text:
# # Partition reduce on GPU

# By using
# [`ReducePartitionBy`](https://juliafolds.github.io/Transducers.jl/dev/reference/manual/#Transducers.ReducePartitionBy),
# per-partition (group) can be computed on GPU in a streaming (single-pass)
# fashion.

using CUDA
using CUDA: @allowscalar

n = 2^24
if has_cuda_gpu()
    xs = CUDA.randn(n)
else
    xs = randn(n)
end
nothing # hide

# `ReducePartitionBy` expects partition to be contiguous; e.g., sorted by the
# key. We will use `floor` as the key. So, plain `sort!` works in this example.

sort!(xs)
nothing # hide

# In GPU, it is convenient to know the output location before computation. So,
# let us build unique index for each partition using `cumsum!`:

function buildindices(f, xs)
    isedge(x, y) = !isequal(f(x), f(y))
    bounds = similar(xs, Bool)
    @views map!(isedge, bounds[2:end], xs[1:end-1], xs[2:end])
    @allowscalar bounds[1] = true
    partitionindices = similar(xs, Int32)
    return cumsum!(partitionindices, bounds)
end

partitionindices_xs = buildindices(floor, xs)
nothing # hide

# ### Counting the size of each partition

import FoldsCUDA  # register the executor
using FLoops
using Transducers

function countparts(partitionindices; ex = nothing)
    nparts = @allowscalar partitionindices[end]
    ys = similar(partitionindices, nparts)

    ## The intra-partition reducing function that reduces each partition to
    ## a 2-tuple of index and count:
    rf_partition = Map(p -> (p, 1))'(ProductRF(right, +))

    index_and_count =
        partitionindices |>
        ReducePartitionBy(
            identity,  # partition by partitionindices
            rf_partition,
            (-1, 0),
        )

    @floop ex for (p, c) in index_and_count
        @inbounds ys[p] = c
    end

    return ys
end

c_xs = countparts(partitionindices_xs)
#-

# ### Computing the average of each partition

function meanparts(xs, partitionindices; ex = nothing)
    nparts = @allowscalar partitionindices[end]
    ys = similar(xs, float(eltype(xs)), nparts)

    ## The intra-partition reducing function that reduces each partition to
    ## a 3-tuple of index, count and sum:
    rf_partition = Map(((i, p),) -> (p, 1, (@inbounds xs[i])))'(ProductRF(right, +, +))

    index_count_and_sum =
        pairs(partitionindices) |>
        ReducePartitionBy(
            ((_, p),) -> p,  # partition by partitionindices
            rf_partition,
            (-1, 0, zero(eltype(ys))),
        )

    @floop ex for (p, c, s) in index_count_and_sum
        @inbounds ys[p] = s / c
    end

    return ys
end

m_xs = meanparts(xs, partitionindices_xs)
#-


text:
module Gherkin

import Base: ==, hash

export Scenario, ScenarioOutline, Feature, FeatureHeader, Given, When, Then

abstract type ParseResult{T} end

struct OKParseResult{T} <: ParseResult{T}
    value::T
end

struct BadParseResult{T} <: ParseResult{T}
    reason::Symbol
    expected::Symbol
    actual::Symbol
end

issuccessful(::OKParseResult{T}) where {T} = true
issuccessful(::BadParseResult{T}) where {T} = false

abstract type ScenarioStep end

function ==(a::T, b::T) where {T <: ScenarioStep}
    a.text == b.text && a.block_text == b.block_text
end
hash(a::T, h::UInt) where {T <: ScenarioStep} = hash((a.text, a.block_text), h)

struct Given <: ScenarioStep
    text::String
    block_text::String

    Given(text::AbstractString; block_text = "") = new(text, block_text)
end
struct When <: ScenarioStep
    text::String
    block_text::String

    When(text::AbstractString; block_text="") = new(text, block_text)
end
struct Then <: ScenarioStep
    text::String
    block_text::String

    Then(text::AbstractString; block_text="") = new(text, block_text)
end

abstract type AbstractScenario end
struct Scenario <: AbstractScenario
    description::String
    tags::Vector{String}
    steps::Vector{ScenarioStep}
end

struct ScenarioOutline <: AbstractScenario
    description::String
    tags::Vector{String}
    steps::Vector{ScenarioStep}
    placeholders::Vector{String}
    examples::AbstractArray
end

struct FeatureHeader
    description::String
    long_description::Vector{String}
    tags::Vector{String}
end

struct Feature
    header::FeatureHeader
    scenarios::Vector{AbstractScenario}
end

mutable struct ByLineParser
    current::String
    rest::Vector{String}
    isempty::Bool

    function ByLineParser(text::String)
        lines = split(text, "\n")
        current = lines[1]
        rest = lines[2:end]
        new(current, rest, false)
    end
end

function consume!(p::ByLineParser)
    if isempty(p.rest)
        p.current = ""
        p.isempty = true
    else
        p.current = p.rest[1]
        p.rest = p.rest[2:end]
    end
end

Base.isempty(p::ByLineParser) = p.isempty

iscurrentlineempty(p::ByLineParser) = strip(p.current) == ""

function parsetags(byline::ByLineParser)
    tags = []
    while !isempty(byline)
        tag_match = matchall(r"(@[^\s]+)", byline.current)
        if isempty(tag_match)
            break
        end

        consume!(byline)
        append!(tags, tag_match)
    end

    tags
end

function parsefeatureheader(byline::ByLineParser) :: ParseResult{FeatureHeader}
    feature_tags = parsetags(byline)

    description_match = match(r"Feature: (?<description>.+)", byline.current)
    if description_match == nothing
        return BadParseResult{FeatureHeader}(:unexpected_construct, :feature, :scenario)
    end
    consume!(byline)

    long_description_lines = []
    while !isempty(byline)
        if iscurrentlineempty(byline)
            consume!(byline)
            break
        end

        push!(long_description_lines, strip(byline.current))
        consume!(byline)
    end

    feature_header = FeatureHeader(description_match[:description],
                                   long_description_lines,
                                   feature_tags)
    return OKParseResult{FeatureHeader}(feature_header)
end

function parseblocktext(byline::ByLineParser)
    consume!(byline)
    block_text_lines = []
    while !isempty(byline)
        line = byline.current
        consume!(byline)
        if ismatch(r"\"\"\"", line)
            break
        end
        push!(block_text_lines, strip(line))
    end
    return OKParseResult{String}(join(block_text_lines, "\n"))
end

function parsescenariosteps(byline::ByLineParser)
    steps = []
    allowed_step_types = Set([Given, When, Then])

    while !isempty(byline)
        if iscurrentlineempty(byline)
            consume!(byline)
            break
        end
        step_match = match(r"(?<step_type>Given|When|Then|And) (?<step_definition>.+)", byline.current)
        block_text_start_match = match(r"\"\"\"", byline.current)
        if step_match == nothing && block_text_start_match == nothing
            return BadParseResult{Vector{ScenarioStep}}(:invalid_step, :step_definition, :invalid_step_definition)
        end

        if block_text_start_match != nothing
            block_text_result = parseblocktext(byline)
            prev_step_type = typeof(steps[end])
            steps[end] = prev_step_type(steps[end].text; block_text=block_text_result.value)
            continue
        end

        step_type = step_match[:step_type]
        step_definition = step_match[:step_definition]
        if step_type == "Given"
            if !(Given in allowed_step_types)
                return BadParseResult{Vector{ScenarioStep}}(:bad_step_order, :NotGiven, :Given)
            end
            step = Given(step_definition)
        elseif step_type == "When"
            if !(When in allowed_step_types)
                return BadParseResult{Vector{ScenarioStep}}(:bad_step_order, :NotWhen, :When)
            end
            step = When(step_definition)
            delete!(allowed_step_types, Given)
        elseif step_type == "Then"
            step = Then(step_definition)
            delete!(allowed_step_types, Given)
            delete!(allowed_step_types, When)
        elseif step_type == "And"
            if isempty(steps)
                return BadParseResult{Vector{ScenarioStep}}(:leading_and, :specific_step, :and_step)
            end
            last_specific_type = typeof(steps[end])
            step = last_specific_type(step_definition)
        end
        push!(steps, step)
        consume!(byline)
    end
    return OKParseResult{Vector{ScenarioStep}}(steps)
end

function parsescenario(byline::ByLineParser)
    tags = parsetags(byline)

    scenario_outline_match = match(r"Scenario Outline: (?<description>.+)", byline.current)
    if scenario_outline_match != nothing
        description = scenario_outline_match[:description]
        consume!(byline)

        # Parse scenario outline steps
        steps_result = parsescenariosteps(byline)
        if !issuccessful(steps_result)
            return steps_result
        end
        steps = steps_result.value

        # Consume the Example: line
        consume!(byline)

        # Consume the placeholders line
        placeholders = matchall(r"(\w+)", byline.current)
        consume!(byline)

        # Parse the examples
        examples = Array{String,2}(length(placeholders), 0)
        while !iscurrentlineempty(byline)
            example = split(strip(byline.current), "|")
            filter!(x -> !isempty(x), example)
            example = map(strip, example)
            consume!(byline)
            examples = [examples example]
        end
        return OKParseResult{ScenarioOutline}(
            ScenarioOutline(description, tags, steps, placeholders, examples))
   end

    scenario_match = match(r"Scenario: (?<description>.+)", byline.current)
    description = scenario_match[:description]
    consume!(byline)

    steps_result = parsescenariosteps(byline)
    if !issuccessful(steps_result)
        return steps_result
    end
    steps = steps_result.value

    OKParseResult{Scenario}(Scenario(description, tags, steps))
end

function parsefeature(text::String) :: ParseResult{Feature}
    byline = ByLineParser(text)

    feature_header_result = parsefeatureheader(byline)
    if !issuccessful(feature_header_result)
        return BadParseResult{Feature}(feature_header_result.reason,
                                       feature_header_result.expected,
                                       feature_header_result.actual)
    end

    scenarios = []
    while !isempty(byline)
        if iscurrentlineempty(byline)
            consume!(byline)
            continue
        end
        scenario_parse_result = parsescenario(byline)
        if issuccessful(scenario_parse_result)
            push!(scenarios, scenario_parse_result.value)
        end
    end

    OKParseResult{Feature}(
        Feature(feature_header_result.value, scenarios))
end

hastag(feature::Feature, tag::String) = tag in feature.header.tags
hastag(scenario::Scenario, tag::String) = tag in scenario.tags

end

ext: jl


text:
using Documenter
using Test

@testset "Manual doctest" begin
    @info "Doctesting Documenter manual"
    doctest(Documenter)

    # Make sure that doctest() fails if there is a manual page with a failing doctest
    # Will need to run it in a Task though, so that we could easily capture the error.
    @info "Doctesting Documenter manual w/ failing doctest"
    tmpfile = joinpath(@__DIR__, "..", "docs", "src", "lib", "internals", "tmpfile.md")
    write(tmpfile, """
    # Temporary source file w/ failing doctest
    ```jldoctest
    julia> 2 + 2
    42
    ```
    """)
    @test isfile(tmpfile)
    # Note: in Julia 1.3 fetch no longer throws the exception direction, but instead
    # wraps it in a TaskFailedException (https://github.com/JuliaLang/julia/pull/32814).
    if isdefined(Base, :TaskFailedException)
        @test_throws TaskFailedException fetch(schedule(Task(() -> doctest(Documenter))))
    else
        @test_throws TestSetException fetch(schedule(Task(() -> doctest(Documenter))))
    end
    println("^^^ Expected error output.")
    rm(tmpfile)
    @test !isfile(tmpfile)
end


text:
"""
    B,T,B_dual,num_it = seysen(H::Array{Td,2}) where Td

Do greedy  Seysen lattice reduction  on the  matrix `H`, returning  `B`, the
reduced lattice basis;  `T` a unimodular matrix that reduces  `H` (i.e. `B =
H*T`); `B_dual`, dual lattice basis (i.e., `B_dual = pinv(B)`); and num_it the number
of iterations (basis updates). See also [`lll`](@ref).

Follows Seysen algorithm in "Lattice Reduction - A Survey with Applications
in Wireless Communications" by D. Wuebben, et al, IEEE Signal Processing
Magazine, 2011.

Technique originally described in "Simultaneous reduction of a lattice basis
and its reciprocal basis" by M. Seysen, Combinatorica, 1993.
http://link.springer.com/article/10.1007%2FBF01202355

# Examples
```jldoctest
julia> H= [1 2; 3 4];B,T = seysen(H); B
2×2 Matrix{Int64}:
 -1  1
  1  1

julia> H= BigFloat.([1.5 2; 3 4]) .+ 2im; B,_= seysen(H); B
2×2 Matrix{Complex{BigFloat}}:
 0.0+1.0im  0.5+0.0im
 0.0+0.0im  1.0+0.0im

```
"""
function seysen(H::Array{Td,2}) where {Td}

if Td<:BigInt || Td<:BigFloat
    Ti= Td<:Complex ? Complex{BigInt} : BigInt
else
    if Td<:Integer
        Ti= Td<:Complex ? Complex{Td} : Td
    else
        Ti= Td<:Complex ? Complex{Int} : Int
    end
end

n,m = size(H); # m-dimensional lattice in an n-dimensional space

# initialization, outputs
B      = copy(H);    # reduced lattice basis
num_it = 0;          # number of iterations
T      = Matrix{Ti}(I, m, m)  # unimodular matrix
A      = (H'*H)*1.0;          # Gram matrix of H
Adual  = inv(A);     # Inverse gram matrix of H
B_dual = H*Adual;    # Dual basis

# calculate update values Λ[s,t] and the corresponding reduction Δ[s,t] in
# Seysen's measure
Λ = zeros(Ti,m,m);
Δ = zeros(Td,m,m)*0.0;

for s = 1:m
    for t = 1:m
        if s != t
	    x = 0.5*(Adual[t,s]/Adual[s,s]-A[t,s]/A[t,t]);
            Λ[s,t] = round(x);
	    AbsΛ   = abs(Λ[s,t])^2;
	    if AbsΛ != 0
	        zw = real(Λ[s,t])*real(x)+imag(Λ[s,t])*imag(x);
	        Δ[s,t]  = Adual[s,s]*A[t,t]*(2*zw-AbsΛ);
	    end
        end
    end
end # - end calculation of Λ and Δ

# find maximum reduction in Seysen's measure (greedy approach)
zw,max_ind = findmax(abs.(Δ[:]));
s,t        = Tuple(CartesianIndices((m,m))[max_ind])

# init loop
do_reduction = true;

if Δ[s,t] == 0  # if no improvement can be achieved
    do_reduction = false;
end

# Lattice reduction loop of subsequent basis updates
while do_reduction

    num_it = num_it + 1;

    B[:,s]  = B[:,s] + Λ[s,t]*B[:,t];    # perform basis update
    T[:,s]  = T[:,s] + Λ[s,t]*T[:,t];    # updater unimodular transformation
    B_dual[:,t] = B_dual[:,t] - Λ[s,t]'*B_dual[:,s];  # update dual basis

    for ind = 1:m
        # update Gram matrix
        if ind != s
            A[s,ind] = A[s,ind]+Λ[s,t]'*A[t,ind];
            A[ind,s] = A[s,ind]';
        else
            A[s,ind] = norm(B[:,s]).^2;
        end
        # update inverse Gram matrix
        if ind != t
            Adual[t,ind] = Adual[t,ind]-Λ[s,t]*Adual[s,ind];
            Adual[ind,t] = Adual[t,ind]';
        else
            Adual[t,ind] = norm(B_dual[:,t])^2;
        end
    end

    # update all possible update values Λ[s,t]
    # and their corresponding reduction Δ[s,t] in Seysen's measure
    for ind1 = 1:m
        for ind2 = 1:m
            if ((ind1==s) | (ind1==t) | (ind2==s) | (ind2==t)) & (ind1!=ind2)
                x = 0.5*(Adual[ind2,ind1]/Adual[ind1,ind1]-A[ind2,ind1]/
                         A[ind2,ind2]);

                # If you get an "InexactError:..." here, try a float type
                # (or float+int type) with more precision. For some reason
                # Seysen is more sensitive than LLL to bit precision.
                Λ[ind1,ind2] = Ti(round(x));

                AbsΛ = abs(Λ[ind1,ind2])^2;
                if AbsΛ != 0
                    zw = real(Λ[ind1,ind2])*real(x)+imag(Λ[ind1,ind2])*imag(x);
                    Δ[ind1,ind2] = Adual[ind1,ind1]*A[ind2,ind2]*(2*zw-AbsΛ);
                else
                    Δ[ind1,ind2] = 0;
                end
            end
        end
    end

    # find maximum reduction in Seysen's measure (greedy approach)
    zw,max_ind = findmax(abs.(Δ[:]));
    s,t = Tuple(CartesianIndices((m,m))[max_ind])

    # if no reduction is possible, exit loop
    if Δ[s,t] == 0
        do_reduction = false;
    end

end # while do_reduction

return B,T,B_dual,num_it
end

max_stars_repo_name: christianpeel/LLLplus

max_issues_repo_name: christianpeel/LLLplus

max_forks_repo_name: christianpeel/LLLplus


text:
module Sugar
# author      : Stefan Harmeling
# date        : 2014-11-16
# description : convenience function, aka syntactic sugar

using Color, Images, Gadfly, DataFrames

import Gadfly.plot
export imagesc, plot, cmblues, cmoranges

## syntactic sugar for Images.jl
using Color, Images
scale0_1(A::Array) = (A - minimum(A))./(maximum(A)-minimum(A))
scale1_255(A::Array) = 254 * scale0_1(A) + 1
cmblues   = colormap("Blues", 255)
cmoranges = colormap("Oranges", 255)
imagesc(x; cm=cmblues, minsize=100, zoom=iceil(minsize ./ [size(x)...])) = ImageCmap(kron(iround(scale1_255(x)), ones(zoom[1],zoom[2])), cm)

##  syntactic sugar for Gadfly.jl
#   plot(v, w)  plots vector w against vector v
#   plot(w)     plots vector w against [1:length(w)]
#   plot(v, A)  plots columns of A against vector v
#   plot(A)     plots columsn of A against [1:size(A,1)]
plot(v::Vector, w::Vector, elements...; mapping...) =
    plot(x=v, y=w, elements...; mapping...)
plot(w::Vector, elements...; mapping...) =
    plot([1:length(w)], w, elements...; mapping...)
plot(v::Vector, A::Matrix, elements...; mapping...) =
    plot(vcat([DataFrame(x=v, y=A[:,i], color="$i")  for i in 1:size(A,2)]...), x="x", y="y", color="color", Geom.line, elements...; mapping...)
plot(A::Matrix, elements...; mapping...) =
    plot([1:size(A,1)], A, elements...; mapping...)

end # module


text:
module Parameters

export Regimen,
       StaticRegimen,
       DynamicRegimen,
       regimen,

       Mean,
       ATE,

       estimand,
       applyparam,
       fluccovar,
       regimens

using ..Qmodels
import ..Common: Parameter, fluccovar
import StatsBase.predict
import Base.==

"""
Represents a particular single time point treatment regimen.
"""
abstract Regimen{T<:AbstractFloat}

"""
A `StaticRegimen` sets treatment to a single value for all observations. Constructed with
`StaticRegimen(a)` where `a` is a floating point 0 or 1.
"""
immutable StaticRegimen{T<:AbstractFloat} <: Regimen{T}
    a::T
    function StaticRegimen(a)
        a == 0 || a == 1 || error(ArgumentError("a should be 0 or 1"))
        new(a)
    end
end
StaticRegimen{T<:AbstractFloat}(a::T) = StaticRegimen{T}(a)
StaticRegimen(a) = StaticRegimen(float(a))

"""
A `DynamicRegimen` sets treatment for each observation to a particular value. Constructed with
`DynamicRegimen(a)` where `a` a vector of floating point is 0 or 1s.
"""
immutable DynamicRegimen{T<:AbstractFloat} <: Regimen{T}
    a::Vector{T}
    function DynamicRegimen(a)
        all(aa -> aa == 0 || aa == 1, a) || error(ArgumentError("Elements of a should be 0 or 1"))
        new(a)
    end
end
DynamicRegimen{T<:AbstractFloat}(a::Vector{T}) = DynamicRegimen{T}(a)
DynamicRegimen(a) = DynamicRegimen(float(a))

==(d1::DynamicRegimen, d2::DynamicRegimen) = d1.a==d2.a

"""
Converts it's argument `a` to a `StaticRegime` if `a` is a scalar, a `DynamicRegime` if `a` is a vector,
or returns `a` if `a` is already a `Regime`.
"""
:regimen

regimen(a::Regimen) = a
regimen(a::Real) = StaticRegimen(float(a))
regimen{T<:Real}(a::Vector{T}) = DynamicRegimen(float(a))

"Returns predicted values from a `Qmodel` for a `StaticRegimen`"
predict{T<:AbstractFloat}(q::Qmodel{T}, d::StaticRegimen{T}) = predict(q, fill(d.a, nobs(q)))
"Returns predicted values from a `Qmodel` for a `DynamicRegimen`"
predict{T<:AbstractFloat}(q::Qmodel{T}, d::DynamicRegimen{T}) = predict(q, d.a, nobs(q))

"""
The `Mean` parameter is \(E_0E_0(Y\mid A=d, W)\).

Under causal assumptions, this can be interpreted as
the mean of the counterfactual outcome under regimen `d`.
"""
type Mean{T<:AbstractFloat} <: Parameter{T}
    d::Regimen{T}
end

Mean(d) = Mean(regimen(d))

"""
The `ATE` parameter is \(E_0[E_0(Y\mid A=d1, W) - E_0(Y\mid A=d0, W)]\).

Under causal assumptions, this can be interpreted as
the difference in mean counterfactual outcome under regimens `d1` and `d0`.  When `d1` is the static regimen setting treatment to 1
and `d0` is the static regimen setting treatment to 0, this is called the averate treatment effect (ATE).
"""
immutable ATE{T<:AbstractFloat} <: Parameter{T}
    d1::Regimen{T}
    d0::Regimen{T}
end

ATE(d1, d0) = ATE(regimen(d1), regimen(d0))
ATE() = ATE(1.0, 0.0)

function fluccovar{T<:AbstractFloat}(param::Mean{T}, a::Vector{T})
    convert(Vector{T}, a .== param.d.a)
end

function fluccovar{T<:AbstractFloat}(param::ATE{T}, a::Vector{T})
    convert(Vector{T}, (a .== param.d1.a) .- (a .== param.d0.a))
end

"""
Computes the parameter estimate and influence curve for a particular `Parmeter` given an estimate of Q and g.

** Arguments **

* `param` - parameter
* `q` - an estimated `Qmodel`
* `g` - vector representing the estimated propensity score for each observation
* `A` - observed treatment vector
* `Y` - observed outcome vector

"""
:applyparam

function applyparam{T<:AbstractFloat}(param::Mean{T}, q::Qmodel{T}, A, Y)
    QnAd = predict(q, param.d)
    psi = mean(QnAd)
    h = weightedcovar(lastfluc(q), A)
    ic = h .* (Y .- QnAd) .+ QnAd #.- psi
    ic .-= mean(ic)
    (psi, ic)
end

function applyparam{T<:AbstractFloat}(param::ATE{T}, q::Qmodel{T}, A, Y)
    QnAd1 = predict(q, param.d1)
    QnAd0 = predict(q, param.d0)
    QnAA = predict(q, A)
    Qndiff = QnAd1 .- QnAd0
    psi = mean(Qndiff)
    h = weightedcovar(lastfluc(q), A)
    ic = h .* (Y .- QnAA) .+ Qndiff # .- psi
    ic .-= mean(ic)
    (psi, ic)
end

estimand(param::Parameter) = "ψ"
estimand(param::Mean) = "E(Y(d))"
estimand(param::ATE) = "ATE"

regimens(param::Mean) = [param.d]
regimens(param::ATE) = [param.d1, param.d0]

end


text:
module pgrad

using Reinforce
using OpenAIGym
using Transformations
using StochasticOptimization
using PenaltyFunctions

using MLPlots; gr(size=(1400,1400), leg=false)

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

# initialize a policy, do the learning, then return the policy
function doit(env = GymEnv("BipedalWalker-v2"))
    s = state(env)
    ns = length(s)
    A = actions(env,s)
    nA = length(A)
    @show s ns

    policy = OnlineActorCritic(s, nA,
        # ϕ = nnet(2ns+nA, 2nA, [100], :softplus, lookback=10000),
        ϕ = resnet(2ns+nA, 2nA, 1, nh=[100], inner_activation=:softplus, lookback=10000),
        penalty = L2Penalty(1e-5),
        γ = 0.995,
        λ = 0.95,
        # αʳ = 0.0001,
        # αᵛ = 0.01,
        # αᵘ = 0.01,
        gaᵛ = OnlineGradAvg(50, lr=0.1, pu=RMSProp()),
        gaᵘ = OnlineGradAvg(50, lr=0.01, pu=RMSProp()),
        # gaʷ = OnlineGradAvg(50, lr=0.5, pu=Adamax())
    )

    # --------------------------------
    # set up the custom visualizations

    # chainplots... put one for each of ϕ/C side-by-side
    ϕ = policy.actor.ϕ
    D = policy.actor.D

    cp_ϕ = ChainPlot(ϕ)
    tpplt = plot(layout=grid(4,2))
    tp_δ = TracePlot(1, sp=tpplt[1], title="delta")
    tp_r̄ = TracePlot(1, sp=tpplt[2], title="r")
    tp_eᵛ = TracePlot(ns, sp=tpplt[3], title="ev")
    tp_v = TracePlot(ns, sp=tpplt[4], title="v")
    tp_eμ = TracePlot(nA, sp=tpplt[5], title="e_mu")
    tp_μ = TracePlot(nA, sp=tpplt[6], title="mu")
    tp_eσ = TracePlot(nA, sp=tpplt[7], title="e_sigma")
    tp_σ = TracePlot(nA, sp=tpplt[8], title="sigma")
    plt = plot(cp_ϕ.plt, tpplt, layout=grid(2,1,heights=[0.3,0.7]))

    # this will be called on every timestep of every episode
    function eachiteration(ep,i)
        @show i, ep.total_reward
        # i<5000 && return

        policy.gaᵛ.lr *= 0.999
        policy.gaᵘ.lr *= 0.999

        # i>=2000 && (policy.αᵘ = 0.0005)
        update!(cp_ϕ)
        for (tp,var) in [(tp_δ, policy.δ), (tp_r̄, mean(policy.r̄)),
                         (tp_eᵛ, policy.eᵛ), (tp_v, policy.v),
                         (tp_eμ, policy.eᵘ[1:nA]), (tp_μ, output_value(ϕ)[1:nA]),
                         (tp_eσ, policy.eᵘ[nA+1:end]), (tp_σ, output_value(ϕ)[nA+1:end]),
                        ]
            push!(tp, i, var)
        end
        gui()
    end

    function renderfunc(ep,i)
        policy.actor.testing = true
        if mod1(ep.niter, 1) == 1
            OpenAIGym.render(env, ep.niter, nothing)
        end
        policy.actor.testing = false
    end

    learn!(policy, Episodes(
        env,
        freq = 1,
        episode_strats = [MaxIter(1000)],
        epoch_strats = [MaxIter(10000), IterFunction(renderfunc, every=5)],
        iter_strats = [IterFunction(eachiteration, every=1000)]
        # append_action = true
    ))

    env, policy
end

end #module


text:
# Test for gaussjacobi().
    tol = 1e-14
    # Test a small n:
    n = 42; a = -.1; b = .3
    x = gaussjacobi(n, -.1, .3)
    w = x[2]; x = x[1]
    @test all(length(x) == n) && all(length(w) == n)
    @test abs(x[37] - 0.912883347814032) < tol
    @test abs(w[37] - 0.046661910947553) < tol

    # Test a larger n (using ASY)
    a = -.7; b = 1.3;
    n = 251;
    x = gaussjacobi(n,-.7,1.3);
    w = x[2]; x = x[1]
    @test all(length(x) == n) && all(length(w) == n)
    @test abs(x[37] + 0.893103435898983) < tol
    @test abs(w[37] - 1.962534523788093e-04) < tol

    # Test n = 1:
    a = 1.0; b = 2.0;
    x = gaussjacobi( 1, a, b );
    @test abs( x[1][1] - (b-a)/(a+b+2) ) < tol
    @test abs( x[2][1] - 2^(a+b+1)*beta(a+1, b+1) ) < tol

    x = gaussjacobi( 1013, .9, -.1 )
    @test abs( x[1][2]  + 0.999986012231899 ) < tol
    @test abs( x[1][13]  + 0.999225722939832 ) < tol
    @test abs( x[2][2] - 9.314674169892358e-05 ) < tol
    @test abs( x[2][13] - 4.654651764553262e-04 ) < tol

    x = gaussjacobi( 10013, .9, -.1 )
    @test abs( x[1][2]  +0.999999856605293 ) < tol
    @test abs( x[1][13]  + 0.999992061552711 ) < tol
    @test abs( x[2][2] - 1.509654630405615e-06 ) < tol
    @test abs( x[2][13] - 7.548275262993863e-06 ) < tol


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

using Markdown
using InteractiveUtils

# ╔═╡ bbdfb85a-309a-11ec-1e68-59d660aab7b3
md"""
# How to Install Genie

Genie development moves fast. Until reaching v1, the recommended installation is to get the latest and greatest Genie, by running off the `master` branch:

```julia
julia> ] # activate package mode
pkg(1.6)> add Genie#master
```

or

```julia
pkg(1.6)> add https://github.com/GenieFramework/Genie.jl
```

Download stable tagged version of Genie

```julia
julia> using Pkg
julia> Pkg.add("Genie")
```

or

```julia
julia> ] # enter package mode
pkg(1.6)> add Genie
```

"""

# ╔═╡ fce3a30a-813f-4875-ba3f-c98008019579
md"""
Genie, just like Julia, uses semantic versioning in the form vX.Y.Z to designate:

- X : major version, introducing breaking changes
- Y : minor version, brings new features, no breaking changes
- Z : patch version, fixes bugs, no new features or breaking changes

---
**HEADS UP**

Pre version 1, changes in Genie's minor version indicate breaking changes. So a new version 0.15 will introduce breaking changes from 0.14. Patch versions indicate non-breaking changes such as new features and patch releases.

---
"""

# ╔═╡ Cell order:
# ╟─bbdfb85a-309a-11ec-1e68-59d660aab7b3
# ╟─fce3a30a-813f-4875-ba3f-c98008019579


text:
# Called in the OrdinaryDiffEQ.__init; All `OrdinaryDiffEqAlgorithm`s have one
function OrdinaryDiffEq.initialize!(integ, cache::GaussianODEFilterCache)
    @assert integ.opts.dense == integ.alg.smooth "`dense` and `smooth` should have the same value! "
    @assert integ.saveiter == 1
    OrdinaryDiffEq.copyat_or_push!(integ.sol.x, integ.saveiter, cache.x)
    OrdinaryDiffEq.copyat_or_push!(integ.sol.pu, integ.saveiter, cache.SolProj*cache.x)
end

"""Perform a step

Not necessarily successful! For that, see `step!(integ)`.

Basically consists of the following steps
- Coordinate change / Predonditioning
- Prediction step
- Measurement: Evaluate f and Jf; Build z, S, H
- Calibration; Adjust prediction / measurement covs if the diffusion model "dynamic"
- Update step
- Error estimation
- Undo the coordinate change / Predonditioning
"""
function OrdinaryDiffEq.perform_step!(integ, cache::GaussianODEFilterCache, repeat_step=false)
    @unpack t, dt = integ
    @unpack d, Proj, SolProj, Precond = integ.cache
    @unpack x, x_pred, u_pred, x_filt, u_filt, err_tmp = integ.cache
    @unpack A, Q = integ.cache

    tnew = t + dt

    # Coordinate change / preconditioning
    P = Precond(dt)
    PI = inv(P)
    x = P * x

    if isdynamic(cache.diffusionmodel)  # Calibrate, then predict cov

        # Predict
        predict_mean!(x_pred, x, A, Q)
        mul!(u_pred, SolProj, PI*x_pred.μ)

        # Measure
        measure!(integ, x_pred, tnew)

        # Estimate diffusion
        integ.cache.diffusion = estimate_diffusion(cache.diffusionmodel, integ)
        # Adjust prediction and measurement
        predict_cov!(x_pred, x, A, apply_diffusion(Q, integ.cache.diffusion))
        copy!(integ.cache.measurement.Σ, X_A_Xt(x_pred.Σ, integ.cache.H))

    else  # Vanilla filtering order: Predict, measure, calibrate

        predict!(x_pred, x, A, Q)
        mul!(u_pred, SolProj, PI*x_pred.μ)
        measure!(integ, x_pred, tnew)
        integ.cache.diffusion = estimate_diffusion(cache.diffusionmodel, integ)

    end

    # Likelihood
    cache.log_likelihood = logpdf(cache.measurement, zeros(d))

    # Update
    x_filt = update!(integ, x_pred)
    mul!(u_filt, SolProj, PI*x_filt.μ)
    integ.u .= u_filt

    # Undo the coordinate change / preconditioning
    copy!(integ.cache.x, PI * x)
    copy!(integ.cache.x_pred, PI * x_pred)
    copy!(integ.cache.x_filt, PI * x_filt)

    # Estimate error for adaptive steps
    if integ.opts.adaptive
        err_est_unscaled = estimate_errors(integ, integ.cache)
        DiffEqBase.calculate_residuals!(
            err_tmp, dt * err_est_unscaled, integ.u, u_filt,
            integ.opts.abstol, integ.opts.reltol, integ.opts.internalnorm, t)
        integ.EEst = integ.opts.internalnorm(err_tmp, t) # scalar

    end
    # stuff that would normally be in apply_step!
    if !integ.opts.adaptive || integ.EEst < one(integ.EEst)
        copy!(integ.cache.x, integ.cache.x_filt)
        integ.sol.log_likelihood += integ.cache.log_likelihood
    end
end

function h!(integ, x_pred, t)
    @unpack f, p, dt = integ
    @unpack u_pred, du, Proj, Precond, measurement = integ.cache
    PI = inv(Precond(dt))
    z = measurement.μ
    E0, E1 = Proj(0), Proj(1)

    u_pred .= E0*PI*x_pred.μ
    IIP = isinplace(integ.f)
    if IIP
        f(du, u_pred, p, t)
    else
        du .= f(u_pred, p, t)
    end
    integ.destats.nf += 1

    z .= E1*PI*x_pred.μ .- du

    return z
end

function H!(integ, x_pred, t)
    @unpack f, p, dt, alg = integ
    @unpack ddu, Proj, Precond, H, u_pred = integ.cache
    E0, E1 = Proj(0), Proj(1)
    PI = inv(Precond(dt))

    if alg isa EK1 || alg isa IEKS
        if alg isa IEKS && !isnothing(alg.linearize_at)
            linearize_at = alg.linearize_at(t).μ
        else
            linearize_at = u_pred
        end

        if isinplace(integ.f)
            f.jac(ddu, linearize_at, p, t)
        else
            ddu .= f.jac(linearize_at, p, t)
            # WIP: Handle Jacobians as OrdinaryDiffEq.jl does
            # J = OrdinaryDiffEq.jacobian((u)-> f(u, p, t), u_pred, integ)
            # @assert J ≈ ddu
        end
        integ.destats.njacs += 1
        mul!(H, (E1 .- ddu * E0), PI)
    else
        mul!(H, E1, PI)
    end

    return H
end

function measure!(integ, x_pred, t)
    @unpack R = integ.cache
    @unpack u_pred, measurement, H = integ.cache

    z, S = measurement.μ, measurement.Σ
    z .= h!(integ, x_pred, t)
    H .= H!(integ, x_pred, t)
    # R .= Diagonal(eps.(z))
    @assert iszero(R)
    copy!(S, X_A_Xt(x_pred.Σ, H))

    return nothing
end

function update!(integ, prediction)
    @unpack measurement, H, R, x_filt = integ.cache
    update!(x_filt, prediction, measurement, H, R)
    assert_nonnegative_diagonal(x_filt.Σ)
    return x_filt
end

function estimate_errors(integ, cache::GaussianODEFilterCache)
    @unpack diffusion, Q, H = integ.cache

    if diffusion isa Real && isinf(diffusion)
        return Inf
    end

    error_estimate = sqrt.(diag(X_A_Xt(apply_diffusion(Q, diffusion), H)))

    return error_estimate
end


text:
import ..UncertainValues:
    AbstractScalarPopulation,
    ConstrainedUncertainScalarPopulation

import StatsBase
import StatsBase: AbstractWeights

"""
    constrain(pop::UncertainScalarPopulation, constraint::SamplingConstraint, n::Int = 30000)

Constrain an `UncertainScalarPopulation` by the given sampling `constraint`. If the
sampling constraint requires resampling to compute, resample `n` times.
"""
function constrain(pop::AbstractScalarPopulation{T, PW}, constraint::SamplingConstraint, n::Int = 30000) where {T, PW}
    # Get the population members and weights that fullfill the constraint and return them as a
    # ConstrainedUncertainScalarPopulation
    truncate(pop, constraint, n)
end

export constrain


text:
heaviside(x::AbstractFloat) = ifelse(x <= 0, zero(x), ifelse(x > 0, one(x), oftype(x,0)))
heaviside(x::Int64) = ifelse(x <= 0, zero(x), ifelse(x > 0, one(x), oftype(x,0)))
heaviside(x::Float64) = ifelse(x <= 0, zero(x), ifelse(x > 0, one(x), oftype(x,0)))

heaviside0(x::AbstractFloat) = ifelse(x < 0, zero(x), ifelse(x >= 0, one(x), oftype(x, 0)))
heaviside0(x::Int64) = ifelse(x < 0, zero(x), ifelse(x >= 0, one(x), oftype(x, 0)))
heaviside0(x::Float64) = ifelse(x < 0, zero(x), ifelse(x >= 0, one(x), oftype(x, 0)))

function heaviside0(v::AbstractArray{Float64})
    m, n = size(v)
    h = zeros(m, n)

    for i in 1:length(v)
      h[i] = heaviside0(v[i])
    end
    return(h)
end

function heaviside0(v::AbstractArray{Int64})
    m, n = size(v)
    h = zeros(m, n)

    for i in 1:length(v)
      h[i] = heaviside0(v[i])
    end
    return(h)
end

function heaviside0(v::Vector{Float64})
    h = zeros(length(v), 1)

    for i in 1:length(v)
      h[i] = heaviside0(v[i])
    end
    return(h)
end

function heaviside0(v::Vector{Int64})
    h = zeros(length(v), 1)

    for i in 1:length(v)
      h[i] = heaviside0(v[i])
    end
    return(h)
end

function heaviside(v::Vector{Float64})
    h = zeros(length(v), 1)

    for i in 1:length(v)
      h[i] = heaviside(v[i])
    end
    return(h)
end

function heaviside(v::AbstractArray{Float64})
    m, n = size(v)
    h = zeros(m, n)

    for i in 1:length(v)
      h[i] = heaviside(v[i])
    end
    return(h)
end

function heaviside(v::AbstractArray{Int64})
    m, n = size(v)
    h = zeros(m, n)

    for i in 1:length(v)
      h[i] = heaviside(v[i])
    end
    return(h)
end

function heaviside(v::Matrix{Int64})
    m, n = size(v)
    h = zeros(m, n)

    for i in 1:length(v)
      h[i] = heaviside(v[i])
    end
    return(h)
end

function heaviside(v::Matrix{Float64})
    m, n = size(v)
    h = zeros(m, n)

    for i in 1:length(v)
      h[i] = heaviside(v[i])
    end
    return(h)
end


text:
using OrdinaryDiffEq, Zygote
using DiffEqSensitivity, Test, ForwardDiff

abstol = 1e-12
reltol = 1e-12
savingtimes = 0.5

# see https://diffeq.sciml.ai/stable/features/callback_functions/#VectorContinuousCallback-Example
function test_vector_continuous_callback(cb,g)

  function f(du, u, p, t)
    du[1] = u[2]
    du[2] = -p[1]
    du[3] = u[4]
    du[4] = 0.0
  end

  u0 = [50.0, 0.0, 0.0, 2.0]
  tspan = (0.0, 10.0)
  p = [9.8, 0.9]
  prob = ODEProblem(f,u0,tspan,p)
  sol = solve(prob, Tsit5(), callback=cb, abstol=abstol, reltol=reltol, saveat=savingtimes)

  du01, dp1 = @time Zygote.gradient(
      (u0, p) -> g(solve(prob, Tsit5(), u0=u0, p=p, callback=cb, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=BacksolveAdjoint())),
      u0, p)

  dstuff = @time ForwardDiff.gradient(
      (θ) -> g(solve(prob, Tsit5(), u0=θ[1:4], p=θ[5:6], callback=cb, abstol=abstol, reltol=reltol, saveat=savingtimes)),
      [u0; p])

  @test du01 ≈ dstuff[1:4]
  @test dp1 ≈ dstuff[5:6]
end

@testset "VectorContinuous callbacks" begin
  @testset "MSE loss function bouncing-ball like" begin
    g(u) = sum((1.0.-u).^2)./2
    function condition(out, u, t, integrator) # Event when event_f(u,t) == 0
      out[1] = u[1]
      out[2] = (u[3]-10.0)u[3]
    end
    @testset "callback with linear affect" begin
      function affect!(integrator, idx)
        if idx == 1
            integrator.u[2] = -integrator.p[2]*integrator.u[2]
        elseif idx == 2
            integrator.u[4] = -integrator.p[2]*integrator.u[4]
        end
      end
      cb = VectorContinuousCallback(condition, affect!, 2)
      test_vector_continuous_callback(cb, g)
    end
  end
end


text:
using Test, DataFrames, Suppressor, HMCResearchRandomWalks.Arlett

@testset "Setters for simulation controls" begin
    set_NUMBER_OF_WALKS()
    @test Arlett.NUMBER_OF_WALKS == 1000
    set_NUMBER_OF_WALKS(10)
    @test Arlett.NUMBER_OF_WALKS == 10

    set_MAX_STEPS_PER_WALK()
    @test Arlett.MAX_STEPS_PER_WALK == 1485604
    set_MAX_STEPS_PER_WALK(2000000)
    @test Arlett.MAX_STEPS_PER_WALK == 2000000

    set_FLOW_BIAS()
    @test Arlett.FLOW_BIAS
    set_FLOW_BIAS(false)
    @test !Arlett.FLOW_BIAS
    set_FLOW_BIAS(true)
    @test Arlett.FLOW_BIAS

    set_CATALASE_ON_WALLS()
    @test Arlett.CATALASE_ON_WALLS
    set_CATALASE_ON_WALLS(false)
    @test !Arlett.CATALASE_ON_WALLS
    set_CATALASE_ON_WALLS(true)
    @test Arlett.CATALASE_ON_WALLS
end

@testset "fast custom arlett_sim" begin
    set_NUMBER_OF_WALKS(10)
    set_MAX_STEPS_PER_WALK(1000000)
    set_FLOW_BIAS(true)
    set_CATALASE_ON_WALLS(true)
    output = @capture_out begin arlett_sim(1234) end
    correct_output = string(
    "############################\n",
    "      Arlett Simulation     \n",
    "############################\n",
    "_________Parameters_________\n",
    "# of trials:\t\t10\n",
    "# of steps (max):\t1000000\n",
    "Step size, water:\t0.1\n",
    "Step size, enz:\t\t0.005\n",
    "Step size, ppd:\t\t0.001\n",
    "Flow bias on:\t\ttrue\n",
    "Catalase on walls:\ttrue\n",
    "Random seed:\t\t1234\n",
    "_________Results____________\n",
    "Side wall collisions:\t0\n",
    "# caught by catalase:\t4\n",
    "# escaped:\t\t2\n",
    "*Sensors from left to right*\n",
    "# in first:\t\t0\n",
    "# in second:\t\t0\n",
    "# in third (spawn):\t2\n",
    "# in fourth:\t\t1\n",
    "# in fifth:\t\t1\n",
    "# unresolved:\t\t0\n",
    "Average steps taken:\t181302\n"
    )
    @test output === correct_output
end

@testset "getdata_arlett_noflow()" begin
    set_NUMBER_OF_WALKS(10)
    set_MAX_STEPS_PER_WALK(1000000)
    output = @capture_out begin
        global gd_arlett_noflow_test_df = getdata_arlett_noflow(1234)
    end
    correct_output = string(
    "############################\n",
    "      Arlett Simulation     \n",
    "############################\n",
    "_________Parameters_________\n",
    "# of trials:\t\t10\n",
    "# of steps (max):\t1000000\n",
    "Step size, water:\t0.1\n",
    "Step size, enz:\t\t0.005\n",
    "Step size, ppd:\t\t0.001\n",
    "Flow bias on:\t\tfalse\n",
    "Catalase on walls:\tfalse\n",
    "Random seed:\t\t1234\n",
    )
    correct_df = DataFrame(
        nth_trial = [n for n in 1:10],
        inner_crosstalk = [100.0, 100.0, 50.0, 50.0, 50.0,
            25.0, 50.0, 50.0, 50.0, 50.0],
        outer_crosstalk = [0.0 for _ in 1:10]
    )
    @test gd_arlett_noflow_test_df == correct_df
    @test output === correct_output
end

@testset "getdata_arlett_animation" begin
    set_FLOW_BIAS(false)
    set_CATALASE_ON_WALLS(false)
    set_MAX_STEPS_PER_WALK(50)
    output = @capture_out begin
        global gd_arlett_ani_df = getdata_arlett_animation(1234)
    end
    correct_output = string(
    "############################\n",
    "      Arlett Simulation     \n",
    "############################\n",
    "_________Parameters_________\n",
    "# of trials:\t\t1\n",
    "# of steps (max):\t50\n",
    "Step size, water:\t0.1\n",
    "Step size, enz:\t\t0.005\n",
    "Step size, ppd:\t\t0.001\n",
    "Flow bias on:\t\tfalse\n",
    "Catalase on walls:\tfalse\n",
    "Random seed:\t\t1234\n",
    "NOTE: # of steps (max) is ignored\n")
    @test output === correct_output
    @test size(gd_arlett_ani_df) === (17885, 2)
    @test names(gd_arlett_ani_df) == ["x_coordinate",
                                        "y_coordinate"]
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
"""
    gerchbergsaxton(Ain, Aout, kw...)

reconstruct the complex signals from their amplitude values based on the Gerchberg-Saxton algorithm.

## References

1. R. W. Gerchberg and W. O. Saxton, "A practical algorithm for the determination of the phase from image and diffraction plane pictures," Optik 35, 237 (1972)
"""
function gerchbergsaxton(Ain, Aout;
        maxiter=100,
        projector=(fft, ifft),
        errorfunc=msd,
        tol=eps(),
        init=randn(ComplexF64, size(Ain))
    )
    Cin = init
    @. Cin = Ain * exp(im * angle(Cin))
    Cout = similar(Aout)
    for i in 1:maxiter
        Cout = projector[1](Cin)
        @. Cout = Aout * exp(im * angle(Cout))
        Cin = projector[2](Cout)
        if errorfunc(abs.(Cin), Ain) < tol
            break
        end
        @. Cin = Ain * exp(im * angle(Cin))
    end
    (Cin, Cout)
end


text:
using Documenter, PregnancyCoagDataAnalysisKit

makedocs(sitename="Documentation",
    pages=[
        "Regression" => "regression.md",
        "Pretreatment" => "pretreatment.md"
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/varnerlab/PregnancyCoagDataAnalysisKit.jl.git",
    devbranch = "main",
    devurl = "dev",
)

max_stars_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_issues_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_forks_repo_name: varnerlab/PregnancyCoagDataAnalysisKit


text:
# Compute C = A * B, using naive matrix-multiplication algorithm,
# with a pre-allocated output array C.  ("!" is a Julia convention
# for functions that modify their arguments.)

function matmul!(C, A, B)
    m,n = size(A)
    n,p = size(B)
    size(C) == (m,p) || error("incorrect dimensions ", size(C), " ≠ $m × $p")
    for i = 1:m
        for k = 1:p
            c = zero(eltype(C))
            for j = 1:n
                @inbounds c += A[i,j] * B[j,k]
            end
            @inbounds C[i,k] = c
        end
    end
    return C
end

# a wrapper that allocates C of an appropriate type
matmul(A, B) = matmul!(Array{promote_type(eltype(A), eltype(B))}(undef,size(A,1), size(B,2)),A, B)

# Example:
# A = rand(5,6)
# B = rand(6,7)
# matmul(A,B)

max_stars_repo_path: Code/Numerical Analysis & Machine Learning/Matrix Multiplication/NaiveMatrixMultiplication.jl

max_stars_repo_name: BambooFlower/Math-Scripts

max_issues_repo_name: BambooFlower/Math-Scripts

max_forks_repo_name: BambooFlower/Math-Scripts


text:
using FastGaussQuadrature

@testset "Gauss–Hermite" begin
    tol = 1e-14

    @testset "Golub-Welsch" begin
        n = 18;
        x,w = gausshermite( n )
        @test isa(x,Vector{Float64})
        @test isa(w,Vector{Float64})
        @test (length(x) == n && length(w) == n)
        @test (dot(w,x) < tol && abs(dot(w,x.^2) - sqrt(pi)/2) < tol)
    end

    @testset "Recurrence" begin
        n = 42
        x,w = gausshermite( n )
        @test isa(x,Vector{Float64})
        @test isa(w,Vector{Float64})
        @test (length(x) == n && length(w) == n)
        @test (dot(w,x) < tol && abs(dot(w,x.^2) - sqrt(pi)/2) < tol)
        @test abs(x[37] - 5.660357581283058) < tol
        @test abs(w[17] - 0.032202101288908) < tol

        f = x -> 1 + x + x^2 + x^3
        @test w'f.(x) ≈ 2.6586807763789717
    end

    @testset "Asymptotics"  begin
        n = 251
        x,w = gausshermite( n )
        @test isa(x,Vector{Float64})
        @test isa(w,Vector{Float64})
        @test (length(x) == n && length(w) == n)
        @test (dot(w,x) < tol && abs(dot(w,x.^2) - sqrt(pi)/2) < 300*tol)
        @test abs(x[37] - -13.292221459334638) < tol
        @test abs(w[123] - 0.117419270715955) < 2*tol

        n = 500
        x,w = gausshermite( n )
        @test isa(x,Vector{Float64})
        @test isa(w,Vector{Float64})
        @test (length(x) == n && length(w) == n)
        @test (dot(w,x) < tol && abs(dot(w,x.^2) - sqrt(pi)/2) < 300*tol)
    end

    @testset "Recurrence" begin
        x = 0.1
        He0 = 1 # Probabilists He_n(x) = H_n(x/sqrt(2))/2^(n/2)
        He1 = x; He2 = x^2-1; He3 = x*(x^2-3); He4 = 3 - 6x^2 + x^4
        @test FastGaussQuadrature.hermpoly_rec(1,x)[1] ≈ He1*exp(-x^2/4)
        @test FastGaussQuadrature.hermpoly_rec(2,x)[1] ≈ He2*exp(-x^2/4)/sqrt(2)
        @test FastGaussQuadrature.hermpoly_rec(3,x)[1] ≈ He3*exp(-x^2/4)/sqrt(2*3)
        @test FastGaussQuadrature.hermpoly_rec(4,x)[1] ≈ He4*exp(-x^2/4)/sqrt(2*3*4)

        @test FastGaussQuadrature.hermpoly_rec(1,60)[1] ≈ 0.0
        @test FastGaussQuadrature.hermpoly_rec(20^2,20)[1] ≈ 0.20019391063012504
        @test FastGaussQuadrature.hermpoly_rec(60^2,60)[1] ≈ 0.07918022667865038
        @test FastGaussQuadrature.hermpoly_rec(60^2,100)[1] ≈ -0.14191347555044895
        @test FastGaussQuadrature.hermpoly_rec(60^2,130)[1] ≈ 1.7334093012299562E-73
        @test FastGaussQuadrature.hermpoly_rec(60^2,200)[1] ≈ 0.0

        @test FastGaussQuadrature.hermpoly_rec(1:9,x) ≈ [FastGaussQuadrature.hermpoly_rec(n,x)[1] for n=1:9]
        @test FastGaussQuadrature.hermpoly_rec(0:20^2,20)[end] ≈ FastGaussQuadrature.hermpoly_rec(20^2,20)[1]
    end

    @testset "Transform" begin
        n = 500
        x,w = FastGaussQuadrature.unweightedgausshermite( n )
        @test w[1] ≠ 0

        V = Array{Float64}(undef,n,n)
        for k=1:n
            V[k,:] = FastGaussQuadrature.hermpoly_rec(0:n-1, sqrt(2)*x[k])
        end

        f = x -> first(FastGaussQuadrature.hermpoly_rec(1, sqrt(2)*x))
        @test V' * (w.* f.(x)) ≈ [0; sqrt(π); zeros(n-2)]
        f = x -> first(FastGaussQuadrature.hermpoly_rec(2, sqrt(2)*x))
        @test V' * (w.* f.(x)) ≈ [0; 0; sqrt(π); zeros(n-3)]
        f = x -> first(FastGaussQuadrature.hermpoly_rec(3, sqrt(2)*x))
        @test V' * (w.* f.(x)) ≈ [0; 0; 0; sqrt(π); zeros(n-4)]
    end
end


text:
s = """
if VERSION < v"0.4-" && VERSION > v"0.1"
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4")
else
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4")
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "Reachable by 0.3")
@test contains(msgs[2].message, "Unreachable by 0.4")
@test contains(msgs[3].message, "Unreachable by 0.3")
@test contains(msgs[4].message, "Reachable by 0.4")

s = """
test() = true
if VERSION < v"0.4-" && test()
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
else
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
end
"""
msgs = lintstr(s)
@test length(msgs) == 4
@test msgs[1].code == :I271
@test contains(msgs[1].message, "Reachable by 0.3")
@test contains(msgs[2].message, "Unreachable by 0.4")
@test contains(msgs[3].message, "Reachable by 0.3")
@test contains(msgs[4].message, "Reachable by 0.4") # we cannot prove unreachable

s = """
test() = true
if VERSION < v"0.4-" || test()
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
else
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
end
"""
msgs = lintstr(s)
@test length(msgs) == 4
@test msgs[1].code == :I271
@test contains(msgs[1].message, "Reachable by 0.3")
@test contains(msgs[2].message, "Reachable by 0.4")
@test contains(msgs[3].message, "Unreachable by 0.3")
@test contains(msgs[4].message, "Reachable by 0.4")

# testing `||`, should be rare
s = """
if VERSION >= v"0.4-" || VERSION < v"0.3"
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
else
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
end
"""
msgs = lintstr(s)
@test length(msgs) == 4
@test msgs[1].code == :I271
@test contains(msgs[1].message, "Unreachable by 0.3")
@test contains(msgs[2].message, "Reachable by 0.4")
@test contains(msgs[3].message, "Reachable by 0.3")
@test contains(msgs[4].message, "Unreachable by 0.4")

s = """
if !(VERSION >= v"0.4-")
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
else
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "Reachable by 0.3")
@test contains(msgs[2].message, "Unreachable by 0.4")
@test contains(msgs[3].message, "Unreachable by 0.3")
@test contains(msgs[4].message, "Reachable by 0.4")


text:
using Unitful
using HDF5
using Mmap
using FortranFiles
using Interpolations
import PhysicalConstants.CODATA2018: h, k_B, R_∞, c_0, m_e, m_u, e, ε_0, a_0
using Transparency
using ProgressMeter
using Base.Threads
@derived_dimension NumberDensity Unitful.𝐋^-3
@derived_dimension PerLength Unitful.𝐋^-1

const Hα = AtomicLine(97489.992u"cm^-1", 82257.172u"cm^-1", 109754.578u"cm^-1",
                    18, 8, 6.411e-01, 1.008 * m_u, 1)
const γ_energy = (h * c_0 / (4 * π * Hα.λ0)) |> u"J"
const unsold_const = const_unsold(Hα)
const quad_stark_const = const_quadratic_stark(Hα)

"""
Reads atmosphere in the input format of MULTI3D. Returns a dictionary.
"""
function read_atmos_multi3d(par_file, atmos_file; dtype=Float32)
    data = Dict()
    # Get parameters and height scale
    fobj = FortranFile(par_file, "r")
    _ = read(fobj, Int32)
    data[:nx] = read(fobj, Int32)
    data[:ny] = read(fobj, Int32)
    data[:nz] = read(fobj, Int32)
    _ = read(fobj, (Float64, data[:nx]))
    _ = read(fobj, (Float64, data[:ny]))
    data[:z] = read(fobj, (Float64, data[:nz]))u"cm"
    close(fobj)
    # Get atmosphere
    fobj = open(atmos_file, "r")
    shape = (data[:nx], data[:ny], data[:nz])
    block_size = data[:nx] * data[:ny] * data[:nz] * sizeof(dtype)
    data[:ne] = Mmap.mmap(fobj, Array{typeof(one(dtype)u"cm^-3"), 3}, shape)
    data[:temperature] = Mmap.mmap(fobj, Array{typeof(one(dtype)u"K"), 3}, shape, block_size)
    data[:vz] = Mmap.mmap(fobj, Array{typeof(one(dtype)u"cm/s"), 3}, shape, block_size * 4)
    close(fobj)
    return data
end

"""
Reads NLTE populations from MULTI3D output
"""
function read_multi3d_pops(pop_file, shape; dtype=Float32)
    fobj = open(pop_file, "r")
    pops = Mmap.mmap(fobj, Array{typeof(one(dtype)u"cm^-3"), 4}, shape) # NLTE pops
    close(fobj)
    return pops
end

"""
Reads NLTE population predictions from SunnyNet. The input file is an HDF5 file
written by SunnyNet.sunnynet_predict_populations().
"""
function read_sunnynet_pops(pop_file)
    pop_nn = h5read(pop_file, "populations")
    tmp = h5readattr(pop_file, "populations/")
    cmass_mean = tmp["cmass_mean"]
    cmass_scale = tmp["cmass_scale"]
    nz = length(cmass_mean)
    pop_nn = permutedims(pop_nn, [2, 1, 3, 4])
    nhydr, ny, nx = size(pop_nn)[end-2:end]
    new_pop = zeros(Float64, (nx, ny, nz, nhydr))
    p = ProgressMeter.Progress(nx * ny * nhydr, desc="Reading populations ")
    Threads.@threads for ij in CartesianIndices((nx, ny, nhydr))
        i, j, l = Tuple(ij)
        interp = LinearInterpolation(cmass_scale, pop_nn[:, l, j, i], extrapolation_bc=Line())
        new_pop[i, j, :, l] .= interp(cmass_mean) * 1.0f-6
        ProgressMeter.next!(p)
    end
    return new_pop * u"cm^-3"
end

"""
Compute continuum extinction.
"""
function αcont(λ::Unitful.Length, temperature::Unitful.Temperature,
               electron_density::NumberDensity, h_ground_density::NumberDensity,
               proton_density::NumberDensity)
    α = Transparency.hminus_ff_stilley(λ, temperature, h_ground_density, electron_density)
    α += Transparency.hminus_bf_geltman(λ, temperature, h_ground_density, electron_density)
    α += hydrogenic_ff(c_0 / λ, temperature, electron_density, proton_density, 1)
    α += h2plus_ff(λ, temperature, h_ground_density, proton_density)
    α += h2plus_bf(λ, temperature, h_ground_density, proton_density)
    α += thomson(electron_density)
    α += rayleigh_h(λ, h_ground_density)
    return α
end

"""
Calculates full spectrum for a given 1D column
"""
function calc_Halpha_1D(
        z_scale::Array{<: Unitful.Length, 1},
        wave::Array{<: Unitful.Length, 1},
        temp::Array{<: Unitful.Temperature, 1},
        ne::Array{<: NumberDensity, 1},
        v_los::Array{<: Unitful.Velocity, 1},
        hpops::Array{<: NumberDensity, 2},
)
    nwave = length(wave)
    nz = length(z_scale)
    intensity = zeros(nwave)u"kW / (m^2 * nm)"
    α_total = zeros(nz)u"m^-1"
    source_function = zeros(nz)u"kW / (m^2 * nm)"
    # Calculate continuum opacity and wavelength-independent quantities
    α_cont = αcont.(Hα.λ0, temp, ne, hpops[:, 1], hpops[:, 6])
    j_cont = α_cont .* (blackbody_λ.(Hα.λ0, temp) .|> u"kW / (m^2 * nm)")
    γ = γ_unsold.(unsold_const, temp, hpops[:, 1]) # van der Waals broadening
    γ .+= 5.701E+08u"s^-1" # from the atom file, not Hα.Aji  # natural broadening
    γ .+= γ_linear_stark.(ne, 3, 2)
    γ .+= γ_quadratic_stark.(ne, temp, stark_constant=quad_stark_const)
    ΔλD = doppler_width.(Hα.λ0, Hα.atom_weight, temp)
    # Calculate line opacity and intensity
    for (i, λ) in enumerate(wave)
        for iz in 1:nz
            # Wavelength-dependent part
            a = damping(γ[iz], λ, ΔλD[iz])
            v = (λ - Hα.λ0 + Hα.λ0 * v_los[iz] / c_0) / ΔλD[iz]
            profile = voigt_profile(a, v, ΔλD[iz])
            # Part that only multiplies by wavelength:
            α_tmp = γ_energy * profile
            j_tmp = α_tmp
            α_tmp *= hpops[iz, 2] * Hα.Bij - hpops[iz, 3] * Hα.Bji
            j_tmp *= hpops[iz, 3] * Hα.Aji
            α_tmp += α_cont[iz]
            j_tmp += j_cont[iz]
            source_function[iz] = j_tmp / α_tmp
            α_total[iz] = α_tmp
        end
        intensity[i] = piecewise_1D_linear(z_scale, α_total, source_function)[1]
    end
    return intensity
end

"""
Performs spectral synthesis for Hα, using NLTE populations from MULTI3D
and predicted NLTE populations from SunnyNet (for comparison). Also needs
an existing 3D model atmosphere in the input format of MULTI3D. File names
and paths can be edited.

FILE_WAVES is an hdf5 file with an array of wavelengths (nm in vacuum)
to compute the line profile.

This example reads the atmosphere from the output of MULTI3D, but this
can be changed. Change the lines of read_atmos_multi3d() to read your
simulation format. 'atmos' needs to be a dictionary with the following keys
for atmospheric arrays :z, :temperature, :ne, :vz (need to have units using
Unitful), plus the auxiliary dimensions :nx, :ny, and :nz.
"""
function do_calc()
    FILE_WAVES = "wavelengths.hdf5"
    M3D_DIR = "/path/to/multi3d/output/"
    SUNNYNET_POPS = "sunnynet_predicted.hdf5"
    OUTPUT = "sunnynet_intensity.hdf5"
    waves = h5read(FILE_WAVES, "wavelength")u"nm"
    nwave = length(waves)
    atmos = read_atmos_multi3d(joinpath(M3D_DIR, "out_par"), joinpath(M3D_DIR, "out_atm"))
    nx, ny, nz = atmos[:nx], atmos[:ny], atmos[:nz]
    lo, hi = (1, nz)  # use to slice the height scale, if needed
    z_scale = atmos[:z][lo:hi] .|> u"m"
    if isfile(OUTPUT)
        println("Intensity already calculated, refusing to overwrite.")
    else
        hpops = read_nn_pops(SUNNYNET_POPS)
        intensity = zeros(nwave, nx, ny)u"kW / (m^2 * nm)"
        p = ProgressMeter.Progress(nx * ny, desc="Calculating spectra ")
        Threads.@threads for ij in CartesianIndices((nx, ny))
            ix, iy = Tuple(ij)
            intensity[:, ix, iy] = calc_Halpha_1D(
                z_scale, waves, atmos[:temperature][ix, iy, lo:hi],
                atmos[:ne][ix, iy, lo:hi], atmos[:vz][ix, iy, lo:hi] .|> u"m/s",
                hpops[ix, iy, lo:hi, :]
            )
            ProgressMeter.next!(p)
        end
        h5write(OUTPUT, "intensity", ustrip(intensity))
    end
end

do_calc()

max_stars_repo_name: tiagopereira/SunnyNet

max_issues_repo_name: tiagopereira/SunnyNet

max_forks_repo_name: tiagopereira/SunnyNet


text:
# Implementes a basis-function (Bayesian) Linear Regressor, where the basis functions are
# provided by an MLP whose parameters are learned via type-2 maximum likelihood
# i.e. maximising the log marginal likelihood of the parameters
# This can be thought of as an MLP where we perform Bayesian inference over only the last
# layer of the network.

# Exact inference is possible in the batch setting. A variational exact that enables mini-
# batching should appear at some point in the future.

# Import a load of things.
using LinearAlgebra, BayesianLinearRegressors, Zygote, Flux, Plots, ProgressMeter, Random
using BayesianLinearRegressors: BayesianLinearRegressor, logpdf, posterior, marginals
using Statistics: mean, std

# Create an MLP that you might have seen in the 90s.
Dlat = 50;
W1, b1 = randn(Dlat, 1), randn(Dlat);
ϕ = Chain(x -> reshape(x, 1, :), x -> tanh.(W1 * x .+ b1))

# Initialise the standard deviation of the observation noise. We will learn this.
logσ = [log(1)]

# Generate toy data.
Ntr = 1000;
x = collect(range(-5.0, 5.0; length=Ntr));
y = sin.(x) + 0.1 * randn(Ntr);

# Construct a simple BLR, with zero mean and identity covariance / precision prior.
blr = BayesianLinearRegressor(zeros(Dlat), Matrix{Float64}(I, Dlat, Dlat))

function nn_blr_training_loop(x, y, pars, Nitr, opt)
    tr_nlml = Vector{Float64}(undef, Nitr)
    p = ProgressMeter.Progress(Nitr)
    for itr in 1:Nitr
        nlml, back = Zygote.forward(Zygote.Params(pars)) do
            -logpdf(blr(ϕ(x), exp(2 * logσ[1])), y)
        end
        tr_nlml[itr] = nlml
        g = back(1.0)
        for par in pars
            Flux.Optimise.update!(opt, par, g[par])
        end
        showvalues = [(:itr, itr), (:nlml, nlml), (:σ_ε, exp(logσ[1]))]
        ProgressMeter.next!(p; showvalues=showvalues)
    end
    return tr_nlml
end

# Gather the parameters and perform 500 iterations of ADAM.
pars = [W1, b1, logσ];
nlmls = nn_blr_training_loop(x, y, pars, 500, ADAM(1e-2, (0.9, 0.999)));

# Compute the posterior BLR, given the optimised network and noise parameters.
blr′ = posterior(blr(ϕ(x), exp(2 * logσ[1])), y);

# Make posterior predictions.
xte = collect(range(-10, 10; length=1000));
ypr_te = marginals(blr′(ϕ(xte), exp(2 * logσ[1])));

# Generate samples from the posterior over the latent function. Note that the observation
# noise is small.
y_samples_te = rand(MersenneTwister(123456), blr′(ϕ(xte), eps()), 25);

# Plot the predicted function vs the ground-truth. Doesn't generalise well, unsurprisingly.
# Does do _really_ well inside the region containing the training data though, and is
# nicely callibrated.
plt = plot();

# Plot the true function with aleatoric uncertainty due to observation noise.
plot!(plt, xte, sin.(xte); linecolor="purple", linewidth=1.5, label="sin");
plot!(plt, xte, sin.(xte) .+ 0.3; linecolor="purple", linewidth=1.5, label="");
plot!(plt, xte, sin.(xte) .- 0.3; linecolor="purple", linewidth=1.5, label="")

# Visualise the posterior marginal uncertainty via 3σ error bars.
m_te, σ_te = mean.(ypr_te), std.(ypr_te)
plot!(
    plt,
    xte,
    [m_te m_te];
    label="",
    fillrange=[m_te .+ 3 .* σ_te, m_te .- 3 .* σ_te],
    linewidth=0,
    fillalpha=0.3,
    fillcolor="red",
);

# Visualise the posterior distribution over the latent function via samples.
plot!(
    plt, xte, y_samples_te[:, 1]; linecolor="blue", linewidth=0.1, label="Posterior Samples"
);
plot!(plt, xte, y_samples_te; linecolor="blue", linewidth=0.1, label="");

# Plot the data.
scatter!(plt, x, y; markercolor="red", label="y", markersize=0.1);
display(plt);

# Compute and display residuals
ypr = marginals(blr′(ϕ(x), exp(2 * logσ[1])));
ε = y .- mean.(ypr);
@show mean(ε), std(ε);
histogram(ε; bins=50)

# Compute and display training curve.
plot(nlmls)


text:
#Tests on the networks containing hill functions.
for rn in reaction_networks_hill
    p_vals = map(amp->rand(1:amp,length(rn.params)),[2,3,6])
    for p in p_vals
        fix_parameters(rn,p)
        fps = steady_states(rn,Vector{Float64}())
        stab = map(fp->stability(fp,rn,Float64.(p)), fps)
        for i=1:length(fps)
            stab[i] ? (tend = 10) : (tend = -10)
            end_point = OrdinaryDiffEq.solve(ODEProblem(rn,fps[i],(0.,tend),p),Rosenbrock23())[end]
            @test maximum(abs.(end_point-fps[i]))<0.0001
        end
        #start_points = map(amp->amp*rand(length(rn.syms)),[1.,10.,100.])
        #for sp in start_points
        #    end_point = OrdinaryDiffEq.solve(ODEProblem(rn,sp,(0.,1000.),p),Rosenbrock23())[end]
        #    @test any(map(fp->maximum(abs.(fp-end_point)),fps) .< 0.001)
        #end
    end
end


text:
# Load file for needed functions
include("../Packages/ReciprocalVariety.jl")

# Basic setup. For ternary quartics: r=2, s=3
    r=2
    s=3
    m=binomial(r+s-1,s-1)
    N=binomial(2r+s-1,2r)
    M=binomial(m + 1, 2)
    @var x[0:2r, 0:2r, 0:2r],y[1:M]

# Define ideal equations
    cat=CatalecticantIndices(r,s,m)
    adj=adjugate(Expression.(vec_to_sym(y,m)));
    eq=CatalecticantSpace(cat,adj,m);

# Set up system
    F = System(eq; variables=y);

# Compute point on the variety
    x₀ = randn(ComplexF64,N)
    @var x[0:2r, 0:2r, 0:2r]
    X = map(ijk -> x[(ijk .+ 1)...], cat)
    y₀ = sym_to_vec(inv(X(variables(X) => x₀)),m)

# Put affine subspace through this point
    V₀ = let
       A₀ = randn(ComplexF64, N, M)
       b₀ = A₀ * y₀
       LinearSubspace(A₀, b₀)
    end

# Perform monodromy
    mres = monodromy_solve(
       F,
       y₀,
       V₀,
       parameter_sampler=_ -> LinearSubspace(
           randn(ComplexF64, N, M),
           randn(ComplexF64, N),
       ),
   )

max_stars_repo_name: LauraBMo/InvertingCats

max_issues_repo_name: LauraBMo/InvertingCats

max_forks_repo_name: LauraBMo/InvertingCats


text:
using NimrodFiles
using Base.Test

@testset "NimrodFiles" begin
    include("./io.jl")
    include("./asciigrid.jl")
end


text: @enum Fruits APPLE BANANA CHERRY

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# f(x,y)=x'-y

FXtransposeMY(x::Float64,y::Float64)=Fxmy(x,y)
FXtransposeMY(x::Float64,y)=Fxmy(x,y)
FXtransposeMY(x,y::Float64)=Fxmy(x',y)
FXtransposeMY(x,y)=Fxmy(x',y)
FXtransposeMY_inplace(value,auxvalue,x,y)=Fxmy_inplace(value,auxvalue,x',y)

function DXtransposeMY(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::Array,y::Array)
    if derivativeIDX==1
        if size(x)==(1,1)
            axpy!(1.0,[sum(grad_c)],grad_n)
        else
            axpy!(1.0,grad_c',grad_n)
        end
    elseif derivativeIDX==2
        if size(y)==(1,1)
            axpy!(-1.0,[sum(grad_c)],grad_n)
        else
            axpy!(-1.0,grad_c,grad_n)
        end
    end
end

if PROC=="GPU"

    function FXtransposeMY_inplace(value::CudaArray,auxvalue,x::CudaArray,y::CudaArray)
        if size(x)==(1,1)
            gfill!(value,x); axpy!(-1.0,y,value)
        elseif size(y)==(1,1)
            gfill!(value,y)
            CUBLAS.geam!('T','N',1.0,x,-1.0,value,value)
        else
            CUBLAS.geam!('T','N',1.0,x,-1.0,y,value)
        end
    end

    function DXtransposeMY(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray,y::CudaArray)
        if derivativeIDX==1
            if size(x)==(1,1)
                axpy!(1.0,sum(grad_c),grad_n)
            else
                tmp=CudaArray(Float64,size(grad_c))
                CUBLAS.geam!('T','N',1.0,grad_c,0.0,tmp,tmp) # tmp=grad_c'
                axpy!(1.0,tmp,grad_n); free(tmp)
            end
        elseif derivativeIDX==2
            if size(y)==(1,1)
                axpy!(-1.0,sum(grad_c),grad_n)
            else
                axpy!(-1.0,grad_c,grad_n)
            end
        end
    end

end

Derivative[FXtransposeMY]=DXtransposeMY
Inplace[FXtransposeMY]=FXtransposeMY_inplace

export FXtransposeMY

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
#=
Euler discovered the remarkable quadratic formula:

n² + n + 41

It turns out that the formula will produce 40 primes for the consecutive values n = 0 to 39. However, when n = 40, 402 + 40 + 41 = 40(40 + 1) + 41 is divisible by 41, and certainly when n = 41, 41² + 41 + 41 is clearly divisible by 41.

The incredible formula  n² − 79n + 1601 was discovered, which produces 80 primes for the consecutive values n = 0 to 79. The product of the coefficients, −79 and 1601, is −126479.

Considering quadratics of the form:

n² + an + b, where |a| < 1000 and |b| < 1000

where |n| is the modulus/absolute value of n
e.g. |11| = 11 and |−4| = 4
Find the product of the coefficients, a and b, for the quadratic expression that produces the maximum number of primes for consecutive values of n, starting with n = 0.
=#

using Primes

function longestprimerun(a,b)
  i = 0
  while isprime((i^2)+(a*i)+b)
    i+= 1
  end
  i
end

function calc()
  (mina, maxa) = (-1000, 1000)
  (largesta, largestb) = (0,0)
  largest = 0
  for a in mina:maxa, b in mina:maxa
    primerun = longestprimerun(a,b)
    if primerun > largest
      (largesta, largestb) = (a,b)
      largest = primerun
    end
  end
  largesta * largestb
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
module PlottingUtils

function reshapedata(A)
        m,n = size(A)
        sz = int(sqrt(m))
        A -= mean(A)

        gridsize = int(ceil(sqrt(n)))
        buffer = 1
        griddata = ones(gridsize*(sz+1)+1,gridsize*(sz+1)+1)
        index = 1
        for i = 1:gridsize
                for j = 1:gridsize
                        if index > n
                                continue
                        end
                        columnlimit = maximum(abs(A[:,index]))
                        griddata[buffer+(i-1)*(sz+buffer)+(1:sz),buffer+(j-1)*(sz+buffer)+(1:sz)] = reshape(A[:,index],sz,sz)/columnlimit
                        index += 1
                end
        end

	return griddata
end

import Plotly
export displaynetwork_plotly
function displaynetwork_plotly(A,filename,username,userkey)
	griddata = reshapedata(A)
	Plotly.signin(username, userkey)
        data = [
          [
            "z" => griddata,
            "colorscale" => "Greys",
            "type" => "heatmap"
          ]
        ]
        layout = [
            "autosize" => false,
            "width" => 500,
            "height"=> 500
        ]
        response = Plotly.plot(data, ["layout" => layout, "filename" => filename, "fileopt" => "overwrite"])
        plot_url = response["url"]
end

import Gadfly
export displaynetwork_gadfly
function displaynetwork_gadfly(A)
	griddata = reshapedata(A)
	Gadfly.spy(A)
end

import Winston
export displaynetwork_winston
function displaynetwork_winston(A)
	griddata = reshapedata(A)
	p = Winston.FramedPlot()
	Winston.colormap("grays")
	Winston.add(p,Winston.imagesc(griddata))
	display(p)
end

end


text:
@testset "CollapsedDim" begin
    using NeuralAttentionlib.Matmul
    x = randn(7,6,5,4,3,2)

    @test collapsed_size(x, 3, 5) == (42, 20, 6)
    @test collapsed_size(x, 2, 3) == (7, 6, 120)
    @test collapsed_size(x, 4, 8) == (210, 24, 1)
    @test noncollapsed_size(x, 3, 5, 1) == (7, 6)
    @test noncollapsed_size(x, 3, 5, 2) == (5, 4)
    @test noncollapsed_size(x, 3, 5, 3) == (3, 2)
    @test noncollapsed_size(x, 2, 3, 1) == (7,)
    @test noncollapsed_size(x, 2, 3, 2) == (6,)
    @test noncollapsed_size(x, 2, 3, 3) == (5, 4, 3, 2)
    @test noncollapsed_size(x, 4, 8, 1) == (7, 6, 5)
    @test noncollapsed_size(x, 4, 8, 2) == (4, 3, 2, 1)
    @test noncollapsed_size(x, 4, 8, 3) == (1,)

    @test reshape(x, (42, 20, 6)) == collapseddim(x, 3, 5)
    @test reshape(x, (7, 6, 120)) == collapseddim(x, 2, 3)
    @test reshape(x, (210, 24, 1)) == collapseddim(x, 4, 8)

    @test reshape(x, (42, 20, 6)) == CollapsedDimArray(x, 3, 5)
    @test reshape(x, (7, 6, 120)) == CollapsedDimArray(x, 2, 3)
    @test reshape(x, (210, 24, 1)) == CollapsedDimArray(x, 4, 8)

    y = rand(5, 4)

    @test collapsed_size(y, 2, 3) == (5, 4, 1)
    @test collapsed_size(y, 3, 5) == (20, 1, 1)
    @test noncollapsed_size(y, 2, 3, 1) == (5,)
    @test noncollapsed_size(y, 3, 5, 1) == (5, 4)
    @test noncollapsed_size(y, 2, 3, 2) == (4,)
    @test noncollapsed_size(y, 3, 5, 2) == (1, 1)
    @test noncollapsed_size(y, 2, 3, 3) == (1,)
    @test noncollapsed_size(y, 3, 5, 3) == (1,)

    @testset "AD" begin
        x = randn(7,6,5,4,3,2)

        test_rrule(parent, CollapsedDimArray(randn(6)))
        test_rrule(parent, CollapsedDimArray(x, 3, 5))
        test_rrule(parent, CollapsedDimArray(x, 2, 3))
        test_rrule(parent, CollapsedDimArray(x, 4, 8))

        test_rrule(CollapsedDimArray, randn(6), (6, 1, 1),
                   static(2), static(3), static(true))
        test_rrule(CollapsedDimArray, x, (42, 20, 6),
                   static(3), static(5), static(false))
        test_rrule(CollapsedDimArray, x, (7, 6, 120),
                   static(2), static(3), static(false))
        test_rrule(CollapsedDimArray, x, (210, 24, 1),
                   static(4), static(8), static(true))

    end

end


text:
using SubsetSelection
using StatsBase, Test

n = 500; p = 1000; k = 10;

indices = sort(sample(1:p, StatsBase.Weights(ones(p)/p), k, replace=false));
w = sample(-1:2:1, k);
X = randn(n,p); Y = X[:,indices]*w;
@time Sparse_Regressor = subsetSelection(OLS(), Constraint(k), Y, X)
for i in 1:k
  @test Sparse_Regressor.indices[i]==indices[i]
end


text:
@latexrecipe function f(
    r::DiffEqBase.AbstractReactionNetwork;
    bracket=false,
    noise=false,
    noise_only=false,
    noise_var=:W,
    env=:align,
)

    env := env
    ## DiffEqBio's chemical arrow notation support is not yet migrated outside
    ## of Latexify.jl. If used, simply pass the reaction network object on and
    ## set the appropriate kwargs.
    if env in [:arrows, :chem, :chemical, :arrow]
        bracket := bracket
        return r
    end

    ### Get the noise terms.
    if noise || noise_only
        noise_matrix = reshape(r.g_func, length(r.syms), :)
        w = [Symbol("d$(noise_var)_$(i)(t)") for _ in 1:size(noise_matrix, 1), i in 1:size(noise_matrix,2)]
        noise_matrix = Expr.(:call, :*, noise_matrix, w)
        noise_rhs = Meta.parse.([join(noise_matrix[i,:], " + ") for i in 1:size(noise_matrix,1)])

        ### Clean out terms which are either zero or multiplied by zero.
        for i in 1:length(noise_rhs)
            filter!(x -> x != 0, noise_rhs[i].args)
            filter!(x -> hasfield(typeof(x), :args) ? x.args[1:2] != [:*, 0] : true, noise_rhs[i].args)
        end
    end

    ### Define the left-hand side
    if noise || noise_only
        if bracket
            lhs = [Meta.parse("d[$x](t)") for x in r.syms]
        else
            lhs = [Meta.parse("d$x(t)") for x in r.syms]
        end
    else
        if bracket
            lhs = [Latexify.LaTeXString("\\frac{d\\left[ $x \\right](t)}{dt}") for x in r.syms]
        else
            lhs = [Latexify.LaTeXString("\\frac{d$x(t)}{dt}") for x in r.syms]
        end
    end

    ### Define the right-hand side
    if noise_only
        rhs = noise_rhs
        separator --> " ∝& "
    elseif noise
        rhs = Expr.(:call, :+, Expr.(:call, :*, r.f_func, :dt), noise_rhs)
    else
        rhs = r.f_func
    end
    bracket && (rhs = Latexify.add_brackets(rhs, r.syms))

    return lhs, rhs
end


text:
## %PREDICTMD_GENERATED_BY%

using PredictMDExtra
PredictMDExtra.import_all()

using PredictMD
PredictMD.import_all()

# PREDICTMD IF INCLUDE TEST STATEMENTS
logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
if isnothing(logger)
    logger_stream = devnull
else
    logger_stream = logger.stream
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### Begin project-specific settings

DIRECTORY_CONTAINING_THIS_FILE = @__DIR__
PROJECT_DIRECTORY = dirname(
    joinpath(splitpath(DIRECTORY_CONTAINING_THIS_FILE)...)
    )
PROJECT_OUTPUT_DIRECTORY = joinpath(
    PROJECT_DIRECTORY,
    "output",
    )
mkpath(PROJECT_OUTPUT_DIRECTORY)
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "data"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "models"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "plots"))

# PREDICTMD IF INCLUDE TEST STATEMENTS
@debug("PROJECT_OUTPUT_DIRECTORY: ", PROJECT_OUTPUT_DIRECTORY,)
if PredictMD.is_travis_ci()
    PredictMD.cache_to_path!(
        ;
        from = ["cpu_examples", "breast_cancer_biopsy", "output",],
        to = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### End project-specific settings

### Begin random forest classifier code

Random.seed!(999)

trainingandtuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_features_df.csv",
    )
trainingandtuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_labels_df.csv",
    )
testing_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_features_df.csv",
    )
testing_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_labels_df.csv",
    )
training_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_features_df.csv",
    )
training_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_labels_df.csv",
    )
tuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_features_df.csv",
    )
tuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_labels_df.csv",
    )
trainingandtuning_features_df = DataFrames.DataFrame(
    FileIO.load(
        trainingandtuning_features_df_filename;
        type_detect_rows = 100,
        )
    )
trainingandtuning_labels_df = DataFrames.DataFrame(
    FileIO.load(
        trainingandtuning_labels_df_filename;
        type_detect_rows = 100,
        )
    )
testing_features_df = DataFrames.DataFrame(
    FileIO.load(
        testing_features_df_filename;
        type_detect_rows = 100,
        )
    )
testing_labels_df = DataFrames.DataFrame(
    FileIO.load(
        testing_labels_df_filename;
        type_detect_rows = 100,
        )
    )
training_features_df = DataFrames.DataFrame(
    FileIO.load(
        training_features_df_filename;
        type_detect_rows = 100,
        )
    )
training_labels_df = DataFrames.DataFrame(
    FileIO.load(
        training_labels_df_filename;
        type_detect_rows = 100,
        )
    )
tuning_features_df = DataFrames.DataFrame(
    FileIO.load(
        tuning_features_df_filename;
        type_detect_rows = 100,
        )
    )
tuning_labels_df = DataFrames.DataFrame(
    FileIO.load(
        tuning_labels_df_filename;
        type_detect_rows = 100,
        )
    )

smoted_training_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "smoted_training_features_df.csv",
    )
smoted_training_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "smoted_training_labels_df.csv",
    )
smoted_training_features_df = DataFrames.DataFrame(
    FileIO.load(
        smoted_training_features_df_filename;
        type_detect_rows = 100,
        )
    )
smoted_training_labels_df = DataFrames.DataFrame(
    FileIO.load(
        smoted_training_labels_df_filename;
        type_detect_rows = 100,
        )
    )

categorical_feature_names_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "categorical_feature_names.jld2",
    )
continuous_feature_names_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "continuous_feature_names.jld2",
    )
categorical_feature_names = FileIO.load(
    categorical_feature_names_filename,
    "categorical_feature_names",
    )
continuous_feature_names = FileIO.load(
    continuous_feature_names_filename,
    "continuous_feature_names",
    )
feature_names = vcat(categorical_feature_names, continuous_feature_names)

single_label_name = :Class
negative_class = "benign"
positive_class = "malignant"
single_label_levels = [negative_class, positive_class]

categorical_label_names = Symbol[single_label_name]
continuous_label_names = Symbol[]
label_names = vcat(categorical_label_names, continuous_label_names)

feature_contrasts = PredictMD.generate_feature_contrasts(
    smoted_training_features_df,
    feature_names,
    )

random_forest_classifier =
    PredictMD.single_labelmulticlassdataframerandomforestclassifier(
        feature_names,
        single_label_name,
        single_label_levels;
        nsubfeatures = 4,
        ntrees = 200,
        package = :DecisionTree,
        name = "Random forest",
        feature_contrasts = feature_contrasts,
        )

PredictMD.fit!(
    random_forest_classifier,
    smoted_training_features_df,
    smoted_training_labels_df,
    )

random_forest_classifier_hist_training =
    PredictMD.plotsinglelabelbinaryclassifierhistogram(
        random_forest_classifier,
        smoted_training_features_df,
        smoted_training_labels_df,
        single_label_name,
        single_label_levels,
        );

# PREDICTMD IF INCLUDE TEST STATEMENTS
filename = string(
    tempname(),
    "_",
    "random_forest_classifier_hist_training",
    ".pdf",
    )
rm(filename; force = true, recursive = true,)
@debug("Attempting to test that the file does not exist...", filename,)
Test.@test(!isfile(filename))
@debug("The file does not exist.", filename, isfile(filename),)
PredictMD.save_plot(filename, random_forest_classifier_hist_training)
if PredictMD.is_force_test_plots()
    @debug("Attempting to test that the file exists...", filename,)
    Test.@test(isfile(filename))
    @debug("The file does exist.", filename, isfile(filename),)
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

display(random_forest_classifier_hist_training)
PredictMD.save_plot(
    joinpath(
        PROJECT_OUTPUT_DIRECTORY,
        "plots",
        "random_forest_classifier_hist_training.pdf",
        ),
    random_forest_classifier_hist_training,
    )

random_forest_classifier_hist_testing =
    PredictMD.plotsinglelabelbinaryclassifierhistogram(
        random_forest_classifier,
        testing_features_df,
        testing_labels_df,
        single_label_name,
        single_label_levels,
        );

# PREDICTMD IF INCLUDE TEST STATEMENTS
filename = string(
    tempname(),
    "_",
    "random_forest_classifier_hist_testing",
    ".pdf",
    )
rm(filename; force = true, recursive = true,)
@debug("Attempting to test that the file does not exist...", filename,)
Test.@test(!isfile(filename))
@debug("The file does not exist.", filename, isfile(filename),)
PredictMD.save_plot(filename, random_forest_classifier_hist_testing)
if PredictMD.is_force_test_plots()
    @debug("Attempting to test that the file exists...", filename,)
    Test.@test(isfile(filename))
    @debug("The file does exist.", filename, isfile(filename),)
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

display(random_forest_classifier_hist_testing)
PredictMD.save_plot(
    joinpath(
        PROJECT_OUTPUT_DIRECTORY,
        "plots",
        "random_forest_classifier_hist_testing.pdf",
        ),
    random_forest_classifier_hist_testing,
    )

show(
    logger_stream, PredictMD.singlelabelbinaryclassificationmetrics(
        random_forest_classifier,
        smoted_training_features_df,
        smoted_training_labels_df,
        single_label_name,
        positive_class;
        sensitivity = 0.95,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

show(
    logger_stream, PredictMD.singlelabelbinaryclassificationmetrics(
        random_forest_classifier,
        testing_features_df,
        testing_labels_df,
        single_label_name,
        positive_class;
        sensitivity = 0.95,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

random_forest_classifier_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "random_forest_classifier.jld2",
    )

PredictMD.save_model(
    random_forest_classifier_filename,
    random_forest_classifier,
    )

# PREDICTMD IF INCLUDE TEST STATEMENTS
random_forest_classifier = nothing
Test.@test isnothing(random_forest_classifier)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### End random forest classifier code

# PREDICTMD IF INCLUDE TEST STATEMENTS
if PredictMD.is_travis_ci()
    PredictMD.path_to_cache!(
        ;
        to = ["cpu_examples", "breast_cancer_biopsy", "output",],
        from = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

## %PREDICTMD_GENERATED_BY%


text:
export sample!
using StatsBase: sample, Weights

"""
    sample!(model::ABM, n [, weight]; kwargs...)

Replace the agents of the `model` with a random sample of the current agents with
size `n`.

Optionally, provide a `weight`: Symbol (agent field) or function (input agent
out put number) to weight the sampling.
This means that the higher the `weight` of the agent, the higher the probability that
this agent will be chosen in the new sampling.

# Keywords
* `replace = true` : whether sampling is performed with replacement, i.e. all agents can
be chosen more than once.
* `rng = GLOBAL_RNG` : a random number generator to perform the sampling with.

See the Wright-Fisher example in the documentation for an application of `sample!`.
"""
function sample!(model::ABM{A, S}, n::Int, weight=nothing; replace=true,
  rng::AbstractRNG=Random.GLOBAL_RNG) where{A, S}

  nagents(model) > 0 || return

  org_ids = collect(keys(model.agents))
  if weight != nothing
    weights = Weights([get_data(a, weight) for a in values(model.agents)])
    newids = sample(rng, org_ids, weights, n, replace=replace)
  else
    newids = sample(rng, org_ids, n, replace=replace)
  end

  n = nextid(model)
  for id in org_ids
    if !in(id, newids)
      kill_agent!(model.agents[id], model)
    else
      noccurances = count(x->x==id, newids)
      for t in 2:noccurances
        newagent = deepcopy(model.agents[id])
        newagent.id = n
        add_agent!(newagent, model)
        n += 1
      end
    end
  end
end


text:
# three dimension vector
import Base:+,-,*,/,^

struct Vec3{T<:Real}
    x::T
    y::T
    z::T
    Vec3{T}(x::T, y::T, z::T) where {T<:Real} = new(x, y, z)
end

Vec3(x::T, y::T, z::T) where {T<:Real} = Vec3{T}(x, y, z)
Vec3(x, y, z) = Vec3(promote(x, y, z)...)
vec3(x, y, z) = Vec3(x, y, z)
vec3_from_spherical(r::Float64, θ::Float64, φ::Float64) = Vec3(
    r*sin(θ)*cos(φ), r*sin(θ)*sin(φ), r*cos(θ)
)

Base.show(io::IO, v::Vec3) = begin
    print(io, "($(v.x), $(v.y), $(v.z))")
end

Base.eltype(::Vec3{T}) where {T<:Real} = T
Base.float(v::Vec3) = Vec3{Float64}(float(v.x), float(v.y), float(v.z))

+(v::Vec3) = v
-(v::Vec3) = Vec3(-v.x, -v.y, -v.z)
+(v1::Vec3, v2::Vec3) = Vec3(v1.x+v2.x, v1.y+v2.y, v1.z+v2.z)
-(v1::Vec3, v2::Vec3) = Vec3(v1.x-v2.x, v1.y-v2.y, v1.z-v2.z)
*(v1::Vec3, v2::Vec3) = v1.x*v2.x + v1.y*v2.y + v1.z*v2.z
*(k::Real, v::Vec3) = Vec3(k*v.x, k*v.y, k*v.z)
*(v::Vec3, k::Real) = k * v
/(v::Vec3, k::Real) = Vec3(v.x/k, v.y/k, v.z/k)
×(v1::Vec3, v2::Vec3) = Vec3(
    v1.y * v2.z - v1.z * v2.y,
    v1.z * v2.x - v1.x * v2.z,
    v1.x * v2.y - v1.y * v2.x
)

Base.abs2(v::Vec3) = v.x^2 + v.y^2 + v.z^2
Base.abs(v::Vec3) = hypot(v.x, v.y, v.z)
Base.length(v::Vec3) = abs(v)
distance(v1::Vec3, v2::Vec3) = abs(v1 - v2)

Base.getindex(v::Vec3, i::Integer) = begin
    i == 1 && return v.x
    i == 2 && return v.y
    i == 3 && return v.z
    error("out of range for Vec3 getindex i = $i")
end


text:
struct tblis_tensor{T}
    type::Int32
    conj::Int32
    attr::ComplexF64
    _data::Ptr{Cvoid}
    ndim::Cuint
    _len::Ptr{Cvoid}
    _stride::Ptr{Cvoid}
end

mutable struct TTensor{T}
    tensor::tblis_tensor{T}
    n::Cuint
    data::Array{T}
    len::Array{Int}
    stride::Array{Int}
end

function TTensor(D)
    TTensor(D,1.0)
end
function TTensor(D,s)
    TTensor{Float64}(D,s)
end
function TTensor{T}(D) where T <: AbstractFloat
    TTensor{T}(D,1.0)
end
"""
    TTensor{T}(D,s) where T <: AbstractFloat

Wraps a Julia array `D` into a TTensor that can be added, multiplied, etc, scaled by s.
"""
function TTensor{T}(D,s) where T <: AbstractFloat
    strides = [1]
    lens = collect(size(D))
    for (i,v) in enumerate(lens[1:end-1])
        push!(strides,v*strides[i])
    end
    n::UInt32 = length(size(D))
    M = TTensor{T}(tblis_tensor{T}(zero(Int32),zero(Int32),
                                   0.0 + 0.0im,
                                   pointer(D),
                                   n,
                                   pointer(lens),
                                   pointer(strides)), n, D, lens, strides)
    if T == Float32
        tblis_init_tensor_scaled = dlsym(tblis,:tblis_init_tensor_scaled_s)
        ccall(tblis_init_tensor_scaled, Cvoid, (Ref{TTensor{T}},Cfloat,Cuint,Ptr{Int},Ptr{T}, Ptr{Int}),
             M, s, n, lens, D, strides)
    elseif T == Float64
        tblis_init_tensor_scaled = dlsym(tblis,:tblis_init_tensor_scaled_d)
        ccall(tblis_init_tensor_scaled, Cvoid, (Ref{TTensor{T}},Cdouble,Cuint,Ptr{Int},Ptr{T}, Ptr{Int}),
             M, s, n, lens, D, strides)
    else
        error("Type $T is not supported by TBLIS :(")
    end
    return M
end


text:
# Mel generalized cepstrum analysis
# re-coded from SPTK

# Overview:
# Tokuda, Keiichi, et al. "Mel-generalized cepstral analysis-a unified approach
# to speech spectral estimation." ICSLP. 1994.
# Details:
# 徳田恵一, et al. "メル一般化ケプストラム分析による音声のスペクトル推定."
# 電子情報通信学会論文誌 A 75.7 (1992): 1124-1134.

# memo(ryuichi):
# Definition of generalized log function:
# sᵧ(ω) =
# (ωᵞ - 1)/γ, 0 < |γ| ≤ 1,
# log(ω), γ = 0.
#
# Definition of inverse generalized log function:
# sᵧ⁻¹(ω) =
# (1 + γω)¹/ᵞ, 0 < |γ| ≤ 1,
# exp(ω), γ = 0.
#
# All-pass filter model:
# z̃⁻¹ = Ψ(z) = (z⁻¹ - α)/(1 - αz⁻¹)
#
# A spectral model based on a mel generalized cepstrum:
# H(z) = sᵧ⁻¹(∑ₘ cₐ,ᵧ(m)z̃⁻ᵐ) =
# (1 + γ∑ₘ cₐ,ᵧ(m)z̃⁻ᵐ)¹/ᵞ, 0 < |γ| ≤ 1,
# exp (∑ₘ cₐ,ᵧ(m) z̃⁻ᵐ), γ = 0.
#
# Another representation with explicit gain `K`:
# H(z) = sᵧ⁻¹(∑ₘ bₐ,ᵧ(m)Φₘ(z)) = K⋅D(z)
# where
# K = sᵧ⁻¹(bₐ,ᵧ(0)),
# D(z) =
#  (1 + γ∑ₘ bₐ,ᵧ'(m)z̃⁻ᵐ)¹/ᵞ, 0 < |γ| ≤ 1,
#  exp (∑ₘ bₐ,ᵧ'(m) z̃⁻ᵐ), γ = 0,
#
# Gain normalization (e.g. b -> b'):
# bₐ,ᵧ'(m) = bₐ,ᵧ(m) / (1 + γbₐ,ᵧ(0))

!isdefined(Base, :FFTW) && using FFTW

function ptrans!{T}(p::AbstractVector{T},  m::Int, α::AbstractFloat)
    d, o = zero(T), zero(T)

    d = p[m+1]
    @inbounds for i=m-1:-1:1
        o = p[i+1] + α * d
        d = p[i+1]
        p[i+1] = o
    end

    o = α * d
    p[1] = (one(T) - α*α) * p[1] + 2o

    p
end

function qtrans!{T}(q::AbstractVector{T},  m::Int, α::AbstractFloat)
    d = q[2]
    @inbounds for i=2:2m
        o = q[i+1] + α*d
        d = q[i+1]
        q[i+1] = o
    end

    q
end

function gain{T}(er::AbstractVector{T}, c::AbstractVector{T}, m::Int,
                 g::AbstractFloat)
    t = zero(T)
    if g != zero(T)
        for i=2:m+1
            @inbounds t += er[i] * c[i]
        end
        return er[1] + g*t
    else
        return er[1]
    end
end

function newton!{T}(c::AbstractVector{T}, # mel-generalized cepstrum stored
                    x::AbstractVector{T}, # modified periodogram
                    order::Int,           # order of cepstrum
                    α::AbstractFloat,     # allpass constant
                    γ::AbstractFloat,     # parameter of generalized log function
                    n::Int,               # the order of recursion
                    iter::Int,            # current iter #
                    y_fft::Array{Complex{T}}, # *length must be equal to length(x)*
                    z_fft::Array{Complex{T}}, # *length must be equal to length(x)*
                    bplan::FFTW.Plan,     # FFTW.Plan used in backward fft
                    cr::Vector{T} = zeros(T, length(x)),
                    pr::Vector{T} = zeros(T, length(x)),
                    rr::Vector{T} = zeros(T, length(x)),
                    ri::Vector{T} = zeros(T, length(x)),
                    qr::Vector{T} = zeros(T, length(x)),
                    qi::Vector{T} = zeros(T, length(x)),
                    Tm::Matrix{T} = Array{T,2}(order, order),
                    Hm::Matrix{T} = Array{T,2}(order, order),
                    Tm_plus_Hm::Matrix{T} = Array{T,2}(order, order),
                    b::Vector{T} = Array{T,1}(order), # right side of equation Ax = b
    )
    @assert length(x) > length(c)
    @assert n < length(x)
    @assert length(y_fft) == length(x)
    @assert length(z_fft) == length(x)

    copy!(cr, 2, c, 2, order)

    if α != zero(T)
        b2c!(view(cr, 1:n+1), cr[1:order+1], -α)
    end

    y = fft(cr)

    γ⁻¹ = one(T)/γ
    if γ == -one(T)
        copy!(pr, 1, x, 1, length(x))
    elseif γ == zero(T)
        for i=1:length(x)
            @inbounds pr[i] = x[i] / exp(2real(y[i]))
        end
    else
        @inbounds for i=1:length(x)
            tr = one(T) + γ*real(y[i])
            ti = γ*imag(y[i])
            trr, tii = tr*tr, ti*ti
            s = trr + tii
            t = x[i] * s^(-γ⁻¹)
            t /= s
            pr[i] = t
            rr[i] = tr * t
            ri[i] = ti * t
            t /= s
            qr[i] = (trr - tii) * t
            s = tr * ti * t
            qi[i] = s + s
        end
    end

    copy!(y_fft, pr)
    A_mul_B!(z_fft, bplan, y_fft)
    scale!(z_fft, 1 / length(y_fft))
    pr = real(z_fft)
    # FFTW.execute(bplan.plan, y_fft, pr)
    # scale!(pr, FFTW.normalization(pr))

    if α != zero(T)
        b2c!(view(pr, 1:2order+1), pr[1:n+1], α)
    end

    if γ == zero(T) || γ == -one(T)
        copy!(qr, 1, pr, 1, 2order+1)
        copy!(rr, 1, pr, 1, order+1)
    else
        for i=1:length(qr)
            @inbounds y_fft[i] = Complex(qr[i], qi[i])
        end

        A_mul_B!(z_fft, bplan, y_fft)
        scale!(z_fft, 1 / length(y_fft))
        qr = real(z_fft)
        # FFTW.execute(bplan.plan, y_fft, qr)
        # scale!(qr, FFTW.normalization(qr))
        for i=1:length(rr)
            @inbounds y_fft[i] = Complex(rr[i], ri[i])
        end
        A_mul_B!(z_fft, bplan, y_fft)
        scale!(z_fft, 1 / length(y_fft))
        rr = real(z_fft)
        # FFTW.execute(bplan.plan, y_fft, rr)
        # scale!(rr, FFTW.normalization(rr))

        if α != zero(T)
            b2c!(view(qr, 1:n+1), qr[1:n+1], α)
            b2c!(view(rr, 1:order+1), rr[1:n+1], α)
        end
    end

    if α != zero(T)
        ptrans!(pr, order, α)
        qtrans!(qr, order, α)
    end

    ϵ = zero(T)
    if γ != -one(T)
        ϵ = gain(rr, c, order, γ)
        c[1] = √ϵ
    end

    if γ == -one(T)
        fill!(qr, zero(T))
    elseif γ != zero(T)
        for i=3:2order+1
            @inbounds qr[i] *= one(T) + γ
        end
    end

    te = view(pr, 1:order)
    fill_toeplitz!(Tm, te)
    he = view(qr, 3:2order+1)
    fill_hankel!(Hm, he)

    for i=1:length(Hm)
        @inbounds Tm_plus_Hm[i] = Hm[i] + Tm[i]
    end

    # Solve Ax = b
    # NOTE: both Tm_plus_Hm and b are overwritten
    copy!(b, 1, rr, 2, order)
    A_ldiv_B!(lufact!(Tm_plus_Hm), b)

    # Add the solution vector
    for i=2:order+1
        @inbounds c[i] += b[i-1]
    end

    if γ == -one(T)
        ϵ = gain(rr, c, order, γ)
        c[1] = √ϵ
    end

    log(ϵ)
end

# mgcepnorm! changes form of mel-generalized cepstrums
# input:
# K, bₐ,ᵧ'(1), ..., bₐ,ᵧ'(m)
#
# output(otype):
# 0: cₐ,ᵧ(0), cₐ,ᵧ(1), ..., cₐ,ᵧ(m)
# 1: bₐ,ᵧ(0), bₐ,ᵧ(1), ..., bₐ,ᵧ(m)
# 2: Kₐ, cₐ,ᵧ'(1), ..., cₐ,ᵧ'(m)
# 3: K, bₐ,ᵧ'(1), ..., bₐ,ᵧ'(m)
# 4: Kₐ, γcₐ,ᵧ'(1), ..., γcₐ,ᵧ'(m)
# 5: K, γbₐ,ᵧ'(1), ..., γbₐ,ᵧ'(m)
#
# For simplicity, we represent bₐ,ᵧ as bᵧ.
function mgcepnorm!{T<:AbstractFloat}(bᵧ′::AbstractVector{T},
                                      α::AbstractFloat,
                                      γ::AbstractFloat,
                                      otype::Int)
    0<=otype<=5 || throw(ArgumentError("0 ≤ otype ≤ 5 are supported"))

    mgc = bᵧ′

    if otype == 0 || otype == 1 || otype == 2 || otype == 4
        # K, bᵧ' -> bᵧ
        ignorm!(mgc, γ)
    end

    if otype == 0 || otype == 2 || otype == 4
        # bᵞ -> cᵞ
        b2mc!(mgc, α)
    end

    if otype == 2 || otype == 4
        # cᵧ -> cᵧ'
        gnorm!(mgc, γ)
    end

    if otype == 4 || otype == 5
        # cᵧ' -> γcᵧ'
        for i=2:length(mgc)
            @inbounds mgc[i] *= γ
        end
    end

    mgc
end

function _mgcep(x::AbstractVector,          # a *windowed* signal
                order=25,                 # order of mgcep
                α=0.35,         # all-pass constant
                γ=0.0;          # parameter of generalized log
                n::Int=length(x)-1,            # order of recursion
                miniter::Int=2,
                maxiter::Int=30,
                criteria::AbstractFloat=0.001, # stopping criteria
                e::Real=zero(eltype(x)),       # floor of
                otype::Int=0,                  # output type
                verbose::Bool=false
    )
    T = eltype(x)
    @assert n < length(x)

    # Periodogram
    periodogram = abs2.(fft(x)) + e

    # Allocate memory
    cr = zeros(T, length(x))
    pr = zeros(T, length(x))
    rr = zeros(T, length(x))
    ri = zeros(T, length(x))
    qr = zeros(T, length(x))
    qi = zeros(T, length(x))
    Tm = Array{T,2}(order, order)
    Hm = Array{T,2}(order, order)
    Tm_plus_Hm = Array{T,2}(order, order)
    b = Array{T,1}(order)

    # FFT workspace
    y = Array{Complex{T},1}(length(x))
    z = Array{Complex{T},1}(length(x))
    # bplan = FFTW.Plan(y, z, 1, FFTW.ESTIMATE, FFTW.NO_TIMELIMIT)
    bplan = plan_bfft(y)

    bᵧ′ = zeros(T, order+1)
    ϵ⁰ = newton!(bᵧ′, periodogram, order, α, -one(T), n, 1, y, z, bplan,
                 cr, pr, rr, ri, qr, qi, Tm, Hm, Tm_plus_Hm, b)

    if γ != -one(T)
        d = Array{T,1}(order+1)
        if α != zero(T)
            ignorm!(bᵧ′, -1.0)
            b2mc!(bᵧ′, α)
            copy!(d, bᵧ′)
            gnorm!(d, -1.0)
        else
            copy!(d, bᵧ′)
        end
        bᵧ′ = gc2gc(d, -1.0, order, γ)

        if α != zero(T)
            ignorm!(bᵧ′, γ)
            mc2b!(bᵧ′, α)
            gnorm!(bᵧ′, γ)
        end
    end

    if γ != -one(T)
        ϵᵗ = ϵ⁰
        for i=1:maxiter
            ϵ = newton!(bᵧ′, periodogram, order, α, γ, n, i, y, z, bplan,
                        cr, pr, rr, ri, qr, qi, Tm, Hm, Tm_plus_Hm, b)
            if i >= miniter
                err = abs((ϵᵗ - ϵ)/ϵ)
                verbose && println("nmse: $err")
                if err < criteria
                    break
                end
            end
            ϵᵗ = ϵ
        end
    end

    mgcepnorm!(bᵧ′, α, γ, otype)
end

function estimate(mgc::MelGeneralizedCepstrum, x::AbstractArray;
                  otype::Int=0,
                  use_sptk::Bool=false,
                  kargs...)
    if otype != 0
        throw(ArgumentError("""otype = 0 is only allowed.
                            Use `_mgcep` function if you need diffrent output types that are used in SPTK.mgcep."""))
    end
    order = param_order(mgc)
    α = allpass_alpha(mgc)
    γ = glog_gamma(mgc)
    normalized = γ == zero(γ) ? true : false
    mgcepfunc = use_sptk ? SPTK.mgcep : _mgcep
    data = mgcepfunc(x, order, α, γ; otype=otype, kargs...)
    SpectralParamState(mgc, data, true, normalized)
end

function mgcep(x::AbstractArray, order=25, α=0.35, γ=0.0; kargs...)
    estimate(MelGeneralizedCepstrum(order, α, γ), x; kargs...)
end


text:
module TestPackage

using Logging
using HDF5

# package code goes here

"""
If running non-interactively, call `exit(expr)`, reporting any exceptions

Use this for deploying command line tools written in julia.
"""
macro run_if_deployed(expr)
    quote
        # Hack to force serialization of package precompilation during package
        # install on 0.4 prerelease.  See
        # https://github.com/JuliaLang/julia/issues/12508
        do_precompile = length(ARGS) > 0 && ARGS[1] == "__precompile__"
        if !do_precompile
            try
                if !isinteractive()
                    blas_set_num_threads(1)
                    exit($(esc(expr)))
                end
            catch err
                bt = catch_backtrace()
                io = IOBuffer()
                showerror(io, err, bt)
                println(STDERR, "ERROR: Caught $(typeof(err)):\n$(takebuf_string(io))")
                exit(1)
            end
        end
    end
end

end # module