text:
module ClusteringGA
# Clustering Genetic Algorithm

# Hruschka, Eduardo & Ebecken, Nelson. (2003). A genetic algorithm for cluster analysis.
# Intell. Data Anal.. 7. 15-25. 10.3233/IDA-2003-7103.

using Random
using Statistics
using LinearAlgebra
using Clustering

# Clustering Genetic Algorithm
export cga, CGAData, CGAResult

# Results Interface

"""
    CGAResult <: ClusteringResult

Contains the results from the computation of [`cga`](@ref).

#### Members
  * `assignments::Vector{Int}` element-to-cluster assignments (`n`)
  * `counts::Vector{Int}`      number of samples assigned to each cluster (`k`)
  * `found_gen::Int`           first generation where the elite was found
  * `total_gen::Int`           total generations the GA has been run
"""
struct CGAResult <: ClusteringResult
    assignments::Vector{Int}    # element-to-cluster assignments (n)
    counts::Vector{Int}         # number of samples assigned to each cluster (k)
    found_gen::Int              # first generation where the elite was found
    total_gen::Int              # total generations the GA has been run
end

"""
    CGAData{S, T<:Real}

Contains the data used in the computation of [`cga`](@ref).

The user does not need to query this object, but can use this as an opaque object. The results of the computation can be obtained from [`CGAResult`](@ref) object.

"""
mutable struct CGAData{S, T <: Real,
                       V <: AbstractVector{S},
                       M <: AbstractMatrix{T}}
    objs::V
    dist::M
    N::Int
    generations::Int
    curr_gen::Int
    population::Matrix{Int}
    scratch::Matrix{Int}
    vpopulation::Vector{SubArray}
    vscratch::Vector{SubArray}

    p::Matrix{Int}
    np::Matrix{Int}
    vp::Vector{SubArray}
    vnp::Vector{SubArray}

    elitev::Float64
    elited::Vector{Int}
    elite_gen::Int
    val::Vector{Float64}
    function CGAData{S, T, V, M}(objs::V,
                                 dist::M,
                                 N::Int,
                                 generations::Int) where {S,
                                                          V <: AbstractVector{S},
                                                          T <: Real,
                                                          M <: AbstractMatrix{T}}
        sz = size(dist)
        @assert sz[1] == sz[2] "The dist matrix should be symmetric"
        population = p  = Matrix{Int}(undef, (sz[1]+1, N))
        scratch =    np = Matrix{Int}(undef, (sz[1]+1, N))
        vpopulation = vp  = [@view  p[:, i] for i = 1:N ]
        vscratch    = vnp = [@view np[:, i] for i = 1:N ]
        curr_gen = elite_gen = 0
        elited = fill(1, sz[1]+1)
        elitev = 0
        val = Vector{Float64}(undef, N)
        return new(objs, dist, N, generations, curr_gen, population,
                   scratch, vpopulation, vscratch, p, np, vp, vnp,
                   elitev, elited, elite_gen, val)
    end
end

# Ensuring silhouttes function does not throw when one cluster is found
# also assigns the silhouette values to zeros to not be picked up for
# selection significantly.
function silhouettes_nothrow(assignments::AbstractVector{<:Integer},
                             dists::AbstractMatrix{T}) where {T <: Real}
    counts = fill(0, maximum(assignments))
    for a in assignments
        counts[a] += 1
    end
    k = length(counts)
    k > 1 && return silhouettes(assignments, counts, dists)
    S = typeof((one(T) + one(T))/2)
    return fill(zero(S), length(assignments))
end

function fitness(idx::Int, data::CGAData)
    nobj = size(data.population, 1) - 1
    vw = @view data.population[1:nobj, idx]
    return fitness(vw, data)
end

function fitness(assignment::AbstractVector{Int}, data::CGAData)
    nobj = size(data.population, 1) - 1
    vw = @view assignment[1:nobj]
    return 1.0 + mean(silhouettes_nothrow(vw, data.dist))
end

"""
    cga(objects::V,
        distances::M=distance_matrix(objs),
        N::Int=length(objs)*20,
        generations::Int=50) where {S, V <: AbstractVector{S},
                                    T <: Real, M <: AbstractMatrix{T}}

Compute the optimal clustering in the data by Genetic Algorithm over the computed mean silhouettes.

  * `objects` the vector of the objects for which clusters are to be computed.
  * `distances` the distance matrix providing the pairwise distances bewtween the `objects`
  * `N` the population size for GA computation
  * `generations` number of generations the GA has to be run

The `fitness` function used is `1+mean(silhouettes())` to ensure positive values.

#### Return Values
The method returns a tuple of ([`CGAData`](@ref), [`CGAResult`](@ref))

#### References
  1. Hruschka, Eduardo & Ebecken, Nelson. (2003). A genetic algorithm for cluster analysis. Intell. Data Anal.. 7. 15-25. 10.3233/IDA-2003-7103.
"""
function cga(objs::V,
             dist::M=distance_matrix(objs),
             N::Int=length(objs)*20,
             generations::Int=50) where {S, V <: AbstractVector{S},
                                         T <: Real, M <: AbstractMatrix{T}}
    data = CGAData{S, T, V, M}(objs, dist, N, generations)
    init_population!(data.population)
    selection!(data)
    nobj = size(data.population, 1) - 1

    N_4  = div(N, 4)
    N_2  = div(N, 2)
    N3_4 = div(3N, 4)
    for gen=1:generations
        data.curr_gen = gen
        for i = 1:N_4
            crossover!(i, N_2-i+1, data)
            mutation_split!(N_2+i, data)
            mutation_merge!(N3_4+i, data)
        end
        selection!(data)
        fval = fitness(1, data)
    end
    elited = data.elited
    assignments = elited[1:nobj]
    k = elited[end]
    counts = Vector{Int}(undef, k)
    for i=1:k
        counts[i] = length(collect(1:nobj)[assignments .== i])
    end
    result = CGAResult(assignments, counts, data.elite_gen, generations)
    return data, result
end

function selection!(data::CGAData)
    nobj = size(data.population, 1) - 1
    population = data.population
    N = data.N
    dist = data.dist
    val = data.val
    elitev, elitei = -Inf, -1
    for i = 1:N
        tval = fitness(i, data)
        elitev < tval && ((elitev, elitei) = (tval, i))
        val[i] = tval
    end
    if elitev > data.elitev
        data.elitev = elitev
        copyto!(data.elited, data.vpopulation[elitei])
        data.elite_gen = data.curr_gen
    end
    for i=2:N
        val[i] += val[i-1]
    end
    p, np, vp, vnp = data.population === data.p ?
        (data.p, data.np, data.vp, data.vnp) :
        (data.np, data.p, data.vnp, data.vp)
    lastv = val[end]
    copyto!(vnp[1], data.elited)
    rn = rand(N)*lastv
    for i = 2:N
        idx = searchsortedfirst(val, rn[i])
        copyto!(vnp[i], vp[idx])
    end
    data.population, data.scratch, data.vpopulation, data.vscratch =
        np, p, vnp, vp
    return
end

@inline function assign_nearest_cluster!(c, objs, nobj, g)
    lg = length(g)
    cgs = Vector{Vector{Float64}}(undef, lg)
    tobjs = similar(objs, nobj)
    for j = 1:lg
        ntobjs = 0
        for i = 1:nobj
            if c[i] == g[j]
                tobjs[ntobjs+=1] = objs[i]
            end
        end
        cgs[j] = centroid(tobjs[1:ntobjs])
    end

    #da = similar(objs[1], Float64)
    for i=1:nobj
        if c[i] == 0
            mind, mini = Inf, 0
            for j = 1:lg
                # copyto!(da, cgs[j])
                da = cgs[j] - objs[i]
                dda = dot(da, da)
                dda < mind && ((mind, mini) = (dda, j))
            end
            c[i] = g[mini]
        end
    end
    return
end

@inline function insert_not_found!(gx, gxc, val)
    found = false
    for j = 1:gxc
        if val == gx[j]
            found = true
            break
        end
    end
    !found && (gx[gxc+=1] = val)
    return gxc
end

@inline function get_xoverchild(a, c, k, l, g)
    gx  = Vector{Int}(undef, k)
    gxc = 0
    for i=1:l
        ai = a[i]
        c[i] = ai in g ? ai : 0
        c[i] == 0 && continue
        gxc = insert_not_found!(gx, gxc, c[i])
    end
    return resize!(gx, gxc)
end

@inline function pre_crossover!(a::AbstractVector{Int}, b::AbstractVector{Int},
                                c::AbstractVector{Int}, d::AbstractVector{Int},
                                g::AbstractVector{Int})
    k1, k2 = a[end], b[end]
    k = max(k1, k2)
    l = lastindex(a) - 1

    g1 = get_xoverchild(a, c, k, l, g)
    g2 = get_xoverchild(b, d, k, l, g1)

    gc,  gd  = Vector{Int}(undef, k), Vector{Int}(undef, k)
    gcn, gdn = 0, 0
    for i=1:l
        if c[i] == 0
            if !(b[i] in g1)
                c[i] = b[i]
                gcn = insert_not_found!(gc, gcn, c[i])
            end
        else
            gcn = insert_not_found!(gc, gcn, c[i])
        end
        if d[i] == 0
            if !(a[i] in g2)
                d[i] = a[i]
                gdn = insert_not_found!(gd, gdn, d[i])
            end
        else
            gdn = insert_not_found!(gd, gdn, d[i])
        end
    end
    return resize!(gc, gcn), resize!(gd, gdn)
end

@inline function crossover!(ia::Int, ib::Int, data::CGAData)
    vpopulation = data.vpopulation
    vscratch    = data.vscratch
    objs       = data.objs
    dist       = data.dist
    nobj = size(data.population, 1) - 1
    a = vpopulation[ia]
    b = vpopulation[ib]
    a == b && return

    c = vscratch[ia]
    d = vscratch[ib]
    k = a[end]
    g = Int[]
    while !(0 < length(g) <= k)
        randsubseq!(g, collect(1:k), 0.5)
    end
    gc, gd = pre_crossover!(a, b, c, d, g)
    assign_nearest_cluster!(c, objs, nobj, gc)
    assign_nearest_cluster!(d, objs, nobj, gd)
    copyto!(a, 1, c, 1, nobj)
    copyto!(b, 1, d, 1, nobj)
    normalize_assignment!(a, c, nobj)
    normalize_assignment!(b, c, nobj)
    return
end

function mutation_split!(ia::Int, data::CGAData{S, T}) where {S, T <: Real}
    vpopulation = data.vpopulation
    vscratch    = data.vscratch
    objs       = data.objs
    nobj = size(data.population, 1) - 1
    a  = vpopulation[ia]
    nc = a[end]
    sc = rand(1:nc)
    ids = Vector{Int}(undef, nobj)
    ln = 0
    for i=1:nobj
        if a[i] == sc
            ids[ln+=1] = i
        end
    end
    resize!(ids, ln)
    cg = centroid(@view objs[ids])
    d = Dict{Int, Float64}()

    l2 = div(ln+1, 2)
    sort!(ids, lt = (x, y) -> begin
          if !haskey(d, x)
          dx = cg - objs[x]
          ddx = dot(dx, dx)
          d[x] = ddx
          end
          if !haskey(d, y)
          dy = cg - objs[y]
          ddy = dot(dy, dy)
          d[y] = ddy
          end
          return d[x] < d[y]
          end
          )
    for ii = l2:ln
        a[ids[ii]] = (nc + 1)
    end
    sa = vscratch[ia]
    normalize_assignment!(a, sa, nobj)
    return
end

function mutation_merge!(ia::Int, data::CGAData{S, T}) where {S, T<:Real}
    vpopulation = data.vpopulation
    vscratch    = data.vscratch
    objs       = data.objs
    dist       = data.dist
    nobj = size(data.population, 1) - 1
    a  = vpopulation[ia]
    nc = a[end]
    if nc > 2
        sc = rand(1:nc)
        cgs = Vector{Vector{Float64}}(undef, nc)
        tobjs = similar(objs, nobj)
        for j = 1:nc
            ntobjs = 0
            for i = 1:nobj
                if a[i] == j
                    tobjs[ntobjs+=1] = objs[i]
                end
            end
            cgs[j] = centroid(@view tobjs[1:ntobjs])
        end

        cg = cgs[sc]
        d = Inf
        j = sc
        for i = 1:nc
            i == sc && continue
            dcg = cgs[i] - cg
            ddcg = dot(dcg, dcg)
            d > ddcg && ((d, j) = (ddcg, i))
        end
        for i = 1:nobj
            a[i] == sc && (a[i] = j)
        end
        sa = vscratch[ia]
        normalize_assignment!(a, sa, nobj)
    end
    return
end

normalize_assignment!(v) =
    normalize_assignment!(v, Vector{Int}(undef, length(v)-1), length(v)-1)

@inline function normalize_assignment!(v, d, nobj)
    fill!(d, 0)
    curr = 0
    for i = 1:nobj
        @inbounds vi = v[i]
        @inbounds d[vi] == 0 && (d[vi] = (curr += 1))
        @inbounds v[i] = d[vi]
    end
    @inbounds v[end] = curr
    return
end

function init_population!(p::Matrix)
    nobj, N = size(p)
    nobj -= 1
    curr = 2
    d = zeros(Int, nobj)
    for i = 1:N
        vw = @view p[:, i]
        rand!(vw, 1:curr)
        normalize_assignment!(vw, d, nobj)
        curr += 1
        curr > nobj && (curr = 2)
    end
    return
end

far_object(objs::AbstractVector{S}) where {T <: Real, S <: AbstractVector{T}} =
    fill(Inf, size(objs[1]))

function centroid(objs::AbstractVector{S}) where {T <: Real, S <: AbstractVector{T}}
    l = length(objs)
    result = fill(0.0, length(objs[1]))
    for i=1:lastindex(objs)
        result .+= objs[i]
    end
    return result /= l
end

function distance_matrix(objs::AbstractVector{S}) where {T <: Real,
                                                         S <: AbstractVector{T}}
    l = length(objs)
    m = Matrix{Float64}(undef, (l, l))
    for i = 1:l
        for j = i:l
            v = objs[i] - objs[j]
            @inbounds m[i, j] = m[j, i] = sqrt(dot(v, v))
        end
    end
    return m
end
end # Module ClusteringGA


text:
function test_sampler_cube(sampler, n_samples; opts=Dict(), add_slice=nothing)
	""" Creates a 3-dimensional unit cube and uses the given sampler
		to sample from it the specified number of times. Then plots the
		sampled distribution
	"""
	A = Array{Float64}(undef, 0, 3)
	b = Array{Float64}(undef, 0, 1)
	p0 = [1.0, 0.5, 0.5]
	v0 = [1.0, 0.0, 0.0]
	p1 = [0, 0.5, 0.5]
	v1 = [-1.0, 0.0, 0.0]
	for (p, v) in zip((p0,p1), (v0,v1))
	for i in 0:2
		p_ = circshift(p, i)
		v_ = circshift(v, i)
		A = vcat(A, v_')
		b = vcat(b, v_' * p_)
	end
	end

	if !isnothing(add_slice)
		v_ = add_slice[1]
		p_ = add_slice[2]
		A = vcat(A, v_')
		b = vcat(b, v_' * p_)
	end

	samples = Array{Float64}(undef, 0, 3)

	x0, r = OPFLearn.chebyshev_center(A, b)

	k = 0
	while k < n_samples
		sample = sampler(A, b, x0'; opts...)'
		samples = vcat(samples, sample)
		k = k + size(sample, 1)
	end
	return samples
end

"""
using PyPlot
const plt = PyPlot
function plot_samples_series(samples, loads=(1, 2, 3); sampler="", case="")

	num_samples = size(samples,1)

    fig = plt.figure()
    ax = fig.add_subplot(projection="3d")

    # Plots with colors indicating sample number
    p = ax.scatter(samples[:, loads[1]],
                   samples[:, loads[2]],
                   samples[:, loads[3]],
                   c=collect(1: size(samples,1)), cmap="plasma")
    fig.colorbar(p, label="Sample Number")
    ax.set_xlabel("Load " * string(loads[1]))
    ax.set_ylabel("Load " * string(loads[2]))
    ax.set_zlabel("Load " * string(loads[3]))
    plt.title(case * " load powers, " * string(num_samples) * " samples, " * sampler * " sampler")
    plt.tight_layout()
    plt.show()
end

function run_test_sampler(num_samples, add_slice=[[1.0,1.0,1.0],[0.5,0.5,0.5]])
	opts = Dict((:method=>"achr"))
	samples = test_sampler_cube(sample_polytope_cprnd, num_samples, opts=opts, add_slice=add_slice)
	plot_samples_series(samples, sampler="cprnd", case="Unit Cube")
end
"""

@testset "sampler uniformity check" begin
    @testset "cprnd unit cube" begin
		num_samples = 10000  # 10,000 samples should take approximately 30 seconds
		opts = Dict((:method=>"achr"))
		samples = test_sampler_cube(sample_polytope_cprnd, num_samples, opts=opts, add_slice=nothing)
		@test size(samples, 1) == num_samples
		@test size(samples, 2) == 3
		@test isapprox(maximum(samples[:,1]), 1; atol = 1e-1)
		@test isapprox(maximum(samples[:,2]), 1; atol = 1e-1)
		@test isapprox(maximum(samples[:,3]), 1; atol = 1e-1)
		test_points = ((1.0, 1.0, 1.0),
					   (0.5, 0.5, 0.5),
					   (0.0, 0.0, 0.0),
					   [circshift([1.0, 0.0, 0.0],i) for i in 1:3]...,
					   [circshift([1.0, 1.0, 0.0],i) for i in 1:3]...,
					   [circshift([1.0, 0.5, 0.5],i) for i in 1:3]...,
					   [circshift([0.0, 0.5, 0.5],i) for i in 1:3]...,
					   )

		for p in test_points  # Test that at least one point sampled is near a test point
		@test isapprox(minimum([sum((sample .- p).^2) for sample in eachrow(samples)]), 0; atol = 5e-2)
		end
    end

    @testset "cprnd unit cube with slice" begin
		num_samples = 10000  # 10,000 samples should take approximately 30 seconds
		opts = Dict((:method=>"achr"))
		normal = [1.0, 1.0, 1.0]
		point = [0.5, 0.5, 0.0]
		slice = [normal, point]
		samples = test_sampler_cube(sample_polytope_cprnd, num_samples, opts=opts, add_slice=slice)
		@test size(samples, 1) == num_samples
		@test size(samples, 2) == 3
		@test isapprox(maximum(samples[:,1]), 1; atol = 1e-1)
		@test isapprox(maximum(samples[:,2]), 1; atol = 1e-1)
		@test isapprox(maximum(samples[:,3]), 1; atol = 1e-1)
		inside_test_points = ((0.0, 0.0, 0.0),
							  (0.5, 0.5, 0.0),
							  (0.25, 0.25, 0.25),
							  [circshift([1.0, 0.0, 0.0],i) for i in 1:3]...,
						      )

		outside_test_points = ((1.0, 1.0, 1.0),
							   (0.5, 0.5, 0.5),
							   [circshift([1.0, 1.0, 0.0],i) for i in 1:3]...,
							   )
		for p in inside_test_points  # Test that at least one point sampled is near a test point
		@test isapprox(minimum([sum((sample .- p).^2) for sample in eachrow(samples)]), 0; atol = 5e-2)
		end
		for p in outside_test_points  # Test that there are no points sampled is near a test point
		@test minimum([sum((sample .- p).^2) for sample in eachrow(samples)]) > 1e-2
		end
    end
end

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

@userplot DistPlot1D

@recipe function f(dp::DistPlot1D; quantiles=[0.25,0.50,0.75])
  # retrieve inputs
  sdata = dp.args[1]
  var   = dp.args[2]

  # fit spatial statistics
  h = EmpiricalHistogram(sdata, var)
  q = quantile(sdata, var, quantiles)
  μ = mean(sdata, var)

  legend --> false
  framestyle --> :box
  grid --> false
  xguide --> var

  # plot histogram
  @series begin
    seriestype --> :step
    seriescolor --> :black
    h
  end

  # plot quantiles
  @series begin
    seriestype --> :vline
    primary --> false
    seriescolor --> :black
    linestyle --> :dash
    q
  end

  # plot mean
  @series begin
    seriestype --> :vline
    primary --> false
    seriescolor --> :green
    [μ]
  end
end


text:
################################################################################
## New storage to reference model
################################################################################
#used
function add_candidate_storage!(ref::Dict{Symbol,<:Any}, data::Dict{String,<:Any})
    for (n, nw_ref) in ref[:nw]
        if haskey(nw_ref, :ne_storage)
            bus_storage_ne = Dict([(i, []) for (i,bus) in nw_ref[:bus]])
            for (i,storage) in nw_ref[:ne_storage]
                push!(bus_storage_ne[storage["storage_bus"]], i)
            end
            nw_ref[:bus_storage_ne] = bus_storage_ne
        end
    end
end

################################################################################
## Model references for distribution
################################################################################

"Like ref_add_ne_branch!, but ne_buspairs are built using calc_buspair_parameters_allbranches"
function ref_add_ne_branch_allbranches!(ref::Dict{Symbol,<:Any}, data::Dict{String,<:Any})
    for (nw, nw_ref) in ref[:nw]
        if !haskey(nw_ref, :ne_branch)
            Memento.error(_LOGGER, "required ne_branch data not found")
        end

        nw_ref[:ne_branch] = Dict(x for x in nw_ref[:ne_branch] if (x.second["br_status"] == 1 && x.second["f_bus"] in keys(nw_ref[:bus]) && x.second["t_bus"] in keys(nw_ref[:bus])))

        nw_ref[:ne_arcs_from] = [(i,branch["f_bus"],branch["t_bus"]) for (i,branch) in nw_ref[:ne_branch]]
        nw_ref[:ne_arcs_to]   = [(i,branch["t_bus"],branch["f_bus"]) for (i,branch) in nw_ref[:ne_branch]]
        nw_ref[:ne_arcs] = [nw_ref[:ne_arcs_from]; nw_ref[:ne_arcs_to]]

        ne_bus_arcs = Dict((i, []) for (i,bus) in nw_ref[:bus])
        for (l,i,j) in nw_ref[:ne_arcs]
            push!(ne_bus_arcs[i], (l,i,j))
        end
        nw_ref[:ne_bus_arcs] = ne_bus_arcs

        if !haskey(nw_ref, :ne_buspairs)
            ismc = haskey(nw_ref, :conductors)
            cid = nw_ref[:conductor_ids]
            nw_ref[:ne_buspairs] = calc_buspair_parameters_allbranches(nw_ref[:bus], nw_ref[:ne_branch], cid, ismc)
        end
    end
end

"""
Add to `ref` the following keys:
- `:frb_branch`: the set of `branch`es whose `f_bus` is the reference bus;
- `:frb_ne_branch`: the set of `ne_branch`es whose `f_bus` is the reference bus.
"""
function ref_add_frb_branch!(ref::Dict{Symbol,Any}, data::Dict{String,<:Any})
    for (nw, nw_ref) in ref[:nw]
        ref_bus_id = first(keys(nw_ref[:ref_buses]))

        frb_branch = Dict{Int,Any}()
        for (i,br) in nw_ref[:branch]
            if br["f_bus"] == ref_bus_id
                frb_branch[i] = br
            end
        end
        nw_ref[:frb_branch] = frb_branch

        if haskey(nw_ref, :ne_branch)
            frb_ne_branch = Dict{Int,Any}()
            for (i,br) in nw_ref[:ne_branch]
                if br["f_bus"] == ref_bus_id
                    frb_ne_branch[i] = br
                end
            end
            nw_ref[:frb_ne_branch] = frb_ne_branch
        end
    end
end

"""
Add to `ref` the following keys:
- `:oltc_branch`: the set of `frb_branch`es that are OLTCs;
- `:oltc_ne_branch`: the set of `frb_ne_branch`es that are OLTCs.
"""
function ref_add_oltc_branch!(ref::Dict{Symbol,Any}, data::Dict{String,<:Any})
    for (nw, nw_ref) in ref[:nw]
        if !haskey(nw_ref, :frb_branch)
            Memento.error(_LOGGER, "ref_add_oltc_branch! must be called after ref_add_frb_branch!")
        end
        oltc_branch = Dict{Int,Any}()
        for (i,br) in nw_ref[:frb_branch]
            if br["transformer"] && haskey(br, "tm_min") && haskey(br, "tm_max") && br["tm_min"] < br["tm_max"]
                oltc_branch[i] = br
            end
        end
        nw_ref[:oltc_branch] = oltc_branch

        if haskey(nw_ref, :frb_ne_branch)
            oltc_ne_branch = Dict{Int,Any}()
            for (i,br) in nw_ref[:ne_branch]
                if br["transformer"] && haskey(br, "tm_min") && haskey(br, "tm_max") && br["tm_min"] < br["tm_max"]
                    oltc_ne_branch[i] = br
                end
            end
            nw_ref[:oltc_ne_branch] = oltc_ne_branch
        end
    end
end

"Like calc_buspair_parameters, but retains indices of all the branches and drops keys that depend on branch"
function calc_buspair_parameters_allbranches(buses, branches, conductor_ids, ismulticondcutor)
    bus_lookup = Dict(bus["index"] => bus for (i,bus) in buses if bus["bus_type"] != 4)

    branch_lookup = Dict(branch["index"] => branch for (i,branch) in branches if branch["br_status"] == 1 && haskey(bus_lookup, branch["f_bus"]) && haskey(bus_lookup, branch["t_bus"]))

    buspair_indexes = Set((branch["f_bus"], branch["t_bus"]) for (i,branch) in branch_lookup)

    bp_branch = Dict((bp, Int[]) for bp in buspair_indexes)

    if ismulticondcutor
        bp_angmin = Dict((bp, [-Inf for c in conductor_ids]) for bp in buspair_indexes)
        bp_angmax = Dict((bp, [ Inf for c in conductor_ids]) for bp in buspair_indexes)
    else
        @assert(length(conductor_ids) == 1)
        bp_angmin = Dict((bp, -Inf) for bp in buspair_indexes)
        bp_angmax = Dict((bp,  Inf) for bp in buspair_indexes)
    end

    for (l,branch) in branch_lookup
        i = branch["f_bus"]
        j = branch["t_bus"]

        if ismulticondcutor
            for c in conductor_ids
                bp_angmin[(i,j)][c] = max(bp_angmin[(i,j)][c], branch["angmin"][c])
                bp_angmax[(i,j)][c] = min(bp_angmax[(i,j)][c], branch["angmax"][c])
            end
        else
            bp_angmin[(i,j)] = max(bp_angmin[(i,j)], branch["angmin"])
            bp_angmax[(i,j)] = min(bp_angmax[(i,j)], branch["angmax"])
        end

        bp_branch[(i,j)] = push!(bp_branch[(i,j)], l)
    end

    buspairs = Dict((i,j) => Dict(
        "branches"=>bp_branch[(i,j)],
        "angmin"=>bp_angmin[(i,j)],
        "angmax"=>bp_angmax[(i,j)],
        "vm_fr_min"=>bus_lookup[i]["vmin"],
        "vm_fr_max"=>bus_lookup[i]["vmax"],
        "vm_to_min"=>bus_lookup[j]["vmin"],
        "vm_to_max"=>bus_lookup[j]["vmax"]
        ) for (i,j) in buspair_indexes
    )

    return buspairs
end


text:
module Images

import Base.Order: Ordering, ForwardOrdering, ReverseOrdering
import Base.Graphics: width, height
importall Base

using Color

# We need a couple of extra features not present in Color
immutable RGB8 <: ColorValue
    r::Uint8
    g::Uint8
    b::Uint8

    function RGB8(r::Real, g::Real, b::Real)
        new(r, g, b)
    end

    RGB8() = RGB8(0,0,0)
end

typemin(::Type{RGB}) = RGB(0,0,0)
typemax(::Type{RGB}) = RGB(1,1,1)

if VERSION.minor < 3
    using Cartesian
else
    using Base.Cartesian
    include("compatibility/forcartesian.jl")
end

# if isdefined(module_parent(Images), :Grid)
#     import ..Grid.restrict
# end

include("core.jl")
include("iterator.jl")
const have_imagemagick = include("ioformats/libmagickwand.jl")
@osx_only include("ioformats/OSXnative.jl")
include("io.jl")
include("scaling.jl")
include("labeledarrays.jl")
include("algorithms.jl")
include("connected.jl")
include("edge.jl")

__init__() = LibMagick.init()

function precompile()
    for T in (Uint8, Uint16, Int, Float32, Float64)
        Tdiv = typeof(one(T)/2)
        for N = 2:3
            precompile(restrict!, (Array{Tdiv, N}, Array{T,N}, Int))
            precompile(imfilter, (Array{T,N}, Array{Float64,N}))
            precompile(imfilter, (Array{T,N}, Array{Float64,2}))
            precompile(imfilter, (Array{T,N}, Array{Float32,N}))
            precompile(imfilter, (Array{T,N}, Array{Float32,2}))
        end
    end
    for T in (Float32, Float64)
        for N = 2:3
            precompile(_imfilter_gaussian!, (Array{T,N}, Vector{Float64}))
            precompile(_imfilter_gaussian!, (Array{T,N}, Vector{Int}))
            precompile(imfilter_gaussian_no_nans!, (Array{T,N}, Vector{Float64}))
            precompile(imfilter_gaussian_no_nans!, (Array{T,N}, Vector{Int}))
            precompile(fft, (Array{T,N},))
        end
    end
    for T in (Complex{Float32}, Complex{Float64})
        for N = 2:3
            precompile(ifft, (Array{T,N},))
        end
    end
end

export # types
    AbstractImage,
    AbstractImageDirect,
    AbstractImageIndexed,
    Image,
    ImageCmap,
    BitShift,
    ClipMin,
    ClipMax,
    ClipMinMax,
    LabeledArray,
    Overlay,
    OverlayImage,
    ScaleAutoMinMax,
    ScaleInfo,
    ScaleMinMax,
    ScaleNone,
    ScaleSigned,
    SliceData,

    # constants
    palette_fire,
    palette_gray32,
    palette_gray64,
    palette_rainbow,

    # core functions
    assert2d,
    assert_scalar_color,
    assert_timedim_last,
    assert_xfirst,
    assert_yfirst,
    colordim,
    colorspace,
    coords,
    coords_spatial,
    data,
    dimindex,
    dimindexes,
    getindexim,
    grayim,
    colorim,
    height,
    isdirect,
    isxfirst,
    isyfirst,
    limits,
    maxabsfinite,
    maxfinite,
    minfinite,
    ncolorelem,
    nimages,
    pixelspacing,
    properties,
    refim,
    rerange!,
    reslice!,
    restrict,
    sdims,
    size_spatial,
    share,
    sliceim,
    spacedirections,
    spatialorder,
    spatialpermutation,
    spatialproperties,
    storageorder,
    subim,
    timedim,
    width,
    widthheight,

    # io functions
    add_image_file_format,
    imread,
    imwrite,
    loadformat,

    # iterator functions
    first_index,
    iterate_spatial,
    parent,

    # color-related functions
    alphaval,
    blueval,
    greenval,
    redval,
    red,
    green,
    blue,
    hsi2rgb,
    imadjustintensity,
    indexedcolor,
    lut,
    ntsc2rgb,
    rgb2gray,
    rgb2hsi,
    rgb2ntsc,
    rgb2ycbcr,
    uint32color,
    uint32color!,
    ycbcr2rgb,

    # Scaling of intensity
    climdefault,
    float32sc,
    float64sc,
    sc,
    scale,
    scaleinfo,
    scaleminmax,
    scalesigned,
    uint8sc,
    uint16sc,
    uint32sc,

    # algorithms
    ando3,
    ando4,
    ando5,
    backdiffx,
    backdiffy,
    dilate,
    erode,
    opening,
    closing,
    forwarddiffx,
    forwarddiffy,
    gaussian2d,
    imaverage,
    imcomplement,
    imdog,
    imedge,
    imfilter,
    imfilter_fft,
    imfilter_gaussian,
    imfilter_gaussian!,
    imfilter_LoG,
    imgaussiannoise,
    imgradients,
    imlaplacian,
    imlineardiffusion,
    imlog,
    imROF,
    imstretch,
#     imthresh,
    label_components,
    label_components!,
    magnitude,
    magnitude_phase,
    ncc,
    orientation,
    padarray,
    phase,
    prewitt,
    sad,
    sadn,
    sobel,
    ssd,
    ssdn,
    thin_edges,
    thin_edges_subpix,
    thin_edges_nonmaxsup,
    thin_edges_nonmaxsup_subpix,

    # phantoms
    shepp_logan

@deprecate cairoRGB     uint32color!
@deprecate refim        getindexim
@deprecate scaledefault climdefault

if VERSION < v"0.3-"
  __init__()
end

end


text:
module HFMod

if length(findin("C:\\Users\\Clinton\\Dropbox\\Projects\\SummerProject17",LOAD_PATH)) == 0
  push!(LOAD_PATH,"C:\\Users\\Clinton\\Dropbox\\Projects\\SummerProject17")
end

#=TODO:
1) IV=#
#2) Run specifications

using DataFrames, Distributions, StatsBase, GZip, JLD, Gadfly, JuMP, CTMod, NLopt
#importall CT

#Set this to an optimal value- I find logical cores *(3/4) works well
BLAS.set_num_threads(Int(round((Sys.CPU_CORES*3÷4)))-1)

#if !isdefined(:constDef)
  const constDef = true
  #other constants here

  const DATA_PATH = pwd() * "\\data"
  const DATA_FILE_NAME = "HFPerformance_1999-2017" #pad this by a year since we are interested in flows
  const YEAR_RANGE = 2000:2016
  const DATA_FILE_NAME_DEAD = DATA_FILE_NAME * "_dead"
  const DATE_FORMAT_STR = "yyyymmdd"
  const NOTICE_PERIOD_ADJ = 0 #days to add to the notice (since data are monthly)
  const MIN_PERIODS_FOR_INCLUSION = 12 #number of periods required for inclusion
  const T_TEST_THRESHOLD = 2.0 #this is a parameter for identificaiton if the flows are lagged
  const PERIODS_SAME_AUM_TO_DELETE = 4 #assume 4 consecutive periods of identical aum implies stale data
  const LAGS_NOTICE_TO_CAPTURE = 2
  const LAGS_PERFORMANCE_TO_CAPTURE = 24
  const MIN_ASSETS = 10.0 #10,000,000
  const START_LAGS_FROM_NOTICE = false #start the performance lags from the notice period (not the redemption date)
  const RUN_LM_INTERACTIONS = false #run the interaction regressions (takes a while)
  const RUN_2SLS_INTERACTIONS =  false
  const RESULTS_PATH = pwd() * "\\results"
  const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
  const HEADER_NAME = "header.tex"
  const USE_AGGRESSIVE_GC = true #this enables very agressive GC which helps memory management at cost of performance
  const R_TEST_PATH = "C:\\Users\\Clinton\\Dropbox\\Projects\\Summer17RTester\\Summer17RTester"
  const MIN_MONTHS_NOTICE_IV = 1 #min months notice to be included in IV analysis
  const MIN_MONTHS_NOTICE_LONG = 1 #min months notice to be included in the long file
  const SIG_REDEMPTION_THRESHOLDS = 0.05:0.05:0.1
  const PERFORMANCE_BOUND = 0.5
  const MIN_FLOW = -0.99

#end

###############HELPER FUNCTIONS

#helper function to extract all positive numbers from a string
#IN: A string
#OUT: A vector of floats, possibly empty
function getPosNumsFromString(s::String)
  len::Int = length(s)
  preParseS::Vector{Char} = collect(s)

  #replace non-numerical values with delimitor Ξ
  for i ∈ 1:len
     preParseS[i] = (isnumber(s[i]) || (
      s[i]=='.' && (i>1 && isnumber(s[i-1]) || (
        i<len && isnumber(s[i+1])))))?s[i]:'Ξ'
  end

  #split the tokens, keeping only numberical strings
  sSplit::Vector{String} = split(preParseS,'Ξ',keep=false)

  #parse and return the values as an array of floats
  return broadcast((x::String)->Float64(parse(x)),sSplit)::Vector{Float64}
end

#=helper function which sorts and splits DF on symbol, extracts all positive
numeric values, averages them, and replaces the column
IN: DataFrame, column symbol, a default value, and whetehr to leave NAs
OUT: None explicit. DataFrame modified with numberic column named col,
old column is named col_str=#
function DFStringtoMeanNum!(DF::DataFrame, col::Symbol;
  default::Float64=0.0, leaveNAs::Bool=false, silent::Bool=true)

  oldCol = Symbol(col, "_str")
  rename!(DF, col, oldCol)
  DF[:,col] = default
  sort!(DF, cols=(oldCol))

  by(DF, oldCol) do DFSub::SubDataFrame
    if ~isna(DFSub[1,oldCol])
      parsed::Vector{Float64} = #parse the string
        getPosNumsFromString(DFSub[1,oldCol])
      if length(parsed)>0 #if string wasn't all garbage
        DFSub[:,col] = mean(parsed)
      elseif !silent#Give a heads up given a garbage incentive fee
        println("WARNING: Could not parse ", col, ". Value: ", DFSub[1, oldCol],
          ", Assuming numeric value of ", default)
      end
    elseif leaveNAs
      DFSub[:,col] = NA
    end
  end
end

#takes a vector and writes a string
function vec2String(v::Vector{T}, delim::W = "")::String where
  {T<:Any, W<:Union{AbstractString,Char}}

  return join(String.(v), delim)
end

#calculates the t statistic of the Pearson Correlation (Lowry 2017)
function tStatOfCor(x::Vector{Float64},y::Vector{Float64})
  r::Float64 = cor(x,y)
  return r/sqrt((1-r^2.0)/(length(x)-2))::Float64
end

#this is a type which stores all of the notice lags
struct NoticeLagSymbols
  noticeLags::Vector{Symbol}
  noticeLLags::Vector{Symbol}
  cumNoticeLags::Vector{Symbol}
  cumNoticeLLags::Vector{Symbol}
end

#Constructor for the above type
#IN: Number of lags
#OUT: A container of symbols with the appropraite lags
function NoticeLagSymbols(lags::Int)::NoticeLagSymbols
  noticeLags::Vector{Symbol} = Vector{Symbol}(lags)
  noticeLLags::Vector{Symbol} = similar(noticeLags)
  cumNoticeLags::Vector{Symbol} = similar(noticeLags)
  cumNoticeLLags::Vector{Symbol} = similar(noticeLags)

  for i ∈ 1:lags
    noticeLags[i] = Symbol("noticeLag$i")
    noticeLLags[i] = Symbol("noticeLLag$i")
    cumNoticeLags[i] = Symbol("cumNoticeLag$i")
    cumNoticeLLags[i] = Symbol("cumNoticeLLag$i")
  end

  return NoticeLagSymbols(noticeLags,noticeLLags,cumNoticeLags,
    cumNoticeLLags)
end

#this is a type which stores all of the performance lags
struct PerformanceLagSymbols
  performanceLags::Vector{Symbol} #regular performance lags and value-weighted if aggregated
end

#Constructor for the above type
#IN: Number of lags
#OUT: A container of symbols with the appropraite lags
function PerformanceLagSymbols(lags::Int)::PerformanceLagSymbols

  performanceLags::Vector{Symbol} = Vector{Symbol}(lags)   #preallocate
  for i ∈ 1:lags
    performanceLags[i] = Symbol("performanceLag$i")
  end

  return PerformanceLagSymbols(performanceLags)
end

#this structure holds symbols for cases where flows msut meet a certain threshold
struct RedemptionThresholds
    negLFlows::Vector{Symbol}
    noticePeriod::Vector{Symbol}
    monthsPastNotice::Vector{Symbol}
end

RedemptionThresholds(rng::Range=SIG_REDEMPTION_THRESHOLDS)::RedemptionThresholds =
    RedemptionThresholds([Symbol("negLFlowsT$(Int(t*10000))") for t ∈ collect(rng)],
    [Symbol("noticePeriodT$(Int(t*10000))") for t ∈ collect(rng)],
    [Symbol("monthsPastNoticeT$(Int(t*10000))") for t ∈ collect(rng)])

contains(e::T where T<:CTExpr, s::String)::Bool  = Base.contains("$e",s)
###########Main Methods
#U for unlagged, L for lagged, C for corrected, R for corrected but with
#a lagged default

function preProcessHFData(;flowsLaggingProcedure::Char = 'U', newBinary::Bool = false,
  dropFirstMonth::Bool = false)

  nPer::Int = 0 #numer of periods
  dtFormat::DateFormat = DateFormat(DATE_FORMAT_STR) ## a date format object

  #this allows us to capture a variable amount of lags
  noticeLS::NoticeLagSymbols = NoticeLagSymbols(LAGS_NOTICE_TO_CAPTURE)
  performanceLS::PerformanceLagSymbols = PerformanceLagSymbols(LAGS_PERFORMANCE_TO_CAPTURE)

  #this makes sure we have a binary of the data (Improves load times 3-4x)
  if !isfile("$DATA_PATH\\$DATA_FILE_NAME.jls.gz") || newBinary

    #extract from the zip file, package into dataframe and serialize into zip
    #Saves a lot of time not to have to do this each time
    gHandleAlive::GZipStream = gzopen("$DATA_PATH\\$DATA_FILE_NAME.gz")
    gHandleDead::GZipStream = gzopen("$DATA_PATH\\$DATA_FILE_NAME_DEAD.gz")
    ogStream::GZipStream = gzopen("$DATA_PATH\\$DATA_FILE_NAME.jls.gz","w")
    serialize(ogStream, [readtable(gHandleAlive); readtable(gHandleDead)])

    close(ogStream)
    close(gHandleAlive)
    close(gHandleDead)

  end

  #load the binary
  igStream::GZipStream= gzopen("$DATA_PATH\\$DATA_FILE_NAME.jls.gz")
  HFData::DataFrame = deserialize(igStream)
  close(igStream)

  println("Initial rows: ", size(HFData,1), ", Initial columns: ", size(HFData,2))
  println("Begining pre-processing...")

  #drop some columns
  HFData = HFData[:, [:fund_id, :main_strategy,
    :sub_strategy, :fund_status, :fund_assets, :firm_assets,
    :advance_notice, :fund_assets_as_of,
    :date, :performance, :nav, :assets]]

  #Drop some null valued fields:
  HFData = dropNullsFromDF(HFData, [:date, :performance, :assets,
  :advance_notice])

  HFData = HFData[HFData[:,:assets] .≥ MIN_ASSETS, :]
  HFData[:,:main_strategy] = Symbol.(HFData[:,:main_strategy])

  #drop unlikely outliers (data quality)
  HFData[:,:performance] ./= 100.0
  HFData = HFData[(HFData[:performance].< PERFORMANCE_BOUND) & (
    HFData[:performance].> (-1.0*PERFORMANCE_BOUND)),:]

  HFData[:,:date] = ((x::Int)->Date("$x",dtFormat)).(HFData[:,:date])

  #sort the data by fund id and date
  sort!(HFData,cols = (:fund_id,:date))

  #clean out some funds with stale aum data
  by(HFData, :fund_id) do HFSub::SubDataFrame
    nPer = size(HFSub,1)

    if nPer > PERIODS_SAME_AUM_TO_DELETE #check if this procedure makes sense
      identicalAUMs::Int = 1
      i::Int =  2
      while i ≤ nPer && identicalAUMs < PERIODS_SAME_AUM_TO_DELETE
        identicalAUMs = (HFSub[i,:assets]==HFSub[i-1,:assets])?identicalAUMs+1:1
        i+=1
      end
      if identicalAUMs == PERIODS_SAME_AUM_TO_DELETE #set assets to NA to flag for deletion
        HFSub[:,:assets] = NA
      end
    end
  end

  HFData =  HFData[~isna(HFData[:,:assets]),:] #delete the identified funds

  #delete funds with gaps in their performance history, or funds with a
  #track record less than six months
  HFData[:,:mcnt]=0
  ctr::Int = 0
  by(HFData, :fund_id) do HFSub::SubDataFrame
    mcnt::Int = size(HFSub,1)
    #now check for track record breaks

    if maximum(HFSub[:,:date]) ≠
      Dates.lastdayofmonth(minimum(HFSub[:,:date])+Dates.Month(mcnt)-Dates.Month(1))
      mcnt = 0

      #=if ctr<10
        showall(HFSub[:,[:fund_id,:date,:performance,:main_strategy, :fund_assets, :fund_status]])
      end
      ctr+=1=#
    end
    HFSub[:,:mcnt] = mcnt
  end
#  println("funds dropped: $ctr")
  HFData=HFData[~(HFData[:,:mcnt].<MIN_PERIODS_FOR_INCLUSION),:]

  #get the months of advance notice
  firstMonthCorrection::Int = dropFirstMonth?:1:0
  HFData[:,:monthsNotice] = broadcast((x::Int)->
    max(0,((x+NOTICE_PERIOD_ADJ)÷30)-(firstMonthCorrection))::Int,HFData[:,:advance_notice])####get rid of MAX

  ###The following sections create nessecary columns with placeholders
  #calculate monthly flows and construct the data set
  HFData[:,:flows] = 0.0 #flows
  HFData[:,:redempNotice] = 0.0 #notifications about flows
  HFData[:,:rollingNotice] = 0.0 #notifications about flows
  HFData[:,:rollingNoticeL] = 0.0 #notifications about flows
  HFData[:,:lFlows]  = 0.0 #will hold log flows
  HFData[:,:redempNoticeDt] = Date(1,1,1)

  #create the appropriate columns for lagged flow notice
  for i::Int ∈ 1:LAGS_NOTICE_TO_CAPTURE
    HFData[:,noticeLS.noticeLags[i]] = 0.0
    HFData[:,noticeLS.noticeLLags[i]] = 0.0
    HFData[:,noticeLS.cumNoticeLags[i]] = 0.0
    HFData[:,noticeLS.cumNoticeLLags[i]] = 0.0
  end

  #create columns for performance lag
  for i::Int ∈ 1:LAGS_PERFORMANCE_TO_CAPTURE
    HFData[:,performanceLS.performanceLags[i]] = 0.0
    HFData[:,performanceLS.performanceLags[i]] = NA
  end

  HFData[:,:isNoticePeriod] = 0.0
  HFData[:, :isRedemptionPeriod] = 0.0
  HFData[:,:isNoticePeriodSingle] = 0.0

  ctr = 0
  ###This is the main section for manipulating data by fund. We will
  #adjust the lagging of assets and assign to appropriate lagged variables
  by(HFData, :fund_id) do HFSub::SubDataFrame
    nPer = HFSub[1,:mcnt]
    useLaggingProc::Bool = false

    #now calculate the flows in the lagged and unlagged flows
    flowsUnlagged::Vector{Float64} =Vector{Float64}(nPer-1)
    flowsLagged::Vector{Float64} =Vector{Float64}(nPer-1)
    for i ∈ 2:nPer

        #back out the flows:
        #F_unlagged=A(t)-A(t-1)*(1+r(t))
        flowsUnlagged[i-1] = HFSub[i,:assets] -
        (HFSub[i-1,:assets])*(1.0+HFSub[i,:performance])

        #F_lagged=A(t)-A(t-1)*(1+r(t-1))
        flowsLagged[i-1]  = HFSub[i,:assets] -
        (HFSub[i-1,:assets])*(1.0+HFSub[i-1,:performance])
    end

    if flowsLaggingProcedure == 'U'  #assume all unlagged
      useLaggingProc = false

    elseif flowsLaggingProcedure == 'L' #assume all lagged
      useLaggingProc = true

      #below assumes a comparison process
    elseif flowsLaggingProcedure == 'C' || flowsLaggingProcedure == 'R'

      #get the differencein performance between lagged and unlagged
      perfDelta::Vector{Float64} = HFSub[2:end,:performance] .-
        HFSub[1:(end-1),:performance]

      #check the correlations
      flowsUnlaggedT::Float64 = tStatOfCor(flowsUnlagged,perfDelta)
      flowsLaggedT::Float64 = tStatOfCor(flowsLagged,perfDelta)

      #='C' Logic: Assume unlagged. If unlagged T test for correlation against
      return deltas is <-2 (or some other threshold), switch to lagged unless
      the lagged T stat is greater than 2
      'R' Logic: Same as C but reverse.=#
      if flowsUnlaggedT > -1.0 * T_TEST_THRESHOLD && flowsLaggedT > 1.0 * T_TEST_THRESHOLD
        useLaggingProc = false
      elseif flowsUnlaggedT < -1.0 * T_TEST_THRESHOLD && flowsLaggedT < 1.0 * T_TEST_THRESHOLD
        useLaggingProc = true
      elseif flowsUnlaggedT * flowsLaggedT ≤ 0.0
        if flowsLaggingProcedure == 'C'
          useLaggingProc = false
        elseif flowsLaggingProcedure == 'R'
          useLaggingProc = true
        end
      end
    end

    #now execute the appropriate lagging procedure
    if useLaggingProc
      HFSub[end,:flows] = NA
      HFSub[1:(end-1),:flows] = flowsLagged
      HFSub[1:(end-1),:lFlows] = flowsLagged ./
        (HFSub[2:end,:assets] .-  flowsLagged)
    else
      HFSub[1,:flows] = NA
      HFSub[1,:lFlows] = NA
      HFSub[2:end,:flows] = flowsUnlagged
      HFSub[2:end,:lFlows] = flowsUnlagged ./
        (HFSub[2:end,:assets] .- flowsUnlagged)
    end

    #now we look up the notice period to determine the appropriate month
    #Only will be for negative flows
    #we will make the table two way, at the expense of redundancy
    for i::Int ∈ 1:nPer
      #note notifications of flows does not include the current month's flows

        if !isna(HFSub[i,:flows]) && i > HFSub[i,:monthsNotice]
          noticeStart::Int = i-HFSub[i,:monthsNotice]+1-firstMonthCorrection
          if HFSub[i,:monthsNotice] ≥ MIN_MONTHS_NOTICE_IV & (
              noticeStart ≤ i) & (noticeStart ≥ 1)
              HFSub[noticeStart : i ,:rollingNotice] += HFSub[i,:flows]
          end

            if HFSub[i,:flows] < 0.0
                HFSub[i,:isNoticePeriod] = 1.0
                HFSub[i,:isNoticePeriodSingle] += 1.0
                HFSub[i,:isRedemptionPeriod] = 1.0

                #record the notification on the notification date
                HFSub[i-HFSub[i,:monthsNotice],:redempNotice] = HFSub[i,:flows]
                HFSub[i-HFSub[i,:monthsNotice],:redempNoticeDt] = HFSub[i,:date]

                #record the redemption notificaion lags and performance lags
                if HFSub[i,:monthsNotice] > 0
                    for l::Int ∈ 1:min(LAGS_NOTICE_TO_CAPTURE, HFSub[i,:monthsNotice])
                        HFSub[i-l, noticeLS.noticeLags[l]] = HFSub[i,:flows]
                        HFSub[i-l, noticeLS.noticeLLags[l]] = HFSub[i,:lFlows]
                    end
                    for l::Int ∈ 1:(min(HFSub[i,:monthsNotice],i-1))
                        HFSub[i-l,:isNoticePeriod] = 1.0
                        HFSub[i-l,:isNoticePeriodSingle] += 1.0
                    end
                end
            end
        end
    end

    #this si where we capture the cumulative notifications of redeumptions from prior months
    for i::Int ∈ 2:nPer

        HFSub[i,:rollingNoticeL] = max(MIN_FLOW,
          isna(HFSub[i,:rollingNotice])?NA:(HFSub[i,:rollingNotice] / HFSub[i,:assets]))

    end

    #now assign to the performance lags
    if START_LAGS_FROM_NOTICE #if true, we start from the end of the notice period
      for i::Int ∈ 2:nPer
        if i - 1 - HFSub[i,:monthsNotice] ≥ 1 #(make sure its positive number)
          for l::Int ∈ 1:min(i-1- HFSub[i,:monthsNotice],LAGS_PERFORMANCE_TO_CAPTURE)
            HFSub[i, performanceLS.performanceLags[l]] = HFSub[i - l, :performance]
          end
        end
      end
    else
      for i::Int ∈ 2:nPer
        for l::Int ∈ 1:min(i-1,LAGS_PERFORMANCE_TO_CAPTURE)
          HFSub[i, performanceLS.performanceLags[l]] = HFSub[i - l, :performance]
        end
      end
    end
  end

    for l::Int ∈ 1:LAGS_NOTICE_TO_CAPTURE
      HFData[HFData[:,:monthsNotice].<MIN_MONTHS_NOTICE_IV,noticeLS.cumNoticeLags[l]] = NA
      HFData[HFData[:,:monthsNotice].<MIN_MONTHS_NOTICE_IV,noticeLS.cumNoticeLags[l]] = NA
    end

    HFData[~isna(HFData[:,:lFlows]),:lFlows] =
      ((x::Float64)->(max(x,MIN_FLOW))).(HFData[~isna(HFData[:,:lFlows]),:lFlows])
    #get a seperate column with net negative flows
    HFData[:,:negFlows] = -1.0*HFData[:,:flows]
    HFData[(~isna(HFData[:,:flows])) & (HFData[:,:flows].≥0.0),:negFlows] = NA

    HFData[:,:negFlowsOrZero] = -1.0*HFData[:,:flows]
    HFData[(~isna(HFData[:,:flows])) & (HFData[:,:flows].≥0.0),:negFlowsOrZero] = 0.0

    HFData[:,:negLFlows] = -1.0*HFData[:,:lFlows]
    HFData[(~isna(HFData[:,:lFlows])) & (HFData[:,:lFlows].≥0.0),:negLFlows] = NA

    HFData[:,:negLFlowsOrZero] = -1.0*HFData[:,:lFlows]
    HFData[(~isna(HFData[:,:flows])) & (HFData[:,:flows].≥0.0),:negLFlowsOrZero] = 0.0

    HFData[:,:negRollingNoticeL] = -1.0*HFData[:,:rollingNoticeL]
    HFData[(~isna(HFData[:,:negRollingNoticeL])) & (HFData[:,:negRollingNoticeL].<0.0),
        :negRollingNoticeL] = NA

    HFData[:,:negRollingNoticeLOrZero] = -1.0*HFData[:,:rollingNoticeL]
    HFData[(~isna(HFData[:,:negRollingNoticeLOrZero])) & (HFData[:,:negRollingNoticeLOrZero].<0.0),
        :negRollingNoticeLOrZero] = 0.0

    HFData[:,:negRollingNoticeLOrZeroSingle] = copy(HFData[:,:negRollingNoticeLOrZero])
    HFData[(HFData[:,:isNoticePeriodSingle] .> 1.0)#= | (HFData[:,:monthsNotice] .!= 1)=#,
      :negRollingNoticeLOrZeroSingle] = NA

    #showall(HFData[HFData[:,:rollingNoticeL] .< -100.0,[:fund_id,:date,:lFlows,:assets, :rollingNoticeL, :negRollingNoticeLOrZero,:negRollingNoticeLOrZeroSingle]] )
  #  println(maximum(HFData[:,:rollingNoticeL]))
  #  println(minimum(HFData[:,:rollingNoticeL]))

    #process the notice periods
    #Premise: Data not relevant if no notice period is present
    HFData[HFData[:,:monthsNotice] .< 1,:isNoticePeriod] = NA
    HFData[HFData[:,:monthsNotice] .< 1, :isNoticePeriodSingle] = NA

    #drop overlaps from the non-overalp specification
    HFData[HFData[:,:isNoticePeriodSingle] .> 1.0, :isNoticePeriodSingle] = NA

    HFData[:,:isNoticePeriod1Mo] = 0.0
    HFData[:,:isNoticePeriod1Mo] .= HFData[:,:isNoticePeriod]
    HFData[HFData[:,:monthsNotice] .!= 1,:isNoticePeriod1Mo]  = NA

    #make the lagged flows positive numbers
    for l::Int ∈ 1:LAGS_NOTICE_TO_CAPTURE, sVec::Vector{Symbol} ∈
        [noticeLS.cumNoticeLags, noticeLS.cumNoticeLLags, noticeLS.noticeLags, noticeLS.noticeLLags]

        HFData[~isna(HFData[:,sVec[l]]),sVec[l]] =
            -1.0*HFData[~isna(HFData[:,sVec[l]]),sVec[l]]
    end

    #showall(HFData[3000:3200, [:fund_id, :performance, :assets, :flows, :lFlows, :negLFlows, :negLLFlows, :negRollingNoticeL]])

  #convert string columns to factors
  HFTypes::Vector{Type} = eltypes(HFData) #all types
  HFNames::Vector{String} = names(HFData) #all names
  for i ∈ 1:size(HFData,2) #for all columns
    if HFTypes[i] <: Union{String,Symbol} && !(typeof(HFData[:,i]) <: PooledDataVector)
      pool!(HFData,Symbol(HFNames[i])) #convert to factor
    end
  end

  #create a factor column for dates
  HFData[:,:fDate] = HFData[:,:date]
  pool!(HFData,:fDate)

  #filter out years that we don't want
  #Note we are filtering here so that the endpoints are better covered
  #with respect to the flows
  HFData = HFData[(((x::Date)->Dates.year(x) ∈ collect(YEAR_RANGE)).(HFData[:,:fDate])),:]

  println(unique(HFData[:,:main_strategy]))

  println("Processing complete.")
  println("Ending rows: ", size(HFData,1), ", Ending columns: ", size(HFData,2))

  #showall(HFData[200:300,[:fund_id,:fDate,:monthsNotice,:lFlows,:isNoticePeriod, :isRedemptionPeriod, :isNoticePeriodSingle]])

  #save the binary
  ogStream = gzopen(
    "$DATA_PATH\\$(DATA_FILE_NAME)_$(flowsLaggingProcedure)_$(dropFirstMonth)_clean.jls.gz","w")
  serialize(ogStream, HFData)
  close(ogStream)

end

#This function aggregates HF data on strategyID
#We will do this in the format of the standard wide file, so that
#we can plug it into the long file creation process
function aggregateOnStrategy(HFData::DataFrame; monthsNotice::T = NA, )::DataFrame where T<:Union{NAtype,Int}

    performanceLS::PerformanceLagSymbols =
        PerformanceLagSymbols(LAGS_PERFORMANCE_TO_CAPTURE)

    ###now calculate the aggregates

    #pre-allocate the aggregate dataframe
    strategies::Vector{Symbol} = unique(HFData[:,:main_strategy])
    dates::Vector{Date} =  unique(HFData[:,:fDate])

    #for each date, have a row for each strategy and a row for the total
    HFAgg::DataFrame = HFData[1:((length(strategies)+1) * length(dates)),:]
    HFAgg[:,:] = NA #blank the rows

    ctr::Int = 0
    for s::Symbol ∈ [strategies; :total], d::Date ∈ dates
        ctr += 1
        HFAgg[ctr,:main_strategy] = s
        HFAgg[ctr,:fDate] = d
        HFAgg[ctr,:monthsNotice] = monthsNotice
    end

    HFAgg[:,:var] = 0.0
    #HFAgg[:,:totNegLFlows] = 0.0
    HFAgg[:,:totFlowsOnNotice] = 0.0
    HFAgg[:,:var] = NA
    #HFAgg[:,:totNegLFlows] = NA
    HFAgg[:,:totFlowsOnNotice] = NA

    #get views for each date
    by(HFData,:fDate) do HFSubPart::SubDataFrame
        HFAggSub::SubDataFrame = view(HFAgg, HFAgg[:,:fDate] .== HFSubPart[1,:fDate]::Date)

        stratTable::Dict = Dict(HFAggSub[i,:main_strategy]=>i for i::Int ∈ 1:length(strategies))

        for s::Symbol ∈ strategies
            if haskey(stratTable,s)
                sInd = stratTable[s]
                HFSub = view(HFSubPart, HFSubPart[:,:main_strategy].==s)

                #now fill in the aggregates
                HFAggSub[sInd,:assets] = sum(dropna(HFSub[:,:assets]))
                HFAggSub[sInd,:flows] = sum(dropna(HFSub[:,:flows]))
                HFAggSub[sInd,:lFlows] = HFAggSub[sInd,:flows] / HFAggSub[sInd,:assets]

                HFAggSub[sInd,:negLFlows] = HFAggSub[sInd,:lFlows]≤0.0?-1.0*HFAggSub[sInd,:lFlows]:NA
              #  HFAggSub[sInd,:totNegLFlows] = sum(dropna(HFSub[:,:negFlows]))
                HFAggSub[sInd,:totFlowsOnNotice] =
                    sum(dropna(HFSub[HFSub[:,:monthsNotice].≥1,:negFlows]))
                #use value weighted or equal weighted method to calculate performance
                HFAggSub[sInd,:performance] =
                    sum(dropna(HFSub[:,:performance].*HFSub[:,:assets])) / HFAggSub[sInd,:assets]
                HFAggSub[sInd,:var] =
                    sum(dropna((HFSub[:,:performance].-HFAggSub[sInd,:performance]).^2 .*
                    HFSub[:,:assets])) / HFAggSub[sInd,:assets]

            end
        end
        sInd = findfirst(HFAggSub[:,:main_strategy], :total)

        #now get the totals
        #HFAggSub[sInd,:totNegLFlows] = sum(dropna(HFSubPart[:,:negFlows]))
        HFAggSub[sInd,:totFlowsOnNotice] =
            sum(dropna(HFSubPart[HFSubPart[:,:monthsNotice].≥1,:negFlows]))
        HFAggSub[sInd,:assets] = sum(dropna(HFSubPart[:,:assets]))
        HFAggSub[sInd,:flows] = sum(dropna(HFSubPart[:,:flows]))
        HFAggSub[sInd,:performance] =
            sum(dropna(HFSubPart[:,:assets] .* HFSubPart[:,:performance])) / HFAggSub[sInd,:assets]
        HFAggSub[sInd,:var] =
            sum(dropna((HFSubPart[:,:performance].-HFAggSub[sInd,:performance]).^2 .*
            HFSubPart[:,:assets])) / HFAggSub[sInd,:assets]

    end

    #now get the lags and assign a "fund" id
    strategyID::Int  = 1
    sort!(HFAgg,cols = [:main_strategy,:fDate])
    for s::Symbol ∈ strategies
        HFAggSub::SubDataFrame = view(HFAgg, HFAgg[:,:main_strategy] .== s)
        HFAggSub[:,:fund_id] = strategyID
        for l::Int ∈ 1:LAGS_PERFORMANCE_TO_CAPTURE
            @simd for i  ∈  2:size(HFAggSub,1)
                if i-l>0
                    HFAggSub[i,performanceLS.performanceLags[l]] = HFAggSub[i-l,:performance]
                end
            end
        end
        strategyID += 1
    end

    return HFAgg
end

#returns a table with the total flows
#High performance requires an initial sort
function aggregateOnDate(HFData::DataFrame, lags::Int=5,
  noFirstNoticeMonth::Bool=false)

  nDates::Int = size(unique(HFData[:,:date]),1)
  noticeLS::NoticeLagSymbols = NoticeLagSymbols(lags)

  #preallocate data frame
  HFByDate::DataFrame =
    DataFrame([Date, Int, Float64, Float64, Float64, Float64,Float64],
    [:date, :numFunds, :totAUM, :totFlows, :totRedemp,
    :totRollingRedemp, :totRedempNotice], nDates)
  HFByDate[:,:date] .= unique(HFData[:,:date])

  #allocate the lag columns
  for s ∈ [noticeLS.noticeLags; noticeLS.cumNoticeLags]

    HFByDate[:,s] = 0.0
  end

  #make a lookup table
  dateTbl = Dict(HFByDate[i,:date] => i for i::Int ∈ 1:nDates)

  #might be a better way to do this with aggregates and joins,
  #although this is reasonably fast with a low number of comparisons
  #Basically we are getting totals for each symbol by a bunch of symbols by date
  by(HFData, :date) do HFSub::SubDataFrame
    dateIdx::Int = dateTbl[HFSub[1,:date]]
    applyToDate::Vector{Bool} = HFByDate[:,:date].==HFSub[1,:date]
    HFByDate[dateIdx,:numFunds]=size(unique(HFSub[:,:fund_id]),1)
    HFByDate[dateIdx,:totAUM]=sum(HFSub[~isna(HFSub[:,:assets]),:assets])
    HFByDate[dateIdx,:totFlows]=sum(HFSub[~isna(HFSub[:,:flows]),:flows])
    HFByDate[dateIdx,:totRedemp]=
      sum(HFSub[broadcast((x)->ifelse(~isna(x)&&x<0,true, false),HFSub[:,:flows]),:flows])
    HFByDate[dateIdx,:totRedempNotice]=
      sum(HFSub[~isna(HFSub[:,:redempNotice]),:redempNotice])
    HFByDate[dateIdx,:totRollingRedemp]=
      sum(HFSub[~isna(HFSub[:,:rollingNotice]),:rollingNotice])
    for s ∈ [noticeLS.noticeLags; noticeLS.cumNoticeLags]
      HFByDate[dateIdx,s]= sum(HFSub[~isna(HFSub[:,s]),s])
    end

  end

  return HFByDate

end

#this executes the linear model tests for a NAV valuation effect
function LM(HFData::DataFrame, titleCaption::String, caption::String;
    flowStr::String = "negLFlows", dependentVar=:performance,
    flowStrDesc::String = "Outflows / AUM")

    performanceLS::PerformanceLagSymbols =
        PerformanceLagSymbols(LAGS_PERFORMANCE_TO_CAPTURE)
        #negLLFlows
    XLMSpecs::Vector{CTExpr} = [parse("$flowStr")]
    XLMNames::Vector{Vector{Symbol}} = [[:intercept; :outflowVar]] #note only non-factor names
    FLMCluster::Vector{CTSym} = [nothing]

    #showall(HFData[1:50, [:fund_id, :main_strategy, :negLFlows, :performance]])

    #####put additional specifications here
    push!(XLMSpecs, parse("$flowStr"))
    push!(XLMNames, [:intercept; :outflowVar])
    push!(FLMCluster,:fDate)

    push!(XLMSpecs, parse("$flowStr+fDate"))
    push!(XLMNames, [:intercept; :outflowVar])
    push!(FLMCluster,:main_strategy)

    push!(XLMSpecs, parse("$flowStr+fDate"))
    push!(XLMNames, [:intercept; :outflowVar])
    push!(FLMCluster,:fund_id)

    push!(XLMSpecs, parse("$flowStr + " *
        "$(vec2String(performanceLS.performanceLags, "+"))"))
    push!(XLMNames, [:intercept; :outflowVar; performanceLS.performanceLags])
    push!(FLMCluster,nothing)

    push!(XLMSpecs, parse("$flowStr + " *
        "$(vec2String(performanceLS.performanceLags, "+"))"))
    push!(XLMNames, [:intercept; :outflowVar; performanceLS.performanceLags])
    push!(FLMCluster,:fDate)

    push!(XLMSpecs, parse("$flowStr + " *
    "$(vec2String(performanceLS.performanceLags, "+")) + fDate"))
    push!(XLMNames, [:intercept; :outflowVar; performanceLS.performanceLags])
    push!(FLMCluster,:main_strategy)

    push!(XLMSpecs, parse("$flowStr + " *
    "$(vec2String(performanceLS.performanceLags, "+")) + fDate"))
    push!(XLMNames, [:intercept; :outflowVar; performanceLS.performanceLags])
    push!(FLMCluster,:fund_id)

    modelsLM::Vector{CTLM} = Vector{CTLM}()

    println("Running LM Spec X: $flowStr Y: performance")
    #run the LM Specs
    for i ∈ 1:length(XLMSpecs)
        push!(modelsLM, CTLM(HFData, XLMSpecs[i], fixedEffectsSym = FLMCluster[i],
            dependentVar, XNames=XLMNames[i], YName = :performance))
    end
    ##########LM IO
    descRowNamesLM::Vector{String} = ["Performance Lags", "Date Fixed Effects",
        "Strategy Fixed Effects","Fund Fixed Effects", "N ('000s)"]

    #need to print the descriptive rows
    descContentLM::Vector{Vector{String}} =
        [Vector{String}(length(modelsLM))for i∈ 1:length(descRowNamesLM)]

    #this builds the descriptive rows. There is Probably a better way to do this,
    #but its fairly specific to the project.
    for i ∈ 1:length(XLMSpecs)
        descContentLM[1][i] = "$(contains(XLMSpecs[i], "performanceLag")?"X":"")"
        descContentLM[2][i] = "$(FLMCluster[i]==:fDate||contains(XLMSpecs[i], "fDate")?"X":"")"
        descContentLM[3][i] = "$(FLMCluster[i]==:main_strategy||contains(XLMSpecs[i], "main_strategy")?"X":"")"
        descContentLM[4][i] = "$(FLMCluster[i]==:fund_id?"X":"")"
        descContentLM[5][i] = "$(round(modelsLM[i].N/1000.0,1))"
    end

    TableTextLM::String = texTable(modelsLM, getModWhiteΣ!, [:outflowVar],
        titleCaption = titleCaption,
        colNames = [["($i)" for i::Int ∈ 1:length(modelsLM)]],
        contentRowNames = ["$flowStrDesc"],
        descRowNames = descRowNamesLM,
        descContent = descContentLM,
        decimalDigits = 2,
        columnSepPt = -15,
        scaling = [100.0],
        clearMem = USE_AGGRESSIVE_GC,
        caption = caption)

    return TableTextLM
end

#this specification runs a conventional panel using dummies for redemption and notice periods
#this executes the linear model tests for a NAV valuation effect
function PM(HFData::DataFrame, titleCaption::String, caption::String;
    focalStr= "isNoticePeriod", dependentVar=:performance, cluster=false,
    coefs::T=nothing)::String where T<:Union{Void,Vector{Float64}}

    XPMSpecs::Vector{CTExpr} = [parse("$focalStr")]
    XPMNames::Vector{Vector{Symbol}} = [[:intercept; :isNoticePeriod]] #note only non-factor names
    FPMCluster::Vector{CTSym} = [nothing]
    #####put additional specifications here

    push!(XPMSpecs, parse("$focalStr"))
    push!(XPMNames, [:intercept; :isNoticePeriod])
    push!(FPMCluster,:fDate)

    push!(XPMSpecs, parse("$focalStr + fDate"))
    push!(XPMNames, [:intercept; :isNoticePeriod])
    push!(FPMCluster,:main_strategy)

    push!(XPMSpecs, parse("$focalStr + fDate"))
    push!(XPMNames, [:intercept; :isNoticePeriod])
    push!(FPMCluster,:fund_id)

    push!(XPMSpecs, parse("$focalStr + isRedemptionPeriod"))
    push!(XPMNames, [:intercept; :isNoticePeriod; :isRedemptionPeriod])
    push!(FPMCluster,nothing)

    push!(XPMSpecs, parse("$focalStr + isRedemptionPeriod"))
    push!(XPMNames, [:intercept; :isNoticePeriod; :isRedemptionPeriod])
    push!(FPMCluster,:fDate)

    push!(XPMSpecs, parse("$focalStr + isRedemptionPeriod + fDate"))
    push!(XPMNames, [:intercept; :isNoticePeriod; :isRedemptionPeriod])
    push!(FPMCluster,:main_strategy)

    push!(XPMSpecs, parse("$focalStr + isRedemptionPeriod + fDate"))
    push!(XPMNames, [:intercept; :isNoticePeriod; :isRedemptionPeriod])
    push!(FPMCluster,:fund_id)

    #####

    if !cluster
        FPMCluster[:] .= nothing
    end

    modelsPM::Vector{CTLM} = Vector{CTLM}()
    #run the PM Specs

    println("Running PM Spec X:$focalStr Y:$dependentVar")
    for i ∈ 1:length(XPMSpecs)
        push!(modelsPM, CTLM(HFData, XPMSpecs[i], fixedEffectsSym = FPMCluster[i],
            dependentVar, XNames=XPMNames[i], YName = :performance))
    end

    ##########PM IO
    descRowNamesPM::Vector{String} = ["Date F.E.",
        "Strategy F.E.", "Fund-level F.E","N ('000s)"]

    #need to print the descriptive rows
    descContentPM::Vector{Vector{String}} =
        [Vector{String}(length(modelsPM))for i∈ 1:length(descRowNamesPM)]

    #this builds the descriptive rows. There is Probably a better way to do this,
    #but its fairly specific to the project.
    for i ∈ 1:length(XPMSpecs)
        descContentPM[1][i] = "$(FPMCluster[i]==:fDate||contains(XPMSpecs[i], "fDate")?"X":"")"
        descContentPM[2][i] = "$(FPMCluster[i]==:main_strategy||contains(XPMSpecs[i], "main_strategy")?"X":"")"
        descContentPM[3][i] = "$(FPMCluster[i]==:fund_id?"X":"")"
        descContentPM[4][i] = "$(round(modelsPM[i].N/1000.0,1))"

    end

    TableTextPM::String = texTable(modelsPM, getModWhiteΣ!, [:isNoticePeriod; :isRedemptionPeriod],
      titleCaption = titleCaption,
      colNames = [["($i)" for i∈1:length(modelsPM)]],
      contentRowNames = ["Notice Period", "Redemption Dt"],
      descRowNames = descRowNamesPM,
      descContent = descContentPM,
      decimalDigits = 3,
      columnSepPt = -27,
      scaling = [100.0,100.0],
      #clearMem = USE_AGGRESSIVE_GC,
      caption = caption)

      #get the coefficients we want
      if coefs!=nothing
        for m::CTLM ∈ modelsPM
          push!(coefs,getTerm(m,:isNoticePeriod))
        end
      end

    return TableTextPM
end

function Summary(HFData::DataFrame)::String
    summaryDecimals::Int = 2

    #list all columns for which we will provide summary stats
    summaryCols::Vector{Symbol} =
        [:assets, :performance, :negLFlows, :monthsNotice]
    summaryNumCols::Int = length(summaryCols)

    #make an index for ease of reference
    summaryDict::Dict = Dict(summaryCols[i] => i for i::Int ∈ 1:summaryNumCols)
    summaryDat::Vector{Vector{Float64}} =
        [HFData[~isna(HFData[:,s]),s] for s::Symbol ∈ summaryCols]

    #scale for readability
    summaryDat[summaryDict[:performance]] .*= 100.0
    summaryDat[summaryDict[:negLFlows]] .*= 100.0

    #elliminate zero outflow values from the vector for summary stat purposes
    filter!((x::Float64)->x>0, summaryDat[summaryDict[:negLFlows]])
    summaryRowNames::Vector{String} =
        ["mean", "median", "\$\\sigma\$","Number of Funds",
        "Number of Points","Months (2000-2016)"]

    summaryNumRows::Int = length(summaryRowNames)
    summaryContent::Vector{Vector{String}} =
        [Vector{String}() for i∈ 1:summaryNumRows]

    #holds the width of each data cell in columns
    summaryWidth::Vector{Vector{Int}} =
        [Vector{Int}() for i∈ 1:summaryNumRows]

    #get the mean
    summaryContent[1] = string.([
        round(mean(summaryDat[i]),summaryDecimals) for i::Int ∈ 1:summaryNumCols])
    summaryWidth[1] = ones(summaryNumCols)

    summaryContent[2] = string.([
        round(median(summaryDat[i]),summaryDecimals) for i::Int ∈ 1:summaryNumCols])
    summaryWidth[2] = ones(summaryNumCols)

    #get the standard deviation
    summaryContent[3] = string.([
        round(std(summaryDat[i]),summaryDecimals) for i::Int ∈ 1:summaryNumCols])
    summaryWidth[3] = ones(summaryNumCols)

    summaryContent[4] = ["$(length(unique(HFData[:,:fund_id]))) " *
        "(Currently Active: $(length(unique(HFData[HFData[:,:fund_status].==
        "Active",:fund_id]))))"]
    summaryWidth[4] = [summaryNumCols]

    summaryContent[5] = ["$(length(HFData[:,:fund_id]))"]
    summaryWidth[5] = [summaryNumCols]

    summaryContent[6] = ["$(length(unique(HFData[:,:fDate])))"]
    summaryWidth[6] = [summaryNumCols]

    TableSummary::String = texTable("Summary of Data",
        """Includes alive and dead funds from 2000 to 2016.
        Initial Data set consisted of 1.6mn rows, reduced by preprocessing to
        approximately 380k.  Pre-processing
        filtered out funds with discontinuous data, monthly performance magnitude >50\\%, AUM < 10mn,
        track record < 12 months, 4 or more identical AUMs in a row (stale data), and missing
        relevant fields. Flows calculated via \$AUM_{t}-(1+r_t)AUM_{t-1}\$ and further
        experessed as a percentage of the relevant AUM. A t-test on
        the correlation of return and flows was used to adjust the timing of reported AUM, with
        a critical value of \$2\\sigma\$ implying an adjustment to timing (main results are robust
        to not making this adjustment.) Seven days were added to the notice period to account for
        informal notice, which was then divided by 30 to calculate months' notice.""", #caption
        [["AUM(\\\$mn)", "Performance \\%", "Redemptions \\%AUM", "Months' Notice"]], #colNames
        Vector{String}(),#contentRowNames
        Vector{Matrix{String}}(), #content
        summaryRowNames, #descRowNames
        summaryContent, #descContent
        Vector{String}(), #notes
        widthDescContent = summaryWidth,
        columnSepPt=5)

    plOutRD = plot(HFData, x=:assets,
        Geom.density,
        Guide.title("Assets Historgram"),
        Guide.xlabel("AUM"),Guide.ylabel("Count"))
      draw(PNG("$RESULTS_PATH\\$(DATA_FILE_NAME)_assets.png", 7.5inch, 6inch), plOutRD)

    return TableSummary
end

function joinDataAgg(HFData::DataFrame, HFDataAgg::DataFrame)::DataFrame
    sort!(HFDataAgg, cols=[:fDate, :main_strategy])

    dateSymTable =  Dict((HFDataAgg[i,:fDate], HFDataAgg[i,:main_strategy]) =>
        HFDataAgg[i,:performance] for i::Int∈1:size(HFDataAgg,1))

    HFData[:,:performanceStrat] = 0.0
    HFData[:,:performanceRel] = 0.0

    HFData[:,:performanceStrat] = NA
    HFData[:,:performanceRel] = NA

    warnCtr::Int = 0
    @simd for i ∈ 1:size(HFData,1)
        if haskey(dateSymTable, (HFData[i,:fDate], HFData[i,:main_strategy]))
            #gets the strategy performance and first difference of the performance
            HFData[i,:performanceStrat] = dateSymTable[(HFData[i,:fDate], HFData[i,:main_strategy])]

            #get the benchmark relative performance and first difference
            HFData[i,:performanceRel] = HFData[i,:performance] -
                HFData[i,:performanceStrat]
        elseif warnCtr≤10
            error("Join Failure on row $i. Row Info:\n $(HFLong[i,[:fDate,
                :performance,:monthsNotice,:main_strategy]])")
            #warnCtr += 1
        end
    end

    return HFData
end

#this function gets the estimate of canceled flows
#NOTE: depending on the implementation of dataframe, it may modify HFData
function getCanceledEstimate!(HFData::DataFrame, titleCaption::String, caption::String;
  focalStr::String="isNoticePeriod", dependentVar=:performance)::Vector{String}

  HFDataModel = HFData[~isna(HFData[:,:isNoticePeriod]),:]
  dates::Vector{Date} = unique(HFData[:,:fDate])
  N::Int = length(dates)

  #first, construct the models for the clustering
  XSpecs::Vector{CTExpr} = [parse("$focalStr + fDate")]
  XNames::Vector{Vector{Symbol}} = [[:intercept; :isNoticePeriod]] #note only non-factor names
  FCluster::Vector{CTSym} = [nothing]
  #####put additional specifications here

  push!(XSpecs, parse("$focalStr + fDate"))
  push!(XNames, [:intercept; :isNoticePeriod])
  push!(FCluster,:main_strategy)

  push!(XSpecs, parse("$focalStr + fDate"))
  push!(XNames, [:intercept; :isNoticePeriod])
  push!(FCluster,:fund_id)

  models::Vector{CTLM} = Vector{CTLM}()

  #keep track of the variable names
  M::Int = length(XSpecs)
  μSyms::Vector{Symbol} = [Symbol("μ$i") for i::Int ∈ 1:M]
  σ2Syms::Vector{Symbol} = [Symbol("σ2$i") for i::Int ∈ 1:M]
  YSyms::Vector{Symbol} = [Symbol("Y$i") for i::Int ∈ 1:M]
  thresholdSyms::Vector{Symbol} = [Symbol("threshold$i") for i::Int ∈ 1:M]
  thresholdMinusμSyms::Vector{Symbol} = [Symbol("thresholdMinusμ$i") for i::Int ∈ 1:M]
  totSubmittedSyms::Vector{Symbol} = [Symbol("totSubmitted$i") for i::Int ∈ 1:M]
  canceledSyms::Vector{Symbol} = [Symbol("canceled$i") for i::Int ∈ 1:M]

  #run the model specs

  for i ∈ 1:M
      println("Running Spec $i for redemption calc")
      push!(models, CTLM(HFDataModel, XSpecs[i], fixedEffectsSym = FCluster[i],
          dependentVar, XNames=XNames[i], YName = :performance))

      #check to make sure we didn't drop a row
      if (models[end].N) != size(HFDataModel[:,:fDate],1)
        error("Size Mismatch spec $i- model rows: $(models[end].N), data rows: $(size(HFData[:,:fDate],1))")
      end
      HFDataModel[:,YSyms[i]] = models[end].Y
  end

  HFFlowDat::DataFrame = DataFrame([Date, Float64, Float64, Float64, Float64],
    [:fDate, :totFlowsOnNotice, :totNegFlows, :assets, :assetsModel], N)

  HFFlowDat[:,:fDate] = dates

  #write the variable column names into the dataframe
  for i::Int ∈ 1:M
    HFFlowDat[:,μSyms[i]]  = 0.0
    HFFlowDat[:,σ2Syms[i]] = 0.0
    HFFlowDat[:,thresholdSyms[i]]  = 0.0
    HFFlowDat[:,thresholdMinusμSyms[i]]  = 0.0
    HFFlowDat[:,totSubmittedSyms[i] ] = 0.0
    HFFlowDat[:,canceledSyms[i]] = 0.0

    HFFlowDat[:,μSyms[i]] = NA
    HFFlowDat[:,σ2Syms[i]] = NA
    HFFlowDat[:,thresholdSyms[i]]  = NA
    HFFlowDat[:,thresholdMinusμSyms[i]]  = NA
    HFFlowDat[:,totSubmittedSyms[i]]  = NA
    HFFlowDat[:,canceledSyms[i]] = NA
  end

  #now aggregate by date from the larger set
  by(HFData, :fDate) do HFSub::SubDataFrame
    r = findfirst(HFFlowDat[:,:fDate], HFSub[1,:fDate])

    #aggregate the flows
    HFFlowDat[r,:totNegFlows] = sum(dropna(HFSub[:,:negFlows]))
    HFFlowDat[r,:totFlowsOnNotice] =
        sum(dropna(HFSub[HFSub[:,:monthsNotice].≥1,:negFlows]))
    HFFlowDat[r,:assets] = sum(dropna(HFSub[:,:assets]))

    #this is the case where we are not demeaning the performance and variance
    #beyond the specific date
    for i::Int ∈ 1:M
      if FCluster[i] == nothing
        HFFlowDat[r,μSyms[i]] =
            sum(dropna(HFSub[:,dependentVar].*HFSub[:,:assets])) / HFFlowDat[r,:assets]
        HFFlowDat[r,σ2Syms[i]] =
          sum(dropna((HFSub[:,dependentVar].-HFFlowDat[r,μSyms[i]]).^2 .*
            HFSub[:,:assets])) / HFFlowDat[r,:assets]
      end
    end
  end

  #where we have the demeaned info
  by(HFDataModel, :fDate) do HFSub::SubDataFrame
    r = findfirst(HFFlowDat[:,:fDate], HFSub[1,:fDate])
    HFFlowDat[r,:assetsModel] = sum(dropna(HFSub[:,:assets]))

    #now aggregate the performance

    for i::Int ∈ 1:M
      if FCluster[i] != nothing
        HFFlowDat[r,μSyms[i]] =
            sum(dropna(HFSub[:,Symbol("Y$i")].*HFSub[:,:assets])) / HFFlowDat[r,:assetsModel]
        HFFlowDat[r,σ2Syms[i]] =
          sum(dropna((HFSub[:,Symbol("Y$i")].-HFFlowDat[r,Symbol("μ$i")]).^2 .*
          HFSub[:,:assets])) / HFFlowDat[r,:assetsModel]
      end
    end
  end

  #this is a check to make sure we aggregated correctly
  #=d::Date = HFFlowDat[5,:fDate]
  assetsChk = sum(HFData[HFData[:,:fDate].==d,:assets])
  μChk = sum(HFData[HFData[:,:fDate].==d,:performance].*
    HFData[HFData[:,:fDate].==d,:assets])/assetsChk
  σChk = sum((HFData[HFData[:,:fDate].==d,:performance]-μChk).^2.0 .*
    HFData[HFData[:,:fDate].==d,:assets])/assetsChk

  assetsChkModel = sum(HFDataModel[HFDataModel[:,:fDate].==d,:assets])
  μChkModel = sum(dropna(HFDataModel[HFDataModel[:,:fDate].==d,Symbol("Y3")].*
    HFDataModel[HFDataModel[:,:fDate].==d,:assets]))/assetsChkModel
  σChkModel = sum((HFDataModel[HFDataModel[:,:fDate].==d,Symbol("Y3")]-μChkModel).^2.0 .*
    HFDataModel[HFDataModel[:,:fDate].==d,:assets])/assetsChkModel

  showall(HFFlowDat[5,:])
  println("\nassets Check: ",assetsChk , " μCheck:", μChk, " σCheck: ",σChk)
  println("assetsModel Check: ",assetsChkModel ," μCheckModel:", μChkModel, " σCheckModel: ",σChkModel)=#

  #define some functions for the solver
  const sqrt2Inv::Float64 = 1/(2.0^.5)
  const sqrt2PiInv::Float64 = 1.0/(2.0*π)^0.5

  #need to define special functions here for the solver
  CTCumNorm(z) = 0.5+0.5*erf(z*sqrt2Inv)
  CTNorm(z) = sqrt2PiInv*exp(-z^2.0/2.0)

  for m ∈ 1:M

    println("Solving Values in Model $m: ")

    for r::Int ∈ 1:N
      #print(HFFlowDat[r,:fDate],", ")

        #form the objective function
      μ::Float64 = HFFlowDat[r,μSyms[m]]
      σ::Float64 = HFFlowDat[r,σ2Syms[m]]^0.5
      β::Float64 = getTerm(models[m],Symbol(focalStr))
      function obj(t)
          z = (t-μ)/σ
          ϵ = β + σ*CTNorm(z)/CTCumNorm(z)
          return ϵ*ϵ
      end

      #build the optimizer
      #Algorithms: :LD_LBFGS-fails, :LN_COBYLA-(too slow),
      #  :LD_SLSQP (works!), :LD_TNEWTON (works!), :LD_VAR2/:LD_VAR1 (fails),
      #  PRAXIS (fails/slow), :LD_MMA (fails/slow),
      mod = Model(solver=NLoptSolver(algorithm=:LD_SLSQP))
      JuMP.register(mod, :CTNorm, 1, CTNorm, autodiff=true)
      JuMP.register(mod, :CTCumNorm, 1, CTCumNorm, autodiff=true)
      @variable(mod, t)

      JuMP.register(mod, :obj, 1, obj, autodiff=true)
      @NLobjective(mod, Min, obj(t))
      output = solve(mod)

      #calculate the formulas
      HFFlowDat[r,thresholdSyms[m]] = getvalue(t)
      HFFlowDat[r,thresholdMinusμSyms[m]] = HFFlowDat[r,thresholdSyms[m]] - μ
      #below is total = notice/Φ(z)
      HFFlowDat[r,totSubmittedSyms[m]]  = HFFlowDat[r,:totFlowsOnNotice]/
          CTCumNorm((HFFlowDat[r,thresholdSyms[m]] - μ)/σ) +
          (HFFlowDat[r,:totNegFlows]- HFFlowDat[r,:totFlowsOnNotice])
      HFFlowDat[r,canceledSyms[m]] = HFFlowDat[r,totSubmittedSyms[m]] -
        HFFlowDat[r,:totNegFlows]
    end

  end

  #aggregate to the annual level
  HFFlowDat[:,:year] = (Dates.year).(HFFlowDat[:,:fDate])
  years::Vector{Float64} = unique(HFFlowDat[:,:year])
  HFFlowDatAnn = aggregate(HFFlowDat[:,[:year; :assets; :totFlowsOnNotice;
    :totNegFlows; thresholdSyms; thresholdMinusμSyms;totSubmittedSyms; canceledSyms]],
    :year, [sum,mean])

  #rename the aggregated columns
  for i∈1:M
    HFFlowDatAnn[:,thresholdSyms[i]] = HFFlowDatAnn[:,Symbol("threshold$(i)_mean")]
    HFFlowDatAnn[:,thresholdMinusμSyms[i]] = HFFlowDatAnn[:,Symbol("thresholdMinusμ$(i)_mean")]
    HFFlowDatAnn[:,totSubmittedSyms[i]] = HFFlowDatAnn[:,Symbol("totSubmitted$(i)_sum")]
    HFFlowDatAnn[:,canceledSyms[i]] = HFFlowDatAnn[:,Symbol("canceled$(i)_sum")]
  end

  HFFlowDatAnn[:,:assets] = HFFlowDatAnn[:,:assets_sum]
  HFFlowDatAnn[:,:totFlowsOnNotice] = HFFlowDatAnn[:,:totFlowsOnNotice_sum]
  HFFlowDatAnn[:,:totNegFlows] = HFFlowDatAnn[:,:totNegFlows_sum]

  HFFlowDatAnn = unique(HFFlowDatAnn[:,[:year; :assets; :totFlowsOnNotice;:totNegFlows;
    thresholdSyms; thresholdMinusμSyms; totSubmittedSyms; canceledSyms]])

  #calculate a couple basic ratios
  cancToAssetsSyms::Vector{Symbol} = [Symbol("cancToAssets$i") for i::Int ∈ 1:M]
  cancToTotalSyms::Vector{Symbol} = [Symbol("cancToTotal$i") for i::Int ∈ 1:M]

  for i::Int ∈ 1:M
    HFFlowDatAnn[:,cancToAssetsSyms[i]] =
      HFFlowDatAnn[:,canceledSyms[i]] ./ HFFlowDatAnn[:,:assets]
    HFFlowDatAnn[:,cancToTotalSyms[i]] =
      HFFlowDatAnn[:,canceledSyms[i]] ./ HFFlowDatAnn[:,totSubmittedSyms[i]]
  end

  #Make the tables

  flowTablesText::Vector{String} =  Vector{String}()
  numYears::Int = size(HFFlowDatAnn,1)

  for i::Int ∈ 1:M
    #need tp build a table
    decimals = 3

    #write the appropriate title
    modelTitle::String = " Spec ($i): \$\\beta=$(round(getTerm(models[i],Symbol(focalStr))*100,2))\$"
    modelTitle = titleCaption*modelTitle
    if FCluster[i]==:main_strategy
      modelTitle = modelTitle*" (Includes Strategy F.E.)"
    elseif FCluster[i]==:fund_id
      modelTitle = modelTitle*" (Includes Fund F.E.)"
    end

    flowContent::Vector{Vector{String}} =
        [Vector{String}() for r∈ 1:numYears]

    #the net threshold is only meaningful with absolute, not relative performance
    if dependentVar == :performance
      colNames::Vector{Vector{String}} =
        [["Avg Threshold", "Avg Threshold - \$\\mu\$", "Canceled",
        "Canceled/Planned (\\%)", "Canceled/Assets (\\%)"]]
        for r ∈ 1:numYears
          flowContent[r] = ([string(round(HFFlowDatAnn[r,thresholdSyms[i]]*100.0,decimals)),
            string(round(HFFlowDatAnn[r,thresholdMinusμSyms[i]]*100.0,decimals)),
            string(Int(round(HFFlowDatAnn[r,canceledSyms[i]]))),
            string(round(HFFlowDatAnn[r,cancToTotalSyms[i]]*100.0,decimals)),
            string(round(HFFlowDatAnn[r,cancToAssetsSyms[i]]*100.0,decimals))])
        end
      else
        colNames =
          [["Avg Threshold", "Canceled",
          "Canceled/Planned (\\%)", "Canceled/Assets (\\%)"]]

        for r ∈ 1:numYears
          flowContent[r] = ([string(round(HFFlowDatAnn[r,thresholdSyms[i]]*100.0,decimals)),
            string(Int(round(HFFlowDatAnn[r,canceledSyms[i]]))),
            string(round(HFFlowDatAnn[r,cancToTotalSyms[i]]*100.0,decimals)),
            string(round(HFFlowDatAnn[r,cancToAssetsSyms[i]]*100.0,decimals))])
        end
      end

    rowNames::Vector{String} = string.(HFFlowDatAnn[:,:year])
    numCols = length(colNames[1])

    #println(i)

    push!(flowTablesText, texTable(modelTitle,
        caption, #caption
        colNames, #colNames
        Vector{String}(),#contentRowNames
        Vector{Matrix{String}}(), #content
        rowNames, #descRowNames
        flowContent, #descContent
        Vector{String}(), #notes
        #widthDescContent = summaryWidth,
        columnSepPt=0))

  end

  return flowTablesText

end

#this is the main method for analysis
function analyzeAggregates(;flowsLaggingProcedure::Char = 'U', dropFirstMonth::Bool = false,
    verbose::Bool = true, runSummary::Bool = true, runLM::Bool = true,
    runPM::Bool = true, process=true, save=false )

    #wrapping this up to save time if the file is long
    if process
        #open the pre-processed data file
        igStream::GZipStream =
            gzopen(
              "$DATA_PATH\\$(DATA_FILE_NAME)_$(flowsLaggingProcedure)_$(dropFirstMonth)_clean.jls.gz")
        HFData::DataFrame = deserialize(igStream)
        close(igStream)

        #get the long files
        HFDataAgg::DataFrame = aggregateOnStrategy(HFData, monthsNotice=MIN_MONTHS_NOTICE_LONG) #get wide strategy aggregate

        HFData = joinDataAgg(HFData,HFDataAgg)

        titleCaption::String = string()
        caption::String = string()

        sort!(HFData,cols = (:date, :fund_id))
        HFByDate::DataFrame = aggregateOnDate(HFData, LAGS_NOTICE_TO_CAPTURE) #get date aggregate

        if save
            ogStream::GZipStream = gzopen(
              "$DATA_PATH\\$(DATA_FILE_NAME)_$(flowsLaggingProcedure)_$(dropFirstMonth)_Proc.jls.gz","w")
            serialize(ogStream, [HFData, HFDataAgg,  HFByDate])
            close(ogStream)
        end
    else # if we arn't going to process, load in the processed files
        igStream =
            gzopen("$DATA_PATH\\$(DATA_FILE_NAME)_$(flowsLaggingProcedure)_$(dropFirstMonth)_Proc.jls.gz")
        HFData, HFDataAgg, HFByDate = deserialize(igStream)
        close(igStream)
    end
    ###################Plots
    Gadfly.push_theme(:dark)
    plOut = plot(HFData, x=:negLFlowsOrZero,  #=shape=:variable, size=[.1mm],=#
        Geom.histogram,
        Guide.title("Fund Flows / AUM"),
        Guide.ylabel("Density"),Guide.xlabel("OutFlows"))
      draw(PNG("$RESULTS_PATH\\$(DATA_FILE_NAME)_Flows.png", 7.5inch, 6inch), plOut)

    ##################OLS Spec

    titleCaption = "Performance vs Redemption Size"
    caption = """ This table illustrates the relationship between redemptions atand
        performance. The key takeway is performance on the redemption date is,
        counterintuitivelly, positively correlated with the size of the redemption.
        Date fixed effects are by month. Performance lags includes the prior
        24 months of performance as a covariate. Data is from alive and dead hedge
        funds 2000-2016. Where fixed effects are present, the standard errors are
        clustered either by fund (4 \\& 8), strategy (3 \\& 7), or date (2 \\& 6).
         Outflows were expressed as a percentage of the
        AUM. A positive value for the coefficient implies a positive relationship
        between performance and redemptions."""

    TableTextLM::String = runLM?LM(HFData, titleCaption, caption,
        flowStr = "negLFlowsOrZero", flowStrDesc = "Outflows / AUM"):""

    titleCaption = "Performance vs Cumulative Rolling Redemption Notice Size"
    caption = """This table reflects the relationship between pending notifications
            of redemptions and performance. To avoide bias in the data,
            the first day Date fixed effects are by month. Performance lags includes the prior
            24 months of performance as a covariate. Data is from alive and dead hedge
            funds 2000-2016. Where fixed effects are present, the standard errors are
            clustered either by fund (4 \\& 8), strategy (3 \\& 7), or date (2 \\& 6).
            Outflows were expressed as a percentage of the
            AUM. A positive value for the coefficient implies a positive relationship
            between performance and redemptions."""

    TableTextLMNotice::String = runLM?LM(HFData, titleCaption, caption,
        flowStr = "negRollingNoticeLOrZero", flowStrDesc = "Rolling Notice / AUM"):""

    #########IV Note: See the postIVO backup for the IV code
    ##############Panel Approach

    coefs::Vector{Float64}=Vector{Float64}()

    titleCaption = "Effect of Notice Period on Performance"
    caption = "This specification uses a dummy for the notice period and a dummy
    for the redemption period. In each specification, funds underperform in the
    notice period. Except with fund-level fixed effects,
    hedge funds realize  further losses on the redemption day. Data is from alive and dead hedge
    funds 2000-2016. Interpret a value of 1.0 for the dummy as corresponding
    to a difference in performance of 1.0\\%. Where fixed effects are present, the standard errors are
    clustered either by fund (4 \\& 8), strategy (3 \\& 7), or date (2 \\& 6)."""
    TableTextPM = runPM?PM(HFData,titleCaption,caption, cluster=true, coefs=coefs):""

    titleCaption = "Effect of Notice Period on Performance (BM-Relative)"
    caption = "This is the same specification as the prior table, only
    performance is net of a strategy-specific value-weighted benchmark. Results
    are generally robust. In each specification, funds underperform in the
    notice period. Except with fund-level fixed effects,
    hedge funds realize  further losses on the redemption day. Data is from alive and dead hedge
    funds 2000-2016. Interpret a value of 1.0 for the dummy as corresponding
    to a difference in performance of 1.0\\%. Where fixed effects are present, the standard errors are
    clustered either by fund (4 \\& 8), strategy (3 \\& 7), or date (2 \\& 6)."""
    TableTextPMRel = runPM?PM(HFData,titleCaption,caption, dependentVar=:performanceRel, cluster=true):""

    titleCaption = "Effect of Notice Period on Performance (No-overlap)"
    caption = "A robustness check on underperformance in the notice period.
    Here, overlapping redemptions are dropped from the data. Again, in each
    specification, funds underperform in the notice period. Except with fund-level fixed effects,
    hedge funds realize  further losses on the redemption day. Data is
    from alive and dead hedge funds 2000-2016. Interpret a value of 1.0 for
    the dummy as corresponding to a difference in performance of 1.0\\%. Where
    fixed effects are present, the standard errors are clustered either by fund
    (4 \\& 8), strategy (3 \\& 7), or date (2 \\& 6)."""

    TableTextPMSingle = runPM?PM(HFData,titleCaption,caption,
        focalStr="isNoticePeriodSingle", cluster=true):""

    titleCaption = "Effect of Notice Period on Performance (No-overlap, BM Relative)"
    caption = "This specification drops overlapping redemptions and uses performance
    net of the value-weighted benchmark. Except with fund-level fixed effects,
    hedge funds realize
    further losses on the redemption day. Data is
    from alive and dead hedge funds 2000-2016. Interpret a value of 1.0 for
    the dummy as corresponding to a difference in performance of 1.0\\%. Where
    fixed effects are present, the standard errors are clustered either by fund
    (4 \\& 8), strategy (3 \\& 7), or date (2 \\& 6)."""
    TableTextPMSingleRel = runPM?PM(HFData,titleCaption,caption,
        dependentVar=:performanceRel,focalStr="isNoticePeriodSingle", cluster=true):""

    titleCaption = "Effect of Notice Period on Performance (1mo Notice Only)"
    caption = "This specification drops overlapping redemptions and limits the
    data set to one months' notice. The purpose is to estimate the amount of
    canceled redemptions. Except with fund-level fixed effects,
    hedge funds realize further losses on the redemption day. Data is
    from alive and dead hedge funds 2000-2016. Interpret a value of 1.0 for
    the dummy as corresponding to a difference in performance of 1.0\\%. Where
    fixed effects are present, the standard errors are clustered either by fund
    (4 \\& 8), strategy (3 \\& 7), or date (2 \\& 6)."""
    TableTextPMSingle1Mo = runPM?PM(HFData,titleCaption,caption,
        dependentVar=:performance,focalStr="isNoticePeriod1Mo", cluster=true):""

    ###################Calculate canceled assets

    titleCaption = "Canceled Flows"
    caption = """This table shows the implied total amount of redemptions which were
      submitted to the hedge fund and canceled due to an improvement in
      performance. Data is from alive and dead hedge funds 2000-2016.  \$\\beta\$ is the impact of a
      notice period on performance as detailed in prior tables. ``Avg threshold''
      corresponds to the calculated level of performance at which an investor
      would have canceled their redemption to realize the \$\\beta\$ value, assuming
      a true \$\\beta\$ value of zero. The average difference between the threshold
      and the monthly return is also presented. ``Canceled''
      is the derived estimate of canceled redemptions based on the threshold value.
      ``Canceled/Total'' is the ratio of canceled redemptions to total submitted
      redemptions, while ``Canceled/Assets'' corresponds to canceled redemptions
      to to total fund assets. """

    flowTablesText = getCanceledEstimate!(HFData,titleCaption,caption)

    titleCaption = "Canceled Flows (Benchmark Relative)"
    caption = """This table shows the implied total amount of redemptions which were
      submitted to the hedge fund and canceled due to an improvement in
      performance. Data is from alive and dead hedge funds 2000-2016, while performance
      used is relative to a value-weighted benchmark.  \$\\beta\$ is the impact of a
      notice period on performance as detailed in prior tables. ``Avg threshold''
      corresponds to the calculated level of performance at which an investor
      would have canceled their redemption to realize the \$\\beta\$ value, assuming
      a true \$\\beta\$ value of zero. The average difference between the threshold
      and the monthly return is also presented. ``Canceled''
      is the derived estimate of canceled redemptions based on the threshold value.
      ``Canceled/Total'' is the ratio of canceled redemptions to total submitted
      redemptions, while ``Canceled/Assets'' corresponds to canceled redemptions
      to to total fund assets. """
    flowTablesText =
      [flowTablesText; getCanceledEstimate!(HFData,titleCaption,caption,dependentVar=:performanceRel)]

    ####################Do the Summary Tables and write all table to a file

    TableSummary::String = runSummary?Summary(HFData::DataFrame):""

    writeTables2File([TableSummary; TableTextLM; TableTextLMNotice;
        TableTextPM; TableTextPMRel; TableTextPMSingle; TableTextPMSingleRel;
        TableTextPMSingle1Mo; flowTablesText],
        HEADER_NAME, FOOTER_NAME, path=RESULTS_PATH,
        outName = "RegTables_$(flowsLaggingProcedure)_$(dropFirstMonth).tex") #cosmetic marker block for summary tables

end

function verify2R(;flowsLaggingProcedure::Char = 'U', dropFirstMonth::Bool=false)

    #read the seria
    igStream::GZipStream =
        gzopen("$DATA_PATH\\$(DATA_FILE_NAME)_$(flowsLaggingProcedure)_$(dropFirstMonth)_clean.jls.gz")
    HFData::DataFrame = deserialize(igStream)
    close(igStream)

    #write the csv
    writetable("$R_TEST_PATH\\RHFData.csv",HFData)

end

@time begin
FLPChar = 'C'
dropFirstMonth = false

#preProcessHFData(flowsLaggingProcedure=FLPChar, newBinary=false, dropFirstMonth = dropFirstMonth)
analyzeAggregates(flowsLaggingProcedure=FLPChar, dropFirstMonth = dropFirstMonth,
  runSummary=false, runLM = false, runPM = false, process=false, save = true)
verify2R(flowsLaggingProcedure=FLPChar, dropFirstMonth =  dropFirstMonth)
end

#=@time for c ∈ ['C', 'R', 'L', 'U']
  PreProcessHFData(flowsLaggingProcedure=c, newBinary=true)
  AnalyzeAggregates(flowsLaggingProcedure=c)
end=#

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
ϵ0 = SolidStateDetectors.ϵ0 * u"F / m"
e = SolidStateDetectors.elementary_charge * u"C"

@testset "Infinite Parallel Plate Capacitor" begin
    sim = Simulation{T}(SSD_examples[:InfiniteParallelPlateCapacitor])
    calculate_electric_potential!(sim, convergence_limit = 1e-6, refinement_limits = [0.2, 0.1], verbose = false)
    calculate_weighting_potential!(sim, 1, convergence_limit = 1e-6, refinement_limits = [0.2, 0.1], verbose = false)
    calculate_weighting_potential!(sim, 2, convergence_limit = 1e-6, refinement_limits = [0.2, 0.1], verbose = false)
    # calculate_electric_field!(sim)
    BV_true = (maximum(broadcast(c -> c.potential, sim.detector.contacts)) - minimum(broadcast(c -> c.potential, sim.detector.contacts)))* u"V"
    # Δd = (sim.detector.contacts[2].decomposed_surfaces[1].loc - sim.detector.contacts[1].decomposed_surfaces[1].loc) * u"m"
    Δd = (sim.detector.contacts[2].geometry.origin[1] - sim.detector.contacts[1].geometry.origin[1]) * u"m"
    Δx = (sim.world.intervals[2].right - sim.world.intervals[2].left) * u"m"
    A = Δx * Δx
    V = A * Δd
    E_true = BV_true / Δd
    ϵr = sim.detector.semiconductor.material.ϵ_r
    W_true = uconvert(u"J", (ϵr * ϵ0 * E_true^2 / 2) * A * Δd)
    C_true = uconvert(u"pF", 2 * W_true / (BV_true^2))
    C_Analytical = [C_true -C_true; -C_true C_true]
    C_ssd = calculate_capacitance_matrix(sim)
    @testset "Capacity" begin
        @test all(isapprox.(C_ssd, C_Analytical, rtol = 0.001))
    end
end

@testset "Two Spehres Capacitor" begin
    sim = Simulation{T}(SSD_examples[:TwoSpheresCapacitor])
    simulate!(sim, refinement_limits = [0.2, 0.1, 0.05, 0.02], convergence_limit = 1e-7)
    C_ssd = calculate_capacitance_matrix(sim)

    R_1 = sim.detector.contacts[1].geometry.r * u"m"
    R_2 = sim.detector.contacts[2].geometry.r * u"m"
    d = abs((sim.detector.contacts[1].geometry.origin - sim.detector.contacts[2].geometry.origin)[3]) * u"m"

    # Analytical solutions for the elements of the capacitance matrix:
    function c_ii(r1, r2, d; n = 10)
        r1 = ustrip(uconvert(u"m", r1))
        r2 = ustrip(uconvert(u"m", r2))
        d  = ustrip(uconvert(u"m", d))
        c = 0
        F = 4π*ϵ0 * 1u"m"
        u = acosh((d^2 - r1^2 - r2^2) / (2*r1*r2))
        for i in 0:n
            c += inv(r1*sinh(i*u) + r2*sinh((i + 1)*u))
        end
        uconvert(u"pF", F * r1 * r2 * sinh(u) * c)
    end
    function c_ij(r1, r2, d; n = 10)
        r1 = ustrip(uconvert(u"m", r1))
        r2 = ustrip(uconvert(u"m", r2))
        d  = ustrip(uconvert(u"m", d))
        c = 0
        F = 4π*ϵ0 * 1u"m"
        u = acosh((d^2 - r1^2 - r2^2) / (2*r1*r2))
        for i in 1:n
            c += inv(sinh(i*u))
        end
        uconvert(u"pF", -F * r1 * r2 * sinh(u) * c / d)
    end
    A_c11 = c_ii(R_1, R_2, d)
    A_c22 = c_ii(R_2, R_1, d)
    A_c12 = c_ij(R_1, R_2, d)
    C_Analytical = [A_c11 A_c12; A_c12 A_c22]
    @test all(isapprox.(C_Analytical, C_ssd, rtol = 0.03))
end

struct DummyImpurityDensity{T} <: SolidStateDetectors.AbstractImpurityDensity{T} end

@testset "InfiniteCoaxialCapacitor" begin
    sim_cyl = Simulation{T}(SSD_examples[:InfiniteCoaxialCapacitor])
    sim_car = Simulation{T}(SSD_examples[:InfiniteCoaxialCapacitorCartesianCoords])
    BV_true = (maximum(broadcast(c -> c.potential, sim_cyl.detector.contacts)) - minimum(broadcast(c -> c.potential, sim_cyl.detector.contacts)))* u"V"
    ϵr = sim_cyl.detector.semiconductor.material.ϵ_r
    R1 = sim_cyl.detector.contacts[1].geometry.r[1][1] * u"m"
    R2 = sim_cyl.detector.contacts[2].geometry.r[1][1] * u"m"
    V1 = sim_cyl.detector.contacts[1].potential * u"V"
    V2 = sim_cyl.detector.contacts[2].potential * u"V"
    L = (sim_cyl.world.intervals[3].right - sim_cyl.world.intervals[3].left) * u"m"
    V = π * R2^2 * L
    intV = (R2^2 - R1^2) * π * L

    calculate_electric_potential!(sim_cyl, grid = Grid(sim_cyl, max_tick_distance = (0.3u"cm", 15u"°", 1u"m")),
        convergence_limit = 1e-6, refinement_limits = [0.2, 0.1, 0.05, 0.01], use_nthreads = 1, verbose = false
    )
    calculate_electric_potential!(sim_car, grid = Grid(sim_car, max_tick_distance = (0.3u"cm", 0.3u"cm", 1u"m")),
        convergence_limit = 1e-6, refinement_limits = [0.2, 0.1, 0.05, 0.01], use_nthreads = 1, verbose = false
    )
    calculate_weighting_potential!(sim_cyl, 1, grid = Grid(sim_cyl, max_tick_distance = (0.3u"cm", 15u"°", 1u"m")),
        convergence_limit = 1e-6, refinement_limits = [0.2, 0.1, 0.05, 0.01], use_nthreads = 1, verbose = false
    )
    calculate_weighting_potential!(sim_car, 1, grid = Grid(sim_car, max_tick_distance = (0.3u"cm", 0.3u"cm", 1u"m")),
        convergence_limit = 1e-3, refinement_limits = [0.2], use_nthreads = 1, verbose = false
    )
    calculate_weighting_potential!(sim_cyl, 2, grid = Grid(sim_cyl, max_tick_distance = (0.3u"cm", 15u"°", 1u"m")),
        convergence_limit = 1e-6, refinement_limits = [0.2, 0.1, 0.05, 0.01], use_nthreads = 1, verbose = false
    )
    calculate_weighting_potential!(sim_car, 2, grid = Grid(sim_car, max_tick_distance = (0.3u"cm", 0.3u"cm", 1u"m")),
        convergence_limit = 1e-3, refinement_limits = [0.2], use_nthreads = 1, verbose = false
    )

    C_true = uconvert(u"pF", 2π * ϵr * ϵ0 / log(R2/R1) * L )
    C_Analytical = [C_true -C_true; -C_true C_true]

    C_cyl_ssd = calculate_capacitance_matrix(sim_cyl)
    C_car_ssd = calculate_capacitance_matrix(sim_car)

    @testset "Capacity" begin
        @test all(isapprox.(C_cyl_ssd, C_Analytical, rtol = 0.01))
        @test all(isapprox.(C_car_ssd, C_Analytical, rtol = 0.06))
    end

    # Add impurity density and compare resulting potential to analytic solution
    Q = uconvert(u"C", -1e7 * e * u"cm^-3" * intV)
    Δz = L #(sim_cyl.electric_potential.grid.axes[3].ticks[end] - sim_cyl.electric_potential.grid.axes[3].ticks[1])* u"m"
    # Δz = 0.04f0 * u"m"
    ρ1 = 4Q / ((R2^4 - R1^4)*Δz*2π)
    function ρ_r(r)
        return ρ1 * r^2
    end
    a = ρ1 / (ϵ0 * ϵr)
    c1 = uconvert(u"V", ((a / 16) * (R2^4 - R1^4) - (V2-V1)) / log(R2/R1))
    c2 = uconvert(u"V", V2 + c1 * log(ustrip(R2)) - a * R2^4 / 16)  # unit: Volt
    function potential_analytic(r)
        return uconvert(u"V", a * r^4 / 16 - c1 * log(ustrip(r)) + c2)
    end
    rs_analytic = range(R1, stop = R2, length = 500);
    pot_analytic = map(r -> potential_analytic(r), rs_analytic);

    function SolidStateDetectors.get_impurity_density(cdm::DummyImpurityDensity{T}, pt::CylindricalPoint{T})::T where {T}
        if ustrip(R1) <= pt[1] <= ustrip(R2)
            return -ustrip(ρ1 / e) * pt[1]^2
        else
            return 0
        end
    end

    function SolidStateDetectors.get_impurity_density(cdm::DummyImpurityDensity{T}, pt::CartesianPoint{T})::T where {T}
        SolidStateDetectors.get_impurity_density(cdm, CylindricalPoint(pt))
    end
    sim_cyl.detector = SolidStateDetector(sim_cyl.detector, DummyImpurityDensity{T}());
    sim_car.detector = SolidStateDetector(sim_car.detector, DummyImpurityDensity{T}());

    calculate_electric_potential!(sim_cyl, grid = Grid(sim_cyl, max_tick_distance = (1u"cm", 15u"°", 1u"m")),
        convergence_limit = 1e-7, refinement_limits = [0.2, 0.1, 0.05], use_nthreads = 1, verbose = false
    )
    calculate_electric_potential!(sim_car, grid = Grid(sim_car, max_tick_distance = (0.5u"cm", 0.5u"cm", 1u"m")),
        convergence_limit = 1e-7, refinement_limits = [0.2, 0.1, 0.05, 0.02], use_nthreads = 1, verbose = false
    )

    idxR1 = searchsortedfirst( sim_cyl.electric_potential.grid.axes[1], ustrip(R1));
    rs = sim_cyl.electric_potential.grid.axes[1].ticks[idxR1:end];
    potential_rms_cyl = sqrt(sum((ustrip.(map(r -> potential_analytic(r * u"m"), rs)) .- sim_cyl.electric_potential.data[idxR1:end, 1, 1]).^2)) / (ustrip(BV_true)*length(rs));
    idxR1 = searchsortedfirst( sim_car.electric_potential.grid.axes[1], ustrip(R1));
    rs = sim_car.electric_potential.grid.axes[1].ticks[idxR1:end];
    idy0 = searchsortedfirst( sim_car.electric_potential.grid.axes[2], 0);
    potential_rms_car = sqrt(sum((ustrip.(map(r -> potential_analytic(r * u"m"), rs)) .- sim_car.electric_potential.data[idxR1:end, idy0, 1]).^2))/ (ustrip(BV_true)*length(rs));

    @testset "Potential RMS" begin
        @test potential_rms_cyl < 1e-3
        @test potential_rms_car < 1e-3
    end
end


text:
"""
OED.jl - Optimal Experimental Design for Inverse Problems with State Constraints

Codes used to generate the results presented in

@article{ruthotto2017optimal,
  title={Optimal Experimental Design for Inverse Problems with State Constraints},
  author={Ruthotto, Lars and Chung, Julianne and Chung, Matthias},
  journal={SIAM Journal on Scientific Computing},
  year={2018}
}

Highlights:
- We solve the OED problem involving constraints on the states
                         (bound and linear equality constraints in inversion)
- This code parallelizes ove the training data using jInv's parallelization.
"""
module OED

	using jInv.ForwardShare
	using jInv.Utils
	using jInv.Mesh
	using jInv.InverseSolve
	using LinearOperators
	using MAT

	# descriptions of OED problems
	include("AbstractOEDParam.jl")
	include("oedTypeA.jl")
	include("oedTypeB.jl")
	include("unconstrainedOED.jl")
	include("eqconstrainedOED.jl")
	include("ineqconstrainedOED.jl")
	include("sensMat.jl")

	# optimization routines
	include("ipqp.jl")
	include("regFun.jl")
	include("models.jl")

	# codes to set up OED problem B
	include("rotation2D.jl")
	include("linearInter.jl")
	include("getLinearInterMatrix.jl")
	include("getLinearInterMatrixTranspose.jl")

	# main drivers
	include("runOEDTypeA.jl")
	include("runOEDTypeB.jl")
end