text:
abstract type AbstractNDSparse end

mutable struct NDSparse{T, D<:Tuple, C<:Columns, V<:AbstractVector} <: AbstractNDSparse
    index::C
    data::V
    _table::IndexedTable

    index_buffer::C
    data_buffer::V
end

function IndexedTable(nds::NDSparse; kwargs...)
    convert(IndexedTable, nds.index, nds.data; kwargs...)
end

convert(::Type{IndexedTable}, nd::NDSparse) = IndexedTable(nd)

"""
    ndsparse(keys, values; kw...)

Construct an NDSparse array with the given `keys` and `values` columns. On construction,
the keys and data are sorted in lexicographic order of the `keys`.

# Keyword Argument Options:

- `agg = nothing` -- Function to aggregate values with duplicate keys.
- `presorted = false` -- Are the key columns already sorted?
- `copy = true` -- Should the columns in `keys` and `values` be copied?
- `chunks = nothing` -- Provide an integer to distribute data into `chunks` chunks.
    - A good choice is `nworkers()` (after `using Distributed`)
    - See also: [`distribute`](@ref)

# Examples:

    x = ndsparse(["a","b"], [3,4])
    keys(x)
    values(x)
    x["a"]

    # Dimensions are named if constructed with a named tuple of columns
    x = ndsparse((index = 1:10,), rand(10))
    x[1]

    # Multiple dimensions by passing a (named) tuple of columns
    x = ndsparse((x = 1:10, y = 1:2:20), rand(10))
    x[1, 1]

    # Value columns can also have names via named tuples
    x = ndsparse(1:10, (x=rand(10), y=rand(10)))
"""
function ndsparse end

function ndsparse(I::Tup, d::Union{Tup, AbstractVector};
                  chunks=nothing, kwargs...)
    if chunks !== nothing
        impl = Val{:distributed}()
    else
        impl = _impl(astuple(I)...)
        if impl === Val{:serial}()
            impl = isa(d, Tup) ?
                _impl(impl, astuple(d)...) : _impl(d)
        end
    end
    ndsparse(impl, I, d; chunks=chunks, kwargs...)
end

function ndsparse(::Val{:serial}, ks::Tup, vs::Union{Tup, AbstractVector};
                  agg=nothing, presorted=false,
                  chunks=nothing, copy=true)

    I = rows(ks)
    d = vs isa Tup ? Columns(vs) : vs

   #if !isempty(filter(x->!isa(x, Int),
   #                   intersect(colnames(I), colnames(d))))
   #    error("All column names, including index and data columns, must be distinct")
   #end
    if !isconcretetype(eltype(d)) || fieldcount(eltype(d)) !== 0
        length(I) == length(d) ||
            error("index and data must have the same number of elements")
    end

    if !presorted && !issorted(I)
        p = sortperm(I)
        I = I[p]
        d = d[p]
    elseif copy
        if agg !== nothing
            iter = igroupreduce(agg, I, d, Base.OneTo(length(I)))
            I, d = collect_columns(iter) |> columns
            agg = nothing
        else
            I = copyto!(similar(I), I)
            d = copyto!(similar(d), d)
        end
    end
    stripnames(x) = isa(x, Columns) ? rows(astuple(columns(x))) : rows((x,))
    _table = convert(IndexedTable, I, stripnames(d); presorted=true, copy=false)
    nd = NDSparse{eltype(d),astuple(eltype(I)),typeof(I),typeof(d)}(
        I, d, _table, similar(I,0), similar(d,0)
    )
    agg===nothing || aggregate!(agg, nd)
    return nd
end

function ndsparse(x::AbstractVector, y; kwargs...)
    ndsparse((x,), y; kwargs...)
end

function ndsparse(x::Tup, y::Columns; kwargs...)
    ndsparse(x, columns(y); kwargs...)
end

function ndsparse(x::Columns, y::AbstractVector; kwargs...)
    ndsparse(columns(x), y; kwargs...)
end

ndsparse(c::Columns{<:Pair}; kwargs...) =
    convert(NDSparse, columns(c).first, columns(c).second; kwargs...)

# backwards compat
NDSparse(idx::Columns, data; kwargs...) = ndsparse(idx, data; kwargs...)

# Easy constructor to create a derivative table
function ndsparse(t::NDSparse; presorted=true, copy=false)
    ndsparse(keys(t), values(t), presorted=presorted, copy=copy)
end

# TableLike API
Base.@pure function colnames(t::NDSparse)
    dnames = colnames(t.data)
    if all(x->isa(x, Integer), dnames)
        dnames = map(x->x+ncols(t.index), dnames)
    end
    (colnames(t.index)..., dnames...,)
end

columns(nd::NDSparse) = concat_tup(columns(nd.index), columns(nd.data))

# IndexedTableLike API

permcache(t::NDSparse) = permcache(t._table)
cacheperm!(t::NDSparse, p) = cacheperm!(t._table, p)

"""
    pkeynames(t::NDSparse)

Names of the primary key columns in `t`.

# Example

    x = ndsparse([1,2],[3,4])
    pkeynames(x)

    x = ndsparse((x=1:10, y=1:2:20), rand(10))
    pkeynames(x)
"""
pkeynames(t::NDSparse) = (dimlabels(t)...,)

# For an NDSparse, valuenames is either a tuple of fieldnames or a
# single name for scalar values
function valuenames(t::NDSparse)
    if isa(values(t), Columns)
        T = eltype(values(t))
        ((ndims(t) .+ (1:fieldcount(eltype(values(t)))))...,)
    else
        ndims(t) + 1
    end
end

"""
`NDSparse(columns...; names=Symbol[...], kwargs...)`

Construct an NDSparse array from columns. The last argument is the data column, and the rest are index columns. The `names` keyword argument optionally specifies names for the index columns (dimensions).
"""
function NDSparse(columns...; names=nothing, rest...)
    keys, data = columns[1:end-1], columns[end]
    ndsparse(Columns(keys, names=names), data; rest...)
end

similar(t::NDSparse) = NDSparse(similar(t.index, 0), similar(t.data, 0))

function copy(t::NDSparse)
    flush!(t)
    NDSparse(copy(t.index), copy(t.data), presorted=true)
end

function (==)(a::NDSparse, b::NDSparse)
    flush!(a); flush!(b)
    return a.index == b.index && a.data == b.data
end

function Base.isequal(a::NDSparse, b::NDSparse)
    flush!(a); flush!(b)
    return isequal(keys(a), keys(b)) && isequal(values(a), values(b))
end

function empty!(t::NDSparse)
    empty!(t.index)
    empty!(t.data)
    empty!(t.index_buffer)
    empty!(t.data_buffer)
    return t
end

_convert(::Type{<:Tuple}, tup::Tuple) = tup
_convert(::Type{T}, tup::Tuple) where {T<:NamedTuple} = T(tup)
convertkey(t::NDSparse{V,K,I}, tup::Tuple) where {V,K,I} = _convert(eltype(I), tup)

ndims(t::NDSparse) = length(fieldarrays(t.index))
length(t::NDSparse) = (flush!(t);length(t.index))
eltype(::Type{NDSparse{T,D,C,V}}) where {T,D,C,V} = T
Base.keytype(::Type{NDSparse{T,D,C,V}}) where {T,D,C,V} = D
Base.keytype(x::NDSparse) = keytype(typeof(x))
dimlabels(::Type{NDSparse{T,D,C,V}}) where {T,D,C,V} = fieldnames(eltype(C))

# Generic ndsparse constructor that also works with distributed
# arrays in JuliaDB

# Keys and Values iterators

keys(t::NDSparse) = t.index
"""
`keys(x::NDSparse[, select::Selection])`

Get the keys of an `NDSparse` object. Same as [`rows`](@ref) but acts only on the index columns of the `NDSparse`.
"""
keys(t::NDSparse, which...) = rows(keys(t), which...)

# works for both IndexedTable and NDSparse
pkeys(t::NDSparse, which...) = keys(t, which...)

values(t::NDSparse) = t.data
"""
`values(x::NDSparse[, select::Selection])`

Get the values of an `NDSparse` object. Same as [`rows`](@ref) but acts only on the value columns of the `NDSparse`.
"""
function values(t::NDSparse, which...)
    if values(t) isa Columns
        rows(values(t), which...)
    else
        if which[1] != 1
            error("column $which not found")
        end
        values(t)
    end
end

## Some array-like API

"""
`dimlabels(t::NDSparse)`

Returns an array of integers or symbols giving the labels for the dimensions of `t`.
`ndims(t) == length(dimlabels(t))`.
"""
dimlabels(t::NDSparse) = dimlabels(typeof(t))

iterate(a::NDSparse) = iterate(a.data)
iterate(a::NDSparse, st) = iterate(a.data, st)

function permutedims(t::NDSparse, p::AbstractVector)
    if !(length(p) == ndims(t) && isperm(p))
        throw(ArgumentError("argument to permutedims must be a valid permutation"))
    end
    flush!(t)
    NDSparse(Columns(columns(t.index)[p]), t.data, copy=true)
end

# showing
function show(io::IO, t::NDSparse{T,D}) where {T,D}
    flush!(t)
    if !(values(t) isa Columns)
        cnames = colnames(keys(t))
        eltypeheader = "$(eltype(t))"
    else
        cnames = colnames(t)
        nf = fieldcount(eltype(t))
        if eltype(t) <: NamedTuple
            eltypeheader = "$(nf) field named tuples"
        else
            eltypeheader = "$(nf)-tuples"
        end
    end
    header = "$(ndims(t))-d NDSparse with $(length(t)) values (" * eltypeheader * "):"
    showtable(io, t; header=header,
              cnames=cnames, divider=length(columns(keys(t))))
end

function showmeta(io, t::NDSparse, cnames)
    nc = length(columns(t))
    nidx = length(columns(keys(t)))
    nkeys = length(columns(values(t)))

    print(io,"    ")
    printstyled(io, "Dimensions", color=:underline)
    metat = Columns(([1:nidx;], [Text(get(cnames, i, "<noname>")) for i in 1:nidx],
                     eltype.([columns(keys(t))...])))
    showtable(io, metat, cnames=["#", "colname", "type"], cstyle=fill(:bold, nc), full=true)
    print(io,"\n    ")
    printstyled(io, "Values", color=:underline)
    if isa(values(t), Columns)
        metat = Columns(([nidx+1:nkeys+nidx;], [Text(get(cnames, i, "<noname>")) for i in nidx+1:nkeys+nidx],
                         eltype.(Any[columns(values(t))...])))
        showtable(io, metat, cnames=["#", "colname", "type"], cstyle=fill(:bold, nc), full=true)
    else
        show(io, eltype(values(t)))
    end
end

@noinline convert(::Type{NDSparse}, @nospecialize(ks), @nospecialize(vs); kwargs...) = ndsparse(ks, vs; kwargs...)
@noinline convert(T::Type{NDSparse}, c::Columns{<:Pair}; kwargs...) = convert(T, columns(c).first, columns(c).second; kwargs...)

# map and convert

"""
    map(f, x::NDSparse; select = values(x))

Apply `f` to every value of `select` selected from `x` (see [`select`](@ref)).

Apply `f` to every data value in `x`. `select` selects fields
passed to `f`. By default, the data values are selected.

If the return value of `f` is a tuple or named tuple the result will contain many data columns.

# Examples

    x = ndsparse((t=[0.01, 0.05],), (x=[1,2], y=[3,4]))

    polar = map(row -> (r = hypot(row.x, row.y), θ = atan(row.y, row.x)), x)

    back2x = map(row -> (x = row.r * cos(row.θ), y = row.r * sin(row.θ)), polar)
"""
function map(f, x::NDSparse; select=x.data)
    ndsparse(copy(x.index), map_rows(f, rows(x, select)), presorted=true, copy=false)
end

# """
# `columns(x::NDSparse, names...)`
#
# Given an NDSparse array with multiple data columns (its data vector is a `Columns` object), return a
# new array with the specified subset of data columns. Data is shared with the original array.
# """
# columns(x::NDSparse, which...) = NDSparse(x.index, Columns(columns(x.data)[[which...]]), presorted=true)

#columns(x::NDSparse, which) = NDSparse(x.index, columns(x.data)[which], presorted=true)

#column(x::NDSparse, which) = columns(x, which)

# NDSparse uses lex order, Base arrays use colex order, so we need to
# reorder the data. transpose and permutedims are used for this.
function convert(::Type{NDSparse}, m::SparseMatrixCSC)
    A = transpose(m)
    nzidx = findall(!iszero, A)
    I,J,V = getindex.(nzidx, 1), getindex.(nzidx, 2), A[nzidx]
    NDSparse(J, I, V, presorted=true)
end

# special method to allow selection on
# ndsparse with repeating names in keys and values
function column(x::NDSparse, which::Integer)
    @assert which > 0
    if which <= ndims(x)
        keys(x, which)
    else
        values(x, which-ndims(x))
    end
end

function convert(::Type{NDSparse}, a::AbstractArray{T}) where T
    n = length(a)
    nd = ndims(a)
    a = permutedims(a, [nd:-1:1;])
    data = reshape(a, (n,))
    idxs = [ Vector{Int}(undef, n) for i = 1:nd ]
    i = 1
    for I in CartesianIndices(size(a))
        for j = 1:nd
            idxs[j][i] = I[j]
        end
        i += 1
    end
    NDSparse(Columns(Tuple(Iterators.reverse(idxs))), data, presorted=true)
end

# aggregation

"""
    aggregate!(f::Function, arr::NDSparse)

Combine adjacent rows with equal indices using the given 2-argument reduction function,
in place.
"""
function aggregate!(f, x::NDSparse)
    idxs, data = x.index, x.data
    n = length(idxs)
    newlen = 0
    for ii in GroupPerm(compact_mem(idxs), Base.OneTo(n))
        newlen += 1
        if newlen != last(ii)
            copyrow!(idxs, newlen, last(ii))
            @inbounds data[newlen] = reduce(f, data[i] for i in ii)
        end
    end
    resize!(idxs, newlen)
    resize!(data, newlen)
    x
end

function subtable(x::NDSparse, idx; presorted=true)
    ndsparse(keys(x)[idx], values(x)[idx], presorted=presorted, copy=false)
end


text:
# read in data from NFL (cut and paste from table produced tab delim data)
#      data from https://www.pro-football-reference.com/years/2009/games.htm
# and output in my format
#
using DataFrames
using CSV
using RatPack

# read in raw TSV file
file = "nfl_2009.txt"
year = 2009 # to 2010
inp_df = CSV.read(file; comment="#", delim='\t', normalizenames=true, )
n = size(inp_df,1)
names(inp_df)
win_tie = sign.( inp_df[:PtsW] .- inp_df[:PtsL] )# Outcome is always 1 or 0, because they order the data Winner,Looser

# output results in my format
out_df = DataFrame(PlayerA=inp_df[:Winner_tie],
                   PlayerB=inp_df[:Loser_tie],
                   Outcome=win_tie,
                   FactorA=zeros(Int,n),
                   FactorB=zeros(Int,n),
                   ScoreA=inp_df[:PtsW],
                   ScoreB=inp_df[:PtsL])
colnames = Symbol.(["Player A", "Player B", "Outcome", "Factor A", "Factor B", "Score A", "Score B"])
names!(out_df, colnames)
out_file = replace(file, ".txt" => ".csv")
CSV.write(out_file, out_df)

# quick test to see if file can be read in OK
(df, player_list) = read_results( out_file )
input_ratings = RatingsList( player_list )
input_competitions = df
ratings = update_ratings(UpdateColley(),
                         input_ratings,
                         input_competitions)
sort(collect(ratings.ratings), by = tuple -> last(tuple), rev=true)

# now do likewise for only the regular season
week = inp_df[:Week]
df = inp_df[ .!(tryparse.(Int, coalesce.(week,"")) .== nothing) , :]
n = size(df,1)
win_tie = sign.( df[:PtsW] .- df[:PtsL] )# Outcome is always 1 or 0, because they order the data Winner,Looser
out_df = DataFrame(PlayerA=df[:Winner_tie],
                   PlayerB=df[:Loser_tie],
                   Outcome=win_tie,
                   FactorA=zeros(Int,n),
                   FactorB=zeros(Int,n),
                   ScoreA=df[:PtsW],
                   ScoreB=df[:PtsL])
colnames = Symbol.(["Player A", "Player B", "Outcome", "Factor A", "Factor B", "Score A", "Score B"])
names!(out_df, colnames)
out_file = replace(file, ".txt" => "_regular.csv")
CSV.write(out_file, out_df)

# now output one file per week (17 weeks in regular season)
weeks = unique( coalesce.(week,"") )
for w in weeks
    if !(w == nothing)
        df = inp_df[ week .== w, :]
        n = size(df,1)
        win_tie = sign.( df[:PtsW] .- df[:PtsL] )# Outcome is always 1 or 0, because they order the data Winner,Looser
        out_df = DataFrame(PlayerA=df[:Winner_tie],
                           PlayerB=df[:Loser_tie],
                           Outcome=win_tie,
                           FactorA=zeros(Int,n),
                           FactorB=zeros(Int,n),
                           ScoreA=df[:PtsW],
                           ScoreB=df[:PtsL])
        colnames = Symbol.(["Player A", "Player B", "Outcome", "Factor A", "Factor B", "Score A", "Score B"])
        names!(out_df, colnames)
        out_file = replace(file, ".txt" => "_week_$(w).csv")
        CSV.write(out_file, out_df)
    end
end

ext: jl


text:
using Distributed
using JLD2
using RCall

@time @everywhere include("src/required.jl")

df = CSV.read("../data/ebd/processed/ebd_warblers_cut.csv", header=true, delim="\t")
head = CSV.read("../data/ebd/ebd_warblers_head.csv", header=true, delim="\t")
prep = CSV.read("../data/ebd/ebd_warblers_prep.csv", header=true, delim="\t")
@load "../data/raw-Y-matrices.jld2" Y Ypred Ytransf inds_pred inds_notpred
@load "../data/pres-abs-ebd.jld2" pres_abs
wc_vars = pmap(worldclim, 1:19)

# Prepare data (arrange values & columns)
df = prepare_ebd_data(df)
head = prepare_ebd_data(head)

# Remove duplicates (BUT not ok for counts, so-so for dates, see group-observation.jl in src/test/)
df_nogroups = filter(x -> ismissing(x[:groupIdentifier]), df)
df_groups = dropmissing(df, :groupIdentifier)
df_groups_unique = unique(df_groups, [:species, :groupIdentifier])
newdf = vcat(df_nogroups, df_groups_unique)

head_nogroups = filter(x -> ismissing(x[:groupIdentifier]), head)
head_groups = dropmissing(head, :groupIdentifier)
head_groups_unique = unique(head_groups, [:species, :groupIdentifier])
newhead = vcat(head_nogroups, head_groups_unique)

$library(xtable)

# Variables base
varbase = select(newhead, [:species, :latitude, :longitude, :observationDate, :allSpeciesReported, :observationCount])
@rput varbase
$varbase$observationDate = as.character(varbase$observationDate)
$xtable(head(varbase), format="latex")

# Variable échant
varechant = select(newhead, [:protocolType, :durationMinutes, :effortDistanceKm, :numberObservers])
@rput varechant
$xtable(head(varechant), format="latex")

speciescounts = by(newdf, :species, nrow)
protocolcounts = by(newdf, :protocolType, nrow)
show(sort(protocolcounts, :x1, rev=true), allrows=true)

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
using HDF5,JLD,KUparser,KUnet,Base.Test
include("pxyequal.jl")

info("Load data")
@date @load "conll07.tst.jld4"
ncpu = 10
nbatch = 10
ndeps = length(deprel)
ft = Flist.acl11eager
pt = ArcEager13
p0 = pt(1,ndeps)
s1 = corpus[1]

info("Let's train a reasonable model")
# net = newnet(relu, flen(p0,s1,ft), 20000, p0.nmove; learningRate=1)
# net[end].f = logp
h0 = 20000
net = [Mmul(h0), Bias(), Relu(),
       Mmul(p0.nmove), Bias(), Logp(),
       LogpLoss()]
@date (p,x,y)=oparse(pt, corpus, ndeps, ft)
@date train(net, hcat(x...), hcat(y...))

info("Single sentence")
@date p1=gparse(pt, [s1], ndeps, ft, net)
@show evalparse(p1, [s1])
@date (p2,x2,y2)=gparse(pt, [s1], ndeps, ft, net; returnxy=true)
@show map(size, (p2,x2,y2))
@show I2 = 1:size(x2[1],2)
@test @show isequal(p2,p1)

info("Multiple sentences")
@date p3=gparse(pt, corpus, ndeps, ft, net)
@show evalparse(p3, corpus)
@test @show isequal(p3[1],p1[1])
@date pxy4=(p4,x4,y4)=gparse(pt, corpus, ndeps, ft, net; returnxy=true)
@show map(size, pxy4)
@test @show isequal(p4, p3)
@test @show isequal(x4[1][:,I2], x2[1])
@test @show isequal(y4[1][:,I2], y2[1])

info("Batch processing")
@date p5=gparse(pt, corpus, ndeps, ft, net; nbatch=nbatch)
@test @show isequal(p5,p3)
@date pxy6=gparse(pt, corpus, ndeps, ft, net; nbatch=nbatch, returnxy=true)
@show map(size, pxy6)
@test @show pxyequal(pxy6, pxy4)

info("Multi-cpu processing")
@date p7=gparse(pt, corpus, ndeps, ft, net; nbatch=nbatch, usepmap=true)
@test @show isequal(p7,pxy6[1])
@date pxy8=gparse(pt, corpus, ndeps, ft, net; nbatch=nbatch, returnxy=true, usepmap=true)
@test @show pxyequal(pxy8, pxy6)


text:
int_rules_4_1_13 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.1.13*(d+e*x)^m*sin(a+b*x+c*x^2)^n =#
    @apply_utils Antiderivative(sin(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => Antiderivative(sin((~b + 2 * ~c * ~x) ^ 2 / (4 * ~c)), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && EqQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative(cos(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => Antiderivative(cos((~b + 2 * ~c * ~x) ^ 2 / (4 * ~c)), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && EqQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative(sin(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => cos(((~b) ^ 2 - 4 * ~a * ~c) / (4 * ~c)) * Antiderivative(sin((~b + 2 * ~c * ~x) ^ 2 / (4 * ~c)), ~x) - sin(((~b) ^ 2 - 4 * ~a * ~c) / (4 * ~c)) * Antiderivative(cos((~b + 2 * ~c * ~x) ^ 2 / (4 * ~c)), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative(cos(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => cos(((~b) ^ 2 - 4 * ~a * ~c) / (4 * ~c)) * Antiderivative(cos((~b + 2 * ~c * ~x) ^ 2 / (4 * ~c)), ~x) + sin(((~b) ^ 2 - 4 * ~a * ~c) / (4 * ~c)) * Antiderivative(sin((~b + 2 * ~c * ~x) ^ 2 / (4 * ~c)), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative(sin(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~n, ~x) => Antiderivative(ExpandTrigReduce(sin(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && IGtQ(~n, 1)
    @apply_utils Antiderivative(cos(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~n, ~x) => Antiderivative(ExpandTrigReduce(cos(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && IGtQ(~n, 1)
    @apply_utils Antiderivative(sin(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(n'), ~x) => Unintegrable(sin(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
    @apply_utils Antiderivative(cos(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(n'), ~x) => Unintegrable(cos(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
    @apply_utils Antiderivative(sin(~v) ^ ~(n'), ~x) => Antiderivative(sin(ExpandToSum(~v, ~x)) ^ ~n, ~x) <-- IGtQ(~n, 0) && (QuadraticQ(~v, ~x) && Not(QuadraticMatchQ(~v, ~x)))
    @apply_utils Antiderivative(cos(~v) ^ ~(n'), ~x) => Antiderivative(cos(ExpandToSum(~v, ~x)) ^ ~n, ~x) <-- IGtQ(~n, 0) && (QuadraticQ(~v, ~x) && Not(QuadraticMatchQ(~v, ~x)))
    @apply_utils Antiderivative((~d + ~(e') * ~x) * sin(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (-(~e) * cos(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && EqQ(2 * ~c * ~d - ~b * ~e, 0)
    @apply_utils Antiderivative((~d + ~(e') * ~x) * cos(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (~e * sin(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && EqQ(2 * ~c * ~d - ~b * ~e, 0)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * sin(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (-(~e) * (~d + ~e * ~x) ^ (~m - 1) * cos(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) + (((~e) ^ 2 * (~m - 1)) / (2 * ~c)) * Antiderivative((~d + ~e * ~x) ^ (~m - 2) * cos(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(2 * ~c * ~d - ~b * ~e, 0) && GtQ(~m, 1))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * cos(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (~e * (~d + ~e * ~x) ^ (~m - 1) * sin(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) - (((~e) ^ 2 * (~m - 1)) / (2 * ~c)) * Antiderivative((~d + ~e * ~x) ^ (~m - 2) * sin(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(2 * ~c * ~d - ~b * ~e, 0) && GtQ(~m, 1))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * sin(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * sin(~a + ~b * ~x + ~c * (~x) ^ 2)) / (~e * (~m + 1)) - ((2 * ~c) / ((~e) ^ 2 * (~m + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m + 2) * cos(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(2 * ~c * ~d - ~b * ~e, 0) && LtQ(~m, -1))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * cos(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * cos(~a + ~b * ~x + ~c * (~x) ^ 2)) / (~e * (~m + 1)) + ((2 * ~c) / ((~e) ^ 2 * (~m + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m + 2) * sin(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(2 * ~c * ~d - ~b * ~e, 0) && LtQ(~m, -1))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) * sin(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (-(~e) * cos(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) + ((2 * ~c * ~d - ~b * ~e) / (2 * ~c)) * Antiderivative(sin(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && NeQ(2 * ~c * ~d - ~b * ~e, 0)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) * cos(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (~e * sin(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) + ((2 * ~c * ~d - ~b * ~e) / (2 * ~c)) * Antiderivative(cos(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && NeQ(2 * ~c * ~d - ~b * ~e, 0)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * sin(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => ((-(~e) * (~d + ~e * ~x) ^ (~m - 1) * cos(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) - ((~b * ~e - 2 * ~c * ~d) / (2 * ~c)) * Antiderivative((~d + ~e * ~x) ^ (~m - 1) * sin(~a + ~b * ~x + ~c * (~x) ^ 2), ~x)) + (((~e) ^ 2 * (~m - 1)) / (2 * ~c)) * Antiderivative((~d + ~e * ~x) ^ (~m - 2) * cos(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ(~b * ~e - 2 * ~c * ~d, 0) && GtQ(~m, 1))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * cos(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => ((~e * (~d + ~e * ~x) ^ (~m - 1) * sin(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) - ((~b * ~e - 2 * ~c * ~d) / (2 * ~c)) * Antiderivative((~d + ~e * ~x) ^ (~m - 1) * cos(~a + ~b * ~x + ~c * (~x) ^ 2), ~x)) - (((~e) ^ 2 * (~m - 1)) / (2 * ~c)) * Antiderivative((~d + ~e * ~x) ^ (~m - 2) * sin(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ(~b * ~e - 2 * ~c * ~d, 0) && GtQ(~m, 1))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * sin(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (((~d + ~e * ~x) ^ (~m + 1) * sin(~a + ~b * ~x + ~c * (~x) ^ 2)) / (~e * (~m + 1)) - ((~b * ~e - 2 * ~c * ~d) / ((~e) ^ 2 * (~m + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m + 1) * cos(~a + ~b * ~x + ~c * (~x) ^ 2), ~x)) - ((2 * ~c) / ((~e) ^ 2 * (~m + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m + 2) * cos(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ(~b * ~e - 2 * ~c * ~d, 0) && LtQ(~m, -1))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * cos(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * cos(~a + ~b * ~x + ~c * (~x) ^ 2)) / (~e * (~m + 1)) + ((~b * ~e - 2 * ~c * ~d) / ((~e) ^ 2 * (~m + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m + 1) * sin(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) + ((2 * ~c) / ((~e) ^ 2 * (~m + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m + 2) * sin(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ(~b * ~e - 2 * ~c * ~d, 0) && LtQ(~m, -1))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * sin(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~n, ~x) => Antiderivative(ExpandTrigReduce((~d + ~e * ~x) ^ ~m, sin(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && IGtQ(~n, 1)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * cos(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~n, ~x) => Antiderivative(ExpandTrigReduce((~d + ~e * ~x) ^ ~m, cos(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && IGtQ(~n, 1)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * sin(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(n'), ~x) => Unintegrable((~d + ~e * ~x) ^ ~m * sin(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * cos(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(n'), ~x) => Unintegrable((~d + ~e * ~x) ^ ~m * cos(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative((~u) ^ ~(m') * sin(~v) ^ ~(n'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * sin(ExpandToSum(~v, ~x)) ^ ~n, ~x) <-- FreeQ(~m, ~x) && (IGtQ(~n, 0) && (LinearQ(~u, ~x) && (QuadraticQ(~v, ~x) && Not(LinearMatchQ(~u, ~x) && QuadraticMatchQ(~v, ~x)))))
    @apply_utils Antiderivative((~u) ^ ~(m') * cos(~v) ^ ~(n'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * cos(ExpandToSum(~v, ~x)) ^ ~n, ~x) <-- FreeQ(~m, ~x) && (IGtQ(~n, 0) && (LinearQ(~u, ~x) && (QuadraticQ(~v, ~x) && Not(LinearMatchQ(~u, ~x) && QuadraticMatchQ(~v, ~x)))))
end

max_stars_repo_path: src/rules/4 Trig functions/4.1 Sine/4.1.13 (d+e x)^m sin(a+b x+c x^2)^n.jl


text:
using Bridge, Distributions
import Bridge: lp

PLOT = false
include("plot.jl")

struct Atan <: ContinuousTimeProcess{Float64}
    α::Float64
    β::Float64
    σ::Float64
    a::Float64
    Atan(a,b,c) = new(a,b,c,c^2)
end

function param(θ, σ)
        θ[1], θ[2], σ[1]
end

Bridge.b(t, x, P::Atan) =  P.α*atan(x) + P.β
phi12(t, x, P::Atan) = [atan(x), 1.]
intercept12(t, x, P::Atan) = 0.

Bridge.σ(t, x, P::Atan) = P.σ
Bridge.a(t, x, P::Atan) = P.a
Bridge.Γ(t, x, P::Atan) = inv(P.a)

function conjugateb(YY, th, xi, P, phif, intc)
        n = length(th)
        G = zeros(n,n)
        mu = zeros(th)
        for Y in YY
            for i in 1:length(Y)-1
                phi = phif(Y[i]..., P)
                Gphi = Bridge.Γ(Y[i]..., P)*phi
                zi = phi*Gphi'
                dy = Y.yy[i+1] - Y.yy[i]
                ds = Y.tt[i+1] - Y.tt[i]
                mu = mu + Gphi*(dy - intc(Y[i]..., P)*ds)
                G[:] = G +  zi * (Y.tt[i+1]-Y.tt[i])
            end
        end #for m
        WW = G .+ diagm(xi)
        WL = transpose(chol(WW))
        th° = WL'\(randn(n)+WL\mu)
end

mcstart(yy) = (zeros(yy), zeros(yy), 0)
function mcnext(mc, x)
    m, m2, n = mc
    delta = x - m
    n = n + 1
    m = m + delta*(1/n)
    m2 = m2 + map((x,y)->x.*y, delta, x - m)
    m, m2, n
end
function mcbandste(mc)
    m, m2, k = mc
    Q = sqrt(2.)*erfinv(0.95)

    ste = eltype(m)[sqrt(v) for v in m2 * (1/(k - 1))]*sqrt(1/k)
    m -Q*ste, m +Q*ste
end
function mcband(mc)
    m, m2, k = mc
    Q = sqrt(2.)*erfinv(0.95)

    std = eltype(m)[sqrt(v) for v in m2 * (1/(k - 1))]
    m-Q*std, m+Q*std
end

function MyProp(u, v, P, proptype=:mbb)
    if proptype == :guip
        #Pt = LinPro(P.α*cos(-P.β/P.α)^2, 0.5*P.α*sin(2P.β/P.α), P.σ)
        Pt = LinPro(P.α*cos(-P.β/P.α)^2, tan(-P.β/P.α), P.σ)
        return GuidedProp(P, u..., v..., Pt)
    else
        return DHBridgeProp(P, u..., v...)
    end
end

############## Configuration ###################################
Random.seed!(10)
K = 100
m = 10 # number of euler steps per segments

simid = 2
propid = 2
proptype = [:mbb,:guip][propid]

simname =["atan", "atantc"][simid]
# simname =["atan", "atantc"][simid] * "$proptype$m"

θ = [-2., 0.]
θtrue=copy(θ)
scaleθ = 0.08*[1., 1.]
σ = [0.75]
scaleσ = [0.1]
σtrue = copy(σ)
Ptrue = Atan(param(θ, σ)...)

n = 100 # number of segments
mextra = 2*div(1000,m) #factor of extra steps for truth
TT = 30.
dt = TT/n/m
tt = range(0., stop=TT, length=n*m+1)
tttrue = range(0., stop=TT, length=n*mextra*m+1)
ttf = tt[1:m:end]

uu = 0.

r = [(:xrange,(-2,2)), (:yrange,(-1,3))]

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

Y = euler(uu, sample(tttrue, Wiener{Float64}()), Ptrue)
Yobs = Y[1:mextra*m:end] #subsample
assert(endof(Yobs) == n+1)

Y2 = SamplePath(tt, copy(Y.yy[1:mextra:end])) #subsample
normalizetime(tt) = (tt - tt[1])/(tt[end]-tt[1])

θs = Float64[]
iter = 0
tts = [collect(tt[m*(i-1)+1:m*i+1]) for i in 1:n]
if simid == 2
    tts = [Bridge.tau(tts[i]) for i in 1:n]
end

BB = SamplePath[Y2[m*(i-1)+1:m*i+1] for i in 1:n]
for i in 1:n
    BB[i].tt[:] = tts[i]
    BB[i].yy[:] =  map(x -> BB[i].yy[1] + x*(BB[i].yy[end] - BB[i].yy[1]),normalize(tts[i]))
end

BBnew = copy(BB)
BBall = vcat([BB[i][1:end-1] for i in 1:n]...)

#
if PLOT
        xr = (5,7)
        plot(Yobs, "o"; xrange=xr, yrange=(-3,3),linewidth=0.1)
        display(oplot(BBall,"b"; xrange=xr, yrange=(-3,3), linewidth=0.7))
end

################### Prior ###################################
conjθs = [1,2] #set of conjugate thetas
phi = phi12
intercept = intercept12
estθ = [true, true] #params to estimate

# arbitrary starting values

σ = [2.]
θ = [-0.1, -0.1]

for i in 1:length(θ) # start with truth for parameters not to be estimated
    if !estθ[i]
        θ[i] = θtrue[i]
        scaleθ[i] = 0.
    end
end

P = Atan(param(θ, σ)...)

# reserve some space
i = 1
ww = Array{Float64,1}(length(m*(i-1)+1:m*(i)+1))
yy = copy(ww)

# Prior
Alpha = 1/500
Beta= 1/500

piσ²(s2) = pdf(InverseGamma(Alpha,Beta), s2[1])
lp(x, si) = Bridge.logpdfnormal(x, si^2)
PiError = InverseGamma(Alpha,Beta)

xi = 1./[5., 5.] #prior prec

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

# Bookkeeping
mkpath(joinpath("output",simname))
try # save cp of this file as documentation
    cp(@__FILE__(), joinpath("output",simname,"$simname.jl"); force=true)
catch
end

open(joinpath("output",simname,"truth.txt"), "w") do f
    println(f, "alpha beta sigma")
    println(f, join(round.([θtrue ; σtrue],3)," "))
end

open(joinpath("output",simname,"params.txt"), "w") do f
    println(f, "n alpha beta sigma")
end

# initialize
bacc = 0
siacc = 0

mc = mcstart(vcat([BB[i][1:end-1] for i in 1:n]...).yy )
mcparams = mcstart([θ ; σ])

Bridge.constdiff(::Atan) = true

perf = @timed while true
    P = Atan(param(θ, σ)...)
    iter += 1

    for i in 1:n-1
        P° = MyProp(Yobs[i], Yobs[i+1], P, proptype)
        B = bridge!(SamplePath(tts[i], yy), sample!(SamplePath(tts[i],ww), Wiener{Float64}()),P°)
        # B = shiftedeulerb!(SamplePath(tts[i], yy), sample!(SamplePath(tts[i],ww), Wiener{Float64}()),P°)
        if iter == 1
             BB[i] = B
        end
        llold = llikelihood(BB[i], P°)
        llnew = llikelihood(B, P°)

        if rand() < exp(llnew - llold)
            bacc += 1
            BB[i] = copy(B)
        end
    end

    # update theta

    if iter % 1 == 0
        #BBall = vcat([BB[i][1:end-1] for i in 1:n]...)
        θ° = θ + (2rand(length(θ)) .- 1).*scaleθ
        Pθ = Atan(param(θ, σ)...)
        Pθ° = Atan(param(θ°, σ)...)
        #ll = girsanov(BBall, Pθ°, Pθ)
        ll = 0.
        for i in 1:n
            ll += girsanov(BB[i], Pθ°, Pθ)
        end
        if rand() < exp(ll)
            θ = θ°
        end
    end

    # update conjugate theta

    θ[conjθs] = conjugateb(BB, θ[conjθs],xi[conjθs], P, phi, intercept)

    # update sigma (and theta)
    if iter % 1 == 0
        σ° = σ .* exp.(scaleσ .* randn(length(σ)))
        θ° = θ #+ (2rand(length(θ)) .- 1).*scaleθ/3
        Pσ = Atan(param(θ, σ)...)
        Pσ° = Atan(param(θ°, σ°)...)
        ll = 0.
        for i in 1:n # reparametrization
            P° = MyProp(Yobs[i], Yobs[i+1], Pσ, proptype)
            P°° = MyProp(Yobs[i], Yobs[i+1], Pσ°, proptype)
            Z = innovations(BB[i], P°)
            BBnew[i] = bridge(Z, P°°)
            # BBnew[i] = eulerb(Z, P°°)
            ll += lptilde(P°°) - lptilde(P°) + llikelihood(BBnew[i], P°°) - llikelihood(BB[i], P°)

        end
      #  println(ll)
        if rand() < exp(ll) * (piσ²(σ°.^2)/piσ²(σ.^2)) * (prod(σ°.^2)/prod(σ.^2))
            σ = σ°
            θ = θ°
            siacc += 1
            for i in 1:n
                BB[i] = BBnew[i]
            end
         #  print("acc")
        end
    end
    open(joinpath("output",simname,"params.txt"), "a") do f; println(f, iter, " ", join(round.([θ ; σ],8)," ")) end
    println(iter, "\t", join(round.([θ; σ./σtrue; 100bacc/iter/n;  100siacc/iter  ],3),"\t"))

    BBall = vcat([BB[i][1:end-1] for i in 1:n]...)

    mc = mcnext(mc, BBall.yy)
    mcparams = mcnext(mcparams, [θ ; σ])

    if iter % 10 == 0
        xr = (5,7)

        if PLOT
            plot(Yobs, "o", "o"; xrange=xr, yrange=(-3,3),linewidth=0.1)
            display(oplot(BBall, "b", "g"; xrange=xr, yrange=(-3,3), linewidth=0.7))
        end
    end
    if iter >= K break end
end

open(joinpath("output",simname,"info.txt"), "w") do f
    println(f, "n $n m $m T $TT A $Alpha B $Beta")
    println(f, "Y0 = $uu")
    println(f, "xi = $xi")
    println(f, "acc = ", 100bacc/iter/n, " (br) ",  100siacc/iter, " (si)")
    println(f, perf)
end

if PLOT

plot(Y, "r","b" ; yrange=(-3,3),linewidth=0.5)
oplot(Yobs,"+r", "+b")
savefig(joinpath("output",simname,"truth.pdf"))

plot(Y, "r","b" ; yrange=(-3,3),linewidth=0.5)
oplot(Yobs,"or","ob";symbolsize=0.3)
mcb = mcband(mc);
oplot(SamplePath(tt, mcb[1]),"r","b";linewidth=0.5)
oplot(SamplePath(tt, mcb[2]),"r","b";linewidth=0.5)
hcat(mcbandste(mcparams)..., [θtrue; σtrue])

savefig(joinpath("output",simname,"band.pdf"))

xr = (6,10)
plot(Y, "r","b" ; xrange=xr,yrange=(-3,3),linewidth=0.5)
oplot(Yobs,"or","ob";xrange=xr,symbolsize=0.3)
mcb = mcband(mc);
oplot(SamplePath(tt, mcb[1]),"r","b";xrange=xr,linewidth=0.5)
oplot(SamplePath(tt, mcb[2]),"r","b";xrange=xr,linewidth=0.5)
savefig(joinpath("output",simname,"bandpart.pdf"))

plot(vcat([BB[i] for i in 1:n]...), "r", "b"; yrange=(-3,3), linewidth=0.5)
oplot(Yobs,"+r", "+b")
savefig(joinpath("output",simname,"sample.pdf"))
end


text:
using YaoArrayRegister
using YaoArrayRegister: matvec

export ControlBlock, control, cnot

struct ControlBlock{N, BT<:AbstractBlock, C, M} <: AbstractContainer{BT, N}
    ctrl_locs::NTuple{C, Int}
    ctrl_config::NTuple{C, Int}
    content::BT
    locs::NTuple{M, Int}
    function ControlBlock{N, BT, C, M}(ctrl_locs, ctrl_config, block, locs) where {N, C, M, BT<:AbstractBlock}
        @assert_locs_safe N (ctrl_locs..., locs...)
        @assert nqubits(block) == M "number of locations doesn't match the size of block"
        @assert block isa AbstractBlock "expect a block, got $(typeof(block))"
        new{N, BT, C, M}(ctrl_locs, ctrl_config, block, locs)
    end
end

"""
    decode_sign(ctrls...)

Decode signs into control sequence on control or inversed control.
"""
decode_sign(ctrls::Int...,) = decode_sign(ctrls)
decode_sign(ctrls::NTuple{N, Int}) where N = tuple(ctrls .|> abs, ctrls .|> sign .|> (x->(1+x)÷2))

function ControlBlock{N}(ctrl_locs::NTuple{C}, ctrl_config::NTuple{C}, block::BT, locs::NTuple{K}) where {N, M, C, K, BT<:AbstractBlock{M}}
    M == K || throw(DimensionMismatch("block position not maching its size!"))
    return ControlBlock{N, BT, C, M}(ctrl_locs, ctrl_config, block, locs)
end

function ControlBlock{N}(ctrl_locs::NTuple{C}, ctrl_config::NTuple{C}, block, locs::NTuple{K}) where {N, M, C, K}
    error("expect a block, got $(typeof(block))")
end

# control bit configs are 1 by default, it use sign to encode control bit code
ControlBlock{N}(ctrl_locs::NTuple{C}, block::AbstractBlock, locs::NTuple) where {N, C} =
    ControlBlock{N}(decode_sign(ctrl_locs)..., block, locs)
ControlBlock{N}(ctrl_locs::NTuple{C}, block::Function, locs::NTuple) where {N, C} =
    ControlBlock{N}(decode_sign(ctrl_locs)..., parse_block(length(locs), block), locs)
ControlBlock{N}(ctrl_locs::NTuple{C}, block, locs::NTuple) where {N, C} =
    ControlBlock{N}(decode_sign(ctrl_locs)..., block, locs) # trigger error

# use pair to represent block under control in a compact way
ControlBlock{N}(ctrl_locs::NTuple{C}, target::Pair) where {N, C} =
    ControlBlock{N}(ctrl_locs, target.second, (target.first...,))

"""
    control(n, ctrl_locs, target)

Return a [`ControlBlock`](@ref) with number of active qubits `n` and control locs
`ctrl_locs`, and control target in `Pair`.

# Example

```jldoctest
julia> control(4, (1, 2), 3=>X)
nqubits: 4
control(1, 2)
└─ (3,) X gate

julia> control(4, 1, 3=>X)
nqubits: 4
control(1)
└─ (3,) X gate
```
"""
control(total::Int, ctrl_locs, target::Pair) = ControlBlock{total}(Tuple(ctrl_locs), target)
control(total::Int, control_location::Int, target::Pair) = control(total, (control_location, ), target)

"""
    control(ctrl_locs, target) -> f(n)

Return a lambda that takes the number of total active qubits as input. See also
[`control`](@ref).

# Example

```jldoctest
julia> control((2, 3), 1=>X)
(n -> control(n, (2, 3), 1 => X gate))

julia> control(2, 1=>X)
(n -> control(n, 2, 1 => X gate))
```
"""
control(ctrl_locs, target::Pair) = @λ(n -> control(n, ctrl_locs, target))
control(control_location::Int, target::Pair) = @λ(n -> control(n, control_location, target))

"""
    control(target) -> f(ctrl_locs)

Return a lambda that takes a `Tuple` of control qubits locs as input. See also
[`control`](@ref).

# Example

```jldoctest
julia> control(1=>X)
(ctrl_locs -> control(ctrl_locs, 1 => X gate))

julia> control((2, 3) => YaoBlocks.ConstGate.CNOT)
(ctrl_locs -> control(ctrl_locs, (2, 3) => CNOT gate))
```
"""
control(target::Pair) = @λ(ctrl_locs -> control(ctrl_locs, target))

"""
    control(ctrl_locs::Int...) -> f(target)

Return a lambda that takes a `Pair` of control target as input.
See also [`control`](@ref).

# Example

```jldoctest
julia> control(1, 2)
(target -> control((1, 2), target))
```
"""
control(ctrl_locs::Int...) = @λ(target -> control(ctrl_locs, target))

"""
    cnot(n, ctrl_locs, location)

Return a speical [`ControlBlock`](@ref), aka CNOT gate with number of active qubits
`n` and locs of control qubits `ctrl_locs`, and `location` of `X` gate.

# Example

```jldoctest
julia> cnot(3, (2, 3), 1)
nqubits: 3
control(2, 3)
└─ (1,) X gate

julia> cnot(2, 1)
(n -> cnot(n, 2, 1))
```
"""
cnot(total::Int, ctrl_locs, locs::Int) = control(total, ctrl_locs, locs=>X)
cnot(ctrl_locs, loc::Int) = @λ(n -> cnot(n, ctrl_locs, loc))

mat(::Type{T}, c::ControlBlock{N, BT, C}) where {T, N, BT, C} = cunmat(N, c.ctrl_locs, c.ctrl_config, mat(T, c.content), c.locs)

function apply!(r::ArrayReg{B, T}, c::ControlBlock) where {B, T}
    _check_size(r, c)
    instruct!(matvec(r.state), mat(T, c.content), c.locs, c.ctrl_locs, c.ctrl_config)
    return r
end

# specialization
for G in [:X, :Y, :Z, :S, :T, :Sdag, :Tdag]
    GT = Expr(:(.), :ConstGate, QuoteNode(Symbol(G, :Gate)))

    @eval function apply!(r::ArrayReg, c::ControlBlock{N, <:$GT}) where N
        _check_size(r, c)
        instruct!(matvec(r.state), Val($(QuoteNode(G))), c.locs, c.ctrl_locs, c.ctrl_config)
        return r
    end
end

PreserveStyle(::ControlBlock) = PreserveAll()

occupied_locs(c::ControlBlock) = (c.ctrl_locs..., map(x->c.locs[x], occupied_locs(c.content))...)
chsubblocks(pb::ControlBlock{N}, blk::AbstractBlock) where {N} = ControlBlock{N}(pb.ctrl_locs, pb.ctrl_config, blk, pb.locs)

# NOTE: ControlBlock will forward parameters directly without loop
cache_key(ctrl::ControlBlock) = cache_key(ctrl.content)

function Base.:(==)(lhs::ControlBlock{N, BT, C, M}, rhs::ControlBlock{N, BT, C, M}) where {BT, N, C, M}
    return (lhs.ctrl_locs == rhs.ctrl_locs) && (lhs.content == rhs.content) && (lhs.locs == rhs.locs)
end

Base.adjoint(blk::ControlBlock{N}) where N = ControlBlock{N}(blk.ctrl_locs, blk.ctrl_config, adjoint(blk.content), blk.locs)

# NOTE: we only copy one hierachy (shallow copy) for each block
function Base.copy(ctrl::ControlBlock{N, BT, C, M}) where {BT, N, C, M}
    return ControlBlock{N, BT, C, M}(ctrl.ctrl_locs, ctrl.ctrl_config, ctrl.content, ctrl.locs)
end

function YaoBase.iscommute(x::ControlBlock{N}, y::ControlBlock{N}) where N
    if x.locs == y.locs && x.ctrl_locs == y.ctrl_locs
        return iscommute(x.content, y.content)
    else
        return iscommute_fallback(x, y)
    end
end


text:
module UnetUtils

export Signals, Recordings

include("signals.jl")
include("recordings.jl")

end # module


text:
@enum DIRECTION RIGHT LEFT DOWN UP
function spiral_matrix(n::Int,m=n)
    spiral = zeros(Int,n,m)
    i = j = counter = 1
    direction = RIGHT
    for counter in 1:n*m
        spiral[i,j] = counter
        t = next_coord(direction,i,j)
        # If the next value is one we haven't set yet, continue on to set it
        if get(spiral,t,1) == 0
            (i,j) = t
        # Otherwise, we have gone out-of-bounds or hit a value we've set in
        # which case we should change direction
        else
            direction = next_dir(direction)
            (i,j) = next_coord(direction,i,j)
        end
    end
    spiral
end

function next_dir(c::DIRECTION)
    c == RIGHT ? DOWN :
    c == DOWN ? LEFT :
    c == LEFT ? UP :
    RIGHT # i.e. c==u
end

function next_coord(c::DIRECTION,i,j)
    c == RIGHT ? (i,j+1) :
    c == LEFT ? (i,j-1) :
    c == DOWN ? (i+1,j) :
    (i-1,j) # c==UP
end


text:
abstract type AbstractBinaryTree end
"""
array-based binary tree
"""
abstract type AbstractBinaryTree_arr<:AbstractBinaryTree end
mutable struct BinaryTree{T}<:AbstractBinaryTree_arr where T
    n::Int
    root::Int
    par::Vector{Int}
    lch::Vector{Int}
    rch::Vector{Int}
    val::Vector{T}
end
function BinaryTree{T}(size::Int,rootval::T)where T
    x=BinaryTree(1,1,zeros(T,size),zeros(T,size),zeros(T,size),Vector{T}(undef,size))
    x.val[1]=rootval
    return x
end

function Base.insert!(tree::BinaryTree,from::Int,value,place::Bool)
    tree.n+=1
    tree.par[tree.n]=from
    if place
        tree.rch[from]=tree.n
    else
        tree.lch[from]=tree.n
    end
    tree.val[tree.n]=value
end

function Base.insert!(tree::BinaryTree,from::Int,value)
    if tree.lch[from]==0
        insert!(tree,from,value,false)
    elseif tree.rch[from]==0
        insert!(tree,from,value,true)
    else
        return false
    end
    return true
end

# for AbstractBinaryTree_arr
function isleaf(tree::Tree,from::Int)::Bool where Tree<:AbstractBinaryTree_arr
    return tree.lch[from]==0&&tree.rch[from]==0
end

function ch(tree::Tree,from::Int,select::Bool)where Tree<:AbstractBinaryTree_arr
    if select return tree.rch[from]
    else return tree.lch[from]
    end
end

function set_ch(tree::Tree,from::Int,select::Bool,value::T)where {Tree<:AbstractBinaryTree_arr,T}
    if select tree.rch[from]=value
    else tree.lch[from]=value
    end
end

function height(tree::Tree,from::Int=tree.root)where Tree<:AbstractBinaryTree_arr
    count=0
    if tree.lch[from]!=0 count=height(tree,tree.lch[from])+1 end
    if tree.rch[from]!=0 count=max(height(tree,tree.rch[from])+1,count) end
    return count
end

function depth(tree::Tree,from::Int)where Tree<:AbstractBinaryTree_arr
    count=0
    while from!=tree.root
        count+=1
        from=tree.par[from]
    end
    return count
end

function left(tree::Tree,from::Int=tree.root)where Tree<:AbstractBinaryTree_arr
    if tree.lch[from]==0
        return from
    else
        return left(tree,tree.lch[from])
    end
end

function right(tree::Tree,from::Int=tree.root)where Tree<:AbstractBinaryTree_arr
    if tree.rch[from]==0
        return from
    else
        return right(tree,tree.rch[from])
    end
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text:
using Random

struct AliasSampler{IndexType <: Integer, ProbabilityParamType <: Real}
    # Some alias_indices will be left uninitialized after the
    # constructor finishes. This is not a bug!
    alias_indices::Vector{IndexType}
    nonalias_probs::Vector{ProbabilityParamType}
end

function AliasSampler(IndexType::Type,
                weights::AbstractVector{ProbabilityParamType}
            ) where{ProbabilityParamType <: Real}
    @assert IndexType <: Integer

    n = length(weights)
    s = sum(weights)
    scf = n/s
    sc_probabilities = (weights*scf)::Vector{ProbabilityParamType}

    smalls = Vector{IndexType}(undef, n)
    smalls_idx = 0
    larges = Vector{IndexType}(undef, n)
    larges_idx = 0

    for idx in 1:n
        if sc_probabilities[idx] >= 1.0
            larges_idx += 1
            larges[larges_idx] = idx
        else
            smalls_idx += 1
            smalls[smalls_idx] = idx
        end
    end

    aliases = Vector{IndexType}(undef, n)

    while larges_idx > 0 && smalls_idx > 0
        sm_p_idx = smalls[smalls_idx]
        lg_p_idx = larges[larges_idx]
        aliases[sm_p_idx] = lg_p_idx

        # This is slightly better numerically than the more obvious: sc_probabilities[lg_p_idx] += sc_probabilities[sm_p_idx] - 1.0
        sc_probabilities[lg_p_idx] -= 1.0
        sc_probabilities[lg_p_idx] += sc_probabilities[sm_p_idx]

        if sc_probabilities[lg_p_idx] < 1.0
            smalls[smalls_idx] = lg_p_idx
            larges_idx -= 1
        else
            smalls_idx -= 1
        end
    end

    while larges_idx > 0
        sc_probabilities[larges[larges_idx]] = 2.0
        larges_idx -= 1
    end

    while smalls_idx > 0
        sc_probabilities[smalls[smalls_idx]] = 2.0
        smalls_idx -= 1
    end

    return AliasSampler(aliases, sc_probabilities)
end

function asample(alias_sampler::AliasSampler{IndexType, ProbabilityParamType}, rng=Random.GLOBAL_RNG)::IndexType where{IndexType <: Integer, ProbabilityParamType <: Real}

    idx = rand(rng, 1:length(alias_sampler.alias_indices))
    if alias_sampler.nonalias_probs[idx] >= 1.0
        return idx
    end
    if rand(rng) < alias_sampler.nonalias_probs[idx]
        return idx
    else
        return alias_sampler.alias_indices[idx]
    end
end

# generalized setup:

function setup_alias_sampler!(
  weights::AbstractVector{T} where T <: Real,
  acceptances::AbstractVector{T} where T <: Real,
  aliases::AbstractVector{Idx},
  smalls::AbstractVector{Idx},
  larges::AbstractVector{Idx},
  wsum::Real = sum(weights)
  ) where Idx <: Integer

  n = length(weights)
  length(acceptances) == length(aliases) == length(smalls) == length(larges) == n || throw(DimensionMismatch("Inconsistent array lengths."))

  scf = n/wsum
  for i = 1:n
    @inbounds acceptances[i] = weights[i] * scf
  end

  smalls_idx = 0
  larges_idx = 0

for idx in 1:n
    if acceptances[idx] >= 1.0
      larges_idx += 1
      larges[larges_idx] = idx
    else
      smalls_idx += 1
      smalls[smalls_idx] = idx
    end
  end

  while larges_idx > 0 && smalls_idx > 0
    sm_p_idx = smalls[smalls_idx]
    lg_p_idx = larges[larges_idx]
    aliases[sm_p_idx] = lg_p_idx

    # This is slightly better numerically than the more obvious: sc_probabilities[lg_p_idx] += sc_probabilities[sm_p_idx] - 1.0
    acceptances[lg_p_idx] -= 1.0
    acceptances[lg_p_idx] += acceptances[sm_p_idx]

    if acceptances[lg_p_idx] < 1.0
      smalls[smalls_idx] = lg_p_idx
      larges_idx -= 1
    else
      smalls_idx -= 1
    end
  end

  while larges_idx > 0
    acceptances[larges[larges_idx]] = 2.0
    larges_idx -= 1
  end

  while smalls_idx > 0
    acceptances[smalls[smalls_idx]] = 2.0
    smalls_idx -= 1
  end
  nothing
end


text:
function secret_handshake(code)
    bits=bitstring(code)
    uncode=reverse(bits)
    ret = String[]
    if uncode[1]=='1'
        push!(ret,"wink")
    end
    if uncode[2]=='1'
        push!(ret,"double blink")
    end
    if(uncode[3]=='1')
        push!(ret,"close your eyes")
    end
    if(uncode[4]=='1')
        push!(ret,"jump")
    end

    if(uncode[5]=='1')
        ret=reverse(ret)
    end
    return ret
end

max_stars_repo_name: stepinski/julia-in

max_issues_repo_name: stepinski/julia-init

max_forks_repo_name: stepinski/julia-init


text:
using GLPlot, GLAbstraction, ModernGL, GeometryTypes

N = 10
r = linspace(0, 6, N)
directions = Vec3f0[(sin(x), cos(y), sin(z)) for x = r,y = r, z = r]
obj = glplot(directions)


text:
using ChainRulesCore
using ChainRulesTestUtils
using BenchmarkTools

function two2three(x1::Float64, x2::Float64)
    return 1.0, 2.0 * x1, 3.0 * x2
end

function ChainRulesCore.frule((Δf, Δx1, Δx2), ::typeof(two2three), x1, x2)
    y = two2three(x1, x2)
    ∂y = Tangent{Tuple{Float64, Float64, Float64}}(
        ZeroTangent(),
        2.0 * Δx1,
        3.0 * Δx2,
    )
    return y, ∂y
end

function ChainRulesCore.rrule(::typeof(two2three), x1, x2)
    y = two2three(x1, x2)
    function two2three_pullback(Ȳ)
        return (NoTangent(), 2.0 * Ȳ[2], 3.0 * Ȳ[3])
    end
    return y, two2three_pullback
end

test_frule(two2three, 3.33, -7.77)

##
function two2three(x::AbstractArray)
    return 1.0, 2.0 * x[1], 3.0 * x[2]
end

function ChainRulesCore.frule((Δf, Δx), ::typeof(two2three), x)
    y = two2three(x)
    ∂y = Tangent{Tuple{Float64, Float64, Float64}}(
        ZeroTangent(),
        2.0 * Δx[1],
        3.0 * Δx[2],
    )
    return y, ∂y
end

function ChainRulesCore.rrule(::typeof(two2three), x)
    y = two2three(x)
    function two2three_pullback(Ȳ)
        return (NoTangent(), [2.0 * Ȳ[2], 3.0 * Ȳ[3]])
    end
    return y, two2three_pullback
end

vec = [3.33, -7.77]
two2three(vec)
test_frule(two2three, vec)
test_rrule(two2three, vec)

##
function ntwo2three(x::AbstractArray)
    return 1.0, 2.0 * x[1]^2, 3.0 * x[2]^3
end

function ChainRulesCore.frule((Δf, Δx), ::typeof(ntwo2three), x)
    y = ntwo2three(x)
    ∂y = Tangent{Tuple{Float64, Float64, Float64}}(
        ZeroTangent(),
        4.0 * x[1] * Δx[1],
        9.0 * x[2]^2 * Δx[2],
    )
    return y, ∂y
end

function ChainRulesCore.rrule(::typeof(ntwo2three), x)
    y = ntwo2three(x)
    function ntwo2three_pullback(Ȳ)
        return (NoTangent(), [4.0 * x[1] * Ȳ[2], 9.0 * x[2]^2 * Ȳ[3]])
    end
    return y, ntwo2three_pullback
end

vec = [3.33, -7.77]
two2three(vec)
test_frule(ntwo2three, vec)
test_rrule(ntwo2three, vec)

##
function lorenz(s; σ = 10, ρ = 28, β = 8 / 3)
    x, y, z = s

    ẋ = σ * (x - y)
    ẏ = -y + (ρ - z) * x
    ż = -β * z + x * y

    return [ẋ, ẏ, ż]
end

using Zygote
tmp = Zygote.jacobian(lorenz, [1.0, 2.0, 3.0])
@show tmp

function lorenz!(ṡ, s; σ = 10, ρ = 28, β = 8 / 3)
    x, y, z = s
    ṡ[1] = -σ * (x - y)
    ṡ[2] = -y + (ρ - z) * x
    ṡ[3] = -β * z + x * y
    return nothing
end

function jlorenz!(ṡ, s, sᵃ; σ = 10, ρ = 28, β = 8 / 3)
    x, y, z = s
    xᵃ, yᵃ, zᵃ = sᵃ
    ṡ[1] = σ * (xᵃ - yᵃ)
    ṡ[2] = -yᵃ + (ρ - zᵃ) * x - z * xᵃ
    ṡ[3] = -β * zᵃ + xᵃ * y + x * yᵃ
    return nothing
end

function jlorenz(s, sᵃ; σ = 10, ρ = 28, β = 8 / 3)
    return Zygote.jacobian(lorenz, s)[1] * sᵃ
end

ṡ = [1.0, 2.0, 3.0]
s = copy(ṡ)
y = copy(ṡ)

@benchmark lorenz!(ṡ, s)
@benchmark lorenz(s)

@benchmark jlorenz!(ṡ, s, y)
@benchmark jlorenz(s, y)

function step!(s, ṡ, s̃, rhs!, Δt)
    rhs!(ṡ, s)
    @. s̃ = s + Δt * ṡ
    @. s = s + Δt * 0.5 * ṡ
    rhs!(ṡ, s̃)
    @. s = s + Δt * 0.5 * ṡ
    return nothing
end

function step(s, rhs, Δt)
    ṡ = rhs(s)
    s̃ = s + Δt * ṡ
    s2 = rhs(s̃)
    return s + Δt * 0.5 * (ṡ + s2)
end

function stepᵃ!(s, ṡ, s̃, sᵃ, i, rhs!, Δt)
    rhs!(ṡ, s[i], sᵃ)
    @. s̃ = s + Δt * ṡ
    @. s = s + Δt * 0.5 * ṡ
    rhs!(ṡ, s[i + 1], s̃)
    @. s = s + Δt * 0.5 * ṡ
    return nothing
end

s̃ = copy(s)
@benchmark step!(s, ṡ, s̃, lorenz!, 0.1)
@benchmark step(s, lorenz, 0.1)

s .= [1.0, 1.0, 1.0]
step(s, lorenz, 0.001)

function evolve(s, N)
    y = s
    for _ in 1:N
        y = step(y, lorenz, 0.001)
    end
    return y
end

function evolve!(s, ṡ, s̃, N)
    for _ in 1:N
        step!(s, ṡ, s̃, lorenz!, 0.001)
    end
    return nothing
end

@benchmark evolve(s, 10)
@benchmark evolve!(s, ṡ, s̃, 10)

max_stars_repo_name: sandreza/GradientLearning

max_issues_repo_name: sandreza/GradientLearning

max_forks_repo_name: sandreza/GradientLearning


text:
export Source, Receiver, Terrain

# ----- #
# Geometry #
# ----- #

"""
    Geometry
All objects describing the model geometry should be a subtype of this.
This includes three objects: Source, Reciver and Terrain.
"""
abstract type Geometry end

# --------------- #
# - Source - #
# --------------- #

"""
    Source
A source of sound which emits sound rays from a point in 2D space.
### Fields:
* `frequency`  - Hz, frequency of source sound emitted from the source.
* `height`     - m, height of the sound source above local ground level.
### Available Constructors:
* `Source(;frequency::Real = 1000.0f0, height::Real = 80.0f0))`
* `Source(frequency, height)`
"""
struct Source{T<:Real} <: Geometry
    frequency::T
    height::T
end

Source(;frequency::Real = 1000.0f0, height::Real = 80.0f0) =
    Source(frequency, height)

# --------------- #
# - Receiver - #
# --------------- #

"""
    Receiver
Multiple receiver in 2D space.
### Fields:
* `depth_point`  - Number of gird points along depth (or height).
* `range_point`  - Number of gird points along propagation range.
* `depth`        - from lowest to highest depth.
* `range`  - from source to receiver.
### Available Constructors:
* `Receiver(;depth_point::Int = 100,
            range_point::Int = 1000,
             depth = Vec2(0.0f0, 1.0f0),
             range = Vec2(0.0f0, 1.0f0))`
* `Receiver(depth_point, range_point, depth,range)`
"""
struct Receiver{T<:AbstractArray} <: Geometry
    depth_point:: Int
    range_point:: Int
    depth::Vec2{T}
    range::Vec2{T}
end

Receiver(;depth_point::Int = 100,
            range_point::Int = 1000,
             depth = Vec2(0.0f0, 1.0f0),
             range = Vec2(0.0f0, 1.0f0)) =
 Receiver(depth_point, range_point, depth,range)

 # --------------- #
 # - Terrain - #
 # --------------- #

 """
     Terrain
 Elevation profile.
 ### Fields:
 * `interp_type`    - Interpolation method 'L' for piecewise linear fit and
                    'C' for a curvilinear fit.
 * `profile`        - x- range in km and y-elevation in m.
 ### Available Constructors:
 * `Terrain(;interp_type::String = "L",
              profile = Vec2(0.0f0, 1.0f0))`
 * `Terrain(interp_type, profile)`
 """
struct Terrain{T<:AbstractArray} <: Geometry
    interp_type::String
    profile::Vec2{T}
end

Terrain(;interp_type::String = "L",
             profile = Vec2(0.0f0, 1.0f0)) =
 Terrain(interp_type, profile)


text:
module MemoryInspectorTest
import MemoryInspector
using Test

mutable struct Obj{T1,T2}
    a::T1
    b::T2
    Obj{A,B}(a::A,b::B) where {A,B} = new{A,B}(a,b)
    Obj{A,B}(a::A) where {A,B} = new{A,B}(a) # Allow partial construction
    Obj{A,B}() where {A,B} = new{A,B}() # Allow partial construction
    Obj(a::A,b::B) where {A,B} = Obj{A,B}(a,b)
end

@testset "self-references" begin
    s = join(rand('a':'z', 1024*1024)); # 1 MiB string
    multistring = Obj(s,s);   # Shared references to large object
    selfie = Obj{Obj,Obj}(multistring);
    selfie.b = selfie;    # Self-reference

    # Not sure how to test an interactive UI like this.
    #MemoryInspector.@inspect selfie
    result = MemoryInspector.MemorySummarySize.summarysize(selfie)

    @test result.size == Base.summarysize(selfie)
    @test length(result.children) == 2

    child_sizes = [v.size for (k,v) in result.children]
    # Test that all the bytes are accounted for: the multistring child and the self-reference.
    @test sum(child_sizes) + sizeof(selfie) == result.size  # allocated children + pointers.
    @test Base.summarysize(selfie.a) in child_sizes
    @test 0 in child_sizes  # selfie.b is a self-reference so adds no size besides its pointers.
end

@testset "handles distinct mutable structs with same value" begin
    double_arr = Obj([1,2], [1,2])

    result = MemoryInspector.MemorySummarySize.summarysize(double_arr)
    @test result.size == Base.summarysize(double_arr)
end
@testset "handles distinct immutable structs with same value" begin
    double_arr = Obj((1,2), (1,2))

    result = MemoryInspector.MemorySummarySize.summarysize(double_arr)
    @test result.size == Base.summarysize(double_arr)
end

@testset "Complex Object" begin
    result = MemoryInspector.MemorySummarySize.summarysize(Base)
    # I'm undercounting some things here... I think I'm being over-generous with self-references,
    # counting things that have the same _value_ but aren't _actually_ self references? I'm not
    # sure though.
    @test_broken result.size == Base.summarysize(Base)
end

end


text:
int_rules_6_1_10 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 6.1.10*(c+d*x)^m*(a+b*sinh)^n =#
    @apply_utils Antiderivative((~u) ^ ~(m') * (~(a') + ~(b') * sinh(~v)) ^ ~(n'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * (~a + ~b * sinh(ExpandToSum(~v, ~x))) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (LinearQ([~u, ~v], ~x) && Not(LinearMatchQ([~u, ~v], ~x)))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~(a') + ~(b') * cosh(~v)) ^ ~(n'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * (~a + ~b * cosh(ExpandToSum(~v, ~x))) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (LinearQ([~u, ~v], ~x) && Not(LinearMatchQ([~u, ~v], ~x)))
end

max_stars_repo_path: src/rules/6 Hyperbolic functions/6.1 Hyperbolic sine/6.1.10 (c+d x)^m (a+b sinh)^n.jl


text:
using Colors, ColorSchemes
using CairoMakie, FileIO
using CSV, DataFrames
using Statistics

imgpath = joinpath("Molenaar", "Figure 4", "Plots")
resultspath = joinpath("Molenaar", "Results")

x = [res["mu"] for res in stats]
a = [res["cat_flux"] for res in statsb]
b = [res["met_flux"] for res in statsb]

fig = Figure()
ax = Axis(fig[1, 1])
scatter1 = scatter!(ax, x, a)
line1 = lines!(ax, x, a)
scatter2 = scatter!(ax, x, b)
line2 = lines!(ax, x, b)

ax.xlabel = "Growth rate"
ax.ylabel = "Fraction of the total substarte flux"
fig

FileIO.save(joinpath(imgpath, "Figure_4B.pdf"), fig)


text:
using DaemonMode
using Documenter

makedocs(;
    modules=[DaemonMode],
    authors="Daniel Molina <dmolina@decsai.ugr.es>",
    repo="https://github.com/dmolina/MoodleQuestions.jl/blob/{commit}{path}#L{line}",
    sitename="DaemonMode.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://dmolina.github.io/DaemonMode.jl",
        assets=String[],
    ),
    pages=[
        "Quick Start" => "index.md",
        "Posibilities" => "posibilities.md",
        "User Guide" => "guide.md",
        "Public API" => "reference.md"
    ],
)

deploydocs(;
    repo="github.com/dmolina/DaemonMode.jl",
)


text:
doc"""
    Gumbel(μ, θ)

The *Gumbel distribution*  with location `μ` and scale `θ` has probability density function

$f(x; \mu, \theta) = \frac{1}{\theta} e^{-(z + e^z)},
\quad \text{ with } z = \frac{x - \mu}{\theta}$

```julia
Gumbel()            # Gumbel distribution with zero location and unit scale, i.e. Gumbel(0.0, 1.0)
Gumbel(u)           # Gumbel distribution with location u and unit scale, i.e. Gumbel(u, 1.0)
Gumbel(u, b)        # Gumbel distribution with location u and scale b

params(d)        # Get the parameters, i.e. (u, b)
location(d)      # Get the location parameter, i.e. u
scale(d)         # Get the scale parameter, i.e. b
```

External links

* [Gumbel distribution on Wikipedia](http://en.wikipedia.org/wiki/Gumbel_distribution)
"""
immutable Gumbel <: ContinuousUnivariateDistribution
    μ::Float64  # location
    θ::Float64  # scale

    Gumbel(μ::Real, θ::Real) = (@check_args(Gumbel, θ > zero(θ)); new(μ, θ))
    Gumbel(μ::Real) = new(μ, 1.0)
    Gumbel() = new(0.0, 1.0)
end

@distr_support Gumbel -Inf Inf

const DoubleExponential = Gumbel

#### Parameters

location(d::Gumbel) = d.μ
scale(d::Gumbel) = d.θ
params(d::Gumbel) = (d.μ, d.θ)

#### Statistics

mean(d::Gumbel) = d.μ + d.θ * 0.57721566490153286

median(d::Gumbel) = d.μ + d.θ * 0.366512920581664327

mode(d::Gumbel) = d.μ

var(d::Gumbel) = 1.6449340668482264 * d.θ^2

skewness(d::Gumbel) = 1.13954709940464866

kurtosis(d::Gumbel) = 2.4

entropy(d::Gumbel) = 1.57721566490153286 + log(d.θ)

#### Evaluation

zval(d::Gumbel, x::Float64) = (x - d.μ) / d.θ
xval(d::Gumbel, z::Float64) = x * d.θ + d.μ

function pdf(d::Gumbel, x::Float64)
    z = zval(d, x)
    exp(-z - exp(-z)) / d.θ
end

function logpdf(d::Gumbel, x::Float64)
    z = zval(d, x)
    - (z + exp(-z) + log(d.θ))
end

cdf(d::Gumbel, x::Float64) = exp(-exp(-zval(d, x)))
logcdf(d::Gumbel, x::Float64) = -exp(-zval(d, x))

quantile(d::Gumbel, p::Float64) = d.μ - d.θ * log(-log(p))

gradlogpdf(d::Gumbel, x::Float64) = - (1.0 + exp((d.μ - x) / d.θ)) / d.θ

#### Sampling

rand(d::Gumbel) = quantile(d, rand())


text:
# include generic methods
include("pred_ared.jl")
include("step_computation.jl")

# # include ARC methods
include("ARC1D/ARC-Cub.jl")
include("ARC1D/ARC-Nwt.jl")
include("ARC1D/ARC-Sec.jl")
include("ARC1D/ARC-SecA.jl")

include("ARC1D/ARC_generic.jl")
include("ARC1D/ARC_direction_computation.jl")

# include Trust Region methods
include("TR1D/TR-Cub.jl")
include("TR1D/TR-Nwt.jl")
include("TR1D/TR-Sec.jl")
include("TR1D/TR-SecA.jl")

include("TR1D/TR_generic.jl")
include("TR1D/TR_direction_computation.jl")

# # include zoom
include("Zoom1D/Trouve_IntervalleA.jl")
include("Zoom1D/zoom_generic.jl")
include("Zoom1D/zoom-Cub.jl")
include("Zoom1D/zoom-Nwt.jl")
include("Zoom1D/zoom-Sec.jl")
include("Zoom1D/zoom-SecA.jl")

# # include bissection
include("Bissection/Trouve_Intervalle.jl")
include("Bissection/Bissect.jl")
include("Bissection/BissectNwt.jl")
include("Bissection/BissectSec.jl")
include("Bissection/BissectSecA.jl")
include("Bissection/BissectCub.jl")


text:
struct AlternatingDescent <: AbstractLoss end

function complete(D::MaskedEDM{T}, loss::AlternatingDescent, d::Int) where T
    D̂, _ = alternating_descent(D, dims=d)

    return D̂
end

# function complete(D::NoisyMaskedEuclideanDistanceMatrix, loss::AlternatingDescent, d::Int)
#     Dm =

# end

"""
    alternating_descent(D::MaskedEDM{T}; dims::Int = 2, max_iterations::Int=50) where T <: Real

Complete an EMD using Alternating Descent.

Authors: Reza Parhizkar and Ivan Dokmanic, 2014. Karel Mundnich (port to Julia), 2019
"""
function alternating_descent(D::MaskedEDM{T}; dims::Int = 2, max_iterations::Int = 50) where T <: Real

    n = size(D, 1)
    L = J(n)
    Dm = D.D .* D.mask # We operate over this matrix

    X₀ = zeros(n, dims) # TODO: Transpose
    X = X₀

    connected_nodes = Array{Array{Int,1},1}(undef, n)
    connectivity_vector = Array{Array{Float64,1},1}(undef, n)

    for i in 1:n
        connected_nodes[i] = findall(x -> x != 0, Dm[i, :])
        connectivity_vector[i] = Dm[i, connected_nodes[i]]
    end

    for iter_ind in 1:max_iterations, ind in 1:n
        X_connected = X[connected_nodes[ind], :]

        for coor_ind in 1:dims
            a = 4 * size(X_connected, 1)
            b = 12 * sum((broadcast(-, X[ind, coor_ind], X_connected[:,coor_ind])))
            c = 4 * sum(sum((repeat(X[ind:ind,:], size(X_connected, 1), 1) .- X_connected).^2, dims=2)
                .+ 2 * (X[ind,coor_ind] .- X_connected[:,coor_ind]).^2 .- connectivity_vector[ind])
            d = 4 * sum((X[ind,coor_ind] .-  X_connected[:,coor_ind]) .* (sum((repeat(X[ind:ind,:], size(X_connected, 1), 1) .- X_connected).^2, dims=2) - connectivity_vector[ind]))

            DeltaX_vec = roots(a, b, c, d)
            if DeltaX_vec != nothing
                # If we have no solutions, we should break and return
                DeltaX_vec = real(DeltaX_vec[abs.(imag(DeltaX_vec)) .< eps()])
            else
                break
            end

            cc = zeros(Float64, length(DeltaX_vec))

            for ii = 1 : length(DeltaX_vec)
                cc[ii] = sum(((X[ind, coor_ind] + DeltaX_vec[ii] .- X_connected[:,coor_ind]).^2
                       + sum((repeat(X[ind:ind,:], size(X_connected, 1), 1) .- X_connected).^2, dims=2)
                       - (X[ind, coor_ind] .- X_connected[:, coor_ind]).^2
                       - connectivity_vector[ind]).^2)
            end

            DeltaX = DeltaX_vec[argmin(cc)]
            X[ind, coor_ind] = X[ind, coor_ind] + DeltaX

            X = L * X
        end
    end

    X = X'L
    D = pairwise(SqEuclidean(), X, dims=2)

    return EuclideanDistanceMatrix(D), Embedding(X)

end

function alternating_descent(D::NoisyMaskedEDM{T}; dims::Int = 2, max_iterations::Int = 50) where T <: Real

end

"""
    roots(a::T, b::T, c::T) where T <: Real

Solve a quadratic equation where a, b, and c are real

    a*x^2 + b*x + c = 0

Authors:  Nam Sun Wang (Matlab version), Karel Mundnich (Julia port)
"""
function roots(a::Real, b::Real, c::Real)
    if a == 0
        if b == 0
            # We have a non-equation therefore, we have no valid solution
            return nothing
        else
            # We have a linear equation with 1 root.
            return -c/b
        end
    else
        # We have a true quadratic equation.  Apply the quadratic formula to find two roots.
        x = zeros(Float64, 2)
        DD = b * b - 4 * a * c
        x[1] = (-b + sqrt(DD))/2/a
        x[2] = (-b - sqrt(DD))/2/a
        return x
    end
end

"""
        roots(a::T, b::T, c::T, d::T) where T <: Real

Solve a cubic equation where a, b, c, and d are real.

    a*x^3 + b*x^2 + c*x + d = 0

# Reference
Tuma, "Engineering Mathematics Handbook", pp, 7, McGraw Hill, 1978.

# Algorithm

        - Step 0: If a is 0, use the quadratic formula to avoid dividing by 0.
        - Step 1: Calculate p and q
                    p = ( 3*c/a - (b/a)**2 ) / 3
                    q = ( 2*(b/a)**3 - 9*b*c/a/a + 27*d/a ) / 27
        - Step 2: Calculate discriminant D
                    D = (p/3)**3 + (q/2)**2
        - Step 3: Depending on the sign of D, we follow different strategy.
                    If D<0, thre distinct real roots.
                    If D=0, three real roots of which at least two are equal.
                    If D>0, one real and two complex roots.
        - Step 3a: For D>0 and D=0,
                    Calculate u and v
                    u = cubic_root(-q/2 + sqrt(D))
                    v = cubic_root(-q/2 - sqrt(D))
                    Find the three transformed roots
                    y1 = u + v
                    y2 = -(u+v)/2 + i (u-v)*sqrt(3)/2
                    y3 = -(u+v)/2 - i (u-v)*sqrt(3)/2
        - Step 3b Alternately, for D<0, a trigonometric formulation is more convenient
                    y1 =  2 * sqrt(|p|/3) * cos(ϕ/3)
                    y2 = -2 * sqrt(|p|/3) * cos((ϕ+pi)/3)
                    y3 = -2 * sqrt(|p|/3) * cos((ϕ-pi)/3)
                    where ϕ = acos(-q/2/sqrt(|p|**3/27))
                                pi  = 3.141592654...
        - Step 4  Finally, find the three roots
                    x = y - b/a/3

Authors: Nam Sun Wang (Matlab version), Karel Mundnich (Julia port)
"""
function roots(a::Real, b::Real, c::Real, d::Real)
    # Step 0: If a is 0 use the quadratic formula
    if a == 0.0
        return roots(b, c, d)
    end

    # Cubic equation with 3 roots
    nroot = 3

    # Step 1: Calculate p and q
    p  = c / a - b * b / a / a / 3.0
    q  = (2.0 * b * b * b / a / a / a - 9.0 * b * c / a / a + 27.0 * d / a) / 27.0

    # Step 2: Calculate DD (discriminant)
    DD = p * p * p / 27.0 + q * q / 4.0

    # Step 3: Branch to different algorithms based on DD
    if DD < 0.0
        # Step 3b
        # 3 real unequal roots -- use the trigonometric formulation
        ϕ = acos(- q / 2.0 / sqrt(abs(p * p * p) / 27.0))
        temp1 = 2.0 * sqrt(abs(p)/3.0)
        y1 =  temp1 * cos(ϕ/3.0)
        y2 = -temp1 * cos((ϕ + pi)/3.0)
        y3 = -temp1 * cos((ϕ - pi)/3.0)
    else
        # Step 3a:
        # 1 real root & 2 conjugate complex roots OR 3 real roots (some are equal)
        temp1 = -q/2.0 + sqrt(DD)
        temp2 = -q/2.0 - sqrt(DD)
        u = abs(temp1)^(1/3)
        v = abs(temp2)^(1/3)
        if (temp1 < 0.0) u = -u end
        if (temp2 < 0.0) v = -v end
        y1  = u + v
        y2r = -(u+v)/2.0
        y2i =  (u-v) * sqrt(3.0)/2.0
    end

    # Step 4: Final transformation
    temp1 = b / a / 3.0
    y1 = y1 - temp1
    if DD < 0.0
        y2 = y2 - temp1
        y3 = y3 - temp1
    else
        y2r = y2r - temp1
    end

    # Assign answers
    x = zeros(Number, 3)

    if DD < 0.0
        x[1] = y1
        x[2] = y2
        x[3] = y3
    elseif DD == 0.0
        x[1] =  y1
        x[2] = y2r
        x[3] = y2r
    else
        x[1] = y1
        x[2] = y2r + y2i * 1im
        x[3] = y2r - y2i * 1im
    end

return x
end


text:
## Algorithm to compute/approximate the MLA values from extended multiroute flow

# Get the slope of the restricted flow function at a given restriction x
function getSlope(g::Network,restriction::Float64,cutedges::Vector{Int})
  slope = 0
  for id in cutedges
    if g.links[id].cap > restriction
      slope += 1
    end
  end
  return slope
end
function getFlowSlope(g::Network,g0::Network,restriction::Float64)
  flow, cutedges = mincut(g)
  slope = getSlope(g0,restriction,cutedges)
  return flow, slope, cutedges
end

# Make a copy of a network with a restricted value on the links capacity
function restriction(x::Float64,g::Network)
  function restrict(a::Arc)
    Arc(a.sym,a.head,min(x,a.cap))
  end
  arcs = map(restrict,g.links)
  new_g = Network(g.source,g.target,g.tails,arcs)
  return new_g
end

function intersection(x1,y1,a1,x2,y2,a2)
  b1 = y1 - a1 * x1
  b2 = y2 - a2 * x2
  x = (b2 - b1)/(a1 - a2)
  y = a1 * x + b1
#  println("$y1 = $a1 * $x1 + $b1")
#  println("$y2 = $a2 * $x2 + $b2")
#  println("x=$x y=$y")
  return (x,y)
end

function auxiliary_points(g::Network)
  λ = round(Int,connectivity(g))
  n = length(g.tails)
  auxpoints = fill((0.,0.),λ+1)
  cuts = Vector{Vector{Int}}()

  resmin, resmax = minmaxCapacity(g)
  resmin = resmin/2.
  gmin = restriction(resmin,g)
#  print(gmin)
  fmin, smin, cutedges = getFlowSlope(gmin,g,resmin)
#  println("smin=",smin,", resmin=",resmin,", fmin=",fmin,", λ=",λ)
  auxpoints[λ + 1 - smin] = (resmin,fmin)
  push!(cuts,cutedges)

  fmax, smax, cutedges = getFlowSlope(g,g,resmax)
  auxpoints[λ + 1 - smax] = (resmax,fmax)
  push!(cuts,cutedges)

  if smin > smax + 1
    queue = [((fmin,smin,resmin),(fmax,smax,resmax))]
  end

  while !isempty(queue)
    (f1,s1,r1),(f2,s2,r2) = pop!(queue)
    res, expectslope = intersection(r1,f1,s1,r2,f2,s2)
    g_res = restriction(res,g)
    f_res, s_res, cutedges = getFlowSlope(g_res,g,res)
    auxpoints[λ + 1 - s_res] = (res,f_res)
    push!(cuts,cutedges)

    if (s1 > s_res + 1) && ((r2,f2) != (res,f_res))
      q = (f1,s1,r1),(f_res,s_res,res)
      push!(queue,q)
    end
    if (s_res > s2 + 1) && ((r1,f1) != (res,f_res))
      q = (f_res,s_res,res),(f2,s2,r2)
      push!(queue,q)
    end
  end
  return auxpoints, cuts
end

function breakingPoints(g::Network)
  auxpoints, cuts = auxiliary_points(g)
  λ = length(auxpoints) - 1
  breakingpoints = Vector{Tuple{Float64,Float64,Int}}()
  left_index = 1

  for (id,p) in enumerate(auxpoints)
    if id == 1
      push!(breakingpoints,(0.,0.,λ))
    else
      p_left = breakingpoints[left_index]
#      println("p=",p)
      if p[1] != 0
        interx, intery = intersection(p_left[1],p_left[2],p_left[3],p[1],p[2], λ + 1 - id)
        push!(breakingpoints,(interx,intery,λ + 1 - id))
        left_index += 1
      else
        push!(breakingpoints,(-1.,-1.,λ + 1 - id))
      end
    end
  end
#  println(auxpoints)
#  println(breakingpoints)
  return breakingpoints
end

function breakingPointsCuts(g::Network)
  auxpoints, cuts = auxiliary_points(g)
  λ = length(auxpoints) - 1
  breakingpoints = Vector{Tuple{Float64,Float64,Int}}()
  left_index = 1

  for (id,p) in enumerate(auxpoints)
    if id == 1
      push!(breakingpoints,(0.,0.,λ))
    else
      p_left = breakingpoints[left_index]
#      println("p=",p)
      if p[1] != 0
        interx, intery = intersection(p_left[1],p_left[2],p_left[3],p[1],p[2], λ + 1 - id)
        push!(breakingpoints,(interx,intery,λ + 1 - id))
        left_index += 1
      else
        push!(breakingpoints,(-1.,-1.,λ + 1 - id))
      end
    end
  end
#  println(auxpoints)
#  println(breakingpoints)
  return breakingpoints, cuts
end


text:
using  Utils

const HALF = 16

function parse_file(f)
    size = fill(2 * HALF, 3)
    re = zeros(Bool, size...)

    for (row, line) in enumerate(eachline(f))
        for (col, c) in enumerate(line)
            re[row+HALF,col+HALF,1+HALF] = c == '#'
        end
    end
    re
end

function load()
    open(joinpath(@__DIR__, "17.txt"), "r") do f
        parse_file(f)
    end
end

function neighbors(cubes, r, c, s)
    count(view(cubes, r-1:r+1, c-1:c+1, s-1:s+1)) - Int(cubes[r,c,s])
end

function neighbors(cubes, r, c, s, d)
    count(view(cubes, r-1:r+1, c-1:c+1, s-1:s+1, d-1:d+1)) - Int(cubes[r,c,s,d])
end

function step(cubes::AbstractArray{Bool,3})
    same_dim = axes(cubes, 1)
    range = firstindex(same_dim)+1:lastindex(same_dim)-1

    prev = copy(cubes)
    for slice in range
        for col in range
            for row in range
                ns = neighbors(prev, row, col, slice)
                if prev[row,col,slice] && !(ns == 2 || ns == 3)
                    cubes[row,col,slice] = false
                elseif !prev[row,col,slice] && ns == 3
                    cubes[row,col,slice] = true
                end
            end
        end
    end
end

function step(cubes::AbstractArray{Bool,4})
    same_dim = axes(cubes, 1)
    range = firstindex(same_dim)+1:lastindex(same_dim)-1

    prev = copy(cubes)
    for dimension in range
        for slice in range
            for col in range
                for row in range
                    ns = neighbors(prev, row, col, slice, dimension)
                    if prev[row,col,slice,dimension] && !(ns == 2 || ns == 3)
                        cubes[row,col,slice,dimension] = false
                    elseif !prev[row,col,slice,dimension] && ns == 3
                        cubes[row,col,slice,dimension] = true
                    end
                end
            end
        end
    end
end

function ex1(cubes, to=6)
    cubes = copy(cubes)
    for i in 1:to
        step(cubes)
    end
    count(cubes)
end

function ex2(cubes, to=6)
    more = zeros(Bool, size(cubes)..., 2 * HALF)
    more[:,:,:,1+HALF] = cubes
    ex1(more, to)
end

cubes = load()

println("Alive after 6 rounds: ", ex1(cubes))
println("Alive after 6 hypercube rounds: ", ex2(cubes))

using Test

@testset "Conway Cubes" begin
    input = """
    .#.
    ..#
    ###
    """

    cubes = parse_file(IOBuffer(input))

    @testset "example 1" begin
        @test ex1(cubes) == 112
    end

    @testset "example 2" begin
        @test ex2(cubes, 0) == 5
        @test ex2(cubes, 1) == 29
        @test ex2(cubes) == 848
    end

    @testset "results" begin
    end
end

max_stars_repo_name: CmdQ/AoC


text:
"""
    read(io::IO, ::Type{T}; compress = :auto) where {T<:MRCData}
    read(fn::AbstractString, ::Type{T}; compress = :auto) where {T<:MRCData}

Read an instance of [`MRCData`](@ref) from an IO stream or existing file.
Use `compress` to specify the decompression with the following options:
- `:auto`: Infer one of the below compression types from the file
- `:gz`: GZ
- `:bz2`: BZ2
- `:xz`: XZ
- `:none`: no compression
"""
read(::Any, ::Type{<:MRCData})
function Base.read(io::IO, ::Type{T}; compress=:auto) where {T<:MRCData}
    if compress == :auto
        compress = checkmagic(io)
    end
    newio = decompressstream(io, compress)
    header = read(newio, MRCHeader)
    extendedheader = read(newio, MRCExtendedHeader; header=header)
    d = MRCData(header, extendedheader)
    read!(newio, d.data)
    close(newio)
    map!(bswaptoh(header.machst), d.data, d.data)
    return d
end
function Base.read(fn::AbstractString, ::Type{T}; compress=:auto) where {T<:MRCData}
    return open(fn; read=true) do io
        return read(io, T; compress=compress)
    end
end

function Base.read(io::IO, ::Type{T}) where {T<:MRCHeader}
    names = fieldnames(T)
    types = T.types
    offsets = fieldoffsets(T)
    bytes = read!(io, Array{UInt8}(undef, HEADER_LENGTH))
    machst = bytes[213:216] # look ahead to machst
    bytes_ptr = pointer(bytes)
    vals = GC.@preserve bytes [
        bytestoentry(names[i], types[i], bytes_ptr + offsets[i]) for i in 1:length(offsets)
    ]
    map!(bswaptoh(machst), vals, vals)
    header = T(vals...)
    return header
end

function Base.read(io::IO, ::Type{T}; length=0, header=nothing) where {T<:MRCExtendedHeader}
    # TODO: use h.exttyp to identify extended header format and parse into a human-readable type
    if header !== nothing && length == 0
        length = header.nsymbt
    end
    data = read!(io, Array{UInt8}(undef, length))
    return T(data)
end

"""
    read_mmap(io::IO, ::Type{MRCData})
    read_mmap(path::AbstractString, ::Type{MRCData})

Read MRC file or stream, using a memory-mapped array to access the data.
"""
function read_mmap(io::IO, ::Type{MRCData})
    head = read(io, MRCHeader)
    exthead = read(io, MRCExtendedHeader; header=head)
    arraytype = Array{MRC.datatype(head),ndims(head)}
    data = Mmap.mmap(io, arraytype, size(head)[1:ndims(head)])
    return MRCData(head, exthead, data)
end
function read_mmap(path::AbstractString, T::Type{MRCData})
    return open(path, "r") do io
        return read_mmap(io, T)
    end
end

"""
    write(io::IO, ::MRCData; compress = :none)
    write(fn::AbstractString, ::MRCData; compress = :auto)

Write an instance of [`MRCData`](@ref) to an IO stream or new file.
Use `compress` to specify the compression with the following options:
- `:auto`: Infer one of the below compression types from the file extension
- `:gz`: GZ
- `:bz2`: BZ2
- `:xz`: XZ
- `:none`: no compression
"""
write(::Any, ::MRCData)
function Base.write(io::IO, d::MRCData; compress=:none, unit_vsize=4096)
    newio = compressstream(io, compress)
    h = header(d)
    sz = write(newio, h)
    sz += write(newio, extendedheader(d))
    T = datatype(h)
    data = parent(d)
    fswap = bswapfromh(h.machst)
    unit_vsize = div(unit_vsize, sizeof(T))
    vlen = length(data)
    vrem = vlen % unit_vsize
    if vrem != 0
        @inbounds @views sz += write(newio, fswap.(T.(data[1:vrem])))
    end
    for i in (vrem + 1):unit_vsize:vlen
        @inbounds @views sz += write(newio, fswap.(T.(data[i:(i + unit_vsize - 1)])))
    end
    close(newio)
    return sz
end
function Base.write(fn::AbstractString, object::T; compress=:auto, unit_vsize=4096) where {T<:Union{MRCData}}
    if compress == :auto
        compress = checkextension(fn)
    end
    return open(fn; write=true) do io
        return write(io, object; compress=compress, unit_vsize=unit_vsize)
    end
end

function Base.write(io::IO, header::MRCHeader)
    names = fieldnames(typeof(header))
    bytes = UInt8[]
    fswap = bswapfromh(header.machst)
    for name in names
        append!(bytes, entrytobytes(name, fswap(getfield(header, name))))
    end
    return write(io, bytes)
end

Base.write(io::IO, eh::MRCExtendedHeader) = write(io, eh.data)


text:
immutable ChineseRestaurantTableProcess <: DiscreteUnivariateDistribution
	n::Int # Number of customers
	alpha::Float64
end

typealias CRT ChineseRestaurantTableProcess

function Distributions.rand(d::ChineseRestaurantTableProcess)
	s = 0
	for t = 1:d.n
		if rand() < d.alpha / (t - 1 + d.alpha)
			s += 1
		end
	end
	return s
end


text:
"""
The Dist object contains information for convolving random variables.
* `lb`: lower bound of density
* `ub`: upper bound of density
* `n_points`: number of discrete points used in transform
* `x`: real values associated with densities
* `density`: densities resulting from Convolution
* `interp_pdf`: an object used to interpolate exact densities from x and density
```@julia
model = Normal(0,1) + Normal(0,1)
convolve!(model)

julia> pdf(model, 0.0)
0.28216441829987216

julia> pdf(Normal(0,sqrt(2)), 0.0)
0.28209479177387814
```
"""
mutable struct Convolution{T}
    model::T
    lb::Float64
    ub::Float64
    n_points::Int64
    x::Vector{Float64}
    density::Vector{Float64}
    interp_pdf::Spline1D
end

function Convolution(;model, lb=-3.0, ub=3.0, n_points = 2^9, x=fill(0.0, 0),
        density=fill(0.0, 0), interp_pdf=Spline1D(zeros(4), zeros(4)))
    Convolution(model, lb, ub, n_points, x, density, interp_pdf)
end

Broadcast.broadcastable(c::Convolution) = Ref(c)

function +(d::Distribution...)
    convolution = Convolution(model=d)
    return convolution
end

function analytic(D::Convolution)
    N = D.n_points;  model = D.model
    lb = D.lb;  ub = D.ub
    i = range(0, length=N); dx = (ub-lb)/N
    x = lb .+ i*dx; dt = 2*π/(N*dx)
    c = -N/2*dt; t = c .+ i*dt; f = fill(1im+1.0, N)
    for m in model
        @. f = f*cf(m, t)
    end
    X = @. exp(-(0+1im)*i*dt*lb)*f
    fft!(X)
    density = @. dt/(2*π)*exp(-(0+1im)*c*x)*X
    # Use abs()to haldle very small negative values, e.g. -2.50e-12
    D.density = abs.(real(density))
    D.x = x
    D.interp_pdf = Spline1D(D.x, D.density)
    return nothing
end

function numeric(D::Convolution)
   n_points = D.n_points;  model = D.model
   lb = D.lb;  ub = D.ub
   Ndist = length(model)
   stepsize = (ub-lb)/n_points
   x = range(lb, stop=ub, length=n_points)
   d = pdf.(model[1], x)
   mn = Int(round(1.5*n_points)) + 1
   st = Int(round(2*n_points))
   mx = n_points ÷ 2
   v = [fill(0.0, mx); d ;fill(0.0, mx)]
   Nv = length(v)
   d1 = pdf.(model[2], x)
   v1 = [fill(0.0, mx); d1; fill(0.0, mx)]
   m = fft(v).*fft(v1)
   cv = real(ifft(m))*Nv*stepsize/(2*n_points)
   cv = abs.(cv[[mn:st;1:mx]])
   cnt = 2
   while cnt < Ndist
       cnt += 1
       v = [fill(0.0, mx); cv; fill(0.0, mx)]
       d1 = pdf.(model[cnt], x)
       v1 = [fill(0.0, mx); d1 ;fill(0.0, mx)]
       m = fft(v).*fft(v1)
       cv = real(ifft(m))*Nv*stepsize/(2*n_points)
       cv = abs.(cv[[mn:st;1:mx]])
    end
    D.density = cv
    D.x = x
    D.interp_pdf = Spline1D(D.x, D.density)
    return nothing
end

"""
Convolves a set of random variables.
```@julia
model = Normal(0,1) + Normal(0,1)
convolve!(model)

julia> pdf(model, 0.0)
0.28216441829987216

julia> pdf(Normal(0,sqrt(2)), 0.0)
0.28209479177387814
```
"""
function convolve!(D::Convolution)
    if all(x->applicable(cf, (x,1.0)...), D.model)
        # Convolve with characteristic function if possible
        analytic(D)
    else
        numeric(D)
    end
    return nothing
end

 function rand(model::Convolution, N::Int64)
     rts = fill(0.0, N)
     dist = model.model
     for d in dist
         rts .+= rand(d, N)
     end
     return rts
 end

Gradient(d::Convolution, x::Float64) = d.interpGrad[x]
pdf(d::Convolution, x::Float64) = abs(d.interp_pdf(x))
logpdf(d::Convolution, x::Float64) = log(pdf(d, x))

"""
Normal approximation to the sum of Uniform rvs
* 'scaling' = (2/3): default scaling used for standard deviation σ
* 'args': a list of NamedTuples for component distributions, e.g. cr =(μ=.05,N=2),...

```@julia

julia> μ,σ = convolve_normal(cr=(μ=.05,N=2),vis=(μ=.085,N=1))
(0.185, 0.021278575558006173)
```
"""
function convolve_normal(;scaling=2/3, args...)
    # uniform factor
    fact1 = 1/12
    μ,σ = 0.0,0.0
    for (k,v) in pairs(args)
        μ += prod(v)
        σ += fact1.*v.N*(scaling*v.μ)^2
    end
    σ = sqrt(σ)
    return μ,σ
end


text:
using GlobalOptimizer, Test

@testset "Solve Ackley Function" begin
    # Define an objective function (using Ackley function here as an example)
    function ackley(x::AbstractVector{T}) where T<:Real
        arg1 = 0.0
        arg2 = 0.0
        dims = length(x)
        @inbounds for i in eachindex(x)
            arg1 += abs2(x[i])
            arg2 += cospi(2.0 * x[i])
        end
        arg1  = 0.2 * sqrt(arg1 / dims)
        arg2 /= dims
        return -20.0 * exp(-arg1) - exp(arg2) + ℯ + 20.0
    end

    nd  = 15
    lb  = ntuple(i -> -32.0, nd) # lower bounds
    ub  = ntuple(i -> +32.0, nd) # upper bounds
    obj = minimizer(nd)          # create an object for the optimization
    minimize!(obj, ackley, lb, ub, avgtimes=3)
    @test obj.xsol ≈ zeros(Float64, nd) atol=1e-12
end


text:
include("utils.jl")

build_so("010")

x = Array(Vector{Float64}, 3)
x[1] = [0.1, 0.2, 0.3]
x[2] = [0.4, 0.5, 0.6]
x[3] = [0.7, 0.8, 0.9]

y = ccall(
    (:foo, "binary/lib010.so"),
    Void,
    (Ptr{Ptr{Cdouble}}, Cint, Cint),
    x, length(x), length(x[1])
)

@assert x[1] == [0.1, 0.2, 0.3] + 1.0
@assert x[2] == [0.4, 0.5, 0.6] + 1.0
@assert x[3] == [0.7, 0.8, 0.9] + 1.0

@assert y === nothing


text:
CONSTRAINT_TYPES = Union{T, Vector{T}} where T <: SamplingConstraint

"""
    ConstrainedValueResampling{N_DATASETS}

Indicates that resampling should be done with constraints on the furnishing distributions/populations.

## Fields

- `constraints`. The constraints for the datasets. The constraints are represented as a tuple of length
    `N_DATASETS`, where the `i`-th tuple element contains the constraints for that dataset.
    Constraints for each dataset must be supplied as either a single sampling constraint,
    or as a vector of sampling constraints with length matching the length of the dataset
    (`Union{SamplingConstraint, Vector{<:SamplingConstraint}}}`). For example, if the `i`-th dataset
    contains 352 observations, then `constraints[i]` must be either a single
    sampling constraint (e.g. `TruncateStd(1.1)`) or a vector of 352 different sampling constraints
    (e.g. `[TruncateStd(1.0 + rand()) for i = 1:352]`).
- `n::Int`. The number of draws.

## Example

Assume we have three collections of uncertain values of, each of length `L = 50`. These should be
resampled `250` times. Before resampling, however, the distributions/populations furnishing
the uncertain values should be truncated:

- For the first collection, truncate each value at `1.5` times its standard deviation
    around its mean. This could simulate measurement errors from an instrument
    that yields stable measurements whose errors are normally distributed, but for
    which we are not interested in outliers or values beyond `1.5` standard devations
    for our analyses.
- For the second collection, truncate each value at the `80`th percentile range.
    This could simulate measurement errors from an instrument that yields stable measurements,
    whose errors are not normally distributed, so that confidence intervals are better
    to use than standard deviations. In this case, we're not interested in outliers,
    and therefore exclude values smaller than the `10`th percentile and larger than the
    `90`th percentile of the data.
- For the third collection, truncate the `i`-th value at an fraction of its standard
    deviation around the mean slightly larger than at the `i-1`-th value, so that
    the standard deviation ranges from `0.5` to `0.5 + L/100`. This could simulate,
    for example, an instrument whose measurement error increases over time.

```julia
L = 50
constraints_d1 = TruncateStd(1.5)
constraints_d2 = TruncateQuantiles(0.1, 0.9)
constraints_d3 = [TruncateStd(0.5 + i/100) for i = 1:L]
```
"""
struct ConstrainedValueResampling{N} <: AbstractUncertainDataResampling
    constraints::Tuple{Vararg{CONSTRAINT_TYPES, N}}
    n::Int
end

# It's tedious for the user to always provide a tuple, so unpack a variable
# number of constraints into a tuple, then call the original constructor.
function ConstrainedValueResampling(constraints::CONSTRAINT_TYPES...) where N
    ConstrainedValueResampling((constraints...,), 1)
end

function ConstrainedValueResampling(n::Int, constraints::CONSTRAINT_TYPES...) where N
    ConstrainedValueResampling((constraints...,), n)
end

Broadcast.broadcastable(c::ConstrainedValueResampling) = Ref(c)

Base.length(c::ConstrainedValueResampling) = length(c.constraints)
Base.firstindex(c::ConstrainedValueResampling) = 1
Base.lastindex(c::ConstrainedValueResampling) = length(c)
Base.getindex(c::ConstrainedValueResampling, i) = c.constraints[i]
Base.iterate(c::ConstrainedValueResampling, state = 1) = iterate(c.constraints, state)
Base.eachindex(c::ConstrainedValueResampling) = Base.OneTo(length(c.constraints))

function Base.show(io::IO, constraints::ConstrainedValueResampling{N_DATASETS}) where {N_DATASETS}
    s = "$(typeof(constraints)) for $N_DATASETS set(s) of value constraints, where n=$(constraints.n)"
    println(io, s)
end

export ConstrainedValueResampling


text:
"""

A package for taking groups of sequences and finding a sequence
that is representitive of the group via the levenshtein distance.

Note that the sequence returned may not be a meaningful sequence, For example, it
would not return a 'correcnt' DNA sequence that would be useful.

"""

module SequenceReduction

export basic_lev, lev_update, lev_trace, lev_distance, lev_create

include("basic-lev-reduction.jl")

end # module


text:
module MLDataStructs

export MLData

#-----------------------------------------------------------------------------# Task
abstract type MLTask end

for T in [:Regression, :Classification, :Clustering]
    @eval begin
        struct $T <: MLTask end
        Base.show(io::IO, t::$T) = print(io, replace(string($T), r".*\." => ""))
    end
end
# struct Regression <: MLTask end
# struct Classification <: MLTask end
# struct Clustering <: MLTask end

#-----------------------------------------------------------------------------# MLData
const AN = Union{AbstractArray,  Nothing}

struct MLData{X <: AN, Y <: AN, W <: AN, T <: MLTask}
    x::X
    y::Y
    w::W
    task::T
    function MLData(x::X, y::Y, w::W, task::T) where {X,Y,W,T}
        nx, ny, nw = _nobs(x), _nobs(y), _nobs(w)
        n = max(nx, ny, nw)
        nx == ny == 0 && error("At least one of x or y must be provided.")
        all(x -> x==n || x==0, (nx, ny, nw)) || error("Components have different nobs.")
        new{X,Y,W,T}(x, y, w, task)
    end
end

MLData(;x=nothing, y=nothing, w=nothing, task=determine_task(x, y)) = MLData(x, y, w, task)

_nobs(x::AbstractArray) = size(x, 1)
_nobs(x::Nothing) = 0

function Base.show(io::IO, o::MLData)
    s = (!isnothing(o.x) && !isnothing(o.y)) ? "Supervised" : "Unsupervised"
    w = !isnothing(o.w) ? "Weighted" : "Un-weighted"
    print(io, "MLData | $s ($(o.task)) | $w\n ▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬▬")
    o.x != nothing && print(io, "\n  > x: ", summary(o.x))
    o.y != nothing && print(io, "\n  > y: ", summary(o.y))
    o.w != nothing && print(io, "\n  > w: ", summary(o.w))
end

xy(o::MLData) = o.x, o.y
xyw(o::MLData) = o.x, o.y, o.w

nobs(o::MLData{<:AbstractArray}) = size(o.x, 1)
nobs(o::MLData{Nothing, <:AbstractArray}) = size(o.y, 1)

npredictors(o::MLData{<:AbstractArray}) = size(o.x, 2)
npredictors(o::MLData{Nothing}) = 0

#-----------------------------------------------------------------------------# determine_task
determine_task(x::AbstractArray, y::Nothing) = Clustering()
function determine_task(x::AbstractArray, y::AbstractArray)
    length(unique(y)) > .1 * _nobs(y) ? Regression() : Classification()
end

end # module


text:
using DataDeps

register(DataDep(
  "thyroid-disease-ann-thyroid",
  "http://archive.ics.uci.edu/ml/datasets/Thyroid+Disease",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/thyroid-disease/ann-train.data",
  "3da53a156bda36cb0c97e9f4b6b111c9226c54c4aa00230de5604b787c47e3a6",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=22,
        feature_indices=1:21,
        delim=' ',
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
struct CodonLookupTable{T}
    table::T
end

const SingleCodonLookup{T} = CodonLookupTable{Vector{T}}
const PairwiseCodonLookup{T} = CodonLookupTable{PairwiseListMatrix{T, true, Vector{T}}}

function SingleCodonLookup{T}(f::Function) where T
    v = SingleCodonLookup{T}(Vector{T}(64))
    cdn = kmer"AAA"
    @inbounds while cdn < kmer"TTT"
        v[cdn] = f(cdn)
        cdn += 1
    end
    v[cdn] = f(cdn)
    return v
end

function PairwiseCodonLookup{T}(f::Function) where T
    v = PairwiseCodonLookup{T}(PairwiseListMatrix(Vector{T}(2080), true))
    ci = kmer"AAA"
    while ci < kmer"TTT"
        cj = ci
        @inbounds while cj < kmer"TTT"
            v[ci, cj] = f(ci, cj)
            cj += 1
        end
        v[ci, cj] = f(ci, cj)
        ci += 1
    end
    return v
end

@inline function Base.getindex(tbl::SingleCodonLookup, cdn::C) where C <: Codon
    return tbl.table[UInt64(cdn) + 1]
end
@inline function Base.setindex!(tbl::SingleCodonLookup{T}, val::T, cdn::C) where {C <: Codon, T}
    return tbl.table[UInt64(cdn) + 1] = val
end

@inline function Base.getindex(tbl::PairwiseCodonLookup, i::C, j::C) where C <: Codon
    return tbl.table[UInt64(i) + 1, UInt64(j) + 1]
end
@inline function Base.setindex!(tbl::PairwiseCodonLookup{T}, val::T, i::C, j::C) where {C <: Codon, T}
    return tbl.table[UInt64(i) + 1, UInt64(j) + 1] = val
end

## To Sort out!

function make_edge_reference(::Type{C}, code::GeneticCode) where C <: Codon
    edges = Vector{Tuple{C, C, Int, Int}}(2016)
    k = 1
    for i in 0x00:0x3F
        for j in (i + 1):0x3F
            x = C(UInt64(i))
            y = C(UInt64(j))
            DS, DN = DS_DN_enumerator(ShortestPath, x, y, code)
            edges[k] = (x, y, Integer(DS), Integer(DN))
            k += 1
        end
    end
    return edges
end

const DSDN_RANK_LOOKUP = [0 2 5 9;
                          1 4 8 0;
                          3 7 0 0;
                          6 0 0 0;]

@inline function rankof(DS::Integer, DN::Integer)
    @inbounds return DSDN_RANK_LOOKUP[DS + 1, DN + 1]
end

struct CodonGraphReference{C <: Codon}
    edges::Vector{Tuple{C, C, Int, Int}}
    permutation::Vector{Int}
    order::Vector{Int}
end

function CodonGraphReference{C}(code::GeneticCode) where C <: Codon
    edgetable = make_edge_reference(C, code)
    edgetable_perm = sortperm(edgetable, by = x -> rankof(x[3], x[4]))
    edgetable_order = zeros(edgetable_perm)
    for i in eachindex(edgetable_perm)
        perm = edgetable_perm[i]
        edgetable_order[perm] = i
    end
    return CodonGraphReference{C}(edgetable, edgetable_perm, edgetable_order)
end

@inline function ij_to_k(i::Integer, j::Integer)
    i, j = ifelse(i > j, (j, i), (i, j))
    nelements = 64
    x = nelements - i
    return div(nelements * (nelements - 1) - (x * (x - 1)), 2) - nelements + j
end

@inline function codons_to_k(i::T, j::T) where T <: Codon
    return ij_to_k(UInt64(i) + 1, UInt64(j) + 1)
end

@inline function lookup_edge(table::CodonGraphReference{T}, i::T, j::T) where T <: Codon
    @inbounds return table.edges[codons_to_k(i, j)]
end

@inline function lookup_edge(table::CodonGraphReference, i::Int)
    @inbounds return table.edges[table.permutation[i]]
end

@inline function lookup_perm(table::CodonGraphReference{T}, i::T, j::T) where T <: Codon
    @inbounds return table.permutation[codons_to_k(i, j)]
end

@inline function lookup_order(table::CodonGraphReference{T}, i::T, j::T) where T <: Codon
    @inbounds return table.order[codons_to_k(i, j)]
end

const DEFAULT_CODON_GRAPH_REFERENCE = CodonGraphReference{Codon{DNA}}(DEFAULT_TRANS)


text:
macro declare_file_product(product_name)
    path_name = Symbol(string(product_name, "_path"))
    return esc(quote
        # These will be filled in by init_file_product()
        $(path_name) = ""
        $(product_name) = ""
    end)
end

macro init_file_product(product_name, product_path)
    path_name = Symbol(string(product_name, "_path"))
    return esc(quote
        # FileProducts are very simple, and we maintain the `_path` suffix version for consistency
        global $(path_name) = joinpath(artifact_dir, $(product_path))
        global $(product_name) = $(path_name)
    end)
end


text:
struct RMSPropIterable{T, F, S, R}
    w0::T
    f::F
    stepsize::S
    alpha::R
    epsilon::R
end

RMSPropIterable(w0, f; stepsize=1e-3, alpha=0.99, epsilon=1e-8) = RMSPropIterable(w0, f, to_iterator(stepsize), alpha, epsilon)

Base.IteratorSize(::Type{<:RMSPropIterable}) = Base.IsInfinite()

mutable struct RMSPropState{T, G, S}
    w::T
    f_w
    grad_f_w::G
    stepsize_iterator::S
    mean_squared_grad::T
end

function Base.iterate(iter::RMSPropIterable)
    w = copy(iter.w0)
    grad_f_w, f_w = gradient(iter.f, w)
    state = RMSPropState(w, f_w, grad_f_w, Iterators.Stateful(iter.stepsize), zero(w))
    return IterationOutput(state.w, state.f_w, state.grad_f_w), state
end

function Base.iterate(iter::RMSPropIterable, state::RMSPropState)
    state.mean_squared_grad .= (1 - iter.alpha) .* (state.grad_f_w .^ 2) .+ iter.alpha .* state.mean_squared_grad
    stepsize = popfirst!(state.stepsize_iterator)
    state.w .-= stepsize .* state.grad_f_w ./ (sqrt.(state.mean_squared_grad) .+ iter.epsilon)
    state.grad_f_w, state.f_w = gradient(iter.f, state.w)
    return IterationOutput(state.w, state.f_w, state.grad_f_w), state
end

struct RMSProp
    kwargs
    RMSProp(; kwargs...) = new(kwargs)
end

(alg::RMSProp)(args...) = RMSPropIterable(args...; alg.kwargs...)