text:
using FastTransforms, Compat
using Compat.Test

import FastTransforms: normalizecolumns!, maxcolnorm

@testset "Fast plan" begin
    m, n = 362, 300

    A = sphrandn(Float64, m, n);
    normalizecolumns!(A);
    B = zero(A);
    C = zero(A);
    D = zero(A);
    E = zero(A);

    @time FP = FastSphericalHarmonicPlan(A; sketch = :none);
    @time SP = SlowSphericalHarmonicPlan(A);

    @time mul!(B, SP, A);
    @time mul!(C, FP, A);
    @time At_mul_B!(D, SP, B);
    @time At_mul_B!(E, FP, C);

    @test maxcolnorm(B-C) < 10eps()*max(m,n)
    @test maxcolnorm(A-D) < 10eps()*max(m,n)
    @test maxcolnorm(A-E) < 10eps()*max(m,n)
    @test maxcolnorm(D-E) < 10eps()*max(m,n)

    println("The forward difference between slow and fast plans: ", maxcolnorm(B-C))
    println("The backward difference between slow plan and original: ", maxcolnorm(A-D))
    println("The backward difference between fast plan and original: ", maxcolnorm(A-E))
    println("The backward difference between slow and fast plans: ", maxcolnorm(D-E))
end


text:
using Test, Random, LinearAlgebra, SparseArrays
using BitBasis

using QuAlgorithmZoo
import QuAlgorithmZoo: _num_grover_step
using Yao

function GroverSearch(oracle, num_bit::Int; psi::DefaultRegister = uniform_state(num_bit))
    it = groveriter(psi, oracle)
    for l_psi in it psi = l_psi end
    return (it.niter, psi)
end

function inference(psi::DefaultRegister, evidense::Vector{Int}, num_iter::Int)
    oracle = inference_oracle(evidense)(nqubits(psi))
    it = groveriter(psi, oracle)
    for l_psi in it psi = l_psi end
    it.niter, psi
end

@testset "oracle" begin
    oracle = inference_oracle([2,-1,3])(3)
    # ≈ method for Identity/PermuteMultiply/Sparse
    # add and mimus between sparse matrices.
    # alway use sorted CSC format.
    v = ones(1<<3)
    v[Int(0b110)+1] *= -1
    @test mat(oracle) ≈ Diagonal(v)
    @test mat(chain(phase(π), Z)) ≈ -mat(Z)
end

@testset "Grover Search" begin
    ####### Construct Grover Search Using Reflection Block
    num_bit = 12
    oracle = inference_oracle(push!(collect(Int, 1:num_bit-1), num_bit))(num_bit)

    niter, psi = GroverSearch(oracle, 12)
    target_state = zeros(1<<num_bit); target_state[end] = 1
    @test isapprox(abs(statevec(psi)'*target_state), 1, atol=1e-3)
end

@testset "groverblock" begin
    psi = uniform_state(5)
    or = inference_oracle(5, [-1,2,5,4,3])
    func_or = or
    gb = groverblock(or, psi)
    gb2 = groverblock(func_or, psi)
    @test apply!(copy(psi), gb) == (for l_psi in groveriter(copy(psi), func_or) psi = l_psi end; psi)
    @test apply!(copy(psi), gb) == apply!(copy(psi), gb2)
end

@testset "test inference" begin
    num_bit = 12
    psi0 = rand_state(num_bit)
    #psi0 = uniform_state(num_bit)
    evidense = [1, 2, 3, 4, 5, 6, 7, 8, 9]
    #evidense = collect(1:num_bit)

    # the desired subspace
    basis = collect(UInt, 0:1<<num_bit-1)
    subinds = [itercontrol(num_bit, abs.(evidense), Int.(evidense.>0))...]

    v_desired = statevec(psi0)[subinds .+ 1]
    p = norm(v_desired)^2
    v_desired[:] ./= sqrt(p)

    # search the subspace
    num_iter = _num_grover_step(p)
    niter, psi = inference(psi0, evidense, num_iter)
    @test isapprox((psi.state[subinds .+ 1]'*v_desired) |> abs2, 1, atol=3e-2)
end

ext: jl


text:
include("generated_defs.jl")
include("helper_functions.jl")

const UA_STRING_NULL = UA_String(0, C_NULL)
const UA_GUID_NULL = UA_Guid(0, 0, 0, Tuple(zeros(UA_Byte, 8)))
const UA_NODEID_NULL = UA_NodeId(
    0,
    UA_NodeIdType(0),
    anonymous_struct_tuple(UInt32(0), fieldtype(UA_NodeId, :identifier))
)
const UA_EXPANDEDNODEID_NULL = UA_ExpandedNodeId(UA_NODEID_NULL, UA_STRING_NULL, 0)

include("types.jl")
include("server.jl")
include("client.jl")
include("init.jl")


text:
using BenchmarkTools

struct MyArray1{N, T}
   A::Array{T, N}
end

Base.length(mA::MyArray1) = length(mA.A)
Base.getindex(mA::MyArray1, I...) = getindex(mA.A, I...)
Base.setindex!(mA::MyArray1, I...) = setindex!(mA.A, I...)

struct MyArray2{N, T}
   A::Array{T, N}
end

Base.length(mA::MyArray2) = length(mA.A)

@propagate_inbounds function Base.getindex(mA::MyArray2, I...)
   @boundscheck checkbounds(A, I...)
   @inbounds mA.A[I...]
end

@propagate_inbounds function Base.setindex!(mA::MyArray2, val, I...)
   @boundscheck checkbounds(A, I...)
   mA.A[I...] = val
end

function dosomething!(A, x)
   @inbounds begin
      A[1] = 1; A[2] = x
      for n = 3:length(A)
         A[n] = 2 * x * A[n-1] - A[n-2]
      end
   end
   return A
end

A = zeros(20)
A1 = MyArray1(A)
A2 = MyArray2(A)

@btime dosomething!($A, 0.123)
@btime dosomething!($A1, 0.123)
@btime dosomething!($A2, 0.123)


text:
# Run code for previous steps
include("neid_solar_1_read.jl")
# Don't import Plots if make_plots set to false
make_plots = isdefined(Main,:make_plots) ? make_plots : true
 if make_plots
   using Plots
 end

order_list_timeseries = RvSpectML.make_order_list_timeseries(solar_data)

order_list_timeseries = RvSpectML.filter_bad_chunks(order_list_timeseries)

RvSpectML.normalize_spectra!(order_list_timeseries,solar_data);

if make_plots
   order_idx = 20:22
   plt = RvSpectML.plot_spectrum_chunks(order_list_timeseries, order_idx)
end

using DataFrames, CSV, Query
 using Statistics

vald_filename = joinpath(ancilary_solar_data_path,"VALD_Fe1_DP_rejectTelluricSlope0.0_badLineFilterESPRESSO-strict-NEID-BIS_overlapcutoff6e-05_depthcutoff0.05_allowBlends0_wavesReiners_depthssolar_nbin1depth0.mas")
 vald_df = RvSpectML.read_mask_vald(vald_filename)
 lambda_range_with_data = (min = maximum(d->minimum(d.λ),solar_data), max = minimum(d->maximum(d.λ),solar_data) )
 orig_chunk_list_df = vald_df |>
    @filter(lambda_range_with_data.min <= _.lambda ) |>
    @filter( _.lambda < lambda_range_with_data.max) |>
#    @filter( _.lambda < 6000.0 ) |>                       # Avoid tellurics at redder wavelengths
#    @filter( _.lambda >6157 || _.lambda < 6155  ) |>   # Avoid "line" w/ large variability
   DataFrame

find_overlapping_chunks(orig_chunk_list_df)
chunk_list_df = RvSpectML.merge_chunks(orig_chunk_list_df)
@assert find_overlapping_chunks(chunk_list_df) == nothing

# TODO:  FIX:  Problem with ESPRESSO mask having lots of lines getting meged into humongous chunks that don't fit into the order!
chunk_list_timeseries = RvSpectML.make_chunk_list_timeseries(solar_data,chunk_list_df)

# Check that no NaN's included
(chunk_list_timeseries, chunk_list_df) = RvSpectML.filter_bad_chunks(chunk_list_timeseries,chunk_list_df)
println(size(chunk_list_df), " vs ", num_chunks(chunk_list_timeseries) )

using Plots
if make_plots
   #order_idx = 20:22
   # plt = RvSpectML.plot_spectrum_chunks(order_list_timeseries, order_idx)
   chunk_idx = 23
   for c in chunk_idx
       t = 1
       #if(sum(chunk_list_df.line_depths[c])<0.25) continue end
       λ_mid = sqrt(chunk_list_df.lambda_hi[c]*chunk_list_df.lambda_lo[c])
       println("c= ",c , " λs= ",chunk_list_df.lambda[c]," depths= ",chunk_list_df.depth[c])
       flush(stdout)
   end
   plt = RvSpectML.plot_spectrum_chunks(chunk_list_timeseries, chunk_idx, plt=plot(), color=:black)
   display(plt)

end


text:
using ApproximateBayesianComputing
const ABC = ApproximateBayesianComputing
using Distributions
using PDMats
import LinearAlgebra.diagm
import ApproximateBayesianComputing.CompositeDistributions.CompositeDist
import Statistics
import SpecialFunctions.erf
import Compat.view # For backward compatability with v0.4

# Set HyperPrior for Population Parameters
param_prior_rate = Gamma(1.0,1.0)
 param_prior_mean  = MvLogNormal(log.([10.0, 10.0]), diagm(log.([10.0,10.0])))
 param_prior_covar_diag = MvLogNormal(zeros(2),ones(2))
 param_prior_cor = MvNormal(zeros(1),ones(1));
 param_prior = CompositeDist( ContinuousDistribution[param_prior_rate,param_prior_mean,param_prior_covar_diag,param_prior_cor] )
 num_star_default = 1000   # How many "stars" to include in universe

 function is_valid_mean_covar(theta::Array{Float64,1})  # Ensure that the proposed covariance matrix is positive definite
  a = theta[4]^2                     # WARNING: Covariance matrix extracted here and when generating data must match
  b = theta[4]*theta[5]*theta[6]
  c = theta[4]*theta[5]*theta[6]
  d = theta[5]^2
  tr_covar = a+d
  det_covar = a*d-b*c
  if (det_covar>0) && (tr_covar>0)
    return true
  else
    return false
  end
 end

function detection_prob_snr_term(snr::Real) # Detection probability as a function of signal-to-noise ratio
  sig =  1.0
  snr_thresh = 7.1
  snr_detect_never =  3.0
  snr_detect_always = 10.0
  sqrt2 = sqrt(2.0)
  pdet = 0.0
  if snr > snr_detect_always
     pdet = 1.0
  elseif snr > snr_detect_never
     pdet = 0.5*( 1+erf((snr - snr_thresh)/(sqrt2*sig)) )
  end
  return pdet
 end

 function detection_prob_geometry_term(period::Real) # Transit probability as a function of orbital geometry
    radius_star = 0.005
    days_in_year = 365.2425
    a = (period/days_in_year)^(2//3)
    prob_geometric = min(1.0,radius_star/a)
 end

 detection_prob_return_one(period::Real, snr::Real) = 1.0
 detection_prob_snr_only(period::Real, snr::Real) = detection_prob_snr_term(snr)
 detection_prob_geometry_only(period::Real, snr::Real) = detection_prob_geometry_term(period)
 detection_prob_snr_and_geometry(period::Real, snr::Real) = detection_prob_geometry_term(period) * detection_prob_snr_term(snr)

 # Code to generate simulated data given array of model parameters
 function gen_number_planets(theta::Array, num_stars::Integer)
    #rate = theta[1]>0 ? sum(rand(Poisson(theta[1]),num_stars)) : 0
    rate = theta[1]^2
    sum(rand(Poisson(rate),num_stars))
 end
 function gen_period_snr(theta::Array, n::Integer)  # Just return "planet" properties
  mean = theta[2:3]
  covar = PDMat([theta[4]^2 theta[4]*theta[5]*theta[6]; theta[4]*theta[5]*theta[6] theta[5]^2])
  period_snr = rand(MvLogNormal(mean,covar),n)
  return period_snr
 end

 function gen_period_snr_detected(theta::Array, num_stars::Integer; detection_prob::Function = detection_prob_snr_and_geometry)  # "Planets" detection efficient depends on SNR and orbital period
  n = gen_number_planets(theta,num_stars)
  period_snr = gen_period_snr(theta,n)
  detected = trues(n)
  for i in 1:n
    period = period_snr[1,i]
    snr = period_snr[2,i]
    prob_detect = detection_prob( period, snr )
    if rand() > prob_detect
      detected[i] = false
    end
  end
  return (period_snr, detected)
 end

 gen_period_snr_detected_all(theta::Array, n::Integer = num_star_default) = gen_period_snr_detected(theta,n,detection_prob=detection_prob_return_one)
 gen_period_snr_detected_snr_only(theta::Array, n::Integer = num_star_default) = gen_period_snr_detected(theta,n,detection_prob=detection_prob_snr_only)
 gen_period_snr_detected_geometry_only(theta::Array, n::Integer = num_star_default) = gen_period_snr_detected(theta,n,detection_prob=detection_prob_geometry_only)
 gen_period_snr_detected_snr_and_geo(theta::Array, n::Integer = num_star_default) = gen_period_snr_detected(theta,n,detection_prob=detection_prob_snr_and_geometry)

function calc_summary_stats_mean_stddev_detected(data::Tuple{Array{Float64,2},BitArray{1}})
  period_snr = data[1]
  detected = data[2]
  @assert size(period_snr,1) == 2
  @assert size(period_snr,2) == length(detected)
  num_detected = length(findall(detected))
  if num_detected<1
    return [0/num_star_default,Inf,Inf,Inf,Inf,Inf]
  end
  mean_log_period = mean(log.(period_snr[1,detected]))
  mean_log_snr = mean(log.(period_snr[2,detected]))
  if num_detected<2
     return [num_detected/num_star_default,mean_log_period,Inf,mean_log_snr,Inf,Inf]
  end
  stddev_log_period = Statistics.stdm(log.(period_snr[1,detected]),mean_log_period)
  stddev_log_snr = Statistics.stdm(log.(period_snr[2,detected]),mean_log_snr)
  covar_cross = sum( (log.(period_snr[1,detected]).-mean_log_period) .* (log.(period_snr[2,detected]).-mean_log_snr) )/(length(findall(detected))-1)
  cor =  covar_cross / ( stddev_log_period*stddev_log_snr*(num_detected-1) )
  #cor_F = atanh(cor)*max(1,sqrt(num_detected-3))
  #c4 = calc_c4(num_detected)
  #stddev_coeff = sqrt((num_detected-1-c4*c4)/num_detected)
  return [num_detected/num_star_default, mean_log_period, stddev_log_period, mean_log_snr, stddev_log_snr, cor]
  #[ sqrt(num_detected)/num_star_default, stddev_log_period/sqrt(num_detected), stddev_coeff*stddev_log_period, stddev_log_snr/sqrt(num_detected), stddev_coeff*stddev_log_snr, 1.0 ]
end

function est_sigma_pred(x::Array{Float64,2},y::Array{Float64,1}) #,index::Integer)
  n = num_star_default*y[1]
  sd_logp = y[3]
  sd_logsnr = y[5]
  return [sqrt(n)/num_star_default, sd_logp/sqrt(n), sd_logp, sd_logsnr/sqrt(n), sd_logsnr, 1.0]
end

function calc_c4(n)
  sqrt(2//(n-1)) * (gamma(n//2)/gamma((n-1)//2))
end

theta_true = [1.0, log(10.0), log(10.0), log(10.0), log(10.0), 0.0] # True Population Parameters
#=
theta_true = fill(NaN,5) #  Generate from hyperprior
 while !is_valid_mean_covar(theta_true)
  theta_true = rand(param_prior)
 end
 theta_true
=#
# Uncomment one of three options below one to pick how data is generated
#gen_data = gen_period_snr_detected_all
#num_star_default = 2000   # How many "planets" to include in universe
gen_data = gen_period_snr_detected_snr_only
num_star_default = 4000   # How many "planets" to include in universe
#gen_data = gen_period_snr_detected_snr_and_geo
#num_star_default = 24000   # How many "planets" to include in universe

#=
ex_data =gen_data(theta_true)
calc_summary_stats_mean_stddev_detected(ex_data)
length(find(ex_data[2]))
=#

# Tell ABC what it needs to know for a simulation
abc_plan = abc_pmc_plan_type(gen_data,calc_summary_stats_mean_stddev_detected,ABC.calc_scaled_dist_max, param_prior;
                             is_valid=is_valid_mean_covar,adaptive_quantiles=true,num_max_attempt=1000,
                             target_epsilon=0.01,epsilon_init=10.0,num_part=50,
   make_proposal_dist = ABC.make_proposal_dist_gaussian_full_covar);

# Generate "true/observed data"
data_true = gen_data(theta_true)   # Draw "real" data from same model as for analysis
ss_true = abc_plan.calc_summary_stats(data_true)
ABC.set_distance_scale(abc_plan, theta_true, ss = ss_true )

#data_star =
#ss_star =
#dist_star = plan.calc_dist(ss_true,plan.calc_summary_stats(plan.gen_data(theta_star)))

#=
pdf(abc_plan.prior,theta_true)
CompositeDistributions.index(abc_plan.prior,1)
CompositeDistributions.index(abc_plan.prior,2)
logpdf(abc_plan.prior.dist[1],theta_true[CompositeDistributions.index(abc_plan.prior,1)])
logpdf(abc_plan.prior.dist[2],theta_true[CompositeDistributions.index(abc_plan.prior,2)])

logpdf(abc_plan.prior.dist[1],theta_true[1:2])
logpdf(abc_plan.prior.dist[2],theta_true[3:5])
println("***************************************************")
=#
# Run ABC simulation
@time pop_out = run_abc(abc_plan,ss_true, verbose=true);

if (Base.find_package("JLD2") != nothing) && (Base.find_package("FileIO") != nothing)

using JLD2, FileIO
save("ex4_out.jld2", "pop_out", pop_out)

end

if Base.find_package("PyPlot") != nothing

using PyPlot
function plot_abc_posterior(theta_plot::Array{Float64,2}, index::Integer)
  lo = minimum(theta_plot[index,:])
  hi = maximum(theta_plot[index,:])
  sigma = 1.0*sqrt(abs(cov(vec(theta_plot[index,:]))))
  x = collect(linspace(lo-3*sigma,hi+3*sigma,100))
  y = zeros(x)
  for i in 1:length(x)
    for j in 1:size(theta_plot,2)
      y[i] += pdf(Normal(theta_plot[index,j],sigma),x[i])/size(theta_plot,2)
    end
  end
  plot(x,y,"-")
end

plot_abc_posterior(pop_out.theta,1)
plot_abc_posterior(pop_out.theta,2)
plot_abc_posterior(pop_out.theta,3)
plot_abc_posterior(pop_out.theta,4)
plot_abc_posterior(pop_out.theta,5)
plot_abc_posterior(pop_out.theta,6)

#=
sampler_plot = abc_plan.make_proposal_dist(pop_out, 1.0) # abc_plan.tau_factor)
num_plot = 100
epsilon = 0.05
theta_plot, dist_plot = generate_abc_sample(abc_plan,sampler_plot, ss_true, epsilon, n=num_plot )
plot_abc_posterior(theta_plot,1)
med_dist = median(dist_plot)
plot_abc_posterior(theta_plot[:,find(dist_plot.<=med_dist)],1)
=#

end


text:
using Dates

## --------------------------------struct Definition------------------------------------------

mutable struct Header
    fileVersion       ::Int64
    dataType1ch       ::String     # 1chのデータ種別
    dataType2ch       ::String     # 2chのデータ種別
    fileHeaderSize    ::Int64      # ファイルヘッダサイズ（コメントのサイズ含む。）
    frameHeaderSize   ::Int64      # フレームヘッダサイズ（コメントのサイズ含む。）
    operatorNameSize  ::Int64
    commentOffsetSize ::Int64
    commentSize       ::Int64
    pixelX            ::Int64
    pixelY            ::Int64
    scanningRangeX    ::Int64
    scanningRangeY    ::Int64
    frameRate         ::Float64
    piezoExtensionX   ::Float64    # Xピエゾ伸び係数[nm/V]
    piezoExtensionY   ::Float64    # Yピエゾ伸び係数[nm/V]
    piezoExtensionZ   ::Float64    # Zピエゾ伸び係数[nm/V]
    piezoGainZ        ::Float64
    adRange           ::Float64    # AD電圧レンジ
    AdResolution      ::Float64    # AD分解能
    isAveraged        ::Bool       # 移動平均化フラグ（tureで移動平均）
    averageWindow     ::Int64      # 1ピクセルに使える最大データ数(最低１)
    day               ::DateTime
    roundingDegreeX   ::Int64
    roundingDegreeY   ::Int64
    maxRangeX         ::Float64
    maxRangeY         ::Float64
    booked1           ::Int64
    booked2           ::Int64
    booked3           ::Int64
    initFrame         ::Int64
    numFrames         ::Int64
    machineId         ::Int64
    fileId            ::Int64
    operatorName      ::Array{Int, 1}
    sensorSensitivity ::Float64       # センサー感度
    phaseSensitivity  ::Float64       # 位相感度
    scannigDirection  ::Int64         # データ取得ステータスコード
    comment           ::Array{Int, 1}
    textEncodingCodepage ::Int64

    Header() = new()
end

struct FrameHeader
    number       ::Int64
    maxData      ::Int64
    minData      ::Int64
    offsetX      ::Int64
    offsetY      ::Int64
    tiltX        ::Float64
    tiltY        ::Float64
    isStimulated ::Bool
    booked1      ::Int64
    booked2      ::Int64
    booked3      ::Int64
    booked4      ::Int64
end

struct Frame
    header   ::FrameHeader
    data     ::Array{Float64, 2}
    subData  ::Union{Array{Float64, 2}, Nothing}
end

# 現在単位はÅ
struct Asd
    header   ::Header
    frames   ::Array{Frame, 1}
end

## --------------------------------struct Definition------------------------------------------
## -----------------------------------read helper---------------------------------------------

function getDataType(data)
    if     data == 0x5054 return "topography"
    elseif data == 0x5245 return "error"
    elseif data == 0x4850 return "phase"
    else return "none" end
end

function getAdRange(data)
    # ADレンジ定数（ユニポーラ0～1V）(使われていないらしい)
    if     data == 0x00000001 @assert false "invalid AdRange"
    # ADレンジ定数（ユニポーラ0～2.5V）(使われていないらしい)
    elseif data == 0x00000002 @assert false "invalid AdRange"
    # ADレンジ定数（ユニポーラ0～5V）(使われていないらしい)
    elseif data == 0x00000004 @assert false "invalid AdRange"
    # ADレンジ定数（バイポーラ±1V）
    elseif data == 0x00010000 return 2.0
    # ADレンジ定数（バイポーラ±2.5V）
    elseif data == 0x00020000 return 5.0
    # ADレンジ定数（バイポーラ±5V）
    elseif data == 0x00040000 return 10.0
    # ADレンジ定数（バイポーラ±80V, データを編集した場合に仮想的にこれを使う。実際にバイポーラ80VでAD変換したわけではない。また、分解能は16ビットにする）
    elseif data == 0x00080000 return 160.0
    # 何も当てはまらない
        else @assert false "invalid AdRange" end
    return nothing
end

function getOperatorName(io, operatorNameSize)
    ret = []
    for i in 1:operatorNameSize
        push!(ret, Int64(read(io, UInt8)))
    end
    return ret
end

function getComment(io, commentOffsetSize, commentSize)
    for i in 1:commentOffsetSize
        read(io, Bool)
    end
    ret = []
    for i in 1:commentSize
        push!(ret, Int64(read(io, UInt8)))
    end
    return ret
end

## -----------------------------------read helper---------------------------------------------

function readDateV0(io::IOStream)
    year              = Int64(read(io, Int16))
    month             = Int64(read(io, UInt8))
    day               = Int64(read(io, UInt8))
    hour              = Int64(read(io, UInt8))
    minute            = Int64(read(io, UInt8))
    second            = Int64(read(io, UInt8))
    return DateTime(year, month, day, hour, minute, second)
end

function readDateV1(io::IOStream)
    year              = Int64(read(io, Int32))
    month             = Int64(read(io, Int32))
    day               = Int64(read(io, Int32))
    hour              = Int64(read(io, Int32))
    minute            = Int64(read(io, Int32))
    second            = Int64(read(io, Int32))
    return DateTime(year, month, day, hour, minute, second)
end

function readHeaderV0(io::IOStream)
    header = Header()
    header.fileVersion       = 0
    header.dataType1ch       = getDataType(read(io, Int16))
    header.dataType2ch       = getDataType(read(io, Int16))
    header.fileHeaderSize    = Int64(read(io, Int32))
    header.frameHeaderSize   = Int64(read(io, Int32))
    header.operatorNameSize  = Int64(read(io, Int32))
    header.commentOffsetSize = Int64(read(io, Int32))
    header.commentSize       = Int64(read(io, Int32))
    header.pixelX            = Int64(read(io, Int16))
    header.pixelY            = Int64(read(io, Int16))
    header.scanningRangeX    = Int64(read(io, Int16))
    header.scanningRangeY    = Int64(read(io, Int16))
    header.frameRate         = Float64(read(io, Float32))
    header.piezoExtensionZ   = Float64(read(io, Float32))
    header.piezoGainZ        = Float64(read(io, Float32))
    header.adRange           = getAdRange(read(io, UInt32))
    header.AdResolution      = (2.0)^Int64(read(io, Int32))
    header.isAveraged        = read(io, Bool)
    header.averageWindow     = Int64(read(io, Int16))
    legacy                   = Int64(read(io, Int32)) # ダミー
    header.day               = readDateV0(io)
    header.roundingDegreeX   = Int64(read(io, UInt8))
    header.maxRangeX         = Float64(read(io, Float32))
    header.maxRangeY         = Float64(read(io, Float32))
    header.booked1           = Int64(read(io, Int32))
    header.booked2           = Int64(read(io, Int32))
    header.booked3           = Int64(read(io, Int32))
    header.initFrame         = Int64(read(io, Int32))
    header.numFrames         = Int64(read(io, Int32))
    header.machineId         = Int64(read(io, Int32))
    header.fileId            = Int64(read(io, Int16))
    header.operatorName      = getOperatorName(io, header.operatorNameSize)
    header.sensorSensitivity = Float64(read(io, Float32))
    header.phaseSensitivity  = Float64(read(io, Float32))
    header.scannigDirection  = Int64(read(io, Int32))
    header.comment           = getComment(io, header.commentOffsetSize, header.commentSize)

    return header
end

function readHeaderV1(io::IOStream)
    header = Header()
    header.fileVersion       = 1
    header.fileHeaderSize    = Int64(read(io, Int32))
    header.frameHeaderSize   = Int64(read(io, Int32))
    header.textEncodingCodepage = Int64(read(io, Int32))
    header.operatorNameSize  = Int64(read(io, Int32))
    header.commentSize       = Int64(read(io, Int32))
    header.dataType1ch       = getDataType(read(io, Int32))
    header.dataType2ch       = getDataType(read(io, Int32))
    header.initFrame         = Int64(read(io, Int32))
    header.numFrames         = Int64(read(io, Int32))
    header.scannigDirection  = Int64(read(io, Int32))
    header.fileId            = Int64(read(io, Int32))
    header.pixelX            = Int64(read(io, Int32))
    header.pixelY            = Int64(read(io, Int32))
    header.scanningRangeX    = Int64(read(io, Int32))
    header.scanningRangeY    = Int64(read(io, Int32))
    header.isAveraged        = read(io, Bool)
    header.averageWindow     = Int64(read(io, Int32))
    header.day               = readDateV1(io)
    header.roundingDegreeX   = Int64(read(io, Int32))
    header.roundingDegreeY   = Int64(read(io, Int32))
    header.frameRate         = Float64(read(io, Float32))
    header.sensorSensitivity = Float64(read(io, Float32))
    header.phaseSensitivity  = Float64(read(io, Float32))
    header.commentOffsetSize = Int64(read(io, Int32))
    header.booked1           = Int64(read(io, Int32))
    header.booked2           = Int64(read(io, Int32))
    header.booked3           = Int64(read(io, Int32))
    header.machineId         = Int64(read(io, Int32))
    header.adRange           = getAdRange(read(io, Int32))
    header.AdResolution      = (2.0)^Int64(read(io, Int32))
    header.maxRangeX         = Float64(read(io, Float32))
    header.maxRangeY         = Float64(read(io, Float32))
    header.piezoExtensionX   = Float64(read(io, Float32))
    header.piezoExtensionY   = Float64(read(io, Float32))
    header.piezoExtensionZ   = Float64(read(io, Float32))
    header.piezoGainZ        = Float64(read(io, Float32))
    header.operatorName      = getOperatorName(io, header.operatorNameSize)
    header.comment           = getComment(io, 0, header.commentSize)

    return header
end

"""
ADのバイナリ－物理量の変換公式は
物理量 = ボードのAD変換レンジ*(サンプリングバイナリデータ)/2^(分解能) - ボードのAD変換レンジの半分。
また、PIDの信号を高さ情報として取り込んでいるので、

バイナリデータが大きい　＝　PIDの出力電圧が高い　＝　試料に対して押し込んでいる　＝　高さが低い

という関係になる。
したがって、高さの最大最小は

高さ最値小 = バイナリ最大値から計算した高さ
高さ最大値 = バイナリ最小値から計算した高さ
"""
function binaryToPhysicalQuantity(data, header, chanelType)
    # 電圧データに変換
    cc = header.adRange / header.AdResolution
    adUiniRange = header.adRange / 2

    # 電圧データを高さor位相データに変換するための乗数。ErrorとPhaseの場合はTopo像とは逆符号になる点に注意。
    multiplier = Float64(0)
    if chanelType == "topography"
        multiplier = header.piezoGainZ * header.piezoExtensionZ
    elseif chanelType == "error"
        multiplier = -1.0 * header.sensorSensitivity
    elseif chanelType == "phase"
        multiplier = if header.phaseSensitivity != 0 phaseSensitivity else -1.0 end
    else
        # ここには来ない
        @assert false "invalid chanelType"
    end

    # nm -> angstrom(要議論)
    unitConversion = 10

    for y in 1:header.pixelY, x in 1:header.pixelX
        data[y, x] = (adUiniRange - data[y, x] * cc) * multiplier * unitConversion
    end
end

function readFrameHeader(io::IOStream)
    number       = Int64(read(io, Int32))
    maxData      = Int64(read(io, Int16))
    minData      = Int64(read(io, Int16))
    offsetX      = Int64(read(io, Int16))
    offsetY      = Int64(read(io, Int16))
    tiltX        = Float64(read(io, Float32))
    tiltY        = Float64(read(io, Float32))
    isStimulated = read(io, Bool)
    booked1      = Int64(read(io, Int8))
    booked2      = Int64(read(io, Int16))
    booked3      = Int64(read(io, Int32))
    booked4      = Int64(read(io, Int32))

    return FrameHeader(number, maxData, minData, offsetX, offsetY, tiltX, tiltY, isStimulated, booked1, booked2, booked3, booked4)
end

function readImage(io::IOStream, header, chanelType)
    data = zeros(header.pixelY, header.pixelX)
    # TODO: 平均化回数が1じゃないことがあるらしい(今は不必要？)
    for y in 1:header.pixelY, x in 1:header.pixelX
        data[y, x] = Int64(read(io, Int16))
    end
    binaryToPhysicalQuantity(data, header, chanelType)
    return data
end

function readFrame(io::IOStream, header)
    frameHeader = readFrameHeader(io)

    data = readImage(io, header, header.dataType1ch)

    if header.dataType2ch == "none"
        return Frame(frameHeader, data, nothing)
    end

    subData = readImage(io, header, header.dataType2ch)

    return Frame(frameHeader, data, subData)
end

function readasd(filePath)
    open(filePath, "r") do io
        headerVersion = Int64(read(io, Int32))
        header = Header()
        if headerVersion == 0
            header = readHeaderV0(io)
        elseif headerVersion == 1
            header = readHeaderV1(io)
        else
            @assert false "can't read v2 file"
        end

        frames = []
        for i in 1:header.numFrames
            push!(frames, readFrame(io, header))
        end

        return Asd(header, frames)
    end
end


text:
immutable AxFile
  version::Uint8
  subsample::Uint8
  fs::Uint32
  nfft::Uint32
  nw::Uint16
  k::Uint16
  pval::Float64
  df::Float64
  data::Array{Float64,2}
end

function ==(x::AxFile, y::AxFile)
  if (x.version != y.version)  return false;  end
  if (x.subsample != y.subsample) return false; end
  if (x.fs != y.fs) return false; end
  if (x.nfft != y.nfft) return false; end
  if (x.nw != y.nw) return false; end
  if (x.k != y.k) return false; end
  if (x.pval != y.pval) return false; end
  if (x.df != y.df) return false; end
  if (x.data != y.data) return false; end
  return true
end

function axread(filename)
  fid=open(filename,"r")

  VERSION=read(fid,Uint8)
  SUBSAMPLE=read(fid,Uint8)
  read(fid,Uint8)
  FS=read(fid,Uint32)
  NFFT=read(fid,Uint32)
  NW=read(fid,Uint16)
  K=read(fid,Uint16)
  PVAL=read(fid,Float64)
  DF=read(fid,Float64)

  first=position(fid)
  seekend(fid)
  skip(fid,-1)
  last=position(fid)
  if read(fid,Uint8)!='Z'  error("bad .ax file");  end
  seek(fid,first)
  data=transpose(read(fid, Float64, (4,int((last-first)/8/4))))

  close(fid)

  AxFile(VERSION,SUBSAMPLE,FS,NFFT,NW,K,PVAL,DF,data)
end

max_stars_repo_name: JaneliaSciComp/Ax

max_issues_repo_name: JaneliaSciComp/Ax

max_forks_repo_name: JaneliaSciComp/Ax


text:
t1 = Time(1,0)
t2 = Time(0, 999_999_999)
t3 = Time(2, 500_000_000)
d1 = Duration(0, 999_999_999)
d2 = Duration(1, 500_000_000)
d3 = Duration(0, 1)

@test t1 == Time(1,0)
@test t1 != t2
@test t1 > t2
@test t1 >= t2
@test t2 < t1
@test t2 <= t1

@test d1 == Duration(0, 999_999_999)
@test d1 != d2
@test d1 < d2
@test d1 <= d2
@test d2 > d1
@test d2 >= d1

@test t1 + d2 == t3
@test t2 + d3 == t1
@test t1 - t2 == d3
@test t1 - d3 == t2
@test d1 + d2 + d3 == Duration(t3.secs, t3.nsecs)
@test d2 - d1 - d3 == Duration(0, 500_000_000)

@test d2*2 == Duration(3,0)
@test 3.0*d2 == Duration(4,500_000_000)

tt = Time(2,0)
@test tt == Time(2.0)
@test convert(Float64,tt) == 2.0
@test to_sec(tt) == 2.0
@test to_nsec(tt) == 2_000_000_000

dt = Duration(3,0)
@test dt == Duration(3.0)
@test convert(Float64,dt) == 3.0
@test to_sec(dt) == 3.0
@test to_nsec(dt) == 3_000_000_000

@test dt + tt == Time(5.0)
@test dt + dt == Duration(6.0)

#PyObject stuff
ptt = convert(PyCall.PyObject, tt)
@test ptt[:secs] == 2
@test ptt[:nsecs] == 0
ptt[:nsecs] = 101
tt2 = convert(Time, ptt)
@test to_nsec(tt2) == 2_000_000_101

pdt = convert(PyCall.PyObject, dt)
@test pdt[:secs] == 3
@test pdt[:nsecs] == 0
pdt[:nsecs] = 202
dt2 = convert(Duration, pdt)
@test to_nsec(dt2) == 3_000_000_202

#rostime and sleeping
t1 = get_rostime()
rossleep(0.5)
t2 = get_rostime()
@test t2 - t1 >= Duration(0.4)

rte = Rate(Duration(0.5))
rossleep(rte)
t1 = RobotOS.now()
rossleep(rte)
t2 = RobotOS.now()
rossleep(rte)
t3 = RobotOS.now()
@test t2 - t1 >= Duration(0.4)
@test t3 - t2 >= Duration(0.4)
@test t3 - t1 >= Duration(0.8)

t1 = get_rostime()
RobotOS.sleep(Duration(0.5))
t2 = get_rostime()
@test t2 - t1 >= Duration(0.4)

RobotOS.sleep(rte)
t1 = get_rostime()
RobotOS.sleep(rte)
t2 = get_rostime()
@test t2 - t1 >= Duration(0.4)


text:
# # Recipe Ingredients Example
# Following example demonstrates prediction of cuisine from set of ingredients.

# ## A gentle introduction to creation of neural networks reflexing structure of JSON documents

# This notebook serves as an introduction to `Mill` and `JsonGrinder` libraries.
# The former provides support for Multi-instance learning problems, their cascades, and their Cartesian product ([see the paper](https://arxiv.org/abs/2105.09107) for theoretical explanation).
# The latter `JsonGrinder` simplifies processing of JSON documents. It allows to infer schema of JSON documents from which it suggests an extractor to convert JSON document to a `Mill` structure.
# `JsonGrinder` defines basic set of "extractors" converting values of keys to numeric representation (matrices) or to convert them to corresponding structures in `Mill`. Naturally, this set of extractors can be extended.
#
# Below, the intended workflow is demonstrated on a simple problem of guessing type of a cuisine from a list of ingrediences.
# The whole dataset and problem description can be found [on this Kaggle page](https://www.kaggle.com/kaggle/recipe-ingredients-dataset/home).
# Note that the goal is not to achieve state of the art, but to demonstrate the workflow.

#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`recipes.ipynb`](@__NBVIEWER_ROOT_URL__/examples/recipes.ipynb)

# todo: dodělat, domigrovat z ipynb a ten pak smazat
# **Caution**
# To decrease the computational load, we keep locally in the repo only a subset of the whole dataset (39774), size of the minibatch, and size of the validation data.
# Of course these numbers are useless in practice, and therefore the resulting accuracy is poor.
# Using all samples (666920), setting minibatch size to 100, and leaving 1000 samples for validation / testing gives you accuracy 0.74 on validation data.
#
#nb # We start by installing JsonGrinder and few other packages we need for the example.
#nb # Julia Ecosystem follows philosophy of many small single-purpose composable packages
#nb # which may be different from e.g. python where we usually use fewer larger packages.
#nb using Pkg
#nb pkg"add JsonGrinder#master Flux Mill#master MLDataPattern Statistics ChainRulesCore JSON"

using JsonGrinder, Flux, Mill, MLDataPattern, Statistics, ChainRulesCore, JSON

# start by loading all samples
#src magic for resolving paths
data_file = "data/recipes.json" #src
data_file = "../../../data/recipes.json" #nb
data_file = "data/recipes.json" #md
data_file = "data/recipes.json" #jl
samples = open(data_file,"r") do fid
	Vector{Dict}(JSON.parse(read(fid, String)))
end
JSON.print(samples[1],2)

# create schema of the JSON
sch = JsonGrinder.schema(samples)

# create extractor and split it into one for loading targets and
# one for loading data, using custom function to set conditions for using n-gram representation
delete!(sch.childs,:id)

extractor = suggestextractor(sch)
extract_data = ExtractDict(deepcopy(extractor.dict))
extract_target = ExtractDict(deepcopy(extractor.dict))
delete!(extract_target.dict, :ingredients)
delete!(extract_data.dict, :cuisine)

extract_data(samples[1])
extract_target(samples[1])[:cuisine]
# we convert JSONs to Datasets
# advised to use all the samples, this is just speedup to demonstrate functionality
data = extract_data.(samples[1:5_000])
data = reduce(catobs, data)
target = extract_target.(samples[1:5_000])
target = reduce(catobs, target)[:cuisine].data

# # Create the model
m = reflectinmodel(sch, extract_data,
	layer -> Dense(layer,20,relu),
	bag -> SegmentedMeanMax(bag),
	fsm = Dict("" => layer -> Dense(layer, size(target, 1))),
)

@non_differentiable getobs(x::DataSubset{<:ProductNode})

# # Train the model
opt = Flux.Optimise.ADAM()
loss(x, y) = Flux.logitcrossentropy(m(x).data, y)
loss(x::DataSubset, y) = loss(getobs(x), y)
loss(xy::Tuple) = loss(xy...)
valdata = data[1:1_000],target[:,1:1_000]
data, target = data[1_001:5_000], target[:,1001:5_000]
# for less recourse-chungry training, we use only part of data for trainng, but it is advised to used all, as i following line:
# data, target = data[1001:nobs(data)], target[:,1001:size(target,2)]

cb = () -> println("accuracy = ",mean(Flux.onecold(m(valdata[1]).data) .== Flux.onecold(valdata[2])))
ps = Flux.params(m)
mean(Flux.onecold(m(data).data) .== Flux.onecold(target))

iterations = 20
minibatchsize = 4000
minibatches = RandomBatches((data, target), size = minibatchsize, count = iterations)

@info "testing the gradient"
loss(first(minibatches))
gs = gradient(() -> loss(first(minibatches)), ps)
Flux.Optimise.update!(opt, ps, gs)

# feel free to train for longer period of time, this example is learns only 20 iterations, so it runs fast
loss(first(minibatches))
gs = gradient(() -> loss(first(minibatches)), ps)
Flux.Optimise.train!(loss, ps, minibatches, opt, cb = Flux.throttle(cb, 2))

# calculate the accuracy
mean(Flux.onecold(m(data).data) .== Flux.onecold(target))


text:
"""
    linsolve(A::AbstractMatrix, b::AbstractVector, [x₀, a₀::Number = 0, a₁::Number = 1];
                kwargs...)
    linsolve(f, b, [x₀, a₀::Number = 0, a₁::Number = 1]; kwargs...)
    # expert version:
    linsolve(f, b, x₀, algorithm, [a₀::Number = 0, a₁::Number = 1])

Compute a solution `x` to the linear system `(a₀ + a₁ * A)*x = b` or
`a₀ * x + a₁ * f(x) = b`, possibly using a starting guess `x₀`. Return the approximate
solution `x` and a `ConvergenceInfo` structure.

### Arguments:

The linear map can be an `AbstractMatrix` (dense or sparse) or a general function or
callable object. If no initial guess is specified, it is chosen as `(zero(a₀)*zero(a₁))*b`
which should generate an object similar to `b` but initialized with zeros. The numbers `a₀`
and `a₁` are optional arguments; they are applied implicitly, i.e. they do not contribute
the computation time of applying the linear map or to the number of operations on vectors of
type `x` and `b`.

### Return values:

The return value is always of the form `x, info = linsolve(...)` with

  - `x`: the approximate solution to the problem, similar type as the right hand side `b`
    but possibly with a different `eltype`

  - `info`: an object of type [`ConvergenceInfo`], which has the following fields

      + `info.converged::Int`: takes value 0 or 1 depending on whether the solution was
        converged up to the requested tolerance
      + `info.residual`: residual `b - f(x)` of the approximate solution `x`
      + `info.normres::Real`: norm of the residual, i.e. `norm(info.residual)`
      + `info.numops::Int`: total number of times that the linear map was applied, i.e. the
        number of times that `f` was called, or a vector was multiplied with `A`
      + `info.numiter::Int`: number of times the Krylov subspace was restarted (see below)

!!! warning "Check for convergence"

    No warning is printed if no converged solution was found, so always check if
    `info.converged == 1`.

### Keyword arguments:

Keyword arguments are given by:

  - `verbosity::Int = 0`: verbosity level, i.e. 0 (no messages), 1 (single message
    at the end), 2 (information after every iteration), 3 (information per Krylov step)
  - `atol::Real`: the requested accuracy, i.e. absolute tolerance, on the norm of the
    residual.
  - `rtol::Real`: the requested accuracy on the norm of the residual, relative to the norm
    of the right hand side `b`.
  - `tol::Real`: the requested accuracy on the norm of the residual that is actually used by
    the algorithm; it defaults to `max(atol, rtol*norm(b))`. So either use `atol` and `rtol`
    or directly use `tol` (in which case the value of `atol` and `rtol` will be ignored).
  - `krylovdim::Integer`: the maximum dimension of the Krylov subspace that will be
    constructed.
  - `maxiter::Integer`: the number of times the Krylov subspace can be rebuilt; see below for
    further details on the algorithms.
  - `orth::Orthogonalizer`: the orthogonalization method to be used, see
    [`Orthogonalizer`](@ref)
  - `issymmetric::Bool`: if the linear map is symmetric, only meaningful if `T<:Real`
  - `ishermitian::Bool`: if the linear map is hermitian
  - `isposdef::Bool`: if the linear map is positive definite

The default values are given by `atol = KrylovDefaults.tol`, `rtol = KrylovDefaults.tol`,
`tol = max(atol, rtol*norm(b))`, `krylovdim = KrylovDefaults.krylovdim`,
`maxiter = KrylovDefaults.maxiter`, `orth = KrylovDefaults.orth`;
see [`KrylovDefaults`](@ref) for details.

The default value for the last three parameters depends on the method. If an
`AbstractMatrix` is used, `issymmetric`, `ishermitian` and `isposdef` are checked for that
matrix, ortherwise the default values are `issymmetric = false`,
`ishermitian = T <: Real && issymmetric` and `isposdef = false`.

### Algorithms

The final (expert) method, without default values and keyword arguments, is the one that is
finally called, and can also be used directly. Here, one specifies the algorithm explicitly.
Currently, only [`CG`](@ref) and [`GMRES`](@ref) are implemented, where `CG` is chosen if
`isposdef == true`. Note that in standard `GMRES` terminology, our parameter `krylovdim` is
referred to as the *restart* parameter, and our `maxiter` parameter counts the number of
outer iterations, i.e. restart cycles. In `CG`, the Krylov subspace is only implicit
because short recurrence relations are being used, and therefore no restarts are required.
Therefore, we pass `krylovdim*maxiter` as the maximal number of CG iterations that can be
used by the `CG` algorithm.
"""
function linsolve end

linsolve(A::AbstractMatrix, b::AbstractVector, a₀::Number = 0, a₁::Number = 1; kwargs...) =
    linsolve(A, b, (zero(a₀) * zero(a₁)) * b, a₀, a₁; kwargs...)

linsolve(f, b, a₀::Number = 0, a₁::Number = 1; kwargs...) =
    linsolve(f, b, (zero(a₀) * zero(a₁)) * b, a₀, a₁; kwargs...)

function linsolve(f, b, x₀, a₀::Number = 0, a₁::Number = 1; kwargs...)
    Tx = promote_type(typeof(x₀))
    Tb = typeof(b)
    Tfx = Core.Compiler.return_type(apply, Tuple{typeof(f),Tx})
    T = promote_type(Core.Compiler.return_type(dot, Tuple{Tb,Tfx}), typeof(a₀), typeof(a₁))
    alg = linselector(f, b, T; kwargs...)
    return linsolve(f, b, x₀, alg, a₀, a₁)
end

function linselector(
    f,
    b,
    T::Type;
    issymmetric::Bool = false,
    ishermitian::Bool = T <: Real && issymmetric,
    isposdef::Bool = false,
    krylovdim::Int = KrylovDefaults.krylovdim,
    maxiter::Int = KrylovDefaults.maxiter,
    rtol::Real = KrylovDefaults.tol,
    atol::Real = KrylovDefaults.tol,
    tol::Real = max(atol, rtol * norm(b)),
    orth = KrylovDefaults.orth,
    verbosity::Int = 0
)
    if (T <: Real && issymmetric) || ishermitian
        isposdef &&
            return CG(maxiter = krylovdim * maxiter, tol = tol, verbosity = verbosity)
        # TODO: implement MINRES for symmetric but not posdef; for now use GRMES
        # return MINRES(krylovdim*maxiter, tol=tol)
        return GMRES(
            krylovdim = krylovdim,
            maxiter = maxiter,
            tol = tol,
            orth = orth,
            verbosity = verbosity
        )
    else
        return GMRES(
            krylovdim = krylovdim,
            maxiter = maxiter,
            tol = tol,
            orth = orth,
            verbosity = verbosity
        )
    end
end
function linselector(
    A::AbstractMatrix,
    b,
    T::Type;
    issymmetric::Bool = T <: Real && LinearAlgebra.issymmetric(A),
    ishermitian::Bool = issymmetric || LinearAlgebra.ishermitian(A),
    isposdef::Bool = ishermitian ? LinearAlgebra.isposdef(A) : false,
    krylovdim::Int = KrylovDefaults.krylovdim,
    maxiter::Int = KrylovDefaults.maxiter,
    rtol::Real = KrylovDefaults.tol,
    atol::Real = KrylovDefaults.tol,
    tol::Real = max(atol, rtol * norm(b)),
    orth = KrylovDefaults.orth,
    verbosity::Int = 0
)
    if (T <: Real && issymmetric) || ishermitian
        isposdef &&
            return CG(maxiter = krylovdim * maxiter, tol = tol, verbosity = verbosity)
        # TODO: implement MINRES for symmetric but not posdef; for now use GRMES
        # return MINRES(krylovdim*maxiter, tol=tol, verbosity = verbosity)
        return GMRES(
            krylovdim = krylovdim,
            maxiter = maxiter,
            tol = tol,
            orth = orth,
            verbosity = verbosity
        )
    else
        return GMRES(
            krylovdim = krylovdim,
            maxiter = maxiter,
            tol = tol,
            orth = orth,
            verbosity = verbosity
        )
    end
end


text:
include("../Tests/Tests.jl")
include("utilities.jl")

using Test, JuMP

@testset "Term" begin
    include("term.jl")
end
@testset "Term fixed" begin
    include("term_fixed.jl")
end
@testset "Quartic constant" begin
    include("quartic_constant.jl")
end
@testset "Quadratic" begin
    include("quadratic.jl")
end
@testset "Quartic ideal" begin
    include("quartic_ideal.jl")
end
@testset "Univariate Sum" begin
    include("univariate_sum.jl")
end
@testset "Rearrangement" begin
    include("rearrangement.jl")
end
@testset "Choi" begin
    include("choi.jl")
end
@testset "Simple matrix" begin
    include("simple_matrix.jl")
end
@testset "Concave then convex cubic" begin
    include("concave_then_convex_cubic.jl")
end
@testset "Horn" begin
    include("horn.jl")
end
@testset "Lyapunov Switched System" begin
    include("lyapunov_switched_system.jl")
end
@testset "Motzkin" begin
    include("motzkin.jl")
end
@testset "BPT12e399" begin
    include("BPT12e399.jl")
end
@testset "Max Cut" begin
    include("maxcut.jl")
end
@testset "Chebyshev" begin
    include("chebyshev.jl")
end
@testset "Quartic comparison" begin
    include("quartic_comparison.jl")
end
@testset "SOSDEMO5" begin
    include("sosdemo5.jl")
end
@testset "SOSDEMO9" begin
    include("sosdemo9.jl")
end
@testset "SOSDEMO10" begin
    include("sosdemo10.jl")
end
@testset "Options Pricing" begin
    include("options_pricing.jl")
end


text:
module VoronoiDelaunay

# Fast, robust 2D Voronoi/Delaunay tesselation
# Implementation follows algorithms described in http://arxiv.org/abs/0901.4107
# and used (for e.g.) in the Illustris Simulation
# http://www.illustris-project.org/
#
# Author: Ariel Keselman (skariel@gmail.com)
# License: MIT
# Bug reportss welcome!

export
	DelaunayTessellation, DelaunayTessellation2D, sizehint!, isexternal,
	min_coord, max_coord, locate, movea, moveb, movec,
	delaunayedges, voronoiedges, voronoiedgeswithoutgenerators,
	start, next, done,
	findindex, push!,
	Point, Point2D, AbstractPoint2D, getx, gety, geta, getb, getc,
	getplotxy

using Compat

using GeometricalPredicates
import GeometricalPredicates
import GeometricalPredicates: geta, getb, getc

import Base: push!, start, done, next, copy, sizehint!
import Colors: RGB, RGBA

const min_coord = GeometricalPredicates.min_coord + eps(Float64)
const max_coord = GeometricalPredicates.max_coord - eps(Float64)

type DelaunayTriangle{T<:AbstractPoint2D} <: AbstractNegativelyOrientedTriangle
  _a::T; _b::T; _c::T
  _bx::Float64; _by::Float64
  _cx::Float64; _cy::Float64
  _px::Float64; _py::Float64
  _pr2::Float64
  _neighbour_a::Int64
  _neighbour_b::Int64
  _neighbour_c::Int64

	function DelaunayTriangle(pa::T, pb::T, pc::T,
			                      na::Int64, nb::Int64, nc::Int64)
      t = new(pa, pb, pc, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, na, nb, nc)
      clean!(t)
      t
  end

	# this constructor is good for making copies
  function DelaunayTriangle(pa::T, pb::T, pc::T,
		                        bx::Float64, by::Float64,
														cx::Float64, cy::Float64,
														px::Float64, py::Float64,
														pr2::Float64,
														na::Int64, nb::Int64, nc::Int64)
		new(pa, pb, pc, bx, by, cx, cy, px, py, pr2, na, nb, nc)
	end
end
function DelaunayTriangle{T<:AbstractPoint2D}(pa::T, pb::T, pc::T,
																							na::Int64, nb::Int64, nc::Int64)
	DelaunayTriangle{T}(pa, pb, pc, na, nb, nc)
end
function DelaunayTriangle{T<:AbstractPoint2D}(pa::T, pb::T, pc::T,
																							bx::Float64, by::Float64,
																							cx::Float64, cy::Float64,
																							px::Float64, py::Float64,
																							pr2::Float64,
																							na::Int64, nb::Int64, nc::Int64)
  DelaunayTriangle{T}(pa, pb, pc, bx, by, cx, cy, px, py, pr2, na, nb, nc)
end

function copy{T<:AbstractPoint2D}(t::DelaunayTriangle{T})
	DelaunayTriangle(
		t._a, t._b, t._c,
		t._bx, t._by,
		t._cx, t._cy,
		t._px, t._py,
		t._pr2,
		t._neighbour_a, t._neighbour_b, t._neighbour_c
	)
end

function isexternal{T<:AbstractPoint2D}(t::DelaunayTriangle{T})
	getx(geta(t)) < min_coord || getx(geta(t)) > max_coord ||
	getx(getb(t)) < min_coord || getx(getb(t)) > max_coord ||
	getx(getc(t)) < min_coord || getx(getc(t)) > max_coord
end

type DelaunayTessellation2D{T<:AbstractPoint2D}
	_trigs::Vector{DelaunayTriangle{T}}
	_last_trig_index::Int64
	_edges_to_check::Vector{Int64}
	_total_points_added::Int64

	function DelaunayTessellation2D(n::Int64 = 100)
		const a = T(GeometricalPredicates.min_coord, GeometricalPredicates.min_coord)
		const b = T(GeometricalPredicates.min_coord, GeometricalPredicates.max_coord)
		const c = T(GeometricalPredicates.max_coord, GeometricalPredicates.min_coord)
		const d = T(GeometricalPredicates.max_coord, GeometricalPredicates.max_coord)
		const t1 = DelaunayTriangle{T}(d,c,b, 2,1,1)
		const t2 = DelaunayTriangle{T}(a,b,c, 3,1,1)
		const t3 = DelaunayTriangle{T}(d,c,b, 2,1,1)
		const _trigs = DelaunayTriangle{T}[t1, t2, t3]
		t = new(_trigs, 3, Int64[], 0)
		sizehint!(t._edges_to_check, 1000)
		sizehint!(t, n)
	end
end
DelaunayTessellation2D(n::Int64) = DelaunayTessellation2D{Point2D}(n)
DelaunayTessellation2D{T<:AbstractPoint2D}(n::Int64, ::T) = DelaunayTessellation2D{T}(n)
DelaunayTessellation(n::Int64=100) = DelaunayTessellation2D(n)

function sizehint!{T<:AbstractPoint2D}(t::DelaunayTessellation2D{T}, n::Int64)
	const required_total_size = 2n + 10
	required_total_size <= length(t._trigs) && return
	sizehint!(t._trigs, required_total_size)
	while length(t._trigs) < required_total_size
		push!(t._trigs, copy(t._trigs[end]))
	end
	t
end

# growing strategy
function sizefit_at_least{T<:AbstractPoint2D}(t::DelaunayTessellation2D{T}, n::Int64)
	const minimal_acceptable_actual_size = 2*n+10
	minimal_acceptable_actual_size <= length(t._trigs) && return
	required_total_size = length(t._trigs)
	while required_total_size < minimal_acceptable_actual_size
		required_total_size += required_total_size >>> 1
	end
	sizehint!(t._trigs, required_total_size)
	while length(t._trigs) < required_total_size
		push!(t._trigs, copy(t._trigs[end]))
	end
	t
end

immutable DelaunayEdge{T<:AbstractPoint2D}
	_a::T
	_b::T
end
geta{T<:AbstractPoint2D}(e::DelaunayEdge{T}) = e._a
getb{T<:AbstractPoint2D}(e::DelaunayEdge{T}) = e._b

immutable VoronoiEdge{T<:AbstractPoint2D}
	_a::Point2D
	_b::Point2D
	_generator_a::T
	_generator_b::T
end
geta{T<:AbstractPoint2D}(e::VoronoiEdge{T}) = e._a
getb{T<:AbstractPoint2D}(e::VoronoiEdge{T}) = e._b
getgena{T<:AbstractPoint2D}(e::VoronoiEdge{T}) = e._generator_a
getgenb{T<:AbstractPoint2D}(e::VoronoiEdge{T}) = e._generator_b

immutable VoronoiEdgeWithoutGenerators
	_a::Point2D
	_b::Point2D
end
geta(e::VoronoiEdgeWithoutGenerators) = e._a
getb(e::VoronoiEdgeWithoutGenerators) = e._b

# TODO: is an iterator faster?
function delaunayedges(t::DelaunayTessellation2D)
	visited = zeros(Bool, t._last_trig_index)
	function delaunayiterator()
		@inbounds for ix in 2:t._last_trig_index
			const tr = t._trigs[ix]
			isexternal(tr) && continue
			visited[ix] && continue
			visited[ix] = true
			const ix_na = tr._neighbour_a
			if !visited[ix_na]
				produce(DelaunayEdge(getb(tr), getc(tr)))
			end
			const ix_nb = tr._neighbour_b
			if !visited[ix_nb]
				produce(DelaunayEdge(geta(tr), getc(tr)))
			end
			const ix_nc = tr._neighbour_c
			if !visited[ix_nc]
				produce(DelaunayEdge(geta(tr), getb(tr)))
			end
		end
	end
	Task(delaunayiterator)
end

# TODO: is an iterator faster?
function voronoiedges(t::DelaunayTessellation2D)
	visited = zeros(Bool, t._last_trig_index)
	visited[1] = true
	function voronoiiterator()
		for ix in 2:t._last_trig_index
			visited[ix] && continue
			const tr = t._trigs[ix]
			visited[ix] = true
			#isexternal(tr) && continue
			const cc = circumcenter(tr)

			const ix_na = tr._neighbour_a
			if !visited[ix_na] #&& !isexternal(t._trigs[ix_na])
				const nb = t._trigs[ix_na]
				produce(VoronoiEdge(cc, circumcenter(nb), getb(tr), getc(tr)))
			end
			const ix_nb = tr._neighbour_b
			if !visited[ix_nb] #&& !isexternal(t._trigs[ix_nb])
				const nb = t._trigs[ix_nb]
				produce(VoronoiEdge(cc, circumcenter(nb), geta(tr), getc(tr)))
			end
			const ix_nc = tr._neighbour_c
			if !visited[ix_nc] #&& !isexternal(t._trigs[ix_nb])
				const nb = t._trigs[ix_nc]
				produce(VoronoiEdge(cc, circumcenter(nb), geta(tr), getb(tr)))
			end
		end
	end
	Task(voronoiiterator)
end

function voronoiedgeswithoutgenerators(t::DelaunayTessellation2D)
	visited = zeros(Bool, t._last_trig_index)
	visited[1] = true
	function voronoiiterator()
		for ix in 2:t._last_trig_index
			visited[ix] && continue
			const tr = t._trigs[ix]
			visited[ix] = true
			#isexternal(tr) && continue
			const cc = circumcenter(tr)

			const ix_na = tr._neighbour_a
			if !visited[ix_na] #&& !isexternal(t._trigs[ix_na])
				const nb = t._trigs[ix_na]
				produce(VoronoiEdgeWithoutGenerators(cc, circumcenter(nb)))
			end
			const ix_nb = tr._neighbour_b
			if !visited[ix_nb] #&& !isexternal(t._trigs[ix_nb])
				const nb = t._trigs[ix_nb]
				produce(VoronoiEdgeWithoutGenerators(cc, circumcenter(nb)))
			end
			const ix_nc = tr._neighbour_c
			if !visited[ix_nc] #&& !isexternal(t._trigs[ix_nc])
				const nb = t._trigs[ix_nc]
				produce(VoronoiEdgeWithoutGenerators(cc, circumcenter(nb)))
			end
		end
	end
	Task(voronoiiterator)
end

type TrigIter
	ix::Int64
end
start(t::DelaunayTessellation2D) = TrigIter(2)
function done(t::DelaunayTessellation2D, it::TrigIter)
	while isexternal(t._trigs[it.ix]) && it.ix <= t._last_trig_index
		it.ix += 1
	end
	it.ix > t._last_trig_index
end
function next(t::DelaunayTessellation2D, it::TrigIter)
	@inbounds const trig = t._trigs[it.ix]
	it.ix += 1
	(trig, it)
end

function findindex{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T}, p::T)
	i::Int64 = tess._last_trig_index
	while true
		@inbounds const w = intriangle(tess._trigs[i], p)
		w > 0 && return i
		@inbounds const tr = tess._trigs[i]
		if w == -1
			i = tr._neighbour_a
		elseif w == -2
			i = tr._neighbour_b
		else
			i = tr._neighbour_c
		end
	end
end

locate{T<:AbstractPoint2D}(t::DelaunayTessellation2D{T}, p::T) = t._trigs[findindex(t, p)]

movea{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
												  trig::DelaunayTriangle{T}) = tess._trigs[trig._neighbour_a]
moveb{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
                          trig::DelaunayTriangle{T}) = tess._trigs[trig._neighbour_b]
movec{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
                          trig::DelaunayTriangle{T}) = tess._trigs[trig._neighbour_c]

function _pushunfixed!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T}, p::T)
	i = findindex(tess, p)
	const ltrigs1::Int64 = tess._last_trig_index+1
	const ltrigs2::Int64 = tess._last_trig_index+2

	@inbounds const t1 = tess._trigs[i]
	old_t1_a = geta(t1)
	seta(t1, p)
	old_t1_b = t1._neighbour_b
	old_t1_c = t1._neighbour_c
	t1._neighbour_b = ltrigs1
	t1._neighbour_c = ltrigs2

	@inbounds const t2 = tess._trigs[ltrigs1]
	setabc(t2, old_t1_a, p, getc(t1))
	t2._neighbour_a = i
	t2._neighbour_b = old_t1_b
	t2._neighbour_c = ltrigs2

	@inbounds const t3 = tess._trigs[ltrigs2]
	setabc(t3, old_t1_a, getb(t1), p)
	t3._neighbour_a = i
	t3._neighbour_b = ltrigs1
	t3._neighbour_c = old_t1_c

	@inbounds const nt2 = tess._trigs[t2._neighbour_b]
	if nt2._neighbour_a==i
		nt2._neighbour_a = ltrigs1
	elseif nt2._neighbour_b==i
		nt2._neighbour_b = ltrigs1
	else
		nt2._neighbour_c = ltrigs1
	end

	@inbounds const nt3 = tess._trigs[t3._neighbour_c]
	if nt3._neighbour_a==i
		nt3._neighbour_a = ltrigs2
	elseif nt3._neighbour_b==i
		nt3._neighbour_b = ltrigs2
	else
		nt3._neighbour_c = ltrigs2
	end

	tess._last_trig_index += 2

	i
end

function _flipa!(tess::DelaunayTessellation2D, ix1::Int64, ix2::Int64)
	@inbounds const ot1 = tess._trigs[ix1]
	@inbounds const ot2 = tess._trigs[ix2]
	if ot2._neighbour_a == ix1
		_flipaa!(tess, ix1, ix2, ot1, ot2)
	elseif ot2._neighbour_b == ix1
		_flipab!(tess, ix1, ix2, ot1, ot2)
	else
		_flipac!(tess, ix1, ix2, ot1, ot2)
	end
end

function _endflipa!(tess::DelaunayTessellation2D,
	                  ix1::Int64, ix2::Int64,
										ot1::DelaunayTriangle, ot2::DelaunayTriangle)
	@inbounds const n1 = tess._trigs[ot1._neighbour_a]
	if n1._neighbour_a==ix2
		n1._neighbour_a = ix1
	elseif n1._neighbour_b==ix2
		n1._neighbour_b = ix1
	else
		n1._neighbour_c = ix1
	end
	@inbounds const n2 = tess._trigs[ot2._neighbour_c]
	if n2._neighbour_a==ix1
		n2._neighbour_a = ix2
	elseif n2._neighbour_b==ix1
		n2._neighbour_b = ix2
	else
		n2._neighbour_c = ix2
	end
end

function _flipaa!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
																			ix1::Int64, ix2::Int64,
																			ot1::DelaunayTriangle{T}, ot2::DelaunayTriangle{T})
	old_ot1_geom_b = getb(ot1)
	setb(ot1, geta(ot2))
	ot1._neighbour_a = ot2._neighbour_c
	old_ot1_neighbour_c = ot1._neighbour_c
	ot1._neighbour_c = ix2

	newc = geta(ot2)
	setabc(ot2, geta(ot1), old_ot1_geom_b, newc)
	ot2._neighbour_a = ot2._neighbour_b
	ot2._neighbour_b = ix1
	ot2._neighbour_c = old_ot1_neighbour_c

	_endflipa!(tess,ix1,ix2,ot1,ot2)
end

function _flipab!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
																			ix1::Int64, ix2::Int64,
																			ot1::DelaunayTriangle{T}, ot2::DelaunayTriangle{T})
	old_ot1_geom_b = getb(ot1)
	setb(ot1, getb(ot2))
	ot1._neighbour_a = ot2._neighbour_a
	old_ot1_neighbour_c = ot1._neighbour_c
	ot1._neighbour_c = ix2

	newc = getb(ot2)
	setabc(ot2, geta(ot1), old_ot1_geom_b, newc)
	ot2._neighbour_a = ot2._neighbour_c
	ot2._neighbour_b = ix1
	ot2._neighbour_c = old_ot1_neighbour_c

	_endflipa!(tess,ix1,ix2,ot1,ot2)
end

function _flipac!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
																			ix1::Int64, ix2::Int64,
																			ot1::DelaunayTriangle{T}, ot2::DelaunayTriangle{T})
	old_ot1_geom_b = getb(ot1)
	setb(ot1, getc(ot2))
	ot1._neighbour_a = ot2._neighbour_b
	old_ot1_neighbour_c = ot1._neighbour_c
	ot1._neighbour_c = ix2

	setab(ot2, geta(ot1), old_ot1_geom_b)
	ot2._neighbour_b = ix1
	ot2._neighbour_c = old_ot1_neighbour_c

	_endflipa!(tess,ix1,ix2,ot1,ot2)
end

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

function _flipb!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
	                                   ix1::Int64, ix2::Int64)
	@inbounds const ot1 = tess._trigs[ix1]
	@inbounds const ot2 = tess._trigs[ix2]
	if ot2._neighbour_a == ix1
		_flipba!(tess, ix1, ix2, ot1, ot2)
	elseif ot2._neighbour_b == ix1
		_flipbb!(tess, ix1, ix2, ot1, ot2)
	else
		_flipbc!(tess, ix1, ix2, ot1, ot2)
	end
end

function _endflipb!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
	                                      ix1::Int64, ix2::Int64,
																				ot1::DelaunayTriangle{T}, ot2::DelaunayTriangle{T})
	@inbounds const n1 = tess._trigs[ot1._neighbour_b]
	if n1._neighbour_a==ix2
		n1._neighbour_a = ix1
	elseif n1._neighbour_b==ix2
		n1._neighbour_b = ix1
	else
		n1._neighbour_c = ix1
	end
	@inbounds const n2 = tess._trigs[ot2._neighbour_a]
	if n2._neighbour_a==ix1
		n2._neighbour_a = ix2
	elseif n2._neighbour_b==ix1
		n2._neighbour_b = ix2
	else
		n2._neighbour_c = ix2
	end
end

function _flipba!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
																		  ix1::Int64, ix2::Int64,
																			ot1::DelaunayTriangle{T}, ot2::DelaunayTriangle{T})
	old_ot1_geom_c = getc(ot1)
	setc(ot1, geta(ot2))
	old_ot1_neighbour_a = ot1._neighbour_a
	ot1._neighbour_a = ix2
	ot1._neighbour_b = ot2._neighbour_c

	setbc(ot2, getb(ot1), old_ot1_geom_c)
	ot2._neighbour_a = old_ot1_neighbour_a
	ot2._neighbour_c = ix1

	_endflipb!(tess,ix1,ix2,ot1,ot2)
end

function _flipbb!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
																			ix1::Int64, ix2::Int64,
																			ot1::DelaunayTriangle{T}, ot2::DelaunayTriangle{T})
	old_ot1_geom_c = getc(ot1)
	setc(ot1, getb(ot2))
	old_ot1_neighbour_a = ot1._neighbour_a
	ot1._neighbour_a = ix2
	ot1._neighbour_b = ot2._neighbour_a

	newa = getb(ot2)
	setabc(ot2, newa, getb(ot1), old_ot1_geom_c)
	ot2._neighbour_a = old_ot1_neighbour_a
	ot2._neighbour_b = ot2._neighbour_c
	ot2._neighbour_c = ix1

	_endflipb!(tess,ix1,ix2,ot1,ot2)
end

function _flipbc!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
																			ix1::Int64, ix2::Int64,
																			ot1::DelaunayTriangle{T}, ot2::DelaunayTriangle{T})
	old_ot1_geom_c = getc(ot1)
	setc(ot1, getc(ot2))
	old_ot1_neighbour_a = ot1._neighbour_a
	ot1._neighbour_a = ix2
	ot1._neighbour_b = ot2._neighbour_b

	newa = getc(ot2)
	setabc(ot2, newa, getb(ot1), old_ot1_geom_c)
	ot2._neighbour_b = ot2._neighbour_a
	ot2._neighbour_a = old_ot1_neighbour_a
	ot2._neighbour_c = ix1

	_endflipb!(tess,ix1,ix2,ot1,ot2)
end

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

function _flipc!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
																		 ix1::Int64, ix2::Int64)
	@inbounds const ot1 = tess._trigs[ix1]
	@inbounds const ot2 = tess._trigs[ix2]
	if ot2._neighbour_a == ix1
		_flipca!(tess, ix1, ix2, ot1, ot2)
	elseif ot2._neighbour_b == ix1
		_flipcb!(tess, ix1, ix2, ot1, ot2)
	else
		_flipcc!(tess, ix1, ix2, ot1, ot2)
	end
end

function _endflipc!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
																				ix1::Int64, ix2::Int64,
																				ot1::DelaunayTriangle{T}, ot2::DelaunayTriangle{T})
	@inbounds const n1 = tess._trigs[ot1._neighbour_c]
	if n1._neighbour_a==ix2
		n1._neighbour_a = ix1
	elseif n1._neighbour_b==ix2
		n1._neighbour_b = ix1
	else
		n1._neighbour_c = ix1
	end
	@inbounds const n2 = tess._trigs[ot2._neighbour_b]
	if n2._neighbour_a==ix1
		n2._neighbour_a = ix2
	elseif n2._neighbour_b==ix1
		n2._neighbour_b = ix2
	else
		n2._neighbour_c = ix2
	end
end

function _flipca!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
																			ix1::Int64, ix2::Int64,
																			ot1::DelaunayTriangle{T}, ot2::DelaunayTriangle{T})
	old_ot1_geom_a = geta(ot1)
	seta(ot1, geta(ot2))
	old_ot1_neighbour_b = ot1._neighbour_b
	ot1._neighbour_b = ix2
	ot1._neighbour_c = ot2._neighbour_c

	newb = geta(ot2)
	setabc(ot2, old_ot1_geom_a, newb, getc(ot1))
	ot2._neighbour_a = ix1
	ot2._neighbour_c = ot2._neighbour_b
	ot2._neighbour_b = old_ot1_neighbour_b

	_endflipc!(tess,ix1,ix2,ot1,ot2)
end

function _flipcb!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
																			ix1::Int64, ix2::Int64,
																			ot1::DelaunayTriangle{T}, ot2::DelaunayTriangle{T})
	old_ot1_geom_a = geta(ot1)
	seta(ot1, getb(ot2))
	old_ot1_neighbour_b = ot1._neighbour_b
	ot1._neighbour_b = ix2
	ot1._neighbour_c = ot2._neighbour_a

	setac(ot2, old_ot1_geom_a, getc(ot1))
	ot2._neighbour_a = ix1
	ot2._neighbour_b = old_ot1_neighbour_b

	_endflipc!(tess,ix1,ix2,ot1,ot2)
end

function _flipcc!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
																			ix1::Int64, ix2::Int64,
																			ot1::DelaunayTriangle{T}, ot2::DelaunayTriangle{T})
	old_ot1_geom_a = geta(ot1)
	seta(ot1, getc(ot2))
	old_ot1_neighbour_b = ot1._neighbour_b
	ot1._neighbour_b = ix2
	ot1._neighbour_c = ot2._neighbour_b

	newb = getc(ot2)
	setabc(ot2, old_ot1_geom_a, newb, getc(ot1))
	ot2._neighbour_c = ot2._neighbour_a
	ot2._neighbour_a = ix1
	ot2._neighbour_b = old_ot1_neighbour_b

	_endflipc!(tess,ix1,ix2,ot1,ot2)
end

function _restoredelaunayhood!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T},
	                                                 ix_trig::Int64)
	@inbounds const center_pt = geta(tess._trigs[ix_trig])

	# `A` - edge
	push!(tess._edges_to_check, ix_trig)
  while length(tess._edges_to_check) > 0
		@inbounds const trix = tess._edges_to_check[end]
		@inbounds const tr_i = tess._trigs[trix]
		@inbounds const nb_a = tr_i._neighbour_a
		if nb_a > 1
			@inbounds const tr_f = tess._trigs[nb_a]
			if incircle(tr_f, center_pt) > 0
				_flipa!(tess, trix, nb_a)
				push!(tess._edges_to_check, nb_a)
				continue
			end
		end
		pop!(tess._edges_to_check)
	end

	# `B` - edge
	push!(tess._edges_to_check, tess._last_trig_index-1)
  while length(tess._edges_to_check) > 0
		@inbounds const trix = tess._edges_to_check[end]
		@inbounds const tr_i = tess._trigs[trix]
		@inbounds const nb_b = tr_i._neighbour_b
		if nb_b > 1
			@inbounds const tr_f = tess._trigs[nb_b]
			if incircle(tr_f, center_pt) > 0
				_flipb!(tess, trix, nb_b)
				push!(tess._edges_to_check, nb_b)
				continue
			end
		end
		pop!(tess._edges_to_check)
	end

	# `C` - edge
	push!(tess._edges_to_check, tess._last_trig_index)
  while length(tess._edges_to_check) > 0
		@inbounds const trix = tess._edges_to_check[end]
		@inbounds const tr_i = tess._trigs[trix]
		@inbounds const nb_c = tr_i._neighbour_c
		if nb_c > 1
			@inbounds const tr_f = tess._trigs[nb_c]
			if incircle(tr_f, center_pt) > 0
				_flipc!(tess, trix, nb_c)
				push!(tess._edges_to_check, nb_c)
				continue
			end
		end
		pop!(tess._edges_to_check)
	end
end

# push a single point. Grows tessellation as needed
function push!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T}, p::T)
	tess._total_points_added += 1
	sizefit_at_least(tess, tess._total_points_added)
	const i = _pushunfixed!(tess, p)
	_restoredelaunayhood!(tess, i)
end

# push an array in given order
function _pushunsorted!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T}, a::Array{T, 1})
	sizehint!(tess, length(a))
	for p in a
		push!(tess, p)
	end
end

# push an array but sort it first for better performance
function push!{T<:AbstractPoint2D}(tess::DelaunayTessellation2D{T}, a::Array{T, 1})
	shuffle!(a)
	mssort!(a)
	_pushunsorted!(tess, a)
end

intensity(c::RGB)  = c.b
intensity(c::RGBA) = c.b
intensity(c) 	     = getfield(c, 1) # Workaround. Gray needs to be imported from images, which would take to long.

# Create DelaunayTessellation with npts points from an image
function from_image(img, npts)

	# placing points in places that represent the image
	pts = Point2D[]
	for i in 1:npts
		x = rand()
		y = rand()
		if intensity(img[@compat Int64(floor(x * size(img)[1])) + 1, @compat Int64(floor(y * size(img)[2])) + 1]) > 0.5
			if rand() < 0.100
				push!(pts, Point2D(1.0 + rand(), 1.0 + rand()))
			end
			continue
		end
		x /= 2
		y /= 2
		push!(pts, Point2D(1.0 + x + 0.3, 2.0 - y*2/3 - 0.3))
	end

	tess = DelaunayTessellation(npts)
	push!(tess, pts)

	tess
end

function getplotxy(edges)
    edges = collect(edges)
    x = Float64[]
    y = Float64[]
    for e in edges
        push!(x, getx(geta(e)))
        push!(x, getx(getb(e)))
        push!(x, NaN)
        push!(y, gety(geta(e)))
        push!(y, gety(getb(e)))
        push!(y, NaN)
    end
    (x, y)
end

end # module VoronoiDelaunay


text:
using BinDeps

@BinDeps.setup

libwsmp = library_dependency("libwsmp")

libdir = joinpath(BinDeps.depsdir(libwsmp),"usr","lib")

provides(SimpleBuild,
         (@build_steps begin
             ChangeDirectory(libdir)
             `make`
         end),libwsmp,os = :Unix)

@BinDeps.install [:libwsmp => :libwsmp]


text:
using GaussianInterpolation, CMBLensing
using GaussianInterpolation: coords, cov!
using BenchmarkTools
using BenchmarkTools: TrialEstimate, TrialJudgement

function print_benchmark(id, results...)
    (id !== nothing) && printstyled(id, "\n", color=:magenta, bold=true)
    for result in results
        show(stdout, "text/plain", result)
    end
    print("\n\n")
end
function print_benchmark(prefix, groups::Vararg{BenchmarkGroup})
    for group_id in sort([key for key in keys(first(groups))])
        id = (prefix === nothing) ? group_id : join([prefix, group_id], '-')
        print_benchmark(id, [group[group_id] for group in groups]...)
    end
end
print_benchmark(groups::Vararg{BenchmarkGroup}) = print_benchmark(nothing, groups...)

macro loop(body)
  quote
    for i = 1:256^2
      $(esc(body))
    end
  end
end

##
suite = BenchmarkGroup()
suite["coords"] = @benchmarkable coords($rand(256, 256))

m = 20
Nside = 32
@unpack f, f̃, ds, n, ϕ = load_sim_dataset(
    θpix=1, Nside=Nside, T=Float64, μKarcminT=6, pol=:I,
	pixel_mask_kwargs=(num_ptsrcs=10, edge_padding_deg=0));

##
X = coords(f)
x = X[1:1, :]
K = make_kernel(rand(Nside, Nside))
out_vec = zeros(m)
out_mat = zeros(m, m)

suite["cov! vector"] = @benchmarkable @loop cov!($out_vec, $K,  $X)
suite["cov! matrix"] = @benchmarkable @loop cov!($out_mat, $K, $X)

# tune!(suite); BenchmarkTools.save("benchmark_params.json", params(suite))
loadparams!(suite, BenchmarkTools.load("benchmark_params.json")[1], :evals, :samples);

results = minimum(run(suite))
##
results_old = BenchmarkTools.load("benchmark.json")[1]
comparison = judge(results, results_old, time_tolerance = 0.05)
print_benchmark(comparison, results)

BenchmarkTools.save("benchmark.json", results)


text:
module Rjasanow

using ..NESSie
using ..NESSie: _etol, cos, cathetus, sign, distance

using LinearAlgebra: norm

export laplacecoll, laplacecoll!

# =========================================================================================
for T in [:ObservationPosition, :InPlane, :InSpace]
    @eval @doc """
        abstract type ObservationPosition end
        struct InPlane <: ObservationPosition end
        stryct InSpace <: ObservationPosition end

    Enum-like representation of the obseration point's position relative to the
    corresponding surface element
    """ $T
end
abstract type ObservationPosition end
struct InPlane <: ObservationPosition end
struct InSpace <: ObservationPosition end

# =========================================================================================
"""
    function laplacepot{T, P <: PotentialType}(
        ptype::Type{P},
        ξ    ::Vector{T},
        elem ::Triangle{T},
        dist ::T
    )

Computes the single or double layer Laplace potential of the given triangle for the given
observation point `ξ`. The latter needs to be projected onto the surface element plane
[[Rja90]](@ref Bibliography).

!!! note
    The result is premultiplied by 4π.

# Arguments
 * `dist` Distance from the original `ξ` to the surface element plane

# Return type
`T`
"""
@inline function laplacepot(
        ptype::Type{P},
        ξ    ::Vector{T},
        elem ::Triangle{T},
        dist ::T
    ) where {T, P <: PotentialType}
    laplacepot(ptype, ξ, elem.v1, elem.v2, elem.normal, dist) +
    laplacepot(ptype, ξ, elem.v2, elem.v3, elem.normal, dist) +
    laplacepot(ptype, ξ, elem.v3, elem.v1, elem.normal, dist)
end

# =========================================================================================
"""
    laplacepot{T, P <: PotentialType}(
        ptype::Type{P},
        ξ::Vector{T},
        x1::Vector{T},
        x2::Vector{T},
        normal::Vector{T},
        dist::T
    )

Computes the single or double layer Laplace potential of the triangle defined by the
observation point `ξ` and two nodes `x1` and `x2` of the surface element. `x2` is required
to be `x1`'s next neighbor in counterclockwise direction. Also, `ξ` needs to be projected
onto the surface element plane  [[Rja90]](@ref Bibliography).

!!! note
    The result is premultiplied by 4π.

# Arguments
 * `normal` Unit normal vector of the surface element
 * `dist` Distance from the original `ξ` to the surface element plane

# Return type
`T`
"""
function laplacepot(
        ptype ::Type{P},
        ξ     ::Vector{T},
        x1    ::Vector{T},
        x2    ::Vector{T},
        normal::Vector{T},
        dist ::T
    ) where {T, P <: PotentialType}
    # Construct triangle sides. Later on, we will use the height h of the triangle at point
    # ξ as well as the angles φ1 abd φ2 between h and the triangle sides extending from ξ.
    u1 = x1 - ξ
    u2 = x2 - ξ
    v  = x2 - x1

    # Compute side lengths
    u1norm = norm(u1)
    u2norm = norm(u2)
    vnorm  = norm(v)

    # Compute the sine of the angle φ1 (φ2) between u1 (u2) and the height, that is, the
    # cosine of the angle θ1 (θ2) between u1 (u2) and v if θ1 (θ2) is an acute or right
    # angle. Otherwise, compute the sine of -φ1 (-φ2), which is again the cosine of θ1
    # (θ2) in these cases.
    # Note that the equation (as given above) uses a polar coordinate system with ξ being
    # the pole and h giving the polar axis. The negative angles are needed whenever the
    # corresponding triangle side lies below the polar axis.
    sinφ1 = max(-one(T), min(one(T), cos(u1, v, u1norm, vnorm)))
    sinφ2 = max(-one(T), min(one(T), cos(u2, v, u2norm, vnorm)))

    # Compute the height of the triangle
    h = cathetus(u1norm, sinφ1)

    # Degenerate triangles
    if max(zero(T), h) < _etol(T) ||      # ξ on the line through v1 and v2 (|φ1| = |φ2| = π/2)
        one(T) - abs(sinφ1) < _etol(T) || # ξ on the line through v1 and v2 (φ1 = π/2)
        one(T) - abs(sinφ2) < _etol(T) || # ξ on the line through v1 and v2 (φ2 = π/2)
        abs(sinφ1 - sinφ2) < _etol(T)     # v1 and v2 on the line through ξ (φ1 = φ2 = 0)
        return zero(T)
    end

    # Since the observation point (projection) lies in the same plane as the surface
    # element, we can decide whether the result of this function is to be added or
    # subtracted from the whole surface triangle's Laplace potential by checking on which
    # side of v the observation point lies. This is equivalent to checking whether the
    # normal of the triangle here and the one of the surface element (which are obviously
    # (anti)parallel) are oriented alike.
    pot = abs(dist) < _etol(T) ?
        laplacepot(ptype, InPlane, sinφ1, sinφ2, h, dist) :
        laplacepot(ptype, InSpace, sinφ1, sinφ2, h, dist)
    sign(u1, u2, normal) * pot
end

# =========================================================================================
"""
    laplacepot{T, P <: PotentialType, O <: ObservationPosition}(
             ::Type{P},
             ::Type{O},
        sinφ1::T,
        sinφ2::T,
        h    ::T,
        d    ::T
    )

Computes the Laplace potential (or its normal derivative) of the triangle with the given
height `h` at the observation point `ξ` (projected onto the surface element plane) and the
sines of the angles ``φ₁`` and ``φ₂`` between `h` and the triangle sides extending from `ξ`
[[Rja90]](@ref Bibliography).

!!! note
    The result is premultiplied by 4π.

# Arguments
 * `d` Distance from the original `ξ` to the surface element plane

# Return type
`T`
"""
@inline function laplacepot(
             ::Type{SingleLayer},
             ::Type{InPlane},
        sinφ1::T,
        sinφ2::T,
        h    ::T,
        d    ::T
    ) where T
    #=
        h/8π * [ln((1 + sin(φ))/(1 - sin(φ)))]   from φ1 to φ2
        = h/8π * [ln((1 + sin(φ2))/(1 - sin(φ2))) - ln((1 + sin(φ1))/(1 - sin(φ1)))]
        = h/8π * ln((1 + sin(φ2))(1 - sin(φ1))/((1 - sin(φ2))(1 + sin(φ1))))

        The equation in the original source misses a factor of 0.5!
    =#
    h * log((1+sinφ2) * (1-sinφ1) / ((1-sinφ2) * (1+sinφ1))) / 2
end

@inline function laplacepot(
             ::Type{SingleLayer},
             ::Type{InSpace},
        sinφ1::T,
        sinφ2::T,
        h    ::T,
        d    ::T
    ) where T
    #=
        h/8π * <1> + d/4π * <2>

        <1>: ln((√(1 - χ² sin²(φ)) + √(1 - χ²) sin(φ))
             / (√(1 - χ² sin²(φ)) - √(1 - χ²) sin(φ)))
             from φ1 to φ2
             = ln((√(1 - χ² sin²(φ2)) + √(1 - χ²) sin(φ2))
               * (√(1 - χ² sin²(φ1)) - √(1 - χ²) sin(φ1))
               / (√(1 - χ² sin²(φ2)) - √(1 - χ²) sin(φ2))
               * √(1 - χ² sin²(φ1)) + √(1 - χ²) sin(φ1))

        <2>: arcsin(χ sin(φ)) - φ   from φ1 to φ2
             = arcsin(χ sin(φ2)) - arcsin(sin(φ2)) - arcsin(χ sin(φ1)) + arcsin(sin(φ1))

        with χ = d / √(d² + h²).
    =#

    d  = abs(d)
    χ2 = d^2 / (d^2 + h^2)
    χ  = √χ2

    # h/8π * <1>
    result = h * log(logterm(χ2, sinφ2) / logterm(χ2, sinφ1)) / 2
    # + d/4π * <2>
    result + d * (asin(χ * sinφ2) - asin(sinφ2) - asin(χ * sinφ1) + asin(sinφ1))
end

@inline function laplacepot(
             ::Type{DoubleLayer},
             ::Type{InPlane},
        sinφ1::T,
        sinφ2::T,
        h    ::T,
        d    ::T
    ) where T
    #=
        1/4π * ∫-1/|ξ-r'|^3 * (ξ-r')⋅n dr'
        = 1/4π * ∫-1/|ξ-r'|^3 * 0 dr'
        = 0
    =#
    zero(T)
end

@inline function laplacepot(
             ::Type{DoubleLayer},
             ::Type{InSpace},
        sinφ1::T,
        sinφ2::T,
        h    ::T,
        d    ::T
    ) where T
    χ  = abs(d) / √(d^2 + h^2)
    sign(d) * (asin(χ * sinφ1) - asin(sinφ1) - asin(χ * sinφ2) + asin(sinφ2))
end

# =========================================================================================
"""
    laplacecoll!{T, P <: PotentialType}(
        ptype   ::Type{P},
        dest    ::DenseArray{T,1},
        elements::Vector{Triangle{T}},
        Ξ       ::Vector{Vector{T}},
        fvals   ::Union{DenseArray{T,1},SubArray{T,1}}
    )

    laplacecoll!{T, P <: PotentialType}(
        ptype   ::Type{P},
        dest    ::DenseArray{T, 2},
        elements::Vector{Triangle{T}},
        Ξ       ::Vector{Vector{T}}
    )

Analytical solution for the single or double layer Laplace potential for a given list of
triangles and observation points `Ξ` [[Rja90]](@ref Bibliography).

The first version of this function uses a vector as destination `dest`, where each element
represents the dot product of the corresponding coefficient matrix row and the `fvals`
vector.

!!! note
    The result is premultiplied by 4π.

# Return type
`Void`
"""
function laplacecoll!(
        ptype   ::Type{P},
        dest    ::DenseArray{T,1},
        elements::Vector{Triangle{T}},
        Ξ       ::Vector{Vector{T}},
        fvals   ::Union{DenseArray{T,1},SubArray{T,1}}
    ) where {T, P <: PotentialType}
    @assert length(dest) == length(Ξ)
    @assert length(fvals) == length(elements)

    @inbounds for eidx in 1:length(elements)
        elem = elements[eidx]
        for oidx in 1:length(Ξ)
            ξ, dist = projectξ(Ξ[oidx], elem)
            dest[oidx] += laplacepot(ptype, ξ, elem, dist) * fvals[eidx]
        end
    end
    nothing
end

function laplacecoll!(
        ptype   ::Type{P},
        dest    ::DenseArray{T, 2},
        elements::Vector{Triangle{T}},
        Ξ       ::Vector{Vector{T}},
    ) where {T, P <: PotentialType}
    @assert size(dest) == (length(Ξ), length(elements))

    @inbounds for eidx in 1:length(elements)
        elem = elements[eidx]
        for oidx in 1:length(Ξ)
            ξ, dist = projectξ(Ξ[oidx], elem)
            dest[oidx, eidx] = laplacepot(ptype, ξ, elem, dist)
        end
    end
    nothing
end

# =========================================================================================
"""
    laplacecoll{T, P <: PotentialType}(
        ptype::Type{P},
        ξ    ::Vector{T},
        elem ::Vector{Triangle{T}}
    )

Analytical solution for the single or double layer Laplace potential for a given triangle
and observation point `ξ` [[Rja90]](@ref Bibliography).

!!! note
    The result is premultiplied by 4π.

# Return type
`Void`
"""
@inline function laplacecoll(
        ptype::Type{P},
        ξ    ::Vector{T},
        elem ::Triangle{T}
    ) where {T, P <: PotentialType}

    ξ_p, dist = projectξ(ξ, elem)
    laplacepot(ptype, ξ_p, elem, dist)
end

# =========================================================================================
"""
    logterm{T}(χ2::T, sinφ::T)

Utility function to compute
```math
\\frac
{\\sqrt{1 - χ² \\sin²(φ)} + \\sqrt{1 - χ² \\sin(φ)}}
{\\sqrt{1 - χ² \\sin²(φ)} - \\sqrt{1 - χ² \\sin(φ)}}
```

# Return type
`T`
"""
@inline function logterm(χ2::T, sinφ::T) where T
    term1 = √(1 - χ2 * sinφ^2)
    term2 = √(1 - χ2) * sinφ
    (term1 + term2) / (term1 - term2)
end

# =========================================================================================
"""
    projectξ{T}(ξ::Vector{T}, elem::Triangle{T})

Projects ξ onto the surface element plane. Returns the projection and the original distance
to the plane.

# Return type
`Tuple{Vector{T}, T}`
"""
function projectξ(ξ::Vector{T}, elem::Triangle{T}) where T
    dist = distance(ξ, elem)
    abs(dist) < _etol(T) && return (ξ, dist)
    # Devectorized version of ξ -= dist * elem.normal
    ([ξ[i] - dist * elem.normal[i] for i in 1:3], dist)
end

end # module


text:
# Similar to SU.jl simplify_rules.jl.

# Set union is commutative and associative like plus, so we can (ab)use the logic for plus.
# Using comma (,) instead of vertical line (|) in P(a | b) because | or my own defined ¦ are considered commutative by SU.

CONDITIONAL_PROBABILITY_RULES = [
    # P(A1 ∪ A2 ∪ A3 ... | B) = P(A1 | B) + P(A2 | B) + P(A3 | B) + ...
    # TODO: Check that sets ~~A are disjoint.
    @acrule P(+(~~A), ~~B) => +([P(a, ~~B...) for a in ~~A]...)
]

SSet = SU.Symbolic{<:Set}
SN = SU.Symbolic{<:Number}

conditional_simplifier() = SU.Chain(CONDITIONAL_PROBABILITY_RULES)

function causal_simplifier()
    return SU.If(SU.istree, SU.Fixpoint(conditional_simplifier()))
end

"""
    identify(G, query)

Given a graph `G` and a `query`, identify the causal effect if possible.
This retuns a probability expression involving only observed quantities
(Pearl, 2009; Corollary 3.4.2).

Unlike Pearl's do-calculus, the implementation in this package is, currently, unlikely to
be complete.
"""
function identify(x)
    simplifiers = [SU.serial_simplifier, causal_simplifier()]
    f = SU.Chain(simplifiers)
    simplified = SU.PassThrough(f)(x)
    return simplified
end


text:
# Multivariate Normal distribution

###########################################################
#
#   Abstract base class for multivariate normal
#
#   Each subtype should provide the following methods:
#
#   - length(d):        vector dimension
#   - mean(d):          the mean vector (in full form)
#   - cov(d):           the covariance matrix (in full form)
#   - invcov(d):        inverse of covariance
#   - logdetcov(d):     log-determinant of covariance
#   - sqmahal(d, x):        Squared Mahalanobis distance to center
#   - sqmahal!(r, d, x):    Squared Mahalanobis distances
#   - gradlogpdf(d, x):     Gradient of logpdf w.r.t. x
#   - _rand!(d, x):         Sample random vector(s)
#
#   Other generic functions will be implemented on top
#   of these core functions.
#
###########################################################

"""

The [Multivariate normal distribution](http://en.wikipedia.org/wiki/Multivariate_normal_distribution)
is a multidimensional generalization of the *normal distribution*. The probability density function of
a d-dimensional multivariate normal distribution with mean vector ``\\boldsymbol{\\mu}`` and
covariance matrix ``\\boldsymbol{\\Sigma}`` is:

```math
f(\\mathbf{x}; \\boldsymbol{\\mu}, \\boldsymbol{\\Sigma}) = \\frac{1}{(2 \\pi)^{d/2} |\\boldsymbol{\\Sigma}|^{1/2}}
\\exp \\left( - \\frac{1}{2} (\\mathbf{x} - \\boldsymbol{\\mu})^T \\Sigma^{-1} (\\mathbf{x} - \\boldsymbol{\\mu}) \\right)
```

We realize that the mean vector and the covariance often have special forms in practice,
which can be exploited to simplify the computation. For example, the mean vector is sometimes
just a zero vector, while the covariance matrix can be a diagonal matrix or even in the form
of ``\\sigma^2 \\mathbf{I}``. To take advantage of such special cases, we introduce a parametric
type `MvNormal`, defined as below, which allows users to specify the special structure of
the mean and covariance.

```julia
struct MvNormal{Cov<:AbstractPDMat,Mean<:AbstractVector} <: AbstractMvNormal
    μ::Mean
    Σ::Cov
end
```

Here, the mean vector can be an instance of any `AbstractVector`. The covariance can be
of any subtype of `AbstractPDMat`. Particularly, one can use `PDMat` for full covariance,
`PDiagMat` for diagonal covariance, and `ScalMat` for the isotropic covariance -- those
in the form of ``\\sigma \\mathbf{I}``. (See the Julia package
[PDMats](https://github.com/lindahua/PDMats.jl) for details).

We also define a set of alias for the types using different combinations of mean vectors and covariance:

```julia
const IsoNormal  = MvNormal{ScalMat,  Vector{Float64}}
const DiagNormal = MvNormal{PDiagMat, Vector{Float64}}
const FullNormal = MvNormal{PDMat,    Vector{Float64}}

const ZeroMeanIsoNormal{Axes}  = MvNormal{ScalMat,  Zeros{Float64,1,Axes}}
const ZeroMeanDiagNormal{Axes} = MvNormal{PDiagMat, Zeros{Float64,1,Axes}}
const ZeroMeanFullNormal{Axes} = MvNormal{PDMat,    Zeros{Float64,1,Axes}}
```

Multivariate normal distributions support affine transformations:
```julia
d = MvNormal(μ, Σ)
c + B * d    # == MvNormal(B * μ + c, B * Σ * B')
dot(b, d)    # == Normal(dot(b, μ), b' * Σ * b)
```
"""
abstract type AbstractMvNormal <: ContinuousMultivariateDistribution end

### Generic methods (for all AbstractMvNormal subtypes)

insupport(d::AbstractMvNormal, x::AbstractVector) =
    length(d) == length(x) && all(isfinite, x)

mode(d::AbstractMvNormal) = mean(d)
modes(d::AbstractMvNormal) = [mean(d)]

"""
    rand(::AbstractRNG, ::Distributions.AbstractMvNormal)

Sample a random vector from the provided multi-variate normal distribution.
"""
rand(::AbstractRNG, ::Distributions.AbstractMvNormal)

function entropy(d::AbstractMvNormal)
    ldcd = logdetcov(d)
    T = typeof(ldcd)
    (length(d) * (T(log2π) + one(T)) + ldcd)/2
end

mvnormal_c0(g::AbstractMvNormal) = -(length(g) * convert(eltype(g), log2π) + logdetcov(g))/2

"""
    invcov(d::AbstractMvNormal)

Return the inversed covariance matrix of d.
"""
invcov(d::AbstractMvNormal)

"""
    logdetcov(d::AbstractMvNormal)

Return the log-determinant value of the covariance matrix.
"""
logdetcov(d::AbstractMvNormal)

"""
    sqmahal(d, x)

Return the squared Mahalanobis distance from x to the center of d, w.r.t. the covariance.
When x is a vector, it returns a scalar value. When x is a matrix, it returns a vector of length size(x,2).

`sqmahal!(r, d, x)` with write the results to a pre-allocated array `r`.
"""
sqmahal(d::AbstractMvNormal, x::AbstractArray)

sqmahal(d::AbstractMvNormal, x::AbstractMatrix) = sqmahal!(Vector{promote_type(partype(d), eltype(x))}(undef, size(x, 2)), d, x)

_logpdf(d::AbstractMvNormal, x::AbstractVector) = mvnormal_c0(d) - sqmahal(d, x)/2

function _logpdf!(r::AbstractArray, d::AbstractMvNormal, x::AbstractMatrix)
    sqmahal!(r, d, x)
    c0 = mvnormal_c0(d)
    for i = 1:size(x, 2)
        @inbounds r[i] = c0 - r[i]/2
    end
    r
end

_pdf!(r::AbstractArray, d::AbstractMvNormal, x::AbstractMatrix) = exp!(_logpdf!(r, d, x))

###########################################################
#
#   MvNormal
#
#   Multivariate normal distribution with mean parameters
#
###########################################################
"""
    MvNormal

Generally, users don't have to worry about these internal details.
We provide a common constructor `MvNormal`, which will construct a distribution of
appropriate type depending on the input arguments.

    MvNormal(sig)

Construct a multivariate normal distribution with zero mean and covariance represented by `sig`.

    MvNormal(mu, sig)

Construct a multivariate normal distribution with mean `mu` and covariance represented by `sig`.

    MvNormal(d, sig)

Construct a multivariate normal distribution of dimension `d`, with zero mean, and an
isotropic covariance matrix corresponding `abs2(sig)*I`.

# Arguments
- `mu::Vector{T<:Real}`: The mean vector.
- `d::Real`: dimension of distribution.
- `sig`: The covariance, which can in of either of the following forms (with `T<:Real`):
    1. subtype of `AbstractPDMat`,
    2. symmetric matrix of type `Matrix{T}`,
    3. vector of type `Vector{T}`: indicating a diagonal covariance as `diagm(abs2(sig))`,
    4. real-valued number: indicating an isotropic covariance matrix corresponding `abs2(sig) * I`.

**Note:** The constructor will choose an appropriate covariance form internally, so that
special structure of the covariance can be exploited.
"""
struct MvNormal{T<:Real,Cov<:AbstractPDMat,Mean<:AbstractVector} <: AbstractMvNormal
    μ::Mean
    Σ::Cov
end

const MultivariateNormal = MvNormal  # for the purpose of backward compatibility

const IsoNormal  = MvNormal{Float64,ScalMat{Float64},Vector{Float64}}
const DiagNormal = MvNormal{Float64,PDiagMat{Float64,Vector{Float64}},Vector{Float64}}
const FullNormal = MvNormal{Float64,PDMat{Float64,Matrix{Float64}},Vector{Float64}}

const ZeroMeanIsoNormal{Axes}  = MvNormal{Float64,ScalMat{Float64},Zeros{Float64,1,Axes}}
const ZeroMeanDiagNormal{Axes} = MvNormal{Float64,PDiagMat{Float64,Vector{Float64}},Zeros{Float64,1,Axes}}
const ZeroMeanFullNormal{Axes} = MvNormal{Float64,PDMat{Float64,Matrix{Float64}},Zeros{Float64,1,Axes}}

### Construction
function MvNormal(μ::AbstractVector{T}, Σ::AbstractPDMat{T}) where {T<:Real}
    dim(Σ) == length(μ) || throw(DimensionMismatch("The dimensions of mu and Sigma are inconsistent."))
    MvNormal{T,typeof(Σ), typeof(μ)}(μ, Σ)
end

function MvNormal(μ::AbstractVector{<:Real}, Σ::AbstractPDMat)
    R = Base.promote_eltype(μ, Σ)
    MvNormal(convert(AbstractArray{R}, μ), convert(AbstractArray{R}, Σ))
end

function MvNormal(μ::AbstractVector, Σ::AbstractPDMat)
    R = Base.promote_eltype(μ, Σ)
    MvNormal(convert(AbstractArray{R}, μ), convert(AbstractArray{R}, Σ))
end

# constructor with general covariance matrix
MvNormal(μ::AbstractVector{<:Real}, Σ::AbstractMatrix{<:Real}) = MvNormal(μ, PDMat(Σ))
MvNormal(μ::AbstractVector{<:Real}, Σ::Diagonal{<:Real}) = MvNormal(μ, PDiagMat(diag(Σ)))
MvNormal(μ::AbstractVector{<:Real}, Σ::UniformScaling{<:Real}) =
    MvNormal(μ, ScalMat(length(μ), Σ.λ))

# constructor with vector of standard deviations
MvNormal(μ::AbstractVector{<:Real}, σ::AbstractVector{<:Real}) = MvNormal(μ, PDiagMat(abs2.(σ)))

# constructor with scalar standard deviation
MvNormal(μ::AbstractVector{<:Real}, σ::Real) = MvNormal(μ, ScalMat(length(μ), abs2(σ)))

# constructor without mean vector
MvNormal(Σ::AbstractVecOrMat{<:Real}) = MvNormal(Zeros{eltype(Σ)}(size(Σ, 1)), Σ)

# special constructor
MvNormal(d::Int, σ::Real) = MvNormal(Zeros{typeof(σ)}(d), σ)

Base.eltype(::Type{<:MvNormal{T}}) where {T} = T

### Conversion
function convert(::Type{MvNormal{T}}, d::MvNormal) where T<:Real
    MvNormal(convert(AbstractArray{T}, d.μ), convert(AbstractArray{T}, d.Σ))
end
function convert(::Type{MvNormal{T}}, μ::AbstractVector, Σ::AbstractPDMat) where T<:Real
    MvNormal(convert(AbstractArray{T}, μ), convert(AbstractArray{T}, Σ))
end

### Show

distrname(d::IsoNormal)  = "IsoNormal"    # Note: IsoNormal, etc are just alias names.
distrname(d::DiagNormal) = "DiagNormal"
distrname(d::FullNormal) = "FullNormal"

distrname(d::ZeroMeanIsoNormal) = "ZeroMeanIsoNormal"
distrname(d::ZeroMeanDiagNormal) = "ZeroMeanDiagNormal"
distrname(d::ZeroMeanFullNormal) = "ZeroMeanFullNormal"

Base.show(io::IO, d::MvNormal) =
    show_multline(io, d, [(:dim, length(d)), (:μ, mean(d)), (:Σ, cov(d))])

### Basic statistics

length(d::MvNormal) = length(d.μ)
mean(d::MvNormal) = d.μ
params(d::MvNormal) = (d.μ, d.Σ)
@inline partype(d::MvNormal{T}) where {T<:Real} = T

var(d::MvNormal) = diag(d.Σ)
cov(d::MvNormal) = Matrix(d.Σ)

invcov(d::MvNormal) = Matrix(inv(d.Σ))
logdetcov(d::MvNormal) = logdet(d.Σ)

### Evaluation

sqmahal(d::MvNormal, x::AbstractVector) = invquad(d.Σ, x .- d.μ)

sqmahal!(r::AbstractVector, d::MvNormal, x::AbstractMatrix) =
    invquad!(r, d.Σ, x .- d.μ)

gradlogpdf(d::MvNormal, x::AbstractVector{<:Real}) = -(d.Σ \ (x .- d.μ))

# Sampling (for GenericMvNormal)

_rand!(rng::AbstractRNG, d::MvNormal, x::VecOrMat) =
    add!(unwhiten!(d.Σ, randn!(rng, x)), d.μ)

# Workaround: randn! only works for Array, but not generally for AbstractArray
function _rand!(rng::AbstractRNG, d::MvNormal, x::AbstractVector)
    for i in eachindex(x)
        @inbounds x[i] = randn(rng, eltype(x))
    end
    add!(unwhiten!(d.Σ, x), d.μ)
end

### Affine transformations

+(d::MvNormal, c::AbstractVector) = MvNormal(d.μ .+ c, d.Σ)

+(c::AbstractVector, d::MvNormal) = d + c

*(B::AbstractMatrix, d::MvNormal) = MvNormal(B * d.μ, X_A_Xt(d.Σ, B))

dot(b::AbstractVector, d::MvNormal) = Normal(dot(d.μ, b), √quad(d.Σ, b))

dot(d::MvNormal, b::AbstractVector) = dot(b, d)

###########################################################
#
#   Estimation of MvNormal
#
###########################################################

### Estimation with known covariance

struct MvNormalKnownCov{Cov<:AbstractPDMat}
    Σ::Cov
end

MvNormalKnownCov(d::Int, σ::Real) = MvNormalKnownCov(ScalMat(d, abs2(Float64(σ))))
MvNormalKnownCov(σ::Vector{Float64}) = MvNormalKnownCov(PDiagMat(abs2.(σ)))
MvNormalKnownCov(Σ::Matrix{Float64}) = MvNormalKnownCov(PDMat(Σ))

length(g::MvNormalKnownCov) = dim(g.Σ)

struct MvNormalKnownCovStats{Cov<:AbstractPDMat}
    invΣ::Cov              # inverse covariance
    sx::Vector{Float64}    # (weighted) sum of vectors
    tw::Float64            # sum of weights
end

function suffstats(g::MvNormalKnownCov{Cov}, x::AbstractMatrix{Float64}) where Cov<:AbstractPDMat
    size(x,1) == length(g) || throw(DimensionMismatch("Invalid argument dimensions."))
    invΣ = inv(g.Σ)
    sx = vec(sum(x, dims=2))
    tw = Float64(size(x, 2))
    MvNormalKnownCovStats{Cov}(invΣ, sx, tw)
end

function suffstats(g::MvNormalKnownCov{Cov}, x::AbstractMatrix{Float64}, w::AbstractVector) where Cov<:AbstractPDMat
    (size(x,1) == length(g) && size(x,2) == length(w)) ||
        throw(DimensionMismatch("Inconsistent argument dimensions."))
    invΣ = inv(g.Σ)
    sx = x * vec(w)
    tw = sum(w)
    MvNormalKnownCovStats{Cov}(invΣ, sx, tw)
end

## MLE estimation with covariance known

fit_mle(g::MvNormalKnownCov{C}, ss::MvNormalKnownCovStats{C}) where {C<:AbstractPDMat} =
    MvNormal(ss.sx * inv(ss.tw), g.Σ)

function fit_mle(g::MvNormalKnownCov, x::AbstractMatrix{Float64})
    d = length(g)
    size(x,1) == d || throw(DimensionMismatch("Invalid argument dimensions."))
    μ = multiply!(vec(sum(x,dims=2)), inv(size(x,2)))
    MvNormal(μ, g.Σ)
end

function fit_mle(g::MvNormalKnownCov, x::AbstractMatrix{Float64}, w::AbstractVector)
    d = length(g)
    (size(x,1) == d && size(x,2) == length(w)) ||
        throw(DimensionMismatch("Inconsistent argument dimensions."))
    μ = BLAS.gemv('N', inv(sum(w)), x, vec(w))
    MvNormal(μ, g.Σ)
end

### Estimation (both mean and cov unknown)

struct MvNormalStats <: SufficientStats
    s::Vector{Float64}  # (weighted) sum of x
    m::Vector{Float64}  # (weighted) mean of x
    s2::Matrix{Float64} # (weighted) sum of (x-μ) * (x-μ)'
    tw::Float64         # total sample weight
end

function suffstats(D::Type{MvNormal}, x::AbstractMatrix{Float64})
    d = size(x, 1)
    n = size(x, 2)
    s = vec(sum(x, dims=2))
    m = s * inv(n)
    z = x .- m
    s2 = z * z'
    MvNormalStats(s, m, s2, Float64(n))
end

function suffstats(D::Type{MvNormal}, x::AbstractMatrix{Float64}, w::AbstractVector)
    d = size(x, 1)
    n = size(x, 2)
    length(w) == n || throw(DimensionMismatch("Inconsistent argument dimensions."))

    tw = sum(w)
    s = x * vec(w)
    m = s * inv(tw)
    z = similar(x)
    for j = 1:n
        xj = view(x,:,j)
        zj = view(z,:,j)
        swj = sqrt(w[j])
        for i = 1:d
            @inbounds zj[i] = swj * (xj[i] - m[i])
        end
    end
    s2 = z * z'
    MvNormalStats(s, m, s2, tw)
end

# Maximum Likelihood Estimation
#
# Specialized algorithms are respectively implemented for
# each kind of covariance
#

fit_mle(D::Type{MvNormal}, ss::MvNormalStats) = fit_mle(FullNormal, ss)
fit_mle(D::Type{MvNormal}, x::AbstractMatrix{Float64}) = fit_mle(FullNormal, x)
fit_mle(D::Type{MvNormal}, x::AbstractMatrix{Float64}, w::AbstractArray{Float64}) = fit_mle(FullNormal, x, w)

fit_mle(D::Type{FullNormal}, ss::MvNormalStats) = MvNormal(ss.m, ss.s2 * inv(ss.tw))

function fit_mle(D::Type{FullNormal}, x::AbstractMatrix{Float64})
    n = size(x, 2)
    mu = vec(mean(x, dims=2))
    z = x .- mu
    C = BLAS.syrk('U', 'N', inv(n), z)
    LinearAlgebra.copytri!(C, 'U')
    MvNormal(mu, PDMat(C))
end

function fit_mle(D::Type{FullNormal}, x::AbstractMatrix{Float64}, w::AbstractVector)
    m = size(x, 1)
    n = size(x, 2)
    length(w) == n || throw(DimensionMismatch("Inconsistent argument dimensions"))

    inv_sw = inv(sum(w))
    mu = BLAS.gemv('N', inv_sw, x, w)

    z = Matrix{Float64}(undef, m, n)
    for j = 1:n
        cj = sqrt(w[j])
        for i = 1:m
            @inbounds z[i,j] = (x[i,j] - mu[i]) * cj
        end
    end
    C = BLAS.syrk('U', 'N', inv_sw, z)
    LinearAlgebra.copytri!(C, 'U')
    MvNormal(mu, PDMat(C))
end

function fit_mle(D::Type{DiagNormal}, x::AbstractMatrix{Float64})
    m = size(x, 1)
    n = size(x, 2)

    mu = vec(mean(x, dims=2))
    va = zeros(Float64, m)
    for j = 1:n
        for i = 1:m
            @inbounds va[i] += abs2(x[i,j] - mu[i])
        end
    end
    multiply!(va, inv(n))
    MvNormal(mu, PDiagMat(va))
end

function fit_mle(D::Type{DiagNormal}, x::AbstractMatrix{Float64}, w::AbstractVector)
    m = size(x, 1)
    n = size(x, 2)
    length(w) == n || throw(DimensionMismatch("Inconsistent argument dimensions"))

    inv_sw = inv(sum(w))
    mu = BLAS.gemv('N', inv_sw, x, w)

    va = zeros(Float64, m)
    for j = 1:n
        @inbounds wj = w[j]
        for i = 1:m
            @inbounds va[i] += abs2(x[i,j] - mu[i]) * wj
        end
    end
    multiply!(va, inv_sw)
    MvNormal(mu, PDiagMat(va))
end

function fit_mle(D::Type{IsoNormal}, x::AbstractMatrix{Float64})
    m = size(x, 1)
    n = size(x, 2)

    mu = vec(mean(x, dims=2))
    va = 0.
    for j = 1:n
        va_j = 0.
        for i = 1:m
            @inbounds va_j += abs2(x[i,j] - mu[i])
        end
        va += va_j
    end
    MvNormal(mu, ScalMat(m, va / (m * n)))
end

function fit_mle(D::Type{IsoNormal}, x::AbstractMatrix{Float64}, w::AbstractVector)
    m = size(x, 1)
    n = size(x, 2)
    length(w) == n || throw(DimensionMismatch("Inconsistent argument dimensions"))

    sw = sum(w)
    inv_sw = 1.0 / sw
    mu = BLAS.gemv('N', inv_sw, x, w)

    va = 0.
    for j = 1:n
        @inbounds wj = w[j]
        va_j = 0.
        for i = 1:m
            @inbounds va_j += abs2(x[i,j] - mu[i]) * wj
        end
        va += va_j
    end
    MvNormal(mu, ScalMat(m, va / (m * sw)))
end


text:
"""
    td_simulation(F; S)

Requires frequency-domain transfer function F(f) and desired time point. F(f)
must be provided as a function of both positive and negative frequencies and the
optional bandwidth must be specified in Hz.
    """
function td_simulation(F,t0;bw=-1e4:1e4)
    f(t) = real(QuadGK.quadgk(f->F(f)*exp(2*pi*im*f*t),bw[1],bw[end])[1])
    return f(t0)
end


text:
using Mimi

@defcomp grosseconomy begin
    YGROSS  = Variable(index=[time])    # Gross output
    K       = Variable(index=[time])    # Capital
    l       = Parameter(index=[time], default=[(1. + 0.015)^t * 6404 for t in 1:20])   # Labor
    tfp     = Parameter(index=[time], default=[(1 + 0.065)^t * 3.57 for t in 1:20])    # Total factor productivity
    s       = Parameter(index=[time], default=(ones(20) * 0.22))   # Savings rate
    depk    = Parameter(default=0.1)    # Depreciation rate on capital - Note that it has no time index
    k0      = Parameter(default=130.0)  # Initial level of capital
    share   = Parameter(default=0.3)    # Capital share

    function run_timestep(p, v, d, t)
        # Define an equation for K
        if is_first(t)
            v.K[t]  = p.k0  # Note the use of v. and p. to distinguish between variables and parameters
        else
            v.K[t]  = (1 - p.depk)^5 * v.K[t-1] + v.YGROSS[t-1] * p.s[t-1] * 5
        end

        # Default values are now provided in @defcomp

        # Define an equation for YGROSS
        v.YGROSS[t] = p.tfp[t] * v.K[t]^p.share * p.l[t]^(1-p.share)
    end
end

@defcomp emissions begin
    E       = Variable(index=[time])    # Total greenhouse gas emissions
    sigma   = Parameter(index=[time], default=[(1. - 0.05)^t * 0.58 for t in 1:20])   # Emissions output ratio
    YGROSS  = Parameter(index=[time])   # Gross output - Note that YGROSS is now a parameter

    function run_timestep(p, v, d, t)
        # Define an eqation for E
        v.E[t] = p.YGROSS[t] * p.sigma[t]   # Note the p. in front of YGROSS
    end
end

model = construct_model()

run(model)

# Show model results
model[:emissions, :E]

# Or, as a DataFrame
getdataframe(model, :emissions, :E)