text:
push!(LOAD_PATH, "../src/")

using fier, Documenter

pages = [
    "Home" => "index.md",
    "Usage" => "usage.md",
    "Configuration" => "config.md"
]

makedocs(;
    modules = [fier],
    authors = "Kel Markert",
    repo = "https://github.com/SERVIR/fier-cli/blob/{commit}{path}#L{line}",
    sitename = "FIER CLI",
    pages = pages,
)

deploydocs(;
    repo = "github.com/SERVIR/fier-cli.git",
    devbranch = "main",
)

max_stars_repo_name: SERVIR/fier-cli

max_issues_repo_name: SERVIR/fier-cli

max_forks_repo_name: SERVIR/fier-cli


text:
# pfk activity function
function pfk_activity(atp,adp,f6p,fbp,amp, params)
	# (alpha,beta,gamma,delta);
	# (0,0,0,0);

	weight1=1;
	topa1=0;
	bottom1=1;

	# (0,0,0,1);
	weight2=atp^2/params["k4"];
	topa2=topa1;
	bottom2=bottom1+weight2;

	# (0,0,1,0);
	weight3=f6p^2/params["k3"];
	topa3=topa2+weight3;
	bottom3=bottom2+weight3;

	# (0,0,1,1);
	weight4=(f6p*atp)^2/(params["fatp"]*params["k3"]*params["k4"]);
	topa4=topa3+weight4;
	bottom4=bottom3+weight4;

	# (0,1,0,0);
	weight5=fbp/params["k2"];
	topa5=topa4;
	bottom5=bottom4+weight5;

	# (0,1,0,1)
	weight6=(fbp*atp^2)/(params["k2"]*params["k4"]*params["fbt"]);
	topa6=topa5;
	bottom6=bottom5+weight6;

	# (0,1,1,0)
	weight7=(fbp*f6p^2)/(params["k2"]*params["k3"]*params["ffbp"]);
	topa7=topa6+weight7;
	bottom7=bottom6+weight7;

	# (0,1,1,1)
	weight8=(fbp*f6p^2*atp^2)/(params["k2"]*params["k3"]*params["k4"]*params["ffbp"]*params["fbt"]*params["fatp"]);
	topa8=topa7+weight8;
	bottom8=bottom7+weight8;

	# (1,0,0,0);
	weight9=amp/params["k1"];
	topa9=topa8;
	bottom9=bottom8+weight9;

	# (1,0,0,1);
	weight10=(amp*atp^2)/(params["k1"]*params["k4"]*params["fmt"]);
	topa10=topa9;
	bottom10=bottom9+weight10;

	# (1,0,1,0);
	weight11=(amp*f6p^2)/(params["k1"]*params["k3"]*params["famp"]);
	topa11=topa10+weight11;
	bottom11=bottom10+weight11;

	# (1,0,1,1);
	weight12=(amp*f6p^2*atp^2)/(params["k1"]*params["k3"]*params["k4"]*params["famp"]*params["fmt"]*params["fatp"]);
	topa12=topa11+weight12;
	bottom12=bottom11+weight12;

	# (1,1,0,0)
	weight13=(amp*fbp)/(params["k1"]*params["k2"]);
	topa13=topa12;
	bottom13=bottom12+weight13;

	# (1,1,0,1);
	weight14=(amp*fbp*atp^2)/(params["k1"]*params["k2"]*params["k4"]*params["fbt"]*params["fmt"]);
	topa14=topa13;
	bottom14=bottom13+weight14;

	# (1,1,1,0) -- the most active state of the enzyme;
	weight15=(amp*fbp*f6p^2)/(params["k1"]*params["k2"]*params["k3"]*params["ffbp"]*params["famp"]);
	topa15=topa14;
	topb=weight15;
	bottom15=bottom14+weight15;

	# (1,1,1,1);
	weight16=(amp*fbp*f6p^2*atp^2)/(params["k1"]*params["k2"]*params["k3"]*params["k4"]*params["ffbp"]*params["famp"]*params["fbt"]*params["fmt"]*params["fatp"]);
	topa16=topa15+weight16;
	bottom16=bottom15+weight16;

	Jpfk= params["vpfk"]*(topb + params["kpfk"]*topa16)/bottom16;
end

# the main ode function
function sys(dy, y, params, t)
    #v, n, c, cer, cam, adp, f6p, fbp = y

    # calc atp/adp ratios
    rad     = sqrt(-4*y[6]^2+(params["atot"]-y[6])^2);
    atp     = (params["atot"]+rad-y[6])/2;
    mgadp   = 0.165*y[6]; # magnesium adp complex
    adp3m   = 0.135*y[6]; # adp for transport into matrix  (löffler petrides s. 238)
    atp4m   = 0.05 *atp;  # atp for transport into matrix (löffler petrides s. 238)
    amp     = y[6]^2/atp; # adenosine monophosphate

    # flux activation functions
    topo    = 0.08+0.89*mgadp^2/params["kdd"]^2+0.16*mgadp/params["kdd"] ;
    bottomo = (1+mgadp/params["kdd"])^2*(1+atp4m/params["ktt"] + adp3m/params["ktd"]) ;
    katpo   = topo/bottomo;                                                 # I_K(ATP) activation function
    minf    = 1/(1+exp((params["vm"]-y[1])/params["sm"]));                  # ca pmca activation function
    ninf    = 1/(1+exp((params["nin"]-y[1])/params["sn"]));                 # rectifing current acti
    qinf    = y[3]^2/(params["kd"]^2+y[3]^2);                               # I_K(Ca) activation functin
    # fluxes

    ik      = params["gk"]*y[2]*(y[1]-params["vk"]);                        # rectifying current
    ikca    = -params["gkca"]*qinf*(params["vk"]-y[1]);                     # ca dependent k current

    ikatp   = params["gkatpbar"]*katpo*(y[1]-params["vk"]);                 # k flux atp dependent
    ica     = params["gca"]*minf*(y[1]-params["vca"]);                      # ca flux for action potentials

    Jer     = params["kserca"]*y[3] - params["pleak"]*(y[4]-y[3]);          # ca flux density across er membrane
    Jm      = params["kuni"]*y[3] - params["knaca"]*(y[5]-y[3]);            # ca flux density across mitochonidium
    Jmem    = -(params["alpha"]/params["vcyt"]*ica + params["kpmca"]*y[3]); # ca flux density (cell membrane)

    # glycolitic oscillations
    Jpfk    = pfk_activity(atp, y[6], y[7], y[8], amp, params)              # pfk activity
    sinfty  = y[5]/(y[5]+params["kCaPDH"]);                                 # Michaelis-Menten function
    Jpdh    = params["vpdh"]*sinfty*sqrt(y[8]);                             # pdh activity

    #save_stats(amp, atp, katpo, minf, ninf, qinf, ik, ikca, Jer, ikatp, Jm, ica, Jmem, Jpfk, sinfty, Jpdh);
    dy[1] = -(ica + ik + ikca + ikatp)/params["Cm"]
    dy[2] = -(y[2]-ninf)/params["taun"]
    dy[3] = params["fca"]*(Jmem - Jm - Jer)
    dy[4] = params["fca"]*params["sigmaer"]*Jer
    dy[5] = params["fca"]*params["sigmam"]*Jm
    dy[6] = (atp-exp((1+2.2 * Jpdh/(Jpdh+0.05)) * (1-y[3]/0.35))*y[6])/params["taua"]
    dy[7] = 0.3*(params["Jgk"]-Jpfk)
    dy[8] = Jpfk-Jpdh/2
end

# the main ode function with constant pfk
function sys_const_pfk(dy, y, params, t)
    #v, n, c, cer, cam, adp, f6p, fbp = y

    # calc atp/adp ratios
    rad     = sqrt(-4*y[6]^2+(params["atot"]-y[6])^2);
    atp     = (params["atot"]+rad-y[6])/2;
    mgadp   = 0.165*y[6]; # magnesium adp complex
    adp3m   = 0.135*y[6]; # adp for transport into matrix  (löffler petrides s. 238)
    atp4m   = 0.05 *atp;  # atp for transport into matrix (löffler petrides s. 238)
    amp     = y[6]^2/atp; # adenosine monophosphate

    # flux activation functions
    topo    = 0.08+0.89*mgadp^2/params["kdd"]^2+0.16*mgadp/params["kdd"] ;
    bottomo = (1+mgadp/params["kdd"])^2*(1+atp4m/params["ktt"] + adp3m/params["ktd"]) ;
    katpo   = topo/bottomo;                                                 # I_K(ATP) activation function
    minf    = 1/(1+exp((params["vm"]-y[1])/params["sm"]));                  # ca pmca activation function
    ninf    = 1/(1+exp((params["nin"]-y[1])/params["sn"]));                 # rectifing current acti
    qinf    = y[3]^2/(params["kd"]^2+y[3]^2);                               # I_K(Ca) activation functin
    # fluxes

    ik      = params["gk"]*y[2]*(y[1]-params["vk"]);                        # rectifying current
    ikca    = -params["gkca"]*qinf*(params["vk"]-y[1]);                     # ca dependent k current
    ikatp   = params["gkatpbar"]*katpo*(y[1]-params["vk"]);                 # k flux atp dependent
    ica     = params["gca"]*minf*(y[1]-params["vca"]);                      # ca flux for action potentials

    Jer     = params["kserca"]*y[3] - params["pleak"]*(y[4]-y[3]);          # ca flux density across er membrane
    Jm      = params["kuni"]*y[3] - params["knaca"]*(y[5]-y[3]);            # ca flux density across mitochonidium
    Jmem    = -(params["alpha"]/params["vcyt"]*ica + params["kpmca"]*y[3]); # ca flux density (cell membrane)

    # glycolitic oscillations
    #Jpfk    = pfk_activity(atp, y[6], y[7], y[8], amp, params)              # pfk activity
    Jpfk    = 0.003
    sinfty  = y[5]/(y[5]+params["kCaPDH"]);                                 # Michaelis-Menten function
    Jpdh    = params["vpdh"]*sinfty*sqrt(y[8]);                             # pdh activity

    #save_stats(amp, atp, katpo, minf, ninf, qinf, ik, ikca, Jer, ikatp, Jm, ica, Jmem, Jpfk, sinfty, Jpdh);
    dy[1] = -(ica + ik + ikca + ikatp)/params["Cm"]
    dy[2] = -(y[2]-ninf)/params["taun"]
    dy[3] = params["fca"]*(Jmem - Jm - Jer)
    dy[4] = params["fca"]*params["sigmaer"]*Jer
    dy[5] = params["fca"]*params["sigmam"]*Jm
    dy[6] = (atp-exp((1+2.2 * Jpdh/(Jpdh+0.05)) * (1-y[3]/0.35))*y[6])/params["taua"]
    dy[7] = 0.3*(params["Jgk"]-Jpfk)
    dy[8] = Jpfk-Jpdh/2
end

# the main ode function with fixed value for katpo variable
function sys_katpo(dy, y, params, t)
    #v, n, c, cer, cam, adp, f6p, fbp = y

    # calc atp/adp ratios
    rad     = sqrt(-4*y[6]^2+(params["atot"]-y[6])^2);
    atp     = (params["atot"]+rad-y[6])/2;
    mgadp   = 0.165*y[6]; # magnesium adp complex
    adp3m   = 0.135*y[6]; # adp for transport into matrix  (löffler petrides s. 238)
    atp4m   = 0.05 *atp;  # atp for transport into matrix (löffler petrides s. 238)
    amp     = y[6]^2/atp; # adenosine monophosphate

    # flux activation functions
    topo    = 0.08+0.89*mgadp^2/params["kdd"]^2+0.16*mgadp/params["kdd"] ;
    bottomo = (1+mgadp/params["kdd"])^2*(1+atp4m/params["ktt"] + adp3m/params["ktd"]) ;
    katpo   = topo/bottomo;                                                 # I_K(ATP) activation function
    minf    = 1/(1+exp((params["vm"]-y[1])/params["sm"]));                  # ca pmca activation function
    ninf    = 1/(1+exp((params["nin"]-y[1])/params["sn"]));                 # rectifing current acti
    qinf    = y[3]^2/(params["kd"]^2+y[3]^2);                               # I_K(Ca) activation functin
    # fluxes
    katpo = 0.009
    ik      = params["gk"]*y[2]*(y[1]-params["vk"]);                        # rectifying current
    ikca    = -params["gkca"]*qinf*(params["vk"]-y[1]);                     # ca dependent k current
    ikatp   = params["gkatpbar"]*katpo*(y[1]-params["vk"]);                 # k flux atp dependent
    ica     = params["gca"]*minf*(y[1]-params["vca"]);                      # ca flux for action potentials

    Jer     = params["kserca"]*y[3] - params["pleak"]*(y[4]-y[3]);          # ca flux density across er membrane
    Jm      = params["kuni"]*y[3] - params["knaca"]*(y[5]-y[3]);            # ca flux density across mitochonidium
    Jmem    = -(params["alpha"]/params["vcyt"]*ica + params["kpmca"]*y[3]); # ca flux density (cell membrane)

    # glycolitic oscillations
    Jpfk    = pfk_activity(atp, y[6], y[7], y[8], amp, params)              # pfk activity
    sinfty  = y[5]/(y[5]+params["kCaPDH"]);                                 # Michaelis-Menten function
    Jpdh    = params["vpdh"]*sinfty*sqrt(y[8]);                             # pdh activity

    #save_stats(amp, atp, katpo, minf, ninf, qinf, ik, ikca, Jer, ikatp, Jm, ica, Jmem, Jpfk, sinfty, Jpdh);
    dy[1] = -(ica + ik + ikca + ikatp)/params["Cm"]
    dy[2] = -(y[2]-ninf)/params["taun"]
    dy[3] = params["fca"]*(Jmem - Jm - Jer)
    dy[4] = params["fca"]*params["sigmaer"]*Jer
    dy[5] = params["fca"]*params["sigmam"]*Jm
    dy[6] = (atp-exp((1+2.2 * Jpdh/(Jpdh+0.05)) * (1-y[3]/0.35))*y[6])/params["taua"]
    dy[7] = 0.3*(params["Jgk"]-Jpfk)
    dy[8] = Jpfk-Jpdh/2
end

ext: jl


text:
using StaticArrays: @SVector, normalize
using Accessors: Accessors, setproperties
export position
export direction

position(o) = o.position
direction(o) = o.direction
set_position(o, pos)  = setproperties(o, (position=pos,))
set_direction(o, dir) = setproperties(o, (direction=dir,))
set_position_direction(o, pos, dir) = setproperties(o, (position=pos, direction=dir))

Accessors.set(o, ::typeof(position), val) = set_position(o, val)
Accessors.set(o, ::typeof(direction), val) = set_direction(o, val)

function direction(p::AbstractParticle)
    @SVector [p.u, p.v, p.w]
end

has_z(::Type{<:AbstractParticle}) = true
has_z(::Type{<:EGSParticle})      = false

function position(p::P, z=nothing) where {P <: AbstractParticle}
    if (z == nothing)
        if has_z(P)
            return @SVector[p.x, p.y, p.z]
        else
            return @SVector[p.x, p.y]
        end
    else
        if has_z(P)
            @argcheck p.z ≈ z
            z = p.z
        end
        return @SVector[p.x, p.y, z]
    end
end

function set_direction(p::AbstractParticle, dir)
    u,v,w = dir
    setproperties(p, (u=u,v=v,w=w))
end

@inline function _xyz(pos, z)
    len = length(pos)
    if len == 3
        return pos
    elseif len == 2
        x,y = pos
        return x,y,z
    else
        throw(ArgumentError("Wrong length $(len))"))
    end
end

function set_position(p::AbstractParticle, pos)
    x,y,z = _xyz(pos, nothing)
    setproperties(p, (x=x,y=y,z=z))
end

function set_position(p::CompressedIAEAParticle, pos)
    set_position(IAEAParticle(p), pos)
end
function set_direction(p::CompressedIAEAParticle, dir)
    set_direction(IAEAParticle(p), dir)
end
function set_position_direction(p::CompressedIAEAParticle, pos, dir)
    set_position_direction(IAEAParticle(p), pos, dir)
end

function set_position_direction(p::AbstractParticle, pos, dir)
    x,y,z = _xyz(pos, nothing)
    u,v,w = dir
    setproperties(p, (x=x,y=y,z=z,u=u,v=v,w=w))
end


text:
@testset "Pairwise" begin
  𝒟 = PointSet(Matrix(1.0I, 3, 3))
  Γ = Variography.pairwise(GaussianVariogram(), 𝒟)
  @test eltype(Γ) == Float64
  @test issymmetric(Γ)

  𝒟 = PointSet(Matrix(1.0f0I, 3, 3))
  Γ_f = Variography.pairwise(GaussianVariogram(range=1f0, sill=1f0, nugget=0f0), 𝒟)
  @test eltype(Γ_f) == Float32
  @test issymmetric(Γ_f)

  𝒟 = CartesianGrid(10, 10)
  Γ = Variography.pairwise(GaussianVariogram(), view(𝒟, 1:5))
  @test size(Γ) == (5, 5)
  @test issymmetric(Γ)
  Γ = Variography.pairwise(GaussianVariogram(), view(𝒟, 1:3), view(𝒟, 7:10))
  @test size(Γ) == (3, 4)
  @test all(Γ .> 0)

  # arbitrary collections
  𝒟 = CartesianGrid(10, 10)
  𝒫 = centroid.(𝒟)
  Γ = Variography.pairwise(GaussianVariogram(), 𝒫)
  @test size(Γ) == (100, 100)
  @test issymmetric(Γ)
  Γ = Variography.pairwise(GaussianVariogram(), view(𝒫, 1:3), view(𝒫, 7:10))
  @test size(Γ) == (3, 4)
  @test all(Γ .> 0)
end


text:
module JsonSerializer

using JSON3

using SearchLight
using SearchLight.Serializer

function SearchLight.Serializer.serialize(value) :: String
  JSON3.write(value)
end

function SearchLight.Serializer.deserialize(::Type{T}, value)::T where {T}
  JSON3.read(value, T)
end

end


text:
#################### Core Model Functionality ####################

#################### Constructors ####################

function Model(; iter::Integer=0, burnin::Integer=0,
               samplers::Vector{Sampler}=Sampler[], nodes...)
  nodedict = Dict{Symbol, Any}()
  for (key, value) in nodes
    isa(value, AbstractDependent) ||
      throw(ArgumentError("nodes are not all Dependent types"))
    node = deepcopy(value)
    node.symbol = key
    nodedict[key] = node
  end
  m = Model(nodedict, Sampler[], ModelState[], iter, burnin, false, false)
  dag = ModelGraph(m)
  dependentkeys = keys(m, :dependent)
  terminalkeys = keys(m, :stochastic)
  for v in vertices(dag.graph)
    vkey = dag.keys[v]
    if vkey in dependentkeys
      m[vkey].targets = intersect(dependentkeys,
                                  gettargets(dag, v, terminalkeys))
    end
  end
  setsamplers!(m, samplers)
end

#################### Indexing ####################

Base.getindex(m::Model, nodekey::Symbol) = m.nodes[nodekey]

function Base.setindex!(m::Model, value, nodekey::Symbol)
  node = m[nodekey]
  if isa(node, AbstractDependent)
    node.value = value
  else
    m.nodes[nodekey] = convert(typeof(node), value)
  end
end

function Base.setindex!(m::Model, values::Dict, nodekeys::Vector{Symbol})
  for key in nodekeys
    m[key] = values[key]
  end
end

function Base.setindex!(m::Model, value, nodekeys::Vector{Symbol})
  length(nodekeys) == 1 || throw(BoundsError())
  m[first(nodekeys)] = value
end

Base.keys(m::Model) = collect(keys(m.nodes))

function Base.keys(m::Model, ntype::Symbol, at...)
  ntype == :block       ? keys_block(m, at...) :
  ntype == :all         ? keys_all(m) :
  ntype == :assigned    ? keys_assigned(m) :
  ntype == :dependent   ? keys_dependent(m) :
  ntype == :independent ? keys_independent(m) :
  ntype == :input       ? keys_independent(m) :
  ntype == :logical     ? keys_logical(m) :
  ntype == :monitor     ? keys_monitor(m) :
  ntype == :output      ? keys_output(m) :
  ntype == :source      ? keys_source(m, at...) :
  ntype == :stochastic  ? keys_stochastic(m) :
  ntype == :target      ? keys_target(m, at...) :
    throw(ArgumentError("unsupported node type $ntype"))
end

function keys_all(m::Model)
  values = Symbol[]
  for key in keys(m)
    node = m[key]
    if isa(node, AbstractDependent)
      push!(values, key)
      append!(values, node.sources)
    end
  end
  unique(values)
end

function keys_assigned(m::Model)
  if m.hasinits
    values = keys(m)
  else
    values = Symbol[]
    for key in keys(m)
      if !isa(m[key], AbstractDependent)
        push!(values, key)
      end
    end
  end
  values
end

function keys_block(m::Model, block::Integer=0)
  block == 0 ? keys_block0(m) : m.samplers[block].params
end

function keys_block0(m::Model)
  values = Symbol[]
  for sampler in m.samplers
    append!(values, sampler.params)
  end
  unique(values)
end

function keys_dependent(m::Model)
  values = Symbol[]
  for key in keys(m)
    if isa(m[key], AbstractDependent)
      push!(values, key)
    end
  end
  intersect(tsort(m), values)
end

function keys_independent(m::Model)
  deps = Symbol[]
  for key in keys(m)
    if isa(m[key], AbstractDependent)
      push!(deps, key)
    end
  end
  setdiff(keys(m, :all), deps)
end

function keys_logical(m::Model)
  values = Symbol[]
  for key in keys(m)
    if isa(m[key], AbstractLogical)
      push!(values, key)
    end
  end
  values
end

function keys_monitor(m::Model)
  values = Symbol[]
  for key in keys(m)
    node = m[key]
    if isa(node, AbstractDependent) && !isempty(node.monitor)
      push!(values, key)
    end
  end
  values
end

function keys_output(m::Model)
  values = Symbol[]
  dag = ModelGraph(m)
  for v in vertices(dag.graph)
    vkey = dag.keys[v]
    if isa(m[vkey], AbstractStochastic) && !any_stochastic(dag, v, m)
      push!(values, vkey)
    end
  end
  values
end

keys_source(m::Model, nodekey::Symbol) = m[nodekey].sources

function keys_source(m::Model, nodekeys::Vector{Symbol})
  values = Symbol[]
  for key in nodekeys
    append!(values, m[key].sources)
  end
  unique(values)
end

function keys_stochastic(m::Model)
  values = Symbol[]
  for key in keys(m)
    if isa(m[key], AbstractStochastic)
      push!(values, key)
    end
  end
  values
end

function keys_target(m::Model, block::Integer=0)
  block == 0 ? keys_target0(m) : m.samplers[block].targets
end

function keys_target0(m::Model)
  values = Symbol[]
  for sampler in m.samplers
    append!(values, sampler.targets)
  end
  intersect(keys(m, :dependent), values)
end

keys_target(m::Model, nodekey::Symbol) = m[nodekey].targets

function keys_target(m::Model, nodekeys::Vector{Symbol})
  values = Symbol[]
  for key in nodekeys
    append!(values, m[key].targets)
  end
  intersect(keys(m, :dependent), values)
end

#################### Display ####################

function Base.show(io::IO, m::Model)
  showf(io, m, Base.show)
end

function showall(io::IO, m::Model)
  showf(io, m, Base.showall)
end

function showf(io::IO, m::Model, f::Function)
  print(io, "Object of type \"$(summary(m))\"\n")
  width = displaysize()[2] - 1
  for node in keys(m)
    print(io, string("-"^width, "\n", node, ":\n"))
    f(io, m[node])
    println(io)
  end
end

#################### Auxiliary Functions ####################

function names(m::Model, monitoronly::Bool)
  values = AbstractString[]
  for key in keys(m, :dependent)
    nodenames = names(m, key)
    v = monitoronly ? nodenames[m[key].monitor] : nodenames
    append!(values, v)
  end
  values
end

function names(m::Model, nodekey::Symbol)
  node = m[nodekey]
  unlist(node, names(node))
end

function names(m::Model, nodekeys::Vector{Symbol})
  values = AbstractString[]
  for key in nodekeys
    append!(values, names(m, key))
  end
  values
end


text:
using AbstractAlgebra

@testset "Groebner modular Integers" begin
    R, (x,) = PolynomialRing(ZZ, ["x"], ordering=:degrevlex)
    @test Groebner.groebner([x]) == [x]
    @test Groebner.groebner([4x]) == [x]
    @test Groebner.groebner([x + 1]) == [x + 1]
    @test Groebner.groebner([4x + 1]) == [4x + 1]
    @test Groebner.groebner([x + 5]) == [x + 5]
    @test Groebner.groebner([x + 2^10]) == [x + 2^10]
    @test Groebner.groebner([x + 2^30]) == [x + 2^30]
    @test Groebner.groebner([x + 2^30 + 3]) == [x + 2^30 + 3]
    @test Groebner.groebner([x + 2^31 - 1]) == [x + 2^31 - 1]
    @test Groebner.groebner([x + 2^31 - 1, x^2]) == [1]
    @test Groebner.groebner([3232323x + 7777]) == [3232323x + 7777]
    @test Groebner.groebner([(2^31-1)x + 1]) == [2147483647x + 1]
    @test Groebner.groebner([x + 2147483647]) == [x + 2147483647]

    # consequtive primes: 2^31-1, 2147483659, 2147483693
    # consequtive primes: 2^30+3, 1073741831, 1073741833
    R, (x, y) = PolynomialRing(ZZ, ["x", "y"], ordering=:degrevlex)
    @test Groebner.groebner([4x + 1, 3y + 6]) == [y + 2, 4x + 1]
    @test Groebner.groebner([3y + 6, 4x + 1]) == [y + 2, 4x + 1]
    @test Groebner.groebner([10x + 10y, 10y + 1]) == [10y + 1, 10x - 1]
    @test Groebner.groebner([(2^31-1)x*y + 2147483659x + 2147483693y]) == [(2^31-1)x*y + 2147483659x + 2147483693y]
    @test Groebner.groebner([ (2^30+3)x*y^2 + 1073741831x^2 + 1073741833y + 2^31-1]) == [(2^30+3)x*y^2 + 1073741831x^2 + 1073741833y + 2^31-1]

    Groebner.groebner([4x + 1, 3y + 6])
end

@testset "Groebner modular Rationals" begin
    R, (x,) = PolynomialRing(QQ, ["x"], ordering=:degrevlex)
    @test Groebner.groebner([x]) == [x]
    @test Groebner.groebner([4x]) == [x]
    @test Groebner.groebner([x + 1]) == [x + 1]
    @test Groebner.groebner([4x + 1]) == [x + 1//4]
    @test Groebner.groebner([x + 5]) == [x + 5]
    @test Groebner.groebner([x + 2^10]) == [x + 2^10]
    @test Groebner.groebner([x + 2^30]) == [x + 2^30]
    @test Groebner.groebner([x + 2^30 + 3]) == [x + 2^30 + 3]
    @test Groebner.groebner([x + 2^31 - 1]) == [x + 2^31 - 1]
    @test Groebner.groebner([x + 2^31 - 1, x^2]) == [1]
    @test Groebner.groebner([(3232323//7777)x + 7777//3232323]) == [x + 60481729//10447911976329]
    @test Groebner.groebner([((2^31-1)//1)x + 1]) == [x + 1//2147483647]
    @test Groebner.groebner([(1//(2^31-1))x + 1]) == [x + 2147483647]
    @test Groebner.groebner([1//(2^30+3)*x^2 + (2^30+3)x + 1//(1073741831)]) == [x^2 + 1152921511049297929x + (2^30+3)//1073741831]

    R, (x, y) = PolynomialRing(QQ, ["x", "y"], ordering=:degrevlex)

    fs = [x, y]
    G = Groebner.groebner(fs)
    @test G ≂ [y, x]

    fs = [5y, 3//4*x]
    G = Groebner.groebner(fs)
    @test G ≂ [y, x]

    fs = [x^2 - (1//6)*y, x*y]
    G = Groebner.groebner(fs)
    @test G ≂ [y^2, x*y, x^2 - (1//6)y]

    fs = [QQ(11, 3)*x^2*y - QQ(2, 4)*x*y - QQ(1, 7)*y,
            QQ(1, 1)*x*y + QQ(7, 13)*x]
    G = Groebner.groebner(fs)
    @test G ≂ [y^2 + 7//13*y,
                x*y + 7//13*x,
                x^2 - 3//22*x + 39//539*y]

    #=
    root = Groebner.change_ordering(Groebner.rootn(3, ground=QQ), :degrevlex)
    gb = Groebner.groebner(root)
    @test Groebner.isgroebner(gb)

    root = Groebner.change_ordering(Groebner.rootn(4, ground=QQ), :degrevlex)
    gb = Groebner.groebner(root)
    @test Groebner.isgroebner(gb)

    root = Groebner.change_ordering(Groebner.rootn(8, ground=QQ), :degrevlex)
    gb = Groebner.groebner(root)
    @test Groebner.isgroebner(gb)
    =#

    #=
    noon = Groebner.change_ordering(Groebner.noon3(ground=QQ), :degrevlex)
    gb = Groebner.groebner(noon)
    @test Groebner.isgroebner(gb)

    noon = Groebner.change_ordering(Groebner.noonn(4, ground=QQ), :degrevlex)
    gb = Groebner.groebner(noon)
    @test Groebner.isgroebner(gb)
    =#
end

@testset "Groebner modular corner cases" begin

    R, (x, y, z, w) = PolynomialRing(QQ, ["x", "y", "z", "w"], ordering=:degrevlex)

    fs = [(12345678//12347)x,
          (222222221111123//2131232232097)y + z]
    G = Groebner.groebner(fs)
    @test G == [y + 2131232232097//222222221111123*z, x]

    fs = [
           (224//225)*x + y,
        x^2 + 1]
    ans = [x + 225//224*y, y^2 + 50176//50625]
    @test Groebner.groebner(fs) == ans

    fs = [
        (12345//12345678)x + y,
        x^2 + 1]
    ans = [x + 4115226//4115*y, y^2 + 16933225//16935085031076]
    @test Groebner.groebner(fs) == ans

    fs = [
        (12345//12345678)x + y,
        x^2 + z^2,
        y^2 + w^2,
        x + 2*y + 3*z + 4*w
    ]
    G = Groebner.groebner(fs)

    @test G ≂ [w^3,
                z*w + 1692834523553//4063974*w^2,
                z^2 - 16935085031076//16933225*w^2,
                y - 12345//4106996*z - 4115//1026749*w,
                x + 6172839//2053498*z + 4115226//1026749*w]

    # what if we are unlucky to start from initial prime in coefficients
    G = [
        (2^30 + 3)x,
        (1//(2^30 + 3))y
    ]
    G = Groebner.groebner(G)
    @test G == [y, x]

    G = [
        (2^31 - 1)*x + y,
        (1//(2^31 - 1))*y + 1
    ]
    G = Groebner.groebner(G)
    @test G == [y + 2147483647, x - 1]

    R, (x, y, z, w) = PolynomialRing(QQ, ["x", "y", "z", "w"], ordering=:lex)
    fs = [(12345678//12347)x,
          (222222221111123//2131232232097)y + z]
    G = Groebner.groebner(fs)
    @test G == [y + 2131232232097//222222221111123*z, x]
end


text:
@testset "System" begin
    ρ = 0.5
    kT = 1.0
    parts = 10
    cutoff = 1.2

    # Without specifying the cutoff
    s = System(ρ, kT, parts, cutoff)
    @test length(s.xpos) == parts
    @test length(s.xpos[1]) == 3

    # Check that the position vectors are different among them
    @test s.xpos[1] != s.xpos[2] && s.xpos[fld(parts, 2)] != s.xpos[parts]
end

@testset "Simulations" begin
    ρ = 0.5
    kT = 1.0
    parts = 10
    cutoff = 1.2
    s = System(ρ, kT, parts, cutoff)
    sim = Simulation(s, NVT(0.35), LennardJones())

    @test sim isa Simulation
end


text:
include("loadAlessioData.jl")
using .loadAlessioData, DelimitedFiles, FortranFiles, Printf, ProgressBars, PyPlot, Peaks, Dierckx, Interpolations

function main(; plotting=false, movie=false, outputObs=true)

	# where is the data
	basepath = "./2011-05-28_Rec78-103_Pace_Apex/"
	inds = getPossibleIndices(basepath)
	ind = rand(inds) # This is just useful for debugging; normally you would choose an index

	# print some diagnostic notices
	print("Loading data for $(ind)...\t")

	# load in some data
	states = loadFile(basepath, ind)

	print("Done.\n")

	# form the space and time masks
	spaceMasks = [spatialMask(state) for state in states]
	timeMasks  = [temporalMask(state) for state in states]
	mutualSpaceMask = mutualMask(spaceMasks)
	mutualTimeMask = mutualMask(timeMasks)

	# determine mutualSpaceMask centroid for square center
	center = spatialMaskCentroid(mutualSpaceMask)
	corners = loadAlessioData.square(theta=0.0, l=3.0, center=center)

	# make some plots
	if plotting
	# sampling position is randomly sampled and checked against mutual space mask
		samplePoint = [rand(1:128), rand(1:128)]
		while mutualSpaceMask[samplePoint[1],samplePoint[2]] == false
			samplePoint = [rand(1:128), rand(1:128)]
		end
		t, V = sampleData(states, samplePoint, mutualTimeMask, mutualSpaceMask)
		fig, axs = plt.subplots(1,3,figsize=(12,3),constrained_layout=true)
		plotTrace!(fig, axs, t, V; threshold=0.25)
		plt.savefig("./trace/$(ind)_trace.svg", bbox_inches="tight")
		plt.close()

		# plotting a single state at a prescribed time
		sampleTime = rand((1:length(mutualTimeMask))[mutualTimeMask])
		fig, axs = plt.subplots(1,2,figsize=(8,3),constrained_layout=true)
		plotState!(fig, axs, sampleTime, states, spaceMasks, timeMasks, [size(state) for state in states], mutualSpaceMask, mutualTimeMask; drawColorbar=true)
		plotSamplingSquare!(fig, axs, corners)
		plt.savefig("./state/$(ind)_state.svg", bbox_inches="tight")
		plt.close()
	end

	# make a movie of the dynamics
	if movie
		recordingMovie(states, spaceMasks, timeMasks, mutualSpaceMask, mutualTimeMask; figname="$(ind)")
	end

	# output observations
	if outputObs
		# read in a single observation file as a template
		obs = readObsFile("./obs/template.dat")

		# generate the splines for the spatiotemporal states
		stateSplines = stateInterpolants(states, mutualTimeMask, mutualSpaceMask)

		# interpolate into stateSplines and mutate the observation template and save to new observation files
		newObsDir = "./obs/$(ind)/"; mkpath(newObsDir);
		sampleTimes = 2.0.*(1:length(mutualTimeMask))[mutualTimeMask]
		generateObservationFileSequence!(obs, stateSplines, newObsDir, "", sampleTimes; center = center)
	end

	return nothing
end

main(;plotting=true, outputObs=false)

max_stars_repo_name: cmarcotte/puttingTheDataInDA

max_issues_repo_name: cmarcotte/puttingTheDataInDA

max_forks_repo_name: cmarcotte/puttingTheDataInDA


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to string processing.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# This was adapted from Julia `escape_string` function. In case of LaTeX, we
# should not escape the sequence `\\`, i.e. the backslash that is already
# escaped.
function _str_latex_escaped(io::IO, s::AbstractString, esc::String = "")
    a = Iterators.Stateful(s)
    for c in a
        if c in esc
            print(io, '\\', c)
        elseif isascii(c)
            c == '\0'          ? print(io, escape_nul(peek(a))) :
            c == '\e'          ? print(io, "\\e") :
          # c == '\\'          ? print(io, "\\\\") :
            '\a' <= c <= '\r'  ? print(io, '\\', "abtnvfr"[Int(c)-6]) :
            isprint(c)         ? print(io, c) :
                                 print(io, "\\x", string(UInt32(c), base = 16, pad = 2))
        elseif !isoverlong(c) && !ismalformed(c)
            isprint(c)         ? print(io, c) :
            c <= '\x7f'        ? print(io, "\\x", string(UInt32(c), base = 16, pad = 2)) :
            c <= '\uffff'      ? print(io, "\\u", string(UInt32(c), base = 16, pad = need_full_hex(peek(a)) ? 4 : 2)) :
                                 print(io, "\\U", string(UInt32(c), base = 16, pad = need_full_hex(peek(a)) ? 8 : 4))
        else # malformed or overlong
            u = bswap(reinterpret(UInt32, c))
            while true
                print(io, "\\x", string(u % UInt8, base = 16, pad = 2))
                (u >>= 8) == 0 && break
            end
        end
    end
end

_str_latex_escaped(s::AbstractString, esc::String = "") =
    sprint(_str_latex_escaped, s, esc, sizehint=lastindex(s))


text:
# Subsetting info-table with quality
function getQScoreOverQ(df::DataFrame, qscore::Int64, totallength::Int64)
	df_subset = df[df.quality .>= qscore, :]
	readNum = size(df_subset, 1)
	baseNum = totalLen(df_subset.length)
	basePerc = round(baseNum * 100 / totallength, digits=2)
	return(readNum, baseNum, "$basePerc%")
end

# Plotting functions
function plotSquiggle(signals::Array{Int16,1})
	gr()
	timepoints = [i for i in 1:length(signals)] ./ 4000
	plot(timepoints, signals, w=3, xlabel = "Time (seconds)", ylabel= "Detected Current (pA)", legend = false)
end

function plotReadQual2IdentScatter(df::DataFrames.DataFrame, outputDir::String)
	gr()
	outputfile = joinpath(outputdir, "read_quality_vs_identity_dot.png")
	plot(df.quality,
		df.identity,
		seriestype = :scatter,
		title = "Read Quality vs Identity",
		xlabel = "Phred Score",
		ylabel = "Identity (%)",
		dpi = 300)
	savefig(outputfile)
end

function plotReadLen2QualScatter(df::DataFrames.DataFrame, N50::Int64, outputDir::String)
	gr()
	outputfile = joinpath(outputDir, "read_length_vs_quality_dot.png")
	if "identity" in names(df)
		plot(df.length,
			df.quality,
			zcolor = df.identity,
			seriestype = :scatter,
			title = "Read Length vs Quality",
			xlabel = "Read Length (bp)",
			ylabel = "Phred Score",
			lab = "Identity (%)",
			dpi = 300)
	else
		plot(df.length,
			df.quality,
			seriestype = :scatter,
			title = "Read Length vs Quality",
			legend = false,
			xlabel = "Read Length (bp)",
			ylabel = "Phred Score",
			dpi = 300)
	end
	vline!([N50], lw = 2, lab = "N50")
	savefig(outputfile)
end

function plotReadLen2QualHistogram2D(df::DataFrames.DataFrame, outputDir::String)
	gr()
	outputfile = joinpath(outputDir, "read_length_vs_quality_histogram.png")
	histogram2d(df.length, df.quality, nbins = 50, dpi=300, xlabel="Read Length (bp)", ylabel="Phred Score", title="Read Length vs Quality")
	savefig(outputfile)
end

function plotReadQual2IdentHistogram2D(df::DataFrames.DataFrame, outputDir::String)
	gr()
	outputfile = joinpath(outputDir, "read_quality_vs_identity_histogram.png")
	histogram2d(df.quality, df.identity, nbins = 50, dpi=300, xlabel="Phred Score", ylabel="Identity (%)", title="Read Quality vs Identity")
	savefig(outputfile)
end

function plotReadLenDist(lengthRecords::Array{Int64,1}, N50::Int64, outputDir::String)
	gr()
	outputfile = joinpath(outputDir, "read_length_distribution_histogram.png")
	plot(lengthRecords,
		 seriestype=:histogram,
		 bins=400,
		 minorticks=true,
		 legend=false,
		 xlabel = "Read Length (bp)",
		 ylabel = "Read Counts",
		 title = "Read Length Distribution",
		 dpi = 300)
	vline!([N50], lw=2, lab="N50")
	savefig(outputfile)
end

function plotGCcontentDist(GCRecords::Array{Float64,1}, avgGCcontent::Float64, outputDir::String)
	transformedGCrecords = GCRecords .* 100
	gr()
	outputfile = joinpath(outputDir, "GC_distribution_histogram.png")
	plot(transformedGCrecords,
		 seriestype=:histogram,
		 bins=400,
		 minorticks=true,
		 legend=false,
		 xlabel = "GC Content (%)",
		 ylabel = "Read Counts",
		 title = "GC-content Distrbution",
		 dpi = 300)
		vline!([avgGCcontent], lw=2, lab="Mean GC-content")
	savefig(outputfile)
end

# Generating statics summary files
function generateStatSummary(df::DataFrames.DataFrame, totalLength::Int64, N50::Int64, outputDir::String)::DataFrames.DataFrame
	stat_summary = describe(df)
	readNum = size(df, 1)
	avg_gc_content = mean_gc_content(df.gc_content)
	readsQ7, basesQ7, percQ7 = getQScoreOverQ(df, 7, totalLength)
	readsQ10, basesQ10, percQ10 = getQScoreOverQ(df, 10, totalLength)
	readsQ12, basesQ12, percQ12 = getQScoreOverQ(df, 12, totalLength)
	readsQ15, basesQ15, percQ15 = getQScoreOverQ(df, 15, totalLength)
	output_df = DataFrame(Property = AbstractString[], Value = Number[])
	if ncol(df) == 4
		meanQual, meanLen, meanIdent = tuple(round.(stat_summary[!,:mean], digits=1)...)
		medianQual, medianLen, medianIdent = tuple(stat_summary[!,:median]...)
		push!(output_df, ("Average Identity", meanIdent))
		push!(output_df, ("Median Identity", medianIdent))
		meanIdenttxt = @sprintf "Mean Identity: %24.1f" meanIdent
		medianIdenttxt = @sprintf "Median Identity: %22.1f" medianIdent
	else
		meanQual, meanLen = tuple(round.(stat_summary[!,:mean], digits=1)...)
		medianQual, medianLen = tuple(stat_summary[!,:median]...)
	end
	push!(output_df, ("Mean Read Length", meanLen))
	push!(output_df, ("Median Read Length", medianLen))
	push!(output_df, ("Mean Read Quality", meanQual))
	push!(output_df, ("Median Read Quality", medianQual))
	push!(output_df, ("Mean GC-Content", avg_gc_content))
	push!(output_df, ("Read N50", N50))
	push!(output_df, ("Read Numbers", readNum))
	push!(output_df, ("Total Bases", totalLength))

	meanQualtxt = @sprintf "Mean Read Quality: %27.1f" meanQual
	meanLentxt = @sprintf "Mean Read Length: %28.1f" meanLen
	medianQualtxt = @sprintf "Median Read Quality: %25.1f" medianQual
	medianLentxt = @sprintf "Median Read Length: %26.1f" medianLen
	readN50txt = @sprintf "Read N50: %36s" format(N50, commas=true)
	readNumtxt = @sprintf "Read Number: %33s" format(readNum, commas=true)
	totalBasestxt = @sprintf "Total Bases: %33s" format(totalLength, commas=true)
	avgGCtxt = @sprintf "Mean GC-Content (%%): %25s" string(avg_gc_content*100, " %")
	readsQ7txt = @sprintf "Reads with quality score > 7: %16s" format(readsQ7, commas=true)
	basesQ7txt = @sprintf "Bases with quality score > 7: %16s" format(basesQ7, commas=true)
	percQ7txt = @sprintf "Percentage of bases (Q>7): %19s" percQ7
	readsQ10txt = @sprintf "Reads with quality score > 10: %15s" format(readsQ10, commas=true)
	basesQ10txt = @sprintf "Bases with quality score > 10: %15s" format(basesQ10, commas=true)
	percQ10txt = @sprintf "Percentage of bases (Q>10): %18s" percQ10
	readsQ12txt = @sprintf "Reads with quality score > 12: %15s" format(readsQ12, commas=true)
	basesQ12txt = @sprintf "Bases with quality score > 12: %15s" format(basesQ12, commas=true)
	percQ12txt = @sprintf "Percentage of bases (Q>12): %18s" percQ12
	readsQ15txt = @sprintf "Reads with quality score > 15: %15s" format(readsQ15, commas=true)
	basesQ15txt = @sprintf "Bases with quality score > 15: %15s" format(basesQ15, commas=true)
	percQ15txt = @sprintf "Percentage of bases (Q>15): %18s" percQ15

	outputfile = joinpath(outputDir, "statistics_summary.txt")
	open(outputfile, "w") do io
		if ncol(df) == 4
			write(io, "$meanQualtxt\n$meanLentxt\n$meanIdenttxt\n$medianQualtxt\n$medianLentxt\n$medianIdenttxt\n$readN50txt\n$avgGCtxt\n$readNumtxt\n$totalBasestxt\n")
		else
			write(io, "$meanQualtxt\n$meanLentxt\n$medianQualtxt\n$medianLentxt\n$readN50txt\n$avgGCtxt\n$readNumtxt\n$totalBasestxt\n")
		end
		write(io, "$readsQ7txt\n$basesQ7txt\n$percQ7txt\n")
		write(io, "$readsQ10txt\n$basesQ10txt\n$percQ10txt\n")
		write(io, "$readsQ12txt\n$basesQ12txt\n$percQ12txt\n")
		write(io, "$readsQ15txt\n$basesQ15txt\n$percQ15txt\n")
	end

	plotReadLen2QualScatter(df, N50, outputDir)
	plotReadLen2QualHistogram2D(df, outputDir)
	plotReadLenDist(df.length, N50, outputDir)
	plotGCcontentDist(df.gc_content, avg_gc_content*100, outputDir)

	if ncol(df) == 4
		plotReadQual2IdentScatter(df, outputDir)
		plotReadQual2IdentHistogram2D(df, outputDir)
	end

	return output_df
end


text:
"""
todot

Create dot file
"""

export todot

import DD: todot
using UUIDs: uuid1

function todot(top::AbstractFaultTreeNode)
    io = IOBuffer()
    visited = Set{AbstractFaultTreeNode}()
    println(io, "digraph { layout=dot; overlap=false; splines=true; node [fontsize=10];")
    _todot!(top, visited, io)
    println(io, "}")
    return String(take!(io))
end

function _todot!(f::FaultTreeEvent, visited::Set{AbstractFaultTreeNode}, io::IO)
#    id = uuid1()
    id = "obj$(objectid(f))"
    println(io, "\"$(id)\" [shape = circle, label = \"$(f.var)\"];")
    id
end

function _todot!(f::AbstractFaultTreeOperation, visited::Set{AbstractFaultTreeNode}, io::IO)
    id = "obj$(objectid(f))"
    (f in visited) && return id
    push!(visited, f)
    println(io, "\"$(id)\" [shape = square, label = \"$(f.op)\"];")
    for x = f.args
        dest = _todot!(x, visited, io)
        println(io, "\"$(id)\" -> \"$(dest)\";")
    end
    id
end

function _todot!(f::FaultTreeKoutofN, visited::Set{AbstractFaultTreeNode}, io::IO)
    id = "obj$(objectid(f))"
    (f in visited) && return id
    push!(visited, f)
    println(io, "\"$(id)\" [shape = square, label = \"$(f.op) $(f.k)\"];")
    for x = f.args
        dest = _todot!(x, visited, io)
        println(io, "\"$(id)\" -> \"$(dest)\";")
    end
    id
end


text:
# Autogenerated wrapper script for Libtiff_jll for i686-w64-mingw32
export libtiff

using JpegTurbo_jll
using Zlib_jll
using Zstd_jll
JLLWrappers.@generate_wrapper_header("Libtiff")
JLLWrappers.@declare_library_product(libtiff, "libtiff-5.dll")
function __init__()
    JLLWrappers.@generate_init_header(JpegTurbo_jll, Zlib_jll, Zstd_jll)
    JLLWrappers.@init_library_product(
        libtiff,
        "bin\\libtiff-5.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# electromagnetics utility functions

Mu_vacuum = 1.25663706212e-6  # vacuum magnetic permeability
Eps_vacuum = 8.8541878128e-12  # vacuum electric permittivity

"""
    Convert between wavelength and frequency
    It should be similar to Matlab function "z_convert_wavelength_freq"
"""
function convert_wavelength_freq(
        wl_or_freq_input::R; input_unit="", output_unit="", decimals_to_round=nothing
    ) where R <: Real
    c = 299792458
    # speed of light
    h_eV = 4.135667662e-15
    # [eV.s]

    input_unit = lowercase(input_unit)
    output_unit = lowercase(output_unit)

    # convert input to Hz
    if input_unit == "hz"
        freq_Hz = wl_or_freq_input

    elseif input_unit == "rad/s"
        freq_Hz = wl_or_freq_input / (2 * pi)

    elseif input_unit == "scufffrequnit"
        freq_radpers = wl_or_freq_input * 3e14
        freq_Hz = freq_radpers / (2 * pi)

    elseif input_unit == "nm"
        wl_m = wl_or_freq_input * 1e-9
        freq_Hz = c / wl_m

    elseif input_unit == "um"
        wl_m = wl_or_freq_input * 1e-6
        freq_Hz = c / wl_m

    elseif input_unit == "mm"
        wl_m = wl_or_freq_input * 1e-3
        freq_Hz = c / wl_m

    elseif input_unit == "cm"
        wl_m = wl_or_freq_input * 1e-2
        freq_Hz = c / wl_m

    elseif input_unit == "m"
        wl_m = wl_or_freq_input
        freq_Hz = c / wl_m

    elseif input_unit == "ev"
        wl_m = c / wl_or_freq_input * h_eV
        freq_Hz = c / wl_m

    elseif input_unit in ["cm^-1", "cm-1"]
        wl_m = 1 / wl_or_freq_input * 1e-2
        freq_Hz = c / wl_m

    else
        throw(DomainError("I don't know how to calculate for this `$input_unit` input unit"))
    end

    # convert Hz into the desired unit
    if output_unit == "hz"
        wl_or_freq_output = freq_Hz

    elseif output_unit == "rad/s"
        wl_or_freq_output = freq_Hz * (2 * pi)

    elseif output_unit == "scufffrequnit"
        wl_or_freq_output = freq_Hz * (2 * pi) / 3e14

    elseif output_unit == "nm"
        wl_m = c / freq_Hz
        wl_or_freq_output = wl_m * 1e9

    elseif output_unit == "um"
        wl_m = c / freq_Hz
        wl_or_freq_output = wl_m * 1e6

    elseif output_unit == "mm"
        wl_m = c / freq_Hz
        wl_or_freq_output = wl_m * 1e3

    elseif output_unit == "cm"
        wl_m = c / freq_Hz
        wl_or_freq_output = wl_m * 1e2

    elseif output_unit == "m"
        wl_m = c / freq_Hz
        wl_or_freq_output = wl_m

    elseif output_unit == "ev"
        wl_m = c / freq_Hz
        wl_or_freq_output = c / wl_m * h_eV

    elseif output_unit in ["cm^-1", "cm-1"]
        wl_m = c / freq_Hz
        wl_or_freq_output = 1 / (wl_m*1e2)

    else
        throw(DomainError("I don't know how to calculate for this `$output_unit` output unit"))
    end

    if decimals_to_round !== nothing
        wl_or_freq_output = round(wl_or_freq_output, digits=Int(decimals_to_round))
    end

    return wl_or_freq_output

end

"""
    It needs one definition for permittivity (either Eps or Eps_r), and one definition for permeability (either Mu or Mu_r)
"""
function get_SpeedOfLight(;Eps = nothing, Mu = nothing, Eps_r = nothing, Mu_r = nothing)
    if (Mu == nothing) && (Mu_r !== nothing)
        Mu = Mu_r * Mu_vacuum
    else
        throw(DomainError("You must input either `Mu` or `Mu_r` kewword arguments."))
    end

    if (Eps == nothing) && (Eps_r !== nothing)
        Eps = Eps_r * Eps_vacuum
    else
        throw(DomainError("You must input either `Eps` or `Eps_r` kewword arguments."))
    end

    return 1 / sqrt(Mu * Eps)
end

"""
    get Refractive Index from electric Permittivity
"""
function get_RefractiveIndex_from_Permittivity(Eps_r)
    return sqrt(Eps_r)
end

"""
It needs one definition for permittivity (either Eps or Eps_r), and one definition for permeability (either Mu or Mu_r)
"""
function get_WaveVector(
        wl_or_freq_input;
        input_unit="",
        Eps = nothing, Mu = nothing,
        Eps_r = nothing, Mu_r = nothing
    )
    Omega = convert_wavelength_freq(
        wl_or_freq_input; input_unit=input_unit, output_unit="rad/s"
    )
    c = get_SpeedOfLight(;Eps = Eps, Mu = Mu, Eps_r = Eps_r, Mu_r = Mu_r)
    k = Omega / c
    return k
end


text:
module train

#Use this to print
#=
using Weave
codeName = "train"
weave(Pkg.dir("$(pwd())\\$(codeName).jl"),
  informat="script",
  out_path = "$(pwd())\\output\\$(codeName)_Appendix.html",
  doctype = "md2html")
=#

#NOTE: To prepare the data, conform missing data to "". That means eliminating
#missing codes "N/A", "missing", and "--"

if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, CSV, StatsBase, GZip, ForwardDiff,# Query,
  Gadfly, NLopt, HCubature, StaticArrays, Measures, Formatting

importall CTIO, CTReg, CTStat, CTVAR, CTNumerical

#############parameters
const VOL_PERIODS = 1

#path constants
const DATA_PATH = pwd() * "\\data" #path where we store the data
const OUTPUT_PATH = pwd() * "\\output"
const DATE_FORMAT = "m/d/yyyy"
const CHINA_DS_NAME = "ChinaDS1.1"
const CHINA_DS_NAME_L = "$(CHINA_DS_NAME)_L"
const CHINA_DS_NAME_W = "$(CHINA_DS_NAME)_W"
const R_TEST_PATH = "C:\\Users\\Clinton Tepper\\Dropbox\\Projects\\InterconnectRTester\\trainTester\\trainTester"
const ACCEPTABLE_RANGE = #use this to check if values are in acceptable ranges
    Dict(:value=>(-10., 10.), :marketCap=>(0.1*10. ^ 9.,Inf))
const DROP_OUTSIDE_RANGE = true
const MIN_SZSECOM_MARKETCAP = 6. * 10. ^ 9.
const KEEP_MISSING_MARKETCAP = false

const MIN_SAMPLES_IN_INDUSTRY = 5
#const MIN_SAMPLES_INTERACTED = 50 #obsolete
const MIN_COMPANY_QUARTERS = 1 #threshold for dropping companies
const SET_MISSING_IF_BELOW = true
const KEEP_LARGEST_WEIGHT_IF_DUP = true # keep the largest weight if multiple shareclasses

const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
const HEADER_NAME = "header.tex"
const DECIMALS = 2
const NUM_FORMAT = "{:.$(DECIMALS)f}"
#const PlotContainer = Union{Plot,Gadfly.Compose.Context,Vector{Plot}} #holds graphs
const OVERRIDE_STAR_STRINGS = ["\\ostar","\\dstar","\\tstar"]
const COLORS = ["DodgerBlue", "OrangeRed", "Green",
  "Brown", "Black", "BlueViolet"]

const PlotContainer = Union{Plot,Gadfly.Compose.Context,Vector{Plot}} #holds graphs

const DATE_SYMBOLS =  Symbol("12/31/2010"), Symbol("3/31/2011"), Symbol("6/30/2011"),
  Symbol("9/30/2011"), Symbol("12/31/2011"), Symbol("3/31/2012"), Symbol("6/30/2012"),
  Symbol("9/30/2012"), Symbol("12/31/2012"), Symbol("3/31/2013"), Symbol("6/30/2013"),
  Symbol("9/30/2013"), Symbol("12/31/2013"), Symbol("3/31/2014"), Symbol("6/30/2014"),
  Symbol("9/30/2014"), Symbol("12/31/2014"), Symbol("3/30/2015"), Symbol("6/30/2015"),
  Symbol("9/30/2015"), Symbol("12/31/2015"), Symbol("3/30/2016"), Symbol("6/30/2016"),
  Symbol("9/30/2016"), Symbol("12/31/2016"), Symbol("3/30/2017"), Symbol("6/30/2017"),
  Symbol("9/30/2017"), Symbol("12/31/2017"), Symbol("3/31/2018")

const LAUNCH_DATE_SH = Date("11/17/2014",  DATE_FORMAT)
const LAUNCH_DATE_SZ = Date("12/5/2016",  DATE_FORMAT)
const FOCAL_DATES = (Date("6/30/2013",  DATE_FORMAT), Date("3/31/2018",  DATE_FORMAT))
const FOCAL_DATES_DD1 = (Date("6/30/2013",  DATE_FORMAT), Date("3/30/2016",  DATE_FORMAT))
const FOCAL_DATES_DD2 = (Date("6/30/2015",  DATE_FORMAT):Date("3/31/2018",  DATE_FORMAT))
const OUTCOMES = ["NORMALIZED_ACCRUALS_CF_METHOD"]

const MFloat64 = Union{Missing, Float64}
const MInt = Union{Missing, Int}
const MBool = Union{Missing, Bool}
const MSymbol = Union{Missing, Symbol}

n2s(x::T where T<:Real) = num2Str(x, 3, Ints=true)

#preprocess rate data
function preProcessChinaDSW()::Nothing
  wideDF::DataFrame = CSV.read("$DATA_PATH\\$CHINA_DS_NAME.csv", strings=:raw)

  stream::IOStream = open("$DATA_PATH\\$CHINA_DS_NAME_W.jls", "w")
  serialize(stream, wideDF)
  close(stream)
  return nothing
end

#to get the df
function getDF(pathStr::String, refreshData::Bool, preProcessFunc::Function)::DataFrame

  #re-process the data if needed or desired
  if !isfile(pathStr) || refreshData
    preProcessFunc()
  end

  #read the binary file
  stream::IOStream = open(pathStr)
  DF::DataFrame = deserialize(stream)
  close(stream)

  return DF
end

function summaryVTable(wideDFV::DataFrame,
    decimals::Int = DECIMALS)::String

  numColumns = 8

  #setup the table headers
  columnNames::Vector{Vector{String}} = [["Shanghai Volume \$\\times 10^{-12}\$", "Shenzhen Volume \$\\times 10^{-12}\$"], ["Shr Vol", "N.B. Buys",
      "N.B. Sells", "Total N.B.", "Shr Vol", "N.B. Buys", "N.B. Sells", "Total N.B."]]
  columnSizes::Vector{Vector{Int}} = [ones(Int, numColumns ÷ 4) .* 4, ones(Int, numColumns)]

  ####################write descriptive rows
  descRowNames::Vector{String} = ["Mean", "Median", "Std Deviation",
    "10\\% Quantile", "90\\% Quantile", "N>0 (Quarters)"]

  #need to print the descriptive rows
  descContent::Vector{Vector{String}} = ((i::Int)->Vector{String}(numColumns)).(1:length(descRowNames))

  #get and scale the summary variables
  series::Vector{Vector{MFloat64}} = [
      wideDFV[:SHCOMP_Index],
      wideDFV[:C1DBTO_Index],
      wideDFV[:C1DSTO_Index],
      wideDFV[:C1DBTO_Index] .+ wideDFV[:C1DSTO_Index],
      wideDFV[:SZCOMP_Index],
      wideDFV[:C2DBTO_Index],
      wideDFV[:C2DSTO_Index],
      wideDFV[:C2DBTO_Index] .+ wideDFV[:C2DSTO_Index]]

  #remove missing
  for i::Int ∈ 1:length(series)
    series[i] = series[i][((x::MFloat64)->!ismissing(x)).(series[i])]
    series[i] .= series[i] .* 10^-9.
  end

  for c::Int ∈ 1:numColumns
    descContent[1][c] = num2Str(mean(series[c]), decimals)
    descContent[2][c] = num2Str(median(series[c]), decimals)
    descContent[3][c] = num2Str(std(series[c]), decimals)
    descContent[4][c] = num2Str(quantile(Vector{Float64}(series[c]),0.1), decimals)
    descContent[5][c] = num2Str(quantile(Vector{Float64}(series[c]),0.9), decimals)
    descContent[6][c] = num2Str(length(series[c]), Ints=true)
  end

  summaryVTable::String = texTable( "Summary of Volume Data",
      """See Tex File""", #caption
      columnNames, #colNames
      Vector{String}(),#contentRowNames
      Vector{Matrix{String}}(), #content
      descRowNames, #descRowNames
      descContent, #descContent
      Vector{String}(),
      widthColNames = columnSizes,
      alignmentColNames = [["c", "c"], ["r", "r", "r", "r|", "r" , "r", "r", "r"]],
      lineSpacer = "\\\\"
    )

    writeTables2File([summaryVTable],
        HEADER_NAME, FOOTER_NAME, path=OUTPUT_PATH,
        outName = "trainVTables.tex")

  return summaryVTable
end

#computes policy outcomes for the different exchanges
function getExchangeVolumeData(longDFV::DataFrame)::DataFrame

  wideDFV::DataFrame = unstack(longDFV, :date, :security, :value)
  N::Int = size(wideDFV,1)

  #sum buy and sell orders

  wideDFV[:C1Vol] = Vector{MFloat64}(wideDFV[:C1DBTO_Index] .+ wideDFV[:C1DSTO_Index])
  wideDFV[:C2Vol] = Vector{MFloat64}(wideDFV[:C2DBTO_Index] .+ wideDFV[:C2DSTO_Index])

  wideDFV[:lC1Vol] = Vector{MFloat64}(wideDFV[:C1DBTO_Index] .- wideDFV[:C1DSTO_Index])
  wideDFV[:lC2Vol] = Vector{MFloat64}(wideDFV[:C2DBTO_Index] .- wideDFV[:C2DSTO_Index])

  #container for the policy variable
  wideDFV[:C1Effect] = Vector{MFloat64}(missings(N))
  wideDFV[:C2Effect] = Vector{MFloat64}(missings(N))

  wideDFV[:lC1Effect] = Vector{MFloat64}(missings(N))
  wideDFV[:lC2Effect] = Vector{MFloat64}(missings(N))

  sort!(wideDFV, :date)

  #this loop computes running totals of volue and sums them
  #the formula is connect period buys + sells / total period volume
  for i::Int ∈ VOL_PERIODS:N
    wideDFV[i, :C1Effect] = sum(skipmissing(wideDFV[(i-VOL_PERIODS+1):i,:C1Vol])) /
      sum(skipmissing(wideDFV[(i-VOL_PERIODS+1):i,:SHCOMP_Index]))

    wideDFV[i, :C2Effect] = sum(skipmissing(wideDFV[(i-VOL_PERIODS+1):i,:C2Vol])) /
      sum(skipmissing(wideDFV[(i-VOL_PERIODS+1):i,:SZCOMP_Index]))

    wideDFV[i, :lC1Effect] = sum(skipmissing(wideDFV[(i-VOL_PERIODS+1):i,:lC1Vol])) /
      sum(skipmissing(wideDFV[(i-VOL_PERIODS+1):i,:SHCOMP_Index]))

    wideDFV[i, :lC2Effect] = sum(skipmissing(wideDFV[(i-VOL_PERIODS+1):i,:lC2Vol])) /
      sum(skipmissing(wideDFV[(i-VOL_PERIODS+1):i,:SZCOMP_Index]))
  end

  return wideDFV
end

#runs any checks they are required on the wideDF post data-quality
function wideDiagnostics(wideDF::DataFrame)::Nothing

  #checks to see if we have cross-listings or duplicates
  countDF = by(wideDF, :security, nrow)
  NSecurities::Int = size(countDF,1)
  securities = countDF[:security]
  countIndex::Dict{Symbol, Int} = Dict(countDF[i,:security]=>countDF[i,:x1] for i::Int ∈ 1:size(countDF,1))

  #check for overloaded securities
  for i::Int ∈ 1:NSecurities
    if countIndex[securities[i]] > 1
      warn("Overloaded security: $(securities[i]) ($(countIndex[securities[i]])")
    end
  end

  return nothing
end

function graphTrends(longDF::DataFrame)::Nothing
  dateFormat::String = "yyyy-mm-dd"

  #get the names of all columns to be NOT unstacked
  #nonDateCols::Vector{Symbol} = setdiff(names(longDF), [variable, value])
  #wideDF::DataFrame = unstack(longDF, nonDateCols, variable, value)

  #set up the plotting environment
  Gadfly.push_theme(:default)
  plotNames::Vector{String} = Vector{String}()
  plots::Vector{PlotContainer} = Vector{PlotContainer}()

  #copy this to transform the data
  outDF::DataFrame = deepcopy(longDF[:,[:oldDate, :exchange, :value]])
  dropmissing!(outDF)
  aggDF::DataFrame = aggregate(outDF, [:exchange, :oldDate], mean)

  #showcols(aggDF)

  #D10::DateTime = Date(2010,12,31)
  #D11::DateTime = Date(2011,12,31)
  #D12::DateTime = Date(2012,12,31)
  D13::DateTime = Date(2013,12,31)
  D14::DateTime = Date(2014,12,31)
  D15::DateTime = Date(2015,12,31)
  D16::DateTime = Date(2016,12,31)
  D17::DateTime = Date(2017,12,31)
  majorDates::Vector{DateTime}= [#=D10, D11, D12, =#D13, D14, D15, D16, D17]
  #majorDates:::Vector{Int} = [2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017]

  push!(plotNames, "Parallel Trends")
  push!(plots, plot(aggDF, x=:oldDate, y=:value_mean, color=:exchange,
      Stat.xticks(majorDates), Geom.line, Guide.xlabel("Year"),
      Guide.ylabel("Average Firm Cashflow Accrual"),
      Scale.x_continuous(labels=x::DateTime->"$(Dates.year(x))"),
      Coord.Cartesian(xmin=Date(2013,6,30), xmax=D17), Guide.colorkey(title="Exchange:"),
      style(key_position=:bottom, guide_title_position=:left),
      xintercept=[LAUNCH_DATE_SH, LAUNCH_DATE_SZ], Geom.vline(style=[:solid, [1mm,1mm]])))

  for i::Int ∈ 1:length(plots)
    draw(SVG("$OUTPUT_PATH\\$(plotNames[i]).svg", 9inch, 7inch), plots[i])
    #display(plots[i])
  end

  return nothing

end

#conforms the industries
function conformIndustries!(longDF::DataFrame;
    variable::Symbol=:date, value::Symbol=:value,
    threshold::Int = MIN_SAMPLES_IN_INDUSTRY, newCol::Symbol = :gics,
    setMissing::Bool = SET_MISSING_IF_BELOW,
    runWideDiagnostics=true)::DataFrame

  dateFormat::String = "yyyy-mm-dd"

  #get the names of all columns to be NOT unstacked
  nonDateCols::Vector{Symbol} = setdiff(names(longDF), [variable, value])
  wideDF::DataFrame = unstack(longDF, nonDateCols, variable, value)

  #convert to a symbol
  wideDF[:,newCol] = Vector{MSymbol}((Symbol).(wideDF[:,:subIndustry]))

  #start at the lowest level of industry classification
  countDF::DataFrame = by(wideDF, newCol, nrow)

  countIndex::Dict{Symbol, Int} = Dict(countDF[i,newCol]=>countDF[i,:x1] for i::Int ∈ 1:size(countDF,1))

  #For each industry, examine the minimum count for subindustries. If its less than desirable
  #then conform the subindustries to industries
  by(wideDF, :industry) do df::SubDataFrame

    subGICS::Vector{Symbol} = unique(df[:,newCol])
    minCount = minimum(((s::Symbol)->countIndex[s]).(subGICS))

    #conform iff the min count is below the desired threshold
    if minCount < threshold
      df[:,newCol] = Symbol(df[1,:industry])
    end
  end

  #repeat for industry groups and industries
  countDF = by(wideDF, newCol, nrow)
  countIndex = Dict(countDF[i,newCol]=>countDF[i,:x1] for i::Int ∈ 1:size(countDF,1))
  by(wideDF, :industryGroup) do df::SubDataFrame

    subGICS::Vector{Symbol} = unique(df[:,newCol])
    minCount = minimum(((s::Symbol)->countIndex[s]).(subGICS))

    #conform iff the min count is below the desired threshold
    if minCount < threshold
      df[:,newCol] = Symbol(df[1,:industryGroup])
    end
  end

  #and again one more time to conform some industry groups to sectors
  countDF = by(wideDF, newCol, nrow)
  countIndex = Dict(countDF[i,newCol]=>countDF[i,:x1] for i::Int ∈ 1:size(countDF,1))
  by(wideDF, :sector) do df::SubDataFrame

    subGICS::Vector{Symbol} = unique(df[:,newCol])
    minCount = minimum(((s::Symbol)->countIndex[s]).(subGICS))

    #conform iff the min count is below the desired threshold
    if minCount < threshold
      df[:,newCol] = Symbol(df[1,:sector])
    end
  end

  #finally, if the threshold is sitll not met, set to missing (if the option is selected)
  if setMissing
    countDF = by(wideDF, newCol, nrow)
    countIndex = Dict(countDF[i,newCol]=>countDF[i,:x1] for i::Int ∈ 1:size(countDF,1))
    by(wideDF, :sector) do df::SubDataFrame

      subGICS::Vector{Symbol} = unique(df[:,newCol])
      minCount = minimum(((s::Symbol)->countIndex[s]).(subGICS))

      #conform iff the min count is below the desired threshold
      if minCount < threshold
        df[:,newCol] = missing
      end
    end
  end

  #reform the long dataframe
  if runWideDiagnostics
    wideDiagnostics(wideDF)
  end

  longDF = melt(wideDF, [nonDateCols; newCol])
  longDF[:,:date] = (s::Symbol->Date(string(s),dateFormat)).(longDF[:,:variable])
  delete!(longDF, :variable)
  longDF = longDF[!(ismissing).(longDF[:,:value]),:] # need to refilter missing variables
  return longDF
end

#conveience function to drop bad rows
function pruneLong!(longDF::DataFrame)::DataFrame

  #make a tracking column and store some convient quantities
  dateFormat::String = DATE_FORMAT
  N::Int = size(longDF,1)
  longDF[:,:keep] = Vector{Bool}(trues(N))
  numQuarters::Int =
      sum((s::Symbol->Date(string(s),dateFormat) ∈ minimum(FOCAL_DATES):maximum(
      FOCAL_DATES)).(DATE_SYMBOLS))

  println(unique(longDF[:date]))
  #now process the firm-level data

  by(longDF, :security) do df::SubDataFrame
    nSub::Int = size(df,1)

    #println("$nSub")
    #println("nSub: $nSub numQuarters:$numQuarters")
    #check for missing values, inconsistent lengths, and stocks switching exchanges
    df[(ismissing).(df[:,:value]) |
      (ismissing).(df[:exchange]) |
      (ismissing).(df[:marketCap]),:keep] = false

    #drop a stock if it has less than two years of data, if its too small,
    # or its listed on multiple exchanges
    if ((sum(df[:, :keep] .== true) < MIN_COMPANY_QUARTERS) ||
        (length(unique(df[:,:exchange])) ≠ 1) ||
        sum((!ismissing).(df[:marketCap])) == 0 || #require at least one marketcap
        (minimum(skipmissing(df[:marketCap])) < 5.) )

      df[:,:keep] = false

    end
  end

  longDF = longDF[longDF[:keep],:]
  ##finally need to check for  stocks listed multiple times

  by(longDF, :description) do dfName::SubDataFrame #for each company name

    by(dfName, :date) do dfUniqueDt #for each date
      if size(dfUniqueDt,1) > 1 #if there is a dup
        if KEEP_LARGEST_WEIGHT_IF_DUP #if true, we keep the largest weighted name
          maxWeight::Float64 = maximum(dfUniqueDt[:,:weight])
          dfUniqueDt[dfUniqueDt[:,:weight] .≠ maxWeight,:keep] = false
        else
          dfUniqueDt[:,:keep] = false
        end
      end
    end
  end

  longDF = longDF[longDF[:, :keep] .== true,:]

  #finally, do the winsorization
  keepSeries::Vector{Bool} = longDF[:,:keep]
  for s::Symbol ∈ [:value, :marketCap]
    series::Vector{MFloat64} = longDF[:,s]

    #μ::Float64 = mean(skipmissing(series))
    #σ::Float64 = std(skipmissing(series))
    lower::Float64 = ACCEPTABLE_RANGE[s][1]
    upper::Float64 = ACCEPTABLE_RANGE[s][2]
    for i::Int ∈ 1:length(series)
      if !ismissing(series[i])
        if series[i] < lower
          if DROP_OUTSIDE_RANGE
            keepSeries[i] = false
          else
            series[i] = lower
          end
        elseif series[i] > upper
          if DROP_OUTSIDE_RANGE
            keepSeries[i] = false
          else
            series[i] = upper
          end
        end
      end
    end

    longDF[:, s] = series
  end

  #use an anonymous function to execute the SZCECOM cutoff
  longDF[(longDF[:, :index] .== :SZSECOM) & ((f::MFloat64->
    ismissing(f) ?  false : (f < MIN_SZSECOM_MARKETCAP)
    ).(longDF[:, :marketCap])), :keep] = false

  longDF[:,:keep] = keepSeries
  longDF = longDF[longDF[:keep] .== true,:]
  delete!(longDF, :keep)

  return longDF
end

#processing activities for the long dataframe
function processLong!(wideDF::DataFrame)::DataFrame
  dateFormat::String = DATE_FORMAT

  #use all columns except the dates as ID columns
  #showcols(wideDF)
  nonDateCols::Vector{Symbol} = setdiff(names(wideDF), DATE_SYMBOLS)
  longDF::DataFrame = melt(wideDF, nonDateCols)

  #convert the newly melted columns to dates
  longDF[:,:date] = (s::Symbol->Date(string(s),dateFormat)).(longDF[:,:variable])
  delete!(longDF, :variable)
  longDF =
    longDF[(d::Date->d ∈ minimum(FOCAL_DATES):maximum(FOCAL_DATES)).(longDF[:,:date]),:] #filter to the correct dates
  #println("names: $(names(longDF))")
  println(unique(longDF[:date]))
  #clean up names
  longDF[:,:security] =
      (s::String->Symbol(replace(strip(s), " ","_"))).(longDF[:,:security])

  #drop and collect non-outcome rows. Start with a flag for outcome variables
  longDF[:,:Outcome] =
      (s::CategoricalString->s ∈ OUTCOMES ? true : false).(longDF[:command])
  longDFV::DataFrame =
      longDF[.!(longDF[:Outcome]),[:date, :security, :value]]
  longDF = longDF[longDF[:Outcome],:]
  delete!(longDF, :Outcome) #no longer need this

  #flag for deleting rows with bad data
  longDF = longDF[(!).((ismissing).(longDF[:exchange])), :]
  longDF[:,:exchange] = Vector{MSymbol}((Symbol).(longDF[:,:exchange]))

  #delete marked rows
  println("Pruning Step 1. Starting rows: $(size(longDF,1))")
  longDF = pruneLong!(longDF)
  println("Pruning complete. Rows: $(size(longDF,1))")
  N::Int = size(longDF,1) #udpate size

  #conform industry codes so each industry gets some variation
  #a broader classificaiton is used for the interacted version
  longDF = conformIndustries!(longDF, runWideDiagnostics=false)
  #longDF = conformIndustries!(longDF, threshold=MIN_SAMPLES_INTERACTED, newCol = :gicsInteracted)

  N = size(longDF,1) #udpate size
  #get wide data for the input, and join to the long DF
  wideDFV::DataFrame = getExchangeVolumeData(longDFV)
  serialize(open("$DATA_PATH\\$(CHINA_DS_NAME)_VW.jls", "w"), wideDFV) #save a copy for later
  longDF = join(longDF,
      wideDFV[[:date, :C1Effect, :C2Effect, :lC1Effect, :lC2Effect]],
      on = :date, kind=:inner)

  summaryVTable(wideDFV)
  #now assign the treatment variables
  #begin by handling missing values and allocating space
  exchanges::Vector{MSymbol} = longDF[:exchange]

  #work with vectors for performance reasons
  C1Effect::Vector{MFloat64} = (x::MFloat64->ismissing(x) ? 0.0 : Float64(x)).(longDF[:, :C1Effect])
  C2Effect::Vector{MFloat64} = (x::MFloat64->ismissing(x) ? 0.0 : Float64(x)).(longDF[:, :C2Effect])
  CEffect::Vector{MFloat64} = Vector{MFloat64}(missings(N))
  lC1Effect::Vector{MFloat64} = (x::MFloat64->ismissing(x) ? 0.0 : Float64(x)).(longDF[:, :lC1Effect])
  lC2Effect::Vector{MFloat64} = (x::MFloat64->ismissing(x) ? 0.0 : Float64(x)).(longDF[:, :lC2Effect])
  lCEffect::Vector{MFloat64} = Vector{MFloat64}(missings(N))
  postTreatment = Vector{MInt}(missings(N))

  for i::Int ∈ 1:N
    if exchanges[i]==:Shanghai
      CEffect[i] = C1Effect[i]
      lCEffect[i] = lC1Effect[i]
    elseif exchanges[i]==:Shenzhen
      CEffect[i] = C2Effect[i]
      lCEffect[i] = lC2Effect[i]
    end

    postTreatment[i] = MInt(abs(CEffect[i]) > 0.0)
  end

  longDF[:, :CEffect] = CEffect #focal policy variable
  longDF[:, :lCEffect] = lCEffect #focal policy variable
  longDF[:, :postTreatment] = postTreatment
  longDF[:, :gicsXDate] =
    Vector{MSymbol}(((s1::MSymbol,s2::Date)->ismissing(s1) ? missing : Symbol(
    s1,s2)).(longDF[:, :gics],longDF[:, :date]))

  longDF[:,:lMarketCap] = Vector{MFloat64}((log).(longDF[:,:marketCap]))
  longDF[:, :oldDate] = Vector{Date}(deepcopy(longDF[:,:date]))

  longDF[:, :DD1] = (d::Date->MInt(
      d ∈ minimum(FOCAL_DATES_DD1):maximum(FOCAL_DATES_DD1))).(longDF[:, :date])
  longDF[:, :DD2] = (d::Date->MInt(
      d ∈ minimum(FOCAL_DATES_DD2):maximum(FOCAL_DATES_DD2))).(longDF[:, :date])

  categorical!(longDF, :security)
  categorical!(longDF, :exchange)
  categorical!(longDF, :gics)
  #  categorical!(longDF, :gicsInteracted)
  categorical!(longDF, :date)
  categorical!(longDF, :gicsXDate)

  #println("unique date= $(length(unique(longDF[:date])))")
  #println("unique industry= $(length(unique(longDF[:gics])))")
  #println("unique interacted = $(length(unique(longDF[:gicsXDate])))")

  println("\nLength: $(size(longDF,1)), Length Unique: $(size(unique(longDF[:security]),1))")
  return longDF
end

function processDF(refreshData::Bool = true)
  wideDF::DataFrame = getDF("$DATA_PATH\\$CHINA_DS_NAME_W.jls", refreshData, preProcessChinaDSW)
  longDF::DataFrame = processLong!(wideDF)

  stream::IOStream = open("$DATA_PATH\\$CHINA_DS_NAME_L.jls", "w")
  serialize(stream, longDF)
  close(stream)
end

function discrete(longDF::DataFrame ;
    YSym::Symbol = :value, focalSpec::String = "postTreatment",
    title::String = "Discrete Specifications",
    focalName=:integrated)::String

    decimals::Int = DECIMALS

  #showcols(longDF[:,[:lMarketCap, :value, :gics, :ey, :exchange, :postTreatment, :date]])

  #println("got here")
  #setup the containers
  YSpecs::Vector{Symbol} = Vector{Symbol}() # holds the LHS
  XSpecs::Vector{CTExpr} = Vector{CTExpr}() #holds the RHS
  XNames::Vector{Vector{Symbol}} = Vector{Vector{Symbol}}() # holds the names of the x variables
  withinSpec::Vector{CTExpr} = Vector{CTExpr}() #holds the within transformation variable
  clusteredSpec::Vector{CTExpr} = Vector{CTExpr}() #holds the clustering variable
  models::Vector{CTLM} = Vector{CTLM}() #holds the regression models

  #showall(describe(longDF))
  #println(longDF[1:5, [:security, :postTreatment, :exchange, :date]])

  #############specifications here

  push!(XSpecs, parse("$focalSpec + date"))
  push!(XNames, [:intercept; focalName])
  push!(YSpecs, YSym)
  push!(withinSpec, :exchange)
  push!(clusteredSpec, :gics) #tag1

  push!(XSpecs, parse("$focalSpec + lMarketCap + date"))
  push!(XNames, [:intercept; focalName; :lMarketCap])
  push!(YSpecs, YSym)
  push!(withinSpec, :exchange)
  push!(clusteredSpec, :gics) #tag1

  push!(XSpecs, parse("$focalSpec + date"))
  push!(XNames, [:intercept; focalName])
  push!(YSpecs, YSym)
  push!(withinSpec, :security)
  push!(clusteredSpec, :security)

  push!(XSpecs, parse("$focalSpec + exchange + date"))
  push!(XNames, [:intercept; focalName])
  push!(YSpecs, YSym)
  push!(withinSpec, :gics)
  push!(clusteredSpec, :gics)

  push!(XSpecs, parse("$focalSpec + lMarketCap + exchange + date"))
  push!(XNames, [:intercept; focalName; :lMarketCap ])
  push!(YSpecs, YSym)
  push!(withinSpec, :gics)
  push!(clusteredSpec, :gics)

  push!(XSpecs, parse("$focalSpec + exchange"))
  push!(XNames, [:intercept; focalName])
  push!(YSpecs, YSym)
  push!(withinSpec, :gicsXDate)
  push!(clusteredSpec, :gics)

  push!(XSpecs, parse("$focalSpec + lMarketCap + exchange"))
  push!(XNames, [:intercept; focalName; :lMarketCap])
  push!(YSpecs, YSym)
  push!(withinSpec, :gicsXDate)
  push!(clusteredSpec, :gics)

  ###################end new regression specs

  NSpecs::Int = length(XSpecs)
  columnNames::Vector{Vector{String}} = [((i::Int)->"($i)").(1:NSpecs)]

  println("Running $title with focal var $focalSpec")
  for i::Int ∈ 1:NSpecs
      #println("spec: $(XSpecs[i])")
      push!(models, CTLM(longDF, XSpecs[i],  YSpecs[i], withinSym = withinSpec[i],
          clusterSym = clusteredSpec[i], XNames=XNames[i], YName = YSpecs[i]))
  end

  ####################write descriptive rows
  descRowNames::Vector{String} = ["Date F.E.",
      "Firm F.E.", "Exchange F.E.", "Industry F.E.", "Industry-Date F.E",
      "\$R^{2}\$ (within)", "N (firm-quarters)"]

  #need to print the descriptive rows
  descContent::Vector{Vector{String}} =
      ((i::Int)->Vector{String}(length(models))).(1:length(descRowNames))

  #this builds the descriptive rows. There is Probably a better way to do this,
  #but its fairly specific to the project.
  for i ∈ 1:length(XSpecs)
    descContent[1][i] = "$(withinSpec[i]==:date ||
        contains(string(XSpecs[i]), "date") ? "X" : "")"
    descContent[2][i] = "$(withinSpec[i]==:security ||
        contains(string(XSpecs[i]), "security") ? "X" : "")"
    descContent[3][i] = "$(withinSpec[i]==:exchange ||
        contains(string(XSpecs[i]), "exchange") ? "X" : "")"
    descContent[4][i] = "$(withinSpec[i]==:gics ? "X" : "")"
    descContent[5][i] = "$(withinSpec[i]==:gicsXDate ? "X" : "")"
    descContent[6][i] =
    "$(num2Str(getR²(models[i], adjusted=false),  decimals))"
    descContent[7][i] = "$(models[i].N)"
  end

  tableText::String = texTable(models,
      getClustered!#=getNeweyWestFunc(5) (lm::CTLM)->getNeweyWestSlow(lm, 5)=#,
      [focalName, :lMarketCap],
      titleCaption = title,
      colNames = columnNames,
      contentRowNames = ["Integrated", "lMarketCap"],
      descRowNames = descRowNames,
      descContent = descContent,
      decimalDigits = 4,
      stars=true,
      starStrings = OVERRIDE_STAR_STRINGS,
      #clearMem = USE_AGGRESSIVE_GC,
      caption = "to be written")

  return tableText
end

function policy(longDF::DataFrame ;
    YSym::Symbol = :value, focalSpec::String = "CEffect",
    title::String = "Policy Specifications",
    focalName=:IMeasure)::String

    decimals::Int = DECIMALS

  #showcols(longDF[:,[:lMarketCap, :value, :gics, :ey, :exchange, :postTreatment, :date]])

  #println("got here")
  #setup the containers
  YSpecs::Vector{Symbol} = Vector{Symbol}() # holds the LHS
  XSpecs::Vector{CTExpr} = Vector{CTExpr}() #holds the RHS
  XNames::Vector{Vector{Symbol}} = Vector{Vector{Symbol}}() # holds the names of the x variables
  withinSpec::Vector{CTExpr} = Vector{CTExpr}() #holds the within transformation variable
  clusteredSpec::Vector{CTExpr} = Vector{CTExpr}() #holds the clustering variable
  models::Vector{CTLM} = Vector{CTLM}() #holds the regression models

  #showall(describe(longDF))
  #println(longDF[1:5, [:security, :postTreatment, :exchange, :date]])

  #############specifications here

  push!(XSpecs, parse("$focalSpec + date"))
  push!(XNames, [:intercept; focalName])
  push!(YSpecs, YSym)
  push!(withinSpec, :exchange)
  push!(clusteredSpec, :gics) #tag1

  push!(XSpecs, parse("$focalSpec + lMarketCap + date"))
  push!(XNames, [:intercept; focalName; :lMarketCap])
  push!(YSpecs, YSym)
  push!(withinSpec, :exchange)
  push!(clusteredSpec, :gics) #tag1

  push!(XSpecs, parse("$focalSpec + date"))
  push!(XNames, [:intercept; focalName])
  push!(YSpecs, YSym)
  push!(withinSpec, :security)
  push!(clusteredSpec, :security)

  push!(XSpecs, parse("$focalSpec + exchange + date"))
  push!(XNames, [:intercept; focalName])
  push!(YSpecs, YSym)
  push!(withinSpec, :gics)
  push!(clusteredSpec, :gics)

  push!(XSpecs, parse("$focalSpec + lMarketCap + exchange + date"))
  push!(XNames, [:intercept; focalName; :lMarketCap])
  push!(YSpecs, YSym)
  push!(withinSpec, :gics)
  push!(clusteredSpec, :gics)

  push!(XSpecs, parse("$focalSpec + exchange"))
  push!(XNames, [:intercept; focalName])
  push!(YSpecs, YSym)
  push!(withinSpec, :gicsXDate)
  push!(clusteredSpec, :gics)

  push!(XSpecs, parse("$focalSpec + lMarketCap + exchange"))
  push!(XNames, [:intercept; focalName; :lMarketCap])
  push!(YSpecs, YSym)
  push!(withinSpec, :gicsXDate)
  push!(clusteredSpec, :gics)

  ###################end new regression specs

  NSpecs::Int = length(XSpecs)
  columnNames::Vector{Vector{String}} = [((i::Int)->"($i)").(1:NSpecs)]

  println("Running $title with focal var $focalSpec")
  for i::Int ∈ 1:NSpecs
      #println("spec: $(XSpecs[i])")
      push!(models, CTLM(longDF, XSpecs[i],  YSpecs[i], withinSym = withinSpec[i],
          clusterSym = clusteredSpec[i], XNames=XNames[i], YName = YSpecs[i]))
  end

  ####################write descriptive rows
  descRowNames::Vector{String} = ["Date F.E.",
      "Firm F.E.", "Exchange F.E.", "Industry F.E.", "Industry-Date F.E",
      "\$R^{2}\$ (within)", "N (firm-quarters)"]

  #need to print the descriptive rows
  descContent::Vector{Vector{String}} =
      ((i::Int)->Vector{String}(length(models))).(1:length(descRowNames))

  #this builds the descriptive rows. There is Probably a better way to do this,
  #but its fairly specific to the project.
  for i ∈ 1:length(XSpecs)
    descContent[1][i] = "$(withinSpec[i]==:date ||
        contains(string(XSpecs[i]), "date") ? "X" : "")"
    descContent[2][i] = "$(withinSpec[i]==:security ||
        contains(string(XSpecs[i]), "security") ? "X" : "")"
    descContent[3][i] = "$(withinSpec[i]==:exchange ||
        contains(string(XSpecs[i]), "exchange") ? "X" : "")"
    descContent[4][i] = "$(withinSpec[i]==:gics ? "X" : "")"
    descContent[5][i] = "$(withinSpec[i]==:gicsXDate ? "X" : "")"
    descContent[6][i] =
    "$(num2Str(getR²(models[i], adjusted=false),  decimals))"
    descContent[7][i] = "$(models[i].N)"
  end

  tableText::String = texTable(models,
      getClustered!#=getNeweyWestFunc(5) (lm::CTLM)->getNeweyWestSlow(lm, 5)=#,
      [focalName, :lMarketCap],
      titleCaption = title,
      colNames = columnNames,
      contentRowNames = ["$focalName", "lMarketCap"],
      descRowNames = descRowNames,
      descContent = descContent,
      decimalDigits = 4,
      stars=true,
      starStrings = OVERRIDE_STAR_STRINGS,
      #clearMem = USE_AGGRESSIVE_GC,
      caption = "to be written")

  return tableText
end

function summary(longDF::DataFrame, decimals::Int = DECIMALS)::String
  numColumns = 6

  #setup the table headers
  columnNames::Vector{Vector{String}} = [["Accrual Ratio", "Market Cap \$\\times 10^{-9}\$"],
      ["Shanghai", "Shenzhen", "Total", "Shanghai", "Shenzhen", "Total"]]
  columnSizes::Vector{Vector{Int}} = [ones(Int, numColumns ÷ 3) .* 3, ones(Int, numColumns)]

  #make views for conveneince
  longDFSH::SubDataFrame = view(longDF, longDF[:,:exchange] .== :Shanghai)
  longDFSZ::SubDataFrame = view(longDF, longDF[:,:exchange] .== :Shenzhen)

  ####################write descriptive rows
  descRowNames::Vector{String} = ["Mean", "Median", "Std Deviation",
    "10\\% Quantile", "90\\% Quantile", "Unique Firms", "N (firm-quarters)"]

  #need to print the descriptive rows
  descContent::Vector{Vector{String}} = ((i::Int)->Vector{String}(numColumns)).(1:length(descRowNames))

  #get and scale the summary variables
  series::Vector{Vector{MFloat64}} = [
      longDFSH[:value], longDFSZ[:value], longDF[:value],
        longDFSH[:marketCap].*10. ^ -9.,
        longDFSZ[:marketCap].*10. ^ -9.,
        longDF[:marketCap].*10. ^ -9.]

  #remove missing
  for i::Int ∈ 1:length(series)
    series[i] = series[i][((x::MFloat64)->!ismissing(x)).(series[i])]
  end

  #get the unique number of firms
  descContent[6] = [num2Str(length(unique(longDFSH[:security])), Ints=true),
      num2Str(length(unique(longDFSZ[:security])), Ints=true),
      num2Str(length(unique(longDF[:security])), Ints=true), "", "", ""]

  for c::Int ∈ 1:numColumns
    descContent[1][c] = num2Str(mean(series[c]), decimals)
    descContent[2][c] = num2Str(median(series[c]), decimals)
    descContent[3][c] = num2Str(std(series[c]), decimals)
    descContent[4][c] = num2Str(quantile(Vector{Float64}(series[c]),0.1), decimals)
    descContent[5][c] = num2Str(quantile(Vector{Float64}(series[c]),0.9), decimals)
    descContent[7][c] = num2Str(length(series[c]), Ints=true)
  end

  summaryTable::String = texTable( "Summary of Data",
      """See Tex File""", #caption
      columnNames, #colNames
      Vector{String}(),#contentRowNames
      Vector{Matrix{String}}(), #content
      descRowNames, #descRowNames
      descContent, #descContent
      Vector{String}(),
      widthColNames = columnSizes,
      alignmentColNames = [["c", "c"], ["r" for i::Int ∈ 1:numColumns]],
      lineSpacer = "\\\\"
    )

  return summaryTable
end

function runRegressions(longDF::DataFrame; runDiscrete::Bool = true, runPolicy::Bool = true,
    runSummary::Bool = true)

  longDF1::DataFrame = longDF[longDF[:DD1] .== 1, :]
  longDF2::DataFrame = longDF[longDF[:DD2] .== 1, :]
  #run the regression with the discrete specification
  discreteTable1::String = runDiscrete ? discrete(longDF1, title="Shanghai Interconnect Discrete") : ""
  discreteTable2::String = runDiscrete ? discrete(longDF2, title="Shenzhen Interconnect Discrete") : ""
  discreteTable::String = runDiscrete ? discrete(longDF, title="Combined Interconnect Discrete") : ""

  policyTableVol1::String = runPolicy ? policy(longDF1, focalSpec="CEffect",
    title="Shanghai Interconnect Volume Policy", focalName=:IMeasureVol) : ""
  policyTableVol2::String = runPolicy ? policy(longDF2, focalSpec="CEffect",
    title="Shenzhen Interconnect Volume Policy", focalName=:IMeasureVol) : ""
  policyTableVol::String = runPolicy ? policy(longDF, focalSpec="CEffect",
    title="Combined Interconnect Volume Policy", focalName=:IMeasureVol) : ""

  policyTableNet1::String = runPolicy ? policy(longDF1, focalSpec="lCEffect",
    title="Shanghai Interconnect Net Flow Policy", focalName=:IMeasureNet) : ""
  policyTableNet2::String = runPolicy ? policy(longDF2, focalSpec="lCEffect",
    title="Shenzhen Interconnect Net Flow Policy", focalName=:IMeasureNet) : ""
  policyTableNet::String = runPolicy ? policy(longDF, focalSpec="lCEffect",
    title="Combined Interconnect Net Flow Policy", focalName=:IMeasureNet) : ""

  #calculate the summary statistics
  summaryTable::String = runSummary ? summary(longDF) : ""
  graphTrends(longDF)

  #write the processed data to a CSV for any tests in R
  CSV.write("$R_TEST_PATH\\longDF.csv", longDF)

  #print tables
  writeTables2File([summaryTable,
      discreteTable1, discreteTable2, discreteTable,
      policyTableVol1, policyTableVol2, policyTableVol,
      policyTableNet1, policyTableNet2, policyTableNet],
      HEADER_NAME, FOOTER_NAME, path=OUTPUT_PATH,
      outName = "trainTables.tex")
end
#main script
function mainScript(;refreshData::Bool = true, processData::Bool =true,
    runDiscrete::Bool = true, runPolicy::Bool = true, runSummary::Bool = false)::Nothing

  longDF::DataFrame = getDF(
      "$DATA_PATH\\$CHINA_DS_NAME_L.jls", processData, ()->processDF(refreshData))

  # run the regressions
  runRegressions(longDF, runDiscrete = runDiscrete, runPolicy = runPolicy, runSummary = runSummary)

  return nothing
end

# Uncomment to run, commented to print the code appendix

@time begin #run main
  mainScript(refreshData = true, processData=true, runDiscrete=true,
      runPolicy=true, runSummary=true)
end

end #end module

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
"""
    function stepBoris(x,p,efield,bfield,q,m,dt)

Performs one complete time step of the Boris algorithm for relativistic charged particles in an electromagnetic field. This version is from Ripperda2017
"""
function stepBoris(x, p, efield, bfield, q, m, dt)
    qdt2m = q * dt / (2.0 * m)

    un = p / m
    γn = γFromP(p, m)
    xnph = x + un * dt / (2.0 * γn)

    eph = efield(xnph)
    um = un + qdt2m * eph

    γm = γFromP(m * um, m)
    t = (qdt2m / γm) * bfield(xnph)
    s = (2.0 / (1 + t ⋅ t)) * t

    up = um + cross((um + cross(um, t)), s)

    unp1 = up + eph * q * dt / (2.0 * m)
    pnew = m * unp1

    γnp1 = γFromP(pnew, m)
    xnew = xnph + (dt / (2.0 * γnp1)) * unp1
    return ArrayPartition(xnew, pnew)
end
export stepBoris

"""
    fBoris!(u::ArrayPartition, params::Dict, t)

Boris stepper wrapper for the discrete solver.
Use with single particle coordinates.

Input:
    - u{ArrayPartition{Number}} ≡ ArrayPartition([position],[momentum])
    - params{Dict} ≡ Dict[
            :t => {Number} current time value,
            :efield => {Function} returns electric field at [x,y,z]
            :bfield => {Function} returns magnetic flux density at [x,y,z]
            :q      => {Number} charge of particle in SI
            :m      => {Number} mass of particle in SI
        ]
    - t{Number} ≡ time at next time step

Returns:
    Particle coordinates at time t in the form
    {ArrayPartition{Number}}([[x,y,z],[px,py,pz]]).
"""
function fBoris!(u::ArrayPartition, params::Dict, t)
    dt = t - params[:t]
    params[:t] = t
    stepBoris(
        u.x[1],
        u.x[2],
        params[:efield],
        params[:bfield],
        params[:q],
        params[:m],
        dt
    )
end

"""
    fBoris!(u::Snapshot, params::Dict, t)

Boris stepper wrapper for the discrete solver.
Use with Snapshot branch of ParticleTree.

Input:
    - u{Snapshot}  ≡ level 2 ParticleTree branch
    - params{Dict} ≡ Dict[
            :t      => {Number} current time value,
            :efield => {Function} returns electric field at [x,y,z]
            :bfield => {Function} returns magnetic flux density at [x,y,z]
        ]
    - t{Number}    ≡ time at next time step

Returns:
    Snapshot at time t
"""
function fBoris!(u::Snapshot, params::Dict, t)
    dt = t - params[:t]
    params[:t] = t
    snapOut = deepcopy(u)
    for species in snapOut.nodes
        q = usistrip(charge(species.id))
        m = usistrip(mass(species.id))
        map(species.nodes) do coords
            stepCoords = stepBoris(
                    coords.values.x[1],
                    coords.values.x[2],
                    params[:efield],
                    params[:bfield],
                    q,
                    m,
                    dt
                )
            coords.values .= stepCoords
        end
    end
    return snapOut
end

export fBoris!

# dictionary of all relativistic integrators availble
rintegrator = Dict(
                :boris => fBoris!
)

export rintegrator


text:
function fractional_trophic_level(N::T) where {T<:UnipartiteNetwork}
  Y = nodiagonal(N)
  producers = keys(filter(spedeg -> spedeg.second == 0, degree(Y; dims=1)))
  sp = shortest_path(Y)
  prod_id = findall(isequal(0), vec(sum(sp; dims=2)))
  return Dict(zip(species(Y; dims=1), maximum(sp[:,prod_id]; dims=2).+1))
end

function trophic_level(N::T) where {T<:UnipartiteNetwork}
  TL = fractional_trophic_level(N)
  Y = nodiagonal(N)
  D = zeros(Float64, size(Y.A))
  ko = degree_out(Y)

  # inner loop to avoid dealing with primary producers
  for i in 1:richness(Y)
    if ko[species(Y)[i]] > 0.0
      D[i,:] = Y[i,:]./ko[species(Y)[i]]
    end
  end

  # return
  return Dict(zip(species(N), 1 .+ D * [TL[s] for s in species(Y)]))
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Test various optimistic bounds for GINI impurity, esp on split groups.
"""
module GINIOptimTest

export gini_optim_test, gini_optim_loop

import Compat.ASCIIString
using Distributions #mode
using Reexport
@reexport using DataFrames
using RLESUtils, MathUtils

function gini_optim_loop(ntrials=100)
  fs = [:evenly, :shortest, :longest, :lowest, :highest,
                       :weighted_lowest, :weighted_highest, :max_p, :nmax_p, :bound1, :bound1multi]
  results = DataFrame([ASCIIString, ASCIIString, Int64, Int64, Float64, Int64, Float64,
                       Float64, Float64, Float64, Float64, Float64, Float64, Float64,
                       Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64],
                      [:v1, :v2, :N, :n1, :g1, :n2, :g2, :max_p1, :max_p2, :nmax_p1, :nmax_p2,
                       fs...], 0)
  for i = 1:ntrials
    n1 = rand(1:10)
    v1 = rand(1:3, n1)
    n2 = rand(1:10)
    v2 = rand(1:3, n2)
    N = rand(1:20)
    trial = loopinner(v1, v2, N)
    push!(results, trial)
  end

  for sym in fs
    results[Symbol(sym, :_ismin)] = map(eachrow(results[:, fs])) do r
      a = convert(Array, r)
      Int64(all(r[sym] .<= a))
    end
  end
  results
end

function loopinner(v1, v2, N)
  n1 = length(v1)
  n2 = length(v2)
  mp1 = maximum(counts(v1) / n1)
  mp2 = maximum(counts(v2) / n2)
  [string(v1), string(v2), N,
   n1, gini_impurity(v1),
   n2, gini_impurity(v2),
   mp1, mp2, n1 * mp1, n2 * mp2,
   gini_impurity(fillevenly(v1, v2, N)...),
   gini_impurity(fill2short(v1, v2, N)...),
   gini_impurity(fill2long(v1, v2, N)...),
   gini_impurity(fill2lowest(v1, v2, N)...),
   gini_impurity(fill2highest(v1, v2, N)...),
   gini_impurity(fill2wlowest(v1, v2, N)...),
   gini_impurity(fill2whighest(v1, v2, N)...),
   gini_impurity(fill2maxp(v1, v2, N)...),
   gini_impurity(fill2nmaxp(v1, v2, N)...),
   gini_impurity(fill2bound1(v1, v2, N)...),
   gini_impurity(fill2bound1multi(v1, v2, N)...)]
end

function gini_optim_test(v1=[1,1,2], v2=[1,1,1,2], N=5)
  @show v1
  @show v2
  @show N
  @show gini_impurity(v1, v2)
  @show gini_impurity(fillevenly(v1, v2, N)...)
  @show gini_impurity(fill2short(v1, v2, N)...)
  @show gini_impurity(fill2long(v1, v2, N)...)
  @show gini_impurity(fill2lowest(v1, v2, N)...)
  @show gini_impurity(fill2highest(v1, v2, N)...)
end

function fillevenly(v1, v2, N)
  a = vcat(v1, fill(mode(v1), ceil(Int, N / 2)))
  b = vcat(v2, fill(mode(v2), floor(Int, N / 2)))
  return a, b
end

function fill2short(v1, v2, N)
  if length(v1) < length(v2)
    v = vcat(v1, fill(mode(v1), N))
    return v, v2
  else
    v = vcat(v2, fill(mode(v2), N))
    return v1, v
  end
end

function fill2long(v1, v2, N)
  if length(v1) > length(v2)
    v = vcat(v1, fill(mode(v1), N))
    return v, v2
  else
    v = vcat(v2, fill(mode(v2), N))
    return v1, v
  end
end

function fill2lowest(v1, v2, N)
  if gini_impurity(v1) < gini_impurity(v2)
    v = vcat(v1, fill(mode(v1), N))
    return v, v2
  else
    v = vcat(v2, fill(mode(v2), N))
    return v1, v
  end
end

function fill2highest(v1, v2, N)
  if gini_impurity(v1) > gini_impurity(v2)
    v = vcat(v1, fill(mode(v1), N))
    return v, v2
  else
    v = vcat(v2, fill(mode(v2), N))
    return v1, v
  end
end

function fill2wlowest(v1, v2, N)
  if length(v1) * gini_impurity(v1) < length(v2) * gini_impurity(v2)
    v = vcat(v1, fill(mode(v1), N))
    return v, v2
  else
    v = vcat(v2, fill(mode(v2), N))
    return v1, v
  end
end

function fill2whighest(v1, v2, N)
  if length(v1) * gini_impurity(v1) > length(v2) * gini_impurity(v2)
    v = vcat(v1, fill(mode(v1), N))
    return v, v2
  else
    v = vcat(v2, fill(mode(v2), N))
    return v1, v
  end
end

function fill2maxp(v1, v2, N)
  if maximum(counts(v1)) > maximum(counts(v2))
    v = vcat(v1, fill(mode(v1), N))
    return v, v2
  else
    v = vcat(v2, fill(mode(v2), N))
    return v1, v
  end
end

function fill2nmaxp(v1, v2, N)
  if length(v1) * maximum(counts(v1)) > length(v2) * maximum(counts(v2))
    v = vcat(v1, fill(mode(v1), N))
    return v, v2
  else
    v = vcat(v2, fill(mode(v2), N))
    return v1, v
  end
end

function fill2bound1(v1, v2, N)
  c1 = counts(v1)
  c2 = counts(v2)
  n1 = sum(c1)
  n2 = sum(c2)
  n1t = maximum(c1)
  n2t = maximum(c2)

  n1f = n1 - n1t
  n2f = n2 - n2t

  b1 = (n1f^2) / (n1 * (n1 + N))
  b2 = (n2f^2) / (n2 * (n2 + N))

  if b1 < b2
    v = vcat(v1, fill(mode(v1), N))
    return v, v2
  else
    v = vcat(v2, fill(mode(v2), N))
    return v1, v
  end
end

function fill2bound1multi(v1, v2, N)
  c1 = counts(v1)
  c2 = counts(v2)
  n1 = sum(c1)
  n2 = sum(c2)

  i1 = indmax(c1)
  i2 = indmax(c2)
  c1[i1] = 0
  c2[i2] = 0

  b1 = (sum(c1)^2 + sumabs2(c1)) / (n1 * (n1 + N))
  b2 = (sum(c2)^2 + sumabs2(c2)) / (n2 * (n2 + N))

  if b1 < b2
    v = vcat(v1, fill(mode(v1), N))
    return v, v2
  else
    v = vcat(v2, fill(mode(v2), N))
    return v1, v
  end
end

end #module


text:
module TotalVariation

using BenchmarkTools
using StructuredOptimization
using AbstractOperators
using Images
using ImageView
using TestImages

function set_up()
	srand(123)
	img = testimage("cameraman")

	X = convert(Array{Float64},img) # convert image to array
	Xt = X.+sqrt(0.006*vecnorm(X,Inf))*randn(size(X)) # add noise
	Xt[Xt .< 0] .= 0. #make sure pixels are in range
	Xt[Xt .> 1] .= 1.

	V = Variation(size(X)) # linear mapping operator
	lambda = 0.07          # level of regularization

	Y = Variable(size(V,1)...) # dual variables
	return V, Y, Xt, X, lambda
end

function run_demo()
	Lf  = 8                # Lipschitz constant
	slv = PANOC(tol = 1e-3, gamma = 1/Lf, adaptive = false)
	setup = set_up()
	@time solve_problem!(slv,setup...)
	return setup
end

function solve_problem!(slv,V, Y, Xt, X, lambda)
	it,  = @minimize ls(-V'*Y+Xt)+conj(lambda*norm(Y,2,1,2)) with slv
	return it
end

function benchmark(;verb = 0, samples = 5, seconds = 100, tol = 1e-3, maxit = 50000 )

	suite = BenchmarkGroup()

	solvers = ["ZeroFPR",
               "PANOC",
               "FPG",
               "PG"]
	slv_opt = ["(verbose = $verb, tol = $tol, gamma = 1/8, maxit = $maxit)",
               "(verbose = $verb, tol = $tol, gamma = 1/8, maxit = $maxit)",
               "(verbose = $verb, tol = $tol, gamma = 1/8, maxit = $maxit)",
               "(verbose = $verb, tol = $tol, gamma = 1/8, maxit = $maxit)"]

	its = Dict([(sol,0.) for sol in solvers])
	for i in eachindex(solvers)

		setup = set_up()
		solver = eval(parse(solvers[i]*slv_opt[i]))

		suite[solvers[i]] =
		@benchmarkable(it = solve_problem!(solver, setup...),
			       setup = (
					it = 0;
					setup = deepcopy($setup);
					solver = deepcopy($solver) ),
			       teardown = (
					  $its[$solvers[$i]] = it;
					  ),
			       evals = 1, samples = samples, seconds = seconds)
	end

	results = run(suite, verbose = (verb != 0))
	println("TotalVariation its")
	println(its)
	return results
end

function show_results(V, Y, Xt, X, lambda)
	Xd = Xt-V'*(~Y)
	imshow(X)
	imshow(Xt)
	imshow(Xd)
end

end


text:
const year = 365.25 * 24. * 60. * 60.

GtCO2_to_ppm(GtCO2) = GtCO2 / (2.13 * (44. /12.))
tCO2_to_ppm(tCO2) = GtCO2_to_ppm(tCO2) * 1.e-9

ppm_to_GtCO2(ppm) = ppm * (2.13 * (44. /12.))
ppm_to_tCO2(ppm) = ppm_to_GtCO2(ppm) * 1.e9


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

using Markdown
using InteractiveUtils

# ╔═╡ e6e24892-e2b7-472f-9a5c-7d1cbe34eeb6
html"<button onclick=present()> Presente </button>"

# ╔═╡ f92f9d30-b368-11eb-2e66-a7bbf5a4a943
md"""
# Plano da disciplina de Linguagens de Programação

## Ementa da Disciplina:

1. Conceitos primários:

	- Algoritmo, programa, linguagem de programação, compilador.

	- Evolução das linguagens de programação.

2. Desenvolvimento de algoritmos: linguagens para algoritmos e refinamento passo a passo.

3. Comandos de uma linguagem procedimental: atribuição, entrada e saída, condicionais, repetitivos e seletivos.

	- Variáveis escalares e estruturadas homogêneas e heterogêneas. Vetores e encadeamento de estruturas.

4. Subprogramação: funções, procedimentos, passagem de parâmetros, recursividade.  Ponteiros.

5. Noções de programação orientada a objetos.

## Resultados de aprendizagem esperados

##### 🥇 Utilizar computadores e softwares no desenvolvimento de atividades de Engenharia Mecânica.

##### 🥈 Desenvolver programas computacionais em linguagens de programações para solução de problema de Engenharia Mecânica.

##### 🥉 Conhecer ferramentas computacionais e aplicativos capazes de auxílio em atividades relacionados à Engenharia.

## Metodologia

### Estratégia

Com base no livro de ELMÔR FILHO, Gabriel et al. $(html"Uma nova sala de aula é possível: aprendizagem ativa na educação em engenharia.") Rio de Janeiro: LTC, 2019.

Será aplicada a **Sala de aula invertida** com utilização da Estratégia *Co-op co-op* que é uma estratégia cooperativa de aprendizagem ativa.

🧑🎓 Os estudantes trabalham juntos, em pequenas equipes (3 ou 4 membros), para investigar um tópico ou resolver um problema, e produzir um produto de grupo que eles compartilham com toda a classe.

🧑🎓 Os estudantes cooperam dentro de suas equipes para produzir algo que irá beneficiar o grande grupo (Estudos de caso que ficarão disponíveis na internet até o fim da disciplina).

👨🔬 Mediação do professor em todas as etapas de desenvolvimento da estratégia.

## Roteiro

#### Momento Pré

🕐 1. O professor anuncia por meio de uma orientação por escrito (Sala):

   - A atividade a ser desenvolvida.
   - O respectivo tema e subtópicos.
   - Uma síntese da estratégia e
   - O cronograma incluindo o período para as avaliações da atividade.

🕑 2. Os estudantes têm a tarefa de ler sobre o assunto e postar no Classroom resumo de uma página do que estudou (Resumo pré) (Extra sala).

## Roteiro

#### Momento Aula

🕒 3. Rápida verificação do entendimento do tema pelo professor (Sala).

🕓 4. Promoção de discussão sobre os interesses no tema que será estudado. O que querem aprender e experienciar durante a unidade? (Sala).

🕔 5. Seleção do tópico que cada grupo irá estudar. Os estudantes discutem seus interesses e selecionam um tópico para seu grupo (Aula).

🕕 6. Cada grupo divide seu tópico em subtópicos para dividir o trabalho dentro do grupo (Aula).

🕖 7. Preparação dos subtópicos pelos membros do grupo (estudo individual).

🕗 8. Os membros do grupo discutem o tema, como um painel de especialistas, apresentando os conceitos e questões chaves (reunião do grupo).

🕘 9. Apresentação formais dos subtópicos em cada grupo (reunião do grupo).

🕙 10. Apresentações de subtópicos pelos grupos (seleção por sorteio a critério do professsor) (sala).

## Roteiro

#### Momento Pós

🕚 11. Preparação das apresentações de cada um dos grupos (reunião do grupo).

🕛 12. Envio via Classroom do relatório da atividade para postagem na internet.

🕐 13. Apresentações de cada um dos grupos (sala).

🕑 14. Avaliações (Extra sala).

🕒 15. Feedback (sala).

## Recursos e Materiais Didáticos

1. Google Classroom;
2. Google Meet para encontros;
3. Jamboard como quadro branco;
4. Google Forms para formulários e avaliações;
5. Google Drive para arquivos (Docs, Slides, Planilhas, Desenhos);
6. Os recurso $(html"Vscode") e $(html"Linguagem Julia") e $(html"Pluto.jl") para Relatórios.
7. Página do professor para registro dos casos.

## Indicações de Referências

### Bibliografia

1. ASCENCIO, Ana Fernanda Gomes; DE CAMPOS, Edilene Aparecida Veneruchi. Fundamentos da programação de computadores. Pearson Educación, 2008.

2. Documentação do [Julia](https://docs.julialang.org/en/v1/)

3. Artigos e Blogs na Internet.

### Bibliografia (PPC):

4. CELES, W. et al. Introdução a Estruturas de Dados. Rio de Janeiro: Campus-Elsevier,2004.

5. MOKARZEL, F.C.; SOMA, N.Y. Introdução à Ciência da Computação. Rio de Janeiro: Campus-Elsevier, 2008.

6. FEOFILOFF, P. Algoritmos em Linguagem C. Rio de Janeiro: Campus-Elsevier, 2009.

### Bibliografia Complementar  (PPC):

1. DROSDEK, A. Estrutura de Dados e Algoritmos em C++. São Paulo: Thomson,2002.

2. MIZRAHI, V.V. Treinamento em Linguagem C. São Paulo: Pearson, 2008.

3. SALIBA, W. L. C. Técnicas de Programação: uma AbordagemAlgorítmica. São Paulo: Makron, 1992.

4. RUMBAUGH, James et al. Object-oriented modeling and design. Englewood Cliffs: Prentice-hall, 1991.

5. ART, Lew; MAUCH, Holger. Introduction to dynamic programming. Springer Berlin Heidelberg, 2007.

## Procedimento de Avaliação da Apredizagem

Avaliação será:

- Contínua
- Individual
- Por atividade a ser desenvolvida

## Procedimento de Avaliação da Apredizagem

### Regras:

- As avaliações serão feitas via Formulário Eletrônico no Google Forms, com excessão do Resumo pré que deve ser entregue via Classroom.
- Só serão aceitos emails da Unifesspa.
- No caso de duplicidade de avaliação será considerada sempre a última avaliação enviada.
- Serão consideradas apenas as avaliações realizadas no período estipulado, ou até quando estiverem disponíveis no formulário.
- Cada avaliação receberá uma pontuação de 0 a 100.
- A nota para o conceito final será a média das notas das atividades propostas.

## Procedimento de Avaliação da Apredizage

### Composição das Notas

A seguir são enumeradas os itens que serão avaliados.

Para cada atividade serão avaliados:

1. Apresentações dos grupos (peso 1).
2. Contribuições individuais (peso 1).
3. Relatório em arquivo html para postagem (peso 1).

## Procedimento de Avaliação da Apredizagem

### Composição das Notas

#### 1. Apresentações dos grupos (peso 1)

Avalições que serão consideradas neste item:

- Avaliação por pares (peso 1)
- Avaliação do professor (peso 1)

## Procedimento de Avaliação da Apredizagem

### Composição das Notas

#### 2. Contribuições individuais (peso 1)

Avaliações que serão consideradas:

- Avaliação do resumo pré (peso 1)
- Avaliação por pares (peso 1)
- Avaliação do professor (peso 1)

## Procedimento de Avaliação da Apredizagem

### Composição das Notas

#### 3. Relatório em arquivo html para postagem (peso 1).

- Autoavaliação (peso 1)
- Avaliação por pares (peso 1)
- Avaliação do professor (peso 1)

"""

# ╔═╡ Cell order:
# ╟─e6e24892-e2b7-472f-9a5c-7d1cbe34eeb6
# ╠═f92f9d30-b368-11eb-2e66-a7bbf5a4a943

max_stars_repo_path: Plano da Disciplina.jl

max_stars_repo_name: jefferds/Programming-languages

max_issues_repo_path: Plano da Disciplina.jl

max_issues_repo_name: jefferds/Programming-languages

max_forks_repo_path: Plano da Disciplina.jl

max_forks_repo_name: jefferds/Programming-languages