text:
macro augmodel()
    @error "You need to provide a list of arguments, please check the documentation on the macro"
end

"""
** Likelihood**
Template file for likelihood creation

```julia
()
```
See all functions you need to implement
---

"""
const InputSymbol = Union{Symbol,Expr,Real}

check_model_name(name::Symbol) = !isnothing(match(r"[^\w]*", string(name)))

function treat_likelihood(likelihood::Expr) end

function treat_params(params) end

function correct_parenthesis(text::AbstractString)
    return replace(text, r"(?=[()*+-.])" => "\\")
end

function check_likelihoodtype(ltype)
    return ltype == :Regression || ltype == :Classification || ltype == :Event
end

const AGP = AugmentedGaussianProcesses

macro augmodel(name::Symbol, likelihoodtype::Symbol, likelihood::Expr, params)
    latent = :x
    @assert occursin(r"p\(y\s*|\s*\w\)", string(likelihood)) "Likelihood should be of the form p(y|x) = C*exp(g(y)*x)*φ(α(y)-β(y)*x+γ(y)*x^2), replacing all functions by 0 if necessary"
    pdf_string = string(likelihood.args[2].args[2])
    # @show pdfexpr = Meta.parse(pdfstring)
    C_string = match(r".*?(?= (\* exp\(.*x\)))", pdf_string).match
    G_string_f =
        match(
            Regex("(?<=$(AGP.correct_parenthesis(C_string)) \\* exp\\().*(?=x\\) \\*)"),
            pdf_string,
        ).match
    G_string = deepcopy(G_string_f)
    while last(G_string) == ' ' || last(G_string) == '*'
        G_string = G_string[1:(end - 1)]
    end
    phi_h_string =
        match(
            Regex("(?<=$(AGP.correct_parenthesis(G_string_f))x\\) \\* ).*"), pdf_string
        ).match
    loc_x² = findfirst("x ^ 2", phi_h_string)
    count_parenthesis = 1
    loc_start = loc_x²[1]
    while count_parenthesis != 0
        loc_start = loc_start - 1
        if phi_h_string[loc_start] == ')'
            count_parenthesis += 1
        elseif phi_h_string[loc_start] == '('
            count_parenthesis -= 1
        end
    end
    h_string = phi_h_string[(loc_start + 1):loc_x²[end]]
    phi_string = phi_h_string[1:loc_start] * "r" * phi_h_string[(loc_x²[end] + 1):end]
    @show alpha_string = match(r"[^(][^-]*", h_string).match[1:(end - 1)]
    gamma_string = match(r"(?<=\+ )[^x]*(?=x \^ 2)", h_string).match
    gamma_string = gamma_string == "" ? "1.0" : gamma_string
    count_parenthesis -= 1
    h_string = phi_h_string[(loc_start + 1):loc_x²[end]]
    phi_string = phi_h_string[1:loc_start] * "r" * phi_h_string[(loc_x²[end] + 1):end]
    @show alpha_string = match(r"[^(][^-]*", h_string).match[1:(end - 1)]
    gamma_string = match(r"(?<=\+ )[^x]*(?=x \^ 2)", h_string).match
    gamma_string = gamma_string == "" ? "1.0" : gamma_string
    while last(gamma_string) == ' ' || last(gamma_string) == '*'
        gamma_string = gamma_string[1:(end - 1)]
    end
    beta_string =
        match(
            Regex("(?<=$(AGP.correct_parenthesis(alpha_string)) -  )[^( x )]*(?= x)"),
            h_string,
        ).match
    while last(beta_string) == ' ' || last(beta_string) == '*'
        beta_string = beta_string[1:(end - 1)]
    end
    return (C_string, G_string, phi_string, alpha_string, beta_string, gamma_string)
    # treat_params(params)
    # C,g,α,β,γ,φ = treat_likelihood(likelihood)
end

macro augmodel(name::Symbol, likelihoodtype::Symbol, likelihood::Expr) end

macro augmodel(
    name,
    ltype,
    C::InputSymbol,
    g::InputSymbol,
    α::InputSymbol,
    β::InputSymbol,
    γ::InputSymbol,
    φ::InputSymbol,
    ∇φ::InputSymbol,
    args...,
)
    add_variables = []
    default_values = []
    for input in args
        @assert input.head == :(=) "Additional variables should be given a default value, for example `b=1`"
        push!(add_variables, input.args[1])
        push!(default_values, input.args[1])
    end
    #### Check args here
    @assert check_model_name(name) "Please only use alphabetic characters for the name of the likelihood"
    @assert check_likelihoodtype(ltype) "Please use a correct likelihood type : Regression, Classification or Event"
    #Find gradient with AD if needed
    return esc(
        generate_likelihood(
            Symbol(name, "Likelihood"), Symbol(ltype, "Likelihood"), C, g, α, β, γ, φ, ∇φ
        ),
    )
end
function generate_likelihood(lname, ltype, C, g, α, β, γ, φ, ∇φ)
    quote
        begin
            using Statistics
            struct $(lname){T<:Real,A<:AbstractVector{T}} <: AGP.$(ltype){T}
                c²::A
                θ::A
                function $(lname){T}() where {T<:Real}
                    return new{T,Vector{T}}()
                end
                function $(lname){T}(c²::A, θ::A) where {T<:Real,A<:AbstractVector{T}}
                    return new{T,A}(c², θ)
                end
            end

            function $(lname)()
                return $(lname){Float64}()
            end

            function AGP.implemented(
                ::$(lname), ::Union{<:AnalyticVI,<:QuadratureVI,<:GibbsSampling}
            )
                return true
            end

            function AGP.init_likelihood(
                likelihood::$(lname){T},
                inference::AbstractInference{T},
                nLatent::Int,
                nSamplesUsed::Int,
            ) where {T}
                if inference isa AnalyticVI || inference isa GibbsSampling
                    $(lname){T}(zeros(T, nSamplesUsed), zeros(T, nSamplesUsed))
                else
                    $(lname){T}()
                end
            end

            function (l::$(lname))(y::Real, f::Real)
                return _gen_C(l) *
                       exp(_gen_g(l, y) * f) *
                       _gen_φ(l, _gen_α(l, y) - _gen_β(l, y) * f + _gen_γ(l, y) * f^2)
            end

            function Distributions.loglikelihood(l::$(lname), y::Real, f::Real)
                return log(_gen_C(l)) +
                       _gen_g(l, y) * f +
                       log(_gen_φ(l, _gen_α(l, y) - _gen_β(l, y) * f + _gen_γ(l, y) * f^2))
            end

            function _gen_C(l::$(lname){T}) where {T<:Real}
                return $(C)
            end

            function _gen_g(l::$(lname), y::T) where {T<:Real}
                return $(g)
            end

            function _gen_α(l::$(lname), y::T) where {T<:Real}
                return $(α)
            end

            function _gen_β(l::$(lname), y::T) where {T<:Real}
                return $(β)
            end

            function _gen_γ(l::$(lname), y::T) where {T<:Real}
                return $(γ)
            end

            function _gen_φ(l::$(lname), r::T) where {T<:Real}
                return $(φ)
            end

            function _gen_∇φ(l::$(lname), r::T) where {T<:Real}
                return $(∇φ)
            end

            function _gen_∇²φ(l::$(lname), r::T) where {T}
                return Zygote.gradient(x -> $(∇φ)(x), r)[1]
            end

            function Base.show(io::IO, model::$(lname){T}) where {T<:Real}
                return print(io, "$(nameof(typeof(model)))")
            end

            function Statistics.var(l::$(lname){T}) where {T<:Real}
                @warn "The variance of the likelihood is not implemented : returnin 0.0"
                return 0.0
            end

            function AGP.compute_proba(
                l::$(lname){T}, μ::AbstractVector{T}, σ²::AbstractVector{T}
            ) where {T<:Real}
                if typeof(l) <: RegressionLikelihood
                    return μ, max.(σ², zero(T)) .+ var(l)
                else
                    N = length(μ)
                    pred = zeros(T, N)
                    sig_pred = zeros(T, N)
                    for i in 1:N
                        x = AGP.pred_nodes .* sqrt(max(σ²[i], zero(T))) .+ μ[i]
                        pred[i] = dot(AGP.pred_weights, AGP.pdf.(l, 1, x))
                        sig_pred[i] =
                            dot(AGP.pred_weights, AGP.pdf.(l, 1, x) .^ 2) - pred[i]^2
                    end
                    return pred, sig_pred
                end
            end

            ### Local Updates Section ###

            function AGP.local_updates!(
                l::$(lname), y::AbstractVector, μ::AbstractVector, diag_cov::AbstractVector
            ) where {T}
                l.c² .=
                    _gen_α.(l, y) - _gen_β.(l, y) .* μ +
                    _gen_γ.(l, y) .* (abs2.(μ) + diag_cov)
                return l.θ .= -_gen_∇φ.(l, l.c²) ./ _gen_φ.(l, l.c²)
            end

            function pω(::$(lname), f)
                #TODO Use the Laplace approx sampler
                @error "You cannot use Gibbs sampling from your likelihood unless you define pω(likelihood,f)"
            end

            function AGP.sample_local!(
                l::$(lname), y::AbstractVector, f::AbstractVector
            ) where {T}
                return set_ω!(
                    l,
                    pω.(
                        l,
                        sqrt.(
                            0.5 * (
                                l,
                                _gen_α.(l, y) - _gen_β.(l, y) .* f +
                                _gen_γ.(l, y) .* (abs2.(f)),
                            ),
                        ),
                    ),
                )
            end

            ### Natural Gradient Section ###

            @inline function AGP.∇E_μ(
                l::$(lname), ::AugmentedGaussianProcesses.AOptimizer, y::AbstractVector
            ) where {T}
                return (_gen_g.(l, y) + l.θ .* _gen_β.(l, y),)
            end
            @inline function AGP.∇E_Σ(
                l::$(lname), ::AugmentedGaussianProcesses.AOptimizer, y::AbstractVector
            ) where {T}
                return (l.θ .* _gen_γ.(l, y),)
            end

            ### ELBO Section ###
            function AGP.expec_loglikelihood(
                l::$(lname),
                i::AnalyticVI,
                y::AbstractVector,
                μ::AbstractVector,
                diag_cov::AbstractVector,
            ) where {T}
                tot = length(y) * log(_gen_C(l))
                tot += dot(_gen_g.(l, y), μ)
                tot +=
                    -(
                        dot(l.θ, _gen_α.(l, y)) - dot(l.θ, _gen_β.(l, y) .* μ) +
                        dot(l.θ, _gen_γ.(l, y) .* (abs2.(μ) + diag_cov))
                    )
                return tot
            end

            function AGP.AugmentedKL(l::$(lname), ::AbstractVector) where {T}
                return -dot(l.c², l.θ) - sum(log, _gen_φ.(l, l.c²))
            end

            ### Gradient Section ###

            @inline function AGP.grad_loglike(
                l::$(lname){T}, y::Real, f::Real
            ) where {T<:Real}
                h² = _gen_α(y) - _gen_β(y) * f + _gen_γ(y) * f^2
                return _gen_g(y) +
                       (-_gen_β(y) + 2 * _gen_γ(y) * f) * _gen_∇φ(h²) / _gen_φ(h²)
            end

            @inline function AGP.hessian_loglike(
                l::$(lname){T}, y::Real, f::Real
            ) where {T<:Real}
                h² = _gen_α(y) - _gen_β(y) * f + _gen_γ(y) * f^2
                φ = _gen_φ(l, h²)
                ∇φ = _gen_∇φ(l, h²)
                ∇²φ = _gen_∇²φ(l, h²)
                return (
                    2 * _gen_γ(y) * ∇φ / φ - ((-_gen_β(y) + 2 * _gen_γ(y) * f) * ∇φ / φ)^2 +
                    (-_gen_β(y) + 2 * _gen_γ(y) * f)^2 * ∇²φ / φ
                )
            end
        end
    end
end


text:
using Polylogarithms
using SpecialFunctions
using Test
using DataFrames, CSV
import Base.MathConstants: π, pi, ℯ, e, γ, eulergamma, catalan, φ, golden
include("test_defs.jl")
include("../src/gamma_derivatives.jl")

@testset "Derivatives of the gamma function at 1.0" begin

    @testset "    throws errors" begin
        @test_throws DomainError D(2, -1, 1.0)
        @test_throws DomainError D(-1, 2, 1.0)
        # @test_throws MethodError stieltjes(1.5)
    end

    @testset "    types" begin
        # @test typeof(stieltjes(1)) == Float64
    end

    @testset "    values" begin
        @test g1(0) ≈ 1.0
        @test g2(0) ≈ 1.0
        @test g1(1) ≈ -γ
        @test g2(1) ≈ -γ
        for n=2:5
            @test g1(n) ≈ g2(n)
        end
    end

end

ext: jl


text:
# ---
# title: 278. First Bad Version
# id: problem278
# author: zhwang
# date: 2022-02-12
# difficulty: Easy
# categories: Binary Search
# link: <https://leetcode.com/problems/first-bad-version/description/>
# hidden: true
# ---
#
# You are a product manager and currently leading a team to develop a new
# product. Unfortunately, the latest version of your product fails the quality
# check. Since each version is developed based on the previous version, all the
# versions after a bad version are also bad.
#
# Suppose you have `n` versions `[1, 2, ..., n]` and you want to find out the
# first bad one, which causes all the following ones to be bad.
#
# You are given an API `bool isBadVersion(version)` which returns whether
# `version` is bad. Implement a function to find the first bad version. You
# should minimize the number of calls to the API.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 5, bad = 4
#     Output: 4
#     Explanation:
#     call isBadVersion(3) -> false
#     call isBadVersion(5) -> true
#     call isBadVersion(4) -> true
#     Then 4 is the first bad version.
#
#
# **Example 2:**
#
#
#
#     Input: n = 1, bad = 1
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= bad <= n <= 231 - 1`
#
#
## @lc code=start
using LeetCode

function first_bad_version(n::Int, is_bad_version::Function)::Int
    left, right = 1, n
    while left <= right
        mid = (left + right) >> 1
        if is_bad_version(mid)
            right = mid -1
        else
            left = mid + 1
        end
    end
    left
end

## @lc code=end

## @lc test=start
@testset "278.first-bad-version.jl" begin
    @test first_bad_version(5, i -> (i>=4)) == 4
    @test first_bad_version(1, i -> (i>=1)) == 1
    @test first_bad_version(2022, i -> (i>=1024)) ==1024
end
## @lc test=end

## @lc info=start
# link: [solution to 278](https://leetcode-cn.com/problems/first-bad-version/solution/278-di-yi-ge-cuo-wu-de-ban-ben-by-rex-rs-92ts/)
## @lc info=end


text:
struct BenchmarkSet
    dynamicalsystem::Function
    n_realizations::Int
    te_func::Function
    tuning_param::Symbol
    tuning_vals::Vector
    which_ts::Vector{Int}
    which_pos::Vector{Int}
    which_lags::Vector{Int}
    vars::TEVars
    ts_lengths::Vector{Int}
    binsizes::Union{Vector{Int}, Vector{Float64}}
end

function summary(b::BenchmarkSet)

end


text:
# ---
# title: 282. Expression Add Operators
# id: problem282
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Divide and Conquer
# link: <https://leetcode.com/problems/expression-add-operators/description/>
# hidden: true
# ---
#
# Given a string that contains only digits `0-9` and a target value, return all
# possibilities to add **binary** operators (not unary) `+`, `-`, or `*` between
# the digits so they evaluate to the target value.
#
# **Example 1:**
#
#
#
#     Input: _num_ = "123", _target_ = 6
#     Output: ["1+2+3", "1*2*3"]
#
#
# **Example 2:**
#
#
#
#     Input: _num_ = "232", _target_ = 8
#     Output: ["2*3+2", "2+3*2"]
#
# **Example 3:**
#
#
#
#     Input: _num_ = "105", _target_ = 5
#     Output: ["1*0+5","10-5"]
#
# **Example 4:**
#
#
#
#     Input: _num_ = "00", _target_ = 0
#     Output: ["0+0", "0-0", "0*0"]
#
#
# **Example 5:**
#
#
#
#     Input: _num_ = "3456237490", _target_ = 9191
#     Output: []
#
#
#
#
# **Constraints:**
#
#   * `0 <= num.length <= 10`
#   * `num` only contain digits.
#
#
## @lc code=start
using LeetCode

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


text:
using Test
using ChowLiuTrees

@testset "ChowLiuTrees Tests" begin
    mi = Float64.([0 1 3 4
          1 0 2 5
          3 2 0 6
          4 5 6 0])

    msts1 = ChowLiuTrees.topk_MST(mi, num_trees=10)

    msts2 = [[(1, 2), (2, 3), (1, 4)],
    [(1, 2), (1, 3), (1, 4)],
    [(1, 2), (2, 3), (2, 4)],
    [(1, 2), (1, 3), (2, 4)],
    [(1, 2), (2, 3), (3, 4)],
    [(2, 3), (1, 3), (1, 4)],
    [(1, 2), (1, 3), (3, 4)],
    [(1, 3), (2, 3), (2, 4)],
    [(2, 3), (1, 4), (2, 4)],
    [(1, 3), (2, 3), (3, 4)]]

    for (a, b) in zip(msts1, msts2)
        @test all(a .== b)
    end
end


text:
##Pakiety

##Zmienne
#Menu z cenami

#Produkty i ich data waznosci

#Początkowy stan produktów

##Rozkłady
#Ilość samochodów na minutę

##Główna funkcja kroku
#jeden krok = jedna minuta

function OneStep()

    #wylosować z rozkładu samochdów na minutę ile samochdów przyjedzie w tym kroku

    #wylosować ilość osób w kazdym samochodzie, 1 samochod = jedna zmienna?

    # wylosować zamówienia dla kazdego samochodu

    #

max_stars_repo_path: Kurczaki i ziemniaki.jl

max_stars_repo_name: myrchas/food-orders-simulation

max_issues_repo_path: Kurczaki i ziemniaki.jl

max_issues_repo_name: myrchas/food-orders-simulation

max_forks_repo_path: Kurczaki i ziemniaki.jl

max_forks_repo_name: myrchas/food-orders-simulation


text:
using DiffEqBase, OrdinaryDiffEq, Catalyst, ReactionNetworkImporters,
      Sundials, Plots, DiffEqDevTools, ODEInterface, ODEInterfaceDiffEq,
      LSODA, TimerOutputs, LinearAlgebra, ModelingToolkit, BenchmarkTools

gr()
datadir  = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr")
const to = TimerOutput()
tf       = 100000.0

# generate ModelingToolkit ODEs
@timeit to "Parse Network" prnbng = loadrxnetwork(BNGNetwork(), joinpath(datadir, "bcr.net"))
rn    = prnbng.rn
@timeit to "Create ODESys" osys = convert(ODESystem, rn)

tspan = (0.,tf)
@timeit to "ODEProb No Jac" oprob = ODEProblem(osys, Float64[], tspan, Float64[])
@timeit to "ODEProb DenseJac" densejacprob = ODEProblem(osys, Float64[], tspan, Float64[], jac=true)

@timeit to "ODEProb SparseJac" sparsejacprob = ODEProblem(osys, Float64[], tspan, Float64[], jac=true, sparse=true)
show(to)

@show numspecies(rn) # Number of ODEs
@show numreactions(rn) # Apprx. number of terms in the ODE
@show length(parameters(rn)) # Number of Parameters

u  = ModelingToolkit.varmap_to_vars(nothing, species(rn); defaults=ModelingToolkit.defaults(rn))
du = copy(u)
p  = ModelingToolkit.varmap_to_vars(nothing, parameters(rn); defaults=ModelingToolkit.defaults(rn))
@timeit to "ODE rhs Eval1" oprob.f(du,u,p,0.)
@timeit to "ODE rhs Eval2" oprob.f(du,u,p,0.)
densejacprob.f(du,u,p,0.)
sparsejacprob.f(du,u,p,0.)

@btime oprob.f($du,$u,$p,0.)

J = zeros(length(u),length(u))
@timeit to "DenseJac Eval1" densejacprob.f.jac(J,u,p,0.)
@timeit to "DenseJac Eval2" densejacprob.f.jac(J,u,p,0.)

Js = similar(sparsejacprob.f.jac_prototype)
@timeit to "SparseJac Eval1" sparsejacprob.f.jac(Js,u,p,0.)
@timeit to "SparseJac Eval2" sparsejacprob.f.jac(Js,u,p,0.)
show(to)

sol = solve(oprob, CVODE_BDF(), saveat=tf/1000., reltol=1e-5, abstol=1e-5)
plot(sol, legend=false, fmt=:png)

@time sol = solve(oprob, CVODE_BDF(), abstol=1/10^12, reltol=1/10^12)
test_sol  = TestSolution(sol)

abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (5:8);
setups = [
          #Dict(:alg=>Rosenbrock23(autodiff=false)),
          Dict(:alg=>TRBDF2(autodiff=false)),
          Dict(:alg=>QNDF(autodiff=false)),
          Dict(:alg=>FBDF(autodiff=false)),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>CVODE_BDF(linear_solver=:LapackDense)),
          #Dict(:alg=>rodas()),
          #Dict(:alg=>radau()),
          #Dict(:alg=>Rodas4(autodiff=false)),
          #Dict(:alg=>Rodas5(autodiff=false)),
          Dict(:alg=>KenCarp4(autodiff=false)),
          Dict(:alg=>KenCarp47(autodiff=false)),
          #Dict(:alg=>RadauIIA5(autodiff=false)),
          #Dict(:alg=>lsoda()),
          ]

wp = WorkPrecisionSet(oprob,abstols,reltols,setups;error_estimate=:l2,
                      saveat=tf/10000.,appxsol=test_sol,maxiters=Int(1e5),numruns=1)
plot(wp)

wp = WorkPrecisionSet(densejacprob,abstols,reltols,setups;error_estimate=:l2,
                      saveat=tf/10000.,appxsol=test_sol,maxiters=Int(1e5),numruns=1)
plot(wp)

setups = [
          #Dict(:alg=>Rosenbrock23(autodiff=false)),
          Dict(:alg=>TRBDF2(autodiff=false)),
          Dict(:alg=>QNDF(autodiff=false)),
          Dict(:alg=>FBDF(autodiff=false)),
          #Dict(:alg=>CVODE_BDF(linear_solver=:KLU)), # Fails!
          #Dict(:alg=>rodas()),
          #Dict(:alg=>radau()),
          #Dict(:alg=>Rodas4(autodiff=false)),
          #Dict(:alg=>Rodas5(autodiff=false)),
          Dict(:alg=>KenCarp4(autodiff=false)),
          Dict(:alg=>KenCarp47(autodiff=false)),
          #Dict(:alg=>RadauIIA5(autodiff=false)),
          #Dict(:alg=>lsoda()),
          ]
wp = WorkPrecisionSet(sparsejacprob,abstols,reltols,setups;error_estimate=:l2,
                      saveat=tf/10000.,appxsol=test_sol,maxiters=Int(1e5),numruns=1)
plot(wp)

using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: SciML/SciMLBenchmarksOutput

max_issues_repo_name: SciML/SciMLBenchmarksOutput

max_forks_repo_name: SciML/SciMLBenchmarksOutput


text:
# matmul.jl: Everything to do with dense matrix multiplication

arithtype(T) = T
arithtype(::Type{Bool}) = Int

# multiply by diagonal matrix as vector
function scale!(C::Matrix, A::Matrix, b::Vector)
    m, n = size(A)
    n==length(b) || throw(DimensionMismatch(""))
    for j = 1:n
        bj = b[j]
        for i = 1:m
            C[i,j] = A[i,j]*bj
        end
    end
    C
end

function scale!(C::Matrix, b::Vector, A::Matrix)
    m, n = size(A)
    m==length(b) || throw(DimensionMismatch(""))
    for j=1:n, i=1:m
        C[i,j] = A[i,j]*b[i]
    end
    C
end
scale(A::Matrix, b::Vector) = scale!(similar(A, promote_type(eltype(A),eltype(b))), A, b)
scale(b::Vector, A::Matrix) = scale!(similar(b, promote_type(eltype(A),eltype(b)), size(A)), b, A)

# Dot products

dot{T<:BlasReal}(x::Vector{T}, y::Vector{T}) = BLAS.dot(x, y)
dot{T<:BlasComplex}(x::Vector{T}, y::Vector{T}) = BLAS.dotc(x, y)
function dot{T<:BlasReal, TI<:Integer}(x::Vector{T}, rx::Union(UnitRange{TI},Range{TI}), y::Vector{T}, ry::Union(UnitRange{TI},Range{TI}))
    length(rx)==length(ry) || throw(DimensionMismatch(""))
    if minimum(rx) < 1 || maximum(rx) > length(x) || minimum(ry) < 1 || maximum(ry) > length(y)
        throw(BoundsError())
    end
    BLAS.dot(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry))
end
function dot{T<:BlasComplex, TI<:Integer}(x::Vector{T}, rx::Union(UnitRange{TI},Range{TI}), y::Vector{T}, ry::Union(UnitRange{TI},Range{TI}))
    length(rx)==length(ry) || throw(DimensionMismatch(""))
    if minimum(rx) < 1 || maximum(rx) > length(x) || minimum(ry) < 1 || maximum(ry) > length(y)
        throw(BoundsError())
    end
    BLAS.dotc(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry))
end
function dot(x::AbstractVector, y::AbstractVector)
    lx = length(x)
    lx==length(y) || throw(DimensionMismatch(""))
    if lx == 0
        return zero(eltype(x))*zero(eltype(y))
    end
    s = conj(x[1])*y[1]
    @inbounds for i = 2:lx
        s += conj(x[i])*y[i]
    end
    s
end
dot(x::Number, y::Number) = conj(x) * y
Ac_mul_B(x::Vector, y::Vector) = [dot(x, y)]
At_mul_B{T<:Real}(x::Vector{T}, y::Vector{T}) = [dot(x, y)]
At_mul_B{T<:BlasComplex}(x::Vector{T}, y::Vector{T}) = [BLAS.dotu(x, y)]

# Matrix-vector multiplication
function (*){T<:BlasFloat,S}(A::StridedMatrix{T}, x::StridedVector{S})
    TS = promote_type(arithtype(T),arithtype(S))
    A_mul_B!(similar(x, TS, size(A,1)), A, convert(AbstractVector{TS}, x))
end
function (*){T,S}(A::AbstractMatrix{T}, x::AbstractVector{S})
    TS = promote_type(arithtype(T),arithtype(S))
    A_mul_B!(similar(x,TS,size(A,1)),A,x)
end
(*)(A::AbstractVector, B::AbstractMatrix) = reshape(A,length(A),1)*B

A_mul_B!{T<:BlasFloat}(y::StridedVector{T}, A::StridedMatrix{T}, x::StridedVector{T}) = gemv!(y, 'N', A, x)
for elty in (Float32,Float64)
    @eval begin
        function A_mul_B!(y::StridedVector{Complex{$elty}}, A::StridedMatrix{Complex{$elty}}, x::StridedVector{$elty})
            Afl = reinterpret($elty,A,(2size(A,1),size(A,2)))
            yfl = reinterpret($elty,y)
            gemv!(yfl,'N',Afl,x)
            return y
        end
    end
end
A_mul_B!(y::StridedVector, A::StridedMatrix, x::StridedVector) = generic_matvecmul!(y, 'N', A, x)

function At_mul_B{T<:BlasFloat,S}(A::StridedMatrix{T}, x::StridedVector{S})
    TS = promote_type(arithtype(T),arithtype(S))
    At_mul_B!(similar(x,TS,size(A,2)), A, convert(AbstractVector{TS}, x))
end
function At_mul_B{T,S}(A::StridedMatrix{T}, x::StridedVector{S})
    TS = promote_type(arithtype(T),arithtype(S))
    At_mul_B!(similar(x,TS,size(A,2)), A, x)
end
At_mul_B!{T<:BlasFloat}(y::StridedVector{T}, A::StridedMatrix{T}, x::StridedVector{T}) = gemv!(y, 'T', A, x)
At_mul_B!(y::StridedVector, A::StridedMatrix, x::StridedVector) = generic_matvecmul!(y, 'T', A, x)

function Ac_mul_B{T<:BlasFloat,S}(A::StridedMatrix{T}, x::StridedVector{S})
    TS = promote_type(arithtype(T),arithtype(S))
    Ac_mul_B!(similar(x,TS,size(A,2)),A,convert(AbstractVector{TS},x))
end
function Ac_mul_B{T,S}(A::StridedMatrix{T}, x::StridedVector{S})
    TS = promote_type(arithtype(T),arithtype(S))
    Ac_mul_B!(similar(x,TS,size(A,2)), A, x)
end

Ac_mul_B!{T<:BlasReal}(y::StridedVector{T}, A::StridedMatrix{T}, x::StridedVector{T}) = At_mul_B!(y, A, x)
Ac_mul_B!{T<:BlasComplex}(y::StridedVector{T}, A::StridedMatrix{T}, x::StridedVector{T}) = gemv!(y, 'C', A, x)
Ac_mul_B!(y::StridedVector, A::StridedMatrix, x::StridedVector) = generic_matvecmul!(y, 'C', A, x)

# Matrix-matrix multiplication

function (*){T,S}(A::StridedMatrix{T}, B::StridedMatrix{S})
    TS = promote_type(arithtype(T),arithtype(S))
    A_mul_B!(similar(B,TS,(size(A,1),size(B,2))),A,B)
end
A_mul_B!{T<:BlasFloat}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = gemm_wrapper!(C, 'N', 'N', A, B)
for elty in (Float32,Float64)
    @eval begin
        function A_mul_B!(C::StridedMatrix{Complex{$elty}}, A::StridedMatrix{Complex{$elty}}, B::StridedMatrix{$elty})
            Afl = reinterpret($elty,A,(2size(A,1),size(A,2)))
            Cfl = reinterpret($elty,C,(2size(C,1),size(C,2)))
            gemm_wrapper!(Cfl,'N','N',Afl,B)
            return C
        end
    end
end
A_mul_B!(C::StridedMatrix, A::StridedMatrix, B::StridedMatrix) = generic_matmatmul!(C, 'N', 'N', A, B)

function At_mul_B{T,S}(A::StridedMatrix{T}, B::StridedMatrix{S})
    TS = promote_type(arithtype(T),arithtype(S))
    At_mul_B!(similar(B,TS,(size(A,2),size(B,2))),A,B)
end
At_mul_B!{T<:BlasFloat}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = is(A,B) ? syrk_wrapper!(C, 'T', A) : gemm_wrapper!(C, 'T', 'N', A, B)
At_mul_B!(C::StridedMatrix, A::StridedMatrix, B::StridedMatrix) = generic_matmatmul!(C, 'T', 'N', A, B)

function A_mul_Bt{T,S}(A::StridedMatrix{T}, B::StridedMatrix{S})
    TS = promote_type(arithtype(T),arithtype(S))
    A_mul_Bt!(similar(B,TS,(size(A,1),size(B,1))),A,B)
end
A_mul_Bt!{T<:BlasFloat}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = is(A,B) ? syrk_wrapper!(C, 'N', A) : gemm_wrapper!(C, 'N', 'T', A, B)
for elty in (Float32,Float64)
    @eval begin
        function A_mul_Bt!(C::StridedMatrix{Complex{$elty}}, A::StridedMatrix{Complex{$elty}}, B::StridedMatrix{$elty})
            Afl = reinterpret($elty,A,(2size(A,1),size(A,2)))
            Cfl = reinterpret($elty,C,(2size(C,1),size(C,2)))
            gemm_wrapper!(Cfl,'N','T',Afl,B)
            return C
        end
    end
end
A_mul_Bt!(C::StridedVecOrMat, A::StridedMatrix, B::StridedMatrix) = generic_matmatmul!(C, 'N', 'T', A, B)

function At_mul_Bt{T,S}(A::StridedMatrix{T}, B::StridedMatrix{S})
    TS = promote_type(arithtype(T),arithtype(S))
    At_mul_Bt!(similar(B,TS,(size(A,2),size(B,1))),A,B)
end
At_mul_Bt!{T<:BlasFloat}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = gemm_wrapper!(C, 'T', 'T', A, B)
At_mul_Bt!(C::StridedMatrix, A::StridedMatrix, B::StridedMatrix) = generic_matmatmul!(C, 'T', 'T', A, B)

Ac_mul_B{T<:BlasReal}(A::StridedMatrix{T}, B::StridedMatrix{T}) = At_mul_B(A, B)
Ac_mul_B!{T<:BlasReal}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = At_mul_B!(C, A, B)
function Ac_mul_B{T,S}(A::StridedMatrix{T}, B::StridedMatrix{S})
    TS = promote_type(arithtype(T),arithtype(S))
    Ac_mul_B!(similar(B,TS,(size(A,2),size(B,2))),A,B)
end
Ac_mul_B!{T<:BlasComplex}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = is(A,B) ? herk_wrapper!(C,'C',A) : gemm_wrapper!(C,'C', 'N', A, B)
Ac_mul_B!(C::StridedMatrix, A::StridedMatrix, B::StridedMatrix) = generic_matmatmul!(C, 'C', 'N', A, B)

A_mul_Bc{T<:BlasFloat,S<:BlasReal}(A::StridedMatrix{T}, B::StridedMatrix{S}) = A_mul_Bt(A, B)
A_mul_Bc!{T<:BlasFloat,S<:BlasReal}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{S}) = A_mul_Bt!(C, A, B)
function A_mul_Bc{T,S}(A::StridedMatrix{T}, B::StridedMatrix{S})
    TS = promote_type(arithtype(T),arithtype(S))
    A_mul_Bc!(similar(B,TS,(size(A,1),size(B,1))),A,B)
end
A_mul_Bc!{T<:BlasComplex}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = is(A,B) ? herk_wrapper!(C, 'N', A) : gemm_wrapper!(C, 'N', 'C', A, B)
A_mul_Bc!(C::StridedMatrix, A::StridedMatrix, B::StridedMatrix) = generic_matmatmul!(C, 'N', 'C', A, B)

Ac_mul_Bc{T,S}(A::StridedMatrix{T}, B::StridedMatrix{S}) = Ac_mul_Bc!(similar(B,promote_type(arithtype(T),arithtype(S)), (size(A,2), size(B,1))), A, B)
Ac_mul_Bc!{T<:BlasFloat}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}) = gemm_wrapper!(C, 'C', 'C', A, B)
Ac_mul_Bc!(C::StridedMatrix, A::StridedMatrix, B::StridedMatrix) = generic_matmatmul!(C, 'C', 'C', A, B)
Ac_mul_Bt{T,S}(A::StridedMatrix{T}, B::StridedMatrix{S}) = Ac_mul_Bt(similar(B, promote_type(arithtype(A),arithtype(B)), (size(A,2), size(B,1))), A, B)
Ac_mul_Bt!(C::StridedMatrix, A::StridedMatrix, B::StridedMatrix) = generic_matmatmul!(C, 'C', 'T', A, B)

# Supporting functions for matrix multiplication

function copytri!(A::StridedMatrix, uplo::Char, conjugate::Bool=false)
    n = chksquare(A)
    @chkuplo
    if uplo == 'U'
        for i = 1:(n-1), j = (i+1):n
            A[j,i] = conjugate ? conj(A[i,j]) : A[i,j]
        end
    elseif uplo == 'L'
        for i = 1:(n-1), j = (i+1):n
            A[i,j] = conjugate ? conj(A[j,i]) : A[j,i]
        end
    else
        throw(ArgumentError("second argument must be 'U' or 'L'"))
    end
    A
end

function gemv!{T<:BlasFloat}(y::StridedVector{T}, tA::Char, A::StridedMatrix{T}, x::StridedVector{T})
    stride(A, 1)==1 || return generic_matvecmul!(y, tA, A, x)

    if tA != 'N'
        (nA, mA) = size(A)
    else
        (mA, nA) = size(A)
    end

    nA==length(x) || throw(DimensionMismatch(""))
    mA==length(y) || throw(DimensionMismatch(""))
    mA == 0 && return zeros(T, 0)
    nA == 0 && return zeros(T, mA)
    return BLAS.gemv!(tA, one(T), A, x, zero(T), y)
end

function syrk_wrapper!{T<:BlasFloat}(C::StridedMatrix{T}, tA::Char, A::StridedMatrix{T})
    nC = chksquare(C)
    if tA == 'T'
        (nA, mA) = size(A)
        tAt = 'N'
    else
        (mA, nA) = size(A)
        tAt = 'T'
    end
    nC == mA || throw(DimensionMismatch("output matrix has size: $(nC), but should have size $(mA)"))
    if mA == 0 || nA == 0; return C; end
    if mA == 2 && nA == 2; return matmul2x2!(C,tA,tAt,A,A); end
    if mA == 3 && nA == 3; return matmul3x3!(C,tA,tAt,A,A); end

    stride(A, 1) == 1 || (return generic_matmatmul!(C, tA, tAt, A, A))
    copytri!(BLAS.syrk!('U', tA, one(T), A, zero(T), C), 'U')
end

function herk_wrapper!{T<:BlasFloat}(C::StridedMatrix{T}, tA::Char, A::StridedMatrix{T})
    nC = chksquare(C)
    if tA == 'C'
        (nA, mA) = size(A)
        tAt = 'N'
    else
        (mA, nA) = size(A)
        tAt = 'C'
    end
    nC == mA || throw(DimensionMismatch("output matrix has size: $(nC), but should have size $(mA)"))
    if mA == 0 || nA == 0; return C; end
    if mA == 2 && nA == 2; return matmul2x2!(C,tA,tAt,A,A); end
    if mA == 3 && nA == 3; return matmul3x3!(C,tA,tAt,A,A); end

    stride(A, 1) == 1 || (return generic_matmatmul!(C,tA, tAt, A, A))

    # Result array does not need to be initialized as long as beta==0
    #    C = Array(T, mA, mA)

    copytri!(BLAS.herk!('U', tA, one(T), A, zero(T), C), 'U', true)
end

function gemm_wrapper{T<:BlasFloat}(tA::Char, tB::Char,
                             A::StridedVecOrMat{T},
                             B::StridedMatrix{T})
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
    C = similar(B, T, mA, nB)
    gemm_wrapper!(C, tA, tB, A, B)
end

function gemm_wrapper!{T<:BlasFloat}(C::StridedVecOrMat{T}, tA::Char, tB::Char,
                             A::StridedVecOrMat{T},
                             B::StridedMatrix{T})
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)

    nA==mB || throw(DimensionMismatch("*"))

    if mA == 0 || nA == 0 || nB == 0; return zeros(T, mA, nB); end
    if mA == 2 && nA == 2 && nB == 2; return matmul2x2!(C,tA,tB,A,B); end
    if mA == 3 && nA == 3 && nB == 3; return matmul3x3!(C,tA,tB,A,B); end

    stride(A, 1)==stride(B, 1)==1 || (return generic_matmatmul!(C, tA, tB, A, B))
    BLAS.gemm!(tA, tB, one(T), A, B, zero(T), C)
end

# blas.jl defines matmul for floats; other integer and mixed precision
# cases are handled here

lapack_size(t::Char, M::AbstractVecOrMat) = (size(M, t=='N' ? 1:2), size(M, t=='N' ? 2:1))

function copy!{R,S}(B::AbstractMatrix{R}, ir_dest::UnitRange{Int}, jr_dest::UnitRange{Int}, tM::Char, M::AbstractMatrix{S}, ir_src::UnitRange{Int}, jr_src::UnitRange{Int})
    if tM == 'N'
        copy!(B, ir_dest, jr_dest, M, ir_src, jr_src)
    else
        Base.copy_transpose!(B, ir_dest, jr_dest, M, jr_src, ir_src)
        tM == 'C' && conj!(B)
    end
end

function copy_transpose!{R,S}(B::AbstractMatrix{R}, ir_dest::UnitRange{Int}, jr_dest::UnitRange{Int}, tM::Char, M::AbstractVecOrMat{S}, ir_src::UnitRange{Int}, jr_src::UnitRange{Int})
    if tM == 'N'
        Base.copy_transpose!(B, ir_dest, jr_dest, M, ir_src, jr_src)
    else
        copy!(B, ir_dest, jr_dest, M, jr_src, ir_src)
        tM == 'C' && conj!(B)
    end
end

# TODO: It will be faster for large matrices to convert to float,
# call BLAS, and convert back to required type.

# NOTE: the generic version is also called as fallback for
#       strides != 1 cases

function generic_matvecmul!{T,S,R}(C::AbstractVector{R}, tA, A::AbstractMatrix{T}, B::AbstractVector{S})
    mB = length(B)
    mA, nA = lapack_size(tA, A)
    mB==nA || throw(DimensionMismatch("*"))
    mA==length(C) || throw(DimensionMismatch("*"))
    z = zero(R)

    Astride = size(A, 1)

    if tA == 'T'  # fastest case
        for k = 1:mA
            aoffs = (k-1)*Astride
            s = z
            for i = 1:nA
                s += A[aoffs+i] * B[i]
            end
            C[k] = s
        end
    elseif tA == 'C'
        for k = 1:mA
            aoffs = (k-1)*Astride
            s = z
            for i = 1:nA
                s += conj(A[aoffs+i]) * B[i]
            end
            C[k] = s
        end
        else # tA == 'N'
        fill!(C, z)
        for k = 1:mB
            aoffs = (k-1)*Astride
            b = B[k]
            for i = 1:mA
                C[i] += A[aoffs+i] * b
            end
        end
    end
    C
end

function generic_matmatmul{T,S}(tA, tB, A::AbstractVecOrMat{T}, B::AbstractMatrix{S})
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
    C = similar(B, promote_type(arithtype(T),arithtype(S)), mA, nB)
    generic_matmatmul!(C, tA, tB, A, B)
end

const tilebufsize = 10800  # Approximately 32k/3
const Abuf = Array(Uint8, tilebufsize)
const Bbuf = Array(Uint8, tilebufsize)
const Cbuf = Array(Uint8, tilebufsize)

function generic_matmatmul!{T,S,R}(C::AbstractVecOrMat{R}, tA, tB, A::AbstractVecOrMat{T}, B::AbstractMatrix{S})
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
    mB==nA || throw(DimensionMismatch("*"))
    if size(C,1) != mA || size(C,2) != nB; throw(DimensionMismatch("*")); end

    if mA == nA == nB == 2; return matmul2x2!(C, tA, tB, A, B); end
    if mA == nA == nB == 3; return matmul3x3!(C, tA, tB, A, B); end

    @inbounds begin
    if isbits(R)
        tile_size = int(ifloor(sqrt(tilebufsize/sizeof(R))))
        sz = (tile_size, tile_size)
        Atile = pointer_to_array(convert(Ptr{R}, pointer(Abuf)), sz)
        Btile = pointer_to_array(convert(Ptr{R}, pointer(Bbuf)), sz)

        z = zero(R)

        if mA < tile_size && nA < tile_size && nB < tile_size
            Base.copy_transpose!(Atile, 1:nA, 1:mA, tA, A, 1:mA, 1:nA)
            copy!(Btile, 1:mB, 1:nB, tB, B, 1:mB, 1:nB)
            for j = 1:nB
                boff = (j-1)*tile_size
                for i = 1:mA
                    aoff = (i-1)*tile_size
                    s = z
                    for k = 1:nA
                        s += Atile[aoff+k] * Btile[boff+k]
                    end
                    C[i,j] = s
                end
            end
        else
            Ctile = pointer_to_array(convert(Ptr{R}, pointer(Cbuf)), sz)
            for jb = 1:tile_size:nB
                jlim = min(jb+tile_size-1,nB)
                jlen = jlim-jb+1
                for ib = 1:tile_size:mA
                    ilim = min(ib+tile_size-1,mA)
                    ilen = ilim-ib+1
                    fill!(Ctile, z)
                    for kb = 1:tile_size:nA
                        klim = min(kb+tile_size-1,mB)
                        klen = klim-kb+1
                        Base.copy_transpose!(Atile, 1:klen, 1:ilen, tA, A, ib:ilim, kb:klim)
                        copy!(Btile, 1:klen, 1:jlen, tB, B, kb:klim, jb:jlim)
                        for j=1:jlen
                            bcoff = (j-1)*tile_size
                            for i = 1:ilen
                                aoff = (i-1)*tile_size
                                s = z
                                for k = 1:klen
                                    s += Atile[aoff+k] * Btile[bcoff+k]
                                end
                                Ctile[bcoff+i] += s
                            end
                        end
                    end
                    copy!(C, ib:ilim, jb:jlim, Ctile, 1:ilen, 1:jlen)
                end
            end
        end
    else
        # Multiplication for non-plain-data uses the naive algorithm
        if tA == 'N'
            if tB == 'N'
                for i = 1:mA, j = 1:nB
                    Ctmp = A[i, 1]*B[1, j]
                    for k = 2:nA
                        Ctmp += A[i, k]*B[k, j]
                    end
                    C[i,j] = Ctmp
                end
            elseif tB == 'T'
                for i = 1:mA, j = 1:nB
                    Ctmp = A[i, 1]*B[j, 1]
                    for k = 2:nA
                        Ctmp += A[i, k]*B[j, k]
                    end
                    C[i,j] = Ctmp
                end
            else
                for i = 1:mA, j = 1:nB
                    Ctmp = A[i, 1]*conj(B[j, 1])
                    for k = 2:nA
                        Ctmp += A[i, k]*conj(B[j, k])
                    end
                    C[i,j] = Ctmp
                end
            end
        elseif tA == 'T'
            if tB == 'N'
                for i = 1:mA, j = 1:nB
                    Ctmp = A[1, i]*B[1, j]
                    for k = 2:nA
                        Ctmp += A[k, i]*B[k, j]
                    end
                    C[i,j] = Ctmp
                end
            elseif tB == 'T'
                for i = 1:mA, j = 1:nB
                    Ctmp = A[1, i]*B[j, 1]
                    for k = 2:nA
                        Ctmp += A[k, i]*B[j, k]
                    end
                    C[i,j] = Ctmp
                end
            else
                for i = 1:mA, j = 1:nB
                    Ctmp = A[1, i]*conj(B[j, 1])
                    for k = 2:nA
                        Ctmp += A[k, i]*conj(B[j, k])
                    end
                    C[i,j] = Ctmp
                end
            end
        else
            if tB == 'N'
                for i = 1:mA, j = 1:nB
                    Ctmp = conj(A[1, i])*B[1, j]
                    for k = 2:nA
                        Ctmp += conj(A[k, i])*B[k, j]
                    end
                    C[i,j] = Ctmp
                end
            elseif tB == 'T'
                for i = 1:mA, j = 1:nB
                    Ctmp = conj(A[1, i])*B[j, 1]
                    for k = 2:nA
                        Ctmp += conj(A[k, i])*B[j, k]
                    end
                    C[i,j] = Ctmp
                end
            else
                for i = 1:mA, j = 1:nB
                    Ctmp = conj(A[1, i])*conj(B[j, 1])
                    for k = 2:nA
                        Ctmp += conj(A[k, i])*conj(B[j, k])
                    end
                    C[i,j] = Ctmp
                end
            end
        end
    end
    end # @inbounds
    C
end

# multiply 2x2 matrices
function matmul2x2{T,S}(tA, tB, A::AbstractMatrix{T}, B::AbstractMatrix{S})
    matmul2x2!(similar(B, promote_type(T,S), 2, 2), tA, tB, A, B)
end

function matmul2x2!{T,S,R}(C::AbstractMatrix{R}, tA, tB, A::AbstractMatrix{T}, B::AbstractMatrix{S})
    if tA == 'T'
        A11 = A[1,1]; A12 = A[2,1]; A21 = A[1,2]; A22 = A[2,2]
    elseif tA == 'C'
        A11 = conj(A[1,1]); A12 = conj(A[2,1]); A21 = conj(A[1,2]); A22 = conj(A[2,2])
    else
        A11 = A[1,1]; A12 = A[1,2]; A21 = A[2,1]; A22 = A[2,2]
    end
    if tB == 'T'
        B11 = B[1,1]; B12 = B[2,1]; B21 = B[1,2]; B22 = B[2,2]
    elseif tB == 'C'
        B11 = conj(B[1,1]); B12 = conj(B[2,1]); B21 = conj(B[1,2]); B22 = conj(B[2,2])
    else
        B11 = B[1,1]; B12 = B[1,2]; B21 = B[2,1]; B22 = B[2,2]
    end
    C[1,1] = A11*B11 + A12*B21
    C[1,2] = A11*B12 + A12*B22
    C[2,1] = A21*B11 + A22*B21
    C[2,2] = A21*B12 + A22*B22
    C
end

# Multiply 3x3 matrices
function matmul3x3{T,S}(tA, tB, A::AbstractMatrix{T}, B::AbstractMatrix{S})
    matmul3x3!(similar(B, promote_type(T,S), 3, 3), tA, tB, A, B)
end

function matmul3x3!{T,S,R}(C::AbstractMatrix{R}, tA, tB, A::AbstractMatrix{T}, B::AbstractMatrix{S})
    if tA == 'T'
        A11 = A[1,1]; A12 = A[2,1]; A13 = A[3,1];
        A21 = A[1,2]; A22 = A[2,2]; A23 = A[3,2];
        A31 = A[1,3]; A32 = A[2,3]; A33 = A[3,3];
    elseif tA == 'C'
        A11 = conj(A[1,1]); A12 = conj(A[2,1]); A13 = conj(A[3,1]);
        A21 = conj(A[1,2]); A22 = conj(A[2,2]); A23 = conj(A[3,2]);
        A31 = conj(A[1,3]); A32 = conj(A[2,3]); A33 = conj(A[3,3]);
    else
        A11 = A[1,1]; A12 = A[1,2]; A13 = A[1,3];
        A21 = A[2,1]; A22 = A[2,2]; A23 = A[2,3];
        A31 = A[3,1]; A32 = A[3,2]; A33 = A[3,3];
    end

    if tB == 'T'
        B11 = B[1,1]; B12 = B[2,1]; B13 = B[3,1];
        B21 = B[1,2]; B22 = B[2,2]; B23 = B[3,2];
        B31 = B[1,3]; B32 = B[2,3]; B33 = B[3,3];
    elseif tB == 'C'
        B11 = conj(B[1,1]); B12 = conj(B[2,1]); B13 = conj(B[3,1]);
        B21 = conj(B[1,2]); B22 = conj(B[2,2]); B23 = conj(B[3,2]);
        B31 = conj(B[1,3]); B32 = conj(B[2,3]); B33 = conj(B[3,3]);
    else
        B11 = B[1,1]; B12 = B[1,2]; B13 = B[1,3];
        B21 = B[2,1]; B22 = B[2,2]; B23 = B[2,3];
        B31 = B[3,1]; B32 = B[3,2]; B33 = B[3,3];
    end

    C[1,1] = A11*B11 + A12*B21 + A13*B31
    C[1,2] = A11*B12 + A12*B22 + A13*B32
    C[1,3] = A11*B13 + A12*B23 + A13*B33

    C[2,1] = A21*B11 + A22*B21 + A23*B31
    C[2,2] = A21*B12 + A22*B22 + A23*B32
    C[2,3] = A21*B13 + A22*B23 + A23*B33

    C[3,1] = A31*B11 + A32*B21 + A33*B31
    C[3,2] = A31*B12 + A32*B22 + A33*B32
    C[3,3] = A31*B13 + A32*B23 + A33*B33

    C
end

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
"""
    Computation{R, T, O, G}

Represents an operation performed over the elements of a field.
"""
struct Computation{R, T, O, G}
      operation :: O
         result :: R
           grid :: G
    return_type :: T
end

"""
    Computation(operation, result; return_type=Array)

Returns a `Computation` representing an `operation` performed over the elements of
`operation.grid` and stored in `result`. `return_type` specifies the output type when the
`Computation` instances is called as a function.
"""
Computation(operation, result; return_type=Array) =
    Computation(operation, result, operation.grid, return_type)

architecture(comp::Computation) = architecture(comp.result)
Base.parent(comp::Computation) = comp # this enables taking a "horizontal average" of a computation

"""
    compute!(computation::Computation)

Perform a `computation`. The result is stored in `computation.result`.
"""
function compute!(computation::Computation)
    arch = architecture(computation.result)
    result_data = data(computation.result)
    @launch(device(arch), config=launch_config(computation.grid, :xyz),
            _compute!(result_data, computation.grid, computation.operation))
    return nothing
end

"""Compute an `operation` over `grid` and store in `result`."""
function _compute!(result, grid, operation)
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds result[i, j, k] = operation[i, j, k]
            end
        end
    end
    return nothing
end

"""
    (computation::Computation)(args...)

Performs the `compute(computation)` and returns the result if `isnothing(return_type)`,
or the result after being converted to `return_type`.
"""
function (computation::Computation)(args...)
    compute!(computation)
    return computation.return_type(interior(computation.result))
end

function (computation::Computation{<:Nothing})(args...)
    compute!(computation)
    return computation.result
end

#####
##### Functionality for using computations with HorizontalAverage
#####

"""
    HorizontalAverage(op::AbstractOperation, result; kwargs...)

Returns the representation of a `HorizontalAverage` over the operation `op`, using
`result` as a temporary array to store the result of `operation` computed on `op.grid`.
"""
function HorizontalAverage(op::AbstractOperation, result; kwargs...)
    computation = Computation(op, result)
    return HorizontalAverage(computation; kwargs...)
end

"""
    HorizontalAverage(op::AbstractOperation, model; kwargs...)

Returns the representation of a `HorizontalAverage` over the operation `op`, using
`model.pressures.pHY′` as a temporary array to store the result of `operation` computed on
`op.grid`.
"""
HorizontalAverage(op::AbstractOperation, model::AbstractModel; kwargs...) =
    HorizontalAverage(op, model.pressures.pHY′; kwargs...)

"""Compute the horizontal average of a computation."""
function run_diagnostic(model, havg::HorizontalAverage{<:Computation})
    compute!(havg.field)
    zero_halo_regions!(parent(havg.field.result), model.grid)
    sum!(havg.result, parent(havg.field.result))
    normalize_horizontal_sum!(havg, model.grid)
    return nothing
end


text:
module TestDecisionTree

using Test
using TSML
using Random
using DataFrames: nrow, DataFrame
using Statistics

function generateXY()
    Random.seed!(123)
    iris = getiris()
    indx = Random.shuffle(1:nrow(iris))
    features=iris[indx,1:4]
    sp = iris[indx,5] |> Vector
    z = rand(nrow(iris))
    (features,sp,z)
end

rmse(x,y)= sqrt(mean((x .- y).^2))
acc(x,y)=score(:accuracy,x,y)

function test_libsvm()
   X,Y,Z = generateXY()

   @test_throws ArgumentError SVMModel("SVCP"; gamma = :auto)
   linear=SVMModel("LinearSVC"; gamma = :auto)
   @test_throws MethodError fit!(linear,X,Y)
   @test_throws MethodError fit(linear,X,Y)

   model = SVMModel("OneClassSVM")
   @test fit_transform!(model,X,Y) |> sum == 133
   @test fit_transform(model,X,Y) |> sum == 133

   learner1 = SVMModel("SVC"; gamma = :auto)
   learner2 = SVMModel("SVC",Dict(:gamma => :auto))

   @test learner1.model[:impl_args] == learner2.model[:impl_args]
   for learner in ["NuSVC","SVC","LinearSVC"]
      model = SVMModel(learner)
      @test crossvalidate(model,X,Y,acc,10,false).mean > 90.0
   end

   for learner in ["NuSVR","EpsilonSVR"]
      model = SVMModel(learner, tolerance = 0.01)
      @test crossvalidate(model,X,Z,rmse,10,false).mean < 0.3
   end
end
@testset "LIBSVM" begin
   Random.seed!(123)
   test_libsvm()
end

end


text:
using HyperdimensionalComputing

# get the different types
println("The time to encode a 4-letter alphabet as bipolar or binary hypervectors, respectively:")
@time bip = encode_alphabet(['A', 'C', 'T', 'G'])
@time bin = encode_alphabet(['A', 'C', 'T', 'G'], vectortype="binary")
println()

# benchmark multiplication
println("The time to combine two bipolar or binary hypervectors, respectively, into a new one that is highly dissimilar:")
@time multiply(bip['A'], bip['C'])
@time multiply(bin['A'], bin['C'])
println()

# benchmark aggregation
println("The time to combine two bipolar or binary hypervectors, respectively, into a new one that is highly similar:")
@time add(bip['A'], bip['C'])
@time add(bin['A'], bin['C'])
println()

# benchmark rotation
println("The time to rotate/permute bipolar or binary hypervectors, respectively:")
@time rotate(bip['A'], 1)
@time rotate(bin['A'], 1)
println()

# benchmark sequence encoding
seq = "ATCGATAGCA"
println("The time to encode a 10-letter sequence as bipolar or binary hypervectors, respectively:")
@time encode_sequence(seq, bip, k=3)
@time encode_sequence(seq, bin, k=3)
println()


text:
__precompile__()

module ExprOptimization

export
        @grammar,
        Grammar,
        RuleNode,
        get_executable,

        optimize,
        ExprOptAlgorithm,
        ExprOptResult,
        get_expr,

        ProbabilisticExprRules,
        PPTs,
        PPT,

        MonteCarlos,
        MonteCarlo,

        GeneticPrograms,
        GeneticProgram,

        GrammaticalEvolutions,
        GrammaticalEvolution,

        CrossEntropys,
        CrossEntropy,

        PIPEs,
        PIPE

using Reexport
@reexport using ExprRules

abstract type ExprOptAlgorithm end

"""
    ExprOptResult

Returned by optimize().  Contains the results of the optimization.
"""
struct ExprOptResult
    tree::RuleNode #best tree
    loss::Float64 #best loss
    expr::Any #best expression
    alg_result::Any #algorithm-specific results
end

"""
    optimize(p::ExprOptAlgorithm, grammar::Grammar, typ::Symbol, loss::Function; kwargs...)

Main entry for expression optimization.  Use concrete ExprOptAlgorithm to specify optimization algorithm. Optimize using grammar and start symbol, typ, and loss function.  Loss function has the form: los::Float64=loss(node::RuleNode).
"""
function optimize end   #implemented by algorithms

"""
    get_expr(result::ExprOptResult)

Returns the expression in the result
"""
get_expr(result::ExprOptResult) = result.expr
get_expr(x::Nothing) = nothing

#############################################################################
# Common base modules
include("../contrib/BoundedPriorityQueues.jl")
using .BoundedPriorityQueues

include("ProbabilisticExprRules/ProbabilisticExprRules.jl")
include("PPT/PPTs.jl")
using .PPTs: PPT

#############################################################################
# Optimization algorithms

include("MonteCarlo/MonteCarlos.jl")
using .MonteCarlos

include("GeneticProgram/GeneticPrograms.jl")
using .GeneticPrograms

include("GrammaticalEvolution/GrammaticalEvolutions.jl")
using .GrammaticalEvolutions

include("CrossEntropy/CrossEntropys.jl")
using .CrossEntropys

include("PIPE/PIPEs.jl")
using .PIPEs

end # module


text:
export Expression, Constant, Numeric, Im, Term
export @vars, Variable
export @term

abstract type Expression end
struct Variable <: Expression
    name::Symbol
end

struct Numeric{T} <: Expression
    value::T
end

struct Constant{sym} <: Expression end
struct Im <: Expression end

Constant(sym::Symbol) = Constant{sym}()
Constant(x::Irrational{sym}) where sym = Constant{sym}()

struct Term{F} <: Expression
    head::F
    args::Vector
end

Term(f::F, xs...) where F = Term(f, collect(xs))

Base.length(t::Term) = length(t.args)
Base.length(::Variable) = 1
Base.length(::Constant) = 1
Base.length(::Im) = 1

Base.:(==)(::Term, ::Term) = false
Base.:(==)(t1::Term{F}, t2::Term{F}) where F = t1.args == t2.args

Base.:(==)(x::Number, y::Numeric) = x == y.value
Base.:(==)(x::Numeric, y::Number) = x.value == y
Base.:(==)(x::Numeric, y::Numeric) = x.value == y.value

Base.Int(x::Numeric) = Int(x.value)

Base.convert(::Type{<:Expression}, x::Number) = Numeric(x)
Base.convert(::Type{<:Expression}, ::Irrational{sym}) where sym = Constant(sym)

macro term(ex)
    term_m(ex)
end
to_expr(x) = x
to_expr(x::Number) = Numeric(x)
to_expr(x::Irrational) = Constant(x)

term_m(x) = :(to_expr($(esc(x))))

function term_m(ex::Expr)
    ex.head === :call || throw(Meta.ParseError("expect function call"))
    :(Term($(term_m.(ex.args)...)))
end

# macro constructors
function varm(xs::Symbol...)
    ts = Expr(:tuple)
    vs = Expr(:tuple)
    for x in xs
        if x in RESERVED_TOKENS
            throw(Meta.ParseError("$x is reserved, use another variable name"))
        end
        push!(vs.args, x)
        push!(ts.args, :(Variable($(Meta.quot(x)))))
    end
    return :($vs = $ts)
end

macro vars(xs::Symbol...)
    return esc(varm(xs...))
end

macro vars(xs...)
    quote
        throw(Meta.ParseError("expect symbols got $xs"))
    end
end


text:
#=
composite_test:
- Julia version:
- Author: Dan
- Date: 2021-06-04
=#
#---------------TEST----------------------
include("../src/CONCEPT/composite.jl")
using Test
@testset "myset" begin
    @testset "under1" begin
        a = [1,2,3]
        b = [1,2,5]
        @test a!=b
    end
    #will open up a new scope and fail.
    #@testset "under2" begin
    #    @test a == [1,2,3]
    #end
end

a_element = Atom()
a_element.name_ = "a"
b_element = Atom()
b_element.name_ = "b"
c_element = Atom()
c_element.name_ = "c"
d_element = Atom()
d_element.name_ = "d"
e_element = Atom()
e_element.name_ = "e"
x_element = Atom()
x_element.name_ = "x"

a = Composite()
b = Composite()
c = Composite()
d = Composite()
e = Composite()
x = Composite()

a.trait_ = a_element
b.trait_ = b_element
c.trait_ = c_element
d.trait_ = d_element
e.trait_ = e_element
x.trait_ = x_element

a.first_child_ = x
a.last_child_ = b

x.first_child_ = missing
x.next_ = b

b.first_child_ = c
b.last_child_ = d

c.first_child_ = e
c.last_child_ = e
c.next_ = d

e.first_child_ = missing
d.first_child_ = missing

for item in a
    #println(item)
end


text:
#!/usr/bin/env julia

# This script generates the table published at https://julialang.org/benchmarks/
# (file _includes/benchmarks.html in the JuliaLang/julialang.github.com repository)

using Compat
import Compat.Statistics
import Compat.Printf

const benchmark_order = [
    "iteration_pi_sum",
    "recursion_fibonacci",
    "recursion_quicksort",
    "parse_integers",
    "print_to_file",
    "matrix_statistics",
    "matrix_multiply",
    "userfunc_mandelbrot",
]

const versions = Dict{String, String}()
const benchmarks = Dict{String, Dict{String, Float64}}()

# Read versions.csv
for line in eachline(ARGS[1])
    lang, version = split(line, ',')
    versions[lang] = version
end

# Read benchmarks.csv
for line in eachline(ARGS[2])
    lang, bench, time_str = split(line, ',')
    times = get!(benchmarks, lang, Dict{String, Float64}())
    times[bench] = parse(Float64, time_str)
end

const labels = Dict{String, String}(
    "c"           => "C"                ,
    "julia"       => "Julia"            ,
    "lua"         => "LuaJIT"           ,
    "fortran"     => "Fortran"          ,
    "java"        => "Java"             ,
    "javascript"  => "JavaScript"       ,
    "matlab"      => "Matlab"           ,
    "python"      => "Python"           ,
    "mathematica" => "Mathematica" ,
    "r"           => "R"                ,
    "octave"      => "Octave"           ,
    "go"          => "Go"               ,
    "rust"        => "Rust"             ,
)

# Produce the sorting order for the list of languages
function lang_by(lang::String)
    # C is placed at the start of the list
    lang == "c" ? -Inf :
    # Julia is sorted immediately after C
    @compat lang == "julia" ? -floatmax() :
    # The rest of the languages are sorted by the geometric mean of their benchmark values
    # See https://en.wikipedia.org/wiki/Geometric_mean#Relationship_with_logarithms for details
    @compat exp(Statistics.mean(log.(collect(values(benchmarks[lang])))))
end

const language_order = sort!(collect(keys(benchmarks)), by=lang_by)

print("""

    <colgroup>
        <col class="name">
        <col class="relative" span="$(length(language_order))">
    </colgroup>

""")

for lang in language_order
    label = get(labels, lang, lang)
    println("            $label	")
end

print("""

""")

for lang in language_order
    version = get(versions, lang, "???")
    println("            $version	")
end

print("""

""")

for benchmark in benchmark_order
    println("
$benchmark	")
    c_time = benchmarks["c"][benchmark]
    for lang in language_order
        rel_time = "n/a"
        if haskey(benchmarks[lang], benchmark)
            @compat rel_time = Printf.@sprintf "%.2f" benchmarks[lang][benchmark]/c_time
        end
        println("            $rel_time	")
    end
    println("
")
end

print("""

""")

max_stars_repo_name: MilesCranmer/Microbenchmarks

max_issues_repo_name: MilesCranmer/Microbenchmarks

max_forks_repo_name: MilesCranmer/Microbenchmarks


text:
# ---
# title: 1675. Minimize Deviation in Array
# id: problem1675
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Heap, Ordered Map
# link: <https://leetcode.com/problems/minimize-deviation-in-array/description/>
# hidden: true
# ---
#
# You are given an array `nums` of `n` positive integers.
#
# You can perform two types of operations on any element of the array any number
# of times:
#
#   * If the element is **even** , **divide** it by `2`.
#     * For example, if the array is `[1,2,3,4]`, then you can do this operation on the last element, and the array will be `[1,2,3, _2_ ].`
#   * If the element is **odd** , **multiply** it by `2`.
#     * For example, if the array is `[1,2,3,4]`, then you can do this operation on the first element, and the array will be `[ _2_ ,2,3,4].`
#
# The **deviation** of the array is the **maximum difference** between any two
# elements in the array.
#
# Return _the **minimum deviation** the array can have after performing some
# number of operations._
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3,4]
#     Output: 1
#     Explanation: You can transform the array to [1,2,3, _2_ ], then to [ _2_ ,2,3,2], then the deviation will be 3 - 2 = 1.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [4,1,5,20,3]
#     Output: 3
#     Explanation: You can transform the array after two operations to [4, _2_ ,5, _5_ ,3], then the deviation will be 5 - 2 = 3.
#
#
# **Example 3:**
#
#
#
#     Input: nums = [2,10,8]
#     Output: 3
#
#
#
#
# **Constraints:**
#
#   * `n == nums.length`
#   * `2 <= n <= 105`
#   * `1 <= nums[i] <= 109`
#
#
## @lc code=start
using LeetCode

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


text:
module MicrobiomePlots

using Microbiome
using RecipesBase
using Colors
using Reexport
using Random
@reexport using StatsPlots

export
    abundanceplot,
    AnnotationBar,
    annotationbar
    # zeroyplot

include("recipes.jl")

end

max_issues_repo_name: BioJulia/MicrobiomePlots

max_forks_repo_name: BioJulia/MicrobiomePlots


text:
export makeTagType,
       op,
       state

function makeTagType(t)
  tag = Tag(t)
  return Val{tag}
end

function _call_op(s::Index,
                  opname::AbstractString;
                  kwargs...)
  use_tag = 0
  nfound = 0
  for n=1:length(tags(s))
    TType = Val{tags(s)[n]}
    if hasmethod(op,Tuple{TType,Index,AbstractString})
      use_tag = n
      nfound += 1
    end
  end
  if nfound == 0
    error("Overload of \"op\" function not found for Index tags $(tags(s))")
  elseif nfound > 1
    error("Multiple tags from $(tags(s)) overload the function \"op\"")
  end
  TType = Val{tags(s)[use_tag]}
  return op(TType(),s,opname;kwargs...)
end

function op(s::Index,
            opname::AbstractString;
            kwargs...)::ITensor

  opname = strip(opname)

  if opname == "Id"
    Op = ITensor(dag(s),s')
    for n=1:dim(s)
      Op[dag(s)(n),s'(n)] = 1.0
    end
    return Op
  end

  # Interpret operator names joined by *
  # as acting sequentially on the same site
  starpos = findfirst("*",opname)
  if !isnothing(starpos)
    op1 = opname[1:starpos.start-1]
    op2 = opname[starpos.start+1:end]
    return multSiteOps(op(s,op1;kwargs...),op(s,op2;kwargs...))
  end

  return _call_op(s,opname;kwargs...)
end

function op(s::Vector{Index},
            opname::AbstractString,
            n::Int;
            kwargs...)::ITensor
  return op(s[n],opname;kwargs...)
end

state(s::Index,n::Integer) = s[n]

function state(s::Index,
               str::String)::IndexVal
  use_tag = 0
  nfound = 0
  for n=1:length(tags(s))
    TType = Val{tags(s)[n]}
    if hasmethod(state,Tuple{TType,AbstractString})
      use_tag = n
      nfound += 1
    end
  end
  if nfound == 0
    error("Overload of \"state\" function not found for Index tags $(tags(s))")
  elseif nfound > 1
    error("Multiple tags from $(tags(s)) overload the function \"state\"")
  end
  TType = Val{tags(s)[use_tag]}
  sn = state(TType(),str)
  return s[sn]
end

function state(sset::Vector{Index},
               j::Integer,
               st)::IndexVal
  return state(sset[j],st)
end


text:
abstract type AbstractAggregateTools end

const ElIndices = Vector{Vector{T}} where {T<:Integer}
const VibIndices = Vector{Vector{Vector{Int64}}} where {T<:Integer}
# TODO: be more specific for T
const Indices = Vector{Vector{Vector{T} where T}}
const IndicesMap = Vector{Vector{T}} where {T<:Integer}
const FCfactorsT = Matrix{T} where {T<:AbstractFloat}
const FCproductT = Matrix{T} where {T<:AbstractFloat}

"""
    vibrationalIndices(aggCore)

Get all electronic indices of the [`Aggregate`](@ref).

# Arguments
* `aggCore`: Instance of [`Aggregate`](@ref).
"""
function electronicIndices(aggCore::AggregateCore)::ElIndices
    # TODO: add typeofs
    vibInds = Array{Array{Int64,1},1}(undef, 0)
    molCount = length(aggCore.molecules)
    currentInds = fill(1, (molCount))
    push!(vibInds, copy(currentInds))
    for indPos = 1:molCount
        currentInds[indPos] += 1
        push!(vibInds, copy(currentInds))
        currentInds[indPos] -= 1
    end
    return vibInds
end

"""
    vibrationalIndices(aggCore)

Get all vibrational indices of the [`Aggregate`](@ref).

"""
function vibrationalIndices(aggCore::AggregateCore)::VibIndices
    # TODO: add typeofs
    NvibMols = getNvib(aggCore)
    molLen = length(aggCore.molecules)
    NvibIndMols = Array{Array{Array{Int64,1},1},1}(undef, 0)
    for mol_i = 1:molLen
        push!(NvibIndMols, vibrationalIndices(NvibMols[mol_i]))
    end
    molNvib = map((vibInds) -> length(vibInds), NvibIndMols)
    molInds = vibrationalIndices(molNvib)
    aggInds = Array{Array{Array{Int64,1},1},1}(undef, 0)
    for molInd in molInds
        aggInd = Array{Array{Int64,1},1}(undef, 0)
        for mol_i = 1:molLen
            push!(aggInd, NvibIndMols[mol_i][molInd[mol_i]])
        end
        push!(aggInds, aggInd)
    end
    return aggInds
end

"""
    getIndices(aggCore)

Get all indices (electronic and vibrational) of the [`Aggregate`](@ref).

# Arguments
* `aggCore`: Instance of [`Aggregate`](@ref).
"""
function getIndices(aggCore::AggregateCore)::Indices
    vibInds = vibrationalIndices(aggCore)
    vibIndsLen = length(vibInds)
    elInds = electronicIndices(aggCore)
    elIndsLen = length(elInds)
    indices = Array{Array{Array{Int64,1} where Int64,1},1}(undef, 0)
    for el_i = 1:elIndsLen
        for vib_i = 1:vibIndsLen
            push!(indices, [elInds[el_i], vibInds[vib_i]])
        end
    end
    return indices
end

"""
getIndicesMap(aggCore)

Get pointers (integers) to the indices of the [`Aggregate`](@ref) separated by
electronic states (e.g. [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]]).

# Arguments
* `agg`: Instance of [`Aggregate`](@ref).
"""
function getIndicesMap(aggCore::AggregateCore, aggIndices::Indices)::IndicesMap
    aggIndLen = length(aggIndices)
    indicesMap = Array{Array{Int64,1},1}(undef, 0)
    elLen = length(aggCore.molecules) + 1
    for el_i = 1:elLen
        push!(indicesMap, Array{Int64,1}(undef, 0))
    end

    for I = 1:aggIndLen
        elind1, vibind1 = aggIndices[I]
        elOrder1 = OpenQuantumSystems.elIndOrder(elind1)
        push!(indicesMap[elOrder1], I)
    end
    return indicesMap
end

"""
    getFranckCondonFactors(aggCore, aggIndices)
    getFranckCondonFactors(aggCore)

Get Frack-Condon factors of the [`Aggregate`](@ref) in a form of matrix.

# Arguments
* `aggCore`: Instance of [`Aggregate`](@ref).
* `aggIndices`: Aggregate indices generated by [`getIndices`](@ref).
"""
function getFranckCondonFactors(aggCore::AggregateCore, aggIndices::Indices)::FCfactorsT
    aggIndLen = length(aggIndices)
    molLen = length(aggCore.molecules)
    FC = zeros(Float64, (aggIndLen, aggIndLen))
    for I = 1:aggIndLen
        elind1, vibind1 = aggIndices[I]
        for J = 1:aggIndLen
            elind2, vibind2 = aggIndices[J]
            if elind1 == elind2
                if vibind1 == vibind2
                    FC[I, J] = 1.0
                end
            else
                fc = 1.0::Float64
                for mi = 1:molLen
                    mol = aggCore.molecules[mi]
                    fc *=
                        getMolStateFC(mol, elind1[mi], vibind1[mi], elind2[mi], vibind2[mi])
                end
                FC[I, J] = fc
            end
        end
    end
    return FC
end

getFranckCondonFactors(aggCore::AggregateCore)::FCfactorsT =
    (aggIndices = getIndices(aggCore);
    getFranckCondonFactors(aggCore, aggIndices))

# TODO: resurrect sparse version
#=
"""
    getFranckCondonFactorsSparse(agg, aggIndices)
    getFranckCondonFactorsSparse(agg)

Sparse version of [`getFranckCondonFactors`](@ref).

"""
function getFranckCondonFactorsSparse(
    agg::C,
    aggIndices::Any
) where {C<:AbstractAggregateCore}
    if aggIndices === nothing
        aggIndices = getIndices(agg)
    end
    aggIndLen = length(aggIndices)
    molLen = length(agg.molecules)
    FC = spzeros(C2, aggIndLen, aggIndLen)
    for I = 1:aggIndLen
        elind1, vibind1 = aggIndices[I]
        for J = 1:aggIndLen
            elind2, vibind2 = aggIndices[J]
            if elind1 == elind2
                if vibind1 == vibind2
                    FC[I, J] = 1.0
                end
            else
                fc = 1.0::C2
                for mi = 1:molLen
                    mol = agg.molecules[mi]
                    fc *=
                        getMolStateFC(mol, elind1[mi], vibind1[mi], elind2[mi], vibind2[mi])
                end
                if fc != 0.0
                    FC[I, J] = fc
                end
            end
        end
    end
    return FC
end

getFranckCondonFactorsSparse(
    agg::C
) where {C<:AbstractAggregateCore} = getFranckCondonFactorsSparse(
    agg::C,
    nothing
)
=#

"""
    getFCProd(agg, FCFact, aggIndices, vibindices)

Get product of Franck-Condon factors. This way the trace over bath will be faster

"""
function getFCproduct(
    aggCore::AggregateCore,
    indices::Indices,
    indicesMap::IndicesMap,
    FCfactors::FCfactorsT,
)::FCproductT
    elLen = length(aggCore.molecules)
    aggIndLen = length(indices)
    vibLen = length(indicesMap[2])
    elLen += 1
    FCproduct = zeros(eltype(FCfactors), aggIndLen, aggIndLen)

    for I = 1:aggIndLen
        elind1, vibind1 = indices[I]
        elOrder1 = OpenQuantumSystems.elIndOrder(elind1)

        for J = 1:aggIndLen
            elind2, vibind2 = indices[J]
            elOrder2 = OpenQuantumSystems.elIndOrder(elind2)
            #=
            K1 = vibindices[elOrder1][1]
            K2 = vibindices[elOrder1][end]
            L1 = vibindices[elOrder2][1]
            L2 = vibindices[elOrder2][end]
            FCProd[I, J] = sum(FCFact[K1:K2, I] .* FCFact[J, L1:L2])
            =#

            for m = 1:vibLen
                # according to quantarhei, trace_over_vibrations()
                K = indicesMap[1][m]
                L = indicesMap[1][m]
                FCproduct[I, J] += FCfactors[K, I] * FCfactors[J, L]
            end
        end
    end
    return FCproduct
end

struct AggregateTools <: AbstractAggregateTools
    elIndices::ElIndices
    vibIndices::VibIndices
    indices::Indices
    indicesMap::IndicesMap
    FCfactors::FCfactorsT
    FCproduct::FCproductT
    basisSystem::GenericBasis
    basisBath::GenericBasis
    basis::GenericBasis
    bSystemSize::Int64
    bBathSize::Int64
    bSize::Int64
    function AggregateTools(
        elIndices::ElIndices,
        vibIndices::VibIndices,
        indices::Indices,
        indicesMap::IndicesMap,
        FCfactors::FCfactorsT,
        FCproduct::FCproductT,
        basisSystem::GenericBasis,
        basisBath::GenericBasis,
        basis::GenericBasis,
        bSystemSize::Int64,
        bBathSize::Int64,
        bSize::Int64,
    )::AggregateTools
        new(
            elIndices,
            vibIndices,
            indices,
            indicesMap,
            FCfactors,
            FCproduct,
            basisSystem,
            basisBath,
            basis,
            bSystemSize,
            bBathSize,
            bSize,
        )
    end
end

function AggregateTools(aggCore::AggregateCore)::AggregateTools
    elIndices = electronicIndices(aggCore)
    vibIndices = vibrationalIndices(aggCore)
    indices = getIndices(aggCore)
    indicesMap = getIndicesMap(aggCore, indices)
    FCfactors = getFranckCondonFactors(aggCore, indices)
    FCproduct = getFCproduct(aggCore, indices, indicesMap, FCfactors)

    bSystemSize = length(elIndices)
    bBathSize = length(vibIndices)
    bSize = length(indices)

    basisSystem = GenericBasis([bSystemSize])
    basisBath = GenericBasis([bBathSize])
    basis = GenericBasis([bSize])

    return AggregateTools(
        elIndices,
        vibIndices,
        indices,
        indicesMap,
        FCfactors,
        FCproduct,
        basisSystem,
        basisBath,
        basis,
        bSystemSize,
        bBathSize,
        bSize,
    )
end

Base.:(==)(x::AggregateTools, y::AggregateTools) =
    x.elIndices == y.elIndices &&
    x.vibIndices == y.vibIndices &&
    x.indices == y.indices &&
    x.FCfactors == y.FCfactors &&
    x.FCproduct == y.FCproduct &&
    x.bSystemSize == y.bSystemSize &&
    x.bBathSize == y.bBathSize &&
    x.bSize == y.bSize &&
    x.basisSystem == y.basisSystem &&
    x.basisBath == y.basisBath &&
    x.basis == y.basis

"""
    take_el_part(A, a, b, indicesMap)

Take electric part specified by electric indices `a` and `b` from the A (type of Array).

"""
function take_el_part(A::Array, a, b, indicesMap)
    a1 = indicesMap[a][1]
    a2 = indicesMap[a][end]
    b1 = indicesMap[b][1]
    b2 = indicesMap[b][end]

    return A[a1:a2, b1:b2]
end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
#Implement a data structure for graphs that allows modification (insertion, deletion).
#It should be possible to store values at edges and nodes.
#It might be easiest to use a dictionary of (node, edgelist) to do this.

using LightGraphs, Gadfly, GraphPlot

graph = SimpleGraph(5)
inp = "x"

function graphGraph(g)
    n=nv(g)
    gplothtml(graph, nodelabel=1:n)
end

function askInput(input="x")
    while input ∉ ["e","v","g"]
        println("Please type v for Verticie, e for Edge, g to Graph")
        input = readline()
	global inp = input
    end
    return  inp
end

function addEdgeNode(graph)
    if inp == "e"
        println("Please enter a edge connection between nodes")
        println("\n", "Edge 1:")
        ed1 = parse(Int64, readline())
        println("\n", "Edge 2:")
        ed2 = parse(Int64, readline())
        add_edge!(graph, ed1,ed2)
    else
       println("Please enter number of vertices to add")
       vert = readline()
       add_vertices!(graph, parse(Int64,vert))
    end
end

while true
    println("Graph Nodes:  ", 1:nv(graph))
    println("Graph Edges:  ", ne(graph))
    askInput()
	println(inp)
    if inp ∈ ["v","e"]
        addEdgeNode(graph)
    else
	graphGraph(graph)
    end
end

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
const adjectives = [
	"groundbreaking"
	"revolutionary"
	"important"
	"novel"
	"fun"
	"interesting"
	"fascinating"
	"exciting"
	"surprising"
	"remarkable"
	"wonderful"
	"stunning"
	"mini"
	"small"
	"tiny"
	"cute"
	"friendly"
	"wild"
]

const nouns = [
	"discovery"
	"experiment"
	"story"
	"journal"
	"notebook"
	"revelation"
	"computation"
	"creation"
	"analysis"
	"invention"
	"blueprint"
	"report"
	"science"
	"magic"
	"program"
	"notes"
	"lecture"
	"theory"
	"proof"
	"conjecture"
]

function cutename()
    titlecase(rand(adjectives)) * " " * rand(nouns)
end

function new_notebooks_directory()
    try
        path = joinpath(first(DEPOT_PATH), "pluto_notebooks")
        if !isdir(path)
            mkdir(path)
        end
        path
    catch
        homedir()
    end
end

function without_dotjl(path)
	extension = let
		parts = split(basename(path), '.')
		length(parts) == 1 ? "" : last(parts)
	end

	if startswith(extension, "jl")
		path[1:end-length(extension)-1]
	else
		path
	end
end

"""
Return `base` * `suffix` if the file does not exist yet.

If it does, return `base * sep * string(n) * suffix`, where `n` is the smallest natural number such that the file is new. (no 0 is not a natural number you snake)
"""
function numbered_until_new(base::AbstractString; sep=" ", suffix=".jl", create_file=true)
    chosen = base * suffix
    n = 1
    while isfile(chosen)
        n += 1
        chosen = base * sep * string(n) * suffix
	end
	if create_file
		touch(chosen)
	end
    chosen
end

"Like `cp` except we create the file manually (to fix permission issues). (It's not plagiarism if you use this function to copy homework.)"
function readwrite(from::AbstractString, to::AbstractString)
    write(to, read(from, String))
end

function tryexpanduser(path)
	try
		expanduser(path)
	catch ex
		path
	end
end

tamepath = abspath ∘ tryexpanduser


text:
using ClimateMachine
using ClimateMachine.ConfigTypes
using ClimateMachine.Mesh.Topologies:
    StackedCubedSphereTopology, equiangular_cubed_sphere_warp, grid1d
using ClimateMachine.Mesh.Grids:
    DiscontinuousSpectralElementGrid,
    VerticalDirection,
    HorizontalDirection,
    EveryDirection,
    min_node_distance
using ClimateMachine.Mesh.Filters
using ClimateMachine.DGMethods: DGModel, init_ode_state, remainder_DGModel
using ClimateMachine.DGMethods.NumericalFluxes:
    RusanovNumericalFlux,
    CentralNumericalFluxGradient,
    CentralNumericalFluxSecondOrder
using ClimateMachine.ODESolvers
using ClimateMachine.SystemSolvers
using ClimateMachine.VTK: writevtk, writepvtu
using ClimateMachine.GenericCallbacks:
    EveryXWallTimeSeconds, EveryXSimulationSteps
using ClimateMachine.Thermodynamics:
    air_density, soundspeed_air, internal_energy, PhaseDry_pT, PhasePartition
using ClimateMachine.TemperatureProfiles: IsothermalProfile
using ClimateMachine.Atmos:
    AtmosPhysics,
    AtmosModel,
    DryModel,
    NoPrecipitation,
    NoRadiation,
    NTracers,
    ConstantDynamicViscosity,
    vars_state,
    Gravity,
    HydrostaticState,
    AtmosAcousticGravityLinearModel,
    AtmosAcousticLinearModel,
    parameter_set
using ClimateMachine.Orientations:
    SphericalOrientation, gravitational_potential, altitude, latitude, longitude
using ClimateMachine.VariableTemplates: flattenednames

using CLIMAParameters
using CLIMAParameters.Planet: planet_radius
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using MPI, Logging, StaticArrays, LinearAlgebra, Printf, Dates, Test

const output_vtk = false

const ntracers = 1

function main()
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    polynomialorder = 5
    numelem_horz = 10
    numelem_vert = 5

    timeend = 60 * 60
    # timeend = 33 * 60 * 60 # Full simulation
    outputtime = 60 * 60

    expected_result = Dict()
    expected_result[Float64, true] = 9.5073337869322578e+13
    expected_result[Float64, false] = 9.5073455070673781e+13

    @testset "acoustic wave" begin
        for FT in (Float64,)# Float32)
            for explicit in (true, false)
                result = test_run(
                    mpicomm,
                    polynomialorder,
                    numelem_horz,
                    numelem_vert,
                    timeend,
                    outputtime,
                    ArrayType,
                    FT,
                    explicit,
                )
                @test result ≈ expected_result[FT, explicit]
            end
        end
    end
end

function test_run(
    mpicomm,
    polynomialorder,
    numelem_horz,
    numelem_vert,
    timeend,
    outputtime,
    ArrayType,
    FT,
    explicit_solve,
)

    setup = AcousticWaveSetup{FT}()

    _planet_radius::FT = planet_radius(param_set)
    vert_range = grid1d(
        _planet_radius,
        FT(_planet_radius + setup.domain_height),
        nelem = numelem_vert,
    )
    topology = StackedCubedSphereTopology(mpicomm, numelem_horz, vert_range)

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = polynomialorder,
        meshwarp = equiangular_cubed_sphere_warp,
    )
    hmnd = min_node_distance(grid, HorizontalDirection())
    vmnd = min_node_distance(grid, VerticalDirection())

    T_profile = IsothermalProfile(param_set, setup.T_ref)
    δ_χ = @SVector [FT(ii) for ii in 1:ntracers]

    fullphysics = AtmosPhysics{FT}(
        param_set;
        ref_state = HydrostaticState(T_profile),
        turbulence = ConstantDynamicViscosity(FT(0)),
        moisture = DryModel(),
        tracers = NTracers{length(δ_χ), FT}(δ_χ),
    )
    fullmodel = AtmosModel{FT}(
        AtmosLESConfigType,
        fullphysics;
        orientation = SphericalOrientation(),
        init_state_prognostic = setup,
        source = (Gravity(),),
    )
    dg = DGModel(
        fullmodel,
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )
    Q = init_ode_state(dg, FT(0))

    # The linear model which contains the fast modes
    # acousticmodel = AtmosAcousticLinearModel(fullmodel)
    acousticmodel = AtmosAcousticGravityLinearModel(fullmodel)

    vacoustic_dg = DGModel(
        acousticmodel,
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient();
        direction = VerticalDirection(),
        state_auxiliary = dg.state_auxiliary,
    )

    # Advection model is the difference between the fullmodel and acousticmodel.
    # This will be handled with explicit substepping (time step in between the
    # vertical and horizontally acoustic models)
    rem_dg = remainder_DGModel(dg, (vacoustic_dg,))

    # determine the time step for the model components
    acoustic_speed = soundspeed_air(param_set, FT(setup.T_ref))
    advection_speed = 1 # What's a reasonable number here?

    vacoustic_dt = vmnd / acoustic_speed
    remainder_dt = min(min(hmnd, vmnd) / advection_speed, hmnd / acoustic_speed)
    odesolver = if explicit_solve
        remainder_dt = 5vacoustic_dt

        nsteps_output = ceil(outputtime / remainder_dt)
        remainder_dt = outputtime / nsteps_output
        nsteps = ceil(Int, timeend / remainder_dt)
        @assert nsteps * remainder_dt ≈ timeend

        vacoustic_solver =
            LSRK54CarpenterKennedy(vacoustic_dg, Q; dt = vacoustic_dt)
        rem_solver = MRIGARKERK45aSandu(
            rem_dg,
            vacoustic_solver,
            Q;
            dt = remainder_dt,
        )

        rem_solver
    else
        vacoustic_dt = 200vacoustic_dt
        element_size = (setup.domain_height / numelem_vert)

        nsteps_output = ceil(outputtime / vacoustic_dt)
        vacoustic_dt = outputtime / nsteps_output
        nsteps = ceil(Int, timeend / vacoustic_dt)
        @assert nsteps * vacoustic_dt ≈ timeend

        rem_solver = LSRK54CarpenterKennedy(rem_dg, Q; dt = remainder_dt)
        vacoustic_solver = MRIGARKESDIRK24LSA(
            vacoustic_dg,
            LinearBackwardEulerSolver(ManyColumnLU(); isadjustable = false),
            rem_solver,
            Q;
            dt = vacoustic_dt,
        )

        vacoustic_solver
    end

    filterorder = 18
    filter = ExponentialFilter(grid, 0, filterorder)
    cbfilter = EveryXSimulationSteps(1) do
        Filters.apply!(Q, :, grid, filter, direction = VerticalDirection())
        nothing
    end

    eng0 = norm(Q)
    @info @sprintf(
        """Starting
           ArrayType       = %s
           FT              = %s
           polynomialorder = %d
           numelem_horz    = %d
           numelem_vert    = %d
           acoustic dt     = %.16e
           remainder_dt    = %.16e
           norm(Q₀)        = %.16e
           """,
        "$ArrayType",
        "$FT",
        polynomialorder,
        numelem_horz,
        numelem_vert,
        vacoustic_dt,
        remainder_dt,
        eng0
    )
    # Set up the information callback
    starttime = Ref(now())
    cbinfo = EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q)
            runtime = Dates.format(
                convert(DateTime, now() - starttime[]),
                dateformat"HH:MM:SS",
            )
            @info @sprintf """Update
                              simtime = %.16e
                              runtime = %s
                              norm(Q) = %.16e
                              """ gettime(odesolver) runtime energy
        end
    end
    callbacks = (cbinfo, cbfilter)

    if output_vtk
        # create vtk dir
        vtkdir =
            "vtk_acousticwave" *
            "_poly$(polynomialorder)_horz$(numelem_horz)_vert$(numelem_vert)" *
            "_$(ArrayType)_$(FT)" *
            "_$(explicit_solve ? "Explicit_Explicit" : "Implicit_Explicit")"
        mkpath(vtkdir)

        vtkstep = 0
        # output initial step
        do_output(mpicomm, vtkdir, vtkstep, dg, Q, fullmodel)

        # setup the output callback
        cbvtk = EveryXSimulationSteps(floor(outputtime / dt)) do
            vtkstep += 1
            Qe = init_ode_state(dg, gettime(odesolver))
            do_output(mpicomm, vtkdir, vtkstep, dg, Q, fullmodel)
        end
        callbacks = (callbacks..., cbvtk)
    end

    solve!(
        Q,
        odesolver;
        numberofsteps = nsteps,
        adjustfinalstep = false,
        callbacks = callbacks,
    )

    # final statistics
    engf = norm(Q)
    @info @sprintf """Finished
    norm(Q)                 = %.16e
    norm(Q) / norm(Q₀)      = %.16e
    norm(Q) - norm(Q₀)      = %.16e
    """ engf engf / eng0 engf - eng0
    engf
end

Base.@kwdef struct AcousticWaveSetup{FT}
    domain_height::FT = 10e3
    T_ref::FT = 300
    α::FT = 3
    γ::FT = 100
    nv::Int = 1
end

function (setup::AcousticWaveSetup)(problem, bl, state, aux, localgeo, t)
    # callable to set initial conditions
    FT = eltype(state)
    param_set = parameter_set(bl)

    λ = longitude(bl, aux)
    φ = latitude(bl, aux)
    z = altitude(bl, aux)

    β = min(FT(1), setup.α * acos(cos(φ) * cos(λ)))
    f = (1 + cos(FT(π) * β)) / 2
    g = sin(setup.nv * FT(π) * z / setup.domain_height)
    Δp = setup.γ * f * g
    p = aux.ref_state.p + Δp

    ts = PhaseDry_pT(param_set, p, setup.T_ref)
    q_pt = PhasePartition(ts)
    e_pot = gravitational_potential(bl.orientation, aux)
    e_int = internal_energy(ts)

    state.ρ = air_density(ts)
    state.ρu = SVector{3, FT}(0, 0, 0)
    state.energy.ρe = state.ρ * (e_int + e_pot)

    state.tracers.ρχ = @SVector [FT(ii) for ii in 1:ntracers]
    nothing
end

function do_output(
    mpicomm,
    vtkdir,
    vtkstep,
    dg,
    Q,
    model,
    testname = "acousticwave",
)
    ## name of the file that this MPI rank will write
    filename = @sprintf(
        "%s/%s_mpirank%04d_step%04d",
        vtkdir,
        testname,
        MPI.Comm_rank(mpicomm),
        vtkstep
    )

    statenames = flattenednames(vars_state(model, Prognostic(), eltype(Q)))
    auxnames = flattenednames(vars_state(model, Auxiliary(), eltype(Q)))
    writevtk(filename, Q, dg, statenames, dg.state_auxiliary, auxnames)

    ## Generate the pvtu file for these vtk files
    if MPI.Comm_rank(mpicomm) == 0
        ## name of the pvtu file
        pvtuprefix = @sprintf("%s/%s_step%04d", vtkdir, testname, vtkstep)

        ## name of each of the ranks vtk files
        prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
            @sprintf("%s_mpirank%04d_step%04d", testname, i - 1, vtkstep)
        end

        writepvtu(pvtuprefix, prefixes, (statenames..., auxnames...), eltype(Q))

        @info "Done writing VTK: $pvtuprefix"
    end
end

main()


text:
# Note: true_density(X::Array{T, 2}) must operate column-wise on X

"""
generate Particles(x, w) ~ true_density: x ~ approx_d, w = true_density(x) ./ approx_d(x)
* ll: log of true_density
"""
function sample_importance(approx_d::Distribution, ll::Function; n=1000, normalized=true, verbose=false)
    x = rand(approx_d, n)
    w_log = ll(x) - logpdf(approx_d, x)
    idx = isfinite(w_log)
    x = ndims(x)==1 ? x[idx] : x[:, idx]
    w_log = w_log[idx]
    if normalized
        w_log = w_log .- maximum(w_log) + 500.0 # prevent overflow
        w = exp(w_log)
        w = w / sum(w)
        idx = w .> 0. # discard particles with zero weight
    else
        w = exp(w_log)
        idx = isfinite(w) & (w .> 0.)
    end
    if verbose
        info(@sprintf("importance samples: %d out of %d valid (%.4f%%), max = %g, sum = %g",
        sum(idx), n, 100.*sum(idx)/n, maximum(w), 1.0))
    end
    w = w[idx]
    x = ndims(x)==1 ? x[idx] : x[:, idx]
    @assert !isempty(x)
    return Particles(x, w)
end

"goodness of fit KL(true||approx_d) as estimated by Particles, the higher the better"
function cross_entropy(approx_d::Distribution, ll::Function; n=1000)
    pt = sample_importance(approx_d, ll; n=n)
    return -dot(pt.w, logpdf(approx_d, pt.x))
end

"goodness of fit KL(true||(1-α) old_approx_d + α new_component)"
function cross_entropy{T<:Real}(old_approx_d::Distribution, new_component::Distribution,
                                α::Union{AbstractFloat, AbstractVector{T}}, ll::Function; n=1000)
    pt = sample_importance(old_approx_d, ll, n=n)
    return cross_entropy(old_approx_d, new_component, pt, α, ll)
end

"version with supplied Particles"
function cross_entropy(old_approx_d::Distribution, new_component::Distribution, pt::Particles,
                       α::AbstractFloat, ll::Function)
    0 <= α <= 1 || error("bad α = $α")
    z1 = pdf(old_approx_d, pt.x)
    z2 = pdf(new_component, pt.x)
    return -dot(pt.w, log((1-α) * z1 + α * z2))
end

"vector version"
function cross_entropy{T<:Real}(old_approx_d::Distribution, new_component::Distribution, pt::Particles,
                                α::AbstractVector{T}, ll::Function)
    return map(b -> cross_entropy(old_approx_d, new_component, pt, b, ll), α)
end

"best possible fitness_score if approx_d = true"
function self_entropy(true_d::Distribution; n=500000)
    x = rand(true_d, n)
    return -mean(logpdf(true_d, x))
end

"""
Compute helper integrals (up to a constant) with supplied samples x ~ p.

* return: ∫ p(x) [L / ((1-α)q + αh)] dx, ∫ p(x) [L(q-h) / ((1-α)q + αh)^2] dx as
          functions of α, both up to the *same* constant exp(-const_factor).
          This is to prevent numerical issue in Newton's method.
"""
function _helpers_dL_d2L_inclusive{T <: Real}(q::Distribution, h::Distribution,
                                             ll::Function, x::Array{T}; const_factor=0.)
    d = length(q)
    z = ll(x) .- const_factor
    # @show const_factor mean(z) maximum(z) minimum(z)
    qq = pdf(q, x)
    hh = pdf(h, x)
    # @show safemean(log(0.5 * qq + 0.5 * hh))
    return α -> safemean(exp(z - log((1-α) * qq + α * hh))), α -> safemean((qq - hh) ./ ((1-α) * qq + α * hh) .* exp(z - log((1-α) * qq + α * hh)))
end

"return approximate dL/dα and d2L/dα2 as a function of α, with n samples"
function dL_d2L_inclusive(q::Distribution, h::Distribution, ll::Function; n=10000)
    samples_from_q = rand(q, n)
    samples_from_h = rand(h, n)
    return dL_d2L_inclusive(q, h, ll, samples_from_q, samples_from_h)
end

"return approximate dL/dα and d2L/dα2 as a function of α, with supplied samples"
function dL_d2L_inclusive{T <: Real}(q::Distribution, h::Distribution, ll::Function,
                      samples_from_q::Array{T}, samples_from_h::Array{T})
    if length(q)==1
        @assert length(samples_from_q) == length(samples_from_h)
    else
        @assert size(samples_from_q) == size(samples_from_h)
    end
    # the same const factor must be used!
    if length(q)==1
        zz = x -> median( ll(x[1:10]) - log(0.5 * pdf(q, x[1:10]) + 0.5 * pdf(h, x[1:10])) )
    else
        zz = x -> median( ll(x[:, 1:10]) - log(0.5 * pdf(q, x[:, 1:10]) + 0.5 * pdf(h, x[:, 1:10])) )
    end
    const_factor = maximum(zz, (samples_from_q, samples_from_h)) # don't take min!

    dL_q, d2L_q = _helpers_dL_d2L_inclusive(q, h, ll, samples_from_q, const_factor=const_factor)
    dL_h, d2L_h = _helpers_dL_d2L_inclusive(q, h, ll, samples_from_h, const_factor=const_factor)
    return α -> dL_q(α) - dL_h(α), α -> d2L_q(α) - d2L_h(α)
end

"search α, with Stochatic Newton's method, to minimize inclusive KL(true||(1-α) old_approx_d + α new_component).
problem is convex."
function search_α(approx_d::Distribution, new_component::Distribution, ll::Function;
                  n=10000, tol=1e-5, max_iter=100, verbose=false, show_plot=false)
    # get stochastic approximates
    approx_dL, approx_d2L = dL_d2L_inclusive(approx_d, new_component, ll, n=n)
    if show_plot
        α_vec = linspace(0.01, 0.99, 20)
        print(lineplot(α_vec, map(approx_dL, α_vec), title="dL/dα"))
        print(lineplot(α_vec, map(approx_d2L, α_vec), title="d2L/dα2"))
    end
    α = 0.01
    iter = 0
    change = 1.0
    while true
        iter += 1
        if iter > max_iter
          warn(@sprintf("change=%.2f > %g, max_iter=%d reached", change, tol, max_iter))
          break
        end
        # Newton
        # @show approx_dL(α), approx_d2L(α), approx_dL(α)/approx_d2L(α)
        α_new = α - approx_dL(α) / approx_d2L(α)
        α_new = min(α_new, 1-1e-3)
        α_new = max(α_new, 1e-3)
        change = abs(α_new - α)
        α = α_new
        if change < tol
            break
        end
    end
    if verbose
        println(@sprintf "α = %5g (final change = %5g < tol = %5g) with %2d iterations." α change tol iter)
    end
    return α
end

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
function Base.summary(io::IO, A::AbstractDimArray{T,N}) where {T,N}
    if N == 0
        print(io, "0-dimensional ")
    elseif N == 1
        print(io, size(A, 1), "-element ")
    else
        print(io, join(size(A), "×"), " ")
    end
    printstyled(io, string(nameof(typeof(A)), "{$T,$N}"); color=:blue)
end

function Base.show(io::IO, mime::MIME"text/plain", A::AbstractDimArray{T,N}) where {T,N}
    lines = 0
    summary(io, A)
    _printname(io, name(A))
    lines += _printdims(io, mime, dims(A))
    !(isempty(dims(A)) || isempty(refdims(A))) && println(io)
    lines += _printrefdims(io, mime, refdims(A))
    println(io)
    ds = displaysize(io)
    ioctx = IOContext(io, :displaysize => (ds[1] - lines, ds[2]))
    _show_array(ioctx, mime, parent(A))
end

function Base.show(io::IO, mime::MIME"text/plain", dim::Dimension)
    get(io, :compact, false) && return _show_compact(io::IO, dim)

    printstyled(io, nameof(typeof(dim)); color=_dimcolor(io))
    if dim isa Dim
        printstyled(io, "{"; color=:red)
        printstyled(io, string(":", name(dim)); color=:yellow)
        printstyled(io, "}"; color=:red)
    else
        if nameof(typeof(dim)) != name(dim)
            _printname(io, name(dim))
        end
    end
    print(io, ": ")
    printstyled(io, "\n  val: ")
    _printdimindex(io, val(dim))

    if !(mode(dim) isa AutoMode)
        printstyled(io, "\n  mode: ")
        show(io, mime, mode(dim))
    end
    if !(metadata(dim) isa NoMetadata)
        printstyled(io, "\n  metadata: ")
        show(io, mime, metadata(dim))
    end
    println(io)
    show(io, mime, typeof(dim))
end

function Base.show(io::IO, mime::MIME"text/plain", metadata::Metadata{N}) where N
    print(io, "Metadata")
    if N !== Nothing
        print(io, "{")
        show(io, N)
        print(io, "}")
    end
    print(io, " of ")
    show(io, mime, val(metadata))
end
Base.show(io::IO, mime::MIME"text/plain", mode::IndexMode) = _printmode(io, mode)

function Base.show(io::IO, mime::MIME"text/plain", mode::AbstractSampled)
    _printmode(io, mode)
    print(io, " - ")
    _printorder(io, mode)
    print(io, " ")
    _printspan(io, mode)
    print(io, " ")
    _printsampling(io, mode)
end

function Base.show(io::IO, mime::MIME"text/plain", mode::AbstractCategorical)
    _printmode(io, mode)
    print(io, " - ")
    _printorder(io, mode)
end

# short printing version for dimensions
function _show_compact(io::IO, dim::Dimension)
    printstyled(io, nameof(typeof(dim)); color=_dimcolor(io))
    if name(dim) != nameof(typeof(dim))
        print(io, " (")
        print(io, name(dim))
        print(io, ")")
    end
    _printdimproperties(io, dim)
end
function _show_compact(io::IO, dim::Dim)
    color = _dimcolor(io)
    printstyled(io, "Dim{"; color=color)
    printstyled(io, string(":", name(dim)); color=:yellow)
    printstyled(io, "}"; color=color)
    _printdimproperties(io, dim)
end

_dimcolor(io) = get(io, :is_ref_dim, false) ? :magenta : :red

function _printname(io::IO, name)
    if !(name == Symbol("") || name isa NoName)
        print(io, " (")
        printstyled(io, string(name); color=:yellow)
        print(io, ")")
    end
end

# Note GeoData uses these
function _printdims(io::IO, mime, dims::Tuple)
    if isempty(dims)
        print(io, ": ")
        return 0
    end
    print(io, " with dimensions: ")
    return _layout_dims(io, mime, dims)
end

function _printrefdims(io::IO, mime, refdims::Tuple)
    if isempty(refdims)
        return 0
    end
    print(io, "and reference dimensions: ")
    ctx = IOContext(io, :is_ref_dim=>true, :show_dim_val=>true)
    lines = _layout_dims(ctx, mime, refdims)
    return lines
end

function _layout_dims(io, mime, dims::Tuple)
    ctx = IOContext(io, :compact=>true)
    if all(m -> m isa NoIndex, mode(dims))
        for d in dims[1:end-1]
            show(ctx, mime, d)
            print(io, ", ")
        end
        show(ctx, mime, dims[end])
        return 0
    else # Dims get a line each
        lines = 1
        for d in dims
            println(io)
            print(io, "  ")
            show(ctx, mime, d)
            lines += 2 # Often they wrap
        end
        return lines
    end
end

function _printdimproperties(io, dim::Dimension)
    mode(dim) isa NoIndex && return nothing
    print(io, ": ")
    _printdimindex(io, val(dim))
    if !(mode(dim) isa AutoMode)
        print(io, " (")
        show(io, MIME"text/plain"(), mode(dim))
        print(io, ")")
    end
    return nothing
end

_printdimindex(io, A) = printstyled(io, A; color=:cyan)
_printdimindex(io, A::AbstractRange) = printstyled(io, A; color=:cyan)
function _printdimindex(io, v::AbstractVector)
    s = string(eltype(v)) * "["
    svals =
    if length(v) > 2 && eltype(v) <: Dates.TimeType
        "$(v[1]), …, $(v[end])"
    elseif length(v) > 5
        "$(v[1]), $(v[2]), …, $(v[end-1]), $(v[end])"
    else
        join((string(va) for va in v), ", ")
    end
    printstyled(io, s * svals * "]"; color=:cyan)
end

_printmode(io, mode) = print(io, nameof(typeof(mode)))
_printorder(io, mode) = print(io, nameof(typeof(order(mode))))
_printspan(io, mode) = print(io, nameof(typeof(span(mode))))
_printsampling(io, mode) = print(io, nameof(typeof(sampling(mode))))

# Thanks to Michael Abbott for the following function
function _show_array(io::IO, mime, A::AbstractArray{T,0}) where T
    print(_ioctx(io, T), "\n", A[])
end
function _show_array(io::IO, mime, A::AbstractArray{T,1}) where T
    Base.print_matrix(_ioctx(io, T), A)
end
function _show_array(io::IO, mime, A::AbstractArray{T,2}) where T
    Base.print_matrix(_ioctx(io, T), A)
end
function _show_array(io::IO, mime, A::AbstractArray{T,N}) where {T,N}
    o = ones(Int, N-2)
    frame = A[:, :, o...]
    onestring = join(o, ", ")
    println(io, "[:, :, $(onestring)]")
    Base.print_matrix(_ioctx(io, T), frame)
    nremaining = prod(size(A,d) for d=3:N) - 1
    nremaining > 0 && print(io, "\n[and ", nremaining," more slices...]")
end

function _ioctx(io, T)
    IOContext(io, :compact=>true, :limit=>true, :typeinfo=>T)
end