text:
import SHA: sha1

"""
	mutable struct Entry <: JSONEntry
		counts::Dict{Any,Int}
		updated::Int
	end

Keeps statistics about scalar values of a one key and also about items inside a key
- `counts` counts how many times given value appeared (at most max_keys() is held)
- `updated` counts how many times the entry was updated
"""
mutable struct Entry{T} <: JSONEntry
	counts::Dict{T,Int}
	updated::Int
end

is_numeric(s::AbstractString, T::Type{<:Number}) = tryparse(T, s) isa Number
is_real_numeric(s::AbstractString, T::Type{<:Real}) = tryparse(T, s) isa Real

Entry(s::T) where {T<:Number} = Entry(Dict{Number,Int}(),0)
function Entry(s::T) where {T<:AbstractString}
	return Entry(Dict{T,Int}(),0)
end

Base.keys(e::Entry) = e.counts |> keys |> collect |> sort
Base.isempty(e::Entry) = false

types(e::Entry) = e.counts |> keys .|> (typeof) |> unique
unify_types(e::Entry) = promote_type(types(e)...)

is_castable(e, T::Type{<:Number}) = unify_types(e) <: T || unify_types(e) <: AbstractString && e.counts |> keys .|> (x->is_numeric(x, T)) |> all
is_intable(e) = is_castable(e, Int32)
is_floatable(e) = is_castable(e, FloatType)
is_numeric_entry(e, T::Type{<:Number}) = unify_types(e) <: T

shorten_if_str(v) = v
shorten_if_str(v::AbstractString) = length(v) > max_len() ? "$(first(v, max_len()))_$(length(v))_$(bytes2hex(sha1(v)))" : v

# Real is the closest supertype of AbstractFloat and Integer
is_numeric_or_numeric_string(e) = is_intable(e) || is_floatable(e) || is_numeric_entry(e, Real)

"""
	function update!(a::Entry, v)

updates the entry when seeing value v
"""
update!(a::Entry{T}, v::Number; path = "") where {T<:Number} = _update!(a, v)
update!(a::Entry{T}, v::AbstractString; path = "") where {T<:AbstractString} = _update!(a, v)
update!(a::Entry{T}, s::AbstractString; path = "") where {T<:Number} = false

function _update!(a::Entry, v)
	v = shorten_if_str(v)
	if length(keys(a.counts)) < max_keys()
		a.counts[v] = get(a.counts,v,0) + 1
		# it's there because otherwise, after filling the count keys not even the existing ones are updates
	elseif haskey(a.counts, v)
		a.counts[v] += 1
	end
	a.updated += 1
	return true
end

function merge(es::Entry...)
	entry_types = es .|> typeof |> unique
	if promote_type(unify_types.(es)...) <: AbstractFloat
		entry_types = [first(filter(x->unify_types(x) <: AbstractFloat, es))]
	elseif promote_type(unify_types.(es)...) <: Integer
		entry_types = [first(filter(x->unify_types(x) <: Integer, es))]
	end
	updates_merged = sum(updated.(es))
	if length(entry_types) > 1
		multi_entry = MultiEntry([], updates_merged)
		other_entries = filter(e->!(e isa MultiEntry), es)
		# merging here also takes care of max_keys(), which is not done for multi-entries
		merged_entries = [merge(filter(x->x isa t, other_entries)...) for t in entry_types]
		multi_entry.childs = merged_entries
		multi_entry
	else
		counts_merged = merge(+, counts.(es)...)
		if length(counts_merged) > max_keys()
			# originally it was this version, but that lead to instable picking of keys based on version of julia, now I enforce sort withing same value
			# counts_merged_list = sort(collect(counts_merged), by=x->x[2], rev=true)
			counts_merged_list = sort(collect(counts_merged), by=x->(x[2], x[1]), rev=true)
			counts_merged = Dict(counts_merged_list[1:max_keys()])
		end
		Entry(counts_merged, updates_merged)
	end
end

function merge_inplace!(e::Entry, es::Entry...)
	es = [e; es...]
	updates_merged = sum(updated.(es))
	counts_merged = merge(+, counts.(es)...)
	if length(counts_merged) > max_keys()
		counts_merged_list = sort(collect(counts_merged), by=x->x[2], rev=true)
		counts_merged = Dict(counts_merged_list[1:max_keys()])
	end
	e.counts = counts_merged
	e.updated = updates_merged
end

function suggestextractor(e::Entry, settings = NamedTuple(); path::String = "", child_less_than_parent = false)
	t = unify_types(e::Entry)
	t == Any && @error "$(path): JSON does not have a fixed type scheme, quitting"

	for (c, ex) in get(settings, :scalar_extractors, default_scalar_extractor())
		c(e) && return ex(e, child_less_than_parent)
	end
end

# todo: here add argument and decide if it should be full or not
function default_scalar_extractor(small_categorical_max_dimension=100, large_categorical_min_occurences=10, large_categorical_max_dimension=min(10000, max_keys()))
	[
	(e -> length(keys(e)) < small_categorical_max_dimension,
		(e, uniontypes) -> ExtractCategorical(keys(e), uniontypes)),
	(e -> is_intable(e),
		(e, uniontypes) -> extractscalar(Int32, e, uniontypes)),
	(e -> is_floatable(e),
	 	(e, uniontypes) -> extractscalar(FloatType, e, uniontypes)),
	# it's important that condition here would be lower than maxkeys
	(e -> (keys_len = length(keys(e)); large_categorical_min_occurences * keys_len <= sum(values(e.counts)) && keys_len < large_categorical_max_dimension && !is_numeric_or_numeric_string(e)),
		(e, uniontypes) -> ExtractCategorical(keys(e), uniontypes)),
	(e -> true,
		(e, uniontypes) -> extractscalar(unify_types(e), e, uniontypes)),]
end

counts(s::T) where {T<:Entry} = s.counts
Base.hash(e::Entry, h::UInt) = hash((e.counts, e.updated), h)
Base.:(==)(e1::Entry, e2::Entry) = e1.updated === e2.updated && e1.counts == e2.counts
sample_synthetic(e::Entry) = first(keys(e.counts))


text:
"""
    ModelBoundaryConditions([FT=Float64;] U = DefaultBoundaryConditions(FT),
                                          V = DefaultBoundaryConditions(FT),
                                          T = DefaultBoundaryConditions(FT),
                                          S = DefaultBoundaryConditions(FT))

Returns a `NamedTuple` of boundary conditions for a `KPP.Model` with solution
fields `U`, `V`, `T`, `S`.

Example
=======

julia> surface_temperature_flux(model) = cos(model.clock.time)

julia> T_bcs = BoundaryConditions(top = FluxBoundaryCondition(surface_flux))

julis> bcs = ModularKPP.ModelBoundaryConditions(T=T_bcs)
"""
function ModelBoundaryConditions(FT=Float64; U = DefaultBoundaryConditions(FT),
                                             V = DefaultBoundaryConditions(FT),
                                             T = DefaultBoundaryConditions(FT),
                                             S = DefaultBoundaryConditions(FT))
    return (U=U, V=V, T=T, S=S)
end

ext: jl


text:
push!(LOAD_PATH, "/Users/Sam/GitHub/KenyaCoV/src")
# import KenyaCoV
using DataFrames,CSV,MAT,Statistics,LinearAlgebra,Optim,Plots,JLD2
gr()
"""
Get data from 2009 and 2019 --- rescale and get movement data from mobile phone movements
"""
T_2009 = readtable("data/2009_population_estimates.csv")
T_2019 = readtable("data/2019_population_estimates.csv")
T = join(T_2009,T_2019,on = :County)
#Rescale the size of the rural and urban population to match 2019 totals
for i = 1:47
    if !ismissing(T.Rural[i])
        T.Rural[i] = round(Int64,T.Rural[i]*T[:total_size][i]/T[:Total][i])
    else
        T.Rural[i] = 0
    end
    if !ismissing(T.Urban[i])
        T.Urban[i] = round(Int64,T.Urban[i]*T[:total_size][i]/T[:Total][i])
    else
        T.Urban[i] = 0
    end
end
CSV.write("data/combined_population_estimates.csv",T)

#Get movement matrix data
# matread("data/movements.mat")
# mv_data = matopen("data/movements.mat")
# mv_matrix = read(mv_data)["movements"]
# close(mv_data)
# CSV.write("data/combined_population_estimates.csv",T)

@load "data/mv_matrix.jld2" mv_matrix;

"""
Idea:
1) Fix the per-capita flux Nairobi <-> Mombassa to 2009 estimates.
2) Construct a gravity model prediction for other moves, varying distance exponent
3) Fit to wider area data using relative entropy
"""
urban_pop_size = [u for u in T.Urban]
x_locations = [x for x in T.x_location]
y_locations = [y for y in T.y_location]
dist_matrix = zeros(47,47) #County level distances
for i = 1:47,j=1:47
    dist_matrix[i,j] = sqrt((x_locations[i] - x_locations[j])^2 + (y_locations[i] - y_locations[j])^2)
end
ind_mombasa = findfirst(T.County .== "Mombasa")
ind_nairobi = findfirst(T.County .== "Nairobi")
normed_movements = copy(mv_matrix) #Data at the wider spatial scale
for j = 1:20
    normed_movements[:,j] = LinearAlgebra.normalize(normed_movements[:,j],1)
end
NtoMmoves = normed_movements[12,4] #Nairobi to Mombasa
MtoNmoves = normed_movements[4,12] #Mombassa to Nairobi

function prediction_moves(dist_matrix,N,α)
    P = similar(dist_matrix)
    d1,d2 = size(dist_matrix)
    for j = 1:d2
        for i = 1:d1
            if j != i
                P[i,j] = N[i]/(dist_matrix[i,j]^α)
            else
                P[i,j] = 0.
            end
        end
        P[:,j] = LinearAlgebra.normalize(P[:,j],1)
        if j == ind_nairobi
            P[ind_mombasa,j] = NtoMmoves
            x = sum(P[:,j]) - NtoMmoves
            for i = 1:d1
                if i != ind_mombasa
                    P[i,j] = P[i,j]*(1-NtoMmoves)/x
                end
            end
        end
        if j == ind_mombasa
            P[ind_nairobi,j] = MtoNmoves
            x = sum(P[:,j]) - MtoNmoves
            for i = 1:d1
                if i != ind_nairobi
                    P[i,j] = P[i,j]*(1-MtoNmoves)/x
                end
            end
        end
    end
    return P
end

function aggregate_movements(P,T)
    n = length(unique(T.wider_area))
    P̂ = zeros(n,n)
    urban_population_size_wa = [sum(T.Urban[findall(T.wider_area .== i)]) for i = 1:20]
    for i = 1:n,j=1:n
        in_i = findall(T.wider_area .== i)
        in_j = findall(T.wider_area .== j)
        for l in in_i,k in in_j
            P̂[i,j] += (T.Urban[k]/urban_population_size_wa[j])*P[l,k]
        end
    end
    P̂[4,12] = MtoNmoves #Set to correct amount
    P̂[12,4] = NtoMmoves
    return P̂
end

function rel_entropy(q,p)#rel entropy of p wrt q
    sum = 0.
    for i=1:length(q)
        if p[i] > 0. && q[i] >0.
            sum += q[i]*log(q[i]/p[i])
        end
    end
    return sum
end

function gravity_model_rel_entropy(mv_data,T,dist_matrix,N,α)
    d1,d2 = size(mv_data) #<--- base distributions
    P = prediction_moves(dist_matrix,N,α) #<--- target distributions
    P̂ = aggregate_movements(P,T)
    sum_rel_entropy = 0.
    for j = 1:d1
        sum_rel_entropy += rel_entropy(P̂[:,j],mv_data[:,j])
    end
    return sum_rel_entropy
end

"""
Optimise the distance exponent
"""
function opt_func(θ)
    α = θ[1]
    error = gravity_model_rel_entropy(normed_movements,T,dist_matrix,T.Urban,α)
end
fit = optimize(opt_func,[0.6],Newton())
α_opt = Optim.minimizer(fit)
rel_min = Optim.minimum(fit)
P_opt = prediction_moves(dist_matrix,urban_pop_size,α_opt[1])

@save "data/optimal_movement_matrix.jld2" P_opt
P[28,30] - NtoMmoves
P[30,28] - MtoNmoves
"""
Plot the total relative entropy
"""
α_range = 0.:0.01:5.
rel_entropy_range = [opt_func([α]) for α in α_range]
plt = plot(α_range,rel_entropy_range,
    lab = "Total relative entropy",lw =3,
    xlabel = "distance exponent",
    title = "Fit of gravity model to aggregated movement data")
scatter!(plt,[α_opt],[rel_min],ms=5, lab = "minimum at alpha = $(round(α_opt[1],digits = 2))")

savefig(plt,"data/rel_entropy.png")
"""
Estimate time spent away by county
"""
median_length_of_trip = 5.
ρ_widerarea = [sum(mv_matrix[:,j])*median_length_of_trip/(1000*365) for j = 1:20]
ρ_county = ones(47)
for i = 1:47
    wider_area = T.wider_area[i]
    ρ_county[i] = ρ_widerarea[wider_area]
end

"""
Calculate optimal transition matrix
"""
T_opt = similar(P_opt)
d1,d2 = size(T_opt)
for i = 1:d1,j=1:d2
    if i == j
        T_opt[i,j] = 1-ρ_county[j]
    else
        T_opt[i,j] = ρ_county[j]*P_opt[i,j]
    end
end
names = T.County
D = Dict("movement_mat" => T_opt,"County_names"=>names)
matwrite("data/mixing_matrix.mat",D)

@save "data/optimal_transition_matrix.jld2" T_opt ρ_county


text:
using ShortStrings, SortingLab, FastGroupBy, SortingAlgorithms

import FastGroupBy: _fastby!
function _fastby!(
    fn::Function,
    byvec::AbstractVector{ShortString{T}},
    valvec::AbstractVector{S}) where {T, S}

    # make structure
    bysc = ShortStrings.size_content.(byvec)
    idx = SortingLab.fsortperm(bysc)
    # sort!(bv, by = x->x[1].size_content, alg = RadixSort)

    byvecv = byvec[idx]
    valvecv = valvec[idx]

    lastby = byvecv[1]
    res = Dict{ShortString{T}, eltype(fn(valvec[1:1]))}()
    start = 1
    @inbounds for i = 2:length(byvec)
        newby = byvecv[i]
        if lastby != newby
            res[lastby] = fn(@view(valvecv[start:i-1]))
            lastby = newby
            start = i
        end
    end

    res[byvecv[end]] = fn(@view(valvecv[start:end]))
    (keys(res) |> collect, values(res))
end

if false
    @time byvec = rand("id".*dec.(1:100,3), 10_000_000) .|> ShortString7;
    @time valvec = rand(length(byvec));
    fn = sum
    T = UInt64
    outType = Float64
    @time _fastby!(sum, byvec, valvec);
    @time fastby(sum, byvec, valvec);

    a = rand([randstring(rand(1:8))  for i = 1:100_000] .|> ShortString15, 10_000_000)
end


text:
#!/usr/bin/env julia

# Accept optional environment-based arguments
host = get(ENV, "PKGSERVER_HOST", "127.0.0.1")
port = parse(Int, get(ENV, "PKGSERVER_PORT", "8000"))
storage_servers = strip.(split(get(ENV, "PKGSERVER_STORAGESERVERS", "https://pkg.julialang.org"), ","))

using PkgServer
empty!(PkgServer.STORAGE_SERVERS)
append!(PkgServer.STORAGE_SERVERS, storage_servers)
PkgServer.start(;host=host, port=port)


text:
using NQCDynamics: RingPolymers
using Distributions: Distributions

struct HarmonicRingPolymer{T,D}
    normal_mode_distribution::D
    normal_mode_transformation::Matrix{T}
end

"""
    HarmonicRingPolymer(ω, β, m, n_beads)

Ring polymer position distribution in a 1D harmonic potential
"""
function HarmonicRingPolymer(ω, β, m, n_beads)
    βₙ = β / n_beads
    ωₙ = 1 / βₙ
    ωₖ = RingPolymers.get_matsubara_frequencies(n_beads, ωₙ)
    σ⁻¹ = sqrt.(βₙ*m .* (ω^2 .+ ωₖ.^2))
    σ = 1 ./ σ⁻¹

    normal_mode_distribution = Distributions.MvNormal(σ)
    normal_mode_transformation = RingPolymers.get_normal_mode_transformation(n_beads)
    HarmonicRingPolymer(normal_mode_distribution, normal_mode_transformation)
end

function Base.rand(rng::AbstractRNG, s::HarmonicRingPolymer)
    R = rand(rng, s.normal_mode_distribution)
    return s.normal_mode_transformation * R
end

function select_item(x::Vector{<:HarmonicRingPolymer{T}}, ::Integer, size::NTuple) where {T}
    length(size) == 3 || throw(error("`size` $size must have 3 dimensions."))
    size[2] == length(x) || throw(error("Sample size does not match distribution."))
    size[3] == length(x[1].normal_mode_distribution) || throw(error("Sample size does not match distribution."))

    out = Array{T,3}(undef, size)
    @views for i in axes(out, 2)
        for j in axes(out, 1)
            out[j,i,:] .= rand(x[i])
        end
    end

    return out
end


text:
using SemiclassicalOrthogonalPolynomials, ClassicalOrthogonalPolynomials, Test
import ClassicalOrthogonalPolynomials: orthogonalityweight, Weighted, associated

@testset "Two Band" begin
    @testset "TwoBandWeight" begin
        w = TwoBandWeight(0.6,1,2,3)
        @test w[0.7] == w[-0.7] == 0.7^6 * (0.7^2-0.6^2)^2 * (1-0.7^2)
        @test_throws BoundsError w[0.4]
        @test_throws BoundsError w[-0.4]
        @test_throws BoundsError w[1.1]
        @test copy(w) == w
        @test_broken sum(w)
    end
    @testset "Chebyshev case" begin
        ρ = 0.5
        T = TwoBandJacobi(ρ, -1/2, -1/2, 1/2)
        @test T[0.6,1:5] ≈ [T[0.6,k] for k=1:5]

        U = TwoBandJacobi(ρ, 1/2, 1/2, -1/2)
        @test U[0.6,1:5] ≈ [U[0.6,k] for k=1:5]

        @test copy(T) == T
        @test U ≠ T
        @test orthogonalityweight(T) == TwoBandWeight(ρ, -1/2, -1/2, 1/2)

        # bug
        @test !issymmetric(jacobimatrix(T)[1:10,1:10])
    end

    @testset "Hilbert" begin
        ρ = 0.5
        w = TwoBandWeight(ρ, -1/2, -1/2, 1/2)
        x = axes(w,1)
        H = inv.(x .- x')
        @test iszero(H*w)
        @test sum(w) == π/2

        T = TwoBandJacobi(ρ, -1/2, -1/2, 1/2)
        Q = associated(T)
        @test (Q[0.6,1:100]' * (Q[:,Base.OneTo(100)] \ exp.(x))) ≈ exp(0.6)
        @test (Q[-0.6,1:100]' * (Q[:,Base.OneTo(100)] \ exp.(x))) ≈ exp(-0.6)

        @test (Q * (Q \ exp.(x)))[0.6] ≈ exp(0.6)
        @test_broken Q \ (H * Weighted(T)) # need to deal with Hcat

        @test_broken H*TwoBandWeight(ρ, 1/2, 1/2, -1/2)
    end
end


text:
immutable Poisson <: DiscreteUnivariateDistribution
    lambda::Float64
    function Poisson(l::Real)
    	if l > 0.0
    		new(float64(l))
    	else
    		error("lambda must be positive")
    	end
    end
end

Poisson() = Poisson(1.0)

@_jl_dist_1p Poisson pois

insupport(d::Poisson, x::Number) = isinteger(x) && 0.0 <= x

function logpdf(d::Poisson, mu::Real, y::Real)
	return ccall((:dpois, Rmath),
		         Float64,
		         (Float64, Float64, Int32),
		         y, mu, 1)
end

mean(d::Poisson) = d.lambda

function mgf(d::Poisson, t::Real)
    l = d.lambda
    return exp(l * (exp(t) - 1.0))
end

function cf(d::Poisson, t::Real)
    l = d.lambda
    return exp(l * (exp(im * t) - 1.0))
end

var(d::Poisson) = d.lambda

# GLM Methods

function devresid(d::Poisson, y::Real, mu::Real, wt::Real)
	return 2.0 * wt * (xlogxdmu(y, mu) - (y - mu))
end

function devresid(d::Poisson, y::Vector{Float64},
	              mu::Vector{Float64}, wt::Vector{Float64})
    [2.0 * wt[i] * (xlogxdmu(y[i], mu[i]) - (y[i] - mu[i])) for i in 1:length(y)]
end

mustart(d::Poisson, y::Real, wt::Real) = y + 0.1

var(d::Poisson, mu::Real) = mu


text:
mutable struct fLDA <: TopicModel
	"Filtered latent Dirichlet allocation."

	K::Int
	M::Int
	V::Int
	N::Vector{Int}
	C::Vector{Int}
	corp::Corpus
	topics::VectorList{Int}
	eta::Float64
	alpha::Vector{Float64}
	kappa::Vector{Float64}
	kappa_old::Vector{Float64}
	kappa_temp::Vector{Float64}
	beta::Matrix{Float64}
	beta_old::Matrix{Float64}
	beta_temp::Matrix{Float64}
	Elogtheta::VectorList{Float64}
	Elogtheta_old::VectorList{Float64}
	gamma::VectorList{Float64}
	tau::VectorList{Float64}
	tau_old::VectorList{Float64}
	phi::MatrixList{Float64}
	elbo::Float64

	function fLDA(corp::Corpus, K::Integer)
		check_corp(corp)
		K > 0 || throw(ArgumentError("Number of topics must be a positive integer."))

		M, V, U = size(corp)
		N = [length(doc) for doc in corp]
		C = [size(doc) for doc in corp]

		topics = [collect(1:V) for _ in 1:K]

		eta = 0.5
		alpha = ones(K)
		kappa = rand(Dirichlet(V, 1.0))
		kappa_old = copy(kappa)
		kappa_temp = zeros(V)
		beta = rand(Dirichlet(V, 1.0), K)'
		beta_old = copy(beta)
		beta_temp = zeros(K, V)
		Elogtheta = [-eulergamma * ones(K) .- digamma(K) for _ in 1:M]
		Elogtheta_old = deepcopy(Elogtheta)
		gamma = [ones(K) for _ in 1:M]
		tau = [fill(eta, N[d]) for d in 1:M]
		tau_old = deepcopy(tau)
		phi = [ones(K, N[d]) / K for d in 1:min(M, 1)]
		elbo = 0

		model = new(K, M, V, N, C, copy(corp), topics, eta, alpha, kappa, kappa_old, kappa_temp, beta, beta_old, beta_temp, Elogtheta, Elogtheta_old, gamma, tau, tau_old, phi, elbo)
		return model
	end
end

function Elogptheta(model::fLDA, d::Int)
	"Compute E_q[log(P(theta))]."

	x = finite(loggamma(sum(model.alpha))) - finite(sum(loggamma.(model.alpha))) + dot(model.alpha .- 1, model.Elogtheta[d])
	return x
end

function Elogpc(model::fLDA, d::Int)
	"Compute E_q[log(P(c))]."

	counts = model.corp[d].counts
	x = log(@boink model.eta^dot(model.tau[d], counts) * (1 - model.eta)^(model.C[d] - dot(model.tau[d], counts)))
	return x
end

function Elogpz(model::fLDA, d::Int)
	"Compute E_q[log(P(z))]."

	counts = model.corp[d].counts
	x = dot(model.phi[1] * counts, model.Elogtheta[d])
	return x
end

function Elogpw(model::fLDA, d::Int)
	"Compute E_q[log(P(w))]."

	terms, counts = model.corp[d].terms, model.corp[d].counts
	x = sum(model.phi[1] .* log.(@boink model.beta[:,terms]) * (counts .* model.tau[d])) + dot(counts .* (1 .- model.tau[d]), log.(@boink model.kappa[terms]))
	return x
end

function Elogqtheta(model::fLDA, d::Int)
	"Compute E_q[log(q(theta))]."

	x = -entropy(Dirichlet(model.gamma[d]))
	return x
end

function Elogqc(model::fLDA, d::Int)
	"Compute E_q[log(q(c))]."

	counts = model.corp[d].counts
	x = -sum([c * entropy(Bernoulli(model.tau[d][n])) for (n, c) in enumerate(counts)])
	return x
end

function Elogqz(model::fLDA, d::Int)
	"Compute E_q[log(q(z))]."

	counts = model.corp[d].counts
	x = -sum([c * entropy(Categorical(model.phi[1][:,n])) for (n, c) in enumerate(counts)])
	return x
end

function update_elbo!(model::fLDA)
	"Update the evidence lower bound."

	model.elbo = 0
	for d in 1:model.M
		terms = model.corp[d].terms
		model.phi[1] = additive_logistic(model.tau_old[d]' .* log.(@boink model.beta_old[:,terms]) .+ model.Elogtheta_old[d], dims=1)
		model.elbo += Elogptheta(model, d) + Elogpc(model, d) + Elogpz(model, d) + Elogpw(model, d) - Elogqtheta(model, d) - Elogqc(model, d) - Elogqz(model, d)
	end

	return model.elbo
end

function update_eta!(model::fLDA)
	"Update eta."
	"Analytic."

	model.eta = sum([dot(model.tau[d], model.corp[d].counts) for d in 1:model.M]) / sum(model.C)
end

function update_alpha!(model::fLDA, niter::Integer, ntol::Real)
	"Update alpha."
	"Interior-point Newton's method with log-barrier and back-tracking line search."

	Elogtheta_sum = sum([model.Elogtheta[d] for d in 1:model.M])

	nu = model.K
	for _ in 1:niter
		rho = 1.0
		alpha_grad = (nu ./ model.alpha + model.M * (digamma(sum(model.alpha)) .- digamma.(model.alpha))) + Elogtheta_sum
		alpha_invhess_diag = -1 ./ (model.M * trigamma.(model.alpha) + nu ./ model.alpha.^2)
		p = (alpha_grad .- dot(alpha_grad, alpha_invhess_diag)  * (model.M * trigamma(sum(model.alpha)) + sum(alpha_invhess_diag))) .* alpha_invhess_diag

		while minimum(model.alpha - rho * p) < 0
			rho *= 0.5
		end
		@finite model.alpha -= rho * p

		if (rho * norm(alpha_grad) < ntol) & (nu / model.K < ntol)
			break
		end
		nu *= 0.5
	end
	@positive model.alpha
end

function update_kappa!(model::fLDA)
	"Reset kappa variables."
	"Analytic."

	model.kappa_old = model.kappa
	model.kappa = model.kappa_temp ./ sum(model.kappa_temp)
	model.kappa_temp = zeros(model.V)
end

function update_kappa!(model::fLDA, d::Int)
	"Update kappa."
	"Analytic."

	terms, counts = model.corp[d].terms, model.corp[d].counts
	model.kappa_temp[terms] += (1 .- model.tau[d]) .* counts
end

function update_beta!(model::fLDA)
	"Reset beta variables."

	model.beta_old = model.beta
	model.beta = model.beta_temp ./ sum(model.beta_temp, dims=2)
	model.beta_temp = zeros(model.K, model.V)
end

function update_beta!(model::fLDA, d::Int)
	"Update beta."
	"Analytic."

	terms, counts = model.corp[d].terms, model.corp[d].counts
	model.beta_temp[:,terms] += model.phi[1] .* (model.tau[d] .* counts)'
end

function update_Elogtheta!(model::fLDA, d::Int)
	"Update E[log(theta)]."
	"Analytic."

	model.Elogtheta_old[d] = model.Elogtheta[d]
	model.Elogtheta[d] = digamma.(model.gamma[d]) .- digamma(sum(model.gamma[d]))
end

function update_gamma!(model::fLDA, d::Int)
	"Update gamma."
	"Analytic."

	counts = model.corp[d].counts
	@positive model.gamma[d] = model.alpha + model.phi[1] * counts
end

function update_tau!(model::fLDA, d::Int)
	"Update tau."
	"Analytic."

	model.tau_old[d] = model.tau[d]

	terms = model.corp[d].terms
	model.tau[d] = model.eta ./ (@boink model.eta .+ (1 - model.eta) * (model.kappa[terms] .* vec(prod(model.beta[:,terms].^-model.phi[1], dims=1))))
end

function update_phi!(model::fLDA, d::Int)
	"Update phi."
	"Analytic."

	terms = model.corp[d].terms
	model.phi[1] = additive_logistic(model.tau[d]' .* log.(@boink model.beta[:,terms]) .+ model.Elogtheta[d], dims=1)
end

function train!(model::fLDA; iter::Integer=150, tol::Real=1.0, niter::Integer=1000, ntol::Real=1/model.K^2, viter::Integer=10, vtol::Real=1/model.K^2, checkelbo::Real=1, printelbo::Bool=true)
	"Coordinate ascent optimization procedure for filtered latent Dirichlet allocation variational Bayes algorithm."

	check_model(model)
	all([tol, ntol, vtol] .>= 0)										|| throw(ArgumentError("Tolerance parameters must be nonnegative."))
	all([iter, niter, viter] .>= 0)										|| throw(ArgumentError("Iteration parameters must be nonnegative."))
	(isa(checkelbo, Integer) & (checkelbo > 0)) | (checkelbo == Inf)	|| throw(ArgumentError("checkelbo parameter must be a positive integer or Inf."))
	all([isempty(doc) for doc in model.corp]) && (iter = 0)
	(checkelbo <= iter) && update_elbo!(model)

	for k in 1:iter
		for d in 1:model.M
			for _ in 1:viter
				update_phi!(model, d)
				update_tau!(model, d)
				update_gamma!(model, d)
				update_Elogtheta!(model, d)
				if norm(model.Elogtheta[d] - model.Elogtheta_old[d]) < vtol
					break
				end
			end
			update_beta!(model, d)
			update_kappa!(model, d)
		end
		update_beta!(model)
		update_kappa!(model)
		update_alpha!(model, niter, ntol)
		#update_eta!(model)

		if check_elbo!(model, checkelbo, printelbo, k, tol)
			break
		end
	end

	model.topics = [reverse(sortperm(vec(model.beta[i,:]))) for i in 1:model.K]
	nothing
end


text:
"""
fit(model, X::AbstractArray{T},Y::AbstractArray{T})
A Partial Least Squares learning algorithm.
"""
function fitting(model,
             X::AbstractArray{T},
             Y::AbstractArray{T}) where T<:AbstractFloat

    fitresult = trainer(model,X,Y)
    report    = nothing
    cache     = nothing

    return (fitresult,cache,report)
end

"""
transform(model, X::AbstractArray{T})
A Partial Least Squares predictor.
"""
function predicting(model,X::AbstractArray{T}) where T<:AbstractFloat
   predictor(model,X)
end


text:
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
using Bootstrap #External packages, must be added
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")

#Esben's path
#cd("$(homedir())/Documents/GitHub/Thesis/Data")
#path = "$(homedir())/Documents/GitHub/Thesis/Data"

#Skipper's path
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/Data"
#HPC path
#path = "/zhome/9f/d/88706/SpecialeCode/Thesis/Data"
#mainData = loadIndexDataNoDur(path)
#fileName = path*"/Results/IndexData/IndexData"
#mainData = loadConcrete(path)
#fileName = path*"/Results/Concrete/Run1/Concrete_Run"
mainData = loadHousingData(path)
#fileName = path*"/Results/HousingData/HousingData"
#mainData = loadCPUData(path)
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results"
fileName = path*"/HousingData/Housing"

#Reset HPC path
#path = "/zhome/9f/d/88706/SpecialeCode/Thesis/ML"
#cd(path)

mainDataArr = Array(mainData)
halfRows = Int64(floor(size(mainDataArr)[1]/2))
train, test = splitDataIn2(mainDataArr, halfRows, size(mainData)[1])
halfRows = Int64(floor(size(test)[1]/2))
test, valid = splitDataIn2(mainDataArr, halfRows, size(test)[1])

colNames = names(mainData)
trainingData = Array(train)
testData = Array(test)
valiData = Array(valid)
nRows = size(trainingData)[1]
nCols = size(trainingData)[2]

#=
mainData, testData = loadElevatorData(path)
dataSize = size(mainData)
colNames = names(mainData)
trainingData = Array(mainData)
testData = Array(testData)
nRows = size(trainingData)[1]
nCols = size(trainingData)[2]
=#

### STAGE 1 ###
println("STAGE 1 INITIATED")
#Define solve
m = JuMP.Model(solver = GurobiSolver(OutputFlag =0 ))

#Add binary variables variables
@variable(m, 0 <= z[1:nCols] <= 1, Bin )

#Calculate highly correlated matix
HC = cor(trainingData)

#Define objective function
@objective(m, Max, sum(z[i] for i=1:length(z)))

#Define max correlation and constraints
rho = 0.8
for i=1:length(z)
	for j=1:length(z)
		if i != j
			if HC[i,j] >= rho
				@constraint(m,z[i]+z[j] <= 1)
			end
		end
	end
end

#Get solution status
status = solve(m);

#Get objective value
println("Objective value kMax: ", getobjectivevalue(m))
kmax = getobjectivevalue(m)

#Get solution value
#zSolved = getvalue(z)
println("STAGE 1 DONE")

### STAGE 2 ###
println("STAGE 2 INITIATED")
println("Standardizing data")

#Split data by X and y
y = trainingData[:,nCols]
X = trainingData[:,1:nCols-1]
yTest = testData[:,nCols]
XTest = testData[:,1:nCols-1]
yVali = valiData[:,nCols]
XVali = valiData[:,1:nCols-1]

# TRANSFORMATIONS ###
X = expandWithTransformations(X)
XTest = expandWithTransformations(XTest)
XVali = expandWithTransformations(XVali)

# Standardize
standX = zScoreByColumn(X)
standY = zscore(y)
standXTest = zScoreByColumn(XTest)
standYTest = zscore(yTest)
standXVali = zScoreByColumn(XVali)
standYVali = zscore(yVali)

# Output for Lasso regression in R
writedlm(fileName*"XTrain.CSV",standX,",")
writedlm(fileName*"YTrain.CSV",standY,",")
writedlm(fileName*"XTest.CSV",standXTest,",")
writedlm(fileName*"YTest.CSV",standYTest,",")
writedlm(fileName*"XVali.CSV",standXVali,",")
writedlm(fileName*"YVali.CSV",standYVali,",")

#Initialise values for check later
bCols = size(X)[2]
kValue = []
RsquaredValue = []
bSolved = []
bestBeta = []
warmStartBeta = []
warmstartZ = []

startIter = 1

bigM = 100

stage2Model = JuMP.Model(solver = GurobiSolver(TimeLimit = 30))
SSTO = sum((standY[i]-mean(standY))^2 for i=1:length(standY))
amountOfGammas = 2
#Spaced between 0 and half the SSTO since this would then get SSTO*absSumOfBeta which would force everything to 0
gammaArray = log10.(logspace(0, log10.(SSTO), amountOfGammas))

bbdata = NodeData[]

function buildStage2(standX, standY, kmax)
	HC = cor(X)
	HCPairCounter = 0
	#println("Building model")
	#Define parameters and model
	stage2Model = JuMP.Model(solver = GurobiSolver(TimeLimit = 20, OutputFlag = 0));
	gamma = 10

	#Define variables
	@variable(stage2Model, b[1:bCols]) #Beta values

	#Define binary variable (5b)
	@variable(stage2Model, 0 <= z[1:bCols] <= 1, Bin )

	@variable(stage2Model, v[1:bCols]) # auxiliary variables for abs

	@variable(stage2Model, T) #First objective term
	@variable(stage2Model, G) #Second objective term

	#Define objective function (5a)
	@objective(stage2Model, Min, T+G)

	#println("Trying to implement new constraint")
	xSquareExpr = @expression(stage2Model, 0*b[1]^2)
	for l = 1:bCols
		coef = 0
		for j = 1:nRows
		   coef += standX[j,l]^(2)
		end
		append!(xSquareExpr, @expression(stage2Model, coef*b[l]^2))
	end
	#println("Implemented x^2")

	ySquaredExpr = @expression(stage2Model, 0*b[1])
	for j = 1:nRows
		append!(ySquaredExpr,@expression(stage2Model, standY[j,1]^2))
	end
	#println("Implemented y^2")

	simpleBetaExpr = @expression(stage2Model, 0*b[1])
	for l = 1:bCols
		coef = 0
		for j = 1:nRows
			coef += -1*2*standX[j, l]*standY[j]
		end
		append!(simpleBetaExpr, @expression(stage2Model, coef*b[l]))
	end
	#println("Implemented simpleBetaExpr")

	crossBetaExpr = @expression(stage2Model, 0*b[1]*b[2])
	iter = 1
	for l = 1:bCols
		for k = (l + 1):bCols
			coef = 0
			for j = 1:nRows
				coef += 2*standX[j,l]*standX[j,k]
			end
			append!(crossBetaExpr, @expression(stage2Model, coef*b[l,1]*b[k,1]))
		end
		#println("Finished loop $l/$bCols")
	end
	#println("Implemented crossBetaExpr")
	totalExpr = @expression(stage2Model, crossBetaExpr+simpleBetaExpr+xSquareExpr+ySquaredExpr)
	@constraint(stage2Model, quadConst, totalExpr <= T)
	#println("Successfully added quadratic constraints")

	#Define constraints (5c)
	@constraint(stage2Model, conBigMN, -1*b .<= bigM*z) #from 0 to bCols -1
	@constraint(stage2Model, conBigMP,  1*b .<= bigM*z) #from bCols to 2bCols-1

	#Second objective term
	@constraint(stage2Model, 1*b .<= v) #from 2bCols to 3bCols-1
	@constraint(stage2Model, -1*b .<= v) #from 3bCols to 4bCols-1
	oneNorm = sum(v[i] for i=1:bCols)

	#gamma[g]*oneNorm <= G ---> -G <= -gamma[g]*oneNorm --> G >= gamma[g]*oneNorm
	g=5
	@constraint(stage2Model, gammaConstr, 1*oneNorm <= G) #4bCols

	#Define kmax constraint (5d)
	@constraint(stage2Model, kMaxConstr, sum(z[j] for j=1:bCols) <= kmax) #4bCols+1

	#Constraint 5f - can only select one of a pair of highly correlated features
	rho = 0.8
	for k=1:bCols
		for j=1:bCols
			if k != j
				if HC[k,j] >= rho
					HCPairCounter += 1
					@constraint(stage2Model,z[k]+z[j] <= 1) #from 4bCols+1 to (4bCols+1+HCPairCounter)
				end
			end
		end
	end

	#Constraint (5g) - only one transformation allowed (x, x^2, log(x) or sqrt(x))
	for j=1:(nCols-1)
		@constraint(stage2Model, z[j]+z[j+(nCols-1)]+z[j+2*(nCols-1)]+z[j+3*(nCols-1)] <= 1) #from (4bCols+1+HCPairCounter) to (4bCols+1+HCPairCounter+nCols)
	end

	#Adding empty cuts to possibly contain cuts. Can generate a maximum of 6 cuts for
	#each stage 3 and default value is 3 iterations between stage 2 and stage 3,
	#so 18 empty constraints are added + 1 for fixing the solution
	nEmptyCuts = 19
	zOnes = ones(bCols)
	for c = 1:nEmptyCuts
		@constraint(stage2Model, sum(zOnes[i]*z[i] for i=1:bCols) <= bCols) #(4bCols+1+HCPairCounter+nCols) to (4bCols+1+HCPairCounter+nCols+18)
	end

	#addinfocallback(stage2Model, infocallback, when = :Intermediate)

	JuMP.build(stage2Model)

	return internalmodel(stage2Model), HCPairCounter
end

function buildLasso(standX, standY)
	#Define parameters and model
	lassoModel = JuMP.Model(solver = GurobiSolver(TimeLimit = 30, OutputFlag = 0));
	gamma = 10

	#Define variables
	@variable(lassoModel, b[1:bCols]) #Beta values

	@variable(lassoModel, v[1:bCols]) # auxiliary variables for abs

	@variable(lassoModel, T) #First objective term
	@variable(lassoModel, G) #Second objective term

	#Define objective function (5a)
	@objective(lassoModel, Min, T+G)

	#println("Trying to implement new constraint")
	xSquareExpr = @expression(lassoModel, 0*b[1]^2)
	for l = 1:bCols
		coef = 0
		for j = 1:nRows
		   coef += standX[j,l]^(2)
		end
		append!(xSquareExpr, @expression(lassoModel, coef*b[l]^2))
	end
	#println("Implemented x^2")

	ySquaredExpr = @expression(lassoModel, 0*b[1])
	for j = 1:nRows
		append!(ySquaredExpr,@expression(lassoModel, standY[j,1]^2))
	end
	#println("Implemented y^2")

	simpleBetaExpr = @expression(lassoModel, 0*b[1])
	for l = 1:bCols
		coef = 0
		for j = 1:nRows
			coef += -1*2*standX[j, l]*standY[j]
		end
		append!(simpleBetaExpr, @expression(lassoModel, coef*b[l]))
	end
	#println("Implemented simpleBetaExpr")

	crossBetaExpr = @expression(lassoModel, 0*b[1]*b[2])
	iter = 1
	for l = 1:bCols
		for k = (l + 1):bCols
			coef = 0
			for j = 1:nRows
				coef += 2*standX[j,l]*standX[j,k]
			end
			append!(crossBetaExpr, @expression(lassoModel, coef*b[l,1]*b[k,1]))
		end
		#println("Finished loop $l/$bCols")
	end
	#println("Implemented crossBetaExpr")
	totalExpr = @expression(lassoModel, crossBetaExpr+simpleBetaExpr+xSquareExpr+ySquaredExpr)
	@constraint(lassoModel, quadConst, totalExpr <= T)
	#println("Successfully added quadratic constraints")

	#Second objective term
	@constraint(lassoModel, 1*b .<= v) #from 2bCols to 3bCols-1
	@constraint(lassoModel, -1*b .<= v) #from 3bCols to 4bCols-1
	oneNorm = sum(v[i] for i=1:bCols)

	#gamma[g]*oneNorm <= G ---> -G <= -gamma[g]*oneNorm --> G >= gamma[g]*oneNorm
	@constraint(lassoModel, gammaConstr, 1*oneNorm <= G) #4bCols

	JuMP.build(lassoModel)

	return internalmodel(lassoModel)
end

function solveLasso(model)
	bestNBeta = zeros(bCols+1,bCols+4)
	for i=1:bCols+1
		bestNBeta[i,1] = i-1
	end

	println("Solving Lasso for all gamma")
	gammaArray = logspace(-1, 3, 100)
	gamma = 0
	tol = 1e-6

	open(fileName*"LassoLog.csv", "w") do f
		allColNames = expandedColNamesToString(colNames)
		write(f, "RsquaredValue, gamma, nNZ, $allColNames\n")
		for j in 1:length(gammaArray)
			gamma = gammaArray[j]
			changeGammaLasso(model, gamma)

			#Gurobi.writeproblem(model, "testproblem.lp")

			#solve problem
			status = Gurobi.optimize!(model)

			#Get solution
			sol = Gurobi.getsolution(model)

			#Get parameters
			bSolved = sol[1:bCols]

			#Shrink values to 0 if within tolerance "tol"
			for i=1:length(bSolved)
				if bSolved[i] < tol
					if bSolved[i] > -tol
						bSolved[i] = 0
					end
				end
			end

			#Out of sample test
			Rsquared = getRSquared(standXVali,standYVali,bSolved)

			#Count number of non zero elements
			numNZ = countnz(bSolved)

			if Rsquared > bestNBeta[numNZ+1,2]
				bestNBeta[numNZ+1,2] = getRSquared(standXVali,standYVali,bSolved)
				bestNBeta[numNZ+1,3] = getRSquared(standXTest,standYTest,bSolved)
				bestNBeta[numNZ+1,4] = gamma
				bestNBeta[numNZ+1,5:bCols+4] = bSolved
			end

			write(f, "$Rsquared, $gamma, $numNZ, $bSolved\n")
			#println("Rsquared =$Rsquared \t gamma =$j \t n NZ=$numNZ")
		end
	end
	f = open(fileName*"LassoBestK.csv", "w")
	write(f, "k,Rsquared Vali, Rsquared Test,gamma,"*expandedColNamesToString(colNames)*"\n")
	writecsv(f,bestNBeta)
	close(f)
	println("Solve all Lasso problems")
	return bestNBeta
end

lassoM = buildLasso(standX, standY)
bestNLASSO = solveLasso(lassoM)

stage2Model, HCPairCounter = buildStage2(standX,standY, kmax)

# best3Beta[:,1] = Rsquared
# best3Beta[:,2] = kmax
# best3Beta[:,3] = gamma
# best3Beta[:,4:bCols+3] = beta

function solveForAllK(model, kmax)
	best3Beta = zeros(3,bCols+3)
	solArr = zeros(kmax*length(gammaArray),3)
	println("Solving for all k and gamma")
	for i in 1:1:kmax
		#println("Calculating warmstart solution")

		#=
		#Calculate warmstart
		warmStartError = 1e6
		for j in 1:5
			warmStartBetaTemp = gradDecent(standX, standY, 30000, 1e-3, i, HC, bSolved)
			tempError = norm(standY- standX*warmStartBetaTemp)^2
			println("Iteration $j of 5: Warmstart error is:", tempError)
			if tempError < warmStartError
				warmStartError = copy(tempError)
				warmStartBeta = copy(warmStartBetaTemp)
				println("Iteration $j of 5: Changed warmstartBeta")
			end
		end

		println("Warmstart error is:", warmStartError)

		#Set warmstart values
		startVar = zeros(length(warmStartBeta)*3+2)
		warmstartZ = zeros(warmStartBeta)
		for j in 1:bCols
			startVar[j] = warmStartBeta[j]
			startVar[j+2*bCols] = abs(warmStartBeta[j])
			if isequal(warmStartBeta[j],0)
				startVar[j+bCols] = 0
			else
				startVar[j+bCols] = 1
				warmstartZ[j] = 1
			end
		end
		startVar[bCols*3+1] = 10000
		startVar[bCols*3+2] = 10000

		Gurobi.setwarmstart!(m2,startVar)
		=#

		#Set kMax rhs constraint
		curUB = Gurobi.getconstrUB(model) #Get current UBbounds
		curUB[bCols*4+2] = i #Change upperbound in current bound vector
		Gurobi.setconstrUB!(model, curUB) #Push bound vector to model
		Gurobi.updatemodel!(model)

		#Set new Big M
		newBigM = 3#tau*norm(warmStartBeta, Inf)
		changeBigM(model,newBigM)
		for j in 1:length(gammaArray)
			changeGamma(model, gammaArray[j])

			#println("Starting to solve stage 2 model with kMax = $i and gamma = $j")

			#Gurobi.writeproblem(model, "testproblem.lp")

			#Solve Stage 2 model
			status = Gurobi.optimize!(model)
			#println("Objective value: ", Gurobi.getobjval(model))

			sol = Gurobi.getsolution(model)

			#printSolutionToCSV("lars.csv", bbdata)
			#Get solution and calculate R^2
			bSolved = sol[1:bCols]
			zSolved = sol[1+bCols:2*bCols]

			#Out of sample test
			Rsquared = getRSquared(standXVali,standYVali,bSolved)

			if Rsquared > best3Beta[1,1] #largest than 3rd largest Rsquared
				if Rsquared > best3Beta[2,1] #largest than 2nd largest Rsquared
					if Rsquared > best3Beta[3,1] #largest than largest Rsquared
						best3Beta[1,:] = best3Beta[2,:]
						best3Beta[2,:] = best3Beta[3,:]
						best3Beta[3,1] = Rsquared
						best3Beta[3,2] = i
						best3Beta[3,3] = j #potentially store actual gamma value
						best3Beta[3,4:bCols+3] = bSolved
					else
						best3Beta[1,:] = best3Beta[2,:]
						best3Beta[2,1] = Rsquared
						best3Beta[2,2] = i
						best3Beta[2,3] = j #potentially store actual gamma value
						best3Beta[2,4:bCols+3] = bSolved
					end
				else
					best3Beta[1,1] = Rsquared
					best3Beta[1,2] = i
					best3Beta[1,3] = j #potentially store actual gamma value
					best3Beta[1,4:bCols+3] = bSolved
				end
			end

			if any(Rsquared .> RsquaredValue) || isempty(RsquaredValue)
				bestBeta = bSolved
			end
			solArr[Int64((i-1)*length(gammaArray)+j),1] = i
			solArr[Int64((i-1)*length(gammaArray)+j),2] = j
			solArr[Int64((i-1)*length(gammaArray)+j),3] = Rsquared
			push!(kValue, (i-1)*length(gammaArray)+j)
			push!(RsquaredValue, Rsquared)
			#println("Rsquared =$Rsquared\t kMax =$i \t gamma =$j")
		end
		#printNonZeroValues(bSolved)
	end

	return best3Beta, solArr, model
end

function solveAndLogForAllK(model, kmax)
	best3Beta = zeros(3,bCols+3)
	solArr = zeros(kmax*length(gammaArray),3)
	println("Solving for all k and gamma")
	open(fileName*"AALRLog.csv", "w") do f
		allColNames = expandedColNamesToString(colNames)
		write(f, "RsquaredValue,gamma,k, $allColNames\n")
		for i in 1:1:kmax
			#println("Calculating warmstart solution")

			#=
			#Calculate warmstart
			warmStartError = 1e6
			for j in 1:5
				warmStartBetaTemp = gradDecent(standX, standY, 30000, 1e-3, i, HC, bSolved)
				tempError = norm(standY- standX*warmStartBetaTemp)^2
				println("Iteration $j of 5: Warmstart error is:", tempError)
				if tempError < warmStartError
					warmStartError = copy(tempError)
					warmStartBeta = copy(warmStartBetaTemp)
					println("Iteration $j of 5: Changed warmstartBeta")
				end
			end

			println("Warmstart error is:", warmStartError)

			#Set warmstart values
			startVar = zeros(length(warmStartBeta)*3+2)
			warmstartZ = zeros(warmStartBeta)
			for j in 1:bCols
				startVar[j] = warmStartBeta[j]
				startVar[j+2*bCols] = abs(warmStartBeta[j])
				if isequal(warmStartBeta[j],0)
					startVar[j+bCols] = 0
				else
					startVar[j+bCols] = 1
					warmstartZ[j] = 1
				end
			end
			startVar[bCols*3+1] = 10000
			startVar[bCols*3+2] = 10000

			Gurobi.setwarmstart!(m2,startVar)
			=#

			#Set kMax rhs constraint
			curUB = Gurobi.getconstrUB(model) #Get current UBbounds
			curUB[bCols*4+2] = i #Change upperbound in current bound vector
			Gurobi.setconstrUB!(model, curUB) #Push bound vector to model
			Gurobi.updatemodel!(model)

			#Set new Big M
			newBigM = 10#tau*norm(warmStartBeta, Inf)
			changeBigM(model,newBigM)
			for j in 1:length(gammaArray)
				gamma=gammaArray[j]
				changeGamma(model, gammaArray[j])

				#println("Starting to solve stage 2 model with kMax = $i and gamma = $j")

				#Gurobi.writeproblem(model, "testproblem.lp")

				#Solve Stage 2 model
				status = Gurobi.optimize!(model)
				#println("Objective value: ", Gurobi.getobjval(model))

				sol = Gurobi.getsolution(model)

				#printSolutionToCSV("lars.csv", bbdata)
				#Get solution and calculate R^2
				bSolved = sol[1:bCols]
				zSolved = sol[1+bCols:2*bCols]

				#Out of sample test
				Rsquared = getRSquared(standXVali,standYVali,bSolved)

				if Rsquared > best3Beta[1,1] #largest than 3rd largest Rsquared
					if Rsquared > best3Beta[2,1] #largest than 2nd largest Rsquared
						if Rsquared > best3Beta[3,1] #largest than largest Rsquared
							best3Beta[1,:] = best3Beta[2,:]
							best3Beta[2,:] = best3Beta[3,:]
							best3Beta[3,1] = Rsquared
							best3Beta[3,2] = i
							best3Beta[3,3] = j #potentially store actual gamma value
							best3Beta[3,4:bCols+3] = bSolved
						else
							best3Beta[1,:] = best3Beta[2,:]
							best3Beta[2,1] = Rsquared
							best3Beta[2,2] = i
							best3Beta[2,3] = j #potentially store actual gamma value
							best3Beta[2,4:bCols+3] = bSolved
						end
					else
						best3Beta[1,1] = Rsquared
						best3Beta[1,2] = i
						best3Beta[1,3] = j #potentially store actual gamma value
						best3Beta[1,4:bCols+3] = bSolved
					end
				end

				if any(Rsquared .> RsquaredValue) || isempty(RsquaredValue)
					bestBeta = bSolved
				end
				solArr[Int64((i-1)*length(gammaArray)+j),1] = i
				solArr[Int64((i-1)*length(gammaArray)+j),2] = j
				solArr[Int64((i-1)*length(gammaArray)+j),3] = Rsquared
				push!(kValue, (i-1)*length(gammaArray)+j)
				push!(RsquaredValue, Rsquared)
				println("Rsquared =$Rsquared\t kMax =$i \t gamma =$j")
				write(f, "$Rsquared,$gamma,$i,$bSolved\n")
			end
			#printNonZeroValues(bSolved)
		end
	end
	return best3Beta, solArr, model
end

#="""
For each of the three beta sets produced, we will test for significance
and condition number of model to see if cuts are necessary
High or low condition number doesn't mean that one correlation matrix is "better"
than the other. All it means is that variables are more correlated or less.
Whether it's good or not depends on the application.
"""=#

function stageThree(best3Beta, X, Y, allCuts)
	#=Condition Number
	  A high condition number indicates a multicollinearity problem. A condition
	  number greater than 15 is usually taken as evidence of
	  multicollinearity and a condition number greater than 30 is
	  usually an instance of severe multicollinearity
	=#

	bCols = size(X)[2]
	nRows = size(X)[1]
	cuts = Matrix(0, bCols+1)
	rowsPerSample = nRows #All of rows in training data to generate beta estimates, but selected with replacement
	totalSamples = 25 #25 different times we will get a beta estimate
	nBoot = 10000

	#For loop start
	for i = 1:size(best3Beta)[1]
		if signifBoolean[i] == 1 #if the previous solution was already completely significant, skip the bootstrapping
			continue
		end
		bestBeta = best3Beta[i,4:bCols+3]
		bestK = best3Beta[i,2]
		bestGamma = best3Beta[i,3]

		xColumns = []
		bSample = Matrix(totalSamples, bCols)

		bZeros = zeros(bCols)
		for j = 1:bCols
			if !isequal(bestBeta[j],0)
				push!(xColumns, j)
			end
		end

		selectedX = X[:,xColumns]
		condNumber = cond(selectedX)
		if condNumber >= 15
			bZeros[xColumns] = 1
			subsetSize = size(xColumns)[1]
			newCut = [bZeros' subsetSize]
			cuts = [cuts; newCut]
			println("A cut based on Condition number = $condNumber has been created from Beta$i")
		end

		# test significance
		bestZ = zeros(bestBeta)
		for l=1:size(bestBeta)[1]
			if bestBeta[l] != 0
				bestZ[l] = 1
			end
		end

		bZeros = zeros(bCols)
		createBetaDistribution(bSample, X, Y, bestK, totalSamples, rowsPerSample,  bestGamma, allCuts, bestZ) #standX, standY, k, sampleSize, rowsPerSample

		confArray99 = createConfidenceIntervalArray(bSample, nBoot, 0.99)
		confArray95 = createConfidenceIntervalArray(bSample, nBoot, 0.95)
		confArray90 = createConfidenceIntervalArray(bSample, nBoot, 0.90)

		significanceResult = testSignificance(confArray99, confArray95, confArray90, bestBeta)
		significanceResultNONSI = [] # = significanceResult[xColumns]
		subsetSize = size(xColumns)[1]
		for n = 1:size(significanceResult)[1]
			for s = 1:subsetSize
				if significanceResult[n] == 0 && xColumns[s] == n
					push!(significanceResultNONSI,xColumns[s])
					println("Parameter $n is selected, but NOT significant")
				elseif significanceResult[n] > 0 && xColumns[s] == n
					println("Parameter $n is significant with ", significanceResult[n])
				end
			end
		end

		if !isempty(significanceResultNONSI)
			bZeros[significanceResultNONSI] = 1
			subsetSize = size(significanceResultNONSI)[1]
			newCut = [bZeros' subsetSize]
			cuts = [cuts; newCut]
			println("A cut based on parameters being non-significant in Beta$i has been created")
		end
		if isempty(significanceResultNONSI)
			signifBoolean[i] = 1
		end
	end
	return cuts
end

bSample = []
allCuts = []
signifBoolean = zeros(3)

stage2Model, HCPairCounter = buildStage2(standX,standY, 2)

#Gurobi.writeproblem(stage2Model, "testproblem1.lp")

best3Beta, solArr, stage2Model = solveAndLogForAllK(stage2Model, 2)

cuts = stageThree(best3Beta, standX, standY, allCuts)

#cuts = stageThree(best3Beta[3,4:bCols+3], best3Beta[3,2], best3Beta[3,3],
# 				  best3Beta[2,4:bCols+3], best3Beta[2,2], best3Beta[2,3],
#				  best3Beta[1,4:bCols+3], best3Beta[1,2], best3Beta[1,3],
#				  standX, standY, allCuts)
#writedlm("betaOut.CSV", cuts,",")

cutCounter = 0
runCount = 1
cutMatrix = []

while !isempty(cuts)
	#Function to add cuts to problem
	#addCuts(stage2Model, cuts)
	preCutCounter = copy(cutCounter)
	allCuts = vcat(allCuts,cuts)
	cutMatrix = cuts
	cutRows = size(cutMatrix)[1]
	cutCols = size(cutMatrix)[2]
	cutCounter += cutRows
	addCuts(stage2Model, cutMatrix, preCutCounter)

	#Resolve problem
	Gurobi.updatemodel!(stage2Model)
	#Gurobi.writeproblem(stage2Model, "testproblem2.lp")
	#println(stage2Model)
	best3Beta, solArr = solveAndLogForAllK(stage2Model, 2)

	#Stage 3
	cuts = stageThree(best3Beta, standX, standY, allCuts)
	println("Finished iteration $runCount")
	runCount += 1
	if runCount == 4
		break
	end
end

best3Beta[1,1] = getRSquared(standXTest, standYTest, best3Beta[1,4:end])
best3Beta[2,1] = getRSquared(standXTest, standYTest, best3Beta[2,4:end])
best3Beta[3,1] = getRSquared(standXTest, standYTest, best3Beta[3,4:end])

best3Beta = cat(2,signifBoolean,best3Beta)

f = open(fileName*"AALRBestK.csv", "w")
write(f, "Significant,Rsquared,k,gamma,"*expandedColNamesToString(colNames)*"\n")
writecsv(f,best3Beta)
close(f)

bSolvedLinReg = inv(standX'*standX)*standX'*standY
standLinReg = getRSquared(standXTest,standYTest,bSolvedLinReg)
f = open(fileName*"standardLinReg.csv", "w")
write(f, "R² \n")
writecsv(f,standLinReg)
close(f)

#println(getRMSE(standXTest, standYTest, best3Beta[1,4:111]))
#println(getRMSE(standXTest, standYTest, best3Beta[2,4:111]))
#println(getRMSE(standXTest, standYTest, best3Beta[3,4:111]))
#Gurobi.writeproblem(stage2Model, "testproblem3.lp")
println("Complete")

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
export neardup

function neardup(X::AbstractVector{T}, epsilon=0.1) where T
    invindex = InvIndex()
    res = KnnResult(1)
    L = zeros(Int, length(X))
    D = zeros(Float64, length(X))
    L[1] = 1
    D[1] = 0.0
    push!(invindex, 1, X[1])
    for i in 2:length(X)
        empty!(res)
        x = X[i]
        search(invindex, cosine_distance, x, res)
        if length(res) == 0 || first(res).dist > epsilon
            push!(invindex, i, x)
            L[i] = i
            D[i] = 0.0
        else
            L[i] = first(res).objID
            D[i] = first(res).dist
        end
    end

    L, D
end


text:
include("todo_app.jl")
app = todo_app(false)
run_server(app)


text:
"""
    log_transformer(x)

Internal function used to get the exponent of x.
Used to select the nearest unit index.
"""
function log_transformer(x)
    return floor(Int, log10(x))
end

"""
    business_round(number::Number, names::Bool)

Internal function to round a number (num) to the nearest IG style.
"""
function business_round(number::Number, names::Bool)
    if number < 1000
        return string(number)
    else
        return extract_identifying_unit(number, names)
    end
end

"""
    number_groups(s::AbstractString)

Internal function to extract the two groups of numbers from a '.' separated string.
"""
function number_groups(s::AbstractString)
    group_1 = match(r"(\d+)", s)[1]
    group_2 = match(r"(\d+$)", s)[1]

    if length(group_2) > 1
        group_2 = group_2[1:2]
    else
        group_2 = group_2[1]
    end
    return group_1, group_2
end

"""
    extract_identifying_unit(num::Number, names::Bool)

Internal function to extract the identifying unit from a number.
"""
function extract_identifying_unit(num::Number, names::Bool)
    log_pos = log_transformer(num)
    div_3 = log_pos ÷ 3

    if names
        unit =  unit_names[div_3]
    else
        unit = units[div_3]
    end

    indentifier = string(num / ( BigInt(10) ^ (3 * div_3) ))
    main_int, rem_float = number_groups(indentifier)

    if names
        return string(main_int, ".", rem_float, " ", unit)
    else
        return string(main_int, ".", rem_float, unit)
    end
end

"""
    round(::Type{IGRound}, x::Number; names::Bool=false)

Function to round a number x to the nearest IG style.

# Arguments
 * x: The number to be rounded.
 * names: If true, the unit will be names, if false, the unit will be the number.

# Example

```julia-repl
julia> round(IGRound, 100_000; names=true)
"100Thousand"
```
"""
function round(::Type{IGRound}, x::Number; names::Bool=false)
    init_round =  Base.round(BigInt, x)
    return business_round(init_round, names)
end


text:
struct Case
	name::Symbol
	key::Symbol
	arg_names::Vector{Symbol}
	arg_typs::Vector{Any}
end
struct MatchType
	name::Symbol
	base::Any
	key::Symbol
	arguments::Symbol
	cstrs::Vector{Case}
end
function parse_constructor(expr)
	@assert(expr.head == :call)
	case_desc = expr.args[1]
	@assert(case_desc.head == :curly)
	case_name = case_desc.args[1]
	case_key = case_desc.args[2]
	arg_typs = Any[]
	arg_names = Symbol[]
	for (index, arg) in enumerate(expr.args[2:end])
		@assert(arg.head == :(::))
		if (length(arg.args) == 1)
			push!(arg_names, Symbol("_" * string(index)))
			push!(arg_typs, arg.args[1])
		elseif (length(arg.args) == 2)
			push!(arg_names, arg.args[1])
			push!(arg_typs, arg.args[2])
		end
	end
	return Case(case_name, case_key.value, arg_names, arg_typs)
end
function parse_constructors(expr)
	@assert(expr.head == :block)
	cstrs = Case[]
	for stmt in expr.args
		if !(stmt isa LineNumberNode || (stmt isa Expr && stmt.head == :call))
			throw("Invalid constructor declaration $stmt")
		end
		if stmt isa LineNumberNode continue end
		push!(cstrs, parse_constructor(stmt))
	end
	return cstrs
end
function parse_definition(expr)
	@assert(expr.head == :(=))
	name_defn = expr.args[1]
	@assert(name_defn.head == :curly)
	name = name_defn.args[1]
	basetype = eval(name_defn.args[2])
	keyname = name_defn.args[3]
	argsname = name_defn.args[4]
	print(expr.args[2])
	cstrs = parse_constructors(expr.args[2])
	return MatchType(name, basetype, keyname, argsname, cstrs)
end
function _matchtype(typedecl, cstrs)
	@assert(typedecl.head == :(<:))
	tyname = typedecl.args[1]
	host_decl = typedecl.args[2]
	@assert(host_decl isa Expr)
	@assert(host_decl.head == :curly)
	hostname = host_decl.args[1]
	tykey = host_decl.args[2]
	tyargs = host_decl.args[3]
	cases = parse_constructors(cstrs)

    expr = Expr(:toplevel)

    constructor_heads = getproperty.(cases, :name)
    case_keys = getproperty.(cases, :key)
    constructor_arguments = getproperty.(cases, :arg_typs)
    constructor_pairs = Pair.(constructor_heads, constructor_arguments)
    key_pairs = Pair.(constructor_heads, case_keys)

    push!(expr.args, quote
    	const $tyname = ASTMatching.UnionType{$hostname, $tykey, $tyargs}
	end)
    push!(expr.args, quote
    	headof(e::$tyname) = e.$tykey
	end)
    push!(expr.args, quote
    	argof(e::$tyname, i) = e.$tyargs[i]
	end)
    push!(expr.args, quote
    	constructors(::Type{$tyname}) = Dict($(constructor_pairs...))
	end)
	push!(expr.args, quote
		constructor_keys(::Type{$tyname}) = Dict($(key_pairs...))
	end)
    return esc(expr)
end
macro matchtype(typedecl, cstrs)
	return _matchtype(typedecl, cstrs)
end

#=
pat ::=
	_ => Str
	x => Var{:x}
	C(pat...) => Cstr{:key, Tuple{pat...}}
	pat || pat => Or{pat, pat}
=#

@generated function _matchfun(receiver::T) where T
end

macro astmatch(discriminant, patterns)
	_compile(patterns)
end


text:
# SAM Reader
# ==========

"""
    SAMReader(input::IO)

Create a data reader of the SAM file format.

# Arguments
* `input`: data source
"""
type SAMReader <: Bio.IO.AbstractReader
    state::Ragel.State
    header::SAMHeader

    function SAMReader(input::BufferedInputStream)
        reader = new(Ragel.State(samparser_start, input), SAMHeader())
        while !eof(input) && BufferedStreams.peek(input) == UInt8('@')
            # NOTE: This reads a header line, not the first record.
            @assert read(reader) == SAMRecord()
        end
        return reader
    end
end

function SAMReader(input::IO)
    return SAMReader(BufferedInputStream(input))
end

function Bio.IO.stream(reader::SAMReader)
    return reader.state.stream
end

function Base.show(io::IO, reader::SAMReader)
    println(io, summary(reader), ":")
      print(io, "  header keys: ", join(keys(reader.header), ", "))
end

function header(reader::SAMReader)
    return reader.header
end

function Base.eltype(::Type{SAMReader})
    return SAMRecord
end


text:
using StreamModels
using Test

@testset "StreamModels tests" begin
    include.(["formula.jl",
              "modelmatrix.jl",
              "summarizers.jl"])
end


text:
# This file was generated, do not modify it. # hide
for i ∈ 1:5, j ∈ 1:5
    print(" ", rpad("*"^i,5), lpad("*"^(6-i),5), j==5 ? "\n" : " "^4)
end


text:
# Posix separator and regex
const POSIX_PATH_SEPARATOR    = "/"
const POSIX_PATH_SEPARATOR_RE = r"/+"
const POSIX_PATH_ABSOLUTE_RE  = r"^/"
const POSIX_PATH_DIRECTORY_RE = r"(?:^|/)\.{0,2}$"
const POSIX_PATH_DIR_SPLITTER = r"^(.*?)(/+)([^/]*)$"
const POSIX_PATH_EXT_SPLITTER = r"^((?:.*/)?(?:\.|[^/\.])[^/]*?)(\.[^/\.]*|)$"

# Windows separator and regex
const WIN_PATH_SEPARATOR    = "\\"
const WIN_PATH_SEPARATOR_RE = r"[/\\]+"
const WIN_PATH_ABSOLUTE_RE  = r"^(?:\w+:)?[/\\]"
const WIN_PATH_DIRECTORY_RE = r"(?:^|[/\\])\.{0,2}$"
const WIN_PATH_DIR_SPLITTER = r"^(.*?)([/\\]+)([^/\\]*)$"
const WIN_PATH_EXT_SPLITTER = r"^((?:.*[/\\])?(?:\.|[^/\\\.])[^/\\]*?)(\.[^/\\\.]*|)$"

# The common READ, WRITE and EXEC constants
# to make working with file permissions easier.
const READ = 0o4
const WRITE = 0o2
const EXEC = 0o1

const USER_COEFF = 0o0100
const GROUP_COEFF = 0o0010
const OTHER_COEFF = 0o0001

# Stardard C style file mode constants taken from cpython's
# stat.py file
const S_IFDIR  = 0o040000     # directory
const S_IFCHR  = 0o020000     # character device
const S_IFBLK  = 0o060000     # block device
const S_IFREG  = 0o100000     # regular file
const S_IFIFO  = 0o010000     # fifo (named pipe)
const S_IFLNK  = 0o120000     # symbolic link
const S_IFSOCK = 0o140000     # socket file

const S_ISUID = 0o4000        # set UID bit
const S_ISGID = 0o2000        # set GID bit
const S_ENFMT = S_ISGID       # file locking enforcement
const S_ISVTX = 0o1000        # sticky bit
const S_IREAD = 0o0400        # Unix V7 synonym for S_IRUSR
const S_IWRITE = 0o0200       # Unix V7 synonym for S_IWUSR
const S_IEXEC = 0o0100        # Unix V7 synonym for S_IXUSR
const S_IRWXU = 0o0700        # mask for owner permissions
const S_IRUSR = 0o0400        # read by owner
const S_IWUSR = 0o0200        # write by owner
const S_IXUSR = 0o0100        # execute by owner
const S_IRWXG = 0o0070        # mask for group permissions
const S_IRGRP = 0o0040        # read by group
const S_IWGRP = 0o0020        # write by group
const S_IXGRP = 0o0010        # execute by group
const S_IRWXO = 0o0007        # mask for others (not in group) permissions
const S_IROTH = 0o0004        # read by others
const S_IWOTH = 0o0002        # write by others
const S_IXOTH = 0o0001        # execute by others

const FILEMODE_TABLE = (
    ((S_IFLNK,         'l'),
     (S_IFREG,         '-'),
     (S_IFBLK,         'b'),
     (S_IFDIR,         'd'),
     (S_IFCHR,         'c'),
     (S_IFIFO,         'p')),

    ((S_IRUSR,         'r'),),
    ((S_IWUSR,         'w'),),
    ((S_IXUSR|S_ISUID, 's'),
     (S_ISUID,         'S'),
     (S_IXUSR,         'x')),

    ((S_IRGRP,         'r'),),
    ((S_IWGRP,         'w'),),
    ((S_IXGRP|S_ISGID, 's'),
     (S_ISGID,         'S'),
     (S_IXGRP,         'x')),

    ((S_IROTH,         'r'),),
    ((S_IWOTH,         'w'),),
    ((S_IXOTH|S_ISVTX, 't'),
     (S_ISVTX,         'T'),
     (S_IXOTH,         'x'))
)

const DATA_SUFFIX = ["", "K", "M", "G", "T", "P", "E", "Z", "Y"]


text:
import Base: ∘

"""
    ∘(f::AbstractGP, g)

Constructs the GP f′ given by f′(x) := f(g(x))
"""
∘(f::AbstractGP, g) = GP(transform, f, g)

μ_p′(::typeof(transform), f::AbstractGP, g) = transform(mean(f), g)
k_p′(::typeof(transform), f::AbstractGP, g) = transform(kernel(f), g)
k_pp′(fp::GP, ::typeof(transform), f::AbstractGP, g) = transform(kernel(fp, f), g, Val(2))
k_p′p(::typeof(transform), f::AbstractGP, g, fp::GP) = transform(kernel(f, fp), g, Val(1))


text:
# Test utililties.jl

using Test, Statistics

println("testing utilites")

tol = 1e-2   # Fairly sloppy, needed for the derivative tests

# nonmaxsuppts - tested in test_cornerfeatures.jl

# derivative 3 5 7
# Generate an image of constant gradient
img = [c for r=1:50, c=1:50]
(gx,gy) = derivative3(img, ("x", "y"))

# Test interior of derivatives matrix matches expected values
@test all(abs.(gx[5:end-5,5:end-5] .- 1) .< tol)
@test all(abs.(gy[5:end-5,5:end-5]) .< tol)

(gx,gy, gxx, gyy, gxy) = derivative5(img, ("x", "y", "xx", "yy", "xy"))
# Test interior of derivatives matrix matches expected values
@test all(abs.(gx[5:end-5,5:end-5] .- 1) .< tol)
@test all(abs.(gy[5:end-5,5:end-5]) .< tol)
@test all(abs.(gxx[5:end-5,5:end-5]) .< tol)
@test all(abs.(gyy[5:end-5,5:end-5]) .< tol)
@test all(abs.(gxy[5:end-5,5:end-5]) .< tol)

(gx,gy, gxx, gyy, gxy) = derivative7(img, ("x", "y", "xx", "yy", "xy"))
# Test interior of derivatives matrix matches expected values
@test all(abs.(gx[5:end-5,5:end-5] .- 1) .< tol)
@test all(abs.(gy[5:end-5,5:end-5]) .< tol)
@test all(abs.(gxx[5:end-5,5:end-5]) .< tol)
@test all(abs.(gyy[5:end-5,5:end-5]) .< tol)
@test all(abs.(gxy[5:end-5,5:end-5]) .< tol)

# gaussfilt
img = zeros(100,100)
img[50,50] = 1
simg = gaussfilt(img,2)
@test abs(sum(simg) - 1) < tol

# dilate1d erode1d
# Generate a 1D signal with a single spike
v = zeros(20)
v[10] = 1
k = 5
vd = dilate1d(v,k)
@test all(vd[8:12] .== 1) && all(vd[1:7] .== 0) && all(vd[13:end] .== 0)

ve = erode1d(vd,k)   # Erode back and check we reconstruct v
@test all(ve .== v)

# imdilate imerode
img = falses(100,100)
img[50,50] = true
dimg = imdilate(img, "rect", [5, 3])
@test all(dimg[47:47,48:52] .== [0 0 0 0 0])
@test all(dimg[50:50,48:52] .== [0 1 1 1 0])
@test all(dimg[53:53,48:52] .== [0 0 0 0 0])
@test all(dimg[47:53,50] .== [0, 1, 1, 1, 1, 1, 0])

eimg = imerode(dimg, "rect", [5, 3])
@test all(eimg[49:51,49:51] .== [0 0 0; 0 1 0; 0 0 0])

# circularstruct

# floatxy

# matchbycorrelation

# grey2census

# imgnormalise
a = rand(5,3)
an = imgnormalise(a)
@test isapprox(minimum(an), 0)
@test isapprox(maximum(an), 1)

reqmean = 5
reqvar = 7
an = imgnormalise(a, reqmean, reqvar)
@test isapprox(mean(an), reqmean)
@test isapprox(var(an), reqvar)

# histtruncate
rows = 100;
cols = 200;
img = rand(rows, cols)
lHistCut = 2
uHistCut = 4
htimg = histtruncate(img, lHistCut, uHistCut)
minval = minimum(htimg)
maxval = maximum(htimg)
# Check that the number of saturated pixels at each extreme is
# approximately lHistCut and uHistCut. Note some integer rounding has
# to occur.
v = findall(abs.(htimg .- minval).<eps())
@test abs(length(v)/(rows*cols) * 100 - lHistCut) < 1

v = findall(abs.(htimg .- maxval).<eps())
@test abs(length(v)/(rows*cols) * 100 - uHistCut) < 1


text:
"""
    Checkpoints

A very minimal module for defining checkpoints or save location in large codebase with
the ability to configure how those checkpoints save data externally
(similar to how Memento.jl works for logging).
"""
module Checkpoints

using AWSS3
using Compat # for contains (julia v1.5)
using ContextVariablesX
using DataStructures: DefaultDict
using FilePathsBase
using FilePathsBase: /, join
using JLSO
using Memento
using OrderedCollections

export checkpoint, with_checkpoint_tags  # creating stuff
export enabled_checkpoints
# indexing stuff
export IndexEntry, index_checkpoint_files, index_files
export checkpoint_name, checkpoint_path, prefixes, tags

const LOGGER = getlogger(@__MODULE__)

__init__() = Memento.register(LOGGER)

include("handler.jl")

const CHECKPOINTS = Dict{String, Union{Nothing, Handler}}()
@contextvar CONTEXT_TAGS::Tuple{Vararg{Pair{Symbol, Any}}} = Tuple{}()

include("session.jl")
include("indexing.jl")
include("deprecated.jl")

"""
    with_checkpoint_tags(f::Function, context_tags::Pair...)
    with_checkpoint_tags(f::Function, context_tags::NamedTuple)

Runs the function `f`, tagging any [`checkpoint`](@ref)s created by `f` with the `context_tags`.
This is normally used via the do-block form:
For example

```julia
with_checkpoint_tags(:foo=>1, :bar=>2) do
    q_out = qux()
    checkpoint("foobar"; :output=q_out)
end
```
This snippet will result in `"foobar"` checkpoint having the `foo=1` and `bar=2` tags, as will any checkpoints created by `qux`().
The context tags are [dynamically scoped](https://en.wikipedia.org/wiki/Scope_(computer_science)#Lexical_scope_vs._dynamic_scope_2) and so are retained through function calls.

Nested contexts (nested `with_checkpoint_tags` calls) are allowed. Duplicate tag names and values are
allowed, including the tags provided directly in the [`checkpoint`](@ref) call.
Duplicate tags are repeated, not overwritten.
"""
function with_checkpoint_tags(f::Function, context_tags::Pair...)
    with_context(f, CONTEXT_TAGS => (CONTEXT_TAGS[]..., context_tags...))
end
with_checkpoint_tags(f::Function, context_tags::NamedTuple) = with_checkpoint_tags(f, pairs(context_tags)...)

"""
    available() -> Vector{String}

Returns a vector of all available (registered) checkpoints.
"""
available() = collect(keys(CHECKPOINTS))

"""
    enabled_checkpoints() -> Vector{String}

Returns a vector of all enabled ([`config`](@ref)ured) checkpoints.
"""
enabled_checkpoints() = filter(k -> CHECKPOINTS[k] !== nothing, available())

"""
    checkpoint([prefix], name, data)
    checkpoint([prefix], name, data::Pair...)
    checkpoint([prefix], name, data::Dict)

Defines a data checkpoint with a specified `label` and values `data`.
By default checkpoints are no-ops and need to be explicitly configured.

    checkpoint(session, data)
    checkpoint(handler, name, data::Dict)

Alternatively, you can also checkpoint with to a session which stages the data to be
commited later by `commit!(session)`.
Explicitly calling checkpoint on a handler is generally not advised, but is an option.
"""
function checkpoint(name::String, data::Dict{Symbol}; tags...)
    checkpoint_deprecation(tags...)
    with_checkpoint_tags(tags...) do
        checkpoint(CHECKPOINTS[name], name, data)
    end
end

function checkpoint(name::String, data::Pair...; tags...)
    checkpoint_deprecation(tags...)
    with_checkpoint_tags(tags...) do
        checkpoint(name, Dict(data...))
    end
end

function checkpoint(name::String, data; tags...)
    checkpoint_deprecation(tags...)
    with_checkpoint_tags(tags...) do
        checkpoint(name, Dict(:data => data))
    end
end

function checkpoint(prefix::Union{Module, String}, name::String, args...; tags...)
    checkpoint_deprecation(tags...)
    with_checkpoint_tags(tags...) do
        checkpoint("$prefix.$name", args...)
    end
end

"""
    config(handler::Handler, labels::Vector{String})
    config(handler::Handler, prefix::String)
    config(labels::Vector{String}, args...; kwargs...)
    config(prefix::String, args...; kwargs...)

Configures the specified checkpoints with a `Handler`.
If the first argument is not a `Handler` then all `args` and `kwargs` are passed to a
`Handler` constructor for you.
"""
function config(handler::Handler, names::Vector{String})
    for n in names
        haskey(CHECKPOINTS, n) || warn(LOGGER, "$n is not a registered checkpoint")
        debug(LOGGER, "Checkpoint $n set to use $(handler)")
        CHECKPOINTS[n] = handler
    end
end

function config(handler::Handler, prefix::Union{Module, String})
    config(handler, filter(l -> startswith(l, prefix), available()))
end

function config(names::Vector{String}, args...; kwargs...)
    config(Handler(args...; kwargs...), names)
end

function config(prefix::Union{Module, String}, args...; kwargs...)
    config(Handler(args...; kwargs...), prefix)
end

"""
    register([prefix], labels)

Registers a checkpoint that may be configured at a later time.
"""
function register end

function register(labels::Vector{String})
    for l in labels
        if haskey(CHECKPOINTS, l)
            warn(LOGGER, "$l has already registered")
        else
            CHECKPOINTS[l] = nothing
        end
    end
end

function register(prefix::Union{Module, String}, labels::Vector{String})
    register(map(l -> join([prefix, l], "."), labels))
end

end  # module


text:
# Logistic regression

using Regression
import Regression: solve, Options

# Data

d = 3
n = 1000

w = randn(d+1)
X = randn(d, n)
y = sign(X'w[1:d] + w[d+1] + 0.2 * randn(n))

# Solve

ret = solve(logisticreg(X, y; bias=1.0),
            reg=SqrL2Reg(1.0e-2),
            options=Options(verbosity=:iter, grtol=1.0e-6 * n))

println()

# Print results

w_e = ret.sol

@printf("corr(truth, estimated) = %.6f\n", dot(w, w_e) / (vecnorm(w) * vecnorm(w_e)))


text:
# Description:
# Dimensions: d
#
# The Michalewicz function has d! local minima, and it is multimodal. The parameter m defines the steepness of they valleys and ridges; a larger m leads to a more difficult search. The recommended value of m is m = 10. The function's two-dimensional form is shown in the plot above.
#
# Input Domain:
# The function is usually evaluated on the hypercube xi ∈ [0, π], for all i = 1, …, d.
#
# Global Minima:
# https://www.sfu.ca/~ssurjano/michal3.png

function michal(xx, m = 10)

    d = length(xx);
    sum = 0;

    for ii = 1:d
       	xi = xx[ii];
       	new = sin(xi) * (sin(ii * xi^2 / pi))^(2 * m);
       	sum  = sum + new;
    end

    y = -sum;

    return [y]
end


text:
using SafeTestsets

module ChargedParticle3dTests

    using Test
    using GeometricIntegrators

    const Δt = 0.01
    const nt = 1

    export test_charged_particle_3d

    function test_charged_particle_3d(ode::ODE)
        int = Integrator(ode, getTableauGLRK(1), Δt)
        sol = integrate(ode, int, nt)
        @test true
    end

    function test_charged_particle_3d(ode::PODE)
        int = Integrator(ode, getTableauIPGLRK(1), Δt)
        sol = integrate(ode, int, nt)
        @test true
    end

    function test_charged_particle_3d(ode::IODE)
        int = Integrator(ode, getTableauVPGLRK(1), Δt)
        sol = integrate(ode, int, nt)
        @test true
    end

end

@safetestset "Charged Particle Dynamics in 3D in Symmetric Magnetic Field                                         " begin

    using ChargedParticleDynamics.ChargedParticle3d.SymmetricField
    using ..ChargedParticle3dTests

    test_charged_particle_3d(SymmetricField.charged_particle_3d_iode())

end

@safetestset "Charged Particle Dynamics in 3D in Singular Magnetic Field                                          " begin

    using ChargedParticleDynamics.ChargedParticle3d.SingularField
    using ..ChargedParticle3dTests

    test_charged_particle_3d(SingularField.charged_particle_3d_iode())

end

@safetestset "Charged Particle Dynamics in 3D in Theta Pinch                                                      " begin

    using ChargedParticleDynamics.ChargedParticle3d.ThetaPinchCanonical
    using ..ChargedParticle3dTests

    test_charged_particle_3d(ThetaPinchCanonical.charged_particle_3d_iode())

end

@safetestset "Charged Particle Dynamics in 3D in Theta Pinch (noncanonical formulation)                           " begin

    using ChargedParticleDynamics.ChargedParticle3d.ThetaPinchNoncanonical
    using ..ChargedParticle3dTests

    test_charged_particle_3d(ThetaPinchNoncanonical.charged_particle_3d_iode())

end

@safetestset "Charged Particle Dynamics in 3D in Tokamak (cartesian coordinates)                                  " begin

    using ChargedParticleDynamics.ChargedParticle3d.TokamakSmallCartesian
    using ..ChargedParticle3dTests

    test_charged_particle_3d(TokamakSmallCartesian.charged_particle_3d_pode())
    test_charged_particle_3d(TokamakSmallCartesian.charged_particle_3d_iode())

end

@safetestset "Charged Particle Dynamics in 3D in Tokamak (cylindrical coordinates)                                " begin

    using ChargedParticleDynamics.ChargedParticle3d.TokamakSmallCylindrical
    using ..ChargedParticle3dTests

    test_charged_particle_3d(TokamakSmallCylindrical.charged_particle_3d_pode())
    test_charged_particle_3d(TokamakSmallCylindrical.charged_particle_3d_iode())

end

@safetestset "Charged Particle Dynamics in 3D in Tokamak (toroidal coordinates)                                   " begin

    using ChargedParticleDynamics.ChargedParticle3d.TokamakSmallToroidal
    using ..ChargedParticle3dTests

    test_charged_particle_3d(TokamakSmallToroidal.charged_particle_3d_pode())
    test_charged_particle_3d(TokamakSmallToroidal.charged_particle_3d_iode())

end

@safetestset "Charged Particle Dynamics in 3D in Tokamak (noncanonical formulation)                               " begin

    using ChargedParticleDynamics.ChargedParticle3d.TokamakSmallNoncanonical
    using ..ChargedParticle3dTests

    test_charged_particle_3d(TokamakSmallNoncanonical.charged_particle_3d_iode())

end


text:
# GFF3 File Format
# ================

module GFF3

using BioGenerics
using Indexes
using FASTX.FASTA #TODO: move responsibility to FASTX.jl.
using TranscodingStreams

import Automa
import Automa.RegExp: @re_str
import Automa.Stream: @mark, @markpos, @relpos, @abspos

import BGZFStreams
import BioGenerics.Exceptions: missingerror
import BioSequences
import GenomicFeatures: GenomicFeatures, Interval, IntervalCollection
import URIParser

include("record.jl")
include("reader.jl")
include("writer.jl")

end # module


text:
using SafeTestsets

@safetestset "SoftmaxTests" begin include("softmaxtests.jl") end
@safetestset "NormalTests" begin include("normaltests.jl") end


text:
# # Solve the heat conduction equation

# Synopsis: Compute the solution of the Poisson equation of heat conduction with a
# nonzero heat source. Quadrilateral four-node elements are used.

# The problem is linear heat conduction equation posed on a bi-unit square,
# solved with Dirichlet boundary conditions around the circumference. Uniform
# nonzero heat generation rate is present. The exact solution is in this way
# manufactured and hence known. That gives us an opportunity to calculate the
# true error.

# The complete code is in the file [`tut_poisson_q4.jl`](tut_poisson_q4.jl).

# The solution will be defined  within a module in order to eliminate conflicts
# with data or functions defined elsewhere.

module tut_poisson_q4

# We'll need some functionality from linear algebra, and the mesh libraries.
# Finally we will need the `Elfel` functionality.
using LinearAlgebra
using MeshCore.Exports
using MeshSteward.Exports
using Elfel.Exports

# This is the top level function.
function run()
    # Input parameters:
    A = 1.0 # length of the side of the square
    kappa =  1.0; # thermal conductivity of the material
    Q = -6.0; # internal heat generation rate
    tempf(x, y) =(1.0 + x^2 + 2.0 * y^2); # the exact distribution of temperature
    N = 1000; # number of element edges along the sides of the square domain

    # Generate the computational mesh.
    mesh = genmesh(A, N)

    # Create the finite element space to represent the temperature solution. The
    # degrees of freedom are real numbers (`Float64`), the quadrilaterals are
    # defined by the mesh, and each of the elements has the continuity ``H
    # ^1``, i. e. both the function values and the derivatives are square
    # integrable.
    Uh = FESpace(Float64, mesh, FEH1_Q4())

    # Apply the essential boundary conditions at the circumference of the square
    # domain. We find the boundary incidence relation (`boundary(mesh)`), and
    # then the list of all vertices connected by the boundary cells. The
    # function `tempf` defines the analytical temperature variation, and hence
    # for each of the vertices `i` on the boundary (they are of manifold
    # dimension  `0`), we set the component of the field (1) to the exact value
    # of the temperature at that location.
    vl = connectedv(boundary(mesh));
    locs = geometry(mesh)
    for i in vl
        setebc!(Uh, 0, i, 1, tempf(locs[i]...))
    end

    # Number the degrees of freedom, both the unknowns and the data
    # (prescribed) degrees of freedom.
    numberdofs!(Uh)
    @show ndofs(Uh), nunknowns(Uh)

    # Assemble the conductivity matrix and the vector of the heat loads. Refer
    # to the definition of this function below.
    K, F = assembleKF(Uh, kappa, Q)

    # This is a vector to hold all degrees of freedom in the system.
    T = fill(0.0, ndofs(Uh))
    # Here we collect the data degrees of freedom (the known values).
    gathersysvec!(T, Uh)

    # The system of linear algebraic equations is solved.
    solve!(T, K, F, nunknowns(Uh))

    # The values of all the degrees of freedom can now be introduced into the
    # finite element space.
    scattersysvec!(Uh, T)

    # Here we associate the values of the finite element space with the entities
    # of the mesh as an attribute.
    makeattribute(Uh, "T", 1)

    # The correctness of the solution is checked by comparing the values at the
    # vertices.
    checkcorrectness(Uh, tempf)

    # The attribute can now be written out for visualization into a VTK file.
    vtkwrite("q4-T", baseincrel(mesh), [(name = "T",)])

    true # return success
end

# The domain is a square, meshed with quadrilateral elements. The function
# `Q4block` creates an incidence relation that defines the quadrilateral
# element shapes by the vertices connected into the shapes. This incidence
# relation is then attached to the mesh and the mesh is returned.
function genmesh(A, N)
    conn = Q4block(A, A, N, N)
    return attach!(Mesh(), conn)
end

# The `assembleKF` function constructs the left-hand side coefficient matrix,
# conductivity matrix, as a sparse matrix, and a vector of the heat loads due
# to the internal heat generation rate `Q`.

# The boundary value problem is expressed in this weak form
# ```math
# \int_{V}(\mathrm{grad}\vartheta)\; \kappa (\mathrm{grad}T
#             )^T\; \mathrm{d} V
#             -\int_{V}  \vartheta Q \; \mathrm{d} V
#              = 0
# ```
# where the test function vanishes on the boundary where the temperature is
# prescribed, ``\vartheta(x) =0``  for  ``x \in{S_1}``
# Substituting ``\vartheta = N_j `` and ``T = \sum_i N_i T_i`` we obtain the
# linear algebraic equations
# ```math
# \sum_i T_i \int_{V} \mathrm{grad}N_j \; \kappa (\mathrm{grad}N_i)^T\; \mathrm{d} V
#  -\int_{V}  N_j Q \; \mathrm{d} V  = 0 , \quad \forall j.
# ```
# The volume element is ``\mathrm{d} V``, which in our case
# becomes ``1.0\times\mathrm{d} S``, since the thickness of the two
# dimensional domain is assumed to be 1.0.

function assembleKF(Uh, kappa, Q)
    # At the top of the `assembleKF` we look at the function `integrate!` to
    # evaluate the weak-form integrals. The key to making this calculation
    # efficient is type stability. All the arguments coming in must have
    # concrete types. This is why the `integrate!` function is an inner
    # function: the function barrier allows for all arguments to be resolved to
    # concrete types.
    function integrate!(am, av, elit, qpit, kappa, Q)
        nedof = ndofsperel(elit)
        # The local assemblers are just like matrices or vectors
        ke = LocalMatrixAssembler(nedof, nedof, 0.0)
        fe = LocalVectorAssembler(nedof, 0.0)
        for el in elit # Loop over all elements
            init!(ke, eldofs(el), eldofs(el)) # zero out elementwise matrix
            init!(fe, eldofs(el)) # and vector
            for qp in qpit # Now loop over the quadrature points
                Jac, J = jacjac(el, qp) # Calculate the Jacobian matrix, Jacobian
                gradN = bfungrad(qp, Jac) # Evaluate the spatial gradients
                JxW = J * weight(qp) # elementary volume
                N = bfun(qp) # Basis function values at the quadrature point
                # This double loop evaluates the elementwise conductivity
                # matrix and the heat load vector precisely as the formula of
                # the weak form  dictates; see above.
                for i in 1:nedof
                    for j in 1:nedof
                        ke[j, i] += dot(gradN[j], gradN[i]) * (kappa * JxW)
                    end
                    fe[j] += N[j] * Q * JxW
                end
            end
            # Assemble the calculated contributions from this element
            assemble!(am, ke)
            assemble!(av, fe)
        end
        return am, av # Return the updated assemblers
    end

    # In the `assembleKF` function we first we create the element iterator. We
    # can go through all the elements that define the domain of integration
    # using this iterator. Each time a new element is accessed, some data are
    # precomputed such as the element degrees of freedom.
    elit = FEIterator(Uh)
    # This is the quadrature point iterator. We know that the elements are
    # quadrilateral, which makes the Gauss integration rule the obvious choice.
    # We also select order 2 for accuracy. Quadrature-point iterators provide
    # access to basis function values and gradients, the Jacobian matrix and
    # the Jacobian determinant, the location of the quadrature point and so
    # on.
    qpit = QPIterator(Uh, (kind = :Gauss, order = 2))
    # Next we create assemblers, one for the sparse system matrix and one for
    # the system vector.
    am = start!(SysmatAssemblerSparse(0.0), ndofs(Uh), ndofs(Uh))
    av = start!(SysvecAssembler(0.0), ndofs(Uh))
    # Now we call the integration function. The assemblers are modified inside
    # this function...
    @time integrate!(am, av, elit, qpit, kappa, Q)
    # ...so that when the integration is done, we can materialize the sparse
    # matrix and the vector and return them.
    return finish!(am), finish!(av)
end

# The linear algebraic system is solved by partitioning. The vector `T` is
# initially all zero, except in the degrees of freedom which are prescribed as
# nonzero. Therefore the product of the conductivity matrix and the vector `T`
# are the heat loads due to nonzero essential boundary conditions. To this we
# add the vector of heat loads due to the internal heat generation rate. The
# submatrix of the heat conduction matrix corresponding to the free degrees of
# freedom (unknowns), `K[1:nu, 1:nu]` is then used to solve for the unknowns `T
# [1:nu]`.
function solve!(T, K, F, nu)
    @time KT = K * T
    @time T[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

# The correctness can be checked in various ways. Here we calculate the mean
# deviation of the calculated temperatures at the nodes relative to the exact
# values of the temperature.
function checkcorrectness(Uh, tempf)
    geom = geometry(Uh.mesh)
    ir = baseincrel(Uh.mesh)
    T = attribute(ir.right, "T")
    std = 0.0
    for i in 1:length(T)
        std += abs(T[i][1] - tempf(geom[i]...))
    end
    @show (std / length(T)) <= 1.0e-9
end

end # module

# The module can now be used.
using .tut_poisson_q4
tut_poisson_q4.run()


text:
const city = Dict(
    "Malmö FF"=>"Malmö",
    "IFK Göteborg"=>"Göteborg",
    "AIK"=>"Stockholm",
    "IFK Norrköping"=>"Norrköping",
    "IF Elfsborg"=>"Borås",
    "Helsingborgs IF"=>"Helsingborg",
    "Djurgårdens IF"=>"Stockholm",
    "Örgryte IS"=>"Göteborg",
    "Halmstads BK"=>"Halmstad",
    "GAIS"=>"Göteborg",
    "Örebro SK"=>"Örebro",
    "Hammarby IF"=>"Stockholm",
    "Östers IF"=>"Växjö",
    "Kalmar FF"=>"Kalmar",
    "Landskrona BoIS"=>"Landskrona",
    "Degerfors IF"=>"Degerfors",
    "Åtvidabergs FF"=>"Åtvidaberg",
    "Sandvikens IF"=>"Sandviken",
    "BK Häcken"=>"Göteborg",
    "Trelleborgs FF"=>"Trelleborg",
    "IK Brage"=>"Borlänge",
    "IK Sleipner"=>"Norrköping",
    "Gefle IF"=>"Gävle",
    "GIF Sundsvall"=>"Sundsvall",
    "IFK Malmö"=>"Malmö",
    "IFK Eskilstuna"=>"Eskilstuna",
    "Jönköpings Södra IF"=>"Jönköping",
    "Västra Frölunda IF"=>"Göteborg",
    "Mjällby AIF"=>"Hällevik",
    "IS Halmia"=>"Halmstad",
    "Gårda BK"=>"Göteborg",
    "IFK Sundsvall"=>"Sundsvall",
    "IF Brommapojkarna"=>"Stockholm",
    "Västerås SK"=>"Västerås",
    "Syrianska FC"=>"Södertälje",
    "Falkenbergs FF"=>"Falkenberg",
    "IK Sirius"=>"Uppsala",
    "Råå IF"=>"Helsingborg",
    "Ljungskile SK"=>"Ljungskile",
    "Östersunds FK"=>"Östersund",
    "Westermalms IF"=>"Stockholm",
    "Umeå FC"=>"Umeå",
    "IFK Uddevalla"=>"Uddevalla",
    "Hallstahammars SK"=>"Hallstahammar",
    "Stattena IF"=>"Helsingborg",
    "Motala AIF"=>"Motala",
    "Redbergslids IK"=>"Göteborg",
    "Ludvika FfI"=>"Ludvika",
    "IK Oddevold"=>"Uddevalla",
    "IFK Luleå"=>"Luleå",
    "IF Saab"=>"Linköping",
    "Reymersholms IK"=>"Stockholm",
    "Norrby IF"=>"Borås",
    "BK Derby"=>"Linköping",
    "Assyriska FF"=>"Södertälje",
    "Brynäs IF"=>"Gävle",
    "Enköpings SK"=>"Enköping",
    "Högadals IS"=>"Karlshamn",
    "Västerås IK"=>"Västerås",
    "IFK Holmsund"=>"Holmsund",
    "Sandvikens AIK"=>"Sandviken",
    "IK City"=>"Eskilstuna",
    "Billingsfors IK"=>"Billingsfors")


text:
function decomp_genocol!{T <: Float}(
    Y     :: DenseVecOrMat{T},
    x     :: BEDFile{T},
    col   :: Int,
    snp   :: Int,
    m     :: T,
    d     :: T;
    quiet :: Bool = true
)
    @inbounds for case = 1:x.geno.n
        t = x.geno[case,snp]
        Y[case,col] = t == ONE8 ? zero(T) : (int2geno(x,t) - m)*d
        quiet || println("Y[$case,$col] = ", Y[case,col])
    end
    return nothing
end

function decomp_genocol!{T <: Float}(
    Y      :: DenseVecOrMat{T},
    x      :: BEDFile{T},
    mask_n :: DenseVector{Int},
    col    :: Int,
    snp    :: Int,
    m      :: T,
    d      :: T;
    quiet  :: Bool = true
)
    @inbounds for case = 1:x.geno.n
        if mask_n[case] == 1
            t = x[case,snp]
            Y[case,col] = t == ONE8 ? zero(T) : (int2geno(x,t) - m)*d
            quiet || println("Y[$case,$col] = ", Y[case,col])
        else
            Y[case,col] = zero(T)
        end
    end
    return nothing
end

function decomp_covarcol!{T <: Float}(
    Y     :: DenseVecOrMat{T},
    x     :: BEDFile{T},
    col   :: Int,
    covar :: Int,
    m     :: T,
    d     :: T;
    quiet :: Bool = true
)
    @inbounds for case = 1:x.geno.n
        Y[case,col] = (x.covar.x[case,(covar-x.geno.p)] - m) * d
        quiet || println("Y[$case,$col] = ", Y[case,col])
    end
end

function decomp_covarcol!{T <: Float}(
    Y      :: DenseVecOrMat{T},
    x      :: BEDFile{T},
    mask_n :: DenseVector{Int},
    col    :: Int,
    covar  :: Int,
    m      :: T,
    d      :: T;
    quiet  :: Bool = true
)
    @inbounds for case = 1:x.geno.n
        if mask_n[case] == 1
            t = x[case,snp]
            Y[case,col] = (x.covar.x[case,(covar-x.geno.p)] - m) * d
            quiet || println("Y[$case,$col] = ", Y[case,col])
        else
            Y[case,col] = zero(T)
        end
    end
end

"""
    decompress_genotypes!(Y, x, [, standardize=true])

Can also be called with a matrix `Y`, in which case all genotypes are decompressed.
Use this function judiciously, since the memory demands from decompressing large portions of `x` can grow quite large.
Use optional argument `standardize` to control standardization of allele dosages.
"""
function decompress_genotypes!{T <: Float}(
    Y       :: DenseMatrix{T},
    x       :: BEDFile{T};
    pids    :: DenseVector{Int} = procs(),
    quiet   :: Bool             = true
)

    # get dimensions of matrix to fill
    (n,p) = size(Y)
    xn,xp = size(x)

    # ensure dimension compatibility
    @assert n == xn "column of Y is not of same length as column of uncompressed x"
    @assert p <= xp "Y has more columns than x"

    ## TODO: parallelize this for SharedArrays
    @inbounds for j = 1:xp
        m = x.means[j]
        s = x.precs[j]
        if j <= x.geno.p
            decomp_genocol!(Y, x, j, j, m, s, quiet=quiet)
        else
            decomp_covarcol!(Y, x, j, j, m, s, quiet=quiet)
        end
    end
    return nothing
end

"""
    decompress_genotypes!(Y, x, idx, means, invstds)

When `Y` is a matrix and `idx` is a `BitArray` or `Int` array that indexes the columns of `x`, then `decompress_genotypes!()` decompresses only a subset of `x`.
"""
function decompress_genotypes!{T <: Float}(
    Y     :: DenseMatrix{T},
    x     :: BEDFile{T},
    idx   :: BitArray{1};
    pids  :: DenseVector{Int} = procs(),
    quiet :: Bool             = true
)

    # get dimensions of matrix to fill
    (n,p)  = size(Y)
    xn, xp = size(x)

    # ensure dimension compatibility
    @assert n == xn "column of Y is not of same length as column of uncompressed x"
    @assert p <= xp "Y has more columns than x"
    @assert sum(idx) <= xp "Vector 'idx' indexes more columns than are available in Y"

    # counter to ensure that we do not attempt to overfill Y
    col = 0

    ## TODO: parallelize this for SharedArrays
    @inbounds for snp = 1:xp

        # use this column?
        if idx[snp]

            # add to counter
            col += 1
            quiet || println("filling current column $current_col with snp $snp")

            # extract column mean, precision
            m = x.means[snp]
            s = x.precs[snp]

            if snp <= x.geno.p
                decomp_genocol!(Y, x, col, snp, m, s, quiet=quiet)
            else
                decomp_covarcol!(Y, x, col, snp, m, s, quiet=quiet)
            end

            # quit when Y is filled
            col == p && return nothing
        end
    end
    return nothing
end

function decompress_genotypes!{T <: Float}(
    Y     :: DenseMatrix{T},
    x     :: BEDFile{T},
    idx   :: DenseVector{Int};
    pids  :: DenseVector{Int} = procs(),
    quiet :: Bool             = true
)

    # get dimensions of matrix to fill
    (n,p) = size(Y)
    xn,xp = size(x)

    # ensure dimension compatibility
    @assert n == xn "column of Y is not of same length as column of uncompressed x"
    @assert p <= xp "Y has more columns than x"
    @assert length(idx) <= p "Vector 'idx' indexes more columns than are available in Y"

    # counter to ensure that we do not attempt to overfill Y
    col = 0

    @inbounds for snp in idx

        # add to counter
        col += 1
        quiet || println("filling current column $current_col with snp $snp")

        # extract column mean, precision
        m = x.means[snp]
        s = x.precs[snp]
        if snp <= x.geno.p
            decomp_genocol!(Y, x, col, snp, m, s, quiet=quiet)
        else
            decomp_covarcol!(Y, x, col, snp, m, s, quiet=quiet)
        end

        # quit when Y is filled
        col == p && return nothing
    end

    return nothing
end

"""
    decompress_genotypes!(Y, x, idx, mask_n [,pids=procs(), quiet=true])

If called with a vector `mask_n` of `0`s and `1`s, then `decompress_genotypes!()` will decompress the columns of `x` indexed by `idx` into `Y` while masking the rows indexed by `mask_n`.
"""
function decompress_genotypes!{T <: Float}(
    Y      :: DenseMatrix{T},
    x      :: BEDFile{T},
    idx    :: BitArray{1},
    mask_n :: DenseVector{Int};
    pids   :: DenseVector{Int} = procs(),
    quiet  :: Bool             = true
)

    # get dimensions of matrix to fill
    (n,p) = size(Y)
    xn,xp = size(x,2)

    # ensure dimension compatibility
    @assert n <= xn "column dimension of of Y exceeds column dimension of uncompressed x"
    @assert n == length(mask_n) "bitmask mask_n indexes different number of columns than column dimension of Y"
    @assert p <= xp "Y has more columns than x"
    @assert sum(idx) <= xp "Vector 'idx' indexes more columns than are available in Y"

    # counter to ensure that we do not attempt to overfill Y
    col = 0

    @inbounds for snp = 1:xp

        # use this column?
        if idx[snp]

            # add to counter
            col += 1
            quiet || println("filling current column $col with snp $snp")

            # extract column mean, precision
            m = x.means[snp]
            s = x.precs[snp]
            if snp <= x.geno.p
                decomp_genocol!(Y, x, mask_n, col, snp, m, s, quiet=quiet)
            else
                decomp_covarcol!(Y, x, mask_n, col, snp, m, s, quiet=quiet)
            end

            # quit when Y is filled
            col == p && return nothing
        end
    end
    return nothing
end

"""
    decompress_genotypes!(y, x, col [, quiet=true])

This function decompresses into vector `y` a standardized column (SNP) of a PLINK BED file housed in `x`.
Missing genotypes are sent to zero.

Arguments:

- `y` is the vector to fill with (standardized) dosages.
- `x` is the BEDfile object that contains the compressed `n` x `p` design matrix.
- `col` is the current predictor to extract.

Optional Arguments:

- `quiet` is a `Bool` to activate output. Defaults to `true` (no output). Used mostly for debugging.
"""
function decompress_genotypes!{T <: Float}(
    y     :: DenseVector{T},
    x     :: BEDFile{T},
    col   :: Int;
    quiet :: Bool = true
)
    m = x.means[col]
    s = x.precs[col]
    if col <= x.geno.p
        decomp_genocol!(y, x, 1, col, m, s, quiet=quiet)
    else
        decomp_covarcol!(y, x, 1, col, m, s, quiet=quiet)
    end
    return nothing
end

max_stars_repo_name: klkeys/PLINK

max_issues_repo_name: klkeys/PLINK


text:
@testset "98.validate-binary-search-tree.jl" begin
    @test is_valid_BST(TreeNode{Int}([5, 1, 4, nothing, nothing, 3, 6])) == false
    @test is_valid_BST(TreeNode{Int}([5, 1, 6])) == true
end


text:
module FishABM

using DataFrames, Distributions, Gadfly, ProgressMeter

export
  # Types
  AgentAssumptions,
  EnvironmentAssumptions,
  StockAssumptions,
  StockDB,

  # Agent functions
  agentdb,
  kill!,
  localmove,
  move!,
  injectagents!,

  # Stock functions
  harvest!,
  ageadults!,

  # Agent-Stock function
  spawn!,
  graduate!,

  # Utilities
  pad_environment!,
  plot,

  # Simulate
  simulate

include("types.jl")
include("utilities.jl")
include("agents.jl")
include("stock.jl")
include("agent_stock_interaction.jl")
include("simulate.jl")

end # module


text:
"""
    NexHotspot(nexSerial, name; pid=pid, cid=cid)

A Nextion Hotspot UI control.
"""
struct NexHotspot <: AbstractNexObject
    _nid::NexID

    touchable::ITouchable

    function NexHotspot(nexSerial::T, name::Name; pid=PageID(), cid=ComponentID()) where {T <: AbstractNexSerial}
        nid = NexID(nexSerial, name, pid, cid)
        new(nid, ITouchable(nid))
    end
end

function setproperty!(obj::NexHotspot, property::Symbol, new_val)
    _setcommonproperty!(obj, property, new_val)
end

function getproperty(obj::NexHotspot, property::Symbol)
    _getcommonproperty(obj, property)
end


text:
using IncrementalInference

global fg = initfg()

global N=100

global doors = [-100.0;0.0;100.0;300.0]'
global cov = [3.0]

global v1 = addVariable!(fg,:x1,ContinuousScalar,N=N)
global f1  = addFactor!(fg,[v1], Prior(Normal(doors, cov)))

# tem = 2.0*randn(1,N)+getVal(v1)+50.0
global v2 = addVariable!(fg,:x2, ContinuousScalar, N=N)
addFactor!(fg,[v1;v2],Odo([50.0]',[2.0]',[1.0]))

global v3=addVariable!(fg,:x3,ContinuousScalar, N=N)
addFactor!(fg,[v2;v3],Odo([50.0]',[4.0]',[1.0]))

global v4=addVariable!(fg,:x4,ContinuousScalar, N=N)
addFactor!(fg,[v3;v4],Odo([50.0]',[2.0]',[1.0]))

global v5=addVariable!(fg,:x5,ContinuousScalar, N=N)
addFactor!(fg,[v4;v5],Odo([50.0]',[2.0]',[1.0]))

global l1=addVariable!(fg, :l1, ContinuousScalar, N=N)
addFactor!(fg, [v3,l1], Ranged([64.0],[0.5],[1.0]))

global l2=addVariable!(fg, :l2, ContinuousScalar, N=N)
addFactor!(fg, [v4,l2], Ranged([16.0],[0.5],[1.0]))

global l3=addVariable!(fg, :l3, ContinuousScalar, N=N)
addFactor!(fg, [v5,l3], Ranged([16.0],[0.5],[1.0]))

# writeGraphPdf(fg);

global tree = wipeBuildNewTree!(fg, drawpdf=true);
# global tree = prepBatchTree!(fg,drawpdf=true);

# run(`evince bt.pdf`)

# do belief propagation inference over tree once
# inferOverTreeR!(fg, tree, N=100)