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

include("../src/rv_model_keplerian_simple.jl")

# 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
   if any(P.<=zero(P)) 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 make_fischer_information_matrix(model::Array{T,2}, Sigma::SigmaT) where {T<:Real, SigmaT<:AbstractPDMat{T}}
   FIM = Xt_invA_X(Sigma,model)
   make_matrix_posdef(FIM)
end

function calc_chisq(data::Vector{T}, model::Array{T,2}, Sigma::SigmaT; FIM::PDMat{T} = make_fischer_information_matrix(model,Sigma), 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} = make_fischer_information_matrix(model,Sigma), calclogprior::Function = x->(0.0,0.0), report_log_MAP::Bool = false ) 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]
  if report_log_MAP     return (param_bf_linalg, sigma_param,log_mode)  end
  LaplaceApprox = 0.5*nparam*log(2pi)-0.5*logdet(FIM)+log_mode + log_jacobian
  return (param_bf_linalg,sigma_param,LaplaceApprox)
end

function compute_max_aposterior_Approx(data::Vector{T}, model::Array{T,2}, Sigma::SigmaT; 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
  return (log_mode, logprior_mode, loglikelihood_mode, log_jacobian)
end

function compute_Laplace_Approx(data::Vector{T}, model::Array{T,2}, Sigma::SigmaT; FIM::PDMat{T} = make_fischer_information_matrix(model,Sigma), calclogprior::Function = x->(0.0,0.0), report_log_MAP::Bool = false ) where {T<:Real, SigmaT<:AbstractPDMat{T} }
  @assert length(data) == size(model,1) == size(Sigma,1)
  nparam = size(model,2)
  nobs = length(data)
#=
  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
=#
  (log_mode, logprior_mode, loglikelihood_mode, log_jacobian) = compute_max_aposterior_Approx(data,model,Sigma,calclogprior=calclogprior)
  if report_log_MAP return log_mode  end

  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)
  j = 1
  while j<=nsamples
    theta = bf .+ sigma*randn(nparam)
	for i in 1:npl
       Ks[j,i] = sqrt.(theta[1+namp_per_pl*(i-1)].^2+theta[2+namp_per_pl*(i-1)].^2)
	   hs[j,i] = theta[3+namp_per_pl*(i-1)] ./ Ks[j,i]
	   ks[j,i] = theta[4+namp_per_pl*(i-1)] ./ Ks[j,i]
	   es[j,i] = sqrt.(hs[j,i].^2.+ks[j,i].^2)
	end
	if all(es[j,:].<1.0)
	     j +=1
	end
  end
  return (Ks,es,hs,ks)
end

function generate_frac_valid_e(bf::Vector, sigma, nsamples::Integer = 100)
  nparam = length(bf)
  const namp_per_pl = 4
  npl = convert(Int64,floor((nparam-1)//namp_per_pl))
  es = generate_Ke_samples(bf,sigma,nsamples)[2]
  num_pass = 0
  for i in 1:size(es,1)
     if all(es[i,:].<1.0) num_pass += 1 end
  end
  num_pass/nsamples
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 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)
    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

function calc_map_no_planet(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}) where T<:Real
  function calc_map_fixed_jitter(sigmaj::Real)
    Sigma = Sigma0 + ScalMat(length(times),sigmaj)
    evid = compute_Laplace_Approx(vels,model,Sigma,report_log_MAP=true)
	evid += logprior_jitter(sigmaj)
    -evid
  end
  model = ones(length(times),1)
  Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.))
  result = optimize(calc_map_fixed_jitter,sigma_j_min,sigma_j_max,rel_tol=1e-6)
  return -result.minimum
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)), report_log_MAP::Bool=false ) 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 = make_fischer_information_matrix(X,Sigma)
     evid = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=calclogprior_linear,report_log_MAP=report_log_MAP)
 	 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 = make_fischer_information_matrix(X,Sigma)
   logp = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=calclogprior_linear,report_log_MAP=report_log_MAP)
   logp += logprior_jitter(sigmaj)+ logprior_periods(periods)
   if !report_log_MAP
      logp += 0.5*(log(2pi)-log(0.5*trace((Sigma \ full(dSdj))^2)))
   end
   #println("# logp = ",logp)
   return logp
end

function make_Sigma_cache(num_Sigmas::Integer, Sigma0::SigmaT ) where {T<:Real, SigmaT<:AbstractPDMat{T}}
   @assert 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(size(Sigma0,1),s),sigmaj_list);
   return (sigmaj_list, weights, SigmaList)
end

function make_Sigma_cache(times::Vector{T}, sigma_obs::Vector{T}, num_Sigmas::Integer ) where T<:Real
   @assert num_Sigmas >= 1
   Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.0))
   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);
   return (sigmaj_list, weights, SigmaList)
end

function compute_marginalized_likelihood_const_sigmalist(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}; Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.0)), num_Sigmas::Integer = 20, Sigma_cache = make_Sigma_cache(num_Sigmas,Sigma0) ) where T<:Real
   X = ones(length(times),1)
   const nparam = size(X,2)
   sigmaj_list = Sigma_cache[1]
   weights = Sigma_cache[2]
   SigmaList = Sigma_cache[3]

	   marginalized_likelihood_fix_periods_sigmaj = map(k->compute_marginalized_likelihood_fix_periods_sigmaj(times, vels, sigma_obs, Float64[], sigmaj_list[k], Sigma=SigmaList[k],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

function compute_marginalized_likelihood_fix_periods_sigmalist(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)), num_Sigmas::Integer = 20, Sigma_cache = make_Sigma_cache(num_Sigmas,Sigma0) ) where T<:Real
   const nparam = size(X,2)
   sigmaj_list = Sigma_cache[1]
   weights = Sigma_cache[2]
   SigmaList = Sigma_cache[3]

	   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

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), FIM = make_fischer_information_matrix(X,Sigma) ) 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)

   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

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(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 = 20, verbose::Integer = 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)
  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);
	 =#
	 (sigmaj_list,weights,SigmaList) = make_Sigma_cache(num_Sigmas,Sigma0)
  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
	   logp = compute_marginalized_likelihood_fix_periods_sigmalist(times, vels, sigma_obs, periods, calclogprior_linear=calclogprior_linear, make_design_matrix=make_design_matrix, X=X, Sigma0=Sigma0,  num_Sigmas=num_Sigmas, Sigma_cache=(sigmaj_list,weights,SigmaList) )
	#=
	   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 verbose >=2 && 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
	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

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, sigmaj::Real = -1.0, num_Sigmas::Integer = 20, verbose::Integer = 0 ) 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_opt = 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_opt) # PDMat(make_Sigma(times,sigma_obs,sigmaj_opt))
  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);
	 =#
	 (sigmaj_list,weights,SigmaList) = make_Sigma_cache(num_Sigmas,Sigma0)
  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
	   if sigmaj < 0
   	      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)
	   else
	      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)
	   end
    else
		  logp = compute_marginalized_likelihood_fix_periods_sigmalist(times, vels, sigma_obs, periods, calclogprior_linear=calclogprior_linear, make_design_matrix=make_design_matrix, X=X, Sigma0=Sigma0,  num_Sigmas=num_Sigmas, Sigma_cache=(sigmaj_list,weights,SigmaList) )
       #=
	   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 verbose >=2 && 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
	period_grid[i] = P
	power_grid[i] = logp

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

# Optimize periods
function compute_marginalized_likelihood_near_periods_laplace(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, report_log_MAP::Bool=false ) 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,report_log_MAP=report_log_MAP)
 	 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,report_log_MAP=report_log_MAP)
 	 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,report_log_MAP=report_log_MAP)
 	 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)
   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,report_log_MAP=report_log_MAP)
   logp += logprior_jitter(sigmaj)+ logprior_periods(periods)

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

   if !report_log_MAP
      # 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))
   else
      log_int_over_sigma_j = 0.
   end

   # Laplace approx to integrate over period
   sigma_P = zeros(length(periods))
   log_int_over_periods = 0
   duration = maximum(times)-minimum(times)
   param_per_pl = (calclogprior_linear == logprior_sinusoids) ? 2 : 4
   for i in 1:length(best_period)
     periods = copy(best_period)

     amp = sqrt(param_bf[1+(i-1)*param_per_pl]^2+param_bf[2+(i-1)*param_per_pl]^2)
     eps = periods[i]/(40*2pi*duration*amp)
	 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(abs(1/d2logpdp2))
     #println("# d^2 logp/dP_",i,"^2 = ",d2logpdp2/log(10))
     log_int_over_periods += 0.5*(log(2pi)-log(abs(d2logpdp2)))
   end

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

   logp += log_int_over_periods # Approx as diagonal
   println("# marginal over P at ",best_period," and sigma_j at ",sigmaj," (laplace): ",logp/log(10))
   if calclogprior_linear == logprior_epicycles
      frac_valid_e = generate_frac_valid_e(param_bf,sigma_param)
      println("# frac w/ valid e = ", frac_valid_e)
      logp += log(frac_valid_e)
   end
   return (param_bf,sigma_param,logp,best_period,sigma_P)
end

function compute_marginalized_likelihood_near_periods_gauss(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, num_Sigmas::Integer = 20 ) 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-4)
	  periods = [result2.minimizer]
   else
	  result2 = optimize(optimize_period_helper,periods,f_tol=1e-4) # abs(result.minimum))
	  periods = result2.minimizer
   end
   #println("# Optimize period: ", result2)
   best_period = copy(periods)
   X = make_design_matrix(times,periods)
  (param_bf,sigma_param,logp) = compute_Laplace_Approx_param_and_integral(vels,X,Sigma,calclogprior=calclogprior_linear)

    (sigmaj_list,weights,SigmaList) = make_Sigma_cache(num_Sigmas,Sigma0)
     logp = compute_marginalized_likelihood_fix_periods_sigmalist(times, vels, sigma_obs, periods, calclogprior_linear=calclogprior_linear, make_design_matrix=make_design_matrix, X=X, Sigma0=Sigma0,  num_Sigmas=num_Sigmas, Sigma_cache=(sigmaj_list,weights,SigmaList) )

   # Laplace approx to integrate over period
   sigma_P = zeros(length(periods))
   log_int_over_periods = 0
   duration = maximum(times)-minimum(times)
   param_per_pl = (calclogprior_linear == logprior_sinusoids) ? 2 : 4
   for i in 1:length(best_period)
     periods = copy(best_period)
	 amp = sqrt(param_bf[1+(i-1)*param_per_pl]^2+param_bf[2+(i-1)*param_per_pl]^2)
     eps = periods[i]/(40*2pi*duration*amp)
	 periods[i] += eps

     X = make_design_matrix(times,periods)
	 logp_hi = compute_marginalized_likelihood_fix_periods_sigmalist(times, vels, sigma_obs, periods, calclogprior_linear=calclogprior_linear, make_design_matrix=make_design_matrix, X=X, Sigma0=Sigma0,  num_Sigmas=num_Sigmas, Sigma_cache=(sigmaj_list,weights,SigmaList) )

	 periods[i] -= 2*eps
     X = make_design_matrix(times,periods)
	 logp_lo = compute_marginalized_likelihood_fix_periods_sigmalist(times, vels, sigma_obs, periods, calclogprior_linear=calclogprior_linear, make_design_matrix=make_design_matrix, X=X, Sigma0=Sigma0,  num_Sigmas=num_Sigmas, Sigma_cache=(sigmaj_list,weights,SigmaList) )

     d2logpdp2 = (logp_hi-2*logp+logp_lo)/eps^2
     sigma_P[i] = sqrt(abs(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," (gauss): ",logp/log(10))
   if calclogprior_linear == logprior_epicycles
      frac_valid_e = generate_frac_valid_e(param_bf,sigma_param)
      println("# frac w/ valid e = ", frac_valid_e)
      logp += log(frac_valid_e)
   end
   return (param_bf,sigma_param,logp,best_period,sigma_P)
end

function find_n_peaks_in_pgram(period::Vector{Float64}, power::Vector{Float64}; num_peaks::Integer = 1, exclude_period_factor::Real = 1.2, exclude_peaks::Vector{Float64} = Float64[] )
  @assert num_peaks>= 1
  peak_periods = zeros(num_peaks)
  peak_logps = zeros(num_peaks)
  idx_active = trues(length(period))
  for i in 1:length(exclude_peaks)
     for j in 1:length(period)
        if exclude_peaks[i]/exclude_period_factor <= period[j] <= exclude_peaks[i]*exclude_period_factor
           idx_active[j] = false
        end # if near an excluded peak
     end # for over periods
  end # for over excluded peaks
  peaks_found = 1
  idx_peak = findmax(power[idx_active])[2]
  peak_periods[peaks_found] = period[idx_active][idx_peak]
  peak_logps[peaks_found] = power[idx_active][idx_peak]
  while peaks_found < num_peaks
     for j in 1:length(period)
        if peak_periods[peaks_found]/exclude_period_factor <= period[j] <= peak_periods[peaks_found]*exclude_period_factor
           idx_active[j] = false
        end # blackout periods near most recent peak
     end # for over periods

     peaks_found += 1
	 idx_peak = findmax(power[idx_active])[2]
     peak_periods[peaks_found] = period[idx_active][idx_peak]
	 peak_logps[peaks_found] = power[idx_active][idx_peak]

	 end # while more peaks to be found
  (peak_periods,peak_logps)
end

function compute_evidences_laplace_circ{T}(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}; num_Sigmas::Integer = 20, samples_per_peak::Integer = 1, Pminsearch::Vector{Float64} = 10.0*ones(3), report_log_MAP::Bool=false, max_num_pl::Integer=3 )
	 #(sigmaj_list,weights,SigmaList) = make_Sigma_cache(num_Sigmas,Sigma0)

runtimes = zeros(max_num_pl+1)
tic()
log10_evidences = zeros(max_num_pl+1)
res0l = compute_marginalized_likelihood_fix_periods(times,vels,sigma_obs,Float64[],report_log_MAP=report_log_MAP)
res0g = compute_marginalized_likelihood_const_sigmalist(times,vels,sigma_obs,num_Sigmas=num_Sigmas)
println("# Evidence for npl=0: ",res0l/log(10),", ", res0g/log(10))
if report_log_MAP
  log10_evidences[1] = res0l/log(10)
else
  log10_evidences[1] = res0g/log(10)
end
runtimes[1] = toq()
if max_num_pl<=0
  return (log10_evidences,[],[res0l],runtimes)
end

tic()
@time res1b = brute_force_over_periods_1d(times,vels,sigma_obs,Pminsearch[1],Pmax,samples_per_peak=samples_per_peak,num_Sigmas=num_Sigmas)
P1 = find_n_peaks_in_pgram(res1b[1],res1b[2])[1][1]

res1l = compute_marginalized_likelihood_near_periods_laplace(times,vels,sigma_obs, [P1],report_log_MAP=report_log_MAP )
res1g = compute_marginalized_likelihood_near_periods_gauss(times,vels,sigma_obs, [P1],num_Sigmas=num_Sigmas )
if report_log_MAP
   log10_evidences[2] = res1l[3]/log(10)
else
   log10_evidences[2] = res1g[3] > res1b[3] ? res1g[3]/log(10) : res1b[3]/log(10)
end
P1 = res1l[4][1]
runtimes[2] = toq()

if max_num_pl<=1
  return (log10_evidences,[P1],[res1l],runtimes)
end

tic()
@time res2b = brute_force_over_periods_multiplanet_1d(times,vels,sigma_obs,[P1],Pminsearch[2],Pmax,samples_per_peak=samples_per_peak,num_Sigmas=num_Sigmas)
println("# About to search for second planet, excluding peak near P=",P1)
P2 = find_n_peaks_in_pgram(res2b[1],res2b[2],exclude_peaks=[P1])[1][1]

res2l = compute_marginalized_likelihood_near_periods_laplace(times,vels,sigma_obs, [P1,P2],report_log_MAP=report_log_MAP )
res2g = compute_marginalized_likelihood_near_periods_gauss(times,vels,sigma_obs, [P1,P2],num_Sigmas=num_Sigmas  )
if report_log_MAP
   log10_evidences[3] = res2l[3]/log(10)
else
   log10_evidences[3] = res2g[3] > res2b[3] ? res2g[3]/log(10) : res2b[3]/log(10)
end
P1 = res2l[4][1]
P2 = res2l[4][2]
runtimes[3] = toq()

if max_num_pl<=2
  return (log10_evidences,[P1,P2],[res1l,res2l],runtimes)
end

tic()
@time res3b = brute_force_over_periods_multiplanet_1d(times,vels,sigma_obs,[P1,P2],Pminsearch[3],Pmax,samples_per_peak=samples_per_peak,num_Sigmas=num_Sigmas)
P3 = find_n_peaks_in_pgram(res3b[1],res3b[2],num_peaks=10,exclude_peaks=[P1,P2])[1][1]

res3l = compute_marginalized_likelihood_near_periods_laplace(times,vels,sigma_obs, [P1,P2,P3],report_log_MAP=report_log_MAP )
res3g = compute_marginalized_likelihood_near_periods_gauss(times,vels,sigma_obs, [P1,P2,P3],num_Sigmas=num_Sigmas  )
if report_log_MAP
   log10_evidences[4] = res3l[3]/log(10)
else
   log10_evidences[4] = res3g[3] > res3b[3] ? res3g[3]/log(10) : res3b[3]/log(10)
end
runtimes[4] = toq()

  return (log10_evidences,[P1,P2,P3],[res1l,res2l,res3l],runtimes)
end

function make_kepler_param_from_linear_output(output, npl::Integer)
  evidence_linear = output[1][npl]
  periods = output[2][1:npl]
  param_linear = output[4][npl]
  sigmaj = output[6][npl]
  chisq_linear = output[10][npl]
  const nparam_per_pl = 5
  const nparam_non_pl = 2
  const nparam_linear_per_pl = 4
  const nparam = npl*nparam_per_pl + nparam_non_pl
  const tref = 300.0
  param = Array{Float64}(nparam)
  for i in 1:npl
      P = periods[i]
	  Ksinl = param_linear[1+(i-1)*nparam_linear_per_pl]
	  Kcosl = param_linear[2+(i-1)*nparam_linear_per_pl]
	  K = sqrt(Ksinl^2+Kcosl^2)
	  wpM = atan2(Ksinl,Kcosl)
	  h = K>0 ? param_linear[3+(i-1)*nparam_linear_per_pl]/K : 0
	  k = K>0 ? param_linear[4+(i-1)*nparam_linear_per_pl]/K : 0
	  w = atan2(h,k)
      param[1+(i-1)*nparam_per_pl] = P
	  param[2+(i-1)*nparam_per_pl] = K > 0 ? K : 0.5
	  param[3+(i-1)*nparam_per_pl] = 1e-8 # h/100
	  param[4+(i-1)*nparam_per_pl] = 1e-4 # k/100
	  param[5+(i-1)*nparam_per_pl] = mod2pi(wpM+2pi*tref/P)
  end
  param[1+npl*nparam_per_pl] = param_linear[1+npl*nparam_linear_per_pl]
  param[2+npl*nparam_per_pl] = sigmaj
  param
end

function optimize_phase{T}(param, times::Vector{T}, vels::Vector{T}, Sigma, plid::Integer)
  const nparam_per_pl = 5
  const nparam_non_pl = 2
  npl = convert(Int64,floor((length(param)-1)//nparam_per_pl))
  function optimize_phase_helper(x::Real)
    param[5+(plid-1)*nparam_per_pl] = x
    -loglikelihood_fixed_jitter(param,times,vels,Sigma)
  end
  result = optimize(optimize_phase_helper,0.0,2pi)
  param[5+(plid-1)*nparam_per_pl] = mod2pi(result.minimizer)
  param
end

function optimize_omega{T}(param, times::Vector{T}, vels::Vector{T}, Sigma, plid::Integer)
  const nparam_per_pl = 5
  const nparam_non_pl = 2
  npl = convert(Int64,floor((length(param)-1)//nparam_per_pl))
  function optimize_omega_helper(x::Vector)
    param[3+(plid-1)*nparam_per_pl:4+(plid-1)*nparam_per_pl] = x
    -loglikelihood_fixed_jitter(param,times,vels,Sigma)
  end
  result = optimize(optimize_omega_helper, param[3+(plid-1)*nparam_per_pl:4+(plid-1)*nparam_per_pl] )
  if sum(result.minimizer.^2)<1
     param[3+(plid-1)*nparam_per_pl:4+(plid-1)*nparam_per_pl] = result.minimizer
  else
     param[3+(plid-1)*nparam_per_pl] = 0.0 #1e-6*param[3+(plid-1)*nparam_per_pl]
	 param[4+(plid-1)*nparam_per_pl] = 0.0 #1e-6*param[4+(plid-1)*nparam_per_pl]
  end
  param
end

function optimize_phases{T}(param, times::Vector{T}, vels::Vector{T}, Sigma)
  const nparam_per_pl = 5
  const nparam_non_pl = 2
  npl = convert(Int64,floor((length(param)-1)//nparam_per_pl))
  for plid in 1:npl
     optimize_phase(param,times,vels,Sigma,plid)
	 println("# param (after phase ",plid," ) = ", param)
	 #optimize_omega(param,times,vels,Sigma,plid)
	 #println("# param = (after omega ",plid," = ", param)
  end
  param
end

function analyze_dataset(filename::String; num_Sigmas::Integer = 20, samples_per_peak::Real = 1.0, Pminsearch::Vector{Float64} = 10.0*ones(3))
  # Read dataset
  data = readdlm(filename);
  times = data[:,1];
  vels = data[:,2];
  sigma_obs = data[:,3];
  nobs = length(vels);

   (evidence_circ_list, P_best, result_laplace_circ_list, runtimes_circ) = compute_evidences_laplace_circ(times,vels,sigma_obs,num_Sigmas=num_Sigmas, Pminsearch=Pminsearch)
   runtimes_kepler = zeros(runtimes_circ)

   param_bf_list = Vector{Vector{Float64}}(length(P_best))
   resid_list = Vector{Vector{Float64}}(length(P_best))
   rms_list = Vector{Float64}(length(P_best))
   chisq_list = Vector{Float64}(length(P_best))
   sigmaj_list = Vector{Float64}(length(P_best))
   sigma_P_list = Vector{Vector{Float64}}(length(P_best))
   sigma_sigmaj_list = Vector{Float64}(length(P_best))
   sigma_param_list = Vector(length(P_best))
   frac_valid_list = Vector{Float64}(length(P_best))
   evidence_kepler_list = Vector{Float64}(length(P_best))
   kepler_fit_list = Vector(length(P_best))
   for np in 1:length(P_best)  # Loop over number of planets to include
      sigma_P_list[np] = result_laplace_circ_list[3][5]
	  tic()
	  X = make_design_matrix_epicycle(times,P_best[1:np])
      sigmaj_opt = calc_opt_sigmaj(times, vels, sigma_obs, X, calclogprior_linear=logprior_epicycles).minimizer
      Sigma = PDMat(make_Sigma(times,sigma_obs,sigmaj_opt))
      param_bf_list[np] = calc_best_fit(vels,X,Sigma)
      predict = X*param_bf_list[np]
      # Estimate Uncertainties
      FIM = make_fischer_information_matrix(X,Sigma)
      sigma_param = inv(FIM).chol[:U]
      frac_valid_e = generate_frac_valid_e(param_bf_list[np],sigma_param)
      # Uncertainty on sigmaj
      dSdj = make_dSigmadsigmaj(times,sigma_obs,sigmaj_opt)
      sigma_sigmaj = 1/sqrt(abs(0.5*trace((Sigma \ full(dSdj))^2)))
      println("# np=", np," sigma_sigmaj = ",sigma_sigmaj,"  frac w/ valid e = ", frac_valid_e)
      evidence_kepler_list[np] = compute_marginalized_likelihood_near_periods_gauss(times,vels,sigma_obs, P_best[1:np],num_Sigmas=num_Sigmas, calclogprior_linear=logprior_epicycles, make_design_matrix=make_design_matrix_epicycle )[3]/log(10)

      resid_list[np] = vels.-predict
      rms_list[np] = std(resid_list[np])
      chisq_list[np] = invquad(Sigma,resid_list[np])
	  sigma_param_list[np] = sigma_param
	  sigmaj_list[np] = sigmaj_opt
	  sigma_sigmaj_list[np] = sigma_sigmaj
	  frac_valid_list[np] = frac_valid_e
	  runtimes_kepler[np] += toq()
  end

   result_linear = (evidence_circ_list, P_best,sigma_P_list, param_bf_list, sigma_param_list, sigmaj_list, sigma_sigmaj_list, resid_list, rms_list, chisq_list, frac_valid_list, runtimes_circ, evidence_kepler_list, runtimes_kepler)

   #=
   Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.0));
   for npl in 1:length(P_best)
     tic()
     param_tmp = make_kepler_param_from_linear_output(result_linear,npl)
     SigmaBig = PDMat(make_Sigma(times,sigma_obs,sigmaj_list[npl]));
     optimize_phases(param_tmp,times,vels,SigmaBig)
     kepler_fit_list[npl] = laplace_approximation_kepler_model_jitter_separate(param_tmp,times,vels,Sigma0)
     evidence_kepler_list[npl] = kepler_fit_list[npl][1]/log(10)
	 println("# logp npl=",npl," (laplace): ",evidence_kepler_list[npl])
     flush(STDOUT)
     evidence_kepler_list[npl] = laplace_approximation_kepler_model_jitter_separate_gauss(param_tmp,times,vels,Sigma0)[3]/log(10)
	 println("# logp npl=",npl," (gauss): ",evidence_kepler_list[npl])
	 flush(STDOUT)
     runtimes_kepler[np] += toq()
  end

   return (evidence_circ_list, P_best,sigma_P_list, param_bf_list, sigma_param_list, sigmaj_list, sigma_sigmaj_list, resid_list, rms_list, chisq_list, frac_valid_list,evidence_kepler_list,runtimes_kepler, kepler_fit_list,runtimes_circ,runtimes_kepler)
   =#
end

function write_evidences_file(filename::String,output; use_method::Integer = 1)
   if use_method==1
      evidence = output[1]
	  runtime = convert(Array{Int64,1},ceil.(output[12]))
   elseif use_method==2
	  evidence = vcat(output[1][1],output[13])
	  runtime = convert(Array{Int64,1},ceil.(output[14]))
   end
   output_circ_data = Array{Any}(4,8)
   output_circ_data[1,1] = 0
   output_circ_data[:,1] = runtime
   output_circ_data[:,2] = collect(0:3)
   for i in 3:8
      output_circ_data[:,i] = evidence
   end
   writedlm(filename, output_circ_data, ", ")
end

function make_param_string(paramid::String, npl::Integer, plid::Integer,bf::Real, sigma::Real; runtime::Integer=0)
	   p50 = bf
	   p02 = bf-2*sigma
	   p16 = bf-sigma
	   p84 = bf+sigma
	   p98 = bf+2*sigma
	   string(runtime,',',paramid,',',npl,',',plid,',',bf,',',p50,',',p02,',',p16,',',p84,',',p98)
end

function write_parameters_file(filename::String,output; use_method::Integer = 1, input_filename::String="../data/rvs_0001.txt", num_pl::Integer=3)
    file = open(filename,"w")
	evidence = output[1]
	runtime = convert(Array{Int64,1},ceil.(output[12]))[num_pl]
    println(file,make_param_string("J",num_pl,0,output[6][num_pl],output[7][num_pl],runtime=runtime))
	(Ks,es,hs,ks) = generate_Ke_samples(output[4][num_pl],output[5][num_pl],1000)
	for plid = 1:num_pl
	   bf = output[2][plid]
	   sigma = output[3][num_pl][plid]
	   output_string = make_param_string("P",num_pl,plid,bf,sigma,runtime=runtime)
	   println(file,output_string)
	   bf = mean(Ks[:,plid])
	   sigma = std(Ks[:,plid])
	   output_string = make_param_string("K",num_pl,plid,bf,sigma,runtime=runtime)
	   println(file,output_string)
	 end
	close(file)
end

#=
for num_pl in 1:3
   write_parameters_file(string("params_",num_pl,"_000",1,".txt"),output1,input_filename="../data/rvs_0001.txt",num_pl=num_pl)
   write_parameters_file(string("params_",num_pl,"_000",2,".txt"),output2,input_filename="../data/rvs_0002.txt",num_pl=num_pl)
   write_parameters_file(string("params_",num_pl,"_000",3,".txt"),output3,input_filename="../data/rvs_0003.txt",num_pl=num_pl)
   write_parameters_file(string("params_",num_pl,"_000",4,".txt"),output4,input_filename="../data/rvs_0004.txt",num_pl=num_pl)
   write_parameters_file(string("params_",num_pl,"_000",5,".txt"),output5,input_filename="../data/rvs_0005.txt",num_pl=num_pl)
   write_parameters_file(string("params_",num_pl,"_000",6,".txt"),output6,input_filename="../data/rvs_0006.txt",num_pl=num_pl)
end
=#

function analyze_dataset(filename::String; num_Sigmas::Integer = 20, samples_per_peak::Real = 1.0, Pminsearch::Vector{Float64} = 10.0*ones(3))
  # Read dataset
  data = readdlm(filename);
  times = data[:,1];
  vels = data[:,2];
  sigma_obs = data[:,3];
  nobs = length(vels);

   (evidence_circ_list, P_best, result_laplace_circ_list, runtimes_circ) = compute_evidences_laplace_circ(times,vels,sigma_obs,num_Sigmas=num_Sigmas, Pminsearch=Pminsearch)
    write_parameters_file("tmp1.out",ans, input_filename=filename, num_pl=1)
	write_parameters_file("tmp2.out",ans, input_filename=filename, num_pl=2)
	write_parameters_file("tmp3.out",ans, input_filename=filename, num_pl=3)
end

function write_evidences_file(filename::String,output; use_method::Integer = 1, input_filename::String="../data/rvs_0001.txt")
   if use_method==1
      evidence = output[1]
	  runtime = convert(Array{Int64,1},ceil.(output[12]))
   elseif use_method==2
	  evidence = vcat(output[1][1],output[13])
	  runtime = convert(Array{Int64,1},ceil.(output[14]))
   elseif use_method==3
       evidence = zeros(4)
	   runtime = zeros(4)
	   data = readdlm(input_filename)
       times = data[:,1];
       vels = data[:,2];
       sigma_obs = data[:,3];
	   Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.0));
	   evidence[1] = NaN
       for np in 1:3
	      tic()
		  param = make_kepler_param_from_linear_output(output,np)
	      SigmaBig = PDMat(make_Sigma(times,sigma_obs,param[end]));
		  optimize_phases(param,times,vels,SigmaBig)
	      param2 = laplace_approximation_kepler_model_jitter_separate(param,times,vels,Sigma0,opt_jitter=false)[2].minimizer
		  evidence[1+np] = laplace_approximation_kepler_model_jitter_separate(param2,times,vels,Sigma0,opt_jitter=true)[1]/log(10)
		  runtime[1+np] = toc()
	   end
   elseif use_method==4
       evidence = zeros(4)
	   runtime = zeros(4)
	   data = readdlm(input_filename)
       times = data[:,1];
       vels = data[:,2];
       sigma_obs = data[:,3];
	   Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.0));
	   evidence[1] = NaN
       for np in 0:3
	      tic()
		  periods = output[2][1:np]
		  if np == 1
		     param = vcat(output[4][np][[1,2,5]],output[6])
		  elseif np == 2
		     idx = vcat(output[4][np][[1,2,5,6,9]],output[6])
		  elseif np == 3
		     idx = vcat(output[4][np][[1,2,5,6,9,10,13]],output[6])
		  end
          if np == 0
 		     evidence[1+np] = calc_map_no_planet(times,vels,sigma_obs)/log(10)
		  else
             evidence[1+np] = compute_marginalized_likelihood_near_periods_laplace(times, vels, sigma_obs, periods,  calclogprior_linear=logprior_sinusoids, make_design_matrix=make_design_matrix_circ, report_log_MAP=true)[3]/log(10) # Actually logMAP
		  end
		  runtime[1+np] = toc()
	   end
   elseif use_method==5
       evidence = zeros(4)
	   runtime = zeros(4)
	   data = readdlm(input_filename)
       times = data[:,1];
       vels = data[:,2];
       sigma_obs = data[:,3];
	   Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.0));
	   evidence[1] = NaN
       for np in 0:3
	      tic()
		  periods = output[2][1:np]
		  if np == 1
		     param = vcat(output[4][np][[1,2,5]],output[6])
		  elseif np == 2
		     idx = vcat(output[4][np][[1,2,5,6,9]],output[6])
		  elseif np == 3
		     idx = vcat(output[4][np][[1,2,5,6,9,10,13]],output[6])
		  end
          if np == 0
 		     evidence[1+np] = calc_map_no_planet(times,vels,sigma_obs)/log(10)
		  else
             evidence[1+np] = compute_marginalized_likelihood_near_periods_laplace(times, vels, sigma_obs, periods,  calclogprior_linear=logprior_epicycles, make_design_matrix=make_design_matrix_epicycle, report_log_MAP=true)[3]/log(10) # Actually logMAP
	      end
		  runtime[1+np] = toc()
	   end
   elseif use_method==6
       evidence = zeros(4)
	   runtime = zeros(4)
	   data = readdlm(input_filename)
       times = data[:,1];
       vels = data[:,2];
       sigma_obs = data[:,3];
	   Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.0));
	   evidence[1] = NaN
       for np in 0:3
	      tic()
		  if np == 0
 		     evidence[1+np] = calc_map_no_planet(times,vels,sigma_obs)/log(10)
		  else
		  param = make_kepler_param_from_linear_output(output,np)
	      SigmaBig = PDMat(make_Sigma(times,sigma_obs,param[end]));
		  optimize_phases(param,times,vels,SigmaBig)
	      #param2 = laplace_approximation_kepler_model_jitter_separate(param,times,vels,Sigma0,opt_jitter=false,report_log_MAP=true)[2].minimizer
		  #evidence[1+np] = laplace_approximation_kepler_model_jitter_separate(param2,times,vels,Sigma0,opt_jitter=true,report_log_MAP=true)[1]/log(10)
		  evidence[1+np] = laplace_approximation_kepler_model_jitter_separate(param,times,vels,Sigma0,opt_jitter=false,report_log_MAP=true)[1]/log(10)
		  end
		  runtime[1+np] = toc()
	   end
   end
   output_circ_data = Array{Any}(4,8)
   output_circ_data[1,1] = 0
   output_circ_data[:,1] = runtime
   output_circ_data[:,2] = collect(0:3)
   for i in 3:8
      output_circ_data[:,i] = evidence
   end
   writedlm(filename, output_circ_data, ", ")
end

#=
include("calc_evidence_laplace_lowe.jl")
output1 = analyze_dataset("../data/rvs_0001.txt",num_Sigmas = 4,samples_per_peak=1,Pminsearch=[30.0,10.0,10.])
write_evidences_file("laplace_linearized_circ/evidences_0001.txt",output1)
write_evidences_file("laplace_linearized_ecc/evidences_0001.txt",output1,use_method=2)
@save "output_kepler_1.jld"

output1 = analyze_dataset("../data/rvs_0001.txt",num_Sigmas = 40,samples_per_peak=4,Pminsearch=[30.0,10.0,10.0])

output1 = analyze_dataset("../data/rvs_0001.txt",num_Sigmas = 4,samples_per_peak=1,Pminsearch=[10.0,10.0,1.2])
output2 = analyze_dataset("../data/rvs_0002.txt",num_Sigmas = 4,samples_per_peak=1,Pminsearch=[10.0,10.0,1.2])
output3 = analyze_dataset("../data/rvs_0003.txt",num_Sigmas = 4,samples_per_peak=1,Pminsearch=[10.0,10.0,1.2])
output4 = analyze_dataset("../data/rvs_0004.txt",num_Sigmas = 4,samples_per_peak=1,Pminsearch=[10.0,10.0,1.2])
output5 = analyze_dataset("../data/rvs_0005.txt",num_Sigmas = 4,samples_per_peak=1,Pminsearch=[10.0,10.0,1.2])
output6 = analyze_dataset("../data/rvs_0006.txt",num_Sigmas = 4,samples_per_peak=1,Pminsearch=[10.0,10.0,1.2])
@save "output_kepler_4.jld"
=#
#=
include("calc_evidence_laplace_lowe.jl")
output1 = analyze_dataset("../data/rvs_0001.txt",num_Sigmas = 40,samples_per_peak=4,Pminsearch=[10.0,10.,1.2])
write_evidences_file("laplace_linearized_circ/evidences_0001.txt",output1)
write_evidences_file("laplace_linearized_ecc/evidences_0001.txt",output1,use_method=2)
@save "output_kepler_40_4_10.jld"
output2 = analyze_dataset("../data/rvs_0002.txt",num_Sigmas = 40,samples_per_peak=4,Pminsearch=[10.0,10.,1.2])
write_evidences_file("laplace_linearized_circ/evidences_0002.txt",output2)
write_evidences_file("laplace_linearized_ecc/evidences_0002.txt",output2,use_method=2)
output3 = analyze_dataset("../data/rvs_0003.txt",num_Sigmas = 40,samples_per_peak=4,Pminsearch=[10.0,10.,1.2])
write_evidences_file("laplace_linearized_circ/evidences_0003.txt",output3)
write_evidences_file("laplace_linearized_ecc/evidences_0003.txt",output3,use_method=2)
output4 = analyze_dataset("../data/rvs_0004.txt",num_Sigmas = 40,samples_per_peak=4,Pminsearch=[10.0,10.,1.2])
write_evidences_file("laplace_linearized_circ/evidences_0004.txt",output4)
write_evidences_file("laplace_linearized_ecc/evidences_0004.txt",output4,use_method=2)
output5 = analyze_dataset("../data/rvs_0005.txt",num_Sigmas = 40,samples_per_peak=4,Pminsearch=[10.0,10.,1.2])
write_evidences_file("laplace_linearized_circ/evidences_0005.txt",output5)
write_evidences_file("laplace_linearized_ecc/evidences_0005.txt",output5,use_method=2)
output6 = analyze_dataset("../data/rvs_0006.txt",num_Sigmas = 40,samples_per_peak=4,Pminsearch=[10.0,10.,1.2])
write_evidences_file("laplace_linearized_circ/evidences_0006.txt",output6)
write_evidences_file("laplace_linearized_ecc/evidences_0006.txt",output6,use_method=2)
@save "output_kepler_40_4_10.jld"
=#

#=
data = readdlm("../data/rvs_0001.txt")
times = data[:,1];
vels = data[:,2];
sigma_obs = data[:,3];
output1 = analyze_dataset("../data/rvs_0001.txt",num_Sigmas = 40,Pminsearch=[30.0,10.0,10.0])
param_tmp = make_kepler_param_from_linear_output(output1,1)
SigmaBig = PDMat(make_Sigma(times,sigma_obs,param_tmp[end]));
optimize_phases(param_tmp,times,vels,SigmaBig)
res1 = laplace_approximation_kepler_model_jitter_separate(param_tmp,times,vels,Sigma0)
res1[1]/log(10)

param_tmp = make_kepler_param_from_linear_output(output1,2)
SigmaBig = PDMat(make_Sigma(times,sigma_obs,param_tmp[end]));
optimize_phases(param_tmp,times,vels,SigmaBig)
res2 = laplace_approximation_kepler_model_jitter_separate(param_tmp,times,vels,Sigma0)
res2[1]/log(10)

param_tmp = make_kepler_param_from_linear_output(output1,3)
SigmaBig = PDMat(make_Sigma(times,sigma_obs,param_tmp[end]));
optimize_phases(param_tmp,times,vels,SigmaBig)
res3 = laplace_approximation_kepler_model_jitter_separate(param_tmp,times,vels,Sigma0)
res3[1]/log(10)

=#

#=
param1a = [42.0,5.0,0.01,0.01,1.,0.]
res = laplace_approximation_kepler_model(param1a,times,vels,SigmaBig,loglikelihood_fixed_jitter)
param1a = res[2].minimizer;
laplace_approximation_kepler_model_jitter_separate(param1a,times,vels,Sigma0)[1]/log(10)
=#

function compute_marginalized_likelihood_cuba(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}, min_param::Vector{T}, max_param::Vector{T}; num_Sigmas::Integer = 20 ) where T<:Real
  @assert length(min_param)==length(max_param)
  log_normalization = 0.0
  Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.))
  Sigma_cache = make_Sigma_cache(num_Sigmas,Sigma0)
  function log_integrand_cuba(x::Vector)
    param = (max_param.-min_param).*x.+min_param

	logp = loglikelihood_marinalize_jitter(param,times,vels,Sigma_cache=Sigma_cache)
	#println("# param= ",param," logp= ",logp)
	logp
  end
  function integrand_cuba(x::Vector, output::Vector)
    output[1] = exp(log_integrand_cuba(x)-log_normalization)
  end

  log_normalization = log_integrand_cuba(0.5*ones(min_param))
  result = suave(integrand_cuba,length(min_param),1,reltol=1e-3,maxevals=10000)
  #result = log(result.integral[1]) + log(prod(max_param.-min_param))
  #result+log_normalization
end

function compute_marginalized_likelihood_cuba_near(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}, param::Vector{T}; num_Sigmas::Integer = 20 ) where T<:Real
  const nparam_per_pl = 5
  const nparam_non_pl = 2
  npl = convert(Int64,floor((length(param)-1)//nparam_per_pl))
  param_min = similar(param); param_max = similar(param);
  for p in 1:npl
    param_min[1+(p-1)*nparam_per_pl] = param[1+(p-1)*nparam_per_pl]*0.98
    param_max[1+(p-1)*nparam_per_pl] = param[1+(p-1)*nparam_per_pl]*1.02
	param_min[2+(p-1)*nparam_per_pl] = param[2+(p-1)*nparam_per_pl]*0.5
    param_max[2+(p-1)*nparam_per_pl] = param[2+(p-1)*nparam_per_pl]*2.00
	param_min[3+(p-1)*nparam_per_pl] = -1.0
    param_max[3+(p-1)*nparam_per_pl] = 1.0
	param_min[4+(p-1)*nparam_per_pl] = -1.0
    param_max[4+(p-1)*nparam_per_pl] = 1.0
	param_min[5+(p-1)*nparam_per_pl] = 0.0
    param_max[5+(p-1)*nparam_per_pl] = 2pi
  end
  param_min[1+npl*nparam_per_pl] = param[1+npl*nparam_per_pl]-2.0
  param_max[1+npl*nparam_per_pl] = param[1+npl*nparam_per_pl]+2.0
  param_min[1+npl*nparam_per_pl] = 0.0
  param_max[1+npl*nparam_per_pl] = 2*param[2+npl*nparam_per_pl]
  res = compute_marginalized_likelihood_cuba(times,vels,sigma_obs,param_min,param_max,num_Sigmas=num_Sigmas)
end

#param_tmp = make_kepler_param_from_linear_output(output1,1)
#compute_marginalized_likelihood_cuba_near(times,vels,sigma_obs,param_tmp,num_Sigmas=3)

#=
include("calc_evidence_laplace_lowe.jl")
using JLD
@load "output_kepler_40_4_10.jld"
@load "output_case_4.jld"

write_evidences_file("map_linearized_circ/evidences_0001.txt",output1,use_method=4,input_filename="../data/rvs_0001.txt")
write_evidences_file("map_linearized_circ/evidences_0002.txt",output2,use_method=4,input_filename="../data/rvs_0002.txt")
write_evidences_file("map_linearized_circ/evidences_0003.txt",output3,use_method=4,input_filename="../data/rvs_0003.txt")
write_evidences_file("map_linearized_circ/evidences_0004.txt",output4,use_method=4,input_filename="../data/rvs_0004.txt")
write_evidences_file("map_linearized_circ/evidences_0005.txt",output5,use_method=4,input_filename="../data/rvs_0005.txt")
write_evidences_file("map_linearized_circ/evidences_0006.txt",output6,use_method=4,input_filename="../data/rvs_0006.txt")
@save "output_method4.jld"

write_evidences_file("map_linearized_ecc/evidences_0001.txt",output1,use_method=5,input_filename="../data/rvs_0001.txt")
write_evidences_file("map_linearized_ecc/evidences_0002.txt",output2,use_method=5,input_filename="../data/rvs_0002.txt")
write_evidences_file("map_linearized_ecc/evidences_0003.txt",output3,use_method=5,input_filename="../data/rvs_0003.txt")
write_evidences_file("map_linearized_ecc/evidences_0004.txt",output4,use_method=5,input_filename="../data/rvs_0004.txt")
write_evidences_file("map_linearized_ecc/evidences_0005.txt",output5,use_method=5,input_filename="../data/rvs_0005.txt")
write_evidences_file("map_linearized_ecc/evidences_0006.txt",output6,use_method=5,input_filename="../data/rvs_0006.txt")
@save "output_method5.jld"

write_evidences_file("map_keplerian/evidences_0001.txt",output1,use_method=6,input_filename="../data/rvs_0001.txt")
write_evidences_file("map_keplerian/evidences_0002.txt",output2,use_method=6,input_filename="../data/rvs_0002.txt")
write_evidences_file("map_keplerian/evidences_0003.txt",output3,use_method=6,input_filename="../data/rvs_0003.txt")
write_evidences_file("map_keplerian/evidences_0004.txt",output4,use_method=6,input_filename="../data/rvs_0004.txt")
write_evidences_file("map_keplerian/evidences_0005.txt",output5,use_method=6,input_filename="../data/rvs_0005.txt")
write_evidences_file("map_keplerian/evidences_0006.txt",output6,use_method=6,input_filename="../data/rvs_0006.txt")
@save "output_method6.jld"

write_evidences_file("laplace_keplerian/evidences_0001.txt",output1,use_method=3,input_filename="../data/rvs_0001.txt")
write_evidences_file("laplace_keplerian/evidences_0002.txt",output2,use_method=3,input_filename="../data/rvs_0002.txt")
write_evidences_file("laplace_keplerian/evidences_0003.txt",output3,use_method=3,input_filename="../data/rvs_0003.txt")
write_evidences_file("laplace_keplerian/evidences_0004.txt",output4,use_method=3,input_filename="../data/rvs_0004.txt")
write_evidences_file("laplace_keplerian/evidences_0005.txt",output5,use_method=3,input_filename="../data/rvs_0005.txt")
write_evidences_file("laplace_keplerian/evidences_0006.txt",output6,use_method=3,input_filename="../data/rvs_0006.txt")
@save "output_method3.jld"

=#

function direct_integration(data::Vector{T}, model::Array{T,2}, Sigma::SigmaT; FIM::PDMat{T} = make_fischer_information_matrix(model,Sigma), calclogprior::Function = x->(0.0,0.0), report_log_MAP::Bool = false ) 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]
  if report_log_MAP     return (param_bf_linalg, sigma_param,log_mode)  end
  LaplaceApprox = 0.5*nparam*log(2pi)-0.5*logdet(FIM)+log_mode + log_jacobian
  return (param_bf_linalg,sigma_param,LaplaceApprox)
end

function compute_marginalized_likelihood_const_gaussgk(times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}; Sigma0 = PDMat(make_Sigma(times,sigma_obs,0.0)) ) where T<:Real
    X = ones(length(times),1)
    const nparam = size(X,2)
    function log_integrand(x)
        Sigma = Sigma0 + ScalMat(length(times),x)
		compute_marginalized_likelihood_fix_periods_sigmaj(times, vels, sigma_obs, Float64[], x, Sigma=Sigma,X=X)+logprior_jitter(x)
	end
    int_norm = log_integrand(0.0)
	function integrand(x)
	   exp(log_integrand(x)-int_norm)
	end
	integral = QuadGK.quadgk(integrand,0.0,2.0,99.0,reltol=1e-14)

    (log(integral[1]) + int_norm, log(integral[2]/integral[1]) )
end

function loglikelihood_zero_pl( param::Vector{T}, times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}; calclogprior::Function = x->(0.0,0.0) ) where T<:Real
  @assert length(vels) == length(sigma_obs) # size(Sigma,1)
  nobs = length(vels)

  predict = ones(nobs)*param[1]
  delta = vels.-predict
  Sigma = make_Sigma(times,sigma_obs,param[2])
  chisq = invquad(Sigma,delta)
  loglikelihood_mode = 0.5*(-chisq-logdet(Sigma)-nobs*log(2pi))
end

function logposterior_zero_pl( param::Vector{T}, times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}) where T<:Real
  logprior_offset(param[1])+logprior_jitter(param[2])+loglikelihood_zero_pl(param,times,vels,sigma_obs)
end

function loglikelihood_zero_pl( param::Vector{T}, data::Vector{T}, Sigma::SigmaT; calclogprior::Function = x->(0.0,0.0) ) where {T<:Real, SigmaT<:AbstractPDMat{T} }
  @assert length(data) == size(Sigma,1)
  nobs = length(data)

  predict = ones(nobs)*param[1]
  delta = data.-predict
  Sigma = make_Sigma(times,sigma_obs,param[2])
  chisq = invquad(Sigma,delta)

  loglikelihood_mode = 0.5*(-chisq-logdet(Sigma)-nobs*log(2pi))
end

function logposterior_zero_pl( param::Vector{T}, data::Vector{T}, Sigma::SigmaT ) where {T<:Real, SigmaT<:AbstractPDMat{T} }
  logprior_offset(param[1])+logprior_jitter(param[2])+loglikelihood_zero_pl(param,data,Sigma)
end

function cuba_integrand(x::Vector{T}, output::Vector{T}) where T<:Real
  param_cuba = similar(x)
  param_cuba[1] = -1000+2000*x[1]
  param_cuba[2] = 0+99*x[2]
  output[1] = exp(logposterior_zero_pl(param_cuba,times,vels,sigma_obs)-int_norm)

end

#=
int_norm = logposterior_zero_pl([0.0,2.0],times,vels,sigma_obs)
diovenne(cuba_integrand,2)
=#

function logposterior_zero_pl_marginalize_offset( sigmaj::T, times::Vector{T}, vels::Vector{T}, sigma_obs::Vector{T}; xmax::T=Cmax, xmin::T=-xmax) where T<:Real
  #@assert length(vels) == length(times) = length(sigma_obs)
  param = zeros(2)

  function log_integrand(x::T)
    param[1] = x
    logprior_offset(param[1])+logprior_jitter(param[2])+loglikelihood_zero_pl(param,times,vels,sigma_obs)
  end

  param[2] = sigmaj
  int_norm = log_integrand(0.)
  function integrand(x::T)
     exp(log_integrand(x)-int_norm)
  end

  log(QuadGK.quadgk(integrand,xmin,xmax)[1])+int_norm
end

#=
data = readdlm("../data/rvs_0001.txt");
times = data[:,1];
vels = data[:,2];
sigma_obs = data[:,3];
Sigma0 = make_Sigma(times,sigma_obs,0.0)
Sigma2 = make_Sigma(times,sigma_obs,2.0)
=#

#=
include("calc_evidence_laplace_lowe.jl")
zero_pl_tests =
	[("rvs_0001.txt", [-0.728509750583953064E+00, 0.208717785347551521E+01]),
	("rvs_0002.txt", [ 0.653665406348636679E+01, 0.171886843613598783E+01]),
	("rvs_0003.txt", [-0.869753669161968901E+01, 0.994907921829712771E+00]),
	("rvs_0004.txt", [-0.658082070168154587E+01, 0.682875011391089748E+00]),
	("rvs_0005.txt", [-0.486556982931438142E+01, 0.930406316677828249E+00]),
	("rvs_0006.txt", [-0.108329221460122653E+02, 0.115487852990739137E+01])]

filename_out = "tmp.out"
file_out = open(filename_out,"w")
for test in zero_pl_tests
  filename = joinpath("..","data",test[1])
  data = readdlm(filename)
  times = data[:,1];
  vels = data[:,2];
  sigma_obs = data[:,3];
  lnlikelihood = loglikelihood_zero_pl(test[2],times,vels,sigma_obs)
  println(file_out,test, ":  ", lnlikelihood)
  test2 = copy(test[2])
  test2[1] = 0.0
  lnlikelihood = loglikelihood_zero_pl(test2,times,vels,sigma_obs)
  println(file_out,(test[1],test2), ":  ", lnlikelihood)
end
close(file_out)

=#

#loglikelihood_zero_pl([-0.639837062706153,0.0],times,vels,sigma_obs)


text:
using Revise
using FwiFlow
using AdFem
using ADCME
using PyPlot
using PyCall
using Statistics
using DelimitedFiles
using ADCMEKit
np = pyimport("numpy")

function visualize_obs(o1)
    close("all")
    plot((0:NT)*Δt, o1[:, 1:m+1])
    xlabel("t")
    ylabel("displacement")
end

function layer_model(m, n)
    z = zeros(4*n*m)
    k = 0
    mm = ones(n)
    mm[n÷3:end] .= 1.5
    for i = 1:n
        for j = 1:m
            for p = 1:2
                for q = 1:2
                    k += 1
                    z[k] = mm[i]
                end
            end
        end
    end
    return z
end

function sin_model(m, n)
end

function visualize_invη(o)
    o = o * 1e12
    Z = zeros(n, m)
    k = 0
    for i = 1:n
        for j = 1:m
            Z[i, j] = mean(o[k+1:k+4])
            k += 4
        end
    end
    pcolormesh(Z)
    axis("scaled")
    xlabel("x")
    ylabel("y")
    gca().invert_yaxis()
    colorbar()
end

ext: jl

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
if basename(pwd()) == "aoc"
    cd("2019/2")
end
include("../IntCodeMachine/src/IntCodeMachine.jl")
using .IntCodeMachine

function part1()
    c = IntCode(readchomp("input.txt"))
    c[1] = 12
    c[2] = 2
    print(execute!(c)[0])
end

function part2()
    c = IntCode(readchomp("input.txt"))
    initmem = copy(c.memory)
    for x in 0:99, y in 0:99
        reset!(c, initmem)
        c[1] = x
        c[2] = y
        if execute!(c)[0] == 19690720
            return c[1] * 100 + c[2]
        end
    end
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
module LibXXX

import Libdl

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("LibXXX was not build properly. Please run Pkg.build(\"LibXXX\").")
end
include(depsjl_path)
# Module initialization function
function __init__()
    check_deps()
end

using CEnum

include("ctypes.jl")
export Ctm, Ctime_t, Cclock_t

include(joinpath(@__DIR__, "..", "gen", "libXXX_common.jl"))
include(joinpath(@__DIR__, "..", "gen", "libXXX_api.jl"))

# export everything
#foreach(names(@__MODULE__, all=true)) do s
#    if startswith(string(s), "SOME_PREFIX")
#        @eval export $s
#    end
#end

end # module


text:
function main()
    N = parse(Int, readline())
    ps = parse.(Int, split(readline()))

    dp = falses(sum(ps) + 1)
    dp[begin] = true # Score zero
    for p in ps
        dp = dp .| vcat(falses(p), dp[begin:end-p])
    end

    println(sum(dp))
end

main()


text:
OTLP_VERSION = v"0.13.0"

using Downloads
using CodecZlib
using Tar
using gRPCClient

download_url = "https://github.com/open-telemetry/opentelemetry-proto/archive/refs/tags/v$OTLP_VERSION.tar.gz"

open(Downloads.download(download_url)) do tar_gz
    tar = GzipDecompressorStream(tar_gz)
    dir = Tar.extract(tar)
    proto_dir = joinpath(dir, "opentelemetry-proto-$OTLP_VERSION")

    trace_proto_file = joinpath(
        proto_dir,
        "opentelemetry",
        "proto",
        "collector",
        "trace",
        "v1",
        "trace_service.proto",
    )

    trace_dir = joinpath(@__DIR__, "..", "src", "trace_generated")

    gRPCClient.generate(trace_proto_file; outdir = trace_dir, includes = [proto_dir])

    metrics_proto_file = joinpath(
        proto_dir,
        "opentelemetry",
        "proto",
        "collector",
        "metrics",
        "v1",
        "metrics_service.proto",
    )

    metrics_dir = joinpath(@__DIR__, "..", "src", "metrics_generated")

    gRPCClient.generate(metrics_proto_file; outdir = metrics_dir, includes = [proto_dir])

    logs_proto_file = joinpath(
        proto_dir,
        "opentelemetry",
        "proto",
        "collector",
        "logs",
        "v1",
        "logs_service.proto",
    )

    logs_dir = joinpath(@__DIR__, "..", "src", "logs_generated")

    gRPCClient.generate(logs_proto_file; outdir = logs_dir, includes = [proto_dir])
end


text:
module TestCauerLowPassFilter

println("\nTestCauerLowPassFilter: Demonstrating the ability to simulate an electrical model translated from Modelica Standar Library")

using Modia
using Modia.Electric

# Desired:
#   using ModiaMath: plot
#
# In order that ModiaMath need not to be defined in the user environment, it is included via Modia:
using Modia.ModiaMath: plot

@model CauerLowPassOPV begin
  # Cauer low pass filter with operational amplifiers
  l1=1.304
  l2=0.8586
  c1=1.072
  c2=1/(1.704992^2*l1)
  c3=1.682
  c4=1/(1.179945^2*l2)
  c5=0.7262

  C1=Capacitor(C=c1 + c2)
  C2=Capacitor(C=c2)
  C3=Capacitor(C=l1)
  C4=Capacitor(C=c4)
  C5=Capacitor(C=c2, v=Float(state=false))
  R1=Resistor(R=1)
  R2=Resistor(R=1)
  R3=Resistor(R=1)
  Op1=IdealOpAmp3Pin()
  G=Ground()
  R4=Resistor(R=-1)
  R5=Resistor(R=-1)
  Op2=IdealOpAmp3Pin()
  Op3=IdealOpAmp3Pin()
  G1=Ground()
  R6=Resistor(R=1)
  R7=Resistor(R=1)
  C6=Capacitor(C=c2 + c3 + c4, v=Float(state=false))
  R8=Resistor(R=-1)
  R9=Resistor(R=-1)
  R10=Resistor(R=1)
  Op4=IdealOpAmp3Pin()
  Op5=IdealOpAmp3Pin()
  C7=Capacitor(C=l2, v=Float(start=0))
  C8=Capacitor(C=c4, v=Float(state=false))
  C9=Capacitor(C=c4 + c5, v=Float(state=false))
  R11=Resistor(R=1)
  n1=Pin()
  n2=Pin()
  n3=Pin()
  n4=Pin()
  n5=Pin()
  p1=Pin()
  n6=Pin()
  n7=Pin()
  n8=Pin()
  p2=Pin()
  out1=Pin()
  p3=Pin()
  n9=Pin()
  n10=Pin()
  n11=Pin()
  G2=Ground()
  n12=Pin()
  n13=Pin()
  p4=Pin()
  n14=Pin()
  G3=Ground()
  G4=Ground()
#  V=StepVoltage(startTime=10, V=1)
  V=ConstantVoltage(V=1)
  Ground1=Ground()
  obs=Float(start=0.0)
@equations begin
  connect(Op1.in_p, G.p)
  connect(G1.p, Op2.in_p)
  connect(R1.n, n1)
  connect(n1, Op1.in_n)
  connect(C2.n, n1)
  connect(R2.n, n2)
  connect(n1, n2)
  connect(n2, n3)
  connect(n3, C1.p)
  connect(n3, R3.p)
  connect(C1.n, n4)
  connect(R3.n, n4)
  connect(n4, Op1.out)
  connect(R4.p, Op1.out)
  connect(C5.p, Op1.out)
  connect(R4.n, n5)
  connect(n5, Op2.in_n)
  connect(C3.p, n5)
  connect(R5.n, n5)
  connect(R5.p, p1)
  connect(C2.p, p1)
  connect(C3.n, n6)
  connect(n6, Op2.out)
  connect(R2.p, n6)
  connect(Op2.out, R7.p)
  connect(R7.n, n7)
  connect(n7, Op3.in_n)
  connect(C5.n, n7)
  connect(R6.n, n8)
  connect(n7, n8)
  connect(C6.p, p2)
  connect(n8, p2)
  connect(C4.n, p2)
  connect(C6.n, Op3.out)
  connect(R9.p, Op3.out)
  connect(Op3.out, out1)
  connect(p1, out1)
  connect(out1, C8.p)
  connect(C4.p, p3)
  connect(p3, R8.p)
  connect(R8.n, n9)
  connect(n9, n10)
  connect(R9.n, n10)
  connect(n10, Op4.in_n)
  connect(n9, C7.p)
  connect(C7.n, n11)
  connect(R6.p, n11)
  connect(n11, Op4.out)
  connect(Op4.out, R10.p)
  connect(G2.p, Op3.in_p)
  connect(R11.n, n12)
  connect(p3, n12)
  connect(C9.n, n13)
  connect(n12, n13)
  connect(n13, Op5.out)
  connect(C9.p, p4)
  connect(R11.p, p4)
  connect(R10.n, n14)
  connect(p4, n14)
  connect(Op5.in_n, n14)
  connect(C8.n, n14)
  connect(Op4.in_p, G3.p)
  connect(Op5.in_p, G4.p)
  connect(V.p, Ground1.p)
  connect(V.n, R1.p)
  der(obs) = 100000*(Op5.out.v-obs)
  end
end

result = checkSimulation(CauerLowPassOPV, 60, "C9.v", -0.5003139583081188, logTranslation=false, removeSingularities=false, logTiming=true)

@model CauerLowPassOPVWithoutNodes begin
  # Cauer low pass filter with operational amplifiers
  l1=1.304
  l2=0.8586
  c1=1.072
  c2=1/(1.704992^2*l1)
  c3=1.682
  c4=1/(1.179945^2*l2)
  c5=0.7262

  C1=Capacitor(C=c1 + c2)
  C2=Capacitor(C=c2)
  C3=Capacitor(C=l1)
  C4=Capacitor(C=c4)
  C5=Capacitor(C=c2, v=Float(state=false))
  R1=Resistor(R=1)
  R2=Resistor(R=1)
  R3=Resistor(R=1)
  Op1=IdealOpAmp3Pin()
  G=Ground()
  R4=Resistor(R=-1)
  R5=Resistor(R=-1)
  Op2=IdealOpAmp3Pin()
  Op3=IdealOpAmp3Pin()
  G1=Ground()
  R6=Resistor(R=1)
  R7=Resistor(R=1)
  C6=Capacitor(C=c2 + c3 + c4, v=Float(state=false))
  R8=Resistor(R=-1)
  R9=Resistor(R=-1)
  R10=Resistor(R=1)
  Op4=IdealOpAmp3Pin()
  Op5=IdealOpAmp3Pin()
  C7=Capacitor(C=l2, v=Float(start=0))
  C8=Capacitor(C=c4, v=Float(state=false))
  C9=Capacitor(C=c4 + c5, v=Float(state=false))
  R11=Resistor(R=1)

  G2=Ground()
  G3=Ground()
  G4=Ground()
#  V=StepVoltage(startTime=1, V=1)
  V=ConstantVoltage(V=1)
  Ground1=Ground()
@equations begin
  connect(Op1.in_p, G.p)
  connect(G1.p, Op2.in_p)
  connect(R1.n, Op1.in_n)
  connect(C2.n, R1.n)
  connect(R1.n, R2.n)
  connect(R2.n, C1.p)
  connect(R2.n, R3.p)
  connect(R3.n, C1.n)
  connect(C1.n, Op1.out)
  connect(R4.p, Op1.out)
  connect(C5.p, Op1.out)
  connect(R4.n, Op2.in_n)
  connect(C3.p, R4.n)
  connect(R5.n, R4.n)
  connect(C2.p, R5.p)
  connect(C3.n, Op2.out)
  connect(R2.p, C3.n)
  connect(Op2.out, R7.p)
  connect(R7.n, Op3.in_n)
  connect(C5.n, R7.n)
  connect(R7.n, R6.n)
  connect(R6.n, C6.p)
  connect(C4.n, C6.p)
  connect(C6.n, Op3.out)
  connect(R9.p, Op3.out)
  connect(R5.p, Op3.out)
  connect(Op3.out, C8.p)
  connect(C4.p, R8.p)
  connect(R9.n, R8.n)
  connect(R8.n, Op4.in_n)
  connect(R8.n, C7.p)
  connect(R6.p, C7.n)
  connect(C7.n, Op4.out)
  connect(Op4.out, R10.p)
  connect(G2.p, Op3.in_p)
  connect(C4.p, R11.n)
  connect(R11.n, C9.n)
  connect(C9.n, Op5.out)
  connect(R11.p, C9.p)
  connect(C9.p, R10.n)
  connect(Op5.in_n, R10.n)
  connect(C8.n, R10.n)
  connect(Op4.in_p, G3.p)
  connect(Op5.in_p, G4.p)
  connect(V.p, Ground1.p)
  connect(V.n, R1.p)
  end
end
#=
result = checkSimulation(CauerLowPassOPVWithoutNodes, 60, "C9.v", -0.5003139583081188)
plot(result, ("C9.v"), figure=20)
result = checkSimulation(CauerLowPassOPVWithoutNodes, 60, "C9.v", -0.5003139583081188, tearing=true)
plot(result, ("C9.v"), figure=20)
=#

end


text:
using .Dagger

HAS_DAGGER[] = true


text:
using LazyArrays, LinearAlgebra
import LazyArrays: @lazymul, @lazylmul, @lazyldiv, materialize!, MemoryLayout, triangulardata, LazyLayout, LazyArrayApplyStyle, UnknownLayout

# used to test general matrix backends
struct MyMatrix{T} <: AbstractMatrix{T}
    A::Matrix{T}
end

MyMatrix{T}(::UndefInitializer, n::Int, m::Int) where T = MyMatrix{T}(Array{T}(undef, n, m))
MyMatrix(A::AbstractMatrix{T}) where T = MyMatrix{T}(Matrix{T}(A))
Base.convert(::Type{MyMatrix{T}}, A::MyMatrix{T}) where T = A
Base.convert(::Type{MyMatrix{T}}, A::MyMatrix) where T = MyMatrix(convert(AbstractArray{T}, A.A))
Base.convert(::Type{MyMatrix}, A::MyMatrix)= A
Base.convert(::Type{AbstractArray{T}}, A::MyMatrix) where T = MyMatrix(convert(AbstractArray{T}, A.A))
Base.convert(::Type{AbstractMatrix{T}}, A::MyMatrix) where T = MyMatrix(convert(AbstractArray{T}, A.A))
Base.convert(::Type{MyMatrix{T}}, A::AbstractArray{T}) where T = MyMatrix{T}(A)
Base.convert(::Type{MyMatrix{T}}, A::AbstractArray) where T = MyMatrix{T}(convert(AbstractArray{T}, A))
Base.convert(::Type{MyMatrix}, A::AbstractArray{T}) where T = MyMatrix{T}(A)
Base.getindex(A::MyMatrix, kj...) = A.A[kj...]
Base.getindex(A::MyMatrix, ::Colon, j::AbstractVector) = MyMatrix(A.A[:,j])
Base.setindex!(A::MyMatrix, v, kj...) = setindex!(A.A, v, kj...)
Base.size(A::MyMatrix) = size(A.A)
Base.similar(::Type{MyMatrix{T}}, m::Int, n::Int) where T = MyMatrix{T}(undef, m, n)
Base.similar(::MyMatrix{T}, m::Int, n::Int) where T = MyMatrix{T}(undef, m, n)
Base.similar(::MyMatrix, ::Type{T}, m::Int, n::Int) where T = MyMatrix{T}(undef, m, n)
LinearAlgebra.factorize(A::MyMatrix) = factorize(A.A)

@lazymul MyMatrix
@lazyldiv MyMatrix

struct MyUpperTriangular{T} <: AbstractMatrix{T}
    A::UpperTriangular{T,Matrix{T}}
end

MyUpperTriangular{T}(::UndefInitializer, n::Int, m::Int) where T = MyUpperTriangular{T}(UpperTriangular(Array{T}(undef, n, m)))
MyUpperTriangular(A::AbstractMatrix{T}) where T = MyUpperTriangular{T}(UpperTriangular(Matrix{T}(A)))
Base.convert(::Type{MyUpperTriangular{T}}, A::MyUpperTriangular{T}) where T = A
Base.convert(::Type{MyUpperTriangular{T}}, A::MyUpperTriangular) where T = MyUpperTriangular(convert(AbstractArray{T}, A.A))
Base.convert(::Type{MyUpperTriangular}, A::MyUpperTriangular)= A
Base.convert(::Type{AbstractArray{T}}, A::MyUpperTriangular) where T = MyUpperTriangular(convert(AbstractArray{T}, A.A))
Base.convert(::Type{AbstractMatrix{T}}, A::MyUpperTriangular) where T = MyUpperTriangular(convert(AbstractArray{T}, A.A))
Base.convert(::Type{MyUpperTriangular{T}}, A::AbstractArray{T}) where T = MyUpperTriangular{T}(A)
Base.convert(::Type{MyUpperTriangular{T}}, A::AbstractArray) where T = MyUpperTriangular{T}(convert(AbstractArray{T}, A))
Base.convert(::Type{MyUpperTriangular}, A::AbstractArray{T}) where T = MyUpperTriangular{T}(A)
Base.getindex(A::MyUpperTriangular, kj...) = A.A[kj...]
Base.getindex(A::MyUpperTriangular, ::Colon, j::AbstractVector) = MyUpperTriangular(A.A[:,j])
Base.setindex!(A::MyUpperTriangular, v, kj...) = setindex!(A.A, v, kj...)
Base.size(A::MyUpperTriangular) = size(A.A)
Base.similar(::Type{MyUpperTriangular{T}}, m::Int, n::Int) where T = MyUpperTriangular{T}(undef, m, n)
Base.similar(::MyUpperTriangular{T}, m::Int, n::Int) where T = MyUpperTriangular{T}(undef, m, n)
Base.similar(::MyUpperTriangular, ::Type{T}, m::Int, n::Int) where T = MyUpperTriangular{T}(undef, m, n)
LinearAlgebra.factorize(A::MyUpperTriangular) = factorize(A.A)

MemoryLayout(::Type{MyUpperTriangular{T}}) where T = MemoryLayout(UpperTriangular{T,Matrix{T}})
triangulardata(A::MyUpperTriangular) = triangulardata(A.A)

@lazylmul MyUpperTriangular

struct MyLazyArray{T,N} <: AbstractArray{T,N}
    data::Array{T,N}
end

Base.size(A::MyLazyArray) = size(A.data)
Base.getindex(A::MyLazyArray, j::Int...) = A.data[j...]
LazyArrays.MemoryLayout(::Type{<:MyLazyArray}) = LazyLayout()
LinearAlgebra.factorize(A::MyLazyArray) = factorize(A.data)

@testset "lazymul/ldiv tests" begin
    @testset "*/" begin
        A = randn(5,5)
        B = randn(5,5)
        x = randn(5)
        @test MyMatrix(A)*x ≈ apply(*,MyMatrix(A),x) ≈ A*x
        @test all(MyMatrix(A)*MyMatrix(A) .=== apply(*,MyMatrix(A),MyMatrix(A)))
        @test all(MyMatrix(A)*A .=== apply(*,MyMatrix(A),A))
        @test all(A*MyMatrix(A) .=== apply(*,A,MyMatrix(A)))
        @test MyMatrix(A)*MyMatrix(A) ≈ MyMatrix(A)*A ≈ A*MyMatrix(A) ≈ A^2

        @test MyMatrix(A)*MyMatrix(A)*MyMatrix(A) ≈ apply(*,MyMatrix(A),MyMatrix(A),MyMatrix(A)) ≈ A^3

        @test all(UpperTriangular(A) * MyMatrix(A) .=== apply(*,UpperTriangular(A), MyMatrix(A)))
        @test all(MyMatrix(A) * UpperTriangular(A) .=== apply(*, MyMatrix(A),UpperTriangular(A)))
        @test all(Diagonal(A) * MyMatrix(A) .=== apply(*,Diagonal(A), MyMatrix(A)))
        @test all(MyMatrix(A) * Diagonal(A) .=== apply(*, MyMatrix(A),Diagonal(A)))

        @test all(MyMatrix(A)' * x .=== apply(*,MyMatrix(A)',x))

        @test all(MyMatrix(A)' * MyMatrix(A)' .=== apply(*,MyMatrix(A)', MyMatrix(A)'))
        @test all(MyMatrix(A)' * A' .=== apply(*,MyMatrix(A)', A'))
        @test all(A' * MyMatrix(A)' .=== apply(*,MyMatrix(A)', MyMatrix(A)'))
        @test all(MyMatrix(A)' * MyMatrix(A) .=== apply(*,MyMatrix(A)', MyMatrix(A)))
        @test all(MyMatrix(A)' * A .=== apply(*,MyMatrix(A)', A))
        @test all(MyMatrix(A) * MyMatrix(A)' .=== apply(*,MyMatrix(A), MyMatrix(A)'))
        @test all(A * MyMatrix(A)' .=== apply(*,A, MyMatrix(A)'))

        @test all(UpperTriangular(A) * MyMatrix(A) .=== apply(*,UpperTriangular(A), MyMatrix(A)))
        @test all(MyMatrix(A) * UpperTriangular(A) .=== apply(*, MyMatrix(A),UpperTriangular(A)))

        @test all(Diagonal(A) * MyMatrix(A)' .=== apply(*,Diagonal(A), MyMatrix(A)'))
        @test all(MyMatrix(A)' * Diagonal(A) .=== apply(*,MyMatrix(A)',Diagonal(A)))
        @test all(UpperTriangular(A) * MyMatrix(A)' .=== apply(*,UpperTriangular(A), MyMatrix(A)'))
        @test all(MyMatrix(A)' * UpperTriangular(A) .=== apply(*,MyMatrix(A)',UpperTriangular(A)))

        @test MyMatrix(A)\x ≈ apply(\,MyMatrix(A),x) ≈ copyto!(similar(x),Ldiv(A,copy(x))) ≈ A\x
        @test eltype(applied(\,MyMatrix(A),x)) == eltype(apply(\,MyMatrix(A),x)) == eltype(MyMatrix(A)\x) == Float64

        @test MyMatrix(A)\MyMatrix(B) ≈ MyMatrix(A)\B ≈ apply(\,MyMatrix(A),B) ≈ copyto!(similar(B),Ldiv(A,copy(B))) ≈ A\B
        @test eltype(applied(\,MyMatrix(A),B)) == eltype(apply(\,MyMatrix(A),B)) == eltype(MyMatrix(A)\B) == Float64

        @test MyMatrix(A) * ApplyArray(exp,B) ≈ apply(*, MyMatrix(A),ApplyArray(exp,B)) ≈ A*exp(B)
        @test ApplyArray(exp,A) * MyMatrix(B)  ≈ apply(*, ApplyArray(exp,A), MyMatrix(B)) ≈ exp(A)*B
        @test ApplyArray(exp,A) * ApplyArray(exp,B) ≈ apply(*, ApplyArray(exp,A),ApplyArray(exp,B)) ≈ exp(A)*exp(B)

        @test MyMatrix(A) * BroadcastArray(exp,B) ≈ apply(*, MyMatrix(A),BroadcastArray(exp,B)) ≈ A*exp.(B)
        @test BroadcastArray(exp,A) * MyMatrix(B)  ≈ apply(*, BroadcastArray(exp,A), MyMatrix(B)) ≈ exp.(A)*B
        @test BroadcastArray(exp,A) * BroadcastArray(exp,B) ≈ apply(*, BroadcastArray(exp,A),BroadcastArray(exp,B)) ≈ exp.(A)*exp.(B)
    end

    @testset "lmul!" begin
        A = randn(5,5)
        B = randn(5,5)
        x = randn(5)

        @test lmul!(MyUpperTriangular(A), copy(x)) ≈ MyUpperTriangular(A) * x
        @test lmul!(MyUpperTriangular(A), copy(B)) ≈ MyUpperTriangular(A) * B

        @test_skip lmul!(MyUpperTriangular(A),view(copy(B),collect(1:5),1:5)) ≈ MyUpperTriangular(A) * B
    end

    @testset "\\" begin
        A = randn(5,5)
        B = randn(5,5)
        x = randn(5)
        @test MyMatrix(A) \ x == apply(\, MyMatrix(A), x)
        @test ldiv!(MyMatrix(A), copy(x)) == materialize!(Ldiv(MyMatrix(A), copy(x)))
        @test MyMatrix(A) \ x ≈ ldiv!(MyMatrix(A), copy(x)) ≈ A\x
        @test MyMatrix(A) \ B == apply(\, MyMatrix(A), B)
        @test ldiv!(MyMatrix(A), copy(B)) == materialize!(Ldiv(MyMatrix(A), copy(B)))
        @test MyMatrix(A) \ B ≈ MyMatrix(A) \ MyMatrix(B) ≈ ldiv!(MyMatrix(A), copy(B)) ≈  A\B
        @test_broken ldiv!(MyMatrix(A), MyMatrix(copy(B))) ≈ A\B

        C = randn(5,3)
        @test all(MyMatrix(C)\x .=== apply(\,MyMatrix(C),x))
        @test MyMatrix(C)\x ≈ C\x
        @test all(MyMatrix(C)\B .=== apply(\,MyMatrix(C),B))
        @test MyMatrix(C)\B ≈ C\B

        @test_throws DimensionMismatch apply(\,MyMatrix(C),randn(4))
        @test_throws DimensionMismatch apply(\,MyMatrix(C),randn(4,3))
    end

    @testset "Lazy" begin
        A = MyLazyArray(randn(2,2))
        B = MyLazyArray(randn(2,2))
        x = MyLazyArray(randn(2))

        @test apply(*,A,x) isa ApplyVector
        @test apply(*,A,Array(x)) isa ApplyVector
        @test apply(*,Array(A),x) isa ApplyVector
        @test apply(*,A,x) ≈ apply(*,Array(A),x) ≈ apply(*,A,Array(x))  ≈ Array(A)*Array(x)

        @test apply(*,A,B) isa ApplyMatrix
        @test apply(*,A,Array(B)) isa ApplyMatrix
        @test apply(*,Array(A),B) isa ApplyMatrix
        @test apply(*,A,B) ≈ apply(*,Array(A),B) ≈ apply(*,A,Array(B))  ≈ Array(A)*Array(B)

        @test apply(\,A,x) isa ApplyVector
        @test apply(\,A,Array(x)) isa ApplyVector
        @test apply(\,Array(A),x) isa ApplyVector
        @test apply(\,A,x) ≈ apply(\,Array(A),x) ≈ apply(\,A,Array(x)) ≈ Array(A)\Array(x)

        @test apply(\,A,B) isa ApplyMatrix
        @test apply(\,A,Array(B)) isa ApplyMatrix
        @test apply(\,Array(A),B) isa ApplyMatrix
        @test apply(\,A,B) ≈ apply(\,Array(A),B) ≈ apply(\,A,Array(B)) ≈ Array(A)\Array(B)

        Ap = applied(*,A,x)
        @test Ap isa Applied{LazyArrayApplyStyle}
        @test copyto!(similar(Ap), Ap) == A*x
        @test copyto!(similar(Ap,BigFloat), Ap) == A*x

        @test MemoryLayout(typeof(Diagonal(x))) isa DiagonalLayout{LazyLayout}
        @test MemoryLayout(typeof(Diagonal(ApplyArray(+,x,x)))) isa DiagonalLayout{LazyLayout}
        @test MemoryLayout(typeof(Diagonal(1:6))) isa DiagonalLayout{UnknownLayout}

        @test MemoryLayout(typeof(A')) isa LazyLayout
        @test MemoryLayout(typeof(transpose(A))) isa LazyLayout
        @test MemoryLayout(typeof(view(A,1:2,1:2))) isa LazyLayout
        @test MemoryLayout(typeof(reshape(A,4))) isa LazyLayout
    end
end


text:
function WriteCSV(A::Assoc,fname)
#Because of potential memory issues, the Assoc will not be converted to dense matrix, but instead directly printed onto the file.

iostream = open(fname,"w")

#First write column
for c = 1:size(A.col,1)
    write(iostream,',') #Separator
    print(iostream,A.col[c])
end

write(iostream,"\n")

valMap = !isempty(A.val) #Check if val needs to be mapped.
#For each row write in row
for r = 1:size(A.row,1)
    print(iostream,A.row[r])
    for c = 1:size(A.col,1)
        write(iostream,',')
        if A.A[r,c] != 0
            if valMap #Mapping needed
                write(iostream,A.val[A.A[r,c]])
            else #Mappping not needed.
                print(iostream,A.A[r,c])
            end
            end
    end
    write(iostream,"\n")
    flush(iostream) #Clear buffer for each line.
end
#Close stream.
close(iostream)
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
using Jive
@useinside module test_stdlib_distributed_sshmanager

using Test
using Distributed
using .Distributed: SSHManager

const test_exeflags = Base.julia_cmd()
const test_exename = popfirst!(test_exeflags.exec)

function addprocs_with_testenv(manager; kwargs...)
    tunnel = false
    sshflags = ``
    max_parallel = 10
    pids = addprocs(manager; exename=test_exename, exeflags=test_exeflags, tunnel=tunnel, sshflags=sshflags, max_parallel=max_parallel, kwargs...)
    pids
end

curr_ws = workers()

machines = [("localhost", 2)]
manager = SSHManager(machines)
ws = try
    addprocs_with_testenv(manager)
catch err
    @test err isa CompositeException
    []
end

using .Distributed: Worker, MsgHeader, CallMsg, send_msg
function myremotecall(f, w::Worker, args...; kwargs...)
    rr = Future(w)
    send_msg(w, MsgHeader(remoteref_id(rr)), CallMsg{:call}(f, args, kwargs))
    return rr
end

@test Distributed.PGRP.workers[1] isa Distributed.LocalProcess

if !isempty(ws)
    wid = first(ws)
    w = Distributed.worker_from_id(wid)
    @test wid == worker_id_from_socket(w.w_stream)
    @test wid == worker_id_from_socket(w.r_stream)
    @test Distributed.map_sock_wrkr[w.r_stream] === w
    @test (fetch ∘ myremotecall)(() -> myid(), w) == wid
end

end # module test_stdlib_distributed_sshmanager


text:
module SSLTestDataFrames
using SimpleSDMLayers
using DataFrames
using Test

temperature = worldclim(1)

df = DataFrame(latitude = [0.0, 1.0], longitude = [30.0, 31.0], values = [42.0, 15.0])

@test eltype(temperature[df]) <: Number

temperature_clip = clip(temperature, df)

@test typeof(temperature_clip) == typeof(temperature)

@test typeof(DataFrame(temperature_clip)) == DataFrame
@test eltype(DataFrame(temperature_clip).values) == eltype(temperature_clip.grid)
@test typeof(DataFrame([temperature_clip, temperature_clip])) == DataFrame
@test eltype(DataFrame([temperature_clip, temperature_clip]).x1) == eltype(temperature_clip.grid)

@test typeof(SimpleSDMPredictor(df, :values, temperature_clip)) <: SimpleSDMLayer
@test typeof(SimpleSDMPredictor(df, :values, temperature_clip)) <: SimpleSDMPredictor

end


text:
using DataFrames, LinearAlgebra, MixedModels, RData, Test

if !@isdefined(dat) || !isa(dat, Dict{Symbol, DataFrame})
    const dat = Dict(Symbol(k) => v for (k, v) in
        load(joinpath(dirname(pathof(MixedModels)), "..", "test", "dat.rda")))
end

@testset "contra" begin
    contra = dat[:Contraception]
    contra[!, :urbdist] = categorical(string.(contra[!, :d], contra[!, :urb]))
    contraform = @formula(use ~ 1+a+abs2(a)+urb+l+(1|urbdist))
    gm0 = fit(MixedModel, contraform, contra, Bernoulli(), fast=true);
    @test gm0.lowerbd == zeros(1)
    @test isapprox(gm0.θ, [0.5720734451352923], atol=0.001)
    @test isapprox(deviance(gm0,true), 2361.657188518064, atol=0.001)
    gm1 = fit(MixedModel, contraform, contra, Bernoulli());
    @test isapprox(gm1.θ, [0.573054], atol=0.005)
    @test lowerbd(gm1) == vcat(fill(-Inf, 7), 0.)
    @test isapprox(deviance(gm1,true), 2361.54575, rtol=0.00001)
    @test isapprox(loglikelihood(gm1), -1180.77288, rtol=0.00001)
    @test dof(gm0) == length(gm0.β) + length(gm0.θ)
    @test nobs(gm0) == 1934
    fit!(gm0, fast=true, nAGQ=7)
    @test isapprox(deviance(gm0), 2360.9838, atol=0.001)
    fit!(gm0, nAGQ=7)
    @test isapprox(deviance(gm0), 2360.8760, atol=0.001)
    @test gm0.β == gm0.beta
    @test gm0.θ == gm0.theta
    @test isnan(gm0.σ)
    @test length(gm0.y) == size(gm0.X, 1)
    @test :θ in propertynames(gm0)
    gm0.β = gm0.beta
    @test gm0.β == gm0.beta
    gm0.θ = gm0.theta
    @test gm0.θ == gm0.theta
    gm0.βθ = vcat(gm0.β, gm0.theta)
    @test gm0.β == gm0.beta
    @test gm0.θ == gm0.theta
    # the next three values are not well defined in the optimization
    #@test isapprox(logdet(gm1), 75.7217, atol=0.1)
    #@test isapprox(sum(abs2, gm1.u[1]), 48.4747, atol=0.1)
    #@test isapprox(sum(gm1.resp.devresid), 2237.349, atol=0.1)
    show(IOBuffer(), gm1)
end

@testset "cbpp" begin
    cbpp = dat[:cbpp]
    cbpp[!, :prop] = cbpp[!, :i] ./ cbpp[!, :s]
    gm2 = fit(MixedModel, @formula(prop ~ 1 + p + (1|h)), cbpp, Binomial(), wts = cbpp[!,:s])
    @test isapprox(deviance(gm2,true), 100.09585619324639, atol=0.0001)
    @test isapprox(sum(abs2, gm2.u[1]), 9.723175126731014, atol=0.0001)
    @test isapprox(logdet(gm2), 16.90099, atol=0.0001)
    @test isapprox(sum(gm2.resp.devresid), 73.47179193718736, atol=0.001)
    @test isapprox(loglikelihood(gm2), -92.02628186555876, atol=0.001)
    @test isnan(sdest(gm2))
    @test varest(gm2) == 1
end

@testset "verbagg" begin
    gm3 = fit(MixedModel, @formula(r2 ~ 1 + a + g + b + s + (1|id)+(1|item)), dat[:VerbAgg],
         Bernoulli())
    @test deviance(gm3) ≈ 8151.40 rtol=1e-5
    @test lowerbd(gm3) == vcat(fill(-Inf, 6), zeros(2))
    @test fitted(gm3) == predict(gm3)
    # these two values are not well defined at the optimum
    @test sum(x -> sum(abs2, x), gm3.u) ≈ 273.31563469936697 rtol=1e-3
    @test sum(gm3.resp.devresid) ≈ 7156.558983084621 rtol=1e-4
end

@testset "grouseticks" begin
    gm4 = fit(MixedModel, @formula(t ~ 1 + y + ch + (1|i) + (1|b) + (1|l)),
              dat[:grouseticks], Poisson(), fast=true)  # fails in pirls! with fast=false
    @test isapprox(deviance(gm4), 851.4046, atol=0.001)
    @test lowerbd(gm4) == vcat(zeros(3))
    # these two values are not well defined at the optimum
    #@test isapprox(sum(x -> sum(abs2, x), gm4.u), 196.8695297987013, atol=0.1)
    #@test isapprox(sum(gm4.resp.devresid), 220.92685781326136, atol=0.1)
end


text:
# Internal facilities for fast random number generation

immutable RandIntSampler  # for generating Int samples in [0, K-1]
    a::Int
    Ku::UInt
    U::UInt

    @compat RandIntSampler(K::Int) = (Ku = UInt(K); new(1, Ku, div(typemax(UInt), Ku) * Ku))
    @compat RandIntSampler(a::Int, b::Int) = (Ku = UInt(b-a+1); new(a, Ku, div(typemax(UInt), Ku) * Ku))
end

function rand(s::RandIntSampler)
    x = rand(UInt)
    while x >= s.U
        x = rand(UInt)
    end
    @compat s.a + Int(rem(x, s.Ku))
end

randi(K::Int) = rand(RandIntSampler(K))
randi(a::Int, b::Int) = rand(RandIntSampler(a, b))