text:
using Pkg
using GLMakie, Test
using GLMakie.FileIO
using GLMakie.AbstractPlotting
using GLMakie.GeometryBasics
using GLMakie.GeometryBasics: origin
using ImageMagick
# ImageIO seems broken on 1.6 ... and there doesn't
# seem to be a clean way anymore to force not to use a loader library?
filter!(x-> x !== :ImageIO, FileIO.sym2saver[:PNG])
filter!(x-> x !== :ImageIO, FileIO.sym2loader[:PNG])

# run the unit test suite
include("unit_tests.jl")

path = normpath(joinpath(dirname(pathof(AbstractPlotting)), "..", "test", "ReferenceTests"))
Pkg.develop(PackageSpec(path = path))
using ReferenceTests
using ReferenceTests: @cell

# Run the AbstractPlotting reference image testsuite
recorded = joinpath(@__DIR__, "recorded")
rm(recorded; force=true, recursive=true); mkdir(recorded)
ReferenceTests.record_tests(recording_dir=recorded)
ReferenceTests.reference_tests(recorded)

# Run the GLMakie specific backend reference tests
empty!(ReferenceTests.DATABASE)
include("glmakie_tests.jl")
recorded = joinpath(@__DIR__, "recorded_glmakie")
rm(recorded; force=true, recursive=true); mkdir(recorded)
ReferenceTests.record_tests(ReferenceTests.DATABASE, recording_dir=recorded)
ref_images = ReferenceTests.download_refimages(; name="glmakie_refimages")
ReferenceTests.reference_tests(recorded; ref_images=ref_images, difference=0.01)

# needs GITHUB_TOKEN to be defined
# ReferenceTests.upload_reference_images()
# Run the below, to generate a html to view all differences:
# recorded, ref_images, scores = ReferenceTests.reference_tests(recorded)
# ReferenceTests.generate_test_summary("preview.html", recorded, ref_images, scores)
# ReferenceTests.upload_reference_images(recorded; name="glmakie_refimages")


text:
module Elements

using FastGaussQuadrature
import Base.Iterators: product

import ..Transforms: TransformChain

export Simplex, Element, Tensor, quadrule

abstract type ReferenceElement end

struct Simplex{NDim} <: ReferenceElement end

Base.ndims(::Simplex{n}) where n = n

quadrule(::Simplex{0}, npts::Int) = zeros(Float64, 1, 0), [1.0]

function quadrule(::Simplex{N}, npts::Int) where N
    (pts, wts) = gausslegendre(npts)
    pts = (pts .+ 1.0) ./ 2
    wts = wts ./ 2

    pts = hcat([collect(k) for k in product((pts for _ in 1:N)...)]...)
    wts = .*((reshape(wts, fill(1, k)..., npts) for k in 0:N-1)...)
    pts', vec(wts)
end

struct Tensor <: ReferenceElement
    terms :: Array{Simplex}
end

Base.ndims(self::Tensor) = sum(ndims(term) for term in self.terms)

struct Element
    reference :: ReferenceElement
    index :: Int
    transform :: Tuple{TransformChain, TransformChain}
end

Element(ref::ReferenceElement, index::Int; transform::TransformChain=(), dimcorr::TransformChain=()) =
    Element(ref, index, (dimcorr, transform))

Base.ndims(self::Element) = ndims(self.reference)

end # module

ext: jl

max_stars_repo_name: TheBB/Jutils

max_issues_repo_name: TheBB/Jutils

max_forks_repo_name: TheBB/Jutils


text:
# TODO:
# mergefiber - this should probably have interpolative, averaging, and fiber cut options/variations
# combinefibers - ditto

###
### Synapse interface
###
"""
    AbstractSynapse{Pre,Pos,AP,A}

"""
abstract type AbstractSynapse{Pre,Pos,AP,A} end

"""
    Synapse
"""
struct Synapse{Pre,Pos,AP,A} <: AbstractSynapse{Pre,Pos,AP,A}
    presynaptic::Pre
    postsynaptic::Pos
    activityparams::AP
    activity::A
end

"""
    presynaptic(s::AbstractSynapse)
"""
function presynaptic(s::AbstractSynapse) end
presynaptic(s::Synapse{Pre,Pos,AP,A}) where {Pre,Pos,AP,A} = getproperty(s, :presynaptic)::Pre

"""
    postsynaptic(s::AbstractSynapse)
"""
function postsynaptic(s::AbstractSynapse) end
postsynaptic(s::Synapse{Pre,Pos,AP,A}) where {Pre,Pos,AP,A} = getproperty(s, :postsynaptic)::Pos

"""
    activityparams(s::AbstractSynapse)
"""
function activityparams(s::AbstractSynapse) end
activityparams(s::Synapse{Pre,Pos,AP,A}) where {Pre,Pos,AP,A} = getproperty(s, :activityparams)::AP

"""
    activity(::AbstractSynapse)
"""
function activity(s::AbstractSynapse) end
activity(s::Synapse{Pre,Pos,AP,A}) where {Pre,Pos,AP,A} = getproperty(s, :activity)::A

"""
    activate!(s::AbstractSynapse)
"""
activate!(s::AbstractSynapse) = activity(s)(presynaptic(s), postsynaptic(s), activityparams(s))

###
### Fiber interface
###

"""
    AbstractFiberBundle

"""
abstract type AbstractFiberBundle{Nu,Syn,T,P} <: AbstractVector{T} end

"""
    FiberBundle

"""
struct FiberBundle{Nu,Syn,T,P<:AbstractVector{T}} <: AbstractFiberBundle{Nu,Syn,T,P}
    nucleus::Nu
    synapse::Syn
    pathway::P
end

"""
    nucleus(f::AbstractFiberBundle) -> AbstractAnatomy
"""
function nucleus(f::AbstractFiberBundle) end
nucleus(f::FiberBundle{Nu,Syn,T,P}) where {Nu,Syn,T,P} = getproperty(f, :nucleus)::Nu

"""
    synapse(f::AbstractFiberBundle) -> AbstractSynapse
"""
function synapse(f::AbstractFiberBundle) end
synapse(f::FiberBundle{Nu,Syn,T,P}) where {Nu,Syn,T,P} = getproperty(f, :synapse)::Syn

"""
    pathway(f::AbstractFiberBundle) -> AbstractSynapse
"""
function pathway(f::AbstractFiberBundle) end
pathway(f::FiberBundle{Nu,Syn,T,P}) where {Nu,Syn,T,P} = getproperty(f, :pathway)::P

"""
    mergefibers(method, f; kwargs...)

Merge fiber bundle segments along its pathway.
"""
function mergefiber(method, f::FiberBundle; kwargs...) end

"""
    combinefibers(method, f...; kwargs)

Merge fiber bundles together into a single fiber bundle.
"""
function combinefibers(method, f::FiberBundle...; kwargs...) end


text:
using FactCheck

include(Pkg.dir("BioFeatures", "src","Features.jl"))

# makedata

df1 = DataFrame( region_id = ["R1","R1","R2"], start = [1,100,1], stop = [50,300,600] )
df2 = DataFrame( region_id = ["R2","R2","R3"], start = [700,1000,60], stop = [800,1200,500] )
df3 = DataFrame( region_id = ["R4","R4","R5"])
df4 = DataFrame( region_id = ["R10","R11","R12"],  start = [700,1000,60], stop = [800,1200,500] )
F = Features(df1, region_id_field = :region_id)
F_keys = collect( keys( F.features))
F_R2_vals = collect(F.features["R2"])

facts("Initialisation Insertions") do
   @fact length(F_keys) => 2
   @fact length(F_R2_vals)=>1
end

add_dataframe_to_features!(F,df2,region_id_field = :region_id)
F_keys = collect( keys( F.features))
F_R2_vals = collect(F.features["R2"])

facts("Appending Features") do
   @fact length(F_keys) => 3
   @fact length(F_R2_vals)=>3
end

# Adds just regions -- even though no start and stop positions
F1 =  Features( df3 )

# Adds region_ids with uninitialised interval trees
# Adds the region into the description
F2 =  Features(df4, region_id_field=:region_id, regions=true )
F2.description
#facts()


text:
# IMC Monitoring
mutable struct IMCMonitor{T,N} <: AbstractMonitor
    # We have one entry in the outer tuple for each socket.
    # Within each socket, we have one entry for each controller.
    # Within each controller, we have one entry for each channel.
    imcs::Record{:socket,T}
    events::NTuple{N, UncoreSelectRegister}
    cleaned::Bool
end

"""
    IMCMonitor(events, socket; [program = true, finalize = true])

Monitor the Integrated Memory Controller (IMC) for `events` on a **single**
selected CPU `socket`. This can gather information such as number of DRAM read and write
operations.  Argument `event` should be a `Tuple` of [`CounterTools.UncoreSelectRegister`](@ref)
and `socket` should be either an `Integer` or `IndexZero`.

If `finalize = true` is passed, a finalizer will be attached to the `IMCMonitor` to clean
up the hardware counter's state.
"""
function IMCMonitor(
        events::NTuple{N, UncoreSelectRegister},
        socket;
        program = true,
        finalize = true,
    ) where {N}
    # Check to see if we already have an active monitor.
    if IMC_RESERVATION[]
        error("An active IMC Monitor already exists!")
    end

    # Get the mapping from socket to bus
    socket_to_bus = findbusses()

    # Memory Controllers
    imcs = ntuple(2) do controller
        # Channels
        return ntuple(3) do channel
            # Get the path for this particular device
            bus = socket_to_bus[socket]
            device = SKYLAKE_IMC_REGISTERS[controller][channel].device
            fn = SKYLAKE_IMC_REGISTERS[controller][channel].fn

            # Construct a monitor
            handle = Handle(bus, device, fn)
            pmu = IMCUncorePMU(handle)
            reset!(pmu)
            return Record{:channel}((pmu,))
        end |> Record{:imc}
    end |> Record{:socket}

    monitor = IMCMonitor(
        imcs,
        events,
        false,
    )

    IMC_RESERVATION[] = true

    # Clean up after ourselves
    finalize && finalizer(cleanup, monitor)

    program && program!(monitor)
    return monitor
end

const IMC_RESERVATION = Ref{Bool}(false)

# Extend `mapleaves` to broadcast across the IMCMonitor
mapleaves(f, monitor::IMCMonitor) = mapleaves(f, monitor.imcs)
program!(monitor::IMCMonitor) = mapleaves(x -> program(x, monitor.events), monitor)
function program(x::IMCUncorePMU, events)
    for (i, evt) in enumerate(events)
        setcontrol!(x, i, evt)
    end
    enable!(x)
    return nothing
end

reset!(monitor::IMCMonitor) = mapleaves(reset!, monitor)
function Base.read(monitor::IMCMonitor{T,N}) where {T,N}
    monitor.cleaned && error("Trying to measure a cleaned IMC Monitor")
    return mapleaves(getallcounters, monitor)
end

function cleanup(monitor::IMCMonitor)
    if !monitor.cleaned
        # Close any open PCI file handles
        reset!(monitor)
        IMC_RESERVATION[] = false
        mapleaves(close, monitor)
        monitor.cleaned = true
    end
    return nothing
end


text:
@doc """Ericsson Model
freq: signal frequency(range from 100 to 2000Mhz)
tyArea: area type (1-rural, 2-suburban e 3-urban).
cityKind : cyte type(1-small, 2-medium e 3-large).
hb: base station's height
hm: mobile's height
"""
type EricssonModel<:Model
  freq::Float64
  txH::Float64
  rxH::Float64
  cityKind::Int
  checkFreqRange::Bool
  EricssonModel() = new(800,50,1.5,CityKind.Medium,true)
end
@doc "Calculate Path loss for distance dist(km) with the Ericsson Model"
function pathloss_op(model::EricssonModel,dist::AbstractFloat)
    # frequency
  f::Float64  = model.freq;
  #Bts and MS Heights
  hm::Float64 = model.rxH
  hb::Float64 = model.txH
   # Checking the frequency range
   if (model.checkFreqRange)
     if (f < 100 || f > 2000)
      error("The frequency range for Ericcson Model is 500MHz-1500Mhz")
     end
   end
   # Pathloss
   g = 44.49*log10(f)-4.78*(log10(f))^2
   a2= 12
   a3= 0.1
  if (model.cityKind==CityKind.Large)
    a0 = 36.2
    a1 = 30.2;
  elseif (model.cityKind==CityKind.Medium)
    a0 = 43.2
    a1 = 68.9
  else
    a0 = 45.9
    a1 = 100.6
  end
  PL=a0+a1*log10(dist)+a2*log10(hb)+a3*(log10(hb))*(log10(dist))-3.2*log10((11.75*hm)^2)+g
  return PL
end


text:
# Test GLS
using PDESolver
using FactCheck
using ODLCommonTools
using ArrayViews
using Input

# insert a command line argument
fname = "input_vals_vortex_GLS.jl"
opts = read_input(fname)
@testset "Check input arguments applied correctly" begin
  @test ( opts["use_GLS"] )== true
  @test ( opts["Q_transpose"] )== true
end

mesh, sbp, eqn, opts = solvePDE(fname)

fill!(eqn.res, 0.0)
fill!(eqn.res_vec,0.0)
# res_0 = zeros(eqn.res_vec)
res_0_norm = physicsRhs(mesh, sbp, eqn, opts, eqn.res_vec, (evalResidual,))

# println("eqn.q = \n", eqn.q)
# println("dxidx = \n", mesh.dxidx)
# println("jacobian determinant = \n", mesh.jac)
# println("coordinates = \n", mesh.coords)
# println("eqn.Axi[:,:,1,1] = \n", eqn.Axi[:,:,1,1])

@testset "\nCheck flux jacobian for conservative variables" begin

  @test isapprox(eqn.Axi[:,:,1,1],   [0.0 0.5 0.0 0.0
                                      0.045125 0.0 0.13435028842 0.2
                                      0.0 -0.335875721 0.0 0.0
                                      0.0 0.7378125 0.0 0.0]) atol=1e-13
  @test isapprox(eqn.Axi[:,:,2,1],   [0.0 0.5 0.0 0.0
                                      0.01128125 0.0 0.06717514421272198 0.2
                                      0.0 -0.167937860531805 0.0 0.0
                                      0.0 0.7378125 0.0 0.0]) atol=1e-13
  @test isapprox(eqn.Axi[:,:,3,1] ,  [0.0 0.5 0.0 0.0
                                     -0.03384375 0.26870057685088766 0.0671751442 0.2
                                      0.05640625 -0.16793786053 0.16793786053 0.0
                                     -0.240235109491 0.71525 0.0225625 0.2351130]) atol-1e-13
  @test isapprox(eqn.Aeta[:,:,1,1],   [0.0 0.0 0.5 0.0
                                       0.0 -0.3358757210636101 0.0 0.0
                                      -0.1805 0.0 -0.5374011537017761 0.2
                                       0.4653138270684989 0.0 0.6475625 -0.47022601]) atol=1e-13
  @test isapprox(eqn.Aeta[:,:,2,1],   [0.0 0.0 0.5 0.0
                                       0.0 -0.16793786053180498 0.0 0.0
                                      -0.045125 0.0 -0.268700576850888 0.2
                                       0.2440242074689959 0.0 0.71525 -0.23511300474452695]) atol=1e-13
  @test isapprox(eqn.Aeta[:,:,3,1],   [0.0 0.0 0.5 0.0
                                       0.05640625 -0.16793786053180487 0.16793786053180484 0.0
                                      -0.03384375 -0.06717514421272192 -0.2687005768508878 0.2
                                       0.24023510949074675 0.0225625 0.71525 -0.2351130047445268]) atol=1e-13

  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      Flux_xi = eqn.Axi[:,:,j,i]*eqn.q[:,j,i]
      Flux_eta = eqn.Aeta[:,:,j,i]*eqn.q[:,j,i]
      @test isapprox( Flux_xi, eqn.flux_parametric[:,j,i,1])
      @test isapprox( Flux_eta, eqn.flux_parametric[:,j,i,2])
    end
  end
end

fill!(eqn.res, 0.0)
fill!(eqn.res_vec,0.0)

parametricFluxJacobian(mesh, sbp, eqn)

GLS(mesh, sbp, eqn)
@testset "Check if GLS residual is getting added to the residual" begin
  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      for k = 1:mesh.numDofPerNode
        @test  abs(eqn.res[k,j,i])  > 0.0
      end
    end
  end
end

# println("Ax = \n", Ax)
# println("\nAy = \n", Ay)

@testset "Check if the Flux Jacobian in the X & Y direction are being computed accrately" begin
  q = ones(Tsol, 4)
  Ax = zeros(Tsol, 4, 4)
  Ay = zeros(Tsol, 4, 4)
  calcFluxJacobian(eqn, q, Ax, Ay)

  @test isapprox(Ax,   [0.0 1.0 0.0  0.0
                       -0.6 1.6 -0.4 0.4
                       -1.0 1.0 1.0 0.0
                       -0.6 0.6 -0.4 1.4]) atol=1e-13
  @test isapprox(Ay,   [0.0 0.0 1.0 0.0
                       -1.0 1.0 1.0 0.0
                       -0.6 -0.4 1.6 0.4
                       -0.6 -0.4 0.6 1.4]) atol=1e-13
end

@testset "Check if eigen value factorization is being computed correctly" begin
  q = Complex{Float64}[2.0,0.0,-1.3435028842544403,2.236964285714286]
  T = zeros(Tsol,4,4)
  Tinv = zeros(T)
  Lambda = eye(T)

  # check Axi factorization is being computed correctly
  dxidx = Tsol[0.5,0.0]
  calcEigenFactorization(eqn, q, dxidx, T, Tinv, Lambda)
  @test isapprox(T,   [1.0 0.0 1.0 1.0
                       0.0 0.0 0.7071067811865476 -0.7071067811865476
                      -0.67175144212722 -1.0 -0.6717514421272202 -0.6717514421
                       0.22562500000000002 0.6717514 1.4756250 1.4756250]) atol=1e-13
  @test isapprox(Tinv,   [0.8195 0.0 -0.5374011537017759 -0.8
                         -0.6717514421272202 0.0 -1.0 0.0
                          0.09025 0.7071067811865475 0.26870057685088794 0.4
                          0.09025 -0.7071067811865475 0.26870057685088794 0.4]) atol=1e-13
  @test isapprox(Lambda,   [0.0 0.0 0.0 0.0
                            0.0 0.0 0.0 0.0
                            0.0 0.0 0.3535533905932738 0.0
                            0.0 0.0 0.0 -0.3535533905932738]) atol=1e-13
  @test isapprox( T*Lambda*Tinv, eqn.Axi[:,:,1,1])

  # Check Aeta factorization in being computed correctly
  dxidx = Tsol[0.0,0.5]
  calcEigenFactorization(eqn, q, dxidx, T, Tinv, Lambda)
  @test isapprox(T,   [1.0 0.0 1.0 1.0
                       0.0 1.0 0.0 0.0
                      -0.6717514421272202 -0.0 0.035355339 -1.3788582233
                       0.225625 0.0 1.0006250000000003 1.950625]) atol=1e-13
  @test isapprox(Tinv,   [0.8195 0.0 -0.5374011537017759 -0.8
                         -0.0 1.0 -0.0 0.0
                          0.56525 0.0 0.9758073580374353 0.4
                         -0.38475 -0.0 -0.43840620433565947 0.4]) atol=1e-13
  @test isapprox(Lambda,    [-0.3358757210636101 0.0 0.0 0.0
                             0.0 -0.3358757210636101 0.0 0.0
                             0.0 0.0 0.01767766952966371 0.0
                             0.0 0.0 0.0 -0.6894291116568838]) atol=1e-13
   @test isapprox( T*Lambda*Tinv, eqn.Aeta[:,:,1,1])
end

@testset "Check if element level product of shape function with nodal flux jacobian is being computed correctly" begin

  ndof = mesh.numDofPerNode
  nnpe = mesh.numNodesPerElement
  endof = nnpe*ndof # degrees of freedom in an element
  Axidxi = zeros(Tsol, endof, endof)
  shapefuncderiv = ones(Tsol, sbp.numnodes, sbp.numnodes, Tdim)
  Axi = ones(Tsol, ndof, ndof, nnpe)
  Aeta = ones(Axi)
  calcAxidxi(Axidxi, shapefuncderiv, Axi, Aeta, ndof, nnpe)

  @test isapprox( Axidxi, 2*ones(Tsol, endof, endof))

end


text:
using Distances
using LinearAlgebra
using SparseArrays
import Base: size
include("../tools/plotting.jl")

#
#
# CONVENTIONS ON x AND y:
# locations x and y will always be arrays, whether it's 1D or 2D
#
# NOTES ABOUT DISTANCES.JL
# - WeightedSqEuclidean([3, 4])(x, y) from Distances.jl will default to SqEuclidean if x and y aren't arrays
#   it will throw an error of x and y are arrays but not the same length as the weights vector
# - one can input single length scale as scalar or array of length 1
#   and multiple length scale as array
# - SqEuclidean and other distance functions from Distances.jl only take in 2-D arrays as inputs
#
# HOW TO DETECT MULTIPLE LENGTH SCALES?
# - when θ is passed as a Real, then distance will default to SqEuclidean (single length scale)
# - when θ is passed as an Array, then distance will default to WeightedSqEuclidean (multiple length scales)
#
#
#
#
@doc raw"""
    AbstractCorrelation

An abstract type for all correlation functions.
"""
abstract type AbstractCorrelation end

(k::AbstractCorrelation)(x::Array, y::Array, θ) = k(distance(k, θ)(x, y),  θ)

"""
Computes matrix of pairwise distances between x and y, which contains points
arranged row-wise
"""
function distancematrix(k::AbstractCorrelation, θ, x, y=x; dims=1)
    ret = Array{Float64}(undef, size(x, dims), size(x, dims))
    computeDists!(ret, k, θ, x, y; dims=dims)
    return ret
end

"""
x, y: arrays of points in R^d
θ: 1D array
"""
function computeDists!(out, k::AbstractCorrelation, θ, x, y=x; dims=1)
    if typeof(θ)<:Array
        @assert max(size(θ, 1), size(θ, 2)) == size(x, 2) == size(y, 2)
    end
    dist = distance(k, θ)
    try
        out .= pairwise!(out, dist, x, y, dims = dims)
    catch(MethodError)
        out .= pairwise!(out, dist, reshape(x,size(x, 1), size(x, 2)), reshape(y, size(y, 1), size(y, 2)), dims = dims)
    end
        return nothing
end

"""
jitter: added to diagonal of kernel matrix to ensure positive-definiteness
dims: 1 if data points are arranged row-wise and 2 if col-wise

"""
function correlation(k::AbstractCorrelation, θ, x; jitter = 1e-12, dims=1)
    ret = Array{Float64}(undef, size(x, dims), size(x, dims))
    correlation!(ret, k, θ, x, jitter = jitter)
    return ret
end

function cross_correlation(k::AbstractCorrelation, θ, x, y; dims=1)
    ret = Array{Float64}(undef, size(x, dims), size(y, dims))
    cross_correlation!(ret, k, θ, x, y;dims=dims)
    return ret
end

function cross_correlation!(out, k::AbstractCorrelation, θ, x, y; dims = 1)
    x = reshape(x, size(x, 1), size(x, 2)); y = reshape(y, size(y, 1), size(y, 2)) #2D arrays
    dist = distance(k, θ)
    pairwise!(out, dist, x, y, dims=dims)
    out .= (τ -> k(τ, θ)).(out)
    return nothing
end

function correlation!(out, k::AbstractCorrelation, θ, x; jitter = 1e-6, dims=1)
    x = reshape(x, size(x, 1), size(x, 2))
    dist = distance(k, θ)
    pairwise!(out, dist, x, dims=dims)
    out .= (τ -> k(τ, θ)).(out)
    if jitter != 0
        out[diagind(out)] .+= jitter
        # out ./= out[1, 1] #covariance must be in [0, 1]
    end
    return nothing
end

@doc raw"""
"""
struct FixedParam{K<:AbstractCorrelation,T} <: AbstractCorrelation
    k::K
    θ::T
end
FixedParam(k, θ...) = FixedParam(k, θ)

(k::FixedParam)(τ) = k.k(τ, k.θ...)
distance(k::FixedParam) = distance(k.k, k.theta...)

@doc raw"""
    Gaussian
"""
struct Gaussian <: AbstractCorrelation end
const RBF = Gaussian
const SqExponential = RBF

#WARNING: these two functions are purely meant to be used in conjunction with (k::AbstractCorrelation)(x::Array, y::Array, θ)
#These are not standalone functions!
(::Gaussian)(τ::Real, θ::Real) = exp(- τ * θ / 2)
(::Gaussian)(τ::Real, ::Union{AbstractVector, Array{T}} where T) = exp(- τ / 2) #theta already taken into account in computation of tau

distance(::Gaussian, θ::Real) = SqEuclidean()
distance(::Gaussian, θ::Union{AbstractVector, Array{T}} where T) = size(θ, 2)>1 ? WeightedSqEuclidean(reshape(θ, length(θ))) : WeightedSqEuclidean(θ)

"""
currently derivative is only implemented for single length scale
"""
function partial_θ(k::Gaussian, τ, θ)
    -τ/2 * k(τ, θ[1])
end
#partial_θ(k::Gaussian, τ, θ...) = ;
#partial_τ(k::Gaussian, τ, θ) = - θ/2 * k(τ, θ)
#partial_θθ(k::RBF, τ, θ) = τ^2/4 * k(τ, θ)
#partial_ττ(k::Gaussian, τ, θ) =  θ^2/4 * k(τ, θ)

struct Spherical <: AbstractCorrelation end
# TODO

struct Matern <: AbstractCorrelation end
# TODO

struct RationalQuadratic <: AbstractCorrelation end
# TODO

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module LossFunctionsModule

import Random: randperm
import LossFunctions: value, AggMode, SupervisedLoss
import ..CoreModule: Options, Dataset, Node, string_tree
import ..EquationUtilsModule: compute_complexity, thermoConstraints
import ..EvaluateEquationModule: eval_tree_array, differentiable_eval_tree_array

# for symbolic constraints
using PyCall
const _sympy = PyNULL()

function loss(
    x::AbstractArray{T}, y::AbstractArray{T}, options::Options{A,B,dA,dB,C,D}
)::T where {T<:Real,A,B,dA,dB,C<:SupervisedLoss,D}
    return value(options.loss, y, x, AggMode.Mean())
end
function loss(
    x::AbstractArray{T}, y::AbstractArray{T}, options::Options{A,B,dA,dB,C,D}
)::T where {T<:Real,A,B,dA,dB,C<:Function,D}
    return sum(options.loss.(x, y)) / length(y)
end

function loss(
    x::AbstractArray{T},
    y::AbstractArray{T},
    w::AbstractArray{T},
    options::Options{A,B,dA,dB,C,D},
)::T where {T<:Real,A,B,dA,dB,C<:SupervisedLoss,D}
    return value(options.loss, y, x, AggMode.WeightedMean(w))
end
function loss(
    x::AbstractArray{T},
    y::AbstractArray{T},
    w::AbstractArray{T},
    options::Options{A,B,dA,dB,C,D},
)::T where {T<:Real,A,B,dA,dB,C<:Function,D}
    return sum(options.loss.(x, y, w)) / sum(w)
end

# Evaluate the loss of a particular expression on the input dataset.
function eval_loss(tree::Node, dataset::Dataset{T}, options::Options)::T where {T<:Real}
    (prediction, completion) = eval_tree_array(tree, dataset.X, options)
    if !completion
        return T(1000000000)
    end

    if dataset.weighted
        return loss(prediction, dataset.y, dataset.weights, options)
    else
        return loss(prediction, dataset.y, options)
    end
end

# Compute a score which includes a complexity penalty in the loss
function loss_to_score(
    loss::T, baseline::T, tree::Node, options::Options
)::T where {T<:Real}
    normalized_loss_term = loss / baseline
    size = compute_complexity(tree, options)
    parsimony_term = size * options.parsimony

    return normalized_loss_term + parsimony_term
end

# Score an equation
function score_func(
    dataset::Dataset{T}, baseline::T, tree::Node, options::Options
)::Tuple{T,T} where {T<:Real}
    result_loss = eval_loss(tree, dataset, options)

    # thermo constraints
    if options.penalties !== nothing
        if ispynull(_sympy)
            println("loading lib")
        end
        # good way to check if lib needs to be loaded from PyCall/serialize.jl
        sympy = ispynull(_sympy) ? copy!(_sympy, pyimport_conda("sympy", "sympy")) : _sympy

        if length(dataset.varMap) >= 1 && ("p" in dataset.varMap)
            # turn tree to string for sympy to use
            # get boolean result for each constraint
            str = string_tree(tree, options)
            expr = sympy.parse_expr(str)
            #expr = sympyVar.parse_expr("p*2 + 1")
            var = sympy.symbols("p")
            results = thermoConstraints(expr, var)
            #results = [true, true, true]
            # for each constraint, apply penalty if expression fails
            for i in 1:length(options.penalties)
                if results[i] == false
                    result_loss *= options.penalties[i]
                end
            end
        else
            println("Need pressure variable \"p\"!")
        end
    end

    score = loss_to_score(result_loss, baseline, tree, options)

    return score, result_loss
end

# Score an equation with a small batch
function score_func_batch(
    dataset::Dataset{T}, baseline::T, tree::Node, options::Options
)::Tuple{T,T} where {T<:Real}
    batch_idx = randperm(dataset.n)[1:(options.batchSize)]
    batch_X = dataset.X[:, batch_idx]
    batch_y = dataset.y[batch_idx]
    (prediction, completion) = eval_tree_array(tree, batch_X, options)
    if !completion
        return T(1000000000), T(1000000000)
    end

    if !dataset.weighted
        result_loss = loss(prediction, batch_y, options)
    else
        batch_w = dataset.weights[batch_idx]
        result_loss = loss(prediction, batch_y, batch_w, options)
    end
    score = loss_to_score(result_loss, baseline, tree, options)
    return score, result_loss
end

end


text:
include("abstractfeatures.jl")
using DataFrames
using ..ChemistryFeaturization: elements
using ..ChemistryFeaturization.AbstractType: AbstractCodec, AbstractAtoms
using ..ChemistryFeaturization.Codec: SimpleCodec
using ..ChemistryFeaturization.Utils.OrbitalFeatureUtils:
    _indexorbital, _econf_to_name, _orbitalsparse
using ..ChemistryFeaturization.Data: valenceshell_conf_df
using SparseArrays

"""
    OrbitalFeatureDescriptor

Describes the orbital configuration of an element.
"""
struct OrbitalFeatureDescriptor <: AbstractEnvironmentFeatureDescriptor
    encoder_decoder::AbstractCodec
    function OrbitalFeatureDescriptor()
        # trim down the DataFrame to have only the Symbol, and valence shell configuration (extracted from "Electronic Structure")
        new(SimpleCodec(default_ofd_encode, default_ofd_decode))
    end
end

function get_value(ofd::OrbitalFeatureDescriptor, el::String)
    @assert el in valenceshell_conf_df[:, :Symbol] "All elements must be valid and accounted for in the periodic table!"

    # pull the value for this element from the "Electronic Structure" column of the orbital data table
    getproperty(
        valenceshell_conf_df[valenceshell_conf_df.Symbol.==el, :][1, :],
        Symbol("Electronic Structure"),
    )
end

get_value(ofd::OrbitalFeatureDescriptor, a::AbstractAtoms) =
    map(e -> get_value(ofd, e), elements(a))

(ofd::OrbitalFeatureDescriptor)(el::String) = get_value(ofd, el)

function output_shape(ofd::OrbitalFeatureDescriptor, sc::SimpleCodec)
    # return the number
end

Base.show(io::IO, efd::OrbitalFeatureDescriptor) =
    print(io, "OrbitalFeatureDescriptor created.")

function Base.show(io::IO, ::MIME"text/plain", efd::OrbitalFeatureDescriptor)
    st = "OrbitalFeatureDescriptor created."
    print(io, st)
end

"""
    default_ofd_encode(ofd::OrbitalFeatureDescriptor, a::AbstractAtoms)

Default encoding scheme for OrbitalFeatureDescriptor, which returns a
sparse matrix, where column_i is a sparse representation of the electronic
configuration of the ith element in the [Atoms](@ref atoms) object.
"""
function default_ofd_encode(elec_configs::Vector{String})
    I = Vector{Int16}()
    J = Vector{Int16}()
    V = Vector{Int16}()
    col = 0
    for s in elec_configs
        # X, Y are equivalent to I, V for a SparseVector
        X, Y = _orbitalsparse(s)
        col += 1    # column number

        for i = 1:length(X)
            append!(I, X[i])
            append!(J, col)
            append!(V, Y[i])
        end
    end
    sparse(I, J, V) # create a SparseMatrix with the valence shell configuration of each element as a column vector
end

"""
    default_ofd_decode(
        ofd::OrbitalFeatureDescriptor,
        encoded_features::SparseArrays.AbstractSparseMatrixCSC{Tv,Ti}
    ) where {Tv,Ti}

Default decoding scheme for OrbitalFeatureDescriptor, which returns the
vector of elements encoded, given the encoded sparse matrix.
"""
function default_ofd_decode(
    encoded_features::SparseArrays.AbstractSparseMatrixCSC{Tv,Ti},
) where {Tv,Ti}
    elements = String[]
    for i = 1:size(encoded_features)[2]
        push!(elements, _econf_to_name(encoded_features[:, i]))
    end
    return elements
end


text:
#!/usr/bin/julia
#
# Helper functions for constructing the Hamiltonian.
#
export sec_dv, energy, ke, pe;

"""
    sec_dv(grid::UniformGrid{S,1}) where S

Construct matrix representation of discretised second order partial derivative
acting on position space.
"""
function sec_dv(grid::UniformGrid{S,1}) where S
  rows = length(grid);
  # Second order method.
  #return LinearAlgebra.SymTridiagonal(repeat([-2.0/grid.Δ^2], rows),
  #                      repeat([1.0/grid.Δ^2], rows-1));

  # Second order method explicitly constructed.
  M = zeros(Float64, rows, rows);
  for i ∈ 1:rows
    # Set incomplete end values to zero to prevent introduction of spurious
    # boundary terms
    if i > 1 && i <= rows-1
      M[i,i-1] = -1;
    #end
    M[i,i] = 2;
    #if i <= rows-1
      M[i,i+1] = -1;
    end
  end
  #M *= -Rational(1/grid.Δ^2);
  M *= -(1/grid.Δ^2);
  return M;

  #=
  # Infinite order expression for 2nd derivative as per
  #J. Chem. Phys. 96, 1982 (1992) Daniel T. Colbert and William H. Miller
  # https://aip-scitation-org.ezphost.dur.ac.uk/doi/10.1063/1.462100
  M = Matrix{Float64}(undef, rows, rows);
  for i ∈ 1:rows
    for j ∈ 1:rows
      if (i == j)
        M[i,i] = π^2/3;
      else
        M[i,j] = (-1)^(i-j) * 2/(i-j)^2;
      end
    end
  end
  M *= -1/grid.Δ^2;
  return M;=#

  # 5 point function
  #M = Matrix{Rational}(zero(Rational), rows, rows);
  #M = zeros(big(Rational), rows, rows);
  #=
  M = zeros(Float64, rows, rows);
  for i ∈ 3:rows-3
    M[i,i-2] = 1/12;
    M[i,i+2] = 1/12;
    M[i,i-1] = -4/3;
    M[i,i+1] = -4/3;
    M[i,i] = 5/2;
    #=
    if i > 1
      M[i,i-1] = -4//3;
      if (i > 2) M[i,i-2] = 1//12; end
    end
    M[i,i] = 5//2;
    if i <= rows-1
      M[i,i+1] = -4//3;
      if (i <= rows-2) M[i,i+2] = 1//12; end
    end=#
  end
  #M *= -Rational(1/grid.Δ^2);
  M *= -(1/grid.Δ^2);
  return M;=#
end

"""
    energy(ψ::WaveFunction{T,1,S,G}, H::SeparableHamiltonian) where
            {T,S,G<:UniformGrid}

Calculate energy of the given wave function with the kinetic energy calculated
in momentum space (uses a fast Fourier transform if wave function is in
configuration space).

See also: [`ke`](@ref), [`pe`](@ref)
"""
function energy(ψ::WaveFunction{T,1,S,G}, H::SeparableHamiltonian) where
            {T,S,G<:UniformGrid}
  return ke(ψ, H) + pe(ψ,H);
end
#In the future this function may implement a different underlying energy calculation.
function energy(ψ::WaveFunction{T,1,S,G}, H::SeparableHamiltonian) where
            {T,S<:ConfigurationSpace,G<:UniformGrid}
  return ke(ψ, H) + pe(ψ, H);
end

"""
    energy(ψ::WaveFunction{T,1,S,G}, H::SeparableHamiltonian,
           g::UniformGrid{S,1}) where {T,S<:ConfigurationSpace,G<:UniformGrid}

Calculate energy of the given wave function with the kinetic energy determined
using matrix multiplication in configuration space.

See also: [`ke`](@ref), [`pe`](@ref)
"""
function energy(ψ::WaveFunction{T,1,S,G}, H::SeparableHamiltonian,
                g::UniformGrid{S,1}) where {T,S<:ConfigurationSpace,G<:UniformGrid}
  return ke(ψ, g) + pe(ψ, H);
end
# TODO: Implement method of finding energy for non-separable Hamiltonians.

"""
    ke(ψ::WaveFunction{T,1,S,G}, H::SeparableHamiltonian) where
            {T,S<:MomentumSpace,G<:UniformGrid}

Compute kinetic energy fully in momentum space.

Uses a one-dimensional wave function defined on a `UniformGrid`.
"""
function ke(ψ::WaveFunction{T,1,S,G}, H::SeparableHamiltonian) where
            {T,S<:MomentumSpace,G<:UniformGrid}
  return real(sum(conj(ψ) .* H.T .* ψ)) * ψ.Δ;
end

"""
    ke(ψ::WaveFunction{T,1,S,G}, g::UniformGrid{S,1}) where
      {T,S<:ConfigurationSpace,G<:UniformGrid}

Compute kinetic energy in configuration space using E_k = \frac{1}{2} p^2
where the wave function is defined on a one-dimensional `UniformGrid`.
"""
function ke(ψ::WaveFunction{T,1,S,G}, g::UniformGrid{S,1}) where
            {T,S<:ConfigurationSpace,G<:UniformGrid}
  return real(ψ' * -0.5sec_dv(g) * ψ) * ψ.Δ;
end

"""
    ke(ψ::WaveFunction{T,1,S,G}, H::SeparableHamiltonian) where
      {T,S<:ConfigurationSpace,G<:UniformGrid}

Compute kinetic energy using a Fourier transform to momentum space and then
using the kinetic component of Hamiltonian.

The wave function will be defined on a one-dimensional `UniformGrid`.
"""
function ke(ψ::WaveFunction{T,1,S,G}, H::SeparableHamiltonian) where
            {T,S<:ConfigurationSpace,G<:UniformGrid}
  return real(sum(conj(ψ) .* FFTW.ifft(H.T .* FFTW.fft(ψ)))) * ψ.Δ;
end
function ke(ψ::WaveFunction{T,1,S,G}, H::SeparableHamiltonian{Q}) where
            {T<:Unitful.Quantity,S<:ConfigurationSpace,G<:UniformGrid,Q<:Unitful.Quantity}
  return ustrip(real(sum(conj(ψ) .* FFTW.ifft(ustrip(H.T) .* FFTW.fft(ψ)))) * ψ.Δ) * unit(Q);
end

"""
    pe(ψ::WaveFunction{T,1,S,G}, H::SeparableHamiltonian) where
            {T,S<:ConfigurationSpace,G<:UniformGrid}

Calculate potential energy of wave function.

Wave function is in configuration space and defined on one-dimensional uniform
grid.
"""
function pe(ψ::WaveFunction{T,1,S,G}, H::SeparableHamiltonian) where
            {T,S<:ConfigurationSpace,G<:UniformGrid}
  # Calculate energy due to spatial potential.
  Ex = real(sum(conj(ψ) .* H.V .* ψ)) * ψ.Δ;
end

# TODO: Implement determination of P.E. for momentum space wave function.


text:
#! format: off
get_variable_multiplier(_, ::Type{PSY.RegulationDevice{PSY.ThermalStandard}}, ::DeviceLimitedRegulation) = NaN
############################ DeltaActivePowerUpVariable, RegulationDevice ###########################

get_variable_binary(::DeltaActivePowerUpVariable, ::Type{<:PSY.RegulationDevice}, ::AbstractRegulationFormulation) = false
get_variable_lower_bound(::DeltaActivePowerUpVariable, ::PSY.RegulationDevice, ::AbstractRegulationFormulation) = 0.0

############################ DeltaActivePowerDownVariable, RegulationDevice ###########################

get_variable_binary(::DeltaActivePowerDownVariable, ::Type{<:PSY.RegulationDevice}, ::AbstractRegulationFormulation) = false
get_variable_lower_bound(::DeltaActivePowerDownVariable, ::PSY.RegulationDevice, ::AbstractRegulationFormulation) = 0.0

############################ AdditionalDeltaActivePowerUpVariable, RegulationDevice ###########################

get_variable_binary(::AdditionalDeltaActivePowerUpVariable, ::Type{<:PSY.RegulationDevice}, ::AbstractRegulationFormulation) = false
get_variable_lower_bound(::AdditionalDeltaActivePowerUpVariable, ::PSY.RegulationDevice, ::AbstractRegulationFormulation) = 0.0

############################ AdditionalDeltaActivePowerDownVariable, RegulationDevice ###########################

get_variable_binary(::AdditionalDeltaActivePowerDownVariable, ::Type{<:PSY.RegulationDevice}, ::AbstractRegulationFormulation) = false
get_variable_lower_bound(::AdditionalDeltaActivePowerDownVariable, ::PSY.RegulationDevice, ::AbstractRegulationFormulation) = 0.0

get_multiplier_value(::ActivePowerTimeSeriesParameter, d::PSY.RegulationDevice, _)  = PSY.get_max_active_power(d)

########################Objective Function##################################################
proportional_cost(::PSY.OperationalCost, ::AdditionalDeltaActivePowerUpVariable, d::PSY.RegulationDevice, ::AbstractRegulationFormulation) = isapprox(PSY.get_participation_factor(d).up, 0.0; atol=1e-2) ? SERVICES_SLACK_COST : 1 / PSY.get_participation_factor(d).up
proportional_cost(::PSY.OperationalCost, ::AdditionalDeltaActivePowerDownVariable, d::PSY.RegulationDevice, ::AbstractRegulationFormulation) = isapprox(PSY.get_participation_factor(d).dn, 0.0; atol=1e-2) ? SERVICES_SLACK_COST : 1 / PSY.get_participation_factor(d).dn

objective_function_multiplier(::VariableType, ::AbstractRegulationFormulation)=OBJECTIVE_FUNCTION_POSITIVE

#! format: on

function get_default_time_series_names(
    ::Type{<:PSY.RegulationDevice{T}},
    ::Type{<:AbstractRegulationFormulation},
) where {T <: PSY.StaticInjection}
    return Dict{Type{<:TimeSeriesParameter}, String}(
        ActivePowerTimeSeriesParameter => "max_active_power",
    )
end

function get_default_attributes(
    ::Type{<:PSY.RegulationDevice{T}},
    ::Type{<:AbstractRegulationFormulation},
) where {T <: PSY.StaticInjection}
    return Dict{String, Any}()
end

function add_constraints!(
    container::OptimizationContainer,
    ::Type{S},
    ::Type{DeltaActivePowerUpVariable},
    devices::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, DeviceLimitedRegulation},
    ::Type{AreaBalancePowerModel},
) where {
    S <: RegulationLimitsConstraint,
    T <: PSY.RegulationDevice{U},
} where {U <: PSY.StaticInjection}
    var_up = get_variable(container, DeltaActivePowerUpVariable(), T)
    var_dn = get_variable(container, DeltaActivePowerDownVariable(), T)
    base_points_param = get_parameter(container, ActivePowerTimeSeriesParameter(), T)
    multiplier = get_multiplier_array(base_points_param)
    base_points = get_parameter_array(base_points_param)

    names = [PSY.get_name(g) for g in devices]
    time_steps = get_time_steps(container)

    container_up =
        add_constraints_container!(container, S(), U, names, time_steps, meta="up")
    container_dn =
        add_constraints_container!(container, S(), U, names, time_steps, meta="dn")

    for d in devices
        name = PSY.get_name(d)
        limits = PSY.get_active_power_limits(d)
        for t in time_steps
            container_up[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                var_up[name, t] <= limits.max - base_points[name, t] * multiplier[name, t]
            )
            container_dn[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                var_dn[name, t] <= base_points[name, t] * multiplier[name, t] - limits.min
            )
        end
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    ::Type{S},
    ::Type{DeltaActivePowerUpVariable},
    devices::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, ReserveLimitedRegulation},
    ::Type{AreaBalancePowerModel},
) where {
    S <: RegulationLimitsConstraint,
    T <: PSY.RegulationDevice{U},
} where {U <: PSY.StaticInjection}
    var_up = get_variable(container, DeltaActivePowerUpVariable(), T)
    var_dn = get_variable(container, DeltaActivePowerDownVariable(), T)

    names = [PSY.get_name(g) for g in devices]
    time_steps = get_time_steps(container)

    container_up =
        add_constraints_container!(container, S(), U, names, time_steps, meta="up")
    container_dn =
        add_constraints_container!(container, S(), U, names, time_steps, meta="dn")

    for d in devices
        name = PSY.get_name(d)
        limit_up = PSY.get_reserve_limit_up(d)
        limit_dn = PSY.get_reserve_limit_dn(d)
        for t in time_steps
            container_up[name, t] =
                JuMP.@constraint(container.JuMPmodel, var_up[name, t] <= limit_up)
            container_dn[name, t] =
                JuMP.@constraint(container.JuMPmodel, var_dn[name, t] <= limit_dn)
        end
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    ::Type{S},
    devices::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, DeviceLimitedRegulation},
    ::Type{AreaBalancePowerModel},
) where {
    S <: RampLimitConstraint,
    T <: PSY.RegulationDevice{U},
} where {U <: PSY.StaticInjection}
    R_up = get_variable(container, DeltaActivePowerUpVariable(), T)
    R_dn = get_variable(container, DeltaActivePowerDownVariable(), T)

    resolution = Dates.value(Dates.Minute(get_resolution(container)))
    names = [PSY.get_name(g) for g in devices]
    time_steps = get_time_steps(container)

    container_up =
        add_constraints_container!(container, S(), U, names, time_steps, meta="up")
    container_dn =
        add_constraints_container!(container, S(), U, names, time_steps, meta="dn")

    for d in devices
        ramp_limits = PSY.get_ramp_limits(d)
        ramp_limits === nothing && continue
        scaling_factor = resolution * SECONDS_IN_MINUTE
        name = PSY.get_name(d)
        for t in time_steps
            container_up[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                R_up[name, t] <= ramp_limits.up * scaling_factor
            )
            container_dn[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                R_dn[name, t] <= ramp_limits.down * scaling_factor
            )
        end
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    ::Type{S},
    devices::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, <:AbstractRegulationFormulation},
    ::Type{AreaBalancePowerModel},
) where {
    S <: ParticipationAssignmentConstraint,
    T <: PSY.RegulationDevice{U},
} where {U <: PSY.StaticInjection}
    time_steps = get_time_steps(container)
    R_up = get_variable(container, DeltaActivePowerUpVariable(), T)
    R_dn = get_variable(container, DeltaActivePowerDownVariable(), T)
    R_up_emergency = get_variable(container, AdditionalDeltaActivePowerUpVariable(), T)
    R_dn_emergency = get_variable(container, AdditionalDeltaActivePowerUpVariable(), T)
    area_reserve_up = get_variable(container, DeltaActivePowerUpVariable(), PSY.Area)
    area_reserve_dn = get_variable(container, DeltaActivePowerDownVariable(), PSY.Area)

    component_names = [PSY.get_name(d) for d in devices]
    participation_assignment_up = add_constraints_container!(
        container,
        S(),
        T,
        component_names,
        time_steps,
        meta="up",
    )
    participation_assignment_dn = add_constraints_container!(
        container,
        S(),
        T,
        component_names,
        time_steps,
        meta="dn",
    )

    expr_up = get_expression(container, EmergencyUp(), PSY.Area)
    expr_dn = get_expression(container, EmergencyDown(), PSY.Area)
    for d in devices
        name = PSY.get_name(d)
        services = PSY.get_services(d)
        if length(services) > 1
            device_agc = (a for a in PSY.get_services(d) if isa(a, PSY.AGC))
            area_name = PSY.get_name.(PSY.get_area.(device_agc))[1]
        else
            device_agc = first(services)
            area_name = PSY.get_name(PSY.get_area(device_agc))
        end
        p_factor = PSY.get_participation_factor(d)
        for t in time_steps
            participation_assignment_up[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                R_up[name, t] ==
                (p_factor.up * area_reserve_up[area_name, t]) + R_up_emergency[name, t]
            )
            participation_assignment_dn[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                R_dn[name, t] ==
                (p_factor.dn * area_reserve_dn[area_name, t]) + R_dn_emergency[name, t]
            )
            JuMP.add_to_expression!(expr_up[area_name, t], -1 * R_up_emergency[name, t])
            JuMP.add_to_expression!(expr_dn[area_name, t], -1 * R_dn_emergency[name, t])
        end
    end

    return
end

function objective_function!(
    container::OptimizationContainer,
    devs::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, U},
) where {
    T <: PSY.RegulationDevice{V},
    U <: AbstractRegulationFormulation,
} where {V <: PSY.StaticInjection}
    add_proportional_cost!(container, AdditionalDeltaActivePowerUpVariable(), devs, U())
    add_proportional_cost!(container, AdditionalDeltaActivePowerDownVariable(), devs, U())
    return
end


text:
#=
hill_climbing_search:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-10-13

Goal is as top of mountain.
Searching to climbing direction.
Cost(distance to Goal) need to be calculated.
Evaluated by Manhattan distance.

Board as 8 puzzle
0 1 2
3 4 5
6 7 8
=#

module HCS
    # adjacent node
    const adjacent = (
        (1, 3),       # 0
        (0, 2, 4),    # 1
        (1, 5),       # 2
        (0, 4, 6),    # 3
        (1, 3, 5, 7), # 4
        (2, 4, 8),    # 5
        (3, 7),       # 6
        (4, 6, 8),    # 7
        (5, 7)        # 8
    )

    # distance to other nodes
    const distance = (
        (),
        (0, 1, 2, 1, 2, 3, 2, 3, 4), # 1
        (1, 0, 1, 2, 1, 2, 3, 2, 3), # 2
        (2, 1, 0, 3, 2, 1, 4, 3, 2), # 3
        (1, 2, 3, 0, 1, 2, 1, 2, 3), # 4
        (2, 1, 2, 1, 0, 1, 2, 1, 2), # 5
        (3, 2, 1, 2, 1, 0, 3, 2, 1), # 6
        (2, 3, 4, 1, 2, 3, 0, 1, 2), # 7
        (3, 2, 3, 2, 1, 2, 1, 0, 1)  # 8
    )

    const goal = [1, 2, 3, 4, 5, 6, 7, 8, 0]

    # distance to goal
    function get_distance(board)
        dist = 0
        for i in 1:9
            panel = board[i]
            if panel != 0
                dist += distance[panel+1][i]
            end
        end
        return dist
    end

    # hill clibming search
    function search(board, space, cost, history)
        # reached goal
        if board == goal
            for h in history
                println(h)
            end
            println("Moves count = $(length(history))")
            return true
        else
            buffer = []
            # define each pattern of next phase
            for i in adjacent[space]
                panel = board[i+1]
                b = [p for p in board]
                b[space] = b[i+1]
                b[i+1] = 0
                # add only new pattern to history
                if (b in history) == false
                    c = cost - distance[panel+1][i+1] + distance[panel+1][space]
                    # phase pattern, adjacent by space, cost
                    push!(buffer, (b, i+1, c))
                end
            end
            # select phase has minimum cost
            sort!(buffer, lt=(x, y)->x[3] < y[3])
            # buf = (phase, space, cost)
            for buf in buffer
                push!(history, buf[1])
                # recursive
                if search(buf[1], buf[2], buf[3], history) == true
                    return true
                end
                pop!(history)
            end
        end
        return false
    end

    function main()
        # initial board
        b = [8, 6, 7, 2, 5, 4, 3, 0, 1]

        @time search(b, findall(x->x == 0, b)[1], get_distance(b), [b])
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .HCS
    HCS.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
if isdir(Pkg.dir("GLMNet"))
    eval(Expr(:import, :GLMNet))
    const glmnet_loaded = true
else
    const glmnet_loaded = false
end
import NarmaxLasso
using Base.Test

## Test utility functions
@testset "Test utility functions" begin
    include("test_util.jl")
end

## Test basis functions
@testset "Test monomials" begin
    include("test_monomials.jl")
end

## Test narmax generated regressors
@testset "Test narmax generated regressors" begin
    include("test_monomials.jl")
end

## Test dynamic system simulation
@testset "Test simulate" begin
    include("test_simulate.jl")
end

## Test pathwise coordinate descent
@testset "Test pathwise coordinate descent algorithm" begin
    if glmnet_loaded
        include("test_pathwise_coordinate_descent.jl")
    end
end

## Test Narmax Lasso estimation
@testset "Test Narmax Lasso estimation" begin
    include("test_narmax_lasso.jl")
end


text:
using RvSpectML

function extract_orders(data::AbstractVector{ST}, pipeline::PipelinePlan; recalc::Bool = false, verbose::Bool = false, orders_to_use = RvSpectML.orders_to_use_default(first(data).inst) ) where { ST<:AbstractSpectra }
   if need_to(pipeline,:extract_orders) || recalc
      if verbose println("# Extracting order list timeseries from spectra.")  end
      @assert !need_to(pipeline,:read_spectra)
      order_list_timeseries = make_order_list_timeseries(data, orders_to_use=orders_to_use)
      order_list_timeseries = filter_bad_chunks(order_list_timeseries,verbose=true)
      #RvSpectML.normalize_spectra!(order_list_timeseries,data);
      set_cache!(pipeline,:extract_orders,order_list_timeseries)
      dont_need_to!(pipeline,:extract_orders)
    end
    if has_cache(pipeline,:extract_orders) return read_cache(pipeline,:extract_orders)
    else   @error("Invalid pipeline state.")          end
end

#=
default_min_chunk_length = 75
default_min_telluric_model = 1.0

function make_chunk_from_existing_chunk(parent::RvSpectML.AbstractChunkOfSpectrum, loc::NamedTuple{(:pixels,:order),Tuple{UnitRange{Int64},Int64}}) # pixels::UnitRange, order::Integer
   λ = parent.λ.parent
   flux = parent.flux.parent
   var = parent.var.parent
   return RvSpectML.ChunkOfSpectrum(λ,flux,var,loc)
end

function find_chunks_avoiding_tellurics(order_list_timeseries::AbstractChunkListTimeseries, chid::Integer; min_chunk_length::Integer = default_min_chunk_length, min_telluric_model::Real = default_min_telluric_model )
   @assert 1 <= chid <= num_chunks(order_list_timeseries)
   @assert 3 <= min_chunk_length <= 128

   order = first(order_list_timeseries.chunk_list).data[chid].flux.indices[2]
   mask_affected_by_tellurics_at_any_time = mapreduce(obsid->(.!(view(order_list_timeseries.metadata[obsid][:tellurics], order_list_timeseries.chunk_list[obsid][chid].flux.indices[1], order_list_timeseries.chunk_list[obsid][chid].flux.indices[2])
                                          .>= min_telluric_model)), +, 1:length(order_list_timeseries) )
   order_length = length(order_list_timeseries.chunk_list[1].data[chid].flux)
   start_looking_at = 1
   chunk_list = NamedTuple{(:pixels,:order),Tuple{UnitRange{Int64},Int64}}[]
   while start_looking_at < order_length
      chunk_begin = findfirst(m->m==0,view(mask_affected_by_tellurics_at_any_time,start_looking_at:order_length) )
      if isnothing(chunk_begin)    # There's nothing left unaffected by tellurics in this order.
         break
      end
      chunk_begin += start_looking_at-1 # offset since didn't start at beginning
      chunk_end = findfirst(m->m!=0,view(mask_affected_by_tellurics_at_any_time,chunk_begin:order_length) )
      if isnothing(chunk_end)       # Don't need to break up into smaller chunks further, at least not on account of tellurics
         chunk_end = order_length
         if chunk_end-chunk_begin >= min_chunk_length  # Chunk big enough to be useful
            push!(chunk_list, ( pixels= chunk_begin:chunk_end, order=order ) )
         end
         break
      else
         chunk_end += chunk_begin-1    # offset since didn't start at beginning
         start_looking_at = chunk_end  # Start looking for next chunk at first pixel flagged as having a telluric
         chunk_end -= 1                # End chunk without tellurics one pixel earlier
         println("# start_looking_at = ", start_looking_at, "  chunk_begin = ", chunk_begin, "  chunk_end = ", chunk_end)
         println(" Found chunk w/ length=",chunk_end-chunk_begin-1)
         if chunk_end-chunk_begin >= min_chunk_length  # Chunk big enough to be useful
            push!(chunk_list, ( pixels= chunk_begin:chunk_end, order=order ))
         end
      end
   end
   #return chunk_list

   return map(obsid->RvSpectML.ChunkList(map(loc->make_chunk_from_existing_chunk(order_list_timeseries.chunk_list[obsid].data[chid], loc), chunk_list)),1:length(order_list_timeseries) )
end

function merge_chunk_lists_from_multiple_orders(vclt::Vector{Vector{ACLT}}) where { ACLT<:AbstractChunkList}
   @assert length(vclt) >= 1
   num_orders = length(vclt)
   num_times = length(first(vclt))
   for t_idx in 1:num_times
      for ord_idx in 2:num_orders
         append!(vclt[1][t_idx], vclt[ord_idx][t_idx])
      end
   end
   return vclt[1]
end

function find_wavelength_ranges_with_tellurics(order_list_timeseries::AbstractChunkListTimeseries, chid::Integer; min_chunk_length::Integer = default_min_chunk_length, min_telluric_model::Real = default_min_telluric_model )
   @assert 1 <= chid <= num_chunks(order_list_timeseries)
   @assert 3 <= min_chunk_length <= 128

   order = get_order_index(first(order_list_timeseries.chunk_list).data[chid]) # .flux.indices[2]

   mask_affected_by_tellurics_at_any_time = mapreduce(obsid->(.!(view(order_list_timeseries.metadata[obsid][:tellurics], get_pixels_range(order_list_timeseries.chunk_list[obsid][chid]), get_order_index(order_list_timeseries.chunk_list[obsid][chid]) )
                                          .>= min_telluric_model)), +, 1:length(order_list_timeseries) )
   order_length = length(order_list_timeseries.chunk_list[1].data[chid].flux)
   start_looking_at = 1
   chunk_list = NamedTuple{(:pixels,:order),Tuple{UnitRange{Int64},Int64}}[]
   while start_looking_at < order_length
      chunk_begin = findfirst(m->m==0,view(mask_affected_by_tellurics_at_any_time,start_looking_at:order_length) )
      if isnothing(chunk_begin)    # There's nothing left unaffected by tellurics in this order.
         break
      end
      chunk_begin += start_looking_at-1 # offset since didn't start at beginning
      chunk_end = findfirst(m->m!=0,view(mask_affected_by_tellurics_at_any_time,chunk_begin:order_length) )
      if isnothing(chunk_end)       # Don't need to break up into smaller chunks further, at least not on account of tellurics
         chunk_end = order_length
         if chunk_end-chunk_begin >= min_chunk_length  # Chunk big enough to be useful
            push!(chunk_list, ( pixels= chunk_begin:chunk_end, order=order ) )
         end
         break
      else
         chunk_end += chunk_begin-1    # offset since didn't start at beginning
         start_looking_at = chunk_end  # Start looking for next chunk at first pixel flagged as having a telluric
         chunk_end -= 1                # End chunk without tellurics one pixel earlier
         println("# start_looking_at = ", start_looking_at, "  chunk_begin = ", chunk_begin, "  chunk_end = ", chunk_end)
         println(" Found chunk w/ length=",chunk_end-chunk_begin-1)
         if chunk_end-chunk_begin >= min_chunk_length  # Chunk big enough to be useful
            push!(chunk_list, ( pixels= chunk_begin:chunk_end, order=order ))
         end
      end
   end
   #return chunk_list

   return map(obsid->RvSpectML.ChunkList(map(loc->make_chunk_from_existing_chunk(order_list_timeseries.chunk_list[obsid].data[chid], loc), chunk_list)),1:length(order_list_timeseries) )
end
=#


text:
"""
    struct ConeMantle{T,TR,TP,D} <: AbstractCurvedSurfacePrimitive{T}

Surface primitive describing the mantle of a [`Cone`](@ref).

## Parametric types
* `T`: Precision type.
* `TR`: Type of the radius `r`.
    * `TR == T`: CylinderMantle (constant radius `r` at all `z`).
    * `TR == Tuple{T, T}`: VaryingCylinderMantle (bottom radius at `r[1]`, top radius at `r[2]`).
* `TP`: Type of the angular range `φ`.
    * `TP == Nothing`: Full 2π Cone.
    * `TP == Tuple{T, T}`: Partial Cone ranging from `φ[1]` to `φ[2]`.
* `D`: Direction in which the normal vector points (`:inwards` or `:outwards`).

## Fields
* `r::TR`: Definition of the radius of the `ConeMantle` (in m).
* `φ::TP`: Range in polar angle `φ` over which the `ConeMantle` extends (in radians).
* `hZ::T`: Half of the height of the `ConeMantle` (in m).
* `origin::CartesianPoint{T}`: Origin of the `Cone` which has this `ConeMantle` as surface.
* `rotation::SMatrix{3,3,T,9}`: Rotation matrix of the `Cone` which has this `ConeMantle` as surface.
"""
@with_kw struct ConeMantle{T,TR,TP,D} <: AbstractCurvedSurfacePrimitive{T}
    r::TR = 1
    φ::TP = nothing
    hZ::T = 1

    # axis::Line{T} = Line{T}(CartesianPoint{T}(zero(T), zero(T), -hZ), CartesianVector{T}(zero(T), zero(T), 2hZ))
    origin::CartesianPoint{T} = zero(CartesianPoint{T})
    rotation::SMatrix{3,3,T,9} = one(SMatrix{3, 3, T, 9})
end

radius_at_z(hZ::T, rBot::T, rTop::T, z::T) where {T} = iszero(hZ) ? rBot : rBot + (hZ+z)*(rTop - rBot)/(2hZ)
radius_at_z(cm::ConeMantle{T,T}, z::T) where {T} = cm.r
radius_at_z(cm::ConeMantle{T,Tuple{T,T}}, z::T) where {T} = radius_at_z(cm.hZ, cm.r[1], cm.r[2], z)

get_φ_limits(cm::ConeMantle{T,<:Any,Tuple{T,T}}) where {T} = cm.φ[1], cm.φ[2]
get_φ_limits(cm::ConeMantle{T,<:Any,Nothing}) where {T} = T(0), T(2π)

function normal(cm::ConeMantle{T,T}, pt::CartesianPoint{T}) where {T}
    pto = _transform_into_object_coordinate_system(pt, cm)
    cyl = CylindricalPoint(pto)
    return CartesianVector(_transform_into_global_coordinate_system(
            CartesianPoint(CylindricalPoint{T}(cyl.r, cyl.φ, zero(T))), cm))
end
function normal(cm::ConeMantle{T,Tuple{T,T},<:Any,:inwards}, pt::CartesianPoint{T}) where {T}
    pto = _transform_into_object_coordinate_system(pt, cm)
    cyl = CylindricalPoint(pto)
    Δr = cm.r[2] - cm.r[1]
    Δz = 2cm.hZ
    return CartesianVector(_transform_into_global_coordinate_system(
            CartesianPoint(CylindricalPoint{T}(-one(T), cyl.φ, Δr / Δz)), cm))
end
function normal(cm::ConeMantle{T,Tuple{T,T},<:Any,:outwards}, pt::CartesianPoint{T}) where {T}
    pto = _transform_into_object_coordinate_system(pt, cm)
    cyl = CylindricalPoint(pto)
    Δr = cm.r[2] - cm.r[1]
    Δz = 2cm.hZ
    return CartesianVector(_transform_into_global_coordinate_system(
            CartesianPoint(CylindricalPoint{T}( one(T), cyl.φ, -Δr / Δz)), cm))
end

function vertices(cm::ConeMantle{T}, n_arc::Int64)::Vector{CartesianPoint{T}} where {T}
    φMin, φMax = get_φ_limits(cm)
    n_arc = _get_n_points_in_arc_φ(cm, n_arc)
    φ = range(φMin, φMax, length = n_arc+1)
    rbot = radius_at_z(cm,-cm.hZ)
    rtop = radius_at_z(cm,cm.hZ)
    botcircle = [_transform_into_global_coordinate_system(CartesianPoint{T}(rbot*cos(φ), rbot*sin(φ), -cm.hZ), cm) for φ in φ]
    topcircle = [_transform_into_global_coordinate_system(CartesianPoint{T}(rtop*cos(φ), rtop*sin(φ), cm.hZ), cm) for φ in φ]
    append!(botcircle, topcircle)
end

function sample(cm::ConeMantle{T}, spacing::T)::Vector{CartesianPoint{T}} where {T}
    φMin, φMax = get_φ_limits(cm)
    Δφ = abs(φMax - φMin)
    full2π = isnothing(cm.φ)
    rbot = radius_at_z(cm,-cm.hZ)
    rtop = radius_at_z(cm,cm.hZ)
    l = hypot(2cm.hZ, rtop - rbot)
    z = range(-cm.hZ, cm.hZ, length = max(2,1+Int(ceil(l/spacing))))
    [_transform_into_global_coordinate_system(CartesianPoint{T}((radius_at_z(cm,z) .* reverse(sincos(φ)))..., z), cm) for z in z for r in [radius_at_z(cm,z)] for φ in (r == 0 ? [φMin] : range(φMin, φMax - full2π*spacing/radius_at_z(cm,z), length = max(2,1+Int(ceil(Δφ*radius_at_z(cm,z)/spacing)))))]
end

function connections(cm::ConeMantle, n_arc::Int64)::Vector{Vector{Int64}}
    n_arc = _get_n_points_in_arc_φ(cm, n_arc)
    [[i,i+1,i+n_arc+2,i+n_arc+1] for i in 1:n_arc]
end

function connections(cm::ConeMantle, n_arc::Int64, n_vert_lines::Int64)::Vector{Vector{Int64}}
    n_arc = _get_n_points_in_arc_φ(cm, n_arc)
    verts = [[i, i + n_arc + 1] for i in _get_vert_lines_range(cm,n_arc,n_vert_lines)]
    circ1 =  [[i, i + 1] for i in 1:n_arc]
    circ2 =  [[i, i + 1] for i in n_arc+2:2*n_arc+1]
    append!(verts, circ1, circ2)
end

get_label_name(::ConeMantle) = "Cone Mantle"

const FullConeMantle{T,D} = ConeMantle{T,Tuple{T,T},Nothing,D} # ugly name but works for now, should just be `ConeMantle`...
const PartialConeMantle{T,D} = ConeMantle{T,Tuple{T,T},Tuple{T,T},D}

extremum(cm::FullConeMantle{T}) where {T} = hypot(cm.hZ, max(cm.r...))

extremum(cm::PartialConeMantle{T}) where {T} = hypot(cm.hZ, max(cm.r...))

function lines(sp::FullConeMantle{T}; n = 2) where {T}
    bot_origin = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T),zero(T), -sp.hZ), sp)
    bot_ellipse = Circle{T}(r = sp.r[1], φ = sp.φ, origin = bot_origin, rotation = sp.rotation)
    top_origin = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T),zero(T), +sp.hZ), sp)
    top_ellipse = Circle{T}(r = sp.r[2], φ = sp.φ, origin = top_origin, rotation = sp.rotation)
    φs = range(T(0), step = T(2π) / n, length = n)
    edges = [ Edge{T}(_transform_into_global_coordinate_system(CartesianPoint(CylindricalPoint{T}(sp.r[1], φ, -sp.hZ)), sp),
                      _transform_into_global_coordinate_system(CartesianPoint(CylindricalPoint{T}(sp.r[2], φ, +sp.hZ)), sp)) for φ in φs ]
    bot_ellipse, top_ellipse, edges
end
function lines(sp::PartialConeMantle{T}; n = 2) where {T}
    bot_origin = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -sp.hZ), sp)
    bot_ellipse = PartialCircle{T}(r = sp.r[1], φ = sp.φ, origin = bot_origin, rotation = sp.rotation)
    top_origin = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +sp.hZ), sp)
    top_ellipse = PartialCircle{T}(r = sp.r[2], φ = sp.φ, origin = top_origin, rotation = sp.rotation)
    φs = range(sp.φ[1], stop = sp.φ[2], length = n)
    edges = [ Edge{T}(_transform_into_global_coordinate_system(CartesianPoint(CylindricalPoint{T}(sp.r[1], φ, -sp.hZ)), sp),
                      _transform_into_global_coordinate_system(CartesianPoint(CylindricalPoint{T}(sp.r[2], φ, +sp.hZ)), sp)) for φ in φs ]
    bot_ellipse, top_ellipse, edges
end

"""
    intersection(cm::ConeMantle{T,Tuple{T,T}}, l::Line{T}) where {T}

Calculates the intersections of a `Line` with a `ConeMantle`.

## Arguments
* `cm::ConeMantle{T,Tuple{T,T}}`: The `ConeMantle`.
* `l::Line{T}`: The `Line`.

!!! note
    The function will always return 2 CartesianPoint's.
    If the line just touches the mantle, the two points will be the same.
    If the line does not touch the mantle at all, the two points will have NaN's as there coordinates.
    If the line crosses the mantle only once, two points will be returned. The two points will be the same point (the intersection).
    If the line lies inside the mantle and is parallel to it. The same point will be returned which is the origin of the line.
"""
function intersection(cm::ConeMantle{T,Tuple{T,T}}, l::Line{T}) where {T}
    obj_l = _transform_into_object_coordinate_system(l, cm) # direction is not normalized

    L1 = obj_l.origin.x
    L2 = obj_l.origin.y
    L3 = obj_l.origin.z
    D1 = obj_l.direction.x
    D2 = obj_l.direction.y
    D3 = obj_l.direction.z

    hZ = cm.hZ
    R0 = cm.r[1]

    S = (cm.r[2] - cm.r[1]) / 2hZ # r slope

    f1 = D1^2 + D2^2 - D3^2*S^2
    # f1 is zero if there is only one intersection.

    λ1, λ2 = if f1 == 0 # one intersection
        term1 = -hZ^2*S^2 - 2hZ*L3*S^2 - 2*hZ*R0*S + L1^2 + L2^2 - L3^2*S^2 - 2L3*R0*S - R0^2
        term2 = L1*sqrt(D3^2*S^2 - D2^2) - D2*L2 + D3*hZ*S^2 + D3*L3*S^2 + D3*R0*S
        λ = term1 / (2term2)
        λ, λ
    else # two or no intersections
        λ = inv(f1)
        term1 = (2D1*L1 + 2D2*L2 - 2D3*hZ*S^2 - 2D3*L3*S^2 - 2D3*R0*S)^2
        term2 = -hZ^2*S^2 - 2hZ*L3*S^2 - 2hZ*R0*S + L1^2 + L2^2 - L3^2*S^2 - 2L3*R0*S - R0^2
        term3 = -D1*L1 - D2*L2 + D3*hZ*S^2 + D3*L3*S^2 + D3*R0*S
        term4 = term1 - 4*f1*term2
        # sq::T = term4 < 0 ? T(NaN) : sqrt(term4)
        sq::T = sqrt(abs(term4))

        # @show term4
        λ1 = λ * (-sq/2 + term3)
        λ2 = λ * (+sq/2 + term3)
        λ1, λ2
    end
    ints1 = obj_l.origin + λ1 * obj_l.direction
    ints2 = obj_l.origin + λ2 * obj_l.direction
    return _transform_into_global_coordinate_system(ints1, cm),
           _transform_into_global_coordinate_system(ints2, cm)
end

# """
#     intersection(cm::ConeMantle{T,T}, l::Line{T}) where {T}

# The function will always return 2 CartesianPoint's.
# If the line just touches the mantle, the two points will be the same.
# If the line does not touch the mantle at all, the two points will have NaN's as there coordinates.
# """
# function intersection(cm::ConeMantle{T,T}, l::Line{T}) where {T}
#     obj_l = _transform_into_object_coordinate_system(l, cm) # direction is not normalized

#     L1 = obj_l.origin.x
#     L2 = obj_l.origin.y
#     L3 = obj_l.origin.z
#     D1 = obj_l.direction.x
#     D2 = obj_l.direction.y
#     D3 = obj_l.direction.z

#     f1 = D1^2 + D2^2
#     λ = inv(f1) # f1 is only 0 if obj_l is parallel to the axis of the cone
#                 # (here eZ -> D1 = D2 = 0)
#                 # We assume here that this is not the case ->
#                 # We check this in choosing the sample / evaluating dimensions in `paint!`
#     hZ = cm.hZ
#     R0 = cm.r

#     term1 = (2D1*L1 + 2D2*L2)^2
#     term2 = L1^2 + L2^2 - R0^2
#     term3 = -D1*L1 - D2*L2
#     term4 = term1 - 4*f1*term2
#     sq::T = term4 < 0 ? T(NaN) : sqrt(term1 - 4*f1*term2) # if this

#     λ1 = λ * (-sq/2 + term3)
#     λ2 = λ * (+sq/2 + term3)

#     ints1 = obj_l.origin + λ1 * obj_l.direction
#     ints2 = obj_l.origin + λ2 * obj_l.direction
#     return _transform_into_global_coordinate_system(ints1, cm),
#            _transform_into_global_coordinate_system(ints2, cm)
# end

# function get_2d_grid_ticks_and_proj(cm::ConeMantle{T}, t) where {T}
#     pts = extreme_points(cm)

#     t_idx_range_x, t_idx_range_y, t_idx_range_z = get_min_max_index_ranges((pts, t))
#     ls = (length(t_idx_range_x), length(t_idx_range_y), length(t_idx_range_z))
#     ls = (
#         ls[1] == 1 ? typemax(eltype(ls)) : ls[1],
#         ls[2] == 1 ? typemax(eltype(ls)) : ls[2],
#         ls[3] == 1 ? typemax(eltype(ls)) : ls[3]
#     )
#     eX = CartesianVector{T}(one(T),zero(T),zero(T))
#     eY = CartesianVector{T}(zero(T),one(T),zero(T))
#     eZ = CartesianVector{T}(zero(T),zero(T),one(T))
#     axis = cm.rotation * CartesianVector{T}(zero(T), zero(T), one(T))

#     proj, t1, t2 = if (axis × eZ != zero(CartesianVector{T}))
#         Val{:xy}(), t_idx_range_x, t_idx_range_y
#     elseif (axis × eY != zero(CartesianVector{T}))
#         Val{:xz}(), t_idx_range_x, t_idx_range_z
#     elseif (axis × eX != zero(CartesianVector{T}))
#         Val{:yz}(), t_idx_range_y, t_idx_range_z
#     else
#         error("Sampling Error. Have to extend cases")
#         # Should never happen. This else case can be removed after testing.
#     end

#     return t1, t2, proj
# end

# function distance_to_surface(pt::AbstractCoordinatePoint{T}, c::ConeMantle{T, <:Any, Nothing, <:Any})::T where {T}
#     pcy = CylindricalPoint(pt)
#     distance_to_line(PlanarPoint{T}(pcy.r,pcy.z), LineSegment(c))
# end

# function distance_to_surface(pt::AbstractCoordinatePoint{T}, c::ConeMantle{T, <:Any, <:AbstractInterval, <:Any})::T where {T}
#     pcy = CylindricalPoint(pt)
#     φMin::T, φMax::T, _ = get_φ_limits(c)
#     if _in_φ(pcy, c.φ)
#         return distance_to_line(PlanarPoint{T}(pcy.r,pcy.z), LineSegment(c))
#     else
#         return distance_to_line(CartesianPoint(pt), LineSegment(c, _φNear(pcy.φ, φMin, φMax)))
#     end
# end


text:
#=
Part of this work is derived from astropy/photutils and astropy/astropy. The relevant derivations
are considered under a BSD 3-clause license. =#

using Statistics
using StatsBase
using Parameters

########################################################################
# Location estimators

"""
    SourceExtractorBackground()

This estimator returns the background of the input using the SourceExtractorBackground algorithm.

The background is calculated using a mode estimator of the form `(2.5 * median) - (1.5 * mean)`.

If `(mean - median) / std > 0.3` then the median is used and if `std = 0` then the mean is used.

# Examples
```jldoctest
julia> data = ones(3, 5);

julia> SourceExtractorBackground()(data)
1.0

julia> SourceExtractorBackground()(data, dims=1)
1×5 Matrix{Float64}:
 1.0  1.0  1.0  1.0  1.0
```
"""
struct SourceExtractorBackground <: LocationEstimator end

""" Utility function for SourceExtractorBackground algorithm"""
function validate_SE(background, _mean, _median, _std)
    _std ≈ 0 && return _mean
    abs(_mean - _median) / _std > 0.3 && return _median
    return background
end

function (::SourceExtractorBackground)(data; dims = :)
    _mean = mean(data, dims = dims)
    _median = median(data, dims = dims)
    _std = std(data, dims = dims)

    background = @. 2.5 * _median - 1.5 * _mean
    return validate_SE.(background, _mean, _median, _std)
end

"""
    MMMBackground(median_factor=3, mean_factor=2)

Estimate the background using a mode estimator of the form `median_factor * median - mean_factor * mean`.
This algorithm is based on the `MMMBackground` routine originally implemented in DAOPHOT. `MMMBackground` uses factors of `median_factor=3` and `mean_factor=2` by default.
This estimator assumes that contaminated sky pixel values overwhelmingly display positive departures from the true value.

# Examples
```jldoctest
julia> x = ones(3, 5);

julia> MMMBackground()(x)
1.0

julia> MMMBackground(median_factor=4, mean_factor=3)(x, dims = 1)
1×5 Matrix{Float64}:
 1.0  1.0  1.0  1.0  1.0
```

# See Also
* [`SourceExtractorBackground`](@ref)
"""
@with_kw struct MMMBackground{T <: Number} <: LocationEstimator
    median_factor::T = 3
    mean_factor::T = 2
end

(alg::MMMBackground)(data; dims = :) = alg.median_factor .* median(data, dims = dims) .- alg.mean_factor .* mean(data, dims = dims)

"""
    BiweightLocationBackground(c = 6.0, M = nothing)

Estimate the background using the robust biweight location statistic.

``\\xi_{biloc}=M + \\frac{\\sum_{|u_i|<1}{(x_i - M)(1 - u_i^2)^2}}{\\sum_{|u_i|<1}{(1-u_i^2)^2}}``

``u_i = \\frac{(x_i - M)}{c\\cdot\\text{MAD}(x)}``

Where ``\\text{MAD}(x)`` is median absolute deviation of `x`.

# Examples
```jldoctest
julia> x = ones(3,5);

julia> BiweightLocationBackground()(x)
1.0

julia> BiweightLocationBackground(c=5.5)(x; dims = 1)
1×5 Matrix{Float64}:
 1.0  1.0  1.0  1.0  1.0
```
"""
@with_kw struct BiweightLocationBackground{T <: Number} <: LocationEstimator
    c::T = 6.0
    M::Union{Nothing,T} = nothing
end

# Consider PR in StatsBase.jl for biweight statistics
function biweight_location(data::AbstractArray, c , M)
    M = M === nothing ? median(data) : M
    MAD = mad(data, normalize = false)
    MAD ≈ 0 && return M
    u = @. (data - M) / (c * MAD)

    num = den = zero(eltype(u))
    for ui in u
        abs(ui) < 1 || continue
        num += ui * (1 - ui^2)^2
        den += (1 - ui^2)^2
    end

    den ≈ 0 && return M
    return M + (c * MAD * num) / den
end

_biweight_location(data::AbstractArray, c, M, ::Colon) = biweight_location(data, c, M)
_biweight_location(data::AbstractArray, c, M, dims) = mapslices(X->biweight_location(X, c, M), data, dims = dims)

(alg::BiweightLocationBackground)(data; dims = :) = _biweight_location(data, alg.c, alg.M, dims)

########################################################################
# RMS Estimators

"""
    StdRMS()

Uses the standard deviation statistic for background RMS estimation.

# Examples
```jldoctest
julia> data = ones(3, 5);

julia> StdRMS()(data)
0.0

julia> StdRMS()(data, dims=1)
1×5 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0
```
"""
struct StdRMS <: RMSEstimator end

(::StdRMS)(data; dims = :) = std(data; corrected = false, dims = dims)

"""
    MADStdRMS()

Uses the standard median absolute deviation (MAD) statistic for background RMS estimation.

This is typically given as

``\\sigma \\approx 1.4826 \\cdot \\text{MAD}``

# Examples
```jldoctest
julia> data = ones(3, 5);

julia> MADStdRMS()(data)
0.0

julia> MADStdRMS()(data, dims=1)
1×5 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0
```
"""
struct MADStdRMS <: RMSEstimator end

_mad(data, ::Colon) = mad(data, normalize=true)
_mad(data, dims) = mapslices(x->mad(x, normalize = true), data; dims = dims)
(::MADStdRMS)(data; dims = :) = _mad(data, dims)

"""
    BiweightScaleRMS(c=9.0, M=nothing)

Uses the robust biweight scale statistic for background RMS estimation.

The biweight scale is the square root of the biweight midvariance. The biweight midvariance uses a tuning constant, `c`, and an optional initial guess of the central value `M`.

``\\zeta^2_{biscl}= \\frac{n\\sum_{|u_i|<1}{(x_i - M)^2(1 - u_i^2)^4}}{\\left[\\sum_{|u_i|<1}{(1-u_i^2)(1-5u_i^2)}\\right]^2}``

``u_i = \\frac{(x_i - M)}{c\\cdot\\text{MAD}(x)}``

Where ``\\text{MAD}(x)`` is median absolute deviation of `x`.

# Examples
```jldoctest
julia> data = ones(3, 5);

julia> BiweightScaleRMS()(data)
0.0

julia> BiweightScaleRMS(c=3.0)(data, dims=1)
1×5 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0
```
"""
@with_kw struct BiweightScaleRMS <: RMSEstimator
    c::Number = 9.0
    M::Union{Nothing,Number} = nothing
end

function biweight_scale(x::AbstractArray{T}, c, M) where T
    length(x) == 1 && return zero(T)
    M = M === nothing ? median(x) : M
    _mad = mad(x, normalize = false)
    # avoid divide by zero error
    _mad = _mad ≈ 0 ? 1 : _mad
    u = @. (x - M) / (c * _mad)

    num = den = zero(eltype(u))
        for ui in u
        abs(ui) < 1 || continue
        num += (c * _mad * ui)^2 * (1 - ui^2)^4
        den += (1 - ui^2) * (1 - 5ui^2)
    end

    return sqrt(length(x) * num) / abs(den)
end

_biweight_scale(x::AbstractArray, c, M, ::Colon) = biweight_scale(x, c, M)
_biweight_scale(x::AbstractArray, c, M, dims) = mapslices(x->biweight_scale(x, c, M), x, dims=dims)
(alg::BiweightScaleRMS)(data; dims = :) = _biweight_scale(data, alg.c, alg.M, dims)


text:
module LinearAlgebraX
using LinearAlgebra, SimplePolynomials, Mods, Permutations

# IntegerX is any sort of real or Gaussian integer
IntegerX = Union{S,Complex{S}} where {S<:Integer}

# RationalX is a Rational or Complex Rational based on integers
RationalX = Union{Rational{S},Complex{Rational{S}}} where {S<:Integer}

TypeX = Union{IntegerX,RationalX,AbstractMod}

function _recip(x::T) where {T<:IntegerX}
    return 1 // x
end
_recip(x) = inv(x)

include("row_ops.jl")
include("detx.jl")
include("cofactor_det.jl")
include("eye.jl")
include("rrefx.jl")
include("invx.jl")
include("rankx.jl")
include("nullspacex.jl")
include("char_poly.jl")
include("homogeneous.jl")
include("perm.jl")

end # module


text:
function basic_lev_multi(dataframe; sequence_id="sequence_id",sequence_event="sequence_event",sequence_order="sequence_order")

    # The lev_master Array is the key object of this funtion. consider it the spine, or memory
    # the Elements are: The Sequence, Lev matrix, Lev Value, Trace
    lev_master = Array[]
    # Elements are: Lev average, and trace average
    lev_history = Array[]

    # ok so first thing is to create a sub dataframe from the dataframe
    gd = groupby(dataframe, Symbol(sequence_id))

    # We need to initalise the master sequence
    master_df = sort!(unique(dataframe[:, [Symbol(sequence_event), Symbol(sequence_order)]]), [Symbol(sequence_order)])
    master_seq=master_df[!,Symbol(sequence_event)]
    # println("Size of master sequence is ",size(master_seq)[1])

    ga=Array[]
    for subdf in gd
        push!(ga,[subdf[!,Symbol(sequence_id)],subdf[!,Symbol(sequence_event)],subdf[!,Symbol(sequence_order)]])
    end
    # println(typeof(ga[1][1]))

    lev_master =  Array{Array}(undef,length(ga),1)
    # lev_master = Matrix{Any}(missing,length(ga),4)
    Threads.@threads  for i in eachindex(ga)
    #     Create the first cut of lev matricies
            lev_matrix = lev_create(master_seq,ga[i][2])
        # Get the lev values
            lev_value = lev_distance(lev_matrix)
        # and do the trace on the lev matrix
            lev_master[i] = [ga[i][2],lev_matrix,lev_value,lev_trace(lev_matrix)]
    end;

    # println(typeof(lev_master))
    # println(length(lev_master))
# Cool so this creates the initial lev master object
# Now we want to begin the ole iteration

# Pull out the values so we can do sum statz
lev_value = Array{Int64}(undef,length(lev_master))
trace =  Array{Array}(undef,length(lev_master))
master_seq_hist = []
push!(master_seq_hist,copy(master_seq))

for i in eachindex(lev_master);
    lev_value[i] = lev_master[i][3]
    trace[i] = lev_master[i][4]
end

# Get the history of the run and add it to the lev hist Array
push!(lev_history,Array{Union{Array,String,Float64}}([mean(lev_value)+1,mean(vcat(trace...),dims=1),trace]))
push!(lev_history,Array{Union{Array,String,Float64}}([mean(lev_value),mean(vcat(trace...),dims=1),trace]))
positions = findall(x -> x == maximum(lev_history[end][2]), lev_history[end][2])

# ok so now iterate over this thing
while (lev_history[end][1] < lev_history[end-1][1])
    # println(lev_history[end][1],"   ",lev_history[end-1][1],)

    to_remove=positions[end][2]-1
    # take positions and remove the last one from the master sequence
    # println(" removing event ",master_seq[to_remove]," position ",to_remove," --- ","Avg Lev:",lev_history[end][1],
    #     " Master Sequence Length:",size(master_seq)[1])

    deleteat!(master_seq,to_remove)
    push!(master_seq_hist,copy(master_seq))

    # init these temp arrays
    lev_value = Array{Float64}(undef,length(lev_master))
    trace =  Array{Array}(undef,length(lev_master))
    # update the lev arrays and calculate new lev value
    Threads.@threads for i in eachindex(lev_master);
        lev_master[i][2] = lev_update(lev_master[i][2],lev_master[i][1],master_seq,to_remove)
        lev_master[i][3] = lev_distance(lev_master[i][2])
        lev_master[i][4] = lev_trace(lev_master[i][2])
        # and create the stats
        # println(typeof(lev_master[i][3]))
        lev_value[i] = lev_master[i][3]
        trace[i] = lev_master[i][4]
    end

    # Get the history of the run and add it to the lev hist Array
    push!(lev_history,Array{Union{Array,String,Float64}}([mean(lev_value),mean(vcat(trace...),dims=1),trace]))
    positions = findall(x -> x == maximum(lev_history[end][2]), lev_history[end][2])

end;

if lev_history[end][1] > lev_history[end-1][1]
    # println("Going back one and ending things")
    master_seq = master_seq_hist[end-1]
end

    Lev_return = [master_seq,lev_history,lev_master]

    return Lev_return

end

# basic_lev_multi(sample_data2, sequence_id=:trace_id, sequence_event=:activities);


text:
function _bounds(
  id::Int64,
  new_state::DiseaseState,
  extents::EventExtents{T},
  obs::EventObservations{T, M},
  events::Events{T},
  network::TransmissionNetwork) where {
    T <: DiseaseStateSequence,
    M <: ILM}

  if new_state == State_E
    lowerbounds = [events.infection[id] - extents.exposure[2]]
    upperbounds = [events.infection[id] - extents.exposure[1]]
    path_to = _pathway_to(id, network, depth = 1)
    if path_to != nothing
      if path_to[2] != nothing
        parent_host = path_to[2]
        push!(lowerbounds, events.infection[parent_host])
        if (T == SEIR) && !isnan(events.removal[parent_host])
          push!(upperbounds, events.removal[parent_host])
        end
      end
    else
      error("Event inconsistent with TransmissionNetwork")
    end
  elseif new_state == State_I
    path_from = _pathway_from(id, network, depth = 1)
    if State_E ∈ T
      lowerbounds = [obs.infection[id] - extents.infection[2]
                     events.exposure[id] + extents.exposure[1]]
      upperbounds = [obs.infection[id] - extents.infection[1]
                     events.exposure[id] + extents.exposure[2]]
      if length(path_from) > 1
        child_hosts = path_from[2:end]
        append!(upperbounds, events.exposure[child_hosts])
      end
    elseif State_E ∉ T
      lowerbounds = [obs.infection[id] - extents.infection[2]]
      upperbounds = [obs.infection[id] - extents.infection[1]]
      path_to = _pathway_to(id, network, depth = 1)
      if length(path_from) > 1
        child_hosts = path_from[2:end]
        append!(upperbounds, events.infection[child_hosts])
      end
      if path_to != nothing
        if path_to[2] != nothing
          parent_host = path_to[2]
          push!(lowerbounds, events.infection[parent_host])
          if (T == SIR) && !isnan(events.removal[parent_host])
            push!(upperbounds, events.removal[parent_host])
          end
        end
      else
        error("Event inconsistent with TransmissionNetwork")
      end
    end
  elseif new_state == State_R
    path_from = _pathway_from(id, network, depth = 1)
    lowerbounds = [obs.removal[id] - extents.removal[2]
                   obs.infection[id]]
    upperbounds = [obs.removal[id] - extents.removal[1]]
    if length(path_from) > 1
      child_hosts = path_from[2:end]
      if T == SEIR
        append!(lowerbounds, events.exposure[child_hosts])
      elseif T == SIR
        append!(lowerbounds, events.infection[child_hosts])
      end
    end
  end
  @debug "Transition of $id into $new_state with bounds: \n  [max($(round.(lowerbounds, digits=3))),\n   min($(round.(upperbounds, digits=3)))]"
  lowerbound = maximum(lowerbounds)
  upperbound = minimum(upperbounds)
  if lowerbound >= upperbound
    @warn "Lower >= upper bound for transition of i=$id into $new_state" lower = lowerbounds upper = upperbounds
  end
  return lowerbound, upperbound
end

function _bounds(
  last_event::Event{T},
  extents::EventExtents{T},
  obs::EventObservations{T, M},
  events::Events{T},
  network::TransmissionNetwork) where {
    T <: DiseaseStateSequence,
    M <: ILM}

  return _bounds(last_event.individual,
                 last_event.new_state,
                 extents, obs, events, network)
end

function _bounds(
  id::Int64,
  new_state::DiseaseState,
  extents::EventExtents{T},
  obs::EventObservations{T, M},
  events::Events{T}) where {
    T <: DiseaseStateSequence,
    M<: ILM}

  if new_state == State_E
    lowerbounds = [events.infection[id] - extents.exposure[2]]
    upperbounds = [events.infection[id] - extents.exposure[1]]
  elseif new_state == State_I
    if State_E ∈ T
      lowerbounds = [obs.infection[id] - extents.infection[2]
                     events.exposure[id] + extents.exposure[1]]
      upperbounds = [obs.infection[id] - extents.infection[1]
                     events.exposure[id] + extents.exposure[2]]
    elseif State_E ∉ T
      lowerbounds = [obs.infection[id] - extents.infection[2]]
      upperbounds = [obs.infection[id] - extents.infection[1]]
    end
  elseif new_state == State_R
    lowerbounds = [obs.removal[id] - extents.removal[2]
                   obs.infection[id]]
    upperbounds = [obs.removal[id] - extents.removal[1]]
  end
  @debug "Transition of $id into $new_state with bounds: \n  [max($(round.(lowerbounds, digits=3))),\n   min($(round.(upperbounds, digits=3)))]"
  lowerbound = maximum(lowerbounds)
  upperbound = minimum(upperbounds)
  if lowerbound >= upperbound
    @warn "Lower >= upper bound for transition of i=$id into $new_state" lower = lowerbounds upper = upperbounds
  end
  return lowerbound, upperbound
end

function _bounds(
  last_event::Event{T},
  extents::EventExtents{T},
  obs::EventObservations{T, M},
  events::Events{T}) where {
    T <: DiseaseStateSequence,
    M <: ILM}

  return _bounds(last_event.individual,
                 last_event.new_state,
                 extents, obs, events)
end


text:
#@testset "TriggerJitter" begin
#    @test 1 == 1
#end


text:
#=
2D advection using higher order FV or structured or unstructured mesh
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("FVadvection2d");

useLog("FVadvection2dlog", level=3)

# Configuration setup
domain(2)
solverType(FV)

use_tri=false;
if use_tri
    timeStepper(RK4)
    # n = 30
    # mesh(QUADMESH, elsperdim=n, bids=4)

    mesh("src/examples/utriangle.msh")  # Using triangles

    add_boundary_ID(2, (x,y) -> (x >= 0.1));
    add_boundary_ID(3, (x,y) -> (y <= 0));
    add_boundary_ID(4, (x,y) -> (y >= 0.3));
    finiteVolumeOrder(2);
else
    timeStepper(RK4)
    n = 10
    mesh(QUADMESH, elsperdim=n, bids=4)
    finiteVolumeOrder(2);
end

# Variables and BCs
u = variable("u", location=CELL)
# boundary(u, 1, FLUX, "(abs(y-0.2) < 0.11) ? sin(2*pi*t)^2 : 0") # x=0
if use_tri
    boundary(u, 1, FLUX, "(abs(y-0.06) < 0.033 && sin(6*pi*t)>0) ? 1 : 0") # x=0
else
    boundary(u, 1, FLUX, "(abs(y-0.2) < 0.11 && sin(6*pi*t)>0) ? 1 : 0") # x=0
end
boundary(u, 2, NO_BC) # x=1
boundary(u, 3, NO_BC) # y=0
boundary(u, 4, NO_BC) # y=1

# Time interval and initial condition
T = 1;
timeInterval(T)
initial(u, "0")

# The flux and source terms of the conservation equation
# F and S in the following equation:
# Dt(int(u dx)) = int(S dx) - int(F.n ds)
if use_tri
    coefficient("a", ["0.1*cos(pi*x/2/0.1)","0.3*sin(pi*x/2/0.1)"], type=VECTOR) # advection velocity
else
    coefficient("a", ["cos(pi*x/2)","sin(pi*x/2)"], type=VECTOR) # advection velocity
end
coefficient("s", ["sin(pi*x)^4 * sin(pi*y)^4"]) # source
# The "upwind" function applies upwinding to the term (a.n)*u with flow velocity a.
# The optional third parameter is for tuning. Default upwind = 0, central = 1. Choose something between these.
flux(u, "upwind(a,u)")
# source(u, "0.1 * s")

# exportCode("fvad2dcode") # uncomment to export generated code to a file
# importCode("fvad2dcode") # uncomment to import code from a file

solve(u)

finalize_finch()

# output to file
# output_values(u, "fvad2d", format="csv");
#output_values(u, "fvad2d", format="vtk");

##### Uncomment below to plot

# xy = Finch.fv_info.cellCenters

# using Plots
# pyplot();
# display(plot(xy[1,:], xy[2,:], u.values[:], st=:surface))

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
module BinaryDecSSA

    """
    Need to more properly document and add ability to plot/output m intstead of N₊.
    """

    export BinDecMod, Ensemble, EnsProb, CentralMoments

    using Distributions, Parameters, StatsBase, LinearAlgebra, Colors

    """
    Define the binary decision model struct for the SSA.
    """
    @with_kw struct BinDecMod
        N::Int64 = 100
        pars::Vector{Float64}; # can be [ϵ,μ], [ϵ1,ϵ2,μ] or [ϵ1,ϵ2,μ1,μ2]
        fixedIC::Bool = true
        p₀::Float64 = 0.5 # initial prob of decided to be on rh food source if fixedIC = false.
    end

    """
    Define the struct for an ensemble of sims.
    """
    @with_kw struct Ensemble
        BD::BinDecMod
        e_size::Int64 = 1 # size of ensemble
        τₘ::Real = 1E3 # max system time
        Δτ::Real = 1.0 # time spacing of storage in trajectory
        sims = SSA(e_size, τₘ, Δτ, BD)
        @assert e_size >= 1
        @assert τₘ > 0
        @assert Δτ > 0 && Δτ < τₘ
    end

    """
    Propensities.
    """
    function props(n::Int64, BD::BinDecMod)
        @unpack pars, N = BD;
        props = zeros(2) # propensity vector.
        if length(pars)==2
            ϵ = pars[1];
            μ = pars[2];
            props[1] = (N-n)*ϵ+μ*n*(N-n)/(N-1);
            props[2] = n*ϵ+μ*n*(N-n)/(N-1);
            return props
        elseif length(pars)==3
            ϵ1 = pars[1];
            ϵ2 = pars[2];
            μ = pars[3];
            props[1] = (N-n)*ϵ1+μ*n*(N-n)/(N-1);
            props[2] = n*ϵ2+μ*n*(N-n)/(N-1);
            return props
        elseif length(pars)==4
            ϵ1 = pars[1];
            ϵ2 = pars[2];
            μ1 = pars[3];
            μ2 = pars[4];
            props[1] = (N-n)*ϵ1+μ1*n*(N-n)/(N-1);
            props[2] = n*ϵ2+μ2*n*(N-n)/(N-1);
            return props
        else
            error("length of pars must be 2, 3 or 4.")
        end
    end

    """
    Mean-field fully connected SSA: species is n. Does much faster mf SSA sims.
    """
    function SSA(ens_its::Int64, τₘ::Real, Δτ::Real, BD::BinDecMod)
        @unpack p₀, fixedIC, N = BD;
        # create the storage array for the ensembles
        ens_store = Array{Array{Real}}(undef, ens_its);
        # assign the initial condition - either fixed or binomially drawn.
        if fixedIC == true
            n = ceil(Int64, N*p₀)
        else
            n = rand(Binomial(N,p₀),1)[1]
        end
        # define the sim times for storage of n
        times = convert(Array{Float64,1},LinRange(τₘ,0.0,floor(Int,τₘ/Δτ)+1))
        # loop over for each sim in the ensemble
        for it in 1:ens_its
            # n at t=0 is either fixed or binomially drawn.
            fixedIC == true ? n = ceil(Int64, N*p₀) : n = rand(Binomial(N,p₀),1)[1]
            # initialise the state vector and sys time
            sim_times = copy(times)
            store = zeros(length(times));
            τ = 0.0;
            m = 1; # counter for updating storage
            while τ < τₘ
                ps = props(n, BD); # get propensities
                f = sum(ps);
                r1, r2 = rand(2); # collect 2 ind rand nums ∈ [0,1]
                u = (1/f)*log(1/r1); # waiting time until next reaction
                k = findfirst(x -> x>=r2*f,cumsum(ps)); # find fired reaction
                # update state vector storage
                while τ+u >= sim_times[length(sim_times)]
                    store[m] = copy(n)
                    pop!(sim_times)
                    m += 1
                    if length(sim_times) == 0
                        break
                    end
                end
                # update the state variable
                k == 1 ? n += 1 : n -= 1;
                # update the time.
                τ += u;
            end
            ens_store[it] = store ./ N; # return the value of m not N₊
        end
        return (reverse(times),ens_store)
    end

    """
    Get the central moments from the ensemble
    """
    function CentralMoments(ens::Ensemble,n::Int64)
        ns = convert(Matrix{Float64},transpose(cat(ens.sims[2]...,dims=2)))
        store = Array{Float64,1}(undef, length(ns[1,:]))
        if n == 1
            for t in 1:length(ns[1,:])
                store[t] = mean(ns[:,t])
            end
        else
            for t in 1:length(ns[1,:])
                store[t] = moment(ns[:,t],n)
            end
        end
        # need to sort out times
        @unpack Δτ = ens;
        times = LinRange(0:Δτ:length(ns[1,:])-1);
        return (times, store)
    end

    """
    Get prob distribution at all times for an ensemble
    """
    function EnsProb(ens::Ensemble, T::Int64)
        @unpack BD = ens;
        @unpack N = BD;
        ns = convert(Matrix{Float64},transpose(cat(ens.sims[2]...,dims=2)))
        mod_bins = LinRange(0.0-(1/N),1.0+(1/N),N+2);
        mid_pts = LinRange(0.0,1.0,N+1);
        bin_vals = normalize(fit(Histogram, ns[:,T], mod_bins), mode=:probability).weights;
        return (mid_pts, bin_vals.*N)
    end

end # module end

max_stars_repo_name: jamesholehouse/AntModels

max_issues_repo_name: jamesholehouse/AntModels

max_forks_repo_name: jamesholehouse/AntModels


text:
using Compat
using Compat.Test
using SemidefiniteModels
import MathProgBase
const MPB = MathProgBase.SolverInterface

function sdtest(solver::MathProgBase.AbstractMathProgSolver; duals=false, tol=1e-6)
    @testset "Testing SDModel with $solver" begin
        m = MPB.ConicModel(solver)
        MPB.loadproblem!(m, "prob.dat-s")
        MPB.optimize!(m)
        @test MPB.status(m) == :Optimal
        @test isapprox(MPB.getobjval(m), 2.75)
        @test norm(MPB.getsolution(m) - [.75, 1.]) < tol
        if duals
            @test norm(MPB.getdual(m) - [0, 0, .125, .125*sqrt(2), .125, 2/3, 0, 0, 0, 0, 0]) < tol
            @test norm(MPB.getvardual(m)) < tol
        end
    end
end


text:
module Mathbox

export bar3d, vec3d, surf
using JSON
using Patchwork

include("patchworkpatch.jl")

const mbjsdir = joinpath(dirname(Base.source_path()), "js")
const mbjs = joinpath(mbjsdir, "mathbox-bundle.js") |> readstring
const mbjlinitjs = joinpath(mbjsdir, "mbjlinit.js") |> readstring

__init__() = begin
  try
      Patchwork.load_js_runtime()
      display("text/html", script_elem(mbjs, mbjlinitjs))
  catch e
      error(e)
  end
end

script_elem(jsstrs...) =
  Elem(:script, join(jsstrs,";\n")) & Dict(:type=>"text/javascript", :charset=>"utf-8")

plotdefaults = Dict(:height=>"500px", :wigglefactor=>0.0)

"""
Returns a div and a script which sets up a scene in the div with axes and a grid
"""
scene_setup_elem(divid, fulldivid, params) = begin
  divsetupjs = """
  var divid = "$divid"
  var fulldivid = "$fulldivid"
  """
  scene_setupjs = joinpath(mbjsdir, "mbscene_axes_setup.js") |> readstring
  [Elem(:div; id=fulldivid, style=Dict(:height=>params[:height])),
   script_elem(divsetupjs, scene_setupjs)]
end

scene_from_template{T<: Real}(
    data_matrices::Vector{Matrix{T}}, template::AbstractString;
    divid=randstring(12), params=plotdefaults, reinit_scene=false) = begin
  params = merge(plotdefaults, params)
  divid = string(divid)
  data_assignments = ""
  for (i,data) in enumerate(data_matrices)
    tdata = data'
    nrows, ncols = size(tdata)
    data_assignments *= """
      var data$i = $(tdata |> vec |> json)
      var drows$i = $nrows
      var dcols$i = $ncols
    """
  end
  template_jsstr = joinpath(mbjsdir, template) |> readstring
  fulldivid = "mathbox-div-"*divid
  pwelems = scene_setup_elem(divid, fulldivid, params)
  variable_assignments = """
$data_assignments
var divid = "$divid"
var mathbox = mbinstances[divid]
var view = mbinstances["view"+divid]
mbjlparams[divid] = $(params |> json)
"""
  script_contents = """
(function mbdiv$divid(){
  $variable_assignments
  $template_jsstr
})()
"""
  push!(pwelems, script_elem(script_contents))
  Elem(:div, pwelems)
end

bar3d{T <: Real}(data::Matrix{T}; divid=randstring(12), params=plotdefaults, reinit_scene=false) =
  scene_from_template(Matrix{T}[data], "mbbar3d_template.js"; divid=divid, params=params, reinit_scene=reinit_scene)

vec3d{T <: Real}(data::Matrix{T}; bases=zeros(data), divid=randstring(12), params=plotdefaults, reinit_scene=false) =
  scene_from_template(Matrix{T}[data, bases], "mbvec3d_template.js"; divid=divid, params=params, reinit_scene=reinit_scene)

surf{T<: Real}(data::Matrix{T}; divid=randstring(12), params=plotdefaults, reinit_scene=false) =
  scene_from_template(Matrix{T}[data], "mbsurf_template.js"; divid=divid, params=params, reinit_scene=reinit_scene)

end


text:
using LinearAlgebraX, LinearAlgebra

import LinearAlgebraX: permanent, detx
import LinearAlgebra: det, svd

permanent(G::SimpleGraph) = permanent(adjacency(G))
det(G::SimpleGraph) = det(adjacency(G))
detx(G::SimpleGraph) = detx(adjacency(G))

export det, detx, permanent

maxeig(A::Matrix) = maximum(eigvals(A))

"""
    maxeig(G::SimpleGraph)

Return the maximum eigenvalue of `G`'s adjacency matrix.
"""
maxeig(G::SimpleGraph) = maxeig(adjacency(G))

export maxeig

svd(G::SimpleGraph) = svdvals(adjacency(G))
export svd


text:
using Dates
using OpenSSL
using Sockets

function test_server()
    x509_certificate = X509Certificate()

    evp_pkey = EvpPKey(rsa_generate_key())
    x509_certificate.public_key = evp_pkey

    x509_name = X509Name()
    add_entry(x509_name, "C", "US")
    add_entry(x509_name, "ST", "Isles of Redmond")
    add_entry(x509_name, "CN", "www.redmond.com")

    adjust(x509_certificate.time_not_before, Second(0))
    adjust(x509_certificate.time_not_after, Year(1))

    x509_certificate.subject_name = x509_name
    x509_certificate.issuer_name = x509_name

    sign_certificate(x509_certificate, evp_pkey)

    server_socket = listen(5000)
    accepted_socket = accept(server_socket)

    # Create and configure server SSLContext.
    ssl_ctx = OpenSSL.SSLContext(OpenSSL.TLSv12ServerMethod())
    _ = OpenSSL.ssl_set_options(ssl_ctx, OpenSSL.SSL_OP_NO_COMPRESSION)

    OpenSSL.ssl_set_ciphersuites(ssl_ctx, "TLS_AES_256_GCM_SHA384:TLS_CHACHA20_POLY1305_SHA256:TLS_AES_128_GCM_SHA256")
    OpenSSL.ssl_use_certificate(ssl_ctx, x509_certificate)
    OpenSSL.ssl_use_private_key(ssl_ctx, evp_pkey)

    ssl_stream = SSLStream(ssl_ctx, accepted_socket, accepted_socket)

    OpenSSL.accept(ssl_stream)

    bytes_available = bytesavailable(ssl_stream)
    request = read(ssl_stream, bytes_available)
    reply = "reply: $(String(request))"

    # eof(ssl_stream) will block

    # Verify the are no more bytes available in the stream.
    @test bytesavailable(ssl_stream) == 0

    write(ssl_stream, reply)

    close(ssl_stream)
    finalize(ssl_ctx)
    return nothing
end

function test_client()
    tcp_stream = connect(5000)

    ssl_ctx = OpenSSL.SSLContext(OpenSSL.TLSv12ClientMethod())
    ssl_options = OpenSSL.ssl_set_options(ssl_ctx, OpenSSL.SSL_OP_NO_COMPRESSION)

    # Create SSL stream.
    ssl_stream = SSLStream(ssl_ctx, tcp_stream, tcp_stream)

    #TODO expose connect
    OpenSSL.connect(ssl_stream)

    # Verify the server certificate.
    x509_server_cert = OpenSSL.get_peer_certificate(ssl_stream)

    @test String(x509_server_cert.issuer_name) == "/C=US/ST=Isles of Redmond/CN=www.redmond.com"
    @test String(x509_server_cert.subject_name) == "/C=US/ST=Isles of Redmond/CN=www.redmond.com"

    request_str = "GET / HTTP/1.1\r\nHost: localhost\r\nUser-Agent: curl\r\nAccept: */*\r\n\r\nRequest_body."

    written = unsafe_write(ssl_stream, pointer(request_str), length(request_str))

    @test length(request_str) == written

    response_str = String(read(ssl_stream))

    @test response_str == "reply: $request_str"

    close(ssl_stream)
    finalize(ssl_ctx)
    return nothing
end

max_stars_repo_name: grlap/OpenSSL

max_issues_repo_name: grlap/OpenSSL

max_forks_repo_name: grlap/OpenSSL


text:
function allocate_array_of_arrays(num, inner_size, inner_type,
                                  constructor = inner_type)
    arrays = Array{inner_type, 1}(num)
    for i in 1:num
        arrays[i] = constructor(inner_size)
    end
    return arrays
end

const UTM = UpperTriangular{Float64, Matrix{Float64}}
make_UTM(args...) = UpperTriangular(Matrix{Float64}(args...))
# TODO: UpperTriangular wastes almost half of memory; fix it


text:
using SafeTestsets

@safetestset "IFGF operator" begin include("ifgfoperator_test.jl") end

max_stars_repo_name: WaveProp/IFGF

max_issues_repo_name: WaveProp/IFGF

max_forks_repo_name: WaveProp/IFGF


text:
module IfElse

ifelse(args...) = Core.ifelse(args...)

end