text:
module Bijectors

using Reexport, Requires
@reexport using Distributions
using StatsFuns
using LinearAlgebra
using MappedArrays
using Roots

export  TransformDistribution,
        PositiveDistribution,
        UnitDistribution,
        SimplexDistribution,
        PDMatDistribution,
        link,
        invlink,
        logpdf_with_trans,
        isclosedform,
        link_jacobian,
        invlink_jacobian,
        transform,
        forward,
        logabsdetjac,
        logabsdetjacinv,
        Bijector,
        ADBijector,
        Inversed,
        Composed,
        compose,
        Stacked,
        stack,
        Identity,
        DistributionBijector,
        bijector,
        transformed,
        TransformedDistribution,
        UnivariateTransformed,
        MultivariateTransformed,
        logpdf_with_jac,
        logpdf_forward,
        PlanarLayer,
        RadialLayer

const DEBUG = Bool(parse(Int, get(ENV, "DEBUG_BIJECTORS", "0")))
_debug(str) = @debug str

_eps(::Type{T}) where {T} = T(eps(T))
_eps(::Type{Real}) = eps(Float64)

#=
  NOTE: Codes below are adapted from
  https://github.com/brian-j-smith/Mamba.jl/blob/master/src/distributions/transformdistribution.jl
  The Mamba.jl package is licensed under the MIT License:
  > Copyright (c) 2014: Brian J Smith and other contributors:
  >
  > https://github.com/brian-j-smith/Mamba.jl/contributors
  >
  > Permission is hereby granted, free of charge, to any person obtaining
  > a copy of this software and associated documentation files (the
  > "Software"), to deal in the Software without restriction, including
  > without limitation the rights to use, copy, modify, merge, publish,
  > distribute, sublicense, and/or sell copies of the Software, and to
  > permit persons to whom the Software is furnished to do so, subject to
  > the following conditions:
  >
  > The above copyright notice and this permission notice shall be
  > included in all copies or substantial portions of the Software.
  >
  > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  > EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  > MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
  > IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
  > CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
  > TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
  > SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#

#############
# a ≦ x ≦ b #
#############

const TransformDistribution{T<:ContinuousUnivariateDistribution} = Union{T, Truncated{T}}
@inline function _clamp(x::T, dist::TransformDistribution) where {T <: Real}
    ϵ = _eps(T)
    bounds = (minimum(dist) + ϵ, maximum(dist) - ϵ)
    clamped_x = ifelse(x < bounds[1], bounds[1], ifelse(x > bounds[2], bounds[2], x))
    DEBUG && _debug("x = $x, bounds = $bounds, clamped_x = $clamped_x")
    return clamped_x
end

link(d::TransformDistribution, x::Real) = _link(d, _clamp(x, d))
function _link(d::TransformDistribution, x::Real)
    a, b = minimum(d), maximum(d)
    lowerbounded, upperbounded = isfinite(a), isfinite(b)
    if lowerbounded && upperbounded
        return StatsFuns.logit((x - a) / (b - a))
    elseif lowerbounded
        return log(x - a)
    elseif upperbounded
        return log(b - x)
    else
        return x
    end
end

invlink(d::TransformDistribution, y::Real) = _clamp(_invlink(d, y), d)
function _invlink(d::TransformDistribution, y::Real)
    a, b = minimum(d), maximum(d)
    lowerbounded, upperbounded = isfinite(a), isfinite(b)
    if lowerbounded && upperbounded
        return (b - a) * StatsFuns.logistic(y) + a
    elseif lowerbounded
        return exp(y) + a
    elseif upperbounded
        return b - exp(y)
    else
        return y
    end
end

function logpdf_with_trans(d::TransformDistribution, x::Real, transform::Bool)
    x = transform ? _clamp(x, d) : x
    return _logpdf_with_trans(d, x, transform)
end
function _logpdf_with_trans(d::TransformDistribution, x::Real, transform::Bool)
    lp = logpdf(d, x)
    if transform
        a, b = minimum(d), maximum(d)
        lowerbounded, upperbounded = isfinite(a), isfinite(b)
        if lowerbounded && upperbounded
            lp += log((x - a) * (b - x) / (b - a))
        elseif lowerbounded
            lp += log(x - a)
        elseif upperbounded
            lp += log(b - x)
        end
    end
    return lp
end

#########
# 0 < x #
#########

const PositiveDistribution = Union{
    BetaPrime, Chi, Chisq, Erlang, Exponential, FDist, Frechet, Gamma, InverseGamma,
    InverseGaussian, Kolmogorov, LogNormal, NoncentralChisq, NoncentralF, Rayleigh, Weibull,
}

_link(d::PositiveDistribution, x::Real) = log(x)
_invlink(d::PositiveDistribution, y::Real) = exp(y)
function _logpdf_with_trans(d::PositiveDistribution, x::Real, transform::Bool)
    return logpdf(d, x) + transform * log(x)
end

#############
# 0 < x < 1 #
#############

const UnitDistribution = Union{Beta, KSOneSided, NoncentralBeta}

_link(d::UnitDistribution, x::Real) = StatsFuns.logit(x)
_invlink(d::UnitDistribution, y::Real) = StatsFuns.logistic(y)
function _logpdf_with_trans(d::UnitDistribution, x::Real, transform::Bool)
    return logpdf(d, x) + transform * log(x * (one(x) - x))
end

###########
# ∑xᵢ = 1 #
###########

const SimplexDistribution = Union{Dirichlet}

_clamp(x, ::SimplexDistribution) = _clamp(x, SimplexBijector())

function link(
    d::SimplexDistribution,
    x::AbstractVecOrMat{<:Real},
    ::Type{Val{proj}} = Val{true}
) where {proj}
    return SimplexBijector{proj}()(x)
end

function link_jacobian(
    d::SimplexDistribution,
    x::AbstractVector{T},
    ::Type{Val{proj}} = Val{true}
) where {T<:Real, proj}
    K = length(x)
    dydxt = similar(x, length(x), length(x))
    @inbounds dydxt .= 0
    ϵ = _eps(T)
    sum_tmp = zero(T)

    @inbounds z = x[1] * (one(T) - 2ϵ) + ϵ # z ∈ [ϵ, 1-ϵ]
    @inbounds dydxt[1,1] = (1/z + 1/(1-z)) * (one(T) - 2ϵ)
    @inbounds @simd for k in 2:(K - 1)
        sum_tmp += x[k - 1]
        # z ∈ [ϵ, 1-ϵ]
        # x[k] = 0 && sum_tmp = 1 -> z ≈ 1
        z = (x[k] + ϵ)*(one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)
        dydxt[k,k] = (1/z + 1/(1-z)) * (one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)
        for i in 1:k-1
            dydxt[i,k] = (1/z + 1/(1-z)) * (x[k] + ϵ)*(one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)^2
        end
    end
    @inbounds sum_tmp += x[K - 1]
    @inbounds if !proj
        @simd for i in 1:K
            dydxt[i,K] = -1
        end
    end

    return UpperTriangular(dydxt)'
end

function invlink(
    d::SimplexDistribution,
    y::AbstractVecOrMat{<:Real},
    ::Type{Val{proj}} = Val{true}
) where {proj}
    return inv(SimplexBijector{proj}())(y)
end

function invlink_jacobian(
    d::SimplexDistribution,
    y::AbstractVector{T},
    ::Type{Val{proj}} = Val{true}
) where {T<:Real, proj}
    K = length(y)
    dxdy = similar(y, length(y), length(y))
    @inbounds dxdy .= 0

    ϵ = _eps(T)
    @inbounds z = StatsFuns.logistic(y[1] - log(T(K - 1)))
    unclamped_x = (z - ϵ) / (one(T) - 2ϵ)
    clamped_x = _clamp(unclamped_x, d)
    @inbounds if unclamped_x == clamped_x
        dxdy[1,1] = z * (1 - z) / (one(T) - 2ϵ)
    end
    sum_tmp = zero(T)
    @inbounds for k = 2:(K - 1)
        z = StatsFuns.logistic(y[k] - log(T(K - k)))
        sum_tmp += clamped_x
        unclamped_x = ((one(T) + ϵ) - sum_tmp) / (one(T) - 2ϵ) * z - ϵ
        clamped_x = _clamp(unclamped_x, d)
        if unclamped_x == clamped_x
            dxdy[k,k] = z * (1 - z) * ((one(T) + ϵ) - sum_tmp) / (one(T) - 2ϵ)
            for i in 1:k-1
                for j in i:k-1
                    dxdy[k,i] += -dxdy[j,i] * z / (one(T) - 2ϵ)
                end
            end
        end
    end
    @inbounds sum_tmp += clamped_x
    @inbounds if proj
    	unclamped_x = one(T) - sum_tmp
        clamped_x = _clamp(unclamped_x, d)
    else
    	unclamped_x = one(T) - sum_tmp - y[K]
        clamped_x = _clamp(unclamped_x, d)
        if unclamped_x == clamped_x
            dxdy[K,K] = -1
        end
    end
    @inbounds if unclamped_x == clamped_x
        for i in 1:K-1
            @simd for j in i:K-1
                dxdy[K,i] += -dxdy[j,i]
            end
        end
    end
    return LowerTriangular(dxdy)
end

function logpdf_with_trans(
    d::SimplexDistribution,
    x::AbstractVector{<:Real},
    transform::Bool,
)
    ϵ = _eps(eltype(x))
    lp = logpdf(d, mappedarray(x -> x + ϵ, x))
    if transform
        lp -= logabsdetjac(bijector(d), x)
    end
    return lp
end

# REVIEW: why do we put this piece of code here?
function logpdf_with_trans(d::Categorical, x::Int)
    return d.p[x] > 0.0 && insupport(d, x) ? log(d.p[x]) : eltype(d.p)(-Inf)
end

###############
# MvLogNormal #
###############

using Distributions: AbstractMvLogNormal

link(d::AbstractMvLogNormal, x::AbstractVector{<:Real}) = log.(x)
invlink(d::AbstractMvLogNormal, y::AbstractVector{<:Real}) = exp.(y)
function logpdf_with_trans(
    d::AbstractMvLogNormal,
    x::AbstractVector{<:Real},
    transform::Bool,
)
    return logpdf(d, x) + transform * sum(log, x)
end

#####################
# Positive definite #
#####################

const PDMatDistribution = Union{InverseWishart, Wishart}

function link(d::PDMatDistribution, X::AbstractMatrix{<:Real})
    Y = Matrix(cholesky(X).L)
    Y[diagind(Y)] .= log.(view(Y, diagind(Y)))
    return Y
end

function invlink(d::PDMatDistribution, Y::AbstractMatrix{<:Real})
    X = copy(Y)
    X[diagind(X)] .= exp.(view(X, diagind(X)))
    return LowerTriangular(X) * LowerTriangular(X)'
end

function logpdf_with_trans(
    d::PDMatDistribution,
    X::AbstractMatrix{<:Real},
    transform::Bool
)
    T = eltype(X)
    Xcf = cholesky(X, check=false)
    if !issuccess(Xcf)
        Xcf = cholesky(X + (eps(T) * norm(X)) * I)
    end
    lp = getlogp(d, Xcf, X)
    if transform && isfinite(lp)
        U = Xcf.U
        @inbounds @simd for i in 1:dim(d)
            lp += (dim(d) - i + 2) * log(U[i, i])
        end
        lp += dim(d) * log(T(2))
    end
    return lp
end
function getlogp(d::Wishart, Xcf, X)
    return 0.5 * ((d.df - (dim(d) + 1)) * logdet(Xcf) - tr(d.S \ X)) - d.c0
end
function getlogp(d::InverseWishart, Xcf, X)
    Ψ = Matrix(d.Ψ)
    return -0.5 * ((d.df + dim(d) + 1) * logdet(Xcf) + tr(Xcf \ Ψ)) - d.c0
end

############################################
# Defaults (assume identity link function) #
############################################

# UnivariateDistributions
using Distributions: UnivariateDistribution

link(d::UnivariateDistribution, x::Real) = x
link(d::UnivariateDistribution, x::AbstractVector{<:Real}) = map(x -> link(d, x), x)

invlink(d::UnivariateDistribution, y::Real) = y
invlink(d::UnivariateDistribution, y::AbstractVector{<:Real}) = map(y -> invlink(d, y), y)

logpdf_with_trans(d::UnivariateDistribution, x::Real, ::Bool) = logpdf(d, x)
function logpdf_with_trans(
    d::UnivariateDistribution,
    x::AbstractVector{<:Real},
    transform::Bool,
)
    return map(x -> logpdf_with_trans(d, x, transform), x)
end

# MultivariateDistributions
using Distributions: MultivariateDistribution

link(d::MultivariateDistribution, x::AbstractVecOrMat{<:Real}) = copy(x)

invlink(d::MultivariateDistribution, y::AbstractVecOrMat{<:Real}) = copy(y)

function logpdf_with_trans(d::MultivariateDistribution, x::AbstractVector{<:Real}, ::Bool)
    return logpdf(d, x)
end
function logpdf_with_trans(
    d::MultivariateDistribution,
    X::AbstractMatrix{<:Real},
    transform::Bool,
)
    return [logpdf_with_trans(d, view(X, :, n), transform) for n in 1:size(X, 2)]
end

# MatrixDistributions
using Distributions: MatrixDistribution

link(d::MatrixDistribution, X::AbstractMatrix{<:Real}) = copy(X)
link(d::MatrixDistribution, X::AbstractVector{<:AbstractMatrix{<:Real}}) = map(x -> link(d, x), X)

invlink(d::MatrixDistribution, Y::AbstractMatrix{<:Real}) = copy(Y)
function invlink(d::MatrixDistribution, Y::AbstractVector{<:AbstractMatrix{<:Real}})
    return map(y -> invlink(d, y), Y)
end

logpdf_with_trans(d::MatrixDistribution, X::AbstractMatrix{<:Real}, ::Bool) = logpdf(d, X)
function logpdf_with_trans(
    d::MatrixDistribution,
    X::AbstractVector{<:AbstractMatrix{<:Real}},
    transform::Bool,
)
    return map(x -> logpdf_with_trans(d, x, transform), X)
end

include("interface.jl")

# optional dependencies
function __init__()
    @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" include("compat/forwarddiff.jl")
    @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" include("compat/tracker.jl")
    @require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" include("compat/zygote.jl")
end

end # module


text:
export SpinHalfSite,
       spinHalfSites

struct SpinHalfSite <: Site
    s::Index
    SpinHalfSite(I::Index) = new(I)
end
SpinHalfSite(n::Int) = SpinHalfSite(Index(2, "Site,S=1/2,n=$n"))

function spinHalfSites(N::Int;kwargs...)::SiteSet
  sites = SiteSet(N)
  for n=1:N
    set(sites,n,SpinHalfSite(n))
  end
  return sites
end

function state(site::SpinHalfSite,
               st::String)::IndexVal
  if st == "Up" || st == "↑"
    return site.s(1)
  elseif st == "Dn" || st == "↓"
    return site.s(2)
  else
    throw(ArgumentError("State string \"$st\" not recognized for SpinHalfSite"))
  end
  return site.s(1)
end

function operator(site::SpinHalfSite,
                  opname::AbstractString)::ITensor
    s = site.s
    sP = prime(site.s)
    Up = s(1)
    UpP = sP(1)
    Dn = s(2)
    DnP = sP(2)

    Op = ITensor(dag(s), s')

    if opname == "S⁺" || opname == "Splus" || opname == "S+"
        Op[Dn, UpP] = 1.
    elseif opname == "S⁻" || opname == "Sminus" || opname == "S-"
        Op[Up, DnP] = 1.
    elseif opname == "Sˣ" || opname == "Sx"
        Op[Up, DnP] = 0.5
        Op[Dn, UpP] = 0.5
    elseif opname == "iSʸ" || opname == "iSy"
        Op[Up, DnP] = -0.5
        Op[Dn, UpP] = 0.5
    elseif opname == "Sʸ" || opname == "Sy"
        Op = complex(Op)
        Op[Up, DnP] = 0.5*im
        Op[Dn, UpP] = -0.5*im
    elseif opname == "Sᶻ" || opname == "Sz"
        Op[Up, UpP] = 0.5
        Op[Dn, DnP] = -0.5
    elseif opname == "projUp"
        Op[Up, UpP] = 1.
    elseif opname == "projDn"
        Op[Dn, DnP] = 1.
    elseif opname == "Up" || opname == "↑"
        pU = ITensor(s)
        pU[Up] = 1.
        return pU
    elseif opname == "Dn" || opname == "↓"
        pD = ITensor(s)
        pD[Dn] = 1.
        return pD
    else
        throw(ArgumentError("Operator name '$opname' not recognized for SpinHalfSite"))
    end
    return Op
end

ext: jl


text:
"""
Frames is a module to work with coordinate systems. The central type is the `Frame` with fields `R` and `t` for the rotation matrix and translation vector respectively. Other inportnat important types include `Point, Points, Plane` and `Line`

All geometric object types have a reference frame associated with them. When a geometric object is created, the reference frame is specified, e.g., `T_RB_T = Frame(T_RB_Tm,"RB","T")` is a transformation between frame `RB` and frame `T`.
When transofrmations are made, the reference frames are updated automatically, e.g., `T_RB_SEAM = T_RB_T*T_T_SEAM`

Usage:
```julia
using Robotlib.Frames
path = Pkg.dir("Robotlib","src","applications","frames/")
add_frame_name!("SEAM","Weld seam frame")
add_frame_name!("TAB","Table frame")

T_RB_Tm    = MAT.matread(path*"T_RB_T.mat")["T_RB_T"]
T_TF_TCPm  = MAT.matread(path*"T_TF_TCP.mat")["T_TF_TCP"]
T_T_TABm   = MAT.matread(path*"T_T_Table.mat")["T_T_Table"]

T_RB_T     = Frame(T_RB_Tm,"RB","T")
T_S_D      = Frame(T_TF_TCPm,"S","D")
T_T_TAB    = Frame(T_T_TABm,"T","TAB")

cloud_seam = readcloud(path*"CloudSeam_edge.txt")
plane_seam = readplane(path*"PlaneSeam_edge.txt")
cloud_seam_projected = project(plane_seam,cloud_seam)
line_seam  = fitline(cloud_seam_projected)

T_T_SEAM      = framefromfeatures(("z+",line_seam),("y-",plane_seam),cloud_seam_projected[1],"SEAM")
T_RB_SEAM     = T_RB_T*T_T_SEAM
T_RB_TAB      = T_RB_T*T_T_TAB
T_TAB_SEAM    = inv(T_T_TAB)*T_T_SEAM

cloud_seam_RB = T_RB_T*cloud_seam
cloud_seam_projected_RB = T_RB_T*cloud_seam_projected
plane_seam_RB = T_RB_T*plane_seam
line_seam_RB  = T_RB_T*line_seam

using Plots
default(markersize=1)
plot(Frame(I4,"RB","U"),200, label=true)
plot!(cloud_seam_RB)
plot!(cloud_seam_projected_RB)
plot!(line_seam_RB,500,label="Line seam")
plot!(plane_seam_RB,200,label="Plane seam")
plot!(T_RB_SEAM,200, label=true)
plot!(T_RB_TAB,200, label=true)

xlabel!("x")
ylabel!("y")
zlabel!("z")

MAT.matwrite(path*"T_TAB_SEAM.mat",["T_TAB_SEAM" => T_TAB_SEAM.T])
MAT.matwrite(path*"T_T_SEAM.mat",["T_T_SEAM" => T_T_SEAM.T])
MAT.matwrite(path*"T_RB_TAB.mat",["T_RB_TAB" => T_RB_TAB.T])
println("Wrote T_TAB_SEAM, T_T_SEAM, T_RB_TAB to files in \$path")
```
"""
module Frames
using StaticArrays, RecipesBase, LinearAlgebra, Statistics, SparseArrays
using ..Robotlib

export Frame, Point, Plane, Points, Line, GeometricObject, add_frame_name!
export readcloud, readTmatrix, readplane, fitline, fitplane, framefromfeatures, project
export plot3Dsmart
export dot

import Base: print, zeros, length, size, getindex, setindex!, convert, promote_rule, push!, show, display, +, *, /,  -, \, inv
import LinearAlgebra: ×, transpose, adjoint, dot, ⋅, det
# using LaTeXStrings
import Robotlib: T2R, T2t, I3, I4

const Vect{Ty} = Union{AbstractVector{Ty}, MVector{3,Ty}}
const Matr{Ty} = Union{AbstractMatrix{Ty}, MMatrix{3,3,Ty}}
const Vec{N,Ty} = MVector{N,Ty} # Translate from FixedSizeArrays
const Mat{N1,N2,Ty} = MMatrix{N1,N2,Ty} # Translate from FixedSizeArrays

# Type definitions --------------------------------------------
# -------------------------------------------------------------
abstract type GeometricObject end

mutable struct Point{Ty} <: GeometricObject
    p::Vec{3,Ty}
    A::String
end
function Point(p::Vect{Ty} = zeros(3), A::String="U") where Ty
    checkframe(A,"U")
    Point{Ty}(p,A)
end

mutable struct Points{Ty} <: GeometricObject
    p::Vector{Point{Ty}}
    A::String
end
function Points(p::AbstractVector{Point{Ty}}, A::String="U") where Ty
    checkframe(A,"U")
    Points{Ty}(p,A)
end

function Points(p::AbstractMatrix{Ty}, A::String="U") where Ty
    checkframe(A,"U")
    N = size(p,1)
    points = [Point{Ty}(squeeze(p[i,:]',2),A) for i in 1:N]
    Points{Ty}(points,A)
end
function Points(A::String="U")
    Points{Float64}(Point{Float64}[],A)
end

push!(points::Points, p::Point) = push!(points.p,p)

getindex(points::Points,i) = points.p[i]
function setindex!(points::Points, p::Point, i)
    points.p[i] = p
end
getindex(point::Point,i) = point.p[i]
function setindex!(point::Point, p::Float64, i)
    point.p[i] = p
end

size(points::Points) = size(points.p)[1]

length(points::Points) = length(points.p)

@static if VERSION >= v"0.7"
    Base.iterate(points::Points, state) = iterate(points.p, state)
    Base.iterate(points::Points) = iterate(points.p)
else
    Base.start(points::Points) = Base.start(points.p)
    Base.next(points::Points, state) = Base.next(points.p,state)
    Base.done(points::Points, state) = Base.done(points.p,state)
end

function convert(::Type{Matrix{Float64}}, p::Points)
    N = length(p)
    points = zeros(N,3)
    for i = 1:N
        points[i,:] = Vector(p[i].p)
    end
    return points
end

function zeros(points::Points)
    p = Point()
    N = length(points)
    newpoints = [p for i in 1:N]
    Points(newpoints,points.A)
end

mutable struct Frame{Ty} <: GeometricObject
    R::Mat{3,3,Ty}
    t::Vec{3,Ty}
    A::String
    B::String
end
function Frame(T::Matr{Ty}, A::String="U",B::String="U") where Ty
    checkframe(A,B)
    f = Frame{Ty}(T2R(T),T2t(T),A,B)
    add_frame!(f)
    f
end
function Frame(R::Matr{Ty}, t::AbstractArray{Ty}, A::String="U", B::String="U") where Ty
    checkframe(A,B)
    f = Frame{Ty}(R,t,A,B)
    add_frame!(f)
    f
end
Frame() = Frame{Float64}(I3,zeros(3),"U","U")

ref_frames = Set{Frame}()
frame_map = Dict{String,String}()
frame_names = Dict("U" => "Undefined", "RB"=>"Robot Base","TF"=>"Tool Flange","TCP"=>"Tool Center Point","T"=>"Tracker","D"=>"Dynamic diods","S"=>"Sensor")
add_frame_name!(key::String,value::String) = setindex!(frame_names,value,key)
frame_name_exists(ref::String) = ref ∈ keys(frame_names)
function add_frame!(f::Frame)
    push!(ref_frames,f)
    setindex!(frame_map,f.B,f.A)
end

function checkframe(A,B)
    if !frame_name_exists(A)
        error("Frame \"$A\" does not exist, add to dictionary using add_ref_frame()")
    elseif !frame_name_exists(B)
        error("Frame \"$B\" does not exist, add to dictionary using add_ref_frame()")
    end
end

checkframes(x::GeometricObject, y::GeometricObject) = x.A != y.A &&  error("The two geometric objects do not have the same reference frame ($x -> $(x.A), $y -> $(y.A))")

function display(f::Frame)
    println(round.(F2T(f),4))
    show(string(frame_names[f.A])*"->"*string(frame_names[f.B]))
end
print(f::Frame) = "\$T_{$(f.A)}^{$(f.B)}\$"
show(f::Frame) = display(f)

mutable struct Plane{Ty} <: GeometricObject
    n::Point{Ty}
    r::Point{Ty}
    A::String
end

Plane(n::Vect{Ty},r::Vect{Ty},A="U") where {Ty} = Plane{Ty}(Point{Ty}(n,A),Point{Ty}(r,A),A)
Plane(points::Points) = fitplane(points)

mutable struct Line{Ty} <: GeometricObject
    v::Point{Ty}
    r::Point{Ty}
    A::String
end
Line(v::Vect{Ty},r::Vect{Ty},A="U") where {Ty} = Line{Ty}(v,r,A)
Line(points::Points) = fitline(points)

# Plot functions ----------------------------------------------
# -------------------------------------------------------------

mat2tup(m) = m[:,1][:], m[:,2][:], m[:,3][:]

@recipe function plotframe(f::Frame, length=1.0)#; annotation=false)
    #Plots a frame using XYZ-RGB
    o = F2t(f)
    x = Rx(f)
    y = Ry(f)
    z = Rz(f)
    seriestype := :path3d
    label     := ""
    @series begin
        data = [o o+x*length]'
        seriescolor := :red
        data |> mat2tup
    end
    @series begin
        data = [o o+y*length]'
        seriescolor := :green
        data |> mat2tup
    end
    # if annotation && f.A != ""
    #     warn("Annotations not currently supported for 3d plots")
    #     # po = o-length/4*(x+y+z)
    #     # annotations := (po[1],po[2],po[3],print(f))
    # end
    # delete!(d,:annotation)
    @series begin
        data = [o o+z*length]'
        seriescolor := :blue
        data |> mat2tup
    end
end

@recipe function plotline(l::Line, length=1.0)#; annotation="")
    coords = [Vector((l.r-length*l.v).p)  Vector((l.r+length*l.v).p)]'
    @series begin
        seriestype := :scatter3d
        ([l.r[1]],[l.r[2]],[l.r[3]])
    end
    # if annotation != ""
    #     warn("Annotations not currently supported for 3d plots")
    #     # po = l.r.p-length/4*normalized(l.r.p)
    #     # annotations := (po[1],po[2],po[3],d[:annotation])
    # end
    # delete!(d,:annotation)
    @series begin
        seriestype := :path3d
        coords |> mat2tup
    end
end

@recipe function plotplane(p::Plane, length=1.0)
    coords = [Vector((p.r).p) Vector((p.r+length*p.n).p)]'
    @series begin
        seriestype := :scatter3d
        markershape --> :utriangle
        ([p.r[1]],[p.r[2]],[p.r[3]])
    end
    @series begin
        seriestype := :path3d
        coords |> mat2tup
    end
    # if label != ""
    #     po = p.r.p-length/4*p.n.p
    #     annotate!(po[1],po[2],po[3],label)
    # end
end

@recipe function plotpoints(p::Points)
    seriestype := :scatter3d
    convert(Matrix{Float64},p) |> mat2tup
end

# Frame functions  --------------------------------------------
# -------------------------------------------------------------

*(p::Point, c::Real) = Point(c.*p.p, p.A)
*(c::Real, p::Point) = *(p, c)
*(p::Point, R::Matr) = error("Post multiplication of ::Point with ::Matrix not supported")
*(R::Matr, p::Point) = Point(R*p.p, p.A)
/(p::Point, c::Real) = Point(p.p./c, p.A)
/(c::Real,p::Point) = /(p,c)
+(p1::Point, p2::Point) = Point(p1.p + p2.p, p1.A)
-(p1::Point, p2::Point) = Point(p1.p - p2.p, p1.A)
+(p1::Vect, p2::Point) = p1 + p2.p
-(p1::Vect, p2::Point) = p1 - p2.p
+(p2::Point, p1::Vect) = Point(p1 + p2.p, p2.A)
-(p2::Point, p1::Vect) = Point(p1 - p2.p, p2.A)
⋅(p1::Point, p2::Point) = p1.p ⋅ p2.p
×(p1::Point, p2::Point) = Point(p1.p × p2.p, p1.A)

transpose(p::Point) = p.p'
adjoint(p::Point) = p.p'
convert(::Type{Vector}, p::Point) = p.p
convert(::Type{Array{T1,n}}, f::Frame{T2}) where {T1,T2,n} = [Matrix{T1}(f.R) Vector{T1}(f.t); 0 0 0 1]
convert(::Type{Matrix}, f::Frame{T2}) where {T2} = convert(Matrix{T2},f)
convert(::Type{Array}, f::Frame{T2}) where {T2} = convert(Matrix{T2},f)
promote_rule(::Type{Array{T1,n}}, f::Type{Frame{T2}} ) where {T1,T2,n} = Array{promote_type(T1,T2),2}

*(x::Frame, y::Frame) = Frame(x.R*y.R, x.t + x.R*y.t , x.A, y.B)
function *(f::Frame, p::Point)
    f.B != p.A &&  error("Reference frames does not match between $f and $p ($(f.B) != $(p.A))")
    pv = f.R*p.p + f.t
    Point(pv,f.A)
end
function *(f::Frame, p::Plane)
    f.B != p.A && error("Reference frames does not match between $f and $p ($(f.B) != $(p.A))")
    n = F2R(f)*p.n
    r = f*p.r
    λ = (r⋅n)/(n⋅n)
    r = λ*n
    Plane(n,r,f.A)
end
function *(f::Frame, l::Line)
    f.B != l.A && error("Reference frames does not match between $f and $l ($(f.B) != $(l.A))")
    v = F2R(f)*l.v
    r2 = f*l.r
    r3 = v×r2
    r = r3×v
    Line(v,r,f.A)
end
function *(f::Frame, points::Points)
    f.B != points.A &&  error("Reference frames does not match between $f and $p ($(f.B) != $(p.A))")
    newpoints = deepcopy(points)
    newpoints.A = f.A
    for i in 1:length(newpoints)
        newpoints[i] = f*newpoints[i]
    end
    newpoints
end
*(p::Point, f::Frame) = error("Post multiplication of ::Point with ::Frame not supported")
tinv(T)       = [T[1:3,1:3]' -T[1:3,1:3]'*T[1:3,4]; 0 0 0 1]
function inv(x::Frame{T}) where T # TODO: can this implementation be faster?
    R = x.R'
    Frame{T}(R, -R*x.t, x.B,x.A)
end
\(x::Frame{T},y::Frame{T}) where {T} = Frame{T}(inv(x)*y,x.B,y.B)
/(x::Frame{T},y::Frame{T}) where {T} = Frame{T}(x*inv(y),x.A,y.A)
F2t(F::Frame) = F.t
F2R(F::Frame) = F.R
F2T(F::Frame) = [Array(F.R) Array(F.t); 0 0 0 1]
Rx(F::Frame)  = F.R[1:3,1]
Ry(F::Frame)  = F.R[1:3,2]
Rz(F::Frame)  = F.R[1:3,3]
det(F::Frame) = det(F.R)
normalized(v::Vect) = v/norm(v)
function normalize!(v::Vect) v/=norm(v); end
normalized(p::Point) = Point(p.p/norm(p.p), p.A)
function normalize!(p::Point) p.p/=norm(p.p); end

function ⊥(R::Matr)
    (U,S,V) = svd(R)
    a = diagm(0=>[1, 1, sign(det(U*V'))])
    return U*a*V'
end

⊥(F::Frame) = Frame(⊥(F2R(F2)),F.t,F.A,F.B)

# Projections and fitting  ------------------------------------
# -------------------------------------------------------------

function center(points::Points{Ty}) where Ty
    N = size(points)
    cog = zeros(Ty,3)
    for i = 1:N
        cog += points[i]
    end
    cog /= N

    cpoints = zeros(Ty,length(points),3)
    for i in 1:N
        cpoints[i,:] = Vector(- cog + points[i])
    end
    return (cog,cpoints)
end

function fitplane(points::Points)
    (cog,cpoints) = center(points)
    (U,S,V) = svd(cpoints)

    n = V[:,3]
    r = n*dot(cog,n)[1]

    @assert abs(norm(n)-1) < 1e-8
    @assert abs(abs(normalized(r-cog)⋅n)) < 1e-14
    @assert abs(abs(n⋅normalized(r))-1) < 1e-14
    @assert norm(r) < norm(cog)

    return Plane(n,r,points.A)
end

function fitline(points::Points)
    (cog,cpoints) = center(points)
    (U,S,V) = svd(cpoints)
    v = V[:,1]
    r = cog-dot(v,cog).*v

    @assert abs(norm(v)-1) < 1e-8
    @assert abs(abs(normalized(r-cog)⋅v) - 1) < 1e-14
    @assert abs(v⋅normalized(r)) < 1e-14
    @assert norm(r) < norm(cog)

    return Line(Point(v,points.A), Point(r,points.A),points.A)
end

function project(plane::Plane, points::Points)
    checkframes(plane,points)
    N = size(points)
    n = plane.n
    P = n.p*n.p'
    ppoints = zeros(points)
    for i = 1:N
        ppoints[i] = points[i] + (- P*points[i] + plane.r)
    end
    return ppoints
end

function project(plane::Plane, line::Line)
    checkframes(plane,line)
    n = plane.n
    n_projplane = normalized(n × line.v)
    P = n.p*n.p'
    vp = normalized(line.v - P*line.v)
    @assert abs(vp⋅n_projplane) < 1e-14

    A = [vp',n', n_projplane']
    b = [0, plane.r⋅n, line.r⋅n_projplane]
    rp = Point(A\b,line.A)
    @assert abs(rp ⋅ vp) < 1e-14
    @assert abs(vp ⋅ n) < 1e-14
    return Line(vp, rp)
end

function project(line::Line, points::Points)
    checkframes(points,line)
    N = size(points)
    v = line.v
    ppoints = zeros(points)
    for i = 1:N
        ppoints[i] = ((points[i]⋅v).*v  + line.r)
    end
    return ppoints
end

function framefromfeatures(feature1, feature2, origin, B)
    checkframes(feature1[2],feature2[2])
    x = y = z = false

    r1 = findfirst(lowercase(feature1[1][1]) .== ['x','y','z'])
    r2 = findfirst(lowercase(feature2[1][1]) .== ['x','y','z'])
    if r1 == 0 || r2 == 0; error("Your feature directions ($(feature1[1]), $(feature2[1])) are not recognized"); end
    r3 = 6-r1-r2
    s1 = feature1[1][2] == '+' ? 1 : -1
    s2 = feature2[1][2] == '+' ? 1 : -1

    R1 = Frames.normalized(Frames.getdir(feature1[2]).p * s1) |> Vector
    R2 = Frames.normalized(Frames.getdir(feature2[2]).p * s2) |> Vector
    R3 = Frames.normalized(R1×R2) |> Vector
    R2 = R3×R1

    _R = zeros(3,3)
    _R[:,r1] = R1
    _R[:,r2] = R2
    _R[:,r3] = R3

    (det(_R) < 0) &&  (_R[:,r3] *= -1)
    @assert abs(det(_R)-1) < 1e-14
    R = Mat{3,3}(_R)

    Frame(R,origin.p,feature1[2].A, B)
end

function getdir(f::GeometricObject)
    if isa(f,Plane)
        return f.n
    elseif isa(f,Line)
        return f.v
    elseif isa(f,Frame)
        return Point(Rx(f),f.A)
    end
end

function readcloud(file)
    coordlines = prepfile(file)
    points = Points("T")
    for line in coordlines
        p = Point(line[2:4],"T")
        push!(points,p)
    end
    println("Read $(length(points)) points into cloud with reference T")
    points
end

function readplane(file)
    coordlines = prepfile(file)
    T = readTmatrix(coordlines)
    n = T[1:3,3]
    o = T[1:3,4]
    r = (n⋅o)*n
    Plane(n,r,"T")
end

function readTmatrix(coordlines)
    T = Matrix(I4)
    for (i,line) in enumerate(coordlines)
        T[i,:] = line
    end
    T
end

function prepfile(file)
    f = open(file)
    text = readlines(f)
    close(f)
    coordlines = Array[]
    for line in text
        stripline = strip(line,[' ','[',']','"',';','\r','\n',' '])
        coords = split(stripline,[',',' '])
        coords = filter(x-> x!="",coords)
        try
            fcoords = parse.(Float64, coords)
            isempty(fcoords) && continue
            push!(coordlines,fcoords)
        catch e
            break
        end
    end
    coordlines
end

end

max_forks_repo_name: baggepinnen/Robotlib


text:
using DataFrames
using ExoplanetsSysSim

function make_actual_stellar_catalog()
  stellar_catalog_file_in = ascii(joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "q1_q17_dr25_stellar.csv"))
  df =  readtable(stellar_catalog_file_in)
  has_mass = ! (isna(df[:mass]) | isna(df[:mass_err1]) | isna(df[:mass_err2]))
  has_radius = ! (isna(df[:radius]) | isna(df[:radius_err1]) | isna(df[:radius_err2]))
  has_teff = ! (isna(df[:teff]) | isna(df[:teff_err1]) | isna(df[:teff_err2]))
  has_logg = ! (isna(df[:logg]) | isna(df[:logg_err1]) | isna(df[:logg_err2]))
  # TODO:  Add additional filters to match stellar catalog of Christian et al
  is_usable = has_mass & has_radius & has_teff & has_logg
  symbols_to_keep = [ :kepid, :mass, :radius, :teff, :logg ]
  df = df[find(is_usable), symbols_to_keep]
  return df
end

# df_star is assumed to have fields kepid, mass and radius for all targets in the survey
function make_actual_planet_candidate_catalog(df_star::DataFrame)
  koi_catalog_file_in = ascii(joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "q1_q17_dr25_koi.csv"))
  (csv_data,csv_header) =  readcsv(koi_catalog_file_in,header=true)
  # Lookup header columns, since DataFrames doesn't like this file
  kepid_idx = findfirst(x->x=="kepid",csv_header)
  koi_period_idx = findfirst(x->x=="koi_period",csv_header)
  koi_time0bk_idx = findfirst(x->x=="koi_time0bk",csv_header)
  koi_depth_idx = findfirst(x->x=="koi_depth",csv_header)
  koi_duration_idx = findfirst(x->x=="koi_duration",csv_header)
  koi_period_err1_idx = findfirst(x->x=="koi_period_err1",csv_header)
  koi_time0bk_err1_idx = findfirst(x->x=="koi_time0bk_err1",csv_header)
  koi_depth_err1_idx = findfirst(x->x=="koi_depth_err1",csv_header)
  koi_duration_err1_idx = findfirst(x->x=="koi_duration_err1",csv_header)
  koi_period_err2_idx = findfirst(x->x=="koi_period_err2",csv_header)
  koi_time0bk_err2_idx = findfirst(x->x=="koi_time0bk_err2",csv_header)
  koi_depth_err2_idx = findfirst(x->x=="koi_depth_err2",csv_header)
  koi_duration_err2_idx = findfirst(x->x=="koi_duration_err2",csv_header)
  koi_disposition_idx = findfirst(x->x=="koi_disposition",csv_header)
  # Choose which KOIs to keep
  idx_keep = (csv_data[:,koi_disposition_idx] .== "CANDIDATE") | (csv_data[:,koi_disposition_idx] .== "CONFIRMED" ) # TODO:  Before publication, check that this is the best disposition field to use
  # TODO:  Add additional filters to match what planet candidates are to be kept Christian et al
  csv_data = csv_data[idx_keep,:]

  output = KeplerObsCatalog([])

  for (j,kepid) in enumerate(df_star[:kepid])
     plids = find(x->x==kepid,csv_data[:,kepid_idx])
	 num_pl = length(plids)
	 if num_pl ==0 continue end
	 target_obs = KeplerTargetObs(num_pl)
	 target_obs.star = ExoplanetsSysSim.StarObs(df_star[j,:radius],df_star[j,:mass])
	 for (plid,i) in enumerate(plids)
	   target_obs.obs[plid] = ExoplanetsSysSim.TransitPlanetObs(csv_data[i,koi_period_idx],csv_data[i,koi_time0bk_idx],csv_data[i,koi_depth_idx],csv_data[i,koi_duration_idx])
	   target_obs.sigma[plid] = ExoplanetsSysSim.TransitPlanetObs((abs(csv_data[i,koi_period_err1_idx])+abs(csv_data[i,koi_period_err2_idx]))/2,(abs(csv_data[i,koi_time0bk_err1_idx])+abs(csv_data[i,koi_time0bk_err2_idx]))/2,(abs(csv_data[i,koi_depth_err1_idx])+abs(csv_data[i,koi_depth_err2_idx]))/2,(abs(csv_data[i,koi_duration_err1_idx])+abs(csv_data[i,koi_duration_err2_idx]))/2)
	   target_obs.prob_detect = ExoplanetsSysSim.OneObserverSystemDetectionProbs( ones(num_pl), ones(num_pl,num_pl), ones(num_pl), fill(Array(Int64,0), 1) )
# ExoplanetsSysSim.OneObserverSystemDetectionProbs(num_pl)
	   push!(output.target,target_obs)
	 end
  end
  return output
end

# Demo of calling the two functions
#df_star = make_actual_stellar_catalog()
#planet_candidate_catalog = make_actual_planet_candidate_catalog(df_star)


text:
export copy_package
"""
    copy_package(model, package)

create a new package based off of model. model must be in the working directory.
"""
function copy_package(model, package)
    if !isdir(model)
        error("cannot find $model")
    elseif isdir(package)
        error("$package already exists")
    else
        @info "copying $model to $package"
        cp(model, package, force = true)
        cd(package)
        @info "removing old git repository"
        rm(".git", recursive = true)
        @info "replacing all occurances of $model with $package"
        foreach(walkdir(pwd())) do dir_tuple
            base_dir, folders, files = dir_tuple
            foreach(files) do file
                new_path = joinpath(base_dir, file)
                write(new_path, join(
                    Generator(readlines(new_path)) do line
                        replace(line, Regex("\\b$model\\b") => package)
                    end,
                    '\n'
                ))
            end
        end
    end

    src_model = joinpath("src", "$model.jl")
    src_package = joinpath("src", "$package.jl")

    if !isfile(src_model)
        @warn "$src_model does not exist"
    else
        @info "renameing $src_model to $src_package"
        mv(src_model, src_package)
    end

    @info "initializing new git repository"
    repo = LibGit2.init(pwd())
    url = "https://github.com/$(settings("username"))/$package.jl.git"
    @info "setting origin to $url"
    LibGit2.set_remote_url(repo, "origin", url)
    @info "initial empty commit"
    LibGit2.commit(repo, "initial empty commit")
    @info "creating gh-pages branch"
    LibGit2.branch!(repo, "gh-pages", set_head=false)
end


text:
"""
    plot_geometry_2D(geometry::Geometry, view::Box, dim::Int64)

Plots a 2D x-y slice of a geometry.

# Arguments
* `geometry::Geometry`: the geometry we want to plot
* `view::Box`: The view box is an axis aligned box that defines where the picture will be taken, with both min and max z dimensions indicating the single z elevation the slice is taken at.
* `dim::Int64`: The dimension is the number of pixels along the x and y axis to use, which determines the resolution of the picture.
"""
function plot_geometry_2D(geometry::Geometry, view::Box, dim::Int64)
    delta_x = (view.upper_right.x - view.lower_left.x) / (dim)
    delta_y = (view.upper_right.y - view.lower_left.y) / (dim)

    x_coords = collect(view.lower_left.x + delta_x/2.0:delta_x:view.upper_right.x - delta_x/2.0)
    y_coords = collect(view.lower_left.y + delta_y/2.0:delta_y:view.upper_right.y - delta_y/2.0)

    pixels = Array{Int64}(undef, dim, dim)

    for i=1:dim
        for j=1:dim
            pixels[i,j] = find_cell_id(Coord(x_coords[i], y_coords[j], view.lower_left.z), geometry)
        end
    end
    pixels = rotl90(pixels)
    colors = Array{RGBA}(undef,0)
    for i=1:length(geometry.cells)
        if (i-1)%4 == 0
            push!(colors, RGBA(1.0, 0.0, 0.0, 1.0))
        elseif (i-1)%4 == 1
            push!(colors, RGBA(0, 1.0, 0, 1.0))
        elseif (i-1)%4 == 2
            push!(colors, RGBA(1.0, 0.0, 1.0, 1.0))
        elseif (i-1)%4 == 3
            push!(colors, RGBA(0, 0, 1.0, 1.0))
        end
        #push!(colors, RGBA(rand(),rand(),rand(),1.0) )
    end
    gradient = cgrad(colors)
    heatmap(x_coords,y_coords,pixels,aspect_ratio=1, color=gradient, leg=false)
end

"""
    plot_cell_2D(geometry::Geometry, view::Box, dim::Int64, cell_id::Int64)

Plots a 2D x-y slice of a geometry, highlighting a specific cell in black.

# Arguments
* `geometry::Geometry`: the geometry we want to plot
* `view::Box`: The view box is an axis aligned box that defines where the picture will be taken, with both min and max z dimensions indicating the single z elevation the slice is taken at.
* `dim::Int64`: The dimension is the number of pixels along the x and y axis to use, which determines the resolution of the picture.
* `cell_id::Int64`: The index of the cell we wish to view
"""
function plot_cell_2D(geometry::Geometry, view::Box, dim::Int64, cell_id::Int64)
    delta_x = (view.upper_right.x - view.lower_left.x) / (dim)
    delta_y = (view.upper_right.y - view.lower_left.y) / (dim)

    x_coords = collect(view.lower_left.x + delta_x/2.0:delta_x:view.upper_right.x - delta_x/2.0)
    y_coords = collect(view.lower_left.y + delta_y/2.0:delta_y:view.upper_right.y - delta_y/2.0)

    pixels = Array{Int64}(undef, dim, dim)

    for i=1:dim
        for j=1:dim
            pixels[i,j] = find_cell_id(Coord(x_coords[i], y_coords[j], view.lower_left.z), geometry)
            if pixels[i, j] == cell_id
                pixels[i, j] = 0
            else
                pixels[i, j] = 1
            end
        end
    end
    pixels = rotl90(pixels)
    colors = Array{RGBA}(undef,0)
    push!(colors, RGBA(0.0, 0.0, 0.0, 1.0))
    push!(colors, RGBA(1.0, 1.0, 1.0, 1.0))
    gradient = cgrad(colors)
    heatmap(x_coords, y_coords, pixels, aspect_ratio=1, color=gradient, leg=false)
end


text:
include("../src/gmm.jl")

# 2D Gaussian mixed model on random data with heatmap.

using Plots

K = 3 # cluster numbers
num_samples = 50
ps, gmm_model, gm_model = create_gmm(K, 2) # prepare model

# prepare data
tmp=ones(2, num_samples)
tmp1=ones(2, num_samples)
tmp1[1, :] .= 0
X = hcat(randn(2, num_samples)/4+tmp1*1.6, randn(2, num_samples)/2, randn(2, num_samples)/4+tmp*1.6)

# learn model params
EM!(ps, X, K, gmm_model, gm_model, 100)

# plot results
x = y = LinRange(minimum(X), maximum(X), 100)
z = Float64[gmm_model(ps, [yi, xi]) for xi = x, yi = y]
p = heatmap(x, y, z, colorbar_title="Likelihood", title = "Gaussian mixture model heatmap")
scatter!(p, X[1, :], X[2, :], label="Data", xlabel="x1", ylabel="x2")

max_stars_repo_name: kozvojtex/AnomalyDetection

max_issues_repo_name: kozvojtex/AnomalyDetection

max_forks_repo_name: kozvojtex/AnomalyDetection


text:
function var(p::AbstractPolynomialLike)
    vars = variables(p)
    if length(vars) == 1
        return vars[1]
    elseif length(vars) == 0
        return nothing
    else
        error("Polynomial should have only one variable")
    end
end

function var(p::AbstractMonomialLike)
    vars = variables(p)
    if length(vars) == 1
        return vars[1]
    elseif length(vars) == 0
        return nothing
    else
        error("Polynomial should have only one variable")
    end
end

var(p::AbstractTermLike) = variables(p)[1]

function var(p)
    vars = get_variables(p)
    if length(vars) == 1
        return vars[1]
    elseif length(vars) == 0
        return nothing
    else
        error("Polynomial should have only one variable")
    end
end

wrap(fun, p::AbstractPolynomialLike) = fun(p)
wrap(fun, p::AbstractPolynomialLike, q::AbstractPolynomialLike) = fun(p, q)

function prewrap(eq)
    x = var(eq)
    x == nothing && return eq * one(𝑦), 𝑦 => x
    p = poly(expand(eq), x => 𝑦)
    p, 𝑦 => x
end

function wrap(fun, eq)
    p, v = prewrap(eq)
    q = fun(p)
    unwrap(q, v)
end

function wrap(fun, eq₁, eq₂)
    x₁ = var(eq₁)
    x₂ = var(eq₂)

    if x₁ == nothing && x₂ == nothing
        return fun(eq₁ * one(𝑦), eq₂ * one(𝑦))
    elseif x₁ == nothing
        p = eq₁ * one(𝑦)
        q = poly(eq₂, x₂ => 𝑦)
        x = x₂
    elseif x₂ == nothing
        p = poly(eq₁, x₁ => 𝑦)
        q = eq₂ * one(𝑦)
        x = x₁
    else
        !isequal(x₁, x₂) && error("incompatible main variables")
        p = poly(eq₁, x₁ => 𝑦)
        q = poly(eq₂, x₂ => 𝑦)
        x = x₁
    end

    unwrap(fun(p, q), 𝑦 => x)
end

unwrap(q::AbstractPolynomial, v) = sym(q, v)
unwrap(q::RationalPoly, v) = sym(numerator(q), 𝑦 => x) / sym(denominator(q), 𝑦 => x)
unwrap(q::FactoredPoly, v) = sym(q, v)
unwrap(q::Tuple, v) = map(x->unwrap(x,v), q)
unwrap(q, v) = q

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

rationalize(p) = wrap(x->polynomial(convert.(Rational{BigInt}, coefficients(x)), terms(x)), p)

function unrationalize(p::AbstractPolynomial)
    t = terms(p)
    c = map(x -> x isa Rational && denominator(x)==1 ? numerator(x) : x, coefficients(p))
    polynomial(c, t)
end

unrationalize(p) = wrap(unrationalize, p)

function gcd_extended(u::AbstractPolynomialLike, v::AbstractPolynomialLike)
    # u = rationalize(u)
    # v = rationalize(v)

    !isequal(var(u), var(v)) && error("incompatible main variable")

    if deg(u) == 0 && deg(v) == 0
        return gcdx(cont(u), cont(v))
    end

    sᵤ = 1
    tᵤ = 0
    sᵥ = 0
    tᵥ = 1

    while !iszero(v)
        q, r = divrem(u, v)
        s, t = sᵤ - q*sᵥ, tᵤ - q*tᵥ
        u, sᵤ, tᵤ = v, sᵥ, tᵥ
        v, sᵥ, tᵥ = r, s, t
    end

    l = leading(u)
    return u / l, sᵤ / l, tᵤ / l
end

Base.gcdx(u, v) = wrap(gcd_extended, u, v)

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

leading(p::AbstractPolynomialLike) = leadingcoefficient(p)
cont(p::AbstractPolynomialLike) = iszero(p) ? 0 : gcd(coefficients(p)...) * sign(leading(p))
prim(p::AbstractPolynomialLike) = polynomial(coefficients(p) .÷ cont(p), terms(p))
derivative(p::AbstractPolynomialLike) = differentiate(p, var(p))
deg(p::Number) = 0
deg(p::AbstractPolynomialLike, x) = x==nothing ? 0 : maxdegree(p, x)
deg(p::AbstractPolynomialLike) = maxdegree(p, var(p))
coef(p::AbstractPolynomialLike) = lcm(denominator.(coefficients(p))...)

# leading(p::AbstractTerm) = leadingcoefficient(p)
# cont(p::AbstractTerm) = gcd(coefficients(p)...) * sign(leading(p))
# prim(p::AbstractTerm) = p ÷ cont(p)

leading(p) = wrap(leading, p)
cont(p) = wrap(cont, p)
prim(p) = wrap(prim, p)
derivative(p) = wrap(derivative, p)
deg(p) = wrap(deg, p)
coef(p) = wrap(coef, p)

function Base.:÷(p::AbstractPolynomialLike, k)
    polynomial(coefficients(p) .÷ k, terms(p))
end

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

function integer_poly_coef(p::Polynomial{true,T}) where T<:Rational
    l = coef(p)
    1//l, polynomial(Int.(numerator.(coefficients(l*p))), terms(p))
end

integer_poly_coef(p::Polynomial{true,T}) where T<:Integer = 1, p
integer_poly_coef(p::Polynomial{true,T}) where T<:Real = 1, integer_poly(p)

function integer_poly(p::Polynomial{true,T}) where T<:Rational
    l = coef(p)
    polynomial(Int.(numerator.(coefficients(l*p))), terms(p))
end

function integer_poly(p::AbstractPolynomialLike)
    polynomial(round.(Int, coefficients(p)), terms(p))
end

function bigint_poly(p::AbstractPolynomialLike)
    p = integer_poly(p)
    polynomial(BigInt.(coefficients(p)), terms(p))
end

integer_poly(p::Polynomial{true,T}) where T<:Integer = p

function to_monic(p::Polynomial{true,T}) where T<:Rational
    c = leading(p)
    n = deg(p)
    d = degree.(terms(p))
    polynomial(coefficients(p) .* c.^((n-1).-d), terms(p)), c
end

function to_monic(p::Polynomial{true,T}) where T<:Integer
    q, c = to_monic(rationalize(p))
    integer_poly(q), Int(numerator(c))
end

function from_monic(p::Polynomial{true,T}, c) where T<:Rational
    n = deg(p)
    d = degree.(terms(p))
    polynomial(coefficients(p) .* (1 .// c.^((n-1).-d)), terms(p))
end

function from_monic(p::Polynomial{true,T}, c) where T<:Integer
    q = from_monic(rationalize(p), Rational(c))
    integer_poly(q)
end

from_monic(p, c) = from_monic(polynomial(p), c)

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

function standard_form(p::Polynomial{true,T}) where T<:Integer
    # q, lc = to_monic(prim(p))
    # undo = s -> from_monic(s, lc)
    # q, undo
    prim(p), x->x
end

function standard_form(p::Polynomial{true,T}) where T<:Rational
    q = integer_poly(p)
    standard_form(q)
end

function standard_form(eq)
    p, v = prewrap(eq)
    q, undo = standard_form(p)
    q, x -> unwrap(undo(x), v)
end


text:
"""
`init_subspec!(m::Model805_alt2)`

Initializes a model subspecification by overwriting parameters from
the original model object with new parameter objects. This function is
called from within the model constructor.
"""
function init_subspec!(m::Model805_alt2)

    if subspec(m) == "ss0"
        # Use σ_r_m as SD for all anticipated shocks without normalizing
        # See measurement equation
        return
    elseif subspec(m) == "ss1"
        # Normalize s.t. sum of variances of anticipated shocks equals the
        # variance of the contemporaneous shock
        # See measurement equation
        return
    elseif subspec(m) =="ss3"
        # Diffuse prior
        ss3!(m)
    elseif subspec(m) == "ss4"
        # Diffuse prior but keep pi* the same
        ss3!(m)
        ss4!(m)
    elseif subspec(m) == "ss20"
        ss20!(m)
    else
        error("This subspec is not defined.")
    end

end

# Diffuse prior
function ss3!(m::Model805_alt2)

    m <= parameter(:α,      0.1687, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),     Normal(0.30, 0.15),         fixed=false,
                   description="α: Capital elasticity in the intermediate goods sector's Cobb-Douglas production function.",
                   tex_label="\\alpha")

    m <= parameter(:ζ_p,   0.7467, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),     Distributions.Uniform(0.0, 1.0),          fixed=false,
                   description="ζ_p: The Calvo parameter. In every period, intermediate goods producers optimize prices with probability (1-ζ_p). With probability ζ_p, prices are adjusted according to a weighted average of the previous period's inflation (π_t1) and steady-state inflation (π_star).",
                   tex_label="\\zeta_p")

    m <= parameter(:ι_p,   0.2684, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(), Distributions.Uniform(0.0, 1.0),         fixed=false,
                   description="ι_p: The weight attributed to last period's inflation in price indexation. (1-ι_p) is the weight attributed to steady-state inflation.",

                   tex_label="\\iota_p")

    m <= parameter(:δ,      0.025,  fixed=true,
                   description="δ: The capital depreciation rate.", tex_label="\\delta" )

    m <= parameter(:Upsilon,  1.000,  (0., 10.),     (1e-5, 0.),      ModelConstructors.Exponential(),    GammaAlt(1., 0.5),          fixed=true,
                   description="Υ: The trend evolution of the price of investment goods relative to consumption goods. Set equal to 1.",
                   tex_label="\\mathcal{\\Upsilon}")

    m <= parameter(:Φ,   1.4933, (1., 10.),     (1.00, 10.00),   ModelConstructors.Exponential(),    Normal(1.25, 0.36),         fixed=false,
                   description="Φ: Fixed costs.",
                   tex_label="\\Phi")

    m <= parameter(:S′′,       3.3543, (-15., 15.),   (-15., 15.),     Untransformed(),  Normal(4., 4.5),            fixed=false,
                   description="S'': The second derivative of households' cost of adjusting investment.", tex_label="S\\prime\\prime")

    m <= parameter(:h,        0.4656, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(), Distributions.Uniform(0.0, 1.0),          fixed=false,
                   description="h: Consumption habit persistence.", tex_label="h")

    m <= parameter(:ppsi,     0.7614, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),  Distributions.Uniform(0.0, 1.0),         fixed=false,
                   description="ppsi: Utilization costs.", tex_label="ppsi")

    m <= parameter(:ν_l,     1.0647, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    Normal(2, 2.25),            fixed=false,
                   description="ν_l: The coefficient of relative risk aversion on the labor term of households' utility function.", tex_label="\\nu_l")

    m <= parameter(:ζ_w,   0.7922, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(), Distributions.Uniform(0.0, 1.0),          fixed=false,
                   description="ζ_w: (1-ζ_w) is the probability with which households can freely choose wages in each period. With probability ζ_w, wages increase at a geometrically weighted average of the steady state rate of wage increases and last period's productivity times last period's inflation.",
                   tex_label="\\zeta_w")

    m <= parameter(:ι_w,   0.5729, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(), Distributions.Uniform(0.0, 1.0),         fixed=false,
                   description="ι_w: No description available.",
                   tex_label="\\iota_w")

    m <= parameter(:λ_w,      1.5000,                                                                               fixed=true,
                   description="λ_w: The wage markup, which affects the elasticity of substitution between differentiated labor services.",
                   tex_label="\\lambda_w")

    m <= parameter(:β, 0.7420, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    GammaAlt(0.25, 0.3),        fixed=false,  scaling = x -> 1/(1 + x/100),
                   description="β: Discount rate.",
                   tex_label="\\beta ")

    m <= parameter(:ψ1,  1.8678, (1e-5, 10.),   (1e-5, 10.00),   ModelConstructors.Exponential(),    Normal(1.5, 0.75),          fixed=false,
                   description="ψ₁: Weight on inflation gap in monetary policy rule.",
                   tex_label="\\psi_1")

    m <= parameter(:ψ2,  0.0715, (-0.5, 0.5),   (-0.5, 0.5),     Untransformed(),  Normal(0.12, 0.15),         fixed=false,
                   description="ψ₂: Weight on output gap in monetary policy rule.",
                   tex_label="\\psi_2")

    m <= parameter(:ψ3, 0.2131, (-0.5, 0.5),   (-0.5, 0.5),     Untransformed(),  Normal(0.12, 0.15),         fixed=false,
                   description="ψ₃: Weight on rate of change of output gap in the monetary policy rule.",
                   tex_label="\\psi_3")
 ## The old prior for π* in SW is GammaAlt(0.62, 0.1) with diffuse GammaAlt(0.62, 0.3) and the old prior here was GammaAlt(0.75, 0.4) so we already had a higher variance, but in keeping with the 'multiply by 3' rule, I'm multiplying by 3 anyways.
    m <= parameter(:π_star,   0.6231, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    GammaAlt(0.75, 1.2),        fixed=false,  scaling = x -> 1 + x/100,
                   description="π_star: The steady-state rate of inflation.",
                   tex_label="\\pi_*")

m <= parameter(:σ_c, 1.5073, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    Normal(1.5, 1.11),          fixed=false,
               description="σ_c: No description available.",
               tex_label="\\sigma_{c}")

m <= parameter(:ρ,      .8519, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(), Distributions.Uniform(0.0, 1.0),        fixed=false,
               description="ρ: The degree of inertia in the monetary policy rule.",
               tex_label="\\rho")

m <= parameter(:ϵ_p,     10.000,                                                                               fixed=true,
               description="ϵ_p: No description available.",
               tex_label="\\varepsilon_{p}")

m <= parameter(:ϵ_w,     10.000,                                                                               fixed=true,
               description="ϵ_w: No description available.",
               tex_label="\\varepsilon_{w}")

# exogenous processes - level
m <= parameter(:γ,      0.3085, (-5.0, 5.0),     (-5., 5.),     Untransformed(), Normal(0.4, 0.3),            fixed=false, scaling = x -> x/100,
               description="γ: The log of the steady-state growth rate of technology.",
               tex_label="\\gamma")

m <= parameter(:Lmean,  -45., (-1000., 1000.), (-1e3, 1e3),   Untransformed(), Normal(-45, 15),   fixed=false,
               description="Lmean: No description available.",
               tex_label="Lmean")

m <= parameter(:g_star,    0.1800,                                                                               fixed=true,
               description="g_star: No description available.",
               tex_label="g_*")

# exogenous processes - autocorrelation
m <= parameter(:ρ_g,      0.9930, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
               description="ρ_g: AR(1) coefficient in the government spending process.",
               tex_label="\\rho_g")

m <= parameter(:ρ_b,      0.8100, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
               description="ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
               tex_label="\\rho_b")

m <= parameter(:ρ_μ,     0.7790, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
               description="ρ_μ: AR(1) coefficient in capital adjustment cost process.",
               tex_label="\\rho_{\\mu}")

m <= parameter(:ρ_z,      0.9676, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
               description="ρ_z: AR(1) coefficient in the technology process.",
               tex_label="\\rho_z")

m <= parameter(:ρ_λ_f,    0.8372, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
               description="ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
               tex_label="\\rho_{\\lambda_f}")

m <= parameter(:ρ_λ_w,    0.9853, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
               description="ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.", # CHECK THIS
               tex_label="\\rho_{\\lambda_w}")

m <= parameter(:ρ_rm,     0.9000, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
               description="ρ_rm: AR(1) coefficient in the monetary policy shock process.", # CHECK THIS
               tex_label="\\rho_{rm}")

m <= parameter(:ρ_π_star,     0.9900, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),
               fixed=true,
               description="ρ_π_star: AR(1) coefficient in the π* process.",
               tex_label="\\rho_{\\pi^*}")
# exogenous processes - standard deviation
m <= parameter(:σ_g,      0.6090, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
               description="σ_g: The standard deviation of the government spending process.",
               tex_label="\\sigma_{g}")

m <= parameter(:σ_b,      0.1818, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
               description="σ_b: No description available.",
               tex_label="\\sigma_{b}")

m <= parameter(:σ_μ,     0.4601, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
               description="σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
               tex_label="\\sigma_{\\mu}")

m <= parameter(:σ_z,      0.4618, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
               description="σ_z: No description available.",
               tex_label="\\sigma_{z}")

m <= parameter(:σ_λ_f,    0.1455, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
               description="σ_λ_f: The mean of the process that generates the price elasticity of the composite good.  Specifically, the elasticity is (1+λ_{f,t})/(λ_{f_t}).",
               tex_label="\\sigma_{\\lambda_f}")

m <= parameter(:σ_λ_w,    0.2089, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
               description="σ_λ_w: No description available.",
               tex_label="\\sigma_{\\lambda_w}")

m <= parameter(:σ_rm,     0.2397, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
               description="σ_rm: No description available.",
               tex_label="\\sigma_{rm}")

m <= parameter(:σ_π_star,      0.0300, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
               description="σ_π_star: The standard deviation of the π* process.",
               tex_label="\\sigma_{\\pi^*}")

m <= parameter(:η_gz,       0.5632, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
               Distributions.Uniform(0.0, 1.0), fixed=false,
               description="η_gz: Correlate g and z shocks.",
               tex_label="\\eta_{gz}")

m <= parameter(:η_λ_f,      0.7652, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
               Distributions.Uniform(0.0, 1.0),         fixed=false,
               description="η_λ_f: No description available.",
               tex_label="\\eta_{\\lambda_f}")

m <= parameter(:η_λ_w,      0.8936, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
               Distributions.Uniform(0.0, 1.0),         fixed=false,
               description="η_λ_w: AR(2) coefficient on wage markup shock process.",
               tex_label="\\eta_{\\lambda_w}")

end

function ss4!(m::Model805_alt2)
    m <= parameter(:π_star,   0.6231, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    GammaAlt(0.75, 0.4),        fixed=false,  scaling = x -> 1 + x/100,
                   description="π_star: The steady-state rate of inflation.",
                   tex_label="\\pi_*")
end

function ss20!(m::Model805_alt2)
    # same specification use for Model1010depo_alt

    m <= parameter(:σ_b_p, sqrt(1/400)/4, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=false,
                   description="σ_b_p: Standard deviation of stationary component of asset preference shifter process.",
                   tex_label="\\sigma_{b^p}")

    m <= parameter(:ρ_z_p,     0.99, (0.0, 1.0),      (0.0, 1.0), ModelConstructors.SquareRoot(),    BetaAlt(0.5, 0.2),
                   fixed=true, description="ρ_z_p: No description available.", tex_label="\\rho_{z^p}")
end


text:
doc"""
Returns a list of pairs (interval, label)
label is 1 if the interval is *excluded* from the setdiff
label is 0 if the interval is included in the setdiff
label is -1 if the intersection of the two intervals was empty
"""
function labelled_setdiff{T}(x::Interval{T}, y::Interval{T})
    intersection = x ∩ y

    isempty(intersection) && return [(x, -1)]
    intersection == x && return [(x, 1)]

    x.lo == intersection.lo && return [(intersection, 1), (Interval(intersection.hi, x.hi), 0)]
    x.hi == intersection.hi && return [(intersection, 1), (Interval(x.lo, intersection.lo), 0)]

    return [(y, 1),
            (Interval(x.lo, y.lo), 0),
            (Interval(y.hi, x.hi), 0)]

end

# function setdiff{N,T}(A::IntervalBox{N,T}, B::IntervalBox{N,T})
#     X = [labelled_setdiff(a,b) for (a, b) in zip(A, B)]
#     lengths = map(length, X)
#     index = ones(N)
#
#     # index[j] represents which set we are looking at in direction j
#``
#
#     while index[1] <= N
#         current_direction = 1
#         current_piece = [ X[1][index[1]] ]
#
#     end
# end
doc"""
    setdiff(A::IntervalBox{N,T}, B::IntervalBox{N,T})

Returns a vector of `IntervalBox`es that are in the set difference `A \ B`,
i.e. the set of `x` that are in `A` but not in `B`.

Algorithm: Start from the total overlap (in all directions);
expand each direction in turn.
"""
function setdiff{N,T}(A::IntervalBox{N,T}, B::IntervalBox{N,T})
    X = [labelled_setdiff(a,b) for (a, b) in zip(A, B)]
    # ordered such that the first in each is the excluded interval

    first = [ i[1] for i in X ]
    labels = [i[2] for i in first]

    any(labels .== -1) && return [A]  # no overlap

    # @assert all(labels .== 1)

    excluded = [i[1] for i in first]

    result_list = IntervalBox{N,T}[]

    for dimension in N:-1:1
        for which in X[dimension][2:end]
            excluded[dimension] = which[1]
            push!(result_list,
                IntervalBox(excluded[1:dimension]..., A[dimension+1:end]...))
        end
    end

    result_list

end

# doc"""
#     setdiff(A::IntervalBox{2,T}, B::IntervalBox{2,T})
#
# Returns a vector of `IntervalBox`es that are in the set difference `A \ B`,
# i.e. the set of `x` that are in `A` but not in `B`.
# """
# function setdiff{T}(A::IntervalBox{2,T}, B::IntervalBox{2,T})
#     X = labelled_setdiff(A[1], B[1])
#     Y = labelled_setdiff(A[2], B[2])
#
#     results_list = typeof(A)[]
#
#     for (x, label) in X
#         label == -1 && return [A]   # intersection in one direction empty, so total intersection empty
#
#         if label == 0
#             push!(results_list, IntervalBox(x, A[2]))
#             continue
#         end
#
#         # label is 1 here, so there is some intersection in the x direction
#         for (y, label) in Y
#             label == -1 && return [A]
#
#             if label == 0
#                 push!(results_list, IntervalBox(x, y))
#                 continue
#             end
#
#             # label == 1:  exclude this box since all labels are 1
#         end
#     end
#
#     return results_list
# end


text:
"""
    RectilinearMesh <: AbstractMesh2D

# Constructors

    RectilinearMesh(domain::RectangleDomain, n1, n2)

Construct a `RectilinearMesh` of equally-spaced `n1` by `n2` elements on `domain`.

    RectilinearMesh(intervalmesh1::IntervalMesh1, intervalmesh2::IntervalMesh2)

Construct the product mesh of `intervalmesh1` and `intervalmesh2`.
"""
struct RectilinearMesh{I1 <: IntervalMesh, I2 <: IntervalMesh} <: AbstractMesh2D
    intervalmesh1::I1
    intervalmesh2::I2
end
RectilinearMesh(domain::RectangleDomain, n1, n2) = RectilinearMesh(
    IntervalMesh(domain.interval1; nelems = n1),
    IntervalMesh(domain.interval2; nelems = n2),
)

function Base.show(io::IO, mesh::RectilinearMesh)
    n1 = nelements(mesh.intervalmesh1)
    n2 = nelements(mesh.intervalmesh2)
    print(io, n1, "×", n2, "-element RectilinearMesh of ", domain(mesh))
end
domain(mesh::RectilinearMesh) =
    RectangleDomain(domain(mesh.intervalmesh1), domain(mesh.intervalmesh2))
nelements(mesh::RectilinearMesh) =
    nelements(mesh.intervalmesh1) * nelements(mesh.intervalmesh2)
function elements(mesh::RectilinearMesh)
    # we use the Base Julia CartesianIndices object to index elements in the mesh
    CartesianIndices((
        elements(mesh.intervalmesh1),
        elements(mesh.intervalmesh2),
    ))
end
function is_boundary_face(mesh::RectilinearMesh, elem::CartesianIndex{2}, face)
    @assert 1 <= face <= 4
    x1, x2 = elem.I
    if face == 1
        return is_boundary_face(mesh.intervalmesh2, x2, 1)
    elseif face == 2
        return is_boundary_face(mesh.intervalmesh1, x1, 2)
    elseif face == 3
        return is_boundary_face(mesh.intervalmesh2, x2, 2)
    else
        return is_boundary_face(mesh.intervalmesh1, x1, 1)
    end
end
function boundary_face_name(
    mesh::RectilinearMesh,
    elem::CartesianIndex{2},
    face,
)
    @assert 1 <= face <= 4
    x1, x2 = elem.I
    if face == 1
        return boundary_face_name(mesh.intervalmesh2, x2, 1)
    elseif face == 2
        return boundary_face_name(mesh.intervalmesh1, x1, 2)
    elseif face == 3
        return boundary_face_name(mesh.intervalmesh2, x2, 2)
    else
        return boundary_face_name(mesh.intervalmesh1, x1, 1)
    end
end

function opposing_face(
    mesh::RectilinearMesh,
    elem::CartesianIndex{2},
    face::Int,
)
    @assert 1 <= face <= 4
    x1, x2 = elem.I
    n1, n2 = size(elements(mesh))
    if face == 1
        return CartesianIndex(x1, mod1(x2 - 1, n2)), 3, true
    elseif face == 2
        return CartesianIndex(mod1(x1 + 1, n1), x2), 4, true
    elseif face == 3
        return CartesianIndex(x1, mod1(x2 + 1, n2)), 1, true
    else
        return CartesianIndex(mod1(x1 - 1, n1), x2), 2, true
    end
end

function coordinates(mesh::RectilinearMesh, elem::CartesianIndex{2}, vert::Int)
    x1, x2 = elem.I
    coord1 = coordinates(mesh.intervalmesh1, x1, vert == 1 || vert == 4 ? 1 : 2)
    coord2 = coordinates(mesh.intervalmesh2, x2, vert == 1 || vert == 2 ? 1 : 2)
    return Geometry.product_coordinates(coord1, coord2)
end
function coordinates(
    mesh::RectilinearMesh,
    elem::CartesianIndex{2},
    (ξ1, ξ2)::NTuple{2},
)
    x1, x2 = elem.I
    coord1 = coordinates(mesh.intervalmesh1, x1, (ξ1,))
    coord2 = coordinates(mesh.intervalmesh2, x2, (ξ2,))
    return Geometry.product_coordinates(coord1, coord2)
end

function containing_element(mesh::RectilinearMesh, coord)
    x1, (ξ1,) =
        containing_element(mesh.intervalmesh1, Geometry.coordinate(coord, 1))
    x2, (ξ2,) =
        containing_element(mesh.intervalmesh2, Geometry.coordinate(coord, 2))
    return CartesianIndex(x1, x2), (ξ1, ξ2)
end


text:
using FixedPointNumbers

include("microcode.jl")

const Fix = Fixed{Int32, 16}
const Fix_frac_bits = 16

int2fix(i::Integer) = Fix(i, 0)
fix2int(x::Fix) = x.i

struct Parameters
	mem_int_bits::Int
	reg_int_bits::Int
	result_int_bits::Int
	frac_bits::Int
end

mutable struct DataPath
	p::Parameters

	result::Fix
	t::Fix
	f1::Fix
	f2::Fix

	positive::Bool # flag stored by some instructions
	overflow::Bool # flag stored by some instructions
	bias::Int # used when C_DITHER is set

	outputs::Vector{Fix}
end
DataPath(p::Parameters) = DataPath(p, 0, 0, 0, 0, true, false, 0, zeros(Fix, 2))

function fix_trunc(x::Fix, int_bits::Int, frac_bits::Int)
	shift = 32 - (Fix_frac_bits + int_bits + 1)
	i::Int32 = x.i << shift
	i >>= shift
	i &= (-1 << (Fix_frac_bits - frac_bits))
	return int2fix(i)
end

mem_trunc(p::Parameters, x::Fix) = fix_trunc(x, p.mem_int_bits, p.frac_bits)
reg_trunc(p::Parameters, x::Fix) = fix_trunc(x, p.reg_int_bits, p.frac_bits)

function get_a(dp::DataPath)
	@assert dp.t == dp.f1 == dp.f2
	return mem_trunc(dp.p, dp.t)
end

fix_round(x, frac_bits) = int2fix(round(Int, x*(1 << frac_bits)) << (Fix_frac_bits - frac_bits))

function madd(p::Parameters, f1::Fix, f2::Fix, t::Fix, bias::Int)
	f1, f2, t = reg_trunc(p, f1), reg_trunc(p, f2), reg_trunc(p, t)
	bias &= (1 << (2*Fix_frac_bits - p.frac_bits)) - 1

	f1i::Int64 = f1.i
	f2i::Int64 = f2.i
	ti::Int64 = t.i

	ti_full = (ti << Fix_frac_bits) + bias
	r = ((f1i*f2i + ti_full) >> Fix_frac_bits)

	return fix_trunc(int2fix(r), p.result_int_bits, p.frac_bits)
end
madd(p::Parameters, f1::Fix, f2::Fix, t::Fix) = madd(p, f1, f2, t, 1 << (2*Fix_frac_bits - (p.frac_bits + 1)))

Base.clamp(p::Parameters, x::Fix) = mem_trunc(p, int2fix(clamp(x.i, -1 << (Fix_frac_bits + p.mem_int_bits), (1 << (Fix_frac_bits + p.mem_int_bits)) - 1)))

# Functions to set result
mem!(  dp::DataPath, x::Fix)       = (dp.result = mem_trunc(dp.p, x))
const!(dp::DataPath, x::Fix)       = (dp.result = reg_trunc(dp.p, x))
madd!( dp::DataPath, dither::Bool) = (dp.result = madd(dp.p, dp.f1, dp.f2, dp.t, dither ? dp.bias : 1 << (2*Fix_frac_bits - (dp.p.frac_bits + 1)) ))
clamp!(dp::DataPath)               = (dp.result = clamp(dp.p, dp.result))
trunc!(dp::DataPath)               = (dp.result = mem_trunc(dp.p, dp.result))

function bit_reverse(x::Integer, num_bits::Integer)
	y = 0
	for k=1:num_bits
		y <<= 1
		y |= x & 1
		x >>= 1
	end
	return y
end

function reg_write!(dp::DataPath, we::Int)
	if (we & SIGN_WE != 0) dp.positive = (dp.f2 >= 0); end
	r = reg_trunc(dp.p, dp.result)
	if (we & T_WE  != 0) dp.t = r; end
	if (we & F1_WE != 0) dp.f1 = r; end
	if (we & F2_WE != 0) dp.f2 = r; end

	if (we & BIAS_WE != 0)
		delta_frac_bits = Fix_frac_bits - dp.p.frac_bits
		bits = fix2int(r) >> delta_frac_bits
		dp.bias = bit_reverse(bits, dp.p.frac_bits) << (2*delta_frac_bits)
	end
end

lfsr16_bit(x::Integer) = xor(x >> 15, x >> 14, x >> 12, x >> 3) & 1

function exec!(dp::DataPath, context::IContext, c::DPControl, mem::Fix)
	# Make sure that these flags are stored in the c.res field for now
	@assert c.reg_we & C_CLAMP == 0

	control = c.res | c.reg_we
	res = control & ((1 << result_sel_bits)-1)

	# Based on original result selector
	result_sel = res < 4 ? res : 0

	const_sel = res
	# TODO: handle res_exp2 here instead
	if control & C_STORED_SIGN != 0
		const_sel = (const_sel & ~1) | (dp.positive ? 0 : 1)
	elseif res == res_lfsr
		x = fix2int(dp.f2) >> (Fix_frac_bits - p.frac_bits)
		const_sel = (const_sel & ~1) | (lfsr16_bit(x) & 1)
	end

	dither = (control & C_DITHER != 0)

	store_overflow = (control & OVERFLOW_WE != 0)
	@assert !(store_overflow && result_sel != res_madd)

	if result_sel == res_mem
		mem!(dp, mem)
	elseif result_sel == res_madd
		madd!(dp, dither)
		if store_overflow; dp.overflow = !(Fix(-1 << p.mem_int_bits) <= dp.result < Fix(1 << p.mem_int_bits)); end
		if (control & C_CLAMP != 0); clamp!(dp); end
	elseif result_sel == res_const
		if     const_sel == res_zero;       const!(dp, Fix(0))
		elseif const_sel == res_maxval;     const!(dp, int2fix((1 << Fix_frac_bits)-1))
		elseif const_sel == res_madd;       const!(dp, Fix(-1)) # paired with res_maxval for clamping and pulse wave
		elseif const_sel == res_one;        const!(dp, Fix(1))
		elseif const_sel == res_m_one;      const!(dp, Fix(-1))
		elseif const_sel == res_two;        const!(dp, Fix(2))
		elseif const_sel == res_m_two;      const!(dp, Fix(-2))
		elseif const_sel == res_lfsr;       const!(dp, Fix(0))
		elseif const_sel == res_eps;        const!(dp, int2fix(1<<(Fix_frac_bits - dp.p.frac_bits)))
		elseif const_sel == res_1_fifth;    const!(dp, fix_round(1/5, dp.p.frac_bits))
		elseif const_sel == res_m_6_fifths; const!(dp, fix_round(-6/5, dp.p.frac_bits))
		elseif const_sel == res_exp2;       const!(dp, context.n <= IC_N_ZERO ? Fix(0) : int2fix(1<<(Fix_frac_bits + context.n)))
		else error(string("Unsupported const_sel = ", const_sel))
		end
	else
		error(string("Unsupported result_sel = ", result_sel))
	end

	if (control & C_TRUNC != 0) trunc!(dp); end

	reg_write!(dp, c.reg_we)

	r = mem_trunc(dp.p, dp.result)
	if (c.reg_we & OUTPUT_WE != 0) dp.outputs[context.n+1] = r; end
	return r
end
exec!(dp::DataPath, c::DPControl) = exec!(dp, c, Fix(0))

function exec!(dp::DataPath, context::IContext, c::CycleControl, imm::Fix, mem::Fix)
	mem_data = c.cn & C_MEM_ADDR != 0 ? mem : imm;
	r = exec!(dp, context, c.cdp, mem_data)
	mem_write_mask = 0
	if (c.cn & C_M_WE != 0) && ((c.cn & C_CONDITIONAL_M_WE == 0) || dp.overflow)
		if (c.cn & C_MEM_ADDR != 0)
			mem_write_mask = 2
			mem = r
		else
			mem_write_mask = 1
			imm = r
		end
	end
	return (imm, mem, mem_write_mask)
end

function exec!(dp::DataPath, context::IContext, cs::Vector{CycleControl}, imm::Fix, mem::Fix)
	for c in cs
		(imm, mem, mem_write_mask) = exec!(dp, context, c, imm, mem)
	end
	return (imm, mem)
end

exec!(dp::DataPath, context::IContext, c, imm::Fix) = ((imm, mem) = exec!(dp, context, c, imm, Fix(0)); imm)
exec!(dp::DataPath, context::IContext, c) =                        (exec!(dp, context, c, Fix(0), Fix(0)); ())

exec!(dp::DataPath, inst::Instruction, imm::Fix, mem::Fix) = exec!(dp, inst.context, opcode_programs[inst.op], imm, mem)
exec!(dp::DataPath, op::Opcode, imm::Fix, mem::Fix) = exec!(dp, Instruction(op), imm, mem)

exec!(dp::DataPath, i::Union{Instruction, Opcode}, imm::Fix) = ((imm, mem) = exec!(dp, i, imm, Fix(0)); imm)
exec!(dp::DataPath, i::Union{Instruction, Opcode}) =                        (exec!(dp, i, Fix(0), Fix(0)); ())

function sina2(dp::DataPath, x)
	exec!(dp, inst_imm, Fix(x))
	exec!(dp, inst_sina2)
	return get_a(dp)
end

function sina4(dp::DataPath, x)
	exec!(dp, inst_imm, Fix(x))
	exec!(dp, inst_sin)
	return get_a(dp)
end

function eval_inst!(dp::DataPath, inst, a, imm)
	exec!(dp, inst_imm, Fix(a))
	exec!(dp, inst, Fix(imm))
	return get_a(dp)
end
eval_inst!(dp::DataPath, inst, a) = eval_inst!(dp, inst, a, Fix(0))

mutable struct DSP
	dp::DataPath
	code::Vector{Instruction}
	data::Vector{Fix}

	pc::Int
end
DSP(dp, code, data) = DSP(dp, code, data, 0)

function step!(dsp::DSP)
	inst = dsp.code[dsp.pc+1]

	for control in opcode_programs[inst.op]
		# read
		(imm, mem) = (dsp.data[dsp.pc+1], dsp.data[inst.context.addr+1])
		# execute
		(imm, mem, mem_write_mask) = exec!(dsp.dp, inst.context, control, imm, mem)
		# write back
		if     mem_write_mask == 1; dsp.data[dsp.pc+1] = imm;
		elseif mem_write_mask == 2; dsp.data[inst.context.addr+1] = mem;
		end
	end

	dsp.pc += 1
	if dsp.pc == length(dsp.code); dsp.pc = 0; end
end

function run_loop!(dsp::DSP)
	last_pc = -1
	while (dsp.pc > last_pc)
		last_pc = dsp.pc
		step!(dsp)
	end
end

function run_loops!(dsp::DSP, output::Vector{Fix}, n::Int)
	for k=1:n
		run_loop!(dsp)
		push!(output, dsp.dp.outputs[1])
	end
	return output
end
run_loops!(dsp::DSP, n::Int) = run_loops!(dsp::DSP, Fix[], n::Int)

function create_dsp(p::Parameters, program)
	code = Vector{Instruction}(undef, length(program))
	data = zeros(Fix, length(program))
	for (k, (c, d)) in enumerate(program)
		if isa(c, Opcode); c = Instruction(c); end
		code[k] = c;
		data[k] = d;
	end

	return DSP(DataPath(p), code, data)
end

max_stars_repo_name: toivoh/fmicrosynth

max_issues_repo_name: toivoh/fmicrosynth

max_forks_repo_name: toivoh/fmicrosynth


text:
function entropy(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)
    V = volume(model, p, T, z; phase=phase, threaded=threaded)
    return VT_entropy(model,V,T,z)
end

function chemical_potential(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)
    V = volume(model, p, T, z; phase=phase, threaded=threaded)
    return VT_chemical_potential(model,V,T,z)
end

function chemical_potential_res(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)
    V = volume(model, p, T, z; phase=phase, threaded=threaded)
    return VT_chemical_potential_res(model,V,T,z)
end

function internal_energy(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)
    V = volume(model, p, T, z; phase=phase, threaded=threaded)
    return VT_internal_energy(model,V,T,z)
end

function enthalpy(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)
    V = volume(model, p, T, z; phase=phase, threaded=threaded)
    return VT_enthalpy(model,V,T,z)
end

function gibbs_free_energy(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)
    V = volume(model, p, T, z; phase=phase, threaded=threaded)
    return VT_gibbs_free_energy(model,V,T,z)
end

function helmholtz_free_energy(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)
    V = volume(model, p, T, z; phase=phase, threaded=threaded)
    return VT_helmholtz_free_energy(model,V,T,z)
end

function isochoric_heat_capacity(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)
    V = volume(model, p, T, z; phase=phase, threaded=threaded)
    return VT_isochoric_heat_capacity(model,V,T,z)
end

function isobaric_heat_capacity(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)
    V = volume(model, p, T, z; phase=phase, threaded=threaded)
    return VT_isobaric_heat_capacity(model,V,T,z)
end

function isothermal_compressibility(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)
    V = volume(model, p, T, z; phase=phase, threaded=threaded)
    return VT_isothermal_compressibility(model,V,T,z)
end

function isentropic_compressibility(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)
    V = volume(model, p, T, z; phase=phase, threaded=threaded)
    return VT_isentropic_compressibility(model,V,T,z)
end

function speed_of_sound(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)
    V = volume(model, p, T, z; phase=phase, threaded=threaded)
    return VT_speed_of_sound(model,V,T,z)
end

function isobaric_expansivity(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)
    V = volume(model, p, T, z; phase=phase, threaded=threaded)
    return VT_isobaric_expansivity(model,V,T,z)
end

function joule_thomson_coefficient(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)
    V = volume(model, p, T, z; phase=phase, threaded=threaded)
    return VT_joule_thomson_coefficient(model,V,T,z)
end

function fugacity_coefficient(model::EoSModel,p,T,z=SA[1.]; phase = :unknown, threaded=true)
    μ_res  = chemical_potential_res(model,p,T,z;phase = phase, threaded=threaded)
    Z      = compressibility_factor(model,p,T,z;phase = phase, threaded=threaded)
    φ_i    = @. exp(μ_res/R̄/T)/Z
    return φ_i
end

function activity_coefficient(model::EoSModel,p,T,z=SA[1.]; phase = :unknown, threaded=true)
    pure   = split_model(model)
    μ_pure = chemical_potential.(pure,p,T;phase = phase, threaded=threaded)
    μ_mixt = chemical_potential(model,p,T,z;phase = phase, threaded=threaded)
    μ_pure = [μ_pure[i][1] for i in 1:length(z)]
    γ_i    = @. exp((μ_mixt-μ_pure)/R̄/T)/z
    return γ_i
end
"""
    compressibility_factor(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)

Calculates the compressibility factor `Z`, defined as:

```julia
Z = p*V(p)/R*T
```
the keywords `phase` and `threaded` are passed to the [volume solver](@ref Clapeyron.volume).
"""
function compressibility_factor(model::EoSModel, p, T, z=SA[1.]; phase = :unknown,threaded=true)
    V = volume(model, p, T, z; phase=phase, threaded=threaded)
    return VT_compressibility_factor(model, V, T, z)
end

function inversion_temperature(model::EoSModel, p, z=SA[1.0])
    T0 = 6.75*T_scale(model,z)
    μⱼₜ(T) = joule_thomson_coefficient(model,p,T,z)
    return Roots.find_zero(μⱼₜ,T0)
end

function molar_density(model::EoSModel,p,T,z=SA[1.0];phase = :unknown,threaded=true)
     V = volume(model,p,T,z;phase=phase,threaded=threaded)
     return sum(z)/V
end

function mass_density(model::EoSModel,p,T,z=SA[1.0];phase = :unknown,threaded=true)
    V = volume(model,p,T,z;phase=phase,threaded=threaded)
    molar_weight = molecular_weight(model,z)
    return molar_weight/V
end

export entropy, chemical_potential, internal_energy, enthalpy, gibbs_free_energy
export helmholtz_free_energy, isochoric_heat_capacity, isobaric_heat_capacity
export isothermal_compressibility, isentropic_compressibility, speed_of_sound
export isobaric_expansivity, joule_thomson_coefficient, compressibility_factor, inversion_temperature
export mass_density,molar_density, activity_coefficient, fugacity_coefficient


text:
reload("OBC.jl")
using Distributions

#srand(1)
blas_set_num_threads(1)

#if 'PARAM' in os.environ:
    #params = yaml.load(os.environ['PARAM'])
#else:
    #params = {}
params = Dict{ASCIIString,Any}()

function setv(p,s,d,conv=null)
    if !(s in keys(p))
    #if !contains(p,s)
        p[s] = d # Possible error TODO str(d)
        return d
    elseif conv != null
        return conv(p[s])
    else
        return p[s]
    end
end

int = x -> round(Int, x)
iters = setv(params, "iters", round(Int, 2e4), int)
num_feat = setv(params, "num_feat", 5, int)
rseed = setv(params, "rseed", rand(UInt), int)
#seed = setv(params, "seed", 1234, int)
seed = setv(params, "seed", rand(UInt), int)
Ntrn = setv(params, "Ntrn", 70, int)
Ntst = setv(params, "Ntst", 500, int)
f_glob = setv(params, "f_glob", 1, int)
subclasses = setv(params, "subclasses", 2, int)
f_het = setv(params, "f_het", 1, int)
f_rand = setv(params, "f_rand", 0, int)
rho = setv(params, "rho", 0.6, float)
f_tot = setv(params, "f_tot", f_glob+f_het*subclasses+f_rand, int)
blocksize = setv(params, "blocksize", 1, int)
mu0 = setv(params, "mu0", 0.2, float)
mu1 = setv(params, "mu1", 0.6, float)
sigma0 = setv(params, "sigma0", 0.1, float)
sigma1 = setv(params, "sigma1", 0.5, float)
lowd = setv(params, "lowd", 9.0, float)
highd = setv(params, "highd", 11.0, float)
numlam = setv(params, "numlam", 20, int)

function gen_data(mu, cov, n)
    lams = rand(MultivariateNormal(mu, cov), n)
    #println("lam stats:")
    #println("$(maximum(lams)) $(minimum(lams)) $(mean(lams)) $(std(lams))")
    ps = zeros(Int, size(lams))
    for i in 1:size(lams,1)
        for j in 1:size(lams,2)
            ps[i,j] = rand(Poisson(rand(Uniform(lowd,highd)) * exp(lams[i,j])))
        end
    end
    lams, ps
end

const kappa = 80
function gen_data_jason(mu, D)
    lmu = zeros(D) .+ mu

    cov = rand(InverseWishart(kappa,eye(D)*(kappa-D-1)))
    #cov = eye(D)*0.5
    #cov[1,2] = -0.1
    #cov[2,1] = -0.1
    lam1, trn_data = gen_data(lmu, cov, Ntrn)
    lam2, tst_data = gen_data(lmu, cov, Ntst)
    return lmu, cov, trn_data', tst_data', lam1, lam2
end

D = num_feat

######################################################################
# Class 0
######################################################################
trumu, trucov, dataa, tst_dataa,lam1,lam2 = gen_data_jason(0.0,D)
#tst_data = rand(510,D) .* 500 .+ 8
#tst_data[:,1] = 0
##tst_data[:,1] = rand(0:1, 510)
#dataa = tst_data[1:10,:]

######################################################################
# Class 1
######################################################################
trumub, trucovb, datab, tst_datab,_,_ = gen_data_jason(-0.0,D)
#tst_datab = rand(510,D) .* 500 .+ 8
#tst_datab[:,1] = 2
##tst_datab[:,1] = rand(0:1, 510)
#datab = tst_datab[1:10,:]

#d=100.0
#d1 = ones(Ntrn)*100.0 .+ randn(Ntrn)*10
#d2 = ones(Ntrn)*100.0 .+ randn(Ntrn)*10
#dataa = iround(dataa .* d1 ./ 100.0)
#datab = iround(datab .* d2 ./ 100.0)
#dmean1,dmean2 = mean(d1),mean(d2)
#dmean = mean([dmean1,dmean2])
d1 = d2 = dmean1 = dmean2 = 100.0
dataa = datab

cls = MPM.mpm_classifier(dataa, datab, burn=1000, thin=50, d1=d1, d2=d2, kappa=float(D+1), usepriors=true)
@time MPM.sample!(cls, 5000)#, verbose=true)

#@show bemc = MPM.bee_e_mc(cls, (dmean1,dmean2),numpts=1000)
#@show bee2s = MPM.bee_moments_2sample(cls, (dmean1,dmean2),numpts=50)
#@show beem1,beem2 = MPM.bee_moments(cls, (dmean1,dmean2),numpts=20)

#Profile.clear()
#@profile bemc = MPM.bee_e_mc(cls, (dmean1,dmean2),numpts=50)
#Profile.print()

#@show bemc = MPM.bee_e_mc(cls, (dmean1,dmean2),numpts=50)
#@show beem1,beem2 = MPM.bee_moments(cls, (dmean1,dmean2))
#@show beem1,beem2 = MPM.bee_moments(cls, (dmean1,dmean2))
#@show MPM.bee_moments_2sample(cls, (dmean1,dmean2),numpts=10)
#@show MPM.bee_moments_2sample(cls, (dmean1,dmean2),numpts=10)

#function project!(obj::MPM.MPMParams, sel)
    #obj.mu = obj.mu[sel]
    #obj.sigma = obj.sigma[sel,sel]
    #obj.sigpre = obj.sigpre[sel,sel]
    #obj.lam = obj.lam[sel,:]
    #return obj
#end

#function project(cls::OBC.BinaryClassifier,sel)
    #subcls = deepcopy(cls)
    #map!(x->project!(x,sel), subcls.mcmc1.db)
    #map!(x->project!(x,sel), subcls.mcmc2.db)
    #project!(subcls.cls1.curr,sel)
    #project!(subcls.cls1.old,sel)
    #project!(subcls.cls2.curr,sel)
    #project!(subcls.cls2.old,sel)
    #return subcls
#end

#for sel = combinations(1:3,2)
    #subcls = project(cls,sel)
    #@show sel
    #@show bemc = MPM.bee_e_mc(subcls, (dmean1,dmean2))
#end

#err = MPM.error_points(cls, [tst_dataa; tst_datab], [zeros(size(tst_dataa,1)), ones(size(tst_datab,1))], dmean=dmean)
#println("holdout error: $err")

#cls = MPM.mpm_classifier(dataa, datab, burn=1000, thin=50, d1=dmean, d2=dmean, usepriors=false)
#@time MPM.sample(cls, 10000)#, verbose=true)
#@show bemc = MPM.bee_e_mc(cls, dmean=dmean)
#@show beem1,beem2 = MPM.bee_moments(cls, dmean=dmean)

#for i=1:2
    #@time bemc = MPM.bee_e_mc(cls, dmean=dmean)
    #@time beem1,beem2 = MPM.bee_moments(cls, dmean=dmean)

    #@show bemc
    #@show beem1, beem2
    #@show beemse = beem2 - beem1^2
#end

#err = MPM.error_points(cls, [tst_dataa; tst_datab], [zeros(size(tst_dataa,1)), ones(size(tst_datab,1))], dmean=dmean)
#println("holdout error: $err")


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(Base.create_expr_cache),String,String,Array{Pair{Base.PkgId,UInt64},1},Base.UUID})
    Base.precompile(Tuple{typeof(Base.require),Module,Symbol})
end


text:
module RefractionCorrection

include("utils.jl")
include("refract.jl")
include("io.jl")

export cart2pol
export pol2cart
export correct!
export correct
export write!

end


text:
@testset "Prefix Parallel Function Tests" begin

  @test isdefined(SetAndForget, :prefix_parallel) == true

end