text:
### A Nbdev script file ###
### Autogenerated file. Don't modify. ###

module Export
#export
import Pluto: Notebook, Cell, load_notebook_nobackup

#export
function processMdChunk(mdChunk::String)
	chopMdMarker(mdChunk)
end

#export
begin
"""
> struct Nucleus--> This is the lowest entity of a code notebook. This type represents the code cell in a Pluto notebook.
> * Fields:
>   * code--> String which makes up a code cell
"""
Base.@kwdef mutable struct Nucleus
	code::String=""
	end

"""
> Nucleus(code)-->Getter for accessing the constructer of the Nucleus type
"""
Nucleus(code) = Nucleus(code = code)
end

#export
begin
"""
> struct Nb--> a notebook having nuclei made up of code cells only.
> * Fields:
>   * nuclei--> Array of Nucleus type
>   * name--> name of notebook
"""
mutable struct Nb
	nuclei::Array{Nucleus,1}
	name::AbstractString
end

"""
> nuclei(nuclei::Array)--> Helper to collect an array of Nucleus.
"""
nuclei(nuclei::Array) = nuclei

"""
> nuclei(nuclei::Array)--> Helper to collect an name of a notebook.
"""
name(name::AbstractString) = name
end

#export
begin
"""
Headers that would appear in a code file
"""
const _header = "### A Nbdev script file ###"
const _sub_header = "### Autogenerated file. Don't modify. ###"
end

#export
"""
This had to be done because having the cell iddelimeter as a constant gives an error during parsing the notebook
"""
cell_id_delimiter=string("#"," ╔═╡ ")

#export
"""
marks the end of a cell
"""
const _cell_suffix = "\n\n"

#export
"""
> struct ScrubbedNotebook--> Represents a notebook from which all but code to be exported are scrubbed off.
> * Fields-->
>   * cells--> Vector of strings.
"""
mutable struct ScrubbedNotebook
	cells::Array{String, 1}
end

#export
"""
> collect_nuclei(notebook::ScrubbedNotebook, marker)--> Reads a notebook, collects code from code cells and returns an array of Nucelus.
"""
function collect_nuclei(notebook::ScrubbedNotebook, marker)
	notebook_cells=notebook.cells
	nuclei=[]

    for i in 1:length(notebook_cells)
    	raw_code=notebook_cells[i]
    	nucleus=Nucleus(raw_code)

    	if startswith(raw_code, marker)
            push!(nuclei,nucleus)
        end

    end
	 nuclei
end

#export
begin
"""
>load_scrubbed_nb(io, path)::ScrubbedNotebook--> These are modification of the Pluto.load_notebook methods. Scrubs the notebook of all stuff but the code to be exported.
"""
function load_scrubbed_nb(io, path)::ScrubbedNotebook
    collected_cells = []

	# ignore first bits of file
    readuntil(io,cell_id_delimiter)

    last_read = ""
    while !eof(io)
        cell_id_str = String(readline(io))
        if cell_id_str == "Cell order:"
            break
        else
            code_raw = String(readuntil(io, cell_id_delimiter))
            # change Windows line endings to Linux
            code_normalised = replace(code_raw, "\r\n" => "\n")
            # remove the cell appendix
            code = code_normalised[1:prevind(code_normalised, end, length(_cell_suffix))]

            read_cell = code
            push!(collected_cells, read_cell)
        end
    end

	ScrubbedNotebook(collected_cells)
end

function load_scrubbed_nb(path::String)::ScrubbedNotebook
    local loaded
    open(path, "r") do io
        loaded = load_scrubbed_nb(io, path)
    end
    loaded
end
end

#export
begin
"""
> load_nb(filename::String, marker::String)--> High level function which calls _load_nb
"""
load_nb(filename::String, marker::String) = _load_nb(filename, marker)

"""
> _load_nb(filename::String, marker::String)--> creates a scrubbed notebook and returns a curated Nb type having code to be exported.
"""
function _load_nb(filename::String, marker::String)
		#notebook=load_notebook_nobackup(filename)
		scrubbedNotebook=load_scrubbed_nb(normpath(joinpath(@__FILE__,"..","..", "nbs", filename)))
		collected_nuclei=nuclei(collect_nuclei(scrubbedNotebook, marker))

		#if marker=="md"
		#	notebook
		#else
		Nb(collected_nuclei, filename)
		#end
end

end

#export
"""
> strip --> Helper fucntion to replace a substring y in a string x with blank.
"""
strip=(x::String,y) -> replace(x, y=>"")

#export
begin
"""
> save_nb(io, notebook)--> Reads the supplied notebook and creates an io and writes stuffs like the module name and the content to the created io.
"""
function save_nb(io, notebook)
    println(io, _header)
    println(io, _sub_header)
    println(io, "")
    println(io, "module $(uppercasefirst(strip(strip(notebook.name, r"[0-9_]"), r".jl")))")
	for nucleus in notebook.nuclei
			println(io, nucleus.code*"\n")
    end
	#notebook
	print(io, "end")
end

"""
> save_nb(notebook::Nb, path::String)--> Creates a file in the supplied path with the name in the NB type.
"""
function save_nb(notebook::Nb, path::String)
	file_name=uppercasefirst(strip(notebook.name, r"[0-9_]"))
	open(joinpath(path, file_name), "w") do io
        save_nb(io, notebook)
    end
end
end

#export
begin
"""
> readfilenames()--> Reads files in the directory and subdirectories in the current path. Reads only the files with ".jl" extension
"""
function readfilenames()
	files=[]
	for file in readdir(normpath(joinpath(@__FILE__,"..","..", "nbs")))
			if endswith(file, ".jl") && !contains(file, "index")
				push!(files,file)
			end
			#if getfile_extension(file)== ".jl"
			#	push!(files,file)
			#end
	end
	files
end

function readfilenames(dir::String)
	files=[]
	for file in readdir(dir)
			if endswith(file, ".jl")
				push!(files,file)
			end
			#if getfile_extension(file)== ".jl"
			#	push!(files,file)
			#end
	end
	files
end
end

#export
"""
> export_file(file::String, path::String, marker::String)--> Loads the file in the supplied path and reads the cells which are marked with "#export". Then saves the notebook in the given path
"""
function export_file(file::String, path::String, marker::String)
	notebook=load_nb(file, "#export")
	save_nb(notebook, path)
end

#export
"""
> export_content(files::AbstractVector, path::String, marker::String)--> maps the `export_file` function to each files
"""
function export_content(files::AbstractVector, path::String, marker::String)
	map(file->export_file(file, path, marker), files)
end

#export
"""
> getfile_extension(filename)--> get the file extensions in the pwd
"""
function getfile_extension(filename)
    return filename[findlast(isequal('.'),filename):end]
end

#export
"""
> notebook2script()--> Export all the code to the ../src path of the project
"""
function notebook2script()
	export_content(readfilenames(), "../src", "#export")
end

end


text:
const DEFAULT_LIN_DISC_FORCE_OUTLIERS_TO_CLOSEST = true

"""
    LinearDiscretizer

A discretizer which encodes (typically) continuous values into discrete bins.
A univariate domain is divided into discrete by edges.
A value V will be encoded into bin B if V ∈ [Bₗ Bᵣ)
(or V ∈ [Bₗ Bᵣ] if B is the rightmost bin)

If force_outliers_to_closest is set to true, then values outside of binedges will be
shunted into the nearest bin. Otherwise an error is thrown.

TODO(tim):
    handle NA
    handle Nullable
    handle NaN
    ability to specify decoding behavior
"""
struct LinearDiscretizer{N<:Real, D<:Integer} <: AbstractDiscretizer{N,D}
    binedges :: Vector{N}   # list of bin edges, sorted smallest to largest
    nbins    :: Int
    i2d      :: Dict{Int,D} # maps bin index to discrete label
    d2i      :: Dict{D,Int} # maps discrete label to bin index

    force_outliers_to_closest :: Bool # if true, real values outside of the bin ranges will be forced to the nearest bin, otherwise will throw an error
end

function LinearDiscretizer( binedges::Vector{N}, i2d::Dict{Int,D};
    force_outliers_to_closest::Bool = DEFAULT_LIN_DISC_FORCE_OUTLIERS_TO_CLOSEST ) where {N<:Real, D<:Integer}

    length(binedges) > 1 || error("bin edges must contain at least 2 values")

    if N <: AbstractFloat # continuous values
        findfirst(i->binedges[i-1] ≥ binedges[i], 2:length(binedges)) == nothing ||
            error("Bin edges must be sorted in increasing order")
    else # for integers, bins of unit width require repeated values
        (findfirst(i->binedges[i-1] ≥ binedges[i], 2:length(binedges)-1) == nothing &&
         findfirst(i->binedges[i-1] > binedges[i], 2:length(binedges)) == nothing ) ||
            error("Bin edges must be sorted in increasing order")
    end

    all(haskey(i2d, i) for i in 1 : length(binedges)-1) || error("i2d must contain all necessary keys")

    d2i = Dict{D,Int}(v => k for (k,v) in i2d)

    return LinearDiscretizer{N,D}(binedges, length(binedges)-1, i2d, d2i, force_outliers_to_closest)
end
function LinearDiscretizer(binedges::Vector{N}, ::Type{D} = Int;
    force_outliers_to_closest::Bool = DEFAULT_LIN_DISC_FORCE_OUTLIERS_TO_CLOSEST ) where {N<:Real, D<:Integer}

    length(binedges) > 1 || error("bin edges must contain at least 2 values")

    if N <: AbstractFloat
        findfirst(i->binedges[i-1] ≥ binedges[i], 2:length(binedges)) == nothing ||
            error("Bin edges must be sorted in increasing order")
    else # for integers, bins of unit width require repeated values
        (findfirst(i->binedges[i-1] ≥ binedges[i], 2:length(binedges)-1) == nothing &&
         findfirst(i->binedges[i-1] > binedges[i], 2:length(binedges)) == nothing ) ||
            error("Bin edges must be sorted in increasing order")
    end

    i2d = Dict{Int,D}(i => convert(D, i) for i in 1 : length(binedges)-1)

    return LinearDiscretizer(binedges,i2d,force_outliers_to_closest=force_outliers_to_closest)
end
function LinearDiscretizer(arr::AbstractVector{N}, ::Type{D}=Int;
    force_outliers_to_closest::Bool = DEFAULT_LIN_DISC_FORCE_OUTLIERS_TO_CLOSEST) where {N<:Real, D<:Integer}
    LinearDiscretizer(convert(Vector{N}, arr), D, force_outliers_to_closest=force_outliers_to_closest)
end

function supports_encoding(ld::LinearDiscretizer{N,D}, x::N) where {N<:Real,D<:Integer}
    if ld.force_outliers_to_closest
        return true
    else
        return ld.binedges[1] ≤ x ≤ ld.binedges[end]
    end
end
supports_decoding(ld::LinearDiscretizer{N,D}, d::D) where {N<:Real,D<:Integer} = 1 ≤ d ≤ ld.nbins

function encode(ld::LinearDiscretizer{N,D}, x::N) where {N<:Real,D<:Integer}
    !isnan(x) || error("cannot encode NaN values")

    if x < ld.binedges[1]
        return ld.force_outliers_to_closest ? ld.i2d[1] : throw(BoundsError())
    elseif x > ld.binedges[end]
        return ld.force_outliers_to_closest ? ld.i2d[ld.nbins] : throw(BoundsError())
    else

        # run bisection search
        binedges = ld.binedges
        a, b = 1, length(binedges)
        va, vb = binedges[a], binedges[b]
        while b-a > 1
            c = div(a+b, 2)
            vc = binedges[c]
            if x < vc
                b, vb = c, vc
            else
                a, va = c, vc
            end
        end
        return ld.i2d[a]
    end
end
encode(ld::LinearDiscretizer{N,D}, x) where {N<:Real,D<:Integer} = encode(ld, convert(N, x))::D
function encode(ld::LinearDiscretizer{N,D}, data::AbstractArray) where {N<:Real,D<:Integer}
    arr = [encode(ld, x) for x in data]
    reshape(arr, size(data))
end

"""
There are several methods for decoding a LinearDiscretizer.
The default is SampleUniform, which uniformly samples from the bin.
SampleBinCenter returns the bin's center value.
"""
function decode(ld::LinearDiscretizer{N,D}, d::D, ::SampleUniform) where {N<:AbstractFloat,D<:Integer}
    ind = ld.d2i[d]
    lo  = ld.binedges[ind]
    hi  = ld.binedges[ind+1]
    convert(N, lo + rand()*(hi-lo))
end
function decode(ld::LinearDiscretizer{N,D}, d::D, ::SampleBinCenter) where {N<:AbstractFloat,D<:Integer}
    ind = ld.d2i[d]
    lo  = ld.binedges[ind]
    hi  = ld.binedges[ind+1]
    convert(N, (hi + lo)/2)
end
decode(ld::LinearDiscretizer{N,D}, d::D) where {N<:AbstractFloat,D<:Integer} = decode(ld, d, SAMPLE_UNIFORM)

function decode(ld::LinearDiscretizer{N,D}, d::D, ::SampleUniform) where {N<:Integer,D<:Integer}
    ind = ld.d2i[d]
    lo  = ld.binedges[ind]
    hi  = ld.binedges[ind+1]
    if hi != ld.binedges[end]
        retval = rand(lo:hi-1)::Int
    else
        retval = rand(lo:hi)::Int
    end
    convert(N, retval)
end
function decode(ld::LinearDiscretizer{N,D}, d::D, ::SampleBinCenter) where {N<:Integer,D<:Integer}
    ind = ld.d2i[d]
    lo = ld.binedges[ind]
    hi = ld.binedges[ind+1]
    convert(N, div(lo+hi,2))
end
decode(ld::LinearDiscretizer{N,D}, d::D) where {N<:Integer,D<:Integer} = decode(ld, d, SAMPLE_UNIFORM)

decode(ld::LinearDiscretizer{N,D}, d::I, method::AbstractSampleMethod=SAMPLE_UNIFORM) where {N<:Real,D<:Integer,I<:Integer} =
    decode(ld, convert(D,d), method)

function decode(ld::LinearDiscretizer{N,D}, data::AbstractArray{D}, ::AbstractSampleMethod=SAMPLE_UNIFORM) where {N<:Real,D<:Integer}
    arr = Vector{N}(undef, length(data))
    for (i,d) in enumerate(data)
        arr[i] = decode(ld, d)
    end
    reshape(arr, size(data))
end

Base.max(ld::LinearDiscretizer) = ld.binedges[end]
Base.min(ld::LinearDiscretizer) = ld.binedges[1]
function Base.extrema(ld::LinearDiscretizer{N,D}) where {N<:Real,D<:Integer}
    lo  = ld.binedges[1]
    hi  = ld.binedges[end]
    return (lo, hi)
end
function Base.extrema(ld::LinearDiscretizer{N,D}, d::D) where {N<:AbstractFloat,D<:Integer}
    ind = ld.d2i[d]
    lo  = ld.binedges[ind]
    hi  = ld.binedges[ind+1]
    return (lo, hi)
end
function Base.extrema(ld::LinearDiscretizer{N,D}, d::D) where {N<:Integer,D<:Integer}
    ind = ld.d2i[d]
    lo  = ld.binedges[ind]
    hi  = ld.binedges[ind+1]
    if hi == ld.binedges[end]
        (lo, hi)
    else
        (lo, hi-1)
    end
end
function totalwidth(ld::LinearDiscretizer)
    ex = extrema(ld)
    ex[2] - ex[1]
end

nlabels(ld::LinearDiscretizer) = ld.nbins
binedges(ld::LinearDiscretizer) = ld.binedges
bincenters(ld::LinearDiscretizer{N,D}) where {N<:AbstractFloat,D<:Integer} = (0.5*(ld.binedges[1:ld.nbins] + ld.binedges[2:end]))::Vector{Float64}
function bincenters(ld::LinearDiscretizer{N,D}) where {N<:Integer,D<:Integer}
    retval = Vector{Float64}(undef, ld.nbins)
    for i = 1 : length(retval)-1
        retval[i] = 0.5(ld.binedges[i+1]-1 + ld.binedges[i])
    end
    retval[end] = 0.5(ld.binedges[end] + ld.binedges[end-1])
    retval
end
binwidth(ld::LinearDiscretizer{N,D}, d::D) where {N<:Real,D<:Integer} = ld.binedges[d+1] - ld.binedges[d]
binwidths(ld::LinearDiscretizer{N,D}) where {N<:Real,D<:Integer} = ld.binedges[2:end] - ld.binedges[1:end-1]

ext: jl


text:
# 16 opcodes.  don't know which one is which.  instructions are of the form
# op A B C
#
# they all write their output to register C.
#
# differ based on
# * operation they perform (add, mul, and, or, assign, greater than, equals)
# * source of input1 (register vs. identity)
# * source of input2 (register vs. identity)

struct Op
    f::Function
    areg::Bool                  # register A or value A
    breg::Bool                  # register B or value B
end

ops = (addr = Op(+, true, true),
       addi = Op(+, true, false),
       mulr = Op(*, true, true),
       muli = Op(*, true, false),
       banr = Op(&, true, true),
       bani = Op(&, true, false),
       borr = Op(|, true, true),
       bori = Op(|, true, false),
       setr = Op((a,b) -> a, true, false),
       seti = Op((a,b) -> a, false, false),
       eqir = Op(isequal, false, true),
       eqri = Op(isequal, true, false),
       eqrr = Op(isequal, true, true),
       gtir = Op((a,b) -> a>b, false, true),
       gtri = Op((a,b) -> a>b, true, false),
       gtrr = Op((a,b) -> a>b, true, true))

function call_op(o::Op, registers, instruction)
    opcode, a, b, c = instruction
    a = o.areg ? registers[a+1] : a
    b = o.breg ? registers[b+1] : b
    return o.f(a,b)
end

function check(o::Op, before, instruction, after)
    c = last(instruction)
    call_op(o, before, instruction) == after[c+1]
end

function count_ops(before, instruction, after)
    matches = 0
    for (name,op) in pairs(ops)
        if check(op, before, instruction, after)
            matches += 1
            println("$name matches")
        end
    end
    return matches
end

input1, input2 = split(read("16.input", String), "\n\n\n");

tupleify(v::Vector) = tuple(v...)
parse_line(s) = tupleify(parse.(Int, match(r"(\d+).*(\d+).*(\d+).*(\d+)", s).captures))

function star1(input1)

    three_ops = 0

    for input in split(input1, "\n\n")
        before, instruction, after = parse_line.(split(input, "\n"))
        if count_ops(before, instruction, after) >= 3
            three_ops += 1
        end
    end

    return three_ops
end

star1(input1)

function op_matches(before, instruction, after)
    matches = Set{Op}()
    for op in ops
        if check(op, before, instruction, after)
            push!(matches, op)
        end
    end
    return matches
end

function star2(input1, input2)
    op_candidates = Dict{Int,Set{Op}}()

    for input in split(input1, "\n\n")
        before, instruction, after = parse_line.(split(input, "\n"))
        matches = op_matches(before, instruction, after)
        intersect!(get!(op_candidates, instruction[1], matches), matches)
    end

    # this isn't enough to uniquely identify.

    opcodes = missings(Op, 16)
    while any(ismissing, opcodes)
        code = findfirst(c->length(c)==1, op_candidates)
        op = first(op_candidates[code])
        opcodes[code+1] = op
        delete!(op_candidates, code)
        for cands in values(op_candidates)
            delete!(cands, op)
        end
    end

    opnames = findfirst.(isequal.(opcodes), Ref(ops))

    registers = [0, 0, 0, 0]
    for line in split(strip(input2), "\n")
        println("Before: $registers")
        inst = parse.(Int, split(line))
        println("$inst ($(opnames[inst[1]+1]))")
        registers[inst[4]+1] = call_op(opcodes[inst[1]+1], registers, inst)
        println("After: $registers\n")
    end

    return op_candidates
end

op_cands = star2(input1, input2)


text:
using Flux, DiffEqFlux, StochasticDiffEq, Plots, DiffEqMonteCarlo

u0 = Float32[2.; 0.]
datasize = 30
tspan = (0.0f0,1.0f0)

function trueODEfunc(du,u,p,t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end
t = range(tspan[1],tspan[2],length=datasize)
mp = Float32[0.2,0.2]
function true_noise_func(du,u,p,t)
    du .= mp.*u
end
prob = SDEProblem(trueODEfunc,true_noise_func,u0,tspan)

# Take a typical sample from the mean
monte_prob = MonteCarloProblem(prob)
monte_sol = solve(monte_prob,SOSRI(),num_monte = 100)
monte_sum = MonteCarloSummary(monte_sol)
sde_data = Array(timeseries_point_mean(monte_sol,t))

dudt = Chain(x -> x.^3,
             Dense(2,50,tanh),
             Dense(50,2))
ps = Flux.params(dudt)
n_sde = x->neural_dmsde(dudt,x,mp,tspan,SOSRI(),saveat=t,reltol=1e-1,abstol=1e-1)

pred = n_sde(u0) # Get the prediction using the correct initial condition

dudt_(u,p,t) = Flux.data(dudt(u))
g(u,p,t) = mp.*u
nprob = SDEProblem(dudt_,g,u0,(0.0f0,1.2f0),nothing)

monte_nprob = MonteCarloProblem(nprob)
monte_nsol = solve(monte_nprob,SOSRI(),num_monte = 100)
monte_nsum = MonteCarloSummary(monte_nsol)
#plot(monte_nsol,color=1,alpha=0.3)
p1 = plot(monte_nsum, title = "Neural SDE: Before Training")
scatter!(p1,t,sde_data',lw=3)

scatter(t,sde_data[1,:],label="data")
scatter!(t,Flux.data(pred[1,:]),label="prediction")

function predict_n_sde()
  n_sde(u0)
end
loss_n_sde1() = sum(abs2,sde_data .- predict_n_sde())
loss_n_sde10() = sum([sum(abs2,sde_data .- predict_n_sde()) for i in 1:10])
Flux.back!(loss_n_sde1())

data = Iterators.repeated((), 10)
opt = ADAM(0.025)
cb = function () #callback function to observe training
  sample = predict_n_sde()
  # loss against current data
  display(sum(abs2,sde_data .- sample))
  # plot current prediction against data
  cur_pred = Flux.data(sample)
  pl = scatter(t,sde_data[1,:],label="data")
  scatter!(pl,t,cur_pred[1,:],label="prediction")
  display(plot(pl))
end

# Display the SDE with the initial parameter values.
cb()

Flux.train!(loss_n_sde1 , ps, Iterators.repeated((), 100), opt, cb = cb)
Flux.train!(loss_n_sde10, ps, Iterators.repeated((), 20), opt, cb = cb)

dudt_(u,p,t) = Flux.data(dudt(u))
g(u,p,t) = mp.*u
nprob = SDEProblem(dudt_,g,u0,(0.0f0,1.2f0),nothing)

monte_nprob = MonteCarloProblem(nprob)
monte_nsol = solve(monte_nprob,SOSRI(),num_monte = 100)
monte_nsum = MonteCarloSummary(monte_nsol)
#plot(monte_nsol,color=1,alpha=0.3)
p2 = plot(monte_nsum, title = "Neural SDE: After Training", xlabel="Time")
scatter!(p2,t,sde_data',lw=3,label=["x" "y" "z" "y"])

plot(p1,p2,layout=(2,1))

savefig("neural_sde.pdf")
savefig("neural_sde.png")

max_stars_repo_name: MikeInnes/zygote-paper

max_issues_repo_name: MikeInnes/zygote-paper

max_forks_repo_name: MikeInnes/zygote-paper


text:
const conemap = Dict("L=" => :Zero, "F" => :Free,
                     "L-" => :NonPos, "L+" => :NonNeg,
                     "Q" => :SOC, "QR" => :SOCRotated,
                     "EXP" => :ExpPrimal, "EXP*" => :ExpDual)
const conemap_rev = Dict(:Zero => "L=", :Free => "F",
                     :NonPos => "L-", :NonNeg => "L+",
                     :SOC => "Q", :SOCRotated => "QR",
                     :ExpPrimal => "EXP", :ExpDual => "EXP*")

function cbfcones_to_mpbcones(c::Vector{Tuple{String,Int}},total)
    i = 1
    mpb_cones = Vector{Tuple{Symbol,Vector{Int}}}(0)

    for (cname,count) in c
        conesymbol = conemap[cname]
        if conesymbol == :ExpPrimal
            @assert count == 3
            indices = i+2:-1:i
        else
            indices = i:(i+count-1)
        end
        push!(mpb_cones, (conesymbol, collect(indices)))
        i += count
    end
    @assert i == total + 1
    return mpb_cones
end

# https://github.com/JuliaLang/julia/issues/13942#issuecomment-217324812
function unzip{T<:Tuple}(A::Array{T})
    res = map(x -> x[], T.parameters)
    res_len = length(res)
    for t in A
        for i in 1:res_len
            push!(res[i], t[i])
        end
    end
    return res
end

psd_len(n) = div(n*(n+1),2)
# returns offset from starting index for (i,j) term in n x n matrix
function idx_to_offset(n,i,j)
    @assert 1 <= i <= n
    @assert 1 <= j <= n
    # upper triangle
    if i > j
        i,j = j,i
    end
    # think row major
    return psd_len(n) - psd_len(n-i+1) + (j-i)
end

function cbftompb(dat::CBFData)
    @assert dat.nvar == (isempty(dat.var) ? 0 : sum(c->c[2],dat.var))
    @assert dat.nconstr == (isempty(dat.con) ? 0 : sum(c->c[2],dat.con))

    c = zeros(dat.nvar)
    for (i,v) in dat.objacoord
        c[i] = v
    end

    var_cones = cbfcones_to_mpbcones(dat.var, dat.nvar)
    con_cones = cbfcones_to_mpbcones(dat.con, dat.nconstr)

    I_A, J_A, V_A = unzip(dat.acoord)
    b = zeros(dat.nconstr)
    for (i,v) in dat.bcoord
        b[i] = v
    end

    psdvarstartidx = Int[]
    for i in 1:length(dat.psdvar)
        if i == 1
            push!(psdvarstartidx,dat.nvar+1)
        else
            push!(psdvarstartidx,psdvarstartidx[i-1] + psd_len(dat.psdvar[i-1]))
        end
        push!(var_cones,(:SDP,psdvarstartidx[i]:psdvarstartidx[i]+psd_len(dat.psdvar[i])-1))
    end
    nvar = (length(dat.psdvar) > 0) ? psdvarstartidx[end] + psd_len(dat.psdvar[end]) - 1 : dat.nvar

    psdconstartidx = Int[]
    for i in 1:length(dat.psdcon)
        if i == 1
            push!(psdconstartidx,dat.nconstr+1)
        else
            push!(psdconstartidx,psdconstartidx[i-1] + psd_len(dat.psdcon[i-1]))
        end
        push!(con_cones,(:SDP,psdconstartidx[i]:psdconstartidx[i]+psd_len(dat.psdcon[i])-1))
    end
    nconstr = (length(dat.psdcon) > 0) ? psdconstartidx[end] + psd_len(dat.psdcon[end]) - 1 : dat.nconstr

    c = [c;zeros(nvar-dat.nvar)]
    for (matidx,i,j,v) in dat.objfcoord
        ix = psdvarstartidx[matidx] + idx_to_offset(dat.psdvar[matidx],i,j)
        @assert c[ix] == 0.0
        scale = (i == j) ? 1.0 : sqrt(2)
        c[ix] = scale*v
    end

    for (conidx,matidx,i,j,v) in dat.fcoord
        ix = psdvarstartidx[matidx] + idx_to_offset(dat.psdvar[matidx],i,j)
        push!(I_A,conidx)
        push!(J_A,ix)
        scale = (i == j) ? 1.0 : sqrt(2)
        push!(V_A,scale*v)
    end

    for (conidx,varidx,i,j,v) in dat.hcoord
        ix = psdconstartidx[conidx] + idx_to_offset(dat.psdcon[conidx],i,j)
        push!(I_A,ix)
        push!(J_A,varidx)
        scale = (i == j) ? 1.0 : sqrt(2)
        push!(V_A,scale*v)
    end

    b = [b;zeros(nconstr-dat.nconstr)]
    for (conidx,i,j,v) in dat.dcoord
        ix = psdconstartidx[conidx] + idx_to_offset(dat.psdcon[conidx],i,j)
        @assert b[ix] == 0.0
        scale = (i == j) ? 1.0 : sqrt(2)
        b[ix] = scale*v
    end

    A = sparse(I_A,J_A,-V_A,nconstr,nvar)

    vartypes = fill(:Cont, nvar)
    vartypes[dat.intlist] = :Int

    return c, A, b, con_cones, var_cones, vartypes, dat.sense, dat.objoffset
end

function mpbtocbf(name, c, A, b, con_cones, var_cones, vartypes, sense=:Min)
    num_scalar_var = 0
    for (cone,idx) in var_cones
        if cone != :SDP
            num_scalar_var += length(idx)
        end
    end
    num_scalar_con = 0
    for (cone,idx) in con_cones
        if cone != :SDP
            num_scalar_con += length(idx)
        end
    end

    # need to shuffle rows and columns to put them in order
    var_idx_old_to_new = zeros(Int,length(c))
    con_idx_old_to_new = zeros(Int,length(b))
    var_idx_new_to_old = zeros(Int,num_scalar_var)
    con_idx_new_to_old = zeros(Int,num_scalar_con)

    # CBF fields
    var = Vector{Tuple{String,Int}}(0)
    con = Vector{Tuple{String,Int}}(0)

    i = 1
    for (cone,idx) in var_cones
        if cone == :ExpPrimal
            @assert all(var_idx_old_to_new[idx] .== 0)
            @assert length(idx) == 3
            # MPB: (x,y,z) : y*exp(x/y) <= z
            # CBF: (z,y,x) : y*exp(x/y) <= z
            var_idx_old_to_new[idx] = i+2:-1:i
            var_idx_new_to_old[i+2:-1:i] = idx
            i += 3
            push!(var, (conemap_rev[cone],length(idx)))
        elseif cone != :SDP
            for k in idx
                var_idx_old_to_new[k] = i
                var_idx_new_to_old[i] = k
                i += 1
            end
            push!(var, (conemap_rev[cone],length(idx)))
        end
    end
    @assert i - 1 == num_scalar_var

    i = 1
    for (cone,idx) in con_cones
        if cone == :ExpPrimal
            @assert all(con_idx_old_to_new[idx] .== 0)
            @assert length(idx) == 3
            con_idx_old_to_new[idx] = i+2:-1:i
            con_idx_new_to_old[i+2:-1:i] = idx
            i += 3
            push!(con, (conemap_rev[cone],length(idx)))
        elseif cone != :SDP
            for k in idx
                @assert con_idx_old_to_new[k] == 0
                con_idx_old_to_new[k] = i
                con_idx_new_to_old[i] = k
                i += 1
            end
            push!(con, (conemap_rev[cone],length(idx)))
        end
    end
    @assert i - 1 == num_scalar_con

    objacoord = collect(zip(findnz(sparse(c[var_idx_new_to_old]))...))
    bcoord = collect(zip(findnz(sparse(b[con_idx_new_to_old]))...))
    # MPB is b - Ax ∈ K, CBF is b + Ax ∈ K
    Acbf = -A[con_idx_new_to_old,var_idx_new_to_old]

    acoord = collect(zip(findnz(Acbf)...))::Vector{Tuple{Int,Int,Float64}}

    intlist = Int[]
    for i in 1:length(vartypes)
        if var_idx_old_to_new[i] == 0 && vartypes[i] != :Cont
            error("CBF format does not support integer restrictions on PSD variables")
        end
        if vartypes[i] == :Cont
        elseif vartypes[i] == :Int
            push!(intlist,var_idx_old_to_new[i])
        elseif vartypes[i] == :Bin
            # TODO: Check if we need to add variable bounds also
            push!(intlist,var_idx_old_to_new[i])
        else
            error("Unrecognized variable category $(vartypes[i])")
        end
    end

    psdvar = Int[]
    psdcon = Int[]

    # Map from MPB linear variable index to (psdvar,i,j)
    psdvar_idx_old_to_new = fill((0,0,0), length(c))
    # Map from MPB linear constraint index to (psdcon,i,j)
    psdcon_idx_old_to_new = fill((0,0,0),length(b))

    for (cone,idx) in var_cones
        if cone == :SDP
            y = length(idx)
            conedim = round(Int, sqrt(0.25 + 2y) - 0.5)
            push!(psdvar, conedim)
            k = 1
            for i in 1:conedim, j in i:conedim
                psdvar_idx_old_to_new[idx[k]] = (length(psdvar), i, j)
                k += 1
            end
            @assert length(idx) == k - 1
        end
    end

    for (cone,idx) in con_cones
        if cone == :SDP
            y = length(idx)
            conedim = round(Int, sqrt(0.25 + 2y) - 0.5)
            push!(psdcon, conedim)
            k = 1
            for i in 1:conedim, j in i:conedim
                psdcon_idx_old_to_new[idx[k]] = (length(psdcon), i, j)
                k += 1
            end
            @assert length(idx) == k - 1
        end
    end

    objfcoord = Vector{Tuple{Int,Int,Int,Float64}}(0)
    for i in 1:length(c)
        if c[i] != 0.0 && psdvar_idx_old_to_new[i] != (0,0,0)
            varidx, vari, varj = psdvar_idx_old_to_new[i]
            scale = (vari == varj) ? 1.0 : sqrt(2)
            push!(objfcoord, (varidx, vari, varj, c[i]/scale))
        end
    end
    dcoord = Vector{Tuple{Int,Int,Int,Float64}}(0)
    for i in 1:length(b)
        if b[i] != 0.0 && psdcon_idx_old_to_new[i] != (0,0,0)
            conidx, coni, conj = psdcon_idx_old_to_new[i]
            scale = (coni == conj) ? 1.0 : sqrt(2)
            push!(dcoord, (conidx, coni, conj, b[i]/scale))
        end
    end

    A_I, A_J, A_V = findnz(A)
    fcoord = Vector{Tuple{Int,Int,Int,Int,Float64}}(0)
    hcoord = Vector{Tuple{Int,Int,Int,Int,Float64}}(0)

    for (i,j,v) in zip(A_I,A_J,A_V)
        if psdvar_idx_old_to_new[j] != (0,0,0)
            if psdcon_idx_old_to_new[i] != (0,0,0)
                error("CBF format does not allow PSD variables to appear in affine expressions defining PSD constraints")
            end
            newrow = con_idx_old_to_new[i]
            @assert newrow != 0
            varidx, vari, varj = psdvar_idx_old_to_new[j]
            scale = (vari == varj) ? 1.0 : sqrt(2)
            push!(fcoord, (newrow, varidx, vari, varj, -v/scale))
        elseif psdcon_idx_old_to_new[i] != (0,0,0)
            newcol = var_idx_old_to_new[j]
            conidx, coni, conj = psdcon_idx_old_to_new[i]
            scale = (coni == conj) ? 1.0 : sqrt(2)
            push!(hcoord, (conidx, newcol, coni, conj, -v/scale))
        end
    end

    return CBFData(name,sense,var,psdvar,con,psdcon,objacoord,objfcoord,0.0,fcoord,acoord,bcoord,hcoord,dcoord,intlist,num_scalar_var,num_scalar_con)
end

# converts an MPB solution to CBF solution
# no transformation needed unless PSD vars present
function mpb_sol_to_cbf(dat::CBFData,x::Vector)
    scalar_solution = x[1:dat.nvar]

    psdvar_solutions = Vector{Matrix{Float64}}(0)
    startidx = dat.nvar+1
    for i in 1:length(dat.psdvar)
        endidx = startidx + psd_len(dat.psdvar[i]) - 1
        svec_solution = x[startidx:endidx]
        push!(psdvar_solutions, make_smat!(Matrix{Float64}(dat.psdvar[i],dat.psdvar[i]), svec_solution))
        startidx = endidx + 1
    end

    return scalar_solution, psdvar_solutions
end

# Copied from Pajarito.jl
function make_smat!(smat::Matrix{Float64}, svec::Vector{Float64})
    dim = size(smat, 1)
    kSD = 1
    for jSD in 1:dim, iSD in jSD:dim
        if jSD == iSD
            smat[iSD, jSD] = svec[kSD]
        else
            smat[iSD, jSD] = smat[jSD, iSD] = (1/sqrt(2)) * svec[kSD]
        end
        kSD += 1
    end
    return smat
end


text:
using Unitful
using UnitfulAtomic
using BesselFunctions
using DifferentialEquations
using LegendrePolynomials
using QuadGK
using Dierckx

const a0 = u"a0_au"
const Eh = u"Eh_au"
const me = u"me_au"
const ħ = u"ħ_au"

kFromE(E,m=1me) = uconvert(a0^-1, sqrt(2*m*E) / ħ)
EFromk(k,m=1me) = uconvert(Eh, ħ^2*k^2/2m )

"""
    DerivFunc(u,u′,r, E,k,potential,l,m=1mₑ)

This is the differential equation derivative function. This function should work
in Unitful quantities.
"""
function DerivFunc(u,u′,r, E,k,potential,l,m=1me)
    # Expressing this as the plane wave expansion with the Y_l0 at the end,
    # gives rise to the sqrt() prefactors.
    jl = sqrt(4pi)*sqrt(2l+1) * 1/k * ricj(l, k*r)

    prefac = 2m/ħ^2
    U = -prefac*potential(r)

    u′′ = @. (l*(l+1)/r^2)*u - ((prefac*E + U)*u + U*jl)

    return u′,u′′
end

"""
    SolveOutwards(rspan, E,k,potential,l ; init_derivs=[0.0,1.0])

Calculates the wavefunction for the off-shell T-matrix problem. Integrates
outwards using `DifferentialEquations` with `rspan`. The ODE derivative is
created from the `potential` given (which should handle `Unitful` quantities).

The incoming wave has a wavenumber of `k` and the energy that the problem is
calculated on is `E`.

TODO: Put the ODE here.

Note: the inhomogeneous part is not included in the solution. i.e. this is
solving for \tilde{f}. Also not that there are two solutions calculated, which
need to be matched to the appropriate boundary conditions using e.g. `ElminateIngoing`.

Returns (r,u,u′)
"""
function SolveOutwards(rspan, E,k,potential,l ; init_derivs=[0.0,1.0])
    # Going to deal with units in the deriv funcs
    # Note: wavefunctions are unitless in the picture with \phi_k = (2\pi)^-3/2 exp(ikr)
    # But u = r*R so there is units for u.

    prob_func = function (f,p,r)
        r = r*a0
        u = f[1,:] .* a0
        u′ = f[2,:]
        du,du′ = DerivFunc(u, u′,r, E, k, potential,l)
        [ustrip.(NoUnits, du)' ; ustrip.(1/a0, du′)']
    end

    # Not sure if I will need to make this a more carefully chosen derivative.
    # At least the value being zero will be necessary.
    u0 = ComplexF64[[0. 0.] ; permutedims(init_derivs)]

    # Needs to include the jl part that is in f, as the BC is for f not tildef.
    jl = sqrt(4pi)*sqrt(2l+1) * 1/k * ricj(l, k*rspan[1])
    u0 .-= [ustrip(a0, jl) ; 0.]

    prob = ODEProblem(prob_func, u0, ustrip.(a0, rspan))
    sol = solve(prob, Tsit5(), abstol=1e-10, reltol=1e-12)

    this_soln = reduce((x,y) -> cat(x,y,dims=3), sol.u)
    r = sol.t .* a0

    u = this_soln[1,:,:] .* a0
    u′ = this_soln[2,:,:]
    return r, u, u′
end

"""
Old function, hasn't been tested since included in the module.
"""
function SolveInwards(r_list, E,k,potential,l, A)

    prob_func = function (f,p,r)
        r = r*a0
        u = f[1] .* a0
        u′ = f[2]
        du,du′ = DerivFunc(u, u′,r, E, k, potential,l)
        [ustrip.(NoUnits, du) ; ustrip.(1/a0, du′)]
    end

    rspan = extrema(r_list)
    rmax = rspan[2]

    κ = kFromE(E)
    u0 = A*[rich₊(l, κ*rmax) ; ustrip(1/a0, drich₊dr(l, κ, rmax))]

    prob = ODEProblem(prob_func, u0, ustrip.(a0, reverse(rspan)))
    sol = solve(prob, Tsit5(), abstol=1e-10, reltol=1e-12)#, callback=DiscreteCallback(RenormCond, RenormAffect!))

    # this_r = filter(x -> rspan[1] <= x <= sol.t[end]*a0, r_list)
    this_r = r_list
    this_soln = sol.(ustrip.(a0,this_r))
    this_soln = reduce(hcat, this_soln)

    r = this_r
    u = this_soln[1,:] .* a0
    u′ = this_soln[2,:]
    # coeffs = MatchSolutionUnitless(u, u′, sol.t[end], ustrip(1/a0,k), l)

    # conv = TransformToPureIngoingOutgoing(coeffs)

    return r, u, u′
end

function MatchSolution(u, u′, r, k, l)
    u = ustrip.(a0,u)
    u′ = u′
    r = ustrip(a0,r)
    k = ustrip(1/a0,k)
    MatchSolutionUnitless(u,u′,r,k,l)
end
function MatchSolution(f, r, k, l)
    r = ustrip(a0,r)
    k = ustrip(1/a0,k)
    MatchSolutionUnitless(f, r, k, l)
end
function MatchSolutionUnitless(u, u′, r, k, l)
    f = [u' ; u′']
    MatchSolutionUnitless(f, r, k, l)
end
function MatchSolutionUnitless(f, r, κ, l)
    M = [rich₊(l,κ*r)    rich₋(l,κ*r) ;
         drich₊dr(l,κ,r) drich₋dr(l,κ,r)]
    coeffs = M\f
end

"""
TransformToPureIngoingOutgoing(C)

Given a coeff matrix `C`, relating two solutions obtained from `SolveOutwards`
to the ingoing/outgoing states, combine the two solutions to satisfy the
boundary conditions.

Note: because of the inhomogeneity, this must satisfy u = A·u1 + (1-A)·u2
instead of the usual freedom of u = A·u1 + B·u2.
"""
function TransformToPureIngoingOutgoing(C)
    x1 = -C[2,2] / (C[2,1] - C[2,2])
    x2 = -C[1,2] / (C[1,1] - C[1,2])

    M = [[x1;(1-x1)] [x2;(1-x2)]]
    return M
end

"""
EliminateIngoing(u, u′, r, k, l)

Matches to the boundary conditions for \tilde{u} using the output from `SolveOutwards`.
"""
function EliminateIngoing(u, u′, r, k, l)
    coeffs = MatchSolution(u[:,end], u′[:,end], r[end], k, l)
    M = TransformToPureIngoingOutgoing(coeffs)
    M = M[:,1:1]
    M = permutedims(M)

    return dropdims(M*u, dims=1),
           dropdims(M*u′, dims=1)
end

"""
    AddBesselPart(r, tildeu, tildeu′, k, l)

Includes the inhomogeneity that was missing from the ODE solution. i.e. converts
from \tilde{u} to u.
"""
function AddBesselPart(r, tildeu, tildeu′, k, l)
    u = sqrt(4pi)*sqrt(2l+1) * 1/k * ricj.(l,k*r) + tildeu
    u′ = sqrt(4pi)*sqrt(2l+1) * 1/k * dricjdr.(l,k,r) + tildeu′

    return u,u′
end

"""
Tmat_l(r,u,potential,p,costh,l)

Once a complete solution for u has been obtained, this function calculates the
`T` matrix element. The `u` going into here must be the full matched solution
and not \tilde{u}.

Note: this was before I considered T_l. This term is instead one of the Legendre
contributions of T(theta) for a specific `costh` and `l`.
"""
function Tmat_l(r,u,potential,p,costh,l)
    # One l component only. costh = p⋅z (or p⋅k more generally)
    V = r -> ustrip(Eh,potential(r*a0))
    p2 = ustrip(1/a0,p)
    r2 = ustrip.(a0,r)

    interp_r = Spline1D(r2,real.(ustrip.(a0,u)),k=1)
    interp_i = Spline1D(r2,imag.(ustrip.(a0,u)),k=1)
    interp = r -> interp_r(r) + im*interp_i(r)
    out, = quadgk(r -> interp(r)*V(r)*ricj(l,p2*r), extrema(r2)...)
    out *= a0*Eh*a0
    # These terms come from the plane wave prefactors for p.
    out *= sqrt(4pi)*sqrt(2l+1) * 1/p * Pl(costh,l)
    # Adjusting for my lack of (2π)^-3/2 in the wavefunctions
    out /= (2π)^3
end

fFromT(T,m=me) = uconvert(a0, -(2pi)^2*m/ħ^2 * T)

"""
TCSandOptTheorem(faa, k)

Compare the total cross section with the optical theorem result, assuming that
`faa` is the forward scattering part (i.e. `faa` = f(θ=0)) and assuming a low
enough energy that the scattering is s-wave.

`faa` can be obtained from `fFromT`.
"""
function TCSandOptTheorem(faa, k)
    DCS = abs2(faa)
    TCS = 4pi * DCS

    opt_theorem = 4pi/k * imag(faa)

    a = sqrt(TCS/4pi)
    return (DCS=DCS, TCS=TCS, opt=opt_theorem, a=a)
end

# |faa|^2 =a^2
# -(2pi)^2*m/\hbar^2 * T = a
# hbr^2 a / (2pi)^2 / m


text:
#=
This file is part of the replication code for: Hasenzagl, T., Pellegrino, F., Reichlin, L., & Ricco, G. (2020). A Model of the Fed's View on Inflation.
Please cite the paper if you are using any part of the code for academic work (including, but not limited to, conference and peer-reviewed papers).
=#

import Pkg;

# Packages used in the estimation
Pkg.add("DataFrames")
Pkg.add("Distributions")
Pkg.add("FileIO")
Pkg.add("JLD")
Pkg.add("XLSX")

# Packages used to make the charts
Pkg.add("Colors")
Pkg.add("CSV")
Pkg.add("ORCA")
Pkg.add("KernelDensity")
Pkg.add("PlotlyJS")


text:
export CVOTDR

mutable struct CVOTDR # Confirmed, vaccinated, tested, deceased, recovered
    confirmed::Int
    vaccinated::Int
    other::Int
    tested::Int
    deceased::Int
    recovered::Int

    CVOTDR(confirmed, vaccinated, other, tested, deceased, recovered) = new(confirmed, vaccinated, other, tested, deceased, recovered)
end

function CVOTDR(;confirmed=0, vaccinated=0, other=0, tested=0, deceased=0, recovered=0)
    return CVOTDR(confirmed, vaccinated, other, tested, deceased, recovered)
end

function CVOTDR(datadict::Dict{String,Int})
    kys = keys(datadict)
    cvotdr = CVOTDR()
    if "confirmed" in kys
        cvotdr.confirmed = datadict["confirmed"]
    end
    if "other" in kys
        cvotdr.other = datadict["other"]
    end
    if "vaccinated" in kys
        cvotdr.vaccinated = datadict["vaccinated"]
    end
    if "tested" in kys
        cvotdr.tested = datadict["tested"]
    end
    if "deceased" in kys
        cvotdr.deceased = datadict["deceased"]
    end
    if "recovered" in kys
        cvotdr.recovered = datadict["recovered"]
    end
    return cvotdr
end

function get_timeseries_statewise(state::String, date::Dates.Date)
    url = timeseries_url
    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj[state_dict[state]]["dates"][string(date)]
        df = DataFrame(date=Dates.Date[], delta7=CVOTDR[], delta=CVOTDR[], total=CVOTDR[])
        push!(df.date, date)
        d7data = datadict["delta7"] # delta 7 data
        tdata = datadict["total"] # total data
        ddata = datadict["delta"] # delta

        d7 = CVOTDR(Dict{String,Int}(d7data))
        d = CVOTDR(Dict{String,Int}(ddata))
        t = CVOTDR(Dict{String,Int}(tdata))

        push!(df.delta7, d7)
        push!(df.delta, d)
        push!(df.total, t)

        return df
    catch e
        return e
    end
end

function get_districts(state::String)
    url = data_today_url
    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        state = state_dict[state]
        datadict = jobj[state]["districts"]
        return keys(datadict)
    catch e
        return e
    end
end

function get_district_current_data(state::String, district::String)
    url = data_today_url
    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        state = state_dict[state]
        datadict = jobj[state]["districts"][district]
        return datadict
    catch e
        return e
    end
end


text:
"""
   run_subnets(
      parent_folder::String,
      objective::String;
      gen_areas::Array=[],
      area_transfer::Array=[],
      gen_zones::Array=[],
      zone_transfer::Array=[],
      enum_branches::Bool=false,
      plot_best_x::Int64=-1,
      scopf::Bool=false,
      ctg_plots::Array{Int64,1}=[0],
      run_fix_f::Bool=false,
      run_indir_PQ::Bool=false,
      print_results::Bool=false
   )

Models and solves an OPF for every network in a directory.

# Arguments
- `parent_folder::String`: a directory containing full network data for one or more networks, each in a folder containing all subnetwork data, *subnetworks.csv*, and *interfaces.csv*
- `obj::String`: the objective function to use, from the following:
   - "mincost": minimize generation cost
   - "areagen": minimize generation in the areas specified in `gen_areas`
   - "zonegen": minimize generation in the zones specified in `gen_zones`
   - "minredispatch": minimize the change in generator dispatch from the initial values defined in the network data
- `gen_areas::Array{Int64,1}`: all areas in which generation should be minimized if `obj=="areagen"`
- `area_transfer::Array{Int64,1}`: two areas with power transfer between them that should be saved and plotted. Results for P, Q, S, and loss are saved for power transfer between the two areas. Must have exactly two elements.
- `gen_zones::Array{Int64,1}`: all zones in which generation should be minimized if `obj=="zonegen"`
- `zone_transfer::Array{Int64,1}`: two zones with power transfer between them that should be saved and plotted. Results for P, Q, S, and loss are saved for power transfer between the two zones. Must have exactly two elements.
- `enum_branches::Bool`: if true, collect results from each folder for plotting bar graphs. This is used when the possible branch upgrades have been enumerated and a comparison is desired. Default false.
- `plot_best_x::Int64`: number of results to plot, sorted from smallest to largest objective. If `plot_best_x` <= 1, the results of all networks which gave feasible solutions are plotted. Default -1.
- `scopf::Bool`: if true, model and solve the N-1 security constrained OPF for each network. Each network folder must contain a contingency specification file (_*.con_) for each subnetwork. Default false.
- `ctg_plots::Array{Int64,1}`: indices of the contingencies to plot. The base case index is 0. Default [0].
"""
function run_subnets(
   parent_folder::String,
   objective::String;
   gen_areas::Array=[],
   area_transfer::Array=[],
   gen_zones::Array=[],
   zone_transfer::Array=[],
   enum_branches::Bool=false,
   plot_best_x::Int64=-1,
   scopf::Bool=false,
   ctg_plots::Array{Int64,1}=[0],
   run_fix_f::Bool=false,
   run_indir_PQ::Bool=false,
   print_results::Bool=false
   )
   folders = [f for f in readdir(parent_folder) if (isdir("$parent_folder/$f"))]
   if scopf
      results_dict = Dict{Int64, Dict{String, Dict{String, Any}}}()
   else
      results_dict = Dict{String, Dict{String, Any}}()
   end
   summary_df = DataFrame()
   folder_split = split(rstrip(parent_folder,'/'),"/")
   folder_split = folder_split[folder_split .!= ""]
   toplevels = folder_split[1:end-2]
   output_folder = join([("$i/") for i in toplevels])*"results/$(folder_split[end])"
   # output_folder = join([(i=="data" ? "results/" : "$i/") for i in split(parent_folder,"/")])
   println("output location: $output_folder")
   if !isdir(output_folder)
      mkpath(output_folder)
   end
   start_vals = Dict{String, Dict}("sn"=>Dict())
   if scopf
      for i in 1:length(folders)
         println("running $parent_folder/$(folders[i])  |  $i/$(length(folders)) folders")
         # full PQ case
         ctg_f_redispatch=false
         ctg_gen_redispatch=false
         ctg_iface_redispatch=true
         fix_post_ctg_pq=false
         direct_pq=true
         master_subnet=1
         suffix=""
         (result, summary_dict, solution_pm) = run_scopf(
            parent_folder,folders[i],
            objective, gen_areas, area_transfer, start_vals, results_dict,
            ctg_f_redispatch,
            ctg_gen_redispatch,
            ctg_iface_redispatch,
            fix_post_ctg_pq,
            direct_pq,
            master_subnet,
            suffix
         )
         if i==1
            start_vals = Dict("sn"=>Dict("0"=>solution_pm["sn"]["0"]))
         end
         res_summary = summary_dict[0]
         res_summary.network = rstrip(folders[i],'/')*suffix
         if i == 1
            summary_df = res_summary
         else
            summary_df = vcat(summary_df, res_summary)
         end

         # fixed post-contingency PQ
         ctg_f_redispatch=false
         ctg_gen_redispatch=false
         ctg_iface_redispatch=false
         fix_post_ctg_pq=true
         direct_pq=true
         master_subnet=1
         suffix="_fixf"
         (result, summary_dict, solution_pm) = run_scopf(
            parent_folder,folders[i],
            objective, gen_areas, area_transfer, start_vals, results_dict,
            ctg_f_redispatch,
            ctg_gen_redispatch,
            ctg_iface_redispatch,
            fix_post_ctg_pq,
            direct_pq,
            master_subnet,
            suffix
         )
         if i==1
            start_vals = Dict("sn"=>Dict("0"=>solution_pm["sn"]["0"]))
         end
         res_summary = summary_dict[0]
         res_summary.network = rstrip(folders[i],'/')*suffix
         if i == 1
            summary_df = res_summary
         else
            summary_df = vcat(summary_df, res_summary)
         end

         # passive post-contingency PQ
         # ctg_f_redispatch=false
         # ctg_gen_redispatch=false
         # ctg_iface_redispatch=false
         # fix_post_ctg_pq=false
         # direct_pq=true
         # master_subnet=1
         # suffix="_indirPQ"
         # (result, summary_dict, solution_pm) = run_scopf(
         #    parent_folder,folders[i],
         #    objective, gen_areas, area_transfer, start_vals, results_dict,
         #    ctg_f_redispatch,
         #    ctg_gen_redispatch,
         #    ctg_iface_redispatch,
         #    fix_post_ctg_pq,
         #    direct_pq,
         #    master_subnet,
         #    suffix
         # )
         # if i==1
         #    start_vals = Dict("sn"=>Dict("0"=>solution_pm["sn"]["0"]))
         # end
         # res_summary = summary_dict[0]
         # res_summary.network = rstrip(folders[i],'/')*suffix
         # if i == 1
         #    summary_df = res_summary
         # else
         #    summary_df = vcat(summary_df, res_summary)
         # end
      end
   else # if not scopf
      for i in 1:length(folders)
         println("running $parent_folder/$(folders[i])  |  $i/$(length(folders)) folders")
         # standard case
         fix_f_override=false
         direct_pq=true
         master_subnet=1
         suffix=""
         (result_temp, res_summary_temp, solution_pm_temp, results_dict_temp) = run_opf(
            parent_folder, folders[i],
            objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
            fix_f_override,
            direct_pq, master_subnet,
            suffix, print_results
         )

         if !(res_summary_temp.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED])
            fix_f_init = false
            noloss_override = create_noloss_dict(solution_pm_temp)

            println("Running OPF again without converter losses.")
            # Now run the opf again without converter losses by applying the noloss_override dictionary
            (result_init, res_summary_init, solution_pm_init, results_dict_init) = run_opf(
               parent_folder, folders[i],
               objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
               fix_f_override,
               direct_pq, master_subnet,
               suffix, print_results,
               override_param=noloss_override
            )
            if !(res_summary_init.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED])
               (result_init, res_summary_init, solution_pm_init, results_dict_init) = run_opf(
                  parent_folder, folders[i],
                  objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
                  true,
                  true, 1,
                  suffix, print_results
               )
               fix_f_init = true
               println("Using fixed frequency solution as initialization, running OPF again.")
            else
               println("Using no-loss solution as initialization, running OPF again with converter losses.")
            end
            # Run the opf again with converter losses and variable frequency, now using the solution from the no-loss or fixed frequency case as initialization

            (result_temp, res_summary_temp, solution_pm_temp, results_dict_temp) = run_opf(
               parent_folder, folders[i],
               objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
               fix_f_override,
               direct_pq, master_subnet,
               suffix, print_results,
               start_vals=solution_pm_init
            )

            if !(res_summary_temp.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED]) && !fix_f_init
               println("Running OPF with fixed frequency for initialization.")
               (result_init, res_summary_init, solution_pm_init, results_dict_init) = run_opf(
                  parent_folder, folders[i],
                  objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
                  true,
                  true, 1,
                  suffix, print_results
               )
               fix_f_init = true
               println("Using fixed frequency solution as initialization, running OPF again.")
               # Run the opf again with converter losses and variable frequency, now using the solution from the no-loss or fixed frequency case as initialization
               (result_temp, res_summary_temp, solution_pm_temp, results_dict_temp) = run_opf(
                  parent_folder, folders[i],
                  objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
                  fix_f_override,
                  direct_pq, master_subnet,
                  suffix, print_results,
                  start_vals=solution_pm_init
               )
               if !(res_summary_temp.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED]) && (res_summary_init.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED])
                  # use fixed frequency solution if it is feasible and the variable frequency solution is not
                  result_temp = result_init
                  res_summary_temp = res_summary_init
                  solution_pm_temp = solution_pm_init
                  results_dict_temp = results_dict_init
               end
            end
         end
         result = result_temp
         res_summary = res_summary_temp
         solution_pm = solution_pm_temp
         results_dict = results_dict_temp

         res_summary.network = rstrip(folders[i],'/')*suffix
         if i == 1
            summary_df = res_summary
         else
            summary_df = vcat(summary_df, res_summary)
         end

         if run_fix_f
            # fixed frequency case
            fix_f_override=true
            direct_pq=true
            master_subnet=1
            suffix="_fixf"
            (result, res_summary, solution_pm, results_dict) = run_opf(
               parent_folder, folders[i],
               objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
               fix_f_override,
               direct_pq, master_subnet,
               suffix, print_results
            )
            res_summary.network = rstrip(folders[i],'/')*suffix
            if i == 1
               summary_df = res_summary
            else
               summary_df = vcat(summary_df, res_summary)
            end
         end

         if run_indir_PQ
            # indirect PQ case
            fix_f_override=false
            direct_pq=false
            master_subnet=1
            suffix="_indirPQ"
            (result_temp, res_summary_temp, solution_pm_temp, results_dict_temp) = run_opf(
               parent_folder, folders[i],
               objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
               fix_f_override,
               direct_pq, master_subnet,
               suffix, print_results
            )

            if !(res_summary_temp.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED])
               fix_f_init = false
               noloss_override = create_noloss_dict(solution_pm_temp)

               println("Running OPF again without converter losses.")
               # Now run the opf again without converter losses by applying the noloss_override dictionary
               (result_init, res_summary_init, solution_pm_init, results_dict_init) = run_opf(
                  parent_folder, folders[i],
                  objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
                  fix_f_override,
                  direct_pq, master_subnet,
                  suffix, print_results,
                  override_param=noloss_override
               )
               if !(res_summary_init.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED])
                  (result_init, res_summary_init, solution_pm_init, results_dict_init) = run_opf(
                     parent_folder, folders[i],
                     objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
                     true,
                     direct_pq, master_subnet,
                     suffix, print_results
                  )
                  fix_f_init = true
                  println("Using fixed frequency solution as initialization, running OPF again.")
               else
                  println("Using no-loss solution as initialization, running OPF again with converter losses.")
               end
               # Run the opf again with converter losses and variable frequency, now using the solution from the no-loss or fixed frequency case as initialization

               (result_temp, res_summary_temp, solution_pm_temp, results_dict_temp) = run_opf(
                  parent_folder, folders[i],
                  objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
                  fix_f_override,
                  direct_pq, master_subnet,
                  suffix, print_results,
                  start_vals=solution_pm_init
               )

               if !(res_summary_temp.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED]) && !fix_f_init
                  println("Running OPF with fixed frequency for initialization.")
                  (result_init, res_summary_init, solution_pm_init, results_dict_init) = run_opf(
                     parent_folder, folders[i],
                     objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
                     true,
                     direct_pq, master_subnet,
                     suffix, print_results
                  )
                  fix_f_init = true
                  println("Using fixed frequency solution as initialization, running OPF again.")
                  # Run the opf again with converter losses and variable frequency, now using the solution from the no-loss or fixed frequency case as initialization
                  (result_temp, res_summary_temp, solution_pm_temp, results_dict_temp) = run_opf(
                     parent_folder, folders[i],
                     objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
                     fix_f_override,
                     direct_pq, master_subnet,
                     suffix, print_results,
                     start_vals=solution_pm_init
                  )
                  if !(res_summary_temp.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED]) && (res_summary_init.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED])
                     # use fixed frequency solution if it is feasible and the variable frequency solution is not
                     result_temp = result_init
                     res_summary_temp = res_summary_init
                     solution_pm_temp = solution_pm_init
                     results_dict_temp = results_dict_init
                  end
               end
            end
            result = result_temp
            res_summary = res_summary_temp
            solution_pm = solution_pm_temp
            results_dict = results_dict_temp

            res_summary.network = rstrip(folders[i],'/')*suffix
            if i == 1
               summary_df = res_summary
            else
               summary_df = vcat(summary_df, res_summary)
            end
         end

      end
   end

   println("Saving summary at $(output_folder)/summary.csv")
   CSV.write("$(output_folder)/summary.csv", summary_df)

   if enum_branches
      if scopf
         for (ctg_idx,ctg_result) in results_dict
            if ctg_idx in ctg_plots
               ctg_out_folder = "$output_folder/ctg$ctg_idx/"
               if !isdir(ctg_out_folder)
                  mkpath(ctg_out_folder)
               end
               collect_subnet_results(ctg_result, objective, gen_areas, ctg_out_folder, plot_best_x)
            end
         end
      else
         collect_subnet_results(results_dict, objective, gen_areas, gen_zones, output_folder, plot_best_x)
      end
   end
end


text:
function chess_n!(r)

    # перемещаем робота в левый нижний угол и запоминаем движение
    sides = []
    while !isborder(r, Sud)
        move!(r, Sud)
        push!(sides, Sud)
    end
    while !isborder(r, West)
        move!(r, West)
        push!(sides, West)
    end

    n = 4 # сторона одного квадрата
    count = 0
    k = 0
    flag = false
    while !flag
        for i in 1:n
            print_line(Ost, n, count)
            if isborder(r, Nord)
                flag = true
                break
            end
            move!(r, Nord)
        end
        if flag
            break
        end
        count += n
    end

    while !isborder(r, Sud)
        move!(r, Sud)
    end
    while !isborder(r, West)
        move!(r, West)
    end
    for i in sides
        move!(r, inverse_side(i))
    end
end

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs


text:
export plot_jacobian_spectrum

# multiply a peak by a number.
 function *(number::Float64, peak::Lorentzian) # multiplication operation
    Lorentzian(peak.ω0, peak.Γ, peak.A*number)
 end

 *(number::Float64, s::JacobianSpectrum) = JacobianSpectrum([number * peak for peak in s.peaks])

function show(io::IO, s::JacobianSpectrum)
    peaks = sort(s.peaks, by=x->x.ω0)
    print(io, "Lorentzian peaks (central frequency ω0, linewidth Γ): \n")
    for peak in peaks
        @printf(io, "%.3e * L(ω0 = %.6e, Γ = %.3e)\n", peak.A, peak.ω0, peak.Γ)
    end
end

function show(io::IO, spectra::Dict{Num,JacobianSpectrum})
    for var in keys(spectra)
        print("VARIABLE: ", var, "\n")
        show(spectra[var])
        print("\n")
    end
end

 "Adds a peak p to JacobianSpectrum s."
 function add_peak(s::JacobianSpectrum, p::Lorentzian)
    JacobianSpectrum(cat(s.peaks, p, dims=1))
 end

 "Adds all peaks from s2 to s1."
function add_peak(s1::JacobianSpectrum, s2::JacobianSpectrum)
    for p in s2.peaks
        s1 = add_peak(s1, p)
    end
    s1
end

"Gives the numerical value of a peak at ω."
evaluate(peak::Lorentzian, ω::Float64) = peak.A / sqrt(( (peak.ω0 - ω)^2 + (peak.Γ)^2 ))

"Gives the numerical value of a JacobianSpectrum at ω"
evaluate(s::JacobianSpectrum, ω::Float64) = sum([evaluate(p, ω) for p in s.peaks])

evaluate(spectra::Dict{Num, JacobianSpectrum}, ω::Float64) = Dict([[var, evaluate(spectra[var], ω)] for var in keys(spectra)])

"Converts an 2-component vector belonging to a VDP variable into the FT of the natural variable."
 function _VDP_to_spectrum(λ, eigvec::Vector; ω::Float64)
    u,v = eigvec
    peak1 = (1 + 2*(imag(u)*real(v) - real(u)*imag(v)) ) * Lorentzian(ω0=ω - imag(λ) , Γ=real(λ))
    peak2 = (1 + 2*(real(u)*imag(v) - real(v)*imag(u)) ) * Lorentzian(ω0=ω + imag(λ) , Γ=real(λ))
    JacobianSpectrum([peak1, peak2])
end

#   Returns the spectra of all variables in `res` for `index` of `branch`.
function JacobianSpectrum(res::Result; index::Int, branch::Int)

    res.classes["stable"][index][branch] || error("\nThe solution is unstable - it has no JacobianSpectrum!\n")

    solution_dict = get_single_solution(res, branch=branch, index=index)

    VDPs = res.problem.eom.variables # fetch the vector of HarmonicVariable
    λs, vs = eigen(res.jacobian(solution_dict))

    solution_dict = get_single_solution(res, index=index, branch=branch)

    # blank JacobianSpectrum for each variable
    all_spectra = Dict{Num, JacobianSpectrum}([[nat_var, JacobianSpectrum([])] for nat_var in getfield.(VDPs, Symbol("natural_variable"))])

    for (j, λ) in enumerate(λs)
        v = vs[:, j]
        # we know how to convert each VDP variable into an excitation -> break down the eigenvector into its constituent VDP pairs
        for (i, VDP_var) in enumerate(VDPs)
            v_pair = v[[2*i-1, 2*i]] # fetch the relevant part of the Jacobian eigenvector
            ωrot = Float64(substitute(VDP_var.ω, solution_dict))
            # v_pair is associated to a natural variable -> this variable gets Lorentzian peaks
            this_spec =  norm(v_pair) * _VDP_to_spectrum(λ, v_pair, ω=ωrot)

            all_spectra[VDP_var.natural_variable] = add_peak(all_spectra[VDP_var.natural_variable], this_spec)
        end
    end
    #_simplify_spectra!(all_spectra) # condense similar peaks
    all_spectra
end

"Condense peaks with similar centers and linewidths in JacobianSpectrum s."
function _simplify_spectrum(s::JacobianSpectrum)
    ω0s = getfield.(s.peaks, Symbol("ω0"))
    new_spectrum = JacobianSpectrum([])
    for ω in unique(ω0s)
        peaks = filter(x -> x.ω0 == ω, s.peaks) # all peaks centered at ω
        length(unique(getfield.(peaks, Symbol("Γ")))) > 1 && return(s)  # there are peaks of multiple linewidths at the same frequency => not simplifiable
        total_height = sqrt(sum(getfield.(peaks, Symbol("A")) .^2 )) # new height is the sum of squares
        new_spectrum = add_peak(new_spectrum, total_height * Lorentzian(ω0=ω,Γ=peaks[1].Γ) )
    end
    new_spectrum
end

"Simplify every JacobianSpectrum in a dictionary assigning variables to spectra."
function _simplify_spectra!(spectra::Dict{Num, JacobianSpectrum})
    for var in keys(spectra)
        spectra[var] = _simplify_spectrum(spectra[var])
    end
end

"""
    plot_jacobian_spectrum(res::Result,
     nat_var::Num;
      Ω_range,
       branch::Int,
        y_offset::String="0.0",
          x_scale=1.0,
           y_scale=1.0,
            logscale=false)

Make a 1D plot of the response spectrum of `res` for the natural variable `nat_var`.
Performs one matrix diagonalization for each element of `Ω_range`.
This method is faster than `plot_response` but results in errors where the noise frequency
is far from the frequency of the harmonic variables.
"""
function plot_jacobian_spectrum(res::Result, nat_var::Num; Ω_range, branch::Int, y_offset::String="0.0", x_scale=1.0, y_scale=1.0, logscale=false)

    length(size(res.solutions)) != 1 && error("1D plots of not-1D datasets are usually a bad idea.")
    stability = classify_branch(res, branch, "stable") # boolean array
    !any(stability) && error("Cannot generate a spectrum - no stable solutions!")

    X = Vector{Float64}(collect(values(res.swept_parameters))[1][stability])

    offset = Vector{Float64}(getindex.(transform_solutions(res, y_offset), branch))[stability]

    # only get spectra of the stable points!
    spectra = [JacobianSpectrum(res, branch=branch, index = i) for i in (1:length(res.solutions))[stability]]
    C = Array{Float64, 2}(undef,  length(Ω_range)-1, length(X)-1)

    bar = Progress(length(CartesianIndices(C)), 1, "Diagonalizing the Jacobian for each solution ... ", 50)
    for ij in CartesianIndices(C)
        C[ij] = abs(evaluate(spectra[ij[2]][nat_var], Ω_range[ij[1]] - offset[ij[2]]))
        next!(bar)
    end

    x_mat = x_scale * hcat([x*ones(length(Ω_range)) for x in X]...)
    y_mat = y_scale * hcat([Ω_range for j=1:length(X)]...)
    C = logscale ? log.(C) : C

    PyPlot.pcolormesh(x_mat, y_mat, C)
    xlabel(Latexify.latexify(string(first(keys(res.swept_parameters)))), fontsize=24);

    y_label = y_offset=="0.0" ? "noise " * latexify("ω") : "noise " * latexify("ω") * " - " * latexify(y_offset)
    ylabel(y_label, fontsize=24, fontname="Times");
    return C
end


text:
import PropertyCheck
using Base.Test
const pc = PropertyCheck

include("strategy_tests.jl")

max_stars_repo_name: cpritcha/PropertyCheck

max_issues_repo_name: cpritcha/PropertyCheck

max_forks_repo_name: cpritcha/PropertyCheck


text:
using ClassicalOrthogonalPolynomials, ContinuumArrays, QuasiArrays, BandedMatrices, Test
import ClassicalOrthogonalPolynomials: Hilbert, StieltjesPoint, ChebyshevInterval, associated, Associated, orthogonalityweight, Weighted, gennormalizedpower, *, dot, PowerLawMatrix, PowKernelPoint, LogKernelPoint
import InfiniteArrays: I

@testset "Associated" begin
    T = ChebyshevT()
    U = ChebyshevU()
    @test associated(T) ≡ U
    @test associated(U) ≡ U
    @test Associated(T)[0.1,1:10] == Associated(U)[0.1,1:10] == U[0.1,1:10]

    P = Legendre()
    Q = Associated(P)
    x = axes(P,1)
    u = Q * (Q \ exp.(x))
    @test u[0.1] ≈ exp(0.1)
    @test grid(Q[:,Base.OneTo(5)]) ≈ eigvals(Matrix(jacobimatrix(Normalized(Q))[1:5,1:5]))

    w = orthogonalityweight(Q)
    @test axes(w,1) == axes(P,1)
    @test sum(w) == 1
end

@testset "Singular integrals" begin
    @testset "weights" begin
        w_T = ChebyshevTWeight()
        w_U = ChebyshevUWeight()
        w_P = LegendreWeight()
        x = axes(w_T,1)
        H = inv.(x .- x')
        @test iszero(H*w_T)
        @test (H*w_U)[0.1] ≈ π/10
        @test (H*w_P)[0.1] ≈ log(1.1) - log(1-0.1)

        w_T = orthogonalityweight(chebyshevt(0..1))
        w_U = orthogonalityweight(chebyshevu(0..1))
        w_P = orthogonalityweight(legendre(0..1))
        x = axes(w_T,1)
        H = inv.(x .- x')
        @test iszero(H*w_T)
        @test (H*w_U)[0.1] ≈ (2*0.1-1)*π
        @test (H*w_P)[0.1] ≈ (log(1+(-0.8)) - log(1-(-0.8)))
    end

    @testset "LogKernelPoint" begin
        wU = Weighted(ChebyshevU())
        x = axes(wU,1)
        z = 0.1+0.2im
        L = log.(abs.(z.-x'))
        @test L isa LogKernelPoint{Float64,ComplexF64,ComplexF64,Float64,ChebyshevInterval{Float64}}

        @testset "Real point" begin
            U = ChebyshevU()
            x = axes(U,1)

            t = 2.0
            @test (log.(abs.(t .- x') )* Weighted(U))[1,1:3] ≈ [1.0362686329607178,-0.4108206734393296, -0.054364775221816465] #mathematica

            t = 0.5
            @test (log.(abs.(t .- x') )* Weighted(U))[1,1:3] ≈ [-1.4814921268505252, -1.308996938995747, 0.19634954084936207] #mathematica

            t = 0.5+0im
            @test (log.(abs.(t .- x') )* Weighted(U))[1,1:3] ≈ [-1.4814921268505252, -1.308996938995747, 0.19634954084936207] #mathematica
        end
    end

    @testset "Stieltjes" begin
        T = Chebyshev()
        wT = ChebyshevWeight() .* T
        x = axes(wT,1)
        z = 0.1+0.2im
        S = inv.(z .- x')
        @test S isa StieltjesPoint{ComplexF64,ComplexF64,Float64,ChebyshevInterval{Float64}}

        @test S * ChebyshevWeight() ≈ π/(sqrt(z-1)sqrt(z+1))
        @test S * JacobiWeight(0.1,0.2) ≈ 0.051643014475741864 - 2.7066092318596726im

        f = wT * [[1,2,3]; zeros(∞)];
        J = T \ (x .* T)

        @test π*((z*I-J) \ f.args[2])[1,1] ≈ (S*f)[1]
        @test π*((z*I-J) \ f.args[2])[1,1] ≈ (S*f.args[1]*f.args[2])[1]

        x = Inclusion(0..1)
        y = 2x .- 1
        wT2 = wT[y,:]
        S = inv.(z .- x')
        f = wT2 * [[1,2,3]; zeros(∞)];

        @test (π/2*(((z-1/2)*I - J/2) \ f.args[2]))[1] ≈ (S*f.args[1]*f.args[2])[1]

        @testset "Real point" begin
            t = 2.0
            T = ChebyshevT()
            U = ChebyshevU()
            x = axes(T,1)
            @test inv.(t .- x') * Weighted(T) ≈ inv.((t+eps()im) .- x') * Weighted(T)
            @test (inv.(t .- x') * Weighted(U))[1:10] ≈ (inv.((t+eps()im) .- x') * Weighted(U))[1:10]

            t = 2
            @test inv.(t .- x') * Weighted(T) ≈ inv.((t+eps()im) .- x') * Weighted(T)
            @test (inv.(t .- x') * Weighted(U))[1:10] ≈ (inv.((t+eps()im) .- x') * Weighted(U))[1:10]

            t = 0.5
            @test (inv.(t .- x') * Weighted(T))[1,1:3] ≈ [0,-π,-π]
            @test (inv.(t .- x') * Weighted(U))[1,1:3] ≈ [π/2,-π/2,-π]

            t = 0.5+0im
            @test (inv.(t .- x') * Weighted(T))[1,1:3] ≈ [0,-π,-π]
            @test (inv.(t .- x') * Weighted(U))[1,1:3] ≈ [π/2,-π/2,-π]
        end
    end

    @testset "Hilbert" begin
        wT = ChebyshevTWeight() .* ChebyshevT()
        wU = ChebyshevUWeight() .* ChebyshevU()
        x = axes(wT,1)
        H = inv.(x .- x')
        @test H isa Hilbert{Float64,ChebyshevInterval{Float64}}

        @testset "weights" begin
            @test H * ChebyshevTWeight() ≡ QuasiZeros{Float64}((x,))
            @test H * ChebyshevUWeight() == π*x
            @test (H * LegendreWeight())[0.1] ≈ log((0.1+1)/(1-0.1))
        end

        @test (Ultraspherical(1) \ (H*wT))[1:10,1:10] == diagm(1 => fill(-π,9))
        @test (Chebyshev() \ (H*wU))[1:10,1:10] == diagm(-1 => fill(1.0π,9))

        # check consistency
        @test (Ultraspherical(1) \ (H*wT) * (wT \ wU))[1:10,1:10] ==
                    ((Ultraspherical(1) \ Chebyshev()) * (Chebyshev() \ (H*wU)))[1:10,1:10]

        @testset "Other axes" begin
            x = Inclusion(0..1)
            y = 2x .- 1
            H = inv.(x .- x')

            wT2 = wT[y,:]
            wU2 = wU[y,:]
            @test (Ultraspherical(1)[y,:]\(H*wT2))[1:10,1:10] == diagm(1 => fill(-π,9))
            @test (Chebyshev()[y,:]\(H*wU2))[1:10,1:10] == diagm(-1 => fill(1.0π,9))
        end

        @testset "Legendre" begin
            P = Legendre()
            x = axes(P,1)
            H = inv.(x .- x')
            Q = H*P
            @test Q[0.1,1:3] ≈ [log(0.1+1)-log(1-0.1), 0.1*(log(0.1+1)-log(1-0.1))-2,-3*0.1 + 1/2*(-1 + 3*0.1^2)*(log(0.1+1)-log(1-0.1))]
            X = jacobimatrix(P)
            @test Q[0.1,1:11]'*X[1:11,1:10] ≈ (0.1 * Array(Q[0.1,1:10])' - [2 zeros(1,9)])
        end

        @testset "mapped" begin
            T = chebyshevt(0..1)
            U = chebyshevu(0..1)
            x = axes(T,1)
            H = inv.(x .- x')
            @test U\H*Weighted(T) isa BandedMatrix
        end
    end

    @testset "Log kernel" begin
        T = Chebyshev()
        wT = ChebyshevWeight() .* Chebyshev()
        x = axes(wT,1)
        L = log.(abs.(x .- x'))
        D = T \ (L * wT)
        @test ((L * wT) * (T \ exp.(x)))[0.] ≈ -2.3347795490945797  # Mathematica

        x = Inclusion(-1..1)
        T = Chebyshev()[1x, :]
        L = log.(abs.(x .- x'))
        wT = (ChebyshevWeight() .* Chebyshev())[1x, :]
        @test (T \ (L*wT))[1:10,1:10] ≈ D[1:10,1:10]

        x = Inclusion(0..1)
        T = Chebyshev()[2x.-1, :]
        wT = (ChebyshevWeight() .* Chebyshev())[2x .- 1, :]
        L = log.(abs.(x .- x'))
        u =  wT * (2 *(T \ exp.(x)))
        @test u[0.1] ≈ exp(0.1)/sqrt(0.1-0.1^2)
        @test (L * u)[0.5] ≈ -7.471469928754152 # Mathematica

        @testset "mapped" begin
            T = chebyshevt(0..1)
            x = axes(T,1)
            L = log.(abs.(x .- x'))
            @test T[0.2,:]'*((T\L*Weighted(T)) * (T\exp.(x))) ≈ -2.9976362326874373 # Mathematica
        end
    end

    @testset "pow kernel" begin
        P = Weighted(Jacobi(0.1,0.2))
        x = axes(P,1)
        S = abs.(x .- x').^0.5
        @test S isa ClassicalOrthogonalPolynomials.PowKernel
        @test_broken S*P
    end

    @testset "Ideal Fluid Flow" begin
        T = ChebyshevT()
        U = ChebyshevU()
        x = axes(U,1)
        H = inv.(x .- x')

        c = exp(0.5im)

        u = Weighted(U) * ((H * Weighted(U)) \ imag(c * x))

        ε  = eps();
        @test (inv.(0.1+ε*im .- x') * u + inv.(0.1-ε*im .- x') * u)/2 ≈ imag(c*0.1)
        @test real(inv.(0.1+ε*im .- x') * u ) ≈ imag(c*0.1)

        v = (s,t) -> (z = (s + im*t); imag(c*z) - real(inv.(z .- x') * u))
        @test v(0.1,0.2) ≈ 0.18496257285081724 # Emperical
    end

    @testset "OffHilbert" begin
        @testset "ChebyshevU" begin
            U = ChebyshevU()
            W = Weighted(U)
            t = axes(U,1)
            x = Inclusion(2..3)
            T = chebyshevt(2..3)
            H = T \ inv.(x .- t') * W;
            @test last(colsupport(H,1)) == 17
            @test last(colsupport(H,6)) ≤ 40
            @test T[2.3,1:100]'*(H * (W \ @.(sqrt(1-t^2)exp(t))))[1:100] ≈ 0.9068295340935111
            @test T[2.3,1:100]' * H[1:100,1:100] ≈ (inv.(2.3 .- t') * W)[:,1:100]
        end

        @testset "mapped" begin
            U = chebyshevu(-1..0)
            W = Weighted(U)
            t = axes(U,1)
            x = Inclusion(2..3)
            T = chebyshevt(2..3)
            H = T \ inv.(x .- t') * W
            N = 100
            @test T[2.3,1:N]' * H[1:N,1:N] ≈ (inv.(2.3 .- t') * W)[:,1:N]

            U = chebyshevu((-2)..(-1))
            W = Weighted(U)
            T = chebyshevt(0..2)
            x = axes(T,1)
            t = axes(W,1)
            H = T \ inv.(x .- t') * W
            @test T[0.5,1:N]'*(H * (W \ @.(sqrt(-1-t)*sqrt(t+2)*exp(t))))[1:N] ≈ 0.047390454610749054
        end
    end

    #################################################
    # ∫f(x)g(x)(t-x)^a dx evaluation where f and g in Legendre
    #################################################

    @testset "Pow kernel" begin
        @testset "Multiplication methods" begin
            P = Normalized(Legendre())
            x = axes(P,1)
            for (a,t) in ((0.1,1.2), (0.5,1.5))
                @test (t.-x).^a isa PowKernelPoint
                w = (t.-x).^a
                @test w .* P isa typeof(P*PowerLawMatrix(P,a,t))
            end
            # some functions
            f = P \ exp.(x.^2)
            g = P \ (sin.(x).*exp.(x.^(2)))
            # some parameters for (t-x)^a
            a = BigFloat("1.23")
            t = BigFloat("1.00001")
            # define powerlaw multiplication
            w = (t.-x).^a

            # check if it can compute the integral correctly
            @test g'*(P'*(w.*P)*f) ≈ -2.656108697646584 # Mathematica
        end
        @testset "Equivalence to multiplication in integer case" begin
            # TODO: overload integer input to make this work
            P = Normalized(Legendre())
            x = axes(P,1)
            a = 1
            t = 1.2
            @test_broken PowerLawMatrix(P,Float64(a),t)[1:20,1:20] ≈ ((t*I-jacobimatrix(P))^a)[1:20,1:20]
            a = 2
            t = 1.0001
            J = ((t*I-jacobimatrix(P)))[1:80,1:80]
            @test_broken PowerLawMatrix(P,BigFloat("$a"),BigFloat("$t"))[1:60,1:60] ≈ (J^2)[1:60,1:60]
        end
        @testset "Cached Legendre power law integral operator" begin
            P = Normalized(Legendre())
            a = 2*rand(1)[1]
            t = 1.0000000001
            Acached = PowerLawMatrix(P,BigFloat("$a"),BigFloat("$t"))
            @test size(Acached) == (∞,∞)
        end
        @testset "PowKernelPoint dot evaluation" begin
            @testset "Set 1" begin
                    P = Normalized(Legendre())
                    x = axes(P,1)
                    f = P \ abs.(π*x.^7)
                    g = P \ (cosh.(x.^3).*exp.(x.^(2)))
                    a = 1.9127
                    t = 1.211
                    w = (BigFloat("$t") .- x).^BigFloat("$a")
                    Pw = P'*(w .* P)
                    @test w isa PowKernelPoint
                    @test Pw[1:20,1:20] ≈ PowerLawMatrix(P,a,t)[1:20,1:20]
                    # this is slower than directly using PowerLawMatrix but it works
                    @test dot(f[1:20],Pw[1:20,1:20],g[1:20]) ≈ 5.082145576355614 # Mathematica
                end
            @testset "Set 2" begin
                P = Normalized(Legendre())
                x = axes(P,1)
                f = P \ exp.(x.^2)
                g = P \ (sin.(x).*exp.(x.^(2)))
                a = 1.23
                t = 1.00001
                W = PowerLawMatrix(P,a,t)
                @test dot(f,W,g) ≈ -2.656108697646584 # Mathematica
            end
            @testset "Set 3" begin
                P = Normalized(Legendre())
                x = axes(P,1)
                t = 1.2
                a = 1.1
                W = PowerLawMatrix(P,a,t)
                f = P \ exp.(x)
                g = P \ exp.(x.^2)
                @test dot(f,W,g) ≈ 2.916955525390389 # Mathematica
            end
            @testset "Set 4" begin
                P = Normalized(Legendre())
                x = axes(P,1)
                t = 1.001
                a = 1.001
                W = PowerLawMatrix(P,a,t)
                f = P \ (sinh.(x).*exp.(x))
                g = P \ cos.(x.^3)
                @test dot(f,W,g) ≈ -0.1249375144525209 # Mathematica
            end
            @testset "More explicit evaluation tests" begin
                # basis
                a = 2.9184
                t = 1.000001
                P = Normalized(Legendre())
                x = axes(P,1)
                # operator
                W = PowerLawMatrix(P,a,t)
                # functions
                f = P \ exp.(x)
                g = P \ sin.(x)
                const1(x) = 1
                onevec = P \ const1.(x)
                # dot() and * methods tests, explicit values via Mathematica
                @test -2.062500116206712 ≈ dot(onevec,W,g)
                @test 2.266485452423447 ≈ dot(onevec,W,f)
                @test -0.954305839543464 ≈ dot(g,W,f)
                @test 1.544769699288028 ≈ dot(f,W,f)
                @test 1.420460011606107 ≈ dot(g,W,g)
            end
        end
        @testset "Tests for -1 < a < 0" begin
            P = Normalized(Legendre())
            x = axes(P,1)
            a = -0.7
            t = 1.271
            # operator
            W = PowerLawMatrix(P,a,t)
            WB = PowerLawMatrix(P,BigFloat("$a"),BigFloat("$t"))
            # functions
            f0 = P \ exp.(2 .*x.^2)
            g0 = P \ sin.(x)
            @test dot(f0,W,g0) ≈  dot(f0,WB,g0) ≈ 1.670106472636101 # Mathematica
            f1 = P \ ((x.^2)./3 .+(x.^3)./3)
            g1 = P \ (x.*exp.(x.^3))
            @test dot(f1,W,g1) ≈ dot(f1,WB,g1) ≈ 0.5362428541997497 # Mathematica
        end
        @testset "Lanczos" begin
            P = Normalized(Legendre())
            x = axes(P,1)
            @time D = ClassicalOrthogonalPolynomials.LanczosData((1.001 .- x).^0.5, P);
            @time ClassicalOrthogonalPolynomials.resizedata!(D,100);
        end
    end
end


text:
ncclose()
v1=NetCDF.open(fn1,"v1",mode=NC_WRITE)
@test typeof(v1)==NetCDF.NcVar{Float64,3}
@test v1[1,1,1]==x1[1,1,1]
@test squeeze(v1[2,:,2],3)==x1[2,:,2]
@test squeeze(v1[1:end,3,1],(2,3))==x1[1:end,3,1]
inds=bitrand(size(x1))
@test v1[inds]==x1[inds]

xnew=rand(size(x1))
v1[:,:,:]=xnew
@test v1==xnew


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/PDEAssembler.jl/blob/master/LICENSE

# Poisson's problem as example

struct Poisson <: FieldProblem end

FEMBase.get_unknown_field_name(::Problem{Poisson}) = "u"

function FEMBase.assemble_elements!(problem::Problem{Poisson},
                                    assembly::Assembly,
                                    elements::Vector{Element{E}},
                                    time::Float64) where E

    bi = BasisInfo(E)
    ndofs = length(bi)
    Ke = zeros(ndofs, ndofs)
    fe = zeros(ndofs)

    for element in elements
        fill!(Ke, 0.0)
        fill!(fe, 0.0)
        for ip in get_integration_points(element)
            J, detJ, N, dN = element_info!(bi, element, ip, time)
            k = 1.0
            if haskey(element, "coefficient")
                k = element("coefficient", ip, time)
            end
            Ke += ip.weight * k*dN'*dN * detJ
            if haskey(element, "source")
                f = element("source", ip, time)
                fe += ip.weight * N'*f * detJ
            end
        end
        gdofs = get_gdofs(problem, element)
        add!(assembly.K, gdofs, gdofs, Ke)
        add!(assembly.f, gdofs, fe)
    end

    return nothing

end

function FEMBase.assemble_elements!(problem::Problem{Poisson},
                                    assembly::Assembly,
                                    elements::Vector{Element{E}},
                                    time::Float64) where E<:Union{Seg2,Seg3}

    bi = BasisInfo(E)
    ndofs = length(bi)
    Ce = zeros(ndofs, ndofs)
    fe = zeros(ndofs)
    ge = zeros(ndofs)

    for element in elements
        fill!(Ce, 0.0)
        fill!(fe, 0.0)
        fill!(ge, 0.0)
        for ip in get_integration_points(element)
            J, detJ, N, dN = element_info!(bi, element, ip, time)
            if haskey(element, "flux")
                f = element("flux", ip, time)
                fe += ip.weight * N'*f * detJ
            end
            if haskey(element, "fixed u")
                f = element("fixed u", ip, time)
                Ce += ip.weight * N'*N * detJ
                ge += ip.weight * N'*f * detJ
            end
        end
        gdofs = get_gdofs(problem, element)
        add!(assembly.C1, gdofs, gdofs, Ce')
        add!(assembly.C2, gdofs, gdofs, Ce)
        add!(assembly.f, gdofs, fe)
        add!(assembly.g, gdofs, ge)
    end

    return nothing

end

"""
    get_unit_square(nel_x=10, nel_y=10; lx=1.0, ly=1.0)

Return elements and boundary elements for a standard test problem, unit square.
"""
function get_unit_square(nel_x=20, nel_y=20; lx=1.0, ly=1.0)

    nnodes_x = nel_x+1
    nnodes_y = nel_y+1
    nnode = nnodes_x*nnodes_y
    nodemap = reshape(1:nnode, nnodes_x, nnodes_y)
    nodes_1 = vec(nodemap[1:nnodes_x-1, 1:nnodes_y-1])
    nodes_2 = vec(nodemap[2:nnodes_x, 1:nnodes_y-1])
    nodes_3 = vec(nodemap[2:nnodes_x, 2:nnodes_y])
    nodes_4 = vec(nodemap[1:nnodes_x-1, 2:nnodes_y])

    X = Dict{Int64, Vector{Float64}}()

    # Create nodes
    nid = 1
    for y in linspace(0, ly, nnodes_y)
        for x in linspace(0, lx, nnodes_x)
            X[nid] = [x, y]
            nid += 1
        end
    end

    # Create elements for volume
    field_elements = []
    for c in zip(nodes_1, nodes_2, nodes_3, nodes_4)
        push!(field_elements, Element(Quad4, collect(c)))
    end

    # add boundary elements to the left side of domain
    boundary_elements = []
    nids = 1:nel_x+1:(nel_x+1)*(nel_y+1) # LEFT and RIGHT
    for c in zip(nids[1:end-1], nids[2:end])
        connectivity = collect(c)
        push!(boundary_elements, Element(Seg2, connectivity))
        push!(boundary_elements, Element(Seg2, connectivity+nel_x))
    end
    nids = 1:nel_x+1 # BOTTOM and TOP
    for c in zip(nids[1:end-1], nids[2:end])
        connectivity = collect(c)
        push!(boundary_elements, Element(Seg2, connectivity))
        push!(boundary_elements, Element(Seg2, connectivity+(nel_x+1)*nel_y))
    end

    update!(field_elements, "geometry", X)
    update!(boundary_elements, "geometry", X)

    return field_elements, boundary_elements
end


text:
"""
    compute_linear_return_array(priceArray::Array{Float64,1})

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function compute_linear_return_array(priceArray::Array{Float64,1})::PSResult

    # initialize -
    linear_return_array = Array{Float64,1}()

    # get the size -
    number_of_time_steps = length(priceArray)
    for time_index = 2:number_of_time_steps
        value = (priceArray[time_index] / priceArray[time_index - 1]) - 1
        push!(linear_return_array, value)
    end

    # return -
    return PSResult(linear_return_array)
end

"""
    compute_linear_return_array(dataTable::DataFrame; key::Symbol = Symbol("adjusted_close"))

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function compute_linear_return_array(dataTable::DataFrame; key::Symbol=Symbol("adjusted_close"))

    # initialize -
    price_array = Float64[]

    # iterate to build the price array -
    (number_of_rows, number_of_cols) = size(dataTable)
    for row_index = 1:number_of_rows
        data_value = dataTable[row_index,key]
        push!(price_array, data_value)
    end

    # return -
    return compute_linear_return_array(price_array)
end

"""
    compute_log_return_array(priceArray::Array{Float64})

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function compute_log_return_array(priceArray::Array{Float64})::PSResult

    # initialize -
    log_return_array = Array{Float64,1}()

    # get the size -
    number_of_time_steps = length(priceArray)
    for time_index = 2:number_of_time_steps
        value = log((priceArray[time_index] / priceArray[time_index - 1]))
        push!(log_return_array, value)
    end

    # return -
    return PSResult(log_return_array)
end

"""
    compute_log_return_array(dataTable::DataFrame; key::Symbol = Symbol("adjusted_close"))

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function compute_log_return_array(dataTable::DataFrame; key::Symbol=Symbol("adjusted_close"))::PSResult

    # initialize -
    price_array = Float64[]

    # iterate to build the price array -
    (number_of_rows, number_of_cols) = size(dataTable)
    for row_index = 1:number_of_rows
        data_value = dataTable[row_index,key]
        push!(price_array, data_value)
    end

    # return -
    return compute_log_return_array(price_array)
end

"""
    compute_return_volatility(dataTable::DataFrame; returnCalcFunction::Function=compute_linear_return_array,
        key::Symbol = Symbol("adjusted_close"))

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function compute_return_volatility(dataTable::DataFrame; returnCalcFunction::Function=compute_linear_return_array,
    key::Symbol=Symbol("adjusted_close"))::PSResult

    # initialize -
    price_array = Float64[]

    # iterate to build the price array -
    (number_of_rows, number_of_cols) = size(dataTable)
    for row_index = 1:number_of_rows
        data_value = dataTable[row_index,key]
        push!(price_array, data_value)
    end

    # we have the difference -
    result = returnCalcFunction(price_array)
    if (isa(result.value, Exception) == true)
        return result
    end
    price_return_array = result.value

    # compute the std -
    volatlity = std(price_return_array)

    # return -
    return PSResult(volatlity)
end

"""
    compute_return_volatility(priceArray::Array{Float64};
        returnCalcFunction::Function=compute_linear_return_array)

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function compute_return_volatility(priceArray::Array{Float64};
    returnCalcFunction::Function=compute_linear_return_array)::PSResult

    # compute the return array -
    result = returnCalcFunction(priceArray)
    if (isa(result.value, Exception) == true)
        return result
    end
    volatlity = std(result.value)

    # return -
    return PSResult(volatlity)
end


text:
# Everything should be the same as usual normals except
# when x is out of boundaries

p = 50
a = zeros(p)
b = ones(p)

mu  = rand(p)
P1  = randn(p, p)
P1 *= P1'
C1  = inv(P1)
C1 += C1'
C1 /= 2

mvg   = MvGaussianCanon(mu, P1)
mvg2  = MvGaussianNatural(P1*mu, -P1)
mvg_d = MvDiagonalGaussian(mu, 2ones(p))
mvgs  = MvGaussianStandard(mu, C1)

tmvg   = TMvGaussianCanon(mu, P1, a, b)
tmvg2  = TMvGaussianNatural(P1*mu, -P1, a, b)
tmvg_d = TMvDiagonalGaussian(mu, 2ones(p), a, b)
tmvgs  = TMvGaussianStandard(mu, C1, a, b)

x_valid = rand(p)

@test loglik(mvg, x_valid) == loglik(tmvg, x_valid)
@test loglik(mvg2, x_valid) == loglik(tmvg2, x_valid)
@test loglik(mvg_d, x_valid) == loglik(tmvg_d, x_valid)
@test loglik(mvgs, x_valid) == loglik(tmvgs, x_valid)

@test gradloglik(mvg, x_valid) == gradloglik(tmvg, x_valid)
@test gradloglik(mvg2, x_valid) == gradloglik(tmvg2, x_valid)
@test gradloglik(mvg_d, x_valid) == gradloglik(tmvg_d, x_valid)
@test gradloglik(mvgs, x_valid) == gradloglik(tmvgs, x_valid)

x_not = -rand(p)

@test loglik(tmvg, x_not) == -Inf
@test loglik(tmvg2, x_not) == -Inf
@test loglik(tmvg_d, x_not) == -Inf
@test loglik(tmvgs, x_not) == -Inf

infvec = Inf*ones(p)

@test gradloglik(tmvg, x_not) == infvec
@test gradloglik(tmvg2, x_not) == infvec
@test gradloglik(tmvg_d, x_not) == infvec
@test gradloglik(tmvgs, x_not) == infvec


text:
"""
    unwrap(A[, wrap_around, seed])

Unwraps the values in A modulo 2π, where A is a 1-, 2-, or 3- dimensional array.

# Arguments
- `A::AbstractArray{T, N}`: Array to unwrap
- `wrap_around::Tuple{Bool}`:  When an element of this tuple is `true`, the
    unwrapping process will consider the edges along the corresponding axis
    of the image to be connected
- `seed::Int`: Unwrapping of 2D or 3D images uses a random initialization. This
    sets the seed of the RNG.
"""
function unwrap end

"""
    unwrap!(A[, wrap_around, seed])

In-place version of unwrap.
"""
function unwrap! end

function unwrap(A::AbstractArray{T, N},
                wrap_around::NTuple{N, Bool}=tuple(zeros(Bool, N)...),
                seed::Int=-1) where {T,N}
    A_copy = copy(A)
    unwrap!(A_copy, wrap_around, seed)
    return A_copy
end

max_forks_repo_name: platawiec/Unwrap


text:
using Geodesy
using Base.Test
using Geodesy: Bounds, onBounds, inBounds, boundaryPoint, NAD27, center

# Construction

# inBounds

function test_bounds{T}(bounds::Bounds{T})

    min_x, min_y, max_x, max_y = bounds.min_x, bounds.min_y, bounds.max_x, bounds.max_y

    @test inBounds(T(min_x, min_y), bounds)
    @test !inBounds(T(min_x - eps(min_x), min_y), bounds)
    @test !inBounds(T(min_x, min_y - eps(min_y)), bounds)

    @test inBounds(T(max_x, max_y), bounds)
    @test !inBounds(T(max_x + eps(max_x), max_y), bounds)
    @test !inBounds(T(max_x, max_y + eps(max_y)), bounds)
end

for bounds in (Bounds{ENU}(1.1, 2.2, 3.3, 4.4),
               Bounds{LLA{NAD27}}(1.1, 2.2, 3.3, 4.4),
               )
    test_bounds(bounds)
end

# boundaryPoint

function test_boundary{T}(bounds::Bounds{T})
    c = center(bounds)
    cx, cy = getX(c), getY(c)

    for _ = 1:1_000
        in_both =    T(cx + rand() - 0.5, cy + rand() - 0.5)
        in_x1 =      T(cx + rand() - 0.5, cy + rand() + 1.0)
        in_x2 =      T(cx + rand() - 0.5, cy - rand() - 1.0)
        in_y =       T(cx + rand() + 1.0, cy + rand() - 0.5)
        in_neither = T(cx + rand() + 1.0, cy + rand() + 1.0)

        @test onBounds(boundaryPoint(in_both, in_x1, bounds), bounds)
        @test onBounds(boundaryPoint(in_x2, in_both, bounds), bounds)
        @test onBounds(boundaryPoint(in_both, in_y, bounds), bounds)
        @test onBounds(boundaryPoint(in_neither, in_both, bounds), bounds)
    end
end

for bounds in (Bounds{LLA{NAD27}}(0, 1, 78, 79),
               Bounds{ENU}(-1, 1, -1, 1))
    test_boundary(bounds)
end


text:
module BLPData

using Dates
using Printf

include("deps.jl")
include("types.jl")
include("results.jl")
include("correlation_id.jl")
include("datetime.jl")
include("c.jl")
include("constant.jl")
include("blpname.jl")
include("session.jl")
include("service.jl")
include("schema.jl")
include("element.jl")
include("request.jl")
include("message.jl")
include("event_queue.jl")
include("print.jl")
include("bdh.jl")
include("bdh_intraday_ticks.jl")
include("bds.jl")
include("bdp.jl")
include("subscription.jl")

function __init__()
    check_deps()
end

@inline function error_check(code::Cint)
    if code != 0
        error(blpapi_getLastErrorDescription(code))
    end

    nothing
end

@inline function error_check(code::Cint, ctx_msg::AbstractString)
    if code != 0
        error("$ctx_msg: $(blpapi_getLastErrorDescription(code))")
    end

    nothing
end

@inline function ptr_check(ptr::Ptr{Cvoid})
    if ptr == C_NULL
        error("Got a NULL pointer.")
    end

    nothing
end

@inline function ptr_check(ptr::Ptr{Cvoid}, ctx_msg::AbstractString)
    if ptr == C_NULL
        error("$ctx_msg: got a NULL pointer.")
    end

    nothing
end

"""
    get_version_info() :: VersionInfo

Returns the version of the shared library for Bloomberg API.
"""
function get_version_info() :: VersionInfo
    major = Ref{Cint}(0)
    minor = Ref{Cint}(0)
    patch = Ref{Cint}(0)
    build = Ref{Cint}(0)
    blpapi_getVersionInfo(major, minor, patch, build)
    return VersionInfo(major[], minor[], patch[], build[])
end

end # module

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
using Base.SparseArrays.CHOLMOD

function chm_analyze_ordering{Tv<:CHOLMOD.VTypes}(A::CHOLMOD.Sparse{Tv}, ordering::Cint, Perm::Vector{Int})
    s = unsafe_load(A.p)
    Parent   = Array(Int,s.nrow)
    Post     = Array(Int,s.nrow)
    ColCount = Array(Int,s.nrow)
    First    = Array(Int,s.nrow)
    Level    = Array(Int,s.nrow)

    ccall((@CHOLMOD.cholmod_name("analyze_ordering",Int),:libcholmod), Cint,
          (Ptr{CHOLMOD.C_Sparse{Tv}}, Cint, Ptr{Int}, Ptr{Int}, Csize_t,
           Ptr{Int}, Ptr{Int}, Ptr{Int},
           Ptr{Int}, Ptr{Int}, Ptr{UInt8}),
           A.p, ordering, Perm, C_NULL, zero(Csize_t), Parent, Post, ColCount, First, Level, CHOLMOD.common())
    Parent, Post, ColCount
end

function pothen_sun{T<:Integer}(par::Vector{T}, post::Vector{T}, colcount::Vector{T})
    n = length(par)
    flag = -ones(T,n)
    ns = n
    for j=post
        mdeg = colcount[j] - 1
        if par[j] != 0 && mdeg == colcount[par[j]] && flag[par[j]] == -1
            # par[j] not assigned to supernode
            ns -= 1
            if flag[j] < 0    # j is a repr. vertex
                flag[par[j]] = j
                flag[j] -= 1
            else              # j is not a repr. vertex
                flag[par[j]] = flag[j]
                flag[flag[j]] -= 1
            end
        end
    end
    ns, flag
end

function chordal_embedding{Tv<:Number,Ti<:Int}(A::SparseMatrixCSC{Tv,Ti}, Perm::Vector{Ti})
    m, n = size(A)
    S = CHOLMOD.Sparse(round(Float64,A))

    F = cholfact(round(Float64,A), shift=n, perm=Perm)
    s = unsafe_load(F.p)
    p = round(Int64,[ unsafe_load(s.Perm, i) for i=1:s.n])

    par, post, colcount = chm_analyze_ordering(S, s.ordering, p)
    ns, flag = pothen_sun(par+1, post+1, colcount)

    # extract cliques
    L = sparse(CHOLMOD.Sparse(F))
    Array{Int,1}[ sort(1+p[L.rowval[L.colptr[i]:L.colptr[i+1]-1]]) for i = find(flag .< 0) ]
end

# If no permutation is specified, use the CHOLMOD ordering
chordal_embedding{Tv<:Number,Ti<:Int}(A::SparseMatrixCSC{Tv,Ti}) = chordal_embedding(A, Array(Int,0))

function clique_index(I::Array{Array{Int,1},1}, S::Array{Int,1}, n::Int)

    J = Int[]
    mask = zeros(Int, n)
    s = sparsevec(S, 1, n)
    for k=1:length(I)
        if length(S) <= length(I[k])
            mask[I[k]] = 1
            if dot(mask, s) == length(S)
                push!(J, k)
            end
            mask[I[k]] = 0
        end
    end
    J
end


text:
module HealthBase

export get_fhir_access_token
export get_fhir_encounter_id
export get_fhir_patient_id
export has_fhir_encounter_id
export has_fhir_patient_id

include("smart_authorization.jl")

end # module