text:
@testset "Constructs_test" begin

    mt_a = Mod("a")
    mt_b = Mod("b")
    mt_c = Mod("c")
    mt_d = Mod("d")
    mt_e = Mod("e")
    # some test constructs
    oc1 = Ordered_Construct([mt_a,mt_d,mt_b])
    oc2 = Ordered_Construct([mt_a,mt_b,mt_e,mt_d])
    unoc1 = Unordered_Construct([mt_e,mt_d,mt_c])
    unoc2 = Unordered_Construct([mt_e,mt_a,mt_d,mt_e])
     @testset "Base" begin
        @test length(oc1) == 3
        @test length(oc2) == 4
        @test length(unoc1) == 3
        @test length(unoc2) == 4

        @test oc1[2] == mt_d
        @test oc2[3] == mt_e

        @test isequal(oc1[1:2],Ordered_Construct([mt_a,mt_d])) == true
        @test isequal(oc1[1:2],Ordered_Construct([mt_a,mt_d])) == true
        @test length(collect(oc1)) == length(oc1)
        @test length(collect(unoc1)) == length(unoc1)
        @test eltype(collect(oc1)) == typeof(mt_a)
    end

    @testset "*_ordered" begin
        @test isequal(mt_a*mt_b,Ordered_Construct([mt_a,mt_b])) == true
        @test isequal(Ordered_Construct([mt_a,mt_b])*mt_c,Ordered_Construct([mt_a,mt_b,mt_c])) == true
        @test isequal(mt_c*Ordered_Construct([mt_a,mt_b]),Ordered_Construct([mt_c,mt_a,mt_b])) == true
        @test isequal((Ordered_Construct([mt_a,mt_b])*Ordered_Construct([mt_d,mt_e])),Ordered_Construct([mt_a,mt_b,mt_d,mt_e])) == true
    end

    @testset "+_unordered" begin
        @test isequal((mt_a+mt_b),Unordered_Construct([mt_a,mt_b])) == true
        @test isequal(Unordered_Construct([mt_a,mt_b])+mt_c,Unordered_Construct([mt_a,mt_b,mt_c])) == true
        @test isequal(mt_c+Unordered_Construct([mt_a,mt_b]),Unordered_Construct([mt_c,mt_a,mt_b])) == true
        @test isequal((Unordered_Construct([mt_a,mt_b])+Unordered_Construct([mt_d,mt_e])),Unordered_Construct([mt_a,mt_b,mt_d,mt_e])) == true
    end

end


text:
l = @lens _[begin]
@test l isa Setfield.DynamicIndexLens
obj = (1,2,3)
@test get(obj, l) == 1
@test set(obj, l, true) == (true,2,3)

l = @lens _[2*begin]
@test l isa Setfield.DynamicIndexLens
obj = (1,2,3)
@test get(obj, l) == 2
@test set(obj, l, true) == (1,true,3)

one = 1
plustwo(x) = x + 2
l = @lens _.a[plustwo(begin) - one].b
obj = (a=(1, (a=10, b=20), 3), b=4)
@test get(obj, l) == 20
@test set(obj, l, true) == (a=(1, (a=10, b=true), 3), b=4)

ext: jl


text:
using AMQPClient

struct ExchangeDeclarationError <: Exception
    exchange_name::String
end

function declare_exchange!(chan, exchange_name, type)
    success = exchange_declare(chan, exchange_name, type)
    success || throw(ExchangeDeclarationError(exchange_name))
    exchange_name
end


text:
# Run this cell to load the graphics packages

using Plots; gr()
using Interact

# Run this cell to use a couple of definitions from the previous Tools notebook

σ(x) = 1 / (1 + exp(-x))
f(x, w) = σ(w * x)

# ------------------------------------------------------------------------------------------
# ## Multiple function parameters
#
# In notebook 6, we saw how we could adjust a parameter to make a curve fit a single data
# point. What if there is more data to fit?
#
# We'll see that as we acquire more data to fit, it can sometimes be useful to add
# complexity to our model via additional parameters.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Adding more data
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Suppose there are now two data points to fit, the previous $(x_0, y_0)$ used in notebook 6
# and also $(x_1, y_1) = (-3, 0.3)$.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
# Make an interactive plot of the function $f_w$ together with the two data points. Can you
# make the graph of $f_w$ pass through *both* data points at the *same* time?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# You should have found that it's actually *impossible* to fit both data points at the same
# time! The best we could do is to *minimise* how far away the function is from the data. To
# do so, we need to somehow balance the distance from each of the two data points.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 2
#
# Play with the slider to find the value $w^*$ of $w$ that you think has the "least error"
# in an intuitive sense.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Defining a loss function
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# How can we *quantify* some kind of overall measure of the distance from *all* of the data?
# Well, we just need to define a new loss function! One way to do so would be to sum up the
# loss functions for each data point, i.e. the sum of the squared vertical distances from
# the graph to each data point:
#
# $$L(w) = [y_0 - f_w(x_0)]^2 + [y_1 - f_w(x_1)]^2.$$
#
# Since the two pieces that we are adding have the same mathematical "shape" or structure,
# we can abbreviate this by writing it as a sum:
#
# $$L(w) = \sum_{i=0}^1 [y_i - f_w(x_i)]^2.$$
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# So now we want to find the value $w^*$ of $w$ that minimizes this new function $L$!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 3
#
# Make an interactive visualization to show the function $f_w$ and to visualize the distance
# from each of the data points.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 4
#
# After playing with this for a while, it is intuitively obvious that we cannot make the
# function pass through both data points for any value of $w$. In other words, our loss
# function, `L(w)`, is never zero.
#
# What is the minimum value of `L(w)` that you can find by altering `w`? What is the
# corresponding value of `w`?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Sums in Julia
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# To generate the above plot we used the `sum` function. `sum` can add together all the
# elements of a collection or range, or it can add together the outputs of applying a
# function to all the elements of a collection or range.
#
# Look up the docs for `sum` via
#
# ```julia
# ?sum
# ```
# if you need more information.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 5
#
# Use `sum` to add together all integers in the range 1 to 16, inclusive. What is the
# result?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 6
#
# What is the sum of the absolute values of all integers between -3 and 3? Use `sum` and the
# `abs` function.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## What does the loss function $L$ look like?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# In our last attempt to minimize `L(w)` by varying `w`, we saw that `L(w)` always seemed to
# be greater than 0.  Is that true? Let's plot it to find out!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 7
#
# Plot the new loss function $L(w)$ as a function of $w$.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Features of the loss function
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The first thing to notice is that $L$ is always positive. Since it is the sum of squares,
# and squares cannot be negative, the sum cannot be negative either!
#
# However, we also see that although $L$ dips close to $0$ for a single, special value $w^*
# \simeq 0.4$, it never actually *reaches* 0! Again we could zoom in on that region of the
# graph to estimate it more precisely.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We might start suspecting that there should be a better way of using the computer to
# minimize $L$ to find the location $w^*$ of the minimum, rather than still doing everything
# by eye. Indeed there is, as we will see in the next two notebooks!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Adding more parameters to the model
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# If we add more parameters to a function, we may be able to improve how it fits to data.
# For example, we could define a new function $g$ with another parameter, a shift or
# **bias**:
#
# $$g(x; w, b) := \sigma(w \, x) + b.$$
# ------------------------------------------------------------------------------------------

g(x, w, b) = σ(w*x) + b

# ------------------------------------------------------------------------------------------
# *Note*: In the last notebook, we added parameters to a sigmoid function to get the form
# $$\sigma(w \, x + b)$$ and here we are working with the form $$\sigma(w \, x) + b$$
# instead. Both of these are valid ways to apply a bias to a function! In machine learning
# terminology, the first corresponds to a neural network with a single neuron, while the
# second can be thought of as a network with two layers.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 8
#
# Make an interactive visualization with two sliders for $w$ and $b$. Play with the sliders
# to try and fit *both* data points at once!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 9
#
# For what values of `w` and `b` does the line pass through both points?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Fitting both data points: a loss function
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Following the procedure that we used when we had a single parameter, we can think of the
# fitting procedure once again as minimizing a suitable *loss function*. The expression for
# the loss function is almost the same, except that now $f$ has two parameters $w$ and $b$,
# so the loss function $L_2$ is itself a function of $w$ *and* $b$:
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# $$L_2(w, b) = \sum_i [y_i - f(x_i; w, b)]^2.$$
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# So we want to minimize this loss function over *both* of the parameters $w$ and $b$! Let's
# plot it.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 10
#
# Define the function `L2` in Julia.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## More data
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# If we add more data, however, we will again not be able to fit all of the data; we will
# only be able to attain a "best fit".
#
# Let's create `xs` and `ys` with some more data:
# ------------------------------------------------------------------------------------------

xs = [2, -3, -1, 1]
ys = [0.8, 0.3, 0.4, 0.4]

# ------------------------------------------------------------------------------------------
# #### Exercise 11
#
# a) Make an interactive visualization of the function $f$ and the data. Try to find the
# values of $w$ and $b$ that give the best fit.
#
# b) Define the loss function and plot it.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 12
#
# We've seen the loss function, $$L_{2D}(w, b) = \sum_i[\mathrm{ys}_i - g(\mathrm{xs}_i, w,
# b)]^2,$$ written in Julia as
#
# ```julia
# L2(w, b) = sum( (ys .- g.(xs, w, b)) .^ 2 )
# ```
#
# a few times now. To ensure you understand what this function is doing, implement your own
# loss function using the commented code below. Do this without using `sum`, `abs2`, or
# broadcasting dot syntax (for example, `.-`). Hint: you'll want to use a `for` loop to do
# this.
# ------------------------------------------------------------------------------------------

# function myL2(w, b)
#     loss = 0.0

#     return loss
# end

# ------------------------------------------------------------------------------------------
# Now that you've defined `L2`, we can plot it using `surface`!
# ------------------------------------------------------------------------------------------

plotlyjs()

ws = -2:0.05:2
bs = -2:0.05:2

surface(ws, bs, L2, alpha=0.8, zlims=(0,3))

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/08. Tools - Multiple function parameters.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
function update_gam!(i::Int, n::Int, j::Int, s::State, c::Constants, d::Data)
  k = s.lam[i][n]
  if k > 0
    z = s.Z[j, k]
    logpriorVec = log.(s.eta[z][i, j, :])
    loglikeVec = logpdf.(Normal.(mus(z, s, c, d), sqrt(s.sig2[i])), s.y_imputed[i][n, j])
    logPostVec = logpriorVec .+ loglikeVec
    s.gam[i][n, j] = MCMC.wsample_logprob(logPostVec)
  else
    s.gam[i][n, j] = 0
  end
end

function update_gam!(s::State, c::Constants, d::Data)
  for i in 1:d.I
    for j in 1:d.J
      for n in 1:d.N[i]
        update_gam!(i, n, j, s, c, d)
      end
    end
  end
end

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch


text:
push!(LOAD_PATH,"../src/")
using Documenter, RDKitMinimalLib

makedocs(
    sitename="RDKitMinimalLib.jl",
    pages=[
        "Home" => "index.md",
        "RDKitMinimalLib" => [
            "I/O" => "io.md",
            "Calculators" => "calculators.md",
            "Standardization" => "standardization.md",
            "Coordinates" => "coordinates.md",
            "Drawing" => "drawing.md",
            "Modification" => "modification.md",
            "Substructure" => "substructure.md",
            "Utils" => "utils.md",
        ],
    ],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    )
)

deploydocs(
    #options
    repo="github.com/eloyfelix/RDKitMinimalLib.jl.git",
    devbranch = "main"
)


text:
module BenchSeq

import Random
using BenchmarkTools
using ParallelRadixSort: msdsort_seq

Random.seed!(1234)

n = 10_000
datasets = [
    # (label, setup)
    ("UInt8 (wide)", rand(UInt8, 100_000)),
    ("UInt64 (wide)", rand(UInt64, 100_000)),
    ("Float64 (wide)", randn(100_000)),
    # `Base` has `Base.Sort.sort_int_range!` for those cases:
    # ("UInt8 (narrow)", rand(UInt8(0):UInt8(9), 100_000)),
    # ("UInt64 (narrow)", rand(UInt64(0):UInt64(9), 100_000)),
    ("Float64 (narrow)", rand(0:0.1:0.9, 100_000)),
    ("String", [join(rand('a':'z', 10)) for _ in 1:10_000]),
    ("Vector{Int}", [rand(0:9, 10) for _ in 1:10_000]),
]
# Poor performance:
# * narrow
# * Vector{Int}

suite = BenchmarkGroup()

for (label, xs) in datasets
    s = suite[label] = BenchmarkGroup()
    s["radix"] = @benchmarkable(msdsort_seq($xs))
    s["base"] = @benchmarkable(sort($xs))
end

end  # module
BenchSeq.suite


text:
module ImPlot

using CImGui
import CImGui: ImVec2, ImVec4

for i in instances(CImGui.ImGuiCond_)
    @eval import CImGui.LibCImGui: $(Symbol(i))
end

# include("ImPlot_jl")

# using .LibCImPlot
using CImGui.LibCImGui
export IMPLOT_AUTO, IMPLOT_AUTO_COL
export ImPlotPoint, ImPlotRange, ImPlotLimits

const ImPlotData = Union{Float32, Float64, Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64}
const IMPLOT_ENUMS = [ImPlotFlags_, ImPlotAxisFlags_, ImPlotCol_, ImPlotStyleVar_, ImPlotMarker_,
          ImPlotColormap_, ImPlotLocation_, ImPlotOrientation_, ImPlotYAxis_]

LibCImPlot = LibCImGui

# SKIPLIST = [
#     :ImPlot_BeginPlot
#     :ImPlot_PlotText
#     :ImPlot_PlotLineG
#     :ImPlot_PlotScatterG
#     :ImPlot_PlotShadedG
#     :ImPlot_PlotBarsG
#     :ImPlot_PlotBarsHG
#     :ImPlot_PlotDigitalG
#     :ImPlot_PlotBars
#     :ImPlot_PlotBarsH
#     :ImPlot_PlotDigital
#     :ImPlot_PlotErrorBars
#     :ImPlot_PlotHeatmap
#     :ImPlot_PlotLine
#     :ImPlot_PlotScatter
#     :ImPlot_PlotStairs
#     :ImPlot_PlotShaded
#     :ImPlot_PlotBars
#     :ImPlot_PlotErrorBars
#     :ImPlot_PlotErrorBarsH
#     :ImPlot_PlotStems
#     :ImPlot_PlotPieChart
#     :ImPlot_PlotDigital
#     :ImPlot_PlotImage
#     :ImPlot_PlotLine
#     :ImPlot_PlotShaded

#     :ImPlot_GetPlotPos
#     :ImPlot_GetPlotSize
#     :ImPlot_SetLegendLocation
#     :ImPlot_GetNextPlotLimits
#     :ImPlot_GetPlotLimits
#     :ImPlot_SetNextPlotLimits
#     :ImPlot_SetPlotLimits
#     :ImPlot_GetPlotMousePos
# ]
# for name in filter(startswith("ImPlot_"), string.(propertynames(LibCImGui)))
#     lhs = Symbol(replace(name, "ImPlot_"=>""))
#     rhs = Symbol(name)
#     if rhs in SKIPLIST
#         continue
#     end
#     @eval $lhs = $rhs
# end

# Export plot flags
for i in IMPLOT_ENUMS
    for j in instances(i)
        @eval export $(Symbol(j))
    end
end

function BeginPlot(title_id::String, x_label = C_NULL, y_label = C_NULL, size::ImVec2 = ImVec2(-1,0);
                    flags = ImPlotFlags_None,
                    x_flags = ImPlotAxisFlags_None,
                    y_flags = ImPlotAxisFlags_None,
                    y2_flags = ImPlotAxisFlags_None,
                    y3_flags = ImPlotAxisFlags_None,
                    y2_label = C_NULL,
                    y3_label = C_NULL,
                    )::Bool

    ImPlot_BeginPlot(title_id, x_label, y_label, size,
                         flags, x_flags, y_flags, y2_flags, y3_flags, y2_label, y3_label)
end

include("constructors.jl")
include("lines.jl")
include("stairs.jl")
include("shaded.jl")
include("scatter.jl")
include("heatmap.jl")
include("digital.jl")
include("barchart.jl")
include("piechart.jl")
include("errorbars.jl")
include("stems.jl")
include("other.jl")
include("util.jl")
include("color.jl")
include("styling.jl")

# SKIPLIST = [
#     :ImPlot_BeginPlot
#     :ImPlot_PlotText
#     :ImPlot_PlotLineG
#     :ImPlot_PlotScatterG
#     :ImPlot_PlotShadedG
#     :ImPlot_PlotBarsG
#     :ImPlot_PlotBarsHG
#     :ImPlot_PlotDigitalG
#     :ImPlot_PlotBars
#     :ImPlot_PlotBarsH
#     :ImPlot_PlotDigital
#     :ImPlot_PlotErrorBars
#     :ImPlot_PlotHeatmap
#     :ImPlot_PlotLine
#     :ImPlot_PlotScatter
#     :ImPlot_PlotStairs
#     :ImPlot_PlotShaded
#     :ImPlot_PlotBars
#     :ImPlot_PlotErrorBars
#     :ImPlot_PlotErrorBarsH
#     :ImPlot_PlotStems
#     :ImPlot_PlotPieChart
#     :ImPlot_PlotDigital
#     :ImPlot_PlotImage
#     :ImPlot_PlotLine
#     :ImPlot_PlotShaded

#     :ImPlot_GetPlotPos
#     :ImPlot_GetPlotSize
#     :ImPlot_SetLegendLocation
#     :ImPlot_GetNextPlotLimits
#     :ImPlot_GetPlotLimits
#     :ImPlot_SetNextPlotLimits
#     :ImPlot_SetPlotLimits
#     :ImPlot_GetPlotMousePos
# ]
for name in filter(startswith("ImPlot_"), string.(propertynames(LibCImGui)))
    lhs = Symbol(replace(name, "ImPlot_"=>""))
    rhs = Symbol(name)
    # if rhs in SKIPLIST
    #     continue
    # end
    try
        @eval $lhs = $rhs
    catch err
        continue
    end
end

end # module


text:
using MRBC
using Test

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


text:
module TestDoctest

using Documenter
using Test
using Transducers

# Workaround: Failed to evaluate `CurrentModule = Transducers` in `@meta` block.
@eval Main import Transducers

@testset "/docs" begin
    doctest(Transducers; manual=true)
end

@testset "/test/doctests" begin
    doctest(joinpath((@__DIR__), "doctests"), Module[])
end

end  # module


text:
# ---
# title: 1011. Capacity To Ship Packages Within D Days
# id: problem1011
# author: Indigo
# date: 2021-01-30
# difficulty: Medium
# categories: Array, Binary Search
# link: <https://leetcode.com/problems/capacity-to-ship-packages-within-d-days/description/>
# hidden: true
# ---
#
# A conveyor belt has packages that must be shipped from one port to another
# within `D` days.
#
# The `i`-th package on the conveyor belt has a weight of `weights[i]`.  Each
# day, we load the ship with packages on the conveyor belt (in the order given
# by `weights`). We may not load more weight than the maximum weight capacity of
# the ship.
#
# Return the least weight capacity of the ship that will result in all the
# packages on the conveyor belt being shipped within `D` days.
#
#
#
# **Example 1:**
#
#
#
#     Input: weights = [1,2,3,4,5,6,7,8,9,10], D = 5
#     Output: 15
#     Explanation:
#     A ship capacity of 15 is the minimum to ship all the packages in 5 days like this:
#     1st day: 1, 2, 3, 4, 5
#     2nd day: 6, 7
#     3rd day: 8
#     4th day: 9
#     5th day: 10
#
#     Note that the cargo must be shipped in the order given, so using a ship of capacity 14 and splitting the packages into parts like (2, 3, 4, 5), (1, 6, 7), (8), (9), (10) is not allowed.
#
#
# **Example 2:**
#
#
#
#     Input: weights = [3,2,2,4,1,4], D = 3
#     Output: 6
#     Explanation:
#     A ship capacity of 6 is the minimum to ship all the packages in 3 days like this:
#     1st day: 3, 2
#     2nd day: 2, 4
#     3rd day: 1, 4
#
#
# **Example 3:**
#
#
#
#     Input: weights = [1,2,3,1,1], D = 4
#     Output: 3
#     Explanation:
#     1st day: 1
#     2nd day: 2
#     3rd day: 3
#     4th day: 1, 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= D <= weights.length <= 50000`
#   * `1 <= weights[i] <= 500`
#
#
## @lc code=start
using LeetCode

function ship_within_days(weights::Vector{Int}, D::Int)
    lo, hi = maximum(weights), sum(weights)
    while lo < hi
        mid = (lo + hi) ÷ 2
        cnt = 1
        tt = 0
        for w in weights
            tt += w
            if tt > mid
                cnt += 1
                tt = w
            end
        end
        cnt > D ? (lo = mid + 1) : (hi = mid)
    end
    return lo
end
## @lc code=end


text:
try
@testset "Lint Helper" begin
    # we must run these scripts explicitly since they are not real packages in standard directories
    include("DEMOMODULE.jl") # this provide the macro that generates functions
    include("DEMOMODULE2.jl") # this uses the first module's macro. It would export the generated functions

    msgs = lintfile("DEMOMODULE2.jl")
    @test_broken isempty(msgs)

    msgs = lintfile("DEMOMODULE3.jl")
    @test_broken msgs[1].code == :E311
    @test_broken contains(msgs[1].message, "cannot find include file")
end
end


text:
module DoubleBLAS
# support for (fairly) efficient Linear Algebra with DoubleFloats

#FIXME: change to explicit lists because namespace pollution is epidemic
using DoubleFloats
using LinearAlgebra

using SIMD
using UnsafeArrays

using Base.Threads

# steal some internals
using LinearAlgebra: lapack_size, BlasInt, checknonsingular, MulAddMul
using LinearAlgebra: _modify!

if VERSION < v"1.2"
    using LinearAlgebra: has_offset_axes
    require_one_based_indexing(x...) =
        has_offset_axes(x...) && throw(ArgumentError("not implemented "
                                                  * "for offset axes"))
else
    using Base: require_one_based_indexing
end

# stuff to extend
import LinearAlgebra: rmul!, lmul!, ldiv!
# see subordinate files for some others

export refinedldiv

# most public functions are defined in LinearAlgebra

################################################################
# SIMD config

# width of SIMD structures to use
const Npref = 8
# Using Npref = 8 on a system w/ 256b SIMD registers imposes a 10% penalty
# for most modest-size Double64 problems.
# There is no obvious penalty for large problems, and this
# presumably gives better performance on 512b systems.
# Kudos to the LLVM people.

################################################################
# Multi-threading internals

# TODO: implement set_num_threads

const mt_thresholds = Dict{Symbol,Any}()

"""
set_mt_threshold(n::Real, problem::Symbol)

Set the size threshold for multi-threading in the DoubleBLAS package to `n`,
for matrix operations of class `problem`.
"""
function set_mt_threshold(n::Real, problem::Symbol)
    if problem ∈ keys(mt_thresholds)
        destref = mt_thresholds[problem]
        destref[] = Float64(n)
    else
        throw(ArgumentError("unrecognized problem key $problem: valid keys are $(keys(mt_thresholds))"))
    end
    nothing
end

function get_mt_threshold(problem::Symbol)
    if problem ∈ keys(mt_thresholds)
        srcref = mt_thresholds[problem]
    else
        throw(ArgumentError("unrecognized problem key $problem: valid keys are $(keys(mt_thresholds))"))
    end
    srcref = mt_thresholds[problem]
    srcref[]
end

################################################################
# SIMD internals

@generated function vgethi(xv::StridedVecOrMat{DoubleFloat{T}},i0,::Type{Vec{N,T}}) where {N,T}
    quote
        $(Expr(:meta, :inline))
        $(Expr(:meta, :inbounds))
        Vec{N,T}(tuple($([:(HI(xv[i0+$i])) for i in 1:N]...)))
    end
end

@generated function vgetlo(xv::StridedVecOrMat{DoubleFloat{T}},i0,::Type{Vec{N,T}}) where {N,T}
    quote
        $(Expr(:meta, :inline))
        $(Expr(:meta, :inbounds))
        Vec{N,T}(tuple($([:(LO(xv[i0+$i])) for i in 1:N]...)))
    end
end

@inline function vputhilo!(xv::StridedVecOrMat{DoubleFloat{T}},i0,
                              zhi::Vec{N,T}, zlo::Vec{N,T}) where {N,T}
    @inbounds for i=1:N
        xv[i0+i] = DoubleFloat{T}((zhi[i],zlo[i]))
    end
end

@generated function vgethire(xv::StridedVecOrMat{Complex{DoubleFloat{T}}},i0,::Type{Vec{N,T}}) where {N,T}
    quote
        $(Expr(:meta, :inline))
        $(Expr(:meta, :inbounds))
        Vec{N,T}(tuple($([:(HI(xv[i0+$i].re)) for i in 1:N]...)))
    end
end

@generated function vgetlore(xv::StridedVecOrMat{Complex{DoubleFloat{T}}},i0,::Type{Vec{N,T}}) where {N,T}
    quote
        $(Expr(:meta, :inline))
        $(Expr(:meta, :inbounds))
        Vec{N,T}(tuple($([:(LO(xv[i0+$i].re)) for i in 1:N]...)))
    end
end

@generated function vgethiim(xv::StridedVecOrMat{Complex{DoubleFloat{T}}},i0,::Type{Vec{N,T}}) where {N,T}
    quote
        $(Expr(:meta, :inline))
        $(Expr(:meta, :inbounds))
        Vec{N,T}(tuple($([:(HI(xv[i0+$i].im)) for i in 1:N]...)))
    end
end

@generated function vgetloim(xv::StridedVecOrMat{Complex{DoubleFloat{T}}},i0,::Type{Vec{N,T}}) where {N,T}
    quote
        $(Expr(:meta, :inline))
        $(Expr(:meta, :inbounds))
        Vec{N,T}(tuple($([:(LO(xv[i0+$i].im)) for i in 1:N]...)))
    end
end

@inline function vputhilo!(xv::StridedVecOrMat{Complex{DoubleFloat{T}}},i0,
                           zrhi::Vec{N,T}, zrlo::Vec{N,T},
                           zihi::Vec{N,T}, zilo::Vec{N,T}) where {N,T}
    @inbounds for i=1:N
        xv[i0+i] = complex(DoubleFloat{T}((zrhi[i],zrlo[i])),
                           DoubleFloat{T}((zihi[i],zilo[i])))
    end
end

include("ops.jl")
include("dots.jl")
include("axpy.jl")
include("level2.jl")
include("givens.jl")

include("gemm.jl")
include("triangular.jl")

include("lu.jl")
include("chol.jl")

include("refine.jl")

end # module

max_stars_repo_name: RalphAS/DoubleBLAS

max_issues_repo_name: RalphAS/DoubleBLAS

max_forks_repo_name: RalphAS/DoubleBLAS


text:
# using
using RCall
using GeoArrays
using Dates
using NCDatasets

R"""
getwd()
# load("dTmax.rda")
"""

file = "data-raw/CN05.1_Rhu_1961_2017_daily_025x025.nc"
gdalinfo(file)

Date(2014,1,29):Dates.Month(1):Date(2014,07,29)
ds = NCDataset(file)
ds.attrib

dates = ds["time"][:]
years = year.(dates)

year_grps = unique(years)
arr = ds

for year in year_grps
    println(year)
end

max_stars_repo_name: CUG-hydro/heatwave

max_issues_repo_name: CUG-hydro/heatwave

max_forks_repo_name: CUG-hydro/heatwave


text:
__precompile__()

module Pluck

  export pluck, pluck!

  function _pluck(cur_array::AbstractArray, remove_item=false)
    cur_length = length(cur_array)
    ( iszero(cur_length) ) && throw(BoundsError())

    cur_index = rand(1:cur_length)
    cur_value = cur_array[cur_index]

    remove_item && deleteat!(cur_array, cur_index)

    cur_value
  end

  function pluck(cur_array::AbstractArray)
    _pluck(cur_array, false)
  end

  function pluck!(cur_array::AbstractArray)
    _pluck(cur_array, true)
  end

  function pluck(cur_iterator::Union{Base.KeyIterator, Base.ValueIterator})
    pluck(collect(cur_iterator))
  end

end


text:
"""
    Actuarial

defines basic types, methods in actuarial science (annuity, life & non-life insurances), including:
1. Types
    1. ::CashFlow, generic type of finite cash flows, defined by pairs of (time, payment)
    2. ::InfCashFlow, generic type of infinite cash flows, equivalent to an infinite annuity
    1. ::StdAnnuity, standard deterministic annuity of discrete time
2. Methods:
    1. AnnuityValue() & PV() & FV(), compute discounted values of a given annuity

There are some non-export members:
1. ::LifeTable, type for life tables
    1. Note: I provide data for China Life Insurance Mortality Table (2000-2003), CL(2000-2003); they are experience life table for life insurance pricing
    1. LT_CL1_03_NonPension::LifeTable, CL(2000-2003) for male,   for non-pension affair
    1. LT_CL2_03_NonPension::LifeTable, CL(2000-2003) for female, for non-pension affair
    1. LT_CL3_03_NonPension::LifeTable, CL(2000-2003) for male,   for pension affair
    1. LT_CL4_03_NonPension::LifeTable, CL(2000-2003) for female, for pension affair

"""
module Actuarial
    import Base  # explicitly import Base to reload basic operators (such as +, -)
    export cashflow, iscashflow, PV, FV, items, cfcompress  # export generic methods
    export StdAnnuity, Payment, CashFlow  # export types

# ==================================
"""
    Payment( Time::Real, Payment::Real )

a synonym of Pair{Float64, Float64}, a "Time => Payment" relationship;
the constructor can automatically convert input args to Float64;
allows finite/infinite Time/Payment; please refer to the doc of CashFlow for more information;
"""
Payment = Pair{Float64, Float64}
# --------------- Atomic Operations
# NOTE: overloads basic Arithmetical operators for Payment
# NOTE: two payments cannot be added, so, no def for + & -
Base.:-( pay::Payment ) = Payment( pay.first, - pay.second )  # but it can be negative!
Base.:*( pay::Payment, k::Real ) = Payment( pay.first, Float64(k * pay.second) )
Base.:*( k::Real, pay::Payment ) = Payment( pay.first, Float64(k * pay.second) )  # exchange-able principle does not always stand, needs to define it manually
Base.:/( pay::Payment, k::Real ) = Payment( pay.first, Float64(pay.second / k) )
Base.:/( k::Real, pay::Payment ) = Payment( pay.first, Float64(k / pay.second) )
# --------------- Type Assertion
Base.isfinite( pay::Payment ) = ( isinf(pay.first) | isinf(pay.second) )  ?  (return false::Bool) : (return true::Bool)
Base.isinf( pay::Payment ) = ( isinf(pay.first) | isinf(pay.second) )  ?  (return true::Bool) : (return false::Bool)
# --------------- Operations
"""
    PV( CF::CashFlow, IntRate::Float64 ; DiscountTo::Union{Int,Float64} = 0.0, IsDiscrete::Bool = false )

computes given payment's present value at a finite moment DiscountTo::Float64;

Receives:
1. CF::Payment, a given payment
2. IntRate::Float64, a constant interest rate to discount at, >= 0
3. DiscountTo::Real, the moment to discount to, 0.0 by default
4. IsDiscrete::Bool, how to deal with the non-integer (digit) part of period, true by default
Returns:
1. Res::Float64, PV of the given payment

Mathematics:
1. if IsDiscrete is false, uses exponential discounting: `` {PV} = \\exp(-rt) ``, where r is IntRate
2. if IsDiscrete is true, uses linear interpolation to discount inner-period (non-integer, digit) part: `` {PV}_{t=1} = {PV}_{t=1.25} / ( 1 + r * 0.25 ) ``, where r is IntRate
3. for infinite time/payment, we have principles (we use (time,payment) to denote a Payment):
    1. we use R to denote finite real number
    1. non-zero R = finite DiscountTo + one of: (R,R)
    2. zero = finite DiscountTo + one of: (Inf,R)
    3. Inf = finite DiscountTo + one of: (-Inf,Inf), (R,Inf)
    5. -Inf = finite DiscountTo + one of: (R,-Inf),(-Inf,-Inf)
    6. NaN = finite DiscountTo + one of: (Inf,Inf),(Inf,-Inf),(-Inf,R)
"""
PV( CF::Payment, IntRate::Float64 ; DiscountTo::Real = 0.0, IsDiscrete::Bool = true ) = begin
    # valid & declare
    isfinite(DiscountTo)  ||  error("expects a finite DiscountTo but received: $(DiscountTo)")
    IntRate >= 0.0  ||  error("expects a non-negative interest rate but received: $(IntRate)")
    # compute
    if isfinite(CF.first) & isfinite(CF.second)
        if IsDiscrete
            # first, discount proportional part to the floor-integer
            TmpVal = CF.second / ( 1.0 + IntRate * ( CF.first - floor(CF.first) ) )
            # then, discount it to floor(DiscountTo) moment
            TmpVal /= ( 1.0 + IntRate ) ^ ( floor(CF.first) - floor(DiscountTo) )
            # finally, discount to proportional DiscountTo moment
            Res = TmpVal * ( 1.0 + IntRate * ( DiscountTo - floor(DiscountTo) ) )
            return Res::Float64
        else
            Res = CF.second * exp( IntRate * ( DiscountTo - CF.first ) )
            return Res::Float64
        end
    elseif CF.first == Inf & isfinite(CF.second)
        return 0.0::Float64
    elseif ( CF.first == -Inf & CF.second == Inf ) | ( isfinite(CF.first) & CF.second == Inf )
        return Inf::Float64
    elseif ( isfinite(CF.first) & CF.second == -Inf ) | ( CF.first == -Inf & CF.second == -Inf )
        return -Inf::Float64
    else
        return NaN::Float64
    end
end # end functio PV
# ==================================
"""
    CashFlow

a synonym of Vector{Pair{Float64,Float64}},
whose element is (Time, Signed Payment);
denotes a series of money flows at different moments;
allows duplicated time/moments unless concise() it;

Mathematics:
1. The domain of Time is: {-Inf,R,Inf}, where R is finite real number
2. The domain of Money is: {-Inf,R,Inf}
3. then we define (mathematically) some operations on a single CashFlow
3. the principles of PV, please refer to the documentation of type Payment
"""
CashFlow = Vector{Payment}
# ----------------- Type Assertion
"""
    iscashflow( Obj::Any )

checks if an instance is one of CashFlow
"""
iscashflow( Obj::Any ) = isa( Obj, CashFlow )
# ----------------- Constructors
cashflow( Time::Real, Cash::Real ) = [Payment( Time, Cash )]
cashflow( Time::Real, Cash::Union{AbstractRange, Vector{T} where T <: Real} ) = begin
    local Res = CashFlow()  # initilize an instance to return
    for x in Cash; push!( Res, Payment( Time, x ) ); end;  # push new (time, payment) pairs to Res::CashFlow
    return Res::CashFlow
end
cashflow( Time::Union{AbstractRange, Vector{T} where T <: Real}, Cash::Real ) = begin
    local Res = CashFlow()
    for x in Time; push!( Res, Payment( x, Cash ) ); end;
    return Res::CashFlow
end
"""
    cashflow( Time::Vector{T} where T <: Real, Payment::Vector{T} where T <: Real )

creates a CashFlow from two vectors: one for time/moments, one for payments;
returns a CashFlow instance;
duplicated time/moments allowed;
"""
cashflow( Time::Union{AbstractRange, Vector{T} where T <: Real}, Cash::Union{AbstractRange, Vector{T} where T <: Real} ) = begin
    # valid & declare
    local T = length(Time)  # size of the cash flow
    T == length(Cash)  ||  error("expects two equal-length vectors but received: ( $(T), $(length(Cash)) )")
    local Res = CashFlow()  # initilize an instance to return
    for t in 1:T
        push!( Res, Payment( Time[t], Cash[t] ) )
    end
    return Res::CashFlow
end  # end function cashflow()
"""
    cashflow( TimeCash::Matrix{T} where T <: Real )

creates a CashFlow from a T*2 matrix (synonym of Array{T,2}),
where the first column is time/moments of each payment,
and the second column is payment amounts;
returns a CashFlow;
"""
cashflow( TimeCash::Matrix{T} where T <: Real ) = begin
    # valid & declare
    local T = size(TimeCash)  # size of the given matrix
    T[2] == 2  ||  error("expects a T*2 matrix but received: $(T)")
    local Res = CashFlow()
    for t in 1:T[1]; push!( Res, Payment( TimeCash[t,1], TimeCash[t,2] ) ); end;
    return Res::CashFlow
end  # end function cashflow()
"""
    cashflow( Time2Cash::Dict{T,T} where T <: Real )

creates a CashFlow from a number-to-number Dict (Hash table),
where the keys are time/moments of each payment,
and the values are payment amounts;
returns a CashFlow;
"""
cashflow( Time2Cash::Dict{T,T} where T <: Real ) = begin
    # valid & declare
    local Res = CashFlow()
    for x in Time2Cash; push!( Res, Payment( x.first, x.second ) ); end;
    return Res::CashFlow
end  # end function cashflow()
# ----------------- Operations
"""
    PV( CF::CashFlow, IntRate::Float64 ; DiscountTo::Union{Int,Float64} = 0.0, IsDiscrete::Bool = true )

computes the present value of a CashFlow at time CountTo::Union{Int,Float64} at a constant non-negative interest rate IntRate::Float64;

Receives:
1. CF::CashFlow, cash flow to discount
2. IntRate::Float64, a constant interest rate to discount, >=0, in digit rather than percentage
3. DiscountTo::Union{Int,Float64}, the moment to discount to
4. IsDiscrete::Bool, to discount as a discrete cash flow or a continuous cash flow, true by default
Returns:
1. Res::Float64, PV of the given cash flow
Depends on:
1. PV()

Note:
1. when IsDiscrete is false, money is discounted as `` PV = \\exp( - r (t-t_0) ) ``
1. when IsDiscrete is true, inner-year part is discounted at linear interpolation, e.g.
`` \\text{PV}_{t=1} = \\text{FV}_{t=1.23} \\cdot \\frac{1}{1+ 10% \\cdot 0.23 } ``
"""
function PV( CF::CashFlow, IntRate::Float64 ; DiscountTo::Real = 0.0, IsDiscrete::Bool = true )
    # NOTE: validations are did inside PV(::Payment,...)
    local Res::Float64 = 0.0  # result
    # compute
    for tmppayment in CF
        Res += PV( tmppayment, IntRate, DiscountTo = DiscountTo, IsDiscrete = IsDiscrete )
    end
    return Res::Float64
end # end function PV
# ----------------- Arithmetical Operations
# NOTE: + & - on CashFlow are set-level operations, not on element (Payment)
# NOTE: but - (get negative sign) works on single Payment
# NOTE: * & / are defined at element-level, not on set-level (CashFlow), please use .* and ./
"""
    Base.:+( cfa::CashFlow, cfb::CashFlow )

overloads addition (+) for CashFlow;
it is a kind of collection-level operation (not defined on element!);
two CashFlow are connected, where payments with different payment amounts but at the same moments are all kept;
returns a new CashFlow instance;
"""
Base.:+( cfa::CashFlow, cfb::CashFlow ) = cat( cfa, cfb , dims = 1)
"""
    Base.:-( cfa::CashFlow, cfb::CashFlow )

overloads minus (-) for CashFlows;
it is a kind of set operation (not defined on element!);
payments in cfb has opposite signs then added to cfa;
returns a new CashFlow instance;
"""
Base.:-( cfa::CashFlow, cfb::CashFlow ) = begin
    return ( cfa + ( .- cfb ) )::CashFlow
end  # end function -
# ----------------- Collection Operations
"""
    Base.keys( cf::CashFlow )

overloads keys() which extracts all time (Payment.first) tags as a vector;
orders are kept;
returns a Vector{Float64};
"""
function Base.keys( cf::CashFlow )
    local Res::Vector{Float64} = Vector{Float64}()
    for x in cf; push!( Res, x.first ); end
    return Res::Vector{Float64}
end  # end function Base.keys
"""
    Base.values( cf::CashFlow )

overloads values() which extracts all payment amount (Payment.second) tags as a vector;
orders are kept;
returns a Vector{Float64};
"""
function Base.values( cf::CashFlow )
    local Res::Vector{Float64} = Vector{Float64}()
    for x in cf; push!( Res, x.second ); end
    return Res::Vector{Float64}
end  # end function Base.values
"""
    items( cf::CashFlow ; sorted = false )

defines a convenient Py-style method to extract the keys and values of a given cash flow at the same time;
returns a Tuple{Vector{Float64},Vector{Float64}},
where the first element is keys vector, the second is values vector;
"""
items( cf::CashFlow ; sorted = false ) = begin
    local Res::Tuple{Vector{Float64},Vector{Float64}} = sorted ? ( keys(sort(cf)), values(sort(cf)) )  :  ( keys(cf), values(cf) )
    return Res::Tuple{Vector{Float64},Vector{Float64}}
end # end function Base.items
"""
    Base.cfcompress( cf::CashFlow )

compresses a cash flow, summing all duplicated records (same time but different payment amount),
and dropping zero payments (payment == 0);
returns a CashFlow with unique members;

Depends on:
1. Base.keys(): overloaded
"""
function cfcompress( cf::CashFlow )
    # declare
    local UniqKeys::Vector{Float64} = unique(keys(cf))  # get unique keys for new cash flow
    # new a new CashFlow with all-zero payment amounts but unique keys
    # NOTE: it can always be safely converted to a Dict because keys are unique now
    # NOTE: when converted to a Dict, we can use Res[keys] to safely index records
    # NOTE: according to our definition, a Dict{Float64,Float64} can always be converted a CashFlow
    local Res = Dict{Float64,Float64}( cashflow( UniqKeys, 0.0 ) )
    # compute, NOTE: because we have 0.0 as bases, now it is easy to just add (possible duplicated) values to their new locations
    for x in cf
        Res[x.first] += x.second
    end
    # finally, converted Res to a CashFlow and return
    return cashflow( Res )::CashFlow
end # end function unique
# -----------------------------

# ==================================
"""
    AbstractAnnuity

a super abstract type for different annuities
"""
abstract type AbstractAnnuity <: Any end
# --------------------
"""
    StdAnnuity

type declare of standard deterministic annuity of discrete time (pays 1 per time, start at time 0 the begin of period 1); it has fields:
1. IsFinite::Bool, if this annuity is a finite annuity, true for finite; true by default
1. IsDue::Bool  # if it is a Annuity Due (pays at the beginning of period); true for Annuity Due, false for Ordinary Annuity (pays at the end of period)
2. T::Int, the duration of this annuity, no less than 1 if IsFinite is true; when IsFinite is false, it will be neglected
3. IntRate::Float64, interest rate of this annuity, positive; discounting rate can be converted from it `` d = \\frac{i}{1+i} ``; ; for a contineous annuity, the force of interest is computed through `` \\delta = \\ln{1+i} ``

Note:
1. deferred & expired annuity can be easily computed through functions, no need to define in a type;
2. you can always construct a generic multi-inyear-payment annuity with standard annuities like constructing a portfolio!
"""
struct StdAnnuity <: AbstractAnnuity
    IsFinite::Bool  # type of annuity, true for finite, false for infinite
    IsDue::Bool  # if it is a Annuity Due (pays at the beginning of period); true for Annuity Due, false for Ordinary Annuity (pays at the end of period)
    T::Int  # period number, only works when IsFinite == true
    IntRate::Float64  # compound interest rate, can be converted to corresponding discounting rate `` d = \\frac{i}{1+i} ``; for a contineous annuity, the force of interest is computed through `` \\delta = \\ln{1+i} ``
    # ----
    """
        StdAnnuity( T::Int, IntRate::Float64 ; IsFinite::Bool = true, IsDue::Bool = false )

    Constructor, where initilize a finite ordinary annuity by default;
    when IsFinite is false, T::Int won't work;
    expects a T >= 1 (when IsFinite is true), and IntRate > 0;
    expects an interest rate in digits but not in percentage, e.g. 3% -> 0.03;
    """
    function StdAnnuity( T::Int, IntRate::Float64 ; IsFinite::Bool = true, IsDue::Bool = false )
        # validation
        IsFinite  &&  ( T >= 1 || error("Annuity type expects a T >= 1 but receives $(T)") )
        IntRate > 0.0  ||  error("Annuity type expects a positive interest rate but receives $(IntRate)")
        # new
        new( IsFinite, IsDue, T, IntRate )
    end # end constructor StdAnnuity
end # end struct StdAnnuity
# ---------------
"""
    AnnuityValue( a::StdAnnuity ; AnnuityStartTime::Int = 0, DiscountTo::Int = 0 )

computes the present value of a defined annuity at a specific time/moment;

Receives:
1. a::StdAnnuity, an StdAnnuity instance
2. AnnuityStartTime::Int, [in Z], an integer, negative or 0 or positive, the time/moment when Annuity starts, equals 0 by default; if > 0, the annuity becomes a Deferred Annuity; if < 0, the function computes the value of this annuity at that moment
3. DiscountTo::Int, [in Z], an interger, indicating what moment to discount the given annuity to; 0 by default
Returns:
1. Res::Float64, the discounted present value of given annuity

Note:
1. this is a generic method to compute the discounted value of an annuity at any time/moment; you may easily compute a Deferred Annuity or an Expired Annuity
"""
function AnnuityValue( a::StdAnnuity ; AnnuityStartTime::Int = 0, DiscountTo::Int = 0 )
    if a.IsFinite  # case: finite annuity
        if a.IsDue # case: finite StdAnnuity Due
            return ( ( 1.0 - 1.0 / (1.0 + a.IntRate) ^ a.T ) / a.IntRate * (1.0 + a.IntRate)^(1.0 - AnnuityStartTime + DiscountTo) )::Float64
        else  # case: finite Ordinary StdAnnuity
            return ( ( 1.0 - 1.0 / (1.0 + a.IntRate) ^ a.T ) / a.IntRate * (1.0 + a.IntRate)^(- AnnuityStartTime + DiscountTo) )::Float64
        end
    else  # case: infinite annuity
        if a.IsDue  # case: infinite StdAnnuity Due
            return ( (1.0 + a.IntRate)^(1.0 - AnnuityStartTime + DiscountTo) / a.IntRate )::Float64
        else  # case: infinite Ordinary StdAnnuity
            return ( (1.0 + a.IntRate)^(- AnnuityStartTime + DiscountTo) / a.IntRate )::Float64
        end
    end
    return -Inf::Float64  # formal return
end
# ----------------
"""
    PV( a::StdAnnuity ; AnnuityStartTime::Int = 0 )

computes the present value of given annuity; the value is that at time 0, the beginning of period 1;

Receives:
1. a::StdAnnuity, an StdAnnuity instance
2. AnnuityStartTime::Int, [in Z], an integer, negative or 0 or positive, the time/moment when Annuity starts, equals 0 by default;
Returns:
1. Res::Float64, the present value of given annuity

Note:
1. you may let AnnuityStartTime > 0 to define a Deferred StdAnnuity
"""
PV( a::StdAnnuity ; AnnuityStartTime::Int = 0 ) = AnnuityValue( a, AnnuityStartTime = AnnuityStartTime, DiscountTo = 0 )
# ----------------
"""
    FV( a::StdAnnuity ; AnnuityStartTime::Int = 0 )

computes the final value of given annuity; the value is that at the moment when the given annuity matures;

Receives:
1. a::StdAnnuity, an StdAnnuity instance
2. AnnuityStartTime::Int, [in Z], an integer, negative or 0 or positive, the time/moment when Annuity starts, equals 0 by default;
Returns:
1. Res::Float64, the present value of given annuity

Note:
1. the returned value is at the moment when the given annuity matures, not at time 0 (unless AnnuityStartTime == -1 * a.T)
2. FV() only works for finite annuity!
"""
FV( a::StdAnnuity ; AnnuityStartTime::Int = 0 ) = a.IsFinite  ?  AnnuityValue( a, AnnuityStartTime = AnnuityStartTime, DiscountTo = AnnuityStartTime + a.T )  :  error("final value is only defined on finite annuity")
# ----------------

# ==================================
# """
#     LifeTable
#
# type/template for life tables, it is defined from age 0 to age Xmax,
# it has members:
# 1. Xmax::Int, the maximum age to live, equals 105 by default
# 1. q::Vector{Float64}, q_{x}, mortality in digit, from age x to age x+1, a vector with length 101 (age 0 to age 100)
# """
# struct LifeTable
#     Xmax::Int  # the maximum age to live, equals 105 by default
#     q::Vector{Float64}  # q_{x}, mortality in digit, from age x to age x+1
#     L::Vector{Float64}  # L_{x},
#
#
# end  # end struct LifeTable
# ---------------

end  # end module Actuarial

max_stars_repo_name: Clpr/NumToolLibs

max_issues_repo_name: Clpr/NumToolLibs

max_forks_repo_name: Clpr/NumToolLibs


text:
mutable struct fixedspending
    #fields
end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
include("euler/euler.jl")
using .BigIntegers: BigInteger
using .Math: fast_max
using .Numbers: digits_sum
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    result = 0
    for a ∈ n - 5:n - 1
        number = BigInteger(a)
        for b ∈ 1:n - 1
            number *= BigInteger(a)
            result = fast_max(result, digits_sum(number))
        end
    end
    return result
end

compute(100)

@benchmark compute(100)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using Printf

using DataFrames
using CSV

using LinearAlgebra
BLAS.set_num_threads(1)

import MathProgBase
const MPB = MathProgBase

import Tulip

global const TIME_LIMIT = 300

# import solvers
include("solvers.jl")

# open-source solvers
const SOLVERS_OS = [sCLP(), sECOS(), sGLPK(), sIPOPT(), sTLP()]
const NAMES_OS   = solver_name.(SOLVERS_OS)

# commercials solvers
const SOLVERS_COMM = [sCPX(), sGRB(), sMSK()]
const NAMES_COMM   = solver_name.(SOLVERS_COMM)

"""
    benchmark(files, solvers names; verbose)

Benchmark the given solvers over the set of instances.
"""
function benchmark(
    dir,
    files,
    solvers;
    verbose=false
)

    names = solver_name.(solvers)

    time_tot = 0.0  # Total solve time

    n_err_read = 0
    n_err_solve = 0
    n_term = 0
    n_success = 0

    results = Dict()

    # Initial log
    verbose && @printf "%12s" ""
    for (s, name) in zip(solvers, names)
        results[name] = Dict{String, Dict}()
        verbose && @printf "%10s" name
    end
    verbose && @printf "\n"

    for f in files
        verbose && @printf "%12s" f

        # Read MPS file
        m, n, c, b, A, collb, colub, ranges = Tulip.readmps(dir * f)

        for (s, name) in zip(solvers, names)

            res = Dict()

            # Instanciate model
            model = MPB.LinearQuadraticModel(s)
            MPB.loadproblem!(model, A, collb, colub, c, b, b, :Min)

            res[:success] = false
            res[:fail] = false
            res[:error] = false
            res[:status] = :Error
            res[:time] = NaN
            res[:error_type] = "None"

            try
                t = @elapsed MPB.optimize!(model)
                s = MPB.status(model)

                res[:time] = t
                res[:status] = "$s"
                # res[:iter] = MPB.getbarrieriter(model)

                if s == :Optimal && t <= TIME_LIMIT
                    # Instance is solved
                    verbose && @printf "%10.3f" t
                    res[:success] = true
                else
                    # Error or timeout
                    verbose && @printf "%10s" "t"
                    res[:time] = Float64(TIME_LIMIT)         # count failure as max time.
                    res[:fail] = true
                end

            catch err
                verbose && @printf "%10s" "f"
                res[:error] = true
                res[:error_type] = "$(typeof(err))"
                res[:time] = Float64(TIME_LIMIT)
            end

            flush(stdout)

            results[name][f] = res

        end

        verbose && @printf "\n"
        flush(stdout)

    end

    return results
end

"""
    extract_results(res)

Parse results and output csv files.
"""
function extract_results(res)

    # Extract solvers
    solvers = collect(keys(res))

    # Extract problem names
    pbnames_raw = collect(keys(res[solvers[1]]))
    pbnames = lowercase.([f[1:end-4] for f in pbnames_raw])
    p = sortperm(pbnames)
    pbnames_raw = pbnames_raw[p]  # sort names
    pbnames = pbnames[p]

    # extract metrics
    metrics = collect(keys(res[solvers[1]][pbnames_raw[1]]))

    # Create data frame with all results
    df = DataFrame()

    df.instance = pbnames  # instance name
    for s in solvers
        for k in metrics
            df[Symbol("$(k)_$(s)")] = [res[s][f][k] for f in pbnames_raw]
        end
    end

    # convert boolean outputs to integers
    for s in solvers
        for k in metrics
            if eltype(df[Symbol("$(k)_$(s)")]) == Bool
                df[Symbol("$(k)_$(s)")] = Int.(df[Symbol("$(k)_$(s)")])
            end
        end
    end

    # Correct reported times
    for s in solvers
        t_ = df[Symbol("time_$s")]
        succ_ = df[Symbol("success_$s")]
        df[Symbol("time_$s")] .= min.(t_ , 300.0) .* succ_ .+ 300.0 .* (1.0 .- succ_)
    end

    return df
end


text:
# TESTED WITH [14b8a8f1] PkgTemplates v0.7.18

module PkgCreate

using Pkg
using Dates
using PkgTemplates

using ..DevEnvConfig: devEnvConfig_src_path

const StrOpt = Union{String,Nothing}
const BoolOpt = Union{Bool,Nothing}

function template(filename::String)
    joinpath(devEnvConfig_src_path(), "pkg", "templates", filename)
end

include("customize.jl")
using ..GitTools
using ..DevEnvConfig: Success, Warning, Error

# DOCUMENTED as `newpkg` in module DevEnvConfig
function create(pkname::String;
    private::Bool,
    dir::StrOpt=nothing,
    docrepo::StrOpt=nothing,
    useextjl::BoolOpt=nothing,
    generalregistry::BoolOpt=nothing,
    license::StrOpt=nothing,
    testing_gitconfig=nothing)

    # Default: a public repo may be able to register latter, a private repo /cannot/
    generalregistry = something(generalregistry, !private)
    # Default: private repo has NO .jl extension, public repo NEED it to register to General
    useextjl = something(useextjl, !private)

    gitconfif = if isnothing(testing_gitconfig); GitTools.getconfig(); else testing_gitconfig end
    tgtdir = something(dir, get(ENV, "JULIA_PKG_DEVDIR", "?"))
    warning = false

    if tgtdir == "?"
        @error "No target directory provided !
        Either define the environment variable JULIA_PKG_DEVDIR,
        or pass the directory with optional argument `dir`."
        return Error
    end

    try
        tgtdir = realpath(tgtdir)
    catch
        @error "The deployment path '$(tgtdir)' does not exist !"
        return Error
    end

    pkpath = joinpath(tgtdir, pkname)
    @info "The package $(pkname) will be created at $(pkpath)"
    if ispath(pkpath)
        @error "The path '$(pkpath)' already exist !"
        return Error
    end

    if generalregistry
        if private
            @error "Defined to both deploy into a private repo and register to the General registry.
            Only public package can be registred in General.
            Set `generalregistry` to `false` or `nothing` (default)."
            return Error
        end

        if !useextjl
            @error "Package was requested or inferred to be published into the General registry.
            The repo must use the .jl extension. See `generalregistry` and `!useextjl`."
            return Error
        end

        if !isnothing(docrepo)
            @error "Package was requested or inferred to be published into the General registry.
            The documentation should be deployed locally in the package repo,
            leave the `docrepo` argument to its default value of `nothing`"
            return Error
        end

        if isnothing(license)
            license = "MIT"
            @warn "Package was requested or inferred to be published into the General registry.
            The license was undefined. Forcing the license to `MIT`.
            An OSI-approved license is required to register the General registry.
            You may choose another OSI-approved license by using the `license` keyword argument."
            warning = true
        end
    end

    if gitconfif.def_branch != "main"
        @error "The local Git default branch should be configured to `main`.
        Execute: git config --global init.defaultBranch main"
        return Error
    end

    if private && isnothing(docrepo)
        @warn "This package is setup to deploy documentation to its own repo.
        This package is also setup to deploy to a private repo. If you have a
        free GitHub account you will NOT be able to enable Pages on this repo.
        You will thus not be able to access the documentation web site at all.
        If you want the repo of the package to be private, then first create
        another public GitHub repo to host the packages documentation and regenerate
        this package to deploy the documentation to the public doc repo."
        warning = true
    end

    pkgconfig = Template(;
        # authors="name <email>"            # Default value from the global Git config (user.name and user.email).
        user=gitconfif.user_github,         # hosting service username, default value from Git config (github.user)
        dir=tgtdir,                         # Directory to place created packages, default to JULIA_PKG_DEVDIR environment variable.
        host="github.com",                  # URL to the code hosting service where packages will reside.
        julia=v"1.6.0",                     # Minimum allowed Julia version for this package

        plugins = [

        !License,                           # no default license, we add it below if needed

        SrcDir(;                             # generate the file in the src directory
            file=template("module.jlt"),     # Custom template, more exemples
        ),

        ProjectFile(; version=v"0.1.0"),     # default Project.toml file

        Tests(;                              # create /test
            file=template("runtests.jlt"),   # Custom template, run doc-tests
            project=true),                   # always use subproject (test specific dependencies of Julia >1.2)

        Readme(; destination="README.md",    # README file. File destination, relative to the repository root.
            file=template("README.md"),      # We use a custome template
            inline_badges=false,             # IGNORED by the custom README template. Badge are never inline with the title.
            badge_order=[                    # Order of the README.md Badges generated by each plugin
                Documenter{GitPagesDoc},     # First the Dev/Stable doc links
                GitHubActions                # then the status of the CI Tests
            ]),

        Git(; name=gitconfif.user_local,     # override Git config user.name, allows for testing
            email=gitconfif.user_email,      # override Git config user.email, allows for testing
            branch="main",                   # The desired name of the repository's default branch.
            ssh=true,                        # Whether or not to use SSH for the remote. If left unset, HTTPS is used.
            jl=useextjl,                     # Whether or not to add a .jl suffix to the remote URL.
            manifest=false,                  # Whether or not to commit Manifest.toml.
            gpgsign=false ),                 # Whether or not to sign commits with your GPG key.

        CompatHelper(;
            file=template("CompatHelper.yml"), # Custom template
            destination="CompatHelper.yml",    # Destination of the workflow file, relative to .github/workflows
            cron="0 0 * * 0"),                 # Run weekly, default is "0 0 * * *" which run daily

        GitHubActions(;                      # If using coverage plugins, don't forget to manually add your API tokens as secrets https://docs.github.com/en/actions/reference/encrypted-secrets#creating-encrypted-secrets
            file=template("Tests.yml"),      # Use our custom template
            destination="Tests.yml",         # Destination of the workflow file, relative to .github/workflows.
            linux=true,                      # Whether or not to run builds on Linux.
            osx=false,                       # Whether or not to run builds on OSX (MacOS).
            windows=false,                   # Whether or not to run builds on Windows.
            x64=true,                        # Whether or not to run builds on 64-bit architecture.
            x86=false,                       # Whether or not to run builds on 32-bit architecture.
            coverage=true,                   # Publish code coverage? A coverage plugin such as Codecov must also be included.
            extra_versions=["1.6"]),         # Julia versions to test, default: ["1.0", "1.6", "nightly"]

        Codecov(;                          # Sets up code coverage submission from CI to Codecov.
            file=nothing),                 # Template file for .codecov.yml, or nothing to create no file.

        # Coveralls(;                      # Sets up code coverage submission from CI to Coveralls.
        #    file=nothing),                # Template file for .coveralls.yml, or nothing to create no file.

        Documenter{GitPagesDoc}(;               # Sets up documentation generation via Documenter.jl; custom GitPages deployment
            make_jl=template("docs_make.jlt"),  # Template file for make.jl.
            index_md=template("docs_index.md"), # Template file for index.md.
            # logo=Logo()),                     # A Logo containing documentation logo information.
            # canonical_url=...                 # This option is IGNORED by the custom GitPagesDoc
            # makedocs_kwargs::Dict{Symbol,Any} # Extra keyword arguments to be inserted into makedocs.
            devbranch="main",                   # Branch that will trigger docs deployment, or `nothing`
            assets=String[]),                   # Extra assets for the generated site.

        ! TagBot,
        # TagBot(;
        #     destination="TagBot.yml",        # Destination of the workflow file, relative to .github/workflows.
        #     trigger="JuliaTagBot",           # Username of the trigger user for custom regsitries.
        #     token=Secret("GITHUB_TOKEN"),    # Name of the token secret to use.
        #     ssh=Secret("DOCUMENTER_KEY"),    # Name of the SSH private key secret to use.
        #     ssh_password=nothing,            # Name of the SSH key password secret to use.
        #     changelog=nothing,               # Custom changelog template.
        #     changelog_ignore=nothing,        # Issue/pull request labels to ignore in the changelog.
        #     gpg=nothing,                     # Name of the GPG private key secret to use.
        #     gpg_password=nothing,            # Name of the GPG private key password secret to use.
        #     registry=nothing,                # Custom registry, in the format owner/repo.
        #     branches=nothing,                # Whether not to enable the branches option.
        #     dispatch=nothing,                # Whether or not to enable the dispatch option.
        #     dispatch_delay=nothing),         # Number of minutes to delay for dispatch events.

        # Logo(;                             # For Documenter{T} above. Logo information for documentation.
        #     light=nothing,                 # Path to a logo file for the light (default) theme.
        #     dark=nothing)                  # Path to a logo file for the dark theme.

        # BlueStyleBadge()                   # Adds a BlueStyle badge to the Readme file, see: https://github.com/invenia/BlueStyle
        # ColPracBadge()                     # Adds a ColPrac badge to the Readme file, see: https://github.com/SciML/ColPrac
        # Develop()                          # Adds generated packages to the current environment by deving them
        # Citation(;                         # Creates a CITATION.bib file for citing package repositories.
            # readme=false)                  # Whether or not to include a section about citing in the README.
        # RegisterAction(;                   # Add a GitHub Actions workflow for registering a package with the General registry
            # destination="register.yml",    # Destination of the workflow file, relative to .github/workflows.
            # prompt="Version?")             # Prompt for workflow dispatch.
    ]) # myTemplate = Template(; ...

    if generalregistry
        push!(pkgconfig.plugins, RegisterAction(;
            file=template("PublicRegister.yml"),
            destination="PublicRegister.yml",
            prompt="Version to register or component to bump"))
    end

    if isnothing(docrepo)
        push!(pkgconfig.plugins, GitPagesDoc(;
            privaterepo    = private,       # true if this package will be upstream-ed to a private GitHub repo
            deploy_gituser = nothing,       # IF nothing: deploy doc in the package repo and in gh-pages branch
            deploy_gitrepo = nothing,       # deploy doc into repo "github/ylvain/PkDoc"
            deploy_branch  = nothing,       # default to "gh-pages" for same repo
            deploy_subdir  = nothing))      # deploys at the root ("/") of the target repo filsystem
    else
        push!(pkgconfig.plugins, GitPagesDoc(;
            privaterepo    = private,       # true if this package will be upstream-ed to a private GitHub repo
            deploy_gituser = pkgconfig.user, # FIXME: assume the docrepo belong to the same GitHub user
            deploy_gitrepo = docrepo,       # deploy doc into repo "github/deploy_gituser/deploy_gitrepo"
            deploy_branch  = nothing,       # default to "main" for same a doc repo
            deploy_subdir  = "<PKG>"))      # "<PKG>" deploys at "/PkgRepoName" (name of the packge repo)
    end

    if !isnothing(license)
        push!(pkgconfig.plugins, License(;
            path=nothing,                  # Path to a custom license file. This keyword takes priority over name.
            name="MIT",                    # See https://github.com/invenia/PkgTemplates.jl/tree/master/templates/licenses
            destination="LICENSE"))        # for the list of open source license that can be given as `name`
    end

    pkgconfig(pkname)

    pkgconfig_env_test(pkpath)
    GitTools.commit(pkpath, "DevEnvConfig newpkg() extensions")

    @info "The package code repo is: $(hostuserrepo_code(pkgconfig, pkname))"
    @info "The package docs repo is: $(hostuserrepo_docs(pkgconfig, pkname))"
    @info "The local package is at $pkpath"
    return warning ? Warning : Success
end

end # module PkgCreate


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Interpreter for input files

"""
    Raytracer.Interpreter

SceneLang interpreter module.
"""
module Interpreter

using Base:
    SizeUnknown

using StaticArrays:
    SVector, @SVector, SOneTo

using Raytracer
using FileIO

import Base:
    show, showerror, print, eof, copy, IteratorSize

export # Interpreter
    SourceLocation,
    TokenValue,
        Keyword, Identifier, LiteralString, LiteralNumber, LiteralSymbol, MathExpression, StopToken,
    Token,
    InputStream,
        open_stream,
    Scene,
        parse_scene,
        parse_variables_from_string,
    InterpreterException

let interpreter_dir = "interpreter"
    include(joinpath(interpreter_dir, "sourcelocation.jl"))
    include(joinpath(interpreter_dir, "token.jl"))
    include(joinpath(interpreter_dir, "exceptions.jl"))
    include(joinpath(interpreter_dir, "utilities.jl"))
    include(joinpath(interpreter_dir, "inputstream.jl"))
    include(joinpath(interpreter_dir, "lexer.jl"))
    include(joinpath(interpreter_dir, "parser.jl"))
end

end # module


text:
# BAM
# ===
#
# The BAM file format.
#
# http://samtools.github.io/hts-specs/SAMv1.pdf

include("bai.jl")
include("record.jl")
include("reader.jl")
include("writer.jl")
include("auxdict.jl")
include("intersect.jl")


text:
function _make_scope_gui()

    popupmenu_scope = Menu()
    popupmenu_voltage=MenuItem("Voltage Gain")
    popupmenu_time=MenuItem("ms/div")
    popupmenu_thres=MenuItem("Threshold")
    popupmenu_signal=MenuItem("Signal")
    push!(popupmenu_scope,popupmenu_voltage)
    push!(popupmenu_scope,popupmenu_time)
    push!(popupmenu_scope,popupmenu_thres)
    push!(popupmenu_scope,popupmenu_signal)

    popupmenu_voltage_select=Menu(popupmenu_voltage)
    if VERSION > v"0.7-"
        scope_v_handles=Array{MenuItemLeaf}(undef,0)
    else
        scope_v_handles=Array{MenuItemLeaf}(0)
    end
    voltage_scales = [1, 50, 100, 200, 500]
    push!(scope_v_handles,MenuItem(string(voltage_scales[1])))
    push!(popupmenu_voltage_select,scope_v_handles[1])
    for i=2:5
        push!(scope_v_handles,MenuItem(scope_v_handles[i-1],string(voltage_scales[i])))
        push!(popupmenu_voltage_select,scope_v_handles[i])
    end

    popupmenu_time_select=Menu(popupmenu_time)
    if VERSION > v"0.7-"
        scope_t_handles=Array{MenuItemLeaf}(undef,0)
    else
        scope_t_handles=Array{MenuItemLeaf}(0)
    end
    time_scales = [1,2,3,4,5] #Need to make this ms / div
    for i=1:length(time_scales)
        time_scales[i]=round(time_scales[i],1)
    end
    push!(scope_t_handles,MenuItem(string(time_scales[1])))
    push!(popupmenu_time_select,scope_t_handles[1])
    for i=2:5
        push!(scope_t_handles,MenuItem(scope_t_handles[i-1],string(time_scales[i])))
        push!(popupmenu_time_select,scope_t_handles[i])
    end

    popupmenu_thres_select=Menu(popupmenu_thres)
    if VERSION > v"0.7-"
        scope_thres_handles=Array{MenuItemLeaf}(undef,0)
    else
        scope_thres_handles=Array{MenuItemLeaf}(0)
    end
    push!(scope_thres_handles,MenuItem("On"))
    push!(popupmenu_thres_select,scope_thres_handles[1])
    push!(scope_thres_handles,MenuItem(scope_thres_handles[1],"Off"))
    push!(popupmenu_thres_select,scope_thres_handles[2])

    popupmenu_signal_select=Menu(popupmenu_signal)
    if VERSION > v"0.7-"
        scope_signal_handles=Array{MenuItemLeaf}(undef,0)
    else
        scope_signal_handles=Array{MenuItemLeaf}(0)
    end
    push!(scope_signal_handles,MenuItem("Spike"))
    push!(popupmenu_signal_select,scope_signal_handles[1])
    push!(scope_signal_handles,MenuItem(scope_signal_handles[1],"LFP"))
    push!(popupmenu_signal_select,scope_signal_handles[2])

    Gtk.showall(popupmenu_scope)

    (scope_t_handles,scope_v_handles,scope_thres_handles,scope_signal_handles,popupmenu_scope)
end

function add_scope_callbacks(v_handles,t_handles,thres_handles,signal_handles,handles)

    for i=1:5
        signal_connect(scope_popup_v_cb,v_handles[i],"activate",Void,(),false,(handles,i-1))
    end

    for i=1:5
        signal_connect(scope_popup_t_cb,t_handles[i],"activate",Void,(),false,(handles,i-1))
    end

    for i=1:2
        signal_connect(scope_popup_thres_cb,thres_handles[i],"activate",Void,(),false,(handles,i-1))
    end

    for i=1:2
        signal_connect(scope_popup_signal_cb,signal_handles[i],"activate",Void,(),false,(handles,i-1))
    end

    nothing
end

function scope_popup_v_cb(widgetptr::Ptr,user_data::Tuple{Gui_Handles,Int64})

    han, event_id = user_data

    if event_id==0
        han.soft.v_div=1.0
    elseif event_id==1
        han.soft.v_div=50.0
    elseif event_id==2
        han.soft.v_div=100.0
    elseif event_id==3
        han.soft.v_div=200.0
    elseif event_id==4
        han.soft.v_div=500.0
    end

    han.soft.v_div /= 1000.0

    nothing
end

function scope_popup_t_cb(widgetptr::Ptr,user_data::Tuple{Gui_Handles,Int64})

    han, event_id = user_data

    #=
    500 pixels across and time boxes are 50 pixels
    So if Sample Rate is 30k, then each box is 50/30000 = 1.67 ms
    20k, 50/20000 = 2.5 ms
    =#

    if event_id==0
        han.soft.t_div=1.0
    elseif event_id==1
        han.soft.t_div=2.0
    elseif event_id==2
        han.soft.t_div=3.0
    elseif event_id==3
        han.soft.t_div=4.0
    elseif event_id==4
        han.soft.t_div=5.0
    end

    nothing
end

function scope_popup_thres_cb(w::Ptr,user_data::Tuple{Gui_Handles,Int64})

    han, event_id = user_data

    if event_id==0
        han.soft.thres_on=true
    else
        han.soft.thres_on=false
    end

    nothing
end

function scope_popup_signal_cb(w::Ptr,user_data::Tuple{Gui_Handles,Int64})

    han, event_id = user_data

    if event_id == 0
        han.soft.signal_type=1
    else
        han.soft.signal_type=2
    end

    nothing
end

function draw_scope(rhd::RHD2000,han::Gui_Handles)

    ctx=Gtk.getgc(han.c)
    myheight=height(ctx)

    if han.soft.draws>9

        startind=SAMPLES_PER_DATA_BLOCK*(han.soft.draws-1)+1
        if han.soft.signal_type==1
            for i=1:SAMPLES_PER_DATA_BLOCK
                han.soft.v[startind]=rhd.v[i,han.sc.spike]
                startind+=1
            end
        else
            for i=1:SAMPLES_PER_DATA_BLOCK
                han.soft.v[startind]=rhd.lfps[i,han.sc.spike]
                startind+=1
            end
        end

        #Paint over old line
        move_to(ctx,1.0,han.soft.last[1])
        for i=2:499
            line_to(ctx,i,han.soft.last[i])
        end

        #Paint over old asterisks
        move_to(ctx,1.0,myheight-4.0)
        line_to(ctx,512.0,myheight-4.0)
        move_to(ctx,1.0,myheight-7.0)
        line_to(ctx,512.0,myheight-7.0)

        set_source(ctx,han.soft.ctx)
        set_line_width(ctx,4.0)
        stroke(ctx)

        s=han.soft.v_div*-1

        #Draw voltage trace from desired channel
        move_to(ctx,1.0,myheight-150.0+han.soft.v[1]*s)
        han.soft.last[1]=myheight-150.0+han.soft.v[1]*s

        t_iter=floor(Int64,han.soft.t_div)

        spike_ind=1

        scope_ind=1+t_iter
        for i=2:512
            if spike_ind > 500
                spike_ind = 500
            end
            y=myheight-150.0+han.soft.v[scope_ind]*s
            line_to(ctx,i,y)
            han.soft.last[i]=y

            if scope_ind>han.soft.spikes[spike_ind]
                han.soft.prev_spikes[spike_ind]=i
                spike_ind+=1
            end

            scope_ind+=t_iter
        end
        set_source_rgb(ctx,1.0,1.0,1.0)
        set_line_width(ctx,0.5)
        stroke(ctx)

        #reset spikes
        for i=1:(spike_ind-1)
            move_to(ctx,han.soft.prev_spikes[i],myheight)
            show_text(ctx,"*")
        end
        han.soft.prev_num_spikes=spike_ind-1
        han.soft.num_spikes=0

        #draw threshold

        if han.soft.thres_on
            plot_thres_scope(han,rhd,ctx)
        end

        han.soft.draws=1
    else
        #copy voltage
        startind=SAMPLES_PER_DATA_BLOCK*(han.soft.draws-1)+1

        if han.soft.signal_type==1
            for i=1:SAMPLES_PER_DATA_BLOCK
                han.soft.v[startind]=rhd.v[i,han.sc.spike]
                startind+=1
            end
        else
            for i=1:SAMPLES_PER_DATA_BLOCK
                han.soft.v[startind]=rhd.lfps[i,han.sc.spike]
                startind+=1
            end
        end

        #get spikes
        for g=1:rhd.nums[han.sc.spike]
            han.soft.num_spikes+=1
            han.soft.spikes[han.soft.num_spikes]=rhd.buf[g,han.sc.spike].inds.start+SAMPLES_PER_DATA_BLOCK*(han.soft.draws-1)+1
        end

        han.soft.draws+=1
    end

    nothing
end

function plot_thres_scope(han,rhd,ctx)

    myheight=height(ctx)

    thres=han.sc.thres*han.soft.v_div

    move_to(ctx,1,myheight-150-thres+2)
    line_to(ctx,500,myheight-150-thres+2)

    move_to(ctx,1,myheight-150-thres-2)
    line_to(ctx,500,myheight-150-thres-2)

    set_line_width(ctx,5.0)
    set_source_rgb(ctx,0.0,0.0,0.0)
    stroke(ctx)

    move_to(ctx,1,myheight-150-thres)
    line_to(ctx,500,myheight-150-thres)
    set_line_width(ctx,1.0)
    set_source_rgb(ctx,1.0,1.0,1.0)
    stroke(ctx)

    nothing
end


text:
immutable NoncentralChisq{T<:Real} <: ContinuousUnivariateDistribution
    ν::T
    λ::T
    function NoncentralChisq(ν::T, λ::T)
        @check_args(NoncentralChisq, ν > zero(ν))
        @check_args(NoncentralChisq, λ >= zero(λ))
    	new(ν, λ)
    end
end

NoncentralChisq{T<:Real}(ν::T, λ::T) = NoncentralChisq{T}(ν, λ)
NoncentralChisq(ν::Real, λ::Real) = NoncentralChisq(promote(ν, λ)...)
NoncentralChisq(ν::Integer, λ::Integer) = NoncentralChisq(Float64(ν), Float64(λ))

@distr_support NoncentralChisq 0.0 Inf

#### Conversions

function convert{T <: Real, S <: Real}(::Type{NoncentralChisq{T}}, ν::S, λ::S)
    NoncentralChisq(T(ν), T(λ))
end
function convert{T <: Real, S <: Real}(::Type{NoncentralChisq{T}}, d::NoncentralChisq{S})
    NoncentralChisq(T(d.ν), T(d.λ))
end

### Parameters

params(d::NoncentralChisq) = (d.ν, d.λ)
@inline partype{T<:Real}(d::NoncentralChisq{T}) = T

### Statistics

mean(d::NoncentralChisq) = d.ν + d.λ
var(d::NoncentralChisq) = 2(d.ν + 2d.λ)
skewness(d::NoncentralChisq) = 2sqrt2*(d.ν + 3d.λ)/sqrt(d.ν + 2d.λ)^3
kurtosis(d::NoncentralChisq) = 12(d.ν + 4d.λ)/(d.ν + 2d.λ)^2

function mgf(d::NoncentralChisq, t::Real)
    exp(d.λ * t/(1 - 2t))*(1 - 2t)^(-d.ν/2)
end

function cf(d::NoncentralChisq, t::Real)
    cis(d.λ * t/(1 - 2im*t))*(1 - 2im*t)^(-d.ν/2)
end

### Evaluation & Sampling

@_delegate_statsfuns NoncentralChisq nchisq ν λ

rand(d::NoncentralChisq) = StatsFuns.RFunctions.nchisqrand(d.ν, d.λ)


text:
using DataFrames, TableOperations, Tables, Random
using TreeParzen
using Recommenders:
    Movielens100k,
    load_dataset,
    ratio_split,
    ImplicitMF,
    evaluate_u2i,
    MeanPrecision,
    MeanRecall,
    MeanNDCG

function main()
    ml100k = Movielens100k()
    download(ml100k)
    rating, _, _ = load_dataset(ml100k)

    Random.seed!(1234)
    train_valid_table, test_table = ratio_split(rating, 0.8)
    train_table, valid_table = ratio_split(train_valid_table, 0.8)

    prec10 = MeanPrecision(10)
    recall10 = MeanRecall(10)
    ndcg10 = MeanNDCG(10)
    metrics = [prec10, recall10, ndcg10]

    space = Dict(
        # :n_epochs => HP.Choice(:n_epochs, [1, 2]),
        :n_epochs => HP.Choice(:n_epochs, [32, 64, 128, 256]),
        :n_negatives => HP.QuantUniform(:n_negatives, 1.0, 16.0, 1.0),
        :learning_rate => HP.LogUniform(:learning_rate, log(1e-3), log(1.0)),
        :log2_dimension => HP.QuantUniform(:log2_dimension, 4.0, 9.0, 1.0),
        :use_bias => HP.Choice(:use_bias, [true, false]),
        :reg_coeff => HP.LogUniform(:reg_coeff, log(1e-3), log(1.0)),
    )

    function invert_output(params)
        @info params
        n_epochs = convert(Int, params[:n_epochs])
        n_negatives = convert(Int, params[:n_negatives])
        learning_rate = params[:learning_rate]
        dimension = convert(Int, 2^params[:log2_dimension])
        use_bias = params[:use_bias]
        reg_coeff = params[:reg_coeff]

        model = ImplicitMF(dimension, use_bias, reg_coeff)

        result = evaluate_u2i(
            model,
            train_table,
            valid_table,
            metrics,
            10,
            col_item = :movieid,
            n_epochs = n_epochs,
            n_negatives = n_negatives,
            learning_rate = learning_rate,
            drop_history = true,
        )
        @info result
        return -result[:ndcg10]
    end

    @info "Tuning start."
    best = fmin(invert_output, space, 100, logging_interval = -1)
    @info best

    @info "Evaluate best model."

    best_model =
        ImplicitMF(convert(Int, 2^best[:log2_dimension]), best[:use_bias], best[:reg_coeff])
    result = evaluate_u2i(
        best_model,
        train_valid_table,
        test_table,
        metrics,
        10,
        col_item = :movieid,
        n_epochs = convert(Int, best[:n_epochs]),
        n_negatives = convert(Int, best[:n_negatives]),
        learning_rate = best[:learning_rate],
        drop_history = true,
    )
    @info result

end

main()


text:
using HMMs
using Brobdingnag
using Dates

hmm = rand( HMMs.HMM{2,HMMs.GenTDist,Brob,Float64,Date} )
size( hmm.inequalityconstraintmatrix )

hmmtype = Models.FittableModel{Date, Float64, HMMs.HMM{statecount, HMMs.GenTDist, Brob, Float64, Date}, typeof(HMMs.em)}

using Models
using BusinessDays
using Dependencies
using Commas
using Distributed

startdate = Date(2004,7,2)
enddate = Date(2005,7,1)
statecount = 2
processes = 20

hmmtype = Models.FittableModel{Date, Float64, HMMs.HMM{statecount, HMMs.GenTDist, Brob, Float64, Date}, typeof(HMMs.em)}
criterion( model ) = HMMs.likelihood(model.model)[end]
modeldates = advancebdays( :USNYSE, startdate, -1 ):Year(1):enddate

    multistartmodeltype = Models.MultiStartModel{Date, Float64, hmmtype, typeof(criterion)}
    fitfunction = FunctionNode(Models.fit)
    singledatemodeltype = Models.FittableModel{Date, Float64, multistartmodeltype, typeof(fitfunction)}
    anmodeltype = Models.ANModel{Date, Float64, Models.RewindableModel{Date, Float64, singledatemodeltype}}
    modeltype = Models.AdaptedModel{Date, Float64, Models.LogReturnModel{Date, anmodeltype}}

    # now update it
    spx = Commas.readcomma( joinpath( homedir(), "data", "SPX" ) )

    index = 2
    dates = Date[]
    prices = Float64[]
    while true
        if spx.time[index-1] <= Time(16) && ( spx.time[index] > Time(16) || spx.date[index] > spx.date[index-1] )
            push!( dates, spx.date[index-1] )
            push!( prices, spx.price[index-1] )
        end
        if spx.date[index] >= startdate
            break
        end
global        index += 1
    end

    if processes > 1
        workers = addprocs(processes)
    end

    model = rand( modeltype, modelperiods=modeldates, seeds=1:100, t=dates[1:1], u=prices[1:1] )
#    model = rand( modeltype, modelperiods=modeldates, seeds=1:2, t=dates[1:1], u=prices[1:1] )
Models.update( model, dates[2:end], prices[2:end], modules=[:HMMs, :Brobdingnag, :Models] )
model.models[1].model.rootmodel.model.model.models[2].model.inequalityconstraintmatrix
msm = model.models[1].model.rootmodel.model.model
msm.models[msm.optimumindex].model.inequalityconstraintmatrix

    if processes > 1
        rmprocs(workers)
    end


text:
"""
    TypedFunction{T}

A callable object that takes an argument of type `T`. This is useful for limiting
dispatch to a function that requires a specific argument type.

See other example at https://github.com/PacktPublishing/Hands-on-Design-Patterns-and-Best-Practices-with-Julia/blob/b20ad0fdfe73b05091fce6e52ca72f03ee0746e2/Chapter12/2_variance.jl#L283
"""
struct TypedFunction{T}
    func::Function
end

# Dispatcher
(tf::TypedFunction{T})(arg::T; kwargs...) where {T} = tf.func(arg; kwargs...)


text:
module TestInteractive
using Debug

println(Debug.UI.helptext)
println()
println("Type an expression to evaluate it in the current scope.")

@debug begin
    function myfunc2(x::Matrix)
        @bp
        for i=1:size(x,1), j=1:size(x,2)
            println("x[$i,$j] = $(x[i,j])");
        end
    end
    myfunc2([1 2; 3 4])
end

end # module


text:
##
## A few handy methods missing in Nemo so far.

"""
$(SIGNATURES)

Iterator for the coefficients and exponent vectors of the given polynomial.

"""
function dissect(p::MPolyElem)
    return zip(coeffs(p), exponent_vectors(p))
end

get_coeffs(p, I=0:length(p)) = [coeff(p, i) for i in I]

# isneg_coeff(p, i) = isless(coeff(p, i), zero(base_ring(p)))
# ispos_coeff(p, i) = isless(zero(base_ring(p)), coeff(p, i))

_ispositive(x) = isless(zero(x), x)
_isnegative(x) = isless(x, zero(x))

# Add generic implementations (I need them for broadcasting)
AA.leading_coefficient(i::Number) = i

normalize(p) = divexact(p, content(p))

"""
$(SIGNATURES)

Given a Julia vector `V` of entries, construct the corresponding AbstractAlgebra.jl one-column matrix over the given ring `R`, assuming all the entries can be coerced into `R`.
"""
function Nemo.matrix(R, V::AbstractVector)
    return Nemo.matrix(R, reshape(V, (:, 1)))
end

"""
$(SIGNATURES)

Given a Julia vector `V` of entries, construct the corresponding AbstractAlgebra.jl diagonal matrix over the given ring `R`, assuming all the entries can be coerced into `R`.
"""
function diagonal(R, V::AbstractVector)
    n = size(V, 1)
    D = Nemo.zero_matrix(R, n, n)
    for i in 1:n D[i,i] = V[i] end
    return D
end

partialdervativeof(i::Integer) = p -> Nemo.derivative(p, i)

"""
$(SIGNATURES)

Returns the Jacobian matrix of a one-column matrix of polynomials `F` with respect to the generators of `R` indexed by `vars`. When `vars` is omitted all the generators of `R` are used.

# Examples
```julia
julia> using Nemo

julia> R, vars = PolynomialRing(ZZ, vcat(["k\$i" for i in 1:5], ["x\$i" for i in 1:4]));

julia> F = [vars[1]*vars[2]*vars[6]-vars[8]; vars[3]*vars[9]+2*vars[7]]
[k1*k2*x1-x3]
[ k3*x4+2*x2]

julia> Jacobian(R, F, 6:9)
[k1*k2  0  -1   0]
[    0  2   0  k3]

julia> Jacobian(R, F)
[k2*x1  k1*x1   0  0  0  k1*k2  0  -1   0]
[    0      0  x4  0  0      0  2   0  k3]
```
"""
function Jacobian(R::MPolyRing, F::AbstractVector, vars=1:nvars(R))
    rows = size(F, 1)
    cols = length(vars)
    J = zeros(R, rows, cols)
    for j in 1:cols
        J[:,j] = partialdervativeof(vars[j]).(F)
    end
    return J
end

"""
$(SIGNATURES)

Returns the `i`-th coefficient of the characteristic polynomial of `M`.
"""
function coeff_charpoly(M::MatElem, i)
    P, _ = PolynomialRing(base_ring(M), "x")
    return coeff(charpoly(P, M), i)
end
# # Examples
# ```jldoctest; setup = :(using SymbolicCRN, Nemo)
# julia> using Nemo


text:
path = dirname(@__FILE__)

# Initialize arguments to function
@load "$path/reference/kalman_filter_args.jld2" y T R C Q Z D E z0 P0

# Kalman Filter (all arguments and no presample)
out = kalman_filter(y, T, R, C, Q, Z, D, E, z0, P0)
@testset "Basic Kalman Filter (all arguments, no presample)" begin
    h5open("$path/reference/kalman_filter_out.h5", "r") do h5
        @test read(h5, "log_likelihood") ≈ sum(out[1])
        @test read(h5, "marginal_loglh") ≈ out[1]
        @test read(h5, "pred")           ≈ out[2]
        @test read(h5, "vpred")          ≈ out[3]
        @test read(h5, "filt")           ≈ out[4]
        @test read(h5, "vfilt")          ≈ out[5]
    end
end

# Kalman Filter (no initial conditions and no presample)
out = kalman_filter(y, T, R, C, Q, Z, D, E)
@testset "Kalman Filter (no initial conditions)" begin
    h5open("$path/reference/kalman_filter_out.h5", "r") do h5
        @test read(h5, "log_likelihood") ≈ sum(out[1])
        @test read(h5, "marginal_loglh") ≈ out[1]
        @test read(h5, "pred")           ≈ out[2]
        @test read(h5, "vpred")          ≈ out[3]
        @test read(h5, "filt")           ≈ out[4]
        @test read(h5, "vfilt")          ≈ out[5]
        @test z0                         ≈ out[6]
        @test P0                         ≈ out[7]
    end
end

# Kalman filter with presample
out = kalman_filter(y, T, R, C, Q, Z, D, E, Nt0=4)
@testset "Kalman Filter (presample)" begin
    h5open("$path/reference/kalman_filter_out_presample.h5", "r") do h5
        @test read(h5, "log_likelihood") ≈ sum(out[1])
        @test read(h5, "marginal_loglh") ≈ out[1]
        @test read(h5, "pred")           ≈ out[2]
        @test read(h5, "vpred")          ≈ out[3]
        @test read(h5, "filt")           ≈ out[4]
        @test read(h5, "vfilt")          ≈ out[5]
    end
end

# Initialize arguments to for multi-regime Kalman Filter (ZLB)
@load "$path/reference/kalman_filter_args_zlb.jld2" y Ts Rs Cs Qs Zs Ds Es regime_inds

out = kalman_filter(regime_inds, y, Ts, Rs, Cs, Qs, Zs, Ds, Es)
@testset "Kalman Filter (Multi-regime/ZLB)" begin
    h5open("$path/reference/kalman_filter_out_zlb.h5", "r") do h5
        @test read(h5, "log_likelihood") ≈ sum(out[1])
        @test read(h5, "marginal_loglh") ≈ out[1]
        @test read(h5, "pred")           ≈ out[2]
        @test read(h5, "vpred")          ≈ out[3]
        @test read(h5, "filt")           ≈ out[4]
        @test read(h5, "vfilt")          ≈ out[5]
    end
end

nothing


text:
function frac_enter_by_parental(ds :: DataSettings)
    # Read ROW totals b/c file is transposed
    load_fct =
        mt -> read_row_totals(raw_entry_gpa_parental(ds; momentType = mt));
    m, ses, cnts = choice_prob_from_xy(load_fct);
    # m, ses, cnts = choice_prob_from_row_total(ds, :fracEnter_pV);
    # m = read_row_totals(raw_file_path(ds, :fracEnter_pV));
    @assert all(m .> 0.0)  &&  all(m .< 1.0)
    @assert all(cnts .> 100)
    return m, ses, cnts
end

## Hours worked, year 1
function work_hours_by_parental(ds :: DataSettings)
    load_fct =
        mt -> read_col_totals(ds, :workTime_qpM, mt);
    m, ses, cnts = load_mean_ses_counts(load_fct);
    @assert all(m .> 300.0)  &&  all(m .< 2200.0)
    @assert all(cnts .> 100)
    return m, ses, cnts
end

# 4y colleges only
function time_to_grad_4y_by_parental(ds :: DataSettings)
    load_fct =
        mt -> read_col_totals(ds, :timeToGrad4y_pV, mt);
    m, ses, cnts = load_mean_ses_counts(load_fct);
    @assert all(m .< 6.0);
    @assert all(cnts .> 100);
    return m, ses, cnts
end

# ------------------------


text:
using aIFileIO
using Test

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


text:
@testset "size" begin
    poollength = 10
    horizon = 6
    datalength = 500
    input, target = get_data(:solar, poollength, datalength, horizon)
    @test size(input) == (137, poollength, 1, datalength)
    @test size(target) == (500,)
    @test_throws ArgumentError get_data(:solar, poollength, 100000, horizon)
end

@testset "data" begin
    poollength = 10
    horizon = 6
    datalength = 500
    for dataset in (:solar, :traffic, :exchange_rate, :electricity)
        input, target = get_data(dataset, poollength, datalength, horizon)
        @test size(input, 1) >= 1
    end
    @test_throws ArgumentError get_data(:mnist, poollength, datalength, horizon)
end

@testset "memory issues" begin
    data = ones(Float32, 100, 10)
    input, target = FluxArchitectures.prepare_data(data, 6, 70, 10, normalise=false)
    @test all(input .== 1.0)
    @test all(target .== 1.0)
end

@testset "Tables.jl interface" begin
    poollength = 10
    horizon = 6
    datalength = 500
    inp = FluxArchitectures.load_data(:solar) |> FluxArchitectures.Tables.table
    @test prepare_data(inp, poollength, datalength, horizon) == get_data(:solar, poollength, datalength, horizon)
end

max_stars_repo_name: sdobber/FluxArchitectures

max_issues_repo_name: sdobber/FluxArchitectures

max_forks_repo_name: sdobber/FluxArchitectures


text:
module TestReferences

using Test
using Mimi

import Mimi: model_params

@defcomp A begin
    p1 = Parameter()
    v1 = Variable(index = [time])
    function run_timestep(p, v, d, t)
        v.v1[t] = gettime(t)
    end
end

@defcomp B begin
    p1 = Parameter()
    p2 = Parameter(index = [time])
end

m = Model()
set_dimension!(m, :time, 10)
refA = add_comp!(m, A, :foo)
refB = add_comp!(m, B)

refA[:p1] = 3   # creates a parameter specific to this component, with name "foo_p1"
@test Mimi.get_model_param_name(m.md, :foo, :p1) == :foo_p1
@test :foo_p1 in keys(model_params(m))
@test Mimi.UnnamedReference(:B, :p1) in Mimi.nothing_params(m.md)

refB[:p1] = 5
@test Mimi.get_model_param_name(m.md, :B, :p1) == :B_p1
@test :B_p1 in keys(model_params(m))

# Use the ComponentReferences to make an internal connection
refB[:p2] = refA[:v1]

run(m)
@test m[:foo, :p1] == 3
@test m[:B, :p1] == 5
@test m[:B, :p2] == collect(1:10)

# Test `connect_param!` methods for ComponentReferences

connect_param!(refB, :p2, refA, :v1)
run(m)
@test m[:B, :p2] == collect(1:10)

@defcomp C begin
    v1 = Parameter(index = [time])
end
refC = add_comp!(m, C)
connect_param!(refC, refA, :v1)
run(m)
@test m[:C, :v1] == collect(1:10)

# test `update_param!` for references
refA[:p1] = 10
run(m)
@test m[:foo, :p1] == 10

end


text:
using Test

using BEAST
using CompScienceMeshes

o, x, y = point(0,0,0), point(1,0,0), point(0,1,0)
G = Mesh([o,x,y],[index(1,2,3)])

@test numvertices(G) == 3
@test numcells(G) == 1

K = MWDoubleLayerTDIO(1.0, 1.0, 0)

X = raviartthomas(G, skeleton(G,1))
Y = buffachristiansen(G, boundary(G))
@test numfunctions(X) == 3
@test numfunctions(Y) == 3

Δt = 20.0
Nt = 5
δ = timebasisdelta(Δt, Nt)
T = timebasiscxd0(Δt, Nt)
T2 = timebasisshiftedlagrange(Δt, Nt, 2)
V = X ⊗ T
W = Y ⊗ T

fr1, store1 = BEAST.allocatestorage(K, W, V,
    Val{:densestorage}, BEAST.LongDelays{:ignore})

BEAST.assemble!(K, W, V, store1)
Z = fr1()
@test all(==(0), Z[:,:,1])

W = X⊗δ
V = Y⊗T2

K = TDMaxwell3D.doublelayer(speedoflight=1.0, numdiffs=1)
fr2, store2 = BEAST.allocatestorage(K, W, V, Val{:densestorage}, BEAST.LongDelays{:ignore})
BEAST.assemble!(K, W, V, store2)
Z2 = fr2()
@test all(==(0), Z2[:,:,1])

γ = geometry(Y)
ch1 = chart(G, first(cells(G)))
qps = quadpoints(ch1, 3)

x = qps[1][1]
w = qps[1][2]

test_els = BEAST.elements(G)
trial_els = BEAST.elements(γ)
time_els = [(0.0, 20.0)]

qs = BEAST.defaultquadstrat(K,X,Y)
qdt = BEAST.quaddata(K, refspace(X), refspace(Y), refspace(T2), test_els, trial_els, nothing, qs)
qrl = BEAST.quadrule(K, refspace(X), refspace(Y), refspace(T2),
    1, test_els[1], 1, trial_els[1], 1, time_els[1], qdt, qs)
z = zeros(Float64, 3, 3, 10)
BEAST.innerintegrals!(z, K, x, refspace(X), refspace(Y), refspace(T2), test_els[1], trial_els[1], (0.0, 20.0), qrl, w)
@test_broken !any(isnan.(z))

verts = trial_els[1].vertices
import WiltonInts84
iG, ivG, igG = WiltonInts84.wiltonints(verts[1], verts[2], verts[3], cartesian(x), 0.0, 20.0, Val{2}, qrl.workspace)
@test_broken !any(isnan(iG))

ctr = WiltonInts84.contour!(verts[1],verts[2],verts[3],cartesian(x),0.0,20.0,qrl.workspace)
iG, ivG, igG = WiltonInts84.wiltonints(ctr,cartesian(x),Val{2})
@test_broken !any(isnan(iG))

import WiltonInts84
a = -0.23570226039551564
b = 1.5700924586837734e-16
p = -6.371760312437358e-16
h = 0.0
using StaticArrays
m = SVector(-0.7071067811865474, 0.7071067811865478, -0.0)
iG, vG = WiltonInts84.segintsg(a, b, p, h, m, Val{2})
@test_broken !any(isnan.(iG))