text:
###########Importing Libraries and Defining Structs############
using  LinearAlgebra, SparseArrays, StaticArrays, BenchmarkTools, Base, Plots

mutable struct Qubit #type for holding qubit information
    bitpos::Int
    vec::Vector{ComplexF64}
end

mutable struct Circuit
    qreg::Array{Qubit,1}
    gates::Vector{SparseMatrixCSC{ComplexF64,Int64}}
end

function Circuit(nqubits::Int)
    G = Vector{SparseMatrixCSC{ComplexF64,Int64}}()
    return Circuit(makeregister(nqubits),G)
end

function makeregister(nqubits::Int)
    #function for generating a quantum register of qubits
    zerovec = [1.0+0.0*im, 0.0+0.0*im]
    return [Qubit(i, copy(zerovec)) for i ∈ 1:nqubits]
end

function Base.max(set::AbstractRange)
    if set[1]<set[end]
        return set[end]
    elseif set[1]>set[end]
        return set[1]
    else
        return set[1]
    end
end

function Base.min(set::AbstractRange)
    if set[1]<set[end]
        return set[1]
    elseif set[1]>set[end]
        return set[end]
    else
        return set[1]
    end
end

######Single qubit gates ###########

function U!(C::Circuit, U::AbstractArray{ComplexF64,2}, bitpos::Int)
    #=Implementation that adds the new matrix into the gates field
    of the circuit=#
    a = length(C.qreg)
    if (bitpos>a) || (bitpos<=0)
        error("The bit specified must be within the register of the circuit!")
    end
    Uf = sparse(U)
    M = similar(Uf)
    I2 = sparse(I*(1.0+0.0*im),2,2)
    for  i ∈ 1:a
        if (i == 1) && (i !== bitpos)
            copy!(M,I2)
        elseif (i == 1) && (i == bitpos)
            copy!(M,Uf)
        else
            if (i == bitpos)
                M = kron(Uf,M)
            else
                M = kron(I2,M)
            end
        end
    end
   return push!(C.gates, M)
end

function U!(C::Circuit, U::AbstractSparseArray{ComplexF64,Int}, bitpos::Int)
    #=Implementation that adds the new matrix into the gates field
    of the circuit=#
    a = length(C.qreg)
    if (bitpos>a) || (bitpos<=0)
        error("The bit specified must be within the register of the circuit!")
    end
    M = similar(U)
    I2 = sparse(I*(1.0+0.0*im),2,2)
    for  i ∈ 1:a
        if (i == 1) && (i !== bitpos)
            copy!(M,I2)
        elseif (i == 1) && (i == bitpos)
            copy!(M,U)
        else
            if (i == bitpos)
                M = kron(U,M)
            else
                M = kron(I2,M)
            end
        end
    end
   return push!(C.gates, M)
end
##Specialized methods for doing multiple gates sequentially
function U!(C::Circuit, U::AbstractArray{ComplexF64,2}, set::AbstractVector{Int})
    #=Implementation that gives user ability to define
    a range of qubits to affect with the specified gate=#
    for i in set
        U!(C,U,i)
    end
end

function U!(C::Circuit, U::AbstractSparseArray{ComplexF64,Int}, set::AbstractVector{Int})
    #=Implementation that gives user ability to define
    a range of qubits to affect with the specified gate=#
    for i in set
        U!(C,U,i)
    end
end

function H!(C::Circuit, bitpos::Int)
    #=Implementation that adds the new matrix into the gates field
    of the circuit=#
    H = sparse([1.0/sqrt(2)+0.0*im  1.0/sqrt(2)+0.0*im
                1.0/sqrt(2)+0.0*im -1.0/sqrt(2)+0.0*im])
    U!(C,H,bitpos)
end

function H!(C::Circuit, set::AbstractVector{Int})
    #=Implementation that gives user ability to define
    a range of qubits to affect with the specified gate=#
    H = sparse([1.0/sqrt(2)+0.0*im  1.0/sqrt(2)+0.0*im
                1.0/sqrt(2)+0.0*im -1.0/sqrt(2)+0.0*im])
    U!(C,H,set)
end

function T!(C::Circuit, bitpos::Int)
    #=Implementation that adds the new matrix into the gates field
    of the circuit=#
    T = sparse([1.0+0.0*im  0.0+0.0*im
                0.0+0.0*im  exp(im*π/4)])
    U!(C,T,bitpos)
end

function T!(C::Circuit, set::AbstractVector{Int})
    #=Implementation that gives user ability to define
    a range of qubits to affect with the specified gate=#
    T = sparse([1.0+0.0*im  0.0+0.0*im
                0.0+0.0*im  exp(im*π/4)])
    U!(C,T,set)
end

######## Pauli Gates ##########
function X!(C::Circuit, bitpos::Int)
    #=Implementation that adds the new matrix into the gates field
    of the circuit=#
    X = sparse([0.0+0.0*im  1.0+0.0*im
                1.0+0.0*im  0.0+0.0*im])
    U!(C,X,bitpos)
end

function X!(C::Circuit, set::AbstractVector{Int})
    #=Implementation that gives user ability to define
    a range of qubits to affect with the specified gate=#
    X = sparse([0.0+0.0*im  1.0+0.0*im
                1.0+0.0*im  0.0+0.0*im])
    U!(C,X,set)
end

function Y!(C::Circuit, bitpos::Int)
    #=Implementation that adds the new matrix into the gates field
    of the circuit=#
    Y = sparse([0.0+0.0*im  0.0-1.0*im
                0.0+1.0*im  0.0+0.0*im])
    U!(C,Y,bitpos)
end

function Y!(C::Circuit, set::AbstractVector{Int})
    #=Implementation that gives user ability to define
    a range of qubits to affect with the specified gate=#
    Y = sparse([0.0+0.0*im  0.0-1.0*im
                0.0+1.0*im  0.0+0.0*im])
    U!(C,Y,set)
end

function Z!(C::Circuit, bitpos::Int)
    #=Implementation that adds the new matrix into the gates field
    of the circuit=#
    Z = sparse([1.0+0.0*im  0.0+0.0*im
                0.0+0.0*im -1.0+0.0*im])
    U!(C,Z,bitpos)
end

function Z!(C::Circuit, set::AbstractVector{Int})
    #=Implementation that gives user ability to define
    a range of qubits to affect with the specified gate=#
    Z = sparse([1.0+0.0*im  0.0+0.0*im
                0.0+0.0*im -1.0+0.0*im])
    U!(C,Z,set)
end

############# Controlled Gates and 2 qubit Gates#############
function CU!(C::Circuit, U::AbstractArray{ComplexF64,2},control::Int,target::Int)
    #Implentation of CU gate using Kronnecker products.
    Uf = sparse(U)

    a = length(C.qreg)

    if (target==control)
        error("Target qubit and control qubit must be Different!")
    elseif ((target>a) || (control>a)) || ((target<=0) || (control<=0))
        error("The target and control qubit must be within the size of the register!")
    end

    M1 = spzeros(ComplexF64,2,2)
    M2 = copy(M1)
    I2 = sparse(I*(1.0+0.0*im),2,2)

    conket1 = sparse([0.0+0.0*im 0.0+0.0*im;
                      0.0+0.0*im 1.0+0.0*im])
    conket0 = sparse([1.0+0.0*im 0.0+0.0*im
                      0.0+0.0*im 0.0+0.0*im])

    for i ∈ 1:a

        if (i==1)
            if (i==control)
                M1 = conket1
            elseif (i==target)
                M1 = Uf
            else
                M1 = I2
            end
        else
            if (i==control)
                M1 = kron(conket1,M1)
            elseif (i==target)
                M1 = kron(Uf,M1)
            else
                M1 = kron(I2,M1)
            end
        end
    end

    for i ∈ 1:a

        if (i==1==control)
            M2 = conket0
        elseif (i==1!==control)
            M2 = I2
        else
            if (i==control)
                M2 = kron(conket0,M2)
            else
                M2 = kron(I2,M2)
            end
        end
    end
    return push!(C.gates,M1+M2)
end

function CU!(C::Circuit, U::AbstractSparseArray{ComplexF64,Int},control::Int,target::Int)
    #Implentation of CU gate using Kronnecker products.
    a = length(C.qreg)

    if (target==control)
        error("Target qubit and control qubit must be Different!")
    elseif ((target>a) || (control>a)) || ((target<=0) || (control<=0))
        error("The target and control qubit must be within the size of the register!")
    end

    M1 = spzeros(ComplexF64,2,2)
    M2 = copy(M1)
    I2 = sparse(I*(1.0+0.0*im),2,2)

    conket1 = sparse([0.0+0.0*im 0.0+0.0*im;
                      0.0+0.0*im 1.0+0.0*im])
    conket0 = sparse([1.0+0.0*im 0.0+0.0*im
                      0.0+0.0*im 0.0+0.0*im])

    for i ∈ 1:a

        if (i==1)
            if (i==control)
                M1 = conket1
            elseif (i==target)
                M1 = U
            else
                M1 = I2
            end
        else
            if (i==control)
                M1 = kron(conket1,M1)
            elseif (i==target)
                M1 = kron(U,M1)
            else
                M1 = kron(I2,M1)
            end
        end
    end

    for i ∈ 1:a

        if (i==1==control)
            M2 = conket0
        elseif (i==1!==control)
            M2 = I2
        else
            if (i==control)
                M2 = kron(conket0,M2)
            else
                M2 = kron(I2,M2)
            end
        end
    end
    return push!(C.gates,M1+M2)
end
###Methods for doing a number of gates in a row
function CU!(C::Circuit,U::AbstractArray{ComplexF64,2},controls::AbstractVector{Int},targets::AbstractVector{Int})
    #convenience method for doing mutliple CU operations at a time.
    a = length(C.qreg)
    b = length(controls)
    if length(targets)!==b
        error("The set of controls and targets must be the same length or the operation is ambiguous.")
    elseif ((max(controls)>a) || (max(targets)>a)) || ((min(controls)<=0) || (min(targets)<=0))
        error("The target and control qubits must all lie within the register!")
    end

    Uf = sparse(U)
    for i in 1:b
        CU!(C,Uf,controls[i],targets[i])
    end
end

function CU!(C::Circuit,U::AbstractArray{ComplexF64,2},control::Int,targets::AbstractVector{Int})
    #convenience method for doing mutliple CU operations at a time.
    a = length(C.qreg)
    if ((control>a) || (max(targets)>a)) || ((control<=0) || (min(targets)<=0))
        error("The target and control qubits must all lie within the register!")
    end

    Uf = sparse(U)
    for i in targets
        CU!(C,Uf,control,i)
    end
end

function CU!(C::Circuit,U::AbstractArray{ComplexF64,2},controls::AbstractVector{Int},target::Int)
    #convenience method for doing mutliple CU operations at a time.
    a = length(C.qreg)
    if ((max(controls)>a) || (target>a)) || ((min(controls)<=0) || (target<=0))
        error("The target and control qubits must all lie within the register!")
    end

    Uf = sparse(U)
    for i in controls
        CU!(C,Uf,i,target)
    end
end

function CU!(C::Circuit,U::AbstractSparseArray{ComplexF64,Int},controls::AbstractVector{Int},targets::AbstractVector{Int})
    #convenience method for doing mutliple CU operations at a time.
    a = length(C.qreg)
    b = length(controls)
    if length(targets)!==b
        error("The set of controls and targets must be the same length or the operation is ambiguous.")
    elseif ((max(controls)>a) || (max(targets)>a)) || ((min(controls)<=0) || (min(targets)<=0))
        error("The target and control qubits must all lie within the register!")
    end

    for i in 1:b
        CU!(C,U,controls[i],targets[i])
    end
end

function CU!(C::Circuit,U::AbstractSparseArray{ComplexF64,Int},control::Int,targets::AbstractVector{Int})
    #convenience method for doing mutliple CU operations at a time.
    a = length(C.qreg)
    if ((control>a) || (max(targets)>a)) || ((control<=0) || (min(targets)<=0))
        error("The target and control qubits must all lie within the register!")
    end

    for i in targets
        CU!(C,U,control,i)
    end
end

function CU!(C::Circuit,U::AbstractSparseArray{ComplexF64,Int},controls::AbstractVector{Int},target::Int)
    #convenience method for doing mutliple CU operations at a time.
    a = length(C.qreg)
    if ((max(controls)>a) || (target>a)) || ((min(controls)<=0) || (target<=0))
        error("The target and control qubits must all lie within the register!")
    end

    for i in controls
        CU!(C,U,i,target)
    end
end

function CX!(C::Circuit,control::Int,target::Int)
    #Implentation of CNOT gate using Kronnecker products.
    X = sparse([0.0+0.0*im 1.0+0.0*im;
                1.0+0.0*im 0.0+0.0*im])

    CU!(C,X,control,target)
end

function CX!(C::Circuit,controls::AbstractVector{Int},targets::AbstractVector{Int})
    #convenience method for doing mutliple CX operations at a time.
    X = sparse([0.0+0.0*im 1.0+0.0*im;
                1.0+0.0*im 0.0+0.0*im])

    CU!(C,X,controls,targets)
end

function CX!(C::Circuit,control::Int,targets::AbstractVector{Int})
    #convenience method for doing mutliple CX operations at a time.
    X = sparse([0.0+0.0*im 1.0+0.0*im;
                1.0+0.0*im 0.0+0.0*im])

    CU!(C,X,control,targets)
end

function CX!(C::Circuit,controls::AbstractVector{Int},target::Int)
    #convenience method for doing mutliple CX operations at a time.
    X = sparse([0.0+0.0*im 1.0+0.0*im;
                1.0+0.0*im 0.0+0.0*im])

    CU!(C,X,controls,target)
end

function CZ!(C::Circuit,control::Int,target::Int)
    #Implentation of CZ gate using Kronnecker products.
    Z = sparse([1.0+0.0*im  0.0+0.0*im;
                0.0+0.0*im -1.0+0.0*im])

    CU!(C,Z,control,target)
end

function CZ!(C::Circuit,controls::AbstractVector{Int},targets::AbstractVector{Int})
    #convenience method for doing mutliple CZ operations at a time.
    Z = sparse([1.0+0.0*im  0.0+0.0*im;
                0.0+0.0*im -1.0+0.0*im])

    CU!(C,Z,controls,targets)
end

function CZ!(C::Circuit,control::Int,targets::AbstractVector{Int})
    #convenience method for doing mutliple CZ operations at a time.
    Z = sparse([1.0+0.0*im  0.0+0.0*im;
                0.0+0.0*im -1.0+0.0*im])

    CU!(C,Z,control,targets)
end

function CZ!(C::Circuit,controls::AbstractVector{Int},target::Int)
    #convenience method for doing mutliple CZ operations at a time.
    Z = sparse([1.0+0.0*im  0.0+0.0*im;
                0.0+0.0*im -1.0+0.0*im])

    CU!(C,Z,controls,target)
end

function Swap!(C::Circuit,b1::Int,b2::Int)
    #quick and dirty implmentation of Swap Gate
    a = length(C.qreg)
    if b1==b2
        error("This ins't really a swap, you just tried to swap a qubit with itself!")
    elseif ((b1>a) || (b2>a)) || ((b1<=0) || (b2<=0))
        error("Both bits must be within the size of the register!")
    end
    CX!(C,b2,b1)
    CX!(C,b1,b2)
    CX!(C,b2,b1)
end

function CCX!(C::Circuit,c1::Int,c2::Int,target::Int)
    #Implementation of the Toffoli gate using Kronnecker products.
    X = sparse([0.0+0.0*im 1.0+0.0*im;
                1.0+0.0*im 0.0+0.0*im])

    a = length(C.qreg)

    if (target==c1) || (target==c2)
        error("Target qubit and control qubits must be Different!")
    elseif (c1==c2)
        error("The control qubits must be different from one another!")
    elseif ((target>a) || (c1>a) || (c2>a)) || ((target<=0) || (c1<=0) ||(c2<=0) )
        error("The target and control qubits must be within the size of the register!")
    end

    M1 = spzeros(ComplexF64,2,2)
    M2 = copy(M1)
    M3 = copy(M1)
    M4 = copy(M1)
    I2 = sparse(I*(1.0+0.0*im),2,2)

    conket1 = sparse([0.0+0.0*im 0.0+0.0*im;
                      0.0+0.0*im 1.0+0.0*im])
    conket0 = sparse([1.0+0.0*im 0.0+0.0*im
                      0.0+0.0*im 0.0+0.0*im])

    for i ∈ 1:a

        if (i==1)
            if (i==c1)||(i==c2)
                M1 = conket1
            elseif (i==target)
                M1 = X
            else
                M1 = I2
            end
        else
            if (i==c1)||(i==c2)
                M1 = kron(conket1,M1)
            elseif (i==target)
                M1 = kron(X,M1)
            else
                M1 = kron(I2,M1)
            end
        end
    end

    for i ∈ 1:a

        if (i==1)
            if (i==c1)||(i==c2)
                M2 = conket0
            else
                M2 = I2
            end
        else
            if (i==c1)||(i==c2)
                M2 = kron(conket0,M2)
            else
                M2 = kron(I2,M2)
            end
        end
    end

    for i ∈ 1:a

        if (i==1)
            if (i==c1)
                M3 = conket0
            elseif (i==c2)
                M3 = conket1
            else
                M3 = I2
            end
        else
            if (i==c1)
                M3 = kron(conket0,M3)
            elseif (i==c2)
                M3 = kron(conket1,M3)
            else
                M3 = kron(I2,M3)
            end
        end
    end

    for i ∈ 1:a

        if (i==1)
            if (i==c1)
                M4 = conket1
            elseif (i==c2)
                M4 = conket0
            else
                M4 = I2
            end
        else
            if (i==c1)
                M4 = kron(conket1,M4)
            elseif (i==c2)
                M4 = kron(conket0,M4)
            else
                M4 = kron(I2,M4)
            end
        end
    end
    return push!(C.gates,M1+M2+M3+M4)

end

########## Running the Circuit ###########
function Run_circ(C::Circuit)

    maxiter = length(C.gates)
    maxbit = length(C.qreg)

    M = copy(C.gates[1])
    ket = copy(C.qreg[1].vec)
    for i ∈ 2:maxbit
    ket = kron(C.qreg[i].vec,ket)
    end

    for i ∈ 2:maxiter
        M .= C.gates[i]*M
    end

    return M*ket
end

function shots(C::Circuit, nshots::Int)

    finalket = Run_circ(C)
    ntrials = rand(nshots)
    ensemble = zeros(Int,length(finalket))
    trials!(ensemble, finalket, ntrials)
end

function trials!(ensemble::AbstractVector{Int}, ket::AbstractVector{ComplexF64}, probs::AbstractVector{Float64})
    nshots = length(probs)

    @inbounds for j ∈ 1:nshots
        local p = probs[j]
        local sumold = 0.0
        local sumnew = 0.0

        @inbounds for i ∈ 1:length(ket)
            sumnew = sumold + real(conj(ket[i])*ket[i])
            if (p<=sumnew)
                ensemble[i] += 1
                break
            end
            sumold = sumnew
        end
    end
    return ensemble
end

function graph(data::AbstractVector)
    pyplot()
    max = length(data)-1
    nqubits = convert(Int,log(2,max+1))
    xs = [string(i,base=2,pad=nqubits) for i ∈ 0:max]

    bar(xs,data, leg=false, title="Simulation of Circuit", xlabel="qubit readout", ylabel="Counts", dpi = 300)
end


text:
@testset "Problems" begin
    @testset "Can create fixed dimensional problem with search space as array" begin
        p = BlackBoxOptim.fixeddim_problem((x)->x; search_space = [(0.0, 1.0)])
        @test typeof(p) === BlackBoxOptim.FixedDimProblem
    end

    @testset "Exception if no valid search space supplied when creating fixed dimensional problem" begin
        @test_throws BlackBoxOptim.fixeddim_problem((x)->x; search_space = (0.0, 1.0))
    end

end

ext: jl


text:
using GradChecker
using Test

@testset "GradChecker.jl" begin
    # Write your tests here.
end


text:
using PDMats  # Positive Definite Matrices
using QuadGK  # 1-d integration
using Cuba    # multi-d integration
using Optim   # Optimizing over non-linear parameters

# Read dataset
data = readdlm("../data/rvs_0001.txt");
times = data[:,1]
vels = data[:,2]
sigma_obs = data[:,3]
nobs = length(vels)

# Set priors
const Cmax = 1000.0
function logprior_offset(param)
  -log(2*Cmax)
end

const sigma_j_min = 0.0
const sigma_j_max = 99.0
const sigma_j_0 = 1.0
function logprior_jitter(sigma_j::Real)
  -log1p(sigma_j/sigma_j_0)-log(log1p(sigma_j_max/sigma_j_0)) # -log(sigma_j_0)
end

const Pmin = 1.0
const Pmax = 1.0e4
function logprior_period(P::Real)
   -log(P)-log(log(Pmax/Pmin))
end
function logprior_periods(P::Vector{T}) where T<:Real
   if length(P) == 0  return 0 end
   -sum(log.(P))-length(P)*log(log(Pmax/Pmin))
end

const Kmin = 0.0
const Kmax = 999.0
const K0 = 1.0
function logprior_amplitude(K::Real)
  logp = -log1p(K/K0)-log(log1p(Kmax/K0)) # -log(K0)
end

function logprior_sinusoids(param)
  const nparam = length(param)
  const namp_per_pl = 2
  npl = convert(Int64,floor((nparam-1)//namp_per_pl))
  logp = 0
  logj = 0
  for i in 0:(npl-1)
    K = sqrt(param[1+namp_per_pl*i]^2+param[2+namp_per_pl*i]^2)
    logp += logprior_amplitude(K)
	logj -= log(K)
  end
  logp -= npl*log(2pi)
  C = param[nparam]
  if !(-Cmax<=C<=Cmax) logp -=Inf
  else  logp -= log(2*Cmax) end
  return (logp,logj)
end

function logprior_epicycles(param)
  nparam = length(param)
  namp_per_pl = 4
  npl = convert(Int64,floor((nparam-1)//namp_per_pl))
  logp = 0
  logj = 0
  for i in 0:(npl-1)
    K  = sqrt(param[1+namp_per_pl*i]^2+param[2+namp_per_pl*i]^2)
	eK = sqrt(param[3+namp_per_pl*i]^2+param[4+namp_per_pl*i]^2)
	if !(0.0<=eK<K) logp -= Inf end
    logp += logprior_amplitude(K)
	logj -= 2*log(K)+log(eK)
  end
  logp -= 2*npl*log(2pi)
  C = param[nparam]
  if !(-Cmax<=C<=Cmax) logp -=Inf
  else  logp -= log(2*Cmax) end
  return (logp,logj)
end

#sincos from LombScargle.jl
# `sincos` is not yet in Julia, but it's available through the math library.  The following
# definitions are from Yichao Yu's code at https://github.com/nacs-lab/yyc-data and are
# likely to be later added to Julia (so they can be removed here).  See also
# https://discourse.julialang.org/t/poor-performance-of-the-cis-function/3402.
if !isdefined(Base, :sincos)
    @inline function sincos(v::Float64)
        Base.llvmcall("""
        %f = bitcast i8 *%1 to void (double, double *, double *)*
        %pres = alloca [2 x double]
        %p1 = getelementptr inbounds [2 x double], [2 x double]* %pres, i64 0, i64 0
        %p2 = getelementptr inbounds [2 x double], [2 x double]* %pres, i64 0, i64 1
        call void %f(double %0, double *nocapture noalias %p1, double *nocapture noalias %p2)
        %res = load [2 x double], [2 x double]* %pres
        ret [2 x double] %res
        """, Tuple{Float64,Float64}, Tuple{Float64,Ptr{Void}}, v,
                      cglobal((:sincos, Base.libm_name)))
    end
    @inline function sincos(v::Float32)
        Base.llvmcall("""
        %f = bitcast i8 *%1 to void (float, float *, float *)*
        %pres = alloca [2 x float]
        %p1 = getelementptr inbounds [2 x float], [2 x float]* %pres, i64 0, i64 0
        %p2 = getelementptr inbounds [2 x float], [2 x float]* %pres, i64 0, i64 1
        call void %f(float %0, float *nocapture noalias %p1, float *nocapture noalias %p2)
        %res = load [2 x float], [2 x float]* %pres
        ret [2 x float] %res
            """, Tuple{Float32,Float32}, Tuple{Float32,Ptr{Void}}, v,
                      cglobal((:sincosf, Base.libm_name)))
    end
    @inline function sincos(v)
        @fastmath (sin(v), cos(v))
    end
end

# Fit planet models
function make_design_matrix_circ(times::Vector{T},periods::Vector{T}) where T<:Real
  const namp_per_pl = 2
  nparam = 1+namp_per_pl*length(periods)
  freqs = 2pi./periods
  X = Array{T}(length(times),nparam)
  for i in 1:length(periods)
     for j in 1:length(times)
       (X[j,1+(i-1)*namp_per_pl],X[j,2+(i-1)*namp_per_pl]) = sincos(freqs[i]*times[j])
	 end
  end
  X[:,nparam] = one(T)
  return X
end

function make_design_matrix_circ_old(times::Vector{T},periods::Vector{T}) where T<:Real
  const namp_per_pl = 2
  nparam = 1+namp_per_pl*length(periods)
  freqs = 2pi./periods
  X = Array{T}(length(times),nparam)
  for i in 1:length(periods)
     X[:,1+(i-1)*namp_per_pl] = sin.(freqs[i]*times)
	 X[:,2+(i-1)*namp_per_pl] = cos.(freqs[i]*times)
  end
  X[:,nparam] = one(T)
  return X
end

function update_design_matrix_circ(X::Array{T,2}, times::Vector{T},periods::Vector{T}, pl_to_update::Integer) where T<:Real
  const namp_per_pl = 2
  nparam = 1+namp_per_pl*length(periods)
  freq = 2pi./periods[pl_to_update]
  @assert size(X) == (length(times),nparam)
  i = pl_to_update
  idx_s = 1+(i-1)*namp_per_pl
  idx_c = 2+(i-1)*namp_per_pl
  for j in 1:length(times)
     (X[j,idx_s],X[j,idx_c]) = sincos(freq*times[j])
  end
  return X
end

function make_design_matrix_epicycle(times::Vector{T},periods::Vector{T}) where T<:Real
  const namp_per_pl = 4
  nparam = 1+namp_per_pl*length(periods)
  freqs = 2pi./periods
  X = Array{T}(length(times),nparam)
  for i in 1:length(periods)
     s = sin.(freqs[i]*times)
	 c = cos.(freqs[i]*times)
	 X[:,1+(i-1)*namp_per_pl] = s
	 X[:,2+(i-1)*namp_per_pl] = c
	 X[:,3+(i-1)*namp_per_pl] = 2*s.*c
	 X[:,4+(i-1)*namp_per_pl] = c.^2.-s.^2
  end
  X[:,nparam] = one(T)
  return X
end

function update_design_matrix_epicycle(X::Array{T,2}, times::Vector{T},periods::Vector{T}, pl_to_update::Integer) where T<:Real
  const namp_per_pl = 4
  nparam = 1+namp_per_pl*length(periods)
  freq = 2pi./periods[pl_to_update]
  @assert size(X) == (length(times),nparam)
  i = pl_to_update
  idx_s = 1+(i-1)*namp_per_pl
  for j in 1:length(times)
	 (s,c) = sincos(freq*times[j])
     X[j,idx_s  ] = s
	 X[j,idx_s+1] = c
	 X[j,idx_s+2] = 2*s.*c
	 X[j,idx_s+3] = c.^2.-s.^2
  end
  return X
end

function calc_best_fit(data::Vector{T}, model::Array{T,2}, Sigma::PDMat{T}) where T<:Real
  inv_chol_Sigma = inv(Sigma.chol[:L])
  param_bestfit = (inv_chol_Sigma*model) \ (inv_chol_Sigma*data)
end

function calc_best_fit(data::Vector{T}, model::Array{T,2}, Sigma::PDiagMat{T}) where T<:Real
  inv_chol_Sigma = sqrt.(Sigma.inv_diag)
  param_bestfit = (inv_chol_Sigma.*model) \ (inv_chol_Sigma.*data)
end

function calc_chisq(data::Vector{T}, model::Array{T,2}, Sigma::SigmaT; FIM::PDMat{T} = PDMat(Xt_invA_X(Sigma,model)), calclogprior::Function = x->(0.0,0.0) ) where {T<:Real, SigmaT<:AbstractPDMat{T} }
  @assert length(data) == size(model,1) == size(Sigma,1)
  nobs = length(data)
  nparam = size(model,2)

  param_bf_linalg = calc_best_fit(data,model,Sigma)
  predict = model*param_bf_linalg
  delta = data.-predict
  chisq = invquad(Sigma,delta) # sum(delta'*Sigma^-1*delta)
end

function compute_Laplace_Approx_param_and_integral(data::Vector{T}, model::Array{T,2}, Sigma::SigmaT; FIM::PDMat{T} = PDMat(Xt_invA_X(Sigma,model)), calclogprior::Function = x->(0.0,0.0) ) where {T<:Real, SigmaT<:AbstractPDMat{T} }
  @assert length(data) == size(model,1) == size(Sigma,1)
  nobs = length(data)
  nparam = size(model,2)

  param_bf_linalg = calc_best_fit(data,model,Sigma)
  predict = model*param_bf_linalg
  delta = data.-predict
  chisq = invquad(Sigma,delta) # sum(delta'*Sigma^-1*delta)

  loglikelihood_mode = 0.5*(-chisq-logdet(Sigma)-nobs*log(2pi))
  (logprior_mode, log_jacobian) = calclogprior(param_bf_linalg)
  log_mode = logprior_mode + loglikelihood_mode

  sigma_param = inv(FIM).chol[:U]
  LaplaceApprox = 0.5*nparam*log(2pi)-0.5*logdet(FIM)+log_mode + log_jacobian
  return (param_bf_linalg,sigma_param,LaplaceApprox)
end

function compute_Laplace_Approx(data::Vector{T}, model::Array{T,2}, Sigma::SigmaT; FIM::PDMat{T} = PDMat(Xt_invA_X(Sigma,model)), calclogprior::Function = x->(0.0,0.0) ) where {T<:Real, SigmaT<:AbstractPDMat{T} }
  @assert length(data) == size(model,1) == size(Sigma,1)
  nobs = length(data)
  nparam = size(model,2)

  param_bf_linalg = calc_best_fit(data,model,Sigma)
  predict = model*param_bf_linalg
  delta = data.-predict
  chisq = invquad(Sigma,delta) # sum(delta'*Sigma^-1*delta)

  loglikelihood_mode = 0.5*(-chisq-logdet(Sigma)-nobs*log(2pi))
  (logprior_mode, log_jacobian) = calclogprior(param_bf_linalg)
  log_mode = logprior_mode + loglikelihood_mode

  LaplaceApprox = 0.5*nparam*log(2pi)-0.5*logdet(FIM)+log_mode + log_jacobian
  return LaplaceApprox
end

function generate_K_samples(bf::Vector, sigma, nsamples::Integer = 1)
  nparam = length(bf)
  const namp_per_pl = 2
  theta = bf .+ sigma*randn(nparam,nsamples)
  #theta = param_bf'.+randn(nsamples,nparam)*sigma_param
  npl = convert(Int64,floor((nparam-1)//namp_per_pl))

  Ks = zeros(nsamples,npl)
  for i in 1:npl
    Ks[:,i] = sqrt.(theta[1+namp_per_pl*(i-1),:].^2+theta[2+namp_per_pl*(i-1),:].^2)
  end
  return Ks
end

function generate_Ke_samples(bf::Vector, sigma, nsamples::Integer = 1)
  nparam = length(bf)
  const namp_per_pl = 4
  theta = bf .+ sigma*randn(nparam,nsamples)
  npl = convert(Int64,floor((nparam-1)//namp_per_pl))

  Ks = zeros(nsamples,npl)
  es = zeros(nsamples,npl)
  hs = zeros(nsamples,npl)
  ks = zeros(nsamples,npl)
  for i in 1:npl
    Ks[:,i] = sqrt.(theta[1+namp_per_pl*(i-1),:].^2+theta[2+namp_per_pl*(i-1),:].^2)
	hs[:,i] = theta[3+namp_per_pl*(i-1),:] ./ Ks[:,i]
	ks[:,i] = theta[4+namp_per_pl*(i-1),:] ./ Ks[:,i]
	es[:,i] = sqrt.(hs[:,i].^2.+ks[:,i].^2)
  end
  return (Ks,es,hs,ks)
end

# Now with covariances
function kernel_quaesiperiodic(dt::Float64)
  alpha_sq     = 3.0
  tau_evolve   = 50.0
  lambda_p     = 0.5
  tau_rotate   = 20.0
  alpha_sq*exp(-0.5*((dt/tau_evolve)^2 + (sin(pi*dt/tau_rotate)/lambda_p)^2) )
end

function make_Sigma(t::Vector{Float64}, sigma_obs::Vector{Float64}, sigma_j::Float64)
  PDMat( [ kernel_quaesiperiodic(t[i]-t[j]) + (i==j ? sigma_obs[i]^2 + sigma_j^2 : 0.0) for i in 1:length(t), j in 1:length(t)] )
end

function make_dSigmadsigmaj(t::Vector{Float64}, sigma_obs::Vector{Float64}, sigma_j::Float64)
  #PDMat( [ 0.0 + (i==j ? 2*sigma_j : 0.0) for i in 1:length(t), j in 1:length(t)] )
  PDiagMat( 2*sigma_j*ones(length(t)) )
  ScalMat(length(t),2*sigma_j)
end

#=
function compute_log_evidence(sigma_j::Real)
  Sigma = make_Sigma(times,sigma_obs,sigma_j)
  log_prior_outside = logprior_jitter(sigma_j) + logprior_periods(periods)
  log_likelihood = compute_Laplace_Approx(vels,X,Sigma,calclogprior=logprior_epicycles)
  log_prior_outside+log_likelihood
end

function compute_evidence(sigma_j::Real; log_normalization::Real=0.0)
  exp(compute_log_evidence(sigma_j)-log_normalization)
end
=#

#=
periods = [42.0]
X = make_design_matrix_epicycle(times,periods);
Sigma = PDMat(make_Sigma(times,sigma_obs,sigma_j_0))
FIM = PDMat(Xt_invA_X(Sigma,X))
(param_bf,sigma_param,evid) = compute_Laplace_Approx_param_and_integral(vels,X,Sigma,FIM=FIM,calclogprior=logprior_epicycles)

(Ks,es) = generate_Ke_samples(param_bf,sigma_param,1000)
=#

#=
lognorm = compute_log_evidence(sigma_j_0)  # can help to prevent underflow
integrand(x) = compute_evidence(x,log_normalization=lognorm)
result = QuadGK.quadgk(integrand,sigma_j_min,sigma_j_max).*(exp(lognorm)/log(10))
=#

function calc_opt_sigmaj(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}, model::Array{T,2}; calclogprior_linear::Function = logprior_sinusoids) where T<:Real
  function calc_laplace_approx_arg_jitter(sigmaj::Real)
    Sigma = Sigma0 + ScalMat(length(times),sigmaj) # PDMat(make_Sigma(times,sigma_obs,sigmaj))
    FIM = PDMat(Xt_invA_X(Sigma,model))
    evid = compute_Laplace_Approx(vels,model,Sigma,FIM=FIM,calclogprior=calclogprior_linear)
	evid += logprior_jitter(sigmaj)
    -evid
  end
  Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.))
  result = optimize(calc_laplace_approx_arg_jitter,sigma_j_min,sigma_j_max,rel_tol=1e-6)
end

#=
periods = [42.0]
X = make_design_matrix_sinusoids(times,periods);
res = calc_opt_sigmaj(times,vels,sigma_obs,X)
dSdj = make_dSigmadsigmaj(times,sigma_obs,res.minimizer)
(compute_log_evidence(res.minimizer)+0.5*(log(2pi)-log(0.5*trace((Sigma \ full(dSdj))^2))))/log(10)
=#

#=
function make_matrix_pd(A::Array{Float64,2}; epsabs::Float64 = 0.0, epsfac::Float64 = 1.0e-6)
  @assert(size(A,1)==size(A,2))
  B = copy(A)
  itt = 1
  while !isposdef(B)
	eigvalB,eigvecB = eig(B)
	neweigval = (epsabs == 0.0) ? epsfac*minimum(eigvalB[eigvalB.>0]) : epsabs
	eigvalB[eigvalB.<0] = neweigval
	B = eigvecB *diagm(eigvalB)*eigvecB'
#	println(itt,": ",B)
	itt +=1
	if itt>size(A,1)
	  error("There's a problem in make_matrix_pd.\n")
  	  break
	end
  end
  return B
end
=#

# Integrate over period
function compute_marginalized_likelihood_fix_periods(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}, periods::Vector{T}; calclogprior_linear::Function = logprior_sinusoids, make_design_matrix::Function = make_design_matrix_circ, X = make_design_matrix(times,periods), Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.0)) ) where T<:Real
   const nparam = size(X,2)

   function optimize_jitter_helper(sigmaj::Real)
     Sigma = Sigma0 + ScalMat(length(times),sigmaj) # PDMat(make_Sigma(times,sigma_obs,sigmaj))
	 FIM = Xt_invA_X(Sigma,X)
	 if !isposdef(FIM)
	   FIM += diagm(1e-3*ones(nparam))
	 end
	 FIM = PDMat(FIM)
     evid = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=calclogprior_linear)
 	 evid += logprior_jitter(sigmaj)
	 -evid
   end
   result = optimize(optimize_jitter_helper,sigma_j_min,sigma_j_max,rel_tol=1e-6)
   #println(result)
   sigmaj = result.minimizer
   #Sigma = PDMat(make_Sigma(times,sigma_obs,sigmaj))
   Sigma = Sigma0 + ScalMat(length(times),sigmaj) # PDMat(make_Sigma(times,sigma_obs,sigmaj))
   dSdj = make_dSigmadsigmaj(times,sigma_obs,sigmaj)
   FIM = Xt_invA_X(Sigma,X)
	 if !isposdef(FIM)
	   FIM += diagm(1e-3*ones(nparam))
	 end
   FIM = PDMat(FIM)

   logp = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=calclogprior_linear)
   logp += logprior_jitter(sigmaj)+ logprior_periods(periods)
   logp += 0.5*(log(2pi)-log(0.5*trace((Sigma \ full(dSdj))^2)))
   #println("# logp = ",logp)
   return logp
end

function compute_marginalized_likelihood_fix_periods_sigmaj(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}, periods::Vector{T}, sigmaj::T; calclogprior_linear::Function = logprior_sinusoids, make_design_matrix::Function = make_design_matrix_circ, Sigma=PDMat(make_Sigma(times,sigma_obs,sigmaj)), X= make_design_matrix(times,periods) ) where T<:Real
   #X = make_design_matrix(times,periods)
   const nparam = size(X,2)
   #Sigma = PDMat(make_Sigma(times,sigma_obs,sigmaj))
   #dSdj = make_dSigmadsigmaj(times,sigma_obs,sigmaj)
   FIM = Xt_invA_X(Sigma,X)
	 if !isposdef(FIM)
	   FIM += diagm(1e-3*ones(nparam))
	 end
   FIM = PDMat(FIM)

   logp = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=calclogprior_linear)
   #logp += logprior_jitter(sigmaj)+ logprior_periods(periods)
   logp += logprior_periods(periods)
   #logp += 0.5*(log(2pi)-log(0.5*trace((Sigma \ full(dSdj))^2)))
   #println("# logp = ",logp)
   return logp
end

function compute_marginalized_likelihood(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}, minperiods::Vector{T}, maxperiods::Vector{T}; calclogprior_linear::Function = logprior_sinusoids, make_design_matrix::Function = make_design_matrix_circ  ) where T<:Real
  @assert length(minperiods)==length(maxperiods)
  function helper(periods::Vector{T})
     logp = compute_marginalized_likelihood_fix_periods(times,vels,sigma_obs,periods,calclogprior_linear=calclogprior_linear,make_design_matrix=make_design_matrix)
	 (logp == -Inf) ? 0. : exp(logp)
  end
  function helper_1d(period::T)
    println("# Evaluating P= ",period)
	flush(STDOUT)
    logp = compute_marginalized_likelihood_fix_periods(times,vels,sigma_obs,[period],calclogprior_linear=calclogprior_linear,make_design_matrix=make_design_matrix)
	println("#      ", logp)
	flush(STDOUT)
	return (logp == -Inf) ? 0. : exp(logp)
  end

  if length(minperiods)!=1
     error("# Sorry, we haven't implemented multiple periods yet")
  end
  result = QuadGK.quadgk(helper_1d,minperiods[1],maxperiods[1],reltol=1e-3)
end

compute_marginalized_likelihood(times,vels,sigma_obs,[11.9],[12.1])./log(10)

#=
compute_marginalized_likelihood(times,vels,sigma_obs,[11.9],[12.1],calclogprior_linear=logprior_epicycles,make_design_matrix=make_design_matrix_epicycle)./log(10)

compute_marginalized_likelihood(times,vels,sigma_obs,[10.],[2400.])./log(10)
=#

function logsumlogs(x::T,y::T) where T<:Real
  (x<=y) ? y + log1p(exp(x-y)) : x + log1p(exp(y-x))
end

#=
function brute_force_over_periods_1d_old(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}, Plo::Real, Phi::Real; samples_per_peak::Real = 4.0, calclogprior_linear::Function = logprior_sinusoids, make_design_matrix::Function = make_design_matrix_circ ) where T<:Real
  min_freq = 1/Phi
  max_freq = 1/Plo
  duration = maximum(times)-minimum(times)
  delta_freq = 1/(samples_per_peak * duration)
  log_delta_freq = log(delta_freq)
  range = min_freq:delta_freq:max_freq
  logintegral = -Inf
  for f in range
    P = 1/f
	logp = compute_marginalized_likelihood_fix_periods(times,vels,sigma_obs,[P], calclogprior_linear=calclogprior_linear, make_design_matrix=make_design_matrix)
	logp -= log(f)
	if logintegral == -Inf
	  logintegral = logp+log_delta_freq
	else
	  logintegral = logsumlogs(logintegral,logp+log_delta_freq)
	end
  end
  logintegral
end

function brute_force_over_periods_1d(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}, Plo::Real, Phi::Real; samples_per_peak::Real = 1.0, calclogprior_linear::Function = logprior_sinusoids, make_design_matrix::Function = make_design_matrix_circ, sigmaj::Real = -1.0 ) where T<:Real
  min_freq = 1/Phi
  max_freq = 1/Plo
  duration = maximum(times)-minimum(times)
  mean_vel = sum(vels./sigma_obs.^2)/sum(sigma_obs.^(-2))
  rms_vel = std((vels.-mean_vel)./sigma_obs)
  num_periods = convert(Int64,ceil(samples_per_peak*2pi*duration*rms_vel*(max_freq-min_freq)))
  println("# num_periods = ",num_periods)
  deltaP_const = 1/(samples_per_peak*2pi*duration*rms_vel)
  P = Plo
  periods = Array(Float64,2*num_periods)
  powers  = Array(Float64,2*num_periods)
  logintegral = -Inf
  i = 0
  while P<=Phi
    deltaP = deltaP_const*P^2
	if i >0 P += deltaP end
    i += 1
	if sigmaj < 0
	   logp = compute_marginalized_likelihood_fix_periods(times,vels,sigma_obs,[P], calclogprior_linear=calclogprior_linear, make_design_matrix=make_design_matrix)
	else
	   logp = compute_marginalized_likelihood_fix_periods_sigmaj(times,vels,sigma_obs,[P],sigmaj, calclogprior_linear=calclogprior_linear, make_design_matrix=make_design_matrix)

	   logp += logprior_jitter(sigmaj)
	end
	logp += log(deltaP/P)
    if i%10 == 0
	  println("#i=",i," P=",P," dP=",deltaP," logp=",logp," logI=",logintegral)
	end
	if logintegral == -Inf
	  logintegral = logp
	else
	  logintegral = logsumlogs(logintegral,logp)
	end
	@assert i<=length(periods)
	periods[i] = P
	powers[i] = logp

  end
  resize!(periods,i)
  resize!(powers,i)
  println("# Last P = ", P)
  return (periods,powers,logintegral)
end

res = brute_force_over_periods_1d(times,vels,sigma_obs,20.0,10000.0,samples_per_peak=1,sigmaj = 1.5)

result = QuadGK.quadgk(x->exp(brute_force_over_periods_1d(times,vels,sigma_obs,20.0,10000.0,samples_per_peak=1,sigmaj=x)),0.0,10.0,reltol=1e-2)

res = brute_force_over_periods_1d(times,vels,sigma_obs,20.0,10000.0,samples_per_peak=1)

=#

using FastGaussQuadrature
num_Sigmas = 40
nodes, weights = gausslegendre( num_Sigmas );
weights *= 0.5;
nodes = 0.5+0.5*nodes;
sigmaj_list = sigma_j_0*(exp.(nodes*log1p(sigma_j_max/sigma_j_0)).-1);
Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.));
SigmaList = map(sigmaj->Sigma0+ScalMat(length(times),sigmaj),sigmaj_list);

marginalized_likelihood_fix_periods_sigmaj = map(i->compute_marginalized_likelihood_fix_periods_sigmaj(times, vels, sigma_obs, [12.1], sigmaj_list[i], Sigma=SigmaList[i]),1:length(sigmaj_list))
int_norm = maximum(marginalized_likelihood_fix_periods_sigmaj);
integral = log(dot( weights, exp.(marginalized_likelihood_fix_periods_sigmaj-int_norm) ))+int_norm
compute_marginalized_likelihood_fix_periods(times, vels, sigma_obs, [12.1] )

function brute_force_over_periods_1d_new(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}, Plo::Real, Phi::Real; samples_per_peak::Real = 1.0, calclogprior_linear::Function = logprior_sinusoids, make_design_matrix::Function = make_design_matrix_circ, sigmaj::Real = -1.0, num_Sigmas::Integer = 40 ) where T<:Real
  min_freq = 1/Phi
  max_freq = 1/Plo
  duration = maximum(times)-minimum(times)
  mean_vel = sum(vels./sigma_obs.^2)/sum(sigma_obs.^(-2))
  rms_vel = std((vels.-mean_vel)./sigma_obs)
  num_periods = convert(Int64,ceil(samples_per_peak*2pi*duration*rms_vel*(max_freq-min_freq)))
  println("# num_periods = ",num_periods)
  periods = zeros(1)
  sigma_eps = 1e-2
  Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.+sigma_eps));
  if sigmaj >= 0.0
     Sigma = Sigma0 + ScalMat(length(times),sigmaj)
  else
     Sigma = Sigma0 + ScalMat(length(times),rms_vel)
  end

  if num_Sigmas > 1
     nodes, weights = gausslegendre( num_Sigmas );
     weights *= 0.5;
     nodes = 0.5+0.5*nodes;
     sigmaj_list = sigma_j_0*(exp.(nodes*log1p(sigma_j_max/sigma_j_0)).-1);
     SigmaList = map(s->Sigma0+ScalMat(length(times),s),sigmaj_list);
  end

  deltaP_const = 1/(samples_per_peak*2pi*duration*rms_vel)
  P = Plo
  period_grid = Array{Float64}(2*num_periods)
  power_grid  = Array{Float64}(2*num_periods)
  logintegral = -Inf
  i = 0
  while P<=Phi
    deltaP = deltaP_const*P^2
	if i >0 P += deltaP end
	periods[1] = P
    i += 1
	X = make_design_matrix(times,periods)
    if num_Sigmas >1
	   marginalized_likelihood_fix_periods_sigmaj = map(k->compute_marginalized_likelihood_fix_periods_sigmaj(times, vels, sigma_obs, periods, sigmaj_list[k], Sigma=SigmaList[k],calclogprior_linear=calclogprior_linear, make_design_matrix=make_design_matrix, X=X),1:length(sigmaj_list))
       int_norm = maximum(marginalized_likelihood_fix_periods_sigmaj);
       logp = log(dot( weights, exp.(marginalized_likelihood_fix_periods_sigmaj-int_norm) ))+int_norm

    else
	   if sigmaj < 0
	      logp = compute_marginalized_likelihood_fix_periods(times,vels,sigma_obs,periods,Sigma0=Sigma0, calclogprior_linear=calclogprior_linear, make_design_matrix=make_design_matrix,X=X)
	   else
	      logp = compute_marginalized_likelihood_fix_periods_sigmaj(times,vels,sigma_obs,periods,sigmaj, Sigma=Sigma, calclogprior_linear=calclogprior_linear, make_design_matrix=make_design_matrix,X=X)
	      logp += logprior_jitter(sigmaj)
	   end

    end
	logp += log(deltaP/P)
    if i%10 == 0
	  println("#i=",i," P=",P," dP=",deltaP," logp=",logp," logI=",logintegral)
	end
	if logintegral == -Inf
	  logintegral = logp
	else
	  logintegral = logsumlogs(logintegral,logp)
	end
	@assert i<=length(period_grid)
	period_grid[i] = P
	power_grid[i] = logp

  end
  resize!(period_grid,i)
  resize!(power_grid,i)
  println("# Last P = ", period_grid[end])
  return (period_grid,power_grid,logintegral)
end

res = brute_force_over_periods_1d_new(times,vels,sigma_obs,20.0,10000.0,samples_per_peak=1,num_Sigmas=1)

function brute_force_over_periods_multiplanet_1d_old(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}, periods_fixed::Vector{T}, Plo::Real, Phi::Real; samples_per_peak::Real = 4.0, calclogprior_linear::Function = logprior_sinusoids, make_design_matrix::Function = make_design_matrix_circ, update_design_matrix::Function = update_design_matrix_circ ) where T<:Real
  periods = zeros(length(periods_fixed)+1)
  periods[1:length(periods_fixed)] = periods_fixed
  periods[end] = (Plo+Phi)/2
  X = make_design_matrix(times,periods)
  min_freq = 1/Phi
  max_freq = 1/Plo
  duration = maximum(times)-minimum(times)
  delta_freq = 1/(samples_per_peak * duration)
  log_delta_freq = log(delta_freq)
  range = min_freq:delta_freq:max_freq
  logintegral = -Inf
  for f in range
    periods[end] =  1/f
	update_design_matrix(X,times,periods,length(periods))
	logp = compute_marginalized_likelihood_fix_periods(times,vels,sigma_obs,periods, calclogprior_linear=calclogprior_linear, X=X)
	logp -= log(f)
	if logintegral == -Inf
	  logintegral = logp+log_delta_freq
	else
	  logintegral = logsumlogs(logintegral,logp+log_delta_freq)
	end
  end
  logintegral
end

function brute_force_over_periods_multiplanet_1d(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}, periods_fixed::Vector{T}, Plo::Real, Phi::Real; samples_per_peak::Real = 2.0, calclogprior_linear::Function = logprior_sinusoids, make_design_matrix::Function = make_design_matrix_circ, update_design_matrix::Function = update_design_matrix_circ, num_Sigmas::Integer = 40 ) where T<:Real
  periods = zeros(length(periods_fixed)+1)
  periods[1:length(periods_fixed)] = periods_fixed

  # Compute how much signal is remaining so can choose appropriate grid density
  X = make_design_matrix(times,periods_fixed)
  sigmaj = calc_opt_sigmaj(times, vels, sigma_obs, X, calclogprior_linear=calclogprior_linear).minimizer
  Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.))
  Sigma = Sigma0 + ScalMat(length(times),sigmaj) # PDMat(make_Sigma(times,sigma_obs,sigmaj))
  param_bf_linalg = calc_best_fit(vels,X,Sigma)
  predict = X*param_bf_linalg

  chisq = invquad(Sigma0,vels.-predict)
  rms_vel = std((vels.-predict)./sigma_obs)
  println("# RMS resid (uncor) = ", rms_vel)
  #rms_vel_eff = sqrt(chisq/length(times))
  #println("# RMS resid (corel) = ", rms_vel_eff)

  periods[end] = (Plo+Phi)/2
  X = make_design_matrix(times,periods)

  if num_Sigmas > 1
     nodes, weights = gausslegendre( num_Sigmas );
     weights *= 0.5;
     nodes = 0.5+0.5*nodes;
     sigmaj_list = sigma_j_0*(exp.(nodes*log1p(sigma_j_max/sigma_j_0)).-1);
     Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.));
     SigmaList = map(sigmaj->Sigma0+ScalMat(length(times),sigmaj),sigmaj_list);
  end

  # Now compute periodogram
  min_freq = 1/Phi
  max_freq = 1/Plo
  duration = maximum(times)-minimum(times)
  num_periods = convert(Int64,ceil(samples_per_peak*2pi*duration*rms_vel*(max_freq-min_freq)))
  println("# num_periods = ",num_periods)
  period_grid = Array{Float64}(2*num_periods)
  power_grid  = Array{Float64}(2*num_periods)
  deltaP_const = 1/(samples_per_peak*2pi*duration*rms_vel)
  P = Plo
  logintegral = -Inf
  i = 0
  while P<=Phi
    deltaP = deltaP_const*P^2
	if i >0 P += deltaP end
    i += 1
    periods[end] =  P
	#println("# size(X) = ",size(X), " len(times)= ", length(times), " len(periods) = ", length(periods))
	update_design_matrix(X,times,periods,length(periods))

    if num_Sigmas == 1
	   #logp = compute_marginalized_likelihood_fix_periods(times,vels,sigma_obs,periods, calclogprior_linear=calclogprior_linear, X=X, Sigma0=Sigma0) # make_design_matrix=make_design_matrix)
	   logp = compute_marginalized_likelihood_fix_periods_sigmaj(times,vels,sigma_obs,periods,sigmaj, calclogprior_linear=calclogprior_linear, X=X, Sigma=Sigma) # make_design_matrix=make_design_matrix)
	else
	   marginalized_likelihood_fix_periods_sigmaj = map(k->compute_marginalized_likelihood_fix_periods_sigmaj(times, vels, sigma_obs, periods, sigmaj_list[k], Sigma=SigmaList[k],calclogprior_linear=calclogprior_linear, make_design_matrix=make_design_matrix, X=X),1:length(sigmaj_list))
       int_norm = maximum(marginalized_likelihood_fix_periods_sigmaj);
       logp = log(dot( weights, exp.(marginalized_likelihood_fix_periods_sigmaj-int_norm) ))+int_norm
	end

	logp += log(deltaP/P)
    if i%100 == 0
	  println("# i=",i," P=",P," dP=",deltaP," logp=",logp," logI=",logintegral)
	end
	if logintegral == -Inf
	  logintegral = logp
	else
	  logintegral = logsumlogs(logintegral,logp)
	end
	period_grid[i] = P
	power_grid[i] = logp

  end
  resize!(period_grid,i)
  resize!(power_grid,i)
  println("# Last P = ", P)
  (period_grid,power_grid,logintegral)
end

# Optimzie periods
function compute_marginalized_likelihood_near_periods(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}, periods::Vector{T}; calclogprior_linear::Function = logprior_sinusoids, make_design_matrix::Function = make_design_matrix_circ ) where T<:Real
   X = make_design_matrix(times,periods)

   function optimize_jitter_helper(sigmaj::Real)
     #println("# making Sigma")
	 Sigma = Sigma0 + ScalMat(length(times),sigmaj) # PDMat(make_Sigma(times,sigma_obs,sigmaj))
	 #println("# making FIM")
     FIM = PDMat(Xt_invA_X(Sigma,X))
     #println("# approximating integral")
	 evid = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=calclogprior_linear)
 	 evid += logprior_jitter(sigmaj)
	 #println("# ",evid/log(10))
     -evid
   end
   Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.))
   result = optimize(optimize_jitter_helper,sigma_j_min,sigma_j_max,rel_tol=1e-6)
   #println("# Optimize jitter: ", result)
   sigmaj = result.minimizer
   Sigma = Sigma0 + ScalMat(length(times),sigmaj) # PDMat(make_Sigma(times,sigma_obs,sigmaj))

   function optimize_period_helper(P::T)
     X = make_design_matrix(times,[P])
     FIM = PDMat(Xt_invA_X(Sigma,X))
     evid = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=calclogprior_linear)
 	 evid += logprior_jitter(sigmaj)
     -evid
   end
   function optimize_period_helper(periods::Vector{T})
     X = make_design_matrix(times,periods)
     FIM = PDMat(Xt_invA_X(Sigma,X))
     evid = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=calclogprior_linear)
 	 evid += logprior_jitter(sigmaj)
     -evid
   end

   if(length(periods)==1)
      result2 = optimize(optimize_period_helper,periods[1]*0.98,periods[1]*1.02,rel_tol=1e-6)
	  periods = [result2.minimizer]
   else
	  result2 = optimize(optimize_period_helper,periods,f_tol=1e-3) # abs(result.minimum))
	  periods = result2.minimizer
   end
   #println("# Optimize period: ", result2)
   best_period = copy(periods)
   #=
   X = make_design_matrix(times,periods)
   result = optimize(optimize_jitter_helper,sigma_j_min,sigma_j_max,rel_tol=1e-6)
   #println("# Optimize jitter: ", result)
   sigmaj = result.minimizer
   #Sigma = PDMat(make_Sigma(times,sigma_obs,sigmaj))
   Sigma = Sigma0 + ScalMat(length(times),sigmaj) # PDMat(make_Sigma(times,sigma_obs,sigmaj))
   result2 = optimize(optimize_period_helper,periods[1]*0.98,periods[1]*1.02,rel_tol=1e-6)
   println("# Optimize period: ", result2)
   best_period = result2.minimizer
   =#

   X = make_design_matrix(times,periods)
   FIM = PDMat(Xt_invA_X(Sigma,X))

   (param_bf,sigma_param,logp) = compute_Laplace_Approx_param_and_integral(vels,X,Sigma,FIM=FIM,calclogprior=calclogprior_linear)
   logp += logprior_jitter(sigmaj)+ logprior_periods(periods)

   #println("# marginal at best sigma_j & P: ",logp/log(10))

   # Laplace approx to integrate over jitter parameter
   dSdj = make_dSigmadsigmaj(times,sigma_obs,sigmaj)
   log_int_over_sigma_j = 0.5*(log(2pi)-log(0.5*trace((Sigma \ full(dSdj))^2)))
   logp += log_int_over_sigma_j

   #println("# marginal over sigma_j at best P: ",logp/log(10))

   # Laplace approx to integrate over period (one peak only)
   sigma_P = zeros(length(periods))
   log_int_over_periods = 0
   for i in 1:length(best_period)
     periods = copy(best_period)
     eps = 1e-5*periods[i]
	 periods[i] += eps

     X = make_design_matrix(times,periods)
     logp_hi = compute_Laplace_Approx(vels,X,Sigma,calclogprior=calclogprior_linear)
     logp_hi += logprior_jitter(sigmaj)+ logprior_periods(periods)
     logp_hi += log_int_over_sigma_j  # Approximates this term as const
	 periods[i] -= 2*eps
     X = make_design_matrix(times,periods)
     logp_lo = compute_Laplace_Approx(vels,X,Sigma,calclogprior=calclogprior_linear)
     logp_lo += logprior_jitter(sigmaj)+ logprior_periods(periods)
     logp_lo += log_int_over_sigma_j # Approximates this term as const
     d2logpdp2 = (logp_hi-2*logp+logp_lo)/eps^2
     sigma_P[i] = sqrt(-1/d2logpdp2)
     println("# d^2 logp/dP_",i,"^2 = ",d2logpdp2/log(10))
     log_int_over_periods += 0.5*(log(2pi)-log(abs(d2logpdp2)))
   end
   logp += log_int_over_periods # Approx as diagonal
   println("# marginal over P at ",best_period," and sigma_j at ",sigmaj," : ",logp/log(10))

   return (param_bf,sigma_param,logp,best_period,sigma_P)
end

res0 = compute_marginalized_likelihood_fix_periods(times,vels,sigma_obs,Float64[])/log(10)
res1  = compute_marginalized_likelihood_near_periods(times,vels,sigma_obs,[42.0])[3] /log(10)
res2  = compute_marginalized_likelihood_near_periods(times,vels,sigma_obs,[42.0,12.1])[3]/log(10)

res1e = compute_marginalized_likelihood_near_periods(times,vels,sigma_obs,[42.0],calclogprior_linear=logprior_epicycles,make_design_matrix=make_design_matrix_epicycle)

res2e = compute_marginalized_likelihood_near_periods(times,vels,sigma_obs,[42.0,12.0],calclogprior_linear=logprior_epicycles,make_design_matrix=make_design_matrix_epicycle)

res2b =  brute_force_over_periods_multiplanet_1d(times,vels,sigma_obs,[42.0777],10.0,2400.0,samples_per_peak=8)

function compute_marginalized_likelihood_cuba(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}, minperiods::Vector{T}, maxperiods::Vector{T}; calclogprior_linear::Function = logprior_sinusoids, make_design_matrix::Function = make_design_matrix_circ  ) where T<:Real
  @assert length(minperiods)==length(maxperiods)
  log_normalization = 0.0
  function helper_1d(period::T) # Note, this exp's unlike multi-d helper
    #println("# Evaluating P= ",period)
	#flush(STDOUT)
    logp = compute_marginalized_likelihood_fix_periods(times,vels,sigma_obs,[period],calclogprior_linear=calclogprior_linear,make_design_matrix=make_design_matrix)
	#println("#      ", logp)
	#flush(STDOUT)
	return (logp == -Inf) ? 0. : exp(logp-log_normalization)
  end
  function helper(periods::Vector{T})
    #println("# Evaluating P= ",periods)
	#flush(STDOUT)
     logp = compute_marginalized_likelihood_fix_periods(times,vels,sigma_obs,periods,calclogprior_linear=calclogprior_linear,make_design_matrix=make_design_matrix)
	#println("#      ", logp)
	#flush(STDOUT)
	return logp
	# (logp == -Inf) ? 0. : exp(logp-log_normalization)
  end
  function integrand_cuba(x::Vector, output::Vector)
    @assert length(output) == 1
	#println("# x=",x)
    periods = (maxperiods.-minperiods).*x.+minperiods
    output[1] = exp(helper(periods)-log_normalization)
	println("# P= ",periods," logp= ",output[1])
  end

  log_normalization = helper(0.5*(minperiods.+maxperiods))
  if length(minperiods)==1
     result = QuadGK.quadgk(helper_1d,minperiods[1],maxperiods[1],reltol=1e-3)
	 result = log(result[1])
     #result = divonne(integrand_cuba,1,1,reltol=1e-3,maxevals=100)*prod(maxperiods.-minperiods) # Problem?
  else

     result = divonne(integrand_cuba,length(minperiods),1,reltol=1e-3,maxevals=1000)
	 result = log(result.integral[1]) + log(prod(maxperiods.-minperiods))
  end
  result+log_normalization
end

res = compute_marginalized_likelihood_cuba(times,vels,sigma_obs,[41.9,11.9],[42.3,12.2])

#=
@time vegas(integrand_cuba,2,1,reltol=1e-3,maxevals=1e5)
@time suave(integrand_cuba,2,1,reltol=1e-3,maxevals=1e5)
@time divonne(integrand_cuba,2,1,reltol=1e-3,maxevals=1e5)
@time cuhre(integrand_cuba,2,1,reltol=1e-3,minevals=1e3,maxevals=1e5)
=#

#=
# Brute force 2-d integration
using Cubature

function log_target(x::Vector)
  offset = x[1]
  sigma_j = x[2]
  logprior = logprior_offset([offset]) + logprior_jitter(sigma_j)
  Sigma = make_Sigma(times,sigma_obs,sigma_j)
  delta = vels.-offset
  chisq = invquad(Sigma,delta)
  loglikelihood = 0.5*(-chisq -logdet(Sigma)-nobs*log(2pi))
  logprior + loglikelihood
end

function integrand(x::Vector; log_normalization::Real=0.0)
  exp(log_target(x)-log_normalization)
end
@time hcubature(integrand, [-Cmax/100,sigma_j_min],[Cmax/100,sigma_j_max],reltol=1e-2)
#@time pcubature(integrand, [-Cmax/100,sigma_j_min],[Cmax/100,sigma_j_max],reltol=1e-2)
=#
#=
using Cuba
function integrand_cuba(x::Vector, output::Vector; log_normalization::Real=0.0)
  xx = zeros(2)
  xx[1] = x[1]*2*Cmax/100-Cmax/100
  xx[2] = x[2]*(sigma_j_max-sigma_j_min)+sigma_j_min
  output[1] = exp(log_target(xx)-log_normalization)*(2*Cmax/100)*(sigma_j_max-sigma_j_min)
end
@time vegas(integrand_cuba,2,1,reltol=1e-3,maxevals=1e5)
@time suave(integrand_cuba,2,1,reltol=1e-3,maxevals=1e5)
@time divonne(integrand_cuba,2,1,reltol=1e-3,maxevals=1e5)
@time cuhre(integrand_cuba,2,1,reltol=1e-3,minevals=1e3,maxevals=1e5)
=#
#=
# Brute Force:  Sequence of Nested 1-d integrations
function compute_log_target_fixed_offset(offset::Real)
  delta = vels.-offset
  function integrand(sigma_j::Real)
    logprior = logprior_offset([offset]) + logprior_jitter(sigma_j)
    Sigma = Sigma0 + ScalMat(length(times),sigmaj) # PDMat(make_Sigma(times,sigma_obs,sigmaj))
    chisq = invquad(Sigma,delta)
    loglikelihood = 0.5*(-chisq -logdet(Sigma)-nobs*log(2pi))
    logtarget = logprior + loglikelihood
	exp(logtarget) # +log(1e210))
  end
  Sigma0 = make_Sigma(times,sigma_obs,0.)
  result = QuadGK.quadgk(integrand,sigma_j_min,sigma_j_0,sigma_j_max,reltol=1e-2)[1]
end
@time result = QuadGK.quadgk(x->begin p=compute_log_target_fixed_offset(x); println(x," ",p); p end,-Cmax/100,Cmax/100,reltol=1e-2)./log(10)

# Brute Force: 2-d integration
function log_target(x::Vector)
  offset = x[1]
  sigma_j = x[2]
  logprior = logprior_offset([offset]) + logprior_jitter(sigma_j)
  Sigma = make_Sigma(times,sigma_obs,sigma_j)
  delta = vels.-offset
  chisq = invquad(Sigma,delta)
  loglikelihood = 0.5*(-chisq -logdet(Sigma)-nobs*log(2pi))
  logprior + loglikelihood
end
=#

#=
# Fit one planet model
periods = [42.0]
nparam = 1+2*length(periods)
freqs = 2pi./periods
X = [sin(freqs[1]*times) cos(freqs[1]*times) ones(times)];
Sigma = PDMat(make_Sigma(times,sigma_obs,sigma_j_0))
FIM = PDMat(Xt_invA_X(Sigma,X))
(param_bf,sigma_param,evid) = compute_Laplace_Approx_param_and_integral(vels,X,Sigma,FIM=FIM,calclogprior=logprior_sinusoids)

function calc_opt_sigmaj(data::Vector{T}, model::Array{T,2}) where T<:Real
  function calc_laplace_approx_arg_jitter(sigmaj::Real)
	Sigma = Sigma0 + ScalMat(length(times),sigmaj) # PDMat(make_Sigma(times,sigma_obs,sigmaj))
    FIM = PDMat(Xt_invA_X(Sigma,X))
    evid = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=logprior_sinusoids)
    evid += logprior_jitter(sigmaj)
	-evid
  end
  Sigma0 = make_Sigma(times,sigma_obs,0.)
  result = optimize(calc_laplace_approx_arg_jitter,sigma_j_min,sigma_j_max,rel_tol=1e-4)
end
res = calc_opt_sigmaj(vels,X)
Sigma = PDMat(make_Sigma(times,sigma_obs,res.minimizer))
dSigmadsigmaj = make_dSigmadsigmaj(times,sigma_obs,res.minimizer)
FIM = PDMat(Xt_invA_X(Sigma,X))
(param_bf,sigma_param,evid) = compute_Laplace_Approx_param_and_integral(vels,X,Sigma,FIM=FIM,calclogprior=logprior_sinusoids)
evid += 0.5*(log(2pi)-log(0.5*trace((Sigma \ full(dSdj))^2)))

Ks = generate_K_samples(param_bf,sigma_param,1000);
m = mean(Ks,1);
sig = std(Ks,1,mean=m);
m-sig,m,m+sig

function profile_sigmaj(data::Vector{T}, model::Array{T,2}) where T<:Real
  function calc_laplace_approx_arg_jitter(sigmaj::Real) #param::Vector{T})
    #sigmaj = param[1]
	Sigma = PDMat(make_Sigma(times,sigma_obs,sigmaj))
    FIM = PDMat(Xt_invA_X(Sigma,X))
    evid = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=logprior_sinusoids)
	evid += logprior_jitter(sigmaj)
    evid
  end
  #result = optimize(calc_laplace_approx_arg_jitter,[sigma_j_0])
  result = optimize(calc_laplace_approx_arg_jitter,sigma_j_min,sigma_j_max,rel_tol=1e-6)
end
res = calc_opt_sigmaj(vels,X)

# Fit two planet model
periods = [42.0,12.0]
X = make_design_matrix(times,periods)
Sigma = PDMat(make_Sigma(times,sigma_obs,0.68442348767604010984))
FIM = PDMat(Xt_invA_X(Sigma,X))
(param_bf,sigma_param,evid) = compute_Laplace_Approx_param_and_integral(vels,X,Sigma,FIM=FIM,calclogprior=logprior_sinusoids)

Ks = generate_K_samples(param_bf,sigma_param,1000)
m = mean(Ks,1)
sig = std(Ks,1,mean=m)
m-sig,m+sig

=#


text:
"""
    PySet{T=Py}([x])

Wraps the Python set `x` (or anything satisfying the set interface) as an `AbstractSet{T}`.

If `x` is not a Python object, it is converted to one using `pyset`.
"""
struct PySet{T} <: AbstractSet{T}
    py :: Py
    PySet{T}(::Val{:new}, py::Py) where {T} = new{T}(py)
end
export PySet

PySet{T}(x=pyset()) where {T} = PySet{T}(Val(:new), ispy(x) ? Py(x) : pyset(x))
PySet(x=pyset()) = PySet{Py}(x)

ispy(::PySet) = true
getpy(x::PySet) = x.py
pydel!(x::PySet) = pydel!(x.py)

pyconvert_rule_set(::Type{T}, x::Py, ::Type{PySet{V}}=Utils._type_ub(T)) where {T<:PySet,V} =
    if PySet{Py} <: T
        pyconvert_return(PySet{Py}(x))
    else
        pyconvert_return(PySet{V}(x))
    end

Base.length(x::PySet) = Int(pylen(x))

Base.isempty(x::PySet) = length(x) == 0

function Base.iterate(x::PySet{T}, it::Py=pyiter(x)) where {T}
    y = pynext(it)
    if ispynull(y)
        pydel!(it)
        return nothing
    else
        return (pyconvert_and_del(T, y), it)
    end
end

function Base.in(v, x::PySet{T}) where {T}
    if v isa T
        return pycontains(x, v)
    else
        r = pyconvert_tryconvert(T, v)
        if pyconvert_isunconverted(r)
            return false
        else
            return pycontains(x, pyconvert_result(T, r))
        end
    end
end

function Base.push!(x::PySet{T}, v) where {T}
    pydel!(@py x.add(@jl(convert(T, v))))
    return x
end

function Base.delete!(x::PySet{T}, v) where {T}
    if v isa T
        pydel!(@py x.discard(v))
    else
        r = pyconvert_tryconvert(T, v)
        if !pyconvert_isunconverted(r)
            pydel!(@py x.discard(@jl pyconvert_result(T, r)))
        end
    end
    return x
end

Base.@propagate_inbounds function Base.pop!(x::PySet{T}) where {T}
    @boundscheck (isempty(x) && throw(ArgumentError("set must be non-empty")))
    return pyconvert_and_del(T, @py x.pop())
end

function Base.pop!(x::PySet, v)
    if v in x
        delete!(x, v)
        return v
    else
        throw(KeyError(v))
    end
end

function Base.pop!(x::PySet, v, d)
    if v in x
        delete!(x, v)
        return v
    else
        return d
    end
end

function Base.empty!(x::PySet)
    pydel!(@py x.clear())
    return x
end

function Base.copy(x::PySet{T}) where {T}
    return PySet{T}(Val(:new), @py x.copy())
end


text:
using Base.Test

# using Lint
# lintpkg( "Discretizers" )
# @test isempty(lintpkg( "Discretizers", returnMsgs=true))

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

using Discretizers

function array_matches{S<:Integer, T<:Integer}(arr::AbstractVector{S}, arr2::AbstractVector{T})
    n = length(arr)
    @assert(length(arr2) == n)
    for i = 1 : n
        if arr[i] != arr2[i]
            return false
        end
    end
    true
end
function array_matches{S<:Real, T<:Real}(arr::AbstractVector{S}, arr2::AbstractVector{T}, abs_tolerance::AbstractFloat=eps(Float64))
    n = length(arr)
    @assert(length(arr2) == n)
    for i = 1 : n
        if abs(arr[i] - arr2[i]) > abs_tolerance
            return false
        end
    end
    true
end

include("test_categorical_discretizer.jl")
include("test_linear_discretizer.jl")
include("test_hybrid_discretizer.jl")

include("test_common.jl")

include("test_disc_uniformwidth.jl")
include("test_disc_uniformcount.jl")
include("test_disc_bayesianblocks.jl")
include("test_disc_MODL.jl")