text:
function tvdrk3!(system::LegSystem,times::CVTimer,n::Int64,splits::Vector{Int64},terms::Vector{Int64},uprox::Float64)
    # allocation for invariant/boundary value RK3 loop
    dW = zeros(2);
    W10 = zeros(2,length(system.branches.ID));
    W20 = zeros(2,length(system.branches.ID));
    W11 = zeros(2,length(system.branches.ID));
    W21 = zeros(2,length(system.branches.ID));
    W12 = zeros(2,length(system.branches.ID));
    W22 = zeros(2,length(system.branches.ID));
    lf = zeros(length(system.branches.ID));
    lb = zeros(length(system.branches.ID));
    ret1 = zeros(2);
    ret2 = zeros(2);
    ret3 = zeros(1);
    ret4 = zeros(1);

    # allocation for interior value RK3 loop
    dU = zeros((length(system.solverparams.colsint))*2);
    Q1 = zeros(length(system.solverparams.colsint),length(system.branches.ID));
    A1 = zeros(length(system.solverparams.colsint),length(system.branches.ID));
    Q2 = zeros(length(system.solverparams.colsint),length(system.branches.ID));
    A2 = zeros(length(system.solverparams.colsint),length(system.branches.ID));
    Acat = zeros(length(system.solverparams.colsint)+2);
    Qcat = zeros(length(system.solverparams.colsint)+2);

    # allocation for WENO3 spatial derivative subroutine
    ISA = zeros(3);
    ISQ = zeros(3);
    PA = Vector{Function}(length(system.solverparams.colsint));
    PQ = Vector{Function}(length(system.solverparams.colsint));
    Aint = Vector{Float64}(length(system.solverparams.colsint));
    Qint = Vector{Float64}(length(system.solverparams.colsint));
    wA = Vector{Float64}(3);
    wQ = Vector{Float64}(3);
    c = Vector{Float64}(3);
    a = Vector{Float64}(2);
    p = Vector{Float64}(2);
    x = Vector{Float64}(2);
    Al = Vector{Float64}(2);
    Ar = Vector{Float64}(2);
    Ql = Vector{Float64}(2);
    Qr = Vector{Float64}(2);
    Fl = Vector{Float64}(2);
    Fr = Vector{Float64}(2);
    Hl = Vector{Float64}(2);
    Hr = Vector{Float64}(2);
    Pi = Vector{Float64}(3);
    σ = Vector{Float64}(2);
    eigs = Vector{Float64}(2);
    abso = Vector{Float64}(2);
    J = zeros(2,2);
    Pret = zeros(1);
    s = zeros(1);

    # first RK step
    # println("First RK:")
    tic();
    for i = 1:length(system.branches.ID)
        # println("First step, artery $i, time step $n.")
        # scalings
        # println("Scalings:")
        As = ((system.solverparams.rho.*system.branches.c0[i][end].^2)./system.branches.beta[i][end]).^2;
        Qs = As.*system.branches.c0[i][end];
        xs = system.branches.lengthinmm[i].*LegHemodynamics.mmTom .- system.branches.k[i];
        ts = As.*xs./Qs;
        Acat .= system.branches.A[i][:,n+1]/As;
        Qcat .= system.branches.Q[i][:,n+1]/Qs;
        # println("Initial area: $(system.branches.A[i][:,n+1]./As)")
        # println("Initial flowrate: $(system.branches.Q[i][:,n+1]./Qs)")
        # invariant step
        # println("Invariant calculation:")
        invariants!(ret1,ret2,ret3,ret4,Acat,Qcat,
            system.branches.c0[i][end],system.branches.beta[i][end],system.solverparams.rho,
            system.solverparams.mu,system.solverparams.diffusioncoeff,Qs,As)
        W10[:,i] .= ret1;
        W20[:,i] .= ret2;
        lf[i] = ret3[1];
        lb[i] = ret4[1];
        # println("Forward invariant at time step $n, artery $i: $(W1n[i]) m/s")
        # println("Invariant update:")
        invariantodes!(dW,W10[:,i],W20[:,i],lf[i],lb[i],Qcat,Acat,
            system.solverparams.rho,system.solverparams.mu,system.solverparams.diffusioncoeff,
            system.branches.k[i]./xs,ts,xs,Qs,As);
        # println("Invariant assignment:")
        if i == 1 # aortic root
            system.branches.W1[i] = (W10[end,i] .+ system.solverparams.h./ts.*dW[1]).*Qs./As;
            system.branches.W2root = (W20[1,i] .+ system.solverparams.h./ts.*dW[2]).*Qs./As;
        elseif isempty(system.branches.children[i]) # terminals
            system.branches.W1end[i] = (W10[end,i] .+ system.solverparams.h./ts.*dW[1]).*Qs./As;
            system.branches.W2[i] = (W20[1,i] .+ system.solverparams.h./ts.*dW[2]).*Qs./As;
        else # all others
            system.branches.W1[i] = (W10[end,i] .+ system.solverparams.h./ts.*dW[1]).*Qs./As;
            system.branches.W2[i] = (W20[1,i] .+ system.solverparams.h./ts.*dW[2]).*Qs./As;
        end
        # primitive step
        # println("WENO3:")
        weno3!(dU,system.branches.k[i]./xs,Acat,
            Qcat,system.branches.beta[i][end],
            system.solverparams.rho,system.solverparams.mu,
            system.solverparams.diffusioncoeff,ts,xs,Qs,As,ISA,ISQ,PA,PQ,wA,wQ,
            c,a,p,x,Al,Ar,Ql,Qr,Hl,Hr,Fl,Fr,Pi,σ,eigs,abso,J,Aint,Qint,Pret,s)
        # println("Intermediate step:")
        A1[:,i] .= (system.branches.A[i][2:end-1,n+1]./As .+
            system.solverparams.h./ts.*dU[system.solverparams.acols[1:end-2]]).*As;
        Q1[:,i] .= (system.branches.Q[i][2:end-1,n+1]./Qs .+
            system.solverparams.h./ts.*dU[system.solverparams.qcols[1:end.-2].-2]).*Qs;
    end
    times.tfd += toq();

    # first set of boundary values
    # println("Split solve:")
    tic();
    for i = 1:length(splits)
        # println("Solving interior boundary, parent ID: $(splits[i])")
        children = system.branches.children[splits[i]];
        Q = zeros(length(children)+1);
        A = zeros(length(children)+1);
        Q[1] = system.branches.Q[splits[i]][system.solverparams.JL,n+1];
        A[1] = system.branches.A[splits[i]][system.solverparams.JL,n+1];
        for j = 1:length(children)
            Q[j+1] = system.branches.Q[children[j]][1,n+1];
            A[j+1] = system.branches.A[children[j]][1,n+1];
        end
        LegHemodynamics.interiorbcs!(system,n,splits[i],Q,A,children);
    end
    times.ts += toq();
    # println("0D-1D coupling:")
    tic();
    LegHemodynamics.coupling!(system,n,system.solverparams.h,terms,uprox);
    times.tc += toq();

    # second RK step
    # println("Second RK:")
    tic()
    for i = 1:length(system.branches.ID)
        # println("Second step, artery $i, time step $n.")
        Acat[1] = system.branches.A[i][1,n+2];
        Acat[2:end-1] .= A1[:,i];
        Acat[end] = system.branches.A[i][end,n+2];
        Qcat[1] = system.branches.Q[i][1,n+2];
        Qcat[2:end-1] .= Q1[:,i];
        Qcat[end] = system.branches.Q[i][end,n+2];
        # scalings
        As = ((system.solverparams.rho.*system.branches.c0[i][end].^2)./system.branches.beta[i][end]).^2;
        Qs = As.*system.branches.c0[i][end];
        xs = system.branches.lengthinmm[i].*LegHemodynamics.mmTom .- system.branches.k[i];
        ts = As.*xs./Qs;
        Acat /= As;
        Qcat /= Qs;
        # println("Area after first step: $(Acat./As)")
        # println("Flowrate after first step: $(Qcat./Qs)")
        # invariant step
        invariants!(ret1,ret2,ret3,ret4,Acat,Qcat,
            system.branches.c0[i][end],system.branches.beta[i][end],system.solverparams.rho,
            system.solverparams.mu,system.solverparams.diffusioncoeff,Qs,As)
        W11[:,i] .= ret1;
        W21[:,i] .= ret2;
        lf[i] = ret3[1];
        lb[i] = ret4[1];
        invariantodes!(dW,W11[:,i],W21[:,i],lf[i],lb[i],Qcat,Acat,
            system.solverparams.rho,system.solverparams.mu,system.solverparams.diffusioncoeff,
            system.branches.k[i]./xs,ts,xs,Qs,As);
        if i == 1 # aortic root
            system.branches.W1[i] = (0.75.*W10[end,i] .+ 0.25.*W11[end,i] .+
                0.25.*system.solverparams.h./ts.*dW[1]).*Qs./As;
            system.branches.W2root = (0.75.*W20[1,i] .+ 0.25.*W21[1,i] .+
                0.25.*system.solverparams.h./ts.*dW[2]).*Qs./As;
        elseif isempty(system.branches.children[i]) # terminals
            system.branches.W1end[i] = (0.75.*W10[end,i] .+ 0.25.*W11[end,i] .+
                0.25.*system.solverparams.h./ts.*dW[1]).*Qs./As;
            system.branches.W2[i] = (0.75.*W20[1,i] .+ 0.25.*W21[1,i] .+
                0.25.*system.solverparams.h./ts.*dW[2]).*Qs./As;
        else # all others
            system.branches.W1[i] = (0.75.*W10[end,i] .+ 0.25.*W11[end,i] .+
                0.25.*system.solverparams.h./ts.*dW[1]).*Qs./As;
            system.branches.W2[i] = (0.75.*W20[1,i] .+ 0.25.*W21[1,i] .+
                0.25.*system.solverparams.h./ts.*dW[2]).*Qs./As;
        end
        # primitive step
        weno3!(dU,system.branches.k[i]./xs,Acat,
            Qcat,system.branches.beta[i][end],
            system.solverparams.rho,system.solverparams.mu,
            system.solverparams.diffusioncoeff,ts,xs,Qs,As,ISA,ISQ,PA,PQ,wA,wQ,
            c,a,p,x,Al,Ar,Ql,Qr,Hl,Hr,Fl,Fr,Pi,σ,eigs,abso,J,Aint,Qint,Pret,s)
        A2[:,i] .= (0.75.*system.branches.A[i][2:end-1,n+1]./As .+ 0.25.*A1[:,i]./As .+
            0.25.*system.solverparams.h./ts.*dU[system.solverparams.acols[1:end-2]]).*As;
        Q2[:,i] .= (0.75.*system.branches.Q[i][2:end-1,n+1]./Qs .+ 0.25.*Q1[:,i]./Qs .+
            0.25.*system.solverparams.h./ts.*dU[system.solverparams.qcols[1:end.-2].-2]).*Qs;
    end
    times.tfd += toq();

    # second set of boundary values
    tic();
    for i = 1:length(splits)
        children = system.branches.children[splits[i]];
        Q = zeros(length(children)+1);
        A = zeros(length(children)+1);
        Q[1] = system.branches.Q[splits[i]][system.solverparams.JL,n+2];
        A[1] = system.branches.A[splits[i]][system.solverparams.JL,n+2];
        for j = 1:length(children)
            Q[j+1] = system.branches.Q[children[j]][1,n+2];
            A[j+1] = system.branches.A[children[j]][1,n+2];
        end
        LegHemodynamics.interiorbcs!(system,n,splits[i],Q,A,children);
    end
    times.ts += toq();
    tic()
    LegHemodynamics.coupling!(system,n,0.5*system.solverparams.h,terms,uprox);
    times.tc += toq();

    # third RK step
    # println("Third RK:")
    tic();
    for i = 1:length(system.branches.ID)
        # println("Third step, artery $i, time step $n.")
        Acat[1] = system.branches.A[i][1,n+2];
        Acat[2:end-1] .= A2[:,i];
        Acat[end] = system.branches.A[i][end,n+2];
        Qcat[1] = system.branches.Q[i][1,n+2];
        Qcat[2:end-1] .= Q2[:,i];
        Qcat[end] = system.branches.Q[i][end,n+2];
        # scalings
        As = ((system.solverparams.rho.*system.branches.c0[i][end].^2)./system.branches.beta[i][end]).^2;
        Qs = As.*system.branches.c0[i][end];
        xs = system.branches.lengthinmm[i].*LegHemodynamics.mmTom .- system.branches.k[i];
        ts = As.*xs./Qs;
        Acat /= As;
        Qcat /= Qs;
        # println("Area after second step: $(Acat./As)")
        # println("Flowrate after second step: $(Qcat./Qs)")
        # invariant step
        invariants!(ret1,ret2,ret3,ret4,Acat,Qcat,
            system.branches.c0[i][end],system.branches.beta[i][end],system.solverparams.rho,
            system.solverparams.mu,system.solverparams.diffusioncoeff,Qs,As)
        W12[:,i] .= ret1;
        W22[:,i] .= ret2;
        lf[i] = ret3[1];
        lb[i] = ret4[1];
        invariantodes!(dW,W12[:,i],W22[:,i],lf[i],lb[i],Qcat,Acat,
            system.solverparams.rho,system.solverparams.mu,system.solverparams.diffusioncoeff,
            system.branches.k[i]./xs,ts,xs,Qs,As);
        if i == 1 # aortic root
            system.branches.W1[i] .= (1/3.*W10[end,i] .+ 2/3.*W12[end,i] .+
                2/3.*system.solverparams.h./ts.*dW[1]).*Qs./As;
            system.branches.W2root = (1/3.*W20[1,i] .+ 2/3.*W22[1,i] .+
                2/3.*system.solverparams.h./ts.*dW[2]).*Qs./As;
        elseif isempty(system.branches.children[i]) # terminals
            system.branches.W1end[i] .= (1/3.*W10[end,i] .+ 2/3.*W12[end,i] .+
                2/3.*system.solverparams.h./ts.*dW[1]).*Qs./As;
            system.branches.W2[i] .= (1/3.*W20[1,i] .+ 2/3.*W22[1,i] .+
                2/3.*system.solverparams.h./ts.*dW[2]).*Qs./As;
        else # all others
            system.branches.W1[i] .= (1/3.*W10[end,i] .+ 2/3.*W12[end,i] .+
                2/3.*system.solverparams.h./ts.*dW[1]).*Qs./As;
            system.branches.W2[i] .= (1/3.*W20[1,i] .+ 2/3.*W22[1,i] .+
                2/3.*system.solverparams.h./ts.*dW[2]).*Qs./As;
        end
        # primitive step
        weno3!(dU,system.branches.k[i]./xs,Acat,
            Qcat,system.branches.beta[i][end],
            system.solverparams.rho,system.solverparams.mu,
            system.solverparams.diffusioncoeff,ts,xs,Qs,As,ISA,ISQ,PA,PQ,wA,wQ,
            c,a,p,x,Al,Ar,Ql,Qr,Hl,Hr,Fl,Fr,Pi,σ,eigs,abso,J,Aint,Qint,Pret,s)
        system.branches.A[i][2:end-1,n+2] .= (1/3.*system.branches.A[i][2:end-1,n+1]./As .+ 2/3.*A2[:,i]./As .+
            2/3.*system.solverparams.h./ts.*dU[system.solverparams.acols[1:end-2]]).*As;
        system.branches.Q[i][2:end-1,n+2] .= (1/3.*system.branches.Q[i][2:end-1,n+1]./Qs .+ 2/3.*Q2[:,i]./Qs .+
            2/3.*system.solverparams.h./ts.*dU[system.solverparams.qcols[1:end.-2].-2]).*Qs;
    end
    times.tfd += toq();

    # third set of boundary values
    # println("Third interior boundaries:")
    tic();
    for i = 1:length(splits)
        # println("Interior boundary allocation:")
        children = system.branches.children[splits[i]];
        Q = zeros(length(children)+1);
        A = zeros(length(children)+1);
        Q[1] = system.branches.Q[splits[i]][system.solverparams.JL,n+2];
        A[1] = system.branches.A[splits[i]][system.solverparams.JL,n+2];
        for j = 1:length(children)
            Q[j+1] = system.branches.Q[children[j]][1,n+2];
            A[j+1] = system.branches.A[children[j]][1,n+2];
        end
        # println("Interior boundary solution:")
        LegHemodynamics.interiorbcs!(system,n,splits[i],Q,A,children);
    end
    times.ts += toq();
    # println("Third coupling:")
    tic();
    LegHemodynamics.coupling!(system,n,system.solverparams.h,terms,uprox);
    times.tc += toq();

end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
#!/usr/bin/env julia

using FileIO
using JLD2
using ArgParse
using CumulantsUpdates
using SymmetricTensors
using Combinatorics
using Distributed
import CumulantsUpdates: cumulants, moment

"""
  cumspeedups(m::Int, n::Vector{Int}, t::Int, tup::Vector{Int}, b::Int)

Returns Vector, a computional speedup of m'th cumulant update of n[i] variate data
"""
function cumspeedups(m::Int, n::Vector{Int}, t::Int, u::Vector{Int}, b::Int, bc::Int = 2)
  compt = zeros(length(u), length(n))
  updt = copy(compt)
  precomp(m)
  for i in 1:length(n)
    println("cumulants calc n = ", n[i])
    X = randn(t, n[i])
    t1 = Float64(time_ns())
    cumulants(X, m, bc)
    compt[:,i] .= Float64(time_ns())-t1
    dm = DataMoments(X, m, b)
    for j in 1:length(u)
      println("update u = ", u[j])
      Xup = rand(u[j], n[i])
      t1 = Float64(time_ns())
      _ = cumulantsupdate!(dm, Xup)
      updt[j,i] = Float64(time_ns()) - t1
    end
  end
  compt, updt
end

"""
  precomp(m::Int)

precompiles updates functions
"""
function precomp(m::Int)
  X = randn(15, 10)
  cumulants(X[1:10,:], m, 2)
  dm = DataMoments(X[1:10,:], m, 4)
  cumulantsupdate!(dm, X[10:15,:])
end

"""
  savecomptime(m::Int, n::Vector{Int}, t::Int, tup::Vector{Int}, b::Int, p::Int)

Saves comptime parameters into a .jld2 file
"""
function savecomptime(m::Int, n::Vector{Int}, t::Int, tup::Vector{Int}, b::Int, p::Int)
  filename = replace("res/$(m)_$(t)_$(n)_$(tup)_$(p).jld2", "["=>"")
  filename = replace(filename, "]"=>"")
  filename = replace(filename, " "=>"")
  compt = Dict{String, Any}()
  cumtime = cumspeedups(m, n, t, tup, b)
  push!(compt, "cumulants" => cumtime[1])
  push!(compt, "cumulants updat" => cumtime[2])
  push!(compt, "tm" => t./(2*tup.+bellnum(m)))
  push!(compt, "t" => t)
  push!(compt, "n" => n)
  push!(compt, "m" => m)
  push!(compt, "tup" => tup)
  push!(compt, "x" => "tup")
  push!(compt, "functions" => [["cumulants", "cumulants updat"]])
  save(filename, compt)
end

"""
  main(args)

Returns plots of the speedup of updates. Takes optional arguments from bash
"""
function main(args)
  s = ArgParseSettings("description")
  @add_arg_table s begin
      "--order", "-d"
        help = "d, the order of cumulant, ndims of cumulant's tensor"
        default = 4
        arg_type = Int
        "--blocksize", "-b"
        help = "the size of blocks of the block structure"
        default = 4
        arg_type = Int
      "--nvar", "-n"
        nargs = '*'
        default = [60]
        help = "n, numbers of marginal variables"
        arg_type = Int
      "--dats", "-t"
        help = "t, numbers of data records"
        default = 500000
        arg_type = Int
      "--updates", "-u"
        help = "u, size of the update"
        nargs = '*'
        default = [12000, 18000, 24000, 30000, 36000]
        arg_type = Int
      "--nprocs", "-p"
        help = "number of processes"
        default = 6
        arg_type = Int
    end
  parsed_args = parse_args(s)
  m = parsed_args["order"]
  n = parsed_args["nvar"]
  t = parsed_args["dats"]
  tup = parsed_args["updates"]
  b = parsed_args["blocksize"]
  p = parsed_args["nprocs"]
  if p > 1
    addprocs(p)
    eval(Expr(:toplevel, :(@everywhere using CumulantsUpdates)))
  end
  println("number of workers = ", nworkers())
  savecomptime(m, n, t, tup, b, p)
end

main(ARGS)

ext: jl


text:
using PointPatterns
using Meshes
using GeoStatsBase
using Test, Random

# list of tests
testfiles = [
  "processes.jl"
]

@testset "PointPatterns.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
"""
    DiscreteMarkovChain{T1,T2}

Discrete-time Markov chain with finite state space.

# Fields
- `π0::T1`: initial state distribution.
- `P::T2`: state transition matrix.
"""
struct DiscreteMarkovChain{T1<:AbstractVector{<:Real},T2<:AbstractMatrix{<:Real}} <:
       AbstractMarkovChain
    π0::T1
    P::T2
end

MeasureTheory.sampletype(::DiscreteMarkovChain) = Vector{Int}

## Access

nb_states(mc::DiscreteMarkovChain) = length(mc.π0)

initial_distribution(mc::DiscreteMarkovChain) = mc.π0

transition_matrix(mc::DiscreteMarkovChain) = mc.P

## Simulation

function Base.rand(rng::AbstractRNG, mc::DiscreteMarkovChain, T::Integer)
    states = Vector{Int}(undef, T)
    states[1] = rand(rng, Dists.Categorical(mc.π0))
    transitions = [Dists.Categorical(mc.P[s, :]) for s = 1:nb_states(mc)]  # TODO: @view
    for t = 2:T
        states[t] = rand(rng, transitions[states[t-1]])
    end
    return states
end

Base.rand(mc::DiscreteMarkovChain, T::Integer) = rand(GLOBAL_RNG, mc, T)

## Logpdf

function MeasureTheory.logdensity(mc::DiscreteMarkovChain, states::Vector{Integer})
    T = length(states)
    l = log(mc.π0[states[1]])
    for t = 2:T
        l += log(mc.P[states[t-1], states[t]])
    end
    return l
end

## Prior

"""
    DiscreteMarkovChainPrior

Define a Dirichlet prior on the initial distribution and on the transitions from each state.

# Fields
- `π0α`: Dirichlet parameter for the initial distribution
- `Pα`: Dirichlet parameters for the transition matrix
"""
struct DiscreteMarkovChainPrior <: AbstractMarkovChainPrior
    π0α::Vector{Float64}
    Pα::Matrix{Float64}
end

## Prior logpdf

function MeasureTheory.logdensity(prior::DiscreteMarkovChainPrior, mc::DiscreteMarkovChain)
    l = logdensity(Dists.Dirichlet(prior.π0α), Dists.Categorical(mc.π0))
    for s = 1:nb_states(mc)
        l += logdensity(Dists.Dirichlet(prior.Pα[s, :]), Dists.Categorical(mc.P[s, :]))
    end
    return l
end

## Fitting

"""
    DiscreteMarkovChainStats

Sufficient statistics for the likelihood of a `DiscreteMarkovChain`.
"""
struct DiscreteMarkovChainStats <: Dists.SufficientStats
    initialization::Vector{Float64}
    transition_count::Matrix{Float64}
end

function Dists.suffstats(::Type{<:DiscreteMarkovChain}, states::Vector{<:Integer})
    S = maximum(states)
    initialization = collect(1:S) .== states[1]
    transition_count = zeros(Int, S, S)
    for t = 1:length(states)-1
        transition_count[states[t], states[t+1]] += 1
    end
    return DiscreteMarkovChainStats(initialization, transition_count)
end

function Dists.suffstats(
    ::Type{<:DiscreteMarkovChain},
    γ::AbstractMatrix,
    ξ::AbstractArray{<:Real,3},
)
    T, S, _ = size(ξ)
    initialization = γ[1, :]
    transition_count = zeros(Float64, S, S)
    for i = 1:S, j = 1:S
        transition_count[i, j] = sum(@view ξ[:, i, j])
    end
    return DiscreteMarkovChainStats(initialization, transition_count)
end

function Dists.suffstats(
    ::Type{<:DiscreteMarkovChain},
    prior::DiscreteMarkovChainPrior,
    args...,
)
    ss = suffstats(DiscreteMarkovChain, args...)
    ss.initialization .+= (prior.π0α .- 1)
    ss.transition_count .+= (prior.Pα .- 1)
    return ss
end

function Dists.fit_mle(::Type{<:DiscreteMarkovChain}, ss::DiscreteMarkovChainStats)
    π0 = ss.initialization
    π0 ./= sum(π0)
    P = ss.transition_count
    P ./= sum(ss.transition_count, dims = 2)
    return DiscreteMarkovChain(π0, P)
end

## Asymptotics

function stationary_distribution(mc::DiscreteMarkovChain)
    π = real.(eigvecs(transition_matrix(mc)')[:, end])
    return π / sum(π)
end


text:
using Stochy

# This is useful for testing inference over models where the number of
# ERP varies across executions.

@pp function f(n)
    factor(n > 3 ? -1 : 0)
    if n < 0
        list()
    else
        flip(0.7) ? 0 .. f(n-1) : list()
    end
end

@pp function g()
    f(5)
end

dist = @pp enum() do
    local x = g()
    x
end

println(dist)
println("\n\n")

dist = @pp mh(1000) do
    local x = g()
    x
end

println(dist)
println("\n\n")


text:
# test script to be run on Travis
# NOTE remove when ContinuousTransformations.jl is registered
Pkg.clone("https://github.com/tpapp/ContinuousTransformations.jl.git")
Pkg.clone(pwd())
Pkg.build("IndirectLikelihood")
Pkg.test("IndirectLikelihood"; coverage=true)


text:
using NCDatasets

abstract type AbstractNcData <: AbstractData end

gettypeparams(data::AbstractData, c::Component) = nothing
axis(data::AbstractData, c::Component) = axis(data, naming(c))
issimple(data::AbstractData, name::Symbol) = false

function components end

struct NcData <: AbstractNcData
    dataset::Dataset
end

function NcData(filename::String)
    ds = Dataset(filename)
    (haskey(ds.attrib, "network_pypsa_version") ? PypsaNcData : NcData)(ds)
end

function load(filename::String)::AbstractData
    if endswith(filename, ".nc")
        Data(fallback=NcData(filename))
    else
        error("file ending of '$filename' not recognized")
    end
end

function Base.get(data::NcData, args...)
    da = data.dataset[String(naming(args...))]
    AxisArray(da[:], (axis(data, n) for n = dimnames(da))...)
end

axis(data::AbstractNcData, n::Symbol) = Axis{n}(nomissing(data.dataset[String(n)][:]))

# Interface for the abstract type AbstractData
#
# * Required methods
#
# Base.get(data, class, param)
# returns an AxisArray of `param` from `class`
#
# devices(data)
# returns Array{DataType} with the described devices
#
# classes(data, ::Type{T}) where T<:Device
# returns Array{Symbol} of the described classes

struct ComponentDesc
    name::Symbol
    componenttype
    data::Dict{Symbol,Union{Scalar,AxisArray}}
end
ComponentDesc(tup) = ComponentDesc(tup...)

struct Data{T} <: AbstractData
    components::Dict{Symbol,ComponentDesc}
    axes::Dict{Symbol,Axis}
    fallback::T

    function Data(components::Dict{Symbol,ComponentDesc}, fallback::Union{Nothing,AbstractData})
        data = new{typeof(fallback)}(components, Dict{Symbol,Axis}(), fallback)
        for cd in values(data.components)
            pushaxes!(data, cd)
        end
        data
    end
end
Data(; components=Dict{Symbol,ComponentDesc}(), fallback=nothing) = Data(components, fallback)
Data(components::Vector{ComponentDesc}; fallback=nothing) = Data(Dict(cd.name => cd for cd in components), fallback)
Data(components::Vararg{Tuple{Symbol,Any,Dict{Symbol}}}; kwargs...) = Data(ComponentDesc.(components); kwargs...)

function pushaxes!(data::AbstractData, ax::Axis; component="unknown", attribute="unknown")
    ## Duck-typing on that data has an `axes` dictionary since we want to use this function for PypsaNcData as well
    axname = axisname(ax)
    if haskey(data.axes, axname)
        if data.axes[axname] != ax
            error("Incompatible axes: The $attribute attribute of component $component has the axis $axname, which we previously encountered with different values:\n\t$(data.axes[axname].val) != $(ax.val) (previous != new)")
        end
    else
        data.axes[axname] = ax
    end
end

function pushaxes!(data::Data, cd::ComponentDesc)
    for (name, attr) in cd.data
        if isa(attr, AxisArray)
            for ax in AxisArrays.axes(attr)
                pushaxes!(data, ax, component=cd.name, attribute=name)
            end
        else
            pushaxes!(data, Axis{name}([name]), component=cd.name, attribute=name)
        end
    end
end

Base.push!(data::Data, cd, ax::Axis{name}) where name = push!(data, cd, Axis{name}(Symbol.(ax.val)))
function Base.push!(data::Data, cd::ComponentDesc, ax::Axis{name,Array{Symbol,1}}) where name
    if name != cd.name
        @warn "Renaming principal axis $(AxisArrays.axisname(ax)) to $(cd.name)"
        ax = Axis{cd.name}(ax.val)
    end
    data.components[cd.name] = cd
    data.axes[cd.name] = ax # We overwrite, since we overwrite the component as well
    pushaxes!(data, cd)
end

_as_axisarrays(d::Dict{Symbol}{T}) where T <: Union{Scalar,AxisArray} = d
Base.push!(data::Data, ::Type{T}, ax::Axis{name}; parameters...) where {name,T <: Component} =
    push!(data, ComponentDesc(name, T, Dict(parameters)), ax)

_components(data) = map(cd->(cd.name => cd.componenttype), values(data.components))
components(data::Data{Nothing}) = _components(data)
components(data::Data{<:AbstractData}) =
    unique(p->p.first, vcat(components(data.fallback), _components(data)))

function axis(data::Data{<:AbstractData}, n::Symbol)
    ret = get(data.axes, n, nothing)
    !isnothing(ret) ? ret : axis(data.fallback, n)
end
axis(data::Data{Nothing}, n::Symbol) = data.axes[n]

function Base.get(data::Data{Nothing}, c::Component, quantity)
    ret = get(data.components[naming(c)].data, quantity, nothing)
    !isnothing(ret) ? ret : getdefault(c, quantity)
end
function Base.get(data::Data{<:AbstractData}, c::Component, quantity)
    componentdesc = get(data.components, naming(c), nothing)
    if !isnothing(componentdesc)
        ret = get(componentdesc.data, quantity, nothing)
        if isnothing(ret)
            ret = get(data.fallback, c, quantity)
        end
    else
        ret = get(data.fallback, c, quantity)
    end

    !isnothing(ret) ? ret : getdefault(c, quantity)
end

store_results!(::AbstractData, model) = error("Not implemented")
store_results!(data::Data{<:AbstractData}, model) = store_results!(data.fallback, model)

issimple(data::Data, name::Symbol) = (ax = axis(data, name); length(ax) == 1 && AxisArrays.axisname(ax) == name)

gettypeparams(data::Data{Nothing}, device::Device) = nothing
gettypeparams(data::Data{<:AbstractData}, device::Device) = gettypeparams(data.fallback, device)

include("pypsa.jl")
include("powersystems.jl")
include("pypsanetwork.jl")

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
"""
`throttle(dt, input::AbstractObservable)`

Throttle a signal to update at most once every `dt` seconds. The throttled signal holds
the last update of the `input` signal during each `dt` second time window.
"""
function throttle(dt, obs::AbstractObservable{T}) where {T}
    throttled = Observable{T}(obs[])
    updatable = Observable(true)
    pending = Ref(false)
    update!() = (throttled[] = obs[]; pending[] = false; updatable[] = false)
    on(updatable) do val
        val ? pending[] && update!() : Timer(t -> updatable[] = true, dt)
    end
    on(obs) do val
        updatable[] ? update!() : (pending[] = true)
    end
    throttled
end


text:
# Functions for negative binomial outcome model in GBM.
#
# Copyright (c) 2020: Jeffrey W. Miller.
# This file is released under the MIT "Expat" License.

module Outcome_NB

export compute_MuWE!, compute_MuWE, loglikelihood, compute_logdispersion_derivatives, update_logdispersions

using Statistics: mean
using SpecialFunctions: logabsgamma, digamma, trigamma
lgamma_(x) = logabsgamma(x)[1]

# Compute Mu, W, and E
function compute_MuWE!(Mu,W,E,Y,logMu,r)
    Mu .= exp.(logMu)
    W .= (r.*Mu)./(r.+Mu)
    E .= (Y - Mu).*(W./Mu)
end

function compute_MuWE(Y,logMu,r)
    I,J = size(Y)
    Mu = zeros(I,J)
    W = zeros(I,J)
    E = zeros(I,J)
    compute_MuWE!(Mu,W,E,Y,logMu,r)
    return Mu,W,E
end

# log-likelihood of the Negative-Binomial distribution with mean Mu and inverse-dispersion r, on data Y.
# loglikelihood(Y,Mu,r) = sum(lgamma_.(r.+Y) .- lgamma_.(Y.+1) .- lgamma_.(r) .+ r.*log.(r) .+ Y.*log.(Mu) .- (r.+Y).*log.(r.+Mu))
# loglikelihood(Y,Mu,r) = sum(lgamma_.(r.+Y) .- lgamma_.(r) .- lgamma_.(Y.+1) .- Y.*log1p.(r./Mu) .- r.*log1p.(Mu./r))
function loglikelihood(Y,Mu,r)
    I,J = size(Y)
    ll = 0.0
    for i = 1:I
        for j = 1:J
            ll += lgamma_(r[i,j]+Y[i,j]) - lgamma_(r[i,j]) - lgamma_(Y[i,j]+1) - Y[i,j]*log1p(r[i,j]/Mu[i,j]) - r[i,j]*log1p(Mu[i,j]/r[i,j])
        end
    end
    return ll
end

digamma_delta(y,r) = ((r > 10^8) ? log1p(y/r) : digamma(y+r)-digamma(r))
trigamma_delta(y,r) = ((r > 10^8) ? -(y/r)/(y+r) : trigamma(y+r)-trigamma(r))

function compute_logdispersion_derivatives(Y,Mu,r)
    I,J = size(Y)
    D1 = zeros(I,J)
    D2 = zeros(I,J)
    for i = 1:I, j = 1:J
        D1[i,j] = -r[i,j]*(digamma_delta(Y[i,j],r[i,j]) - log1p(Mu[i,j]/r[i,j]) - (Y[i,j]-Mu[i,j])/(r[i,j]+Mu[i,j]))

        if isinf(Mu[i,j]^2)
            D2[i,j] = -D1[i,j] + (r[i,j]^2)*trigamma_delta(Y[i,j],r[i,j]) + r[i,j]
        else
            D2[i,j] = -D1[i,j] + (r[i,j]^2)*trigamma_delta(Y[i,j],r[i,j]) + (Y[i,j] + Mu[i,j]^2/r[i,j]) / (1.0 + Mu[i,j]/r[i,j])^2
        end
    end
    return D1,D2
end

function update_logdispersions(S,T,omega,Y,Mu,r,p,max_step,b_S,b_T,update_S,update_T,verbose)
    I,J = size(Y)

    if update_S
    # Update S
    D1,D2 = compute_logdispersion_derivatives(Y,Mu,r)
    g = vec(sum(D1; dims=2)) .- p.lambda_S*(S .- p.mu_S)
    h = vec(sum(D2; dims=2)) .- p.lambda_S
    delta = (g./h).*(h.<0) .+ (-g).*(h.>=0)
    S = S - delta.*min.(b_S./abs.(delta), 1)
    if any(abs.(delta) .> b_S) && verbose; println(" max step size enforced in S[i] update for one or more i."); end
    b_S = [(abs(delta[i]) > b_S[i] ?  b_S[i]/2 : max_step) for i = 1:I]
    shift = log(mean(exp.(S)))
    S = S .- shift  # Project onto constrained space
    omega = omega + shift  # Compensate to preserve likelihood
    r = exp.(-(S.+T'.+omega))
    end

    if update_T
    # Update T
    D1,D2 = compute_logdispersion_derivatives(Y,Mu,r)
    g = vec(sum(D1; dims=1)) .- p.lambda_T*(T .- p.mu_T)
    h = vec(sum(D2; dims=1)) .- p.lambda_T
    delta = (g./h).*(h.<0) .+ (-g).*(h.>=0)
    T = T - delta.*min.(b_T./abs.(delta), 1)
    if any(abs.(delta) .> b_T) && verbose; println(" max step size enforced in T[j] update for one or more j."); end
    b_T = [(abs(delta[j]) > b_T[j] ?  b_T[j]/2 : max_step) for j = 1:J]
    shift = log(mean(exp.(T)))
    T = T .- shift  # Project onto constrained space
    omega = omega + shift  # Compensate to preserve likelihood
    r = exp.(-(S.+T'.+omega))
    end

    return S,T,omega,r,b_S,b_T
end

end # module Outcome_NB


text:
using jlmie
using Test

@testset "jlmie" begin
    @test jlmie_Qsca(1,1) == [0.0]
    @test abs.(jlmie_Qsca(4,1) - [6.062]) < [0.01]
end

max_stars_repo_name: Hinamoooon/jlmie

max_issues_repo_name: Hinamoooon/jlmie

max_forks_repo_name: Hinamoooon/jlmie


text:
const spectroscopic = "spdfghiklmnoqrtuv"
spectroscopic_label(ℓ) =
    ℓ + 1 ≤ length(spectroscopic) ? spectroscopic[ℓ+1] : "[$(ℓ)]"

# Nicer string representation for rationals
rs(r::Number) = "$(r)"
rs(r::Rational) = denominator(r) == 1 ? "$(numerator(r))" : "$(numerator(r))/$(denominator(r))"


text:
using Revise
using AdFem
using PyPlot
using LinearAlgebra
using Statistics
using MAT

mmesh = Mesh(50, 50, 1/50, degree=2)

left = bcnode((x,y)->x<1e-5, mmesh)
right = bcedge((x1,y1,x2,y2)->(x1>0.049-1e-5) && (x2>0.049-1e-5), mmesh)

t1 = eval_f_on_boundary_edge((x,y)->1.0e-4, right, mmesh)
t2 = eval_f_on_boundary_edge((x,y)->0.0, right, mmesh)
rhs = compute_fem_traction_term(t1, t2, right, mmesh)

E = 1.5
ν = 0.3
D = compute_plane_stress_matrix(E*ones(get_ngauss(mmesh)), ν*ones(get_ngauss(mmesh)))
K = compute_fem_stiffness_matrix(D, mmesh)

bdval = [eval_f_on_boundary_node((x,y)->0.0, left, mmesh);
        eval_f_on_boundary_node((x,y)->0.0, left, mmesh)]
DOF = [left;left .+ mmesh.ndof]
K, rhs = impose_Dirichlet_boundary_conditions(K, rhs, DOF, bdval)
u = K\rhs
sess = Session(); init(sess)
S = run(sess, u)

matwrite("fenics/data1.mat", Dict("u"=>S))

close("all")
figure(figsize=(20, 5))
subplot(131)
visualize_scalar_on_fem_points(S[1:mmesh.nnode], mmesh)
title("x displacement")
subplot(132)
visualize_scalar_on_fem_points(S[1+mmesh.ndof:mmesh.nnode + mmesh.ndof], mmesh)
title("y displacement")
subplot(133)
Dval = run(sess, D)
visualize_von_mises_stress(Dval, S, mmesh)
title("von Mises Stress")
savefig("fenics/fwd1.png")


text:
"""
    module IdentityMatrix

Methods for `Diagonal`, `FillArrays.Fill`, and `FillArrays.Eye`.
"""
module IdentityMatrix

using LinearAlgebra, FillArrays

using  FillArrays: AbstractFill, ZerosVecOrMat, getindex_value

using Base: promote_op, has_offset_axes

import Base: any, all, inv, permutedims, imag, one, isone, zero, oneunit,
       sqrt, sum, prod, first, last, minimum, maximum, extrema,
       kron, Matrix

import Base: *, /, +, -, ^, ==, \

import LinearAlgebra: diag, triu, triu!, tril, tril!, eigmin, eigmax,
       norm, normp, norm1, norm2, normInf, normMinusInf, opnorm, isposdef,
       det, logdet

# FillArrays
export Eye, Fill, Ones, Zeros

# LinearAlgebra
export Diagonal

include("diagonal.jl")

# FIXME: replace checkuniquedim with something more efficient (from LinearAlgebra or Base)
checkuniquedim(D::LinearAlgebra.Diagonal) = size(D, 1)
checkuniquedim(AV::AbstractVector) = size(AV, 1)

function checkuniquedim(A::AbstractMatrix)
    m, n = size(A)
    m == n || throw(DimensionMismatch("matrix is not square: dimensions are $(size(A))"))
    return m
end

function checkuniquedim(A, B)
    m1 = checkuniquedim(A)
    m2 = checkuniquedim(B)
    m1 == m2 || throw(DimensionMismatch("matrix A has dimensions $(size(A)). matirix B has dimensions $(size(B))"))
    return m1
end

# Following are functions that I determined will *not* benefit from
# a specialized method.
# The fallback method for (to AbstractArray) is already efficient for Eye for:
# length, firstindex, lastindex, axes, checkbounds

# The fallback method (to Diagonal) is already efficient for Eye for:
# getindex, factorize, real, float, ishermitian, issymmetric, isdiag, istriu,
# istril, log, most trig functions

# The default method for Diagonal returns Diagonal with a dense diagonal

# Test these well!
one(DF::Diagonal{T,V}) where {T, V <: AbstractFill} = Eye{T}(size(DF, 1))
one(AF::AbstractFill{T, 2, <: Any}) where T = Eye{T}(size(AF, 1))

#one(IM::Eye) = IM

imag(IM::Eye{T}) where T = Diagonal(Zeros{real(T)}(size(IM, 1)))
oneunit(IM::Eye) = one(IM)
zero(IM::Eye{T}) where T = Diagonal(Zeros{T}(size(IM, 1)))

isposdef(::Eye) = true

# Return a Vector to agree with other `diag` methods
diag(IM::Eye{T}) where T = ones(T, size(IM, 1))

## Reduction

sum(f::Function, IM::Eye{T}) where T = (m = size(IM, 1); (m * (m - 1)) * f(zero(T)) + m * f(one(T)))
sum(IM::Eye{T}) where T = convert(T, size(IM, 1))

function sum(f::Function, x::Fill)
    dims = size(x)
    return prod(dims) * f(FillArrays.getindex_value(x))
end
sum(x::Fill) = sum(identity, x)

prod(f, IM::Eye{T}) where T = (m = size(IM, 1); f(zero(T))^(m * (m - 1)) * f(one(T))^m)
prod(IM::Eye{T}) where T = size(IM, 1) > 1 ? zero(T) : one(T)

function prod(f::Function, x::Fill)
    dims = size(x)
    return f(FillArrays.getindex_value(x))^prod(dims)
end
prod(x::Fill) = prod(identity, x)

norm2(IM::Eye{T}) where T = sqrt(T(size(IM, 1)))
norm1(IM::Eye{T}) where T = T(size(IM, 1))
normInf(IM::Eye{T}) where T = one(T)

function normMinusInf(IM::Eye{T}) where T
    m = size(IM, 1)
    m == 1 && return one(T)
    return zero(T)
end

normp(IM::Eye{T}, p::Real) where T = (m = size(IM, 1); return iszero(p) ? float(T(size(IM, 1))) : T(m)^(1 / p))

# We have to reproduce this, because LinearAlgebra does not split norm0 into a function.
# So, we are unable to replace it.
# We could file an issue.
function norm(IM::Eye{T}, p::Real=2) where T
    isempty(IM) && return float(zero(T))
    if p == 2
        return norm2(IM)
    elseif p == 1
        return norm1(IM)
    elseif p == Inf
        return normInf(IM)
    elseif p == 0
        return float(T(size(IM, 1)))
    elseif p == -Inf
        return normMinusInf(IM)
    else
        return normp(IM, p)
    end
end

LinearAlgebra.opnorm(IM::Eye{T}, p::Real=1) where T = (isempty(IM) ? zero(T) : one(T)) |> float

for f in (:first, :last)
    @eval begin
        function (Base.$f)(IM::Eye{T}) where T
            size(IM, 1) == 0 && throw(BoundsError("0-element $(typeof(IM))", 0))
            return  one(T)
        end
    end
end

function minimum(IM::Eye{T}) where T
    m = size(IM, 1)
    m > 1 && return zero(T)
    m < 1 && return minimum(T[]) # Error
    return one(T)
end

function maximum(IM::Eye{T}) where T
    size(IM, 1) == 0 && return maximum(T[])
    return one(T)
end

extrema(IM::Eye) = (minimum(IM), maximum(IM)) # FIXME: implement extrema(IM, dims = dims)

# StatsBase.mean is sum/length and therefore is efficient.
# StatsBase.mean(IM::Eye{T}) where T = (m = size(IM, 1); convert(T, m * (m - 1)))
function median end
function IdentityMatrix.median(IM::Eye{T}) where T
    m = size(IM, 1)
    m == 0 && throw(ArgumentError("median of an empty Eye is undefined"))
    Tout = promote_op(/, T, T)
    m == 1 && return one(Tout)
    m == 2 && return Tout(1) / Tout(2)
    return zero(Tout)
end

for f in (:permutedims, :triu, :triu!, :tril, :tril!, :inv)
    @eval ($f)(IM::Eye) = IM
end

function inv(DF::Diagonal{<:Any, Tf}) where {Tf <: Fill}
    value = getindex_value(DF.diag)
    return Diagonal(Fill(inv(value), size(DF, 1)))
end

# This agrees with diag for dense matrices in that any integer `k` is allowed.
# It disagrees with diag for Diagonal, which throws an error when `abs(k)` is too large.
function LinearAlgebra.diag(IM::Eye{T}, k::Integer) where T
    k == 0 && return LinearAlgebra.diag(IM)
    m = size(IM, 1)
    (k > m || k < -m) && return T[]
    return zeros(T, m - abs(k))
end

det(::Eye{T}) where T = one(T)
logdet(::Eye{T}) where T = log(one(T))
sqrt(IM::Eye) = IM

# There are a few trig functions `f` for which `f(IM)` cannot be computed
# at compile time. This is one way to solve this problem.
const _cschval = csch(1.0)
_mycsch(::Union{Type{Float64}, Type{Int}}) = _cschval
_mycsch(::Type{T}) where T = csch(one(T))
Base.csch(IM::Eye) = LinearAlgebra.Diagonal(Fill(_mycsch(eltype(IM)), size(IM, 1)))

==(x::Fill, y::Fill) = FillArrays.axes(x) == FillArrays.axes(y) &&
    FillArrays.getindex_value(x) == FillArrays.getindex_value(y)

==(IMa::Eye, IMb::Eye) = size(IMa, 1) == size(IMb, 1)

^(IM::Eye, p::Integer) = IM
/(AM::AbstractMatrix, IM::Eye) = IM * AM
/(AM::Eye, IM::Eye) = IM * AM
/(AM::DenseMatrix, IM::Eye) = IM * AM

function *(IM::Eye{T}, AV::AbstractVector{V}) where {T, V}
    _ = checkuniquedim(IM, AV)
    return convert(Vector{promote_op(*, T, V)}, AV)
end

const AbstractFillVecOrMat{V} = Union{AbstractFill{V, 1}, AbstractFill{V, 2}}

# Handles in particular Eye(n) * Ones(n), Eye(n) * Ones(n, n)
function *(IM::Eye{T}, AF::AT) where {T, AT <: AbstractFillVecOrMat{V}} where {V}
    Tout = typeof(one(T) * one(V))
    return convert(AbstractFill{Tout}, AF)
end

# Needed for disambiguation
*(IM::Eye{T}, Z::ZerosVecOrMat{V}) where {T, V} = Zeros{promote_op(*, T, V)}(size(Z))

function *(IM::Eye{T}, AM::AbstractMatrix{V}) where {T,V}
    junk = checkuniquedim(IM, AM)
    return convert(Matrix{promote_op(*, T, V)}, AM)
end

*(AM::AbstractMatrix{T}, IM::Eye{V}) where {T, V} = IM * AM
*(AM::Diagonal, IM::Eye) = IM * AM
*(IMa::Eye{T}, IMb::Eye{V}) where {T, V} = Eye{promote_op(*, T, V)}(size(IMa, 1))

# Kron with first argment either Diagonal or Eye
for (diagonaltype, A_jj, outputelement) in ((:(A::Diagonal{T}), :(A_jj = A[j, j]), :(A_jj * B[l,k]) ),
                        ( :(A::Eye{T}), :(nothing), :(B[l,k])))
    @eval begin
        function (kron)($diagonaltype, B::AbstractMatrix{S}) where {T<:Number, S<:Number}
            @assert ! Base.has_offset_axes(B)
            (mA, nA) = size(A); (mB, nB) = size(B)
            R = zeros(promote_op(*, T, S), mA * mB, nA * nB)
            m = 1
            for j = 1:nA
                $A_jj
                for k = 1:nB
                    for l = 1:mB
                        R[m] = $outputelement
                        m += 1
                    end
                    m += (nA - 1) * mB
                end
                m += mB
            end
            return R
        end
    end
end

# It is not clear if there is any advantage over kron(::AbstractMatrix, ::Diagonal)
function kron(A::AbstractMatrix{T}, B::Eye{S}) where {T<:Number, S<:Number}
    @assert ! Base.has_offset_axes(A)
    (mA, nA) = size(A); (mB, nB) = size(B)
    R = zeros(promote_op(*, T, S), mA * mB, nA * nB)
    m = 1
    for j = 1:nA
        for l = 1:mB
            for k = 1:mA
                R[m] = A[k,j]
                m += nB
            end
            m += 1
        end
        m -= nB
    end
    return R
end

kron(a::Eye{T}, b::Eye{T}) where {T<:Number} = Eye{T}(size(a, 1) * size(b, 1))
kron(a::Eye{T}, b::Eye{V}) where {T<:Number, V<:Number} = Eye{promote_op(*, T, V)}((size(a, 1) * size(b, 1)))

LinearAlgebra.eigvals(IM::Eye{T}) where T = diag(IM)
LinearAlgebra.eigvecs(IM::Eye) = IM # method for Diagonal returns a material matrix
LinearAlgebra.eigen(IM::Eye) = LinearAlgebra.Eigen(LinearAlgebra.eigvals(IM), LinearAlgebra.eigvecs(IM))

for f in (:eigmin, :eigmax)
    @eval ($f)(M::Diagonal{T,Tf}) where {T, Tf <: AbstractFill} = getindex_value(M.diag)
end

materialize(IM::Eye) = Matrix{eltype(IM)}(I, size(IM)) # This is a bit slower
Matrix(IM::Eye) = materialize(IM)

# For Eye{T}, the fallback method only materializes the diagonal.
Base.copymutable(IM::Eye) = Diagonal(ones(eltype(IM), size(IM, 1)))

+(IM::Eye{T}, s::UniformScaling) where T = Diagonal(Fill(one(T) + s.λ, size(IM, 1)))
+(s::UniformScaling, IM::Eye) = IM + s
-(IM::Eye{T}, s::UniformScaling) where T = Diagonal(Fill(one(T) - s.λ, size(IM, 1)))
-(s::UniformScaling, IM::Eye{T}) where T = Diagonal(Fill(s.λ - one(T), size(IM, 1)))

# Diagonal is already efficient. But, we use `Eye` to remove fatal method ambiguity intrduced
# by the methods below.
\(IMa::Eye{T}, IMb::Eye{V}) where {T, V} = IMa * IMb
\(AM::AbstractMatrix{T}, IM::Eye{V}) where {T, V} = convert(AbstractMatrix{promote_op(*, T, V)}, inv(AM))
/(IM::Eye, AM::AbstractMatrix) = \(AM, IM)
\(IM::Eye{V}, AM::AbstractMatrix{T}) where {T, V} = convert(AbstractMatrix{promote_op(*, T, V)}, AM)

\(IM::Eye, s::UniformScaling) = s.λ * IM
\(s::UniformScaling, IM::Eye) =  IM / s.λ

end # module


text:
struct FastSphericalHarmonicPlan{T} <: SphericalHarmonicPlan{T}
    RP::RotationPlan{T}
    BF::Vector{Butterfly{T}}
    p1::NormalizedLegendreToChebyshevPlan{T}
    p2::NormalizedLegendre1ToChebyshev2Plan{T}
    p1inv::ChebyshevToNormalizedLegendrePlan{T}
    p2inv::Chebyshev2ToNormalizedLegendre1Plan{T}
    B::Matrix{T}
end

function FastSphericalHarmonicPlan(A::Matrix{T}, L::Int; opts...) where T
    M, N = size(A)
    n = (N+1)÷2
    RP = RotationPlan(T, M-1)
    p1 = plan_normleg2cheb(A)
    p2 = plan_normleg12cheb2(A)
    p1inv = plan_cheb2normleg(A)
    p2inv = plan_cheb22normleg1(A)
    B = zero(A)
    Ce = Matrix{T}(I, M, M)
    Co = Matrix{T}(I, M, M)
    BF = Vector{Butterfly{T}}(undef, n-2)
    P = Progress(n-2, 0.1, "Pre-computing...", 43)
    for j = 1:2:n-2
        mul!(Ce, RP.layers[j])
        BF[j] = Butterfly(Ce, L; isorthogonal = true, opts...)
        next!(P)
    end
    for j = 2:2:n-2
        mul!(Co, RP.layers[j])
        BF[j] = Butterfly(Co, L; isorthogonal = true, opts...)
        next!(P)
    end
    FastSphericalHarmonicPlan(RP, BF, p1, p2, p1inv, p2inv, B)
end

FastSphericalHarmonicPlan(A::Matrix; opts...) = FastSphericalHarmonicPlan(A, floor(Int, log2(size(A, 1)+1)-6); opts...)

if VERSION ≥ v"0.7-"
    adjoint(P::FastSphericalHarmonicPlan) = Adjoint(P)
    transpose(P::FastSphericalHarmonicPlan) = Transpose(P)
end

function LAmul!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix)
    RP, BF, p1, p2, B = FP.RP, FP.BF, FP.p1, FP.p2, FP.B
    fill!(B, zero(eltype(B)))
    M, N = size(X)
    copyto!(B, 1, X, 1, 3M)
    @stepthreads for J = 2:N÷2
        mul_col_J!(B, BF[J-1], X, 2J)
        2J < N && mul_col_J!(B, BF[J-1], X, 2J+1)
    end

    mul_col_J!!(Y, p1, B, 1)
    @stepthreads for J = 2:4:N
        mul_col_J!!(Y, p2, B, J)
        J < N && mul_col_J!!(Y, p2, B, J+1)
    end
    @stepthreads for J = 4:4:N
        mul_col_J!!(Y, p1, B, J)
        J < N && mul_col_J!!(Y, p1, B, J+1)
    end
    Y
end
if VERSION < v"0.7-"
    function Base.At_mul_B!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix)
        RP, BF, p1inv, p2inv, B = FP.RP, FP.BF, FP.p1inv, FP.p2inv, FP.B
        copyto!(B, X)
        M, N = size(X)
        mul_col_J!!(Y, p1inv, B, 1)
        @stepthreads for J = 2:4:N
            mul_col_J!!(Y, p2inv, B, J)
            J < N && mul_col_J!!(Y, p2inv, B, J+1)
        end
        @stepthreads for J = 4:4:N
            mul_col_J!!(Y, p1inv, B, J)
            J < N && mul_col_J!!(Y, p1inv, B, J+1)
        end

        copyto!(B, Y)
        @stepthreads for J = 2:N÷2
            At_mul_B_col_J!(Y, BF[J-1], B, 2J)
            2J < N && At_mul_B_col_J!(Y, BF[J-1], B, 2J+1)
        end
        sph_zero_spurious_modes!(Y)
    end

    Base.Ac_mul_B!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix) = At_mul_B!(Y, FP, X)
else
    function LinearAlgebra.mul!(Y::Matrix, FPt::Transpose{T,<:FastSphericalHarmonicPlan}, X::Matrix) where T
        FP = parent(FPt)
        RP, BF, p1inv, p2inv, B = FP.RP, FP.BF, FP.p1inv, FP.p2inv, FP.B
        copyto!(B, X)
        M, N = size(X)
        mul_col_J!!(Y, p1inv, B, 1)
        @stepthreads for J = 2:4:N
            mul_col_J!!(Y, p2inv, B, J)
            J < N && mul_col_J!!(Y, p2inv, B, J+1)
        end
        @stepthreads for J = 4:4:N
            mul_col_J!!(Y, p1inv, B, J)
            J < N && mul_col_J!!(Y, p1inv, B, J+1)
        end

        copyto!(B, Y)
        @stepthreads for J = 2:N÷2
            At_mul_B_col_J!(Y, BF[J-1], B, 2J)
            2J < N && At_mul_B_col_J!(Y, BF[J-1], B, 2J+1)
        end
        sph_zero_spurious_modes!(Y)
    end

    LinearAlgebra.mul!(Y::Matrix, FPc::Adjoint{T,<:FastSphericalHarmonicPlan}, X::Matrix) where T =
        mul!(Y, transpose(parent(FPc)), X)
end

allranks(FP::FastSphericalHarmonicPlan) = mapreduce(allranks,vcat,FP.BF)


text:
module MonteCarloTreeSearch

@with_kw struct ExpectimaxOptions{S<:Function,V<:Function,Q<:Function,N<:Function,SA<:Function,R<:AbstractRNG}
    iterations::Int64   = 100
    maxdepth::Int64     = typemax(Int)
    rng::R              = Random.GLOBAL_RNG
    select::S           = MaxUCB(1.0)
    leafvalue::V        = RandomRollout(rng)
    initq::Q            = (m, s, a)->0.0
    initn::N            = (m, s, a)->0
    selectaction::SA    = maxtried
    keepstats::Bool     = false
end

inittree(options, model, state) = ExpectimaxTree{statetype(model), actiontype(model)}(options.iterations)

function simulate!(tree::ExpectimaxTree,
                   opt::ExpectimaxOptions,
                   m::MDP,
                   n::StateNode,
                   stepstogo::Int)

    # once depth is zero return
    if stepstogo == 0 || isterminal(m, state(node))
        return 0.0
    end

    # pick action
    sanode = opt.select(n)
    said = sanode.id

    # transition to a new state
    sp, r = generate_sr(m, state(n), action(sanode), rng)

    spid = get(tree.statemap, sp, 0)
    if spid == 0
        spn = insertnode!(tree, planner, sp)
        spid = spn.id
        q = r + discount(mdp) * opt.leafvalue(m, sp, stepstogo-1)
    else
        q = r + discount(mdp) * simulate!(tree, opt, m, StateNode(tree, spid), stepstogo-1)
    end

    if opt.keepstats
        recordvisit!(tree, said, spid)
    end

    tree.totaln[n.id] += 1
    tree.n[said] += 1
    tree.q[said] += (q - tree.q[said]) / tree.n[said] # moving average of Q value
    return q
end

maketree(options, model, state)

search(options, model, state)

struct RandomRollout{R<:AbstractRNG}
    rng::R
end

function (rr::RandomRollout)(m, s, stepstogo)
    sim = RolloutSimulator(rr.rng, stepstogo)
    policy = RandomPolicy(rr.rng, m, NothingUpdater())
    return POMDPs.simulate(sim, mdp, policy, s)
end

end # module


text:
fname = "../../../externals/Core.Math.Samples/data/Pejcic_318.csv"
data = readdlm(fname, separator = ',')

print(data)
# using DataFrames, CSV, Statistics, StatsBase
# data = CSV.read("Pejcic_318.csv", header=true, delim=',', copycols = true)

# println((data),"\n")


text:
using Transducers
# using Folds
using Profile
using DataStructures

function slurp_file(filename)
    lines = readlines(filename)
    template = collect(lines[1])
    rules = split.(lines[3:end], " -> ")
    template, rules
end

template, rules = slurp_file("14eg.txt")

function rules_to_f(rules)
    mapping = DefaultDict(Dict)
    for (pair, to_insert) in rules
        left, right = collect(pair)
        mapping[left][right] = only(to_insert)
    end
    function (left, right)
        if haskey(mapping, left) && haskey(mapping[left], right)
            mapping[left][right]
        else
        end
    end
end

f = rules_to_f(rules)

function insertions(s, rule_f)
    result = s |> Consecutive(2, 1) |> MapCat(rule_f) |> collect
    push!(result, s[length(s)])
    result
end

function insertions2(s, rule_f)
    left = s[1]
    result = Vector()
    for right in s[2:end]
        push!(result, left)
        to_insert = rule_f(left, right)
        if !isnothing(to_insert)
            push!(result, to_insert)
        end
        left = right
    end
    push!(result, s[end])
    result
end

insertions2(template, f)
insertions2(template, f) == collect("NCNBCHB")

function solve(filename, steps)
    template, rules = slurp_file(filename)
    rule_f = rules_to_f(rules)
    for c = 1:steps
        template = insertions2(template, rule_f)
    end
    vals = values(counter(template))
    max(vals...) - min(vals...)
end

solve1(filename) = solve(filename, 10)

solve1("14.txt") == 2975

to_pair_counts(template) = counter(collect(Consecutive(2,1)(template)))

function to_letter_counts(counts)
    result = counter(Char)
    for ((l,r), count) in counts
        result[l] += count
        result[r] += count
    end
    map!(x -> iseven(x) ? x/2 : (x+1)/2, values(result))
    result
end

function solve2(filename, num_steps)
    template, rules = slurp_file(filename)
    rule_f = rules_to_f(rules)
    counts = to_pair_counts(template)
    for c = 1:num_steps
        new_counts = counter(Tuple{Char, Char})
        for ((l,r), count) in counts
            new = rule_f(l,r)
            new_counts[(l, new)] += count
            new_counts[(new, r)] += count
        end
        counts = new_counts
    end
    letter_counts = to_letter_counts(counts)
    vals = values(letter_counts)
    max(vals...) - min(vals...)
end

solve2("14.txt", 10) == 2975
solve2("14.txt", 40)


text:
"""
    Econometrics

    Econometrics in Julia.
"""
module Econometrics
    using Base.Iterators: flatten
    using LinearAlgebra: bunchkaufman, bunchkaufman!, cholesky!, diag, diagm, Diagonal,
                         Hermitian, I, LowerTriangular, qr, UpperTriangular
    using DataFrames: AbstractDataFrame, categorical!, DataFrame, dropmissing,
                      # CategoricalArrays
                      AbstractCategoricalVector, categorical, levels, isordered
    using Distributions: cdf, ccdf, Categorical, Chisq, FDist, logpdf, Logistic, Normal,
                         pdf, TDist,
                         # Statistics
                         mean, quantile, var,
                         # StatsFuns
                         softmax
    using FillArrays: Fill
    using Optim: hessian!, optimize, minimizer, TwiceDifferentiable,
                 # ForwardDiff
                 ForwardDiff.Dual, ForwardDiff.value
    using Parameters: @unpack
    using Printf: @sprintf
    using StatsModels: AbstractContrasts, AbstractTerm, apply_schema, CategoricalTerm,
                       ConstantTerm, ContrastsMatrix, DummyCoding, @formula, FormulaTerm,
                       FunctionTerm, InterceptTerm, MatrixTerm, modelcols, schema, terms,
                       termvars,
                       # StatsBase
                       aic, aicc, bic, harmmean, FrequencyWeights, CoefTable,
                       ConvergenceException
    import Base: show
    import StatsModels: hasintercept, implicit_intercept,
                        # StatsBase
                        coef, coefnames, coeftable, confint, deviance, islinear, nulldeviance,
                        loglikelihood, nullloglikelihood, score, nobs, dof, dof_residual,
                        mss, rss, informationmatrix, vcov, stderror, weights, isfitted, fit,
                        fit!, r2, adjr2, fitted, response, meanresponse, modelmatrix, leverage,
                        residuals, predict, predict!, dof_residual, RegressionModel,
                        StatsBase.params

    foreach(file -> include(joinpath(dirname(@__DIR__), "src", "$file.jl")),
            ["structs", "transformations", "formula", "solvers", "main", "statsbase", "wald"])

    export categorical!, DataFrame, @formula, DummyCoding, aic, aicc, bic, coef, coefnames,
           coeftable, confint, deviance, islinear, nulldeviance, loglikelihood,
           nullloglikelihood, score, nobs, dof, mss, rss, informationmatrix, vcov, stderror,
           weights, isfitted, fit, fit!, r2, adjr2, fitted, response, meanresponse,
           modelmatrix, leverage, residuals, predict, predict!, dof_residual,
           params, EconometricModel, absorb, between, PID, TID
end


text:
import Pkg
if VERSION >= v"1.1.0"
    function pkg_version(_module::Module)
        Pkg.Types.read_package(joinpath(dirname(dirname(pathof(_module))), "Project.toml")).version
    end
else
    function pkg_version(_module::Module)
        pkg_id = Base.PkgId(_module)
        env = Pkg.Types.Context().env
        pkg_info = Pkg.Types.manifest_info(env, pkg_id.uuid)
        return VersionNumber(pkg_info["version"])
    end
end

if pkg_version(CUDAnative) >= v"2.1.0"
    for op = [:add, # :sub, :max, :min,
              ]
        @eval function $(Symbol("scatter_", op, "!"))(ys::CuMatrix{T}, us::CuArray{T}, xs::CuArray{Int}) where T
            function kernel!(ys, us, xs)
                li = threadIdx().y + (blockIdx().y - 1) * blockDim().y
                i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

                @inbounds if li <= length(xs) && i <= size(ys, 1)
                    ind = Tuple(CartesianIndices(xs)[li])
                    CUDAnative.$(Symbol("atomic_", op, "!"))(
                        pointer(ys,
                                Base._to_linear_index(ys,
                                                      i, xs[li]
                                                      )
                                ),
                        us[i, ind...]
                    )
                end

                return
            end

            max_threads = 256
            thread_x = min(max_threads, size(ys, 1))
            thread_y = min(max_threads ÷ thread_x, length(xs))
            threads = (thread_x, thread_y)
            blocks = ceil.(Int, (size(ys, 1), length(xs)) ./ threads)
            @cuda blocks=blocks threads=threads kernel!(ys, us, xs)
            return ys
        end

        @eval function $(Symbol("scatter_", op, "!"))(ys::CuArray{T}, us::CuArray{T}, xs::CuArray{<:Tuple}) where T
            function kernel!(ys, us, xs)
                li = threadIdx().y + (blockIdx().y - 1) * blockDim().y
                i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

                @inbounds if li <= length(xs) && i <= size(ys, 1)
                    ind = Tuple(CartesianIndices(xs)[li])
                    CUDAnative.$(Symbol("atomic_", op, "!"))(
                        pointer(ys,
                                Base._to_linear_index(ys,
                                                      i, xs[li]...
                                                      )
                                ),
                        us[i, ind...]
                    )
                end

                return
            end

            max_threads = 256
            thread_x = min(max_threads, size(ys, 1))
            thread_y = min(max_threads ÷ thread_x, length(xs))
            threads = (thread_x, thread_y)
            blocks = ceil.(Int, (size(ys, 1), length(xs)) ./ threads)
            @cuda blocks=blocks threads=threads kernel!(ys, us, xs)
            return ys
        end
    end
else
    function scatter_add!(ys::CuMatrix{T}, us::CuArray{T}, xs::CuArray{Int}) where T
        function kernel!(ys::CuDeviceArray{T}, us::CuDeviceArray{T}, xs)
            xi = threadIdx().x + (blockIdx().x - 1) * blockDim().x

            if xi <= size(ys, 1)
                @inbounds for i = 1:length(xs)
                    ind = Tuple(CartesianIndices(xs)[i])
                    ys[xi, xs[i]] += us[xi, ind...]
                end
            end

            return
        end

        max_threads = 256
        threads = max_threads
        blocks = ceil(Int, size(ys, 1) / threads)
        CuArrays.@cuda blocks=blocks threads=threads kernel!(ys, us, xs)
        return ys
    end

    function scatter_add!(ys::CuArray{T}, us::CuArray{T}, xs::CuArray{<:Tuple}) where T
        function kernel!(ys::CuDeviceArray{T}, us::CuDeviceArray{T}, xs)
            xi = threadIdx().x + (blockIdx().x - 1) * blockDim().x

            if xi <= size(ys, 1)
                @inbounds for i = 1:length(xs)
                    ind = Tuple(CartesianIndices(xs)[i])
                    ys[xi, xs[i]...] += us[xi, ind...]
                end
            end

            return
        end

        max_threads = 256
        threads = max_threads
        blocks = ceil(Int, size(ys, 1) / threads)
        CuArrays.@cuda blocks=blocks threads=threads kernel!(ys, us, xs)
        return ys
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "CMake"
version = v"3.17.1"

# Collection of sources required to build CMake
sources = [
    ArchiveSource("https://github.com/Kitware/CMake/releases/download/v$(version)/cmake-$(version).tar.gz",
                  "3aa9114485da39cbd9665a0bfe986894a282d5f0882b1dea960a739496620727"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/cmake-*/

cmake -DCMAKE_INSTALL_PREFIX=$prefix

make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.
# CMake is in C++ and it exports the C++ string ABIs, but when compiling it with
# the C++03 string ABI it seems to ignore our request, so let's just build for
# the C++11 string ABI.
platforms = [
    Platform("i686", "linux"; libc="glibc", cxxstring_abi = "cxx11"),
    Platform("x86_64", "linux"; libc="glibc", cxxstring_abi = "cxx11"),
    Platform("x86_64", "linux"; libc="musl", cxxstring_abi = "cxx11"),
]

# platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("cmake", :cmake),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("OpenSSL_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: vchuravy/Yggdrasil

max_issues_repo_name: vchuravy/Yggdrasil

max_forks_repo_name: vchuravy/Yggdrasil


text:
#=
Problem 4
A palindromic number reads the same both ways. The largest palindrome made from the product of two 2-digit numbers is 9009 = 91 × 99.
Find the largest palindrome made from the product of two 3-digit numbers.
=#

function isPal(x)
  s = string(x)
  l = length(s)
  for i=1:l
    if (s[i] != s[l-i+1])
      return false
    end
  end
  return true
end

maxPal = 0
p = 0
for i=1:999, j=1:999
  p = i*j
  if (isPal(p) && p > maxPal)
    maxPal = p
  end
end

println("The largest palindrome is: $(maxPal)")

max_stars_repo_name: mebrunet/project-euler

max_issues_repo_name: mebrunet/project-euler

max_forks_repo_name: mebrunet/project-euler


text:
function lagrange_poly_basis(nodes, evalPt)
    order = length(nodes)
    poly_coefs = zeros(order)
    for n in 1:order
        notn = [1:n-1;n+1:order]
        poly_coefs[n] = prod(evalPt.-nodes[notn])./prod(nodes[n].-nodes[notn])
    end
    return poly_coefs
end

function lagrange_interpolant()

end
function gauss_lobatto_quadrature(fun::Function,a::Float64,b::Float64,n_evals::Int)
    (b<=a)&&throw(DomainError("must have a<b"))
    (n_evals<1)&&throw(DomainError("n_evals must be > 0"))
    if n_evals > 1
        # the GL nodes
        # nodes in [-1,1]
        nodes = Jacobi.zglj(n_evals, 0, 0)
        weights = Jacobi.wglj(nodes,0,0)*(b-a)/2
        # shift nodes to [a,b]
        nodes *= 0.5*(b-a)
        nodes .+= 0.5*(a+b)
    else
        nodes = 0.5*(a+b)
        weights = (b-a)
    end
    fun_vals = fun.(nodes)
    quad = LinearAlgebra.dot(fun_vals,weights)
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
module Lasso
module Util
    # Extract fields from object into function locals
    # See https://github.com/JuliaLang/julia/issues/9755
    macro extractfields(from, fields...)
        esc(Expr(:block, [:($(fields[i]) = $(from).$(fields[i])) for i = 1:length(fields)]...))
    end
    export @extractfields
end

include("FusedLasso.jl")
include("TrendFiltering.jl")

using Reexport, LinearAlgebra, SparseArrays, Random, .Util, MLBase
import Random: Sampler
@reexport using GLM, Distributions, .FusedLassoMod, .TrendFiltering
using GLM: FPVector
export RegularizationPath, LassoPath, GammaLassoPath, NaiveCoordinateDescent,
       CovarianceCoordinateDescent, fit, fit!, coef, predict,
       minAICc, hasintercept, dof, aicc, distfun, linkfun, cross_validate_path

## HELPERS FOR SPARSE COEFFICIENTS

struct SparseCoefficients{T} <: AbstractVector{T}
    coef::Vector{T}              # Individual coefficient values
    coef2predictor::Vector{Int}  # Mapping from indices in coef to indices in original X
    predictor2coef::Vector{Int}  # Mapping from indices in original X to indices in coef

    SparseCoefficients{T}(n::Int) where {T} = new(T[], Int[], zeros(Int, n))
end

function LinearAlgebra.mul!(out::Vector, X::Matrix, coef::SparseCoefficients{T}) where T
    fill!(out, zero(eltype(out)))
    @inbounds for icoef = 1:nnz(coef)
        ipred = coef.coef2predictor[icoef]
        c = coef.coef[icoef]
        @simd for i = 1:size(X, 1)
            out[i] += c*X[i, ipred]
        end
    end
    out
end

function LinearAlgebra.mul!(out::Vector, X::SparseMatrixCSC, coef::SparseCoefficients{T}) where T
    @extractfields X colptr rowval nzval
    fill!(out, zero(eltype(out)))
    @inbounds for icoef = 1:nnz(coef)
        ipred = coef.coef2predictor[icoef]
        c = coef.coef[icoef]
        @simd for i = colptr[ipred]:colptr[ipred+1]-1
            out[rowval[i]] += c*nzval[i]
        end
    end
    out
end

function LinearAlgebra.dot(x::Vector{T}, coef::SparseCoefficients{T}) where T
    v = 0.0
    @inbounds @simd for icoef = 1:nnz(coef)
        v += x[coef.coef2predictor[icoef]]*coef.coef[icoef]
    end
    v
end

Base.size(x::SparseCoefficients) = (length(x.predictor2coef),)
SparseArrays.nnz(x::SparseCoefficients) = length(x.coef)
Base.getindex(x::SparseCoefficients{T}, ipred::Int) where {T} =
    x.predictor2coef[ipred] == 0 ? zero(T) : x.coef[x.predictor2coef[ipred]]

function Base.setindex!(A::Matrix{T}, coef::SparseCoefficients, rg::UnitRange{Int}, i::Int) where T
    A[:, i] = zero(T)
    for icoef = 1:nnz(coef)
        A[rg[coef.coef2predictor[icoef]], i] = coef.coef[icoef]
    end
    A
end

function Base.copyto!(x::SparseCoefficients, y::SparseCoefficients)
    length(x) == length(y) || throw(DimensionMismatch())
    n = length(y.coef)
    resize!(x.coef, n)
    resize!(x.coef2predictor, n)
    copyto!(x.coef, y.coef)
    copyto!(x.coef2predictor, y.coef2predictor)
    copyto!(x.predictor2coef, y.predictor2coef)
    x
end

# Add a new coefficient to x, returning its index in x.coef
function addcoef!(x::SparseCoefficients{T}, ipred::Int) where T
    push!(x.coef, zero(T))
    push!(x.coef2predictor, ipred)
    coefindex = nnz(x)
    x.predictor2coef[ipred] = coefindex
end

# Add newcoef to column i of coefs
function addcoefs!(coefs::SparseMatrixCSC, newcoef::SparseCoefficients, i::Int)
    n = nnz(coefs)
    nzval = coefs.nzval
    rowval = coefs.rowval
    resize!(nzval, n+nnz(newcoef))
    resize!(rowval, n+nnz(newcoef))
    @inbounds for ipred = 1:length(newcoef.predictor2coef)
        icoef = newcoef.predictor2coef[ipred]
        if icoef != 0
            cval = newcoef.coef[icoef]
            if cval != 0
                n += 1
                nzval[n] = cval
                rowval[n] = ipred
            end
        end
    end
    resize!(nzval, n)
    resize!(rowval, n)
    coefs.colptr[i+1:end] .= n+1
end

## COEFFICIENT ITERATION IN SEQUENTIAL OR RANDOM ORDER
struct RandomCoefficientIterator
    rng::MersenneTwister
    rg::Sampler
    coeforder::Vector{Int}
end
const RANDOMIZE_DEFAULT = true

function RandomCoefficientIterator()
    rng = MersenneTwister(1337)
    RandomCoefficientIterator(rng, Sampler(rng, 1:2), Int[])
end

const CoefficientIterator = Union{UnitRange{Int},RandomCoefficientIterator}

# Iterate over coefficients in random order
function Base.iterate(x::RandomCoefficientIterator)
    if !isempty(x.coeforder)
        @inbounds for i = length(x.coeforder):-1:2
            j = rand(x.rng, x.rg)
            x.coeforder[i], x.coeforder[j] = x.coeforder[j], x.coeforder[i]
        end
        x.coeforder[1], 2
    else
        nothing
    end
end

Base.iterate(x::RandomCoefficientIterator, i) = (i > length(x.coeforder)) ? nothing : (x.coeforder[i], i += 1)

# Add an additional coefficient and return a new CoefficientIterator
function addcoef(x::RandomCoefficientIterator, icoef::Int)
    push!(x.coeforder, icoef)
    RandomCoefficientIterator(x.rng, Sampler(x.rng, 1:length(x.coeforder)), x.coeforder)
end
addcoef(x::UnitRange{Int}, icoef::Int) = 1:length(x)+1

abstract type RegularizationPath{S<:Union{LinearModel,GeneralizedLinearModel},T} <: RegressionModel end
## LASSO PATH

mutable struct LassoPath{S<:Union{LinearModel,GeneralizedLinearModel},T} <: RegularizationPath{S,T}
    m::S
    nulldev::T                    # null deviance
    nullb0::T                     # intercept of null model, if one was fit
    λ::Vector{T}                  # shrinkage parameters
    autoλ::Bool                   # whether λ is automatically determined
    Xnorm::Vector{T}              # original squared norms of columns of X before standardization
    pct_dev::Vector{T}            # percent deviance explained by each model
    coefs::SparseMatrixCSC{T,Int} # model coefficients
    b0::Vector{T}                 # model intercepts
    niter::Int                    # number of coordinate descent iterations

    LassoPath{S,T}(m, nulldev::T, nullb0::T, λ::Vector{T}, autoλ::Bool, Xnorm::Vector{T}) where {S,T} =
        new(m, nulldev, nullb0, λ, autoλ, Xnorm)
end

function Base.show(io::IO, path::RegularizationPath)
    prefix = isa(path.m, GeneralizedLinearModel) ? string(typeof(distfun(path)).name.name, " ") : ""
    pathsize = size(path)
    println(io, "$(prefix)$(typeof(path).name.name) ($(pathsize[2])) solutions for $(pathsize[1]) predictors in $(path.niter) iterations):")

    if isdefined(path, :coefs)
        coefs = path.coefs
        ncoefs = zeros(Int, size(coefs, 2))
        for i = 1:size(coefs, 2)-1
            ncoefs[i] = coefs.colptr[i+1] - coefs.colptr[i]
        end
        ncoefs[end] = nnz(coefs) - coefs.colptr[size(coefs, 2)] + 1
        show(io, CoefTable(Union{Vector{Int},Vector{Float64}}[path.λ, path.pct_dev, ncoefs], ["λ", "pct_dev", "ncoefs"], []))
    else
        print(io, "    (not fit)")
    end
end

## MODEL CONSTRUCTION

# Controls early stopping criteria with automatic λ
const MIN_DEV_FRAC_DIFF = 1e-5
const MAX_DEV_FRAC = 0.999

# Compute automatic λ values based on λmax and λminratio
function computeλ(λmax, λminratio, α, nλ)
    λmax /= α
    logλmax = log(λmax)
    exp.(range(logλmax, stop=logλmax + log(λminratio), length=nλ))
end

# Compute automatic λ values based on X'y and λminratio without weights ω
function computeλ(Xy::Vector{T}, λminratio, α, nλ, ω::Nothing) where T
    λmax = zero(T)
    for i = eachindex(Xy)
        x = abs(Xy[i])
        if x > λmax
            λmax = x
        end
    end
    computeλ(λmax, λminratio, α, nλ)
end

# Compute automatic λ values based on X'y and λminratio with specific weights ω
function computeλ(Xy::Vector{T}, λminratio, α, nλ, ω::Vector) where T
    λmax = zero(T)
    for i = eachindex(Xy)
        if ω[i] > 0
            x = abs(Xy[i]) / ω[i]
            if x > λmax
                λmax = x
            end
        end
    end
    computeλ(λmax, λminratio, α, nλ)
end

function weightedresiduals(model::LinearModel)
    r = model.rr
    y = r.y
    mu = r.mu
    if isempty(r.wts)
        y - mu
    else
        wts = r.wts
        resid = similar(y)
        @simd for i = eachindex(resid,y,mu,wts)
            @inbounds resid[i] = (y[i] - mu[i]) * wts[i]
        end
        resid
    end
end

# NOTE: mutates nullmodel
maxλ!(nullmodel, λ::Vector, X::AbstractMatrix{T}, λminratio, α, nλ, ω) where T =
    convert(Vector{T}, λ)

function maxλ!(nullmodel::LinearModel, λ::Nothing,
    X::AbstractMatrix{T}, λminratio, α, nλ, ω) where T

    # muscratch = nullmodel.rr.mu.*nullmodel.rr.wts
    Xy = X'*weightedresiduals(nullmodel)

    computeλ(Xy, λminratio, α, nλ, ω)
end

function maxλ!(nullmodel::GeneralizedLinearModel, λ::Nothing,
    X::AbstractMatrix{T}, λminratio, α, nλ, ω) where T

    Xy = X'*broadcast!(*, nullmodel.rr.wrkresid, nullmodel.rr.wrkresid, nullmodel.rr.wrkwt)
    computeλ(Xy, λminratio, α, nλ, ω)
end

# rescales A so that it sums to base
rescale(A, base) = A * (base / sum(A))

# null model's X includes unpenalized (ω_j==0) variables
nullX(X::AbstractMatrix{T}, intercept::Bool, ω::Nothing) where T =
    ones(T, size(X,1), ifelse(intercept, 1, 0))

function nullX(X::AbstractMatrix{T}, intercept::Bool, ω::Vector) where T
    ixunpenalized = findall(iszero, ω)
    [ones(T, size(X,1), ifelse(intercept, 1, 0)) view(X, :, ixunpenalized)]
end

# version for Linear models
function build_model(X::AbstractMatrix{T}, y::FPVector, d::Normal, l::IdentityLink,
                     lp::LinPred, λminratio::Real, λ::Union{Vector,Nothing},
                     wts::Union{FPVector,Nothing}, offset::Vector, α::Real, nλ::Int,
                     ω::Union{Vector,Nothing}, intercept::Bool, irls_tol::Real, dofit::Bool) where T
    mu = isempty(offset) ? y : y + offset
    nullmodel = LinearModel(LmResp{typeof(y)}(fill!(similar(y), 0), offset, wts, y),
                            GLM.cholpred(nullX(X, intercept, ω), false))
    fit!(nullmodel)
    nulldev = deviance(nullmodel)
    nullb0 = intercept ? coef(nullmodel)[1] : zero(T)

    λ = maxλ!(nullmodel, λ, X, λminratio, α, nλ, ω)

    # mu is just a placeholder here, setting it to nullmodel.mu causes problems
    model = LinearModel(LmResp{typeof(y)}(mu, offset, wts, y), lp)

    (model, nulldev, nullb0, λ)
end

# version for GLMs
function build_model(X::AbstractMatrix{T}, y::FPVector, d::UnivariateDistribution, l::Link,
                     lp::LinPred, λminratio::Real, λ::Union{Vector,Nothing},
                     wts::Union{FPVector,Nothing}, offset::Vector, α::Real, nλ::Int,
                     ω::Union{Vector, Nothing}, intercept::Bool, irls_tol::Real, dofit::Bool) where T
    # Fit to find null deviance
    # Maybe we should reuse this GlmResp object?
    nullmodel = fit(GeneralizedLinearModel, nullX(X, intercept, ω), y, d, l;
                    wts=wts, offset=offset, convTol=irls_tol, dofit=dofit)
    nulldev = deviance(nullmodel)
    nullb0 = intercept ? coef(nullmodel)[1] : zero(T)

    λ = maxλ!(nullmodel, λ, X, λminratio, α, nλ, ω)

    rr = GlmResp(y, d, l, offset, wts)
    model = GeneralizedLinearModel(rr, lp, false)

    (model, nulldev, nullb0, λ)
end

defaultalgorithm(d::Normal, l::IdentityLink, n::Int, p::Int) = p > 5n ? NaiveCoordinateDescent : CovarianceCoordinateDescent
defaultalgorithm(d::UnivariateDistribution, l::Link, n::Int, p::Int) = NaiveCoordinateDescent

# following glmnet rescale penalty factors to sum to the number of coefficients
initpenaltyfactor(penalty_factor::Nothing,p::Int,::Bool) = nothing
initpenaltyfactor(penalty_factor::Vector,p::Int,standardizeω::Bool) =
    standardizeω ? rescale(penalty_factor, p) : penalty_factor

# Standardize predictors if requested
function standardizeX(X::AbstractMatrix{T}, standardize::Bool) where T
    if standardize
        Xnorm = vec(convert(Matrix{T},std(X; dims=1, corrected=false)))
        if any(x -> x == zero(T), Xnorm)
            @warn("""One of the predicators (columns of X) is a constant, so it can not be standardized.
                  To include a constant predicator set standardize = false and intercept = false""")
        end
        for i = 1:length(Xnorm)
            @inbounds Xnorm[i] = 1/Xnorm[i]
        end
        X = X .* transpose(Xnorm)
    else
        Xnorm = T[]
    end

    X, Xnorm
end

function StatsBase.fit(::Type{LassoPath},
                       X::AbstractMatrix{T}, y::V, d::UnivariateDistribution=Normal(),
                       l::Link=canonicallink(d);
                       wts::Union{FPVector,Nothing}=ones(T, length(y)),
                       offset::V=similar(y, 0),
                       α::Number=one(eltype(y)), nλ::Int=100,
                       λminratio::Number=ifelse(size(X, 1) < size(X, 2), 0.01, 1e-4),
                       λ::Union{Vector,Nothing}=nothing, standardize::Bool=true,
                       intercept::Bool=true,
                       algorithm::Type=defaultalgorithm(d, l, size(X, 1), size(X, 2)),
                       dofit::Bool=true,
                       irls_tol::Real=1e-7, randomize::Bool=RANDOMIZE_DEFAULT,
                       maxncoef::Int=min(size(X, 2), 2*size(X, 1)),
                       penalty_factor::Union{Vector,Nothing}=nothing,
                       standardizeω::Bool=true,
                       fitargs...) where {T<:AbstractFloat,V<:FPVector}
    size(X, 1) == size(y, 1) || DimensionMismatch("number of rows in X and y must match")
    n = length(y)
    length(wts) == n || error("length(wts) = $(length(wts)) should be 0 or $n")

    X, Xnorm = standardizeX(X, standardize)

    # Lasso initialization
    α = convert(T, α)
    λminratio = convert(T, λminratio)
    coefitr = randomize ? RandomCoefficientIterator() : (1:0)

    # penalty_factor (ω) defaults to a vector of ones
    ω = initpenaltyfactor(penalty_factor, size(X, 2), standardizeω)

    cd = algorithm{T,intercept,typeof(X),typeof(coefitr),typeof(ω)}(X, α, maxncoef, 1e-7, coefitr, ω)

    # GLM response initialization
    autoλ = λ == nothing
    model, nulldev, nullb0, λ = build_model(X, y, d, l, cd, λminratio, λ, wts .* T(1/sum(wts)),
                                            Vector{T}(offset), α, nλ, ω, intercept, irls_tol, dofit)

    # Fit path
    path = LassoPath{typeof(model),T}(model, nulldev, nullb0, λ, autoλ, Xnorm)
    if dofit
        fit!(path; irls_tol=irls_tol, fitargs...)
    end
    path
end

StatsBase.nobs(path::RegularizationPath) = length(path.m.rr.y)

dispersion_parameter(path::RegularizationPath) = GLM.dispersion_parameter(distfun(path))

function StatsBase.loglikelihood(path::RegularizationPath{S,T}) where {S<:LinearModel,T}
    n = nobs(path)
    -0.5.*n.*log.(deviance(path))
end

function StatsBase.loglikelihood(path::RegularizationPath{S,T}) where {S,T}
    n = nobs(path)
    -0.5.*n.*deviance(path)
end

"""
    dof(path::RegularizationPath)

Approximates the degrees-of-freedom in each segment of the path as the number of non zero coefficients
plus a dispersion parameter when appropriate.
Note that for GammaLassoPath this may be a crude approximation, as gamlr does this differently.
"""
function StatsBase.dof(path::RegularizationPath)
    nλ = length(path.λ)
    βs = coef(path)
    df = zeros(Int,nλ)
    for s=1:nλ
        df[s] = sum(βs[:,s].!=0)
    end

    if dispersion_parameter(path)
        # add one for dispersion_parameter
        df.+=1
    end

    df
end

function infocrit(d::T,l::F,n,k) where {T,F}
    if d + one(T) > n
        floatmax(F)
    else
        -2l + (k*d + k*d*(d+1)/(n-d-1))
    end
end

function StatsBase.aicc(path::RegularizationPath;k=2)
    dfs = dof(path)
    ls = loglikelihood(path)
    n = nobs(path)
    broadcast((d,l)->infocrit(d,l,n,k), dfs, ls)
end

minAIC(path::RegularizationPath)=argmin(aic(path))
minAICc(path::RegularizationPath;k=2)=argmin(aicc(path;k=k))
minBIC(path::RegularizationPath)=argmin(bic(path))

hasintercept(path::RegularizationPath) = hasintercept(path.m.pp)

"""
size(path) returns (p,nλ) where p is the number of coefficients (including
any intercept) and nλ is the number of path segments.
If model was only initialized but not fit, returns (p,1).
"""
function Base.size(path::RegularizationPath)
  if isdefined(path,:coefs)
    p,nλ = size(path.coefs)
  else
    X = path.m.pp.X
    p = size(X,2)
    nλ = 1
  end

  if hasintercept(path)
    p += 1
  end

  p,nλ
end

"""
    coef(path::RegularizationPath; select=:AICc)

Returns a p by nλ coefficient array where p is the number of
coefficients (including any intercept) and nλ is the number of path segments,
or a selected segment's coefficients.

If model was only initialized but not fit, returns a p vector of zeros.
Consistent with StatsBase.coef, if the model has an intercept it is included.

# Example:
```julia
  m = fit(LassoPath,X,y)
  coef(m; select=:CVmin)
```

# Keywords
- `select=:all` returns a p by nλ matrix of coefficients
- `select=:AIC` selects the AIC minimizing segment
- `select=:AICc` selects the corrected AIC minimizing segment
- `select=:BIC` selects the BIC minimizing segment
- `select=:CVmin` selects the segment that minimizing out-of-sample mean squared
    error from cross-validation with `nCVfolds` random folds.
- `select=:CV1se` selects the segment whose average OOS deviance is no more than
    1 standard error away from the one minimizing out-of-sample mean squared
    error from cross-validation with `nCVfolds` random folds.
- `nCVfolds=10` number of cross-validation folds
- `kwargs...` are passed to [`minAICc(::RegularizationPath)`](@ref) or to
    [`cross_validate_path(::RegularizationPath)`](@ref)
"""
function StatsBase.coef(path::RegularizationPath; select=:all, nCVfolds=10, kwargs...)
    if !isdefined(path,:coefs)
        X = path.m.pp.X
        p,nλ = size(path)
        if select == :all
            return zeros(eltype(X),p,nλ)
        else
            return zeros(eltype(X),p)
        end
    end

    if select == :all
        if hasintercept(path)
            vcat(path.b0',path.coefs)
        else
            path.coefs
        end
    elseif select in [:AICc, :AIC, :BIC]
        if select == :AICc
            seg = minAICc(path; kwargs...)
        elseif select == :AIC
            seg = minAIC(path)
        elseif select == :BIC
            seg = minBIC(path)
        end

        if hasintercept(path)
            vec(vcat(path.b0[seg],path.coefs[:,seg]))
        else
            path.coefs[:,seg]
        end
    elseif select == :CVmin || select == :CV1se
        gen = Kfold(length(path.m.rr.y),nCVfolds)
        segCV = cross_validate_path(path;gen=gen,select=select, kwargs...)
        if hasintercept(path)
            vec(vcat(path.b0[segCV],path.coefs[:,segCV]))
        else
            path.coefs[:,segCV]
        end
    else
        error("unknown selector $select")
    end
end

"link of underlying GLM"
GLM.linkfun(path::RegularizationPath{M}) where {M<:LinearModel} = IdentityLink()
GLM.linkfun(path::RegularizationPath{GeneralizedLinearModel{GlmResp{V,D,L},L2}}) where {V<:FPVector,D<:UnivariateDistribution,L<:Link,L2<:GLM.LinPred} = L()

## Prediction function for GLMs
function StatsBase.predict(path::RegularizationPath, newX::AbstractMatrix{T}; offset::FPVector=T[], select=:all) where {T<:AbstractFloat}
    # add an interecept to newX if the model has one
    if hasintercept(path)
        newX = [ones(eltype(newX),size(newX,1),1) newX]
    end

    # calculate etas for each obs x segment
    eta = newX * coef(path;select=select)

    # get model
    mm = path.m

    # adjust for any offset
    if length(mm.rr.offset) > 0
        length(offset) == size(newX, 1) ||
            throw(ArgumentError("fit with offset, so `offset` kw arg must be an offset of length `size(newX, 1)`"))
        broadcast!(+, eta, eta, offset)
    else
        length(offset) > 0 && throw(ArgumentError("fit without offset, so value of `offset` kw arg does not make sense"))
    end

    if typeof(mm) <: LinearModel
        eta
    else
        # invert all etas to mus
        μ(η) = linkinv(linkfun(path), η)
        map(μ, eta)
    end
end

"distribution of underlying GLM"
distfun(path::RegularizationPath{M}) where {M<:LinearModel} = Normal()
distfun(path::RegularizationPath) = path.m.rr.d

"deviance at each segment of the path for the fitted model and data"
StatsBase.deviance(path::RegularizationPath) = (1 .- path.pct_dev) .* path.nulldev # * nobs(path)

"""
deviance at each segement of the path for (potentially new) data X and y
select=:all or :AICc like in coef()
"""
function StatsBase.deviance(path::RegularizationPath, X::AbstractMatrix{T}, y::V;
                    offset::FPVector=T[], select=:all,
                    wts::FPVector=ones(T, length(y))) where {T<:AbstractFloat,V<:FPVector}
    μ = predict(path, X; offset=offset, select=select)
    deviance(path, y, μ, wts)
end

"""
deviance at each segment of the path for (potentially new) y and predicted values μ
"""
function StatsBase.deviance(path::RegularizationPath, y::V, μ::AbstractArray{T},
                wts::FPVector=ones(T, length(y))) where {T<:AbstractFloat,V<:FPVector}
    # get model specs from path
    d = distfun(path)

    # rescale weights
    wts = wts .* convert(T, 1/sum(wts))

    # closure for deviance of a single observation
    dev(ys,μs,ws) = ws * devresid(d, ys, μs)

    # deviances of all obs x segment
    devresidv = broadcast(dev,y,μ,wts)

    # deviance is just their sum
    if size(μ,2) > 1
        vec(sum(devresidv,dims=1))
    else
        sum(devresidv)
    end
end

include("coordinate_descent.jl")
include("gammalasso.jl")
include("cross_validation.jl")

end


text:
@testset "Controller: Objective" begin
    s = get_simulation("quadruped", "flat_2D_lc", "flat")
    model = s.model
    env = s.env

    H = 10
    obj = ContactImplicitMPC.TrackingObjective(model, env, H)
    @test length(obj.q) == H
    @test length(obj.u) == H
    @test length(obj.γ) == H
    @test length(obj.b) == H
end


text:
macro define(name, definition)
    quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end

alloc_x(x::Number) = typeof(x)(NaN)
alloc_f(x::Number) = real(typeof(x))(NaN)
alloc_d(x::Number) = typeof(x)(NaN)

alloc_x(x::AbstractArray) = eltype(x)(NaN) .* x
alloc_f(x::AbstractArray) = real(eltype(x))(NaN)
alloc_g(x::AbstractArray) = eltype(x)(NaN) .* x
alloc_h(x::AbstractArray) = eltype(x)(NaN) .* x*x'
alloc_j(x::AbstractArray, f::AbstractArray) = eltype(x)(NaN) .* vec(f) .* vec(x)'

function L2norm(x)
    local l2::eltype(x) = 0
    for xᵢ in x
        l2 += xᵢ^2
    end
    l2
end

function l2norm(x)
    sqrt(L2norm(x))
end

function maxnorm(x)
    local r² = zero(eltype(x))
    @inbounds for xᵢ in x
        r² = max(r², xᵢ^2)
    end
    sqrt(r²)
end

function outer!(O, x, y)
    @assert axes(O,1) == axes(x,1)
    @assert axes(O,2) == axes(y,1)
    @inbounds @simd for i in axes(O, 1)
        for j in axes(O, 2)
            O[i,j] = x[i] * y[j]
        end
    end

end