text:
@testset "IEEE1788.a" begin

    @test interval(-Inf,Inf) === entireinterval()

end

@testset "IEEE1788.b" begin

    @test @interval("[1.2345]") === Interval(0x1.3C083126E978DP+0, 0x1.3C083126E978EP+0)

    @test @interval("[1,+infinity]") === Interval(1.0, Inf)

    @test @decorated("[1,1e3]_com") === DecoratedInterval(Interval(1.0, 1000.0), com)

    @test @decorated("[1,1E3]_COM") === DecoratedInterval(Interval(1.0, 1000.0), com)

end

@testset "IEEE1788.c" begin

    @test @interval("[1.e-3, 1.1e-3]") === Interval(0x4.189374BC6A7ECP-12, 0x4.816F0068DB8BCP-12)

    @test @interval("[-0x1.3p-1, 2/3]") === Interval(-0x9.8000000000000P-4, +0xA.AAAAAAAAAAAB0P-4)

    @test @interval("[3.56]") === Interval(0x3.8F5C28F5C28F4P+0, 0x3.8F5C28F5C28F6P+0)

    @test @interval("3.56?1") === Interval(0x3.8CCCCCCCCCCCCP+0, 0x3.91EB851EB8520P+0)

    @test @interval("3.56?1e2") === Interval(355.0, 357.0)

    @test @interval("3.560?2") === Interval(0x3.8ED916872B020P+0, 0x3.8FDF3B645A1CCP+0)

    @test @interval("3.56?") === Interval(0x3.8E147AE147AE0P+0, 0x3.90A3D70A3D70CP+0)

    @test @interval("3.560?2u") === Interval(0x3.8F5C28F5C28F4P+0, 0x3.8FDF3B645A1CCP+0)

    @test @interval("-10?") === Interval(-10.5, -9.5)

    @test @interval("-10?u") === Interval(-10.0, -9.5)

    @test @interval("-10?12") === Interval(-22.0, 2.0)

end

@testset "IEEE1788.d" begin

    @test @interval("[1.234e5,Inf]") === Interval(123400.0, Inf)

    @test @interval("3.1416?1") === Interval(0x3.24395810624DCP+0, 0x3.24467381D7DC0P+0)

    @test @interval("[Empty]") === emptyinterval()

end

@testset "IEEE1788.e" begin

    @test isnai(DecoratedInterval(2,1))

end

@testset "IEEE1788.e" begin

    @test @decorated("[ ]") === DecoratedInterval(emptyinterval(), trv)

    @test @decorated("[entire]") === DecoratedInterval(Interval(-Inf, +Inf), dac)

    @test @decorated("[1.e-3, 1.1e-3]") === DecoratedInterval(Interval(0x4.189374BC6A7ECP-12, 0x4.816F0068DB8BCP-12), com)

    @test @decorated("[-Inf, 2/3]") === DecoratedInterval(Interval(-Inf, +0xA.AAAAAAAAAAAB0P-4), dac)

    @test @decorated("[0x1.3p-1,]") === DecoratedInterval(Interval(0x1.3p-1, Inf), dac)

    @test @decorated("[,]") === DecoratedInterval(entireinterval(), dac)

    @test @decorated("3.56?1") === DecoratedInterval(Interval(0x3.8CCCCCCCCCCCCP+0, 0x3.91EB851EB8520P+0), com)

    @test @decorated("3.56?1e2") === DecoratedInterval(Interval(355.0, 357.0), com)

    @test @decorated("3.560?2") === DecoratedInterval(Interval(0x3.8ED916872B020P+0, 0x3.8FDF3B645A1CCP+0), com)

    @test @decorated("3.56?") === DecoratedInterval(Interval(0x3.8E147AE147AE0P+0, 0x3.90A3D70A3D70CP+0), com)

    @test @decorated("3.560?2u") === DecoratedInterval(Interval(0x3.8F5C28F5C28F4P+0, 0x3.8FDF3B645A1CCP+0), com)

    @test @decorated("-10?") === DecoratedInterval(Interval(-10.5, -9.5), com)

    @test @decorated("-10?u") === DecoratedInterval(Interval(-10.0, -9.5), com)

    @test @decorated("-10?12") === DecoratedInterval(Interval(-22.0, 2.0), com)

    @test @decorated("-10??u") === DecoratedInterval(Interval(-10.0, Inf), dac)

    @test @decorated("-10??") === DecoratedInterval(Interval(-Inf, Inf), dac)

    @test isnai(@decorated("[nai]"))

    @test @decorated("3.56?1_def") === DecoratedInterval(Interval(0x3.8CCCCCCCCCCCCP+0, 0x3.91EB851EB8520P+0), def)

end

@testset "IEEE1788.f" begin

    @test @interval("[]") === emptyinterval()

    @test @interval("[empty]") === emptyinterval()

    @test @interval("[ empty ]") === emptyinterval()

    @test @interval("[,]") === entireinterval()

    @test @interval("[ entire ]") === entireinterval()

end


text:
using AugmentedGaussianProcesses.KernelModule
using Distances, LinearAlgebra
using BenchmarkTools
using Random
if !@isdefined(SUITE)
    const SUITE = BenchmarkGroup(["Kernel"])
end

paramfile = "params/kernel.json"
dim = 50
Random.seed!(1234)
N1 = 1000; N2 = 500;
X = rand(Float64,N1,dim)
Y = rand(Float64,N2,dim)
KXY = rand(Float64,N1,N2)
KX = rand(Float64,N1,N1)
sKX = Symmetric(rand(Float64,N1,N1))
kX = rand(Float64,N1)

kernelnames = ["Matern","RBF"]
kerneltypes = ["ARD","ISO"]
kernels=Dict{String,Dict{String,Kernel}}()
for k in kernelnames
    kernels[k] = Dict{String,Kernel}()
    SUITE["Kernel"][k] = BenchmarkGroup(kerneltypes)
    for kt in kerneltypes
        SUITE["Kernel"][k][kt] = BenchmarkGroup(["XY","XYinplace","X","Xinplace","diagX","diagXinplace","dXY","dX","ddiagX"])
        kernels[k][kt] = eval(Meta.parse(k*"Kernel("*(kt == "ARD" ? "[2.0]" : "2.0" )*",variance=10.0,dim=dim)"))
    end
end

for k in kernelnames
    for kt in kerneltypes
        SUITE["Kernel"][k][kt]["XY"] = @benchmarkable kernelmatrix($X,$Y,$(kernels[k][kt]))
        SUITE["Kernel"][k][kt]["XYinplace"] = @benchmarkable kernelmatrix!(KXY,$X,$Y,$(kernels[k][kt])) setup=(KXY=copy($KXY))
        SUITE["Kernel"][k][kt]["X"] = @benchmarkable kernelmatrix($X,$(kernels[k][kt]))
        SUITE["Kernel"][k][kt]["Xinplace"] = @benchmarkable kernelmatrix!(KX,$X,$(kernels[k][kt])) setup=(KX=copy($KX))
        SUITE["Kernel"][k][kt]["diagX"] = @benchmarkable kerneldiagmatrix($X,$(kernels[k][kt]))
        SUITE["Kernel"][k][kt]["diagXinplace"] = @benchmarkable kerneldiagmatrix!(kX,$X,$(kernels[k][kt])) setup=(kX=copy($kX))
        SUITE["Kernel"][k][kt]["dXY"] = @benchmarkable kernelderivativematrix($X,$Y,$(kernels[k][kt]))
        # SUITE["Kernel"][k][kt]["dXY_K"] = @benchmarkable kernelderivativematrix_K($X,$Y,KXY,$(kernels[k][kt])) setup=(KXY=copy($KXY))
        SUITE["Kernel"][k][kt]["dX"] = @benchmarkable kernelderivativematrix($X,$(kernels[k][kt]))
        # SUITE["Kernel"][k][kt]["dX_K"] = @benchmarkable kernelderivativematrix_K($X,sKX,$(kernels[k][kt])) setup=(sKX=copy($sKX))
        SUITE["Kernel"][k][kt]["ddiagX"] = @benchmarkable kernelderivativediagmatrix($X,$(kernels[k][kt]))
    end
end

# if isfile(paramfile)
#     loadparams!(suite,BenchmarkTools.load(paramfile)[1])
# else
#     tune!(suite,verbose=true)
#     BenchmarkTools.save(paramfile,params(suite))
# end
#
# results = run(suite,verbose=true,seconds=10)
# save_target = "results/kernel_"*("$(now())"[1:10])
# i = 1
# while isfile(save_target*"_$(i).json")
#     global i += 1
# end
# BenchmarkTools.save(save_target*"_$(i).json",results)

ext: jl


text:
## dsolve

## A type akin to SymFunction but with the ability to keep track of derivative
type SymFunction <: SymPy.SymbolicObject
    u::PyCall.PyObject
    n::Int
end

function PyCall.PyObject(f::SymFunction)
    f.n == 0 && return f.u
    __x__ = symbols("__x__")
    diff(f(__x__), __x__, f.n).x
end

"""

Create a symbolic function. These can be used for specifying differential equations.
For these objects we can specify derivatives with the transpose
operator (e.g., `u''`) as opposed to, say `diff(u(x), x, 2)`.

Example:
```
u = SymFunction("u")
u'
```

Alternatively, we can pass `symfunction` to the `cls` argument of
`symbols`. This provides a convenient way to create more than one
symbolic function per call.

```
F, G = symbols("F,G", cls=symfunction)
```

"""
function SymFunction{T<:AbstractString}(x::T)
    u = sympy[:Function](x)
    SymFunction(u, 0)
end

IVPSolution{T<:AbstractString}(x::T) = SymFunction(x)
export IVPSolution
## Need to deprecate this...
@deprecate IVPSolution(x::AbstractString) SymFunction(x)

symfunction(x) = SymFunction(x) # for use with symbols("F", cls=symfunction)

## macro to create functions
"""

Macro to create multiple symbolic functions in the Main workspace.

Example:
```
@symfuns Vd Vq id iq
@syms Rs Lq=>"\lambda_q"  Ld=>"lambda_d" FLq FLd We Wm Np FLf t

basicDQEquations = [
    -Vq(t) + Rs*iq(t) + We*Ld*id(t) + We*FLf + Lq*iq'(t)
    -Vd(t) + Rs*id(t) - We*Lq*iq(t) + Ld*id'(t) ]
```

Thanks to `@alhirzel` for the contribution.
"""
macro symfuns(x...)
    q = Expr(:block)
    fs = []    # running list of functions created
    for s in x
        if isa(s, Expr) && s.head == :(=>) # named function
            push!(fs, s.args[1])
            push!(q.args, Expr(:(=), s.args[1], Expr(:call, :SymFunction, s.args[2])))
        elseif isa(s, Symbol)
            push!(fs, s)
            push!(q.args, Expr(:(=), s, Expr(:call, :SymFunction, string(s))))
        else
            throw(AssertionError("@fns expected a list of symbols"))
        end
    end
    push!(q.args, Expr(:tuple, fs...)) # return all of the functions we created
    eval(Main, q)
end

# some display of objects
Base.display(u::SymFunction) = println("$(string(Sym(u.u)))" * repeat("'", u.n))

## override call is in call.jl

## rather than use `diff(u(x),x,1)` we can use `u'(x)`
function Base.ctranspose(x::SymFunction)
    SymFunction(x.u, x.n + 1)
end

"""

Solve an ordinary differential equation.

Examples:

```
f = SymFunction("f")
x = Sym("x")
dsolve(diff(f(x), x) + f(x)) ## solve f'(x) + f(x) = 0
dsolve(diff(f(x), x, x) + f(x)) ## solve f''(x) + f(x) = 0
```

References: [SymPy Docs](http://docs.sympy.org/0.7.5/modules/solvers/ode.html#ode-docs)

In addition, we add a `Julia`n interface to solve initial-value
problems. (Julia `v"0.4.0"+`) For this, the calling pattern is

`dsolve(eqn::Sym, var::Sym, args::Tuple...; kwargs...)` where

* `eqn` the equation. Can be written as `u'(x) - u(x)` where `u` is of type `SymFunction`
* `var` the variable of the equation. Typically `free_symbols(x)[1]`, but still, for now, is specified
* `args::Tuple...` Specification of initial values in style of `(u, x0, y0)` or `(u', x0, y0)`.
* `kwargs...` are passed onto `dsolve`

Example:

```
u = SymFunction("u")
@vars x
dsolve(u'(x) - 2u(x))                  # u(x) = C_1 e^(2x)
dsolve(u'(x) - 2u(x), x, (u, 0 , 1)) # u(x) = e^(2x)
dsolve(u''(x) - 2u(x), x, (u, 0, 1), (u', 0, 2))  ## some expression
```

This is a simple wrapper for: 1) substituting in the values for the
initial conditions 2) calling `solve` for the constants `C1`, `C2`,
... 3) putting these values back into the equations. If this automated
approach breaks down, then those steps can be done manually starting
with the output of `dsolve` without the initial conditions.

"""

dsolve(ex::Sym;kwargs...) = sympy_meth(:dsolve, ex; kwargs...)
dsolve(exs::Vector{Sym};kwargs...) = sympy_meth(:dsolve, exs; kwargs...)
dsolve(exs::Vector{Sym}, fx::Sym; kwargs...) = sympy_meth(:dsolve, exs, fx; kwargs...)
## Note, dsolve(ex, var::Sym; kwargs...) is depracated so that initial value problems can be specifiedx

## The dsolve function is not great for initial value problems
## The `ics` argument seems to only work with power series solutions
##
## This adds the ability to more naturally specify the equations.
##
## There are issues with the number of solutions that need ironing out.
function dsolve(eqn::Sym, var::Sym, args::Tuple...; kwargs...)

    if length(args) == 0
        throw(ArgumentError("""Some initial value specification is needed.
Specifying the function, as in `dsolve(ex, f(x))`, is deprecated.
Use `sympy_meth(:dsolve, ex, f(x); kwargs...)` directly for that underlying interface.
"""))
    end

    out = dsolve(eqn; kwargs...)
    ord = sympy_meth(:ode_order, eqn, var)

    ## `out` may be an array of solutions. If so we do each one.
    ## we want to use an array for output only if needed
    if !isa(out, Array)
        return _solve_ivp(out, var, args,ord)
    else
        output = Sym[]
        for o in out
            a = _solve_ivp(o, var, args,ord)
            a != nothing && push!(output, a)
        end
        return length(output) == 1 ? output[1] : output
    end
end

## Helper.
## out is an equation in var with constants. Args are intial conditions
## Return `nothing` if initial condition is not satisfied (found by `solve`)
function _solve_ivp(out, var, args, o)
    eqns = Sym[rhs(diff(out, var, f.n))(var=>x0) - y0 for (f, x0, y0) in args]
    sols = solve(eqns, Sym["C$i" for i in 1:o])
    if length(sols) == 0
       return nothing
    end

    ## massage output
    ## Might have more than one solution, though unlikely. But if we substitute a variable
    ## for y0 we will get an array back from solve which may have length 1.
    if isa(sols, Array)
        if length(sols) == 1
            sols = sols[1]
        else
            return [out([Pair(k,v) for (k,v) in sol]...) for sol in sols]
        end
    end

    out([Pair(k,v) for (k,v) in sols]...)
end

export SymFunction, symfunction, dsolve


text:
# TODO temporary fix, delete later
function _nullStringsFix!{T<:AbstractString}(v::NullableVector{T})
    for i ∈ 1:length(v)
        if v.isnull[i]
            v[i] = ""
        end
    end
end

# TODO this  is a temporary fix for a bug in Feather
function _nullStringsFix!(data::DataFrame)
    for (name, dtype) ∈ zip(names(data), eltype.(eltypes(data)))
        if dtype <: AbstractString
            _nullStringsFix!(data[name])
        end
    end
end

# TODO: this will be changed once Feather.jl gets updated
"""
    featherWrite(filename, df[, overwrite=false])

A wrapper for writing dataframes to feather files.  To be used while Feather.jl package
is in development.

If `overwrite`, this will delete the existing file first (an extra step taken to avoid some
strange bugs).
"""
function featherWrite(filename::AbstractString, df::DataFrames.DataFrame;
                      overwrite::Bool=false)::Void
    if isfile(filename)
        if !overwrite
            error("File already exists.  Use overwrite=true.")
        end
        rm(filename)
    end
    Feather.write(filename, df)
    nothing
end

function featherWrite(filename::AbstractString, dt::DataFrame;
                      overwrite::Bool=false)::Void
    dt = copyColumns(dt)
    _nullStringsFix!(dt)  # TODO this is for fixing a bug in Feather, remove when fixed
    featherWrite(filename, convert(DataFrames.DataFrame, dt), overwrite=overwrite)
end
export featherWrite

# TODO these will also get changed once Feather gets updated
# .... other parts of thsi might change in 0.6
"""
    convertWeakRefStrings(df)
    convertWeakRefStrings!(df)

Converts all columns with eltype `Nullable{WeakRefString}` to have eltype `Nullable{String}`.
`WeakRefString` is a special type of string used by the feather package to improve deserialization
performance.

Note that this will no longer be necessary in Julia 0.6.
"""
function convertWeakRefStrings(df::AbstractDataFrame)
    odf = DataFrame()
    for col ∈ names(df)
        if eltype(df[col]) <: Nullable{Feather.WeakRefString{UInt8}}
            odf[col] = convert(NullableVector{String}, df[col])
        else
            odf[col] = df[col]
        end
    end
    odf
end
export convertWeakRefStrings
function convertWeakRefStrings!(df::AbstractDataFrame)
    for col ∈ names(df)
        if eltype(df[col]) <: Nullable{Feather.WeakRefString{UInt8}}
            df[col] = convert(NullableVector{String}, df[col])
        end
    end
    df
end
export convertWeakRefStrings!

# TODO: this will also get changed once Feather.jl gets updated
"""
    featherRead(filename[; convert_strings=true])

A wrapper for reading dataframes which are saved in feather files.  The purpose of this
wrapper is primarily for converting `WeakRefString` to `String`.  This will no longer
be necessary in Julia 0.6.
"""
function featherRead(filename::AbstractString; convert_strings::Bool=true)::DataFrame
    df = convert(DataFrame, Feather.read(filename))  # this shouldn't be needed now
    # df = convert(DataFrame, Feather.read(filename))
    if convert_strings
        convertWeakRefStrings!(df)
    end
    df
end
export featherRead


text:
using Test
using Elemental

for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
    A = Elemental.DistMatrix(T)
    Elemental.gaussian!(A, 10)
    @test isa(similar(A), Elemental.DistMatrix{T})
    @test isa(similar(A, Elemental.ElInt), Elemental.DistMatrix{Elemental.ElInt})
    @test isa(similar(A, Float32), Elemental.DistMatrix{Float32})
    @test isa(similar(A, Elemental.ElInt, (5,5)), Elemental.DistMatrix{Elemental.ElInt})
    @test isa(similar(A, Float32, (5,5)), Elemental.DistMatrix{Float32})
end


text:
"""
`Mesh{N,V,F} <: AbstractMesh{N,V,F}`

`N`-dimensional surface-mesh with faces of type `F` and `N`-dimensional vertices of type `V`

Mesh(vertices::Vector{<:AbstractPoint}, faces::Vector{<:AbstractPolygon})
"""
struct Mesh{N,V,F} <: AbstractMesh{N,V,F}
  vertices :: OffsetVector{V}
  faces :: Vector{F}
end

"""
`TriangularMesh{N,V}`

`N`-dimensional surface-mesh with triangular faces and `N`-dimensional vertices of type `V`
"""
const TriangleMesh{N,V} = Mesh{N,V,Triangle}

Mesh(vertices::Vector{V}, faces::Vector{F}) where {N,V<:AbstractPoint{N},F} =
  Mesh{N,V,F}(zerobased(vertices), faces)

face(mesh::Mesh, i) = mesh.faces[i]
vertex(mesh::Mesh, i) = mesh.vertices[i]

Base.length(mesh::Mesh) = length(mesh.faces)
Base.size(mesh::Mesh) = size(mesh.faces)
Base.getindex(mesh::Mesh, i) = face(mesh, i)


text:
# AS clause for a common table expression.

mutable struct CTEClause <: AbstractSQLClause
    over::Union{SQLClause, Nothing}
    name::Symbol
    columns::Union{Vector{Symbol}, Nothing}
    materialized::Union{Bool, Nothing}

    CTEClause(;
              over = nothing,
              name::Union{Symbol, AbstractString},
              columns::Union{AbstractVector{<:Union{Symbol, AbstractString}}, Nothing} = nothing,
              materialized = nothing) =
        new(over,
            Symbol(name),
            !(columns isa Union{Vector{Symbol}, Nothing}) ?
                Symbol[Symbol(col) for col in columns] : columns,
            materialized)
end

CTEClause(name; over = nothing, columns = nothing, materialized = nothing) =
    CTEClause(over = over, name = name, columns = columns, materialized = materialized)

"""
    CTE(; over = nothing, name, columns = nothing, materialized = nothing)
    CTE(name; over = nothing, columns = nothing, materialized = nothing)

An `AS` clause for a common table expression.

# Examples

```jldoctest
julia> c = FROM(:condition_occurrence) |>
           WHERE(OP("=", :condition_concept_id, 320128)) |>
           SELECT(:person_id) |>
           CTE(:essential_hypertension);

julia> print(render(c))
"essential_hypertension" AS (
  SELECT "person_id"
  FROM "condition_occurrence"
  WHERE ("condition_concept_id" = 320128)
)
```
"""
CTE(args...; kws...) =
    CTEClause(args...; kws...) |> SQLClause

dissect(scr::Symbol, ::typeof(CTE), pats::Vector{Any}) =
    dissect(scr, CTEClause, pats)

function PrettyPrinting.quoteof(c::CTEClause, ctx::QuoteContext)
    ex = Expr(:call, nameof(CTE), quoteof(c.name))
    if c.columns !== nothing
        push!(ex.args, Expr(:kw, :columns, Expr(:vect, quoteof(c.columns, ctx)...)))
    end
    if c.materialized !== nothing
        push!(ex.args, Expr(:kw, :materialized, c.materialized))
    end
    if c.over !== nothing
        ex = Expr(:call, :|>, quoteof(c.over, ctx), ex)
    end
    ex
end

rebase(c::CTEClause, c′) =
    CTEClause(over = rebase(c.over, c′), name = c.name, columns = c.columns, materialized = c.materialized)


text:
# this implements a very simple spiking layer mechanism - it simply emits a single spike when input > 0
# common use case: drives actual biophysical layers in a heterogenous network with a source of spikes
# For now, this completely ignores incoming synaptic inputs (efficiency reasons)

function reset!(neuron::SimpleSpikingNeuron)
    neuron.lastSpike = -Inf
    neuron.didSpike = false
end

# type signature of update function should match others for operator overloading
# this is because simplespiking might not be a source layer (for whatever reason)
function update!(neuron::SimpleSpikingNeuron, doSpike::Bool, t::Float64)
	neuron.didSpike = doSpike
	if (doSpike)
		neuron.lastSpike = t
	end
end


text:
#Import matlab files and saving them as .JLD2 files!

#Kristen R. Hunter-Cevera, Marine Biological Laboratory, Spring 2020

using MAT
using Dates
using JLD2
include("/home/kristen/Documents/useful_julia_scripts/save_by_regex.jl")

file = matopen("/mnt/Lab_data/MVCO/FCB/Syn_and_MVCO_packaged_data/syndata_current.mat")
varnames = names(file)

#assign variables one by one:
#mu_avg=read(file, "mu_avg") # note that this does NOT introduce a variable ``varname`` into scope

# @eval $(Symbol(string)) = $value or @eval $(Symbol(string)) = value.

## Success! Finally!
#to_add=[] #does this need to allocated first?
for var in varnames

    if occursin("time",var) #need to deal with time differently as these are matlab dates/times
        @show var

        ex = :(temp = read(file,$var)) #import variable as temporary
        eval(ex)

        global to_add = temp .- temp[1,1] #reset to 0 days first entry should be Jan-1-03!
        @show typeof(to_add)
        v=Symbol(var)
        ex=:($v = DateTime(2003,1,1) + map(x->Dates.Day(x), to_add)) #sweet! Success!
        eval(ex)

    else
        v=Symbol(var)
        ex = :($v = read(file,$var))
        #@show ex
        eval(ex)
    end
end

close(file)
##

#saving with regex! thanks to Dehann!
save("/home/kristen/Documents/V6V8_analysis/analysis_products/mvco_syndata.jld2", [r"syn",r"mu",r"net",r"loss",r"PE",r"SSC",r"vol",r"CHL",r"allgrowthrates",r"allMR",r"modelallmatdate",r"allmatdat",r"FCB"])
#*syn* *mu* *net* *loss* *PE* *SSC* *vol* *CHL* allgrowthrates allMR modelallmatdate allmatdat* FCBnum
#older, default save:
#@save("/home/kristen/Documents/V6V8_analysis/scripts/analysis_products/mvco_syndata.jld2")

## and the environmental data!

file = matopen("/mnt/Lab_data/MVCO/FCB/Syn_and_MVCO_packaged_data/mvco_envdata_current.mat")
varnames = names(file)

#to_add=[] #does this need to allocated first?
for var in varnames

    if occursin("time",var)#need to deal with time differently as these are matlab dates/times
        @show var

        ex = :(temp = read(file,$var)) #import variable as temporary
        eval(ex)
        global to_add = temp .- temp[1,1] #reset to 0 days first entry should be Jan-1-03!
        # @show typeof(to_add)
        v=Symbol(var)

        uu=findall(map(x->isnan(x),temp))
        if isempty(uu)

            to_add = map(x->Dates.Day(x), to_add)
            @show typeof(to_add)
            ex=:($v = DateTime(2003,1,1) + to_add) #sweet! Success!
            @show ex
            eval(ex)

        end
    else
        v=Symbol(var)
        ex = :($v = read(file,$var))
        #@show ex
        eval(ex)
    end
end

close(file)
##
save("/home/kristen/Documents/V6V8_analysis/analysis_products/mvco_envdata.jld2",[r"light", r"Tcorr", r"corr_years", r"Tbeam_corr", r"dy_", r"_month nutyears", r"nutyrdy", r"time_nut", r"SiOH", r"NH4", r"NO3", r"PO4"])
#*light* *Tcorr* corr_years Tbeam_corr dy_*  *_month nutyears nutyrdy time_nut* SiOH_wk NH4_wk NO3_wk PO4_wk SiOH_*_mn NH4_*_mn NO3_*_mn PO4_*_mn'
#@save "/home/kristen/Documents/V6V8_analysis/scripts/analysis_products/mvco_envdata.jld2"

#EVAL AND EXPRESSIONS:
# var="mu_avg"
# v=Symbol(var)
# ex = :($v = read(file,$var))
# eval(ex)

# a=1
# b=3
# ex=:(a+b)
# eval(ex)

#we need to work with time....so, easy way?
#731582.0 = Jan 1 2003 > work it!
# to_add = time_syn .- time_syn[1,1]
#
# map(x->Dates.Day(x), to_add)
# a= DateTime(2003,1,1) + map(x->Dates.Day(x), to_add) #sweet! Success!

## to match mvco event numbers:
file = matopen("/home/kristen/Documents/seasons_of_syn_paper/matlab_scripts/nut_data_reps.mat")
varnames = names(file)
for var in varnames
    v=Symbol(var)
    ex = :($v = read(file,$var))
    #@show ex
    eval(ex)
end
close(file)
##

@save "/home/kristen/Documents/V6V8_analysis/scripts/analysis_products/mvco_nut_data.jld2" MVCO_nut_reps header_nut

### plotting....
# gr()
# scatter(ydmu, mu_avg)
# xaxis!("Year day",(0, 366),0:100:360,font(12, "Times"))
# # xlabel = "my label",
# #    xlims = (0,10),
# #    xticks = 0:0.5:10,
# #    xscale = :log,
# #    xflip = true,
# #    xtickfont = font(20, "Courier"))
#
# plot(time_syn,daily_syn, yscale = :log10, size=(600,200), framestyle =:box, legend=:outertopright)

max_stars_repo_name: hsosik/NES-LTER

max_issues_repo_name: hsosik/NES-LTER

max_forks_repo_name: hsosik/NES-LTER


text:
# function to parse the time-series data from the given folder and output a dictionary
# implementing the data into the  original RTS_GMLC network.
#
# function parameters:
######################################################################################
# month_st: the month to start
# day_st: the day to start
# period_st: the hour to start
# month_en: the month that ends parsing
# day_en: the month that ends parsing
# period_en: the month that ends parsing (here the periods need to be chosen corresponding to the data set use in the future)
# folder: the folder containing the time-sereis data and subnetworks.csv(For example I set it to be time_series)

function hour_parse(month_st, day_st, period_st, month_en, day_en, period_en, folder)

    # set the time range
    start_t = DateTime(2020,month_st,day_st,period_st)
    end_t = DateTime(2020,month_en,day_en,period_en)
    range_t = start_t:Hour(1):end_t

    # direct to the current directory
    # println("$(@__DIR__)")
    cur_dir_split = split(rstrip(folder, '/'), "/")
    cur_dir_split = cur_dir_split[cur_dir_split .!= ""]
    cur_dir_new = cur_dir_split[1:end-3]
    # cur_dir = join([("$i/") for i in cur_dir_new])
    cur_dir = pwd()

    # for the base structure of network, parse all the subnetwork files
    subnet_file = CSV.read("$(folder)/subnetworks.csv", DataFrame)
    subnet_file_1 = subnet_file[!, "file"][1]
    case = Dict()
    for subnet in eachrow(subnet_file)
        case[subnet.file] = parse_file("$(folder)/$(subnet.file)")
    end

    # Set the generator minimum active power limits to zero for surrent stage (release the constraints)
    for (subnet_file, subnet) in case
       for (gen_idx, gen) in subnet["gen"]
          gen["pmin"] = 0
       end
    end

    # read the time-series data from the folder
    hour_csp_ori = CSV.read("$(cur_dir)/time_series/CSP/DAY_AHEAD_Natural_Inflow.csv", DataFrame)
    hour_hydro_ori = CSV.read("$(cur_dir)/time_series/Hydro/DAY_AHEAD_hydro.csv", DataFrame)
    hour_load_ori = CSV.read("$(cur_dir)/time_series/Load/DAY_AHEAD_regional_Load.csv", DataFrame)
    hour_pv_ori = CSV.read("$(cur_dir)/time_series/PV/DAY_AHEAD_pv.csv", DataFrame)

    # here didn't include reserve data for now (now optimizations on different time steps are independent)
    hour_rtpv_ori = CSV.read("$(cur_dir)/time_series/RTPV/DAY_AHEAD_rtpv.csv", DataFrame)
    hour_wind_ori = CSV.read("$(cur_dir)/time_series/WIND/DAY_AHEAD_wind.csv", DataFrame)

    hour_wind = filter(row -> DateTime(row.Year, row.Month, row.Day, row.Period) in range_t, hour_wind_ori)
    hour_pv = filter(row -> DateTime(row.Year, row.Month, row.Day, row.Period) in range_t, hour_pv_ori)
    hour_csp = filter(row -> DateTime(row.Year, row.Month, row.Day, row.Period) in range_t, hour_csp_ori)
    hour_rtpv = filter(row -> DateTime(row.Year, row.Month, row.Day, row.Period) in range_t, hour_rtpv_ori)
    hour_load = filter(row -> DateTime(row.Year, row.Month, row.Day, row.Period) in range_t, hour_load_ori)
    hour_hydro = filter(row -> DateTime(row.Year, row.Month, row.Day, row.Period) in range_t, hour_hydro_ori)

    wind_names = names(hour_wind)[5:end]
    pv_names = names(hour_pv)[5:end]
    csp_names = names(hour_csp)[5:end]
    rtpv_names = names(hour_rtpv)[5:end]
    #load_names = names(hour_load)[5:end]
    hydro_names = names(hour_hydro)[5:end]
    #input_csp = filter()

    # set the load in the three areas
    total_load = Dict(1=>0.0, 2=>0.0, 3=>0.0)
    for (id, load) in case[subnet_file_1]["load"]
        bus= case[subnet_file_1]["bus"]["$(load["load_bus"])"]
        total_load[bus["area"]] += load["pd"]*case[subnet_file_1]["baseMVA"]
    end

    hour_data = Dict()
    # made a copy of original network at each time step and make modifications using time-series data.
    for i in 1:size(hour_load, 1)
        load_data = hour_load[i,:]
        for subnet in eachrow(subnet_file)
            # file name and index of current subnetwork file
            name = subnet.file
            hour_data[i] = deepcopy(case)
            hour_data[i][name]["year"] = load_data.Year
            hour_data[i][name]["month"] = load_data.Month
            hour_data[i][name]["day"] = load_data.Day
            hour_data[i][name]["period"] = load_data.Period

        end
        hour_data[i][subnet_file_1]["ori_load"] = Dict(1=>load_data[Symbol(1)], 2=>load_data[Symbol(2)], 3=>load_data[Symbol(3)])    end

    # here distribute the power over the 3 areas
    for (id, name) in hour_data
        #for (ind, network) in name
            for (load_id, load) in name[subnet_file_1]["load"]
                #if haskey(network, "ori_load")
                    area = name[subnet_file_1]["bus"]["$(load["load_bus"])"]["area"]
                    scale = name[subnet_file_1]["ori_load"][area]/total_load[area]
                    load["pd"] = load["pd"]*scale
                    load["qd"] = load["qd"]*scale
                #end
            end
        #end
    end

    # here change all the data of the points according to time-serires data of PV,
    for (id, name) in hour_data
        #for (ind, network) in name
            for (gen_id, gen) in name[subnet_file_1]["gen"]
                #for subnet in eachrow(subnet_file)
                    #if haskey(gen, "name")
                        if gen["name"] in wind_names
                            gen["pmax"] = hour_wind[id, gen["name"]]/case[subnet_file_1]["baseMVA"]
                        end
                        if gen["name"] in pv_names
                            gen["pmax"] = hour_pv[id, gen["name"]]/case[subnet_file_1]["baseMVA"]
                        end
                        if gen["name"] in csp_names
                            gen["pmax"] = hour_csp[id, gen["name"]]/case[subnet_file_1]["baseMVA"]
                        end
                        if gen["name"] in rtpv_names
                            gen["pmax"] = hour_rtpv[id, gen["name"]]/case[subnet_file_1]["baseMVA"]
                        end
                    #else
                    #    gen["pmax"] =
                    #end
                #end
            end
        #end
    end

    # write the data to a specified JSON file for debug
    string = JSON.json(hour_data)
    open("$(folder)/RTS_GMLC_($(month_st).$(day_st).$(period_st)__$(month_en).$(day_en).$(period_en)).json", "w") do f
        write(f, string)
    end

    time_size = size(hour_load, 1)
    return hour_data
end


text:
"""
    AbstractMultivector{T <: Real} <: Number

An abstract low-dimensional multivector.
"""
abstract type AbstractMultivector{T <: Real} <: Number end

function isreal(z::AbstractMultivector)
    iszero(unreal(z))
end

function asarray(z::AbstractMultivector)
    vcat(real(z), unreal(z))
end

function iszero(z::AbstractMultivector)
    iszero(asarray(z))
end

(+)(z::AbstractMultivector{T}) where {T <: Real} = z

"""
    abs(z::AbstractMultivector)

The absolute value of the real part of z.
"""
abs(z::AbstractMultivector) = abs(real(z))

"""
    abs2(z::AbstractMultivector)

The square of the real part of z.
Also known as the quadrance of z.
"""
abs2(z::AbstractMultivector) = (real(z))^2

"""
    iszerodivisor(z::AbstractMultivector)

Returns true if the real part of z is zero.
Note that it follows that also ``z * z = 0``.
That is, if z is a zero-divisor, then it is also nilpotent.
"""
iszerodivisor(z::AbstractMultivector) = iszero(real(z))

"""
    commutator(x::AbstractMultivector, y::AbstractMultivector)
Measures the failure of commutativity of multivector multiplication.
"""
function commutator(x::AbstractMultivector, y::AbstractMultivector)
    (x * y) - (y * x)
end

"""
    associator(x::AbstractMultivector, y::AbstractMultivector, z::AbstractMultivector)

Measures the failure of associativity of multivector multiplication.
"""
function associator(x::AbstractMultivector, y::AbstractMultivector, z::AbstractMultivector)
    ((x * y) * z) - (x * (y * z))
end

"""
    jacobiator(x::AbstractMultivector, y::AbstractMultivector, z::AbstractMultivector)

Measures the failure of multivector multiplication to satisfy the Jacobi identity.
"""
function jacobiator(x::AbstractMultivector, y::AbstractMultivector, z::AbstractMultivector)
    (x * (y * z)) + (y * (z * x)) + (z * (x * y))
end

"""
    inv(z::AbstractMultivector)

Inverse of a multivector. An error occurs if argument is a zero divisor.
"""
function inv(z::AbstractMultivector)
    if iszerodivisor(z)
        error(ZeroDivisorInverse)
    end

    conj(z) / abs2(z)
end

function (/)(a::Real, z::AbstractMultivector)
    if iszerodivisor(z)
        error(ZeroDivisorDenominator)
    end

    a * inv(z)
end

function (\)(z::AbstractMultivector, a::Real)
    if iszerodivisor(z)
        error(ZeroDivisorDenominator)
    end

    inv(z) * a
end

"""
    crossratioL(w::AbstractMultivector,
               x::AbstractMultivector,
               y::AbstractMultivector,
               z::AbstractMultivector)
The left cross-ratio:
```julia
    ((x-y) \\ (w-y)) * ((w-z) \\ (x-z))
```
The left cross-ratio is invariant under right Möbius transformations.
"""
function crossratioL(w::AbstractMultivector,
                     x::AbstractMultivector,
                     y::AbstractMultivector,
                     z::AbstractMultivector)
    ((x - y) \ (w - y)) * ((w - z) \ (x - z))
end

"""
    crossratioR(w::AbstractMultivector,
                x::AbstractMultivector,
                y::AbstractMultivector,
                z::AbstractMultivector)
The right cross-ratio:
```julia
    ((w-y) / (x-y)) * ((x-z) / (w-z))
```
The right cross-ratio is invariant under left Möbius transformations.
"""
function crossratioR(w::AbstractMultivector,
                     x::AbstractMultivector,
                     y::AbstractMultivector,
                     z::AbstractMultivector)
    ((w - y) / (x - y)) * ((x - z) / (w - z))
end

"""
    möbiusL(z::AbstractMultivector,
           a::AbstractMultivector,
           b::AbstractMultivector,
           c::AbstractMultivector,
           d::AbstractMultivector)
    möbiusL(z::AbstractMultivector, a::Real, b::Real, c::Real, d::Real)
The left Möbius transformation:
```julia
    ((z * c) + d) \\ ((z * a) + b)
```
This transformation is also know as a fractional linear transformation.
"""
function möbiusL(z::AbstractMultivector,
                a::AbstractMultivector,
                b::AbstractMultivector,
                c::AbstractMultivector,
                d::AbstractMultivector)
    ((z * c) + d) \ ((z * a) + b)
end

function möbiusL(z::AbstractMultivector, a::Real, b::Real, c::Real, d::Real)
    ((c * z) + d) \ ((a * z) + b)
end

"""
    möbiusR(z::AbstractMultivector,
           a::AbstractMultivector,
           b::AbstractMultivector,
           c::AbstractMultivector,
           d::AbstractMultivector)
    möbiusR(z::AbstractMultivector, a::Real, b::Real, c::Real, d::Real)
The right Möbius transformation:
```julia
    ((a * z) + b) / ((c * z) + d)
```
This transformation is also know as a fractional linear transformation.
"""
function möbiusR(z::AbstractMultivector,
                a::AbstractMultivector,
                b::AbstractMultivector,
                c::AbstractMultivector,
                d::AbstractMultivector)
    ((a * z) + b) / ((c * z) + d)
end

function möbiusR(z::AbstractMultivector, a::Real, b::Real, c::Real, d::Real)
    ((a * z) + b) / ((c * z) + d)
end


text:
# This file defines the Prerequisite type, used to define all types of
#   prerequisites.

# It is important to note that prerequisites that define an inequality relation,
#   the relation is ALWAYS stated as
#     "prereqValue prereqRelation value_in_personnel_record".
# For example, if something requires a minimum age of 21.0, the prerequisite is
#   written as
#     "21.0 <= :age"

# The Prerequisite type requires no extra types.
requiredTypes = [ ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export Prerequisite
"""
This type defines a prerequisite on a single attribute of a personnel member.

The type contains the following fields:
* `attr::String`: the attribute upon which the prerequisite is defined.
* `val::Union{String,Vector{String}}`: the value with which the the attribute is
compared.
* `rel::Function`: the relation between the attribute and the value. Valid
relations are `==`, `!=`, `∈`, and `∉`.

The prerequisites are always of the form `attr rel val`.
"""
type Prerequisite

    attr::String
    val::Union{String, Vector{String}}
    rel::Function

    # Constructors.
    function Prerequisite( attr::String, rel::Function, val::String )

        if rel ∉ [ ==, != ]
            error( "A String value permits only == and != as operators for a Prerequisite." )
        end  # if rel ∉ [ ==, != ]

        newPrereq = new()
        newPrereq.attr = attr
        newPrereq.val = val
        newPrereq.rel = rel
        return newPrereq

    end  # Prerequisite( attr, rel, val )

    function Prerequisite( attr::String, rel::Function, valList::Vector{String} )

        if rel ∉ [ ∈, ∉ ]
            error( "A Vector{String} value permits only ∈ and ∉ as operators for a Prerequisite." )
        end  # if rel ∉ [ ==, != ]

        newPrereq = new()
        newPrereq.attr = attr
        newPrereq.val = valList
        newPrereq.rel = rel
        return newPrereq

    end  # Prerequisite( attr, rel, valList )

end  # type Prerequisite


text:
using DSP:conv; using BenchmarkTools

function boltzconvolve_closure(mod, t, dt, ϵdot)
  mod_vector = mod(t)
  out = conv(mod_vector, ϵdot)
end

function objectivefunction_closure(params, grad, modulus, t, dt, ϵdot)
  mod = (x -> modulus(x, params))
  convolved = boltzconvolve_closure(mod, t, dt, ϵdot)
end

function boltzconvolve_sendparams(mod, params, t, dt, ϵdot)
  mod_vector = mod(t, params)
  out = conv(mod_vector, ϵdot)
end

function objectivefunction_sendparams(params, grad, mod, t, dt, ϵdot)
  convolved = boltzconvolve_sendparams(mod, params, t, dt, ϵdot)
end

function SLS2_gmod(t, params)
  params[1] .+ params[2]*exp.(-t/params[3]) .+ params[4]*exp.(-t/params[5])
end

dt = 0.1
t = Vector{Float64}(0.0:dt:1000.0)
ϵdot = ones(length(t))
params = ones(5)
@btime objectivefunction_closure($params, [0.0], $SLS2_gmod, $t, $dt, $ϵdot);
@btime objectivefunction_sendparams($params, [0.0], $SLS2_gmod, $t, $dt, $ϵdot);

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

#=
---------------------------------------------------------------
versioninfo()
Julia Version 1.3.1
Commit 2d5741174c (2019-12-30 21:36 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: AMD Athlon(tm) II X4 640 Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, amdfam10)
Environment:
  JULIA_EDITOR = "atom.exe"  -a
  JULIA_NUM_THREADS = 4
---------------------------------------------------------------

The execution of this file leads to the following error message:

Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks.
Exception: EXCEPTION_ACCESS_VIOLATION at 0x2b33f08d -- fmpz_set at /workspace/srcdir/flint2/fmpz\set.c:23
in expression starting at C:\Users\Home\JuliaProjects\IntegerSequences.jl\src\SetPartitionsTypeM.jl:539
fmpz_set at /workspace/srcdir/flint2/fmpz\set.c:20
coeff at C:\Users\Home\.julia\packages\Nemo\LBWHq\src\flint\fmpz_poly.jl:38
unknown function (ip: 0000000018263017)
#7 at .\none:0 [inlined]
iterate at .\generator.jl:47 [inlined]
collect at .\array.jl:622
SetPartitionsOfTypeM at C:\Users\Home\JuliaProjects\IntegerSequences.jl\src\Polynomials.jl:56
L036040 at C:\Users\Home\JuliaProjects\IntegerSequences.jl\src\SetPartitionsTypeM.jl:394 [inlined]
macro expansion at C:\Users\Home\JuliaProjects\IntegerSequences.jl\src\SetPartitionsTypeM.jl:456 [inlined]
macro expansion at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\Test\src\Test.jl:1107 [inlined]
test at C:\Users\Home\JuliaProjects\IntegerSequences.jl\src\SetPartitionsTypeM.jl:455
main at C:\Users\Home\JuliaProjects\IntegerSequences.jl\src\SetPartitionsTypeM.jl:534
unknown function (ip: 0000000018233FA0)
_jl_invoke at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\gf.c:2141 [inlined]
jl_apply_generic at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\gf.c:2305
jl_apply at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\julia.h:1631 [inlined]
do_call at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:328
eval_value at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:417
eval_stmt_value at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:368 [inlined]
eval_body at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:764
jl_interpret_toplevel_thunk_callback at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:888
unknown function (ip: FFFFFFFFFFFFFFFE)
unknown function (ip: 000000002934610F)
unknown function (ip: FFFFFFFFFFFFFFFF)
jl_toplevel_eval_flex at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:814
jl_eval_module_expr at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:181
jl_toplevel_eval_flex at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:640
jl_parse_eval_all at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\ast.c:873
jl_load at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:878 [inlined]
jl_load_ at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:885
include at .\boot.jl:328 [inlined]
include_relative at .\loading.jl:1105
include at .\Base.jl:31
_jl_invoke at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\gf.c:2135 [inlined]
jl_apply_generic at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\gf.c:2305
jl_apply at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\julia.h:1631 [inlined]
do_call at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:328
eval_value at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:417
eval_stmt_value at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:368 [inlined]
eval_body at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:764
eval_body at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:705
jl_interpret_toplevel_thunk_callback at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:888
unknown function (ip: FFFFFFFFFFFFFFFE)
unknown function (ip: 0000000010DE970F)
unknown function (ip: 0000000000000002)
jl_toplevel_eval_flex at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:814
jl_toplevel_eval at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:823 [inlined]
jl_toplevel_eval_in at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:843
eval at .\boot.jl:330 [inlined]
eval at .\client.jl:425
_jl_invoke at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\gf.c:2135 [inlined]
jl_apply_generic at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\gf.c:2305
top-level scope at .\none:3
jl_toplevel_eval_flex at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:808
jl_toplevel_eval_flex at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:764
jl_toplevel_eval at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:823 [inlined]
jl_toplevel_eval_in at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:843
eval at .\boot.jl:330
_jl_invoke at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\gf.c:2135 [inlined]
jl_apply_generic at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\gf.c:2305
exec_options at .\client.jl:263
_start at .\client.jl:460
jfptr__start_2086 at C:\Users\Home\AppData\Local\Julia-1.3.1\lib\julia\sys.dll (unknown line)
_jl_invoke at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\gf.c:2135 [inlined]
jl_apply_generic at /cygdrive/d/buildbot/worker/package_win64/build/src/cygdrive/d/buildbot/worker/package_win64/build/src\gf.c:2305
unknown function (ip: 00000000004017DD)
unknown function (ip: 0000000000401BC5)
unknown function (ip: 00000000004013DE)
unknown function (ip: 000000000040151A)
BaseThreadInitThunk at C:\WINDOWS\System32\KERNEL32.DLL (unknown line)
RtlUserThreadStart at C:\WINDOWS\SYSTEM32\ntdll.dll (unknown line)
Allocations: 10894447 (Pool: 10893053; Big: 1394); GC: 12
ERROR: LoadError: Failed to precompile SetPartitionsTypeM [top-level] to C:\Users\Home\.julia\compiled\v1.3\SetPartitionsTypeM.ji.
Stacktrace:
 [1] compilecache(::Base.PkgId, ::String) at .\loading.jl:1283
 [2] _require(::Base.PkgId) at .\loading.jl:1024
 [3] require(::Base.PkgId) at .\loading.jl:922
 [4] require(::Module, ::Symbol) at .\loading.jl:917
in expression starting at C:\Users\Home\JuliaProjects\IntegerSequences.jl\src\TestSetPartitions.jl:8
=#

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module TestSetPartitions

using Nemo, SetPartitionsTypeM, OrderedSetPartitionsTypeM, SeqUtils, OEISUtils

function Test()

    rowsum(f, m, len) = [sum(f(m, k)) for k ∈ 0:len-1]
    altrowsum(f, m, len) = [sum([(iseven(j) ? -c : c) for (j, c) in enumerate(f(m, k))]) for k ∈ 0:len-1]

    for m ∈ 0:4
        println("\n++++++++++++++++++++++++++\n")
        println("\nShapePartitions m = $m")
        for n ∈ 0:6
            ShapePartitionsOfTypeM(m, n) |> Println
        end
        println()
        seq = ShapePartitionsOfTypeM(m, 5)
        oeis_search(seq, 2)
        println()
        println("Row sums:")
        rsum = rowsum(ShapePartitionsOfTypeM, m, 6)
        oeis_search(rsum, 2)
        println()
        println("Alternating row sums:")
        asum = altrowsum(ShapePartitionsOfTypeM, m, 6)
        oeis_search(asum, 2)
        println()

        println()
        println("\nSetPartitions m = $m")
        for n ∈ 0:6
            SetPartitionsOfTypeM(m, n) |> Println
        end
        println()
        seq = SetPartitionsOfTypeM(m, 5)
        oeis_search(seq, 2)
        println()
        println("Row sums:")
        rsum = rowsum(SetPartitionsOfTypeM, m, 6)
        oeis_search(rsum, 2)
        println()
        println("Alternating row sums:")
        asum = altrowsum(SetPartitionsOfTypeM, m, 6)
        oeis_search(asum, 2)
        println()
    end

    for m ∈ 0:4
        println("\n++++++++++++++++++++++++++\n")
        println("\nO-ShapePartitions m = $m")
        for n ∈ 0:6
            #OrderedShapePartitionsOfTypeM(m, n) |> Println
            OrderedShapePartitions(m, n) |> Println
        end
        println()
        seq = OrderedShapePartitions(m, 5)
        oeis_search(seq, 2)
        println()
        println("Row sums:")
        rsum = rowsum(OrderedShapePartitions, m, 6)
        oeis_search(rsum, 2)
        println()
        println("Alternating row sums:")
        asum = altrowsum(OrderedShapePartitions, m, 6)
        oeis_search(asum, 2)
        println()

        println("\nO-SetPartitions m = $m")
        for n ∈ 0:6
            #OrderedSetPartitionsOfTypeM(m, n) |> Println
            OrderedSetPartitions(m, n) |> Println
        end
        println()
        seq = OrderedSetPartitions(m, 5)
        oeis_search(seq, 2)
        println()
        println("Row sums:")
        rsum = rowsum(OrderedSetPartitions, m, 6)
        oeis_search(rsum, 2)
        println()
        println("Alternating row sums:")
        asum = altrowsum(OrderedSetPartitions, m, 6)
        oeis_search(asum, 2)
        println()
    end

end

Test()

end # module

#=

('--', 0, '---------------------')

by Shape    A000012
[0] [1]
[1] [1]
[2] [1, 1]
[3] [1, 1, 1]
[4] [1, 1, 1, 1, 1]
[5] [1, 1, 1, 1, 1, 1, 1]
[6] [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

by Length   A008284
[0] [1]
[1] [0, 1]
[2] [0, 1, 1]
[3] [0, 1, 1, 1]
[4] [0, 1, 2, 1, 1]
[5] [0, 1, 2, 2, 1, 1]
[6] [0, 1, 3, 3, 2, 1, 1]

[1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 2, 1, 1, 0, 1, 2, 2, 1, 1]

RowSum      A000041 [1, 1, 2, 3, 5, 7, 11]
AlterRowSum A081362 [1, -1, 0, -1, 1, -1, 1]
Diagonal    A000012 [1, 1, 1, 1, 1, 1, 1]
Central     A000041 [1, 1, 2, 3, 5, 7]

('--', 1, '---------------------')

by Shape    A036040
[0] [1]
[1] [1]
[2] [1, 1]
[3] [1, 3, 1]
[4] [1, 4, 3, 6, 1]
[5] [1, 5, 10, 10, 15, 10, 1]
[6] [1, 6, 15, 10, 15, 60, 15, 20, 45, 15, 1]

[1, 1, 1, 1, 1, 3, 1, 1, 4, 3, 6, 1]

by Length   A048993
[0] [1]
[1] [0, 1]
[2] [0, 1, 1]
[3] [0, 1, 3, 1]
[4] [0, 1, 7, 6, 1]
[5] [0, 1, 15, 25, 10, 1]
[6] [0, 1, 31, 90, 65, 15, 1]

by Recur/k!
[0] [1]
[1] [0, 1]
[2] [0, 1, 1]
[3] [0, 1, 3, 1]
[4] [0, 1, 7, 6, 1]
[5] [0, 1, 15, 25, 10, 1]
[6] [0, 1, 31, 90, 65, 15, 1]

by Recurrence
[0] [1]
[1] [0, 1]
[2] [0, 1, 2]
[3] [0, 1, 6, 6]
[4] [0, 1, 14, 36, 24]
[5] [0, 1, 30, 150, 240, 120]
[6] [0, 1, 62, 540, 1560, 1800, 720]

[1, 0, 1, 0, 1, 1, 0, 1, 3, 1, 0, 1, 7, 6, 1, 0, 1, 15, 25, 10, 1]

RowSum      A000110 [1, 1, 2, 5, 15, 52, 203]
AlterRowSum A000587 [1, -1, 0, 1, 1, -2, -9]
Diagonal    A000012 [1, 1, 1, 1, 1, 1, 1]
Central     A007820 [1, 1, 7, 90, 1701, 42525]

('--', 2, '---------------------')

by Shape    A257490
[0] [1]
[1] [1]
[2] [1, 3]
[3] [1, 15, 15]
[4] [1, 28, 35, 210, 105]
[5] [1, 45, 210, 630, 1575, 3150, 945]
[6] [1, 66, 495, 462, 1485, 13860, 5775, 13860, 51975, 51975, 10395]

[1, 1, 1, 3, 1, 15, 15, 1, 28, 35, 210, 105]

by Length   A156289
[0] [1]
[1] [0, 1]
[2] [0, 1, 3]
[3] [0, 1, 15, 15]
[4] [0, 1, 63, 210, 105]
[5] [0, 1, 255, 2205, 3150, 945]
[6] [0, 1, 1023, 21120, 65835, 51975, 10395]

by Recur/k!
[0] [1]
[1] [0, 1]
[2] [0, 1, 3]
[3] [0, 1, 15, 15]
[4] [0, 1, 63, 210, 105]
[5] [0, 1, 255, 2205, 3150, 945]
[6] [0, 1, 1023, 21120, 65835, 51975, 10395]

by Recurrence
[0] [1]
[1] [0, 1]
[2] [0, 1, 6]
[3] [0, 1, 30, 90]
[4] [0, 1, 126, 1260, 2520]
[5] [0, 1, 510, 13230, 75600, 113400]
[6] [0, 1, 2046, 126720, 1580040, 6237000, 7484400]

[1, 0, 1, 0, 1, 3, 0, 1, 15, 15, 0, 1, 63, 210, 105, 0, 1, 255, 2205, 3150, 945]

RowSum      A005046 [1, 1, 4, 31, 379, 6556, 150349]
AlterRowSum A260884 [1, -1, 2, -1, -43, 254, 4157]
Diagonal    A001147 [1, 1, 3, 15, 105, 945, 10395]
Central     A327416 [1, 1, 63, 21120, 20585565, 44025570225]

('--', 3, '---------------------')

by Shape    A327003
[0] [1]
[1] [1]
[2] [1, 10]
[3] [1, 84, 280]
[4] [1, 220, 462, 9240, 15400]
[5] [1, 455, 5005, 50050, 210210, 1401400, 1401400]
[6] [1, 816, 18564, 24310, 185640, 4084080, 2858856, 13613600, 85765680, 285885600, 190590400]

[1, 1, 1, 10, 1, 84, 280, 1, 220, 462, 9240, 15400]

by Length   A291451
[0] [1]
[1] [0, 1]
[2] [0, 1, 10]
[3] [0, 1, 84, 280]
[4] [0, 1, 682, 9240, 15400]
[5] [0, 1, 5460, 260260, 1401400, 1401400]
[6] [0, 1, 43690, 7128576, 99379280, 285885600, 190590400]

by Recur/k!
[0] [1]
[1] [0, 1]
[2] [0, 1, 10]
[3] [0, 1, 84, 280]
[4] [0, 1, 682, 9240, 15400]
[5] [0, 1, 5460, 260260, 1401400, 1401400]
[6] [0, 1, 43690, 7128576, 99379280, 285885600, 190590400]

by Recurrence
[0] [1]
[1] [0, 1]
[2] [0, 1, 20]
[3] [0, 1, 168, 1680]
[4] [0, 1, 1364, 55440, 369600]
[5] [0, 1, 10920, 1561560, 33633600, 168168000]
[6] [0, 1, 87380, 42771456, 2385102720, 34306272000, 137225088000]

[1, 0, 1, 0, 1, 10, 0, 1, 84, 280, 0, 1, 682, 9240, 15400, 0, 1, 5460, 260260, 1401400, 1401400]

RowSum      A291973 [1, 1, 11, 365, 25323, 3068521, 583027547]
AlterRowSum A291974 [1, -1, 9, -197, 6841, -254801, -3000807]
Diagonal    A025035 [1, 1, 10, 280, 15400, 1401400, 190590400]
Central     A327417 [1, 1, 682, 7128576, 429120851544, 94066556834970720]

('--', 4, '---------------------')

by Shape    A327004
[0] [1]
[1] [1]
[2] [1, 35]
[3] [1, 495, 5775]
[4] [1, 1820, 6435, 450450, 2627625]
[5] [1, 4845, 125970, 4408950, 31177575, 727476750, 2546168625]
[6] [1, 10626, 735471, 1352078, 25741485, 1338557220, 1577585295, 15616500900, 165646455975, 1932541986375, 4509264634875]

[1, 1, 1, 35, 1, 495, 5775, 1, 1820, 6435, 450450, 2627625]

by Length   A291452
[0] [1]
[1] [0, 1]
[2] [0, 1, 35]
[3] [0, 1, 495, 5775]
[4] [0, 1, 8255, 450450, 2627625]
[5] [0, 1, 130815, 35586525, 727476750, 2546168625]
[6] [0, 1, 2098175, 2941884000, 181262956875, 1932541986375, 4509264634875]

by Recur/k!
[0] [1]
[1] [0, 1]
[2] [0, 1, 35]
[3] [0, 1, 495, 5775]
[4] [0, 1, 8255, 450450, 2627625]
[5] [0, 1, 130815, 35586525, 727476750, 2546168625]
[6] [0, 1, 2098175, 2941884000, 181262956875, 1932541986375, 4509264634875]

by Recurrence
[0] [1]
[1] [0, 1]
[2] [0, 1, 70]
[3] [0, 1, 990, 34650]
[4] [0, 1, 16510, 2702700, 63063000]
[5] [0, 1, 261630, 213519150, 17459442000, 305540235000]
[6] [0, 1, 4196350, 17651304000, 4350310965000, 231905038365000, 3246670537110000]

[1, 0, 1, 0, 1, 35, 0, 1, 495, 5775, 0, 1, 8255, 450450, 2627625, 0, 1, 130815, 35586525, 727476750, 2546168625]

RowSum      A291975 [1, 1, 36, 6271, 3086331, 3309362716, 6626013560301]
AlterRowSum A291976 [1, -1, 34, -5281, 2185429, -1854147586, 2755045819549]
Diagonal    A025036 [1, 1, 35, 5775, 2627625, 2546168625, 4509264634875]
Central     A327418 [1, 1, 8255, 2941884000, 11957867341948125, 294040106448733743008625]

=#

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
function POMDPs.solve{S,A}(solver::RobustMCTSSolver, rmdp::RobustMDP{S,A})
    if isa(solver.rollout_solver, Solver)
        rollout_policy = solve(solver.rollout_solver, solver.rollout_nature)
    else
        rollout_policy = solver.rollout_solver
    end
    return RobustMCTSPlanner(solver, rmdp, Dict{S,RobustStateNode{S,A}}(), rollout_policy)
end

function POMDPs.action{S,A}(p::RobustMCTSPlanner{S,A}, s::S, a::A=create_action(p.rmdp))
    # This function calls simulate and chooses the approximate best action from the reward approximations
    # XXX do we need to make a copy of the state here?
    if p.solver.reset_tree_at_every_decision
        p.tree = Dict{S,RobustStateNode{S,A}}()
    end
    for i = 1:p.solver.n_iterations
        simulate(p, deepcopy(s), p.solver.depth)
    end
    snode = p.tree[s]
    best_Q = -Inf
    local best_a
    for (a, sanode) in snode.A
        if sanode.Q > best_Q
            best_Q = sanode.Q
            best_a = a
        end
    end
    # XXX some publications say to choose action that has been visited the most
    return best_a # choose action with highest approximate value
end

function simulate{S,A}(p::RobustMCTSPlanner{S,A}, s::S, d::Int)
    # TODO: reimplement this as a loop instead of a recursion?

    # This function returns the reward for one iteration of MCTS
    if d == 0 || isterminal(p.rmdp, s)
        return 0.0 # XXX is this right or should it be a rollout?
    end
    if !haskey(p.tree,s) # if state is not yet explored, add it to the set of states, perform a rollout
        p.tree[s] = RobustStateNode{S,A}() # TODO: Mechanism to set N0
        p.tree[s].N += 1
        return estimate_value(p,s,d)
    end

    snode = p.tree[s] # save current state node so we do not have to iterate through map many times
    snode.N += 1

    # action progressive widening
    if length(snode.A) <= p.solver.k_action*snode.N^p.solver.alpha_action # criterion for new action generation
        a = next_action(p.solver.action_generator, p.rmdp, s, snode) # action generation step
        if !haskey(snode.A,a) # make sure we haven't already tried this action
            snode.A[a] = RobustActionNode{S}() # TODO: Mechanism to set N0, Q0
        end
    end

    best_UCB = -Inf
    local a
    sN = snode.N
    for (act, sanode) in snode.A
        if sN == 1 && sanode.N == 0
            UCB = sanode.Q
        else
            c = p.solver.c # for clarity
            UCB = sanode.Q + c*sqrt(log(sN)/sanode.N)
        end
        @assert !isnan(UCB)
        @assert !isequal(UCB, -Inf)
        if UCB > best_UCB
            best_UCB = UCB
            a = act
        end
    end

    sanode = snode.A[a]
    sanode.N += 1

    # nature's progressive widening
    if length(sanode.nature_nodes) <= p.solver.k_nature*sanode.N^p.solver.alpha_nature
        tau = next_model(p.solver.model_generator, p.rmdp, s, a, snode, sanode)
        if !haskey(sanode.nature_nodes, tau)
            sanode.nature_nodes[tau] = NatureNode{S}()
        end
    end

    # choose worst nature
    worst_UCB = Inf
    local tau
    saN = sanode.N
    for (mdp, nnode) in sanode.nature_nodes
        if saN == 1 && nnode.N == 0
            UCB = nnode.Q
        else
            c = p.solver.c_nature # for clarity
            UCB = nnode.Q - c*sqrt(log(saN)/nnode.N)
        end
        @assert !isnan(UCB)
        @assert !isequal(UCB, Inf)
        if UCB < worst_UCB
            worst_UCB = UCB
            tau = mdp
        end
    end

    nnode = sanode.nature_nodes[tau]

    # state progressive widening
    if length(nnode.transitions) <= p.solver.k_state*nnode.N^p.solver.alpha_state # criterion for new transition state consideration
        sp, r = generate_sr(tau, s, a, p.solver.rng) # choose a new state and get reward

        if !haskey(nnode.transitions,sp) # if transition state not yet explored, add to set and update reward
            nnode.transitions[sp] = TransitionNode() # TODO: mechanism for assigning N0
            nnode.transitions[sp].R = r
        end
        nnode.transitions[sp].N += 1

    else # sample from transition states proportional to their occurence in the past
        # warn("sampling states: |V|=$(length(nnode.V)), N=$(nnode.N)")
        total_N = reduce(add_N, 0, values(nnode.transitions))
        rn = rand(p.solver.rng, 1:total_N) # this is where Jon's bug was (I think)
        cnt = 0
        local sp, tnode
        for (sp,tnode) in nnode.transitions
            cnt += tnode.N
            if rn <= cnt
                break
            end
        end

        r = tnode.R
    end

    q = r + discount(p.rmdp)*simulate(p,sp,d-1)

    nnode.N += 1

    sanode.Q += (q - sanode.Q)/sanode.N
    nnode.Q += (q - nnode.Q)/nnode.N

    return q
end

"""
Add the N's of two sas nodes - for use in reduce
"""
add_N(a::TransitionNode, b::TransitionNode) = a.N + b.N
add_N(a::Int, b::TransitionNode) = a + b.N

"""
Estimate the value of being at state s

By default, run a rollout simulation.
"""
function estimate_value(p::RobustMCTSPlanner, s, d::Int)
    rollout(p, s, d)
end

function rollout(p::RobustMCTSPlanner, s, d::Int)
    sim = RolloutSimulator(rng=p.solver.rng, max_steps=d) # TODO(?) add a mechanism to customize this
    POMDPs.simulate(sim, p.solver.rollout_nature, p.rollout_policy, s)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule alsa_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("alsa")
JLLWrappers.@generate_main_file("alsa", UUID("45378030-f8ea-5b20-a7c7-1a9d95efb90e"))
end  # module alsa_jll


text:
"""
Represent a string in terms of the counts of the ASCII lowercase letters in it.
"""
struct LetterCounts
    counts::MVector{26, Int}

    function LetterCounts()
        counts = fill(0, MVector{26, Int})
        new(counts)
    end
end

"""
Read from LetterCounts using `letter_counts['a']`.
"""
function Base.getindex(letter_counts::LetterCounts, c::Char)
    letter_counts.counts[to_index(c)]
end

"""
Write to LetterCounts using `letter_counts['a'] = 3`.
"""
function Base.setindex!(letter_counts::LetterCounts, cnt::Int, c::Char)
    letter_counts.counts[to_index(c)] = cnt
    cnt
end


text:
"""
Functions to estimate Q|_{U^0}. See our paper for details.
"""

using LinearAlgebra
using SparseArrays
using FastRounding
using ValidatedNumerics

export norms_of_powers, refine_norms_of_powers, norms_of_powers_dfly, norms_of_powers_trivial, norms_of_powers_from_coarser_grid

"""
Returns the maximum number of (structural) nonzeros in a row of A
"""
function max_nonzeros_per_row(A::SparseMatrixCSC)
    rows = rowvals(A)
    m, n = size(A)
    nonzeros_in_each_row = zeros(eltype(rows), m)
    for i in rows
        nonzeros_in_each_row[i] += 1
    end
    return maximum(nonzeros_in_each_row)
end

"""
γₙ constants for floating point error estimation, as in [Higham, Accuracy and Stability of Numerical Algorithms]
"""
function gamma(T, n::Integer)
    u = eps(T)
    nu = u ⊗₊ T(n) #TODO: in theory, this should be rounded up/down. In practice, all integers up to 2^53 or so fit in Float64, so it won't be needed.
    return nu ⊘₊ (one(T) ⊖₋ nu)
end

using ProgressMeter
"""
Estimates the norms ||Q||, ||Q^2||, ... ||Q^m|| on U^0.

U is the matrix [ones(1,n-1); -I_(n-1,n-1)]. It is currently assumed that
f*U==0 (i.e., all elements of f are equal).

f must be an interval vector.

The following constants may be specified as keyword arguments:

normQ, normE, normv0, normEF, normIEF, normN

otherwise they are computed (which may be slower).

e and f must be specified in case is_integral_preserving==false
In case is_integral_preserving is true, they may be specified but they are then ignored.

Implementation note: currently we perform this computation one column at a time,
to be able to scale (slowly) to cases with large size; for moderate sizes, it would
indeed be better to do the computation all columns at the same time, in BLAS level 3.
"""
function norms_of_powers(N::Type{<:NormKind}, m::Integer, Q::DiscretizedOperator, f::AbstractArray;
        normv0::Real=-1., #used as "missing" value
        normQ::Real=-1.,
        normE::Real=-1.,
        normEF::Real=-1.,
        normIEF::Real=-1.,
        normN::Real=-1.)

    @assert eltype(f) <: Interval
    T = typeof(zero(eltype(Q.L)).hi) # gets "Float64" from Q.L
    n = size(Q.L, 1)
    M = mid.(Q.L)
    R = radius.(Q.L)
    δ = opnormbound(N, R)
    γz = gamma(T, max_nonzeros_per_row(Q.L))
    γn = gamma(T, n+3) # not n+2 like in the paper, because we wish to allow for f to be the result of rounding
    ϵ = zero(T)

    nrmM = opnormbound(N, M)

    # precompute norms
    if !is_integral_preserving(Q)
        if normE == -1.
            normE = opnormbound(N, Q.e)
        end
        if normEF == -1.
            normEF = opnormbound(N, Q.e*f)
        end
        if normIEF == -1.
            normIEF =  opnormbound(N, [Matrix(UniformScaling{Float64}(1),n,n) Q.e*f])
        end
        if normN == -1.
            normN = opnormbound(N, Matrix(UniformScaling{Float64}(1),n,n) - Q.e*f)
        end
    end

    if normQ == -1.
        if is_integral_preserving(Q)
            normQ = nrmM ⊕₊ δ
        else
            defect = opnormbound(N, Q.w)
            normQ = nrmM ⊕₊ δ ⊕₊ normE ⊗₊ defect
        end
    end

    # initialize normcachers
    normcachers = [NormCacher{N}(n) for j in 1:m]

    midf = map(mid, f)

    # main loop

    v = zeros(T, n)
    @showprogress for j = 1:n-1
        v .= zero(T) # TODO: check for type stability in cases with unusual types
        v[1] = one(T) # TODO: in full generality, this should contain entries of f rather than ±1
        v[j+1] = -one(T)
        if normv0 == -1.
            nrmv = opnormbound(N, v)
        else
            nrmv = normv0
        end
        ϵ = 0.
        nrmw = nrmv # we assume that initial vectors are already integral-preserving
        for k = 1:m
            w = M * v
            if is_integral_preserving(Q)
                v = w
                ϵ = (γz ⊗₊ nrmM ⊕₊ δ) ⊗₊ nrmv ⊕₊ normQ ⊗₊ ϵ
            else
                v = w - Q.e * (midf*w)  # TODO: we are currently assuming that f is not too large, to estimate the error (the result of only one floating point operation)
                new_nrmw = opnormbound(N, w)
                ϵ = γn ⊗₊ normIEF ⊗₊ (new_nrmw ⊕₊ normEF ⊗₊ nrmw) ⊕₊ normN ⊗₊ (γz ⊗₊ nrmM ⊕₊ δ) ⊗₊ nrmv ⊕₊ normQ ⊗₊ ϵ
                nrmw = new_nrmw
            end
            nrmv = opnormbound(N, v)
            add_column!(normcachers[k], v, ϵ) #TODO: Could pass and reuse nrmv in the case of norm-1
        end
    end
    return map(get_norm, normcachers)
end

"""
Array of "trivial" bounds for the powers of a DiscretizedOperator (on the whole space)
coming from from ||Q^k|| ≤ ||Q||^k
"""
function norms_of_powers_trivial(N::Type{<:NormKind}, Q::DiscretizedOperator, m::Integer)
    norms = fill(NaN, m)
    norms[1] = opnormbound(N, Q)
    for i = 2:m
        norms[i] = norms[i-1] ⊗₊ norms[1]
    end
    return norms
end

"""
Arrays of bounds to ||Q^k||_{w → s} = sup_{||f||_w=1} ||Q^k f||_s
and to ||Q^k||_{w}
coming theoretically from iterated DFLY inequalities (the "small matrix method").

Returns two arrays (strongs, norms) of length m:
strongs[k] bounds ||Q^k f||_s, norms[k] bounds ||Q^k f||)
"""
function norms_of_powers_dfly(Bas::Basis, D::Dynamic, m)
    A, B = dfly(strong_norm(Bas), aux_norm(Bas), D)
    Eh = BasisDefinition.aux_normalized_projection_error(Bas)
    M₁n = BasisDefinition.strong_weak_bound(Bas)
    M₂ = BasisDefinition.aux_weak_bound(Bas)
    S₁, S₂ = BasisDefinition.weak_by_strong_and_aux_bound(Bas)

    norms = fill(NaN, m)
    strongs = fill(NaN, m)

    v = Array{Float64}([M₁n; M₂])
    # We evaluate [S₁ S₂] * ([1 0; Eh 1]*[A B; 0 1])^k * [M₁n; M₂] (with correct rounding)
    for k = 1:m
        # invariant: v[1] bounds ||Q^kf||_s for ||f||_w=1
        # v[2] bounds |||Q^kf||| for ||f||_w=1
        v[1] = A ⊗₊ v[1] ⊕₊ B ⊗₊ v[2]
        v[2] = Eh ⊗₊ v[1] ⊕₊ v[2]
        strongs[k] = v[1]
        norms[k] = S₁ ⊗₊ v[1] ⊕₊ S₂ ⊗₊ v[2]
    end
    return strongs, norms
end

"""
Compute better and/or more estimates of power norms using
the fact that ||Q^{k+h}|| ≤ ||Q^k|| * ||Q^h||.
This uses multiplicativity, so it will not work for mixed norms,
e.g., ||Q^k||_{s → w}, or ||M^k|_{U^0}||
(unless M preserves U^0, which is the case for Q|_{U^0}).
"""
function refine_norms_of_powers(norms::Vector, m)
    better_norms = fill(NaN, m)
    better_norms[1] = norms[1]
    for k in 2:m
        better_norms[k] = minimum(better_norms[i] ⊗₊ better_norms[k-i] for i in 1:k-1)
        if k <= length(norms)
            better_norms[k] = min(better_norms[k], norms[k])
        end
    end
    return better_norms
end
refine_norms_of_powers(norms::Vector) = refine_norms_of_powers(norms, length(norms))

"""
Estimate norms of powers from those on a coarser grid (see paper for details)
"""
function norms_of_powers_from_coarser_grid(fine_basis::Basis, coarse_basis::Basis, D::Dynamic, coarse_norms::Vector, normQ::Real)
    if !BasisDefinition.is_refinement(fine_basis, coarse_basis)
        @error "The fine basis is not a refinement of the coarse basis"
    end
    m = length(coarse_norms)
    fine_norms = fill(NaN, m)
    (strongs, norms) = norms_of_powers_dfly(fine_basis, D, m)

    # adds a 0th element to strongs
    strongs0(k::Integer) = k==0 ? BasisDefinition.strong_weak_bound(fine_basis) : strongs[k]
    coarse_norms0(k::Integer) = k==0 ? 1. : coarse_norms[k]

    Kh =  BasisDefinition.weak_projection_error(coarse_basis)
    for k in 1:m
		temp = 0.
		for k in 0:m-1
			temp = temp ⊕₊ coarse_norms0(m-1-k) ⊗₊ (normQ ⊗₊ strongs0(k) ⊕₊ strongs0(k+1))
		end
		fine_norms[k] = coarse_norms[k] ⊕₊ 2. ⊗₊ Kh ⊗₊ temp
	end
    return fine_norms
end

"""
Estimate ||I + Q + Q^2 + … || (infinite sum) using a list of computed norm bounds norms[k] ≥ ||Q^k||.
"""
function infinite_sum_norms(norms::Vector)
    m = length(norms)
    if norms[m] >= 1
        @error "The last norm is >= 1: the bounds are not sufficient to show that the sum converges"
    end
    S = 1.
    for k = 1:m-1
        S = S ⊕₊ norms[k]
    end
    S = S ⊘₊ (1. ⊖₋ norms[m])
    return S
end


text:
#=
 Train a simple network to predict formation energy per atom (downloaded from Materials Project).
=#
using CSV, DataFrames
using Random, Statistics
using Flux
using Flux: @epochs
using ChemistryFeaturization
using AtomicGraphNets
using Serialization

cd(@__DIR__)

println("Setting things up...")

# where to find the inputs
csv_path = "qm9.csv"
xyz_dir = "xyz/"
graph_dir = "graphs/"

# data-related options
num_pts = 111
train_frac = 0.8 # what fraction for training?
num_epochs = 40 # how many epochs to train?
num_train = Int32(round(train_frac * num_pts))
num_test = num_pts - num_train
prop = :u0 # internal energy at 0K
info = DataFrame!(CSV.File(csv_path))
y = Array(Float32.(info[!, Symbol(prop)]))

# atom featurization, pretty arbitrary choices for now
features = Symbol.(["Group", "Row", "Block", "Atomic mass", "Atomic radius", "X"])
num_bins = [18, 9, 4, 16, 10, 10]
num_features = sum(num_bins) # we'll use this later
logspaced = [false, false, false, true, true, false]
# returns actual vectors (in a dict with keys of elements) plus Vector of AtomFeat objects describing featurization metadata
atom_feature_vecs, featurization = make_feature_vectors(features, nbins=num_bins, logspaced=logspaced)

# model hyperparameters – keeping it pretty simple for now
num_conv = 5 # how many convolutional layers?
atom_fea_len = 80
crys_fea_len = 32 # length of crystal feature vector after pooling (keep node dimension constant for now)
pool_type = "max"
num_hidden_layers = 2 # how many fully-connected layers after convolution and pooling?
opt = ADAM(0.003) # optimizer

# shuffle data and pick out subset
indices = shuffle(1:size(info,1))[1:num_pts]
info = info[indices,:]
output = y[indices]

# next, make graphs and build input features (matrices of dimension (# features, # nodes))
println("Building graphs and feature vectors from structures...")
inputs = AtomGraph[]

# build the graphs
build_graphs_batch(xyz_dir; output_folder = graph_dir)

for r in eachrow(info)
    graph_path = joinpath(graph_dir, string(r.mol_id, ".jls"))
    gr = deserialize(graph_path)
    feature_mat = hcat([atom_feature_vecs[e] for e in gr.elements]...)
    add_features!(gr, feature_mat, featurization)
    push!(inputs, gr)
end

# pick out train/test sets
println("Dividing into train/test sets...")
train_output = output[1:num_train]
test_output = output[num_train+1:end]
train_input = inputs[1:num_train]
test_input = inputs[num_train+1:end]
train_data = zip(train_input, train_output)

# build the model
println("Building the network...")
model = Xie_model(num_features, num_conv=num_conv, atom_conv_feature_length=atom_fea_len, pool_type=pool_type, pooled_feature_length=crys_fea_len, num_hidden_layers=num_hidden_layers)

# define loss function
loss(x,y) = Flux.Losses.mse(model(x), y)
# and a callback to see training progress
evalcb() = @show(mean(loss.(test_input, test_output)))
evalcb()

# train
println("Training!")
@epochs num_epochs Flux.train!(loss, params(model), train_data, opt, cb = Flux.throttle(evalcb, 10))


text:
# ---
# title: 722. Remove Comments
# id: problem722
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/remove-comments/description/>
# hidden: true
# ---
#
# Given a C++ program, remove comments from it. The program `source` is an array
# where `source[i]` is the `i`-th line of the source code. This represents the
# result of splitting the original source code string by the newline character
# `\n`.
#
# In C++, there are two types of comments, line comments, and block comments.
#
# The string `//` denotes a line comment, which represents that it and rest of
# the characters to the right of it in the same line should be ignored.
#
# The string `/*` denotes a block comment, which represents that all characters
# until the next (non-overlapping) occurrence of `*/` should be ignored. (Here,
# occurrences happen in reading order: line by line from left to right.) To be
# clear, the string `/*/` does not yet end the block comment, as the ending
# would be overlapping the beginning.
#
# The first effective comment takes precedence over others: if the string `//`
# occurs in a block comment, it is ignored. Similarly, if the string `/*` occurs
# in a line or block comment, it is also ignored.
#
# If a certain line of code is empty after removing comments, you must not
# output that line: each string in the answer list will be non-empty.
#
# There will be no control characters, single quote, or double quote characters.
# For example, `source = "string s = "/* Not a comment. */";"` will not be a
# test case. (Also, nothing else such as defines or macros will interfere with
# the comments.)
#
# It is guaranteed that every open block comment will eventually be closed, so
# `/*` outside of a line or block comment always starts a new comment.
#
# Finally, implicit newline characters can be deleted by block comments. Please
# see the examples below for details.
#
# After removing the comments from the source code, return the source code in
# the same format.
#
# **Example 1:**
#
#
#
#     Input:
#     source = ["/*Test program */", "int main()", "{ ", "  // variable declaration ", "int a, b, c;", "/* This is a test", "   multiline  ", "   comment for ", "   testing */", "a = b + c;", "}"]
#
#     The line by line code is visualized as below:
#     /*Test program */
#     int main()
#     {
#       // variable declaration
#     int a, b, c;
#     /* This is a test
#        multiline
#        comment for
#        testing */
#     a = b + c;
#     }
#
#     Output: ["int main()","{ ","  ","int a, b, c;","a = b + c;","}"]
#
#     The line by line code is visualized as below:
#     int main()
#     {
#
#     int a, b, c;
#     a = b + c;
#     }
#
#     Explanation:
#     The string /* denotes a block comment, including line 1 and lines 6-9. The string // denotes line 4 as comments.
#
#
# **Example 2:**
#
#
#
#     Input:
#     source = ["a/*comment", "line", "more_comment*/b"]
#     Output: ["ab"]
#     Explanation: The original source string is "a/*comment **\n** line **\n** more_comment*/b", where we have bolded the newline characters.  After deletion, the _implicit_ newline characters are deleted, leaving the string "ab", which when delimited by newline characters becomes ["ab"].
#
#
# **Note:**
#
# * The length of `source` is in the range `[1, 100]`.
# * The length of `source[i]` is in the range `[0, 80]`.
# * Every open block comment is eventually closed.
# * There are no single-quote, double-quote, or control characters in the source code.
#
#
## @lc code=start
using LeetCode

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


text:
# data structure for a story
struct Story
    id::Int64
    creator_username::String
    creation_time::Int64
    title::String
    score::Int64
    comment_count::Int64
    scrape_timestamp::Float64
end

# convenience constructor with named parameters
Story(;
    id::Int64, creator_username::String, creation_time::Int64, title::String, score::Int64, comment_count::Int64, scrape_timestamp::Float64
) = Story(
    id, creator_username, creation_time, title, score, comment_count, scrape_timestamp
)

# default constructor for convenient filtering of faulty stories
Story() = Story(id = -1, creator_username = "NA", creation_time = -1, title = "NA", score = -1, comment_count = -1, scrape_timestamp = -1)

# data structure for a topstories list
struct TopStoriesList
    itemlist::Array{Story}
    slice_timestamp::Float64
end

# get top stories from hacker news (every `interval_minutes` for a total of `n_timeslices` iterations)
function scrape_topstories(; n_timeslices::Int, interval_minutes::Int, base_filename::String)
    (interval_minutes >= 1) || throw(DomainError(interval_minutes, "sampling rate cannot be less than 1 minute"))
    if !(".tmp" in readdir())
        mkdir(".tmp")
    end
    df_list = []
    n = 0
    for i in 1:n_timeslices
        t = Timer(interval_minutes * 60)
        current_time = Dates.now()
        @info "Request at $current_time"
        ts_df = log_topstories()
        push!(df_list, deepcopy(ts_df))
        if ((i % 60) == 0) || (i == n_timeslices)
            n += 1
            filename = base_filename * "_$n" * ".feather"
            write_tmp_stories(df_list, filename)
            df_list = []
        end
        @info "Logged time slice. Waiting for next request."
        wait(t)
    end
    df = aggregate_data()
    rm(".tmp", recursive = true)
    return df
end

function write_tmp_stories(df_list, filename)
    topstories = reduce(vcat, df_list)
    Feather.write(joinpath(".tmp", filename), topstories)
    @info "Saved current data and emptied cache."
    return nothing
end

function aggregate_data()
    filenames = joinpath.(".tmp", readdir(".tmp"))
    df_list = asyncmap(Feather.read, filenames, ntasks = 20)
    df = reduce(vcat, df_list)
    return df
end

# log topstories with timestamp of extraction time
function log_topstories()
    topstories = get_topstories()
    topstories_df = DataFrames.DataFrame(topstories.itemlist)
    topstories_df[!, :slice_timestamp] .= topstories.slice_timestamp
    topstories_df[!, :rank_at_timestamp] = 1:length(topstories)
    return topstories_df
end

Base.length(tsl::TopStoriesList) = length(tsl.itemlist)

# get topstories
function get_topstories()
    identifier = "https://hacker-news.firebaseio.com/v0/topstories.json"
    slice_timestamp = Dates.datetime2unix(Dates.now(UTC))
    response = HTTP.get(identifier)
    response_body = String(response.body)
    topstories_id_list = JSON.parse(response_body)
    topstories = asyncmap(get_story, topstories_id_list, ntasks = 20)
    return TopStoriesList(topstories, slice_timestamp)
end

# get story from id (factory-like)
function get_story(story_id::Int)
    identifier = "https://hacker-news.firebaseio.com/v0/item/$story_id" * ".json"
    scrape_timestamp = Dates.datetime2unix(Dates.now(UTC))
    response = HTTP.get(identifier)
    response_body = String(response.body)
    story = JSON.parse(response_body)
    element_keys = keys(story)
    return Story(
        id = "id" in element_keys ? story["id"] : -1,
        creator_username = "by" in element_keys ? story["by"] : "NA",
        creation_time = "time" in element_keys ? story["time"] : -1,
        title = "title" in element_keys ? story["title"] : "NA",
        score = "score" in element_keys ? story["score"] : -1,
        comment_count = "descendants" in element_keys ? story["descendants"] : -1,
        scrape_timestamp = scrape_timestamp
    )
end

# get new stories (IDs) from hacker news (every `interval_minutes` for a total of `n_timeslices` iterations)
function scrape_newstories(; n_timeslices::Int, interval_minutes::Int)
    (interval_minutes >= 1) || throw(DomainError(interval_minutes, "sampling rate cannot be less than 1 minute"))
    df_list = []
    identifier = "https://hacker-news.firebaseio.com/v0/newstories.json"
    for i in 1:n_timeslices
        t = Timer(interval_minutes * 60)
        timestamp = Dates.now(UTC)
        response = HTTP.get(identifier)
        response_body = String(response.body)
        newstories = Array{Int64}(JSON.parse(response_body))
        newstories_df = DataFrames.DataFrame(
            story_id = newstories, timestamp = Dates.datetime2unix(timestamp), rank = 1:length(newstories)
        )
        # newstories_df[!, :timestamp] .= timestamp

        push!(df_list, deepcopy(newstories_df))
        @info "Logged new stories at $timestamp"
        wait(t)
    end
    newstories = reduce(vcat, df_list)
    return newstories
end

# POTENTIAL TO DO: subscribe to updates

max_stars_repo_name: JohannesNakayama/hacker-news-analysis

max_issues_repo_name: JohannesNakayama/hacker-news-analysis

max_forks_repo_name: JohannesNakayama/hacker-news-analysis


text:
#=
# Support
GitHub issues will be used for tracking requests and bugs. For questions go to [forums.xilinx.com](http://forums.xilinx.com/).

# License
Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License.
You may obtain a copy of the License at [http://www.apache.org/licenses/LICENSE-2.0]( http://www.apache.org/licenses/LICENSE-2.0 )

Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.

XD0xx | © Copyright 2021 Xilinx, Inc.

=#

# Ref: https://docs.juliadsp.org/latest/filters/#Filters-filter-design-and-filtering

# coefficient matrix written out as 1d array

#using DSP
#using PyPlot
#using DoubleFloats
#using Printf

using DSP
using FFTW
using Noise

# --- begin user parameters

f1 = 2.0e6              # frequency of s1
f2 = 18.0e6             # frequency of s2
fs = 100.0e6            # sampling frequency
N = 256                 # no. of samples for impulse response
a1 = 0.5                # amplitude of s1
a2 = 0.5                # amplitude of s2

# --- end user parameters

n = 0 : (N - 1)         # sampling index
Ts = 1/fs               # sampling period

s1 = a1 * sin.(2 * pi * f1 * n * Ts)
s2 = a2 * sin.(2 * pi * f2 * n * Ts)
s = add_gauss(s1 + s2, 0.2)

xfactor = 16
s_pad = zeros(N * xfactor)
s_pad[1:N] = s          # zero padding for FFT

win = kaiser(N * xfactor, 3, padding=0, zerophase=false)

S = (fft(win .* s_pad)) |> fftshift
freqs = fftfreq(length(s_pad), fs) |> fftshift
PS = 20*log10.(abs.(S)) # power spectrum

close("all")

fig1, ax1 = subplots()
ax1.plot(n, s)
ax1.grid("on")
ax1.set_title("Noisy Input Signal")
ax1.set_xlabel("Sampling Index")
ax1.set_ylabel("Amplitude")
ax1.set_xlim([0, N])
ax1.set_ylim([-1.5, 1.5])

fig2, ax2 = subplots()
ax2.plot(freqs[Int32(N * xfactor/2) : end], PS[Int32(N * xfactor/2) : end])
ax2.grid("on")
ax2.set_title("Power Spectrum")
ax2.set_xlabel("Frequency (Hz)")
ax2.set_ylabel("Amplitude (dB)")
ax2.set_xlim([0, fs/2])
ax2.set_ylim([-60, 0])

fid = open("two_freqs.dat", "w")
for i = 1 : N
    @printf(fid, "%f\n", s[i])
end # for i = 1 : N
close(fid)

max_stars_repo_name: jlamperez/Vitis-Tutorials

max_issues_repo_name: jlamperez/Vitis-Tutorials

max_forks_repo_name: jlamperez/Vitis-Tutorials


text:
using Distributed
addprocs(15)
@everywhere include("rijke_tangent_state_estimation.jl")
filenames = Array{String,1}(undef, nworkers())
for i = 1:nworkers()
	filenames[i] = string("../data/rijke_exp9_", string(i), ".jld")
end
pmap(assimilate_parameter_and_trajectory, filenames)

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies