text:
@with_kw struct CSetup{S,V}
    a::S = 0.5u"cm^2"
    A::S = 71.0u"cm^2"
    fixed_ΔV::V = nothing
end

function specific_capacitance(C, porosity, folder, a, A)
    porosities = sort!(parse.(Int, dirs_in_folder(folder, false)))

    df = CSV.read("$folder/weights.csv")
    mass = (df[!,:after].*u"g" .- df[!,:before].*u"g") .* a ./ A .|> u"μg"
    mass_dict = Dict(porosities.=>mass[1:length(porosities)])

    # @show mass_dict
    C / mass_dict[porosity]
end

function energy(C, ΔV)
    return (C*ΔV^2)/2
end

power(E, Δt) = E / Δt

function add_report!(df, cr)
    add_aux!(df, cr)

    push!(df[!, :C], ustrip(u"μF", cr.C))
    push!(df[!, :C_specific], ustrip(u"F/g", cr.C_specific))
    push!(df[!, :E], ustrip(u"W*hr", cr.E))
    push!(df[!, :E_specific], ustrip(u"W*hr/kg", cr.E_specific))
    push!(df[!, :P], ustrip(u"W", cr.P))
    push!(df[!, :P_specific], ustrip(u"W/kg", cr.P_specific))

    return nothing
end

function add_report!(result_df, datafile, folder, setup)
    if endswith(datafile.filename, "_D")
        df = read_file(datafile)
        cr = CDCapacitanceReport(datafile, df, folder, setup)
        add_report!(result_df, cr)
    end
end

function add_report!(result_df, datafile, quadrant, folder, setup)
    df = read_file(datafile, 3, false)
    cr = CVCapacitanceReport(datafile, df, quadrant, folder, setup)
    add_report!(result_df, cr)
end

function compute_capacitances(folder; cv_setup=[CSetup(),CSetup()], cd_setup=CSetup())
    data = find_files(folder)
    processed_data = find_files(folder, exclude_with=r"!", select_with=".dat", rename=false)

    grouped = groupbyfolder(data["C&D"])
    processed_grouped = groupbyfolder(processed_data["CV"])
    @assert keys(grouped) == keys(processed_grouped)

    cd_units = cd_report_units()
    cv_units = cv_report_units()

    for f in keys(grouped)
        cd_datafiles = grouped[f]
        cv_datafiles = processed_grouped[f]

        cd_report = cd_result()
        cv_report4 = cv_result()
        cv_report2 = cv_result()

        for datafile in cd_datafiles
            add_report!(cd_report, datafile, folder, cd_setup)
        end

        for datafile in cv_datafiles
            add_report!(cv_report4, datafile, 4, folder, cv_setup[1])
            add_report!(cv_report2, datafile, 2, folder, cv_setup[2])
        end

        rename_cols!(df) = rename!(df, Dict(
            "E_specific"=>"Energy density",
            "P_specific"=>"Power density",
            "C"=>"C_abs",
            "C_specific"=>"C",
            "Δt"=>replace_unicode("Δt"),
            "ΔV"=>replace_unicode("ΔV")))
        rename_cols!(cd_report)
        rename_cols!(cv_report4)
        rename_cols!(cv_report2)
        sort!(cd_report, :I)
        sort!(cv_report4, :scan_rate)
        sort!(cv_report2, :scan_rate)

        comment(df) = join(repeat([f*replace_powers(" mA cm^-2")], size(df,2)), ',')

        write_file(
            cd_report,
            (cd_units, comment(cd_report)),
            joinpath(folder, f, "cd_capacitances.csv"),
            ',')
        write_file(
            cv_report4,
            (cv_units, comment(cv_report4)),
            joinpath(folder, f, "cv_capacitances4.csv"),
            ',')
        write_file(
            cv_report2,
            (cv_units, comment(cv_report2)),
            joinpath(folder, f, "cv_capacitances2.csv"),
            ',')
    end

    return nothing
end


text:
#CALIBRATION FILE
# Parameters used into the processing + Calibrations protocol functions
#---------------------------------------------------------------------

#List of Calibration Values, allow the soft to do EVERYTHING.
Cal = Dict{String,Int16}()
#Format : "CalibrationVar" => Value(Int16)

#coefR,coefG,coefB and bias are neural parameter to filter solder grey.
#density,order are paramters for the lowpass noise filter (after the RGBtoWHITE processing).
#spray is the amount of decimation to use while doing the AprxConv processing.
#match is the percentage needed to be considerer as a footprint after the Conv processing.
#t_OFF is the Temperature Offset for the (un)soldering iron.

#Default Calibration Values.
DefaultCal = Dict{String,Int16}("coefR"=>12,"coefG"=>12,"coefB"=>12,"bias"=>1,
"density"=>51,"order"=>3,"spray"=>2,"match"=>80,"t_OFF"=>0)

"""
	\tLoadCal(CalibrationFile::String)
	\t\tReturn Boolean
Check if there is EBoAR Calibration File corresponding to the specified argument.\n
If YES, the function will gather the parameters and put them into a "Cal" dictionnary variable.
If NO, the function will create a file corresponding to the specified argument and fill it with Default values.
Then, put theses values into the "Cal" variable.
"""
function LoadCal(adr::String)
	#search Calibration File ("Calibration.txt")
	#If here, fill Cal with all the written values.
	#If missing, create it with default value into the current directory and fill Cal with default values.

	if(isfile(adr))
		println("Calibration file found.")
		open(adr,"r") do c_file
			while ((!eof(c_file))&&isreadable(c_file))
				(k,v) = split(readline(c_file))
				try
					Cal[k] = parse(Int16,v)
					println(k*" => "*v)
				catch
					println("\t -> Wrong number format (need Int16)")
					return false
				end
			end
		end
		return (keys(Cal)==keys(DefaultCal))
	else
		println("Calibration file no found.")
		println("Creating with default values...")
		open(adr,"w+") do c_file
			for k in keys(DefaultCal)
				try
					write(c_file,k*"\t\t"*string(DefaultCal[k])*"\n")
				catch
					println("Cannot write into calibration file !")
					return false
				end
			end
		end
		return true
	end
end

ext: jl


text:
#=##############################################################################
# DESCRIPTION
    Methods for manipulating airfoil geometries.
# AUTHORSHIP
  * Author    : Eduardo J Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Jul 2018
  * License   : MIT License
=###############################################################################

"Receives a .dat file as pulled from airfoiltools.com containing the x and y
contour coordinates of an airfoil, and returns arrays x and y."
function readcontour(file_name; header_len=1, delim=" ", path="", output="arrays")

    x, y = Float64[], Float64[]

    open(joinpath(path,file_name)) do f
        for (i,line) in enumerate(eachline(f))

            # Ignores header
            if i<=header_len
                nothing
                # Parses each line
            else
                this_x, this_y = split(line, delim; keepempty=false)
                push!(x, Base.parse(Float64, this_x))
                push!(y, Base.parse(Float64, this_y))
            end

        end
    end

    if output=="arrays"
        return x,y
    elseif output=="matrix"
        return hcat(x,y)
    else
        error("Invalid `output` argument $(output).")
    end
end
#TODO: This is type unstable, since it returns either a vector or a matrix.  Should
# we make this type stable?

function rediscretize_airfoil(x::AbstractVector{T}, y::AbstractVector{T},
    upperNDIVS::multidisctype, lowerNDIVS::multidisctype; spl_s::Real=0.00001,
    spl_k::Integer=4, verify_spline::Bool=true, pltargs...) where{T<:Real}

    # Separate upper and lower sides to make the contour injective in x
    upper, lower = splitcontour(x, y)

    # Parameterize both sides independently
    fun_upper = parameterize(upper[1], upper[2], zeros(eltype(upper[1]), size(upper[1])); inj_var=1,
                                                      s=spl_s, kspl=spl_k)
    fun_lower = parameterize(lower[1], lower[2], zeros(eltype(lower[1]), size(lower[1])); inj_var=1,
                                                      s=spl_s, kspl=spl_k)

    # Discretizes both sides
    new_upper = multidiscretize(fun_upper, 0, 1, upperNDIVS)
    new_lower = multidiscretize(fun_lower, 0, 1, lowerNDIVS)

    # Merges sides back together
    points = vcat(reverse(new_upper), new_lower[2:end])
    new_x = [p[1] for p in points]
    new_y = [p[2] for p in points]

    # Plots
    if verify_spline
        if isdefined(Main, :PyPlot)
            plot_airfoil(x, y; label="Original", style="--^k", alpha=0.5, pltargs...)
            plot_airfoil(new_x, new_y; label="Parameterized", style=":.b", pltargs...)
        else
            @eval Main import PyPlot
            Base.invokelatest(plot_airfoil, x, y; label="Original", style="--^k", alpha=0.5, pltargs...)
            Base.invokelatest(plot_airfoil, new_x, new_y; label="Parameterized", style=":.b", pltargs...)
        end
    end

    return new_x, new_y
end

"""
    Returns the y-coordinates of upper and lower surfaces of a 4-digit NACA
    airfoil. x input and y outputs are normalized by the chord.

    * First digit describing maximum camber as percentage of the chord.
    * Second digit describing the distance of maximum camber from the airfoil
            leading edge in tens of percents of the chord.
    * Last two digits describing maximum thickness of the airfoil as percent of
            the chord.

    https://en.wikipedia.org/wiki/NACA_airfoil#Four-digit_series
"""
function naca4digits(d1, d2, d34, xs)

    m = d1/100
    p = d2/10
    t = d34/100

    yt = 5*t * (0.2969*sqrt.(xs) .- 0.126*xs .- 0.3516*xs.^2 .+ 0.2843*xs.^3 .- 0.1015*xs.^4)

    if d1==0 && d2==0
        return xs, xs, yt, -yt

    elseif d1==0 || d2==0
        error("Invalid airfoil definition d1=$d1, d2=$d2")

    else

        yc = similar(xs)
        thetas = similar(xs)

        for (i,x) in enumerate(xs)
            if x<0 || x>1
                error("Got invalid x-coordinate: $x")
            elseif x <= p
                yc[i] = m/p^2 * (2*p*x - x^2)
                thetas[i] = atan(2*m/p^2 * (p-x))
            else
                yc[i] = m/(1-p)^2 * ((1-2*p) + 2*p*x - x^2)
                thetas[i] = atan(2*m/(1-p)^2 * (p-x))
            end
        end

        xU = xs .- yt.*sin.(thetas)
        xL = xs .+ yt.*sin.(thetas)
        yU = yc .+ yt.*cos.(thetas)
        yL = yc .- yt.*cos.(thetas)

        return xU, xL, yU, yL
    end
end

function naca4digits(d1::Int64, d2::Int64, d34::Int64, n::Int64, r::Float64;
                                                                    off::Real=0)
    xs = [float(x) for x in discretize(x->x, off, 1, n, r; central=true)]
    xs[end] = 1.0
    xU, xL, yU, yL = naca4digits(d1, d2, d34, xs)
    xU[end] = 1.0
    xL[end] = 1.0

    if off==0
        return vcat(reverse(xU), xL[2:end]), vcat(reverse(yU), yL[2:end])
    else
        return vcat(reverse(xU), xL), vcat(reverse(yU), yL)
    end
end


text:
using MLJ, StableRNGs
import DataFrames
Ridge = @load RidgeRegressor pkg=MultivariateStats

rng = StableRNG(551234) # for reproducibility

x1 = rand(rng, 300)
x2 = rand(rng, 300)
x3 = rand(rng, 300)
y = exp.(x1 - x2 -2x3 + 0.1*rand(rng, 300))

X = DataFrames.DataFrame(x1=x1, x2=x2, x3=x3)
first(X, 3) |> pretty

test, train = partition(eachindex(y), 0.8);

Xs = source(X)
ys = source(y)

stand = machine(Standardizer(), Xs)
W = transform(stand, Xs)

fit!(W, rows=train);

W()             # transforms all data
W(rows=test, )  # transforms only test data
W(X[3:4, :])    # transforms specific data

box_model = UnivariateBoxCoxTransformer()
box = machine(box_model, ys)
z = transform(box, ys)

ridge_model = Ridge(lambda=0.1)
ridge = machine(ridge_model, W, z)
ẑ = predict(ridge, W)

ŷ = inverse_transform(box, ẑ)

fit!(ŷ, rows=train);

rms(y[test], ŷ(rows=test))

ridge_model.lambda = 5.0;

fit!(ŷ, rows=train)
rms(y[test], ŷ(rows=test))

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
export CIFAR100
module CIFAR100
    using DataDeps
    using BinDeps
    using FixedPointNumbers
    using ..MLDatasets: bytes_to_type, datafile, download_dep, download_docstring
    using ..CIFAR10: convert2image, convert2features

    export

        classnames_coarse,
        classnames_fine,

        traintensor,
        testtensor,

        trainlabels,
        testlabels,

        traindata,
        testdata,

        convert2image,
        convert2features,

        download

    const DEPNAME = "CIFAR100"
    const TRAINSET_FILENAME = joinpath("cifar-100-binary", "train.bin")
    const TESTSET_FILENAME  = joinpath("cifar-100-binary", "test.bin")
    const COARSE_FILENAME = joinpath("cifar-100-binary", "coarse_label_names.txt")
    const FINE_FILENAME = joinpath("cifar-100-binary", "fine_label_names.txt")

    const TRAINSET_SIZE = 50_000
    const TESTSET_SIZE  = 10_000

    """
        download([dir]; [i_accept_the_terms_of_use])

    Trigger the (interactive) download of the full dataset into
    "`dir`". If no `dir` is provided the dataset will be
    downloaded into "~/.julia/datadeps/$DEPNAME".

    This function will display an interactive dialog unless
    either the keyword parameter `i_accept_the_terms_of_use` or
    the environment variable `DATADEPS_ALWAYS_ACCEPT` is set to
    `true`. Note that using the data responsibly and respecting
    copyright/terms-of-use remains your responsibility.
    """
    download(args...; kw...) = download_dep(DEPNAME, args...; kw...)

    include(joinpath("Reader","Reader.jl"))
    include("interface.jl")

    function __init__()
        register(DataDep(
            DEPNAME,
            """
            Dataset: The CIFAR-100 dataset
            Authors: Alex Krizhevsky, Vinod Nair, Geoffrey Hinton
            Website: https://www.cs.toronto.edu/~kriz/cifar.html
            Reference: https://www.cs.toronto.edu/~kriz/learning-features-2009-TR.pdf

            [Krizhevsky, 2009]
                Alex Krizhevsky.
                "Learning Multiple Layers of Features from Tiny Images",
                Tech Report, 2009.

            The CIFAR-100 dataset is a labeled subsets of the 80
            million tiny images dataset. It consists of 60000
            32x32 colour images in 100 classes. Specifically, it
            has 100 classes containing 600 images each. There are
            500 training images and 100 testing images per class.
            The 100 classes in the CIFAR-100 are grouped into 20
            superclasses. Each image comes with a "fine" label
            (the class to which it belongs) and a "coarse" label
            (the superclass to which it belongs).

            The compressed archive file that contains the
            complete dataset is available for download at the
            offical website linked above; specifically the binary
            version for C programs. Note that using the data
            responsibly and respecting copyright remains your
            responsibility. The authors of CIFAR-100 aren't
            really explicit about any terms of use, so please
            read the website to make sure you want to download
            the dataset.
            """,
            "https://www.cs.toronto.edu/~kriz/cifar-100-binary.tar.gz",
            "58a81ae192c23a4be8b1804d68e518ed807d710a4eb253b1f2a199162a40d8ec",
            post_fetch_method = file -> (run(BinDeps.unpack_cmd(file,dirname(file), ".gz", ".tar")); rm(file))
        ))
    end
end


text:
## --- Test bilinear exponential function

    @test bilinear_exponential(66.02, [1.219384666893631, 66.0846112095295, 0.10822726922099106, -1.277645745450562, 0.6971298010153942]) ≈ 0.15736531960405353
    @test bilinear_exponential([66.02, 66.08, 66.15], [1.219384666893631, 66.0846112095295, 0.10822726922099106, -1.277645745450562, 0.6971298010153942]) ≈ [0.1573653196040532, 4.614683898246965, 8.99724709626313e-8]
    @test bilinear_exponential_ll(66.02, [1.219384666893631, 66.0846112095295, 0.10822726922099106, -1.277645745450562, 0.6971298010153942]) ≈ -1.8491853001597212

## --


text:
using Test,RSCG
using SparseArrays
using LinearMaps

function test1(n,M)
    σ = zeros(ComplexF64,M)
    σmin = -10.0*im
    σmax = 10.0*im
    for i=1:M
        σ[i] = (i-1)*(σmax-σmin)/(M-1) + σmin
    end
    A = make_mat(n)
    i = 1
    j = 1
    Gij = greensfunctions(i,j,σ,A)

    Ii = spzeros(ComplexF64,n,n)
    for ii=1:n
        Ii[ii,ii] = σ[1]
    end
    G = inv(Matrix(Ii-A))

    @test abs(G[1,1]-Gij[1]) < 1e-7
end

function test1(n,M)
    σ = zeros(ComplexF64,M)
    σmin = -10.0*im
    σmax = 10.0*im
    for i=1:M
        σ[i] = (i-1)*(σmax-σmin)/(M-1) + σmin
    end
    A = make_mat(n)
    i = 1
    j = 1
    Gij = greensfunctions(i,j,σ,A)

    Ii = spzeros(ComplexF64,n,n)
    for ii=1:n
        Ii[ii,ii] = σ[1]
    end
    G = inv(Matrix(Ii-A))
    println("Residual ",abs(G[1,1]-Gij[1]))
    @test abs(G[1,1]-Gij[1]) < 1e-7
end

function test2(n,M)
    σ = zeros(ComplexF64,M)
    σmin = -10.0*im
    σmax = 10.0*im
    for i=1:M
        σ[i] = (i-1)*(σmax-σmin)/(M-1) + σmin
    end
    A2 = make_matc(n)
    i = 1
    j = 1
    Gij = greensfunctions(i,j,σ,A2)

    Ii = spzeros(ComplexF64,2n,2n)
    for ii=1:2n
        Ii[ii,ii] = σ[1]
    end
    G = inv(Matrix(Ii-A2))
    println("Residual ",abs(G[1,1]-Gij[1]))
    @test abs(G[1,1]-Gij[1]) < 1e-7
end

function test3(n,M)
    σ = zeros(ComplexF64,M)
    σmin = -4.0*im+0.2
    σmax = 4.0*im -0.3
    for i=1:M
        σ[i] = (i-1)*(σmax-σmin)/(M-1) + σmin
    end
    A2 = make_matc(n)
    vec_i = [1,4,div(n,2),n]
    j = div(n,2)
    Gij = greensfunctions(vec_i,j,σ,A2)

    Ii = spzeros(ComplexF64,2n,2n)
    for ii=1:2n
        Ii[ii,ii] = σ[div(M,2)]
    end
    G = inv(Matrix(Ii-A2))

    for ii = 1:length(vec_i)
        i = vec_i[ii]
        println("Residual ",abs(G[i,j]-Gij[div(M,2),ii]) )
        @test abs(G[i,j]-Gij[div(M,2),ii]) < 1e-7
    end
end

function set_diff(v)
    function calc_diff!(y::AbstractVector, x::AbstractVector)
        n = length(x)
        length(y) == n || throw(DimensionMismatch())
        μ = -1.5
        for i=1:n
            dx = 1
            jp = i+dx
            jp += ifelse(jp > n,-n,0) #+1方向
            dx = -1
            jm = i+dx
            jm += ifelse(jm < 1,n,0) #-1方向
            y[i] = v*(x[jp]+x[jm])-μ*x[i]
        end

        return y
    end
    (y,x) -> calc_diff!(y,x)
end

function test4(n,M)
    σ = zeros(ComplexF64,M)
    σmin = -4.0*im-1.2
    σmax = 4.0*im +4.3
    for i=1:M
        σ[i] = (i-1)*(σmax-σmin)/(M-1) + σmin
    end
    A = set_diff(-1.0)
    D = LinearMap(A, n; ismutating=true,issymmetric=true)
    A2 = make_mat(n)

    vec_i = [2,4,div(n,2),n]
    j = n
    Gij = greensfunctions(vec_i,j,σ,D)

    Ii = spzeros(ComplexF64,n,n)
    for ii=1:n
        Ii[ii,ii] = σ[div(M,2)]
    end
    G = inv(Matrix(Ii-A2))

    for ii = 1:length(vec_i)
        i = vec_i[ii]
        println("Residual ",abs(G[i,j]-Gij[div(M,2),ii]) )
        @test abs(G[i,j]-Gij[div(M,2),ii]) < 1e-7
    end
end

function make_mat(n)
    A = spzeros(Float64,n,n)
    t = -1.0
    μ = -1.5
    for i=1:n
        dx = 1
        jp = i+dx
        jp += ifelse(jp > n,-n,0) #+1方向
        dx = -1
        jm = i+dx
        jm += ifelse(jm < 1,n,0) #-1方向
        A[i,jp] = t
        A[i,i] = -μ
        A[i,jm] = t
    end
    return A
end

function make_matc(n)
    A = make_mat(n)
    A2 = spzeros(ComplexF64,2n,2n)

    Delta = spzeros(ComplexF64,n,n)
    phi = rand(n)
    for i=1:n
        Delta[i,i] = exp(im*phi[i])
    end

    A2[1:n,1:n] = A[1:n,1:n]
    A2[1+n:2n,1+n:2n] = -conj(A[1:n,1:n])
    A2[1:n,1+n:2n] = Delta[1:n,1:n]
    A2[1+n:2n,1:n] = Delta'[1:n,1:n]

    return A2
end

@testset "SparseArrays" begin
    @testset "real" begin
        M = 100
        σ = zeros(M)
        n = 10
        println("n = $(n)")
        @time test1(n,M)
        n = 100
        println("n = $(n)")
        @time test1(n,M)
        n=  1000
        println("n = $(n)")
        @time test1(n,M)
    end

    @testset "complex" begin
        M = 100
        σ = zeros(M)
        n = 5
        println("n = $(2n)")
        @time  test2(n,M)
        n = 50
        println("n = $(2n)")
        @time  test2(n,M)
        n=  500
        println("n = $(2n)")
        @time  test2(n,M)
    end

    @testset "complex: multi values" begin
        M = 100
        σ = zeros(M)
        n = 5
        println("n = $(2n)")
        @time  test3(n,M)
        n = 50
        println("n = $(2n)")
        @time  test3(n,M)
        n=  500
        println("n = $(2n)")
        @time  test3(n,M)
    end
end

@testset "Linearmaps" begin
    @testset "complex: multi values" begin
        M = 100
        σ = zeros(M)
        n = 10
        println("n = $(n)")
        @time  test4(n,M)
        n = 100
        println("n = $(n)")
        @time  test4(n,M)
        n=  1000
        println("n = $(n)")
        @time  test4(n,M)
    end
end


text:
facts("sensekeys") do
    const mock_db = DB(joinpath(dirname(@__FILE__), "mock_db"))

    lem = mock_db["section", 'n']
    ss = synsets(mock_db, lem)
    @fact sensekeys(mock_db, lem) --> ["section%1:15:01::"]
    @fact sensekey(mock_db, ss[1], lem) --> "section%1:15:01::"
end


text:
"variable controling a linear genetor responce "
function variable_response_delta(pm::_PM.AbstractPowerModel; nw::Int=nw_id_default, report::Bool=true)
    delta = var(pm, nw)[:delta] = @variable(pm.model,
        base_name="$(nw)_delta",
        start = 0.0
    )

    if report
        sol(pm, nw)[:delta] = delta
    end
end

"generates variables for both `active` and `reactive` bus deltas"
function variable_bus_delta_abs(pm::_PM.AbstractPowerModel; kwargs...)
    variable_bus_delta_abs_power_real(pm; kwargs...)
    variable_bus_delta_abs_power_imaginary(pm; kwargs...)
end

""
function variable_bus_delta_abs_power_real(pm::_PM.AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    p_delta_abs = var(pm, nw)[:p_delta_abs] = @variable(pm.model,
        [i in ids(pm, :bus)], base_name="$(nw)_p_delta_abs",
        start = 0.0
    )

    if bounded
        for (i, bus) in ref(pm, nw, :bus)
            JuMP.set_lower_bound(p_delta_abs[i], 0.0)
            JuMP.set_upper_bound(p_delta_abs[i], 0.5)
        end
    end

    report && _PM.sol_component_value(pm, nw, :bus, :p_delta_abs, ids(pm, nw, :bus), p_delta_abs)
end

""
function variable_bus_delta_abs_power_imaginary(pm::_PM.AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
     q_delta_abs = var(pm, nw)[:q_delta_abs] = @variable(pm.model,
        [i in ids(pm, :bus)], base_name="$(nw)_q_delta_abs",
        start = 0.0
    )

    if bounded
        for (i, bus) in ref(pm, nw, :bus)
            JuMP.set_lower_bound(q_delta_abs[i], 0.0)
            JuMP.set_upper_bound(q_delta_abs[i], 0.5)
        end
    end

    report && _PM.sol_component_value(pm, nw, :bus, :q_delta_abs, ids(pm, nw, :bus), q_delta_abs)
end

""
function variable_shunt_admittance_imaginary(pm::_PM.AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    bs = var(pm, nw)[:bs] = @variable(pm.model,
        [i in ids(pm, nw, :shunt_var)], base_name="$(nw)_bs",
        start = _PM.comp_start_value(ref(pm, nw, :shunt, i), "bs_start")
    )

    if bounded
        for i in ids(pm, nw, :shunt_var)
            shunt = ref(pm, nw, :shunt, i)
            JuMP.set_lower_bound(bs[i], shunt["bmin"])
            JuMP.set_upper_bound(bs[i], shunt["bmax"])
        end
    end

    report && _PM.sol_component_value(pm, nw, :shunt, :bs, ids(pm, nw, :shunt_var), bs)
end

""
function variable_shunt_admittance_imaginary(pm::_PM.AbstractWModels; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    bs = var(pm, nw)[:bs] = @variable(pm.model,
        [i in ids(pm, nw, :shunt_var)], base_name="$(nw)_bs",
        start = _PM.comp_start_value(ref(pm, nw, :shunt, i), "bs_start")
    )

    wbs = var(pm, nw)[:wbs] = @variable(pm.model,
        [i in ids(pm, nw, :shunt_var)], base_name="$(nw)_wbs",
        start = 0.0
    )

    if bounded
        for i in ids(pm, nw, :shunt_var)
            shunt = ref(pm, nw, :shunt, i)
            JuMP.set_lower_bound(bs[i], shunt["bmin"])
            JuMP.set_upper_bound(bs[i], shunt["bmax"])
        end
    end

    report && _PM.sol_component_value(pm, nw, :shunt, :bs, ids(pm, nw, :shunt_var), bs)
    report && _PM.sol_component_value(pm, nw, :shunt, :wbs, ids(pm, nw, :shunt_var), wbs)
end

function variable_branch_power_slack(pm::_PM.AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    sm_slack = var(pm, nw)[:sm_slack] = JuMP.@variable(pm.model,
        [l in ids(pm, nw, :branch_sm_active)], base_name="$(nw)_sm_slack",
        start = _PM.comp_start_value(ref(pm, nw, :branch, l), "sm_slack_start")
    )

    if bounded
        for (l,branch) in ref(pm, nw, :branch_sm_active)
            JuMP.set_lower_bound(sm_slack[l], 0.0)
        end
    end

    report && _PM.sol_component_value(pm, nw, :branch, :sm_slack, ids(pm, nw, :branch_sm_active), sm_slack)
end

function variable_bus_voltage_magnitude_delta(pm::_PM.AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    vvm_delta = var(pm, nw)[:vvm_delta] = JuMP.@variable(pm.model,
        [i in ids(pm, nw, :bus)], base_name="$(nw)_vvm_delta",
        start = _PM.comp_start_value(ref(pm, nw, :bus, i), "vvm_delta_start", 0.1)
    )

    report && _PM.sol_component_value(pm, nw, :bus, :vvm_delta, ids(pm, nw, :bus), vvm_delta)
end

function variable_gen_power_real_delta(pm::_PM.AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    pg_delta = var(pm, nw)[:pg_delta] = JuMP.@variable(pm.model,
        [i in ids(pm, nw, :gen)], base_name="$(nw)_pg_delta",
        start = _PM.comp_start_value(ref(pm, nw, :gen, i), "pg_delta_start", 0.1)
    )

    report && _PM.sol_component_value(pm, nw, :gen, :pg_delta, ids(pm, nw, :gen), pg_delta)
end

function variable_branch_contigency_power_violation(pm::_PM.AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    branch_cont_flow_vio = var(pm, nw)[:branch_cont_flow_vio] = JuMP.@variable(pm.model,
        [i in 1:length(ref(pm, :branch_flow_cuts))], base_name="$(nw)_branch_cont_flow_vio",
        #start = _PM.comp_start_value(ref(pm, nw, :bus, i), "cont_branch_vio_start")
    )

    if bounded
        for i in 1:length(ref(pm, :branch_flow_cuts))
            JuMP.set_lower_bound(branch_cont_flow_vio[i], 0.0)
        end
    end

    #report && _PM.sol_component_value(pm, nw, :gen, :pg_delta, ids(pm, nw, :gen), pg_delta)
end

function variable_gen_contigency_power_violation(pm::_PM.AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    gen_cont_flow_vio = var(pm, nw)[:gen_cont_flow_vio] = JuMP.@variable(pm.model,
        [i in 1:length(ref(pm, :gen_flow_cuts))], base_name="$(nw)_gen_cont_flow_vio",
        #start = _PM.comp_start_value(ref(pm, nw, :bus, i), "cont_branch_vio_start")
    )

    if bounded
        for i in 1:length(ref(pm, :gen_flow_cuts))
            JuMP.set_lower_bound(gen_cont_flow_vio[i], 0.0)
        end
    end

    #report && _PM.sol_component_value(pm, nw, :gen, :pg_delta, ids(pm, nw, :gen), pg_delta)
end

function variable_gen_contigency_capacity_violation(pm::_PM.AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    gen_cont_cap_vio = var(pm, nw)[:gen_cont_cap_vio] = JuMP.@variable(pm.model,
        [i in 1:length(ref(pm, :gen_contingencies))], base_name="$(nw)_gen_cont_cap_vio",
        #start = _PM.comp_start_value(ref(pm, nw, :bus, i), "cont_branch_vio_start")
    )

    if bounded
        for i in 1:length(ref(pm, :gen_contingencies))
            JuMP.set_lower_bound(gen_cont_cap_vio[i], 0.0)
        end
    end

    #report && _PM.sol_component_value(pm, nw, :gen, :pg_delta, ids(pm, nw, :gen), pg_delta)
end


text:
"""
cond(m::MixedModel)

Return a vector of condition numbers of the λ matrices for the random-effects terms
"""
LinearAlgebra.cond(m::MixedModel) = cond.(m.λ)

function σs(m::MixedModel)
    σ = dispersion(m)
    NamedTuple{fnames(m)}(((σs(t, σ) for t in m.reterms)...,))
end

function σρs(m::MixedModel)
    σ = dispersion(m)
    NamedTuple{fnames(m)}(((σρs(t, σ) for t in m.reterms)...,))
end

"""
    vcov(m::LinearMixedModel)

Returns the variance-covariance matrix of the fixed effects.
If `corr=true`, then correlation of fixed effects is returned instead.
"""
function StatsBase.vcov(m::MixedModel; corr=false)
    Xtrm = fetrm(m)
    iperm = invperm(Xtrm.piv)
    p = length(iperm)
    r = Xtrm.rank
    Linv = inv(feL(m))
    T = eltype(Linv)
    permvcov = dispersion(m, true) * (Linv'Linv)
    if p == Xtrm.rank
        vv = permvcov[iperm, iperm]
    else
        covmat = fill(zero(T) / zero(T), (p, p))
        for j = 1:r, i = 1:r
            covmat[i, j] = permvcov[i, j]
        end
        vv = covmat[iperm, iperm]
    end

    corr ?  StatsBase.cov2cor!(vv, stderror(m)) : vv
end


text:
using Morsel
using Meddle
using JSON
using WebSockets
using REPL
using Gadfly
using Compose
using Color
import YAML

f = open("settings.yml")
settings = YAML.load(f)
close(f)

const dark_theme =
    Gadfly.Theme(color("steel blue"),      # default_color
          0.6mm,                    # default_point_size
          0.2mm,                    # line_width
          color("#000000"),         # panel_fill
          color("#4d87c7"),         # panel_stroke
          color("#333333"),         # grid_color
          color("#f0f0f0"),         # grid_color_focused
          0.2mm,                    # grid line width
          Gadfly.default_font_desc,        # minor_label_font
          9pt,                      # minor_label_font_size
          color("#dddddd"),         # minor_label_color
          Gadfly.default_font_desc,        # major_label_font
          11pt,                     # major_label_font_size
          color("#dddddd"),         # major_label_color
          Gadfly.default_font_desc,        # point_label_font
          8pt,                      # point_label_font_size
          color("#dddddd"),         # point_label_color
          0.0mm,                    # bar_spacing
          1mm,                      # boxplot_spacing
          0.3mm,                    # highlight_width
          Gadfly.default_highlight_color,  # highlight_color
          Gadfly.default_middle_color,     # middle_color
          1000,                     # label_placement_iterations
          10.0,                     # label_out_of_bounds_penalty
          0.5,                      # label_hidden_penalty
          0.2)                      # label_visibility_flip_pr

include("msg_types.jl")

app = Morsel.app()

import Base: repl_show

function repl_show(io::WebSocket,result,user_id)
	message = JSON.to_json([MSG_OUTPUT_EVAL_RESULT,user_id,sprint(repl_show,result)])
	for sock in clients
		write(sock,message)
	end
end
function repl_show(io::WebSocket,p::Union(Plot,Canvas),user_id)
	p.theme = dark_theme
	div_id = "plot_$(dec(uint16(rand(Uint64))))"
	out = IOBuffer()
	draw(D3(out,6inch,6inch; emit_on_finish=false),p)
	message = JSON.to_json([MSG_OUTPUT_HTML,user_id,"""

	"""]) #"
	for sock in clients
		isopen(sock) && write(sock,message)
	end
end

repl_channel, response_channel = RemoteRef(),RemoteRef()
REPL.start_repl_backend(repl_channel,response_channel)

clients = Array(WebSocket,0)

websocket(app,"/") do req, res
	socket = res.res
	push!(clients,socket)
	while true
        msg = read(socket)
        julia_msg = JSON.parse(bytestring(msg))
        msg_type = julia_msg[1]
        if msg_type == MSG_INPUT_EVAL
        	username = julia_msg[2]
        	user_id = julia_msg[3]
        	ast = Base.parse_input_line(julia_msg[4])
        	put(repl_channel, (ast,1))
        	@async begin
        		message = JSON.to_json([MSG_OUTPUT_EVAL_INPUT,user_id,username,julia_msg[4]])
        		for sock in clients
        			isopen(sock) && write(sock, message)
        		end
        		(val, bt) = take(response_channel)
        		repl_show(socket,val,user_id)
        	end
        end
    end
    splice!(clients,findfirst(clients,socket))
end

start(app, parseint(setting["worker_port"]))


text:
using LoopVectorization: Instruction

"""
Contains a sequence of instructions and dependencies for each instruction.
Dependencies are vectors of integers, indicating all the arguments for the instruction.
Arguments ≥ 1 are intermediaries.
Argument 0 is the gradient with respect to the output
Arguments -N,...,-1 are the inputs
Arguments -2N,...,-N-1 are the gradients with respect to inputs, for use with updating operations (zero-initialized is common case)

Sections are:
1: End first pass block
2: End shared second pass block
2+1,...,2+n,...,2+N: End block for argument n

Returns:
1: return value of function
1+1,...,1+n,...,1+N: return value for gradient with respect to argument n

A few examples:

    exp(x)
instructions = [   :exp,     :fmadd    ]
dependencies = [ [ -1 ], [ 1,0,-2 ] ]
sections = [ 1:1, 2:1, 2:2 ]
returns  = [ 1, 2 ]

Means that we use the following instruction sequence for calculating the value and derivative of `y = exp(x)`:
# definitions
%-1 = x
%0 = ∂y
# first pass
%1 = exp(%-1)
# shared
# calc ∂x
%2 = %1 * %0

then `y = %1` and `∂x = %2`.

    atan(x, y)
instructions = [   atan,  abs2, abs2,   +,     /,     -,    *,      *   ]
dependencies = [ [-2,-1], [-2], [-1], [1,2], [0,4], [-2], [7,5], [-1,5] ]
sections = [ 1:1, 2:5, 6:7, 8:8 ]
returns = [ 1, 7, 8 ]

This corresponds to the following instruction sequence for `z = atan(x,y)`:
# definitions
%-2 = x
%-1 = y
%0 = ∂z
# first pass
%1 = atan(%-2, %-1)
# shared
%2 = abs2(%-2)
%3 = abs2(%-1)
%4 = %2 + %3
%5 = %0 / %4
# calc ∂x
%6 = - %-2
%7 = %6 * %5
# calc ∂y
%8 = %-1 * %5

With
z = %1
∂x = %7
∂y = %8

    sincos(x)
instructions = [ sincos, first,  last,    -,     *,     *,    tuple   ]
dependencies = [  [-1],  [ 1 ], [ 1 ], [ 2 ], [3,0], [4,0], [3, 4] ]
sections = [ 1:1, 2:1, 2:7 ]
returns = [ 1, 7 ]

For `y = sincos(x)`
# definitions
%-1 = x
%0 = ∂y
# first pass
%1 = sincos(x)
# shared
# calc ∂x
%2 = first(%1)
%3 = last(%1)
%4 = -%2
%5 = %3 * %0
%6 = %4 * %0
%7 = tuple( %5, %6 )

"""
struct DiffRule
    instructions::Vector{Instruction}
    dependencies::Vector{Vector{Int}}
    sections::Vector{UnitRange{Int}}
    returns::Vector{Int}
    num_parents::Int
end

LoopVectorization.instruction(dr::DiffRule, i::Int) = dr.instructions[i]
dependencies(dr::DiffRule, i::Int) = dr.dependencies[i]

struct InstructionArgs
    instr::Instruction
    num_args::Int
end
# function InstructionArgs(instr::Instruction, args)
    # uargs = zero(UInt)
    # for b ∈ args
        # uargs <<= 1
        # uargs |= b
    # end
    # InstructionArgs(instr, uargs, length(args))
# end
# InstructionArgs(instr::Symbol, num_args) = InstructionArgs(Instruction(instr), args)
num_parents(dr::DiffRule) = dr.num_parents
section_two_returns(dr::DiffRule) = length(dr.sections) == length(dr.returns)
function Base.hash(ia::InstructionArgs, h::UInt)
    hash(ia.instr, hash(ia.num_args, h))
end
function Base.:(==)(ia1::InstructionArgs, ia2::InstructionArgs)
    instr1 = ia1.instr; instr2 = ia2.instr
    # Rough order of likelihood they differ
    instr1.instr === instr2.instr && ia1.num_args == ia2.num_args && instr1.mod === instr2.mod
end

const DERIVATIVERULES = Dict{InstructionArgs,DiffRule}()
# Sometimes we may want to drop additions and set equal to a value rather than update.
# """
# Drop last argument to the function
# """
# const NOADDITION = Dict{Instruction,Instruction}(
#     Instruction(:+) => Instruction(:identity),
#     Instruction(:vadd) => Instruction(:identity),
#     Instruction(:vnsub) => Instruction(:vsub),
#     Instruction(:muladd) => Instruction(:*),
#     Instruction(:vmuladd) => Instruction(:vmul),
#     Instruction(:vfma) => Instruction(:vmul),
#     Instruction(:vfmadd) => Instruction(:vmul),
#     Instruction(:vfmadd_fast) => Instruction(:vmul),
#     Instruction(:vfnmadd) => Instruction(:vnmul),
#     Instruction(:vfnmadd_fast) => Instruction(:vnmul),
#     Instruction(:vmullog2add) => Instruction(:vmullog2),
#     Instruction(:vmullog10add) => Instruction(:vmullog10),
#     Instruction(:vdivlog2add) => Instruction(:vdivlog2),
#     Instruction(:vdivlog10add) => Instruction(:vdivlog10)
# )
# To support mutating and non-mutating APIs.
# Generated code should be SAC-style
const MAKEUPDATING = Dict{Instruction,Instruction}(
    Instruction(:identity) => Instruction(:vadd!),
    Instruction(:vsub) => Instruction(:vsub!),
    Instruction(:*) => Instruction(:vfmadd!),
    Instruction(:vmul) => Instruction(:vfmadd!),
    Instruction(:vnmul) => Instruction(:vfnmadd!),
    Instruction(:vmullog2) => Instruction(:vmullog2add!),
    Instruction(:vmullog10) => Instruction(:vmullog10add!),
    Instruction(:vdivlog2) => Instruction(:vdivlog2add!),
    Instruction(:vdivlog10) => Instruction(:vdivlog10add!)
)
const MAKEINPLACE = Dict{Instruction,Instruction}(
    Instruction(:*) => Instruction(:mul!),
    Instruction(:vmul) => Instruction(:mul!),
    Instruction(:vnmul) => Instruction(:nmul!)
)

for arity ∈ 2:8
    let sections = [ i > 1 ? (i:i) : (i+1):1 for i ∈ 0:arity+1 ], returns = collect(1:arity+1), deps = [ iszero(i) ? collect(-arity:-1) : [0] for i in 0:arity ]
        for add ∈ [:+, :vadd, :evadd]
            instructions = Instruction[ iszero(i) ? Instruction(add) : Instruction(:identity) for i ∈ 0:arity ]
            DERIVATIVERULES[InstructionArgs(add, arity)] = DiffRule(
                instructions, deps, sections, returns, arity
            )
        end
    end
    # let sections = [ 1:1, 2:1 ], returns = [ iszero(i) ? 1 : i - 1 - 2arity for i ∈ 0:arity ]
    #     DERIVATIVERULES[InstructionArgs(:logdensity, arity)] = DiffRule(
    #         [ :∂logdensity! ], deps, sections, returns, arity
    #     )
    # end
end
let deps = [ [-2, -1], [-1], [0,2], [-2], [4,0] ], sections = [ 1:1, 2:1, 2:3, 4:5], returns = [ 1, 3, 5 ]
    for mul ∈ [:*, :vmul, :evmul]
        DERIVATIVERULES[InstructionArgs(mul, 2)] = DiffRule(
            Instruction[ mul, :adjoint, :vmul, :adjoint, :vmul ],
            deps, sections, returns, 2
        )
    end
end
let deps = [ [ -2,-1], [0], [0] ], sections = [1:1,2:1,2:2,3:3], returns = [1,2,3]
    for sub ∈ [:-, :vsub, :evsub]
        instructions = Instruction[sub, :identity, :vsub]
        DERIVATIVERULES[InstructionArgs(sub,2)] = DiffRule(
            instructions, deps, sections, returns, 2
        )
    end
end
let deps = [ [-2,-1], [-1], [0,2], [1], [4,3,-3] ], sections = [ 1:1, 2:3, 4:3, 4:5 ], returns = [1,3,5]
    for div ∈ [:/, :vfdiv]
        instructions = [ div, :adjoint, :vfdiv, :adjoint, :vnmul ]
        DERIVATIVERULES[InstructionArgs(div,2)] = DiffRule(
            instructions, deps, sections, returns, 2
        )
    end
end

# fmadd for the ∂ updates so compiler can elide zero-initialization
let deps = [ [-3, -2, -1], [-2], [0,2], [-3], [4,0] ], sections = [ 1:1, 2:1, 2:3, 4:5, 6:5 ], returns = [ 1, 3, 5, 0 ]
    for fma_instr ∈ [:muladd, :fma, :vmuladd, :vfma, :vfmadd, :vfmadd_fast]
       # fma_instr2 = :vfmadd_fast
        DERIVATIVERULES[InstructionArgs(fma_instr, 3)] = DiffRule(
            Instruction[ fma_instr, :adjoint, :vmul, :adjoint, :vmul ],
            deps, sections, returns, 3
        )
    end
    for vfnmadd_instr ∈ [:vfnmadd, :vfnmadd_fast]
        # fnmaadd_instr2 = :vfnmadd_fast
        DERIVATIVERULES[InstructionArgs(vfnmadd_instr, 3)] = DiffRule(
            Instruction[ vfnmadd_instr, :adjoint, :vnmul, :adjoint, :vnmul ],
            deps, sections, returns, 3
        )
    end
end
let deps = [ [-3, -2, -1], [-2], [0,2], [-3], [4,0], [0] ], sections = [ 1:1, 2:1, 2:3, 4:5, 6:6 ], returns = [ 1, 3, 5, 6 ]
    for vfmsub_instr ∈ [:vfmsub, :vfmsub_fast]
        # fmadd = vfmsub_instr === :vfmsub ? :vfmadd : :vfmadd_fast
        DERIVATIVERULES[InstructionArgs(vfmsub_instr, 3)] = DiffRule(
            Instruction[ vfmsub_instr, :adjoint, :vmul, :adjoint, :vmul, :vsub ],
            deps, sections, returns, 3
        )
    end
    for vfnmsub_instr ∈ [:vfnmsub, :vfnmsub_fast]
        # fnmadd = vfnmsub_instr === :vfnmsub ? :vfnmadd : :vfnmadd_fast
        # fnmadd = :vfnmadd_fast
        DERIVATIVERULES[InstructionArgs(vfnmsub_instr, 3)] = DiffRule(
            Instruction[ vfnmsub_instr, :adjoint, :vnmul, :adjoint, :vnmul, :vsub ],
            deps, sections, returns, 3
        )
    end
end # let

DERIVATIVERULES[InstructionArgs(:exp, 1)] = DiffRule(
    Instruction[ :exp, :vmul ],
    [ [-1], [1,0] ],
    [ 1:1, 2:1, 2:2 ],
    [ 1, 2 ],
    1
)
DERIVATIVERULES[InstructionArgs(:expm1, 1)] = DiffRule(
    Instruction[ :expm1, :vfmadd ],
    [ [-1], [1], [2,0,0] ],
    [ 1:1, 2:1, 2:2 ],
    [ 1, 2 ],
    1
)
DERIVATIVERULES[InstructionArgs(:exp2, 1)] = DiffRule(
    Instruction[ :exp, :vmullog2, :vmul ],
    [ [-1], [1], [2,0] ],
    [ 1:1, 2:1, 2:3 ],
    [ 1, 3 ],
    1
)
DERIVATIVERULES[InstructionArgs(:exp10, 1)] = DiffRule(
    Instruction[ :exp, :vmullog10, :vmul ],
    [ [-1], [1], [2,0] ],
    [ 1:1, 2:1, 2:3 ],
    [ 1, 3 ],
    1
)

for sq ∈ [:sqrt, :vsqrt]
    DERIVATIVERULES[InstructionArgs(sq, 1)] = DiffRule(
        Instruction[ sq, :vadd, :vfdiv ],
        [ [-1], [1,1], [0,2] ],
        [ 1:1, 2:1, 2:3 ],
        [ 1, 3 ],
        1
    )
end
for cb ∈ [:cbrt, :cbrt_fast]
    DERIVATIVERULES[InstructionArgs(cb,1)] = DiffRule(
        Instruction[ cb, :vabs2, :vmul3, :vfdiv ],
        [ [-1], [1], [2], [0,3] ],
        [ 1:1, 2:1, 2:4 ],
        [ 1, 4 ],
        1
    )
end

for ln ∈ [:log, :log_fast, :vlog]
    DERIVATIVERULES[InstructionArgs(ln,1)] = DiffRule(
        Instruction[ ln, :vfdiv ],
        [ [-1], [0,-1] ],
        [ 1:1, 2:1, 2:2 ],
        [ 1, 2 ],
        1
    )
end
for (ln, comb) ∈ [(:log2,Instruction(:vmullog2)), (:log10,Instruction(:vmullog10))]
    DERIVATIVERULES[InstructionArgs(ln,1)] = DiffRule(
        Instruction[ ln, :vfdiv, comb ],
        [ [-1], [0,-1], [2] ],
        [ 1:1, 2:1, 2:3 ],
        [ 1, 3 ],
        1
    )
end

DERIVATIVERULES[InstructionArgs(:log1p,1)] = DiffRule(
    Instruction[ :log1p, :vadd1, :vfdiv ],
    [ [-1], [-1], [0, 2] ],
    [ 1:1, 2:1, 2:3 ],
    [ 1, 3 ],
    1
)

for p ∈ [:^, :pow]
    DERIVATIVERULES[InstructionArgs(p,2)] = DiffRule(
        Instruction[ :log, :vmul, :exp, :vfdiv, :vmul, :vmul ],
        [ [-2], [1,-1], [2], [-1,-2], [4,0], [3,0] ],
        [ 1:3, 4:3, 4:5, 6:6 ],
        [ 3, 5, 6 ],
        2
    )
end

DERIVATIVERULES[InstructionArgs(:sin,1)] = DiffRule(
    Instruction[ :sincos, :first, :last ],
    [ [-1], [1], [1] ],
    [ 1:3, 4:3, 4:3 ],
    [ 2, 3 ],
    1
)
DERIVATIVERULES[InstructionArgs(:cos,1)] = DiffRule(
    Instruction[ :sincos, :first, :last, :vsub ],
    [ [-1], [1], [1], [2] ],
    [ 1:3, 4:3, 4:4 ],
    [ 3, 4 ],
    1
)

DERIVATIVERULES[InstructionArgs(:tan,1)] = DiffRule(
    Instruction[ :tan, Instruction(:vfmaddaddone), :vmul ],
    [ [-1], [1], [2,0] ],
    [ 1:1, 2:1, 2:3 ],
    [ 1, 3 ],
    1
)
DERIVATIVERULES[InstructionArgs(:cot,1)] = DiffRule(
    Instruction[ :cot, Instruction(:vfmaddaddone), :vnmul ],
    [ [-1], [1], [2,0] ],
    [ 1:1, 2:1, 2:3 ],
    [ 1, 3 ],
    1
)

DERIVATIVERULES[InstructionArgs(Instruction(:DistributionParameters,:constrain),3)] = DiffRule(
    # Args are θ, descript
    # caller should calc gep(θ, offset)
    Instruction[ :constrain_pullback, :first, :last, :constrain_reverse! ],
    [ [-3, -2, -1], [1], [1], [0, 3, -1] ],
    [ 1:2, 3:2, 3:4, 5:4, 5:4 ],
    [ 2, 4, 0, 0 ],
    3
)

DERIVATIVERULES[InstructionArgs(:identity,1)] = DiffRule(
    Instruction[ :identity, :identity ],
    [ [-1], [0] ],
    [ 1:1, 2:1, 2:2 ],
    [ 1, 2 ],
    1
)

for f ∈ [:abs2, :vabs2]
    DERIVATIVERULES[InstructionArgs(f,1)] = DiffRule(
        Instruction[ f, :vadd, :vmul ],
        [ [-1], [-1,-1], [2,0] ],
        [ 1:1, 2:1, 2:3 ],
        [ 1, 3 ],
        1
    )
end
for at ∈ [:atan, :atan_fast]
    DERIVATIVERULES[InstructionArgs(at, 2)] = DiffRule(
        Instruction[ at, :vabs2, :vfmadd_fast, :vfdiv, :vnmul, :vmul ],
        [ [-2,-1], [ -2 ], [-1,-1,1], [0,3], [-2,4], [-1,4] ],
        [ 1:1, 2:4, 5:5, 6:6 ],
        [ 1, 5, 6 ],
        2
    )
end

@inline callunthunk(f::F, x) where {F} = f(unthunk(x))
const FALLBACK_RULES = Vector{DiffRule}(undef, 8);
let getters = [:second, :third, :fourth, :fifth, :sixth, :seventh, :eighth, :ninth]
    for i ∈ 1:8
        # Special case `rrule` in handling to add func?
        FALLBACK_RULES[i] = DiffRule(
            Instruction[:rrule, :first, :last, :callunthunk, (getters[1:i])...],
            [ collect(-i:-1), [1], [1], [3, 0], ([4] for _ ∈ 1:i)... ],
            [ 1:2, 3:4, (j:j for j ∈ 5:i+4)... ],
            [ 2, 4, (j for j ∈ 5:i+4)... ],
            i
        )
    end
end


text:
using UUIDs
is_ch_type(::Val{:LowCardinality})  = true
can_be_nullable(::Val{:LowCardinality}) = false

# Need to read additional keys.
# Additional keys are stored before indexes as value N and N keys
# after them.
const lc_has_additional_keys_bit = 1 << 9
# Need to update dictionary.
# It means that previous granule has different dictionary.
const lc_need_update_dictionary = 1 << 10

const lc_serialization_type = lc_has_additional_keys_bit | lc_need_update_dictionary

const lc_index_int_types = [:UInt8, :UInt16, :UInt32, :UInt64]

categorical_vector_type(::Type{Vector{T}}) where {T} =
                                                 CategoricalVector{T}
categorical_vector_type(::Type{CategoricalVector{T}}) where {T} =
                                                CategoricalVector{T}

result_type(::Val{:LowCardinality}, nested)  =
                             categorical_vector_type(result_type(nested))

function read_state_prefix(sock::ClickHouseSock, ::Val{:LowCardinality}, nested::TypeAst)
    ver = chread(sock, UInt64) # KeysSerializationVersion
    ver == 1 || error("unsupported LC serialization version: $(ver)")
    return ver
end

function write_state_prefix(sock::ClickHouseSock, ::Val{:LowCardinality}, nested::TypeAst)
    # KeysSerializationVersion. See ClickHouse docs.
    chwrite(sock, Int64(1))
end

function make_result(index::Vector{T}, keys, is_nullable) where {T}

    result = is_nullable ?
            CategoricalVector{Union{T, Missing}}(undef, 0, levels = index)  :
            CategoricalVector{T}(undef, 0, levels = index)
    result.refs = keys
    return result
end

function make_result(index::CategoricalVector{T}, keys, is_nullable) where {T}

    result = is_nullable ?
            CategoricalVector{Union{T, Missing}}(undef, 0, levels = get.(index))  :
            CategoricalVector{T}(undef, 0, levels = get.(index))
    result.refs = keys
    return result
end

function read_col_data(sock::ClickHouseSock, num_rows::VarUInt,
                                            ::Val{:LowCardinality}, nested::TypeAst)

    UInt64(num_rows)  == 0 && return read_col_data(sock, num_rows, nested)

    is_nested_nullable = (nested.name == :Nullable)
    notnullable_nested = is_nested_nullable ? nested.args[1] : nested

    serialization_type = chread(sock, UInt64)
    int_type = serialization_type & 0xf

    index_size = chread(sock, UInt64)
    index = read_col_data(sock, VarUInt(index_size), notnullable_nested)
    is_nested_nullable && (index = index[2:end])

    keys_size = chread(sock, UInt64)
    keys = read_col_data(sock, VarUInt(keys_size), Val(lc_index_int_types[int_type + 1]))

    (nested.name != :Nullable) && (keys .= keys .+ 1)

    return make_result(index, keys, nested.name == :Nullable)
end

unmissing_type(::Type{Union{Missing, T}}) where {T} = T
function write_col_data(sock::ClickHouseSock,
                                data::AbstractCategoricalVector{T},
                                ::Val{:LowCardinality}, nested::TypeAst) where {T}

    is_nested_nullable = (nested.name == :Nullable)
    notnullable_nested = is_nested_nullable ? nested.args[1] : nested

    isempty(data) && return

    int_type = floor(Int, log2(length(levels(data))) / 2)

    serialization_type = lc_serialization_type | int_type
    chwrite(sock, serialization_type)

    index = is_nested_nullable ?
                    vcat(missing_replacement(unmissing_type(T)), levels(data)) :
                    levels(data)

    chwrite(sock, length(index))
    write_col_data(sock, index, notnullable_nested)

    chwrite(sock, length(data))

    #In c++ indexes started from 0, in case of nullable nested 0 means null and
    # it's ok, but if nested not nullable we must sub 1 from index
    keys = is_nested_nullable ? data.refs : data.refs .- 1
    write_col_data(sock, keys, Val(lc_index_int_types[int_type + 1]))
end

function write_col_data(sock::ClickHouseSock,
                                data::AbstractVector{T},
                                v::Val{:LowCardinality}, nested::TypeAst) where {T}
    write_col_data(sock, CategoricalVector{T}(data), v, nested)
end


text:
__precompile__()

module Showoff

using Dates

export showoff

# suppress compile errors when there isn't a grisu_ccall macro
macro grisu_ccall(x, mode, ndigits)
    quote end
end

function grisu(v::AbstractFloat, mode, requested_digits)
    return tuple(Base.Grisu.grisu(v, mode, requested_digits)..., Base.Grisu.DIGITS)
end

# Fallback
function showoff(xs::AbstractArray, style=:none)
    result = Vector{String}(undef, length(xs))
    buf = IOBuffer()
    for (i, x) in enumerate(xs)
        show(buf, x)
        result[i] = String(take!(buf))
    end

    return result
end

# Floating-point

function concrete_minimum(xs)
    if isempty(xs)
        throw(ArgumentError("argument must not be empty"))
    end

    x_min = first(xs)
    for x in xs
        if isa(x, AbstractFloat) && isfinite(x)
            x_min = x
            break
        end
    end

    for x in xs
        if isa(x, AbstractFloat) && isfinite(x) && x < x_min
            x_min = x
        end
    end
    return x_min
end

function concrete_maximum(xs)
    if isempty(xs)
        throw(ArgumentError("argument must not be empty"))
    end

    x_max = first(xs)
    for x in xs
        if isa(x, AbstractFloat) && isfinite(x)
            x_max = x
            break
        end
    end

    for x in xs
        if isa(x, AbstractFloat) && isfinite(x) && x > x_max
            x_max = x
        end
    end
    return x_max
end

function plain_precision_heuristic(xs::AbstractArray{<:AbstractFloat})
    ys = filter(isfinite, xs)
    precision = 0
    for y in ys
        len, point, neg, digits = grisu(convert(Float32, y), Base.Grisu.SHORTEST, 0)
        precision = max(precision, len - point)
    end
    return max(precision, 0)
end

function scientific_precision_heuristic(xs::AbstractArray{<:AbstractFloat})
    ys = [x == 0.0 ? 0.0 : x / 10.0^floor(log10(abs(x)))
          for x in xs if isfinite(x)]
    return plain_precision_heuristic(ys) + 1
end

function showoff(xs::AbstractArray{<:AbstractFloat}, style=:auto)
    x_min = concrete_minimum(xs)
    x_max = concrete_maximum(xs)
    x_min = Float64(x_min)
    x_max = Float64(x_max)

    if !isfinite(x_min) || !isfinite(x_max)
        throw(ArgumentError("At least one finite value must be provided to formatter."))
    end

    if style == :auto
        if x_max != x_min && abs(log10(x_max - x_min)) > 4
            style = :scientific
        else
            style = :plain
        end
    end

    if style == :plain
        precision = plain_precision_heuristic(xs)
        return String[format_fixed(x, precision) for x in xs]
    elseif style == :scientific
        precision = scientific_precision_heuristic(xs)
        return String[format_fixed_scientific(x, precision, false)
                      for x in xs]
    elseif style == :engineering
        precision = scientific_precision_heuristic(xs)
        return String[format_fixed_scientific(x, precision, true)
                      for x in xs]
    else
        throw(ArgumentError("$(style) is not a recongnized number format"))
    end
end

# Print a floating point number at fixed precision. Pretty much equivalent to
# @sprintf("%0.$(precision)f", x), without the macro issues.
function format_fixed(x::AbstractFloat, precision::Integer)
    @assert precision >= 0

    if x == Inf
        return "∞"
    elseif x == -Inf
        return "-∞"
    elseif isnan(x)
        return "NaN"
    end

    len, point, neg, digits = grisu(x, Base.Grisu.FIXED, precision)

    buf = IOBuffer()
    if x < 0
        print(buf, '-')
    end

    for c in digits[1:min(point, len)]
        print(buf, convert(Char, c))
    end

    if point > len
        for _ in len:point-1
            print(buf, '0')
        end
    elseif point < len
        if point <= 0
            print(buf, '0')
        end
        print(buf, '.')
        if point < 0
            for _ in 1:-point
                print(buf, '0')
            end
            for c in digits[1:len]
                print(buf, convert(Char, c))
            end
        else
            for c in digits[point+1:len]
                print(buf, convert(Char, c))
            end
        end
    end

    trailing_zeros = precision - max(0, len - point)
    if trailing_zeros > 0 && point >= len
        print(buf, '.')
    end

    for _ in 1:trailing_zeros
        print(buf, '0')
    end

    String(take!(buf))
end

const superscript_numerals = ['⁰', '¹', '²', '³', '⁴', '⁵', '⁶', '⁷', '⁸', '⁹']

# Print a floating point number in scientific notation at fixed precision. Sort of equivalent
# to @sprintf("%0.$(precision)e", x), but prettier printing.
function format_fixed_scientific(x::AbstractFloat, precision::Integer,
                                 engineering::Bool)
    if x == 0.0
        return "0"
    elseif x == Inf
        return "∞"
    elseif x == -Inf
        return "-∞"
    elseif isnan(x)
        return "NaN"
    end

    mag = floor(Int, log10(abs(x)))
    if mag < 0
        grisu_precision = precision + abs(round(Int, mag))
    else
        grisu_precision = precision
    end

    len, point, neg, digits = grisu((x / 10.0^mag), Base.Grisu.FIXED, grisu_precision)
    point += mag

    @assert len > 0

    buf = IOBuffer()
    if x < 0
        print(buf, '-')
    end

    print(buf, convert(Char, digits[1]))
    nextdigit = 2
    if engineering
        while (point - 1) % 3 != 0
            if nextdigit <= len
                print(buf, convert(Char, digits[nextdigit]))
            else
                print(buf, '0')
            end
            nextdigit += 1
            point -= 1
        end
    end

    if precision > 1
        print(buf, '.')
    end

    for i in nextdigit:len
        print(buf, convert(Char, digits[i]))
    end

    for i in (len+1):precision
        print(buf, '0')
    end

    print(buf, "×10")
    for c in string(point - 1)
        if '0' <= c <= '9'
            print(buf, superscript_numerals[c - '0' + 1])
        elseif c == '-'
            print(buf, '⁻')
        end
    end

    return String(take!(buf))
end

function showoff(ds::AbstractArray{T}, style=:none) where T<:Union{Date,DateTime}
    years = Set()
    months = Set()
    days = Set()
    hours = Set()
    minutes = Set()
    seconds = Set()
    for d in ds
        push!(years, Dates.year(d))
        push!(months, Dates.month(d))
        push!(days, Dates.day(d))
        push!(hours, Dates.hour(d))
        push!(minutes, Dates.minute(d))
        push!(seconds, Dates.second(d))
    end
    all_same_year         = length(years)   == 1
    all_one_month         = length(months)  == 1 && 1 in months
    all_one_day           = length(days)    == 1 && 1 in days
    all_zero_hour         = length(hours)   == 1 && 0 in hours
    all_zero_minute       = length(minutes) == 1 && 0 in minutes
    all_zero_seconds      = length(minutes) == 1 && 0 in minutes
    all_zero_milliseconds = length(minutes) == 1 && 0 in minutes

    # first label format
    label_months = false
    label_days = false
    f1 = "u d, yyyy"
    f2 = ""
    if !all_zero_seconds
        f2 = "HH:MM:SS.sss"
    elseif !all_zero_seconds
        f2 = "HH:MM:SS"
    elseif !all_zero_hour || !all_zero_minute
        f2 = "HH:MM"
    else
        if !all_one_day
            first_label_format = "u d yyyy"
        elseif !all_one_month
            first_label_format = "u yyyy"
        elseif !all_one_day
            first_label_format = "yyyy"
        end
    end
    if f2 != ""
        first_label_format = string(f1, " ", f2)
    else
        first_label_format = f1
    end

    labels = Vector{String}(undef, length(ds))
    labels[1] = Dates.format(ds[1], first_label_format)
    d_last = ds[1]
    for (i, d) in enumerate(ds[2:end])
        if Dates.year(d) != Dates.year(d_last)
            if all_one_day && all_one_month
                f1 = "yyyy"
            elseif all_one_day && !all_one_month
                f1 = "u yyyy"
            else
                f1 = "u d, yyyy"
            end
        elseif Dates.month(d) != Dates.month(d_last)
            f1 = all_one_day ? "u" : "u d"
        elseif Dates.day(d) != Dates.day(d_last)
            f1 = "d"
        else
            f1 = ""
        end

        if f2 != ""
            f = string(f1, " ", f2)
        elseif f1 != ""
            f = f1
        else
            f = first_label_format
        end

        labels[i+1] = Dates.format(d, f)
        d_last = d
    end

    return labels
end

end # module


text:
"""
    BlockDescription

Description of blocks of `A` and `L` in a [`LinearMixedModel`](@ref)

## Fields
* `blknms`: Vector{String} of block names
* `blkrows`: Vector{Int} of the number of rows in each block
* `ALtypes`: Matrix{String} of datatypes for blocks in `A` and `L`.

When a block in `L` is the same type as the corresponding block in `A`, it is
described with a single name, such as `Dense`.  When the types differ the entry
in `ALtypes` is of the form `Diag/Dense`, as determined by a `shorttype` method.
"""
struct BlockDescription
    blknms::Vector{String}
    blkrows::Vector{Int}
    ALtypes::Matrix{String}
end
function BlockDescription(m::LinearMixedModel)
    A = m.A
    L = m.L
    blknms = push!(string.([fnames(m)...]), "fixed")
    k = length(blknms)
    ALtypes = fill("", k, k)
    Ltypes = fill(Nothing, k, k)
    for i in 1:k, j in 1:i
        ALtypes[i, j] = shorttype(A[Block(i,j)], L[Block(i,j)])
    end
    BlockDescription(
        blknms,
        [size(A[Block(i, 1)], 1) for i in 1:k],
        ALtypes,
     )
end

BlockDescription(m::GeneralizedLinearMixedModel) = BlockDescription(m.LMM)

shorttype(::UniformBlockDiagonal,::UniformBlockDiagonal) = "BlkDiag"
shorttype(::UniformBlockDiagonal,::Matrix) = "BlkDiag/Dense"
shorttype(::SparseMatrixCSC,::BlockedSparse) = "Sparse"
shorttype(::Diagonal,::Diagonal) = "Diagonal"
shorttype(::Diagonal,::Matrix) = "Diag/Dense"
shorttype(::Matrix,::Matrix) = "Dense"
shorttype(::SparseMatrixCSC,::SparseMatrixCSC) = "Sparse"
shorttype(::SparseMatrixCSC,::Matrix) = "Sparse/Dense"

function Base.show(io::IO, ::MIME"text/plain", b::BlockDescription)
    rowwidth = max(maximum(ndigits, b.blkrows) + 1, 5)
    colwidth = max(maximum(textwidth, b.blknms) + 1, 14)
    print(io, rpad("rows:", rowwidth))
    println(io, cpad.(b.blknms, colwidth)...)
    for (i, r) in enumerate(b.blkrows)
        print(io, lpad(string(r, ':'), rowwidth))
        for j in 1:i
            print(io, cpad(b.ALtypes[i, j], colwidth))
        end
        println(io)
    end
end
Base.show(io::IO, b::BlockDescription) = show(io, MIME"text/plain"(), b)

@deprecate describeblocks(io, m) show(io, BlockDescription(m))
@deprecate describeblocks(m) show(stdout, BlockDescription(m))


text:
using TomoForward
using Images
using Plots
using XfromProjections
using StaticArrays
using PyCall
using Logging

include("../tv_primaldual_flow.jl")
include("../optical_flow.jl")
include("../../phantoms/sperm_phantom.jl")

H,W = 609, 609
detmin, detmax = -38.0, 38.0
function radon_operator(height, width, detcount, angles)
    proj_geom = ProjGeom(0.5, detcount, angles)
    A = fp_op_parallel2d_line(proj_geom, height, width, detmin,detmax, detmin,detmax)
    return A
end
grid = collect(detmin:0.1:detmax)
r(s) = 1.0

images, tracks = get_sperm_phantom(301,r,grid)
frames = images[:,:,collect(271:274)]
@info size(frames)

w_tv = 0.01
w_flow  = 0.001
c=10.0

detcount = Int(floor(H*1.4))

cwd = @__DIR__
path = normpath(joinpath(@__DIR__, "result"))
!isdir(path) && mkdir(path)
cd(path)

#4 angles evenly spaced between 0 and pi
angles = [rand(0.0:0.001:π, 3), rand(0.0:0.001:π, 3), rand(0.0:0.001:π, 3), rand(0.0:0.001:π, 3)]#collect(range(π/2,3*π/2,length=ang_nr+1))[1:end-1]
@info angles
#match size of input image (generating data)
As = map(t -> radon_operator(size(frames[:,:,1],1),size(frames[:,:,1],2),detcount, angles[t]),1:size(frames)[3])
bs = zeros(size(As[1])[1],size(frames)[3])
map(t -> bs[:,t] = As[t]*vec(frames[:,:,t]), 1:size(frames)[3])

u0s = zeros(H,W,size(frames)[3])

As  = map(t -> radon_operator(H,W, detcount, angles[t]),1:size(u0s)[3])

@info "Reconstruction using tv regularization frame by frame"
niter=450
us_tv = zeros(H,W,size(frames)[3])
for t = 1:size(frames)[3]
    A = As[t]
    p = bs[:,t]
    u0 = u0s[:,:,t]
    us_tv[:,:,t] .= recon2d_tv_primaldual!(us_tv[:,:,t], A, p, niter, w_tv, c)
end

img2 = us_tv[:,:,2]#zeros(128,128)
img1 = us_tv[:,:,1]
H, W = size(img2)
#img2[30:60,40:50] .= 1.0

#flow = zeros(H, W, 2)
#flow[:,:, 1] .= -13.0
#flow[:,:, 2] .= 10.0
#flow[:,:, 1] +=rand(H,W)*4
#flow[:,:, 2] +=rand(H,W)*4

#Wop = compute_warping_operator(flow)
#img1_ = Wop * vec(img2)
#img1 = reshape(img1_, H, W)

us = zeros(H,W,2)
us[:,:,2] = img2
us[:,:,1] = img1

v = get_flows(us)

W_list = mapslices(f -> compute_warping_operator(f), v,dims=[1,2,3])
@info size(W_list)
#t = 1
#Wuv .= mul!(Wuv, W_list[t], vec(us[:,:,t+1])) .- vec(us[:,:,t])

#flow_x, flow_y, im2Warped = py"py_flow"(img1,img2)
#v = zeros(H,W,2)
#v[:,:,1] = flow_x
#v[:,:,2] = flow_y#W_list = mapslices(f -> compute_warping_operator(f), v, dims=[1,2,3])
#Woptical = W_list[1]#compute_warping_operator(v)
#img3_ = Woptical * vec(img2)
#img3 = reshape(img3_,H, W)
#W_list = mapslices(f -> compute_warping_operator(f), v,dims=[1,2,3])

t = 1
_ubar_1 = view(us,:,:,t+1)
_ubar = view(us,:,:,t)
Wu = W_list[t]*(collect(Iterators.flatten(_ubar_1))) - (collect(Iterators.flatten(_ubar)))
p3_ascent = reshape(Wu, H, W)

p3_adjoint = reshape(W_list[t]'*vec(p3_ascent),H,W)-p3_ascent

img3 = max.(us[:,:,1].-p3_adjoint, 0.0)

p1 = plot(Gray.(img1[:,:,1]), aspect_ratio=:equal, framestyle=:none, title="img1")
p2 = plot(Gray.(img2[:,:,1]), aspect_ratio=:equal, framestyle=:none, title="img2")
p3 = plot(Gray.(img3), aspect_ratio=:equal, framestyle=:none, title="rec")
p4 = plot(Gray.(p3_ascent), aspect_ratio=:equal, framestyle=:none, title="ascent")
p5 = plot(Gray.(p3_adjoint), aspect_ratio=:equal, framestyle=:none, title="adjoint")
l = @layout [a b c d e]
plot(p1, p2, p3, p4, p5, layout = l, size=(2000,600), linewidth=5)
#XX, YY = repeat(collect(1:W)', H, 1), repeat(collect(1:H), 1, W)
#p1 = quiver([(XX...)...][1:20:end], [(YY...)...][1:20:end], quiver=([(v[:,:,1]...)...][1:20:end]*0.001, [(v[:,:,2]...)...][1:20:end]*0.001), title="v")
#p2 = quiver([(XX...)...][1:20:end], [(YY...)...][1:20:end], quiver=([(flow[:,:,1]...)...][1:20:end]*0.1, [(flow[:,:,2]...)...][1:20:end]*0.1), title="flow")
#l = @layout [a b]
#plot(p1, p2, layout = l, size=(2000,600), linewidth=5)

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
using DataStructures

# Returns all neighboring positions of 'pos'.
function getNbors(M::Matrix{Int}, pos::CartesianIndex{2}) :: Vector{CartesianIndex{2}}
  s = size(M)
  nbors :: Vector{CartesianIndex{2}} = []
  (pos[1] != 1   ) && push!(nbors, pos + CartesianIndex(-1,0))
  (pos[1] != s[1]) && push!(nbors, pos + CartesianIndex(1,0))
  (pos[2] != 1   ) && push!(nbors, pos + CartesianIndex(0,-1))
  (pos[2] != s[2]) && push!(nbors, pos + CartesianIndex(0,1))
  return nbors
end

# Calculates all low points for a given height map.
function lowPts(M::Matrix{Int}) :: Vector{CartesianIndex{2}}
  s = size(M)
  lpts :: Vector{CartesianIndex{2}} = []
  for i in 1:s[1]
    for j in 1:s[2]
      pos = CartesianIndex(i,j)
      nbors = getNbors(M,pos)
      (all(map( n -> M[n] > M[pos], nbors))) && push!(lpts,pos)
    end
  end
  return lpts
end

# Computes risk level of a position.
riskLevel(M::Matrix{Int}, pos::CartesianIndex{2}) :: Int = M[pos] + 1
riskLevel(M::Matrix{Int}, pos::Vector{CartesianIndex{2}}) :: Int = sum(map(p -> riskLevel(M,p), pos))

# Calculates all basins, given the set of low points.
function basins(M::Matrix{Int}, lpts::Vector{CartesianIndex{2}}) :: Vector{Int}
  visited :: Matrix{Int} = -ones(size(M))

  nBasin = 1
  sizes :: Vector{Int} = []
  for p in lpts
    if visited[p] == -1 # If new bassin
      append!(sizes,0)
      q :: Queue{CartesianIndex{2}} = Queue{CartesianIndex{2}}()
      enqueue!(q,p)
      while !isempty(q)
        next = dequeue!(q)
        if visited[next] == -1 && M[next] != 9
          visited[next] = nBasin
          sizes[nBasin] += 1
          map(n -> enqueue!(q,n), getNbors(M,next))
        end
      end
      nBasin += 1
    end
  end

  return sort(sizes,rev=true)
end

function readInput() :: Matrix{Int}
  f = open("./inputs/day9.txt")
  data :: Matrix{Int} = zeros(100,100)

  i = 1
  for line in eachline(f)
    data[i,:] = map(x -> parse(Int,x), collect(line))
    i += 1
  end
  return data
end

M = readInput()

println("Problem 1: ")
lpts = lowPts(M)
res = riskLevel(M,lpts)
println(res)

sizes = basins(M,lpts)
println("Problem 2: ")
println(sizes)
println(sizes[1]*sizes[2]*sizes[3])


text:
# Benchmark some basic operations of `Table`

# Choice of benchmarks drawn from: https://blog.thedataincubator.com/2018/01/pandas-vs-postgresql/

using MinimumViableTables
using BenchmarkTools
using SplitApplyCombine
using DataFrames

include("create_dataset.jl")

function printtime(t)
    if t < 10^3
        return "$t ns"
    elseif t < 10^6
        t = 1e-3 * t
        return "$t μs"
    elseif t < 10^9
        t = 1e-6 * t
        return "$t ms"
    else
        t = 1e-9 * t
        return "$t s"
    end
end

ns = [10, 100, 1000, 10_000, 100_000, 1_000_000, 10_000_000]
n_ns = length(ns)

select_times = zeros(n_ns)
select_times_df = zeros(n_ns)

filter_times = zeros(n_ns)
filter_times_sorted = zeros(n_ns)
filter_times_primary_sorted = zeros(n_ns)
filter_times_hashed = zeros(n_ns)
filter_times_df = zeros(n_ns)

groupreduce_times = zeros(n_ns)
groupreduce_times_df = zeros(n_ns)

join_times_unique_sorted = zeros(n_ns)
join_times_unique_hashed = zeros(n_ns)
join_times_sorted = zeros(n_ns)
join_times_hashed = zeros(n_ns)
join_times_merge = zeros(n_ns)
join_times_merge_unique1 = zeros(n_ns)
join_times_merge_unique2 = zeros(n_ns)
join_times_merge_unique12 = zeros(n_ns)
join_times_df = zeros(n_ns)

index_times_sort_section = zeros(n_ns)
index_times_uniquesort_section = zeros(n_ns)
index_times_hash_section = zeros(n_ns)
index_times_uniquehash_section = zeros(n_ns)

index_times_sort_id = zeros(n_ns)
index_times_uniquesort_id = zeros(n_ns)
index_times_sort_id2 = zeros(n_ns)
index_times_uniquesort_id2 = zeros(n_ns)
index_times_hash_id2 = zeros(n_ns)
index_times_uniquehash_id2 = zeros(n_ns)

function select_column(t)
    Project{(:score_1,)}()(t)
end

function grouping_df(df_in)
    by(df_in, :section) do df
         DataFrame(mean_1 = mean(df[:score_1]), max_2 = maximum(df[:score_2]))
    end
end

for (i, n) in enumerate(ns)
    println("n = $n")
    (t_a, t_b) = create_dataset(n)
    df_a = DataFrame(; columns(t_a)...)
    df_b = DataFrame(; columns(t_b)...)
    rename!(df_b, :id2 => :id)

    index_times_sort_section[i] = median((@benchmark accelerate($t_a, $(SortIndex{(:section,)}))).times) * 1e-9
    println("Sort by section: t = $(printtime(index_times_sort_section[i]*1e9))")
    index_times_hash_section[i] = median(@benchmark(accelerate($t_a, $(HashIndex{(:section,)}))).times) * 1e-9
    println("Hash by section: t = $(printtime(index_times_hash_section[i]*1e9))")

    index_times_sort_id[i] = median(@benchmark(accelerate($t_a, $(SortIndex{(:id,)}))).times) * 1e-9
    println("Sort by id: t = $(printtime(index_times_sort_id[i]*1e9))")
    index_times_uniquesort_id[i] = median(@benchmark(accelerate($t_a, $(UniqueSortIndex{(:id,)}))).times) * 1e-9
    println("Unique sort by id: t = $(printtime(index_times_uniquesort_id[i]*1e9))")
    index_times_sort_id2[i] = median(@benchmark(accelerate($t_b, $(SortIndex{(:id2,)}))).times) * 1e-9
    println("Sort by id2: t = $(printtime(index_times_sort_id2[i]*1e9))")
    index_times_uniquesort_id2[i] = median(@benchmark(accelerate($t_b, $(UniqueSortIndex{(:id2,)}))).times) * 1e-9
    println("Unique sort by id2: t = $(printtime(index_times_uniquesort_id2[i]*1e9))")
    index_times_hash_id2[i] = median((@benchmark accelerate($t_b, $(HashIndex{(:id2,)}))).times) * 1e-9
    println("Hash by id2: t = $(printtime(index_times_hash_id2[i]*1e9))")
    index_times_uniquehash_id2[i] = median(@benchmark(accelerate($t_b, $(UniqueHashIndex{(:id2,)}))).times) * 1e-9
    println("Unique hash by id2: t = $(printtime(index_times_uniquehash_id2[i]*1e9))")

    bm = @benchmark select_column($t_a)
    time_ns = median(bm.times)
    select_times[i] = time_ns * 1e-9
    println("Select column: t = $(printtime(time_ns))")

    global bm = @benchmark getindex($(df_a), $[:score_1])
    time_ns = median(bm.times)
    select_times_df[i] = time_ns * 1e-9
    println("Select column (DataFrame): t = $(printtime(time_ns))")

    perm = sortperm(columns(t_a).section)
    t_a_sorted = Table(columns(t_a), (SortIndex{(:section,)}(perm),))
    tmp = t_a[perm]
    t_a_primary_sorted = Table(columns(tmp), (SortIndex{(:section,)}(Base.OneTo(n)),))
    t_a_hashed = accelerate(t_a, HashIndex{(:section,)})

    n = length(filter(IsEqual(section = 'A'), t_a))
    global bm = @benchmark filter($(IsEqual(section = 'A')), $t_a)
    time_ns = median(bm.times)
    filter_times[i] = time_ns * 1e-9
    println("Filter column: t = $(printtime(time_ns)) ($n rows)")

    n = length(filter(IsEqual(section = 'A'), t_a_sorted))
    global bm = @benchmark filter($(IsEqual(section = 'A')), $t_a_sorted)
    time_ns = median(bm.times)
    filter_times_sorted[i] = time_ns * 1e-9
    println("Filter column (sort index): t = $(printtime(time_ns)) ($n rows)")

    n = length(filter(IsEqual(section = 'A'), t_a_primary_sorted))
    global bm = @benchmark filter($(IsEqual(section = 'A')), $t_a_primary_sorted)
    time_ns = median(bm.times)
    filter_times_primary_sorted[i] = time_ns * 1e-9
    println("Filter column (primary sort index): t = $(printtime(time_ns)) ($n rows)")

    n = length(filter(IsEqual(section = 'A'), t_a_hashed))
    global bm = @benchmark filter($(IsEqual(section = 'A')), $t_a_hashed)
    time_ns = median(bm.times)
    filter_times_hashed[i] = time_ns * 1e-9
    println("Filter column (hash index): t = $(printtime(time_ns)) ($n rows)")

    n = size(filter(row -> isequal(row[:section], 'A'), df_a), 1)
    global bm = @benchmark filter($(row -> isequal(row[:section], 'A')), $df_a)
    time_ns = median(bm.times)
    filter_times[i] = time_ns * 1e-9
    println("Filter column (DataFrame): t = $(printtime(time_ns)) ($n rows)")

    n = length(groupreduce(Project{(:section,)}(), identity, (cum, row) -> (cum[1] + row.score_1, max(cum[2], row.score_2)), (0.0, 0.0), Project{(:section, :score_1, :score_2)}()(t_a)))
    bm = @benchmark groupreduce($(Project{(:section,)}()), $identity, $((cum, row) -> (cum[1] + row.score_1, max(cum[2], row.score_2))), $((0.0, 0.0)), $(Project{(:section, :score_1, :score_2)}()(t_a)))
    time_ns = median(bm.times)
    groupreduce_times[i] = time_ns * 1e-9
    println("Group reduction: t = $(printtime(time_ns)) ($n groups)")

    n = size(grouping_df(df_a), 1)
    global bm = @benchmark grouping_df($df_a)
    time_ns = median(bm.times)
    groupreduce_times_df[i] = time_ns * 1e-9
    println("Group reduction (DataFrame): t = $(printtime(time_ns)) ($n groups)")

    t_b_sorted = accelerate(t_b, SortIndex{(:id2,)})
    t_product_sorted = ProductTable(t_a, t_b_sorted)
    n = length(filter(Equals{(:id, :id2)}(), t_product_sorted))
    bm = @benchmark filter($(Equals{(:id, :id2)}()), $t_product_sorted)
    time_ns = median(bm.times)
    join_times_sorted[i] = time_ns * 1e-9
    println("Join (sort): t = $(printtime(time_ns)) ($n rows)")

    t_b_unique_sorted = accelerate(t_b, UniqueSortIndex{(:id2,)})
    t_product_unique_sorted = ProductTable(t_a, t_b_unique_sorted)
    n = length(filter(Equals{(:id, :id2)}(), t_product_unique_sorted))
    global bm = @benchmark filter($(Equals{(:id, :id2)}()), $t_product_unique_sorted)
    time_ns = median(bm.times)
    join_times_unique_sorted[i] = time_ns * 1e-9
    println("Join (unique sort): t = $(printtime(time_ns)) ($n rows)")

    t_b_hashed = accelerate(t_b, HashIndex{(:id2,)})
    t_product_hashed = ProductTable(t_a, t_b_hashed)
    n = length(filter(Equals{(:id, :id2)}(), t_product_hashed))
    global bm = @benchmark filter($(Equals{(:id, :id2)}()), $t_product_hashed)
    time_ns = median(bm.times)
    join_times_hashed[i] = time_ns * 1e-9
    println("Join (hash): t = $(printtime(time_ns)) ($n rows)")

    t_b_unique_hashed = accelerate(t_b, UniqueHashIndex{(:id2,)})
    t_product_unique_hashed = ProductTable(t_a, t_b_unique_hashed)
    n = length(filter(Equals{(:id, :id2)}(), t_product_unique_hashed))
    global bm = @benchmark filter($(Equals{(:id, :id2)}()), $t_product_unique_hashed)
    time_ns = median(bm.times)
    join_times_unique_hashed[i] = time_ns * 1e-9
    println("Join (unique hash): t = $(printtime(time_ns)) ($n rows)")

    t_a_merge = accelerate(t_a, SortIndex{(:id,)})
    t_product_merge = ProductTable(t_a_merge, t_b_sorted)
    n = length(filter(Equals{(:id, :id2)}(), t_product_merge))
    global bm = @benchmark filter($(Equals{(:id, :id2)}()), $t_product_merge)
    time_ns = median(bm.times)
    join_times_merge[i] = time_ns * 1e-9
    println("Join (merge): t = $(printtime(time_ns)) ($n rows)")

    t_a_merge_unique = accelerate(t_a, UniqueSortIndex{(:id,)})
    t_product_merge_unique1 = ProductTable(t_a_merge_unique, t_b_sorted)
    n = length(filter(Equals{(:id, :id2)}(), t_product_merge_unique1))
    global bm = @benchmark filter($(Equals{(:id, :id2)}()), $t_product_merge_unique1)
    time_ns = median(bm.times)
    join_times_merge_unique1[i] = time_ns * 1e-9
    println("Join (merge unique/non-unique): t = $(printtime(time_ns)) ($n rows)")

    t_product_merge_unique2 = ProductTable(t_a_merge, t_b_unique_sorted)
    n = length(filter(Equals{(:id, :id2)}(), t_product_merge_unique2))
    global bm = @benchmark filter($(Equals{(:id, :id2)}()), $t_product_merge_unique2)
    time_ns = median(bm.times)
    join_times_merge_unique2[i] = time_ns * 1e-9
    println("Join (merge non-unique/unique): t = $(printtime(time_ns)) ($n rows)")

    t_product_merge_unique12 = ProductTable(t_a_merge_unique, t_b_unique_sorted)
    n = length(filter(Equals{(:id, :id2)}(), t_product_merge_unique12))
    global bm = @benchmark filter($(Equals{(:id, :id2)}()), $t_product_merge_unique12)
    time_ns = median(bm.times)
    join_times_merge_unique12[i] = time_ns * 1e-9
    println("Join (merge unique/unique): t = $(printtime(time_ns)) ($n rows)")

    n = size(join(df_a, df_b, on = :id), 1)
    global bm = @benchmark join($df_a, $df_b, on = :id)
    time_ns = median(bm.times)
    join_times_df[i] = time_ns * 1e-9
    println("Join (DataFrame): t = $(printtime(time_ns)) ($n rows)")

    println()
end

#using UnicodePlots

#plot = UnicodePlots.scatterplot(log10.(ns), log10.(filter_times))
#UnicodePlots.scatterplot!(plot, log10.(ns), log10.(filter_times_sorted))
#UnicodePlots.scatterplot!(plot, log10.(ns), log10.(filter_times_hashed))


text:
module Fields

import ..slab
import ..column
import ..DataLayouts: DataLayouts, AbstractData, DataStyle
import ..Spaces: Spaces, AbstractSpace
import ..Geometry: Geometry, Cartesian12Vector
import ..RecursiveApply
import ..Topologies

import LinearAlgebra

"""
    Field(values, space)

A set of `values` defined at each point of a `space`.
"""
struct Field{V <: AbstractData, S <: AbstractSpace}
    values::V
    space::S
    # add metadata/attributes?
    function Field{V, S}(values::V, space::S) where {V, S}
        #TODOneed to enforce that the data size matches the space
        return new{V, S}(values, space)
    end
end
Field(values::V, space::S) where {V <: AbstractData, S <: AbstractSpace} =
    Field{V, S}(values, space)

Field(::Type{T}, space::S) where {T, S <: AbstractSpace} =
    Field(similar(Spaces.coordinates_data(space), T), space)

# Spectral Element Field
const SpectralElementField2D{V, S} =
    Field{V, S} where {V <: AbstractData, S <: Spaces.SpectralElementSpace2D}

const SpectralElementField1D{V, S} =
    Field{V, S} where {V <: AbstractData, S <: Spaces.SpectralElementSpace1D}

# Finite Difference Field
const FiniteDifferenceField{V, S} =
    Field{V, S} where {V <: AbstractData, S <: Spaces.FiniteDifferenceSpace}

const FaceFiniteDifferenceField{V, S} =
    Field{V, S} where {V <: AbstractData, S <: Spaces.FaceFiniteDifferenceSpace}

const CenterFiniteDifferenceField{V, S} = Field{
    V,
    S,
} where {V <: AbstractData, S <: Spaces.CenterFiniteDifferenceSpace}

# Extruded Fields
const ExtrudedFiniteDifferenceField{V, S} = Field{
    V,
    S,
} where {V <: AbstractData, S <: Spaces.ExtrudedFiniteDifferenceSpace}

const FaceExtrudedFiniteDifferenceField{V, S} = Field{
    V,
    S,
} where {V <: AbstractData, S <: Spaces.FaceExtrudedFiniteDifferenceSpace}

const CenterExtrudedFiniteDifferenceField{V, S} = Field{
    V,
    S,
} where {V <: AbstractData, S <: Spaces.CenterExtrudedFiniteDifferenceSpace}

Base.propertynames(field::Field) = propertynames(getfield(field, :values))
@inline field_values(field::Field) = getfield(field, :values)

# Define the axes field to be the todata(bc) of the return field
@inline Base.axes(field::Field) = getfield(field, :space)

# need to define twice to avoid ambiguities
@inline Base.getproperty(field::Field, name::Symbol) = Field(
    DataLayouts._getproperty(field_values(field), Val{name}()),
    axes(field),
)

@inline Base.getproperty(field::Field, name::Integer) =
    Field(getproperty(field_values(field), name), axes(field))

Base.eltype(field::Field) = eltype(field_values(field))
Base.parent(field::Field) = parent(field_values(field))

# to play nice with DifferentialEquations; may want to revisit this
# https://github.com/SciML/SciMLBase.jl/blob/697bd0c0c7365e77fa311f2d32eade70f43a8d50/src/solutions/ode_solutions.jl#L31
Base.size(field::Field) = ()
Base.length(field::Fields.Field) = 1

slab(field::Field, inds...) =
    Field(slab(field_values(field), inds...), slab(axes(field), inds...))

column(field::Field, inds...) =
    Field(column(field_values(field), inds...), column(axes(field), inds...))

const SlabField{V, S} =
    Field{V, S} where {V <: AbstractData, S <: Spaces.SpectralElementSpaceSlab}

const SlabField1D{V, S} = Field{
    V,
    S,
} where {V <: AbstractData, S <: Spaces.SpectralElementSpaceSlab1D}

const SlabField2D{V, S} = Field{
    V,
    S,
} where {V <: AbstractData, S <: Spaces.SpectralElementSpaceSlab2D}

Topologies.nlocalelems(field::Field) = Topologies.nlocalelems(axes(field))

# nice printing
# follow x-array like printing?
# repl: #https://earth-env-data-science.github.io/lectures/xarray/xarray.html
# html: https://unidata.github.io/MetPy/latest/tutorials/xarray_tutorial.html
function Base.show(io::IO, field::Field)
    print(io, eltype(field), "-valued Field:")
    _show_compact_field(io, field, "  ", true)
    # print(io, "\non ", axes(field)) # TODO: write a better space print
end
function _show_compact_field(io, field, prefix, isfirst = false)
    #print(io, prefix1)
    if eltype(field) <: Number
        if isfirst
            print(io, "\n", prefix)
        end
        print(
            IOContext(io, :compact => true, :limit => true),
            vec(parent(field)),
        )
    else
        names = propertynames(field)
        for name in names
            subfield = getproperty(field, name)
            if sizeof(eltype(subfield)) == 0
                continue
            end
            print(io, "\n", prefix)
            print(io, name, ": ")
            _show_compact_field(io, getproperty(field, name), prefix * "  ")
        end
    end
end

# https://github.com/gridap/Gridap.jl/blob/master/src/Fields/DiffOperators.jl#L5
# https://github.com/gridap/Gridap.jl/blob/master/src/Fields/FieldsInterfaces.jl#L70

Base.similar(field::Field, ::Type{Eltype}) where {Eltype} =
    Field(Eltype, axes(field))
Base.similar(field::Field) = Field(similar(field_values(field)), axes(field))

# fields on different spaces
function Base.similar(field::Field, space_to::AbstractSpace)
    similar(field, space_to, eltype(field))
end
function Base.similar(
    field::Field,
    space_to::AbstractSpace,
    ::Type{Eltype},
) where {Eltype}
    Field(Eltype, space_to)
end

Base.copy(field::Field) = Field(copy(field_values(field)), axes(field))

function Base.copyto!(dest::Field{V, M}, src::Field{V, M}) where {V, M}
    @assert axes(dest) == axes(src)
    copyto!(field_values(dest), field_values(src))
    return dest
end

function Base.zeros(::Type{FT}, space::AbstractSpace) where {FT}
    field = Field(FT, space)
    data = parent(field)
    fill!(data, zero(eltype(data)))
    return field
end
Base.zeros(space::AbstractSpace) = zeros(Spaces.undertype(space), space)

function Base.ones(::Type{FT}, space::AbstractSpace) where {FT}
    field = Field(FT, space)
    data = parent(field)
    fill!(data, one(eltype(data)))
    return field
end
Base.ones(space::AbstractSpace) = ones(Spaces.undertype(space), space)

function Base.zero(field::Field)
    zfield = similar(field)
    zarray = parent(zfield)
    fill!(zarray, zero(eltype(zarray)))
    return zfield
end

"""
    coordinate_field(space::AbstractSpace)

Construct a `Field` of the coordinates of the space.
"""
coordinate_field(space::AbstractSpace) =
    Field(Spaces.coordinates_data(space), space)
coordinate_field(field::Field) = coordinate_field(axes(field))

"""
    local_geometry_field(space::AbstractSpace)

Construct a `Field` of the `LocalGeometry` of the space.
"""
local_geometry_field(space::AbstractSpace) =
    Field(Spaces.local_geometry_data(space), space)
local_geometry_field(field::Field) = local_geometry_field(axes(field))

include("broadcast.jl")
include("mapreduce.jl")
include("compat_diffeq.jl")
include("fieldvector.jl")

function interpcoord(elemrange, x::Real)
    n = length(elemrange) - 1
    z = x == elemrange[end] ? n : searchsortedlast(elemrange, x) # element index
    @assert 1 <= z <= n
    lo = elemrange[z]
    hi = elemrange[z + 1]
    # Find ξ ∈ [-1,1] such that
    # x = (1-ξ)/2 * lo + (1+ξ)/2 * hi
    #   = (lo + hi) / 2 + ξ * (hi - lo) / 2
    ξ = (2x - (lo + hi)) / (hi - lo)
    return z, ξ
end

"""
    Spaces.variational_solve!(field)

Divide `field` by the mass matrix.
"""
function Spaces.variational_solve!(field::Field)
    Spaces.variational_solve!(field_values(field), axes(field))
    return field
end

function Spaces.horizontal_dss!(field::Field)
    Spaces.horizontal_dss!(field_values(field), axes(field))
    return field
end
function Spaces.weighted_dss!(field::Field)
    Spaces.weighted_dss!(field_values(field), axes(field))
    return field
end

end # module


text:
type InvalidFormatError <: Exception
    msg
    InvalidFormatError(msg = "") = new(msg)
end

type NoDataException <: Exception end


text:
@testset "Constraints" begin
    T = 5
    num_state = 2
    num_action = 1
    num_parameter = 0
    dim_x = [num_state for t = 1:T]
    dim_u = [num_action for t = 1:T-1]
    dim_w = [num_parameter for t = 1:T]
    x = [rand(dim_x[t]) for t = 1:T]
    u = [[rand(dim_u[t]) for t = 1:T-1]..., zeros(0)]
    w = [rand(dim_w[t]) for t = 1:T]

    ct = (x, u, w) -> [-ones(num_state) - x; x - ones(num_state)]
    cT = (x, u, w) -> x

    cont = Constraint(ct, num_state, num_action, indices_inequality=collect(1:2num_state), num_parameter=num_parameter)
    conT = Constraint(cT, num_state, 0, indices_inequality=collect(1:num_state), num_parameter=num_parameter)

    constraints = [[cont for t = 1:T-1]..., conT]

    nct = 2 * num_state
    ncT = num_state
    ct0 = zeros(nct)
    cT0 = zeros(ncT)
    cont.evaluate(ct0, x[1], u[1], w[1])
    conT.evaluate(cT0, x[T], u[T], w[T])
    @test norm(ct0 - [-ones(num_state) - x[1]; x[1] - ones(num_state)]) < 1.0e-8
    @test norm(cT0 - x[T]) < 1.0e-8

    cc = [[zeros(nct) for t = 1:T-1]..., zeros(ncT)]
    IterativeLQR.constraint!(cc, constraints, x, u, w)

    @test norm(vcat(cc...) - vcat([ct(x[t], u[t], w[t]) for t = 1:T-1]..., cT(x[T], u[T], w[T]))) < 1.0e-8

    jx = [[zeros(nct, num_state) for t = 1:T-1]..., zeros(ncT, num_state)]
    ju = [[zeros(nct, num_action) for t = 1:T-1]..., zeros(ncT, num_action)]
    IterativeLQR.jacobian!(jx, ju, constraints, x, u, w)

    for t = 1:T-1
        @test norm(jx[t] - ForwardDiff.jacobian(x -> ct(x, u[t], w[t]), x[t])) < 1.0e-8
        @test norm(ju[t] - ForwardDiff.jacobian(u -> ct(x[t], u, w[t]), u[t])) < 1.0e-8
    end
    @test norm(jx[T] - ForwardDiff.jacobian(x -> cT(x, u[T], w[T]), x[T])) < 1.0e-8
end


text:
"""
    FullyObservablePOMDP(mdp)

Turn `MDP` `mdp` into a `POMDP` where the observations are the states of the MDP.
"""
struct FullyObservablePOMDP{M,S,A} <: POMDP{S,A,S}
    mdp::M
end

function FullyObservablePOMDP(m::MDP)
    return FullyObservablePOMDP{typeof(m), statetype(m), actiontype(m)}(m)
end

mdptype(::Type{FullyObservablePOMDP{M,S,A}}) where {M,S,A} = M

function POMDPs.DDNStructure(::Type{M}) where M <: FullyObservablePOMDP
    MM = mdptype(M)
    add_obsnode(DDNStructure(MM))
end

add_obsnode(ddn) = add_node(ddn, :o, FunctionDDNNode((m,sp)->sp), (:sp,)) # for ::DDNStructure, but this is not declared yet POMDPs in v0.7.3

POMDPs.observations(pomdp::FullyObservablePOMDP) = states(pomdp.mdp)
POMDPs.obsindex(pomdp::FullyObservablePOMDP{S, A}, o::S) where {S, A} = stateindex(pomdp.mdp, o)

POMDPs.convert_o(T::Type{V}, o, pomdp::FullyObservablePOMDP) where {V<:AbstractArray} = convert_s(T, s, pomdp.mdp)
POMDPs.convert_o(T::Type{S}, vec::V, pomdp::FullyObservablePOMDP) where {S,V<:AbstractArray} = convert_s(T, vec, pomdp.mdp)

POMDPs.gen(::DDNNode{:o}, m::FullyObservablePOMDP, sp, rng) = sp

function POMDPs.observation(pomdp::FullyObservablePOMDP, a, sp)
    return Deterministic(sp)
end

function POMDPs.observation(pomdp::FullyObservablePOMDP, s, a, sp)
    return Deterministic(sp)
end

# inherit other function from the MDP type

POMDPs.states(pomdp::FullyObservablePOMDP) = states(pomdp.mdp)
POMDPs.actions(pomdp::FullyObservablePOMDP) = actions(pomdp.mdp)
POMDPs.transition(pomdp::FullyObservablePOMDP, s, a) = transition(pomdp.mdp, s, a)
POMDPs.initialstate_distribution(pomdp::FullyObservablePOMDP) = initialstate_distribution(pomdp.mdp)
POMDPs.initialstate(pomdp::FullyObservablePOMDP, rng::AbstractRNG) = initialstate(pomdp.mdp, rng)
POMDPs.isterminal(pomdp::FullyObservablePOMDP, s) = isterminal(pomdp.mdp, s)
POMDPs.discount(pomdp::FullyObservablePOMDP) = discount(pomdp.mdp)
POMDPs.stateindex(pomdp::FullyObservablePOMDP, s) = stateindex(pomdp.mdp, s)
POMDPs.actionindex(pomdp::FullyObservablePOMDP, a) = actionindex(pomdp.mdp, a)
POMDPs.convert_s(T::Type{V}, s, pomdp::FullyObservablePOMDP) where V<:AbstractArray = convert_s(T, s, pomdp.mdp)
POMDPs.convert_s(T::Type{S}, vec::V, pomdp::FullyObservablePOMDP) where {S,V<:AbstractArray} = convert_s(T, vec, pomdp.mdp)
POMDPs.convert_a(T::Type{V}, a, pomdp::FullyObservablePOMDP) where V<:AbstractArray = convert_a(T, a, pomdp.mdp)
POMDPs.convert_a(T::Type{A}, vec::V, pomdp::FullyObservablePOMDP) where {A,V<:AbstractArray} = convert_a(T, vec, pomdp.mdp)

POMDPs.gen(d::DDNNode, m::FullyObservablePOMDP, args...) = gen(d, m.mdp, args...)
POMDPs.gen(m::FullyObservablePOMDP, s, a, rng) = gen(m.mdp, s, a, rng)
POMDPs.reward(pomdp::FullyObservablePOMDP, s, a) = reward(pomdp.mdp, s, a)

# deprecated in POMDPs v0.8
add_obsnode(ddn::POMDPs.DDNStructureV7{(:s,:a,:sp,:r)}) = POMDPs.DDNStructureV7{(:s,:a,:sp,:o,:r)}()
add_obsnode(ddn::POMDPs.DDNStructureV7) = error("FullyObservablePOMDP only supports MDPs with the standard DDN Structure (DDNStructureV7{(:s,:a,:sp,:r)}) with POMDPs v0.7.")

POMDPs.generate_s(pomdp::FullyObservablePOMDP, s, a, rng::AbstractRNG) = generate_s(pomdp.mdp, s, a, rng)
POMDPs.generate_sr(pomdp::FullyObservablePOMDP, s, a, rng::AbstractRNG) = generate_sr(pomdp.mdp, s, a, rng)
POMDPs.n_actions(pomdp::FullyObservablePOMDP) = n_actions(pomdp.mdp)
POMDPs.n_states(pomdp::FullyObservablePOMDP) = n_states(pomdp.mdp)
function POMDPs.generate_o(pomdp::FullyObservablePOMDP, s, rng::AbstractRNG)
    return s
end


text:
export generate_noise_mpi

import CorrNoise
using Random
using FITSIO

try
    import MPI
catch
end

"""
This structure holds information about a certain chunk of data produced by a certain MPI process to send to another MPI process.

Field         | Type  | Meaning
:------------ |:----- |:--------------------------------------------
`pol_number`  | Int   | ID number of the polarimeter
`first_idx`   | Int   | Index of the first element of the chunk
`last_idx`    | Int   | Index of the last element of the chunk
`dest_rank`   | Int   | Number of destination rank

# Example
Suppose that the following array of data: [1, 10, 20, 3, 4, 5, 7]
is measured by the polarimeter number 2.
If the chunk [20, 3, 4] needs to be sent to rank number 3, the corresponding structure will be:

chunk_to_send(2, 3, 5, 3)
"""
struct chunk_to_send
    pol_number::Int
    first_idx::Int
    last_idx::Int
    dest_rank::Int
end

function generate_noisechunks_to_send(chunks, baseline_length_s, fsamp_hz, total_number_of_polarimeters)

    commsize = length(chunks)
    noisechunks_to_send = []

    for cur_rank in 0:(commsize - 1) #loop on ranks
        pol_cur_rank  = cur_rank + 1 #polarimeter assigned to cur_rank for noise production
        chunk_cur_rank = []

        #Only the first "total_number_of_polarimeters" ranks will generate noise, so they will have something to send.
        (0 <= cur_rank < total_number_of_polarimeters) || continue

        for i in 1:length(chunks) #loop on chunks
            rank_it = i - 1
            for j in 1:length(chunks[i]) #loop on detectors per rank
                if chunks[i][j].pol_number == pol_cur_rank #if cur_rank is going to produce noise for the "correct" polarimeter
                    first_idx = round(Int64, (chunks[i][j].first_idx - 1) * baseline_length_s * fsamp_hz + 1)
                    last_idx = round(Int64, first_idx + chunks[i][j].num_of_elements * baseline_length_s * fsamp_hz - 1)
                    dest_rank = rank_it
                    chunk = chunk_to_send(pol_cur_rank, first_idx, last_idx, dest_rank)
                    chunk_cur_rank = append!(chunk_cur_rank, [chunk])
                end
            end
        end
        noisechunks_to_send = append!(noisechunks_to_send, [chunk_cur_rank])
    end

    noisechunks_to_send
end

"""
    generate_noise_mpi(chunks, baselines_per_process, baseline_length_s, total_time_s, fsamp_hz, sigma_k, fknee_hz, slope, rank; comm = nothing, input_seed = missing)

This MPI based function can be used to generate white noise and 1/f noise when the simulation is splitted in various MPI processes.

Input parameters:

- `chunks`: data chunks, usually obtained using [`split_tod_mpi`](@ref)
- `baselines_per_process`: array containing the number of 1/f baselines to simulate for each process. Use [`split_into_n`](@ref) (see below)
- `baseline_length_s`: the length (in seconds) of each 1/f noise baseline
- `total_time_s`: total duration (in seconds) of data acquisition
- `fsamp_hz`: the sampling frequency (in Hz)
- `sigma_k`: an array containing the RMS of the temperature (in K) for
  each polarimeter. Usually `sigma_k[i] = tsys_k[i] / sqrt(β_hz[i] * τ_s`
- `fknee_hz`: an array containing the knee frequency (in Hz) for each polarimeter
- `slope`: an array containing the oof slope for each polarimeter

Keywords:

- `rank`: the rank of current process
- `comm`: the MPI communicator (defaults to `nothing`)
- `input_seed`: seed for the pseudo-random generator (defaults to
  `missing`; in this case a random number will be used).

Return value:

Noise TOD for the current rank.

```julia
# Number of baselines in each polarimeter
baselines_per_pol = Int64(total_time_s / baseline_length_s)
split_into_n(num_of_polarimeters * baselines_per_pol, num_of_MPI_proc)
```

The function behaves differently according to the number of processes used:

1. If the number of processes is greater than the number of polarimeters:
    the first "number of polarimeters" processes will simultaneously produce the noise, one process for one polarimeter
    (rank 0 for pol.1, rank 1 for pol.2 and so on).
    After that, the noise is redistributed between the ranks according to the input chunk structure.

2. If the number of processes is lesser than the number of polarimeters,
    Each partial noise TOD is generated by rank 0 and then sent to the correspondent rank.

N.B. if you want to use this function without MPI, remember to put rank = 0 and comm = nothing

"""
function generate_noise_mpi(
    chunks,
    baselines_per_process,
    baseline_length_s,
    total_time_s,
    fsamp_hz,
    sigma_k,
    fknee_hz,
    slope;
    rank = 0,
    comm = nothing,
    input_seed = missing,
)
    maxrand = 10_000_000

    total_number_of_polarimeters = length(sigma_k)
    commsize = length(chunks)

    father_seed = ismissing(input_seed) ? rand(1:maxrand) : input_seed
    Random.seed!(father_seed)

    # Create random seeds for each polarimeters, and make sure
    # that there are no duplicated seeds (otherwise two polarimeters
    # would have the same noise realization!)
    seeds = Array{Int64}(undef, total_number_of_polarimeters)
    visited_seeds = Set(Int64[])
    for i in 1:total_number_of_polarimeters
        while true
            new_seed = rand(1:maxrand)
            if ! (new_seed in visited_seeds)
                seeds[i] = new_seed
                push!(visited_seeds, new_seed)
                break
            end
        end
    end

    ################ IF MORE PROCESSES THAN POLARIMETERS ################
    if(commsize >= total_number_of_polarimeters)

        noisechunks_to_send = generate_noisechunks_to_send(chunks, baseline_length_s, fsamp_hz, total_number_of_polarimeters)

        ##### NOISE PRODUCTION #####
        pol_number = rank + 1
        pol_noise = Float64[]
        noise = Array{Float64}(undef, round(Int, baselines_per_process[rank + 1] * baseline_length_s * fsamp_hz))

        if  0 <= rank < total_number_of_polarimeters  #if first "total_number_of_polarimeters" ranks
            rng = CorrNoise.OofRNG(Random.MersenneTwister(seeds[pol_number]), -abs(slope[pol_number]), 1.15e-5, fknee_hz[pol_number], fsamp_hz)
            samples_per_pol = total_time_s * fsamp_hz
            pol_noise = Float64[CorrNoise.randoof(rng) * sigma_k[pol_number] for i in 1:(samples_per_pol)]
        end

        comm != nothing && MPI.Barrier(comm)

        ##### SEND AND RECEIVE #####

        start_idx = 1
        end_idx = 0

        for i in 1:length(noisechunks_to_send)  #loop on ranks
            cur_rank = i - 1

            for j in 1:length(noisechunks_to_send[i]) #loop on the noise chunks produced by cur_rank
                start_idx_noise = noisechunks_to_send[i][j].first_idx
                end_idx_noise = noisechunks_to_send[i][j].last_idx
                dest_rank = noisechunks_to_send[i][j].dest_rank

                cur_num_noise_samples = end_idx_noise - start_idx_noise + 1

                #if rank has to send noise (since it produced it)
                if cur_rank == rank

                    #if the destination rank is itself, keep noise chunk
                    if dest_rank == rank
                        end_idx = start_idx + cur_num_noise_samples - 1
                        noise[start_idx:end_idx] = pol_noise[start_idx_noise:end_idx_noise]
                        start_idx += cur_num_noise_samples

                    #otherwise send to the destination rank
                    else
                        cur_noise_send = pol_noise[start_idx_noise:end_idx_noise]
                        MPI.Send(cur_noise_send, dest_rank, 0, comm)
                    end

                #if rank is a destination rank for someone
                elseif dest_rank == rank

                    cur_noise_rec = Array{Float64}(undef, cur_num_noise_samples)
                    end_idx = start_idx + cur_num_noise_samples - 1
                    MPI.Recv!(cur_noise_rec, cur_rank, 0, comm)
                    noise[start_idx:end_idx] = cur_noise_rec
                    start_idx += cur_num_noise_samples
                end

            end

        end

###################################################################

    ################IF LESS PROCESSES THAN POLARIMETERS################
    else
        if rank == 0
            previous_detector = chunks[1][1].pol_number
            noise = Float64[]  #noise rank 0

            rng = CorrNoise.OofRNG(Random.MersenneTwister(seeds[1]), -abs(slope[1]), 1.15e-5, fknee_hz[1], fsamp_hz)

            for i in 1:length(chunks)   #loop on ranks
                num_noise_samples = baselines_per_process[i] * baseline_length_s * fsamp_hz

                for j in 1:length(chunks[i]) #loop on detectors per rank
                    cur_detector = chunks[i][j].pol_number

                    if cur_detector != previous_detector  #if new detector generate new noise
                        rng = CorrNoise.OofRNG(Random.MersenneTwister(seeds[cur_detector]), -abs(slope[cur_detector]), 1.15e-5, fknee_hz[cur_detector], fsamp_hz)
                    end
                    cur_num_noise_samples = chunks[i][j].num_of_elements * baseline_length_s * fsamp_hz
                    cur_noise = Float64[CorrNoise.randoof(rng) * sigma_k[cur_detector] for i in 1:(cur_num_noise_samples)]
                    previous_detector =  cur_detector

                    if (i > 1)  #send to other ranks, not rank 0
                        MPI.Send(cur_noise, i - 1, 0, comm)
                    else
                        noise = append!(noise, cur_noise)
                    end
                end
            end
        else
            noise = Float64[]
            for j in 1:length(chunks[rank + 1])
                cur_noise = Array{Float64}(undef, chunks[rank + 1][j].num_of_elements * baseline_length_s * fsamp_hz)
                MPI.Recv!(cur_noise, 0, 0, comm)
                noise = append!(noise, cur_noise)
            end
        end

    end

    noise
end