text:
module EnsemblePTSampler

using Distributed

using Statistics

function propose(ps, qs)
    nd, nw, nt = size(ps)

    ps_out = zeros(nd, nw, nt)
    zs = zeros(nw, nt)

    for i in 1:nw
        for j in 1:nt
            zs[i,j] = exp(log(0.5) + (log(2.0)-log(0.5))*rand())
            ii = rand(1:nw)
            ps_out[:,i,j] = qs[:,ii,j] + zs[i,j]*(ps[:,i,j] - qs[:,ii,j])
        end
    end

    ps_out, zs
end

function lnlike_lnprior(x, ll, lp)
    p = lp(x)
    if p == -Inf
        (-Inf, -Inf)
    else
        (ll(x), p)
    end
end

function lnlikes_lnpriors(ps, ll, lp)
    nd, nw, nt = size(ps)
    if length(workers()) > 1
        parr = Array{Float64,1}[]
        for j in 1:nt
            for i in 1:nw
                push!(parr, ps[:,i,j])
            end
        end
        lls_lps = pmap(x -> lnlike_lnprior(x, ll, lp), parr, batch_size=div(nw*nt,(4*length(workers()))))
        lls = reshape(Float64[l for (l,p) in lls_lps], (nw, nt))
        lps = reshape(Float64[p for (l,p) in lls_lps], (nw, nt))
        (lls, lps)
    else
        lls = zeros(nw, nt)
        lps = zeros(nw, nt)
        for i in 1:nw
            for j in 1:nt
                lps[i,j] = lp(ps[:,i,j])
                if lps[i,j] == -Inf
                    lls[i,j] = -Inf
                else
                    lls[i,j] = ll(ps[:,i,j])
                end
            end
        end
        (lls, lps)
    end
end

function update_half(ps, llps, lpps, qs, ll, lp, betas)
    nd, nw, nt = size(ps)

    prop_ps, zs = propose(ps, qs)
    prop_llps, prop_lpps = lnlikes_lnpriors(prop_ps, ll, lp)

    new_ps = zeros(size(ps)...)
    new_llps = zeros(size(llps)...)
    new_lpps = zeros(size(lpps)...)

    for i in 1:nw
        for j in 1:nt
            lpacc = betas[j]*(prop_llps[i,j] - llps[i,j]) + prop_lpps[i,j] - lpps[i,j] + nd*log(zs[i,j])
            if log(rand()) < lpacc
                new_ps[:,i,j] = prop_ps[:,i,j]
                new_llps[i,j] = prop_llps[i,j]
                new_lpps[i,j] = prop_lpps[i,j]
            else
                new_ps[:,i,j] = ps[:,i,j]
                new_llps[i,j] = llps[i,j]
                new_lpps[i,j] = lpps[i,j]
            end
        end
    end

    (new_ps, new_llps, new_lpps)
end

"""
    mcmc_step(pts, lnlikes, lnpriors, loglike, logprior, betas)

Returns `(new_pts, new_lnlikes, new_lnpriors)` from a single MCMC step operating
on the given ensemble.

# Arguments
* `pts` a `(ndim, nwalkers, ntemps)` array of the current state of the chain.
* `lnlikes` a `(nwalkers, ntemps)` array giving the current log-likelihoods.
* `lnpriors` a `(nwalkers, ntemps)` array giving the current log-priors.
* `loglike` a function that gives the log-likelihood for parameter values.
* `logprior` a function that gives the log-prior for parameter values.
* `betas` a `(ntemps,)` array with the current inverse temperatures.

"""
function mcmc_step(pts, lnlikes, lnpriors, loglike, logprior, betas)
    nd, nw, nt = size(pts)

    @assert nw%2==0 "must have even number of walkers"

    ihalf = div(nw, 2)

    ps = pts[:,1:ihalf,:]
    llps = lnlikes[1:ihalf,:]
    lpps = lnpriors[1:ihalf,:]

    qs = pts[:,ihalf+1:end,:]
    llqs = lnlikes[ihalf+1:end,:]
    lpqs = lnpriors[ihalf+1:end,:]

    new_ps, new_llps, new_lpps = update_half(ps, llps, lpps, qs, loglike, logprior, betas)
    new_qs, new_llqs, new_lpqs = update_half(qs, llqs, lpqs, new_ps, loglike, logprior, betas)

    new_pts = cat(new_ps, new_qs, dims=2)
    new_lnlikes = cat(new_llps, new_llqs, dims=1)
    new_lnpriors = cat(new_lpps, new_lpqs, dims=1)
    (new_pts, new_lnlikes, new_lnpriors)
end

"""
    tswap_step(pts, lnlikes, lnpriors, betas)

Implements the parallel-tempering swap step between ensembles of different
temperatures.  Arguments are as for `mcmc_step`.

"""
function tswap_step(pts, lnlikes, lnpriors, betas)
    nd, nw, nt = size(pts)

    new_pts = copy(pts)
    new_lnlikes = copy(lnlikes)
    new_lnpriors = copy(lnpriors)
    tswaps = zeros(Int, nt-1)

    for k in nt:-1:2
        for j in 1:nw
            jj = rand(1:nw)

            llhigh = new_lnlikes[j,k]
            lllow = new_lnlikes[jj, k-1]

            bhigh = betas[k]
            blow = betas[k-1]
            db = bhigh-blow

            lpacc = db*(lllow - llhigh)

            if log(rand()) < lpacc
                phigh = new_pts[:,j,k]
                lphigh = new_lnpriors[j,k]

                new_pts[:,j,k] = new_pts[:,jj,k-1]
                new_lnlikes[j,k] = new_lnlikes[jj,k-1]
                new_lnpriors[j,k] = new_lnpriors[jj,k-1]

                new_pts[:,jj,k-1] = phigh
                new_lnlikes[jj,k-1] = llhigh
                new_lnpriors[jj,k-1] = lphigh

                tswaps[k-1] += 1
            end # If not accepted, do nothing
        end
    end

    new_pts, new_lnlikes, new_lnpriors, tswaps
end

"""
    tevolve(swapfraction, betas, tc)

Returns `new_betas`, which are evolved according to a rule to equalise
temperature swap rates over a timescale `tc`.

The rule is formulated in terms of a quantity `x[i]`, which is the logit
transform of `beta[i]` relative to its neighbours (recall that `beta` is the
*inverse* temperature, so decreases as `i` increases):

    x[i] = log(beta[i] - beta[i+1]) - log(beta[i-1] - beta[i])

`x` maps the range between `beta[i+1]` and `beta[i-1]` to `-Inf` to `Inf`.  We
evolve `x` via

    x_new[i] = x[i] + 2.0*(swapfraction[i] - swapfraction[i-1])/(tc*(swapfraction[i] + swapfraction[i-1]))

where `swapfraction[i]` measures the fraction of accepted swaps between
temperature `i` and temperature `i+1` (i.e. between chain `i` and the
next-highest temperature chain) and `tc` is a "time constant" that controls the
(exponential) rate of convergence of `x`.

This is similar to the evolution rule in [Vousden, Farr, & Mandel
(2016)](https://ui.adsabs.harvard.edu/#abs/2016MNRAS.455.1919V/abstract).

"""
function tevolve(swapfraction, betas, tc)
    new_betas = copy(betas)

    for i in 2:size(betas, 1)-1
        if swapfraction[i] == 0 && swapfraction[i-1] == 0
            # Do nothing---we don't know which way to move!
        else
            x = log(betas[i] - betas[i+1]) - log(betas[i-1] - betas[i])
            xnew = x + 2.0*(swapfraction[i] - swapfraction[i-1])/(tc*(swapfraction[i] + swapfraction[i-1]))
            exnew = exp(xnew)
            new_betas[i] = (new_betas[i+1] + exnew*new_betas[i-1])/(1.0 + exnew)
        end
    end

    new_betas
end

"""
    run_mcmc(pts, loglike, logprior, betas, nstep; [thin=1])

Run a MCMC simulation of length `nstep`.

The simulation will consist of a "burnin" phase and a "sampling" phase, each
involving half `nstep` iterations.  During "burnin", the temperatures of the
chain will adapt with a time constant of 1/10 of the burnin length.  During
sampling the temperatures are fixed to preserve detailed balance.

The function will return `(chain, chainloglike, chainlogprior, new_betas)` where
each returned chain value has an extra dimension appended counting steps of the
chain (so `chain` is of shape `(ndim, nwalkers, ntemp, nstep)`, for example).

# Arguments
* `pts` a `(ndim, nwalkers, ntemp)` array with the current state.
* `loglike` a function giving the log-likelihood for parameters.
* `logprior` a function giving the log-prior for parameters.
* `betas` a `(ntemp,)` array of inverse temperatures.
* `thin=1` save the output every `thin` steps.

"""
function run_mcmc(pts, loglike, logprior, betas, nstep; thin=1)
    nd, nw, nt = size(pts)
    lnlikes, lnpriors = lnlikes_lnpriors(pts, loglike, logprior)

    nhalf = div(nstep, 2)

    for i in 1:nhalf
        pts, lnlikes, lnpriors = mcmc_step(pts, lnlikes, lnpriors, loglike, logprior, betas)
        pts, lnlikes, lnpriors, ntswap = tswap_step(pts, lnlikes, lnpriors, betas)
        swapfraction = ntswap / nw
        tc = 10.0 + i/10.0 # tc > 10 always, but asymptotes to "current length / 10"
        betas = tevolve(swapfraction, betas, tc)
    end

    nsave = div(nhalf, thin)
    chain = zeros(nd, nw, nt, nsave)
    chainlnlike = zeros(nw, nt, nsave)
    chainlnprior = zeros(nw, nt, nsave)
    isave = 1

    for i in 1:nhalf
        pts, lnlikes, lnpriors = mcmc_step(pts, lnlikes, lnpriors, loglike, logprior, betas)
        pts, lnlikes, lnpriors, _ = tswap_step(pts, lnlikes, lnpriors, betas)

        if i % thin == 0
            chain[:,:,:,isave] = pts
            chainlnlike[:,:,isave] = lnlikes
            chainlnprior[:,:,isave] = lnpriors
            isave += 1
        end
    end

    chain, chainlnlike, chainlnprior, betas
end

function trapz(ys, xs)
    sum(0.5*diff(xs).*(ys[1:end-1] + ys[2:end]))
end

"""
    lnZ(lnlikes, betas)

Returns `(ln_Z, delta_ln_Z)` estimating the log-evidence and uncertainty using
thermodynamic integration.

# Arguments
* `lnlikes` an `(nwalkers, ntemps, nsteps)` array of log-likelihoods.
* `betas` an `(ntemps,)` array of the corresponding beta values.

"""
function lnZ(lnlikes, betas)
    mean_lnlikes = vec(mean(lnlikes, dims=(1, 3)))

    bs_thin = [betas[1]]
    lnls_thin = [mean_lnlikes[1]]
    for i in 2:2:size(betas,1)-1
        push!(bs_thin, betas[i])
        push!(lnls_thin, mean_lnlikes[i])
    end
    push!(bs_thin, betas[end])
    push!(lnls_thin, mean_lnlikes[end])

    lnZ = -trapz(mean_lnlikes, betas)
    lnZ2 = -trapz(lnls_thin, bs_thin)

    lnZ, abs(lnZ-lnZ2)
end

end


text:
using Ladder
using Base.Test

ext: jl


text:
using JuMP
N = 6;
m = Model();				# Use default solver
@defVar(m, x[1:N], Bin);		# Define array variable to hold results
profit  =  [ 5, 3, 2, 7, 4, 4 ];	# Profit vector of size N
weight = [ 2, 8, 4, 2, 5, 6 ];		# Weights vector of size

maxcap = 12;
@setObjective(m, Max, dot(profit, x));
@addConstraint(m, dot(weight, x) <= maxcap);

status = solve(m);			# Solve problem using MIP solver
@printf "Value of the objective is %.1f" getObjectiveValue(m);

println("Solution is:")
for i = 1:N
    print("x[$i] = ", getValue(x[i]));
    println(", p[$i]/w[$i] = ", profit[i]/weight[i]);
end

max_stars_repo_path: Module 3/Chapter06/knapsack.jl

max_issues_repo_path: Module 3/Chapter06/knapsack.jl

max_forks_repo_path: Module 3/Chapter06/knapsack.jl


text:
#  MIT License
#
#  Copyright (c) 2020 Michael J Simms. All rights reserved.
#
#  Permission is hereby granted, free of charge, to any person obtaining a copy
#  of this software and associated documentation files (the "Software"), to deal
#  in the Software without restriction, including without limitation the rights
#  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#  copies of the Software, and to permit persons to whom the Software is
#  furnished to do so, subject to the following conditions:
#
#  The above copyright notice and this permission notice shall be included in all
#  copies or substantial portions of the Software.
#
#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
#  SOFTWARE.

include("Distance.jl")
include("Peaks.jl")
include("Powers.jl")
include("Signals.jl")

using Pkg
Pkg.add("ArgParse")
using ArgParse
Pkg.add("CSV")
using CSV

function read_accelerometer_csv(filename::String)
    data = []
    data = CSV.read(filename)
    ts = data[1]
    x = data[2]
    y = data[3]
    z = data[4]
    ts, x, y, z
end

function distance_tests()
	println("Distance Tests:")
	println("---------------")

	distance = Distance.hamming_distance("1011101", "1001001")
	println("Hamming Distance: ", distance)
	@assert distance == 2

	distance = Distance.levenshtein_distance("foo", "foobar")
	println("Levenshtein Distance: ", distance)
	@assert distance == 3
end

function signals_tests()
	println("Signals Tests:");
	println("--------------");

	indata = [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0 ]
    windowsize = 2
    outdata = Signals.smooth(indata, windowsize)

    println(outdata)
end

function power_tests()
	println("Power Tests:")
	println("------------")

	nearest = Powers.nearest_power_of_2(63)
	println("Nearest power of 2 for 63 is ", nearest)
	@assert nearest == 64
end

function peak_finding_tests(data)
	println("Peak Finding Tests:")
	println("-------------------")

    peaks = Peaks.find_peaks_over_stddev(data)
    println(peaks)
end

# Parses the command line arguments
function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--csv"
            help = "A CSV file with accelerometer data"
            arg_type = String
            default = "../data/10_pullups.csv"
    end

    return parse_args(s)
end

# Run all the unit tests
parsed_args = parse_commandline()
ts, x, y, z = read_accelerometer_csv(parsed_args["csv"])
distance_tests()
println("")
signals_tests()
println("")
power_tests()
println("")
peak_finding_tests(x)
println("")

max_stars_repo_name: msimms/LibMath

max_issues_repo_name: msimms/LibMath

max_forks_repo_name: msimms/LibMath


text:
module VDPTag2

# package code goes here
using POMDPs
using StaticArrays
using Parameters
using Plots
using Distributions
using POMDPModelTools
using POMDPPolicies
using POMDPSimulators
using ParticleFilters
using Random
using LinearAlgebra

const Vec2 = SVector{2, Float64}
const Vec8 = SVector{8, Float64}

# importall POMDPs
import Base: rand, eltype, convert
import MCTS: next_action, n_children
import ParticleFilters: obs_weight
import POMDPs: actions, isterminal

export
    TagState,
    TagAction,
    VDPTagMDP,
    VDPTagPOMDP,
    Vec2,

    AODiscreteVDPTagPOMDP,
    ADiscreteVDPTagPOMDP,
    cproblem,

    convert_s,
    convert_a,
    convert_o,

    ToNextML,
    ToNextMLSolver,
    NextMLFirst,
    ManageUncertainty,
    CardinalBarriers,
    mdp
    isterminal

struct TagState
    agent::Vec2
    target::Vec2
end

struct TagAction
    look::Bool
    angle::Float64
end

@with_kw struct VDPTagMDP{B} <: MDP{TagState, Float64}
    mu::Float64          = 2.0
    agent_speed::Float64 = 1.0
    dt::Float64          = 0.1
    step_size::Float64   = 0.5
    tag_radius::Float64  = 0.1
    tag_reward::Float64  = 100.0
    step_cost::Float64   = 1.0
    pos_std::Float64     = 0.05
    barriers::B          = nothing
    tag_terminate::Bool  = true
    discount::Float64    = 0.95
end

@with_kw struct VDPTagPOMDP{B} <: POMDP{TagState, TagAction, Vec8}
    mdp::VDPTagMDP{B}           = VDPTagMDP()
    meas_cost::Float64          = 5.0
    active_meas_std::Float64    = 0.1
    meas_std::Float64           = 5.0
end

const VDPTagProblem = Union{VDPTagMDP,VDPTagPOMDP}
mdp(p::VDPTagMDP) = p::VDPTagMDP
mdp(p::VDPTagPOMDP) = p.mdp::VDPTagMDP

function next_ml_target(p::VDPTagMDP, pos::Vec2)
    steps = round(Int, p.step_size/p.dt)
    for i in 1:steps
        pos = rk4step(p, pos)
    end
    return pos
end
next_ml_target(p::VDPTagMDP, pos::AbstractVector) = next_ml_target(p, convert(Vec2, pos))

function POMDPs.transition(pp::VDPTagProblem, s::TagState, a::Float64)
    ImplicitDistribution(pp, s, a) do pp, s, a, rng
        p = mdp(pp)
        targ = next_ml_target(p, s.target) + p.pos_std*SVector(randn(rng), randn(rng))
        agent = barrier_stop(p.barriers, s.agent, p.agent_speed*p.step_size*SVector(cos(a), sin(a)))
        return TagState(agent, targ)
    end
end

POMDPs.reward(pp::VDPTagPOMDP, s::TagState, a::Float64, sp::TagState) = reward(pp.mdp, s, a, sp)
function POMDPs.reward(p::VDPTagMDP, s::TagState, a::Float64, sp::TagState)
    if norm(sp.agent-sp.target) < p.tag_radius
        return p.tag_reward
    else
        return -p.step_cost
    end
end

POMDPs.discount(pp::VDPTagProblem) = mdp(pp).discount
POMDPs.isterminal(pp::VDPTagProblem, s::TagState) = mdp(pp).tag_terminate && norm(s.agent-s.target) < mdp(pp).tag_radius

struct AngleSpace end
rand(rng::AbstractRNG, ::AngleSpace) = 2*pi*rand(rng)
POMDPs.actions(::VDPTagMDP) = AngleSpace()

POMDPs.transition(p::VDPTagPOMDP, s::TagState, a::TagAction) = transition(p, s, a.angle)

struct POVDPTagActionSpace end
rand(rng::AbstractRNG, ::POVDPTagActionSpace) = TagAction(rand(rng, Bool), 2*pi*rand(rng))
POMDPs.actions(::VDPTagPOMDP) = POVDPTagActionSpace()

function POMDPs.reward(p::VDPTagPOMDP, s::TagState, a::TagAction, sp::TagState)
    return reward(mdp(p), s, a.angle, sp) - a.look*p.meas_cost
end

#=
Beam | covers (deg)
-------------------
1    | (0,45]
2    | (45,90]
etc.
=#

struct BeamDist
    abeam::Int
    an::Normal{Float64}
    n::Normal{Float64}
end

function rand(rng::AbstractRNG, d::BeamDist)
    o = MVector{8, Float64}(undef)
    for i in 1:length(o)
        if i == d.abeam
            o[i] = rand(rng, d.an)
        else
            o[i] = rand(rng, d.n)
        end
    end
    return SVector(o)
end

function POMDPs.pdf(d::BeamDist, o::Vec8)
    p = 1.0
    for i in 1:length(o)
        if i == d.abeam
            p *= POMDPs.pdf(d.an, o[i])
        else
            p *= POMDPs.pdf(d.n, o[i])
        end
    end
    return p
end

function active_beam(rel_pos::Vec2)
    angle = atan(rel_pos[2], rel_pos[1])
    while angle <= 0.0
        angle += 2*pi
    end
    bm = ceil(Int, 8*angle/(2*pi))
    return clamp(bm, 1, 8)
end

function POMDPs.observation(p::VDPTagPOMDP, a::TagAction, sp::TagState)
    rel_pos = sp.target - sp.agent
    dist = norm(rel_pos)
    abeam = active_beam(rel_pos)
    if a.look
        an = Normal(dist, p.active_meas_std)
    else
        an = Normal(dist, p.meas_std)
    end
    n = Normal(1.0, p.meas_std)
    BeamDist(abeam, an, n)
end

POMDPs.observation(p::VDPTagPOMDP, a::Float64, sp::TagState) = observation(p, TagAction(false, a), sp)

include("rk4.jl")
include("barriers.jl")
include("initial.jl")
include("discretized.jl")
include("visualization.jl")
include("heuristics.jl")

function POMDPModelTools.gbmdp_handle_terminal(pomdp::VDPTagPOMDP, updater::Updater, b::ParticleCollection, s, a, rng)
    return ParticleCollection([s])
end

end # module


text:
import SortingAlgorithms: uint_mapping, RadixSort
import Base.Ordering

function load_uint(::Type{T}, s::String, skipbytes) where T <: Unsigned
    ss = sizeof(s)
    if skipbytes > ss
        return T(0)
    elseif skipbytes + sizeof(T) > ss
        s1 = ((s |> pointer |> Ptr{T}) + skipbytes) |> unsafe_load |> ntoh
        extra_bits_to_shift = 8(skipbytes + sizeof(T) - ss)
        return (s1 >> extra_bits_to_shift) << extra_bits_to_shift
    else
        return s1 = ((s |> pointer |> Ptr{T}) + skipbytes) |> unsafe_load |> ntoh
    end
end

load_uint(s::String, skipbytes) = load_uint(UInt, s, skipbytes)

function load_uint(sptr::Ptr{T}, ss, skipbytes) where T <: Unsigned
    if skipbytes > ss
        return T(0)
    elseif skipbytes + sizeof(T) > ss
        s1 = ((sptr |> Ptr{T}) + skipbytes) |> unsafe_load |> ntoh
        extra_bits_to_shift = 8(skipbytes + sizeof(T) - ss)
        return (s1 >> extra_bits_to_shift) << extra_bits_to_shift
    else
        return s1 = ((sptr |> Ptr{T}) + skipbytes) |> unsafe_load |> ntoh
    end
end

# function load_uint_shift(s::String, skipbytes)
#     s1 = ((((s |> pointer |> Ptr{UInt64}) + skipbytes) |> unsafe_load) << 32) >> 32
#     Base.zext_int(UInt64, s1) |> ntoh
# end

# timing is similar
# @time load_uint.(svec, 4);
# @time load_uint_shift.(svec, 4);

# @time radixsort(svec); #20
# @time sort!(svec, by = x -> load_uint(x, 4)); #61

# @time svec_sorted = radixsort(svec);
# print(issorted(svec_sorted))

"""
    fsortperm(svec)

Faster sortperm for string vectors
"""
function fsortperm(svec::AbstractVector{String}, ::Type{T} = UInt) where T<:Unsigned
    strlen = maximum(sizeof, svec)
    strlen = max(strlen-sizeof(T), 0)

    l = length(svec)

    pairs = Vector{Tuple{T, UInt32}}(undef, l)
    for (i, svec1) = enumerate(svec)
        pairs[i] = tuple(load_uint(T, svec1, strlen), i)
    end

    sort!(pairs, by=x->x[1], alg = RadixSort)

    while strlen > 0
        strlen = max(strlen-sizeof(T), 0)

        for (i, p) in enumerate(pairs)
            p2 = p[2]
            pairs[i] = tuple(load_uint(T, svec[p2], strlen), p2)
        end

        sort!(pairs, by=x->x[1], alg = RadixSort)
    end
    return [x[2] for x in pairs]
end

# function fsortperm3(svec::AbstractVector{String}, ::Type{T} = UInt) where T<:Unsigned
#     strlen = maximum(sizeof, svec)
#     strlen = max(strlen-sizeof(T), 0)

#     l = length(svec)

#     pairs = Vector{Pair{T, UInt32}}(l)
#     for (i, svec1) = enumerate(svec)
#         pairs[i] = Pair(load_uint(T, svec1, strlen), i)
#     end

#     sort!(pairs, by=x->x.first, alg = RadixSort)

#     while strlen > 0
#         strlen = max(strlen-4, 0)

#         for (i, svec1) = enumerate(svec)
#             pairs[i] = Pair(load_uint(T, svec1, strlen), pairs[2][2])
#         end
#         sort!(pairs, by=x->x.first, alg = RadixSort)
#     end
#     return [x.second for x in pairs]
# end


text:
module LoopVectorization

using Static: StaticInt, gt
using VectorizationBase, SLEEFPirates, UnPack, OffsetArrays
using VectorizationBase: register_size, register_count, cache_linesize, cache_size, has_opmask_registers,
    mask, pick_vector_width, MM, AbstractMask, data, grouped_strided_pointer,
    maybestaticlength, maybestaticsize, staticm1, staticp1, staticmul, vzero,
    maybestaticrange, offsetprecalc, lazymul,
    maybestaticfirst, maybestaticlast, scalar_less, scalar_greaterequal, gep, gesp, NativeTypes, #llvmptr,
    vfmadd, vfmsub, vfnmadd, vfnmsub, vfmadd_fast, vfmsub_fast, vfnmadd_fast, vfnmsub_fast, vfmadd231, vfmsub231, vfnmadd231, vfnmsub231,
    vfma_fast, vmuladd_fast, vdiv_fast, vadd_fast, vsub_fast, vmul_fast,
    relu, stridedpointer, StridedPointer, StridedBitPointer, AbstractStridedPointer, _vload, _vstore!,
    reduced_add, reduced_prod, reduce_to_add, reduce_to_prod, reduced_max, reduced_min, reduce_to_max, reduce_to_min,
    reduced_all, reduced_any, reduce_to_all, reduce_to_any,
    vsum, vprod, vmaximum, vminimum, vany, vall, unwrap, Unroll, VecUnroll,
    preserve_buffer, zero_vecunroll, vbroadcast_vecunroll, _vzero, _vbroadcast,
    contract_add, collapse_add,
    contract_mul, collapse_mul,
    contract_max, collapse_max,
    contract_min, collapse_min,
    contract_and, collapse_and,
    contract_or,  collapse_or,
    num_threads, num_cores,
    max_mask

using IfElse: ifelse

using ThreadingUtilities, CheapThreads
using SLEEFPirates: pow
using Base.Broadcast: Broadcasted, DefaultArrayStyle
using LinearAlgebra: Adjoint, Transpose, Diagonal
using Base.Meta: isexpr
using DocStringExtensions
import LinearAlgebra # for check_args

using Base.FastMath: add_fast, sub_fast, mul_fast, div_fast, inv_fast, abs2_fast, rem_fast, max_fast, min_fast, log_fast, log2_fast, log10_fast

using ArrayInterface
using ArrayInterface: OptionallyStaticUnitRange, OptionallyStaticRange, Zero, One, StaticBool, True, False, reduce_tup, indices, UpTri, LoTri
@static if VERSION ≥ v"1.6.0-rc1" #TODO: delete `else` when dropping 1.5 support
    using ArrayInterface: static_step
else # Julia 1.5 did not define `step` on CartesianIndices
    @inline static_step(x) = ArrayInterface.static_step
    static_step(x::CartesianIndices) = CartesianIndex(map(step, x.indices))
end
const Static = ArrayInterface.StaticInt

using Requires

export LowDimArray, stridedpointer, indices,
    @avx, @avxt, @_avx, *ˡ, _avx_!,
    vmap, vmap!, vmapt, vmapt!, vmapnt, vmapnt!, vmapntt, vmapntt!,
    tanh_fast, sigmoid_fast,
    vfilter, vfilter!, vmapreduce, vreduce

const VECTORWIDTHSYMBOL, ELTYPESYMBOL, MASKSYMBOL = Symbol("##Wvecwidth##"), Symbol("##Tloopeltype##"), Symbol("##mask##")

include("vectorizationbase_compat/contract_pass.jl")
include("vectorizationbase_compat/subsetview.jl")
include("closeopen.jl")
include("getconstindexes.jl")
include("predicates.jl")
include("simdfunctionals/map.jl")
include("simdfunctionals/filter.jl")
include("modeling/costs.jl")
include("modeling/operations.jl")
include("modeling/graphs.jl")
include("codegen/operation_evaluation_order.jl")
include("parse/memory_ops_common.jl")
include("parse/add_loads.jl")
include("parse/add_stores.jl")
include("parse/add_compute.jl")
include("parse/add_constants.jl")
include("parse/add_ifelse.jl")
include("modeling/determinestrategy.jl")
include("codegen/line_number_nodes.jl")
include("codegen/loopstartstopmanager.jl")
include("codegen/lower_compute.jl")
include("codegen/lower_constant.jl")
include("codegen/lower_memory_common.jl")
include("codegen/lower_load.jl")
include("codegen/lower_store.jl")
include("codegen/lowering.jl")
include("codegen/split_loops.jl")
include("codegen/lower_threads.jl")
include("condense_loopset.jl")
include("reconstruct_loopset.jl")
include("constructors.jl")
include("user_api_conveniences.jl")
include("simdfunctionals/mapreduce.jl")
include("broadcast.jl")

"""
LoopVectorization provides macros and functions that combine SIMD vectorization and
loop-reordering so as to improve performance:

- [`@avx`](@ref): transform `for`-loops and broadcasting
- [`@_avx`](@ref): similar to `@avx` but does not use type information
- [`vmap`](@ref) and `vmap!`: vectorized version of `map` and `map!`
- [`vmapnt`](@ref) and `vmapnt!`: non-temporal variants of `vmap` and `vmap!`
- [`vmapntt`](@ref) and `vmapntt!`: threaded variants of `vmapnt` and `vmapnt!`
- [`vfilter`](@ref) and `vfilter!`: vectorized versions of `filter` and `filter!`
"""
LoopVectorization

include("precompile.jl")
_precompile_()

# import ChainRulesCore, ForwardDiff
# include("vmap_grad.jl")
function __init__()
    @require ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" begin
        @require ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" include("simdfunctionals/vmap_grad.jl")
    end
    @require SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" begin
        @eval SpecialFunctions.erf(x::VectorizationBase.AbstractSIMD) = VectorizationBase.verf(float(x))
    end
end

end # module


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.0
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # $\chi^2$ distribution
#
#
# The Chi squared distribution $\chi^2(k)$ with some integer $k$ degrees of freedom is the sum of squares of $k$ independent standard Normal variates.
#

# +
using Distributions
using Plots
using SpecialFunctions

using LaTeXStrings
# -

?Chisq

# # Prepare standard normal distributions

d = Normal(0,1)

# ## Plot histogram

sample = rand(d, 10000)
p = plot(xlim=[-3, 3])
histogram!(sample, normalize=:pdf, label="hist")
plot!(x->pdf(d, x),label="pdf")

# # Prepare $\Gamma$ function

?gamma

Γ(x) = gamma(x)

function T(n, z)
    if z > 0
        a = z^(n/2 - 1)*exp(-z/2)
        b = 2^(n/2) * Γ(n/2)
        return a/b
    else
        return 0.
    end
end

# # Generate $\chi^2$ distribution with ...

# ## degree = 1

# +
n_trial = 100000
X = rand(d, n_trial)
X² = X .^ 2.
p = plot()
histogram!(p, X², normalize=:pdf, label=L"X^2")
x = range(0.01,2,length=100)
y = T.(1, x)
plot!(p, x, y, xlim=[0,2], ylim=[0,4], label=L"T_1")

y = pdf.(Chisq(1), x)
plot!(p, x, y, line=:dash, label=L"\chi^2(1)")

# -

# ## degree = 2

# +
n_trial = 100000
X = rand(d, n_trial)
Y = rand(d, n_trial)
X² = X .^ 2
Y² = Y .^ 2
p=plot()

histogram!(p, X².+Y², normalize=:pdf, label=L"X^2+Y^2+Z^2")
x = range(0.01,2,length=100)
y = T.(2, x)
plot!(p, x, y,xlim=[0,5], ylim=[0,0.5], label=L"T_3")

χ² = Chisq(2)

y = pdf.(χ², x)
χ² = Chisq(2)
plot!(p, x, y, label=L"\chi^2(2)",line=:dash)
# -

# ## degree = 3
#

# +
n_trial = 100000
X = rand(d, n_trial)
Y = rand(d, n_trial)
Z = rand(d, n_trial)
X² = X .^ 2
Y² = Y .^ 2
Z² = Z .^ 2
p = plot()
histogram!(p, X².+Y².+Z², normalize=:pdf, label=L"X^2+Y^2+Z^2")
χ² = Chisq(3)

x = range(0.01,2,length=100)
y = T.(3, x)
plot!(p, x, y,xlim=[0,5], ylim=[0,0.5], label=L"T_3")
y = pdf.(χ², x)

plot!(p, x, y, label=L"\chi^2(3)",line=:dash)


text:
@testset "community" begin

    #Create community
    @test_nowarn begin
        m = @agent(
            0,

            id::Identity,
            [v,l,i]::Local,
            g::Global,
            ga::GlobalArray,

            Equation = v̇ = 34*dt ,
            UpdateLocal = l += 1,
            UpdateGlobal = begin
                g += Normal(1,2) + Normal(1,2) + Uniform(1,2)
                ga[3,2] += 3
            end,
            UpdateInteraction = i += 1,
            UpdateLocalInteraction = i += 1
        )

        m = compile(m)

        com = Community(m)
        com = Community(m,N=30)
    end

    #Get variables
    @test_nowarn begin
        m = @agent(
            0,

            [id1,id2,id3]::Identity,
            [l1,l2,l3]::Local,
            [g1,g2,g3]::Global,
            [ga1,ga2,ga3]::GlobalArray,
        )

        m = compile(m)

        com = Community(m,N=10)

        com.declaredSymbols_
        com.t
        com.N
        com.id1
        com.id2
        com.id3
        com.l1
        com.l2
        com.l3
        com.g1
        com.g2
        com.g3
        com.ga1
        com.ga2
        com.ga3
    end

    #Get variables
    @test_nowarn begin
        m = @agent(
            0,

            [id1,id2,id3]::Identity,
            [l1,l2,l3]::Local,
            [g1,g2,g3]::Global,
            [ga1,ga2,ga3]::GlobalArray,
        )

        m = compile(m)

        com = Community(m,N=10)

        com.id1[1] = 1
        com.id2[1] = 2
        com.id3[1] = 3
        com.l1[1] = 1
        com.l2[1] = 2
        com.l3[1] = 3
        com.ga1 = ones(2,2)
        com.ga2 = 2*ones(2,2)
        com.ga3 = 3*ones(2,2)
        com.g1 = 1
        com.g2 = 2
        com.g3 = 3

        if com.id1[1] != 1 error() end
        if com.id2[1] != 2 error() end
        if com.id3[1] != 3 error() end
        if com.l1[1] != 1. error() end
        if com.l2[1] != 2. error() end
        if com.l3[1] != 3. error() end
        if com.ga1[1,1] != 1. error() end
        if com.ga2[1,1] != 2. error() end
        if com.ga3[1,1] != 3. error() end
        if com.g1 != 1. error() end
        if com.g2 != 2. error() end
        if com.g3 != 3. error() end

        com.id1 .= 1
        com.id2 .= 2
        com.id3 .= 3
        com.l1 .= 1
        com.l2 .= 2
        com.l3 .= 3
        com.ga1 .= 4
        com.ga2 .= 5
        com.ga3 .= 6

        if com.id1 != 1*ones(com.N) error() end
        if com.id2 != 2*ones(com.N) error() end
        if com.id3 != 3*ones(com.N) error() end
        if com.l1 != 1. *ones(com.N) error() end
        if com.l2 != 2. *ones(com.N) error() end
        if com.l3 != 3. *ones(com.N) error() end
        if com.ga1 != 4. *ones(2,2) error() end
        if com.ga2 != 5. *ones(2,2) error() end
        if com.ga3 != 6. *ones(2,2) error() end
    end

    @test_throws ErrorException begin
        m = @agent(0,id::Identity)
        m = compile(m)
        com = Community(m,N=10)

        com.id = 1
    end

    @test_throws MethodError begin
        m = @agent(0,id::Identity)
        m = compile(m)
        com = Community(m,N=10)

        com.id = [1]
    end

    @test_throws ErrorException begin
        m = @agent(0,id::Local)
        m = compile(m)
        com = Community(m,N=10)

        com.id = 1
    end

    @test_throws MethodError begin
        m = @agent(0,id::Local)
        m = compile(m)
        com = Community(m,N=10)

        com.id = [1]
    end

    @test_throws ErrorException begin
        m = @agent(0,ga::GlobalArray)
        m = compile(m)
        com = Community(m,N=10)

        com.ga = 1
    end

end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
function random_walk(tuning_run::Run)
    initial_cost = tuning_run.starting_cost
    x            = deepcopy(tuning_run.starting_point)
    name         = "Random Walk"
    cost_calls   = 0
    neighbor!(x)
    f_x = @fetch (tuning_run.measurement_method(tuning_run, x))
    cost_calls += tuning_run.cost_evaluations
    return Result(name, tuning_run.starting_point, x, f_x, 1,
                  1, cost_calls, false)
end


text: isautodifferentiable(alg::DEAlgorithm) = false


text:
#=
Created on Sunday 5 January 2019
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Build all the notebooks and stuff. Based on a file by Bram De Jaegher.
=#

t₀ = time()

using Weave

dirs = [#="00.Introduction",
        "01.Brackets",
        "02.Quadratic",
        "03.Unconstrained",
        "04.Constrained",
        "05.AutoDiff",
        "06.OptimalTransport",
        "07.MST",
        "08.ShortestPath",=#
        "09.NP-Complete",
        "11.TSP"
        ]

function cleanTemps(filename, dir; exts=[".out", ".log", ".aux"])
  for ext in exts
    rm(dir * filename * ext)
  end
end

for dir in dirs
  dir = "chapters/$dir/"
  files = readdir(dir)
  filter!(fn -> occursin(".jmd", fn), files)
  for filename in files
    filename = filename[1:end-4]
    println("Building $(filename)...")
    weave(dir * filename * ".jmd"; doctype="md2pdf")
    cleanTemps(filename, dir)
    convert_doc(dir * filename * ".jmd", dir * filename * ".ipynb")
  end
end

# FIGURES

include("scripts/maketestfigs.jl")
include("scripts/monge.jl")
include("scripts/unconstrained.jl")
include("scripts/colortransfer.jl")

tₑ = time()

t = Int(round(tₑ - t₀))

println("Building took $(div(t, 60))m$(t % 60)s")

max_stars_repo_name: dimiboeckaerts/STMO

max_issues_repo_name: dimiboeckaerts/STMO

max_forks_repo_name: dimiboeckaerts/STMO


text:
using ITensors
using ITensorGPU

import ITensors: space, state, op

space(::SiteType"Qubit") = 2
state(::SiteType"Qubit", ::StateName"0") = 1
state(::SiteType"Qubit", ::StateName"1") = 2

op_matrix(s::String) = op_matrix(OpName(s))

op_matrix(::OpName"Id") = [
  1 0
  0 1
]

op_matrix(::OpName"I") = op_matrix("Id")

op_matrix(::OpName"X") = [
  0 1
  1 0
]

function op_matrix(on::OpName, s::Index...; kwargs...)
  rs = reverse(s)
  return itensor(op_matrix(on; kwargs...), prime.(rs)..., dag.(rs)...)
end

op_matrix(gn::String, s::Index...; kwargs...) = op_matrix(OpName(gn), s...; kwargs...)

function op_matrix(gn::String, s::Vector{<:Index}, ns::Int...; kwargs...)
  return op_matrix(OpName(gn), s[[ns...]]...; kwargs...)
end

op(gn::OpName, ::SiteType"Qubit", s::Index...; kwargs...) = op_matrix(gn, s...; kwargs...)

N = 10
s = siteinds("Qubit", N)
X = cu.(ops(s, [("X", n) for n in 1:N]))

initstate = fill("0", N)

ψ0 = productCuMPS(s, initstate)

gates = [X[n] for n in 1:2:N]
ψ = apply(gates, ψ0)


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# Maximal Coverage Location Problem applied to creating a nested compliance table

@testset "nested mclp" begin
	## tiny
	pointDemands = [1.0]
	coverMatrix = fill(true, 1, 1)
	weights = [1.0]
	results = Dict()
	@test solveNestedMclp(pointDemands, coverMatrix, weights; results = results) == [1]

	## small
	pointDemands = Float[5,3,2,1]
	coverMatrix = convert(Array{Bool,2}, [ # coverMatrix[i,j] = true if facility location i can cover point j
		1 1 0 0
		0 1 1 0
		1 0 0 1
	])

	weights = [1,1,1]/3
	@test solveNestedMclp(pointDemands, coverMatrix, weights, results = results) == [1,2,3]
	results[:objVal] = sum(weights .* [5+3, 5+3+2, 5+3+2+1])

	weights = [1,4,1]/6
	@test solveNestedMclp(pointDemands, coverMatrix, weights, results = results) == [3,2,1]
	results[:objVal] = sum(weights .* [5+1, 5+3+2+1, 5+3+2+1])
end


text:
using ..Helper: init

function simulate!(x, dt, sigma, rng)
    sgm = sigma * sqrt(dt)
    steps, = size(x)
    for t = 2:steps
        x[t] = (1 - dt) * x[t-1] + sgm * randn(rng)
    end
end

function prepare(name)
    res = init(name)
    rng = MersenneTwister(res[:seed])
    function run()
        simulate!(res[:x], res[:dt], res[:sigma], rng)
        res
    end
end

max_stars_repo_name: tkf/comparatist

max_issues_repo_name: tkf/comparatist

max_forks_repo_name: tkf/comparatist


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## initdefs.jl - initialization and runtime management definitions

"""
    PROGRAM_FILE

A string containing the script name passed to Julia from the command line. Note that the
script name remains unchanged from within included files. Alternatively see
[`@__FILE__`](@ref).
"""
global PROGRAM_FILE = ""

"""
    ARGS

An array of the command line arguments passed to Julia, as strings.
"""
const ARGS = String[]

"""
    exit(code=0)

Stop the program with an exit code. The default exit code is zero, indicating that the
program completed successfully. In an interactive session, `exit()` can be called with
the keyboard shortcut `^D`.
"""
exit(n) = ccall(:jl_exit, Cvoid, (Int32,), n)
exit() = exit(0)

const roottask = current_task()

is_interactive = false

"""
    isinteractive() -> Bool

Determine whether Julia is running an interactive session.
"""
isinteractive() = (is_interactive::Bool)

## package depots (registries, packages, environments) ##

const DEPOT_PATH = String[]

function init_depot_path(BINDIR::String = Sys.BINDIR)
    if haskey(ENV, "JULIA_DEPOT_PATH")
        depots = split(ENV["JULIA_DEPOT_PATH"], Sys.iswindows() ? ';' : ':')
        append!(empty!(DEPOT_PATH), map(expanduser, depots))
    else
        push!(empty!(DEPOT_PATH), joinpath(homedir(), ".julia"))
        push!(DEPOT_PATH, abspath(BINDIR, "..", "local", "share", "julia"))
        push!(DEPOT_PATH, abspath(BINDIR, "..", "share", "julia"))
    end
end

## load-path types ##

abstract type AbstractEnv end

struct CurrentEnv <: AbstractEnv
    create::Bool
    CurrentEnv(; create::Bool=false) = new(create)
end

struct NamedEnv <: AbstractEnv
    name::String
    create::Bool
    NamedEnv(name::String; create::Bool=false) = new(name, create)
end

function show(io::IO, env::CurrentEnv)
    print(io, CurrentEnv, "(")
    env.create && print(io, "create=true")
    print(io, ")")
end

function show(io::IO, env::NamedEnv)
    print(io, NamedEnv, "(", repr(env.name))
    env.create && print(io, ", create=true")
    print(io, ")")
end

function parse_env(env::Union{String,SubString{String}})
    isempty(env) && return Any[]
    env == "@" && return CurrentEnv()
    env == "@!" && return CurrentEnv(create=true)
    if env[1] == '@'
        create = env[2] == '!'
        name = env[2+create:end]
        name = replace(name, '#' => VERSION.major, count=1)
        name = replace(name, '#' => VERSION.minor, count=1)
        name = replace(name, '#' => VERSION.patch, count=1)
        return NamedEnv(name, create=create)
    end
    return env # literal path
end

"""
    LOAD_PATH

An array of paths as strings or custom loader objects for the `require`
function and `using` and `import` statements to consider when loading
code.
"""
const LOAD_PATH = Any[]

function parse_load_path(str::String)
    envs = Any[split(str, Sys.iswindows() ? ';' : ':');]
    for (i, env) in enumerate(envs)
        if '|' in env
            envs[i] = Any[parse_env(e) for e in split(env, '|')]
        else
            envs[i] = parse_env(env)
        end
    end
    return envs
end

const default_named = parse_load_path("@v#.#.#|@v#.#|@v#|@default|@!v#.#")

function init_load_path(BINDIR::String = Sys.BINDIR)
    vers = "v$(VERSION.major).$(VERSION.minor)"
    stdlib = abspath(BINDIR, "..", "share", "julia", "stdlib", vers)
    load_path = Base.creating_sysimg           ? Any[stdlib] :
                haskey(ENV, "JULIA_LOAD_PATH") ? parse_load_path(ENV["JULIA_LOAD_PATH"]) :
                                                 Any[CurrentEnv(); default_named; stdlib]
    append!(empty!(LOAD_PATH), load_path)
end

const atexit_hooks = []

"""
    atexit(f)

Register a zero-argument function `f()` to be called at process exit. `atexit()` hooks are
called in last in first out (LIFO) order and run before object finalizers.
"""
atexit(f::Function) = (pushfirst!(atexit_hooks, f); nothing)

function _atexit()
    for f in atexit_hooks
        try
            f()
        catch err
            show(stderr, err)
            println(stderr)
        end
    end
end

max_stars_repo_name: ninjin/julia

max_issues_repo_name: ninjin/julia

max_forks_repo_name: ninjin/julia


text:
#P6
#1920 1080
#255
#<frame RGB>
header = zeros(UInt8, 128)
read!(stdin, header)
println(header)
# ffmpeg.exe -i foobar.mp4 -f image2pipe -vcodec ppm pipe:1 | julia ppmhead.jl


text:
module problem126

export Cuboid

description = """
Cuboid layers

The minimum number of cubes to cover every visible face on a cuboid measuring 3 x 2 x 1 is twenty-two.
If we then add a second layer to this solid it would require forty-six cubes to cover every visible face,
the third layer would require seventy-eight cubes, and the fourth layer would require one-hundred and eighteen cubes to cover every visible face.

However, the first layer on a cuboid measuring 5 x 1 x 1 also requires twenty-two cubes;
similarly the first layer on cuboids measuring 5 x 3 x 1, 7 x 2 x 1, and 11 x 1 x 1 all contain forty-six cubes.

We shall define C(n) to represent the number of cuboids that contain n cubes in one of its layers.
So C(22) = 2, C(46) = 4, C(78) = 5, and C(118) = 8.
It turns out that 154 is the least value of n for which C(n) = 10.

Find the least value of n for which C(n) = 1000.
"""

# A cube with length, width and height
immutable Cuboid
  l :: Int
  w :: Int
  h :: Int
end

# calculates the middle section of the n'th layer
# equal to (the perimeter of original cube + diagonals) * height or original cube
middle(c :: Cuboid, n :: Int) = (2 * (c.l + c.w) + 4 * (n - 1)) * c.h

# Recurrance relation for calculating the top section of the n'th layer
# For n = 1, just the area of the top of original cube
# otherwise, equal to the top of the previous layer + a one cube thick rim from the middle of previous layer
top(c :: Cuboid, n :: Int) =
  if (n == 1)
    return (c.l * c.w)
  else
    return (top(c, n-1) + div(middle(c, n-1), c.h))
  end

# Calculate the n'th layer, top and bottom are symetric + middle section
layer(c :: Cuboid, n :: Int) = 2 * top(c, n) + middle(c, n)

# Calculate the size of all layers of cuboid up to a maximum
# re-uses the results of previous iterations to efficiently calculate the next
function layers(c :: Cuboid, maxsize :: Int)
  n = 1
  t = top(c, n)
  m = middle(c, n)
  x = 0
  result = Int[]
  while true
    x = 2 * t + m
    if x > maxsize break end
    push!(result, x)
    t = t + div(m, c.h)
    n += 1
    m = middle(c, n)
  end
  return result
end

# Search for a solution x, to C(x) == N, for x in 1:guess
# The lower bound is given by layer 1 of cuboid (Length, 1, 1)
# Once a solution is found < lower bound, search can stop
function search(N :: Int, guess :: Int)
  counts = zeros(Int, guess)
  l = 1
  lowerbound = layer(Cuboid(l, 1, 1), 1)
  while lowerbound <= guess
    for i = 1:lowerbound-1
      if counts[i] == N
        return i
      end
    end
    w = 1
    while w <= l && layer(Cuboid(l, w, 1), 1) <= guess
      h = 1
      while h <= w && layer(Cuboid(l, w, h), 1) <= guess
        sizes = layers(Cuboid(l,w,h), guess)
        counts[sizes] += 1
        h += 1
      end
      w += 1
    end
    l += 1
    lowerbound = layer(Cuboid(l, 1, 1), 1)
  end
  println("no solution found < $guess")
  return -1
end

using Base.Test

@test layer(Cuboid(3, 2, 1), 1) == 22
@test layer(Cuboid(3, 2, 1), 2) == 46
@test layer(Cuboid(3, 2, 1), 3) == 78
@test layer(Cuboid(3, 2, 1), 4) == 118
@test layer(Cuboid(5, 1, 1), 1) == 22
@test layer(Cuboid(5, 3, 1), 1) == 46
@test layer(Cuboid(7, 2, 1), 1) == 46
@test layer(Cuboid(11, 1, 1), 1) == 46
@test layers(Cuboid(3,2,1), 78) == [22, 46, 78]
@test search(10, 200) == 154

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
using Documenter, Literate
import LocalStore

ENV["JULIA_DEBUG"] = "Documenter,Literate,LocalStore"

const literate_dir = joinpath(@__DIR__, "src/literate")

function clear_md_files(dir::String)
    for file in readdir(dir; join=true)
        if endswith(file, ".md")
            rm(file)
        end
    end
end

clear_md_files(literate_dir)

for file in readdir(literate_dir; join=true)
    if endswith(file, ".jl")
        Literate.markdown(file, literate_dir)
    end
end

makedocs(
    modules = [LocalStore],
    sitename = "LocalStore.jl",
    pages = [
        "Home" => "index.md",
        "Examples" => [
            "Save / load" => "literate/io.md",
        ],
        "Reference" => "reference.md"
    ],
    strict = true
)

clear_md_files(literate_dir)

deploydocs(
    repo = "github.com/cossio/LocalStore.jl.git",
    devbranch = "master"
)


text:
module TailRiskScenGen

# using Distributions, SkewDist
using Distributions
using Cones
using Gurobi
using Clustering
using DataStructures
using StatsFuns

import Base: ∈, length
import PDMats: dim

export aggregate_scenarios, Cone, FiniteCone, PolyhedralCone, project, AbstractRiskRegion, EllipticalRiskRegion, MonotonicEllipticalRiskRegion, transformed_size, in_RiskRegion, aggregation_sampling, nonrisk_clustering, dim, cone_from_constraints, quota_cone, MonotonicRiskRegion, ExactMonotonicRiskRegion, SurvivorApproximator

const VecF64 = Vector{Float64}
const MatF64 = Matrix{Float64}

include("utilities.jl")
include("RiskRegion.jl")
include("monotonic_utils.jl")
include("monotonic_risk_regions.jl")
include("nonrisk_clustering.jl")
include("BruteForceRiskRegion.jl")
include("ScenarioReduction.jl")

end # module


text:
module StencilFunctorTest

using StencilFunctor
using Base.Test

function main()
    # Array
    let
        diff(c, bm, bp) = bp - bm
        bc(c, dir) = c
        @test stencil_fmap(diff, bc, Int[], 0, 1) == Int[]
        @test stencil_fmap(diff, bc, [1], 0, 2) == [2]
        @test stencil_fmap(diff, bc, [1,2], 0, 3) == [2,2]
        @test stencil_fmap(diff, bc, [1,2,3,4,5], 0, 6) == [2,2,2,2,2]
    end
end

end


text:
module TestFormula
    using Base.Test
    using DataFrames

    # TODO:
    # - grouped variables in formulas with interactions
    # - is it fast?  Can expand() handle DataFrames?
    # - deal with intercepts
    # - implement ^2 for datavector's
    # - support more transformations with I()?

    ## Formula parsing
    import DataFrames.Terms

    ## totally empty
    t = Terms(Formula(nothing, 0))
    @test t.response == false
    @test t.intercept == false
    @test t.terms == []
    @test t.eterms == []

    ## empty RHS
    t = Terms(y ~ 0)
    @test t.intercept == false
    @test t.terms == []
    @test t.eterms == [:y]
    t = Terms(y ~ -1)
    @test t.intercept == false
    @test t.terms == []

    ## intercept-only
    t = Terms(y ~ 1)
    @test t.response == true
    @test t.intercept == true
    @test t.terms == []
    @test t.eterms == [:y]

    ## terms add
    t = Terms(y ~ 1 + x1 + x2)
    @test t.intercept == true
    @test t.terms == [:x1, :x2]
    @test t.eterms == [:y, :x1, :x2]

    ## implicit intercept behavior:
    t = Terms(y ~ x1 + x2)
    @test t.intercept == true
    @test t.terms == [:x1, :x2]
    @test t.eterms == [:y, :x1, :x2]

    ## no intercept
    t = Terms(y ~ 0 + x1 + x2)
    @test t.intercept == false
    @test t.terms == [:x1, :x2]

    t = Terms(y ~ -1 + x1 + x2)
    @test t.intercept == false
    @test t.terms == [:x1, :x2]

    t = Terms(y ~ x1 & x2)
    @test t.terms == [:(x1 & x2)]
    @test t.eterms == [:y, :x1, :x2]

    ## `*` expansion
    t = Terms(y ~ x1 * x2)
    @test t.terms == [:x1, :x2, :(x1 & x2)]
    @test t.eterms == [:y, :x1, :x2]

    ## associative rule:
    ## +
    t = Terms(y ~ x1 + x2 + x3)
    @test t.terms == [:x1, :x2, :x3]

    ## &
    t = Terms(y ~ x1 & x2 & x3)
    @test t.terms == [:((&)(x1, x2, x3))]
    @test t.eterms == [:y, :x1, :x2, :x3]

    ## distributive property of + and &
    t = Terms(y ~ x1 & (x2 + x3))
    @test t.terms == [:(x1&x2), :(x1&x3)]

    ## FAILS: ordering of expanded interaction terms is wrong
    ## (only has an observable effect when both terms are categorical and
    ## produce multiple model matrix columns that are multiplied together...)
    ##
    ## t = Terms(y ~ (x2 + x3) & x1)
    ## @test t.terms == [:(x2&x1), :(x3&x1)]

    ## three-way *
    t = Terms(y ~ x1 * x2 * x3)
    @test t.terms == [:x1, :x2, :x3,
                      :(x1&x2), :(x1&x3), :(x2&x3),
                      :((&)(x1, x2, x3))]
    @test t.eterms == [:y, :x1, :x2, :x3]

    ## Interactions with `1` reduce to main effect.  All fail at the moment.
    ## t = Terms(y ~ 1 & x1)
    ## @test t.terms == [:x1]              # == [:(1 & x1)]
    ## @test t.eterms == [:y, :x1]

    ## t = Terms(y ~ (1 + x1) & x2)
    ## @test t.terms == [:x2, :(x1&x2)]    # == [:(1 & x1)]
    ## @test t.eterms == [:y, :x1, :x2]

    ## Tests for constructing ModelFrame and ModelMatrix

    d = DataFrame()
    d[:y] = [1:4;]
    d[:x1] = [5:8;]
    d[:x2] = [9:12;]
    d[:x3] = [13:16;]
    d[:x4] = [17:20;]

    x1 = [5.:8;]
    x2 = [9.:12;]
    x3 = [13.:16;]
    x4 = [17.:20;]
    f = y ~ x1 + x2
    mf = ModelFrame(f, d)
    ## @test mm.response_colnames == ["y"] # nope: no response_colnames
    @test coefnames(mf) == ["(Intercept)","x1","x2"]
    ## @test model_response(mf) == transpose([1. 2 3 4]) # fails: Int64 vs. Float64
    mm = ModelMatrix(mf)
    @test mm.m[:,1] == ones(4)
    @test mm.m[:,2:3] == [x1 x2]

    #test_group("expanding a PooledVec into a design matrix of indicators for each dummy variable")

    d[:x1p] = PooledDataArray(d[:x1])
    mf = ModelFrame(y ~ x1p, d)
    mm = ModelMatrix(mf)

    @test mm.m[:,2] == [0, 1., 0, 0]
    @test mm.m[:,3] == [0, 0, 1., 0]
    @test mm.m[:,4] == [0, 0, 0, 1.]
    @test coefnames(mf)[2:end] == ["x1p - 6", "x1p - 7", "x1p - 8"]

    #test_group("create a design matrix from interactions from two DataFrames")
    ## this was removed in commit dead4562506badd7e84a2367086f5753fa49bb6a

    ## b = DataFrame()
    ## b["x2"] = DataVector(x2)
    ## df = interaction_design_matrix(a,b)
    ## @test df[:,1] == DataVector([0, 10., 0, 0])
    ## @test df[:,2] == DataVector([0, 0, 11., 0])
    ## @test df[:,3] == DataVector([0, 0, 0, 12.])

    #test_group("expanding an singleton expression/symbol into a DataFrame")
    ## generalized expand was dropped, too
    ## df = deepcopy(d)
    ## r = expand(:x2, df)
    ## @test isa(r, DataFrame)
    ## @test r[:,1] == DataVector([9,10,11,12])  # TODO: test float vs int return

    ## df = deepcopy(d)
    ## ex = :(log(x2))
    ## r = expand(ex, df)
    ## @test isa(r, DataFrame)
    ## @test r[:,1] == DataVector(log([9,10,11,12]))

    # ex = :(x1 & x2)
    # r = expand(ex, df)
    # @test isa(r, DataFrame)
    # @test ncol(r) == 1
    # @test r[:,1] == DataArray([45, 60, 77, 96])

    ## r = expand(:(x1 + x2), df)
    ## @test isa(r, DataFrame)
    ## @test ncol(r) == 2
    ## @test r[:,1] == DataVector(df["x1"])
    ## @test r[:,2] == DataVector(df["x2"])

    ## df["x1"] = PooledDataArray(x1)
    ## r = expand(:x1, df)
    ## @test isa(r, DataFrame)
    ## @test ncol(r) == 3
    ## @test r == expand(PooledDataArray(x1), "x1", DataFrame())

    ## r = expand(:(x1 + x2), df)
    ## @test isa(r, DataFrame)
    ## @test ncol(r) == 4
    ## @test r[:,1:3] == expand(PooledDataArray(x1), "x1", DataFrame())
    ## @test r[:,4] == DataVector(df["x2"])

    ## df["x2"] = PooledDataArray(x2)
    ## r = expand(:(x1 + x2), df)
    ## @test isa(r, DataFrame)
    ## @test ncol(r) == 6
    ## @test r[:,1:3] == expand(PooledDataArray(x1), "x1", DataFrame())
    ## @test r[:,4:6] == expand(PooledDataArray(x2), "x2", DataFrame())

    #test_group("Creating a model matrix using full formulas: y ~ x1 + x2, etc")

    df = deepcopy(d)
    f = y ~ x1 & x2
    mf = ModelFrame(f, df)
    mm = ModelMatrix(mf)
    @test mm.m == [ones(4) x1.*x2]

    f = y ~ x1 * x2
    mf = ModelFrame(f, df)
    mm = ModelMatrix(mf)
    @test mm.m == [ones(4) x1 x2 x1.*x2]

    df[:x1] = PooledDataArray(x1)
    x1e = [[0, 1, 0, 0] [0, 0, 1, 0] [0, 0, 0, 1]]
    f = y ~ x1 * x2
    mf = ModelFrame(f, df)
    mm = ModelMatrix(mf)
    @test mm.m == [ones(4) x1e x2 [0, 10, 0, 0] [0, 0, 11, 0] [0, 0, 0, 12]]

    #test_group("Basic transformations")

    ## the log(x2) appears to be broken (again, in The Great Purge, by the
    ## removal of the x -> with(d, x) in the ModelFrame constructor)
    ## df = deepcopy(d)
    ## f = y ~ x1 + log(x2)
    ## mf = ModelFrame(f, df)
    ## mm = ModelMatrix(mf)
    ## @test mm.m == [ones(4) x1 log(x2)]

    ## df = deepcopy(d)
    ## df["x1"] = PooledDataArray([5:8])
    ## f = Formula(:(y ~ x1 * (log(x2) + x3)))
    ## mf = ModelFrame(f, df)
    ## mm = ModelMatrix(mf)
    ## @test mm.model_colnames == [
    ##  "(Intercept)"
    ##  "x1:6"
    ##  "x1:7"
    ##  "x1:8"
    ##  "log(x2)"
    ##  "x3"
    ##  "x1:6&log(x2)"
    ##  "x1:6&x3"
    ##  "x1:7&log(x2)"
    ##  "x1:7&x3"
    ##  "x1:8&log(x2)"
    ##  "x1:8&x3" ]

    #test_group("Model frame response variables")

    ## also does not work, not sure what used to happen but now it seems that
    ## the LHS is assumed to be just a single symbol that indexes a DF column

    ## f = x1 + x2 ~ y + x3
    ## mf = ModelFrame(f, d)
    ## @test mf.y_indexes == [1, 2]
    ## @test isequal(mf.formula.lhs, [:(x1 + x2)])
    ## @test isequal(mf.formula.rhs, [:(y + x3)])

    # unique_symbol tests
    #@test unique_symbols(:(x1 + x2)) ==  {"x1"=>:x1, "x2"=>:x2}
    #unique_symbols(:(y ~ x1 + x2 + x3))

    # additional tests from Tom
    y = [1., 2, 3, 4]
    mf = ModelFrame(y ~ x2, d)
    mm = ModelMatrix(mf)
    @test mm.m == [ones(4) x2]
    ## @test model_response(mf) == y''     # fails: Int64 vs. Float64

    df = deepcopy(d)
    df[:x1] = PooledDataArray(df[:x1])

    f = y ~ x2 + x3 + x3*x2
    mm = ModelMatrix(ModelFrame(f, df))
    @test mm.m == [ones(4) x2 x3 x2.*x3]
    mm = ModelMatrix(ModelFrame(y ~ x3*x2 + x2 + x3, df))
    @test mm.m == [ones(4) x3 x2 x2.*x3]
    mm = ModelMatrix(ModelFrame(y ~ x1 + x2 + x3 + x4, df))
    @test mm.m[:,2] == [0, 1., 0, 0]
    @test mm.m[:,3] == [0, 0, 1., 0]
    @test mm.m[:,4] == [0, 0, 0, 1.]
    @test mm.m[:,5] == x2
    @test mm.m[:,6] == x3
    @test mm.m[:,7] == x4

    mm = ModelMatrix(ModelFrame(y ~ x2 + x3 + x4, df))
    @test mm.m == [ones(4) x2 x3 x4]
    mm = ModelMatrix(ModelFrame(y ~ x2 + x2, df))
    @test mm.m == [ones(4) x2]
    mm = ModelMatrix(ModelFrame(y ~ x2*x3 + x2&x3, df))
    @test mm.m == [ones(4) x2 x3 x2.*x3]
    mm = ModelMatrix(ModelFrame(y ~ x2*x3*x4, df))
    @test mm.m == [ones(4) x2 x3 x4 x2.*x3 x2.*x4 x3.*x4 x2.*x3.*x4]
    mm = ModelMatrix(ModelFrame(y ~ x2&x3 + x2*x3, df))
    @test mm.m == [ones(4) x2 x3 x2.*x3]

    f = y ~ x2 & x3 & x4
    mf = ModelFrame(f, df)
    mm = ModelMatrix(mf)
    @test mm.m == [ones(4) x2.*x3.*x4]

    f = y ~ x1 & x2 & x3
    mf = ModelFrame(f, df)
    mm = ModelMatrix(mf)
    @test mm.m[:, 2:end] == vcat(zeros(1, 3), diagm(x2[2:end].*x3[2:end]))

    #test_group("Column groups in formulas")
    ## set_group was removed in The Great Purge (55e47cd)

    ## set_group(d, "odd_predictors", ["x1","x3"])
    ## @test expand(:odd_predictors, d) == d["odd_predictors"]
    ## mf = ModelFrame(Formula(:(y ~ odd_predictors)), d)
    ## @test mf.df[:,1] == d["y"]
    ## @test mf.df[:,2] == d["x1"]
    ## @test mf.df[:,3] == d["x3"]
    ## @test ncol(mf.df) == 3
    ## mf = ModelFrame(Formula(:(y ~ odd_predictors * x2)), d)
    ## mm = ModelMatrix(mf)
    ## @test mm.model == [ones(4) x1 x3 x2 x1.*x2 x3.*x2]

    ## Interactions between three PDA columns
    df = DataFrame(y=1:27,
                   x1 = PooledDataArray(vec([x for x in 1:3, y in 4:6, z in 7:9])),
                   x2 = PooledDataArray(vec([y for x in 1:3, y in 4:6, z in 7:9])),
                   x3 = PooledDataArray(vec([z for x in 1:3, y in 4:6, z in 7:9])))
    f = y ~ x1 & x2 & x3
    mf = ModelFrame(f, df)
    @test coefnames(mf)[2:end] ==
        vec([string("x1 - ", x, " & x2 - ", y, " & x3 - ", z) for
             x in 2:3,
             y in 5:6,
             z in 8:9])

    mm = ModelMatrix(mf)
    @test mm.m[:,2] == 0. + (df[:x1] .== 2) .* (df[:x2] .== 5) .* (df[:x3].==8)
    @test mm.m[:,3] == 0. + (df[:x1] .== 3) .* (df[:x2] .== 5) .* (df[:x3].==8)
    @test mm.m[:,4] == 0. + (df[:x1] .== 2) .* (df[:x2] .== 6) .* (df[:x3].==8)
    @test mm.m[:,5] == 0. + (df[:x1] .== 3) .* (df[:x2] .== 6) .* (df[:x3].==8)
    @test mm.m[:,6] == 0. + (df[:x1] .== 2) .* (df[:x2] .== 5) .* (df[:x3].==9)
    @test mm.m[:,7] == 0. + (df[:x1] .== 3) .* (df[:x2] .== 5) .* (df[:x3].==9)
    @test mm.m[:,8] == 0. + (df[:x1] .== 2) .* (df[:x2] .== 6) .* (df[:x3].==9)
    @test mm.m[:,9] == 0. + (df[:x1] .== 3) .* (df[:x2] .== 6) .* (df[:x3].==9)

    ## Distributive property of :& over :+
    df = deepcopy(d)
    f = y ~ (x1+x2) & (x3+x4)
    mf = ModelFrame(f, df)
    mm = ModelMatrix(mf)
    @test mm.m == hcat(ones(4), x1.*x3, x1.*x4, x2.*x3, x2.*x4)

    ## Condensing nested :+ calls
    f = y ~ x1 + (x2 + (x3 + x4))
    @test ModelMatrix(ModelFrame(f, df)).m == hcat(ones(4), x1, x2, x3, x4)

    ## Extra levels in categorical column
    mf_full = ModelFrame(y ~ x1p, d)
    mm_full = ModelMatrix(mf_full)
    @test size(mm_full) == (4,4)

    mf_sub = ModelFrame(y ~ x1p, d[2:4, :])
    mm_sub = ModelMatrix(mf_sub)
    ## should have only three rows, and only three columns (intercept plus two
    ## levels of factor)
    @test size(mm_sub) == (3,3)

    ## Missing data
    d[:x1m] = @data [5, 6, NA, 7]
    mf = ModelFrame(y ~ x1m, d)
    mm = ModelMatrix(mf)
    @test mm.m[:, 2] == d[complete_cases(d), :x1m]

end


text:
using Random
using Statistics

export Spot, SpotDynamics, modulate

"""
Siderial rotation period at the equator
"""
const PROT_SUN = 24.5
const OMEGA_SUN = 2 * π / PROT_SUN / 86400

struct Spot{T <: Number}
    nday::Int # day index
    lat::T
    lon::T
    Φmax::T
    area_max::T
end

struct SpotDynamics{T <: Number}
    spots::Vector{<:Spot}
    duration::T
    inclination::T
    ω::T
    Δω::T
    equatorial_period::T
    τ_emergence::T
    τ_decay::T
end

"""
    SpotDynamics(spots::AbstractVector{Spot};
        duration = maximum([s.nday for s in spots]),
        alpha_med = 0.0001,
        inclination = asin(rand()),
        ω = 1.0,
        Δω = 0.2,
        τ_decay = 5.0,
        threshold = 0.1)

A container for the dynamics of starspots.

# Parameters
* spots - The list of `Spot`s evolved over time
* duration - The length of the evolution time
* alpha_med - An activation parameter for the magnetic flux
* inclination - Inclination of the star from our line of sight in radians
* ω - Rotational velocity of the star in solar units
* Δω - change in rotational velocity over the time in solar untis
* τ_decay - The decay timescale
* threshold - The threshold in magnetic flux for filtering starspots
"""
function SpotDynamics(spots::AbstractVector{Spot};
    duration = maximum([s.nday for s in spots]),
    alpha_med = 0.0001,
    inclination = asin(rand()),
    ω = 1.0,
    Δω = 0.2,
    τ_decay = 5.0,
    threshold = 0.1)

    ω *= OMEGA_SUN
    Δω *= OMEGA_SUN
    eq_per = 2π / ω / 86400
    τ_emergence = max(2.0, eq_per * τ_decay / 5)
    τ_decay *= eq_per

    filt(s) = (s.nday < duration) && (s.Φmax > threshold)
    spots = filter(filt, spots)
    Φmax = [s.Φmax for s in spots]
    area_max = alpha_med .* Φmax ./ median(Φmax)
    spots = [Spot(s.nday, s.lat, s.lon, s.Φmax, a) for (s, a) in zip(spots, area_max)]

    return SpotDynamics(spots, promote(duration, inclination, ω, Δω, eq_per, τ_emergence, τ_decay)...)
end

Base.broadcastable(r::SpotDynamics) = Ref(r)
Base.length(r::SpotDynamics) = length(r.spots)
Base.size(r::SpotDynamics) = size(r.spots)
Base.size(r::SpotDynamics, i) = size(r.spots, i)

function Base.show(io::IO, sd::SpotDynamics)
    out = """
    $(typeof(sd))
      nspots: $(length(sd))
      duration: $(sd.duration)
      inclination: $(sd.inclination)
      ω: $(sd.ω)
      Δω: $(sd.Δω)
      equatorial_period: $(sd.equatorial_period)
      τ_emergence: $(sd.τ_emergence)
      τ_decay: $(sd.τ_decay)"""
    print(io, out)
end
"""
    diffrot(ω₀, Δω, lat)

Default differental rotation function

Returns angular velocity as a function of latitude [0°, 90°]
"""
@inline diffrot(ω₀, Δω, lat) = ω₀ - Δω * sin(lat)^2

function modulate(spot::Spot, t,
    τ_emergence, τ_decay,
    ω, Δω, inclination, diffrot_func = diffrot)
    # Get spot area
    tt = t - spot.nday
    area = spot.area_max
    timescale = tt < 0 ? τ_emergence : τ_decay

    area *= exp(-(tt / timescale)^2 / 2)

    # rotation rate
    phase = diffrot(ω, Δω, spot.lat) * t * 86400 + spot.lon

    # foreshortening
    cos_beta = cos(inclination) * sin(spot.lat) + sin(inclination) * cos(spot.lat) * cos(phase)

    return -area * max(cos_beta, 0.0)
end

"""
    modulate(::SpotDynamics, time)

Modulate the flux due to starspots at the given timestep in days.
"""
function modulate(r::SpotDynamics, t)
    dFlux = sum(modulate.(r.spots, t, r.τ_emergence, r.τ_decay, r.ω, r.Δω, r.inclination))
    return dFlux
end

max_stars_repo_name: zclaytor/butterpy

max_issues_repo_name: zclaytor/butterpy

max_forks_repo_name: zclaytor/butterpy


text:
module Got
export demo
using Parameters

@with_kw mutable struct Cols
  x  = (all  = [], nums = [], syms = [])
  y = (all  = [], nums = [], syms = [],
           goals= [22], klass= nothing)
  all  = []
  nums = []
  syms = (10, 20,   30, 40)
  name=""
  pos=0
end

z=Cols(name="jane", pos=2)
#show(z)
#w = Cols(),
#dump(w)
#dump(nothing)

function num(s)
  v = tryparse(Float64,s)
  v == nothing ? s : v
end
function dd()
open("sherlock-holmes.txt") do f
   line = 1
   while !eof(f)
     x = readline(f)
     println("$line $x")
     line += 1
   end
 end
end

#[take(PSerie(2),10)...]
#
#struct PSerie
#           p
#       end
#
#
#function Base.iterate(ps::PSerie, (s,n) = (0,1))
#    s += 1 / (n^ps.p)
#           s, (s, n+1)
#       end
#

demo = begin
  open("got.jl") do file
      n=1
      old=""
      function prep(s)
        if sizeof(s) > 0
          println( map(num, split(s,",")))end end
      for line in eachline(file)
          line = replace(line, r"([ \t\n]|#.*)" => "")
          if sizeof(line) == 0  continue end
          print(line[end])
          if line[end] == ','
            old = old * line
          else
            prep(old*line)
            old=""
          end
      end
      prep(old)
  end
end
end

max_stars_repo_name: timm/misc

max_issues_repo_name: timm/misc

max_forks_repo_name: timm/misc


text:
module {{ class }}

  {{fields}}

end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.7.0-beta3
#     language: julia
#     name: julia-1.7
# ---

# %%
using Random, Statistics, Plots

function diff!(X)
    @inbounds for i in reverse(keys(X)[begin:end-1])
        X[i+1] -= X[i]
    end
    X
end

function randsimplex!(rng::AbstractRNG, X)
    X[end] = one(eltype(X))
    Y = @view X[begin:end-1]
    rand!(rng, Y)
    sort!(Y)
    X .*= length(X)
    diff!(X)
end
randsimplex(rng::AbstractRNG, N, T=Float64) = randsimplex!(rng, Vector{T}(undef, N))
randsimplex!(X) = randsimplex!(Random.default_rng(), X)
randsimplex(N, T=Float64) = randsimplex!(Vector{T}(undef, N))

n = 10^4
A = randsimplex(n)
@show mean(A), var(A)
histogram(A; norm=true, alpha=0.3, label="randsimplex($n)")
plot!(x->exp(-x), 0, 6; xlim=(-0.1, 6), lw=2, label="Exponential(1) dist")

# %%
struct Simplex{T} N::Int end
Simplex(N) = Simplex{Float64}(N)

function Random.rand!(rng::AbstractRNG,
        X::AbstractVector{T}, d::Random.SamplerTrivial{Simplex{T}}) where T
    @assert length(X) == d[].N
    randsimplex!(rng, X)
end
Base.rand(rng::AbstractRNG, d::Random.SamplerTrivial{Simplex{T}}) where T =
    rand!(rng, Vector{T}(undef, d[].N), d)

# %%
using BenchmarkTools

X = zeros(10)
d = Simplex(10)
@btime rand!($X, $d)

# %%
rand(Simplex(10))

# %%
using Random, LinearAlgebra, Statistics, Plots

function randsphere!(rng::AbstractRNG, X)
    randn!(X)
    X .*= √(length(X)) / norm(X)
end
randsphere(rng::AbstractRNG, N, T=Float64) = randsphere!(rng, Vector{T}(undef, N))
randsphere!(X) = randsphere!(Random.default_rng(), X)
randsphere(N, T=Float64) = randsphere!(Vector{T}(undef, N))

n = 10^4
A = randsphere(n)
@show mean(A), var(A)
histogram(A; norm=true, alpha=0.3, xlim=(-5, 5), label="randsphere($n)")
plot!(x->exp(-x^2/2)/√(2π), -5, 5; lw=2, label="standard normsl dist")

# %%

max_stars_repo_path: 0013/rand of uniform dist on simplex.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0013/rand of uniform dist on simplex.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0013/rand of uniform dist on simplex.jl

max_forks_repo_name: genkuroki/public


text:
module CustomCollections

export
    ConstVector,
    ConstDict,
    NullDict,
    CacheElement,
    AbstractIndexDict,
    IndexDict,
    CacheIndexDict,
    SegmentedVector,
    SegmentedBlockDiagonalMatrix,
    UnorderedPair

export
    foreach_with_extra_args,
    isdirty,
    segments,
    ranges

using TypeSortedCollections
using DocStringExtensions
using LinearAlgebra
using UnsafeArrays

using Base: @propagate_inbounds

include("foreach_with_extra_args.jl")
include("ConstVector.jl")
include("ConstDict.jl")
include("NullDict.jl")
include("CacheElement.jl")
include("IndexDict.jl")
include("SegmentedVector.jl")
include("SegmentedBlockDiagonalMatrix.jl")
include("UnorderedPair.jl")

end # module


text:
"""
    Planet(optimizer, eager::Bool)

Planet integrates a SAT solver (`PicoSAT.jl`) to find an activation pattern that maps a feasible input to an infeasible output.

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: hyperrectangle or hpolytope
3. Output: halfspace

# Return
`BasicResult`

# Method
Binary search of activations (0/1) and pruning by optimization. Our implementation is non eager.
- `optimizer` default `GLPKSolverMIP()`;
- `eager` default `false`;

# Property
Sound and complete.

# Reference
[R. Ehlers, "Formal Verification of Piece-Wise Linear Feed-Forward Neural Networks,"
in *International Symposium on Automated Technology for Verification and Analysis*, 2017.](https://arxiv.org/abs/1705.01320)

[https://github.com/progirep/planet](https://github.com/progirep/planet)
"""
@with_kw struct Planet
    optimizer::AbstractMathProgSolver
    eager::Bool                        = false
end

Planet(x::AbstractMathProgSolver) = Planet(optimizer = x)

function solve(solver::Planet, problem::Problem)
    @assert ~solver.eager "Eager implementation not supported yet"
    # Refine bounds. The bounds are values after activation
    status, bounds = tighten_bounds(problem, solver.optimizer)
    status == :Optimal || return CounterExampleResult(:SAT)
    ψ = init_ψ(problem.network, bounds)
    δ = PicoSAT.solve(ψ)
    opt = solver.optimizer
    # Main loop to compute the SAT problem
    while δ != :unsatisfiable
        status, conflict = elastic_filtering(problem, δ, bounds, opt)
        status == :Infeasible || return CounterExampleResult(:UNSAT, conflict)
        append!(ψ, Any[conflict])
        δ = PicoSAT.solve(ψ)
    end
    return CounterExampleResult(:SAT)
end

function init_ψ(nnet::Network, bounds::Vector{Hyperrectangle})
    ψ = Vector{Vector{Int64}}(undef, 0)
    index = 0
    for i in 1:length(bounds)-1
        before_act_bound = linear_transformation(nnet.layers[i], bounds[i])
        lower = low(before_act_bound)
        upper = high(before_act_bound)
        for j in 1:length(lower)
            index += 1
            lower[j] > 0 && push!(ψ, [index])
            upper[j] < 0 && push!(ψ, [-index])
            lower[j] <= 0 <= upper[j] && push!(ψ, [index, -index])
        end
    end
    return ψ
end

function elastic_filtering(problem::Problem, δ::Vector{Vector{Bool}}, bounds::Vector{Hyperrectangle}, optimizer::AbstractMathProgSolver)
    model = Model(solver = optimizer)
    neurons = init_neurons(model, problem.network)
    add_set_constraint!(model, problem.input, first(neurons))
    add_complementary_set_constraint!(model, problem.output, last(neurons))
    encode_Δ_lp!(model, problem.network, bounds, neurons)
    slack = encode_slack_lp!(model, problem.network, neurons, δ)
    min_sum!(model, slack)
    conflict = Vector{Int64}()
    act = get_activation(problem.network, bounds)
    while true
        status = solve(model, suppress_warnings = true)
        status == :Optimal || return (:Infeasible, conflict)
        (m, index) = max_slack(getvalue(slack), act)
        m > 0.0 || return (:Feasible, getvalue(neurons[1]))
        # activated neurons get a factor of (-1)
        coeff = δ[index[1]][index[2]] ? -1 : 1
        node = coeff * get_node_id(problem.network, index)
        push!(conflict, node)
        @constraint(model, slack[index[1]][index[2]] == 0.0)
    end
end

elastic_filtering(problem::Problem, list::Vector{Int64}, bounds::Vector{Hyperrectangle}, optimizer::AbstractMathProgSolver) = elastic_filtering(problem, get_assignment(problem.network, list), bounds, optimizer)

function max_slack(x::Vector{Vector{Float64}}, act)
    m = 0.0
    index = (0, 0)
    for i in 1:length(x)
        for j in 1:length(x[i])
            if x[i][j] > m && act[i][j] == 0 # Only return undetermined nodes
                m = x[i][j]
                index = (i, j)
            end
        end
    end
    return (m, index)
end

function tighten_bounds(problem::Problem, optimizer::AbstractMathProgSolver)
    bounds = get_bounds(problem)
    model = Model(solver = optimizer)
    neurons = init_neurons(model, problem.network)
    add_set_constraint!(model, problem.input, first(neurons))
    add_complementary_set_constraint!(model, problem.output, last(neurons))
    encode_Δ_lp!(model, problem.network, bounds, neurons)

    o = min_sum!(model, neurons)
    status = solve(model, suppress_warnings = true)
    status == :Optimal || return (:Infeasible, bounds)
    lower = getvalue(neurons)

    o = max_sum!(model, neurons)
    status = solve(model, suppress_warnings = true)
    status == :Optimal || return (:Infeasible, bounds)
    upper = getvalue(neurons)

    new_bounds = Vector{Hyperrectangle}(undef, length(neurons))
    for i in 1:length(neurons)
        new_bounds[i] = Hyperrectangle(low = lower[i], high = upper[i])
    end
    return (:Optimal, new_bounds)
end

function get_assignment(nnet::Network, list::Vector{Int64})
    p = Vector{Vector{Bool}}(undef, length(nnet.layers))
    n = 0
    for (i, layer) in enumerate(nnet.layers)
        p[i] = zeros(Bool, length(layer.bias))
        for j in 1:length(p[i])
            if list[n+j] > 0
                p[i][j] = true
            end
        end
        n += length(p[i])
    end
    return p
end

function get_node_id(nnet::Network, x::Tuple{Int64, Int64})
    n = 0
    for i in 1:x[1]-1
        n += length(nnet.layers[i].bias)
    end
    return n + x[2]
end

function get_node_id(nnet::Network, n::Int64)
    i = 0
    j = n
    while j > length(nnet.layers[i+1].bias)
        i += 1
        j = j - length(nnet.layers[i].bias)
    end
    return (i, j)
end

# function init_ψ(p_list::Vector{Int64})
#     ψ = Vector{Vector{Int64}}(undef, length(p_list))
#     for i in 1:length(p_list)
#         if p_list[i] == 0
#             ψ[i] = [i, -i]
#         else
#             ψ[i] = [p_list[i]*i]
#         end
#     end
#     return ψ
# end

# function get_list(p::Vector{Vector{Int64}})
#     list = Vector{Int64}()
#     for i in 1:length(p)
#         for j in 1:length(p[i])
#             append!(list, Int64[p[i][j]])
#         end
#     end
#     return list
# end

# To be implemented
# function infer_node_phases(nnet::Network, p::Vector{Vector{Int64}}, bounds::Vector{Hyperrectangle})
#     extra = Vector{Vector{Int64}}()
#     return extra
# end

# function encode_partial_assignment(model::Model, nnet::Network, p::Vector{Vector{Int64}}, neurons, slack::Bool)
#     if slack
#         slack_var = Vector{Vector{Variable}}(undef, length(nnet.layers))
#         sum_slack = 0.0
#     end

#     for (i, layer) in enumerate(nnet.layers)
#         (W, b, act) = (layer.weights, layer.bias, layer.activation)
#         before_act = W * neurons[i] + b
#         if slack
#             slack_var[i] = @variable(model, [1:length(b)])
#         end
#         for j in 1:length(b)
#             if p[i][j] != 0
#                 if slack
#                     sum_slack += slack_var[i][j]
#                     if p[i][j] == 1
#                         @constraint(model, neurons[i+1][j] == before_act[j] + slack_var[i][j])
#                         @constraint(model, before_act[j] + slack_var[i][j] >= 0.0)
#                     else
#                         @constraint(model, neurons[i+1][j] == 0.0)
#                         @constraint(model, 0.0 >= before_act[j] - slack_var[i][j])
#                     end
#                 else
#                     if p[i][j] == 1
#                         @constraint(model, neurons[i+1][j] <= before_act[j])
#                     else
#                         @constraint(model, neurons[i+1][j] == 0.0)
#                         @constraint(model, 0.0 >= before_act[j])
#                     end
#                 end
#             end
#         end
#     end
#     if slack
#         return slack_var, sum_slack
#     else
#         return nothing
#     end
# end

# function get_tight_clause(problem::Problem, p::Vector{Vector{Int64}}, bounds::Vector{Hyperrectangle})
#     model = Model(solver = optimizer)

#     neurons = init_neurons(solver, model, problem.network)
#     add_set_constraint!(model, problem.input, first(neurons))
#     add_complementary_output_constraint(model, problem.output, last(neurons))
#     encode_Δ_lp!(model, problem.network, bounds, neurons)
#     encode_partial_assignment(model, problem.network, p, neurons, false)

#     o = 0.0
#     for i in 1:length(p)
#         for j in 1:length(p[i])
#             if p[i][j] == 0
#                 o += neurons[i+1][j]
#             end
#         end
#     end

#     @objective(model, Min, o)
#     status = solve(model, suppress_warnings = true)
#     if status != :Optimal
#         return (:Infeasible, [])
#     end

#     v = getvalue(neurons)
#     tight = Vector{Int64}()
#     complete = :Complete
#     for i in 1:length(p)
#         for j in 1:length(p[i])
#             if p[i][j] == 0
#                 if v[i+1][j] > 0
#                     append!(tight, Any[get_node_id(problem.network, (i,j))])
#                 else
#                     complete = :Incomplete
#                 end
#             end
#         end
#     end
#     return (complete, tight)
# end


text:
# Generates a batch of plots

using TrajOptSOCPs

include("plotConstraintViolation.jl")
include("plotObjective.jl")
include("plotTrajectory.jl")
include("plotTrajectory3D.jl")
include("plotThrust.jl")

function batchPlot(trajStates, cM::constraintManager, nDim::Int64,
                            penalty::Float64 = 1.0, maxThrust::Float64 = 0.0,
                            maxAngleDeg::Float64 = 0.0)
    ptList = [getParseTrajectory(traj, nDim) for traj in trajStates]
    # println("ptList: ")
    # display(ptList)

    ptrajStart = ptList[1]
    ptrajLast = ptList[end]

    if nDim == 2
        pltTraj = plotTrajPos2D_Multiple(ptList)
        display(pltTraj)
    elseif nDim == 3
        pltTraj = plotTrajPos3D_Multiple(ptList)
        display(pltTraj)
    end

    pltCV, pltCV2 = plotConstraintViolation(cM, trajStates, penalty)
    display(pltCV)
    display(pltCV2)

    pltObj = plotObjective(costFun, trajStates)
    display(pltObj)

    # plotSVUTime_Simple(ptList[1])
    if nDim == 2
        plts, pltv, pltu = plotSVUTime_StartEnd([ptrajStart, ptrajLast])
    elseif nDim == 3
        plts, pltv, pltu = plotSVUTime3D_StartEnd([ptrajStart, ptrajLast])
    end

    pltMagT = plotMagThrust_StartEnd(ptrajStart.uList, ptrajLast.uList,
                                     maxThrust)
    display(pltMagT)

    pltAngleT = plotAngle3D_StartEnd(ptrajStart.uList, ptrajLast.uList,
                                     maxAngleDeg)
    display(pltAngleT)

    pltAngleTNoLine = plotAngle3D_StartEnd(ptrajStart.uList, ptrajLast.uList)
    display(pltAngleTNoLine)

    pltList = [pltTraj, pltCV, pltCV2, pltObj, plts, pltv, pltu,
                pltMagT, pltAngleT, pltAngleTNoLine]

    pltLabels = ["Trajectory", "ConstraintViolation", "DynamicsViolation",
                    "Objective", "PositionTime", "VelocityTime",
                    "ControlsTime", "ThrustMagTime", "ThrustAngleTime",
                    "ThrustAngleTimeNoLine"]

    # Return the result as a dictionary
    return Dict(zip(pltLabels, pltList))
end

function saveBulk(pltDict, header::String)

    for (index, value) in pltDict
        println("Saving: $index $value")
        savefig(value, header * index)
    end

end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
"""
    core_periphery_deg(g)

A simple degree-based core-periphery detection algorithm (see [Lip](http://arxiv.org/abs/1102.5511)).
Returns the vertex assignments (1 for core and 2 for periphery).
"""
function core_periphery_deg(g::Graph)
    degs = degree(g)
    p = sortperm(degs, rev=true)
    s = sum(degs) / 2.
    sbest = +Inf
    kbest = 0
    for k = 1:nv(g)-1
        s = s + k  - 1 - degree(g, p[k])
        if s < sbest
            sbest = s
            kbest = k
        end
    end
    c = 2 + zeros(Int, nv(g))
    c[p[1:kbest]] = 1
    c
end


text:
"""
    TypeJ(var_1::TypeK)

Lorem ipsum dolor sit amet.
"""
struct TypeJ

  var_1::TypeK

end


text:
# Objective
function JuMP.set_objective_function(m::RobustModel, f::JuMP.AbstractJuMPScalar)
    return m.uncertain_objective_function = f
end
function JuMP.set_objective_function(m::RobustModel, f::Real)
    return m.uncertain_objective_function = JuMP.GenericAffExpr{Float64,UncertainVariableRef}(f)
end
JuMP.objective_sense(model::RobustModel) = model.objectivesense
function JuMP.set_objective_sense(model::RobustModel, sense)
    return model.objectivesense = sense
end
JuMP.objective_function_type(model::RobustModel) = typeof(model.uncertain_objective_function)
JuMP.objective_function(model::RobustModel) = model.uncertain_objective_function
function JuMP.objective_function(model::RobustModel, FT::Type)
    # InexactError should be thrown, this is needed in `objective.jl`
    if !(model.uncertain_objective_function isa FT)
        throw(
            InexactError(
                :objective_function,
                FT,
                typeof(model.uncertain_objective_function),
            ),
        )
    end
    return model.uncertain_objective_function::FT
end


text:
#######################################################################################
# GPU-crutch solve method
# Makes the simplest possible method for GPU-compatibility
# Out of place only
#######################################################################################
struct GPUSimpleTsit5 <: DiffEqBase.AbstractODEAlgorithm end
export GPUSimpleTsit5

@muladd function DiffEqBase.solve(prob::ODEProblem,
                          alg::GPUSimpleTsit5;
                          dt = 0.1f0)
  @assert !isinplace(prob)
  u0 = prob.u0
  tspan = prob.tspan
  f = prob.f
  p = prob.p
  t = tspan[1]
  tf = prob.tspan[2]
  ts = tspan[1]:dt:tspan[2]
  us = MVector{length(ts),typeof(u0)}(undef)
  us[1] = u0
  u = u0
  k7 = f(u, p, t)

  cs, as, btildes, rs = _build_atsit5_caches(eltype(u0))
  c1, c2, c3, c4, c5, c6 = cs
  a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
  a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = as
  btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = btildes

  # FSAL
  for i in 2:length(ts)
      uprev = u; k1 = k7; t = ts[i]
      tmp = uprev+dt*a21*k1
      k2 = f(tmp, p, t+c1*dt)
      tmp = uprev+dt*(a31*k1+a32*k2)
      k3 = f(tmp, p, t+c2*dt)
      tmp = uprev+dt*(a41*k1+a42*k2+a43*k3)
      k4 = f(tmp, p, t+c3*dt)
      tmp = uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)
      k5 = f(tmp, p, t+c4*dt)
      tmp = uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)
      k6 = f(tmp, p, t+dt)
      u = uprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6)
      k7 = f(u, p, t+dt)
      us[i] = u
  end

  sol = DiffEqBase.build_solution(prob,alg,ts,SArray(us),
                                  k = nothing, destats = nothing,
                                  calculate_error = false)
  DiffEqBase.has_analytic(prob.f) && DiffEqBase.calculate_solution_errors!(sol;timeseries_errors=true,dense_errors=false)
  sol
end

#######################################################################################
# GPU-crutch solve method
# Makes the simplest possible adaptive method for GPU-compatibility
# Out of place only
#######################################################################################
struct GPUSimpleATsit5 end
export GPUSimpleATsit5

@muladd function DiffEqBase.solve(prob::ODEProblem,
                          alg::GPUSimpleATsit5;
                          dt = 0.1f0,saveat = nothing,
                          save_everystep = true,
                          abstol = 1f-6, reltol = 1f-3)
  @assert !isinplace(prob)
  u0 = prob.u0
  tspan = prob.tspan
  f = prob.f
  p = prob.p

  t = tspan[1]
  tf = prob.tspan[2]

  if saveat === nothing
    ts = Vector{eltype(dt)}(undef,1)
    ts[1] = prob.tspan[1]
    us = Vector{typeof(u0)}(undef,0)
    push!(us,recursivecopy(u0))
  else
    ts = saveat
    cur_t = 1
    us = MVector{length(ts),typeof(u0)}(undef)
    if prob.tspan[1] == ts[1]
      cur_t += 1
      us[1] = u0
    end
  end

  u = u0
  qold = qoldinit
  k7 = f(u, p, t)

  cs, as, btildes, rs = _build_atsit5_caches(eltype(u0))
  c1, c2, c3, c4, c5, c6 = cs
  a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
  a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = as
  btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = btildes

  # FSAL
  while t < tspan[2]
      uprev = u
      k1 = k7
      EEst = Inf

      while EEst > 1
        dt < 1e-14 && error("dt<dtmin")

        tmp = uprev+dt*a21*k1
        k2 = f(tmp, p, t+c1*dt)
        tmp = uprev+dt*(a31*k1+a32*k2)
        k3 = f(tmp, p, t+c2*dt)
        tmp = uprev+dt*(a41*k1+a42*k2+a43*k3)
        k4 = f(tmp, p, t+c3*dt)
        tmp = uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)
        k5 = f(tmp, p, t+c4*dt)
        tmp = uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)
        k6 = f(tmp, p, t+dt)
        u = uprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6)
        k7 = f(u, p, t+dt)

        tmp = dt*(btilde1*k1+btilde2*k2+btilde3*k3+btilde4*k4+
                     btilde5*k5+btilde6*k6+btilde7*k7)
        tmp = tmp./(abstol.+max.(abs.(uprev),abs.(u))*reltol)
        EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t)

        if iszero(EEst)
          q = inv(qmax)
        else
          @fastmath q11 = EEst^beta1
          @fastmath q = q11/(qold^beta2)
        end

        if EEst > 1
          dt = dt/min(inv(qmin),q11/gamma)
        else # EEst <= 1
          @fastmath q = max(inv(qmax),min(inv(qmin),q/gamma))
          qold = max(EEst,qoldinit)
          dtold = dt
          dt = dt/q #dtnew
          dt = min(abs(dt),abs(tf-t-dtold))
          told = t
          if (tf - t - dtold) < 1e-14
            t = tf
          else
            t += dtold
          end

          if saveat === nothing && save_everystep
            push!(us,recursivecopy(u))
            push!(ts,t)
          else saveat !== nothing
            while cur_t <= length(ts) && ts[cur_t] <= t
              savet = ts[cur_t]
              θ = (savet - told)/dtold
              b1θ, b2θ, b3θ, b4θ, b5θ, b6θ, b7θ = bθs(rs, θ)
              us[cur_t] = uprev + dtold*(
                   b1θ*k1 + b2θ*k2 + b3θ*k3 + b4θ*k4 + b5θ*k5 + b6θ*k6 + b7θ*k7)
              cur_t += 1
            end
          end

        end
      end
  end

  if saveat === nothing && !save_everystep
      push!(us,u)
      push!(ts,t)
  end
  sol = DiffEqBase.build_solution(prob,alg,ts,us,
                                  calculate_error = false)
  DiffEqBase.has_analytic(prob.f) && DiffEqBase.calculate_solution_errors!(sol;timeseries_errors=true,dense_errors=false)
  sol
end


text:
pcre_option =
    Either(
        # with_name(:MARK, "mark" => UInt32(0)),
        # with_name(:aftertext, "aftertext" => UInt32(0)),
        with_name(:DUPNAMES, "dupnames" => Base.PCRE.DUPNAMES),
        # with_name(:no_start_optimize, "no_start_optimize" => UInt32(0)),
        # with_name(:subject_literal, "subject_literal" => UInt32(0)),
        # "jitstack=256" => UInt32(0),
        with_name(:EXTENDED_MORE, "xx" => Base.PCRE.EXTENDED_MORE),
        with_name(:CASELESS, 'i' => Base.PCRE.CASELESS),
        with_name(:MULTILINE, 'm' => Base.PCRE.MULTILINE),
        with_name(:NO_AUTO_CAPTURE, 'n' => Base.PCRE.NO_AUTO_CAPTURE),
        with_name(:UNGREEDY, 'U' => Base.PCRE.UNGREEDY),
        with_name(:DUPNAMES, 'J' => Base.PCRE.DUPNAMES),
        with_name(:DOTALL, 's' => Base.PCRE.DOTALL),
        with_name(:EXTENDED, 'x' => Base.PCRE.EXTENDED),
        # 'g' => UInt32(0),
        with_name(:BINCODE, 'B' => UInt32(0)), # bincode
        with_name(:INFO, 'I' => UInt32(0)) # info
    );

splat_or(v) = (isempty(v) ? 0x00000000 : (|(v...)))::UInt32
pcre_options = Repeat(splat_or,map(IndexAt(1),Sequence(pcre_option,Optional(','))))

pcre_options_parser=map(IndexAt(2),Sequence(AtStart(),pcre_options,AtEnd()))

"""
    parse_options(options::AbstractString)

Return PCRE option mask parsed from `options`.

Parser for `flags` in [`@re_str`](@ref).

```jldoctest
julia> CombinedParsers.Regexp.pcre_options_parser
🗄 Sequence[2]
├─ ^ AtStart
├─ 🗄+? Sequence[1] |> Repeat |> map(splat_or)|0
│  ├─ |🗄... Either
│  │  ├─ dupnames  => 0x00000040 |> with_name(:DUPNAMES)
│  │  ├─ xx  => 0x01000000 |> with_name(:EXTENDED_MORE)
│  │  ├─ i  => 0x00000008 |> with_name(:CASELESS)
│  │  ├─ m  => 0x00000400 |> with_name(:MULTILINE)
│  │  ├─ n  => 0x00002000 |> with_name(:NO_AUTO_CAPTURE)
│  │  ├─ U  => 0x00040000 |> with_name(:UNGREEDY)
│  │  ├─ J  => 0x00000040 |> with_name(:DUPNAMES)
│  │  ├─ s  => 0x00000020 |> with_name(:DOTALL)
│  │  ├─ x  => 0x00000080 |> with_name(:EXTENDED)
│  │  ├─ B  => 0x00000000 |> with_name(:BINCODE)
│  │  └─ I  => 0x00000000 |> with_name(:INFO)
│  └─ ,? |missing
└─ \$ AtEnd
::UInt32

```
"""
function parse_options(options::AbstractString)
    flags = tryparse(pcre_options_parser,options)
    if flags === nothing
        throw(UnsupportedError("options $options"))
    else
        flags
    end
end

function print_opts(io,opts)
    if (opts & Base.PCRE.CASELESS ) != 0; print(io, 'i'); end
    if (opts & Base.PCRE.MULTILINE) != 0; print(io, 'm'); end
    if (opts & Base.PCRE.DOTALL   ) != 0; print(io, 's'); end
    if (opts & Base.PCRE.EXTENDED_MORE ) != 0; print(io, "xx"); end
    if (opts & Base.PCRE.EXTENDED ) != 0; print(io, 'x'); end
    if (opts & Base.PCRE.DUPNAMES   ) != 0; print(io, 'J'); end
    if (opts & Base.PCRE.UNGREEDY ) != 0; print(io, 'U'); end
end
options_string(flags) =
    let sio = IOBuffer()
        print_opts(sio, flags)
        String(take!(sio))
    end

"""
A lazy element transformation type (e.g. AbstractString),
`getindex` wraps elements in `with_options(flags,...)`.

With parsing options

TODO: make flags a transformation function?
"""
struct WithOptions{S} <: AbstractString
    x::S
    flags::UInt32
    function WithOptions(x,flags)
        @assert !isa(x,WithOptions)
        new{typeof(x)}(x,flags)
    end
end
export flags
flags(x::WithOptions) = x.flags
flags(x) = UInt32(0)
regex_string(x::WithOptions) = WithOptions(regex_escape(x.x),x.flags)

"""
    with_options(options::AbstractString,str::AbstractString)

"""
with_options(flags::UInt32,x::WithOptions) =
    with_options(flags|x.flags,x.x)
with_options(flags::UInt32,x) =
    flags == 0 ? x : WithOptions(x,flags)
with_options(options::AbstractString,str::AbstractString) =
    with_options(parse_options(options),str)
with_options(set_flags::UInt32, unset_flags::UInt32,x::AbstractString) =
    with_options(set_flags,x)
with_options(set_flags::UInt32, unset_flags::UInt32,x::WithOptions) =
    with_options(set_flags | ( x.flags & ~unset_flags ),x.x)

"""
    convert(::Type{AbstractToken},x::Union{AbstractString,Char})

A [`ConstantParser`](@ref) matching `x`.
"""
Base.convert(::Type{AbstractToken},x::WithOptions{Char}) =
    if !iszero(x.flags & Base.PCRE.CASELESS)
        CharIn(lowercase(x.x),uppercase(x.x))
    else
        convert(AbstractToken,x.x)
    end

import Base: Regex
Base.Regex(x::WithOptions) =
    Regex(x.x, options_string(x.flags))
Base.show(io::IO, x::WithOptions) =
    print(io,x.x)
Base.print(io::IO, x::WithOptions{Char}) =
    print(io,x.x)

Base.getindex(x::WithOptions,i::Integer) =
    with_options(x.flags,(getindex(x.x,i)))
Base.iterate(x::WithOptions{<:AbstractString}) =
    let n = iterate(x.x)
        n===nothing ? nothing : ( with_options(x.flags,n[1]),n[2] )
    end
Base.iterate(x::WithOptions{<:AbstractString},i::Integer) =
    let n = iterate(x.x,i)
        n===nothing ? nothing : ( with_options(x.flags,n[1]),n[2] )
    end

Base.isless(x::WithOptions{Char},y) = isless(x.x,y)
(==)(x::WithOptions{Char},y) = x.x==y
Base.SubString(x::WithOptions,start::Int,stop::Int) = with_options(x.flags,SubString(x.x,start,stop))
Base.length(x::WithOptions) =
    length(x.x)
Base.lastindex(x::WithOptions) =
    lastindex(x.x)
Base.firstindex(x::WithOptions) =
    firstindex(x.x)
Base.prevind(x::WithOptions,i::Int,n::Int) =
    prevind(x.x,i,n)
Base.nextind(x::WithOptions,i::Int,n::Int) =
    nextind(x.x,i,n)
Base.prevind(x::WithOptions,i::Int) =
    prevind(x.x,i)
Base.nextind(x::WithOptions,i::Int) =
    nextind(x.x,i)
Base.ncodeunits(x::WithOptions) =
    ncodeunits(x.x)

import ..CombinedParsers: MatchState
@inline function _iterate(p::WithOptions{Char}, sequence, till, posi, next_i, state::Nothing, nc=0)
    @inbounds sc,j=iterate(sequence,posi)
    if ismatch(p,sc)
        j, MatchState()
    else
        nothing
    end
end

import ..CombinedParsers: ismatch, _ismatch
function ismatch(c::WithOptions{Char},p)::Bool
    _ismatch(c.x,p)
end

function _ismatch(c,p::WithOptions{Char})::Bool
    if !iszero(p.flags & Base.PCRE.CASELESS)
        _ismatch(lowercase(c),lowercase(p.x))
    else
        _ismatch(c,p.x)
    end
end

function Base.convert(::Type{CombinedParser},x::Char)
    CharIn(x)
end

import Base: convert
function Base.convert(::Type{Char},y::WithOptions{Char})
    if !iszero(y.flags & Base.PCRE.CASELESS)
        lowercase(y.x)
    else
        y.x
    end
end
function Base.convert(::Type{CombinedParser},x::WithOptions{Char})
    if x.flags & Base.PCRE.CASELESS > 0
        CharIn(unique([lowercase(x.x),uppercase(x.x)])...)
    else
        CharIn(x.x)
    end
end

revert(x::WithOptions) =
    WithOptions(revert(x.x),x.flags)

export set_options, with_options, on_options, map
"""
A wrapper matching the inner parser on `with_options(set_flags, unset_flags, sequence)`.
"""
struct ParserOptions{P,S,T} <: WrappedParser{P,S,T}
    parser::P
    set_flags::UInt32
    unset_flags::UInt32
    ParserOptions(parser,set::UInt32,unset::UInt32) =
        new{typeof(parser),state_type(parser),result_type(parser)}(parser,set,unset)
end
deepmap_parser(f::Function,mem::AbstractDict,x::ParserOptions,a...; kw...) =
    get!(mem,x) do
        ParserOptions(
            deepmap_parser(f,mem,x.parser,a...; kw...),
            x.set_flags,x.unset_flags)
    end

set_options(set::UInt32,unset::UInt32,p) =
    ParserOptions(parser(p),set,unset)

set_options(set::UInt32,parser) =
    set_options(set,UInt32(0),parser)

function regex_prefix(x::ParserOptions)
    "(?" * options_string(x.set_flags) *
    if x.unset_flags!=0
        "-" * options_string(x.unset_flags)
    else
        ""
    end
end
function regex_suffix(x::ParserOptions)
    ")"
end
function print_constructor(io::IO,x::ParserOptions)
    print_constructor(io,x.parser)
    print(io, " |> set_options")
end

@inline Base.get(parser::ParserOptions, sequence, till, after, i, state) =
    get(parser.parser,
        with_options(parser.set_flags,parser.unset_flags,sequence),
        till, after, i, state)

@inline function _iterate(parser::ParserOptions, sequence, till, posi, next_i, state)
    _iterate(parser.parser,
             with_options(parser.set_flags,parser.unset_flags,sequence),
             till, posi, next_i, state)
end

"""
Lazy wrapper for a sequence, masking elements in `getindex` with MatchingNever if any of `flags` are not set.

TODO: make flags a filter function?
"""
struct FilterOptions{S}
    x::S
    flags::UInt32
end
import Base: Regex
Base.show(io::IO, x::FilterOptions) =
    print(io,"\"",x.x,"\"[",options_string(x.flags),"]")

Base.getindex(x::FilterOptions,i) =
    if_options(x.flags,x.x[i])
if_options(flags::UInt32,x::Char) =
    iszero(flags) ? x : MatchingNever{Char}()
if_options(flags::UInt32,x::WithOptions{Char}) =
    if (flags & x.flags) == flags
        x.x
    else
        MatchingNever{Char}()
    end

if_options(flags::UInt32,x::FilterOptions) =
    FilterOptions(x.x,flags | x.flags)
if_options(flags::UInt32,x) =
    iszero(flags) ? x : FilterOptions(x,flags)

Base.length(x::FilterOptions) =
    length(x.x)
Base.lastindex(x::FilterOptions) =
    lastindex(x.x)
Base.firstindex(x::FilterOptions) =
    firstindex(x.x)
Base.prevind(x::FilterOptions,i::Integer,n::Integer) =
    prevind(x.x,i,n)
Base.nextind(x::FilterOptions,i::Integer,n::Integer) =
    nextind(x.x,i,n)
Base.prevind(x::FilterOptions,i::Integer) =
    prevind(x.x,i)
Base.nextind(x::FilterOptions,i::Integer) =
    nextind(x.x,i)
Base.ncodeunits(x::FilterOptions) =
    ncodeunits(x.x)
Base.iterate(x::FilterOptions,a...) =
    let n = iterate(x.x,a...)
        n===nothing ? nothing : if_options(x.flags,tuple_pos(n)),tuple_state(n)
    end

Base.convert(::Type{Union{Char,CharIn}},x::FilterOptions{Char}) = x.x

function Base.in(x::FilterOptions,set)
    if !iszero(x.flags & Base.PCRE.CASELESS)
        lowercase(x.x) in set || uppercase(x.x) in set
    else
        x.x in set
    end
end

export on_options
"""
Parser wrapper sequence with `if_options`.
"""
struct OnOptionsParser{P,S,T} <: WrappedParser{P,S,T}
    parser::P
    flags::UInt32
    OnOptionsParser(parser,flags::UInt32) =
        new{typeof(parser),state_type(parser),result_type(parser)}(parser,flags)
end

function print_constructor(io::IO, x::OnOptionsParser)
    print_constructor(io,x.parser)
    print(io," |> on_options(\"$(options_string(x.flags))\")")
end

deepmap_parser(f::Function,mem::AbstractDict,x::OnOptionsParser,a...; kw...) =
    get!(mem,x) do
        OnOptionsParser(
            deepmap_parser(f,mem,x.parser,a...; kw...),
            x.flags)
    end

on_options(flags::Integer,p) =
    OnOptionsParser(parser(p),UInt32(flags))

@inline function _iterate(parser::OnOptionsParser, sequence, till, posi, next_i, state)
    _iterate(parser.parser,
             (if_options(parser.flags,sequence)), till, posi, next_i, state)
end

"""
    @pcre_testset

Define `@syntax pcre_test` and `@syntax pcre_tests` for parsing unit test output of the PCRE library.
The parser is used for testing `CombinedParser` and benchmarking against `Regex`.
"""
macro pcre_tests()
    esc(quote
        ## not handled in escaped_character, but backreference, if a capture with number (in decimal) is defined
        charparser = Either(
            Sequence(
                '\\',
                CombinedParsers.Regexp.integer_base(8,1,3)
            ) do v
            Char(v[2])
            end,
            CombinedParsers.Regexp.escaped_character,
            AnyChar())
        ## compile
        parse(charparser,"a");

        unescaped=map(Repeat_until(
            AnyChar(), Sequence(Repeat(' '),'\n');
            wrap=JoinSubstring)) do v
        ## join Chars after unescaping
        join(parse(Repeat(charparser),v))
        end;
        parse(unescaped,"abc\n");
        comment_or_empty = Repeat(
            JoinSubstring(Either(
                Sequence(
                    CombinedParsers.Regexp.at_linestart,
                    '#',Repeat_until(AnyChar(),'\n')),
                Sequence(
                    CombinedParsers.Regexp.at_linestart,
                    Repeat_until(
                        CombinedParsers.Regexp.whitespace_char,'\n')))));
        parse(comment_or_empty,"# test");

        #    @test parse(unescaped,"A\\123B\n") == "ASB"
        @syntax pcre_test = begin
        match_test = Sequence(Repeat1(' '),
                              :sequence => unescaped,
                              :expect => Repeat(Sequence(
                                  Repeat(' '),
                                  :i => Either(CombinedParsers.Regexp.integer,"MK"),':',
                                  Repeat(' '),
                                  :result => unescaped))
                              );

        Sequence(
            :pattern => after(CharIn("/'\""),Any) do s
            Repeat_until(
                AnyChar(),
                Sequence(3, NegativeLookbehind('\\'),
                         s, Repeat_until(
                             AnyChar(),
                             Sequence(Repeat(
                                 CombinedParsers.Regexp.whitespace_char), '\n'),
                             wrap=JoinSubstring)),
                true; wrap=JoinSubstring)
            end,
            :test => Repeat(match_test),
            :tests_nomatch => Optional(
                Sequence(
                    2,
                    Optional("\\= Expect no match",
                             Repeat_until(AnyChar(), '\n'; wrap=JoinSubstring)),
                    Repeat(Sequence(2,
                                    Repeat1(' '),
                                    unescaped,
                                    Optional(Sequence("No match",
                                                      Repeat_until(AnyChar(), '\n'; wrap=JoinSubstring)))
                                    ))))
        );
        end;

        function Base.show(io::IO, x::result_type(pcre_test))
        print(io, "Pattern: ")
        printstyled(io,"r(e)\"$(x.pattern[1])\"$(x.pattern[2])\n", color=:underline)
        println(io, "Test Examples:")
        for (i,t) in enumerate(x.test)
        println(io, "   $i. $(t.sequence)")
        end
        println(io, "Not Examples:")
        for (i,t) in enumerate(x.tests_nomatch)
        println(io, "   $i. $t")
        end
        end

        @syntax pcre_tests = Sequence(
            Repeat(Sequence(comment_or_empty,
                            pcre_test)),
            comment_or_empty,
            AtEnd());
        end)
end