text:
## ExoplanetsSysSim/src/transit_detection_model.jl
## (c) 2015 Eric B. Ford

# Several functions below based on https://github.com/christopherburke/KeplerPORTs/blob/master/KeplerPORTs_utils.py
# That follows the procedure outlined in Burke et al.(2015).
# However we don't currently interpolate the CDDP or mesthreshold to the relevant duration
import SpecialFunctions.lgamma

function real_log_choose(m::Float64, n::Float64)::Float64
  lgamma(m+1)-lgamma(n+1)-lgamma(m-n+1.0)
end

function real_binom(k::Float64, BigM::Float64, f::Float64)::Float64
    F1 = real_log_choose(BigM,k)
    F2 = k*log(f)
    F3 = (BigM-k)*log(1.0-f)
    x = exp(F1+F2+F3)
    return x
end

function kepler_window_function_binomial_model(exp_num_transits_no_gaps::Float64, duty_cycle::Float64; min_transits::Float64 = 3.0)::Float64
  if exp_num_transits_no_gaps < min_transits
     return 0.0
  else
     return max(1.0 - real_binom(min_transits,exp_num_transits_no_gaps,duty_cycle), 0.0)
  end
end

function kepler_window_function_binomial_model(t::KeplerTarget, exp_num_transits_no_gaps::Float64, period::Float64, duration::Float64; min_transits::Float64 = 3.0)::Float64
   kepler_window_function_binomial_model(exp_num_transits_no_gaps, t.duty_cycle, min_transits=min_transits)
end

function kepler_window_function_dr25_model(t::KeplerTarget, exp_num_transits_no_gaps::Float64, period::Float64, duration::Float64)::Float64
   ExoplanetsSysSim.WindowFunction.eval_window_function(t.window_function_id, Duration=duration, Period=period)
end

#kepler_window_function = kepler_window_function_binomial_model
kepler_window_function = kepler_window_function_dr25_model

# See Christiansen et al. (2015)  This assumes a linear limbdarkening coefficient of 0.6
function frac_depth_to_tps_depth(frac_depth::Float64)
    alp = 1.0874
    bet = 1.0187
    REALDEPTH2TPSSQUARE = 1.0  # WARNING: Waiting for this to be confirmed
    k = sqrt(frac_depth)
    tps_depth = min( (alp-bet*k) * frac_depth* REALDEPTH2TPSSQUARE, 1.0)   # NOTE: I added the max based on common sense
    return tps_depth::Float64
end

function detection_efficiency_theory(mes::Float64; min_pdet_nonzero::Float64 = 0.0)
   muoffset =  0.0
   sig =  1.0
   mesthresh = 7.1
   if mes > (9.0 - mesthresh - muoffset)
      return 1.0
   else
      pdet = 0.5 + 0.5*erf((mes - mesthresh - muoffset) / sqrt(2.0*sig*sig))
      pdet = pdet >= min_pdet_nonzero ? pdet : 0.0
      return pdet
   end
end

function detection_efficiency_fressin2013(mes::Float64)
    mesmin =  6.0
    mesmax =  16.0
   if mes <= mesmin
      return 0.0
   elseif mes >= mesmax
      return 1.0
   else
     return (mes - mesmin) / (mesmax - mesmin)
  end
end

function detection_efficiency_christiansen2015(mes::Float64; mes_threshold::Float64 = 7.1, min_pdet_nonzero::Float64 = 0.0)
    a =  4.65  # from code for detection_efficiency(...) at https://github.com/christopherburke/KeplerPORTs/blob/master/KeplerPORTs_utils.py
    b =  0.98
   # a =  4.35 # from arxiv abstract. Informal testing showed it didn't matter
   # b =  1.05
   usemes = max(0.0,mes - 7.1 - (mes_threshold - 7.1))
   pdet = cdf(Gamma(a,b), usemes)
   pdet = pdet >= min_pdet_nonzero ? pdet : 0.0
   return pdet
end

function detection_efficiency_dr25_simple(mes::Float64; min_pdet_nonzero::Float64 = 0.0)::Float64
   a = 30.87  # from pg 16 of https://exoplanetarchive.ipac.caltech.edu/docs/KSCI-19110-001.pdf
   b = 0.271
   c = 0.940
   dist = Gamma(a,b)
   pdet::Float64 = c*cdf(dist, mes)::Float64
   pdet = pdet >= min_pdet_nonzero ? pdet : 0.0
   return pdet
end

# From Mulders et al. 2018 (arXiv 1805.08211)
# WARNING: Assumes inputs are in R_sol and days
function vetting_efficiency_dr25_mulders(R_p::Real, P::Real)
    c = 0.63
    a_R = 0.19
    P_break = 53.
    a_P = -0.07
    b_P = -0.39
    pvet = c*(R_p/earth_radius)^a_R
    if P < P_break
        pvet *= (P/P_break)^a_P
    else
        pvet *= (P/P_break)^b_P
    end
    return pvet
end

# WARNING: Hardcoded choice of transit detection efficiency here for speed and so as to not have it hardcoded in multiple places
#detection_efficiency_model = detection_efficiency_christiansen2015
detection_efficiency_model = detection_efficiency_dr25_simple

# Resume code original to SysSim

function interpolate_cdpp_to_duration(t::KeplerTarget, duration::Float64)::Float64
   duration_in_hours = duration *24.0
   dur_idx = searchsortedlast(cdpp_durations,duration_in_hours)   # cdpp_durations is defined in constants.jl
   if dur_idx <= 0
      cdpp = t.cdpp[1,1]
   elseif dur_idx==length(cdpp_durations) && (duration_in_hours > cdpp_durations[end]) # Should be 15 cdpp_durations.
      cdpp = t.cdpp[length(cdpp_durations),1]
   else
      w = ((duration_in_hours)-cdpp_durations[dur_idx]) / (cdpp_durations[dur_idx+1]-cdpp_durations[dur_idx])
      cdpp = w*t.cdpp[dur_idx+1,1] + (1-w)*t.cdpp[dur_idx,1]
   end
   return cdpp
end

function calc_snr_if_transit(t::KeplerTarget, depth::Float64, duration::Float64, cdpp::Float64, sim_param::SimParam; num_transit::Float64 = 1)
   depth_tps = frac_depth_to_tps_depth(depth)                  # TODO SCI:  WARNING: Hardcoded this conversion.  Remove once depth calculated using limb darkening model
   snr = depth_tps*sqrt(num_transit*duration*LC_rate)/cdpp     # WARNING: Assumes measurement uncertainties are uncorrelated & CDPP based on LC
end

function calc_snr_if_transit_central(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  depth = calc_transit_depth(t,s,p)
  duration_central = calc_transit_duration_central(t,s,p)
  cdpp = interpolate_cdpp_to_duration(t, duration_central)
  num_transit = calc_expected_num_transits(t,s,p,sim_param)
  calc_snr_if_transit(t,depth,duration_central,cdpp, sim_param,num_transit=num_transit)
end

function calc_prob_detect_if_transit(t::KeplerTarget, snr::Float64, period::Float64, duration::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  min_pdet_nonzero = 1.0e-4                                                # TODO OPT: Consider raising threshold to prevent a plethora of planets that are very unlikely to be detected due to using 0.0 or other small value here
  wf = kepler_window_function(t, num_transit, period, duration)
  return wf*detection_efficiency_model(snr, min_pdet_nonzero=min_pdet_nonzero)
end

function calc_prob_detect_if_transit(t::KeplerTarget, depth::Float64, period::Float64, duration::Float64, cdpp::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  snr = calc_snr_if_transit(t,depth,duration,cdpp, sim_param, num_transit=num_transit)
  return calc_prob_detect_if_transit(t, snr, period, duration, sim_param, num_transit=num_transit)
end

function calc_prob_detect_if_transit_central(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  period = t.sys[s].orbit[p].P
  depth = calc_transit_depth(t,s,p)
  duration_central = calc_transit_duration_central(t,s,p)
  cdpp = interpolate_cdpp_to_duration(t, duration_central)
  ntr = calc_expected_num_transits(t,s,p,sim_param)
  calc_prob_detect_if_transit(t,depth,period,duration_central,cdpp, sim_param, num_transit=ntr)
end

function calc_prob_detect_if_transit_with_actual_b(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  period = t.sys[s].orbit[p].P
  size_ratio = t.sys[s].planet[p].radius/t.sys[s].star.radius
  depth = calc_transit_depth(t,s,p)
  duration = calc_transit_duration(t,s,p)
  cdpp = interpolate_cdpp_to_duration(t, duration)
  ntr = calc_expected_num_transits(t,s,p,sim_param)
  calc_prob_detect_if_transit(t,depth,period,duration,cdpp, sim_param, num_transit=ntr)
end

# Compute probability of detection if we average over impact parameters b~U[0,1)
function calc_ave_prob_detect_if_transit_from_snr(t::KeplerTarget, snr_central::Float64, period::Float64, duration_central::Float64, size_ratio::Float64, cdpp_central::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  min_pdet_nonzero = 1.0e-4
  wf = kepler_window_function(t, num_transit, period, duration_central)

  detection_efficiency_central = detection_efficiency_model(snr_central, min_pdet_nonzero=min_pdet_nonzero)
  if wf*detection_efficiency_central <= min_pdet_nonzero
     return 0.
  end

  # Breaking integral into two sections [0,1-b_boundary) and [1-b_boundary,1], so need at least 5 points to evaluate integral via trapezoid rule
  num_impact_param_low_b =  7                            # Number of points to evaluate integral over [0,1-b_boundary) via trapezoid rule
  num_impact_param_high_b = 5 # (size_ratio<=0.05) ? 5 : 11  # Number of points to evaluate integral over [1-b_boudnary,1) via trapezoid rule.  If using 2*size_ratio for bondary for small planets, then keep this odd, so one point lands on 1-size_ratio.
  @assert(num_impact_param_low_b >= 5)
  @assert(num_impact_param_high_b >= 3)
  num_impact_param = num_impact_param_low_b+num_impact_param_high_b-1 # One point is shared
  b_boundary = (size_ratio <= 0.15) ? 2*size_ratio : min(max(0.3,size_ratio),0.5)
  b = Array{Float64}(undef,num_impact_param)
  weight = Array{Float64}(undef,num_impact_param)
  b[1:num_impact_param_low_b] = range(0.0,stop=1-b_boundary,length=num_impact_param_low_b)
  b[num_impact_param_low_b:num_impact_param] .= range(1-b_boundary,stop=1.0,length=num_impact_param_high_b)
  weight[1:num_impact_param_low_b] .= (1-b_boundary)/(num_impact_param_low_b-1)  # Points for first integral
  weight[1] *= 0.5                        # Lower endpoint of first integral
  weight[num_impact_param_low_b] *= 0.5   # Upper endpoint of first integral
  weight[num_impact_param_low_b] += 0.5*(b_boundary)/(num_impact_param_high_b-1) # Also lower endpoint of second integral
  weight[(num_impact_param_low_b+1):num_impact_param] .= b_boundary/(num_impact_param_high_b-1)
  weight[num_impact_param] *= 0.5         # Upper endpoint of second integral
  @assert isapprox(sum(weight),1.0)

  function integrand(b::Float64)::Float64
     depth_factor = calc_depth_correction_for_grazing_transit(b,size_ratio)
     duration_factor = calc_transit_duration_factor_for_impact_parameter_b(b,size_ratio)
     cdpp = interpolate_cdpp_to_duration(t,duration_central*duration_factor)
     snr_factor = depth_factor*sqrt(duration_factor)*(cdpp_central/cdpp)
     detection_efficiency_model(snr_central*snr_factor, min_pdet_nonzero=min_pdet_nonzero)
  end

  ave_detection_efficiency = sum(weight .* map(integrand,b)::Vector{Float64} )

  return wf*ave_detection_efficiency
end

function calc_ave_prob_detect_if_transit(t::KeplerTarget, depth::Float64, period::Float64, duration_central::Float64, size_ratio::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  cdpp_central = interpolate_cdpp_to_duration(t, duration_central)
  snr_central = calc_snr_if_transit(t,depth,duration_central,cdpp_central, sim_param, num_transit=num_transit)
  return calc_ave_prob_detect_if_transit_from_snr(t, snr_central, period, duration_central, size_ratio, cdpp_central, sim_param, num_transit=num_transit)
end

function calc_ave_prob_detect_if_transit(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  size_ratio = t.sys[s].planet[p].radius/t.sys[s].star.radius
  depth = calc_transit_depth(t,s,p)
  period = t.sys[s].orbit[p].P
  duration_central = calc_transit_duration_central(t,s,p)
  ntr = calc_expected_num_transits(t,s,p,sim_param)
  calc_ave_prob_detect_if_transit(t,depth,period,duration_central, size_ratio, sim_param, num_transit=ntr)
end


text:
using Sundials, Compat

## f routine. Compute function f(t,y).

function f(t, y_nv, ydot_nv, user_data)
    y = convert(Vector, y_nv)
    ydot = convert(Vector, ydot_nv)
    ydot[1] = -0.04*y[1] + 1.0e4*y[2]*y[3]
    ydot[3] = 3.0e7*y[2]*y[2]
    ydot[2] = -ydot[1] - ydot[3]
    return Sundials.CV_SUCCESS
end

neq = 3

t0 = 0.0
t1 = 0.4
tmult = 10.0
nout = 12
y0 = [1.0, 0.0, 0.0]
reltol = 1e-4
abstol = 1e-11
userdata = nothing
h0 = 1e-4 * reltol

mem_ptr = Sundials.ARKodeCreate()
arkode_mem = Sundials.Handle(mem_ptr)
Sundials.@checkflag Sundials.ARKodeInit(arkode_mem, C_NULL, f, t0, y0)
Sundials.@checkflag Sundials.ARKodeSetInitStep(arkode_mem, h0)
Sundials.@checkflag Sundials.ARKodeSetMaxErrTestFails(arkode_mem, 20)
Sundials.@checkflag Sundials.ARKodeSetMaxNonlinIters(arkode_mem, 8)
Sundials.@checkflag Sundials.ARKodeSetNonlinConvCoef(arkode_mem, 1.e-7)
Sundials.@checkflag Sundials.ARKodeSetMaxNumSteps(arkode_mem, 100000)
Sundials.@checkflag Sundials.ARKodeSetPredictorMethod(arkode_mem, 1)

Sundials.@checkflag Sundials.ARKodeSStolerances(arkode_mem, reltol, abstol)
A = Sundials.SUNDenseMatrix(neq,neq)
LS = Sundials.SUNDenseLinearSolver(y0,A)
Sundials.@checkflag Sundials.ARKDlsSetLinearSolver(arkode_mem, LS, A)

iout = 0
tout = t1
t = [t0]

while iout < nout
    y = similar(y0)
    flag = Sundials.ARKode(arkode_mem, tout, y, t, Sundials.ARK_NORMAL)
    println("T=", tout, ", Y=", y)
    iout += 1
    tout *= tmult
end

#Sundials.ARKodeFree(mem_ptr)
#Sundials.SUNLinSolFree_Dense(LS)
#Sundials.SUNMatDestroy_Dense(A)

ext: jl


text:
using LinearAlgebra
using SparseArrays
using Plots

"This routine generates a gif movie showing how the heat equation rapidly dissipates
high frequencies."

# m = number of grid points
function build_FD(m,f,α,β)
    x = LinRange(-1,1,m+2)
    xint = x[2:end-1]
    h = x[2]-x[1]
    A = (1/h^2) * spdiagm(0=>2*ones(m),-1=>-ones(m-1),1=>-ones(m-1))

    function F(t)
        b = f.(xint,t)
        b[1] += α(t)/h^2
        b[m] += β(t)/h^2
        return b
    end
    return A,x,xint,h,F
end

m = 200
T = .025

dt = .005/m

A,x,xint,h,F = build_FD(m,(x,t)->0,t->0,t->0)

u0(x) = sin(pi*x) + .1*sin(2*pi*x) - .33*sin(3*pi*x) + 10*sin(8*pi*x)
u = u0.(xint)
Nsteps = ceil(Int,T/dt)
dt = T/Nsteps
@gif for i = 1:Nsteps
    tprev = (i-1)*dt
    t = i*dt
    u .= (I + .5*dt*A)\((I - .5*dt*A)*u + .5*dt*(F(t)+F(tprev)))
    if i%10==0
        plot(xint,u,mark=:dot,leg=false)
        println("on timestep $i out of $Nsteps")
    end
end every 10


text:
"""
    initialise_updates!(store, ::updateInfo)
provides the time = 0 values of relevant quantities (e.g. loss, grad) in store. these are used for the time = 1 update.
time = 1 quantities are calculated after the time=1 update.

When writing a new initialise_updates! method for a different type of ::updateInfo, make sure to initialise (in store),
exactly those fields which are automatically updated by the do_updates! method
"""
function initialise_updates!(store::Dict{Symbol,Any}, upd::fixedNormFluctuatingNoisyGD)

    ## initialise storage holders
    data = get_training_batch(store[:trainingInfo])
    store = get_loss_and_grad!(store, data) #sets loss and grad
    store[:sys_noise], store[:intr_noise], store[:update] = make_some_model_templates(store[:lossInfo],3).value

    ## initialise time zero values for storage holders
    # set noise
    for s in (store[:sys_noise],store[:intr_noise])
        avoid_activations(s)
    end

    g_norm2, sys_noise_norm2, intr_noise_norm2 = Float32(0),Float32(0),Float32(0) #squared quantities
    for (k,p) in  enumerate(store[:params])
        if length(size(p)) > 1 #avoid activation layers
            g_norm2 += sum(store[:grad][p].^2)
            sys_noise_norm2 += sum(store[:sys_noise][k].^2)
            intr_noise_norm2 += sum(store[:intr_noise][k].^2)
        end
    end

    store[:grad_norm] = sqrt(Tracker.data(g_norm2))
    store[:sys_noise_norm] = sqrt(Tracker.data(sys_noise_norm2))
    store[:intr_noise_norm] = sqrt(Tracker.data(intr_noise_norm2))

"""
Want γ1, γ2:
    γ1^2 + γ2^2 = upd.norm^2 :=m^2
    γ1/γ2 = upd.sn_ratio := r^2
Soln:
γ1 = sqrt(m^2(1 - 1/(1+r^2))) = m*sqrt(1 - (1/(1+r^2)))
γ2 = sqrt(m^2/1+r^2)
"""
    store[:γ1] = upd.norm*sqrt(      1 - ( 1/(1+upd.sn_ratio^2) )      )
    store[:γ2] = upd.norm*sqrt(      1/(1 + upd.sn_ratio^2)    )

    return store
end

function initialise_updates!(store::Dict{Symbol,Any}, upd::vanillaGD)
    data = get_training_batch(store[:trainingInfo])
    store = get_loss_and_grad!(store, data)
    store[:update] = make_some_model_templates(store[:lossInfo],1).value[1]
    return store
end

function initialise_updates!(store::Dict{Symbol,Any}, upd::beelineLearning)
    data = get_training_batch(store[:trainingInfo])
    store[:loss] = get_loss(store[:lossInfo], data)
    store[:intr_noise], store[:update] = make_some_model_templates(store[:lossInfo],2).value
    avoid_activations(store[:intr_noise])
    store[:intr_noise_norm] = sqrt(sum(sum(el.^2) for el in store[:intr_noise]))
    store[:teacher_array] = [Tracker.data(el) for el in params(get_teacher(store[:lossInfo]))]
    store[:offset] = store[:teacher_array] - [Tracker.data(el) for el in store[:params]]
    store[:offset_norm] = sqrt(sum(sum(el.^2) for el in store[:offset]))

    return store
end

"""
    do_updates(store, ::updateInfo)

mutates values in the dictionary store, using loss function in ::lossInfo and update rule in ::updateInfo
"""
function do_updates!(store::Dict{Symbol,Any}, upd::fixedNormFluctuatingNoisyGD)
    # for fixed Norm noisy gradient descent with an intrinsic plasticity term
    ## seed noise

# multiply things by approprite norms: find multipliers
    sys_noise_multiplier = store[:γ2]/store[:sys_noise_norm]
    grad_multiplier = -store[:γ1]/store[:grad_norm]
    # signal_proportion = upd.sn_ratio/(1 + upd.sn_ratio) # in terms of signal:noise ratio
    # sys_noise_multiplier = sqrt((upd.norm^2)*signal_proportion)/(store[:sys_noise_norm])
    # grad_multiplier = -sqrt((upd.norm^2)*(1-signal_proportion))/(store[:grad_norm])
    intr_noise_multiplier = upd.intrinsic_norm/store[:intr_noise_norm]

# update parameters
    for (k,p) in enumerate(store[:params])
        if length(size(p)) > 1 #avoid activation layers
            store[:update][k] = Tracker.data(grad_multiplier*store[:grad][p]) .+
                                    sys_noise_multiplier*store[:sys_noise][k] .+
                                    intr_noise_multiplier*store[:intr_noise][k]
            Tracker.update!(p,store[:update][k])
        end
    end
    ## get loss and gradient
    data = get_training_batch(store[:trainingInfo])
    store = get_loss_and_grad!(store, data)
    store[:loss] = Tracker.data(store[:loss])

    ## set noise for next iteration
    for s in (store[:sys_noise],store[:intr_noise])
        avoid_activations(s)
    end

    ### calculate norms of things for next iteration
    g_norm2, sys_noise_norm2, intr_noise_norm2 = Float32(0),Float32(0),Float32(0) #squared quantities
    for (k,p) in  enumerate(store[:params])
        if length(size(p)) > 1 #avoid activation layers
            g_norm2 += sum(store[:grad][p].^2)
            sys_noise_norm2 += sum(store[:sys_noise][k].^2)
            intr_noise_norm2 += sum(store[:intr_noise][k].^2)
        end
    end
    store[:grad_norm] = sqrt(Tracker.data(g_norm2))
    store[:sys_noise_norm] = sqrt(Tracker.data(sys_noise_norm2))
    store[:intr_noise_norm] = sqrt(Tracker.data(intr_noise_norm2))
    return store
end

function do_updates!(store::Dict{Symbol,Any}, upd::vanillaGD)
    # update parameters
    for (k,p) in enumerate(store[:params])
        if length(size(p)) > 1 #avoid activation layer
            store[:update][k] = Tracker.data(-upd.learning_rate*store[:grad][p])
            Tracker.update!(p, store[:update][k])
            # Tracker.update!(p, -upd.learning_rate*store[:grad][p])
            # store[:update][k] = Tracker.data(-upd.learning_rate*store[:grad][p])
        end
    end

    # recalculate loss and grads at new parameter value
    data = get_training_batch(store[:trainingInfo])
    store = get_loss_and_grad!(store, data)
end

function do_updates!(store::Dict{Symbol,Any}, upd::beelineLearning)

    for (k,p) in enumerate(store[:params])
        if length(size(p)) > 1
            store[:update][k] = (upd.step_size/store[:offset_norm])*store[:offset][k] +
            (upd.intrinsic_norm/store[:intr_noise_norm])*store[:intr_noise][k]
            Tracker.update!(p, store[:update][k])

            ## update offset for next iteration
            store[:offset][k] = store[:teacher_array][k] - Tracker.data(p)
        end
    end

    avoid_activations(store[:intr_noise]) # update noise for next iter
    store[:intr_noise_norm] = sqrt(sum(sum(el.^2) for el in store[:intr_noise]))
    store[:offset_norm] = sqrt(sum(sum(el.^2) for el in store[:offset]))
    data = get_training_batch(store[:trainingInfo])
    store[:loss] = get_loss(store[:lossInfo], data)
    return store
end

function get_test_loss!(::updateInfo, store)
    if store[:am_i_updated][:test_loss] == true
        return store
    end
    store[:test_loss] = get_loss(store[:lossInfo], get_test_batch(store[:trainingInfo]))
    store[:am_i_updated][:test_loss] = true
    return store
end

function get_grad_norm!(::updateInfo, store)
    if store[:am_i_updated][:grad_norm] == true
        return store
    end
    store[:grad_norm] = Tracker.data(sqrt(sum(sum(store[:grad][p].^2) for p in store[:params])))
    store[:am_i_updated][:grad_norm] = true
    return store
end

function get_update_norm!(::updateInfo, store)
    if store[:am_i_updated][:update_norm] == true
        return store
    end
    store[:update_norm] = sqrt(sum(sum(el.^2) for el in store[:update]))
    store[:am_i_updated][:update_norm] = true
    return store
end

function get_test_grad!(::updateInfo, store)
    if store[:am_i_updated][:test_grad] == true
        return store
    end
    store = get_test_loss!(store[:updateInfo], store)
    store[:test_grad] = Tracker.gradient( () -> store[:test_loss], store[:params]; nest=true)
    store[:am_i_updated][:test_grad] = true
    return store
end

function get_full_training_loss!(::updateInfo, store)
    if store[:am_i_updated][:full_training_loss] == true
        return store
    end
    store[:full_training_loss] = get_loss(store[:lossInfo], store[:trainingInfo].training_data)
    store[:am_i_updated][:full_training_loss] = true
    return store
end

function get_full_training_grad!(::updateInfo, store)
    if store[:am_i_updated][:full_training_grad] == true
        return store
    end
    store = get_full_training_loss!(store[:updateInfo], store)
    store[:full_training_grad] = Tracker.gradient( () -> store[:full_training_loss], store[:params] )
    store[:am_i_updated][:full_training_grad] = true
    return store
end

function get_full_training_grad_norm!(::updateInfo, store)
    if store[:am_i_updated][:full_training_grad_norm] == true
        return store
    end
    store = get_full_training_grad!(store[:updateInfo], store)
    store[:full_training_grad_norm] = sqrt(sum(sum(store[:full_training_grad][p].^2) for p in store[:params]))
    store[:am_i_updated][:full_training_grad_norm] = true
    return store
end

function get_test_grad_norm!(::updateInfo, store)
    if store[:am_i_updated][:test_grad_norm] == true
        return store
    end
        store = get_test_grad!(store[:updateInfo], store)
        store[:test_grad_norm] = sqrt(sum(sum(store[:test_grad][p].^2) for p in store[:params]))
        store[:am_i_updated][:test_grad_norm] = true
    return store
end

"""
    get_local_task_difficulty(store::Dict{Symbol,Any}, tr::trainingInfo, li::lossInfo)
note that this function requires that store[:update_norm] is updated in the do_updates!() function, not post-hoc
"""
function get_local_task_difficulty!(::updateInfo, store)
    tr = store[:trainingInfo]
    li = store[:lossInfo]
    store = get_test_grad!(store[:updateInfo], store)
    store = get_test_grad_norm!(store[:updateInfo], store)
    store = get_update_norm!(store[:updateInfo], store)

    #  Hessian*update_direction
    store[:extra_tplt] = Tracker.gradient(() -> sum(sum(sum(store[:test_grad][p].*store[:update][k])) for (k,p) in enumerate(store[:params])), store[:params])

    store[:local_task_difficulty] = Tracker.data(
                        (1/(2*store[:test_grad_norm]*store[:update_norm]^2))*
                        (sum(sum(sum(store[:extra_tplt][p].*store[:update][k]))
                                        for (k,p) in enumerate(store[:params]))) )
    store[:am_i_updated][:local_task_difficulty] = true
    return store
end

function get_training_test_grad_correlation!(::updateInfo, store)
    # store = get_grad_norm!(store[:updateInfo], store)
    store = get_test_grad_norm!(store[:updateInfo], store) # also gets test grad
    store = get_full_training_grad_norm!(store[:updateInfo], store)
    temp = 0.
    for p in store[:params]
        temp += dot(store[:full_training_grad][p], store[:test_grad][p])
    end
    store[:training_test_grad_correlation] = Tracker.data(temp/(store[:full_training_grad_norm]*store[:test_grad_norm]) )
    store[:am_i_updated][:training_test_grad_correlation] = true
    return store
end

"""
solution to double updating:
1. make updated = make_update_checker part of store
2. if local_task_difficulty comes first: when it evaluates store[:test_loss], change updated(:test_loss) to true
3. if it comes after test_loss, don't need to calculate test loss

so an if statement in local_task_difficulty
"""


text:
module RuntimeGeneratedFunctions

using ExprTools, Serialization, SHA

const function_cache = Dict{Tuple,Expr}()
struct RuntimeGeneratedFunction{uuid,argnames}
    function RuntimeGeneratedFunction(ex)
        def = splitdef(ex)
        args, body = normalize_args(def[:args]), def[:body]
        uuid = expr2bytes(body)
        function_cache[uuid] = body
        new{uuid,Tuple(args)}()
    end
end
(f::RuntimeGeneratedFunction)(args::Vararg{Any,N}) where N = generated_callfunc(f, args...)

@inline @generated function generated_callfunc(f::RuntimeGeneratedFunction{uuid,argnames},__args...) where {uuid,argnames}
    setup = (:($(argnames[i]) = @inbounds __args[$i]) for i in 1:length(argnames))
    quote
        $(setup...)
        $(function_cache[uuid])
    end
end

export RuntimeGeneratedFunction

###
### Utilities
###
normalize_args(args::Vector) = map(normalize_args, args)
normalize_args(arg::Symbol) = arg
function normalize_args(arg::Expr)
    arg.head === :(::) || error("argument malformed. Got $arg")
    arg.args[1]
end

function expr2bytes(ex)
    io = IOBuffer()
    Serialization.serialize(io, ex)
    return Tuple(sha512(take!(io)))
end

end


text:
"
An iid mean zero normal dataset.

##Fields

- `yobs` The observed data.
- `nobs` Number of observations.
"
type IIDNormalDetails <: ApplicationDetails
    yobs::Array{Float64, 1}
    nobs::Int
end

"
A `IIDNormalParticles` variable stores the particles used in an SMC algorithm for the iid normal application.

##Fields

- `nparticles` How many particles are used.
- `scale` The scale parameter (must be >=0).
- `precision` The precision parameter i.e. `scale^-2.0`
- `u` Underlying uniform draws. `u[i]` is a vector of draws for the ith particle, controlling `y[i]`.
- `y` The simulated observation. `y[i]` is the simulated data for the ith particle.
- `accrate` Acceptance rate in the most recent iteration (fixed to 1.0 here)
- `iteration` How many SMC iterations have been performed.
- `reachedBound` Whether the lower bound on log likelihood has been reached.
- `schedule` Schedule of epsilon values used (or, for an adaptive algorithm, schedule so far)
- `maxz` Maximum slice sampling step size in previous iteration
"
type IIDNormalParticles <: Particles
    nparticles::Int
    scale::Float64
    precision::Float64
    u::Array{Array{Float64, 1}, 1}
    y::Array{Array{Float64, 1}, 1}
    accrate::Float64
    llest::Float64
    ε::Float64
    iteration::Int
    reachedBound::Bool
    schedule::Array{Float64, 1}
    maxz::Float64
end

"
Converts a vector of uniform draws to a vector of normal draws with the specified scale
"
function u2y(u::Array{Float64, 1}, scale::Float64)
    return quantile(Normal(), u) .* scale
end

function Particles(a::IIDNormalDetails, nparticles::Int, parameters::Array{Float64, 1})
    u = Array{Float64, 1}[rand(a.nobs) for i in 1:nparticles]
    y = Array{Float64, 1}[u2y(uvec, parameters[1]) for uvec in u]
    IIDNormalParticles(nparticles, parameters[1], parameters[1]^-2.0, u, y, 1.0, 0.0, Inf, 0, false, Float64[], 0.5)
end

"
Returns Euclidean distance between data vector `y` and observations.
"
function getDistance(a::IIDNormalDetails, y::Array{Float64, 1})
    sqrt(sum((y-a.yobs) .^ 2))
end

"
Returns vector of Euclidean distances between each particle and the observations.
"
function getDistances(a::IIDNormalDetails, p::IIDNormalParticles)
    Float64[getDistance(a, p.y[i]) for i in 1:p.nparticles]
end

"
Returns `NormalParticles` variable which is a copy of `p` except that its ith particle is a copy of the `indices[i]`th old particle.
"
function shuffle!(p::IIDNormalParticles, indices::Array{Int, 1})
    p.u = deepcopy(p.u[indices])
    p.y = deepcopy(p.y[indices])
    return nothing
end

"
Performs a slice sampling update on all particles.
"
function updateParticles!(a::IIDNormalDetails, p::IIDNormalParticles, θ::Array{Float64, 1})
    function dosim(u::Array{Float64, 1})
        y = u2y(u, θ[1])
        d = getDistance(a, y)
        return (y, d)
    end
    z = zeros(p.nparticles)
    for i in 1:p.nparticles
        p.u[i], p.y[i], z[i] = slicesample(p.u[i], dosim, p.ε, w=2.0*p.maxz)
    end
    p.maxz = maximum(abs(z))
    return nothing
end


text:
export zoo_sdd, zoo_sdd_file,
    SddFormat, SddVtreeFormat

struct SddFormat <: FileFormat end

const SddVtreeFormat = Tuple{SddFormat,VtreeFormat}
Tuple{SddFormat,VtreeFormat}() = (SddFormat(),VtreeFormat())

##############################################
# Read SDDs
##############################################

zoo_sdd_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/sdds/$name"

"""
    zoo_sdd(name)

Loads SDD file with given name from model zoo. See https://github.com/UCLA-StarAI/Circuit-Model-Zoo.
"""
zoo_sdd(name) =
    read(zoo_sdd_file(name), LogicCircuit, SddFormat())

const sdd_grammar = raw"""
    start: header (_NL node)+ _NL?

    header : "sdd" _WS INT

    node : "F" _WS INT -> false_node
         | "T" _WS INT -> true_node
         | "L" _WS INT _WS INT _WS SIGNED_INT -> literal_node
         | "D" _WS INT _WS INT _WS INT _WS elems  -> decision_node

    elems : elem (_WS elem)*
    elem : INT _WS INT

    %import common.INT
    %import common.SIGNED_INT
    %import common.WS_INLINE -> _WS
    %import common.NEWLINE -> _NL
    """ * dimacs_comments

const sdd_parser = Lark(sdd_grammar)

abstract type SddParse <: JuiceTransformer end

@inline_rule header(t::SddParse, x) =
    Base.parse(Int,x)

@rule start(t::SddParse, x) = begin
    @assert sdd_num_nodes_leafs(x[end]) == x[1]
    x[end]
end

@rule elem(t::SddParse, x) =
    [t.nodes[x[1]], t.nodes[x[2]]]

@rule elems(t::SddParse, x) =
    Array(x)

#  parse unstructured
struct PlainSddParse <: SddParse
    nodes::Dict{String,PlainLogicCircuit}
    PlainSddParse() = new(Dict{String,PlainLogicCircuit}())
end

@rule literal_node(t::PlainSddParse, x) =
    t.nodes[x[1]] = PlainLiteralNode(Base.parse(Lit,x[3]))

@rule false_node(t::PlainSddParse, x) =
    t.nodes[x[1]] = PlainConstantNode(false)

@rule true_node(t::PlainSddParse, x) =
    t.nodes[x[1]] = PlainConstantNode(true)

@rule decision_node(t::PlainSddParse,x) = begin
    @assert length(x[4]) == Base.parse(Int,x[3])
    elems = map(x[4]) do elem
        Plain⋀Node(elem)
    end
    t.nodes[x[1]] = Plain⋁Node(elems)
end

function Base.parse(::Type{PlainLogicCircuit}, str, ::SddFormat)
    ast = Lerche.parse(sdd_parser, str)
    Lerche.transform(PlainSddParse(), ast)
end

Base.read(io::IO, ::Type{PlainLogicCircuit}, ::SddFormat) =
    parse(PlainLogicCircuit, read(io, String), SddFormat())

#  parse structured
struct StructSddParse <: SddParse
    id2vtree::Dict{String,<:Vtree}
    nodes::Dict{String,StructLogicCircuit}
    StructSddParse(id2vtree) =
        new(id2vtree,Dict{String,StructLogicCircuit}())
end

@rule literal_node(t::StructSddParse, x) = begin
    lit = Base.parse(Lit,x[3])
    vtree = t.id2vtree[x[2]]
    t.nodes[x[1]] = PlainStructLiteralNode(lit, vtree)
end

@rule false_node(t::StructSddParse, x) =
    t.nodes[x[1]] = PlainStructConstantNode(false)

@rule true_node(t::StructSddParse, x) =
    t.nodes[x[1]] = PlainStructConstantNode(true)

@rule decision_node(t::StructSddParse,x) = begin
    @assert length(x[4]) == Base.parse(Int,x[3])
    vtree = t.id2vtree[x[2]]
    elems = map(x[4]) do elem
        PlainStruct⋀Node(elem[1], elem[2], vtree)
    end
    t.nodes[x[1]] = PlainStruct⋁Node(elems, vtree)
end

function Base.parse(::Type{PlainStructLogicCircuit}, str::AbstractString, ::SddFormat, id2vtree)
    ast = Lerche.parse(sdd_parser, str)
    Lerche.transform(StructSddParse(id2vtree), ast)
end

function Base.parse(::Type{PlainStructLogicCircuit}, strings, format::SddVtreeFormat)
    id2vtree = parse(Dict{String,Vtree}, strings[2], format[2])
    parse(PlainStructLogicCircuit, strings[1], format[1], id2vtree)
end

Base.read(io::IO, ::Type{PlainStructLogicCircuit}, ::SddFormat, id2vtree) =
    parse(PlainStructLogicCircuit, read(io, String), SddFormat(), id2vtree)

function Base.read(ios::Tuple{IO,IO}, ::Type{PlainStructLogicCircuit}, ::SddVtreeFormat)
    circuit_str = read(ios[1], String)
    vtree_str = read(ios[2], String)
    parse(PlainStructLogicCircuit, (circuit_str,vtree_str), SddVtreeFormat())
end

#  parse as SDD

function Base.parse(::Type{Sdd}, str::AbstractString, ::SddFormat, id2mgr::Dict{String,<:SddMgr})
    ast = Lerche.parse(sdd_parser, str)
    # create a structured logic circuit for Sdd Mgr
    struct_circuit = Lerche.transform(StructSddParse(id2mgr), ast)
    # turn logic circuit into Sdd
    compile(vtree(struct_circuit)::SddMgr, struct_circuit)::Sdd
end

function Base.parse(::Type{Sdd}, strings, format::SddVtreeFormat)
    id2mgr = parse(Dict{String,SddMgr}, strings[2], format[2])
    parse(Sdd, strings[1], format[1], id2mgr)
end

Base.read(io::IO, ::Type{Sdd}, ::SddFormat, id2mgr::Dict{String,<:SddMgr}) =
    parse(Sdd, read(io, String), SddFormat(), id2mgr)

function Base.read(ios::Tuple{IO,IO}, ::Type{Sdd}, ::SddVtreeFormat)
    circuit_str = read(ios[1], String)
    vtree_str = read(ios[2], String)
    parse(Sdd, (circuit_str,vtree_str), SddVtreeFormat())
end

##############################################
# Write SDDs
##############################################

const SDD_FORMAT = """c this file was saved by LogicCircuits.jl
c ids of sdd nodes start at 0
c sdd nodes appear bottom-up, children before parents
c
c file syntax:
c sdd count-of-sdd-nodes
c F id-of-false-sdd-node
c T id-of-true-sdd-node
c L id-of-literal-sdd-node id-of-vtree literal
c D id-of-decomposition-sdd-node id-of-vtree number-of-elements {id-of-prime id-of-sub}*
c"""

function Base.write(io::IO, sdd::Sdd, ::SddFormat, vtree2id::Function = (x -> 0))

    id = -1

    println(io, SDD_FORMAT)
    println(io, "sdd $(sdd_num_nodes_leafs(sdd))")

    f_con(n) = begin
        nid = id += 1
        sign = isfalse(n) ? "F" : "T"
        println(io, "$sign $nid")
        nid
    end

    f_lit(n) = begin
        nid = id += 1
        println(io, "L $nid $(vtree2id(mgr(n))) $(literal(n))")
        nid
    end

    f_a(n, ids) = tuple(ids...)

    f_o(n, ids) = begin
        nid = id += 1
        print(io, "D $nid $(vtree2id(mgr(n))) $(length(ids))")
        for el in ids
            print(io, " $(el[1]) $(el[2])")
        end
        println(io)
        nid
    end

    foldup_aggregate(sdd, f_con, f_lit, f_a, f_o, Union{Int, Tuple{Int,Int}})

    nothing
end

function Base.write(ios::Tuple{IO,IO}, sdd::Sdd, format::SddVtreeFormat)
    vtree2id = write(ios[2], mgr(sdd), format[2])
    write(ios[1], sdd, format[1], i -> vtree2id[i])
end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
"""
Convenience alias for a matrix of a sub-type of `Real`
"""
const MT = AbstractMatrix{<:Real}

"""
Convenience alias for a vector of a sub-type of `Real`
"""
const VT = AbstractVector{<:Real}

"""
A bilevel linear optimization problem of the form:
```
min cx^T * x + cy^T * y
s.t. G x + H y <= q
     x_j ∈ [xl_j,xu_j]
     x_j ∈ ℤ ∀ j ∈ Jx
     y ∈ arg min {
        d^T * y + x^T * F * y
        s.t. A x + B y <= b
             y_j ∈ [yl_j,yu_j]
        }
```
Note that integer variables are allowed at the upper level.
"""
mutable struct BilevelLP{V<:VT,M<:MT,MQ<:MT,VB<:AbstractVector{Bool}}
    cx::V
    cy::V
    G::M
    H::M
    q::V
    d::V
    A::M
    B::M
    b::V
    nu::Int
    nl::Int
    mu::Int
    ml::Int
    xl::V
    xu::V
    yl::VB # if lowerbound yi >= 0, otherwise yi >= -∞. For other bounds, add a constraint
    yu::V
    Jx::Vector{Int} # ∀ j ∈ Jx, x[j] is integer
    F::MQ

    function BilevelLP(cx::V,
                       cy::V,
                       G::M,
                       H::M,
                       q::V,
                       d::V,
                       A::M,
                       B::M,
                       b::V,
                       Jx::Vector{Int} = Int[],
                       F::MQ = zeros(length(cx),length(cy)); ylowerbound = true) where {V<:VT,M<:MT,MQ<:MT}
        nu = length(cx)
        nl = length(cy)
        nl == length(d) || DimensionMismatch("Objectives")
        mu = length(q)
        ml = length(b)
        size(A) == (ml,nu) && size(B) == (ml,nl) || DimensionMismatch("Lower constraints")
        size(G) == (mu,nu) && size(H) == (mu,nl) || DimensionMismatch("Higher constraints")
        xl = zeros(nu)
        xu = Inf64 .* ones(nu)
        yl = zeros(Bool, nl) .+ ylowerbound # by default, y >= 0
        yu = Inf64 .* ones(nl)
        size(F) == (nu,nl) || DimensionMismatch("Quadratic constraint")
        new{V,M,MQ,Vector{Bool}}(cx,cy,G,H,q,d,A,B,b,nu,nl,mu,ml,xl,xu,yl,yu,Jx,F)
    end
end

"""
VariableType enum distinguishing upper- and lower-level
variables for setting upper and lower bounds
"""
@enum VariableType begin
    lower
    upper
end

"""
Set a lower bound on a lower or higher variable of `bp` depending on `vartype`
"""
function JuMP.set_lower_bound(bp::BilevelLP, vartype::VariableType, j::Integer, v::T) where {T<:Real}
    if vartype == lower::VariableType
        if v ≈ 0
            bp.yl[j] = true
        elseif isinf(v) && v < 0.
            bp.yl[j] = false
        else
            throw(DomainError(v,"The lower bound on a lower variable can only be 0 or -∞"))
        end
    else
        bp.xl[j] = v
    end
    return nothing
end

"""
Set an upper bound on a lower or higher variable of `bp` depending on `vartype`
Since we need the dual for lower-level upper bounds, this case is registered as
a constraint.
"""
function JuMP.setupperbound(bp::BilevelLP, vartype::VariableType, j::Integer, v::T) where {T<:Real}
    if vartype == upper::VariableType
        bp.xu[j] = v
    else
        row = zeros(1,bp.nl); row[1,j] = 1.
        bp.B = vcat(bp.B, row)
        bp.A = vcat(bp.A, zeros(1,bp.nu))
        push!(bp.b, v)
        bp.ml = bp.ml + 1
    end
end


text:
# This file is a part of UpROOT.jl, licensed under the MIT License (MIT).

using UpROOT
using Test

using ArraysOfArrays
using TypedTables

@testset "testfiles" begin
    @testset "leaflist" begin
        file = TFile(UpROOT.testfiles["leaflist"])
        @test "tree;1" in keys(file)
        tree = file["tree"]

        @test tree[1] isa NamedTuple
        @test tree[1].leaflist isa NamedTuple
        @test tree[1:5] isa Table
        @test tree[1:5].leaflist isa Table
        @test tree[:] isa Table
        @test tree[:].leaflist isa Table

        @test tree.leaflist[1] isa NamedTuple
        @test tree.leaflist[1:5] isa Table
        @test tree.leaflist[:] isa Table

        @test tree[1].leaflist == tree.leaflist[1]
    end

    @testset "HZZ" begin
        file = TFile(UpROOT.testfiles["HZZ"])
        @test "events;1" in keys(file)
        tree = file["events"]

        @test tree[1] isa NamedTuple
        @test tree[1:5] isa Table
        @test tree[:] isa Table

        @test tree.Jet_E[1] isa AbstractVector
        @test tree.Jet_E[1:5] isa VectorOfVectors
        @test tree.Jet_E[:] isa VectorOfVectors

        @test copy(tree.Jet_E) == tree.Jet_E[:]
        @test Array(tree.Jet_E) == tree.Jet_E[:]
        @test convert(Array, tree.Jet_E) == tree.Jet_E[:]

        @test Table(tree) == tree[:]

        @test tree[1:5, (:Jet_Px, :Jet_Py, :Jet_Pz)] isa Table{<:NamedTuple{(:Jet_Px, :Jet_Py, :Jet_Pz)}}
        @test tree[:, [:Jet_Px, :Jet_Py, :Jet_Pz]] isa Table{<:NamedTuple{(:Jet_Px, :Jet_Py, :Jet_Pz)}}
        @test tree[[2, 4, 7, 11], ["Jet_Px", "Jet_Py", "Jet_Pz"]] isa Table{<:NamedTuple{(:Jet_Px, :Jet_Py, :Jet_Pz)}}
    end
end


text:
"""
    intervals(a, [b]; endincluded=false, step=one(eltype(a)))

Return a vector of ranges from `a[i]:step:a[i+1]` or `a[i]:step:b[j]`, if `b` is given.

All ranges returned will have a non-zero length and will not overlap (except by one `step` if endincluded is `true`). Values in `a` and/or `b` will be skipped when required to enforce this behavior.

If `endincluded == true`, ranges will include the `a[i+1]` or `b[j]` value, otherwise, ranges will end one `step` before that value.

"""
function intervals(a::AbstractVector{T}; endincluded=false, step=one(T)) where T
    length(a) > 1 || throw(ArgumentError("a is not long enough to create any intervals"))
    ranges = Vector{typeof(first(a):step:a[2])}()
    for i in 1:(length(a)-1)
        if endincluded
            rg = a[i]:step:a[i+1]
        else
            rg = (a[i]:step:a[i+1])[1:end-1]
        end
        push!(ranges, rg)
    end

    return ranges
end

function intervals(
    a::AbstractVector{T}, b::AbstractVector{T}; endincluded=false, step=one(T)
) where T
    ranges = Vector{typeof(first(a):step:first(b))}()

    firstbi = searchsortedfirst(b, first(a))

    bi = firstbi
    ai = firstindex(a)
    while bi < lastindex(b) && ai < lastindex(a)
        bi = searchsortedfirst(b, a[ai], bi, lastindex(b), Base.Order.Forward)
        ai = searchsortedlast(a, b[bi], ai, lastindex(a), Base.Order.Forward)
        if endincluded
            rg = a[ai]:step:b[bi]
        else
            rg = (a[ai]:step:b[bi])[1:end-1]
        end
        push!(ranges, rg)

        bi += 1
        ai += 1
    end

    return ranges
end

function _rotating_diff(vectors::Vararg{<:AbstractVector{T},N};
    i=argmin(first.(vectors))
) where {T,N}
    diffed = ntuple(x -> similar(Vector{T}, (0,)), N)
    iterators::Vector{Union{Nothing,Tuple{T,Any}}} = collect(iterate.(vectors))

    nextvector = i
    currvector = nextvector
    nextvector = mod1(nextvector+1, N)

    while !isnothing(iterators[nextvector])
        cval, cstate = iterators[currvector]
        nval, nstate = iterators[nextvector]
        iterators[currvector] = iterate(vectors[currvector], cstate)

        push!(diffed[currvector], nval - cval)
        currvector = nextvector
        nextvector = mod1(nextvector+1, N)
    end

    return diffed
end

function interleave(vectors::Vararg{<:AbstractVector{T},N}) where {T,N}
    minlen = minimum(length.(vectors))
    arr = Vector{T}(undef, N*minlen)
    iterators::Vector{Union{Nothing,Tuple{T,Any}}} = collect(iterate.(vectors))

    currvector = 1
    i = 0

    while !isnothing(iterators[currvector])
        cval, cstate = iterators[currvector]
        iterators[currvector] = iterate(vectors[currvector], cstate)

        arr[i += 1] = cval
        currvector = mod1(currvector+1, N)
    end

    return arr
end


text:
function new_f(x::Int, y::Int)
    if y == 0
        error("division by zero is not defined")
    else
        a = x / y
    end
    return a
end

max_stars_repo_name: Mikejmnez/MyExample

max_issues_repo_name: Mikejmnez/MyExample

max_forks_repo_name: Mikejmnez/MyExample


text:
"""
# module Namelists

# Examples

```jldoctest
julia>
```
"""
module Namelists

export NamelistVariable,
    Namelist,
    CONTROL_NAMELIST

struct NamelistVariable{T}
    name::String
    in_namelist::String
    value::T
end

NamelistVariable(name::String, in_namelist::String, value::T) where T = NamelistVariable{T}(name, in_namelist, value)
NamelistVariable{T}(name::String, in_namelist::String) where T = value::T -> NamelistVariable(name, in_namelist, value)

struct Namelist
    name::String
    values::Dict
end

const CONTROL_NAMELIST_VARIABLES = [
    ("calculation", String),
    ("title", String),
    ("verbosity", String),
    ("restart_mode", String),
    ("wf_collect", Bool),
    ("nstep", Int),
    ("iprint", Int),
    ("tstress", Bool),
    ("tprnfor", Bool),
    ("dt", Float64),
    ("outdir", String),
    ("wfcdir", String),
    ("prefix", String),
    ("lkpoint_dir", Bool),
    ("max_seconds", Float64),
    ("etot_conv_thr", Float64),
    ("forc_conv_thr", Float64),
    ("disk_io", String),
    ("pseudo_dir", String),
    ("tefield", Bool),
    ("dipfield", Bool),
    ("lelfield", Bool),
    ("nberrycyc", Int),
    ("lorbm", Bool),
    ("lberry", Bool),
    ("gdir", Int),
    ("nppstr", Int),
    ("lfcpopt", Bool),
    ("gate", Bool)
]

const CONTROL_NAMELIST = Namelist("CONTROL", Dict(name => NamelistVariable{T}(name, "CONTROL") for (name, T) in CONTROL_NAMELIST_VARIABLES))

end


text:
using Base.Threads: @spawn, @threads
using LinearAlgebra: norm
using StaticArrays
using Test
using FLoops

# Parallel code for sum(norm.( x - y for x in x, y in y ))
function sumd_threads(x,y,nbatches,aux)
    @assert length(x)%nbatches == 0
    batchsize = length(x) ÷ nbatches
    aux .= 0.
    @threads for ibatch in 1:nbatches
        ifirst = (ibatch-1)*batchsize + 1
        ilast = ibatch*batchsize
        begin
            acc = 0.
            for i in ifirst:ilast
                for j in 1:length(y)
                    @inbounds acc += norm(y[j]-x[i])
                end
            end
            aux[ibatch] = acc
        end
    end
    return sum(aux)
end

function sumd_spawn(x,y,nbatches,aux)
    @assert length(x)%nbatches == 0
    batchsize = length(x) ÷ nbatches
    aux .= 0
    @sync for ibatch in 1:nbatches
        ifirst = (ibatch-1)*batchsize + 1
        ilast = ibatch*batchsize
        @spawn begin
            acc = 0.
            for i in ifirst:ilast
                for j in 1:length(y)
                    @inbounds acc += norm(y[j]-x[i])
                end
            end
            aux[ibatch] = acc
        end
    end
    return sum(aux)
end

function sumd_floop(x,y,nbatches)
    @floop ThreadedEx(basesize = div(length(x),nbatches)) for i in 1:length(x)
        for j in 1:length(y)
            d = norm(y[j]-x[i])
            @reduce( dsum = 0. + d )
        end
    end
    return dsum
end

function test(nbatches)
    x = [ rand(SVector{3,Float64}) for _ in 1:6400 ]
    y = [ rand(SVector{3,Float64}) for _ in 1:100 ]
    aux = zeros(nbatches)
    simple = sum(norm.( x - y for x in x, y in y ))
    @test sumd_threads(x,y,nbatches,aux) ≈ simple
    @test sumd_spawn(x,y,nbatches,aux) ≈ simple
    @test sumd_floop(x,y,nbatches) ≈ simple
    print("sumd threads: "), @btime sumd_threads($x,$y,$nbatches,$aux)
    print("sumd spawn:   "), @btime sumd_spawn($x,$y,$nbatches,$aux)
    print("floops:       "), @btime sumd_floop($x,$y,$nbatches)
    nothing
end


text:
import LinearAlgebra
import MLBase

"""
"""
is_square(m::AbstractMatrix) = size(m, 1) == size(m, 2)

"""
"""
function cohen_kappa(contingency_table::AbstractMatrix)
    if !is_square(contingency_table)
        error("contingency_table must be a square matrix")
    end
    numclasses = size(contingency_table, 1)
    if numclasses < 2
        error("number of classes must be >= 2")
    end
    totalnumcases = sum(contingency_table)
    # p_o = "relative observed agreement among raters (identical to accuracy)"
    p_o = sum(LinearAlgebra.diag(contingency_table))/totalnumcases
    # n_k1 = number of times rater 1 predicted class k
    # n_k2 = number of times rater 2 predicted class k
    # rater 1 = rows
    # rater 2 = columns
    n_k1 = [sum( contingency_table[k,:] ) for k = 1:numclasses]
    n_k2 = [sum( contingency_table[:,k] ) for k = 1:numclasses]
    # p_e = "hypothetical probability of chance agreement"
    p_e = sum(n_k1 .* n_k2)/(totalnumcases^2)
    kappa = (p_o - p_e)/(1 - p_e)
    return kappa
end

"""
"""
function compute_contingency_table(y1::AbstractVector, y2::AbstractVector)
    classes = sort(unique(vcat(y1, y2)); rev = false,)
    numclasses = length(classes)
    if numclasses < 2
        error("number of classes must be >= 2")
    end
    return contingency_table
end

"""
"""
function cohen_kappa(y1::AbstractVector, y2::AbstractVector)
    contingency_table = compute_contingency_table(y1, y2)
    result = cohen_kappa(contingency_table)
    return result
end

"""
"""
function compute_contingency_table(rocnums::MLBase.ROCNums)
    # we will arbitrarily set rows = predicted, columns = true/gold
    contingency_table = [rocnums.tp rocnums.fp; rocnums.fn rocnums.tp]
    return contingency_table
end

"""
"""
function cohen_kappa(rocnums::MLBase.ROCNums)
    contingency_table = compute_contingency_table(rocnums)
    result = cohen_kappa(contingency_table)
    return result
end


text:
@testset "Greedy" begin
    seed!(42)
    N = 20
    D = 3
    noise = 0.01
    nInd = 10
    kernel = SqExponentialKernel()
    X = ColVecs(rand(D, N))
    y = rand(N)
    alg = Greedy(nInd, N)
    @test repr(alg) == "Greedy Selection of Inducing Points"
    Z = inducingpoints(alg, X; y=y, kernel=kernel, noise=noise)
    @test length(Z) == nInd
    @test_throws ArgumentError Greedy(-1, 2)
    @test_throws ArgumentError Greedy(10, -3)
    @test_throws ArgumentError inducingpoints(
        alg, X; y=rand(N + 1), kernel=kernel, noise=noise
    )
    @test_throws ArgumentError inducingpoints(alg, X; y=y, kernel=kernel, noise=-1)
    @test_throws ErrorException inducingpoints(
        Greedy(N + 1, N), X; y=y, kernel=kernel, noise=noise
    )
    test_Zalg(Greedy(nInd, N), N; kernel=kernel, noise=noise, y=y)
end


text:
module VisualRegressionTests

# Gtk is optional since it's only used for the GUI popup
try
    @eval using Gtk
catch err
    warn("Gtk not loaded. err: $err")
end

import Images

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

type VisualTest
    testFilename::AbstractString

    # function taking (output_filename, args...; kw...) and writing to output_filename
    testFunction::Function
    args
    kw

    referenceFilename::AbstractString
end

function VisualTest(testFunction::Function, referenceFilename::AbstractString, args...; kw...)
    VisualTest(tempname() * ".png", testFunction, args, kw, referenceFilename)
end

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

@enum VisualTestStatus EXACT_MATCH CLOSE_MATCH DOES_NOT_MATCH PROCESSING_ERROR

type VisualTestResult
    testFilename::AbstractString
    testImage
    referenceFilename::AbstractString
    referenceImage
    status::VisualTestStatus
    diff::Float64
    err
end

Base.success(result::VisualTestResult) = result.status in (EXACT_MATCH, CLOSE_MATCH)

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

export
    compare_images,
    test_images,
    VisualTest,
    VisualTestStatus,
    VisualTestResult,
    EXACT_MATCH,
    CLOSE_MATCH,
    DOES_NOT_MATCH,
    PROCESSING_ERROR

include("gui.jl")
include("imgcomp.jl")

end # module


text:
using Plots

f(x) = 1/sqrt(1-x^2)

x = -0.99:0.01:0.99

savefig("plot_function.png")


text:
"""
    constraining(f; onget=true, onset=true)

Create a lens to impose constraints by a callable `f`.

* The callable `f` must be idempotent; i.e., `f ∘ f = identity`.

* If the original object already satisfies the constraint
  (i.e. `f(obj) == obj`), `onget=false` can be passed safely to skip
  calling `f` during `get`.

# Examples
```jldoctest
julia> using Kaleido, Setfield

julia> obj = (a = 1, b = 1);

julia> constraint = constraining() do obj
           @set obj.b = obj.a
       end
constraining(#1)

julia> lens = constraint ∘ @lens _.a
constraining(#1) ∘ (@lens _.a)

julia> get(obj, lens)
1

julia> set(obj, lens, 2)
(a = 2, b = 2)
```

`constraining` is useful when combined with [`@batchlens`](@ref) or
[`MultiLens`](@ref):

```jldoctest
julia> using Kaleido, Setfield

julia> obj = (a = 1, b = 2, c = 3);

julia> constraint = constraining() do obj
           @set obj.c = obj.a + obj.b
       end;

julia> lens = constraint ∘ MultiLens((
           (@lens _.a),
           (@lens _.b),
       ));

julia> get(obj, lens)
(1, 2)

julia> set(obj, lens, (100, 20))
(a = 100, b = 20, c = 120)
```
"""
constraining(f; onget=true, onset=true) =
    ConstrainingLens(prefer_singleton_callable(f), asval(onget), asval(onset))

asval(v::Val) = v
asval(v) = Val(v)

const ValBool = Union{Val{true}, Val{false}}

struct ConstrainingLens{TF, TG <: ValBool, TS <: ValBool} <: Lens
    f::TF
    onget::TG
    onset::TS
end

Setfield.get(obj, lens::ConstrainingLens) =
    lens.onget === Val(true) ? lens.f(obj) : obj

Setfield.set(::Any, lens::ConstrainingLens, obj) =
    lens.onset === Val(true) ? lens.f(obj) : obj

function Base.show(io::IO, lens::ConstrainingLens)
    # TODO: handle module prefix
    print(io, "constraining(")
    print(io, lens.f)
    default = constraining(lens.f)
    if lens !== default
        # TODO: handle compat
        prefix = "; "
        if lens.onget !== default.onget
            print(io, prefix, "onget=", lens.onget === Val(true))
            prefix = ", "
        end
        if lens.onset !== default.onset
            print(io, prefix, "onset=", lens.onset === Val(true))
        end
    end
    print(io, ")")
end
Base.show(io::IO, ::MIME"text/plain", lens::ConstrainingLens) = show(io, lens)


text:
# ---
# title: 258. Add Digits
# id: problem258
# author: zhwang
# date: 2022-03-03
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/add-digits/description/>
# hidden: true
# ---
#
# Given a non-negative integer `num`, repeatedly add all its digits until the
# result has only one digit.
#
# **Example:**
#
#
#
#     Input: 38
#     Output: 2
#     Explanation: The process is like: 3 + 8 = 11, 1 + 1 = 2.
#                  Since 2 has only one digit, return it.
#
#
# **Follow up:**
# Could you do it without any loop/recursion in O(1) runtime?
#
#
## @lc code=start
using LeetCode

function add_digits(num::Int)::Int
    num <= 9 && return num
    return iszero(num % 9) ? 9 : num % 9
end
## @lc code=end

## @lc test=start
@testset "258.add-digits.jl" begin
    @test all(add_digits(i) == i for i in 0:9)
    @test add_digits(18) == 9
    @test add_digits(10) == 1
    @test add_digits(38) == 2
end
## @lc test=end

## @lc info=start
# link: [solution to 258]()
## @lc info=end


text:
const libwci = find_library(["libwrapclang",],[Pkg.dir("Clang", "deps", "usr", "lib"),])
@assert(libwci != "", "Failed to find required library libwrapclang. Try re-running the package script using Pkg.build(\"Clang\")")


text:
module Day05

import Base.show, Base.string
import ..data_dir # from parent module

export part1, part2

struct Point
    x::Int64
    y::Int64
end

function string(point::Point)
    "$(point.x),$(point.y)"
end

function show(io::IO, point::Point)
    println(io, string(point))
end

function Point(point_defn::AbstractString)
    x_string, y_string = split(point_defn, ",")
    x, y = parse(Int64, x_string), parse(Int64, y_string)
    Point(x, y)
end

struct Line
    a::Point
    b::Point
end

function string(line::Line)
    "$(string(line.a)) -> $(string(line.b))"
end

function show(io::IO, line::Line)
    println(io, string(line))
end

function Line(line_defn::AbstractString)
    point_defns = split(line_defn, " -> ")
    a_string, b_string = point_defns[1], point_defns[2]
    a, b = Point(a_string), Point(b_string)
    Line(a, b)
end

is_horizontal(line::Line) = (line.a.y == line.b.y)
is_vertical(line::Line) = (line.a.x == line.b.x)

"""
Julia sees no points for the range `5:1` since `5 > 1`. This function makes sure
the step is correct (either 1 or -1) so that `5:1 == [5, 4, 3, 2, 1]`.
"""
function sequence(i, j)
    if i > j
        range(i, j, step = -1)
    else
        range(i, j, step = 1)
    end
end

"""
Also known as _gradient_
"""
slope(line) = (line.b.y - line.a.y) / (line.b.x - line.a.x)

"""
Returns a dictionary indexed by points and with values the number of occurrences of
each points. Points that aren't on any of the given lines will not appear in the
dictionary at all.

Lines are either horizontal, vertical, or diagonal at a 45° angle. In the diagonal
case, the slope is ±1, so we don't need to worry about non-integral coordinates.

The vertical line case has to be handled separately because the slope is NaN due to
a division by zero
"""
function points_on_line(line::Line)
    if is_vertical(line)
        return ([Point(line.a.x, j) for j in sequence(line.a.y, line.b.y)])
    end

    slope_of_line = slope(line)
    points = Vector{Point}()
    for i in sequence(0, line.b.x - line.a.x)
        # Move 1 along and 1 up/down each time
        point = Point(
            Int(line.a.x + i),
            Int(line.a.y + (i * slope_of_line))
        )
        push!(points, point)
    end
    points
end

"""
Returns a dictionary indexed by points and with values the number of occurrences of
each points. Points that aren't on any of the given lines will not appear in the
dictionary at all.
"""
function points_on_lines(lines::Vector{Line})
    occurrences = Dict{Point,Int64}()
    for line in lines
        for point in points_on_line(line)
            occurrences[point] = get(occurrences, point, 0) + 1
        end
    end
    occurrences
end

function count_multiple_occurrences(lines::Vector{Line})
    occurrences = points_on_lines(lines)
    length([p for p in keys(occurrences) if occurrences[p] > 1])
end

# Parse the input
const data_file = joinpath(data_dir, "day05")
const input = readlines(data_file) # read into an array, each line is an element
const lines = [Line(i) for i in input]

function part1(input = input)
    lines = [Line(i) for i in input]
    no_diagonals = [line for line in lines if is_horizontal(line) || is_vertical(line)]
    count_multiple_occurrences(no_diagonals)
end
function part2(input = input)
    lines = [Line(i) for i in input]
    count_multiple_occurrences(lines)
end

end # module


text:
export alphaFilter, alphaSimplex, delaunayTriangulation
#===============================================================================
#
#	src/alpha_complex.jl
#
#	This File Contains:
#
#	 - alphaFilter(
#			V::Matrix,
#			DT = Array{Int64,1}[];
#			digits=64
#		)::DataStructures.SortedMultiDict{}
#
#	 - alphaSimplex(
#			V::Matrix,
#			filtration::DataStructures.SortedMultiDict{},
#			α_threshold::Float64
#		)::Array{Array{Array{Int64,1},1},1}
#
#	 - delaunayTriangulation(
#			V::Matrix
#		)::Array{Array{Int64,1},1}
#
===============================================================================#

"""
	alphaFilter(
		V::Matrix, DT = Array{Int64,1}[];
		digits=64
	)::DataStructures.SortedMultiDict{}

Return ordered collection of pairs `(alpha charatteristic, complex)`.

This method evaluates the `α`-filter over the sites `S`.
If a Delaunay Triangulation `DT` is not specified than it is evaluated
via `AlphaStructures.delaunayTriangulation()`.

# Examples
```jldoctest
julia> V = [1. 2. 1. 2.; 0. 0. 1. 2. ];

julia> AlphaStructures.alphaFilter(V)
SortedMultiDict(Base.Order.ForwardOrdering(),
	0.0 => [1],
	0.0 => [2],
	0.0 => [3],
	0.0 => [4],
	0.5 => [1, 2],
	0.5 => [1, 3],
	0.70710678118655 => [2, 3],
	0.70710678118655 => [3, 4],
	0.70710678118655 => [1, 2, 3],
	1.0 => [2, 4],
	1.0 => [2, 3, 4]
)

```
"""
function alphaFilter(
		V::Matrix,
		DT = Array{Int64,1}[];
		digits=64
	)::DataStructures.SortedDict{}

	dim = size(V, 1)
	filtration = DataStructures.SortedDict{Array{Int64,1},Float64}()

	# 1 - Each point => alpha_char = 0.
	for i = 1 : size(V, 2)
		insert!(filtration, [i], 0.)
	end

	# 2 - Delaunay triangulation of ``V``
	if isempty(DT)
		DT = AlphaStructures.delaunayTriangulation(V)
	end

	n_upsimplex = length(DT)

	# 3 - process all upper simplex
	ind = 1
	for upper_simplex in DT
		if ind % 500000 == 0
			println(ind," simplices processed of ", n_upsimplex)
		end
		AlphaStructures.processuppersimplex(V,upper_simplex,filtration; digits = digits)
		ind = ind + 1
	end

	return filtration
end

"""
	processuppersimplex(
		V::Matrix,
		up_simplex::Array{Int64,1},
		filtration::DataStructures.SortedDict{};
		digits=64
		)

Process the upper simplex.
"""
function processuppersimplex(
		V::Matrix,
		up_simplex::Array{Int64,1},
		filtration::DataStructures.SortedDict{};
		digits=64
		)

	α_char = AlphaStructures.findRadius(V[:, up_simplex], digits=digits);
	insert!(filtration, up_simplex, α_char)

	d = length(up_simplex)-1
	if d > 1
		# It gives back combinations in natural order
		newsimplex = collect(Combinatorics.combinations(up_simplex,d))
		for lowsimplex in newsimplex
			AlphaStructures.processlowsimplex(V, up_simplex, lowsimplex, filtration; digits=digits)
		end
	end
end

"""
	processlowsimplex(
		V::Matrix,
		up_simplex::Array{Int64,1},
		lowsimplex::Array{Int64,1},
		filtration::DataStructures.SortedDict{};
		digits=64
		)

Process the lower simplex knowing the upper.
"""
function processlowsimplex(
	V::Matrix,
	up_simplex::Array{Int64,1},
	lowsimplex::Array{Int64,1},
	filtration::DataStructures.SortedDict{};
	digits=64)

	α_char = AlphaStructures.findRadius(V[:, lowsimplex], digits=digits)
	point = V[:, setdiff(up_simplex, lowsimplex)]

	if AlphaStructures.vertexInCircumball(V[:, lowsimplex], α_char, point)
		filtration[lowsimplex] = filtration[up_simplex]

	elseif !haskey(filtration, lowsimplex)
		filtration[lowsimplex] = α_char

	end

	d = length(lowsimplex)-1
	if d > 1
		# It gives back combinations in natural order
		newsimplex = collect(Combinatorics.combinations(lowsimplex,d))
		for simplex in newsimplex
			 AlphaStructures.processlowsimplex(V, lowsimplex, simplex, filtration, digits=digits)
		end
	end
end

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

"""
	alphaSimplex(
		V::Matrix,
		filtration::DataStructures.SortedMultiDict{},
		α_threshold::Float64
	)::Array{Array{Array{Int64,1},1},1}

Return collection of all `d`-simplex, for `d ∈ [0,dimension]`,
	with characteristic α less than a given value `α_threshold`.
"""
function alphaSimplex(
		V::Matrix,
		filtration::DataStructures.SortedDict{},
		α_threshold::Float64
	)::Array{Array{Array{Int64,1},1},1}

	dim = size(V, 1)
	# [VV, EV, FV, ...]
	simplexCollection = [ Array{Array{Int64,1},1}() for i = 1 : dim+1 ]

	for (k, v) in filtration
        if v <= α_threshold
        	push!(simplexCollection[length(k)], k)
    	end
    end

	sort!.(simplexCollection)

	return simplexCollection
end


text:
_ft(id, day) = isodd(id + day)
function run!(::Experiment{:Amodio20_exp1}, models)
    results = DataFrame(id = Int[], action = String[], foodtype = String[], A = Int[], B = Int[], C = Int[])
    for id in eachindex(models)
        m = models[id]
        trays = [Tray(i) for i in 1:3] # A, B, C
        for day in 1:9
            remove!(m, Any)
            wait!(m, 12u"hr")
            add!(m, trays[_abc(id, day)])
            if _ft(id, day)
                add!(m, PowderedKibble, 50)
            else
                add!(m, PowderedPeanut, 50)
            end
            wait!(m, 15u"minute")
            remove!(m, Any)
            add!(m, MaintenanceDiet)
            if day == 9
                wait!(m, 3u"hr")
                remove!(m, Any)
                add!(m, trays)
                add!(m, Kibble, 25)
                add!(m, Peanut, 25)
                wait!(m, 15u"minute")
                remove!(m, Any)
                push!(results, [id; "cache"; "peanut"; countcache.(trays, Peanut)])
                push!(results, [id; "cache"; "kibble"; countcache.(trays, Kibble)])
            else
                wait!(m, 3u"hr")
                remove!(m, Any)
                add!(m, trays)
                wait!(m, 15u"minute")
                add!(m, MaintenanceDiet)
                wait!(m, 8.5u"hr")
            end
        end
    end
    results
end

# 1. Rome 2. Quito 3. Lisbon 4. Caracas 5. Wellington 6. Washington
_nf(id, day) = isodd(id + day)
_abc(id, day) = (day + div(id - 1, 2) - 1) % 3 + 1
function run!(::Experiment{:Amodio20_exp2}, models)
    results = DataFrame(id = Int[], action = String[], A = Int[], B = Int[], C = Int[])
    for id in eachindex(models)
        m = models[id]
        trays = [Tray(i) for i in 1:3] # A, B, C
        for day in 1:9
            remove!(m, Any)
            wait!(m, 12u"hr")
            add!(m, trays[_abc(id, day)])
            if _nf(id, day)
                add!(m, PowderedKibble, 50)
            end
            wait!(m, 15u"minute")
            remove!(m, Any)
            add!(m, MaintenanceDiet)
            if day == 9
                wait!(m, 3u"hr")
                remove!(m, Any)
                add!(m, trays)
                add!(m, Kibble, 50)
                wait!(m, 15u"minute")
                remove!(m, Any)
                push!(results, [id; "cache"; countcache.(trays)])
            else
                wait!(m, 3u"hr")
                remove!(m, Any)
                add!(m, trays)
                wait!(m, 15u"minute")
                add!(m, MaintenanceDiet)
                wait!(m, 8.5u"hr")
            end
        end
    end
    results
end


text:
#!/usr/bin/env julia

include("HHChannels.jl")

using JSON
using Unitful
using Unitful.DefaultSymbols
using Main.HHChannels

scale(quantity, unit) = uconvert(NoUnits, quantity/unit)

radius = 18.8µm/2
area = 4*pi*radius^2

stim = Stim(10ms, 100ms, 0.1nA/area)
ts, vs = run_hh(100ms, stim=stim, sample_dt=0.025ms)

trace = Dict(
    :name => "membrane voltage",
    :sim => "numeric",
    :model => "soma",
    :units => "mV",
    :data => Dict(
        :time => scale.(ts, 1ms),
        Symbol("soma.mid") => scale.(vs, 1mV)
    )
)

println(JSON.json([trace]))

max_stars_repo_name: kanzl/arbor

max_issues_repo_name: kanzl/arbor

max_forks_repo_name: kanzl/arbor


text:
# ---
# title: 1605. Find Valid Matrix Given Row and Column Sums
# id: problem1605
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/find-valid-matrix-given-row-and-column-sums/description/>
# hidden: true
# ---
#
# You are given two arrays `rowSum` and `colSum` of non-negative integers where
# `rowSum[i]` is the sum of the elements in the `ith` row and `colSum[j]` is the
# sum of the elements of the `jth` column of a 2D matrix. In other words, you do
# not know the elements of the matrix, but you do know the sums of each row and
# column.
#
# Find any matrix of **non-negative** integers of size `rowSum.length x
# colSum.length` that satisfies the `rowSum` and `colSum` requirements.
#
# Return _a 2D array representing **any** matrix that fulfills the
# requirements_. It's guaranteed that **at least one** matrix that fulfills the
# requirements exists.
#
#
#
# **Example 1:**
#
#
#
#     Input: rowSum = [3,8], colSum = [4,7]
#     Output: [[3,0],
#              [1,7]]
#     Explanation:
#     0th row: 3 + 0 = 0 == rowSum[0]
#     1st row: 1 + 7 = 8 == rowSum[1]
#     0th column: 3 + 1 = 4 == colSum[0]
#     1st column: 0 + 7 = 7 == colSum[1]
#     The row and column sums match, and all matrix elements are non-negative.
#     Another possible matrix is: [[1,2],
#                                  [3,5]]
#
#
# **Example 2:**
#
#
#
#     Input: rowSum = [5,7,10], colSum = [8,6,8]
#     Output: [[0,5,0],
#              [6,1,0],
#              [2,0,8]]
#
#
# **Example 3:**
#
#
#
#     Input: rowSum = [14,9], colSum = [6,9,8]
#     Output: [[0,9,5],
#              [6,0,3]]
#
#
# **Example 4:**
#
#
#
#     Input: rowSum = [1,0], colSum = [1]
#     Output: [[1],
#              [0]]
#
#
# **Example 5:**
#
#
#
#     Input: rowSum = [0], colSum = [0]
#     Output: [[0]]
#
#
#
#
# **Constraints:**
#
#   * `1 <= rowSum.length, colSum.length <= 500`
#   * `0 <= rowSum[i], colSum[i] <= 108`
#   * `sum(rows) == sum(columns)`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#!/usr/bin/env julia

using Cumulants
using SymmetricTensors
using CumulantsUpdates
using JLD2
using FileIO
using ArgParse
using Distributed

function comptime(dm::DataMoments{Float64}, Xup::Matrix{Float64})
  t = time_ns()
  _ = cumulantsupdate!(dm, Xup)
  Float64(time_ns()-t)/1.0e9
end

function precomp(m::Int)
  X = randn(15, 10)
  dm = DataMoments(X[1:10,:], m, 4)
  cumulantsupdate!(dm, X[10:15,:])
end

function savect(u::Vector{Int}, p::Int, n::Int, m::Int, b::Int)
  comptimes = zeros(p, length(u))
  X = randn(maximum(u)+10, n)
  precomp(m)
  addprocs(p)
  eval(Expr(:toplevel, :(@everywhere using CumulantsUpdates)))
  dm = DataMoments(X, m, b)
  for i in 1:p
    rmprocs(procs()[2:end])
    addprocs(i)
    eval(Expr(:toplevel, :(@everywhere using CumulantsUpdates)))
    println("number of workers = ", nworkers())
    for k in 1:length(u)
      Xup = randn(u[k], n)
      comptimes[i, k] = comptime(dm, Xup)
      println("u = ", u[k])
    end
  end
  filename = replace("res/$(m)_$(u)_$(n)_$(b)_nprocs.jld2", "["=>"")
  filename = replace(filename, "]"=>"")
  filename = replace(filename, " "=>"")
  compt = Dict{String, Any}("cumulants"=> comptimes)
  push!(compt, "t" => u)
  push!(compt, "n" => n)
  push!(compt, "m" => m)
  push!(compt, "x" => "number of procs")
  push!(compt, "number of procs" => [collect(1:p)...])
  push!(compt, "functions" => [["cumulants"]])
  save(filename, compt)
end

function main(args)
  s = ArgParseSettings("description")
  @add_arg_table s begin
      "--order", "-d"
        help = "d, the order of cumulant, ndims of cumulant's tensor"
        default = 4
        arg_type = Int
      "--nvar", "-n"
        default = 48
        help = "n, numbers of marginal variables"
        arg_type = Int
      "--tup", "-u"
        help = "u, numbers of data updates"
        nargs = '*'
        default = [10000, 20000]
        arg_type = Int
      "--blocksize", "-b"
        help = "the size of blocks of the block structure"
        default = 4
        arg_type = Int
      "--nprocs", "-p"
        help = "maximal number of processes"
        default = 6
        arg_type = Int
    end
  parsed_args = parse_args(s)
  m = parsed_args["order"]
  n = parsed_args["nvar"]
  u = parsed_args["tup"]
  p = parsed_args["nprocs"]
  b = parsed_args["blocksize"]
  savect(u, p, n, m, b)
end

main(ARGS)


text:
using Random: AbstractRNG
using Distributions: Univariate, Sampleable
using HDF5: h5open
using ComponentArrays: ComponentVector

"""
    DynamicalDistribution(velocity, position, size)

Sampleable that provides positions and velocities sampled from a variety of distributions.

# Arguments

**`velocity`** and **`position`** provide the velocities and positions and can be any type for
which `select_item` has been implemented.

**`size`** denotes to the size of the samples that should be produced.

# Keywords

# Example

```jldoctest; setup = :(using Random; Random.seed!(1))
using NQCDynamics.InitialConditions
using Distributions

d = DynamicalDistribution(5.0, Normal(), (1, 1))
rand(d)

# output

ComponentVector{Float64}(v = [5.0;;], r = [0.5314767537831963;;])
```
"""
struct DynamicalDistribution{V,R,S} <: NuclearDistribution
    velocity::V
    position::R
    size::NTuple{S,Int}
end

Base.size(s::DynamicalDistribution) = s.size
Base.length(s::DynamicalDistribution) = lastindex(s)

function Base.rand(rng::AbstractRNG, s::DynamicalDistribution)
    i = rand(rng, 1:lastindex(s))
    return getindex(s, i)
end

Base.firstindex(s::DynamicalDistribution) = 1
function Base.lastindex(s::DynamicalDistribution)::Int
    if s.velocity isa Vector
        return length(s.velocity)
    elseif s.position isa Vector
        return length(s.position)
    elseif s.velocity isa String
        h5open(s.velocity, "r") do fid
            return length(fid)
        end
    elseif s.position isa String
        h5open(s.position, "r") do fid
            return length(fid)
        end
    else
        return 1 end
end

function Base.getindex(s::DynamicalDistribution, i::Integer)
    return ComponentVector(
        v=select_item(s.velocity, i, s.size),
        r=select_item(s.position, i, s.size)
    )
end
Base.getindex(s::DynamicalDistribution, I) = [s[i] for i in I]

# Indexed selections
select_item(x::Vector{<:AbstractArray}, i::Integer, ::NTuple) = x[i]
select_item(x::Vector{<:Number}, i::Integer, size::NTuple) = fill(x[i], size)
function select_item(x::AbstractString, i::Integer, ::NTuple)
    h5open(x, "r") do fid
        return read(fid, string(i))
    end
end

# Sampled selection
select_item(x::Sampleable{Univariate}, ::Integer, size::NTuple) = rand(x, size)
function select_item(x::Vector{<:Sampleable{Univariate}}, ::Integer, size::NTuple)
    if length(x) == size[2]
        if length(size) == 3
            tmpsize = (size[1], 1, size[3])
        elseif length(size) == 2
            tmpsize = (size[1], 1)
        end
        a = rand.(x, Ref(tmpsize))
        return cat(a...; dims=2)
    else
        return throw(ErrorException(
            "Distribution size does not match sample size. natoms != length(distribution)"
            ))
    end
end

# Deterministic selections
select_item(x::Real, ::Integer, size::NTuple) = fill(x, size)
select_item(x::Matrix, ::Integer, ::NTuple) = x
select_item(x::AbstractArray{T,3}, ::Integer, ::NTuple) where T = x

function Base.show(io::IO, s::DynamicalDistribution)
    print(io, "DynamicalDistribution with size: ", size(s))
end

function Base.show(io::IO, ::MIME"text/plain", s::DynamicalDistribution)
    print(io, "DynamicalDistribution:\n  ", "size: ", size(s))
end

function write_hdf5(filename, data)
    h5open(filename, "w") do fid
        for (i,d) in enumerate(data)
            fid[string(i)] = d
        end
    end
end

function Base.iterate(s::DynamicalDistribution, state=1)
    state > lastindex(s) ? nothing : (s[state], state+1)
end


text:
module RulesNormalMeanVarianceVarTest

using Test
using ReactiveMP
using Random

import DomainSets
import ReactiveMP: @test_rules

@testset "rules:NormalMeanVariance:var" begin

    @testset "Belief Propagation: (m_out::PointMass, m_μ::UnivariateNormalDistributionsFamily)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:v, Marginalisation) [
            (input = (m_out = PointMass(2.0), m_μ = NormalMeanVariance(0.0, 1.0)),  output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(1.0 + x) / 2 - 4 / 2x)),
            (input = (m_out = PointMass(0.5), m_μ = NormalMeanVariance(1.0, 0.5)),  output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(0.5 + x) / 2 - 0.25 / 2x)),
            (input = (m_out = PointMass(-3.5), m_μ = NormalMeanVariance(0.5, 2.0)), output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.0 + x) / 2 - 16 / 2x)),
            (input = (m_out = PointMass(3.5), m_μ = NormalMeanVariance(-0.5, 2.0)), output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.0 + x) / 2 - 16 / 2x)),
        ]

    end

    @testset "Belief Propagation: (m_out::UnivariateNormalDistributionsFamily, m_μ::PointMass)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:v, Marginalisation) [
            (input = (m_out = NormalMeanVariance(0.0, 1.0), m_μ = PointMass(2.0)),  output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(1.0 + x) / 2 - 4 / 2x)),
            (input = (m_out = NormalMeanVariance(1.0, 0.5), m_μ = PointMass(0.5)),  output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(0.5 + x) / 2 - 0.25 / 2x)),
            (input = (m_out = NormalMeanVariance(0.5, 2.0), m_μ = PointMass(-3.5)), output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.0 + x) / 2 - 16 / 2x)),
            (input = (m_out = NormalMeanVariance(-0.5, 2.0), m_μ = PointMass(3.5)), output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.0 + x) / 2 - 16 / 2x)),
        ]

    end

    @testset "Belief Propagation: (m_out::UnivariateNormalDistributionsFamily, m_μ::UnivariateNormalDistributionsFamily)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:v, Marginalisation) [
            (input = (m_out = NormalMeanVariance(0.0, 1.0), m_μ = NormalMeanVariance(-2.0, 1.0)),  output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.0 + x) / 2 - 4 / 2x)),
            (input = (m_out = NormalMeanVariance(1.0, 0.5), m_μ = NormalMeanVariance(0.0, 2.0)),   output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.5 + x) / 2 - 1 / 2x)),
            (input = (m_out = NormalMeanVariance(0.5, 2.0), m_μ = NormalMeanVariance(1.5, 0.5)),   output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.5 + x) / 2 - 1 / 2x)),
            (input = (m_out = NormalMeanVariance(-0.5, 2.0), m_μ = NormalMeanVariance(-1.5, 0.5)), output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.5 + x) / 2 - 1 / 2x)),
        ]

    end

end

end


text:
# parameter p in geometric probability distribution of PMD
const PMDpparam = 0.3

# constant C in geometric probability distribution of PMD
const PMDconstant = 0.01

# True biological polymorphism between the ancient individual and the ref_aligned seq_aligned
const polymorphism_ancient = 0.001

# True biological polymorphism between the contaminants and the ref_aligned seq_aligned
const polymorphism_contamination = 0.001

function phred2prob(Q)
    q = Q % Int
    return 10.0^(-q / 10.0)
end

function phreds2probs(Qs)
    return phred2prob.(Qs)
end

function damage_model_modern(z)
    return 0.001
end

function damage_model_ancient(z, p = PMDpparam, C = PMDconstant)
    return Dz = p * (1 - p)^(z - 1) + C
end

function L_match(z, damage_model, quality, polymorphism)
    P_damage = damage_model(z)
    P_error = phred2prob(quality) / 3
    P_poly = polymorphism
    P_match =
        (1.0 - P_damage) * (1.0 - P_error) * (1.0 - P_poly) +
        (P_damage * P_error * (1.0 - P_poly)) +
        (P_error * P_poly * (1.0 - P_damage))
    return P_match
end

function L_mismatch(z, damage_model, quality, polymorphism)
    P_match = L_match(z, damage_model, quality, polymorphism)
    P_mismatch = 1 - P_match
    return P_mismatch
end

function compute_PMD_score(
    seq_aligned::LongDNASeq,
    ref_aligned::LongDNASeq,
    qual_aligned::Vector{UInt8},
    max_position::Int = -1,
)::Float64

    if max_position < 1
        max_position = length(seq_aligned)
    end

    L_D = 1.0
    L_M = 1.0

    z = 1
    for (s, r, q) in zip(seq_aligned, ref_aligned, qual_aligned)

        if s == DNA_N | r == DNA_N
            continue
        end

        if r == DNA_C
            if s == DNA_T
                L_D *= L_mismatch(z, damage_model_ancient, q, polymorphism_ancient)
                L_M *= L_mismatch(z, damage_model_modern, q, polymorphism_contamination)
            elseif s == DNA_C
                L_D *= L_match(z, damage_model_ancient, q, polymorphism_ancient)
                L_M *= L_match(z, damage_model_modern, q, polymorphism_contamination)
            end

        elseif r == DNA_G

            if s == DNA_A
                L_D *= L_mismatch(z, damage_model_ancient, q, polymorphism_ancient)
                L_M *= L_mismatch(z, damage_model_modern, q, polymorphism_contamination)
            elseif s == DNA_G
                L_D *= L_match(z, damage_model_ancient, q, polymorphism_ancient)
                L_M *= L_match(z, damage_model_modern, q, polymorphism_contamination)
            end
        end

        z += 1

        if z > max_position
            break
        end

    end

    PMD = log(L_D / L_M)
    return PMD

end

max_stars_repo_name: ChristianMichelsen/reference-free

max_issues_repo_name: ChristianMichelsen/reference-free

max_forks_repo_name: ChristianMichelsen/reference-free


text:
module NetSci02

using LightGraphs
using Distributions
using Loess

include("utils.jl")
include("tests.jl")

"""
    samplepercolations(measure::Function, graph::Graph, node_distributions::Function,
                       repetitions::Integer, stepsize = 0.1)

Sample `repetitions` node removals on `graph`, over deletion rates `0.0:stepsize:1.0`, calculating
`measure` each time.  The probability of deleting a nodes is given by `node_distributions(φ)`
for each φ (so `node_distributions` must have type Float -> Vector{Bernoulli}).
"""
function samplepercolations(measure::Function,
                            graph::Graph,
                            node_distributions::Function,
                            repetitions::Integer,
                            stepsize = 0.05,
                            result_type = typeof(measure(Graph(1, 0))))
    const vertices = nv(graph)
    const range = 0.0:stepsize:1.0
    curves = SharedArray(result_type, repetitions, length(range))

    pmap(1:repetitions) do r
        for (i, φ) in enumerate(range)
            victim = copy(graph)
            distributions = node_distributions(φ)

            for v = vertices:-1:1
                if rand(distributions[v]) == 1
                    rem_vertex!(victim, v)
                end
            end

            curves[r, i] = measure(victim)
        end
    end

    return curves
end

"""
    informed_probabilities(graph::Graph, weights::Vector{Float64}, infos::Vector{Symbol})

Construct a paramtrized vector of distributions over nodes which integrates knowledge about local
clustering, degree, and second-order neighbourhood, taking into account `weights`.  `infos` can
be used to choose which information to include.
"""
function informed_probabilities(graph::Graph,
                                weights::Vector{Float64},
                                infos::Vector{Symbol})::Function
    const maxdegrees = nv(graph) * (nv(graph) + 1) / 2

    # wherever meaningful, there is a normalized (*_n) and an unnormalized variant
    const information_functions = Dict(
        :degree => g -> degree(g),
        :degree_n => g -> degree(g) / maxdegrees,
        :inv_degree => g -> 1 ./ (degree(g) + eps()),
        :inv_degree_n => g -> (1 ./ (degree(g) + eps())) / maxdegrees,
        :local_clustering => g -> local_clustering_coefficient(g),
        :neighborhood2 => g -> map(v -> length(neighborhood(g, v, 2)), vertices(g)),
        :neighborhood2_n => g -> map(v -> length(neighborhood(g, v, 2)), vertices(g)) / maxdegrees,
        :neighborhood3 => g -> map(v -> length(neighborhood(g, v, 3)), vertices(g)),
        :neighborhood3_n => g -> map(v -> length(neighborhood(g, v, 3)), vertices(g)) / maxdegrees
    )

    information = Array(Float64, nv(graph), length(infos))
    for (i, info) in enumerate(infos)
        information[:, i] = information_functions[info](graph)
    end

    return makedistributions(information * weights)
end

"""
    uninformed_probabilities(graph::Graph)

Like `informed_probabilities`, but just use the same distribution based on `φ` for every node.
"""
function uninformed_probabilities(graph::Graph)
    return φ -> Bernoulli.(fill(φ, nv(graph)))
end

"""
    makedistributions(logits::Vector{Float64})::Function

Convert a vector of weight values for every node (`logits`) into a parametrized Bernoulli
distribution for each node, such that the expected fraction of nodes after using this for
percolation depends on the parameter (but the individual probabilities may differ).
"""
function makedistributions(logits::AbstractArray{Float64, 1})::Function
    probabilities = logistic.(logits)
    return function (φ)
        scaled_probabilities = probabilities .* (φ / mean(probabilities))
        Bernoulli.(clamp(scaled_probabilities, 0.0, 1.0)) # because rounding errors
    end
end

"""
    smoothcurves(curves::AbstractArray{Float64, 2})

Return the smoothed average of the given `curves` (using Loess regression).
"""
function smoothcurves{T<:Number}(curves::AbstractArray{T, 2}, use_loess = true)
    average_curve = squeeze(mapslices(mean, curves, 1), 1)

    if use_loess
        stepsize = 1 / (size(curves)[end] - 1)
        range = 0.0:stepsize:1.0
        return predict(loess(range, average_curve), range)
    else
        return average_curve
    end
end

"""
    normalizecurves(curves::AbstractArray{Float64, 2})

Normalize all curves to the range [0, 1] by their maximum.
"""
normalizecurves{T<:Number}(curves::AbstractArray{T, 2}) = curves ./ mapslices(maximum, curves, 2)

"""
    percolation_area{T<:Number}(curves::AbstractArray{T, 2}, smooth = true)

Estimate the area under the sampled percolation curves by averaging and calculating the area.
If `smooth` is true, uses Loess regression and numerical integration, otherwise simply sum up and
weight by the stepsize.
"""
function percolation_area{T<:Number}(curves::AbstractArray{T, 2}, smooth = true)
    normalized_curves = normalizecurves(curves)
    samples = size(curves)[1]
    stepsize = 1 / (size(curves)[end] - 1)

    if smooth
        # numerical integration over a local regression curve
        average_curve = squeeze(mapslices(mean, normalized_curves, 1), 1)
        smoothed_curve = loess(0.0:stepsize:1.0, average_curve)
        return quadgk(x -> predict(smoothed_curve, x), 0.0, 1.0)[1]
    else
        # just take the average of everything...
        return sum(normalized_curves) * stepsize / samples
    end
end

"""
    max_component_size(g::Graph)

Calculate the size of the larges connected component of `g`.
"""
max_component_size(g::Graph) = reduce(max, 0, length.(connected_components(g)))

"""
    percolation_loss(graph::Graph, infos::Vector{Symbol}, samples, stepsize)

Construct a loss function which samples the expected area of the percolation curve of `graph`,
using the parameter as weightings for local information.
"""
function percolation_loss(graph::Graph, infos::Vector{Symbol}, samples = 10, stepsize = 0.05,
                          measure = max_component_size)
    function loss(parameter)
        dist = informed_probabilities(graph, parameter, infos)
        curves = samplepercolations(measure, graph, dist, samples, stepsize)
        return percolation_area(curves)
    end
end

####################################################################################################
# TRAINING AND TESTING
####################################################################################################

function log_annealing(;kwargs...)
    args = Dict(kwargs)
    # push!(error_trace, args[:f_new])

    if args[:steps] % 100 == 0
        println(args[:steps], " steps, T = ", args[:T])
        println("\tParameters: ", args[:x_new])
        println("\tError: ", args[:f_new])
    end
end

function trainandtest(name::String, train_graph::Graph, test_graph::Graph,
                      results_file::String, variant::String;
                      mode = "a", parameters...)
    for (s, v) in parameters
        @eval $s = $v
    end

    println("Training on ", name, " using ", variant, "...")
    trained_param = annealing(percolation_loss(train_graph, infos, samples, stepsize),
                              fill(0.0, length(infos)),
                              init_temp, temp_decay, exploration_size,
                              (s, v) -> s > iterations;
                              debug_callback = log_annealing)
    println("Learned parameters (", name, "): ", trained_param, "\n")

    baseline_dist = uninformed_probabilities(test_graph)
    optimized_dist = informed_probabilities(test_graph, trained_param, infos)

    baseline_samples = samplepercolations(max_component_size, test_graph,
                                          baseline_dist, samples, stepsize)
    optimized_samples = samplepercolations(max_component_size, test_graph,
                                           optimized_dist, samples, stepsize)

    baseline_curve = normalizecurves(baseline_samples)
    optimized_curve = normalizecurves(optimized_samples)

    range = 0.0:stepsize:1.0
    open(results_file, mode) do f
        for (x, y, y2) in zip(range,
                              smoothcurves(baseline_curve, false),
                              smoothcurves(baseline_curve, true))
            write(f, "$name $variant baseline $x $y $y2\n")
        end

        for (x, y, y2) in zip(range,
                          smoothcurves(optimized_curve, false),
                          smoothcurves(optimized_curve, true))
            write(f, "$name $variant optimized $x $y $y2\n")
        end
    end
end

function main()
    fb1 = readnetwork("../data/facebook1.txt")
    fb2 = readnetwork("../data/facebook2.txt")
    austin = readnetwork("../data/austin.txt")
    philadelphia = readnetwork("../data/philadelphia.txt")
    random1 = Graph(2000, 5000)
    random2 = Graph(2000, 5000)

    examples = ["facebook" => (fb1, fb2),
                "traffic" => (austin, philadelphia),
                "random" => (random1, random2)]

    info_variants = Dict("d" => [:degree],
                         "i" => [:inv_degree],
                         "l" => [:local_clustering],
                         "n" => [:neighborhood2],
                         "dl" => [:degree, :local_clustering],
                         "dln" => [:degree, :local_clustering, :neighborhood2])

    info_variants_n = Dict("i" => [:inv_degree_n],
                           "d" => [:degree_n],
                           "l" => [:local_clustering],
                           "n" => [:neighborhood2_n],
                           "dl" => [:degree_n, :local_clustering],
                           "dln" => [:degree_n, :local_clustering, :neighborhood2_n])

    training_parameters = Dict(
        :stepsize => 0.05,
        :samples => 10,
        :init_temp => 10.0,
        :temp_decay => 0.9,
        :exploration_size => 0.05,
        :iterations => 400)

    for (variant, infos) in info_variants
        for (name, (train, test)) in examples
            trainandtest(name, train, test,
                         "../evaluation/results1.txt", variant;
                         mode = "a", infos = infos, training_parameters...)
        end
    end

    # this turns out to be worse...
    # for (variant, infos) in info_variants_n
    #     for (name, (train, test)) in examples
    #         trainandtest(name, train, test,
    #                      "../evaluation/results2.txt", variant;
    #                      mode = "a", infos = infos, training_parameters...)
    #     end
    # end
end

end


text:
using BenchmarkTools
using ConvexBodyProximityQueries
using CurveProximityQueries
using Plots
using StaticArrays
import Random: seed!

pgfplots()

function randomtraj(num)
    Bs = []
    xax = collect(0.0:0.2:1.0)'
    i1 = @SVector [-0.2, 0.5]
    i2 = @SVector [ 0.0, 0.5]
    e1 = @SVector [ 1.0, 0.5]
    e2 = @SVector [ 1.2, 0.5]
    for i = 1:num
        m1 = @SVector rand(2)
        m2 = @SVector rand(2)
        b = Bernstein(i1, i2, m1, m2, e1, e2)
        push!(Bs, b)
    end
    return Bs
end

import CurveProximityQueries:tolerance_verification
function tolerance_verification(obs::Array, traj::Curve, Δ::Real)
    ret = true
    for o in obs
        ret = ret && tolerance_verification(o, traj, Δ)
    end
    return ret
end

import CurveProximityQueries:collision_detection
function collision_detection(obs::Array, traj::Curve)
    ret = false
    for o in obs
        ret = ret || collision_detection(o, traj)
    end
    return ret
end

function counttv(obs, trajs)
    cnt = 0.0
    for t in trajs
        cnt += tolerance_verification(obs, t, 0.03)
    end
    return cnt
end

seed!(22)
obs = [randpoly([0.35, 0.6]; n=20, scale=0.2),
       randpoly([0.65, 0.4]; n=20, scale=0.2)]
trajs = randomtraj(1000)

#  @btime counttv($obs, $trajs)

plot()
plot!.(obs,
       fillalpha=1.0,
       fillcolor=Gray(0.7),
      )
cnt = zeros(3)
for t in trajs
    if tolerance_verification(obs, t, 0.03)
        cnt[1] += 1
        plot!(t,
              linecolor = RGBA(77/255,141/225,205/255,1.0),
              linewidth = 1.2)
    elseif !collision_detection(obs, t)
        cnt[2] += 1
        plot!(t,
              linecolor = RGBA(231/255,95/225,97/255,0.3),
              linewidth = 1.2)
    else
        cnt[3] += 1
        plot!(t,
              linecolor = GrayA(0.5, 0.05),
              linewidth = 0.5)
    end
end
plot!(xlim=(-0.3, 1.3), ylim=(0.0, 1.0))
savefig("out/fig12.pdf")


text:
function is_between(value::String, low::Int, high::Int)::Bool
    return low <= parse(Int, value) <= high
end

function is_valid_height(value::String)::Bool
    if endswith(value, "cm")
        return is_between(value[1:end - 2], 150, 193)
    elseif endswith(value, "in")
        return is_between(value[1:end - 2], 59, 76)
    else
        return false
    end
end

function check_constraints(fields::Dict{String,String})::Bool
    required_fields = Set{String}(("byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"))
    if !(issubset(required_fields, keys(fields)))
        return false
    elseif !is_between(fields["byr"], 1920, 2002)
        return false
    elseif !is_between(fields["iyr"], 2010, 2020)
        return false
    elseif !is_between(fields["eyr"], 2020, 2030)
        return false
    elseif !is_valid_height(fields["hgt"])
        return false
    elseif !occursin(r"^#[0-9a-f]{6}$", fields["hcl"])
        return false
    elseif !(fields["ecl"] in Set{String}(("amb", "blu", "brn", "gry", "grn", "hzl", "oth")))
        return false
    elseif !occursin(r"^[0-9]{9}$", fields["pid"])
        return false
    end
    return true
end

function run(s::String)::Int
    # Your code here
    passports::Array{String} = split(s, "\n\n")
    counter::Int = 0
    for passport_str in passports
        tags::Array{String} = split(replace(passport_str, '\n' => ' '), ' ')
        fields = Dict{String,String}(map(x -> split(x, ":"), tags))
        if check_constraints(fields)
            counter += 1
        end
    end
    return counter
end

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

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()


text:
using Compat, FastTransforms, LowRankApprox
using Compat.Test
import FastTransforms: Cnλ, Λ, lambertw, Cnαβ, Anαβ, pochhammer
import FastTransforms: clenshawcurtisnodes, clenshawcurtisweights, fejernodes1, fejerweights1, fejernodes2, fejerweights2
import FastTransforms: chebyshevmoments1, chebyshevmoments2, chebyshevjacobimoments1, chebyshevjacobimoments2, chebyshevlogmoments1, chebyshevlogmoments2
import Compat: range

@testset "Special functions" begin
    n = 0:1000_000
    λ = 0.123
    @time FastTransforms.Cnλ(n,λ)
    @time FastTransforms.Cnλ(n,λ)

    x =range(0,stop=20,length=81)
    @test norm((Λ.(x) .- Λ.(big.(x)))./Λ.(x),Inf) < 2eps()
    @test norm((lambertw.(x) .- lambertw.(big.(x)))./max.(lambertw.(x),1), Inf) < 2eps()

    x = 0:0.5:10_000
    λ₁,λ₂ = 0.125,0.875
    @test norm((Λ.(x,λ₁,λ₂) .- Λ.(big.(x),big(λ₁),big(λ₂)))./Λ.(big.(x),big(λ₁),big(λ₂)),Inf) < 4eps()
    λ₁,λ₂ = 1//3,2//3
    @test norm((Λ.(x,Float64(λ₁),Float64(λ₂)) .- Λ.(big.(x),big(λ₁),big(λ₂))) ./ Λ.(big.(x),big(λ₁),big(λ₂)),Inf) < 4eps()

    n = 0:1000
    α = 0.125
    β = 0.375

    @time FastTransforms.Cnαβ.(n,α,β);
    @test norm(FastTransforms.Cnαβ.(n,α,β) ./ FastTransforms.Cnαβ.(n,big(α),big(β)) .- 1,Inf) < 3eps()

    @time FastTransforms.Anαβ.(n,α,β);
    @test norm(FastTransforms.Anαβ.(n,α,β) ./ FastTransforms.Anαβ.(n,big(α),big(β)) .- 1,Inf) < 4eps()

    @test pochhammer(2,1) == pochhammer(2,1.0) == pochhammer(2.0,1) == 2
    @test pochhammer(1.1,2.2) ≈ gamma(3.3)/gamma(1.1)
    @test pochhammer(-2,1) == pochhammer(-2,1.0) == pochhammer(-2.0,1) == -2
end

@testset "Fejer and Clenshaw--Curtis quadrature" begin
    N = 20
    f(x) = exp(x)

    x = clenshawcurtisnodes(Float64, N)
    μ = chebyshevmoments1(Float64, N)
    w = clenshawcurtisweights(μ)
    @test norm(dot(f.(x),w)-2sinh(1)) ≤ 4eps()

    μ = chebyshevjacobimoments1(Float64, N, 0.25, 0.35)
    w = clenshawcurtisweights(μ)
    @test norm(dot(f.(x),w)-2.0351088204147243) ≤ 4eps()

    μ = chebyshevlogmoments1(Float64, N)
    w = clenshawcurtisweights(μ)
    @test norm(sum(w./(x .- 3)) - π^2/12) ≤ 4eps()

    x = fejernodes1(Float64, N)
    μ = chebyshevmoments1(Float64, N)
    w = fejerweights1(μ)
    @test norm(dot(f.(x),w)-2sinh(1)) ≤ 4eps()

    μ = chebyshevjacobimoments1(Float64, N, 0.25, 0.35)
    w = fejerweights1(μ)
    @test norm(dot(f.(x),w)-2.0351088204147243) ≤ 4eps()

    μ = chebyshevlogmoments1(Float64, N)
    w = fejerweights1(μ)
    @test norm(sum(w./(x .- 3)) - π^2/12) ≤ 4eps()

    x = fejernodes2(Float64, N)
    μ = chebyshevmoments2(Float64, N)
    w = fejerweights2(μ)
    @test norm(dot(f.(x),w)-2sinh(1)) ≤ 4eps()

    μ = chebyshevjacobimoments2(Float64, N, 0.25, 0.35)
    w = fejerweights2(μ)
    @test norm(dot(f.(x),w)-2.0351088204147243) ≤ 4eps()

    μ = chebyshevlogmoments2(Float64, N)
    w = fejerweights2(μ)
    @test norm(sum(w./(x .- 3)) - π^2/12) ≤ 4eps()
end

if VERSION < v"0.7-"
    @testset "Allocation-free ID matrix-vector products" begin
        for T in (Float64, ComplexF64)
            r = rand(T)
            A = idfact([r/(i+j-1) for i in 1:200, j = 1:50])
            P = A[:P]
            k, n = size(A)

            x = rand(T, n)
            y = zeros(T, k)

            @test norm(FastTransforms.mul!(y, A, P, x, 1, 1) - A*x) < 10eps()*norm(A*x)

            x = rand(T, k)
            y = zeros(T, n)

            @test norm(FastTransforms.At_mul_B!(y, A, P, x, 1, 1) - At_mul_B(A,x), Inf) < 10eps()

            fill!(y, zero(T))
            @test norm(FastTransforms.Ac_mul_B!(y, A, P, x, 1, 1) - A'x, Inf) < 10eps()
        end
    end
else
    @testset "Allocation-free ID matrix-vector products" begin
        for T in (Float64, ComplexF64)
            r = rand(T)
            A = idfact([r/(i+j-1) for i in 1:200, j = 1:50])
            P = A[:P]
            k, n = size(A)

            x = rand(T, n)
            y = zeros(T, k)

            @test norm(FastTransforms.mul!(y, A, P, x, 1, 1) - A*x) < 10eps()*norm(A*x)

            x = rand(T, k)
            y = zeros(T, n)

            @test norm(FastTransforms.At_mul_B!(y, A, P, x, 1, 1) - transpose(A)*x, Inf) < 10eps()

            fill!(y, zero(T))
            @test norm(FastTransforms.Ac_mul_B!(y, A, P, x, 1, 1) - A'x, Inf) < 10eps()
        end
    end
end


text:
# Calculate structure factor
# special case: only for 1 species with Z = 1
function structure_factor( Xpos::Array{Float64,2}, G::Array{Float64,2} )
    Ng = size(G)[2]
    Na = size(Xpos)[2]
    Sf = zeros(ComplexF64,Ng)
    for ia = 1:Na
        for ig = 1:Ng
            GX = Xpos[1,ia]*G[1,ig] + Xpos[2,ia]*G[2,ig] + Xpos[3,ia]*G[3,ig]
            Sf[ig] = Sf[ig] + cos(GX) - im*sin(GX)
        end
    end
    return Sf
end


text:
using ScrapeSEC
using Documenter

Documenter.makedocs(
    modules = [ScrapeSEC],
    sitename = "ScrapeSEC.jl",
)