text:
## ExoplanetsSysSim/examples/hsu_etal_2018/christiansen_func.jl
## (c) 2018 Danley C. Hsu & Eric B. Ford
# Collection of functions specific to estimating Q1-Q16 FGK
#   planet candidate occurrence rates over a 2D period-radius grid

using ExoplanetsSysSim
using StatsFuns
using JLD
using CSV
using DataFrames
using Distributions

## simulation_parameters
macro isdefinedlocal(var)
    quote
        try
            $(esc(var))
            true
        catch err
            isa(err, UndefVarError) ? false : rethrow(err)
        end
    end
end

function setup_sim_param_christiansen(args::Vector{String} = Array{String}(0) )   # allow this to take a list of parameter (e.g., from command line)
    sim_param = ExoplanetsSysSim.SimParam()

    add_param_fixed(sim_param,"max_tranets_in_sys",7)
    add_param_fixed(sim_param,"generate_star",ExoplanetsSysSim.generate_star_dumb)
    add_param_fixed(sim_param,"generate_planetary_system", ExoplanetsSysSim.generate_planetary_system_uncorrelated_incl)
    add_param_fixed(sim_param,"generate_kepler_target",ExoplanetsSysSim.generate_kepler_target_from_table)
    add_param_fixed(sim_param,"star_table_setup",setup_star_table_christiansen)
    add_param_fixed(sim_param,"stellar_catalog","q1_q16_christiansen.jld")
    add_param_fixed(sim_param,"generate_num_planets",generate_num_planets_christiansen)
    add_param_fixed(sim_param,"generate_planet_mass_from_radius",ExoplanetsSysSim.generate_planet_mass_from_radius_powerlaw)
    add_param_fixed(sim_param,"mr_power_index",2.0)
    add_param_fixed(sim_param,"mr_const",1.0)
    add_param_fixed(sim_param,"generate_period_and_sizes", generate_period_and_sizes_christiansen)
    #p_lim_arr_num = [0.5, 1.25, 2.5, 5., 10., 20., 40., 80., 160., 320.]
    #r_lim_arr_num = [0.5, 0.75, 1., 1.25, 1.5, 1.75, 2., 2.5, 3., 4., 6., 8., 12., 16.]
    #p_dim = length(p_lim_arr_num)-1
    #r_dim = length(r_lim_arr_num)-1
    #rate_tab_init = reshape(fill(1.0, p_dim*r_dim)*0.01,(r_dim,p_dim))
    #add_param_fixed(sim_param, "p_lim_arr", p_lim_arr_num)
    #add_param_fixed(sim_param, "r_lim_arr", r_lim_arr_num*ExoplanetsSysSim.earth_radius)
    #add_param_active(sim_param,"obs_par", rate_tab_init)
    add_param_fixed(sim_param,"generate_e_omega",ExoplanetsSysSim.generate_e_omega_rayleigh)
    add_param_fixed(sim_param,"sigma_hk",0.03)
    add_param_fixed(sim_param,"sigma_incl",2.0)   # degrees
    add_param_fixed(sim_param,"calc_target_obs_sky_ave",ExoplanetsSysSim.calc_target_obs_sky_ave)
    add_param_fixed(sim_param,"calc_target_obs_single_obs",ExoplanetsSysSim.calc_target_obs_single_obs)
    add_param_fixed(sim_param,"transit_noise_model",ExoplanetsSysSim.transit_noise_model_diagonal)
    return sim_param
end

function set_test_param(sim_param_closure::SimParam)
    @eval(include(joinpath(pwd(),"param.in")))

    if @isdefinedlocal(stellar_catalog)
        @assert (typeof(stellar_catalog) == String)
        add_param_fixed(sim_param_closure,"stellar_catalog",stellar_catalog)
    end
    if @isdefinedlocal(koi_catalog)
        @assert (typeof(koi_catalog) == String)
        add_param_fixed(sim_param_closure,"koi_catalog",koi_catalog)
    end

    if @isdefinedlocal(num_targ_sim)
        @assert (typeof(num_targ_sim) == Int)
        add_param_fixed(sim_param_closure,"num_targets_sim_pass_one",num_targ_sim)
    end

    @assert (typeof(p_bin_lim) == Array{Float64,1})
    add_param_fixed(sim_param_closure, "p_lim_arr", p_bin_lim)

    @assert (typeof(r_bin_lim) == Array{Float64,1})
    add_param_fixed(sim_param_closure, "r_lim_arr", r_bin_lim*ExoplanetsSysSim.earth_radius)

    p_dim = length(get_any(sim_param_closure, "p_lim_arr", Array{Float64,1}))-1
    r_dim = length(get_any(sim_param_closure, "r_lim_arr", Array{Float64,1}))-1
    n_bin = p_dim*r_dim

    if @isdefinedlocal(rate_init)
        if typeof(rate_init) <: Real
            @assert (rate_init >= 0.0)
            rate_init = fill(rate_init, n_bin)
        end

        @assert (ndims(rate_init) <= 2)
        if ndims(rate_init) == 1
            @assert (length(rate_init) == n_bin)
            rate_tab_init = reshape(rate_init*0.01, (r_dim, p_dim))
        else
            @assert (size(rate_init) == (r_dim, p_dim))
            rate_tab_init = rate_init*0.01
        end
        add_param_active(sim_param_closure, "obs_par", rate_tab_init)
    else
        rate_init = fill(1.0, n_bin)
        rate_tab_init = reshape(rate_init*0.01, (r_dim, p_dim))
        add_param_active(sim_param_closure, "obs_par", rate_tab_init)
    end

    return sim_param_closure
end

## planetary_system
function generate_num_planets_christiansen(s::Star, sim_param::SimParam)
  const max_tranets_in_sys::Int64 = get_int(sim_param,"max_tranets_in_sys")
  rate_tab::Array{Float64,2} = get_any(sim_param, "obs_par", Array{Float64,2})
  lambda = sum_kbn(rate_tab)
  #println("# lambda= ", lambda)
  ExoplanetsSysSim.generate_num_planets_poisson(lambda,max_tranets_in_sys)
end

function generate_period_and_sizes_christiansen(s::Star, sim_param::SimParam; num_pl::Integer = 1)
  rate_tab::Array{Float64,2} = get_any(sim_param, "obs_par", Array{Float64,2})

  limitP::Array{Float64,1} = get_any(sim_param, "p_lim_arr", Array{Float64,1})
  limitRp::Array{Float64,1} = get_any(sim_param, "r_lim_arr", Array{Float64,1})

  @assert ((length(limitP)-1) == size(rate_tab, 2))
  @assert ((length(limitRp)-1) == size(rate_tab, 1))
  Plist = []
  Rplist = []
  rate_tab_1d = reshape(rate_tab,length(rate_tab))
  #logmaxcuml = logsumexp(rate_tab_1d)
  #cuml = cumsum_kbn(exp(rate_tab_1d-logmaxcuml))
  maxcuml = sum(rate_tab_1d)
  cuml = cumsum_kbn(rate_tab_1d/maxcuml)

  for n in 1:num_pl
    rollp = Base.rand()
    idx = findfirst(x -> x > rollp, cuml)
    i_idx = (idx-1)%size(rate_tab,1)+1
    j_idx = floor(Int64,(idx-1)//size(rate_tab,1))+1
    ### TODO: Keep uniform log sampling here?
    Rp = exp(Base.rand()*(log(limitRp[i_idx+1])-log(limitRp[i_idx]))+log(limitRp[i_idx]))
    P = exp(Base.rand()*(log(limitP[j_idx+1])-log(limitP[j_idx]))+log(limitP[j_idx]))
    push!(Plist, P)
    push!(Rplist, Rp)
  end
  return Plist, Rplist
end

## stellar_table
function setup_christiansen(sim_param::SimParam; force_reread::Bool = false)
  #global df
  if haskey(sim_param,"read_stellar_catalog") && !force_reread
     return ExoplanetsSysSim.StellarTable.df
     #return data
  end
  stellar_catalog_filename = convert(String,joinpath(Pkg.dir("ExoplanetsSysSim"), "data", convert(String,get(sim_param,"stellar_catalog","q1_q17_dr25_stellar.csv")) ) )
  df = setup_christiansen(stellar_catalog_filename)
  add_param_fixed(sim_param,"read_stellar_catalog",true)
  add_param_fixed(sim_param,"num_kepler_targets",StellarTable.num_usable_in_star_table())
  if !haskey(sim_param.param,"num_targets_sim_pass_one")
      add_param_fixed(sim_param_closure,"num_targets_sim_pass_one", StellarTable.num_usable_in_star_table())
  end
  StellarTable.set_star_table(df)
  return df
end

function setup_christiansen(filename::String; force_reread::Bool = false)
  #df = ExoplanetsSysSim.StellarTable.df
  if ismatch(r".jld$",filename)
  try
    data = load(filename)
    df::DataFrame = data["stellar_catalog"]
  catch
    error(string("# Failed to read stellar catalog >",filename,"< in jld format."))
  end
  else
  try
    #df = CSV.read(filename,nullable=true)
    df = CSV.read(filename, allowmissing=:all)
  catch
    error(string("# Failed to read stellar catalog >",filename,"< in ascii format."))
  end

  has_mass = .! (ismissing.(df[:mass]) .| ismissing.(df[:mass_err1]) .| ismissing.(df[:mass_err2]))
  has_radius = .! (ismissing.(df[:radius]) .| ismissing.(df[:radius_err1]) .| ismissing.(df[:radius_err2]))
  has_dens = .! (ismissing.(df[:dens]) .| ismissing.(df[:dens_err1]) .| ismissing.(df[:dens_err2]))
  has_rest = .! (ismissing.(df[:rrmscdpp04p5]) .| ismissing.(df[:dataspan]) .| ismissing.(df[:dutycycle]))
  in_Q1Q12 = []
  for x in df[:st_quarters]
    subx = string(x)
    subx = ("0"^(17-length(subx)))*subx
    indQ = search(subx, '1')
    if ((indQ < 1) | (indQ > 12))
      push!(in_Q1Q12, false)
    else
      push!(in_Q1Q12, true)
    end
  end
  is_FGK = []
  for x in 1:length(df[:teff])
    if ((df[x,:teff] > 4000.0) & (df[x,:teff] < 7000.0) & (df[x,:logg] > 4.0))
      push!(is_FGK, true)
    else
      push!(is_FGK, false)
    end
  end
  is_usable = has_radius .& is_FGK .& has_mass .& has_rest .& has_dens
  if contains(filename,"q1_q16_stellar.csv")
    is_usable = is_usable .& in_Q1Q12
  end
  # See options at: http://exoplanetarchive.ipac.caltech.edu/docs/API_keplerstellar_columns.html
  symbols_to_keep = [ :kepid, :mass, :mass_err1, :mass_err2, :radius, :radius_err1, :radius_err2, :dens, :dens_err1, :dens_err2, :rrmscdpp01p5, :rrmscdpp02p0, :rrmscdpp02p5, :rrmscdpp03p0, :rrmscdpp03p5, :rrmscdpp04p5, :rrmscdpp05p0, :rrmscdpp06p0, :rrmscdpp07p5, :rrmscdpp09p0, :rrmscdpp10p5, :rrmscdpp12p0, :rrmscdpp12p5, :rrmscdpp15p0, :dataspan, :dutycycle ]

  delete!(df, [~(x in symbols_to_keep) for x in names(df)])    # delete columns that we won't be using anyway
  usable = find(is_usable)
  df = df[usable, symbols_to_keep]
  tmp_df = DataFrame()
  for col in names(df)
      tmp_df[col] = collect(skipmissing(df[col]))
  end
  df = tmp_df
  StellarTable.set_star_table(df)
  end
#    println("# Removing stars observed <5 quarters (not done in original Hsu et al paper).")
#    df[:wf_id] = map(x->ExoplanetsSysSim.WindowFunction.get_window_function_id(x,use_default_for_unknown=false),df[:kepid])
#    obs_5q = df[:wf_id].!=-1
#    df = df[obs_5q,keys(df.colindex)]
#    StellarTable.set_star_table(df)
  return df
end

setup_star_table_christiansen(sim_param::SimParam; force_reread::Bool = false) = setup_christiansen(sim_param, force_reread=force_reread)
setup_star_table_christiansen(filename::String; force_reread::Bool = false) = setup_christiansen(filename, force_reread=force_reread)

## summary_statistics
function calc_summary_stats_sim_pass_one_binned_rates(cat_obs::KeplerObsCatalog, cat_phys::KeplerPhysicalCatalog, param::SimParam )      # Version for simulated data, since includes cat_phys
  ssd = Dict{String,Any}()
  cache = Dict{String,Any}()

  max_tranets_in_sys = get_int(param,"max_tranets_in_sys")    # Demo that simulation parameters can specify how to evalute models, too
  @assert max_tranets_in_sys >= 1
  idx_tranets = find(x::KeplerTargetObs-> length(x.obs) > 0, cat_obs.target)::Array{Int64,1}             # Find indices of systems with at least 1 tranet = potentially detectable transiting planet

  # Count total number of tranets and compile indices for N-tranet systems
  num_tranets = 0
  idx_n_tranets = Vector{Int64}[ Int64[] for m = 1:max_tranets_in_sys]
  for n in 1:max_tranets_in_sys-1
    idx_n_tranets[n] = find(x::KeplerTargetObs-> length(x.obs) == n, cat_obs.target[idx_tranets] )
    num_tranets += n*length(idx_n_tranets[n])
  end
  idx_n_tranets[max_tranets_in_sys] = find(x::KeplerTargetObs-> length(x.obs) >= max_tranets_in_sys, cat_obs.target[idx_tranets] )

  num_tranets += max_tranets_in_sys*length(idx_n_tranets[max_tranets_in_sys])  # WARNING: this means we need to ignore planets w/ indices > max_tranets_in_sys
  num_tranets  = convert(Int64,num_tranets)            # TODO OPT: Figure out why isn't this already an Int.  I may be doing something that prevents some optimizations

  cache["num_tranets"] = num_tranets
  cache["idx_tranets"] = idx_tranets                                   # We can save lists of indices to summary stats for pass 2, even though we won't use these for computing a distance or probability
  #cache["idx_n_tranets"] = idx_n_tranets

  expected_num_detect = 0.0
  expected_num_sys_n_tranets = zeros(max_tranets_in_sys)
  period_list = zeros(num_tranets)
  weight_list = zeros(num_tranets)
  radius_list = zeros(num_tranets)

  n = 1    # tranet id
  for i in idx_tranets
    for j in 1:num_planets(cat_obs.target[i])
      p_tr_and_det = ExoplanetsSysSim.prob_detect(cat_obs.target[i].prob_detect,j)
      expected_num_detect += p_tr_and_det
      (s,p) = cat_obs.target[i].phys_id[j]

      period_list[n] = cat_phys.target[i].sys[s].orbit[p].P
      weight_list[n] = p_tr_and_det
      radius_list[n] = cat_phys.target[i].sys[s].planet[p].radius
      n = n+1
    end
    for k in 1:max_tranets_in_sys
      expected_num_sys_n_tranets[k] += ExoplanetsSysSim.prob_detect_n_planets(cat_obs.target[i].prob_detect,k)
    end
  end
  ssd["expected planets detected"] = expected_num_detect
  ssd["num_sys_tranets"] = expected_num_sys_n_tranets
  ssd["num targets"] = get_int(param,"num_targets_sim_pass_one")
  #println("expected planets = ",expected_num_detect,", num_sys_tranets = ",expected_num_sys_n_tranets,", num targets = ",ssd["num targets"])

  limitP::Array{Float64,1} = get_any(param, "p_lim_arr", Array{Float64,1})
  limitRp::Array{Float64,1} = get_any(param, "r_lim_arr", Array{Float64,1})

  np_bin = zeros((length(limitP)-1) * (length(limitRp)-1))
  np_bin_idx = 1
  for i in 1:(length(limitP)-1)
    P_match = find(x -> ((x > limitP[i]) && (x < limitP[i+1])), period_list)
    for j in 1:(length(limitRp)-1)
      R_match = find(x -> ((x > limitRp[j]) && (x < limitRp[j+1])), radius_list)

      bin_match = intersect(P_match, R_match)

      np_bin[np_bin_idx] = sum(weight_list[bin_match])
      np_bin_idx += 1
    end
  end

  #ssd["expected planets detected"] = sum(np_bin)
  ssd["expected planets table"] = np_bin

  return CatalogSummaryStatistics(ssd, cache)
end

function calc_summary_stats_obs_binned_rates(cat_obs::KeplerObsCatalog, param::SimParam; trueobs_cat::Bool = false)
  ssd = Dict{String,Any}()
  cache = Dict{String,Any}()

  if !trueobs_cat
    ssd["num targets"] = get_int(param,"num_targets_sim_pass_one")
  else
    ssd["num targets"] = get_int(param,"num_kepler_targets")
  end

  max_tranets_in_sys = get_int(param,"max_tranets_in_sys")    # Demo that simulation parameters can specify how to evalute models, too
  @assert max_tranets_in_sys >= 1
  idx_tranets = find(x::KeplerTargetObs-> length(x.obs) > 0, cat_obs.target)::Array{Int64,1}             # Find indices of systems with at least 1 tranet = potentially detectable transiting planet

  # Count total number of tranets and compile indices for N-tranet systems
  num_tranets = 0
  idx_n_tranets = Vector{Int64}[ Int64[] for m = 1:max_tranets_in_sys]
  for n in 1:max_tranets_in_sys-1
    idx_n_tranets[n] = find(x::KeplerTargetObs-> length(x.obs) == n, cat_obs.target[idx_tranets] )
    num_tranets += n*length(idx_n_tranets[n])
  end
  idx_n_tranets[max_tranets_in_sys] = find(x::KeplerTargetObs-> length(x.obs) >= max_tranets_in_sys, cat_obs.target[idx_tranets] )

  num_tranets += max_tranets_in_sys*length(idx_n_tranets[max_tranets_in_sys])  # WARNING: this means we need to ignore planets w/ indices > max_tranets_in_sys
  if ( length( find(x::KeplerTargetObs-> length(x.obs) > max_tranets_in_sys, cat_obs.target[idx_tranets] ) ) > 0)   # Make sure max_tranets_in_sys is at least big enough for observed systems
    warn("Observational data has more transiting planets in one systems than max_tranets_in_sys allows.")
  end
  num_tranets  = convert(Int64,num_tranets)            # TODO OPT: Figure out why isn't this already an Int.  I may be doing something that prevents some optimizations

  num_sys_tranets = zeros(max_tranets_in_sys)                           # Since observed data, don't need to calculate probabilities.
  for n in 1:max_tranets_in_sys                                         # Make histogram of N-tranet systems
    num_sys_tranets[n] = length(idx_n_tranets[n])
  end
  ssd["num_sys_tranets"] = num_sys_tranets
  ssd["planets detected"] = num_tranets

  period_list = zeros(num_tranets)
  weight_list = zeros(num_tranets)
  radius_list = zeros(num_tranets)

  n = 1    # tranet id
  for i in idx_tranets
    for j in 1:num_planets(cat_obs.target[i])
      period_list[n] = cat_obs.target[i].obs[j].period
      weight_list[n] = 1.0
      radius_list[n] = sqrt(cat_obs.target[i].obs[j].depth)*cat_obs.target[i].star.radius
      n = n+1
    end
  end

  limitP::Array{Float64,1} = get_any(param, "p_lim_arr", Array{Float64,1})
  limitRp::Array{Float64,1} = get_any(param, "r_lim_arr", Array{Float64,1})

  np_bin = zeros((length(limitP)-1) * (length(limitRp)-1))
  np_bin_idx = 1
  for i in 1:(length(limitP)-1)
    P_match = find(x -> ((x > limitP[i]) && (x < limitP[i+1])), period_list)
    for j in 1:(length(limitRp)-1)
      R_match = find(x -> ((x > limitRp[j]) && (x < limitRp[j+1])), radius_list)

      bin_match = intersect(P_match, R_match)

      np_bin[np_bin_idx] = sum(weight_list[bin_match])
      np_bin_idx += 1
    end
  end

  #ssd["planets detected"] = sum(np_bin)
  ssd["planets table"] = np_bin

  return CatalogSummaryStatistics(ssd, cache)
end

## abc_distance
function calc_distance_vector_binned(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, pass::Int64, sim_param::SimParam ; verbose::Bool = false)
  d = Array{Float64}(0)
  if pass == 1
    if verbose
      println("# Summary 1, pass 1: ",summary1)
      println("# Summary 2, pass 1: ",summary2)
    end
    d = zeros(3)
    # Since observed and simulated catalogs can have different summary statistics for the number of planets, prefer detections if avaliable (e.g., after pass2), otherwise use expected (e.g., from pass 1)
    #np1 = haskey(summary1.stat,"planets detected") ? summary1.stat["planets detected"] : summary1.stat["expected planets detected"]
    #np2 = haskey(summary2.stat,"planets detected") ? summary2.stat["planets detected"] : summary2.stat["expected planets detected"]
    #d[1] = dist_L1_abs(np1/summary1.stat["num targets"],np2/summary2.stat["num targets"])    #  Normalize so different statistics weighted appropriately and not dominated by this one
    #println("np1 = ",np1,", np2 = ",np2)
    #println("np1 (normalized) = ",np1/summary1.stat["num targets"],", np2 (normalized) = ",np2/summary2.stat["num targets"],", d[1] = ",d[1])

    np1 = haskey(summary1.stat,"planets table") ? summary1.stat["planets table"] : summary1.stat["expected planets table"]
    np2 = haskey(summary2.stat,"planets table") ? summary2.stat["planets table"] : summary2.stat["expected planets table"]

    np_bin = zeros(length(np1))
    for n in 1:length(np1)
        #np_bin[n] = dist_L1_abs(np1[n]/summary1.stat["num targets"], np2[n]/summary2.stat["num targets"])
        np_bin[n] = dist_L2_abs(np1[n]/summary1.stat["num targets"], np2[n]/summary2.stat["num targets"])

      #println("True # [Bin ", n,"] = ",np1[n],", Expected # [Bin ", n,"] = ",np2[n])
    end
      #d[1] = maximum(np_bin)
      d[1] = sum(np_bin)
    else
    println("# calc_distance_vector_demo doesn't know what to do for pass= ", pass)
  end
  return d
end

## eval_model
function test_christiansen()
  global sim_param_closure = setup_sim_param_christiansen()
  cat_phys = generate_kepler_physical_catalog(sim_param_closure)
  cat_obs = observe_kepler_targets_single_obs(cat_phys,sim_param_closure)
  global summary_stat_ref_closure = calc_summary_stats_obs_demo(cat_obs,sim_param_closure)
  global cat_phys_try_closure  = generate_christiansen_catalog(sim_param_closure)
  global cat_obs_try_closure  = observe_kepler_targets_sky_avg(cat_phys_try_closure,sim_param_closure)
  global summary_stat_try_closure  = calc_summary_stats_sim_pass_one_demo(cat_obs_try_closure,cat_phys_try_closure,sim_param_closure)
  summary_stat_try_closure   = calc_summary_stats_sim_pass_two_demo(cat_obs_try_closure,cat_phys_try_closure,summary_stat_try_closure,sim_param_closure)
  param_guess = make_vector_of_sim_param(sim_xparam_closure)
  evaluate_model_scalar_ret( param_guess)
end

## inverse_detection & simple bayesian
function inv_det(cat_obs::KeplerObsCatalog, param::SimParam)
    num_targ = ExoplanetsSysSim.StellarTable.num_usable_in_star_table()

    limitP::Array{Float64,1} = get_any(param, "p_lim_arr", Array{Float64,1})
    limitRp::Array{Float64,1} = get_any(param, "r_lim_arr", Array{Float64,1})

    println("------------------------------")
    cnt_bin, np_bin = cnt_np_bin(cat_obs, param)
    println("------------------------------")

    println("Inverse Detection Rates:")
    for i in 1:(length(limitP)-1)
        for j in 1:(length(limitRp)-1)
            rate_f = np_bin[(i-1)*(length(limitRp)-1) + j]/num_targ*100.
            if cnt_bin[(i-1)*(length(limitRp)-1) + j] > 0.
                println(rate_f,
                        " +/- ", rate_f/sqrt(cnt_bin[(i-1)*(length(limitRp)-1) + j]), " %")
            else
                println(rate_f,
                        " +/- N/A %")
            end
        end
    end
    println()
end

function simp_bayes(cat_obs::KeplerObsCatalog, param::SimParam)
    num_targ = ExoplanetsSysSim.StellarTable.num_usable_in_star_table()

    limitP::Array{Float64,1} = get_any(param, "p_lim_arr", Array{Float64,1})
    limitRp::Array{Float64,1} = get_any(param, "r_lim_arr", Array{Float64,1})

    println("------------------------------")
    cnt_bin, np_bin = cnt_np_bin(cat_obs, param)
    println("------------------------------")
    ess_bin = stellar_ess(param)
    println("------------------------------")

    println("Simple Bayesian Rates:")
    for i in 1:(length(limitP)-1)
        for j in 1:(length(limitRp)-1)
            rate_f = (1.0+cnt_bin[(i-1)*(length(limitRp)-1) + j])/(1.0+ess_bin[(i-1)*(length(limitRp)-1) + j])*100.
            up_quant = quantile(Gamma(1.0+cnt_bin[(i-1)*(length(limitRp)-1) + j], 1.0/(1.0+ess_bin[(i-1)*(length(limitRp)-1) + j])), 0.8413)*100.
            low_quant = quantile(Gamma(1.0+cnt_bin[(i-1)*(length(limitRp)-1) + j], 1.0/(1.0+ess_bin[(i-1)*(length(limitRp)-1) + j])), 0.1587)*100.
            println(rate_f,
                    " + ", up_quant - rate_f,
                    " - ", rate_f - low_quant, " %")
        end
    end
    println()
end

function inv_det_simp_bayes(cat_obs::KeplerObsCatalog, param::SimParam)
    num_targ = ExoplanetsSysSim.StellarTable.num_usable_in_star_table()

    limitP::Array{Float64,1} = get_any(param, "p_lim_arr", Array{Float64,1})
    limitRp::Array{Float64,1} = get_any(param, "r_lim_arr", Array{Float64,1})

    println("------------------------------")
    cnt_bin, np_bin = cnt_np_bin(cat_obs, param)
    println("------------------------------")
    ess_bin = stellar_ess(param)
    println("------------------------------")

    println("Inverse Detection Rates:")
    for i in 1:(length(limitP)-1)
        for j in 1:(length(limitRp)-1)
            rate_f = np_bin[(i-1)*(length(limitRp)-1) + j]/num_targ*100.
            if cnt_bin[(i-1)*(length(limitRp)-1) + j] > 0.
                println(rate_f,
                        " +/- ", rate_f/sqrt(cnt_bin[(i-1)*(length(limitRp)-1) + j]), " %")
            else
                println(rate_f,
                        " +/- N/A %")
            end
        end
    end

    println()
    println("Simple Bayesian Rates:")
    for i in 1:(length(limitP)-1)
        for j in 1:(length(limitRp)-1)
            rate_f = (1.0+cnt_bin[(i-1)*(length(limitRp)-1) + j])/(1.0+ess_bin[(i-1)*(length(limitRp)-1) + j])*100.
            up_quant = quantile(Gamma(1.0+cnt_bin[(i-1)*(length(limitRp)-1) + j], 1.0/(1.0+ess_bin[(i-1)*(length(limitRp)-1) + j])), 0.8413)*100.
            low_quant = quantile(Gamma(1.0+cnt_bin[(i-1)*(length(limitRp)-1) + j], 1.0/(1.0+ess_bin[(i-1)*(length(limitRp)-1) + j])), 0.1587)*100.
            println(rate_f,
                    " + ", up_quant - rate_f,
                    " - ", rate_f - low_quant, " %")
        end
    end
    println()
end

## cnt_bin & np_bin (inverse detection & simple bayesian)
function cnt_np_bin(cat_obs::KeplerObsCatalog, param::SimParam, verbose::Bool = true)
    num_targ = ExoplanetsSysSim.StellarTable.num_usable_in_star_table()
    idx_tranets = find(x::KeplerTargetObs-> length(x.obs) > 0, cat_obs.target)::Array{Int64,1}

    limitP::Array{Float64,1} = get_any(param, "p_lim_arr", Array{Float64,1})
    limitRp::Array{Float64,1} = get_any(param, "r_lim_arr", Array{Float64,1})

    np_bin = zeros((length(limitP)-1) * (length(limitRp)-1))
    cnt_bin = zeros((length(limitP)-1) * (length(limitRp)-1))
    pl_idx = 1

    println("Calculating completeness for each planet...")
    for i in idx_tranets
        for j in 1:num_planets(cat_obs.target[i])
            pper = cat_obs.target[i].obs[j].period
            prad = sqrt(cat_obs.target[i].obs[j].depth)*cat_obs.target[i].star.radius

            pbin = findfirst(x -> ((pper > limitP[x]) && (pper < limitP[x+1])), collect(1:(length(limitP)-1)))
            rbin = findfirst(x -> ((prad > limitRp[x]) && (prad < limitRp[x+1])), collect(1:(length(limitRp)-1)))
            if (pbin > 0 && rbin > 0)
                cnt_bin[(pbin-1)*(length(limitRp)-1) + rbin] += 1
                pgeo = ExoplanetsSysSim.calc_transit_prob_single_planet_approx(pper, cat_obs.target[i].star.radius, cat_obs.target[i].star.mass)
	        pdet = 0.0
	        for star_id in 1:num_targ
	            star = SingleStar(ExoplanetsSysSim.StellarTable.star_table(star_id,:radius),ExoplanetsSysSim.StellarTable.star_table(star_id,:mass),1.0, star_id)
	            cdpp = 1.0e-6 * ExoplanetsSysSim.StellarTable.star_table(star_id, :rrmscdpp04p5) * sqrt(4.5/24.0 / ExoplanetsSysSim.LC_duration )
	            contam = 0.0
	            data_span = ExoplanetsSysSim.StellarTable.star_table(star_id, :dataspan)
	            duty_cycle = ExoplanetsSysSim.StellarTable.star_table(star_id, :dutycycle)
	            pl_arr = Array{Planet}( 1)
	            orbit_arr = Array{Orbit}( 1)
                    incl = acos(Base.rand()*star.radius*ExoplanetsSysSim.rsol_in_au/ExoplanetsSysSim.semimajor_axis(pper, star.mass))
	            orbit_arr[1] = Orbit(pper, 0., incl, 0., 0., Base.rand()*2.*pi)
	            pl_arr[1] = Planet(prad, 1.0e-6)
	            kep_targ = KeplerTarget([PlanetarySystem(star, pl_arr, orbit_arr)], fill(cdpp,ExoplanetsSysSim.num_cdpp_timescales,ExoplanetsSysSim.num_quarters),contam,data_span,duty_cycle)

	            duration_central = ExoplanetsSysSim.calc_transit_duration(kep_targ,1,1)
	            if duration_central <= 0.
	                continue
	            end
	            ntr = ExoplanetsSysSim.calc_expected_num_transits(kep_targ, 1, 1, param)
	            depth = ExoplanetsSysSim.calc_transit_depth(kep_targ,1,1)
                    snr_central = ExoplanetsSysSim.calc_snr_if_transit(kep_targ, depth, duration_central, param, num_transit=ntr)
	            pdet += ExoplanetsSysSim.calc_prob_detect_if_transit(kep_targ, snr_central, param, num_transit=ntr)
	        end
                np_bin[(pbin-1)*(length(limitRp)-1) + rbin] += 1.0/pgeo/(pdet/num_targ)
                if verbose
	            println("Planet ",pl_idx," => Bin ", (pbin-1)*(length(limitRp)-1) + rbin, ", C = ", 1.0/pgeo/(pdet/num_targ))
                end
	        pl_idx += 1
            end
        end
    end
    return cnt_bin, np_bin
end

## stellar catalog ess (simple bayesian)
function stellar_ess(param::SimParam, verbose::Bool = true)
  num_realiz = 100
  num_targ = ExoplanetsSysSim.StellarTable.num_usable_in_star_table()

  limitP::Array{Float64,1} = get_any(param, "p_lim_arr", Array{Float64,1})
  limitRp::Array{Float64,1} = get_any(param, "r_lim_arr", Array{Float64,1})

  ess_bin = zeros((length(limitP)-1) * (length(limitRp)-1))

  println(string("Stellar ESS calculation beginning..."))
  for star_id in 1:num_targ
    star = SingleStar(ExoplanetsSysSim.StellarTable.star_table(star_id,:radius),ExoplanetsSysSim.StellarTable.star_table(star_id,:mass),1.0, star_id)
    cdpp = 1.0e-6 * ExoplanetsSysSim.StellarTable.star_table(star_id, :rrmscdpp04p5) * sqrt(4.5/24.0 / ExoplanetsSysSim.LC_duration )
    contam = 0.0
    data_span = ExoplanetsSysSim.StellarTable.star_table(star_id, :dataspan)
    duty_cycle = ExoplanetsSysSim.StellarTable.star_table(star_id, :dutycycle)

    for i_idx in 1:(length(limitP)-1)
      for j_idx in 1:(length(limitRp)-1)
        temp_bin = 0.0
        for n_test in 1:num_realiz
          pper = exp(Base.rand()*(log(limitP[i_idx+1])-log(limitP[i_idx]))+log(limitP[i_idx]))
	  prad = exp(Base.rand()*(log(limitRp[j_idx+1])-log(limitRp[j_idx]))+log(limitRp[j_idx]))

          pgeo = ExoplanetsSysSim.calc_transit_prob_single_planet_approx(pper, star.radius, star.mass)
	  pdet = 0.0

	  pl_arr = Array{Planet}(1)
	  orbit_arr = Array{Orbit}(1)
          incl = acos(Base.rand()*star.radius*ExoplanetsSysSim.rsol_in_au/ExoplanetsSysSim.semimajor_axis(pper, star.mass))
	  orbit_arr[1] = Orbit(pper, 0., incl, 0., 0., Base.rand()*2.*pi)
	  pl_arr[1] = Planet(prad, 1.0e-6)
	  kep_targ = KeplerTarget([PlanetarySystem(star, pl_arr, orbit_arr)], fill(cdpp,ExoplanetsSysSim.num_cdpp_timescales,ExoplanetsSysSim.num_quarters),contam,data_span,duty_cycle)

	  duration_central = ExoplanetsSysSim.calc_transit_duration(kep_targ,1,1)
	  if duration_central <= 0.
	    continue
	  end
	  ntr = ExoplanetsSysSim.calc_expected_num_transits(kep_targ, 1, 1, param)
	  depth = ExoplanetsSysSim.calc_transit_depth(kep_targ,1,1)
          snr_central = ExoplanetsSysSim.calc_snr_if_transit(kep_targ, depth, duration_central, param, num_transit=ntr)
	  pdet = ExoplanetsSysSim.calc_prob_detect_if_transit(kep_targ, snr_central, param, num_transit=ntr)

	  temp_bin += (pgeo*pdet)
        end
        ess_bin[(i_idx-1)*(length(limitRp)-1) + j_idx] += temp_bin/num_realiz
      end
    end
    if verbose && rem(star_id, 10000) == 0.
      println(string("Star #", star_id, " finished"))
    end
  end

  if verbose
      println("")
      for i in 1:(length(limitP)-1)
          for j in 1:(length(limitRp)-1)
              println("Period limits: ", limitP[i:i+1], " / Radius limits: ", limitRp[j:j+1]/ExoplanetsSysSim.earth_radius, " / Stellar ESS = ", ess_bin[(i-1)*(length(limitRp)-1) + j])
          end
      end
  end
  return ess_bin
end


text:
fsabs(x) = sum(abs, x)
fsum_abs_and_sq(x) = (sum(abs, x), sum(abs2, x))

@testset "minimization_problem()" begin

@testset "1-dimensional, single-objective sum(abs, x)" begin
    p = minimization_problem(fsabs, "sumabs", (-1.0, 1.0), 1, 0.0)

    @test fitness_scheme(p) == MinimizingFitnessScheme
    @test numobjectives(p) == 1
    @test numdims(p) == 1
    @test search_space(p) == RectSearchSpace(1, (-1.0, 1.0))

    @test fitness([0.0], p) == 0.0
    @test fitness([1.2], p) == 1.2
    @test fitness([-1.9], p) == 1.9
end

@testset "3-dimensional, single-objective sum(abs, x)" begin
    p = minimization_problem(fsabs, "sumabs", (-1.0, 1.0), 3, 0.0)

    @test fitness_scheme(p) == MinimizingFitnessScheme
    @test numobjectives(p) == 1
    @test numdims(p) == 3
    @test search_space(p) == RectSearchSpace(3, (-1.0, 1.0))

    @test fitness([0.0, 1.0, 2.0], p) == 3.0
    @test fitness([-1.0, 1.0, 2.0], p) == 4.0
end

@testset "1-dimensional, multi-objective sumabs_sumsq" begin
    ss = RectSearchSpace(1)
    p = FunctionBasedProblem(fsum_abs_and_sq, "sumabs_sumsq",
                             ParetoFitnessScheme{2}(), ss, (0.0, 0.0))

    @test numdims(p) == 1
    @test numobjectives(p) == 2
    @test search_space(p) == ss

    @test fitness([0.0], p) == (0.0, 0.0)
    @test fitness([1.2], p) == (1.2, 1.2^2)
    @test fitness([-1.9], p) == (1.9, 1.9^2)
end

@testset "MinimizationProblemFamily" begin
    pfam = MinimizationProblemFamily(fsabs, "sumabs", (0.0, 1.0), 0.0)
    p1 = instantiate(pfam, 1)

    @test numobjectives(p1) == 1
    @test numdims(p1) == 1

    @test fitness([0.0], p1) == 0.0
    @test fitness([1.2], p1) == 1.2
    @test fitness([-1.9], p1) == 1.9

    p3 = instantiate(pfam, 3)

    @test numobjectives(p3) == 1
    @test numdims(p3) == 3

    @test fitness([0.0, 1.0, 2.0], p3) == 3.0
    @test fitness([-1.0, 1.0, 2.0], p3) == 4.0
end

end

@testset "ShiftedAndBiasedProblem" begin

@testset "Only x shifted 1-dim" begin
    p = minimization_problem(fsabs, "sumabs", (-1.0, 1.0), 1, 0.0)
    sp = BlackBoxOptim.ShiftedAndBiasedProblem(p, xshift=[0.5])

    @test BlackBoxOptim.orig_problem(sp) === p
    @test numobjectives(sp) == 1
    @test numdims(sp) == 1
    @test search_space(sp) == RectSearchSpace(1, (-1.0, 1.0))
    xs = sp.xshift
    @test xs == [0.5]
    @test sp.fitshift == 0.0

    @test fitness([0.0], sp) == 0.5
    @test fitness([1.2], sp) == 0.7
    @test fitness([-1.9], sp) == 2.4
end

@testset "Shifted and biased 2-dim" begin
    ss = RectSearchSpace(2, (-0.5, 1.0))
    p = minimization_problem(fsabs, "sumabs", (-0.5, 1.0), 2, 0.0)
    sp = BlackBoxOptim.ShiftedAndBiasedProblem(p; fitshift = 1.3)

    @test BlackBoxOptim.orig_problem(sp) === p
    @test numobjectives(sp) == 1
    @test numdims(sp) == 2
    @test sp.xshift == [0.0, 0.0]
    @test sp.fitshift == 1.3
    @test search_space(sp) == RectSearchSpace(2, (-0.5, 1.0))

    xs = sp.xshift
    @test fitness([0.0, 1.0], sp) == 1.0 + 1.3
    @test fitness([1.2, -1.3], sp) == 2.5 + 1.3
    @test fitness([-1.9, 0.0], sp) == 1.9 + 1.3
end

@testset "Within ftol" begin
    subp = minimization_problem(fsabs, "sumabs", (-1.0, 1.0), 2, 0.0)

    @test fitness_is_within_ftol(subp, 0.1, 0.2)
    @test !fitness_is_within_ftol(subp, 0.1, 0.09)
end

end

ext: jl


text:
using Test
using LinearAlgebra
using Random
using KrylovKit

const n = 10
const N = 100

const η₀   = 1/sqrt(2) # conservative choice, probably 1/2 is sufficient
const cgs = ClassicalGramSchmidt()
const mgs = ModifiedGramSchmidt()
const cgs2 = ClassicalGramSchmidt2()
const mgs2 = ModifiedGramSchmidt2()
const cgsr = ClassicalGramSchmidtIR(η₀)
const mgsr = ModifiedGramSchmidtIR(η₀)

Random.seed!(1234567);

include("linalg.jl")
include("factorize.jl")
include("eigsolve.jl")
include("schursolve.jl")
include("svdsolve.jl")
include("linsolve.jl")
include("exponentiate.jl")
include("recursivevec.jl")


text:
########################################################################
# Atoms of the trace-norm ball.
########################################################################
"""
    TraceBallAtom(n, v)

Atom in the TraceBall of dimension n×n.
"""
struct TraceBallAtom{T1<:Int64, T2<:Vector{Float64}} <: AbstractAtom
    n::T1
    v::T2
    function TraceBallAtom(v::Vector{Float64})
        n = length(v)
        new{Int64, Vector{Float64}}(n, v)
    end
end

"""
    *(M::AbstractLinearOp, a::TraceBallAtom) -> Y

Multiply a trace-norm atom by a linear map, and return the matrix `Y`.
"""
Base.:(*)(M::LinearOp, a::TraceBallAtom{Int64, Vector{Float64}}) = M*vec(a.v*a.v')
LinearAlgebra.mul!(Ma::Vector{Float64}, M::LinearOp, a::TraceBallAtom{Int64, Vector{Float64}}) = mul!(Ma, M, vec(a.v*a.v'))
Base.vec(a::TraceBallAtom{Int64, Vector{Float64}}) = vec(a.v*a.v')
Base.length(a::TraceBallAtom{Int64, Vector{Float64}}) = a.n * a.n
Base.size(a::TraceBallAtom{Int64, Vector{Float64}}) = a.n, a.n

"""
    TraceBall(n)

Atomic set defined by trace-norm ball in `n by n` matrices.
"""
struct TraceBall{T<:Int64} <: AbstractAtomicSet
    n::T
    maxrank::T
    function TraceBall(n::Int64, maxrank::Int64)
        n ≥ 1 || throw(DomainError(n,"n must be a positive integer"))
        n ≥ maxrank ≥ 0 || throw(DomainError(maxrank,"maxrank must be ≥ 0"))
        new{Int64}(n, maxrank)
    end
end
TraceBall(n::Int64; maxrank=n) = TraceBall(n, maxrank)

"""
    gauge(A::TraceBall, x::Matrix)

Gives the trace-norm of `x`.
"""
gauge(::TraceBall{Int64}, x::Matrix{Float64}) = tr(x)
gauge(A::TraceBall{Int64}, x::Vector{Float64}) = gauge(A, reshape(x, A.n, A.n))

"""
    support(A::TraceBall, z::AbstractMatrix)

Gives the largest singular value of `z`.
"""
support(::TraceBall{Int64}, z::AbstractMatrix{Float64}) = real(eigs(z, nev = 1, which=:LR)[1][1])
support(A::TraceBall{Int64}, z::AbstractVector{Float64}) = support(A, reshape(z, A.n, A.n))

"""
    expose!(A::TraceBall, z::Matrix, a::TraceBallAtom)

Gives the top eigen vector of `z`.
"""
function expose!(::TraceBall{Int64}, z::AbstractMatrix{Float64}, a::TraceBallAtom{Int64, Vector{Float64}})
    ~, v = eigs(z, nev = 1, which=:LR)
    a.v .= v[:]
    return nothing
end
expose!(A::TraceBall{Int64}, z::AbstractVector{Float64}, a::TraceBallAtom{Int64, Vector{Float64}}) = expose!(A, reshape(z, A.n, A.n), a)

function expose(A::TraceBall{Int64}, z::AbstractMatrix{Float64})
    n = A.n
    a = TraceBallAtom(zeros(n))
    expose!(A, z, a)
    return a
end
expose(A::TraceBall{Int64}, z::AbstractVector{Float64}) = expose(A, reshape(z, A.n, A.n))

Base.length(A::TraceBall{Int64}) = A.n*A.n
Base.size(A::TraceBall{Int64}) = A.n, A.n
rank(A::TraceBall{Int64}) = A.maxrank
atom_name(A::TraceBall{Int64}) = "trace-norm ball"
atom_description(A::TraceBall{Int64}) = "{ vv^T | v ∈ ℝⁿ, ||v||₂ ≤ 1 }"
atom_parameters(A::TraceBall{Int64}) = "n = $(A.n), maxrank = $(A.maxrank)"

########################################################################
# Face of the trace-norm ball.
########################################################################

mutable struct TraceBallFace{T1<:Int64, T2<:Matrix{Float64}} <: AbstractFace
    n::T1
    r::T1
    V::T2
    function TraceBallFace(V::Matrix{Float64})
        n, r = size(V)
        new{Int64, Matrix{Float64}}(n, r, V)
    end
end

"""
    face(A, z) -> F

Return a face `F` of the atomic set `A` exposed by the vector `z`. The dimension of
the exposed face is limited by `k=maxrank(A)`. If the dimension is being limited
by `k`, then the routine returns a set of rank `k` atoms that define the subset
of the exposed face.
"""
function face!(A::TraceBall{Int64}, z::AbstractMatrix{Float64}, F::TraceBallFace{Int64, Matrix{Float64}})
    r = A.maxrank
    if r == 0
        V = zeros(A.n, 1)
    elseif r < min(A.n) - 1
        ~, V = eigs(z, nev=r, which=:LR)
        F.V .= V
    else
        V = eigvecs(z)
        F.V .= V[:, 1:r]
    end
    return nothing
end
face!(A::TraceBall{Int64}, z::AbstractVector{Float64}, F::TraceBallFace{Int64, Matrix{Float64}}) = face!(A, reshape(z, A.n, A.n), F)

function face(A::TraceBall{Int64}, z::AbstractMatrix{Float64})
    n, r = A.n, max(A.maxrank, 1)
    V = zeros(n, r)
    F = TraceBallFace(V)
    face!(A, z, F)
    return F
end
face(A::TraceBall{Int64}, z::AbstractVector{Float64}) = face(A, reshape(z, A.n, A.n))

"""
Given a face and a set of weights, reveal the corresponding
point on the face.
"""
Base.:(*)(F::TraceBallFace{Int64, Matrix{Float64}}, S::Matrix{Float64}) = vec(F.U*S*F.U')
Base.:(*)(F::TraceBallFace{Int64, Matrix{Float64}}, c::Vector{Float64}) = F*reshape(c, F.r, F.r)

"""
    *(M::LinearOp, F::TraceBallFace)

Return a LinearMap `L` whose forward and adjoing operators
are defined as

    L*p = M*vec(V*matrix(p)*V')
    L'*q = vec(V'*matrix(M'*q)*V)

where the `n`-by-`k` matrice `V` is the left
eigen vectors that define the face `F` of the trace-norm ball.
"""
function Base.:(*)(M::LinearOp, F::TraceBallFace{Int64, Matrix{Float64}})
    n, r = F.n, F.r
    f = p->M*vec(F.V*reshape(p, r, r)*F.V')
    fc = q->vec(F.V'*reshape(M'*q, n, n)*F.V)
    return LinearMap(f, fc, size(M, 1), r*r)
end

Base.:(*)(λ::Real, F::TraceBallFace{Int64, Matrix{Float64}}) = λ*vec(F)
Base.length(F::TraceBallFace{Int64, Matrix{Float64}}) = F.n*F.n
Base.size(F::TraceBallFace{Int64, Matrix{Float64}}) = F.n, F.n
rank(F::TraceBallFace{Int64, Matrix{Float64}}) = F.r^2
vec(F::TraceBallFace{Int64, Matrix{Float64}}) = I(F.n*F.n)*F
face_name(F::TraceBallFace{Int64, Matrix{Float64}}) = "face of trace-norm ball"
face_parameters(F::TraceBallFace{Int64, Matrix{Float64}}) = "n = $(F.n), rank = $(rank(F))"


text:
export Ball

mutable struct Ball{N,N²} <: Shape{N,N²}
    c::SVector{N,Float64}  # center of ball
    r::Float64  # radius
    Ball{N,N²}(c,r) where {N,N²} = new(c,r)  # suppress default outer constructor
end

Ball(c::SVector{N,<:Real}, r::Real) where {N} = Ball{N,N*N}(c, r)
Ball(c::AbstractVector{<:Real}, r::Real) = (N = length(c); Ball(SVector{N}(c), r))

Base.:(==)(s1::Ball, s2::Ball) = s1.c==s2.c && s1.r==s2.r
Base.isapprox(s1::Ball, s2::Ball) = s1.c≈s2.c && s1.r≈s2.r
Base.hash(s::Ball, h::UInt) = hash(s.c, hash(s.r, hash(:Ball, h)))

level(x::SVector{N,<:Real}, s::Ball{N}) where {N} = √(sum(abs2, x - s.c) / s.r^2) - 1.0

function surfpt_nearby(x::SVector{N,<:Real}, s::Ball{N}) where {N}
    nout = x==s.c ? SVector(ntuple(k -> k==1 ? 1.0 : 0.0, Val(N))) :  # nout = e₁ for x == s.c
                    normalize(x-s.c)
    return s.c+s.r*nout, nout
end

bounds(s::Ball) = (s.c.-s.r, s.c.+s.r)


text:
function param_values(f::AbstractParameterizedFunction)
  [getfield(f,s) for s in f.params]
end

num_params(f::AbstractParameterizedFunction) = length(f.params)

# Fallbacks

param_values(f) = nothing
num_params(f) = 0

function problem_new_parameters(prob::ODEProblem,p)
  f = (t,u,du) -> prob.f(t,u,p,du)
  uEltype = eltype(p)
  u0 = [uEltype(prob.u0[i]) for i in 1:length(prob.u0)]
  tspan = (uEltype(prob.tspan[1]),uEltype(prob.tspan[2]))
  ODEProblem(f,u0,tspan)
end
param_values(prob::ODEProblem) = param_values(prob.f)
num_params(prob::ODEProblem) = num_params(prob.f)

function problem_new_parameters(prob::DAEProblem,p)
  f = (t,u,du,resid) -> prob.f(t,u,p,du,resid)
  uEltype = eltype(p)
  u0 = [uEltype(prob.u0[i]) for i in 1:length(prob.u0)]
  tspan = (uEltype(prob.tspan[1]),uEltype(prob.tspan[2]))
  DAEProblem(f,u0,tspan)
end
param_values(prob::DAEProblem) = param_values(prob.f)
num_params(prob::DAEProblem) = num_params(prob.f)

function problem_new_parameters(prob::ConstantLagDDEProblem,p)
  f = (t,u,h,du) -> prob.f(t,u,h,p,du)
  uEltype = eltype(p)
  u0 = [uEltype(prob.u0[i]) for i in 1:length(prob.u0)]
  tspan = (uEltype(prob.tspan[1]),uEltype(prob.tspan[2]))
  DAEProblem(f,prob.h,u0,prob.lags,tspan)
end
param_values(prob::ConstantLagDDEProblem) = param_values(prob.f)
num_params(prob::ConstantLagDDEProblem) = num_params(prob.f)

function problem_new_parameters(prob::SDEProblem,p)
  fpars = num_params(prob.f)
  if fpars > 0
    f = (t,u,du) -> prob.f(t,u,@view(p[1:fpars]),du)
    if num_params(prob.g) > 0
      g = (t,u,du) -> prob.g(t,u,@view(p[(fpars+1):end]),du)
    else
      g = prob.g
    end
  else
    f = prob.f
    g = (t,u,du) -> prob.g(t,u,p,du)
  end
  uEltype = eltype(p)
  u0 = [uEltype(prob.u0[i]) for i in 1:length(prob.u0)]
  tspan = (uEltype(prob.tspan[1]),uEltype(prob.tspan[2]))
  SDEProblem(f,g,u0,tspan)
end
param_values(prob::SDEProblem) = (A = [param_values(prob.f) ; param_values(prob.g)]; [A.!=nothing])
num_params(prob::SDEProblem) = num_params(prob.f) + num_params(prob.g)

function problem_new_parameters(prob::MonteCarloProblem,p)
  probpars = num_params(prob.prob)
  pfpars = num_params(prob.prob_func)
  opars = num_params(prob.output_func)
  if probpars > 0
    tmp_prob = problem_new_parameters(prob.prob,@view(p[1:probpars]))
    if pfpars > 0
      prob_func = (prob,i) -> prob.prob_func(prob,i,@view(p[(probpars+1):(probpars+pfpars)]))
      opars>0 ? output_func = prob.output_func(sol,@view(p[(probpars+pfpars+1):end])) : output_func = prob.output_func
    else
      prob_func = prob.prob_func
      opars>0 ? output_func = prob.output_func(sol,@view(p[(probpars+pfpars+1):end])) : output_func = prob.output_func
    end
  else
    tmp_prob = prob.prob
    if pfpars > 0
      prob_func = (prob,i) -> prob.prob_func(prob,i,@view(p[1:pfpars]))
      opars>0 ? output_func = prob.output_func(sol,@view(p[(pfpars+1):end])) : output_func = prob.output_func
    else
      prob_func = prob.prob_func
      output_func = prob.output_func(sol,p)
    end
  end
  MonteCarloProblem(tmp_prob,prob_func,output_func)
end

param_values(prob::MonteCarloProblem) = (A = [param_values(prob.prob) ; param_values(prob.prob_func) ; param_values(prob.output_func)]; [A.!=nothing])
num_params(prob::MonteCarloProblem) = num_params(prob.prob) + num_params(prob.prob_func) + num_params(prob.output_func)


text:
function rand(rng::AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{DoubleFloat{T}}}) where {T<:IEEEFloat}
    hi, lo  = rand(rng, T, 2)
    if hi === zero(T)
        if lo === zero(T)
            return zero(DoubleFloat(T))
        end
        hi, lo = lo, hi
    end
    frlo, xplo  = frexp(lo)
    xplo = Base.exponent(hi) - min(1, fld(xplo,4)) - abs(Base.exponent(eps(hi)))
    lo = ldexp(frlo, xplo)
    lo = rand(rng, Bool) ? lo : -lo

    DoubleFloat(hi, lo)
end

function rand(rng::AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{Complex{DoubleFloat{T}}}}) where {T<:IEEEFloat}
    re = rand(rng, DoubleFloat{T})
    im = rand(rng, DoubleFloat{T})
    return Complex{DoubleFloat{T}}(re, im)
end

function randpm(::Type{DoubleFloat{T}}) where {T<:IEEEFloat}
    r = rand(DoubleFloat{T})
    r = rand(Bool) ? r : -r
    return r
end

function randpm(rng::MersenneTwister, ::Type{DoubleFloat{T}}) where {T<:IEEEFloat}
    r = rand(rng, DoubleFloat{T})
    r = rand(rng, Bool) ? r : -r
    return r
end

function randpm(::Type{Complex{DoubleFloat{T}}}) where {T<:IEEEFloat}
    re = randpm(DoubleFloat{T})
    im = randpm(DoubleFloat{T})
    return Complex{DoubleFloat{T}}(re, im)
end

function randpm(rng::MersenneTwister, ::Type{Complex{DoubleFloat{T}}}) where {T<:IEEEFloat}
    re = randpm(rng, DoubleFloat{T})
    im = randpm(rng, DoubleFloat{T})
    return Complex{DoubleFloat{T}}(re, im)
end

function randpm(::Type{DoubleFloat{T}}, n::Int) where {T<:IEEEFloat}
    rs = rand(DoubleFloat{T}, n)
    sgns = ones(n) .- (2 .* rand(Bool, n))
    return rs .* sgns
end

function randpm(::Type{Complex{DoubleFloat{T}}}, n::Int) where {T<:IEEEFloat}
    res = rand(DoubleFloat{T}, n)
    ims = rand(DoubleFloat{T}, n)
    o = ones(n)
    rsgns = o .- (2 .* rand(Bool,n))
    isgns = o .- (2 .* rand(Bool,n))
    res = res .* rsgns
    ims = ims .* isgns
    return map((re,im)->Complex{DoubleFloat{T}}(re,im), res, ims)
end

function randpm(rng::MersenneTwister, ::Type{DoubleFloat{T}}, n::Int) where {T<:IEEEFloat}
    rs = rand(rng, DoubleFloat{T}, n)
    sgns = ones(n) .- (2 .* rand(rng, Bool, n))
    return rs .* sgns
end

function randpm(rng::MersenneTwister, ::Type{Complex{DoubleFloat{T}}}, n::Int) where {T<:IEEEFloat}
    res = rand(rng, DoubleFloat{T}, n)
    ims = rand(rng, DoubleFloat{T}, n)
    o = ones(n)
    rsgns = o .- (2 .* rand(rng, Bool,n))
    isgns = o .- (2 .* rand(rng, Bool,n))
    res = res .* rsgns
    ims = ims .* isgns
    return map((re,im)->Complex{DoubleFloat{T}}(re,im), res, ims)
end

# normal variates

function randn(rng::AbstractRNG, ::Type{DoubleFloat{T}}) where {T<:IEEEFloat}
    urand1, urand2 = rand(rng, DoubleFloat{T}, 2)
    urand1 = urand1 + urand1 - 1
    urand2 = urand2 + urand2 - 1
    s = urand1*urand1 + urand2*urand2

    while s >= 1 || s === 0
        urand1, urand2 = rand(rng, DoubleFloat{T}, 2)
        urand1 = urand1 + urand1 - 1
        urand2 = urand2 + urand2 - 1
        s = urand1*urand1 + urand2*urand2
    end

    s = sqrt( -log(s) / s )
    return (urand1 + urand2) * s
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "PulseAudio"
version = v"14.2.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://freedesktop.org/software/pulseaudio/releases/pulseaudio-14.2.tar.gz", "902dd1928801bb5dc7b121754aa4110ce55768b5dff94a700e7bd58d3f597970")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
apk add gettext
apk add glib
apk add orc-compiler
apk add perl-xml-parser
apk add bash-completion
# make sure meson can find everything
sed -i -e "s~c_args = .*~c_args = ['-I${includedir}', '-L${libdir}']~" ${MESON_TARGET_TOOLCHAIN}
# I guess pulseaudio doesn't set install_rpath correctly?
find pulseaudio-* -type f | xargs sed -i "s~install_rpath : privlibdir~install_rpath : '\$ORIGIN/pulseaudio'~"
# For some reason, librt fails to get linked correctly, so add a flag
sed -i -e "s~c_link_args = .*~c_link_args = ['-lrt']~" ${MESON_TARGET_TOOLCHAIN}
cd pulseaudio-*
# Disable ffast-math; I repented
sed -i -e "s/link_args : \['-ffast-math'],//" src/daemon/meson.build
# pulseaudio seems to check for iconv_open but use libiconv_open?
sed -i -e "s/cc.has_function('iconv_open')/cc.has_function('libiconv_open')/" meson.build
# Force meson to use some libraries
if [[ "${target}" == powerpc64le-* ]]; then     sed -i -e "s~'sys/capability.h',~~"  meson.build; fi
mkdir build
cd build
# I can't figure out how to build tdb, see https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/887
# BlueZ is in progress, see https://github.com/JuliaPackaging/Yggdrasil/pull/2582
meson ..  -Ddatabase="gdbm" -Dbluez5="false" --cross-file=${MESON_TARGET_TOOLCHAIN}
ninja
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(Sys.islinux, supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libpulse-simple", :libpulse_simple),
    ExecutableProduct("pasuspender", :pasuspender),
    ExecutableProduct("pacmd", :pacmd),
    ExecutableProduct("pactl", :pactl),
    ExecutableProduct("pulseaudio", :pulseaudio),
    LibraryProduct("libpulse", :libpulse),
    LibraryProduct("libpulse-mainloop-glib", :libpulse_mainloop_glib),
    ExecutableProduct("pacat", :pacat)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="alsa_jll", uuid="45378030-f8ea-5b20-a7c7-1a9d95efb90e"))
    Dependency(PackageSpec(name="Check_jll", uuid="491db154-c145-5abe-9c32-446728d60cce"))
    Dependency(PackageSpec(name="Dbus_jll", uuid="ee1fde0b-3d02-5ea6-8484-8dfef6360eab"))
    Dependency(PackageSpec(name="eudev_jll", uuid="35ca27e7-8b34-5b7f-bca9-bdc33f59eb06"))
    Dependency(PackageSpec(name="FFTW_jll", uuid="f5851436-0d7a-5f13-b9de-f02708fd171a"))
    Dependency(PackageSpec(name="Gdbm_jll", uuid="54ca2031-c8dd-5cab-9ed4-295edde1660f"))
    Dependency(PackageSpec(name="Gettext_jll", uuid="78b55507-aeef-58d4-861c-77aaff3498b1"))
    Dependency(PackageSpec(name="Glib_jll", uuid="7746bdde-850d-59dc-9ae8-88ece973131d"))
    Dependency(PackageSpec(name="GStreamer_jll", uuid="aaaaf01e-2457-52c6-9fe8-886f7267d736"))
    Dependency(PackageSpec(name="libsndfile_jll", uuid="5bf562c0-5a39-5b4f-b979-f64ac885830c"))
    Dependency(PackageSpec(name="libasyncns_jll", uuid="ed080073-db63-57db-a029-74e11ae80737"))
    Dependency(PackageSpec(name="libcap_jll", uuid="eef66a8b-8d7a-5724-a8d2-7c31ae1e29ed"))
    Dependency(PackageSpec(name="Libiconv_jll", uuid="94ce4f54-9a6c-5748-9c1c-f9c7231a4531"))
    Dependency(PackageSpec(name="Libtool_jll", uuid="a76c16ae-fb8f-5ff0-8826-da3b7a640f0b"))
    Dependency(PackageSpec(name="OpenSSL_jll", uuid="458c3c95-2e84-50aa-8efc-19380b2a3a95"))
    Dependency(PackageSpec(name="ORC_jll", uuid="fb41591b-4dee-5dae-bf56-d83afd04fbc0"))
    Dependency(PackageSpec(name="SBC_jll", uuid="da37f231-8920-5702-a09a-bdd970cb6ddc"))
    Dependency(PackageSpec(name="SoXResampler_jll", uuid="fbe68eb6-6641-54c6-99e3-f7c7c4d73a57"))
    Dependency(PackageSpec(name="SpeexDSP_jll", uuid="f2f9631b-9a4e-5b48-9975-88f638ec36a7"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
# Should not be exposed to the user.
const StringRef = MLIR.API.MlirStringRef
create_string_ref(str::String) = MLIR.API.mlirStringRefCreateFromCString(str)
StringRef(str::String) = create_string_ref(str)

# Unwraps API structs to ptr field.
unwrap(o) = getfield(o, :ptr)
unwrap(ptr::P) where P <: Ptr = ptr


text:
using Base.Threads

neighborhood_count_matrix(data::Union{BmmData, T} where T <: AbstractDataFrame, k::Int; kwargs...) =
    neighborhood_count_matrix(position_data(data), composition_data(data), k; kwargs...)

function neighborhood_count_matrix(pos_data::Matrix{<:Real}, genes::Vector{Int}, k::Int; confidences::Union{Vector{Float64}, Nothing}=nothing,
    n_genes::Int=maximum(genes), normalize_by_dist::Bool=true, normalize::Bool=true)
    if k < 3
        @warn "Too small value of k: $k. Setting it to 3."
        k = 3
    end

    if !normalize
        normalize_by_dist = false
    end

    if (confidences !== nothing) && normalize
        @warn "`confidences` are not supported with `normalize=true`. Ignoring them."
        confidences = nothing
    end

    k = min(k, size(pos_data, 2))

    neighbors, dists = knn(KDTree(pos_data), pos_data, k, true);

    n_cm = zeros(((normalize || (confidences !== nothing)) ? Float64 : Int), n_genes, size(pos_data, 2));

    if !normalize_by_dist
        @threads for (i,ids) in collect(enumerate(neighbors))
            if confidences !== nothing
                count_array!(view(n_cm, :, i), genes[ids], confidences[ids])
            elseif !normalize
                count_array!(view(n_cm, :, i), genes[ids])
            else
                prob_array!(view(n_cm, :, i), genes[ids])
            end
        end

        return n_cm
    end

    # normalize_by_dist doesn't affect result much, but neither it slows the work down. In theory, should work better for border cases.
    med_closest_dist = median([d[d .> 1e-15][1] for d in dists if any(d .> 1e-15)]) # account for problems with points with duplicating coordinates
    @threads for (i,(ids, dists)) in collect(enumerate(zip(neighbors, dists)))
        c_probs = view(n_cm, :, i)
        for (gene, dist) in zip(genes[ids], dists)
            c_probs[gene] += 1 / max(dist, med_closest_dist)
        end
    end

    return n_cm ./ sum(n_cm, dims=1);
end

truncate_pca(pca::MultivariateStats.PCA, outdim::Int) =
    MultivariateStats.PCA(pca.mean, pca.proj[:, 1:outdim], pca.prinvars[1:outdim], pca.tvar)

function gene_composition_transformation(count_matrix::Matrix{Float64}, confidence::Vector{Float64}=ones(size(count_matrix, 2));
        sample_size::Int=10000, seed::Int=42, method::Symbol=:umap, return_all::Bool=false, n_pcs::Int=15, out_dim::Int=3, kwargs...)
    sample_size = min(sample_size, size(count_matrix, 2))
    Random.seed!(seed)

    if n_pcs <= 3
        n_pcs = 3
        method = :pca
    end

    pca = fit(MultivariateStats.PCA, count_matrix, maxoutdim=n_pcs)
    pca_res = transform(truncate_pca(pca, 2), count_matrix);
    sample_ids = select_ids_uniformly(pca_res', confidence, n=sample_size)

    count_matrix_sample = count_matrix[:,sample_ids]

    emb = nothing
    if method == :umap
        emb = fit(UmapFit, count_matrix_sample, pca; n_components=out_dim, kwargs...);
    elseif method == :pca
        emb = truncate_pca(pca, 3);
    else
        error("Unknown method: '$method'")
    end

    if return_all
        return (emb, sample_ids, pca)
    end

    return emb
end

gene_composition_colors(count_matrix::Matrix{Float64}, transformation::UmapFit; kwargs...) =
    gene_composition_colors!(MultivariateStats.transform(transformation, count_matrix); kwargs...)

gene_composition_colors(embedding::Matrix{Float64}; kwargs...) =
    gene_composition_colors!(deepcopy(embedding); kwargs...)

function gene_composition_colors!(embedding::Matrix{Float64}; lrange::Tuple{<:Real, <:Real}=(10, 90), log_colors::Bool=false, trim_frac::Float64=0.0125)
    @assert size(embedding, 1) == 3 "Color embedding must have exactly 3 rows"

    embedding .-= mapslices(x -> quantile(x, trim_frac), embedding, dims=2)
    embedding .= max.(embedding, 0.0)
    max_val = quantile(vec(embedding), 1.0 - trim_frac)
    embedding ./= max_val;
    embedding .= min.(embedding, 1.0)

    if log_colors
        embedding .= log10.(embedding .+ max(quantile(vec(embedding), 0.05), 1e-3))
        embedding .-= minimum(embedding, dims=2)
        embedding ./= maximum(embedding, dims=2)
    end

    embedding[1,:] .*= lrange[2] - lrange[1]
    embedding[1,:] .+= lrange[1]

    embedding[2:3,:] .-= 0.5
    embedding[2:3,:] .*= 200

    return vec(mapslices(col -> Colors.Lab(col...), embedding, dims=1))
end

function gene_composition_colors(df_spatial::T where T <: AbstractDataFrame, k::Int; kwargs...)
    neighb_cm = neighborhood_count_matrix(df_spatial, k);
    confidence = (:confidence in propertynames(df_spatial)) ? df_spatial.confidence : ones(size(df_spatial, 1))
    transformation = gene_composition_transformation(neighb_cm, confidence)
    return gene_composition_colors(neighb_cm, transformation; kwargs...)
end


text:
export find_next_delim

"""
    Use: find_next_delim(x::AbstractVector, i::Int, probe_length::Int)

From starting index `i`, return the index of the next element in the vector that holds a different value.

# EXAMPLE
```
julia> x = vec([1 1 1 2 2 2 2 3 3 3 3]);

julia> find_next_delim(x, 1, 1)
4

julia> find_next_delim(x, 4, 1)
8
```
"""
function find_next_delim(x::AbstractVector, i::Int, probe_length::Int)

    # Setup
    j_prev = i
    j = j_prev + probe_length
    undershoot_counter = 0
    end_counter = 0
    n = length(x)
    probed_end_flag = false

    # Evaluate probe and update
    while true

        # Check inbounds, and update values
        if j <= n
            val = x[j_prev]
            probed_val = x[j]

            # Undershoot delim
            if val == probed_val
                j_prev = j
                j += probe_length
                undershoot_counter += 1

                # If keep undershooting, extend probe unless at end of x
                if (undershoot_counter >= 10) & ~probed_end_flag
                    probe_length *= 10
                    undershoot_counter = 0
                end

            # Delim found
            elseif j_prev + 1 == j
                return j

            # Overshoot delim
            else
                j = j_prev
                probe_length = max(1,floor(Int, probe_length/2))
            end
        # Probed out of bounds
        else
            j = n
            probe_length = 1
            probed_end_flag = true
            (end_counter>10) ? (return j) : (end_counter += 1)
        end
    end
end


text:
export SD_AIO_clockIOconfig, SD_AIO_clockGetFrequency, SD_AIO_clockGetSyncFrequency,
	SD_AIO_clockSetFrequency, SD_AIO_clockResetPhase, SD_AIO_triggerIOconfig,
	SD_AIO_triggerIOwrite, SD_AIO_triggerIOread
# Clock
## int __cdecl SD_AIO_clockIOconfig(int moduleID, int port, int clockConfig);
SD_AIO_clockIOconfig(moduleID::Integer, port::Integer, clockConfig::Integer) =
	ccall((:SD_AIO_clockIOconfig, klib), Cint, (Cint, Cint, Cint),
		moduleID, port, clockConfig)

## double SD_AIO_clockGetFrequency(int moduleID, int port);
SD_AIO_clockGetFrequency(moduleID::Integer, port::Integer) =
	ccall((:SD_AIO_clockGetFrequency, klib), Cdouble, (Cint, Cint),
		moduleID, port)

## double SD_AIO_clockGetSyncFrequency(int moduleID, int port);
SD_AIO_clockGetSyncFrequency(moduleID::Integer, port::Integer) =
	ccall((:SD_AIO_clockGetSyncFrequency, klib), Cdouble, (Cint, Cint),
		moduleID, port)

## double SD_AIO_clockSetFrequency(int moduleID, int port, double frequency, int mode SD_DEFAULT_CLOCK_FREQ_MODE);
SD_AIO_clockSetFrequency(moduleID::Integer, port::Integer, frequency::Real,
	mode::Integer) =
	ccall((:SD_AIO_clockSetFrequency, klib), Cdouble, (Cint, Cint, Cdouble, Cint),
		moduleID, port, frequency, mode)

## int SD_AIO_clockResetPhase(int moduleID, int port, int triggerBehavior, int PXItrigger, double skew SD_DEFAULT_CLOCK_SKEW);
SD_AIO_clockResetPhase(moduleID::Integer, port::Integer, triggerBehavior::Integer,
	PXItrigger::Integer, skew::Real) =
	ccall((:SD_AIO_clockResetPhase, klib), Cint,
		(Cint, Cint, Cint, Cint, Cdouble),
		moduleID, port, triggerBehavior, PXItrigger, skew)

# Trigger
## int SD_AIO_triggerIOconfig(int moduleID, int port, int direction);
SD_AIO_triggerIOconfig(moduleID::Integer, port::Integer, direction::Integer) =
	ccall((:SD_AIO_triggerIOconfig, klib), Cint, (Cint, Cint, Cint),
		moduleID, port, direction)

## int SD_AIO_triggerIOwrite(int moduleID, int port, int value, int syncMode);
SD_AIO_triggerIOwrite(moduleID::Integer, port::Integer, value::Integer,
	syncMode::Integer) =
	ccall((:SD_AIO_triggerIOwrite, klib), Cint, (Cint, Cint, Cint, Cint),
		moduleID, port, value, syncMode)

## int SD_AIO_triggerIOread(int moduleID, int port);
SD_AIO_triggerIOread(moduleID::Integer, port::Integer) =
	ccall((:SD_AIO_triggerIOread, klib), Cint, (Cint, Cint), moduleID, port)


text:
function slurp_grid(filename)
    permutedims(hcat(map(collect, readlines(filename))...))
end

eg = slurp_grid("3eg.txt")
input = slurp_grid("3a.txt")

flip(bitchar) = bitchar == '1' ? '0' : '1'

gamma(col) = length(filter(==('1'), col)) / length(col) >= 0.5 ? '1' : '0'

function gamma_rate(m)
    map(gamma, eachcol(m))
end

function to_decimal(bits)
    bitstring = string(bits...)
    parse(Int, bitstring, base = 2)
end

function solve1(m)
    g = gamma_rate(m)
    epsilon_rate = map(flip, g)
    to_decimal(g) * to_decimal(epsilon_rate)
end

solve1(eg)

solve1(input)

function select_max(m, idx)
    max = gamma(m[:, idx])
    permutedims(hcat([row for row in eachrow(m) if row[idx] == max]...))
end

function select_min(m, idx)
    min = flip(gamma(m[:, idx]))
    permutedims(hcat([row for row in eachrow(m) if row[idx] == min]...))
end

function oxygen(m)
    s = size(m, 2)
    for idx = 1:s
        m = select_max(m, idx)
        if size(m, 1) == 1
            return m
        end
    end
end

function co2(m)
    s = size(m, 2)
    for idx = 1:s
        m = select_min(m, idx)
        if size(m, 1) == 1
            return m
        end
    end
end

function solve2(m)
    to_decimal(oxygen(m)) * to_decimal(co2(m))
end

solve2(eg)
solve2(input) == 1613181


text:
using Documenter, ExtractMacro

makedocs()

deploydocs(
    deps   = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
    repo   = "github.com/carlobaldassi/ExtractMacro.jl.git",
    julia  = "release"
)


text:
using MultiModalMuSig
using Test

K = [2, 3]
α = [0.1, 0.1]
X = Vector{Matrix{Int}}[
    # document 1
    [
        # modality 1
        [
            1 5;
            2 8
        ],
        # modality 2
        [
            1 2;
            2 5
        ]
    ],
    # document 2
    [
        # modality 1
        [
            3 4;
            4 9
        ],
        # modality 2
        [
            3 4;
            4 6
        ]
    ],
]

@testset "constructor" begin
    model = MultiModalMuSig.MMCTM(K, α, X)
    @test model.D == 2
    @test model.N == [[13, 7], [13, 10]]
    @test model.M == 2
    @test model.V == [4, 4]

    @test length(model.μ) == 5
    @test size(model.Σ) == (5, 5)
    @test size(model.invΣ) == (5, 5)

    @test length(model.ζ) == 2
    @test length(model.ζ[1]) == 2
    @test sum(model.θ[1][1], dims=1) ≈ ones(2)'
    @test length(model.λ[1]) == 5
    @test model.ν[1] == ones(5)

    @test length(model.γ) == 2
    @test length(model.γ[1]) == 2
    @test length(model.γ[1][1]) == 4
    @test length(model.γ[1][2]) == 4
    @test all(model.γ[1][2] .> 0) == true
end

@testset "calculate_Ndivζ" begin
    model = MultiModalMuSig.MMCTM(K, α, X)
    model.ζ = [[2, 3], [4,5]]

    Ndivζ = [
        sum(X[1][1][:, 2]) / model.ζ[1][1],
        sum(X[1][1][:, 2]) / model.ζ[1][1],
        sum(X[1][2][:, 2]) / model.ζ[1][2],
        sum(X[1][2][:, 2]) / model.ζ[1][2],
        sum(X[1][2][:, 2]) / model.ζ[1][2]
    ]

    @test MultiModalMuSig.calculate_Ndivζ(model, 1) ≈ Ndivζ
end

@testset "calculate_sumθ" begin
    model = MultiModalMuSig.MMCTM(K, α, X)
    θ = Matrix{Float64}[]
    push!(θ, [0.4 0.1; 0.6 0.9])
    push!(θ, [0.3 0.4; 0.3 0.5; 0.4 0.1])
    model.θ[1] = θ

    sumθ = [
        X[1][1][1, 2] * θ[1][1, 1] + X[1][1][2, 2] * θ[1][1, 2],
        X[1][1][1, 2] * θ[1][2, 1] + X[1][1][2, 2] * θ[1][2, 2],
        X[1][2][1, 2] * θ[2][1, 1] + X[1][2][2, 2] * θ[2][1, 2],
        X[1][2][1, 2] * θ[2][2, 1] + X[1][2][2, 2] * θ[2][2, 2],
        X[1][2][1, 2] * θ[2][3, 1] + X[1][2][2, 2] * θ[2][3, 2],
    ]

    @test MultiModalMuSig.calculate_sumθ(model, 1) ≈ sumθ
end

@testset "λ_objective" begin
    model = MultiModalMuSig.MMCTM(K, α, X)

    λ = Float64[1, 2, 3, 4, 1]
    model.λ[1] = copy(λ)

    MultiModalMuSig.update_λ!(model, 1)
    @test model.λ[1] ≉ λ
    @test isnan.(model.λ[1]) == fill(false, sum(K))
end

function calc_ν_L(invΣ, λ, ν, ζ, X)
    L = (
        -0.5 * tr(diagm(0 => ν) * invΣ) -
        sum(X[1][1][:, 2]) / ζ[1] *
            (exp(λ[1] + 0.5ν[1]) + exp(λ[2] + 0.5ν[2])) -
        sum(X[1][2][:, 2]) / ζ[2] *
            (exp(λ[3] + 0.5ν[3]) + exp(λ[4] + 0.5ν[4]) + exp(λ[5] + 0.5ν[5])) +
        0.5 * (log(ν[1]) + log(ν[2]) + log(ν[3]) + log(ν[4]) + log(ν[5]))
    )
    return L
end

function calc_ν_grad(invΣ, λ, ν, ζ, X)
    grad = -0.5 * diag(invΣ)
    grad[1] += (
        -sum(X[1][1][:, 2]) / (2 * ζ[1]) * exp(λ[1] + 0.5ν[1]) + 1 / (2ν[1])
    )
    grad[2] += (
        -sum(X[1][1][:, 2]) / (2 * ζ[1]) * exp(λ[2] + 0.5ν[2]) + 1 / (2ν[2])
    )
    grad[3] += (
        -sum(X[1][2][:, 2]) / (2 * ζ[2]) * exp(λ[3] + 0.5ν[3]) + 1 / (2ν[3])
    )
    grad[4] += (
        -sum(X[1][2][:, 2]) / (2 * ζ[2]) * exp(λ[4] + 0.5ν[4]) + 1 / (2ν[4])
    )
    grad[5] += (
        -sum(X[1][2][:, 2]) / (2 * ζ[2]) * exp(λ[5] + 0.5ν[5]) + 1 / (2ν[5])
    )
    return grad
end

@testset "ν_objective" begin
    model = MultiModalMuSig.MMCTM(K, α, X)

    μ = Float64[1, 1, 2, 2, 1]
    invΣ = Matrix{Float64}(I, sum(K), sum(K))
    λ = Float64[1, 2, 3, 4, 1]
    ν = Float64[1, 1, 1, 2, 1]
    ζ = Float64[2, 1]

    ∇ν = Array{Float64}(undef, sum(K))
    Ndivζ = vcat([fill(model.N[1][m] / ζ[m], K[m]) for m in 1:model.M]...)
    objective = MultiModalMuSig.ν_objective(ν, ∇ν, λ, Ndivζ, μ, invΣ)
    @test objective ≈ calc_ν_L(invΣ, λ, ν, ζ, X)
    @test ∇ν ≈ calc_ν_grad(invΣ, λ, ν, ζ, X)

    model.μ = μ
    model.λ[1] = λ
    model.ν[1] = ν
    model.ζ[1] = ζ
    MultiModalMuSig.update_ν!(model, 1)
    @test (model.ν[1] .> 0.0) == fill(true, sum(K))
end

@testset "update_ζ!" begin
    model = MultiModalMuSig.MMCTM(K, α, X)
    model.λ = [[1,2,3,4,1], [2,3,1,4,2]]
    model.ν = [[1,1,1,2,1], [1,3,1,2,1]]
    MultiModalMuSig.update_ζ!(model, 1)

    ζ = Float64[exp(1.5) + exp(2.5), exp(3.5) + exp(5) + exp(1.5)]
    @test model.ζ[1] ≈ ζ
end

@testset "update_θ!" begin
    model = MultiModalMuSig.MMCTM(K, α, X)
    model.λ = [[1,2,3,4,1], [2,3,1,4,2]]
    model.γ = [
        [
            [1, 2, 2, 6],
            [2, 3, 1, 2]
        ],
        [
            [1, 2, 3, 4],
            [2, 1, 2, 6],
            [1, 1, 3, 1]
        ]
    ]
    MultiModalMuSig.update_Elnϕ!(model)
    MultiModalMuSig.update_θ!(model, 1)

    θ = Array{Float64}(undef, 2, 2)
    θ[1, 1] = exp(1 + digamma(1) - digamma(11))
    θ[2, 1] = exp(2 + digamma(2) - digamma(8))
    θ[1, 2] = exp(1 + digamma(2) - digamma(11))
    θ[2, 2] = exp(2 + digamma(3) - digamma(8))
    θ[:, 1] ./= sum(θ[:, 1])
    θ[:, 2] ./= sum(θ[:, 2])

    @test sum(model.θ[1][1], dims=1) ≈ ones(size(X[1][1])[1])'
    @test model.θ[1][1] ≈ θ
    @test any(model.θ[1][1] .< 0.0) == false

    MultiModalMuSig.update_θ!(model, 2)
    θ = Array{Float64}(undef, 3, 2)
    θ[1, 1] = exp(1 + digamma(3) - digamma(10))
    θ[2, 1] = exp(4 + digamma(2) - digamma(11))
    θ[3, 1] = exp(2 + digamma(3) - digamma(6))
    θ[1, 2] = exp(1 + digamma(4) - digamma(10))
    θ[2, 2] = exp(4 + digamma(6) - digamma(11))
    θ[3, 2] = exp(2 + digamma(1) - digamma(6))
    θ[:, 1] ./= sum(θ[:, 1])
    θ[:, 2] ./= sum(θ[:, 2])

    @test model.θ[2][2] ≈ θ
end

@testset "update_μ!" begin
    model = MultiModalMuSig.MMCTM(K, α, X)
    model.λ = [[1,2,3,4,1], [2,3,1,4,2]]

    MultiModalMuSig.update_μ!(model)

    @test model.μ ≈ [1.5, 2.5, 2.0, 4.0, 1.5]
end

@testset "update_Σ!" begin
    model = MultiModalMuSig.MMCTM(K, α, X)
    model.λ = [[1,2,3,4,1], [2,3,1,4,2]]
    model.ν = [[1,1,1,2,1], [1,3,1,2,1]]
    model.μ = [1, 1, 2, 2, 1]

    MultiModalMuSig.update_Σ!(model)

    diff1 = model.λ[1] .- model.μ
    diff2 = model.λ[2] .- model.μ
    Σ = 0.5 * (
        diagm(0 => model.ν[1]) + diagm(0 => model.ν[2]) +
        (diff1 * diff1') + (diff2 * diff2')
    )
    @test model.Σ ≈ Σ
    @test model.invΣ ≈ inv(Σ)
end

@testset "update_γ!" begin
    model = MultiModalMuSig.MMCTM(K, α, X)
    model.θ[1][1] = [0.4 0.1; 0.6 0.9]
    model.θ[2][1] = [0.3 0.5; 0.7 0.5]
    model.θ[1][2] = [0.2 0.6; 0.7 0.3; 0.1 0.1]
    model.θ[2][2] = [0.1 0.3; 0.7 0.5; 0.2 0.2]
    MultiModalMuSig.update_γ!(model)

    γ1 = [0.1 + 5 * 0.4, 0.1 + 8 * 0.1, 0.1 + 4 * 0.3, 0.1 + 9 * 0.5]
    γ2 = [0.1 + 5 * 0.6, 0.1 + 8 * 0.9, 0.1 + 4 * 0.7, 0.1 + 9 * 0.5]
    @test model.γ[1][1] ≈ γ1
    @test model.γ[1][2] ≈ γ2

    γ1 = [0.1 + 2 * 0.2, 0.1 + 5 * 0.6, 0.1 + 4 * 0.1, 0.1 + 6 * 0.3]
    γ2 = [0.1 + 2 * 0.7, 0.1 + 5 * 0.3, 0.1 + 4 * 0.7, 0.1 + 6 * 0.5]
    γ3 = [0.1 + 2 * 0.1, 0.1 + 5 * 0.1, 0.1 + 4 * 0.2, 0.1 + 6 * 0.2]
    @test model.γ[2][1] ≈ γ1
    @test model.γ[2][2] ≈ γ2
    @test model.γ[2][3] ≈ γ3
end

@testset "update_Elnϕ!" begin
    model = MultiModalMuSig.MMCTM(K, α, X)
    model.γ[1][1] .= [1, 2, 1, 3]

    MultiModalMuSig.update_Elnϕ!(model)

    @test model.Elnϕ[1][1][1] ≈ digamma(1) - digamma(7)
end

@testset "update_α!" begin
    model = MultiModalMuSig.MMCTM(K, α, X)

    sum_Elnϕ = sum(model.Elnϕ[1][1]) + sum(model.Elnϕ[1][2])
    L = K[1] * (lgamma(4α[1]) - 4lgamma(α[1])) + α[1] * sum_Elnϕ
    grad = 4K[1] * (digamma(4α[1]) - digamma(α[1])) + sum_Elnϕ

    ∇α = Array{Float64}(undef, 1)
    @test MultiModalMuSig.α_objective(model.α[1:1], ∇α, sum_Elnϕ, K[1], 4) ≈ L
    @test ∇α[1] ≈ grad

    sum_Elnϕ = (
        sum(model.Elnϕ[2][1]) + sum(model.Elnϕ[2][2]) + sum(model.Elnϕ[2][3])
    )
    L_before = (
        K[2] * (lgamma(4model.α[2]) - 4lgamma(model.α[2])) +
        model.α[2] * sum_Elnϕ
    )
    MultiModalMuSig.update_α!(model)
    L_after = (
        K[2] * (lgamma(4model.α[2]) - 4lgamma(model.α[2])) +
        model.α[2] * sum_Elnϕ
    )
    @test model.α ≉ α
    @test L_after > L_before
end

#@testset "calculate_ElnPϕ" begin
    #model = MultiModalMuSig.MMCTM(K, α, X)
    #ElnPϕ = 0.0
    #@pending @test MultiModalMuSig.calculate_ElnPϕ(model) ≈ ElnPϕ
#end

#@testset "calculate_ElnPη" begin
    #model = MultiModalMuSig.MMCTM(K, α, X)
    #ElnPη = 0.0
    #@pending @test MultiModalMuSig.calculate_ElnPη(model) ≈ ElnPη
#end

#@testset "calculate_ElnPZ" begin
    #model = MultiModalMuSig.MMCTM(K, α, X)
    #ElnPZ = 0.0
    #@pending @test MultiModalMuSig.calculate_ElnPZ(model) ≈ ElnPZ
#end

#@testset "calculate_ElnPX" begin
    #model = MultiModalMuSig.MMCTM(K, α, X)
    #ElnPX = 0.0
    #@pending @test MultiModalMuSig.calculate_ElnPX(model) ≈ ElnPX
#end

#@testset "calculate_ElnQϕ" begin
    #model = MultiModalMuSig.MMCTM(K, α, X)
    #ElnQϕ = 0.0
    #@pending @test MultiModalMuSig.calculate_ElnQϕ(model) ≈ ElnQϕ
#end

#@testset "calculate_ElnQη" begin
    #model = MultiModalMuSig.MMCTM(K, α, X)
    #ElnQη = 0.0
    #@pending @test MultiModalMuSig.calculate_ElnQη(model) ≈ ElnQη
#end

#@testset "calculate_ElnQZ" begin
    #model = MultiModalMuSig.MMCTM(K, α, X)
    #ElnQZ = 0.0
    #@pending @test MultiModalMuSig.calculate_ElnQZ(model) ≈ ElnQZ
#end

@testset "calculate_elbo" begin
    model = MultiModalMuSig.MMCTM(K, α, X)
    @test MultiModalMuSig.calculate_elbo(model) <= 0.0
end

@testset "fit" begin
    model = MultiModalMuSig.MMCTM(K, α, X)
    ll = MultiModalMuSig.fit!(model, maxiter=1, verbose=false)
    @test length(ll) == 1
    @test length(ll[1]) == 2
end

@testset "loglikelihoods" begin
    η = [[1.0, 2.0], [2.0, 3.0]]
    props = [exp.(η[d]) ./ sum(exp.(η[d])) for d in 1:2]

    γ = [
        [1, 2, 1, 3],
        [1, 1, 2, 4]
    ]
    ϕ = [γ[k] ./ sum(γ[k]) for k in 1:2]

    Xm1 = [X[d][1] for d in 1:2]

    sum_ll = [
        (
            Xm1[1][1, 2] * log(props[1][1] * ϕ[1][1] + props[1][2] * ϕ[2][1]) +
            Xm1[1][2, 2] * log(props[1][1] * ϕ[1][2] + props[1][2] * ϕ[2][2])
        ),
        (
            Xm1[2][1, 2] * log(props[2][1] * ϕ[1][3] + props[2][2] * ϕ[2][3]) +
            Xm1[2][2, 2] * log(props[2][1] * ϕ[1][4] + props[2][2] * ϕ[2][4])
        )
    ]
    N = [sum(X[d][1][:, 2]) for d in 1:2]

    res = MultiModalMuSig.calculate_docmodality_loglikelihood(
        Xm1[1], props[1], ϕ
    )
    @test res ≈ sum_ll[1] / N[1]

    res = MultiModalMuSig.calculate_modality_loglikelihood(Xm1, props, ϕ)
    @test res ≈ sum(sum_ll) / sum(N)

    model = MultiModalMuSig.MMCTM(K, α, X)
    model.props[1][1] = props[1]
    model.props[2][1] = props[2]
    model.ϕ[1] = deepcopy(ϕ)

    res = MultiModalMuSig.calculate_loglikelihoods(X, model)
    @test res[1] ≈ sum(sum_ll) / sum(N)
end


text:
@doc """
A class defining a 2D elliptical gaussian
""" ->
type Gaussian
	# All of the parameters
	mass::Float64
	q::Float64
	xc::Float64
	yc::Float64
	width::Float64
	theta::Float64

	# Derived parameters
	qinv::Float64
	width2inv::Float64
	sin_theta::Float64
	cos_theta::Float64
	coeff::Float64
end

@doc """
Constructor taking a vector of five parameters
""" ->
function Gaussian(params::Vector{Float64})
	gaussian = Gaussian(params[1], params[2], params[3], params[4], params[5],
							params[6], 0.0, 0.0, 0.0, 0.0, 0.0)
	calculate_derived!(gaussian)
	return gaussian
end

@doc """
Calculate the derived parameters
""" ->
function calculate_derived!(gaussian::Gaussian)
	gaussian.qinv = 1.0/gaussian.q
	gaussian.width2inv = 1.0/gaussian.width^2
	gaussian.sin_theta = sin(gaussian.theta)
	gaussian.cos_theta = cos(gaussian.theta)
	gaussian.coeff = gaussian.mass*gaussian.width2inv/(2*pi)
end

@doc """
Evaluate the gaussian at the position (x, y)
""" ->
function evaluate(gaussian::Gaussian, x::Float64, y::Float64)
	xx =  (x - gaussian.xc)*gaussian.cos_theta +
			(y - gaussian.yc)*gaussian.sin_theta
	yy = -(x - gaussian.xc)*gaussian.sin_theta +
			(y - gaussian.yc)*gaussian.cos_theta
	rsq = gaussian.q*xx^2 + yy^2*gaussian.qinv
	f = gaussian.coeff*exp(-0.5*rsq*gaussian.width2inv)
	return f
end

max_stars_repo_name: eggplantbren/MogTrack

max_issues_repo_name: eggplantbren/MogTrack

max_forks_repo_name: eggplantbren/MogTrack


text:
using Test, LazySets, Reachability, MathematicalSystems, MathematicalPredicates,
      LinearAlgebra, SparseArrays, IntervalMatrices

# fix namespace conflicts with MathematicalSystems
using LazySets: LinearMap, AffineMap, ResetMap, Interval

include("Systems/alltests.jl")
include("ReachSets/alltests.jl")
include("Reachability/alltests.jl")
include("Properties/alltests.jl")


text:
module JuliaAcademyData

import Pkg

export datapath, activate

const COURSE = Ref("")

"""
    datapath(paths...)

Like [`joinpath`](@ref), but specifically for loading supplimentary data and scripts for JuliaAcademy notebooks.

The course must first be [`activate`](@ref)'d before calling this function.
"""
function datapath(paths...; course=COURSE[])
    isempty(course) && error("no course has been activated; be sure to execute the first cell of the notebook to initialize the lecture")
    return joinpath(dirname(@__DIR__), "courses", course, paths...)
end

"""
    activate(course)

Activate the Pkg environment for a JuliaAcademy course and set up the package for loading data for that course.
"""
function activate(course)
    COURSE[] = course
    isdir(datapath()) || error("course $course does not exist")
    Pkg.activate(datapath())
    Pkg.instantiate()
    return nothing
end

end # module


text:
module TestShow
    using Compat, Compat.Test, DataFrames

    # In the future newline character \n should be added to this test case
    df = DataFrame(A = Int64[1:4;], B = ["x\"", "∀ε>0: x+ε>x", "z\$", "ABC"],
                   C = Float32[1.0, 2.0, 3.0, 4.0])

    refstr = """
    4×3 $DataFrame
    │ Row │ A │ B           │ C   │
    ├─────┼───┼─────────────┼─────┤
    │ 1   │ 1 │ x\"          │ 1.0 │
    │ 2   │ 2 │ ∀ε>0: x+ε>x │ 2.0 │
    │ 3   │ 3 │ z\$          │ 3.0 │
    │ 4   │ 4 │ ABC         │ 4.0 │"""

    for f in [show, showall], allcols in [true, false]
        io = IOBuffer()
        f(io, df, allcols)
        str = String(take!(io))
        @test str == refstr
    end

    refstr = """
    4×3 $DataFrame
    │ Col # │ Name │ Eltype  │ Missing │
    ├───────┼──────┼─────────┼─────────┤
    │ 1     │ A    │ Int64   │ 0       │
    │ 2     │ B    │ String  │ 0       │
    │ 3     │ C    │ Float32 │ 0       │"""
    for a in [true, false]
        io = IOBuffer()
        showcols(io, df, a, false)
        str = String(take!(io))
        @test str == refstr
    end

    refstr = """
    4×3 $DataFrame
    │ Col # │ Name │ Eltype  │ Missing │ Values      │
    ├───────┼──────┼─────────┼─────────┼─────────────┤
    │ 1     │ A    │ Int64   │ 0       │ 1  …  4     │
    │ 2     │ B    │ String  │ 0       │ x\"  …  ABC  │
    │ 3     │ C    │ Float32 │ 0       │ 1.0  …  4.0 │"""
    for a in [true, false]
        io = IOBuffer()
        showcols(io, df, a, true)
        str = String(take!(io))
        @test str == refstr
    end

    srand(1)
    df_big = DataFrame(rand(25,5))

    io = IOContext(IOBuffer(), :displaysize=>(10,40))
    show(io, df_big)
    str = String(take!(io.io))
    @test str == """
    25×5 $DataFrame. Omitted printing of 2 columns
    │ Row │ x1       │ x2       │ x3       │
    ├─────┼──────────┼──────────┼──────────┤
    │ 1   │ 0.236033 │ 0.644883 │ 0.440897 │
    ⋮
    │ 24  │ 0.278582 │ 0.241591 │ 0.990741 │
    │ 25  │ 0.751313 │ 0.884837 │ 0.550334 │"""

    io = IOContext(IOBuffer(), :displaysize=>(10,40))
    show(io, df_big, true)
    str = String(take!(io.io))
    @test str == """
    25×5 $DataFrame
    │ Row │ x1       │ x2       │ x3       │
    ├─────┼──────────┼──────────┼──────────┤
    │ 1   │ 0.236033 │ 0.644883 │ 0.440897 │
    ⋮
    │ 24  │ 0.278582 │ 0.241591 │ 0.990741 │
    │ 25  │ 0.751313 │ 0.884837 │ 0.550334 │

    │ Row │ x4       │ x5       │
    ├─────┼──────────┼──────────┤
    │ 1   │ 0.580782 │ 0.138763 │
    ⋮
    │ 24  │ 0.762276 │ 0.755415 │
    │ 25  │ 0.339081 │ 0.649056 │"""

    io = IOContext(IOBuffer(), :displaysize=>(10,40))
    showall(io, df_big)
    str = String(take!(io.io))
    @test str == """
    25×5 $DataFrame
    │ Row │ x1         │ x2        │ x3        │ x4        │ x5        │
    ├─────┼────────────┼───────────┼───────────┼───────────┼───────────┤
    │ 1   │ 0.236033   │ 0.644883  │ 0.440897  │ 0.580782  │ 0.138763  │
    │ 2   │ 0.346517   │ 0.0778264 │ 0.404673  │ 0.768359  │ 0.456446  │
    │ 3   │ 0.312707   │ 0.848185  │ 0.736787  │ 0.519525  │ 0.739918  │
    │ 4   │ 0.00790928 │ 0.0856352 │ 0.953803  │ 0.514863  │ 0.816004  │
    │ 5   │ 0.488613   │ 0.553206  │ 0.0951856 │ 0.998136  │ 0.114529  │
    │ 6   │ 0.210968   │ 0.46335   │ 0.519675  │ 0.603682  │ 0.748928  │
    │ 7   │ 0.951916   │ 0.185821  │ 0.0135403 │ 0.758775  │ 0.878108  │
    │ 8   │ 0.999905   │ 0.111981  │ 0.303399  │ 0.590953  │ 0.930481  │
    │ 9   │ 0.251662   │ 0.976312  │ 0.702557  │ 0.722086  │ 0.896291  │
    │ 10  │ 0.986666   │ 0.0516146 │ 0.596537  │ 0.953207  │ 0.663145  │
    │ 11  │ 0.555751   │ 0.53803   │ 0.638935  │ 0.384411  │ 0.472799  │
    │ 12  │ 0.437108   │ 0.455692  │ 0.872347  │ 0.320011  │ 0.880525  │
    │ 13  │ 0.424718   │ 0.279395  │ 0.548635  │ 0.865625  │ 0.0141033 │
    │ 14  │ 0.773223   │ 0.178246  │ 0.262992  │ 0.45457   │ 0.502774  │
    │ 15  │ 0.28119    │ 0.548983  │ 0.526443  │ 0.420287  │ 0.224851  │
    │ 16  │ 0.209472   │ 0.370971  │ 0.465019  │ 0.225151  │ 0.287858  │
    │ 17  │ 0.251379   │ 0.894166  │ 0.275519  │ 0.286169  │ 0.104033  │
    │ 18  │ 0.0203749  │ 0.648054  │ 0.461823  │ 0.309144  │ 0.475749  │
    │ 19  │ 0.287702   │ 0.417039  │ 0.951861  │ 0.170391  │ 0.416681  │
    │ 20  │ 0.859512   │ 0.144566  │ 0.288737  │ 0.147162  │ 0.521387  │
    │ 21  │ 0.0769509  │ 0.622403  │ 0.661232  │ 0.230063  │ 0.908499  │
    │ 22  │ 0.640396   │ 0.872334  │ 0.194568  │ 0.0929292 │ 0.102832  │
    │ 23  │ 0.873544   │ 0.524975  │ 0.393193  │ 0.681415  │ 0.670421  │
    │ 24  │ 0.278582   │ 0.241591  │ 0.990741  │ 0.762276  │ 0.755415  │
    │ 25  │ 0.751313   │ 0.884837  │ 0.550334  │ 0.339081  │ 0.649056  │"""

    io = IOContext(IOBuffer(), :displaysize=>(10,40))
    showall(io, df_big, false)
    str = String(take!(io.io))
    @test str == """
    25×5 $DataFrame. Omitted printing of 3 columns
    │ Row │ x1         │ x2        │
    ├─────┼────────────┼───────────┤
    │ 1   │ 0.236033   │ 0.644883  │
    │ 2   │ 0.346517   │ 0.0778264 │
    │ 3   │ 0.312707   │ 0.848185  │
    │ 4   │ 0.00790928 │ 0.0856352 │
    │ 5   │ 0.488613   │ 0.553206  │
    │ 6   │ 0.210968   │ 0.46335   │
    │ 7   │ 0.951916   │ 0.185821  │
    │ 8   │ 0.999905   │ 0.111981  │
    │ 9   │ 0.251662   │ 0.976312  │
    │ 10  │ 0.986666   │ 0.0516146 │
    │ 11  │ 0.555751   │ 0.53803   │
    │ 12  │ 0.437108   │ 0.455692  │
    │ 13  │ 0.424718   │ 0.279395  │
    │ 14  │ 0.773223   │ 0.178246  │
    │ 15  │ 0.28119    │ 0.548983  │
    │ 16  │ 0.209472   │ 0.370971  │
    │ 17  │ 0.251379   │ 0.894166  │
    │ 18  │ 0.0203749  │ 0.648054  │
    │ 19  │ 0.287702   │ 0.417039  │
    │ 20  │ 0.859512   │ 0.144566  │
    │ 21  │ 0.0769509  │ 0.622403  │
    │ 22  │ 0.640396   │ 0.872334  │
    │ 23  │ 0.873544   │ 0.524975  │
    │ 24  │ 0.278582   │ 0.241591  │
    │ 25  │ 0.751313   │ 0.884837  │"""

    io = IOContext(IOBuffer(), :displaysize=>(10,40))
    showcols(io, df_big, false, false)
    str = String(take!(io.io))
    @test str == """
    25×5 $DataFrame
    │ Col # │ Name │ Eltype  │ Missing │
    ├───────┼──────┼─────────┼─────────┤
    │ 1     │ x1   │ Float64 │ 0       │
    │ 2     │ x2   │ Float64 │ 0       │
    │ 3     │ x3   │ Float64 │ 0       │
    │ 4     │ x4   │ Float64 │ 0       │
    │ 5     │ x5   │ Float64 │ 0       │"""

    io = IOContext(IOBuffer(), :displaysize=>(10,40))
    showcols(io, df_big, true, false)
    str = String(take!(io.io))
    @test str == """
    25×5 $DataFrame
    │ Col # │ Name │ Eltype  │ Missing │
    ├───────┼──────┼─────────┼─────────┤
    │ 1     │ x1   │ Float64 │ 0       │
    │ 2     │ x2   │ Float64 │ 0       │
    │ 3     │ x3   │ Float64 │ 0       │
    │ 4     │ x4   │ Float64 │ 0       │
    │ 5     │ x5   │ Float64 │ 0       │"""

    io = IOContext(IOBuffer(), :displaysize=>(10,40))
    showcols(io, df_big, false, true)
    str = String(take!(io.io))
    @test str == """
    25×5 $DataFrame
    │ Col # │ Name │ Eltype  │ Missing │
    ├───────┼──────┼─────────┼─────────┤
    │ 1     │ x1   │ Float64 │ 0       │
    │ 2     │ x2   │ Float64 │ 0       │
    │ 3     │ x3   │ Float64 │ 0       │
    │ 4     │ x4   │ Float64 │ 0       │
    │ 5     │ x5   │ Float64 │ 0       │

    │ Col # │ Values                │
    ├───────┼───────────────────────┤
    │ 1     │ 0.236033  …  0.751313 │
    │ 2     │ 0.644883  …  0.884837 │
    │ 3     │ 0.440897  …  0.550334 │
    │ 4     │ 0.580782  …  0.339081 │
    │ 5     │ 0.138763  …  0.649056 │"""

    io = IOContext(IOBuffer(), :displaysize=>(10,40))
    showcols(io, df_big, true, true)
    str = String(take!(io.io))
    @test str == """
    25×5 $DataFrame
    │ Col # │ Name │ Eltype  │ Missing │ Values                │
    ├───────┼──────┼─────────┼─────────┼───────────────────────┤
    │ 1     │ x1   │ Float64 │ 0       │ 0.236033  …  0.751313 │
    │ 2     │ x2   │ Float64 │ 0       │ 0.644883  …  0.884837 │
    │ 3     │ x3   │ Float64 │ 0       │ 0.440897  …  0.550334 │
    │ 4     │ x4   │ Float64 │ 0       │ 0.580782  …  0.339081 │
    │ 5     │ x5   │ Float64 │ 0       │ 0.138763  …  0.649056 │"""

    io = IOBuffer()
    df_small = DataFrame(transpose([1.0:5.0;]))
    showcols(io, df_small)
    str = String(take!(io))
    @test str == """
    1×5 $DataFrame
    │ Col # │ Name │ Eltype  │ Missing │ Values │
    ├───────┼──────┼─────────┼─────────┼────────┤
    │ 1     │ x1   │ Float64 │ 0       │ 1.0    │
    │ 2     │ x2   │ Float64 │ 0       │ 2.0    │
    │ 3     │ x3   │ Float64 │ 0       │ 3.0    │
    │ 4     │ x4   │ Float64 │ 0       │ 4.0    │
    │ 5     │ x5   │ Float64 │ 0       │ 5.0    │"""

    io = IOBuffer()
    df_min = DataFrame(rand(0,5))
    showcols(io, df_min)
    str = String(take!(io))
    @test str == """
    0×5 $DataFrame
    │ Col # │ Name │ Eltype  │ Missing │
    ├───────┼──────┼─────────┼─────────┤
    │ 1     │ x1   │ Float64 │ 0       │
    │ 2     │ x2   │ Float64 │ 0       │
    │ 3     │ x3   │ Float64 │ 0       │
    │ 4     │ x4   │ Float64 │ 0       │
    │ 5     │ x5   │ Float64 │ 0       │"""

    subdf = view(df, [2, 3]) # df[df[:A] .> 1.0, :]
    show(io, subdf)
    show(io, subdf, true)
    showall(io, subdf)
    showall(io, subdf, false)

    dfvec = DataFrame[df for _=1:3]
    show(io, dfvec)
    showall(io, dfvec)

    gd = groupby(df, :A)
    show(io, gd)
    showall(io, gd)

    dfr = DataFrameRow(df, 1)
    show(io, dfr)

    df = DataFrame(A = Vector{String}(3))

    A = DataFrames.StackedVector(Any[[1, 2, 3], [4, 5, 6], [7, 8, 9]])
    show(io, A)
    A = DataFrames.RepeatedVector([1, 2, 3], 5, 1)
    show(io, A)
    A = DataFrames.RepeatedVector([1, 2, 3], 1, 5)
    show(io, A)

    #Test show output for REPL and similar
    df = DataFrame(Fish = ["Suzy", "Amir"], Mass = [1.5, missing])
    if VERSION >= v"0.7.0-DEV.3077"
        @test sprint(show, df, context=:color=>true) == """
        2×2 DataFrame
        │ Row │ Fish │ Mass    │
        ├─────┼──────┼─────────┤
        │ 1   │ Suzy │ 1.5     │
        │ 2   │ Amir │ \e[90mmissing\e[39m │"""
    else
        io = IOBuffer()
        show(io, df)
        str = String(take!(io))
        @test str == (Base.have_color ? """
        2×2 $DataFrame
        │ Row │ Fish │ Mass    │
        ├─────┼──────┼─────────┤
        │ 1   │ Suzy │ 1.5     │
        │ 2   │ Amir │ \e[90mmissing\e[39m │""" : """
        2×2 $DataFrame
        │ Row │ Fish │ Mass    │
        ├─────┼──────┼─────────┤
        │ 1   │ Suzy │ 1.5     │
        │ 2   │ Amir │ missing │""")
    end

    io = IOBuffer()
    showcols(io, df)
    str = String(take!(io))
    if VERSION ≥ v"0.7.0-DEV.2762"
        ref = """
        2×2 $DataFrame
        │ Col # │ Name │ Eltype                  │ Missing │ Values          │
        ├───────┼──────┼─────────────────────────┼─────────┼─────────────────┤
        │ 1     │ Fish │ String                  │ 0       │ Suzy  …  Amir   │
        │ 2     │ Mass │ $(Union{Missing,Float64}) │ 1       │ 1.5  …  missing │"""
    else
        ref = """
        2×2 $DataFrame
        │ Col # │ Name │ Eltype                           │ Missing │ Values          │
        ├───────┼──────┼──────────────────────────────────┼─────────┼─────────────────┤
        │ 1     │ Fish │ String                           │ 0       │ Suzy  …  Amir   │
        │ 2     │ Mass │ Union{Float64, $Missing} │ 1       │ 1.5  …  missing │"""
    end
    @test str == ref

    # Test showing missing
    df = DataFrame(A = [:Symbol, missing, :missing],
                   B = [missing, "String", "missing"],
                   C = [:missing, "missing", missing])
    if VERSION ≥ v"0.7.0-DEV.2762"
        @test sprint(show, df, context=:color=>true) == """
        3×3 DataFrame
        │ Row │ A       │ B       │ C       │
        ├─────┼─────────┼─────────┼─────────┤
        │ 1   │ Symbol  │ \e[90mmissing\e[39m │ missing │
        │ 2   │ \e[90mmissing\e[39m │ String  │ missing │
        │ 3   │ missing │ missing │ \e[90mmissing\e[39m │"""
    else
        io = IOBuffer()
        show(io, df)
        str = String(take!(io))
        @test str == (Base.have_color ? """
        3×3 $DataFrame
        │ Row │ A       │ B       │ C       │
        ├─────┼─────────┼─────────┼─────────┤
        │ 1   │ Symbol  │ \e[90mmissing\e[39m │ missing │
        │ 2   │ \e[90mmissing\e[39m │ String  │ missing │
        │ 3   │ missing │ missing │ \e[90mmissing\e[39m │""" : """
        3×3 $DataFrame
        │ Row │ A       │ B       │ C       │
        ├─────┼─────────┼─────────┼─────────┤
        │ 1   │ Symbol  │ missing │ missing │
        │ 2   │ missing │ String  │ missing │
        │ 3   │ missing │ missing │ missing │""")
    end

    io = IOBuffer()
    showcols(io, df)
    str = String(take!(io))
    if VERSION ≥ v"0.7.0-DEV.2762"
        ref = """
        3×3 DataFrame
        │ Col # │ Name │ Eltype                 │ Missing │ Values              │
        ├───────┼──────┼────────────────────────┼─────────┼─────────────────────┤
        │ 1     │ A    │ Union{Missing, Symbol} │ 1       │ Symbol  …  missing  │
        │ 2     │ B    │ Union{Missing, String} │ 1       │ missing  …  missing │
        │ 3     │ C    │ Any                    │ 1       │ missing  …  missing │"""
    elseif VERSION ≥ v"0.7.0-DEV.2762"
        ref = """
        3×3 $DataFrame
        │ Col # │ Name │ Eltype                 │ Missing │
        ├───────┼──────┼────────────────────────┼─────────┤
        │ 1     │ A    │ $(Union{Missing,Symbol}) │ 1       │
        │ 2     │ B    │ $(Union{Missing,String}) │ 1       │
        │ 3     │ C    │ Any                    │ 1       │

        │ Col # │ Values              │
        ├───────┼─────────────────────┤
        │ 1     │ Symbol  …  missing  │
        │ 2     │ missing  …  missing │
        │ 3     │ missing  …  missing │"""
    else
        ref = """
        3×3 $DataFrame
        │ Col # │ Name │ Eltype                          │ Missing │
        ├───────┼──────┼─────────────────────────────────┼─────────┤
        │ 1     │ A    │ $(Union{Missing,Symbol}) │ 1       │
        │ 2     │ B    │ $(Union{Missing,String}) │ 1       │
        │ 3     │ C    │ Any                             │ 1       │

        │ Col # │ Values              │
        ├───────┼─────────────────────┤
        │ 1     │ Symbol  …  missing  │
        │ 2     │ missing  …  missing │
        │ 3     │ missing  …  missing │"""
    end
    @test str == ref

    # Test computing width for Array{String} columns
    df = DataFrame(Any[["a"]], [:x])
    io = IOBuffer()
    show(io, df)
    str = String(take!(io))
    @test str == """
    1×1 $DataFrame
    │ Row │ x │
    ├─────┼───┤
    │ 1   │ a │"""

    # Test escape characters
    df = DataFrame(a = ["1\n1", "2\t2", "3\r3", "4\$4", "5\"5", "6\\6"])
    io = IOBuffer()
    show(io, df)
    str = String(take!(io))
    @test str == """
    6×1 $DataFrame
    │ Row │ a    │
    ├─────┼──────┤
    │ 1   │ 1\\n1 │
    │ 2   │ 2\\t2 │
    │ 3   │ 3\\r3 │
    │ 4   │ 4\$4  │
    │ 5   │ 5\"5  │
    │ 6   │ 6\\\\6 │"""

end


text:
using HypothesisTests, Base.Test

@test SWCoeffs(3).A == [sqrt(2.0)/2.0]
@test SWCoeffs(3).N == 3

S = SWCoeffs(10)
@test S.N == 10
@test length(S) == 10
@test endof(S) == 10
@test length(S.A) == 5

@test S.A[4] == S[4] == -S[7]
@test S.A[2] == S[2] == -S[9]

S = SWCoeffs(11)
@test S.N == 11
@test length(S) == 11
@test endof(S) == 11
@test length(S.A) == 5

@test S[6] == 0.0
@test S.A[5] == S[5] == -S[7]
@test S.A[3] == S[3] == -S[9]
@test S.A[1] == S[1] == -S[11]

#anti-symmetry
S = SWCoeffs(20)
@test all([S[i] == -S[end-i+1] for i in 1:length(S)])

# Values obtained by calling `_swilk` fortran subroutine directly.

SWc10 = SWCoeffs(10)
a = SWc10.A .- [0.573715, 0.32897, 0.214349, 0.122791, 0.0400887]
@test norm(a,1) ≈ 0.0 atol=S.N*eps(Float32)

SWc20 = SWCoeffs(20)
b = SWc20.A .- [0.473371, 0.32174, 0.255663, 0.208297, 0.16864, 0.133584, 0.101474, 0.0712893, 0.0423232, 0.0140351]
@test norm(b,1) ≈ 0.0 atol=S.N*eps(Float32)

srand(1)
X = sort(rand(Float64, 10))
@test HypothesisTests.swstat(X, SWc10) ≈ 0.8453033192812377 atol=1.2e-8
@test pvalue(HypothesisTests.swstat(X, SWc10), SWc10) ≈ 0.05106360837817192 atol=4.9e-8

# **Worked Example** (Section 4) from
# PATRICK ROYSTON Approximating the Shapiro-Wilk W-test for non-normality
# *Statistics and Computing* (1992) **2**, 117-119

X = [48.4, 49.0, 59.5, 59.6, 60.7, 88.8, 98.2, 109.4, 169.1, 227.1]
SWc = SWCoeffs(length(X))
@test all(abs.(SWc.A .- [0.5737, 0.3290, 0.2143, 0.1228, 0.0401]) .< 5.0e-5)
W = HypothesisTests.swstat(X, SWc)
@test W ≈ 0.8078 atol=2.9e-5
pval = pvalue(W, SWc)
@test pval ≈ 0.018 atol=4.7e-5

t = ShapiroWilkTest(X)
@test t.W == W
@test pvalue(t) == pval
@test t.N1 == length(X)
@test t.SWc.N == length(X)
@test length(t.SWc) == length(X)

@test pvalue(t) == pval


text:
push!(LOAD_PATH, ".")
using Statistics
using Points

function mean_var(ps :: Array{P}) where {P<:Point}
    m = mean(ps)
    ps2 = map(p->abs2(p-m),ps)
    v = sum(ps2)/(length(ps2)-1)
    return abs(m),v
end

function randomwalk(::Type{P}, num::Integer, nsteps::Integer) where {P<:Point}
    means, vars = zeros(nsteps+1), zeros(nsteps+1)
    ps = zeros(P,num)
    means[1], vars[1] = mean_var(ps)
    for t in 0:nsteps
        map!(update, ps, ps)
        means[t+1], vars[t+1] = mean_var(ps)
    end
    means, vars
end

max_stars_repo_name: isuke/JuliaBook-Samples

max_issues_repo_name: isuke/JuliaBook-Samples

max_forks_repo_name: isuke/JuliaBook-Samples


text:
using BlockBandedMatrices, ArrayLayouts, Test
import BlockBandedMatrices: BandedBlockBandedRowMajor, BandedBlockBandedRows, BandedBlockBandedColumns,BlockBandedRows, BlockBandedColumns

@testset "Adj/Trans" begin
    @testset "BandedBlockBanded" begin
        A = BandedBlockBandedMatrix(randn(ComplexF64,10,14), 1:4,2:5, (1,2), (2,1))

        @test MemoryLayout(transpose(A)) isa BandedBlockBandedRowMajor
        @test MemoryLayout(A') isa BandedBlockBandedRows
        @test MemoryLayout(transpose(A)') isa BandedBlockBandedColumns

        @test A'[Block(1,1)] == A[Block(1,1)]'
        @test A'[Block(2,3)] == A[Block(3,2)]'
        @test transpose(A)[Block(1,1)] == transpose(A[Block(1,1)])
        @test transpose(A)[Block(2,3)] == transpose(A[Block(3,2)])

        @test blockbandwidths(A') == blockbandwidths(transpose(A)) == (2,1)
        @test subblockbandwidths(A') == subblockbandwidths(transpose(A)) == (1,2)

        @test BandedBlockBandedMatrix(A') == A'
        @test BandedBlockBandedMatrix(transpose(A)) == transpose(A)
    end

    @testset "BlockBanded" begin
        A = BlockBandedMatrix(randn(ComplexF64,10,14), 1:4,2:5, (1,2))

        @test MemoryLayout(transpose(A)) isa BlockBandedRows
        @test MemoryLayout(A') isa BlockBandedRows
        @test MemoryLayout(transpose(A)') isa BlockBandedColumns

        @test A'[Block(1,1)] == A[Block(1,1)]'
        @test A'[Block(2,3)] == A[Block(3,2)]'
        @test transpose(A)[Block(1,1)] == transpose(A[Block(1,1)])
        @test transpose(A)[Block(2,3)] == transpose(A[Block(3,2)])

        @test blockbandwidths(A') == blockbandwidths(transpose(A)) == (2,1)

        @test BlockBandedMatrix(A') == A'
        @test BlockBandedMatrix(transpose(A)) == transpose(A)
    end
end


text:
function soi_distance(t, orb::Orbit, bd::CelestialBody)
    if orb.primary != bd.orbit.primary
        throw(ArgumentError("orbit and body do not share a primary body"))
    end
    return norm(time_orbital_position(t, orb) - time_orbital_position(t, bd.orbit)) - bd.SoI
end

"""
    nextorb = propagate_across_soi(t, orb)

Searches for SoI changes starting at time `t`, and returns the orbit after the SoI
patch if one occurs.

For periodic (elliptical) orbits, the search begins at time `t` and ends at a full period
beyond time `t`. If no SoI change occurs, the original orbit is returned.

For hyperbolic orbits, the search begins at time `t` and ends at the time of SoI ejection.
If the orbit's primary body has no SoI (i.e. it is a Star), the search ends at the time when
the orbit is gauranteed to have passed beyond the outermost satellite of the star.
"""
function propagate_across_soi(t, orb::Orbit)
    if orb.e < 1
        endtime = orb.period
    else
        if typeof(orb.primary) == Star
            f_sat = max(orb.primary.satellite_bodies)
            endtime = true_to_time(orbital_angle(f_sat.a*(1+f_sat.e)+f_sat.SoI, orb), orb)
        else
            endtime = ejection_time(orb)
        end
    end

    intersects = Tuple{CelestialBody, Float64}[]
    for bd in orb.primary.satellite_bodies
        #TO DO: ensure that earliest intersection is found, if multiple exist
        obj_fun(t) = soi_distance(t, orb, bd)
        res = optimize(obj_fun, t, endtime, Brent())
        if Optim.minimum(res) < 0
            push!(intersects, (bd, Optim.minimizer(res)))
            endtime = Optim.minimizer(res)
        end
    end

    if !isempty(intersects)
        intsct_bd, mindist_time = intersects[end]
        obj_fun2(t) = abs(soi_distance(t, orb, intsct_bd))
        res2 = optimize(obj_fun2, t, mindist_time, Brent())
        intsct_time = Optim.minimizer(res2)
        r̄ = time_orbital_position(intsct_time, orb) - time_orbital_position(intsct_time, intsct_bd.orbit)
        v̄ = time_orbital_velocity(intsct_time, orb) - time_orbital_velocity(intsct_time, intsct_bd.orbit)
        @show r̄, v̄, intsct_time
    elseif orb.e > 1
        intsct_bd = orb.primary.orbit.primary
        intsct_time = endtime
        r̄ = time_orbital_position(intsct_time, orb) + time_orbital_position(intsct_time, orb.primary.orbit)
        v̄ = time_orbital_velocity(intsct_time, orb) + time_orbital_velocity(intsct_time, orb.primary.orbit)
    else
        return orb
    end
    return Orbit(intsct_time, r̄, v̄, intsct_bd)
end


text:
function scatter(
    X::Vector{Float64}, Y::Vector{Float64},
    SZ::Vector{Float64}, COL::Vector{String},
    STYLE::Symbol,
    FILL::Bool,
    OPACITY::Real=1
    )
    f = nothing
    if FILL
        f = svg_points
    else
        if STYLE==:circle
            f = svg_circles
        elseif STYLE==:dot
            f = svg_points
        else
            f = svg_circles
        end
    end
    lsvg = [f(collect(zip(X[COL.==color],Y[COL.==color])), SZ[COL.==color], color, OPACITY) for color ∈ COL]
    return vcat(lsvg...)
end

function scatter(
    X::Vector{Float64}, Y::Vector{Float64},
    SZ::Float64, COL::Vector{String},
    STYLE::Symbol,
    FILL::Bool,
    OPACITY::Real=1
    )
    scatter(X, Y, [SZ for i=1:length(X)], COL, STYLE, FILL, OPACITY)
end

function scatter(
    X::Vector{Float64}, Y::Vector{Float64},
    SZ, COL::String,
    STYLE::Symbol,
    FILL::Bool,
    OPACITY::Real=1
    )
    scatter(X, Y, SZ, [COL for i=1:length(X)], STYLE, FILL, OPACITY)
end

max_stars_repo_name: algorithmx/PlotSVG

max_issues_repo_name: algorithmx/PlotSVG

max_forks_repo_name: algorithmx/PlotSVG