text:
module JLRasterization

using Colors, ColorVectorSpace, StaticArrays
using GeometryTypes, Interpolations
using ..Visualize: DepthBuffer, ColorBuffer, Area, Scene, AbstractWindow, Window, Canvas, JLCanvas, AbstractCanvas
import ..Visualize: rasterizer

using FieldTraits
using FieldTraits: ComposableLike, @composed, Partial
import FieldTraits: convertfor, default
import Interpolations: interpolate

@composed type JLWindow <: AbstractWindow
    <: Window
    Canvas
end

default(::Type{Canvas}, window::Partial{JLWindow}) = JLCanvas(window.val)

@inline function edge_function(a, b, c)
    (c[1] - a[1]) * (b[2] - a[2]) - (c[2] - a[2]) * (b[1] - a[1])
end
function src_alpha{T <: Colorant}(c::T)
    a = alpha(c)
    a == 0.0 ? zero(T) : c ./ a
end
one_minus_alpha{T <: Colorant}(c::T) = one(T) .- src_alpha(c)
blend(source, dest, src_func, dest_func) = clamp01(src_func(source) .+ dest_func(dest))
ColorTypes.alpha(x::StaticVector) = x[4]
function standard_transparency{T}(source, dest::T)
    (alpha(source) .* source) .+ ((one(eltype(T)) - alpha(source)) .* dest)
end

type FixedGeomView{GeomOut, VT}
    buffer::Vector{GeomOut}
    view::VT
    idx::Int
end

function FixedGeomView(T, max_primitives, primitive_in, primitive_out)
    buffer = Vector{Tuple{Point4f0, T}}(max_primitives)
    # TODO implement primitive_in and out correctly
    # this is for triangle_strip and 4 max_primitives
    if max_primitives != 4 || primitive_out == :triangle_strip
        error("Not implemented for max_primitives == $max_primitives and primitive_out == $primitive_out.")
    end
    geometry_view = if primitive_in == :triangle_strip
        view(buffer, [Face(1, 2, 3), Face(3, 2, 4)])
    else
        error("$primitive_in not supported. Only :triangle_strip supported right now")
    end
    FixedGeomView(buffer, geometry_view, 1)
end

function reset!(A::FixedGeomView)
    A.idx = 1
end
function Base.push!(A::FixedGeomView, element)
    if A.idx > length(A.buffer)
        error("Emit called more often than max_primitives. max_primitives: $(length(A.buffer))")
    end
    A.buffer[A.idx] = element
    A.idx += 1
    return
end

immutable JLRasterizer{Vertex, Args, FragN, VS, FS, GS, GV, EF}
    vertexshader::VS
    fragmentshader::FS

    geometryshader::GS
    geometry_view::GV
    emit::EF
end
function (::Type{JLRasterizer{Vertex, Args, FragN}}){Vertex, Args, FragN, VS, FS, GS, GV, EF}(
        vertexshader::VS,
        fragmentshader::FS,
        geometryshader::GS,
        geometry_view::GV,
        emit::EF
    )
    JLRasterizer{Vertex, Args, FragN, VS, FS, GS, GV, EF}(
        vertexshader,
        fragmentshader,
        geometryshader,
        geometry_view,
        emit
    )
end

function geometry_return_type(vertex_array, vertexshader, geometryshader, uniforms)
    typ = Any
    emit_t{T}(position, ::T) = (typ = T)
    face1 = first(vertex_array)
    vertex_stage = map(reverse(face1)) do f
        vertexshader(f, uniforms...)
    end
    geometryshader(emit_t, vertex_stage, uniforms...) # figure out typ
    typ
end

arglength{T <: Tuple}(::Type{T}) = length(T.parameters)
arglength{T <: AbstractArray}(::Type{T}) = 1
arglength{T}(::Type{T}) = nfields(T)

function rasterizer(
        window::AbstractWindow,
        vertexarray::AbstractArray,
        uniforms::Tuple,
        vertexshader::Function,
        fragmentshader::Function;
        geometryshader = nothing,
        max_primitives = 4,
        primitive_in = :points,
        primitive_out = :triangle_strip,
    )

    emit, geometry_view = nothing, nothing
    fragment_in_ndim = if geometryshader != nothing
        T = geometry_return_type(vertexarray, vertexshader, geometryshader, uniforms)
        geometry_view = FixedGeomView(T, max_primitives, primitive_in, primitive_out)
        emit = (position, fragment_args) -> push!(geometry_view, (position, fragment_args))
        arglength(T)
    else
        # when we don't have a geometry shader, vertex shader will feed fragment shader
        T = Base.Core.Inference.return_type(vertexshader, Tuple{eltype(vertexarray), map(typeof, uniforms)...})
        if T <: Tuple
            # TODO error handling
            arglength(T.parameters[2])
        else # if not a tuple, vertex shader doesn't pass any arguments to fragment shader
            0
        end
    end

    raster = JLRasterizer{eltype(vertexarray), typeof(uniforms), fragment_in_ndim}(
        vertexshader,
        fragmentshader,
        geometryshader,
        geometry_view,
        emit
    )
    raster, (vertexarray, uniforms)
end

Base.@pure Next{N}(::Val{N}) = Val{N - 1}()
@inline function interpolate{N, T}(bary, face::NTuple{N, T}, vn::Val{0}, aggregate)
    if T <: Tuple
        aggregate
    else
        T(aggregate...)
    end
end
@inline function interpolate{N}(bary, face, vn::Val{N}, aggregate = ())
    @inbounds begin
        res = (
            bary[1] * getfield(face[1], N) .+
            bary[2] * getfield(face[2], N) .+
            bary[3] * getfield(face[3], N)
        )
    end
    interpolate(bary, face, Next(vn), (res, aggregate...))
end

broadcastmin(a, b) = min.(a, b)
broadcastmax(a, b) = max.(a, b)

function clip2pixel_space(position, resolution)
    clipspace = position / position[4]
    p = clipspace[Vec(1, 2)]
    (((p + 1f0) / 2f0) .* (resolution - 1f0)) + 1f0
end

function (r::JLRasterizer{Vert, Args, FragN}){Vert, Args, FragN}(
        canvas, vertex_array::AbstractArray{Vert}, uniforms::Args
    )
    framebuffers = canvas[ColorBuffer]; depthbuffer = canvas[DepthBuffer]
    resolution = Vec2f0(size(framebuffers[1]))
    # hoisting out functions... Seems to help inference a bit. Or not?
    vshader = r.vertexshader
    gshader = r.geometryshader
    fshader = r.fragmentshader
    FragNVal = Val{FragN}()
    fragments_drawn = 0
    for face in vertex_array
        vertex_stage = map(reverse(face)) do f
            vshader(f, uniforms...)
        end
        geom_stage = if isa(r.geometryshader, Void)
            (vertex_stage,)
        else
            reset!(r.geometry_view)
            gshader(r.emit, vertex_stage, uniforms...)
            r.geometry_view.view
        end
        for geom_face in geom_stage
            fdepth = map(geom_face) do vert
                fv = first(vert)
                p = clip2pixel_space(fv, resolution)
                p[Vec(1, 2)], p[3]
            end
            f = map(first, fdepth)
            depths = map(last, fdepth)
            vertex_out = map(last, geom_face)
            # Bounding rectangle
            mini = max.(reduce(broadcastmin, f), 1f0)
            maxi = min.(reduce(broadcastmax, f), resolution)
            area = edge_function(f[1], f[2], f[3])
            for y = mini[2]:maxi[2]
                for x = mini[1]:maxi[1]
                    p = Vec(x, y)
                    w = Vec(
                        edge_function(f[2], f[3], p),
                        edge_function(f[3], f[1], p),
                        edge_function(f[1], f[2], p)
                    )
                    yi, xi = round(Int, y), round(Int, x)
                    if all(w-> w <= 0f0, w) && checkbounds(Bool, framebuffers[1], yi, xi)
                        bary = w / area
                        depth = bary[1] * depths[1] + bary[2] * depths[2] + bary[3] * depths[3]

                        if depth <= depthbuffer[yi, xi]
                            depthbuffer[yi, xi] = depth
                            fragment_in = interpolate(bary, vertex_out, FragNVal)
                            fragment_out = fshader(fragment_in, uniforms...)
                            for i = eachindex(fragment_out)
                                src_color = framebuffers[i][yi, xi]
                                dest_color = fragment_out[i]
                                fragments_drawn += 1
                                framebuffers[i][yi, xi] = standard_transparency(
                                    src_color,
                                    RGBA{Float32}(dest_color[1], dest_color[2], dest_color[3], dest_color[4])
                                )
                            end
                        end
                    end
                end
            end
        end
    end
    println("fragments drawn: ", fragments_drawn)
    return
end

end


text:
import AWSLambdaRuntimeInterfaceClient
const ric = AWSLambdaRuntimeInterfaceClient

import JSON

# Make handler
function handler(event, context)
    return "Hello World! $(JSON.json(event))"
end

# Pass handler to RIC and start
ric.start(handler)

ext: jl


text:
using OffsetArrays
using PencilArrays: Permutation, PermutedLinearIndices, PermutedCartesianIndices
using Test

@testset "Permuted indices" begin
    Aoff = OffsetArray(rand(3, 4, 5), -2, -1, -3)
    lin = LinearIndices(Aoff)
    cart = CartesianIndices(Aoff)

    perm = Permutation(3, 1, 2)
    plin = @inferred PermutedLinearIndices(lin, perm)
    pcart = @inferred PermutedCartesianIndices(cart, perm)

    for f in (size, axes)
        @test f(plin) == perm \ f(lin)
        @test f(pcart) == perm \ f(cart)
    end

    for n in LinearIndices(lin)
        @test lin[n] == plin[n]
        @test cart[n] == perm * pcart[n]
    end

    for I in CartesianIndices(lin)
        @test lin[I] == plin[perm \ I]
    end

    for J in CartesianIndices(plin)
        @test lin[perm * J] == plin[J]
    end

    # Iterate over permuted indices
    for n in plin
        @test n == lin[n] == plin[n]
        @test cart[n] == perm * pcart[n]
    end

    for J in pcart
        @test pcart[J] == J
    end

    for (I, J) in zip(cart, pcart)
        @test I == perm * J
    end
end


text:
function test_poly_mi(s)
    global solver = s

    test_dir = dirname(@__FILE__)

    @testset "Nonconvex MIPolynomial Models" begin end
end


text:
function get_random(seed,sequence,distribution)
    Random.seed!(seed)
    return last(rand(distribution,sequence))
end
# function get_random(seed,sequence,list)
#     Random.seed!(seed)
#     return last(rand(list,sequence))
# end
function get_random(seed,sequence)
    return get_random(seed,sequence,Uniform())
end

"""
Helper function that checks if x value is in vector V
"""
is_in(x::Int64,V) = x in V

"""
    Checks if a directory exists, otherwise create it.
"""
function check_create_dir!(dir_name)
    tdir = split(dir_name,"/")
    dirc = ""

    for tk=1:length(tdir)
        dirc *= tdir[tk] * "/"
        if !isdir(dirc)
            mkdir(dirc)
        end
    end
end

"""
Helper function to create a csv file template out of a given config object passed
"""
function create_csv_template(file_name,config)
    open(file_name * ".csv","w") do io
        writedlm(io,vcat(reshape(collect(keys(config)),1,length(config)),reshape(collect(values(config)),1,length(config))),";")
    end
end

"""
Find paths such as the network is covered according to passed coverage.
"""
function find_paths_by_seed(seed,g::SimpleGraph,coverage::Float64)
    Random.seed!(seed)
    cc = closeness_centrality(g)
    cci = sort([ (i,cc[i]) for i=1:length(cc) ],by=last,rev=true)

    pending = cci
    cp =[]

        while length(pending) >= (1- coverage) * nv(g)
            pending_i = [ first(p) for p in pending]

            #node with the most closeness_centrality
            s = first(first(pending))

            #max distance to any other pending node
            ds = gdistances(g,s)
            sds = sort([(i,ds[i]) for i=1:length(ds) if i in pending_i ],by=last,rev=true)
            d = first(first(sds))

            #shortest path between these two nodes
            sp = first(yen_k_shortest_paths(g,s,d).paths)
            push!(cp, sp)

            #remove nodes in the shortest path from pending list
            pending_i = collect(setdiff([ first(p) for p in pending ],Set(sp)))
            # println(pending_i)
            pending = [ p for p in pending  if first(p) in pending_i ]
        end
    return cp
end


text:
# script for Fig.8(b)(c), Fig.9, Fig.10(b)(c), Fig.11

using MultiscaleGraphSignalTransforms, LightGraphs, Plots; gr(dpi = 200)

## Build weighted graph
G, L, X = SunFlowerGraph(N = 400); N = nv(G)
𝛌, 𝚽 = eigen(Matrix(L))
sgn = (maximum(𝚽, dims = 1)[:] .> -minimum(𝚽, dims = 1)[:]) .* 2 .- 1
𝚽 = 𝚽 * Diagonal(sgn)
Q = incidence_matrix(G; oriented = true)
W = 1.0 * adjacency_matrix(G)
edge_weight = [e.weight for e in edges(G)]

## Build Dual Graph by DAG metric
distDAG = eigDAG_Distance(𝚽, Q, N; edge_weight = edge_weight)
Gstar_Sig = dualgraph(distDAG)
G_Sig = GraphSig(W, xy = X)
GP_dual = partition_tree_fiedler(Gstar_Sig; swapRegion = false)
GP_primal = pairclustering(𝚽, GP_dual)

@time VM_NGWP = vm_ngwp(𝚽, GP_dual) #54.986524 seconds (1.19 M allocations: 25.127 GiB, 2.36% gc time)
@time PC_NGWP = pc_ngwp(𝚽, GP_dual, GP_primal) #0.611176 seconds (225.41 k allocations: 844.488 MiB, 11.71% gc time)

#################### Fig. 8(b) barbara eye graph signal
using MAT
f = matread(joinpath(@__DIR__, "../datasets",
            "sunflower_barbara_voronoi.mat"))["f_eye_voronoi"]
G_Sig.f = reshape(f, (N, 1))
scatter_gplot(X; marker = f, ms = LinRange(4.0, 14.0, N), c = :greys);
plt = plot!(xlim = [-1.2,1.2], ylim = [-1.2,1.2], frame = :none)
# savefig(plt, joinpath(@__DIR__, "../paperfigs/SunFlower_barbara_feye.png"))

#################### Fig. 8(c) barbara eye relative l2 approximation error by various methods
DVEC = getall_expansioncoeffs(G_Sig, GP_dual, VM_NGWP, PC_NGWP, 𝚽)
approx_error_plot(DVEC);
plt = plot!(xguidefontsize = 16, yguidefontsize = 16, legendfontsize = 12)
# savefig(plt, joinpath(@__DIR__, "../paperfigs/SunFlower_barbara_feye_DAG_approx.png"))

#################### Fig. 9 barbara eye 16 most important VM-NGWP vectors (ignore the DC vector)
dmatrix_VM = ngwp_analysis(G_Sig, VM_NGWP)
dvec_vm_ngwp, BS_vm_ngwp = ngwp_bestbasis(dmatrix_VM, GP_dual)
important_idx = sortperm(dvec_vm_ngwp[:].^2; rev = true)
for i in 2:17
    dr, dc = BS_vm_ngwp.levlist[important_idx[i]]
    w = VM_NGWP[dr, dc, :]
    println("(j, k, l) = ", NGWP_jkl(GP_dual, dr, dc))
    scatter_gplot(X; marker = w, ms = LinRange(4.0, 14.0, N), c = :greys)
    plt = plot!(xlim = [-1.2, 1.2], ylim = [-1.2, 1.2], frame = :none,
                cbar = false, clims = (-0.15, 0.15))
    # savefig(plt, joinpath(@__DIR__,
    #     "../paperfigs/SunFlower_barbara_feye_DAG_VM_NGW_important_basis_vector$(lpad(i,2,"0")).png"))
end

#################### Fig. 10(b) barbara pants graph signal
f = matread(joinpath(@__DIR__, "../datasets",
            "sunflower_barbara_voronoi.mat"))["f_trouser_voronoi"]
scatter_gplot(X; marker = f, ms = LinRange(4.0, 14.0, N), c = :greys);
plt = plot!(xlim = [-1.2,1.2], ylim = [-1.2,1.2], frame = :none)
# savefig(plt, joinpath(@__DIR__, "../paperfigs/SunFlower_barbara_ftrouser.png"))

#################### Fig. 10(c) barbara eye relative l2 approximation error by various methods
G_Sig.f = reshape(f, (N, 1))
DVEC = getall_expansioncoeffs(G_Sig, GP_dual, VM_NGWP, PC_NGWP, 𝚽)
approx_error_plot(DVEC);
plt = plot!(xguidefontsize = 16, yguidefontsize = 16, legendfontsize = 12)
# savefig(plt, joinpath(@__DIR__, "../paperfigs/SunFlower_barbara_ftrouser_DAG_approx.png"))

#################### Fig. 11 barbara pants 16 most important VM-NGWP vectors (ignore the DC vector)
dmatrix_VM = ngwp_analysis(G_Sig, VM_NGWP)
dvec_vm_ngwp, BS_vm_ngwp = ngwp_bestbasis(dmatrix_VM, GP_dual)
important_idx = sortperm(dvec_vm_ngwp[:].^2; rev = true)
for i in 2:17
    dr, dc = BS_vm_ngwp.levlist[important_idx[i]]
    w = VM_NGWP[dr, dc, :]
    println("(j, k, l) = ", NGWP_jkl(GP_dual, dr, dc))
    scatter_gplot(X; marker = w, ms = LinRange(4.0, 14.0, N), c = :greys)
    plt = plot!(xlim = [-1.2, 1.2], ylim = [-1.2, 1.2], frame = :none,
                cbar = false, clims = (-0.15, 0.15))
    # savefig(plt, joinpath(@__DIR__,
    #     "../paperfigs/SunFlower_barbara_ftrouser_DAG_VM_NGW_important_basis_vector$(lpad(i,2,"0")).png"))
end


text:
# Custom function to allow for batch evaluation outside of julia
# or to exploit parallelized objective functions
# (included from "AbstractObjectiveInterface")
struct BatchObjectiveFunction{F<:Function} <: Function
    function_handle :: F
end

# evaluation of a BatchObjectiveFunction
function (objf::BatchObjectiveFunction)(x :: Vec)
    return ensure_vec(objf.function_handle(x))
end

# overload broadcasting for BatchObjectiveFunction's
# that are assumed to handle arrays themselves
function Broadcast.broadcasted( objf::BatchObjectiveFunction, X :: VecVec)
    ensure_vec.( objf.function_handle( X ) )
end

function _new_batch( func1 :: F, func2 :: T ) where{ F <:Function ,T <: Function }
    return BatchObjectiveFunction(
        function( x :: Union{ VecVec, Vec } )
            if isa( x, Vec )
                [ func1(x); func2(x) ]
            else
                f1 = func1.(x)
                f2 = func2.(x)
                [ [f1[i];f2[i] ] for i = eachindex(x) ]
            end
        end
    )
end

combine( func1 :: Function, func2 :: BatchObjectiveFunction ) = _new_batch( func1, func2)
combine( func1 :: BatchObjectiveFunction, func2 :: Function) = _new_batch( func1, func2 )
combine( func1 :: BatchObjectiveFunction, func2 :: BatchObjectiveFunction ) = _new_batch( func1, func2 )

"Get a new function function handle stacking the output of `func1` and `func2`."
function combine(func1 :: F, func2 :: T ) where{ F <:Function ,T <: Function }
    return function( x :: Vec )
        return vcat( func1(x), func2(x) )
    end
end


text:
module IntroToJulia

function find_easter()

println("""
Welcome to the IntroToJulia Easter Egg hunt!

If you found this, the you have not successfully checked out
the Easter branch of the repository. Please go through the
materials again to find out how to change the branch of
of a repository, or ask for help!
""")
end

export find_easter

end

max_issues_repo_name: BenjaminBorn/IntroToJulia

max_forks_repo_name: BenjaminBorn/IntroToJulia


text:
using Bennu
using WriteVTK

grid = brickgrid(LobattoCell(5, 4), (-1:1//8:1, -1:1//8:1);
                 ordering=HilbertOrdering())

vtk_grid("grid", grid) do vtk
    vtk["CellNumber"] = 1:length(grid)
end


text:
using RNGPool
using RData
using CoupledConditionalSMC
import Statistics: mean, var

dataJLS = RData.load("demo/ar1data.RData")
ys = vec(dataJLS["observations"])

include("../demo/lgModel.jl")

# they start at time 0 with observations from time 1,2,...
# so we have an initial distribution of N(0,1+0.9^2)
theta = LGTheta(0.9, 1.0, 1.0, 1.0, 0.0, 1.81)
model = makeLGModel(theta, ys)
lM = LinearGaussian.makelM(theta)
ko = kalman(theta, ys)

setRNGs(12345)

function runDemoBS(n::Int64, N::Int64, m::Int64, maxcost::Int64)
    maxiter = Int(round(maxcost/N))
    results = CoupledConditionalSMC.couplingTimes(model, lM, N, n, m,
                                                 :BS, true, true, maxiter)
    results[results .== 0] .= maxiter
    results.*N
end

function runDemoAS(n::Int64, N::Float64, m::Int64, maxcost::Int64)
    N_ = Int(ceil(N*n))
    maxiter = Int(round(maxcost/N_))
    results = CoupledConditionalSMC.couplingTimes(model, lM, N_, n, m,
                                                  :AS, true, true, maxiter)
    results[results .== 0] .= maxiter
    results.*N_
end

function runDemoAT(n::Int64, N::Float64, m::Int64, maxcost::Int64)
  # Scale number of particles proportional to length
  N_ = Int(ceil(N*n))
  maxiter = Int(round(maxcost/N_))
  results = CoupledConditionalSMC.couplingTimes(model, N_, n, m,
                                                true, true, maxiter)
  results[results .== 0] .= maxiter
  results.*N_
end

# Number of replications
Nrepl = 10_000
# Time series lengths
ts = collect(100:50:300)
# Number of particles
Ns = 2 .^ (2:6)

BS_cost = zeros(length(ts), length(Ns), Nrepl)
AS_cost = deepcopy(BS_cost)

#Threads.@threads
for i in length(ts):-1:1
  t = ts[i]
  for j in 1:length(Ns)
    N = Ns[j]
    # BS estimator with given number of particles:
    BS_cost[i,j,:] = runDemoBS(t, N, Nrepl, 50000)
    # AT estimator with t*N/10 particles:
    #AT_cost[i,j,:] = runDemoAT(t, N/2.0, Nrepl, 50000)
    AS_cost[i,j,:] = runDemoAS(t, N/5.0, Nrepl, 50000)
  end
end

using JLD2
@save "paper/output/lgDemoCompareOptimalCost_results.jld2" BS_cost AS_cost ts Nrepl Ns


text:
#Load a (trained) neural network from an HDF5 file
function nnLoad(fileName::String)
  #read from file
  f = h5open(fileName);
  nNeurons = read(f,"nNeurons");
  w        = read(f,"w");
  b        = read(f,"b");
  close(f);

  #build the neural net
  nn = neuralNet(nNeurons);

  #assign the saved weights and biases
  nn.w = w;
  nn.b = b;

  return nn;
end

max_stars_repo_name: krometis/SimpleNeuralNet

max_issues_repo_name: krometis/SimpleNeuralNet

max_forks_repo_name: krometis/SimpleNeuralNet


text:
"""
    AnalyticalNonlinearShrinkage

Analytical nonlinear shrinkage estimator. See docs for
`analytical_nonlinear_shrinkage` for details.
"""
struct AnalyticalNonlinearShrinkage{TEigen<:Union{Eigen,Nothing}} <: CovarianceEstimator
    decomp::TEigen
    corrected::Bool
    function AnalyticalNonlinearShrinkage(;corrected=false)
        new{Nothing}(nothing, corrected)
    end
    function AnalyticalNonlinearShrinkage(decomp::Eigen; corrected=false)
        new{Eigen}(decomp, corrected)
    end
end

const SQRT5  = sqrt(5.0)
const INVPI  = 1.0 / π
const EPAN_1 = 3.0/(4.0 * SQRT5)
const EPAN_2 = EPAN_1 * INVPI
const EPAN_3 = 0.3 * INVPI

"""
    epanechnikov(x)

Return the Epanechnikov kernel evaluated at `x`.
"""
epanechnikov(x::Real) = EPAN_1 * max(0.0, 1.0 - x^2/5.0)

"""
    epnanechnikov_HT(x)

Return the Hilbert Transform of the Epanechnikov kernel evaluated at `x`
if `|x|≂̸√5`.
"""
function epanechnikov_HT1(x::Real)
    -EPAN_3 * x + EPAN_2 * (1.0 - x^2/5.0) * log(abs((SQRT5 - x)/(SQRT5 + x)))
end

"""
    epnanechnikov_HT2(x)
Return the Hilbert Transform of the Epanechnikov kernel evaluated at `x`
if `|x|=√5`.
"""
epanechnikov_HT2(x::Real) = -EPAN_3*x

"""
    analytical_nonlinear_shrinkage(S, n, p; decomp)

Internal implementation of the analytical nonlinear shrinkage. The
implementation is inspired from the Matlab code given in section C of
Olivier Ledoit and Michael Wolf's paper "Analytical Nonlinear
Shrinkage of Large-Dimensional Covariance Matrices". (Nov 2018)
http://www.econ.uzh.ch/static/wp/econwp264.pdf
"""
function analytical_nonlinear_shrinkage(S::AbstractMatrix{<:Real},
                                        n::Int, p::Int, est_mean::Bool;
                                        decomp::Union{Nothing,Eigen}=nothing)

    η = ifelse(p < n, n, n - Int(est_mean)) # effective sample size
    # sample eigenvalues sorted in ascending order and eigenvectors
    F    = isa(decomp, Nothing) ? eigen(S) : decomp # O(p^3)
    perm = sortperm(F.values)
    sample_perm = @view perm[max(1, (p - η) + 1):p]
    λ    = @view F.values[sample_perm]
    U    = F.vectors[:, perm]

    # dominant cost forming of S or eigen(S) --> O(max{np^2, p^3})

    # compute analytical nonlinear shrinkage kernel formula
    L = repeat(λ, outer=(1, min(p, η)))

    # Equation (4.9)
    h = η^(-1//3)
    H = h * L'
    x = (L .- L') ./ H

    # additional useful definitions
    γ  = p/η
    πλ = π * λ

    # Equation (4.7) in http://www.econ.uzh.ch/static/wp/econwp264.pdf
    f̃ = vec(mean(epanechnikov.(x) ./ H, dims=2))

    # Equation (4.8)
    Hf̃_tmp = epanechnikov_HT1.(x)
    # if by any chance there are x such that |x| ≈ √5 ...
    mask = (@. abs(x) ≈ SQRT5)
    any(mask) && (Hf̃_tmp[mask] = epanechnikov_HT2.(@view x[mask]))
    Hf̃ = vec(mean(Hf̃_tmp ./ H, dims=2))

    # dominant cost up to here: elementwise ops on x --> O(max{p^2, η^2})

    if p <= η
        # Equation (4.3)
        πγλ = γ * πλ
        denom = @. (πγλ * f̃)^2 + (1.0 - γ - πγλ * Hf̃)^2
        d̃ = λ ./ denom
    else
        # Equation (C.8)
        hs5  = h * SQRT5
        Hf̃0  = (0.3/h^2 + 0.75/hs5 * (1.0 - 0.2/h^2) * log((1+hs5)/(1-hs5)))
        Hf̃0 *= mean(1.0 ./ πλ)
        # Equation (C.5)
        d̃0 = INVPI / ((γ - 1.0) * Hf̃0)
        # Eq. (C.4)
        d̃1 = @. 1.0 / (π * πλ * (f̃^2 + Hf̃^2))
        d̃  = [fill(d̃0, (p - η, 1)); d̃1]
    end

    # Equation (4.4)
    return Symmetric(U*(d̃ .* U'))
end

"""
    cov(ans::AnalyticalNonlinearShrinkage, X; dims=1, mean=nothing)

Nonlinear covariance estimator derived from the sample covariance estimator `S`
and its eigenvalue decomposition (which can be given through `decomp`).
See Ledoit and Wolf's paper http://www.econ.uzh.ch/static/wp/econwp264.pdf
The keyword `mean` can be `nothing` (centering via estimated mean),
zero (no centering) or a provided vector. In the first case, a rank-1
modification is applied and therefore the effective sample size is decreased
by one (see `analytical_nonlinear_shrinkage`). In the latter two case the mean
cannot have been estimated on the data (otherwise the effective sample size
will be 1 larger than it should be resulting in numerical instabilities).
If you are unsure, use either `nothing` or provide an explicit (non-estimated)
vector (possibly a zero vector) and avoid the use of `mean=0`.

* Time complexity (including formation of `S`)
    - (p<n): O(np^2 + n^2) with moderate constant
    - (p>n): O(p^3) with low constant (dominated by eigendecomposition of S)
"""
function cov(ans::AnalyticalNonlinearShrinkage, X::AbstractMatrix{<:Real};
             dims::Int=1, mean=nothing)

    @assert dims ∈ [1, 2] "Argument dims can only be 1 or 2 (given: $dims)"

    szx    = size(X)
    (n, p) = ifelse(dims==1, szx, reverse(szx))

    # explained in the paper there must be at least 12 samples
    (n < 12) && throw(ArgumentError("The number of samples `n` must be at " *
                                    "least 12 (given: $n)."))

    S  = cov(SimpleCovariance(corrected=ans.corrected), X; dims=dims, mean=mean)
    return analytical_nonlinear_shrinkage(S, n, p, mean === nothing;
                decomp=ans.decomp)
end


text:
module qTab

function xcodegen(labels)

  ob = "(";
  cb = ")";
  xstr = "" * ob;
  for (i, j) in enumerate(labels)
    xstr = xstr * j * " = datar[$(i)], ";
    end;
  xstr = xstr[1:end-2] * cb;
  return xstr;
  end;

function qfunc(header)

  dataexec = xcodegen(header);
  myf = eval(Meta.parse(" function (datar) \n" * dataexec * "\n end;"));
  return myf;
  end;

function quickTable(data, xfun)

  ret = [ xfun(data[i, :]) for i in 1:size(data,1)];
  return ret;
  end;

function quickTable(data)

  x = qfunc(data[1, :]);
  return Base.invokelatest(quickTable, data, x);
  end;

export quickTable, qfunc;

end # module


text:
# Part 1
program = readlines("day14_input.txt")
getval(line) = last(split(line, " = "))
getmemloc(line) = parse(Int64, first(split(line, " = "))[5:end-1])
applymask(mask, val) = parse(Int64, join([m == 'X' ? v : m for (m, v) in zip(mask, bitstring(val)[end-35:end])]), base=2)

function run(prog)
    mask = repeat("0", 36)
    mem = Dict()
    for line in prog
        if startswith(line, "mask")
            mask = getval(line)
        else
            mem[getmemloc(line)] = applymask(mask, parse(Int64, getval(line)))
        end
    end
    mem
end
println("Part 1: sum = $(sum(values(run(program))))")

# Part 2
function masklocs(mask, loc)
    idx = 1
    locformat = ""
    for (m, v) in zip(mask, bitstring(loc)[end-35:end])
        if m == 'X'
            locformat *= "[$idx]"
            idx += 1
        elseif m == '1'
            locformat *= '1'
        else
            locformat *= v
        end
    end

    if idx == 1 return [parse(Int64, locformat, base=2)] end
    xcount = idx - 1

    locs = []
    for i in 0:2^xcount-1
        replaced = locformat
        for (j, x) in enumerate(collect(bitstring(i)[end-(xcount-1):end]))
            replaced = replace(replaced, "[$j]" => "$x")
        end
        push!(locs, replaced)
    end
    locs
end

function run2(prog)
    mask = repeat("0", 36)
    mem = Dict()
    for line in prog
        if startswith(line, "mask")
            mask = getval(line)
        else
            loc = getmemloc(line)
            for ml in masklocs(mask, loc)
                mem[ml] = parse(Int64, getval(line))
            end
        end
    end
    mem
end
println("Part 2: sum = $(sum(values(run2(program))))")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
# v0.6.0

function printlyrics()
    const gifts = split("""
    A partridge in a pear tree.
    Two turtle doves
    Three french hens
    Four calling birds
    Five golden rings
    Six geese a-laying
    Seven swans a-swimming
    Eight maids a-milking
    Nine ladies dancing
    Ten lords a-leaping
    Eleven pipers piping
    Twelve drummers drumming
    """, '\n')
    const days = split("""
    first second third fourth fifth
    sixth seventh eighth ninth tenth
    eleventh twelfth""")
    for (n, day) in enumerate(days)
        g = gifts[n:-1:1]
        print("\nOn the $day day of Christmas\nMy true love gave to me:\n")
        if n == 1
            print(join(g[1:end], '\n'), '\n')
        else
            print(join(g[1:end-1], '\n'), " and\n", g[end], '\n')
        end
    end
end

printlyrics()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
abstract type AbstractThreadPool end

@deprecate pforeach(pool, fn::Function, itr) tforeach(pool, fn::Function, itr)
@deprecate pforeach(fn::Function, pool, itr) tforeach(fn::Function, pool, itr)
@deprecate pmap(pool, fn::Function, itr) tmap(pool, fn::Function, itr)
@deprecate pmap(fn::Function, pool, itr) tmap(fn::Function, pool, itr)

_detect_type(fn, itr) = eltype(map(fn, empty(itr)))
_detect_type(fn, itrs::Tuple) = eltype(map(fn, [empty(x) for x in itrs]...))

"""
    tforeach(pool, fn::Function, itr)
    tforeach(fn::Function, pool, itr)

Mimics `Base.foreach`, but launches the function evaluations onto the provided
pool to assign the tasks.

# Example
```
julia> pool = pwith(ThreadPools.LoggedQueuePool(1,2)) do pool
         tforeach(pool,  x -> println((x,Threads.threadid())), 1:8)
       end;
(2, 2)
(1, 1)
(3, 2)
(5, 2)
(4, 1)
(6, 2)
(7, 1)
(8, 2)

julia> plot(pool)
```
"""
function tforeach(pool, fn::Function, itr::AbstractVector)
    tmap(pool, fn, itr)
    nothing
end

tforeach(fn::Function, pool, itr) = tforeach(pool, fn, itr::AbstractVector)
tforeach(pool, fn::Function, itrs...) = tforeach(pool, (x) -> fn(x...), zip(itrs...))
tforeach(fn::Function, pool, itrs...) = tforeach(pool, (x) -> fn(x...), zip(itrs...))

# function tmap(pool::AbstractThreadPool, fn::Function, itr::AbstractVector)
#     error("Not Implemented")
# end

"""
    tmap(pool, fn::Function, itr)
    tmap(fn::Function, pool, itr)

Mimics `Base.map`, but launches the function evaluations onto the provided
pool to assign the tasks.

# Example
```
julia> pool = pwith(ThreadPools.LoggedQueuePool(1,2)) do pool
         tmap(pool, 1:8) do x
           println((x,Threads.threadid()))
         end
       end;
(2, 2)
(1, 1)
(3, 2)
(4, 1)
(5, 2)
(6, 1)
(7, 2)
(8, 1)

julia> plot(pool)
```
"""
tmap(fn::Function, pool, itr) = tmap(pool, fn, itr)
tmap(pool, fn::Function, itrs...) = tmap(pool, (x) -> fn(x...), zip(itrs...))
tmap(fn::Function, pool, itrs...) = tmap(pool, (x) -> fn(x...), zip(itrs...))

"""
    pwith(fn::Function, pool) -> pool

Apply the functon `fn` to the provided pool and close the pool.  Returns the
closed pool for any desired analysis or plotting.

# Example
```
julia> pwith(ThreadPools.QueuePool(1,2)) do pool
         tforeach(pool, x -> println((x,Threads.threadid())), 1:8)
       end;
(2, 2)
(1, 1)
(3, 2)
(4, 1)
(5, 2)
(6, 1)
(7, 2)
(8, 1)
```
Note in the above example, only two threads were used, as set by the
`QueuePool` setting.
"""
function pwith(fn::Function, pool)
    fn(pool)
    close(pool)
    pool
end

"""
    @pthreads pool

Mimic the `Base.Threads.@threads` macro, but uses the provided pool to
assign the tasks.

# Example
```
julia> pwith(ThreadPools.QueuePool(1,2)) do pool
         @pthreads pool for x in 1:8
           println((x,Threads.threadid()))
         end
       end;
(2, 2)
(3, 2)
(1, 1)
(4, 2)
(5, 1)
(6, 2)
(8, 2)
(7, 1)
```
"""
macro pthreads(pool, args...)
    na = length(args)
    if na != 1
        throw(ArgumentError("wrong number of arguments in @pthreads"))
    end
    ex = args[1]
    if !isa(ex, Expr)
        throw(ArgumentError("need an expression argument to @pthreads"))
    end
    if ex.head === :for
        if ex.args[1] isa Expr && ex.args[1].head === :(=)
            index = ex.args[1].args[1]
            range = ex.args[1].args[2]
            body = ex.args[2]
            return quote
                tforeach($(esc(pool)), $(esc(range))) do $(esc(index))
                    $(esc(body))
                end
            end
        else
            throw(ArgumentError("nested outer loops are not currently supported by @pthreads"))
        end
    else
        throw(ArgumentError("unrecognized argument to @pthreads"))
    end
end

function Base.finalize(pool::AbstractThreadPool)
    close(pool)
end

"""
    Base.close(pool::AbstractThreadPool)

Closes the pool, shuts down any handlers and finalizes any logging activities.
"""
function Base.close(pool::AbstractThreadPool)
    nothing
end


text:
using LinearAlgebraicRepresentation, ViewerGL, SparseArrays
Lar = LinearAlgebraicRepresentation; GL = ViewerGL

# using Revise
using OhMyREPL

# 3D Boolean example generation
#-------------------------------------------------------------------------------
n,m = 1,1
V,(VV,EV,FV) = Lar.cuboidGrid([n,m],true)
square = V,FV,EV
q = sqrt(2)
V, (VV,EV,FV) = Lar.cuboid([-q,-q], true, [q,q])

twosquares = Lar.Struct([ #(V,FV,EV),
    Lar.Struct([ Lar.t(-sqrt(2)/4, -sqrt(2)/2 ), Lar.r(pi/4), square ]),
    Lar.Struct([ Lar.t( sqrt(2)/4, -sqrt(2)/2 ), Lar.r(pi/4), square ])
    ])

V,FV,EV = Lar.struct2lar(twosquares)
GL.VIEW([ GL.GLGrid(V,FV), GL.GLFrame2 ]);
W,EW = Lar.fragmentlines((V,EV))
V,FVs,EVs = Lar.arrange2D(W,EW)
GL.VIEW(GL.GLExplode(V,FVs[1:3],1.,1.,1.,99,1));

circle = Lar.circle(1.06066/2)()

n,m = 1,1
V,(VV,EV,FV) = Lar.cuboidGrid([n,m],true)
square = V,EV
q = sqrt(2)
assembly = Lar.Struct([ #(V,FV,EV),
    # Lar.Struct([ Lar.t(-q/4, 0 ), Lar.r(pi/4), circle ]),
    Lar.Struct([ Lar.t( -q/4, -q/2 ), Lar.r(pi/4), square ])
    Lar.Struct([ Lar.t( q/4, -q/2 ), Lar.r(pi/4), square ])
])
V,EV = Lar.struct2lar(assembly)
GL.VIEW([ GL.GLGrid(V,EV, GL.COLORS[1],1), GL.GLFrame2 ]);

# boolmatrix = Lar.booleanops2d((V,EV))
W, copEV, copFE, boolmatrix = Lar.bool2d(assembly)
@show Matrix(boolmatrix)
@show Matrix(copFE)
A = boolmatrix[:,1]
B = boolmatrix[:,2]
AorB = A .| B
AandB = A .& B
AxorB = AorB .& (.! AandB)
expr = .&(AorB,(.!AandB)) # xor in prefix polac form

# TODO new method for "bool2d()" in src/Bool2d.jl (see template)
# W, copEV, copFE, boolmat, result = Lar.bool2d( expr, assembly )

union = Matrix(copFE)' * Int.(AorB)
intersection = Matrix(copFE)' * Int.(AandB)
xor = Matrix(copFE)' * AxorB

V = convert(Lar.Points,W')
EV = Lar.cop2lar(copEV)
EVor = [ev for (k,ev) in enumerate(EV) if abs(union[k])==1 ]
EVand = [ev for (k,ev) in enumerate(EV) if abs(intersection[k])==1 ]
EVxor = [ev for (k,ev) in enumerate(EV) if abs(xor[k])==1 ]
GL.VIEW([ GL.GLGrid(V,EVor, GL.COLORS[1],1), GL.GLFrame2 ]);
GL.VIEW([ GL.GLGrid(V,EVand, GL.COLORS[1],1), GL.GLFrame2 ]);
GL.VIEW([ GL.GLGrid(V,EVxor, GL.COLORS[1],1), GL.GLFrame2 ]);

VV = [[k] for k in 1:size(V,2)];
model = (V,[VV, EV])
GL.VIEW( GL.numbering(.5)( model,GL.COLORS[1],0.1 ) );


text:
function test_linop()
  (nrow, ncol) = (10, 6)
  ϵ = eps(Float64)
  rtol = sqrt(ϵ)
  A1 = simple_matrix(ComplexF64, nrow, ncol)

  @testset ExtendedTestSet "Basic operations" begin
    for op in (
      LinearOperator(A1),
      LinearOperator(A1')',
      transpose(LinearOperator(transpose(A1))),
      conj(LinearOperator(conj(A1))),
      # PreallocatedLinearOperator(A1),
    )
      show(op)

      @testset "Data type" begin
        @test eltype(op) == eltype(A1)
        @test !isreal(op)
      end

      @testset "Size" begin
        @test(size(op) == (nrow, ncol))
        @test(shape(op) == (nrow, ncol))
        @test(size(op, 1) == nrow)
        @test(size(op, 2) == ncol)
        @test_throws LinearOperatorException size(op, 3)
        @test_throws LinearOperatorException op * ones(ncol + 1)
      end

      @testset "Boolean operators" begin
        @test(symmetric(op) == false)
        @test(hermitian(op) == false)
      end

      @testset "Full" begin
        @test(norm(A1 - Matrix(op)) <= ϵ * norm(A1))
      end

      @testset "Unary +." begin
        @test(norm(A1 - Matrix(+op)) <= ϵ * norm(A1))
      end
    end

    @testset "LinearOperator(Matrix)" begin
      A2 = simple_sparse_matrix(ComplexF64, nrow, ncol)
      for A in (A1, A2)
        op = LinearOperator(A)
        @test(op.nrow == nrow)
        @test(op.ncol == ncol)

        @test(norm(A - Matrix(op)) <= rtol * norm(A))
        @test(norm(transpose(A) - Matrix(transpose(op))) <= rtol * norm(A))
        @test(norm(A' - Matrix(op')) <= rtol * norm(A))

        v = simple_vector(ComplexF64, ncol)
        @test(norm(A * v - op * v) <= rtol * norm(v))

        u = simple_vector(ComplexF64, nrow)
        @test(norm(transpose(A) * u - transpose(op) * u) <= rtol * norm(u))
        @test(norm(A' * u - op' * u) <= rtol * norm(u))
      end
    end

    # @testset "Preallocated LinearOperator(Matrix)" begin
    #   A2 = simple_matrix(Float64, nrow, ncol)
    #   for A in (A1, A2)
    #     T = eltype(A)
    #     Mv = zeros(T, nrow)
    #     Mtu = zeros(T, ncol)
    #     Maw = zeros(T, ncol)
    #     op = PreallocatedLinearOperator(Mv, Mtu, Maw, A)
    #     v = simple_vector(T, ncol)
    #     u = simple_vector(T, nrow)
    #     @test norm(op * v - A * v) <= rtol * norm(A)
    #     @test norm(transpose(op) * u - transpose(A) * u) <= rtol * norm(A)
    #     @test norm(op' * u - A' * u) <= rtol * norm(A)

    #     al = @allocated op * v
    #     @test al == 0

    #     op = PreallocatedLinearOperator(A)
    #     v = simple_vector(T, ncol)
    #     u = simple_vector(T, nrow)
    #     @test norm(op * v - A * v) <= rtol * norm(A)
    #     @test norm(transpose(op) * u - transpose(A) * u) <= rtol * norm(A)
    #     @test norm(op' * u - A' * u) <= rtol * norm(A)
    #   end

    #   for B in (simple_matrix(Float64, nrow, nrow), simple_sparse_matrix(Float64, nrow, nrow))
    #     for A in (SymTridiagonal(Symmetric(B)), Symmetric(B))
    #       v = simple_vector(Float64, nrow)

    #       Mv = zeros(nrow)
    #       op = PreallocatedLinearOperator(Mv, A)
    #       @test norm(op * v - A * v) <= rtol * norm(A)
    #       @test norm(transpose(op) * v - transpose(A) * v) <= rtol * norm(A)
    #       @test norm(op' * v - A' * v) <= rtol * norm(A)

    #       op = PreallocatedLinearOperator(A)
    #       @test norm(op * v - A * v) <= rtol * norm(A)
    #       @test norm(transpose(op) * v - transpose(A) * v) <= rtol * norm(A)
    #       @test norm(op' * v - A' * v) <= rtol * norm(A)
    #     end
    #   end

    #   v = simple_vector(Float64, nrow)

    #   A = simple_matrix(ComplexF64, nrow, nrow)
    #   # Symmetric and Hermitian - actually a Real matrix, but with Complex type.
    #   # This tests a specific condition located in PreallocatedLinearOperators.jl
    #   # when eltype(A) is not Real but A is symmetric and hermitian.
    #   A = A + A'
    #   A = A + transpose(A)
    #   op = PreallocatedLinearOperator(A, symmetric = true, hermitian = true)
    #   @test norm(op * v - A * v) <= rtol * norm(A)
    #   @test norm(transpose(op) * v - transpose(A) * v) <= rtol * norm(A)
    #   @test norm(op' * v - A' * v) <= rtol * norm(A)

    #   A = simple_matrix(ComplexF64, nrow, nrow)
    #   A = A + A' # Hermitian
    #   op = PreallocatedLinearOperator(A, symmetric = false, hermitian = true)
    #   @test norm(op * v - A * v) <= rtol * norm(A)
    #   @test norm(transpose(op) * v - transpose(A) * v) <= rtol * norm(A)
    #   @test norm(op' * v - A' * v) <= rtol * norm(A)

    #   A = simple_matrix(ComplexF64, nrow, nrow)
    #   A = A + transpose(A) # Symmetric
    #   op = PreallocatedLinearOperator(A, symmetric = true, hermitian = false)
    #   @test norm(op * v - A * v) <= rtol * norm(A)
    #   @test norm(transpose(op) * v - transpose(A) * v) <= rtol * norm(A)
    #   @test norm(op' * v - A' * v) <= rtol * norm(A)

    #   A = simple_matrix(ComplexF64, nrow, nrow)
    #   A = Hermitian(A)
    #   op = PreallocatedLinearOperator(A)
    #   @test norm(op * v - A * v) <= rtol * norm(A)
    #   @test norm(transpose(op) * v - transpose(A) * v) <= rtol * norm(A)
    #   @test norm(op' * v - A' * v) <= rtol * norm(A)
    # end

    @testset "Basic arithmetic operations" begin
      B1 = simple_matrix(ComplexF64, nrow, ncol)

      for q in (+, -)
        C = q(A1, B1)
        opC = q(LinearOperator(A1), LinearOperator(B1))
        v = simple_vector(ComplexF64, ncol)
        @test(norm(opC * v - C * v) <= rtol * norm(v))
        u = simple_vector(ComplexF64, nrow)
        @test(norm(transpose(opC) * u - transpose(C) * u) <= rtol * norm(u))
        @test(norm(opC' * u - C' * u) <= rtol * norm(u))

        opC = q(A1, LinearOperator(B1))
        @test(norm(opC * v - C * v) <= rtol * norm(v))
        @test(norm(transpose(opC) * u - transpose(C) * u) <= rtol * norm(u))
        @test(norm(opC' * u - C' * u) <= rtol * norm(u))

        opC = q(LinearOperator(A1), B1)
        @test(norm(opC * v - C * v) <= rtol * norm(v))
        @test(norm(transpose(opC) * u - transpose(C) * u) <= rtol * norm(u))
        @test(norm(opC' * u - C' * u) <= rtol * norm(u))
      end
    end

    B2 = simple_matrix(ComplexF64, ncol, ncol + 1)
    @testset "Operator ± scalar" begin
      opC = LinearOperator(A1) + 2.12345
      @test(norm(A1 .+ 2.12345 - Matrix(opC)) <= rtol * norm(A1 .+ 2.12345))

      opC = 2.12345 + LinearOperator(A1)
      @test(norm(A1 .+ 2.12345 - Matrix(opC)) <= rtol * norm(A1 .+ 2.12345))

      opC = LinearOperator(A1) - 2.12345
      @test(norm((A1 .- 2.12345) - Matrix(opC)) <= rtol * norm(A1 .- 2.12345))

      opC = 2.12345 - LinearOperator(A1)
      @test(norm((2.12345 .- A1) - Matrix(opC)) <= rtol * norm(2.12345 .- A1))

      C = A1 * B2
      opC = LinearOperator(A1) * LinearOperator(B2)
      v = simple_vector(ComplexF64, ncol + 1)
      @test(norm(opC * v - C * v) <= rtol * norm(v))
      u = simple_vector(ComplexF64, nrow)
      @test(norm(transpose(opC) * u - transpose(C) * u) <= rtol * norm(u))
      @test(norm(opC' * u - C' * u) <= rtol * norm(u))

      @test_throws LinearOperatorException LinearOperator(A1) + LinearOperator(B2)
      @test_throws LinearOperatorException LinearOperator(B2) * LinearOperator(A1)
    end

    A1B2 = A1 * B2
    @testset "Matrix × operator" begin
      opC = A1 * LinearOperator(B2)
      @test(norm(A1B2 - Matrix(opC)) <= rtol * norm(A1B2))
    end

    @testset "Operator × matrix" begin
      opC = LinearOperator(A1) * B2
      @test(norm(A1B2 - Matrix(opC)) <= rtol * norm(A1B2))
    end

    AA1 = 2.12345 * A1
    @testset "Scalar × operator" begin
      opC = 2.12345 * LinearOperator(A1)
      @test(norm(AA1 - Matrix(opC)) <= rtol * norm(AA1))
    end

    @testset "Operator × scalar" begin
      opC = LinearOperator(A1) * 2.12345
      @test(norm(AA1 - Matrix(opC)) <= rtol * norm(AA1))
    end
  end

  @testset ExtendedTestSet "Basic operators" begin
    @testset "Identity" begin
      opI = opEye(nrow)
      v = simple_vector(ComplexF64, nrow)
      @test(abs(norm(opI * v - v)) <= ϵ * norm(v))
      @test(abs(norm(transpose(opI) * v - v)) <= ϵ * norm(v))
      @test(abs(norm(opI' * v - v)) <= ϵ * norm(v))
      @test(norm(Matrix(opI) - Matrix(1.0I, nrow, nrow)) <= ϵ * norm(Matrix(1.0I, nrow, nrow)))

      w = opI * v
      w[1] = -1.0
      @test v[1] != w[1]

      opI = opEye(nrow, ncol)
      v = simple_vector(ComplexF64, ncol)
      v0 = [v; zeros(nrow - ncol)]
      vu = [v; ones(nrow - ncol)]
      @test(abs(norm(opI * v - v0)) <= ϵ * norm(v))
      @test(abs(norm(transpose(opI) * vu - v)) <= ϵ * norm(v))
      @test(abs(norm(opI' * vu - v)) <= ϵ * norm(v))
      @test(norm(Matrix(opI) - Matrix(1.0I, nrow, ncol)) <= ϵ * norm(Matrix(1.0I, nrow, ncol)))

      opI = opEye(ncol, nrow)
      @test(abs(norm(opI * vu - v)) <= ϵ * norm(v))
      @test(abs(norm(transpose(opI) * v - v0)) <= ϵ * norm(v))
      @test(abs(norm(opI' * v - v0)) <= ϵ * norm(v))
      @test(norm(Matrix(opI) - Matrix(1.0I, ncol, nrow)) <= ϵ * norm(Matrix(1.0I, ncol, nrow)))
    end

    @testset "Identity (non-convertible to matrix)" begin
      op = opEye()

      v = simple_vector(Float64, 5)
      w = op * v
      @test w === v
      w = v * op

      @test w === v
      A2 = op * A1
      @test A2 === A1
      A2 = A1 * op
      @test A2 === A1

      T1 = LinearOperator(A1)
      T2 = op * T1
      @test T2 === T1
      T2 = T1 * op
      @test T2 === T1

      op2 = opEye()
      @test op === op2
      @test op === op * op2 === op2 * op
    end

    @testset "Ones" begin
      E = opOnes(nrow, ncol)
      v = simple_vector(ComplexF64, nrow)
      u = simple_vector(ComplexF64, ncol)
      @test(norm(E * u - sum(u) * ones(nrow)) <= rtol * norm(u))
      @test(norm(transpose(E) * v - sum(v) * ones(ncol)) <= rtol * norm(v))
      @test(norm(E' * v - sum(v) * ones(ncol)) <= rtol * norm(v))
    end

    @testset "Zeros" begin
      O = opZeros(nrow, ncol)
      v = simple_vector(ComplexF64, nrow)
      u = simple_vector(ComplexF64, ncol)
      @test(norm(O * u) <= ϵ)
      @test(norm(transpose(O) * v) <= ϵ)
      @test(norm(O' * v) <= ϵ)
    end

    @testset "Diagonal" begin
      v = simple_vector(ComplexF64, nrow)
      D = opDiagonal(v)
      u = simple_vector(ComplexF64, nrow)
      @test(norm(D * u - v .* u) <= ϵ * norm(u))
      @test(norm(transpose(D) * u - v .* u) <= ϵ * norm(u))
      @test(norm(D' * u - conj(v) .* u) <= ϵ * norm(u))
    end

    @testset "Rectangular diagonal" begin
      nmin = min(nrow, ncol)
      nmax = max(nrow, ncol)
      A = zeros(ComplexF64, nmax, nmin)
      v = simple_vector(ComplexF64, nmin)
      for i = 1:nmin
        A[i, i] = v[i]
      end
      D = opDiagonal(nmax, nmin, v)
      u = simple_vector(ComplexF64, nmin)
      @test(norm(A * u - D * u) <= ϵ * norm(u))
      w = simple_vector(ComplexF64, nmax)
      @test(norm(transpose(A) * w - transpose(D) * w) <= ϵ * norm(w))
      @test(norm(A' * w - D' * w) <= ϵ * norm(w))

      A = zeros(ComplexF64, nmin, nmax)
      for i = 1:nmin
        A[i, i] = v[i]
      end
      D = opDiagonal(nmin, nmax, v)
      @test(norm(A * w - D * w) <= ϵ * norm(w))
      @test(norm(transpose(A) * u - transpose(D) * u) <= ϵ * norm(u))
      @test(norm(A' * u - D' * u) <= ϵ * norm(u))
    end

    @testset "posdef" begin
      v = simple_vector(ComplexF64, nrow)
      H = opHouseholder(v)
      Λ = collect(eltype(v), 1:nrow)
      op = H * opDiagonal(Λ) * H'
      @test check_positive_definite(op)
      @test check_positive_definite(op, semi = true)
      Λ = collect(eltype(v), 0:(nrow - 1))
      op = H * opDiagonal(Λ) * H'
      @test check_positive_definite(op, semi = true)
      A = Matrix(op)
      @test check_positive_definite(A, semi = true)
    end

    @testset "Hermitian" begin
      A = simple_matrix(ComplexF64, nrow, nrow)
      d = real.(diag(A))
      A = tril(A, -1)
      C = A + A' + diagm(0 => d)
      H = opHermitian(d, A)
      v = simple_vector(ComplexF64, nrow)
      @test(norm(H * v - C * v) <= rtol * norm(v))
      @test(norm(transpose(H) * v - transpose(C) * v) <= rtol * norm(v))
      @test(norm(H' * v - C * v) <= rtol * norm(v))

      @test(!check_hermitian(LinearOperator(A - A')))
      @test(!check_positive_definite(LinearOperator(-A' * A)))

      C = simple_matrix(ComplexF64, nrow, nrow, symmetric = true)
      H = opHermitian(C)
      v = simple_vector(ComplexF64, nrow)
      @test(norm(H * v - C * v) <= rtol * norm(v))
      @test(norm(transpose(H) * v - transpose(C) * v) <= rtol * norm(v))
      @test(norm(H' * v - C * v) <= rtol * norm(v))
    end

    @testset "Transpose and adjoint" begin
      A = simple_matrix(ComplexF64, nrow, nrow)
      v = simple_vector(ComplexF64, nrow)

      op = LinearOperator(nrow, nrow, false, false,
        (res, v, α, β) -> mul!(res, A, v, α, β),
        (res, v, α, β) -> mul!(res, transpose(A), v, α, β),
        nothing)
      @test(norm(transpose(A) * v - transpose(op) * v) <= rtol * norm(v))
      @test(norm(adjoint(A) * v - adjoint(op) * v) <= rtol * norm(v))
      @test(norm(A * v - transpose(transpose(op)) * v) <= rtol * norm(v))
      @test(norm(A * v - adjoint(adjoint(op)) * v) <= rtol * norm(v))
      @test(norm(conj.(A) * v - transpose(adjoint(op)) * v) <= rtol * norm(v))
      @test(norm(conj.(A) * v - adjoint(transpose(op)) * v) <= rtol * norm(v))

      op = LinearOperator(nrow, nrow, false, false,
        (res, v, α, β) -> mul!(res, A, v, α, β),
        nothing,
        (res, v, α, β) -> mul!(res, adjoint(A), v, α, β))
      @test(norm(transpose(A) * v - transpose(op) * v) <= rtol * norm(v))
      @test(norm(adjoint(A) * v - adjoint(op) * v) <= rtol * norm(v))
      @test(norm(A * v - transpose(transpose(op)) * v) <= rtol * norm(v))
      @test(norm(A * v - adjoint(adjoint(op)) * v) <= rtol * norm(v))
      @test(norm(conj.(A) * v - transpose(adjoint(op)) * v) <= rtol * norm(v))
      @test(norm(conj.(A) * v - adjoint(transpose(op)) * v) <= rtol * norm(v))
    end

    @testset "Integer" begin
      A = round.(Int, simple_matrix(Float64, nrow, nrow))
      op = LinearOperator(A)
      @test check_ctranspose(op)
      @test check_hermitian(op + op')
      @test check_positive_definite(op * op')
    end

    @testset "Restriction and Extension" begin
      n = 10
      J = [1; 2; 4; 7]
      r = 3:6
      s = 1:2:7
      k = 4
      v = simple_vector(Float64, nrow)

      for idx in (J, r, s, Colon(), k)
        P = opRestriction(idx, n)
        Z = opExtension(idx, n)

        # 1d slices are different; in Julia, v[idx] is a scalar
        w = v[idx]
        typeof(idx) <: Number && (w = [w])
        vz = zeros(n)
        vz[idx] = v[idx]

        @test P * v == w
        @test P' * w == vz
        @test Z * w == vz
        @test Z' * v == w
        @test (P * Z) * w == w
        @test (Z * P) * v == vz
      end
    end
  end

  @testset ExtendedTestSet "Linear system operators" begin
    A = simple_matrix(ComplexF64, nrow, nrow)
    v = simple_vector(Float64, nrow)

    @testset "Inverse" begin
      Ainv = opInverse(A)
      @test(norm(A \ v - Ainv * v) <= rtol * norm(v))
      @test(norm(transpose(A) \ v - transpose(Ainv) * v) <= rtol * norm(v))
      @test(norm(A' \ v - Ainv' * v) <= rtol * norm(v))
    end

    @testset "Cholesky and LDL" begin
      B = A' * A
      Binv = opCholesky(B)  #, check=true);
      @test(norm(B \ v - Binv * v) <= rtol * norm(v))
      @test(norm(transpose(B) \ v - transpose(Binv) * v) <= rtol * norm(v))
      @test(norm(B' \ v - Binv' * v) <= rtol * norm(v))

      @test_throws LinearOperatorException opCholesky(simple_matrix(Float64, 3, 5))
      @test_throws LinearOperatorException opCholesky(simple_matrix(Float64, 5, 5), check = true)

      # Test Cholesky operator on SQD matrix.
      A = simple_matrix(Float64, 3, 3, symmetric = true)
      B = simple_matrix(Float64, 2, 3)
      C = simple_matrix(Float64, 2, 2, symmetric = true)
      K = Symmetric([A B'; B -C])

      # Dense Cholesky should throw an exception.
      @test_throws LinearAlgebra.PosDefException opCholesky(K)

      # Compute the LDL' factorization.
      LDL = opLDL(sparse(K))
      e = ones(size(K, 1))
      @test(norm(LDL * (K * e) - e) < rtol * norm(e))
    end

    @testset "Householder" begin
      v = simple_vector(ComplexF64, nrow)
      H = opHouseholder(v)
      u = simple_vector(ComplexF64, nrow)
      @test(norm(H * u - (u - 2 * dot(v, u) * v)) <= rtol * norm(u))
      @test(norm(transpose(H) * u - (u - 2 * dot(conj(v), u) * conj(v))) <= rtol * norm(u))
      @test(norm(H' * u - (u - 2 * dot(v, u) * v)) <= rtol * norm(u))
    end
  end

  @testset ExtendedTestSet "Inference" begin
    function test_func(res)
      res .= 1.0 .+ im * 1.0
    end
    op = LinearOperator(5, 3, false, false, (res, p, α, β) -> test_func(res))
    @test eltype(op) == ComplexF64
    v = rand(5)
    @test_throws LinearOperatorException transpose(op) * v  # cannot be inferred
    @test_throws LinearOperatorException op' * v            # cannot be inferred

    op2 = conj(op)
    @test(norm(Matrix(op2) - conj(Matrix(op))) <= ϵ * norm(Matrix(op)))

    A = simple_matrix(ComplexF64, 5, 3)
    op = LinearOperator(A)
    @test(check_ctranspose(A))
    @test(check_ctranspose(op))
    @test_throws LinearOperatorException opCholesky(A)  # Shape mismatch

    A = simple_matrix(ComplexF64, 5, 5)
    @test_throws LinearOperatorException opCholesky(A, check = true)  # Not Hermitian / positive definite
    @test_throws LinearOperatorException opCholesky(-A' * A, check = true)  # Not positive definite

    # Adjoint of a symmetric non-hermitian
    A = simple_matrix(ComplexF64, 3, 3)
    A = A + transpose(A)
    op = LinearOperator(3, 3, true, false, (res, v, α, β) -> mul!(res, A, v))
    v = rand(3)
    @test op' * v ≈ A' * v
  end

  @testset ExtendedTestSet "Type specific operator" begin
    function prod(res, v, α, β)
      res[1] = v[1] + v[2]
      res[2] = v[2]
    end
    function ctprod(res, v, α, β)
      res[1] = v[1]
      res[2] = v[1] + v[2]
    end
    op = LinearOperator(2, 2, false, false, prod, nothing, ctprod)
    @test eltype(op) == Complex{Float64}
    for T in (Complex{Float64}, Complex{Float32}, BigFloat, Float64, Float32, Float16, Int32)
      op = LinearOperator(T, 2, 2, false, false, prod, nothing, ctprod)
      w = ones(T, 2)
      @test eltype(op) == T
      @test op * w == T[2; 1]
      @test eltype(op * w) == T
    end

    A = [im 1.0; 0.0 1.0]
    function prod!(res, v, α, β)
      mul!(res, A, v)
    end
    function tprod!(res, u, α, β)
      mul!(res, transpose(A), u)
    end
    function ctprod!(res, w, α, β)
      mul!(res, A', w)
    end
    opC = LinearOperator(2, 2, false, false, prod!, tprod!, ctprod!)
    v = simple_vector(ComplexF64, 2)
    @test A == Matrix(opC)
    opF = LinearOperator(Float64, 2, 2, false, false, prod!, tprod!, ctprod!) # The type is a lie
    @test eltype(opF) == Float64
    @test_throws InexactError Matrix(opF) # changed here TypeError to InexactError
  end

  # Issue #80
  @testset ExtendedTestSet "Test mul!" begin
    A = [1.0 1.0; 1.0 0.0]
    op = LinearOperator(A)
    y = zeros(2)
    x = ones(2)
    mul!(y, op, x)
    @test y == [2.0; 1.0]
  end

  # Issue #107
  @testset ExtendedTestSet "Unary and scalar operations on Adjoint and Transpose operators" begin
    op = LinearOperator(rand(5, 3))
    for adjtrans in [adjoint, transpose]
      @test Matrix(adjtrans(-op)) == Matrix(-adjtrans(op))
      @test Matrix(adjtrans(2 * op)) == Matrix(2 * adjtrans(op))
    end
  end

  # Issue #109
  @testset ExtendedTestSet "Sum with Adjoint and Transpose" begin
    A = rand(3, 3) + im * rand(3, 3)
    opA = LinearOperator(A)
    for adjtrans in [adjoint, transpose]
      @test Matrix(adjtrans(opA) + opA) == Matrix(opA + adjtrans(opA)) == A + adjtrans(A)
      @test Matrix(adjtrans(opA) + A) == Matrix(A + adjtrans(opA)) == A + adjtrans(A)
    end
  end

  # Issue #109
  @testset ExtendedTestSet "Cat with Adjoint and Transpose" begin
    A = rand(3, 3) + im * rand(3, 3)
    opA = LinearOperator(A)
    for adjtrans in [adjoint, transpose]
      @test Matrix([adjtrans(opA) opA]) == Matrix([adjtrans(A) A])
      @test Matrix([adjtrans(opA); opA]) == Matrix([adjtrans(A); A])
      @test Matrix([adjtrans(opA); adjtrans(opA)]) == Matrix([adjtrans(A); adjtrans(A)])
      @test Matrix([opA adjtrans(opA)]) == Matrix([A adjtrans(A)])
      @test Matrix([opA; adjtrans(opA)]) == Matrix([A; adjtrans(A)])
      @test Matrix([adjtrans(opA) A]) == Matrix([adjtrans(A) A])
      @test Matrix([adjtrans(opA); A]) == Matrix([adjtrans(A); A])
      @test Matrix([A adjtrans(opA)]) == Matrix([A adjtrans(A)])
      @test Matrix([A; adjtrans(opA)]) == Matrix([A; adjtrans(A)])
      @test Matrix([adjtrans(opA) opA; opA adjtrans(opA)]) == Matrix([adjtrans(A) A; A adjtrans(A)])
    end
  end

  @testset ExtendedTestSet "Counters" begin
    op = LinearOperator(rand(3, 4) + im * rand(3, 4))
    @test nprod(op) == 0
    @test ntprod(op) == 0
    @test nctprod(op) == 0
    nprods = 5
    ntprods = 4
    nctprods = 7
    for _ = 1:nprods
      op * rand(4)
    end
    for _ = 1:ntprods
      transpose(op) * rand(3)
    end
    for _ = 1:nctprods
      op' * rand(3)
    end
    @test nprod(op) == nprods
    @test ntprod(op) == ntprods
    @test nctprod(op) == nctprods
    for _ = 1:nprods
      conj(op) * rand(4)
    end
    @test nprod(op) == 2 * nprods

    opᵀ = transpose(op)
    @test nprod(opᵀ) == ntprod(op)
    @test ntprod(opᵀ) == nprod(op)
    @test nctprod(opᵀ) == nprod(op)

    opᴴ = op'
    @test nprod(opᴴ) == nctprod(op)
    @test ntprod(opᴴ) == nprod(op)
    @test nctprod(opᴴ) == nprod(op)

    reset!(op)
    @test nprod(op) == 0
    @test ntprod(op) == 0
    @test nctprod(op) == 0
  end

  @testset ExtendedTestSet "Timers" begin
    op = LinearOperator(rand(3, 4) + im * rand(3, 4))
    top = TimedLinearOperator(op)
    nprods = 5
    ntprods = 4
    nctprods = 7
    for _ = 1:nprods
      op * rand(4)
    end
    for _ = 1:ntprods
      transpose(op) * rand(3)
    end
    for _ = 1:nctprods
      op' * rand(3)
    end
    for fn ∈ (
      :size,
      :shape,
      :symmetric,
      :issymmetric,
      :hermitian,
      :ishermitian,
      :nprod,
      :ntprod,
      :nctprod,
    )
      @eval begin
        @test $fn($top) == $fn($top.op)
      end
    end

    reset!(op)
    reset!(top)

    top2 = TimedLinearOperator(op')  # the same as top'
    nrow, ncol = size(op)
    u = rand(nrow) + im * rand(nrow)
    @test all(top' * u .== top2 * u)
    v = rand(ncol) + im * rand(ncol)
    @test all(top * v .== top2' * v)

    top3 = TimedLinearOperator(transpose(op))  # the same as transpose(top)
    nrow, ncol = size(op)
    u = rand(nrow) + im * rand(nrow)
    @test all(transpose(top) * u .== top3 * u)
    v = rand(ncol) + im * rand(ncol)
    @test all(top * v .== transpose(top3) * v)
  end

  @testset ExtendedTestSet "BlockDiagonal" begin
    A = rand(3, 4) + im * rand(3, 4)
    B = rand(3, 3) + im * rand(3, 3)
    C = rand(4, 2) + im * rand(4, 2)
    D = [
      A zeros(3, 3) zeros(3, 2)
      zeros(3, 4) B zeros(3, 2)
      zeros(4, 4) zeros(4, 3) C
    ]
    M = BlockDiagonalOperator(LinearOperator.((A, B, C))...)
    @test size(M, 1) == size(A, 1) + size(B, 1) + size(C, 1)
    @test size(M, 2) == size(A, 2) + size(B, 2) + size(C, 2)
    @test norm(Matrix(M) - D) ≤ sqrt(eps()) * norm(D)
    @test norm(Matrix(transpose(M)) - transpose(D)) ≤ sqrt(eps()) * norm(D)
    @test norm(Matrix(M') - D') ≤ sqrt(eps()) * norm(D)
  end

  # Issue #139
  @testset ExtendedTestSet "Matrix-vector products with SparseMatrix and SparseVector" begin
    A = sprand(10, 10, 0.2)
    b = sprand(10, 0.2)
    opA = LinearOperator(A)
    @test opA * b == A * b
    @test transpose(opA) * b == transpose(A) * b
    @test adjoint(opA) * b == adjoint(A) * b
  end
end

test_linop()


text:
using MultivariateStats
using LinearAlgebra
using Test
using StableRNGs
import Statistics: mean, cov, var
import SparseArrays
import StatsBase

@testset "Probabilistic PCA" begin

    rng = StableRNG(34568)

    ## PCA with zero mean

    X = randn(rng, 5, 10)
    Y = randn(rng, 3, 10)

    W = qr(randn(rng, 5, 5)).Q[:, 1:3]
    σ² = 0.1
    M = PPCA(Float64[], W, σ²)

    @test size(M) == (5,3)
    @test mean(M) == zeros(5)
    @test loadings(M) == W
    @test var(M) == σ²

    T = inv(W'*W .+ σ²*Matrix(I, 3, 3))*W'
    @test predict(M, X[:,1]) ≈ T * X[:,1]
    @test predict(M, X) ≈ T * X

    R = W*inv(W'W)*(W'W .+ σ²*Matrix(I, 3, 3))
    @test reconstruct(M, Y[:,1]) ≈ R * Y[:,1]
    @test reconstruct(M, Y) ≈ R * Y

    ## PCA with non-zero mean

    mval = rand(rng, 5)
    M = PPCA(mval, W, σ²)

    @test size(M) == (5,3)
    @test mean(M) == mval
    @test loadings(M) == W
    @test var(M) == σ²

    @test predict(M, X[:,1]) ≈ T * (X[:,1] .- mval)
    @test predict(M, X) ≈ T * (X .- mval)

    @test reconstruct(M, Y[:,1]) ≈ R * Y[:,1] .+ mval
    @test reconstruct(M, Y) ≈ R * Y .+ mval

    ## prepare training data

    d = 5
    n = 1000

    R = collect(qr(randn(rng, d, d)).Q)
    @test R'R ≈ Matrix(I, 5, 5)
    rmul!(R, Diagonal(sqrt.([0.5, 0.3, 0.1, 0.05, 0.05])))

    X = R'randn(rng, 5, n) .+ randn(rng, 5)
    mval = vec(mean(X, dims=2))
    Z = X .- mval

    M0 = fit(PCA, X; mean=mval, maxoutdim = 4)

    ## ppcaml (default)

    M = fit(PPCA, X)
    P = projection(M)
    W = loadings(M)

    @test size(M) == (5,4)
    @test mean(M) == mval
    @test P'P ≈ Matrix(I, 4, 4)
    @test reconstruct(M, predict(M, X)) ≈ reconstruct(M0, predict(M0, X))

    M = fit(PPCA, X; mean=mval)
    @test loadings(M) ≈ W

    M = fit(PPCA, Z; mean=0)
    @test loadings(M) ≈ W

    M = fit(PPCA, X; maxoutdim=3)
    P = projection(M)
    W = loadings(M)

    @test size(M) == (5,3)
    @test P'P ≈ Matrix(I, 3, 3)

    # ppcaem

    M = fit(PPCA, X; method=:em)
    P = projection(M)
    W = loadings(M)

    @test size(M) == (5,4)
    @test mean(M) == mval
    @test P'P ≈ Matrix(I, 4, 4)
    @test all(isapprox.(reconstruct(M, predict(M, X)), reconstruct(M0, predict(M0, X)), atol=1e-2))

    M = fit(PPCA, X; method=:em, mean=mval)
    @test loadings(M) ≈ W

    M = fit(PPCA, Z; method=:em, mean=0)
    @test loadings(M) ≈ W

    M = fit(PPCA, X; method=:em, maxoutdim=3)
    P = projection(M)

    @test size(M) == (5,3)
    @test P'P ≈ Matrix(I, 3, 3)

    @test_throws StatsBase.ConvergenceException fit(PPCA, X; method=:em, maxiter=1)

    # bayespca
    M0 = fit(PCA, X; mean=mval, maxoutdim = 3)

    M = fit(PPCA, X; method=:bayes)
    P = projection(M)
    W = loadings(M)

    @test size(M) == (5,3)
    @test mean(M) == mval
    @test P'P ≈ Matrix(I, 3, 3)
    @test reconstruct(M, predict(M, X)) ≈ reconstruct(M0, predict(M0, X))

    M = fit(PPCA, X; method=:bayes, mean=mval)
    @test loadings(M) ≈ W

    M = fit(PPCA, Z; method=:bayes, mean=0)
    @test loadings(M) ≈ W

    M = fit(PPCA, X; method=:em, maxoutdim=2)
    P = projection(M)

    @test size(M) == (5,2)
    @test P'P ≈ Matrix(I, 2, 2)

    @test_throws StatsBase.ConvergenceException fit(PPCA, X; method=:bayes, maxiter=1)

    # Different data types
    # --------------------
    X = randn(rng, Float64, 5, 10)
    XX = convert.(Float32, X)

    Y = randn(rng, Float64, 1, 10)
    YY = convert.(Float32, Y)

    for method in (:bayes, :em)
        M = fit(PPCA, X ; maxoutdim=1, method=method)
        MM = fit(PPCA, XX ; maxoutdim=1, method=method)

        # mix types
        predict(M, XX)
        predict(MM, X)
        reconstruct(M, YY)
        reconstruct(MM, Y)

        # type consistency
        for func in (mean, projection, var, loadings)
            @test eltype(func(M)) == Float64
            @test eltype(func(MM)) == Float32
        end
    end

    # views
    X = randn(rng, 5, 200)
    M = fit(PPCA, view(X, :, 1:100), maxoutdim=3)
    M = fit(PPCA, view(X, :, 1:100), maxoutdim=3, method=:em)
    M = fit(PPCA, view(X, :, 1:100), maxoutdim=3, method=:bayes)
    # sparse
    @test_throws AssertionError fit(PPCA, SparseArrays.sprandn(rng, 100d, n, 0.6))

end


text:
# # Examples

## load packages
using BenchmarkTools
using WaterPhysics
FT = Float32;
#------------------------------------------------------------------------------

## define the variables
rand_r = (rand(FT) + 20) * FT(1e-6);
rand_T = rand(FT) + 298;
rand_α = rand(FT) * 50;
rand_Ψ = rand(FT) - 3;

gas_air = TraceGasAir{FT}();
gas_CO₂ = TraceGasCO₂{FT}();
#------------------------------------------------------------------------------

# ## Capillary pressure
@info "Benchmarking capillary_pressure functions...";
@btime capillary_pressure($rand_r, $rand_T);
@btime capillary_pressure($rand_r, $rand_T, $rand_α);
#------------------------------------------------------------------------------

# ## Diffusive coefficient
@info "Benchmarking diffusive_coefficient functions...";
@btime diffusive_coefficient($rand_T, $gas_CO₂, $gas_air);
@btime relative_diffusive_coefficient($rand_T);
#------------------------------------------------------------------------------

# ## Latent heat of evaporation
@info "Benchmarking latent_heat_vapor functions...";
@btime latent_heat_vapor($rand_T);
#------------------------------------------------------------------------------

# ## Saturation vapor pressure
@info "Benchmarking saturation_vapor_pressure functions...";
@btime pressure_correction($rand_T, $rand_Ψ);
@btime saturation_vapor_pressure($rand_T);
@btime saturation_vapor_pressure($rand_T, $rand_Ψ);
@btime saturation_vapor_pressure_slope($rand_T);
@btime saturation_vapor_pressure_slope($rand_T, $rand_Ψ);
#------------------------------------------------------------------------------

# ## Surface tension
@info "Benchmarking surface_tension functions...";
@btime surface_tension($rand_T);
@btime relative_surface_tension($rand_T);
#------------------------------------------------------------------------------

# ## Viscosity
@info "Benchmarking viscosity functions...";
@btime viscosity($rand_T);
@btime relative_viscosity($rand_T);
#------------------------------------------------------------------------------


text:
println("Compile libraries on main processor..."); flush(stdout)
include("imports.jl")
println("Finished loading libraries."); flush(stdout)

using Distributed
rmprocs(workers())
addprocs(20)
println("Compile libraries on workers ..."); flush(stdout)
@everywhere include("imports.jl")
println("Finished loading libraries."); flush(stdout)

sims = dict_list(Dict(:beta => collect(0:1),
                      :K => collect(2:6),
                      :snum => collect(1:4)))

# Run simulations in parallel.
@time res = pmap(run, sims, on_error=identity)

# Print results status.
foreach(z -> println("$(z[1]) => $(z[2])"), zip(sims, res))

# Post process
pp_res = pmap(postprocess, sims)
# postprocess(Dict(:beta=>1, :K=>5, :snum=>1))

# Compute BF
for d in dict_list(Dict(:K => collect(2:6), :snum => collect(1:4)))
  d0 = copy(d); d0[:beta] = 0
  d1 = copy(d); d1[:beta] = 1
  compute_bf(d0, d1)
end

# TODO: Combine results
for d in dict_list(Dict(:snum => collect(1:4)))
  combine_results(d)
end

# Model comparisons

# Send results to S3
# Util.s3sync(from="$(Info.resultsdir_simstudy)/$(simname)",
#             to="$(Info.awsbucket_simstudy)/$(simname)",
#             tags=`--exclude '*.nfs'`)


text:
using SparseArrays
using LinearAlgebra
using Random
using Clp
using MathOptInterface
const MOI = MathOptInterface

using TimerOutputs

struct RandomLP
    rows::Int
    cols::Int
    dens::Float64
end

function generate_moi_problem(model, At, b, c;
    var_bounds = false, scalar = true)
    cols, rows = size(At)
    x = MOI.add_variables(model, cols)
    A_cols = rowvals(At)
    A_vals = nonzeros(At)
    if var_bounds
        for col in 1:cols
            MOI.add_constraint(model, MOI.SingleVariable(x[col]),
                MOI.LessThan(10.0))
            MOI.add_constraint(model, MOI.SingleVariable(x[col]),
                MOI.GreaterThan(-10.0))
        end
    end
    if scalar
        for row in 1:rows
            MOI.add_constraint(model, MOI.ScalarAffineFunction(
                [MOI.ScalarAffineTerm(A_vals[i], x[A_cols[i]]) for i in nzrange(At, row)], 0.0),
                MOI.LessThan(b[row]))
        end
    else
        for row in 1:rows
            MOI.add_constraint(model, MOI.VectorAffineFunction(
                [MOI.VectorAffineTerm(1,
                    MOI.ScalarAffineTerm(A_vals[i], x[A_cols[i]])
                ) for i in nzrange(At, row)], [-b[row]]),
                MOI.Nonpositives(1))
        end
    end
    objective = MOI.ScalarAffineFunction(
        [MOI.ScalarAffineTerm(c[i], x[i]) for i in findall(!iszero, c)],
            0.0)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objective)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    return x
end

function random_data(seed, data)

    rows = data.rows
    cols = data.cols
    density = data.dens

    p_neg_element = 0.0 # 0.25

    rng = Random.MersenneTwister(seed)

    f_A(r, n) = ifelse.(rand(r, n) .> p_neg_element, 1, -1) .* (15 .+ 30 .* rand(r, n))

    At = sprand(rng, cols, rows, density, f_A)
    b = 50 * rand(rng, rows)

    # not using signs now
    sign = ifelse.(rand(rng, rows) .> 0.2, 'L', 'G')

    f_c(r, n) = 20 .* 2 .* (rand(r, n) .- 0.5)
    c = sprand(rng, cols, 0.5, f_c)

    return At, b, c
end

function bridged_cache_and_solver()
    model = MOI.Bridges.full_bridge_optimizer(MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        MOI.Utilities.MANUAL), Float64)
    clp = Clp.Optimizer()
    MOI.set(clp, MOI.Silent(), true)
    return model, clp
end
function cache_and_solver()
    model = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        MOI.Utilities.MANUAL)
    clp = Clp.Optimizer()
    MOI.set(clp, MOI.Silent(), true)
    return model, clp
end
function bridged_cached_solver()
    model = MOI.Bridges.full_bridge_optimizer(MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        Clp.Optimizer()), Float64)
    MOI.set(model, MOI.Silent(), true)
    return model
end
function cached_solver()
    model = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        Clp.Optimizer())
    MOI.set(model, MOI.Silent(), true)
    return model
end

function time_build_and_solve(to_build, to_solve, At, b, c, scalar = true)
    @timeit "build" x = generate_moi_problem(to_build, At, b, c, scalar = scalar)
    if to_build !== to_solve
        @timeit "copy" MOI.copy_to(to_solve, to_build, copy_names = false)
    end
    @time @timeit "opt" MOI.optimize!(to_solve)
    MOI.get(to_solve, MOI.ObjectiveValue())
    val = MOI.get(to_solve, MOI.SolveTime())
    println(val)
end

function solve_clp(seed, data; time_limit_sec=Inf)

    reset_timer!()

    for i in 1:seed

        @timeit "data"  At, b, c = random_data(i, data)

        bridged_cache, pure_solver = bridged_cache_and_solver()
        @timeit "bc + s" time_build_and_solve(bridged_cache, pure_solver, At, b, c)

        cache, pure_solver2 = cache_and_solver()
        @timeit "c + s" time_build_and_solve(cache, pure_solver2, At, b, c)

        full_solver = bridged_cached_solver()
        @timeit "bcs" time_build_and_solve(full_solver, full_solver, At, b, c)

        full_solver = bridged_cached_solver()
        @timeit "bcs + v" time_build_and_solve(full_solver, full_solver, At, b, c, false)

        cache_solver = cached_solver()
        @timeit "cs" time_build_and_solve(cache_solver, cache_solver, At, b, c)

    end

    print_timer()

end

solve_clp(10, RandomLP(10000, 20000, 0.01); time_limit_sec=5)


text:
# types.jl
using Dates:
    Date,
    DateTime,
    Period,
    Year,
    Month,
    Day,
    TimeType,
    UTInstant,
    value

struct QDate <: TimeType
    instant::UTInstant{Day}
    QDate(instant::UTInstant{Day}) = new(_check_instant(instant))
end
@inline function _check_instant(instant::UTInstant{Day})
    val = instant.periods.value
    if !(FIRST_VALUE <= val <= LAST_VALUE)
        throw(ArgumentError("Instant value: $val out of range ($FIRST_VALUE:$LAST_VALUE)"))
    end
    instant
end

@inline QDate(year::Integer, month::Integer=1, day::Integer=1) = QDate(year, month, false, day)
function QDate(year::Integer, month::Integer, leap::Bool, day::Integer)
    qdinfo = QREF.rqref_strict(year, month, leap, day)
    QDate(UTD(qdinfo.j - DAYS_OFFSET))
end
@inline _ci(x) = convert(Int, x)
@inline QDate(y,m=1,d=1) = QDate(_ci(y), _ci(m), _ci(d))
@inline QDate(y,m,l::Bool,d) = QDate(_ci(y), _ci(m), l, _ci(d))
@inline QDate(qdt::QDate) = qdt

@inline QDate(y::Year, m::Month=Month(1), d::Day=Day(1)) = QDate(value(y), value(m), false, value(d))
@inline QDate(y::Year, m::Month, l::Bool, d::Day=Day(1)) = QDate(value(y), value(m), l, value(d))

function QDate(periods::Union{Period,Bool}...)
    y = Year(0); m = Month(1); l = false; d = Day(1)
    _isyearspecified = false
    for p in periods
        isa(p, Year) && (_isyearspecified = true; y = p::Year)
        isa(p, Month) && (m = p::Month)
        isa(p, Bool) && (l = p::Bool)
        isa(p, Day) && (d = p::Day)
    end
    if !_isyearspecified
        throw(ArgumentError("Dates.Year must be specified"))
    end
    return QDate(y, m, l, d)
end

Base.eps(::QDate) = Day(1)

Base.typemax(::Union{QDate,Type{QDate}}) = QDate(UTD(LAST_VALUE))
Base.typemin(::Union{QDate,Type{QDate}}) = QDate(UTD(FIRST_VALUE))

Base.isless(x::QDate, y::QDate) = isless(value(x), value(y))

Base.promote_rule(::Type{QDate}, x::Type{Date}) = Date
Base.promote_rule(::Type{QDate}, x::Type{DateTime}) = DateTime

# for convenience
function QREF.qref(qdt::QDate)
    QREF.qref(value(qdt) + DAYS_OFFSET)
end


text:
struct SBM
    w::Matrix{Float64}      # Matrix of probabilities
    dist::String            # Distribution (either "poisson" or "bernoulli")
    n_communities::Int      # Number of communities/clusters
    variant::String

    # Constructor function for struct SBM
    function SBM(w::Matrix{Float64}, dist::String, variant::String)
        # Check if w is a square matrix
        if size(w,1) != size(w,2)
            throw(ArgumentError("Matrix w must be a square matrix."))
        end
        # Check if elements in w are strictly positive
        if any(w .< 0)
            throw(DomainError("All elements in matrix of probabilities w must be strictly positive. Got w = $w"))
        end
        # Check if matrix is symmetric
        if ~isapprox(w, w', rtol=1e-4) # TODO: define global variable for rtol
            throw(ArgumentError("Matrix w must be symmetric. Got w = $w"))
        end
        q = size(w,1)

        # Check input argument 'distribution'
        if ~(dist in ["poisson","bernoulli"])
            throw(ArgumentError(string("Invalid distribution: ", dist)))
        end
        if dist == "bernoulli"
            # Check if w is a matrix of 0-1 values
            if any(w .> 1)
                throw(DomainError("When using a bernoulli distribution all elements in matrix of probabilities w must be between 0 and 1. Got w = $w"))
            end
        end
        if ~(variant in ["standard","degree-corrected"])
            throw(ArgumentError("Invalid SBM variant: $variant"))
        end
        return new(w, dist, q, variant)
    end

    function SBM(w::Matrix{Float64}, dist::String)
        variant = "standard"
        return SBM(w, dist, variant)
    end
end

function generate(sbm::SBM, n_per_community::Vector{Int}, seed::Union{Nothing,Int}=nothing)::Dataset
    """
    Generates a graph from the Stochastic Block Model.

    Parameters
    ----------
    sbm : SBM
        An SBM structure.
    n_per_community : Vector{Int}
        Number of nodes in each community.
    seed : Union{Nothing,Int}
        Seed for the random number generator.

    Returns
    -------
    adj_matrix : Matrix{Int}
        Adjacency matrix of the generated graph
    """
    if sbm.variant != "standard"
        throw("not implemented")
    end

    # Check input argument 'n_per_community'
    if any(n_per_community .<= 0)
        throw(DomainError(string("Invalid value in argument n_per_community: ",
                                 "the number of nodes in each community must be strictly positive.")))
    end
    # Check if n_per_community matches the dimensions of w
    if ~(length(n_per_community) == sbm.n_communities)
        throw(ArgumentError("Number of communities in n_per_community must match the number of communities in SBM"))
    end

    # Set the random seed if one was provided
    if seed != nothing
        Random.seed!(seed)
    end
    # Define number of vertices, initialize adjacency matrix and community labels
    n_vertices = sum(n_per_community);
    adj_matrix = zeros(Int, n_vertices,n_vertices);
    community_labels = convert(Array{Int}, vcat([i*ones(ni) for (i,ni) in enumerate(n_per_community)]...));
    # Construct adjacency matrix
    for i=1:n_vertices
        for j=1:i
            proba = sbm.w[community_labels[i], community_labels[j]];
            if sbm.dist == "poisson"
                if i == j
                    # Note that this implies that the expected number of self-edges at a vertex in group r is 1/2 * ωrr
                    # because of the factor of two in the definition of the diagonal elements of the adjacency matrix.
                    proba = proba / 2 # "The factor of half is included solely because it makes the algebra easier"
                end
                pois = Distributions.Poisson(proba);
                adj_matrix[i,j] = rand(pois);
            elseif sbm.dist == "bernoulli"
                if i == j
                    adj_matrix[i,j] = 0;
                else
                    bernoulli = Bernoulli(proba);
                    adj_matrix[i,j] = rand(bernoulli);
                end
            end
        end
    end
    # "We have adopted the common convention that a self-edge is represented by A[i,j] = 2
    # (and not 1 as one might at first imagine)"
    adj_matrix += transpose(adj_matrix) #- Diagonal(adj_matrix)

    # Instantiate a dataset
    n_communities = length(n_per_community)
    dataset = Dataset(adj_matrix, n_communities)
    return dataset
end

function generate(probability_matrix::Matrix{Float64}, n_per_community::Vector{Int};
                            distribution::String="poisson", seed::Union{Nothing, Int}=nothing)::Dataset
    """
    Generates a graph from the Stochastic Block Model.

    Parameters
    ----------
    probability_matrix : Matrix{Float64}
        Matrix of expected number of edges under Poisson distribution.
    n_per_community : Vector{Int}
        Number of nodes in each community.
    distribution : String
        Distribution to be used for generating edges (either "poisson" or "bernoulli").
    seed : Union{Nothing,Int}
        Seed for the random number generator.

    Returns
    -------
    dataset : Dataset
        Dataset representing a graph.
    """
    sbm = SBM(probability_matrix, distribution)
    return generate(sbm, n_per_community, seed)
end

function getObjectiveConstant(dataset::Dataset)::Float64
    """ Calculate constant term of the objective function """
    A = dataset.A
    n = dataset.n
    m = dataset.m
    k = dataset.k

    Cnst = 0
    for i=1:n, j=1:(i-1)
        if A[i,j] == 0
            Cnst += -1
        else
            Cnst += ( A[i,j] * log(k[i]*k[j] / (2*m)) - log(factorial(A[i,j])) )
        end
    end

    for i=1:n
        if A[i,i] == 0
            Cnst += -1
        else
            Cnst += ( 0.5 * A[i,i] * log((0.5*k[i]^2)/(2*m)) - log(0.5*factorial(A[i,i])))
        end
    end
    return Cnst
end

function getWUpperBound(dataset::Dataset)::Float64
    # Upper bound on w
    n = dataset.n
    m = dataset.m
    A = dataset.A
    k = dataset.k

    rho = 0
    for i=1:n, j=1:n
        if A[i,j] != 0
            rho = maximum([rho, (A[i,j] / (k[i]*k[j]))])
        end
    end
    w_upper = 2 * m * rho
    return w_upper
end

function calculateObjective(dataset::Dataset, w::Matrix{Float64}, z::Matrix{Int})::Float64
    """ Calculates the maximum log-likelihood value

    Parameters
    ----------
    dataset : Dataset
        Dataset representing an observed graph.
    sbm : SBM
        An Stochastic Block Model
    z : Matrix{Int}
        Matrix of assignments of nodes to groups

    Returns
    -------
    L : Float64
        Maximum log-likelihood value
    """
    A = dataset.A
    n = dataset.n
    m = dataset.m
    q = dataset.n_communities
    k = dataset.k
    if ~((n,q) == size(z))
        throw(ArgumentError("The dimensions of the assignment matrix do not match the dimensions of the dataset."))
    end

    # Get objective constant
    Cnst = getObjectiveConstant(dataset)

    # Calculate objective value
    L = 0
    for g=1:q, h=1:q, i=1:n, j=1:n
        if z[i,g] * z[j,h] == 1
            if A[i,j] != 0
                L += 0.5*(- A[i,j]*log(w[g,h]) + ((k[i]*k[j])/(2*m)) * w[g,h])
            else
                L += 0.5*( ((k[i]*k[j])/(2*m)) * w[g,h] )
            end
        end
    end
    L = L - Cnst
    return L
end

max_stars_repo_name: vidalt/Exact-SBM

max_issues_repo_name: vidalt/Exact-SBM

max_forks_repo_name: vidalt/Exact-SBM


text:
using DelayEmbeddings, DynamicalSystemsBase
using Test
using Random
import Peaks

@testset "Pecora" begin
# %% Generate data
lor = Systems.lorenz([0.0;1.0;0.0];ρ=60)
data = trajectory(lor, 200; Δt=0.02, Ttr = 10)
metric = Chebyshev()

UNDERSAMPLING = false

@testset "Pecora univariate" begin
# %% Timeseries case

s = data[:, 2] # input timeseries = first entry of lorenz
optimal_τ = estimate_delay(s, "mi_min")
Tmax = 100
K = 14
samplesize = 1

# using PyPlot
# pygui(true)
# figure()
# ax1 = subplot(211)
# ylabel("⟨ε★⟩")
# ax2 = subplot(212)
# xlabel("τ (index units)")
# ylabel("⟨Γ⟩")

τs = (0,)
Random.seed!(123)
es_ref, Γs = pecora(s, τs; delays = 0:Tmax, w = optimal_τ, samplesize = samplesize, K = K, metric = metric, undersampling = UNDERSAMPLING)
(max1,_) = Peaks.findmaxima(vec(es_ref))
@test optimal_τ - 2 ≤ max1[1]-1 ≤ optimal_τ + 2
maxi = maximum(es_ref)
@test maxi ≤ 1.3
# ax1.plot(es_ref, label = "τs = $(τs)")
# ax2.plot(Γs)

τs = (0, max1[1]-1,)
Random.seed!(123)
es, Γs = pecora(s, τs; delays = 0:Tmax, w = optimal_τ, samplesize = samplesize, K = K, metric = metric, undersampling = UNDERSAMPLING)
(max2,_) = Peaks.findmaxima(vec(es))
(min1,_) = Peaks.findminima(vec(es))
@test 1 ≤ max2[1]-1 ≤ 4
@test optimal_τ - 2 ≤ min1[1]-1 ≤ optimal_τ + 2
# ax1.plot(es, label = "τs = $(τs)")
# ax2.plot(Γs)

τs = (0, max1[1]-1, max2[1]-1)
Random.seed!(123)
es, Γs = pecora(s, τs; delays = 0:Tmax, w = optimal_τ, samplesize = samplesize, K = K, metric = metric, undersampling = UNDERSAMPLING)
(min2,_) = Peaks.findminima(vec(es))
@test max2[1]-2 ≤ min2[1]-1 ≤ max2[1]
# ax1.plot(es, label = "τs = $(τs)")
# ax2.plot(Γs)
# ax1.legend()
# ax1.set_title("lorenz data")

end

@testset "Pecora multivariate" begin
## %% Trajectory case
s = Dataset(data)
optimal_τ = estimate_delay(s[:,2], "mi_min")
Tmax = 100
K = 14
samplesize = 1

js = (2,)
τs = (0,)
Random.seed!(123)
es_ref, Γs = pecora(s[:,2], τs; delays = 0:Tmax, w = optimal_τ, samplesize = samplesize, K = K, metric = metric, undersampling = UNDERSAMPLING)
Random.seed!(123)
es, Γs = pecora(s, τs, js; delays = 0:Tmax, samplesize = samplesize, w = optimal_τ, K = K, metric = metric, undersampling = UNDERSAMPLING)

@test round.(es[:,2], digits = 4) == round.(vec(es_ref),digits = 4)
(x_maxi,_) = Peaks.findmaxima(es[:,1])
@test 8 ≤ x_maxi[1]-1 ≤ 10
(z_maxi,_) = Peaks.findmaxima(es[:,3])
@test 13 ≤ z_maxi[1]-1 ≤ 15

# using PyPlot
# pygui(true)
# figure()
# subplot(211)
# for i in 1:3
#     plot(es[:, i], label = "τs = $(τs), js = $(js), J = $(i)")
# end
# ylabel("⟨ε★⟩")
# title("lorenz system, multiple timeseries")
# legend()
# grid()
# subplot(212)
# plot(es_ref, label = "τs = $(τs), js = $(js), J = 2")
# ylabel("⟨ε★⟩-ref")
# xlabel("τ (index units)")
# grid()

end
end


text:
# Helper functions

function generic_LL_V(data::Array{Int,2}, lfun)

    data_stats = get_data_stats(data)
    LL = 0.
    all_signals = SortedDict(Dict{Int,Step_by_step_signals}())

    for subj_id in data_stats.all_subject_ids
        subj_data = filter_by_person(data, subj_id)

        LL_subj, signals_subj = lfun(subj_id, subj_data)

        LL += LL_subj
        all_signals[signals_subj.subject_id] = signals_subj
    end
    LL, all_signals
end

function get_MAP(samples, top_n)
    idx_perm = sortperm(samples[:,end],rev=true)

    top_n = idx_perm[1:top_n]
    top_n_samples = samples[top_n,:]
    top_MAP = top_n_samples[1,:]

    top_n_avgMAP = mean(top_n_samples, dims=1)
    top_n_stdMAP = std(top_n_samples, dims=1)
    top_n_min = minimum(top_n_samples, dims=1)
    top_n_max = maximum(top_n_samples, dims=1)
    top_MAP, top_n_avgMAP, top_n_min, top_n_max, top_n_stdMAP, top_n_samples
end

function get_idx_marginal_posterior(mcmc_samples, idx_param_of_interest, param_Map_values, cond_width)
    # identifies the samples in mcmc_samples which are closest around the
    # marginal posterior.
    # for each parameter i (except the parameter of interest), only the
    # samples within param_Map_values[i] +- cond_width/2. are kept.
    # Used mostly for visualization purproses
    nr_params = length(param_Map_values)
    pram_indexes = setdiff(1:nr_params, idx_param_of_interest)
    m,n = size(mcmc_samples)

    idxs = trues(m)
    for i in pram_indexes
        parval_min = param_Map_values[i] - cond_width/2.
        parval_max = param_Map_values[i] + cond_width/2.
        idx_param = (mcmc_samples[:,i] .>= parval_min) .& (mcmc_samples[:,i] .<= parval_max)
        idxs .&= idx_param
    end
    return idxs
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
data for electric tariff in reopt model
    can be defined using custom rates or URDB rate
"""
# TODO function for creating BAU inputs: don't need for tariff b/c Tech repeats no longer needed?
struct ElectricTariff
    energy_rates::Array{Float64,1}

    monthly_demand_rates::Array{Float64,1}
    time_steps_monthly::Array{Array{Int64,1},1}  # length = 0 or 12

    tou_demand_rates::Array{Float64,1}
    tou_demand_ratchet_timesteps::Array{Array{Int64,1},1}  # length = n_tou_demand_ratchets

    fixed_monthly_charge::Float64
    annual_min_charge::Float64
    min_monthly_charge::Float64

    export_rates::DenseAxisArray{Array{Float64,1}}
    export_bins::Array{Symbol,1}
end

"""
    ElectricTariff

    function ElectricTariff(;
        urdb_label::String="",
        urdb_response::Dict=Dict(),
        year::Int=2020,
        time_steps_per_hour::Int=1,
        NEM::Bool=false,
        wholesale_rate::T=nothing,
        monthly_energy_rates::Array=[],
        monthly_demand_rates::Array=[],
    ) where {T <: Union{Nothing, Int, Float64, Array}, S <: Union{Nothing, Int, Float64, Array}}

"""
function ElectricTariff(;
    urdb_label::String="",
    urdb_response::Dict=Dict(),
    year::Int=2020,
    time_steps_per_hour::Int=1,
    NEM::Bool=false,
    wholesale_rate::T=nothing,
    monthly_energy_rates::Array=[],
    monthly_demand_rates::Array=[],
    ) where {T <: Union{Nothing, Int, Float64, Array}, S <: Union{Nothing, Int, Float64, Array}}

    nem_rate = [100.0 for _ in 1:8760*time_steps_per_hour]

    u = nothing
    if !isempty(urdb_label)

        u = URDBrate(urdb_label, year, time_steps_per_hour=time_steps_per_hour)

    elseif !isempty(urdb_response)

        u = URDBrate(urdb_response, year, time_steps_per_hour=time_steps_per_hour)

    elseif !isempty(monthly_energy_rates) && !isempty(monthly_demand_rates)

        invalid_args = String[]
        if !(length(monthly_energy_rates) == 12)
            push!(invalid_args, "length(monthly_energy_rates) must equal, got length $(length(monthly_energy_rates))")
        end
        if !(length(monthly_demand_rates) == 12)
            push!(invalid_args, "length(monthly_demand_rates) must equal, got length $(length(monthly_demand_rates))")
        end
        if length(invalid_args) > 0
            error("Invalid argument values: $(invalid_args)")
        end

        tou_demand_rates = Float64[]
        tou_demand_ratchet_timesteps = []
        demand_rates_monthly = monthly_demand_rates
        time_steps_monthly = get_monthly_timesteps(year, time_steps_per_hour=time_steps_per_hour)
        energy_rates = Real[]
        for m in 1:12
            append!(energy_rates, [monthly_energy_rates[m] for ts in time_steps_monthly[m]])
        end

        fixed_monthly_charge = 0.0
        annual_min_charge = 0.0
        min_monthly_charge = 0.0

        if NEM
            nem_rate = [-0.999 * x for x in energy_rates]
        end
    else
        error("Creating ElectricTariff requires at least urdb_label or monthly rates.")
    end

    if !isnothing(u)

        if NEM
            t = get_tier_with_lowest_energy_rate(u)
            nem_rate = [-0.999 * x for x in u.energy_rates[t,:]]
        end

        energy_rates, monthly_demand_rates, tou_demand_rates = remove_tiers_from_urdb_rate(u)
        time_steps_monthly = Array[]
        if !isempty(u.monthly_demand_rates)
            time_steps_monthly =
                get_monthly_timesteps(year, time_steps_per_hour=time_steps_per_hour)
        end

        tou_demand_ratchet_timesteps = u.tou_demand_ratchet_timesteps
        fixed_monthly_charge = u.fixed_monthly_charge
        annual_min_charge = u.annual_min_charge
        min_monthly_charge = u.min_monthly_charge
    end

    #= export_rates
    3 "tiers": 1. NEM (Net Energy Metering), 2. WHL (Wholesale), 3. CUR (Curtail)
    - if NEM then set ExportRate[:Nem, :] to energy_rate[tier_with_lowest_energy_rate, :]
        - otherwise set to 100 dollars/kWh
    - user can provide either scalar wholesale rate or vector of timesteps,
        - otherwise set to 100 dollars/kWh
    - curtail cost set to zero by default, but can be specified same as wholesale rate
    =#
    whl_rate = create_export_rate(wholesale_rate, length(energy_rates), time_steps_per_hour)

    if !NEM & (sum(whl_rate) >= 0)
        export_rates = DenseAxisArray{Array{Float64,1}}(undef, [])
        export_bins = Symbol[]
    elseif !NEM
        export_bins = [:WHL]
        export_rates = DenseAxisArray([whl_rate], export_bins)
    elseif (sum(whl_rate) >= 0)
        export_bins = [:NEM]
        export_rates = DenseAxisArray([nem_rate], export_bins)
    else
        export_bins = [:NEM, :WHL]  # TODO add :EXC to align with REopt Lite API
        export_rates = DenseAxisArray([nem_rate, whl_rate], export_bins)
    end

    ElectricTariff(
        energy_rates,
        monthly_demand_rates,
        time_steps_monthly,
        tou_demand_rates,
        tou_demand_ratchet_timesteps,
        fixed_monthly_charge,
        annual_min_charge,
        min_monthly_charge,
        export_rates,
        export_bins
    )
end

function get_tier_with_lowest_energy_rate(u::URDBrate)
    """
    ExportRate should be lowest energy cost for tiered rates.
    Otherwise, ExportRate can be > FuelRate, which leads REopt to export all PV energy produced.
    """
    tier_with_lowest_energy_cost = 1
    if length(u.energy_tier_limits) > 1
        annual_energy_charge_sums = Float64[]
        for etier in u.energy_rates
            push!(annual_energy_charge_sums, sum(etier))
        end
        tier_with_lowest_energy_cost =
            findall(annual_energy_charge_sums .== minimum(annual_energy_charge_sums))[1]
    end
    return tier_with_lowest_energy_cost
end

# TODO: dispatch custom rates based on options: TOU, monthly, etc.
function CustomRate(flat_energy::Real, flat_demand::Real=0) end

"""
    function create_export_rate(e::Nothing, N::Int, ts_per_hour::Int=1)
No export rate provided by user: set to 100 dollars/kWh for all time
"""
function create_export_rate(e::Nothing, N::Int, ts_per_hour::Int=1)
    [100 for _ in range(1, stop=N) for ts in 1:ts_per_hour]
end

"""
    function create_export_rate(e::T, N::Int, ts_per_hour::Int=1) where T<:Real
Case for scaler export rate provided -> convert to array of timesteps
"""
function create_export_rate(e::T, N::Int, ts_per_hour::Int=1) where T<:Real
    [float(-1*e) for ts in range(1, stop=N) for ts_per_hour::Int=1]
end

"""
    function create_export_rate(e::Array{<:Real, 1}, N::Int, ts_per_hour::Int=1)

Check length of e and upsample if length(e) != N
"""
function create_export_rate(e::Array{<:Real, 1}, N::Int, ts_per_hour::Int=1)
    Ne = length(e)
    if Ne != Int(N/ts_per_hour) || Ne != N
        @error "Export rates do not have correct number of entries. Must be $(N) or $(Int(N/ts_per_hour))."
    end
    if Ne != N  # upsample
        export_rates = [-1*x for x in e for ts in 1:ts_per_hour]
    else
        export_rates = -1*e
    end
    return export_rates
end

"""
    get_monthly_timesteps(year::Int; time_steps_per_hour=1)

return Array{Array{Int64,1},1}, size = (12,)
"""
function get_monthly_timesteps(year::Int; time_steps_per_hour=1)
    a = Array[]
    i = 1
    for m in range(1, stop=12)
        n_days = daysinmonth(Date(string(year) * "-" * string(m)))
        stop = n_days * 24 * time_steps_per_hour + i - 1
        if m == 2 && isleapyear(year)
            stop -= 24 * time_steps_per_hour  # TODO support extra day in leap years?
        end
        steps = [step for step in range(i, stop=stop)]
        append!(a, [steps])
        i = stop + 1
    end
    return a
end

# TODO use this function only for URDBrate
function remove_tiers_from_urdb_rate(u::URDBrate)
    # tariff args: have to validate that there are no tiers
    if length(u.energy_tier_limits) > 1
        @warn "Energy rate contains tiers. Using the first tier!"
    end
    elec_rates = vec(u.energy_rates[1,:])

    if u.n_monthly_demand_tiers > 1
        @warn "Monthly demand rate contains tiers. Using the last tier!"
    end
    if u.n_monthly_demand_tiers > 0
        demand_rates_monthly = vec(u.monthly_demand_rates[:,u.n_monthly_demand_tiers])
    else
        demand_rates_monthly = vec(u.monthly_demand_rates)  # 0×0 Array{Float64,2}
    end

    if u.n_tou_demand_tiers > 1
        @warn "TOU demand rate contains tiers. Using the last tier!"
    end
    if u.n_tou_demand_tiers > 0
        demand_rates = vec(u.tou_demand_rates[:,u.n_tou_demand_tiers])
    else
        demand_rates = vec(u.tou_demand_rates)
    end

    return elec_rates, demand_rates_monthly, demand_rates
end

max_stars_repo_name: NLaws/REoptLite

max_issues_repo_name: NLaws/REoptLite

max_forks_repo_name: NLaws/REoptLite