text:
################################################################################
#
# TWRI module
#
################################################################################

module TWRIdual

using LinearAlgebra
using JUDI.TimeModeling
using PyCall
using Distributed
using FFTW, DSP
using Optim

const R = Float32; export R
const C = ComplexF32; export C
const SPACE_ORDER = 8
const NB = 100

include("utils.jl")
include("gendata.jl")
include("data_filter.jl")
include("FWIFun.jl")
include("TWRIdualFun.jl")

end


text:
struct CPEGDensityParameters
    a::Float64
    b::Float64
    c::Float64
    d::Float64
    e::Float64
    f::Float64
    g::Float64
    h::Float64
    i::Float64
    function CPEGDensityParameters()
        # default values from a marsgram sample
        p = [ -4.001833776317166
              -0.015696377977412
              -0.129412788692219
               0.000199820058253
               0.010377518080309
               0.000043652882189
              -0.000539682362508
              -0.000000205086106
               0.000000874980179]
        new(p...)
    end
end

function density(p::CPEGDensityParameters, h::T) where T
   h = h / 1000

   # clamp in a forward diff friendly way
   if h > 125.0
       h = one(h)*125.0
   elseif h < 0.2
       h = one(h)*0.2
   end
   num = @evalpoly(h, p.a, p.c, p.e, p.g, p.i)
   den = @evalpoly(h, 1.0, p.b, p.d, p.f, p.h)
   exp(num/den)
end

ext: jl


text:
# WaveSim.jl, wave propagation simulator.
# Wave propagation computation.
# Christophe MEYER, 2016-2019
module WaveSim

using Parameters
using StaticArrays
using Statistics

export WaveSimParameters

@enum ApodizationShape Rect Hann

# Configuration
@with_kw struct WaveSimParameters
  # Transmit center frequency.
  tx_frequency::Float64 = 3_000_000.0  # [Hz]
  # Length of the transmit pulse.
  pulse_cycles::Float64 = 2  # [cycles]
  # Speed of sound.
  c::Float64 = 1540.0  # [m/s]
  # Resolution at which to divide the simulation time span.
  temporal_res::Float64 = 0.1e-6  # [s]
  # Time span to simulate.
  end_simulation_time::Float64 = 56.0e-6  # [s] starts at 0 s
  # Spatial resolution of the simulation.
  spatial_res::SVector{2} = @SVector [256, 512]  # [pixels]
  # Field of view, x then z. x=0 centered on aperture center, z=0 at aperture plane.
  fov::SVector{2} = @SVector [4e-2, 8e-2]  # [m]
  # Physical length of the transducer array.
  transducer_array_size::Float64 = 0.03  # [m]
  # Spacing between physical elements of transducer array.
  transducer_pitch::Float64 = 208e-6  # [m]
  # Active aperture size of transducer (centered)
  aperture_size::Float64 = 0.02  # [m]
  # Distance from transducer to focus point, Inf means plane wave.
  focus_depth::Float64 = 0.03  # [m]
  # Steering angle in azimuth.
  steer_angle::Float64 = 10.0  # [deg]
  # Shape of the transmit pulse.
  pulse_shape_func = phase -> cos(phase)
  # Shape of the transmit aperture apodization.
  apodization_shape::ApodizationShape = Rect
  # Display range.
  dbrange::Float64 = 40
  # Plot orientation: beam is horizontal or vertical.
  orientation::Symbol = :horizontal
end

# compute dependent parameters given global configuration
function init(trans_delays, sim_params)
  @unpack tx_frequency, transducer_pitch, spatial_res, temporal_res, end_simulation_time, fov, pulse_cycles, apodization_shape = sim_params
  n_transducers = length(trans_delays)
  x_transducers = [transducer_pitch * itrans for itrans in 1:n_transducers]  # [m]
  x_transducers .+= fov[1] / 2 - mean(extrema(x_transducers))  # center around FOV in x
  time_vec = collect(0.0:temporal_res:end_simulation_time)
  pulse_length = pulse_cycles / tx_frequency
  if apodization_shape == Rect
      apodization_vec = ones(n_transducers)
  else
      apodization_vec = [(sin(i*pi/(n_transducers-1)))^2 for i in 1:n_transducers]
  end
  return x_transducers, time_vec, pulse_length, apodization_vec
end

# simulate one time step of wave propagation
function simulate_one_time_step!(image, t, trans_delays, x_transducers, pulse_length, tx_frequency, c, spatial_res, fov, pulse_shape_func, apodization_vec)
  # println("simulate_one_time_step")
  image_pitch = fov ./ spatial_res
  transducers_that_are_firing = findall(trans_delays .>= 0)
  trans_coord = zeros(MVector{2})  # [m] space
  pix_coord = zeros(MVector{2})  # [m] space
  one_over_c = 1.0 / c  # For computation speed improvement [s/m]

  for y in 1:spatial_res[2]
    pix_coord[2] = y * image_pitch[2]  # [m]
    for x in 1:spatial_res[1]
      pix_coord[1] = x * image_pitch[1]  # [m]

      # for transducers that will fire...
      @inbounds for i_trans in transducers_that_are_firing
        xt = x_transducers[i_trans]  # index
        trans_delay = trans_delays[i_trans]
        trans_coord[1] = xt
        trans_coord[2] = 0.0
        # from pixel to transducer
        dist_to_transducer = sqrt((trans_coord[1] - pix_coord[1])^2 +
                                  (trans_coord[2] - pix_coord[2])^2)
        # wave is spreading spherically in space, energy spreading loss goes with r^2, amp with r
        # see: https://ccrma.stanford.edu/~jos/pasp/Spherical_Waves_Point_Source.html
        wave_spreading_factor = 1 / dist_to_transducer
        time_to_transducer = dist_to_transducer * one_over_c
        time_to_reach = time_to_transducer + trans_delay
        # if the transducer wave reached this pixel...
        if time_to_reach <= t <= time_to_reach + pulse_length
          amp = pulse_shape_func((t - time_to_reach) * tx_frequency * 2pi) * wave_spreading_factor
        else
          amp = 0.0
        end
        amp *= apodization_vec[i_trans]
        # wave interference
        image[x,y] += amp
        # TODO: maybe scale by pi/4 (Kirchhoff's formula)
      end
    end
  end
end

# run the simulation time steps
function wavesim(trans_delays, sim_params)
  @unpack tx_frequency, pulse_cycles, spatial_res = sim_params
  @unpack c, fov, pulse_shape_func = sim_params
  x_transducers, tvec, pulse_length, apodization_vec = init(trans_delays, sim_params)

  images = zeros(Float32, (spatial_res[1], spatial_res[2], length(tvec)))

  Threads.@threads for i_time in 1:length(tvec)
    t = tvec[i_time]
    image = view(images, :, :, i_time)
    simulate_one_time_step!(image, t, trans_delays, x_transducers, pulse_length, tx_frequency, c, spatial_res, fov, pulse_shape_func, apodization_vec)
  end

  return images
end

# Get the beam profile spatial map and transmit time of beam energy, where each pixel indicates the maximum energy that was seen at that place, and at what time.
function beam_energy_map_and_transmit_time_map(images, sim_params)
    @unpack temporal_res = sim_params

    maxval, linindices = findmax(images, dims=3)

    beam_energy_map = dropdims(maxval, dims=3)

    transmit_time_map = similar(beam_energy_map)

    for linind in eachindex(linindices)
        x, y, t = Tuple(CartesianIndices(images)[linindices[linind]])
        transmit_time_map[linind] = t * temporal_res
    end

    return beam_energy_map, transmit_time_map
end

# Compute transmit time delays for transducer elements given aperture_size [m] (impacts how many elements are firing) and focus depth [m] and steer angle [deg] (both impact amount of delay).
# Note: Steer angle is about center of active aperture.
# Since current implementation places the active aperture about the center of the physical aperture, steer angle is also about center of physical aperture.
# Focus equation adapted from: Tumsys, 2014, http://dx.doi.org/10.5755/j01.eee.20.3.3638
# Steer equation adapted from: Ramm, 1983, http://dx.doi.org/10.1109/TBME.1983.325149
function delays_from_focus_and_steer(sim_params)
    @unpack focus_depth, steer_angle, aperture_size, c, transducer_pitch, fov = sim_params

    num_elements = round(Int, aperture_size / transducer_pitch)

    x_transducers = [transducer_pitch * ielem for ielem in 1:num_elements]  # [m]
    x_transducers .+= fov[1] / 2 - mean(extrema(x_transducers))  # center around FOV in x

    trans_delays = zeros(num_elements)

    # Focus
    focus_coord = [0.0, focus_depth]
    if Inf == focus_depth
        # Focus at infinity is considered planewave. Zero delay is fine.
    else
        # Define variables as in the reference paper (see above).
        l0 = focus_depth
        A = aperture_size
        e = transducer_pitch
        n = num_elements
        v1 = c
        # Compute dependent variables.
        l02 = l0^2
        Ae = A - e
        first_part = sqrt(l02 + (Ae / 2) ^ 2)  # longest dist from focus to elem (i.e. to the edge)

        # for transducers that will fire...
        for (i_elem, x_elem) = enumerate(x_transducers)
            k = i_elem
            second_part = sqrt(l02 + (Ae * abs(n - 2k + 1) / (2n - 1)) ^ 2)  # offset_dist
            trans_delays[i_elem] = (first_part - second_part) / v1
        end
    end

    # Steer
    # Define variables as in the reference paper (see above).
    θ = deg2rad(steer_angle)
    d = transducer_pitch
    # n * d = x_elem
    # for transducers that will fire...
    for (i_elem, x_elem) = enumerate(x_transducers)
        trans_delays[i_elem] += x_elem / c * sin(θ)
    end

    # Causal delays
    trans_delays .-= minimum(trans_delays)

    return trans_delays
end

end  # module WaveSim


text:
using Documenter, SubspaceInference

makedocs(
    modules = [SubspaceInference],
    format = Documenter.HTML(),
    sitename = "SubspaceInference.jl",
    doctest = true,
    pages = [
        "Home" => "index.md",
        "Subspace Construction" => "subspace.md",
        "Inference" => "inference.md",
        "NN Example" => "nn_example.md",
        "Neural ODE Example" => "node_example.md",
    ]
)

deploydocs(
    repo = "github.com/efmanu/SubspaceInference.jl.git",
)


text:
if !isdefined(:PDMats)      using PDMats      end # Positive Definite Matrices
if !isdefined(:Optim)       using Optim       end # Optimizing over non-linear parameters
if !isdefined(:DiffBase)    using DiffBase    end # For accessing results of autodiff
if !isdefined(:ForwardDiff) using ForwardDiff end # For computing derivatives analytically
if !isdefined(:QuadGK)      using QuadGK      end # For 1-d integration

const default_tol_kepler_eqn = 1.e-8

function calc_rv_pal_one_planet{T,A<:AbstractArray{T,1}}( theta::A, time::Float64; tol::Real = default_tol_kepler_eqn, tref::Real = 300.0  )
  P = theta[1]
  K = theta[2]
  h = theta[3]
  k = theta[4]
  wpM0 = theta[5]
  ecc = sqrt(h*h+k*k)
  w = atan2(h,k)
  M0 = wpM0-w-2pi*tref/P
  n = 2pi/P
  M = time*n-M0
  #M = mod2pi(M)
  Mdiv2pi = floor(M/(2pi))
  M -= 2pi*Mdiv2pi
  lambda = w+M
  #E = ecc_anom_itterative_laguerre(M,ecc,tol=tol)   # Kepler equation solve integrated into this function due to limitations of ForwardDiff
  if(M<zero(M)) M += 2pi end
  E = (M<pi) ? M + 0.85*ecc : M - 0.85*ecc;
  const max_its_laguerre = 200
  local E_old
  for i in 1:max_its_laguerre
       E_old = E
       es = ecc*sin(E)
       ec = ecc*cos(E)
	   F = (E-es)-M
       Fp = 1.0-ec
       Fpp = es
       const n = 5
       root = sqrt(abs((n-1)*((n-1)*Fp*Fp-n*F*Fpp)))
       denom = Fp>zero(Fp) ? Fp+root : Fp-root
       E = E-n*F/denom

       if abs(E-E_old)<tol break end
  end
  @assert abs(E-E_old)<tol

  c = cos(lambda+ecc*sin(E))
  s = sin(lambda+ecc*sin(E))
  if ecc >= 1.0
    println("# ERROR in calc_rv_pal_one_planet: ecc>=1.0:  ",theta)
  end
  @assert(0.0<=ecc<1.0)
  j = sqrt((1.0-ecc)*(1.0+ecc))
  p, q = (ecc == 0.0) ? (zero(T), zero(T)) : (ecc*sin(E), ecc*cos(E))
  a = K/(n/sqrt((1.0-ecc)*(1.0+ecc)))
  zdot = a*n/(1.0-q)*( cos(lambda+p)-k*q/(1.0+j) )
end

function calc_model_rv{T,A<:AbstractArray{T,1}}( theta::A, time::Float64; tol::Real = default_tol_kepler_eqn )
  const nparam_per_pl = 5
  const nparam_non_pl = 1
  num_pl = convert(Int64,floor((length(theta)-nparam_non_pl)//nparam_per_pl))
  v = zero(T)
  for p in 1:num_pl
     P = theta[1+(p-1)*nparam_per_pl]
     K = theta[2+(p-1)*nparam_per_pl]
     h = theta[3+(p-1)*nparam_per_pl]
	 k = theta[4+(p-1)*nparam_per_pl]
	 ecc = h*h+k*k
	 if P>zero(P) && K>=zero(K) && ecc<one(ecc)
	    #v += calc_rv_pal_one_planet(theta[1+(p-1)*nparam_per_pl:p*nparam_per_pl],time)
		v += calc_rv_pal_one_planet(view(theta,1+(p-1)*nparam_per_pl:p*nparam_per_pl),time)
     end
  end
  rv_offset = theta[num_pl*nparam_per_pl+1]
  v += rv_offset
  return v
end

function calc_model_rv{T,A<:AbstractArray{T,1}}( theta::A, times::Array{Float64}; tol::Real = default_tol_kepler_eqn )
  v = map(t->calc_model_rv(theta,t),times)
end

#=
time = 0.0
f(x) = calc_rv_pal_one_planet(x,time)
param = [12.0,5.0,0.1,0.1,pi/4]
f(param)
ForwardDiff.gradient(f,param)
cfg = ForwardDiff.GradientConfig(f, param)
ForwardDiff.gradient(f,param,cfg)
result = DiffBase.GradientResult(param)
ForwardDiff.gradient!(result,f,param,cfg)
@time ForwardDiff.gradient(f,param)
@time ForwardDiff.gradient(f,param,cfg)
@time ForwardDiff.gradient!(result,f,param,cfg)
f(x) = calc_model_rv(x,time)
param = [12.0,5.0,0.1,0.1,pi/4,  30.0,2.0,0.1,0.1,pi/4,  0.0, 0.0]
f(param)

times = collect(linspace(0.0,100.0,10));
f(x) = calc_model_rv(x,times)
ForwardDiff.gradient(f,param)
ForwardDiff.jacobian(f,param)
=#

function calc_model_rv_deltas{T}(theta::Array{T,1}, times::Array{Float64}, vels::Array{Float64})
  const nparam_per_pl = 5
  const nparam_non_pl = 1
  num_pl = convert(Int64,floor((length(theta)-nparam_non_pl)//nparam_per_pl))
  logp = zero(T)
  for p in 1:num_pl
     P = theta[1+(p-1)*nparam_per_pl]
	 if P<0.0
	   logp -= 1000*length(times)*P^2
	 end
     K = theta[2+(p-1)*nparam_per_pl]
	 if K<0.0
	   logp -= 1000*length(times)*K^2
	 end
	 h = theta[3+(p-1)*nparam_per_pl]
	 k = theta[4+(p-1)*nparam_per_pl]
	 ecc = h*h+k*k
	 if ecc>=1.0
	   logp -= 1000*length(times)*(ecc-1)^2
	 end
  end
  delta = calc_model_rv(view(theta,1:num_pl*nparam_per_pl+1),times).-vels
  (delta,logp)
end
function loglikelihood_fixed_jitter_from_delta_Sigma{T,PDMatT<:AbstractPDMat}(delta::Array{T,1}, Sigma::PDMatT)
  chisq = dot(delta, Sigma \ delta)
  -0.5*(chisq + logdet(Sigma) +length(delta)*log(2pi) )
end
function loglikelihood_fixed_jitter{T,PDMatT<:AbstractPDMat}(theta::Array{T,1}, times::Array{Float64}, vels::Array{Float64}, Sigma::PDMatT)
   (delta,logp) = calc_model_rv_deltas(theta,times,vels)
   logp += loglikelihood_fixed_jitter_from_delta_Sigma(delta,Sigma)
end

function loglikelihood_marinalize_jitter{T}(theta::Array{T,1}, times::Array{Float64}, vels::Array{Float64}; sigma_obs::Array{Float64}=ones(length(vels)), num_Sigmas::Integer=20, Sigma_cache=make_Sigma_cache(num_Sigmas,PDMat(make_Sigma(times,sigma_obs,0.))) )
	(sigmaj_list, weights, SigmaList) = Sigma_cache
	(delta,logp) = calc_model_rv_deltas(theta,times,vels)
	marginalized_likelihood_fix_periods_sigmaj = map(k-> loglikelihood_fixed_jitter_from_delta_Sigma(delta,SigmaList[k]),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 loglikelihood_add_whitenoise{T,PDMatT<:AbstractPDMat}(theta::Array{T,1}, times::Array{Float64}, vels::Array{Float64}, Sigma::PDMatT)
  const nparam_per_pl = 5
  const nparam_non_pl = 2
  num_pl = convert(Int64,floor((length(theta)-nparam_non_pl)//nparam_per_pl))
  logp = 0
  for p in 1:num_pl
     P = theta[1+(p-1)*nparam_per_pl]
	 if P<0.0
	   logp -= 1000*length(times)*P^2
	 end
     K = theta[2+(p-1)*nparam_per_pl]
	 if K<0.0
	   logp -= 1000*length(times)*K^2
	 end
	 h = theta[3+(p-1)*nparam_per_pl]
	 k = theta[4+(p-1)*nparam_per_pl]
	 ecc = h*h+k*k
	 if ecc>=1.0
	   logp -= 1000*length(times)*(ecc-1)^2
	 end
  end
  delta = calc_model_rv(view(theta,1:(length(theta)-1)),times).-vels
  sigmaj = theta[end]
  SigmaNew = Sigma + PDiagMat(sigmaj^2*ones(length(times)))
  chisq = dot(delta, SigmaNew \ delta)
  logp -= 0.5*(chisq + logdet(SigmaNew) + length(times)*log(2pi) )
  logp += logprior_kepler(theta)
  return logp
end
=#

#=
# 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(sigma_j_0)-log(log1p(sigma_j_max/sigma_j_0))
end

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

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

=#

arg_2d(h::Real,k::Real,sigma::Real) = exp(-0.5*(h^2+k^2)/sigma^2)
arg_1d(h::Real,sigma::Real) = QuadGK.quadgk(k->arg_2d(h,k,sigma),-sqrt(1-h^2),sqrt(1-h^2))[1]
compute_prior_ecc_norm(sigma::Real) = 1/QuadGK.quadgk(h->arg_1d(h,sigma),-1.0,1.0)[1]
sigma_ecc = 0.2
logprior_eccentricity_normalization = log(compute_prior_ecc_norm(sigma_ecc))

function logprior_eccentricty{T,A<:AbstractArray{T,1}}( theta::A )
  # logp = 0 # for e~Uniform[0,1)
  # For e~ Rayleigh(sigma_ecc) Truncated to unit disk
  h = theta[3]
  k = theta[4]
  logp = -0.5*(h^2+k^2)/sigma_ecc^2+ logprior_eccentricity_normalization
end

function logprior_eccentricty{T}( ecc::T )
  # logp = 0 # for e~Uniform[0,1)
  # For e~ Rayleigh(sigma_ecc) Truncated to unit disk
  logp = -0.5*(ecc^2)/sigma_ecc^2+ logprior_eccentricity_normalization +log(2pi)  # remove omega prior
end

function logprior_kepler{T,A<:AbstractArray{T,1}}( theta::A )
  const nparam_per_pl = 5
  const nparam_non_pl = 1
  num_pl = convert(Int64,floor((length(theta)-nparam_non_pl)//nparam_per_pl))
  @assert(num_pl>=0)
  logp = -2*num_pl*log(2pi)
  for p in 0:(num_pl-1)
    logp += logprior_period(theta[1+p*nparam_per_pl])
    logp += logprior_amplitude(theta[2+p*nparam_per_pl])
    logp += logprior_eccentricity(view(theta,1+p*nparam_per_pl:(p+1)*nparam_per_pl))
  end
  logp  += logprior_offset(theta[num_pl*nparam_per_pl+1])
  if length(theta) == num_pl*nparam_per_pl+2
    logp  += logprior_jitter(theta[num_pl*nparam_per_pl+2])
  end
  return logp
end

#=
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)
  ScalMat(length(t),2*sigma_j)
end

=#

function make_matrix_posdef_dumb(X::Array{T,2}) where T<:Real
   eps = 1e-9
   while !isposdef(X)
      X = X + ScalMat(size(X,1),eps)
	  eps *= 10
   end
   PDMat(X)
end

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

function make_matrix_posdef_softabs(hessian::Array{T,2}; alpha::Float64 = 1000.0) where T<:Real
  lam, Q = eig(-hessian)
  lam_twig = lam./tanh.(alpha*lam)
  H_twig_try = full(Hermitian(Q*diagm(lam_twig)*Q'))
  if isposdef(H_twig_try)
     H_twig = PDMat(H_twig_try)
  else
     println("# Problem in softabs: ",hessian)
	 println("# eigvals_orig =", lam)
	 println("# eigvals_new = ",eigvals(H_twig_try))
	 println("# diff = ",H_twig_try.-hessian)
     H_twig = PDMat(H_twig_try)
  end
  return H_twig
end
make_matrix_posdef = make_matrix_posdef_softabs
#make_matrix_posdef = make_matrix_pd_eigval

function laplace_approximation_kepler_model{T,PDMatT<:AbstractPDMat}(theta::Array{T,1}, times::Array{Float64}, vels::Array{Float64}, Sigma::PDMatT, loglikelihood::Function; report_log_MAP::Bool = false, verbose::Integer = 3 )
    f(x) = - ( loglikelihood(x, times, vels, Sigma)+logprior_kepler(theta) )
	td = TwiceDifferentiable(f, theta; autodiff = :forward)
	if verbose >= 2
	   optim_opts =  Optim.Options(show_trace=true,f_tol=1e-4,time_limit=60)
	else
	   optim_opts =  Optim.Options(show_trace=false,f_tol=1e-4,time_limit=60)
	end
    res_opt = optimize(td,theta,optim_opts)
    bestfit = res_opt.minimizer
	logp_max = -res_opt.minimum
	if verbose >= 1
	   println("# best fit: ",logp_max/log(10), " : ",bestfit)
	end
    const nparam_per_pl = 5
    const nparam_non_pl = 1
    num_pl = convert(Int64,floor((length(theta)-nparam_non_pl)//nparam_per_pl))
    const nparam_no_jitter =  num_pl*nparam_per_pl+nparam_non_pl
	if report_log_MAP
	   if length(theta) > nparam_no_jitter
	      sigmaj = theta[end]
	      logp_max += logprior_jitter(sigmaj)
	   end
	   return (logp_max,res_opt)
	end
    hess_cfg = ForwardDiff.HessianConfig(f, view(theta,1:nparam_no_jitter))
	hess_result = DiffBase.HessianResult(view(theta,1:nparam_no_jitter))
    #ForwardDiff.hessian!(hess_result,f,view(bestfit,1:nparam_no_jitter))
	ForwardDiff.hessian!(hess_result,f,bestfit[1:nparam_no_jitter])

	logdet_hess = logdet(make_matrix_posdef(DiffBase.hessian(hess_result)))
    if verbose >= 3
 	  println("# hessian = ",DiffBase.hessian(hess_result))
	  println("# logdet(hessian) = ",logdet_hess)
	end
	logmarginal = 0.5*(length(view(theta,1:nparam_no_jitter))*log(2pi)-logdet_hess)+logp_max
	return (logmarginal,res_opt)
end

function laplace_approximation_kepler_model_jitter_separate{T,PDMatT<:AbstractPDMat}(theta::Array{T,1}, times::Array{Float64}, vels::Array{Float64}, Sigma::PDMatT; sigmaj_min::Real = 0.0, sigmaj_max::Real = 10.0, opt_jitter::Bool=true, report_log_MAP::Bool = false)
                                                 # delta_sigmaj_factor::Real = 0.1 )
    rms = std(calc_model_rv(theta,times).-vels)
	if opt_jitter
    SigmaTmp = Sigma + PDiagMat(rms^2*ones(length(times)))
	res_opt = laplace_approximation_kepler_model(theta,times,vels,SigmaTmp,loglikelihood_fixed_jitter,report_log_MAP=report_log_MAP)[2]
	theta = res_opt.minimizer
	function helper(sigmaj)
     SigmaTmp = Sigma + PDiagMat(sigmaj^2*ones(length(times)))
	 -(logprior_jitter(sigmaj)+laplace_approximation_kepler_model(theta,times,vels,SigmaTmp,loglikelihood_fixed_jitter,report_log_MAP=report_log_MAP)[1])
	 #-(laplace_approximation_kepler_model(theta,times,vels,SigmaTmp,loglikelihood_fixed_jitter)[1])
	end
	res_opt = optimize(helper,sigmaj_min,sigmaj_max,show_trace=false)
	sigmaj = res_opt.minimizer
	println("# Best sigma_j = ",sigmaj)
	SigmaTmp = Sigma + PDiagMat(sigmaj^2*ones(length(times)))
	else
	   sigmaj = theta[end]
	   SigmaTmp = Sigma
    end
	(logp, res_opt) = laplace_approximation_kepler_model(theta,times,vels,SigmaTmp,loglikelihood_fixed_jitter,report_log_MAP=report_log_MAP)
    #logp += logprior_jitter(sigmaj)
	if report_log_MAP
	  return (logp, res_opt)
	end
	theta = res_opt.minimizer
	dSigmadJ = diagm(2*sigmaj*ones(length(times)))
	d2logpdJ2 = 0.5*trace((SigmaTmp \ dSigmadJ)^2)
	# println("# logdetjitter (ana) = ",logdetjitter)
    #=
    delta_sigmaj = delta_sigmaj_factor*sigmaj
	SigmaTmp = Sigma + PDiagMat((sigmaj+delta_sigmaj)^2*ones(length(times)))
	logp_hi = laplace_approximation_kepler_model(theta,times,vels,SigmaTmp,loglikelihood_fixed_jitter)[1]
	logp_hi += logprior_jitter(sigmaj+delta_sigmaj)
	SigmaTmp = Sigma + PDiagMat((sigmaj-delta_sigmaj)^2*ones(length(times)))
    logp_lo = laplace_approximation_kepler_model(theta,times,vels,SigmaTmp,loglikelihood_fixed_jitter)[1]
    logp_lo += logprior_jitter(sigmaj-delta_sigmaj)
	d2logpdJ2 = (logp_hi-2*logp+logp_lo)/delta_sigmaj^2
	println("# logp = ",logp)
	println("# logp_hi = ",logp_hi)
	println("# logp_lo = ",logp_lo)
    println("# d^2 logp/dsigmaj^2 = ",d2logpdJ2/log(10))
	println("# logdetjitter (num) = ",log(abs(d2logpdJ2)))
	=#
    log_int_over_sigma = logp + 0.5*(log(2pi)-log(abs(d2logpdJ2)))
	(log_int_over_sigma,res_opt)
end

function laplace_approximation_kepler_model_jitter_separate_gauss{T,PDMatT<:AbstractPDMat}(theta::Array{T,1}, times::Array{Float64}, vels::Array{Float64}, Sigma0::PDMatT, num_Sigmas::Integer = 30; Sigma_cache = make_Sigma_cache(num_Sigmas,Sigma0) )

	(sigmaj_list, weights, SigmaList) = Sigma_cache
	 X = ones(length(times),1)
     const nparam = size(X,2)

	   marginalized_likelihood_fix_periods_sigmaj = map(k->laplace_approximation_kepler_model(theta,times, vels, SigmaList[k],loglikelihood_fixed_jitter)[1],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

	#(log_int_over_sigma,res_opt)
end

#=
param = [12.0,15.0,0.1,0.1,pi/4,  30.0,10.0,-0.1,-0.1,pi*3/4,  0.0, 1.0]
nobs = 100
times = collect(linspace(0.0,100.0,nobs));
vels = calc_model_rv(param,times) + 1.0*randn(size(times))
Sigma = PDiagMat(ones(length(vels)))
loglikelihood_add_whitenoise(param, times, vels, Sigma)

laplace_approximation_kepler_model(param,times,vels,Sigma,loglikelihood_add_whitenoise)/log(10)
=#

#=
f(x) = -loglikelihood_add_whitenoise(param, times, vels, Sigma)
f(param)

grad_result = DiffBase.GradientResult(param);
hess_result = DiffBase.HessianResult(param);
grad_cfg = ForwardDiff.GradientConfig(f, param);
hess_cfg = ForwardDiff.HessianConfig(f, param);
ForwardDiff.gradient!(grad_result,f,param)
ForwardDiff.hessian!(hess_result,f,param)

@time f(param)
@time ForwardDiff.gradient!(grad_result,f,param)
@time ForwardDiff.hessian!(hess_result,f,param)

optim_opts = Optim.Options(show_trace=true,f_tol=1e-4);
optimize(f,param,optim_opts)

#od = OnceDifferentiable(f, param; autodiff = :forward)
#@time res_opt = optimize(od,param,optim_opts)

td = TwiceDifferentiable(f, param; autodiff = :forward)
@time res_opt = optimize(td,param,optim_opts)
@time ForwardDiff.hessian!(hess_result,f,res_opt.minimizer)
DiffBase.hessian(hess_result)
DiffBase.gradient(hess_result)
DiffBase.value(hess_result)

=#


text:
using PolygonIO
using Documenter

DocMeta.setdocmeta!(PolygonIO, :DocTestSetup, :(using PolygonIO); recursive=true)

makedocs(;
    modules=[PolygonIO],
    authors="PyDataBlog <bbrenyah@gmail.com> and contributors",
    repo="https://github.com/PyDataBlog/PolygonIO.jl/blob/{commit}{path}#{line}",
    sitename="PolygonIO.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://PyDataBlog.github.io/PolygonIO.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Reference" => "reference.md",
        "Stocks" => "stocks.md",
        "Crypto" => "crypto.md",
        "Forex" => "forex.md",
        "Streaming" => "streaming.md"
    ],
)

deploydocs(;
    repo="github.com/PyDataBlog/PolygonIO.jl",
)


text:
module TestRuns

using Test
using RLEVectors

@testset begin

    # numruns single
    @test numruns([1, 2, 2, 3, 3]) == 3
    @test numruns([1, 1, 2, 2, 3]) == 3

    # numruns pair
    @test numruns(["a", "b", "c", "d", "d", "e", "f", "g"], [3, 6, 9, 12, 15, 18, 21, 24]) == 7 # Run of 2 runvalues
    @test numruns(["a", "b", "c", "d", "d", "d", "e", "f"], [3, 6, 9, 12, 15, 18, 21, 24]) == 6 # Run of 3 runvalues
    @test numruns(["a", "b", "c", "d", "e", "f", "g"], [3, 6, 9, 12, 12, 15, 17]) == 6 # run of 2 runends
    @test numruns(["a", "b", "c", "d", "e", "f"], [3, 6, 12, 12, 12, 15]) == 4 # run of 3 runends
    @test numruns([1, 2, 3], [1, 2, 2]) == 2
    @test_throws ArgumentError numruns([1, 2, 3], [3, 2, 4])

    # one vector
    @test ree([1, 1, 2, 2, 3, 3, 4, 4, 4]) == ([1, 2, 3, 4], [2, 4, 6, 9])
    @test ree([1, 1, 2, 2, 3, 3, 4, 4, 7]) == ([1, 2, 3, 4, 7], [2, 4, 6, 8, 9])
    @test ree([0, 1, 2, 2, 3, 3, 4, 4, 7]) == ([0, 1, 2, 3, 4, 7], [1, 2, 4, 6, 8, 9])

    # Clean up an RLEVector
    @test ree(["a", "b", "c", "d", "d", "e", "f", "g"], [3, 6, 9, 12, 15, 18, 21, 24]) ==
          (["a", "b", "c", "d", "e", "f", "g"], [3, 6, 9, 15, 18, 21, 24]) # Run of 2 runvalues
    @test ree(["a", "b", "c", "d", "d", "d", "e", "f"], [3, 6, 9, 12, 15, 18, 21, 24]) ==
          (["a", "b", "c", "d", "e", "f"], [3, 6, 9, 18, 21, 24]) # Run of 3 runvalues
    @test ree(["a", "b", "c", "d", "e", "f", "g"], [3, 6, 9, 12, 12, 15, 17]) ==
          (["a", "b", "c", "d", "f", "g"], [3, 6, 9, 12, 15, 17]) # run of 2 runends
    @test ree(["a", "b", "c", "d", "e", "f"], [3, 6, 12, 12, 12, 15]) ==
          (["a", "b", "c", "f"], [3, 6, 12, 15]) # run of 3 runends
    @test ree([1, 2, 3], [1, 2, 2]) == ([1, 2], [1, 2])

    # two vectors
    @test ree([9], [1]) == ([9], [1])
    @test ree([4, 3], [2, 4]) == ([4, 3], [2, 4])
    @test ree([4, 3], [0, 4]) == ([3], [4])

    # scalar
    @test ree(5) == ([5], [1])
    @test ree([5]) == ([5], [1])

    # inverse_ree
    @test inverse_ree([4, 5, 6], [2, 5, 7]) == [4, 4, 5, 5, 5, 6, 6]
    @test inverse_ree(Real[], Int64[]) == Real[]
    @test_throws ArgumentError inverse_ree([], [1])

end # testset

end # module


text:
function build_arg_table()
    settings = ArgParseSettings(
        description="\033[32mNeuralNet\033[0m\n\n\n\n" *
        "\033[31m Simple Neural Net Framework for Binding Site Prediction t\033[0m"
    )

    @add_arg_table settings begin
        "predict"
            action = :command
            help = "run simple predictor"
        "simple_encoder"
            action = :command
            help = "run 8x3x8 encododer"
        "train_model"
            action = :command
            help = "train the neural net on sequence data"
    end

    @add_arg_table settings["simple_encoder"] begin
        "-N"
            help = "Number of iterations to train neural net"
            default = 1000
            arg_type = Int
        "-a"
            help = "Learning Rate"
            default = 1e-3
            arg_type = Float64
        "output_file"
            help = "File to output results to [.CSV, .TSV, etc]"
            required = true
    end

    settings["simple_encoder"].description = "Runs the simple 8x3x8 predictor. "

    @add_arg_table settings["train_model"] begin
        "-N"
            help = "Number of iterations to train neural net"
            default = 1000
            arg_type = Int
        "-a"
            help = "Learning Rate"
            default = 1e-3
            arg_type = Float64
        "--hl_size"
            help = "hidden layer size"
            arg_type = Int
            default = 3
        "-w"
            help = "Choose whether or not to write output"
            default = true
        "-b"
            help = "Class inbalance parameter"
            arg_type = Int
            default = 1
        "-C"
            help = "Number of cross-validation rounds"
            arg_type = Int
            default = 10
        "pos_seq"
            help = "data set of labeled sites"
            required = true
        "total_seq"
            help = "Fasta file for predictions"
            required = true
        "output_file"
            help = "File to output results to [.CSV, .TSV, etc]"
            required = true
    end

    settings["train_model"].description = "Trains the neural net for testing purposes."

    @add_arg_table settings["predict"] begin
        "-N"
            help = "Number of iterations to train neural net"
            default = 1000
            arg_type = Int
        "-a"
            help = "Learning Rate"
            default = 1e-3
            arg_type = Float64
        "--hl_size"
            help = "hidden layer size"
            arg_type = Int
            default = 3
        "-w"
            help = "Choose whether or not to write output"
            default = true
        "-b"
            help = "Class inbalance parameter"
            arg_type = Int
            default = 1
        "-C"
            help = "Number of cross-validation rounds"
            arg_type = Int
            default = 10
        "pos_seq"
            help = "data set of labeled sites"
            required = true
        "total_seq"
            help = "Fasta file for predictions"
            required = true
        "predict_seq"
            help = "Contains the sequences you wish to predict"
            required = true
        "output_file"
            help = "File to output results to [.CSV, .TSV, etc]"
            required = true
    end

    settings["train_model"].description = "Trains the neural net for testing purposes."

    if typeof(Base.source_dir()) != Void
        settings.epilog = readstring(normpath(joinpath(Base.source_dir(),"..","LICENSE.md")))
    end

    return settings
end


text:
module AlgebraicDependencies

export @deps, dependencies

using LinearAlgebra
using SymEngine
using SymPy
using Nemo
using Hecke
using Singular
using Polynomials.PolyCompat

import Base: length, iterate, collect

include("ContinuedFractions.jl")

"""
    @deps r1 r2 r3 ...

Compute a basis for the algebraic dependencies among the ``r_i^n``.

# Examples
```julia-repl
julia> @deps 2 1/2
1-element Array{Expr,1}:
 :(v1 * v2 - 1)

julia> @deps (1-sqrt(5))/2 (1+sqrt(5))/2
1-element Array{Expr,1}:
 :(v1 ^ 2 * v2 ^ 2 - 1)
```
"""
macro deps(args...)
    dependencies(collect(args))
end

function dependencies(::Type{sideal}, roots::Vector{Basic}; variables=Basic[])
    if length(roots) < 2
        return nothing
    end
    if isempty(variables)
        variables = [Basic("v$i") for i in 1:length(roots)]
    elseif length(variables) != length(roots)
        throw("Number of variables does not match number exponentials. Got $(length(variables)), need $(length(roots))")
    end
    lattice = findrelations(roots)
    @debug "Lattice" lattice
    if isempty(lattice)
        return nothing
    end
    return ideal(lattice, variables)
end

dependencies(::Type{sideal}, roots; variables=[]) =
    dependencies(sideal, map(Basic ∘ string, roots); variables = map(Basic ∘ string, variables))

function dependencies(::Type{Basic}, roots::Vector{Basic}; variables=Basic[])
    ideal = dependencies(sideal, roots; variables = variables)
    if ideal === nothing
        return Basic[]
    end
    [Basic(string(ideal[i])) for i in 1:Singular.ngens(ideal)]
end

function dependencies(roots; variables=[])
    ideal = dependencies(sideal, map(Basic ∘ string, roots); variables = map(Basic ∘ string, variables))
    if ideal === nothing
        return Expr[]
    end
    [Meta.parse(string(ideal[i])) for i in 1:Singular.ngens(ideal)]
end

is_rational(x::Basic) = SymEngine.BasicType(x) isa SymEngine.BasicType{Val{:Integer}} || SymEngine.BasicType(x) isa SymEngine.BasicType{Val{:Rational}}

function minpoly(r::Basic, z::Symbol)
    s = sympify(z)
    p = sympy.minpoly(sympify(string(r)), s)
    cs = sympy.Poly(p, s).coeffs()
    bs = map(Basic ∘ string, cs)
    Basic(string(p)), PolyCompat.Poly(bs)
end

function common_number_field(roots::Vector{Basic})
    v = :a
    minpolys = Basic[]
    ps = PolyCompat.Poly[]
    for r in roots
        mp, p = minpoly(r, v)
        if is_rational(r)
            push!(minpolys, r)
        else
            push!(minpolys, mp)
        end
        push!(ps, p)
    end
    indices = findall(is_rational, minpolys)
    ratroots = [minpolys[i] for i in indices]
    deleteat!(minpolys, indices)
    mps = [convert(Expr, b) for b in minpolys]
    if isempty(mps)
        return 1, roots, ps
    end
    deleteat!(roots, indices)

    uniquemps = Base.unique(mps)
    mpsindices = [findfirst(x->x==m, uniquemps) for m in mps]
    R, _ = Nemo.PolynomialRing(Nemo.QQ, string(v))
    K, g = number_field(map(R, uniquemps))
    @debug "Number field" K, mpsindices
    S, mS = simple_extension(K)
    @debug "Simple extension" S
    gs = [SymEngine.lambdify(Basic(replace(string(mS\(x)), "//"=>"/"))) for x in g]
    rs = [gs[j](roots[i]) for (i, j) in enumerate(mpsindices)]
    for (i, ind) in enumerate(indices)
        splice!(rs, ind:ind-1, ratroots[i])
    end
    Hecke.degree(S), rs, ps
end

function masser_bound(degree::Int, roots::Vector{Basic}, mpolys::Vector{PolyCompat.Poly})
    k = length(roots)
    # assume all roots belong to the same number field
    d = degree
    # h = maximum height of the a[i]. The height of an algebraic number is the
    # sum of the degree and the binary length of all coefficients in the
    # defining equation over Q
    h = 0
    for p in mpolys
        cs = filter(x->!iszero(x), PolyCompat.coeffs(p))
        h0 = ceil(PolyCompat.degree(p) + sum([log(abs(c)) for c in cs]))
        if h0 > h
            h = h0
        end
    end
    @debug "Parameters for Masser bound" d h k
    return ceil(d^2 * (4*h*k*d* (log(2+d)/log(log(2+d)))^3)^(k-1) + 1)
end

gensym_unhashed(s::Symbol) = Symbol(replace(string(gensym(s)), "#"=>""))

function ideal(m::Matrix{BigInt}, x::Array{Basic,1})
    if isempty(m)
        return []
    end

    y = [Basic(gensym_unhashed(:vv)) for i in 1:length(x)]
    base = Basic[]
    for i in 1:nrows(m)
        r = 1
        for j in 1:ncols(m)
            exp = m[i,j]
            v = exp > 0 ? x[j] : y[j]
            r *= v^abs(exp)
        end
        base = [base; r - 1]
    end

    inv = [xi*yi - 1 for (xi,yi) in zip(x,y)]
    base = [base; inv]

    R, g = Singular.PolynomialRing(Nemo.QQ, map(string, [y; x]))
    sbasis = [R(convert(Expr, b)) for b in base]
    I = Singular.Ideal(R, sbasis)
    return Singular.eliminate(I, g[1:length(y)]...)
end

function findrelations(roots::Vector{Basic})
    # first treat zeros in the root list
    zpos = findall(iszero, roots)
    if length(zpos) == length(roots)
        return []
    end
    if !isempty(zpos)
        b = findrelations(filter(!iszero, roots))
        for i in zpos
            b = hcat(b[:,1:i-1], zeros(nrows(b)), b[:,i:end])
        end
        return b
    end

    degree, roots, mpolys = common_number_field(roots)
    @debug "Common number field" degree roots mpolys
    rs = [convert(Complex{BigFloat}, x) for x in roots]
    @debug "Complex roots" rs

    relog = [real(log(x)) for x in rs]
    imlog = [imag(log(x)) for x in rs]
    imlog = [imlog; 2*BigFloat(pi)]

    # Replace implicit zeros by explicit ones
    for i in 1:length(rs)
        z = rs[i]

        # abs(z) == 1
        if abs(abs(z) - 1) < 0.1 && abs(z) == 1
            relog[i] = 0
        end

        # z is real and z >= 0
        if isreal(z) && isreal(sqrt(z))
            imlog[i] = 0
        else
            # TODO: try harder: If[ Element[RootReduce[z], Reals] && Element[Sqrt[RootReduce[z]], Reals], imLog[[i]] = 0 ];
        end
    end
    @debug "reLog and imLog" relog imlog

    # Compute a bound for the coefficients of the generators
    bound = Int(masser_bound(degree, roots, mpolys))
    @debug "Masser bound" bound

    m = eye(Int, length(roots)) # identity matrix

    # Successively refine the approximation until only valid generators are returned
    level = Int(ceil(log2(bound) + 1))
    while prod(Bool[check_relation(roots, m[i,:]) for i in 1:size(m)[1]]) == 0
        m1 = getbasis(relog, level, bound)
        m2 = getbasis(imlog, level, bound)
        @debug "" m1 m2[:,1:end-1]
        m = z_module_intersect(m1, m2[:,1:end-1])
        level = level + 1
    end

    return m
end

hnf_with_transform(m::Matrix{Int}) = Matrix{Int}.(Nemo.hnf_with_transform(matrix(FlintZZ, m)))

function lattice_divide(l::Matrix{Int}, d::Int)
    n = ncols(l)
    return z_module_intersect(l, d*eye(Int, n)) / d
end

function check_relation(a::Vector{Basic}, e::Vector{<:Integer})
    return SymEngine.expand(prod([ax^ex for (ax, ex) in zip(a,e)])) == 1
end

function convergent(x, n)
    cf = ContinuedFraction(x)
    co = convergents(cf)
    return iterate(co, n)[1]
end

nrows(a::Matrix{<:Any}) = size(a)[1]
ncols(a::Matrix{<:Any}) = size(a)[2]

function getbasis(l::Vector{BigFloat}, level::Int, bound::Int)
    n = length(l)

    # First treat zero elements in l as special case
    zpos = findall(iszero, l)
    if length(zpos) == length(l)
        return eye(Int, n)
    end

    if length(zpos) > 0
        ll = deleteat!(copy(l), zpos)
        b = getbasis(ll, level, bound) # basis for nonzero numbers
        zvec = zeros(nrows(b), 1)
        # Insert new dimensions
        for pos in zpos
            b = hcat(b[:,1:pos-1], zvec, b[:,pos:end])
        end
        # Add unit vectors at the zero positions
        b = vcat(b, eye(Int, n)[zpos,:])
        return Matrix{Int}(b)
    end

    # Now we assume that l is a list of nonzero real numbers
    c1 = [convergent(x, level) for x in l]
    c2 = [convergent(x, level+1) for x in l]

    e = [1//denominator(x1*x2) for (x1, x2) in zip(c1, c2)]
    # Refine the approximation such that all errors are smaller than the smallest
    # element of l in absolute value
    lev = level + 1
    while length(filter(x -> maximum(e) >= abs(x), c1)) > 0 && lev < level + 5
        ex = findall(x->x==maximum(e), e) # indices with greates error

        lev++
        for i in 1:length(ex)
            j = ex[i]
            c1[j] = c2[j]
            c2[j] = convergent(l[j], lev)
            e[j] = c1[j] == l[j] ? 0 : 1/denominator(c1[j]*c2[j])
        end
    end

    # now: max e[i] < min |c1[i]|
    # this bound guarantees that generators whose norm is at most bound will
    # appear in the LLL-reduced basis
    minc1 = minimum([abs(c) for c in c1])
    maxe = maximum([abs(c) for c in e])
    d = BigInt(ceil(2^((length(c1) - 1)/2)*bound/(minc1 - maxe)))
    identity = eye(Int, n)
    row = c1 * d
    b = hcat(identity, row)
    b = lll(b)
    # Vectors whose right hand side is greater than the errors permit can be
    #   discarded; they cannot correspond to integer relations.
    # b = vcat([b[i,:] for i in 1:nrows(b) ])
    # TODO: find better way to filter rows
    res = Matrix{Rational{Int}}(undef, 0,n+1)
    for i in 1:nrows(b)
        if (abs(b[i,:][end]) <= d*sum(abs(bb*ee) for (bb, ee) in zip(b[i,1:n],e)))
            res = vcat(res, transpose(b[i,:]))
        end
    end

    # all remaining vectors are returned as candidates
    # TODO: result should be integer matrix?
    return Matrix{Int}(res[:,1:end-1])
end

function z_nullspace(matrix::Matrix{Int})
    h, t = hnf_with_transform(matrix)
    t = t * -1
    # println("HNF: $(h) | $(t)")

    # kernel is generated by the rows of t that correspond to zero rows in h
    zvec = zeros(size(h, 2))

    # TODO: find better way to filter zero vectors
    res = Matrix{Int}(undef, 0, ncols(t))
    for i in 1:nrows(t)
        if iszero(h[i,:])
            res = vcat(res, transpose(t[i,:]))
        end
    end
    return res
    # return [t[i,:] for i in 1:size(h, 1) if h[i,:] == zvec]
end

function z_module_intersect(base1::Matrix{Int}, base2::Matrix{Int})
    if isempty(base1) || isempty(base2)
        return []
    end

    sol = z_nullspace(vcat(base1, base2))

    if isempty(sol)
        return []
    end

    m1 = transpose(base1)
    m2 = transpose(sol[:, 1:nrows(base1)])
    return lll(collect(transpose(m1 * m2)))
end

function clear_denom(a::Matrix{Rational{BigInt}})
    d = lcm(denominator.(a))
    return a*d, d
end

import Nemo.lll

function lll(a::Matrix{Rational{BigInt}})
    m, d = clear_denom(a)
    m = numerator.(m)
    m = Matrix{BigInt}(Nemo.lll(matrix(FlintZZ, m)))
    return m // d
end

lll(m::Matrix{Int}) = Matrix{BigInt}(Nemo.lll(matrix(FlintZZ, m)))

eye(::Type{T}, n) where {T} = Matrix{T}(I, n, n)

function (R::FmpqPolyRing)(p::Expr)
    v = gen(R)
    vs = [:($(Symbol(string(v))) = $v)]
    q = quote
        let $(vs...)
            $(p)
        end
    end
    eval(q)
end

function (R::Singular.PolyRing)(p::Expr)
    vs = gens(R)
    vs = [:($(Symbol(string(v))) = $v) for v in vs]
    q = quote
        let $(vs...)
            $(p)
        end
    end
    eval(q)
end

function (K::NfAbsNS)(a::fmpq_poly)
    q, w = divrem(a, K.pol)
    z = NfAbsNSElem(w)
    z.parent = K
    return z
end

end # module


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 360367c3-a17d-4c04-b8eb-e2f0366baa86
import Pkg;

# ╔═╡ a69bf247-acbc-4b4d-889d-3f52ff046ef4
Pkg.activate("/store/users/ehinger/unfoldjl_dev/")

# ╔═╡ cf3dcd54-cfa0-11eb-2e99-756f9df71a4d
using Revise,Unfold,Plots,BSplines

# ╔═╡ a9e6971f-8a50-4564-ace5-1fabe65b1522
import PlutoUI

# ╔═╡ 8a62a281-d91d-4e70-b301-4ebafdd58b38

# ╔═╡ 539da607-f382-4eaa-93ad-31c906aa67ee

# ╔═╡ 17943406-ac61-49bd-bc8d-7d3f62441b8e

# ╔═╡ 89eace50-c0f5-4951-821f-cf010cdb9726
@bind nsplines PlutoUI.Slider(2:10; default=8, show_value=true)

# ╔═╡ cb62fb85-364d-4142-b28b-c8b0871b4ccc
@bind srate PlutoUI.Slider(1:50; default=10, show_value=true)

# ╔═╡ 8496a69b-20bc-45c2-92cb-179d0ad7127c
@bind τ1 PlutoUI.Slider(-5:0.1:2; default=-.5, show_value=true)

# ╔═╡ c8285ed7-bffb-4047-89ff-22a664c19a78
@bind τ2 PlutoUI.Slider(-2:0.1:5; default=2, show_value=true)

# ╔═╡ 9df0d50d-2a1a-4431-9725-12178c824d62
τ = [τ1,τ2]

# ╔═╡ 94a258ea-febd-4793-89c4-d3755e82c48c
begin
	function splinebasis(τ,sfreq,nsplines,name::String)
		    τ = Unfold.round_times(τ,sfreq)
    		times =range(τ[1],stop=τ[2],step=1 ./sfreq)
		    kernel=e->splinekernel(e,times,nsplines)
    		type = "splinebasis"

    		shiftOnset = Int64(floor(τ[1] * sfreq))

			return Unfold.BasisFunction(kernel,times,times,type,name,shiftOnset)
	end

	function spl_breakpoints(times,nsplines)
		# calculate the breakpoints, evenly spaced
		return collect(range(minimum(times),stop=maximum(times),length=nsplines))
	end
	function splinekernel(e,times,nsplines)
			breakpoints = spl_breakpoints(times,nsplines)
	   	 	basis = BSplineBasis(4,breakpoints)  # 4= cubic
			return Unfold.splFunction(times,basis)
	end
	spl=splinebasis(τ,srate,nsplines,"test")

end

# ╔═╡ 3f612e26-aba8-46ed-af19-67cd1be672fa
plot(spl.times,spl.kernel(0))

# ╔═╡ d375b64c-0e64-4b39-8c3f-5a14d365877d

# ╔═╡ 7086e800-d51a-4fe3-b33f-309acd676ae7

# ╔═╡ 4d6fc512-205e-4e25-b3b5-b7981d9dec75

# ╔═╡ ba9bafc2-9733-43d0-80cd-75b2637ba4a7

# ╔═╡ Cell order:
# ╠═360367c3-a17d-4c04-b8eb-e2f0366baa86
# ╠═a9e6971f-8a50-4564-ace5-1fabe65b1522
# ╠═a69bf247-acbc-4b4d-889d-3f52ff046ef4
# ╠═8a62a281-d91d-4e70-b301-4ebafdd58b38
# ╠═539da607-f382-4eaa-93ad-31c906aa67ee
# ╠═17943406-ac61-49bd-bc8d-7d3f62441b8e
# ╠═94a258ea-febd-4793-89c4-d3755e82c48c
# ╠═3f612e26-aba8-46ed-af19-67cd1be672fa
# ╠═89eace50-c0f5-4951-821f-cf010cdb9726
# ╠═cb62fb85-364d-4142-b28b-c8b0871b4ccc
# ╠═8496a69b-20bc-45c2-92cb-179d0ad7127c
# ╠═c8285ed7-bffb-4047-89ff-22a664c19a78
# ╠═9df0d50d-2a1a-4431-9725-12178c824d62
# ╠═d375b64c-0e64-4b39-8c3f-5a14d365877d
# ╠═cf3dcd54-cfa0-11eb-2e99-756f9df71a4d
# ╠═7086e800-d51a-4fe3-b33f-309acd676ae7
# ╠═4d6fc512-205e-4e25-b3b5-b7981d9dec75
# ╠═ba9bafc2-9733-43d0-80cd-75b2637ba4a7

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
function register_datadep()
    register(
        DataDep(
            "impute-v1.0.0",
            "Datasets for testing and demonstrating Impute.jl",
            "https://invenia-public-datasets.s3.amazonaws.com/Impute/v1.0.0/datasets.tar.gz",
            "938b3705752eb73141476a2abc7a36cfdaba9ec45f99f0796f44e0870e006e1c",
            post_fetch_method=unpack,
        )
    )
end

function datasets()
    dep = datadep"impute-v1.0.0/data/"

    # Only select paths containing a data.x file
    selected = Iterators.filter(walkdir(dep)) do (root, dirs, files)
        any(f -> first(splitext(f)) == "data", files)
    end

    # Return just the root path with the data dep path part removed
    return [first(t)[length(dep)+2:end] for t in selected]
end

function dataset(name)
    dep = @datadep_str joinpath("impute-v1.0.0/data", name)
    files = readdir(dep)
    idx = findfirst(f -> first(splitext(f)) == "data", files)
    idx === nothing && throw(ArguementError("No data file found for $name."))
    fullpath = joinpath(dep, files[idx])
    ext = splitext(fullpath)[end]

    # This is necessary because CSV isn't registered in FileIO
    if ext == ".csv"
        return CSV.File(fullpath)
    elseif ext == ".bson"
        return BSON.load(fullpath)
    else
        throw(ArgumentError("Unsupported file type $ext."))
    end
end


text:
#= convex.jl
Philip Barrett, pobarrett@gmail.com
01jul2016, Chicago

Tests the (two) convex hull routines to guarantee consistency
=#

Z = [ 1 -1; 1 1; -1 1; 0 -1 ; -1 -1 ]
    # Add a redundant point
a = Polygon( pts=Z ) ;
b=grahamScan(a)

c = Point{2, Float64}[Point(Z[i, 1], Z[i, 2]) for i=1:size(Z)[1]]
d = CHull2D.convexhull(c)

Y = [ 1 -1; 1 1; -1 1; 0 -1 + 1e-06 ; -1 -1 ]
e = Polygon( pts=Y ) ;
f = grahamScan(e)

g = Point{2, Float64}[Point(Y[i, 1], Y[i, 2]) for i=1:size(Y)[1]]
h = CHull2D.convexhull(g)


text:
struct BinaryConfusionMatrix{T<:Number} <: AbstractArray{T,2}
    true_positives::T
    false_positives::T
    false_negatives::T
    true_negatives::T
end

(::Type{BinaryConfusionMatrix{T}})() where {T} = BinaryConfusionMatrix{T}(0,0,0,0)
BinaryConfusionMatrix() = BinaryConfusionMatrix{Int}()
function BinaryConfusionMatrix(A::AbstractArray)
    @assert length(A) == 4
    BinaryConfusionMatrix{Int}(Tuple(A)...)
end

function Base.:(+)(a::BinaryConfusionMatrix{T},
                   b::BinaryConfusionMatrix{S}) where {T,S}
    BinaryConfusionMatrix{promote_type(T,S)}(
        a.true_positives  + b.true_positives,
        a.false_positives + b.false_positives,
        a.false_negatives + b.false_negatives,
        a.true_negatives  + b.true_negatives,
    )
end

@inline Core.Tuple(c::BinaryConfusionMatrix) = (c.true_positives, c.false_positives, c.false_negatives, c.true_negatives)
@inline Base.size(c::BinaryConfusionMatrix) = (2, 2)
function Base.getindex(c::BinaryConfusionMatrix, i::Int, j::Int)
    ifelse((i == 1) & (j == 1), c.true_positives,
    ifelse((i == 2) & (j == 1), c.false_positives,
    ifelse((i == 1) & (j == 2), c.false_negatives, c.true_negatives)))
end

true_positives(c::BinaryConfusionMatrix)  = c.true_positives
false_positives(c::BinaryConfusionMatrix) = c.false_positives
true_negatives(c::BinaryConfusionMatrix)  = c.true_negatives
false_negatives(c::BinaryConfusionMatrix) = c.false_negatives
condition_positive(c::BinaryConfusionMatrix) =
    true_positives(c) + false_negatives(c)
condition_negative(c::BinaryConfusionMatrix) =
    true_negatives(c) + false_positives(c)
LearnBase.nobs(c::BinaryConfusionMatrix) =
    condition_positive(c) + condition_negative(c)
predicted_positive(c::BinaryConfusionMatrix) =
    true_positives(c) + false_positives(c)
predicted_negative(c::BinaryConfusionMatrix) =
    true_negatives(c) + false_negatives(c)
correctly_classified(c::BinaryConfusionMatrix) =
    true_positives(c) + true_negatives(c)
incorrectly_classified(c::BinaryConfusionMatrix) =
    false_positives(c) + false_negatives(c)
prevalence(c::BinaryConfusionMatrix) =
    condition_positive(c) / nobs(c)

function Base.show(io::IO, ::MIME"text/plain", c::BinaryConfusionMatrix)
    println(io, summary(c), ":")
    tp = string(true_positives(c))
    fp = string(false_positives(c))
    fn = string(false_negatives(c))
    tn = string(true_negatives(c))
    len_p = max(maximum(map(length, (tp, fp))), 3)
    len_n = max(maximum(map(length, (fn, tn))), 3)
    tp = lpad(tp, len_p); fp = lpad(fp, len_p)
    fn = lpad(fn, len_n); tn = lpad(tn, len_n)
    pad = "  "
    println(io, pad, " ", "Predicted")
    println(io, pad, "  ", lpad("+",len_p), "   ", lpad("-",len_n))
    println(io, pad, "┌", repeat("─",len_p+2), "┬", repeat("─",len_n+2), "┐")
    println(io, pad, "│ ", tp, " │ ", fn, " │ +")
    println(io, pad, "├", repeat("─",len_p+2), "┼", repeat("─",len_n+2), "┤   Actual")
    println(io, pad, "│ ", fp, " │ ", tn, " │ -")
    println(io, pad, "└", repeat("─",len_p+2), "┴", repeat("─",len_n+2), "┘")
    println(io, pad, "   f1_score: ", round(f1_score(c),digits=4))
    println(io, pad, "   accuracy: ", round(accuracy(c),digits=4))
    println(io, pad, "  precision: ", round(precision_score(c),digits=4))
    println(io, pad, "sensitivity: ", round(sensitivity(c),digits=4))
    print(io, pad,   "specificity: ", round(specificity(c),digits=4))
end

# --------------------------------------------------------------------

function confusions(targets::AbstractArray{T},
                    outputs::AbstractArray,
                    encoding::AbstractVector) where T
    length(encoding) == 2 || error("Multiclass not yet implemented")
    confusions(targets, outputs, LabelEnc.NativeLabels{T,2}(encoding))
end

# Generic fallback. Tries to infer label encoding
function confusions(targets, outputs)
    encoding = _labelenc(targets, outputs)
    nlabel(encoding) == 2 || error("Multiclass not yet implemented")
    confusions(targets, outputs, encoding)
end

function confusions(targets::AbstractArray,
                    outputs::AbstractArray,
                    encoding::BinaryLabelEncoding)
    @_dimcheck length(targets) == length(outputs)
    tp::Int = 0; tn::Int = 0;
    fp::Int = 0; fn::Int = 0;
    @inbounds for I in eachindex(targets, outputs)
        target = targets[I]
        output = outputs[I]
        tp +=  true_positives(target, output, encoding)
        fp += false_positives(target, output, encoding)
        fn += false_negatives(target, output, encoding)
        tn +=  true_negatives(target, output, encoding)
    end
    BinaryConfusionMatrix{Int}(tp, fp, fn, tn)
end

function confusions(targets::AbstractArray,
                    outputs::AbstractArray,
                    encoding::BinaryLabelEncoding,
                    thresholds::AbstractVector)
    @_dimcheck length(targets) == length(outputs)
    curve = fill(BinaryConfusionMatrix{Int}(), length(thresholds))
    @inbounds for I in eachindex(targets, outputs)
        target = targets[I]
        output = outputs[I]
        for (j, th) in enumerate(thresholds)
            pred = convertlabel(encoding, classify(output, th), LabelEnc.ZeroOne())
            curve[j] += BinaryConfusionMatrix{Int}(
                 true_positives(target, pred, encoding),
                false_positives(target, pred, encoding),
                false_negatives(target, pred, encoding),
                 true_negatives(target, pred, encoding),
            )
        end
    end
    curve
end

# --------------------------------------------------------------------

struct ROCCurve{T<:Number,F<:Number,V<:AbstractVector{F}} <: AbstractArray{Tuple{F,BinaryConfusionMatrix{T}},1}
    thresholds::V
    confusions::Vector{BinaryConfusionMatrix{T}}
end

function ROCCurve(thresholds::AbstractVector, confusions::Vector{BinaryConfusionMatrix{T}}) where T
    @_dimcheck length(thresholds) == length(confusions)
    ROCCurve{T,eltype(thresholds),typeof(thresholds)}(thresholds, confusions)
end

function Base.:(+)(r1::ROCCurve, r2::ROCCurve)
    @assert r1.thresholds == r2.thresholds
    nm = r1.confusions .+ r2.confusions
    ROCCurve(r1.thresholds, nm)
end

Base.size(r::ROCCurve) = size(r.confusions)
Base.@propagate_inbounds function Base.getindex(r::ROCCurve, I::Int...)
    (r.thresholds[I...], r.confusions[I...])
end
Base.@propagate_inbounds function Base.setindex!(r::ROCCurve, tup::Tuple, I::Int...)
    r.thresholds[I...] = tup[1]
    r.confusions[I...] = tup[2]
end
Base.@propagate_inbounds function Base.setindex!(r::ROCCurve, v::BinaryConfusionMatrix, I::Int...)
    r.confusions[I...] = v
end

for fun in (:true_positives,  :true_negatives,
            :false_positives, :false_negatives,
            :condition_positive, :condition_negative,
            :correctly_classified, :incorrectly_classified,
            :prevalence,
            :predicted_positive, :predicted_negative,
            :true_positive_rate, :false_positive_rate,
            :true_negative_rate, :false_negative_rate,
            :positive_predictive_value, :negative_predictive_value,
            :false_discovery_rate, :false_omission_rate,
            :accuracy, :f_score, :f1_score,
            :positive_likelihood_ratio, :negative_likelihood_ratio,
            :diagnostic_odds_ratio)
    @eval ($fun)(r::ROCCurve) = Base.broadcast(($fun), r.confusions)
end

confusions(r::ROCCurve) = r.confusions

function confusions_at_sensitivity(r::ROCCurve, at::Number)
    0 <= at <= 1 || throw(ArgumentError("parameter \"at\" has to be in [0, 1]"))
    idx = findfirst(s -> s >= at, sensitivity(r))
    r[idx][2]
end

function specificity_at_sensitivity(r::ROCCurve, at::Number)
    specificity(confusions_at_sensitivity(r, at))
end

function precision_at_sensitivity(r::ROCCurve, at::Number)
    precision_score(confusions_at_sensitivity(r, at))
end

function accuracy_at_sensitivity(r::ROCCurve, at::Number)
    accuracy(confusions_at_sensitivity(r, at))
end

function Base.show(io::IO, r::ROCCurve)
    print(io, length(r), "-element ", typeof(r).name, " (auc: ", round(auc(r),digits=5), ")")
end

function Base.show(io::IO, ::MIME"text/plain", r::ROCCurve)
    println(io, summary(r), ":")
    fpr = false_positive_rate(r)
    tpr = true_positive_rate(r)
    p = lineplot(fpr, tpr, height=10, width=21, xlim=[0,1], ylim=[0,1], canvas=DotCanvas)
    annotate!(p, :r, 1, "auc: $(round(auc_from_rates(fpr, tpr),digits=5))")
    lineplot!(p, 0, 1, color=:white)
    xlabel!(p, "FPR")
    ylabel!(p, "TPR")
    print(io, p)
end

# --------------------------------------------------------------------

auc(r::ROCCurve) = auc_from_rates(false_positive_rate(r), true_positive_rate(r))
auc(args...; kw...) = auc(roc(args...; kw...))

function auc_from_rates(fpr::AbstractVector{T}, tpr::AbstractVector{R}) where {T,R}
    O = promote_type(T,R)
    @_dimcheck length(fpr) == length(tpr)
    dfpr = diff(fpr)
    dtpr = diff(tpr)
    a = zero(O); b = zero(O)
    @inbounds for i in 1:length(fpr)-1
        a += tpr[i] * dfpr[i]
        b += dtpr[i] * dfpr[i]
    end
    a + b / 2
end

# --------------------------------------------------------------------

@inline _thresholds(th::AbstractVector) = th
@inline _thresholds(th::Number) = range(1, stop=0, length=th)

function roc(targets::AbstractArray,
             outputs::AbstractArray,
             labels::AbstractVector;
             thresholds = 100)
    encoding = LabelEnc.NativeLabels(labels)
    nlabel(encoding) == 2 || error("Multiclass not yet implemented")
    roc(targets, outputs, encoding; thresholds=thresholds)
end

function roc(targets::AbstractArray,
             outputs::AbstractArray,
             encoding::BinaryLabelEncoding = labelenc(targets);
             thresholds = 100)
    th = _thresholds(thresholds)
    curve = confusions(targets, outputs, encoding, th)
    ROCCurve(th, curve)
end


text:
load("debug.jl")

module TestRecode
import Base.*
import Debug.*

code = :(function f(n)
    x = 0
    for k=1:n
        x += k
        println("k = ", k)
    end
    x = x*x
    let y=x
        z=2y
    end
    try
        error()
        x=1
    catch e
        n=2
        x
    end
end)

println(getdefs(code))

println()
println(code_debug(code))

end


text:
"""
    get_indices(arr::Containers.SparseAxisArray)

Get indices in SparseAxisArray.

    get_indices(arr)

Get indices in Array or DenseAxisArray.
"""
get_indices(arr::Containers.SparseAxisArray) = keys(arr.data)
get_indices(arr) = Iterators.product(axes(arr)...)

get_reform_param(param::Missing, args...; constr) =
    apply_interval_arithmetic(constr)
get_reform_param(param::Number, args...; kwargs...) = param
get_reform_param(param::Union{Vector,Tuple}, i::Int, args...; kwargs...) =
    get_reform_param(param[i], args...; kwargs...)
function get_reform_param(param::Dict, args...; kwargs...)
    arg_list = [arg for arg in args if !ismissing(arg)] #remove mising if j or k are missing
    get_reform_param(param[arg_list...]; kwargs...)
end

"""
    gen_constraint_name(constr)

Generate constraint name for a constraint to be split (Interval or EqualTo).
"""
function gen_constraint_name(constr)
    constr_name = name.(constr)
    if any(isempty.(constr_name))
        constr_name = gensym("constraint")
    elseif !isa(constr_name, String)
        c_names = union(first.(split.(constr_name,"[")))
        if length(c_names) == 1
            constr_name = c_names[1]
        else
            constr_name = gensym("constraint")
        end
    end

    return Symbol("$(constr_name)_split")
end

function replace_Symvars!(expr, model; logical_proposition = false)
    #replace JuMP variables with symbolic variables
    name = join(split(string(expr)," "))
    var = variable_by_name(model, name)
    if !isnothing(var)
        logical_proposition && @assert is_binary(var) "Only binary variables are allowed in $expr."
        return Symbol(name)
    end
    if expr isa Expr
        for i in eachindex(expr.args)
            expr.args[i] = replace_Symvars!(expr.args[i], model)
        end
    end

    return expr
end

function replace_JuMPvars!(expr, model)
    #replace symbolic variables and any matching expressions with JuMP variables
    var = variable_by_name(model, string(expr))
    !isnothing(var) && return var
    if expr isa Expr #run recursion
        for i in eachindex(expr.args)
            expr.args[i] = replace_JuMPvars!(expr.args[i], model)
        end
    end

    return expr
end

function replace_operators!(expr)
    #replace operators with their symbol. NOTE: Is this still needed for the hull_reformulation! of nl constraints? (check this)
    if expr isa Expr #run recursion
        for i in eachindex(expr.args)
            expr.args[i] = replace_operators!(expr.args[i])
        end
    elseif expr isa Function #replace Function with its symbol
        return Symbol(expr)
    end

    return expr
end

function replace_intevals!(expr, intervals)
    #replace variables with their intervals
    expr_str = replace(string(expr), ", " => ",") #remove any blank space after index commas
    if expr_str in keys(intervals) #check if expression is one of the model variables in the intervals dict
        return intervals[expr_str] #replace expression with interval
    elseif expr isa Expr
        if length(expr.args) == 1 #run recursive relation on the leaf node on expression tree
            expr.args[i] = replace_intevals!(expr.args[i], intervals)
        else #run recursive relation on each internal node of the expression tree, but skip the first element, which will always be the operator (this will avoid issues if the user creates a model variable called exp)
            for i in 2:length(expr.args)
                expr.args[i] = replace_intevals!(expr.args[i], intervals)
            end
        end
    end

    return expr
end

function symbolic_variable(var_ref)
    var_sym = Symbol(var_ref)
    return eval(:(Symbolics.@variables($var_sym)[1]))
end

function name_disaggregated_variable(var_ref, bin_var, i)
    #get disaggregated variable reference
    if occursin("[", string(var_ref))
        var_name_i = replace(string(var_ref), "[" => "_$bin_var$i[")
    else
        var_name_i = "$(var_ref)_$bin_var$i"
    end

    return var_name_i
end

function name_split_constraint(con_name, side)
    #get disaggregated variable reference
    if occursin("[", string(con_name))
        con_name = replace(string(con_name), "]" => ",$side]")
    else
        con_name = "$(con_name)[$side]"
    end

    return con_name
end


text:
using SparseArrays

n = 2 .^ [6 6 7 8 9 10 11 12 13 14 15 16 17 18]

K = 8

MaxGB = 2
MaxGF = 4 * 1.6

m = K .* n
sparse_gbytes = 9 * 8 * m ./ 1e9
sparse_flops = zeros(1,length(n))
sparse_gflops = zeros(1,length(n))
sparse_time = zeros(1,length(n))

for i = 1:length(n)

    global sparse_gflops, sparse_time

    ii = round.(Int,floor.(rand(m[i]) .* n[i]) .+ 1)
    jj = round.(Int,floor.(rand(m[i]) .* n[i]) .+ 1)
    A  = sparse(ii , jj , 1.0, n[i], n[i])
    ii = round.(Int,floor.(rand(m[i]) .* n[i]) .+ 1)
    jj = round.(Int,floor.(rand(m[i]) .* n[i]) .+ 1)
    B  = sparse(ii,jj,1.0,n[i],n[i])

    sparse_time[i] = @elapsed C = A*B

    sparse_flops[i] = 2*sum(C)
    sparse_gflops[i] = sparse_flops[i] ./ sparse_time[i] ./ 1e9

    print("Time: ", string(sparse_time[i]))
    print(", GFlops: ", string(sparse_gflops[i]))
    println(", GBytes: ", string(sparse_gbytes[i]))

end


text:
using IntervalArithmetic
using Base.Filesystem
using Glob

output_file_path = joinpath(dirname(dirname(pathof(IntervalArithmetic))), "test", "ITF1788_tests")
input_file_path = joinpath("output", "julia", "Base.Test", "ValidatedNumerics")
files = glob("*.jl", input_file_path)

for file in files
    bname = basename(file)
    filepath = joinpath(input_file_path, bname)
    ofilepath = joinpath(output_file_path, bname)
    cp(filepath, ofilepath; force=true)
end


text:
function init_setup()
    param = :complexConjug
    θ₀ = [10.0, -5.0, 5.0, 0.0, 3.0]
    P˟ = BSI.FitzhughDiffusion(param, θ₀...)
    obs = [[0.0, 0.0], [-0.3, -0.2], [0.4, 0.2]]
    tt = [0.0, 1.0, 1.5]
    P̃ = [BSI.FitzhughDiffusionAux(param, θ₀..., tt[1], obs[1], tt[2], obs[2]),
         BSI.FitzhughDiffusionAux(param, θ₀..., tt[2], obs[2], tt[3], obs[3])]
    setup = BSI.MCMCSetup(P˟, P̃, BSI.PartObs())
    (setup = setup, θ = θ₀, trgt = P˟, obs = obs, tt = tt, aux = P̃)
end

@testset "acceptance tracker" begin
    setup = init_setup().setup
    updt_coord = ((1,2),(2,3))
    BSI.set_transition_kernels!(setup, nothing, nothing, true, updt_coord,
                            (BSI.MetropolisHastingsUpdt(),BSI.MetropolisHastingsUpdt()))
    at = BSI.AccptTracker(setup)

    @testset "initialisation" begin
        @test at.accpt_imp == 0
        @test at.prop_imp == 0
        @test at.accpt_updt == [0, 0]
        @test at.prop_updt == [0, 0]
        @test at.updt_len == length(updt_coord)
    end

    accept_reject = [true, false, false, true, true, false]

    for ar in accept_reject
        BSI.update!(at, BSI.ParamUpdate(), 1, ar)
        BSI.update!(at, BSI.Imputation(), ar)
    end

    @testset "testing update! (1/2)" begin
        accept_rate = sum(accept_reject)/length(accept_reject)
        @test BSI.accpt_rate(at, BSI.ParamUpdate())[1] == accept_rate
        @test isnan(BSI.accpt_rate(at, BSI.ParamUpdate())[2])
        @test BSI.accpt_rate(at, BSI.Imputation()) == accept_rate
    end

    for ar in accept_reject BSI.update!(at, BSI.ParamUpdate(), 2, ar) end

    @testset "testing update! (2/2)" begin
        accept_rate = sum(accept_reject)/length(accept_reject)
        @test BSI.accpt_rate(at, BSI.ParamUpdate())[1] == accept_rate
        @test BSI.accpt_rate(at, BSI.ParamUpdate())[2] == accept_rate
        @test BSI.accpt_rate(at, BSI.Imputation()) == accept_rate
    end

end

@testset "parameter history" begin
    out = init_setup()
    setup, θ = out.setup, out.θ
    updt_coord = ((1,2),(2,3))
    BSI.set_transition_kernels!(setup, nothing, nothing, true, updt_coord,
                            (BSI.MetropolisHastingsUpdt(),BSI.MetropolisHastingsUpdt()))

    num_mcmc_steps = 1000
    warm_up = 50
    BSI.set_mcmc_params!(setup, num_mcmc_steps, nothing, nothing, nothing, warm_up)
    ph = BSI.ParamHistory(setup)

    _foo(x::BSI.ParamHistory{T}) where T = T

    @testset "initialisation" begin
        @test eltype(ph.θ_chain) == typeof(θ) == _foo(ph)
        @test length(ph.θ_chain) == length(updt_coord)*(num_mcmc_steps-warm_up)+1
        @test ph.θ_chain[1] == last(ph) == θ
        @test ph.counter == 1
    end

    N = 10
    chain = [rand(5) for i in 1:N]
    for i in 1:N BSI.update!(ph, chain[i]) end

    @testset "updating chain" begin
        @test ph.θ_chain[1] == θ
        @test all([ph.θ_chain[i+1] == chain[i] for i in 1:N])
        @test BSI.last(ph) == ph.θ_chain[N+1] == chain[N]
        @test BSI.last(ph) != θ
    end
end

@testset "action tracker" begin
    setup = init_setup().setup
    save_iter = 10
    verb_iter = 3
    warm_up = 50
    BSI.set_mcmc_params!(setup, nothing, save_iter, verb_iter, nothing, warm_up)
    param_updt = true
    BSI.set_transition_kernels!(setup, nothing, nothing, param_updt)
    at = BSI.ActionTracker(setup)

    @testset "initialisation" begin
        @test at.save_iter == save_iter
        @test at.verb_iter == verb_iter
        @test at.warm_up == warm_up
        @test at.param_updt == param_updt
    end

    @testset "correct acting" begin
        @test !BSI.act(BSI.SavePath(), at, 1)
        @test !BSI.act(BSI.SavePath(), at, 10)
        @test BSI.act(BSI.SavePath(), at, 60)
        @test !BSI.act(BSI.SavePath(), at, 51)
        @test !BSI.act(BSI.SavePath(), at, 61)
        @test BSI.act(BSI.SavePath(), at, 40000)

        @test !BSI.act(BSI.Verbose(), at, 1)
        @test BSI.act(BSI.Verbose(), at, 3)
        @test !BSI.act(BSI.Verbose(), at, 10)
        @test BSI.act(BSI.Verbose(), at, 30)

        @test !BSI.act(BSI.ParamUpdate(), at, 1)
        @test !BSI.act(BSI.ParamUpdate(), at, 10)
        @test !BSI.act(BSI.ParamUpdate(), at, 50)
        @test BSI.act(BSI.ParamUpdate(), at, 51)
        @test BSI.act(BSI.ParamUpdate(), at, 1000)
    end
end

@testset "workspace" begin
    out = init_setup()
    setup, θ, trgt, obs, tt = out.setup, out.θ, out.trgt, out.obs, out.tt
    aux = out.aux

    L = [1. 0.; 0. 1.]
    Σ = [0.5 0.0; 0.0 1.0]
    BSI.set_observations!(setup, [L, L], [Σ, 2*Σ], obs, tt)

    dt = 0.01
    τ(t₀,T) = (x) ->  t₀ + (x-t₀) * (2-(x-t₀)/(T-t₀))
    BSI.set_imputation_grid!(setup, dt, τ)

    t_kernels = [BSI.RandomWalk([0.002, 0.1], [true, true]),
                 BSI.RandomWalk([0.2, 1.0], [false, true])]
    ρ = 0.5
    param_updt = true
    updt_coord = (Val((true,true,false)),
                  Val((false,true,true)))
    updt_type=(BSI.MetropolisHastingsUpdt(),
               BSI.ConjugateUpdt())
    BSI.set_transition_kernels!(setup, t_kernels, ρ, param_updt, updt_coord,
                            updt_type)

    priors = BSI.Priors((BSI.ImproperPrior(), BSI.ImproperPrior()))
    x0_prior = BSI.KnownStartingPt(obs[1])
    BSI.set_priors!(setup, priors, x0_prior)

    num_mcmc_steps = 100
    BSI.set_mcmc_params!(setup, num_mcmc_steps)

    BSI.set_solver!(setup)

    BSI.initialise!(Float64, setup)

    ws = BSI.Workspace(setup).workspace
    @testset "initialisation" begin
        @test typeof(ws.Wnr) == Wiener{Float64}
        @test length(ws.XX) == length(ws.XXᵒ) == 2
        @test length(ws.WWᵒ) == length(ws.WW) == 2
        @test all([ws.XX[i].tt == ws.XXᵒ[i].tt for i in 1:2])
        @test all([ws.XX[i].yy == ws.XXᵒ[i].yy for i in 1:2])
        @test all([ws.WW[i].tt == ws.WWᵒ[i].tt for i in 1:2])
        @test all([ws.WW[i].yy == ws.WWᵒ[i].yy for i in 1:2])
        @test eltype(eltype(ws.WW)) == Float64
        @test eltype(eltype(ws.XX)) == SArray{Tuple{2},Float64,1,2}
        @test ws.WW[1].tt[1] == ws.XX[1].tt[1] == tt[1]
        @test ws.WW[2].tt[1] == ws.XX[2].tt[1] == tt[2]
        @test typeof(ws.P[1].Target) <: BSI.FitzhughDiffusion
        @test typeof(ws.P[1].Pt) <: BSI.FitzhughDiffusionAux
        @test all([ws.P[i].tt == ws.WW[i].tt == ws.XX[i].tt for i in 1:2])
        @test ws.recompute_ODEs == [true, true]
        @test ws.ρ == [[0.5]]
    end

    # NOTE this constructor has been deprecated
    #ws2 = Workspace(ws, 0.25)
    #@testset "copy constructor" begin
    #    @test ws2.ρ == 0.25
    #    @test ws.ρ == 0.5
    #    @test ws.Wnr == ws2.Wnr
    #    @test ws.XX == ws2.XX
    #    @test ws.WW == ws2.WW
    #    @test ws.P == ws.P
    #end

    @testset "action determination" begin
        @test !BSI.act(BSI.SavePath(), ws, 1)
        @test !BSI.act(BSI.SavePath(), ws, 10)
        @test !BSI.act(BSI.SavePath(), ws, 60)
        @test !BSI.act(BSI.SavePath(), ws, 51)
        @test !BSI.act(BSI.SavePath(), ws, 61)
        @test !BSI.act(BSI.SavePath(), ws, 40000)

        @test !BSI.act(BSI.Verbose(), ws, 1)
        @test !BSI.act(BSI.Verbose(), ws, 3)
        @test !BSI.act(BSI.Verbose(), ws, 10)
        @test !BSI.act(BSI.Verbose(), ws, 30)

        @test BSI.act(BSI.ParamUpdate(), ws, 1)
        @test BSI.act(BSI.ParamUpdate(), ws, 10)
        @test BSI.act(BSI.ParamUpdate(), ws, 50)
        @test BSI.act(BSI.ParamUpdate(), ws, 51)
        @test BSI.act(BSI.ParamUpdate(), ws, 1000)
    end
end

@testset "definition of gibbs sweep" begin
    out = init_setup()
    setup, obs, tt = out.setup, out.obs, out.tt
    L = [1. 0.; 0. 1.]
    Σ = [0.5 0.0; 0.0 1.0]
    BSI.set_observations!(setup, [L, L], [Σ, 2*Σ], obs, tt)

    dt = 0.01
    τ(t₀,T) = (x) ->  t₀ + (x-t₀) * (2-(x-t₀)/(T-t₀))
    BSI.set_imputation_grid!(setup, dt, τ)

    t_kernels = [BSI.RandomWalk([0.002, 0.1], [true, true]),
                 BSI.RandomWalk([0.2, 1.0], [false, true])]
    ρ = 0.5
    param_updt = true
    updt_coord = (Val((true,true,false)),
                  Val((false,true,true)))
    updt_type=(BSI.MetropolisHastingsUpdt(),
               BSI.ConjugateUpdt())
    BSI.set_transition_kernels!(setup, t_kernels, ρ, param_updt, updt_coord,
                            updt_type)

    priors = BSI.Priors((BSI.ImproperPrior(), BSI.ImproperPrior()))
    x0_prior = BSI.KnownStartingPt(obs[1])
    BSI.set_priors!(setup, priors, x0_prior)

    num_mcmc_steps = 100
    BSI.set_mcmc_params!(setup, num_mcmc_steps)

    BSI.set_solver!(setup)

    BSI.initialise!(Float64, setup)
    gs = BSI.GibbsDefn(setup)
    @testset "initialisation" begin
        @test length(gs) == 2
        @test all([gs[i] == gs.updates[i] for i in 1:2])
        @test all([gs[i].updt_type == updt_type[i] for i in 1:2])
        @test all([gs[i].updt_coord == updt_coord[i] for i in 1:2])
        @test all([gs[i].priors == priors[i] for i in 1:2])
        @test gs[1].recompute_ODEs && gs[2].recompute_ODEs
    end
end


text:
import LatticeBoltzmann: LatticeBoltzmannModel, ProcessingMethod

@testset "TRT with τ_a = τ_s is identical to SRT" begin
    τ = 0.8

    q = D2Q9()
    N = 10
    f_in = zeros(N, N, length(q.weights))
    nx, ny, nf = size(f_in)
    @inbounds for x in 1:nx, y in 1:ny
        @inbounds for f_idx in 1:nf
            f_in[x, y, f_idx] = q.weights[f_idx]
        end
    end
    f_srt = similar(f_in)
    srt = SRT(τ)
    collide!(srt, q, f_in, f_srt)

    f_trt = similar(f_in)
    trt = LatticeBoltzmann.TRT(τ, τ)
    collide!(trt, q, f_in, f_trt)

    @test f_srt ≈ f_trt
    @testset "Relaxation time $τ" for τ in [0.51, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1]
        cs² = q.speed_of_sound_squared
        ν = (τ - 0.5) / cs²

        problem = PoiseuilleFlow(ν, 1, static = false)

        process_method = ProcessingMethod(problem, false, 10)
        force_field = LatticeBoltzmann.has_external_force(problem) ?
            (x_idx, y_idx, t) -> lattice_force(problem, x_idx, y_idx, t) : nothing

        srt_model = LatticeBoltzmannModel(
            problem,
            q,
            collision_model = SRT(τ),
            process_method = process_method,
        )
        trt_model = LatticeBoltzmannModel(
            problem,
            q,
            collision_model = LatticeBoltzmann.TRT(τ, τ),
            process_method = process_method,
        )

        srt_result = simulate(srt_model, 0:10)
        trt_result = simulate(trt_model, 0:10)

        nx, ny, nf = size(srt_result.f_stream)
        @inbounds for x in 1:nx, y in 1:ny, f_idx in 1:nf
            @test srt_result.f_stream[x, y, f_idx] .≈ trt_result.f_stream[x, y, f_idx]
        end
    end
end

@testset "MRT with τ_i = τ is identical to SRT" begin
    τ = 0.8

    q = D2Q9()
    N = 4
    f_in = zeros(N, N, length(q.weights))
    nx, ny, nf = size(f_in)
    @inbounds for x in 1:nx, y in 1:ny
        @inbounds for f_idx in 1:nf
            f_in[x, y, f_idx] = q.weights[f_idx]
        end
    end
    f_srt = similar(f_in)
    srt = SRT(τ)
    collide!(srt, q, f_in, f_srt)

    f_mrt = similar(f_in)
    mrt = LatticeBoltzmann.MRT(q, τ)
    collide!(mrt, q, f_in, f_mrt)

    @test f_srt ≈ f_mrt
    @testset "Relaxation time $τ" for τ in [0.51, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1]
        cs² = q.speed_of_sound_squared
        ν = (τ - 0.5) / cs²

        problem = PoiseuilleFlow(ν, 1, static = false)

        process_method = ProcessingMethod(problem, false, 10)
        force_field = LatticeBoltzmann.has_external_force(problem) ?
            (x_idx, y_idx, t) -> lattice_force(problem, x_idx, y_idx, t) : nothing

        srt_model = LatticeBoltzmannModel(
            problem,
            q,
            collision_model = SRT(τ),
            process_method = process_method,
        )
        mrt_model = LatticeBoltzmannModel(
            problem,
            q,
            collision_model = LatticeBoltzmann.MRT(q, τ),
            process_method = process_method,
        )

        srt_result = simulate(srt_model, 0:10)
        mrt_result = simulate(mrt_model, 0:10)

        nx, ny, nf = size(srt_result.f_stream)
        @inbounds for x in 1:nx, y in 1:ny, f_idx in 1:nf
            if τ == 1.0
                @test srt_result.f_stream[x, y, f_idx] .≈ mrt_result.f_stream[x, y, f_idx]
            else
                @test_broken srt_result.f_stream[x, y, f_idx] .≈
                             mrt_result.f_stream[x, y, f_idx]
            end
        end
    end
end


text:
using Gadfly, DataArrays, RDatasets

plot(dataset("datasets", "iris"),
     layer(x=:SepalLength, y=:SepalWidth, Geom.point),
     layer([sin, cos], 0, 25))


text:
module ExampleEffect

using ..Ahorn, Maple

@mapdef Effect "ExampleMod/ExampleEffect" Example(only::String="*", exclude::String="")
#= evaluates to:
    const Example = Effect{Symbol("YourMod/YourEffect")}
    Example(only::String = "*", exclude::String = "") =
        Effect{Symbol("YourMod/YourEffect")}("YourMod/YourEffect", only = only, exclude = exclude)
=#

@mapdef Effect "ExampleMod/ExampleEffect2" Example2(only::String="*", exclude::String="", tag::String="", flag::String="", notFlag::String="")

const placements = Example
# `placements` can also be an array to register multiple effects:
#const placements = [Example, Example2]

Ahorn.canFgBg(effect::Example) = true, true

# The order in which attributes will be displayed when editing.
# Attributes that aren't listed are appended in alphabetical order.
editingOptions(entity::Maple.Effect) = Dict{String, Any}()

# Attributes that should not be displayed when editing.
# `multiple` is true when more than one entity is being edited simultaneously.
editingIgnored(entity::Maple.Effect) = String[]

# Attributes that should not be displayed when editing.
# `multiple` is true when more than one entity is being edited simultaneously.
editingOrder(entity::Maple.Effect) = String[
    "name", "only", "exclude", "tag",
    "flag", "notflag"
]

end


text:
using MLDataPattern
using DataFrames

function crossvalidate(data,target,createmodel,loss,bs=100,n=10000)
	function onerun(xTrn,xVal)
		m = createmodel()
		opt = Flux.Optimise.ADAM(params(m))
		Flux.train!((xx...) -> loss(m(getobs(xx[1])),getobs(xx[2])), RandomBatches(xTrn,bs,n), opt)
		(mean(Flux.argmax(m(getobs(xTrn[1]))) .!= Flux.argmax(getobs(xTrn[2]))),
		mean(Flux.argmax(m(getobs(xVal[1]))) .!= Flux.argmax(getobs(xVal[2]))))
	end
	map(x -> onerun(x[1],x[2]), kfolds(shuffleobs((data,target)),5))
end

function gridsearch(cr,parameters)
	errs = map(p -> (p,cr(p)),parameters)
	i = indmin(map(v -> mean(v[2]),errs))
	return(errs[i][1],errs[2][2],errs)
end

function gridsearch(cr,parameters,names::Vector{Symbol},ofname::String = "")
	df = []
	for p in parameters
		dff = DataFrame([cr(p)],[:errs])
		dff[:fold] = 1:size(dff,1)
		foreach(pn -> dff[pn[1]] = pn[2],zip(names,p))
		push!(df,dff)
		!isempty(ofname) && CSV.write(ofname, dff, append = true, header = true)
		display(dff)
	end
	vcat(df...)
end
gridsearch(cr,parameters,names::Vector{String},ofname::String = "") = gridsearch(cr,parameters,Symbol.(names),ofname)

max_stars_repo_name: pevnak/FluxExtensions

max_issues_repo_name: pevnak/FluxExtensions

max_forks_repo_name: pevnak/FluxExtensions


text:
module LAModelData

    using DataFrames
    using CSV
    using DataFramesMeta
    using Query
    using Logging

    using GlobalDecls
    using Utils
    using ONSCodes

    export is_aggregate, get_pay_class_from_region_code
    export get_staying_put_rates, get_target_la, do_exit_rate_query
    export UDATA, get_18s_level_from_doe, load_one

    function load_one( name :: AbstractString ) :: DataFrame
        df = CSV.File(
            name,
            delim=',',
            missingstrings=["x","","-", "c",".", "..", "*"],
            types=make_type_block(6:1000)) |> DataFrame
        lcnames = Symbol.(lowercase.(string.(names(df))))
        rename!(df,lcnames)
        df
    end

    #
    # DFE Raw data on care movements 2015-2018; see
    # only 2017 and 2018, 2019 have Staying Put data
    #
    function get_underlying_data()::NamedTuple
        (
        c17_18_2017=load_one(
                DATADIR*"underlying_data/2017/SFR50_CareLeavers17182017.csv" ),
        c19_21_2017 = load_one(
                DATADIR*"underlying_data/2017/SFR50_CareLeavers19to212017.csv" ),
        c17_18_2018 = load_one(
                DATADIR*"underlying_data/2018/CareLeavers17182018_amended.csv" ),
        c19_21_2018 = load_one(
                DATADIR*"underlying_data/2018/CareLeavers_Acc_StayPut19to212018.csv" ),
        c17_18_2019 = load_one(
                DATADIR*"underlying_data/2019/Care_Leavers_17182019.csv" ),
        c19_21_2019 = load_one(
                DATADIR*"underlying_data/2019/CareLeavers_Acc_StayPut19to212019.csv" )
        )
    end

    UDATA = get_underlying_data()

    """
    get the number of 18s leaving care into staying put from the DoE underlying data as an
    alternative to a calculation from OFSTED data and ageing rates
    """
    function get_18s_level_from_doe( target :: AbstractString )
        which = UDATA.c17_18_2019.new_geog_code .== target
        #cl_stayput_18 Care leavers in the year ending 31 March [201X] aged 18 who ceased to be looked after from a foster placement on their 18th birthday
        v = UDATA.c17_18_2019[which,:cl_stayput_18][1]
        #if( ismissing( v ))
    #        v=-1
        #end
        return v

        q = @from i in UDATA.c17_18_2019 begin
            @where i.new_geog_code == target
            @select i.cl_stayput_18
            @collect
        end
        q[1].value
    end

    """
    extract the number leaving care and the number staying put
    for the given target (LA, country, etc.)
    age (18,19,20) and year (2017,2018,2019)
    """
    function do_exit_rate_query(
        year   :: Integer,
        target :: AbstractString,
        age    :: Integer ) :: AbstractArray
        as = age == 18 ? "17_18" : "19_21"
        # something like :c17_18_2019
        field = Symbol( "c"*as*"_"*"$year" )
        f1 = Symbol( "cl_stayput_$age" )
        f2 = Symbol( "cl_stayput_ffc_$age" )
        @info  "target $target age $age f1 $f1 f2 $f2 field $field"
        q = @from i in UDATA[field] begin
            @where i.new_geog_code == target
            @select i[f1], i[f2]
            @collect
        end
        # @assert length( q ) == 1 no - 3 londons ffs we use London not Inner/Outer London
        @info  "q1 $q[1]"
        q[1]
    end

    """
     return ONS code of either the LA, its region, or England, depending on `AggLevel`
    """
    function get_target_la(
        lacode :: AbstractString,
        agglev :: AggLevel ) :: AbstractString
        target = lacode
        if agglev == regional
            target = region_code_from_code( lacode )
        elseif agglev == national
            target = "E92000001" # england
        end
        target
    end

    """
    return the exit rates for the given LA - a named tuple a18,a19,a20
      will be the regional one if the calculation is not possible for an LA
      rates are proportions for that age - so 0.5,0.5,0.5 means 50% stay 1 year, 25% 2 12.5% 3 and so on
    """
    function get_staying_put_rates(
        lacode    :: AbstractString,
        agglev    :: AggLevel,
        poolyears :: Bool  ) :: Vector
        # poolyears is false by defailt
        # @info  "get_exit_rates for lacode $lacode agglev $agglev" ;
        avprop = zeros(3)
        target = get_target_la( lacode, agglev ) # just swaps region or national code if the agglevel<>council
        ages = [18,19,20]
        years = poolyears ? [2017,2018,2019] : [GlobalDecls.SIMULATION_YEAR]
        for year in years
            i = 0
            numreachingage=zeros(3)
            numstaying=zeros(3)
            prop=zeros(3)
            for age in ages
                i+=1
                q = do_exit_rate_query( year, target, age )
                if zeroormissing( q )
                    @assert agglev == local_authority
                    target = get_target_la( lacode, regional )
                    @info  "zeros detected! lacode $lacode new target $target $year age $age " ;
                    q = do_exit_rate_query( year, target, age )
                end
                @assert (! zeroormissing( q ))
                numreachingage[i] += get(q[1],0)
                numstaying[i] += get(q[2],0)
            end
            prop = zeros(3)
            #
            # FIXME this goes wrong if 1 year is region and the
            # other is LA
            #
            @info  "got exit data numstaying = $numstaying numreachingage $numreachingage"

            # prop[1] = numstaying[1]/numreachingage[1]
            # prop[2] = numstaying[2]/numstaying[1]
            # prop[3] = numstaying[3]/numstaying[2]
            # correct way, I think - prop is scale
            prop = numstaying./numreachingage
            prop[2] /= prop[1]
            prop[3] /= (prop[2]*prop[1])
            # @info  "prop $prop " ;
            # @info  "numstaying $numstaying numreachingage $numreachingage"
            avprop .+= prop
        end
        avprop /= length( years )
        # t = (a18=avprop[1],a19=avprop[2],a20=avprop[3])
        return avprop
    end

    # df meta example:
    # @linq las |> where( :LAD19NM.=="City of London" ) |> select( :LAD19CD )

    function get_pay_class_from_region_code( regcode :: AbstractString ) :: Region
        # E12000001  #	A 	North East England
        # E12000002  # 	B 	North West England
        # E12000003  #	D 	Yorkshire and the Humber
        # E12000004  #	E 	East Midlands
        # E12000005  # 	F 	West Midlands
        # E12000006  # 	G 	East of England
        if regcode == "E12000007"
            r = london # 	H 	London
        elseif regcode == "E12000008"
            r = se
            # J 	South East England
        else
            r = rest_of_england
        end
        # E12000009 	K 	South West England
        r
    end

    function is_aggregate( name )
        return name in [
            "East of England",
            "England",
            "North East England",
            "North West England",
            "Yorkshire and The Humber",
            "East Midlands",
            "West Midlands",
            "London",
            "Inner London",
            "Outer London",
            "South East England",
            "South West England",
            "London Tri-borough" ]
    end

end # LAModelData


text:
module CryptoSignatures

using CryptoGroups
using Random

rng() = RandomDevice()

# const _default_rng = Ref{RandomDevice}()
# function __init__()
#     _default_rng[] = RandomDevice()
# end

# default_rng() = _default_rng[]

abstract type AbstractSignature end

id(s::AbstractSignature) = s.pubkey ### One can overwrite this as one wishes

"""
Never use this for deployed application!!!
"""
verify(data,hashnum::UInt64) = hash(data)==hashnum
verify(data,s::AbstractSignature,G::AbstractGroup) = verify(s,G) && verify(data,s.hash)

abstract type AbstractSigner end
id(s::AbstractSigner) = s.pubkey

abstract type AbstractEncryptionKey end
abstract type AbstractDecryptionKey end

# The singer does store group to prevent the private key to be accidentally missused with different group which would lower the security. And the singature would be useless in such case.
struct Signer{K,P} <: AbstractSigner
    privkey::K
    pubkey::P
    G::AbstractGroup
end

function Signer(G::AbstractGroup;rng=rng())
    x = rand(1:order(G))
    y = binary(G^x)
    Signer(x,y,G)
end

### There are many different ways one could sing stuff.
# include("rsasignatures.jl")
include("dsasignatures.jl")

export verify, rsasign, id, Signer, DSASignature

end # module