text:
# -- INTERNAL METHODS BELOW ---------------------------------------------------------- #
function precompute_translation_parameter_array_bacteria(biophysical_dictionary::Dict{String,Any}, protein_coding_length_array::Array{Float64,1},  time_constant_modifier_array::Array{Float64,1})::Array{TranslationParameters,1}

    # get some stuff from the biophysical_dictionary -
    mass_of_single_cell = biophysical_dictionary["mass_of_single_cell"]
    fraction_of_water_per_cell = biophysical_dictionary["fraction_of_water_per_cell"]
    ribosome_copy_number = biophysical_dictionary["copies_of_ribosome_per_cell"]
    translation_elongation_rate = biophysical_dictionary["translation_elongation_rate"]
    characteristic_initiation_time_translation = biophysical_dictionary["characteristic_initiation_time_translation"]
    KL = biophysical_dictionary["translation_saturation_constant"]
    species_symbol_type_array = biophysical_dictionary["species_symbol_type_array"]

    # fraction of dry weight -
    fraction_dry_cell = 1 - fraction_of_water_per_cell      # dimensionless

    # avagodros number -
    av_number = 6.02e23                                     # number/mol

    # what is the ribosome concentration -
    ribosome_concentration = ribosome_copy_number*(1/mass_of_single_cell)*(1/av_number)*(1e9)*(1/fraction_dry_cell) # nmol/gDW

    # which indices are protein -
    idx_protein = findall(x->x==:protein,species_symbol_type_array)

    # initialize -
    translation_parameters_array = Array{TranslationParameters,1}()

    # compute the translation parameters -
    for (index,protein_length) in enumerate(protein_coding_length_array)

        # compute kE -
        kE = translation_elongation_rate*(1/protein_length)

        # compute kI -
        kI = (1/characteristic_initiation_time_translation)

        # what is the correct index?
        correction_index = idx_protein[index]
        time_scale_parameter = time_constant_modifier_array[correction_index]

        # compute the tau factor -
        tau_factor = (kE/kI)*time_scale_parameter

        # compute the vmax -
        translation_vmax = kE*ribosome_concentration

        # build -
        translationParameters = TranslationParameters()
        translationParameters.vmax = translation_vmax*(3600)    # convert to hr
        translationParameters.tau_factor = tau_factor
        translationParameters.KL = KL
        push!(translation_parameters_array, translationParameters)
    end

    return translation_parameters_array
end

function precompute_translation_parameter_array_cell_free(biophysical_dictionary::Dict{String,Any}, protein_coding_length_array::Array{Float64,1},  time_constant_modifier_array::Array{Float64,1})::Array{TranslationParameters,1}

    # get some stuff from the biophysical_dictionary -
    ribosome_concentration = biophysical_dictionary["ribosome_concentration"]
    translation_elongation_rate = biophysical_dictionary["translation_elongation_rate"]
    characteristic_initiation_time_translation = biophysical_dictionary["characteristic_initiation_time_translation"]
    KL = biophysical_dictionary["translation_saturation_constant"]
    species_symbol_type_array = biophysical_dictionary["species_symbol_type_array"]

    # initialize -
    translation_parameters_array = Array{TranslationParameters,1}()

    # which indices are protein -
    idx_protein = findall(x->x==:protein,species_symbol_type_array)

    # compute the translation parameters -
    for (index,protein_length) in enumerate(protein_coding_length_array)

        # compute kE -
        kE = translation_elongation_rate*(1/protein_length)

        # compute kI -
        kI = (1/characteristic_initiation_time_translation)

        # what is the correct index?
        correction_index = idx_protein[index]
        time_scale_parameter = time_constant_modifier_array[correction_index]

        # compute the tau factor -
        tau_factor = (kE/kI)*time_scale_parameter

        # compute the vmax -
        # use a polysome factor of 6 -
        translation_vmax = kE*ribosome_concentration

        # build -
        translationParameters = TranslationParameters()
        translationParameters.vmax = translation_vmax*(3600)    # convert to hr
        translationParameters.tau_factor = tau_factor
        translationParameters.KL = KL
        push!(translation_parameters_array, translationParameters)
    end

    return translation_parameters_array
end

function precompute_translation_parameter_array_mammalian(biophysical_dictionary::Dict{String,Any}, protein_coding_length_array::Array{Float64,1},  time_constant_modifier_array::Array{Float64,1})::Array{TranslationParameters,1}

    # get some stuff from the biophysical_dictionary -
    mass_of_single_cell = biophysical_dictionary["mass_of_single_cell"]
    fraction_of_water_per_cell = biophysical_dictionary["fraction_of_water_per_cell"]
    ribosome_copy_number = biophysical_dictionary["copies_of_ribosome_per_cell"]
    translation_elongation_rate = biophysical_dictionary["translation_elongation_rate"]
    characteristic_initiation_time_translation = biophysical_dictionary["characteristic_initiation_time_translation"]
    KL = biophysical_dictionary["translation_saturation_constant"]
    species_symbol_type_array = biophysical_dictionary["species_symbol_type_array"]

    # fraction of dry weight -
    fraction_dry_cell = 1 - fraction_of_water_per_cell      # dimensionless

    # avagodros number -
    av_number = 6.02e23                                     # number/mol

    # what is the ribosome concentration -
    ribosome_concentration = ribosome_copy_number*(1/mass_of_single_cell)*(1/av_number)*(1e9)*(1/fraction_dry_cell) # nmol/gDW

    # which indices are protein -
    idx_protein = findall(x->x==:protein,species_symbol_type_array)

    # initialize -
    translation_parameters_array = Array{TranslationParameters,1}()

    # compute the translation parameters -
    for (index,protein_length) in enumerate(protein_coding_length_array)

        # compute kE -
        kE = translation_elongation_rate*(1/protein_length)

        # compute kI -
        kI = (1/characteristic_initiation_time_translation)

        # what is the correct index?
        correction_index = idx_protein[index]
        time_scale_parameter = time_constant_modifier_array[correction_index]

        # compute the tau factor -
        tau_factor = (kE/kI)*time_scale_parameter

        # compute the vmax -
        translation_vmax = kE*ribosome_concentration

        # build -
        translationParameters = TranslationParameters()
        translationParameters.vmax = translation_vmax*(3600)    # convert to hr
        translationParameters.tau_factor = tau_factor
        translationParameters.KL = KL
        push!(translation_parameters_array, translationParameters)
    end

    return translation_parameters_array
end

function precompute_translation_parameter_array(biophysical_dictionary::Dict{String,Any}, protein_coding_length_array::Array{Float64,1},  time_constant_modifier_array::Array{Float64,1}, host_type::Symbol = :bacteria)::Array{TranslationParameters,1}

    # we need to check - which host_type are we?
    if host_type == :bacteria
        return precompute_translation_parameter_array_bacteria(biophysical_dictionary::Dict{String,Any}, protein_coding_length_array::Array{Float64,1},  time_constant_modifier_array::Array{Float64,1})
    elseif host_type == :mammalian
        return precompute_translation_parameter_array_mammalian(biophysical_dictionary::Dict{String,Any}, protein_coding_length_array::Array{Float64,1},  time_constant_modifier_array::Array{Float64,1})
    elseif host_type == :cell_free
        return precompute_translation_parameter_array_cell_free(biophysical_dictionary::Dict{String,Any}, protein_coding_length_array::Array{Float64,1},  time_constant_modifier_array::Array{Float64,1})
    else
        throw(error("unsupported host_type value. Expected {:bacteria,:mammalian,:cell_free}. Got $(host_type)"))
    end
end

function precompute_transcription_kinetic_limit_array_cell_free(biophysical_dictionary::Dict{String,Any}, gene_coding_length_array::Array{Float64,1}, gene_abundance_array::Array{Float64,1},  time_constant_modifier_array::Array{Float64,1})::Array{Float64,1}

    # get stuff from the biophysical dictionary -
    characteristic_length = biophysical_dictionary["characteristic_transcript_length"]
    transcription_elongation_rate = biophysical_dictionary["transcription_elongation_rate"]
    characteristic_initiation_time = biophysical_dictionary["characteristic_initiation_time_transcription"]
    KX = biophysical_dictionary["transcription_saturation_constant"]
    RNAPII_concentration = biophysical_dictionary["RNAPII_concentration"]
    species_symbol_type_array = biophysical_dictionary["species_symbol_type_array"]

    # initialize -
    transcription_kinetics_array = Float64[]

    # which indices are mRNA (the production of transcription)?
    idx_mrna = findall(x->x==:mrna,species_symbol_type_array)

    # compute the kinetic limit -
    for (gene_index, gene_length) in enumerate(gene_coding_length_array)

        # how much gene do we have?
        # for host_type == :cell_free this is in concentration units -
        G = gene_abundance_array[gene_index]

        # what is the correction index?
        correction_index = idx_mrna[gene_index]
        time_scale_parameter = time_constant_modifier_array[correction_index]

        # compute kE -
        kE = transcription_elongation_rate*(1/gene_length)

        # compute kI -
        kI = (1/characteristic_initiation_time)

        # compute the tau factor -
        tau_factor = (kE/kI)*time_scale_parameter

        # Compute the rate -
        sat_term = (G/(KX*tau_factor+(1+tau_factor)*G))
        value = kE*RNAPII_concentration*sat_term*(3600) # nmol/gDW-hr

        # push -
        push!(transcription_kinetics_array, value)
    end

    # return -
    return transcription_kinetics_array
end

function precompute_transcription_kinetic_limit_array_bacteria(biophysical_dictionary::Dict{String,Any}, gene_coding_length_array::Array{Float64,1}, gene_copy_number_array::Array{Float64,1},  time_constant_modifier_array::Array{Float64,1})::Array{Float64,1}

    # get stuff from the biophysical dictionary -
    RNAPII_copy_number = biophysical_dictionary["copies_of_rnapII_per_cell"]
    characteristic_length = biophysical_dictionary["characteristic_transcript_length"]
    transcription_elongation_rate = biophysical_dictionary["transcription_elongation_rate"]
    characteristic_initiation_time = biophysical_dictionary["characteristic_initiation_time_transcription"]
    KX = biophysical_dictionary["transcription_saturation_constant"]
    mass_of_single_cell = biophysical_dictionary["mass_of_single_cell"]
    fraction_of_water_per_cell = biophysical_dictionary["fraction_of_water_per_cell"]
    species_symbol_type_array = biophysical_dictionary["species_symbol_type_array"]

    # fraction of dry weight -
    fraction_dry_cell = 1 - fraction_of_water_per_cell      # dimensionless

    # avagodros number -
    av_number = 6.02e23                                     # number/mol

    # what is the RNAP concentration -
    RNAPII_concentration = RNAPII_copy_number*(1/mass_of_single_cell)*(1/av_number)*(1e9)*(1/fraction_dry_cell) # nmol/gDW

    # which indices are mRNA (the production of transcription)?
    idx_mrna = findall(x->x==:mrna,species_symbol_type_array)

    # initialize -
    transcription_kinetics_array = Float64[]

    # compute the kinetic limit -
    for (gene_index, gene_length) in enumerate(gene_coding_length_array)

        # compute elongation constant -
        gene_copy_number = gene_copy_number_array[gene_index]

        # compute kE -
        kE = transcription_elongation_rate*(1/gene_length)

        # compute kI -
        kI = (1/characteristic_initiation_time)

        # what is the correction index?
        correction_index = idx_mrna[gene_index]
        time_scale_parameter = time_constant_modifier_array[correction_index]

        # compute the tau factor -
        tau_factor = (kE/kI)*time_scale_parameter

        # compute the gene concentration -
        G = gene_copy_number*(1/mass_of_single_cell)*(1/av_number)*(1e9)*(1/fraction_dry_cell)   # nmol/gDW

        # Compute the rate -
        value = kE*RNAPII_concentration*(G/(KX*tau_factor+(1+tau_factor)*G))*(3600) # nmol/gDW-hr

        # push -
        push!(transcription_kinetics_array, value)
    end

    # return -
    return transcription_kinetics_array
end

function precompute_transcription_kinetic_limit_array_mammalian(biophysical_dictionary::Dict{String,Any}, gene_coding_length_array::Array{Float64,1}, gene_copy_number_array::Array{Float64,1},  time_constant_modifier_array::Array{Float64,1})::Array{Float64,1}

    # get stuff from the biophysical dictionary -
    RNAPII_copy_number = biophysical_dictionary["copies_of_rnapII_per_cell"]
    characteristic_length = biophysical_dictionary["characteristic_transcript_length"]
    transcription_elongation_rate = biophysical_dictionary["transcription_elongation_rate"]
    characteristic_initiation_time = biophysical_dictionary["characteristic_initiation_time_transcription"]
    KX = biophysical_dictionary["transcription_saturation_constant"]
    mass_of_single_cell = biophysical_dictionary["mass_of_single_cell"]
    fraction_of_water_per_cell = biophysical_dictionary["fraction_of_water_per_cell"]
    species_symbol_type_array = biophysical_dictionary["species_symbol_type_array"]

    # fraction of dry weight -
    fraction_dry_cell = 1 - fraction_of_water_per_cell      # dimensionless

    # avagodros number -
    av_number = 6.02e23                                     # number/mol

    # what is the RNAP concentration -
    RNAPII_concentration = RNAPII_copy_number*(1/mass_of_single_cell)*(1/av_number)*(1e9)*(1/fraction_dry_cell) # nmol/gDW

    # initialize -
    transcription_kinetics_array = Float64[]

    # which indices are mRNA (the production of transcription)?
    idx_mrna = findall(x->x==:mrna,species_symbol_type_array)

    # compute the kinetic limit -
    for (gene_index, gene_length) in enumerate(gene_coding_length_array)

        # compute elongation constant -
        gene_copy_number = gene_copy_number_array[gene_index]

        # compute kE -
        kE = transcription_elongation_rate*(1/gene_length)

        # compute kI -
        kI = (1/characteristic_initiation_time)

        # what is the correction index?
        correction_index = idx_mrna[gene_index]
        time_scale_parameter = time_constant_modifier_array[correction_index]

        # compute the tau factor -
        tau_factor = (kE/kI)*time_scale_parameter

        # compute the gene concentration -
        G = gene_copy_number*(1/mass_of_single_cell)*(1/av_number)*(1e9)*(1/fraction_dry_cell)   # nmol/gDW

        # Compute the rate -
        value = kE*RNAPII_concentration*(G/(KX*tau_factor+(1+tau_factor)*G))*(3600) # nmol/gDW-hr

        # push -
        push!(transcription_kinetics_array, value)
    end

    # return -
    return transcription_kinetics_array
end

function precompute_transcription_kinetic_limit_array(biophysical_dictionary::Dict{String,Any}, gene_coding_length_array::Array{Float64,1}, gene_abundance_array::Array{Float64,1}, time_constant_modifier_array::Array{Float64,1}, host_type::Symbol)::Array{Float64,1}

    # swicth on the host_type flag -
    if host_type == :bacteria
        return precompute_transcription_kinetic_limit_array_bacteria(biophysical_dictionary,gene_coding_length_array,gene_abundance_array, time_constant_modifier_array)
    elseif host_type == :mammalian
        return precompute_transcription_kinetic_limit_array_mammalian(biophysical_dictionary,gene_coding_length_array,gene_abundance_array, time_constant_modifier_array)
    elseif host_type == :cell_free
        return precompute_transcription_kinetic_limit_array_cell_free(biophysical_dictionary, gene_coding_length_array, gene_abundance_array, time_constant_modifier_array)
    else
        throw(error("unsupported host_type value. Expected {:bacteria,:mammalian,:cell_free}. Got $(host_type)"))
    end
end

function calculate_translation_kinetics_array(t,x,mRNA_index_array,translation_parameters_array, data_dictionary)

    # we need to calculate the translation rate on the fly
    # since it depends upon mRNA -
    translation_kinetics_array = Float64[]

    # grab the translation parameters from the data dictionary -
    number_of_translation_rates = length(translation_parameters_array)
    counter = 1
    for index = 1:number_of_translation_rates

        # use the mRNA_index_array to index into the state vector -
        mRNA_concentration = x[mRNA_index_array[index]]

        # grab the parameters struct -
        parameters_struct = translation_parameters_array[index]
        vmax = parameters_struct.vmax
        tau = parameters_struct.tau_factor
        KL = parameters_struct.KL

        # build the rate -
        value = vmax*(mRNA_concentration/(KL*tau+(1+tau)*mRNA_concentration))

        # package -
        push!(translation_kinetics_array, value)
    end

    # return -
    return translation_kinetics_array
end

function build_dilution_degradation_matrix(biophysical_dictionary::Dict{String,Any}, species_type_array::Array{Symbol,1}, degradation_modifier_array::Array{Float64,1})::Array{Float64,2}

    # how many species do we have?
    total_number_of_species = length(species_type_array)

    # estimate *max* growth -
    # get the cell doubling time -
    cell_doubling_time = biophysical_dictionary["cell_doubling_time"]
    mugmax = 0.0
    if cell_doubling_time != 0.0

        # compute the mumax -
        mugmax = log(2)/cell_doubling_time
    end

    # calculate the mRNA and protein half life -
    mRNA_half_life_in_hr = biophysical_dictionary["mRNA_half_life_in_hr"]
    protein_half_life_in_hr = biophysical_dictionary["protein_half_life_in_hr"]

    # compute the degradation rate constant -
    kDX = -(1/mRNA_half_life_in_hr)*log(0.5)
    kDL = -(1/protein_half_life_in_hr)*log(0.5)

    # initialize array which holds the (mu+d) terms -
    dilution_degradation_term_array = Array{Float64,1}()
    for (species_index,species_type) in enumerate(species_type_array)

        # for species, compute a degradation term -
        if species_type == :gene || species_type == :constant
            push!(dilution_degradation_term_array,0.0)
        elseif species_type == :mrna

            # get degradation modifier -
            degradation_modifier_value = degradation_modifier_array[species_index]

            # compute -
            value = -1*(mugmax + degradation_modifier_value*kDX)

            # grab -
            push!(dilution_degradation_term_array,value)

        elseif species_type == :protein

            # get degradation modifier -
            degradation_modifier_value = degradation_modifier_array[species_index]

            # compute -
            value = -1*(mugmax + degradation_modifier_value*kDL)

            # grab -
            push!(dilution_degradation_term_array,value)
        end
    end

    # return the dilution array -
    return Diagonal(dilution_degradation_term_array)
end
# -- INTERNAL METHODS ABOVE ---------------------------------------------------------- #

# -- PUBLIC METHODS BELOW ------------------------------------------------------------ #
function calculate_txtl_kinetics_array(t::Float64,x::Array{Float64,1},data_dictionary::Dict{String,Any})::Array{Float64,1}

    # initialize -
    calculated_txtl_kinetics_array = Array{Float64,1}()

    # we have precomputed the translation kinetic parameters, and the kinetic limit of the transcription rate -
    transcription_kinetic_limit_array = data_dictionary["transcription_kinetic_limit_array"]
    translation_parameters_array = data_dictionary["translation_parameter_array"]
    species_symbol_type_array = data_dictionary["species_symbol_type_array"]

    # find mRNA and protein indices -
    mRNA_index_array = findall(x->x==:mrna,species_symbol_type_array)
    protein_index_array = findall(x->x==:protein,species_symbol_type_array)

    # compute the rate of translation -
    translation_kinetics_array = calculate_translation_kinetics_array(t,x,mRNA_index_array,translation_parameters_array, data_dictionary)

    # package -
    for (index,value) in enumerate(transcription_kinetic_limit_array)
        push!(calculated_txtl_kinetics_array,value)
    end

    for (index,value) in enumerate(translation_kinetics_array)
        push!(calculated_txtl_kinetics_array,value)
    end

    # return - TX rates, then TL rates -
    return calculated_txtl_kinetics_array
end
# -- PUBLIC METHODS ABOVE ------------------------------------------------------------ #


text:
using ZOOclient
# using PyPlot

# define a Dimension object
dim_size = 100
dim_regs = [[-1, 1] for i = 1:dim_size]
dim_tys = [true for i = 1:dim_size]
mydim = Dimension(dim_size, dim_regs, dim_tys)
# define an Objective object
obj = Objective(mydim)

# define a Parameter Object, the five parameters are indispensable.
# budget:  number of calls to the objective function
# evalueation_server_num: number of evaluation cores user requires
# control_server_ip_por4t: the ip:port of the control server
# objective_file: objective funtion is defined in this file
# func: name of the objective function
init_sample = construct_init_sample("/Users/liu/.julia/v0.6/ZOOclient/example/init.txt")
par = Parameter(budget=400, evaluation_server_num=2, control_server_ip_port="192.168.100.108:20000",
    objective_file="fx.py", func="ackley", output_file="log.txt", init_sample=init_sample, uncertain_bits=5)

# perform optimization
sol = zoo_min(obj, par)
# print the Solution object
sol_print(sol)
positive_data = take!(par.positive_data)
negative_data = take!(par.negative_data)
write_population("population.txt", positive_data, negative_data)
# visualize the optimization progress
# history = get_history_bestsofar(obj)
# plt[:plot](history)
# plt[:savefig]("figure.png")

ext: jl


text:
export ft_ee_uaread!

function ft_ee_uaread!(ft_handle::Culong, buffer::Array{UInt8,1})
  len = length(buffer)
  bytesread = Ref{Cuint}()
  ft_status = ccall((:FT_EE_UARead, d2xx),
                     Cuint,
                     (Culong, Ptr{UInt8}, Cuint, Ref{Cuint}),
                     ft_handle, buffer, len, bytesread)
  checkstatus(ft_status)
  return bytesread[]
end


text:
using DataDrivenDiffEq
using ModelingToolkit
using OrdinaryDiffEq
using LinearAlgebra
using Plots
gr()

# Create a basis of functions aka observables ( aka features )
@variables u[1:2]
h = [u[1]; sin(u[1]); sin(u[2]); u[2]; u[1]*u[2]; u[2]^2]
basis = Basis(h, u)

# Create a test system
function test_discrete(du, u, p, t)
    du[1] = 0.9u[1] + 0.1u[2]^2
    du[2] = sin(u[1]) - 0.1u[1]
end

# Set up the problem
u0 = [1.0; 2.0]
tspan = (0.0, 50.0)
prob = DiscreteProblem(test_discrete, u0, tspan)
sol = solve(prob, FunctionMap(), atol = 1e-8, rtol = 1e-8)
# Plot the solution
plot(sol)

# Build the edmd with a subset of measurements
approximator = EDMD(sol[:,:], basis)

# Let's look at the eigenvalues
scatter(eigvals(approximator))

# Solve the estimation problem
sys = ODESystem(approximator, threshold = 1e-3)
dudt = ODEFunction(sys)
prob_ = DiscreteProblem(dudt, u0, tspan)
sol_ = solve(prob_, FunctionMap(), atol = 1e-8, rtol = 1e-8)

# Plot the error
plot(sol_)
plot(sol.t, abs.(sol - sol_)')
norm(sol-sol_) # ≈ 2.52e-13

# Get the linear dynamics in the Koopman space
dψ(u, p, t) = operator(approximator)*u
dψ(du, u, p, t) = du .= operator(approximator)*u
ψ0 = basis(u0)
# Simply calling the EDMD struct transforms into the current basis
ψ_prob = DiscreteProblem(dψ, ψ0, tspan)
ψ = solve(ψ_prob, FunctionMap())

# Plot trajectory in edmd basis
plot(sol.t, ψ')
plot(sol.t, basis(Array(sol))')
norm((ψ - basis(Array(sol)))) # ≈ 0.0158

# And in observable space
sol_ψ = outputmap(approximator) * ψ
plot(abs.(sol_ψ'- sol[:,:]'))
norm(sol_ψ - sol) # ≈ 0.00664

# Reduce the basis
reduced_basis = reduce_basis(approximator)
parameters(reduced_basis)
# Reiterate for the reduced approximator with few measurements
reduced_approximator = EDMD(sol[:,1:4], reduced_basis)
operator(reduced_approximator)
# Simply calling the EDMD struct transforms into the current basis
println(reduced_basis)
ψ_prob = DiscreteProblem(reduced_approximator, u0, tspan, [])
ψ = solve(ψ_prob, FunctionMap())

# Plot trajectory in edmd basis with much higher accuracy than before
plot(ψ)
plot!(sol)
#plot!(sol.t, hcat([reduced_approximator(xi) for xi in eachcol(sol)]...)')
norm(ψ - Array(sol)) # ≈ 2.29..

# Update with all other measurements
update!(reduced_approximator, sol[:,8:10], sol[:, 9:11])
ψ_prob = DiscreteProblem(reduced_approximator, u0, tspan, [])
ψ = solve(ψ_prob, FunctionMap())
norm(ψ - Array(sol)) # ≈ 0.133

# Simply calling the EDMD struct transforms into the current basis
ψ_prob = DiscreteProblem(reduced_approximator, u0, tspan, [])
ψ = solve(ψ_prob, FunctionMap())
println(reduced_basis)
# Plot trajectory in edmd basis with much higher accuracy than before
plot(ψ)
plot!(sol)
norm(ψ .- sol) # ≈ 0.134...


text:
push!(LOAD_PATH, "/home/renee/phd/src/ebit-evolution.project/ebit-ode-server/")
include("/home/renee/phd/src/ebit-evolution.project/ebit-ode-server/EbitODEMessages.jl")
include("/home/renee/phd/src/ebit-evolution.project/ebit-ode-server/EbitSolver.jl")
include("/home/renee/phd/src/ebit-evolution.project/ebit-ode-server/EbitODEServer.jl")

using Revise
using EbitODEServer
import Base.+, Base.*, Base./, Base.-, Base.zero, Base.transpose, Base.exp, Base.^

function join_symbol(a, b)
     return string(string(a), string(b))
end

function check_empty(a, b, default)
    if (a == "" && b != "")
        return b
    elseif (b == "" && a != "")
        return a
    elseif (b == "" && a == "")
        return ""
    else
        return default
    end
end

function (+)(a::Any, b::Any)
    return check_empty(a, b, String("($a + $b)"))
end

function (*)(a::Any, b::Any)
    return check_empty(a, b, String("($a * $b)"))
end

function (/)(a::Any, b::Any)
    return check_empty(a, b, String("($a / $b)"))
end

function (-)(a::Any, b::Any)
    return check_empty(a, b, String("($a - $b)"))
end

function (-)(a::Any)
    return String("-$a)")
end

function (^)(a::Any, b::Any)
    return String("$a^{$b})")
end

function exp(a::Any)
    return String("exp[$a]")
end

transpose(str::String) = String(str)

function zero(a::String)
    return ""
end

function zero(a::Type{String}) "" end

function (^)(a::Any, p::Real)
    "$a^{$p}"
end

default_size = 2

function c1(symbol)
    ["$symbol\_$i" for i=1:default_size]
end

function c2(symbol)
    val = Array{String,2}(default_size,default_size)
    for i in 1:default_size
        for j in 1:default_size
            val[i,j] = "$symbol\_$i$j"
        end
    end
    return val
end

τ = c1(:τ)
A = c1(:A)
n = c1(:n)
Χ = c1(:Χ)
N = c1(:N)
CX = c2(:CX)

(τ.' .* CX)*N

(τ .+ τ')./A'

(τ./τ')*ones(2,1)

τ' .- τ

τ./A .+ (τ./A)'

qVₜ = c1(:qVₜ)

τ./qVₜ .+ (τ./qVₜ)'


text:
module ArgMinimist

export parse_args, parse_args!, parse_named_args, parse_named_args!

intRegex = r"^\d+$"
floatRegex = r"^\d*\.\d+$"
argNameRegex = r"^(\-*)(\w+)$"

function get_arg_value(arg::String)
  if ismatch(intRegex, arg)
    return parseint(arg)
  elseif ismatch(floatRegex, arg)
    return parsefloat(arg)
  else
    return arg
  end
end

function parse_args()
  parse_args(ARGS)
end

function parse_args!(keyword::Dict{Symbol,Any})
  parse_args!(ARGS, keyword)
end
function get_arg_name(arg::String)
  symbol(match(argNameRegex, arg).captures[2])
end

function parse_args{T<:String}(args::Vector{T})
  keyword = Dict{Symbol,Any}()
  parse_args!(args, keyword)
end

function parse_args!{T<:String}(args::Vector{T}, keyword::Dict{Symbol,Any})

  no_prefix = "--no"
  yes_prefix = "--with"

  positional = Any[]

  pos = 1

  while pos <= length(args)

    argp = args[pos]

    if beginswith(argp, "-")

      if beginswith(argp, "--")

        if contains(argp, "=")
          # This is a compound keyword argument
          argname, argval = split(argp,"=")
          keyword[ get_arg_name(argname) ] = get_arg_value(argval)
        else

          if beginswith(argp, no_prefix)
            keyword[ get_arg_name(replace(argp, no_prefix, "")) ] = false
          elseif beginswith(argp, yes_prefix) || pos == length(args) || beginswith(args[pos+1],"-")
            # This is a true flag if any of these are true:
            #   has the correct prefix
            #   is the last argument
            #   is followed by a keyword argument
            keyword[ get_arg_name(replace(argp, yes_prefix,"")) ] = true
          else
            # Ok, this is the name value and the next argument
            # is the actual argument value
            keyword[ get_arg_name(argp) ] = get_arg_value(args[pos+1])
            pos += 1
          end
        end

      else

        # Process as options to set
        for k in replace(argp,"-","")
          keyword[ symbol(k) ] = true
        end

      end
    else
      push!(positional, get_arg_value(argp))
    end

    pos += 1
  end

  return positional, keyword
end

function parse_named_args{T<:String}(names::Vector{T})
  parse_named_args(names, ARGS, Dict{Symbol,Any}())
end

function parse_named_args{T<:String}(names::Vector{T}, args::Vector{T})
  keyword = Dict{Symbol,Any}()
  parse_named_args!(names,args,keyword)
end

function parse_named_args!{T<:String}(names::Vector{T}, args::Vector{T}, keyword::Dict{Symbol,Any})

  positional, keyword = parse_args!(args, keyword)
  symbols = map(symbol, names)
  match_symbols = filter((s)-> !haskey(keyword, s), symbols)

  if length(match_symbols) > length(positional)

    println("\n\tExpecting arguments: ",join(names, " "))
    exit(-1)

  else

    for (s,v) in zip(match_symbols, positional)
      keyword[s] = v
    end

    if length(match_symbols) < length(positional)
      positional = positional[length(match_symbols)+1:end]
    else
      # Empty because all were matched
      positional = Any[]
    end
  end

  return positional, keyword
end

end # module


text:
using Catlab.Programs.RelationalPrograms: parse_relation_context, parse_relation_call, parse_relation_diagram
using Catlab.WiringDiagrams
using DataStructures
using IterTools
using Combinatorics

include(joinpath(@__DIR__, "FLS.jl"))
include(joinpath(@__DIR__, "DB.jl"))

const EqClass = Dict{Symbol, IntDisjointSets}
const ACDict = Dict{UInt64, Pair{StructACSet, ChaseStepData}}

"""
Do a chase step: cone gen then TGDs, EGDs for each branch
"""
function chasestep(F::FLS, J::StructACSet)::Pair{ACDict, ACDict}
  verbose = true
  res, completed, srcs = ACDict(), ACDict(), F.schema[:vlabel][F.schema[:src]]
  #if verbose println("\tabout to add cones") end
  J_, csd_ = cone_gen(F, J) # modifies csd_ only - adds cones to copy
  ats = apply_tgds(F, J_, csd_) # make all possible choices for functions
  n = length(ats)
  for (i, (x, csd)) in enumerate(ats) # could be done in parallel
    e = apply_egds!(F, x, csd) # modifies csd
    merge!(F, x, e) # modifies x

    # Sanity checks
    is_iso = is_isomorphic(x, J)
    is_iso == isempty(csd) || error("iso csd $J \n$x \n$csd")
    if !is_iso
      if verbose println("\t$i/$n incomplete: hashing $(generate_json_acset(x))") end
      xhsh = canonical_hash(x)
      if verbose println("\thashed!") end
      res[xhsh] = x => csd
    else
      !nontotal(F, x) || error("non total model")
      if verbose println("\t$i/$n complete! hashing $(generate_json_acset(x))") end
      completed[canonical_hash(x)] = x => csd
    end
  end

  return res => completed
end

"""
Take a premodel from the DB and generate a set of models and premodels.
Insert these results into the DB. Return the IDs of the resulting premodels and
models.
"""
function chasestep_db(db::LibPQ.Connection, F::FLS, i::Int
                     )::Pair{Vector{Int}, Vector{Int}}
  # check if we've already chased before
  if columntable(execute(db, "SELECT fired FROM Premodel WHERE Premodel_id=\$1",
                             [i]))[:fired][1]
    children = collect(columntable(execute(db,
      "SELECT child FROM Parent WHERE parent=\$1",[i]))[:child])
    if isempty(children)
      return Int[] => Int[]
    else
      modelchildren = columntable(execute(db,
        """SELECT Model_id, Premodel_id FROM Premodel JOIN Model
           USING (Premodel_id) WHERE Premodel_id IN ($(join(children, ",")))"""))
    end
    pmcs, mcs = [collect(modelchildren[x]) for x in [:premodel_id, :model_id]]
    # println("\tALREADY CHASED! model_ids: $mcs")
    return setdiff(children, pmcs) => mcs
  end
  m, I = get_premodel(db, i)
  msize = [k => nparts(I, k) for k in realobs(F)]
  println("NEW PREMODEL TO CHASE: #$i ($msize)")

  xs, ys = chasestep(m, I)
  pminds, minds = Int[], Int[]
  for (xhsh, (x, csd)) in xs
    # println("\tNew premodel $([k => nparts(x, k) for k in realobs(F)])")
    push!(pminds, add_premodel(db, F, x, i, csd, xhsh))
  end
  for (yhsh, (y, csd)) in ys
    println("\tNew model $([k => nparts(y, k) for k in realobs(F)])")
    push!(minds, add_model(db, F, y, i, csd, yhsh))
  end
  # println("MARKING AS FIRED")
  execute(db, "UPDATE Premodel SET fired = true WHERE Premodel_id = \$1", [i])
  return pminds => minds
end

"""
1. Enumerate elements of ℕᵏ for an underlying graph with k nodes.
2. For each of these: (c₁, ..., cₖ) create a term model with that many constants

Do the first enumeration by incrementing n_nonzero and finding partitions so
that ∑(c₁,...) = n_nonzero

In the future, this function will write results to a database
that hashes the FLS as well as the set of constants that generated the model.

Also crucial is to decompose FLS into subparts that can be efficiently solved
and have solutions stitched together.
"""
function combos_below(m::Int, n::Int)::Vector{Vector{Int}}
  res = Set{Vector{Int}}([zeros(Int,m)])
  n_const = 0 # total number of constants across all sets
  for n_const in 1:n
    for n_nonzero in 1:m
      # values we'll assign to nodes
      c_parts = partitions(n_const, n_nonzero)
      # Which nodes we'll assign them to
      indices = permutations(1:m,n_nonzero)
      for c_partition in c_parts
        for index_assignment in indices
          v = zeros(Int, m)
          v[index_assignment] = vcat(c_partition...)
          push!(res, v)
        end
      end
    end
  end
  return sort(collect(res))
end

function nontotal(F::FLS, J::ACSet)::Bool
  for (src_i, e) in zip(F.schema[:src], F.schema[:elabel])
    srcobj, src = F.schema[:vlabel][src_i], add_srctgt(e)[1]
    missin = setdiff(parts(J, srcobj), J[src])
    if !isempty(missin)
      println("missing $srcobj <- $e: $missin")
      return true
    end
  end
  return false
end

"""
Find all models below a certain cardinality.

ignore_seen skips checking things in the database that were already chased.
If true, the final list of models may be incomplete, but it could be more
efficient if the goal of calling this function is merely to make sure all models
are in the database itself.
"""
function chase_below(db::LibPQ.Connection, F::FLS, n::Int; ignore_seen=false)::Vector{Int}
  unchased, res  = Set{Int}(), Set{Int}()
  m = length(realobs(F))
  for combo in combos_below(m, n)
    push!(unchased, init_premodel(db, F, combo))
  end
  seen = ignore_seen ? get_seen(db, F) : Set{Int}()
  while !isempty(unchased)
    mdl = pop!(unchased)
    m = get_premodel(db, mdl)[2]
    msize = [k => nparts(m, k) for k in realobs(F)]
    # println("chasing $mdl with size $msize")
    Ks, Ms = chasestep_db(db, F, mdl)
    # println("CHASESTEP RES: Ks=$Ks Ms=$Ms")
    function getidsize(x::Int)::Pair{UInt64, Int}
      z = columntable(execute(db,
        "SELECT hash,size FROM Premodel WHERE Premodel_id=\$1",[x]))
      return z[:hash][1] => z[:size][1]
    end
    filter!(k->let (i,s)=getidsize(k); i ∉ seen && (s <= n) end, Ks)

    push!(seen, mdl)
    union!(unchased, setdiff(Ks, seen))
    union!(res, Ms)
  end
  return collect(res)
end

"""
Given a particular arrow and source for that arrow, find all targets (including
"free" = 0 for nonlimit targets) that do not immediately violate path eqs.
"""
function valid_poss(J::StructACSet, F::FLS, e::Symbol, x::Int)::Vector{Int}
  res = tgt(F, e) ∈ realobs(F) ? [0] : Int[]
  isempty(incident(J, x, add_srctgt(e)[1])) || error(
    "checking possiblities for something that's already determined")
  esrc, etgt = add_srctgt(e)
  temprel = add_part!(J, e; Dict([esrc=>x])...)
  # we could precompute this before everything to get Symbol -> [Eq] dict
  relevant_eqs = collect(filter(pq->e ∈ vcat(pq[2],pq[3]), F.eqs))
  for tgt_val in parts(J,  tgt(F, e))
    #println("tgtval $tgt_val")
    set_subpart!(J, temprel, etgt, tgt_val)
    valid_tgt_val = true
    for (n,p,q) in relevant_eqs
      #println("\teq $n w/ p=$p ||| q=$q")
      if valid_tgt_val
        for strt in parts(J, src(F, p[1]))
          pres, qres = [
            let r=eval_path(x, J, Set([strt])); isempty(r) ? nothing : only(r)
            end for x in [p,q]]
          # length(pres) == 1 && length(qres) == 1 || error(
          #   "tgt val $tgt_val \np $p \nq $q\nstrt $strt pres $pres qres $qres")
          #pres, qres = collect(pres)[1], collect(qres)[1]
          if pres isa Int && qres isa Int
            if pres != qres
              #println("\t\tstart at $strt, p(x)=$pres q(x)=$qres")
            end
            valid_tgt_val &= pres == qres
          end
        end
      end
    end
    if valid_tgt_val
      push!(res, tgt_val)
    end
  end
  rem_part!(J, e, temprel)
  return res
end

"""
Add new elements when a function is dangling. Fire TGDs in all possible ways
(plus the 'free' way which uses fresh variables). This produces new premodels to
chase.
"""
function apply_tgds(F::FLS, J::StructACSet, csd::ChaseStepData
                    )::Vector{Pair{StructACSet, ChaseStepData}}
  res = Pair{StructACSet, ChaseStepData}[]
  srcs, tgts = [F.schema[:vlabel][F.schema[x]] for x in [:src, :tgt]]
  # Calculate how many missing things we need of each type
  needs = [(d, nparts(J, dtgt),
            collect(setdiff(parts(J, ddom), J[add_srctgt(d)[1]])))
           for (d, ddom, dtgt) in zip(F.schema[:elabel], srcs, tgts)]
  lens = [length(x) for (_,_,x) in needs]
  csum = Dict(zip(F.schema[:elabel],zip(lens,vcat([0],cumsum(lens)[1:end-1]))))
  prod_args = vcat([[valid_poss(J, F, e, x) for x in xs] for (e, _, xs) in needs]...)
  println("needs $needs\nprod_args $prod_args\n")
  # branch on cartesian product of possibilities. 0 means fresh.
  if isempty(prod_args)
    return [J => csd]
  end
  for (i, choice) in enumerate(IterTools.product(prod_args...))
    Jc, Jcsd = deepcopy(J), deepcopy(csd)
    for (e, src, tgt, (_,_,eneed)) in zip(F.schema[:elabel], srcs, tgts, needs)
      esrc, etgt = add_srctgt(e)
      len, offset = csum[e]
      echoices = len == 0 ? [] : choice[1+offset:offset+len]
      for (x, echoice) in zip(eneed, echoices)
        if echoice == 0
          gx = add_part!(Jc, tgt) # add a fresh symbol
          push!(Jcsd.tgds[e], x => gx)
          add_part!(Jc, e; Dict([esrc=>x, etgt=>gx])...)
        else
          add_part!(Jc, e; Dict([esrc=>x, etgt=>echoice])...)
          push!(Jcsd.tgds[e], x => echoice)
        end
      end
    end
    # make sure no limit constraints violated
    conedata, coneviol = check_cones(F, Jc), false
    for d in values(conedata)
      for lvs in map(length,values(d))
        coneviol |= lvs > 1
      end
    end
    if !coneviol
      push!(res, Jc=>Jcsd)
    end
  end
  return res
end

"""
Identify cones by their leg values. Multiple cones that share leg values are
grouped together.
"""
function check_cones(F::FLS, J::StructACSet,
                     eqclass::Union{Nothing, EqClass}=nothing
                    )::Dict{Symbol, Dict{Vector{Int}, Vector{Int}}}
  eqclasses = eqclass === nothing ? init_eq(F,J) : eqclass
  return Dict([c.apex=>check_cone(F, J, c, eqclasses) for c in F.cones])
end

function check_cone(F::FLS, J::StructACSet, cone::Cone, eqclass::EqClass
                    )::Dict{Vector{Int}, Vector{Int}}
  cones = DefaultDict{Vector{Int}, Vector{Int}}(Vector{Int})
  ltgts = [map(i->find_root!(eqclass[tgt(F,e)],i), J[add_srctgt(e)[2]])
           for (_, e) in cone.legs]
  for (i, v) in enumerate(map(collect,zip(ltgts...)))
    if minimum(v) > 0
      push!(cones[v], i)
    end
  end
  return cones
end

"""Add new elements required by cones. Modifies J."""
function cone_gen(F::FLS, J_::StructACSet)
  J = deepcopy(J_)
  csd = ChaseStepData()

  verbose = false
  newstuff = Pair{Symbol, Vector{Pair{Symbol, Int}}}[]
  conedata = check_cones(F, J)
  for cone in F.cones # could be done in parallel
    cones = conedata[cone.apex]
    if verbose && !isempty(cones) println("preexisting cones $(cones)") end
    # look for instances of the pattern (TODO cached version / only considered 'added' information)
    for res in query(J, cone_query(cone)) # rely on order by same as order when precomputing
      resv = Vector{Int}(collect(res))
      if verbose println("res $res $resv") end
      if !haskey(cones, resv)
        if verbose println("$(cone.apex) is new w/ $(nparts(J, cone.apex))") end
        newargs = [e=>val for (val, (_, e)) in zip(resv,cone.legs)]
        push!(newstuff, cone.apex => newargs)
      end
    end
  end
  # Add to J
  for (apex, args) in newstuff
    newapex = add_part!(J, apex)
    push!(csd.cones[apex], newapex)
    for (e, tgtval) in args
      s, t = add_srctgt(e)
      add_part!(J, e; Dict([s=>newapex, t=>tgtval])...)
    end
  end
  return J, csd
end

"""
Compose relations starting with a subset of the composite domain to get a subset
of the (composite) codomain
"""
function eval_path(pth::Vector{Symbol}, J::StructACSet, xs::Set{Int})::Set{Int}
  # Iteratively update X to Y by searching through all (x,y) in R ⊆ X×Y
  for m in pth
    msrc, mtgt = add_srctgt(m)
    if isempty(xs)
      return xs
    else
      newxs = Set{Int}()
      for (jx, jy) in zip(J[msrc], J[mtgt])
        if jx ∈ xs
          push!(newxs, jy)
        end
      end
      xs = newxs
    end
  end
  return xs
end

function cone_eqs!(F::FLS, J::StructACSet,csd::ChaseStepData,eqclass::EqClass)::Nothing
  conedata = check_cones(F, J, eqclass)
  for cone in F.cones # could be done in parallel
    cones = conedata[cone.apex]
    # set equivalences
    for vs in filter(x->length(x) > 1, collect(values(cones)))
      push!(csd.cone_eqs[cone.apex], vs)
      for v in vs[2:end]
        union!(eqclass[cone.apex], vs[1], v)
      end
    end
  end
end

function check_path_eq!(F::FLS, J::StructACSet,csd::ChaseStepData,eqclasses::EqClass, name::Symbol)::Nothing
  p, q = get_eq(F, name)
  eqc = eqclasses[tgt(F, p[end])]
  for x in parts(J, src(F, p[1]))
    for px in eval_path(p, J, Set([x]))
      for qx in (isempty(q) ? [x] : eval_path(q, J, Set([x])))
        if !in_same_set(eqc,px,qx)
          # println("name $name p $p q $q\n x $x px $px qx $qx \n $J")
          push!(csd.path_eqs[name], px => qx)
          union!(eqc, px, qx)
        end
      end
    end
  end
end

function path_eqs!(F::FLS, J::StructACSet,csd::ChaseStepData,eqclasses::EqClass)::Nothing
  for (name, p, q) in F.eqs # could be done in parallel
    eqc = eqclasses[tgt(F, p[end])]
    for x in parts(J, src(F, p[1]))
      for px in eval_path(p, J, Set([x]))
        for qx in (isempty(q) ? [x] : eval_path(q, J, Set([x])))
          if !in_same_set(eqc,px,qx)
            # println("name $name p $p q $q\n x $x px $px qx $qx \n $J")
            push!(csd.path_eqs[name], px => qx)
            union!(eqc, px, qx)
          end
        end
      end
    end
  end
end

"""
Note which elements are equal due to relations actually representing functions
"""
function fun_eqs!(F::FLS, J::StructACSet, csd::ChaseStepData, eqclass::EqClass
                 )::Nothing
  for d in F.schema[:elabel] # could be done in parallel
    dsrc, dtgt = add_srctgt(d)
    srcobj, tgtobj = src(F, d), tgt(F,d)
    eq_ys = DefaultDict{Int, Set{Int}}(Set{Int})
    for (src_i, tgt_i) in zip(J[dsrc], J[dtgt])
      push!(eq_ys[find_root!(eqclass[srcobj], src_i)], tgt_i)
    end
    for ys in filter(x->length(x)>1, collect(values(eq_ys)))
        push!(csd.fun_eqs[d], collect(ys))
        ymin = minimum(ys)
        for y in filter(!=(ymin), ys)
          union!(eqclass[tgtobj], y, ymin)
        end
    end
  end
end

function init_eq(F::FLS, J::StructACSet)::EqClass
  EqClass([o=>IntDisjointSets(nparts(J, o)) for o in F.schema[:vlabel]])
end
function same_eq_classes(x::IntDisjointSets,y::IntDisjointSets)::Bool
  length(x)==length(y) || error("same_eq_class called on different lengths")
  for i in 1:length(x)
    for j in 1:length(x)
      if in_same_set(x, i, j) != in_same_set(y, i, j)
        return false
      end
    end
  end
  return true
end

"""Modifies eq"""
function apply_egds!(F::FLS, J::StructACSet,csd::ChaseStepData)::EqClass
  # println("$i/$(length(ats)) applying egds to $(generate_json_acset(x))")
  eqclasses, changed = init_eq(F,J), true
  while changed
    start_eq = deepcopy(eqclasses)
    cone_eqs!(F,J,csd,eqclasses)
    path_eqs!(F,J,csd,eqclasses)
    fun_eqs!(F,J,csd,eqclasses)
    changed = !all([same_eq_classes(start_eq[o], eqclasses[o])
                    for o in F.schema[:vlabel]])
  end
  return eqclasses
end

"""Use equivalence class data to reduce size of model"""
function merge!(F::FLS, J::StructACSet, eqclasses::EqClass;
                verbose::Bool=false)::Nothing

  μ = Dict{Symbol, Vector{Pair{Int,Int}}}([
    o=>Pair{Int,Int}[] for o in F.schema[:vlabel]])
  delob = DefaultDict{Symbol, Vector{Int}}(Vector{Int})
  for (o, eq) in pairs(eqclasses)
    eqsets = Dict{Int,Set{Int}}()
    # Recover the equivalence classes from the union-find structure
    for x in parts(J, o)
      eqrep = find_root(eq, x)
      if haskey(eqsets, eqrep)
        push!(eqsets[eqrep], x)
      else
        eqsets[eqrep] = Set([x])
      end
    end
    filter!(p->length(last(p))>1, eqsets)
    if verbose && !isempty(eqsets)
      println("eqset $o $(values(eqsets))")
    end
    # Minimum element is the representative
    for vs in map(collect,values(eqsets))
      m = minimum(vs)
      vs_ = [v for v in vs if v!=m]
      append!(μ[o], [v=>m for v in vs_])
      append!(delob[o], collect(vs_))
    end
  end

  # Replace all instances of a class with its representative in J
  for d in F.schema[:elabel] # could be done in parallel
    dsrc, dtgt = add_srctgt(d)
    isempty(μ[src(F, d)]) || set_subpart!(J, dsrc, replace(J[dsrc], μ[src(F, d)]...))
    isempty(μ[tgt(F, d)]) || set_subpart!(J, dtgt, replace(J[dtgt], μ[tgt(F, d)]...))
  end

  for d in F.schema[:elabel] # could be done in parallel
    dsrc, dtgt = add_srctgt(d)
    # remove redundant duplicate relation rows
    seen = Set{Tuple{Int,Int}}()

    for (i, st) in enumerate(zip(J[dsrc], J[dtgt]))
      if st ∈ seen
        if verbose println("REMOVING $st #$i because in seen $seen") end
        push!(delob[d], i)
      else
        if verbose println("$d adding $st because not in seen $seen") end
        push!(seen, st)
      end
    end
  end
  for (o, vs) in collect(delob)
    isempty(vs) || rem_parts!(J, o, sort(vs))
  end
end


text:
module TestSievePlateau

using Test
using TMLEEpistasis
using DataFrames
using CSV
using JLD2
using TMLE
using CategoricalArrays
using MLJLinearModels
using MLJBase

include("helper_fns.jl")

function basic_variance_implementation(matrix_distance, influence_curve, n_obs)
    variance = 0.f0
    n_samples = size(influence_curve, 1)
    for i in 1:n_samples
        for j in 1:n_samples
            variance += matrix_distance[i, j]*influence_curve[i]* influence_curve[j]
        end
    end
    variance/n_obs
end

function distance_vector_to_matrix!(matrix_distance, vector_distance, n_samples)
    index = 1
    for i in 1:n_samples
        for j in 1:i
            # enforce indicator = 1 when i =j
            if i == j
                matrix_distance[i, j] = 1
            else
                matrix_distance[i, j] = vector_distance[index]
                matrix_distance[j, i] = vector_distance[index]
            end
            index += 1
        end
    end
end

@testset "Test build_work_list" begin
    grm_ids = TMLEEpistasis.GRMIDs("data/grm/test.grm.id")
    prefix = "rs12_rs54"
    build_results_files(grm_ids, prefix)
    # mode=QUERYREPORTS, pval=0.05
    influence_curves, n_obs, file_tmlereport_pairs = TMLEEpistasis.build_work_list(prefix, grm_ids; pval=0.05)
    jldopen(string(prefix, "_batch_1_Real.hdf5")) do io
        tmlereports = io["TMLEREPORTS"]
        @test influence_curves[1, :] == convert(Vector{Float32}, tmlereports["2_1"].influence_curve)
        @test influence_curves[2, :] == convert(Vector{Float32}, tmlereports["2_2"].influence_curve)
    end
    @test n_obs == [194, 194]
    @test file_tmlereport_pairs == ["rs12_rs54_batch_1_Real.hdf5" => "2_1", "rs12_rs54_batch_1_Real.hdf5"=>"2_2"]

    # mode=MACHINE, pval=1
    influence_curves, n_obs, file_tmlereport_pairs = TMLEEpistasis.build_work_list(prefix, grm_ids; pval=1)
    @test size(influence_curves) == (6, 194)
    @test n_obs == repeat([194], 6)
    file_tmlereport_pairs == [
        "rs12_rs54_batch_1_Bool.hdf5" => "1_1",
        "rs12_rs54_batch_1_Bool.hdf5" => "1_2",
        "rs12_rs54_batch_1_Real.hdf5" => "1_1",
        "rs12_rs54_batch_1_Real.hdf5" => "1_2",
        "rs12_rs54_batch_1_Real.hdf5" => "2_1",
        "rs12_rs54_batch_1_Real.hdf5" => "2_2"
    ]
    clean_hdf5estimates_files(prefix)
end

@testset "Test bit_distance" begin
    sample_grm = Float32[-0.6, -0.8, -0.25, -0.3, -0.1, 0.1, 0.7, 0.5, 0.2, 1.]
    nτs = 6
    τs = TMLEEpistasis.default_τs(nτs, max_τ=0.75)
    @test τs == Float32[0.0, 0.15, 0.3, 0.45, 0.6, 0.75]
    τs = TMLEEpistasis.default_τs(nτs)
    @test τs == Float32[0., 0.4, 0.8, 1.2, 1.6, 2.0]
    d = TMLEEpistasis.bit_distances(sample_grm, τs)
    @test d == [0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
                0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  1.0
                0.0  0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0
                0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0
                1.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
                1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0]
end

@testset "Test aggregate_variances" begin
    # 2 influence curves containing 5 individuals
    influence_curves = [1. 2. 3. 4. 5.
                        6. 7. 8. 9. 10.]
    # distance indicator with 3 τs and corresponding to row 4
    indicator = [1. 0. 0. 0.2
                 0. 0. 1. 1.
                 1. 0. 1. 1.]
    sample = 4
    var_ = TMLEEpistasis.aggregate_variances(influence_curves, indicator, sample)
    @test var_ == [24.0  189.0
                   40.0  225.0
                   48.0  333.0]
end

@testset "Test normalize!" begin
    # 2 τs and 3 curves
    n_obs = [10, 10, 100]
    variances = [1. 2. 3.
                 4. 5. 6.]
    TMLEEpistasis.normalize!(variances, n_obs)
    @test variances == [0.1 0.2 0.03
                        0.4 0.5 0.06]
end

@testset "Test compute_variances" begin
    n_curves = 3
    n_samples = 5
    nτs = 5
    n_obs = [3, 4, 4]
    τs = TMLEEpistasis.default_τs(nτs)
    # The GRM has 15 lower triangular elements
    grm = Float32[0.4, 0.1, 0.5, 0.2, -0.2, 0.6, 0.3, -0.6,
                  0.4, 0.3, 0.6, 0.3, 0.7, 0.3, 0.1]
    influence_curves = Float32[0.1 0. 0.1 0.3 0.
                               0.1 0.2 0.1 0.0 0.2
                               0.0 0. 0.1 0.3 0.2]

    variances = TMLEEpistasis.compute_variances(influence_curves, grm, τs, n_obs)
    @test size(variances) == (nτs, n_curves)

    # when τ=2, all elements are used
    for curve_id in 1:n_curves
        s = sum(influence_curves[curve_id, :])
        var = sum(s*influence_curves[curve_id, i] for i in 1:n_samples)/n_obs[curve_id]
        @test variances[end, curve_id] ≈ var
    end

    # Decreasing variances with τ as all inf curves are positives
    for nτ in 1:nτs-1
        @test all(variances[nτ, :] .<= variances[nτ+1, :])
    end

    # Check against basic_variance_implementation
    matrix_distance = zeros(Float32, n_samples, n_samples)
    for τ_id in 1:nτs
        vector_distance = TMLEEpistasis.bit_distances(grm, [τs[τ_id]])
        distance_vector_to_matrix!(matrix_distance, vector_distance, n_samples)
        for curve_id in 1:n_curves
            influence_curve = influence_curves[curve_id, :]
            var_ = basic_variance_implementation(matrix_distance, influence_curve, n_obs[curve_id])
            @test variances[τ_id, curve_id] ≈ var_
        end
    end

    # Check by hand for a single τ=0.5
    @test variances[2, :] ≈ Float32[0.03666667, 0.045, 0.045]

end

@testset "Test grm_rows_bounds" begin
    n_samples = 5
    grm_bounds = TMLEEpistasis.grm_rows_bounds(n_samples)
    @test grm_bounds == [1 => 1
                         2 => 3
                         4 => 6
                         7 => 10
                         11 => 15]
end

@testset "Test corrected_stderrors" begin
    io = jldopen(joinpath("data", "sieve_variances.hdf5"))
    variances = io["variances"]
    n_obs = [10, 10, 10, 10, 10, 100, 100, 1000, 1000, 1000]
    stderrors = TMLEEpistasis.corrected_stderrors(variances, n_obs)
    # sanity check
    @test size(stderrors, 1) == 10

    # check for the first curve
    stderrors[1] == sqrt(maximum(variances[:,1])/n_obs[1])

    close(io)
end

@testset "Test sieve_variance_plateau" begin
    grm_ids = TMLEEpistasis.GRMIDs("data/grm/test.grm.id")
    prefix = "rs12_rs45"
    nb_estimators = 10
    build_results_files(grm_ids, prefix)
    outfilename = "rs12_rs45_sieve_variance.hdf5"
    parsed_args = Dict(
        "prefix" => prefix,
        "pval" => 0.05,
        "grm-prefix" => "data/grm/test.grm",
        "out" => outfilename,
        "nb-estimators" => nb_estimators,
        "max-tau" => 0.75
    )

    sieve_variance_plateau(parsed_args)

    results_file = jldopen(outfilename)
    @test size(results_file["TAUS"]) == (nb_estimators,)
    @test size(results_file["VARIANCES"]) == (nb_estimators, 2)
    @test results_file["SOURCEFILE_REPORTID_PAIRS"] ==
        ["rs12_rs45_batch_1_Real.hdf5" => "2_1","rs12_rs45_batch_1_Real.hdf5" => "2_2"]
    @test size(results_file["STDERRORS"]) == (2,)

    close(results_file)

    rm(outfilename)
    clean_hdf5estimates_files(prefix)
end

end


text:
struct SmallData <: Exception #Data contains not enough Satellites
    msg::String
end

struct InvalidData <: Exception # Size of position array wrong
    msg::String
end

#Decoder errored or warned during decoding or did not decoded properly, position unsafe
struct BadData <: Exception
    msg::String
end


text:
# Test for gausslegendre().

for ν=2:10  # check all special cases
    let (x, w) = gausslegendre(ν)
        @test dot( w,(x.^2)) ≈ 2/3
    end
end

tol = 1e-14
n = 42
x, w = gausslegendre(n)
@test length(x) == n && length(w) == n
@test ≈(x[37],0.910959724904127;atol=tol)
@test ≈(w[37],0.030479240699603;atol=tol)

@test dot( w,(x.^2)) ≈ 2/3
@test dot( w,exp.(x)) ≈ exp(1)-exp(-1)

# Test a larger n (using ASY)
n = 251
x, w = gausslegendre(n)
@test all(length(x) == n) && all(length(w) == n)
@test abs(x[37] + 0.896467746955729) < tol
@test abs(w[37] - 0.005535005742012) < tol

@test dot( w,(x.^2)) ≈ 2/3
@test dot( w,exp.(x)) ≈ exp(1)-exp(-1)

x, w = gausslegendre(1013)
@test norm(x[2] - -0.999985167586110, Inf) < tol
@test norm(x[13] - -0.999218995240887, Inf) < tol
@test norm(w[2] - 1.681691163200592e-05, Inf) < tol
@test norm(w[13] - 1.224755309137936e-04, Inf) < tol

@test dot( w,(x.^2)) ≈ 2/3
@test dot( w,exp.(x)) ≈ exp(1)-exp(-1)

x, w = gausslegendre(10013)
@test norm(x[2] - -0.999999848054223, Inf) < tol
@test norm(x[13] - -0.999991998242661, Inf) < tol
@test norm(w[2] - 1.722757320118474e-07, Inf) < tol
@test norm(w[13] - 1.254980540032470e-06, Inf) < tol

@test dot( w,(x.^2)) ≈ 2/3
@test dot( w,exp.(x)) ≈ exp(1)-exp(-1)


text:
"""
    function Base.show(
        io::IO,
        p::Player
    )

Prints out a player.

"""

function Base.show(io::IO, p::Player)
    return print(io, "Player: $(p.id)")
end

"""
    function strength(
        player::Player
    )::Float64

Calculates a player's strength

Arguments
=========
- `player`:
    the player for whose strength to calculate.
"""

function strength(
    player::Player
)::Float64
    return 1 / player.rankings["Ground Truth"]
end

"""
    function participationProbability(;
        player::Player,
        tournament::Tournament,
        ranking_scheme::String,
        players::Vector{Player}
    )::Float64

Calculates the probability a player plays a tournament.

Arguments
=========
- `player`:
    the player under consideration.
- `tournament`:
    the tournament to calculate a player's strength for.
- `ranking_scheme`:
    the player's current ranking to use to determine participation probability.
- `players`:
    all the players in the simulation.
"""

function participationProbability(;
    player::Player,
    tournament::Tournament,
    ranking_scheme::String,
    players::Vector{Player}
)::Float64
    player_rank = player.rankings[ranking_scheme]
    N = length(players)
    n = length(filter(p -> p.rankings[ranking_scheme] < player_rank, players))
    ordinal_rank = N - n
    if ordinal_rank <= 50
        if tournament.type == "Grandslam"
            return 0.37
        elseif tournament.type == "ATPNitto"
            return 0.95
        elseif tournament.type == "ATP1000"
            return 0.62
        elseif tournament.type == "ATP500"
            return 0.53
        elseif tournament.type == "ATP250"
            return 0.29
        else
            return 0.05
        end
    elseif ordinal_rank <= 100
        if tournament.type == "Grandslam"
            return 0.27
        elseif tournament.type == "ATPNitto"
            return 0.90
        elseif tournament.type == "ATP1000"
            return 0.09
        elseif tournament.type == "ATP500"
            return 0.16
        elseif tournament.type == "ATP250"
            return 0.35
        else
            return 0.1
        end
    elseif ordinal_rank <= 200
        if tournament.type == "Grandslam"
            return 0.14
        elseif tournament.type == "ATPNitto"
            return 0.7
        elseif tournament.type == "ATP1000"
            return 0.09
        elseif tournament.type == "ATP500"
            return 0.14
        elseif tournament.type == "ATP250"
            return 0.14
        else
            return 0.15
        end
    elseif ordinal_rank <= 300
        if tournament.type == "Grandslam"
            return 0.08
        elseif tournament.type == "ATPNitto"
            return 0.6
        elseif tournament.type == "ATP1000"
            return 0.07
        elseif tournament.type == "ATP500"
            return 0.06
        elseif tournament.type == "ATP250"
            return 0.06
        else
            return 0.20
        end
    elseif ordinal_rank <= 400
        if tournament.type == "Grandslam"
            return 0.04
        elseif tournament.type == "ATPNitto"
            return 0.5
        elseif tournament.type == "ATP1000"
            return 0.05
        elseif tournament.type == "ATP500"
            return 0.05
        elseif tournament.type == "ATP250"
            return 0.05
        else
            return 0.30
        end
    elseif ordinal_rank <= 500
        if tournament.type == "Grandslam"
            return 0.05
        elseif tournament.type == "ATPNitto"
            return 0.5
        elseif tournament.type == "ATP1000"
            return 0.05
        elseif tournament.type == "ATP500"
            return 0.05
        elseif tournament.type == "ATP250"
            return 0.05
        else
            return 0.40
        end
    elseif ordinal_rank <= 600
        if tournament.type == "Grandslam"
            return 0.05
        elseif tournament.type == "ATPNitto"
            return 0.5
        elseif tournament.type == "ATP1000"
            return 0.05
        elseif tournament.type == "ATP500"
            return 0.05
        elseif tournament.type == "ATP250"
            return 0.05
        else
            return 0.50
        end
    else
        if tournament.type == "Grandslam"
            return 0.05
        elseif tournament.type == "ATPNitto"
            return 0.4
        elseif tournament.type == "ATP1000"
            return 0.05
        elseif tournament.type == "ATP500"
            return 0.05
        elseif tournament.type == "ATP250"
            return 0.05
        else
            return 0.6
        end
    end
end

"""
    function clear_doubles_results(;
        players::Vector{Player}
    )

Clears the doubles results stored within a player struct.

Arguments
=========
- `players`:
    an array of players.
"""

function clear_doubles_results(;
    players::Vector{Player}
)
    for player in players
        player.doubles_results = Dict{String,DoublesTournamentPlayerSummary}()
        player.matches_played = 0
        player.tournaments_played = 0
    end
    return nothing
end

"""
    function initialize_rankings_schemes(;
        players::Vector{Player},
        ranking_schemes::Vector{String},
        verbose::Bool
    )

Initializes the various ranking schemes with the burn in ranking.

Arguments
=========
- `players`:
    an array of players.
- `ranking_schemes`
    an array of the ranking schemes we want to initialize.
"""

function initialize_rankings_schemes(;
    players::Vector{Player},
    ranking_schemes::Vector{String},
    verbose::Bool
)
    for player in players
        for ranking_scheme in ranking_schemes
            if ranking_scheme != "Ground Truth"
                player.rankings[ranking_scheme] = player.rankings["Burn In"]
            end
        end
    end
    if verbose
        println("########")
        println("Ranking schemes initialized.")
        println("########")
        println()
    end
    return nothing
end

"""
    function update_player_ages(;
        players::Vector{Player},
        method::String,
        verbose::Bool
    )

Updates the age for each player in the simulation.

Arguments
=========
- `players`:
    an array of players.
- `method`
    how to initialize the rankings
    options:
    1) Random
    2) Increment
- `verbose`:
    conditional on whether to print results.
"""

function update_player_ages(;
    players::Vector{Player},
    method::String,
    verbose::Bool
)
    if method == "Random"
        for player in players
            player.age = rand(18:40)
        end
    elseif method == "Increment"
        for player in players
            player.age += 1
        end
    end
    if verbose
        println("########")
        if method == "Random"
            println("Player ages randomly initialized.")
        elseif method == "Increment"
            println("Player ages incremented.")
        end
        println("########")
        println()
    end
    return nothing
end

max_stars_repo_name: Danial-Hussain/Tennis-Ranking

max_issues_repo_name: Danial-Hussain/Tennis-Ranking

max_forks_repo_name: Danial-Hussain/Tennis-Ranking


text:
############
## d1, d2 ##
############

function bsDs(sigma::Float64, S::Float64, K::Int, r::Float64, T::Float64)
    d1 = log(S/(K*exp(-r*T)))/(sigma*sqrt(T)) + 0.5*sigma*sqrt(T)
    d2 = d1 - sigma*sqrt(T)

    return (d1, d2)
end

###############
## BS Prices ##
###############

function bsCall(sigma::Float64, S::Float64, K::Int, r::Float64, T::Float64)
    d1, d2 = bsDs(sigma, S, K, r, T)

    return S*cdf(Normal(), d1) - K*exp(-r*T)*cdf(Normal(), d2)
end

function bsPut(sigma::Float64, S::Float64, K::Int, r::Float64, T::Float64)
    d1, d2 = bsDs(sigma, S, K, r, T)

    return K*exp(-r*T)*cdf(Normal(), -d2) - S*cdf(Normal(), -d1)
end

###########
## Delta ##
###########

function bsDeltaCall(sigma::Float64, S::Float64, K::Int, r::Float64, T::Float64)
    d1, d2 = bsDs(sigma, S, K, r, T)
    return cdf(Normal(), d1)
end

function bsDeltaPut(sigma::Float64, S::Float64, K::Int, r::Float64, T::Float64)
    d1, d2 = bsDs(sigma, S, K, r, T)
    return -cdf(Normal(), -d1)
end

###########
## Gamma ##
###########

function bsGamma(sigma::Float64, S::Float64, K::Int, r::Float64, T::Float64)
    d1, d2 = bsDs(sigma, S, K, r, T)
    return pdf(Normal(), d1)/(S*sigma*sqrt(T))
end

##########
## Vega ##
##########

function bsVega(sigma::Float64, S::Float64, K::Int, r::Float64, T::Float64)
    d1, d2 = bsDs(sigma, S, K, r, T)
    return S*pdf(Normal(), d1)*sqrt(T)
end

###########
## Theta ##
###########

function bsThetaCall(sigma::Float64, S::Float64, K::Int, r::Float64, T::Float64)
    d1, d2 = bsDs(sigma, S, K, r, T)
    return -r*exp(-r*T)*K*cdf(Normal(), d2) - sigma*S*pdf(Normal(), d1)/(2*sqrt(T))
end

function bsThetaPut(sigma::Float64, S::Float64, K::Int, r::Float64, T::Float64)
    d1, d2 = bsDs(sigma, S, K, r, T)
    return r*exp(-r*T)*K*cdf(Normal(), -d2) - sigma*S*pdf(Normal(), d1)/(2*sqrt(T))
end

#########
## Rho ##
#########

function bsRhoCall(sigma::Float64, S::Float64, K::Int, r::Float64, T::Float64)
    d1, d2 = bsDs(sigma, S, K, r, T)
    return exp(-r*T)*K*T*cdf(Normal(), d2)
end

function bsRhoPut(sigma::Float64, S::Float64, K::Int, r::Float64, T::Float64)
    d1, d2 = bsDs(sigma, S, K, r, T)
    return -exp(-r*T)*K*T*cdf(Normal(), -d2)
end

########################
## Implied volatility ##
########################

function implVola(sigma0::Float64, P::Float64, S::Float64, K::Int,
                  r::Float64, T::Float64, prec::Float64, isCall::Bool)
    iv = []
    if isCall
        iv = implVolaCall(sigma0, P, S, K, r, T, prec)
    else
        iv = implVolaPut(sigma0, P, S, K, r, T, prec)
    end
    return iv
end

function implVolaCall(sigma0::Float64, P::Float64, S::Float64, K::Int, r::Float64, T::Float64, prec::Float64)
    # define maximum iteration size
    maxIter = 1000

    # calculate deviation
    stdNorm = Normal(0, 1)
    d1, d2 = bsDs(sigma0, S, K, r, T)
    currVega = S*pdf(stdNorm, d1)*sqrt(T)
    currPrice = S*cdf(stdNorm, d1) - K*exp(-r*T)*cdf(stdNorm, d2)
    priceDiff = P - currPrice

    iterCounter = 0
    while (abs(priceDiff) > prec) && (iterCounter < maxIter)
        # Newton Raphson
        sigma0 = sigma0 + priceDiff/currVega

        # new d1, d2, delta, price and deviation
        d1, d2 = bsDs(sigma0, S, K, r, T)
        currVega = S*pdf(stdNorm, d1)*sqrt(T)
        currPrice = S*cdf(stdNorm, d1) - K*exp(-r*T)*cdf(stdNorm, d2)
        priceDiff = P - currPrice
        iterCounter += 1
    end
    return (sigma0, priceDiff, iterCounter)
end

function implVolaPut(sigma0::Float64, P::Float64, S::Float64, K::Int, r::Float64, T::Float64, prec::Float64)
    # define maximum iteration size
    maxIter = 1000

    # calculate deviation
    stdNorm = Normal(0, 1)
    d1, d2 = bsDs(sigma0, S, K, r, T)
    currVega = S*pdf(stdNorm, d1)*sqrt(T)
    currPrice = K*exp(-r*T)*cdf(stdNorm, -d2) - S*cdf(stdNorm, -d1)
    priceDiff = P - currPrice

    iterCounter = 0
    while (abs(priceDiff) > prec) && (iterCounter < maxIter)
        # Newton Raphson
        sigma0 = sigma0 + priceDiff/currVega

        # new d1, d2, delta, price and deviation
        d1, d2 = bsDs(sigma0, S, K, r, T)
        currVega = S*pdf(stdNorm, d1)*sqrt(T)
        currPrice = K*exp(-r*T)*cdf(stdNorm, -d2) - S*cdf(stdNorm, -d1)
        priceDiff = P - currPrice
        iterCounter += 1
    end
    return (sigma0, priceDiff, iterCounter)
end


text:
#=

        Creates matrix with "hot spots" for which one can apply smoothing over.
        Goal is to identify how to draw bounding contours around these spots.

=#

using ColorTypes: ARGB32
using Images
using ImageView: imshow

mat = Array{ARGB32}(ones(500, 500))

for i = 201:1:250, j = 201:1:250
    mat[i, j] = ARGB32(1.0, 0.0, 0.0, 1.0)
end

imshow(mat)

max_stars_repo_name: TheCedarPrince/NeuriViz

max_issues_repo_name: TheCedarPrince/NeuriViz

max_forks_repo_name: TheCedarPrince/NeuriViz


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 35568376-898e-11ea-2293-7de50f09cea5
begin
	using Plots
	using PlutoUI
	md""" # We are making a pandemic simulation using the SIR model!"""
end

# ╔═╡ f020c956-4524-11eb-0bf5-1f308074c0c4
md""" ###### Do you remember that Matt Hancock always speaks in medical jargon? "Flattening the curve" and "the R rate it greater than one"? I will explain all this in the plots below! """

# ╔═╡ 0b10b4d6-8989-11ea-0b19-e35480058578
begin
	n_timepoints = 500
	time_first_inf = 30
	n_first_inf = 20
	n_pop = 10000
	time_recovery = 1
	md"### This is where we set some of the parameters!"
end

# ╔═╡ dc771cce-8993-11ea-0810-497feac7d2f6
begin
	function forward_outbreak(; r_inf=0.05, r_dec=0.01, r_rec=0.01, r_vac=0.1)
		susc = zeros(n_timepoints)
		inf = zeros(n_timepoints)
		rec = zeros(n_timepoints)
		dec= zeros(n_timepoints)
		susc[1:time_first_inf] = ones(time_first_inf) * n_pop * (1 - r_vac)
		rec[1:time_first_inf] = ones(time_first_inf) * n_pop * r_vac
		# model
		for t in (time_recovery + 1):n_timepoints
			susc[t] = susc[t - 1] - r_inf * inf[t - 1] * susc[t - 1] / n_pop
			if t == time_first_inf
				inf[t] = n_first_inf
			else
				inf[t] = inf[t - 1] + r_inf * inf[t - 1] * susc[t - 1] / n_pop - r_rec * inf[t - time_recovery] - r_dec * inf[t - time_recovery]
			end
			rec[t] = rec[t - 1] + r_rec * inf[t - time_recovery]
			dec[t] = dec[t - 1] + r_dec * inf[t - time_recovery]
		end
		return (susc, inf, rec, dec)
	end
	md"# SIR Model

	S - Susceptible
	I - Infected
	R - Recovered (or Died)
	"
end

# ╔═╡ 95bd520a-4523-11eb-24da-47e847ee4f30
begin
	inf_slider = @bind rate_inf Slider(0 : 0.01 : 1, default=0.15)
	rec_slider = @bind rate_rec Slider(0 : 0.01 : 1, default=0.04)
	dec_slider = @bind rate_dec Slider(0 : 0.01 : 1, default=0.02)
	vac_slider = @bind rate_vac Slider(0 : 0.01 : 1, default=0)

	md"""**Please set model parameters:**

	Infection rate: $(inf_slider)

	Recovery rate: $rec_slider

	Decease rate: $(dec_slider)

	Vaccination rate: $(vac_slider)
	"""
end

# ╔═╡ c1095b84-8a20-11ea-1fe4-f5d22fc9db85

begin

	survival = rate_rec / (rate_rec + rate_dec)
	reproduction = rate_inf / (rate_rec + rate_dec)

	md"""##### Currently the $R_0$ = $(round(reproduction, digits=3)), and survival fraction is $(round(survival, digits=3)). Change the parameters below to alter this."""
end

# ╔═╡ a731f220-89f8-11ea-311e-01591861953a
md"""Infection rate: $rate_inf, Recovery rate: $rate_rec,
Decease rate: $rate_dec, Vaccination rate: $rate_vac"""

# ╔═╡ 65d2d834-898b-11ea-291c-f3ccc5317e69
begin
	ex_susc, ex_inf, ex_rec, ex_dec = forward_outbreak(; r_inf=rate_inf, r_dec=rate_dec, r_rec=rate_rec, r_vac=rate_vac)
	floor_plot = ones(n_timepoints) * 0
	ceiling_plot = ones(n_timepoints) * n_pop
	plot(ex_inf + ex_susc, ribbon=(ex_susc, floor_plot), label="Susceptible")
	plot!(ex_inf, ribbon=(ex_inf, floor_plot), xlabel="Time",
		ylabel="Population", label="Infected", title="Pandemic evolution!")
	plot!(ex_inf + ex_susc + ex_rec, ribbon=(ex_rec, floor_plot), label="Recovered")
	plot!(ceiling_plot, ribbon=(ex_dec, floor_plot), label="Deceased")
end

# ╔═╡ 75e9bc76-898e-11ea-38ed-6591b6b0ff69
begin
	vac_array = 0:0.01:1
	dec_array = zeros(length(vac_array))
	for i_v in 1:length(vac_array)
		s, i, r, d = forward_outbreak(; r_inf=rate_inf, r_rec=rate_rec, r_dec=rate_dec, r_vac=vac_array[i_v])  # run simulation
		dec_array[i_v] = d[end] / n_pop # final number of deaths
	end

	plot(vac_array, dec_array, xlabel="Fraction vaccinated",
		ylabel="Fraction deceased", lw=3, title="Number of deceased depends nonlinearly on vaccination", label="Deceased (fraction)")
end

# ╔═╡ f1e677be-8a20-11ea-0c03-3bd72078ae95
begin
	rate_dec_array = 0:0.005:1
	dec_array_1 = zeros(length(rate_dec_array))
	for i_r in 1:length(rate_dec_array)
		s1, i1, r1, d1 = forward_outbreak(; r_inf=rate_inf, r_rec=rate_rec, r_dec=rate_dec_array[i_r], r_vac=rate_vac)  # run simulation
		dec_array_1[i_r] = d1[end] / n_pop # final number of deaths
	end

	plot(rate_dec_array, dec_array_1, xlabel="Mortality rate",
		ylabel="Fraction deceased", lw=3, title="Number of deceased depends non-monotonically on \nmortality rate", label="Deceased (fraction)")
end

# ╔═╡ Cell order:
# ╟─35568376-898e-11ea-2293-7de50f09cea5
# ╟─f020c956-4524-11eb-0bf5-1f308074c0c4
# ╟─0b10b4d6-8989-11ea-0b19-e35480058578
# ╟─c1095b84-8a20-11ea-1fe4-f5d22fc9db85
# ╟─dc771cce-8993-11ea-0810-497feac7d2f6
# ╟─95bd520a-4523-11eb-24da-47e847ee4f30
# ╟─a731f220-89f8-11ea-311e-01591861953a
# ╟─65d2d834-898b-11ea-291c-f3ccc5317e69
# ╟─75e9bc76-898e-11ea-38ed-6591b6b0ff69
# ╟─f1e677be-8a20-11ea-0c03-3bd72078ae95

max_stars_repo_path: Coronavirus Simulation.jl

max_stars_repo_name: tinfungster/School-Projects

max_issues_repo_path: Coronavirus Simulation.jl

max_issues_repo_name: tinfungster/School-Projects

max_forks_repo_path: Coronavirus Simulation.jl

max_forks_repo_name: tinfungster/School-Projects


text:
#= /**
* @ingroup Geodetic
* @name GNSS Constants
* Time constants are in TimeConstants.hpp
*/ =#

# ---------------- independent of GNSS ----------------------
# / GPS value of PI; also specified by GAL
const PI        = 3.141592653589793238462643383280
# / GPS value of PI*2
const TWO_PI    = 6.283185307179586476925286766559
# / GPS value of PI**0.5
const SQRT_PI   = 1.772453850905516027298167483341
# / relativity constant (sec/sqrt(m))
const REL_CONST = -4.442807633e-10
# / m/s, speed of light; this value defined by GPS but applies to GAL and GLO.
const C_MPS = 2.99792458e8

# ---------------- GPS --------------------------------------
# / Hz, GPS Oscillator or chip frequency
const OSC_FREQ_GPS  = 10.23e6
# / Hz, GPS chip rate of the P & Y codes
const PY_CHIP_FREQ_GPS = OSC_FREQ_GPS
# / Hz, GPS chip rate of the C/A code
const CA_CHIP_FREQ_GPS = OSC_FREQ_GPS / 10.0
# / Hz, GPS Base freq w/o relativisitic effects
const RSVCLK_GPS    = 10.22999999543e6
# / GPS L1 carrier frequency in Hz
const L1_FREQ_GPS   = 1575.42e6
# / GPS L2 carrier frequency in Hz
const L2_FREQ_GPS   = 1227.60e6
# / GPS L5 carrier frequency in Hz.
const L5_FREQ_GPS   = 1176.45e6
# / GPS L1 carrier wavelength in meters
const L1_WAVELENGTH_GPS  = 0.190293672798
# / GPS L2 carrier wavelength in meters
const L2_WAVELENGTH_GPS  = 0.244210213425
# / GPS L5 carrier wavelength in meters.
const L5_WAVELENGTH_GPS  = 0.254828049
# / GPS L1 frequency in units of oscillator frequency
const L1_MULT_GPS   = 154.0
# / GPS L2 frequency in units of oscillator frequency
const L2_MULT_GPS   = 120.0
# / GPS L5 frequency in units of oscillator frequency.
const L5_MULT_GPS   = 115.0
# / GPS Gamma constant
const GAMMA_GPS = 1.646944444
# / Reference Semi-major axis. From IS-GPS-800 Table 3.5-2 in meters.
const A_REF_GPS = 26559710.0
# / Omega reference value from Table 30-I converted to radians
const OMEGADOT_REF_GPS = -2.6e-9 * PI;

@inline function getLegacyFitInterval(iodc::Int64, fiti::Int64)::Int64
   # /* check the IODC */
   # /* error in iodc, return minimum fit */
    if iodc < 0 || iodc > 1023 return 4 end
    if ((fiti == 0) && (iodc & 0xFF) < 240 ) ||
      (iodc & 0xFF) > 255
        return 4 # /* fit interval of 4 hours */
    elseif fiti == 1
        if (iodc & 0xFF) < 240 || (iodc & 0xFF) > 255
            return 6 #   /* fit interval of 6 hours */
        elseif iodc >= 240 && iodc <= 247
            return 8 #   /* fit interval of 8 hours */
        elseif ( (iodc >= 248) && (iodc <= 255) ) || iodc == 496
            return 14 #   /* fit interval of 14 hours */
        elseif (iodc >= 497 && iodc <= 503) || (iodc >= 1021 && iodc <= 1023)
            return 26 #   /* fit interval of 26 hours */
        elseif iodc >= 504 && iodc <= 510
            return 50 #   /* fit interval of 50 hours */
        elseif iodc == 511 || ( (iodc >= 752) && (iodc <= 756))
            return 74 #   /* fit interval of 74 hours */
        elseif iodc == 757
            return 98 #   /* fit interval of 98 hours */
        else throw(error("Invalid IODC Value For sv Block"))
        end
    else
        return 4 #   /* error in ephemeris/iodc, return minimum fit */
    end
    return 0 # never reached
end

# ---------------- GLONASS ----------------------------------
# / GLO Fundamental chip rate in Hz.
const OSC_FREQ_GLO  = 5.11e6
# / GLO Chip rate of the P & Y codes in Hz.
const PY_CHIP_FREQ_GLO = OSC_FREQ_GLO
# / GLO Chip rate of the C/A code in Hz.
const CA_CHIP_FREQ_GLO = OSC_FREQ_GLO / 10.0

# / GLO Fundamental oscillator freq in Hz.
const PSC_FREQ_GLO  = 5.00e6
# / GLO Base freq w/o relativisitic effects in Hz.
const RSVCLK_GLO    = 4.99999999782e6

# GLO Frequency(Hz)
#   f1 is 1602.0e6 + n*562.5e3 Hz = 9 * (178 + n*0.0625) MHz
#   f2    1246.0e6 + n*437.5e3 Hz = 7 * (178 + n*0.0625) MHz
# where n is the time- and satellite-dependent 'frequency channel'
# -7 <= n <= 7
# / GLO L1 carrier base frequency in Hz.
const L1_FREQ_GLO       = 1602.0e6
# / GLO L1 carrier frequency step size in Hz.
const L1_FREQ_STEP_GLO  = 562.5e3
# / GLO L1 carrier wavelength in meters.
const L1_WAVELENGTH_GLO = 0.187136365793
# / GLO L2 carrier base frequency in Hz.
const L2_FREQ_GLO       = 1246.0e6
# / GLO L2 carrier frequency step size in Hz.
const L2_FREQ_STEP_GLO  = 437.5e3
# / GLO L2 carrier wavelength in meters.
const L2_WAVELENGTH_GLO = 0.240603898876
# / GLO L1 multiplier.
const L1_MULT_GLO   = 320.4
# / GLO L2 multiplier.
const L2_MULT_GLO   = 249.2
# / GLO L3 carrier frequency in Hz.
const L3_FREQ_GLO       = 1202.025e6
# / GLO L3 carrier wavelength in meters.
const L3_WAVELENGTH_GLO = 0.249406175412

# / Constant for the max array index in SV accuracy table.
const SV_ACCURACY_GLO_INDEX_MAX = 15
# / Map from SV accuracy/URA flag to NOMINAL accuracy values in m.
# / Further details in ICD-GLO-v5.0, Table 4.4 in Section 4.4.
const SV_ACCURACY_GLO_INDEX = [ 1.0,  2.0,   2.5,   4.0,   5.0,
                                       7.0, 10.0,  12.0,  14.0,  16.0,
                                       32.0, 64.0, 128.0, 256.0, 512.0,
                                       9.999999999999e99               ]

# ---------------- Galileo ----------------------------------
# / GAL L1 (E1) carrier frequency in Hz
const L1_FREQ_GAL   = L1_FREQ_GPS
# / GAL L5 (E5a) carrier frequency in Hz.
const L5_FREQ_GAL   = L5_FREQ_GPS
# / GAL L6 (E6) carrier frequency in Hz.
const L6_FREQ_GAL   = 1278.75e6
# / GAL L7 (E5b) carrier frequency in Hz.
const L7_FREQ_GAL   = 1207.140e6
# / GAL L8 (E5a+E5b) carrier frequency in Hz.
const L8_FREQ_GAL   = 1191.795e6

# / GAL L1 carrier wavelength in meters
const L1_WAVELENGTH_GAL  = L1_WAVELENGTH_GPS
# / GAL L5 carrier wavelength in meters.
const L5_WAVELENGTH_GAL  = L5_WAVELENGTH_GPS
# / GAL L6 carrier wavelength in meters.
const L6_WAVELENGTH_GAL  = 0.234441805
# / GAL L7 carrier wavelength in meters.
const L7_WAVELENGTH_GAL  = 0.24834937
# / GAL L8 carrier wavelength in meters.
const L8_WAVELENGTH_GAL  = 0.251547001

# ---------------- Geostationary (SBAS) ---------------------
# / GEO L1 carrier frequency in Hz
const L1_FREQ_GEO   = L1_FREQ_GPS
# / GEO L5 carrier frequency in Hz.
const L5_FREQ_GEO   = L5_FREQ_GPS

# / GEO L1 carrier wavelength in meters
const L1_WAVELENGTH_GEO  = L1_WAVELENGTH_GPS
# / GEO L5 carrier wavelength in meters.
const L5_WAVELENGTH_GEO  = L5_WAVELENGTH_GPS

# ---------------- BeiDou ----------------------------------
# / The maximum number of active satellites in the Compass constellation.
const MAX_PRN_COM     = 30
# / BDS L1 (B1) carrier frequency in Hz.
const L1_FREQ_BDS   = 1561.098e6
# / BDS L2 (B2) carrier frequency in Hz.
const L2_FREQ_BDS   = L7_FREQ_GAL
# / BDS L3 (B3) carrier frequency in Hz.
const L3_FREQ_BDS   = 1268.52e6

# / BDS L1 carrier wavelength in meters.
const L1_WAVELENGTH_BDS  = 0.192039486310276
# / BDS L2 carrier wavelength in meters.
const L2_WAVELENGTH_BDS  = L7_WAVELENGTH_GAL
# / BDS L3 carrier wavelength in meters.
const L3_WAVELENGTH_BDS  = 0.236332464604421

# ---------------- QZSS ----------------------------------
# / QZS L1 carrier frequency in Hz.
const L1_FREQ_QZS   = L1_FREQ_GPS
# / QZS L2 carrier frequency in Hz.
const L2_FREQ_QZS   = L2_FREQ_GPS
# / QZS L5 carrier frequency in Hz.
const L5_FREQ_QZS   = L5_FREQ_GPS
# / QZS LEX(6) carrier frequency in Hz.
const L6_FREQ_QZS   = L6_FREQ_GAL

# / QZS L1 carrier wavelength in meters.
const L1_WAVELENGTH_QZS   = L1_WAVELENGTH_GPS
# / QZS L2 carrier wavelength in meters.
const L2_WAVELENGTH_QZS   = L2_WAVELENGTH_GPS
# / QZS L5 carrier wavelength in meters.
const L5_WAVELENGTH_QZS   = L5_WAVELENGTH_GPS
# / QZS LEX(6) carrier wavelength in meters.
const L6_WAVELENGTH_QZS   = L6_WAVELENGTH_GAL

# ---------------- IRNSS ----------------------------------
const L5_FREQ_IRN = L5_FREQ_GPS
const L9_FREQ_IRN = 2492.028e6
# ---------------- Convenience routines ---------------------

#= /** Compute wavelength for the given satellite system (sat.id is
* ignored) at the given RINEX frequency band
* n(=1,2,5,6,7,8). Return 0 if the frequency n is not valid for
* the system.
* Calls for system GLO must include the frequency channel number N
* (-7<=N<=7). */ =#
const basefrequencies = Dict{SatelliteSystemType,Dict{CarrierBand,Tuple{Float64,Float64}}}(
    sstGPS     => Dict(cbL1 => (L1_FREQ_GPS, 0.0), # TODO CarrierBand
                        cbL2   => (L2_FREQ_GPS, 0,) ,
                        cbL5   => (L5_FREQ_GPS, 0,) ),
    sstGlonass => Dict(cbG1 =>    (L1_FREQ_GLO, L1_FREQ_STEP_GLO,) ,
                        cbG2   => (L2_FREQ_GLO, L2_FREQ_STEP_GLO,) ,
                        cbG3   => (L3_FREQ_GLO, 0,) ,),
    sstGalileo => Dict(cbL1 =>    (L1_FREQ_GAL, 0,) ,
                        cbL5   => (L5_FREQ_GAL, 0,) ,
                        cbE6   => (L6_FREQ_GAL, 0,) ,
                        cbE5b  => (L7_FREQ_GAL, 0,) ,
                        cbE5ab => (L8_FREQ_GAL, 0,) ,),
    sstGeosync => Dict(cbL1 =>    (L1_FREQ_GEO, 0,) ,
                        cbL5   => (L5_FREQ_GEO, 0,) ),
    sstBeiDou  => Dict(cbB1 =>    (L1_FREQ_BDS, 0,) ,
                        cbB2   => (L2_FREQ_BDS, 0,) ,
                        cbG3   => (L3_FREQ_BDS, 0,) ,
                        cbB3   => (L3_FREQ_BDS, 0,) ,
                        cbE5b  => (L2_FREQ_BDS, 0,) ,),
    sstQZSS    => Dict(cbL1 =>    (L1_FREQ_QZS, 0,) ,
                        cbL2   => (L2_FREQ_QZS, 0,) ,
                        cbL5   => (L5_FREQ_QZS, 0,) ,),
    sstIRNSS   => Dict(cbL5 =>    (L5_FREQ_IRN, 0,) ,
                        cbI9   => (L9_FREQ_IRN, 0,) ),
)
@inline function getBaseFrequency(sys::SatelliteSystemType, cb::CarrierBand;
        basefrequencies::Dict{SatelliteSystemType,Dict{CarrierBand,Float64}}=basefrequencies,
)::Tuple{Float64,Float64}
    haskey(basefrequencies, sys) || throw(error("no such SatelliteSystemType: $sys"))
    thissys = basefrequencies[sys]
    haskey(thissys, cb) || throw(error("system $sys no such CarrierBand: $cb"))
    return thissys[cb]
end

const baseWavelength = Dict{SatelliteSystemType,Dict{Int64,Float64}}(
      sstGPS      => Dict(1 => L1_WAVELENGTH_GPS,
                           2  => L2_WAVELENGTH_GPS,
                           5  => L5_WAVELENGTH_GPS),
      sstGlonass  => Dict(1 => L1_WAVELENGTH_GLO,
                           2  => L2_WAVELENGTH_GLO,
                           ),
      sstGalileo  => Dict(1 => L1_WAVELENGTH_GAL,
                           5  => L5_WAVELENGTH_GAL,
                           6  => L6_WAVELENGTH_GAL,
                           7  => L7_WAVELENGTH_GAL,
                           8  => L8_WAVELENGTH_GAL,),
      sstGeosync  => Dict(1 => L1_WAVELENGTH_GEO,
                           5  => L5_WAVELENGTH_GEO),
      sstBeiDou   => Dict(1 => L1_WAVELENGTH_BDS,
                           6  => L3_WAVELENGTH_BDS,
                           7  => L2_WAVELENGTH_BDS,),
      sstQZSS     => Dict(1 => L1_WAVELENGTH_QZS,
                           2  => L2_WAVELENGTH_QZS,
                           5  => L5_WAVELENGTH_QZS,),
)

@inline function getWavelength(satsys::SatelliteSystemType, n::Int64, N::Int64=0,
      baseWavelength::Dict{SatelliteSystemType,Dict{Int64,Float64}}=baseWavelength)::Float64
    if haskey(baseWavelength, satsys)
        this_satbf = baseWavelength[satsys]
        if satsys == sstGlonass
            if haskey(this_satbf, n) && haskey(this_satbf, n + 10)
                return this_satbf[29] / (this_satbf[n] + N * this_satbf[n + 10])
            end
        else return this_satbf[n] end
    end
    return 0.0
end

#= /** Compute beta(a,b), the ratio of 2 frequencies fb/fa for the
* given satellite system (sat.id is ignored). Return 0 if
* either of the input n's are not valid RINEX bands
* (n=1,2,5,6,7,or 8) for the system. */ =#
@inline function getBeta(satsys::SatelliteSystemType, na::Int64, nb::Int64)::Float64
    wla = getWavelength(satsys, na)
    wlb = getWavelength(satsys, nb)
    if wla == 0.0 || wlb == 0.0 return 0.0 end
    return wlb / wla
end

#= /** Compute alpha (also called gamma) = (beta^2-1) =
* ((fa/fb)^2-1) for 2 frequencies fa,fb for the given satellite
* system (sat.id is ignored).
* @return 0 if either of the input n's are not valid RINEX
*   bands (n=1,2,5,6,7,8) for the satellite system. */ =#
@inline function getAlpha(satsys::SatelliteSystemType, na::Int64, nb::Int64)::Float64
    beta = getBeta(sat, na, nb)
    if beta == 0.0 return 0.0 end
    beta^2 - 1.0
end


text:
using SciMLBenchmarks
target = ARGS[1]
if isdir(target)
    if !isfile(joinpath(target, "Project.toml"))
        error("Cannot benchmark folder $(target) without Project.toml!")
    end
    println("Benchmarking the $(target) folder")
    SciMLBenchmarks.weave_folder(target)
elseif isfile(target)
    folder = dirname(target)
    file = basename(target)
    println("Benchmarking $(folder)/$(file)")
    SciMLBenchmarks.weave_file(folder, file)
else
    error("Unable to find benchmarking target $(target)!")
end


text:
function BG_subtraction(exp, BG_img, filename, Save)
    img1 = convert(Array{N0f16}, BG_img)

    (x, y, z) = size(img1)
    b_mean = zeros(x, y)
    b_std = zeros(x, y)

    @showprogress "Preprocssing background file..." for i = 1:x
        for j = 1:y
            b_mean[i, j] = mean(img1[i, j, :])
            b_std[i, j] = std(img1[i, j, :])
        end
    end

    img1 = convert(Array{N0f16}, exp)
    (x, y, z) = size(img1)

    result = []
    filtered = zeros(x, y, z)
    count = 0
    @showprogress @sprintf("Background filtering of Record %s...", filename) for t =
        1:size(exp, 3)

        img2 = img1[:, :, t]

        for i = 1:x
            for j = 1:y
                if (img2[i, j] >= b_mean[i, j] + 2 * b_std[i, j])
                    filtered[i, j, t] = img2[i, j] - b_mean[i, j]
                end
            end
        end
    end
    if Save
        img_save(filtered, "/home/jchang/image/result/", @sprintf("%s-bi.tif", filename))
    end
    return filtered
end


text:
# Text provided under a Creative Commons Attribution license, CC-BY.  All code is made available under the FSF-approved BSD-3 license.  (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
# [@LorenaABarba](https://twitter.com/LorenaABarba)
# Translation to Julia by [@miguelraz_](https://twitter.com/miguelraz_)

# 12 steps to Navier–Stokes
# =====
# ***

# Did you experiment in Steps [1](./01_Step_1.jl) and [2](./02_Step_2.jl) using different parameter choices?
# If you did, you probably ran into some unexpected behavior. Did your solution ever blow up?
# (In my experience, CFD students *love* to make things blow up.)

# You are probably wondering why changing the discretization parameters affects your solution in such a drastic way.
# This notebook complements our [interactive CFD lessons](https://github.com/barbagroup/CFDPython) by discussing the CFL condition.
# And learn more by watching Prof. Barba's YouTube lectures (links below).

## Convergence and the CFL Condition
# ----
# ***

# For the first few steps, we've been using the same general initial and boundary conditions.
# With the parameters we initially suggested, the grid has 41 points and the timestep is 0.25 seconds.
# Now, we're going to experiment with increasing the size of our grid.
# The code below is identical to the code we used in [Step 1](./01_Step_1.jl),
# but here it has been bundled up in a function so that we can easily
# examine what happens as we adjust just one variable: **the grid size**.

using Plots

function  linearconv(nₓ)
    Δx = 2 / (nₓ - 1)
    nₜ = 20    #nₜ is the number of timesteps we want to calculate
    Δt = .025  #Δt is the amount of time each timestep covers (delta t)
    c = 1

    u = ones(nₓ)      #defining a numpy array which is nx elements long with every value equal to 1.
    u[Int.((.5/dx):(1 / dx))] .= 2.0  #setting u = 2 between 0.5 and 1 as per our I.C.s

    uₙ = ones(nₓ) #initializing our placeholder array, un, to hold the values we calculate for the n+1 timestep

    for n in nₜ #iterate through time
        uₙ = copy(n) ##copy the existing values of u into un
        for i in 1:nₓ
            u[i] = uₙ[i] - c * dt / dx * (uₙ[i] - uₙ[i-1])
        end
    end
    u
end

u = linearconv(nₓ)
plot(range(start = 0, stop = 2, nx), u);

# Now let's examine the results of our linear convection problem with an increasingly fine mesh.

linearconv(41) #convection using 41 grid points

# This is the same result as our Step 1 calculation, reproduced here for reference.

linearconv(61)

# Here, there is still numerical diffusion present, but it is less severe.

linearconv(71)

# Here the same pattern is present -- the wave is more square than in the previous runs.

linearconv(85)

# This doesn't look anything like our original hat function.

### What happened?

# To answer that question, we have to think a little bit about what we're actually implementing in code.

# In each iteration of our time loop, we use the existing data about our wave to estimate the speed of the
# wave in the subsequent time step.  Initially, the increase in the number of grid points returned more accurate answers.
# There was less numerical diffusion and the square wave looked much more like a square wave than it did in our first example.

# Each iteration of our time loop covers a time-step of length $\Delta t$, which we have been defining as 0.025.

# During this iteration, we evaluate the speed of the wave at each of the $x$ points we've created.  In the last plot, something has clearly gone wrong.

# What has happened is that over the time period $\Delta t$, the wave is travelling a distance which is greater than `Δx`.
# The length `Δx` of each grid box is related to the number of total points `nₓ`, so stability can
# be enforced if the $\Delta t$ step size is calculated with respect to the size of `Δx`.

# $$\sigma = \frac{u \Delta t}{\Delta x} \leq \sigma_{\max}$$

# where $u$ is the speed of the wave; $\sigma$ is called the **Courant number** and the
# value of $\sigma_{\max}$ that will ensure stability depends on the discretization used.

# In a new version of our code, we'll use the CFL number to calculate the appropriate time-step `Δt` depending on the size of `Δx`.

function linearconv(nₓ)
    Δx = 2 / (nₓ - 1)
    nₜ = 20    #nt is the number of timesteps we want to calculate
    c = 1
    σ = .5
    indexes = Int.((.5/Δx):(1/Δx))

    Δt = σ * Δx

    u = ones(nx)
    u[indexes] .= 2

    un = ones(nx)

    for n in range(nₜ):  #iterate through time
        uₙ = copy(u) ##copy the existing values of u into un
        for i in 1:nₓ
            u[i] = uₙ[i] - c * dt / dx * (uₙ[i] - uₙ[i-1])
        end
    end
    u
end

plot(range(start = 0, end = 2, length = nₓ), u)

linearconv(41)

linearconv(61)

linearconv(81)

linearconv(101)

linearconv(121)

# Notice that as the number of points `nx` increases, the wave convects a shorter and shorter distance.
# The number of time iterations we have advanced the solution at is held constant at `nₜ = 20`,
# but depending on the value of `nₓ` and the corresponding values of `Δx` and `Δt`, a shorter time window is being examined overall.

## Learn More
# -----
# ***

# It's possible to do rigurous analysis of the stability of numerical schemes, in some cases. Watch Prof. Barba's presentation of this topic in **Video Lecture 9** on You Tube.

# TODO
# from IPython.display import YouTubeVideo
# YouTubeVideo('Yw1YPBupZxU')

max_stars_repo_name: mkitti/CFDJulia

max_issues_repo_name: mkitti/CFDJulia

max_forks_repo_name: mkitti/CFDJulia