text:
module config

# This is the main configuration file for the whole of the kokako-detect project.
# Any constants and parameters that need to be changed must be changed here only.
# Likewise the addition of any such parameters should also be added here.
# Author: Tim Evans
# Date: 30/07/14

#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
# paths and constants pertaining to extracting features
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------

export wavdir, traindir, speclength, path_rel, test_size, specwidth, labelleddir,datafile,
ntrees, ratio, subsampleratio, modeldir, specdir,boxwidth,boxheight,thresh,kernwidth,kernsigma,
minarea

#the following paths are all relative to the project directory and will be concatenated
#To use absolute paths sets path_rel = false
path_rel = true
wavdir = "data/wavs/"
traindir = "data/training/"
specdir = "data/spectrograms/"
labelleddir = "data/labelled/"
datafile = "data.csv"

#------------------------------------------------------------------------------
#                           filter parameters
#------------------------------------------------------------------------------
fc = 1250
npoles = 4
speclength = 10 #seconds
specwidth = 930 #setting the width of spectrograms produced

#------------------------------------------------------------------------------
#                       feature extraction paramaters
#------------------------------------------------------------------------------
boxwidth = 3
boxheight = 10
# This determines how much we will subsample the noise and labelled noise in our
# dataset. E.g. if we have 30 observations of bird call pixels, and subsampleratio = 0.5
# then we will add 15 labelled noise, and 15 unlabelled noise to the dataset
subsampleratio = 1

#------------------------------------------------------------------------------
#                          Random Forest parameters
#------------------------------------------------------------------------------
modeldir = "data/model/"
ntrees = 10;
ratio = 0.2;

#------------------------------------------------------------------------------
#                prediction parameters (post-processing)
#------------------------------------------------------------------------------
thresh = 0.6
kernwidth = 5 # Gaussian kernel parameters
kernsigma = 2
minarea = 16 # minimum sized syllable to consider as a bird

end

max_stars_repo_name: kokakodetect/kokako-detection

max_issues_repo_name: kokakodetect/kokako-detection

max_forks_repo_name: kokakodetect/kokako-detection


text:
using Documenter, Literate
ENV["PLOTS_DEFAULT_BACKEND"] = "GR"
ENV["GKSwstype"]="nul"
using Plots
using Plots.PlotMeasures: mm

default(
    fontfamily = "Computer Modern", linewidth=1.5,
    titlefontsize=(16+8), guidefontsize=(11+5),
    tickfontsize=(8+4), legendfontsize=(8+4),
    left_margin=5mm, right_margin=5mm)

src = joinpath(@__DIR__, "src")
lit = joinpath(@__DIR__, "src")

config = Dict(
    "credit" => false,  # credit is configured to render in Documenter instead
    "repo_root_url"=> "https://github.com/simonsobs/PSpipe/blob/master/project/Planck_cov/",
)

nonexecute_config = copy(config)
nonexecute_config["codefence"] = "```julia" => "```"
execution_exclusion = []

mkpath(joinpath(@__DIR__, "build"))
cp(joinpath(@__DIR__, "..", "src", "util.jl"), joinpath(@__DIR__, "build", "util.jl"); force=true)
cp(joinpath(@__DIR__, "..", "example.toml"), joinpath(@__DIR__, "build", "example.toml"); force=true)
cp(joinpath(@__DIR__, "..", "input"), joinpath(@__DIR__, "input"); force=true)

for (root, _, files) ∈ walkdir(joinpath(@__DIR__, "../src")), file ∈ files
    splitext(file)[2] == ".jl" || continue
    ipath = joinpath(root, file)
    opath = lit
    if file ∉ execution_exclusion
        Literate.markdown(ipath, opath; config=config)
    else
        @warn("not executing unless forced for ", file)
        Literate.markdown(ipath, opath; config=nonexecute_config)
    end
end

## manually clean and add back files
rm(joinpath(@__DIR__, "build"); force=true, recursive=true)
mkpath(joinpath(@__DIR__, "build"))
cp(joinpath(@__DIR__, "..", "src", "util.jl"), joinpath(@__DIR__, "build", "util.jl"); force=true)
cp(joinpath(@__DIR__, "..", "example.toml"), joinpath(@__DIR__, "build", "example.toml"); force=true)
cp(joinpath(@__DIR__, "..", "input"), joinpath(@__DIR__, "input"); force=true)

makedocs(
    sitename = "PSPipe Planck",
    modules = Module[],
    clean=false,

    format=Documenter.HTML(
        prettyurls=get(ENV, "CI", nothing) == "true",
        assets=["assets/so.css"],
    ),
    pages = [
        "Introduction" => "index.md",
        "setup.jl" => "setup.md",
        "util.jl" => "util.md",
        "rawspectra.jl" => "rawspectra.md",
        "fitnoisemodel.jl" => "fitnoisemodel.md",
        "signalsim.jl" => "signalsim.md",
        "corrections.jl" => "corrections.md",
        "whitenoise.jl" => "whitenoise.md",
        "covmat.jl" => "covmat.md",
        "slurmgen.jl" => "slurmgen.md",

        # "Pipeline" => "pipeline.md"
        ]
    )

ext: jl

max_forks_repo_name: sgiardie/PSpipe


text:
using Documenter, WiSER

makedocs(
    format = Documenter.HTML(),
    sitename = "WiSER.jl",
    authors = "Hua Zhou, Chris German",
    clean = true,
    debug = true,
    pages = [
        "index.md",
        "model_fitting.md",
        "simulation.md"
    ]
)

deploydocs(
    repo   = "github.com/OpenMendel/WiSER.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
function is_at_widest_option_type(x)
    types = Set()
    for v in x
        union!(types, Set((typeof(v),)))
    end
    lt = length(types)
    lt > 2 && return false
    lt > 1 && return ((Nothing in types) || (Missing in types) ? true : false)
    true
end
export is_at_most_option_type

function get_types_from_iterable(x)
    types = Set()
    for v in x
        union!(types, Set((typeof(v),)))
    end
    collect(types)
end

function narrow_column_types!(data)
    for k in keys(data)
        types = get_types_from_iterable(data[k])
        TYPE = length(types) > 1 ? Union{types...} : types[1]
        data[k] = convert(Vector{TYPE}, data[k])
    end
end

to_df(data::Dict{String,Any}) = DataFrame(Any[values(data)...], [Symbol(k) for k in keys(data)])

function get_data(results, column_names)
    data = Dict{String, Any}()
    for column_name in column_names
        data[column_name] = []
        cnlp = column_name * "_logprob"
        data[cnlp] = []
        for (ix, trace) in enumerate(results.traces)
            kt = keys(trace)
            if column_name in kt
                push!(data[column_name], trace[column_name].value)
                push!(data[cnlp], trace[column_name].logprob_sum)
            else
                push!(data[column_name], missing)
                push!(data[cnlp], missing)
            end
        end
    end
    data
end

function to_df(f; params = (), num_iterations = 100, method::InferenceType = FORWARD, inference_params = Dict())
    check_sf(f)
    inference_params["num_iterations"] = num_iterations
    results = inference(f, method; params = params, inference_params = inference_params)
    column_names = addresses(results)
    types_cnames = get_types_from_iterable(column_names)
    !((length(types_cnames) == 1) && (types_cnames[1] == String)) && error("$f does not have all string addresses")
    data = get_data(results, column_names)
    !all(map(v -> is_at_widest_option_type(v), values(data))) && error("Each column type is not typed at widest Option[T]")
    narrow_column_types!(data)
    to_df(data)
end
export to_df

abstract type TabularModel{F} end

struct DataFrameModel{F} <: TabularModel{F}
    f::F
    table::DataFrame
end
function DataFrameModel(f::F; params = (), num_iterations = 100, method::InferenceType = FORWARD, inference_params = Dict()) where F
    df = to_df(f; params = params, num_iterations = num_iterations, method = method, inference_params = inference_params)
    DataFrameModel(f, df)
end

@doc raw"""
    struct SQLModel{F} <: TabularModel{F}
        f::F
        db::SQLite.DB
    end

A representation of the model `f` as a SQL table.
"""
struct SQLModel{F} <: TabularModel{F}
    f::F
    db::SQLite.DB
end

@doc raw"""
    function SQLModel(f::F; params = (), num_iterations = 100, method::InferenceType = FORWARD, inference_params = Dict(), name = nothing) where F

Create a representation of the model `f` as a database table which can be queried using SQL (SQLite, specifically).
This method works by first conducting inference on the model (including forward sampling if the user just wants to query the prior),
and then collecting sampled values in a table. The resulting table is titled `<name>_results` and includes a column for each address
(titled `a`) that occurs in the union of all sampled traces and a corresponding column for each address's log probability (titled `a_logprob`).

This methood *does* work with open-universe structured models; traces (rows in the table) that do not include address `a` have a
`NULL` in columns `a` and `a_logprob` in that row.

+ `params`: additional parameters to pass to the model during inference
+ `num_iterations`: number of iterations of inference to complete; the meaning is dependent on the inference algorithm used
+ `method`: an `InferenceType`, defaults to `Forward()` for forward sampling.
+ `inference_params`: a dict of parameters to pass to the inference method.
+ `name`: a name to use instead of the model's name, which is used by default.
"""
function SQLModel(f::F; params = (), num_iterations = 100, method::InferenceType = FORWARD, inference_params = Dict(), name = nothing) where F
    df = to_df(f; params = params, num_iterations = num_iterations, method = method, inference_params = inference_params)
    db = SQLite.DB()
    name === nothing && (name = "$(f)")
    tablename = df |> SQLite.load!(db, "$(name)_results")
    @info "Created table $tablename"
    SQLModel(f, db)
end
export DataFrameModel, SQLModel

@doc raw"""
    function query(sqm::SQLModel, q::String)

Execute the SQLite query `q` against the `SQLModel` and return results as a `DataFrame`.
"""
function query(sqm::SQLModel, q::String)
    DBInterface.execute(sqm.db, q) |> DataFrame
end
export query


text:
using Plots
using Formatting

function logistic(r,x)
    r*x*(1.0-x)
end

function logisticmap(;x0=0.2,r=4,n=10)
    results = [(x0,logistic(r,x0))];
    for i in 1:(n-1)
      xold = results[end][2]
      xnew = logistic(r,xold)
      append!(results,[(xold,xold)])
      append!(results,[(xold,xnew)])
  end
  results
end

function vary(dname;x0_vals=[],r_vals=[], nvals=[])
 if !isdir(dname)
  mkdir(dname)
 end
 fno = 1
 for x0 in x0_vals
  for r in r_vals
      for n in nvals
          fname = format(dname*"/frame{:05d}.png",fno)
          fno = fno + 1
          simres = logisticmap(x0=x0, r=r, n=n)
          xvals = [x for (x, y) in simres]
          yvals = [y for (x, y) in simres]
          plot(xvals,yvals,legend=false,xaxis=false,yaxis=false,xlim=(0,1),ylim=(0,1));
          savefig(fname)
        end
    end
end
end

# name/number each iteration then plot 1, 1+2, 1+2+3 ...

# To make frames varying r, specify x0_vals as an array with a single value
# along with whatever values of r you want to loop through
vary("frames_r";x0_vals=[0.5],r_vals=0:0.1:5, n=[100])

# To make frames varying x0, specify r_vals as an array with a single value
# along with whatever values of x0 you want to loop through
vary("frames_x0";x0_vals=0:0.05:1,r_vals=[4.0], n=[100])

vary("frames_n";x0_vals=[0.2],r_vals=[4], nvals = 1:500)

r_vals = 4.1
vary("loop3";x0_vals=[(r_vals+sqrt(r_vals-4)*sqrt(r_vals))/(2*r_vals)],r_vals=[4], nvals=1:10)

vary("loop2";x0_vals=[(r_vals-sqrt(r_vals-4)*sqrt(r_vals))/(2*r_vals)],r_vals=[4], nvals=1:10)

ffmpeg -i frame%05d.png -vf "fps=12,split[s0][s1];[s0]palettegen[p];[s1][p]paletteuse" -loop 0 output.gif


text:
_string(::EickerHuberWhite) = "Eicker White Huber variance"

Base.@kwdef struct NaiveLocalLinearRD{K,B,V<:VarianceEstimator} <: SharpRD
    kernel::K
    bandwidth::B
    variance::V = EickerHuberWhite()
end

Base.@kwdef struct FittedLocalLinearRD{R,F,K,B,S,T,C}
    rdd_setting::R
    fitted_lm::F
    fitted_kernel::K
    fitted_bandwidth::B
    data_subset::S
    tau_est::T
    se_est::T
    coeftable::C
end

linearweights(fitted::FittedLocalLinearRD) = linearweights(fitted.fitted_lm)

function linearweights(fitted_lm::RegressionModel; idx = 2)
    wts = fitted_lm.model.rr.wts
    γs = (fitted_lm.model.pp.chol\fitted_lm.model.pp.X'*Diagonal(wts))[idx, :]
    γs
end

function var(::EickerHuberWhite, fitted_lm::RegressionModel)
    γs = linearweights(fitted_lm)
    dot(abs2.(γs), abs2.(residuals(fitted_lm)))
end

fit(method::SharpRD, ZsR::RunningVariable, Y) = fit(method, RDData(Y, ZsR))

function fit(method::NaiveLocalLinearRD, rddata::RDData; level=0.95)
    c = rddata.cutoff
    kernel = method.kernel
    variance = method.variance

    h = bandwidth(method.bandwidth, kernel, rddata)
    fitted_kernel = setbandwidth(kernel, h)

    lb = support(fitted_kernel).lb
    ub = support(fitted_kernel).ub

    new_support = Interval(lb + c, ub + c)

    rddata_filt = rddata[new_support]
    wts = weights(fitted_kernel, rddata_filt.ZsR)

    fitted_lm = fit(LinearModel, @formula(Ys ~ Ws * ZsC), rddata_filt, wts = wts)

    tau_est = coef(fitted_lm)[2]
    se_est = sqrt(var(variance, fitted_lm))
    γs = linearweights(fitted_lm)

    z = tau_est/se_est
    pval = 1-cdf(Normal(), abs(z))
    z_quantile = quantile(Normal(), 1-level/2)
    ci = Interval(tau_est - se_est*z_quantile, tau_est + se_est*z_quantile)
    # as in GLM.jl
    levstr = isinteger(level*100) ? string(Integer(level*100)) : string(level*100)

    res = [h tau_est se_est "unaccounted" z pval leftendpoint(ci) rightendpoint(ci)]
    colnms = ["h"; "τ̂"; "se"; "bias"; "z"; "p-val"; "Lower $levstr%"; "Upper $levstr%"]
    rownms = ["Sharp RD estimand"]
    coeftbl = CoefTable(res, colnms, rownms, 6, 5)

    FittedLocalLinearRD(
        rdd_setting = method,
        fitted_lm = fitted_lm,
        fitted_kernel = fitted_kernel,
        fitted_bandwidth = h,
        data_subset = rddata_filt,
        tau_est = tau_est,
        se_est = se_est,
        coeftable = coeftbl,
    )
end

function Base.show(io::IO, rdd_fit::RegressionDiscontinuity.FittedLocalLinearRD)
    println(io, "Local linear regression for regression discontinuity design")
    println(io, "       ⋅⋅⋅⋅ " * "Naive inference (not accounting for bias)")
    println(io, "       ⋅⋅⋅⋅ " * _string(rdd_fit.rdd_setting.kernel))
    println(io, "       ⋅⋅⋅⋅ " * _string(rdd_fit.rdd_setting.bandwidth))
    println(io, "       ⋅⋅⋅⋅ " * _string(rdd_fit.rdd_setting.variance))
    Base.show(io, rdd_fit.coeftable)
end

@recipe function f(rdd_fit::FittedLocalLinearRD; show_local_support = false)
    sorted_subset = rdd_fit.data_subset |> DataFrame |> df -> sort(df, [:Zs])
    sorted_preds = predict(rdd_fit.fitted_lm, sorted_subset)
    if show_local_support
        lb = support(rdd_fit.fitted_kernel).lb
        ub = support(rdd_fit.fitted_kernel).ub
        @series begin
            label := nothing
            seriestype := :vline
            linecolor := :lightgrey
            inestyle := :dot
            linewidth := 2
            [lb; ub]
        end
    end
    seriestype = :path
    linecolor --> :lightblue
    linewidth --> 2
    label --> "Local linear fit"
    sorted_subset.Zs, sorted_preds
end


text:
module HMCResearchRandomWalks

module ViewOut
using DataFrames, CSV, Dates
include(string(@__DIR__, "/view_out/current.jl"))
include(string(@__DIR__, "/view_out/data.jl"))
include(string(@__DIR__, "/view_out/plot_noflow.jl"))
export plotdata_noflow,
       savedata,
       savetofolder,
       loaddata,
       Current, current_df, current_seed, current_path, current_file
end

module Enz
using Random, DataFrames, Reexport, ProgressMeter
@reexport using HMCResearchRandomWalks.ViewOut
include(string(@__DIR__, "/step_size/PARAMETERS_enz.jl"))
include(string(@__DIR__, "/step_size/locations_enz.jl"))
include(string(@__DIR__, "/walk_logic/walk_logic.jl"))
include(string(@__DIR__, "/arlett/flow_arlett.jl"))
include(string(@__DIR__, "/arlett/spawn_arlett.jl"))
include(string(@__DIR__, "/step_size/enz_sim.jl"))
export enz_sim,
       getdata_enzstepsize,
       set_NUMBER_OF_WALKS,
       set_MAX_STEPS_PER_WALK,
       set_THICK_ENZ,
       set_ENZ_STEP_SIZE
end

module PPD
using Random, DataFrames, Reexport, ProgressMeter
@reexport using HMCResearchRandomWalks.ViewOut
include(string(@__DIR__, "/step_size/PARAMETERS_ppd.jl"))
include(string(@__DIR__, "/step_size/locations_ppd.jl"))
include(string(@__DIR__, "/walk_logic/walk_logic.jl"))
include(string(@__DIR__, "/arlett/flow_arlett.jl"))
include(string(@__DIR__, "/step_size/ppd_sim.jl"))
export ppd_sim,
       getdata_ppdstepsize,
       set_NUMBER_OF_WALKS,
       set_MAX_STEPS_PER_WALK,
       set_PPD_ON,
       set_PPD_STEP_SIZE
end

module Arlett
using Random, DataFrames, Reexport, ProgressMeter
@reexport using HMCResearchRandomWalks.ViewOut
include(string(@__DIR__, "/arlett/arlett.jl"))
include(string(@__DIR__, "/walk_logic/walk_logic.jl"))
export arlett_sim,
       getdata_arlett_animation,
       getdata_arlett_noflow,
       set_NUMBER_OF_WALKS,
       set_MAX_STEPS_PER_WALK,
       set_FLOW_BIAS,
       set_CATALASE_ON_WALLS
end

module Quinto
using Random, DataFrames, Reexport, ProgressMeter
@reexport using HMCResearchRandomWalks.ViewOut
include(string(@__DIR__, "/quinto/quinto_sim.jl"))
export quinto_sim,
       set_NUMBER_OF_WALKS,
       set_MAX_STEPS_PER_WALK
end

end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.1
#   kernelspec:
#     display_name: Julia 1.1.0
#     language: julia
#     name: julia-1.1
# ---

# # Streamfunction For Vertically Integrated Transport
#
# This notebook proceeds in several steps:
#
# - read `TrspX,TrspY` transport vector field from file and compute its convergence.
# - split `TrspX,TrspY` into rotational and divergent components (i.e., `Helmholtz Decomposition`).
# - compute streamfunction from the non-divergent component.
#
# _Note: this requires the `MeshArrays.jl v0.2.3` (upcoming release)_

# +
using MeshArrays, Plots

include(joinpath(dirname(pathof(MeshArrays)),"Plots.jl"))
include("prepare_transports.jl")

#1) get grid
!isdir("GRID_LLC90") ? run(`git clone https://github.com/gaelforget/GRID_LLC90`) : nothing
mygrid=GridSpec("LLC90"); GridVariables=GridLoad(mygrid);

#2) get transport
(TrspX, TrspY, TauX, TauY, SSH)=trsp_read(mygrid,"GRID_LLC90/")

#3) compute convergence
TrspCon=convergence(TrspX,TrspY);

#4) apply mask & avoid isolated Canyons / singular matrices
msk=1.0 .+ 0.0 * mask(view(GridVariables["hFacC"],:,1),NaN,0.0)
TrspCon=msk*TrspCon;
# -

# ### Helmholtz Decomposition

# +
#scalar potential
TrspPot=ScalarPotential(TrspCon)

#plot result
heatmap(TrspPot,clims=(-5e5,5e5))

#Divergent vector field
(TrspXdiv,TrspYdiv)=gradient(TrspPot,GridVariables)
TrspXdiv=TrspXdiv.*GridVariables["DXC"]
TrspYdiv=TrspYdiv.*GridVariables["DYC"];

# Check that divergences match:
tmpCon=convergence(TrspXdiv,TrspYdiv)
heatmap(tmpCon-TrspCon) #TrspCon-tmpCon should be << TrspCon

#Rotational vector field
TrspXrot = TrspX-TrspXdiv
TrspYrot = TrspY-TrspYdiv

#vector Potential
TrspPsi=VectorPotential(TrspX,TrspY,GridVariables)
heatmap(1e-6*msk*TrspPsi)
# -

max_stars_repo_name: lmilechin/MeshArrayNotebooks

max_issues_repo_name: lmilechin/MeshArrayNotebooks

max_forks_repo_name: lmilechin/MeshArrayNotebooks


text:
using Random
using Ranges

indices = randperm(MersenneTwister(1234), 8000)[1:20]
# rates = [2000, 4000, 6000, 8000, 10000, 14000, 18000, 22000, 26000, 30000, 34000, 38000, 42000, 46000, 50000]
rates = [8000]
d_range = [0.1, 0.5, 1, 2, 4, 8, 12]
d_range = convert(Array{Float64,1}, d_range)
max_iter = 10
mode = "untargeted"  # "untargeted" or "targeted"
optim = "joint"  # "univariate" or "joint"
save_imp = false
save_normacc_graph = true
test_name = "joint_norm_acc"  # univariate_norm_acc

fw_arguments = (
    line_search=FrankWolfe.MonotonousNonConvexStepSize(),
    max_iteration=max_iter,
    print_iter=max_iter / 10,
    verbose=true,
    #lazy=true,
)


text:
import Printf
import JSON
import Healpix
import Quaternions
import LinearAlgebra: dot, ×
import Base: show

export MINUTES_PER_DAY,
       DAYS_PER_YEAR,
       rpm2angfreq,
       angfreq2rpm,
       period2rpm,
       rpm2period,
       ScanningStrategy,
       update_scanning_strategy,
       load_scanning_strategy,
       to_dict,
       save,
       time2pointing,
       time2pointing!,
       genpointings,
       genpointings!

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

rpm2angfreq(rpm) = 2π * rpm / 60
angfreq2rpm(ω) = ω / 2π * 60

@doc raw"""
    rpm2angfreq(rpm)
    angfreq2rpm(ω)

Convert rotations per minute into angular frequency 2πν (in Hertz), and vice
versa.
"""
rpm2angfreq, angfreq2rpm

period2rpm(p) = 60 / p
rpm2period(rpm) = 60 / rpm

@doc raw"""
    period2rpm(p)
    rpm2perriod(rpm)

Convert a period (time span) in seconds into a number of rotations per minute,
and vice versa.
"""
period2rpm, rpm2period

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

const MINUTES_PER_DAY = 60 * 24
const DAYS_PER_YEAR = 365.25

@doc raw"""
The structure `ScanningStrategy` encodes the information needed to build a
set of pointing directions. It contains the following fields:

- `omega_spin_hz`: angular speed of the rotation around the spin axis (equal to 2πν)
- `omega_prec_hz`: angular speed of the rotation around the precession axis (equal to 2πν)
- `omega_year_hz`: angular speed of the rotation around the Elliptic axis (equal to 2πν)
- `omega_hwp_hz`: angular speed of the rotation of the half-wave plate (equal to 2πν)
- `spinsunang_rad`: angle between the spin axis and the Sun-Earth direction
- `borespinang_rad`: angle between the boresight direction and the spin axis
- `q1`, `q3`: quaternions used to generate the pointings

Each field has its measure unit appended to the name. For instance, field
`spinsunang_rad` must be expressed in radians.

"""
struct ScanningStrategy
    omega_spin_hz::Float64
    omega_prec_hz::Float64
    omega_year_hz::Float64
    omega_hwp_hz::Float64
    spinsunang_rad::Float64
    borespinang_rad::Float64
    # First quaternion used in the rotation
    q1::Quaternion
    # Third quaternion used in the rotation
    q3::Quaternion

    ScanningStrategy(; spin_rpm = 0,
        prec_rpm = 0,
        yearly_rpm = 1  / (MINUTES_PER_DAY * DAYS_PER_YEAR),
        hwp_rpm = 0,
        spinsunang_rad = deg2rad(45.0),
        borespinang_rad = deg2rad(50.0),) = new(rpm2angfreq(spin_rpm),
        rpm2angfreq(prec_rpm),
        rpm2angfreq(yearly_rpm),
        rpm2angfreq(hwp_rpm),
        spinsunang_rad,
        borespinang_rad,
        compute_q1(borespinang_rad),
        compute_q3(spinsunang_rad))

    ScanningStrategy(io::IO) = load_scanning_strategy(io)
    ScanningStrategy(filename::AbstractString) = load_scanning_strategy(filename)
end

compute_q1(borespinang_rad) = qrotation_y(borespinang_rad)
compute_q3(spinsunang_rad) = qrotation_y(π / 2 - spinsunang_rad)

@doc raw"""
    update_scanning_strategy(sstr::ScanningStrategy)

Update the internal fields of a `ScanningStrategy` object. If you change any of the
fields in a `ScanningStrategy` object after it has been created using the constructors,
call this function
before using one of the functions `time2pointing`, `time2pointing!`, `genpointings`,
and `genpointings!`, as they rely on a number of internal parameters that need to be
updated.

```julia
sstr = ScanningStrategy()
# ... use sstr ...

sstr.borespinang_rad *= 0.5
update_scanning_strategy(sstr)
```

"""
function update_scanning_strategy(sstr::ScanningStrategy)
    sstr.q1 = compute_q1(sstr.borespinang_rad)
    sstr.q3 = compute_q3(sstr.spinsunang_rad)
end

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

function load_scanning_strategy(io::IO)
    data = JSON.Parser.parse(io)

    sstr_data = data["scanning_strategy"]

    ScanningStrategy(spin_rpm = sstr_data["spin_rpm"],
        prec_rpm = sstr_data["prec_rpm"],
        yearly_rpm = sstr_data["yearly_rpm"],
        hwp_rpm = sstr_data["hwp_rpm"],
        spinsunang_rad = sstr_data["spinsunang_rad"],
        borespinang_rad = sstr_data["borespinang_rad"])
end

function load_scanning_strategy(filename::AbstractString)
    open(filename) do inpf
        load_scanning_strategy(inpf)
    end
end

@doc raw"""
    load_scanning_strategy(io::IO) -> ScanningStrategy
    load_scanning_strategy(filename::AbstractString) -> ScanningStrategy

Create a `ScanningStrategy` object from the definition found in the JSON file
`io`, or from the JSON file with path `filename`. See also
[`load_scanning_strategy`](@ref).
"""
load_scanning_strategy

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

@doc raw"""
    to_dict(sstr::ScanningStrategy) -> Dict{String, Any}

Convert a scanning strategy into a dictionary suitable to be serialized using
JSON or any other structured format. See also [`save`](@ref).
"""
function to_dict(sstr::ScanningStrategy)
    Dict("scanning_strategy" => Dict("spin_rpm" => angfreq2rpm(sstr.omega_spin_hz),
        "prec_rpm" => angfreq2rpm(sstr.omega_prec_hz),
        "yearly_rpm" => angfreq2rpm(sstr.omega_year_hz),
        "hwp_rpm" => angfreq2rpm(sstr.omega_hwp_hz),
        "spinsunang_rad" => sstr.spinsunang_rad,
        "borespinang_rad" => sstr.borespinang_rad))
end

function save(io::IO, sstr::ScanningStrategy)
    print(io, JSON.json(to_dict(sstr), 4))
end

function save(filename::AbstractString, sstr::ScanningStrategy)
    open(filename, "w") do outf
        save(outf, sstr)
    end
end

@doc raw"""
    save(io::IO, sstr::ScanningStrategy)
    save(filename::AbstractString, sstr::ScanningStrategy)

Write a definition of the scanning strategy in a self-contained JSON file.
You can reload this definition using one of the constructors of
[`ScanningStrategy`](@ref).
"""
save

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

function show(io::IO, sstr::ScanningStrategy)
    Printf.@printf(io, """Scanning strategy:
    Spin angular velocity.................................... %g rot/s
    Precession angular velocity.............................. %g rot/s
    Yearly angular velocity around the Sun................... %g rot/s
    Half-wave plate angular velocity......................... %g rot/s
    Angle between the spin axis and the Sun-Earth direction.. %f°
    Angle between the boresight direction and the spin axis.. %f°
""",
        sstr.omega_spin_hz,
        sstr.omega_prec_hz,
        sstr.omega_year_hz,
        sstr.omega_hwp_hz,
        rad2deg(sstr.spinsunang_rad),
        rad2deg(sstr.borespinang_rad))
end

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

function time2pointing!(sstr::ScanningStrategy, time_s, beam_dir, polangle_rad, resultvec)
    curpolang = mod2pi(polangle_rad + 4 * sstr.omega_hwp_hz * time_s)
    # The polarization vector lies on the XY plane; if polangle_rad=0 then
    # the vector points along the X direction at t=0.
    poldir = StaticArrays.SVector(cos(curpolang), sin(curpolang), 0.0)

    q2 = qrotation_z(sstr.omega_spin_hz * time_s)
    q4 = qrotation_x(sstr.omega_prec_hz * time_s)
    q5 = qrotation_z(sstr.omega_year_hz * time_s)

    qtot = q5 * (q4 * (sstr.q3 * (q2 * sstr.q1)))
    rot = rotationmatrix_normalized(qtot)
    # Direction in the sky of the beam main axis
    resultvec[4:6] = rot * beam_dir
    # Direction in the sky of the beam polarization axis
    poldir = rot * poldir

    # The North for a vector v is just -dv/dθ, as θ is the
    # colatitude and moves along the meridian
    (θ, ϕ) = Healpix.vec2ang(resultvec[4:6]...)
    northdir = StaticArrays.SVector(-cos(θ) * cos(ϕ), -cos(θ) * sin(ϕ), sin(θ))

    cosψ = clamp(dot(northdir, poldir), -1, 1)
    crosspr = northdir × poldir
    sinψ = clamp(sqrt(dot(crosspr, crosspr)), -1, 1)
    resultvec[3] = atan(cosψ, sinψ)

    resultvec[1], resultvec[2] = θ, ϕ
end

function time2pointing(sstr::ScanningStrategy, time_s, beam_dir, polangle_rad)
    resultvec = Float64[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    time2pointing!(sstr, time_s, beam_dir, polangle_rad, resultvec)
    resultvec
end

@doc raw"""
    time2pointing!(sstr::ScanningStrategy, time_s, beam_dir, polangle_rad, resultvec)
    time2pointing(sstr::ScanningStrategy, time_s, beam_dir, polangle_rad)

Calculate the pointing direction of a beam along the direction `beam_dir`, with
a detector sensitive to the polarization along the angle `polangle_rad`. The
result is saved in `resultvec` for `time2pointing!`, and it is the return value
of `time2pointing`; it is a 6-element array containing the following fields:

1. The colatitude (in radians) of the point in the sky
2. The longitude (in radians) of the point in the sky
3. The polarization angle, in the reference frame of the sky
4. The X component of the normalized pointing vector
5. The Y component
6. The Z component

Fields #4, #5, #6 are redundant, as they can be derived from the colatitude
(field #1) and longitude (field #2). They are returned as the code already
computes them.

The vector `beam_dir` and the angle `polangle_rad` must be expressed in the
reference system of the focal plane. If `polangle_rad == 0`, the detector
measures polarization along the x axis of the focal plane. The normal direction
to the focal plane is along the z axis; thus, the boresight director is such
that `beam_dir = [0., 0., 1.]`.

"""
time2pointing!, time2pointing

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

function genpointings!(sstr::ScanningStrategy, timerange_s, beam_dir, polangle_rad, result)
    @assert size(result)[1] == length(timerange_s)

    @inbounds for (idx, t) in enumerate(timerange_s)
        time2pointing!(sstr, t, beam_dir, polangle_rad, view(result, idx, :))
    end
end

function genpointings(sstr::ScanningStrategy, timerange_s, beam_dir, polangle_rad)
    result = Array{Float64}(undef, length(timerange_s), 6)
    genpointings!(sstr, timerange_s, beam_dir, polangle_rad, result)

    result
end

@doc raw"""
    genpointings!(sstr::ScanningStrategy, timerange_s, beam_dir, polangle_rad, result)
    genpointings(sstr::ScanningStrategy, timerange_s, beam_dir, polangle_rad)

Generate a set of pointing directions and angles for a given orientation
`beam_dir` (a 3-element vector) of the boresight beam, assuming the scanning
strategy in `sstr`. The pointing directions are calculated over all the elements
of the list `timerange_s`.

The two functions only differ in the way the result is returned to the caller.
Function `genpointings` returns a N×6 matrix containing the following fields:

1. The colatitude (in radians)
2. The longitude (in radians)
3. The polarization angle (in radians)
4. The X component of the one-length pointing vector
5. The Y component
6. The Z component

Function `genpointings!` works like `genpointings`, but it accept a
pre-allocated matrix as input (the `result` parameter) and will save the result
in it. The matrix must have two dimensions with size `(N, 6)` at least.

Both functions are simple iterators wrapping [`time2pointing!`](@ref) and
[`time2pointing`](@ref).

"""
genpointings!, genpointings


text:
module PathDistribution

# package code goes here

include("misc.jl")
include("monte_carlo.jl")
include("path_enumeration.jl")

export
	PathSample,
	monte_carlo_path_sampling,
	monte_carlo_path_number,
	estimate_cumulative_count,
	PathEnum,
	path_enumeration,
	actual_cumulative_count

end # module


text:
module Doug

" Compute Shannon entropy from a probability distribution "
function entropy(dist)
    -reduce(+, 0, map(xlog2x, dist))
end

function xlog2x(x)
    x * log2(x)
end

# using Gadfly
# Plot the entropy in the case of two possibilities with probabilities p and (1-p)
# from mathematical theory of communication page 50.
# plot([p -> Doug.entropy([p, 1-p])], 0.0, 1.0)

end

max_stars_repo_name: drichardson/examples

max_issues_repo_name: drichardson/examples

max_forks_repo_name: drichardson/examples


text:
include("../Models21.jl")
import Knet
using Knet.Layers21, Knet.Ops21
using Knet.Train20: param
using Knet.Ops20: softmax
using Knet.KnetArrays: KnetArray
using CUDA: CuArray
using PyCall, AutoGrad, SHA, Tar, DelimitedFiles
using ImageCore, ImageTransformations, FileIO, Artifacts, Base.Threads, Base.Iterators

# pytorch models are saved in e.g. /home/dyuret/.cache/torch/hub/checkpoints/resnext50_32x4d-7cdf4587.pth
torch = pyimport("torch")
nn = pyimport("torch.nn")
models = pyimport("torchvision.models")
transforms = pyimport("torchvision.transforms")
PIL = pyimport("PIL")

t2a(x) = x.cpu().detach().numpy()
chkparams(a,b)=((pa,pb)=params.((a,b)); length(pa)==length(pb) && all(isapprox.(pa,pb)))
rdims(a)=permutedims(a, ((ndims(a):-1:1)...,))
typename(p::PyObject) = pytypeof(p).__name__

"""
    torchimport("resnet18", "resnet18", ResNet)
    torchimport("mobilenet_v2", "mobilenet_v2_100_224_pt", MobileNetV2)
    torchimport("mobilenet_v3_large", "mobilenet_v3_large_100_224_pt", MobileNetV3)
"""
function torchimport(model, file=nothing, testmodel=nothing)
    global pm, am, km, cm, dm, ax, ay
    @assert haskey(models, model)
    pm = getproperty(models,model)(pretrained=true).eval()
    px = randn(Float32, 224, 224, 3, 1)
    py = px |> rdims |> torch.tensor |> pm |> t2a |> rdims

    T = Float32 ## use Float64 for @gcheck

    save_type = Knet.array_type[]
    Knet.array_type[] = Array{T}
    am = torch2knet(pm)
    ax = Knet.atype(px)
    ay = am(ax)
    @show ay ≈ py
    #ap = Param(ax)
    #@show @gcheck am(ap) (nsample=3,)

    Knet.array_type[] = KnetArray{T}
    km = torch2knet(pm)
    @show chkparams(km,am)
    kx = Knet.atype(px)
    ky = km(kx)
    @show ky ≈ py
    #kp = Param(kx)
    #@show @gcheck km(kp) (nsample=3,)

    Knet.array_type[] = CuArray{T}
    cm = torch2knet(pm)
    @show chkparams(cm,am)
    cx = Knet.atype(px)
    cy = cm(cx)
    @show cy ≈ py
    #cp = Param(cx)
    #@show @gcheck cm(cp) (nsample=3,)

    if file !== nothing
        @info "Saving $(file).jld2"
        saveweights("$(file).jld2", cm; atype=CuArray)

        if testmodel !== nothing
            @info "Loading $(file).jld2"
            dm = testmodel(file; pretrained=false)
            setweights!(dm, "$(file).jld2")
            @show chkparams(dm,am)
            dx = Knet.atype(px)
            dy = dm(dx)
            @show dy ≈ py
            #dp = Param(dx)
            #@show @gcheck dm(dp) (nsample=3,)
        end

        run(`tar cf $(file).tar $(file).jld2`)
        sha1 = Tar.tree_hash("$(file).tar")
        @info "git-tree-sha1 = \"$sha1\""

        run(`gzip $(file).tar`)
        sha2 = open("$(file).tar.gz") do f; bytes2hex(sha256(f)); end
        @info "sha256 = \"$sha2\""
    end

    Knet.array_type[] = save_type
    am
end

function torch2knet(p::PyObject; o...)
    converter = Symbol(pytypeof(p).__name__ * "_pt")
    if isdefined(@__MODULE__, converter)
        eval(converter)(p; o...)
    else
        @warn "$converter not defined"
    end
end

### ResNet

function ResNet_pt(p::PyObject)
    layers = (p.layer1, p.layer2, p.layer3, p.layer4)
    s = Block(ResNetInput_pt(p.conv1, p.bn1); name="ResNet$(length.(layers))")
    for (i, layer) in enumerate(layers)
        blocks = Block(; name="Layer$i")
        for block in layer
            push!(blocks, ResNetBlock_pt(block))
        end
        push!(s, blocks)
    end
    push!(s, ResNetOutput_pt(p.fc))
end

function ResNetBlock_pt(p::PyObject)
    i, layers = 1, []
    while haskey(p, "conv$i")
        convi = getproperty(p, "conv$i")
        bni = getproperty(p, "bn$i")
        push!(layers, Conv2d_pt(convi; normalization=torch2knet(bni), activation=relu))
        i += 1
    end
    layers[end].activation = nothing
    r = (p.downsample === nothing ? identity :
         Conv2d_pt(p.downsample[1]; normalization=torch2knet(p.downsample[2])))
    Add(Block(layers...), r; activation=relu)
end

function ResNetInput_pt(conv1::PyObject, bn1::PyObject)
    Block(
        Conv2d_pt(conv1; normalization=torch2knet(bn1), activation=relu),
        Op(pool; window=3, stride=2, padding=1);
        name = "Input"
    )
end

function ResNetOutput_pt(fc::PyObject)
    w = param(t2a(fc.weight))
    bias = param(t2a(fc.bias))
    Block(
        Op(pool; op=mean, window=typemax(Int)),
        reshape2d,
        Linear(w; bias);
        name = "Output"
    )
end

function Conv2d_pt(p::PyObject; normalization=nothing, activation=nothing)
    w = param(permutedims(t2a(p.weight), (4,3,2,1)))
    bias = (p.bias === nothing ? nothing :
            param(reshape(t2a(p.bias), (1,1,:,1))))
    Conv(w; bias, normalization, activation,
         p.padding, p.stride, p.dilation, p.groups)
end

function BatchNorm2d_pt(b::PyObject)
    bnweight(x) = param(reshape(t2a(x), (1,1,:,1)))
    BatchNorm(
        ; use_estimates = nothing,
        update = b.momentum,
        mean = bnweight(b.running_mean).value,
        var = bnweight(b.running_var).value,
        bias = bnweight(b.bias),
        scale = bnweight(b.weight),
        epsilon = b.eps,
    )
end

function BatchNorm2d_gcheck(b::PyObject) # use this for @gcheck, defaults will be fixed when weights loaded from file
    bnweight(x) = param(reshape(t2a(x), (1,1,:,1)))
    BatchNorm(
        ; use_estimates = true, ### nothing,
        update = 0, ### b.momentum,
        mean = bnweight(b.running_mean), ### .value,
        var = bnweight(b.running_var), ### .value,
        bias = bnweight(b.bias),
        scale = bnweight(b.weight),
        epsilon = b.eps,
    )
end

##### MobileNetV2 ##################################################################

function adaptive_avg_pool2d(x)
    y = pool(x; op=mean, window=typemax(Int))
    reshape2d(y)
end

function MobileNetV2_pt(p::PyObject)
    s = Block()
    for l in p.features
        push!(s, torch2knet(l))
    end
    push!(s, adaptive_avg_pool2d)
    dropout = p.classifier[1].p
    push!(s, torch2knet(p.classifier[2]; dropout))
    return s
end

function MobileNetV3_pt(p::PyObject)
    s = Block()
    for l in p.features
        push!(s, torch2knet(l))
    end
    push!(s, torch2knet(p.avgpool))
    push!(s, torch2knet(p.classifier[1]; activation=torch2knet(p.classifier[2])))
    dropout = p.classifier[3].p
    push!(s, torch2knet(p.classifier[4]; dropout))
    return s
end

function Dropout_pt(p::PyObject)
    Op(dropout, p.p)
end

function Linear_pt(p::PyObject; dropout=0, activation=nothing)
    w = param(t2a(p.weight))
    bias = p.bias === nothing ? nothing : param(t2a(p.bias))
    Linear(w; bias, dropout, activation)
end

function ConvBNReLU_pt(p::PyObject)
    activation = torch2knet(p[3])
    normalization = torch2knet(p[2])
    Conv2d_pt(p[1]; activation, normalization)
end

function ReLU6_pt(p::PyObject)
    Op(relu; max_value=6)
end

function InvertedResidual_pt(p::PyObject)
    haskey(p, :conv) ? InvertedResidualV2(p) : InvertedResidualV3(p)
end

function InvertedResidualV2(p::PyObject)
    s = Block()
    for l in p.conv
        if l == p.conv[end]
            @assert typename(l) == "BatchNorm2d"
            s[end].normalization = torch2knet(l)
        else
            push!(s, torch2knet(l))
        end
    end
    if p.use_res_connect # s[2].stride == (1,1) && size(s[1].w,3) == size(s[end].w,4)
        return Add(s, identity)
    else
        return s
    end
end

function InvertedResidualV3(p::PyObject)
    s = Block()
    for l in p.block
        push!(s, torch2knet(l))
    end
    if p.use_res_connect # s[2].stride == (1,1) && size(s[1].w,3) == size(s[end].w,4)
        return Add(s, identity)
    else
        return s
    end
end

function SqueezeExcitation_pt(p::PyObject)
    Mul(Block(
        torch2knet(p.fc1; activation=relu),
        torch2knet(p.fc2; activation=hardsigmoid),
    ), identity)
end

function ConvBNActivation_pt(p::PyObject)
    activation = torch2knet(p[3])
    normalization = torch2knet(p[2])
    torch2knet(p[1]; activation, normalization)
end

function AdaptiveAvgPool2d_pt(p::PyObject)
    return adaptive_avg_pool2d
end

function Hardswish_pt(p::PyObject)
    return hardswish
end

function ReLU_pt(p::PyObject)
    return relu
end

function Identity_pt(p::PyObject)
    return identity
end


text:
# # Instructions for Assignment02
#
# ## Introduction
#
# In this assignment,
# you will start to put the pieces together.
# Note: this file is runnable in its current state,
# but is incomplete.
# You can run the file from the command line in script mode,
# copy code into the REPL in interactive mode,
# or use the VS Code julia extension to run individual lines.

# The following code is used for set up;
# there is no need to change it.
# If you _do_ change it to figure out what it's doing,
# make sure you change it back or your assignment may not work correctly

## To view additional instructions for this assignment, #jl
## see https://wellesley-bisc195.github.io/assignments/Assignment02 #jl

## To view a rendered version of this document, #jl
## see https://wellesley-bisc195.github.io/assignments/Assignment02/#assignment02_code #jl

using Random
Random.seed!(42)

function generate_sequence(len)
    seq = join(rand(['A','C','G','T'], len))
    println("Your sequence is: ", seq)
    return seq
end

my_seq = generate_sequence(20)

# ## Question 1
#
# The code above generates a random 20nt DNA sequence,
# and assigns it to the variable `my_seq`.
# What kind of data type is `my_seq` (eg `Float64`, `Int64`, `String`, or something else)?
# Assign the variable `question1` to the correct type.
#
# Hint1: rather than guess, you can just use the `typeof()` function.

question1 = ""

# Hint2: If you're having trouble,
# note that your answer should not contain quotes.
# Eg. if you think that `seq` has the type `MyType`,
# your code should read `question1 = MyType` and not `question1 = "MyType"`

# ## Question 2
#
# The function below has a documentation string or "doc string",
# which explains what the function should do.
# Many julia functions have doc strings,
# which makes it easy to get help right from the REPL
#
# Open a julia REPL and type `?`.
# You will see the `julia>` prompt change to `help?>`.
# Now type `println` and hit enter.
# The docstring of `println` is printed
# and you're returned to a julia prompt
#
# Evaluate the following function, including the docstring, in the REPL
# (that is, copy from the triple quotes down to the "end"
# and paste it into the REPL, then hit enter)

"""
    question2(sequence)

Checks if `sequence` is a String.

Example
≡≡≡≡≡≡≡≡≡

    julia> question2("hello")
    true

    julia> question2(1001)
    false

    julia> if question2("I'm a string!")
               println("yeah, that was a string")
           end
    yeah, that was a string
"""
function question2()
    ## put your code here
end

# At the moment, this function doesn't do what it says it should.
# Actually, it doesn't do anything.
# Fix the function so that it takes one argument
# and returns `true` if the argument is a String,
# and returns `false` otherwise.
# You should be able to run the examples in the docstring
# and get the correct answer,
# and to run this function using `my_seq` as the argument
# (it should return `true`).
#
# Hint: You haven't explicitly encountered a way to check if a type is another type.
# Learning to search for answers is a key programming skill!
# Try searching "check if type is string julia".
# Typically, when I see search results, I prioritize in this order:
# 1. stackoverflow.com
# 2. docs.julialang.org
# 3. discourse.julialang.org
# 4. github.com

# ## Question 3
#
# Let's try something a bit more complicated.
# The `question3()` function has a doc string,
# but is also incomplete.
# I've writtten some code to get you started,
# complete it so that the examples work as expected in the docstring.

"""
    question3(sequence)

Calculates the GC ratio of a DNA sequence.
The GC ratio is the total number of G and C bases divided by the total length of the sequence.
For more info about GC content, see here: https://en.wikipedia.org/wiki/GC-content

Example
≡≡≡≡≡≡≡≡≡≡

    julia> question3("AATG")
    0.25

    julia> question3("CCCGG")
    1.0

    julia> question3("ATTA")
    0.0
"""
function question3(sequence)
    ## throw an error if the string contains anything other than ACGT
    if any(c-> !in(c, ['A','C','G','T']), sequence)
        throw(ArgumentError("Sequence must only contain ACGT"))
    end

    ## change line to assign `seqlength` to the length of `sequence` instead of `1`
    ## If you're stuck, search for "length of string julia"
    seqlength = 1

    ## count the number of G's
    gs = count(==('G'), sequence)
    ## count the number of C's
    cs = count(==('C'), sequence)

    return gs + cs / seqlength # something is wrong with this line...
end

# ## Question 4
#
# This process should be familiar by now.
#
# Hint: you do not need to re-write code if you already have it.
# That is, don't copy any code from `question3`,
# just call it!

"""
    question4(sequence)

Calculates the GC content of a DNA sequence
and prints it to the screen.

Example
≡≡≡≡≡≡≡≡≡≡

    julia> question4("AATC")
    Sequence:
    AATC
    GC Content:
    0.25

    julia> question4("CCCGG")
    Sequence:
    CCCGG
    GC Content:
    1.0

    julia> question4("CCCGT"); # note the semicolon
    Sequence:
    CCCGG
    GC Content:
    0.8
"""
function question4(sequence)
    ## Your code here
end

## Question 5
#
# The following are related to questions in the exercises
# from chapter 2 of _Think Julia_.
# The format of the varables below refer to specific exercise questions.
# For example, `ce_2_1_1` would refer to chapter excercise 2-2, question 1.
#
# > We’ve seen that `n = 42` is legal. What about `42 = n`?
#
# Assign the variable to `true` if legal, `false` if illegal

ce_2_2_1 = Bool

# > How about `x = y = 1`?
#
# Assign the variable to `true` if legal, `false` if illegal

ce_2_2_2 = Bool

# > The volume of a sphere with radius `r` is `4/3πr^3`.
# What is the volume of a sphere with radius `5`?
#
# Assign the variable to the correct answer

ce_2_3_1 = Float64

# > Suppose the cover price of a book is \$24.95,
# but bookstores get a 40 % discount.
# Shipping costs \$3 for the first copy and 75 cents for each additional copy.
# What is the total wholesale cost for 60 copies?
#
# Assign the variable to the correct answer.

ce_2_3_2 = Float64

## Question 6
#
# Write a function that calculates the cost of books,
# given some list price, discount rate, and number.

"""
    bookprice(list, discount, count)

Calculates the price of books, given some list price, discount rate,
and number of books ordered.

Example
≡≡≡≡≡≡≡≡≡≡

    julia> bookprice(24.95, 0.4, 60)
    # this should be the same answer you got for `c_2_3_2`

    julia> bookprice(24.95, 0.4, 1)
    17.97

    julia> bookprice(24.95, 0.4, 2)
    33.69

    julia> bookprice(1, 0, 1)
    4.0
"""
function bookprice(list, discount, count)
    ## your code here
end


text:
# --------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
function get_env(t::TF, pulse::Pulse{TE, TF}) where {TE, TF}
	@unpack tL, tR, T, env_type, ncyc = pulse
	if t < tL || tR < t
		return 0
	end
	if env_type == "sin2"
		return sin((π * t) / (ncyc * T))^2
	elseif env_type == "sin4"
		return sin((π * t) / (ncyc * T))^4
	else env_type == "flat" || env_type == "delta"
		return 1.0
	end
end
function get_env(t::TVF, pulse::Pulse) where {TVF <: AbstractVector{<:AbstractFloat}}
	env = zeros(eltype(t), size(t))
	for it in 1:length(t)
		env[it] = get_env(t[it], pulse)
	end
	return env
end
## electric field
function getEt(t::TF, pulse::Pulse{Et, TF}) where TF
	@unpack E, ω, ϕ = pulse
	return E * get_env(t, pulse) * sin(ω * t + ϕ)
end
function getEt(t::TVF, pulse::Pulse{Et}) where {TVF <: AbstractVector{<:AbstractFloat}}
	Et = zeros(eltype(t), length(t))
	for it in 1:length(t)
		Et[it] = getEt(t[it], pulse)
	end
	return Et
end
function getEt(t::TF, laser::Vector{Pulse{Et}}) where TF <: AbstractFloat
	Et = zero(t)
	for p in laser
		Et += getEt(t, p)
	end
	return Et
end
function getEt(t::TVF, laser::Vector{Pulse{Et}})  where {TVF <: AbstractVector{<:AbstractFloat}}
	Et = zeros(eltype(t), length(t))
	for it in 1:length(t)
		Et[it] = getEt(t[it], laser)
	end
	return Et
end
## vector potential
function getAt(t::TF, pulse::Pulse{Et, TF},
		tprev::TF=0., Atprev::TF=0.) where TF
	@unpack tL, tR, T = pulse
	if t < tL || tR < t
		return 0
	end
	nt = 1
	while (t - tprev) / nt > T / 100000
		nt *= 10
	end
	dt = (t - tprev) / nt
	At = Atprev
	for it in 1:nt
		t0 = tprev + (it - 1) * dt
		t1 = t0 + dt / 2
		t2 = t0 + dt
		At += (-dt / 6) * (getEt(t0, pulse) + 4getEt(t1, pulse) + getEt(t2, pulse))
	end
	return At
end
function getAt(t::TVF, pulse::Pulse{Et}) where {TVF <: AbstractVector{<:AbstractFloat}}
	At = zeros(eltype(t), size(t))
	At[begin] = getAt(t[begin], pulse)
	for it in 2:length(t)
		At[it] = getAt(t[it], pulse, t[it - 1], At[it - 1])
	end
	return At
end
function getAt(t::TF, laser::Vector{Pulse{Et}},
		tprev::TVF=zeros(eltype(t), size(laser)),
		Atprev::TVF=zeros(eltype(t), size(laser))) where {TF <: AbstractFloat,TVF <: AbstractVector{<:AbstractFloat}}
	At = zero(t)
	for ip in 1:length(laser)
		At += getAt(t, laser[ip], tprev[ip], Atprev[ip])
	end
	return At
end
function getAt(t::TVF, laser::Vector{Pulse{Et}}) where {TVF <: AbstractVector{<:AbstractFloat}}
	At = zeros(eltype(t), size(t))
	for p in laser
		Attmp = getAt(t, p)
		@. At += Attmp
	end
	return At
end


text:
using CPUOptimizedConvNet
using Documenter

makedocs(;
    modules=[CPUOptimizedConvNet],
    authors="Chris Elrod <elrodc@gmail.com> and contributors",
    repo="https://github.com/chriselrod/CPUOptimizedConvNet.jl/blob/{commit}{path}#L{line}",
    sitename="CPUOptimizedConvNet.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)


text:
using Distributed, Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
if haskey(ENV, "FOOD_CACHING_PARALLEL")
    nprocs() < Sys.CPU_THREADS && addprocs(Sys.CPU_THREADS - nprocs(),
                                           exeflags = "--project",
                                           dir = joinpath(@__DIR__, ".."))
end
using FoodCachingModels, DataFrames, Distances, Statistics
import FoodCachingExperiments: bsave, target
import FoodCachingFitting: run_experiments, ALLEXPERIMENTS, DATADIR, _logp_hat, simname

function bootstrap(f, x; N = 100)
    result = [f(x[rand(1:length(x), length(x))]) for _ in 1:N]
    std(result)
end
results = DataFrame(model = [], experiment = [], logp_hat = [],
                    logp_hat_std = [], id = [],
                    avg = [], best = [], best_seed = [])
rev = "b15ffc2"
for model in [Baseline, MotivationalControl, EpisodicLikeMemory,
              PlasticCaching, ReplayAndPlan]
    for e in ALLEXPERIMENTS
        for id in 1:4
            println("Running $model, $e, $id.")
            tmp = run_experiments(simname(model, [e], id, rev), [e], 20000)
            metric = Distances.Euclidean()
            x = target(e)
            d = metric.(Ref(x), tmp.results)
            s = sortperm(d)
            d = d[s]
            f = d -> _logp_hat(d, x, 5, metric)
            l = f(d)
            st = bootstrap(f, d)
            push!(results, [Symbol(model), e, l, st, id, mean(tmp.results),
                            tmp.results[s[1]], tmp.seed[s[1]]])
        end
    end
end
bsave(joinpath(DATADIR, "run_indi_$rev"), Dict(:results => results))


text:
module TestModels

using Test
using TMLEEpistasis
using MLJ

@testset "Test InteractionTransformer" begin
    X = (rs1234=[1, 2, 3], rs455=[4, 5, 6], rs4489=[7, 8, 9], rstoto=[1, 2, 3])
    t = TMLEEpistasis.InteractionTransformer(r"^rs[0-9]+")
    mach = machine(t, X)
    fit!(mach)
    Xt = MLJ.transform(mach, X)

    @test Xt == (
        rs1234 = [1, 2, 3],
        rs455 = [4, 5, 6],
        rs4489 = [7, 8, 9],
        rstoto = [1, 2, 3],
        rs1234_rs455 = [4.0, 10.0, 18.0],
        rs1234_rs4489 = [7.0, 16.0, 27.0],
        rs455_rs4489 = [28.0, 40.0, 54.0]
    )
    @test mach.fitresult.ninter == 3
    @test mach.fitresult.interaction_pairs == [:rs1234 => :rs455, :rs1234 => :rs4489, :rs455 => :rs4489]

end

@testset "Test InteractionLM" begin
    n = 100
    X = (rs1234=rand(n), rs455=rand(n), rstoto=rand(n))
    # Classifier
    y = categorical(rand([0,1], n))
    mach = machine(
        TMLEEpistasis.InteractionLMClassifier(),
        X,
        y
    )
    fit!(mach)
    fp = fitted_params(mach)
    @test fp.interaction_transformer.fitresult.interaction_pairs == [:rs1234 => :rs455]
    @test predict(mach) isa MLJ.UnivariateFiniteVector{Multiclass{2}, Int64, UInt32, Float64}
    # Regressor
    y = rand(n)
    mach = machine(
        TMLEEpistasis.InteractionLMRegressor(),
        X,
        y
    )
    fit!(mach)
    fp = fitted_params(mach)
    @test fp.interaction_transformer.fitresult.interaction_pairs == [:rs1234 => :rs455]
    @test predict(mach) isa Vector{Float64}

    @test target_scitype(TMLEEpistasis.InteractionLMRegressor()) == AbstractVector{<:MLJ.Continuous}
    @test target_scitype(TMLEEpistasis.InteractionLMClassifier()) == AbstractVector{<:MLJ.Finite}

end

end
true


text:
# ---
# title: 779. K-th Symbol in Grammar
# id: problem779
# author: Indigo
# date: 2021-01-28
# difficulty: Medium
# categories: Recursion
# link: <https://leetcode.com/problems/k-th-symbol-in-grammar/description/>
# hidden: true
# ---
#
# On the first row, we write a `0`. Now in every subsequent row, we look at the
# previous row and replace each occurrence of `0` with `01`, and each occurrence
# of `1` with `10`.
#
# Given row `N` and index `K`, return the `K`-th indexed symbol in row `N`. (The
# values of `K` are 1-indexed.) (1 indexed).
#
#
#
#     **Examples:**
#     Input: N = 1, K = 1
#     Output: 0
#
#     Input: N = 2, K = 1
#     Output: 0
#
#     Input: N = 2, K = 2
#     Output: 1
#
#     Input: N = 4, K = 5
#     Output: 1
#
#     Explanation:
#     row 1: 0
#     row 2: 01
#     row 3: 0110
#     row 4: 01101001
#
#
# **Note:**
#
#   1. `N` will be an integer in the range `[1, 30]`.
#   2. `K` will be an integer in the range `[1, 2^(N-1)]`.
#
#
## @lc code=start
using LeetCode

function kth_grammar(N::Int, K::Int)
    res = false
    for i in 1:N
        iseven(K) && (res = !res)
        K = (K + 1) >> 1
    end
    res
end
## @lc code=end


text:
# warn("Must remove optional dev/ISAMRemoteSolve.jl dependency here.")
# include("dev/ISAMRemoteSolve.jl")

immutable RangeAzimuthElevation
  range::Float64
  azimuth::Float64
  elevation::Union{Void,Float64}
end

function convert{T <: CoordinateTransformations.AffineMap}(::Type{T}, x::SE3)
  q = convert(TransformUtils.Quaternion, x.R)
  Translation(x.t...) ∘ LinearMap( Rotations.Quat(q.s, q.v...) )
end
function convert{T <: CoordinateTransformations.AffineMap{Rotations.Quat{Float64}}}(::Type{SE3}, x::T)
  SE3(x.v[1:3], TransformUtils.Quaternion(x.m.w, [x.m.x,x.m.y,x.m.z]) )
end

function convert(::Type{SE3}, t::Tuple{Symbol, Vector{Float64}})
  if t[1]==:XYZqWXYZ
    return SE3(t[2][1:3],TransformUtils.Quaternion(t[2][4],t[2][5:7]))
  else
    error("Unknown conversion type $(t[1])")
  end
end

function convert(::Type{RangeAzimuthElevation}, val::Tuple{Symbol, Vector{Float64}})
  if val[1] == :rangeazimuth
    return RangeAzimuthElevation(val[2][1],val[2][2],nothing)
  elseif val[1] == :rangeazimuthelevation
    return RangeAzimuthElevation(val[2][1],val[2][2],val[2][3])
  else
    error("Unknown conversion from $(val[1]) to RangeAzimuthElevation.")
  end
end

# should be deprecated or indicated more clearly
lsrBR(a) = [a[2,:];a[1,:]]';

function veePose3(s::SE3)
  TransformUtils.veeEuler(s)
end
function veePose(s::SE3)
  TransformUtils.veeEuler(s)
end

function \(s::SE3, wTr::CTs.Translation)
  bTr = s.R.R'*(wTr.v-s.t)
  Dtr = bTr
  range = norm(Dtr)
  azi = atan2(Dtr[2], Dtr[1])
  elev = atan2(Dtr[3], Dtr[1])
  RangeAzimuthElevation(range, azi, elev)
end

"""
    getRangeKDEMax2D(fgl::FactorGraph, vsym1::Symbol, vsym2::Symbol)

Calculate the cartesian distance between two vertices in the graph using their symbol name, and by maximum belief point.
"""
function getRangeKDEMax2D(fgl::FactorGraph, vsym1::Symbol, vsym2::Symbol)
  x1 = getKDEMax(getVertKDE(fgl, vsym1))
  x2 = getKDEMax(getVertKDE(fgl, vsym2))
  norm(x1[1:2]-x2[1:2])
end

function measureMeanDist(fg::FactorGraph, a::T, b::T) where {T <: AbstractString}
    #bearrang!(residual::Array{Float64,1}, Z::Array{Float64,1}, X::Array{Float64,1}, L::Array{Float64,1})
    res = zeros(2)
    A = getVal(fg,a)
    B = getVal(fg,b)
    Ax = Base.mean(vec(A[1,:]))
    Ay = Base.mean(vec(A[2,:]))
    Bx = Base.mean(vec(B[1,:]))
    By = Base.mean(vec(B[2,:]))
    dx = Bx - Ax
    dy = By - Ay
    b = atan2(dy,dx)
    r = sqrt(dx^2 + dy^2)
    return r, b
end

function predictBodyBR(fg::FactorGraph, a::T, b::T) where {T <: AbstractString}
  res = zeros(2)
  A = getVal(fg,a)
  B = getVal(fg,b)
  Ax = Base.mean(vec(A[1,:]))
  Ay = Base.mean(vec(A[2,:]))
  Ath = Base.mean(vec(A[3,:]))
  Bx = Base.mean(vec(B[1,:]))
  By = Base.mean(vec(B[2,:]))
  wL = SE2([Bx;By;0.0])
  wBb = SE2([Ax;Ay;Ath])
  bL = se2vee((wBb \ eye(3)) * wL)
  dx = bL[1] - 0.0
  dy = bL[2] - 0.0
  b = (atan2(dy,dx))
  r = sqrt(dx^2 + dy^2)
  return b, r
end

function getNextLbl(fgl::FactorGraph, chr)
  # TODO convert this to use a double lookup
  warn("getNextLbl(::FactorGraph..) to be deprecated, use getlastpose/landm(::SLAMWrapper..) instead.")
  max = -1
  maxid = -1
  for vid in fgl.IDs
  # for v in fgl.v #fgl.g.vertices # fgl.v
      v = (vid[2], fgl.g.vertices[vid[2]])
      if v[2].attributes["label"][1] == chr
        val = parse(Int,v[2].attributes["label"][2:end])
        if max < val
          max = val
          maxid = v[1]
        end
      end
  end
  if maxid != -1
    v = getVert(fgl,maxid)
    X = getVal(v)
    return v, X, Symbol(string(chr,max+1))
  else
    return nothing, nothing, Symbol(string(chr,max+1)) # Union{}
  end
end

function getLastPose(fgl::FactorGraph)
  return getNextLbl(fgl, 'x')
end
getLastPose2D(fgl::FactorGraph) = getLastPose(fgl)

function getlastpose(slam::SLAMWrapper)
  error("getlastpose -- Not implemented yet")
end

function getLastLandm2D(fgl::FactorGraph)
  return getNextLbl(fgl, 'l')
end

function odomKDE(p1,dx,cov)
  X = getPoints(p1)
  sig = diag(cov)
  RES = zeros(size(X))
  # increases the number of particles based on the number of modes in the measurement Z
  for i in 1:size(X,2)
      ent = [randn()*sig[1]; randn()*sig[2]; randn()*sig[3]]
      RES[:,i] = addPose2Pose2(X[:,i], dx + ent)
  end
  return kde!(RES, "lcv")
end

"""
    addOdoFG!(fg, name, DX, cov, , N=100, ready=1)

Create a new variable node and insert odometry constraint factor between
which will automatically increment latest pose symbol x<k+1> for new node new node and
constraint factor are returned as a tuple.

"""
function addOdoFG!(
        fg::FactorGraph,
        n::Symbol,
        DX::Array{Float64,1},
        cov::Array{Float64,2};
        N::Int=0,
        ready::Int=1,
        labels::Vector{<:AbstractString}=String[]  )
    #
    prev, X, nextn = getLastPose2D(fg)
    r,c = size(X)
    if N==0
      N = c
    end
    sig = diag(cov)
    XnextInit = zeros(r,c)
    # increases the number of particles based on the number of modes in the measurement Z
    for i in 1:c
        ent = [randn()*sig[1]; randn()*sig[2]; randn()*sig[3]]
        XnextInit[:,i] = addPose2Pose2(X[:,i], DX + ent)
    end

    v = addNode!(fg, n, XnextInit, cov, N=N, ready=ready, labels=labels)
    pp = Pose2Pose2(vectoarr2(DX), cov, [1.0]) #[prev;v],
    f = addFactor!(fg, [prev;v], pp, ready=ready)
    infor = inv(cov^2)
    # addOdoRemote(prev.index,v.index,DX,infor) # this is for remote factor graph ref parametric solution -- skipped internally by global flag variable
    return v, f
end

function addOdoFG!(
        fg::FactorGraph,
        n::T,
        Z::Pose3Pose3;
        N::Int=0,
        ready::Int=1,
        labels::Vector{T}=String[]  ) where {T <: AbstractString}
  #
  DX=Z.μ
  cov=Z.Σ.mat

  addOdoFG(fg, n, DX, cov, N=N, ready=ready, labels=labels)
end

"""
    addOdoFG!(fg, odo, N=100, ready=1)

Create a new variable node and insert odometry constraint factor between
which will automatically increment latest pose symbol x<k+1> for new node new node and
constraint factor are returned as a tuple.

"""
function addOdoFG!(
        fgl::FactorGraph,
        odo::PP;
        N::Int=0,
        ready::Int=1,
        labels::Vector{T}=String[] ) where {PP <: RoME.BetweenPoses, T <: AbstractString}
    #
    vprev, X, nextn = getLastPose(fgl)
    if N==0
      N = size(X,2)
    end
    vnext = addNode!(fgl, nextn, X⊕odo, ones(1,1), N=N, ready=ready, labels=labels)
    fact = addFactor!(fgl, [vprev;vnext], odo)
    return vnext, fact
end

function initfg(;sessionname="NA")
  fgl = emptyFactorGraph()
  fgl.sessionname=sessionname
  # registerCallback!(fgl, RoME.getSample) # RoME.evalPotention
  return fgl
end

function initFactorGraph!(fg::FactorGraph;
      P0::Union{Array{Float64,2},Void}=nothing,
      init::Union{Vector{Float64},Void}=nothing,
      N::Int=100,
      lbl::Symbol=:x0,
      ready::Int=1,
      firstPoseType=Pose2,
      labels::Vector{<:AbstractString}=String[])
  #
  nodesymbols = Symbol[]
  if firstPoseType == Pose2
      init = init!=nothing ? init : zeros(3)
      P0 = P0!=nothing ? P0 : diagm([0.03;0.03;0.001])
      init = vectoarr2(init)
      addNode!(fg,lbl,Pose2,N=N,autoinit=true,ready=ready,labels=String["VARIABLE"; labels] )
      push!(nodesymbols, lbl)
      # v1 = addNode!(fg, lbl, init, P0, N=N, ready=ready, labels=labels)
      fctVert = addFactor!(fg, [lbl;], PriorPose2(init, P0,  [1.0]), ready=ready, labels=String["FACTOR"; labels]) #[v1],
      push!(nodesymbols, fctVert.label)
  end
  if firstPoseType == Pose3
      init = init!=nothing ? init : zeros(6)
      P0 = P0!=nothing ? P0 : diagm([0.03;0.03;0.03;0.001;0.001;0.001])
      addNode!(fg,lbl,Pose2,N=N,autoinit=true,ready=ready,labels=String["VARIABLE"; labels] )
      push!(nodesymbols, lbl)
      # v1 = addNode!(fg, lbl, init, P0, N=N, ready=ready, labels=labels)
      fctVert = addFactor!(fg, [lbl;], PriorPose3(MvNormal(init, P0)), ready=ready, labels=String["FACTOR"; labels]) #[v1],
      push!(nodesymbols, fctVert.label)
  end
  warn("initFactorGraph! -- returns variable and factor symbols ")
  return nodesymbols
end

function newLandm!(fg::FactorGraph, lm::T, wPos::Array{Float64,2}, sig::Array{Float64,2};
                  N::Int=100, ready::Int=1, labels::Vector{T}=String[]) where {T <: AbstractString}

    # TODO -- need to confirm this function is updating the correct memory location. v should be pointing into graph
    vert=addNode!(fg, Symbol(lm), wPos, sig, N=N, ready=ready, labels=labels)

    vert.attributes["age"] = 0
    vert.attributes["maxage"] = 0
    vert.attributes["numposes"] = 0
    updateFullVert!(fg, vert)

    println("newLandm! -- added $(lm)")
    return vert
end

function updateLandmAge(vlm::Graphs.ExVertex, pose::T) where {T <: AbstractString}
  error("still working here")
end

function addBRFG!(fg::FactorGraph,
      pose::T,
      lm::T,
      br::Array{Float64,1},
      cov::Array{Float64,2};
      ready::Int=1  ) where {T <: AbstractString}
  #
  vps = getVert(fg,pose)
  vlm = getVert(fg,lm)
  testlbl = vps.label*vlm.label
  for nei in getOutNeighbors(fg, vlm)
    if nei.label == testlbl
      # TODO -- makes function call brittle
      warn("We already have $(testlbl), skipping this constraint")
      return nothing
    end
  end
  @show keys(vlm.attributes)
  np = vlm.attributes["numposes"]
  la = vlm.attributes["age"]
  nage = parse(Int,pose[2:end])
  vlm.attributes["numposes"] = np+1
  vlm.attributes["age"] = ((la*np)+nage)/(np+1)
  vlm.attributes["maxage"] = nage
  updateFullVert!(fg, vlm)

  # @show br, cov
  pbr = Pose2DPoint2DBearingRange{Normal, Normal}(Normal(br[1], cov[1,1]), Normal(br[2],  cov[2,2]))
  # pbr = Pose2DPoint2DBearingRange((br')',  cov,  [1.0])
  f = addFactor!(fg, [vps;vlm], pbr, ready=ready ) #[vps;vlm],

  # only used for max likelihood unimodal tests.
  u, P = pol2cart(br[[2;1]], diag(cov))
  infor = inv(P^2)
  # addLandmMeasRemote(vps.index,vlm.index,u,infor) # for iSAM1 remote solution as reference
  return f
end

function addMMBRFG!(fg::FactorGraph, pose::T,
                  lm::Array{T,1}, br::Array{Float64,1},
                  cov::Array{Float64,2}; w=[0.5;0.5], ready::Int=1) where {T <: AbstractString}
    #
    vps = getVert(fg,pose)
    vlm1 = getVert(fg,lm[1])
    vlm2 = getVert(fg,lm[2])

    pbr = Pose2DPoint2DBearingRange(vectoarr2(br),  cov,  w) #[vps;vlm1;vlm2],
    f = addFactor!(fg, [vps;vlm1;vlm2], pbr, ready=ready )
    return f
end

function projNewLandmPoints(vps::Graphs.ExVertex, br::Array{Float64,1}, cov::Array{Float64,2})
    # TODO -- convert to use Distributions and common projection function
    Xps = getVal(vps)
    lmPts = zeros(2,size(Xps,2))
    for i in 1:size(Xps,2)
        ent = [cov[1,1]*randn(); cov[2,2]*randn()]
        init = vec(Xps[1:2,i])+randn(2)
        lmPts[:,i] = solveLandm(br + ent, vec(Xps[:,i]), init)
    end
    return lmPts
end

function projNewLandm!(fg::FactorGraph, pose::T, lm::T, br::Array{Float64,1}, cov::Array{Float64,2};
                        addfactor=true, N::Int=100, ready::Int=1,
                        labels::Vector{T}=String[]) where {T <: AbstractString}
    #
    vps = getVert(fg,pose)

    lmPts = projNewLandmPoints(vps, br, cov)
    vlm = newLandm!(fg, lm, lmPts, cov, N=N, ready=ready, labels=labels) # cov should not be required here
    if addfactor
      fbr = addBRFG!(fg, pose, lm, br, cov, ready=ready)
      return vlm, fbr
    end
    return vlm
end

function calcIntersectVols(fgl::FactorGraph, predLm::BallTreeDensity;
                          currage=0, maxdeltaage=Inf)
    # all landmarks of interest
    xx,ll = ls(fgl)
    # output result
    rr = Dict{String, RemoteRef}()
    fetchlist = String[]
    iv = Dict{String, Float64}()
    for l in ll
      pvlm = getVert(fgl,l)
      # TODO -- can be improved via query in DB case
      if currage - pvlm.attributes["maxage"] < maxdeltaage
        p = getVertKDE(fgl, l)
        rr[l] = remotecall(uppA(), intersIntgAppxIS, p,predLm)
        push!(fetchlist, l)
      else
        println("calcIntersectVols -- ignoring $(l) because maxdeltaage exceeded")
        iv[l] = 0
      end
    end
    max = 0
    maxl = String("")
    for l in fetchlist #ll
      # p = getVertKDE(fgl, l)
      # tv = intersIntgAppxIS(p,predLm)
      # iv[l] = tv
      tv = fetch(rr[l])
      iv[l] = tv
      if max < tv  max = tv; maxl = l; end
    end
    return iv, maxl
end

function maxIvWithoutID(ivs::Dict{String, Float64}, l::T) where {T <: AbstractString}
  max = 0
  maxl = String("")
  for i in ivs
    if max < i[2] && i[1] != l;  max = i[2]; maxl = i[1]; end
  end
  return maxl
end

# binary tests to distinguish how to automatically add a landmark to the existing factor graph
function doAutoEvalTests(fgl::FactorGraph, ivs::Dict{T, Float64}, maxl::T, lmid::Int, lmindx::Int) where {T <: AbstractString}
  maxAnyval = maxl != String("") ? ivs[maxl] : 0.0
  # maxid = fgl.IDs[maxl]
  lmidSugg = lmid != -1 # a landmark ID has been suggested
  maxAnyExists = maxAnyval > 0.03 # there is notable intersection with previous landm
  lmIDExists = haskey(fgl.v, lmid) # suggested lmid already in fgl
  newlmindx = lmindx
  if lmIDExists
    lmSuggLbl = String(getVert(fgl,lmid).label) # TODO -- wasteful
  else
    newlmindx = lmindx + 1
    lmSuggLbl = String(string('l',newlmindx))
  end
  maxl2 = lmIDExists ? maxIvWithoutID(ivs, lmSuggLbl) : String("")
  maxl2Exists = lmIDExists ? (maxl2 != "" ? ivs[maxl2] > 0.03 : false) : false # there is notable intersection with previous landm
  intgLmIDExists = lmIDExists ? ivs[lmSuggLbl] > 0.03 : false

  return lmidSugg, maxAnyExists, maxl2Exists, maxl2, lmIDExists, intgLmIDExists, lmSuggLbl, newlmindx
end

function evalAutoCases!(fgl::FactorGraph, lmid::Int, ivs::Dict{T, Float64}, maxl::T,
                        pose::T, lmPts::Array{Float64,2}, br::Array{Float64,1}, cov::Array{Float64,2}, lmindx::Int;
                        N::Int=100, ready::Int=1 ) where {T <: AbstractString}
  lmidSugg, maxAnyExists, maxl2Exists, maxl2, lmIDExists, intgLmIDExists, lmSuggLbl, newlmindx = doAutoEvalTests(fgl,ivs,maxl,lmid, lmindx)

  println("evalAutoCases -- found=$(lmidSugg), $(maxAnyExists), $(maxl2Exists), $(lmIDExists), $(intgLmIDExists)")

  vlm = Union{}; fbr = Union{};
  if (!lmidSugg && !maxAnyExists)
    #new landmark and UniBR constraint
    v,L,lm = getLastLandm2D(fgl)
    vlm = newLandm!(fgl, lm, lmPts, cov, N=N,ready=ready)
    fbr = addBRFG!(fgl, pose, lm, br, cov, ready=ready)
  elseif !lmidSugg && maxAnyExists
    # add UniBR to best match maxl
    vlm = getVert(fgl,maxl)
    fbr = addBRFG!(fgl, pose, maxl, br, cov, ready=ready)
  elseif lmidSugg && !maxl2Exists && !lmIDExists
    #add new landmark and add UniBR to suggested lmid
    vlm = newLandm!(fgl, lmSuggLbl, lmPts, cov, N=N, ready=ready)
    fbr = addBRFG!(fgl, pose, lmSuggLbl, br, cov, ready=ready)
  elseif lmidSugg && !maxl2Exists && lmIDExists && intgLmIDExists
    # doesn't self intesect with existing lmid, add UniBR to lmid
    vlm = getVert(fgl, lmid)
    fbr = addBRFG!(fgl, pose, lmSuggLbl, br, cov, ready=ready)
  elseif lmidSugg && maxl2Exists && !lmIDExists
    # add new landmark and add MMBR to both maxl and lmid
    vlm = newLandm!(fgl, lmSuggLbl, lmPts, cov, N=N, ready=ready)
    addMMBRFG!(fgl, pose, [maxl2;lmSuggLbl], br, cov, ready=ready)
  elseif lmidSugg && maxl2Exists && lmIDExists && intgLmIDExists
    # obvious case, add MMBR to both maxl and lmid. Double intersect might be the same thing
    println("evalAutoCases! -- obvious case is happening")
    addMMBRFG!(fgl, pose, [maxl2;lmSuggLbl], br, cov, ready=ready)
    vlm = getVert(fgl,lmSuggLbl)
  elseif lmidSugg && maxl2Exists && lmIDExists && !intgLmIDExists
    # odd case, does not intersect with suggestion, but does with some previous landm
    # add MMBR
    warn("evalAutoCases! -- no self intersect with suggested $(lmSuggLbl) detected")
    addMMBRFG!(fgl, pose, [maxl;lmSuggLbl], br, cov, ready=ready)
    vlm = getVert(fgl,lmSuggLbl)
  elseif lmidSugg && !maxl2Exists && lmIDExists && !intgLmIDExists
  #   # landm exists but no intersection with existing or suggested lmid
  #   # may suggest some error
    warn("evalAutoCases! -- no intersect with suggested $(lmSuggLbl) or map detected, adding  new landmark MM constraint incase")
    v,L,lm = getLastLandm2D(fgl)
    vlm = newLandm!(fgl, lm, lmPts, cov, N=N, ready=ready)
    addMMBRFG!(fgl, pose, [lm; lmSuggLbl], br, cov, ready=ready)
  else
    error("evalAutoCases! -- unknown case encountered, can reduce to this error to a warning and ignore user request")
  end

  return vlm, fbr, newlmindx
end

function addAutoLandmBR!(fgl::FactorGraph, pose::T, lmid::Int, br::Array{Float64,1}, cov::Array{Float64,2}, lmindx::Int;
                      N::Int=100, ready::Int=1 ) where {T <: AbstractString}
    vps = getVert(fgl, pose)
    lmPts = projNewLandmPoints(vps, br, cov)
    lmkde = kde!(lmPts)
    currage = parse(Int, pose[2:end])
    ivs, maxl = calcIntersectVols(fgl, lmkde, currage=currage,maxdeltaage=10)

    # There are 8 cases of interest
    vlm, fbr, newlmindx = evalAutoCases!(fgl, lmid, ivs, maxl,pose,lmPts, br,cov,lmindx,N=N,ready=ready)

    return vlm, fbr, newlmindx
end

function malahanobisBR(measA, preA, cov::Array{Float64,2})
    # measure landmark with noise
    res = measA - preA
    mala2 = Union{}
    #Malahanobis distance
    if false
      lambda = cov \ eye(2)
      mala2 = res' * lambda * res
    else
      mala2 = res' * (cov \ res)
    end

    mala = sqrt(mala2)
    return mala
end

# ------------------------------------
# Transfered from IncrementalInference

function get2DSamples(fg::FactorGraph,
      sym;
      from::Int=0, to::Int=999999999,
      minnei::Int=0,
      api::DataLayerAPI=IncrementalInference.localapi  )
  #
  X = Array{Float64,1}()
  Y = Array{Float64,1}()

  # if sym = 'l', ignore single measurement landmarks

  for id in fg.IDs
    vertlbl = string(id[1])
    if vertlbl[1] == sym
      val = parse(Int,vertlbl[2:end])
      if from <= val && val <= to
        if length( getOutNeighbors(fg, id[2], api=api, needdata=true ) ) >= minnei
          # if length(out_neighbors(fg.v[id[2]],fg.g)) >= minnei
          X=[X; vec(getVal(fg,id[2], api=api)[1,:]) ]
          Y=[Y; vec(getVal(fg,id[2], api=api)[2,:]) ]
          # X=[X;vec(fg.v[id[2]].attributes["val"][1,:])]
          # Y=[Y;vec(fg.v[id[2]].attributes["val"][2,:])]
        end
      end
    end
  end
  return X,Y
end

function getAll2D(fg, sym; minnei::Int=0, api::DataLayerAPI=IncrementalInference.localapi)
  warn("getAll2D deprecated, use get2DSamples instead")
  return get2DSamples(fg, sym, minnei=minnei, api=api)
end

function get2DSampleMeans(fg::FactorGraph,
      sym;
      from::Int=0, to::Int=9999999999,
      minnei::Int=0,
      api::DataLayerAPI=IncrementalInference.localapi  )
  #
  X = Array{Float64,1}()
  Y = Array{Float64,1}()
  Th = Array{Float64,1}()
  LB = String[]

  # if sym = 'l', ignore single measurement landmarks
  allIDs = Array{Int,1}()
  for id in fg.IDs
    vertlbl = string(id[1])
    if vertlbl[1] == sym
      val = parse(Int,vertlbl[2:end])
      if from <= val && val <= to
        if length( getOutNeighbors(fg, id[2], api=api, needdata=true) ) >= minnei
          push!(allIDs, val)
        end
      end
    end
  end
  allIDs = sort(allIDs)

  for id in allIDs
    X=[X;Base.mean( vec( getVal(fg,string(sym,id), api=api)[1,:] ) )]
    Y=[Y;Base.mean( vec( getVal(fg, string(sym,id), api=api)[2,:] ) )]
    if sym == 'x'
      Th=[Th;Base.mean( vec( getVal(fg, string(sym,id), api=api)[3,:] ) )]
    end
    push!(LB, string(sym,id))
  end
  return X,Y,Th,LB
end

#draw landmark positions
function getAll2DMeans(fg, sym, api::DataLayerAPI=IncrementalInference.localapi)
  return get2DSampleMeans(fg, sym, api=api)
end

function getAll2DPoses(fg::FactorGraph, api::DataLayerAPI=IncrementalInference.localapi)
    return getAll2DSamples(fg, 'x', api=api)
end

function get2DPoseSamples(fg::FactorGraph; from::Int=0, to::Int=999999999, api::DataLayerAPI=IncrementalInference.localapi)
  return get2DSamples(fg, 'x'; from=from, to=to, api=api)
end

function get2DPoseMeans(fg::FactorGraph; from::Int=0, to::Int=999999999, api::DataLayerAPI=IncrementalInference.localapi)
  return get2DSampleMeans(fg, 'x', from=from, to=to, api=api)
end

function get2DPoseMax(fgl::FactorGraph;
            from::Int=-99999999999, to::Int=9999999999, api::DataLayerAPI=IncrementalInference.localapi )
  xLB,ll = ls(fgl) # TODO add: from, to, special option 'x'
  X = Array{Float64,1}()
  Y = Array{Float64,1}()
  Th = Array{Float64,1}()
  LB = String[]
  for slbl in xLB
    lbl = string(slbl)
    if from <= parse(Int,lbl[2:end]) <=to
      mv = getKDEMax(getVertKDE(fgl,slbl))
      push!(X,mv[1])
      push!(Y,mv[2])
      push!(Th,mv[3])
      push!(LB, lbl)
    end
  end
  return X, Y, Th, LB
end

function getAll2DLandmarks(fg::FactorGraph, minnei::Int=0, api::DataLayerAPI=IncrementalInference.localapi)
    return getAll2DSamples(fg, 'l', minnei=minnei, api=api)
end

function get2DLandmSamples(fg::FactorGraph; from::Int=0, to::Int=999999999, minnei::Int=0, api::DataLayerAPI=IncrementalInference.localapi)
  return get2DSamples(fg, 'l', from=from, to=to, minnei=minnei, api=api)
end

function get2DLandmMeans(fg::FactorGraph; from::Int=0, to::Int=999999999, minnei::Int=0, api::DataLayerAPI=IncrementalInference.localapi)
  return get2DSampleMeans(fg, 'l', from=from, to=to, minnei=minnei, api=api)
end

function removeKeysFromArr(fgl::FactorGraph, torm::Array{Int,1}, lbl::Array{String,1})
  retlbs = String[]
  for i in 1:length(lbl)
    id = parse(Int,lbl[i][2:end])
    if findfirst(torm,id) == 0
      push!(retlbs, lbl[i])
    else
      println("removeKeysFromArr -- skipping $(lbl[i]), id=$(id)")
    end
  end
  return retlbs
end

function get2DLandmMax(fgl::FactorGraph;
                from::Int=-99999999999,
                to::Int=9999999999,
                showmm=false, MM=Union{},
                api::DataLayerAPI=IncrementalInference.localapi  )
  #
  xLB,lLB = ls(fgl) # TODO add: from, to, special option 'x'
  if !showmm lLB = removeKeysFromArr(fgl, collect(keys(MM)), lLB); end
  X = Array{Float64,1}()
  Y = Array{Float64,1}()
  Th = Array{Float64,1}()
  LB = String[]
  for lb in lLB
    @show lb
    @show lbl = string(lb)
    if from <= parse(Int,lbl[2:end]) <=to
      mv = getKDEMax(getVertKDE(fgl,lb, api=api))
      push!(X,mv[1])
      push!(Y,mv[2])
      push!(LB, lbl)
    end
  end
  return X, Y, Th, LB
end

# convenience function to add DIDSON sonar constraints to graph
function addLinearArrayConstraint(fgl::FactorGraph,
      rangebearing::Union{Tuple{Float64, Float64}, Vector{Float64}},
      pose::Symbol,
      landm::Symbol ;
      rangecov::Float64=3e-4,
      bearingcov::Float64=3e-4 )

  #

  cl = LinearRangeBearingElevation((rangebearing[1],rangecov),(rangebearing[2],bearingcov))
  if !haskey(fgl.IDs, landm)
    pts = getVal(fgl, pose) + cl
    N = size(pts,2)
    vl1 = addNode!(fgl, landm,  pts,  N=N)
    println("Automatically added $(landm) to the factor graph")
  end
  addFactor!(fgl, [getVert(fgl, pose); getVert(fgl, landm)], cl)
  nothing
end

function addSoftEqualityPoint2D(fgl::FactorGraph, l1::Symbol, l2::Symbol;
    dist=MvNormal([0.0;0.0],eye(2)), ready::Int=1 )
  pp = Point2DPoint2D(dist)
  addFactor!(fgl, [l1,l2], pp, ready=ready)
end

#
# function +(p1::Rigid6DOF, p2::Rigid6DOF)
#     t = rotate(p1.rot, p2.trl) + p1.trl
#     return Rigid6DOF(p1.rot*p2.rot,t)
# end
#
# function +(p1::PoseSE3, p2::PoseSE3)
#     return PoseSE3(p1.utime+p2.utime,
#                    "+",
#                    p1.mu+p2.mu,
#                    p1.cov + p2.cov,
#                    Dict())
# end
#
# function makePose3(q::Quaternion=Quaternion(1.,zeros(3)), t::Array{Float64,1}=[0.,0,0]; ut=0, name="pose", cov=0.1*eye(6))
#     return PoseSE3(ut,
#                    name,
#                    Rigid6DOF( q, t ),
#                    cov,
#                    Dict())
# end

# function makeDynPose3(q::Quaternion=Quaternion(1.,zeros(3)), t::Array{Float64,1}=[0.,0,0]; ut=0, name="pose", cov=0.1*eye(18))
#     return DynPose3(ut,
#                     name,
#                     Dynamic6DOF( q, t, [0.,0,0], [0.,0,0] ),
#                     IMUComp([0.,0,0],[0.,0,0]),
#                     cov,
#                     Dict())
# end
#
# function wTo(p::PoseSE3)
#     T = eye(4)
#     T[1:3,1:3] =  convert(SO3, p.mu.rot).R
#     T[1:3,4] = p.mu.trl
#     return T
# end
#
# function oTw(p::PoseSE3)
#     return wTo(p) \ eye(4)
# end


text:
struct EnumerableIterable{T,S} <: Enumerable
    source::S
end

function query(source)
    IteratorInterfaceExtensions.isiterable(source) || error()
    typed_source = IteratorInterfaceExtensions.getiterator(source)
	T = eltype(typed_source)
    S = typeof(typed_source)

    source_enumerable = EnumerableIterable{T,S}(typed_source)

    return source_enumerable
end

Base.IteratorSize(::Type{EnumerableIterable{T,S}}) where {T,S} = haslength(S)

Base.eltype(::Type{EnumerableIterable{T,S}}) where {T,S} = T

Base.length(iter::EnumerableIterable{T,S}) where {T,S} = length(iter.source)

function Base.iterate(iter::EnumerableIterable{T,S}, state...) where {T,S}
    return iterate(iter.source, state...)
end


text:
using Test, GenericLinearAlgebra, LinearAlgebra

@testset "The General eigenvalue problem" begin

cplxord = t -> (real(t), imag(t))

@testset "General eigen problem with n=$n and element type=$T" for
    n in (10, 23, 100),
    T in (Float64, Complex{Float64})

    A = randn(T, n, n)
    vGLA    = GenericLinearAlgebra._eigvals!(copy(A))
    vLAPACK = eigvals(A)
    vBig    = eigvals(big.(A)) # not defined in LinearAlgebra so will dispatch to the version in GenericLinearAlgebra
    @test sort(vGLA, by = cplxord) ≈ sort(vLAPACK, by = cplxord)
    @test sort(vGLA, by = cplxord) ≈ sort(complex(eltype(A)).(vBig), by = cplxord)
end

@testset "make sure that solver doesn't hang" begin
    for i in 1:1000
        A = randn(8, 8)
        sort(abs.(GenericLinearAlgebra._eigvals!(copy(A)))) ≈ sort(abs.(eigvals(A)))
    end
end

@testset "Convergence in corner cases. Issue 29." begin
    function H(n::Int)
        H = zeros(2n, 2n)
        for i = 1 : 2 : 2n
            H[i, i+1] = 1
            H[i+1, i] = 1
        end
        H
    end

    function E(n::Int)
        E = zeros(2n, 2n)
        for i = 1 : (n - 1)
            E[2i + 1, 2i] = 1
        end
        E[1, 2n] = 1
        E
    end

    my_matrix(n::Int, η::Float64 = 1e-9) = H(n) .+ η .* E(n)

    A = my_matrix(4, 1e-3);
    @test sort(GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(A))), by = t -> (real(t), imag(t))) ≈ sort(eigvals(A), by = t -> (real(t), imag(t)))
end

@testset "Extract Schur factor" begin
    A = randn(5, 5)
    @test sum(eigvals(schur(A).T))     ≈ sum(eigvals(Float64.(schur(big.(A)).T)))
    @test sum(eigvals(schur(A).Schur)) ≈ sum(eigvals(Float64.(schur(big.(A)).Schur)))
end

end


text:
using Documenter
using GraphsMatching
import Graphs

makedocs(
    modules = [GraphsMatching],
    format = :html,
    sitename = "GraphsMatching",
	pages = Any[
		"Getting started"         => "index.md",
    ]
)

deploydocs(
    deps        = nothing,
    make        = nothing,
    repo        = "github.com/JuliaGraphs/GraphsMatching.jl.git",
    target      = "build",
    julia       = "0.6",
    osname      = "linux"
)