text:
type AmericanExercise <: EarlyExercise
  dates::Vector{Date}
end

type BermudanExercise <: EarlyExercise
  dates::Vector{Date}
end

type EuropeanExercise <: Exercise
  dates::Vector{Date}
end

EuropeanExercise(d::Date) = EuropeanExercise([d])
AmericanExercise(d1::Date, d2::Date) = AmericanExercise(Date[d1, d2])

type RebatedExercise{E <: Exercise, C <: BusinessCalendar, B <: BusinessDayConvention} <: Exercise
  exercise::E
  rebate::Float64
  rebateSettlementDays::Int
  rebatePaymentCalendar::C
  rebatePaymentConvention::B

  RebatedExercise(exercise::E, rebate::Float64 = 0.0, rebateSettlementDays::Int = 0, rebatePaymentCalendar::C = NullCalendar(),
                rebatePaymentConvention::B = Following()) = new{E, C, B}(exercise, rebate, rebateSettlementDays, rebatePaymentCalendar, rebatePaymentConvention)
end


text:
using WizardsChess
using Test

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

ext: jl


text:
using SignalDecomposition, Dates, Random, Plots
using ColorSchemes, Colors

pyplot()

Random.seed!(41516)
y = Date(2001):Day(1):Date(2025)
dy = dayofyear.(y)
cy =  @. 4 + 7.2cos(2π*dy/365.26) + 5.6cos(4π*dy/365.26 + 3π/5)
r0 = randn(length(dy))/2
sy = cy .+ r0

t = collect(1:length(y)) ./ 365.26 # true time in years
x, r = decompose(t, sy, Sinusoidal([1.0, 2.0]))

is = collect(0:0.02:1.0)
s1 = x .+ 6
s2 = r .- 10

cg1 = ColorScheme([parse(RGB, "black"), parse(RGB, "#745bb1")])
cg2 = ColorScheme([parse(RGB, "black"), parse(RGB, "#1aa0a4")])

append!(is, ones(100))
p1 = plot(t, sy, label = "input", color = "black")
annotate!(p1, [(1, 22.5, text("SignalDecomposition.jl", 48; family = "Bahnschrift"))])
ylims!(p1, -12, 25)
xlims!(0, 2)
plot!(p1, t, s1, label = "decomposed", alpha = 1.0)
plot!(p1; size = (1200, 600), showaxis = false, grid = false, legend=false)

# %%
anim = @animate for i in is
    p1 = plot(t, sy, label = "input", color = "black")
    x1 = @. i*s1 + (1-i)*sy
    x2 = @. i*s2 + (1-i)*sy
    plot!(p1, t, x1, label = "decomposed", alpha = 1.0, color = get(cg1, i))
    plot!(p1, t, x2, color = get(cg2, i))
    xlims!(0, 2)
    annotate!(p1, [(1, 22.5, text("SignalDecomposition.jl", 48; family = "Bahnschrift"))])
    ylims!(p1, -12, 25)
    plot!(p1; size = (1200, 600), showaxis = false, grid = false, legend=false)
end
gif(anim, "signaldecomposition.gif", fps=30)


text:
function fill_history(obj::BulletinAFilesLoader, arg0::IERSConventions_NutationCorrectionConverter, arg1::SortedSet)
    return jcall(obj, "fillHistory", void, (IERSConventions_NutationCorrectionConverter, SortedSet), arg0, arg1)
end


text:
import Base:haslength
import Serialization
import Knet

import Photon.Layers: Layer
import Photon.Losses: Loss

export Workout, saveworkout, loadworkout, predict, train!, hasmetric,
        freeze!, unfreeze!, validate, gradients, stop

# Callback niceties from Flux.jl
call_fn(f, xs...) = f(xs...)
runall(f) = f
runall(fs::AbstractVector) = () -> foreach(call_fn, fs)

"""
The Workout keeps track of the progress of the training session. At least a model
and a loss function needs to be provided. Optional an optimizer and one or more
metrics can be specified.

If no optimizer is specified, SGD will be used. If no metrics are provided, only
the loss during training and validation will be registered (:loss and :val_loss).

The provided mover will move data to the correct device. See also SmartMover. If
no mover is required, you can provide: (x) -> x or simple identity

# Usage

```julia
workout = Workout(model, L1Loss())

workout = Workout(model, CrossEntropy(), Adam())

workout = Workout(model, HingeLoss(); acc=BinaryAccuracy())

workout = Workout(model, L1Loss(), mover=identity)
```
"""
mutable struct Workout
    model::Layer
    loss::Union{Loss,Function}
    opt
    metrics::Vector{Pair}
    history::IdDict{Symbol,MetricStore}
    steps::Int
    epochs::Int
    mover::Union{Mover,Function}

    function Workout(model::Layer, loss::Union{Loss,Function},
                     opt=Knet.SGD(); mover=SmartMover(), metrics...)
        new(model, loss, opt, collect(metrics), IdDict(), 0, 0, mover)
    end
end

"""
Save a workout to a file. This will save all the state that is captured in the workout
and enables to continue at a later stage using the loadWorkout function. Under the hood
this function uses Julia serialization.
"""
function saveworkout(workout::Workout, filename="workout_$(workout.steps).dat")::String
    Serialization.serialize(filename, workout)
    return filename
end

"""
Load a workout from a file and return the initialized Workout.

# Usage

```julia
workout = loadWorkout("workout_1000.dat")
train!(workout, mydata)
```
"""
function loadworkout(filename)::Workout
    workout = Serialization.deserialize(filename)
    return workout
end

struct StopException <: Exception end

"""
Stop a training session. Typically invoked by a callback function that detects
that the training is not progressing anymore.

If this function is called outside the scope of a trianing session, an exception is thrown.
"""
stop(workout::Workout, reason::String) = throw(StopException())

"""
Checks if the workout has any recorded values for a certain metric. The provided
metricname has to be fully qualified.

# Usage

```julia
if hasmetric(workout, :val_loss) ...
```

"""
hasmetric(workout::Workout, metricname::Symbol)::Bool = haskey(workout.history, metricname)

"""
Update the workout history with a single metric value.
"""
function updatemetric!(workout::Workout, metricname::Symbol, value)
    e = get!(workout.history, metricname, SmartReducer())
    update!(e, workout.steps, value)
end

"""
Invoke the configured metrics. The loss metric will always be logged and available.
Metrics are stored in the history attribute of the workout.
"""
function updatemetrics!(workout::Workout, loss::Number, y, y_pred, phase=:train)

    # First register the loss
    metricname = getmetricname(:loss, phase)
    e = get!(workout.history, metricname, SmartReducer())
    update!(e, workout.steps, loss)

    # And now run and register any additional metrics
    for (name,fn) in workout.metrics
        try
            metricname = getmetricname(name, phase)
            metricvalue = fn(y_pred, y)
            updatemetric!(workout, metricname, metricvalue)
        catch
            @warn "Failed executing metric." metricname maxlog=1
        end
    end
    return loss
end

"""
Get the metric value for a fully qualified metric name and a certain step. If
step is not provided the last step will be used. If no value is found the passed
function will not be invoked.

# Usage

```julia
getmetricvalue(workout, :val_loss) do value
    println("validation loss", value)
end
```
"""
function getmetricvalue(f::Function, workout::Workout, metricname::Symbol, step=workout.steps)
    if haskey(workout.history, metricname)
        m = workout.history[metricname]
        value = get(m.state, step, nothing)
        value !== nothing && f(value)
    end
end

"""
Get the history of a metric. The provided metric has the fully qualified name
and the returned value is a tuple of steps and values.

```julia
h = history(workout, :val_loss)
# returns e.g ([1000, 2000, 3000, 4000], [0.81, 0.73, 0.64, 0.61])

h = history(workout, :loss)
```

"""
function history(workout::Workout, metric::Symbol)::Tuple
      h = workout.history[metric].state
      steps = sort(collect(keys(h)))
      return (steps, [h[step] for step in steps])
end

"""
Utility function to calculate the gradients. Useful when checking for vanishing or
exploding gradients. The returned value is a Vector of (Param, Gradient).
"""
function gradients(workout::Workout, minibatch)
    x, y = workout.mover(minibatch)
    J = Knet.@diff begin
        y_pred = workout.model(x)
        workout.loss(y_pred, y)
    end
    gradients = []
    for p in Knet.params(J)
        p.opt === :NOUPDATE && continue
        if p.opt === nothing; p.opt = Knet.clone(workout.opt); end
        push!(gradients, (p, Knet.grad(J,p)))
    end
    return gradients
end

"""
Freeze a parameter so it no longer will be updated during training.
"""
freeze!(p::Knet.Param) = p.opt = :NOUPDATE

"""
Unfreeze a parameter so it will be updated again during training.
"""
unfreeze!(p::Knet.Param)= p.opt = nothing

"""
Perform the back propagation and update of weights in a single go.
"""
function back!(J, opt)
    for p in Knet.params(J)
        p.opt === :NOUPDATE && continue
        if p.opt === nothing; p.opt = Knet.clone(opt); end
        Knet.update!(p, Knet.grad(J,p))
    end
end

"""
Take a single step in updating the weights of a model. This function
will be invoked from train! to do the actual learning.

For a minibatch (x,y) of data, the folowing sequence will be executed:

1. perform the forward pass
2. calculate the loss
3. update and remember the metrics, if any
4. do the backpropagation and update the weights
"""
function step!(workout::Workout, minibatch; zerograd=true)
    workout.steps += 1
    x, y = workout.mover(minibatch)
    J = Knet.@diff begin
        y_pred = workout.model(x)
        loss = workout.loss(y_pred, y)
        updatemetrics!(workout, Knet.value(loss), y, y_pred)
        loss
    end
    back!(J, workout.opt)
end

"""
Predict a sample, either a single value or a batch. Compared to invoking
the model directory with model(x), predit takes care of:

- Moving the data to the GPU if required.
- Shaping the data into a batch (controlled by makebatch parameter)

# Usage

```julia
x = randn(Float32, 224, 224, 3)
predict(workout, x)
```
"""
function predict(workout, x; makebatch=true)
    x = makebatch ? addlast(x) : x
    y = workout.model(workout.mover(x))
    makebatch ? droplast(y) : y
end

"""
Validate a minibatch and calculate the loss and metrics. Typically this function
is called from the train! method. But if required can also be invoked directly.
"""
function validate(workout::Workout, minibatch)
    x, y = workout.mover(minibatch)
    y_pred = workout.model(x)
    loss = workout.loss(y_pred, y)
    updatemetrics!(workout, loss, y, y_pred, :valid)
end

"""
Train the model based on a supervised dataset and for a number
of epochs. train! can be called multiple times and will continue
to train where is left of last time.

By default the train! function will try to ensure the data is of the right
type (e.g. Float32) and on the right device (e.g. GPU) before feeding it to
the model.

# Usage

```julia
train!(workout, traindata)
train!(workout, traindata, testdata, epochs=50)
```

"""
function train!(workout::Workout, data, validation=nothing;
    epochs=1, cb = nothing)

    cb = isnothing(cb) ? Photon.Callbacks.ConsoleMeter() : cb
    cb = runall(cb)

    for epoch in 1:epochs
        workout.epochs += 1
        d = data isa Function ? data() : data
        for minibatch in d
            step!(workout, minibatch)
            cb(workout, :train)
        end

        if validation !== nothing
            d = validation isa Function ? validation() : validation
            for minibatch in d
                validate(workout, minibatch)
            end
        end

        updatemetric!(workout, :epoch, workout.epochs)

        try
            cb(workout, :valid)
        catch ex
            ex isa StopException ? break : rethrow(ex)
        end
    end

end

@debug "Loaded Training module"


text:
# Decompressor Codec
# ==================

struct ZstdDecompressor <: TranscodingStreams.Codec
    dstream::DStream
end

function Base.show(io::IO, codec::ZstdDecompressor)
    print(io, summary(codec), "()")
end

"""
    ZstdDecompressor()

Create a new zstd decompression codec.
"""
function ZstdDecompressor()
    return ZstdDecompressor(DStream())
end

const ZstdDecompressorStream{S} = TranscodingStream{ZstdDecompressor,S} where S<:IO

"""
    ZstdDecompressorStream(stream::IO; kwargs...)

Create a new zstd decompression stream (`kwargs` are passed to `TranscodingStream`).
"""
function ZstdDecompressorStream(stream::IO; kwargs...)
    return TranscodingStream(ZstdDecompressor(), stream; kwargs...)
end

# Methods
# -------

function TranscodingStreams.initialize(codec::ZstdDecompressor)
    code = initialize!(codec.dstream)
    if iserror(code)
        zstderror(codec.dstream, code)
    end
    return
end

function TranscodingStreams.finalize(codec::ZstdDecompressor)
    if codec.dstream.ptr != C_NULL
        code = free!(codec.dstream)
        if iserror(code)
            zstderror(codec.dstream, code)
        end
        codec.dstream.ptr = C_NULL
    end
    return
end

function TranscodingStreams.startproc(codec::ZstdDecompressor, mode::Symbol, error::Error)
    code = reset!(codec.dstream)
    if iserror(code)
        error[] = ErrorException("zstd error")
        return :error
    end
    return :ok
end

function TranscodingStreams.process(codec::ZstdDecompressor, input::Memory, output::Memory, error::Error)
    dstream = codec.dstream
    dstream.ibuffer.src = input.ptr
    dstream.ibuffer.size = input.size
    dstream.ibuffer.pos = 0
    dstream.obuffer.dst = output.ptr
    dstream.obuffer.size = output.size
    dstream.obuffer.pos = 0
    code = decompress!(dstream)
    Δin = Int(dstream.ibuffer.pos)
    Δout = Int(dstream.obuffer.pos)
    if iserror(code)
        error[] = ErrorException("zstd error")
        return Δin, Δout, :error
    else
        return Δin, Δout, code == 0 ? :end : :ok
    end
end

function TranscodingStreams.expectedsize(codec::ZstdDecompressor, input::Memory)
    ret = find_decompressed_size(input.ptr, input.size)
    if ret == ZSTD_CONTENTSIZE_ERROR
        # something is bad, but we ignore it here
        return Int(input.size)
    elseif ret == ZSTD_CONTENTSIZE_UNKNOWN
        # random guess
        return Int(input.size * 2)
    else
        # exact size
        return Int(ret)
    end
end


text:
module LogicalKoans

function not(x)
end

function equals(x, y)
end

function not_equals(x, y)
end

function less_than(x, y)
end

function greater_than(x, y)
end

function less_than_or_equal_to(x, y)
end

function greater_than_or_equal_to(x, y)
end

end


text:
abstract type MeanFunction end

"""
    ZeroMean{T<:Real} <: MeanFunction

Returns `zero(T)` everywhere.
"""
struct ZeroMean{T<:Real} <: MeanFunction end

Base.map(::ZeroMean{T}, x::AbstractVector) where T = zeros(T, length(x))

ZeroMean() = ZeroMean{Float64}()

"""
    ConstMean{T<:Real} <: MeanFunction

Returns `c` everywhere.
"""
struct ConstMean{T<:Real} <: MeanFunction
    c::T
end

Base.map(m::ConstMean, x::AbstractVector) = fill(m.c, length(x))

"""
    CustomMean{Tf} <: MeanFunction

A wrapper around whatever unary function you fancy. Must be able to be mapped over an
`AbstractVector` of inputs.
"""
struct CustomMean{Tf} <: MeanFunction
    f::Tf
end

Base.map(f::CustomMean, x::AbstractVector) = map(f.f, x)


text:
using DataFrames: Highlighter
using StatsBase: self_avoid_sample!
# using ScottishTaxBenefitModel
using Test
using ScottishTaxBenefitModel

include( "testutils.jl")

include( "test_utils_tests.jl")

include( "randoms_tests.jl")

include( "general_tests.jl")

include( "results_tests.jl")

include( "equivence_scale_tests.jl")

include( "historic_benefits_tests.jl")

include("non_means_tested_bens_tests.jl")

include("benefit_generosity_tests.jl")

include( "income_tax_tests.jl")
include( "matching_tests.jl" )
include( "parameter_tests.jl")
include( "ni_tests.jl")

include( "legacy_mt_tests.jl")
include( "complete_mt_bens_tests.jl")

include( "complete_calc_tests.jl")
include( "uprating_tests.jl")
include( "social_security_age_tests.jl")
include( "minimum_wage_tests.jl")
include( "local_level_calculations_tests.jl" )
include( "scottish_benefits_tests.jl" )
include( "ubi_tests.jl")
include( "output_tests.jl")

include( "vs_policy_in_practice_tests.jl")
include( "vs_age_uk_tests.jl")

# These will only run if datasets are locally installed
if IS_LOCAL
	# These will only run if datasets are locally installed
    include( "household_tests.jl")
	include( "weighting_tests.jl")
    include( "uc_transition_tests.jl")
    include( "simple_runner_tests.jl")
    include( "equaliser_tests.jl")
end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
 #=

Implementations of pcg that use BLAS when they can.
Computes error in matrix norm, if correct left hand side is present.

Started by Dan Spielman.
Contributors: Serban Stan

=#

"""
`pcg(mat, b, pre; tol, maxits, maxtime, verbose)` solves a symmetric linear system using preconditioner `pre`.
`pre` can be a function or a matrix.  If a matrix, a function to solve it is created with cholFact.
`tol` is set to 1e-6 by default,
`maxits` defaults to Inf
`maxtime` defaults to Inf.  It measures seconds.
`verbose` defaults to false
`lhs` defaults to the all zeros vector
"""

function pcgMatNorm{Tv}(mat, b::Array{Tv,1}, pre::Function, lhs::Array{Tv,1}; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)
    return pcgBLASMatNorm(mat, b, pre, lhs, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)
end

#====================================
   PCG Code
=====================================#

# uses BLAS.  As fast as Matlab's pcg.
# set to use similar paramaters
function pcgBLASMatNorm{Tval}(mat, b::Array{Tval,1}, pre, lhs::Array{Tval,1}; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)

    n = size(mat,2)

    x = zeros(Tval,n)

    r = copy(b)
    z = pre(r)
    p = copy(z)

    rho = dot(r, z) # BLAS.dot does not seem faster

    t1 = time()

    itcnt = 0
    while itcnt < maxits
        itcnt = itcnt+1

        q = mat*p

        al = rho/dot(p, q)

        BLAS.axpy!(al,p,x)  # x = x + al * p
        BLAS.axpy!(-al,q,r)  # r -= al*q

        # This will take some time. Can we do something faster?
        if sqrt((lhs - x)' * mat * (lhs - x))[1] < tol
            break
        end

        # here is the top of the code in numerical templates

        z = pre(r)

        oldrho = rho
        rho = dot(z, r)

        # the following would have higher accuracy
        #       rho = sum(r.^2)

        beta = rho/oldrho

        p = z + beta*p
        # BLAS.scal!(n,beta,p,1) # p *= beta
        # BLAS.axpy!(1.0,z,p) # p += z

        if (time() - t1) > maxtime
            if verbose
                println("PCG stopped at maxtime.")
            end
            break
        end

      end

    if verbose
        println("PCG stopped after: ", itcnt, " iterations with relative error ", sqrt((lhs - x)' * mat * (lhs - x))[1], ".")
    end

    return x
end


text:
function _get_extension(file_name)
    ext = splitext(file_name)[end][2:end]
    ext = ext == "hdf5" ? "h5" : ext
    !(ext in ["csv", "h5"]) && error("File type $ext is not supported.")
    ext
end

_check_file(file_path) = !isfile(file_path) && error("File: $file_path does not exist.")
_get_function(operation, ext) = getfield(DM, Symbol(operation * ext))

function convert_recursive_array(a)
    if eltype(a) <: AbstractArray
        data = [convert_recursive_array(x) for x in a]
        d = ndims(data[1]) + 1
        cat(data..., dims=d)
    else
        a
    end
end

"""
    parse_string_to_num(d)

Parse the strings in a dictionary `d` to numbers.
```
"""
function parse_string_to_num(d)
    res = []
    for (k, v) in d
        if isa(v, String)
            pv = eval(Meta.parse(v))
        else
            pv = v
        end
        push!(res, k => pv)
    end
    Dict(res...)
end


text:
using Statistics
using MATLAB
include("atn_module.jl")
# using .time_series

N = 10_000          # number of nodes
m = 30               # number of active links per active node
mu = 0.1            # probability I -> S
tmax = 1000          # maximum duration of the simulation
Ni = Int(0.1 * N)   # initial number of infects
gamma = -2.22       # exponent of the power law
NTrials = 5         # number of trials for each value of lambda
tolerance = 1e-15   # termination tolerance
TWindow = 1000      # time window for evaluating steady state
TMinConvergence = 500

# health state of agents S=0 I=1

# diffferent lambda
# lambda=0.1:0.1:1;
# IAvgLambda=zeros(length(lambda))
# beta=zeros(length(lambda))

lambda = 0.2

ITrend = atn_module.time_series(Ni, N, m, gamma, lambda, mu, tolerance, TWindow,
                                TMinConvergence, tmax);

mat"plot($ITrend)"


text:
using NeRCA
using Optim

track = NeRCA.Track([0, 0, 0.2], [0, 200, 300])
ccalc = NeRCA.make_cherenkov_calculator(track)
hit_positions = [Position(0, 0, 0),
                 Position(20, 50, 100),
                 Position(100, -40, 400),
                 Position(23, 56, 50)]
hit_positions = [Position(0, 0, z) for z in 1:20:800]
times = ccalc.(hit_positions)

function make_quality_function(positions, times)
    function quality_function(params)
        x, y, z, dx, dy, dz = params
        track = NeRCA.Track([dx, dy, dz], [x, y, z])
        ccalc = NeRCA.make_cherenkov_calculator(track)
        expected_times = ccalc.(positions)
        return sum((times - expected_times).^2)
    end
    return quality_function
end

qfunc = make_quality_function(hit_positions, times)
qfunc([100,200,300,0,1,0])

function make_single_du_quality_function(positions, times)
    function quality_function(params)
        z, d, θ = params
        track = NeRCA.Track([0, sin(θ), cos(θ)], [0, d, z])
        ccalc = NeRCA.make_cherenkov_calculator(track)
        expected_times = ccalc.(positions)
        return sum((times - expected_times).^2)
    end
    return quality_function
end

qfunc = make_single_du_quality_function(hit_positions, times)
qfunc([300, 200, 0.2])


text:
using SymPy

# National Income Model

eq1 =  SymFunction("eq1")
eq2 =  SymFunction("eq2")

Y, C, b, I₀, G₀, a = symbols("Y C b I₀ G₀, a")

eq1 = Eq(Y - C, I₀ + G₀)
eq2 = Eq(-b*Y + C, a)

display(eq1)
display(eq2)

row1 = [1 -1 I₀ + G₀]
row2 = [-b 1 a]

system  = sympy.Matrix((row1, row2))

display(system)

solution = linsolve(system, (Y, C))
display(solution)

# Demand and Supply Model

eq1 = SymFunction("eq1")
eq2 = SymFunction("eq2")
eq3 = SymFunction("eq3")

Qᵈ, Qˢ, a, b, P, f, Y, c, d, g, Pᵢ = symbols("Qᵈ Qˢ a b P f Y c d g Pᵢ")

eq1 = Eq(Qᵈ + b*P, a + f*Y)
eq2 = Eq(Qˢ + d*s, c - g*Pᵢ)
eq3 = Eq(Qᵈ - Qˢ)

row1 = [1 0 b a+f*Y]
row2 = [0 1 -d c-g*Pᵢ]
row3 = [1 -1 0 0]

system = sympy.Matrix(((1, 0, b, a+f*Y), (0, 1, -d, c-g*Pᵢ), (1, -1, 0, 0)))

solution = linsolve(system, (Qᵈ, Qˢ, P))

display(solution)

sympy.init_printing(use_latex = "mathjax")


text:
using DistributionsAD

using Turing, Random, Test, LinearAlgebra
using Turing: Variational
using Turing.Variational: TruncatedADAGrad, DecayedADAGrad, AdvancedVI

include("../test_utils/AllUtils.jl")

@testset "advi.jl" begin
    @turing_testset "advi constructor" begin
        Random.seed!(0)
        N = 500

        s1 = ADVI()
        q = vi(gdemo_default, s1)
        c1 = rand(q, N)
    end
    @numerical_testset "advi inference" begin
        @testset for opt in [TruncatedADAGrad(), DecayedADAGrad()]
            Random.seed!(1)
            N = 500

            alg = ADVI(10, 5000)
            q = vi(gdemo_default, alg; optimizer = opt)
            samples = transpose(rand(q, N))
            chn = Chains(reshape(samples, size(samples)..., 1), ["s", "m"])

            # TODO: uhmm, seems like a large `eps` here...
            check_gdemo(chn, atol = 0.5)
        end
    end

    @turing_testset "advi different interfaces" begin
        Random.seed!(1234)

        target = MvNormal(ones(2))
        logπ(z) = logpdf(target, z)
        advi = ADVI(10, 1000)

        # Using a function z ↦ q(⋅∣z)
        getq(θ) = TuringDiagMvNormal(θ[1:2], exp.(θ[3:4]))
        q = vi(logπ, advi, getq, randn(4))

        xs = rand(target, 10)
        @test mean(abs2, logpdf(q, xs) - logpdf(target, xs)) ≤ 0.07

        # OR: implement `update` and pass a `Distribution`
        function AdvancedVI.update(d::TuringDiagMvNormal, θ::AbstractArray{<:Real})
            return TuringDiagMvNormal(θ[1:length(q)], exp.(θ[length(q) + 1:end]))
        end

        q0 = TuringDiagMvNormal(zeros(2), ones(2))
        q = vi(logπ, advi, q0, randn(4))

        xs = rand(target, 10)
        @test mean(abs2, logpdf(q, xs) - logpdf(target, xs)) ≤ 0.05
    end
end


text:
"""
    EffectSizes

A Julia package for effect size measures.
"""
module EffectSizes

export
    EffectSize,
    CohenD,
    HedgeG,
    GlassΔ,
    effectsize,
    confint,
    lower,
    upper,
    quantile

import Statistics: mean, std
import Distributions
import Distributions: quantile, Normal
import HypothesisTests
import HypothesisTests: confint
import Base.Grisu: PRECISION

include("confint.jl")
include("effectsize.jl")

"""
    const EffectSize = CohenD

See [`CohenD`](@ref).
"""
const EffectSize = CohenD

end # module


text:
using GeometricIntegratorsDiffEq
using Test

using DiffEqProblemLibrary
using GeometricIntegrators

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_2Dlinear

prob = prob_ode_2Dlinear

sol = solve(prob,GIEuler(),dt=0.1)
sol = solve(prob,GIMidpoint(),dt=0.1)
sol = solve(prob,GIHeun(),dt=0.1)
sol = solve(prob,GIKutta(),dt=0.1)
sol = solve(prob,GIERK4(),dt=0.1)
sol = solve(prob,GIERK438(),dt=0.1)
sol = solve(prob,GIImplicitEuler(),dt=0.1)
sol = solve(prob,GIImplicitMidpoint(),dt=0.1)
sol = solve(prob,GIRadIIA2(),dt=0.1)
sol = solve(prob,GIRadIIA3(),dt=0.1)
sol = solve(prob,GISRK3(),dt=0.1)
sol = solve(prob,GIGLRK(2),dt=0.1)

u0 = zeros(2)
v0 = ones(2)
f2 = function (dv,v,u,t)
  dv .= -u
end
function (::typeof(f2))(::Type{Val{:analytic}}, y0, p, x)
  u0, v0 = y0
  ArrayPartition(u0*cos(x) + v0*sin(x), -u0*sin(x) + v0*cos(x))
end
prob = SecondOrderODEProblem(f2,u0,v0,(0.0,5.0))

sol = solve(prob,GISymplecticEulerA(),dt=0.1)
sol = solve(prob,GISymplecticEulerB(),dt=0.1)
sol = solve(prob,GILobattoIIIAIIIB2(),dt=0.1)
sol = solve(prob,GILobattoIIIBIIIA2(),dt=0.1)


text:
export Physics
export addRigidBox!, addStaticBox!, addTriggerBox!, removePhysicsObject!
export other

other(e::Entity, col::TS_CollisionEvent) = (e.id == col.id1) ? col.id2 : col.id1

mutable struct PhysicsObjectInfo
  hx::AbstractFloat
  hy::AbstractFloat
  hz::AbstractFloat
  m::AbstractFloat
  isKinematic::Bool
  mx::AbstractFloat
  my::AbstractFloat
  mz::AbstractFloat
  PhysicsObjectInfo(hx = 0, hy = 0, hz = 0, m = 0, isKinematic = false, mx = 0, my = 0, mz = 0) = new(hx, hy, hz, m, isKinematic, mx, my, mz)
end

mutable struct Physics
  ids::Dict{Number, PhysicsObjectInfo}
  Physics() = new(Dict{Number, PhysicsObjectInfo}())
end

function addRigidBox!(e, hx, hy, hz, m, px, py, pz, isKinematic = false, mx = 0, my = 0, mz = 0)
  phys().ids[e.id] = PhysicsObjectInfo(hx, hy, hz, m, isKinematic, mx, my, mz)
  TS_BtAddRigidBox(e.id, hx + mx, hy + my, hz + mz, m, px, py, pz, isKinematic)
end

function addStaticBox!(e, hx, hy, hz, px, py, pz, mx = 0, my = 0, mz = 0)
  phys().ids[e.id] = PhysicsObjectInfo(hx, hy, hz, 0.0, false, mx, my, mz)
  TS_BtAddStaticBox(e.id, hx + mx, hy + my, hz + mz, px, py, pz)
end

function addTriggerBox!(e, hx, hy, hz, px, py, pz, mx = 0, my = 0, mz = 0)
  phys().ids[e.id] = PhysicsObjectInfo(hx, hy, hz, 1.0, false, mx, my, mz)
  TS_BtAddTriggerBox(e.id, hx + mx, hy + my, hz + mz, px, py, pz)
end

function removePhysicsObject!(e)
  delete!(phys().ids, e.id)
  TS_BtRemovePhysicsObject(e.id)
end

function handleMessage!(p::Physics, m::TICK)
  Log.@debug "Physics tick"
  TS_BtStepSimulation()
  for (id, pinfo) in p.ids
    pos = TS_BtGetPosition(id)
    # TODO: this is a bad hack for working specifically with rectangular objects, needs improvement
    getEntityById(id).pos = XYZ(pos.x - pinfo.hx, pos.y - pinfo.hy, pos.z - pinfo.hz)
  end
  while true
    col = TS_BtGetNextCollision()
    if col.id1 == -1 && col.id2 == -1 break end
    Log.@debug "entities $(col.id1) and $(col.id2) have a collision event"
    e1 = getEntityById(col.id1)
    e2 = getEntityById(col.id2)
    handleMessage!(e1, col)
    handleMessage!(e2, col)
  end
end

function awake!(p::Physics)
  Log.@debug "Physics awake!"
  listenFor(p, TICK)
end

function shutdown!(p::Physics)
  Log.@debug "Physics shutdown!"
  unlistenFrom(p, TICK)
end


text:
using Plots
using Plots.PlotMeasures
using Colors
using StatsPlots

include("collect_data.jl")

function plot_D()

    mcolor = ARGB(186/255,12/255,47/255, 0.8)
    pcolor = ARGB(39/255, 56/255, 150/255, 0.75)
    ocolor = :darkolivegreen

    mkl   = 1000 .* raw_average(:fermi, "MKL")
    oct   = 1000 .* raw_average(:fermi, "Octavian")
    psi4  = 1000 .* raw_average(:psi4)

    nao = [58, 48, 104, 114, 264, 251, 321, 68, 96, 148, 228, 208, 264, 275, 321, 86, 138, 143, 147, 228, 275, 256]
    perm = sortperm(nao)
    nao = nao[perm]

    p1 = groupedbar([oct mkl psi4], c = [ocolor mcolor pcolor], legend=:topright,
                label=["Fermi - Octavian" "Fermi - MKL" "Psi4"], dpi=1200, orientation = :v,
                bar_width=0.8,xticks=([2, 6, 12, 18, 22],["S2","S6", "S12","S18","S22"]), framestyle=:box)

    ylabel!("Time (seconds)")
    xlabel!("Dimer")
    p1
    p2 = scatter(nao, mkl[perm] ./ psi4[perm], label="", yticks=[0.5, 1.0, 1.5], ylims=[0.3,1.8],
                color=:black, ylabel="Relative time\n(Fermi/Psi4)", framestyle=:box, marker=:diamond, markercolor=mcolor,
                xlabel="Number of basis functions", legend=:bottomleft)

    scatter!(p2, nao, oct[perm] ./ psi4[perm], label="", markercolor = ocolor, marker=:diamond)
    hline!(p2, [1], color=:gray, linestyle=:dot, label="")

    plot(p1, p2, layout = grid(2,1, heights=[0.7, 0.3]), leftmargin=0.8cm)
end

function plot_E()

    println("Plot (T) data")
    mcolor = ARGB(186/255,12/255,47/255, 0.8)
    pcolor = ARGB(39/255, 56/255, 150/255, 0.75)
    ocolor = :palegreen4

    # MKL is actually tblis in this case... ops

    mkl   = 1000 .* raw_average(:fermi, "MKL")
    oct   = 1000 .* raw_average(:fermi, "Octavian")
    psi4  = 1000 .* raw_average(:psi4)

    nao = [58, 48, 104, 114, 264, 251, 321, 68, 96, 148, 228, 208, 264, 275, 321, 86, 138, 143, 147, 228, 275, 256]

    perm = sortperm(nao)
    nao = nao[perm]

    p1 = groupedbar([oct mkl psi4], c = [ocolor mcolor pcolor], legend=:topright,
                label=["Fermi - Octavian" "Fermi - TBLIS" "Psi4"], dpi=1200, orientation = :v,
                yaxis=:log,
                ylim = (10^1,10^9),
                bar_width=0.8,xticks=([2, 6, 12, 18, 22],["S2","S6", "S12","S18","S22"]), framestyle=:box)

    ylabel!("Time (miliseconds)")
    xlabel!("Dimer")
    p1
    p2 = scatter(nao, mkl[perm] ./ psi4[perm], label="TBLIS", yticks=[0.5, 1.0, 1.5], ylims=[0.3,1.8],
                color=:black, ylabel="Relative time\n(Fermi/Psi4)", framestyle=:box, marker=:diamond, markercolor=mcolor,
                xlabel="Number of basis functions", legend=:bottomleft, legendfontsize=6, legendfonthalign=:left)

    scatter!(p2, nao, oct[perm] ./ psi4[perm], label="Octavian", markercolor = ocolor, marker=:diamond, legend=:topright)
    hline!(p2, [1], color=:gray, linestyle=:dot, label="")

    plot(p1, p2, layout = grid(2,1, heights=[0.7, 0.3]), leftmargin=0.8cm, dpi=1200)

end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
struct LCMControlReceiver
    result::DynamicsResult{Float64, Float64}
    tprev::Base.RefValue{Float64}
    τprev::Vector{Float64}
    robot_info::HumanoidRobotInfo{Float64}
    lcm::LCM
    robot_state_channel::String
    robot_state_msg::robot_state_t
    atlas_command_msg::atlas_command_t
    encodebuffer::FastWriteBuffer
    decodebuffer::FastReadBuffer
    new_command::Base.RefValue{Bool}
end

# TODO: reduce cost of encoding

function LCMControlReceiver(robot_info::HumanoidRobotInfo;
        robot_state_channel::String = "EST_ROBOT_STATE",
        robot_command_channel::String = "ROBOT_COMMAND")
    mechanism = robot_info.mechanism
    result = DynamicsResult(mechanism)
    tprev = Ref(0.)
    τprev = zeros(num_velocities(mechanism))
    lcm = LCM()
    robot_state_msg = robot_state_t()
    atlas_command_msg = atlas_command_t()
    encodebuffer = FastWriteBuffer()
    decodebuffer = FastReadBuffer()
    receiver = LCMControlReceiver(result, tprev, τprev, robot_info, lcm, robot_state_channel,
        robot_state_msg, atlas_command_msg, encodebuffer, decodebuffer, Ref(false))

    for joint in tree_joints(mechanism)
        num_velocities(joint) == 1 || continue
        push!(robot_state_msg.joint_name, joint.name)
    end
    robot_state_msg.num_joints = length(robot_state_msg.joint_name)
    resize!(robot_state_msg)

    let receiver = receiver
        sub = subscribe(lcm, robot_command_channel, (channel, data) -> handle_robot_command_msg(receiver, data))
        set_queue_capacity(sub, 1)
    end

    receiver
end

function initialize!(receiver::LCMControlReceiver)
    receiver.tprev[] = 0
    receiver.τprev .= 0
    nothing
end

function handle_robot_command_msg(receiver::LCMControlReceiver, data::AbstractVector{UInt8})
    setdata!(receiver.decodebuffer, data)
    msg = receiver.atlas_command_msg
    decode!(msg, receiver.decodebuffer)
    receiver.new_command[] = true
end

function publish_robot_state(receiver::LCMControlReceiver, t::Number, state::MechanismState)
    set!(receiver.robot_state_msg, receiver.result, receiver.robot_info, receiver.τprev, t, state)
    encode(receiver.encodebuffer, receiver.robot_state_msg)
    bytes = take!(receiver.encodebuffer)
    publish(receiver.lcm, receiver.robot_state_channel, bytes)
end

struct NoCommandError <: Exception end
Base.showerror(io::IO, e::NoCommandError) = print(io, "Didn't receive a command message.")

function (receiver::LCMControlReceiver)(τ::AbstractVector, t::Number, state::MechanismState)
    # send state info on first tick to get things started
    if t == 0 # TODO: consider creating a separate flag in the controller for doing this
        publish_robot_state(receiver, t, state)
        sleep(0.5); publish_robot_state(receiver, t, state) # send a second time to work around a bug in the controller
    end

    # process command
    LCMCore.handle(receiver.lcm, Second(1))
    receiver.new_command[] || throw(NoCommandError())
    set!(τ, receiver.tprev[] - t, state, receiver.τprev, receiver.atlas_command_msg, receiver.robot_info)
    receiver.new_command[] = false
    receiver.tprev[] = t
    receiver.τprev .= τ

    # send state info
    publish_robot_state(receiver, t, state)
    τ
end

function RigidBodySim.PeriodicController(Δt::Number, receiver::LCMControlReceiver) # TODO: switch order?
    initialize = let receiver = receiver
        (c, t, u, integrator) -> initialize!(receiver)
    end
    PeriodicController(zero(receiver.τprev), Δt, receiver; initialize = initialize)
end


text:
"""
    base_method(M::AbstractOptBuffer)

For a `Wrapper` object, return the wrapped buffer. For an `OptBuffer` object, return the
    object itself.
"""
function base_method(::AbstractOptBuffer) end

argumentvec(M::Wrapper) = argumentvec(base_method(M))
gradientvec(M::Wrapper) = gradientvec(base_method(M))
step_origin(M::Wrapper) = step_origin(base_method(M))
fnval(M::Wrapper) = fnval(base_method(M))
fnval_origin(M::Wrapper) = fnval_origin(base_method(M))

__descent_dir!(M::Wrapper) = __descent_dir!(base_method(M))

"""
    callfn!(fdf, M::AbstractOptBuffer, x, α, d)

Return the value of `fdf(x + α * d)` overwriting the internal buffers of `M` as needed.
    `fdf(x, g)` must return a tuple (f(x), ∇f(x)) and, if `g` is mutable, overwrite
    it with the gradient (see [`optimize!`](@ref)).
"""
@inline callfn!(fdf, M::Wrapper, x, α, d) = callfn!(fdf, base_method(M), x, α, d)

function init!(fdf, M::Wrapper, args...; kw...)
    init!(fdf, base_method(M), args...; kw...)
    return
end

"""
    reset!(M::AbstractOptBuffer, args...; kwargs...)

Reset the solver parameters to the default (or to specific value -- see the documentation
    for the specific types).

Each method has to implement a parameter-free `reset!(M)` method.
"""
function reset!(M::Wrapper, args...; kw...)
    reset!(base_method(M), args...; kw...)
    return
end

"""
    step!(fdf, M::AbstractOptBuffer; kw...)

Make one iteration of optimization routine from the current state of `M` using `fdf` as the
    objective function and modifying the internal buffers as needed.
"""
step!(fn::F, M::Wrapper; kw...) where {F} = step!(fn, base_method(M); kw...)

"""
    stopcond(M::AbstractOptBuffer)

Decide if the optimization should be stopped from the state of `M`.
"""
stopcond(M::Wrapper) = stopcond(base_method(M))
@inline stopcond(M::OptBuffer) = false

"""
    conv_success(M::AbstractOptBuffer)

Decide if the optimization has converged from the state of `M` (useful to distinguish
    between stopping by convergence and stopping by iterations count when the convergence
    has been achieved on the limiting iteration).
"""
conv_success(M::Wrapper) = conv_success(base_method(M))
@inline conv_success(M::OptBuffer) = false

iter_count(M::Wrapper) = iter_count(base_method(M))
call_count(M::Wrapper) = call_count(base_method(M))

iter_count(M::OptBuffer) = -1
call_count(M::OptBuffer) = -1

"""
    convstat(M::Wrapper)

For a converged state, return the statistics in the form of `NamedTuple`
  `(converged = true/false, argument, gradient, iterations, calls)`.
  Negative values of `calls` or `iterations` mean that the number has not been tracked.
"""
function convstat(M::Wrapper)
    converged = conv_success(M)
    argument = argumentvec(M)
    gradient = gradientvec(M)
    iterations = iter_count(M)
    calls = call_count(M)
    @logmsg OptLogLevel """

    ==FINAL STATISTICS==
    # converged: $(converged)
    # final gradient: $(gradient)
    # final argument: $(argument)
    # previous argument: $(step_origin(M))
    # final func value: $(fnval(M))
    # previous func value: $(fnval_origin(M))
    # number of iterations: $(iterations)
    # number of function calls: $(calls)
    """
    return (;
        converged,
        argument,
        gradient,
        iterations,
        calls,
    )
end

"""
    BasicConvergenceStats

Wrapper type to provide basic stop conditions: magnitude of gradient is less than the
    specified value, objective function call count exceeds threshold, iteration count
    exceeds threshold.
"""
mutable struct BasicConvergenceStats{T<:AbstractOptBuffer,F} <: Wrapper
    descent::T
    convcond::F
    converged::Bool
    call_limit::Int
    call_count::Int
    iter_limit::Int
    iter_count::Int

    function BasicConvergenceStats(
        M::T;
        convcond::F,
        call_limit=typemax(Int),
        iter_limit=typemax(Int),
    ) where {T<:AbstractOptBuffer, F<:Base.Callable}
        return new{T, F}(M, convcond, false, call_limit, 0, iter_limit, 0)
    end
end

base_method(M::BasicConvergenceStats) = M.descent

function init!(fdf, M::BasicConvergenceStats, args...; kw...)
    M.call_count = 0
    M.iter_count = 0
    M.converged = false
    init!(fdf, M.descent, args...; kw...)
    return
end

function reset!(M::BasicConvergenceStats, args...; kw...)
    M.call_count = 0
    M.iter_count = 0
    M.converged = false
    reset!(M.descent, args...; kw...)
    return
end

function callfn!(fdf::F, M::BasicConvergenceStats, x, α, d) where {F}
    fg = callfn!(fdf, M.descent, x, α, d)
    M.call_count += 1
    return fg
end

function step!(fdf::F, M::BasicConvergenceStats; kw...) where {F}
    s = step!(fdf, M.descent; kw...)
    M.iter_count += 1
    return s
end

function stopcond(M::BasicConvergenceStats)
    base = base_method(M)

    x, xpre = argumentvec(base), step_origin(base)
    y, ypre = fnval(base), fnval_origin(base)
    g = gradientvec(base)
    M.converged = M.convcond(x, xpre, y, ypre, g)

    if M.converged
        return true
    elseif M.call_count >= M.call_limit
        return true
    elseif M.iter_count >= M.iter_limit
        return true
    else
        return stopcond(base)
    end
end

function conv_success(M::BasicConvergenceStats)
    base = base_method(M)
    if M.converged
        return true
    else
        return conv_success(base)
    end
end

@inline iter_count(M::BasicConvergenceStats) = M.iter_count
@inline call_count(M::BasicConvergenceStats) = M.call_count

"""
    ConstrainStepSize

Wrapper type to limit step sizes attempted in optimization,
given a function `(origin, direction) -> max_step`.
"""
struct ConstrainStepSize{F, T<:AbstractOptBuffer} <: Wrapper
    descent::T
    constraint::F

    function ConstrainStepSize(constr::F, optbuf::T) where {F,T<:AbstractOptBuffer}
        return new{F,T}(optbuf, constr)
    end
end

ConstrainStepSize(M::AbstractOptBuffer) = ConstrainStepSize(infstep, M)

base_method(M::ConstrainStepSize) = M.descent

function init!(fdf, M::ConstrainStepSize, args...; kw...)
    init!(fdf, M.descent, args...; constrain_step=M.constraint, kw...)
    return
end

step!(optfn!, M::ConstrainStepSize) = step!(optfn!, M.descent; constrain_step=M.constraint)

struct OptFunc{F<:Base.Callable,M<:AbstractOptBuffer}<:Function
    fdf::F
    buffer::M
end

function (optfn::OptFunc)(x, α, d)
    y, g = callfn!(optfn.fdf, optfn.buffer, x, α, d)
    @logmsg OptLogLevel "$(join(x .+ α .* d, ' ')) $y $(join(g, ' '))"
    y, g
end


text:
#! /usr/bin/julia

# Rosetta Code, Sorting algorithms/Bead sort

function beadsort{T<:Integer}(a::Array{T,1})
    (lo, hi) = extrema(a)
    if lo < 1
        throw(DomainError())
    end
    hi = convert(Int, hi)
    len = length(a)
    abacus = falses(len, hi)
    for (i, v) in enumerate(a)
           abacus[i,1:v] = true
    end
    for i in 1:hi
        v = sum(abacus[:,i])
        if v < len
            abacus[1:end-v,i] = false
            abacus[end-v+1:end,i] = true
        end
    end
    return T[sum(abacus[i,:]) for i in 1:len]
end

a = Uint8[rand(1:typemax(Uint8)) for i in 1:20]
println("Sort of Unsigned Bytes:")
println("    Before Sort:")
println("    ", a)
a = beadsort(a)
println("\n    After Sort:")
println("    ", a, "\n")

a = [rand(1:2^10) for i in 1:20]
println("Sort of Integers:")
println("    Before Sort:")
println("    ", a)
a = beadsort(a)
println("\n    After Sort:")
println("    ", a)

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using Documenter, Transformers

using Transformers.Basic
using Transformers.Basic: MultiheadAttention
using Transformers.Pretrain
using Transformers.Stacks
using Transformers.GenerativePreTrain
using Transformers.BidirectionalEncoder

makedocs(sitename="Transformers.jl",
         pages = Any[
           "Home" => "index.md",
           "Tutorial" => "tutorial.md",
           "Basic" => "basic.md",
           "Stacks" => "stacks.md",
           "Pretrain" => "pretrain.md",
           "Models" => [
             "GPT" => "gpt.md",
             "BERT" => "bert.md",
           ],
           "Datasets" => "datasets.md"
         ],
         )

deploydocs(
    repo = "github.com/chengchingwen/Transformers.jl.git",
)


text:
struct InteriorMap
    var
    pde
    I
    lower
    upper
end

#to get an equal mapping, you want to associate every equation to a unique dependent variable that it's solving for
# this is tearing
# so there's u and v
# first equation has second derivative of u so we assign that one to u
# second equation, assume it has second derivative of both
# then we assign v to it because u is already assigned somewhere else.
# and use the interior based on the assignment

function InteriorMap(pdes, boundarymap, s::DiscreteSpace{N,M}) where {N,M}
    @assert length(pdes) == M "There must be the same number of equations and unknowns, got $(length(pdes)) equations and $(M) unknowns"
    m = buildmatrix(pdes, s)
    varmap = Dict(build_variable_mapping(m, s.ū, pdes))

    # Determine the interiors for each pde
    vlower = []
    vupper = []

    interior = map(pdes) do pde
        u = varmap[pde]
        boundaries = reduce(vcat, collect(values(boundarymap[operation(u)])))
        n = ndims(u, s)
        lower = zeros(Int, n)
        upper = zeros(Int, n)
        # Determine thec number of points to remove from each end of the domain for each dimension
        for b in boundaries
            #@show b
            clip_interior!!(lower, upper, s, b)
        end
        push!(vlower, pde => lower)
        push!(vupper, pde => upper)
        args = remove(arguments(u), s.time)
        # Don't update this x2i, it is correct.
        pde => s.Igrid[u][[(1+lower[x2i(s, u, x)]:length(s.grid[x])-upper[x2i(s, u, x)]) for x in args]...]
    end
    pdemap = [k.second => k.first for k in varmap]
    return InteriorMap(varmap, Dict(pdemap), Dict(interior), Dict(vlower), Dict(vupper))
end

function buildmatrix(pdes, s::DiscreteSpace{N,M}) where {N,M}
    m = zeros(Int, M, M)
    elegiblevars = [getvars(pde, s) for pde in pdes]
    u2i = Dict([u => k for (k, u) in enumerate(s.ū)])
    #@show elegiblevars, s.ū
    for (i, varmap) in enumerate(elegiblevars)
        for var in keys(varmap)
            m[i, u2i[var]] = varmap[var]
        end
    end
    return m
end

function build_variable_mapping(m, vars, pdes)
    notzero(x) = x > 0 ? 1 : 0
    varpdemap = []
    N = length(pdes)
    rows = sum(m, dims=2)
    cols = sum(m, dims=1)
    i = findfirst(isequal(0), rows)
    j = findfirst(isequal(0), cols)
    @assert i === nothing "Equation $(pdes[i[1]]) is not an equation for any of the dependent variables."
    @assert j === nothing "Variable $(vars[j[2]]) does not appear in any equation, therefore cannot be solved for"
    for k in 1:N
        # Check if any of the pdes only have one valid variable
        m_ones = notzero.(m)
        cols = sum(m_ones, dims=1)
        j = findfirst(isequal(1), cols)
        if j !== nothing
            j = j[2]
            for i in 1:N
                if m[i, j] > 0
                    push!(varpdemap, pdes[i] => vars[j])
                    m[i, :] .= 0
                    m[:, j] .= 0
                    break
                end
            end
            continue
        end
        # Check if any of the variables only have one valid pde
        rows = sum(m_ones, dims=2)
        i = findfirst(isequal(1), rows)
        if i !== nothing
            i = i[1]
            for j in 1:N
                if m[i, j] > 0
                    push!(varpdemap, pdes[i] => vars[j])
                    m[i, :] .= 0
                    m[:, j] .= 0
                    break
                end
            end
            continue
        end
        # Check if any of the variables have more than one valid pde, and pick one
        I = findmax(m)[2]
        i = I[1]
        j = I[2]
        push!(varpdemap, pdes[i] => vars[j])
        m[i, :] .= 0
        m[:, j] .= 0
    end
    @assert length(varpdemap) == N "Could not map all PDEs to variables to solve for, the system is unbalanced."
    return varpdemap
end

@inline function split(children)
    count = [child[1] for child in children]
    # if there are more than one depvars at that order, include all of them
    is = findall(a -> a == maximum(count), count)
    vars = unique(reduce(vcat, [child[2] for child in children[is]]))
    return (sum(count), vars)
end
"""
Creates a ranking of the variables in the term, based on their derivative order.
The heuristic that should work is, if there's a time derivative then use that variable, otherwise use the highest derivative for that variable. If there are two with the highest derivative, pick first from the list that hasn't been chosen for another equation
"""
function get_ranking!(varmap, term, x, s)
    if !istree(term)
        return (0, [])
    end
    S = Symbolics
    SU = SymbolicUtils
    #@show term
    if findfirst(isequal(term), s.ū) !== nothing
        if varmap[term] < 1
            varmap[term] = 1
        end
        return (1, [term])
    else
        op = SU.operation(term)
        children = map(arg -> get_ranking!(varmap, arg, x, s), SU.arguments(term))
        count, vars = split(children)
        if op isa Differential && isequal(op.x, x)
            for var in vars
                if varmap[var] < count + 1
                    varmap[var] = count + 1
                end
            end
            return (1 + count, vars)
        end
        return (count, vars)
    end
end

function getvars(pde, s)
    ct = 0
    ut = []
    # Create ranking for each variable
    varmap = Dict([u => 0 for u in s.ū])
    if s.time !== nothing
        l = get_ranking!(varmap, pde.lhs, s.time, s)
        r = get_ranking!(varmap, pde.rhs, s.time, s)
        for u in s.ū
            if varmap[u] > 0
                varmap[u] += div(typemax(Int) + 1, 2) #Massively weight derivatives in time
            end
        end
    end
    for x in s.x̄
        l = get_ranking!(varmap, pde.lhs, x, s)
        r = get_ranking!(varmap, pde.rhs, x, s)
        split([l, r])
    end
    return varmap
end


text:
"""
    AdamsBashforthTimeStepper(float_type, arch, grid, tracers, χ=0.125;
                              Gⁿ = TendencyFields(arch, grid, tracers),
                              G⁻ = TendencyFields(arch, grid, tracers))

Return an AdamsBashforthTimeStepper object with tendency fields on `arch` and
`grid` with AB2 parameter `χ`. The tendency fields can be specified via optional
kwargs.
"""
struct AdamsBashforthTimeStepper{T, TG} <: AbstractTimeStepper
     χ :: T
    Gⁿ :: TG
    G⁻ :: TG
end

AdamsBashforthTimeStepper(float_type, arch, grid, tracers, χ=0.125;
    Gⁿ = TendencyFields(arch, grid, tracers),
    G⁻ = TendencyFields(arch, grid, tracers)
    ) = AdamsBashforthTimeStepper{float_type, typeof(Gⁿ)}(χ, Gⁿ, G⁻)

#####
##### Time steppping
#####

"""
    time_step!(model::IncompressibleModel{<:AdamsBashforthTimeStepper}, Δt; euler=false)

Step forward `model` one time step `Δt` with a 2nd-order Adams-Bashforth method and
pressure-correction substep. Setting `euler=true` will take a forward Euler time step.
"""
function time_step!(model::IncompressibleModel{<:AdamsBashforthTimeStepper}, Δt; euler=false)
    χ = ifelse(euler, convert(eltype(model.grid), -0.5), model.timestepper.χ)

    # Convert NamedTuples of Fields to NamedTuples of OffsetArrays
    velocities, tracers, pressures, diffusivities, Gⁿ, G⁻ =
        datatuples(model.velocities, model.tracers, model.pressures, model.diffusivities,
                   model.timestepper.Gⁿ, model.timestepper.G⁻)

    ab2_store_previous_source_terms!(G⁻, model.architecture, model.grid, Gⁿ)

    calculate_explicit_substep!(Gⁿ, velocities, tracers, pressures, diffusivities, model)

    ab2_update_source_terms!(Gⁿ, model.architecture, model.grid, χ, G⁻)

    calculate_pressure_correction!(pressures.pNHS, Δt, Gⁿ, velocities, model)

    complete_pressure_correction_step!(velocities, Δt, tracers, pressures, Gⁿ, model)

    model.clock.time += Δt
    model.clock.iteration += 1

    return nothing
end

#####
##### Adams-Bashforth-specific kernels
#####

""" Store previous source terms before updating them. """
function ab2_store_previous_source_terms!(G⁻, arch, grid, Gⁿ)

    # Velocity fields
    @launch(device(arch), config=launch_config(grid, :xyz),
            ab2_store_previous_velocity_source_terms!(G⁻, grid, Gⁿ))

    # Tracer fields
    for i in 4:length(G⁻)
        @inbounds Gc⁻ = G⁻[i]
        @inbounds Gcⁿ = Gⁿ[i]
        @launch(device(arch), config=launch_config(grid, :xyz),
                ab2_store_previous_tracer_source_term!(Gc⁻, grid, Gcⁿ))
    end

    return nothing
end

""" Store previous source terms for `u`, `v`, and `w` before updating them. """
function ab2_store_previous_velocity_source_terms!(G⁻, grid, Gⁿ)
    @loop_xyz i j k grid begin
        @inbounds G⁻.u[i, j, k] = Gⁿ.u[i, j, k]
        @inbounds G⁻.v[i, j, k] = Gⁿ.v[i, j, k]
        @inbounds G⁻.w[i, j, k] = Gⁿ.w[i, j, k]
    end
    return nothing
end

""" Store previous source terms for a tracer before updating them. """
function ab2_store_previous_tracer_source_term!(Gc⁻, grid, Gcⁿ)
    @loop_xyz i j k grid begin
        @inbounds Gc⁻[i, j, k] = Gcⁿ[i, j, k]
    end
    return nothing
end

"""
Evaluate the right-hand-side terms for velocity fields and tracer fields
at time step n+½ using a weighted 2nd-order Adams-Bashforth method.
"""
function ab2_update_source_terms!(Gⁿ, arch, grid, χ, G⁻)
    # Velocity fields
    @launch(device(arch), config=launch_config(grid, :xyz),
            ab2_update_velocity_source_terms!(Gⁿ, grid, χ, G⁻))

    # Tracer fields
    for i in 4:length(Gⁿ)
        @inbounds Gcⁿ = Gⁿ[i]
        @inbounds Gc⁻ = G⁻[i]
        @launch(device(arch), config=launch_config(grid, :xyz),
                ab2_update_tracer_source_term!(Gcⁿ, grid, χ, Gc⁻))
    end

    return nothing
end

"""
Evaluate the right-hand-side terms at time step n+½ using a weighted 2nd-order
Adams-Bashforth method

    `G^{n+½} = (3/2 + χ)G^{n} - (1/2 + χ)G^{n-1}`
"""
function ab2_update_velocity_source_terms!(Gⁿ, grid::AbstractGrid{FT}, χ, G⁻) where FT
    @loop_xyz i j k grid begin
        @inbounds Gⁿ.u[i, j, k] = (FT(1.5) + χ) * Gⁿ.u[i, j, k] - (FT(0.5) + χ) * G⁻.u[i, j, k]
        @inbounds Gⁿ.v[i, j, k] = (FT(1.5) + χ) * Gⁿ.v[i, j, k] - (FT(0.5) + χ) * G⁻.v[i, j, k]
        @inbounds Gⁿ.w[i, j, k] = (FT(1.5) + χ) * Gⁿ.w[i, j, k] - (FT(0.5) + χ) * G⁻.w[i, j, k]
    end

    return nothing
end

"""
Evaluate the right-hand-side terms at time step n+½ using a weighted 2nd-order
Adams-Bashforth method

    `G^{n+½} = (3/2 + χ)G^{n} - (1/2 + χ)G^{n-1}`
"""
function ab2_update_tracer_source_term!(Gcⁿ, grid::AbstractGrid{FT}, χ, Gc⁻) where FT
    @loop_xyz i j k grid begin
        @inbounds Gcⁿ[i, j, k] = (FT(1.5) + χ) * Gcⁿ[i, j, k] - (FT(0.5) + χ) * Gc⁻[i, j, k]
    end
    return nothing
end


text:
using Test

include("game.jl")

# This exercise isn't about mathematical methods to determine the distribution of a random sample,
# therefore we only test that all possible values are returned at least once.

@testset "fog" begin
    @test Set(is_foggy() for _ in 1:1000) == Set([true, false])
end

@testset "distracted dog" begin
    @test Set(is_dog_distracted() for _ in 1:1000) == Set([true, false])
end

@testset "loot purse" begin
    @test Set(loot() for _ in 1:1000) == Set(3:13)
end

@testset "loot crate" begin
    crate = Set(["Cabbage", "Daring Dagger", "Sneaky Shoes"])
    @test Set(loot(crate) for _ in 1:1000) == crate
end


text:
"""
	remove(strains::Array{<:AbstractStrain}, names, [field=:strain]; verbose=false)

Remove strains for which `s.data[field]` is in `names`.
"""
function remove(strains::Array{<:AbstractStrain}, names, field=:strain; verbose=false)
	out = Array{eltype(strains),1}(undef, length(strains))
	idx = 1
	for s in strains
		if !in(s[field], names)
			out[idx] = s
			idx += 1
		end
	end
	verbose && println("Filtered $(length(strains) - idx + 1) strains")
	deleteat!(out, idx:length(strains))
end
"""
	remove!(strains::Array{<:AbstractStrain}, names, [field=:strain]; verbose=false)

Remove strains for which `s.data[field]` is in `names`.
"""
function remove!(strains::Array{<:AbstractStrain}, names, field=:strain; verbose=false)
	indices = Int64[]
	for (i,s) in enumerate(strains)
		if in(s[field], names)
			push!(indices, i)
		end
	end
	verbose && println("Filtered $(length(indices)) strains")
	deleteat!(strains, indices)
end

"""
	gapfilter(s::Strain; threshold=0.1)

Return true if `s.seq` has less than `threshold * length(s.seq)` gaps.
"""
gapfilter(s::Strain; threshold=0.05) = countgaps(s.seq) < threshold * length(s.seq)
"""
	gapfilter(s::ArtificialStrain; threshold=0.)
"""
gapfilter(s::ArtificialStrain; threshold=0.)  = true
# gapfilter(threshold::Float64) = s->gapfilter(s, threshold)

hasdate(x::AbstractStrain) = (!ismissing(get(x.data, "date", missing)) || !ismissing(get(x.data, :date, missing)))

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


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

using Markdown
using InteractiveUtils

# ╔═╡ 3244173c-e227-11eb-39eb-93a74dce1c9e
using PlutoUI, PDFIO, Taro, WordTokenizers, DeepDiffs

# ╔═╡ 0f5786cb-6430-4a14-bc05-bb127f8b73df
md"""
This notebook does not work!
"""

# ╔═╡ c3353ab9-0ae1-4785-b561-631945710a35
#Taro.init()

# ╔═╡ 74f3335a-db62-49ba-8a62-84b6eebed5a2
function getPDFText(src)
    doc = pdDocOpen(src)
    docinfo = pdDocGetInfo(doc)
	result = IOBuffer()
	npage = pdDocGetPageCount(doc)
	for i=1:npage
		page = pdDocGetPage(doc, i)
		pdPageExtractText(result, page)
	end
    pdDocClose(doc)
    return docinfo, result #
end

# ╔═╡ 9b83224c-1529-4ab7-9094-ba3b7b3fe3ba
document1 = "some.pdf"

# ╔═╡ 9318c6dc-dded-403e-8e03-f994c86027c9
document2 = "other.pdf"

# ╔═╡ cf87add1-239a-4917-b721-8482fae365f2
function readpdf(name)
	docinfo, iobuffer = getPDFText(name)
	(docinfo, String(take!(iobuffer)))
end

# ╔═╡ c71c6c25-c00a-42fc-9df7-5a28984bb69f
info1, text1 = readpdf(document1)

# ╔═╡ cb3b3dca-4712-4292-8e54-0f23dc934e38
info2, text2 = readpdf(document2)

# ╔═╡ e4055352-93a6-4b3c-9485-9869800a2ea6
tokenized_text1 = tokenize.(split_sentences(text1));

# ╔═╡ 0870f966-1242-47b0-aba9-e4f6a1f67727
tokenized_text1

# ╔═╡ eb80a7b3-de56-463b-b889-844417aefba1
tokenized_text2 = tokenize.(split_sentences(text2));

# ╔═╡ f992f24e-b6ce-43bd-a4c3-bf83dba86649
tokenized_text2

# ╔═╡ a2a5d048-6440-4d50-b81b-4a98cf08989d
deepdiff(tokenized_text1, tokenized_text2)

# ╔═╡ 71165030-3fe6-4879-89f7-d3f1ce1904b7

# ╔═╡ c102e8e6-2d8f-4582-b85c-30fdd7029be6
with_terminal() do
	print(text1)
end

# ╔═╡ 50f96bd0-21c9-4f7b-8b84-4b1425648395
with_terminal() do
	print(text2)
end

# ╔═╡ 79538987-e2bc-442c-bc75-4e9576c0ef73

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
PDFIO = "4d0d745f-9d9a-592e-8d18-1ad8a0f42b92"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Taro = "61d0e4fa-4e73-5030-88a9-ae4c27b203dd"
WordTokenizers = "796a5d58-b03d-544a-977e-18100b691f6e"

[compat]
DeepDiffs = "~1.2.0"
PDFIO = "~0.1.12"
PlutoUI = "~0.7.9"
Taro = "~0.8.3"
WordTokenizers = "~0.5.6"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[AbstractTrees]]
deps = ["Markdown", "Test"]
git-tree-sha1 = "6621d9645702c1c4e6970cc6a3eae440c768000b"
uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
version = "0.2.1"

[[AdobeGlyphList]]
deps = ["DelimitedFiles", "Pkg", "Test"]
git-tree-sha1 = "dce62e14d1eb385ef241a05094a7bc9511a276dc"
uuid = "715cd884-8f22-55db-b077-449d1f7eef73"
version = "0.1.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[BinDeps]]
deps = ["Compat", "Libdl", "SHA", "URIParser"]
git-tree-sha1 = "12093ca6cdd0ee547c39b1870e0c9c3f154d9ca9"
uuid = "9e28174c-4ba2-5203-b857-d8d62c4213ee"
version = "0.8.10"

[[BinaryProvider]]
deps = ["Libdl", "Logging", "SHA"]
git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058"
uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
version = "0.5.10"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "b0b7e8a0d054fada22b64095b46469627a138943"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "2.2.1"

[[DataAPI]]
git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.7.0"

[[DataDeps]]
deps = ["BinaryProvider", "HTTP", "Libdl", "Reexport", "SHA", "p7zip_jll"]
git-tree-sha1 = "4f0e41ff461d42cfc62ff0de4f1cd44c6e6b3771"
uuid = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
version = "0.7.7"

[[DataStructures]]
deps = ["InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "88d48e133e6d3dd68183309877eac74393daa7eb"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.17.20"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DeepDiffs]]
git-tree-sha1 = "9824894295b62a6a4ab6adf1c7bf337b3a9ca34c"
uuid = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
version = "1.2.0"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[HTML_Entities]]
deps = ["StrTables"]
git-tree-sha1 = "c4144ed3bc5f67f595622ad03c0e39fa6c70ccc7"
uuid = "7693890a-d069-55fe-a829-b4a6d304f0ee"
version = "1.0.1"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "c6a1fff2fd4b1da29d3dccaffb1e1001244d844e"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.12"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.1"

[[JavaCall]]
deps = ["DataStructures", "Dates", "Libdl", "WinReg"]
git-tree-sha1 = "75ec422f88cca591b379a68f85dd24dc5c037291"
uuid = "494afd89-becb-516b-aafa-70d2670c0337"
version = "0.7.7"

[[LabelNumerals]]
deps = ["RomanNumerals", "Test"]
git-tree-sha1 = "03f5757e29f13d76a93416c0b913c19bf8f9baac"
uuid = "2e475f56-4567-5c41-94ff-3f69a29378f9"
version = "0.1.0"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "15003dcb7d8db3c6c857fda14891a539a8f2705a"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.10+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[PDFIO]]
deps = ["AbstractTrees", "AdobeGlyphList", "BinDeps", "Dates", "DelimitedFiles", "LabelNumerals", "Libdl", "LinearAlgebra", "OpenSSL_jll", "Pkg", "Printf", "Rectangle", "RomanNumerals", "Zlib_jll"]
git-tree-sha1 = "9e793da6e602ba4440caf087adcf890485a4706f"
uuid = "4d0d745f-9d9a-592e-8d18-1ad8a0f42b92"
version = "0.1.12"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Primes]]
git-tree-sha1 = "afccf037da52fa596223e5a0e331ff752e0e845c"
uuid = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
version = "0.5.0"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Rectangle]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "2bf927306e130e4718449a8b2c72a6b71afe96e0"
uuid = "9a9db56c-8f71-5460-add5-cb450131785e"
version = "0.1.2"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[RomanNumerals]]
deps = ["Compat", "Primes"]
git-tree-sha1 = "06e621d4a9eda8d09eab8aade41e85cb93cadae9"
uuid = "37834d88-8936-577c-80c9-1066ecf66832"
version = "0.3.1"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StrTables]]
deps = ["Dates"]
git-tree-sha1 = "5998faae8c6308acc25c25896562a1e66a3bb038"
uuid = "9700d1a9-a7c8-5760-9816-a99fda30bb8f"
version = "1.0.1"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "8ed4a3ea724dac32670b062be3ef1c1de6773ae8"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.4.4"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Taro]]
deps = ["Dates", "JavaCall", "Tables"]
git-tree-sha1 = "fad39fbbe8a3eb412c87d9ff85e24392b694aabc"
uuid = "61d0e4fa-4e73-5030-88a9-ae4c27b203dd"
version = "0.8.3"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[URIParser]]
deps = ["Unicode"]
git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d"
uuid = "30578b45-9adc-5946-b283-645ec420af67"
version = "0.4.1"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[WinReg]]
deps = ["Test"]
git-tree-sha1 = "808380e0a0483e134081cc54150be4177959b5f4"
uuid = "1b915085-20d7-51cf-bf83-8f477d6f5128"
version = "0.3.1"

[[WordTokenizers]]
deps = ["DataDeps", "HTML_Entities", "StrTables", "Unicode"]
git-tree-sha1 = "01dd4068c638da2431269f49a5964bf42ff6c9d2"
uuid = "796a5d58-b03d-544a-977e-18100b691f6e"
version = "0.5.6"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─0f5786cb-6430-4a14-bc05-bb127f8b73df
# ╠═3244173c-e227-11eb-39eb-93a74dce1c9e
# ╠═c3353ab9-0ae1-4785-b561-631945710a35
# ╠═74f3335a-db62-49ba-8a62-84b6eebed5a2
# ╠═9b83224c-1529-4ab7-9094-ba3b7b3fe3ba
# ╠═9318c6dc-dded-403e-8e03-f994c86027c9
# ╠═cf87add1-239a-4917-b721-8482fae365f2
# ╠═c71c6c25-c00a-42fc-9df7-5a28984bb69f
# ╠═cb3b3dca-4712-4292-8e54-0f23dc934e38
# ╠═e4055352-93a6-4b3c-9485-9869800a2ea6
# ╠═0870f966-1242-47b0-aba9-e4f6a1f67727
# ╠═eb80a7b3-de56-463b-b889-844417aefba1
# ╠═f992f24e-b6ce-43bd-a4c3-bf83dba86649
# ╠═a2a5d048-6440-4d50-b81b-4a98cf08989d
# ╠═71165030-3fe6-4879-89f7-d3f1ce1904b7
# ╠═c102e8e6-2d8f-4582-b85c-30fdd7029be6
# ╠═50f96bd0-21c9-4f7b-8b84-4b1425648395
# ╠═79538987-e2bc-442c-bc75-4e9576c0ef73
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: StatisticalMice/julia-tutorials

max_issues_repo_name: StatisticalMice/julia-tutorials

max_forks_repo_name: StatisticalMice/julia-tutorials


text:
abstract type AbstractProduction end

"""
    Production(y_P_V, j_P_mai, n_N_P, w_P)

$(FIELDDOCTABLE)
"""
@columns struct Production{MoMo,MoMoD,GMo} <: AbstractProduction
    # Field        | Default | Unit            | Bounds       | Log | Description
    y_P_V::MoMo    | 0.02    | mol*mol^-1      | (0.0,1.0)    | _   | "Product formation linked to growth"
    j_P_mai::MoMoD | 0.001   | mol*mol^-1*d^-1 | (0.0,0.1)    | _   | "Product formation linked to maintenance"
    n_N_P::MoMo    | 0.1     | mol*mol^-1      | (0.0, 1.0)   | _   | "N/C in product (wood)"
    w_P::GMo       | 25.0    | g*mol^-1        | (10.0, 40.0) | _   | "Mol-weight of shoot product (wood)"
end

for fn in fieldnames(Production)
    @eval $fn(p::Production) = p.$fn
end

growth_production!(o, growth) = growth_production!(production_pars(o), o, growth)
growth_production!(p::Production, o, growth) = flux(o)[:P,:gro] = growth * p.y_P_V
growth_production!(p::Nothing, o, growth) = zero(growth)

maintenance_production!(o, u) = maintenance_production!(production_pars(o), o, u)
maintenance_production!(p::Production, o, u) = flux(o)[:P,:mai] = p.j_P_mai * tempcorrection(o) * u[:V]
maintenance_production!(p::Nothing, o, u) = zero(eltype(flux(o)))


text:
module BenchMissingDot

using BenchmarkTools
using Transducers

function manual_missing_dot(x, y)
    acc = zero(eltype(x))
    @inbounds for i = eachindex(x)
        xi = x[i]  # this helps type inference
        yi = y[i]
        # Nested if helps type inference:
        if !(xi isa Missing)
            if !(yi isa Missing)
                acc += xi * yi
            end
        end
    end
    return acc
end

function equiv_missing_dot(x, y)
    acc = zero(eltype(x))
    @inbounds for i = eachindex(x)
        xy = x[i] * y[i]
        if !(xy isa Missing)
            acc += xy
        end
    end
    return acc
end

function naive_missing_dot(x, y)
    acc = zero(eltype(x))
    @inbounds for i = eachindex(x)
        if !(x[i] isa Missing) && !(y[i] isa Missing)
            acc += x[i] * y[i]
        end
    end
    return acc
end

random_missings(n::Int, th=2) = [abs(x) > th ? missing : x for x in randn(n)]

suite = BenchmarkGroup()

n = 10^3
let xs = random_missings(n)
    ys = random_missings(n)
    @assert manual_missing_dot(xs, ys) == foldl(
        +,
        eduction(
            OfType(Tuple{Vararg{Number}}) |> MapSplat(*),
            zip(xs, ys)))
end

suite["xf_nota"] = @benchmarkable(
    foldl(+, ed),
    setup=(ed = eduction(
        MapSplat(*) |> NotA(Missing),
        zip(random_missings.(($n, $n))...))))

suite["xf"] = @benchmarkable(
    foldl(+, ed),
    setup=(ed = eduction(
        OfType(Tuple{Vararg{Number}}) |> MapSplat(*),
        zip(random_missings.(($n, $n))...))))

# This is a bit "cheating" since it's using non-public API.  It is
# just to show the lower-bound of Transducers.jl runtime:
rf_nota = reducingfunction(
    MapSplat(*) |> NotA(Missing),
    +;
    # simd = true,
)
suite["rf_nota"] = @benchmarkable(
    transduce($rf_nota, 0.0, zs),
    setup=(zs = zip(random_missings.(($n, $n))...)))

rf_oftype = reducingfunction(
    OfType(Tuple{Vararg{Number}}) |> MapSplat(*),
    +;
    # simd = true,
)
suite["rf"] = @benchmarkable(
    transduce($rf_oftype, 0.0, zs),
    setup=(zs = zip(random_missings.(($n, $n))...)))

suite["man"] = @benchmarkable(
    manual_missing_dot(xs, ys),
    setup=((xs, ys) = random_missings.(($n, $n))))
suite["equiv"] = @benchmarkable(
    equiv_missing_dot(xs, ys),
    setup=((xs, ys) = random_missings.(($n, $n))))
suite["naive"] = @benchmarkable(
    naive_missing_dot(xs, ys),
    setup=((xs, ys) = random_missings.(($n, $n))))

end  # module
BenchMissingDot.suite


text:
# Diffusion maps
# --------------
# Diffusion maps,
# Coifman, R. & Lafon, S., Applied and Computational Harmonic Analysis, Elsevier, 2006, 21, 5-30

#### DiffMap type
immutable DiffMap{T <: Real} <: SpectralResult
    t::Int
    ɛ::Float64
    K::AbstractMatrix{T}
    proj::Projection{T}

    DiffMap{T}(t::Int, ɛ::T, K::AbstractMatrix{T}, proj::Projection{T}) = new(t, ɛ, K, proj)
end

## properties
outdim(M::DiffMap) = size(M.proj, 1)
projection(M::DiffMap) = M.proj
kernel(M::DiffMap) = M.K

## show & dump
function show(io::IO, M::DiffMap)
    print(io, "Diffusion Maps(outdim = $(outdim(M)), t = $(M.t), ɛ = $(M.ɛ))")
end

function dump(io::IO, M::DiffMap)
    show(io, M)
    println(io, "kernel: ")
    Base.showarray(io, M.K, header=false, repr=false)
    println(io)
    println(io, "projection:")
    Base.showarray(io, projection(M), header=false, repr=false)
end

## interface functions
function transform{T<:Real}(::Type{DiffMap}, X::DenseMatrix{T}; d::Int=2, t::Int=1, ɛ::T=1.0)
    transform!(fit(UnitRangeTransform, X), X)

    sumX = sum(X.^ 2, 1)
    K = exp(( sumX' .+ sumX .- 2*At_mul_B(X,X) ) ./ ɛ)

    p = sum(K, 1)'
    K ./= ((p * p') .^ t)
    p = sqrt(sum(K, 1))'
    K ./= (p * p')

    U, S, V = svd(K, thin=false)
    U ./= U[:,1]
    Y = U[:,2:(d+1)]

    return DiffMap{T}(t, ɛ, K, Y')
end


text:
"""A type for all the tensors and tensor like objects.

The rank means the number of indices. For example AbstractTensor{0} is a
scalar, and AbstractTensor{2} a matrix.
"""
abstract type AbstractTensor{rank} <: Node end

"""A type for the tensors which are terminal nodes, i.e not created by
operations.
"""
abstract type TerminalTensor{rank} <: AbstractTensor{rank} end

"""A union type for everything that can be treated as a scalar."""
Scalar = Union{Number, AbstractTensor{0}}

"""Tensors of which we only know in which vector spaces their indices are."""
struct VariableTensor{rank} <: TerminalTensor{rank}
    name::String
    shape::Tuple{Vararg{AbstractVectorSpace}}
    children::Tuple{}
    freeindices::Tuple{}
    function VariableTensor(name::String, shape::Vararg{AbstractVectorSpace})
        new{length(shape)}(name, shape, (), ())
    end
end

function checktensordimensions(x::AbstractArray, Vs::Vararg{AbstractVectorSpace})
    if size(x) == (1,) && length(Vs) == 0
    elseif ndims(x) != length(Vs)
        throw(DomainError((x, Vs) , string(x, " does not  fit into the tensor space ", Vs)))
    else
        for i in 1:ndims(x)
            if size(x)[i] != dim(Vs[i])
                throw(DomainError(Vs[i], string("Dimension ",i,  " of ", x," does not match with vector space rank")))
            end
        end
    end
end

"""Tensor created from a multidimensional array."""
struct Tensor{T, rank} <: TerminalTensor{rank}
    value::Array{T}
    shape::Tuple{Vararg{AbstractVectorSpace}}
    children::Tuple{}
    freeindices::Tuple{}
    # TODO Check x consists of scalars, if possible
    function Tensor(x::Array{T}, Vs::Vararg{AbstractVectorSpace}) where T
        checktensordimensions(x, Vs...)
        new{T, length(Vs)}(x, Vs, (), ())
    end
end

struct DeltaTensor{rank} <: TerminalTensor{rank}
    shape::Tuple{Vararg{AbstractVectorSpace}}
    children::Tuple{}
    freeindices::Tuple{}
    function DeltaTensor(Vs::Vararg{AbstractVectorSpace})
        if length(Vs) == 1
            throw(DomainError(Vs, "Cannot have delta of dimension 1"))
        end
        for V in Vs[2:end]
            if V != Vs[1] && V != Vs[1]'
                throw(DomainError(Vs,
                "The shape of a delta can contain only a single vector space and its dual."))
            end
        end
        new{length(Vs)}(Vs)
    end
end

struct ZeroTensor{rank} <: TerminalTensor{rank}
    shape::Tuple{Vararg{AbstractVectorSpace}}
    children::Tuple{}
    freeindices::Tuple{}
    function ZeroTensor(Vs::Vararg{AbstractVectorSpace})
        new{length(Vs)}(Vs, (), ())
    end
end

"""A tensor where every entry is of the same value."""
struct ConstantTensor{T, rank} <: TerminalTensor{rank}
    shape::Tuple{Vararg{AbstractVectorSpace}}
    children::Tuple{}
    freeindices::Tuple{}
    value::T
    function ConstantTensor(value::T, Vs::Vararg{AbstractVectorSpace}) where (T <: Scalar)
        if value == 0
            return ZeroTensor(Vs...)
        end
        new{T, length(Vs)}(Vs, (), (), value)
    end
end

"""A convenience function. Allows calling Tensor on any scalar."""
Tensor(x::AbstractTensor{0}) = x
ConstantTensor(x::AbstractTensor{0}) = x
"""Turns non-tensor scalar into a tensor."""
Tensor(x::T) where (T <: Scalar) = ConstantTensor(x)

function printtensor(io, s::String, A::AbstractTensor)
    print(io, typeof(A), ", ", s, "shape: ")
    if A.shape == ()
        print(io, "()")
    else
        print(io, A.shape[1])
        for V in A.shape[2:end]
            print(io, "⊗ ",V)
        end
    end
end

Base.show(io::IO, A::Union{Tensor, ConstantTensor}) = printtensor(io, string(A.value, ", "), A)

Base.show(io::IO, A::VariableTensor) = printtensor(io, string(A.name, ", "), A)

Base.show(io::IO, A::Union{DeltaTensor, ZeroTensor}) = printtensor(io, "", A)

max_stars_repo_name: Sohva/TensorDSL

max_issues_repo_name: Sohva/TensorDSL

max_forks_repo_name: Sohva/TensorDSL


text:
using .Flux

# Network -> Flux

activation(x) = GeneralAct(x)

activation(::typeof(identity)) = Id()
activation(::typeof(relu)) = ReLU()

layer(x) = error("Can't use $x as a layer.")

if isdefined(Flux, :Tracker)
    layer(d::Dense) = Layer(Flux.Tracker.data(d.W), Flux.Tracker.data(d.b), activation(d.σ))
else
    layer(d::Dense) = Layer(d.W, d.b, activation(d.σ))
end

network(c::Chain) = Network([layer.(c.layers)...])

Problem(c::Chain, input::AbstractPolytope, output::AbstractPolytope) =
  Problem(network(c), input, output)

# Flux -> Network

_flux(::ReLU) = relu
_flux(::Id) = identity
_flux(f::GeneralAct) = f.f

_flux(m::Layer) = Dense(m.weights, m.bias, _flux(m.activation))
_flux(m::Network) = Chain(_flux.(m.layers)...)

Flux.Chain(m::Network) = _flux(m)


text:
# sparseDirichletMixture.jl

export SparseDirMix, logSDMweights,
    logMarginal, rand;

struct SparseDirMix
  α::Union{Float64, Vector{Float64}}
  β::Float64
  SparseDirMix(α, β) = all(α .> 0.0) && β > 1.0 ? new(α, β) :
    error("Invalid parameter values.")
end

"""
    logSDMweights(d::SparseDirMix)

  Calculate the log of mixture weights for the SDM prior.
### Example
```julia
    α = exp.(rand(5))
    β = 2.0
    d = SparseDirMix(α, β)
    logSDMweights(d)
```
"""
function logSDMweights(d::SparseDirMix)
  K = length(d.α)
  X = reshape(repeat(d.α, inner=K), (K,K)) + Diagonal(d.β*ones(K))
  lgX = loggamma.(X)
  lpg = reshape(sum(lgX, dims=2), K)
  lpg_denom = logsumexp(lpg)

  lw = lpg .- lpg_denom
  lw
end

"""
    logMarginal(x, d::SparseDirMix)

Calculate the log of the SDM-multinomial compound (marginal) probability mass function.

"""
function logMarginal(prior::SparseDirMix, x::Vector{Int})
    lwSDM = logSDMweights(prior)

    K = length(prior.α)
    A = reshape(repeat(prior.α, inner=K), (K,K)) + Diagonal(prior.β*ones(K))
    AX = A + reshape(repeat(x, inner=K), (K,K))

    lnum = [ lmvbeta(AX[k,:]) for k in 1:K ]
    ldenom = [ lmvbeta(A[k,:]) for k in 1:K ]

    lv = lwSDM .+ lnum .- ldenom

    logsumexp( lv )
end

"""
    rand(d::SparseDirMix[, logscale=FALSE])

  Draw from sparse Dirichlet mixture: p(Θ) ∝ Dir(α)⋅∑Θ^β

"""
function Base.rand(d::SparseDirMix; logout::Bool=false)
  K = length(d.α)
  X = reshape(repeat(d.α, inner=K), (K,K)) + Diagonal(d.β*ones(K))
  lgX = loggamma.(X)
  lpg = reshape(sum(lgX, dims=2), K)
  lpg_denom = logsumexp(lpg)

  lw = lpg .- lpg_denom
  z = StatsBase.sample(Weights( exp.(lw) ))

  rDirichlet(X[z,:], logout=logout)
end


text:
using VariationalImaging.OpDenoise

###################
# Learning function
###################
function tv_learning_function(x,data)
    u = denoise(data[2],x)
    cost = 0.5*norm₂²(u-data[1])
    grad = gradient(x,u,data[1])
    return u,cost,grad
end

function denoise(data,x::Real)
    denoise_params = (
        ρ = 0,
        α = x,
		op = FwdGradientOp(),
        # PDPS
        τ₀ = 5,
        σ₀ = 0.99/5,
        accel = true,
        save_results = false,
        maxiter = 1000,
        verbose_iter = 1001,
        save_iterations = false
    )
    st_opt, iterate_opt = initialise_visualisation(false)
    opt_img = op_denoise_pdps(data; iterate=iterate_opt, params=denoise_params)
    finalise_visualisation(st_opt)
    return opt_img
end

function denoise(data,x::AbstractArray)
	p = PatchOp(x,data[:,:,1]) # Adjust parameter size
	x̄ = zeros(p.size_out)
	inplace!(x̄,p,x)
    denoise_params = (
        ρ = 0,
        α = x̄,
		op = FwdGradientOp(),
        # PDPS
        τ₀ = 5,
        σ₀ = 0.99/5,
        accel = true,
        save_results = false,
        maxiter = 1000,
        verbose_iter = 1001,
        save_iterations = false
    )
    st_opt, iterate_opt = initialise_visualisation(false)

    opt_img = op_denoise_pdps(data; iterate=iterate_opt, params=denoise_params)
    finalise_visualisation(st_opt)
    return opt_img
end

function gradient(α,u::AbstractArray{T,3},ū::AbstractArray{T,3}) where T

	M,N,O = size(u)
	grad = 0
	for i = 1:O
		u1 = @view u[:,:,i]
		u2 = @view ū[:,:,i]
		g = gradient(α,u1,u2)
		grad += g
	end
	return grad
end

function gradient(α::Real,u::AbstractArray{T,2},ū::AbstractArray{T,2}) where T
	u = Float64.(Gray{Float64}.(u))
	ū = Float64.(Gray{Float64}.(ū))
	# Obtain Active and inactive sets
	n = size(u,1)

	# Generate centered gradient matrix
	G = createDivMatrix(n)
	Gu = G*u[:]
	nGu = xi(Gu)
	act = nGu .< 1e-12
	inact = 1 .- act
	Act = spdiagm(0=>act)
	Inact = spdiagm(0=>inact)

	# Vector with grad norm in inactive components and one in the active
	den = Inact*nGu+act
	Den = spdiagm(0=>1 ./den)

	# prod KuKuᵗ/norm³
	prodKuKu = prodesc(Gu ./den.^3,Gu)

	Adj = [spdiagm(0=>ones(n^2)) G';
			Act*G+Inact*α*(prodKuKu-Den)*G Inact+eps()*Act]

	Track=[(u[:]-ū[:]);zeros(2*n^2)]
	mult = Adj\Track
	p = @view mult[1:n^2]
	return -p'*(G'*Inact*Den*Gu)
end

function gradient(α::AbstractArray,u::AbstractArray{T,3},ū::AbstractArray{T,3}) where T

	M,N,O = size(u)
	m,n = size(α)
	p = PatchOp(α,u[:,:,1]) # Adjust parameter size
	ᾱ = zeros(p.size_out)
	inplace!(ᾱ,p,α)
	grad = zeros(size(ᾱ[:]))
	for i = 1:O
		u1 = @view u[:,:,i]
		u2 = @view ū[:,:,i]
		g = gradient(ᾱ,u1,u2)
		grad .+= g
	end
	grad = reshape(grad,M,N)
	grad₊ = zeros(size(α))
	inplace!(grad₊,p',grad)
	return grad₊
end

function gradient(α::AbstractArray,u::AbstractArray{T,2},ū::AbstractArray{T,2}) where T
	u = Float64.(Gray{Float64}.(u))
	ū = Float64.(Gray{Float64}.(ū))
	# Obtain Active and inactive sets
	n = size(u,1)
	M,N = size(α)

	# Generate centered gradient matrix
	G = createDivMatrix(n)
	Gu = G*u[:]
	nGu = xi(Gu)
	act = nGu .< 1e-12
	inact = 1 .- act
	Act = spdiagm(0=>act)
	Inact = spdiagm(0=>inact)

	# Vector with grad norm in inactive components and one in the active
	den = Inact*nGu+act
	Den = spdiagm(0=>1 ./den)

	# prod KuKuᵗ/norm³
	prodKuKu = prodesc(Gu ./den.^3,Gu)
	Adj = [spdiagm(0=>ones(n^2)) G';
			Act*G+Inact*spdiagm(0=>[α[:];α[:]])*(prodKuKu-Den)*G Inact+eps()*Act]

	Track=[(u[:]-ū[:]);zeros(2*n^2)]
	mult = Adj\Track
	p = @view mult[1:n^2]
	return -spdiagm(0=>p[:])*(G'*Inact*Den*Gu)
end

max_stars_repo_name: dvillacis/BPLDenoising

max_issues_repo_name: dvillacis/BPLDenoising

max_forks_repo_name: dvillacis/BPLDenoising


text:
## Order downvalues by specificity, or canonically.
## Probably impossible, in general. We start with
## getting some common cases correct.
## But, we get many simple cases wrong.
## Associating Julia data structures with downvalues would be
## more efficient.

# @sym ClearAll(g)
# @sym g(x_^2) := x
# @sym g(x_Integer^2) := x
# @sym g(z^2) := x
# dvs = @sym DownValues(g)

function has_pattern(mx::Mxpr)
    for i in 1:length(mx)
        res = has_pattern(mx[i])
        res[1] == true && return res
    end
    if symname(mhead(mx)) == :Pattern
        return (true, mx[2])
    end
    return (false,false)
end

function has_pattern(x)
    return (false,false)
end

#  a == f(blah,blah) , b == f(foo,barblah), (or they can be at a deeper level)
function isless_patterns(a::Mxpr, b::Mxpr)
    n = min(length(a),length(b))
    for i in 1:n
        (haspata,blanka) = has_pattern(a)
        (haspatb,blankb) = has_pattern(b)
        res = haspatb && (! haspata)
        if haspatb && (! haspata)  return true end
        length(blanka) > length(blankb) && return true
    end
    return false
end


text:
scaling = true
maxiter = 1450

#
# Plain implementation, no Stopping, simplified Armijo (Wolfe) linesearch
#

using Logging
println("\n bfgs with Armijo line search,  ")

logger = Logging.ConsoleLogger(stderr,Logging.Warn)

Logging.with_logger(logger) do
    # comment the definition of B₀ and in the call to have "regular" BFGS
    # use B₀ to hage L-BFGS
    #B₀ =  InverseLBFGSOperator(Float64, n, mem=mem, scaling=scaling)

    @time  iter, fopt, gopt = bfgs(nlp, maxiter = maxiter, scaling = scaling)# , B₀ = B₀);
    @show iter, (gopt), fopt
end

@show nlp.counters

# to compare plain BFGS as implemented by L-BFGS
mem = 401

using LinearOperators
println("\n L-bfgs with Armijo line search, mem = $mem ")

logger = Logging.ConsoleLogger(stderr,Logging.Warn)
Logging.with_logger(logger) do
    reset!(nlp)
    # comment the definition of B₀ and in the call to have "regular" BFGS
    # use B₀ to hage L-BFGS
    B₀ =  InverseLBFGSOperator(Float64, n, mem=mem, scaling=scaling)

    @time  iter, fopt, gopt = bfgs(nlp, maxiter = maxiter, scaling = scaling , B₀ = B₀);
    @show iter, (gopt), fopt
end

@show nlp.counters

#
# Stopping implementation, simplified Armijo (Wolfe) linesearch
#

using Stopping

println("\n bfgs Stopping,  ")
stp = NLPStopping(nlp, NLPAtX(nlp.meta.x0)  )
stp.meta.optimality_check = unconstrained_check
stp.meta.max_iter = maxiter
stp.meta.rtol = 0
using LSDescent

let stp = stp; #nlp = nlp
    logger = Logging.ConsoleLogger(stderr,Logging.Warn)
    Logging.with_logger(logger) do # execute twice to eliminate compilation time randomness
        reset!(nlp)
        # comment the definition of B₀ and in the call to have "regular" BFGS
        # use B₀ to have L-BFGS
        #B₀ =  InverseLBFGSOperator(Float64, n, mem, scaling=scaling)

        @time  stp, B = bfgs_Stop(nlp, scaling = scaling, stp = stp) #, B₀ = B₀);
        @show stp.meta.nb_of_stop, norm(stp.current_state.gx), stp.current_state.fx
    end
end

@show nlp.counters

#
# Stopping implementation, line search by 1D enhanced interval reduction
#

using OneDmin

println("\n bfgs Stopping line search,  ")

let stp = stp; #nlp = nlp
logger = Logging.ConsoleLogger(stderr,Logging.Warn)
    Logging.with_logger(logger) do # execute twice to eliminate compilation time randomness
        reset!(nlp)
        reinit!(stp)
        # comment the definition of B₀ and in the call to have "regular" BFGS
        # use B₀ to have L-BFGS
        #B₀ =  InverseLBFGSOperator(Float64, n, mem, scaling=scaling)
        lslog = Logging.ConsoleLogger(stderr,Logging.Debug)

        @time  stp, B = bfgs_StopLS(nlp, stp = stp, scaling = scaling, # B₀ = B₀,
                                              # line search options below
                                              pick_in = pick_ins2N,
                                              best = false,
                                              strongWolfe = false
                                              #, LS_logger = lslog
                                              )
        @show stp.meta.nb_of_stop, norm(stp.current_state.gx), stp.current_state.fx
    end
end

@show nlp.counters

;

max_stars_repo_name: tmigot/LSDescent

max_issues_repo_name: tmigot/LSDescent

max_forks_repo_name: tmigot/LSDescent


text:
using MetacommunityDynamics
using DynamicGrids
using DynamicGrids: CellRule, SetCellRule, applyrule,Rule,Ruleset;
using Test

test_list = [
  "./constructors/speciespools.jl",
  "./constructors/mechanisms/colonization.jl",

]

global test_n
global error_bit
test_n = 1
error_bit = false

for this_test in test_list
  try
    include(this_test)
    println("[TEST $(lpad(test_n,2))] \033[1m\033[32mPASS\033[0m $(this_test)")
  catch e
    global error_bit = true
    println("[TEST $(lpad(test_n,2))] \033[1m\033[31mFAIL\033[0m $(this_test)")
    showerror(stdout, e, backtrace())
    println()
    throw("TEST FAILED")
  end
  global test_n += 1
end

if error_bit
  throw("Tests failed")
end


text:
# GridArray

export GridArray
struct GridArray{S,T,D} <: AbstractGridComponent{S,T,D}
    domain::Domain{S,D}
    box::BBox{D}
    cell_centered::SVector{D,Bool}
    data::AbstractArray{T,D}
end

function invariant(ga::GridArray)
    size(bbox) == size(data) || return false
    return true
end

domain(ga::GridArray) = ga.domain
bbox(ga::GridArray) = ga.box
box(ga::GridArray) = ga.box

Base.size(ga::GridArray) = size(bbox(ga))
Base.length(ga::GridArray) = length(bbox(ga))

function Base.similar(ga::GridArray{S,T,D}, element_type::Type{R}=T) where {S,T,D,R}
    return GridArray{S,R,D}(ga.domain, ga.box, ga.cell_centered, similar(ga.data, R))
end

function Base.copy(ga::GridArray{S,T,D}) where {S,T,D}
    return GridArray{S,T,D}(ga.domain, ga.box, ga.cell_centered, copy(ga.data))
end

function Base.fill!(dest_ga::GridArray{S,R,D}, a) where {S,R,D}
    fill!(dest_ga.data, a)
    return dest_ga
end

function Base.map(f, ga::GridArray{S,T,D} where {T}, gas::(GridArray{S,T,D} where {T})...) where {S,D}
    resdata = map(f, ga.data, (ga.data for ga in gas)...)
    R = eltype(resdata)
    return GridArray{S,R,D}(ga.domain, ga.box, ga.cell_centered, resdata)
end

function Base.map!(f, dest_ga::GridArray{S,R,D}, ga::GridArray{S,T,D} where {T}, gas::(GridArray{S,T,D} where {T})...) where {S,R,D}
    map!(f, dest_ga.data, ga.data, (ga.data for ga in gas)...)
    return dest_ga
end

function Base.mapreduce(f, op, ga::GridArray{S,T,D} where {T}, gas::(GridArray{S,T,D} where {T})...; init) where {S,D}
    res = mapreduce(f, op, ga.data, (ga.data for ga in gas)...; init=init)
    return res
end

function Base.reduce(op, ga::GridArray; init)
    res = reduce(op, ga.data; init=init)
    return res
end

function Random.rand!(rng::AbstractRNG, dest_ga::GridArray)
    rand!(rng, dest_ga.data)
    return dest_ga
end
Random.rand(rng::AbstractRNG, ga::GridArray) = rand!(rng, similar(ga))

Random.rand!(dest_ga::GridArray) = rand!(Random.GLOBAL_RNG, dest_ga)
Random.rand(ga::GridArray) = rand(Random.GLOBAL_RNG, ga)

function Base.:(==)(ga1::GridArray{S,T1,D}, ga2::GridArray{S,T2,D}) where {S,T1,T2,D}
    return ga1.domain == ga2.domain && ga1.box == ga2.box && ga1.cell_centered == ga2.cell_centered && ga1.data == ga2.data
end

function find_rtol(::Type{T1}, ::Type{T2}) where {T1,T2}
    T = typeof(zero(T1) + zero(T2))
    T = eltype(T)               # needed if T is an array type
    T <: Union{Integer,Rational} && return 0
    return sqrt(eps(T))
end

function Base.isapprox(ga1::GridArray{S,T1,D}, ga2::GridArray{S,T2,D}; atol::Real=0, rtol::Real=atol > 0 ? 0 : find_rtol(T1, T2),
                       nans::Bool=false) where {S,T1,T2,D}
    return ga1.domain == ga2.domain &&
           ga1.box == ga2.box &&
           ga1.cell_centered == ga2.cell_centered &&
           isapprox(ga1.data, ga2.data; atol=atol, rtol=rtol, nans=nans)
end

export iota!
function iota!(dest_ga::GridArray)
    for i in CartesianIndices(Tuple(size(dest_ga.box)))
        dest_ga.data[i] = SVector(Tuple(i)) + min(dest_ga.box) .- 1
    end
    return dest_ga
end

export iota
iota(ga::GridArray{S,T,D}) where {S,T,D} = iota!(similar(ga, SVector{D,Int}))

export coordinate
function coordinate(ga::GridArray{S,X,D}, pos::SVector{D,I}) where {S,X,I,D}
    imin = min(ga.box) + SVector{D,S}(ga.cell_centered) / 2
    imax = max(ga.box) - SVector{D,S}(ga.cell_centered;) / 2
    return linterp(imin, min(ga.domain), imax, max(ga.domain), pos)::SVector{D,S}
end

export coordinates!
function coordinates!(dest_ga::GridArray)
    return map!(i -> coordinate(dest_ga, i), dest_ga, iota(dest_ga))
end

export coordinates
function coordinates(ga::GridArray{S,T,D}) where {S,T,D}
    return map(i -> coordinate(ga, i), iota(ga))
end

export index
function index(ga::GridArray{S,T,D}, coord::SVector{D,S}) where {S,T,D}
    imin = min(ga.box) + SVector{D,S}(ga.cell_centered) / 2
    imax = max(ga.box) - SVector{D,S}(ga.cell_centered;) / 2
    return linterp(min(ga.domain), imin, max(ga.domain), imax, coord)::SVector{D,S}
end

export zero!
zero!(dest_ga::GridArray) = map!(x -> zero(x), dest_ga)
Base.zero(ga::GridArray) = map(x -> zero(eltype(ga)), ga)

function Base.:+(ga1::GridArray{S,T1,D}, ga2s::GridArray{S,T2,D}...) where {S,T1,T2,D}
    return map(+, ga1, ga2s...)
end

Base.:-(ga::GridArray) = map(-, ga)
function Base.:-(ga1::GridArray{S,T1,D}, ga2::GridArray{S,T2,D}) where {S,T1,T2,D}
    return map(-, ga1, ga2)
end

Base.:*(α::Number, ga::GridArray) = map(x -> α * x, ga)
Base.:*(ga::GridArray, α::Number) = map(x -> x * α, ga)
Base.:\(α::Number, ga::GridArray) = map(x -> α \ x, ga)
Base.:/(ga::GridArray, α::Number) = map(x -> x / α, ga)

Base.:*(α::T, ga::GridArray{S,T} where {S}) where {T} = map(x -> α * x, ga)
Base.:*(ga::GridArray{S,T} where {S}, α::T) where {T} = map(x -> x * α, ga)
Base.:\(α::T, ga::GridArray{S,T} where {S}) where {T} = map(x -> α \ x, ga)
Base.:/(ga::GridArray{S,T} where {S}, α::T) where {T} = map(x -> x / α, ga)

function Base.getindex(ga::GridArray{S,T,D}, i::SVector{D,Int}) where {S,T,D}
    arrmin = map(a -> a[begin], SVector(axes(ga.data)))
    return ga.data[CartesianIndex(Tuple(arrmin + i - min(ga.box)))]::T
end
Base.getindex(ga::GridArray, i::Tuple) = ga[SVector(i)]
Base.getindex(ga::GridArray, i::CartesianIndex) = ga[Tuple(i)]
Base.getindex(ga::GridArray, i::Int...) = ga[i]

function evaluate(ga::GridArray{S,T,D}, x::SVector{D,S}) where {S,T,D}
    @assert all(min(ga.domain) .≤ x .≤ max(ga.domain))
    i = map(s -> round(Int, s), index(ga, x))
    @assert all(ga.box.min .≤ i .≤ ga.box.max)
    # Return the value from the nearest point (don't interpolate)
    ga[i]
end

export load_adios
function load_adios(filename::AbstractString, varname::AbstractString)
    file = adios_open_serial(filename, mode_read)
    varref = adios_get(file, varname)
    @assert varref ≢ nothing
    varpath = replace(varname, r"/[^/]*$" => s"")
    gridGlobalOffset = adios_attribute_data(file, "$varpath/gridGlobalOffset")
    @assert gridGlobalOffset ≢ nothing
    gridSpacing = adios_attribute_data(file, "$varpath/gridSpacing")
    @assert gridSpacing ≢ nothing
    # TODO: ensure that geometry == "cartesian"
    # TODO: ensure that dataOrder == "C"
    position = adios_attribute_data(file, "$varname/position")
    @assert position ≢ nothing
    vardata = fetch(varref)
    close(file)

    @assert vardata isa AbstractArray
    D = ndims(vardata)
    T = eltype(vardata)
    @assert gridGlobalOffset isa AbstractVector
    @assert length(gridGlobalOffset) == D
    S = eltype(gridGlobalOffset)
    @assert gridSpacing isa AbstractVector
    @assert length(gridSpacing) == D
    @assert eltype(gridSpacing) == S
    @assert position isa AbstractVector
    @assert length(position) == D
    @assert eltype(position) <: AbstractFloat

    vardomain_lo = SVector{D,S}(gridGlobalOffset...)
    vardomain_delta = SVector{D,S}(gridSpacing...)
    vardomain_hi = vardomain_lo + (SVector{D,Int}(size(vardata)) .- 1) .* vardomain_delta
    vardomain = Domain{S,D}(vardomain_lo, vardomain_hi)
    varbbox = BBox{D}(fill(1, SVector{D,Int}), SVector{D,Int}(size(vardata)))
    var_cell_centered = SVector{D}(position...) .≠ 0
    return GridArray{S,T,D}(vardomain, varbbox, var_cell_centered, vardata)
end


text:
module TransposedVectors

import Base: transpose, ctranspose, conj, size, length, similar, norm, getindex,
       setindex!, convert, map, broadcast, linearindexing, LinearFast, promote_op,
       LinAlg, SparseArrays, hcat, typed_hcat, vcat, typed_vcat, hvcat_fill, promote_typeof

import Base: *, A_mul_Bt, At_mul_B, At_mul_Bt, A_mul_Bc, Ac_mul_B, Ac_mul_Bc
import Base: /, A_rdiv_Bt, At_rdiv_B, At_rdiv_Bt, A_rdiv_Bc, Ac_rdiv_B, Ac_rdiv_Bc
import Base: \, A_ldiv_Bt, At_ldiv_B, At_ldiv_Bt, A_ldiv_Bc, Ac_ldiv_B, Ac_ldiv_Bc

import Base: @propagate_inbounds, @pure

export TransposedVector, transpose_type

include("TransposedVector.jl")
include("hcat.jl")
include("mul.jl")
include("ldiv.jl")
include("rdiv.jl")
include("diagonal.jl")
include("triangular.jl")
include("sparse.jl")
include("qr.jl")
include("norm.jl")

end # module


text:
using HeadlessChromium
using Base.Test

# open dummy target ('about:'' page)

tg1 = Target("about:")
sleep(5)
resp = send(tg1, "Browser.getVersion")

@test haskey(resp, "result")
@test haskey(resp["result"], "protocolVersion")
@test resp["result"]["protocolVersion"] == "1.2"

close(tg1)

# open file target

src = joinpath(dirname(@__FILE__), "example.html")
tg2 = Target("file://$src")
sleep(5)

plotfile = tempname()
send(tg2, "Page.printToPDF", format="A4") do resp
    open(plotfile, "w") do io
        write(io, base64decode(resp["result"]["data"]))
    end
end
sleep(5)  # give some time for the file to be written

@test isfile(plotfile)
@test stat(plotfile).size > 5000

# navigate to another URL

resp = send(tg2, "Page.navigate", url="https://www.yahoo.com")
sleep(5)
@test haskey(resp, "result")

resp = send(tg2, "DOM.getDocument")
@test haskey(resp, "result")
@test haskey(resp["result"], "root")
@test haskey(resp["result"]["root"], "baseURL")

close(tg2)

stopChromium()