text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAppsV1beta2StatefulSet <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiAppsV1beta2StatefulSetSpec } # spec name: spec
    status::Any # spec type: Union{ Nothing, IoK8sApiAppsV1beta2StatefulSetStatus } # spec name: status

    function IoK8sApiAppsV1beta2StatefulSet(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1beta2StatefulSet, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiAppsV1beta2StatefulSet, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiAppsV1beta2StatefulSet, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiAppsV1beta2StatefulSet, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        validate_property(IoK8sApiAppsV1beta2StatefulSet, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApiAppsV1beta2StatefulSet

const _property_map_IoK8sApiAppsV1beta2StatefulSet = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApiAppsV1beta2StatefulSet = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiAppsV1beta2StatefulSetSpec", Symbol("status")=>"IoK8sApiAppsV1beta2StatefulSetStatus")
Base.propertynames(::Type{ IoK8sApiAppsV1beta2StatefulSet }) = collect(keys(_property_map_IoK8sApiAppsV1beta2StatefulSet))
Swagger.property_type(::Type{ IoK8sApiAppsV1beta2StatefulSet }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAppsV1beta2StatefulSet[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1beta2StatefulSet }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1beta2StatefulSet[property_name]

function check_required(o::IoK8sApiAppsV1beta2StatefulSet)
    true
end

function validate_property(::Type{ IoK8sApiAppsV1beta2StatefulSet }, name::Symbol, val)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

module DrillHoles

using CSV
using DataFrames
using StatsBase: mean, weights

include("definitions.jl")
include("compositing.jl")
include("desurvey.jl")
include("mergetables.jl")
include("validations.jl")

export
  composite,
  drillhole,
  exportwarns,
  Collar,
  Interval,
  Survey
end

ext: jl


text:
module FreqTables
    using CategoricalArrays
    using Tables
    using NamedArrays

    include("freqtable.jl")

    export freqtable, proptable, prop
end # module


text:
using RCall

include("ProblemGenerator.jl")
include("VirusCallBacks.jl")

prob = ModelSetup(:Virus,:Stochastic)

sol = solve(prob,CVODE_BDF(linear_solver=:GMRES),saveat=0.1,callback=cb)

#-----------IFN Animation-----------

#Get the max value for IFNβ concentration for whole simulation
Maxval = maximum(sol[:,:,7,:])
#Loop through time points to make an animation
anim = @animate for i = tspan[1]:0.1:tspan[2]
    heatmap(sol(i)[:,:,7],clims=(0.0,Maxval),
            title="Time = " * string(i) * " hrs")
end

gif(anim,"./Figures/IFNAnimation.gif")

#-----------Cell States Plots-----------
saveTimePoints = range(tspan[1],tspan[2],step=0.1)
allCellStates = zeros(Int64,length(saveTimePoints),3)

for k=1:length(saveTimePoints)
    allCellStates[k,:] = cellStates(saveTimePoints[k],prob.p)
end

plot(saveTimePoints,allCellStates./nCells,framestyle = :box,linewidth=2)

#-----------IRNb State Plots-----------

allStateData = DataFrame()
for (i,st) in enumerate(statesNames)
    allStateData[!,Symbol(st)] = vec(sol[:,:,i,:])
end

timeLength = length(sol.t)
allStateData.Cell = repeat(1:nCells,timeLength)
allStateData.Time = repeat(sol.t,inner=nCells)
allStateData.Infected = repeat(vec(prob.u0[:,:,2] .> 0.0),timeLength)

@rput allStateData
R"""
library(ggplot2)
library(ggpubr)

stateAve = aggregate(allStateData[,1:14], list(allStateData$Infected,allStateData$Time), mean)
colnames(stateAve)[1] <- "Cell.State"
colnames(stateAve)[2] <- "Time"

logic<- unlist(lapply(stateAve$Cell.State, function(x) x == TRUE))
stateAve$Cell.State[which(logic == TRUE)] <- "Primary"
stateAve$Cell.State[which(logic == FALSE)] <- "Secondary"
low = aggregate(allStateData[,1:14], list(allStateData$Infected,allStateData$Time), FUN = 'quantile',probs=0.05)
high = aggregate(allStateData[,1:14], list(allStateData$Infected,allStateData$Time), FUN = 'quantile',probs=0.95)

commonFigureOptions <- list( scale_x_continuous(breaks=seq(0, 48, 12)),
  theme_pubr(border=TRUE),
  ylab("nM"),
  xlab("Time (hours)")
  theme(plot.title = element_text(hjust = 0.5),aspect.ratio = 1))

p1 <- ggplot(stateAve) + geom_line(aes(y=DNA, x=Time, group=Cell.State, color = Cell.State))+
  geom_ribbon(aes(ymin=low$DNA, ymax=high$DNA, x=Time,group=Cell.State,fill = Cell.State), alpha = 0.3) +
  ggtitle("DNA") +
  commonFigureOptions

p2 <- ggplot(stateAve) + geom_line(aes(y=IFNb, x=Time, group=Cell.State, color = Cell.State))+
  geom_ribbon(aes(ymin=low$IFNb, ymax=high$IFNb, x=Time,group=Cell.State,fill = Cell.State), alpha = 0.3) +
  ggtitle("IFNb") +
  commonFigureOptions

p3 <- ggplot(stateAve) + geom_line(aes(y=IRF7, x=Time, group=Cell.State, color = Cell.State))+
  geom_ribbon(aes(ymin=low$IRF7, ymax=high$IRF7, x=Time,group=Cell.State,fill = Cell.State), alpha = 0.3) +
  ggtitle("IRF7") +
  commonFigureOptions

p4 <- ggplot(stateAve) + geom_line(aes(y=TREX1, x=Time, group=Cell.State, color = Cell.State))+
  geom_ribbon(aes(ymin=low$TREX1, ymax=high$TREX1, x=Time,group=Cell.State,fill = Cell.State), alpha = 0.3) +
  ggtitle("TREX1") +
  commonFigureOptions

figure <- ggarrange(p1, p2, p3, p4,
                    labels = c("A", "B", "C", "D"),
                    common.legend = TRUE, legend = "right",
                    align = "hv",
                    ncol = 2, nrow = 2)

ggsave("./Figures/Figure2Stoch.pdf")
"""

#-----------HeatMap Cell States -----------

function CellStateGrid(θ::ParContainer,t::Float64)
  grid = zeros(Int64,size(θ.cellsInfected))
  grid[θ.cellsInfected .<= t] .= 1
  grid[θ.cellsDead .<= t] .= 2

  return grid
end

anim = @animate for i = tspan[1]:0.1:tspan[2]
    heatmap(CellStateGrid(prob.p,i),clims=(0.0,2.0),
            title="Time = " * string(i) * " hrs")
end

gif(anim,"./Figures/CellStateGrid.gif")


text:
struct CFL
    prefactor::Float64
end

function calcSingleTimeStep(cell::PrimitiveVars,
                            dx::Float64,
                            ss::Float64)::Float64
    return dx/(ss+abs(cell.velocity))
end

function calcTimeStep(tsf::CFL,
                      state::HydroSnapshot,
                      eos,
                      cached)
    dx_list = diff(state.grid)
    dt_list = [calcSingleTimeStep(cell, dx, ss)
               for (cell, dx, ss)
               in zip(state.cells, dx_list, cached["sound_speeds"])]
    return minimum(dt_list)*tsf.prefactor
end

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
# Autogenerated wrapper script for MicrosoftMPI_jll for i686-w64-mingw32
export generate_compile_time_mpi_constants, libload_time_mpi_constants, libmpi, mpiexec

JLLWrappers.@generate_wrapper_header("MicrosoftMPI")
JLLWrappers.@declare_executable_product(generate_compile_time_mpi_constants)
JLLWrappers.@declare_library_product(libload_time_mpi_constants, "libload_time_mpi_constants.dll")
JLLWrappers.@declare_library_product(libmpi, "msmpi.dll")
JLLWrappers.@declare_executable_product(mpiexec)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        generate_compile_time_mpi_constants,
        "bin\\generate_compile_time_mpi_constants.exe",
    )

    JLLWrappers.@init_library_product(
        libload_time_mpi_constants,
        "bin\\libload_time_mpi_constants.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmpi,
        "bin\\msmpi.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        mpiexec,
        "bin\\mpiexec.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct Sampler
    mVkDevice::vk.VkDevice
    mInfo::SamplerCreateInfo
    mHandle::vk.VkSampler

    function Sampler(vkDevice::vk.VkDevice, info::SamplerCreateInfo)
        handleRef = Ref{vk.VkSampler}()
        if vk.vkCreateSampler(vkDevice, info.handleRef(), C_NULL, handleRef) != vk.VK_SUCCESS
            error("Failed to create sampler.")
        end
        this = new(vkDevice, info, handleRef[])
        return this
    end
end
@class Sampler

function destroy(this::Sampler)
    vk.vkDestroySampler(this.mVkDevice, this.mHandle, C_NULL)
end

function handle(this::Sampler)::vk.VkSampler
    return this.mHandle
end


text: using Stheno: ConstantMean, ZeroMean, ConstantKernel, ZeroKernel, pairwise, AV, AM, GPC


text:
//Test 3 in Julia

function a()
	x = 1
		if ~= x 1 then
			print(0)
		else
			print(1)
		end
end

max_stars_repo_name: tylerholmes/Interpreter

max_issues_repo_name: tylerholmes/Interpreter

max_forks_repo_name: tylerholmes/Interpreter


text:
using Plots

function create_plot(x,y)
	plt = plot([x],[y],title="Error function",label=["Error"])
	sleep(0.005)
	return plt
end

function update_plot(plt,x,y)
    push!(plt, 1, x, 1.0*y)
    display(plt)
    sleep(0.005)
end


text:
module Util

import Base.==

mutable struct BinMaxHeap{TD<:AbstractVector{Int}, TI<:Integer}
    data::TD
    ind::TI
end

BinMaxHeap(cap::T) where {T <: Integer} = BinMaxHeap(Vector{Int}(undef, cap), 0)

==(bh1::BinMaxHeap, bh2::BinMaxHeap) = bh1.data[1:bh1.ind] == bh2.data[1:bh2.ind]

Base.length(bh::BinMaxHeap) = bh.ind

Base.getindex(bh::BinMaxHeap, ind::Int) = bh.data[ind]

function top(bh::BinMaxHeap)
    if bh.ind == 0
        throw(BoundsError(
            "attempt to access top of 0-element BinMaxheap"
            )
        )
    end
    return bh.data[1]
end

# BEGIN: From DataStructures.jl/src/heaps/binary_heap.jl

function push!(bh::BinMaxHeap, r::Int)
    bh.data[bh.ind+1] = r
    bh.ind += 1
    _heap_bubble_up!(bh)
end

function pop!(bh::BinMaxHeap)
    max_r = bh.data[1]
    bh.data[1] = bh.data[bh.ind]
    bh.ind -= 1
    if bh.ind != 0
        _heap_bubble_down!(bh)
    end
    return max_r
end

function replace_least!(bh::BinMaxHeap, r::Int)
    pop!(bh)
    push!(bh, r)
end

function _heap_bubble_up!(bh::BinMaxHeap)
    i::Int = bh.ind
    @inbounds r = bh.data[i]

    while i > 1  # nd is not root
        p = i >> 1
        @inbounds rp = bh.data[p]

        if r > rp
            # move parent downward
            @inbounds bh.data[i] = rp
            i = p
        else
            break
        end
    end

    if i != bh.ind
        @inbounds bh.data[i] = r
    end
end

function _heap_bubble_down!(bh::BinMaxHeap)
    i = 1
    @inbounds r = bh.data[i]
    swapped = true
    n = length(bh)
    last_parent = n >> 1

    while swapped && i <= last_parent
        lc = i << 1
        if lc < n   # contains both left and right children
            rc = lc + 1
            @inbounds lr = bh.data[lc]
            @inbounds rr = bh.data[rc]
            if rr > lr #compare(comp, rv, lv)
                if rr > r  #compare(comp, rv, v)
                    @inbounds bh.data[i] = rr
                    i = rc
                else
                    swapped = false
                end
            else
                if lr > r
                    @inbounds bh.data[i] = lr
                    i = lc
                else
                    swapped = false
                end
            end
        else        # contains only left child
            @inbounds lr = bh.data[lc]
            if lr > r
                @inbounds bh.data[i] = lr
                i = lc
            else
                swapped = false
            end
        end
    end

    bh.data[i] = r
end

# END: From DataStructures.jl/src/heaps/binary_heap.jl

end  # module


text:
using OrdinaryDiffEq, Test, LinearAlgebra, Statistics

@testset "Mass Matrix Accuracy Tests" begin

  # create mass matrix problems
  function make_mm_probs(mm_A, ::Type{Val{iip}}) where iip
    # iip
    mm_b = vec(sum(mm_A; dims=2))
    function mm_f(du,u,p,t)
      mul!(du,mm_A,u)
      du .+= t * mm_b
      nothing
    end
    mm_g(du,u,p,t) = (@. du = u + t; nothing)

    # oop
    mm_f(u,p,t) = mm_A * (u .+ t)
    mm_g(u,p,t) = u .+ t

    mm_analytic(u0, p, t) = @. 2 * u0 * exp(t) - t - 1

    u0 = ones(3)
    tspan = (0.0, 1.0)

    prob = ODEProblem(ODEFunction{iip,true}(mm_f, analytic=mm_analytic, mass_matrix=mm_A), u0, tspan)
    prob2 = ODEProblem(ODEFunction{iip,true}(mm_g, analytic=mm_analytic), u0, tspan)

    prob, prob2
  end

  # test each method for exactness
  for iip in (false, true)
    mm_A = Float64[-2 1 4; 4 -2 1; 2 1 3]
    prob, prob2 = make_mm_probs(mm_A, Val{iip})

    sol = solve(prob,  ImplicitEuler(),dt=1/10,adaptive=false)
    sol2 = solve(prob2,ImplicitEuler(),dt=1/10,adaptive=false)

    @test norm(sol .- sol2) ≈ 0 atol=1e-7

    sol = solve(prob,  RadauIIA5(),dt=1/10,adaptive=false)
    sol2 = solve(prob2,RadauIIA5(),dt=1/10,adaptive=false)

    @test norm(sol .- sol2) ≈ 0 atol=1e-12

    sol = solve(prob,  ImplicitMidpoint(extrapolant = :constant),dt=1/10)
    sol2 = solve(prob2,ImplicitMidpoint(extrapolant = :constant),dt=1/10)

    if iip
      sol = solve(prob  ,Rosenbrock23())
      sol2 = solve(prob2,Rosenbrock23())

      @test norm(sol .- sol2) ≈ 0 atol=1e-11

      sol = solve(prob, Rosenbrock32())
      sol2 = solve(prob2,Rosenbrock32())

      @test norm(sol .- sol2) ≈ 0 atol=1e-11

      sol = solve(prob,  ROS3P())
      sol2 = solve(prob2,ROS3P())

      @test norm(sol .- sol2) ≈ 0 atol=1e-11

      sol = solve(prob,  Rodas3())
      sol2 = solve(prob2,Rodas3())

      @test norm(sol .- sol2) ≈ 0 atol=1e-11

      sol = solve(prob,  RosShamp4())
      sol2 = solve(prob2,RosShamp4())

      @test norm(sol .- sol2) ≈ 0 atol=1e-10

      sol = solve(prob,  Rodas4())
      sol2 = solve(prob2,Rodas4())

      @test norm(sol .- sol2) ≈ 0 atol=1e-9

      sol = solve(prob,  Rodas5())
      sol2 = solve(prob2,Rodas5())

      @test norm(sol .- sol2) ≈ 0 atol=1e-7
    end
  end

  # test functional iteration
  for iip in (false, true)
    prob, prob2 = make_mm_probs(Matrix{Float64}(1.01I, 3, 3), Val{iip})

    sol = solve(prob,ImplicitEuler(
                          nlsolve=NLFunctional(tol=1e-7)),dt=1/10,adaptive=false)
    sol2 = solve(prob2,ImplicitEuler(nlsolve=NLFunctional(tol=1e-7)),dt=1/10,adaptive=false)
    @test norm(sol .- sol2) ≈ 0 atol=1e-7

    sol = solve(prob, ImplicitMidpoint(extrapolant = :constant,
                          nlsolve=NLFunctional(tol=1e-7)),dt=1/10)
    sol2 = solve(prob2,ImplicitMidpoint(extrapolant = :constant, nlsolve=NLFunctional(tol=1e-7)),dt=1/10)
    @test norm(sol .- sol2) ≈ 0 atol=1e-7

    sol = solve(prob,ImplicitEuler(nlsolve=NLAnderson()),dt=1/10,adaptive=false)
    sol2 = solve(prob2,ImplicitEuler(nlsolve=NLAnderson()),dt=1/10,adaptive=false)
    @test norm(sol .- sol2) ≈ 0 atol=1e-7
    @test norm(sol[end] .- sol2[end]) ≈ 0 atol=1e-7

    sol = solve(prob, ImplicitMidpoint(extrapolant = :constant, nlsolve=NLAnderson(tol=1e-6)),dt=1/10)
    sol2 = solve(prob2,ImplicitMidpoint(extrapolant = :constant, nlsolve=NLAnderson(tol=1e-6)),dt=1/10)
    @test norm(sol .- sol2) ≈ 0 atol=1e-7
    @test norm(sol[end] .- sol2[end]) ≈ 0 atol=1e-7
  end
end

# Singular mass matrices

@testset "Mass Matrix Tests with Singular Mass Matrices" begin
  function f!(du, u, p, t)
      du[1] = u[2]
      du[2] = u[2] - 1.
      return
  end

  u0 = [0.,1.]
  tspan = (0.0, 1.0)

  M = fill(0., 2,2)
  M[1,1] = 1.

  m_ode_prob = ODEProblem(ODEFunction(f!;mass_matrix=M), u0, tspan)
  @test_nowarn sol = solve(m_ode_prob, Rosenbrock23())

  M = [0.637947  0.637947
       0.637947  0.637947]

  inv(M) # not caught as singular

  function f2!(du, u, p, t)
      du[1] = u[2]
      du[2] = u[1]
      return
  end
  u0 = zeros(2)

  m_ode_prob = ODEProblem(ODEFunction(f2!;mass_matrix=M), u0, tspan)
  @test_nowarn sol = solve(m_ode_prob, Rosenbrock23())
end


text:
function properties_solve(n)

  reps = parse(Int64, ENV["LAMP_REPS"])

  B = randn(n, 200)

  # SPD
  A = randn(n,n)
  A = Matrix(Symmetric(A, :L)) + I*n
  @test issymmetric(A)
  @test isposdef(A)
  Benchmarker.add_data(csv, "solve_spd", Benchmarker.measure(reps, solve, A, B))
  @test isposdef(A) && issymmetric(A)

  # Symmetric
  A = randn(n,n)
  A = Matrix(Symmetric(A, :L))
  A[1, 1] = -1.0
  Benchmarker.add_data(csv, "solve_sym", Benchmarker.measure(reps, solve, A, B))
  @test issymmetric(A) && !isposdef(A)

  # Triangular
  A = randn(n,n)
  A = Matrix(LowerTriangular(A))
  Benchmarker.add_data(csv, "solve_tri", Benchmarker.measure(reps, solve, A, B))
  @test istril(A)

  # Diagonal
  A = randn(n,n)
  A = Matrix(Diagonal(A))
  Benchmarker.add_data(csv, "solve_dia", Benchmarker.measure(reps, solve, A, B))
  @test isdiag(A)
end

function solve(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2})
  result = zeros(size(B))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result = A\B

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end


text:
# ODE23S  Solve stiff systems based on a modified Rosenbrock triple
# (also used by MATLAB's ODE23s); see Sec. 4.1 in
#
# [SR97] L.F. Shampine and M.W. Reichelt: "The MATLAB ODE Suite," SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22
#
# supports keywords: points = :all | :specified (using dense output)
#                    jacobian = G(t,y)::Function | nothing (FD)

# Internal

immutable ModifiedRosenbrockStepper{T<:Number} <: AbstractStepper
    d :: T
    e32 :: T

    function ModifiedRosenbrockStepper()
        d   = T(1/(2 + sqrt(2)))
        e32 = T(6 + sqrt(2))
        new(d,e32)
    end
end

# TODO: is this correct?
order(ModifiedRosenbrockStepper) = 2

# define the set of ODE problems with which this stepper can work
solve(ode :: ExplicitODE, stepper :: ModifiedRosenbrockStepper, options :: Options) =
    Solver{ModifiedRosenbrockStepper}(ode,stepper,options)

# lower level interface (iterator)

"""
The state for the Rosenbrock stepper

- step:  Last successful step
- F1,F2: Work arrays for storing the intermediate values of y'
- J:     Jacobian
- iters: Number of successful steps made

"""
type RosenbrockState{T,S} <: AbstractState
    step ::Step{T,S}
    dt   ::T
    F1   ::S
    F2   ::S
    J # :: ?
    iters::Int
end

# for debugging
function show(io::IO, state :: RosenbrockState)
    show(io,state.step)
    println("dt =$(state.dt)")
    println("F1 =$(state.F1)")
    println("F2 =$(state.F2)")
    println("J  =$(state.J)")
end

function start(s :: Solver{ModifiedRosenbrockStepper})
    t  = s.ode.t0
    dt = s.options.initstep
    y  = s.ode.y0
    dy = zero(y)

    J  = Array(eltype(y),length(y),length(y))

    step  = Step(t,deepcopy(y),deepcopy(dy))
    state = RosenbrockState(step,
                            dt,
                            zero(y), # F1
                            zero(y), # F2
                            J,       # J
                            0)       # iters
    # initialize the derivative and the Jacobian
    s.ode.F!(t,y,step.dy)
    s.ode.jac!(t,y,state.J)

    return state
end

function next(s     :: Solver{ModifiedRosenbrockStepper},
              state :: RosenbrockState)

    stepper = s.stepper
    ode     = s.ode
    step    = state.step
    opts    = s.options

    F1, F2, J = state.F1, state.F2, state.J

    t, dt, y, dy = step.t, state.dt, step.y, step.dy
    # F!, jac! = ode.F!, ode.jac!
    d, e32 = stepper.d, stepper.e32

    F0 = dy

    while true

        state.iters += 1
        if state.iters > s.options.maxiters
            return ((step.t,step.y), state)
        end

        # trim the step size to match the bounds of integration
        dt = min(s.options.tstop-t,dt)

        # TODO: this should go to a specialized function for type stabilty sake
        # maybe make W a part of ExplicitODE?  Same for tder below?
        if size(J,1) == 1
            W = one(J) - dt*d*J
        else
            # note: if there is a mass matrix M on the lhs of the ODE, i.e.,
            #   M * dy/dt = F(t,y)
            # we can simply replace eye(J) by M in the following expression
            # (see Sec. 5 in [SR97])

            W = lufact( eye(J) - dt*d*J )
        end

        # Approximate time-derivative of F, we are using F1 as a
        # temporary array
        ode.F!(t+dt/100,y,F1)
        tder = 100*d*(F1-F0)

        # modified Rosenbrock formula
        # TODO: allocate some temporary space for these variables
        k1 = W \ (F0 + tder)
        ode.F!(t+dt/2, y+dt*k1/2, F1)
        k2 = W \ (F1 - k1) + k1
        ynew = y + dt*k2
        ode.F!(t+dt,   ynew,      F2)
        k3 = W \ (F2 - e32*(k2 - F1) - 2*(k1 - F0) + tder )

        delta = max(opts.reltol*max(opts.norm(y),
                                    opts.norm(ynew)),
                    opts.abstol) # allowable error

        err = (dt/6)*opts.norm(k1 - 2*k2 + k3)/delta # error estimate

        # upon a failed step decrease the step size
        dtnew = min(opts.maxstep,
                    dt*0.8*err^(-1/3) )

        # check if the new solution is acceptable
        if err <= 1

            # update the state and return
            step.t     = t+dt
            state.dt   = dtnew
            step.y[:]  = ynew
            step.dy[:] = F2
            ode.jac!(step.t,step.y,J)

            return ((step.t,step.y), state)
        else
            # continue with the decreased time step
            dt = dtnew
        end

    end

    return tout, yout

end


text:
# A collection of date and time functions.
#
# Some of these have/should be changed to use the built-in Dates module.
#
# The functions which use Julian Day Numbers are valid only for positive values,
# i.e., for dates after -4712 (4713BC).
#
# Unless otherwise specified, Julian Day Numbers are fractional values.
#
# Numeric years use the astronomical convention of a year 0: 0 = 1BC, -1 = 2BC,
# etc.
#
# Numeric months are 1-based: Jan = 1 Dec = 12.
#
# Numeric days are the same as the calendar value.
#
# Reference: Jean Meeus, _Astronomical Algorithms_, second edition, 1998, Willmann-Bell, Inc.

export
    is_leap_year,
    cal_to_jd,
    date_to_jd,
    jd_to_date,
    jd_to_cal,
    cal_to_jd_m,
    cal_to_day_of_year,
    day_of_year_to_cal,
    pesach,
    jewish_new_year,
    moslem_to_christian,
    christian_to_moslem,
    easter,
    jd_to_day_of_week,
    is_dst,
    jd_to_jcent,
    lt_to_str,
    ut_to_lt

"""
    is_leap_year(yr, gregorian=true)

Parameters:
yr        : year
gregorian : If True, use Gregorian calendar, else use Julian calendar (default: True)

Return:
True is this is a leap year, else False.

Return True if this is a leap year in the Julian or Gregorian calendars
"""
function is_leap_year(yr, gregorian=true)
    yr =  floor(Integer, yr)
    if gregorian
        return (yr % 4 == 0) && ((yr % 100 != 0) || (yr % 400 == 0))
    else
        return yr % 4 == 0
    end
end

"""
    cal_to_jd(yr, mo, day, gregorian=true)

Convert a date in the Julian or Gregorian calendars to the Julian Day Number

    Parameters:
        yr        : year
        mo        : month
        day       : day, may be fractional
        gregorian : if True, use Gregorian calendar, else use Julian calendar (default: True)

    Return:
        jd        : julian day number

presumably returns midday time, since Julian days start at noon?

"""
function cal_to_jd(yr, mo, day, gregorian=true)
    a =  floor(Integer, (14 - mo)/12)

    y = yr + 4800 - a
    m = mo + 12 * a - 3

    jdn = day +  floor(Integer, ((153 * m + 2) / 5)) + 365 * y +  floor(Integer, (y / 4)) - 32083.5
    if !gregorian
    	return jdn
    else
    	return jdn -  floor(Integer, (y/100)) +  floor(Integer, y/400) + 38
    end
end

"""
    date_to_jd(yr, mo, d, h, m, s, gregorian=true)

be careful as they don't take into account the current Daylight Savings Time.
"""
function date_to_jd(yr, mo, d, h, m, s, gregorian=true)
   return(cal_to_jd(yr, mo, d, gregorian) + hms_to_fday(h,m,s))
end

"""
    jd_to_date(jd, gregorian=true)

be careful as they don't take into account the current Daylight Savings Time.
"""
function jd_to_date(jd, gregorian=true)
    (y, month, df) = jd_to_cal(jd, gregorian)
    (f, d) = modf(df)
    (h, m, s) = fday_to_hms(f)
    return (y, month, floor(Integer, d), h, m, s)
end

"""
    jd_to_cal(jd, gregorian=true)

Convert a Julian day number to a date in the Julian or Gregorian calendars.

Parameters:
    jd        : Julian Day number
    gregorian : if True, use Gregorian calendar, else use Julian calendar (default: True)

Return:
    year
    month
    day (may be fractional)

Return a tuple (year, month, day).
"""
function jd_to_cal(jd, gregorian=true)
    (F, Z) = modf(jd + 0.5)
    if gregorian
        alpha =  floor(Integer, ((Z - 1867216.25) / 36524.25))
        A = Z + 1 + alpha -  floor(Integer, alpha / 4)
    else
        A = Z
    end

    B = A + 1524
    C =  floor(Integer, (B - 122.1) / 365.25)
    D =  floor(Integer, 365.25 * C)
    E =  floor(Integer, (B - D) / 30.6001)
    day = B - D -  floor(Integer, 30.6001 * E) + F
    if E < 14
        mo = E - 1
    else
        mo = E - 13
    end

    if mo > 2
        yr = C - 4716
    else
        yr = C - 4715
    end

    return (yr, mo, day)
end

"""
    cal_to_jd_m(yr, mo, day, gregorian=true)

Meeus' implementation of cal_to_jd - which gives the same results
"""
function cal_to_jd_m(yr, mo, day, gregorian=true)
    if mo <= 2
        yr = yr-1
        mo = mo+12
    end
    if gregorian
        A = floor(Integer, yr / 100)
        B = 2 - A + floor(Integer, A / 4)
    else
        B = 0
    end
    return  floor(Integer, 365.25 * (yr + 4716)) +  floor(Integer, 30.6001 * (mo + 1)) + day + B - 1524.5
end

"""
    cal_to_day_of_year(yr, mo, dy, gregorian=true)

Convert a date in the Julian or Gregorian calendars to day of the year (Meeus 7.1).

Parameters:
    yr        : year
    mo        : month
    day       : day
    gregorian : if True, use Gregorian calendar, else use Julian calendar (default: True)

Return:
    day number : 1 = Jan 1...365 (or 366 for leap years) = Dec 31.
"""
function cal_to_day_of_year(yr, mo, dy, gregorian=true)
    if gregorian == false
        gregorian = true
    end
    if is_leap_year(yr, gregorian)
        K = 1
    else
        K = 2
    end
    dy =  floor(Integer, dy)
    return  floor(Integer, 275 * mo / 9.0) - (K *  floor(Integer, (mo + 9) / 12.0)) + dy - 30
end

"""
    day_of_year_to_cal(yr, N, gregorian=true)

Convert a day of year number to a month and day in the Julian or Gregorian calendars.

Parameters:
    yr        : year
    N        : day of year, 1..365 (or 366 for leap years)
    gregorian : if True, use Gregorian calendar, else use Julian calendar (default: True)

Return:
    month
    day
"""
function day_of_year_to_cal(yr, N, gregorian=true)
    if is_leap_year(yr, gregorian)
        K = 1
    else
        K = 2
    end
    if (N < 32)
        mo = 1
    else
        mo =  floor(Integer, 9 * (K+N) / 275.0 + 0.98)
    end
    dy =  floor(Integer, N -  floor(Integer, 275 * mo / 9.0) + K *  floor(Integer, (mo + 9) / 12.0) + 30)
    return (mo, dy)
end

"""
    pesach(yr, gregorian=true)

Return the date of Jewish Pesach  for a year in the Julian or Gregorian calendars.

Parameters:
    yr        : year (in the Jewish calendar)
    gregorian : if True, use Gregorian calendar, else use Julian calendar (default: True)

Return:
    month
    day

Meeus - chapter 9
"""
function pesach(yr, gregorian=true)
    C =  floor(Integer, yr/100)
    S = 0
    if gregorian
        S =  floor(Integer, (3*C-5)/4)
    end
    A = yr - 3760
    a = (12 * yr + 12) % 19
    b = yr % 4
    Q = -1.904412361576 + 1.554241796621 * a + 0.25 * b -0.003177794022 * yr + S
    j = ( floor(Integer, Q) + 3 * yr + 5 * b + 2 - S) % 7
    r = Q- floor(Integer, Q)
    if j==2 || j==4 || j==6
        D =  floor(Integer, Q) + 23
    elseif j==1 && a > 6 && r >= 0.63287037
        D =  floor(Integer, Q)+24
    elseif j==0 && a > 11 && r >= 0.897723765
        D =  floor(Integer, Q) + 23
    else
        D =  floor(Integer, Q) + 22
    end

    if D > 31
        return (4, D-31)
    else
        return (3, D)
    end
end

"""
    jewish_new_year(yr)

"""
function jewish_new_year(yr)
    gr_yr = yr-1-3760
    (m, d) = pesach(yr-1)
    jd = cal_to_jd(gr_yr, m, d) + 316
    return jd_to_cal(Float64)
end

"""
    moslem_to_christian(h, m, d)

Converts a date in the moslem calendar to a date in the Christian calendar
Meeus - chapter 9

Input: h - moslem yearhe christian year
       m - moslem month number
       d - moslem day
Returns: yr, month, day in the christian calendar

This function will return meaningless results if called with a date earlier than
622, July 16th
"""
function moslem_to_christian(h, m, d)
    N = d +  floor(Integer, 29.5001*(m-1)+0.99)
    Q  =  floor(Integer, h/30)
    R  = h%30
    A =  floor(Integer, (11 * R + 3)/30)
    W = 404*Q + 354*R+208 + A
    Q1 =  floor(Integer, W/1461)
    Q2 = W % 1461
    G  = 621 + 4 *  floor(Integer, 7 * Q + Q1)
    K =  floor(Integer, Q2/365.2422)
    E =  floor(Integer, K*365.2422)
    J = Q2-E+N-1
    X = G+K

    if J > 366 && X%4 == 0
        X = X+1
        J = J - 366
    elseif J > 365 && X%4 > 0
        X = X+1
        J = J - 365
    end

    # X is the year, J is the day in the julian calendar
    # Convert to MJD
    jd =  floor(Integer, 365.25*(X-1)) + 1721423 + J
    alpha =  floor(Integer, (jd-1867216.25)/36524.25)
    if jd < 2299161
        # Before 1582, Oct 15th
        bbeta = jd
    else
        bbeta = jd+ 1 + alpha -  floor(Integer, (alpha/4))
    end
    b = bbeta+1524
    c =  floor(Integer, (b-122.1)/365.25)
    d =  floor(Integer, 365.25*c)
    e1 =  floor(Integer, (b-d)/30.6001)
    D = b-d- floor(Integer, 30.6001*e1)
    if e1 < 14
        M = e1 - 1 else M = e1 - 13
    end
    if M > 2
        X = c - 4716 else X = c - 4715
    end
    return (X, M, D)
end

"""
    christian_to_moslem(X, M, D, gregorian=true)

"""
function christian_to_moslem(X, M, D, gregorian=true)
    # first convert the gregorian date to a julian date
    if gregorian
        if M < 3
            X = X-1
            M = M+12
        end
        alpha =  floor(Integer, X/100)
        bbeta  = 2-alpha+ floor(Integer, alpha/4)
        b =  floor(Integer, 365.25 * X) +  floor(Integer, 30.6001*(M+1)) + D + 1722519 + bbeta
        c =  floor(Integer, (b-122.1)/365.25)
        d =  floor(Integer, 365.25 * c)
        e1 =  floor(Integer, (b - d)/30.6001)
        D = b-d- floor(Integer, 30.6001 * e1)
        if e1 < 14
            M = e1 - 1 else M = e1 - 13
        end
        if M > 2
            X = c-4716
        else
            X = c-4715
        end
    end
    if X % 4 == 0
        W = 1
    else
        W = 2
    end
    N =  floor(Integer, 275 * M / 9) - W *  floor(Integer, (M + 9) / 12) + D - 30
    A = X - 623
    B =  floor(Integer, A/4)
    C = A % 4
    C1 = 365.2501 * C
    C2 =  floor(Integer, C1)
    if C1 - C2 > 0.5
        C2 = C2 + 1
    end
    Dp = 1461 * B + 170 + C2
    Q =  floor(Integer, Dp/10631)
    R = Dp % 10631
    J =  floor(Integer, R/354)
    K = R % 354
    O =  floor(Integer, (11 * J + 14)/30)
    H = 30 * Q + J + 1
    JJ = K - O + N - 1

    if JJ > 354
        CL = H % 30
        DL = (11 * CL + 3) % 30
        if DL < 19
           JJ = JJ - 354
           H=H + 1
        else
           JJ = JJ - 355
           H=H + 1
        end
        if JJ == 0
           JJ = 355
           H = H - 1
        end
    end
    S =  floor(Integer, (JJ - 1)/29.5)
    m = 1 + S
    d =  floor(Integer, JJ - 29.5 * S)
    if JJ == 355
        m = 12
        d = 30
    end
    return (H, m, d)
end

"""
    easter(yr, gregorian=true)

Return the date of Western ecclesiastical Easter for a year in the Julian or Gregorian calendars.

Parameters:
    yr        : year
    gregorian : if True, use Gregorian calendar, else use Julian calendar (default: True)

Return:
    month
    day

Julian formula found in Meeus.

Gregorian formula taken from Wikipedia (for the formula published in Nature, April 20th, 1876)
"""
function easter(yr, gregorian=true)
    yr =  floor(Integer, yr)
    if gregorian
        a = yr % 19
        b =  floor(Integer, yr / 100)
        c = yr % 100
        d =  floor(Integer, b / 4)
        e1 = b % 4
        f =  floor(Integer, (b + 8) / 25)
        g =  floor(Integer, (b - f + 1) / 3)
        h = (19 * a + b - d - g + 15) % 30
        i =  floor(Integer, c / 4)
        k = c % 4
        l = (32 + 2 * e1 + 2 * i - h - k) % 7
        m =  floor(Integer, (a + 11 * h + 22 * l) / 451)
        tmp = h + l - 7 * m + 114
    else
        a = yr % 4
        b = yr % 7
        c = yr % 19
        d = (19 * c + 15) % 30
        e1 = (2 * a + 4 * b - d + 34) % 7
        tmp = d + e1 + 114
    end
    mo =  floor(Integer, tmp / 31)
    dy = (tmp % 31) + 1
    return (mo, dy)
end

"""
    jd_to_day_of_week(jd::Float64)

Return the day of week for a Julian Day Number.

The Julian Day Number must be for 0h UT.

Parameters:
    jd : Julian Day number

Return:
    day of week: 0 = Sunday...6 = Saturday.

"""
function jd_to_day_of_week(jd)
    i = jd + 1.5
    return  floor(Integer, i) % 7
end

"""
    is_dst(jd)

Is this instant within the Daylight Savings Time period

# for Europe: last sunday of march 1:00UTC to last sunday of october 1:00UTC

if daylight_timezone_name is false, the function always returns false.

Parameters:
    jd : Julian Day number representing an instant in Universal Time

Return:
    True if Daylight Savings Time is in effect, False otherwise.
"""
function is_dst(jd::Float64)
    global daylight_timezone_name
    if daylight_timezone_name == ""
        return false
    end
    #  What year is this?
    (yr, mon, day) = jd_to_cal(jd)
    # Last day of march
    starting = cal_to_jd(yr, 3, 31)
    # Back to last sunday
    dow = jd_to_day_of_week(starting)
    if dow >=0
        starting = starting - dow
    end
    # Advance to 1AM
    starting = starting + 1.0 / 24
    # Before the beginning of the period ?
    if jd < starting
        return false
    end
    # Last day in October
    stop = cal_to_jd(yr, 10, 31)
    # Backup to the last Sunday
    dow = jd_to_day_of_week(stop)
    stop = stop - dow
    #  Advance to 1AM
    stop = stop + 1.0 / 24
    # Before the end of the period ?
    if jd < stop
        return true
    end
    # After the end of the period
    return false
end

"""
    jd_to_jcent(jd)

Return the number of Julian centuries since J2000.0

Parameters:
    jd : Julian Day number

Return:
    Julian centuries
"""
function jd_to_jcent(jd)
    return (jd - 2451545.0) / 36525.0
end

"""
    lt_to_str(jd, zone="", level="second")

This now uses Julia Dates module. Timezones are currently not done - everything's in UTC.

Convert time in Julian Days to a formatted string.

    The general format is:

        YYYY-MMM-DD HH:MM:SS ZZZ

    Truncate the time value to seconds, minutes, hours or days as
    indicated. if level = "day", don't print the time zone string.

    Pass an empty string ("", the default) for zone if you want to do
    your own zone formatting in the calling module.

    Parameters:
        jd    : Julian Day number
        zone  : Time zone string (default = "")
        level : "day" or "hour" or "minute" or "second" (default = "second")

    Return:
        formatted date/time string

"""
function lt_to_str(jd, zone="", level="second")
    (yr, mon, day) = jd_to_cal(jd)
    (fday, iday) = modf(day)
    iday =  floor(Integer, iday)
    (hr, mn, second) = fday_to_hms(fday)
    second =  floor(Integer, second)

    if level == "second"
        return Dates.format(DateTime(yr, mon, iday, hr, mn, second), Dates.DateFormat("e, dd u yyyy HH:MM:SS"))
    end

    if level == "minute"
        return Dates.format(DateTime(yr, mon, iday, hr, mn), Dates.DateFormat("e, dd u yyyy HH:MM"))
    end

    if level == "day"
        return Dates.format(DateTime(yr, mon, iday), Dates.DateFormat("e, dd u yyyy"))
    end

    error("Unknown time level = $level")
end

"""
    ut_to_lt(jd)

Convert universal time in Julian Days to a time.

Include Daylight Savings Time offset, if any.

Parameters:
    jd : Julian Day number, universal time

Return:
    Julian Day number, time
    zone string of the zone used for the conversion

not yet doing time zones... :(

"""
function ut_to_lt(jd)
    global daylight_timezone_name, daylight_timezone_offset, standard_timezone_name, standard_timezone_offset
    if is_dst(jd)
        zone   = daylight_timezone_name
        offset = daylight_timezone_offset
    else
        zone   = standard_timezone_name
        offset = standard_timezone_offset
    end
    return (jd - offset/24, zone)
end


text:
@testset "Topology" begin
  @testset "FullTopology" begin
    # 2 triangles
    elems = connect.([(1,2,3),(4,3,2)])
    t = FullTopology(elems)
    @test paramdim(t) == 2
    @test connec4elem(t, 1) == (1,2,3)
    @test connec4elem(t, 2) == (4,3,2)
    @test nvertices(t) == 4
    @test nelements(t) == 2
    @test nfaces(t, 2) == 2
    @test nfaces(t, 1) == 0
    @test nfaces(t, 0) == 4

    # 2 triangles + 2 quadrangles
    elems = connect.([(1,2,6,5),(2,4,6),(4,3,5,6),(1,5,3)])
    t = FullTopology(elems)
    @test connec4elem(t, 1) == (1,2,6,5)
    @test connec4elem(t, 2) == (2,4,6)
    @test connec4elem(t, 3) == (4,3,5,6)
    @test connec4elem(t, 4) == (1,5,3)
    @test nvertices(t) == 6
    @test nelements(t) == 4
    @test nfaces(t, 2) == 4
    @test nfaces(t, 1) == 0
    @test nfaces(t, 0) == 6

    # 1 triangle + 3 quadrangles + 1 triangle hole
    elems = connect.([(1,2,6,5),(2,4,7,6),(4,3,7),(3,1,5,7)])
    t = FullTopology(elems)
    @test connec4elem(t, 1) == (1,2,6,5)
    @test connec4elem(t, 2) == (2,4,7,6)
    @test connec4elem(t, 3) == (4,3,7)
    @test connec4elem(t, 4) == (3,1,5,7)
    @test nvertices(t) == 7
    @test nelements(t) == 4
    @test nfaces(t, 2) == 4
    @test nfaces(t, 1) == 0
    @test nfaces(t, 0) == 7
  end

  @testset "GridTopology" begin
    t = GridTopology(3)
    @test paramdim(t) == 1
    @test size(t) == (3,)
    @test rank2type(t, 1) == Segment
    @test_throws AssertionError rank2type(t, 2)
    @test elem2cart(t, 1) == (1,)
    @test elem2cart(t, 2) == (2,)
    @test elem2cart(t, 3) == (3,)
    @test cart2corner(t, 1) == 1
    @test cart2corner(t, 2) == 2
    @test cart2corner(t, 3) == 3
    @test elem2corner(t, 1) == 1
    @test elem2corner(t, 2) == 2
    @test elem2corner(t, 3) == 3
    @test corner2elem(t, 1) == 1
    @test corner2elem(t, 2) == 2
    @test corner2elem(t, 3) == 3
    @test nelements(t) == 3
    @test nfacets(t) == 4
    @test nvertices(t) == 4
    @test nfaces(t, 1) == 3
    @test nfaces(t, 0) == 4
    @test element(t, 1) == connect((1,2))
    @test element(t, 2) == connect((2,3))
    @test element(t, 3) == connect((3,4))
    @test faces(t, 1) == elements(t)

    t = GridTopology(3, 4)
    @test paramdim(t) == 2
    @test size(t) == (3, 4)
    @test rank2type(t, 1) == Segment
    @test rank2type(t, 2) == Quadrangle
    @test_throws AssertionError rank2type(t, 3)
    @test elem2cart(t, 1) == (1, 1)
    @test elem2cart(t, 2) == (2, 1)
    @test elem2cart(t, 3) == (3, 1)
    @test elem2cart(t, 4) == (1, 2)
    @test elem2cart(t, 5) == (2, 2)
    @test elem2cart(t, 6) == (3, 2)
    @test elem2cart(t, 7) == (1, 3)
    @test elem2cart(t, 8) == (2, 3)
    @test elem2cart(t, 9) == (3, 3)
    @test elem2cart(t, 10) == (1, 4)
    @test elem2cart(t, 11) == (2, 4)
    @test elem2cart(t, 12) == (3, 4)
    @test cart2corner(t, 1, 1) == 1
    @test cart2corner(t, 2, 1) == 2
    @test cart2corner(t, 3, 1) == 3
    @test cart2corner(t, 1, 2) == 5
    @test cart2corner(t, 2, 2) == 6
    @test cart2corner(t, 3, 2) == 7
    @test cart2corner(t, 1, 3) == 9
    @test cart2corner(t, 2, 3) == 10
    @test cart2corner(t, 3, 3) == 11
    @test cart2corner(t, 1, 4) == 13
    @test cart2corner(t, 2, 4) == 14
    @test cart2corner(t, 3, 4) == 15
    @test elem2corner(t, 1) == 1
    @test elem2corner(t, 2) == 2
    @test elem2corner(t, 3) == 3
    @test elem2corner(t, 4) == 5
    @test elem2corner(t, 5) == 6
    @test elem2corner(t, 6) == 7
    @test elem2corner(t, 7) == 9
    @test elem2corner(t, 8) == 10
    @test elem2corner(t, 9) == 11
    @test elem2corner(t, 10) == 13
    @test elem2corner(t, 11) == 14
    @test elem2corner(t, 12) == 15
    @test corner2elem(t, 1) == 1
    @test corner2elem(t, 2) == 2
    @test corner2elem(t, 3) == 3
    @test corner2elem(t, 5) == 4
    @test corner2elem(t, 6) == 5
    @test corner2elem(t, 7) == 6
    @test corner2elem(t, 9) == 7
    @test corner2elem(t, 10) == 8
    @test corner2elem(t, 11) == 9
    @test corner2elem(t, 13) == 10
    @test corner2elem(t, 14) == 11
    @test corner2elem(t, 15) == 12
    @test nelements(t) == 12
    @test nfacets(t) == 31
    @test nvertices(t) == 20
    @test nfaces(t, 2) == 12
    @test nfaces(t, 1) == 31
    @test nfaces(t, 0) == 20
    @test element(t, 1) == connect((1,2,6,5))
    @test element(t, 5) == connect((6,7,11,10))
    @test faces(t, 2) == elements(t)
    @test facet.(Ref(t), 1:31) == connect.([(1, 2), (2, 3), (3, 4), (5, 6), (6, 7),
                                            (7, 8), (9, 10), (10, 11), (11, 12), (13, 14),
                                            (14, 15), (15, 16), (17, 18), (18, 19), (19, 20),
                                            (1, 5), (2, 6), (3, 7), (4, 8), (5, 9),
                                            (6, 10), (7, 11), (8, 12), (9, 13), (10, 14),
                                            (11, 15), (12, 16), (13, 17), (14, 18), (15, 19),
                                            (16, 20)])

    t = GridTopology(3, 4, 2)
    @test paramdim(t) == 3
    @test size(t) == (3, 4, 2)
    @test rank2type(t, 1) == Segment
    @test rank2type(t, 2) == Quadrangle
    @test rank2type(t, 3) == Hexahedron
    @test_throws AssertionError rank2type(t, 4)
    @test elem2cart(t, 1) == (1, 1, 1)
    @test elem2cart(t, 2) == (2, 1, 1)
    @test elem2cart(t, 3) == (3, 1, 1)
    @test elem2cart(t, 4) == (1, 2, 1)
    @test elem2cart(t, 5) == (2, 2, 1)
    @test elem2cart(t, 6) == (3, 2, 1)
    @test elem2cart(t, 7) == (1, 3, 1)
    @test elem2cart(t, 8) == (2, 3, 1)
    @test elem2cart(t, 9) == (3, 3, 1)
    @test elem2cart(t, 10) == (1, 4, 1)
    @test elem2cart(t, 11) == (2, 4, 1)
    @test elem2cart(t, 12) == (3, 4, 1)
    @test elem2cart(t, 13) == (1, 1, 2)
    @test elem2cart(t, 14) == (2, 1, 2)
    @test elem2cart(t, 15) == (3, 1, 2)
    @test elem2cart(t, 16) == (1, 2, 2)
    @test elem2cart(t, 17) == (2, 2, 2)
    @test elem2cart(t, 18) == (3, 2, 2)
    @test elem2cart(t, 19) == (1, 3, 2)
    @test elem2cart(t, 20) == (2, 3, 2)
    @test elem2cart(t, 21) == (3, 3, 2)
    @test elem2cart(t, 22) == (1, 4, 2)
    @test elem2cart(t, 23) == (2, 4, 2)
    @test elem2cart(t, 24) == (3, 4, 2)
    @test cart2corner(t, 1, 1, 1) == 1
    @test cart2corner(t, 2, 1, 1) == 2
    @test cart2corner(t, 3, 1, 1) == 3
    @test cart2corner(t, 1, 2, 1) == 5
    @test cart2corner(t, 2, 2, 1) == 6
    @test cart2corner(t, 3, 2, 1) == 7
    @test cart2corner(t, 1, 3, 1) == 9
    @test cart2corner(t, 2, 3, 1) == 10
    @test cart2corner(t, 3, 3, 1) == 11
    @test cart2corner(t, 1, 4, 1) == 13
    @test cart2corner(t, 2, 4, 1) == 14
    @test cart2corner(t, 3, 4, 1) == 15
    @test cart2corner(t, 1, 1, 2) == 21
    @test cart2corner(t, 2, 1, 2) == 22
    @test cart2corner(t, 3, 1, 2) == 23
    @test cart2corner(t, 1, 2, 2) == 25
    @test cart2corner(t, 2, 2, 2) == 26
    @test cart2corner(t, 3, 2, 2) == 27
    @test cart2corner(t, 1, 3, 2) == 29
    @test cart2corner(t, 2, 3, 2) == 30
    @test cart2corner(t, 3, 3, 2) == 31
    @test cart2corner(t, 1, 4, 2) == 33
    @test cart2corner(t, 2, 4, 2) == 34
    @test cart2corner(t, 3, 4, 2) == 35
    @test elem2corner(t, 1) == 1
    @test elem2corner(t, 2) == 2
    @test elem2corner(t, 3) == 3
    @test elem2corner(t, 4) == 5
    @test elem2corner(t, 5) == 6
    @test elem2corner(t, 6) == 7
    @test elem2corner(t, 7) == 9
    @test elem2corner(t, 8) == 10
    @test elem2corner(t, 9) == 11
    @test elem2corner(t, 10) == 13
    @test elem2corner(t, 11) == 14
    @test elem2corner(t, 12) == 15
    @test elem2corner(t, 13) == 21
    @test elem2corner(t, 14) == 22
    @test elem2corner(t, 15) == 23
    @test elem2corner(t, 16) == 25
    @test elem2corner(t, 17) == 26
    @test elem2corner(t, 18) == 27
    @test elem2corner(t, 19) == 29
    @test elem2corner(t, 20) == 30
    @test elem2corner(t, 21) == 31
    @test elem2corner(t, 22) == 33
    @test elem2corner(t, 23) == 34
    @test elem2corner(t, 24) == 35
    @test corner2elem(t, 1) == 1
    @test corner2elem(t, 2) == 2
    @test corner2elem(t, 3) == 3
    @test corner2elem(t, 5) == 4
    @test corner2elem(t, 6) == 5
    @test corner2elem(t, 7) == 6
    @test corner2elem(t, 9) == 7
    @test corner2elem(t, 10) == 8
    @test corner2elem(t, 11) == 9
    @test corner2elem(t, 13) == 10
    @test corner2elem(t, 14) == 11
    @test corner2elem(t, 15) == 12
    @test corner2elem(t, 21) == 13
    @test corner2elem(t, 22) == 14
    @test corner2elem(t, 23) == 15
    @test corner2elem(t, 25) == 16
    @test corner2elem(t, 26) == 17
    @test corner2elem(t, 27) == 18
    @test corner2elem(t, 29) == 19
    @test corner2elem(t, 30) == 20
    @test corner2elem(t, 31) == 21
    @test corner2elem(t, 33) == 22
    @test corner2elem(t, 34) == 23
    @test corner2elem(t, 35) == 24
    @test nelements(t) == 24
    @test nfacets(t) == 3*24 + 3*4 + 4*2 + 3*2
    @test nvertices(t) == 60
    @test nfaces(t, 3) == 24
    @test nfaces(t, 2) == 3*24 + 3*4 + 4*2 + 3*2
    @test nfaces(t, 0) == 60
    @test element(t, 1) == connect((1,2,6,5,21,22,26,25), Hexahedron)
    @test element(t, 5) == connect((6,7,11,10,26,27,31,30), Hexahedron)
    @test faces(t, 3) == elements(t)
  end

  @testset "HalfEdgeTopology" begin
    function test_halfedge(elems, topology)
      @test nelements(topology) == length(elems)
      for e in 1:nelements(topology)
        he = half4elem(topology, e)
        inds = indices(elems[e])
        @test he.elem == e
        @test he.head ∈ inds
      end
    end

    # 2 triangles as a list of half-edges
    h1  = HalfEdge(1, 1)
    h2  = HalfEdge(2, nothing)
    h3  = HalfEdge(2, 1)
    h4  = HalfEdge(3, 2)
    h5  = HalfEdge(3, 1)
    h6  = HalfEdge(1, nothing)
    h7  = HalfEdge(2, 2)
    h8  = HalfEdge(4, nothing)
    h9  = HalfEdge(4, 2)
    h10 = HalfEdge(3, nothing)
    h1.half =  h2; h2.half  = h1
    h3.half =  h4; h4.half  = h3
    h5.half =  h6; h6.half  = h5
    h7.half =  h8; h8.half  = h7
    h9.half = h10; h10.half = h9
    h1.prev = h5;  h1.next = h3
    h3.prev = h1;  h3.next = h5
    h4.prev = h9;  h4.next = h7
    h5.prev = h3;  h5.next = h1
    h7.prev = h4;  h7.next = h9
    h9.prev = h7;  h9.next = h4
    halves  = [(h1,h2),(h3,h4),(h5,h6),(h7,h8),(h9,h10)]
    struc = HalfEdgeTopology(halves)
    @test half4elem(struc, 1) == h1
    @test half4elem(struc, 2) == h4
    @test half4vert(struc, 1) == h1
    @test half4vert(struc, 2) == h3
    @test half4vert(struc, 3) == h4
    @test half4vert(struc, 4) == h9
    @test edge4pair(struc, (1,2)) == 1
    @test edge4pair(struc, (2,1)) == 1
    @test edge4pair(struc, (2,3)) == 2
    @test edge4pair(struc, (3,2)) == 2
    @test edge4pair(struc, (3,1)) == 3
    @test edge4pair(struc, (1,3)) == 3
    @test edge4pair(struc, (2,4)) == 4
    @test edge4pair(struc, (4,2)) == 4
    @test edge4pair(struc, (4,3)) == 5
    @test edge4pair(struc, (3,4)) == 5

    # 2 triangles
    elems = connect.([(1,2,3),(4,3,2)])
    t = HalfEdgeTopology(elems)
    @test paramdim(t) == 2
    @test nelements(t) == 2
    @test nfacets(t) == 5
    @test nvertices(t) == 4
    @test nfaces(t, 2) == 2
    @test nfaces(t, 1) == 5
    @test nfaces(t, 0) == 4
    test_halfedge(elems, t)

    # 2 triangles + 2 quadrangles
    elems = connect.([(1,2,6,5),(2,4,6),(4,3,5,6),(1,5,3)])
    t = HalfEdgeTopology(elems)
    @test paramdim(t) == 2
    @test nelements(t) == 4
    @test nfacets(t) == 9
    @test nvertices(t) == 6
    @test nfaces(t, 2) == 4
    @test nfaces(t, 1) == 9
    @test nfaces(t, 0) == 6
    test_halfedge(elems, t)

    # 1 triangle + 3 quadrangles + 1 triangle hole
    elems = connect.([(1,2,6,5),(2,4,7,6),(4,3,7),(3,1,5,7)])
    t = HalfEdgeTopology(elems)
    @test paramdim(t) == 2
    @test nelements(t) == 4
    @test nfacets(t) == 11
    @test nvertices(t) == 7
    @test nfaces(t, 2) == 4
    @test nfaces(t, 1) == 11
    @test nfaces(t, 0) == 7
    test_halfedge(elems, t)
  end
end


text:
"""
This function calculates the annualized capital, queue and fixed O&M costs
"""
function calculate_annualized_costs(project::P,
                                    sz::Int64,
                                    queue_cost::Vector{Float64}) where P <: Project{<: BuildPhase}

    finance = get_finance_data(project)
    r = get_discount_rate(finance)

    adjusted_investment_cost = get_effective_investment_cost(finance) * (1 + r) ^ (get_lag_time(finance) - 1)

    adjusted_queue_cost = 0

    for q in 1:length(queue_cost)
        adjusted_queue_cost += queue_cost[q] * (1 + r) ^ (get_lag_time(finance) + (length(queue_cost) - q))
    end

    α =  r / (1 - (1 + r) ^ (-get_capex_years(finance)))

    annualized_cap_cost = fill(α * (adjusted_investment_cost + adjusted_queue_cost), sz)

    annual_OM_cost = fill(get_fixed_OM_cost(finance), sz)

    total_annual_cost = annualized_cap_cost + annual_OM_cost

    return total_annual_cost

end

"""
This function returns the net present value of the existing and planned projects:
NPV for these buildphases doesn't include investment and queue costs.
"""
function calculate_npv(project::P,
                       scenario_name::String,
                       iteration_year::Int64,
                       yearly_horizon::Int64,
                       queue_time::Int64,
                       queue_cost::Vector{Float64},
                       update_start_year::Int64) where P <: Union{Project{Existing}, Project{Planned}}

    project_life_end = get_end_life_year(project)
    finance = get_finance_data(project)
    npv_years = project_life_end - iteration_year + 1
    discount_rate = get_discount_rate(finance)

    OM_cost_array = create_OMcost_array(finance, scenario_name, project_life_end, update_start_year, iteration_year)

    npv_cost_array = [OM_cost_array[y]
                     for y in iteration_year:project_life_end]

    npv_profit_array = [sum(get_scenario_profit(finance)[scenario_name][iteration_year][:, y])
                       for y in iteration_year:project_life_end]

    npv = sum(((1 / (1 + discount_rate)) ^ y) * (npv_profit_array[y] - npv_cost_array[y])
              for y in 1:npv_years)

    return npv
end

"""
This function returns the net present value of the queue and option projects:
NPV for these buildphases considers annualized capital and queue costs.
"""
function calculate_npv(project::P,
                       scenario_name::String,
                       iteration_year::Int64,
                       yearly_horizon::Int64,
                       queue_time::Int64,
                       queue_cost::Vector{Float64},
                       update_start_year::Int64) where P <: Union{Project{Option}, Project{Queue}}

    project_life_end = get_end_life_year(project)
    finance = get_finance_data(project)
    npv_years = min(project_life_end, yearly_horizon + iteration_year - 1)
    discount_rate = get_discount_rate(finance)

    queue_cost_array = create_queue_cost_array(size(get_scenario_profit(finance)[scenario_name][iteration_year], 2),
                                              get_decision_year(project),
                                              queue_time,
                                              queue_cost)

    npv_profit_array = [sum(get_scenario_profit(finance)[scenario_name][iteration_year][:, y])
                        for y in get_construction_year(project):npv_years]

    npv_cost_array = calculate_annualized_costs(project, length(npv_profit_array), queue_cost)

    npv = sum(((1 / (1 + discount_rate)) ^ (y + get_construction_year(project) - iteration_year )) * (npv_profit_array[y] - npv_cost_array[y])
              for y in 1:length(npv_profit_array))

    return npv
end

"""
This function updates the NPV of the project.
"""
function update_project_npv!(project::P,
                               scenario_data::Vector{Scenario},
                               market_prices::MarketPrices,
                               rep_hour_weight::Vector{Float64},
                               iteration_year::Int64,
                               yearly_horizon::Int64,
                               queue_cost::Vector{Float64},
                               capacity_forward_years::Int64,
                               solver::JuMP.MOI.OptimizerWithAttributes) where P <: Project{<: BuildPhase}

    price_years = find_price_years(get_construction_year(project),
                                   get_end_life_year(project),
                                   iteration_year,
                                   yearly_horizon)

    update_years = find_update_years(get_construction_year(project),
                                     get_end_life_year(project),
                                     iteration_year,
                                     yearly_horizon)

    update_expected_profit!(project,
                            scenario_data,
                            market_prices,
                            price_years,
                            update_years,
                            rep_hour_weight,
                            queue_cost,
                            capacity_forward_years,
                            iteration_year,
                            solver)
    expected_npv = 0.0

    finance_data = get_finance_data(project)

    for scenario in scenario_data
        scenario_name = get_name(scenario)
        scenario_npv = calculate_npv(project,
                            scenario_name,
                            iteration_year,
                            yearly_horizon,
                            length(queue_cost),
                            queue_cost,
                            update_years[:start_year])

        if -1 <= scenario_npv / get_maxcap(project) < 0   # Converting values from very small negative values to 0 to avoid rounding off issues
            scenario_npv = 0.0
        end
        set_scenario_npv!(finance_data, scenario_name, iteration_year, scenario_npv)
        expected_npv += scenario_npv * get_probability(scenario)
    end

    set_expected_npv!(finance_data, iteration_year, expected_npv)

    return
end

"""
This function does nothing if the project is retired.
Returns nothing.
"""
function update_project_npv!(project::P,
                               scenario_data::Vector{Scenario},
                               market_prices::MarketPrices,
                               rep_hour_weight::Vector{Float64},
                               iteration_year::Int64,
                               yearly_horizon::Int64,
                               queue_cost::Vector{Float64},
                               capacity_forward_years::Int64,
                               solver::JuMP.MOI.OptimizerWithAttributes ) where P <: Project{Retired}
    return
end


text:
module Test_gsdec

using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

println("Test_gsdec")
@testset "gsdec " begin

sys = rdss(0,0,0);

for job in ("finite", "infinite", "stable", "unstable")
    @time sys1, sys2 = gsdec(sys,job = job)
    @test iszero(sys-sys1-sys2)
end

fast = true; Ty = Complex{Float64}; Ty = Float64
m = 2; n = 5; p = 3;

for fast in (true, false)
# random examples
for Ty in (Float64, Complex{Float64})

# continuous, standard
sys = rss(n,p,m,T = Ty,disc=false);
@time sys1, sys2 = gsdec(sys,fast = fast, job = "finite", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)
@time sys1, sys2 = gsdec(sys,fast = fast, job = "infinite", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)
@time sys1, sys2 = gsdec(sys,fast = fast, job = "stable", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)    &&
      all(real(gpole(sys1)) .< 0) && all(real(gpole(sys2)) .> 0)
@time sys1, sys2 = gsdec(sys,fast = fast, job = "unstable", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)   &&
      all(real(gpole(sys2)) .< 0) && all(real(gpole(sys1)) .> 0)

# discrete, standard
sys = rss(n,p,m,T = Ty,disc=true);
@time sys1, sys2 = gsdec(sys,fast = fast, job = "finite", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)
@time sys1, sys2 = gsdec(sys,fast = fast, job = "infinite", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)
@time sys1, sys2 = gsdec(sys,fast = fast, job = "stable", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)    &&
      all(abs.(gpole(sys1)) .< 1) && all(abs.(gpole(sys2)) .> 1)
@time sys1, sys2 = gsdec(sys,fast = fast, job = "unstable", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)   &&
      all(abs.(gpole(sys2)) .< 1) && all(abs.(gpole(sys1)) .> 0)

# continuous, descriptor, no infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false);
@time sys1, sys2 = gsdec(sys,fast = fast, job = "finite", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)
@time sys1, sys2 = gsdec(sys,fast = fast, job = "infinite", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)
@time sys1, sys2 = gsdec(sys,fast = fast, job = "stable", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)    &&
      all(real(gpole(sys1)) .< 0) && all(real(gpole(sys2)) .> 0)
@time sys1, sys2 = gsdec(sys,fast = fast, job = "unstable", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)   &&
      all(real(gpole(sys2)) .< 0) && all(real(gpole(sys1)) .> 0)

# discrete, descriptor, no infinite eigenvalues
sys = rdss(n,p,m,T = Ty, disc=true);
@time sys1, sys2 = gsdec(sys,fast = fast, job = "finite", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)
@time sys1, sys2 = gsdec(sys,fast = fast, job = "infinite", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)
@time sys1, sys2 = gsdec(sys,fast = fast, job = "stable", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)    &&
      all(abs.(gpole(sys1)) .< 1) && all(abs.(gpole(sys2)) .> 1)
@time sys1, sys2 = gsdec(sys,fast = fast, job = "unstable", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)   &&
      all(abs.(gpole(sys2)) .< 1) && all(abs.(gpole(sys1)) .> 0)

# continuous, descriptor, infinite poles
sys = rdss(n,p,m,T = Ty, disc=false,id=[3*ones(Int,1);2*ones(Int,1)]);
@time sys1, sys2 = gsdec(sys,fast = fast, job = "finite", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)   &&
      all(isinf.(gpole(sys2,atol=1.e-7))) && all(isfinite.(gpole(sys1,atol=1.e-7)))
@time sys1, sys2 = gsdec(sys,fast = fast, job = "infinite", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)   &&
      all(isinf.(gpole(sys1,atol=1.e-7))) && all(isfinite.(gpole(sys2,atol=1.e-7)))
@time sys1, sys2 = gsdec(sys,fast = fast, job = "stable", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)    &&
      all(real(gpole(sys1)) .< 0) && all(real(gpole(sys2,atol=1.e-7)) .> 0)
@time sys1, sys2 = gsdec(sys,fast = fast, job = "unstable", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)   &&
      all(real(gpole(sys2)) .< 0) && all(real(gpole(sys1,atol=1.e-7)) .> 0)

# discrete, descriptor, infinite poles
sys = rdss(n,p,m,T = Ty, disc=true,id=[3*ones(Int,1);2*ones(Int,1)]);
@time sys1, sys2 = gsdec(sys,fast = fast, job = "finite", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)   &&
      all(isinf.(gpole(sys2,atol=1.e-7))) && all(isfinite.(gpole(sys1,atol=1.e-7)))
@time sys1, sys2 = gsdec(sys,fast = fast, job = "infinite", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)    &&
      all(isinf.(gpole(sys1,atol=1.e-7))) && all(isfinite.(gpole(sys2,atol=1.e-7)))
@time sys1, sys2 = gsdec(sys,fast = fast, job = "stable", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)    &&
      all(abs.(gpole(sys1)) .< 1) && all(abs.(gpole(sys2,atol=1.e-7)) .> 1)
@time sys1, sys2 = gsdec(sys,fast = fast, job = "unstable", atol = 1.e-7);
@test iszero(sys-sys1-sys2, atol=1.e-7)   &&
      all(abs.(gpole(sys2)) .< 1) && all(abs.(gpole(sys1,atol=1.e-7)) .> 1)

end
end

end # gsdec

end # module


text:
using SmartShockFinder

datatyp = "CNN"
data_input_size = 4
n_meshes = 1
polydeg = collect(2:7)

generate_traindataset2d(datatyp, n_meshes, polydeg, data_input_size)


text:
import Pkg; Pkg.activate(".")

const USE_GPU = false
using Test
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2);
else
    @init_parallel_stencil(Threads, Float64, 2);
end

_viscosity(x, y) = exp(1 - ( (1-x) + y*(1-y) ))

@parallel function viscosity!(η, x, y)
    @all(η) = _viscosity(@idx_x(x), @idx_y(y))
    return
end

@parallel function viscosity2!(η, x, y)
    @all(η) = _viscosity(@all(x), @all(y))
    return
end

viscosity(x, y) = [ _viscosity(xi, yi) for xi in x, yi in y]

function testme()

    nx = ny = 64
    x = y = LinRange(0.0, 1.0, nx)
    η1 = @zeros(nx, ny)
    @parallel viscosity!(η1, x, y)
    η2 = [ _viscosity(xi, yi) for xi in x, yi in y]

    return η1 == η2
end

@test  testme() == true


text:
using BenchmarkTools
using CatViews
import Base.Cartesian: @ntuple

A = randn(5000,5000)
x = CatView( @ntuple 500 (n)->randn(10) )

@benchmark ($A)*($x)


text:
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

# Method (1) real
t, n, k=10, 20, 10
Xset = [genDataMatrix(t, n) for i = 1:k]
Xfixed=randn(t, n)./1
for i=1:length(Xset) Xset[i]+=Xfixed end
C1=Hermitian( mean((X'*X)/t for X∈Xset) )
C2=Hermitian( mean((X*X')/n for X∈Xset) )
Xbar=mean(Xset)
c=cstp(Xbar, C1, C2; simple=true)
@test c.F[1]'*C2*c.F[1]≈I
@test c.F[2]'*C1*c.F[2]≈I
Z=c.F[1]'*Xbar*c.F[2]
n=minimum(size(Z))
@test norm(Z[1:n, 1:n]-Diagonal(Z[1:n, 1:n]))+1. ≈ 1.
cX=cstp(Xset; simple=true)
@test c==cX

# Method (1) complex
t, n, k=10, 20, 10
Xcset = [genDataMatrix(ComplexF64, t, n) for i = 1:k]
Xcfixed=randn(ComplexF64, t, n)./1
for i=1:length(Xcset) Xcset[i]+=Xcfixed end
C1c=Hermitian( mean((Xc'*Xc)/t for Xc∈Xcset) )
C2c=Hermitian( mean((Xc*Xc')/n for Xc∈Xcset) )
Xcbar=mean(Xcset)
cc=cstp(Xcbar, C1c, C2c; simple=true)
@test cc.F[1]'*C2c*cc.F[1]≈I
@test cc.F[2]'*C1c*cc.F[2]≈I
Zc=cc.F[1]'*Xcbar*cc.F[2]
n=minimum(size(Zc))
@test norm(Zc[1:n, 1:n]-Diagonal(Zc[1:n, 1:n]))+1. ≈ 1.
cXc=cstp(Xcset; simple=true)
@test cc==cXc

# Method (2) real
c=cstp(Xset)

# ... selecting subspace dimension allowing an explained variance = 0.9
c=cstp(Xset; eVar=0.9)

# ... giving weights `w` to the covariance matrices
c=cstp(Xset; w=abs2.(randn(k)), eVar=0.9)

# ... subtracting the means
c=cstp(Xset; meanXd₁=nothing, meanXd₂=nothing, w=abs2.(randn(k)), eVar=0.9)

# explained variance
c.eVar

# name of the filter
c.name

using Plots
# plot the original covariance matrices and the transformed counterpart
c=cstp(Xset)

C1Max=maximum(abs.(C1));
 h1 = heatmap(C1, clim=(-C1Max, C1Max), title="C1", yflip=true, c=:bluesreds);
 D1=c.F[1]'*C2*c.F[1];
 D1Max=maximum(abs.(D1));
 h2 = heatmap(D1, clim=(0, D1Max), title="F[1]'*C2*F[1]", yflip=true, c=:amp);
 C2Max=maximum(abs.(C2));
 h3 = heatmap(C2, clim=(-C2Max, C2Max), title="C2", yflip=true, c=:bluesreds);
 D2=c.F[2]'*C1*c.F[2];
 D2Max=maximum(abs.(D2));
 h4 = heatmap(D2, clim=(0, D2Max), title="F[2]'*C1*F[2]", yflip=true, c=:amp);

XbarMax=maximum(abs.(Xbar));
 h5 = heatmap(Xbar, clim=(-XbarMax, XbarMax), title="Xbar", yflip=true, c=:bluesreds);
 DX=c.F[1]'*Xbar*c.F[2];
 DXMax=maximum(abs.(DX));
 h6 = heatmap(DX, clim=(0, DXMax), title="F[1]'*Xbar*F[2]", yflip=true, c=:amp);
 📈=plot(h1, h3, h5, h2, h4, h6, size=(800,400))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigCSTP.png")

# Method (2) complex
cc=cstp(Xcset)

# ... selecting subspace dimension allowing an explained variance = 0.9
cc=cstp(Xcset; eVar=0.9)

# ... giving weights `w` to the covariance matrices
cc=cstp(Xcset; w=abs2.(randn(k)), eVar=0.9)

# ... subtracting the mean
cc=cstp(Xcset; meanXd₁=nothing, meanXd₂=nothing,
        w=abs2.(randn(k)), eVar=0.9)

# explained variance
c.eVar

# name of the filter
c.name


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

function tri_lu(u,num,lu)
    tri_lu_ = load_op_and_grad("./build/libTriLu","tri_lu")
    u,num,lu = convert_to_tensor([u,num,lu], [Float64,Int64,Int64])
    tri_lu_(u,num,lu)
end

lu = 1
num = 15
m = 5
n = 10
u = rand(1,m, n)
ref = zeros(size(u,1), m, n)
for i = 1:size(u,1)
    if lu==0
        ref[i,:,:] = triu(u[i,:,:], num)
    else
        ref[i,:,:] = tril(u[i,:,:], num)
    end
end
out = tri_lu(u,num,lu)
sess = Session(); init(sess)
@show run(sess, out)- ref

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(m)
    return sum(tri_lu(m,num,lu)^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand(10, m, n))
v_ = rand(10, m, n)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
module BehaviorTree
using AbstractTrees

abstract type BT end

struct Sequence <: BT
    tasks
    name::String
end
struct Selector <: BT
    tasks
    name::String
end
Sequence(tasks) = Sequence(tasks, "")
Selector(tasks) = Selector(tasks, "")
function run_task(task::Function)
    @debug("RUNNING task $(task)")
    result = task()
    return result, result
end
function run_task(task::Function, state)
    @debug("RUNNING task $(task)")
    result = task(state)
    return result, result
end
function run_task(task::BT)
    tick(task)
end

function run_task(task::BT, state)
    tick(task, state)
end

function format(task::BT)
    task.name
end
function format(task::Function)
    string(task)
end

function sequence(tree::Sequence, task_runner)
    @debug("RUNNING sequence $(tree.name)")
    results = []
    for task in tree.tasks
        result, status = task_runner(task)
        push!(results, status)
        if result == :running
            @info("sequence $(tree.name) running at $(format(task))")
            return :running, results
        end
        if result == :failure
            @info("sequence $(tree.name) failed at $(format(task))")
            return :failure, results
        end
    end
    return :success, results
end

function tick(tree::Sequence)
    task_runner(x) = run_task(x)
    sequence(tree, task_runner)
end

function tick(tree::Sequence, state)
    task_runner(x) = run_task(x, state)
    sequence(tree, task_runner)
end

function selector(tree::Selector, task_runner)
    @debug("RUNNING selector $(tree.name)")
    results = []
    for task in tree.tasks
        result, status = task_runner(task)
        push!(results, status)
        if result == :running
            @info("selector $(tree.name) running at $(format(task))")
            return :running, results
        end
        if result == :success
            @info("selector $(tree.name) succeeded at $(format(task))")
            return :success, results
        end
    end
    return :failure, results
end
function tick(tree::Selector)
    task_runner(x) = run_task(x)
    selector(tree, task_runner)
end

function tick(tree::Selector, state)
    task_runner(x) = run_task(x, state)
    selector(tree, task_runner)
end

include("abstractrees.jl")
include("visualization.jl")

export tick, Sequence, Selector
end # module


text:
using Pkg
Pkg.update()
Pkg.add("Plots")


text:
using HypothesisTests, Test
using HypothesisTests: default_tail

@testset "Binomial" begin
t = BinomialTest(26, 78)
@test pvalue(t) ≈ 0.004334880883507431
@test pvalue(t, tail=:left) ≈ 0.002167440441753716
@test pvalue(t, tail=:right) ≈ 0.9989844298129187
@test default_tail(t) == :both
@test_ci_approx confint(t) (0.23058523962930383, 0.4491666887959782)
@test_ci_approx confint(t, tail=:left) (0.0, 0.4313047758370174)
@test_ci_approx confint(t, tail=:right) (0.2451709633730693, 1.0)
@test_ci_approx confint(t, method=:wald) (0.22871819521037956, 0.43794847145628707)
@test_ci_approx confint(t, tail=:left, method=:wald) (0.0, 0.42112912485444692)
@test_ci_approx confint(t, tail=:right, method=:wald) (0.24553754181221971, 1.0)
@test_ci_approx confint(t, method=:wilson) (0.23872670036358601, 0.44358590287381217)
@test_ci_approx confint(t, tail=:left, method=:wilson) (0.0, 0.42541288951088108)
@test_ci_approx confint(t, tail=:right, method=:wilson) (0.25242832328277831, 1.0)
@test_ci_approx confint(t, method=:jeffrey) (0.23626570247518358, 0.44251318323879296)
@test_ci_approx confint(t, tail=:left, method=:jeffrey) (0.0, 0.42466492683653623)
@test_ci_approx confint(t, tail=:right, method=:jeffrey) (0.25098836986261724, 1.0)
@test_ci_approx confint(t, method=:agresti_coull) (0.2384423809121706, 0.44387022232522744)
@test_ci_approx confint(t, tail=:left, method=:agresti_coull) (0.0, 0.42558712894362222)
@test_ci_approx confint(t, tail=:right, method=:agresti_coull) (0.25225408385003706, 1.0)
@test_ci_approx confint(t, method=:arcsine) (0.23366209634204066,0.44117918327686334)
@test_ci_approx confint(t, tail=:left, method=:arcsine) (0.0,0.4235046425920888)
@test_ci_approx confint(t, tail=:right, method=:arcsine) (0.2489264087216164,1.0)

show(IOBuffer(), t)

t = BinomialTest([trues(6); falses(3)])
@test pvalue(t) ≈ 0.5078125000000002
@test pvalue(t, tail=:left) ≈ 0.91015625
@test pvalue(t, tail=:right) ≈ 0.2539062500000001
@test_ci_approx confint(t) (0.2992950562085405, 0.9251453685803082)
@test_ci_approx confint(t, tail=:left) (0.0, 0.9022531865607242)
@test_ci_approx confint(t, tail=:right) (0.3449413659437032, 1.0)
@test_ci_approx confint(t, method=:wald) (0.35868803903340479, 0.97464529429992841)
@test_ci_approx confint(t, tail=:left, method=:wald) (0.0,0.92513047859481645)
@test_ci_approx confint(t, tail=:right, method=:wald) (0.40820285473851681,1.0)
@test_ci_approx confint(t, method=:wilson) (0.35420213558039609,0.87941618161308899)
@test_ci_approx confint(t, tail=:left, method=:wilson) (0.0,0.85802909820500495)
@test_ci_approx confint(t, tail=:right, method=:wilson) (0.39825972868840931,1.0)
@test_ci_approx confint(t, method=:jeffrey) (0.34779179347226591,0.89578677833922582)
@test_ci_approx confint(t, tail=:left, method=:jeffrey) (0.0,0.86830830610561005)
@test_ci_approx confint(t, tail=:right, method=:jeffrey) (0.39604343455469687,1.0)
@test_ci_approx confint(t, method=:agresti_coull) (0.350905767251112,0.88271254994237336)
@test_ci_approx confint(t, tail=:left, method=:agresti_coull) (0.0,0.86049746046629294)
@test_ci_approx confint(t, tail=:right, method=:agresti_coull) (0.39579136642712159,1.0)
@test_ci_approx confint(t, method=:arcsine) (0.345812446615087,0.9188773496172281)
@test_ci_approx confint(t, tail=:left, method=:arcsine) (0.0,0.8879439981269358)
@test_ci_approx confint(t, tail=:right, method=:arcsine) (0.3965293068864491,1.0)
show(IOBuffer(), t)

x = [55, 58, 61, 61, 62, 62, 62, 63, 63, 64, 66, 68, 68, 69, 69, 69, 70, 71, 72, 72]
@test pvalue(SignTest(x)) ≈ 1.907348632812499e-6
@test pvalue(SignTest(x), tail=:left) ≈ 1.0
@test pvalue(SignTest(x), tail=:right) ≈ 9.536743164062495e-7
@test pvalue(SignTest(x, 70)) ≈ 0.004425048828125003
@test pvalue(SignTest(x, 70), tail=:left) ≈ 0.0022125244140625013
@test pvalue(SignTest(x, 70), tail=:right) ≈ 0.9996356964111328
@test default_tail(SignTest(x)) == :both
@test_ci_approx confint(SignTest(x, 70)) (62, 69)
@test_ci_approx confint(SignTest(x, 70), 0.0002) (61, 71)
show(IOBuffer(), SignTest(x, 70))

x = [9, 2, 7, 5]
y = [7, 2, 6, 4]
@test pvalue(SignTest(x, y)) ≈ 0.25
@test_ci_approx confint(SignTest(x, y)) (0, 2)
show(IOBuffer(), SignTest(x, y))

# www.stat.umn.edu/geyer/old03/5102/notes/rank.pdf
x = [-4.7, 3.7, 22.4, 13.6, 8.7, 9.1, -7.8, 10.8, 15.6, 23.5, 14.4, 20.2, 6.5, 10.1, -6.9]
@test pvalue(SignTest(x)) ≈ 0.03515625
@test pvalue(SignTest(x), tail=:left) ≈ 0.996307373046875
@test pvalue(SignTest(x), tail=:right) ≈ 0.017578125000000007
@test_ci_approx confint(SignTest(x), 1-0.99995) (-7.8, 23.5)
@test_ci_approx confint(SignTest(x), 1-0.9999) (-6.9, 22.4)
@test_ci_approx confint(SignTest(x), 1-0.993) (-4.7, 20.2)
@test_ci_approx confint(SignTest(x), 1-0.965) (3.7, 15.6)
@test_ci_approx confint(SignTest(x), 1-0.882) (6.5, 14.4)
@test_ci_approx confint(SignTest(x), 1-0.7) (8.7, 13.6)
@test_ci_approx confint(SignTest(x), 1-0.6) (9.1, 10.8)
show(IOBuffer(), SignTest(x))
end


text:
using Base.Test
using MIPVerify: setloglevel!
using MIPVerify: remove_cached_models
using MIPVerify: get_max_index, get_norm
using JuMP

isdefined(:TestHelpers) || include("TestHelpers.jl")
using TestHelpers: get_new_model

@testset "MIPVerify" begin
    MIPVerify.setloglevel!("info")
    MIPVerify.remove_cached_models()

    include("integration.jl")
    include("net_components.jl")
    include("utils.jl")
    include("models.jl")
    include("batch_processing_helpers.jl")

    @testset "get_max_index" begin
        @test_throws MethodError get_max_index([])
        @test get_max_index([3]) == 1
        @test get_max_index([3, 1, 4]) == 3
        @test get_max_index([3, 1, 4, 1, 5, 9, 2]) == 6
    end

    @testset "get_norm" begin
        @testset "real-valued arrays" begin
            xs = [1, -2, 3]
            @test get_norm(1, xs) == 6
            @test get_norm(2, xs) == sqrt(14)
            @test get_norm(Inf, xs) == 3
            @test_throws DomainError get_norm(3, xs)
        end
        @testset "variable-valued arrays" begin
            @testset "l1" begin
                m = get_new_model()
                x1 = @variable(m, lowerbound=1, upperbound=5)
                x2 = @variable(m, lowerbound=-8, upperbound=-2)
                x3 = @variable(m, lowerbound=3, upperbound=10)
                xs = [x1, x2, x3]
                n_1 = get_norm(1, xs)
                n_2 = get_norm(2, xs)
                n_inf = get_norm(Inf, xs)

                @objective(m, Min, n_1)
                solve(m)
                @test getobjectivevalue(m)≈6

                if Pkg.installed("Gurobi") != nothing
                    # Skip these tests if Gurobi is not installed.
                    # Cbc does not solve problems with quadratic objectives
                    @objective(m, Min, n_2)
                    solve(m)
                    @test getobjectivevalue(m)≈14
                end

                @objective(m, Min, n_inf)
                solve(m)
                @test getobjectivevalue(m)≈3

                @test_throws DomainError get_norm(3, xs)
            end
        end
    end
end


text:
function search(s::AbstractString, c::Char)
    result = findfirst(isequal(c), s)
    result != nothing ? result : 0
end

function inqtext(x, y, s, wc=false)
    if length(s) >= 2 && s[1] == '$' && s[end] == '$'
        tbx, tby = GR.inqmathtex(x, y, s[2:end-1])
    elseif search(s, '\\') != 0 || search(s, '_') != 0 || search(s, '^') != 0
        tbx, tby = GR.inqtextext(x, y, s)
    else
        tbx, tby = GR.inqtext(x, y, s)
    end
    if wc
        for i = 1:4
            tbx[i], tby[i] = GR.ndctowc(tbx[i], tby[i])
        end
    end
    tbx, tby
end

function stringsize(s, wc=false)
    tbx, tby = inqtext(0, 0, s, wc)
    w = tbx[2] - tbx[1]
    h = tby[4] - tby[1]
    w, h
end

function text(x, y, s, wc = false)
    if wc
        win = GR.inqwindow()
        vp = GR.inqviewport()
        x = (vp[2] - vp[1]) / (win[2] - win[1]) * x + (vp[1] * win[2] - vp[2] * win[1]) / (win[2] - win[1])
        y = (vp[4] - vp[3]) / (win[4] - win[3]) * y + (vp[3] * win[4] - vp[4] * win[3]) / (win[4] - win[3])
    end
    if length(s) >= 2 && s[1] == '$' && s[end] == '$'
        GR.mathtex(x, y, s[2:end-1])
    elseif search(s, '\\') != 0 || search(s, '_') != 0 || search(s, '^') != 0
        GR.textext(x, y, s)
    else
        GR.text(x, y, s)
    end
end


text:
push!(LOAD_PATH, "C:\\Users\\Zhang\\Documents\\GitHub\\pcreo\\src")

using LinearAlgebra: I, ×
using Random: shuffle!
using UUIDs

using MultiFloats
using NearestNeighbors: KDTree, knn
using ProgressMeter: @showprogress

using DZLinearAlgebra
using RieszEnergy

use_clean_multifloat_arithmetic()
setprecision(64 * 3)

function split_into_runs(items::AbstractArray{T,N};
                         by=identity, threshold=zero(T)) where {T,N}
    runs = Vector{T}[]
    @inbounds if !isempty(items)
        current_run = T[items[1]]
        current_score = by(items[1])
        for i = 2 : length(items)
            item = items[i]
            next_score = by(item)
            if abs(next_score - current_score) <= threshold
                push!(current_run, item)
            else
                push!(runs, current_run)
                current_run = T[item]
            end
            current_score = next_score
        end
        push!(runs, current_run)
    end
    return runs
end

################################################################################

const NUM_BITS = 90
const EPSILON = inv(Float64x3(2)^NUM_BITS)

################################################################################

function candidate_isometries(p1, p2, q1, q2)
    return [
        hcat(q1, q2, q1 × q2) / hcat(p1, p2, p1 × p2),
        hcat(q1, q2, q1 × q2) / hcat(p1, p2, p2 × p1),
        hcat(q1, q2, q1 × q2) / hcat(p2, p1, p1 × p2),
        hcat(q1, q2, q1 × q2) / hcat(p2, p1, p2 × p1),
    ]
end

function riesz_energy_distances!(
        result::AbstractVector{Tuple{T,Int,Int}},
        points::AbstractMatrix{T}) where {T<:Real}
    dim, num_points = size(points)
    energy = zero(T)
    k = 0
    @inbounds for i = 1 : num_points
        for j = i+1 : num_points
            norm_sq = zero(T)
            @simd ivdep for d = 1 : dim
                temp = points[d,i] - points[d,j]
                norm_sq += temp * temp
            end
            k += 1
            dist = renormalize(unsafe_sqrt(norm_sq))
            energy += inv(dist)
            result[k] = (dist, i, j)
        end
    end
    return energy
end

function isocmp(x::Tuple{T,Int,Int}, y::Tuple{T,Int,Int}) where {T<:Real}
    @inbounds return MultiFloats._lt(x[1], y[1])
end

function isisomorphic(p::Matrix{T}, q::Matrix{T};
                      verbose::Bool=false) where {T <: Real}
    if size(p) != size(q)
        if verbose
            println("DIFFERENT SIZES")
        end
        return false
    end
    dim, num_points = size(p)
    num_pairs = div(num_points * (num_points - 1), 2)
    threshold = EPSILON * dim * num_points
    p_dist = Vector{Tuple{T,Int,Int}}(undef, num_pairs)
    q_dist = Vector{Tuple{T,Int,Int}}(undef, num_pairs)
    Ep = riesz_energy_distances!(p_dist, p)
    Eq = riesz_energy_distances!(q_dist, q)
    if !(abs((Ep - Eq) / (Ep + Eq)) <= threshold)
        if verbose
            println("DIFFERENT ENERGIES: ", abs((Ep - Eq) / (Ep + Eq)))
        end
        return false
    end
    sort!(p_dist; lt=isocmp)
    sort!(q_dist; lt=isocmp)
    @inbounds for i = 1 : num_pairs
        if !(abs(p_dist[i][1] - q_dist[i][1]) <= threshold)
            if verbose
                println("DIFFERENT DISTANCES: ",
                        abs(p_dist[i][1] - q_dist[i][1]))
            end
            return false
        end
    end
    p_runs = split_into_runs(p_dist; by=x->x[1], threshold=threshold)
    q_runs = split_into_runs(q_dist; by=x->x[1], threshold=threshold)
    filter!(run -> !(abs(run[end][1] - 2) <= threshold), p_runs)
    filter!(run -> !(abs(run[end][1] - 2) <= threshold), q_runs)
    if length.(p_runs) != length.(q_runs)
        if verbose
            println("DIFFERENT BARCODES")
        end
        return false
    end
    _, run_idx = findmin(length.(p_runs))
    _, pi, pj = p_runs[run_idx][1]
    for (_, qi, qj) in q_runs[run_idx]
        for R in candidate_isometries(p[:,pi], p[:,pj], q[:,qi], q[:,qj])
            @assert maximum(abs.(R' * R - I)) < threshold
            p_tree = KDTree(R * p)
            idxs, dists = knn(p_tree, q, 1)
            idxs = vcat(idxs...)
            dists = vcat(dists...)
            @assert length(idxs) == num_points
            @assert length(dists) == num_points
            if issetequal(idxs, 1:num_points) && (maximum(dists) < threshold)
                return true
            end
        end
    end
    if verbose
        println("DID NOT FIND ISOMETRY")
    end
    return false
end

################################################################################

const SOURCE_DIR = "C:/Users/Zhang/PCREO/PCREO"
const DATABASE_DIR = "C:/PCREO-DB"
const DATABASE_CACHE = Dict{Int,Vector{Tuple{Float64x3,String}}}()

function build_database_cache(num_points::Int)
    @assert !haskey(DATABASE_CACHE, num_points)
    println("Building database cache for size: ", num_points)
    identifier = "PCREO-" * lpad(num_points, 8, '0')
    cache_entry = Tuple{Float64x3,String}[]
    for foldername in readdir(DATABASE_DIR)
        if occursin(identifier, foldername)
            folderpath = joinpath(DATABASE_DIR, foldername)
            @assert isdir(folderpath)
            exemplarpath = joinpath(DATABASE_DIR,
                foldername, first(readdir(folderpath)))
            @assert isfile(exemplarpath)
            exemplar = Float64x3.(BigFloat.(readlines(exemplarpath)))
            @assert length(exemplar) % 3 == 0
            num_points = div(length(exemplar), 3)
            exemplar = reshape(exemplar, (3, num_points))
            @assert norm(riesz_gradient(exemplar)) <= EPSILON * length(exemplar)
            energy = riesz_energy(exemplar)
            push!(cache_entry, (energy, folderpath))
        end
    end
    DATABASE_CACHE[num_points] = cache_entry
end

function find_in_database(points::Matrix{Float64x3})
    dim, num_points = size(points)
    threshold = EPSILON * dim * num_points
    Ep = riesz_energy(points)
    if !haskey(DATABASE_CACHE, num_points)
        build_database_cache(num_points)
    end
    @assert haskey(DATABASE_CACHE, num_points)
    for (Eq, folderpath) in DATABASE_CACHE[num_points]
        if abs((Ep - Eq) / (Ep + Eq)) <= threshold
            @assert isdir(folderpath)
            exemplarpath = joinpath(folderpath, first(readdir(folderpath)))
            @assert isfile(exemplarpath)
            exemplar = Float64x3.(BigFloat.(readlines(exemplarpath)))
            @assert length(exemplar) % 3 == 0
            exemplar = reshape(exemplar, (3, div(length(exemplar), 3)))
            if isisomorphic(points, exemplar)
                return folderpath
            end
        end
    end
    return nothing
end

function add_to_database_entry(sourcepath::String, folderpath::String)
    mv(sourcepath, joinpath(folderpath, basename(sourcepath)))
end

function create_database_entry(points::Matrix{Float64x3}, sourcepath::String)
    dim, num_points = size(points)
    @assert haskey(DATABASE_CACHE, num_points)
    energy = riesz_energy(points)
    foldername = "PCREO-$(lpad(num_points, 8, '0'))-$(uppercase(string(uuid4())))"
    folderpath = joinpath(DATABASE_DIR, foldername)
    push!(DATABASE_CACHE[num_points], (energy, folderpath))
    mkdir(folderpath)
    add_to_database_entry(sourcepath, folderpath)
end

function main()
    for filename in readdir(SOURCE_DIR)
        path = joinpath(SOURCE_DIR, filename)
        if isfile(path)
            println("Reading: ", path)
            lines = readlines(path)
            @assert length(lines) % 3 == 0
            num_points = div(length(lines), 3)
            points = reshape(Float64x3.(BigFloat.(lines)), (3, num_points))
            if norm(riesz_gradient(points)) <= EPSILON * length(points)
                folderpath = find_in_database(points)
                if isnothing(folderpath)
                    println("Creating new database entry.")
                    create_database_entry(points, path)
                else
                    println("Matched existing database entry: ", folderpath)
                    add_to_database_entry(path, folderpath)
                end
            else
                println("WARNING: Gradient norm exceeds threshold.")
            end
        end
    end
end

# function readpoints(path)
#     lines = readlines(path)
#     @assert length(lines) % 3 == 0
#     reshape(Float64x3.(BigFloat.(lines)),
#             (3, div(length(lines), 3)))
# end

# function main()
#     for foldername in readdir(DATABASE_DIR)
#         println("Verifying directory: ", foldername)
#         folderpath = joinpath(DATABASE_DIR, foldername)
#         @assert isdir(folderpath)
#         exemplarpaths = [joinpath(folderpath, exemplarname)
#                          for exemplarname in readdir(folderpath)]
#         points = readpoints.(exemplarpaths)
#         sizes = unique!(size.(points))
#         @assert length(sizes) == 1
#         dim, num_points = sizes[1]
#         threshold = EPSILON * dim * num_points
#         Emin, Emax = extrema(riesz_energy.(points))
#         @assert (Emax - Emin) / (Emax + Emin) < threshold
#         for i = 1 : length(points)
#             for j = i+1 : length(points)
#                 if !isisomorphic(points[i], points[j]; verbose=true)
#                     println("ERROR:")
#                     println("    ", exemplarpaths[i])
#                     println("    ", exemplarpaths[j])
#                     println(isisomorphic(readpoints(exemplarpaths[i]), readpoints(exemplarpaths[j])))
#                 end
#             end
#         end
#     end
# end

main()


text:
by_la_sys_iteration_and_year = main_results |>
    @filter( _.year > 2018 && _.year < 2025 ) |>
    @groupby( [_.ccode,  _.year, _.sysno, _.iteration ] ) |>
    @map({
        year=first(_.year),
        ccode=first(_.ccode),
        sysno=first(_.sysno),
        iteration=first(_.iteration),
        cnt=length( _ ),
        income=sum( _.income_recieved )*52.0,
        contribs = sum( _.contributions_from_yp )*52.0,
        payments = sum( _.payments_from_la )*52.0
        }

        ) |>
    @groupby( [_.ccode,  _.year, _.sysno] ) |>
    @map(
            {
               year   = first(_.year),
               ccode  = first(_.ccode),
               sysno  = first(_.sysno),

               avg_income=mean( _.income  ),
               min_income=minimum( _.income ),
               max_income=maximum( _.income ),
               pct_10_income=quantile( _.income, [0.10] )[1],
               pct_25_income=quantile( _.income, [0.25] )[1],
               pct_75_income=quantile( _.income, [0.75] )[1],
               pct_90_income=quantile( _.income, [0.90] )[1],

               avg_contribs = mean( _.contribs ),
               min_contribs = minimum( _.contribs ),
               max_contribs = maximum( _.contribs ),
               pct_10_contribs=quantile( _.contribs, [0.10] )[1],
               pct_25_contribs=quantile( _.contribs, [0.25] )[1],
               pct_75_contribs=quantile( _.contribs, [0.75] )[1],
               pct_90_contribs=quantile( _.contribs, [0.90] )[1],

               avg_payments = mean( _.payments  ),
               min_payments = minimum( _.payments ),
               max_payments = maximum( _.payments ),
               pct_10_payments=quantile( _.payments, [0.10] )[1],
               pct_25_payments=quantile( _.payments, [0.25] )[1],
               pct_75_payments=quantile( _.payments, [0.75] )[1],
               pct_90_payments=quantile( _.payments, [0.90] )[1]
            }
        ) |>
    @orderby( [_.year,_.ccode,_.sysno] ) |>
    DataFrame


text:
using Test, bad

p0  = .2
α = .05
p1 = .4
β = .2
mcr = .3
pmp = .7

prior = (0.2*Beta(1, 1) + 0.8*Beta(mean = .35, sd = .1)) <= pmp

problem = Problem(
    minimise(
        SampleSize(prior)
    ),
    Power(prior  | p0) <= α,
    Power(prior >= mcr)  >= 1 - β
)
design = optimise(problem; verbosity = 0)

𝚾 = sample_space(design)

ordering = EstimatorOrdering(MaximumLikelihoodEstimator(); orientation = :superiority)
pval     = PValue(ordering, design, p0)
pvals    = pval.(𝚾[:,1], 𝚾[:,2])

@test all(
    pvals .≈ p_value.(𝚾[:,1], 𝚾[:,2], p0, ordering, design)
)

prob0 = pmf.(𝚾[:,2], n2.(design, 𝚾[:,1]), p0) .* pmf.(𝚾[:,1], n1(design), p0)

function f(p)
    inds = pvals .<= p
    min(1, max(0, sum(prob0[inds])))
end

p = 0:.01:1
@test all(f.(p) .<= p)


text:
"""
    AbstractParametersBinary <: AbstractParameters
"""
abstract type AbstractParametersBinary <: AbstractParameters end

include("1pl.jl")
include("2pl.jl")
include("3pl.jl")
include("Npl.jl")


text:
################################################################################
# String management for Speckles.jl
################################################################################

"""
    makeName(params::Dict; excpt::Dict = Dict())

Returns a string based on parameters in the simulation dictionary.
"""
function makeName(params::Dict; excpt::Dict = Dict())
    date = today()
    mm = length(month(date)) == 1 ? string(0,month(date)) : month(date)
    dd = length(day(date)) == 1 ? string(0,day(date)) : day(date)
    out = string(year(date),mm,dd,"_")
    for key in keys(params)
        keyname = key
        val = params[key]
        if isa(val,Array)
            val = length(val)
            keyname = string("len-",keyname)
        end
        out = string(out,keyname,"=",val,"_")
    end
    return out
end

export makeName

"""
    function paramTable(params::Dict)

Returns a markdown table with given parameters
"""
function paramTable(params::Dict,names::Dict = Dict())
    out = "| Name | Value |\n|:---:|:---:|\n"
    for (key,value) in params
        paramName = string(key)
        if key in keys(names)
            paramName = string(names[key])
        end
        out = string(out,"|$paramName|$value|\n")
    end
    return out

end

export paramTable


text:
#=
An irrational decimal fraction is created by concatenating the positive integers:

0.123456789101112131415161718192021...

It can be seen that the 12th digit of the fractional part is 1.

If dn represents the nth digit of the fractional part, find the value of the following expression.

d1 × d10 × d100 × d1000 × d10000 × d100000 × d1000000
=#
function champernowne_digit(n)
  index, digit = 1, 1
  while index < n
    d = reverse(digits(digit))
    nextindex = index + length(d)
    if nextindex > n
      return d[(n+1)-index]
    end
    digit += 1
    index += length(d)
  end
  reverse(digits(digit))[1]
end

function calc()
  arr = Int64[]
  for i in [1, 10, 100, 1_000, 10_000, 100_000, 1_000_000]
    push!(arr, champernowne_digit(i))
  end
  prod(arr)
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
# Duffing resonance curves
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot

function amplitude_curve(d,a,ωvec,ntrans,nattr,nres,y_init)
    nω = length(ωvec) # no. of driving frequencies (control parameter)
    amp_vec = zeros(nω)    # oscillation amplitudes (to be computed)
    ds = Systems.duffing(; d, f = a)

    for iω = 1:nω
        ω = ωvec[iω]
        set_parameter!(ds, 1, ω) # sets the ω parameter
        period = 2π/ω
        Δt = period/nres
        tr = trajectory(ds, nattr*period, y_init; Ttr = ntrans*period)
        amp_vec[iω] = maximum(u[1] for u in tr)
        y_init = tr[end]
    end
    return amp_vec
end

@time begin
# parameters and initial values
d = 0.1
omega_min = 0.02
omega_max = 3.3
nomega = 1500 # 1500
ωvec = [LinRange(omega_min, omega_max, nomega); LinRange(omega_max, omega_min, nomega)]
avec = [1.1, 0.2, 0.04]
ccol = [COLORS[3], COLORS[2], COLORS[1]]
ntrans = 400   # no. of transient periods
nattr  = 2     # no. of periods on the attractor
nres   = 300   # resolution
amps = []

# compute resonance curves
for ia = 1:length(avec)
    a = avec[ia]
    y_init = [0.1, 0.1]
    amp_vec = amplitude_curve(d,a,ωvec,ntrans,nattr,nres,y_init)
    push!(amps, amp_vec)
end
end

# %% plot
fig = figure(figsize=(0.55*figx,figy))
ax = gca()
ax.set_xlabel(L"\omega")
ax.set_ylabel("\$x_{\\mathrm{max}}\$")

for ia = 1:length(avec)
    a = avec[ia]
    ax.scatter(ωvec,amps[ia],s=2,c=ccol[ia], label = "\$a=$(avec[ia])\$")
end
ax.set_xlim([0,omega_max])
ax.legend(markerscale = 10)
fig.tight_layout(pad=0.3)

wsave(plotsdir("9", "duffing_resonance"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook