text:
module TimeSeriesA1P2

#Use this to print
#=
weave(Pkg.dir("$(pwd())\\FamaMacbeth.jl"),
  informat="script",
  out_path = "$(pwd())\\FamaMacbethOUT.html",
  doctype = "md2html")
=#

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

using DataFrames, Distributions, StatsBase, GZip, JLD,
  Gadfly, CTIO, CTReg, JuMP, NLopt, HCubature, StaticArrays, CTNumerical,
  ParallelAccelerator, Measures, Formatting

const DATA_PATH = pwd() * "\\data" #path where we store the data
const OUTPUT_PATH = pwd() * "\\output"
const YAHOO_DATE_FORMAT = "m/d/yyyy"
const SP500_NAME = "GSPC"
const DBV_NAME = "DBV"

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} #holds graphs
const OVERRIDE_STAR_STRINGS = ["\\ostar","\\dstar","\\tstar"]

function getNormalData(N::Int; μ::Float64=0., σ²::Float64=1.)::Vector{Float64}
  dist::Distribution = Normal(μ,σ²^0.5)

  return rand(dist, N)
end

function getJumpData(N::Int; μ::Float64=0., σ²::Float64=1.,
  p::Float64=0.5, μⱼ::Float64 = 0.0, σ²ⱼ::Float64 = 1.)::Vector{Float64}

  normDist::Distribution = Normal(μ,σ²^0.5)
  bernDist::Distribution = Bernoulli(p)
  jumpDist::Distribution = Normal(μⱼ,σ²ⱼ^0.5)

  println("kurtosis: ", kurtosis(rand(normDist, N) .+ rand(bernDist, N) .* rand(jumpDist, N)))

  return rand(normDist, N) .+ rand(bernDist, N) .* rand(jumpDist, N)
end

jumpVar(; μ::Float64=0., σ²::Float64=1.,
  p::Float64=0.5, μⱼ::Float64 = 0.0, σ²ⱼ::Float64 = 1.) = σ²+p*(-1*(p-1)*μⱼ^2+σ²ⱼ)

jumpSkew(; μ::Float64=0., σ²::Float64=1.,
  p::Float64=0.5, μⱼ::Float64 = 0.0, σ²ⱼ::Float64 = 1.) =
    ((-1+p)*p*μⱼ*((-1+2*p)*μⱼ^2-3*σ²ⱼ))/(σ²+p*(-1*(-1+p)*μⱼ^2+σ²ⱼ))^(1.5)

jumpKurt(; μ::Float64=0., σ²::Float64=1.,
  p::Float64=0.5, μⱼ::Float64 = 0.0, σ²ⱼ::Float64 = 1.) =
    ((-1+p)*p*((-1-6*(-1+p)*p)*μⱼ^4+6*(-1+2*p)*μⱼ^2*σ²ⱼ-3*σ²ⱼ^2))/(σ²+p*(-1*(p-1)*μⱼ^2+σ²ⱼ))^2

formNum(x::T where T<:Real)::String = format(NUM_FORMAT, x)

#Calculates the summary statistics
function summaryA1P1(normalReturns::Vector{Float64}, jumpReturns::Vector{Float64})
  α::Float64 = 0.05

  #label the table
  colNames::Vector{Vector{String}} =
    [["mean", "\$\\sigma\$", "skewness", "excess kurtosis"]]
  numCols = length(colNames[end])

  rowNames::Vector{String} = ["Normal", "Normal Jumps"]
  numRows::Int = length(rowNames)

  NNormal::Int = length(normalReturns)
  NJumps::Int = length(jumpReturns)

  data::Vector{Vector{Float64}} = [normalReturns, jumpReturns]

  content::Vector{Vector{String}} = [Vector{String}(numCols) for i::Int ∈ 1:numRows]

  #run the tests and record the results
  for r::Int ∈ 1:numRows
    content[r][1] = formNum(mean(data[r]))
    content[r][2] = formNum(std(data[r]))
    content[r][3] = formNum(skewness(data[r]))
    content[r][4] = formNum(kurtosis(data[r]))
  end

  #write the table
  summaryTable::String = texTable( "Normal and Normal Jump Distributions",
    """See Tex File""", #caption
    colNames, #colNames
    Vector{String}(),#contentRowNames
    Vector{Matrix{String}}(), #content
    rowNames, #descRowNames
    content, #descContent
    Vector{String}(), #notes
  )

  return summaryTable
end

#plots the descriptive graphs
function visualizeA1P1(normalReturns::Vector{Float64}, jumpReturns::Vector{Float64})

  #convenience containers
  plots::Vector{PlotContainer} = Vector{PlotContainer}()
  plotNames::Vector{String} = Vector{String}()

  Gadfly.push_theme(:default)
  DF::DataFrame = DataFrame(N = 1:length(normalReturns),
    normalReturns = normalReturns, jumpReturns = jumpReturns)
  outDF::DataFrame = melt(DF, :N) #put data in long form

  push!(plotNames, "Histogram P1")
  push!(plots, plot(outDF, x=:value, color=:variable,
      Guide.title("Histograms"),
      Guide.xlabel("Log Return"),Guide.ylabel("Density"),
      Coord.Cartesian(xmin=-0.2, xmax=0.2),
      Geom.histogram(bincount=20, density=true, position=:dodge)))

  push!(plotNames, "Time Series P1")
  push!(plots, vstack(
        plot(DF, y=:normalReturns, x=:N, Geom.line,
          Guide.xlabel("Year"),Guide.ylabel("Return"), Guide.title("Normal")),
        plot(DF, y=:jumpReturns, x=:N, Geom.line,
          Guide.xlabel("Year"),Guide.ylabel("Return"), Guide.title("Normal with Normal Jumps"))))

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

end

function A1P1Script(; refreshData::Bool = false)
  N::Int = 600
  normalReturns = getNormalData(N, μ=0.008, σ²=0.063^2)
  jumpReturns = getJumpData(N, μ=0.012, σ²=0.05^2, p=0.15, μⱼ=-0.03, σ²ⱼ=0.1^2)
  visualizeA1P1(normalReturns, jumpReturns)
  summaryTable::String = summaryA1P1(normalReturns, jumpReturns)

  println("Jump Var: $(jumpVar(μ=0.012, σ²=0.05^2, p=0.15, μⱼ=-0.03, σ²ⱼ=0.1^2))")
  println("Jump Skew: $(jumpSkew(μ=0.012, σ²=0.05^2, p=0.15, μⱼ=-0.03, σ²ⱼ=0.1^2))")
  println("Jump Kurt: $(jumpKurt(μ=0.012, σ²=0.05^2, p=0.15, μⱼ=-0.03, σ²ⱼ=0.1^2))")

  writeTables2File([summaryTable],
    HEADER_NAME, FOOTER_NAME, path=OUTPUT_PATH,
    outName = "A1P1 Summary.tex")

end

@time begin
#  A1P1Script()
end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
import Base:
  # types
  Array,
  CartesianIndices,
  Vector,
  NTuple,
  Tuple,
  # symbols
  +,
  -,
  *,
  ^,
  /,
  ==,
  <,
  >,
  !,
  # functions
  adjoint,
  allunique,
  axes,
  complex,
  convert,
  copy,
  copyto!,
  deepcopy,
  deleteat!,
  eachindex,
  eltype,
  fill!,
  filter,
  filter!,
  findall,
  findfirst,
  getindex,
  hash,
  intersect,
  intersect!,
  isapprox,
  isassigned,
  isempty,
  isless,
  iterate,
  keys,
  lastindex,
  length,
  map,
  map!,
  ndims,
  permutedims,
  promote_rule,
  push!,
  resize!,
  setdiff,
  setdiff!,
  setindex!,
  show,
  similar,
  size,
  summary,
  truncate,
  zero,
  # macros
  @propagate_inbounds

import Base.Broadcast:
  # types
  AbstractArrayStyle,
  Broadcasted,
  BroadcastStyle,
  DefaultArrayStyle,
  Style,
  # functions
  _broadcast_getindex,
  broadcasted,
  broadcastable,
  instantiate

import ITensors.ContractionSequenceOptimization:
  contraction_cost, optimal_contraction_sequence

import HDF5: read, write

import LinearAlgebra:
  axpby!,
  axpy!,
  dot,
  eigen,
  exp,
  factorize,
  lmul!,
  mul!,
  norm,
  normalize!,
  qr,
  rmul!,
  svd,
  tr

using ITensors.NDTensors:
  EmptyNumber,
  _Tuple,
  _NTuple,
  blas_get_num_threads,
  eachblock,
  eachdiagblock,
  fill!!,
  randn!!,
  timer

import ITensors.NDTensors:
  # Modules
  Strided, # to control threading
  # Types
  AliasStyle,
  AllowAlias,
  NeverAlias,
  # Methods
  array,
  blockdim,
  blockoffsets,
  contract,
  dense,
  denseblocks,
  dim,
  dims,
  disable_tblis,
  eachnzblock,
  enable_tblis,
  ind,
  inds,
  insertblock!,
  insert_diag_blocks!,
  matrix,
  maxdim,
  mindim,
  nblocks,
  nnz,
  nnzblocks,
  nzblock,
  nzblocks,
  outer,
  permuteblocks,
  polar,
  scale!,
  setblockdim!,
  setinds,
  setstorage,
  sim,
  storage,
  sum,
  tensor,
  truncate!,
  using_tblis,
  vector,
  # Deprecated
  addblock!,
  store

import Random: randn!

using SerializedElementArrays: SerializedElementVector

const DiskVector{T} = SerializedElementVector{T}

import SerializedElementArrays: disk

ext: jl


text:
include("./util.jl")

inputs = util.input_lines(3)

digits_at_pos(inputs::Array{String}, pos::Int)::Array{Int} =
    map(l -> parse(Int, l[pos]), inputs)

most_common_bin(inputs::Array{Int}) =
    if sum(inputs) >= length(inputs) / 2
        1
    else
        0
    end

least_common_bin(inputs::Array{Int}) =
    if most_common_bin(inputs) == 1
        0
    else
        1
    end

gamma = 0
epsilon = 0

for i = 1:length(inputs[1])
    g_value = most_common_bin(digits_at_pos(inputs, i))
    e_value = least_common_bin(digits_at_pos(inputs, i))
    global gamma += (g_value << (length(inputs[1]) - i))
    global epsilon += (e_value << (length(inputs[1]) - i))
end

println("Part 1: $(gamma * epsilon)")

function find_rating(inputs::Array{String}, f::Function, curr_pos::Int)::Int
    if length(inputs) == 1
        result = 0
        for (idx, i) in enumerate(inputs[1])
            result += parse(Int, i) << (length(inputs[1]) - idx)
        end
        result
    elseif length(inputs) == 0
        throw("No inputs remaining.")
    else
        target_digit = string(f(digits_at_pos(inputs, curr_pos)))[1]
        find_rating(filter(s -> s[curr_pos] == target_digit, inputs), f, curr_pos + 1)
    end
end

ox_rating = find_rating(inputs, most_common_bin, 1)
co2_rating = find_rating(inputs, least_common_bin, 1)

println("Part 2: $(ox_rating * co2_rating)")


text:
using Statistics, StatsBase
using PyCall
xr = pyimport("xarray")
using Blosc
using Elefridge
using JLD
using ZfpCompression

path = "/Users/milan/cams/gridded"
filelist = filter(x->endswith(x,".grib"),readdir(path))

n = length(filelist)
varnames = fill("",n)
methods = ["Lin24","Log16","R16+Blosc","R14+Blosc","R12+Blosc",
            "R16+LZ4HC","R14+LZ4HC","R12+LZ4HC",
            "ZFP15","ZFP12","ZFP9"]
stats = ["cf","meanerror","absmedian","abs90%","absmax","decmedian","dec90%","decmax"]

E = fill(0f0,n,length(methods),length(stats))

function statscalc(X0::Array{T,N},Xa::Array{T,N}) where {T,N}
    Xamean = mean(Xa)
    X0mean = mean(X0)
    meanerror = abs(X0mean-Xamean)/X0mean       # normalised

    abserr = abs.(X0-Xa)/X0mean                 # normalised
    abserr = sort(vec(abserr))
    abs50 = T(quantile(abserr,.5,sorted=true))
    abs90 = T(quantile(abserr,.9,sorted=true))
    absmax = abserr[end]

    # exception cases for decimal error
    Xa0 = Xa .== zero(T)
    X00 = X0 .== zero(T)

    Xas = signbit.(Xa)
    X0s = signbit.(X0)

    decerr = abs.(log10.(abs.(X0./Xa)))

    decerr[Xas .⊻ X0s] .= Inf        # signs don't match
    decerr[Xa0 .⊻ X00] .= Inf        # one is zero the other one not
    decerr[Xa0 .& X00] .= zero(T)    # both zero, set decimal error to zero

    decerr = sort(vec(decerr))
    dec50 = T(quantile(decerr,.5,sorted=true))
    dec90 = T(quantile(decerr,.9,sorted=true))
    decmax = decerr[end]

    return [meanerror,abs50,abs90,absmax,dec50,dec90,decmax]
end

for (i,file) in enumerate(filelist)
    varname = split(split(file,"_")[end],".")[1]
    varnames[i] = varname
    println("---")
    println("Reading $varname")
    X = xr.open_dataarray(joinpath(path,file),engine="cfgrib").data

    # transpose array
    X = permutedims(X,[3,2,1])  # unravels as longitude, latitude, level

    negatives = any(X .< 0f0)  # check for negative entries

    ## LinQuant24
    E[i,1,1] = 32/24     # compression factor
    Xc = Array{Float32}(LinQuant24Array(X))
    E[i,1,2:end] = statscalc(X,Xc)
    println("LinQuant24 finished.")

    ## LogQuant16
    E[i,2,1] = 32/16
    if negatives       # don't compress
        E[i,2,2:end] .= 0
    else
        Xc = Array{Float32}(LogQuant16Array(X))
        E[i,2,2:end] = statscalc(X,Xc)
    end
    println("LogQuant16 finished.")

    ## Round+lossless
    for (j,r) in enumerate([7,5,3])
        Xc = round(X,r)         # RN16,14,12
        E[i,2+j,2:end] = statscalc(X,Xc)
        E[i,5+j,2:end] = E[i,2+j,2:end]

        Blosc.set_compressor("blosclz")
        E[i,2+j,1] = sizeof(X)/sizeof(compress(Xc,level=5))
        Blosc.set_compressor("lz4hc")
        E[i,5+j,1] = sizeof(X)/sizeof(compress(Xc,level=5))
        println("RN$(r+9)+Blosc/LZ4HC finished.")
    end

    ## Zfp
    for (j,r) in enumerate([15,12,9])
        Xc = zfp_compress(X,precision=r)
        E[i,8+j,1] = sizeof(X)/sizeof(Xc)
        X2 = zfp_decompress(Xc)
        E[i,8+j,2:end] = statscalc(X,X2)
        println("ZFP$r finished.")
    end

    @save "/Users/milan/cams/error/linlogroundzfp_all.jld" varnames methods stats E
end


text:
module CommonFunctions

    export nbasis, nnodes, nodes, order, degree

    nbasis() = nothing
    nnodes() = nothing
    nodes()  = nothing
    order()  = nothing
    degree() = nothing

    export evaluate, evaluate!

    evaluate() = nothing
    evaluate!() = nothing

    export derivative, integral

    derivative() = nothing
    integral() = nothing

    export write_to_hdf5

    write_to_hdf5() = nothing

    export reset!, cut_periodic_solution!

    reset!() = nothing
    cut_periodic_solution!() = nothing

    export eachsample, eachtimestep

    eachsample() = nothing
    eachtimestep() = nothing

end


text:
# the fully specified version
# TODO: use a layered transform parameter to determine if we're returning a st or thinst instead
@doc """
    layers = stParallel(m, Xlength; CWTType=morl, subsampling = 2,
                    outputSize=[1,1], varargs...)
 """
function stParallel(m, Xlength; CWTType = morl, subsampling = 2,
    outputSize = [1, 1], varargs...)
    CWTType = makeTuple(m + 1, CWTType)
    subsampling = makeTuple(m + 1, subsampling)
    pairedArgs = processArgs(m + 1, varargs)
    shears = map(x -> ContinuousWavelets.wavelet(x[1]; x[2]...), zip(CWTType, pairedArgs))
    # @info "Treating as a 1D Signal. Vector Lengths: $Xlength nScales:" *
    #         "$nScales subsampling: $subsampling"
    stParallel{typeof(shears[1]),1,m,typeof(subsampling),typeof(outputSize)}((Xlength,), shears, subsampling, outputSize)
end


text:
# FASTQ Writer
# ============

struct Writer <: BioGenerics.IO.AbstractWriter
    output::IO
    quality_header::Bool
end

function BioGenerics.IO.stream(writer::Writer)
    return writer.output
end

"""
    FASTQ.Writer(output::IO; quality_header=false)

Create a data writer of the FASTQ file format.

# Arguments
* `output`: data sink
* `quality_header=false`: output the title line at the third line just after '+'
"""
function Writer(output::IO; quality_header::Bool=false)
    return Writer(output, quality_header)
end

function Base.write(writer::Writer, record::Record)
    checkfilled(record)
    output = writer.output
    n = 0
    # sequence
    n += write(output, '@')
    n += unsafe_write(output, pointer(record.data, first(record.identifier)), length(record.identifier))
    if hasdescription(record)
        n += write(output, ' ')
        n += unsafe_write(output, pointer(record.data, first(record.description)), length(record.description))
    end
    n += write(output, '\n')
    n += unsafe_write(output, pointer(record.data, first(record.sequence)), length(record.sequence))
    n += write(output, '\n')
    # quality
    n += write(output, '+')
    if writer.quality_header
        n += unsafe_write(output, pointer(record.data, first(record.identifier)), length(record.identifier))
        if hasdescription(record)
            n += write(output, ' ')
            n += unsafe_write(output, pointer(record.data, first(record.description)), length(record.description))
        end
    end
    n += write(output, '\n')
    n += unsafe_write(output, pointer(record.data, first(record.quality)), length(record.quality))
    n += write(output, '\n')
    return n
end


text:
@testset "compare" begin

    @test ArbNumerics.Mag(2.0) == ArbNumerics.Mag(2.0)
    @test ArbNumerics.Mag(1.0) != ArbNumerics.Mag(2.0)
    @test ArbNumerics.Mag(2.0) > ArbNumerics.Mag(1.0)
    @test ArbNumerics.Mag(1.0) < ArbNumerics.Mag(2.0)
    @test ArbNumerics.Mag(2.0) >= ArbNumerics.Mag(1.0)
    @test ArbNumerics.Mag(1.0) <= ArbNumerics.Mag(2.0)

    @test ArbFloat(2.0) > 1
    @test ArbFloat(2.0) >= 1
    @test 1 < ArbFloat(2.0)
    @test 1 <= ArbFloat(2.0)

    @test ArbFloat(2.0) > 1.0
    @test 1.0 < ArbFloat(2.0)
    @test ArbFloat(2.0) >= 1.0
    @test 1.0 <= ArbFloat(2.0)

    @test ArbFloat(2.0) == ArbFloat(2.0)
    @test ArbFloat(1.0) != ArbFloat(2.0)
    @test ArbFloat(2.0) > ArbFloat(1.0)
    @test ArbFloat(1.0) < ArbFloat(2.0)
    @test ArbFloat(2.0) >= ArbFloat(1.0)
    @test ArbFloat(1.0) <= ArbFloat(2.0)

    @test ArbFloat(2.0) == ArbReal(2.0)
    @test ArbFloat(1.0) != ArbReal(2.0)
    @test ArbFloat(2.0) > ArbReal(1.0)
    @test ArbFloat(1.0) < ArbReal(2.0)
    @test ArbFloat(2.0) >= ArbReal(1.0)
    @test ArbFloat(1.0) <= ArbReal(2.0)

    @test ArbReal(2.0) == ArbFloat(2.0)
    @test ArbReal(1.0) != ArbFloat(2.0)
    @test ArbReal(2.0) > ArbFloat(1.0)
    @test ArbReal(1.0) < ArbFloat(2.0)
    @test ArbReal(2.0) >= ArbFloat(1.0)
    @test ArbReal(1.0) <= ArbFloat(2.0)

    @test ArbReal(2.0) > 1
    @test ArbReal(2.0) >= 1
    @test 1 < ArbReal(2.0)
    @test 1 <= ArbReal(2.0)

    @test ArbReal(2.0) > 1.0
    @test 1.0 < ArbReal(2.0)
    @test ArbReal(2.0) >= 1.0
    @test 1.0 <= ArbReal(2.0)

    @test ArbReal(2.0) == ArbReal(2.0)
    @test ArbReal(1.0) != ArbReal(2.0)
    @test ArbReal(2.0) > ArbReal(1.0)
    @test ArbReal(1.0) < ArbReal(2.0)
    @test ArbReal(2.0) >= ArbReal(1.0)
    @test ArbReal(1.0) <= ArbReal(2.0)

end


text:
using
  PyPlot,
  Printf,
  FourierFlows

using FFTW: rfft

import GeophysicalFlows.TwoDTurb

 dev = CPU()   # Device (CPU/GPU)

  nx = 256     # Resolution
  Lx = 2π      # Domain size
   ν = 1e-6    # Viscosity
  nν = 1       # Order of (hyper-)viscosity. nν=1 means Laplacian
  dt = 0.1     # Timestep
nint = 200     # Number of steps between plots
ntot = 10nint  # Number of total timesteps

# Define problem
prob = TwoDTurb.Problem(nx=nx, Lx=Lx, ν=ν, nν=nν, dt=dt, stepper="FilteredRK4", dev=dev)
TwoDTurb.set_zeta!(prob, ArrayType(dev)(rand(nx, nx)))

cl, vs, gr = prob.clock, prob.vars, prob.grid
x, y = gridpoints(gr)

"Plot the vorticity of the twodturb problem `prob`."
function makeplot!(ax, prob)
  sca(ax)
  cla()
  pcolormesh(x, y, vs.zeta)
  title("Vorticity")
  xlabel(L"x")
  ylabel(L"y")
  pause(0.01)
  nothing
end

# Step forward
fig1, ax = subplots(figsize=(8, 8))

while cl.step < ntot
  @time begin
    stepforward!(prob, nint)
    @printf("step: %04d, t: %6.1f", cl.step, cl.t)
  end

  TwoDTurb.updatevars!(prob)
  makeplot!(ax, prob)
end

# Plot the radial energy spectrum
E  = @. 0.5*(vs.u^2 + vs.v^2) # energy density
Eh = rfft(E)
kr, Ehr = FourierFlows.radialspectrum(Eh, gr, refinement=1)

fig2, axs = subplots(ncols=2, figsize=(8, 4))

sca(axs[1])
pcolormesh(x, y, vs.zeta)
xlabel(L"x")
ylabel(L"y")
title("Vorticity")

sca(axs[2])
plot(kr, abs.(Ehr))
xlabel(L"k_r")
ylabel(L"\int | \hat{E} | k_r \mathrm{d} k_{\theta}")
title("Radial energy spectrum")

xlim(0, nx/4)
axs[2].set_yscale("log")

tight_layout(w_pad=0.1)
show()


text:
using Mmap
using JLD

#%%
x = rand(100, 100)
data_dict = Dict("x" => x)

save("files/test-file.jld", data_dict)

function split_file(;filePath::String, fileSizeMb::Real, folderPath="default")
    if folderPath=="default"
        folderPath=dirname(filePath)
    end

    fileName=basename(filePath)

    Kbyte = 1024
    Mbyte = 1024 * 1024
    maxFileSize = fileSizeMb * Mbyte
    nfiles = ceil(Int, bytes / maxFileSize)

    splits_indexes = collect(ceil.(Int, range(0, bytes, length = (files + 1))))

    splits = [file[(splits_indexes[i]+1):(splits_indexes[i+1])] for i = 1:(length(splits_indexes)-1)]
    let s = Mmap.open("$folderPath/$fileName.0", "w+")
        write(s, Int(length(splits_indexes)-1))
        close(s)
    end

    for (i, split) in enumerate(splits)
        let s = Mmap.open("$folderPath/$fileName.$i", "w+")
            write(s, split)
            close(s)
        end
    end
    close(s)
end

function recompose_file(;filePath::String, folderPath="default")
    if folderPath=="default"
        folderPath=dirname(filePath)
    end

    fileName=basename(filePath)

    s0 = Mmap.open("$folderPath/$fileName.0", "r")

    nfiles = read(s0, Int)
    close(s0)

    restored = Vector{UInt8}[]

    for i = 1:nfiles
        let s = Mmap.open("$folderPath/$fileName.$i", "r")
            file = read(s)
            push!(restored, file)
            close(s)
        end
    end

    let s = Mmap.open("$folderPath/$fileName", "w+")
        file = vcat(restored...)
        write(s, file)
        close(s)
    end
    @info "$fileName recomposed"
    return
end
# %%
split_file(filePath="files/test-file.jld", folderPath="files2", fileSizeMb=0.01)
recompose_file(filePath="files2/test-file.jld")
# %%

data_dict2 = load("files/test-file_restored.jld")
# data_dict2 = load("test-file.jld")
data_dict2["x"]==x

basename("files/test-file_restored.jld")

#%%
struct Person
    height::Float64
    weight::Float64
    BMI::Float64

    function Person(heightInMeters::Float64, weightInKilos::Float64)
        BMI = weightInKilos / heightInMeters^2
        new(heightInMeters, weightInKilos, BMI)
    end

    function Person(;heightInMeters::Float64, weightInKilos::Float64)
        BMI = weightInKilos / heightInMeters^2
        new(heightInMeters, weightInKilos, BMI)
    end
end

bob = Person(heightInMeters=1.84, weightInKilos=84.0)
carl = Person(1.84, 84.0)

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
@testset "call_tetgen.jl" begin

end


text:
using LinearAlgebra, Plots, Statistics, TaylorSeries, Printf, FileIO
# if !(pwd() in LOAD_PATH) push!(LOAD_PATH, pwd()) end
# using TaylorBigF

path1=realpath(dirname(@__FILE__)*"/..")
include(string(path1,"/src/TaylorBigF.jl"))

# 3.7 System Identification

p=2000;setprecision(p); setrounding(BigFloat, RoundUp);
L=BigFloat("1"); n=BigInt(100); iN=convert(Int64,n);
dx=BigFloat("1")*L/(n-1); x=convert(Array{BigFloat},0:dx:(n-1)*dx)
I=zeros(BigFloat, n, n);for i=1:n I[i,i]=BigFloat("1") end

V=TaylorBigF.calc_Vandermonde(x,n,n)
dV=TaylorBigF.calc_Vandermonde_dot(x,n,n)
ddV=TaylorBigF.calc_Vandermonde_dot_dot(x,n,n)
dddV=TaylorBigF.calc_Vandermonde_dot_dot_dot(x,n,n)
ddddV=TaylorBigF.calc_Vandermonde_dot_dot_dot_dot(x,n,n)

# x=time (t in manuscript) s=space dds=a=F/m
s=x.^BigInt(2)
# plot(x,s)

a=V\s
# plot(a)

# s_dot=s[2:end]-s[1:end-1]
# s_dot_dot=s_dot[2:end]-s_dot[1:end-1]
# VV=[ones(n-2) s[2:end-1] s_dot[1:end-1] s_dot_dot]
# aa=VV\ones(n-2)
# plot(aa)

# ss=V*a;ssd=dV*a;ssdd=ddV*a;
# VV=[ss ssd ssdd ss.^2 ss.*ssd ss.*ssdd ssd.^2 ssd.*ssdd ssdd.^2]
VV=[V*a dV*a ddV*a]
b=VV\ones(n)
# plot(b)

# xi=x[end]+dx:dx:x[end]+100dx
xi=x[end]+BigFloat("9525000000000")dx:dx:x[end]+BigFloat("9525000000010")*dx
Vi=TaylorBigF.calc_Vandermonde(xi,length(xi),n)
IVi=TaylorBigF.calc_Vandermonde_int(xi,length(xi),n)
IIVi=TaylorBigF.calc_Vandermonde_int_int(xi,length(xi),n)
IV=TaylorBigF.calc_Vandermonde_int(x,length(x),n)
IIV=TaylorBigF.calc_Vandermonde_int_int(x,length(x),n)

# t2=(x.^BigInt(2))/BigInt(2)
# mat1=(-t2+b[1]*IIV*a+b[2]*IV*a+b[3]*V*a)
# cc=[s[1] BigInt(1);s[end] BigInt(1)]\[mat1[1,:];mat1[end,:]]
# s_comp=(t2+cc[1].*x+cc[2].*ones(BigFloat,length(x))-b[1]*IIV*a-b[2]*IV*a)/b[3]
# # plot(s_comp)#confirmation
# t2=(xi.^BigInt(2))/BigInt(2)
# s_comp=(t2+cc[1].*xi+cc[2].*ones(BigFloat,length(xi))-b[1]*IIVi*a-b[2]*IVi*a)/b[3] # Equation 5
# plot(x,s,legend=:topleft,xaxis="t",yaxis="s",label="given s")
# plot!(xi,s_comp,label=string("identified s, p=",p))
# # savefig("syst_id.pdf")

t2=(xi.^BigInt(2))/BigInt(2)
s_comp=(t2+-b[1]*IIVi*a-b[2]*IVi*a)/b[3] # Equation 5
# plot(x,s,legend=:topleft,xaxis="t",yaxis="s",label="given s")
# plot!(xi,s_comp,label=string("identified s, p=",p))
# savefig("syst_id.pdf")

s_analyt=xi.^BigInt(2)
error1=s_comp-s_analyt
iok=abs.(error1).<1
xi[iok][end]


text:
"""
Print out a dashed line of certain length.

Arguments
---------
length:Integer
    length of the dashed line

marker="-":AbstractString
    (optional) maker of the dashed line
"""
function print_dashed_line(length::Integer, marker::AbstractString="-")
    if length > 0
        print("-")
        print_dashed_line(length-1)
    elseif length == 0
        print("\n")
    else
        print("")
    end
end


text:
## using_glmakie
using GLMakie

## glmakie_lineplot
f, ax, l = lines(1:10, 1:10)
mktempdir() do path
    save(joinpath(path, "test.png"), f)
end


text:
function default_fmt(T::Type{S}) where {S}
    error("Format for $T is not defined")
end
"""
Usage is: @fmt T = fmtstring

Available types and their default formats are:
    1. Real            "{:.3g}"
    2. Int             "{:,n}"
    3. Bool            "{:}"
    4. AbstractString  "{:}"
"""
macro fmt(ex)
    msg = """
    Usage is: @fmt T = fmtstring

    Available types and their default formats are:
        1. Real            "{:.3g}"
        2. Int             "{:,n}"
        3. Bool            "{:}"
        4. AbstractString  "{:}"
    """
    @assert(ex.head == :(=), msg)
    @assert(length(ex.args) == 2, msg)
    @assert(ex.args[1] isa Symbol, msg)
    @assert(isa(ex.args[2], String), msg)
    @assert(ex.args[1] in [:Real, :Int, :Bool, :AbstractString, :Missing], msg)

    ex1 = ex.args[1]
    ex2 = ex.args[2]

    q = quote
        TexTables.default_fmt(T::Type{S}) where {S <: $ex1} = $ex2
    end
    return q
end

@fmt Real = "{:.3g}"
@fmt Int  = "{:,n}"
@fmt Bool = "{:}"
@fmt AbstractString = "{:}"
@fmt Missing = ""

const _fmt_spec_gG = r"[gG]"
const _fmt_spec_g  = r"[g]"
const _fmt_spec_G  = r"[G]"

function fixed_or_scientific(val, format)
    if occursin(_fmt_spec_gG, format)
        if val isa Integer
            return replace(format, _fmt_spec_gG => "n")
        else
            mag = log(abs(val))/log(10)
            if  (-Inf < mag <= -3) | (mag >= 5)
                r = "e"
            else
                r = "f"
            end

            if occursin(_fmt_spec_g, format)
                return replace(format, _fmt_spec_g => lowercase(r))
            else
                return replace(format, _fmt_spec_G => uppercase(r))
            end
        end
    end
    return format
end

abstract type FormattedNumber{T} end

mutable struct FNum{T} <: FormattedNumber{T}
    val::T
    star::Int
    format::String
    function FNum(val::T, star::Int, format::String) where T
        return new{T}(val, star, fixed_or_scientific(val, format))
    end
end

==(x1::FNum, x2::FNum) =    x1.val      == x2.val &&
                            x1.format   == x2.format &&
                            x1.star     == x2.star

mutable struct FNumSE{T} <: FormattedNumber{T}
    val::T
    se::Float64
    star::Int
    format::String
    format_se::String

    function FNumSE(val::T, se::Float64, star::Int, format::String,
                    format_se::String) where T
        return new{T}(val, se, star, fixed_or_scientific(val, format),
                      fixed_or_scientific(se, format_se))
    end
end

==(x1::FNumSE, x2::FNumSE) =   x1.val       == x2.val       &&
                               x1.se        == x2.se        &&
                               x1.star      == x2.star      &&
                               x1.format    == x2.format    &&
                               x1.format_se == x2.format_se

function FormattedNumber(val::T, format::String=default_fmt(T)) where T
    return FNum(val, 0, format)
end

function FormattedNumber(val::T, se::S,
             format::String=default_fmt(T),
             format_se::String=default_fmt(S)) where
             {T<:AbstractFloat, S <: AbstractFloat}
    se2 = Float64(se)
    newval, newse = promote(val, se)
    return FNumSE(newval, newse, 0, format, format_se)
end

function FormattedNumber(val::T, se::S,
                         format::String=default_fmt(T),
                         format_se::String=default_fmt(S)) where
                         {T, S<:AbstractFloat}
    se2 = Float64(se)
    @assert(isnan(se), "Cannot have non-NaN Standard Errors for $T")
    return FNumSE(val, se, 0, format, format_se)
end

FormattedNumber(x::FormattedNumber) = x

# Unpack Tuples of Floats for precision
FormattedNumber(x::Tuple{T1, T2}) where {T1<: AbstractFloat,
                                         T2<:AbstractFloat} = begin
    return FormattedNumber(x[1], x[2])
end

Base.show(io::IO, x::FNum) = print(io, format(x.format, x.val))
Base.show(io::IO, x::FNumSE)= begin
    if isnan(x.se)
        print(io, value(x))
    else
        str = string(value(x), " ", se(x))
        print(io, str)
    end
end

Base.convert(::Type{FNumSE}, x::FNum) = FormattedNumber(x.val, NaN)
Base.promote_rule(::Type{FNumSE{T}}, ::Type{FNum{S}}) where {T,S} = FNumSE

value(x::FormattedNumber)   = format(x.format, x.val)
se(x::FNumSE)               = format("($(x.format_se))", x.se)
se(x::FNum)                 = ""

function star!(x::FormattedNumber, i::Int)
    x.star = i
    return x
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -31.730360167469108 <= q <= 25.84895899492546)

                     add_NL_constraint(m, :(log(1 + exp(0.0048766473277672695 + 0.6370542772950833*log(1 + exp(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5]))) + 0.7379502189935647*log(1 + exp(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5]))) + -0.017781498731094914*log(1 + exp(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5]))) + -0.30964053180650675*log(1 + exp(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5]))) + -0.19693223852539354*log(1 + exp(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))))) + log(1 + exp(0.6750746912592471 + -0.6689814871888631*log(1 + exp(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5]))) + 0.08174772592079638*log(1 + exp(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5]))) + 0.43040165086652493*log(1 + exp(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5]))) + 0.6360324545158482*log(1 + exp(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5]))) + 0.2590621611090005*log(1 + exp(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))))) + log(1 + exp(-0.5603004507305789 + -0.6080760485894823*log(1 + exp(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5]))) + -0.17827205367530619*log(1 + exp(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5]))) + -0.3024929990799108*log(1 + exp(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5]))) + 0.2876553916986979*log(1 + exp(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5]))) + 0.1394466763243627*log(1 + exp(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))))) + log(1 + exp(-0.3104079613918116 + -0.8967757209559282*log(1 + exp(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5]))) + 0.9554659562201504*log(1 + exp(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5]))) + 0.4284452595185697*log(1 + exp(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5]))) + -0.8546417933790886*log(1 + exp(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5]))) + -0.9550224035691546*log(1 + exp(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))))) + log(1 + exp(0.11524778075403264 + -0.08842997013200282*log(1 + exp(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5]))) + -0.2446780712327863*log(1 + exp(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5]))) + 0.5457935761049009*log(1 + exp(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5]))) + -0.1673283141057702*log(1 + exp(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5]))) + 0.6343389483483679*log(1 + exp(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module TestPlotter

using ..Example: schedule_timings
using MemoryAccesses
using Plots
using Test

function shows_something(mime, obj)
    io = IOBuffer()
    origin = position(io)
    show(io, mime, obj)
    return position(io) > origin
end

function test_show()
    table = MemoryAccesses.table(schedule_timings(100))
    @test shows_something("text/plain", table)
end

function test_plot()
    @testset for transform in [:collapse, :pack, :none]
        check_plot(transform)
    end
end

function check_plot(transform::Symbol)
    table = MemoryAccesses.table(schedule_timings(100); transform = transform)
    @test shows_something("image/png", plot(table))
    @test shows_something("image/png", MemoryAccesses.location_time(table))
    @test shows_something("text/plain", MemoryAccesses.location_time(table))
end

end  # module


text:
"""
    LinkedList

Linked list abstract supertype.
"""
abstract type LinkedList{T} end

"""
    LLNode

Linked list node supertype.
"""
abstract type LLNode{T} end

"""
    NValueType

Node value type.
"""
NValueType = Union{T, Nothing} where T

"""
SLLNode

Singly linked list node. Holds pointer to next node and some value.
"""
mutable struct SLLNode{T} <: LLNode{T}
    next::SLLNode{T}
    value::NValueType{T}
    SLLNode{T}() where T = (x = new(); x.next = x; x.value = nothing; x)
    SLLNode{T}(y::T) where T = (x = new(); x.next = x; x.value = y; x)
    SLLNode{T}(n::SLLNode{T}, y::T) where T = new(n, y)
end

"""
    DLLNode

Double linked list node. Holds pointers to previous and next nodes and some value.
"""
mutable struct DLLNode{T} <: LLNode{T}
    prev::DLLNode{T}
    next::DLLNode{T}
    value::NValueType{T}
    DLLNode{T}() where T = (x = new(); x.prev = x; x.next = x; x.value = nothing; x)
    DLLNode{T}(y::T) where T = (x = new(); x.prev = x; x.next = x; x.value = y; x)
    DLLNode{T}(p::DLLNode{T}, n::DLLNode{T}, y::T) where T = new(p, n, y)
end

"""
    SLinkedList

A singly linked list implementation.
"""
mutable struct SLinkedList{T} <: LinkedList{T}
    head::SLLNode{T}
    tail::SLLNode{T}
    size::Int
    #SLinkedList{T}() where T = (x = new(); x.size = 0; x)
    SLinkedList{T}() where T = (x = SLLNode{T}(); new(x, x, 0))
end

"""
    DLinkedList

A double linked list implementation.
"""
mutable struct DLinkedList{T} <: LinkedList{T}
    head::DLLNode{T}
    tail::DLLNode{T}
    size::Int
    #DLinkedList{T}() where T = (x = new(); x.size = 0; x)
    DLinkedList{T}() where T = (x = DLLNode{T}(); new(x, x, 0))
end

"""
    addlast!(l::SLinkedList{T}, x::T)

Add node to last position in linked list. Time complexity O(1).
"""
function addlast!(l::SLinkedList{T}, x::T) where T
    l.tail.next.value = x
    l.tail = l.tail.next
    l.tail.next = SLLNode{T}()
    l.size += 1

    return nothing
end

"""
    addlast!(l::DLinkedList{T}, x::T)

Add node to last position in linked list. Time complexity O(1).
"""
function addlast!(l::DLinkedList{T}, x::T) where T
    l.tail.next.value = x
    l.tail = l.tail.next
    l.tail.next = DLLNode{T}()
    l.tail.next.prev = l.tail
    l.size += 1

    return nothing
end

"""
    addfirst!(l::SLinkedList{T}, x::T)

Add node to first position in linked list. Time complexity O(1).
"""
function addfirst!(l::SLinkedList{T}, x::T) where T
    n = l.head
    l.head = SLLNode{T}(x)
    l.head.next = n
    if l.size == 0 # do this only if called first on empty list
        l.tail = l.head
    end
    l.size += 1

    return nothing
end

"""
    addfirst!(l::DLinkedList{T}, x::T)

Add node to first position in linked list. Time complexity O(1).
"""
function addfirst!(l::DLinkedList{T}, x::T) where T
    n = l.head
    l.head = DLLNode{T}(x)
    l.head.next = n
    l.head.next.prev = l.head
    if l.size == 0 # do this only if called first on empty list
        l.tail = l.head
    end
    l.size += 1

    return nothing
end

"""
    add!(l::SLinkedList{T}, x::T, pos::Int)

Add node to linked list at some position. Time complexity O(n).
"""
function add!(l::SLinkedList{T}, x::T, pos::Int) where T
    if pos <= 1
        addfirst!(l, x)
    elseif pos > l.size
        addlast!(l, x)
    else
        n = l.head
        i = 1
        while i < pos
            n = n.next
            i += 1
        end
        n.next = SLLNode{T}(n.next, n.value)
        n.value = x
        l.size += 1
    end

    return nothing
end

"""
    add!(l::DLinkedList{T}, x::T, pos::Int)

Add node to linked list at some position. Time complexity O(n).
"""
function add!(l::DLinkedList{T}, x::T, pos::Int) where T
    if pos <= 1
        addfirst!(l, x)
    elseif pos > l.size
        addlast!(l, x)
    else
        n = l.head
        i = 1
        while i < pos
            n = n.next
            i += 1
        end
        n.next = DLLNode{T}(n, n.next, n.value)
        n.value = x
        n.next.next.prev = n.next
        l.size += 1
    end

    return nothing
end

"""
    removelast!(l::SLinkedList)

Remove last node from linked list. Time complexity O(n).
"""
function removelast!(l::SLinkedList)
    x = l.tail.value
    n = l.head

    i = 1
    while i < l.size - 1
        n = n.next
        i += 1
    end
    empty = l.tail.next

    if l.size == 1
        # all this necessary?
        l.head = empty
        l.head.next = l.head
        l.tail = empty
        l.tail.next = l.tail
    else
        l.tail = n
        l.tail.next = empty
    end
    l.size -= 1

    return x
end

"""
    removelast!(l::DLinkedList)

Remove last node from linked list. Time complexity O(1).
"""
function removelast!(l::DLinkedList)
    x = l.tail.value
    empty = l.tail.next

    if l.size == 1
        # all this necessary?
        l.head = empty
        l.head.prev = l.head
        l.head.next = l.head
        l.tail = empty
        l.tail.prev = l.tail
        l.tail.next = l.tail
    else
        l.tail = l.tail.prev
        l.tail.next = empty
        l.tail.next.prev = l.tail
    end
    l.size -= 1

    return x
end

"""
    removefirst!(l::SLinkedList)

Remove first node from linked list. Time complexity O(1).
"""
function removefirst!(l::SLinkedList)
    x = l.head.value
    l.head = l.head.next
    if l.size == 1
        l.tail = l.head
    end
    l.size -= 1

    return x
end

"""
    removefirst!(l::DLinkedList)

Remove first node from linked list. Time complexity O(1).
"""
function removefirst!(l::DLinkedList)
    x = l.head.value
    l.head = l.head.next
    l.head.prev = l.head
    if l.size == 1
        l.tail = l.head
    end
    l.size -= 1

    return x
end

"""
    remove!(l::SLinkedList, pos::Int)

Remove node from linked list. Time complexity O(n).
"""
function remove!(l::SLinkedList, pos::Int)
    x = 0

    if pos <= 1
        x = removefirst!(l)
    elseif pos > l.size
        x = removelast!(l)
    else
        n = l.head

        i = 1
        while i < l.size-2
            n = n.next
            i += 1
        end
        x = n.next.value
        n.next = n.next.next

        l.size -= 1
    end

    return x
end

"""
    remove!(l::DLinkedList, pos::Int)

Remove node from linked list. Time complexity O(n).
"""
function remove!(l::DLinkedList, pos::Int)
    x = 0

    if pos <= 1
        x = removefirst!(l)
    elseif pos > l.size
        x = removelast!(l)
    else
        n = l.head

        i = 1
        while i < l.size - 1
            n = n.next
            i += 1
        end
        x = n.value
        n.prev.next = n.next
        n.next.prev = n.prev

        l.size -= 1
    end

    return x
end

"""
    add!(l::LinkedList{T}, x::T)

Add node to last position in linked list. Time complexity O(1).
"""
add!(l::LinkedList{T}, x::T) where T = addlast!(l, x)

"""
    remove!(l::LinkedList)

 Remove first node from linked list. Time complexity O(1).
"""
remove!(l::LinkedList) = removefirst!(l)

"""
    peekfirst(l::LinkedList)

Lookup first element in list.
"""
peekfirst(l::LinkedList) = l.head.value

"""
    peeklast(l::LinkedList)

Lookup last element in list.
"""
peeklast(l::LinkedList) = l.tail.value

"""
    length(l)

Get length of linked list.
"""
length(l::LinkedList) = l.size

"""
    size(l)

Get length of linked list.
"""
size(l::LinkedList) = length(l)

"""
    :(==)

Test list equality. O(n)
"""
function (==)(x::LinkedList, y::LinkedList)
    collect(x) == collect(y)
end
# (==)(x::LinkedList, y::LinkedList) = (x.head == y.head) && (x.tail == y.tail) && (x.size == y.size) # this is not complete

"""
    show(l)

Display linked list types.
"""
function show(io::IO, l::SLinkedList{T}) where T
    println(io, "$(typeof(l))")
    print(io, " Head: ")

    n = l.head
    while n != n.next
        print(io, " $(n.value)  -> ")
        n = n.next
    end
    println(io, " $(n.value)")
    println(io, " Tail:  $(l.tail.value)")
    print(io, " Size:  $(l.size)")
end

function show(io::IO, l::DLinkedList{T}) where T
    println(io, "$(typeof(l))")
    print(io, " Head: ")

    n = l.head
    while n != n.next
        print(io, " $(n.value)  <-> ")
        n = n.next
    end
    println(io, " $(n.value)")
    println(io, " Tail:  $(l.tail.value)")
    print(io, " Size:  $(l.size)")
end

"""
    iterate(l, n)

Iterate over a linked list.
"""
iterate(l::LinkedList, n = l.head) = n == n.next ? nothing : (n.value, n.next)

"""
    eltype(type)

Type of elements generated by iteration over a linked list.
"""
eltype(::LinkedList{T}) where T = T


text:
# Contains code helpful for using transient decay models

struct DecayParameter
    # Generic type for configuring sub-module
    val::AbstractFloat
    fixed::Bool
    glob::Bool
end

function param_vector(pmain::AbstractVector{<:AbstractFloat},
    fxparam::AbstractVector{<:Bool},
    gparam::AbstractVector{<:Bool})

    # Converts guesses to DecayParameters

    # Ensure sizes match
    if length(pmain) != length(fxparam) ||
        length(fxparam) != length(gparam)
        error("Parameter vectors must have same length ($(length(pmain))")
    end
    npar = length(pmain)
    param = Vector{DecayParameter}(undef,npar)
    # Fill in vector
    for i = 1:npar
        param[i] = DecayParameter(pmain[i],fxparam[i],gparam[i])
    end
    return param
end

struct DecayComponent
    # Stores amplitude-weighted exponential component
    α::DecayParameter
    τ::DecayParameter
end

struct TcspcFitModel
    # Configuration for fitting data
    τ::AbstractVector{DecayComponent}
    θ::AbstractVector{DecayComponent}
    Q::DecayParameter
    S::DecayParameter
    V::DecayParameter
    Z::DecayParameter
    function TcspcFitModel(pmain::AbstractVector{DecayParameter},
        ntau::Integer,nani::Integer)
        # Builds fit model using parameter vector

        # Ensure that dimensions match specs
        if length(pmain) != ((2 * ntau) + (2 * nani) + 4)
            error("Parameter vector must be $(((2 * ntau) + (2 * nani) + 4)) elements long")
        end
        if ntau < 1
            error("Must provide at least one lifetime")
        end
        if nani < 0
            error("Number of anisotropies must be nonnegative")
        end
        α = pmain[1:2:(2 * ntau)]
        τ = pmain[2:2:(2 * ntau)]
        τvec = Vector{DecayComponent}(undef,ntau)
        for i = 1:ntau
            τvec[i] = DecayComponent(α[i],τ[i])
        end
        if nani > 0
            sa = (2 * ntau) + 1
            ea = sa + (2 * nani)
            β = pmain[sa:2:ea]
            θ = pmain[(sa + 1):2:ea]
            θvec = Vector{DecayComponent}(undef,nani)
            for i = 1:nani
                θvec[i] = DecayComponent(β[i],θ[i])
            end
        else
            θvec = Vector{DecayComponent}(undef,0)
        end
        Q = pmain[(end - 3)]
        S = pmain[(end - 2)]
        V = pmain[(end - 1)]
        Z = pmain[end]
        new(τvec,θvec,Q,S,V,Z)
    end
end

function TcspcFitModel(pmain::AbstractVector{<:AbstractFloat},
    fxparam::AbstractVector{Bool},
    gparam::AbstractVector{Bool},
    ntau::Integer,nani::Integer,Q::Bool)
    # Wrapper for TcspcFitModel construction
    param = param_vector(pmain,fxparam,gparam)
    if !Q
        # Add a temporary, fixed, non-global color shift with val = 0
        q = DecayParameter(0,true,false)
        np = length(param) + 1
        ptemp = Vector{DecayParameter}(undef,np)
        ind = vcat(1:(np - 4),(np - 2):np)
        ptemp[ind] = param
        ptemp[(np - 3)] = q
        param = ptemp
    end
    TcspcFitModel(param,ntau,nani)
end

nlcomp(mod::TcspcFitModel) = length(mod.τ)

nacomp(mod::TcspcFitModel) = length(mod.θ)

npar(mod::TcspcFitModel) = 2 * (nlcomp(mod) + nacomp(mod)) + 3

function npar(arr::AbstractVector{TcspcFitModel})
    # Vectorized npar
    np = zeros(Int,length(arr))
    modcount = 1
    for mod ∈ arr
        np[modcount] = npar(mod)
        modcount += 1
    end
    return np
end

function param_buffer(mod::TcspcFitModel)
    # Groups parameters for easy access in fitting

    pmain = zeros(npar(mod),3)
    # Column 1: Values
    # Column 2: Fixed/free (specified as 1)
    # Column 3: Global/non-global (specified as 1)

    # Start with tau
    τᵢ = 1
    for τ ∈ mod.τ
        pmain[2 * τᵢ - 1,:] = [τ.α.val,τ.α.fixed,τ.α.glob]
        pmain[2 * τᵢ,:] = [τ.τ.val,τ.τ.fixed,τ.τ.glob]
        τᵢ += 1
    end
    if nacomp(mod) > 0
        θᵢ = τᵢ
        for θ ∈ mod.θ
            pmain[2 * θᵢ - 1,:] = [θ.α.val,θ.α.fixed,θ.α.glob]
            pmain[2 * θᵢ,:] = [θ.τ.val,θ.τ.fixed,θ.τ.glob]
            θᵢ += 1
        end
    end
    pmain[end - 2,:] = [mod.S.val,mod.S.fixed,mod.S.glob]
    pmain[end - 1,:] = [mod.V.val,mod.V.fixed,mod.V.glob]
    pmain[end,:] = [mod.Z.val,mod.Z.fixed,mod.Z.glob]
    return pmain
end

function isunique(arr::AbstractVector{TcspcFitModel})
    # Determines whether the models are the same

    # Let's assume they are at first
    for i = 2:length(arr)
        modₐ = arr[i - 1]
        modᵦ = arr[i]
        # First, check if length of parameters is the same
        if ncomp(modₐ) != ncomp(modᵦ)
            return false
        end
        # Now, check lifetimes for consistency
        for j = 1:length(modₐ.τ)
            # Amplitude check
            if modₐ.τ[j].α.fixed != modᵦ.τ[j].α.fixed ||
                modₐ.τ[j].α.glob != modᵦ.τ[j].α.glob
                return false
            end
            # Lifetime check
            if modₐ.τ[j].τ.fixed != modᵦ.τ[j].τ.fixed ||
                modₐ.τ[j].τ.glob != modᵦ.τ[j].τ.glob
                return false
            end
        end
        # Now, check anisotropies for consistency
        for j = 1:length(modₐ.θ)
            # Amplitude check
            if modₐ.θ[j].α.fixed != modᵦ.θ[j].α.fixed ||
                modₐ.θ[j].α.glob != modᵦ.θ[j].α.glob
                return false
            end
            # Lifetime check
            if modₐ.θ[j].τ.fixed != modᵦ.θ[j].τ.fixed ||
                modₐ.θ[j].τ.glob != modᵦ.θ[j].τ.glob
                return false
            end
        end
        # Finally, check the others
        if modₐ.Q.glob != modᵦ.Q.glob
            return false
        end
        if modₐ.S.fixed != modᵦ.S.fixed ||
            modₐ.S.glob != modᵦ.S.glob
            return false
        end
        if modₐ.V.fixed != modᵦ.V.fixed ||
            modₐ.V.glob != modᵦ.V.glob
            return false
        end
        if modₐ.Z.fixed != modᵦ.Z.fixed ||
            modₐ.Z.glob != modᵦ.Z.glob
            return false
        end
    end
    return true
end

function isnonnegative(model::TcspcFitModel)
    # Checks whether linear coefficients are nonnegative
    tf = true
    for τ ∈ model.τ
        if τ.α.val < 0
            tf = false
            return tf
        end
    end
    for θ ∈ model.θ
        if θ.α.val < 0
            tf = false
            return tf
        end
    end
    tf = tf && model.S.val >= 0 && model.V.val >= 0 && model.Z.val >= 0
end


text:
push!(LOAD_PATH, abspath(@__DIR__,".."))

using ThinningAndShift ; global const T = ThinningAndShift

using Test
using Plots; using NamedColors; theme(:dark)
using LinearAlgebra
using Distributions

using Random
Random.seed!(0)

##

function crosscov_an(t)
  t < 0 && return 0.0
  τ = antikernels[1].τ
  p = probs[1]*antiprobs[1]
  return - rate_parent*p*exp(-t/τ)/τ#*0.8
end

_ = let dtcov = 20E-3,
  Tcov = 1.0
  (covtimes,covs) =  T.covariance_density_ij(trains...,dtcov,Tcov)
  plt = plot()
  scatter!(plt,covtimes,covs;label="numerical",leg=false,xlabel="time (s)",
    ylabel="cross covariance density")
  ts = range(-Tcov,Tcov;length=200)
  plot!(plt,ts,crosscov_an.(ts);leg=false, linewidth=2 )
  plt
end

##

##

_ = let plt = plot(),
  tlims = (0,2.0),
  msize = 34,
  trplot1 = filter(t-> tlims[1]< t<tlims[2],trains[1])
  trplot2 = filter(t-> tlims[1]< t<tlims[2],trains[2])
  scatter!(plt,trplot1,fill(1.0,length(trplot1));
    marker=:vline,markersize=msize)
  scatter!(plt,trplot2,fill(1.05,length(trplot2));marker=:vline,
    markersize=msize)
  plot!(plt;ylims=(0.95,1.1),leg=false)
end

##
_ = let plt=plot()
end

##

_ = let x=range(0,10.;length=150)
  d = Exponential(inv(3.0))
  plot(x,cdf.(d,x);ylims=(0,1.1),xlims=(0,10),leg=false,linewidth=2)
end

_ = let d = Exponential(0.1)
  cdf(d, 0.8725210739840403 - 0.17 )
end

##

markings = [ [1,2],]
markings_probs = [1.0]
const mySigma = let σ1 = 0.17,σ2 = 0.3, ρ = 0.8
    [ σ1^2 ρ*σ1*σ2   ; ρ*σ1*σ2  σ2^2 ]
end
jitters = [T.JitterDistribution( MultivariateNormal(mySigma)) ,]

const rate_parent = 50.0
const Tend = 20_000.0

gtas_test = T.GTAS(rate_parent,markings,markings_probs,jitters)
trains,t_parent,attr = T.make_samples_with_parent(gtas_test,20_000.0)

##
const dtcov = 0.25
const Tcov = 5.0
timescov,covboth = T.covariance_density_numerical(trains,dtcov,Tcov)

@show rate_parent/dtcov;
@show covboth[1,1,1] ;

std_diff = sqrt(mySigma[1,1]+mySigma[2,2] - 2*mySigma[1,2])
density_an = Normal(0.0,std_diff)
density_an_vals = rate_parent .* pdf.(density_an,timescov)
density_num =  @. 0.5(covboth[:,1,2]+covboth[:,2,1])
plot(timescov,[density_an_vals density_num])

##

const dT = 25.0
times,rates  = T.bin_and_rates(trains,dT,Tend)
@show mean(rates;dims=2)
covnum = cov(rates;dims=2) * dT

rate_parent

##

struct booh{N,NT<:NTuple{N,Number}}
  vals::NT
end

b = booh((1.,4,Float32(3.3)))


text:
facts("Bimodal Cauchy Distributions") do

  context("sample bimodal cauchy with truncation on one side") do
    bc = BlackBoxOptim.bimodal_cauchy(0.65, 0.1, 1.0, 0.1)

    for(i in 1:100)
      v = BlackBoxOptim.sample_bimodal_cauchy(bc; truncateAbove1 = true, truncateBelow0 = false)

      @fact v <= 1.0 => true
      @fact v > 0.0 => true # Very unlikely to be 0.0!!?
    end
  end

  context("sample bimodal cauchy with truncation on both sides") do
    bc = BlackBoxOptim.bimodal_cauchy(0.1, 0.1, 0.95, 0.1)

    for(i in 1:100)
      v = BlackBoxOptim.sample_bimodal_cauchy(bc; truncateAbove1 = true, truncateBelow0 = true)

      @fact v <= 1.0 => true
      @fact v >= 0.0 => true
    end
  end

end


text:
using ContinuousBenchmarks, FileIO, JLD2

# const ldastandata = [
# Dict(
#   "K" => 2,
#   "V" => 5,
#   "M" => 25,
#   "N" => 262,
#   "w" => [4, 3, 5, 4, 3, 3, 3, 3, 3, 4, 5, 3, 4, 4, 5,
# 3, 4, 4, 4, 3, 5, 4, 5, 2, 3, 3, 1, 5, 5, 1, 4,
# 3, 1, 2, 5, 4, 4, 3, 5, 4, 2, 4, 5, 3, 4, 1, 4,
# 4, 3, 2, 1, 2, 1, 2, 2, 2, 1, 2, 2, 3, 1, 2, 2,
# 4, 4, 5, 4, 5, 5, 4, 3, 5, 4, 4, 4, 2, 2, 1, 1,
# 2, 1, 3, 1, 2, 1, 1, 1, 3, 2, 3, 3, 5, 4, 5, 4,
# 3, 5, 4, 2, 2, 2, 1, 3, 2, 1, 3, 1, 3, 1, 1, 2,
# 1, 2, 2, 4, 4, 4, 5, 5, 4, 4, 5, 4, 3, 3, 3, 1,
# 3, 3, 4, 2, 1, 3, 4, 4, 5, 4, 4, 4, 3, 4, 3, 4,
# 5, 1, 2, 1, 3, 2, 1, 1, 2, 3, 3, 3, 3, 4, 1, 4,
# 4, 4, 4, 3, 4, 4, 1, 2, 2, 3, 3, 1, 1, 4, 1, 3,
# 1, 5, 3, 2, 2, 1, 1, 2, 3, 3, 4, 4, 5, 3, 4, 3,
# 1, 5, 5, 5, 3, 3, 4, 5, 3, 3, 3, 2, 3, 1, 3, 3,
# 1, 3, 1, 5, 5, 5, 2, 2, 3, 3, 3, 1, 1, 5, 5, 5,
# 3, 1, 5, 4, 1, 3, 3, 3, 3, 4, 2, 5, 1, 3, 5, 2,
# 5, 5, 2, 1, 3, 3, 5, 3, 5, 3, 3, 5, 1, 2, 2, 1,
# 1, 2, 1, 2, 3, 1, 1],
#   "doc" => [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2,
# 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
# 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4,
# 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
# 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8,
# 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9,
# 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
# 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12,
# 12, 12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 13, 14,
# 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 15, 15, 15,
# 15, 15, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16,
# 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 17, 17,
# 17, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 19, 19,
# 19, 19, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 20,
# 20, 20, 20, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21,
# 22, 22, 22, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23,
# 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 24, 24, 24,
# 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24,
# 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25],
#   "alpha" => [0.5, 0.5],
#   "beta" => [0.2, 0.2, 0.2, 0.2, 0.2]
# )
# ]

const ldastandata = load(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "ldastandata.jld2"))["data"]
# const ldastandata = load(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "ldastandataV20K2M25L100.data"))["data"]
# const ldastandata = load(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "ldastandataV100K5M10L1000.data"))["data"]

# using PyCall
# @pyimport pickle
# f = open("/afs/inf.ed.ac.uk/user/s16/s1672897/fun/20news-sub.pickle", "r")
# data = pickle.loads(pybytes(read(f)))
# const ldastandata = [data]


text:
using CoordinateSystems
using Test

@testset "CoordinateSystems.jl" begin
    # Write your tests here.
end


text:
module VirusLoadCurve

    export fitVLF,
           LogVirusLoadFunction,
           fitVTM,
           LogViralTargetModel,
           VirusLoadData,
           Boxplots
    using OrdinaryDiffEq
    using Optim, LsqFit
    using Plots, LaTeXStrings
    using DelimitedFiles, DataFrames
    using Statistics
    using Turing
    using Printf
    using ProgressMeter

    import Base.summary

    struct VirusLoadData{T<:Real}
        t::Vector{T}
        v::Vector{T}
    end

    @recipe function f(data::VirusLoadData)
        tmin, tmax = extrema(data.t)
        vmin, vmax = extrema(data.v)
        x := data.t
        y := data.v
        seriestype := :scatter
        xaxis := ("Time (days)", (0.0, tmax), font(14))
        yaxis := (L"\log\,V(t)", font(14))
        label := "Data"
        grid := :none
        ()
    end

    include("virus_load_function.jl")
    include("virus_target_model.jl")

end #module


text:
using MakieLayout
using Makie

begin
    scene = Scene(resolution = (1000, 1000));
    screen = display(scene)
    campixel!(scene);

    maingl = GridLayout(scene, alignmode = Outside(30))

    for i in 1:7
        maingl[i, 1:end] = LAxis(scene)
        sleep(0.3)
        maingl[1:end, i + 1] = LAxis(scene)
        sleep(0.3)
    end

    # maingl[2:end-1, 0] = LAxis(scene)

end


text:
## Predicted number of links in a food web using the flexible links model

function links_muphi(mu::Float64, phi::Float64, S::Int64, n::Int64)
    # mu and phi: parameters of the flexible links model
    # S: number of species
    # n: number of simulation

    # simulate total number of links for specific mu and phi

    # parameters alpha and beta
    α = mu*phi
    β = (1.0-mu)*phi

    # beta-binomial distribution
    bb = BetaBinomial(S^2-S+1, α, β)

    # simulate total number of links
    links_rdm = rand(bb, n) .+ S .- 1

    return links_rdm
end

function predict_L(S)
    # S: number of species

    # simulate total number of links for the entire joint posterior distribution

    # empty vector
    links_rdm_post = Vector{Int64}()

    # simulate over the joint posterior distribution
    for i in 1:length(mu_post)
        append!(links_rdm_post, links_muphi(mu_post[i], phi_post[i], S, 1000))
    end

    return shuffle(links_rdm_post)
end

max_stars_repo_name: FrancisBanville/Trophic-METE

max_issues_repo_name: FrancisBanville/Trophic-METE

max_forks_repo_name: FrancisBanville/Trophic-METE


text:
# #############################################################################
# #############################################################################
#Combined computation and analysis of rays in cavity

#This file contains the main function for running the ray solver on a cavity for multiple initial conditions, and storing the desired results on disk.

#   resultsdir = run_rays(bnd::Boundary,idx::RefractiveIndex;kwargs...)
#   This function runs the ray solvers and analysis on a given cavity with boundary <bnd> and index distribution <idx> on a number of initial conditions. The keyword arguments allow control over the following parameters regarding the initial conditions, the solver settings and results recording:
#
#>> Initial conditions
#   pssinitarray::Array{Float64,2} = Array(Float64,0,2)
#       array of initial conditions as (theta,chi) points on Poincare Surface of Section plot. Replaced with autogenerated points if empty.
#   initarray::Array{Float64,2} = Array(Float64,0,3)
#       array of initial conditions as (r,theta,phi) ray parameters. Replaced with equivalent values from <pssinitarray> if empty.
#
#>> Initial conditions autogeneration parameters
#   symmetry::Int64 = 1
#       order of rotational symmetry of cavity to assume
#   thetares::Int64 = 25
#       number of theta values to sample
#   chires::Int64= 25
#       number of chi values to sample
#
#>> Solver parameters
#   maxpathlength::Float64 = 200.0
#       maximum pathlength before halting solvers
#   maxbounces::Int64 = 500
#       maximum number of bounces before halting solvers
#   relativetolerance::Float64 = 1.0e-12
#       relative tolerance limits for solvers
#   absolutetolerance::Float64 = 1.0e-12
#       relative tolerance limits for solvers
#
#>> Results root directory
#   resultsroot::AbstractString = "."
#       root directory to store result subdirectories
#
#>> Results to record
#   record_pssinitarray::Bool = false
#   record_initarray::Bool = true
#   record_initialconditions::Bool = false
#   record_raypath::Bool = false
#   record_bounceindices::Bool = false
#   record_bouncepoints::Bool = true
#   record_cavityimage::Bool = true
#   record_rayimage::Bool = true
#   record_pathlengths::Bool = true
#   record_actions::Bool = true
#   record_modeproperties::Bool = false
#   record_farfield::Bool = false)
#       record value of parameter/result to disk

# #############################################################################
# #############################################################################
#Initiate

#Dependencies
#require("util.jl")
#require("boundary.jl")
#require("refractive_index.jl")
#require("solver.jl")
#require("io.jl")
#require("ray_plot.jl")
#require("ray_analysis.jl")

# #############################################################################
# #############################################################################
#Initial conditions generation

#Use PSS-spanning initial conditions
#Automatically generates an array of (theta,chi) pairs as rows, that covers as much of the Poincare Surface of Section plot up to a specified rotational symmetry.
function gen_pssinitarray(thetares::Int64=25,chires::Int64=25,symmetry::Int64=1)
    #Choose total number of points (see explanation below)
    interval::Int64 = ceil(Int64,sqrt(thetares*chires))
    N::Int64 = interval^2+3
    while gcd(interval,N) != 1
        N += 1
    end

    #Avoid repetition of values by generating N values of theta and N value of
    #chi, then permute them relative to each other so that theta and chi values
    #are matched to form points on a lattice spanning the PSS.
    #generate theta values spanning angle up to cavity rotational symmetry
    theta::Array{Float64,1} = linspace(0,2*pi/symmetry,N+1)[1:N]
    #generate chi values excluding 0 and pi/2
    chi::Array{Float64,1} = linspace(0,(1-1/chires)*pi/2,N+1)[2:end]

    #In order to generate a lattice that is as evenly spread out as possible, can do
    #this by stepping through appropriate intervals in mod <N>. By having N be in the
    #form M^2 + P, where gcd(M,P) = 1, P < M, then lattice is fairly even with the P
    #smallest values spread out as much as possible, in order to sample the low
    #chi behavior at P points.
    return [theta chi[1+mod([(interval*(0:N-1))...],N)]]::Array{Float64,2}
end

#Converts a (theta,chi)-based set initial conditions to a (r,theta,phi)-based one
function pss2raycoord(bnd::Boundary,theta::Float64,chi::Float64)
    r::Float64,alpha::Float64 = rsys(bnd,theta)
    phi::Float64 = alpha + chi
    return (r,theta,phi)
end

#Converts a (r,theta,phi)-based initial conditions array from a given (theta,chi)-based array.
function pss2initarray(bnd::Boundary,pssinitarray::Array{Float64,2})
    count::Int64 = size(pssinitarray,1)
    initA::Array{Float64,2} = Array(Float64,count,3)
    for i=1:count
        initA[i,1],initA[i,2],initA[i,3] = pss2raycoord(bnd,pssinitarray[i,1],pssinitarray[i,2])
    end
    return initA
end

#Automatically generates an array of (r,theta,phi) as rows, that covers as much of the Poincare Surface of Section plot up to a specified rotational symmetry.
function gen_initarray(bnd::Boundary,thetares::Int64=25,chires::Int64=25,symmetry::Int64=1)
    pssinitarray::Array{Float64,2} = gen_pssinitarray(thetares,chires,symmetry)
    return pss2initarray(bnd,pssinitarray)
end

# #############################################################################
# #############################################################################
#Main ray simulation

#Run params hash
#Chunk all the simulation-relevant parameters together
#Then get a hash to use as a checksum
function get_runparamshash( initarray::Array{Float64,2},
                            maxpathlength::Float64,maxbounces::Int64,
                            relativetolerance::Float64,absolutetolerance::Float64,
                            record_pssinitarray::Bool,record_initarray::Bool,
                            record_initialconditions::Bool,record_raypath::Bool,
                            record_bounceindices::Bool,record_bouncepoints::Bool,
                            record_cavityimage::Bool,record_rayimage::Bool,
                            record_pathlengths::Bool,record_actions::Bool,
                            record_modeproperties::Bool,record_farfield::Bool)
    run_params = cell(3)
    run_params[1] = initarray
    run_params[2] = Number[maxpathlength,maxbounces,relativetolerance,absolutetolerance]
    run_params[3] = Bool[record_cavityimage,record_pssinitarray,record_initarray,record_initialconditions,record_raypath,record_bounceindices,record_bouncepoints,record_cavityimage,record_rayimage,record_pathlengths,record_actions,record_modeproperties,record_farfield]
    return hash(run_params)
end

#Actual run
function run_rays(  #Cavity parameters
                    bnd::Boundary,idx::RefractiveIndex;

                    #Initial conditions
                    pssinitarray::Array{Float64,2} = Array(Float64,0,2),
                    initarray::Array{Float64,2} = Array(Float64,0,3),

                    #Initial conditions generator parameters
                    symmetry::Int64 = 1,
                    thetares::Int64 = 25,
                    chires::Int64 = 25,

                    #Solver parameters
                    maxpathlength::Float64 = 200.0,
                    maxbounces::Int64 = 500,
                    relativetolerance::Float64 = 1.0e-12,
                    absolutetolerance::Float64 = 1.0e-12,

                    #Results root directory
                    resultsroot::AbstractString = ".",

                    #Results to record
                    record_pssinitarray::Bool = false,
                    record_initarray::Bool = true,
                    record_initialconditions::Bool = false,
                    record_raypath::Bool = false,
                    record_bounceindices::Bool = false,
                    record_bouncepoints::Bool = true,
                    record_cavityimage::Bool = true,
                    record_rayimage::Bool = true,
                    record_pathlengths::Bool = true,
                    record_actions::Bool = true,
                    record_modeproperties::Bool = false,
                    record_farfield::Bool = false)

    tic()

    #Autogenerate initarray
    if isempty(initarray)
        if isempty(pssinitarray)
            pssinitarray = gen_pssinitarray(thetares,chires,symmetry)
        end
        initarray = pss2initarray(bnd,pssinitarray)
    end

    #Get results directory
    const run_params_hash::UInt64 =
        get_runparamshash(  initarray,
                            maxpathlength,maxbounces,
                            relativetolerance,absolutetolerance,
                            record_pssinitarray,record_initarray,
                            record_initialconditions,record_raypath,
                            record_bounceindices,record_bouncepoints,
                            record_cavityimage,record_rayimage,
                            record_pathlengths,record_actions,
                            record_modeproperties,record_farfield)
    const resultsdir::AbstractString = getresultsdir(run_params_hash,bnd,idx,resultsroot)
    mkpath(resultsdir)

    #Write cavity image and get plot range
    if record_cavityimage
        range::Array{Float64,1} = writecavityimg(resultsdir,bnd,idx)
    else
        #Set plot range to default
        cavityplot = prepareplot()
        plotbnd(bnd,axes=cavityplot)
        range = [plt.axis()...]
    end
    #Initialize ray plot line object to nothing
    rayline = nothing

    #Write initial conditions arrays
    if record_initarray
        writeresults(resultsdir,0,"initarray",initarray)
    end
    if record_pssinitarray
        writeresults(resultsdir,0,"pssinitarray",pssinitarray)
    end

    #Iterate through various initial conditions in parallel
    const num_of_runs::Int64 = size(initarray,1)
    @sync begin
        @parallel for resultid = 1:num_of_runs

            #Run ray tracer
            (raypath::Array{Float64,2},
             bounceindices::Array{Int64,1},
             bouncepoints::Array{Float64,2}) =
                rayevolve_gsl(
                    bnd, idx,
                    initarray[resultid,:][:];
                    tmax = maxpathlength,
                    bouncemax = maxbounces,
                    reltol = relativetolerance,
                    abstol = absolutetolerance)

            #Analyze results
            if record_pathlengths || record_actions
                #Compute
                pathlengths,actions =
                    getpathintegrals(raypath,idx,bounceindices)
            end

            #Write results
            if record_initialconditions
                writeresults(resultsdir,resultid,
                             "initialconditions",initarray[resultid,:][:])
            end
            if record_raypath
                writeresults(resultsdir,resultid,"raypath",raypath)
            end
            if record_bounceindices
                writeresults(resultsdir,resultid,"bounceindices",bounceindices)
            end
            if record_rayimage
                if rayline == nothing
                    rayline = writerayimg(resultsdir,resultid,raypath,range)
                else
                    writerayimg(resultsdir,resultid,raypath,rayline)
                end
            end
            if record_bouncepoints
                writeresults(resultsdir,resultid,"bouncepoints",bouncepoints)
            end
            if record_pathlengths
                writeresults(resultsdir,resultid,"pathlengths",pathlengths)
            end
            if record_actions
                writeresults(resultsdir,resultid,"actions",actions)
            end

            #Mode detected! (not implemented yet)
            if false
                if record_modeproperties

                end
                if record_farfield

                end
            end

            #Display progress
            println("Completed #$resultid of $(num_of_runs) run(s).")
        end
    end

    #End
    plt[:close]()
    println("Cavity analyzed in $(toq()) seconds.")

    return resultsdir
end


text:
using Flux
using FluxModels
using Test

@testset "Vision Models" begin
    include("VisionModel/alexnet.jl")
    include("VisionModel/resnet.jl")
    include("VisionModel/vgg.jl")
    include("VisionModel/googlenet.jl")
end


text:
include("plotHelpers.jl")
Egrid, PAgrid = logrange(10,1000,21), 6:4:90

# themis_lolat = load_resultant_matrix("210429_themis_lolat", "results/themis_lolat_m1pro", "210429_themis_lolat_96600", "setupasrun.conf", 161);
# themis_midlat = load_resultant_matrix("210429_themis_midlat", "results/themis_midlat_m1pro", "210429_themis_midlat_96600", "setupasrun.conf", 161);
# themis_hilat = load_resultant_matrix("210429_themis_hilat", "results/themis_hilat_m1pro", "210429_themis_hilat_96600", "setupasrun.conf", 161);

# themis_1pkt = load_resultant_matrix("210429_themis_1pkt", "results/210429_themis_1pkt2", "210429_themis_1pkt2_96600", "setupasrun.conf", 161);
# themis_30pkt = load_resultant_matrix("210429_themis_30pkt", "results/210429_themis_30pkt2", "210429_themis_30pkt2_96600", "setupasrun.conf", 161);
# themis_100pkt = load_resultant_matrix("210429_themis_100pkt", "results/210429_themis_100pkt2", "210429_themis_100pkt2_96600", "setupasrun.conf", 161);

# mms_lowerlat = load_resultant_matrix("200922_mms_lowerlat", "results/mms_lowerlat_m1pro", "200922_mms_lowerlat_96600", "setupasrun.conf", 161);
# mms_lolat = load_resultant_matrix("200922_mms_lolat", "results/mms_lolat_m1pro", "200922_mms_lolat_96600", "setupasrun.conf", 161);
# mms_midlat = load_resultant_matrix("200922_mms_midlat", "results/mms_midlat_m1pro", "200922_mms_midlat_96600", "setupasrun.conf", 161);
# mms_hilat = load_resultant_matrix("200922_mms_hilat", "results/mms_hilat_m1pro", "200922_mms_hilat_96600", "setupasrun.conf", 161);

# @time @save "results/210429_themis_lolat.jld2" themis_lolat
# @time @save "results/210429_themis_midlat.jld2" themis_midlat
# @time @save "results/210429_themis_hilat.jld2" themis_hilat

# @time @save "results/210429_themis_1pkt.jld2" themis_1pkt
# @time @save "results/210429_themis_30pkt.jld2" themis_30pkt
# @time @save "results/210429_themis_100pkt.jld2" themis_100pkt

# @time @save "results/200922_mms_lowerlat.jld2" mms_lowerlat;
# @time @save "results/200922_mms_lolat.jld2" mms_lolat
# @time @save "results/200922_mms_midlat.jld2" mms_midlat;
# @time @save "results/200922_mms_hilat.jld2" mms_hilat;

@time @load "results/210429_themis_hilat.jld2" themis_hilat;
@time @load "results/210429_themis_midlat.jld2" themis_midlat;
@time @load "results/210429_themis_lolat.jld2" themis_lolat;

@time @load "results/210429_themis_1pkt.jld2" themis_1pkt
@time @load "results/210429_themis_30pkt.jld2" themis_30pkt
@time @load "results/210429_themis_100pkt.jld2" themis_100pkt

@time @load "results/200922_mms_lowerlat.jld2" mms_lowerlat;
@time @load "results/200922_mms_lolat.jld2" mms_lolat;
@time @load "results/200922_mms_midlat.jld2" mms_midlat;
@time @load "results/200922_mms_hilat.jld2" mms_hilat;
@info "Loaded result matrices."

@time equatorial_fluxes_042921 = calc_equatorial_fluxes(themis_hilat, f0_042921);
@time themis_hilat_042921  = export_results("042921_themis_hilat",  calc_prec_flux(themis_hilat, 10,f0_042921,.02));
@time themis_midlat_042921 = export_results("042921_themis_midlat", calc_prec_flux(themis_midlat,10,f0_042921,.04));
@time themis_lolat_042921  = export_results("042921_themis_lolat",  calc_prec_flux(themis_lolat, 10,f0_042921,5*.06));
@time themis_1pkt_042921   = export_results("042921_themis_1pkt",   calc_prec_flux(themis_1pkt, 10,f0_042921,.02));
@time themis_30pkt_042921  = export_results("042921_themis_30pkt",  calc_prec_flux(themis_30pkt,10,f0_042921,.02));
@time themis_100pkt_042921 = export_results("042921_themis_100pkt", calc_prec_flux(themis_100pkt,10,f0_042921,.02));
@time elfin_measurements_042921, elfin_error_042921 = extract_idl_csv("042921_time.csv", "042921_prec.csv",
                                                "042921_precerror.csv", "ebins.csv", # csvs containing ELFIN measurements
                                                DateTime(2021,4,29,3,14,45),DateTime(2021,4,29,3,15,0)) # time to sample from ELFIN measurements
@time @save "results/210429_lat_storage.jld2" equatorial_fluxes_042921 themis_hilat_042921 themis_midlat_042921 themis_lolat_042921;
@time @save "results/210429_pkt_storage.jld2" equatorial_fluxes_042921 themis_1pkt_042921 themis_30pkt_042921 themis_100pkt_042921;

@time equatorial_fluxes_092220 = calc_equatorial_fluxes(mms_lolat, f0_092220);
@time mms_lowerlat_092220  = export_results("092220_mms_lowerlat",  calc_prec_flux(mms_lowerlat, 10,f0_092220,0.06));
@time mms_lolat_092220  = export_results("092220_mms_lolat",  calc_prec_flux(mms_lolat, 10,f0_092220,0.06));
@time mms_midlat_092220 = export_results("092220_mms_midlat", calc_prec_flux(mms_midlat,10,f0_092220,0.02));
@time mms_hilat_092220  = export_results("092220_mms_hilat",  calc_prec_flux(mms_hilat, 10,f0_092220,0.02));
@time @save "results/200922_lat_storage.jld2" equatorial_fluxes_092220 mms_lowerlat_092220 mms_lolat_092220 mms_midlat_092220 mms_hilat_092220;

equatorial_fluxes_092220, elfin_measurements_092220, prec_flux_timeseries_092220 = generate_flux_comparison(mms_midlat,
                                                                10, f0_092220, 1.2,     # timebins, dist_func, whistler occurence rate
                                                                Egrid, PAgrid, # Ebins and PA bins to use
                                                                "092220_time.csv", "092220_prec.csv", "ebins.csv", # csvs containing ELFIN measurements
                                                                DateTime(2020,9,22,9,16,15), DateTime(2020,9,22,9,16,50)) # time to sample from ELFIN measurements

# Plots for 9/22/20
plot(Egrid, equatorial_fluxes_092220, label = "Equatorial Flux", color = bipride_blue, linewidth=2, markershape=:circle);
plot!(Egrid, prec_flux_timeseries_092220[2:end-1], color = bipride_lavender, alpha = .5, label=false, markershape=:x);
plot!(Egrid, prec_flux_timeseries_092220[1], color = bipride_lavender, alpha = .5, label="Modelled Precipitating Flux", markershape=:x);
plot!(elfin_measurements_092220, label = "ELFIN Measured Precipitating Flux", color = bipride_pink, linewidth=3, markershape=:dtriangle);
plot!(ylim =(1e1,1e9), xlim=(50,800), yscale=:log10);
plot!(xlabel=L"\mathrm{Energy\ (keV)}", ylabel=L"\mathrm{Flux\ (1/cm^{2}/s/sr/MeV)}", title=L"\mathrm{Flux\ Comparison\ of\ Precipitating\ Particles\ on\ 09/22/20}",
xtickfontsize=12, ytickfontsize=12, xguidefontsize=16, yguidefontsize=16, legendfontsize=10, titlefontsize=16);
plot1 = plot!(dpi = 300,size=(800,450), margin=3mm, bottom_margin=4mm)
#savefig(plot1, "092220_flux_comparison.pdf")
equatorial_fluxes_102720, elfin_measurements_102720, prec_flux_timeseries_102720 = generate_flux_comparison(themis_hilat,
                                                                10, f0_102720, .05,     # timebins, dist_func, whistler occurence rate
                                                                Egrid, PAgrid, # Ebins and PA bins to use
                                                                "102720_time.csv", "102720_prec.csv", "ebins.csv", # csvs containing ELFIN measurements
                                                                DateTime(2020,10,27,10,34,7),DateTime(2020,10,27,10,34,40)) # time to sample from ELFIN measurements

# Plots for 10/27/20
plot(Egrid, equatorial_fluxes_102720, label = "Equatorial Flux", color = bipride_blue, linewidth=2, markershape=:circle);
plot!(Egrid, prec_flux_timeseries_102720[2:end-1], color = bipride_lavender, alpha = .5, label=false, markershape=:x);
plot!(Egrid, prec_flux_timeseries_102720[1], color = bipride_lavender, alpha = .5, label="Modelled Precipitating Flux", markershape=:x);
plot!(elfin_measurements_102720, label = "ELFIN Measured Precipitating Flux", color = bipride_pink, linewidth=3, markershape=:dtriangle);
plot!(ylim =(1e4,1e11), xlim=(50,800), yscale=:log10);
plot!(xlabel=L"\mathrm{Energy\ (keV)}", ylabel=L"\mathrm{Flux\ (1/cm^{2}/s/sr/MeV)}", title=L"\mathrm{Flux\ Comparison\ of\ Precipitating\ Particles\ on\ 10/27/20}",
xtickfontsize=12, ytickfontsize=12, xguidefontsize=16, yguidefontsize=16, legendfontsize=10, titlefontsize=16);
plot2 = plot!(dpi = 300,size=(800,450), margin=3mm, bottom_margin=4mm)
#savefig(plot2, "102720_flux_comparison.pdf")

equatorial_fluxes_042921, elfin_measurements_042921, prec_flux_timeseries_042921 = generate_flux_comparison(themis_midlat,
                                                                10, f0_042921, .05,     # timebins, dist_func, whistler occurence rate
                                                                Egrid, PAgrid, # Ebins and PA bins to use
                                                                "042921_time.csv", "042921_prec.csv", "ebins.csv", # csvs containing ELFIN measurements
                                                                DateTime(2021,4,29,3,14,45),DateTime(2021,4,29,3,15,0)) # time to sample from ELFIN measurements

# Plots for 04/29/21
plot(Egrid, equatorial_fluxes_042921, label = "Equatorial Flux", color = bipride_blue, linewidth=2, markershape=:circle);
plot!(Egrid, prec_flux_timeseries_042921[2:end-1], color = bipride_lavender, alpha = .5, label=false, markershape=:x);
plot!(Egrid, prec_flux_timeseries_042921[1], color = bipride_lavender, alpha = .5, label="Modelled Precipitating Flux", markershape=:x);
plot!(elfin_measurements_042921, label = "ELFIN Measured Precipitating Flux", color = bipride_pink, linewidth=3, markershape=:dtriangle);
plot!(ylim =(1e2,1e9), xlim=(50,800), yscale=:log10);
plot!(xlabel=L"\mathrm{Energy\ (keV)}", ylabel=L"\mathrm{Flux\ (1/cm^{2}/s/sr/MeV)}", title=L"\mathrm{Flux\ Comparison\ of\ Precipitating\ Particles\ on\ 4/29/21}",
xtickfontsize=12, ytickfontsize=12, xguidefontsize=16, yguidefontsize=16, legendfontsize=10, titlefontsize=16);
plot3 = plot!(dpi = 300,size=(800,450), margin=3mm, bottom_margin=4mm)
# savefig(plot3, "042921_flux_comparison.pdf")

mms_hilat = export_results("200922_mms_hilat", equatorial_fluxes_092220, prec_flux_timeseries_092220)
plot(Egrid, mms_hilat.precipitating_fluxes_mean, yerror=(mms_hilat.precipitating_fluxes_minus, mms_hilat.precipitating_fluxes_plus), ylim =(1e2,1e9), xlim=(50,800), yscale=:log10)
plot!(Egrid, mms_hilat.equatorial_fluxes)
@save "200922_data_storage.jld2" mms_midlat

themis_midlat = export_results("210429_themis_midlat", equatorial_fluxes_042921, prec_flux_timeseries_042921)
plot(Egrid, themis_midlat.precipitating_fluxes_mean, yerror=(themis_midlat.precipitating_fluxes_minus, themis_midlat.precipitating_fluxes_plus), ylim =(1e2,1e9), xlim=(50,800), yscale=:log10)
plot!(Egrid, themis_midlat.equatorial_fluxes)
@save "210429_data_storage.jld2" themis_hilat_100wav themis_hilat_nowav themis_hilat themis_midlat themis_lolat

trapped_fluxes = generate_trapped_psd(f0_092220, 1.0)
simulated_ratio = [prec_flux_timeseries[i]./((trapped_fluxes[i]).-prec_flux_timeseries[i]) for i in 1:length(trapped_fluxes)]
plot(Egrid, simulated_ratio, xlim=(50,1000),legend=false, xscale=:log10)
elfin_prec = extract_idl_csv("092220_time.csv", "092220_prec.csv", "ebins.csv", DateTime(2020,9,22,9,16,15), DateTime(2020,9,22,9,16,50))
elfin_trap = extract_idl_csv("092220_time.csv", "092220_perp.csv", "ebins.csv", DateTime(2020,9,22,9,16,15), DateTime(2020,9,22,9,16,50))
plot!(elfin_prec[1], elfin_prec[2] ./ elfin_trap[2], xlim=(50,1000), yscale=:log10)

trapped_fluxes = generate_trapped_psd(f0_102720, 0.1)
simulated_ratio = [prec_flux_timeseries[i]./((trapped_fluxes[i]./45.).-prec_flux_timeseries[i]) for i in 1:length(trapped_fluxes)]
plot(Egrid, simulated_ratio, xlim=(50,1000),legend=false, xscale=:log10)
elfin_prec = extract_idl_csv("102720_time.csv", "102720_prec.csv", "ebins.csv", DateTime(2020,10,27,10,34,7),DateTime(2020,10,27,10,34,40))
elfin_trap = extract_idl_csv("102720_time.csv", "102720_perp.csv", "ebins.csv", DateTime(2020,10,27,10,34,7),DateTime(2020,10,27,10,34,40))
plot!(elfin_prec[1], elfin_prec[2] ./ elfin_trap[2], xlim=(50,1000), yscale=:log10)

plot(Egrid,binned_psd_prec_timeseries[14].*psd_0i.*Egrid./1000,xlim=(Egrid[1], Egrid[end]), xscale=:log10)
extract_idl_csv("time.csv", "prec.csv", "ebins.csv", DateTime(2020,9,22,9,16,30),DateTime(2020,9,22,9,16,35))
plot(ebins,flux, xlim = (ebins[1], 1e3), xscale =:log10, ylim = (1, 1e8), yscale = :log10)

plot!(xlim=(10,1000), ylim=(1e2,1e11), yscale=:log10)

animateNewPSD("21000_animation.gif", Egrid, PAgrid)

plot_numLostParticles = plot(lostParticles[:,1]*Re/(c), lostParticles[:,2], legend=false,
    xlim = [0,maximum(maximum(allT))*Re/(c)], ylim = [0,numParticles],
    title = "Number of Lost Particles", xlabel="Time (s)", ylabel="Number of Particles")
savefig(plot_numLostParticles, string("particleLosses.png"))

plot3 = plot(
    (1000*allT*Re)./(c), allPA, title = "Pitch Angles",
    legend=false,xlabel="Time (ms)", ylabel="Pitch Angle (deg)",
    xlim = [0,50], ylim = [0,90]);
plot4 = plot(
    (1000*allT*Re)./(c), allE, title = "Energies",
    legend=false,xlabel="Time (ms)", ylabel="Energy (keV)",
    xlim = [0,50], ylim = [0,3000]);
bigplot = plot(plot3,plot4, dpi = 300, layout = (2,1))

animatePAD("63000PADevolution.gif", 2000, 5)
animateESD("4200ESDevolution.gif", 4200, 50)
animatePSD("63000PSDevolution2.gif", 5, 50, 450)

# makes
allPrecip, indexArray = precipitatingParticles(themis_100pkt, 10);
animatePrecipitatingParticles("test_precipanimation.gif", themis_100pkt, 10, 0.015)

initial, final = 1, 5
Egrid = logrange(10,1000,20)
PAgrid = 2:4:90
f, psd_init, psd_final = recalcDistFunc(Ematrix,PAmatrix,initial,final,f0_102720, Egrid, PAgrid, 1.0);
checkDistFunc(f, psd_init, psd_final, initial, final, Egrid, PAgrid)
savefig(plot_numLostParticles, string("particleLosses.png"))

# function f0(E, alpha)
#     Estar = 30; # kev
#     return (E/Estar * (sin(deg2rad(alpha)) - sin(deg2rad(lossConeAngle))))
# end

# check ICs
# plot(1:numParticles,sort(Ematrix[1,:]))
# plot(1:numParticles,sort(PAmatrix[1,:]))

# energyDistributionofPrecipiatingParticles
# energy and time phase space dist at bin of time
# can use large bin of time
# 1 simulation time = 0.1 s
# each bin = 10 simulation time

# trajectory of conjunctions for mms and elfin
# plot whistler wave spectrum
# elfin data

# plot trajectory
trajectoryChecking(45)
trajectoryTracing([7300,7309,6900,6905,6910,6918,6927], 10)

animDec = 1; # make a png for animation every 10 points
animScale = 10; # i.e. animscale = 10 means every 10 seconds in animation is 1 second of simulation time (increase for longer animation)
anim = @animate for i in eachindex(tVec)
    p1 = plot!(tVec[1:i], PAmatrix[:,j][1:i], label = "")
    scatter!((tVec[i], PAmatrix[:,j][i]), color = 1, label = "")

    # annotate!(20, 0.1*maxFraction, text("t = $(round(tVec[indexArray[i]]*Re*L/(c),digits=4)) s"), :left)
end every animDec
savename = string("results/","trajectories")
gif(anim, savename, fps = (length(tVec)/animDec)/(animScale*endTime*Re*L/(c)))

max_stars_repo_name: ethantsai/nlwhistlers

max_issues_repo_name: ethantsai/nlwhistlers

max_forks_repo_name: ethantsai/nlwhistlers


text:
using ProgressMeter

"""
    dimension_test(numSamples,numFeatures,blockSize)

    # Description
    Throws an error if the dimensions of the particular problem instance do not
    make sense.
"""
function dimension_test(numSamples::Integer, numFeatures::Integer, blockSize::Integer)
    # PCA of a single vector makes no sense...
    if numSamples == 1
        error("StreamingPCA.jl: <<Dataset consists of a single sample.>>")
    end

    # PCA of samples with a single feature...
    if numFeatures == 1
        error("StreamingPCA.jl: <<Dataset consists of samples with a single feature.>>")
    end

    # Require that the size of the dataset be divisible by the block size
    if (numSamples % blockSize) != 0
        error("StreamingPCA.jl:  <<Dataset size ($(numSamples)) must be divisible by block size ($(blockSize)).>>")
    end
end

"""
    block_ranges(numSamples::Integer,blockSize::Integer)

    # Description
    Return a UnitRange{Int64} type for the dataset indicies of the
    specified batch.
"""
function block_range(numSamples::Integer,blockSize::Integer, batchNumber::Integer)
    loIdx = blockSize*(batchNumber-1) + 1
    hiIdx = loIdx + blockSize - 1

    return loIdx:hiIdx
end

"""
    Stream over samples, update the Q/R as needed if we have enough data to
    warrant a full batch
"""
# function streaming_pca_process_sample!(x::Array{Float64},
#                                        Q::Array{Float64}, R::Array{Float64},
#                                        S::Array{Float64}, sBQ::Array{Float64},
#                                        batchSize::Integer,sampleNumber::Integer)

#     gemm!('T','N',1/batchSize,x,sBQ[sample,m:],1.0,S)
# end

"""
    streaming_pca_process_batch!(batch::Array{Float64},
                                 Q::Array{Float64},R::Array{Float64})

    # Description
    Given a current estimate of the QR decomposition of the dataset (the
    passed QR values), update this estimate given a new batch of data. This
    function will modify the factorization in place. To do this in place, we
    will have to use `qrfact!` and the factorization object instead of `qr`.
"""
function streaming_pca_process_batch!(batch::Array{Float64},
                                      Q::Array{Float64},R::Array{Float64})
    numFeatures,batchSize = size(batch)
    rank = size(Q,2)
    S = zeros(numFeatures,rank)
    scale = 1/batchSize

    sBQ = scale.*(batch*Q)

    for sample in 1:batchSize
        gemm!('T','N',1.0,batch[sample,:],sBQ[sample,:],1.0,S)
    end

    Qloc,Rloc = qr(S)
end

"""
   streaming_pca(X::Array{Float64}, blockSize::Integer)

   # Description
   Given a dataset X, where samples are arranged as rows, and a valid block
   size, return the principal components of the dataset as estimated by the
   'Block-Stochastic Orthogonal Iteration' given in Alg. 1 of
   (Mitliagkas, Caramanis & Jain, 2013).
"""
function streaming_pca(X::Array{Float64}, blockSize::Integer; rank::Integer = 1)
    numSamples, numFeatures = size(X)
    dimension_test(numSamples,numFeatures,blockSize)

    numBlocks = convert(Integer,numSamples / blockSize)

    # Initial Parameters
    Q, R = qr(randn(numFeatures,rank))
    scale = 1 / blockSize

    S = []  # Declaring for access outside block
    @showprogress 0.001 "Streaming batches..." 20 for block in 1:numBlocks
        S = zeros(numFeatures,rank)
        # Get indicies for this block
        blockRange = block_range(numSamples,blockSize,block)
        sXQ = scale.*(X[blockRange,:]*Q)

        for sample in blockRange
            # Over samples, the core procedure is (assuming row
            # vector samples)...
            #     S = (a'aQ)/s + (b'bQ)/s + (c'cQ)/s + (d'dQ)/s + ...
            # So, we should be able to put together this operation with
            # just gemm, now that we have precomputed all (aQ)/s, (bQ)/s, ...
            gemm!('T','N',1.0,X[sample,:],sXQ[sample-blockRange[1]+1,:],1.0,S)
        end

        Q, R = qr(S)
    end

    # Now get the data projection
    # XQ = (Q\(X'))'
    XQ = X*Q;
    # What is the order of the components (via norm)?
    sq_pc_energy = vec(sum(XQ.^2,1));
    # Sort components by energy
    Q = Q[:,sortperm(sq_pc_energy;rev=true)]
    sort!(sq_pc_energy;rev=true)

    return Q,XQ,sq_pc_energy
end


text:
module ModelingToolkit

using DiffEqBase, Distributed
using StaticArrays, LinearAlgebra, SparseArrays, LabelledArrays
using Latexify, Unitful, ArrayInterface
using MacroTools
using UnPack: @unpack
using DiffEqJump
using DataStructures: OrderedDict, OrderedSet
using SpecialFunctions, NaNMath
using Base.Threads
import MacroTools: splitdef, combinedef, postwalk, striplines
import GeneralizedGenerated
using DocStringExtensions
using Base: RefValue

using RecursiveArrayTools

import SymbolicUtils
import SymbolicUtils: to_symbolic, FnType, @rule, Rewriters

import LightGraphs: SimpleDiGraph, add_edge!

import TreeViews

using Requires

"""
$(TYPEDEF)

Base type for a symbolic expression.
"""
abstract type Expression <: Number end

"""
$(TYPEDEF)

TODO
"""
abstract type AbstractSystem end
abstract type AbstractODESystem <: AbstractSystem end

Base.promote_rule(::Type{<:Number},::Type{<:Expression}) = Expression
Base.zero(::Type{<:Expression}) = Constant(0)
Base.zero(::Expression) = Constant(0)
Base.one(::Type{<:Expression}) = Constant(1)
Base.one(::Expression) = Constant(1)
Base.oneunit(::Expression) = Constant(1)
Base.oneunit(::Type{<:Expression}) = Constant(1)

"""
$(TYPEDSIGNATURES)

Get the set of independent variables for the given system.
"""
function independent_variables end

"""
$(TYPEDSIGNATURES)

Get the set of states for the given system.
"""
function states end

"""
$(TYPEDSIGNATURES)

Get the set of parameters variables for the given system.
"""
function parameters end

include("variables.jl")
include("context_dsl.jl")
include("operations.jl")
include("differentials.jl")

function Base.convert(::Type{Variable},x::Operation)
    if x.op isa Variable
        x.op
    elseif x.op isa Differential
        var = x.args[1].op
        rename(var,Symbol(var.name,:ˍ,x.op.x))
    else
        throw(error("This Operation is not a Variable"))
    end
end

include("equations.jl")
include("function_registration.jl")
include("simplify.jl")
include("utils.jl")
include("linearity.jl")
include("direct.jl")
include("domains.jl")

include("systems/abstractsystem.jl")

include("systems/diffeqs/odesystem.jl")
include("systems/diffeqs/sdesystem.jl")
include("systems/diffeqs/abstractodesystem.jl")
include("systems/diffeqs/first_order_transform.jl")
include("systems/diffeqs/modelingtoolkitize.jl")
include("systems/diffeqs/validation.jl")

include("systems/jumps/jumpsystem.jl")

include("systems/nonlinear/nonlinearsystem.jl")

include("systems/optimization/optimizationsystem.jl")

include("systems/pde/pdesystem.jl")

include("systems/reaction/reactionsystem.jl")
include("systems/dependency_graphs.jl")

include("latexify_recipes.jl")
include("build_function.jl")

export ODESystem, ODEFunction, ODEFunctionExpr, ODEProblemExpr
export SDESystem, SDEFunction, SDEFunctionExpr, SDESystemExpr
export JumpSystem
export ODEProblem, SDEProblem
export NonlinearProblem, NonlinearProblemExpr
export OptimizationProblem, OptimizationProblemExpr
export SteadyStateProblem, SteadyStateProblemExpr
export JumpProblem, DiscreteProblem
export NonlinearSystem, OptimizationSystem
export ode_order_lowering
export PDESystem
export Reaction, ReactionSystem, ismassaction
export Differential, expand_derivatives, @derivatives
export IntervalDomain, ProductDomain, ⊗, CircleDomain
export Equation, ConstrainedEquation
export Operation, Expression, Variable
export independent_variable, states, parameters, equations

export calculate_jacobian, generate_jacobian, generate_function
export calculate_tgrad, generate_tgrad
export calculate_gradient, generate_gradient
export calculate_factorized_W, generate_factorized_W
export calculate_hessian, generate_hessian
export calculate_massmatrix, generate_diffusion_function

export BipartiteGraph, equation_dependencies, variable_dependencies
export eqeq_dependencies, varvar_dependencies
export asgraph, asdigraph

export simplified_expr, rename, get_variables
export simplify, substitute
export build_function
export @register
export modelingtoolkitize
export @variables, @parameters

const HAS_DAGGER = Ref{Bool}(false)
function __init__()
    @require Dagger="d58978e5-989f-55fb-8d15-ea34adc7bf54" include("dagger.jl")
end

end # module


text:
using BenchmarkTools, SpecialFunctions, FiniteDifferences, ForwardDiff, StaticArrays
include("shared.jl")

const FDFAST(fn, x) = (fn(x+h)-fn(x))/h
const FD2 = central_fdm(2, 1)

dbesselk_fdfast(v, x) = FDFAST(_v->besselk(_v, x), v)
dbesselk_fd2(v, x)    = FD2(_v->besselk(_v, x), v)
dbesselk_ad(v, x)     = ForwardDiff.derivative(_v->_besselk(_v, x), v)

#=
# And note zero allocations for the AD version. Pretty good. And faster by
# significant margins---like a factor of five---than even the most reckless FD.
for (v, x) in Iterators.product((0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 2.25, 3.0, 4.25),
                                 (1e-4, 0.25, 2.5, 5.0, 7.5, 12.0, 25.0))
  println("(v,x) = ($v, $x):")
  print("FDFAST:")
  @btime dbesselk_fdfast($v, $x)
  print("FD2:   ")
  @btime dbesselk_fd2($v, $x)
  print("AD:    ")
  @btime dbesselk_ad($v, $x)
  print("\n\n")
end
=#

# matern function:
const oo = @SVector ones(2)
const zz = @SVector zeros(2)

@inline pv(v) = @SVector [1.1, 1.1, v]
function matern_d3_ad(v)
  ForwardDiff.derivative(_v->BesselK.matern(oo, zz, pv(_v)), v)
end

function matern_d3_fd(v)
  FDFAST(_v->BesselK.matern(oo, zz, pv(_v)), v)
end


text:
module SimpleANOVA

#=

Workflow:
results = anova() # get table for all effects and detect interactions
plot(results)     # assess whether interactions are trivial or large

# if interactions are trivial, can immediately do posthoc test of choice on the independent main effects
tukey(results)

# if interactions are nontrivial, do subanovas with one less factor, at each level of the removed factor, with α/k for k levels
subanova(results)

# finally, do pairwise (e.g. tukey) tests among levels of one factor
tukey(results, [1,2])
=#

using CategoricalArrays
using Combinatorics
using Distributions
using InvertedIndices
using Requires
using Statistics

include("data/AnovaEffect.jl")
include("data/AnovaValue.jl")
include("data/AnovaFactor.jl")
include("data/AnovaResult.jl")
include("data/AnovaData.jl")
include("data/AnovaContrastResult.jl")
include("data/FactorType.jl")
include("anova.jl")
include("pretests.jl")
include("contrasts.jl")

function __init__()
    @require DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
        include("anova_dataframes.jl")
        include("pretests_dataframes.jl")
    end
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("anova_plots.jl")
end

export anova, ftest, plot, levene, contrast, differencecontrasts, repeatedcontrasts, simplecontrasts
export AnovaEffect, AnovaValue, AnovaFactor, AnovaResult, AnovaData, AnovaContrastResult, AnovaContrastResults
export FactorType, fixed, random, nested, subject, block

end