text:
"""
    c2k(temp)

Convert temperature in Celsius to Kelvin.
"""
function c2k(temp)
    return temp + 273.15
end

"""
    k2c(temp)

Convert temperature in Kelvin to Celsius.
"""
function k2c(temp)
    return temp - 273.15
end

ext: jl


text:
#
# Newton's method
#
# Reference:
# Nicholas J. Higham, Computing the Polar Decomposition ---with Applications,
# SIAM J. Sci. Statist. Comput. Vol. 7, Num 4 (1986) pp. 1160-1174.
#

type NewtonAlg{T} <: PolarAlg
    maxiter::Int        # maximum number of iterations.
    scale::Bool         # whether to scale Newton iteration.
    verbose::Bool       # whether to show procedural information
    tol::T        # tolerance for convergence
    scale_tol::T  # tolerance for acceleration scaling

    function NewtonAlg(;maxiter::Integer=100,
                       verbose::Bool=false,
                       scale::Bool=true,
                       tol::Real=cbrt(eps(T)),
                       scale_tol::Real=eps(T)^(1/4))
        maxiter > 1 || error("maxiter must be greater than 1.")
        tol > 0 || error("tol must be positive.")
        scale_tol > 0 || error("scale_tol must be positive.")

        new(maxiter,
            scale,
            verbose,
            tol,
            scale_tol)
    end
end

function solve!{T}(alg::NewtonAlg{T},
                   X::Matrix{T}, U::Matrix{T}, H::Matrix{T})
    common_iter_scal!(NewtonUpdater(alg.scale, alg.scale_tol), X, U, H, alg.maxiter, alg.verbose, alg.tol)
end

type NewtonUpdater{T} <: PolarUpdater
    scale::Bool
    scale_tol::T
end

function update_U!{T}(upd::NewtonUpdater, U::Matrix{T})
    scale = upd.scale
    Uinv = Array(T, size(U))
    Uinvt = Array(T, size(U))
    copy!(Uinv, inv(U))

    # 1, Inf-norm scaling
    if scale
        g = (norm(Uinv,1) * norm(Uinv, Inf) / (norm(U,1) * norm(U, Inf)) )^(1/4)
    else
        g = one(T)
    end
    transpose!(Uinvt, Uinv)
    copy!(U, (g * U + Uinvt /g) / convert(T, 2))

end

# Newton Schulz algrithm

type NewtonSchulzAlg{T} <: PolarAlg
    maxiter::Int     # maximum number of iterations.
    verbose::Bool    # whether to show procedural information
    tol::T     # convergence tolerance.

    function NewtonSchulzAlg(;maxiter::Integer=100,
                             verbose::Bool=false,
                             tol::Real=cbrt(eps(T)))
        maxiter > 1 || error("maxiter must be greater than 1.")
        tol > 0 || error("tol must be positive.")

        new(int(maxiter),
            verbose,
            tol)
    end
end

function solve!{T}(alg::NewtonSchulzAlg,
                   X::Matrix{T}, U::Matrix{T}, H::Matrix{T})

    # Newton Schulz converge quadratically if norm(X) < sqrt(3)
    norm(X) < convert(T, sqrt(3)) || throw(ArgumentError("The norm of the input matrix must be smaller than sqrt(3)."))

    common_iter!(NewtonSchulzUpdater(), X, U, H, alg.maxiter, alg.verbose, alg.tol)
end

immutable NewtonSchulzUpdater <: PolarUpdater end

function update_U!{T}(upd::NewtonSchulzUpdater, U::Matrix{T})
    UtU = Array(T, size(U))
    At_mul_B!(UtU, U, U)
    copy!(U, 0.5*U*(3*I - UtU))
end


text:
function ImageSegmentation(
		blocks::Tuple{<:Image{N},<:Mask{N}},
        data=nothing;
		size=ntuple(i -> 128, N),
		aug_projections=DataAugmentation.Identity(),
		aug_image=DataAugmentation.Identity(),
        C=RGB{N0f8},
        computestats=false,
	) where N
	return SupervisedTask(
		blocks,
		(
			ProjectiveTransforms(size; augmentations=aug_projections),
			getimagepreprocessing(data, computestats; C=C, augmentations=aug_image),
        	OneHot()
		)
	)
end

"""
    ImageSegmentation(size, classes; kwargs...)

Learning task for image segmentation. Images are
resized to `size` and a class is predicted for every pixel.

## Keyword arguments

- `computestats = false`: Whether to compute image statistics on dataset `data` or use
    default ImageNet stats.
- `aug_projections = `[`DataAugmentation.Identity`](#): augmentation to apply during
  [`ProjectiveTransforms`](#) (resizing and cropping)
- `aug_image = `[`DataAugmentation.Identity`](#): pixel-level augmentation to apply during
  [`ImagePreprocessing`](#)
- `C = RGB{N0f8}`: Color type images are converted to before further processing. Use `Gray{N0f8}`
    for grayscale images.
"""
function ImageSegmentation(size::NTuple{N,Int}, classes::AbstractVector; kwargs...) where N
    blocks = (Image{N}(), Mask{N}(classes))
    return ImageSegmentation(blocks; size=size, kwargs...)
end

registerlearningtask!(FASTAI_METHOD_REGISTRY, ImageSegmentation, (Image, Mask))

# ## Tests

@testset "ImageSegmentation [task]" begin
    @testset "2D" begin
        task = ImageSegmentation((16, 16), 1:4)
        testencoding(getencodings(task), getblocks(task).sample)
        FastAI.checktask_core(task)
        @test_nowarn tasklossfn(task)
        @test_nowarn taskmodel(task, Models.xresnet18())
        @testset "Show backends" begin
            @testset "ShowText" begin
                FastAI.test_task_show(task, ShowText(Base.DevNull()))
            end
        end
    end
    @testset "3D" begin
        task = SupervisedTask(
            (Image{3}(), Mask{3}(1:4)),
            (
                ProjectiveTransforms((16, 16, 16), inferencefactor=8),
                ImagePreprocessing(),
                FastAI.OneHot()
            )
        )
        testencoding(getencodings(task), getblocks(task).sample)
        FastAI.checktask_core(task)
        @test_nowarn tasklossfn(task)
    end

end


text:
module MimiDICE2016R2

using Mimi
using XLSX: readxlsx

include("helpers.jl")
include("parameters.jl")

include("marginaldamage.jl")

include("components/totalfactorproductivity_component.jl")
include("components/grosseconomy_component.jl")
include("components/emissions_component.jl")
include("components/co2cycle_component.jl")
include("components/radiativeforcing_component.jl")
include("components/climatedynamics_component.jl")
include("components/damages_component.jl")
include("components/neteconomy_component.jl")
include("components/welfare_component.jl")

export constructdice, getdiceexcel

const model_years = 2015:5:2510

function constructdice(params)

    m = Model()
    set_dimension!(m, :time, model_years)

    #--------------------------------------------------------------------------
    # Add components in order
    #--------------------------------------------------------------------------

    add_comp!(m, totalfactorproductivity, :totalfactorproductivity)
    add_comp!(m, grosseconomy, :grosseconomy)
    add_comp!(m, emissions, :emissions)
    add_comp!(m, co2cycle, :co2cycle)
    add_comp!(m, radiativeforcing, :radiativeforcing)
    add_comp!(m, climatedynamics, :climatedynamics)
    add_comp!(m, damages, :damages)
    add_comp!(m, neteconomy, :neteconomy)
    add_comp!(m, welfare, :welfare)

    #--------------------------------------------------------------------------
    # Make internal parameter connections
    #--------------------------------------------------------------------------

    # Socioeconomics
    connect_param!(m, :grosseconomy, :AL, :totalfactorproductivity, :AL)
    connect_param!(m, :grosseconomy, :I, :neteconomy, :I)
    connect_param!(m, :emissions, :YGROSS, :grosseconomy, :YGROSS)

    # Climate
    connect_param!(m, :co2cycle, :E, :emissions, :E)
    connect_param!(m, :radiativeforcing, :MAT, :co2cycle, :MAT)
    connect_param!(m, :climatedynamics, :FORC, :radiativeforcing, :FORC)

    # Damages
    connect_param!(m, :damages, :TATM, :climatedynamics, :TATM)
    connect_param!(m, :damages, :YGROSS, :grosseconomy, :YGROSS)
    connect_param!(m, :neteconomy, :YGROSS, :grosseconomy, :YGROSS)
    connect_param!(m, :neteconomy, :DAMAGES, :damages, :DAMAGES)
	connect_param!(m, :neteconomy, :SIGMA, :emissions, :SIGMA)
    connect_param!(m, :welfare, :CPC, :neteconomy, :CPC)

    #--------------------------------------------------------------------------
    # Set external parameter values
    #--------------------------------------------------------------------------
    # Set unshared parameters - name is a Tuple{Symbol, Symbol} of (component_name, param_name)
    for (name, value) in params[:unshared]
        update_param!(m, name[1], name[2], value)
    end

    # Set shared parameters - name is a Symbol representing the param_name, here
    # we will create a shared model parameter with the same name as the component
    # parameter and then connect our component parameters to this shared model parameter
    add_shared_param!(m, :fco22x, params[:shared][:fco22x]) #Forcings of equilibrium CO2 doubling (Wm-2)
    connect_param!(m, :climatedynamics, :fco22x, :fco22x)
    connect_param!(m, :radiativeforcing, :fco22x, :fco22x)

    add_shared_param!(m, :l, params[:shared][:l], dims = [:time]) #Level of population and labor (millions)
    connect_param!(m, :grosseconomy, :l, :l)
    connect_param!(m, :neteconomy, :l, :l)
    connect_param!(m, :welfare, :l, :l)

    add_shared_param!(m, :MIU, params[:shared][:MIU], dims = [:time]) #Optimized emission control rate results from DICE2016R (base case)
    connect_param!(m, :neteconomy, :MIU, :MIU)
    connect_param!(m, :emissions, :MIU, :MIU)

    return m

end

function getdiceexcel(;datafile = joinpath(dirname(@__FILE__), "..", "data", "DICE2016R-090916ap-v2-REVISEDtoR2.xlsm"))
    params = getdice2016r2excelparameters(datafile)

    m = constructdice(params)

    return m
end

# get_model function for standard Mimi API: use the Excel version
get_model = getdiceexcel

include("mcs.jl")

end # module


text:
function leftorthoQR_MPO!(mpoin)
    n = length(mpoin);
    mpoin[1] = permutedims(mpoin[1],[1,2,4,3])
    for i = 1:n-1
        (dl,d,_,dr) = size(mpoin[i])
        bi = reshape(mpoin[i],(dl*d^2,dr))
        q,r = qr(bi);
        mpoin[i] = permutedims(reshape(q,dl,d,d,size(q,2)),[1,2,4,3])
        mo = mpoin[i+1]
        @tensor mo[-1,-2,-4,-3] := r[-1,1]*mo[1,-2,-3,-4]
        mpoin[i+1] = mo
    end
#    (dl,d,_,dr) = size(mpoin[n])
#    mo = reshape(mpoin[n],(dl,d^2,dr))
    mpoin[n] = permutedims(mpoin[n],[1,2,4,3])
#    @tensor no = real(scalar(mo[1,2,3]*conj(mo[1,2,3])))
#    return mpoin, no;
    return mpoin
end

"""
leftorthonormalizationQR!(A)

Left-orthonormalizes the MPS A in place and returns the norm
"""
function leftorthonormalizationQR!(mpsin)

    no = leftorthogonalizeQR!(mpsin)
    mpsin[end] = mpsin[end]./sqrt(no)

    return no

end

"""
leftorthogonalizeQR!(A)

Left-orthogonalizes the MPS A in place and calculate the norm, but doesn't actually normalize
"""
function leftorthogonalizeQR!(mpsin)

    n = length(mpsin)

    for i = 1:(n - 1)

        mpsin[i], r = leftorth(mpsin[i])
        mpsin[i + 1] = prod_LR(r, mpsin[i + 1])

    end

    no = sum(abs2.(mpsin[n][:]))

end

"""
rightorthonormalizationQR(A)

Right-orthonormalizes the MPS A in place and returns the norm
"""
function rightorthonormalizationQR!(mpsin)

    no = rightorthogonalizeQR!(mpsin)
    mpsin[1] = mpsin[1]./sqrt(no)

    return no

end
"""
rightorthogonalizeQR(A)

Right-orthogonalizes the MPS A in place and returns the norm
"""
function rightorthogonalizeQR!(mpsin)

    n = length(mpsin)

    for i = n:-1:2

        mpsin[i], r = rightorth(mpsin[i])
        mpsin[i - 1] = prod_LR(mpsin[i - 1], r)

    end

    no = sum(abs2.(mpsin[1][:]))

end

"""
function to normalize left orthogonal MPS
"""

function normalize_lo!(mps)

    no = sum(abs2.(mps[end][:]))
	mps[end] = mps[end]./sqrt(no)
	return no

end

function normalize_lo(mps)

    no = sum(abs2.(mps[end][:]))
    mps[end] = mps[end]./sqrt(no)
    return mps

end

"""
function to get norm of left orthogonal MPS
"""

function norm_lo(mps)

    no = sum(abs2.(mps[end][:]))

end

"""
function to normalize left orthogonal MPS
"""

function normalize_ro!(mps)

    no = sum(abs2(mps[1][:]))
	mps[1] = mps[1]./sqrt(no)
	return no

end

"""
function to get norm of right orthogonal MPS
"""

function norm_ro(mps)

    no = sum(abs2(mps[1][:]))

end

"""
Separate a MPS tensor A into a left orthonormal tensor AL and a bond
factor C
"""

function leftorth(A)

    Dl, d, Dr = size(A)
    Q, C = qr(reshape(A, d*Dl, Dr))
    Q = reshape(Q, Dl, d, size(Q, 2))
    return Q, C

end

"""
Separate a MPS tensor A into a left orthonormal tensor AL and a bond
factor C in place
"""

function leftorth!(A)

    Dl, d, Dr = size(A)
    A = reshape(A, d*Dl, Dr)
    Q = qrq!(A)
    reshape(Q, Dl, d, size(Q, 2))

end

function qrq_julia!(A::Matrix)

    QR = qrfact!(A)
    Array(QR[:Q])

end

function qrq!(A::Matrix)

    LAPACK.orgqr!(LAPACK.geqrf!(A)...)

end

function qrq!(Agpu::AbstractMatrix)

    A = collect(Agpu)
    n = min(size(A)...)
    LAPACK.orgqr!(LAPACK.geqrf!(A)...)
    copy!(Agpu, A)
    if n < size(A, 2)
        Agpu[:, 1:n]
    else
        Agpu
    end

end

"""
Separate a MPS tensor A into a right orthonormal tensor AR and a bond
factor C
"""
function rightorth(A)

    Dl, d, Dr = size(A)
    Q, C = qr(reshape(A, Dl, d*Dr).')
    Q = reshape(Q.', size(Q, 2), d, Dr)
    return Q, C.'

end

function rightorth_lq(A)

    Dl, d, Dr = size(A)
    C, Q = lq(reshape(A, Dl, d*Dr))
    Q = reshape(Q, size(Q, 1), d, Dr)
    return Q, C

end

"""
Separate a MPS tensor A into a right orthonormal tensor AR and a bond
factor C in place
"""

function rightorth!(A)

    Dl, d, Dr = size(A)
    A = reshape(A, Dl, d*Dr).'
    Q = qrq!(A)
    reshape(Q.', size(Q, 2), d, Dr)

end

function rightorth_lq!(A::Array{<:Number, 3})

    Dl, d, Dr = size(A)
    A = reshape(A, Dl, d*Dr)
    Q = lqq!(A)
    reshape(Q, size(Q, 1), d, Dr)

end

function lqq_julia!(A::Matrix)

    LQ = lqfact!(A)
    Array(LQ[:Q])

end

function lqq!(A::Matrix)

    LAPACK.orglq!(LAPACK.gelqf!(A)...)

end

function lqq!(Agpu::AbstractMatrix)

    A = collect(Agpu)
    m = min(size(A)...)
    LAPACK.orglq!(LAPACK.gelqf!(A)...)
    if m < size(A, 1)
        Agpu = Agpu[1:m, :]
        copy!(Agpu, A)
    else
        copy!(Agpu, A)
    end

end


text:
function chineseremainder(n::Array, a::Array)
    Π = prod(n)
    mod(sum(ai * invmod(Π ÷ ni, ni) * Π ÷ ni for (ni, ai) in zip(n, a)), Π)
end

@show chineseremainder([3, 5, 7], [2, 3, 2])

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using CSV

f = CSV.File(joinpath("D:\\", "xwd-customers.csv"))


text:
module LogisticOptTools

using LinearAlgebra, SparseArrays, Random
using Statistics
using Optim

include("utils.jl")
include("loggers.jl")
include("io/libsvm_parser.jl")

include("dataset.jl")
include("loss.jl")
include("penalty.jl")

include("logistic/dense_logistic.jl")
include("logistic/sparse_logistic.jl")
include("logistic/dual_logistic.jl")

include("model.jl")
include("optim.jl")

end # module


text:
using FileIO, NRRD, ImageCore

Alast = nothing
for (A, name) in ((rand(UInt8, 3, 5), "uint8"),
                  (reinterpret(N0f8, rand(UInt8, 3, 5)), "u8"),
                  (rand(Gray{N0f8}, 3, 5), "gray_u8"),
                  (rand(RGB{N0f8}, 3, 5), "rgb_u8"),
                  (rand(Gray{N0f16}, 3, 5), "gray_u16"),
                  (reinterpret(Gray{N4f12}, rand(0x0000:0x0fff, 3, 5)), "gray_u12"))
    fn = joinpath(dirname(@__FILE__), "io", "eltype_$name.nrrd")
    save(fn, A)
    Alast = A
end

# Now manually tweak for the "hex" version
fn = joinpath(dirname(@__FILE__), "io", "eltype_gray_u12_hex.nrrd")
save(fn, Alast, props=Dict("sample units" => "gray 0x0fff"))


text:
module LAR
# basic implementation of lar-core.jl module

using JSON
using PyCall
@pyimport larlib as p

export p,larview,viewexploded,viewLarIndices,lar2hpc,cellComplex,
json2larmodel,rebase,rotate,scalate,scale
export lines2lar,larFromLines,larModelCheck,boxBucketing,lar2boxes,
randomLines,vcode,lineIntersection,lineBucketIntersect,segmentIntersect,
biconnectedComponents
export boundary,boundaryOp
export Surface,Volume,chainAreas

using IntervalTrees

# root of LAR type hierarchy
abstract LarModel

# top LAR subtypes
abstract LarGeometry <: LarModel
abstract CellComplex <: LarModel
abstract ChainComplex <: LarModel

# concrete type for models
type Model <: LarModel
	Verts
	Lar
end

# `Cells` type is a list (of lists of integers)
typealias Cells Array{Any,1}

# concrete type `LAR` is a quadruple of `Cells`
type Lar <: CellComplex
    VV::Cells
    EV::Cells
    FV::Cells
    CV::Cells
end

# constructor of an empty quadruple of lists
function Lar()
	Lar(Any[],Any[],Any[],Any[])
end

# constructor of a `Lar` instance (some fields possibly empty)
function Lar(bases::Array{Any,1})
	Lar( vv::Array{Any,2}, ev::Array{Any,2}, fv::Array{Any,2}, cv::Array{Any,2} )
end

# `Basis` type as alias for Sparse binay Matrix by Columns
typealias Basis SparseMatrixCSC

# `ChainBases` is a quadruple (just 3D, for now) of `Basis` binary matrices
type ChainBases <: ChainComplex
    M0::Basis
    M1::Basis
    M2::Basis
    M3::Basis
end

# constructor of an empty `ChainBases`
function ChainBases()
	m0 = cellComplex()
	m1 = cellComplex()
	m2 = cellComplex()
	m3 = cellComplex()
	ChainBases(m0,m1,m2,m3)
end

# constructor of a Basis instance from its cells
function cellComplex(cells::Cells)
	I,J,V = Int[],Int[],Int[]
	for (i,row) in enumerate(cells)
		for (j,col) in enumerate(row)
			push!(I,i); push!(J,col); push!(V,1)
		end
	end
	sparse(I,J,V)
end

function cellComplex(cells::Array{Array{Int64,1},1})
	I,J,V = Int[],Int[],Int[]
	m,n = length(cells),maximum(vcat(cells...))
	for j=1:m
		for i=1:length(cells[j])
			push!(I,cells[j][i]); push!(J,j); push!(V,1)
		end
	end
	sparse(I,J,V)
end

# constructor of a Basis instance from its cells
function cellComplex(cells::Array{Any,2})
	I,J,V = Int[],Int[],Int[]
	m,n = size(cells)
	for i=1:m
		for j=1:n
			push!(I,cells[i,j]); push!(J,j); push!(V,1)
		end
	end
	sparse(I,J,V)
end

# constructor of a Basis instance from its cells
function cellComplex(cells::Array{Int64,2})
	I,J,V = Int[],Int[],Int[]
	m,n = size(cells)
	for i=1:m
		for j=1:n
			push!(I,cells[i,j]); push!(J,j); push!(V,1)
		end
	end
	sparse(I,J,V)
end

# constructor of an empty Basis instance
function cellComplex()
	sparse(Int[],Int[],Int[])
end

# alias type used for boundary and coboundary operators
typealias Map SparseMatrixCSC

# concrete type `CCMap` (for Chain-Cochain Map) for chain complexes (3D)
type CCMap <: ChainComplex
    B1::Map
    B2::Map
    B3::Map
end

# constructor of empty `CCMap`
function CCMap()
	CCMap(sparse(Int[],Int[],Int[]))
end

# alias `Verts` for array of vertices in any embedding space (nD)
typealias Verts Array{Float64,2}

# Constructor of vertices from a Lar dictionary
function Verts(lar::Dict)
	verts = get(lar, "V", 0)
	if verts != 0
		if length(size(verts))==1
			verts = map(Float64, reduce(hcat,verts))
			vs = Verts(verts)
		else
			println("\nerror: vertex input")
		end
	end
end

# conversion from a JSON dictionary to a concrete `Model` type,
# i.e. to a pair of (Verts, Lar) types
function json2larmodel(strng::AbstractString)
	lardict = JSON.parse(strng)
	vs = Verts(lardict)
	verts = get(lardict, "VV", 0)
	edges = get(lardict, "EV", 0)
	faces = get(lardict, "FV", 0)
	cells = get(lardict, "CV", 0)
	lar = Lar()
	if verts != 0 lar.VV = verts end
	if edges != 0 lar.EV = edges end
	if faces != 0 lar.FV = faces end
	if cells != 0 lar.CV = cells end
	Model(vs,lar)
end

# used to convert a 2-array to a list of lists
function listOfList(cellArray::Array{Any,2})
	out = []
	for k=1:size(EV1)[1]
		cell = copy(cellArray[k,:])
		cell = cell[cell .!= ""]
		push!(out, cell)
	end
	out
end

# transform a 0-based array to a 1-based array
function rebase(faces::Array{Any,1},shift=+1)
	Any[Int[v+shift for v in face] for face in faces]
end

# transform a 0-based array to a 1-based array
function rebase(faces::Array{Any,2},shift=+1)
	faces + shift
end

# transform a `larlib` pair (varts,bases) into a Julia triple
# of instances of (Verts,Lar,ChainComplex) types
function importmodel(larmodel)
	if typeof(larmodel) == Tuple{Array{Any,2},Array{Any,1}}
		verts,bases = larmodel
		n = length(verts)
		chainbases = ChainBases()
		cells = Lar()
		v =  map(Float64,verts)'
		if length(bases) in (3,4)
			vv = rebase(bases[1]')
			ev = rebase(bases[2]')
			fv = rebase(bases[3]')
			cells.VV = Vector{Int}[vv[:,i] for i=1:size(vv,2)]
			cells.EV = Vector{Int}[ev[:,i] for i=1:size(ev,2)]
			cells.FV = Vector{Int}[fv[:,i] for i=1:size(fv,2)]
			chainbases.M0 = cellComplex(vv)
			chainbases.M1 = cellComplex(ev)
			chainbases.M2 = cellComplex(fv)
		end
		if length(bases) == 4
			cv = rebase(larmodel[2][4]')
			cells.CV = Vector{Int}[cv[:,i] for i=1:size(cv,2)]
			chainbases.M3 = cellComplex(cv)
		end
		v,cells,chainbases
	end
end

function lar2hpc(V::Array{Float64,2}, EV::Array{Any,1})
	EV = map(Int32, EV)
	lar2hpc(V,EV)
end

function lar2hpc(V::Array{Float64,2}, EV::Array{Array{Int64,1},1})
	EV = hcat(EV...)
	lar2hpc(V,EV)
end

# function lar2hpc(V::Array{Float64,2}, EV::Array{Array{Int64,1},1})
# 	EV = map(Any, EV)
# 	lar2hpc(V,EV)
# end

function lar2hpc(V::Array{Float64,2}, EV::Array{Int32,2})
	a,b = PyObject(V'), PyObject(EV')
	verts = PyObject(a[:tolist]())
	cells = PyObject(b[:tolist]())
	p.MKPOL([verts,cells,1])
end

function lar2hpc(V::Array{Float64,2}, EV::Array{Int64,2})
	EV = map(Int32, EV)
	a,b = PyObject(V'), PyObject(EV')
	verts = PyObject(a[:tolist]())
	cells = PyObject(b[:tolist]())
	p.MKPOL([verts,cells,1])
end

# visualize an HPC value from a Julia pair (Verts,Cells)
function larview(V::Array{Float64,2}, EV::Array{Int64,2})
	p.VIEW(lar2hpc(V,EV))
end

# visualize an HPC value from a Julia pair (Verts,Cells)
function larview(V::Array{Float64,2}, EV::Array{Any,2})
	EV = map(Int64, EV)
	larview(V,EV)
end

# visualize an HPC value from a Julia pair (Verts,Cells)
function larview(V::Array{Any,2}, EV::Array{Any,1})
	a,b = PyObject(V'), PyObject(EV)
	p.VIEW(p.MKPOL([a,b,1]))
end

# visualize an HPC value from a Julia pair (Verts,Cells)
function larview(V::Array{Float64,2}, EV::Array{Any,1})
	p.VIEW(lar2hpc(V,EV))
end

# visualize an HPC value from a Julia pair (Verts,Cells)
function larview(V::Array{Float64,2}, EV::Array{Array{Int64,1},1})
	a,b = PyObject(V'), PyObject(EV)
	verts = PyObject(a[:tolist]())
	cells = b
	p.VIEW(p.MKPOL([verts,cells,1]))
end

# visualize an HPC value from a Julia pair (Verts,Cells)
function larview(V::Array{Any,2}, EV::Array{Any,2})
	a,b = PyObject(V'), PyObject(EV')
	verts = a
	cells = b
	p.VIEW(p.MKPOL([verts,cells,1]))
end

# generate the triple of scaling factors needed by `larlib` explode
function scalingargs(scaleargs)
	if length(scaleargs)==1
		sx = sy = sz = scaleargs[1]
	elseif length(scaleargs)==3
		sx, sy, sz = scaleargs
	else
		println("error: wrong number of scaling arguments")
	end
	sx, sy, sz
end

# visualise an exploded `larlib` pair from a Julia pair (Verts,Cells)
function viewexploded(v::Array{Float64,2}, fv::Array{Array{Int64,1},1}, scaleargs=1.2)
	sx, sy, sz = scalingargs(scaleargs)
	a = PyObject(v')
	b = PyObject(Array{Any,1}[fv[k]-1 for k in 1:length(fv)])
	verts = PyObject(a[:tolist]())
	p.VIEW(p.EXPLODE(sx,sy,sz)(p.MKPOLS((verts,b))))
end

# visualise an exploded `larlib` pair from a Julia pair (Verts,Cells)
function viewexploded(v::Array{Float64,2}, fv::Array{Any,1}, scaleargs=1.2)
	sx, sy, sz = scalingargs(scaleargs)
	a = PyObject(v')
	b = PyObject(Array{Any,1}[fv[k]-1 for k in 1:length(fv)])
	verts = PyObject(a[:tolist]())
	p.VIEW(p.EXPLODE(sx,sy,sz)(p.MKPOLS((verts,b))))
end

# visualise an exploded `larlib` pair from a Julia pair (Verts,Cells)
function viewexploded(v::Array{Any,2}, fv::Array{Any,1}, scaleargs=1.2)
	sx, sy, sz = scalingargs(scaleargs)
	verts = PyObject(v')
	if typeof(verts)==PyObject
		a = verts
	else
		a = PyObject(verts[:tolist]())
	end
	b = PyObject(Array{Any,1}[fv[k]-1 for k in 1:length(fv)])
	p.VIEW(p.EXPLODE(sx,sy,sz)(p.MKPOLS((a,b))))
end

# visualise an exploded `larlib` pair from a Julia pair (Verts,Cells)
function viewexploded(v::Array{Int64,2}, fv::Array{Int64,2}, scaleargs=(1.2,))
	v = map(Float64,v)
	viewexploded(v, fv, scaleargs)
end

# visualise an exploded `larlib` pair from a Julia pair (Verts,Cells)
function viewexploded(v::Array{Float64,2}, fv::Array{Int64,2}, scaleargs=(1.2,))
	sx, sy, sz = scalingargs(scaleargs)
	a,b = PyObject(v'), PyObject(fv'-1)
	verts = PyObject(a[:tolist]())
	cells = PyObject(b[:tolist]())
	p.VIEW(p.EXPLODE(sx,sy,sz)(p.MKPOLS((verts,cells))))
end

# visualise an exploded `larlib` pair from a Julia pair (Verts,Cells)
function viewexploded(v::Array{Float64,2}, fv::Array{Any,2}, scaleargs=(1.2,))
	fv = map(Int64, fv)
	viewexploded(v, fv, scaleargs)
end

# visualise an exploded `larlib` pair from a Julia pair (Verts,Cells)
function viewexploded(v::Array{Any,2}, fv::Array{Any,2}, scaleargs=(1.2,))
	sx, sy, sz = scalingargs(scaleargs)
	fv = map(Int,fv)-1
	cells = [(fv[:,k]) for k in 1:size(fv,2)]
	a,b = PyObject(v'), PyObject(cells)
	p.VIEW(p.EXPLODE(sx,sy,sz)(p.MKPOLS((a,b))))
end

function viewexploded(v::Array{Any,2}, fv::Array{Int64,2}, scaleargs=(1.2,))
	sx, sy, sz = scalingargs(scaleargs)
	fv = fv-1
	cells = [(fv[:,k]) for k in 1:size(fv,2)]
	a,b = PyObject(v'), PyObject(cells)
	p.VIEW(p.EXPLODE(sx,sy,sz)(p.MKPOLS((a,b))))
end

# embed vertices array in a higher dimensional space, adding `ndims` zero coords
function embed(verts::Verts; ndims=1)
	vs = copy(verts)
	thetype = typeof(vs[1,1])
	added = zeros(thetype, (ndims, size(vs)[2]))
	vs = [vs ; added]
	out = Verts(vs)
end

# translate the columns of `V` matrix by sum with `t` vector
function translate( t, V )
	broadcast(+,t,V)
end

# scale the columns of `V` matrix by product times `s` vector
function scale( s, V )
	broadcast(*,s,V)
end

# rotate the columns of `V` matrix by properly using the `args` parameters
function rotate(args,V)
	n = length(args)
	if n == 1 # rotation in 2D
		angle = args[1]; Cos = cos(angle); Sin = sin(angle)
		mat = eye(2)
		mat[1,1] = Cos;    mat[1,2] = -Sin;
		mat[2,1] = Sin;    mat[2,2] = Cos;
	elseif n == 3 # rotation in 3D
		mat = eye(3)
		args = PyObject(PyObject(args)[:tolist]())
		angle = p.VECTNORM(args); axis = p.UNITVECT(args)
		Cos = cos(angle); Sin = sin(angle)
		# elementary rotations (in 3D)
		if axis[2]==axis[3]==0.0    # rotation about x
			mat[2,2] = Cos;    mat[2,3] = -Sin;
			mat[3,2] = Sin;    mat[3,3] = Cos;
		elseif axis[1]==axis[3]==0.0    # rotation about y
			mat[1,1] = Cos;    mat[1,3] = Sin;
			mat[3,1] = -Sin;    mat[3,3] = Cos;
		elseif axis[1]==axis[2]==0.0    # rotation about z
			mat[1,1] = Cos;    mat[1,2] = -Sin;
			mat[2,1] = Sin;    mat[2,2] = Cos;
		# general rotations (in 3D)
		else  # general 3D rotation (Rodrigues' rotation formula)
			I = eye(3) ; u = axis;
			Ux = Array([
				 0        -u[3]      u[2];
				 u[3]        0      -u[1];
				-u[2]      u[1]        0 ])
			UU = Array([
				 u[1]*u[1]    u[1]*u[2]    u[1]*u[3];
				 u[2]*u[1]    u[2]*u[2]    u[2]*u[3];
				 u[3]*u[1]    u[3]*u[2]    u[3]*u[3] ])
			mat = Cos*I + Sin*Ux + (1.0-Cos)*UU
		end
	end
	mat*V
end

# Identify -0.0 values with 0.0 values
function zeros!(W)
	for k=1:length(W)
	   if W[k] == -0.0 W[k] = 0.0 end
	end
	W
end

# vectorized parametric toroidal surface
function toroidalmap(V,r,R,m=40,n=80,angle1=2pi,angle2=2pi)
	V = map(Float64,V)
	V = (scale(V, [angle1/m angle2/n]))
	x =  (R .+ r.*cos(V[:,1])) .* cos(V[:,2])
	y =  (R .+ r.*cos(V[:,1])) .* sin(V[:,2])
	z = -r .* sin(V[:,1])
	[x y z]'
end

# valid (closed, i.e. no boundary) toroidal surface
function toroidalsurface(r,R,m=10,n=20,angle1=2pi,angle2=2pi)
	larmodel = p.larCuboids((m,n),true)
	verts,cells,chainbases = importmodel(larmodel)
	V = toroidalmap(verts',r,R,m,n)
	W,newcells,oldindex = larvalidate(V,cells)
	Z = map(Float64,hcat([V[:,k] for k in oldindex]...))
	Z,newcells
end

# topological validation of cellular complexes
function larvalidate(V,lar,prec=10^15) # prec for precision
	W = zeros!(map(round, V*prec)/prec)
	vdict = Dict([ ("$(W[:,k])",k) for k=1:size(W,2) ])
	verts = Dict(zip( keys(vdict), 1:length(vdict) ))
	newindex = convertindex(W,verts,prec)
	oldindex = invertindex(newindex)
	vs = [eval(parse(key))' for key in keys(vdict)]
	W = vcat(vs...)'
	newcells = Lar()
	if lar.FV != []
		FW = relink(lar.FV,newindex)
		newcells.FV = [ FW[:,k][:] for k=1:size(FW,2) ]
	end
	if lar.EV != []
		EW = relink(lar.EV,newindex)
		newcells.EV = [ EW[:,k][:] for k=1:size(EW,2) ]
	end
	if lar.VV != []
		VW = relink(lar.VV,newindex)
		newcells.VV = [ VW[:,k][:] for k=1:size(VW,2) ]
	end
	W,newcells,oldindex
end

# relinking of chain basis (d-cells, named FV) after topology validation
# including removal of multiple vertices (TODO ...)
# and
function relink(basis,newindex)
	m0,n0 = length(basis),length(basis[1])
	FV = reshape(vcat(basis...),n0,m0)
	FW = round(Int64, zeros(size(FV)))
	for k in 1:size(FV,2)
		for h in 1:size(FV,1)
			FW[h,k] = newindex[FV[h,k]]
		end
		FW[:,k] = sort(FW[:,k])
	end
	FW
end

# mapping between vertex indexing in dictionary and old indexing
function convertindex(W,vdict,precision)
	[vdict[ string(zeros!(map(round,W[:,k]*precision)/precision)) ] for k=1:size(W,2)]
end

function invertindex(index)
	n = max(index...)
	out = round(Int,zeros(n))
	for k=1:length(index)
		out[index[k]] = k
	end
	out
end

# Printing a Lar model (V,EV)
function printLar(modelname, model)
	V,EV = model

	function writefile(filename,data::Matrix)
		out_file = open(filename, "w")
		n,m = size(data)
		# file records printed by column
		for h=1:m
			for k=1:n-1
				print(out_file, data[k,h], ", ")
			end
			println(out_file, data[end,h])
		end
		close(out_file)
	end

	# vertex file
	filename = join([modelname,".V"])
	writefile(filename,V)
	# cell file
	filename = join([modelname,".EV"])
	writefile(filename,EV)
end

# Reading a Lar model (V,EV)
function readLar(modelname)
	filename = join([modelname,".V"])
	W = readcsv(filename)
	filename = join([modelname,".EV"])
	EW = readcsv(filename,Int64)
	W',EW'
end

function cols2any(EW)
	EZ = Any[]
	[push!(EZ, EW[:,k]) for k=1:size(EW,2)]
	EZ
end

# Produce an indexed view of Lar model (W,EW)
function viewLarIndices(W,EW,unit=1.0)
	EW = cols2any(EW)
	function convertData(W,EW)
		# Shifted vertices
		Z = zeros(length(W)+size(W,1))
		for k=1:length(W)
			Z[k + size(W,1)] = W[k]
		end
		Z = reshape(Z,size(W,1),size(W,2)+1)
		# shifted edges
		EZ = Any[[1,1]]
		for k=1:length(EW)
			push!(EZ, EW[k]+1)
		end
		Z,EZ
	end
	V,EV = convertData(W,EW)
	submodel = lar2hpc(V,hcat(EV...))
	VV = PyObject(Any[Any[k] for k in  0:size(V,2)-1])
	EZ = map(Array{Int32},EV-1)
	EV = PyObject(Any[PyObject(EZ[e])[:tolist]() for e=1:length(EZ) ])
	V = PyObject(Any[PyObject(V[:,v])[:tolist]() for v=1:size(V,2) ])
	scale = unit*maximum(p.SIZE(Any[1,2])(submodel))/6
	hpc = p.larModelNumbering(1,1,1)(V,PyObject([VV,EV]),submodel,scale)
	p.VIEW(hpc)
end

# Produce an indexed view of Lar model (W,EW,FW)
function viewLarIndices(W::Array{Float64,2},EW::Array{Int64,2},FW::Array{Int64,2},unit=1.0)
	EV = [EW[:,k] for k=1:size(FW,2)]
	return viewLarIndices(W,EV,FW,unit)
end

# Produce an indexed view of Lar model (W,EW,FW)
function viewLarIndices(W::Array{Float64,2},EW::Array{Int64,2},FW::Array{Array{Int64,1},1},
						unit=1.0)
	EW = cols2any(EW)
	function convertData(W,EW)
		# Shifted vertices
		Z = zeros(length(W)+size(W,1))
		for k=1:length(W)
			Z[k + size(W,1)] = W[k]
		end
		Z = reshape(Z,size(W,1),size(W,2)+1)
		# shifted edges
		EZ = Any[[1,1]]
		for k=1:length(EW)
			push!(EZ, EW[k]+1)
		end
		Z,EZ
	end
	V,EV = convertData(W,EW)
	submodel = lar2hpc(V,hcat(EV...))
	VV = PyObject(Any[Any[k] for k in  0:size(V,2)-1])
	EZ = map(Array{Int32},EV-1)
	EV = PyObject(Any[PyObject(EZ[e])[:tolist]() for e=1:length(EZ) ])
	FV = PyObject(Any[PyObject(FW[f])[:tolist]() for f=1:length(FW) ])
	V = PyObject(Any[PyObject(V[:,v])[:tolist]() for v=1:size(V,2) ])
	Scale = unit*maximum(p.SIZE(Any[1,2])(submodel))/6
	hpc = p.larModelNumbering(1,1,1)(V,PyObject([VV,EV,FV]),submodel,Scale)
	p.VIEW(hpc)
end

include("inters.jl")
include("boundary.jl")
include("integr.jl")

# end # module LAR
end


text:
module Data

using Random: shuffle!
using Base: @propagate_inbounds

include("dataloader.jl")
export DataLoader

## TODO: remove in v0.13 ##############
module MNIST
images() = error("Dataset is deprecated, use MLDatasets.jl instead.")
labels() = error("Dataset is deprecated, use MLDatasets.jl instead.")
end
module Iris
features() = error("Dataset is deprecated, use MLDatasets.jl instead.")
labels() = error("Dataset is deprecated, use MLDatasets.jl instead.")
end
module FashionMNIST
images() = error("Dataset is deprecated, use MLDatasets.jl instead.")
labels() = error("Dataset is deprecated, use MLDatasets.jl instead.")
end
module CMUDict
phones() = error("Dataset is deprecated, use MLDatasets.jl instead.")
symbols() = error("Dataset is deprecated, use MLDatasets.jl instead.")
rawdict() = error("Dataset is deprecated, use MLDatasets.jl instead.")
cmudict() = error("Dataset is deprecated, use MLDatasets.jl instead.")
end
module Sentiment
train() = error("Dataset is deprecated, use MLDatasets.jl instead.")
test() = error("Dataset is deprecated, use MLDatasets.jl instead.")
dev() = error("Dataset is deprecated, use MLDatasets.jl instead.")
end

export MNIST, Iris, FashionMNIST, CMUDict, Sentiment
#########################################

end#module


text:
n = rand(1:10)
Z = Companion(randn(n))

#Special properties
@test_approx_eq full(inv(Z)) inv(full(Z))

#Matvec product
b = randn(n)
@test_approx_eq Z*b full(Z)*b

m = rand(1:10)
A = randn(m, n)
@test_approx_eq A*Z A*full(Z)


text:
using Base: OneTo

struct tmpvars where {V0<:AbstractVector}
    u::V0
    v::V0
    gradis::Matrix{T}
end

function simloglik!(grad, info_matrix, tmpvars, model, theta, data, dograd)

    check_theta_ok(model, theta, data) # check that royalty model is OK
    compute_xbetas!(tmpvars, data) # does x'β for royalty & production data

    LL = zero(eltype(theta))
    for indiv in data
        gradi = _gradi(tmpvars, i)
        LL += simloglik_i!(gradi, tmpvars, model, theta, indiv, dograd)
    end

    if dograd
        gradis = _gradis(tmpvars)
        sum!(reshape(gradis, :, 1), )
        mul!(info_matrix, gradis', gradis)
    end

    return LL
end

function simloglik_i!(grad::AbstractVector, tmpvars, model, theta::AbstractVector, data, dograd::Bool)

    dograd && zero_gradm!(tmpvars)

    # royalty rate
    rc = RoyaltyTmpVar
    simloglik_royalty!(rc, royalty(model), theta, dograd)

    # production
    plc = ProductionLikelihoodComputations(data)
    simloglik_produce!(LL, produce(model), theta, plc, _ψ2(tmpvars))

    # drilling
    solve_value_function(_vftmp(tmpvars), drilling(model), theta, itype(data), dograd)
    simloglik_threaded!(grad, tmpvars, model, thet, data, dograd)

    # log likelihood & gradient
    if !dograd
        LL = logsumexp(llm) - log(M)
    else
        LL = logsumexp_and_softmax!(llm) - log(M)

        # royalty rate
        update_grad_royalty!(grad, rc, royalty(model), theta)

        # likelihood
        pgc = ProductionGradientComputations
        update_grad_produce!(grad, produce(model), theta, plc, pgc)

        # drilling
        mul!(drilling_gradient, gradm, llm)
        grad .+= gradtmps[1]
    end

    isfinite(LL) || @warn "infinite simulated log likelihood" #  for unit $i, which had $(j1length(data,i)) leases. SLLi = $lli.  θ[rng_d] = $(θ[rng_d])"
    return LL
end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))

@static if Sys.isapple()
    # download binaries from LLVM
    bin_prefix = "http://releases.llvm.org/6.0.0"
    url = "$bin_prefix/clang+llvm-6.0.0-x86_64-apple-darwin.tar.xz"
    tarball_hash = "0ef8e99e9c9b262a53ab8f2821e2391d041615dd3f3ff36fdf5370916b0f4268"
    # install
    mac_prefix = Prefix(joinpath(prefix, "clang+llvm-6.0.0-x86_64-apple-darwin"))
    products = [
        LibraryProduct(mac_prefix, String["libclang"], :libclang),
        ExecutableProduct(mac_prefix, "llvm-config", :llvm_config),
    ]
    install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    # write out a deps.jl file that will contain mappings for our products
    write_deps_file(joinpath(@__DIR__, "deps.jl"), products)
end

@static if Sys.islinux()
    if !haskey(ENV, "LLVM_CONFIG")
        @info "LLVM_CONFIG is not defined, looking for system llvm-config"
        try
            llvmcfg = readchomp(`which llvm-config`)
            isfile(llvmcfg) && (ENV["LLVM_CONFIG"] = llvmcfg;)
        catch err
            error("no system llvm-config: $err")
        end
    end
    LLVM_CONFIG = ENV["LLVM_CONFIG"]
    llvm_libdir = readchomp(`$(LLVM_CONFIG) --libdir`)
    if !ispath(llvm_libdir)
        error("llvm-config returned non-existent libdir path: '$(llvm_libdir)")
    end

    const libclang = joinpath(llvm_libdir, "libclang.so")

    open(joinpath(@__DIR__, "deps.jl"), "w") do f
        println(f, """
            import Libdl
            const libclang = "$libclang"
            const LLVM_CONFIG = "$LLVM_CONFIG"
            function check_deps()
                global libclang
                if !isfile(libclang)
                    error("\$(libclang) does not exist, Please re-run Pkg.build(\\"Clang\\"), and restart Julia.")
                end
                hdl = Libdl.dlopen_e(libclang)
                if hdl == C_NULL
                    error("\$(libclang) cannot be opened, Please re-run Pkg.build(\\"Clang\\"), and restart Julia.")
                end
            end
        """)
    end
    @info "build success!"
end

@static if Sys.iswindows()
    import Libdl

    const LLVM_WINDOWS_DEFAULT = joinpath(homedir(), "..", "..", "Program Files", "LLVM", "bin") |> normpath
    if !haskey(ENV, "LLVM_WINDOWS")
        @info "LLVM_WINDOWS is not defined, searching default directory"
        lib = Libdl.find_library(["libclang"], [LLVM_WINDOWS_DEFAULT])
        if lib == ""
            error("couldn't find libclang.dll in the default directory $LLVM_WINDOWS_DEFAULT.\n
                 Please download and install LLVM from http://releases.llvm.org/download.html and then run Pkg.build(\"Clang\").")
        else
            ENV["LLVM_WINDOWS"] = LLVM_WINDOWS_DEFAULT
        end
    end
    const libclang = Libdl.dlpath(joinpath(ENV["LLVM_WINDOWS"], "libclang.dll")) |> x->replace(x, "\\"=>"\\\\")
    open(joinpath(@__DIR__, "deps.jl"), "w") do f
        println(f, """
            import Libdl
            const libclang = "$libclang"
            function check_deps()
                global libclang
                if !isfile(libclang)
                    error("\$(libclang) does not exist, Please re-run Pkg.build(\\"Clang\\"), and restart Julia.")
                end
                hdl = Libdl.dlopen_e(libclang)
                if hdl == C_NULL
                    error("\$(libclang) cannot be opened, Please re-run Pkg.build(\\"Clang\\"), and restart Julia.")
                end
            end
        """)
    end
end


text:
module StatsBase
    import Base: length, isempty, eltype, values, sum, mean, show, quantile
    import Base: rand, rand!
    import Base.LinAlg: BlasReal

    export

    # common
    WeightVec, weights,

    # means
    geomean,     # geometric mean
    harmmean,    # harmonic mean
    trimmean,    # trimmed mean
    wmean,       # weighted mean

    # scalar_stats
    skewness,   # (standardized) skewness
    kurtosis,   # (excessive) kurtosis
    variation,  # ratio of standard deviation to mean
    sem,        # standard error of the mean, i.e. sqrt(var / n)
    mad,        # median absolute deviation
    middle,     # the mean of two real numbers
    midrange,   # the mean of minimum and maximum
    range,      # the difference between maximum and minimum
    percentile, # quantile using percentage (instead of fraction) as argument
    iqr,        # interquatile range
    nquantile,  # quantiles at [0:n]/n
    mode,       # find a mode from data
    modes,      # find all modes from data
    summarystats,   # summary statistics
    describe,       # print the summary statistics

    # counts
    addcounts!,     # add counts to an accumulating array or map
    counts,         # count integer values in given arrays
    proportions,    # proportions of integer values in given arrays
                    # (normalized version of counts)
    countmap,       # count distinct values and return a map
    proportionmap,  # proportions of distinct values returned as a map

    # ranking
    ordinalrank,    # ordinal ranking ("1234" ranking)
    competerank,    # competition ranking ("1 2 2 4" ranking)
    denserank,      # dense ranking ("1 2 2 3" ranking)
    tiedrank,       # tied ranking ("1 2.5 2.5 4" ranking)

    # corr
    autocov!, autocov,      # auto covariance
    autocor!, autocor,      # auto correlation
    crosscov!, crosscov,    # cross covariance
    crosscor!, crosscor,    # cross correlation
    pacf!, pacf,  # partial auto-correlation

    # rankcorr
    corspearman,       # spearman's rank correlation
    corkendall,        # kendall's rank correlation

    # sampling
    samplepair,     # draw a pair of distinct elements
    sample,         # sampling from a population
    wsample,        # sampling from a population with weights

    # empirical
    ecdf,           # empirical cumulative distribution function
    kde,            # kernel density estimation

    # misc
    rle,            # run-length encoding
    inverse_rle,    # inverse run-length encoding
    indexmap,       # construct a map from element to index
    levelsmap,      # construct a map from n unique elements to [1, ..., n]
    findat,         # find the position within a for elements in b
    indicatormat,   # construct indicator matrix

    # statistical models
    StatisticalModel,
    RegressionModel,

    coef,
    coeftable,
    confint,
    deviance,
    loglikelihood,
    nobs,
    stderr,
    vcov,
    predict,
    predict!,
    residuals,
    model_response

    # source files

    include("common.jl")
    include("misc.jl")
    include("means.jl")
    include("scalarstats.jl")
    include("counts.jl")
    include("ranking.jl")
    include("toeplitzsolvers.jl")
    include("corr.jl")
    include("rankcorr.jl")
    include("rand.jl")
    include("sampling.jl")
    include("empirical.jl")

    include("statmodels.jl")

    include("deprecates.jl")

end # module


text:
using ImageMagick, Documenter, ESDL, ESDLPlots, Cairo

newcubedir = mktempdir()
ESDLdir(newcubedir)
# Download Cube subset
if !isempty(ESDL.ESDLDefaults.cubedir[])
  c = S3Cube()
  csa = c[
    region = "South America",
    var = ["country_mask","c_emissions","gross", "net_ecosystem", "air_temperature_2m", "terrestrial_ecosystem", "land_surface_temperature"],
    time = 2003:2006
  ]
  saveCube(csa,"southamericacube", chunksize=(90,90,92,1))
  ESDL.ESDLDefaults.cubedir[] = joinpath(newcubedir,"southamericacube")
end

exampledir = joinpath(@__DIR__,"src","examples")
allex = map(readdir(exampledir)) do fname
  n = splitext(fname)[1]
  n => joinpath("examples",fname)
end

makedocs(
    modules = [ESDL, ESDLPlots],
    clean   = true,
    format   = Documenter.HTML(),
    sitename = "ESDL.jl",
    authors = "Fabian Gans",
    pages    = [ # Compat: `Any` for 0.4 compat
        "Home" => "index.md",
        "Examples" => allex,
        "Manual" => Any[
            "cube_access.md",
            "analysis.md",
            "plotting.md",
            "iotools.md"
        ],
        "Other functions" => "./lib/misc.md"
        ]
)

deploydocs(
    #deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/esa-esdl/ESDL.jl.git",
)


text:
abstract type AbstractPort end
abstract type AbstractSignalLine{Tv} end

abstract type AbstractBlock end
abstract type AbstractFunctionBlock <: AbstractBlock end
abstract type AbstractArithmeticBlock <: AbstractBlock end


text:
include("gumbel_softmax.jl")

using ReverseDiff

#= and(x, y) = Int.((x + y) .== 2) =#
#= or(x, y) = Int.((x + y) .>= 1) =#
#= xor(x, y) = Int.((x .* y) .== 0) =#
and(x, y) = x+y
or(x, y) = x-y
xor(x, y) = x*y

gt(X) = and(xor(X[1],X[5]), or(X[2],X[4])) #(1,3,2,4,8,5,7)

const N = 6
X = hcat([parse.(Int,vcat(bin(i,N)...)) for i=0:2^N-1]...)' #all combos of bits up to N

function f1(ps1, ps2, ps3, ps4, ps5, ps6, ps7, x)
    T = 1.0
    logits = log.(ps1)
    w = gumbel_softmax(logits, T, true)
    return w' * [and(f2(ps2,ps4,ps5,x),f3(ps3,ps6,ps7,x)), or(f2(ps2,ps4,ps5,x),f3(ps3,ps6,ps7,x)), xor(f2(ps2,ps4,ps5,x),f3(ps3,ps6,ps7,x)), x[1], x[2], x[3], x[4], x[5], x[6]]
end
function f2(ps2, ps4, ps5, x)
    T = 1.0
    logits = log.(ps2)
    w = gumbel_softmax(logits, T, true)
    return w' * [and(f4(ps4,x),f5(ps5,x)), or(f4(ps4,x),f5(ps5,x)), xor(f4(ps4,x),f5(ps5,x)), x[1], x[2], x[3], x[4], x[5], x[6]]
end
function f3(ps3, ps6, ps7, x)
    T = 1.0
    logits = log.(ps3)
    w = gumbel_softmax(logits, T, true)
    return w' * [and(f6(ps6,x),f7(ps7,x)), or(f6(ps6,x),f7(ps7,x)), xor(f6(ps6,x),f7(ps7,x)), x[1], x[2], x[3], x[4], x[5], x[6]]
end
function f4(ps4, x)
    T = 1.0
    logits = log.(ps4)
    w = gumbel_softmax(logits, T, true)
    return w' * [x[1], x[2], x[3], x[4], x[5], x[6]]
end
function f5(ps5, x)
    T = 1.0
    logits = log.(ps5)
    w = gumbel_softmax(logits, T, true)
    return w' * [x[1], x[2], x[3], x[4], x[5], x[6]]
end
function f6(ps6, x)
    T = 1.0
    logits = log.(ps6)
    w = gumbel_softmax(logits, T, true)
    return w' * [x[1], x[2], x[3], x[4], x[5], x[6]]
end
function f7(ps7, x)
    T = 1.0
    logits = log.(ps7)
    w = gumbel_softmax(logits, T, true)
    return w' * [x[1], x[2], x[3], x[4], x[5], x[6]]
end
function objective(ps1, ps2, ps3, ps4, ps5, ps6, ps7, x, y_target)
    abs2.(y_target - f1(ps1, ps2, ps3, ps4, ps5, ps6, ps7, x))
end

Y = [gt(X[i,:]) for i=1:size(X,1)]
ps1 = ones(9) ./ 9
ps2 = deepcopy(ps1)
ps3 = deepcopy(ps1)
ps4 = ones(6) ./ 6
ps5 = deepcopy(ps4)
ps6 = deepcopy(ps4)
ps7 = deepcopy(ps4)
#= j = 20 =#
#= f1(ps1, ps2, ps3, ps4, ps5, ps6, ps7, X[j,:]) =#
#= g = ReverseDiff.gradient(objective, (ps1, ps2, ps3, ps4, ps5, ps6, ps7, X[j,:], [Y[j]])) =#

function learn(α=0.1, N=5)
    Y = [gt(X[i,:]) for i=1:size(X,1)]
    ps1 = ones(9) ./ 9
    ps2 = deepcopy(ps1)
    ps3 = deepcopy(ps1)
    ps4 = ones(6) ./ 6
    ps5 = deepcopy(ps4)
    ps6 = deepcopy(ps4)
    ps7 = deepcopy(ps4)

    for i = 1:N
        println("i=$i")
        @show ps1
        @show ps2
        @show ps3
        @show ps4
        @show ps5
        @show ps6
        @show ps7
        for j = 1:size(X,1)
            g1,g2,g3,g4,g5,g6,g7,_,_ = ReverseDiff.gradient(objective, (ps1, ps2, ps3, ps4, ps5, ps6, ps7, X[j,:], [Y[j]]))
            #= @show ps1 =#
            #= @show ps2 =#
            #= @show ps3 =#
            #= @show ps4 =#
            #= @show ps5 =#
            #= @show ps6 =#
            #= @show ps7 =#
            #= @show g1 =#
            #= @show g2 =#
            #= @show g3 =#
            #= @show g4 =#
            #= @show g5 =#
            #= @show g6 =#
            #= @show g7 =#
            ps1 = update(ps1, g1, α)
            ps2 = update(ps2, g2, α)
            ps3 = update(ps3, g3, α)
            ps4 = update(ps4, g4, α)
            ps5 = update(ps5, g5, α)
            ps6 = update(ps6, g6, α)
            ps7 = update(ps7, g7, α)
        end
    end
    return indmax(ps1), indmax(ps2), indmax(ps3), indmax(ps4), indmax(ps5), indmax(ps6), indmax(ps7)
end
function update(ps, g, α)
    if !any(isnan.(g))
        ps -= α * g
    end
    ps = max.(min.(ps, 1.0), 1e-12)
    normalize!(ps,1)
    ps
end

#=  =#
#= function test(n_seeds=50) =#
#=     result = zeros(length(OP), length(OP)) #truth x result =#
#=     for (i, gt) in enumerate(OP) =#
#=         for s=1:n_seeds =#
#=             srand(s) =#
#=             imax = learn(gt) =#
#=             result[i, imax] += 1 =#
#=         end =#
#=     end =#
#=     result =#
#= end =#

#Issues:
#1. y_hard: tell ReverseDiff to hard forward, soft back
#2. variable temperature: tell ReverseDiff that T is a constant
#3. ReverseDiff throws and error when batch_size > 1


text:
function moi_sensorloc(optimizer, seed, n; verbose = false, test = false, scalar = false)

    rng = Random.MersenneTwister(seed)
    MOI.empty!(optimizer)
    if test
        @test MOI.is_empty(optimizer)
    end
    # Generate randomized problem data
    m, x_true, a, d, d_bar = sensorloc_data(seed, n)

    # Decision variable
    nvars = ProxSDP.sympackedlen(n + 2)
    X = MOI.add_variables(optimizer, nvars)
    Xsq = Matrix{MOI.VariableIndex}(undef, n + 2, n + 2)
    ProxSDP.ivech!(Xsq, X)
    Xsq = Matrix(LinearAlgebra.Symmetric(Xsq, :U))
    vov = MOI.VectorOfVariables(X)
    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(n + 2))

    # Constraint with distances from anchors to sensors
    for j in 1:n
        for k in 1:m
            if scalar
                MOI.add_constraint(optimizer,
                    MOI.ScalarAffineFunction([
                            MOI.ScalarAffineTerm(a[k][1]*a[k][1], Xsq[1,1]),
                            MOI.ScalarAffineTerm(a[k][2]*a[k][2], Xsq[2,2]),
                            MOI.ScalarAffineTerm(-2 * a[k][1], Xsq[1, j+2]),
                            MOI.ScalarAffineTerm(-2 * a[k][2], Xsq[2, j+2]),
                            MOI.ScalarAffineTerm(1.0, Xsq[j+2, j+2]),
                        ], 0.0), MOI.EqualTo(d_bar[k, j]^2))
            else
                MOI.add_constraint(optimizer,
                MOI.VectorAffineFunction(
                    MOI.VectorAffineTerm.([1], [
                        MOI.ScalarAffineTerm(a[k][1]*a[k][1], Xsq[1,1]),
                        MOI.ScalarAffineTerm(a[k][2]*a[k][2], Xsq[2,2]),
                        MOI.ScalarAffineTerm(-2 * a[k][1], Xsq[1, j+2]),
                        MOI.ScalarAffineTerm(-2 * a[k][2], Xsq[2, j+2]),
                        MOI.ScalarAffineTerm(1.0, Xsq[j+2, j+2]),
                    ]), -[d_bar[k, j]^2]), MOI.Zeros(1))
            end
        end
    end

    # Constraint with distances from sensors to sensors
    count, count_all = 0, 0
    for i in 1:n
        for j in 1:i - 1
            count_all += 1
            if Random.rand(rng) > 0.9
                count += 1
                if scalar
                    MOI.add_constraint(optimizer,
                        MOI.ScalarAffineFunction([
                            MOI.ScalarAffineTerm(1.0, Xsq[i+2,i+2] ),
                            MOI.ScalarAffineTerm(1.0, Xsq[j+2,j+2] ),
                            MOI.ScalarAffineTerm(-2.0, Xsq[i+2,j+2]),
                        ], 0.0), MOI.EqualTo(d[i, j]^2))
                else
                    MOI.add_constraint(optimizer,
                        MOI.VectorAffineFunction(
                            MOI.VectorAffineTerm.([1],[
                                MOI.ScalarAffineTerm(1.0, Xsq[i+2,i+2] ),
                                MOI.ScalarAffineTerm(1.0, Xsq[j+2,j+2] ),
                                MOI.ScalarAffineTerm(-2.0, Xsq[i+2,j+2]),
                            ]), -[d[i, j]^2]), MOI.Zeros(1))
                end
            end
        end
    end
    if verbose
        @show count_all, count
    end
    if scalar
        MOI.add_constraint(optimizer, MOI.SingleVariable(Xsq[1, 1]), MOI.EqualTo(1.0))
        MOI.add_constraint(optimizer, MOI.SingleVariable(Xsq[1, 2]), MOI.EqualTo(0.0))
        MOI.add_constraint(optimizer, MOI.SingleVariable(Xsq[2, 1]), MOI.EqualTo(0.0))
        MOI.add_constraint(optimizer, MOI.SingleVariable(Xsq[2, 2]), MOI.EqualTo(1.0))
    else
        MOI.add_constraint(optimizer, MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.([1],[MOI.ScalarAffineTerm(1.0, Xsq[1, 1])]), [-1.0]), MOI.Zeros(1))
        MOI.add_constraint(optimizer, MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.([1],[MOI.ScalarAffineTerm(1.0, Xsq[1, 2])]), [ 0.0]), MOI.Zeros(1))
        MOI.add_constraint(optimizer, MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.([1],[MOI.ScalarAffineTerm(1.0, Xsq[2, 1])]), [ 0.0]), MOI.Zeros(1))
        MOI.add_constraint(optimizer, MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.([1],[MOI.ScalarAffineTerm(1.0, Xsq[2, 2])]), [-1.0]), MOI.Zeros(1))
    end

    objf_t = [MOI.ScalarAffineTerm(0.0, Xsq[1, 1])]
    if false
        objf_t = [MOI.ScalarAffineTerm(1.0, Xsq[i,i]) for i in 1:n + 2]
    end
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(objf_t, 0.0))

    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(optimizer)

    objval = MOI.get(optimizer, MOI.ObjectiveValue())

    stime = -1.0
    try
        stime = MOI.get(optimizer, MOI.SolveTimeSec())
    catch
        println("could not query time")
    end

    Xsq_s = MOI.get.(optimizer, MOI.VariablePrimal(), Xsq)

    verbose && sensorloc_eval(n, m, x_true, Xsq_s)

    rank = -1
    status = 0
    if MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
        status = 1
    end
    return (objval, stime, rank, status)
end


text:
# Test file to check the projections
# workspace()
# include("../src/Solver.jl")
using Test, COSMOTestUtils, COSMO.Projections, LinearAlgebra, SparseArrays

rng = Random.MersenneTwister(1234)

nn = 1

@testset "Test projection onto zero cone" begin
        x = rand(rng, rand(rng, 1:100))
        Projections.zeroCone!(x, 1, length(x))
        @test x == zeros(size(x, 1))
end

@testset "Test projection onto nonNegativeOrthant" begin
    for iii = 1:nn
        dim = rand(rng, 1:100)
        x = rand(rng, -1:0.001:1, dim)
        Projections.nonNegativeOrthant!(x, 1, length(x))
        @test minimum(x) >= 0
    end
end

@testset "Test projection onto Box" begin
    for iii = 1:nn
        dim = rand(rng, 1:100)
        vMin = -1000
        vMax = 1000
        x = vMin + rand(rng, dim, 1) * (vMax - vMin)
        u = (vMin / 5) + rand(rng, 1e-6:1e-3:1e6, dim ,1) * ((vMax - vMin) / 5)
        l = -u
        xProj = Projections.box(x, l, u)

        projection_failed = false
        for kkk = 1:dim
            if xProj[kkk] < l[kkk] || xProj[kkk] > u[kkk]
                projection_failed = true
            end
        end
        @test projection_failed == false
    end
end

@testset "Test projection onto second-order cone" begin

    for iii = 1:nn
        # create random vector

        dim = rand(rng, 1:100)
        vMin = -2
        vMax = 2
        x = vMin + rand(rng, dim, 1) * (vMax - vMin)
        t = rand(rng)
        x = [x; t]
        # project matrix onto positive semidefinite cone
        xNew = Projections.secondOrderCone(x)

        # check if relation holds {(t,x) | ||x||_2 <= t}
        @test (norm(xNew[2:end], 2) - xNew[1]) < 1e-14

    end
end

@testset "Test projection onto pos-def cone" begin

  for iii = 1:nn

    # create random matrix
    dim = rand(rng, 1:100)
    A = rand(rng, dim, dim)
    A = full(Symmetric(A))
    # project matrix onto positive semidefinite cone
    a = vec(A)
    X = Projections.sdcone(a)

    # check positive semidefiniteness
    @test minimum(eig(reshape(X, dim, dim))[1]) >= -1e-10

  end

end
nothing


text:
struct GraphConv{A <: AbstractMatrix, B}
    weight1::A
    weight2::A
    bias::B
    σ
end

function GraphConv(ch::Pair{Int, Int}, σ=identity; init=Flux.glorot_uniform, bias::Bool=true, T::DataType=Float32)
    in, out = ch
    W1 = init(out, in)
    W2 = init(out, in)
    b = bias ? T.(init(ch[2])) : zeros(T, ch[2])
    return GraphConv(W1, W2, b, σ)
end

Flux.@functor GraphConv

function (g::GraphConv)(fgs::BatchedFeaturedGraph{Float32})
    A, X = fgs.graph, fgs.nf
    return BatchedFeaturedGraph{Float32}(
        fgs.graph;
        nf = g.σ.(g.weight1 ⊠ X .+ g.weight2 ⊠ X ⊠ A .+ g.bias),
        ef = fgs.ef,
        gf = fgs.gf
    )
end

function (g::GraphConv)(fg::FeaturedGraph)
    A, X = fg.graph, fg.nf
    return FeaturedGraph(
        fg.graph,
        g.σ.(g.weight1 * X .+ g.weight2 * X * A .+ g.bias),
        fg.ef,
        fg.gf,
        fg.directed
    )
end


text:
export eye

"""
`eye(T,n)` returns an `n`-by-`n` identity matrix whose entries
are numbers of type `T` (`Float64` by default).
"""
function eye(T,n)
    A = zeros(T,n,n)
    for j=1:n
        A[j,j] = 1
    end
    return A
end

eye(n) = eye(Float64,n)


text:
using Plots
using Roots

function plot_dielectric_functions()
    χ0_metal(q, ktf=1) = -ktf^2 / (4π)
    χ0_dielectric(q, εr, C₀=1-εr, ktf=1) = C₀*q^2 / (4π * (1 - C₀*q^2/ktf^2))
    ε(χ0, q) = (1 - (4π)/q^2 * χ0(q))

    kerker(q, ktf=1) = q^2 / (ktf^2 + q^2)

    xs = range(1e-5, 1.5, length=10000)
    p = plot(xlims=(0, Inf), )
    plot!(p, x -> ε(χ0_metal, x),                    xs, label="Al",   ls=:solid, lw=2)
    plot!(p, x -> ε(q -> χ0_dielectric(q, 14.9), x), xs, label="GaAs", ls=:dash, lw=2)
    plot!(p, x -> ε(q -> χ0_dielectric(q,  2.4), x), xs, label="SiO₂", ls=:dashdot, lw=2)
    ylims!(p, (0, 16))
    xaxis!(p, "q")
    yaxis!(p, "ε(q)")

    p
end

function intersection_kerker_preconditioned()
    # Find intersection between Kerker-preconditioned and not preconditioned

    function ratio(qmin; qmax=100, prec=q -> 1)
        a = prec(qmin) * ε(q -> χ0_dielectric(q, 14.9), qmin)
        b = prec(qmax) * ε(q -> χ0_dielectric(q, 14.9), qmax)
        max(a, b) / min(a, b)
    end
    p = plot(xlims=(0, Inf), )
    plot!(p, q -> ratio(q),              xs, ls=:solid, label="None")
    plot!(p, q -> ratio(q, prec=kerker), xs, ls=:dash,  label="Kerker")
    ylims!(p, (0, 16))
    xaxis!(p, "q")
    yaxis!(p, "κ(q)")
    # savefig(p, "prec_dielectric_functions.pdf")

    qcut = Roots.find_zero(x -> ratio(x) - ratio(x, prec=kerker), [0.05, 0.1])
    println("qcut = $qcut")
    println("L    = $(2π / qcut)")
end


text:
##
## Data Structures
##

# Each separator is associated with a tree like
#
#       1
#       |
#       1
#     / | \
#    1  2  3
#   / \ |  |
#   1 2 3  4
#
# The ids are increasing left-right

# A node a the tree
mutable struct NodeIT
    i::Int64                        # Id
    c::Array{NodeIT,1}              # Childrens
    p::Union{Nothing, NodeIT}       # Parent
    function NodeIT(i::Int64)
        this = new()
        this.i = i
        this.c = NodeIT[]
        this.p = nothing
        this
    end
end
Base.show(io::IO, n::NodeIT) = print(io,n.i)

# The tree
mutable struct IT
    nit::Array{Array{NodeIT,1},1}
    function IT(i::Int64)
        this = new()
        this.nit = [ [ NodeIT(i) ] ]
        this
    end

    function IT()
        this = new()
        this.nit = [ [] ]
        this
    end
end

IP = NamedTuple{(:p, :l, :r),Tuple{Int64,Int64,Int64}} # Parent/Left/Right

function lt_ip(x::IP, y::IP)
    return      (x.p <  y.p) ||
              ( (x.p == y.p) && (x.l <  y.l) ) ||
              ( (x.p == y.p) && (x.l == y.l) && (x.r <  y.r) )
end

ord = Base.ord(isless,identity,false,Base.Forward)

const IP0 = (p=1,l=1,r=1)

# Partition the graph represented by the sparse matrix A
# from level 1 to maxLevel
#
# This performs a 'modified' nested dissection ordering, i.e., a ND ordering but with 'matching' partitions.
#
# By convention, level=1 are the leaf separators and level=maxLevel is the root separator.
#
# Each ND separator is further divided in a hierarchy of clusters
# Separator at level 1 <= l <= maxLevel consists of a l-levels hierarchy of nested clusters.
# Hierarchy is encoded into a integer-tree representing the nested hierarchy
#
# For instance, a given separator's (lvl=3, sep=?) dofs (degrees of freedom, i.e., original unknowns) are numbered
#       subseps[dofs[lvl][sep]] = [1 1 2 3 3 4 5 5 5] (in practice, they are not necessary sorted)
# with a tree stored in
#       hrch[lvl][sep] =        1
#                             /   \
#                            1     2
#                          / | \   | \
#                         1  2  3  4  5
# indicating a 3-levels hierarchy and the way the merging should be performed
# In particular, the hierarchy is
#                       [1 1 2 3 3 4 5 5 5] 5 partitions
#                       [1 1 1 1 1 2 2 2 2] 2 partitions
#                       [1 1 1 1 1 1 1 1 1] 1 partition - root
#
# Outputs
#       * subseps       subseps[i] = 1 ... #partitions, indicates which cluster within its ND separator node i belongs to
#       * hrch          hrch[l][s] a l-levels hierarchy for separator (l,s)
#       * dofs          dofs[l][s] the dofs of separator (l,s)
function mnd(A::SparseMatrixCSC{Tv,Ti}, maxLevel ; verbose::Bool=false) where {Tv, Ti}

    @assert maxLevel >= 1
    if verbose
        println("Algebraic partitionning of graph with $(size(A,1)) vertices, $maxLevel levels")
    end

    n = size(A,1)
    parts = Array{IP,1}(undef, n) # (parent, left, right) for each dof
    fill!(parts,IP0)

    # Prepare stuff to separate
    tosep  = [ collect(Cint, 1:n) ]

    # Temp arrays for Metis
    colptr = Vector{Cint}(undef, n+1)
    rowval = Vector{Cint}(undef, nnz(A))
    pids   = Vector{Cint}(undef, n)

    # Output structures
    dofs    = Vector{Vector{Vector{Int64}}}(undef, maxLevel)
    hrch    = Vector{Vector{IT}}(undef, maxLevel)

    # Do the work
    for lvl = 1:maxLevel

        dofs[lvl] = Vector{Vector{Int64}}(undef, 2^(lvl-1))
        hrch[lvl] = Vector{IT}(undef, 2^(lvl-1))
        for sep = 1:2^(lvl-1)
            hrch[lvl][sep] = IT(1)
        end
        tosep_next = Vector{Vector{Int64}}(undef, 2^lvl)

        if verbose println("Level $lvl, partitioning") end

        for (sep, list) in enumerate(tosep)

            # What is our separator -unique- id
            id = 2^(lvl-1) + (sep - 1)

            # Do the ordering, i.e., cut A[list,list]
            if lvl == maxLevel # Leaf level
                pids[1:length(list)] .= 2
            elseif length(list) <= 1 # Really isn't anything to separate
                pids[1:length(list)] .= 2
            else
                fill_nodiag_list_sorted!(A, colptr, rowval, list)
                vertex_sep_bissect_fast!(length(list), colptr, rowval, pids)
            end

            # Update parts, so that each node has a proper (left, right) pair
            for (i, side) in zip(list, pids)
                p = parts[i]
                if side == 2 && p.r == id && p.l == id # Separator, brand new
                    p = (p=1, l=2*id, r=2*id+1)
                elseif side == 0 # In the left interior
                    if p.l == id && p.r == id # Interior, brand new
                        p = (p=1, l=2*id, r=2*id)
                    elseif p.l == id # Upper level separator, on its left
                        p = (p=p.p, l=2*id, r=p.r)
                    elseif p.r == id # Upper level separator, on its right
                        p = (p=p.p, l=p.l, r=2*id)
                    end
                elseif side == 1 # In the right interior
                    if p.l == id && p.r == id # Interior, brand new
                        p = (p=1, l=2*id+1, r=2*id+1)
                    elseif p.l == id # Upper level separator, on its left
                        p = (p=p.p, l=2*id+1, r=p.r)
                    elseif p.r == id # Upper level separator, on its right
                        p = (p=p.p, l=p.l, r=2*id+1)
                    end
                end
                parts[i] = p
            end

            # Find new sep & create initial hierarchy
            mid = findall(i -> parts[i] == (p=1, l=2*id, r=2*id+1), list)
            dofs[lvl][sep] = list[mid]

            # Define new interiors + boundaries
            left  = findall(i -> parts[i].l == 2*id  , list)
            right = findall(i -> parts[i].r == 2*id+1, list)
            tosep_next[2*sep-1] = list[left]
            tosep_next[2*sep  ] = list[right]
        end

        # Go through edges and update hierarchy & parents
        if lvl < maxLevel
            for lvl2 = 1:lvl
                for sep2 = 1:2^(lvl2-1)
                    vals  = parts[dofs[lvl2][sep2]]
                    lr_srt = sortperm(vals, lt=lt_ip) # Sort by parent first - arbitrary after
                    vals_srt = vals[lr_srt]
                    dofs_srt = dofs[lvl2][sep2][lr_srt]
                    u_ptr = sorted_to_ptr(vals_srt)
                    # Create hierarchy
                    it = hrch[lvl2][sep2]
                    push!(it.nit, NodeIT[]) # Create new level
                    if length(dofs_srt) == 0
                        @assert length(it.nit[end-1]) == 1
                        x = NodeIT(1)
                        x.p = it.nit[end-1][1]
                        push!(it.nit[end-1][1].c,x)
                        push!(it.nit[end], x)
                    else
                        for p = 1:length(u_ptr)-1
                            dofs_p = dofs_srt[u_ptr[p]:u_ptr[p+1]-1] # The dofs - this is non-empty by definition of unique
                            old_p  = parts[dofs_p[1]].p # All same since we sorted by parent
                            left   = parts[dofs_p[1]].l # ..
                            right  = parts[dofs_p[1]].r # ..
                            # Update partition
                            for i in dofs_p
                                parts[i] = (p=p, l=left, r=right)
                            end
                            # Create node & update hierarchy
                            x     = NodeIT(p)
                            x.p   = it.nit[end-1][old_p]
                            push!(it.nit[end-1][old_p].c, x)
                            push!(it.nit[end], x)
                        end
                    end
                end
            end
        end

        # Next clusters to consider
        tosep = tosep_next
    end

    subseps = map(x -> x.p, parts)
    reverse!(hrch)
    reverse!(dofs)

    return subseps, hrch, dofs

end


text:
"""
    MonteCarloModel(core, dates, paths)

A `MonteCarloModel` is the result of a simulation of a series of asset prices.

* `core`: a reference `CoreModel`
* `dates`: an `AbstractVector{Date}`
* `paths`: a matrix of the scenario paths: the rows are the scenarios, and the columns are the values at each date in `dates`.
"""
struct MonteCarloModel{C,D,T} <: AbstractModel
    core::C
    dates::D
    paths::Matrix{T}
end

"""
    MonteCarloScenario(core, dates, path)

A `MonteCarloScenario` is a single simulation scenario of a `MonteCarloModel`.

* `core`: a reference `CoreModel`
* `dates`: an `AbstractVector{Date}`
* `paths`: an `AbstractVector` of the values at each date in `dates`.
"""
struct MonteCarloScenario{C,D,S} <: AbstractModel
    core::C
    dates::D
    path::S
end

numeraire(m::Union{MonteCarloModel, MonteCarloScenario}) = numeraire(m.core)
startdate(m::Union{MonteCarloModel, MonteCarloScenario}) = startdate(m.core)

yearfractionto(m::Union{MonteCarloModel, MonteCarloScenario}, dt::Date) =
    yearfractionto(m.core,  dt)

# the value of currency is the same under every scenario
value(m::MonteCarloModel, c::WhenAt{Receive{T}}) where {T} =
    value(m.core,c)
value(m::MonteCarloScenario, c::WhenAt{Receive{T}}) where {T} =
    value(m.core,c)

struct ScenarioIterator{M<:MonteCarloModel}
    m::M
    n::Int
end

"""
    scenarios(m::MonteCarloModel)

Returns an iterator over each `MonteCarloScenario` in `m`.
"""
scenarios(m::MonteCarloModel) = ScenarioIterator(m, size(m.paths, 1))

Base.length(sc::ScenarioIterator) = sc.n
Base.iterate(sc::ScenarioIterator, i::Int=1) = i > sc.n ? nothing :
    (MonteCarloScenario(sc.m.core, sc.m.dates, view(sc.m.paths, i, :)), i+1)

"""
    date2index(m::Union{MonteCarloScenario, MonteCarloModel}, dt::Date)

Returns the index of `dt` in the path(s) of `m`.
"""
function date2index(m::Union{MonteCarloScenario, MonteCarloModel}, dt::Date)
    ii = searchsorted(m.dates, dt)
    isempty(ii) && throw(DomainError(dt))
    return ii[1]
end

index2date(m::Union{MonteCarloScenario, MonteCarloModel}, i) = m.dates[i]

function forwardprice(m::MonteCarloScenario, s::SingleStock, dt::Date)
    valueat(m,s,date2index(m,dt))
end
function forwardprice(m::MonteCarloModel, s::SingleStock, dt::Date)
    mean(forwardprice(ms, s, dt) for ms in scenarios(m))
end

valueat(m::MonteCarloScenario, ::SingleStock, i::Int) = m.path[i]
valueat(m::MonteCarloScenario, ::SingleStock, i::Int, ::Type{Dual}) =
    Dual(m.path[i], 1)

# fallback
function value(m::MonteCarloModel, c::WhenAt)
    N = date2index(m, maturitydate(c))
    discount(m.core.yieldcurve, maturitydate(c)) * mean(valueat(ms, c.c, N) for ms in scenarios(m))
end

"""
    montecarlo(m::GeomBMModel, dates, npaths)

Sample `npaths` Monte Carlo paths of the model `m`, at time `dates`.
"""
function montecarlo(m::GeomBMModel{CoreModel{T,R,Q}, V}, dates::StepRange{Date}, npaths::Integer) where {T,R,Q,V}
    σ = m.volatility
    S = typeof(m.core.yieldcurve.rate)
    Xt = Array{promote_type(T,V,S)}(undef, length(dates), npaths)
    Δt = yearfraction(daycount(m.core.yieldcurve), step(dates))
    df = discount(m.core.carrycurve, Δt) / discount(m.core.yieldcurve, Δt)
    for i = 1:npaths
        x = value(m, SingleStock())
        for (j, dt) in enumerate(dates)
            if j == 1
                Δt1 = yearfraction(daycount(m.core.yieldcurve), startdate(m), first(dates))
                df1 = Δt1 == 0 ? 1.0 : discount(m.core.carrycurve, Δt1) / discount(m.core.yieldcurve, Δt1)
                x *= df1 * exp(-σ^2*Δt1/2 + σ*sqrt(Δt1)*randn())
            else
                x *= df * exp(-σ^2*Δt/2 + σ*sqrt(Δt)*randn())
            end
            Xt[j,i] = x
        end
    end
    MonteCarloModel(m.core, dates, copy(transpose(Xt)))
end

function value(m::GeomBMModel, c::Contract, ::Type{MonteCarloModel}, dates::StepRange{Date}, npaths::Integer)
    mcm = montecarlo(m, dates, npaths)
    value(mcm, c)
end

value(m::GeomBMModel, c::Contract, ::Type{MonteCarloModel}, npaths::Integer) =
    value(m, c, MonteCarloModel, startdate(m):Day(1):maturitydate(c), npaths)


text:
# matrix-form of spin-1/2 single spin operators,
#   in Zeeman basis.

### for understanding.
function getsingleI⁺()::Matrix{Float64}
    return [0 1; 0 0]
end

function getsingleI⁻()::Matrix{Float64}
    return [0 0; 1 0]
end

function getsingleIy()::Matrix{Complex{Float64}}
    return 1/(2*im) .* [0 1; -1 0]
end

function getsinglehalfId()::Matrix{Float64}
    return 0.5 .* [1 0; 0 1]
end

### actual use.

function getsingleIx()::Matrix{Float64}
    return 0.5 .* [0 1; 1 0]
end

function getsingleIynoim()::Matrix{Float64}
    return 0.5 .* [0 1; -1 0]
end

function getsingleIz()::Matrix{Float64}
    return 0.5 .* [1 0; 0 -1]
end

function getsingleId()::Matrix{Float64}
    return [1 0; 0 1]
end

#### product operators.

function getmultiIz(   Id::Matrix{T},
                        j::Int,
                        N::Int) where T <: Real

    # prepare base operators.
    # prepare base operators.
    operators = Vector{Matrix{T}}(undef, N)
    fill!(operators, Id)

    operators[j] = getsingleIz()

    # get product operator.
    #out = 2^(N-1) .* Kronecker.kronecker(operators...)
    out = Kronecker.kronecker(operators...)

    return out
end

function getmultiIjIk(   Id,
                        Ix,
                        Iy_no_im,
                        Iz,
                        j::Int,
                        k::Int,
                        N::Int,
                        c::T) where T <: Real

    ### x.
    # prepare base operators.
    x_operators = Vector{Matrix{T}}(undef, N)
    fill!(x_operators, Id)

    x_operators[j] = c .* Ix
    x_operators[k] = Ix

    # get product operator.
    x_out = Kronecker.kronecker(x_operators...)

    ### z.
    # prepare base operators.
    z_operators = Vector{Matrix{T}}(undef, N)
    fill!(z_operators, Id)

    z_operators[j] = c .* Iz
    z_operators[k] = Iz

    # get product operator.
    z_out = Kronecker.kronecker(z_operators...)

    ### y.
    # prepare base operators.
    y_operators = Vector{Matrix{T}}(undef, N)
    fill!(y_operators, Id)

    y_operators[j] = c .* Iy_no_im
    y_operators[k] = Iy_no_im

    # get product operator.
    y_out = -Kronecker.kronecker(y_operators...)

    out = x_out + y_out + z_out

    return out
end

function getmultiIjIk0(   Id::Matrix{T},
                        Ix::Matrix{T},
                        Iy_no_im::Matrix{T},
                        Iz::Matrix{T},
                        j::Int,
                        k::Int,
                        N::Int) where T <: Real

    ### x.
    # prepare base operators.
    x_operators = Vector{Matrix{T}}(undef, N)
    fill!(x_operators, Id)

    x_operators[j] = Ix
    x_operators[k] = Ix

    # get product operator.
    x_out = Kronecker.kronecker(x_operators...)

    ### z.
    # prepare base operators.
    z_operators = Vector{Matrix{T}}(undef, N)
    fill!(z_operators, Id)

    z_operators[j] = Iz
    z_operators[k] = Iz

    # get product operator.
    z_out = Kronecker.kronecker(z_operators...)

    ### y.
    # prepare base operators.
    y_operators = Vector{Matrix{T}}(undef, N)
    fill!(y_operators, Id)

    y_operators[j] = Iy_no_im
    y_operators[k] = Iy_no_im

    # get product operator.
    y_out = -Kronecker.kronecker(y_operators...)

    out = x_out + y_out + z_out

    return out
end


text:
"""
This defines abstract method type
    and implements functions for abstract mehtod.
"""
abstract type AbstractMethod end

# This returns the internal JuMP.Model in BundleModel.
get_jump_model(method::AbstractMethod)::JuMP.Model = get_model(get_model(method))

set_optimizer(method::AbstractMethod, optimizer) = set_optimizer(get_model(method), optimizer)

# This returns solution.
function get_solution(method::AbstractMethod) end

# This returns objective value.
get_objective_value(method::AbstractMethod)::Float64 = Inf

# This sets the termination tolerance.
function set_bundle_tolerance!(method::AbstractMethod, tol::Float64) end

# This sets the termination wrt objective value
set_obj_limit(method::AbstractMethod, bound::Float64) = set_parameter(method.params, "obj_limit", bound)

# This builds the initial bundle model.
function build_bundle_model!(method::AbstractMethod)
    add_variables!(method)
    add_objective_function!(method)
    add_constraints!(method)
end

# This creates variables to the bundle model.
function add_variables!(method::AbstractMethod)
    bundle = get_model(method)
    model = get_model(bundle)
    @variable(model, x[i=1:bundle.n])
    @variable(model, θ[j=1:bundle.ncuts_per_iter])
end

# This creates an objective function to the bundle model.
function add_objective_function!(method::AbstractMethod)
    bundle = get_model(method)
    model = get_model(bundle)
    θ = model[:θ]
    @objective(model, Min, sum(θ[j] for j = 1:bundle.ncuts_per_iter))
end

# This creates constraints to the bundle model.
function add_constraints!(method::AbstractMethod) end

# This implements the algorithmic steps.
function run!(method::AbstractMethod)
    add_initial_bundles!(method)
    update_iteration!(method)
    solve_bundle_model!(method)
    display_info!(method)
    while !termination_test(method)
        evaluate_functions!(method)
        update_bundles!(method)
        update_iteration!(method)
        solve_bundle_model!(method)
        display_info!(method)
    end
end

# This implements adding initial bundle to the model.
function add_initial_bundles!(method::AbstractMethod)
    evaluate_functions!(method)
    add_bundles!(method)
end

# This implements the bundle model solution step.
function solve_bundle_model!(method::AbstractMethod)
    model = get_model(method)
    solve_model!(model)
    collect_model_solution!(method)
end

# This may collect solutions from the bundle model.
function collect_model_solution!(method::AbstractMethod) end

# Should the method terminate?
termination_test(method::AbstractMethod) = true

# This should call BundleModel.evaluate_f function
# to get any necessary information (e.g., function value and subgradient).
function evaluate_functions!(method::AbstractMethod) end

# This updates the bundle pool by removing and/or adding bundle objects.
function update_bundles!(method::AbstractMethod)
    purge_bundles!(method)
    add_bundles!(method)
end
function purge_bundles!(method::AbstractMethod) end
function add_bundles!(method::AbstractMethod) end

# This updates any information for the current iteration (e.g., iteration count).
function update_iteration!(method::AbstractMethod) end

# This displays iteration information.
function display_info!(method::AbstractMethod) end


text:
export iauPv2s
"""
Convert position/velocity from Cartesian to spherical coordinates.

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

Status:  vector/matrix support function.

Given:
    pv       double[2][3]  pv-vector

Returned:
    theta    double        longitude angle (radians)
    phi      double        latitude angle (radians)
    r        double        radial distance
    td       double        rate of change of theta
    pd       double        rate of change of phi
    rd       double        rate of change of r

Notes:

    1. If the position part of pv is null, theta, phi, td and pd
    are indeterminate.  This is handled by extrapolating the
    position through unit time by using the velocity part of
    pv.  This moves the origin without changing the direction
    of the velocity component.  If the position and velocity
    components of pv are both null, zeroes are returned for all
    six results.

    2. If the position is a pole, theta, td and pd are indeterminate.
    In such cases zeroes are returned for all three.

This revision:  2013 June 18

SOFA release 2018-01-30

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

# void iauPv2s(double pv[2][3],
#              double *theta, double *phi, double *r,
#              double *td, double *pd, double *rd)

function iauPv2s(pv::AbstractMatrix{<:Real})
    # Preallocate return values
    ref_theta = Ref{Float64}(0.0)
    ref_phi = Ref{Float64}(0.0)
    ref_r = Ref{Float64}(0.0)
    ref_td = Ref{Float64}(0.0)
    ref_pd = Ref{Float64}(0.0)
    ref_rd = Ref{Float64}(0.0)

    ccall((:iauPv2s, libsofa_c), Cvoid,
          (Ptr{Cdouble},
          Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
          Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
           convert(Array{Float64, 2}, pv'),
           ref_theta, ref_phi, ref_r, ref_td, ref_pd, ref_rd)

    return ref_theta[], ref_phi[], ref_r[], ref_td[], ref_pd[], ref_rd[]
 end


text:
using MOCNeutronTransport
@testset "Abaqus" begin
    @testset "c5g7 pin - triangles" begin
        filepath = "./mesh/mesh_files/c5g7_UO2_pin.inp"
        ref_points = [Point_2D(0.0, 0.0),
                      Point_2D(1.26, 0.0),
                      Point_2D(0.0, 1.26),
                      Point_2D(1.26, 1.26),
                      Point_2D(1.17, 0.63),
                      Point_2D(0.63, 0.0),
                      Point_2D(0.0, 0.63),
                      Point_2D(1.26, 0.63),
                      Point_2D(0.63, 1.26),
                      Point_2D(0.96668449300372, 1.0521890005327),
                      Point_2D(0.50983869566359, 1.1564610725782),
                      Point_2D(0.14347681133269, 0.86429721912348),
                      Point_2D(0.14347681133269, 0.39570278087652),
                      Point_2D(0.50983869566359, 0.10353892742181),
                      Point_2D(0.96668449300372, 0.20781099946726),
                      Point_2D(0.63, 0.63)]
        ref_faces = [
                    SVector(5, 12, 16, 11),
                    SVector(5, 13, 16, 12),
                    SVector(5, 14, 16, 13),
                    SVector(5, 11, 16, 10),
                    SVector(5, 10, 16, 5),
                    SVector(5, 15, 16, 14),
                    SVector(5, 5, 16, 15),
                    SVector(5, 7, 1, 13),
                    SVector(5, 7, 12, 3),
                    SVector(5, 1, 6, 14),
                    SVector(5, 1, 14, 13),
                    SVector(5, 6, 2, 15),
                    SVector(5, 15, 2, 8),
                    SVector(5, 3, 11, 9),
                    SVector(5, 3, 12, 11),
                    SVector(5, 10, 8, 4),
                    SVector(5, 9, 10, 4),
                    SVector(5, 10, 5, 8),
                    SVector(5, 15, 8, 5),
                    SVector(5, 14, 6, 15),
                    SVector(5, 7, 13, 12),
                    SVector(5, 11, 10, 9)
                   ]
        ref_face_sets = Dict{String, Set{Int64}}()
        ref_face_sets["MATERIAL_UO2-3.3"] = Set([1, 2, 3, 4, 5, 6, 7])
        ref_face_sets["GRID_L1_1_1"] = Set([1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
                                            11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22])
        ref_face_sets["MATERIAL_WATER"] = Set([8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
                                               18, 19, 20, 21, 22])

        # read_abaqus
        mesh = read_abaqus_2d(filepath)
        @test mesh.points == ref_points
        @test mesh.faces == ref_faces
        @test mesh.name == "c5g7_UO2_pin"
        @test mesh.face_sets == ref_face_sets
    end
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
using StokesParticles, Test

@test_throws DimensionMismatch LineWall([Point2D(0.0,0.0)],0.1)

ArcWall(Point2D(1.0,0.0),Point2D(0.0,0.0),Point2D(0.0,1.0),0.1)

@test_throws MethodError ArcWall(Point2D(1.0,0.0),Point2D(0.0,0.0),Point2D(0.0,2.0),0.1)


text:
function ECCO_path_etc(sol0::String,calc::String,nam::String)
    sol="ECCOv4"*sol0*"_analysis"

    if sol0=="r1"||sol0=="r2"
        nt=240
    elseif sol0=="r3"
        nt=288
    elseif sol0=="r4"
        nt=312
    elseif sol0=="r5"
        nt=336
    end

    if sol0!=="r5"
        pth_in="ECCOv4"*sol0*"/nctiles_monthly/"
    else
        pth_in="ECCOv4"*sol0*"/diags/"
    end
    list_steps=list_time_steps(pth_in)

    pth_out=joinpath("ECCO_diags",sol)

    if sum(calc.==("overturn","MHT","trsp"))==0
        pth_tmp=joinpath(pth_out,nam*"_"*calc)
    else
        pth_tmp=joinpath(pth_out,calc)
    end

    return pth_in,pth_out,pth_tmp,nt,list_steps
end

#STATE/state_3d_set1.0000241020.meta
#    'THETA   ' 'SALT    ' 'DRHODR  '
#TRSP/trsp_3d_set1.0000241020.meta
#    'UVELMASS' 'VVELMASS' 'WVELMASS' 'GM_PsiX ' 'GM_PsiY '
#TRSP/trsp_3d_set3.0000241020.meta
#    'DFxE_TH ' 'DFyE_TH ' 'ADVx_TH ' 'ADVy_TH ' 'DFxE_SLT' 'DFyE_SLT' 'ADVx_SLT' 'ADVy_SLT'

function list_time_steps(pth_in)
    if isdir(joinpath(pth_in,"STATE"))
        list=readdir(joinpath(pth_in,"STATE"))
        list=list[findall(occursin.(Ref("state_3d_set1"),list))]
        list=list[findall(occursin.(Ref("data"),list))]
        [list[i][15:end] for i in 1:length(list)]
    else
        list=[]
    end
    return list
end

function get_nr_nt(pth_in,nam)
    nct_path=joinpath(pth_in,nam)
    lst=readdir(nct_path)
    lst=lst[findall(occursin.(Ref(".nc"),lst))]
    fil1=joinpath(nct_path,lst[1])
    ds=NCTiles.NCDatasets.Dataset(fil1)
    siz=size(ds[nam])
    siz[end-1],siz[end]
end

nansum(x) = sum(filter(!isnan,x))
nansum(x,y) = mapslices(nansum,x,dims=y)

function GridLoad_Main()

    γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
    nr=50
    XC=GridLoadVar("XC",γ)
    YC=GridLoadVar("YC",γ)
    RAC=GridLoadVar("RAC",γ)
    hFacC=GridLoadVar("hFacC",γ)
    hFacW=GridLoadVar("hFacW",γ)
    hFacS=GridLoadVar("hFacS",γ)
    DRF=GridLoadVar("DRF",γ)
    DXC=GridLoadVar("DXC",γ)
    DYC=GridLoadVar("DYC",γ)
    DXG=GridLoadVar("DXG",γ)
    DYG=GridLoadVar("DYG",γ)

    mskC=hFacC./hFacC

    tmp=hFacC./hFacC
    tmp=[nansum(tmp[i,j].*RAC[i]) for j in 1:nr, i in eachindex(RAC)]
    tot_RAC=nansum(tmp,2)

    tmp=[nansum(hFacC[i,j].*RAC[i].*DRF[j]) for j in 1:nr, i in eachindex(RAC)]
    tot_VOL=nansum(tmp,2)

    Γ=(XC=XC,YC=YC,RAC=RAC,DXG=DXG,DYG=DYG,DXC=DXC,DYC=DYC,DRF=DRF,
        hFacC=hFacC,hFacW=hFacW,hFacS=hFacS,
        mskC=mskC,tot_RAC=tot_RAC,tot_VOL=tot_VOL)

    return γ,Γ
end

## generic read function

function read_monthly(sol,nam,t,list_steps)
    if nam=="SSH"
        read_monthly_SSH(sol,nam,t,list_steps)
    elseif nam=="MHT"
        read_monthly_MHT(sol,nam,t,list_steps)
    elseif nam=="BSF"
        read_monthly_BSF(sol,nam,t,list_steps)
    else
        read_monthly_default(sol,nam,t,list_steps)
    end
end

###

function read_monthly_SSH(sol,nam,t,list_steps)
    ETAN=read_monthly_default(sol,"ETAN",t,list_steps)
    sIceLoad=read_monthly_default(sol,"sIceLoad",t,list_steps)
    (ETAN+sIceLoad/1029.0)*Γ.mskC[:,1]
end

function read_monthly_MHT(sol,nam,t,list_steps)
    U=read_monthly(sol,"ADVx_TH",t,list_steps)
    V=read_monthly(sol,"ADVy_TH",t,list_steps)
    U=U+read_monthly(sol,"DFxE_TH",t,list_steps)
    V=V+read_monthly(sol,"DFyE_TH",t,list_steps)

    [U[i][findall(isnan.(U[i]))].=0.0 for i in eachindex(U)]
    [V[i][findall(isnan.(V[i]))].=0.0 for i in eachindex(V)]

    Tx=0.0*U[:,1]
    Ty=0.0*V[:,1]
    [Tx=Tx+U[:,z] for z=1:nr]
    [Ty=Ty+V[:,z] for z=1:nr]

    return Tx,Ty
end

function read_monthly_BSF(sol,nam,t,list_steps)
    U=read_monthly(sol,"UVELMASS",t,list_steps)
    V=read_monthly(sol,"VVELMASS",t,list_steps)
    (Utr,Vtr)=UVtoTransport(U,V,Γ)

    nz=size(Γ.hFacC,2)
    μ=Γ.mskC[:,1]
    Tx=0.0*Utr[:,1]
	Ty=0.0*Vtr[:,1]
	for z=1:nz
		Tx=Tx+Utr[:,z]
		Ty=Ty+Vtr[:,z]
	end

    #convergence & land mask
    TrspCon=μ.*convergence(Tx,Ty)

    #scalar potential
    TrspPot=ScalarPotential(TrspCon)

    #Divergent transport component
    (TxD,TyD)=gradient(TrspPot,Γ)
    TxD=TxD.*Γ.DXC
    TyD=TyD.*Γ.DYC

    #Rotational transport component
    TxR = Tx-TxD
    TyR = Ty-TyD

    #vector Potential
    TrspPsi=VectorPotential(TxR,TyR,Γ)

    return TrspPsi
end

###

function read_monthly_default(sol,nam,t,list)
    var_list3d=("THETA","SALT","UVELMASS","VVELMASS",
                "ADVx_TH","ADVy_TH","DFxE_TH","DFyE_TH")
    mdsio_list3d=("STATE/state_3d_set1","STATE/state_3d_set1",
        "TRSP/trsp_3d_set1","TRSP/trsp_3d_set1","TRSP/trsp_3d_set2",
        "TRSP/trsp_3d_set2","TRSP/trsp_3d_set2","TRSP/trsp_3d_set2")

    var_list2d=("MXLDEPTH","SIarea","sIceLoad","ETAN")
    mdsio_list2d=("STATE/state_2d_set1","STATE/state_2d_set1",
                  "STATE/state_2d_set1","STATE/state_2d_set1")

    if (sol=="ECCOv4r1_analysis")||(sol=="ECCOv4r2_analysis")||(sol=="ECCOv4r3_analysis")
        nct_path=joinpath(pth_in,nam)
        if sum(var_list3d.==nam)==1
            tmp=read_nctiles(nct_path,nam,γ,I=(:,:,:,t))
        else
            tmp=read_nctiles(nct_path,nam,γ,I=(:,:,t))
        end
    elseif (sol=="ECCOv4r4_analysis")
        y0=Int(floor((t-1)/12))+1992
        m0=mod1(t,12)
        nct_path=joinpath(pth_in,nam,string(y0))
        m0<10 ? fil=nam*"_$(y0)_0$(m0).nc" : fil=nam*"_$(y0)_$(m0).nc"
        tmp0=NCTiles.NCDatasets.Dataset(joinpath(nct_path,fil))[nam]
        til0=Tiles(γ,90,90)
        if sum(var_list3d.==nam)==1
            tmp=MeshArray(γ,γ.ioPrec,nr)
            for i in 1:13, k in 1:50
              ff=til0[i].face
              ii=collect(til0[i].i)
              jj=collect(til0[i].j)
              tmp[ff,k][ii,jj]=tmp0[:,:,i,k,1]
            end
            tmp
        else
            tmp=MeshArray(γ,γ.ioPrec)
            for i in 1:13
              ff=til0[i].face
              ii=collect(til0[i].i)
              jj=collect(til0[i].j)
              tmp[ff][ii,jj]=tmp0[:,:,i,1]
            end
            tmp
        end
    else
      if !isempty(findall(var_list3d.==nam))
        ii=findall(var_list3d.==nam)[1];
        fil=mdsio_list3d[ii]
        meta=read_meta(joinpath(pth_in,fil*".0000241020.meta"))
        kk=findall(vec(meta.fldList).==nam)[1]
        tmp=read_mdsio(joinpath(pth_in,fil*list[t][14:end]))[:,:,:,kk]
        tmp=Γ.mskC*read(tmp,γ)
      else
        ii=findall(var_list2d.==nam)[1];
        fil=mdsio_list2d[ii]
        meta=read_meta(joinpath(pth_in,fil*".0000241020.meta"))
        kk=findall(vec(meta.fldList).==nam)[1]
        tmp=read_mdsio(joinpath(pth_in,fil*list[t][14:end]))[:,:,kk]
        tmp=Γ.mskC[:,1]*read(tmp,Γ.XC)
      end
    end
end

import Base:push!

function push!(allcalc::Vector{String},allnam::Vector{String},allkk::Vector{Int};
    calc="unknown",nam="unknown",kk=1)
    push!(allcalc,calc)
    push!(allnam,nam)
    push!(allkk,kk)
end

function ECCO_standard_list_toml(fil)

    allcalc=String[]
    allnam=String[]
    allkk=Int[]

    push!(allcalc,allnam,allkk;calc="trsp")
    push!(allcalc,allnam,allkk;calc="MHT")
    push!(allcalc,allnam,allkk;calc="zonmean2d",nam="SIarea")
    push!(allcalc,allnam,allkk;calc="zonmean2d",nam="MXLDEPTH")
    push!(allcalc,allnam,allkk;calc="zonmean2d",nam="SSH")
    push!(allcalc,allnam,allkk;calc="zonmean",nam="THETA")
    push!(allcalc,allnam,allkk;calc="glo2d",nam="THETA")
    push!(allcalc,allnam,allkk;calc="glo3d",nam="THETA")
    push!(allcalc,allnam,allkk;calc="zonmean",nam="SALT")
    push!(allcalc,allnam,allkk;calc="glo2d",nam="SALT")
    push!(allcalc,allnam,allkk;calc="glo3d",nam="SALT")
    push!(allcalc,allnam,allkk;calc="overturn")
    [push!(allcalc,allnam,allkk;calc="clim",nam="THETA",kk=kk) for kk in [1 10 20 29 38 44]]
    [push!(allcalc,allnam,allkk;calc="clim",nam="SALT",kk=kk) for kk in [1 10 20 29 38 44]]
    push!(allcalc,allnam,allkk;calc="clim",nam="BSF")
    push!(allcalc,allnam,allkk;calc="clim",nam="SSH")
    push!(allcalc,allnam,allkk;calc="clim",nam="MXLDEPTH")
    push!(allcalc,allnam,allkk;calc="clim",nam="SIarea")

    tmp1=Dict("calc"=>allcalc,"nam"=>allnam,"kk"=>allkk)
    open(fil, "w") do io
        TOML.print(io, tmp1)
    end

end

##

function ECCO_transport_lines()
    lonPairs=[]
    latPairs=[]
    namPairs=[]

    push!(lonPairs,[-173 -164]); push!(latPairs,[65.5 65.5]); push!(namPairs,"Bering Strait");
    push!(lonPairs,[-5 -5]); push!(latPairs,[34 40]); push!(namPairs,"Gibraltar");
    push!(lonPairs,[-81 -77]); push!(latPairs,[28 26]); push!(namPairs,"Florida Strait");
    push!(lonPairs,[-81 -79]); push!(latPairs,[28 22]); push!(namPairs,"Florida Strait W1");
    push!(lonPairs,[-76 -76]); push!(latPairs,[21 8]); push!(namPairs,"Florida Strait S1");
    push!(lonPairs,[-77 -77]); push!(latPairs,[26 24]); push!(namPairs,"Florida Strait E1");
    push!(lonPairs,[-77 -77]); push!(latPairs,[24 22]); push!(namPairs,"Florida Strait E2");
    push!(lonPairs,[-65 -50]); push!(latPairs,[66 66]); push!(namPairs,"Davis Strait");
    push!(lonPairs,[-35 -20]); push!(latPairs,[67 65]); push!(namPairs,"Denmark Strait");
    push!(lonPairs,[-16 -7]); push!(latPairs,[65 62.5]); push!(namPairs,"Iceland Faroe");
    push!(lonPairs,[-6.5 -4]); push!(latPairs,[62.5 57]); push!(namPairs,"Faroe Scotland");
    push!(lonPairs,[-4 8]); push!(latPairs,[57 62]); push!(namPairs,"Scotland Norway");
    push!(lonPairs,[-68 -63]); push!(latPairs,[-54 -66]); push!(namPairs,"Drake Passage");
    push!(lonPairs,[103 103]); push!(latPairs,[4 -1]); push!(namPairs,"Indonesia W1");
    push!(lonPairs,[104 109]); push!(latPairs,[-3 -8]); push!(namPairs,"Indonesia W2");
    push!(lonPairs,[113 118]); push!(latPairs,[-8.5 -8.5]); push!(namPairs,"Indonesia W3");
    push!(lonPairs,[118 127 ]); push!(latPairs,[-8.5 -15]); push!(namPairs,"Indonesia W4");
    push!(lonPairs,[127 127]); push!(latPairs,[-25 -68]); push!(namPairs,"Australia Antarctica");
    push!(lonPairs,[38 46]); push!(latPairs,[-10 -22]); push!(namPairs,"Madagascar Channel");
    push!(lonPairs,[46 46]); push!(latPairs,[-22 -69]); push!(namPairs,"Madagascar Antarctica");
    push!(lonPairs,[20 20]); push!(latPairs,[-30 -69.5]); push!(namPairs,"South Africa Antarctica");
    push!(lonPairs,[-76 -72]); push!(latPairs,[21 18.5]); push!(namPairs,"Florida Strait E3");
    push!(lonPairs,[-72 -72]); push!(latPairs,[18.5 10]); push!(namPairs,"Florida Strait E4");

    lonPairs,latPairs,namPairs
end

function ECCO_transport_lines(pth_trsp)
    mkdir(pth_trsp)
    lonPairs,latPairs,namPairs=ECCO_transport_lines()
    for ii in 1:length(lonPairs)
        lons=Float64.(lonPairs[ii])
        lats=Float64.(latPairs[ii])
        name=namPairs[ii]
        Trsct=Transect(name,lons,lats,Γ)
        jldsave(joinpath(pth_trsp,"$(Trsct.name).jld2"),
            tabC=Trsct.tabC,tabW=Trsct.tabW,tabS=Trsct.tabS);
    end
    return true
end

function reload_transport_lines(pth_trsp)
    list_trsp=readdir(pth_trsp)
    ntr=length(list_trsp)
    TR=[load(joinpath(pth_trsp,list_trsp[itr])) for itr in 1:ntr]
    return list_trsp,MeshArrays.Dict_to_NamedTuple.(TR),ntr
end


text:
# Tests for Revenue General Direct Approach Models
@testset "RevenueRussellDEAModel" begin

    # Test using Book data
    X = [1; 1; 1; 1; 1; 1; 1; 1];
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5];
    P = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

    # Revenue GDA
    revenuegda = dearevenuegda(X, Y, P, :ERG)

    @test typeof(revenuegda) == RevenueGDADEAModel

    @test nobs(revenuegda) == 8
    @test ninputs(revenuegda) == 1
    @test noutputs(revenuegda) == 2
    @test ismonetary(revenuegda) == false

    @test efficiency(revenuegda) == efficiency(revenuegda, :Economic)
    @test efficiency(revenuegda, :Economic)   ≈ [0.0; 0.25; 0.25; 13/15; 4/3; 1.0; 11/24; 2.5]
    @test efficiency(revenuegda, :Technical)  ≈ [0.0; 0.0; 0.0; 13/15; 4/3; 0.5; 11/24; 37/18]
    @test efficiency(revenuegda, :Allocative) ≈ [0.0; 0.25; 0.25; 0.0; 0.0; 0.5; 0.0; 4/9]
    @test normfactor(revenuegda) ≈ [14.0; 8.0; 8.0; 90/13; 6.0; 4.0; 96/11; 3.0] atol = 1e-5

    # With monetary = true
    revenuegdamon = dearevenuegda(X, Y, P, :ERG, monetary = true)

    @test ismonetary(revenuegdamon) == true

    @test efficiency(revenuegdamon) ≈ efficiency(revenuegda) .* normfactor(revenuegda)
    @test efficiency(revenuegdamon, :Technical) ≈ efficiency(revenuegda, :Technical) .* normfactor(revenuegda)
    @test efficiency(revenuegdamon, :Allocative) ≈ efficiency(revenuegda, :Allocative) .* normfactor(revenuegda)

    # With CRS (Same result in this example)
    revenugdacrs = dearevenuegda(X, Y, P, :ERG, rts = :CRS)

    @test efficiency(revenugdacrs, :Economic)   ≈ [0.0; 0.25; 0.25; 13/15; 4/3; 1.0; 11/24; 2.5]
    @test efficiency(revenugdacrs, :Technical)  ≈ [0.0; 0.0; 0.0; 13/15; 4/3; 0.5; 11/24; 37/18]
    @test efficiency(revenugdacrs, :Allocative) ≈ [0.0; 0.25; 0.25; 0.0; 0.0; 0.5; 0.0; 4/9]

    # Print
    show(IOBuffer(), revenuegda)

    # Test errors
    @test_throws DimensionMismatch dearevenuegda([1; 2 ; 3], [4 ; 5], [1; 1; 1], :ERG) #  Different number of observations
    @test_throws DimensionMismatch dearevenuegda([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], :ERG) # Different number of observation in prices
    @test_throws DimensionMismatch dearevenuegda([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], :ERG) # Different number of output prices and outputs
    @test_throws ArgumentError dearevenuegda([1; 2; 3], [4; 5; 6], [4; 5; 6], :Error) # Invalid associated efficiency measure
    @test_throws ArgumentError dearevenuegda([1; 2; 3], [4; 5; 6], [1; 2; 3], :ERG, rts = :Error) # Invalid returns to scale

end


text:
@testset "Filter time series" begin
	sig = ones(15);
	imp = ones(5)./5;
	filtsig = mmconv(sig,imp);
	@test sum(filter(!isnan,filtsig)) ≈ 11.
	for i in [1,2,14,15]
		@test isnan(filtsig[i])
	end
	@test length(filtsig) == 15

	invec = collect(1.:1:17.);
	invec[[6,10,11,12,14,17]] .= NaN;
	start1,stop1 = findblocks(invec);
	@test start1 == [1,7,13,15];
	@test stop1 == [5,9,13,16];
	# Add NaN to the first position
	invec[1] = NaN;
	invec = vcat(invec,[18]);
	start2,stop2 = findblocks(invec);
	@test start2 == [2,7,13,15,18];
	@test stop2 == [5,9,13,16,18];
	# Test Vector input without NaNs
	start3,stop3 = findblocks([1.,2.,3.,4.,5.]);
	@test start3 == [1];
	@test stop3 == [5]

	invec = collect(1.:1:17.);
	invec[[6,10,11,12,14,17]] .= NaN;
	start1,stop1 = findnanblocks(invec);
	@test start1 == [6,10,14,17];
	@test stop1 == [6,12,14,17];
	# Add NaN to the first & remove from last position
	invec[1] = NaN;
	invec = vcat(invec,[18]);
	start2,stop2 = findnanblocks(invec);
	@test start2 == [1,6,10,14,17];
	@test stop2 == [1,6,12,14,17];
	# Test Vector input without NaNs
	start3,stop3 = findnanblocks([1.,2.,3.,4.,5.]);
	@test isempty(start3)
	@test isempty(stop3)

	sig = vcat(ones(10),NaN,ones(12) .+ 1);
	out = filtblocks(sig,ones(3)./3);
	for i in [1,10,11,12,length(sig)]
		@test isnan(out[i])
	end
	@test sum(filter(!isnan,out)) ≈ 8+10*2.0
end

@testset "De-mean & -trend" begin
	o = [-1.,NaN,0.,1.];
	s = o .+ 1.234
	sf = demean(s)
	for (i,v) in enumerate(sf)
		if i != 2
			@test v ≈ o[i]
		else
			@test isnan(v)
		end
	end

	x = collect(1.:1:10);
	y = ones(length(x));
	yc0 = detrend(x,y,deg=0);
	@test isapprox(sum(yc0),0.0)

	yc1 = detrend(x,x .* 0.5,deg=1);
	@test isapprox(sum(yc1),0.0,atol=1e-10)

	o = [-1.,NaN,0.,1.];
	sf = defirst(o)
	for (i,v) in enumerate(sf)
		if i != 2
			@test v ≈ o[i]-o[1]
		else
			@test isnan(v)
		end
	end
end