text:
using DataFrames, RDatasets, CSV, JLD, HDF5, Query, Match
using Plots, PyPlot, GR, StatsBase, Distributions
using GLM, HypothesisTests, KernelDensity, Clustering

using Statistics, Printf, Dates, DelimitedFiles

 write(stdout, "Help ME!!!\n")

write(stderr,"Error messages should come here.")

#=
read(stdin) will not work in a Jupyter notebook

It is possible to use a terminal

Use File/Open and choose New->Terminal from the pulldowmn menu
Run the Julia REPL

TIP:
On my system (OSX) I have setup the following symbolic link in my
~/bin folder (which is on the executable path)

julia -> /Applications/Julia-1.0.app/Contents/Resources/julia/bin/julia

Try some of the following:- (CR => carriage return, ^D => control-D)

read(stdin,Char)        - type a letter + CR
read(stdin,String)      - type a string + ^D (CR is part of the string)

read(stdin,Int32)       - type a number + CR
read(stdin,UInt32)      - type same number + CR
readline(stdin)         -   "    "    "       "

# However readline will work in Julia
a = readline()

# Because readline() works in Jupyter this is OK.

function getInt()
    s = chomp(readline())
    try
        parse(eltype(1),s)
    catch ex
        println(ex, "\nCan't convert $s to an integer")
    end
end

getInt()

# As is this

function getFloat()
    s = chomp(readline())
    try
        parse(eltype(1.1),s)
    catch ex
        println(ex, "\nCan't convert $s to a float")
    end
end

getFloat()

bb = Array{UInt8}(undef,0)     # This will not work in a notebook

#=

Opening a file has to main syntax:

open(filename::AbstractString; keywords...)
open(filename::AbstractString, [mode::AbstractString])

If neither keywords of a mode string are given then the default is to open the file for a read

=#

# Discuss diffence between s1*s2 vs String(s1,s2)

flip = ENV["HOME"]*"/PacktPub/Alice/jabberwocky.txt";
isfile(flip)

# Classic paradigm
# This prints the first 4 lines of the Jabberwocky
# We need the let/end because of the new scoping rules
let
    fin = open(flip)
    k = 0;
    while !eof(fin)
        ln = readline(fin)
        k = k+1
        (k > 4) ? break : println(ln)
    end
    close(fin)
end

# Because file processing is a common process Julia has an alternative syntax
# Note jabber4.txt is just the first verse of the Jabberwocky, save us having to count
# and break to reduce the output

flp4 = ENV["HOME"]*"/PacktPub/Alice/jabber4.txt";

open(flp4) do pn4
  while !eof(pn4)
     println(readline(pn4))
  end
end

# The file handle is not defined after the loop, so does not need to be closed.
julia> eof(pn4)

#=
Conveniently we can read all a file in to a variable without explicitly opening it
Note that at the data is returned as a byte array (the file may be a binary)
=#
# For a text file we need to convert it

jb4 = String(read(flp4))

# NB: This needs to be split of "\n"s to reproduce the familiar first verse

# There is an additional format:  open(function, filename)

capitalize(f::IOStream) = chomp(uppercase(String(read(f))));
split(open(capitalize, flp4),"\n")

# We need the chomp() to remove the "\n" in the last line,
# otherwise a 5-element array would be returned with the last element "".

# In the previous chapter we used Perl to reverse a lines of a poem.
# Here is a native Julia version

open( flp4) do pn4
    while !eof(pn4)
        println(reverse(chomp(readline(pn4))))
    end
end

# Again we need to chomp the line, reverse it and use println
# Otherwise the "\n" would come at the front

# We can also emulate the 'wc' or or wc() approaches from previously
# In this routine only the words or counted
# Actually this is the most difficult, counting lines and characters are easy

# We will split on white space and also pinctation marks

const PUNCTS =  [' ','\n','\t','-','.',',',':',';','!','?','\'','"'];

# The routine returns a hash (Dict) with the words found as keys and the count as the values,

function wordcount(text)
  wds = split(lowercase(text), PUNCTS; keepempty = false)
  d = Dict()
  for w = wds
     d[w] = get(d,w,0)+1
  end
  return d
end

# Test it on our 4 line Jabberwocky
wordcount(String(read(flp4)))

# Most words only occur once but 'and' & 'the' have a higher frequency

# Wordcount returns a dictionary of the words in a file
#
# So if we collect the values (i.e the counts) and sum them ...
# ... this gives a total for the file.

using Printf
AliceDir = ENV["HOME"]*"/PacktPub/Alice/";  # Note the convenient trailing '/'

# Filter to look just a the '.txt' files in the Alice directory
# We can collect all the value sof the Dict in an array and sum it for the total i the file
#=
for fname in filter!(r"\.txt$", readdir(AliceDir))
  open(AliceDir*fname) do f
    n = sum(collect(values(wordcount(String(read(f))))))
    @printf "%s: %d\n" fname n
  end
end
=#

for fname in readdir(AliceDir)
  if match(r"\.txt$", fname) != nothing
    open(AliceDir*fname) do f
      n = sum(collect(values(wordcount(String(read(f))))))
      @printf "%s: %d\n" fname n
    end
  end
end

# Look at some of the characters in the Hunting of the Snark
# All of the ones in the boat begin with a 'B'

snarkDict = wordcount(String(read(AliceDir*"hunting-the-snark.txt")))

wds = ["baker","banker","barrister","beaver","bellman","boots","butcher"];
for w in wds
  @printf "%12s  => %4d\n" w snarkDict[w]
end

# Set the correct directory and open the Juliaset
# Read the first line for the "magic" number

cd(ENV["HOME"]*"/PacktPub/Chp06");
img = open("Files/juliaset.pgm");
magic = chomp(readline(img))

#=
Check is is a PGM file
If so pick up the next line to get the image size
Open an 'output' file
This will be another PGM file of the same size ...
... so we can write 'magic' and 'params' to the output
=#
if magic == "P5"
  out = open("Files/jsetinvert.pgm", "w");
  println(out, magic);
  params = chomp(readline(img));  # => "800 400 255"
  println(out, params);
# Params splits into strings, we need integers
  (wd, ht, pmax) = parse.(Int64,split(params));
# Create a byte array and read ALL the image data in one call.
  np = wd*ht;
  buf = Array{UInt8,1}(undef,np)
  readbytes!(img, buf, np);
# Invert the gra
  bufX = [UInt8(255 - buf[i]) for i = 1:np]
  write(out,bufX))
  close(out);
else
  error("Not a NetPBM grayscale file")
end

close(img)

# Define a pseudocolor filter
#
function pseudocolor(pix)
  if pix < 64
    pr = UInt8(0); pg = UInt8(0); pb = UInt8(4*pix)
  elseif pix < 128
    pr = UInt8(0); pg = UInt8(min(4*(pix - 64),255)); pb = UInt8(255)
  elseif pix < 192
    pr = UInt8(0); pg = UInt8(255); pb = UInt8(min(4*(192 - pix),255))
  else
    pr = UInt8(min(4*(pix - 192),255))
    pg = UInt8(min(4*(256 - pix),255))
    pb = UInt8(0)
  end
  return (pr, pg, pb)
end

pseudocolor(0xa1)

cd(ENV["HOME"]*"/PacktPub/Chp06");
img = open("Files/juliaset.pgm");
magic = chomp(readline(img))

#=
Same logical as before
Except:-
1. Changing magic number from P5 => P6
2. Applying the pseudocolor filter
3. Writing 3 bytes for each (r,g,b)
=#
if magic == "P5"
  out = open("Files/jsetcolor.ppm", "w");
  println(out, "P6");
  params = chomp(readline(img));  # => "800 400 255"
  println(out, params);
  (wd, ht, pmax) = parse.(Int64,split(params));
  np = wd*ht;
  buf = Array{UInt8,1}(undef,np)
  readbytes!(img, buf, np);
  for j = 1:np
    (r,g,b) = pseudocolor(buf[j]);
    write(out,UInt8(r)); write(out,UInt8(g)); write(out,UInt8(b))
  end
  close(out);
else
  error("Not a NetPBM grayscale file")
end

close(out)
close(img)

# Start by working with some Apple stock prices
using CSV, Statistics, Printf

cd(ENV["HOME"]*"/PacktPub/Chp06/")
aaplcsv = "Files/aapl.csv"; isfile(aaplcsv)

# The CSV.File routine returns the schema
# Notice the unions as the data may have missing values
aapl = CSV.File(aaplcsv)

# Look at the fields comprising the datastructure.
fieldnames(typeof(aapl))

# We are at the first position
aapl.currentrow

# Print the first 5 rows
k = 0;
for r in aapl
    @printf "%s : %.2f\n" r.Date (r.Close - r.Open)
    global k = k + 1
    if k > 5 break end
end

# Check that we are now at the 6th row.
aapl.currentrow

# It is possible to 'pipe' the CSV.File structure to a Dataframe
# Useful to sort the dataframe in place
#
using DataFrames
df = aapl |> DataFrame

sort!(df)

# Dataframes can be queried with the "Queryverse"
# Of which more later
using Query, Dates

x = @from i in df begin
    @where i.Date >= Date(2013,12,20)
    @select {i.Date, i.Open, i.High, i.Low, i.Close}
    @collect DataFrame
end

using DelimitedFiles

cd(ENV["HOME"]*"/PacktPub/Chp06/")
ukhptsv = "Files/UKH-Prcs.tsv"; isfile(ukhptsv)

(ukhpData,ukhpHead) = readdlm(ukhptsv, '\t'; header=true)
ukhpHead

(ukd1, ukd2) = size(ukhpData)

typeof(ukhpData)

using PyPlot

t = collect(1:ukd1);

for i in 1:10
  plot(t,ukhpData[:,i])
end

aaplcsv = "Files/aapl.csv"; isfile(aaplcsv)

(aaplData,aaplHead) = readdlm(aaplcsv, ','; header=true)
aaplHead

aaplData[1:10, 1:6]

using Dates, PyPlot

d0 = Date("2000-01-01")
aaplDate  = reverse(Date.(aaplData[:,1]))
aaplOpen = reverse(Float64.(aaplData[:,2]))
aaplClose = reverse(Float64.(aaplData[:,5]))

const NAAPL = length(aaplDate)
aapl_days = zeros(Int64, NAAPL)

d0 = aaplDate[1]
dt = [(aaplDate[i] - d0).value for i = 1:NAAPL]

plot(dt, aaplClose)
title("Apple Stock - Closing Price")

aaplDifs = [(aaplClose[i] - aaplOpen[i]) for i = 1:NAAPL]
plot(dt,aaplDifs)
title("Apple Stock - Daily Change")

using JLD, HDF5

cd(ENV["HOME"]*"/PacktPub/Chp06/")
h5file = "Files/mydata.h5"
aa = [u + v*rand() for u = 0.5:0.5:10.0, v = 0.5:0.5:6.0]

h5open(h5file, "w") do f
    write(f, "aa", aa)
end

# Alternatively, we can say "@write h5file aa" or else h5write(h5file,"aa",aa)

# This can be read back without the h5open statement (similar to h5write)
# Also we can create a slice of the file on-the-fly

bb = h5read(h5file, "aa", (2:3:14, 4:3:10))

aapljld = "Files/aapl.jld";
rm(aapljld, force=true)

jldopen(aapljld, "w") do fid
  write(fid,"aaplDate",aaplDate)
  write(fid,"aaplClose",aaplClose)
  write(fid,"aaplDifs",aaplDifs)
end

isfile(aapljld)

fid = jldopen(aapljld, "r")

aaDate = read(fid, "aaplDate")
aaDifs = read(fid, "aaplDifs")
close(fid)

fid = jldopen("Files/myaapl.jld", "w")
g = g_create(fid, "aapl")
g["aaplDate"]  = aaplDate
g["aaplClose"] = aaplClose
g["aaplDifs"]  = aaplDifs

dump(fid)

g = fid["aapl"]
adf = g["aaplDifs"]
close(fid)

# Alternatively: adf = fid["aapl/aaplDifs"]

using LightXML
xdoc = parse_file("Files/books.xml");

xtop = LightXML.root(xdoc);
println(LightXML.name(xtop));

xdoc

using Printf

for c in child_nodes(xtop)
    if is_elementnode(c)
        e = XMLElement(c)
        t = find_element(e, "title")
        title = content(t)
        genre = attribute(t, "genre")
        @printf "%30s -:- %s\n" title genre
    end
end

using Dates
for c in child_nodes(xtop)
   if is_elementnode(c)
      e = XMLElement(c)
      t = find_element(e, "title")
      genre = attribute(t, "genre")
      if genre == "Computing"
         a = find_element(e,"author")
         p = find_element(e,"price")
         curr = attribute(p, "currency")
         d = find_element(e,"publish_date")
         dc = DateTime(content(d))
         ds = string(day(dc)," ",monthname(dc)," ",year(dc))
         desc = find_element(e,"description")
         println("Title:      ", content(t))
         println("Author:     " ,content(a))
         println("Date:       " ,ds)
         println("Price:      " ,p ," (", curr, ")")
         println(content(desc),"\n");
      end
   end
end

# https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/00Index.html
# http://juliastats.github.io/StatsBase.jl/latest

using DataFrames, RDatasets

RDatasets.packages()

RDatasets.datasets("MASS")

using Statistics, StatsBase

quakes = dataset("datasets", "quakes");

quakes[1:5,:]

mags  = Float64.(quakes[:Mag]);
depth = Float64.(quakes[:Depth]);

cor(mags,depth)

using Plots
gr

scatter(depth,mags)

describe(mags)

mags = Float64.(quakes[:Mag]);
(m1, m2, m3, m4) = map(x -> round(x,digits=4),
              [mean(mags), std(mags), skewness(mags), kurtosis(mags)]);

using Printf
@printf "Mean: %.4f\nStdV: %.4f\nSkew: %.4f\nKurt: %.4f\n" m1 m2 m3 m4

histmag = fit(Histogram, mags, 4.0:0.1:6.6, closed=:left)

plot(histmag)

using RDatasets

mlmf = dataset("mlmRev","Gcsemv");
size(mlmf)

mlmf[1:5,:]

describe(mlmf)

writtenF = collect(skipmissing(mlmf[mlmf.Gender .== "F", :Written]));
writtenM = collect(skipmissing(mlmf[mlmf.Gender .== "M", :Written]));

(μWM, μWF) = round.((mean(writtenM), mean(writtenF)), digits=3)
(48.286, 45.005)

(σWM, σWF) = (σWM, σWF) = round.((std(writtenM), std(writtenF)),digits=3)
(12.905, 13.535)

(nWM, nWF) = (length(writtenM), length(writtenF))
(706, 997)

σW = sqrt((σWM*σWM)/(nWM - 1) + (σWF*σWF)/(nWF - 1))
tt = round(abs(σWM - σWF)/σW , digits=4)

# p ~ 0.33; 95% ~ 0.06, 90% ~ 0.13
0.9726

using RDatasets, KernelDensity

mlmf = dataset("mlmRev", "Gcsemv");

df = mlmf[completecases(mlmf[[:Written, :Course]]), :]

macro F64(sym)
    quote
      Float64.(skipmissing(Array($sym)))
    end
end

dc = @F64 df[:Course];
dW = @F64 df[:Written];

kdc = kde(dc)

summarystats(dc)

kdw = kde(dw)

summarystats(dw)

using PyPlot

PyPlot.plot(kdc.x, kdc.density)
PyPlot.plot(kdw.x, kdw.density, linestyle="--")

for subdf in groupby(df, :School)
  (size(subdf)[1] > 40) &&
    let
     sch = subdf[:School][1]
     msw = mean(subdf[:Written])
     msc = mean(subdf[:Course])
     @printf "%10s : %8.4f %8.4f\n" sch msw msc
    end
end

using HypothesisTests

df68107 = mlmf[mlmf[:School] .== "68107", :];
df68107cc = df68107[completecases(df68107[[:Written, :Course]]), :];

df68411 = mlmf[mlmf[:School] .== "68411", :];
df68411cc = df68411[completecases(df68411[[:Written, :Course]]), :];

df68107wri = @F64 df68107cc[:Written];
df68411wri = @F64 df68411cc[:Written];

UnequalVarianceTTest(df68107wri, df68411wri)

df68107cou = @F64 df68107cc[:Course];
df68411cou = @F64 df68411cc[:Course];

UnequalVarianceTTest(df68107cou, df68411cou)

using GLM

dw68411ss =  sort(sample(df68411wri,50));
dw68107ss =  sort(sample(df68107wri,50));
dc68411ss =  sort(sample(df68411cou,50));
dc68107ss =  sort(sample(df68107cou,50));

cor(dw68107ss, dw68411ss)

cor(dc68107ss, dc68411ss)

dwf = convert(DataFrame, hcat(dw68107ss, dw68411ss))
names!(dwf, [:s68107, :s68411])

dwf[1:5,:]

dcf = convert(DataFrame, hcat(dc68107ss, dc68411ss))
names!(dcf, [:s68107, :s68411])

dcf[1:5, :]

lm1 = fit(LinearModel, @formula(s68107 ~ s68411), dwf)

lm2 = fit(LinearModel, @formula(s68107 ~ s68411), dcf)

using TimeSeries, MarketData

# Two years of AAPL data (for 2000 / 2001)

ohlcv[1:5]

#=
struct TimeArray{T,N,D<:TimeType,A<:AbstractArray{T,N}} <: AbstractTimeSeries{T,N,D}
    timestamp::Vector{D}
    values::A   # an AbstractArray{T,N}
    colnames::Vector{Symbol}
    meta::ANY
end

=#

meta(ohlc)    # ohlc does not include the volume (v)
colnames(ohlc)
timestamp(ohlc)[end-2:end]
values(ohlc)[end-2:end]
cl[1:3]   # Also op, hi, lo

using PyPlot

PyPlot.plot(timestamp(cl),values(cl))

# Marketdata used to provide routines for getting data from Yahoo (no-defunct)
# Also from FRED (Federal Reserve Bank of St Louis), for US economic data
# see: https://fred.stlouisfed.org
#
# As with the ohlc dataset fred() returns a default series, the US CPI
#
cpi =  fred()
vcat(cpi[1:3],cpi[end-2:end])

dsg10 = fred("DGS10")
vcat(dsg10[1:3],dsg10[end-2:end])

t = collect(1:length(dsg10))
y = values(dsg10)
PyPlot.plot(t,y)

# Temporal.jl is an alternative timeseries package, maintained by JuliaQuant
# and integrating with some packages from the group, such as Quandl and Indicators
# which we will meet later in the book
#
# * Motivated by the `xts` package in R and the `pandas` package in Python *
#
using Temporal, MarketData

md = MarketData   # Ensure no name clash on 'ohlc'
aaplts = TS(values(md.ohlc),timestamp(md.ohlc),colnames(md.ohlc))

using Printf
aapl_spread = aaplts[:High] - aaplts[:Low]
@printf "Average spread : %.3f\n" round(sum(aapl_spread)/length(aapl_spread), digits=3)

describe(aapl_spread.values[:,1])

plot(aapl_spread.index,aapl_spread.values[:,1])


text:
include("activate.jl")

import LinearAlgebra

n = Int(1e5)
k = 1000

xpi = rand(n);
total = sum(xpi);
const xp = xpi ./ total;

f(x) = LinearAlgebra.norm(x - xp)^2

function grad!(storage, x)
    @. storage = 2 * (x - xp)
end

lmo = FrankWolfe.KSparseLMO(40, 1.0);
x00 = FrankWolfe.compute_extreme_point(lmo, zeros(n));

FrankWolfe.benchmark_oracles(x -> f(x), (str, x) -> grad!(str, x), () -> randn(n), lmo; k=100)

println("\n==> Short Step rule - if you know L.\n")

x0 = deepcopy(x00)

@time x, v, primal, dual_gap, trajectorySs = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
);

println("\n==> Adaptive if you do not know L.\n")

x0 = deepcopy(x00)

@time x, v, primal, dual_gap, trajectoryAda = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
);

println("\n==> Agnostic if function is too expensive for adaptive.\n")

x0 = deepcopy(x00)

@time x, v, primal, dual_gap, trajectoryAg = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Agnostic(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
);

data = [trajectorySs, trajectoryAda, trajectoryAg]
label = ["short step", "adaptive", "agnostic"]

FrankWolfe.plot_trajectories(data, label, xscalelog=true)

ext: jl


text:
"""
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
run_powerflow!(sys)
# Passing NLsolve arguments
run_powerflow!(sys, method=:newton)
# Using Finite Differences
run_powerflow!(sys, finite_diff=true)
```
"""
function run_powerflow!(system::PSY.System; finite_diff=false, kwargs...)
    #Save per-unit flag
    settings_unit_cache = deepcopy(system.units_settings.unit_system)
    #Work in System per unit
    PSY.set_units_base_system!(system, "SYSTEM_BASE")
    res = _run_powerflow(system, finite_diff; kwargs...)
    if res.f_converged
        write_powerflow_solution!(system, res.zero)
        @info("PowerFlow solve converged, the results have been stored in the system")
        #Restore original per unit base
        PSY.set_units_base_system!(system, settings_unit_cache)
        return res.f_converged
    end
    @error("The powerflow solver returned convergence = $(res.f_converged)")
    PSY.set_units_base_system!(system, settings_unit_cache)
    return res.f_converged
end

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

## Examples

```julia
res = run_powerflow(sys)
# Passing NLsolve arguments
res = run_powerflow(sys, method=:newton)
# Using Finite Differences
res = run_powerflow(sys, finite_diff=true)
```
"""
function run_powerflow(system::PSY.System; finite_diff=false, kwargs...)
    #Save per-unit flag
    settings_unit_cache = deepcopy(system.units_settings.unit_system)
    #Work in System per unit
    PSY.set_units_base_system!(system, "SYSTEM_BASE")
    res = _run_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.zero)
        #Restore original per unit base
        PSY.set_units_base_system!(system, settings_unit_cache)
        return df_results
    end
    @error("The powerflow solver returned convergence = $(res.f_converged)")
    PSY.set_units_base_system!(system, settings_unit_cache)
    return res.f_converged
end

function _run_powerflow(system::PSY.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(PSY.get_components(PSY.Bus, system)), by=x -> PSY.get_number(x))
    N_BUS = length(buses)

    # assumes the ordering in YPSY.Bus is the same as in the buses.
    Yb = PSY.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 = PSY.get_bustype(buses[ix_t])
            #Set to 0.0 only on connected buses
            if b == PSY.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 == PSY.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 == PSY.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 =
        PSY.get_components(PSY.StaticInjection, system, d -> !isa(d, PSY.ElectricLoad))

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

        P_LOAD_BUS[ix], Q_LOAD_BUS[ix] = _get_load_data(system, b)
        if b.bustype == PSY.BusTypes.REF
            injection_components = PSY.get_components(
                PSY.StaticInjection,
                system,
                d -> PSY.get_available(d) && PSY.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] = PSY.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 == PSY.BusTypes.PV
            x0[state_variable_count] = Q_GEN_BUS[ix]
            x0[state_variable_count + 1] = bus_angle
            state_variable_count += 2
        elseif b.bustype == PSY.BusTypes.PQ
            x0[state_variable_count] = bus_voltage_magnitude
            x0[state_variable_count + 1] = bus_angle
            state_variable_count += 2
        else
            throw(ArgumentError("PSY.Bustype not recognized"))
        end
    end

    @assert state_variable_count - 1 == N_BUS * 2
    bus_types = PSY.get_bustype.(buses)
    function pf!(F::Vector{Float64}, X::Vector{Float64})
        for (ix, b) in enumerate(bus_types)
            if b == PSY.BusTypes.REF
                # When bustype == REFERENCE PSY.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 == PSY.BusTypes.PV
                # When bustype == PV PSY.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 == PSY.BusTypes.PQ
                # When bustype == PQ PSY.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 == PSY.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 == PSY.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 == PSY.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 = PSY.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:
abstract type AbstractServer end
abstract type AbstractFile end

mutable struct File
    path::String
    source::String
    cst::EXPR
    root::Union{Nothing,File}
    server
end

mutable struct FileServer <: AbstractServer
    files::Dict{String,File}
    roots::Set{File}
    symbolserver::Dict{String,SymbolServer.ModuleStore}
end

# Interface spec.
# AbstractServer :-> (has/canload/load/set/get)file, getsymbolserver
# AbstractFile :-> (get/set)path, (get/set)root, (get/set)cst, scopepass, (get/set)server

hasfile(server::FileServer, path::String) = haskey(server.files, path)
canloadfile(server, path) = isfile(path)
function setfile(server::FileServer, path::String, file::File)
    server.files[path] = file
end
getfile(server::FileServer, path::String) = server.files[path]
function loadfile(server::FileServer, path::String)
    source = read(path, String)
    cst = CSTParser.parse(source, true)
    f = File(path, source, cst, nothing, server)
    setroot(f, f)
    setfile(server, path, f)
    return getfile(server, path)
end
getsymbolserver(server::FileServer) = server.symbolserver

function scopepass(file)
    server = file.server
    setscope!(getcst(file), Scope(nothing, Dict(), Dict{String,Any}("Base" => getsymbolserver(server)["Base"], "Core" => getsymbolserver(server)["Core"]), false))
    state = State(file, scopeof(getcst(file)), false, false, false, EXPR[], server)
    state(getcst(file))
    for uref in state.urefs
        s = retrieve_delayed_scope(uref)
        if s !== nothing
            resolve_ref(uref, s, state)
        end
    end
end

getpath(file::File) = file.path
function setpath(file::File, path)
    file.path = path
    return file
end

getroot(file::File) = file.root
function setroot(file::File, root::File)
    file.root = root
    return file
end

getcst(file::File) = file.cst
function setcst(file::File, cst::EXPR)
    file.cst = cst
    return file
end

getserver(file::File) = file.server
function setserver(file::File, server::FileServer)
    file.server = server
    return file
end

function Base.display(f::File)
    println(f.path)
end

function Base.display(s::FileServer)
    n = length(s.files)
    println(n, "-file Server")
    cnt = 0
    for (p,f) in s.files
        cnt+=1
        println(" ", p)
        cnt > 10 && break
    end
end

function get_path(x::EXPR)
    if typof(x) === Call && length(x.args) == 4
        parg = x.args[3]
        if CSTParser.is_lit_string(parg)
            path = CSTParser.str_value(parg)
            return normpath(path)
        elseif typof(parg) === Call && typof(parg.args[1]) === IDENTIFIER && CSTParser.str_value(parg.args[1]) == "joinpath"
            path = ""
            for i = 2:length(parg.args)
                arg = parg.args[i]
                if typof(arg) === PUNCTUATION
                elseif CSTParser.is_lit_string(arg)
                    path = string(path, "/", CSTParser.str_value(arg))
                else
                    return ""
                end
            end
            return normpath(path)
        end
    end
    return ""
end


text:
using DashBootstrapComponents

jumbotron = html_div(
    dbc_container(
        [
            html_h1("Jumbotron", className = "display-3"),
            html_p(
                "Use Containers to create a jumbotron to call attention to " *
                "featured content or information.",
                className = "lead",
            ),
            html_hr(className = "my-2"),
            html_p(
                "Use utility classes for typography and spacing to suit the " *
                "larger container.",
            ),
            html_p(dbc_button("Learn more", color = "primary"), className = "lead"),
        ],
        fluid = true,
        className = "py-3",
    ),
    className = "p-3 bg-light rounded-3",
)

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
module RocketIgnoreOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: ignore()" begin

    println("Testing: operator ignore()")

    run_proxyshowcheck("Ignore", ignore(0))

    run_testset([
        (
            source = from(1:3) |> ignore(-1),
            values = @ts([ 1, 2, 3, c ])
        ),
        (
            source = from(1:3) |> ignore(0),
            values = @ts([ 1, 2, 3, c ])
        ),
        (
            source = from(1:3) |> ignore(1),
            values = @ts([ 2, 3, c ])
        ),
        (
            source = from(1:3) |> ignore(2),
            values = @ts([ 3, c ])
        ),
        (
            source = from(1:3) |> ignore(3),
            values = @ts(c)
        ),
        (
            source = from(1:3) |> ignore(4),
            values = @ts(c)
        ),
        (
            source = from(1:3) |> async(0) |> ignore(2),
            values = @ts([ 3 ] ~ c)
        ),
        (
            source = completed(Int) |> ignore(2),
            values = @ts(c)
        ),
        (
            source = never(Int) |> ignore(2),
            values = @ts()
        ),
        (
            source = faulted("e") |> ignore(2),
            values = @ts(e("e"))
        )
    ])

end

end


text:
"""
    plotDSM(DSM,label)

Plot the binary DSM.

# Arguments
- `DSM`: adjacency matrix.
- `label`: label for elements in this adjacency matrix.

# Examples

```jldoctest
Julia> DSM = [0 1 0 0 0 1 0;
            0 0 0 1 0 0 0;
            1 0 0 0 0 0 1;
            0 0 0 0 1 0 0;
            0 1 0 0 0 1 0;
            0 0 1 0 0 0 0;
            1 0 0 0 1 0 0];
julia> label = ["A", "B", "C", "D", "E", "F", "G"];
julia> plotDSM(DSM,label)
"hugahuga"
```

"""

function plotDSM(DSM,label)
    (DSMr,DSMc) = size(DSM);

    @svg begin

        CellSizeRow = min(50,450/(DSMr+5));
        CellSizeColumn = min(50,450/(DSMc+5));
        t = Table(DSMr+1,DSMc+1, CellSizeRow, CellSizeColumn);

        sethue("black")

        for r in 2 : DSMr+1
            for c in 2 : DSMc+1
                box(t[r,c], t.colwidths[t.currentcol], t.rowheights[t.currentrow], :stroke);
                if r == c
                    sethue("grey")
                    box(t[r,c], t.colwidths[t.currentcol], t.rowheights[t.currentrow], :fill);
                    sethue("white")
                    textfit(label[r-1], BoundingBox(box(t[r,c],CellSizeRow, CellSizeColumn)))
                    sethue("black")
                end

                if DSM[r-1,c-1] != 0
                    circle(t[r,c], CellSizeColumn/3, :fill)
                end
            end
        end

        for r in 1 : DSMr
            textfit(label[r], BoundingBox(box(t[r+1,1],CellSizeRow, CellSizeColumn)))
        end

        for c in 1 : DSMc
            textfit(label[c], BoundingBox(box(t[1,c+1],CellSizeRow, CellSizeColumn)))
        end

    end
end


text:
module Coord

using Gadfly
using Compose

import Gadfly.Maybe
export cartesian

# Cartesian coordinates with position given by the x and y (and similar)
# aesthetics.
#
# Args:
#   xvars: Aesthetics to consider when choosing x bounds.
#   yvars: Aesthetics to consider when choosing y bounds.
#   xmin, xmax, ymin, ymax:
#     Force a particular x or y bound, rather than trying to choose it
#     from the data.
#
#
type Cartesian <: Gadfly.CoordinateElement
    xvars::Vector{Symbol}
    yvars::Vector{Symbol}
    xmin
    xmax
    ymin
    ymax

    function Cartesian(
            xvars=[:x, :x_viewmin, :x_viewmax, :x_min, :x_max],
            yvars=[:y, :y_viewmin, :y_viewmax, :y_min, :y_max, :middle,
                   :lower_hinge, :upper_hinge, :lower_fence, :upper_fence, :outliers];
            xmin=nothing, xmax=nothing, ymin=nothing, ymax=nothing)
        new(xvars, yvars, xmin, xmax, ymin, ymax)
    end
end

const cartesian = Cartesian()

# Produce a canvas with suitable cartesian coordinates.
#
# Args:
#   coord: cartesian coordinate instance.
#
# Returns:
#   A compose Canvas.
#
function apply_coordinate(coord::Cartesian, aess::Gadfly.Aesthetics...)
    xmin = Inf
    xmax = -Inf
    for var in coord.xvars
        for aes in aess
            if getfield(aes, var) === nothing
                continue
            end

            for val in getfield(aes, var)
                if val < xmin
                    xmin = val
                end

                if val > xmax
                    xmax = val
                end
            end
        end
    end

    ymin = Inf
    ymax = -Inf
    for var in coord.yvars
        for aes in aess
            if getfield(aes, var) === nothing
                continue
            end

            # Outliers is an odd aesthetic that needs special treatment.
            if var == :outliers
                for vals in aes.outliers
                    for val in vals
                        if val < ymin
                            ymin = val
                        end

                        if val > ymax
                            ymax = val
                        end
                    end
                end

                continue
            end

            for val in getfield(aes, var)
                if val < ymin
                    ymin = val
                end

                if val > ymax
                    ymax = val
                end
            end
        end
    end

    xmin = coord.xmin === nothing ? xmin : coord.xmin
    xmax = coord.xmax === nothing ? xmax : coord.xmax
    ymin = coord.ymin === nothing ? ymin : coord.ymin
    ymax = coord.ymax === nothing ? ymax : coord.ymax

    # A bit of kludge. We need to extend a bit to be able to fit bars when
    # using discrete scales. TODO: Think more carefully about this. Is there a
    # way for the geometry to let the coordinates know that a little extra room
    # is needed to draw everything?
    if all([aes.x === nothing || typeof(aes.x) == Array{Int64, 1} for aes in aess])
        xmin -= 0.5
        xmax += 0.5
        xpadding = 0
    else
        xpadding = 0.03 * (xmax - xmin)
    end

    if all([aes.y === nothing || typeof(aes.y) == Array{Int64, 1} for aes in aess])
        ymin -= 0.5
        ymax += 0.5
        ypadding = 0
    else
        ypadding = 0.03 * (ymax - ymin)
    end

    width  = xmax - xmin + 2xpadding
    height = ymax - ymin + 2ypadding

    canvas(unit_box=Units(xmin - xpadding, ymax + ypadding, width, -height))
end

end # module Coord


text:
module TestMath

using Test
using RLEVectors

@testset begin

    x = RLEVector(["a","b","c","b"],[2,4,6,10])
    m = countmap(x)
    @test collect(m) == ["c" => 2, "b" => 6, "a" => 2]
    x = RLEVector([4,5,6,5],cumsum([10,5,3,6]))
    @test mode(x) == 5

end # testset

end # module


text:
abstract type AbstractAggregation end
abstract type AbstractAggregator end
# ------------------------------------------------------------
include("cut_collection.jl")
include("distance_measures.jl")
include("selection_rules.jl")
include("cluster_rules.jl")
include("lock_rules.jl")
include("no_aggregation.jl")
include("partial_aggregation.jl")
include("dynamic_aggregation.jl")
include("cluster_aggregation.jl")
include("hybrid_aggregation.jl")


text:
using FANN
using Base.Test

# test simple linear problem for which we know the solution
net = ANN([2, 1], [:linear], b=0.1)

X = [ 0.0 1.0 0.0 ;
	  0.0 0.0 1.0 ]
y = [ 0.0 1.0 1.0 ]

dset = DataSet(X, y)

# train
train!(net, dset,
	   max_epochs=1000,
	   epochs_between_reports=0,
	   desired_error=1e-8)

# there is somethign strange here, since the network stops training
# but the mean square error is not less than the desired error. Where
# is this bug from?
out = predict(net, X)
@test sumabs2(out - y) < 100*1e-8


text:
using TerminalLoggers
using Test

using Logging:
    LogLevel, BelowMinLevel, Debug, Info, Warn, Error, AboveMaxLevel,
    shouldlog, handle_message, min_enabled_level, catch_exceptions

include("TerminalLogger.jl")
include("StickyMessages.jl")


text:
# Check that package loads
using StatsKit

using Pkg

# TODO: run all tests in the same process after 'using StatsKit'
#       to check that they pass when all packages are loaded at the same time
Pkg.test(["Bootstrap", "CategoricalArrays", "Clustering", "CSV", "DataFrames",
          "Distances", "Distributions", "GLM", "HypothesisTests", "KernelDensity",
          "Loess", "MultivariateStats", "StatsBase", "TimeSeries"])


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    bspkntins(d::I,c::Matrix{F},k::Vector{F},u::Vector{F}
             ) where {I<:Integer,F<:AbstractFloat}

Insert knots into a univariate B-Spline.

# Arguments:
- `d`: Degree of the B-Spline.
- `c`: Control points, matrix of size `(dim,nc)`.
- `k`: Knot sequence, vector of size `nk`.
- `u`: Row vector of knots to be inserted, size `nu`.

# Output:
- `ic`: Control points of the new B-Spline, of size `(dim,nc+nu)`.
- `ik`: Knot vector of the new B-Spline, of size `(nk+nu)`.

# Examples:
```julia
julia> ic,ik = bspkntins(d,c,k,u)
```

# Reference:
- Algorithm A5.4 from 'The NURBS BOOK' pg164 ([1997Piegl](@cite))
"""
function bspkntins(d::I,c::Matrix{F},k::Vector{F},u::Vector{F}
                   ) where {I<:Integer,F<:AbstractFloat}
# TODO: the original C function passed c, k, ic, ik and u as reference
# TODO: this function has optimization potential

mc, nc = size(c);
nu = length(u);
nk = length(k);

ic = zeros(F,(mc,nc+nu));
ik = zeros(F,nk+nu);

n = nc - 1;
r = nu - 1;

m = n + d + 1;
a = findspan(n, d, u[1], k);
b = findspan(n, d, u[r+1], k);
b += 1;

for q in 0:mc-1
  for j in 0:a-d; ic[q+1,j+1] = c[q+1,j+1]; end # for j
  for j in b-1:n; ic[q+1,j+r+2] = c[q+1,j+1]; end # for j
end # for q

for j in 0:a; ik[j+1] = k[j+1]; end # for j
for j in b+d:m; ik[j+r+2] = k[j+1]; end # for j

i = b + d - 1;
s = b + d + r;

for j in r:-1:0
  while u[j+1] <= k[i+1] && i > a
    for q in 0:mc-1
      ic[q+1,s-d] = c[q+1,i-d];
    end # for q
    ik[s+1] = k[i+1];
    s -= 1;
    i -= 1;
  end # for q

  for q in 0:mc-1
    ic[q+1,s-d] = ic[q+1,s-d+1];
  end # for q

  for l in 1:d
    ind = s - d + l;
    alfa = ik[s+l+1] - u[j+1];
    if abs(alfa) == 0
      for q in 0:mc-1
        ic[q+1,ind] = ic[q+1,ind+1];
      end
    else
      alfa /= (ik[s+l+1] - k[i-d+l+1]);
      for q=0:mc-1
        ic[q+1,ind] = alfa*ic[q+1,ind] + (1-alfa)*ic[q+1,ind+1];
      end # for q
    end # if
  end # for l

  ik[s+1] = u[j+1];
  s -= 1;
end # for j
return ic, ik;
end # bspkntins


text:
export # White matter
       CorpusCallosum,
       ForcepsMajor,
       ForcepsMinor,
       CingulateCingulum,  # TODO I dislike the label cingulate cingulum very much
       ParahippocampalCingulum,
       Fornix,
       InferiorFrontalOccipitalFasciculus,
       InferiorLongitudinalFasciculus,
       CorticospinalTract,
       SuperiorLongitudinalFasciculus,
       ArcuateFasciculus,
       AnteriorThalamicRadiation,
       Uncinate

#include("neurotypes.jl")
include("cerebellum.jl")
include("cingulate.jl")
include("fibers.jl")
include("temporal_lobe.jl")
include("frontal_lobe.jl")
include("occipital_lobe.jl")
include("parietal_lobe.jl")
include("diencephalon.jl")


text:
# # Visualization of Frank-Wolfe running on a 2-dimensional polytope

# This example provides an intuitive view of the Frank-Wolfe algorithm
# by running it on a polyhedral set with a quadratic function.
# The Linear Minimization Oracle (LMO) corresponds to a call to a generic simplex solver from `MathOptInterface.jl` (MOI).

# ## Import and setup

# We first import the necessary packages, including Polyhedra to visualize the feasible set.

using LinearAlgebra
using FrankWolfe

import MathOptInterface
const MOI = MathOptInterface
using GLPK

using Polyhedra
using Plots

# We can then define the objective function,
# here the squared distance to a point in the place, and its in-place gradient.

n = 2
y = [3.2, 0.5]

function f(x)
    return 1 / 2 * norm(x - y)^2
end
function grad!(storage, x)
    @. storage = x - y
end

# ## Custom callback
#
# FrankWolfe.jl lets users define custom callbacks to record information about each iteration.
# In that case, the callback will copy the current iterate `x`, the current vertex `v`, and the current step size `gamma`
# to an array thanks to a closure.
# We then declare the array and the callback over this array.
# Each iteration will then push to this array.

function build_callback(trajectory_arr)
    return function callback(state)
        return push!(trajectory_arr, (copy(state.x), copy(state.v), state.gamma))
    end
end

iterates_information_vector = []
callback = build_callback(iterates_information_vector)

# ## Creating the Linear Minimization Oracle
# The LMO is defined as a call to a linear optimization solver, each iteration resets the objective and calls the solver.
# The linear constraints must be defined only once at the beginning and remain identical along iterations.
# We use here MathOptInterface directly but the constraints could also be defined with JuMP or Convex.jl.

o = GLPK.Optimizer()
x = MOI.add_variables(o, n)

## −x + y ≤ 2
c1 = MOI.add_constraint(
    o,
    -1.0x[1] + x[2],
    MOI.LessThan(2.0),
)

## x + 2 y ≤ 4
c2 = MOI.add_constraint(
    o,
    x[1] + 2.0x[2],
    MOI.LessThan(4.0),
)

## −2 x − y ≤ 1
c3 = MOI.add_constraint(
    o,
    -2.0x[1] - x[2],
    MOI.LessThan(1.0),
)

## x − 2 y ≤ 2
c4 = MOI.add_constraint(
    o,
    x[1] - 2.0x[2],
    MOI.LessThan(2.0),
)

## x ≤ 2
c5 = MOI.add_constraint(
    o,
    x[1] + 0.0x[2],
    MOI.LessThan(2.0),
)

# The LMO is then built by wrapping the current MOI optimizer

lmo_moi = FrankWolfe.MathOptLMO(o)

# ## Calling Frank-Wolfe
# We can now compute an initial starting point from any direction
# and call the Frank-Wolfe algorithm.
# Note that we copy `x0` before passing it to the algorithm because it is modified in-place by `frank_wolfe`.

x0 = FrankWolfe.compute_extreme_point(lmo_moi, zeros(n))

xfinal, vfinal, primal_value, dual_gap, traj_data = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo_moi,
    copy(x0),
    line_search = FrankWolfe.Adaptive(),
    max_iteration = 10,
    epsilon=1e-8,
    callback=callback,
    verbose=true,
    print_iter=1,
)

# We now collect the iterates and vertices across iterations.

iterates = Vector{Vector{Float64}}()
push!(iterates, x0)
vertices = Vector{Vector{Float64}}()
for s in iterates_information_vector
    push!(iterates, s[1])
    push!(vertices, s[2])
end

# ## Plotting the algorithm run

# We define another method for `f` adapted to plot its contours.

function f(x1, x2)
    x = [x1, x2]
    return f(x)
end

xlist = collect(range(-1, 3, step = 0.2))
ylist = collect(range(-1, 3, step = 0.2))

X = repeat(reshape(xlist, 1, :), length(ylist), 1)
Y = repeat(ylist, 1, length(xlist))

# The feasible space is represented using Polyhedra.

h = HalfSpace([-1, 1], 2) ∩
    HalfSpace([1, 2], 4) ∩
    HalfSpace([-2, -1], 1) ∩
    HalfSpace([1, -2], 2) ∩
    HalfSpace([1, 0], 2)

p = polyhedron(h)

p1 = contour(xlist, ylist, f, fill = true, line_smoothing = 0.85)
plot(p1, opacity = 0.5)
plot!(p, ratio = :equal, opacity = 0.5, label = "feasible region", framestyle = :zerolines, legend = true, color=:blue);

# Finally, we add all iterates and vertices to the plot.
colors = ["gold", "purple", "darkorange2", "firebrick3"]
iterates = unique!(iterates)
for i in 1:3
    scatter!([iterates[i][1]], [iterates[i][2]], label = string("x_", i - 1), markersize = 6, color = colors[i])
end
scatter!([last(iterates)[1]], [last(iterates)[2]], label = string("x_", length(iterates) - 1), markersize = 6, color = last(colors))

# plot chosen vertices
scatter!([vertices[1][1]], [vertices[1][2]], m = :diamond, markersize = 6, color = colors[1], label = "v_1")
scatter!([vertices[2][1]], [vertices[2][2]], m = :diamond, markersize = 6, color = colors[2], label = "v_2", legend = :outerleft, colorbar = true)


text:
"""
struct containing grids used in a simulation

# Examples

```jldoctest
julia> using JultraDark

julia> len = 1;

julia> resol = 16;

julia> PencilGrids(len, resol);

```
"""
struct PencilGrids
    "Real space distance array"
    dist
    "Fourier space postition array"
    k
    "Fourier space postition array for use with `rfft`"
    rk
    "ψ field"
    ψx
    "ψ field in Fourier space"
    ψk
    "density field ρ"
    ρx
    "density field ρ in Fourier space"
    ρk
    "gravitational potential field Φ"
    Φx
    "gravitational potential field Φ in fourier space"
    Φk
    fft_plan
    rfft_plan

    function PencilGrids(dist, k, rk, ψx, ψk, ρx, ρk, Φx, Φk, fft_plan, rfft_plan)
        n_dims = 3
        resol_tuple = size_global(dist)
        resol_tuple_realfft = (size_global(dist)[1] ÷ 2 + 1, size_global(dist)[2], size_global(dist)[3])

        for var in [dist, k, rk, ψx, ψk, ρx, ρk, Φx, Φk]
            @assert(ndims(var) == n_dims)
        end

        for var in [dist, k, ψx, ψk, ρx, Φx]
            @assert(size_global(var) == resol_tuple)
        end

        for var in [rk, ρk, Φk]
            @assert(size_global(var) == resol_tuple_realfft)
        end

        new(dist, k, rk, ψx, ψk, ρx, ρk, Φx, Φk, fft_plan, rfft_plan)
    end
end

"""
    PencilGrids(length::Real, resol::Integer)

Constructor for `PencilGrids`

Create an empty grid with length `length` and resolution `resol`.  Uses `PencilFFTs` to create `PencilArrays`.

# Examples

```jldoctest
julia> using JultraDark

julia> PencilGrids(1.0, 64);

```
"""
function PencilGrids(length::Real, resol::Integer)::PencilGrids

    resol_tuple = (resol, resol, resol)
    resol_tuple_realfft = (resol ÷ 2 + 1, resol, resol)

    # MPI topology information
    comm = MPI.COMM_WORLD  # we assume MPI.Comm_size(comm) == 12
    Nproc = MPI.Comm_size(comm)

    # Let MPI_Dims_create choose the decomposition.
    proc_dims = let pdims = zeros(Int, 2)
        MPI.Dims_create!(Nproc, pdims)
        pdims[1], pdims[2]
    end

    FFTW.set_num_threads(Threads.nthreads())
    # Plan a 3D complex-to-complex (c2c) FFT.
    fft_plan = PencilFFTPlan(resol_tuple, Transforms.FFT(), proc_dims, comm)

    # Allocate ψx and ψk arrays
    ψx = allocate_input(fft_plan)
    ψk = allocate_output(fft_plan)

    # Plan a 3D real-to-complex (r2c) FFT.
    rfft_plan = PencilFFTPlan(resol_tuple, Transforms.RFFT(), proc_dims, comm)

    ρx = allocate_input(rfft_plan)
    ρk = allocate_output(rfft_plan)

    Φx = allocate_input(rfft_plan)
    Φk = allocate_output(rfft_plan)

    dist = allocate_input(rfft_plan)
    dist_glob = global_view(dist)
    _dist = dist_array(length, resol)
    for I in CartesianIndices(dist_glob)
        dist_glob[I] = _dist[I]
    end

    k = similar(ψk, Float64)
    k_glob = global_view(k)
    _k = k_norm((length, length, length), (resol, resol, resol))
    for I in CartesianIndices(k_glob)
        k_glob[I] = _k[I]
    end

    rk = similar(ρk, Float64)
    rk_glob = global_view(rk)
    _rk = rk_norm((length, length, length), (resol, resol, resol))
    for I in CartesianIndices(rk_glob)
        rk_glob[I] = _rk[I]
    end

    PencilGrids(
        dist,
        k, rk,
        ψx, ψk,
        ρx, ρk,
        Φx, Φk,
        fft_plan, rfft_plan,
    )
end

"""
    PencilGrids(ψx::Array{Complex{Float64}}, length::Real)

Constructor for `PencilGrids`

Create a grid with given ψ field, length `length` and resolution inferred
from `ψx`

# Examples

Can be contructed from a ψ field and box length,
```jldoctest
julia> using JultraDark

julia> ψ = zeros(Complex{Float64}, 16, 16, 16);

julia> len = 1;

julia> PencilGrids(ψ, len);

```
"""
function PencilGrids(ψx::Array{Complex{Float64}}, length::Real)::PencilGrids
    @assert(
        ndims(ψx) == 3,
        "Invalid ψ: only three dimensions supported"
    )
    @assert(
        size(ψx, 1) == size(ψx, 2) == size(ψx, 3),
        "Invalid ψ: heterogenous resolutions"
    )

    resol = size(ψx, 1)
    grids = PencilGrids(length, resol)

    ψx_glob = global_view(grids.ψx)
    for I in CartesianIndices(ψx_glob)
        ψx_glob[I] = ψx[I]
    end

    grids.ρx .= abs2.(grids.ψx)
    a_init = 1  # TODO

    grids.ρk .= grids.rfft_plan * grids.ρx
    grids.Φk .= -4 * π * grids.ρk ./ (a_init * grids.rk.^2)
    grids.Φk[1, 1, 1] = 0
    grids.Φx .= grids.rfft_plan \ grids.Φk

    grids

end


text:
"""
    add!!(x, y)

Returns `x+y`, potentially mutating `x` in-place to hold this value.
This avoids allocations when `x` can be mutated in this way.
"""
add!!(x, y) = x + y

"""
    add!!(x, t::ImplacableThunk)

The specialization of `add!!` for [`InplaceableThunk`](@ref) promises to only call
`t.add!` on `x` if `x` is suitably mutable; otherwise it will be out of place.
"""
function add!!(x, t::InplaceableThunk)
    return if is_inplaceable_destination(x)
        if !debug_mode()
            t.add!(x)
        else
            debug_add!(x, t)
        end
    else
        x + t
    end
end

function add!!(x::AbstractArray{<:Any, N}, y::AbstractArray{<:Any, N}) where N
    return if is_inplaceable_destination(x)
        x .+= y
    else
        x + y
    end
end

"""
    is_inplaceable_destination(x) -> Bool

Returns true if `x` is suitable for for storing inplace accumulation of gradients.
For arrays this boils down `x .= y` if will work to mutate `x`, if `y` is an appropriate
differential.
Wrapper array types do not need to overload this if they overload `Base.parent`, and are
`is_inplaceable_destination` if and only if their parent array is.
Other types should overload this, as it defaults to `false`.
"""
is_inplaceable_destination(::Any) = false
is_inplaceable_destination(::Array) = true
is_inplaceable_destination(::SparseVector) = true
is_inplaceable_destination(::SparseMatrixCSC) = true
is_inplaceable_destination(::BitArray) = true
function is_inplaceable_destination(x::AbstractArray)
    p = parent(x)
    p === x && return false  # no parent
    # basically all wrapper types delegate `setindex!` to their `parent` after some
    # processing and so are mutable if their `parent` is.
    return is_inplaceable_destination(p)
end

# Hermitian and Symmetric are too fussy to deal with right now
# https://github.com/JuliaLang/julia/issues/38056
# TODO: https://github.com/JuliaDiff/ChainRulesCore.jl/issues/236
is_inplaceable_destination(::LinearAlgebra.Hermitian) = false
is_inplaceable_destination(::LinearAlgebra.Symmetric) = false

function debug_add!(accumuland, t::InplaceableThunk)
    returned_value = t.add!(accumuland)
    if returned_value !== accumuland
        throw(BadInplaceException(t, accumuland, returned_value))
    end
    return returned_value
end

struct BadInplaceException <: Exception
    ithunk::InplaceableThunk
    accumuland
    returned_value
end

function Base.showerror(io::IO, err::BadInplaceException)
    println(io, "`add!!(accumuland, ithunk))` did not return an updated accumuland.")
    println(io, "ithunk = $(err.ithunk)")
    println(io, "accumuland = $(err.accumuland)")
    println(io, "returned_value = $(err.returned_value)")

    if err.accumuland == err.returned_value
        println(
            io,
            "Which in this case happenned to be equal. But they are not the same object."
        )
    end
end


text:
using Test, OrdinaryDiffEq, HyperbolicDiffEq, DiffEqCallbacks

function runtest()
    xmin = -1.0
    xmax = +0.5
    p = 1
    N = 2^8
    tspan = (0., 1)
    ode_solver = OrdinaryDiffEq.Euler()
    dt(p,N) = 1/((p^2+1)*N)
    u01func(x) = x < 0 ? 1.5 : -1.725862; u02func(x) = x < 0 ? 0.0 : 1.276293
    u0func(x) = KeyfitzKranzerVar(u01func(x), u02func(x))

    mesh = UniformPeriodicMesh1D(xmin, xmax, N)
    basis = LobattoLegendre(p)
    balance_law = KeyfitzKranzer()
    fvol = EnergyConservativeFlux()
    fnum = FluxPlusDissipation(EnergyConservativeFlux(), LocalLaxFriedrichsDissipation())
    semidisc = UniformPeriodicFluxDiffDisc1D(balance_law, mesh, basis, fvol, fnum)
    ode = semidiscretise(semidisc, u0func, tspan)

    save_func = (u,t,integrator) -> integrate(u->IntegralQuantitiesKeyfitzKranzer(u,balance_law), u, integrator.f.f)
    saved_values = SavedValues(Float64, IntegralQuantitiesKeyfitzKranzer{Float64})
    saving = SavingCallback(save_func, saved_values, saveat=range(tspan..., length=10^3))
    stepsize = StepsizeLimiter((u,p,t)->max_dt(t,u,semidisc), safety_factor=1.0)

    sol = solve(ode, ode_solver, dt=dt(p,N), save_everystep=false, callback=CallbackSet(stepsize,saving))

    u1 = map(q->q.u1, saving.affect!.saved_values.saveval)
    u2 = map(q->q.u2, saving.affect!.saved_values.saveval)
    entropy = map(q->q.entropy, saving.affect!.saved_values.saveval)
    @test u1[1] ≈ u1[end]
    @test u2[1] ≈ u2[end]
    @test entropy[1] >= entropy[end]

    nothing
end

runtest()