text:
"""
	fillnans(datavec,maxgap)

Find and (try to) remove all NaNs in the input data vector (`datavec`) via
linear interpolation

**Input:**
`datavec`: vector to be corrected for NaNs. It is assumed that this
	vector corresponds to equally sampled time vector! Use the `time2regular`
	function to ensure regular sampling, or check the sampling using `isregular`
	function.
`maxgap`: parameter controlling the maximum gap to be filled. Corresponds to
	index range.

**Output:**
corrected vector

**Example**
```
# NaN will be replaced by interpolated value
data = DataFrame(Temp=[10.,11.,12.,14.,NaN, 14.,],
	Grav=[1.,2.,NaN,NaN,NaN,6.],
	datetime=collect(DateTime(2010,1,1,1):Dates.Hour(1):DateTime(2010,1,1,6)));
out = fillnans(data[:Temp],2);

# To find indices corresponding to corrected values just search for difference
# `fillnans` was called:
corrindex = findall(isnan.(data[:Temp]) .& !isnan.(out)); # will return [5]

# NaN will NOT be replaced as the missing window is too long (3>2)
data[:Grav] = fillnans(data[:Grav],2);
```
"""
function fillnans(datavec::Vector{Float64},maxgap::Int)
   nanlines = .!isnan.(datavec);
   dataout = copy(datavec);
   for i in 2:length(datavec)-1
	   if nanlines[i]==false # only for NaNs
		   # Get maximum possible data range (the maxgap length will be checked inside lininterp fce)
		   starti = i <= maxgap ? 1 : i - maxgap
		   stopi  = i >= length(datavec)-maxgap ? length(datavec) : i + maxgap
		   # Interpolate only if valid data exist at both sides of the NaN
		   if any(nanlines[starti:i]) && any(nanlines[i:stopi])
			   dataout[i] = lininterp(datavec[starti:stopi],i-starti+1,maxgap);
		   end
	   end
   end
   return dataout
end

"""
	replacenans!(datain,replaceby)

Replace all NaNs by given value

**Input:**
`datain`: DataFrame to be corrected (all Float64 columns)
`replaceby`: replace NaNs by this value

**Output:**
corrected vector

**Example**
```
# NaN will be replaced by 0.0
datain = DataFrame(Temp=[10.,11.,12.,14.,NaN, 14.,],
	Grav=[1.,2.,NaN,NaN,NaN,6.],
	datetime=collect(DateTime(2010,1,1,1):Dates.Hour(1):DateTime(2010,1,1,6)));
out = replacenans!(datain,0.0);
```
"""
function replacenans!(datain::DataFrame,replaceby::Float64)
	for i in names(datain)
		if eltype(datain[i]) == Float64
			for j in 1:length(datain[i])
				if isnan(datain[i][j])
					datain[i][j] = replaceby;
				end
			end
		end
	end
end

"""
Function to convert NA to NaNs (if (el)type=Float64)
"""
function na2nan!(datain::DataFrame)
	for i in names(datain)
		if eltype(datain[i]) == Union{Float64, Missings.Missing}
			datain[i] = collect(Missings.replace(datain[i],NaN))
		end
	end
end

"""
	missing2nan(datain)

Convert vector or DataFrame of type missing to Float64 or DateTime (missing replaced by NaNs)

**Input:**
`datain`: Vector/DataFrame to be converted to Float64 or DateTime Vector/DataFrame

**Output:**
corrected Vector/DataFrame

**Example**
```
df = DataFrame(datetime=[DateTime(2000),DateTime(2001)],
			datetime2 = [DateTime(2002),missing],
			grav=[missing,3.],
			pres=[1,missing],
			temp=[missing,missing]);
df[:datetime3] = convert(Array{Union{Missing, DateTime},1},df[:datetime])

@test missing2nan(df[:datetime]) == [DateTime(2000),DateTime(2001)]
@test missing2nan(df[:datetime3]) == [DateTime(2000),DateTime(2001)]
@test missing2nan(df[:datetime2])[1] == DateTime(2002)
@test ismissing(missing2nan(df[:datetime2])[2])
@test missing2nan(df[:grav])[2] == 3.
@test missing2nan(df[:pres])[1] == 1.
@test isnan(missing2nan(df[:grav])[1])
@test isnan(missing2nan(df[:pres])[2])
@test isnan(missing2nan(df[:temp])[1])
@test isnan(missing2nan(df[:temp])[1])

# input should stay unchanged
@test ismissing(df[:grav][1])
@test df[:pres][1] == 1
@test eltype(typeof(df[:temp])) == Missing

# test for dataframe input
dfout = missing2nan(df)

@test missing2nan(dfout[:datetime]) == [DateTime(2000),DateTime(2001)]
@test missing2nan(dfout[:datetime3]) == [DateTime(2000),DateTime(2001)]
@test missing2nan(dfout[:datetime2])[1] == DateTime(2002)
@test ismissing(missing2nan(dfout[:datetime2])[2])
@test missing2nan(dfout[:grav])[2] == 3.
@test missing2nan(dfout[:pres])[1] == 1.
@test isnan(missing2nan(dfout[:grav])[1])
@test isnan(missing2nan(dfout[:pres])[2])
@test isnan(missing2nan(dfout[:temp])[1])
@test isnan(missing2nan(dfout[:temp])[1])

# input should stay unchanged
@test ismissing(df[:grav][1])
@test df[:pres][1] == 1
@test eltype(typeof(df[:temp])) == Missing
```
"""
function missing2nan(datain::Vector)::Vector
	outtype = eltype(typeof(datain));
	if outtype == DateTime || outtype == Int ||
		outtype == Float64 || outtype == String  || outtype == Symbol
		return datain
	elseif outtype == Missing
		return zeros(Float64,length(datain)).+NaN
	elseif !any(ismissing.(datain))
		return convert(Vector{outtype.b},datain)
	elseif outtype.b == DateTime
		return datain # cannot convert DateTime
	else
		out = zeros(Float64,length(datain)).+NaN;
		for i in 1:length(datain)
			if !ismissing(datain[i])
				out[i] = datain[i];
			end
		end
		return out
	end
end
function missing2nan(datain::DataFrame)::DataFrame
	dataout = copy(datain);
	for i in names(datain)
		dataout[i] = missing2nan(datain[i]);
	end
	return dataout
end

"""
Auxiliary function for linear interpolation
"""
function lininterp(datain::Vector{Float64},ind::Int,maxgap::Int)
	xl,xr,yl,yr = prepwindow(datain,ind);
	if xr - xl <= maxgap+1 # make sure the window is not too long
		return (yr-yl)/(xr-xl)*(ind-xl) + yl; # (slope)*distance + offset
	else
		return NaN;
	end
end
"""
Auxiliary function to prepare data for interpolation (remove NaNs and convert
	input to x,y coordinates)
"""
function prepwindow(datain::Vector{Float64},ind::Int)
	li = datain[1:ind]; # values left of NaN
	ri = datain[ind:end]; # values right of NaN
	rl = .!isnan.(li); # find all valid values on left
	rr = .!isnan.(ri); # ------------------------ right
	xl = 1:ind |> x -> x[rl][end]; # get only left x coodinate closest to NaN
	xr = ind:length(datain) |> x -> x[rr][1];#right--------------------
	yl = li[rl][end]; # y coordinate on left
	yr = ri[rr][1];	  # ----------------right
	return xl,xr,yl,yr
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SpatialPartition(object, subsets, [metadata])

A partition of a spatial `object` into `subsets`.
Optionally, save `metadata` as a dictionary.
"""
struct SpatialPartition{O}
  object::O
  subsets::Vector{Vector{Int}}
  metadata::Dict
end

SpatialPartition(object, subsets, metadata=Dict()) =
  SpatialPartition{typeof(object)}(object, subsets, metadata)

"""
    subsets(partition)

Return the subsets of indices in spatial object
that make up the `partition`.
"""
subsets(partition::SpatialPartition) = partition.subsets

"""
    metadata(partition)

Return the metadata dictionary saved in the partition.
"""
metadata(partition::SpatialPartition) = partition.metadata

# --------------
# INDEXABLE API
# --------------

Base.iterate(partition::SpatialPartition, state=1) =
  state > length(partition) ? nothing : (partition[state], state + 1)

Base.length(partition::SpatialPartition) = length(partition.subsets)

Base.getindex(partition::SpatialPartition, ind::Int) =
  view(partition.object, partition.subsets[ind])

Base.getindex(partition::SpatialPartition, inds::AbstractVector{Int}) =
  [getindex(partition, ind) for ind in inds]

# ------------
# IO methods
# ------------
function Base.show(io::IO, partition::SpatialPartition)
  nsubsets = length(partition.subsets)
  print(io, "$nsubsets SpatialPartition")
end

function Base.show(io::IO, ::MIME"text/plain", partition::SpatialPartition)
  subs = partition.subsets
  meta = partition.metadata
  lines = ["  └─$(length(sub))" for sub in subs]
  lines = length(lines) > 11 ? [lines[1:5]; ["  ⋮"]; lines[end-4:end]] : lines
  println(io, partition)
  println(io, "  N° points")
  print(io, join(lines, "\n"))
  if !isempty(meta)
    print(io, "\n  metadata: ", join(keys(meta), ", "))
  end
end

"""
    PartitionMethod

A method for partitioning spatial objects.
"""
abstract type PartitionMethod end

"""
    partition(object, method)

Partition `object` with partition `method`.
"""
function partition end

"""
    PredicatePartitionMethod

A method for partitioning spatial objects with predicate functions.
"""
abstract type PredicatePartitionMethod <: PartitionMethod end

function partition(object, method::PredicatePartitionMethod)
  subsets = Vector{Vector{Int}}()
  for i in randperm(nelms(object))
    inserted = false
    for subset in subsets
      j = subset[1]
      if method(i, j)
        push!(subset, i)
        inserted = true
        break
      end
    end

    if !inserted
      push!(subsets, [i])
    end
  end

  SpatialPartition(object, subsets)
end

"""
    SPredicatePartitionMethod

A method for partitioning spatial objects with spatial predicate functions.
"""
abstract type SPredicatePartitionMethod <: PartitionMethod end

function partition(object, method::SPredicatePartitionMethod)
  N = ncoords(object)
  T = coordtype(object)

  # pre-allocate memory for coordinates
  x = MVector{N,T}(undef)
  y = MVector{N,T}(undef)

  subsets = Vector{Vector{Int}}()
  for i in randperm(nelms(object))
    coordinates!(x, object, i)

    inserted = false
    for subset in subsets
      coordinates!(y, object, subset[1])

      if method(x, y)
        push!(subset, i)
        inserted = true
        break
      end
    end

    if !inserted
      push!(subsets, [i])
    end
  end

  SpatialPartition(object, subsets)
end

#------------------
# IMPLEMENTATIONS
#------------------
include("partitioning/random.jl")
include("partitioning/fraction.jl")
include("partitioning/slic.jl")
include("partitioning/block.jl")
include("partitioning/bisect_point.jl")
include("partitioning/bisect_fraction.jl")
include("partitioning/ball.jl")
include("partitioning/plane.jl")
include("partitioning/direction.jl")
include("partitioning/predicate.jl")
include("partitioning/spatial_predicate.jl")
include("partitioning/variable.jl")
include("partitioning/product.jl")
include("partitioning/hierarchical.jl")

ext: jl


text: print("Hello World\n")

max_stars_repo_name: MarshMue/math-tweets

max_issues_repo_name: MarshMue/math-tweets

max_forks_repo_name: MarshMue/math-tweets


text:
export StochEnKF, SeqStochEnKF

"""
$(TYPEDEF)

A structure for the stochastic ensemble Kalman filter (sEnKF)

References:

Evensen, G. (1994). Sequential data assimilation with a nonlinear quasi‐geostrophic model using Monte Carlo methods to forecast error statistics. Journal of Geophysical Research: Oceans, 99(C5), 10143-10162.
## Fields
$(TYPEDFIELDS)
"""

struct StochEnKF<:SeqFilter
    "Filter function"
    G::Function

    "Standard deviations of the measurement noise distribution"
    ϵy::AdditiveInflation

    "Time step dynamic"
    Δtdyn::Float64

    "Time step observation"
    Δtobs::Float64

	"Boolean: is the covariance matrix localized"
	islocal::Bool

    "Boolean: is state vector filtered"
    isfiltered::Bool
end

function StochEnKF(G::Function, ϵy::AdditiveInflation,
    Δtdyn, Δtobs; islocal = false, isfiltered = false)
    @assert norm(mod(Δtobs, Δtdyn))<1e-12 "Δtobs should be an integer multiple of Δtdyn"

    return StochEnKF(G, ϵy, Δtdyn, Δtobs, islocal, isfiltered)
end

# If no filtering function is provided, use the identity in the constructor.
function StochEnKF(ϵy::AdditiveInflation,
    Δtdyn, Δtobs; islocal = false)
    @assert norm(mod(Δtobs, Δtdyn))<1e-12 "Δtobs should be an integer multiple of Δtdyn"

    return StochEnKF(x -> x, ϵy, Δtdyn, Δtobs, islocal, false)
end

function Base.show(io::IO, enkf::StochEnKF)
	println(io,"Stochastic EnKF  with filtered = $(enkf.isfiltered)")
end

function (enkf::StochEnKF)(X, ystar::Array{Float64,1}, t::Float64; laplace::Bool=false)

    Ny = size(ystar,1)
    Nx = size(X,1)-Ny
    Ne = size(X, 2)

    meas  = viewmeas(X,Ny,Nx)
    state = viewstate(X,Ny,Nx)

    x̄ = copy(mean(X[Ny+1:Ny+Nx,:],dims=2)[:,1])
    Xf = copy(X[Ny+1:Ny+Nx,:])
    Xf .-= x̄
	rmul!(Xf, 1/sqrt(Ne-1))

	# Need the covariance to perform the localisation
	u = zeros(Ny, Ne)
	if laplace == false
		u .= enkf.ϵy.σ*randn(Ny, Ne) .+ enkf.ϵy.m
	else
		u .= sqrt(2.0)*enkf.ϵy.σ*rand(Laplace(),(Ny, Ne)) .+ enkf.ϵy.m
	end
	# ū = mean(u,dims=2)[:,1]
	Yf = copy(meas .- mean(meas; dims=2)[:,1])

	rmul!(Yf,1/sqrt(Ne-1))

	"Analysis step with representers, Evensen, Leeuwen et al. 1998"

	b = (Yf*Yf' + enkf.ϵy.Σ) \ (ystar .+ (u - meas))

	Bᵀb = (Xf*Yf')*b

	state .+= Bᵀb

	return X
end

function (enkf::StochEnKF)(X, ystar::Array{Float64,1}, ȳf::Array{Float64,1}; laplace::Bool=false)
    Ny = size(ystar,1)
    Nx = size(X,1)-Ny
    Ne = size(X, 2)

    meas  = viewmeas(X,Ny,Nx)
    state = viewstate(X,Ny,Nx)

    x̄ = copy(mean(X[Ny+1:Ny+Nx,:],dims=2)[:,1])

    Xf = copy(X[Ny+1:Ny+Nx,:])
    Xf .-= x̄

	# Need the covariance to perform the localisation
	u = zeros(Ny, Ne)
	if laplace == false
		u .= enkf.ϵy.σ*randn(Ny, Ne) .+ enkf.ϵy.m
	else
		u .= sqrt(2.0)*enkf.ϵy.σ*rand(Laplace(),(Ny, Ne)) .+ enkf.ϵy.m
	end
	ū = mean(u,dims=2)[:,1]

	Yf = copy(meas) .- ȳf

	"Analysis step with representers, Evensen, Leeuwen et al. 1998"
		b = (Yf*transpose(Yf) + (Ne-1)*enkf.ϵy.Σ)\(ystar .+ (u - meas))
		Bᵀb = (Xf*Yf')*b
		# @show norm(Bᵀb)
		state .+= Bᵀb

	return X
end

##### Sequential and localized ensemble Kalman filter

struct SeqStochEnKF<:SeqFilter
    "Filter function"
    G::Function

	"Local observation function"
	h::Function

	"Multiplicative inflation"
	β::MultiplicativeInflation

    "Standard deviations of the measurement noise distribution"
    ϵy::AdditiveInflation

	"Observation dimension"
	Ny::Int64

	"State dimension"
	Nx::Int64

    "Time step dynamic"
    Δtdyn::Float64

    "Time step observation"
    Δtobs::Float64

	"Index of measurement"
	idx::Array{Int64,2}
	# idx contains the dictionnary of the mapping
	# First line contains the range of integer 1:Ny
	# Second line contains the associated indice of each measurement

	# Define cache
	"Cache for the inference"
	cache::Array{Float64,2}

	"Boolean: is state vector filtered"
    isfiltered::Bool

	"Boolean: is assimilation localized"
	islocalized::Bool
end

function SeqStochEnKF(G::Function, h::Function, β::Float64, ϵy::AdditiveInflation,
	 Ny, Nx, Ne, Δtdyn, Δtobs, idx::Array{Int64,2}; isfiltered = false, islocalized = true)
    @assert norm(mod(Δtobs, Δtdyn))<1e-12 "Δtobs should be an integer multiple of Δtdyn"

    return SeqStochEnKF(G, h, MultiplicativeInflation(β), ϵy, Ny, Nx, Δtdyn, Δtobs, idx, zeros(Nx+1, Ne), isfiltered, islocalized)
end

function Base.show(io::IO, enkf::SeqStochEnKF)
	println(io,"Sequential Stochastic EnKF  with filtered = $(enkf.isfiltered)")
end

# """
#     Define action of SeqStocEnKF on EnsembleStateMeas
# """
# Baptista, Spantini, Marzouk 2019

# Version with localization
function (enkf::SeqStochEnKF)(X, ystar, t, Loc::Localization)
	h = enkf.h

	Nx = enkf.Nx
	Ny = enkf.Ny
	Na = Nx+1
	NxX, Ne = size(X)

	cache = enkf.cache
	fill!(cache, 0.0)

	# Localize covariance
	locXY = Locgaspari((Nx, Ny), Loc.L, Loc.Gxy)

	# Sequential assimilation of the observations
	@inbounds for i=1:Ny
		idx1, idx2 = enkf.idx[:,i]
		ylocal = ystar[idx1]

		# Inflate state
		cache[2:end,:] .= X[Ny+1:Ny+Nx,:]
		Aβ = enkf.β
		Aβ(cache)

		# Generate samples from local likelihood with inflated ensemble
		@inbounds for i=1:Ne
			col = view(X, Ny+1:Ny+Nx, i)
			cache[1,i] = h(col, idx2, t) + enkf.ϵy.m[idx1] + dot(enkf.ϵy.σ[idx1,:], randn(Ny))
		end

		XYf = copy(cache) .- mean(cache; dims = 2)[:,1]
		rmul!(XYf, 1/sqrt(Ne-1))

		Yf = view(XYf,1,:)
		Xf = view(XYf,2:Nx+1, :)

		#Construct the observation with the un-inflated measurement
		@inbounds for i=1:Ne
			col = view(X, Ny+1:Ny+Nx, i)
			cache[1,i] = h(col, idx2, t) + enkf.ϵy.m[idx1] + dot(enkf.ϵy.σ[idx1,:], randn(Ny))
		end

		# Since Yf is a vector (Yf Yf') is reduced to a scalar
		"Analysis step with representers, Evensen, Leeuwen et al. 1998"
		K = view(locXY,:,idx1).* (Xf*Yf)
		# @show K./dot(Yf,Yf)
		# K ./= dot(Yf, Yf)
		# Bᵀb = K*(ylocal .- view(enkf.ensa.S,1,:))
		#In-place analysis step with gemm!
		BLAS.gemm!('N', 'T', 1/dot(Yf,Yf), K, ylocal .- view(cache,1,:), 1.0, view(X, Ny+1:Ny+Nx, :))
	end
	return X
end

# """
#     Define action of StocEnKF on EnsembleStateMeas
# """
# Bocquet Data assimilation p.160 Chapter 6 Stochastic EnKF
# function (enkf::StochEnKF)(ens::EnsembleStateMeas{Nx, Ny, Ne}, ystar::Array{Float64,1}, t::Float64; laplace::Bool=false) where {Nx, Ny, Ne}
#
# 	Xf = deviation(ens.state).S
# 	rmul!(Xf, 1/sqrt(Ne-1))
#
# 	# Need the covariance to perform the localisation
# 	u = zeros(Ny, Ne)
# 	if laplace == false
# 	u .= enkf.ϵy.σ*randn(Ny, Ne) .+ enkf.ϵy.m
# 	else
# 	u .= sqrt(2.0)*enkf.ϵy.σ*rand(Laplace(),(Ny, Ne)) .+ enkf.ϵy.m
# 	end
# 	# ū = mean(u,dims=2)[:,1]
# 	Yf = deviation(ens.meas).S
# 	rmul!(Yf,1/sqrt(Ne-1))
#
# 	"Analysis step with representers, Evensen, Leeuwen et al. 1998"
# 	b = (Yf*Yf' + enkf.ϵy.Σ) \ (ystar .+ (u - ens.meas.S))
#
# 	Bᵀb = (Xf*Yf')*b
#
# 	ens.state.S .+= Bᵀb
#
# 	return ens
# end
#
# # Bocquet Data assimilation p.160 Chapter 6 Stochastic EnKF
# function (enkf::StochEnKF)(ens::EnsembleStateMeas{Nx, Ny, Ne}, ystar::Array{Float64,1}, t::Float64; laplace::Bool=false) where {Nx, Ny, Ne}
#
# 	Xf = deviation(ens.state).S
# 	rmul!(Xf, 1/sqrt(Ne-1))
#
# 	# Need the covariance to perform the localisation
# 	u = zeros(Ny, Ne)
# 	if laplace == false
# 	u .= enkf.ϵy.σ*randn(Ny, Ne) .+ enkf.ϵy.m
# 	else
# 	u .= sqrt(2.0)*enkf.ϵy.σ*rand(Laplace(),(Ny, Ne)) .+ enkf.ϵy.m
# 	end
# 	# ū = mean(u,dims=2)[:,1]
# 	Yf = deviation(ens.meas).S
# 	rmul!(Yf,1/sqrt(Ne-1))
#
# 	"Analysis step with representers, Evensen, Leeuwen et al. 1998"
# 	b = (Yf*Yf' + enkf.ϵy.Σ) \ (ystar .+ (u - ens.meas.S))
#
# 	Bᵀb = (Xf*Yf')*b
#
# 	ens.state.S .+= Bᵀb
#
# 	return ens
# end

# # Version with localization
# function (enkf::StochEnKF)(ens::EnsembleStateMeas{Nx, Ny, Ne}, ystar::Array{Float64,1}, t::Float64, Loc::Localization) where {Nx, Ny, Ne}
#
# 	Xf = deviation(ens.state).S
# 	rmul!(Xf, 1/sqrt(Ne-1))
#
# 	# Localize covariance
# 	locXY = Locgaspari((Nx, Ny), Loc.L, Loc.Gxy)
# 	locYY = Locgaspari((Ny, Ny), Loc.L, Loc.Gyy)
#
# 	# Need the covariance ot peform the localisation
# 	u = zeros(Ny, Ne)
# 	u .= enkf.ϵy.σ*randn(Ny, Ne) .+ enkf.ϵy.m
# 	ū = mean(u,dims=2)[:,1]
#
# 	Yf = deviation(ens.meas).S
# 	rmul!(Yf,1/sqrt(Ne-1))
#
# 	"Analysis step with representers, Evensen, Leeuwen et al. 1998"
# 	b = (locYY .* (Yf*Yf') + enkf.ϵy.Σ) \ ((ū + ystar) .- (u + ens.meas.S))
# 	Bᵀb = (locXY .* (Xf*Yf'))*b
#
#
# 	@inbounds for i = 1:Ne
# 		Bᵀbi = view(Bᵀb,:,i)
# 		ens.state.S[:,i] .+= Bᵀbi
# 	end
# 	return ens
# end


text:
import Rates
using Test

@testset "Continuous Compounding" begin
    r = 0.1
    t = 1.5

    @test Rates.erf(Rates.Continuous, r, t) ≈ exp(r*t)
    @test Rates.discountfactor(Rates.Continuous, r, t) ≈ 1.0/exp(r*t)
    @test Rates.erf_to_rate(Rates.Continuous, exp(r*t), t) ≈ r
    @test Rates.discountfactor_to_rate(Rates.Continuous, 1.0/exp(r*t), t) ≈ r
    @test_throws AssertionError Rates.discountfactor(Rates.Continuous, r, -0.1)
    @test_throws AssertionError Rates.discountfactor_to_rate(Rates.Continuous, 1.0, -0.1)
end

@testset "Simple Compounding" begin
    r = 0.1
    t = 1.5

    @test Rates.erf(Rates.Simple, r, t) ≈ 1+r*t
    @test Rates.discountfactor(Rates.Simple, r, t) ≈ 1.0/(1+r*t)
    @test Rates.erf_to_rate(Rates.Simple, 1+r*t, t) ≈ r
    @test Rates.discountfactor_to_rate(Rates.Simple, 1.0/(1+r*t), t) ≈ r
    @test_throws AssertionError Rates.discountfactor(Rates.Simple, r, -0.1)
    @test_throws AssertionError Rates.discountfactor_to_rate(Rates.Simple, 1.0, -0.1)
end

@testset "Exponential Compounding" begin
    r = 0.1
    t = 1.5

    @test Rates.erf(Rates.Exponential, r, t) ≈ (1+r)^t
    @test Rates.discountfactor(Rates.Exponential, r, t) ≈ 1.0/((1+r)^t)
    @test Rates.erf_to_rate(Rates.Exponential, (1+r)^t, t) ≈ r
    @test Rates.discountfactor_to_rate(Rates.Exponential, 1.0/((1+r)^t), t) ≈ r
    @test_throws AssertionError Rates.discountfactor(Rates.Exponential, r, -0.1)
    @test_throws AssertionError Rates.discountfactor_to_rate(Rates.Exponential, 1.0, -0.1)
end


text:
@testset "Faces" begin
    @testset "Tetrahedron" begin
        for T in Floats
            P₁ = Point{3,T}(0, 0, 0)
            P₂ = Point{3,T}(1, 0, 0)
            P₃ = Point{3,T}(0, 1, 0)
            P₄ = Point{3,T}(0, 0, 1)
            tet = Tetrahedron(P₁, P₂, P₃, P₄)
            tet_faces = faces(tet)
            @test length(tet_faces) == 4
            @test tet_faces[1] == Triangle(P₁, P₂, P₃)
            @test tet_faces[2] == Triangle(P₁, P₂, P₄)
            @test tet_faces[3] == Triangle(P₂, P₃, P₄)
            @test tet_faces[4] == Triangle(P₃, P₁, P₄)
        end
    end

    @testset "Hexahedron" begin
        for T in Floats
            P₁ = Point{3,T}(0, 0, 0)
            P₂ = Point{3,T}(1, 0, 0)
            P₃ = Point{3,T}(1, 1, 0)
            P₄ = Point{3,T}(0, 1, 0)
            P₅ = Point{3,T}(0, 0, 1)
            P₆ = Point{3,T}(1, 0, 1)
            P₇ = Point{3,T}(1, 1, 1)
            P₈ = Point{3,T}(0, 1, 1)
            hex = Hexahedron(P₁, P₂, P₃, P₄, P₅, P₆, P₇, P₈)
            v = (P₁, P₂, P₃, P₄, P₅, P₆, P₇, P₈)
            hex_faces = faces(hex)
            @test length(hex_faces) == 6
            @test hex_faces[1] == Quadrilateral(v[1], v[2], v[3], v[4])
            @test hex_faces[2] == Quadrilateral(v[5], v[6], v[7], v[8])
            @test hex_faces[3] == Quadrilateral(v[1], v[2], v[6], v[5])
            @test hex_faces[4] == Quadrilateral(v[2], v[3], v[7], v[6])
            @test hex_faces[5] == Quadrilateral(v[3], v[4], v[8], v[7])
            @test hex_faces[6] == Quadrilateral(v[4], v[1], v[5], v[8])
        end
    end

    @testset "QuadraticTetrahedron" begin
        for T in Floats
            P₁  = Point{3,T}(0, 0, 0)
            P₂  = Point{3,T}(1, 0, 0)
            P₃  = Point{3,T}(0, 1, 0)
            P₄  = Point{3,T}(0, 0, 1)
            P₅  = Point{3,T}(1//2,    0,    0)
            P₆  = Point{3,T}(1//2, 1//2,    0)
            P₇  = Point{3,T}(0,    1//2,    0)
            P₈  = Point{3,T}(0,       0, 1//2)
            P₉  = Point{3,T}(1//2,    0, 1//2)
            P₁₀ = Point{3,T}(0,    1//2, 1//2)
            tet = QuadraticTetrahedron(P₁, P₂, P₃, P₄, P₅, P₆, P₇, P₈, P₉, P₁₀)
            v = (P₁, P₂, P₃, P₄, P₅, P₆, P₇, P₈, P₉, P₁₀)
            tet_faces = faces(tet)
            @test length(tet_faces) == 4
            @test tet_faces[1] == QuadraticTriangle(v[1], v[2], v[3], v[5],  v[6],  v[7])
            @test tet_faces[2] == QuadraticTriangle(v[1], v[2], v[4], v[5],  v[9],  v[8])
            @test tet_faces[3] == QuadraticTriangle(v[2], v[3], v[4], v[6], v[10],  v[9])
            @test tet_faces[4] == QuadraticTriangle(v[3], v[1], v[4], v[7],  v[8], v[10])
        end
    end

    @testset "QuadraticHexahedron" begin
        for T in Floats
            P₁  = Point{3,T}(   0,    0,    0)
            P₂  = Point{3,T}(   1,    0,    0)
            P₃  = Point{3,T}(   1,    1,    0)
            P₄  = Point{3,T}(   0,    1,    0)
            P₅  = Point{3,T}(   0,    0,    1)
            P₆  = Point{3,T}(   1,    0,    1)
            P₇  = Point{3,T}(   1,    1,    1)
            P₈  = Point{3,T}(   0,    1,    1)
            P₉  = Point{3,T}(1//2,    0,    0)
            P₁₀ = Point{3,T}(   1, 1//2,    0)
            P₁₁ = Point{3,T}(1//2,    1,    0)
            P₁₂ = Point{3,T}(   0, 1//2,    0)
            P₁₃ = Point{3,T}(1//2,    0,    1)
            P₁₄ = Point{3,T}(   1, 1//2,    1)
            P₁₅ = Point{3,T}(1//2,    1,    1)
            P₁₆ = Point{3,T}(   0, 1//2,    1)
            P₁₇ = Point{3,T}(   0,    0, 1//2)
            P₁₈ = Point{3,T}(   1,    0, 1//2)
            P₁₉ = Point{3,T}(   1,    1, 1//2)
            P₂₀ = Point{3,T}(   0,    1, 1//2)
            hex = QuadraticHexahedron(P₁,  P₂,  P₃,  P₄,  P₅,  P₆,  P₇,
                                      P₈,  P₉,  P₁₀, P₁₁, P₁₂, P₁₃, P₁₄,
                                      P₁₅, P₁₆, P₁₇, P₁₈, P₁₉, P₂₀)
            v = (P₁,  P₂,  P₃,  P₄,  P₅,  P₆,  P₇,
                 P₈,  P₉,  P₁₀, P₁₁, P₁₂, P₁₃, P₁₄,
                 P₁₅, P₁₆, P₁₇, P₁₈, P₁₉, P₂₀)
            hex_faces = faces(hex)
            @test length(hex_faces) == 6
            @test hex_faces[1] == QuadraticQuadrilateral(v[1], v[2], v[3], v[4], v[ 9], v[10], v[11], v[12])
            @test hex_faces[2] == QuadraticQuadrilateral(v[5], v[6], v[7], v[8], v[13], v[14], v[15], v[16])
            @test hex_faces[3] == QuadraticQuadrilateral(v[1], v[2], v[6], v[5], v[ 9], v[18], v[13], v[17])
            @test hex_faces[4] == QuadraticQuadrilateral(v[2], v[3], v[7], v[6], v[10], v[19], v[14], v[18])
            @test hex_faces[5] == QuadraticQuadrilateral(v[3], v[4], v[8], v[7], v[11], v[20], v[15], v[19])
            @test hex_faces[6] == QuadraticQuadrilateral(v[4], v[1], v[5], v[8], v[12], v[17], v[16], v[20])
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
if !@isdefined TEST
    using Revise
end
using Mitosis
using Random, Test, LinearAlgebra, Statistics
using Soss
using StatsBase
DIAG = true
if DIAG
    using MCMCChains
end

N = 50_000 # samples

# Define some vectors and matrices
ξ0 = [1., 0.]
x = [1.2, 0.1]
P0 = Matrix(1.0I, 2, 2)

Φ = [0.8 0.5; -0.1 0.8]
β = [0.1, 0.2]
Q = [0.2 0.0; 0.0 1.0]

H = [1.0 0.0]
R = Matrix(1.0I, 1, 1)

# and a nonlinear vector function
f(x) = [0.8 0.5; -0.1 0.8]*[atan(x[1]), atan(x[2])] + [0.1, 0.2]

# We define the equivalent of the Soss model
(m = @model ξ0, H, f, P0, R, Q begin
    a ~ MvNormal(ξ0, P0) # priortransition
    b ~ MvNormal(H*a, R) # partial observation
    c ~ MvNormal(f(a), Q) # nonlineartransition
    d ~ MvNormal(H*c, R) # partial observation
    e ~ MvNormal(c, P0) # full observation
    # return b, d, e
end)

# Define some transition kernels.

# We use AffineMap, LinearMap and ConstantMap
# For example
@test AffineMap(Φ, β)(x) == Φ*x + β

# Prior
# a ~ N(ξ0, P0)
k_0a = Mitosis.kernel(Gaussian; μ=LinearMap(I(2)), Σ=ConstantMap(P0))

# Partial observation
# b ~ N(H*a, R)
k_ab = Mitosis.kernel(Gaussian; μ=LinearMap(H), Σ=ConstantMap(R))

# Nonlinear transition and linear approximation
# c ~ N(f(a), Q)
k_ac = Mitosis.kernel(Gaussian; μ=f, Σ=ConstantMap(Q))
k_ac_approx = Mitosis.kernel(Gaussian; μ=AffineMap(Φ, β), Σ=ConstantMap(Q))

# d ~ N(H*c, R)
k_cd = Mitosis.kernel(Gaussian; μ=LinearMap(H), Σ=ConstantMap(R))

# e ~ N(c, P0)
k_ce = Mitosis.kernel(Gaussian; μ=LinearMap(I(2)), Σ=ConstantMap(P0))

# And a deterministic copy kernel
cp2 = Mitosis.Copy{2}() # copy kernel
@test cp2(1.1) == (1.1, 1.1)

# Forward sample a Bayes net # rand(m(ξ0=ξ0))
a = rand(k_0a(ξ0))
b = rand(k_ab(a))
c = rand(k_ac(a)) # forward model with nonlinear transition
d = rand(k_cd(c))
e = rand(k_ce(c))

# We actually want to write down the
# forward sample with explicit copies
# so every state is only used as input of a single kernel
a = rand(k_0a(ξ0))
a1, a2 = cp2(a)
b = rand(k_ab(a1))
c = rand(k_ac(a2)) # forward model with nonlinear transition
c1, c2 = cp2(c)
d = rand(k_cd(c1))
e = rand(k_ce(c2))

# sample posterior with dynamicHMC via Soss for comparison
# monkey patch for xform
Soss.xform(d::MvNormal, _data) = Soss.as(Array, Soss.asℝ, size(d))
@time posterior = dynamicHMC(m(ξ0=ξ0, H=H, f=f, P0=P0, R=R, Q=Q), (b=b, d=d, e=e,), N)

# observations are b d e (all leaves)
mc2, pc2 = backward(BFFG(), k_ce, e; unfused=true) # child of copy
mc1, pc1 = backward(BFFG(), k_cd, d; unfused=true) # child of copy
mc, pc = backward(BFFG(), cp2, pc1, pc2) # copy
ma2, pa2 = backward(BFFG(), k_ac_approx, pc; unfused=true) # # child of copy
ma1, pa1 = backward(BFFG(), k_ab, b; unfused=true) # child of copy
ma, pa = backward(BFFG(), cp2, pa1, pa2) # copy
m0, p0 = backward(BFFG(), k_0a, pa) # not a child of a copy

samples = []
lls = []
@time for i in 1:N
    a = rand(forward(BFFG(), k_0a, m0, weighted(ξ0)))
    a1, a2 = rand(forward(BFFG(), cp2, ma, a))
    b = rand(forward(BFFG(), k_ab, ma1, a1))
    c = rand(forward(BFFG(), (k_ac, k_ac_approx), ma2, a2)) # forward model with nonlinear transition
    c1, c2 = rand(forward(BFFG(), cp2, mc, c))
    d = rand(forward(BFFG(), k_cd, mc1, c1))
    e = rand(forward(BFFG(), k_ce, mc2, c2))

    # Sample and weight (sum over the weights of the leafs.)
    push!(samples, (a=a.x,b=b.x,c=c.x,d=d.x,e=e.x))
    push!(lls, b.ll + d.ll + e.ll)

end

# Importance weights
w = exp.(lls)/mean(exp.(lls))

# Estimate Mitosis
â1 = mean(getfield.(samples, :a).*w)
v1 = mean(Mitosis.outer.(getfield.(samples, :a).-Ref(â1)).*w)

# Estimate Soss
â2 = mean(getfield.(posterior, :a))
v2 = cov(getfield.(posterior, :a))

if DIAG
    # Kish's Effective Sample Size
    ess1 = sum(w).^2/sum(w.^2)
    ess2 = MCMCChains.ess(Chains(reshape(first.(getfield.(posterior, :a)), :,1,1)))
end

[â1 â2]


text:
#=
This file defines functions that simply return all timeseries of velocity
and timing from the different sources we have. E.g. PG13, Uwe recordings,
Drum playalongs, drum tapping, etc.
Any extra options (like e.g. using only melody, etc.), are passed
as keyword arguments.

All functions always return (mtds, vels)
=#
include("drums_processing.jl")

function uwe_data(;subset = "all")
    uwedir = datadir("exp_raw", "uwe")
    mtds = []
    vels = []
    for file in readdir(uwedir)
        midi = readMIDIFile(joinpath(uwedir, file))
        piano = getnotes(midi, 2)
        d = estimate_delay(piano, 0:1//3:1)
        piano = translate(piano, -round(Int, d))

        asr, σasr = Jazz.average_swing_ratio(piano, "swung8s")

        x = asr/(asr + 1)
        grid = [0, x/2, x, 1]

        if subset == "all"
            qnotes = quantize(piano, grid)
            class = classify(piano, grid)
            # in the following do not use notes in the x/2 bin; outliers
            a = findall(!isequal(2), class)

            m = positions(piano)[a] .- positions(quantize(piano, grid))[a]
            v = velocities(piano)[a]
        elseif subset == "melody"
            # only highest pitch notes
            fvel(notes) = ((m, i) = findmax(pitches(notes)); notes[i].velocity)
            fpos(notes) = ((m, i) = findmax(pitches(notes)); notes[i].position)
            t, v = timeseries(piano, fvel, grid)

            tt, p = timeseries(piano, fpos, grid)
            nomiss = @. !ismissing.(v)

            m = @. round(Int, p[nomiss] .- tt[nomiss])
            v = @. round(Int, v[nomiss])
        end
        push!(vels, v); push!(mtds, m)
    end
    return mtds, vels
end

function playalongs_data(; subset = "all")
    pdir = datadir("exp_raw", "playalongs")
    mtds = []
    vels = []
    for file in readdir(pdir)
        grid = file[3] == 'j' ? (0:1//3:1) : (0:1//4:1)
        midi = readMIDIFile(joinpath(pdir, file))
        drums = Drums.getnotes_td50(midi.tracks[2], midi.tpq)
        drums = drums_cleanup(drums)

        d = estimate_delay_recursive(drums, grid, 5)
        drums = translate(drums, -d)

        if subset == "all"
            drums = drums
        elseif subset == "lead"
            drums = filterpitches(drums, [0x1a, 0x2e, 0x35, 0x33])
        elseif subset == "snare"
            drums = filterpitches(drums, [0x26, 0x28])
        elseif subset == "kick"
            drums = filterpitches(drums, [0x24])
        elseif subset == "quarters"
            class = classify(drums, grid)
            keep = [c ∈ (1, length(grid)) ? true : false for c in class]
            drums = drums[keep]
        elseif subset == "aligned"
            channeled = channel_direct(drums)
            aligned = align_notes(channeled)
            drums = combine(aligned)
        end
        drums = remove_peaking_velocities!(drums)
        qnotes = quantize(drums, grid)
        m = positions(drums) .- positions(qnotes)
        v = velocities(drums)
        length(v) == 0 && continue
        push!(vels, v); push!(mtds, m)
    end
    return mtds, vels
end

function tapping_data(;subset = "all")
    cordir = datadir("exp_raw", "tapping")
    mtds = []
    vels = []
    for file in readdir(cordir)
        midi = readMIDIFile(joinpath(cordir, file))
        hits = getnotes(midi)
        grid = 0:1
        d = round(Int, mean(relpos(hits, grid)))
        hits = translate(hits, -d)
        drum_q = quantize(hits, 0:1)
        mtd = positions(hits) .- positions(drum_q)
        vel = velocities(hits)
        push!(mtds,mtd)
        push!(vels,vel)
    end
    return mtds, vels
end

function find_first(notes)
  tmp, first_notes = Notes(), Notes()
  for (index,qn) in enumerate(quantize(notes,0:1)[1:end-1])
    push!(tmp, notes[index])
    if quantize(notes,0:1)[index+1].position != qn.position
      push!(first_notes, tmp[findmax(pitches(tmp))[2]])
      tmp = Notes()
    end
  end
  return first_notes
end

"""
    pgmusic_data(bool, interval) → mtd, vels

extracts the MTDs and velocities for all the pieces in the pgmusic datadir.

  input :
    - first : boolean argument that determines if we take the whole chords into account or just the top notes that corresponds to the melody line
    - interval : the interval into which we consider a note to be a quarter note. Since the recordings are complex, it is hqrd to find q single grid that describes everything for quantization.
      so every notes which is in [960 - lower_limit, upper_limit] ticks is considered a quarter note.
  returns :
    - mtd : a vector of vector containing the MTDs for each piece.
    - vel : a vector of vector containing the vels for each piece.
"""
function pgmusic_data(;subset = "all")
    if subset == "all"
        first = false
    elseif subset == "melody"
        first = true
    end
    boundary = [-100,125]

    pgdir = datadir("exp_raw","pgmusic")
    MTDs, vels = [], []
    for file in readdir(pgdir)
        notes = getnotes(readMIDIFile(joinpath(pgdir, file)),1)
        tmp, quarter_notes = Notes(), Notes()
        lower_limit, upper_limit = boundary[1], boundary[2]
        #Taking all the notes in the interval - limit: limit
        for n in notes
            if mod(n.position,960) >= 960 + lower_limit || mod(n.position,960) < upper_limit
                push!(tmp,n)
            end
        end
        if first
            quarter_notes = find_first(tmp)
        elseif first == false
            quarter_notes = tmp
        end
        # moduloing the positions, and substracting tpq when applicable to get mtdS
        # can't use '''quantize''' since we don't have a clear grid.
        pos = mod.(positions(quarter_notes),960)
        for i in 1:length(pos)
            if pos[i] >= 960 + lower_limit
            pos[i] -= 960
            end
        end
        push!(MTDs, pos)
        push!(vels, velocities(quarter_notes))
    end
    return MTDs, vels
end

"""
    find_top_notes(notes)
finds the top notes (gets rid of chords, keeps melody) on a grid [0:1//6:1].
"""
function find_top_notes(notes)
  tmp, first_notes = Notes(), Notes()
  for (index,qn) in enumerate(quantize(notes,collect(0:1//6:1))[1:end-1])
    push!(tmp, notes[index])
    if quantize(notes,collect(0:1//6:1))[index+1].position != qn.position
      push!(first_notes, tmp[findmax(pitches(tmp))[2]])
      tmp = Notes()
    end
  end
  return first_notes
end

function get_dataset(c::Dict)
    dataset = c[:dataset]
    subset = get(c, :subset, "all")
    get_dataset(dataset, subset)
end
function get_dataset(dataset, subset = "all")
    if dataset == "uwe"
        mtds, vels = uwe_data(;subset=subset)
    elseif dataset == "playalongs"
        mtds, vels = playalongs_data(;subset=subset)
    elseif dataset == "tapping"
        mtds, vels = tapping_data(;subset=subset)
    elseif dataset == "pgmusic"
        mtds, vels = pgmusic_data(;subset=subset)
    end
    return mtds, vels
end

get_dataset_names(c::Dict) = get_dataset_names(c[:dataset])
function get_dataset_names(dataset::String)
    dir = datadir("exp_raw", dataset)
    return [f[1:end-4] for f in readdir(dir)]
end

max_stars_repo_name: Datseris/VelocityMTD

max_issues_repo_name: Datseris/VelocityMTD

max_forks_repo_name: Datseris/VelocityMTD


text:
const Void = Cvoid
const NULL = C_NULL

#=
# These are not used anymore, I think.
const DOUBLE_CLASS = 1
const SINGLE_CLASS = 2
const INT64_CLASS  = 3
const UINT64_CLASS = 4
const INT32_CLASS  = 5
const UINT32_CLASS = 6
const INT16_CLASS  = 7
const UINT16_CLASS = 8
const INT8_CLASS   = 9
const UINT8_CLASS  = 10 =#

struct GMT_OPTION			# Structure for a single GMT command option
	option::UInt8			# 1-char command line -<option> (e.g. D in -D) identifying the option (* if file)
	arg::Cstring			# If not NULL, contains the argument for this option
	next::Ptr{GMT_OPTION}
	previous::Ptr{GMT_OPTION}
end

mutable struct GMT_PEN		# Structure to hold pen attributes
	width::Cdouble			# In points
	offset::Cdouble			# In points
	rgb::NTuple{4,Cdouble}	# RGB color of pen + Transparency 0-1 [0 = opaque] */
	style::NTuple{128,UInt8}
	# For line modifications
	mode::UInt32			# Line-type: PSL_LINEAR [0; default] or PSL_BEZIER [1]
	cptmode::UInt32		# How a cpt affects pens and fills: 0-none, 1=use CPT for line, 2 = use CPT for fill, 3 = both
	# end::NTuple{2,GMT_LINE_END}
	end_::NTuple{2,Ptr{Cvoid}}		# This is a dangereous thing. If accessed, will crash Julia
	GMT_PEN(width, offset, rgb, style, mode, cptmode, end_) = new(width, offset, rgb, style, mode, cptmode, end_)
	GMT_PEN() = new(0.0, 0.0, (0.0, 0.0, 0.0, 0.0), map(UInt8, (repeat('\0', 128)...,)), 0, 0, (pointer([0]), pointer([0])))
end

mutable struct GMT_GRID_HEADER
	n_columns::UInt32
	n_rows::UInt32
	registration::UInt32
	wesn::NTuple{4,Cdouble}
	z_min::Cdouble
	z_max::Cdouble
	inc::NTuple{2,Cdouble}
	z_scale_factor::Cdouble
	z_add_offset::Cdouble
	x_unit::NTuple{80,UInt8}
	y_unit::NTuple{80,UInt8}
	z_unit::NTuple{80,UInt8}
	title::NTuple{80,UInt8}
	command::NTuple{320,UInt8}
	remark::NTuple{160,UInt8}
	# Items not stored in the data file for grids but explicitly used in macros computing node numbers
	nm::Csize_t 			# Number of data items in this grid (n_columns * n_rows) [padding is excluded]
	size::Csize_t 			# Actual number of items required to hold this grid (= mx * my)
	bits::UInt32
	complex_mode::UInt32
	type::UInt32
	n_bands::UInt32
	mx::UInt32 				# Actual dimensions of the grid in memory, allowing for the padding
	my::UInt32
	pad::NTuple{4,UInt32}
	mem_layout::NTuple{4,UInt8}
	nan_value::Cfloat
	xy_off::Cdouble
	ProjRefPROJ4::Ptr{UInt8}
	ProjRefWKT::Ptr{UInt8}
	ProjRefEPSG::Cint
	hidden::Ptr{Cvoid}		# Lower-level information for GMT use only
end

mutable struct GMT_GRID
	header::Ptr{GMT_GRID_HEADER}
	data::Ptr{Cfloat}
	x::Ptr{Cdouble}
	y::Ptr{Cdouble}
	hidden::Ptr{Cvoid}
end

mutable struct GMT_CUBE
	header::Ptr{GMT_GRID_HEADER}
	data::Ptr{Cfloat}
	x::Ptr{Cdouble}
	y::Ptr{Cdouble}
	hidden::Ptr{Cvoid}
	mode::UInt32
	z_range::NTuple{2,Cdouble}
	z_inc::Cdouble
	z::Ptr{Cdouble}
	name::NTuple{80,UInt8}
	units::NTuple{80,UInt8}
end

mutable struct GMT_DATASEGMENT
	n_rows::UInt64
	n_columns::UInt64
	min::Ptr{Cdouble}
	max::Ptr{Cdouble}
	data::Ptr{Ptr{Cdouble}}
	label::Ptr{UInt8}
	header::Ptr{UInt8}
	text::Ptr{Ptr{UInt8}}
	hidden::Ptr{Cvoid}
end

mutable struct GMT_DATATABLE
	n_headers::UInt32
	n_columns::UInt64
	n_segments::UInt64
	n_records::UInt64
	min::Ptr{Cdouble}
	max::Ptr{Cdouble}
	header::Ptr{Ptr{UInt8}}
	segment::Ptr{Ptr{GMT_DATASEGMENT}}
	hidden::Ptr{Cvoid}
end

mutable struct GMT_DATASET
	n_tables::UInt64
	n_columns::UInt64
	n_segments::UInt64
	n_records::UInt64
	min::Ptr{Cdouble}
	max::Ptr{Cdouble}
	table::Ptr{Ptr{GMT_DATATABLE}}
	type_::UInt32
	geometry::UInt32
	ProjRefPROJ4::Ptr{UInt8}
	ProjRefWKT::Ptr{UInt8}
	ProjRefEPSG::Cint
	hidden::Ptr{Cvoid}
end

struct GMT_FILL
	rgb::NTuple{4,Cdouble}
	f_rgb::NTuple{4,Cdouble}
	b_rgb::NTuple{4,Cdouble}
	use_pattern::Bool
	pattern_no::Int32
	dpi::UInt32
	pattern::NTuple{256,UInt8}		# was char pattern[GMT_BUFSIZ];
end

struct GMT_LUT
	z_low::Cdouble
	z_high::Cdouble
	i_dz::Cdouble
	rgb_low::NTuple{4,Cdouble}
	rgb_high::NTuple{4,Cdouble}
	rgb_diff::NTuple{4,Cdouble}
	hsv_low::NTuple{4,Cdouble}
	hsv_high::NTuple{4,Cdouble}
	hsv_diff::NTuple{4,Cdouble}
	annot::UInt32
	skip::UInt32
	fill::Ptr{GMT_FILL}
	label::Ptr{UInt8}
	key::Ptr{UInt8}
end

struct GMT_BFN
	rgb::NTuple{4,Cdouble}
	hsv::NTuple{4,Cdouble}
	skip::UInt32
	fill::Ptr{GMT_FILL}
end

mutable struct GMT_PALETTE
	data::Ptr{GMT_LUT}
	bfn::NTuple{3,GMT_BFN}
	n_headers::UInt32
	n_colors::UInt32
	mode::UInt32
	model::UInt32
	is_wrapping::UInt32
	is_gray::UInt32
	is_bw::UInt32
	is_continuous::UInt32
	has_pattern::UInt32
	has_hinge::UInt32
	has_range::UInt32
	categorical::UInt32
	minmax::NTuple{2,Cdouble}
	hinge::Cdouble
	wrap_length::Cdouble
	header::Ptr{Ptr{UInt8}}
	hidden::Ptr{Cvoid}
end

mutable struct GMT_IMAGE
	type::UInt32
	colormap::Ptr{Cint}
	n_indexed_colors::Cint
	header::Ptr{GMT_GRID_HEADER}
	data::Ptr{Cvoid}
	alpha::Ptr{Cuchar}
	color_interp::Ptr{UInt8}
	x::Ptr{Cdouble}
	y::Ptr{Cdouble}
	hidden::Ptr{Cvoid}
end

mutable struct GMT_POSTSCRIPT
	n_bytes::Csize_t
	mode::UInt32
	n_headers::UInt32
	data::Ptr{UInt8}
	header::Ptr{Ptr{UInt8}}
	hidden::Ptr{Cvoid}
end

struct GMT_UNIVECTOR
	uc1::Ptr{UInt8}
	sc1::Ptr{Int8}
	ui2::Ptr{UInt16}
	si2::Ptr{Int16}
	ui4::Ptr{UInt32}
	si4::Ptr{Int32}
	ui8::Ptr{UInt64}
	si8::Ptr{Int64}
	f4::Ptr{Float32}
	f8::Ptr{Float64}
end

struct GMT_VECTOR
	n_columns::UInt64
	n_rows::UInt64
	n_headers::UInt32
	registration::UInt32
	type::Ptr{UInt32}
	range::NTuple{2,Cdouble}
	# data::Ptr{GMT_UNIVECTOR}
	data::Ptr{Ptr{Cvoid}}
	text::Ptr{Ptr{UInt8}};			# Pointer to optional array of strings [NULL]
	command::NTuple{320,UInt8}
	remark::NTuple{160,UInt8}
	header::Ptr{Ptr{UInt8}};		# Array with all Vector header records, if any)
	ProjRefPROJ4::Ptr{UInt8}
	ProjRefWKT::Ptr{UInt8}
	ProjRefEPSG::Cint
	hidden::Ptr{Cvoid}
end

mutable struct GMT_MATRIX
	n_rows::UInt64
	n_columns::UInt64
	n_layers::UInt64
	n_headers::UInt32
	shape::UInt32
	registration::UInt32
	dim::Csize_t
	size::Csize_t
	type::UInt32
	range::NTuple{6,Cdouble}
	inc::NTuple{3,Cdouble}
#	data::GMT_UNIVECTOR
	data::Ptr{Cvoid}
	text::Ptr{Ptr{UInt8}}
	command::NTuple{320,UInt8}
	remark::NTuple{160,UInt8}
	header::Ptr{Ptr{UInt8}};		# Array with all Matrix header records, if any)
	ProjRefPROJ4::Ptr{UInt8}
	ProjRefWKT::Ptr{UInt8}
	ProjRefEPSG::Cint
	hidden::Ptr{Cvoid}
end

# Virtual file name length
VF_LEN = (GMTver < v"6.1") ? 16 : 32

mutable struct GMT_RESOURCE
	family::UInt32          # GMT data family, i.e., GMT_IS_DATASET, GMT_IS_GRID, etc.
	geometry::UInt32        # One of the recognized GMT geometries
	direction::UInt32       # Either GMT_IN or GMT_OUT
	option::Ptr{GMT_OPTION} # Pointer to the corresponding module option
	name::NTuple{VF_LEN,UInt8}  # Object ID returned by GMT_Register_IO
	pos::Cint               # Corresponding index into external object in|out arrays
	mode::Cint              # Either primary (0) or secondary (1) resource
	object::Ptr{Cvoid}      # Pointer to the actual GMT object
end

struct GMTAPI_DATA_OBJECT
	# Information for each input or output data entity, including information
	# needed while reading/writing from a table (file or array)
	rec::UInt64                 # Current rec to read [GMT_DATASET to/from MATRIX/VECTOR only]
	n_rows::UInt64              # Number or rows in this array [GMT_DATASET to/from MATRIX/VETOR only]
	n_columns::UInt64			# Number of columns to process in this dataset [GMT_DATASET only]
	n_expected_fields::UInt64	# Number of expected columns for this dataset [GMT_DATASET only]
	n_alloc::Csize_t			# Number of items allocated so far if writing to memory
	ID::UInt32					# Unique identifier which is >= 0
	alloc_level::UInt32			# Nested module level when object was allocated
	status::UInt32				# 0 when first registered, 1 after reading/writing has started, 2 when finished
	selected::Cint				# true if requested by current module, false otherwise
	close_file::Cint			# true if we opened source as a file and thus need to close it when done
	region::Cint				# true if wesn was passed, false otherwise
	no_longer_owner::Cint		# true if the data pointed to by the object was passed on to another object
	messenger::Cint				# true for output objects passed from the outside to receive data from GMT. If true we destroy data pointer before writing
	alloc_mode::UInt32			# GMT_ALLOCATED_{BY_GMT|EXTERNALLY}
	direction::UInt32			# GMT_IN or GMT_OUT
	family::UInt32				# One of GMT_IS_{DATASET|TEXTSET|CPT|IMAGE|GRID|MATRIX|VECTOR|COORD}
	actual_family::UInt32		# May be GMT_IS_MATRIX|VECTOR when one of the others are created via those
	method::UInt32              # One of GMT_IS_{FILE,STREAM,FDESC,DUPLICATE,REFERENCE} or sum with enum GMT_enum_via (GMT_VIA_{NONE,VECTOR,MATRIX,OUTPUT}); using unsigned type because sum exceeds enum GMT_enum_method
	geometry::UInt32			# One of GMT_IS_{POINT|LINE|POLY|PLP|SURFACE|NONE}
	wesn::NTuple{4,Cdouble}	# Grid domain limits
	resource::Ptr{Cvoid}			# Points to registered filename, memory location, etc., where data can be obtained from with GMT_Get_Data.
	data::Ptr{Cvoid}				# Points to GMT object that was read from a resource
	# FILE *fp;					# Pointer to source/destination stream [For rec-by-rec procession, NULL if memory location]
	fp::Ptr{Cvoid}				# Pointer to source/destination stream [For rec-by-rec procession, NULL if memory location]
	filename::Ptr{UInt8}		# Filename, stream, of file handle (otherwise NULL)
	# Cvoid *(*import) (struct GMT_CTRL *, FILE *, uint64_t *, int *);	# Pointer to input function (for DATASET/TEXTSET only)
	ifun::Ptr{Cvoid} 			# Pointer to input function (for DATASET/TEXTSET only)
	# Start of temporary variables for API debug - They are only set when building GMT with /DEBUG
	G::Ptr{Cvoid}				# struct GMT_GRID *G;
	D::Ptr{Cvoid}				# struct GMT_DATASET *D;
	T::Ptr{Cvoid}				# struct GMT_TEXTSET *T;
	C::Ptr{Cvoid}				# struct GMT_PALETTE *C;
	M::Ptr{Cvoid}				# struct GMT_MATRIX *M;
	V::Ptr{Cvoid}				# struct GMT_VECTOR *V;
	I::Ptr{Cvoid}				# struct GMT_IMAGE *I;
end

struct Gmt_libinfo
	name::Ptr{UInt8}	# Library tag name [without leading "lib" and extension], e.g. "gmt", "gmtsuppl" */
	path::Ptr{UInt8}	# Full path to library as given in GMT_CUSTOM_LIBS */
	skip::Ptr{Bool}		# true if we tried to open it and it was not available the first time */
	handle::Ptr{Cvoid}	# Handle to the shared library, returned by dlopen or dlopen_special */
end

if (GMTver >= v"6.1")
	mutable struct OGR_FEATURES
		n_rows::Cint
		n_cols::Cint
		n_layers::Cint
		n_filled::Cint
		is3D::Cint
		np::Cuint
		att_number::Cint
		n_islands::Cint				# Number of islands in a polygon
		name::Ptr{UInt8}
		wkt::Ptr{UInt8}
		proj4::Ptr{UInt8}
		type::Ptr{UInt8}            # Geometry type. E.g. Point, Polygon or LineString
		att_names::Ptr{Ptr{UInt8}}  # Names of the attributes of a Feature
		att_values::Ptr{Ptr{UInt8}} # Values of the attributes of a Feature as strings
		att_types::Ptr{Cint}
		islands::Ptr{Cint}
		BoundingBox::NTuple{6,Cdouble}
		BBgeom::Ptr{Cdouble};       # Not currently assigned (would be the BoundingBox of each individual geometry)
		x::Ptr{Cdouble}
		y::Ptr{Cdouble}
		z::Ptr{Cdouble}
	end
else
	mutable struct OGR_FEATURES
		n_rows::Cint
		n_cols::Cint
		n_layers::Cint
		n_filled::Cint
		is3D::Cint
		np::Cuint
		att_number::Cint
		name::Ptr{UInt8}
		wkt::Ptr{UInt8}
		proj4::Ptr{UInt8}
		type::Ptr{UInt8}            # Geometry type. E.g. Point, Polygon or LineString
		att_names::Ptr{Ptr{UInt8}}  # Names of the attributes of a Feature
		att_values::Ptr{Ptr{UInt8}} # Values of the attributes of a Feature as strings
		att_types::Ptr{Cint}
		islands::Ptr{Cint}
		BoundingBox::NTuple{6,Cdouble}
		BBgeom::Ptr{Cdouble};       # Not currently assigned (would be the BoundingBox of each individual geometry)
		x::Ptr{Cdouble}
		y::Ptr{Cdouble}
		z::Ptr{Cdouble}
	end
end

struct OGRREAD_CTRL
	info::Cint					# If != 0 gmt_ogrread will return only Info
	layer::Cint					# If > 0 will return only data from that layer. Use negative to return all layers
	name::Ptr{UInt8}			# Vector file name
	region::NTuple{6,Cdouble}	# For when a sub-region is required
end

struct COMMON_R
	active::NTuple{4,Bool}
	oblique::Bool
	registration::Csize_t
	row_oder::Cint
	mode::Cuint
	wesn::NTuple{6,Cdouble}
	wesn_orig::NTuple{4,Cdouble}
	inc::NTuple{2,Cdouble}
	string::NTuple{256,UInt8}
end

#=
mutable struct GMTAPI_CTRL
	# Master controller which holds all GMT API related information at run-time for a single session.
	# Users can run several GMT sessions concurrently; each session requires its own structure.
	# Use GMTAPI_Create_Session to initialize a new session and GMTAPI_Destroy_Session to end it.

	current_rec::NTuple{2,UInt64}	# Current record number >= 0 in the combined virtual dataset (in and out)
	n_objects::UInt32			# Number of currently active input and output data objects
	unique_ID::UInt32			# Used to create unique IDs for duration of session
	session_ID::UInt32			# ID of this session
	unique_var_ID::UInt32		# Used to create unique object IDs (grid,dataset, etc) for duration of session
	current_item::NTuple{2,Cint}	# Array number of current dataset being processed (in and out)
	pad::UInt32					# Session default for number of rows/cols padding for grids [2]
	external::UInt32			# 1 if called via external API (MATLAB, Python) [0]
	runmode::UInt32				# nonzero for GMT modern runmode [0 = classic]
	shape::Cint                 # GMT_IS_COL_FORMAT (1) if column-major (MATLAB, Fortran), GMT_IS_ROW_FORMAT (0) if row-major
	leave_grid_scaled::UInt32	# 1 if we dont want to unpack a grid after we packed it for writing [0]
	n_cores::UInt32             # Number of available cores on this system
	verbose::UInt32             # Used until GMT is set up
	registered::NTuple{2,Bool}	# true if at least one source/destination has been registered (in and out)
	io_enabled::NTuple{2,Bool}	# true if access has been allowed (in and out)
	module_input::Bool          # true when we are about to read inputs to the module (command line) */
	n_objects_alloc::Csize_t	# Allocation counter for data objects
	error::Int32				# Error code from latest API call [GMT_OK]
	last_error::Int32			# Error code from previous API call [GMT_OK]
	shelf::Int32				# Place to pass hidden values within API
	io_mode::NTuple{2,UInt32}	# 1 if access as set, 0 if record-by-record
	PPID::Cint                  # The Process ID of the parent (e.g., shell) or the external caller
	#GMT::Ptr{GMT_CTRL}			# Key structure with low-level GMT internal parameters
	GMT::Ptr{Cvoid}				# Maybe one day. Till than just keep it as Cvoid
	object::Ptr{Ptr{GMTAPI_DATA_OBJECT}}	# List of registered data objects
	session_tag::Ptr{UInt8}		# Name tag for this session (or NULL)
	tmp_dir::Ptr{UInt8}         # System tmp_dir (NULL if not found)
	gwf_dir::Ptr{UInt8}         # GMT WorkFlow dir (NULL if not running in modern mode)
	internal::Bool				# true if session was initiated by gmt.c
	deep_debug::Bool			# temporary for debugging
	#int (*print_func) (FILE *, const char *);	# Pointer to fprintf function (may be reset by external APIs like MEX)
	pf::Ptr{Cvoid}				# Don't know what to put here, so ley it be *void
	do_not_exit::UInt32			# 0 by default, mieaning it is OK to call exit  (may be reset by external APIs like MEX to call return instead)
	lib::Ptr{Gmt_libinfo}		# List of shared libs to consider
	n_shared_libs::UInt32		# How many in lib
end =#


text:
using FileIO
include("julia2ims.jl")

field = ["s1-6"  "s3-4"  "s7-12" "s8-3"  "s8-4"  "s12-1" "s17-9" "s14-7" "s14-5" "s15-11"    "s19-9" "s20-7" "s21-1" "s26-2" "s31-1" "s33-8" "s37-7" "s37-8" "s30-7" "s40-3" ]
data_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_hamamatsu-60X-TIRF/20200316";
output_dir = "$(data_dir)_visualization"

try
    mkdir("$output_dir")
catch
end

function max_projection(cell)
    println("Processing $cell")
    #raw_img = load(File(format"TIFF", "$data_dir/HE7-11-1-80uw-const_3_$pos.ome.tiff"));
    raw_img = loadims2("$(data_dir)_result/output/d16$cell.ims");
    img_size = size(raw_img);
    img_max = zeros(N0f16, img_size[1], img_size[2], Int(img_size[3]/20));
    for i in 1:Int(img_size[3]/20)
        img_max[:, :, i] = maximum(raw_img[:, :, (i-1)*20+1:i*20], dims=3);
    end
    save("$output_dir/$cell.tiff", img_max);
end

for cell in field
    max_projection(cell)
end


text:
struct ZSpec{Tscale<:Dict}
  originalms::MeasureSpec
  Zms::MeasureSpec
  scale::Tscale
end

#standardizes the data to make it easier to process
#WARNING- some field specific code here
function condition!(panel::AbstractDataFrame, ms::MeasureSpec)

  #scale the fields by the standard deviation
  #IMPORTANT: We need to use the same std dev on the lags in order to maintain integrity

  #assemble the new names and replace characteristics with standardized versions
  standardnames::Dict{NSymbol,NSymbol} = Dict{NSymbol,NSymbol}()
  standardnames[nothing]=nothing

  #create a dict to store the scale values
  scale::Dict{Symbol, Float64} = Dict()

  #start by standardizing the characteristics
  #NOTE: each characteristic weight and the lags, and controls are standardized by a single value
  #to avoid distorting changes in weights (so each characteristic uses a single standardizer)
  #QUESTION- could I use a seperate stnadardization value for the controls?
  Zξs::Vector{Characteristic} = Vector{Characteristic}()
  for ξ ∈ ms.ξs

    #IMPORTANT: use same scaling for lagged variables as unlagged equivelents
    standardnames[ξ.Fξ] = Symbol(:Z, ξ.Fξ)
    standardnames[ξ.X[:Fw]] = Symbol(:Z, ξ.X[:Fw])
    σ::Float64 = std(skipmissing(panel[!, ξ.X[:Fw]])) #use the clean panel to compute sds

    #now apply the same scaling to the lagged equivelents
    scale[standardnames[ξ.X[:Fw]]] = σ #do it this way to make sure we got the correct standard name
    for s ∈ [values(ξ.X)...; values(ms.FWbykt[ξ.Fξ])...;]
      standardnames[s] = Symbol(:Z, s)
      panel[!,standardnames[s]] = panel[!, s] ./ σ
      scale[standardnames[s]] = σ #record the scale
    end

    #construct the standardized characteristic
    #println("ZX: $(Dict(s=>standardnames[ξ.X[s]] for s ∈ keys(ξ.X)))")
    push!(Zξs, Characteristic(ξ, standardnames[ξ.Fξ],
      Dict(s=>standardnames[ξ.X[s]] for s ∈ keys(ξ.X))))
  end

  additionalfields::Vector{Symbol} = [ms.Fvol; ms.FW; ms.FWbyt]
  if !(ms.FDvol===nothing)
    push!(additionalfields, ms.FDvol)
  end

  #now standardize all other fields
  for F ∈ additionalfields
    standardnames[F] = Symbol(:Z, F)
    σ::Float64 = std(skipmissing(panel[!, F]))
    if σ + 0.01 ≈ 0.01 #check if σ is approaximately 0
      σ = 1.0#if so, leave the scale alone
    end
    scale[standardnames[F]] = σ
    panel[!, standardnames[F]] = panel[!, F] ./ σ
  end

  #create new names for the standardized controls
  ZFWbykt::Dict{Symbol, Vector{Symbol}} = Dict(standardnames[Fξ]=>(
    FW::Symbol->standardnames[FW]).(ms.FWbykt[Fξ])
    for Fξ ∈ keys(ms.FWbykt))
  Zcontrolnames = (;
    FW = (F->Symbol(:Z, F)).(ms.FW),
    FWbyt = (F->Symbol(:Z, F)).(ms.FWbyt),
    FWbykt = ZFWbykt,
  )

  #construct the standardized measurespec
  Zms = MeasureSpec(Zξs...,
    Fvol=standardnames[ms.Fvol],
    FDvol=standardnames[ms.FDvol],
    controlnames=Zcontrolnames)

  #check on the name integrity
  panelnames = propertynames(panel)
  for f ∈ [weightfields.(ms.ξs)...;(Fξ->ms.FWbykt[Fξ]).(ms.Fξs)...; additionalfields]
    (f ∈ panelnames) || error("$f not found in panel names")
    (Symbol(:Z, f) ∈ panelnames) || error("$(Symbol(:Z, f)) not found in panel names")
  end

  zs = ZSpec(ms, Zms, scale)

  #record the scaling values
  b = IOBuffer()
  #write(b, "Scaling info for $scale:")
  for p ∈ collect(scale)
    write(b, "$(p[1]) σ:\t $(p[2])\n")
  end
  standardizationinfo = String(take!(b))
  open("$(PARAM[:testpath])\\scaling\\$(iterresultsname()).txt", "w+") do f
    write(f, standardizationinfo)
  end

  return zs
end

function destandardize!(results::DataFrame, zs::ZSpec; intercept::Bool = true)
  #coefficients = intercept ? Fcoefficients : Fcoefficients0

  #this assumes the other controls are projected as opposed to calculated directly
  Fcoefficients(ms::MeasureSpec) = [ms.Fξs;#= ms.FW;=#]

  Fcoefs::Vector{Symbol} = Fcoefficients(zs.originalms)
  FZcoefs::Vector{Symbol} = Fcoefficients(zs.Zms)
  Fresults::Vector{Symbol} = propertynames(results)

  #make sure we have the expected coefficients
  if length(Fcoefs) ≠ length(FZcoefs) || length(setdiff(FZcoefs, Fresults)) ≠ 0
    error("Inconsistentcies found between Zspec and measure spec, and results:
      Fcoefs: $Fcoefs\nFZcoefs: $FZcoefs\nFresults: $Fresults")
  end

  #make sure the names are what we expect
  originalnames::Dict = Dict()
  for (Fcoef, FZcoef) ∈ zip(Fcoefs, FZcoefs)
    if FZcoef == Symbol(:Z, Fcoef)
      originalnames[FZcoef] = Fcoef
    else
      error("Unexpected field names found in Zms: \n Fcoefs: $Fcoefs \nFZcoefs: $FZcoefs")
    end
  end
  Symbol(:Z, zs.originalms.Fvol) == zs.Zms.Fvol || "Fvol and ZFvol have inconsistent names"
  originalnames[zs.Zms.Fvol] = zs.originalms.Fvol

  #now do the scaling
  @info "Beginning linear descaling of result fields $FZcoefs"
  #=NOTE: since V_nt/σᵥ =  sum_k 1/σ_kw * |ZA_tk * w_ntk - LZA_t * Lw_t*r_nt| + ε
    we have A_tk = σᵥ/σ_kw * ZA_tk=#

  σᵥ::Float64 = zs.scale[zs.Zms.Fvol]
  for (Fcoef, FZcoef) ∈ zip(Fcoefs, FZcoefs)
    resultcol::Vector = results[!, FZcoef]

    FZw::Symbol = zs.Zms[FZcoef].X[:Fw]
    σ_kw::Float64 = zs.scale[FZw]
    results[!, Fcoef] = resultcol .* σᵥ ./ σ_kw
  end

  standardizenames!(results, zs)
  return nothing

end

function standardizenames!(results::DataFrame, zs::ZSpec)
  Fcoefs::Vector{Symbol} = zs.originalms.Fξs
  FZcoefs::Vector{Symbol} = zs.Zms.Fξs
  FGcoefs::Vector{Symbol} = (s->Symbol(:G_,s)).(Fcoefs)

  @assert setdiff([:date; Fcoefs; FZcoefs; FGcoefs], propertynames(results)) |> isempty
  #setdiff(FGcoefs, propertynames(results)) |> isempty #"

#    Invalid setdiff: FGcoefs: $FGcoefs propertynames(results): $(propertynames(results))"
  #@assert setdiff(propertynames(results), [:date; Fcoefs; FZcoefs; FGcoefs]) |> isempty

  rename!(results, Dict(s=>Symbol(:A_,s) for s ∈ Fcoefs))
  rename!(results, Dict(s=>Symbol(:Z_,s) for s ∈ FZcoefs))

  return results
end

#some re-usable code regarding a variety of items needed to run the optimizatinos
function modelparts(panel::AbstractDataFrame,
    ::Type{T}=throw("3rd arg Type{T} is required!"),
    ::Type{TV}=throw("4th arg Type{TV} is required!"),
    ::Type{TM}=throw("5th arg Type{TM} is required!");
    zs::Union{ZSpec, Nothing} = nothing,
    ms::MeasureSpec = zs ≡ nothing ? throw("Need either zs::ZSpec or ms::MeasureSpec") : zs.Zms,
    weightdata::NamedTuple = error("weight data are required"),
    Fvol::Symbol = throw("Fvol is required"),
    FW::Vector{Symbol} = Vector{Symbol}(),
    FWgroup::Vector{Symbol} = Vector{Symbol}(),
    placebovol::Bool = PARAM[:placebovol],) where {
      T<:Real, TV<:AbstractVector{T}, TM<:AbstractMatrix{T}}

  #get some convenience quantities
  dims = getdims(panel, ms)
  tidx::Dict = gettidx(panel, ms)

  #get the time index for each date
  essentialfields::Vector{Symbol} = [[collect(values(weightdata))...;]...;
    FW...; FWgroup...; Fvol]
  cleanpanel::SubDataFrame = view(panel, completecases(panel, essentialfields), :)

  dat::Dict = Dict{Symbol, Any}()
  for k ∈ keys(weightdata)
    dat[k] = TM(reduce(hcat,(f->(T).(cleanpanel[!, f])).(weightdata[k])))
  end

  #test routine
  if placebovol
    dat[:Fvol] = (T).(cleanpanel[!, Fvol]) |> Random.shuffle |> TV
  else
    dat[:Fvol] = (T).(cleanpanel[!, Fvol]) |> TV
  end

  #acquire controls data. Creates an empty matrix if no controls
  emptymat = Matrix{T}(undef, size(cleanpanel,1), 0)
  dat[:FW] = reduce(hcat, (f->(T).(cleanpanel[!, f])).(FW), init=emptymat) |>TM
  dat[:FWgroup] = reduce(hcat, (f->(T).(cleanpanel[!, f])).(FWgroup), init=emptymat) |>TM
  dat[:permno] = cleanpanel[!, :permno] |> deepcopy
  dat[:date] = cleanpanel.date |> deepcopy

  ts::Vector{Int} = (d::Date->tidx[d]).(cleanpanel[!,:date])
  @assert issorted(ts)
  #V::TV = (T).(cleanpanel[!,ms.Fvol]) |> fluxgpu

  #now get the scaling values
  #Ascale = Vector{T}(undef, dims.K)
  #@assert length(Ascale) == size(dat[:Fw],2) #integrity check
  #FZcoefs::Vector{Symbol} = Fcoefficients0(zs.Zms)
  #@assert length(FZcoefs)  == length(Ascale)

  #=
  σᵥ::Float64 = zs.scale[zs.Zms.Fvol]
  for (i,FZcoef) ∈ enumerate(FZcoefs) #loop through each coefficient and identify the scaling value
    FZw::Symbol = zs.Zms[FZcoef].X[:Fw] #pull the standard weighting name (this is keyed to the scale)
    σ_kw::Float64 = zs.scale[FZw]
    Ascale[i] = σᵥ / σ_kw
  end

  @info "Scaling Values: (σᵥ, $σᵥ), $((zip((FZcoef->zs.Zms[FZcoef].X[:Fw]).(FZcoefs), Ascale)...))"
  =#
  return (dims=dims, tidx=tidx, cleanpanel=cleanpanel,
    dat=dat, ts=ts#=, Ascale=Ascale=#)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Trust-region
# ------------------------------------------------------------
const TRConstraint = CI{VectorAffineDecisionFunction{Float64}, MOI.NormInfinityCone}
@with_kw mutable struct TRData{T <: AbstractFloat}
    Q̃::T = 1e10
    Δ::MOI.VariableIndex = MOI.VariableIndex(0)
    constraint::TRConstraint = TRConstraint(0)
    cΔ::Int = 0
    incumbent::Int = 1
    major_iterations::Int = 0
    minor_iterations::Int = 0
end

@with_kw mutable struct TRParameters{T <: AbstractFloat}
    γ::T = 1e-4
    Δ::T = 1.0
    Δ̅::T = 1000.0
end

"""
    TrustRegion

Functor object for using trust-region regularization in an L-shaped algorithm. Create by supplying a [`TR`](@ref) object through `regularize` in `LShaped.Optimizer` or by setting the [`Regularizer`](@ref) attribute.

...
# Parameters
- `γ::T = 1e-4`: Relative tolerance for deciding if a minor iterate should be accepted as a new major iterate.
- `Δ::AbstractFloat = 1.0`: Initial size of ∞-norm trust-region.
- `Δ̅::AbstractFloat = 1000.0`: Maximum size of ∞-norm trust-region.
...
"""
struct TrustRegion{T <: AbstractFloat, A <: AbstractVector} <: AbstractRegularization
    data::TRData{T}
    parameters::TRParameters{T}

    decisions::DecisionMap
    projection_targets::Vector{MOI.VariableIndex}
    ξ::Vector{Decision{T}}

    Q̃_history::A
    Δ_history::A
    incumbents::Vector{Int}

    function TrustRegion(decisions::DecisionMap, ξ₀::AbstractVector; kw...)
        T = promote_type(eltype(ξ₀), Float32)
        A = Vector{T}
        ξ = map(ξ₀) do val
            KnownDecision(val, T)
        end
        return new{T, A}(TRData{T}(),
                         TRParameters{T}(; kw...),
                         decisions,
                         Vector{MOI.VariableIndex}(undef, length(ξ₀)),
                         ξ,
                         A(),
                         A(),
                         Vector{Int}())
    end
end

function initialize_regularization!(lshaped::AbstractLShaped, tr::TrustRegion{T}) where T <: AbstractFloat
    n = length(tr.ξ) + 1
    # Sense-correct initial objective
    sense = MOI.get(lshaped.master, MOI.ObjectiveSense())
    coeff = sense == MOI.MIN_SENSE ? 1.0 : -1.0
    tr.data.Q̃ *= coeff
    # Add projection targets
    add_projection_targets!(tr, lshaped.master)
    # Add trust region
    name = string(:Δ)
    trust_region = KnownDecision(tr.parameters.Δ, T)
    set = SingleDecisionSet(1, trust_region, NoSpecifiedConstraint(), false)
    tr.data.Δ, _ = MOI.add_constrained_variable(lshaped.master, set)
    set_decision!(tr.decisions, tr.data.Δ, trust_region)
    MOI.set(lshaped.master, MOI.VariableName(), tr.data.Δ, name)
    x = VectorOfDecisions(all_decisions(tr.decisions))
    ξ = VectorOfDecisions(tr.projection_targets)
    Δ = SingleDecision(tr.data.Δ)
    # Add trust-region constraint
    f = MOIU.operate(vcat, T, Δ, x) -
        MOIU.operate(vcat, T, zero(tr.parameters.Δ), ξ)
    tr.data.constraint =
        MOI.add_constraint(lshaped.master, f,
                           MOI.NormInfinityCone(n))
    return nothing
end

function restore_regularized_master!(lshaped::AbstractLShaped, tr::TrustRegion)
    # Delete trust region constraint
    if !iszero(tr.data.constraint.value)
        MOI.delete(lshaped.master, tr.data.constraint)
        tr.data.constraint = TRConstraint(0)
    end
    # Delete trust region
    if !iszero(tr.data.Δ.value)
        MOI.delete(lshaped.master, tr.data.Δ)
        StochasticPrograms.remove_decision!(tr.decisions, tr.data.Δ)
        tr.data.Δ = MOI.VariableIndex(0)
    end
    # Delete projection targets
    for var in tr.projection_targets
        MOI.delete(lshaped.master, var)
        StochasticPrograms.remove_decision!(tr.decisions, var)
    end
    empty!(tr.projection_targets)
    return nothing
end

function filter_variables!(tr::TrustRegion, list::Vector{MOI.VariableIndex})
    # Filter projection targets
    filter!(vi -> !(vi in tr.projection_targets), list)
    # Filter Δ
    i = something(findfirst(isequal(tr.data.Δ), list), 0)
    if !iszero(i)
        MOI.deleteat!(list, i)
    end
    return nothing
end

function filter_constraints!(tr::TrustRegion, list::Vector{<:CI})
    # Filter trust-region constraint
    i = something(findfirst(isequal(tr.data.constraint), list), 0)
    if !iszero(i)
        MOI.deleteat!(list, i)
    end
    return nothing
end

function log_regularization!(lshaped::AbstractLShaped, tr::TrustRegion)
    @unpack Q̃, Δ, incumbent = tr.data
    push!(tr.Q̃_history, Q̃)
    push!(tr.Δ_history, StochasticPrograms.decision(tr.decisions, Δ).value)
    push!(tr.incumbents, incumbent)
    return nothing
end

function log_regularization!(lshaped::AbstractLShaped, t::Integer, tr::TrustRegion)
    @unpack Q̃,Δ,incumbent = tr.data
    tr.Q̃_history[t] = Q̃
    tr.Δ_history[t] = StochasticPrograms.decision(tr.decisions, Δ).value
    tr.incumbents[t] = incumbent
    return nothing
end

function take_step!(lshaped::AbstractLShaped, tr::TrustRegion)
    @unpack Q,θ = lshaped.data
    @unpack τ = lshaped.parameters
    @unpack Q̃ = tr.data
    @unpack γ = tr.parameters
    need_update = false
    t = timestamp(lshaped)
    Q̃t = incumbent_objective(lshaped, t, tr)
    sense = MOI.get(lshaped.master, MOI.ObjectiveSense())
    coeff = sense == MOI.MIN_SENSE ? 1.0 : -1.0
    if coeff*Q + τ <= coeff*Q̃ && (tr.data.major_iterations == 1 || coeff*Q <= coeff*Q̃t - γ*abs(Q̃t-θ))
        need_update = true
        enlarge_trustregion!(lshaped, tr)
        tr.data.cΔ = 0
        x = current_decision(lshaped)
        for i in eachindex(tr.ξ)
            tr.ξ[i].value = x[i]
        end
        tr.data.Q̃ = Q
        tr.data.incumbent = timestamp(lshaped)
        tr.data.major_iterations += 1
    else
        need_update = reduce_trustregion!(lshaped, tr)
        tr.data.minor_iterations += 1
    end
    if need_update
        update_trustregion!(lshaped, tr)
    end
    return nothing
end

function process_cut!(lshaped::AbstractLShaped, cut::HyperPlane{FeasibilityCut}, tr::TrustRegion)
    @unpack τ = lshaped.parameters
    if !satisfied(cut, decision(lshaped), τ)
        # Project decision to ensure prevent master infeasibility
        A = [I cut.δQ; cut.δQ' 0*I]
        b = [zeros(length(tr.ξ)); -gap(cut, decision(lshaped))]
        t = A\b
        for i in eachindex(tr.ξ)
            tr.ξ[i].value += t[i]
        end
        update_trustregion!(lshaped, tr)
    end
    return nothing
end

function update_trustregion!(lshaped::AbstractLShaped, tr::TrustRegion)
    @unpack Δ = tr.data
    # Update projection targets
    for vi in tr.projection_targets
        ci = CI{MOI.SingleVariable,SingleDecisionSet{Float64}}(vi.value)
        MOI.modify(lshaped.master,
                   ci,
                   KnownValuesChange())
    end
    # Update trust-region
    ci = CI{MOI.SingleVariable,SingleDecisionSet{Float64}}(Δ.value)
    MOI.modify(lshaped.master,
               ci,
               KnownValuesChange())
    return nothing
end

function enlarge_trustregion!(lshaped::AbstractLShaped, tr::TrustRegion)
    @unpack Q,θ = lshaped.data
    @unpack τ, = lshaped.parameters
    @unpack Δ = tr.data
    @unpack Δ̅ = tr.parameters
    t = timestamp(lshaped)
    Δ̃ = incumbent_trustregion(lshaped, t, tr)
    ξ = incumbent_decision(lshaped, t, tr)
    Q̃ = incumbent_objective(lshaped, t, tr)
    sense = MOI.get(lshaped.master, MOI.ObjectiveSense())
    coeff = sense == MOI.MIN_SENSE ? 1.0 : -1.0
    if coeff*(Q̃ - Q) >= 0.5*coeff*(Q̃ - θ) && abs(norm(ξ - lshaped.x, Inf) - Δ̃) <= τ
        # Enlarge the trust-region radius
        Δ = StochasticPrograms.decision(tr.decisions, tr.data.Δ)
        Δ.value = max(Δ.value, min(Δ̅, 2 * Δ̃))
        return true
    else
        return false
    end
end

function reduce_trustregion!(lshaped::AbstractLShaped, tr::TrustRegion)
    @unpack Q,θ = lshaped.data
    @unpack Q̃,Δ,cΔ = tr.data
    t = timestamp(lshaped)
    Δ̃ = incumbent_trustregion(lshaped, t, tr)
    Q̃ = incumbent_objective(lshaped, t, tr)
    ρ = min(1, Δ̃)*(Q-Q̃)/(Q̃-θ)
    if ρ > 0
        tr.data.cΔ += 1
    end
    if ρ > 3 || (cΔ >= 3 && 1 < ρ <= 3)
        # Reduce the trust-region radius
        tr.data.cΔ = 0
        Δ = StochasticPrograms.decision(tr.decisions, tr.data.Δ)
        Δ.value = min(Δ.value, (1/min(ρ,4))*Δ̃)
        return true
    else
        return false
    end
end

# API
# ------------------------------------------------------------
"""
    TR

Factory object for [`TrustRegion`](@ref). Pass to `regularize` in `LShaped.Optimizer` or set the [`Regularizer`](@ref) attribute.. Equivalent factory calls: `TR`, `WithTR`, `TrustRegion`, `WithTrustRegion`. See ?TrustRegion for parameter descriptions.

"""
mutable struct TR <: AbstractRegularizer
    parameters::TRParameters{Float64}
end
TR(; kw...) = TR(TRParameters(; kw...))
WithTR(; kw...) = TR(TRParameters(; kw...))
TrustRegion(; kw...) = TR(TRParameters(; kw...))
WithTrustRegion(; kw...) = TR(TRParameters(; kw...))

function (tr::TR)(decisions::DecisionMap, x::AbstractVector)
    return TrustRegion(decisions, x; type2dict(tr.parameters)...)
end

function str(::TR)
    return "L-shaped using trust-region"
end


text:
"""
Naive index type for storing vectors. It is a wrapper
around a vector of embeddings and performs brute search using
the cosine similarity between vectors.
"""
struct NaiveIndex{E, A<:Vector{<:AbstractVector{E}}} <: Garamond.AbstractIndex
    data::A
end

NaiveIndex(data, args...; kwargs...) =
    NaiveIndex([data[:, i] for i in 1:size(data, 2)])  # args, kwargs are ignored

# Nearest neighbor search method
function knn_search(index::NaiveIndex{E},
                    point::AbstractVector,
                    k::Int,
                    keep::Vector{Int}=collect(1:length(index));
                    kwargs...
                   ) where {E<:AbstractFloat}
    n = length(keep)
    _k = min(n, k)
    scores = zeros(E, n)
    _point = densify(point)'
    # Cosine similarity
    @inbounds for i in eachindex(keep)
        scores[i] = _point * index.data[i]
    end
    #scores = index.data[:, keep]' * __densify(point)
    idxs = sortperm(scores, rev=true)[1:_k]
    return keep[idxs], one(E) .- scores[idxs]
end

# Length method
length(index::NaiveIndex) = length(index.data)

# push!, pushfirst!, pop!, popfirst!, delete_from_index!
Base.push!(index::NaiveIndex, point) = push!(index.data, point)

Base.pushfirst!(index::NaiveIndex, point) = pushfirst!(index.data, point)

Base.pop!(index::NaiveIndex) = pop!(index.data)

Base.popfirst!(index::NaiveIndex) = popfirst!(index.data)

delete_from_index!(index::NaiveIndex, pos) = begin
    deleteat!(index.data, pos)
    nothing
end


text:
using Documenter, VIEstimators

makedocs(;
    modules=[VIEstimators],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/willtebbutt/VIEstimators.jl/blob/{commit}{path}#L{line}",
    sitename="VIEstimators.jl",
    authors="Will Tebbutt <wt0881@my.bristol.ac.uk>",
    assets=String[],
)

deploydocs(;
    repo="github.com/willtebbutt/VIEstimators.jl",
)


text:
module AARMBEM

############################################################################################
module MolModule

using LinearAlgebra;
using DelimitedFiles;
using SharedArrays;

export OneMol, OneMol_WithPhonons, OneMol_NoPhonons, MolSystem, readMolSystem
## export extractKIconst!, extractRequiredMolDataVec, extractOptionalMolDataVec
## export extractRequiredMolDataMat, extractBescale, extractGeneralMolData
## export addKIconsttoKIendatoms!
include("MolData.jl");

end
############################################################################################

############################################################################################
module GenericMatrixVectorUtils

using LinearAlgebra;

include("GenericMatrixVectorUtils.jl");
export addtodiag!, addtodiag, subtractfromdiag!, subtractfromdiag, mulMatDiagleft!
export mulMatDiagleft, mulMatDiagright!, mulMatDiagright, vecNto3N, sympart!, sympart
export hermpart!, hermpart, ahermpart!, ahermpart

end
############################################################################################

############################################################################################
module PeriodicModule

using LinearAlgebra;
using DelimitedFiles;

include("PeriodicData.jl");
export PeriodicData, readPeriodicData, constructlatticefromreciprocal!
export constructreciprocalfromlattice!

end
############################################################################################

############################################################################################
module GFGG

using LinearAlgebra;
using SpecialFunctions;

using ..GenericMatrixVectorUtils;
export addtodiag!, addtodiag, subtractfromdiag!, subtractfromdiag, mulMatDiagleft!
export mulMatDiagleft, mulMatDiagright!, mulMatDiagright, vecNto3N, sympart!, sympart
export hermpart!, hermpart, ahermpart!, ahermpart

include("E1func.jl");
include("GFVACGG.jl");
export c, GFVACGG!, GFVACGG_sca!, GFVACGG0!, GFVACGG0_sca!, GFVACGGcoincident!
export GFVACGG1Ewald!, GFVACGG2Ewald!, GFVACGGEwaldSR, GFVACGGEwaldLR1D, GFVACGGEwaldLR2D
export GFVACGG, GFVACGG_sca, GFVACGG0, GFVACGG0_sca, GFVACGGcoincident, GFVACGG1Ewald
export GFVACGG2Ewald

include("GFPECGG.jl");
export GFPECGG!, GFPECGG0!, GFPECGG1Ewald!, GFPECGG2Ewald!, GFPECGG, GFPECGG0
export GFPECGG1Ewald, GFPECGG2Ewald

end
############################################################################################

############################################################################################
module MolTrans

using LinearAlgebra;
using DelimitedFiles;
using SharedArrays;

using ..MolModule;
export OneMol, OneMol_WithPhonons, OneMol_NoPhonons, MolSystem, readMolSystem

include("MolTransUtils.jl");
export findCOM, translateCOM!, rotatecoords!, translateCOM, rotatecoords, rotate3Nmat!
export rotate3Nmat

include("TransData.jl");
export TransData, readTransData, initializeAtomPos!, transformAtomPos!, untransformAtomPos!

end
############################################################################################

############################################################################################
module MolAlpha

using LinearAlgebra;
using DelimitedFiles;
using SharedArrays;

using ..MolModule;
export OneMol, OneMol_WithPhonons, OneMol_NoPhonons, MolSystem, readMolSystem

using ..GenericMatrixVectorUtils;
export addtodiag!, addtodiag, subtractfromdiag!, subtractfromdiag, mulMatDiagleft!
export mulMatDiagleft, mulMatDiagright!, mulMatDiagright, vecNto3N, sympart!, sympart
export hermpart!, hermpart, ahermpart!, ahermpart

using ..PeriodicModule;
export PeriodicData, readPeriodicData, constructlatticefromreciprocal!
export constructreciprocalfromlattice!

include("MolAlpha.jl");
export constructMolAlpha!, KIk

end
############################################################################################

############################################################################################
module MolEM

using LinearAlgebra;
using DelimitedFiles;
using SharedArrays;
using Distributed;

using ..MolModule;
export OneMol, OneMol_WithPhonons, OneMol_NoPhonons, MolSystem, readMolSystem

include("MolTransUtils.jl");
export findCOM, translateCOM!, rotatecoords!, translateCOM, rotatecoords, rotate3Nmat!
export rotate3Nmat

include("TransData.jl");
export TransData, readTransData, initializeAtomPos!, transformAtomPos!, untransformAtomPos!

using ..GenericMatrixVectorUtils;
export addtodiag!, addtodiag, subtractfromdiag!, subtractfromdiag, mulMatDiagleft!
export mulMatDiagleft, mulMatDiagright!, mulMatDiagright, vecNto3N, sympart!, sympart
export hermpart!, hermpart, ahermpart!, ahermpart

using ..PeriodicModule;
export PeriodicData, readPeriodicData, constructlatticefromreciprocal!
export constructreciprocalfromlattice!

include("MolAlpha.jl");
export constructMolAlpha!, KIk

using ..GFGG;
export c, GFVACGG!, GFVACGG_sca!, GFVACGG0!, GFVACGG0_sca!, GFVACGGcoincident!
export GFVACGG1Ewald!, GFVACGG2Ewald!, GFVACGGEwaldSR, GFVACGGEwaldLR1D, GFVACGGEwaldLR2D
export GFVACGG, GFVACGG_sca, GFVACGG0, GFVACGG0_sca, GFVACGGcoincident, GFVACGG1Ewald
export GFVACGG2Ewald
export GFPECGG!, GFPECGG0!, GFPECGG1Ewald!, GFPECGG2Ewald!, GFPECGG, GFPECGG0
export GFPECGG1Ewald, GFPECGG2Ewald

include("MolAssembleGFGG.jl");
export assemble1MolGvacinf!, assemble1MolGvacinfdiag!, subtract1MolalphainvGvac!
export assembleGscaMolBlock!, assembleGvacMolBlock!

end
############################################################################################

############################################################################################
module AARMBEMCMDARGS

using LinearAlgebra;
using DelimitedFiles;
using Printf;

include("AARMBEM_ReadCMDARGS.jl");
export readRequiredFilenameARG, readOptionalFilenameARG, readOptionalBoolcondARG
export readOptionalFloatType, readfreqlist!, readklist, readGenvstr

end
############################################################################################

############################################################################################
module AARMBEMvdW

using LinearAlgebra;
using DelimitedFiles;
using SharedArrays;
using Distributed;
using Printf;

using ..MolEM;
using ..AARMBEMCMDARGS;

include("vdWenergy.jl");
export vdWenergy

end
############################################################################################

############################################################################################
module AARMBEMheat

using LinearAlgebra;
using DelimitedFiles;
using SharedArrays;
using Distributed;
using Printf;

using ..MolEM;
using ..AARMBEMCMDARGS;

include("heattransfer.jl");
export heattransfer

end
############################################################################################

end


text:
function get_test_roadway()
    roadway = Roadway()

    seg1 = RoadSegment(1,
        [Lane(LaneTag(1,1),
            [CurvePt(VecSE2(-2.0,0.0,0.0), 0.0),
             CurvePt(VecSE2(-1.0,0.0,0.0), 1.0)]),
         Lane(LaneTag(1,2),
            [CurvePt(VecSE2(-2.0,1.0,0.0), 0.0),
             CurvePt(VecSE2(-1.0,1.0,0.0), 1.0)]),
         Lane(LaneTag(1,3),
            [CurvePt(VecSE2(-2.0,2.0,0.0), 0.0),
             CurvePt(VecSE2(-1.0,2.0,0.0), 1.0)])
        ])

    seg2 = RoadSegment(2,
        [Lane(LaneTag(2,1),
            [CurvePt(VecSE2(0.0,0.0,0.0), 0.0),
             CurvePt(VecSE2(1.0,0.0,0.0), 1.0),
             CurvePt(VecSE2(2.0,0.0,0.0), 2.0),
             CurvePt(VecSE2(3.0,0.0,0.0), 3.0)]),
         Lane(LaneTag(2,2),
            [CurvePt(VecSE2(0.0,1.0,0.0), 0.0),
             CurvePt(VecSE2(1.0,1.0,0.0), 1.0),
             CurvePt(VecSE2(2.0,1.0,0.0), 2.0),
             CurvePt(VecSE2(3.0,1.0,0.0), 3.0)]),
         Lane(LaneTag(2,3),
            [CurvePt(VecSE2(0.0,2.0,0.0), 0.0),
             CurvePt(VecSE2(1.0,2.0,0.0), 1.0),
             CurvePt(VecSE2(2.0,2.0,0.0), 2.0),
             CurvePt(VecSE2(3.0,2.0,0.0), 3.0)])
         ])

    seg3 = RoadSegment(3,
        [Lane(LaneTag(3,1),
            [CurvePt(VecSE2(4.0,0.0,0.0), 0.0),
             CurvePt(VecSE2(5.0,0.0,0.0), 1.0)]),
         Lane(LaneTag(3,2),
            [CurvePt(VecSE2(4.0,1.0,0.0), 0.0),
             CurvePt(VecSE2(5.0,1.0,0.0), 1.0)]),
         Lane(LaneTag(3,3),
            [CurvePt(VecSE2(4.0,2.0,0.0), 0.0),
             CurvePt(VecSE2(5.0,2.0,0.0), 1.0)])
        ])

    for i in 1:3
        connect!(seg1.lanes[i], seg2.lanes[i])
        connect!(seg2.lanes[i], seg3.lanes[i])
    end

    push!(roadway.segments, seg1)
    push!(roadway.segments, seg2)
    push!(roadway.segments, seg3)

    roadway
end

let
    curve = get_test_curve1()

    lanetag = LaneTag(1,1)
    lane = Lane(lanetag, curve)
    show(IOBuffer(), lanetag)

    @test !has_next(lane)
    @test !has_prev(lane)

    show(IOBuffer(), RoadIndex(CurveIndex(1,0.0), lanetag))

    roadway = get_test_roadway()

    lane = roadway[LaneTag(1,1)]
    @test has_next(lane)
    @test !has_prev(lane)
    @test lane.exits[1].target == RoadIndex(CurveIndex(1,0.0), LaneTag(2,1))

    lane = roadway[LaneTag(2,1)]
    @test has_next(lane)
    @test has_prev(lane)
    @test lane.exits[1].target == RoadIndex(CurveIndex(1,0.0), LaneTag(3,1))
    @test lane.entrances[1].target == RoadIndex(CurveIndex(1,1.0), LaneTag(1,1))

    res = proj(VecSE2(1.0,0.0,0.0), lane, roadway)
    @test res.curveproj.ind == CurveIndex(2, 0.0)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == lane.tag

    res = proj(VecSE2(1.0,0.0,0.0), lane, roadway, move_along_curves=false)
    @test res.curveproj.ind == CurveIndex(2, 0.0)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == lane.tag

    res = proj(VecSE2(1.5,0.25,0.1), lane, roadway)
    @test res.curveproj.ind == CurveIndex(2, 0.5)
    @test isapprox(res.curveproj.t, 0.25)
    @test isapprox(res.curveproj.ϕ, 0.1)
    @test res.tag == lane.tag

    res = proj(VecSE2(0.0,0.0,0.0), lane, roadway)
    @test res.curveproj.ind == CurveIndex(1, 0.0)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == lane.tag

    res = proj(VecSE2(-0.75,0.0,0.0), lane, roadway)
    @test res.curveproj.ind == CurveIndex(0, 0.25)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == lane.tag

    res = proj(VecSE2(-1.75,0.0,0.0), lane, roadway)
    @test res.curveproj.ind == CurveIndex(1, 0.25)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == prev_lane(lane, roadway).tag

    res = proj(VecSE2(-1.75,0.0,0.0), lane, roadway, move_along_curves=false)
    @test res.curveproj.ind == CurveIndex(0, 0.0)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == lane.tag

    res = proj(VecSE2(4.25,0.2,0.1), lane, roadway)
    @test res.curveproj.ind == CurveIndex(1, 0.25)
    @test isapprox(res.curveproj.t, 0.2)
    @test isapprox(res.curveproj.ϕ, 0.1)
    @test res.tag == next_lane(lane, roadway).tag

    res = proj(VecSE2(4.25,0.2,0.1), lane, roadway, move_along_curves=false)
    @test res.curveproj.ind == CurveIndex(4, 1.0)
    @test isapprox(res.curveproj.t, 0.2)
    @test isapprox(res.curveproj.ϕ, 0.1)
    @test res.tag == lane.tag

    res = proj(VecSE2(3.25,0.2,0.1), lane, roadway)
    @test res.curveproj.ind == CurveIndex(4, 0.25)
    @test isapprox(res.curveproj.t, 0.2)
    @test isapprox(res.curveproj.ϕ, 0.1)
    @test res.tag == lane.tag

    res = proj(VecSE2(4.25,0.2,0.1), prev_lane(lane, roadway), roadway)
    @test res.curveproj.ind == CurveIndex(1, 0.25)
    @test isapprox(res.curveproj.t, 0.2)
    @test isapprox(res.curveproj.ϕ, 0.1)
    @test res.tag == next_lane(lane, roadway).tag

    res = proj(VecSE2(-0.75,0.0,0.0), next_lane(lane, roadway), roadway)
    @test res.curveproj.ind == CurveIndex(0, 0.25)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == lane.tag

    ####

    seg = roadway[2]

    res = proj(VecSE2(1.0,0.0,0.0), seg, roadway)
    @test res.curveproj.ind == CurveIndex(2, 0.0)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == lane.tag

    res = proj(VecSE2(1.5,0.25,0.1), seg, roadway)
    @test res.curveproj.ind == CurveIndex(2, 0.5)
    @test isapprox(res.curveproj.t, 0.25)
    @test isapprox(res.curveproj.ϕ, 0.1)
    @test res.tag == lane.tag

    res = proj(VecSE2(0.0,0.0,0.0), seg, roadway)
    @test res.curveproj.ind == CurveIndex(1, 0.0)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == lane.tag

    res = proj(VecSE2(-0.75,0.0,0.0), seg, roadway)
    @test res.curveproj.ind == CurveIndex(0, 0.25)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == lane.tag

    res = proj(VecSE2(-1.75,0.0,0.0), seg, roadway)
    @test res.curveproj.ind == CurveIndex(1, 0.25)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == prev_lane(lane, roadway).tag

    res = proj(VecSE2(4.25,0.2,0.1), seg, roadway)
    @test res.curveproj.ind == CurveIndex(1, 0.25)
    @test isapprox(res.curveproj.t, 0.2)
    @test isapprox(res.curveproj.ϕ, 0.1)
    @test res.tag == next_lane(lane, roadway).tag

    res = proj(VecSE2(3.25,0.2,0.1), seg, roadway)
    @test res.curveproj.ind == CurveIndex(4, 0.25)
    @test isapprox(res.curveproj.t, 0.2)
    @test isapprox(res.curveproj.ϕ, 0.1)
    @test res.tag == lane.tag

    res = proj(VecSE2(4.25,0.2,0.1), roadway[1], roadway)
    @test res.curveproj.ind == CurveIndex(1, 0.25)
    @test isapprox(res.curveproj.t, 0.2)
    @test isapprox(res.curveproj.ϕ, 0.1)
    @test res.tag == next_lane(lane, roadway).tag

    res = proj(VecSE2(-0.75,0.0,0.0), roadway[3], roadway)
    @test res.curveproj.ind == CurveIndex(0, 0.25)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == lane.tag

    res = proj(VecSE2(1.0,1.0,0.0), seg, roadway)
    @test res.curveproj.ind == CurveIndex(2, 0.0)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == LaneTag(2,2)

    res = proj(VecSE2(1.0,2.0,0.0), seg, roadway)
    @test res.curveproj.ind == CurveIndex(2, 0.0)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == LaneTag(2,3)

    ###

    res = proj(VecSE2(1.0,0.0,0.0), roadway)
    @test res.curveproj.ind == CurveIndex(2, 0.0)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == lane.tag

    res = proj(VecSE2(1.5,0.25,0.1), roadway)
    @test res.curveproj.ind == CurveIndex(2, 0.5)
    @test isapprox(res.curveproj.t, 0.25)
    @test isapprox(res.curveproj.ϕ, 0.1)
    @test res.tag == lane.tag

    res = proj(VecSE2(0.0,0.0,0.0), roadway)
    @test res.curveproj.ind == CurveIndex(1, 0.0)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == lane.tag

    res = proj(VecSE2(-0.75,0.0,0.0), roadway)
    @test res.curveproj.ind == CurveIndex(2, 0.25)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == prev_lane(lane, roadway).tag

    res = proj(VecSE2(-1.75,0.0,0.0), roadway)
    @test res.curveproj.ind == CurveIndex(1, 0.25)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == prev_lane(lane, roadway).tag

    res = proj(VecSE2(4.25,0.2,0.1), roadway)
    @test res.curveproj.ind == CurveIndex(1, 0.25)
    @test isapprox(res.curveproj.t, 0.2)
    @test isapprox(res.curveproj.ϕ, 0.1)
    @test res.tag == next_lane(lane, roadway).tag

    res = proj(VecSE2(3.25,0.2,0.1), roadway)
    @test res.curveproj.ind == CurveIndex(4, 0.25)
    @test isapprox(res.curveproj.t, 0.2)
    @test isapprox(res.curveproj.ϕ, 0.1)
    @test res.tag == lane.tag

    res = proj(VecSE2(4.25,0.2,0.1), roadway[1], roadway)
    @test res.curveproj.ind == CurveIndex(1, 0.25)
    @test isapprox(res.curveproj.t, 0.2)
    @test isapprox(res.curveproj.ϕ, 0.1)
    @test res.tag == next_lane(lane, roadway).tag

    res = proj(VecSE2(-0.75,0.0,0.0), roadway[3], roadway)
    @test res.curveproj.ind == CurveIndex(0, 0.25)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == lane.tag

    res = proj(VecSE2(1.0,1.0,0.0), roadway)
    @test res.curveproj.ind == CurveIndex(2, 0.0)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == LaneTag(2,2)

    res = proj(VecSE2(1.0,2.0,0.0), roadway)
    @test res.curveproj.ind == CurveIndex(2, 0.0)
    @test isapprox(res.curveproj.t, 0.0)
    @test isapprox(res.curveproj.ϕ, 0.0)
    @test res.tag == LaneTag(2,3)

    ###

    roadind_0 = RoadIndex(CurveIndex(1,0), LaneTag(2,1))
    roadind = move_along(roadind_0, roadway, 0.0)
    @test roadind == roadind_0

    roadind = move_along(roadind_0, roadway, 1.0)
    @test roadind == RoadIndex(CurveIndex(2,0), LaneTag(2,1))

    roadind = move_along(roadind_0, roadway, 1.25)
    @test roadind == RoadIndex(CurveIndex(2,0.25), LaneTag(2,1))

    roadind = move_along(roadind_0, roadway, 3.0)
    @test roadind == RoadIndex(CurveIndex(3,1.0), LaneTag(2,1))

    roadind = move_along(roadind_0, roadway, 4.0)
    @test roadind == RoadIndex(CurveIndex(1,0.0), LaneTag(3,1))

    roadind = move_along(roadind_0, roadway, 4.5)
    @test roadind == RoadIndex(CurveIndex(1,0.5), LaneTag(3,1))

    roadind = move_along(roadind_0, roadway, 3.75)
    @test roadind == RoadIndex(CurveIndex(0,0.75), LaneTag(3,1))

    roadind = move_along(roadind_0, roadway, -1.0)
    @test roadind == RoadIndex(CurveIndex(0,0.0), LaneTag(2,1))

    roadind = move_along(roadind_0, roadway, -1.75)
    @test roadind == RoadIndex(CurveIndex(1,0.25), LaneTag(1,1))

    roadind = move_along(roadind_0, roadway, -0.75)
    @test roadind == RoadIndex(CurveIndex(0,0.25), LaneTag(2,1))

    roadind = move_along(roadind_0, roadway, -Inf)
    @test roadind == RoadIndex(CurveIndex(1,0.0), LaneTag(1,1))

    roadind = move_along(roadind_0, roadway, Inf)
    @test roadind == RoadIndex(CurveIndex(1,1.0), LaneTag(3,1))

    ####

    @test n_lanes_right(roadway[LaneTag(2,1)], roadway) == 0
    @test n_lanes_right(roadway[LaneTag(2,2)], roadway) == 1
    @test n_lanes_right(roadway[LaneTag(2,3)], roadway) == 2
    @test n_lanes_left(roadway[LaneTag(2,1)], roadway) == 2
    @test n_lanes_left(roadway[LaneTag(2,2)], roadway) == 1
    @test n_lanes_left(roadway[LaneTag(2,3)], roadway) == 0

end

let
    roadway = Roadway()

    seg1 = RoadSegment(1,
        [Lane(LaneTag(1,1),
            [CurvePt(VecSE2(0.0,0.0,0.0), 0.0),
             CurvePt(VecSE2(1.0,0.0,0.0), 1.0)],
             width=1.0,
             boundary_left = LaneBoundary(:solid, :white),
             boundary_right = LaneBoundary(:broken, :yellow)),
         Lane(LaneTag(1,2),
            [CurvePt(VecSE2(0.0,1.0,0.0), 0.0),
             CurvePt(VecSE2(1.0,1.0,0.0), 1.0)],
             width=2.0,
             boundary_left = LaneBoundary(:double, :white))],
         )

    seg2 = RoadSegment(2,
        [Lane(LaneTag(2,1),
            [CurvePt(VecSE2(4.0,0.0,0.0), 0.0),
             CurvePt(VecSE2(5.0,0.0,0.0), 1.0)]),
         Lane(LaneTag(2,2),
            [CurvePt(VecSE2(4.0,1.0,0.0), 0.0),
             CurvePt(VecSE2(5.0,1.0,0.0), 1.0)])],
         )

    connect!(seg1.lanes[1], seg2.lanes[2])

    push!(roadway.segments, seg1)
    push!(roadway.segments, seg2)

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

    path, io = mktemp()
    write(io, MIME"text/plain"(), roadway)
    close(io)

    lines = open(readlines, path)

    for (line_orig, line_test) in zip(lines,
            [
             "ROADWAY",
             "2",
             "1",
             "   2",
             "   1",
             "       1.000",
             "       -Inf Inf",
             "       solid white",
             "       broken yellow",
             "       1",
             "           D (1 1.000000) 1 0.000000 2 2",
             "       2",
             "           (0.0000 0.0000 0.000000) 0.0000 NaN NaN",
             "           (1.0000 0.0000 0.000000) 1.0000 NaN NaN",
             "   2",
             "       2.000",
             "       -Inf Inf",
             "       double white",
             "       unknown unknown",
             "       0",
             "       2",
             "           (0.0000 1.0000 0.000000) 0.0000 NaN NaN",
             "           (1.0000 1.0000 0.000000) 1.0000 NaN NaN",
             "2",
             "   2",
             "   1",
             "       3.000",
             "       -Inf Inf",
             "       unknown unknown",
             "       unknown unknown",
             "       0",
             "       2",
             "           (4.0000 0.0000 0.000000) 0.0000 NaN NaN",
             "           (5.0000 0.0000 0.000000) 1.0000 NaN NaN",
             "   2",
             "       3.000",
             "       -Inf Inf",
             "       unknown unknown",
             "       unknown unknown",
             "       1",
             "           U (1 0.000000) 1 1.000000 1 1",
             "       2",
             "           (4.0000 1.0000 0.000000) 0.0000 NaN NaN",
             "           (5.0000 1.0000 0.000000) 1.0000 NaN NaN",
             ]
        )
        @test strip(line_orig) == strip(line_test)
    end

    io = open(path)
    roadway2 = read(io, MIME"text/plain"(), Roadway)
    close(io)
    rm(path)

    @test length(roadway.segments) == length(roadway2.segments)
    for (seg1, seg2) in zip(roadway.segments, roadway2.segments)
        @test seg1.id == seg2.id
        @test length(seg1.lanes) == length(seg2.lanes)
        for (lane1, lane2) in zip(seg1.lanes, seg2.lanes)
            @test lane1.tag == lane2.tag
            @test isapprox(lane1.width, lane2.width, atol=1e-3)
            @test isapprox(lane1.speed_limit.lo, lane2.speed_limit.lo, atol=1e-3)
            @test isapprox(lane1.speed_limit.hi, lane2.speed_limit.hi, atol=1e-3)
            @test lane1.boundary_left == lane2.boundary_left
            @test lane1.boundary_right == lane2.boundary_right
            for (conn1, conn2) in zip(lane1.exits, lane2.exits)
                @test conn1.downstream == conn2.downstream
                @test conn1.mylane.i == conn2.mylane.i
                @test isapprox(conn1.mylane.t, conn2.mylane.t, atol=1e-3)
                @test conn1.target.tag == conn2.target.tag
                @test conn1.target.ind.i == conn2.target.ind.i
                @test isapprox(conn1.target.ind.t, conn2.target.ind.t, atol=1e-3)
            end
            for (conn1, conn2) in zip(lane1.entrances, lane2.entrances)
                @test conn1.downstream == conn2.downstream
                @test conn1.mylane.i == conn2.mylane.i
                @test isapprox(conn1.mylane.t, conn2.mylane.t, atol=1e-3)
                @test conn1.target.tag == conn2.target.tag
                @test conn1.target.ind.i == conn2.target.ind.i
                @test isapprox(conn1.target.ind.t, conn2.target.ind.t, atol=1e-3)
            end
            for (pt1, pt2) in zip(lane1.curve, lane2.curve)
                @test abs2(convert(VecE2, pt1.pos) - convert(VecE2, pt2.pos)) < 0.01
                @test angledist(pt1.pos.θ, pt2.pos.θ) < 1e-5
                @test isnan(pt1.s) || isapprox(pt1.s, pt2.s, atol=1e-3)
                @test isnan(pt1.k) || isapprox(pt1.k, pt2.k, atol=1e-6)
                @test isnan(pt1.kd) || isapprox(pt1.kd, pt2.kd, atol=1e-6)
            end
        end
    end
end

let
    lc = LaneConnection(true, CurveIndex(1,0.0), RoadIndex(CurveIndex(2,0.0), LaneTag(2,2)))
    show(IOBuffer(), lc)
end