text:
#shape example.
import Pkg
Pkg.activate(".")
using WavesAndEigenvalues.Helmholtz

## configure set-up
mesh=Mesh("./examples/tutorials/Rijke_mm.msh",scale=0.001)
case="Rijke_test"
h=1E-9
mesh=octosplit(mesh);case*="_fine"  #activate this line to refine the mesh

dscrp=Dict()
dscrp["Interior"]=(:interior,())
dscrp["Outlet"]=(:admittance, (:Y,1E15))
hot=false
hot=true#activate this line for the active solution
if !hot
    c=ones(length(mesh.tetrahedra))*347.0
    case*="_cold"
    init=510.0
else
    dscrp["Outlet"]=(:admittance, (:Y,1E15)) #specify outlet BC
    γ=1.4 #ratio of specific heats
    ρ=1.225 #density at the reference location upstream to the flame in kg/m^3
    Tu=300.0    #K unburnt gas temperature
    Tb=1200.0    #K burnt gas temperature
    P0=101325.0 # ambient pressure in Pa
    A=pi*0.025^2 # cross sectional area of the tube
    Q02U0=P0*(Tb/Tu-1)*A*γ/(γ-1) #the ratio of mean heat release to mean velocity Q02U0
    x_ref=[0.0; 0.0; -0.00101] #reference point
    n_ref=[0.0; 0.0; 1.00] #directional unit vector of reference velocity
    n=1 #interaction index
    τ=0.001 #time delay
    ref_idx =  find_tetrahedron_containing_point(mesh,x_ref)
    dscrp["Flame"]=(:flame,(γ,ρ,Q02U0,ref_idx,x_ref,n_ref,:n,:τ,n,τ)) #flame dynamics
    #dscrp["Flame"]=(:fancyflame,(γ,ρ,Q02U0,x_ref,n_ref,[:n1,:n2],[:τ1,:τ2],[:a1,:a2],[1.0,.5],[.001,0.0005],[.01,.01]))
    #unify!(mesh,"Flame2","Flame")
    #D["Flame"]=(:fancyflame,(γ,ρ,Q02U0,x_ref,n_ref,:n1,:τ1,:a1,1.0,.001,.01,),)
    #D["Flame2"]=(:fancyflame,(γ,ρ,Q02U0,x_ref,n_ref,:n1,:τ1,:a1,.5,.0005,.01,),)
    R=287.05 # J/(kg*K) specific gas constant (air)
    speedofsound(x,y,z) = z<0. ? sqrt(γ*R*Tu) : sqrt(γ*R*Tb)
    c=generate_field(mesh,speedofsound)
    case*="_hot"
    init=700.0
end

##
L=discretize(mesh, dscrp, c)

##
sol, n, flag = householder(L,init*2*pi,maxiter=14, tol=1E-11,output = true, n_eig_val=3)

## prepare shape
# identify surface points
surface_points, tri_mask, tet_mask =get_surface_points(mesh)
# compute normal vectors on surface triangles
normal_vectors=get_normal_vectors(mesh)
## DA approach
sens=discrete_adjoint_shape_sensitivity(mesh,dscrp,c,surface_points,tri_mask,tet_mask,L,sol,h=h)
fd_sens=forward_finite_differences_shape_sensitivity(mesh,dscrp,c,surface_points,tri_mask,tet_mask,L,sol,h=h)

## Postprocessing
#normalize point sensitivity with directed surface triangle area
normed_sens=normalize_sensitivity(surface_points,normal_vectors,tri_mask,sens)
# boundary-mass-based-normalizations
nsens = bound_mass_normalize(surface_points,normal_vectors,tri_mask,mesh,sens)
# compute sensitivity in unit normal direction
normal_sens = normal_sensitivity(normal_vectors, normed_sens)

## save data
DD=Dict()
DD["real normed_sens1"]=real.(normed_sens[1,:])
DD["real normed_sens2"]=real.(normed_sens[2,:])
DD["real normed_sens3"]=real.(normed_sens[3,:])
DD["real normal_sens"]=real.(normal_sens)
DD["imag normed_sens1"]=imag.(normed_sens[1,:])
DD["imag normed_sens2"]=imag.(normed_sens[2,:])
DD["imag normed_sens3"]=imag.(normed_sens[3,:])
DD["imag normal_sens"]=imag.(normal_sens)
vtk_write_tri(case*"_tri",mesh,DD)

##
mode="[$(string(round((sol.params[:ω]/2/pi),digits=2)))]Hz"
data=Dict()
data["abs_p"*mode]=abs.(sol.v)./maximum(abs.(sol.v))
data["phase_p"*mode]=angle.(sol.v)
data["abs_p_adj"*mode]=abs.(sol.v_adj)./maximum(abs.(sol.v_adj))
data["phase_p_adj"*mode]=angle.(sol.v_adj)
data["real DA x"*mode]=real.(sens[1,:])
data["real DA y"*mode]=real.(sens[2,:])
data["real DA z"*mode]=real.(sens[3,:])
data["imag DA x"*mode]=imag.(sens[1,:])
data["imag DA y"*mode]=imag.(sens[2,:])
data["imag DA z"*mode]=imag.(sens[3,:])
data["real nDA x"*mode]=real.(nsens[1,:])
data["real nDA y"*mode]=real.(nsens[2,:])
data["real nDA z"*mode]=real.(nsens[3,:])
data["imag nDA x"*mode]=imag.(nsens[1,:])
data["imag nDA y"*mode]=imag.(nsens[2,:])
data["imag nDA z"*mode]=imag.(nsens[3,:])
#data["finite difference"*mode]=real.(sens_FD)
#data["central finite difference"]=real.(sens_CFD)
#data["diff DA -- FD"*mode]=abs.((sens.-sens_FD))
vtk_write(case, mesh, data)
## FD approach
#fd_sens=forward_finite_differences_shape_sensitivity(mesh,dscrp,c,surface_points,tri_mask,tet_mask,L,sol,h=h)
## write output
mode="[$(string(round((sol.params[:ω]/2/pi),digits=2)))]Hz"
data=Dict()
data["abs_p"*mode]=abs.(sol.v)./maximum(abs.(sol.v))
data["phase_p"*mode]=angle.(sol.v)
data["abs_p_adj"*mode]=abs.(sol.v_adj)./maximum(abs.(sol.v_adj))
data["phase_p_adj"*mode]=angle.(sol.v_adj)
data["real DA x"*mode]=real.(fd_sens[1,:])
data["real DA y"*mode]=real.(fd_sens[2,:])
data["real DA z"*mode]=real.(fd_sens[3,:])
data["imag DA x"*mode]=imag.(fd_sens[1,:])
data["imag DA y"*mode]=imag.(fd_sens[2,:])
data["imag DA z"*mode]=imag.(fd_sens[3,:])
#data["finite difference"*mode]=real.(fd_sens_FD)
#data["central finite difference"]=real.(sens_CFD)
#data["diff DA -- FD"*mode]=abs.((sens.-sens_FD))
vtk_write(case*"_fd", mesh, data)
println("Done writing!")


text:
struct JAXTAMError <: Exception
    msg   :: String
    step  :: Symbol
    super :: Union{Nothing, Exception}

    JAXTAMError(msg::String, step::Symbol) = new(msg, step, nothing)

    JAXTAMError(msg::String, step::Symbol, super::Exception) = new(msg, step, super)
end

ext: jl


text:
# Julia version: LTS (v1.0.3)
using Test

graph = Dict()
graph["start"] = Dict()
graph["start"]["a"] = 6
graph["start"]["b"] = 2

graph["a"] = Dict()
graph["a"]["fin"] = 1

graph["b"] = Dict()
graph["b"]["a"] = 3
graph["b"]["fin"] = 5

graph["fin"] = Dict()

# the costs table
infinity = Inf
costs = Dict()
costs["a"] = 6
costs["b"] = 2
costs["fin"] = infinity

# the parents table
parents = Dict()
parents["a"] = "start"
parents["b"] = "start"
parents["fin"] = nothing

processed = []

function find_lowest_cost_node(costs)
    lowest_cost = Inf
    lowest_cost_node = nothing
    for node in keys(costs)
        cost = costs[node]
        if cost < lowest_cost && node ∉ processed
            lowest_cost = cost
            lowest_cost_node = node
        end
    end
    return lowest_cost_node
end

node = find_lowest_cost_node(costs)
while node != nothing
    global node
    cost = costs[node]
    neighbors = graph[node]
    for n in keys(neighbors)
        new_cost = cost + neighbors[n]
        if costs[n] > new_cost
            costs[n] = new_cost
            parents[n] = node
        end
    end
    push!(processed,node)
    node = find_lowest_cost_node(costs)
end

println("Cost from the start to each node:")
println(costs)
@test costs == Dict("fin"=>6,"b"=>2,"a"=>5)


text:
function Player1(x::Int,y::Int, num_food_types::Int, starting_food::Int)
    # Yes these types are necessary
    Player{num_food_types}(
        (x,y),
        NTuple{num_food_types, Int}(starting_food for i in 1:num_food_types),
        false,
        starting_food)
end


text:
"""
    pkg ReDiverse
---

# General description
A package to organize and analysis ReDiverse data. Below are the `ToDo`s

# Urgent
- `Start imputation asap`

# General roadmap
1. Genotype collection, and convert to `plink` format.
2. [x] Genotype data clean
3. [x] Determine of SNP chip platforms, and final SNP set.
4. [x] Imputation → multi-breed GRM
   - [x] GRM validation
5. GBLUP
6. Cross-validation strategy

"""
module ReDiverse
using Test
# , Printf, DataFrames, CSV
# using Statistics, LinearAlgebra
################################################################################

export sandbox_ReDiverse

work_dir = pwd()                # all data and results are relative to work_dir
rdDir = begin                   # package directory
    t = splitdir(pathof(ReDiverse))[1]
    l = findlast('/', t) - 1
    t[1:l]
end
rdBin = joinpath(rdDir, "bin")
rdCpp = joinpath(rdDir, "src/cpp")

beagle = joinpath(rdBin, "beagle.jar")
plink = joinpath(rdBin, "plink")

#include("work-flow.jl")
#
## raw genotype orgnization
#include("genotype/prepare-maps.jl")
#include("genotype/raw-genotype-org.jl")
#include("genotype/norsk-genotypes.jl")
#include("genotype/quality-control.jl")
#include("genotype/imputation.jl")
#
## about the genomic relationship matrix
#include("grm/calc-grm.jl")
#include("grm/check-grm.jl")
#
## some utilities
#include("utils/makefile.jl")
#include("utils/commonFunc.jl")

#include("utils/abg-cmds.jl")
end # module


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

#using ExoplanetsSysSim
#using DataFrames
#using JLD

type KeplerPhysicalCatalog
  target::Array{KeplerTarget,1}
end
#KeplerPhysicalCatalog() = KeplerPhysicalCatalog([])

type KeplerObsCatalog
  target::Array{KeplerTargetObs,1}
end
#KeplerObsCatalog() = KeplerObsCatalog([])

function generate_kepler_physical_catalog(sim_param::SimParam)
   if haskey(sim_param,"stellar_catalog")
      star_tab_func = get_function(sim_param, "star_table_setup")
      star_tab_func(sim_param)
   end
   num_sys = get_int(sim_param,"num_targets_sim_pass_one")
   generate_kepler_target = get_function(sim_param,"generate_kepler_target")
   target_list = Array{KeplerTarget}(num_sys)
   map!(generate_kepler_target, target_list, fill(sim_param,num_sys) )
  return KeplerPhysicalCatalog(target_list)
end

function observe_kepler_targets_sky_avg(input::KeplerPhysicalCatalog, sim_param::SimParam )
  calc_target_obs = get_function(sim_param,"calc_target_obs_sky_ave")
  return observe_kepler_targets(calc_target_obs, input, sim_param)
end

function observe_kepler_targets_single_obs(input::KeplerPhysicalCatalog, sim_param::SimParam )
  calc_target_obs = get_function(sim_param,"calc_target_obs_single_obs")
  return observe_kepler_targets(calc_target_obs, input, sim_param)
end

function observe_kepler_targets(calc_target_obs::Function, input::KeplerPhysicalCatalog, sim_param::SimParam )
  #calc_target_obs = get_function(sim_param,"calc_target_obs_sky_ave")
  #calc_target_obs = get_function(sim_param,"calc_target_obs_single_obs")
  output = KeplerObsCatalog([])
  if haskey(sim_param,"mem_kepler_target_obs")
     output.target = get(sim_param,"mem_kepler_target_obs",Array{KeplerTargetObs}(0) )
  end
  num_targets_sim_pass_one = get_int(sim_param,"num_targets_sim_pass_one")
  if length(output.target) < num_targets_sim_pass_one
     output.target = Array{KeplerTargetObs}(num_targets_sim_pass_one)
  end
  #output.target = Array{KeplerTargetObs}(length(input.target) )  # Replaced to reduce memory allocation
  map!(x::KeplerTarget->calc_target_obs(x,sim_param)::KeplerTargetObs, output.target, input.target)
  return output
end

function select_targets_one_obs(ps::PlanetarySystemAbstract)
 for pl in 1:length(ps.orbit)
   ecc::Float64 = ps.orbit[pl].ecc
   incl::Float64 = ps.orbit[pl].incl
   a::Float64 = semimajor_axis(ps,pl)
   Rstar::Float64 = rsol_in_au*ps.star.radius
   if (Rstar > (a*(1-ecc)*(1+ecc))/(1+ecc*sin(ps.orbit[pl].omega))*cos(incl))
     return true
   end
 end
 return false
end

function generate_obs_targets(cat_phys::KeplerPhysicalCatalog, sim_param::SimParam )
  for t in 1:length(cat_phys.target)
    for ps in 1:length(cat_phys.target[t].sys)
      kep_targ = cat_phys.target[t].sys[ps]
      for pl in length(kep_targ.orbit):-1:1
        ecc::Float64 = kep_targ.orbit[pl].ecc
	incl::Float64 = kep_targ.orbit[pl].incl
   	a::Float64 = semimajor_axis(kep_targ,pl)
   	Rstar::Float64 = rsol_in_au*kep_targ.star.radius

   	if (Rstar < (a*(1-ecc)*(1+ecc))/(1+ecc*sin(kep_targ.orbit[pl].omega))*cos(incl)) || (rand() > calc_prob_detect_if_transit(cat_phys.target[t], 1, pl, sim_param))
    	  splice!(cat_phys.target[t].sys[ps].orbit, pl)
	  splice!(cat_phys.target[t].sys[ps].planet, pl)
     	end
      end
    end
  end
  return cat_phys
end

function simulated_read_kepler_observations(sim_param::SimParam ) # TODO SCI:  IMPORTANT:  Eventually, replace this with a function to read data from input file (see koi_table.jl)
   if haskey(sim_param,"stellar_catalog")
      star_tab_func = get_function(sim_param, "star_table_setup")
      star_tab_func(sim_param)
   end
   num_sys = get_int(sim_param,"num_kepler_targets")
   generate_kepler_target = get_function(sim_param,"generate_kepler_target")
   target_list = Array{KeplerTarget}(num_sys)
   map!(x->generate_kepler_target(sim_param), target_list, 1:num_sys )

   cat_phys_cut = generate_obs_targets(KeplerPhysicalCatalog(target_list), sim_param)
   calc_target_obs = get_function(sim_param,"calc_target_obs_single_obs")
   output = KeplerObsCatalog([])
   output.target = map(x::KeplerTarget->calc_target_obs(x,sim_param)::KeplerTargetObs, cat_phys_cut.target)
   return output
end

# df_star is assumed to have fields kepid, mass and radius for all targets in the survey
function setup_actual_planet_candidate_catalog(df_star::DataFrame, sim_param::SimParam)
  add_param_fixed(sim_param,"num_kepler_targets",num_usable_in_star_table())  # For "observed" catalog
  koi_catalog_file_in = convert(ASCIIString,joinpath(Pkg.dir("ExoplanetsSysSim"), "data", convert(ASCIIString,get(sim_param,"koi_catalog","q1_q17_dr25_koi.csv")) ) )
  (csv_data,csv_header) =  readcsv(koi_catalog_file_in,header=true)
  # Lookup header columns, since DataFrames doesn't like this file
  kepid_idx = findfirst(x->x=="kepid",csv_header)
  koi_period_idx = findfirst(x->x=="koi_period",csv_header)
  koi_time0bk_idx = findfirst(x->x=="koi_time0bk",csv_header)
  koi_depth_idx = findfirst(x->x=="koi_depth",csv_header)
  koi_duration_idx = findfirst(x->x=="koi_duration",csv_header)
  koi_ror_idx = findfirst(x->x=="koi_ror",csv_header)
  koi_period_err1_idx = findfirst(x->x=="koi_period_err1",csv_header)
  koi_time0bk_err1_idx = findfirst(x->x=="koi_time0bk_err1",csv_header)
  koi_depth_err1_idx = findfirst(x->x=="koi_depth_err1",csv_header)
  koi_duration_err1_idx = findfirst(x->x=="koi_duration_err1",csv_header)
  koi_period_err2_idx = findfirst(x->x=="koi_period_err2",csv_header)
  koi_time0bk_err2_idx = findfirst(x->x=="koi_time0bk_err2",csv_header)
  koi_depth_err2_idx = findfirst(x->x=="koi_depth_err2",csv_header)
  koi_duration_err2_idx = findfirst(x->x=="koi_duration_err2",csv_header)
  koi_disposition_idx = findfirst(x->x=="koi_disposition",csv_header)
  koi_pdisposition_idx = findfirst(x->x=="koi_pdisposition",csv_header)
  # Choose which KOIs to keep
  #is_cand = (csv_data[:,koi_disposition_idx] .== "CONFIRMED") | (csv_data[:,koi_disposition_idx] .== "CANDIDATE")
  is_cand = (csv_data[:,koi_pdisposition_idx] .== "CANDIDATE")

  idx_keep = is_cand & !isna.(csv_data[:,koi_ror_idx]) & ([typeof(x) for x in csv_data[:,koi_ror_idx]] .== Float64)
  idx_keep = idx_keep & !isna.(csv_data[:,koi_period_err1_idx]) & ([typeof(x) for x in csv_data[:,koi_period_err1_idx]] .== Float64) # DR25 catalog missing uncertainties for some candidates
  csv_data = csv_data[idx_keep,:]

  output = KeplerObsCatalog([])

  for (j,kepid) in enumerate(df_star[:kepid])
     plids = find(x->x==kepid,csv_data[:,kepid_idx])
	 num_pl = length(plids)
	 if num_pl ==0 continue end
	 target_obs = KeplerTargetObs(num_pl)
	 target_obs.star = ExoplanetsSysSim.StarObs(df_star[j,:radius],df_star[j,:mass])
	 for (plid,i) in enumerate(plids)
	   target_obs.obs[plid] = ExoplanetsSysSim.TransitPlanetObs(csv_data[i,koi_period_idx],csv_data[i,koi_time0bk_idx],csv_data[i,koi_depth_idx]/1.0e6,csv_data[i,koi_duration_idx])
           target_obs.sigma[plid] = ExoplanetsSysSim.TransitPlanetObs((abs(csv_data[i,koi_period_err1_idx])+abs(csv_data[i,koi_period_err2_idx]))/2,(abs(csv_data[i,koi_time0bk_err1_idx])+abs(csv_data[i,koi_time0bk_err2_idx]))/2,(abs(csv_data[i,koi_depth_err1_idx]/1.0e6)+abs(csv_data[i,koi_depth_err2_idx]/1.0e6))/2,(abs(csv_data[i,koi_duration_err1_idx])+abs(csv_data[i,koi_duration_err2_idx]))/2)
	   target_obs.prob_detect = ExoplanetsSysSim.OneObserverSystemDetectionProbs(num_pl)
	 end
      push!(output.target,target_obs)
  end
  return output
end

function setup_actual_planet_candidate_catalog_csv(df_star::DataFrame, sim_param::SimParam)
  add_param_fixed(sim_param,"num_kepler_targets",num_usable_in_star_table())  # For "observed" catalog
  koi_catalog_file_in = convert(ASCIIString,joinpath(Pkg.dir("ExoplanetsSysSim"), "data", convert(ASCIIString,get(sim_param,"koi_catalog","q1_q17_dr25_koi.csv")) ) )
  (csv_data,csv_header) =  readcsv(koi_catalog_file_in,header=true)
  # Lookup header columns, since DataFrames doesn't like this file
  kepid_idx = findfirst(x->x=="kepid",csv_header)
  koi_period_idx = findfirst(x->x=="koi_period",csv_header)
  koi_time0bk_idx = findfirst(x->x=="koi_time0bk",csv_header)
  koi_depth_idx = findfirst(x->x=="koi_depth",csv_header)
  koi_duration_idx = findfirst(x->x=="koi_duration",csv_header)
  koi_ror_idx = findfirst(x->x=="koi_ror",csv_header)
  koi_period_err1_idx = findfirst(x->x=="koi_period_err1",csv_header)
  koi_time0bk_err1_idx = findfirst(x->x=="koi_time0bk_err1",csv_header)
  koi_depth_err1_idx = findfirst(x->x=="koi_depth_err1",csv_header)
  koi_duration_err1_idx = findfirst(x->x=="koi_duration_err1",csv_header)
  koi_period_err2_idx = findfirst(x->x=="koi_period_err2",csv_header)
  koi_time0bk_err2_idx = findfirst(x->x=="koi_time0bk_err2",csv_header)
  koi_depth_err2_idx = findfirst(x->x=="koi_depth_err2",csv_header)
  koi_duration_err2_idx = findfirst(x->x=="koi_duration_err2",csv_header)
  koi_disposition_idx = findfirst(x->x=="koi_disposition",csv_header)
  koi_pdisposition_idx = findfirst(x->x=="koi_pdisposition",csv_header)

  koi_subset = fill(false, length(csv_data[:,kepid_idx]))

  if haskey(sim_param, "koi_subset_csv")
    subset_df = readtable(convert(ASCIIString,get(sim_param,"koi_subset_csv", "christiansen_kov.csv")), header=true, separator=' ')

    col_idx = findfirst(x->x==string(names(subset_df)[1]),csv_header)
    for n in 1:length(subset_df[:,1])
      subset_colnum = 1
      subset_entry = find(x->x==subset_df[n,1], csv_data[:,col_idx])
      # println("Initial cut: ", subset_entry)
      while (length(subset_entry) > 1) & (subset_colnum < length(names(subset_df)))
        subset_colnum += 1
        col_idx = findfirst(x->x==string(names(subset_df)[subset_colnum]),csv_header)
        subsubset = find(x->round(x*10.)==round(subset_df[n,subset_colnum]*10.), csv_data[subset_entry,col_idx])
	# println("Extra cut: ", subset_df[n,subset_colnum], " / ", csv_data[subset_entry,col_idx], " = ", subsubset)
	subset_entry = subset_entry[subsubset]
      end
      if subset_colnum > 1
 	col_idx = findfirst(x->x==string(names(subset_df)[1]),csv_header)
      end
      if length(subset_entry) > 1
	cand_sub = find(x->x == "CANDIDATE",csv_data[subset_entry,koi_pdisposition_idx])
	subset_entry = subset_entry[cand_sub]
	if length(subset_entry) > 1
          println("Multiple planets found in final cut: ", subset_df[n,1])
        end
      end
      if length(subset_entry) < 1
        println("No planets found in final cut: ", subset_df[n,:])
      end
      koi_subset[subset_entry] = true
    end
  else
    println("No KOI csv file given!")
    koi_subset = fill(true, length(csv_data[:,kepid_idx]))
  end

  csv_data = csv_data[koi_subset,:]

  output = KeplerObsCatalog([])

  tot_plan = count(x->x, koi_subset)
  for (j,kepid) in enumerate(df_star[:kepid])
     plids = find(x->x==kepid,csv_data[:,kepid_idx])
     num_pl = length(plids)
     tot_plan = tot_plan - num_pl
     if num_pl ==0 continue end
     target_obs = KeplerTargetObs(num_pl)
     target_obs.star = ExoplanetsSysSim.StarObs(df_star[j,:radius],df_star[j,:mass])
     for (plid,i) in enumerate(plids)
       target_obs.obs[plid] = ExoplanetsSysSim.TransitPlanetObs(csv_data[i,koi_period_idx],csv_data[i,koi_time0bk_idx],csv_data[i,koi_depth_idx]/1.0e6,csv_data[i,koi_duration_idx])
       target_obs.sigma[plid] = ExoplanetsSysSim.TransitPlanetObs((abs(csv_data[i,koi_period_err1_idx])+abs(csv_data[i,koi_period_err2_idx]))/2,(abs(csv_data[i,koi_time0bk_err1_idx])+abs(csv_data[i,koi_time0bk_err2_idx]))/2,(abs(csv_data[i,koi_depth_err1_idx]/1.0e6)+abs(csv_data[i,koi_depth_err2_idx]/1.0e6))/2,(abs(csv_data[i,koi_duration_err1_idx])+abs(csv_data[i,koi_duration_err2_idx]))/2)
       target_obs.prob_detect = ExoplanetsSysSim.OneObserverSystemDetectionProbs(num_pl)
# ExoplanetsSysSim.OneObserverSystemDetectionProbs( ones(num_pl), ones(num_pl,num_pl), ones(num_pl), fill(Array{Int64}(0), 1) )
     end
     push!(output.target,target_obs)
  end
  println("Number of planet candidates with no matching star in table: ", tot_plan)
  return output
end

# Two functions below were just for debugging purposes
function calc_snr_list(cat::KeplerPhysicalCatalog, sim_param::SimParam)
  snrlist = Array{Float64}(0)
  for t in 1:length(cat.target)
    for p in 1:length(cat.target[t].sys[1].planet)
      snr = calc_snr_if_transit(cat.target[t],1,p,sim_param)
      if snr>0.0
        push!(snrlist,snr)
      end
    end
  end
  snrlist[find(x->x>7.1,snrlist)]
end

function calc_prob_detect_list(cat::KeplerPhysicalCatalog, sim_param::SimParam)
  pdetectlist = Array{Float64}(0)
  for t in 1:length(cat.target)
    for p in 1:length(cat.target[t].sys[1].planet)
      pdet = calc_prob_detect_if_transit(cat.target[t],1,p,sim_param)
      if pdet>0.0
        push!(pdetectlist,pdet)
      end
    end
  end
  idx = find(x->x>0.0,pdetectlist)
  pdetectlist[idx]
end

function test_catalog_constructors(sim_param::SimParam)
  cat_phys = generate_kepler_physical_catalog(sim_param)::KeplerPhysicalCatalog
  id = findfirst( x->num_planets(x)>=1 , cat_phys.target)   # fast forward to first target that has some planets
  @assert(length(id)>=1)
  semimajor_axis(cat_phys.target[id].sys[1],1)
  pdetlist = calc_prob_detect_list(cat_phys,sim_param)
  calc_target_obs_single_obs(cat_phys.target[id],sim_param)
  calc_target_obs_sky_ave(cat_phys.target[id],sim_param)
  @assert( length(cat_phys.target[id].sys[1].planet)  == num_planets(cat_phys.target[id]) )
  cat_obs = simulated_read_kepler_observations(sim_param)
  return (cat_phys, cat_obs)
end


text:
include("./preprocess.jl")
using BenchmarkTools, Statistics
using Random
## Modify for different problem size
inpGrid = 4
numOfPoints = 1<<17
numOfQueries = 1<<17
Points = 0
Queries = 0
dimensions = 3
if length(ARGS) == 1
    println("Running for .csv files")
    global inpGrid
    global numOfPoints
    global numOfQueries
    global Points
    global Queries
    global dimensions
    inpGrid = parse(Int64, ARGS[1])
    Points = ReadCSV("./pts.csv")
    Queries = ReadCSV("./qrs.csv")
    numOfPoints = size(Points)[1]
    numOfQueries = size(Queries)[1]
    dimensions = size(Points)[2]
elseif length(ARGS) == 3
    println("Running for random generated dataset")
    global rseed
    global numOfPoints
    global numOfQueries
    global inpGrid
    global Points
    global Queries
    global dimensions
    Random.seed!(parse(Int64, ARGS[3]))
    numOfPoints = 1<< parse(Int64, ARGS[1])
    numOfQueries = 1<<parse(Int64, ARGS[1])
    inpGrid = parse(Int64, ARGS[2])
    Points = rand(Float32, numOfPoints ,dimensions)
    Queries = rand(Float32, numOfQueries, dimensions)

end

numOfGrids = 1<<inpGrid

BlockOfPoint = AssignPointsToBlock(Points, numOfGrids, dimensions)
BlockOfQuery = AssignPointsToBlock(Queries, numOfGrids, dimensions)

PointsPerBlock, IntegralPointsPerBlock = CountPointsPerBlock(Points, numOfGrids, dimensions)
QueriesPerBlock, IntegralQueriesPerBlock = CountPointsPerBlock(Queries, numOfGrids, dimensions)

OrderedPoints = ReorderPointsByBlock(Points, BlockOfPoint)
OrderedQueries = ReorderPointsByBlock(Queries, BlockOfQuery)
println("Num pts: ", numOfPoints, " Gridsize: ", numOfGrids)

############################################################
#   Uncomment the implementations you want to benchmark.   #
############################################################
include("kernel_simple.jl")

println("Running for simple version: ")
gpu_idxs, gpu_dists = cuda_knn_simple(OrderedPoints, OrderedQueries,PointsPerBlock,
    QueriesPerBlock, IntegralPointsPerBlock, IntegralQueriesPerBlock,numOfPoints,
    numOfQueries, numOfGrids, dimensions)
#
println(gpu_idxs[1:10])
println(gpu_dists[1:10])
##############################################
include("kernel_simple_view_function.jl")

println("Running for simple version with @view and distance calcs in function: ")
gpu_idxs, gpu_dists = cuda_knn_simple_view_function(OrderedPoints, OrderedQueries,PointsPerBlock,
    QueriesPerBlock, IntegralPointsPerBlock, IntegralQueriesPerBlock,numOfPoints,
    numOfQueries, numOfGrids, dimensions)

println(gpu_idxs[1:10])
println(gpu_dists[1:10])
###############################################
include("kernel_with_skip.jl")

println("Running for version with skip: ")
gpu_idxs, gpu_dists = cuda_knn_with_skip(OrderedPoints, OrderedQueries,PointsPerBlock,
    QueriesPerBlock, IntegralPointsPerBlock, IntegralQueriesPerBlock,numOfPoints,
    numOfQueries, numOfGrids, dimensions)

println(gpu_idxs[1:10])
println(gpu_dists[1:10])
###############################################
include("kernel_with_skip_view_function.jl")
println("Running for version with skip and @view and distance calcs in function: ")
gpu_idxs, gpu_dists = cuda_knn_with_skip_view_function(OrderedPoints, OrderedQueries,PointsPerBlock,
    QueriesPerBlock, IntegralPointsPerBlock, IntegralQueriesPerBlock,numOfPoints,
    numOfQueries, numOfGrids, dimensions)

println(gpu_idxs[1:10])
println(gpu_dists[1:10])
###############################################
include("kernel_simple_no_types.jl")
println("Running for version with no types, skip and @view and distance calcs in function: ")
gpu_idxs, gpu_dists = cuda_knn_simple_no_types_function(OrderedPoints, OrderedQueries,PointsPerBlock,
    QueriesPerBlock, IntegralPointsPerBlock, IntegralQueriesPerBlock,numOfPoints,
    numOfQueries, numOfGrids, dimensions)

println(gpu_idxs[1:10])
println(gpu_dists[1:10])


text:
# This file holds the functions to save the configuration of a manpower
#   simulation to an SQLite database.

export saveSimConfigToDatabase

"""
```
saveSimConfigToDatabase( mpSim::ManpowerSimulation,
                         dbName::String,
                         configName::String = "config" )
```
This function stores the configuration parameters of the manpower simulation
`mpSim` in the database with filename `dbName`, in the table with name
`configName`. If the database filename doesn't have the proper extension,
`.sqlite`, it will be appended to the name.

This function returns `nothing`.
"""
function saveSimConfigToDatabase( mpSim::ManpowerSimulation, dbName::String,
    configName::String = "config" )::Void

    tmpDBname = ( dbName == "" ) || endswith( dbName, ".sqlite" ) ?
        dbName : dbName * ".sqlite"
    configDB = SQLite.DB( tmpDBname )

    createConfigTable( configName, configDB )

    command = "BEGIN TRANSACTION"
    SQLite.execute!( configDB, command )

    readGeneralParsFromSim( mpSim, configDB, configName )
    readAttributesFromSim( mpSim, configDB, configName )
    readStatesFromSim( mpSim, configDB, configName )
    readTransitionsFromSim( mpSim, configDB, configName )
    readRecruitmentFromSim( mpSim, configDB, configName )
    readAttritionFromSim( mpSim, configDB, configName )
    readRetirementFromSim( mpSim, configDB, configName )

    command = "COMMIT"
    SQLite.execute!( configDB, command )

    return

end  # saveSimConfigToDatabase( mpSim, dbName, configName )

"""
```
saveSimConfigToDatabase( mpSim::ManpowerSimulation )
```
This function stores the configuration parameters of the manpower simulation
`mpSim` in table `config` of the database defined in the simulation.

This function returns `nothing`.
"""
function saveSimConfigToDatabase( mpSim::ManpowerSimulation )::Void

    saveSimConfigToDatabase( mpSim, mpSim.dbName )

end  # saveSimConfigToDatabase( mpSim::ManpowerSimulation )

# ==============================================================================
# Non-exported methods.
# ==============================================================================

"""
```
readGeneralParsFromSim( mpSim::ManpowerSimulation,
                        configDB::SQLite.DB,
                        configName::String )
```
This function reads the general simulation parameters from the manpower
simulation `mpSim` and writes them to the SQLite database `configDB` in the
table with name `configName`.

This function returns `nothing`.
"""
function readGeneralParsFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    # Name of the simulation database file.
    command = "INSERT INTO config
        (parName, parType, strPar) VALUES
        ('dbName', 'General', '$(mpSim.dbName)')"
    SQLite.execute!( configDB, command )

    # Name of the simulation.
    simName = replace( mpSim.personnelDBname, "Personnel_", "", 1 )
    command = "INSERT INTO config
        (parName, parType, strPar) VALUES
        ('simName', 'General', '$simName')"
    SQLite.execute!( configDB, command )

    # Personnel cap.
    command = "INSERT INTO config
        (parName, parType, intPar) VALUES
        ('dbName', 'General', '$(mpSim.personnelTarget)')"
    SQLite.execute!( configDB, command )

    # Simulation length.
    simLength = mpSim.simLength / 12
    command = "INSERT INTO config
        (parName, parType, realPar) VALUES
        ('dbName', 'General', '$simLength')"
    SQLite.execute!( configDB, command )

    # Database commits.
    dbCommits = Int( mpSim.simLength / mpSim.commitFrequency )
    command = "INSERT INTO config
        (parName, parType, intPar) VALUES
        ('dbName', 'General', '$dbCommits')"
    SQLite.execute!( configDB, command )

    return

end  # readGeneralParsFromSim( mpSim, configDB, configName )

"""
```
readAttributesFromSim( mpSim::ManpowerSimulation,
                       configDB::SQLite.DB,
                       configName::String )
```
This function reads the personnel attributes from the manpower simulation
`mpSim` and writes them to the SQLite database `configDB` in the table with name
`configName`.

This function returns `nothing`.
"""
function readAttributesFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    # Read every single attribute.
    foreach( attr -> saveAttrToDatabase( attr, configDB, configName ),
        vcat( mpSim.initAttrList, mpSim.otherAttrList ) )

    return

end  # readAttributesFromSim( mpSim, configDB, configName )

"""
```
readStatesFromSim( mpSim::ManpowerSimulation,
                   configDB::SQLite.DB,
                   configName::String )
```
This function reads the personnel states from the manpower simulation `mpSim`
and writes them to the SQLite database `configDB` in the table with name
`configName`.

This function returns `nothing`.
"""
function readStatesFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    # Read every single state.
    foreach( state -> saveStateToDatabase( state, state.attrScheme.name,
        configDB, configName ), keys( mpSim.initStateList ) )
    foreach( state -> saveStateToDatabase( state, state.attrScheme.name,
        configDB, configName ), keys( mpSim.otherStateList ) )

    return

end  # readStatesFromSim( mpSim, configDB, configName )

"""
```
readTransitionsFromSim( mpSim::ManpowerSimulation,
                        configDB::SQLite.DB,
                        configName::String )
```
This function reads the state transitions from the manpower simulation `mpSim`
and writes them to the SQLite database `configDB` in the table with name
`configName`.

This function returns `nothing`.
"""
function readTransitionsFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    # Read every single transition.
    for state in keys( mpSim.initStateList )
        foreach( trans -> saveTransitionToDatabase( trans, state.name,
            trans.endState.name, configDB, configName ),
            mpSim.initStateList[ state ] )
    end  # for state in keys( mpSim.initStateList )

    for state in keys( mpSim.otherStateList )
        foreach( trans -> saveTransitionToDatabase( trans, state.name,
            trans.endState.name, configDB, configName ),
            mpSim.otherStateList[ state ] )
    end  # for state in keys( mpSim.otherStateList )

    return

end  # readStatesFromSim( mpSim, configDB, configName )

"""
```
readRecruitmentFromSim( mpSim::ManpowerSimulation,
                        configDB::SQLite.DB,
                        configName::String )
```
This function reads the recruitment schemes from the manpower simulation
`mpSim` and writes them to the SQLite database `configDB` in the table with name
`configName`.

This function returns `nothing`.
"""
function readRecruitmentFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    # Read every single recruitment scheme.
    foreach( recScheme -> saveRecruitmentToDatabase( recScheme, configDB,
        configName ), mpSim.recruitmentSchemes )

    return

end  # readRecruitmentFromSim( mpSim, configDB, configName )

"""
```
readAttritionFromSim( mpSim::ManpowerSimulation,
                      configDB::SQLite.DB,
                      configName::String )
```
This function reads the attrition parameters from the manpower simulation
`mpSim` and writes them to the SQLite database `configDB` in the table with name
`configName`.

This function returns `nothing`.
"""
function readAttritionFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    attrScheme = mpSim.defaultAttritionScheme
    saveAttritionToDatabase( attrScheme, configDB, configName )
    foreach( attrName -> saveAttritionToDatabase(
        mpSim.attritionSchemes[ attrName ], configDB, configName),
        keys( mpSim.attritionSchemes ) )
    return

end  # readAttritionFromSim( mpSim, configDB, configName )

"""
```
readRetirementFromSim( mpSim::ManpowerSimulation,
                       configDB::SQLite.DB,
                       configName::String )
```
This function reads the retirement parameters from the manpower simulation
`mpSim` and writes them to the SQLite database `configDB` in the table with name
`configName`.

This function returns `nothing`.
"""
function readRetirementFromSim( mpSim::ManpowerSimulation,
    configDB::SQLite.DB, configName::String )::Void

    # Only process if there is a retirement scheme.
    if mpSim.retirementScheme === nothing
        return
    end  # if mpSim.retirementScheme === nothing

    retScheme = mpSim.retirementScheme

    command = "INSERT INTO $configName
        (parName, parType, boolPar, strPar) VALUES
        ('Retirement', 'Retirement', '$(retScheme.isEither)',
            '$(retScheme.maxCareerLength),$(retScheme.retireAge),$(retScheme.retireFreq),$(retScheme.retireOffset)')"
    SQLite.execute!( configDB, command )

    return

end  # readRetirementFromSim( mpSim, configDB, configName )


text:
# functions related to gamma distribution

import .RFunctions:
    gammapdf,
    gammalogpdf,
    gammacdf,
    gammaccdf,
    gammalogcdf,
    gammalogccdf,
    gammainvcdf,
    gammainvccdf,
    gammainvlogcdf,
    gammainvlogccdf


text:
"""
For an initial guess from the index of a Fold bifurcation point located in ContResult.specialpoint, returns a point which will be refined using `newtonFold`.
"""
function FoldPoint(br::AbstractBranchResult, index::Int)
	bptype = br.specialpoint[index].type
	@assert bptype == :bp || bptype == :nd || bptype == :fold "This should be a Fold / BP point"
	specialpoint = br.specialpoint[index]
	return BorderedArray(_copy(specialpoint.x), specialpoint.param)
end
####################################################################################################
# constructors
FoldProblemMinimallyAugmented(F, J, Ja, d2F, lens::Lens, a, b, issymmetric::Bool, linsolve::AbstractLinearSolver, linbdsolver = MatrixBLS()) = FoldProblemMinimallyAugmented(F, J, Ja, d2F, lens, a, b, 0*a, linsolve, linsolve, linbdsolver, linbdsolver, issymmetric)

FoldProblemMinimallyAugmented(F, J, Ja, d2F, lens::Lens, a, b, linsolve::AbstractLinearSolver, linbdsolver = MatrixBLS()) = FoldProblemMinimallyAugmented(F, J, Ja, d2F, lens, a, b, false, linsolve, linbdsolver)

FoldProblemMinimallyAugmented(F, J, Ja, lens::Lens, a, b, issymmetric::Bool, linsolve::AbstractLinearSolver, linbdsolver = MatrixBLS()) = FoldProblemMinimallyAugmented(F, J, Ja, nothing, lens, a, b, 0*a, linsolve, linbdsolver, issymmetric)

FoldProblemMinimallyAugmented(F, J, Ja, lens::Lens, a, b, linsolve::AbstractLinearSolver, linbdsolver = MatrixBLS()) = FoldProblemMinimallyAugmented(F, J, Ja, lens, a, b, false, linsolve, linbdsolver)
####################################################################################################
getVec(x, ::FoldProblemMinimallyAugmented) = getVec(x)
getP(x, ::FoldProblemMinimallyAugmented) = getP(x)

function (fp::FoldProblemMinimallyAugmented)(x, p::T, params) where T
	# These are the equations of the minimally augmented (MA) formulation of the Fold bifurcation point
	# input:
	# - x guess for the point at which the jacobian is singular
	# - p guess for the parameter value `<: Real` at which the jacobian is singular
	# The jacobian of the MA problem is solved with a BLS method
	a = fp.a
	b = fp.b
	# update parameter
	par = set(params, fp.lens, p)
	# ┌      ┐┌  ┐ ┌ ┐
	# │ J  a ││v │=│0│
	# │ b  0 ││σ1│ │1│
	# └      ┘└  ┘ └ ┘
	# In the notations of Govaerts 2000, a = w, b = v
	# Thus, b should be a null vector of J
	#       a should be a null vector of J'
	# we solve Jv + a σ1 = 0 with <b, v> = n
	# the solution is v = -σ1 J\a with σ1 = -n/<b, J^{-1}a>
	n = T(1)
	J = fp.J(x, par)
	σ1 = fp.linbdsolver(J, a, b, T(0), fp.zero, n)[2]
	return fp.F(x, par), σ1
end

# this function encodes the functional
function (foldpb::FoldProblemMinimallyAugmented)(x::BorderedArray, params)
	res = foldpb(x.u, x.p, params)
	return BorderedArray(res[1], res[2])
end

# Struct to invert the jacobian of the fold MA problem.
struct FoldLinearSolverMinAug <: AbstractLinearSolver; end

function foldMALinearSolver(x, p::T, pb::FoldProblemMinimallyAugmented, par,
							rhsu, rhsp;
							debugArray = nothing) where T
	################################################################################################
	# debugArray is used as a temp to be filled with values used for debugging. If debugArray = nothing, then no debugging mode is entered. If it is AbstractArray, then it is used
	################################################################################################
	# Recall that the functional we want to solve is [F(x,p), σ(x,p)] where σ(x,p) is computed in the above function.
	# The jacobian has to be passed as a tuple as Jac_fold_MA(u0, pb::FoldProblemMinimallyAugmented) = (return (u0, pb, d2F::Bool))
	# The Jacobian J of the vector field is expressed at (x, p)
	# We solve here Jfold⋅res = rhs := [rhsu, rhsp]
	# The Jacobian expression of the Fold problem is
	#           ┌         ┐
	#  Jfold =  │  J  dpF │
	#           │ σx   σp │
	#           └         ┘
	# where σx := ∂_xσ and σp := ∂_pσ
	# We recall the expression of σx = -< w, d2F(x,p)[v, x2]> where (w, σ2) is solution of J'w + b σ2 = 0 with <a, w> = n
	########################## Extraction of function names ########################################
	F = pb.F
	J = pb.J

	d2F = pb.d2F
	a = pb.a
	b = pb.b

	# parameter axis
	lens = pb.lens
	# update parameter
	par0 = set(par, lens, p)

	# we define the following jacobian. It is used at least 3 times below. This avoids doing 3 times the (possibly) costly building of J(x, p)
	J_at_xp = J(x, par0)

	# we do the following in order to avoid computing J_at_xp twice in case pb.Jadjoint is not provided
	if issymmetric(pb)
		JAd_at_xp = J_at_xp
	else
		JAd_at_xp = hasAdjoint(pb) ? pb.Jᵗ(x, par0) : transpose(J_at_xp)
	end

	# normalization
	n = T(1)

	# we solve Jv + a σ1 = 0 with <b, v> = n
	# the solution is v = -σ1 J\a with σ1 = -n/<b, J\a>
	v, σ1, _, itv = pb.linbdsolver(J_at_xp, a, b, T(0), pb.zero, n)

	# we solve J'w + b σ2 = 0 with <a, w> = n
	# the solution is w = -σ2 J'\b with σ2 = -n/<a, J'\b>
		w, σ2, _, itw = pb.linbdsolver(JAd_at_xp, b, a, T(0), pb.zero, n)

	δ = T(1e-8)
	ϵ1, ϵ2, ϵ3 = T(δ), T(δ), T(δ)
	################### computation of σx σp ####################
	################### and inversion of Jfold ####################
	dpF = minus(F(x, set(par, lens, p + ϵ1)), F(x, set(par, lens, p - ϵ1))); rmul!(dpF, T(1) / T(2ϵ1))
	dJvdp = minus(apply(J(x, set(par, lens, p + ϵ3)), v), apply(J(x, set(par, lens, p - ϵ3)), v)); rmul!(dJvdp, T(1) / T(2ϵ3))
	σp = -dot(w, dJvdp) / n

	if hasHessian(pb) == false
		# We invert the jacobian of the Fold problem when the Hessian of x -> F(x, p) is not known analytically.
		# apply Jacobian adjoint
		u1 = applyJacobian(pb, x + ϵ2 * v, par0, w, true)
		u2 = apply(JAd_at_xp, w)
		σx = minus(u2, u1); rmul!(σx, 1 / ϵ2)
		########## Resolution of the bordered linear system ########
		# we invert Jfold
		dX, dsig, flag, it = pb.linbdsolver(J_at_xp, dpF, σx, σp, rhsu, rhsp)
	else
		# We invert the jacobian of the Fold problem when the Hessian of x -> F(x, p) is known analytically.
		# we solve it here instead of calling linearBorderedSolver because this removes the need to pass the linear form associated to σx
		# !!! Carefull, this method makes the linear system singular
		x1, x2, _, it = pb.linsolver(J_at_xp, rhsu, dpF)

		d2Fv = d2F(x, par0, x1, v)
		σx1 = -dot(w, d2Fv ) / n

		copyto!(d2Fv, d2F(x, par0, x2, v))
		σx2 = -dot(w, d2Fv ) / n

		dsig = (rhsp - σx1) / (σp - σx2)

		# dX = x1 .- dsig .* x2
		dX = _copy(x1); axpy!(-dsig, x2, dX)
	end

	if debugArray isa AbstractArray
		debugArray .= [J(x, par0) dpF ; σx' σp]
	end

	return dX, dsig, true, sum(it) + sum(itv) + sum(itw)
end

function (foldl::FoldLinearSolverMinAug)(Jfold, du::BorderedArray{vectype, T}; debugArray = nothing, kwargs...) where {vectype, T}
	# kwargs is used by AbstractLinearSolver
	out =  foldMALinearSolver((Jfold.x).u,
				 (Jfold.x).p,
				 Jfold.fldpb,
				 Jfold.params,
				 du.u, du.p;
				 debugArray = debugArray)
	# this type annotation enforces type stability
	return BorderedArray{vectype, T}(out[1], out[2]), out[3], out[4]
end

################################################################################################### Newton / Continuation functions
"""
$(SIGNATURES)

This function turns an initial guess for a Fold point into a solution to the Fold problem based on a Minimally Augmented formulation. The arguments are as follows
- `F   = (x, p) -> F(x, p)` where `p` is a set of parameters.
- `dF  = (x, p) -> d_xF(x, p)` associated jacobian
- `foldpointguess` initial guess (x_0, p_0) for the Fold point. It should be a `BorderedArray` as returned by the function `FoldPoint`
- `par` parameters used for the vector field
- `lens` parameter axis used to locate the Fold point.
- `eigenvec` guess for the 0 eigenvector
- `eigenvec_ad` guess for the 0 adjoint eigenvector
- `options::NewtonPar` options for the Newton-Krylov algorithm, see [`NewtonPar`](@ref).

# Optional arguments:
- `issymmetric` whether the Jacobian is Symmetric
- `Jᵗ = (x, p) -> transpose(d_xF(x, p))` jacobian adjoint, it should be implemented in an efficient manner. For matrix-free methods, `transpose` is not readily available and the user must provide a dedicated method. In the case of sparse based jacobian, `Jᵗ` should not be passed as it is computed internally more efficiently, i.e. it avoids recomputing the jacobian as it would be if you pass `Jᵗ = (x, p) -> transpose(dF(x, p))`
- `d2F = (x, p, v1, v2) ->  d2F(x, p, v1, v2)` a bilinear operator representing the hessian of `F`. It has to provide an expression for `d2F(x, p)[v1, v2]`.
- `normN = norm`
- `bdlinsolver` bordered linear solver for the constraint equation
- `kwargs` keywords arguments to be passed to the regular Newton-Krylov solver

# Simplified call
Simplified call to refine an initial guess for a Fold point. More precisely, the call is as follows

	newtonFold(F, J, br::AbstractBranchResult, ind_fold::Int; options = br.contparams.newtonOptions, kwargs...)

where the optional argument `Jᵗ` is the jacobian transpose and the Hessian is `d2F`. The parameters / options are as usual except that you have to pass the branch `br` from the result of a call to `continuation` with detection of bifurcations enabled and `index` is the index of bifurcation point in `br` you want to refine. You can pass newton parameters different from the ones stored in `br` by using the argument `options`.

!!! tip "Jacobian transpose"
    The adjoint of the jacobian `J` is computed internally when `Jᵗ = nothing` by using `transpose(J)` which works fine when `J` is an `AbstractArray`. In this case, do not pass the jacobian adjoint like `Jᵗ = (x, p) -> transpose(d_xF(x, p))` otherwise the jacobian will be computed twice!

!!! tip "ODE problems"
    For ODE problems, it is more efficient to pass the Bordered Linear Solver using the option `bdlinsolver = MatrixBLS()`
"""
function newtonFold(F, J,
				foldpointguess, par,
				lens::Lens,
				eigenvec, eigenvec_ad,
				options::NewtonPar;
				normN = norm,
				issymmetric::Bool = false,
				Jᵗ = nothing,
				d2F = nothing,
				bdlinsolver::AbstractBorderedLinearSolver = MatrixBLS(),
				kwargs...)

	foldproblem = FoldProblemMinimallyAugmented(
		F, J, Jᵗ, d2F,
		lens,
		_copy(eigenvec),
		_copy(eigenvec_ad), issymmetric,
		options.linsolver,
		# do not change linear solver if user provides it
		@set bdlinsolver.solver = (isnothing(bdlinsolver.solver) ? options.linsolver : bdlinsolver.solver))

	# Jacobian for the Fold problem
	Jac_fold_MA = (x, param) -> (x = x, params = param, fldpb = foldproblem)

	opt_fold = @set options.linsolver = FoldLinearSolverMinAug()

	# solve the Fold equations
	return newton(foldproblem, Jac_fold_MA, foldpointguess, par, opt_fold; normN = normN, kwargs...)..., foldproblem
end

function newtonFold(F, J,
				br::AbstractBranchResult, ind_fold::Int;
				normN = norm,
				issymmetric::Bool = false,
				Jᵗ = nothing,
				d2F = nothing,
				options = br.contparams.newtonOptions,
				nev = br.contparams.nev,
				startWithEigen = false,
				bdlinsolver::AbstractBorderedLinearSolver = MatrixBLS(),
				kwargs...)
	foldpointguess = FoldPoint(br, ind_fold)
	bifpt = br.specialpoint[ind_fold]
	eigenvec = bifpt.τ.u; rmul!(eigenvec, 1/normN(eigenvec))
	eigenvec_ad = _copy(eigenvec)

	if startWithEigen
		λ = zero(getvectoreltype(br))
		p = bifpt.param
		parbif = setParam(br, p)

		# jacobian at bifurcation point
		L = J(bifpt.x, parbif)

		# computation of zero eigenvector
		ζstar, = getAdjointBasis(L, λ, br.contparams.newtonOptions.eigsolver; nev = nev, verbose = false)
		eigenvec .= real.(ζstar)

		# computation of adjoint eigenvector
		_Jt = isnothing(Jᵗ) ? adjoint(L) : Jᵗ(bifpt.x, parbif)
		ζstar, = getAdjointBasis(_Jt, λ, br.contparams.newtonOptions.eigsolver; nev = nev, verbose = false)
		eigenvec_ad .= real.(ζstar)
		rmul!(eigenvec_ad, 1/normN(eigenvec_ad))
	end

	# solve the Fold equations
	return newtonFold(F, J, foldpointguess, br.params, br.lens, eigenvec, eigenvec_ad, options; issymmetric = issymmetric, Jᵗ = Jᵗ, d2F = d2F, normN = normN, bdlinsolver = bdlinsolver, kwargs...)
end

"""
$(SIGNATURES)

Codim 2 continuation of Fold points. This function turns an initial guess for a Fold point into a curve of Fold points based on a Minimally Augmented formulation. The arguments are as follows
- `F = (x, p) ->	F(x, p)` where `p` is a set of parameters
- `J = (x, p) -> d_xF(x, p)` associated jacobian
- `foldpointguess` initial guess (x_0, p1_0) for the Fold point. It should be a `BorderedArray` as returned by the function `FoldPoint`
- `par` set of parameters
- `lens1` parameter axis for parameter 1
- `lens2` parameter axis for parameter 2
- `eigenvec` guess for the 0 eigenvector at p1_0
- `eigenvec_ad` guess for the 0 adjoint eigenvector
- `options_cont` arguments to be passed to the regular [`continuation`](@ref)

# Optional arguments:
- `issymmetric` whether the Jacobian is Symmetric
- `Jᵗ = (x, p) -> transpose(d_xF(x, p))` associated jacobian transpose
- `d2F = (x, p, v1, v2) -> d2F(x, p, v1, v2)` this is the hessian of `F` computed at `(x, p)` and evaluated at `(v1, v2)`.
- `bdlinsolver` bordered linear solver for the constraint equation
- `updateMinAugEveryStep` update vectors `a, b` in Minimally Formulation every `updateMinAugEveryStep` steps
- `computeEigenElements = false` whether to compute eigenelements. If `options_cont.detecttEvent>0`, it allows the detection of ZH points.
- `kwargs` keywords arguments to be passed to the regular [`continuation`](@ref)

# Simplified call

	continuationFold(F, J, br::AbstractBranchResult, ind_fold::Int64, lens2::Lens, options_cont::ContinuationPar ; kwargs...)

where the parameters are as above except that you have to pass the branch `br` from the result of a call to `continuation` with detection of bifurcations enabled and `index` is the index of Fold point in `br` that you want to continue.

!!! tip "Jacobian tranpose"
    The adjoint of the jacobian `J` is computed internally when `Jᵗ = nothing` by using `transpose(J)` which works fine when `J` is an `AbstractArray`. In this case, do not pass the jacobian adjoint like `Jᵗ = (x, p) -> transpose(d_xF(x, p))` otherwise the jacobian would be computed twice!

!!! tip "ODE problems"
    For ODE problems, it is more efficient to pass the Bordered Linear Solver using the option `bdlinsolver = MatrixBLS()`

!!! tip "Detection of Bogdanov-Takens and Cusp bifurcations"
    In order to trigger the detection, pass `detectEvent = 1,2` in `options_cont`.
"""
function continuationFold(F, J,
				foldpointguess::BorderedArray{vectype, T}, par,
				lens1::Lens, lens2::Lens,
				eigenvec, eigenvec_ad,
				options_cont::ContinuationPar ;
				normC = norm,
				updateMinAugEveryStep = 0,
				issymmetric::Bool = false,
				Jᵗ = nothing,
				d2F = nothing,
				bdlinsolver::AbstractBorderedLinearSolver = MatrixBLS(),
			 	computeEigenElements = false,
				kwargs...) where {T, vectype}
	@assert lens1 != lens2 "Please choose 2 different parameters. You only passed $lens1"

	# options for the Newton Solver inheritated from the ones the user provided
	options_newton = options_cont.newtonOptions

	foldPb = FoldProblemMinimallyAugmented(
			F, J, Jᵗ, d2F,
			lens1,
			_copy(eigenvec),
			_copy(eigenvec_ad), issymmetric,
			options_newton.linsolver,
			# do not change linear solver if user provides it
			@set bdlinsolver.solver = (isnothing(bdlinsolver.solver) ? options_newton.linsolver : bdlinsolver.solver))

	# Jacobian for the Fold problem
	Jac_fold_MA = (x, param) -> (x = x, params = param, fldpb = foldPb)

	opt_fold_cont = @set options_cont.newtonOptions.linsolver = FoldLinearSolverMinAug()

	# this functions allows to tackle the case where the two parameters have the same name
	lenses = getLensSymbol(lens1, lens2)

	# global variables to save call back
	BT::T = one(T)
	CP::T = one(T)
	ZH::Int = 1

	# this function is used as a Finalizer
	# it is called to update the Minimally Augmented problem
	# by updating the vectors a, b
	function updateMinAugFold(z, tau, step, contResult; kUP...)
		# we first check that the continuation step was successful
		# if not, we do not update the problem with bad information!
		success = get(kUP, :state, nothing).isconverged
		(~modCounter(step, updateMinAugEveryStep) || success == false) && return true

		x = getVec(z.u)	# fold point
		p1 = getP(z.u)	# first parameter
		p2 = z.p	# second parameter
		newpar = set(par, lens1, p1)
		newpar = set(newpar, lens2, p2)

		a = foldPb.a
		b = foldPb.b

		# expression of the jacobian
		J_at_xp = foldPb.J(x, newpar)

		# compute new b
		newb = foldPb.linbdsolver(J_at_xp, a, b, T(0), foldPb.zero, T(1))[1]

		# compute new a
		if foldPb.issymmetric
			JAd_at_xp = J_at_xp
		else
			JAd_at_xp = hasAdjoint(foldPb) ? foldPb.Jᵗ(x, newpar) : transpose(J_at_xp)
		end
		newa = foldPb.linbdsolver(JAd_at_xp, b, a, T(0), foldPb.zero, T(1))[1]

		copyto!(foldPb.a, newa); rmul!(foldPb.a, 1/normC(newa))
		# do not normalize with dot(newb, foldPb.a), it prevents from BT detection
		copyto!(foldPb.b, newb); rmul!(foldPb.b, 1/normC(newb))

		# call the user-passed finalizer
		finaliseUser = get(kwargs, :finaliseSolution, nothing)
		if isnothing(finaliseUser) == false
			return finaliseUser(z, tau, step, contResult; prob = foldPb, kUP...)
		end
		return true
	end

	function testForBT_CP(iter, state)
		z = getx(state)
		x = getVec(z)		# fold point
		p1 = getP(z)		# first parameter
		p2 = getp(state)	# second parameter
		newpar = set(par, lens1, p1)
		newpar = set(newpar, lens2, p2)

		a = foldPb.a
		b = foldPb.b

		# expression of the jacobian
		J_at_xp = foldPb.J(x, newpar)

		# compute new b
		ζ = foldPb.linbdsolver(J_at_xp, a, b, T(0), foldPb.zero, T(1))[1]
		ζ ./= norm(ζ)

		# compute new a
		JAd_at_xp = hasAdjoint(foldPb) ? foldPb.Jᵗ(x, newpar) : transpose(J_at_xp)
		ζstar = foldPb.linbdsolver(JAd_at_xp, b, a, T(0), foldPb.zero, T(1))[1]
		ζstar ./= norm(ζstar)

		BT = dot(ζstar, ζ)
		CP = getP(state.τ)

		return BT, CP
	end

	function testForZH(iter, state)
		# printstyled(color = :magenta, "\n******* Entree dans Fold-Event-Dis\n")
		if isnothing(state.eigvals)
			ZH = 1
		else
			ϵ = iter.contParams.precisionStability
			ρ = minimum(abs ∘ real, state.eigvals)
			ZH = mapreduce(x -> ((real(x) > ρ) & (imag(x) > ϵ)), +, state.eigvals)
		end
		return ZH
	end

	# the following allows to append information specific to the codim 2 continuation to the user data
	_printsol = get(kwargs, :recordFromSolution, nothing)
	_printsol2 = isnothing(_printsol) ?
		(u, p; kw...) -> (zip(lenses, (getP(u), p))..., BT = BT) :
		(u, p; kw...) -> (; namedprintsol(_printsol(u, p;kw...))..., zip(lenses, (getP(u), p))..., BT = BT,)

	# eigen solver
	eigsolver = FoldEigsolver(getsolver(opt_fold_cont.newtonOptions.eigsolver))

	# solve the Fold equations
	br, u, tau = continuation(
		foldPb, Jac_fold_MA,
		foldpointguess, par, lens2,
		(@set opt_fold_cont.newtonOptions.eigsolver = eigsolver);
		linearAlgo = BorderingBLS(solver = opt_fold_cont.newtonOptions.linsolver, checkPrecision = false),
		kwargs...,
		normC = normC,
		recordFromSolution = _printsol2,
		finaliseSolution = updateMinAugFold,
		event = PairOfEvents(ContinuousEvent(2, testForBT_CP, computeEigenElements, ("bt", "cusp"), 0), DiscreteEvent(1, testForZH, false, ("zh",)))
		)
		@assert ~isnothing(br) "Empty branch!"
	return correctBifurcation(setproperties(br; type = :FoldCodim2, functional = foldPb)), u, tau
end

function continuationFold(F, J,
				br::AbstractBranchResult, ind_fold::Int,
				lens2::Lens,
				options_cont::ContinuationPar = br.contparams ;
				normC = norm,
				issymmetric::Bool = false,
				Jᵗ = nothing,
				d2F = nothing,
				nev = br.contparams.nev,
				startWithEigen = false,
				kwargs...)
	foldpointguess = FoldPoint(br, ind_fold)
	bifpt = br.specialpoint[ind_fold]
	eigenvec = bifpt.τ.u; rmul!(eigenvec, 1/norm(eigenvec))
	eigenvec_ad = _copy(eigenvec)

	p = bifpt.param
	parbif = setParam(br, p)

	if startWithEigen
		# jacobian at bifurcation point
		L = J(bifpt.x, parbif)

		# computation of adjoint eigenvalue
		eigenvec .= real.(	geteigenvector(options_cont.newtonOptions.eigsolver ,br.eig[bifpt.idx].eigenvec, bifpt.ind_ev))
		rmul!(eigenvec, 1/normC(eigenvec))

		# jacobian adjoint at bifurcation point
		_Jt = isnothing(Jᵗ) ? transpose(L) : Jᵗ(bifpt.x, parbif)

		ζstar, λstar = getAdjointBasis(_Jt, 0, br.contparams.newtonOptions.eigsolver; nev = nev, verbose = options_cont.newtonOptions.verbose)
		eigenvec_ad = real.(ζstar)
		rmul!(eigenvec_ad, 1/dot(eigenvec, eigenvec_ad))
	end

	return continuationFold(F, J, foldpointguess, parbif, br.lens, lens2, eigenvec, eigenvec_ad, options_cont ; issymmetric = issymmetric, Jᵗ = Jᵗ, d2F = d2F, normC = normC, kwargs...)
end

# structure to compute eigen-elements along branch of Fold points
struct FoldEigsolver{S} <: AbstractCodim2EigenSolver
	eigsolver::S
end

function (eig::FoldEigsolver)(Jma, nev; kwargs...)
	n = min(nev, length(Jma.x.u))
	J = Jma.fldpb.J(Jma.x.u, set(Jma.params, Jma.fldpb.lens, Jma.x.p))
	eigenelts = eig.eigsolver(J, n; kwargs...)
	return eigenelts
end

geteigenvector(eig::FoldEigsolver, vectors, i::Int) = geteigenvector(eig.eigsolver, vectors, i)


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

using Markdown
using InteractiveUtils

# ╔═╡ 72f2dc98-7923-11eb-3bf1-8d83a92ec5dc
begin
	using Markdown
	using InteractiveUtils
end

# ╔═╡ 443591a0-7925-11eb-198a-7557d7efb6b8
using Images

# ╔═╡ 75c648f4-7925-11eb-26fe-35998b00111f
using JImages

# ╔═╡ 839e25e4-7925-11eb-011c-09c7d65161de
md"reading Lena"

# ╔═╡ 92f5c24c-7925-11eb-1308-591e390d97ae
begin
	path = string(JImages.datadir, "/", "lena.h5")
	println("path to Lena is $path")
	lena = JImages.load_h5(path, "lena")
	Gray.(lena)
end

# ╔═╡ d46760b4-7925-11eb-21e4-db7c070b7c4c

# ╔═╡ b440870c-7925-11eb-3cb2-035863e491e2

# ╔═╡ Cell order:
# ╠═72f2dc98-7923-11eb-3bf1-8d83a92ec5dc
# ╠═443591a0-7925-11eb-198a-7557d7efb6b8
# ╠═75c648f4-7925-11eb-26fe-35998b00111f
# ╠═839e25e4-7925-11eb-011c-09c7d65161de
# ╠═92f5c24c-7925-11eb-1308-591e390d97ae
# ╠═d46760b4-7925-11eb-21e4-db7c070b7c4c
# ╠═b440870c-7925-11eb-3cb2-035863e491e2

max_stars_repo_name: jjgomezcadenas/JImages

max_issues_repo_name: jjgomezcadenas/JImages

max_forks_repo_name: jjgomezcadenas/JImages


text:
input_nums = zeros(Bool, 2021)

open("day01/day01.in") do io
    for line in readlines(io)
        num = parse(Int32, line)
        input_nums[parse(Int32, line)] = true
    end
end

for num in 1:1011
    if input_nums[num] && input_nums[2020-num]
        println(num * (2020-num))
    end
end

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code


text:
using JuMP, Base.Test, AmplNLWriter

## Solve test problem with non-linear binary variables
 #
 #  min   100 * (x2 - (0.5 + x1) ^ 2) ^ 2 + (1 - x1) ^ 2
 #  s.t.  x1, x2 binary
 #
 #  The solution is (0, 0).
 ##

# solver = AmplNLSolver(Ipopt.amplexe, ["print_level=0"])

@testset "example: jump_nonlinearbinary" begin
    m = Model(solver=solver)
    @variable(m, x[1:2], Bin)

    # Set some non-binary bounds on x1 and x2. These should be ignored.
    # The optimal solution if x is Int is (1, 2) which is allowed by these bounds
    setupperbound(x[1], 2)
    setupperbound(x[2], 2)

    @NLobjective(m, Min, 100*(x[2] - (0.5 + x[1])^2)^2 + (1 - x[1])^2)

    @test solve(m) == :Optimal

    if contains(getsolvername(solver), "ipopt")
        # Ipopt solves the relaxation
        @test isapprox(getvalue(x), [0.501245, 1.0], atol=1e-6)
        @test isapprox(getobjectivevalue(m), 0.249377, atol=1e-6)
    else
        @test getvalue(x)[:] == [0.0, 0.0]
        @test getobjectivevalue(m) == 7.25
    end
end


text:
#
# Operations on matrices using *random* data.
# Algorithms from "Finding structure with randomness: Probabilistic
# algorithms for constructing approximate matrix decompositions"
# by Nathan Halko, Per-Gunnar Martinsson, Joel A. Tropp.
# See http://arxiv.org/abs/0909.4061.

"""
    find_range_randomly(A, num_lhs, [q])

Find range of matrix `A` (with size m,n) using `num_lhs` random vectors and with
`q` power iterations (default q=1). Returns range matrix Q of size (m,num_lhs).

Based on Halko Martinsson & Tropp Algorithm 4.3.
"""
function find_range_randomly(A::AbstractMatrix, num_lhs::Integer, q=1)
    m,n = size(A)
    Omega = randn(n, num_lhs)
    Y = A*Omega
    for _ in 1:q
        Y = A' * Y
        Y = A * Y
    end
    Q,R = qr(Y)
    Q
end

"""
    find_svd_randomly(A, num_lhs, [q])

Compute a randomized SVD of `A` as *approximately* the matrix `U*W*V'` with
rank=`num_lhs`, where W is diagonal. Returns `(U, diag(W), V)` (note: *not* the
transpose of V). U and V will have only `num_lhs` columns, and similarly
`diag(W)` will be of length `num_lhs`. Uses `q` (default q=2) power iterations.

Based on Halko Martinsson & Tropp Algorithm 5.1.
"""
function find_svd_randomly(A::AbstractMatrix, num_lhs::Integer, q=2)
    Q = find_range_randomly(A, num_lhs, q)
    B = Q' * A
    u_b,w,v = svd(B)
    u = Q*u_b
    u,w,v
end


text:
"""
Author: Michael Palumbo
Created: December 2019
Contact: mlp95@psu.edu
Based on code by Alex Wise (aw@psu.edu)
Refactored and optimized by Eric Ford
"""

default_measure_width_at_frac_depth = 0.5
default_frac_of_width_to_fit = 0.5
default_init_guess_ccf_σ = 2000.0   # m/s

# Abstract types
""" Abstract type for functors to estimate the raidal velocitiy from a CCF and its velocity grid.  """
abstract type AbstractMeasureRvFromCCF end

"""
Estimate RV based on centroid velocity of the CCF.
Inputs:
* vels: Array of velocites where CCF was evaluated.
* ccf:  Array of values of CCF
"""
function (::AbstractMeasureRvFromCCF)(vels::A1, ccf::A2 ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1} } end

# Functions to be exported

"""
    `measure_rv_from_ccf(vels, ccf, [ccf_var]; alg )`
Return estimated RV based on the CCF using specified algorithm object.
Inputs:
* `vels`: Array of velocites where CCF was evaluated.
* `ccf`:  Array of values of CCF
Optional Arguements:
* `alg`: Functor specifying how to measure the RV and it's uncertainty from the CCF.  Options include:
MeasureRvFromCCFGaussian (default), MeasureRvFromCCFQuadratic, MeasureRvFromCCFCentroid, and MeasureRvFromMinCCF.
"""
function measure_rv_from_ccf  end

function measure_rv_from_ccf(vels::A1, ccf::A2, ; alg::AlgT=MeasureRvFromCCFGaussian() )  where {
                T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}, AlgT<:AbstractMeasureRvFromCCF }
    return alg(vels,ccf)
end

function measure_rv_from_ccf(vels::A1, ccf::A2, ccf_var::A3 ; alg::AlgT=MeasureRvFromCCFGaussian() )  where {
                T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,1}, AlgT<:AbstractMeasureRvFromCCF }
    return alg(vels,ccf,ccf_var)
end

function measure_rv_from_ccf(vels::A1, ccf::A2, ccf_var::A3 ; alg::AlgT=MeasureRvFromCCFGaussian() )  where {
                T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,2}, AlgT<:AbstractMeasureRvFromCCF }
    return alg(vels,ccf,ccf_var)
end

function measure_rv_from_ccf(vels::A1, ccf::A2, ccf_var::A3, covar_base::A4 ; alg::AlgT=MeasureRvFromCCFGaussian() )  where {
                T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,1}, T4<:Real, A4<:AbstractArray{T4,2}, AlgT<:AbstractMeasureRvFromCCF }
    return alg(vels,ccf,ccf_var,covar_base)
end

"""   `measure_rvs_from_ccf(vels, ccf, [ccf_var]; alg )`
At each time, return the estimated radial velocities based on the CCFs using the specified algorithm.
Inputs:
* `vels`: Array of velocites where CCF was evaluated.
* `ccf`:  Array of values of CCF
Optional Arguements:
* `alg`: Functor specifying how to measure the RV and its uncertainty from the CCF.  Options include:
MeasureRvFromCCFGaussian (default), MeasureRvFromCCFQuadratic, MeasureRvFromCCFCentroid, and MeasureRvFromMinCCF.
"""
function measure_rvs_from_ccf end

function measure_rvs_from_ccf(vels::A1, ccf::A2; alg::AlgT=MeasureRvFromCCFGaussian() )  where {
                T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,2}, AlgT<:AbstractMeasureRvFromCCF }
    nobs = size(ccf,2)
    rvs = zeros(nobs)
    σ_rvs = zeros(nobs)
    for i in 1:nobs
        result = measure_rv_from_ccf(vels, view(ccf,:,i), alg=alg)
        rvs[i] = result.rv
        σ_rvs[i] = result.σ_rv
    end
    return (rvs=rvs, σ_rvs=σ_rvs)
end

function measure_rvs_from_ccf(vels::A1, ccf::A2, ccf_var::A3; alg::AlgT=MeasureRvFromCCFGaussian() )  where {
                T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,2}, T3<:Real, A3<:AbstractArray{T3,2}, AlgT<:AbstractMeasureRvFromCCF }
    nobs = size(ccf,2)
    rvs = zeros(nobs)
    σ_rvs = zeros(nobs)
    for i in 1:nobs
        result = measure_rv_from_ccf(vels, view(ccf,:,i), view(ccf_var,:,i), alg=alg)
        rvs[i] = result.rv
        σ_rvs[i] = result.σ_rv
    end
    return (rvs=rvs, σ_rvs=σ_rvs)
end

function measure_rvs_from_ccf(vels::A1, ccf::A2, ccf_covar::A3; alg::AlgT=MeasureRvFromCCFGaussian() )  where {
                T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,2}, T3<:Real, A3<:AbstractArray{T3,3}, AlgT<:AbstractMeasureRvFromCCF }
    nobs = size(ccf,2)
    rvs = zeros(nobs)
    σ_rvs = zeros(nobs)
    for i in 1:nobs
        result = measure_rv_from_ccf(vels, view(ccf,:,i), view(ccf_covar,:,:,i), alg=alg)
        rvs[i] = result.rv
        σ_rvs[i] = result.σ_rv
    end
    return (rvs=rvs, σ_rvs=σ_rvs)
end

function measure_rvs_from_ccf(vels::A1, ccf::A2, ccf_var::A3, covar_base::A4; alg::AlgT )  where {
                T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,2}, T3<:Real, A3<:AbstractArray{T3,2}, T4<:Real, A4<:AbstractArray{T4,2}, AlgT<:AbstractMeasureRvFromCCF }
    nobs = size(ccf,2)
    rvs = zeros(nobs)
    σ_rvs = zeros(nobs)
    for i in 1:nobs
        result = measure_rv_from_ccf(vels, view(ccf,:,i), view(ccf_var,:,i), covar_base, alg=alg)
        rvs[i] = result.rv
        σ_rvs[i] = result.σ_rv
    end
    return (rvs=rvs, σ_rvs=σ_rvs)
end

# Utility functions

"""  `est_full_width(vels, ccf; measure_width_at_frac_depth = 0.5 )`
Return rough estimate of a ccf (or line) full width at the specified fractional depth (i.e., fraction of total depth).
This is based on the velocities of the first/last pixel to have a value less than the target value, with no interpolation.
Assumes vels is sorted.  Depth is measured assuming to the maximum value of ccf represents the continuum.
Could be improved for noisy data and segments with a steep slope due to the blaze or neighboring lines.
Inputs:
* `vels`: Array of velocites where CCF was evaluated.
* `ccf`:  Array of values of CCF
Optional Arguements:
* `measure_width_at_frac_depth`: What fractional CCF depth should be used for defining the width (0.5)
"""
function est_full_width(vels::A1, ccf::A2; measure_width_at_frac_depth::Real = default_measure_width_at_frac_depth ) where  {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1} }
    minccf, maxccf = extrema(ccf)
    depth = maxccf - minccf
    target_val = minccf + (1-measure_width_at_frac_depth) * depth
    ind1 = findfirst(ccf .<= target_val)
    ind2 = findlast(ccf .<= target_val)
    if isnothing(ind1) || isnothing(ind2)
        println("ccf = ",ccf)
        println("minccf= ", minccf, " maxccf= ", maxccf, " depth= ", depth, " measure_width_at_frac_depth= ", measure_width_at_frac_depth, " targetval= ",target_val, " ind1= ", ind1, " ind2= ", ind2)
        @error "est_full_width failed."
    end
    return vels[ind2] - vels[ind1]
end

"""  `find_idx_at_and_around_minimum(vels, ccf; frac_of_width, measure_width_at_frac_depth )`
Return the a pair with the index of the lowest value of ccf and a range of pixels surrounding it.
The range is based on finding the pixels with velocities nearest to the
Assumes vels is sorted.
Inputs:
* `vels`: Array of velocites where CCF was evaluated.
* `ccf`:  Array of values of CCF
Optional Arguements:
* `frac_of_width_to_fit`: How large of a range of velocities should be included in the fit (0.5)
* `measure_width_at_frac_depth`: What fractional CCF depth should be used for defining the width (0.5)
"""
function find_idx_at_and_around_minimum(vels::A1, ccf::A2; frac_of_width_to_fit::Real = default_frac_of_width_to_fit, measure_width_at_frac_depth::Real = default_measure_width_at_frac_depth) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1} }
    # do a prelim fit to get the width
    full_width = est_full_width(vels, ccf, measure_width_at_frac_depth=measure_width_at_frac_depth)

    # find the min and fit only that
    amin = argmin(ccf)
    if amin == 1 || amin==length(vels)
        offset = max(1,floor(Int64,length(vels)//4))
        amin = argmin(view(ccf,offset:(length(vels)-offset)))
        amin += offset-1
    end
    lend = vels[amin] - frac_of_width_to_fit * full_width
    rend = vels[amin] + frac_of_width_to_fit * full_width
    # get the indices
    lind = searchsortednearest(view(vels,1:amin), lend)
    rind = amin + searchsortednearest(view(vels,(amin+1):length(vels)), rend)
    inds = lind:rind

    return (amin, inds)
end


text:
push!(LOAD_PATH,"../src/")

using Splittings
using Documenter
using Literate
using Plots # to not capture precompilation output

examples = [
"Vlasov-Ampere"  => "examples/vlasov-ampere",
"Vlasov-Poisson" => "examples/vlasov-poisson",
"Bump On Tail"   => "examples/bump_on_tail",
"Rotation 2D"    => "examples/rotation2d_bsl",
"Vlasov-HMF"     => "examples/vlasov-hmf"
]

# generate examples

for example in examples

    EXAMPLE     = joinpath(@__DIR__, "..",  string(example[2],".jl"))
    DOC_OUTPUT  = joinpath(@__DIR__, "src", "examples")
    NB_OUTPUT   = joinpath(@__DIR__, "src", "notebooks")

    Literate.markdown(EXAMPLE, DOC_OUTPUT)
    Literate.notebook(EXAMPLE, NB_OUTPUT, execute=false)

end

makedocs(modules=[Splittings],
         doctest = false,
         format = :html,
         sitename = "Splittings.jl",
         pages = ["Introduction"    => "index.md",
                  "Semi-Lagrangian" => "bsl.md",
		      "Advection functions" => "advections.md",
                  "Examples" => [string(example[2],".md") for example in examples],
                  "User Documentation" => [
                      "How to Contribute" => "contributing.md"],
		          "Contents" => "contents.md"])

deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/pnavaro/Splittings.jl.git",
 )


text:
include("hold_action_policy.jl")
include("running_avg_policy.jl")

include("simple.jl")
include("sac.jl")
include("ppo.jl")
include("ddpg.jl")

max_stars_repo_name: albheim/ServiceMeshControl

max_issues_repo_name: albheim/ServiceMeshControl

max_forks_repo_name: albheim/ServiceMeshControl


text:
x = [1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6]
y = log(x) + sqrt(x)
gradients = "Not-Supplied"

using SchumakerSpline

########################
# Linear Extrapolation #
extrapolation = "Linear"
Spline, DerivSpline, SecondDerivSpline = schumaker(x,y, gradients, extrapolation)

#Pkg.add("PyPlot")
using PyPlot
#Pkg.update()
xArray = linspace(-5, 10, 100)
Spp = map(x -> DerivSpline(x), xArray)
plt = plot(xArray, Spp, color="green", linewidth=2.0, linestyle="--")
plt
Spp = map(x -> Spline(x), xArray)
plt = plot(x, y, color="blue", linewidth=2.0, linestyle="-")
plt
plt = plot(xArray, Spp, color="red", linewidth=2.0, linestyle="--")
plt
Spp = map(x -> DerivSpline(x), xArray)
plt = plot(xArray, Spp, color="green", linewidth=2.0, linestyle="--")
plt
Spp = map(x -> SecondDerivSpline(x), xArray)
plt = plot(xArray, Spp, color="Black", linewidth=2.0, linestyle="--")
plt

##########################
# Constant Extrapolation #
extrapolation = "Constant"
Spline, DerivSpline, SecondDerivSpline = schumaker(x,y, gradients, extrapolation)

#Pkg.add("PyPlot")
using PyPlot
#Pkg.update()
xArray = linspace(-5, 10, 100)
Spp = map(x -> Spline(x), xArray)
plt = plot(x, y, color="blue", linewidth=2.0, linestyle="-")
plt
plt = plot(xArray, Spp, color="red", linewidth=2.0, linestyle="--")
plt

##########################
# Curve    Extrapolation #
Spline, DerivSpline, SecondDerivSpline = schumaker(x,y, gradients)

#Pkg.add("PyPlot")
using PyPlot
#Pkg.update()
xArray = linspace(-5, 10, 100)
Spp = map(x -> Spline(x), xArray)
plt = plot(x, y, color="blue", linewidth=2.0, linestyle="-")
plt
plt = plot(xArray, Spp, color="red", linewidth=2.0, linestyle="--")
plt


text:
using Flux, Test
using Flux: maxpool, meanpool
using Flux: gradient

@testset "Pooling" begin
  x = randn(Float32, 10, 10, 3, 2)
  y = randn(Float32, 20, 20, 3, 2)
  ampx = AdaptiveMaxPool((5,5))
  @test ampx(x) == maxpool(x, PoolDims(x, 2))
  ampx = AdaptiveMeanPool((5,5))
  @test ampx(x) == meanpool(x, PoolDims(x, 2))
  ampy = AdaptiveMaxPool((10, 5))
  @test ampy(y) == maxpool(y, PoolDims(y, (2, 4)))
  ampy = AdaptiveMeanPool((10, 5))
  @test ampy(y) == meanpool(y, PoolDims(y, (2, 4)))
  gmp = GlobalMaxPool()
  @test size(gmp(x)) == (1, 1, 3, 2)
  gmp = GlobalMeanPool()
  @test size(gmp(x)) == (1, 1, 3, 2)
  mp = MaxPool((2, 2))
  @test mp(x) == maxpool(x, PoolDims(x, 2))
  mp = MeanPool((2, 2))
  @test mp(x) == meanpool(x, PoolDims(x, 2))
end

@testset "CNN" begin
  r = zeros(Float32, 28, 28, 1, 5)
  m = Chain(
    Conv((2, 2), 1=>16, relu),
    MaxPool((2,2)),
    Conv((2, 2), 16=>8, relu),
    MaxPool((2,2)),
    x -> reshape(x, :, size(x, 4)),
    Dense(288, 10), softmax)

  @test size(m(r)) == (10, 5)

  # Test bias switch
  bias = Conv(ones(Float32, 2, 2, 1, 3), ones(Float32, 3))
  ip = zeros(Float32, 28,28,1,1)

  op = bias(ip)
  @test sum(op) == prod(size(op))

  bias = Conv((2,2), 1=>3, bias = Flux.Zeros())
  op = bias(ip)
  @test sum(op) ≈ 0.f0
  gs = gradient(() -> sum(bias(ip)), Flux.params(bias))
  @test gs[bias.bias] == nothing

  # Train w/o bias and make sure no convergence happens
  # when only bias can be converged
  bias = Conv((2, 2), 1=>3, bias = Flux.Zeros());
  ip = zeros(Float32, 28,28,1,1)
  op = zeros(Float32, 27,27,3,1) .+ 2.f0
  opt = Descent()

  for _ = 1:10^3
    gs = gradient(params(bias)) do
      Flux.Losses.mse(bias(ip), op)
    end
    Flux.Optimise.update!(opt, params(bias), gs)
  end

  @test Flux.Losses.mse(bias(ip), op) ≈ 4.f0
end

@testset "asymmetric padding" begin
  r = ones(Float32, 28, 28, 1, 1)
  m = Conv((3, 3), 1=>1, relu; pad=(0,1,1,2))
  m.weight[:] .= 1.0
  m.bias[:] .= 0.0
  y_hat = m(r)[:,:,1,1]
  @test size(y_hat) == (27, 29)
  @test y_hat[1, 1] ≈ 6.0
  @test y_hat[2, 2] ≈ 9.0
  @test y_hat[end, 1] ≈ 4.0
  @test y_hat[1, end] ≈ 3.0
  @test y_hat[1, end-1] ≈ 6.0
  @test y_hat[end, end] ≈ 2.0
end

@testset "Depthwise Conv" begin
  r = zeros(Float32, 28, 28, 3, 5)
  m1 = DepthwiseConv((2, 2), 3=>15)
  @test size(m1(r), 3) == 15

  m3 = DepthwiseConv((2, 3), 3=>9)
  @test size(m3(r), 3) == 9

  # Test that we cannot ask for non-integer multiplication factors
  @test_throws AssertionError DepthwiseConv((2,2), 3=>10)
end

@testset "ConvTranspose" begin
  x = zeros(Float32, 28, 28, 1, 1)
  y = Conv((3,3), 1 => 1)(x)
  x_hat = ConvTranspose((3, 3), 1 => 1)(y)
  @test size(x_hat) == size(x)

  m = ConvTranspose((3,3), 1=>1)
  # Test that the gradient call does not throw: #900
  @test gradient(()->sum(m(x)), params(m)) isa Flux.Zygote.Grads
end

@testset "CrossCor" begin
  x = rand(Float32, 28, 28, 1, 1)
  w = rand(2,2,1,1)
  y = CrossCor(w, [0.0])

  @test isapprox(sum(w .* x[1:2, 1:2, :, :]), y(x)[1, 1, 1, 1], rtol=1e-7)

  r = zeros(Float32, 28, 28, 1, 5)
  m = Chain(
    CrossCor((2, 2), 1=>16, relu),
    MaxPool((2,2)),
    CrossCor((2, 2), 16=>8, relu),
    MaxPool((2,2)),
    x -> reshape(x, :, size(x, 4)),
    Dense(288, 10), softmax)

  @test size(m(r)) == (10, 5)
  @test y(x) != Conv(w, [0.0])(x)
  @test CrossCor(w[end:-1:1, end:-1:1, :, :], [0.0])(x) == Conv(w, [0.0])(x)
end

@testset "Conv with non quadratic window #700" begin
  data = zeros(Float32, 7,7,1,1)
  data[4,4,1,1] = 1

  l = Conv((3,3), 1=>1)
  expected = zeros(eltype(l.weight),5,5,1,1)
  expected[2:end-1,2:end-1,1,1] = l.weight
  @test expected ≈ l(data)

  l = Conv((3,1), 1=>1)
  expected = zeros(eltype(l.weight),5,7,1,1)
  expected[2:end-1,4,1,1] = l.weight
  @test expected ≈ l(data)

  l = Conv((1,3), 1=>1)
  expected = zeros(eltype(l.weight),7,5,1,1)
  expected[4,2:end-1,1,1] = l.weight
  @test expected ≈ l(data)

  @test begin
    # we test that the next expression does not throw
    randn(Float32, 10,10,1,1) |> Conv((6,1), 1=>1, Flux.σ)
    true
  end
end

@testset "conv output dimensions" begin
  m = Conv((3, 3), 3 => 16)
  @test Flux.outdims(m, (10, 10)) == (8, 8)
  m = Conv((3, 3), 3 => 16; stride = 2)
  @test Flux.outdims(m, (5, 5)) == (2, 2)
  m = Conv((3, 3), 3 => 16; stride = 2, pad = 3)
  @test Flux.outdims(m, (5, 5)) == (5, 5)
  m = Conv((3, 3), 3 => 16; stride = 2, pad = 3, dilation = 2)
  @test Flux.outdims(m, (5, 5)) == (4, 4)

  m = ConvTranspose((3, 3), 3 => 16)
  @test Flux.outdims(m, (8, 8)) == (10, 10)
  m = ConvTranspose((3, 3), 3 => 16; stride = 2)
  @test Flux.outdims(m, (2, 2)) == (5, 5)
  m = ConvTranspose((3, 3), 3 => 16; stride = 2, pad = 3)
  @test Flux.outdims(m, (5, 5)) == (5, 5)
  m = ConvTranspose((3, 3), 3 => 16; stride = 2, pad = 3, dilation = 2)
  @test Flux.outdims(m, (4, 4)) == (5, 5)

  m = DepthwiseConv((3, 3), 3 => 6)
  @test Flux.outdims(m, (10, 10)) == (8, 8)
  m = DepthwiseConv((3, 3), 3 => 6; stride = 2)
  @test Flux.outdims(m, (5, 5)) == (2, 2)
  m = DepthwiseConv((3, 3), 3 => 6; stride = 2, pad = 3)
  @test Flux.outdims(m, (5, 5)) == (5, 5)
  m = DepthwiseConv((3, 3), 3 => 6; stride = 2, pad = 3, dilation = 2)
  @test Flux.outdims(m, (5, 5)) == (4, 4)

  m = CrossCor((3, 3), 3 => 16)
  @test Flux.outdims(m, (10, 10)) == (8, 8)
  m = CrossCor((3, 3), 3 => 16; stride = 2)
  @test Flux.outdims(m, (5, 5)) == (2, 2)
  m = CrossCor((3, 3), 3 => 16; stride = 2, pad = 3)
  @test Flux.outdims(m, (5, 5)) == (5, 5)
  m = CrossCor((3, 3), 3 => 16; stride = 2, pad = 3, dilation = 2)
  @test Flux.outdims(m, (5, 5)) == (4, 4)

  m = MaxPool((2, 2))
  @test Flux.outdims(m, (10, 10)) == (5, 5)
  m = MaxPool((2, 2); stride = 1)
  @test Flux.outdims(m, (5, 5)) == (4, 4)
  m = MaxPool((2, 2); stride = 2, pad = 3)
  @test Flux.outdims(m, (5, 5)) == (5, 5)

  m = MeanPool((2, 2))
  @test Flux.outdims(m, (10, 10)) == (5, 5)
  m = MeanPool((2, 2); stride = 1)
  @test Flux.outdims(m, (5, 5)) == (4, 4)
  m = MeanPool((2, 2); stride = 2, pad = 3)
  @test Flux.outdims(m, (5, 5)) == (5, 5)
end

@testset "$ltype SamePad kernelsize $k" for ltype in (Conv, ConvTranspose, DepthwiseConv, CrossCor), k in ( (1,), (2,), (3,), (4,5), (6,7,8))
  data = ones(Float32, (k .+ 3)..., 1,1)
  l = ltype(k, 1=>1, pad=SamePad())
  @test size(l(data)) == size(data)

  l = ltype(k, 1=>1, pad=SamePad(), dilation = k .÷ 2)
  @test size(l(data)) == size(data)

  stride = 3
  l = ltype(k, 1=>1, pad=SamePad(), stride = stride)
  if ltype == ConvTranspose
    @test size(l(data))[1:end-2] == stride .* size(data)[1:end-2] .- stride .+ 1
  else
    @test size(l(data))[1:end-2] == ceil.(Int, size(data)[1:end-2] ./ stride)
  end
end

@testset "$ltype SamePad windowsize $k" for ltype in (MeanPool, MaxPool), k in ( (1,), (2,), (3,), (4,5), (6,7,8))
  data = ones(Float32, (k .+ 3)..., 1,1)

  l = ltype(k, pad=SamePad())
  @test size(l(data))[1:end-2] == ceil.(Int, size(data)[1:end-2] ./ k)
end


text:
# Diagonal Approximation algorithm

using LinearAlgebra
using SparseArrays
using DataStructures
using Preconditioners
using IncompleteLU
using Statistics
using IterativeSolvers

export diagonalapprox, tr_inv, diag_inv

## pchip start
struct _pchip
    N :: Integer
    xs :: Array{Float64}
    ys :: Array{Float64}
    ds :: Array{Float64}
end

ϕ(t) = 3t^2 - 2t^3
ψ(t) = t^3 - t^2

function _interp(pchip :: _pchip, x :: Number)
    i = _pchip_index(pchip, x)
    x1, x2 = pchip.xs[i:i+1]
    y1, y2 = pchip.ys[i:i+1]
    d1, d2 = pchip.ds[i:i+1]
    h = x2 - x1

    (y1 * ϕ((x2-x)/h)
     + y2 * ϕ((x-x1)/h)
     - d1*h * ψ((x2-x)/h)
     + d2*h * ψ((x-x1)/h))
end

function _pchip_index(pchip :: _pchip, x)
    N = pchip.N
    if N < 200  # Approximate performance cross-over on my old intel i7-3517U
        i = _pchip_index_linear_search(pchip, x)
    else
        i = _pchip_index_bisectional_search(pchip, x)
    end
    if i == N
        # Treat right endpoint as part of rightmost interval
        @assert (x ≈ pchip.xs[N])
        i = N-1
    end
    i
end

function _pchip_index_linear_search(pchip :: _pchip, x)
    xmin = pchip.xs[1]
    @assert (x >= xmin)

    i = 1
    N = pchip.N
    while i < N  &&  x >= pchip.xs[i+1]
        i = i + 1
    end
    i
end

function _pchip_index_bisectional_search(pchip :: _pchip, x)
    N = pchip.N
    imin, imax = 1, N
    xmin = pchip.xs[imin]
    xmax = pchip.xs[imax]
    @assert (x >= xmin && x <= xmax)

    i = imin + div(imax - imin + 1, 2)
    while imin < imax
        if x < pchip.xs[i]
            imax = i - 1
        elseif x >= pchip.xs[i+1]
            imin = i + 1
        else
            break
        end
        i = imin + div(imax - imin + 1, 2)
    end
    i
end

function _initial_ds_scipy(xs, ys)
    h(i) = xs[i+1]-xs[i]
    Δ(i) = (ys[i+1]-ys[i]) / h(i)

    N = length(xs)
    ds = similar(xs)
    if N == 2
        ds[:] .= Δ(1)
    else
        Δl = Δ(1)
        hl = h(1)
        for i ∈ 2:N-1
            Δr = Δ(i)
            hr = h(i)
            if sign(Δl) != sign(Δr) || Δl ≈ 0.0 || Δr ≈ 0.0
                ds[i] = 0.0
            else
                wl = 2hl + hr
                wr = hl + 2hr
                axx = (wl + wr)
                bxx = (wl/Δl + wr/Δr)
                ds[i] =  axx / bxx
            end
            Δl = Δr
            hl = hr
        end
        ds[1] = _edge_derivative(h(1), h(2), Δ(1), Δ(2))
        ds[N] = _edge_derivative(h(N-1), h(N-2), Δ(N-1), Δ(N-2))
    end
    ds
end

function _edge_derivative(h1, h2, Δ1, Δ2)
    d = ((2h1 + h2)*Δ1 - h2*Δ2) / (h1 + h2)
    if sign(d) != sign(Δ1)
        d = 0.0
    elseif sign(Δ1) != sign(Δ2)  &&  abs(d) > abs(3Δ1)
        d = 3Δ1
    end
    d
end

function interpolate(xs, ys)
    xs_ = [x for x ∈ xs]
    ys_ = [y for y ∈ ys]
    _assert_xs_ys(xs_, ys_)
    ds = _initial_ds_scipy(xs_, ys_)
    pchip = _pchip(length(xs_), xs_, ys_, ds)

    x -> _interp(pchip, x)
end

function _assert_xs_ys(xs, ys)
    N = length(xs)
    @assert (N > 1)
    @assert (N == length(ys))
    assert_monotonic_increase(xs)
end

function assert_monotonic_increase(xs)
    foldl((a,b) -> (@assert (a < b); b), xs)
end
## pcip end

# find approximate inverse
function z_approx_chol(A::AbstractMatrix, v)
    P = CholeskyPreconditioner(A)
    return (P.L \ v) , P
end
function z_approx_ilu(A::AbstractMatrix, v)
    P = ilu(A)
    L = P.L + I
    U = P.U
    return (U \ (L \ v)), P
end
function z_approx_amg(A::AbstractMatrix, v)
    P = AMGPreconditioner(A)
    z = P \ v
    return z, P
end

# Minimize this function for value of t
# i < t < j, t ∈ 𝐙
function t_bisection(M, i, j)
    vals = []
    for t in i+1:j-1
        push!(vals, abs((M[i] - M[j])*(i - j) - (M[i] - M[t])*(i-t) - (M[t] - M[j])*(t-j)))
    end
    return argmin(vals)+i
end

# Create new index in the middle of the longest L-R range
function middle_index!(S, Q = Nothing, Mₒ = Nothing)
    temp = 0
    L = 0
    R = 0
    interval = 0
    for i in S
        if temp != 0
            if interval < i - temp
                interval = i - temp
                L = temp
                R = i
            end
            temp = i
        else
            temp = i
        end
    end
    t = Int64(ceil((L + R)/2))
    push!(S, t)
    if Q != Nothing
        delete!(Q, (L, R))
        if t - L > 1
            enqueue!(Q, (L, t), abs(sum(Mₒ[L:t]) - (Mₒ[L] + Mₒ[t])*((t-L))/2))
        end
        if R - t > 1
            enqueue!(Q, (t, R), abs(sum(Mₒ[t:R]) - (Mₒ[t] + Mₒ[R])*((R-t))/2))
        end
    end
end

function point_identification(M, maxPts)
    N = size(M, 1)
    Mₒ = sort(M)
    J = sortperm(M)
    numSamples = 1
    initErr = abs(sum(Mₒ) - (Mₒ[1] + Mₒ[N])*(N/2))
    tempErr = initErr
    Sₒ = SortedSet{Int64}()
    push!(Sₒ, 1)
    push!(Sₒ, N)
    Q = PriorityQueue{Tuple{Int64, Int64}, Float64}(Base.Order.Reverse)
    enqueue!(Q, (1, N), initErr)
    while numSamples < maxPts && peek(Q)[2] > 0.001*initErr
        # pop interval (L, R) with the largest error from Q
        interval = dequeue!(Q)
        L, R = interval[1], interval[2]

        # for t = L+1:R-1 find bisection index t
        t = t_bisection(Mₒ, L, R)

        # Add t to Sₒ and increase numsamples
        if t - L >= 1 && R - t >= 1
            push!(Sₒ, t)
            numSamples = numSamples + 1
        end

        # Push interval (L, t) and (t, R) with temperror into Q
        if t - L > 1
            enqueue!(Q, (L, t), abs(sum(Mₒ[L:t]) - (Mₒ[L] + Mₒ[t])*((t-L)/2)))
        end
        if R - t > 1
            enqueue!(Q, (t, R), abs(sum(Mₒ[t:R]) - (Mₒ[t] + Mₒ[R])*((R-t)/2)))
        end

        # insert midpoint into largest interval after every 5 samples
        if numSamples%5 == 0
            middle_index!(Sₒ, Q, Mₒ)
            numSamples = numSamples + 1
        end
    end

    # insert midpoints in the largest intervals
    while numSamples < maxPts
        middle_index!(Sₒ)
        numSamples = numSamples + 1
    end

    ## unique Mₒ values and Sₒ indices
    temp = Set()
    S = []
    for i in Sₒ
        if (Mₒ[i] in temp) == false
            push!(S, J[i])
            push!(temp, Mₒ[i])
        end
    end
    return S
end

function basis_vector!(e, k)
    e .= Ref(0)
    e[k] = one(eltype(e))
end

function linreg(X, Y)
    n = size(Y, 1)
    sX = sum(X)
    sY = sum(Y)
    sXY = sum(X.*Y)
    sX2 = sum(X.^2)

    c = (sY*sX2 - sX*sXY)/(n*sX2 - sX^2)
    b = (n*sXY - sX*sY)/(n*sX2 - sX^2)
    return b, c
end

function linear_model(S, D, M, n)
    b, c = linreg(M[S], D)
    Tf = zero(eltype(D))
    Ms = Vector{eltype(D)}(undef, n)
    for i in 1:n
        Ms[i] = b * M[i] + c
        Tf += b * M[i] + c
    end
    return Tf, Ms
end

function pchip_iterpolation(S, D, M, n)
    pchip = interpolate(M[S], D)
    Tf = zero(eltype(D))
    Ms = Vector{eltype(D)}(undef, n)
    for i in 1:n
        Ms[i] = pchip(M[i])
        Tf += pchip(M[i])
    end
    return Tf, Ms
end

function diagonalapprox(A::AbstractMatrix, n::Int64, p::Int64, pc, model)
    # Compute M = diag(approximation of A⁻¹)
    v = Matrix{eltype(A)}(undef, size(A,1), n)
    rademacherDistribution!(v)
    if pc == "ilu"
        Z, pl = z_approx_ilu(A, v)
    elseif pc == "amg"
        Z, pl = z_approx_chol(A, v)
    elseif pc == "cheby"
        M = chebydiagonal(A, 4, 6)
        pl = ilu(A)
    else
        Z, pl = z_approx_chol(A, v)
    end

    if pc != "cheby"
        M = vec(mean(v .* Z , dims=2))
    end

    # Compute fitting sample S, set of k indices
    S = point_identification(M, p)

    # Solve for D(i) = e'(i) A⁻¹ e(i) for i = 1 ... k
    e = zeros(size(A, 1))
    D = eltype(A)[]
    if pc == "cheby"
        lmin, lmax = lczeigen(A, invfun, rademacherDistribution!)
        for i in 1:size(S, 1)
            basis_vector!(e,S[i])
            push!(D, e' * chebyshev(A, e, lmin, lmax, Pl = pl))
        end
    else
        for i in 1:size(S, 1)
            basis_vector!(e, S[i])
            push!(D, e' * cg(A, e, Pl = pl))
        end
    end

    # Obtain fitting model f(M) ≈ D by fitting f(M(S)) to D(S)
    Tf = zero(eltype(A))
    Ms = Vector{eltype(A)}(undef, size(A, 1))
    if model == "pchip"
        Tf, Ms = pchip_iterpolation(S, D, M, size(A, 1))
    else
        Tf, Ms = linear_model(S, D, M, size(A, 1))
    end
    return (Tf, Ms)
end
"""
    tr_inv(A::AbstractMatrix, n::Int64, p::Int64; pc = "chol", model = "linear")

Diagonal Approximation algorithm for calculating trace of the matrix inverse.

# Arguments
- `A` : Symmetric Positive Definite Matrix
- `n` : Probing vector count for initial approximation.
- `p` : Sample Points count for interpolation
- `pc` : Preconditioner used for initial approximation and cg ("chol" - Incomplete Cholesky, "ilu" - IncompleteLU, "amg" - AlgebraicMultigrid, "cheby" - Chebyshev Approximation). Default = "chol"
- `model` : Fitting model used for calculation of trace ("linear" - Linear Regression, "pchip" - PCHIP interpolation). Default = "linear".
"""
function tr_inv(A::AbstractMatrix, n::Int64, p::Int64; pc = "chol", model = "linear")
    return diagonalapprox(A, n, p, pc, model)[1];
end

"""
    diag_inv(A::AbstractMatrix, n::Int64, p::Int64; pc = "chol", model = "linear")

Diagonal Approximation algorithm for the matrix inverse.

# Arguments
- `A` : Symmetric Positive Definite Matrix
- `n` : Probing vector count for initial approximation.
- `p` : Sample Points count for interpolation
- `pc` : Preconditioner used for initial approximation and cg ("chol" - Incomplete Cholesky, "ilu" - IncompleteLU, "amg" - AlgebraicMultigrid, "cheby" - Chebyshev Approximation). Default = "chol"
- `model` : Fitting model used for calculation of trace ("linear" - Linear Regression, "pchip" - PCHIP interpolation). Default = "linear".
"""
function diag_inv(A::AbstractMatrix, n::Int64, p::Int64; pc = "chol", model = "linear")
    return diagonalapprox(A, n, p, pc, model)[2];
end


text:
using OneDimensionalCorrelation
##
println("Calculate CI by fist calculating CI on smaller element blocks")
println("and then combining lowest 2 CI vectors to form a new CI calculation.")
println("The result is then compared to full CI calculation.")
##
b4 = BasisLobatto(-12, 12, 4, 24)
b6 = BasisLobatto(-12, 12, 6, 16)
b8 = BasisLobatto(-12, 12, 8, 12)

# Nuclear potential
V(x) = -2( exp(-0.5(x-2)^2) + exp(-0.5(x+2)^2) )

##
println("\n *** Hartree Fock calculations ***\n")

E4, Eorb4, orbitals4 = solve_hartree_fock(b4, V)
E6, Eorb6, orbitals6 = solve_hartree_fock(b6, V)
E8, Eorb8, orbitals8 = solve_hartree_fock(b8, V)

##

println("\n *** CI in blocks calculations ***\n")

eb4 = block_ci(b4, orbitals4, V; nstates=2)
eb6 = block_ci(b6, orbitals6, V; nstates=2)
eb8 = block_ci(b8, orbitals8, V; nstates=2)

##
println("\n *** Full CI calculation ***\n")
ci = full_ci(b6, orbitals6, V)

##

E_corr = ci["energies"][1] - E6

Eb4_corr = eb4["energy"][1] - E4
Eb6_corr = eb6["energy"][1] - E6
Eb8_corr = eb8["energy"][1] - E8

E2 = ci["energies"][2]

E2_err4 = abs( ( eb4["energy"][2] - E2 ) / E2 )
E2_err6 = abs( ( eb6["energy"][2] - E2 ) / E2 )
E2_err8 = abs( ( eb8["energy"][2] - E2 ) / E2 )

##
println("\nHartree-Fock energy: ", E6)
println("Total correlation energy (ground state): ", E_corr, "\n")

println("Relative correlation energy based on number of elements when calculating with blocks:\n")
println("4 elements ", round(100 * (Eb4_corr/E_corr); digits=2),"% of total correlation")
println("6 elements ", round(100 * (Eb6_corr/E_corr); digits=2),"% of total correlation")
println("8 elements ", round(100 * (Eb8_corr/E_corr); digits=2),"% of total correlation")

println("\n***\n")

println("First excited state energy FCI: ", E2,"\n")
println("Excited state energies when using element blocks:")
println("4 elements: ", eb4["energy"][2], "  relative accuracy: ", round(E2_err4; sigdigits=2))
println("6 elements: ", eb6["energy"][2], "  relative accuracy: ", round(E2_err6; sigdigits=2))
println("8 elements: ", eb8["energy"][2], "  relative accuracy: ", round(E2_err8; sigdigits=2))

max_stars_repo_name: tjjarvinen/OneDimensionalCorrelation

max_issues_repo_name: tjjarvinen/OneDimensionalCorrelation

max_forks_repo_name: tjjarvinen/OneDimensionalCorrelation


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1ObjectFieldSelector <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    fieldPath::Any # spec type: Union{ Nothing, String } # spec name: fieldPath

    function IoK8sApiCoreV1ObjectFieldSelector(;apiVersion=nothing, fieldPath=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ObjectFieldSelector, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiCoreV1ObjectFieldSelector, Symbol("fieldPath"), fieldPath)
        setfield!(o, Symbol("fieldPath"), fieldPath)
        o
    end
end # type IoK8sApiCoreV1ObjectFieldSelector

const _property_map_IoK8sApiCoreV1ObjectFieldSelector = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("fieldPath")=>Symbol("fieldPath"))
const _property_types_IoK8sApiCoreV1ObjectFieldSelector = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("fieldPath")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1ObjectFieldSelector }) = collect(keys(_property_map_IoK8sApiCoreV1ObjectFieldSelector))
Swagger.property_type(::Type{ IoK8sApiCoreV1ObjectFieldSelector }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ObjectFieldSelector[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ObjectFieldSelector }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ObjectFieldSelector[property_name]

function check_required(o::IoK8sApiCoreV1ObjectFieldSelector)
    (getproperty(o, Symbol("fieldPath")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ObjectFieldSelector }, name::Symbol, val)
end


text:
module BasicInterpolators

using Base: RefValue
using LinearAlgebra: mul!, dot
using Memoize

include("base.jl")
include("polynomial.jl")
include("spline.jl")
include("chebyshev.jl")
include("scattered.jl")

end


text:
using Combinatorics

numbers = open("xmas_cypher.txt") do f
    parse.(Int, readlines(f))
end

function find_weak(numbers::Array{Int64, 1}, preamble::Int64)::Int64
    for i in (preamble+1):length(numbers)
        found = false
        n = numbers[i]

        for pair in combinations(numbers[(i-preamble):(i-1)], 2)
            if sum(pair) == n
                found = true
                break
            end
        end

        if !found
            return n
        end
    end

    return 0
end

function find_contiguous_set(numbers::Array{Int64,1}, sum_to::Int64)::Array{Int64, 1}

    for i in 1:(length(numbers)-1)
        for j in (i+1):length(numbers)
            if sum(numbers[i:j]) == sum_to
                return numbers[i:j]
            end
        end
    end

    return []
end

weak = find_weak(numbers, 25)
println(weak)
set = find_contiguous_set(numbers, weak)
println(minimum(set) + maximum(set))


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

using Markdown
using InteractiveUtils

# ╔═╡ 078beb8e-3003-11eb-2024-832e809585c7
md"The goal of this notebook is to simultaneously learn how to program with julia while learning the material for AP Calculus, AP Chemistry, and computer science. In each section contains another notebook where all the notebooks for each subject can be accessed. To reach this, click on the linked here below. "

# ╔═╡ 12f4af38-3003-11eb-192f-6f16e8048709
md"Calculus notebooks are [here](/open?path=./calculus/intro.jl)"

# ╔═╡ 3c91c06a-3003-11eb-2eae-07b252c6c016
md"Chemistry notebooks are [here](/open?path=./chemistry/intro.jl)"

# ╔═╡ 49a7f12a-3003-11eb-01c7-efc80a813ddb
md"Computer Science notebooks are [here](/open?path=./julia/intro.jl)"

# ╔═╡ 8dfe479c-3004-11eb-3d02-3f83cb926b36
md"Add a new notebook [here](/new)"

# ╔═╡ Cell order:
# ╟─078beb8e-3003-11eb-2024-832e809585c7
# ╟─12f4af38-3003-11eb-192f-6f16e8048709
# ╟─3c91c06a-3003-11eb-2eae-07b252c6c016
# ╟─49a7f12a-3003-11eb-01c7-efc80a813ddb
# ╟─8dfe479c-3004-11eb-3d02-3f83cb926b36


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct NetworkInterfaceIPConfigurationListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{NetworkInterfaceIPConfiguration} } # spec name: value
    nextLink::Any # spec type: Union{ Nothing, String } # spec name: nextLink

    function NetworkInterfaceIPConfigurationListResult(;value=nothing, nextLink=nothing)
        o = new()
        validate_property(NetworkInterfaceIPConfigurationListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(NetworkInterfaceIPConfigurationListResult, Symbol("nextLink"), nextLink)
        setfield!(o, Symbol("nextLink"), nextLink)
        o
    end
end # type NetworkInterfaceIPConfigurationListResult

const _property_map_NetworkInterfaceIPConfigurationListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"), Symbol("nextLink")=>Symbol("nextLink"))
const _property_types_NetworkInterfaceIPConfigurationListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{NetworkInterfaceIPConfiguration}", Symbol("nextLink")=>"String")
Base.propertynames(::Type{ NetworkInterfaceIPConfigurationListResult }) = collect(keys(_property_map_NetworkInterfaceIPConfigurationListResult))
Swagger.property_type(::Type{ NetworkInterfaceIPConfigurationListResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NetworkInterfaceIPConfigurationListResult[name]))}
Swagger.field_name(::Type{ NetworkInterfaceIPConfigurationListResult }, property_name::Symbol) =  _property_map_NetworkInterfaceIPConfigurationListResult[property_name]

function check_required(o::NetworkInterfaceIPConfigurationListResult)
    true
end

function validate_property(::Type{ NetworkInterfaceIPConfigurationListResult }, name::Symbol, val)
end


text:
using Autotest

# We need this import since there are conflicting describe methods in DataFrames and StatsBase...
import Autotest: describe

if length(ARGS) > 0 && ARGS[1] == "continuous"
  Autotest.autorun("PerfTest", "src", "test")
else
  Autotest.run("PerfTest", "src", "test")
end