text:
const CONSTANT_PREFIX = "parameters"
const SWITCHER_PREFIX = "events_active"
const SWITCHER_SAVE_PREFIX = "events_save"
const SAVEAT_HEADER = Symbol("saveat[]")
const TSPAN_HEADER = Symbol("tspan")
const OBSERVABLES_HEADER = Symbol("observables[]")

# general interface
"""
    Scenario(model::Model;
      measurements::Vector{AbstractMeasurementPoint}=AbstractMeasurementPoint[],
      observables::Union{Nothing,Vector{Symbol}}=nothing,
      parameters::Vector{Pair{Symbol,Float64}} = Pair{Symbol,Float64}[],
      events_active::Union{Nothing, Vector{Pair{Symbol,Bool}}} = Pair{Symbol,Bool}[],
      events_save::Union{Tuple,Vector{Pair{Symbol, Tuple{Bool, Bool}}}} = (true,true),
      saveat::Union{Nothing,AbstractVector} = nothing,
      tspan::Union{Nothing,Tuple} = nothing,
      save_scope::Bool = true,
    )

Builds simulation scenario of type [`Scenario`](@ref)

Example: `Scenario(model; tspan = (0., 200.), saveat = [0.0, 150., 250.])`

Arguments:

- `model` : model of type [`Model`](@ref)
- `measurements` : `Vector` of measurements. Default is empty `Vector{AbstractMeasurementPoint}`
- `observables` : names of output observables. Overwrites default model's values. Default is `nothing`
- `parameters` : `Vector` of `Pair`s containing constants' names and values. Overwrites default model's values. Default is empty vector.
- `events_active` : `Vector` of `Pair`s containing events' names and true/false values. Overwrites default model's values. Default is empty `Vector{Pair}`
- `events_save` : `Tuple` or `Vector{Tuple}` marking whether to save solution before and after event. Default is `(true,true)` for all events
- `saveat` : time points, where solution should be saved. Default `nothing` values stands for saving solution at timepoints reached by the solver
- `tspan` : time span for the ODE problem
- `save_scope` : should scope be saved together with solution. Default is `true`
"""
function Scenario(
  model::Model;
  measurements::Vector{AbstractMeasurementPoint} = AbstractMeasurementPoint[],
  observables::Union{Nothing,Vector{Symbol}} = nothing,
  kwargs... # all arguments of build_ode_problem()
)
  # ODE problem
  prob = build_ode_problem(
    model;
    observables_ = observables,
    kwargs...
  )

  return Scenario(model.init_func, prob, measurements)
end

# Scenario struct method
function add_scenarios!(p, scen::Vector{P}) where P <: Pair
  [push!(scenarios(p), s) for s in scen]
  return nothing
end

# CSV methods
function add_scenarios!(
  platform::Platform,
  vector::DataFrames.DataFrameRows;
  subset::AbstractVector{P} = Pair{Symbol, Symbol}[]
) where P <: Pair{Symbol, Symbol}
  selected_rows = _subset(vector, subset)

  for row in selected_rows
    _add_scenario!(platform, row)
  end

  return nothing
end

# DataFrame methods
"""
    add_scenarios!(
      platform::Platform,
      df::DataFrame;
      subset::AbstractVector{P} = Pair{Symbol, Symbol}[]
    ) where P <: Pair{Symbol, Symbol}

Adds a new `Scenario` to the `Platform`

Arguments:

- `platform` : platform of [`Platform`](@ref) type
- `df` : `DataFrame` with scenarios setup, typically obtained with [`read_scenarios`](@ref) function
- `subset` : subset of scenarios which will be added to the `platform`. Default `Pair{Symbol, Symbol}[]` adds all scenarios from the `df`
"""
function add_scenarios!(
  platform::Platform,
  df::DataFrame;
  kwargs...
)
  add_scenarios!(platform, eachrow(df); kwargs...)
end

# private function to add one scenario row into platform

function _add_scenario!(platform::Platform, row::Any) # maybe not any
  _id = row[:id]
  model_name = get(row, :model, :nameless)

  if !haskey(platform.models, model_name)
    @warn "Lost model $(model_name). Scenario $_id has been skipped."
    return nothing # BRAKE
  else
    model = platform.models[model_name]
  end

  # iterate through parameters
  _parameters = Pair{Symbol,Float64}[]
  _events_active = Pair{Symbol,Bool}[]
  _events_save = Pair{Symbol, Tuple{Bool, Bool}}[]
  for key in keys(row)
    if !ismissing(row[key])
      splitted_key = split(string(key), ".")
      if splitted_key[1] == CONSTANT_PREFIX
        push!(_parameters, Symbol(splitted_key[2]) => row[key])
      elseif splitted_key[1] == SWITCHER_PREFIX
        push!(_events_active, Symbol(splitted_key[2]) => bool(row[key]))
      elseif splitted_key[1] == SWITCHER_SAVE_PREFIX
        save_evt_vec = split(row[key], ";")
        @assert length(save_evt_vec) == 2 "Events saving setup accepts two values (e.g. true;false). Check the scenarios table."
        push!(_events_save, Symbol(splitted_key[2]) => (bool(save_evt_vec[1]), bool(save_evt_vec[2])))
      end
    end
  end

  if haskey(row, SAVEAT_HEADER) && !ismissing(row[SAVEAT_HEADER])
    save_times = row[SAVEAT_HEADER]
    _saveat = typeof(save_times) <: Number ? [Float64(save_times)] : parse.(Float64, split(save_times, ";"))
  else
    _saveat = nothing
  end

  if haskey(row, TSPAN_HEADER) && !ismissing(row[TSPAN_HEADER])
    _tspan = (0., row[TSPAN_HEADER])
  else
    _tspan = nothing
  end

  if haskey(row, OBSERVABLES_HEADER) && !ismissing(row[OBSERVABLES_HEADER])
    observables_str = split(row[OBSERVABLES_HEADER], ";")
    _observables = Symbol.(observables_str)
  else
    _observables = nothing
  end

  scenario = Scenario(
    model;
    parameters = _parameters,
    events_active = _events_active,
    events_save = _events_save,
    saveat = _saveat,
    tspan = _tspan,
    observables = _observables
  )

  push!(platform.scenarios, _id => scenario)
end

# helper to read from csv and xlsx

function read_scenarios_csv(filepath::String; kwargs...)
  df = DataFrame(CSV.File(
    filepath,
    types = Dict(:id => Symbol, :model=>Symbol, :tspan => Float64);
    kwargs...)
  )
  assert_scenarios(df)

  return df
end

function read_scenarios_xlsx(filepath::String, sheet=1; kwargs...)
  df = DataFrame(XLSX.readtable(filepath, sheet,infer_eltypes=true)...)
  assert_scenarios(df)

  df[!,:id] .= Symbol.(df[!,:id])
  "tspan" in names(df) && (df[!,:tspan] .= float64.(df[!,:tspan]))
  "model" in names(df) && (df[!,:model] .= Symbol.(df[!,:model]))
  return df
end

"""
    read_scenarios(filepath::String, sheet=1; kwargs...)

Reads table file with scenarios to `DataFrame`

Arguments:

- `filepath` : path to table file. Supports ".csv" and ".xlsx" files
- `sheet` : number of sheet in case of ".xlsx" file. Default is `1`
- kwargs : other arguments supported by `CSV.File`
"""
function read_scenarios(filepath::String, sheet=1; kwargs...)
  ext = splitext(filepath)[end]

  if ext == ".csv"
    df = read_scenarios_csv(filepath; kwargs...)
  elseif ext == ".xlsx"
    df = read_scenarios_xlsx(filepath, sheet)
  else
    error("Extension $ext is not supported.")
  end
  return df
end

function assert_scenarios(df)
  names_df = names(df)
  for f in ["id"]
    @assert f ∈ names_df "Required column name $f is not found in measurements table."
  end
  return nothing
end


text:
## Dependencies

#To install Laplacians you need to put in Julia's Prompt "add Laplacians#master"
#For some reason Pkg.add() doesn't work with this.
#Pkg.add("MATLAB")
#using Laplacians
#using Pkg
#include(string(Pkg.dir("Laplacians") , "/src/matlabSolvers.jl"))
getlagged(x) = [NaN; x[1:(end - 1)]]

using DocStringExtensions

#Defining types and structures
abstract type AbstractLLSAlgorithm end
abstract type AbstractGraphLLSAlgorithm <: AbstractLLSAlgorithm end  # i.e. requires graph laplacian
struct CMGPreconditionedLLS <: AbstractGraphLLSAlgorithm  end
struct AMGPreconditionedLLS <: AbstractGraphLLSAlgorithm  end
struct DirectLLS <: AbstractLLSAlgorithm  end  # no graph required

abstract type AbstractLeverageAlgorithm end

"""
$(TYPEDEF)

Data type to pass to VCHDFESettings type, to indicate Exact algorithm
"""
struct ExactAlgorithm <: AbstractLeverageAlgorithm  end

"""
$(TYPEDEF)

Data type to pass to VCHDFESettings type, to indicate JLA algorithm

### Fields

* `num_simulations`: number of simulations in estimation. If num_simulations = 0, defaults to 100 * log(#total fixed effect)"
"""
@with_kw struct JLAAlgorithm <: AbstractLeverageAlgorithm
    num_simulations::Int64 = 0
end

"""
$(TYPEDEF)

The VCHDFESettings type is to pass information to methods regarding which algorithm to use.

### Fields

* `cg_maxiter`: maximum number of iterations (default = 300)
* `leave_out_level`: leave-out level (default = match)
* `leverage_algorithm`: which type of algorithm to use (default = JLAAlgorithm())
* `first_id_effects`: includes first id effects. At this version it is required to include the first_id_effects. (default = true)
* `cov_effects`: includes covariance of first-second id effects. At this version it is required to include the cov_effects. (default = true)
* `print_level`: prints the state of the program in std output. If print_level = 0, the app prints nothing in the std output. (default = 1)
* `first_id_display_small`: name of the first id in lower cases (default = person)
* `first_id_display`: name of the first id (default = Person)
* `second_id_display_small`: name of the second id in lower cases (default = firm)
* `second_id_display`: name of the second id (default = Firm)
* `outcome_id_display_small`: name of the observation id in lower cases (default = wage)
* `outcome_id_display`: name of the observation id (default = Wage)
"""
@with_kw struct VCHDFESettings{LeverageAlgorithm}
    cg_maxiter::Int64 = 300
    leverage_algorithm::LeverageAlgorithm = JLAAlgorithm()
    leave_out_level::String = "match"
    first_id_effects::Bool = true
    cov_effects::Bool = true
    print_level::Int64 = 1
    first_id_display_small::String = "person"
    first_id_display::String = "Person"
    second_id_display_small::String = "firm"
    second_id_display::String = "Firm"
    outcome_id_display_small::String = "wage"
    outcome_id_display::String = "Wage"
end

## FUNCTIONS FROM THE VECTORIZEDROUTINE.JL PACKAGE

function accumarray(subs, val, sz=(maximum(subs),))
    A = zeros(eltype(val), sz...)
    for i = 1:length(val)
        @inbounds A[subs[i]] += val[i]
    end
    A
end

function accumarray(subs, val::Number, sz=(maximum(subs),))
    A = zeros(eltype(val), sz...)
    for i = 1:length(subs)
        @inbounds A[subs[i]] += val[i]
    end
    A
end

## Methods and Types Definitions

#1) Finds AKM largest connected set
"""
$(SIGNATURES)

Returns a tuple of observation belonging to the largest connected set with the corresponding identifiers and outcomes. This requires to have the data sorted by first identifier, and time period (e.g. we sort by worked id and year). This is also the set where we can run AKM models with the original data.

### Arguments
* `y`: outcome (e.g. log wage)
* `first_id`: first identifier (e.g. worker id)
* `second_id`: second identifier (e.g. firm id)
* `settings`: settings based on data type `VCHDFESettings`. Please see the reference provided below.
"""
function find_connected_set(y, first_idvar, second_idvar, settings)

    #Observation identifier to join later the FE
    obs_id = collect(1:size(y,1))

    seconds = unique(sort(second_idvar))
    second_id = indexin(second_idvar, seconds)

    firsts = unique(sort(first_idvar))
    first_id = indexin(first_idvar, firsts)

    #Save ids just in case
    second_id_old=second_idvar;
    first_id_old=first_idvar;

    n_seconds=length(seconds)
    n_firsts=length(firsts)
    second_id=second_id.+n_firsts

    graph_size=length(seconds)+length(firsts)
    G=Graph(graph_size)
    for i in 1:size(y,1)
        add_edge!(G, second_id[i], first_id[i])
    end

    cs=connected_components(G)
    (settings.print_level > 1) && println("Largest connected set has been found")
    pos = indexin(  [  maximum(size.(cs,1))] , size.(cs,1) )[1]
    lcs=cs[pos]

    connected_seconds=lcs[lcs.>n_firsts]

    sel=findall(in(connected_seconds),second_id)

    obs_id = [obs_id[x] for x in sel ]
    yvec = [y[x] for x in sel ]
    second_id = [second_idvar[x] for x in sel ]
    first_id = [first_idvar[x] for x in sel ]

    #Relabel seconds
    seconds = unique(sort(second_id))
    second_id = indexin(second_id, seconds)

    #Relabel firsts
    first_ids = unique(sort(first_id))
    first_id = indexin(first_id, first_ids)

    return (obs_id = obs_id , y = yvec , first_id = first_id, second_id = second_id )

end

#2) Pruning and finding Leave-Out Largest connected set
"""
$(SIGNATURES)

This function prunes the dataset from articulation points. If the first identifier is worker id it means that it prunes workers that would disconnect the graph if they were dropped.

### Arguments
* `yvec`: outcome (e.g. log wage)
* `first_idvar`: first identifier (e.g. worker id)
* `second_idvar`: second identifier (e.g. firm id)
* `obs_id`: observation identifier.
* `settings`: settings based on data type `VCHDFESettings`. Please see the reference provided below.
"""
function prunning_connected_set(yvec, first_idvar, second_idvar, obs_id, settings)

    seconds = unique(sort(second_idvar))
    second_id = indexin(second_idvar, seconds)
    firsts = unique(sort(first_idvar))
    first_id = indexin(first_idvar, firsts)

    nbadfirsts=1
    while nbadfirsts>0

        n_seconds=length(seconds)
        n_firsts=length(firsts)
        second_id=second_id.+n_firsts

        graph_size=n_firsts + n_seconds
        lcs_graph=Graph(graph_size)
            for i in 1:size(yvec,1)
                add_edge!(lcs_graph, second_id[i], first_id[i])
            end

        #get articulation vertex
        artic_vertex = articulation(lcs_graph)

        sel=findall(x->x==nothing, indexin(first_id,artic_vertex))

        bad_firsts = artic_vertex[artic_vertex.<=n_firsts]
        nbadfirsts = size(bad_firsts,1)

        settings.print_level > 1 && println("Number of $(settings.first_id_display_small) that are articulation points: ", nbadfirsts)

        #Restrict the sample
        yvec = [yvec[x] for x in sel ]
        second_id = [second_id[x] for x in sel ]
        first_id = [first_id[x] for x in sel ]
        obs_id = [obs_id[x] for x in sel ]

        #Relabel seconds
        seconds = unique(sort(second_id))
        second_id = indexin(second_id, seconds)

        #Relabel firsts
        first_ids = unique(sort(first_id))
        first_id = indexin(first_id, first_ids)

        n_firsts=maximum(first_id)
        n_seconds=maximum(second_id);

        second_id = second_id .+ n_firsts

        #Constructing new Graph
        G=Graph(n_firsts+n_seconds)
        Nprimeprime = size(sel,1)
        for i in 1:Nprimeprime
            add_edge!(G, second_id[i], first_id[i])
        end

        #Find Largest Connected Set (LCS)
        cs=connected_components(G)

        pos = indexin(  [  maximum(size.(cs,1))] , size.(cs,1) )[1]
        lcs=cs[pos]

        connected_seconds=lcs[lcs.>n_firsts]

        sel=findall(in(connected_seconds),second_id)

        obs_id = [obs_id[x] for x in sel ]
        yvec = [yvec[x] for x in sel ]
        second_id = [second_id[x] for x in sel ]
        first_id = [first_id[x] for x in sel ]

        #Relabel seconds
        seconds = unique(sort(second_id))
        second_id = indexin(second_id, seconds)

        #Relabel firsts
        first_ids = unique(sort(first_id))
        first_id = indexin(first_id, first_ids)

    end

    return (obs_id = obs_id , y = yvec , first_id = first_id, second_id = second_id )
end

#3) Drops Single Observations
"""
$(SIGNATURES)

This function drops observations that correspond with first identifiers with a single observation. For example, if first identifier is worker id, it will drop observations for workers that only appear once in the data.

### Arguments
* `yvec`: outcome (e.g. log wage)
* `first_idvar`: first identifier (e.g. worker id)
* `second_idvar`: second identifier (e.g. firm id)
* `obs_id`: observation identifier.
"""
function drop_single_obs(yvec, first_idvar, second_idvar,obs_id)

    seconds = unique(sort(second_idvar))
    second_id = indexin(second_idvar, seconds)
    firsts = unique(sort(first_idvar))
    first_id = indexin(first_idvar, firsts)

    T=accumarray(first_id,ones(length(first_id)))
    T = [T[x] for x in first_id]
    sel = T.>1
    sel = findall(x->x==true,sel )

    obs_id = [obs_id[x] for x in sel ]
    yvec = [yvec[x] for x in sel ]
    second_id = [second_id[x] for x in sel ]
    first_id = [first_id[x] for x in sel ]

    #Relabel seconds
    seconds = unique(sort(second_id))
    second_id = indexin(second_id, seconds)

    #Relabel firsts
    first_ids = unique(sort(first_id))
    first_id = indexin(first_id, first_ids)

    return (obs_id = obs_id , y = yvec , first_id = first_id, second_id = second_id )
end

#4) Compute Movers
"""
$(SIGNATURES)

Returns a vector that indicates whether the `first_id` (e.g. worker) is a mover across `second_id` (e.g. firms), as well as a vector with the number of periods that each `first_id` appears.

### Arguments
* `first_id`: first identifier (e.g. worker id)
* `second_id`: second identifier (e.g. firm id)
"""
function compute_movers(first_id,second_id)

    gcs = [NaN; first_id[1:end-1]]
    gcs = first_id.!=gcs

    lagsecond_id=[NaN; second_id[1:end-1]]
    for x in findall(x->x ==true , gcs)
        lagsecond_id[x] = NaN
    end

    stayer=(second_id.==lagsecond_id)
    for x in findall(x->x ==true , gcs)
        stayer[x] = true
    end

    stayer=Int.(stayer)
    stayer=accumarray(first_id,stayer)
    T=accumarray(first_id,ones(length(first_id)))
    stayer=T.==stayer
    movers=stayer.==false

    movers = [movers[x] for x in first_id]
    T = [T[x] for x in first_id]

    return (movers = movers, T = T)
end

#5) Creates match first_id using second_id id
"""
$(SIGNATURES)

Computes a match identifier for every combination of first and second identifier. For example, this can be the match identifier of worker-firm combinations.

### Arguments
* `first_id`: first identifier (e.g. worker id)
* `second_id`: second identifier (e.g. firm id)
"""
function compute_matchid(second_id,first_id)
    match_id2 = string.(first_id).*string.("+",second_id)
    match_id2 = indexin(match_id2, unique(match_id2))

    return match_id2
end

#6) Compute Leave-Out Connected Set : No worker is an articulation vertex and no single observations
"""
$(SIGNATURES)

Returns a tuple with the observation number of the original dataset that belongs to the Leave-out connected set as described in Kline,Saggio, Solvesten. It also provides the corresponding outcome and identifiers in this connected set.

### Arguments
* `y`: outcome vector
* `first_id`: first identifier (e.g. worker id)
* `second_id`: second identifier (e.g. firm id)
* `settings`: settings based on `VCHDFESettings`
* `controls`: at this version only `controls=nothing` is supported.
"""
function get_leave_one_out_set(y, first_id, second_id, settings, controls)
    # @assert settings.first_id_effects == true && settings.cov_effects == true

    # compute y, id firmid, controls, settings
    # compute y, first_id second_id, controls, settings
    (settings.print_level > 1) && println("\nFinding the leave-one-out connected set")
    @unpack obs_id,  y  , first_id , second_id  = find_connected_set(y,first_id,second_id,settings)
    @unpack obs_id,  y  , first_id , second_id  = prunning_connected_set(y,first_id,second_id, obs_id,settings)
    @unpack obs_id,  y  , first_id , second_id  = drop_single_obs(y,first_id,second_id, obs_id)
    controls == nothing ? nothing : controls = controls[obs_id,:]

    return (obs = obs_id, y = y, first_id = first_id, second_id = second_id, controls = controls)
end

#7) Main Routine Function : Performs Leave Out Estimation and Inference
"""
$(SIGNATURES)

Returns a tuple with the observation number of the original dataset that belongs to the Leave-out connected set as described in Kline,Saggio, Solvesten. It also provides the corresponding outcome and identifiers in this connected set.

### Arguments
* `y`: outcome vector
* `first_id`: first identifier (e.g. worker id)
* `second_id`: second identifier (e.g. firm id)
* `controls`: covariates that will be partialled out from outcome before it performs KSS.
* `do_lincom`: boolean that indicates whether it runs inference.
* `Z_lincom`: matrix of covariates to be used in lincom inference.
* `lincom_labels`: vector of labels of the columns of Z_lincom.
* `settings`: settings based on `VCHDFESettings`
* `controls`: at this version only `controls=nothing` is supported.
"""
function leave_out_KSS(y,first_id,second_id;controls = nothing, do_lincom = false, Z_lincom = nothing , lincom_labels = nothing, settings = VCHDFESettings())

    algo = typeof(settings.leverage_algorithm) == JLAAlgorithm ? "JLA" : "Exact"
    println("Running KSS correction with the following options")
    println("Leave Out Strategy : Leave $(settings.leave_out_level) out")
    if algo == "JLA"
        simul = settings.leverage_algorithm.num_simulations == 0 ? 200 : settings.leverage_algorithm.num_simulations
        println("Algorithm for Computation of Statistical Leverages : $(algo) with $(simul) simulations")
    else
        println("Algorithm for Computation of Statistical Leverages : $(algo) \n")
    end

    first_id_old = first_id
    second_id_old = second_id
    y_untransformed = y

    #Compute Leave Out Connected Set
    @unpack obs,  y  , first_id , second_id, controls = get_leave_one_out_set(y, first_id, second_id, settings, controls)

    if settings.print_level>1
        num_movers = length(unique(compute_movers(first_id,second_id).movers .* first_id)) - 1

        summary = """
        Information on Leave Out Connected Sample \n
        Number of observations: $(length(obs))
        Number of $(settings.first_id_display_small)s: $(maximum(first_id))
        Number of $(settings.second_id_display_small)s: $(maximum(second_id))
        Number of Movers : $(num_movers)
        Mean of $(settings.outcome_id_display): $(mean(y))
        Variance of $(settings.outcome_id_display): $(var(y))
        """
        println(summary)
    end

    #Residualize outcome variable
    if controls != nothing
        println("\nPartialling out controls from $(settings.outcome_id_display)...")
        NT = size(y,1)
        J = maximum(second_id)
        N = maximum(first_id)
        K = size(controls,2)
        nparameters = N + J + K

        D = sparse(collect(1:NT),first_id,1)
        F = sparse(collect(1:NT),second_id,1)
        S= sparse(1.0I, J-1, J-1)
        S=vcat(S,sparse(-zeros(1,J-1)))
        X = hcat(D, -F*S, controls)

        #My best shot is to wrap AMG as LinearOperator
        buff = zeros(size(X,2))
        xx = X'*X
        P = AmgOperator(ruge_stuben(xx),buff)
        xy=X'*y
        beta, stats = Krylov.cg(xx,[xy...]; M = P , rtol = 1e-6, itmax = 300)

        y=y-X[:,N+J:end]*beta[N+J:end]
        controls = nothing
        println("Partialling out completed.\n")
    end

    @unpack θ_first, θ_second, θCOV, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov, y, X, sigma_i = leave_out_estimation(y,first_id,second_id,controls,settings)

    if do_lincom == true
        #Subset to the Leave Out Sample
        Z_lincom = Z_lincom[obs,:]
        F = sparse(collect(1:length(second_id)),second_id,1)
        J = size(F,2)
        S = sparse(1.0I, J-1, J-1)
        S = vcat(S,sparse(-zeros(1,J-1)))
        Transform = hcat(spzeros(length(second_id),maximum(first_id)), -F*S )

        #Collapse and reweight to person-year observations
        match_id = compute_matchid(second_id, first_id)
        Z_lincom_col = ones(size(Z_lincom,1),1)
        for i = 1:size(Z_lincom,2)
            Z_lincom_col = hcat(Z_lincom_col,(transform(groupby(DataFrame(z = Z_lincom[:,i], match_id = match_id), :match_id), :z => mean  => :z_py).z_py))
        end

        #Run Inference
        println("\nRegressing the $(settings.second_id_display_small) effects on observables Z.")
        @unpack test_statistic, linear_combination , SE_linear_combination_KSS, SE_naive = lincom_KSS(y,X, Z_lincom_col, Transform, sigma_i; lincom_labels)
    end

    return (θ_first = θ_first, θ_second = θ_second , θCOV = θCOV)
end

"""
$(SIGNATURES)

Computes variance of errors for stayers when leaving out a match.

### Arguments
* `y`: outcome variable
* `first_id`: first identifier (e.g. worker id)
* `second_id`: second identifier (e.g. firm id)
* `weight`: spell length of every match.
* `b`: fixed effects coefficients vector
"""
function sigma_for_stayers(y,first_id, second_id, weight, b)
    #Go back to person-year space
    first_id_weighted = vcat(fill.(first_id, weight)...)
    second_id_weighted = vcat(fill.(second_id, weight)...)

    #Compute Pii for Stayers = inverse of # of obs
    T = accumarray(first_id_weighted, 1)
    Pii = 1 ./T
    Mii = 1 .- Pii[first_id_weighted]

    #Compute OLS residual
    NT = size(y,1)
    D = sparse(collect(1:NT),first_id_weighted, 1)
    F = sparse(collect(1:NT),second_id_weighted, 1)
    J = size(F,2)
    S = sparse(1.0I, J-1, J-1)
    S = vcat(S,sparse(-zeros(1,J-1)))
    X = hcat(D, -F*S)

    eta = y - X*b
    eta_h = eta./Mii
    sigma_stayers = (y .- mean(y)).*eta_h

    #Collapse to match
    match_id = compute_matchid(second_id_weighted,first_id_weighted)
    sigma_stayers = (combine(groupby(DataFrame(sigma = sigma_stayers , match_id = match_id), :match_id), :sigma => mean).sigma_mean)

    return sigma_stayers
end

"""
$(SIGNATURES)

Computes a KSS quadratic form to correct bias.

### Arguments
* `sigma_i`: variance estimator of observation i.
* `A_1`: left matrix of quadratic form
* `A_2`: right matrix of quadratic form
* `beta`: fixed effects coefficients vector
* `Bii`: Bii correction elements.
"""
function kss_quadratic_form(sigma_i, A_1, A_2, beta, Bii)
    right                               = A_2*beta
    left                                = A_1*beta
    theta                               = cov(left,right)
    theta                               = theta[1]
    dof                                 = size(left,1)-1
    theta_KSS                           = theta-(1/dof)*sum(Bii.*sigma_i)
end

"""
$(SIGNATURES)

Returns the bias-corrected components, the vector of coefficients, the corresponding fixed effects for every observation, and the diagonal matrices containing the Pii and Biis.

### Arguments
* `y`: outcome vector
* `first_id`: first identifier (e.g. worker id)
* `second_id`: second identifier (e.g. firm id)
* `settings`: settings based on `VCHDFESettings`
* `controls`: matrix of control variables. At this version it doesn't work properly for very large datasets.
"""
function leave_out_estimation(y,first_id,second_id,controls,settings)

    #Create matrices for computations
    NT = size(y,1)
    J = maximum(second_id)
    N = maximum(first_id)
    K = controls == nothing ? 0 : size(controls,2)
    nparameters = N + J + K

    match_id = compute_matchid(first_id, second_id)

    #first (worker) Dummies
    D = sparse(collect(1:NT),first_id,1)

    #second (firm) Dummies
    F = sparse(collect(1:NT),second_id,1)

    # N+J x N+J-1 restriction matrix
    S = sparse(1.0I, J-1, J-1)
    S = vcat(S,sparse(-zeros(1,J-1)))

    X = hcat(D, -F*S)

    #SET DIMENSIONS
    n=size(X,1)

    # PART 1: ESTIMATE HIGH DIMENSIONAL MODEL
    xx=X'*X
    xy=X'*y
    compute_sol = approxcholSolver(xx;verbose = settings.print_level > 0)
    beta = compute_sol([xy...];verbose=false)
    eta=y-X*beta

    pe=D * beta[1:N]
    fe=F*S * beta[N+1:N+J-1]

    println("\n","Plug-in Variance Components:")

    var_den = var(y)
    σ2_ψ_AKM = var(fe)
    println("Plug-in Variance of $(settings.second_id_display_small) Effects: ", σ2_ψ_AKM )
    σ2_α_AKM = var(pe)
    println("Plug-in Variance of $(settings.first_id_display_small) Effects: ", σ2_α_AKM )
    σ2_ψα_AKM = cov(pe,-fe)
    println("Plug-in Covariance of $(settings.first_id_display_small)-$(settings.second_id_display_small) Effects: ", σ2_ψα_AKM)
    println("Correlation of $(settings.first_id_display_small)-$(settings.second_id_display_small) Effects: ", σ2_ψα_AKM/(sqrt(σ2_ψ_AKM)*sqrt(σ2_α_AKM)))
    println("Fraction of Variance explained by  $(settings.first_id_display_small)-$(settings.second_id_display_small) Effects: ", (σ2_ψ_AKM+2*σ2_ψα_AKM+σ2_ψα_AKM)/var_den, "\n" )

    #Part 2: Collapse & Reweight (if needed)
    weight = ones(NT,1)
    y_py = y
    if settings.leave_out_level == "match"
        #Compute Weights
        weight = accumarray(match_id, 1)

        first_id_weighted = first_id
        second_id_weighted = second_id

        #Collapse at match level averages
        first_id = Int.(combine(groupby(DataFrame(first_id = first_id, match_id = match_id), :match_id), :first_id => mean).first_id_mean)
        second_id = Int.(combine(groupby(DataFrame(second_id = second_id, match_id = match_id), :match_id), :second_id => mean).second_id_mean)
        y = (combine(groupby(DataFrame(y = y, match_id = match_id), :match_id), :y => mean).y_mean)

        #At person-year level (like a Matlab's repelem)
        first_id_weighted = vcat(fill.(first_id,weight)...)
        second_id_weighted = vcat(fill.(second_id,weight)...)

        #Build Design
        NT = size(y,1)
        J = maximum(second_id)
        N = maximum(first_id)
        nparameters = N + J + K
        D = sparse(collect(1:NT),first_id,1)
        F = sparse(collect(1:NT),second_id,1)
        S = sparse(1.0I, J-1, J-1)
        S = vcat(S,sparse(-zeros(1,J-1)))

        X = hcat(D, -F*S)
        Dvar = hcat( sparse(collect(1:length(match_id)),first_id_weighted,1) , spzeros(length(match_id),J-1))
        Fvar = hcat(spzeros(length(match_id),N), -1*sparse(collect(1:length(match_id)), second_id_weighted,1  )*S )

        #Weighting
        weight_mat = sparse(collect(1:NT), collect(1:NT), weight.^(0.5) , NT, NT )
        X = weight_mat * X
        y = weight_mat * y

    else
        #Build Design
        NT = size(y,1)
        J = maximum(second_id)
        N = maximum(first_id)
        nparameters = N + J + K
        D = sparse(collect(1:NT),first_id,1)
        F = sparse(collect(1:NT),second_id,1)
        S= sparse(1.0I, J-1, J-1)
        S=vcat(S,sparse(-zeros(1,J-1)))

        X = hcat(D, -F*S)
        Dvar = hcat(  D, spzeros(NT,J-1) )
        Fvar = hcat(spzeros(NT,N), -F*S )

    end

    #Part 3: Compute Pii, Bii
    @unpack Pii , Mii  , correction_JLA , Bii_first , Bii_second , Bii_cov = leverages(settings.leverage_algorithm, X, Dvar, Fvar, settings)

    (settings.print_level > 1) && println("Pii and Bii have been computed.")

    #Compute Leave-out residual
    xx=X'*X
    xy=X'*y
    compute_sol = approxcholSolver(xx;verbose = settings.print_level > 0)
    beta = compute_sol([xy...];verbose=false)
    eta=y-X*beta

    eta_h = eta ./ Mii
    sigma_i = ( ( y .- mean(y) ) .* eta_h ) .* correction_JLA

    if settings.leave_out_level == "match"
        T = accumarray(first_id,1)
        stayers = T.==1
        stayers = [stayers[j] for j in first_id]

        sigma_stayers = sigma_for_stayers(y_py, first_id, second_id, weight, beta)
        sigma_i[stayers] .= [sigma_stayers[j] for j in findall(x->x==true,stayers )]
    end

    #Compute bias corrected variance comp of second (Firm) Effects
    θ_second = kss_quadratic_form(sigma_i, Fvar, Fvar, beta, Bii_second)

    θ_first = settings.first_id_effects==true  ? kss_quadratic_form(sigma_i, Dvar, Dvar, beta, Bii_first) : nothing

    θCOV = settings.cov_effects==true ? kss_quadratic_form(sigma_i, Fvar, Dvar, beta, Bii_cov) : nothing

    #Pii = diag(Pii)

    #Bii_second = diag(Bii_second)

    #Bii_first = settings.first_id_effects == true ? diag(Bii_first) : nothing

    #Bii_cov = settings.cov_effects == true ? diag(Bii_cov) : nothing

    #TODO print estimates
    if settings.print_level > 0
        println("Bias-Corrected Variance Components:")
        println("Bias-Corrected Variance of $(settings.second_id_display_small) Effects: ", θ_second)
        (settings.first_id_effects > 0) && println("Bias-Corrected Variance of $(settings.first_id_display_small) Effects: ", θ_first)
        (settings.cov_effects > 0) && println("Bias-Corrected Covariance of $(settings.first_id_display_small)-$(settings.second_id_display_small) Effects: ", θCOV)
        θ_first = θ_first == nothing ? nothing : abs(θ_first) #This is only required for the very small dataset that does the precompilation
        θ_second = θ_second == nothing ? nothing : abs(θ_second) #This is only required for the very small dataset that does the precompilation
        (settings.cov_effects > 0) && (settings.first_id_effects > 0) && println("Bias-Corrected Correlation of $(settings.first_id_display_small)-$(settings.second_id_display_small) Effects: ", θCOV/(sqrt(θ_first)*sqrt(θ_second)) )
        (settings.cov_effects > 0) && (settings.first_id_effects > 0) && println("Bias-Corrected Fraction of Variance explained by  $(settings.first_id_display_small)-$(settings.second_id_display_small) Effects: ", (θ_second+2*θCOV+θ_first)/var_den)
    end

    return (θ_first = θ_first, θ_second = θ_second, θCOV = θCOV, β = beta, Dalpha = pe, Fpsi = fe, Pii = Pii, Bii_first = Bii_first,
            Bii_second = Bii_second, Bii_cov = Bii_cov, y = y , X = X, sigma_i = sigma_i )
end

"""
$(SIGNATURES)

This function computes the diagonal matrices containing Pii and Bii under Exact Algorithm. See appendix in KSS for more information.

### Arguments
* `lev`: an instance of Exact algorithm structure.
* `X`: the design matrix in the linear model.
* `Dvar`: matrix with incidence information of first_id
* `Fvar`: matrix with incidence information of second_id
* `settings`: settings based on data type `VCHDFESettings`. Please see the reference provided below.
"""
function leverages(lev::ExactAlgorithm, X,Dvar,Fvar, settings)

    M = size(X,1)

    #Define solver
    S_xx = X'*X

    # Create the solvers
    ldli, la = computeLDLinv(S_xx)
    buffs = zeros(size(la)[1],Threads.nthreads())
    compute_sol = []
    for i in 1:Threads.nthreads()
        P = approxcholOperator(ldli,buffs[:,i])
        push!(compute_sol,approxcholSolver(P,la))
    end

    #Initialize output
    Pii = zeros(M)
    Bii_second = zeros(M)
    Bii_cov= settings.cov_effects ==true ? zeros(M) : nothing
    Bii_first= settings.first_id_effects == true ? zeros(M) : nothing

    Threads.@threads for i=1:M

        #Only one inversion needed for exact alg
        zexact = compute_sol[Threads.threadid()]( [X[i,:]...] ; verbose=false)

        #Compute Pii
        Pii[i] = X[i,:]'*zexact

        #Compute Bii
        COV = cov(Fvar*zexact,Fvar*zexact)
        Bii_second[i] = COV[1]*(size(Fvar,1)-1)

        if Bii_first != nothing
            COV = cov(Fvar*zexact,Dvar*zexact)
            Bii_first[i] = COV[1]*(size(Fvar,1)-1)
        end

        if Bii_cov != nothing
            COV = cov(Dvar*zexact,Dvar*zexact)
            Bii_cov[i] = COV[1]*(size(Dvar,1)-1)
        end

    end

    #Censor
    #Pii[ findall(Pii.>=0.99)] .= 0.99

    correction_JLA = 1
    Mii = 1 .- Pii

    return (Pii = Pii , Mii = Mii , correction_JLA = correction_JLA, Bii_first = Bii_first , Bii_second = Bii_second , Bii_cov = Bii_cov)
end

"""
$(SIGNATURES)

This function computes the diagonal matrices containing Pii and Bii under Johnson-Linderstrauss Algorithm. See appendix in KSS for more information.

### Arguments
* `lev`: an instance of JLA algorithm structure.
* `X`: the design matrix in the linear model.
* `Dvar`: matrix with incidence information of first_id
* `Fvar`: matrix with incidence information of second_id
* `settings`: settings based on data type `VCHDFESettings`. Please see the reference provided below.
"""
function leverages(lev::JLAAlgorithm, X,Dvar,Fvar, settings)

    NT = size(X,1)
    FE = size(X,2)

    p = lev.num_simulations == 0 ? 200 : lev.num_simulations

    #Define solver
    S_xx = X'*X

    # Create the solvers
    ldli, la = computeLDLinv(S_xx)
    buffs = zeros(size(la)[1],Threads.nthreads())
    compute_sol = []
    for i in 1:Threads.nthreads()
        P = approxcholOperator(ldli,buffs[:,i])
        push!(compute_sol,approxcholSolver(P,la))
    end

    #Clear that memory
    ldli = nothing
    la = nothing
    buffs = nothing
    S_xx = nothing

    #Pre-allocate Rademacher Seeds
    mts = MersenneTwister.(1:Threads.nthreads())

    #Initialize output
    Pii=zeros(NT)
    Bii_second=zeros(NT)
    Bii_cov= settings.cov_effects ==true ? zeros(NT) : nothing
    Bii_first= settings.first_id_effects == true ? zeros(NT) : nothing

    Mii = zeros(NT)
    Pii_sq = zeros(NT)
    Mii_sq = zeros(NT)
    Pii_Mii = zeros(NT)

    Threads.@threads for i=1:p

        rademach =  rand(mts[Threads.threadid()],1,NT) .> 0.5
        rademach  = rademach  .- (rademach .== 0)
        rademach  = rademach / sqrt(p)

        Z= compute_sol[Threads.threadid()]( [rademach*X...] ; verbose=false)
        Z = X*Z

        #Auxiliaries for Non-linear Correction

        aux = Z.^2 / p
        Pii = Pii .+ aux

        aux = Z.^4 / p
        Pii_sq = Pii_sq .+ aux

        aux			= ((rademach' - Z).^2)/p
		Mii			= Mii .+ aux
		aux			= ((rademach' - Z).^4)/p
		Mii_sq		= Mii_sq .+ aux

        Pii_Mii		= Pii_Mii .+ ((Z.^2).*((rademach' - Z).^2) )/p

        #Demeaned Rademacher
        rademach =  rand(mts[Threads.threadid()],1,size(Fvar,1)) .> 0.5
        rademach  = rademach  .- (rademach .== 0)
        rademach  = rademach /sqrt(p)
        rademach = rademach .- mean(rademach)

        Z = compute_sol[Threads.threadid()]( [rademach*Fvar...] ; verbose=false)
        Z = X*Z
        Bii_second = Bii_second .+ (Z.*Z)

        if settings.first_id_effects == true | settings.cov_effects == true
            Z_pe = compute_sol[Threads.threadid()]( [rademach*Dvar...] ; verbose=false)
            Z_pe = X*Z_pe

            if settings.first_id_effects == true
                Bii_first = Bii_first .+ (Z_pe.*Z_pe)
            end

            if settings.cov_effects == true
                Bii_cov = Bii_cov .+ (Z_pe.*Z)
            end

        end
    end

    #Censor
    #Pii[ findall(Pii.>=0.99)] .= 0.99

    #Account for Non-linear Bias
    Pii = Pii ./ (Pii + Mii)
    Mii = 1 .- Pii
    Vi = (1/p)*((Mii.^2).*Pii_sq+(Pii.^2).*Mii_sq-2*Mii.*Pii.*Pii_Mii)
    Bi = (1/p)*(Mii.*Pii_sq-Pii.*Mii_sq+2*(Mii-Pii).*Pii_Mii)
    correction_JLA = (1 .- Vi./(Mii.^2)+Bi./Mii)

    return (Pii = Pii , Mii = Mii , correction_JLA = correction_JLA, Bii_first = Bii_first , Bii_second = Bii_second , Bii_cov = Bii_cov)
end

"""
$(SIGNATURES)

This function regresses fixed effects based onto some observables. See appendix in KSS for more information.

### Arguments
* `y`: outcome variable.
* `X`: the design matrix in the linear model.
* `Z`: matrix of observables to use in regression.
* `Transform`: matrix to compute fixed effects (e.g. Transform = [0 F] recovers second fixed effects).
* `sigma_i`: estimate of the unbiased variance of observation i.
* `lincom_labels`: labels of the columns of Z.
* `settings`: settings based on data type `VCHDFESettings`. Please see the reference provided below.
"""
function lincom_KSS(y,X, Z, Transform, sigma_i; lincom_labels = nothing)
    #lincom_KSS(y,X, Z, Transform, sigma_i, labels; joint_test =false, joint_test_regressors = nothing, nsim = 10000)
    #regressors is a vector of strings
    #if joint_test_regressors != nothing
    #    positionvec = indexin(joint_test_regressors, labels) .+ 1
    #end

    #PART 1: COMPUTE FE
    n = size(y,1)
    r = size(Z,2)

    xx=X'*X
    xy=X'*y
    compute_sol = approxcholSolver(xx;verbose = false)
    beta = compute_sol([xy...];verbose=false)

    #PART 2: SET UP MATRIX FOR SANDWICH FORMULA
    wy = Transform*beta
    zz = Z'*Z
    numerator=Z\wy
    sigma_i  = sparse(collect(1:n),collect(1:n),[sigma_i ...],n,n)
    sigma_i_naive = sparse(collect(1:n),collect(1:n),[(y-X*beta).^2 ...],n,n)

    #PART 3: COMPUTE STATS
    denominator  = zeros(r,1)
    denominator_naive   = zeros(r,1)

    for q=1:r
        v=sparse([q],[1.0],[1.0],r,1)
        v=zz\[v...]
        v=Z*v
        v=Transform'*v

        right = compute_sol(v;verbose=false)

        left=right'

        denominator[q]=left*(X'*sigma_i*X)*right
        denominator_naive[q]=left*(X'*sigma_i_naive*X)*right
    end
    test_statistic=numerator./(sqrt.(denominator))

    #PART 4: REPORT
    println("Inference on Linear Combinations:")
    if lincom_labels == nothing
        for q=2:r
            if q <= r
                println("\nCoefficient of Column ", q-1,": ", numerator[q] )
                println("Traditional HC Standard Error of Column ", q-1,": ", sqrt(denominator_naive[q]) )
                println("KSS Standard Error of Column ", q-1,": ", sqrt(denominator[q]) )
                println("T-statistic of Column ", q-1, ": ", test_statistic[q])
            end
        end
    else
        for q=2:r
            tell_me = lincom_labels[q-1]
            println("\nCoefficient on ", tell_me,": ", numerator[q] )
            println("Traditional HC Standard Error on ", tell_me,": ", sqrt(denominator_naive[q]) )
            println("KSS Standard Error on ", tell_me,": ", sqrt(denominator[q]) )
            println("T-statistic on ", tell_me,": ", test_statistic[q])
        end
    end

    test_statistic = test_statistic[2:end]
    linear_combination = numerator[2:end]
    SE_linear_combination_KSS = sqrt.(denominator[2:end])
    SE_naive = sqrt.(denominator[2:end])
    return (test_statistic = test_statistic , linear_combination = linear_combination, SE_linear_combination_KSS = SE_linear_combination_KSS, SE_naive =SE_naive)
end

#    # PART 5: Joint-test. Quadratic form beta'*A*beta  MAYBE FOR FUTURE VERSION!
#    if joint_test ==false &  (joint_test_regressors != nothing )
#        println("Joint test will not be computed. You need to set the argument option joint_test to true.")
#    end
#
#    if joint_test == true
#
#        if joint_test_regressors  == nothing
#            restrict=sparse(collect(1:r-1),collect(2:r),1.0,r-1,r)
#        elseif joint_test_regressors  != nothing
#            restrict = sparse( collect(1:length(positionvec)), positionvec, 1.0, length(positionvec), r)
#        end
#
#        v=restrict*(zz\(Z'*Transform))
#        v=v'
#        #v=sparse(v) #ldiv doesn't work for sparse RHS
#        r=size(v,2)
#
#        #Auxiliary
#        aux=xx\v[:,:]
#        opt_weight=v'*aux
#        opt_weight=opt_weight^(-1)
#        opt_weight=(1/r)*(opt_weight+opt_weight')/2
#
#        #Eigenvalues, eigenvectors, and relevant components
#        lambda , Qtilde = eigs( v*opt_weight*v', xx; nev=r,ritzvec=true)
#        lambda = Real.(lambda)
#        Qtilde = Real.(Qtilde)
#        #lambdaS, QtildeS = eigs(v*opt_weight*v', xx; nev=1,which=:SM,ritzvec=true)
#        #lambdaS = [lambdaL; lambdaS]
#        #Qtilde = hcat(QtildeL,QtildeS)
#
#        W=X*Qtilde
#        V_b=W'*sigma_i*W
#        V_b = (1/2)*(V_b + V_b')
#
#        #Now focus on obtaining matrix Lambda_B with the A test associated with a joint hypothesis testing.
#        Bii=opt_weight^(0.5)*aux';
#        Bii=Bii*X'
#        Bii=Bii'
#        Bii = 0.5*(Bii[rows,:].*Bii[columns,:] + Bii[columns,:].*Bii[rows,:])
#        Bii = sum(Bii,dims=2)[:]
#        Lambda_B=sparse(rows,columns,Bii,n,n)
#
#        #Leave Out Joint-Statistic
#        stat=(v'*β)'*opt_weight*(v'*β)-y'*Lambda_B*eta_h
#
#        #Now simulate critical values under the null.
#        mu=zeros(r)
#        sigma = V_b
#        b_sim = MvNormal(mu,sigma)
#        b_sim = rand(b_sim, nsim)
#
#        #theta_star_sim=sum(lambda'.*(b_sim.^2 - diag(V_b)'),2)
#        theta_star_sim = sum(lambda.*b_sim.^2 .- lambda.* diag(V_b),dims=1)
#        pvalue=mean(theta_star_sim.>stat)
#
#        #Report
#        println("Joint-Test Statistic: ", stat)
#        println("p-value: ", pvalue)
#
#    end

ext: jl


text:
# Code for Video 2.3:
# The "iris" dataset contains 150 measurements on 3 species
# of Iris flowers, namely I. setosa, I. versicolor and I. virginica.
# The columns 1 to 4 respectively contain the:
# sepal length - sepal width -petal length - petal width.
# The species is specified in the last column.
# The dataset can be found in the file iris.txt

module IrisTools

### exports:
export Iris
export species
export read_file, counting_species, length_matrix, calculations_setosa, to_string

### module variables:
###################################################
fname = "iris.txt"
species = ["I. setosa", "I. versicolor", "I. virginica"]

irism = [] # will contain the measurements of the dataset, each array element is a line containing 4 measurements
iris = [] # temporary array  to contain the names of iris specimens
lengths = [] # matrix with the length measurements
species_freq = Dict{ASCIIString, Int64}()

### types:
###################################################
type Iris
    sepal_length::Float64
    sepal_width::Float64
    petal_length::Float64
    petal_width::Float64
    species::ASCIIString
end

type Setosa
    # some fields
end

type Versicolor
    # some fields
end

type Virginica
    # some fields
end

IrisU = Union{Setosa, Versicolor, Virginica}

### functions:
###################################################
"executed when the module is loaded"
function __init__()
    read_file()
end

"reads in the file iris.txt with iris measurements"
function read_file()
    open(fname) do file
        global irism = readlines(file)
    end
    # println(irism)
end

"gives a count of measurements for each species"
function counting_species()
    # gather names of species in iris array:
    for line in irism
        # println(line)
        # the species name starts from column 17:
        species = line[17:end]
        # removes the trailing \r\n:
        species = strip(species)
        push!(iris, species)
    end
    # count species in dictionary species-freq:
    for species in iris
        haskey(species_freq, species) ?
            species_freq[species] += 1 :
            species_freq[species] = 1
    end
    # sort the the keys and print out the frequencies:
    println("Species : frequency \n")
    species = sort!(collect(keys(species_freq)))
    for sp in species
        println("$sp : $(species_freq[sp])")
    end
end

" extracts a matrix with the length measurements"
function length_matrix()
    # extract the measurements as numbers in a matrix lengths:
    lengths_temp = Array(Float64, 1, 4)
    for line in irism
        # The columns 1 to 16 respectively contain the lengths of 1 measurement:
        mlengths = line[1:16] # measurement lengths
        arr_temp = split(mlengths) # split on spaces to an array of strings
        # convert to an array of Float64 values (a column vector 4 x 1)
        arr_mlength = [parse(Float64, str) for str in arr_temp]
        # transpose arr_mlength to obtain a 1 x 4 array
        # and append it vertically to lengths
        lengths_temp = vcat(lengths_temp, arr_mlength')
        # we must indicate that we mean the variable lengths
        # defined in global scope by prepending it with global
    end
    # cut off the first row which contained arbitrary values
    global lengths = lengths_temp[2:end, :]
end

" performs some calculations on I. setosa measurements"
function calculations_setosa()
    # extract I. setosa measurements:
    lengths_setosa = lengths[1:50, :]
    # do some calculations:
    println("sum of each column")
    println(sum(lengths_setosa, 1))
    println("mean of the features:")
    println(mean(lengths_setosa, 1))
    println("smallest values of the features:")
    println(sort(lengths_setosa, 1)[1, :])
    println("largest values of the features:")
    println(sort(lengths_setosa, 1)[end, :])
end

function show(i::Iris)
    println("This $(i.species) has measurements: ")
    println("Sepal length: $(i.sepal_length)")
    println("Sepal width: $(i.sepal_width)")
    println("Petal length: $(i.petal_length)")
    println("Petal width: $(i.petal_width)")
end

function to_string(i::Iris)
    "This $(i.species) has measurements: \n" *
    "\tSepal length: $(i.sepal_length)\n" *
    "\tSepal width: $(i.sepal_width)\n" *
    "\tPetal length: $(i.petal_length)\n" *
    "\tPetal width: $(i.petal_width)\n"
end

end  # module end

max_stars_repo_path: Section 2/iris_tools.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 2/iris_tools.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 2/iris_tools.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
using Plots
using Optim

# Calculate dose-response
# flux reaction:
function reaction(conc,Km,Vmax)
    10 .^Vmax.*conc./(10 .^Km .+ conc)
end

Km = 0
Vmax = 1
conc = collect(0:.01:2)

flux = reaction(conc, Km, Vmax)
plot(conc,flux)

# simulate data with 7% std
conc_sim = [.1,.5,1,2]
flux_sim = reaction(conc_sim,  Km, Vmax)
flux_sim = randn(length(flux_sim))*0.07.*flux_sim+flux_sim
std_sim = flux_sim * .07

scatter!(conc_sim,flux_sim,yerror=std_sim,grid=false,linealpha = 0)

# negative log-likelihood for simulated data
function L(conc,flux,std,p)
    X = reaction(conc,p[1],p[2])
    LL = -sum(-length(X)*log(2π)/2 .- length(X).*log.(abs.(std)) .- (sum((flux .- X).^2) ./ (2std.^2)))
end

# maximum likelihood fit
init = [1.0,1.0]        # start parameters Km,Vmax
optimum = optimize(p -> L(conc_sim,flux_sim,std_sim,p),init,Newton())

### We will continue with the profile likelihood on July, 10th


text:
using KernelDensity, AbstractPlotting, Makie
import AbstractPlotting: convert_arguments, PlotFunc, plottype
to_tuple(t::Tuple) = t
to_tuple(t) = (t,)

function convert_arguments(P::PlotFunc, f::Function, args...; kwargs...)
    tmp = f(args...; Iterators.filter(t -> last(t) != automatic, kwargs)...) |> to_tuple
    convert_arguments(P, tmp...)
end

# remove convert_arguments(P, f, x) = (x, f.(x))
function convert_arguments(P::PlotFunc, d::KernelDensity.UnivariateKDE)
    ptype = plottype(P, Lines) # choose the more concrete one
    ptype => convert_arguments(ptype, d.x, d.density)
end

function convert_arguments(P::Type{<: Combined{T}}, d::KernelDensity.BivariateKDE) where T
    ptype = plottype(P, Heatmap)
    ptype => convert_arguments(ptype, d.x, d.y, d.density)
end

plot(kde, rand(100)) #line plot
scatter(kde, rand(100), markersize = 0.02, color = (:black, 0.1))

plot(kde, rand(100, 2)) #heatmap
surface(kde, rand(100, 2))


text:
# Submit coverage to Coveralls.io
export Coveralls
"""
Coverage.Coveralls Module

This module provides functionality to push coverage information to the Coveralls
web service. It exports the `submit` and `submit_local` methods.
"""
module Coveralls
    using Coverage
    using CoverageTools
    using HTTP
    using JSON
    using LibGit2
    using MbedTLS

    export submit, submit_local

    #=
    JSON structure for Coveralls
    Accessed 2015/07/24:
    https://coveralls.zendesk.com/hc/en-us/articles/201774865-API-Introduction
    {
      "service_job_id": "1234567890",
      "service_name": "travis-ci",
      "source_files": [
        {
          "name": "example.rb",
          "source": "def four\n  4\nend",
          "coverage": [null, 1, null]
        },
        {
          "name": "lib/two.rb",
          "source": "def seven\n  eight\n  nine\nend",
          "coverage": [null, 1, 0, null]
        }
      ]
    }
    =#

    # to_json
    # Convert a FileCoverage instance to its Coveralls JSON representation
    function to_json(fc::FileCoverage)
        name = Sys.iswindows() ? replace(fc.filename, '\\' => '/') : fc.filename
        return Dict("name"          => name,
                    "source_digest" => digest(MD_MD5, fc.source, "secret"),
                    "coverage"      => fc.coverage)
    end

    # Format the body argument to HTTP.post
    function makebody(data::Dict)
        iodata = IOBuffer(JSON.json(data))
        json_file = HTTP.Multipart("json_file", iodata, "application/json")
        return Dict("json_file" => json_file)
    end

    """
        submit(fcs::Vector{FileCoverage}; kwargs...)

    Take a vector of file coverage results (produced by `process_folder`),
    and submits them to Coveralls. Assumes that this code is being run
    on TravisCI, AppVeyor or Jenkins. If running locally, use `submit_local`.
    """
    function submit(fcs::Vector{FileCoverage}; kwargs...)
        data = prepare_request(fcs, false)
        post_request(data)
    end

    function prepare_request(fcs::Vector{FileCoverage}, local_env::Bool, git_info=query_git_info)
        data = Dict{String,Any}("source_files" => map(to_json, fcs))

        if local_env
            # Attempt to parse git info via git_info, unless the user explicitly disables it by setting git_info to nothing
            data["service_name"] = "local"
            data["git"] = parse_git_info(git_info)
        elseif lowercase(get(ENV, "APPVEYOR", "false")) == "true"
            data["service_job_number"] = ENV["APPVEYOR_BUILD_NUMBER"]
            data["service_job_id"] = ENV["APPVEYOR_BUILD_ID"]
            data["service_name"] = "appveyor"
            appveyor_pr = get(ENV, "APPVEYOR_PULL_REQUEST_NUMBER", "")
            isempty(appveyor_pr) || (data["service_pull_request"] = appveyor_pr)
        elseif lowercase(get(ENV, "TRAVIS", "false")) == "true"
            data["service_number"] = ENV["TRAVIS_BUILD_NUMBER"]
            data["service_job_id"] = ENV["TRAVIS_JOB_ID"]
            data["service_name"] = "travis-ci"
            travis_pr = get(ENV, "TRAVIS_PULL_REQUEST", "")
            isempty(travis_pr) || (data["service_pull_request"] = travis_pr)
        elseif lowercase(get(ENV, "JENKINS", "false")) == "true"
            data["service_job_id"] = ENV["BUILD_ID"]
            data["service_name"] = "jenkins-ci"
            data["git"] = parse_git_info(git_info)

            # get the name of the branch if not a pull request
            if get(ENV, "CI_PULL_REQUEST", "false") == "false"
                data["git"]["branch"] = split(ENV["GIT_BRANCH"], "/")[2]
            end
        elseif haskey(ENV, "GITHUB_ACTION")
            data["service_job_id"] = ENV["GITHUB_RUN_ID"]
            data["service_name"] = "github"
            data["git"] = parse_git_info(git_info)

            event_path = open(JSON.Parser.parse, ENV["GITHUB_EVENT_PATH"])
            github_pr_info = get(event_path, "pull_request", Dict())
            github_pr = get(github_pr_info, "number", "")
            isempty(github_pr) || (data["service_pull_request"] = github_pr)
        else
            data["git"] = parse_git_info(git_info)
        end

        service_name = get(ENV, "COVERALLS_SERVICE_NAME", "")
        isempty(service_name) || (data["service_name"] = service_name)

        service_number = get(ENV, "COVERALLS_SERVICE_NUMBER", "")
        isempty(service_number) || (data["service_number"] = service_number)

        service_job_number = get(ENV, "COVERALLS_SERVICE_JOB_NUMBER", "")
        isempty(service_job_number) || (data["service_job_number"] = service_job_number)

        jobid = get(ENV, "COVERALLS_SERVICE_JOB_ID", "")
        isempty(jobid) || (data["service_job_id"] = jobid)

        flag_name = get(ENV, "COVERALLS_FLAG_NAME", "")
        isempty(flag_name) || (data["flag_name"] = flag_name)

        ci_pr = get(ENV, "COVERALLS_PULL_REQUEST", "")
        isempty(ci_pr) || (data["service_pull_request"] = ci_pr)

        if !haskey(data, "service_name")
            error("No compatible CI platform detected")
        end

        data = add_repo_token(data, local_env)
        if get(ENV, "COVERALLS_PARALLEL", "false") == "true"
            data["parallel"] = "true"
        end
        return data
    end

    function parse_git_info(git_info::Function)
        result = nothing
        try
            result = git_info()
        catch ex
            @warn "Parse of git information failed" exception=e, catch_backtrace()
        end
        return result
    end

    parse_git_info(git_info::Dict) = git_info

    # query_git_info
    # Pulls information about the repository that isn't available if we
    # are running somewhere other than TravisCI
    function query_git_info(dir=pwd())
        repo            = LibGit2.GitRepoExt(dir)
        head            = LibGit2.head(repo)
        head_cmt        = LibGit2.peel(head)
        head_oid        = LibGit2.GitHash(head_cmt)
        commit_sha      = string(head_oid)
        author_name     = string(LibGit2.author(head_cmt).name)
        author_email    = string(LibGit2.author(head_cmt).email)
        committer_name  = string(LibGit2.committer(head_cmt).name)
        committer_email = string(LibGit2.committer(head_cmt).email)
        message         = LibGit2.message(head_cmt)
        remote_name     = "origin"
        branch          = LibGit2.shortname(head)

        # determine remote url, but only if repo is not in detached state
        remote = ""
        if branch != "HEAD"
            LibGit2.with(LibGit2.get(LibGit2.GitRemote, repo, remote_name)) do rmt
                remote = LibGit2.url(rmt)
            end
        end
        LibGit2.close(repo)

        return Dict(
            "branch"    => branch,
            "remotes"   => [
                Dict(
                    "name"  => remote_name,
                    "url"   => remote
                )
            ],
            "head" => Dict(
                "id" => commit_sha,
                "author_name"       => author_name,
                "author_email"      => author_email,
                "committer_name"    => committer_name,
                "committer_email"   => committer_email,
                "message"           => message
            )
        )
    end

    """
        submit_local(fcs::Vector{FileCoverage}, git_info=query_git_info; kwargs...)

    Take a `Vector` of file coverage results (produced by `process_folder`),
    and submits them to Coveralls. For submissions not from CI.

    git_info can be either a `Dict` or a function that returns a `Dict`.
    """
    function submit_local(fcs::Vector{FileCoverage}, git_info=query_git_info; kwargs...)
        data = prepare_request(fcs, true, git_info)
        post_request(data)
    end

    # posts the actual request given the data
    function post_request(data)
        @info "Submitting data to Coveralls..."
        coveralls_url = get(ENV, "COVERALLS_URL", "https://coveralls.io/api/v1/jobs")
        req = HTTP.post(coveralls_url, HTTP.Form(makebody(data)))
        @debug "Result of submission:\n" * String(req)
        nothing
    end

    # adds the repo token to the data
    function add_repo_token(data, local_submission)
        repo_token =
                get(ENV, "COVERALLS_TOKEN") do
                    get(ENV, "REPO_TOKEN") do #backward compatibility
                        # error unless we are on Travis
                        if local_submission || (data["service_name"] != "travis-ci")
                            error("Coveralls submission requires a COVERALLS_TOKEN environment variable")
                        end
                    end
                end
        if repo_token !== nothing
            data["repo_token"] = repo_token
        end
        return data
    end

end  # module Coveralls


text:
# version 0.1

module EcoEvo3D

using StatsBase

function getSampleFromArray(items)
  return sample(items)
end

function readLocations(filename)
	dat = readdlm(filename,'\t');
	@inbounds w = hcat(dat[:,2],dat[:,3]);
	w;
end

function readDistanceMatrix(filename)
	dat = readdlm(filename,' ');
	dat;
end

function changediagonal!(MA,nrows,value)
	@inbounds for(i in 1:nrows)
		@inbounds MA[i,i]=value;
	end
end

function createMRM(MRM,Sti,Sp,ts)
	MRM = [Sti Sp 1 ts];
	MRM;
end

function createMC(MC,Sti,Sp,ts)
	MC = [Sti Sp 1 ts];
	return MC;
end

function createMA(MA,Sti,Stj,Sp,anaG)
	MA = [Sti Stj Sp anaG];
	return MA;
end

function updateMC(MC)
	if length(MC)>0
		MC[:,3] = MC[:,3].+1;#event: -1 for extinction; 1 for speciation
	end
	return MC;
end

function updateMA(MA)
	if length(MA)>0
		MA[:,4] = MA[:,4] .- 1;
	end
	return MA;
end

function checkIfThereIsMRM(MRM,Sti,Sp,ts)
	if length(MRM)==0
		MRM = createMRM(MRM,Sti,Sp,ts);
	else
		pos = find( (MRM[:,1].==Sti) & (MRM[:,2].==Sp))
		if length(pos)==0
			MRM = cat(1,MRM,[Sti Sp 1 ts]);
		end
	end
	return MRM;
end

function checkIfThereIsMC(MC,Sti,Sp,ts)
	if length(MC)==0
		MC = createMC(MC,Sti,Sp,ts);
	else
		pos = find( (MC[:,1].==Sti) & (MC[:,2].==Sp))#position in the matrix MA referred to the presence of individuals of species 'Sp' coming from site 'Stj' to site 'Sti'
		if length(pos)==0
			MC = cat(1,MC,[Sti Sp 1 ts]);
		end
	end
	return MC;
end

function GetGammaRichness(R,S)
	richnessspeciesR = [];
	#%gamma richness
	for (i in 1:S)
		AR = sort(R[i])
		richnessspeciesR = [richnessspeciesR; unique(AR)];
	end;
	extantR = length(unique(sort(richnessspeciesR)));
	gamma = extantR;

	return gamma;
end

function OutputPerGeneration(outputfilepergen,ri,cost,J,G,S,k,anaG,retG,mr,ml,v,lakesIsolation,lakesArea,gamma,alpharich,SpecANA,SpecCLA,SpecMR,DispersalRich)
	for i in 1:S
		writedlm(outputfilepergen, [ri cost J G k anaG retG mr ml v gamma i lakesIsolation[i] lakesArea[i] alpharich[i] SpecANA[i] SpecCLA[i] SpecMR[i] DispersalRich[i]],' ');
	end
	flush(outputfilepergen);#To print in the output file each realization
	return;
end

function printPhylogeny(phylogenyfile,old,new,ts,ri)
	writedlm(phylogenyfile,[ri old new ts],' ');
        flush(phylogenyfile);
end

function UAB(MA,MC)
	MC = updateMC(MC);
	MA = updateMA(MA);
	return MC,MA;
end

function UAC(MA)
	MA = updateMA(MA);
	return MA;
end

function UARM(MA,MC)
	MC = updateMC(MC);
	MA = updateMA(MA);
	return MC,MA;
end

function UALM(MA,MC,R,Sti,Stj,Sp,anaG,retG,ts)
	MC = updateMC(MC);
	MA = updateMA(MA);

	if length(MA)==0 #Si no hay MA
  	MA = createMA(MA,Sti,Stj,Sp,anaG);#Create MA
	else#Si hay MA
		pos = find( (MA[:,1].==Sti) & (MA[:,2].==Stj) & (MA[:,3].==Sp))#position in the matrix MA referred to the presence of individuals of species 'Sp' coming from site 'Stj' to site 'Sti'
		if length(pos)==0 #No hay la linea
			indalive = length(find(R[Sti].==Sp))#Hay individuos de la specie sp vivos en el sitio Sti
			if (indalive == 0) #Checking if there are individuals of species 'Sp' alive in site 'Sti'
#        			println("<<<<<<<<<<< A D D    R O W   T O    A N A G E N E S I S   <<<<<<<<<<<<<<<<<");
				MA = cat(1,MA,[Sti Stj Sp anaG]);#Create the row in MA matrix, starting by 1
      			end
		else #Ya hay la linea
#			println("<<<<<<<<<<< R E T A R D    A N A G E N E S I S   <<<<<<<<<<<<<<<<<");
			MA[pos,4] = MA[pos,4] .+ retG;#retards anagenesis by increasing the remaining
		end
	end
  	return MC,MA,R;
end

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

function checkAna(MA,R,anaG,lastspecies,listofanagenesis,ts,phylogenyfile,ri)
	pos = [];
	if length(MA)>0
		pos = find(MA[:,4] .== 0)#find the anagenesis events
		if length(pos)>0#is there was an anagenesis event
		  @inbounds for (p in pos)
  			Sti=round(Integer,MA[p,1]);#pos represents the rows of the matrix. p is one row. MA[p,1] is the target site
  			Stj=round(Integer,MA[p,2]);#pos represents the rows of the matrix. p is one row. MA[p,2] is the source site
		  	Sp=round(Integer,MA[p,3]);#pos represents the rows of the matrix. p is one row. MA[p,3] is the ancient species
			  MA, R, lastspecies,listofanagenesis = AnagenesisSpeciation(MA,R,Sti,Stj,Sp,lastspecies,listofanagenesis,ts,phylogenyfile,ri);#speciation in target site
#			  MA, R, lastspecies,listofanagenesis = MA, R, lastspecies, listofanagenesis;
		  end
	  end
  end
  return MA,R,lastspecies,listofanagenesis;
end

function AnagenesisSpeciation(MA,R,Sti,Stj,Sp,lastspecies,listofanagenesis,ts,phylogenyfile,ri)
	newspeciesAna = lastspecies + 1;#the id of the new species
	ancientindividuals = find( R[Sti].==Sp )#the position of all the individuals of the ancient species 'Sp' in the target site
	R[Sti][ancientindividuals] = newspeciesAna;#the speciation itself: all the individuals of former species 'Sp' in the target site are now from a new species 'newspeciesAna'
        println("new ANAgenesis Speciation : Species ",newspeciesAna," - time: ",ts);
 	pos = find( (MA[:,1].==Sti) & (MA[:,3].==Sp))#position in the matrix MA referred to the presence of individuals of species 'Sp' in site 'Sti'
	MA = MA[1:size(MA,1).!=pos,:];#Borra la linea 'pos' de la matriz MA!!
	printPhylogeny(phylogenyfile,Sp,newspeciesAna,ts,ri);

  	if length(listofanagenesis)>0
		listofanagenesis = cat(1,listofanagenesis,[Sti newspeciesAna]);
	else
		listofanagenesis = cat(1,[Sti newspeciesAna]);
	end

	return MA,R,newspeciesAna,listofanagenesis;
end

function CladogenesisEvent(MC,R,Sti,Individual,lastspecies,ts,phylogenyfile,ri)
	newspeciesClado = lastspecies + 1;
	printPhylogeny(phylogenyfile,R[Sti][Individual],newspeciesClado,ts,ri);
        println("new CLADOgenesis Speciation : Species ",newspeciesClado," - time: ",ts);
   	R[Sti][Individual] = newspeciesClado;
	MC = checkIfThereIsMC(MC,Sti,newspeciesClado,ts);
	return MC,R,newspeciesClado;
end

function LocalMigrationEvent(R,KillHab,MigrantHab,KillInd,Dc,Ji,S)
	MigrantSpecies = -1;

	allpos = find(Dc[KillHab,:] .>= MigrantHab);#All the sites at a distance lower than the threshold 'MigrantHab'
	kr = minimum(allpos[find(allpos .!= KillHab)]);
	MigrantInd = rand(1:Ji[kr]);
	MigrantSpecies = R[kr][MigrantInd];
#        println("LOCAL Migration");
	R[KillHab][KillInd] = R[kr][MigrantInd];
	return kr,R,MigrantSpecies;
end;

function RegionalMigrationEvent(MRM,R,Sti,Individual,ts,lastspecies)
  newspeciesMR = lastspecies + 1;
  R[Sti][Individual] = newspeciesMR;
  println("new RegionalMig Speciation: Species ",newspeciesMR," - time: ",ts)
  MRM = checkIfThereIsMRM(MRM,Sti,newspeciesMR,ts);

  return MRM,R,newspeciesMR;
end;

function BirthEvent(R,BirthLocal,KillInd,KillHab)
	R[KillHab][KillInd] = R[KillHab][BirthLocal];
	return R;
end

function calculateSpeciationMA(listofanagenesis,R)
	S = length(R);#number of sites equals the number of rows in R
	speciatedMA = zeros(S);
	if (length(listofanagenesis) > 0)
		@inbounds for i in 1:S
			@inbounds speciesR = unique(sort(R[i]))';
			pos = find(listofanagenesis[:,1].== i);
			if length(pos)>0
				speciesMA = listofanagenesis[pos,2];
				speciatedMA[i] = length(setdiff(speciesMA,setdiff(speciesMA,speciesR))');
			end
		end
	end
	return speciatedMA;
end

function calculateSpeciationMC(MC,R,S,k,Ji)
	speciatedMC = zeros(S);
	if (length(MC) > 0)
		@inbounds for i in 1:S
		@inbounds speciesR = unique(sort(R[i]))';
			pos = find( (MC[:,1].== i) & (MC[:,3].>= k) );#Only the species that are in the system for more than k time steps (Protracted)
			if length(pos)>0
				speciesMC = unique(sort(MC[pos,2]));
				speciatedMC[i] = length(setdiff(speciesMC,setdiff(speciesMC,speciesR))');
			end
		end
	end
	return speciatedMC;
end

function calculateSpeciationMR(MRM,R,S,Ji)
	speciatedMRM = zeros(S);
	if (length(MRM) > 0)
		@inbounds for i in 1:S
      @inbounds speciesR = unique(sort(R[i]))';
      pos = find(MRM[:,1].== i);#check for site i
      if length(pos)>0
        speciesMRM = unique(sort(MRM[pos,2]));
        @inbounds speciatedMRM[i] = length(setdiff(speciesMRM,setdiff(speciesMRM,speciesR))');
      end
		end
	end
#	writedlm("MRMmatrix.dat",MRM,' ');
	return speciatedMRM;
end

function richnessanalysis!(S,R,Ji,richnessspeciesR,alpharich)
	#%gamma richness
	@inbounds for (i in 1:S)
		@inbounds AR = sort(R[i]);
		richnessspeciesR = [richnessspeciesR; unique(AR)];
		@inbounds alpharich[i] = length(unique(AR));
	end;
	gamma = length(unique(sort(richnessspeciesR)));
	return gamma,alpharich;
end

function normalizeDc!(Dc,S)
  for(i in 1:S)
    Dc[i,:] = Dc[i,:]/maximum(Dc[i,:]);
  end
end

function repeatval(x,N)
  return map(v -> v=x,1:N);
end

#Each site starts with one different species
function initialPopulation(R,Ji)
  lastspecies = 0;
  for(iJ in Ji)#iterating over columns
     lastspecies = lastspecies+1;
     popi = repeatval(lastspecies,iJ);
     push!(R,popi);
  end
  return R,lastspecies;
end

function dynamic(seed,nreal,Gmax,J,v,mr,ml,anaG,retG,distmatfile,verticesdata,model)
	DI = readDistanceMatrix(distmatfile);#the location of the points of the landscape.
	Dc = cumsum(DI,2);#Should we divide
  dT = sum(DI,2);
	const S = length(DI[:,1]);#Number of sites
  normalizeDc!(Dc,S);
#### To get the cost
        substring = distmatfile[search(distmatfile,'_')+1:end];
        modelstop = search(substring,'_');
        coststop = search(substring,'t');
	const cost = float(substring[coststop+1:modelstop-1]);

	sitesdata = readdlm(verticesdata,' ',header=true)[1];#Proportion of individuals of each site
#	Pj = ones(S);#Proportion of individuals of each site - if an input file is not defined, the proportion is the same for each site
        lakesVolume = sitesdata[:,2];
        lakesArea = sitesdata[:,3];
	Pj = lakesVolume;#In case we define different carrying capacities for each site we consider the volume of the lakes (third column)
	mins = find(sitesdata.==minimum(sitesdata[:,1]));#the first column represents the height of the lakes
	entrypoint = mins[rand(1:length(mins))];
	t=1;#Sites have different sizes and are located at different height.
	Ji=round(Integer,J * Pj/sum(Pj));

	outputfilepergen = open(string("RichnessPerGen_AnaG_",anaG,"_retG_",retG,"_cost_",cost,"_MR_",signif(mr,3),"_VR_",signif(v,3),".txt"),"a")
	writedlm(outputfilepergen, ["Real Cost J G Gi anaG retG mr ml v gamma Site dT lakeArea alpharich SpecANA SpecCLA SpecMR DispersalRich"]);

	outputfile = open(string("RichnessPerSite_AnaG_",anaG,"_retG_",retG,"_cost_",cost,"_MR_",signif(mr,3),"_VR_",signif(v,3),".txt"),"a")
	writedlm(outputfile,["Real Cost Model J G anaG retG Site lakeArea Ji dT mr ml v gamma alpharich SpecANA SpecCLA SpecMR DispersalRich"]);

	phylogenyfile = open(string("Phylogeny_AnaG_",anaG,"_retG_",retG,"_cost_",cost,"_MR_",signif(mr,3),"_VR_",signif(v,3),".txt"),"a")
	writedlm(phylogenyfile,["Repl Ancestral Derived Age"]);

  @inbounds for (ri in 1:nreal)#realizations
 	  lastspecies = 0;
  	ts=0;
    R=[];
    R,lastspecies = initialPopulation(R,Ji)
	  MA = Array(Number,0,0);#Matrix to calculate Anagenesis speciation (MA = [Sti, Stj, S, C])
	  MC = Array(Number,0,0);#Matrix to control Cladogenesis speciation (MC = [Sti, S, E])
	  MRM = Array(Number,0,0);#Matrix to control events of Regional Migration (MRM = [Sti, S, ts])
		srand(seed+(7*ri));
		listofanagenesis = Array(Number,0,0);

		#%Resources
#		G = rand(15:Gmax);#Minimum of 15 generations
		G = Gmax;#Minimum of 15 generations

		@inbounds for (k = 1:G)#%population-metapopulation-metacommunity dynamics (not-tracking multitrophic metacommunity dynamics!)
			@inbounds for (j = 1:J)#For each individual in each site
				ts = ts+1;#Time step increases
				mvb = rand();
		       		#Demography - Resources
    				KillHab = rand(1:S);#which site to kill
				KillInd = rand(1:Ji[KillHab]);#which individual to kill
        			KillIndEntryPoint = rand(1:Ji[entrypoint]);
				MigrantHab = rand()*maximum(Dc[KillHab,:]);
				BirthLocal = getSampleFromArray(1:length(R[KillHab]));#which individual to born

				if mvb <= ml;#Local Migration event
					kr,R,MigrantSpecies = LocalMigrationEvent(R,KillHab,MigrantHab,KillInd,Dc,Ji,S);
					MC,MA,R = UALM(MA,MC,R,KillHab,kr,MigrantSpecies,anaG,retG,ts);#Update Anagenesis after Local Migration event
				elseif (ml < mvb <= (ml+mr));#Regional Migration event
					MRM,R,lastspecies = RegionalMigrationEvent(MRM,R,entrypoint,KillIndEntryPoint,ts,lastspecies);#Speciation through the entry point
					MC,MA = UARM(MA,MC);#Update Anagenesis after Regional Migration event
				elseif ((ml+mr) < mvb <= (ml+mr+v))
					if(v > 0)#we only simulate Cladogenesis when the probability is higher than 0
						MC,R,lastspecies = CladogenesisEvent(MC,R,KillHab,KillInd,lastspecies,ts,phylogenyfile,ri);
						MA = UAC(MA);#Update Anagenesis after Cladogenesis Speciation
					end
				else #Birth event
					R = BirthEvent(R,BirthLocal,KillInd,KillHab);#Birth event
					MC,MA = UAB(MA,MC);#Update MA after Birth event
		       		end;
####ANAGENESIS
				lengthlistbefore = length(listofanagenesis);#Size of Anagenesis Matrix before checking for Anagenesis
				MA,R,lastspecies,listofanagenesis = checkAna(MA,R,anaG,lastspecies,listofanagenesis,ts,phylogenyfile,ri); #After the update of matrix MA, we check for events of Anagenesis
				lengthlistafter = length(listofanagenesis);#Size of Anagenesis Matrix after checking for Anagenesis
####CLADOGENESIS
#				if ((lengthlistafter - lengthlistbefore) > 0)#If there is at least one Anagenesis Speciation Event
#					mvc = rand();
#					if (mvc <= v);#Probability to occur a Cladogenesis Speciation event
#						allInd = find(R[St,:].==lastspecies);#All individuals from the new species
#						if (length(allInd)>0)
#							KillInd = allInd[rand(1:length(allInd))];#Randomly choose an individual from the new species
#							MC,R,lastspecies = CladogenesisEvent(MC,R,KillHab,KillInd,lastspecies,ts,phylogenyfile,ri);#Cladogenesis event
#							MA = UAC(MA);#Update Anagenesis after Cladogenesis Speciation event
#						end #if allInd
#					end #if mvc
#				end #if lengthlistafter
			end;#end S*Ji
			richnessspeciesR = [];
			alpharich = zeros(S);
			gamma,alpharich = richnessanalysis!(S,R,Ji,richnessspeciesR,alpharich);
			SpecANA = calculateSpeciationMA(listofanagenesis,R);
			SpecCLA = calculateSpeciationMC(MC,R,S,anaG,Ji);
			SpecMR = calculateSpeciationMR(MRM,R,S,Ji);
			DispersalRich = alpharich - (SpecANA + SpecCLA + SpecMR);
			OutputPerGeneration(outputfilepergen,ri,cost,J,G,S,k,anaG,retG,mr,ml,v,dT,lakesArea,gamma,alpharich,SpecANA,SpecCLA,SpecMR,DispersalRich);
		end;#end Gmax

		#To analyze the resulting richness
		richnessspeciesR = [];
		alpharich = zeros(S);
		gamma,alpharich = richnessanalysis!(S,R,Ji,richnessspeciesR,alpharich);
		SpecANA = calculateSpeciationMA(listofanagenesis,R);
		SpecCLA = calculateSpeciationMC(MC,R,S,anaG,Ji);
		SpecMR = calculateSpeciationMR(MRM,R,S,Ji);
		DispersalRich = alpharich - (SpecANA + SpecCLA + SpecMR);
#		matrixfile = open("Rmatrix.txt","a");
		for i in 1:S
			writedlm(outputfile,[ri cost model J G anaG retG i lakesArea[i] Ji[i] dT[i] mr ml v gamma alpharich[i] SpecANA[i] SpecCLA[i] SpecMR[i] DispersalRich[i]],' ');
#                        for j in 1:length(R[i])
#	                        print(matrixfile,R[i][j],' ');
#			end
#	                print(matrixfile,'\n');
		end
		flush(outputfile);
#		ANAfile = open("ANAmatrix.txt","a");
#		MRMfile = open("MRMmatrix.txt","a");
#		rowsANA = size(listofanagenesis,1);rowsMRM = size(MRM,1);
#		colsANA = size(listofanagenesis,2);colsMRM = size(MRM,2);
#		for i in 1:rowsANA #number of rows
#			for j in 1:colsANA
#	                        print(ANAfile,listofanagenesis[i,j],' ');
#			end
#	                print(ANAfile,'\n');
#		end
#		for i in 1:rowsMRM #number of rows
#			for j in 1:colsMRM
#	                        print(MRMfile,MRM[i,j],' ');
#			end
#	                print(MRMfile,'\n');
#		end

	end#%ri
	close(outputfile);
	close(outputfilepergen);
	close(phylogenyfile);
#	close(logfile);
end#ecoevo3d function

end #module


text:
import PredictMDAPI

# types

const AbstractFittable = PredictMDAPI.AbstractFittable
const AbstractEstimator = PredictMDAPI.AbstractEstimator
const AbstractPipeline = PredictMDAPI.AbstractPipeline
const AbstractTransformer = PredictMDAPI.AbstractTransformer

const AbstractFeatureContrasts = PredictMDAPI.AbstractFeatureContrasts
const AbstractNonExistentFeatureContrasts = PredictMDAPI.AbstractNonExistentFeatureContrasts

const AbstractNonExistentUnderlyingObject = PredictMDAPI.AbstractNonExistentUnderlyingObject

const AbstractBackend = PredictMDAPI.AbstractBackend

const AbstractPlot = PredictMDAPI.AbstractPlot

# traits

# const TargetStyle = PredictMDAPI.TargetStyle
# const UnknownTargetStyle = PredictMDAPI.UnknownTargetStyle
# const MixedTargetStyle = PredictMDAPI.MixedTargetStyle
# const Regression = PredictMDAPI.Regression
# const Classification{N} = PredictMDAPI.Classification{N}
# const BinaryClassification = PredictMDAPI.BinaryClassification

# functions

const fit! = PredictMDAPI.fit!
const get_history = PredictMDAPI.get_history
const get_underlying = PredictMDAPI.get_underlying
const parse_functions! = PredictMDAPI.parse_functions!
const predict = PredictMDAPI.predict
const predict_proba = PredictMDAPI.predict_proba
const set_feature_contrasts! = PredictMDAPI.set_feature_contrasts!
const set_max_epochs! = PredictMDAPI.set_max_epochs!
const transform = PredictMDAPI.transform

const accuracy = PredictMDAPI.accuracy
const auprc = PredictMDAPI.auprc
const aurocc = PredictMDAPI.aurocc
const binary_brier_score = PredictMDAPI.binary_brier_score
const cohen_kappa = PredictMDAPI.cohen_kappa
const f1score = PredictMDAPI.f1score
const false_negative_rate = PredictMDAPI.false_negative_rate
const false_positive_rate = PredictMDAPI.false_positive_rate
const fbetascore = PredictMDAPI.fbetascore
const mean_squared_error = PredictMDAPI.mean_squared_error
const negative_predictive_value = PredictMDAPI.negative_predictive_value
const positive_predictive_value = PredictMDAPI.positive_predictive_value
const prcurve = PredictMDAPI.prcurve
const precision = PredictMDAPI.precision
const r2score = PredictMDAPI.r2score
const recall = PredictMDAPI.recall
const roccurve = PredictMDAPI.roccurve
const root_mean_squared_error = PredictMDAPI.root_mean_squared_error
const sensitivity = PredictMDAPI.sensitivity
const specificity = PredictMDAPI.specificity
const true_negative_rate = PredictMDAPI.true_negative_rate
const true_positive_rate = PredictMDAPI.true_positive_rate

const plotlearningcurve = PredictMDAPI.plotlearningcurve
const plotprcurve = PredictMDAPI.plotprcurve
const plotroccurve = PredictMDAPI.plotroccurve

const load_model = PredictMDAPI.load_model
const save_model = PredictMDAPI.save_model
const save_plot = PredictMDAPI.save_plot


text:
doc"""
    rsi(ta, n=14; wilder=false)

Relative Strength Index

```math
    RSI = \frac{EMA(Up, n)}{EMA(Up, n) + EMA(Dn, n)}
```
"""
function rsi(ta::TimeArray, n::Int=14; wilder::Bool=false)
    ret = diff(ta.values)
    ups = zeros(size(ta.values, 1) - 1, size(ta.values, 2))
    dns = zeros(size(ta.values, 1) - 1, size(ta.values, 2))

    @inbounds for i in 1:size(ta.values, 1) - 1
        for j in 1:size(ta.values, 2)
            if ret[i, j] >= 0
                ups[i, j] += ret[i, j]
            else
                dns[i, j] += ret[i, j]
            end
        end
    end

    upsema = ema(ups, n, wilder=wilder)
    dnsema = abs.(ema(dns, n, wilder=wilder))
    rs     = upsema ./ dnsema

    res  = @. 100 - (100 / (1 + rs))

    cname   = String[]
    cols    = colnames(ta)
    for c in 1:length(cols)
        push!(cname, string(cols[c], "_rsi_", n))
    end

    TimeArray(ta.timestamp[n+1:end], res, cname, ta.meta)
end

doc"""
    cci(ohlc, ma=20, c=0.015)

Commodity Channel Index

```math
    CCI = \frac{P_{typical} - SMA(P_{typical})}{c \times \sigma(P_{typical})}
```

**Reference**

- https://en.wikipedia.org/wiki/Commodity_channel_index

"""
function cci(ohlc::TimeArray, ma::Int=20, c::AbstractFloat=0.015)
    pt = typical(ohlc)
    cci = safediv.((pt .- sma(pt, ma)), moving(mean_abs_dev, pt, ma)) ./ c
    rename(cci, "cci")
end

doc"""
    chaikinoscillator(ohlcv, fast=3, slow=10; h="High", l="Low", c="Close")

**Chaikin Oscillator**

Developed by Marc Chaikin

**Formula**

```math
    Chaikin\ OSC = EMA(ADL, fast) - EMA(ADL, slow)
```

where the [`adl`](@ref) is the Accumulation/Distribution Line.

**Reference**

- [StockCharts]
  (http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:chaikin_oscillator)
"""
function chaikinoscillator(ohlcv::TimeArray, fast::Integer=3, slow::Integer=10;
                           h="High", l="Low", c="Close")
    _adl = adl(ohlcv, h=h, l=l, c=c)
    rename(ema(_adl, fast) .- ema(_adl, slow), ["chaikinoscillator"])
end

doc"""
    macd(ta, fast=12, slow=26, signal=9; wilder=false)

Moving Average Convergence / Divergence

```math
    \begin{align*}
        MACD Bar & = DIF - DEM \\
        DIF & = EMA(P_{close}, fast) - EMA(P_{close}, slow) \\
        DEM & = EMA(DIF, 9) \tag{signal}
    \end{align*}
```

**Return**:

`TimeArray` with 3 columns `["macd", "dif", "signal"]`.

If the input is a multi-column `TimeArray`, the new column names will be
`["A_macd", "B_macd", "A_dif", "B_dif", "A_signal", "B_signal"]`.

"""
function macd(ta::TimeArray{T,N},
              fast::Int=12, slow::Int=26, signal::Int=9;
              wilder::Bool=false) where {T,N}
    dif = ema(ta, fast, wilder=wilder) .- ema(ta, slow, wilder=wilder)
    sig = ema(dif, signal, wilder=wilder)
    osc = dif .- sig

    cols = ["macd", "dif", "signal"]
    new_cols = (N > 1) ? gen_colnames(ta.colnames, cols) : cols

    merge(merge(osc, dif), sig, colnames=new_cols)
end

doc"""
    roc(ta, n)

**Rate of Change**

**Formula**:

```math
    roc = \frac{close_{t} - close_{t-n}}{close_{t-n}}
```

**Reference**:

- [Wikipedia](https://en.wikipedia.org/wiki/Momentum_(technical_analysis))

"""
function roc(ta::TimeArray, n::Integer)
    prev = lag(ta, n)
    rename((ta .- prev) ./ prev, ["$c\_roc_$n" for c ∈ ta.colnames])
end

doc"""
    aroon(ohlc, n=25; h="High", l="Low")

**Aroon Oscillator**

**Formula**

```math
    \begin{align*}
        up   & = \frac{\mathop{argmax}(High_{t-n} \dots High_t)}{n} \times 100 \\
        down & = \frac{\mathop{argmin}(Low_{t-n} \dots Low_t)}{n} \times 100 \\
        osc  & = up - down
    \end{align*}
```

**Reference**

- [StockCharts]
  (http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:aroon_oscillator)
"""
function aroon(ohlc::TimeArray, n::Integer=25; h="High", l="Low")
    up = rename(moving(indmax, ohlc[h], n) ./ n .* 100, "up")
    dn = rename(moving(indmin, ohlc[l], n) ./ n .* 100, "dn")
    osc = rename(up .- dn, "osc")

    merge(merge(up, dn), osc)
end

doc"""
    adx(ohlc, n=14; h="High", l="Low", c="Close")

**Average Directional Movement Index**

Developed by J. Welles Wilder.
This Implementation follows StockCharts.

**Reference**

- [wikipedia]
  (https://en.wikipedia.org/wiki/Average_directional_movement_index)

- [StockCharts]
  (http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:average_directional_index_adx)
"""
function adx(ohlc::TimeArray, n::Integer=14; h="High", l="Low", c="Close")
    dm = relu(merge(ohlc[h] .- lag(ohlc[h]), lag(ohlc[l]) .- ohlc[l]))
    dm.values[findmin(dm.values, 2)[2]] = 0

    dm = wilder_smooth(dm, n)

    tr = wilder_smooth(truerange(ohlc, h=h, l=l, c=c), n)

    di = (dm ./ tr) .* 100
    di = rename(di, ["+di", "-di"])

    dx = @. abs((di["+di"] - di["-di"]) / (di["+di"] + di["-di"])) * 100
    adx = wilder_smooth(dx, n, dx=true)

    rename(merge(adx, merge(dx, di)), ["adx", "dx", "+di", "-di"])
end

doc"""
    stochasticoscillator(ohlc, n=14, fast_d=3, slow_d=3; h="High", l="Low", c="Close")

**Stochastic Oscillator**

A.k.a *%K%D*, or *KD*

**Parameter**

- `n`: period of fast(raw) `%K`

- `fast_d`: MA period of fast `%D`

- `slow_d`: MA period of slow `%D`

**Formula**

```math
    \begin{align*}
        fast\ \%K & = \frac{Close_t - \max(High_{t-n}, \dots, High_t)}
            {\max(High_{t-n}, \dots, High_t) - \min(Low_{t-n}, \dots, Low_t)}
            \times 100 \\
        fast\ \%D & = SMA(fast\ \%K) \\
        slow\ \%D & = SMA(fast\ \%D)
    \end{align*}
```

**Reference**

- [Wikipedia]
  (https://en.wikipedia.org/wiki/Stochastic_oscillator)

- [FMLabs]
  (http://www.fmlabs.com/reference/default.htm?url=StochasticOscillator.htm)
"""
function stochasticoscillator(ohlc::TimeArray, n::Integer=14, fast_d::Integer=3,
                              slow_d::Integer=3; h="High", l="Low", c="Close")
    high = moving(maximum, ohlc[h], n)
    low = moving(minimum, ohlc[l], n)
    fast_k = rename((ohlc[c] .- low) ./ (high .- low) .* 100, "fast_k")
    fast_d = rename(sma(fast_k, fast_d), "fast_d")
    slow_d = rename(sma(fast_d, slow_d), "slow_d")
    merge(merge(fast_k, fast_d), slow_d)
end


text:
module ImageSegmentation

import Base: show

using LinearAlgebra, Statistics
using Images, DataStructures, StaticArrays, ImageFiltering, LightGraphs, SimpleWeightedGraphs, RegionTrees, Distances, StaticArrays, Clustering
import Clustering: kmeans, fuzzy_cmeans

# For efficient hashing of CartesianIndex
if !isdefined(Base.IteratorsMD, :cartindexhash_seed)
    const cartindexhash_seed = UInt == UInt64 ? 0xd60ca92f8284b8b0 : 0xf2ea7c2e
    function Base.hash(ci::CartesianIndex, h::UInt)
        h += cartindexhash_seed
        for i in ci.I
            h = hash(i, h)
        end
        return h
    end
end

include("core.jl")
include("region_growing.jl")
include("felzenszwalb.jl")
include("fast_scanning.jl")
include("watershed.jl")
include("region_merging.jl")
include("meanshift.jl")
include("clustering.jl")

export
    #accessor methods
    labels_map,
    segment_labels,
    segment_pixel_count,
    segment_mean,

    # methods
    seeded_region_growing,
    unseeded_region_growing,
    felzenszwalb,
    fast_scanning,
    watershed,
    hmin_transform,
    region_adjacency_graph,
    rem_segment,
    rem_segment!,
    prune_segments,
    region_tree,
    region_splitting,
    meanshift,
    kmeans,
    fuzzy_cmeans,

    # types
    SegmentedImage,
    ImageEdge

end # module


text:
# define your working directory path here

using Pkg
Pkg.add("DataFrames")
Pkg.add("GLM")
using DataFrames, LinearAlgebra, Statistics, GLM
data = DataFrame(Array{Union{Missing,Float64}}(randn(1000,202)))
data[1:end,202]=1
headers = [ :y ; [ Symbol("x$i") for i = 1:size(data,2) - 1 ] ]
names!(data, headers)

cols = [:x1,:x2,:x3,:x4,:x5,:x6,:x7,:x8,:x9,:x10,:x11,:x12,:x13,:x14,:x15,:x16,:x17,:x18,:x19,:x20,:x21,:x22,:x23,:x24,:x25,:x26,:x27,:x28,:x29,:x30,:x31,:x32,:x33,:x34,:x35,:x36,:x37,:x38,:x39,:x40,:x41,:x42,:x43,:x44,:x45,:x46,:x47,:x48,:x49,:x50,:x51,:x52,:x53,:x54,:x55,:x56,:x57,:x58,:x59,:x60,:x61,:x62,:x63,:x64,:x65,:x66,:x67,:x68,:x69,:x70,:x71,:x72,:x73,:x74,:x75,:x76,:x77,:x78,:x79,:x80,:x81,:x82,:x83,:x84,:x85,:x86,:x87,:x88,:x89,:x90,:x91,:x92,:x93,:x94,:x95,:x96,:x97,:x98,:x99,:x100,:x101,:x102,:x103,:x104,:x105,:x106,:x107,:x108,:x109,:x110,:x111,:x112,:x113,:x114,:x115,:x116,:x117,:x118,:x119,:x120,:x121,:x122,:x123,:x124,:x125,:x126,:x127,:x128,:x129,:x130,:x131,:x132,:x133,:x134,:x135,:x136,:x137,:x138,:x139,:x140,:x141,:x142,:x143,:x144,:x145,:x146,:x147,:x148,:x149,:x150,:x151,:x152,:x153,:x154,:x155,:x156,:x157,:x158,:x159,:x160,:x161,:x162,:x163,:x164,:x165,:x166,:x167,:x168,:x169,:x170,:x171,:x172,:x173,:x174,:x175,:x176,:x177,:x178,:x179,:x180,:x181,:x182,:x183,:x184,:x185,:x186,:x187,:x188,:x189,:x190,:x191,:x192,:x193,:x194,:x195,:x196,:x197,:x198,:x199,:x200,:x201]
depvar_data = convert(Array{Float64},@view(data[1:end, 1]))
expvars_data = convert(Array{Float64}, @view(data[1:end, cols]))
nobs = size(depvar_data, 1)
ncoef = size(expvars_data, 2)

include("qrreg2.jl")

#1
b=@elapsed qrreg(data)

#2 with new dataset (operation were chached)
data2 = DataFrame(Array{Union{Missing,Float64}}(randn(1000,202)))
data2[1:end,202]=1
headers = [ :y ; [ Symbol("x$i") for i = 1:size(data2,2) - 1 ] ]
names!(data2, headers)

b=@elapsed qrreg(data2)

#3
data3 = DataFrame(Array{Union{Missing,Float64}}(randn(1000,202)))
data3[1:end,202]=1
headers = [ :y ; [ Symbol("x$i") for i = 1:size(data3,2) - 1 ] ]
names!(data3, headers)

b=@elapsed qrreg(data3)

#4
data4 = DataFrame(Array{Union{Missing,Float64}}(randn(1000,202)))
data4[1:end,202]=1
headers = [ :y ; [ Symbol("x$i") for i = 1:size(data4,2) - 1 ] ]
names!(data4, headers)

b=@elapsed qrreg(data4)

#5
data5 = DataFrame(Array{Union{Missing,Float64}}(randn(1000,202)))
data5[1:end,202]=1
headers = [ :y ; [ Symbol("x$i") for i = 1:size(data5,2) - 1 ] ]
names!(data5, headers)

b=@elapsed qrreg(data5)

################################################################
#GML regression

#1
data5 = DataFrame(Array{Union{Missing,Float64}}(randn(1000,202)))
data5[1:end,202]=1
headers = [ :y ; [ Symbol("x$i") for i = 1:size(data5,2) - 1 ] ]
names!(data5, headers)

@elapsed ols = lm(@formula(y ~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13+x14+
x15+x16+x17+x18+x19+x20+x21+x22+x23+x24+x25+x26+x27+x28+x29+x30+x31+x32+x33+x34
+x35+x36+x37+x38+x39+x40+x41+x42+x43+x44+x45+x46+x47+x48+x49+x50+x51+x52+x53+
x54+x55+x56+x57+x58+x59+x60+x61+x62+x63+x64+x65+x66+x67+x68+x69+x70+x71+x72+x73
+x74+x75+x76+x77+x78+x79+x80+x81+x82+x83+x84+x85+x86+x87+x88+x89+x90+x91+x92+x93
+x94+x95+x96+x97+x98+x99+x100+x101+x102+x103+x104+x105+x106+x107+x108+x109+x110
+x111+x112+x113+x114+x115+x116+x117+x118+x119+x120+x121+x122+x123+x124+x125+x126
+x127+x128+x129+x130+x131+x132+x133+x134+x135+x136+x137+x138+x139+x140+x141+x142
+x143+x144+x145+x146+x147+x148+x149+x150+x151+x152+x153+x154+x155+x156+x157+x158
+x159+x160+x161+x162+x163+x164+x165+x166+x167+x168+x169+x170+x171+x172+x173+x174
+x175+x176+x177+x178+x179+x180+x181+x182+x183+x184+x185+x186+x187+x188+x189+x190
+x191+x192+x193+x194+x195+x196+x197+x198+x199+x200), data5)

#2 with new dataset (operation were chached)
data5 = DataFrame(Array{Union{Missing,Float64}}(randn(1000,202)))
data5[1:end,202]=1
headers = [ :y ; [ Symbol("x$i") for i = 1:size(data5,2) - 1 ] ]
names!(data5, headers)

@elapsed ols = lm(@formula(y ~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13+x14+
x15+x16+x17+x18+x19+x20+x21+x22+x23+x24+x25+x26+x27+x28+x29+x30+x31+x32+x33+x34
+x35+x36+x37+x38+x39+x40+x41+x42+x43+x44+x45+x46+x47+x48+x49+x50+x51+x52+x53+
x54+x55+x56+x57+x58+x59+x60+x61+x62+x63+x64+x65+x66+x67+x68+x69+x70+x71+x72+x73
+x74+x75+x76+x77+x78+x79+x80+x81+x82+x83+x84+x85+x86+x87+x88+x89+x90+x91+x92+x93
+x94+x95+x96+x97+x98+x99+x100+x101+x102+x103+x104+x105+x106+x107+x108+x109+x110
+x111+x112+x113+x114+x115+x116+x117+x118+x119+x120+x121+x122+x123+x124+x125+x126
+x127+x128+x129+x130+x131+x132+x133+x134+x135+x136+x137+x138+x139+x140+x141+x142
+x143+x144+x145+x146+x147+x148+x149+x150+x151+x152+x153+x154+x155+x156+x157+x158
+x159+x160+x161+x162+x163+x164+x165+x166+x167+x168+x169+x170+x171+x172+x173+x174
+x175+x176+x177+x178+x179+x180+x181+x182+x183+x184+x185+x186+x187+x188+x189+x190
+x191+x192+x193+x194+x195+x196+x197+x198+x199+x200), data5)

#3
data5 = DataFrame(Array{Union{Missing,Float64}}(randn(1000,202)))
data5[1:end,202]=1
headers = [ :y ; [ Symbol("x$i") for i = 1:size(data5,2) - 1 ] ]
names!(data5, headers)

@elapsed ols = lm(@formula(y ~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13+x14+
x15+x16+x17+x18+x19+x20+x21+x22+x23+x24+x25+x26+x27+x28+x29+x30+x31+x32+x33+x34
+x35+x36+x37+x38+x39+x40+x41+x42+x43+x44+x45+x46+x47+x48+x49+x50+x51+x52+x53+
x54+x55+x56+x57+x58+x59+x60+x61+x62+x63+x64+x65+x66+x67+x68+x69+x70+x71+x72+x73
+x74+x75+x76+x77+x78+x79+x80+x81+x82+x83+x84+x85+x86+x87+x88+x89+x90+x91+x92+x93
+x94+x95+x96+x97+x98+x99+x100+x101+x102+x103+x104+x105+x106+x107+x108+x109+x110
+x111+x112+x113+x114+x115+x116+x117+x118+x119+x120+x121+x122+x123+x124+x125+x126
+x127+x128+x129+x130+x131+x132+x133+x134+x135+x136+x137+x138+x139+x140+x141+x142
+x143+x144+x145+x146+x147+x148+x149+x150+x151+x152+x153+x154+x155+x156+x157+x158
+x159+x160+x161+x162+x163+x164+x165+x166+x167+x168+x169+x170+x171+x172+x173+x174
+x175+x176+x177+x178+x179+x180+x181+x182+x183+x184+x185+x186+x187+x188+x189+x190
+x191+x192+x193+x194+x195+x196+x197+x198+x199+x200), data5)

#4
data5 = DataFrame(Array{Union{Missing,Float64}}(randn(1000,202)))
data5[1:end,202]=1
headers = [ :y ; [ Symbol("x$i") for i = 1:size(data5,2) - 1 ] ]
names!(data5, headers)

@elapsed ols = lm(@formula(y ~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13+x14+
x15+x16+x17+x18+x19+x20+x21+x22+x23+x24+x25+x26+x27+x28+x29+x30+x31+x32+x33+x34
+x35+x36+x37+x38+x39+x40+x41+x42+x43+x44+x45+x46+x47+x48+x49+x50+x51+x52+x53+
x54+x55+x56+x57+x58+x59+x60+x61+x62+x63+x64+x65+x66+x67+x68+x69+x70+x71+x72+x73
+x74+x75+x76+x77+x78+x79+x80+x81+x82+x83+x84+x85+x86+x87+x88+x89+x90+x91+x92+x93
+x94+x95+x96+x97+x98+x99+x100+x101+x102+x103+x104+x105+x106+x107+x108+x109+x110
+x111+x112+x113+x114+x115+x116+x117+x118+x119+x120+x121+x122+x123+x124+x125+x126
+x127+x128+x129+x130+x131+x132+x133+x134+x135+x136+x137+x138+x139+x140+x141+x142
+x143+x144+x145+x146+x147+x148+x149+x150+x151+x152+x153+x154+x155+x156+x157+x158
+x159+x160+x161+x162+x163+x164+x165+x166+x167+x168+x169+x170+x171+x172+x173+x174
+x175+x176+x177+x178+x179+x180+x181+x182+x183+x184+x185+x186+x187+x188+x189+x190
+x191+x192+x193+x194+x195+x196+x197+x198+x199+x200), data5)


text:
using SnoopCompile

### Log the compiles
# This only needs to be run once (to generate "/tmp/colortypes_compiles.csv")

SnoopCompile.@snoop "/tmp/colortypes_compiles.csv" begin
    include(Pkg.dir("ColorTypes", "test","runtests.jl"))
end

### Parse the compiles and generate precompilation scripts
# This can be run repeatedly

# IMPORTANT: we must have the module(s) defined for the parcelation
# step, otherwise we will get no precompiles for the Colortypes module
using ColorTypes

data = SnoopCompile.read("/tmp/colortypes_compiles.csv")

# Use these two lines if you want to create precompile functions for
# individual packages
pc, discards = SnoopCompile.parcel(data[end:-1:1,2])
SnoopCompile.write("/tmp/precompile", pc)

# Use these two lines if you want to add to your userimg.jl
pc = SnoopCompile.format_userimg(data[end:-1:1,2])
SnoopCompile.write("/tmp/userimg_ColorTypes.jl", pc)

function notisempty(filename, minlength=1)
    @test isfile(filename)
    @test length(readlines(filename)) >= minlength
end

notisempty("/tmp/precompile/precompile_ColorTypes.jl", 100)
notisempty("/tmp/precompile/precompile_FixedPointNumbers.jl", 2)
notisempty("/tmp/userimg_ColorTypes.jl", 100)