text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	hydro_inter_period_linkage(EP::Model, inputs::Dict)
This function creates variables and constraints enabling modeling of long duration storage resources when modeling representative time periods.
**Storage inventory balance at beginning of each representative period**
The constraints in this section are used to approximate the behavior of long-duration energy storage technologies when approximating annual grid operations by modeling operations over representative periods. Previously, the state of charge balance for storage (as defined in ```storage_all()```) assumed that state of charge at the beginning and end of each representative period has to be the same. In other words, the amount of energy built up or consumed by storage technology $o$ in zone $z$ over the representative period $m$, $\Delta Q_{o,z,m} = 0$. This assumption implicitly excludes the possibility of transferring energy from one representative period to the other which could be cost-optimal when the capital cost of energy storage capacity is relatively small. To model long-duration energy storage using representative periods, we replace the state of charge equation, such that the first term on the right hand side accounts for change in storage inventory associated with representative period $m$ ($\Delta Q_{o,z,m}$), which could be positive (net accumulation) or negative (net reduction).
```math
\begin{aligned}
& \Gamma_{o,z,(m-1)\times \tau^{period}+1 } =\left(1-\eta_{o,z}^{loss}\right)\times \left(\Gamma_{o,z,m\times \tau^{period}} -\Delta Q_{o,z,m}\right) -  \\
& \frac{1}{\eta_{o,z}^{discharge}}\Theta_{o,z,(m-1)\times \tau^{period}+1} + \eta_{o,z}^{charge}\Pi_{o,z,(m-1)\times \tau^{period}+1} \quad \forall o \in \mathcal{O}^{LDES}, z \in \mathcal{Z}, m \in \mathcal{M}
\end{aligned}
```
By definition $\mathcal{T}^{start}=\{\left(m-1\right) \times \tau^{period}+1 | m \in \mathcal{M}\}$, which implies that this constraint is defined for all values of $t \in T^{start}$.
**Storage inventory change input periods**
We need additional variables and constraints to approximate energy exchange between representative periods, while accounting for their chronological occurence in the original input time series data and the possibility that two representative periods may not be adjacent to each other (see Figure below). To implement this, we introduce a new variable $Q_{o,z, n}$ that models inventory of storage technology $o \in O$ in zone $z$ in each input period $n \in \mathcal{N}$. Additionally we define a function mapping, $f: n \rightarrow m$, that uniquely maps each input period $n$ to its corresponding representative period $m$. This mapping is available as an output of the process used to identify representative periods (E.g. k-means clustering [Mallapragada et al., 2018](https://www.sciencedirect.com/science/article/pii/S0360544218315238?casa_token=I-6GVNMtAVIAAAAA:G8LFXFqXxRGrXHtrzmiIGm02BusIUmm83zKh8xf1BXY81-dTnA9p2YI1NnGuzlYBXsxK12by)).
![Modeling inter-period energy exchange via long-duration storage when using representative period temporal resolution to approximate annual grid operations](assets/LDES_approach.png)
*Figure. Modeling inter-period energy exchange via long-duration storage when using representative period temporal resolution to approximate annual grid operations*
The following two equations define the storage inventory at the beginning of each input period $n+1$ as the sum of storage inventory at begining of previous input period $n$ plus change in storage inventory for that period. The latter is approximated by the change in storage inventory in the corresponding representative period, identified per the mapping $f(n)$.  The second constraint relates the storage level of the last input period, $|N|$, with the storage level at the beginning of the first input period. Finally, if the input period is also a representative period, then a third constraint enforces that initial storage level estimated by the intra-period storage balance constraint should equal the initial storage level estimated from the inter-period storage balance constraints. Note that $|N|$ refers to the last modeled period.
```math
\begin{aligned}
& Q_{o,z,n+1} = Q_{o,z,n} + \Delta Q_{o,z,f(n)}
\quad \forall  o \in \mathcal{O}^{LDES}, z \in \mathcal{Z}, n \in \mathcal{N}\setminus\{|N|\}
\end{aligned}
```
```math
\begin{aligned}
& Q_{o,z,1} = Q_{o,z,|N|} + \Delta Q_{o,z,f(|N|)}
\quad \forall  o \in \mathcal{O}^{LDES}, z \in \mathcal{Z}, n = |N|
\end{aligned}
```
```math
\begin{aligned}
& Q_{o,z,n} =\Gamma_{o,z,f(n)\times \tau^{period}} - \Delta Q_{o,z,m}
\quad \forall  o \in \mathcal{O}^{LDES}, z \in \mathcal{Z}, n \in   \mathcal{N}^{rep},
\end{aligned}
```
Finally, the next constraint enforces that the initial storage level for each input period $n$ must be less than the installed energy capacity limit. This constraint ensures that installed energy storage capacity is consistent with the state of charge during both the operational time periods $t$ during each sample period $m$ as well as at the start of each chronologically ordered input period $n$ in the full annual time series.
```math
\begin{aligned}
    Q_{o,z,n} \leq \Delta^{total, energy}_{o,z}
\quad \forall n \in \mathcal{N}, o \in \mathcal{O}^{LDES}
\end{aligned}
```
"""
function hydro_inter_period_linkage(EP::Model, inputs::Dict)

	println("Long Duration Storage Module for Hydro Reservoir")

	dfGen = inputs["dfGen"]

	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones
	REP_PERIOD = inputs["REP_PERIOD"]     # Number of representative periods

	STOR_HYDRO_LONG_DURATION = inputs["STOR_HYDRO_LONG_DURATION"]

	START_SUBPERIODS = inputs["START_SUBPERIODS"]

	hours_per_subperiod = inputs["hours_per_subperiod"] #total number of hours per subperiod

	dfPeriodMap = inputs["Period_Map"] # Dataframe that maps modeled periods to representative periods
	NPeriods = size(inputs["Period_Map"])[1] # Number of modeled periods

	MODELED_PERIODS_INDEX = 1:NPeriods
	REP_PERIODS_INDEX = MODELED_PERIODS_INDEX[dfPeriodMap[!,:Rep_Period] .== MODELED_PERIODS_INDEX]

	### Variables ###

	# Variables to define inter-period energy transferred between modeled periods

	# State of charge of storage at beginning of each modeled period n
	@variable(EP, vSOC_HYDROw[y in STOR_HYDRO_LONG_DURATION, n in MODELED_PERIODS_INDEX] >= 0)

	# Build up in storage inventory over each representative period w
	# Build up inventory can be positive or negative
	@variable(EP, vdSOC_HYDRO[y in STOR_HYDRO_LONG_DURATION, w=1:REP_PERIOD])

	### Constraints ###

	# Links last time step with first time step, ensuring position in hour 1 is within eligible change from final hour position
	# Modified initial state of storage for long-duration storage - initialize wth value carried over from last period
	# Alternative to cSoCBalStart constraint which is included when not modeling operations wrapping and long duration storage
	# Note: tw_min = hours_per_subperiod*(w-1)+1; tw_max = hours_per_subperiod*w
	@constraint(EP, cSoCBalLongDurationStorageStart_H[w=1:REP_PERIOD, y in STOR_HYDRO_LONG_DURATION],
				    EP[:vS_HYDRO][y,hours_per_subperiod*(w-1)+1] == (EP[:vS_HYDRO][y,hours_per_subperiod*w]-vdSOC_HYDRO[y,w])-(1/dfGen[!,:Eff_Down][y]*EP[:vP][y,hours_per_subperiod*(w-1)+1])-EP[:vSPILL][y,hours_per_subperiod*(w-1)+1]+inputs["pP_Max"][y,hours_per_subperiod*(w-1)+1]*EP[:eTotalCap][y])

	# Storage at beginning of period w = storage at beginning of period w-1 + storage built up in period w (after n representative periods)
	## Multiply storage build up term from prior period with corresponding weight
	@constraint(EP, cSoCBalLongDurationStorageInterior_H[y in STOR_HYDRO_LONG_DURATION, r in MODELED_PERIODS_INDEX[1:(end-1)]],
					vSOC_HYDROw[y,r+1] == vSOC_HYDROw[y,r] + vdSOC_HYDRO[y,dfPeriodMap[!,:Rep_Period_Index][r]])

	## Last period is linked to first period
	@constraint(EP, cSoCBalLongDurationStorageEnd_H[y in STOR_HYDRO_LONG_DURATION, r in MODELED_PERIODS_INDEX[end]],
					vSOC_HYDROw[y,1] == vSOC_HYDROw[y,r] + vdSOC_HYDRO[y,dfPeriodMap[!,:Rep_Period_Index][r]])

	# Storage at beginning of each modeled period cannot exceed installed energy capacity
	@constraint(EP, cSoCBalLongDurationStorageUpper_H[y in STOR_HYDRO_LONG_DURATION, r in MODELED_PERIODS_INDEX],
					vSOC_HYDROw[y,r] <= dfGen[!,:Hydro_Energy_to_Power_Ratio][y]*EP[:eTotalCap][y])

	# Initial storage level for representative periods must also adhere to sub-period storage inventory balance
	# Initial storage = Final storage - change in storage inventory across representative period
	@constraint(EP, cSoCBalLongDurationStorageSub_H[y in STOR_HYDRO_LONG_DURATION, r in REP_PERIODS_INDEX],
					vSOC_HYDROw[y,r] == EP[:vS_HYDRO][y,hours_per_subperiod*dfPeriodMap[!,:Rep_Period_Index][r]] - vdSOC_HYDRO[y,dfPeriodMap[!,:Rep_Period_Index][r]])

	return EP
end

ext: jl

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
@testset "Test DIRECTIONAL.S2 on short predefined sequence" begin
    a = [1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0]

    corr = mean(Directional.s2(a, 1; len = 12))
    @test corr ≈ [7, 4, 3, 2, 1, 1, 2, 3, 2, 2, 1, 0] ./ [12 - x for x in 0:11]

    corr = mean(Directional.s2(a, 0; len = 12))
    @test corr ≈ [5, 2, 2, 1, 1, 1, 1, 1, 0, 1, 0, 0] ./ [12 - x for x in 0:11]

    corr = mean(Directional.s2(a, 1; len = 12, periodic = true))
    @test corr ≈ [7, 4, 4, 4, 3, 4, 4, 4, 3, 4, 4, 4] ./ 12

    corr = mean(Directional.s2(a, 0; len = 12, periodic = true))
    @test corr ≈ [5, 2, 2, 2, 1, 2, 2, 2, 1, 2, 2, 2] ./ 12
end

@testset "Test DIRECTIONAL.L2 on short predefined sequences" begin
    a = [1, 0, 1, 0, 1, 1, 1]

    corr = mean(Directional.l2(a, 1; len = 7))
    @test corr ≈ [5, 2, 1, 0, 0, 0, 0] ./ [7 - x for x in 0:6]

    corr = mean(Directional.l2(a, 0; len = 7))
    @test corr ≈ [2, 0, 0, 0, 0, 0, 0] ./ [7 - x for x in 0:6]

    corr = mean(Directional.l2(a, 1; len = 7, periodic = true))
    @test corr ≈ [5, 3, 2, 1, 0, 0, 0] ./ 7

    a = [1, 1, 0, 1, 0]
    @test Directional.l2(a, 1; len = 5).success ==
        Directional.l2(a, 1; len = 5, periodic = true).success

    a = [0, 1, 1, 0, 1]
    @test Directional.l2(a, 1; len = 5).success ==
        Directional.l2(a, 1; len = 5, periodic = true).success
end

@testset "Test on array with all elements equal to 1" begin
    a = fill(1, (50, 50, 50))
    for p in (false, true)
        for func in (Directional.s2, Directional.l2, Directional.c2)
            @test all(x -> x ≈ 1,
                      mean(func(a, 1; periodic = p, directions = known_directions)))
        end
    end

    @test length(pore_size(a, 0).weights) == 0
end


text:
import Base.show

###############
# Bytevector
###############

# Long form, as in display(h) or REPL: h enter
function Base.show(io::IO, ::MIME"text/plain", v::Bytevector)
    ioc = IOContext(io)
    n = length(v.data)
    println(ioc, "  ", n, "-element Bytevector{UInt8}:")
    print(ioc, "\t")
    if get(ioc, :color, false)
        for i = 1:n
            printstyled(ioc, string(i)[end], "¹²³⁴ ⁵⁶⁷⁸", color = :green, bold = false)
        end
    else
        for i = 1:n
            print(ioc, string(i)[end], "¹²³⁴ ⁵⁶⁷⁸")
        end
    end
    println(ioc)
    print(ioc, "\t")
    for i = 1:n
        bs = string(v.data[i], base=2, pad = 8)
        print(ioc, " ", bs[1:4], " ", bs[5:8])
    end
    print(ioc," ")
end

# Short form, as in print(stdout, h)
function Base.show(io::IO, v::Bytevector)
    ioc = IOContext(io)
    n = length(v.data)
    print(ioc, "\t")
    if get(ioc, :color, false)
        for i = 1:n
            printstyled(ioc, string(i)[end], "¹²³⁴ ⁵⁶⁷⁸", color = :green, bold = false)
        end
    else
        for i = 1:n
            print(ioc, string(i)[end], "¹²³⁴ ⁵⁶⁷⁸")
        end
    end
    println(ioc)
    print(ioc, "\t")
    for i = 1:n
        bs = string(v.data[i], base=2, pad = 8)
        print(ioc, " ", bs[1:4], " ", bs[5:8])
    end
    print(ioc," ")
end
# Standard output is fine

# The default Juno / Atom display works nicely with standard output
Base.show(io::IO, ::MIME"application/prs.juno.inline", v::Bytevector) = Base.show(io, v)

###############
# DevState
###############

# Long form, as in display(h) or REPL: h enter
function Base.show(io::IO, ::MIME"text/plain", ds::DevState)
    print(io, "DevState:")
    show(io, MIME("text/plain"), ds.values)
    print(io, "\n\t timestamp = ", ds.timestamp, ", proceed = ", ds.proceed)
end

# Short form, as in print(stdout, h)
function Base.show(io::IO, ds::DevState)
    print(io, "DevState(")
    show(io, ds.values)
    print(io, ", ", ds.timestamp, ", ", ds.proceed, ")")
end

# The default Juno / Atom display works nicely with standard output
Base.show(io::IO, ::MIME"application/prs.juno.inline", ds::DevState) = Base.show(io, v)

###############
# Utilty
###############
function bytestring(b::UInt8)
    io = IOBuffer()
    bs = string(b, base=2, pad = 8)
    print(io, bs[1:4], " ", bs[5:8])
    String(take!(io))
end


text:
using GNSSTools
using JLD
using Statistics
using PyPlot
pygui(true)
include("ch4_data_plot_functions.jl")
include("ch4_simulation_function.jl")

directory = "/home/sjbilardi/Dropbox/Apps/Overleaf/MastersThesis/"
folder = "figures/l1ca_leo_sat_sim/"

# Iridium
a = Rₑ+780*1000;  # meters
plane_num = 6;
sat_per_plane = 11;
incl = 86;
ΔΩ = 360/plane_num/2;
Δf_per_plane = 360/sat_per_plane/2;  # degrees
eop = get_eop();
CN0 = 45
h_parms = [1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15] ./ 16.3676e6^2  # Rakon IT5300B TCXO
trackresults, data, doppler_curve, doppler_t, sat_elevations,
sat_azimuth = demo(a, plane_num, sat_per_plane; ΔΩ=ΔΩ, eop=eop, t_length=10*60,
                   include_phase_noise=true, include_databits=false,
                   h_parms=h_parms, CN0=CN0, showplot=false);

save(string(directory, folder, "ch4_l1ca_leo_sat_data_no_databits.jld"),
     "trackresults", trackresults,
     "data", data,
     "doppler_curve", doppler_curve,
     "doppler_t", doppler_t,
     "sat_elevations", sat_elevations,
     "sat_azimuth", sat_azimuth,
     "a", a,
     "plane_num", plane_num,
     "sat_per_plane", sat_per_plane,
     "incl", incl,
     "d_omega", ΔΩ,
     "Δf_per_plane", Δf_per_plane,
     "h_parms", h_parms,
     "CN0", CN0)

trackresults, data, doppler_curve, doppler_t, sat_elevations, sat_azimuth,
a, plane_num, sat_per_plane, incl, ΔΩ,Δf_per_plane = load(string(directory, folder, "ch4_l1ca_leo_sat_data_no_databits.jld"),
                                                          "trackresults",
                                                          "data",
                                                          "doppler_curve",
                                                          "doppler_t",
                                                          "sat_elevations",
                                                          "sat_azimuth",
                                                          "a",
                                                          "plane_num",
                                                          "sat_per_plane",
                                                          "incl",
                                                          "d_omega",
                                                          "Δf_per_plane")

file_name = string(directory, folder, "ch4_l1ca_leo_sat_example_no_databits.png")
phi0 = π/4
doppler_rate = diff(doppler_curve)
doppler_rate = [doppler_rate[1]; doppler_rate]
n0 = trackresults.data_init_code_chip
cn0 = 45
Tsys = 535
B = 2*l1ca_chipping_rate
A = sqrt(2*GNSSTools.k*Tsys)*10^(cn0/20)
σ = sqrt(GNSSTools.k*B*Tsys)
PS = floor(Int, 5e6*1e-3)*A^2
PN = σ^2
snr = calc_snr(PS, PN)
plot_example_track_results(trackresults; f_if=data.f_if,
                           code_length=l1ca_code_length,
                           save_to_file=file_name,
                           bottom=0.075, top=0.95, left=0.125, right=0.875,
                           hspace=0.75, wspace=0.35, marker_size=1,
                           line_width=0.8, use_ylims=false)#,
                         #   truth_t=t, truth_doppler=doppler,
                         #   truth_doppler_rate=doppler_rate,
                         #   truth_code_phase=code_phase,
                         #   truth_snr=snr, truth_phi=phi)


text:
import Base: convert
import Base: ==

export CalcFactor

"""
$TYPEDEF

User factor interface method for computing the residual values of factors.

Notes
- Also see #467 on API consolidation

```juila
function (cf::CalcFactor{<:LinearRelative})(res::AbstractVector{<:Real}, z, xi, xj)
  cf.metadata.variablelist
  cf.metadata.targetvariable
  cf.metadata.usercache
  generic on-manifold residual function

  return distance(z, distance(xj, xi))
end
```

DevNotes
- Follow the Github project in IIF to better consolidate CCW FMD CPT CF CFM

Related

[`CalcFactorMahalanobis`](@ref), [`CommonConvWrapper`](@ref), [`FactorMetadata`](@ref), [`ConvPerThread`](@ref)
"""
struct CalcFactor{T <: AbstractFactor, M, P <: Union{<:Tuple,Nothing,AbstractVector}, X}
  # the interface compliant user object functor containing the data and logic
  factor::T
  # the metadata to be passed to the user residual function
  metadata::M
  # what is the sample (particle) id for which the residual is being calculated
  _sampleIdx::Int
  # legacy support when concerned with how many measurement tuple elements are used by user
  _measCount::Int
  # legacy suport for measurement sample values of old functor residual functions
  _legacyMeas::P
  # legacy support for variable values old functor residual functions
  _legacyParams::X
end

"""
    $TYPEDEF

Internal parametric extension to [`CalcFactor`](@ref) used for buffering measurement and calculating Mahalanobis distance

Related

[`CalcFactor`](@ref)
"""
struct CalcFactorMahalanobis{CF<:CalcFactor, S, N}
  calcfactor!::CF
  varOrder::Vector{Symbol}
  meas::NTuple{N, <:AbstractVector{Float64}}
  iΣ::NTuple{N, Matrix{Float64}}
  specialAlg::S
end

abstract type _AbstractThreadModel end

"""
$(TYPEDEF)
"""
struct SingleThreaded <: _AbstractThreadModel
end
"""
$(TYPEDEF)
"""
struct MultiThreaded <: _AbstractThreadModel
end

"""
$(TYPEDEF)

Notes
- type-stable `usercache`, see https://github.com/JuliaRobotics/IncrementalInference.jl/issues/783#issuecomment-665080114

DevNotes
- TODO why not just a NamedTuple? Perhaps part of #467
- Follow the Github project in IIF to better consolidate CCW FMD CPT CF CFM
- TODO standardize -- #927, #1025, #784, #692, #640
- TODO make immutable #825
"""
mutable struct FactorMetadata{FV<:AbstractVector{<:DFGVariable},
                              VL<:AbstractVector{Symbol},
                              AR, # <:AbstractVector{<:AbstractArray}
                              CD}
  # full list of Vector{DFGVariable} connected to the factor
  fullvariables::FV # Vector{<:DFGVariable}
  #TODO full variable can perhaps replace this
  variablelist::VL # Vector{Symbol}
  # TODO consolidate, same as ARR used in CCW,
  arrRef::AR # Vector{Matrix{Float64}}
  # label of which variable is being solved for
  solvefor::Symbol
  # for type specific user data, see (? #784)
  cachedata::CD
end

"""
$(TYPEDEF)

DevNotes
- FIXME consolidate with CalcFactor
- TODO consolidate with CCW, FMd, CalcFactor
- TODO consider renaming `.p` to `.decisionDims`
  - TODO `.decisionDims::DD where DD <: Union{<:AbstractVector{Int},Colon}` -- ensure type stability
- TODO figure out if we want static parameter THRID
- TODO make static params {XDIM, ZDIM, P}
- TODO make immutable
"""
mutable struct ConvPerThread{R,F<:FactorMetadata,P}
  thrid_::Int
  # the actual particle being solved at this moment
  particleidx::Int
  # additional/optional data passed to user function
  factormetadata::F
  # subsection indices to select which params should be used for this hypothesis evaluation
  activehypo::Vector{Int}
  # slight numerical perturbation for degenerate solver cases such as division by zero
  perturb::Vector{Float64}
  # working memory location for optimization routines on target decision variables
  X::Vector{P}
  # working memory to store residual for optimization routines
  res::R # was Vector{Float64}
end

"""
$(TYPEDEF)

Main factor memory container used during inference operations -- i.e. values specific to one complete convolution operation

Notes
- CCW does not get serialized / persisted
- At writing, the assumption is there is just one CCW per factor
- Any multithreaded design needs to happens as sub-constainers inside CCW or otherwise, to carry separate memory.
- Since #467, `CalcFactor` is the only type 'seen by user' during `getSample` or function residual calculations `(cf::CalcFactor{<:MyFactor})`, s.t. `MyFactor <: AbstractRelative___`
- There also exists a `CalcFactorMahalanobis` for parameteric computations using as much of the same mechanics as possible.

DevNotes
- Follow the Github project in IIF to better consolidate CCW FMD CPT CF CFM

Related

[`CalcFactor`](@ref), [`CalcFactorMahalanobis`](@ref), [`FactorMetadata`](@ref), [`ConvPerThread`](@ref)
"""
mutable struct CommonConvWrapper{ T<:AbstractFactor,
                                  H<:Union{Nothing, Distributions.Categorical},
                                  C<:Union{Nothing, Vector{Int}},
                                  NTP <: NamedTuple,
                                  G,
                                  MT} <: FactorOperationalMemory
  #
  ### Values consistent across all threads during approx convolution
  usrfnc!::T # user factor / function
  # general setup
  xDim::Int
  zDim::Int
  # is this a partial constraint as defined by the existance of factor field `.partial::Tuple`
  partial::Bool
  # multi hypothesis settings
  hypotheses::H
  # categorical to select which hypothesis is being considered during convolution operation
  certainhypo::C
  nullhypo::Float64
  # parameters passed to each hypothesis evaluation event on user function, #1321
  params::NTP
  # which index is being solved for in params?
  varidx::Int
  # FIXME make type stable, JT should now be type stable if rest works
  #   user defined measurement values for each approxConv operation
  measurement::Vector{MT}
  # TODO refactor and deprecate this old approach, Union{Type{SingleThreaded}, Type{MultiThreaded}}
  threadmodel::Type{<:_AbstractThreadModel}
  # FIXME, deprecate for only `readonly` and build CalcFactor objects on stack instead
    ## will be obsolete: particular convolution computation values per particle idx (varies by thread)
  cpt::Vector{<:ConvPerThread}
  # inflationSpread
  inflation::Float64
  # Which dimensions does this factor influence.  Sensitive (mutable) to both which 'solvefor index' variable and whether the factor is partial dimension
  partialDims::Vector{<:Integer}
  # variable types for points in params
  vartypes::Vector{DataType}
  # experimental feature to embed gradient calcs with ccw
  _gradients::G
end

#


text:
using GadgetGalaxies
using Documenter

DocMeta.setdocmeta!(GadgetGalaxies, :DocTestSetup, :(using GadgetGalaxies); recursive=true)

makedocs(;
    modules=[GadgetGalaxies],
    authors="Lucas Valenzuela <lucasvalenzuela@users.noreply.github.com>",
    repo="https://github.com/lucasvalenzuela/GadgetGalaxies.jl/blob/{commit}{path}#{line}",
    sitename="GadgetGalaxies.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://lucasvalenzuela.github.io/GadgetGalaxies.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Kinematic Parameters" => "kinematic_parameters.md",
        "API" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/lucasvalenzuela/GadgetGalaxies.jl",
)


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

using Markdown
using InteractiveUtils

# ╔═╡ c7cd922c-1b76-11eb-2860-2bdefdcbf8b2
begin
	import Pkg
	Pkg.activate(".")

	using SimpleValueGraphs

	# For creating reproducable random matrices
	using Random: MersenneTwister
	using SparseArrays

	# We are going to test our new packages with some functions
	# from Graphs.jl and GraphsRecipes
	using Graphs
	using GraphRecipes, Plots
	pyplot() # gr backend has issues showing edge labels so we use pyplot instead
end

# ╔═╡ 9ce46c16-2085-11eb-076f-256b70803020
md"""
# Creating a new value graph type

This notebook show how one create a new graph type that subclasses `AbstractValGraph`.

The new graph type will be a wrapper type around a matrix `A`, so that the graph has
a directed edge from vertex `i` to vertex `j` with weight `A[i, i]` if that entry is not zero.

We start by importing a few packages that we are going to use later:
"""

# ╔═╡ daf7b178-2087-11eb-23cb-6d4f29edc176
md"""
The definition of `AbstractValGraph` is
```julia
abstract type AbstractValGraph{V<:Integer, V_VALS, E_VALS, G_VALS} <: AbstractGraph{V} end
```
where
* `V` is the type of the vertices
* `V_VALS` are the types of the vertex values
* `E_VALS` are the types of the edge values
* `G_VALS` are the types of the edge values

As this type of `Graphs.AbstractGraph`, it can also be used with Graphs functions as long as we correctly implement the `AbstractGraph` interface.

In our case the vertex type will be `Int` as this type is also used for indexing rows and columns in a matrix.

We don't have any graph or vertex values, so use the empty tuple type for `V_VALS` and `G_VALS`.

For the edge values we have a single type (the type of the values in our matrix). We also want to give the name `:weight` to these values. This can be done by using the named tuple type `NamedTuple{(:weight,), Tuple{T}}` where `T` is the type of the matrix values.

Then our definition is:
"""

# ╔═╡ 9259d936-1b79-11eb-15cc-abaa483a07c0
struct GraphView{T, M<:AbstractMatrix{<:T}} <: SimpleValueGraphs.AbstractValGraph{
		Int,
		Tuple{},
		NamedTuple{(:weight,), Tuple{T}},
		Tuple{}
	}

	matrix::M
end

# ╔═╡ ec4e6c94-2082-11eb-3d82-d799a91863a7
md"""

## Implementing the AbstractValGraph interface

To be able to do anything meaningful with that graph, we have to implement the interface for `Graphs.AbstractGraph` as well as the interface for `SimpleGraphs.AbstractValGraph`. Luckily, SimpleValueGraphs provides already some sensible defaults (that can be overriden for performance or other reasons) for a lot of Graphs.jl functions so that what we will have to implement is less than is usually required by Graphs.jl. Nevertheless we need to implement the following functions:

* `nv(::GraphView)`
* `is_directed(::Type{<:GraphView})`
* `has_edge(::GraphView, s, d)`
* `zero(::Type{<:GraphView{T, M}) where {T, M}`

And because our graphs has edge values, we also need to implement
* `get_edgeval(::GraphView, s, d, :)`

"""

# ╔═╡ 3f533e8e-208e-11eb-090a-6543f23256a2
md"""
##### nv
`nv(g)` should return the number of vertices in our graph. We could make sure in the constructor that the wrapped matrix is square, but instead we just the minimum of the number of rows and columns as the number of vertices so that any additional values in the matrix are simply ignored.
"""

# ╔═╡ 5e162db8-208e-11eb-09e0-5f14ac2f631c
SimpleValueGraphs.nv(g::GraphView) = minimum(size(g.matrix))

# ╔═╡ 5923245a-208e-11eb-0dd9-0fc7c7344c84
md"""
##### is_directed

`is_directed` returns whether the graph is directed or not. This should be a property of the graph type itself, regardless of the values in that graph. As our matrix is not necessarily symmetric (although we could ensure that with a check in the constructor), our graph will be directed.
"""

# ╔═╡ bb9187f8-1b7f-11eb-15c5-d7884b4cca33
SimpleValueGraphs.is_directed(::Type{<:GraphView}) = true

# ╔═╡ a3b59cc2-208f-11eb-32d0-238d11964bf9
md"""
##### has_edge
`has_edge(g, u, v)` should return true if there is an edge between the vertices `u` and `v`. This should be the case when the entry in the wrapped matrix is zero. We can check that by using the `iszero` function. This also means that our graph type does not work for value types that do not have `iszero` defined.

We also verify that `u` and `v` are actually vertices of our graph. For this we use the `vertices` function that alreay has a default implementation through the implementation of `nv`.
"""

# ╔═╡ e28e02c8-1b7f-11eb-2949-1fc7e27f6cc2
function SimpleValueGraphs.has_edge(g::GraphView, u, v)
	u ∉ vertices(g) && return false
	v ∉ vertices(g) && return false

	return !iszero(g.matrix[u, v])
end

# ╔═╡ 10913788-2091-11eb-31a9-8dc219364439
md"""
##### get_edgeval
`get_edgeval(g, u, v, key::Integer)` should return the edge value associated with the edge between `u` and `v` for a  specific key. We only have to implement this method for integer keys. As we only have a single value for each edge, we accept all values for the key - a more sophisticated approach would be to throw an error if the key is not correct.

In our case we do not explicetly verify that the edge exist altought we could do so and throw a specialised error message. But this is not part of the interface and functions that deal with value graphs should not relay on that.
"""

# ╔═╡ 5fce373e-56a0-11eb-015f-e39d395121d4
SimpleValueGraphs.get_edgeval(g::GraphView, u, v, key::Integer) = weight=g.matrix[u, v]

# ╔═╡ ed96451c-5596-11eb-1c39-8b28c0e661af
md"""
##### zero
This function is a bit of an anomaly and should in my opinion not be part of the Graphs.jl interface. It is also rarely used in Graphs.jl, so it might not be a very big issue to omit it. Nevertheless we implement it here for the sake of completeness.

`zero(G)` should create a graph with zero vertices, given a graph type. We do this here by trying to create a matrix of size (0, 0) and the correct matrix type, and then wrap a `GraphView` around it:
"""

# ╔═╡ d594af1c-5690-11eb-3b9d-4f3e1f6384ee
function SimpleValueGraphs.zero(::Type{<:GraphView{T, M}}) where {T, M}
	matrix = convert(M, Matrix{T}(undef, 0, 0))
	return GraphView{T, M}(matrix)
end

# ╔═╡ 455072ea-5692-11eb-1a6b-23f3cfac6a16
zero(GraphView{String, Matrix{String}})

# ╔═╡ a36d032a-5692-11eb-1b38-3772dea5f067
zero(GraphView{Float64, SparseMatrixCSC{Float64}})

# ╔═╡ 45463dee-56a0-11eb-2181-8f0089913ef5
md"## Using our new graph type"

# ╔═╡ 7757542e-2092-11eb-20b4-ff4234b49d32
md"""
Now we are ready to create `GraphView` of a matrix and verify that it indeed works.

We start by creating a sparse 5x5 matrix where each entry is _not_ zero with probability 0.3. For better reproducibility we used a `MersenneTwister` random generator with a fixed seed.

"""

# ╔═╡ 9ef5fe50-1b7d-11eb-1ad4-115e1adeeebe
A = sprand(MersenneTwister(1), Float64, 5, 5, 0.3)

# ╔═╡ 339348dc-2093-11eb-3896-d1e25cbfa433
md"""
Then we wrap this matrix in our `GraphView`. As `A` has 11 non-zero entries, our directed graph should also have 11 edges.
"""

# ╔═╡ bee02bdc-1b7d-11eb-2f87-bb568e3e6b7c
g = GraphView(A)

# ╔═╡ 8d4cc722-2093-11eb-0b98-8741fa511c84
md"Let's verify that our graph has the correct number of edges"

# ╔═╡ bc09149c-2093-11eb-3d18-1b90f05311a8
ne(g)

# ╔═╡ dad3ce14-2093-11eb-076d-e1c784c74c8c
md"We can also look at the list of edges in our graph:"

# ╔═╡ e8b78692-2093-11eb-233b-63061135f033
edges(g) |> collect

# ╔═╡ ecaf3358-2093-11eb-33a9-a1404f9571b0
md"""
This seems to work as it should. By using the `weights` function, we can also get a matrix of the weights in our graph:
"""

# ╔═╡ 4c9dc1f8-2094-11eb-0660-75c6a339a717
weights(g)

# ╔═╡ 598a9a44-2094-11eb-0518-9f2b052d0140
md"We can verify that these are the same values as in our wrapped matrix:"

# ╔═╡ 6f99a4f6-2094-11eb-3987-3d079af43bd8
weights(g) == A

# ╔═╡ 7a02edde-56a0-11eb-28cd-f15320b018ed
md"""
We should also quickly check that `zero` works and then forget about it:
"""

# ╔═╡ 90c970ba-56a0-11eb-039b-45777dda3357
zero(typeof(g))

# ╔═╡ 75a1a2a6-2094-11eb-2989-d5fa04fb8a33
md"""
### Plotting

Using [GraphRecipes.jl](https://github.com/JuliaPlots/GraphRecipes.jl) we can now also plot our graph:
"""

# ╔═╡ a6f47b56-2094-11eb-05ce-a336e6dd028b
graphplot(g)

# ╔═╡ d035e7d4-2094-11eb-2121-3b675521b9db
md"""
Let's add some vertex labels and the eddge weights so that is easier to verify that this corresponds to our matrix:
"""

# ╔═╡ ff2d75f6-2045-11eb-3d9c-218b3f5adfb3
graphplot(g,
	names=vertices(g),
	edgelabel=trunc.(weights(g), digits=2), # truncate weights for better readability
	nodecolor=:lightgreen,
	nodesize=0.2
)

# ╔═╡ 9a2272c8-2096-11eb-349d-91f34702dd1f
md"""
### Shortest paths

Let's use some functions from Graphs.jl to find the shortest path from vertex `2` to `5`:
"""

# ╔═╡ 1147357a-2095-11eb-224e-61ad831db786
path = enumerate_paths(dijkstra_shortest_paths(g, 2), 5)

# ╔═╡ f24f1466-2098-11eb-11ad-479e850f798d
md"""
We want to visualize this path with GraphRecipes. For that we create a `Dict` that assigns to each edge (a tuple of source and destination of that edge) a color.
"""

# ╔═╡ be361882-2097-11eb-3003-a3503eefca35
begin
	edgecolor = Dict()
	# Set all edgecolors to black
	for e in edges(g)
		edgecolor[(src(e), dst(e))] = :black
	end
	# Set the edgecolor for edges on the path to red
	for i in 1:length(path)-1
		edgecolor[(path[i], path[i+1])] = :red
	end
end

# ╔═╡ a33449d8-2099-11eb-3740-a91b14a285f6
md"""
Now we can pass `edgecolor` to the plotting function:
"""

# ╔═╡ 4abed890-2099-11eb-2287-6b1dce2dd4f2
graphplot(g,
	names=vertices(g),
	edgelabel=trunc.(weights(g), digits=2),
	edgecolor=edgecolor,
	nodecolor=:lightgreen,
	nodesize=0.2
)

# ╔═╡ e8acb5ca-1b7d-11eb-3f4e-4f4dfb631517

# ╔═╡ 4daf3fde-1bcb-11eb-2eec-ebe2c0c0f11d

# ╔═╡ a3d03dde-1bcb-11eb-37ff-b558b6ecde1a

# ╔═╡ 6f14c1fa-1bcb-11eb-1ea3-5547297c6fa0

# ╔═╡ Cell order:
# ╟─9ce46c16-2085-11eb-076f-256b70803020
# ╠═c7cd922c-1b76-11eb-2860-2bdefdcbf8b2
# ╟─daf7b178-2087-11eb-23cb-6d4f29edc176
# ╠═9259d936-1b79-11eb-15cc-abaa483a07c0
# ╠═ec4e6c94-2082-11eb-3d82-d799a91863a7
# ╟─3f533e8e-208e-11eb-090a-6543f23256a2
# ╠═5e162db8-208e-11eb-09e0-5f14ac2f631c
# ╟─5923245a-208e-11eb-0dd9-0fc7c7344c84
# ╠═bb9187f8-1b7f-11eb-15c5-d7884b4cca33
# ╟─a3b59cc2-208f-11eb-32d0-238d11964bf9
# ╠═e28e02c8-1b7f-11eb-2949-1fc7e27f6cc2
# ╟─10913788-2091-11eb-31a9-8dc219364439
# ╠═5fce373e-56a0-11eb-015f-e39d395121d4
# ╟─ed96451c-5596-11eb-1c39-8b28c0e661af
# ╠═d594af1c-5690-11eb-3b9d-4f3e1f6384ee
# ╠═455072ea-5692-11eb-1a6b-23f3cfac6a16
# ╠═a36d032a-5692-11eb-1b38-3772dea5f067
# ╟─45463dee-56a0-11eb-2181-8f0089913ef5
# ╟─7757542e-2092-11eb-20b4-ff4234b49d32
# ╠═9ef5fe50-1b7d-11eb-1ad4-115e1adeeebe
# ╟─339348dc-2093-11eb-3896-d1e25cbfa433
# ╠═bee02bdc-1b7d-11eb-2f87-bb568e3e6b7c
# ╟─8d4cc722-2093-11eb-0b98-8741fa511c84
# ╠═bc09149c-2093-11eb-3d18-1b90f05311a8
# ╟─dad3ce14-2093-11eb-076d-e1c784c74c8c
# ╠═e8b78692-2093-11eb-233b-63061135f033
# ╟─ecaf3358-2093-11eb-33a9-a1404f9571b0
# ╠═4c9dc1f8-2094-11eb-0660-75c6a339a717
# ╟─598a9a44-2094-11eb-0518-9f2b052d0140
# ╠═6f99a4f6-2094-11eb-3987-3d079af43bd8
# ╟─7a02edde-56a0-11eb-28cd-f15320b018ed
# ╠═90c970ba-56a0-11eb-039b-45777dda3357
# ╟─75a1a2a6-2094-11eb-2989-d5fa04fb8a33
# ╠═a6f47b56-2094-11eb-05ce-a336e6dd028b
# ╟─d035e7d4-2094-11eb-2121-3b675521b9db
# ╠═ff2d75f6-2045-11eb-3d9c-218b3f5adfb3
# ╟─9a2272c8-2096-11eb-349d-91f34702dd1f
# ╠═1147357a-2095-11eb-224e-61ad831db786
# ╟─f24f1466-2098-11eb-11ad-479e850f798d
# ╠═be361882-2097-11eb-3003-a3503eefca35
# ╟─a33449d8-2099-11eb-3740-a91b14a285f6
# ╠═4abed890-2099-11eb-2287-6b1dce2dd4f2
# ╟─e8acb5ca-1b7d-11eb-3f4e-4f4dfb631517
# ╟─4daf3fde-1bcb-11eb-2eec-ebe2c0c0f11d
# ╟─a3d03dde-1bcb-11eb-37ff-b558b6ecde1a
# ╟─6f14c1fa-1bcb-11eb-1ea3-5547297c6fa0


text:
using MaximumSubMatrix
using Test

@testset "MaximumSubMatrix.jl" begin
    include("histogram.jl")
    include("submatrix.jl")
end


text:
# MathOptInterface objective

MOI.get(model::Optimizer, ::MOI.ObjectiveSense) = model.sense

# Objective definition.
function add_objective!(model::Optimizer, objective::MOI.ScalarQuadraticFunction)
    eidx = rhp_add_equ(model.ctx)
    # We parse the expression passed in arguments.
    qvidx1, qvidx2, qcoefs = canonical_quadratic_reduction(objective)
    lvidx, lcoefs = canonical_linear_reduction(objective)
    rhp_equ_add_quadratic(model.ctx, eidx, qvidx1, qvidx2, qcoefs)
    avar = _ensure_avar(model)
    rhp_avar_set(avar, lvidx)
    rhp_equ_add_linear_chk(model.ctx, eidx, avar, lcoefs)
    reshop_set_rhs(model.ctx, eidx, objective.constant)
    rhp_set_objeqn(model.ctx, eidx)
    return
end

function add_objective!(model::Optimizer, objective::MOI.ScalarAffineFunction)
    eidx = rhp_add_equ(model.ctx)
    lvidx, lcoefs = canonical_linear_reduction(objective)
    avar = _ensure_avar(model)
    rhp_avar_set(avar, lvidx)
    rhp_equ_add_linear(model.ctx, eidx, avar, lcoefs)
    reshop_set_rhs(model.ctx, eidx, objective.constant)
    rhp_set_objeqn(model.ctx, eidx)
    return
end

function add_objective!(model::Optimizer, var::MOI.SingleVariable)
    check_inbounds(model, var)
    # TODO(Xhub) because of ovf, we need to always add an equation
#    rhp_set_objvar(model.ctx, var.variable.value - 1)
    eidx = rhp_add_equ(model.ctx)
    avar = _ensure_avar(model)
    rhp_avar_set(avar, var.variable.value - 1)
    rhp_equ_add_linear(model.ctx, eidx, avar, [1.,])
    rhp_set_objeqn(model.ctx, eidx)
    return
end

function MOI.set(model::Optimizer, ::MOI.ObjectiveFunction,
                 func::Union{MOI.SingleVariable, MOI.ScalarAffineFunction,
                             MOI.ScalarQuadraticFunction})
    add_objective!(model, func)
    return
end

function MOI.set(model::Optimizer, ::MOI.ObjectiveSense,
                 sense::MOI.OptimizationSense)
    model.sense = sense
    reshop_set_objsense(model.ctx, sense)
end


text:
# calculates column max of a large matrix stored in HDFS
# julia -p 3
# include("test_maxval.jl")

using Elly
using HadoopBlocks

const INP = "hdfs://tan@localhost:9000/maxvalinp.csv"

if myid() == 1
    function gendata(M, N)
        open(HDFSFile(INP), "w") do f
            for m in 1:M
                write(f, join(map(string, rand(N)), ","), "\n")
            end
        end
    end
    gendata(10^5, 5)
end

# split script here and uncomment below lines to run data generation and computation separately
#using Elly
#using HadoopBlocks
#const INP = "hdfs://tan@localhost:9000/maxvalinp.csv"

findrow(r::HdfsBlockReader, iter_status) = HadoopBlocks.find_rec(r, iter_status, Vector)

maprow(rec) = [tuple([parse(Float64, x) for x in rec]...)]

function collectrow(results, rec)
    isempty(rec) && (return results)
    (results == nothing) && (return rec)
    tuple([max(results[x],rec[x]) for x in 1:length(results)]...)
end

function reducerow(reduced, results...)
    (nothing == reduced) && (reduced = zeros(Float64, length(results[1])))
    for res in results
        (nothing == res) && continue
        for x in 1:length(res)
            reduced[x] = max(reduced[x], res[x])
        end
    end
    reduced
end

function wait_results(j_mon)
    loopstatus = true
    while(loopstatus)
        sleep(5)
        jstatus,jstatusinfo = status(j_mon,true)
        ((jstatus == "error") || (jstatus == "complete")) && (loopstatus = false)
        (jstatus == "running") && println("$(j_mon): $(jstatusinfo)% complete...")
    end
    wait(j_mon)
    println("time taken (total time, wait time, run time): $(times(j_mon))")
    println("")
end

j = dmapreduce(MRHdfsFileInput([INP], findrow), maprow, collectrow, reducerow)
wait_results(j)
println(results(j))

max_stars_repo_name: tanmaykm/HadoopBlocks

max_issues_repo_name: tanmaykm/HadoopBlocks

max_forks_repo_name: tanmaykm/HadoopBlocks


text:
"""
    estimate_volatility(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts),
                        method::Symbol = :two_scales_volatility;
                        time_grid::Union{Missing,Dict} = missing ,
                        fixed_spacing::Union{Missing,Dict,<:Real} = missing,
                        use_all_obs::Bool = false, rough_guess_number_of_intervals::Integer = 5,
                        num_grids::Real = default_num_grids(ts))

This is a convenience wrapper for the two volatility estimation techniques included in this package.
### General Inputs
* `ts` - The tick data.
* `assets` - What assets from ts that you want to estimate the covariance for.
* `method` - The method can be `:simple_volatility` (for the simple volatility method) or `:two_scales_volatility` (for the two scales volatility method)
#### Inputs only used in `:simple_volatility` method.
* `time_grid` - The grid with which to calculate returns. If missing one is generated with a fixed spacing (if that is provided) or a default spacing.
* `fixed_spacing` - A spacing used to calculate a time grid. Not used if a `time_grid` is input or if `use_all_obs = true`.
* `use_all_obs` - Use all observations to estimate volatilities. Not used if a `time_grid` is provided.
* `rough_guess_number_of_intervals` - A rough number of intervals to calculate a default spacing. Not used if a `time_grid` or `fixed_spacing` is provided or if `use_all_obs = true`.
#### Inputs only used in `:two_scales_volatility` method.
* `num_grids` - Number of grids used in order in two scales estimation.
### Returns
* A `Dict` with estimated volatilities for each asset.
"""
function estimate_volatility(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts),
                             method::Symbol = :two_scales_volatility;
                             time_grid::Union{Missing,Dict} = missing ,
                             fixed_spacing::Union{Missing,Dict,<:Real} = missing,
                             use_all_obs::Bool = false, rough_guess_number_of_intervals::Integer = 5,
                             num_grids::Real = default_num_grids(ts))
    if method == :simple_volatility
        return simple_volatility(ts, assets;
                                   time_grid = time_grid, fixed_spacing = fixed_spacing,
                                   use_all_obs = use_all_obs, rough_guess_number_of_intervals = rough_guess_number_of_intervals)
    elseif method == :two_scales_volatility
        return two_scales_volatility(ts, assets; num_grids = num_grids)[1]
    else
        error("The volatilty method chosen must be either :simple_volatility or :two_scales_volatility")
    end
end

"""
    estimate_microstructure_noise(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts);
                                  num_grids::Real = default_num_grids(ts))

This estimates microstructure noise with the two_scales_volatility method.
### Inputs
* `ts` - The tick data.
* `assets` - What assets from ts that you want to estimate the covariance for.
* `num_grids` - Number of grids used in order in two scales estimation.
### Returns
* A `Dict` with estimated microstructure noise variances for each asset.
"""
function estimate_microstructure_noise(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts);
                             num_grids::Real = default_num_grids(ts))
    return two_scales_volatility(ts, assets; num_grids = num_grids)[2]
end

"""
    estimate_covariance(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts),
                        method::Symbol = :preaveraged_covariance;
                        regularisation::Union{Missing,Symbol} = :default, regularisation_params::Dict = Dict(),
                        only_regulise_if_not_PSD::Bool = false,
                        time_grid::Union{Missing,Vector} = missing,
                        fixed_spacing::Union{Missing,<:Real} = missing, refresh_times::Bool = false,
                        rough_guess_number_of_intervals::Integer = 5,
                        kernel::HFC_Kernel{<:Real} = parzen,
                        H::Real = kernel.c_star * ( mean(map(a -> length(ts.groupingrows[a]), assets))   )^0.6,
                        m::Integer = 2, # BNHLS parameters, numJ::Integer = 100, num_blocks::Integer = 10,
                        block_width::Real = (maximum(ts.df[:,ts.time]) - minimum(ts.df[:,ts.time])) / num_blocks,
                        microstructure_noise_var::Dict{Symbol,<:Real} = two_scales_volatility(ts, assets)[2],
                        theta::Real = 0.15, g::NamedTuple = g,
                        equalweight::Bool = false, num_grids::Real = default_num_grids(ts))

This is a convenience wrapper for the regularisation techniques.
### General Inputs
* `ts` - The tick data.
* `assets` - What assets from ts that you want to estimate the covariance for.
* `method`  - The method you want to use. This can be `:simple_covariance`, `:bnhls_covariance`, `:spectral_covariance`, `:preaveraged_covariance` or `:two_scales_covariance`.
* `regularisation` - The regularisation method to use. This can be `:identity_regularisation`, `:eigenvalue_clean`, `:nearest_correlation_matrix` or `:nearest_psd_matrix`. You can also choose `:covariance_default` (which is `:nearest_psd_matrix`) or  `:correlation_default` (which is `:nearest_correlation_matrix`). If missing then the default regularisation method for your chosen covariance estimation method will be used.
* `regularisation_params` - Keyword arguments that will be used by your chosen regularisation method.
* `only_regulise_if_not_PSD` - Should the resultant matrix only be regularised if it is not psd.
#### Inputs only used in `:simple_covariance` method.
* `time_grid` - The grid with which to calculate returns (`:simple_covariance` method only).
* `fixed_spacing` - A spacing used to calculate a time grid. Not used if `refresh_times=true` (`:simple_covariance` method only).
* `refresh_times` - Should refresh times be used to estimate covariance (`:simple_covariance` method only).
* `rough_guess_number_of_intervals` - A rough number of intervals to calculate a default spacing. Not used if a `time_grid` or `fixed_spacing` is provided or if `refresh_times=true` (`:simple_covariance` method only).
#### Inputs only used in `:bnhls_covariance` method.
* `kernel` - The kernel used. See the bnhls paper for details. (`:bnhls_covariance` method only)
* `H` - The number of lags/leads used in estimation. See the bnhls paper for details. (`:bnhls_covariance` method only)
* `m` - The number of end returns to average. (`:bnhls_covariance` method only)
#### Inputs only used in `:spectral_covariance` method.
* `numJ` - The number of J values. See the paper for details (`:spectral_covariance` method only).
* `num_blocks` - The number of blocks to split the time frame into. See the preaveraging paper for details (`:spectral_covariance` method only).
* `block_width` - The width of each block to split the time frame into (`:spectral_covariance` method only).
* `microstructure_noise_var` - Estimates of microstructure noise variance for each asset (`:spectral_covariance` method only).
#### Inputs only used in `:preaveraged_covariance` method.
* `theta` - A theta value. See paper for details (`:preaveraged_covariance` method only).
* `g` - A tuple containing a preaveraging method (with name "f") and a ψ value. See paper for details (`:preaveraged_covariance` method only).
#### Inputs only used in `:two_scales_covariance` method.
* `equalweight` - Should we use equal weight for the two different linear combinations of assets. If false then an optimal weight is calculated (from volatilities) (`:two_scales_covariance` method only).
* `num_grids` - Number of grids used in order in two scales estimation (`:two_scales_covariance` method only).
### Returns
* A `CovarianceMatrix`
"""
function estimate_covariance(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts), method::Symbol = :preaveraged_covariance;
                             regularisation::Union{Missing,Symbol} = :default, regularisation_params::Dict = Dict(),
                             only_regulise_if_not_PSD::Bool = false,
                             time_grid::Union{Missing,Vector} = missing,
                             fixed_spacing::Union{Missing,<:Real} = missing, refresh_times::Bool = false, rough_guess_number_of_intervals::Integer = 5, # General Inputs
                             kernel::HFC_Kernel{<:Real} = parzen, H::Real = kernel.c_star * ( mean(map(a -> length(ts.groupingrows[a]), assets))   )^0.6, m::Integer = 2, # BNHLS parameters
                             numJ::Integer = 100, num_blocks::Integer = 10, block_width::Real = (maximum(ts.df[:,ts.time]) - minimum(ts.df[:,ts.time])) / num_blocks, microstructure_noise_var::Dict{Symbol,<:Real} = two_scales_volatility(ts, assets)[2], # Spectral Covariance parameters
                             theta::Real = 0.15, g::NamedTuple = g, # Preaveraging
                             equalweight::Bool = false, num_grids::Real = default_num_grids(ts)) # Two Scales parameters
    if (ismissing(regularisation) == false) && (regularisation == :default)
        regularisation = (method == :two_scales_covariance) ? :correlation_default : :covariance_default
    end

    if method == :simple_covariance
        return simple_covariance(ts, assets; regularisation = regularisation, only_regulise_if_not_PSD = only_regulise_if_not_PSD,
                                   time_grid = time_grid,
                                   fixed_spacing = fixed_spacing, refresh_times = refresh_times, rough_guess_number_of_intervals = rough_guess_number_of_intervals)
    elseif method == :bnhls_covariance
        return bnhls_covariance(ts, assets; regularisation = regularisation,
                                  only_regulise_if_not_PSD = only_regulise_if_not_PSD, kernel = kernel, H = H,
                                  m = m)
    elseif method == :spectral_covariance
        return spectral_covariance(ts, assets; regularisation = regularisation,
                                     only_regulise_if_not_PSD = only_regulise_if_not_PSD, numJ = numJ, num_blocks = num_blocks, block_width = block_width,
                                     microstructure_noise_var = microstructure_noise_var)
    elseif method == :preaveraged_covariance
        return preaveraged_covariance(ts, assets; regularisation = regularisation,
                                     only_regulise_if_not_PSD = only_regulise_if_not_PSD, theta = theta, g = g)
    elseif method == :two_scales_covariance
        return two_scales_covariance(ts, assets; regularisation = regularisation, only_regulise_if_not_PSD = only_regulise_if_not_PSD,
                                    equalweight = equalweight, num_grids = num_grids)
    else
        error("The covariance method chosen must be :simple_covariance, :bnhls_covariance, :spectral_covariance, :preaveraged_covariance or :two_scales_covariance")
    end
end

# The Hermitian version
"""
    regularise(mat::Hermitian, ts::SortedDataFrame,  mat_labels::Vector, method::Symbol = :correlation_default;
               spacing::Union{Missing,<:Real} = missing,
               weighting_matrix = Diagonal(eltype(mat).(I(size(mat)[1]))),
               doDykstra = true, stop_at_first_correlation_matrix = true, max_iterates = 1000)

This is a convenience wrapper for the regularisation techniques.
### General Inputs
* `mat` - The matrix you want to regularise.
* `ts` - The tick data.
* `mat_labels` - The name of the assets for each row/column of the matrix.
* `method`  - The method you want to use. This can be `:identity_regularisation`, `:eigenvalue_clean`, `:nearest_correlation_matrix` or `:nearest_psd_matrix`. You can also choose `:covariance_default` (which is `:nearest_psd_matrix`) or  `:correlation_default` (which is `:nearest_correlation_matrix`).
#### Inputs only used in `:identity_regularisation` method.
* `spacing` - The interval spacing used in choosing an identity weight (`identity_regularisation` method only).
#### Inputs only used in `:nearest_correlation_matrix` method.
* `weighting_matrix` - The weighting matrix used to calculate the nearest psd matrix (`:nearest_correlation_matrix` method only).
* `doDykstra` - Should a Dykstra correction be applied (`:nearest_correlation_matrix` method only).
* `stop_at_first_correlation_matrix` - Should we stop at first valid correlation matrix (`:nearest_correlation_matrix` method only).
* `max_iterates` - Maximum number of iterates (`:nearest_correlation_matrix` method only).
### Returns
* A `Hermitian`

    regularise(covariance_matrix::CovarianceMatrix, ts::SortedDataFrame, method::Symbol = :nearest_correlation_matrix;
               apply_to_covariance::Bool = true,
               spacing::Union{Missing,<:Real} = missing,
               weighting_matrix = Diagonal(eltype(covariance_matrix.correlation).(I(size(covariance_matrix.correlation)[1]))),
               doDykstra = true, stop_at_first_correlation_matrix = true, max_iterates = 1000)

This is a convenience wrapper for the regularisation techniques.
### General Inputs
* `covariance_matrix` - The matrix you want to regularise.
* `ts` - The tick data.
* `method`  - The method you want to use. This can be `:identity_regularisation`, `:eigenvalue_clean`, `:nearest_correlation_matrix` or `:nearest_psd_matrix`. You can also choose `:covariance_default` (which is `:nearest_psd_matrix`) or  `:correlation_default` (which is `:nearest_correlation_matrix`).
* `apply_to_covariance` - Should regularisation be applied to the covariance matrix. If false it is applied to the correlation matrix.
#### Inputs only used in `:identity_regularisation` method.
* `spacing` - The interval spacing used in choosing an identity weight (`identity_regularisation` method only).
#### Inputs only used in `:nearest_correlation_matrix` method.
* `weighting_matrix` - The weighting matrix used to calculate the nearest psd matrix (`:nearest_correlation_matrix` method only).
* `doDykstra` - Should a Dykstra correction be applied (`:nearest_correlation_matrix` method only).
* `stop_at_first_correlation_matrix` - Should we stop at first valid correlation matrix (`:nearest_correlation_matrix` method only).
* `max_iterates` - Maximum number of iterates (`:nearest_correlation_matrix` method only).
### Returns
* A `CovarianceMatrix`
"""
function regularise(mat::Hermitian, ts::SortedDataFrame,  mat_labels::Vector, method::Symbol = :correlation_default;
                    spacing::Union{Missing,<:Real} = missing,
                    weighting_matrix = Diagonal(eltype(mat).(I(size(mat)[1]))),
                    doDykstra = true, stop_at_first_correlation_matrix = true, max_iterates = 1000)
    if method == :covariance_default
        method = :nearest_psd_matrix
    elseif method == :correlation_default
        method = :nearest_correlation_matrix
    end
    if method == :identity_regularisation
        return identity_regularisation(mat, ts, mat_labels; spacing = spacing)
    elseif method == :eigenvalue_clean
        return eigenvalue_clean(mat, ts)
    elseif method == :nearest_correlation_matrix
        return nearest_correlation_matrix(mat; weighting_matrix = weighting_matrix, doDykstra = doDykstra,
                                          stop_at_first_correlation_matrix = stop_at_first_correlation_matrix, max_iterates = max_iterates)
    elseif method == :nearest_psd_matrix
        return nearest_psd_matrix(mat)
    else
        error("The covariance method chosen must be :identity_regularisation, :eigenvalue_clean, :nearest_correlation_matrix or :nearest_psd_matrix. You can also choose :covariance_default (which is :nearest_psd_matrix) or  :correlation_default (which is :nearest_correlation_matrix).")
    end
end
function regularise(covariance_matrix::CovarianceMatrix, ts::SortedDataFrame, method::Symbol = :nearest_correlation_matrix;
                    apply_to_covariance::Bool = true,
                    spacing::Union{Missing,<:Real} = missing,
                    weighting_matrix = Diagonal(eltype(covariance_matrix.correlation).(I(size(covariance_matrix.correlation)[1]))),
                    doDykstra = true, stop_at_first_correlation_matrix = true, max_iterates = 1000)
    if method == :covariance_default
        method = :nearest_psd_matrix
    elseif method == :correlation_default
        method = :nearest_correlation_matrix
    end
    if method == :identity_regularisation
        return identity_regularisation(covariance_matrix, ts; spacing = spacing, apply_to_covariance = apply_to_covariance)
    elseif method == :eigenvalue_clean
        return eigenvalue_clean(covariance_matrix, ts; apply_to_covariance = apply_to_covariance)
    elseif method == :nearest_correlation_matrix
        return nearest_correlation_matrix(covariance_matrix, ts; weighting_matrix = weighting_matrix, doDykstra = doDykstra,
                                          stop_at_first_correlation_matrix = stop_at_first_correlation_matrix, max_iterates = max_iterates)
    elseif method == :nearest_psd_matrix
        return nearest_psd_matrix(covariance_matrix; apply_to_covariance = apply_to_covariance)
    else
        error("The covariance method chosen must be :identity_regularisation, :eigenvalue_clean, :nearest_correlation_matrix or :nearest_psd_matrix. You can also choose :covariance_default (which is :nearest_psd_matrix) or  :correlation_default (which is :nearest_correlation_matrix).")
    end
end


text:
include("libavdevice_h.jl")

include("avdevice.jl")
#include("version.jl")


text:
module DPFEHM

import AlgebraicMultigrid
import Calculus
import ChainRulesCore
import ForwardDiff
import Interpolations
import NLsolve
import NonlinearEquations
import SparseArrays
import StaticArrays

include("grid.jl")
include("groundwater.jl")
include("RelPerm.jl")
include("richards.jl")

end


text:
# testing discretization procedures
# example is based on Walter Gautschi "Orthogonal Polynomials: Computation and Approximation"
# examples 2.36

using PolyChaos, Test
import LinearAlgebra: norm

nodes = [40, 80, 160, 320]
tol = 1e-14

@time @testset "Stieltjes procedure" begin
    for n in nodes
        myfile = open("fej$n.txt")
        βref = parse.(Float64,readlines(myfile))
        βcom = stieltjes(n,fejer(n)...)[2]
        @test isapprox(norm(βref-βcom,Inf),0.;atol=tol)
        close(myfile)
    end
end

@time @testset "Lanczos procedure" begin
    for n in nodes
        myfile = open("fej$n.txt")
        βref = parse.(Float64,readlines(myfile))
        βcom = lanczos(n,fejer(n)...)[2]
        @test isapprox(norm(βref-βcom,Inf),0.;atol=tol)
        close(myfile)
    end
end


text:
@testset "sparse_variational" begin
    @testset "AbstractGPs interface" begin
        rng = MersenneTwister(123456)
        N_cond = 5
        M = 4
        N_a = 6
        N_b = 7

        # Specify prior.
        f = GP(Matern32Kernel())

        # Sample from prior.
        x = collect(range(-1.0, 1.0; length=N_cond))
        fx = f(x, 1e-3)
        y = rand(rng, fx)

        # Specify inducing variables.
        z = range(-1.0, 1.0; length=M)
        fz = f(z, 1e-6)

        # Construct approximate posterior.
        q_Centered = optimal_variational_posterior(fz, fx, y)
        approx_Centered = SparseVariationalApproximation(Centered(), fz, q_Centered)
        f_approx_post_Centered = posterior(approx_Centered)

        # Check that approximate posterior is self-consistent.
        a = collect(range(-1.0, 1.0; length=N_a))
        b = randn(rng, N_b)
        TestUtils.test_internal_abstractgps_interface(rng, f_approx_post_Centered, a, b)

        @testset "NonCentered" begin

            # Construct optimal approximate posterior.
            q = optimal_variational_posterior(fz, fx, y)
            Cuu = cholesky(Symmetric(cov(fz)))
            q_ε = MvNormal(
                Cuu.L \ (mean(q) - mean(fz)), Symmetric((Cuu.L \ cov(q)) / Cuu.U)
            )

            # Check that q_ε has been properly constructed.
            @test mean(q) ≈ mean(fz) + Cuu.L * mean(q_ε)
            @test cov(q) ≈ Cuu.L * cov(q_ε) * Cuu.U

            # Construct equivalent approximate posteriors.
            approx_non_Centered = SparseVariationalApproximation(NonCentered(), fz, q_ε)
            f_approx_post_non_Centered = posterior(approx_non_Centered)
            TestUtils.test_internal_abstractgps_interface(
                rng, f_approx_post_non_Centered, a, b
            )

            # Unit-test kl_term.
            @test isapprox(
                ApproximateGPs.kl_term(approx_non_Centered, f_approx_post_non_Centered),
                ApproximateGPs.kl_term(approx_Centered, f_approx_post_Centered);
                rtol=1e-5,
            )

            # Verify that the non-centered approximate posterior agrees with centered.
            @test mean(f_approx_post_non_Centered, a) ≈ mean(f_approx_post_Centered, a)
            @test cov(f_approx_post_non_Centered, a, b) ≈ cov(f_approx_post_Centered, a, b)
            @test elbo(approx_non_Centered, fx, y) ≈ elbo(approx_Centered, fx, y)
        end
    end

    @testset "elbo" begin
        rng, N = MersenneTwister(654321), 20
        x = rand(rng, N) * 10
        y = sin.(x) + 0.9 * cos.(x * 1.6) + 0.4 * rand(rng, N)
        z = x[1:5]

        kernel = make_kernel([0.2, 0.6])
        f = GP(kernel)
        fx = f(x, 0.1)
        fz = f(z)
        q_ex = optimal_variational_posterior(fz, fx, y)

        sva = SparseVariationalApproximation(fz, q_ex)
        @test elbo(sva, fx, y) isa Real
        @test elbo(sva, fx, y) ≤ logpdf(fx, y)

        fx_bad = f(x, fill(0.1, N))
        @test_throws ErrorException elbo(sva, fx_bad, y)

        lf = LatentGP(f, GaussianLikelihood(0.1), 1e-18)
        lfx = lf(x)

        @test elbo(sva, lfx, y) ≈ elbo(sva, fx, y) atol = 1e-10
    end

    @testset "equivalences" begin
        rng, N = MersenneTwister(654321), 20
        x = rand(rng, N) * 10
        y = sin.(x) + 0.9 * cos.(x * 1.6) + 0.4 * rand(rng, N)

        z = copy(x) # Set inducing inputs == training inputs

        k_init = [0.2, 0.6] # initial kernel parameters
        lik_noise = 0.1 # The (fixed) Gaussian likelihood noise

        @testset "exact posterior" begin
            # There is a closed form optimal solution for the variational posterior
            # q(u) (e.g. # https://krasserm.github.io/2020/12/12/gaussian-processes-sparse/
            # equations (11) & (12)). The SVGP posterior with this optimal q(u)
            # should therefore be equivalent to the sparse GP (Titsias) posterior
            # and exact GP regression (when z == x).

            kernel = make_kernel(k_init)
            f = GP(kernel)
            fx = f(x, lik_noise)
            fz = f(z)
            q_ex = optimal_variational_posterior(fz, fx, y)

            gpr_post = posterior(fx, y) # Exact GP regression
            vfe_post = posterior(VFE(fz), fx, y) # Titsias posterior
            svgp_post = posterior(SparseVariationalApproximation(Centered(), fz, q_ex)) # Hensman (2013) exact posterior

            @test mean(gpr_post, x) ≈ mean(svgp_post, x) atol = 1e-10
            @test cov(gpr_post, x) ≈ cov(svgp_post, x) atol = 1e-10

            @test mean(vfe_post, x) ≈ mean(svgp_post, x) atol = 1e-10
            @test cov(vfe_post, x) ≈ cov(svgp_post, x) atol = 1e-10

            @test(
                isapprox(
                    elbo(SparseVariationalApproximation(Centered(), fz, q_ex), fx, y),
                    logpdf(fx, y);
                    atol=1e-6,
                )
            )
        end

        @testset "optimised posterior" begin
            jitter = 1e-5

            make_gp(kernel) = GP(kernel)

            # SVGP model
            struct SVGPModel
                k # kernel parameters
                z # inducing points
                m # variational mean
                A # variational covariance sqrt (Σ = A'A)
            end
            Flux.@functor SVGPModel (m, A) # Only train the variational parameters

            function construct_parts(m::SVGPModel, x)
                f = make_gp(make_kernel(m.k))
                fx = f(x, lik_noise)
                fz = f(m.z, jitter)

                S = PDMat(Cholesky(LowerTriangular(m.A)))
                q = MvNormal(m.m, S)
                return SparseVariationalApproximation(fz, q), fx
            end

            # initialise the variational parameters
            m, A = zeros(N), Matrix{Float64}(I, N, N)
            svgp_model = SVGPModel(copy(k_init), copy(z), m, A)
            function svgp_loss(x, y)
                approx, fx = construct_parts(svgp_model, x)
                return -elbo(approx, fx, y)
            end

            # Train the SVGP model
            data = [(x, y)]
            opt = ADAM(0.001)

            svgp_ps = Flux.params(svgp_model)

            # Optimise q(u)
            Flux.train!((x, y) -> svgp_loss(x, y), svgp_ps, ncycle(data, 20000), opt)

            ## construct the posteriors
            f_gpr = make_gp(make_kernel(k_init))
            gpr_post = posterior(f_gpr(x, lik_noise), y)

            svgp_post = posterior(first(construct_parts(svgp_model, x)))

            ## FIFTH - test equivalences
            @test all(isapprox.(mean(gpr_post, x), mean(svgp_post, x), atol=1e-4))
            @test all(isapprox.(cov(gpr_post, x), cov(svgp_post, x), atol=1e-4))
        end
    end
end


text:
Array{T}(::UndefInitializer, ::Infinity) where T = throw(ArgumentError("Cannot create infinite Array"))

Array{T}(::UndefInitializer, ::Infinity, ::Infinity) where T = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::Infinity, ::Integer) where T = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::Integer, ::Infinity) where T = throw(ArgumentError("Cannot create infinite Array"))

Matrix{T}(::UndefInitializer, ::Infinity, ::Infinity) where T = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::Infinity, ::Integer) where T = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::Integer, ::Infinity) where T = throw(ArgumentError("Cannot create infinite Array"))

Vector{T}(::UndefInitializer, ::Infinity) where T = throw(ArgumentError("Cannot create infinite Array"))

# This gets called when infinit number of columns
print_matrix_row(io::IO,
        X::AbstractVecOrMat, A::Vector,
        i::Integer, cols::AbstractVector{<:Infinity}, sep::AbstractString) = nothing

print_matrix_vdots(io::IO, vdots::AbstractString,
        A::Vector, sep::AbstractString, M::Integer, ::NotANumber) = nothing

# Avoid infinite loops on maximum
Base.mapreduce_impl(f, op, A::AbstractArray, ifirst::Integer, ::Infinity) =
    throw(ArgumentError("Cannot call mapreduce on an infinite length $(typeof(A))"))

function show_delim_array(io::IO, itr::AbstractArray, op, delim, cl,
                          delim_one, i1, ::Infinity)
    print(io, op)
    l = 20
    if !show_circular(io, itr)
        recur_io = IOContext(io, :SHOWN_SET => itr)
        if !haskey(io, :compact)
          recur_io = IOContext(recur_io, :compact => true)
        end
        first = true
        i = i1
        if 20 >= i1
          while true
              if !isassigned(itr, i)
                  print(io, undef_ref_str)
              else
                  x = itr[i]
                  show(recur_io, x)
              end
              i += 1
              if i > l
                  print(io, delim)
                  print(io, ' ')
                  print(io, '…')
                  delim_one && first && print(io, delim)
                  break
              end
              first = false
              print(io, delim)
              print(io, ' ')
          end
        end
    end
    print(io, cl)
end

#####
# FillArrays
#####

# Lazy Broadacasting
for typ in (:Ones, :Zeros, :Fill)
    @eval begin
        BroadcastStyle(::Type{$typ{T,N,NTuple{N,<:OneToInf}}}) where {T,N} = LazyArrayStyle{N}()
        BroadcastStyle(::Type{$typ{T,2,<:Tuple{<:Any,<:OneToInf}}}) where {T} = LazyArrayStyle{2}()
        BroadcastStyle(::Type{$typ{T,2,<:Tuple{<:OneToInf,<:Any}}}) where {T} = LazyArrayStyle{2}()
    end
end

BroadcastStyle(::Type{Eye{T,OneToInf{I}}}) where {T,I} = LazyArrayStyle{2}()

#####
# Diagonal
#####

BroadcastStyle(::Type{<:Diagonal{<:Any,<:AbstractInfUnitRange}}) = LazyArrayStyle{2}()


text:
#=
run stepper benchmarks
see stepper/README.md
=#

using Test
using Printf
import DataFrames
import CSV
import Hypatia
import Hypatia.Solvers

include(joinpath(@__DIR__, "../../examples/Examples.jl"))
using Main.Examples

# path to write results DataFrame to CSV
results_path = joinpath(mkpath(joinpath(@__DIR__, "raw")), "bench.csv")

# script verbosity
script_verbose = false

# default options to solvers
default_options = (
    verbose = false,
    # verbose = true,
    default_tol_relax = 10,
    iter_limit = 10000,
    time_limit = 10000,
    )

# stepper option sets to run
porc = Solvers.PredOrCentStepper{Float64}
comb = Solvers.CombinedStepper{Float64}
stepper_options = [
    "basic" => porc(use_adjustment = false, use_curve_search = false,
        use_max_prox = false, prox_bound = 0.2844),
    "prox" => porc(use_adjustment = false, use_curve_search = false),
    "toa" => porc(use_adjustment = true, use_curve_search = false),
    "curve" => porc(use_adjustment = true, use_curve_search = true),
    "comb" => comb(shift_sched = 0),
    "prox_val2" => porc(use_adjustment = false, use_curve_search = false,
        prox_bound = 0.2844),
    "prox_val3" => porc(use_adjustment = false, use_curve_search = false,
        prox_bound = 0.5),
    "prox_val4" => porc(use_adjustment = false, use_curve_search = false,
        prox_bound = 0.9),
    "prox_val5" => porc(use_adjustment = false, use_curve_search = false,
        prox_bound = 0.9999),
    "comb_val2" => comb(shift_sched = 0, prox_bound = 0.2844),
    "comb_val3" => comb(shift_sched = 0, prox_bound = 0.5),
    "comb_val4" => comb(shift_sched = 0, prox_bound = 0.9),
    "comb_val5" => comb(shift_sched = 0, prox_bound = 0.9999),
    ]

# instance sets and real types to run and corresponding time limits (seconds)
inst_sets = [
    # "minimal",
    # "fast",
    "compile",
    "various",
    ]

time_all = time()

@testset "examples tests" begin
test_insts = Examples.get_test_instances()
perf = Examples.setup_benchmark_dataframe()
CSV.write(results_path, perf)

@testset "$mod, $ex" for (mod, mod_insts) in test_insts,
    (ex, (ex_type, ex_insts)) in mod_insts
@testset "$inst_set" for inst_set in inst_sets
    if inst_set == "compile"
        haskey(ex_insts, "various") || continue
        ex_insts["compile"] = ex_insts["various"]
    end
    haskey(ex_insts, inst_set) || continue
    inst_subset = ex_insts[inst_set]
    isempty(inst_subset) && continue

    for (step_name, stepper) in stepper_options
        info_perf = (; inst_set, :example => ex, :model_type => mod,
            :real_T => Float64, :solver_options => (step_name,))
        new_default_options = (; default_options..., stepper = stepper)

        str = "$mod $ex $inst_set $step_name"
        println("\nstarting $str")
        @testset "$str" begin
            Examples.run_instance_set(inst_subset, ex_type{Float64}, info_perf,
                new_default_options, script_verbose, perf, results_path)
        end
    end

end
end

println("\n")
DataFrames.show(perf, allrows = true, allcols = true)
println("\n")
flush(stdout); flush(stderr)
end

@printf("\nbenchmarks total time: %8.2e seconds\n\n", time() - time_all)
;


text:
# ---
# title: 1592. Rearrange Spaces Between Words
# id: problem1592
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/rearrange-spaces-between-words/description/>
# hidden: true
# ---
#
# You are given a string `text` of words that are placed among some number of
# spaces. Each word consists of one or more lowercase English letters and are
# separated by at least one space. It's guaranteed that `text` **contains at
# least one word**.
#
# Rearrange the spaces so that there is an **equal** number of spaces between
# every pair of adjacent words and that number is **maximized**. If you cannot
# redistribute all the spaces equally, place the **extra spaces at the end** ,
# meaning the returned string should be the same length as `text`.
#
# Return _the string after rearranging the spaces_.
#
#
#
# **Example 1:**
#
#
#
#     Input: text = "  this   is  a sentence "
#     Output: "this   is   a   sentence"
#     Explanation: There are a total of 9 spaces and 4 words. We can evenly divide the 9 spaces between the words: 9 / (4-1) = 3 spaces.
#
#
# **Example 2:**
#
#
#
#     Input: text = " practice   makes   perfect"
#     Output: "practice   makes   perfect "
#     Explanation:  There are a total of 7 spaces and 3 words. 7 / (3-1) = 3 spaces plus 1 extra space. We place this extra space at the end of the string.
#
#
# **Example 3:**
#
#
#
#     Input: text = "hello   world"
#     Output: "hello   world"
#
#
# **Example 4:**
#
#
#
#     Input: text = "  walks  udp package   into  bar a"
#     Output: "walks  udp  package  into  bar  a "
#
#
# **Example 5:**
#
#
#
#     Input: text = "a"
#     Output: "a"
#
#
#
#
# **Constraints:**
#
#   * `1 <= text.length <= 100`
#   * `text` consists of lowercase English letters and `' '`.
#   * `text` contains at least one word.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# TODO disabled

function blockdiag(::Type{<:ValueGraph{V, E_VAL}}, iter::AbstractGraph{<:Integer}...) where {V, E_VAL}
    n::V = V(0)
    for g in iter
        n += nv(g)
        # TODO check for overflow
    end

    resultg = ValueGraph(n, E_VAL)

    # TODO this is not very efficient
    Δ::V = zero(V)
    for g in iter
        w = weights(g)
        for u in vertices(g)
            for v in neighbors(g, u)
                add_edge!(resultg, V(u) + Δ, V(v) + Δ, convert(E_VAL, w[u, v]))
            end
        end
        Δ += nv(g)
    end

    return resultg
end

function blockdiag(::Type{<:ValueDiGraph{V, E_VAL}}, iter::AbstractGraph{<:Integer}...) where {V, E_VAL}
    n::V = V(0)
    for g in iter
        n += nv(g)
        # TODO check for overflow
    end

    resultg = ValueDiGraph(n, E_VAL)

    # TODO this is not very efficient
    Δ::V = zero(V)
    for g in iter
        w = weights(g)
        for u in vertices(g)
            for v in outneighbors(g, u)
                add_edge!(resultg, V(u) + Δ, V(v) + Δ, convert(E_VAL, w[u, v]))
            end
        end
        Δ += nv(g)
    end

    return resultg
end

blockdiag(g::AbstractValueGraph, iter::AbstractGraph...) = blockdiag(typeof(g), g, iter...)

# ==== weight related operators =====

# TODO docstring
# TODO maybe better pass a ValueEdge to f
function transform_edgevals!(f::Function, g::ValueGraph)
    V = eltype(g)
    E_VAL = edgeval_type(g)
    n = nv(g)
    fadjlist = g.fadjlist
    edgevals = g.edgevals
    for i = 1:n
        list_i = fadjlist[i]
        len = length(list_i)
        index::Int = 1
        while index <= len && list_i[index] <= i
            j = list_i[index]
            new_val = f(V(i), V(j), value_for_index(edgevals, E_VAL, i, index))
            set_value_for_index!(edgevals, V(i), index, new_val)

            index2 = searchsortedfirst(fadjlist[j], i)
            set_value_for_index!(edgevals, V(j), index2, new_val)

            index += 1
        end
    end
end

function transform_edgevals!(f::Function, g::ValueGraph{V, E_VAL, <: TupleOrNamedTuple}, key) where {V, E_VAL}
    n = nv(g)
    fadjlist = g.fadjlist
    edgevals = g.edgevals
    for i = 1:n
        list_i = fadjlist[i]
        len = length(list_i)
        index::Int = 1
        while index <= len && list_i[index] <= i
            j = list_i[index]
            new_val = f(V(i), V(j), value_for_index(edgevals, E_VAL, i, index, key))
            set_value_for_index!(edgevals, i, index, key, new_val)

            index2 = searchsortedfirst(fadjlist[j], i)
            set_value_for_index!(edgevals, j, index2, key, new_val)

            index += 1
        end
    end
end

function transform_edgevals!(f::Function, g::OutValueDiGraph)
    V = eltype(g)
    E_VAL = edgeval_type(g)
    n = nv(g)
    fadjlist = g.fadjlist
    edgevals = g.edgevals
    for i = 1:n
        list_i = fadjlist[i]
        len = length(list_i)
        index::Int = 1
        while index <= len
            j = list_i[index]
            new_val = f(V(i), V(j), value_for_index(edgevals, E_VAL, i, index))
            set_value_for_index!(edgevals, V(i), index, new_val)

            index += 1
        end
    end
end

function transform_edgevals!(f::Function, g::OutValueDiGraph{V, E_VAL, <: TupleOrNamedTuple}, key) where {V, E_VAL}
    n = nv(g)
    fadjlist = g.fadjlist
    edgevals = g.edgevals
    for i = 1:n
        list_i = fadjlist[i]
        len = length(list_i)
        index::Int = 1
        while index <= len
            j = list_i[index]
            new_val = f(V(i), V(j), value_for_index(edgevals, E_VAL, i, index, key))
            set_value_for_index!(edgevals, i, index, key, new_val)

            index += 1
        end
    end
end


text:
abstract type PFSoftware end

"Organon software corrections"
mutable struct Organon <: PFSoftware end

"Anarede software corrections"
mutable struct ANAREDE <: PFSoftware end

function generators_from_bus(pm_data::Dict, bus::Int; filters::Vector = [])
    filters = vcat(gen -> gen["gen_bus"] == bus, filters)
    return findall(
            x -> (
                all([f(x) for f in filters])
            ),
            pm_data["gen"]
        )
end

function load_from_bus(pm_data::Dict, bus::Int; filters::Vector = [])
    filters = vcat(load -> load["load_bus"] == bus, filters)
    return findall(
            x -> (
                all([f(x) for f in filters])
            ),
            pm_data["load"]
        )
end

function sum_generators_power_and_turn_off(pm_data::Dict, gen_keys::Vector)
    Pg = 0.0
    Qg = 0.0
    for (i, key) in enumerate(gen_keys)
        gen = pm_data["gen"][key]
        Pg += gen["pg"]
        Qg += gen["qg"]
        gen["gen_status"] = 0
    end
    return Pg, Qg
end

function _fix_shunt_voltage_bounds(shunt::Dict, pm_data::Dict)
    shunt["control_data"]["vmmin"] = pm_data["bus"]["$(shunt["shunt_bus"])"]["vm"]
    shunt["control_data"]["vmmax"] = pm_data["bus"]["$(shunt["shunt_bus"])"]["vm"]
end

function _pwf2pm_corrections!(pm_data::Dict, pwf_data::Dict, software::PFSoftware)
    _pwf2pm_corrections_shunt!(pm_data, software)
    _pwf2pm_corrections_gen!(pm_data, pwf_data, software)
    _pwf2pm_corrections_PV!(pm_data, pwf_data, software)
    _pwf2pm_corrections_PQ!(pm_data, software)

    return
end

function _set_controlled_bus_voltage_bounds!(pm_data::Dict, shunt_control_data::Dict)
    ctrl_bus = shunt_control_data["controlled_bus"]
    pm_data["bus"]["$ctrl_bus"]["control_data"]["vmmin"] = pop!(shunt_control_data, "vmmin")
    pm_data["bus"]["$ctrl_bus"]["control_data"]["vmmax"] = pop!(shunt_control_data, "vmmax")
end

function _correct_pwf_network_data(pm_data::Dict)
    mva_base = pm_data["baseMVA"]

    rescale        = x -> x/mva_base

    if haskey(pm_data, "shunt")
        for (i, shunt) in pm_data["shunt"]
            if haskey(shunt, "control_data")
                PowerModels._apply_func!(shunt["control_data"], "bsmin", rescale)
                PowerModels._apply_func!(shunt["control_data"], "bsmax", rescale)

                _set_controlled_bus_voltage_bounds!(pm_data, shunt["control_data"])
            end
        end
    end

end


text:
@testset "Dispatcher" begin
    dispatcher = Vulkan.APIDispatcher()
    @test !haskey(dispatcher.instance_table, :vkEnumerateDeviceExtensionProperties)
    Vulkan.add_fptr!(dispatcher, instance, :vkEnumerateDeviceExtensionProperties)
    @test haskey(dispatcher.instance_table, :vkEnumerateDeviceExtensionProperties)
    @test !isnothing(enumerate_instance_extension_properties(Vulkan.get_fptr(dispatcher, nothing, :vkEnumerateInstanceExtensionProperties)))
    @test !isnothing(enumerate_device_extension_properties(device.physical_device, Vulkan.get_fptr(dispatcher, instance, :vkEnumerateDeviceExtensionProperties)))
    @test Fence(device, Vulkan.get_fptr(dispatcher, device, :vkCreateFence), Vulkan.get_fptr(dispatcher, device, :vkDestroyFence)) isa Fence
    if WITH_DEBUG
        @test Vulkan.add_fptr!(dispatcher, instance, :vkCreateDebugUtilsMessengerEXT) ≠ C_NULL
    end
    @test_throws ErrorException Vulkan.add_fptr!(dispatcher, device, :vkCmdPushDescriptorSetWithTemplateKHR)
end

GC.gc()


text:
"""
    increment_bins!(h, data)

Increments the bin counts for a 1D histogram `h` using the `data`. Elements of `data` that are outside the range of
the histogram's bins will NOT be filtered out, they will be considered as members of the closest bin.
"""
increment_bins!(h, data) = increment_bins!(BinSearchAlgorithm(h), HistogramParallelization(h), h, data)

"""
    increment_bins!(h, data1, data2)

Increments the bin counts for a 2D histogram `h` using the data `data1` and `data2`. Elements of `data` that are
outside the range of the histogram's bins will NOT be filtered out, they will be considered as members of the closest
bin.
"""
increment_bins!(h, data1, data2) = increment_bins!(BinSearchAlgorithm(h), HistogramParallelization(h), h, data1, data2)

"""
    get_weights(h)

Returns the weights array. All histograms must implement this.
"""
get_weights

"""
    get_subweights(h)

Returns the subweights array. All histograms implementing SIMD and PrivateThreads parallelization must implement this.
"""
get_subweights

"""
    counts(h)

Returns the bin counts of the histogram `h`. All histograms must implement this.
"""
counts(h) = get_weights(h)

"""
    zero!(h)

Sets all bin counts of the histogram `h` to zero. All histograms must implement this.
"""
function zero!(h)
    get_weights(h) .= 0
    get_subweights(h) .= 0
    nothing
end

function increment_bins!(::BinSearchAlgorithm, ::NoParallelization, h, data)
    for d in data
        @inbounds i = bin_search(h, 1, d)
        @inbounds get_weights(h)[i] += 1
    end

    nothing
end

function increment_bins!(
    ::BinSearchAlgorithm,
    ::NoParallelization,
    h,
    data1::Union{AbstractVector,AbstractMatrix},
    data2::Union{AbstractVector,AbstractMatrix},
)
    for c = 1:size(data1, 2)
        for r = 1:size(data1, 1)
            @inbounds tx = data1[r, c]
            @inbounds ty = data2[r, c]
            @inbounds ix = bin_search(h, 1, tx)
            @inbounds iy = bin_search(h, 2, ty)
            @inbounds get_weights(h)[ix, iy] += 1
        end
    end

    nothing
end

# This implementation is slower than NoParallelization in my benchmarks.
# TODO: Match performance with https://github.com/pcjung/AHTL
function increment_bins!(::Arithmetic, ::SIMD, h, data::Union{AbstractVector,AbstractMatrix})
    rows = size(data, 1)
    align_rows = floor(Int, rows / 3)

    w = get_subweights(h)
    for c = 1:size(data, 2)
        r = 1

        while r < align_rows
            @turbo for i = 0:2
                @inbounds tx = data[r+i, c]
                @inbounds w[bin_search(h, 1, tx), i+1] += 1
            end
            r += 3
        end

        for r2 = r:rows
            @inbounds tx = data[r2, c]
            @inbounds w[bin_search(h, 1, tx), 1] += 1
        end
    end

    sum_subweights!(h)

    nothing
end

# This implementation is slower than NoParallelization in my benchmarks.
# TODO: Match performance with https://github.com/pcjung/AHTL
function increment_bins!(
    ::Arithmetic,
    ::SIMD,
    h,
    img1::Union{AbstractVector,AbstractMatrix},
    img2::Union{AbstractVector,AbstractMatrix},
)
    rows = size(img1, 1)
    align_rows = floor(Int, rows / 4)

    w = get_subweights(h)
    for c = 1:size(img1, 2)
        r = 1

        while r < align_rows
            @turbo for i = 0:3
                @inbounds tx = img1[r+i, c]
                @inbounds ty = img2[r+i, c]
                ix = bin_search(h, 1, tx)
                iy = bin_search(h, 2, ty)
                @inbounds w[ix, iy, i+1] += 1
            end
            r += 4
        end

        for r2 = r:rows
            @inbounds tx = img1[r2, c]
            @inbounds ty = img2[r2, c]
            ix = bin_search(h, 1, tx)
            iy = bin_search(h, 2, ty)
            @inbounds w[ix, iy, 1] += 1
        end
    end

    sum_subweights!(h)

    nothing
end

function increment_bins!(::BinSearchAlgorithm, ::PrivateThreads, h, data::Union{AbstractVector,AbstractMatrix})
    nthreads = size(get_subweights(h), 2)
    Threads.@threads for thread_idx = 1:nthreads
        if length(data) <= nthreads
            if thread_idx == 1
                start_idx = 1
                end_idx = length(data)
            else
                # Make an empty range
                start_idx = 1
                end_idx = 0
            end
        else
            step = length(data) ÷ nthreads
            start_idx = (thread_idx - 1) * step + 1
            end_idx = step * thread_idx
            if end_idx + step > length(data) - (length(data) % nthreads)
                end_idx = length(data)
            end
        end

        for i = start_idx:end_idx
            @inbounds i = bin_search(h, 1, data[i])
            @inbounds get_subweights(h)[i, thread_idx] += 1
        end
    end

    sum!(h.weights, get_subweights(h))

    nothing
end

function increment_bins!(
    ::BinSearchAlgorithm,
    ::PrivateThreads,
    h,
    data1::Union{AbstractVector,AbstractMatrix},
    data2::Union{AbstractVector,AbstractMatrix},
)
    nthreads = size(get_subweights(h), 3)
    Threads.@threads for thread_idx = 1:nthreads
        if length(data1) <= nthreads
            if thread_idx == 1
                start_idx = 1
                end_idx = length(data1)
            else
                # Make an empty range
                start_idx = 1
                end_idx = 0
            end
        else
            step = length(data1) ÷ nthreads
            start_idx = (thread_idx - 1) * step + 1
            end_idx = step * thread_idx
            if end_idx + step > length(data1) - (length(data1) % nthreads)
                end_idx = length(data1)
            end
        end

        for i = start_idx:end_idx
            @inbounds ix = bin_search(h, 1, data1[i])
            @inbounds iy = bin_search(h, 2, data2[i])
            @inbounds get_subweights(h)[ix, iy, thread_idx] += 1
        end
    end

    sum!(h.weights, get_subweights(h))

    nothing
end

function increment_bins!(::HashFunction, ::NoParallelization, h, data::AbstractVector{<:AbstractString})
    t = get_table(h)
    for d in data
        t[d] = get!(t, d, 0) + 1
    end

    nothing
end

function increment_bins!(::HashFunction, ::PrivateThreads, h, data::AbstractVector{<:AbstractString})
    nthreads = length(get_subtable(h))
    Threads.@threads for thread_idx = 1:nthreads
        if length(data) <= nthreads
            if thread_idx == 1
                start_idx = 1
                end_idx = length(data)
            else
                # Make an empty range
                start_idx = 1
                end_idx = 0
            end
        else
            step = length(data) ÷ nthreads
            start_idx = (thread_idx - 1) * step + 1
            end_idx = step * thread_idx
            if end_idx + step > length(data) - (length(data) % nthreads)
                end_idx = length(data)
            end
        end

        st = get_subtable(h)[thread_idx]
        for i = start_idx:end_idx
            @inbounds d = data[i]
            st[d] = get!(st, d, 0) + 1
        end
    end

    t = get_table(h)
    for thread_idx = 1:nthreads
        st = get_subtable(h)[thread_idx]
        for (k, v) in pairs(st)
            t[k] = get!(t, k, 0) + v
        end
    end

    nothing
end

sum_subweights!(h) = sum!(get_weights(h), get_subweights(h))


text:
module KerrQuasinormalModes

using LinearAlgebra
using NLsolve
using Statistics
using CSV
using StaticArrays
using PyCall

include("AngularFunctions.jl")
include("ParameterFunctions.jl")
include("ContinuedFractionFunctions.jl")
include("SolverFunctions.jl")
include("Schwarzschild.jl")
include("ModeCalculator.jl")
include("QnmRotationSeries.jl")
include("SpinWeightedSphericalLookup.jl")
include("Interface.jl")
include("ConfluentHeun.jl")
include("LinearCombinations.jl")
include("Derivative.jl")

#Ψ = qnmfunction(s=-2,l=2,n=2,m=-2,a=0.4)
#Ψ(2.0)
function asad()
    print("blah blah blo new yesterday heres somthing new")
end

export SpinWeightedSpherical
export HeunConfluentRadial
export qnmfunction
export QuasinormalModeFunction
export asad
export ComplexPlot

#For debugging
export HeunConfluentLocalExpansion, ComputeSeries, qnm, Custom
#export importqnm, qnmfunctionnew
end


text:
using Test
using Quantica

@testset "Quantica.jl" begin
    include("test_lattice.jl")
    include("test_model.jl")
    include("test_ket.jl")
    include("test_hamiltonian.jl")
    include("test_mesh.jl")
    include("test_bandstructure.jl")
    include("test_greens.jl")
    include("test_KPM.jl")
end

max_stars_repo_name: pablosanjose/Quantica

max_issues_repo_name: pablosanjose/Quantica

max_forks_repo_name: pablosanjose/Quantica


text:
using Knet
libknet8handle = Libdl.dlopen(Knet.libknet8)

SIZE = 100000
ITER = 100000
x32 = KnetArray(rand(Float32,SIZE))
y32 = similar(x32)
x64 = KnetArray(rand(Float64,SIZE))
y64 = similar(x64)

function cuda1test(fname, jname=fname, o...)
    println(fname)
    fcpu = eval(parse(jname))
    f32 = Libdl.dlsym(libknet8handle, fname*"_32")
    @time cuda1rep(f32,x32,y32)
    isapprox(Array(y32),fcpu(Array(x32))) || warn("$fname 32")
    f64 = Libdl.dlsym(libknet8handle, fname*"_64")
    @time cuda1rep(f64,x64,y64)
    isapprox(Array(y64),fcpu(Array(x64))) || warn("$fname 64")
end

function cuda1rep{T}(f,x::KnetArray{T},y::KnetArray{T})
    n = Cint(length(y))
    for i=1:ITER
        ccall(f,Void,(Cint,Ptr{T},Ptr{T}),n,x,y)
    end
    Knet.@cuda(cudart,cudaDeviceSynchronize,())
    Knet.@cuda(cudart,cudaGetLastError,())
end

for f in Knet.cuda1
    isa(f,Tuple) || (f=(f,))
    cuda1test(f...)
end


text:
include("../Tests/Tests.jl")
include("utilities.jl")

using Test, JuMP

@testset "ZeroPolynomial" begin
    include("zero_polynomial.jl")
end
@testset "ZeroPolynomialInFixedVariablesSet" begin
    include("zero_polynomial_in_fixed_variables_set.jl")
end
@testset "ZeroPolynomialInAlgebraicSet" begin
    include("zero_polynomial_in_algebraic_set.jl")
end
@testset "PlusMinus" begin
    include("plus_minus.jl")
end


text:
using VariableTypes
using Base.Test

@testset "Single value dispatch" begin
    @test vartype(1) === Discrete
    @test vartype(1.) === Continuous
    @test vartype(true) === Categorical
    @test vartype("a") === Categorical
end

@testset "Vector dispatch" begin
    @test vartype(repeat([1], inner=10)) === Discrete
    @test vartype(repeat([2.0], inner=15)) === Continuous
    a = Array{Any, 1}(); append!(a, 1); append!(a, 2); append!(a, 3.5)
    @test vartype(a, 2) === Discrete
    @test vartype([1, 2, 3, "a"]) === Categorical
    @test vartype([1, 2, 3.]) === Continuous
end

@testset "Matrix dispatch" begin
    a = Matrix{Any}(10, 3)
    a[:, 1] = rand(10)
    a[:, 2] = repeat([1, 2], inner=5)
    a[:, 3] = repeat(["A", "B"], inner=5)

    @test all(vartype(a, 1, 10) .=== [Continuous, Discrete, Categorical])
    @test all(vartype(a, 2) .=== repeat([Categorical], inner=10))
end


text:
using CompScienceMeshes
using Test
using SauterSchwabQuadrature

pI = point(1,5,3)
pII = point(2,5,3)
pIII = point(7,1,0)
pIV = point(5,1,-3)
pV = point(0,0,0)

Sourcechart = simplex(pI,pIII,pII)
Testchart = simplex(pI,pIV,pV)

Accuracy = 12
cv = CommonVertex(SauterSchwabQuadrature._legendre(Accuracy,0.0,1.0))

function integrand(x,y)
			return(((x-pI)'*(y-pV))*exp(-im*1*norm(x-y))/(4pi*norm(x-y)))
 end

 function INTEGRAND(û,v̂)
 	n1 = neighborhood(Testchart, û)
 	n2 = neighborhood(Sourcechart, v̂)
 	x = cartesian(n1)
 	y = cartesian(n2)
 	output = integrand(x,y)*jacobian(n1)*jacobian(n2)
  	return(output)
 end

 result = sauterschwab_parameterized(INTEGRAND, cv)-
 			verifintegral2(Sourcechart, Testchart, integrand, Accuracy)

@test norm(result) < 1.e-3

kernel(x,y) = 1/norm(cartesian(x)-cartesian(y))

t1 = simplex(
    @SVector[0.180878, -0.941848, -0.283207],
    @SVector[0.0, -0.980785, -0.19509],
    @SVector[0.0, -0.92388, -0.382683])
t2 = simplex(
    @SVector[0.180878, -0.941848, -0.283207],
    @SVector[0.373086, -0.881524, -0.289348],
    @SVector[0.294908, -0.944921, -0.141962])

@test indexin(t1.vertices, t2.vertices) == [1, nothing, nothing]

rt = RTRefSpace{Float64}()
igd = generate_integrand_uv(kernel, rt, rt, t1, t2)

i5 = sauterschwab_parameterized(igd, CommonVertex(SauterSchwabQuadrature._legendre(5,0.0,1.0)))
i10 = sauterschwab_parameterized(igd, CommonVertex(SauterSchwabQuadrature._legendre(10,0.0,1.0)))

# brute numerical approach
q1 = quadpoints(t1, 10)
q2 = quadpoints(t2, 10)

M = N = numfunctions(rt)
iref = zero(i5)
for (x,w1) in q1
    f = rt(x)
    for (y,w2) in q2
        g = rt(y)
        G = kernel(x,y)
        ds = w1*w2
        global iref += SMatrix{M,N}([dot(f[i][1], G*g[j][1])*ds for i=1:M, j=1:N])
    end
end

include(joinpath(dirname(@__FILE__,),"numquad.jl"))
ibf = numquad(kernel, rt, rt, t1, t2, zero(i5))

@test i5  ≈ iref atol=1e-7
@test i10 ≈ iref atol=1e-7
@test i10 ≈ ibf  atol=1e-8


text:
using JuliaDB, IndexedTables, OnlineStats

using StatsPlots, StaticArrays   # Note the name change from Stat[s]Plots

cd("/Users/malcolm/PacktPub/Chp09")

# Load the data
flights = loadtable("Files/hhflights.csv");

# Define a filter (1st day os each month)

filter(i -> (i.Month == 1) && (i.DayofMonth == 1), flights)

# Select all flights Departure/Arrival/FlightNum
select(flights, (:DepTime, :ArrTime, :FlightNum))

filter(i -> i.DepDelay > 60, select(flights, (:UniqueCarrier, :DepDelay)))

# Apply function to each row, using 'map'

speedy = map(i -> i.Distance / i.AirTime * 60, flights)

# Add speed to the dataset
flights = pushcol(flights, :Speedy, speedy)

using ShiftedArrays
y = groupby(length, flights, :Month)
lengths = columns(y, :length)
pushcol(y, :change, lengths .- lag(lengths))

# Apply several operations
import Lazy
Lazy.@as x flights begin
    select(x, (:UniqueCarrier, :DepDelay))
    filter(i -> i.DepDelay > 60, x)
end

# Performance

using BenchmarkTools
sortedflights = reindex(flights, :Dest)

# println("Presorted timing:")
# @benchmark groupby(NamedTuple(avg_delay = mean.(dropmissing(sortedflights), sortedflights, select = :ArrDelay)))

# Grouping

using StatsBase, Nullables
fc = filter(t->!isnull(t.DepDelay), flights)

gfc = groupby(fc, :UniqueCarrier, select = (:Month, :DayofMonth, :DepDelay), flatten = true) do dd
    rks = ordinalrank(column(dd, :DepDelay), rev = true)
    sort(dd[rks .<= 2], by =  i -> i.DepDelay, rev = true)
end

#=
groupby(fc, :UniqueCarrier, select = (:Month, :DayofMonth, :DepDelay), flatten = true) do dd
    select(dd, 1:2, by = i -> i.DepDelay, rev = true)
end
=#

# Visualization

using StatsPlots
gr(fmt = :png) # choose the fast GR backend and set format to png: svg would probably crash with so many points
@df flights scatter(:DepDelay, :ArrDelay, group = :Distance .> 1000,  fmt = :png, layout = 2, legend = :topleft)

# using Pkg

# path = Pkg.dir("JuliaDB", "test", "sample")
# path = joinpath(dirname(pathof(JuliaDB)), "..","test","sample")
path =  joinpath(homedir(), "PacktPub","Chp09","Files","Stocks")

### sampledata = loadfiles(path, indexcols=["date", "ticker"])
stockdata = loadndsparse(path, indexcols=["date", "ticker"])

using Dates
stockdata[Date("2010-06-01"), "GOOGL"]

stockdata[Date("2012-01"):Dates.Month(1):Date("2014-12"), ["GOOGL", "KO"]]

filter(x ->  x.close >= 100.0 && x.close <= 140.0, stockdata[:, "GS"])

### select(stockdata, 1=>Dates.ismonday, 2=>x->startswith(x, "G"))
goFri = filter((1=>Dates.isfriday, 2=>x->startswith(x, "GO")), stockdata)

# Select just the GOOGLE stocks
googl = stockdata[:, ["GOOGL"]]

# Create a set of spread values (i.e High - Low)
spread = map(x -> x.high - x.low, googl)

# Compute the average spread
round(reduce(+,(mean.(spread)))/length(spread), digits=4)

# More realistic is the average gain (i.e. Open - Close)
gain = map(x -> x.open - x.close, googl)
round(reduce(+,(mean.(gain)))/length(gain), digits=4)

# Note that dataset is only the FIRST day in the month so the value is not necessarily

# pushcol(googl, :gain, gain)

# sratio = (Open - Close)/(High - Low ) ϵ [-1.0,1.0]
ρ = map(x -> (x.open - x.close) / (x.high - x.low), googl)
round(reduce(+,(mean.(ρ)))/length(ρ),digits=4)

# If this is Weiner (Brownian) process the value should be about 0.5

using OnlineStats

reduce(Mean(), googl; select = :close)

groupreduce(Mean(),stockdata,:ticker; select=:close)

using StatsPlots

@df stockdata plot(:date, :close, group=:ticker, layout = 4, legend = :topleft)

spread[:,:]

volumes = map(x -> x.volume, stockdata)

describe(volumes)

# volumes = map(pick(:volume), stockdata)

reduce(+, map(x -> x.volume, stockdata))

# Does not work in Jupyter
#=
@everywhere function agg_ohlcv(x, y)
    Named_Tuple(
        open = x.open,
        high = maximum(x.high, y.high),
        low = minimum(x.low, y.low),
        close = y.close,
        volume = x.volume + y.volume,
    )
end

reducedim(agg_ohlcv, sampledata, 1)
=#

using MySQL, IndexedTables

conn = MySQL.connect("localhost", "root", "LQaLxxy6",db="Chinook",
            unix_socket="/Applications/MAMP/tmp/mysql/mysql.sock");

sql = "select a.FirstName,a.LastName, sum(b.Total) as 'TotalSpent'";
sql *= " from Customer a";
sql *= " join Invoice b on a.CustomerId = b.CustomerId";
sql *= " group by a.FirstName, a.LastName";

tbl = MySQL.query(conn,sql) |> table

t = table(NamedTuple(x = randn(100), y = randn(100)))
@df t scatter(:x, :y)

path = Pkg.dir(homedir(), "PacktPub", "Chp09")

cd(path)
# diamonds = loadtable("datasets/diamonds.csv"; indexcols = [:carat])
diamonds = loadndsparse("Files/diamonds.csv"; indexcols = [:carat])

t = table(1:100, rand(Bool, 100), randn(100))

reduce(Mean(), t; select = 3)

groupreduce(Mean(), t, 2; select=3)