text:
# Modifier and sampler types

"""
`AbstractSampler` type is responsible for picking samples (voxels or
pixels) during annealing process.

See also: [`UniformSampler`](@ref), [`InterfaceSampler`](@ref),
[`DPNSampler`](@ref).
"""
abstract type AbstractSampler end

"""
    UniformSampler()

Create a sampler which picks samples randomly with uniform
distribution.

See also: [`AbstractSampler`](@ref), [`InterfaceSampler`](@ref),
[`DPNSampler`](@ref).
"""
struct UniformSampler <: AbstractSampler end

"""
    InterfaceSampler()

Create a sampler which picks samples randomly from an interface
between phases.

See also: [`UniformSampler`](@ref), [`AbstractSampler`](@ref),
[`DPNSampler`](@ref).
"""
struct InterfaceSampler <: AbstractSampler end

# Here and later: DPN = different phase neighbors
struct DPNSampler <: AbstractSampler
    histogram :: Vector{Int}
    α         :: Float64
end

@doc raw"""
    DPNSampler(array, α)

Create a sampler which picks a random sample with $n$ different phase
neighbors with probability  $p \propto \alpha^n N(n)$ where $N(n)$ is
a number of pixels having $n$ different phase neighbors in the system.

See also: [`UniformSampler`](@ref), [`InterfaceSampler`](@ref),
[`AbstractSampler`](@ref).
"""
DPNSampler(array :: AbstractArray, α :: Float64) =
    DPNSampler(dpn_histogram(array), α)

"""
`AbstractModifier` type is responsible for how configuration update
steps are performed during an annealing process.

See also: [`Flipper`](@ref), [`Swapper`](@ref).
"""
abstract type AbstractModifier end

"""
    Flipper(sampler)

Create a modifier which changes phase of a sample taken from a binary
array.

See also: [`AbstractModifier`](@ref), [`Swapper`](@ref).
"""
struct Flipper{T <: AbstractSampler} <: AbstractModifier
    sampler :: T
end

"""
    Swapper(sampler)

Create a modifier which takes two samples having different phases and
swaps phases of these two samples.

See also: [`Flipper`](@ref), [`AbstractModifier`](@ref).
"""
struct Swapper{T <: AbstractSampler} <: AbstractModifier
    sampler :: T
end

struct BatchRollbackToken{T <: AbstractRollbackToken} <: AbstractRollbackToken
    tokens :: Vector{T}
end

"""
    BatchModifier(modifier, mods)

Create a modifier which performs `mods` modifications with a modifier
`modifier`. These modifications are then either accepted alltogether
or rejected alltogether.
"""
struct BatchModifier{T <: AbstractModifier} <: AbstractModifier
    modifier :: T
    mods     :: Int
end

# Interface

"""
    sample(array, sampler)

Get a random sample from the array. This function must return an
index, e.g. `CartesianIndex`.
"""
function sample end

"""
    update_pre!(array, index, sampler)

Update the state of sampler before the value `array[index]` is
changed. Should only be implemented for stateful samplers.
"""
function update_pre! end

"""
    update_post!(array, index, sampler)

Update the state of sampler after the value `array[index]` is
changed. Should only be implemented for stateful samplers.
"""
function update_post! end

"""
    modify!(array, modifier)

Randomly modify the array. This function must return a state used in
`reject!` to bring the array to a previous state if this modification
is rejected.
"""
function modify! end

"""
    reject!(array, modifier, state)

Bring the array back to the previous state.
"""
function reject! end

# Methods

# Stateless samples do not have to implement this
update_pre!( :: Any, :: Any, :: AbstractSampler) = nothing
update_post!(:: Any, :: Any, :: AbstractSampler) = nothing

sample(array :: AbstractArray, :: UniformSampler) =
    array |> CartesianIndices |> rand

function sample(array :: AbstractArray, :: InterfaceSampler)
    indices = CartesianIndices(array)

    while true
        rndidx = rand(indices)
        ray = RandomLineIterator(rndidx)
        start_phase = array[rndidx]

        cont(x :: CartesianIndex) =
            checkbounds(Bool, array, x) && array[x] == start_phase
        index = dropwhile(cont, ray) |> first

        if checkbounds(Bool, array, index)
            return index
        end
    end
end

# That horrible stateful sampler
function count_dpn(array :: AbstractArray{T, N},
                   idx   :: CartesianIndex{N}) where {T, N}
    indices = CartesianIndices(array)
    fidx, lidx = first(indices), last(indices)
    uidx = oneunit(idx)
    elt = array[idx]
    return sum(array[idx] != elt for idx in max(idx - uidx, fidx):min(idx + uidx, lidx))
end

function dpn_histogram(array :: AbstractArray)
    dims = ndims(array)
    hist = zeros(Int, 3^dims)

    for idx in CartesianIndices(array)
        different = count_dpn(array, idx)
        hist[different+1] += 1
    end

    return hist
end

function update_pdf(histogram :: Vector{Int}, α :: Float64)
    prob = [p*α^n for (p, n) in zip(histogram, countfrom(0))]
    prob = prob ./ sum(prob)
    return prob
end

function neighbors_for_update(histogram :: Vector{Int}, α :: Float64)
    prob = update_pdf(histogram, α)
    for idx in 2:length(prob)
        prob[idx] += prob[idx-1]
    end

    r = rand(Float64)
    return findfirst(x -> r < x, prob) - 1
end

function update_histogram!(array   :: AbstractArray{T, N},
                           index   :: CartesianIndex{N},
                           sampler :: DPNSampler,
                           fn      :: Function) where {T, N}
    histogram = sampler.histogram

    indices = CartesianIndices(array)
    fidx, lidx = first(indices), last(indices)
    uidx = oneunit(fidx)

    for idx in max(index - uidx, fidx):min(index + uidx, lidx)
        n_neighbor = count_dpn(array, idx)
        histogram[n_neighbor + 1] = fn(histogram[n_neighbor + 1])
    end
end

update_pre!(array   :: AbstractArray{T, N},
            index   :: CartesianIndex{N},
            sampler :: DPNSampler) where {T, N} =
                update_histogram!(array, index, sampler, x -> x - 1)

update_post!(array   :: AbstractArray{T, N},
             index   :: CartesianIndex{N},
             sampler :: DPNSampler) where {T, N} =
                 update_histogram!(array, index, sampler, x -> x + 1)

function sample(array :: AbstractArray, sampler :: DPNSampler)
    histogram = sampler.histogram
    α         = sampler.α

    n = neighbors_for_update(histogram, α)
    indices = CartesianIndices(array)
    @assert histogram[n + 1] > 0

    while true
        index = rand(indices)
        if n == count_dpn(array, index)
            return index
        end
    end
end

# Update correlation functions and sampler state
function update_all!(array   :: AbstractArray{T, N},
                     sampler :: AbstractSampler,
                     index   :: CartesianIndex{N},
                     val) where {T, N}
    update_pre!(array, index, sampler)
    token = update_corrfns!(array, val, index)
    update_post!(array, index, sampler)

    return token
end

# Rollback correlation functions and sampler state
function reject_all!(array   :: AbstractArray{T, N},
                     sampler :: AbstractSampler,
                     token   :: AbstractRollbackToken) where {T, N}
    update_pre!(array, token.index, sampler)
    rollback!(array, token)
    update_post!(array, token.index, sampler)

    return nothing
end

function modify!(array :: AbstractArray, modifier :: Flipper)
    sampler = modifier.sampler
    index   = sample(array, modifier.sampler)

    return update_all!(array, sampler, index, 1 - array[index])
end

reject!(array    :: AbstractArray,
        modifier :: Flipper,
        state    :: AbstractRollbackToken) =
            reject_all!(array, modifier.sampler, state)

function modify!(array :: AbstractArray, modifier :: Swapper)
    sampler = modifier.sampler
    index1 = sample(array, sampler)
    val1 = array[index1]

    while true
        index2 = sample(array, modifier.sampler)
        val2 = array[index2]
        if val1 ≠ val2
            token1 = update_all!(array, sampler, index1, val2)
            token2 = update_all!(array, sampler, index2, val1)
            return token1, token2
        end
    end
end

function reject!(array    :: AbstractArray,
                 modifier :: Swapper,
                 state    :: NTuple{2, AbstractRollbackToken})
    token1, token2 = state
    sampler = modifier.sampler

    reject_all!(array, sampler, token2)
    reject_all!(array, sampler, token1)
end

function modify!(array :: AbstractArray, modifier :: BatchModifier)
    m = modifier.modifier
    n = modifier.mods

    tokens = [modify!(array, m) for k in 1:n]
    return BatchRollbackToken{tokens |> first |> typeof}(tokens)
end

function reject!(array    :: AbstractArray,
                 modifier :: BatchModifier,
                 state    :: BatchRollbackToken)
    foreach(reverse(state.tokens)) do token
        reject!(array, modifier.modifier, token)
    end

    return nothing
end


text:
#name: Template
#description: Following template calculates number of cells in table
#language: julia
#tags: template, demo
#sample: cars.csv
#input: dataframe table [Data table]
#output: int count [Number of cells in table]
count = size(table)[1] * size(table)[2]

ext: jl

max_stars_repo_name: SoftwareCountry/datagrok-ai-public

max_issues_repo_name: SoftwareCountry/datagrok-ai-public

max_forks_repo_name: SoftwareCountry/datagrok-ai-public


text:
using JSON, HTTP
using HTTP.URIs: URI
import HTTP: get
import Base: write

export write, query

"""
    InfluxServer

Object that represents an InfluxDB server endpoint, with optional basic authentication
"""
struct InfluxServer
    # HTTP API endpoints
    endpoint::URI

    # Optional authentication stuffage
    username::Union{String,Nothing}
    password::Union{String,Nothing}

    # Build a server object that we can use in queries from now on
    function InfluxServer(address::AbstractString; username=nothing, password=nothing)
        # If there wasn't a schema defined (we only recognize http/https), default to http
        if match(r"^https?://", address) == nothing
            uri = URI("http://$address")
        else
            uri = URI(address)
        end

        # If we didn't get an explicit port, default to 8086
        if uri.port == 0
            uri =  URI(uri.scheme, uri.host, 8086, uri.path)
        end

        # URIs are the new hotness
        return new(uri, username, password)
    end
end

# Add authentication to a query dict, if we need to
function authenticate!(server::InfluxServer, query::Dict)
    if server.username != nothing && server.password != nothing
        query["u"] = server.username
        query["p"] = server.password
    end
end

function query(server::InfluxServer, query_data::Dict; type::Symbol = :get)
    # Add authentication to a query dict
    authenticate!(server, query_data)

    if type == :get
        response = HTTP.get(string(server.endpoint, "/query"); query=query_data)
    else
        response = HTTP.post(string(server.endpoint, "/query"); query=query_data)
    end
    response_data = String(response.body)
    if response.status != 200
        error(response_data)
    end

    # Grab JSON data
    response = JSON.parse(response_data)

    if !haskey(response, "results") || !haskey(response["results"][1], "series")
        return
    end

    # Convert the JSON object's values into a DataFrame
    function series_df(series_dict)
        df = DataFrame()
        for name_idx in 1:length(series_dict["columns"])
            df[!,Symbol(series_dict["columns"][name_idx])] = [x[name_idx] for x in series_dict["values"]]
        end
        return df
    end

    # Return a DataFrame for each measurement we requested
    return series_df.(response["results"][1]["series"])
end

# Shortcut for the common case of just a single "q" data element
query(server::InfluxServer, data::AbstractString; kwargs...) = query(server, Dict("q" => data); kwargs...)

# Shortcut for Query objects
query(server::InfluxServer, db::AbstractString, q::Query) = query(server, Dict("q" => payload(q), "db" => db))

function write(server::InfluxServer, db::AbstractString, query_data::Dict, payload::AbstractString)
    # Add authentication to a data dict
    authenticate!(server, query_data)

    response = HTTP.post(string(server.endpoint, "/write"), [], payload; query=query_data)
    response_data = String(response.body)
    if response.status != 204
        error(response_data)
    end
    return nothing
end

function write(server::InfluxServer, db::AbstractString, measurements::Vector{Measurement}; precision = :nanosecond)
    if isempty(measurements)
        return
    end

    # Start by building our query dict, pointing at a particular database and timestamp precision
    prec_dict = Dict(
        :nanosecond => "ns", :ns => "ns", "ns" => "ns",
        :microsecond => "us", :us => "us", "us" => "us",
        :millisecond => "ms", :ms => "ms", "ms" => "ms",
        :second => "s", :s => "s", "s" => "s",
        :minute => "m", :m => "m", "m" => "m",
        :hour => "h", :h => "h", "h" => "h",
    )
    query_data = Dict(
        "db"=>db,
        "precision"=>prec_dict[precision],
    )

    # Next, build payload
    buff = IOBuffer()
    for m in measurements
        payload!(buff, m)
        write(buff, "\n")
    end

    return write(server, db, query_data, String(take!(buff)))
end


text:
# # This is a title

# [This is a description](@ref custom_id)
# # [This is another title](@id custom_id)

# This isn't a description

x = 1+1


text:
using Plots
using GeometryTypes
gr()t
plot()

### SANITY CHECK
x = 1:10; y = rand(10,2) # 2 columns means two lines
plot(x,y)

#rect = HyperRectangle(Vec(0.5, 0.5), Vec(1., 1))

#module foo
#    struct Rect
#        NE::GeometryTypes.Point2f0
#        SW::GeometryTypes.Point2f0
#    end
#end

@recipe function f(r::HyperRectangle)
    points = decompose(Point{2,Float64}, r)
    rectpoints = points[[1,2,4,3],:]
    xs = [pt[1] for pt in rectpoints];
    ys = [pt[2] for pt in rectpoints];
    seriestype := :shape
	s = Shape(xs[:], ys[:])
end
sizePlot = (400,400)
r = HyperRectangle(Vec(1,1), Vec(2,4))
fig = plot(opacity=0.5, legend=false, size=sizePlot, yaxis=( (0,10), 0:1:10), xaxis=( (0,10), 0:1:10))
plot(r, opacity=0.5, legend=false, size=sizePlot, yaxis=( (0,10), 0:1:10), xaxis=( (0,10), 0:1:10))
plot!(r)


text:
#==
This is the modified ZDT example taken from
Esmaile Khorram, K Khaledian, and M Khaledyan.
A numerical method for
constructing the pareto front of multi-objective optimization problems.
Journal of Computational and Applied Mathematics, 261:158–171, 2014.

This is an example of a nonconvex Pareto front, where the weighted sum method
would not work.
==#
using JuMP, MultiJuMP, Ipopt

m = multi_model(solver=IpoptSolver())
n = 30

l = -ones(n); l[1] = 0
u = ones(n)
@variable(m, l[i] <= x[i=1:n] <= u[i])
@NLexpression(m, f1, x[1])
@NLexpression(m, g, 1 + 9*sum{x[j]^2, j=2:n}/(n-1))
@NLexpression(m, h, 1 - (f1/g)^2)
@NLexpression(m, f2, g*h)

obj1 = SingleObjective(f1, sense = :Min)

obj2 = SingleObjective(f2, sense = :Min)

multim = get_multidata(m)
multim.objectives = [obj1, obj2]
multim.pointsperdim = 30
solve(m, method = NBI())

plot(multim)


text:
module NUTSGibbs

using LinearAlgebra, Statistics
using TransformVariables, LogDensityProblems, DynamicHMC, ForwardDiff, Parameters
using Random: Random

import DynamicHMC: NUTS_init, NUTS_transition, get_final_ϵ, adapting_ϵ, adapt_stepsize
import LogDensityProblems: ADgradient

abstract type AbstractProblem end

export NUTSGibbsSampler

mutable struct NUTSGibbsSampler{Tp<:AbstractProblem, Tv, Tf, TA, TApar, Tfunc<:Function}
    iter::Int
    warmup::Int
    nuts::NUTS{Tv,Tf}
    problem::Tp
    transformation::Tfunc
    eps::Float64
    A_params::TApar
    A::TA
    sample::Vector{NUTS_Transition{Tv,Tf}}
    tune_steps::Vector{Int}
end

function NUTSGibbsSampler(nuts::NUTS{Tv,Tf}, problem, transformation, eps, warmup, A_params, A) where {Tv,Tf}
    sample = Vector{NUTS_Transition{Tv,Tf}}(undef, 0)
    tune_steps = floor.(Int, [0.15, 0.2, 0.3, 0.5, 0.9].*warmup)
    @info tune_steps
    NUTSGibbsSampler(0, warmup, nuts, problem, transformation, eps, A_params, A, sample, tune_steps)
end

function init(problem::Tp, transformation::Tf, warmup::Int) where {Tp <: AbstractProblem, Tf <: Function}
    P = TransformedLogDensity(transformation(problem), problem)
    nuts = nuts_try_init(ADgradient(:ForwardDiff, P))
    A_params, A = adapting_ϵ(nuts.ϵ)
    NUTSGibbsSampler(nuts, problem, transformation, nuts.ϵ, warmup, A_params, A)
end

function nuts_try_init(P)
    try
        nuts = NUTS_init(Random.GLOBAL_RNG, P; report=ReportSilent())
        @info "NUTS initial stepsize eps=$(nuts.ϵ)!"
        return nuts
    catch ex
        nuts = NUTS_init(Random.GLOBAL_RNG, ϵ=0.1, P)
        @warn "NUTS using fixed step size (eps=$(nuts.ϵ))!", ex
        return nuts
    end
end

function tune_cov(nuts::NUTS, A, sample; regularize=5.0)
    @unpack rng, H, max_depth, report = nuts

    Σ = DynamicHMC.sample_cov(sample)
    Σ += (DynamicHMC.UniformScaling(max(1e-3, median(diag(Σ))))-Σ) * regularize/length(sample)
    κ = DynamicHMC.GaussianKE(Symmetric(Σ))

    NUTS(rng, DynamicHMC.Hamiltonian(H.ℓ, κ), sample[end].q, get_final_ϵ(A), max_depth, report)
end

function step(x::NUTSGibbsSampler, current_pos)
    @unpack rng, H, q, ϵ, max_depth, report = x.nuts

    res = current_pos
    ac = 0.0
    x.iter += 1

    try
        pos = first(inverse.(Ref(x.transformation(x.problem)), [current_pos]))
        trans = NUTS_transition(rng, H, pos, x.eps, max_depth)
        res = first(transform.(Ref(x.transformation(x.problem)), [get_position(trans)]))
        ac = get_acceptance_rate(trans)

        if x.iter <= x.warmup
            x.A = adapt_stepsize(x.A_params, x.A, trans.a)
            x.eps = get_final_ϵ(x.A)

            push!(x.sample, trans)

            if x.iter ∈ x.tune_steps ##[100, 150, 250, 450, 850]
                x.nuts = tune_cov(x.nuts, x.A, x.sample)
                resize!(x.sample, 0)
            end

            if x.iter == x.warmup
                @info "Finish adaptation"
                @info "Final stepsize: $(x.eps)"
                @info "Mean acceptance rate at current transition: $(ac)"
            end
        end
    catch ex
        @warn "Exception in NUTS_transition: ", ex
    end

    return res, ac
end

end


text:
# Fetch a bunch of precomputed Hadamard matrices from neilsloane.com/hadamard/
# and convert them into BitArrays.  Then print the sizes and bitarray 64-bit
# chunks to give us a compact representation of the data for inlining.

const sloane = "http://neilsloane.com/hadamard/had." # URL prefix

# There are multiple inequivalent Hadamard matrices for many sizes;
# we just pick one, semi-arbitrarily.

const files = [
# powers of 2 (and 1) are handled specially
#               1 => "1",
#               2 => "2",
               12 => "12",
               20 => "20.1",
               28 => "28.pal2",
               36 => "36.pal2",
               44 => "44.pal",
               52 => "52.will",
               60 => "60.pal",
               68 => "68.pal",
               72 => "72.pal",
               76 => "76.pal2",
               84 => "84.pal",
               88 => "88.tpal",
               92 => "92.will",
               96 => "96.tpal",
               100 => "100.will",
               104 => "104.pal",
               108 => "108.pal",
               116 => "116.will",
               120 => "120.tpal",
               124 => "124.pal2",
               132 => "132.pal",
               136 => "136.tpal",
               140 => "140.pal",
               148 => "148.pal2",
               152 => "152.pal",
               156 => "156.will",
               164 => "164.pal",
               168 => "168.pal",
               172 => "172.will",
               180 => "180.pal",
               184 => "184.twill",
               188 => "188.tur",
               196 => "196.pal2",
               204 => "204.pal2",
               212 => "212.pal",
               216 => "216.tpal",
               220 => "220.pal2",
               228 => "228.pal",
               232 => "232.twill",
               236 => "236.od",
               244 => "244.will",
               248 => "248.twill",
               252 => "252.pal",
               ]

# the inverse of the printsigns function: parse +-+-... data (Sloane's format)
# into a bitarray.
function parsesigns(io::IO)
    lines = BitVector[]
    linelen = -1
    for s in eachline(io)
        s = chomp(s)
        line = BitVector(length(s))
        for (i,c) in enumerate(s)
            line[i] = c=='+' ? true : c=='-' ? false : error("invalid char $c")
        end
        if linelen == -1
            linelen = length(line)
        elseif linelen != length(line)
            error("line ", length(lines)+1, " has length ", length(line),
                  " != ", linelen)
        end
        push!(lines, line)
        length(lines) == linelen && break
    end
    B = BitArray(length(lines), linelen)
    for i = 1:length(lines)
        for j = 1:linelen
            B[i,j] = lines[i][j]
        end
    end
    B
end
parsesigns(s::String) = parsesigns(IOBuffer(s))

open("cache.dat", "w") do io
    f = tempname()
    for k in sort!(collect(keys(files)))
        println("FETCHING order $k...\n")
        B = try
            download(sloane*files[k]*".txt", f)
            open(parsesigns, f, "r")
        catch
            println("# ERROR for size $k: ", sloane*files[k]*".txt", " -> ", f)
            rethrow()
        end
        size(B) != (k,k) && error("size $k mismatch for ", files[k])
        write(io, hton(int64(k)), map(hton, B.chunks))
    end

    # size 428 is in a different format, grr
    k = 428
    println("FETCHING order $k...\n")
    B = try
        download(sloane*"428.txt", f)
        s = readall(f)
        s = replace(s, r" +1", "+")
        s = replace(s, r" *-1", "-")
        parsesigns(s)
    catch
        println("# ERROR for size $k: ", sloane*"428.txt", " -> ", f)
        rethrow()
    end
    size(B) != (k,k) && error("size $k mismatch for 428.txt")
    write(io, hton(int64(k)), map(hton, B.chunks))

    # others not on Sloane's page can be found at http://www.iasri.res.in/webhadamard/
    #= ... the site seems to be broken (truncated output) at the moment ...
    for k in [260,268,276,284,288,292,300,308,316,324,332,340,348,356,364,372,380,388,396,404,412,420,436,444,452,460,468,476,484,492,500,508,516,520,524,532,536,540,548,552,556,564,568,572,576,580,584,588,596,600,604,612,616,620,628,632,636,644,648,652,660,664,676,680,684,692,696,700,708,712,724,728,732,740,744,748,756,760,764,772,776,780,788,792,796,804,808,812,820,824,828,836,840,844,852,860,868,872,884,888,900,904,908,916,920,924,932,936,940,948,952,956,964,968,972,980,984,988,996,1000]
        B = try
            download("http://www.iasri.res.in/webhadamard/Validate.dll?MfcISAPICommand=Hadamard&Order=$k&Normalization=3&output=2", f)
            s = readall(f)
            s = replace(s, r"</p><p> *", "\n")
            s = replace(s, r"<html>.*<p> *"s, "")
            parsesigns(s)
        catch
            println("# ERROR for size $k from www.iasri.res.in/webhadamard")
            rethrow()
        end
        size(B) != (k,k) && error("size $k mismatch for webhadamard")
    end
    write(io, hton(int64(k)), map(hton, B.chunks))
    =#
end


text:
#Sets
set_generators = 1:2
set_nodes = 1:3
set_nodes_ref=3
set_nodes_noref=1:2
set_nodes_g= 1:2
set_nodes_d= 2:3
set_demands = 1:2

#Parameters
generators=[100      12;
            80       20]

demands=[100  40;
         50   35]

capacity=1
varcost=2
load_node=1
utility_node=2

network=[ 1  2  100 500;
          1  3   40 500;
          2  3  100 500;
          2  1  100 500;
          3  1  100 500;
          3  2  100 500]

start_node=network[:,1]
end_node=network[:,2]
F_max=network[:,3]
B=network[:,4]

n_node=3
n_link=6

set_links=1:n_link
#create a n*m Tuple Array with ‘undefined’ contents
links = Array{Tuple{Int64, Int64},2}(undef, n_link, 1)

for i=1:n_link
  links[i] = (start_node[i], end_node[i])
end

links_rev = Array{Tuple{Int64, Int64},2}(undef, n_link, 1)

for i=1:n_link
  links_rev[i] = (end_node[i], start_node[i])
end

F_max_dict=Dict()
B_dict=Dict()

for i=1:n_link
  F_max_dict[(start_node[i], end_node[i])] = F_max[i]
  B_dict[(start_node[i], end_node[i])] = B[i]
end

MapG=[(1,1), (2,2), (0,0)]
MapD=[(1,2),(2,3),(0,0)]


text:
#= ======================
   LINEAR REGRESSOR (OLS)
   ====================== =#

"""
Standard linear regression model with objective function

``|Xθ - y|₂²/2``

## Parameters

$TYPEDFIELDS

$(example_docstring("LinearRegressor"))
"""
@with_kw_noshow mutable struct LinearRegressor <: MMI.Deterministic
    "whether to fit the intercept or not."
    fit_intercept::Bool    = true
    "type of solver to use (if `nothing` the default is used). The solver is
    Cholesky by default but can be Conjugate-Gradient as well. See `?Analytical`
    for more information."
    solver::Option{Solver} = nothing
end

glr(m::LinearRegressor) = LinearRegression(fit_intercept=m.fit_intercept)

descr(::Type{LinearRegressor}) = "Regression with objective function ``|Xθ - y|₂²/2``."

#= ===============
   RIDGE REGRESSOR
   =============== =#

"""
Ridge regression model with objective function

``|Xθ - y|₂²/2 + n⋅λ|θ|₂²/2``

where ``n`` is the number of samples `size(X, 1)`.
With `scale_penalty_with_samples = false` the objective function is
``|Xθ - y|₂²/2 + λ|θ|₂²/2``.

## Parameters

$TYPEDFIELDS

$(example_docstring("RidgeRegressor"))
"""
@with_kw_noshow mutable struct RidgeRegressor <: MMI.Deterministic
    "strength of the L2 regularisation."
    lambda::Real             = 1.0
    "whether to fit the intercept or not."
    fit_intercept::Bool      = true
    "whether to penalize the intercept."
    penalize_intercept::Bool = false
    "whether to scale the penalty with the number of samples."
    scale_penalty_with_samples::Bool = true
    "type of solver to use (if `nothing` the default is used). The
     solver is Cholesky by default but can be Conjugate-Gradient as
     well. See `?Analytical` for more information."
    solver::Option{Solver}   = nothing
end

glr(m::RidgeRegressor) =
    RidgeRegression(m.lambda,
                    fit_intercept=m.fit_intercept,
                    penalize_intercept=m.penalize_intercept,
                    scale_penalty_with_samples=m.scale_penalty_with_samples)

descr(::Type{RidgeRegressor}) = "Regression with objective function ``|Xθ - y|₂²/2 + λ|θ|₂²/2``."

#= ===============
   LASSO REGRESSOR
   =============== =#

"""
Lasso regression model with objective function

``|Xθ - y|₂²/2 + n⋅λ|θ|₁``

where ``n`` is the number of samples `size(X, 1)`.
With `scale_penalty_with_samples = false` the objective function is
``|Xθ - y|₂²/2 + λ|θ|₁``

## Parameters

$TYPEDFIELDS

$(example_docstring("LassoRegressor"))
"""
@with_kw_noshow mutable struct LassoRegressor <: MMI.Deterministic
    "strength of the L1 regularisation."
    lambda::Real             = 1.0
    "whether to fit the intercept or not."
    fit_intercept::Bool      = true
    "whether to penalize the intercept."
    penalize_intercept::Bool = false
    "whether to scale the penalty with the number of samples."
    scale_penalty_with_samples::Bool = true
    "type of solver to use (if `nothing` the default is used). Either `FISTA` or
    `ISTA` can be used (proximal methods, with/without acceleration)."
    solver::Option{Solver}   = nothing
end

glr(m::LassoRegressor) =
    LassoRegression(m.lambda,
                    fit_intercept=m.fit_intercept,
                    penalize_intercept=m.penalize_intercept,
                    scale_penalty_with_samples=m.scale_penalty_with_samples)

descr(::Type{LassoRegressor}) = "Regression with objective function ``|Xθ - y|₂²/2 + λ|θ|₁``."

#= =====================
   ELASTIC NET REGRESSOR
   ===================== =#

"""
Elastic net regression model with objective function

``|Xθ - y|₂²/2 + n⋅λ|θ|₂²/2 + n⋅γ|θ|₁``

where ``n`` is the number of samples `size(X, 1)`.
With `scale_penalty_with_samples = false` the objective function is
``|Xθ - y|₂²/2 + λ|θ|₂²/2 + γ|θ|₁``

## Parameters

$TYPEDFIELDS

$(example_docstring("ElasticNetRegressor"))
"""
@with_kw_noshow mutable struct ElasticNetRegressor <: MMI.Deterministic
    "strength of the L2 regularisation."
    lambda::Real             = 1.0
    "strength of the L1 regularisation."
    gamma::Real              = 0.0
    "whether to fit the intercept or not."
    fit_intercept::Bool      = true
    "whether to penalize the intercept."
    penalize_intercept::Bool = false
    "whether to scale the penalty with the number of samples."
    scale_penalty_with_samples::Bool = true
    "type of solver to use (if `nothing` the default is used). Either `FISTA` or
    `ISTA` can be used (proximal methods, with/without acceleration)."
    solver::Option{Solver}   = nothing
end

glr(m::ElasticNetRegressor) =
    ElasticNetRegression(m.lambda, m.gamma,
                         fit_intercept=m.fit_intercept,
                         penalize_intercept=m.penalize_intercept,
                         scale_penalty_with_samples=m.scale_penalty_with_samples)

descr(::Type{ElasticNetRegressor}) = "Regression with objective function ``|Xθ - y|₂²/2 + λ|θ|₂²/2 + γ|θ|₁``."

#= ==========================
   ROBUST REGRESSOR (General)
   ========================== =#

"""
Robust regression model with objective function

``∑ρ(Xθ - y) + n⋅λ|θ|₂² + n⋅γ|θ|₁``

where ``ρ`` is a robust loss function (e.g. the Huber function) and
``n`` is the number of samples `size(X, 1)`.
With `scale_penalty_with_samples = false` the objective function is
``∑ρ(Xθ - y) + λ|θ|₂² + γ|θ|₁``.

## Parameters

$TYPEDFIELDS

$(example_docstring("RobustRegressor"))
"""
@with_kw_noshow mutable struct RobustRegressor <: MMI.Deterministic
    "the type of robust loss to use (see `HuberRho`, `TalwarRho`, ...)"
    rho::RobustRho           = HuberRho(0.1)
    "strength of the regulariser if `penalty` is `:l2` or `:l1`.
    Strength of the L2 regulariser if `penalty` is `:en`."
    lambda::Real             = 1.0
    "strength of the L1 regulariser if `penalty` is `:en`."
    gamma::Real              = 0.0
    "the penalty to use, either `:l2`, `:l1`, `:en` (elastic net) or `:none`."
    penalty::SymStr          = :l2
    "whether to fit the intercept or not."
    fit_intercept::Bool      = true
    "whether to penalize the intercept."
    penalize_intercept::Bool = false
    "whether to scale the penalty with the number of samples."
    scale_penalty_with_samples::Bool = true
    "type of solver to use, default if `nothing`."
    solver::Option{Solver}   = nothing
end

glr(m::RobustRegressor) =
    RobustRegression(m.rho, m.lambda, m.gamma;
                     penalty=Symbol(m.penalty),
                     fit_intercept=m.fit_intercept,
                     penalize_intercept=m.penalize_intercept,
                     scale_penalty_with_samples=m.scale_penalty_with_samples)

descr(::Type{RobustRegressor}) = "Robust regression with objective ``∑ρ(Xθ - y) + λ|θ|₂² + γ|θ|₁`` for a given robust `ρ`."

#= ===============
   HUBER REGRESSOR
   =============== =#

"""
Huber Regression is the same as `RobustRegressor` but with the robust loss
set to `HuberRho`.

## Parameters

$TYPEDFIELDS

$(example_docstring("HuberRegressor"))
"""
@with_kw_noshow mutable struct HuberRegressor <: MMI.Deterministic
    "parametrises the `HuberRho` function (radius of the ball within which the loss
is a quadratic loss)"
    delta::Real              = 0.5
    "strength of the regulariser if `penalty` is `:l2` or `:l1`.
    Strength of the L2 regulariser if `penalty` is `:en`."
    lambda::Real             = 1.0
    "strength of the L1 regulariser if `penalty` is `:en`."
    gamma::Real              = 0.0
    "the penalty to use, either `:l2`, `:l1`, `:en` (elastic net) or `:none`."
    penalty::SymStr          = :l2
    "whether to fit the intercept or not."
    fit_intercept::Bool      = true
    "whether to penalize the intercept."
    penalize_intercept::Bool = false
    "whether to scale the penalty with the number of samples."
    scale_penalty_with_samples::Bool = true
    "type of solver to use, default if `nothing`."
    solver::Option{Solver}   = nothing
end

glr(m::HuberRegressor) =
    HuberRegression(m.delta, m.lambda, m.gamma;
                    penalty=Symbol(m.penalty),
                    fit_intercept=m.fit_intercept,
                    penalize_intercept=m.penalize_intercept,
                    scale_penalty_with_samples=m.scale_penalty_with_samples)

descr(::Type{HuberRegressor}) = "Robust regression with objective ``∑ρ(Xθ - y) + λ|θ|₂² + γ|θ|₁`` where `ρ` is the Huber Loss."

#= ==================
   QUANTILE REGRESSOR
   ================== =#

"""
Quantile Regression is the same as `RobustRegressor` but with the robust
loss set to `QuantileRho`.

## Parameters

$TYPEDFIELDS

$(example_docstring("QuantileRegressor"))
"""
@with_kw_noshow mutable struct QuantileRegressor <: MMI.Deterministic
    "parametrises the `QuantileRho` function (indicating the  quantile to use
with default `0.5` for the median regression)"
    delta::Real              = 0.5
    "strength of the regulariser if `penalty` is `:l2` or `:l1`.
    Strength of the L2 regulariser if `penalty` is `:en`."
    lambda::Real             = 1.0
    "strength of the L1 regulariser if `penalty` is `:en`."
    gamma::Real              = 0.0
    "the penalty to use, either `:l2`, `:l1`, `:en` (elastic net) or `:none`."
    penalty::SymStr          = :l2
    "whether to fit the intercept or not."
    fit_intercept::Bool      = true
    "whether to penalize the intercept."
    penalize_intercept::Bool = false
    "whether to scale the penalty with the number of samples."
    scale_penalty_with_samples::Bool = true
    "type of solver to use, default if `nothing`."
    solver::Option{Solver}   = nothing
end

glr(m::QuantileRegressor) =
    QuantileRegression(m.delta, m.lambda, m.gamma;
                       penalty=Symbol(m.penalty),
                       fit_intercept=m.fit_intercept,
                       penalize_intercept=m.penalize_intercept,
                       scale_penalty_with_samples=m.scale_penalty_with_samples)

descr(::Type{QuantileRegressor}) = "Robust regression with objective ``∑ρ(Xθ - y) + λ|θ|₂² + γ|θ|₁`` where `ρ` is the Quantile Loss."

#= ==================================
   LEAST ABSOLUTE DEVIATION REGRESSOR
   ================================== =#

"""
Least Absolute Deviation regression with with objective function

``∑ρ(Xθ - y) + n⋅λ|θ|₂² + n⋅γ|θ|₁``

where ``ρ`` is the absolute loss and
``n`` is the number of samples `size(X, 1)`.
With `scale_penalty_with_samples = false` the objective function is
``∑ρ(Xθ - y) + λ|θ|₂² + γ|θ|₁``

See also `RobustRegressor`.

## Parameters

$TYPEDFIELDS

$(example_docstring("LADRegressor"))
"""
@with_kw_noshow mutable struct LADRegressor <: MMI.Deterministic
    "strength of the regulariser if `penalty` is `:l2` or `:l1`.
    Strength of the L2 regulariser if `penalty` is `:en`."
    lambda::Real             = 1.0
    "strength of the L1 regulariser if `penalty` is `:en`."
    gamma::Real              = 0.0
    "the penalty to use, either `:l2`, `:l1`, `:en` (elastic net) or `:none`."
    penalty::SymStr          = :l2
    "whether to fit the intercept or not."
    fit_intercept::Bool      = true
    "whether to penalize the intercept."
    penalize_intercept::Bool = false
    "whether to scale the penalty with the number of samples."
    scale_penalty_with_samples::Bool = true
    "type of solver to use, default if `nothing`."
    solver::Option{Solver}   = nothing
end

glr(m::LADRegressor) =
    LADRegression(m.lambda, m.gamma;
                  penalty=Symbol(m.penalty),
                  fit_intercept=m.fit_intercept,
                  penalize_intercept=m.penalize_intercept,
                  scale_penalty_with_samples=m.scale_penalty_with_samples)

descr(::Type{LADRegressor}) = "Robust regression with objective ``∑ρ(Xθ - y) + λ|θ|₂² + γ|θ|₁`` where `ρ` is the Absolute Loss."


text:
using Pkg
#  If NetABM is not available then it should be installed via
#  Uncomment if not installed
#  using NetABM
#  ] add https://github.com/ollin18/NetABM.jl

function useit(list::Array{Symbol})
    installed = [key for key in keys(Pkg.installed())]
    strpackages = @. string(list)
    uninstalled = setdiff(strpackages,installed)

    map(Pkg.add,uninstalled)
    for package ∈ list
        @eval using $package
    end
end

thep = [:LightGraphs, :StatsBase, :Plots, :GraphPlot, :Statistics, :LaTeXStrings, :Bootstrap, :NetABM]
useit(thep)

global num_agents = 100
global n_boot = 1000
global cil = 0.95
global j
global adoptstep = (0,[0.1])
global maxI = Array{Float64}(undef,0)
global sumI = []

@time let
    for copera in 0.0:0.5:1
        for adap in 0.0:0.5:1
            global meantotI = Array{Float64}(undef,0)
            global Imean = Array{Float64}(undef,0)
            global Ilow = Array{Float64}(undef,0)
            global Ihigh = Array{Float64}(undef,0)
            global Itotall = Array{Float64}(undef,0)
            global Iall = Array{Array}(undef,0)
            global Sall = Array{Array}(undef,0)
            global Rall = Array{Array}(undef,0)
            for a in 1:5#100
                g = barabasi_albert(num_agents,4)
                agents = [Agent(i) for i = 1:num_agents];
                    init_demographics!(agents;states=["S","I"],initial=[1.0,0.0]);
                    pat_0 = rand(1:nv(g))
                    agents[pat_0].state = "I"
                #  init_demographics!(agents;states=["S","I"],initial=[0.9,0.1]);
                set_coop_agents!(agents;p_cop=copera);
                set_adapt_agents!(agents;p_cop=adap);
                map(x -> assign_contacts!(g,x), agents);
                num_S = []
                num_I = []
                num_R = []
                global oldS = ([ag.state for ag in agents if ag.state=="S"] |> length)/num_agents
                global oldI = ([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents
                global totI = ([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents
                for i in 1:100
                    push!(num_S, ([ag.state for ag in agents if ag.state=="S"] |> length)/num_agents)
                    push!(num_I, ([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
                    push!(num_R, ([ag.state for ag in agents if ag.state=="R"] |> length)/num_agents)

                    next_state!(agents;fun=SI_attitude!,inf_prob=0.1, rec_prob=0.03, R=true);
                    update_state!(agents);
                    update_effect_given_distance!(agents,g,adoptstep[1],1,adoptstep[2]);

                    newI = oldS - num_S[end]
                    global totI = totI + newI
                    oldS = num_S[end]

                end
                push!(Iall,num_I)
                push!(Itotall,totI)
                GC.gc()
            end
            for step in eachindex(Iall[1])
                bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
                bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
                bci1 = confint(bs1, BasicConfInt(cil))
                bci2 = confint(bs2, BasicConfInt(cil))
                push!(Imean,bci1[1][1])
                push!(Ilow,bci2[1][2])
                push!(Ihigh,bci2[1][3])
            end
            bs1 = bootstrap(mean, Itotall,BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            push!(sumI,bci1[1][1])
            push!(maxI,maximum(Imean))
        end
        GC.gc()
    end
end
p
savefig(p,"figs/sis_attitudes_01.png")

theme(:default)

plt = heatmap(
        xlabel = "Cooperative Agents",
        ylabel = "Adaptive Agents",
        title = "Maximum Infections",
        guidefontsize = 16,
        tickfontsize = 10,
        xlims = (0,11)
    )
xs = [string(i) for i = 0:0.1:1]
ys = [string(i) for i = 0:0.1:1]
z = reshape(maxI, (11, 11))'
plt = heatmap!(xs, ys, z, aspect_ratio = 1,c = :algae,dpi=300)
plt
png(plt,"figs/heatmap_scaled_maximum.png")
savefig(plt,"figs/heatmap_scaled_maximum.pdf")

plt = heatmap(
        xlabel = "Cooperative Agents",
        ylabel = "Adaptive Agents",
        title = "Total Infections",
        guidefontsize = 16,
        tickfontsize = 10,
        xlims = (0,11)
    )
xs = [string(i) for i = 0:0.1:1]
ys = [string(i) for i = 0:0.1:1]
z = reshape(sumI, (11, 11))'
plt = heatmap!(xs, ys, z, aspect_ratio = 1,c = :algae,dpi=300)
plt
png(plt,"figs/heatmap_scaled_total.png")
savefig(plt,"figs/heatmap_scaled_total.pdf")

################################################################
################################################################
################################################################
################################################################
################################################################
# The same as before but instead of looking at neighbor infections
# we look use a mean cooperation effect by distance

global num_agents = 100
global n_boot = 1000
global cil = 0.95
global j
#  global muladol = [(0,[0.1]), (1,repeat([0.1],2)),(2,repeat([0.1],3)),(3,repeat([0.1],4))]
global muladol = [(i,repeat([0.5],i+1)) for i in 1:10]

global adoptstep = (1,[0.1, 0.1])
global maxI = Array{Float64}(undef,0)
global sumI = []

@time let
    for adoptstep in muladol
    for copera in 0.0:0.5:1
        for adap in 0.0:0.5:1
            global meantotI = Array{Float64}(undef,0)
            global Imean = Array{Float64}(undef,0)
            global Ilow = Array{Float64}(undef,0)
            global Ihigh = Array{Float64}(undef,0)
            global Itotall = Array{Float64}(undef,0)
            global Iall = Array{Array}(undef,0)
            global Sall = Array{Array}(undef,0)
            global Rall = Array{Array}(undef,0)
            for a in 1:5
                g = barabasi_albert(num_agents,9)
                agents = [Agent(i) for i = 1:num_agents];
                init_demographics!(agents;states=["S","I"],initial=[1.0,0.0]);
                pat_0 = rand(1:nv(g))
                agents[pat_0].state = "I"
                #  init_demographics!(agents;states=["S","I"],initial=[0.9,0.1]);
                set_coop_agents!(agents;p_cop=copera);
                set_adapt_agents!(agents;p_cop=adap);
                map(x -> assign_contacts!(g,x), agents);
                num_S = []
                num_I = []
                num_R = []
                global oldS = ([ag.state for ag in agents if ag.state=="S"] |> length)/num_agents
                global oldI = ([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents
                global totI = ([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents
                for i in 1:100
                    push!(num_S, ([ag.state for ag in agents if ag.state=="S"] |> length)/num_agents)
                    push!(num_I, ([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
                    push!(num_R, ([ag.state for ag in agents if ag.state=="R"] |> length)/num_agents)

                    next_state!(agents;fun=SI_attitude!,inf_prob=0.1, rec_prob=0.03, R=true);
                    update_state!(agents);
                    update_effect_given_distance_coop!(agents,g,adoptstep[1],1,adoptstep[2]);

                    newI = oldS - num_S[end]
                    global totI = totI + newI
                    oldS = num_S[end]

                end
                push!(Iall,num_I)
                push!(Itotall,totI)
                GC.gc()
            end
            for step in eachindex(Iall[1])
                bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
                bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
                bci1 = confint(bs1, BasicConfInt(cil))
                bci2 = confint(bs2, BasicConfInt(cil))
                push!(Imean,bci1[1][1])
                push!(Ilow,bci2[1][2])
                push!(Ihigh,bci2[1][3])
            end
            bs1 = bootstrap(mean, Itotall,BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            push!(sumI,bci1[1][1])
            push!(maxI,maximum(Imean))
        end
        GC.gc()
    end
    end
end

theme(:default)
using DelimitedFiles

plt = heatmap(
        xlabel = "Cooperative Agents",
        ylabel = "Adaptive Agents",
        title = "Maximum Infections",
        guidefontsize = 16,
        tickfontsize = 10,
        xlims = (0,11)
    )
xs = [string(i) for i = 0:0.5:1]
ys = [string(i) for i = 0:0.5:1]
z = reshape(maxI, (11, 11, 11))
z = reshape(maxI, (3,3,10))
z1 = z[:,:,4]
plt = heatmap!(xs, ys, z1, aspect_ratio = 1, clim=(0,1), c = :algae,dpi=300)
plt
png(plt,"figs/heatmap_scaled_maximum_coop.png")
savefig(plt,"figs/heatmap_scaled_maximum_coop.pdf")
writedlm("data/max_inf_coop_neigh.csv",maxI)

plt = heatmap(
        xlabel = "Cooperative Agents",
        ylabel = "Adaptive Agents",
        title = "Total Infections",
        guidefontsize = 16,
        tickfontsize = 10,
        xlims = (0,11)
    )
xs = [string(i) for i = 0:0.1:1]
ys = [string(i) for i = 0:0.1:1]
z = reshape(sumI, (11, 11))
plt = heatmap!(xs, ys, z, aspect_ratio = 1, clim=(0,1), c = :algae,dpi=300)
plt
png(plt,"figs/heatmap_scaled_total_coop.png")
savefig(plt,"figs/heatmap_scaled_total_coop.pdf")
writedlm("data/tot_inf_coop_neigh.csv",sumI)


text:
using Test, Statistics, Random, LinearAlgebra, SparseArrays
using DiffEqBiological, DiffEqBase
using OrdinaryDiffEq, StochasticDiffEq, DiffEqJump, SteadyStateDiffEq

# full macro tests
@time begin
  @time @testset "Model Macro" begin include("make_model_test.jl") end
  @time @testset "Gillespie Tests" begin include("gillespie.jl") end
  @time @testset "Test Solvers" begin include("solver_test.jl") end
  @time @testset "Higher Order" begin include("higher_order_reactions.jl") end
  @time @testset "Additional Functions" begin include("func_test.jl") end
  @time @testset "Steady state solver" begin include("steady_state.jl") end
  @time @testset "Mass Action Jumps" begin include("mass_act_jump_tests.jl") end
  @time @testset "Other Tests" begin include("misc_tests.jl") end
  @time @testset "Network query tests" begin include("networkquery_test.jl") end
end

# min macro tests
@time begin
  @time @testset "Model Macro (Min)" begin include("make_model_test_min.jl") end
  @time @testset "Gillespie Tests (Min)" begin include("gillespie_min.jl") end
  @time @testset "Test Solvers (Min)" begin include("solver_test_min.jl") end
  @time @testset "Higher Order (Min)" begin include("higher_order_reactions_min.jl") end
  @time @testset "Additional Functions (Min)" begin include("func_test_min.jl") end
  @time @testset "Steady state solver (Min)" begin include("steady_state_min.jl") end
  @time @testset "Mass Action Jumps (Min)" begin include("mass_act_jump_tests_min.jl") end
end

# tests that handle both macros
@time begin
  @time @testset "Discrete Problem" begin include("discreteproblem_test.jl") end
  @time @testset "Add Reactions API" begin include("addreactions_test.jl") end
end


text:
module SensitiveFloats

export SensitiveFloat,
       zero_sensitive_float,
       clear!,
       multiply_sfs!,
       add_scaled_sfs!,
       combine_sfs!,
       add_sources_sf!,
       set_hess!

abstract ParamSet

"""
A function value and its derivative with respect to its arguments.

Attributes:
  v:  The value
  d:  The derivative with respect to each variable in
      P-dimensional VariationalParams for each of S celestial objects
      in a local_P x local_S matrix.
  h:  The second derivative with respect to each variational parameter,
      in the same format as d.  This is used for the full Hessian
      with respect to all the sources.
"""
type SensitiveFloat{ParamType <: ParamSet, NumType <: Number}
    # Actually a single value, but an Array to avoid memory allocation
    v::Vector{NumType}

    # local_P x local_S matrix of gradients
    d::Matrix{NumType}

    # h is ordered so that p changes fastest.  For example, the indices
    # of a column of h correspond to the indices of d's stacked columns.
    h::Matrix{NumType} # (local_P * local_S) x (local_P * local_S)
    # ids::ParamType
end

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

"""
Set a SensitiveFloat's hessian term, maintaining symmetry.
"""
function set_hess!{ParamType <: ParamSet, NumType <: Number}(
    sf::SensitiveFloat{ParamType, NumType},
    i::Int, j::Int, v::NumType)
  i != j ?
    sf.h[i, j] = sf.h[j, i] = v:
    sf.h[i, j] = v

    true # Set definite return type
end

function zero_sensitive_float{ParamType <: ParamSet}(
  ::Type{ParamType}, NumType::DataType, local_S::Int)
    local_P = length(ParamType)

    v = zeros(NumType, 1)
    d = zeros(NumType, local_P, local_S)
    h = zeros(NumType, local_P * local_S, local_P * local_S)
    SensitiveFloat{ParamType, NumType}(v, d, h)
end

function zero_sensitive_float{ParamType <: ParamSet}(
    ::Type{ParamType}, NumType::DataType)
        # Default to a single source.
        zero_sensitive_float(ParamType, NumType, 1)
end

# If no type is specified, default to using Float64.
function zero_sensitive_float{ParamType <: ParamSet}(
            param_arg::Type{ParamType}, local_S::Int)
    zero_sensitive_float(param_arg, Float64, local_S)
end

function zero_sensitive_float{ParamType <: ParamSet}(
            param_arg::Type{ParamType})
    zero_sensitive_float(param_arg, Float64, 1)
end

function clear!{ParamType <: ParamSet, NumType <: Number}(
            sp::SensitiveFloat{ParamType, NumType})
    clear!(sp, true)
end

function clear!{ParamType <: ParamSet, NumType <: Number}(
                sp::SensitiveFloat{ParamType, NumType}, clear_hessian::Bool)
    fill!(sp.v, zero(NumType))
    fill!(sp.d, zero(NumType))
    if clear_hessian
      fill!(sp.h, zero(NumType))
    end

    true # Set definite return type
end

"""
Factor out the hessian part of combine_sfs!.
"""
function combine_sfs_hessian!{ParamType <: ParamSet,
                              T1 <: Number, T2 <: Number, T3 <: Number}(
        sf1::SensitiveFloat{ParamType, T1},
        sf2::SensitiveFloat{ParamType, T1},
        sf_result::SensitiveFloat{ParamType, T1},
        g_d::Vector{T2}, g_h::Matrix{T3})

    p1, p2 = size(sf_result.h)
    @assert size(sf_result.h) == size(sf1.h) == size(sf2.h)
    @assert p1 == p2 == prod(size(sf1.d)) == prod(size(sf2.d))
    for ind2 = 1:p2
        sf11_factor = g_h[1, 1] * sf1.d[ind2] + g_h[1, 2] * sf2.d[ind2]
        sf21_factor = g_h[1, 2] * sf1.d[ind2] + g_h[2, 2] * sf2.d[ind2]

        @inbounds for ind1 = 1:ind2
            sf_result.h[ind1, ind2] =
                g_d[1] * sf1.h[ind1, ind2] +
                g_d[2] * sf2.h[ind1, ind2] +
                sf11_factor * sf1.d[ind1] +
                sf21_factor * sf2.d[ind1]
            sf_result.h[ind2, ind1] = sf_result.h[ind1, ind2]
        end
    end

    true # Set definite return type
end

"""
Updates sf_result in place with g(sf1, sf2), where
g_d = (g_1, g_2) is the gradient of g and
g_h = (g_11, g_12; g_12, g_22) is the hessian of g,
each evaluated at (sf1, sf2).

The result is stored in sf_result.  The order is done in such a way that
it can overwrite sf1 or sf2 and still be accurate.
"""
function combine_sfs!{ParamType <: ParamSet,
                      T1 <: Number, T2 <: Number, T3 <: Number}(
        sf1::SensitiveFloat{ParamType, T1},
        sf2::SensitiveFloat{ParamType, T1},
        sf_result::SensitiveFloat{ParamType, T1},
        v::T1, g_d::Vector{T2}, g_h::Matrix{T3};
        calculate_hessian::Bool=true)

    # TODO: time consuming **************

    # TODO: this line is allocating a lot of memory and I don't know why.
    # Commenting this line out attributes the same allocation to the next line.
    # Is memory being allocated lazily or misattributed?
    @assert g_h[1, 2] == g_h[2, 1]

    # You have to do this in the right order to not overwrite needed terms.
    if calculate_hessian
        combine_sfs_hessian!(sf1, sf2, sf_result, g_d, g_h)
    end

    for ind in eachindex(sf_result.d)
        sf_result.d[ind] = g_d[1] * sf1.d[ind] + g_d[2] * sf2.d[ind]
    end

    sf_result.v[1] = v

    true # Set definite return type
end

"""
Updates sf1 in place with g(sf1, sf2), where
g_d = (g_1, g_2) is the gradient of g and
g_h = (g_11, g_12; g_12, g_22) is the hessian of g,
each evaluated at (sf1, sf2).

The result is stored in sf1.
"""
function combine_sfs!{ParamType <: ParamSet,
                                            T1 <: Number, T2 <: Number, T3 <: Number}(
        sf1::SensitiveFloat{ParamType, T1},
        sf2::SensitiveFloat{ParamType, T1},
        v::T1, g_d::Vector{T2}, g_h::Matrix{T3};
        calculate_hessian::Bool=true)

    combine_sfs!(sf1, sf2, sf1, v, g_d, g_h, calculate_hessian=calculate_hessian)

    true # Set definite return type
end

# Decalare outside to avoid allocating memory.
const multiply_sfs_hess = Float64[0 1; 1 0]

"""
Updates sf1 in place with sf1 * sf2.
"""
function multiply_sfs!{ParamType <: ParamSet, NumType <: Number}(
        sf1::SensitiveFloat{ParamType, NumType},
        sf2::SensitiveFloat{ParamType, NumType},
        calculate_hessian::Bool=true)

    v = sf1.v[1] * sf2.v[1]
    g_d = NumType[sf2.v[1], sf1.v[1]]
    #const g_h = NumType[0 1; 1 0]

    combine_sfs!(sf1, sf2, v, g_d, multiply_sfs_hess,
                             calculate_hessian=calculate_hessian)

    true # Set definite return type
end

"""
Update sf1 in place with (sf1 + scale * sf2).
"""
function add_scaled_sfs!{ParamType <: ParamSet, NumType <: Number}(
        sf1::SensitiveFloat{ParamType, NumType},
        sf2::SensitiveFloat{ParamType, NumType},
        scale::AbstractFloat, calculate_hessian::Bool)

    sf1.v[1] = sf1.v[1] + scale * sf2.v[1]

    @inbounds for i in eachindex(sf1.d)
        sf1.d[i] = sf1.d[i] + scale * sf2.d[i]
    end

    if calculate_hessian
        p1, p2 = size(sf1.h)
        @assert (p1, p2) == size(sf2.h)
        @inbounds for ind2=1:p2, ind1=1:ind2
            sf1.h[ind1, ind2] += scale * sf2.h[ind1, ind2]
            sf1.h[ind2, ind1] = sf1.h[ind1, ind2]
        end
    end

    true # Set definite return type
end

"""
Adds sf2_s to sf1, where sf1 is sensitive to multiple sources and sf2_s is only
sensitive to source s.
"""
function add_sources_sf!{ParamType <: ParamSet, NumType <: Number}(
        sf_all::SensitiveFloat{ParamType, NumType},
        sf_s::SensitiveFloat{ParamType, NumType},
        s::Int, calculate_hessian::Bool)

    # TODO: This line, too, allocates a lot of memory.    Why?
    sf_all.v[1] = sf_all.v[1] + sf_s.v[1]

    # TODO: time consuming **************
    P = length(ParamType)
    Ph = size(sf_all.h)[1]

    @assert Ph == prod(size(sf_all.d))
    @assert Ph == size(sf_all.h)[2]
    @assert Ph >= P * s
    @assert size(sf_s.d) == (P, 1)
    @assert size(sf_s.h) == (P, P)

    @inbounds for s_ind1 in 1:P
        s_all_ind1 = P * (s - 1) + s_ind1
        sf_all.d[s_all_ind1] = sf_all.d[s_all_ind1] + sf_s.d[s_ind1]
        if calculate_hessian
            @inbounds for s_ind2 in 1:s_ind1
                s_all_ind2 = P * (s - 1) + s_ind2
                sf_all.h[s_all_ind2, s_all_ind1] =
                    sf_all.h[s_all_ind2, s_all_ind1] + sf_s.h[s_ind2, s_ind1]
                # TODO: move outside the loop?
                sf_all.h[s_all_ind1, s_all_ind2] = sf_all.h[s_all_ind2, s_all_ind1]
            end
        end
    end

    true # Set definite return type
end

end


text:
using Bijectors
using Compat: only
using Distributions
using LinearAlgebra
using Optim
using ParameterHandling
using PDMats
using Test
using Zygote

using ParameterHandling: value
using ParameterHandling.TestUtils: test_flatten_interface, test_parameter_interface

const tuple_infers = VERSION < v"1.5" ? false : true

@testset "ParameterHandling.jl" begin
    include("flatten.jl")
    include("parameters.jl")
end


text:
export MvNormalMeanPrecision

import Distributions: logdetcov, distrname, sqmahal, sqmahal!, AbstractMvNormal
import LinearAlgebra: diag, Diagonal, dot
import Base: ndims, precision, length, size, prod

struct MvNormalMeanPrecision{ T <: Real, M <: AbstractVector{T}, P <: AbstractMatrix{T} } <: AbstractMvNormal
    μ :: M
    Λ :: P
end

function MvNormalMeanPrecision(μ::AbstractVector{ <: Real }, Λ::AbstractMatrix{ <: Real })
    T = promote_type(eltype(μ), eltype(Λ))
    return MvNormalMeanPrecision(convert(AbstractArray{T}, μ), convert(AbstractArray{T}, Λ))
end

function MvNormalMeanPrecision(μ::AbstractVector{ <: Integer}, Λ::AbstractMatrix{ <: Integer })
    return MvNormalMeanPrecision(float.(μ), float.(Λ))
end

function MvNormalMeanPrecision(μ::AbstractVector, λ::AbstractVector)
    return MvNormalMeanPrecision(μ, matrix_from_diagonal(promote_type(eltype(μ), eltype(λ)), λ))
end

function MvNormalMeanPrecision(μ::AbstractVector{T}) where T
    return MvNormalMeanPrecision(μ, convert(AbstractArray{T}, ones(length(μ))))
end

Distributions.distrname(::MvNormalMeanPrecision) = "MvNormalMeanPrecision"

weightedmean(dist::MvNormalMeanPrecision) = precision(dist) * mean(dist)

Distributions.mean(dist::MvNormalMeanPrecision)      = dist.μ
Distributions.mode(dist::MvNormalMeanPrecision)      = mean(dist)
Distributions.var(dist::MvNormalMeanPrecision)       = diag(cov(dist))
Distributions.cov(dist::MvNormalMeanPrecision)       = cholinv(dist.Λ)
Distributions.invcov(dist::MvNormalMeanPrecision)    = dist.Λ
Distributions.std(dist::MvNormalMeanPrecision)       = cholsqrt(cov(dist))
Distributions.logdetcov(dist::MvNormalMeanPrecision) = -logdet(invcov(dist))

Distributions.sqmahal(dist::MvNormalMeanPrecision, x::AbstractVector) = sqmahal!(similar(x), dist, x)

function Distributions.sqmahal!(r, dist::MvNormalMeanPrecision, x::AbstractVector)
    μ = mean(dist)
    for i in 1:length(r)
        @inbounds r[i] = μ[i] - x[i]
    end
    return xT_A_x(r, invcov(dist))
end

Base.eltype(::MvNormalMeanPrecision{T})     where T = T
Base.precision(dist::MvNormalMeanPrecision)         = invcov(dist)
Base.length(dist::MvNormalMeanPrecision)            = length(mean(dist))
Base.ndims(dist::MvNormalMeanPrecision)             = length(dist)
Base.size(dist::MvNormalMeanPrecision)              = (length(dist), )

Base.convert(::Type{ <: MvNormalMeanPrecision }, μ::AbstractVector, Λ::AbstractMatrix) = MvNormalMeanPrecision(μ, Λ)

function Base.convert(::Type{ <: MvNormalMeanPrecision{T} }, μ::AbstractVector, Λ::AbstractMatrix) where { T <: Real }
    MvNormalMeanPrecision(convert(AbstractArray{T}, μ), convert(AbstractArray{T}, Λ))
end

vague(::Type{ <: MvNormalMeanPrecision }, dims::Int) = MvNormalMeanPrecision(zeros(dims), fill(tiny, dims))

prod_analytical_rule(::Type{ <: MvNormalMeanPrecision }, ::Type{ <: MvNormalMeanPrecision }) = ProdAnalyticalRuleAvailable()

function Base.prod(::ProdPreserveType, left::MvNormalMeanPrecision, right::MvNormalMeanPrecision)
    Λ = invcov(left) + invcov(right)
    μ = cholinv(Λ) * (invcov(left) * mean(left) + invcov(right) * mean(right))
    return MvNormalMeanPrecision(μ, Λ)
end

function Base.prod(::ProdAnalytical, left::MvNormalMeanPrecision, right::MvNormalMeanPrecision)
    W = precision(left) + precision(right)
    xi = precision(left)*mean(left) + precision(right)*mean(right)
    return MvNormalWeightedMeanPrecision(xi, W)
end


text:
fruits = "🍌 🍎 🍐"

index = firstindex(fruits)
while index <= sizeof(fruits)
    letter = fruits[index]
    println(letter)
    global index = nextind(fruits, index)
end

println(uppercase("abcdefg"))
println('a' in "abcdef")

# This is not valid:
# println("a" in "abcdef")

# These are also valid:
println(in('a', "abcdef"))
println(!in('z', "abcdef"))

# Regex
function usesonly(word, available)
    r = Regex("[^$(available)]")
    !occursin(r, word)
end

println(usesonly("banana", "abn"))
println(usesonly("bananas", "abn"))

m1 = match(r"[^abn]", "banana")
m2 = match(r"[^abn]", "bananas")
println(m1)
println(m2)

println(fieldnames(Regex))
println(fieldnames(RegexMatch))
println(m2.match)
println(m2.offset)


text:
@testset "Simulator: utilities" begin
    # dimensions
    nz = RoboDojo.num_var(biped)
    nθ = RoboDojo.num_data(biped, nf=length(friction_coefficients(biped)))

    @test nz == 57
    @test nθ == 33

    # trajectory
    traj = RoboDojo.Trajectory(biped, 5)
    grad = RoboDojo.GradientTrajectory(biped, 5)
    @test length(traj.q) == 7
    @test length(traj.v) == 6
    @test length(traj.u) == 5
    traj.q[1] .= rand(biped.nq)
    RoboDojo.reset!(traj)
    @test all(traj.q[1] .== 0.0)
    @test length(grad.∂q3∂q1) == 5
    grad.∂q3∂q1[1] .= rand(biped.nq, biped.nq)
    RoboDojo.reset!(grad)
    @test all(grad.∂q3∂q1[1] .== 0.0)

    # disturbances
    dist = RoboDojo.empty_disturbances(biped)
    @test length(dist.w) == biped.nw
    @test length(RoboDojo.disturbances(dist, traj, 1)) == biped.nw

    # policy
    p = RoboDojo.empty_policy(biped)
    @test length(RoboDojo.policy(p, traj, 1)) == biped.nu

    # indices

    idx_z = RoboDojo.indices_z(biped)
    @test idx_z.q == collect(1:9)
    @test idx_z.γ == collect(10:17)
    @test idx_z.sγ == collect(18:25)
    @test idx_z.ψ == collect(26:33)
    @test idx_z.b == collect(34:41)
    @test idx_z.sψ == collect(42:49)
    @test idx_z.sb == collect(50:57)

    idx_θ = RoboDojo.indices_θ(biped, nf=length(friction_coefficients(biped)))
    @test idx_θ.q1 == collect(1:9)
    @test idx_θ.q2 == collect(10:18)
    @test idx_θ.u == collect(19:24)
    @test idx_θ.w == collect(25:24)
    @test idx_θ.f == collect(25:32)
    @test idx_θ.h == collect(33:33)

    idx = RoboDojo.IndicesOptimization()
    @test idx.nz == 0
    @test idx.nΔ == 0

    idx_opt = RoboDojo.indices_optimization(biped)
    @test idx_opt.nz == nz
    @test idx_opt.nΔ == nz
    @test idx_opt.ortz == [[10, 11, 12, 13, 14, 15, 16, 17], [18, 19, 20, 21, 22, 23, 24, 25]]
    @test idx_opt.ortΔ == [[10, 11, 12, 13, 14, 15, 16, 17], [18, 19, 20, 21, 22, 23, 24, 25]]
    @test idx_opt.socz == [[[26, 34], [42, 50]],
                            [[27, 35], [43, 51]],
                            [[28, 36], [44, 52]],
                            [[29, 37], [45, 53]],
                            [[30, 38], [46, 54]],
                            [[31, 39], [47, 55]],
                            [[32, 40], [48, 56]],
                            [[33, 41], [49, 57]]]
    @test idx_opt.socΔ == [[[26, 34], [42, 50]],
                            [[27, 35], [43, 51]],
                            [[28, 36], [44, 52]],
                            [[29, 37], [45, 53]],
                            [[30, 38], [46, 54]],
                            [[31, 39], [47, 55]],
                            [[32, 40], [48, 56]],
                            [[33, 41], [49, 57]]]
    @test idx_opt.equr == collect(1:33)
    @test idx_opt.ortr == collect(34:41)
    @test idx_opt.socr == collect(42:57)
    @test idx_opt.socri == [[42, 43],
                            [44, 45],
                            [46, 47],
                            [48, 49],
                            [50, 51],
                            [52, 53],
                            [54, 55],
                            [56, 57]]
    @test idx_opt.bil == collect(34:57)

    # initialize
    z0 = zeros(nz)
    θ0 = zeros(nθ)
    q0 = rand(biped.nq)
    u0 = rand(biped.nu)
    w0 = rand(biped.nw)
    f0 = friction_coefficients(biped)
    h0 = [0.235]
    RoboDojo.initialize_z!(z0, biped, idx_z, q0)
    @test norm(z0 - [q0; 1.0 * ones(3 * biped.nc); 0.1 * ones(biped.nc); 1.0 * ones(biped.nc); 0.1 * ones(biped.nc)]) < 1.0e-8

    RoboDojo.initialize_θ!(θ0, biped, idx_θ, q0, q0, u0, w0, f0, h0)
    @test norm(θ0 - [q0; q0; u0; w0; f0; h0]) < 1.0e-8

    # floating base dimensions
    @test RoboDojo.floating_base_dim(biped) == 3
end


text:
using GnuplotSimple, Gnuplot
using Test


text:
# ---
# title: 547. Friend Circles
# id: problem547
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Union Find
# link: <https://leetcode.com/problems/friend-circles/description/>
# hidden: true
# ---
#
# There are **N** students in a class. Some of them are friends, while some are
# not. Their friendship is transitive in nature. For example, if A is a
# **direct** friend of B, and B is a **direct** friend of C, then A is an
# **indirect** friend of C. And we defined a friend circle is a group of
# students who are direct or indirect friends.
#
# Given a **N*N** matrix **M** representing the friend relationship between
# students in the class. If M[i][j] = 1, then the ith and jth students are
# **direct** friends with each other, otherwise not. And you have to output the
# total number of friend circles among all the students.
#
# **Example 1:**
#
#
#
#     Input:
#     [[1,1,0],
#      [1,1,0],
#      [0,0,1]]
#     Output: 2
#     Explanation: The 0th and 1st students are direct friends, so they are in a friend circle.
#     The 2nd student himself is in a friend circle. So return 2.
#
#
#
#
# **Example 2:**
#
#
#
#     Input:
#     [[1,1,0],
#      [1,1,1],
#      [0,1,1]]
#     Output: 1
#     Explanation: The 0th and 1st students are direct friends, the 1st and 2nd students are direct friends,
#     so the 0th and 2nd students are indirect friends. All of them are in the same friend circle, so return 1.
#
#
#
#
#
# **Constraints:**
#
#   * `1 <= N <= 200`
#   * `M[i][i] == 1`
#   * `M[i][j] == M[j][i]`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    InflationConstant <: InflationFunction
    InflationConstant(c)

Metodología de inflación constante con valor interanual `c`.
"""
struct InflationConstant <: InflationFunction
    c::Float32
end
InflationConstant() = InflationConstant(1)

# Método para obtener la variación interanual constante igual a c
function (inflfn::InflationConstant)(cs::CountryStructure)
    t = infl_periods(cs::CountryStructure)
    fill(eltype(cs)(inflfn.c), t)
end

# Variación intermensual correspondiente a la variación interanual constante igual a c
function (inflfn::InflationConstant)(base::VarCPIBase{T}) where T
    v = 100 * ((inflfn.c / 100 + 1)^(1 // 12) - 1)
    fill(T(v), periods(base))
end

# Nombre de la medida
measure_name(inflfn::InflationConstant) = "Variación interanual constante igual a " * string(round(inflfn.c, digits=2))


text:
# Day 17: Trick Shot
function load(::Day{17}, path)
    m = match(r"target area: x=(-?\d*)\.\.(-?\d*), y=(-?\d*)\.\.(-?\d*)", readline(path))
    return parse.(Int, m.captures)
end

# Part 1 – highest shot
# For this task, the trajectory in x doesn't matter as the probe will always fall
# vertically into the target area. We just brute-force compute all trajectories
# that end in the target area and return the highest point across all of these trajectories.
function solve1(::Day{17}, data; vmax=150)
    _, _, ymin, ymax = data

    ytrajectories = filter(
        trajectory -> in(last(trajectory), ymin:ymax),
        [cumsum(ẏ₀:-1:ymin) for ẏ₀ in 1:vmax],
    )
    return maximum(maximum.(ytrajectories))
end

# Part 2 – count all viable trajectories
function solve2(::Day{17}, data; vmax=150)
    xmin, xmax, ymin, ymax = data
    xrange = xmin:xmax
    yrange = ymin:ymax
    count = 0

    # List of viable trajectories in x
    xtrajectories = filter(
        trajectory -> any(in(xrange), trajectory), [cumsum(ẋ₀:-1:0) for ẋ₀ in 1:xmax]
    )

    for xtraj in xtrajectories
        # There are two types of xtrajectories:
        # a) those that have zero velocity in x when reaching the target
        # b) those that pass through the target

        steps_until_target = findall(in(xrange), xtraj)
        if length(xtraj) >= 2 && in(xtraj[end], xrange) && xtraj[end] == xtraj[end - 1] # case a
            ytrajectories = filter(
                trajectory -> any(in(yrange), trajectory[first(steps_until_target):end]),
                [cumsum(ẏ₀:-1:ymin) for ẏ₀ in ymin:vmax],
            )
            count += length(ytrajectories)
        else # case b
            vels = Vector{Int}()
            for step in steps_until_target
                ytrajectories = filter(
                    trajectory -> in(last(trajectory), yrange),
                    [cumsum(range(; start=ẏ₀, step=-1, length=step)) for ẏ₀ in ymin:vmax],
                )
                append!(vels, first.(ytrajectories))
            end
            count += length(unique!(vels))
        end
    end
    return count
end

testresult1(::Day{17}) = 45
testresult2(::Day{17}) = 112


text: @test problem2(4e6) == 4613732

max_stars_repo_name: halleysfifthinc/ProjectEuler

max_issues_repo_name: halleysfifthinc/ProjectEuler

max_forks_repo_name: halleysfifthinc/ProjectEuler


text:
# GPUArrays.jl interface

#
# Device functionality
#

## execution

struct CuArrayBackend <: AbstractGPUBackend end

struct CuKernelContext <: AbstractKernelContext end

@inline function GPUArrays.launch_heuristic(::CuArrayBackend, f::F, args::Vararg{Any,N};
                                            maximize_blocksize=false) where {F,N}
    kernel_args = map(cudaconvert, args)
    kernel_tt = Tuple{CuKernelContext, map(Core.Typeof, kernel_args)...}
    kernel = cufunction(f, kernel_tt)
    if maximize_blocksize
        # some kernels benefit (algorithmically) from a large block size
        launch_configuration(kernel.fun)
    else
        # otherwise, using huge blocks often hurts performance: even though it maximizes
        # occupancy, we'd rather have a couple of blocks to switch between.
        launch_configuration(kernel.fun; max_threads=256)
    end
end

@inline function GPUArrays.gpu_call(::CuArrayBackend, f::F, args::TT, threads::Int,
                                    blocks::Int; name::Union{String,Nothing}) where {F,TT}
    @cuda threads=threads blocks=blocks name=name f(CuKernelContext(), args...)
end

## on-device

# indexing

GPUArrays.blockidx(ctx::CuKernelContext) = blockIdx().x
GPUArrays.blockdim(ctx::CuKernelContext) = blockDim().x
GPUArrays.threadidx(ctx::CuKernelContext) = threadIdx().x
GPUArrays.griddim(ctx::CuKernelContext) = gridDim().x

# math

@inline GPUArrays.cos(ctx::CuKernelContext, x) = cos(x)
@inline GPUArrays.sin(ctx::CuKernelContext, x) = sin(x)
@inline GPUArrays.sqrt(ctx::CuKernelContext, x) = sqrt(x)
@inline GPUArrays.log(ctx::CuKernelContext, x) = log(x)

# memory

@inline function GPUArrays.LocalMemory(::CuKernelContext, ::Type{T}, ::Val{dims}, ::Val{id}
                                      ) where {T, dims, id}
    ptr = CUDA._shmem(Val(id), T, Val(prod(dims)))
    CuDeviceArray(dims, reinterpret(LLVMPtr{T, AS.Shared}, ptr))
end

# synchronization

@inline GPUArrays.synchronize_threads(::CuKernelContext) = sync_threads()

#
# Host abstractions
#

GPUArrays.device(A::CuArray) = device(CuCurrentContext())

GPUArrays.backend(::Type{<:CuArray}) = CuArrayBackend()

GPUArrays.unsafe_reinterpret(::Type{T}, A::CuArray, size::NTuple{N, Integer}) where {T, N} =
  CuArray{T,N}(convert(CuPtr{T}, pointer(A)), size, A)


text:
using PurlinLine
using StructuresKit

design_code = "AISI S100-16 ASD"

# Define the properties of each purlin segment along the line.
                #length, section_properties, material_properties
segments = [(23.0*12, 1, 1),
                (2.0*12, 2, 1),
                (2.0*12, 2, 1),
                (23.0*12, 1, 1)]

# Define the purlin spacing.
spacing = 60;  #in.

# Define the roof slope.
roof_slope = 0.0;   #degrees

#Define the purlin cross-section type.
cross_section_dimensions = [("Z", 0.059, 0.91, 2.5, 8.0, 2.5, 0.91, -50.0, 0.0, 90.0, 0.0, -50.0, 3*0.059, 3*0.059, 3*0.059, 3*0.059),
                                        ("Z", 0.059*2, 0.91, 2.5, 8.0, 2.5, 0.91, -50.0, 0.0, 90.0, 0.0, -50.0, 3*0.059, 3*0.059, 3*0.059, 3*0.059)]

material_properties = [(29500.0, 0.30, 55.0, 70.0),
                              (29500.0, 0.30, 55.0, 70.0)]

#type="screw-fastened", thickness, fastener spacing, fastener diameter, fastener_shear_strength

#type="standing seam", thickness, clip spacing, clip stiffness
deck_details = ("screw-fastened", 0.0179, 12.0, 0.212, 2.50)

deck_material_properties = (29500.0, 0.30, 55.0, 70.0)

frame_flange_width = 24.0

support_locations = [0.0, 25.0*12, 50.0*12]

purlin_frame_connections = "bottom flange connection"
# purlin_frame_connections = "anti-roll clip"

bridging_locations =[0.0, 10.0*12, 50.0*12]

#Calculate purlin line design variables from user inputs and store them in the data structure.
purlin_line = PurlinLine.build(design_code, segments, spacing, roof_slope, cross_section_dimensions, material_properties, deck_details, deck_material_properties, frame_flange_width, support_locations, purlin_frame_connections, bridging_locations);
# purlin_line.applied_pressure = 10 / 1000 / 144 #psf to kips/in^2

purlin_line.loading_direction = "gravity"

#Perform a test to collapse.
purlin_line = PurlinLine.test(purlin_line);

purlin_line.applied_pressure*1000*144

using Plots
plot(purlin_line.model.z, purlin_line.internal_forces.T)

max_stars_repo_name: runtosolve/PurlinLine

max_issues_repo_name: runtosolve/PurlinLine

max_forks_repo_name: runtosolve/PurlinLine


text:
@testset "152.maximum-product-subarray.jl" begin
    @test max_product([2,3,-2,4]) == 6
    @test max_product([-2,0,-1]) == 0
    @test max_product([-3, 1, -2, 4, 0, 2, 0, 3, 3, 3, 3, 2, 2, 3, -2, 0, 4, -1, 1, 1, 1, 2, -2, 0, -2, 4, 1, 2, -3, 3]) == 972
end


text:
@testset "7.5 Inverse hyperbolic secant" begin
    include("7.5.1 u (a+b arcsech(c x))^n.jl")
    include("7.5.2 Inverse hyperbolic secant functions.jl")
end

max_stars_repo_path: test/7 Inverse hyperbolic functions/7.5 Inverse hyperbolic secant/.jl


text:
using Test

@testset "Geo.jl" begin
    include("tests.jl")
end


text:
"""
    ESw1(P₀::Vector,P₁::Vector)

Compute the effect size `w` for two sets of `k` (number of cells) probabilities.
`P₀` is the vector of `k` proportions under the null hypothesis
and `P₁` is the vector of `k` proportions under the alternative hypothesis.
"""
function ESw1(P0::Vector,P1::Vector)
    return sqrt(sum((P1-P0).^2./P0))
end

max_stars_repo_name: mwsohn/pwr

max_issues_repo_name: mwsohn/pwr

max_forks_repo_name: mwsohn/pwr


text:
using ZygoteRules
using ZygoteRules: @adjoint

# From the flux file, this one replaces
# @grad -> @adjoint
# .data ->
# ::TrackedArray -> ::AbstractArray
# Delete every line containing track(
# Delete all docstrings -- @adjoint dislikes them
# Delete inv!! whatever that was.

@adjoint function exp_(A::AbstractArray)
    expA = exp_(A)
    expA_copy = copy(expA) # ensures that mutating output won't damage the gradient
    expA, Δ -> ( scale!(expA_copy, Δ) ,) # and we can then mutate the copy here... for 1st deriv?
end
@adjoint function exp!(A::AbstractArray)
    expA = exp!(A)
    expA_copy = copy(expA)
    expA, Δ -> ( scale!(expA_copy, Δ) ,)
end
@adjoint function exp!!(A::AbstractArray)
    expA = exp!(A)
    expA, Δ -> ( scale!(Δ, expA) ,)
end

@adjoint function exp_(name::Symbol, A::AbstractArray)
    expA = exp_(name, A)
    expA_copy = copy_(:exp_copy, expA) # opts in but has a distinct name? NOT safe
    expA, Δ -> (nothing, scale!(expA_copy, Δ) ,)
end

@adjoint log_(A::AbstractArray) =
    log_(A), Δ -> ( iscale_(Δ, A) ,)
@adjoint function log!(A::AbstractArray)
    invA = inv_(A)
    log!(A), Δ -> ( scale!(invA, Δ) ,)
end
@adjoint function log!!(A::AbstractArray)
    invA = inv_(A)
    logA = log!(A)
    logA, Δ -> ( scale!(Δ, invA) ,)
end

@adjoint function log_(name::Symbol, A::AbstractArray)
    invA = inv_(:log_copy, A)
    logA = log_(name, A)
    logA, Δ -> (nothing, scale!(invA, Δ) ,)
end

@adjoint scale_(A::AbstractArray, b::Number) = scale_(A, b), Δ -> (scale_(Δ, b), nothing)
@adjoint scale!(A::AbstractArray, b::Number) = scale!(A, b), Δ -> (scale_(Δ, b), nothing)
@adjoint scale!!(A::AbstractArray, b::Number) = scale!(A, b), Δ -> (scale!(Δ, b), nothing)

# @adjoint scale_(A::AbstractArray, b::TrackedReal) = scale_(A, b), Δ -> (scale_(Δ, b), dot(A,Δ))
# @adjoint scale!(A::AbstractArray, b::TrackedReal) = @error "hmm"

@adjoint scale_(A::AbstractArray, B::Union{Array, RVector}) = scale_(A, B), Δ -> (scale_(Δ, B), nothing)
@adjoint scale!(A::AbstractArray, B::Union{Array, RVector}) = scale!(A, B), Δ -> (scale_(Δ, B), nothing)
@adjoint scale!!(A::AbstractArray, B::Union{Array, RVector}) = scale!(A, B), Δ -> (scale!(Δ, B), nothing)

@adjoint  scale_(A::AbstractArray, B::AbstractArray) =
    scale_(A, B), Δ ->  (scale_(Δ, B), sum!(similar(B), scale_(Δ,A)) )
@adjoint function scale!(A::AbstractArray, B::AbstractArray)
    Ac = copy(A)
    scale!(A, B), Δ -> (scale_(Δ, B), sum!(similar(B), scale!(Ac,Δ)) )
end
@adjoint function scale!!(A::AbstractArray, B::AbstractArray)
    Ac = copy(A)
    scale!(A, B), function(Δ)
        ∇B = sum!(similar(B), scale!(Ac,Δ))
        (scale!(Δ, B),  ∇B)
    end
end

@adjoint  scale_(A::Array, B::AbstractArray) =
    scale_(A, B), Δ ->  (nothing, sum!(similar(B), scale_(Δ,A)) )
@adjoint function scale!(A::Array, B::AbstractArray)
    Ac = copy(A)
    scale!(A, B), Δ -> (nothing, sum!(similar(B), scale!(Ac,Δ)) )
end

@adjoint iscale_(A::AbstractArray, b::Number) = iscale_(A, b), Δ -> (iscale_(Δ, b), nothing)
@adjoint iscale!(A::AbstractArray, b::Number) = iscale!(A, b), Δ -> (iscale!(Δ, b), nothing)

# @adjoint iscale_(A::AbstractArray, b::TrackedReal) = iscale_(A, b), Δ -> (iscale_(Δ, b), nothing)
# @adjoint iscale!(A::AbstractArray, b::TrackedReal) = iscale!(A, b), Δ -> (iscale!(Δ, b), nothing)

@adjoint iscale_(A::AbstractArray, B::Union{Array, RVector}) = iscale_(A, B), Δ -> (iscale_(Δ, B), nothing)
@adjoint iscale!(A::AbstractArray, B::Union{Array, RVector}) = iscale!(A, B), Δ -> (iscale!(Δ, B), nothing)
@adjoint iscale!!(A::AbstractArray, B::Union{Array, RVector}) = scale!(A, inv!(B)), Δ -> (scale!(Δ, B), nothing)

@adjoint iscale_(A::AbstractArray, B::AbstractArray) = iscale_(A, B), Δ -> (iscale_(Δ, B), sum!(similar(B), -1 .* Δ .* A ./ B .^2))
# @adjoint function iscale!(A::AbstractArray, B::AbstractArray)
#     Ac = copy(A)
#     iscale!(A, B), Δ -> (scale(Δ, B), sum!(similar(B), scale!(Ac,Δ)) )
# end

@adjoint iscale_(A::Array, B::AbstractArray) = iscale_(A, B), Δ -> (nothing, sum!(similar(B), -1 .* Δ .* A ./ B .^2))

# @adjoint inv!(A::AbstractArray, b::Number=1) = inv_(A, b), Δ -> (-1 .* Δ .* b .* A .^ (-2) , nothing) # one copy
@adjoint function inv_(A::AbstractArray, b=1)
    invA = inv_(A, b)
    invA_copy = copy(invA) # don't invert twice?
    invA, Δ -> (invA_copy .= -1 .* Δ .* b .* invA_copy .^ 2 , nothing) # total one copy
end
# @adjoint inv!(A::AbstractArray, b::Number=1) = inv!(A, b), Δ -> (-1 .* Δ .* b .* A .^ 2 , nothing) # one copy
@adjoint function inv!(A::AbstractArray, b=1)
    invA = inv!(A, b)
    invA_copy = copy(invA) # keep gradient calc safe from downstream mutations
    invA, Δ -> (invA_copy .= -1 .* Δ .* b .* invA_copy .^ 2 , nothing) # total one copy
end


text:
# some kernels
epanechnikov(x::Tp) where {Tp<:AbstractFloat} = Tp( (abs(x) ≤  1.0) ? 0.75 * (1.0-x*x) : 0.0 )
gaussian(x::Tp) where {Tp<:AbstractFloat} = Tp( exp( -0.5 * x * x) * 0.3989422804014327 )
triweight(x::Tp) where {Tp<:AbstractFloat} = Tp( (abs(x) ≤ 1.0) ? 1.09375 * (1.0-x*x)^3  : 0.0 )
uniform(x::Tp) where {Tp<:AbstractFloat} = Tp( (abs(x) ≤ 1.0) ? 0.5 : 0.0 )
cosinus(x::Tp) where {Tp<:AbstractFloat} = Tp( (abs(x) ≤ 1.0) ?   0.7853981633974483    * cos( 1.5707963267948966 * x ) : 0.0 )
quartic(x::Tp) where {Tp<:AbstractFloat} = Tp( (abs(x) ≤ 1.0) ? 0.9375 * (1.0 - x*x)^2 : 0.0 )


text:
function isMandel(rep::Int, lim::Int64, x::Int64, y::Int64)
  zr::Int64 = 0
  zi::Int64 = 0
  for i in 0:(rep-1)
    zrNext = div(zr*zr,lim) - div(zi*zi,lim) + x
    ziNext = div(zr*zi*2,lim) + y
    dist2 = div(zrNext*zrNext,lim) + div(ziNext*ziNext,lim)
    if 4*lim < dist2
      return false
    end
    zr=zrNext
    zi=ziNext
  end
  return true
end

function countMandel( rep, div )
  lim::Int64 = 2^div
  count::Int64 = 0
  for y::Int64 in -lim:lim
    for x::Int64 in -lim:lim
      if (isMandel( rep, lim, x, y))
        count+=1
      end
    end
  end
  return count
end

function main()
  rep = length(ARGS)<1 ? 100 : tryparse(Int, ARGS[1])
  div = length(ARGS)<2 ? 8 : tryparse(Int, ARGS[2])
  println( "rep=", rep, "  div=", div )
  println( "result=", countMandel(rep,div))
end

main()


text:
struct Butterfly{T} <: Factorization{T}
    columns::Vector{Matrix{T}}
    factors::Vector{Vector{IDPackedV{T}}}
    permutations::Vector{Vector{ColumnPermutation}}
    indices::Vector{Vector{Int}}
    temp1::ThreadSafeVector{T}
    temp2::ThreadSafeVector{T}
    temp3::ThreadSafeVector{T}
    temp4::ThreadSafeVector{T}
end

function size(B::Butterfly, dim::Integer)
    if dim == 1
        m = 0
        for i = 1:length(B.columns)
            m += size(B.columns[i], 1)
        end
        return m
    elseif dim == 2
        n = 0
        for j = 1:length(B.factors[1])
            n += size(B.factors[1][j], 2)
        end
        return n
    else
        return 1
    end
end

size(B::Butterfly) = size(B, 1), size(B, 2)

function Butterfly(A::AbstractMatrix{T}, L::Int; isorthogonal::Bool = false, opts...) where T
    m, n = size(A)
    tL = 2^L

    LRAOpts = LRAOptions(T; opts...)
    LRAOpts.rtol = eps(real(T))*max(m, n)

    columns = Vector{Matrix{T}}(undef, tL)
    factors = Vector{Vector{IDPackedV{T}}}(undef, L+1)
    permutations = Vector{Vector{ColumnPermutation}}(undef, L+1)
    indices = Vector{Vector{Int}}(undef, L+1)
    cs = Vector{Vector{Vector{Int}}}(undef, L+1)

    factors[1] = Vector{IDPackedV{T}}(undef, tL)
    permutations[1] = Vector{ColumnPermutation}(undef, tL)
    indices[1] = Vector{Int}(undef, tL+1)
    cs[1] = Vector{Vector{Int}}(undef, tL)

    ninds = range(1, stop=n+1, length=tL+1)
    indices[1][1] = 1
    for j = 1:tL
        nl = round(Int, ninds[j])
        nu = round(Int, ninds[j+1]) - 1
        nd = nu-nl+1
        if isorthogonal
            factors[1][j] = IDPackedV{T}(collect(1:nd),Int[],Array{T}(undef,nd,0))
        else
            factors[1][j] = idfact!(A[:,nl:nu], LRAOpts)
        end
        permutations[1][j] = factors[1][j][:P]
        indices[1][j+1] = indices[1][j] + size(factors[1][j], 1)
        cs[1][j] = factors[1][j].sk .+ nl .- 1
    end

    ii, jj = 2, (tL>>1)
    for l = 2:L+1
        factors[l] = Vector{IDPackedV{T}}(undef, tL)
        permutations[l] = Vector{ColumnPermutation}(undef, tL)
        indices[l] = Vector{Int}(undef, tL+1)
        cs[l] = Vector{Vector{Int}}(undef, tL)

        ctr = 0
        minds = range(1, stop=m+1, length=ii+1)
        indices[l][1] = 1
        for i = 1:ii
            shft = 2jj*div(ctr,2jj)
            ml = round(Int, minds[i])
            mu = round(Int, minds[i+1]) - 1
            for j = 1:jj
                cols = vcat(cs[l-1][2j-1+shft],cs[l-1][2j+shft])
                lc = length(cols)
                Av = A[ml:mu,cols]
                if maximum(abs, Av) < floatmin(real(T))/eps(real(T))
                    factors[l][j+ctr] = IDPackedV{T}(Int[], collect(1:lc), Array{T}(undef,0,lc))
                else
                    LRAOpts.rtol = eps(real(T))*max(mu-ml+1, lc)
                    factors[l][j+ctr] = idfact!(Av, LRAOpts)
                end
                permutations[l][j+ctr] = factors[l][j+ctr][:P]
                indices[l][j+ctr+1] = indices[l][j+ctr] + size(factors[l][j+ctr], 1)
                cs[l][j+ctr] = cols[factors[l][j+ctr].sk]
            end
            ctr += jj
        end
        ii <<= 1
        jj >>= 1
    end

    minds = range(1, stop=m+1, length=tL+1)
    for i = 1:tL
        ml = round(Int, minds[i])
        mu = round(Int, minds[i+1]) - 1
        columns[i] = A[ml:mu, cs[L+1][i]]
    end

    kk = sumkmax(indices)

    Butterfly(columns, factors, permutations, indices, threadsafezeros(T, kk), threadsafezeros(T, kk), threadsafezeros(T, kk), threadsafezeros(T, kk))
end

if VERSION ≥ v"0.7-"
    LinearAlgebra.adjoint(B::Butterfly) = Adjoint(B)
    LinearAlgebra.transpose(B::Butterfly) = Transpose(B)
end

function sumkmax(indices::Vector{Vector{Int}})
    ret = 0
    @inbounds for j = 1:length(indices)
        ret = max(ret, indices[j][end])
    end
    ret
end

#### Helper

function rowperm!(fwd::Bool, x::AbstractVecOrMat, p::Vector{Int}, jstart::Int)
    n = length(p)
    jshift = jstart-1
    rmul!(p, -1)
    @inbounds if (fwd)
        for i = 1:n
            p[i] > 0 && continue
            j    =    i
            p[j] = -p[j]
            k    =  p[j]
            while p[k] < 0
                x[jshift+j], x[jshift+k] = x[jshift+k], x[jshift+j]
                j    =    k
                p[j] = -p[j]
                k    =  p[j]
            end
        end
    else
        for i = 1:n
            p[i] > 0 && continue
            p[i] = -p[i]
            j    =  p[i]
            while p[j] < 0
                x[jshift+i], x[jshift+j] = x[jshift+j], x[jshift+i]
                p[j] = -p[j]
                j    =  p[j]
            end
        end
    end
    x
end

function rowperm!(fwd::Bool, y::AbstractVector, x::AbstractVector, p::Vector{Int}, jstart::Int)
    n = length(p)
    jshift = jstart-1
    @inbounds if (fwd)
        @simd for i = 1:n
            y[jshift+i] = x[jshift+p[i]]
        end
    else
        @simd for i = 1:n
            y[jshift+p[i]] = x[jshift+i]
        end
    end
    y
end

## ColumnPermutation
mul!(A::ColPerm, B::AbstractVecOrMat, jstart::Int) = rowperm!(false, B, A.p, jstart)
At_mul_B!(A::ColPerm, B::AbstractVecOrMat, jstart::Int) = rowperm!(true, B, A.p, jstart)
Ac_mul_B!(A::ColPerm, B::AbstractVecOrMat, jstart::Int) = At_mul_B!(A, B, jstart)

mul!(y::AbstractVector, A::ColPerm, x::AbstractVector, jstart::Int) = rowperm!(false, y, x, A.p, jstart)
At_mul_B!(y::AbstractVector, A::ColPerm, x::AbstractVector, jstart::Int) = rowperm!(true, y, x, A.p, jstart)
Ac_mul_B!(y::AbstractVector, A::ColPerm, x::AbstractVector, jstart::Int) = At_mul_B!(y, x, A, jstart)

# Fast mul!, At_mul_B!, and Ac_mul_B! for an ID. These overwrite the output.

function mul!(y::AbstractVecOrMat{T}, A::IDPackedV{T}, P::ColumnPermutation, x::AbstractVecOrMat{T}, istart::Int, jstart::Int) where {T}
    k, n = size(A)
    At_mul_B!(P, x, jstart)
    copyto!(y, istart, x, jstart, k)
    mul!(y, A.T, x, istart, jstart+k)
    mul!(P, x, jstart)
    y
end

function mul!(y::AbstractVector{T}, A::IDPackedV{T}, P::ColumnPermutation, x::AbstractVector{T}, temp::AbstractVector{T}, istart::Int, jstart::Int) where {T}
    k, n = size(A)
    At_mul_B!(temp, P, x, jstart)
    copyto!(y, istart, temp, jstart, k)
    mul!(y, A.T, temp, istart, jstart+k)
    y
end

### mul!, At_mul_B!, and  Ac_mul_B! for a Butterfly factorization.
mul!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = mul_col_J!(u, B, b, 1)

for f! in (:At_mul_B!, :Ac_mul_B!)
    @eval begin
        function $f!(y::AbstractVecOrMat{T}, A::IDPackedV{T}, P::ColumnPermutation, x::AbstractVecOrMat{T}, istart::Int, jstart::Int) where {T}
            k, n = size(A)
            copyto!(y, istart, x, jstart, k)
            $f!(y, A.T, x, istart+k, jstart)
            mul!(P, y, istart)
            y
        end

        function $f!(y::AbstractVector{T}, A::IDPackedV{T}, P::ColumnPermutation, x::AbstractVector{T}, temp::AbstractVector{T}, istart::Int, jstart::Int) where {T}
            k, n = size(A)
            copyto!(temp, istart, x, jstart, k)
            $f!(temp, A.T, x, istart+k, jstart)
            mul!(y, P, temp, istart)
            y
        end
    end
end

if VERSION < v"0.7-"
    Base.A_mul_B!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = mul_col_J!(u, B, b, 1)
    Base.At_mul_B!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = At_mul_B_col_J!(u, B, b, 1)
    Base.Ac_mul_B!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = Ac_mul_B_col_J!(u, B, b, 1)
else
    LinearAlgebra.mul!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = mul_col_J!(u, B, b, 1)
    LinearAlgebra.mul!(u::Vector{T}, Bt::Transpose{T,Butterfly{T}}, b::Vector{T}) where T = At_mul_B_col_J!(u, parent(Bt), b, 1)
    LinearAlgebra.mul!(u::Vector{T}, Bc::Adjoint{T,Butterfly{T}}, b::Vector{T}) where T = Ac_mul_B_col_J!(u, parent(Bc), b, 1)
end

function mul_col_J!(u::VecOrMat{T}, B::Butterfly{T}, b::VecOrMat{T}, J::Int) where T
    L = length(B.factors) - 1
    tL = 2^L

    M = size(b, 1)

    COLSHIFT = M*(J-1)

    temp1 = B.temp1
    temp2 = B.temp2
    temp3 = B.temp3
    fill!(temp1, zero(T))
    fill!(temp2, zero(T))

    factors = B.factors[1]
    permutations = B.permutations[1]
    inds = B.indices[1]
    nu = 0
    for j = 1:tL
        nl = nu+1
        nu += size(factors[j], 2)
        mul!(temp1, factors[j], permutations[j], b, inds[j], nl+COLSHIFT)
    end

    ii, jj = 2, (tL>>1)
    for l = 2:L+1
        factors = B.factors[l]
        permutations = B.permutations[l]
        indsout = B.indices[l]
        indsin = B.indices[l-1]
        ctr = 0
        for i = 1:ii
            shft = 2jj*div(ctr,2jj)
            for j = 1:jj
                mul!(temp2, factors[j+ctr], permutations[j+ctr], temp1, temp3, indsout[j+ctr], indsin[2j+shft-1])
            end
            ctr += jj
        end
        temp1, temp2 = temp2, fill!(temp1, zero(T))
        ii <<= 1
        jj >>= 1
    end

    columns = B.columns
    inds = B.indices[L+1]
    mu = 0
    for i = 1:tL
        ml = mu+1
        mu += size(columns[i], 1)
        mul!(u, columns[i], temp1, ml+COLSHIFT, inds[i])
    end

    u
end

for f! in (:At_mul_B!,:Ac_mul_B!)
    f_col_J! = Meta.parse(string(f!)[1:end-1]*"_col_J!")
    @eval begin
        function $f_col_J!(u::VecOrMat{T}, B::Butterfly{T}, b::VecOrMat{T}, J::Int) where T
            L = length(B.factors) - 1
            tL = 2^L

            M = size(b, 1)

            COLSHIFT = M*(J-1)

            temp1 = B.temp1
            temp2 = B.temp2
            temp3 = B.temp3
            temp4 = B.temp4
            fill!(temp1, zero(T))
            fill!(temp2, zero(T))
            fill!(temp3, zero(T))

            columns = B.columns
            inds = B.indices[L+1]
            mu = 0
            for i = 1:tL
                ml = mu+1
                mu += size(columns[i], 1)
                $f!(temp1, columns[i], b, inds[i], ml+COLSHIFT)
            end

            ii, jj = tL, 1
            for l = L+1:-1:2
                factors = B.factors[l]
                permutations = B.permutations[l]
                indsout = B.indices[l-1]
                indsin = B.indices[l]
                ctr = 0
                for i = 1:ii
                    shft = 2jj*div(ctr,2jj)
                    fill!(temp4, zero(T))
                    for j = 1:jj
                        $f!(temp3, factors[j+ctr], permutations[j+ctr], temp1, temp4, indsout[2j+shft-1], indsin[j+ctr])
                        addtemp3totemp2!(temp2, temp3, indsout[2j+shft-1], indsout[2j+shft+1]-1)
                    end
                    ctr += jj
                end
                temp1, temp2 = temp2, fill!(temp1, zero(T))
                ii >>= 1
                jj <<= 1
            end

            factors = B.factors[1]
            permutations = B.permutations[1]
            inds = B.indices[1]
            nu = 0
            for j = 1:tL
                nl = nu+1
                nu += size(factors[j], 2)
                $f!(u, factors[j], permutations[j], temp1, nl+COLSHIFT, inds[j])
            end

            u
        end
    end
end

function addtemp3totemp2!(temp2::AbstractVector, temp3::AbstractVector, i1::Int, i2::Int)
    z = zero(eltype(temp3))
    @inbounds @simd for i = i1:i2
        temp2[i] += temp3[i]
        temp3[i] = z
    end
    temp2
end

function allranks(B::Butterfly)
    L = length(B.factors)-1
    tL = 2^L
    ret = zeros(Int, tL, L+1)
    @inbounds for l = 1:L+1
        for j = 1:tL
            ret[j,l] = size(B.factors[l][j], 1)
        end
    end

    ret
end