text:
import LinearAlgebra.dot
function dot(u::AbstractLinearOperator, v::Union{AbstractVector, AbstractLinearOperator})
  if size(u, 1) == 1
    vu = Matrix(u')[:]
  elseif size(u, 2) == 1
    vu = Matrix(u)[:]
  else
    throw("Wrong dimension")
  end
  if size(v, 1) == 1
    vv = Matrix(v')[:]
  elseif size(v, 2) == 1
    vv = typeof(v) <: AbstractVector ? v : Matrix(v)[:]
  else
    throw("Wrong dimension")
  end
  if length(vu) != length(vv)
    throw("Wrong dimension")
  end
  return dot(vu, vv)
end

import Base.-
function -(u::AbstractVector, v::AbstractLinearOperator)
  vv = zeros(length(u))
  if size(v, 1) == 1
    vv = Matrix(v')[:]
  elseif size(v, 2) == 1
    vv = Matrix(v)[:]
  else
    throw("Wrong dimension")
  end
  if length(u) != length(vv)
    throw("Wrong dimension")
  end

  return u - vv
end


text:
module PkgTest

include("constrsparsereg_test.jl")
include("classopath_test.jl")
include("admm_test.jl")
include("genlasso_test.jl")

end

ext: jl


text:
# stochastic-optimization.jl
#
# Wrapper methods for stochastic optimization problems.

module StochOpt

using LinearAlgebra;
using Distributions;

include("linreg.jl");
include("poissonreg.jl");
include("phase-retrieval.jl");
include("utilities.jl");
include("logreg.jl");

# Update types are
#
#  :proximal  - Full proximal update
#  :truncated  - Truncated update
#  :sgm  -  Stochastic gradient method (linear approximation) update
#
# The problem types are
#
#  :linreg - Linear regression
#  :logreg - Logistic regression
#  :poisson - Poisson regression
#  :phase - Phase retrieval

# obj_gaps =
#   MultiStepsizeExperiment(problem_type, update_type; n_tests, n_sample,
#                           dim, epsilons, stepsizes, maxiter, condition,
#                           noise)
#
# Performs a number of experiments to determine how sensitive problems
# of type problem_type are to initial stepsizes for the given
# update_type. Returns a tensor of size
#
#  num_computed_objectives -by- length(stepsizes) -by- n_tests
#
# where the entry [i,j,k] in the tensor corresponds to the obective
# value gap at "iteration" i of the given method, during experiment k,
# using initial stepsize stepsizes[j].
#
# The input data is as follows:
#
#  problem_type - one of the problem types above
#  update_type - type of optimizer: :sgm, :proximal, :truncated.
#  n_tests - number of experiments to run
#  n_sample - sample size for individual experiments
#  dim - dimension of problem to solve
#        Note that for matrix completion, the dimension is taken to be
#        the *rank* of the desired solution.
#  stepsizes - Different stepsizes to use
#  maxiter - maximum number of iterations
#  noise - magnitude of noise to add
#  condition - condition number of data matrix
function MultiStepsizeExperiment(
  problem_type::Symbol = :linreg,
  update_type::Symbol = :proximal;
  n_tests::Int64 = 20,
  n_sample::Int64 = 400,
  dim::Int64 = 40,
  stepsizes = logspace(-1.5, log10(100), 13),
  step_pow::Float64 = 0.6,
  maxiter::Int64 = 20 * n_sample,
  compute_every::Int64 = floor(Int64, maxiter / 100),
  condition::Float64 = 1.0,
  noise::Float64 = 0.0)

  # First, get the matrices of different times to epsilon.
  num_stepsizes = length(stepsizes);
  objectives = zeros(floor(Int64, maxiter / compute_every) + 1,
                     num_stepsizes, n_tests);
  # Now iterate over number of tests.
  for test_ind = 1:n_tests
    println("*** Experiment ", test_ind, " of ", n_tests, " ***");

    opt_objval = 0.0;
    (A, b, x_opt) = GenerateData(problem_type, n_sample, dim, noise, condition);
    minimal_value = FindMinimalValue(problem_type, A, b, quiet=true);
    for step_ind = 1:length(stepsizes)
      init_step = stepsizes[step_ind];
      println("\tInit step = ", init_step);
      x_init = zeros(size(A, 2));
      if (problem_type == :phase)
        x_init = randn(size(A, 2)) / 2;
      end
      obj_vals = StochasticProxPoint(A, b;
                                     problem_type = problem_type,
                                     update_type = update_type,
                                     init_stepsize = init_step,
                                     step_pow = step_pow,
                                     maxiter = maxiter,
                                     compute_every = compute_every,
                                     x_init = x_init);
      # remove NaNs
      obj_vals[isnan.(obj_vals)] .= Inf;
      objectives[:, step_ind, test_ind] = obj_vals .- minimal_value;
    end
  end
  return objectives;
end

# obj_gaps =
#   MultiMinibatchExperiment(problem_type, update_type; n_tests, n_sample,
#                           dim, epsilons, stepsize, maxiter, condition,
#                           noise, minibatch_sizes)
#
# Performs a number of experiments to determine how sensitive problems
# of type problem_type are to initial minibatch size for the given
# update_type. Returns a tensor of size
#
#  num_computed_objectives -by- length(minibatch_sizes) -by- n_tests
#
# where the entry [i,j,k] in the tensor corresponds to the obective
# value gap at "iteration" i of the given method, during experiment k,
# using minibatch size minibatch_sizes[j].
#
# The input data is as follows:
#
#  problem_type - one of the problem types above
#  update_type - type of optimizer: :sgm, :proximal, :truncated.
#  n_tests - number of experiments to run
#  n_sample - sample size for individual experiments
#  dim - dimension of problem to solve
#        Note that for matrix completion, the dimension is taken to be
#        the *rank* of the desired solution.
#  stepsize - Initial stepsize to use
#  maxiter - maximum number of iterations
#  noise - magnitude of noise to add
#  condition - condition number of data matrix
function MultiMinibatchExperiment(
  problem_type::Symbol = :linreg,
  update_type::Symbol = :proximal;
  n_tests::Int64 = 20,
  n_sample::Int64 = 400,
  dim::Int64 = 40,
  stepsize::Float64 = 1.0,
  minibatch_sizes::Array{Int64} = [1, 2, 4, 8, 16, 32, 64],
  step_pow::Float64 = 0.6,
  maxiter::Int64 = 20 * n_sample,
  compute_every::Int64 = floor(Int64, maxiter / 100),
  condition::Float64 = 1.0,
  noise::Float64 = 0.0)

  # First, get the matrices of different times to epsilon.
  num_minibatches = length(minibatch_sizes);
  objectives = zeros(floor(Int64, maxiter / compute_every) + 1,
                     num_minibatches, n_tests);
  # Now iterate over number of tests.
  for test_ind = 1:n_tests
    println("*** Experiment ", test_ind, " of ", n_tests, " ***");

    opt_objval = 0.0;
    (A, b, x_opt) = GenerateData(problem_type, n_sample, dim, noise, condition);
    minimal_value = FindMinimalValue(problem_type, A, b, quiet=true);
    for minibatch_ind = 1:length(minibatch_sizes)
      x_init = zeros(size(A, 2));
      if (problem_type == :phase)
        x_init = randn(size(A, 2)) / 2;
      end
      minibatch_size = minibatch_sizes[minibatch_ind];
      println("Minibatch size = ", minibatch_size);
      init_stepsize = stepsize * sqrt(minibatch_size);
      obj_vals = StochasticProxPoint(A, b;
                                     problem_type = problem_type,
                                     update_type = update_type,
                                     init_stepsize = init_stepsize,
                                     minibatch_size = minibatch_size,
                                     step_pow = step_pow,
                                     maxiter = maxiter,
                                     compute_every = compute_every,
                                     x_init = x_init);
      # remove NaNs
      obj_vals[isnan.(obj_vals)] .= Inf;
      objectives[:, minibatch_ind, test_ind] = obj_vals .- minimal_value;
    end
  end
  return objectives;
end

# objectives = StochasticProxPoint(A, b;
#                                  problem_type, update_type,
#                                  init_stepsize, step_pow,
#                                  maxiter,
#                                  compute_every, x_init, minibatch_size)
#
# Performs a stochastic model-based minimization on data (A, b),
# beginning from initial point x_init.
#
# Returns a vector of (estimated) *gaps* to optimality in the
# objective.
function StochasticProxPoint(A::Matrix{Float64}, b::Vector;
                             problem_type::Symbol = :linreg,
                             update_type::Symbol = :proximal,
                             init_stepsize::Float64 = 1.0,
                             step_pow::Float64 = .5,
                             maxiter::Int64 = 20 * size(A, 1),
                             compute_every::Int64 = floor(Int64,
                                                          maxiter / 100),
                             x_init::Vector{Float64} = zeros(size(A, 2)),
                             minibatch_size::Int64 = 1)
  # First, get the optimizers and updaters
  (objective, updater) = GetUpdateAndObjective(problem_type, update_type);
  (mm, nn) = size(A);
  if (length(b) != mm)
    error("Mismatched data sizes");
  end
  # Now, run the stochastic model-based minimization method.
  objectives = zeros(floor(Int64, maxiter / compute_every) + 1);
  x = x_init;
  objectives[1] = objective(A, b, x);
  num_objs_computed = 1;
  ind = 1;
  for iter = 1:maxiter
    # Choose a random index
    stepsize = init_stepsize / (iter^step_pow);
    if (minibatch_size > 1)
      inds = ceil.(Int64, mm * rand(minibatch_size));
      x = updater(@view(A[inds, :]), @view(b[inds]), x, stepsize);
    else
      ind = ceil(Int64, mm * rand())
      x = updater(A[ind, :], b[ind], x, stepsize);
    end
    if (mod(iter, compute_every) == 0)
      num_objs_computed += 1;
      objectives[num_objs_computed] = objective(A, b, x);
    end
  end
  return objectives;
end

# (obj, updater) = GetUpdateAndObjective(problem_type::Symbol,
#                                        update_type::Symbol)
#
# Given the symbols for problem type and update type, returns an
# objective function and an update function. The semantics of each of
# these is as follows. See above for lists of problem types and update types.
#
# The obj function takes as arguments (A, b, x), where A is an m-by-n
# data matrix, b is an m-dimensional label (or measurement or target)
# vector, and x is an n-dimensional parameter being optimized over.
# For example, if the problem_type is :linreg, then we would have
#
#  obj(A, b, x) = norm(A * x - b)^2 / (2 * m)
#
# where (m, n) = size(A).
#
# The updater function takes four arguments and computes a model-based
# minimization update. In particular, if the stochastic optimization
# problem being solved is to minimize
#
#  f(x) = (1/m) sum_{i=1}^m F(x; A[i, :], b[i])
#
# the updater takes as an argument (a, b, x_0, alpha), where a is a
# single data vector (i.e. one row of the matrix A), b a single label,
# x the current point, and alpha the stepsize, and returns
#
#  x = updater(a, b, x_0, alpha)
#    = argmin_x { F_{x_0}(x; a, b) + norm(x - x_0)^2 / (2 * alpha) }
#
# This is the model based update.
function GetUpdateAndObjective(problem_type::Symbol,
                               update_type::Symbol)
  objective = LinearRegressionObj;
  updater = ZeroUpdate;
  if (problem_type == :linreg)
    objective = LinearRegressionObj;
    if (update_type == :proximal)
      updater = ProximalUpdateLinreg;
    elseif (update_type == :truncated)
      updater = TruncatedUpdateLinreg;
    elseif (update_type == :sgm)
      updater = SGUpdateLinreg;
    else
      error("Unknown update type ", update_type);
    end
  elseif (problem_type == :poisson)
    objective = PoissonObj;
    if (update_type == :proximal)
      updater = ProximalUpdatePoisson;
    elseif (update_type == :truncated)
      updater = TruncatedUpdatePoisson;
    elseif (update_type == :sgm)
      updater = SGUpdatePoisson;
    else
      error("Unknown update type ", update_type);
    end
  elseif (problem_type == :phase)
    objective = PhaseRetrievalObj;
    if (update_type == :proximal)
      updater = ProximalUpdatePhaseRetrieval;
    elseif (update_type == :truncated)
      updater = TruncatedUpdatePhaseRetrieval;
    elseif (update_type == :sgm)
      updater = SGUpdatePhaseRetrieval;
    else
      error("Unknown update type ", update_type);
    end
  elseif (problem_type == :logreg)
    objective = LogisticRegressionObj;
    if (update_type == :proximal)
      error("No proximal update for logistic regression");
    elseif (update_type == :truncated)
      updater = TruncatedUpdateLogistic;
    elseif (update_type == :sgm)
      updater = SGUpdateLogistic;
    else
      error("Unknown update type ", update_type);
    end
  else
    error("Unknown problem type ", problem_type);
  end
  return (objective, updater);
end

# (A, b, x) = GenerateData(problem_type::Symbol,
#                          mm::Int64, nn::Int64,
#                          noise::Float64, condition::Float64)
#
# Generates data matrix A, target vector b, and optimal point x for
# the given problem type. The data matrix A is of size m-by-n, b of
# size m, and x of size n.
function GenerateData(problem_type::Symbol,
                      mm::Int64, nn::Int64,
                      noise::Float64, condition::Float64)
  if (problem_type == :linreg)
    return GenerateLinearRegressionData(mm, nn,
                                        condition_number = condition,
                                        noise = noise);
  elseif (problem_type == :poisson)
    return GeneratePoissonData(mm, nn, condition_number = condition);
  elseif (problem_type == :phase)
    return GeneratePhaseRetrievalData(mm, nn, condition_number = condition);
  elseif (problem_type == :logreg)
    return GenerateLogisticData(mm, nn, condition_number = condition,
                                snr = min(1 / noise, 100.0));
  end
  error("Data generation not implemented for problem type ", problem_type);
  return 0;
end

# Default zero objective
function ZeroObjective(A::Matrix{Float64}, b::Real,
                       x::Vector{Float64})
  return 0;
end

# Default zero update just returns input x.
function ZeroUpdate(a::Vector{Float64}, b::Real,
                    x::Vector{Float64}, stepsize::Float64)
  return x;
end

# FindMinimalValue(problem_type::Symbol,
#                  A::Matrix{Float64}, b::Vector)
#
# Performs a gradient descent method with backtracking line search
# (unless there is a known closed-form solution) to give the minimal value
# of the specified problem.
function FindMinimalValue(problem_type::Symbol,
                          A::Matrix{Float64}, b::Vector;
                          quiet::Bool = true)
  (mm, nn) = size(A);
  if (!quiet)
    println("Finding minimal value for problem of type ", problem_type);
  end
  if (problem_type == :linreg)
    if (mm <= nn)
      return 0;
    end
    x_opt = A\b;
    return norm(A * x_opt - b)^2 / (2 * mm);
  elseif (problem_type == :poisson)
    # Do a gradient method with backtracking line search
    return GradientDescentToValue(A, b,
                                  x -> PoissonObj(A, b, x),
                                  x -> FullPoissonGradient(A, b, x));
  elseif (problem_type == :logreg)
    # Do a gradient method with backtracking line search
    return GradientDescentToValue(A, b,
                                  x -> LogisticRegressionObj(A, b, x),
                                  x -> LogisticRegressionGrad(A, b, x));
  elseif (problem_type == :phase || problem_type == :matrix)
    # The minimal value for these problems--as there is a true,
    # consistent solution--is 0.
    return 0.0;
  else
    error("Unknown problem type ", problem_type);
  end
end

# obj = GradientDescentToValue(A::Matrix{Float64}, b::Vector,
#                              F::Function, Grad::Function;
#                              x_init = zeros(size(A, 2)),
#                              tol::Float64 = 1e-5)
#
# Performs gradient descent with a backtracking (Armijo) line search
# on the given objective F with gradient Grad until the gradient norm
# is less than tol.  Returns the value of the function F at termination.
function GradientDescentToValue(A::Matrix{Float64}, b::Vector,
                                F::Function, Grad::Function;
                                x_init = zeros(size(A, 2)),
                                tol::Float64 = 1e-5,
                                quiet::Bool = true)
  x = x_init;
  curr_obj = F(x);
  alpha = .25;  # Multiplier for sufficient decrease
  beta = .5;  # Multiplier for stepsizes
  stepsize = 1.0;
  grad_norm = Inf;
  curr_iter = 0;
  while (grad_norm > tol)
    g = Grad(x);
    grad_norm = norm(g);
    x_next = x - stepsize * g;
    guess_next_obj = curr_obj + alpha * dot(g, x_next - x);
    next_obj = F(x_next);
    while (guess_next_obj < next_obj)
      stepsize = stepsize * beta;
      x_next = x - stepsize * g;
      guess_next_obj = curr_obj + alpha * dot(g, x_next - x);
      next_obj = F(x_next);
    end
    curr_obj = next_obj;
    x = x_next;
    stepsize = min(stepsize / beta, 1.0);
    if (mod(curr_iter, 100) == 0 && !quiet)
      println("\tGradient descent: iteration ", curr_iter,
              ". Gradient norm: ", grad_norm);
    end
    curr_iter += 1;
  end
  return curr_obj;
end

end  # module StochOpt

max_stars_repo_path: Machine Learning Summer School 2019 (London, UK)/tutorials/optimization/stochastic-optimization.jl


text:
number_training_data = 6

r = rand()
δ = 1e-2
number_holes = 5
lb = 0.1 # for a blue wavelength this corresponds to a minimum feature of ≈40 nm
ub = 0.95 - lb # this assumes a period of 0.95 λ

# create parameters
ps_array = create_parameters(number_holes, number_training_data, δ, lb, ub)

# reference simulation
refractive_indexes_ref = ones(3) * 1.45
ref_local_field = get_local_field(Float64[0.5], refractive_indexes=refractive_indexes_ref)[3]

# simulations
refractive_indexes_sim = Float64[1.0, 1.0, 1.45]

sim_local_fields = ComplexF32[get_local_field(ps_array[it_ps, :], refractive_indexes=refractive_indexes_sim)[3] for it_ps=1:size(ps_array)[1]]

sim_local_fields2 = zeros(ComplexF64, (number_holes+1)*number_training_data)
@inbounds Base.Threads.@threads for it_ps=1:size(ps_array)[1]
    sim_local_fields2[it_ps] = get_local_field(ps_array[it_ps, :], refractive_indexes=refractive_indexes_sim)[3]
end

@test all(sim_local_fields .≈ sim_local_fields2)


text:
# This script runs the tests for all models that are known to us. For
# registered models it will use the latest tagged version, for all other
# models it will use the latest `master` branch version.
#
# The script assumes that Mimi is deved into ~/.julia/dev/Mimi, and will
# run the tests with the version of Mimi that is in that location.
#
# To run the script, simply call it via
#
#   julia --color=yes test_all_models.jl
#

using Pkg
pkg_that_errored = []

# first set of packages to test
packages_to_test = [
    "MimiDICE2010" => ("https://github.com/anthofflab/MimiDICE2010.jl", "master"),
    "MimiDICE2013" => ("https://github.com/anthofflab/MimiDICE2013.jl", "master"),
    "MimiDICE2016" => ("https://github.com/AlexandrePavlov/MimiDICE2016.jl", "master"),
    "MimiDICE2016R2" => ("https://github.com/anthofflab/MimiDICE2016R2.jl", "master"),
    "MimiRICE2010" => ("https://github.com/anthofflab/MimiRICE2010.jl", "master"),
    "MimiFUND" => ("https://github.com/fund-model/MimiFUND.jl", "master"),
    "MimiPAGE2009" => ("https://github.com/anthofflab/MimiPAGE2009.jl", "master"),
    "MimiPAGE2020" => ("https://github.com/lrennels/MimiPAGE2020.jl", "master"),
    "MimiSNEASY" => ("https://github.com/anthofflab/MimiSNEASY.jl", "master"),
    "MimiFAIR" => ("https://github.com/anthofflab/MimiFAIR.jl", "master"),
    "MimiMAGICC" => ("https://github.com/anthofflab/MimiMAGICC.jl", "master"),
    "MimiHector" => ("https://github.com/anthofflab/MimiHector.jl", "master")
]

mktempdir() do folder_name

    Pkg.activate(folder_name)

    Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "..")))

    Pkg.add([i isa Pair ? PackageSpec(url=i[2][1], rev=i[2][2]) : PackageSpec(i) for i in packages_to_test])

    Pkg.resolve()

    for pkg_info in packages_to_test
        pkg = pkg_info isa Pair ? pkg_info[1] : pkg_info
        @info "Now testing $pkg."
        try
            Pkg.test(PackageSpec(pkg))
        catch err
            push!(pkg_that_errored, pkg)
        end
    end

    println()
    println()
    println()

    println("The following packages errored:")
    for p in pkg_that_errored
        println(p)
    end
end

# test separately because needs MimiFUND 3.8.6
packages_to_test = [
   "MimiIWG" => ("https://github.com/rffscghg/MimiIWG.jl", "master")
]

mktempdir() do folder_name

    Pkg.activate(folder_name)

    Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "..")))

    Pkg.add([i isa Pair ? PackageSpec(url=i[2][1], rev=i[2][2]) : PackageSpec(i) for i in packages_to_test])

    Pkg.resolve()

    for pkg_info in packages_to_test
        pkg = pkg_info isa Pair ? pkg_info[1] : pkg_info
        @info "Now testing $pkg."
        try
            Pkg.test(PackageSpec(pkg))
        catch err
            push!(pkg_that_errored, pkg)
        end
    end
end

println()
println()
println()

println("The following packages errored:")
for p in pkg_that_errored
    println(p)
end


text:
using AbstractOperators
using LinearAlgebra, FFTW, DSP, SparseArrays, RecursiveArrayTools
using Printf
using Random
using Test

include("utils.jl")
Random.seed!(0)

verb = true

@testset "AbstractOperators" begin

@testset "Linear operators" begin
  include("test_linear_operators.jl")
end

@testset "Non-Linear operators" begin
  include("test_nonlinear_operators.jl")
end

@testset "Linear Calculus rules" begin
  include("test_linear_operators_calculus.jl")
end

@testset "Nonlinear Calculus rules" begin
  include("test_nonlinear_operators_calculus.jl")
end

@testset "L-BFGS" begin
  include("test_lbfgs.jl")
end

@testset "Syntax shorthands" begin
  include("test_syntax.jl")
end

end


text:
module LBFGSB

using L_BFGS_B_jll

include("subroutine.jl")
export setulb, timer

include("wrapper.jl")
export L_BFGS_B,lbfgsb

end # module


text:
using FISH_MDS

mds_main()


text:
using HPAT
using MPI

@acc @inline function cndf2( In::Array{Float64,1} )
    Out = 0.5 .+ 0.5 .* erf(0.707106781 .* In);
    return Out;
end

@acc hpat function blackscholes(iterations)

    sptprice   = Float64[ 42.0 for i=1:iterations ]
    strike     = Float64[ 40.0 + (i / iterations) for i=1:iterations ]
    rate       = Float64[ 0.5 for i=1:iterations ]
    volatility = Float64[ 0.2 for i=1:iterations ]
    time       = Float64[ 0.5 for i=1:iterations ]

    logterm = log10(sptprice ./ strike);
    powterm = .5 .* volatility .* volatility;
    den = volatility .* sqrt(time);
    d1 = (((rate .+ powterm) .* time) .+ logterm) ./ den;
    d2 = d1 .- den;
    NofXd1 = cndf2(d1);
    NofXd2 = cndf2(d2);
    futureValue = strike .* exp(- rate .* time);
    c1 = futureValue .* NofXd2
    Call = sptprice .* NofXd1 .- c1;
    Put  = Call .- futureValue .+ sptprice;
    sum(Put)
end

n = 10^8

if (length(ARGS) > 0)
	n = parse(Int, ARGS[1])
end

rank = MPI.Comm_rank(MPI.COMM_WORLD)
pes = MPI.Comm_size(MPI.COMM_WORLD)

checksum = blackscholes(10000)

MPI.Barrier(MPI.COMM_WORLD)
t1 = time_ns()
checksum = blackscholes(n)
t2 = time_ns()

if rank==0
	println("nodes: ", pes, "\ntime: ", (t2-t1)/1.0e9, "\nchecksum: ", checksum)
end

MPI.Barrier(MPI.COMM_WORLD)


text:
using FastAI, FluxMPI

FluxMPI.Init()

data, blocks = loaddataset("imagenette2-160", (Image, Label))
data_ddp = DistributedDataContainer(data)

method = ImageClassificationSingle(blocks)
learner = methodlearner(method, data_ddp; callbacks=[ToGPU(), Metrics(accuracy)],
                        optimizer=DistributedOptimiser(ADAM()))

fitonecycle!(learner, 10)


text:
import Knet.Ops20: softmax, logsoftmax
using Knet.Ops20: _softmax, _logsoftmax
using Knet.KnetArrays: DevArray
using AutoGrad: AutoGrad, @primitive1, Value
using CUDA.CUDNN: CUDNN, unsafe_cudnnSoftmaxForward, unsafe_cudnnSoftmaxBackward #, handle
using CUDA.CUDNN: CUDNN_SOFTMAX_FAST, CUDNN_SOFTMAX_ACCURATE, CUDNN_SOFTMAX_LOG, CUDNN_SOFTMAX_MODE_INSTANCE
#include("cudnn_retry.jl") # @cudnn_retry

function logsoftmax(x::Union{DevArray,Value{<:DevArray}}; dims=:)
    if dims isa Colon || (ndims(x) === 1 && (1 in dims)) # TODO: compare speed to fallback when n>20000
        _x = reshape(x, (1,1,length(x),1))
        _y = _cudnnSoftmaxForward(_x, algo=CUDNN_SOFTMAX_LOG)
        y = reshape(_y, size(x))
    elseif unique(dims) == [1]
        _x = reshape(x, (1,1,size(x,1),:))
        _y = _cudnnSoftmaxForward(_x, algo=CUDNN_SOFTMAX_LOG)
        y = reshape(_y, size(x))
    elseif ndims(x) == 2 && unique(dims) == [2]
        y = permutedims(logsoftmax(permutedims(x); dims=1))
    else
        y = _logsoftmax(x; dims=dims)
    end
    return y
end

function softmax(x::Union{DevArray,Value{<:DevArray}}; dims=:)
    if dims isa Colon || (ndims(x) === 1 && (1 in dims)) # TODO: compare speed to fallback when n>20000
        _x = reshape(x, (1,1,length(x),1))
        _y = _cudnnSoftmaxForward(_x, algo=CUDNN_SOFTMAX_ACCURATE)
        y = reshape(_y, size(x))
    elseif unique(dims) == [1]
        _x = reshape(x, (1,1,size(x,1),:))
        _y = _cudnnSoftmaxForward(_x, algo=CUDNN_SOFTMAX_ACCURATE)
        y = reshape(_y, size(x))
    elseif ndims(x) == 2 && unique(dims) == [2]
        y = permutedims(softmax(permutedims(x); dims=1))
    else
        y = _softmax(x; dims=dims)
    end
    return y
end

function _cudnnSoftmaxForward(x::R; algo) where {T,R<:DevArray{T}}
    mode = CUDNN_SOFTMAX_MODE_INSTANCE
    y = similar(x)
    @cudnn_retry unsafe_cudnnSoftmaxForward(CUDNN.handle(), algo, mode, Ref(T(1)), TD4(x), x, Ref(T(0)), TD4(y), y)
    return y
end

function _cudnnSoftmaxBackward(y::R, dy::R; algo) where {T,R<:DevArray{T}}
    mode = CUDNN_SOFTMAX_MODE_INSTANCE
    dx = similar(y)
    @cudnn_retry unsafe_cudnnSoftmaxBackward(CUDNN.handle(), algo, mode, Ref(T(1)), TD4(y), y, TD4(dy), dy, Ref(T(0)), TD4(dx), dx)
    return dx
end

@primitive1 _cudnnSoftmaxForward(x;o...),dy,y _cudnnSoftmaxBackward(y,dy;o...)
@primitive1 _cudnnSoftmaxBackward(y,dy;o...),ddx,dx csb1(y,dy,dx,ddx;o...) csb2(y,dy,dx,ddx;o...)

function csb1(y,dy,dx,ddx;algo)
    cdim = ndims(y) - 1
    dims = ((1:cdim)...,) # (mode == 0 ? ((1:cdim)...,) : (cdim,))
    if algo==CUDNN_SOFTMAX_FAST || algo==CUDNN_SOFTMAX_ACCURATE
        ddx .* dy - dy .* sum(y .* ddx, dims=dims) - ddx .* sum(y .* dy, dims=dims)
    elseif algo==CUDNN_SOFTMAX_LOG
        -ddx .* exp.(y) .* sum(dy,dims=dims)
    else
        error("Unknown algo: $algo")
    end
end

function csb2(y,dy,dx,ddx;algo)
    cdim = ndims(y) - 1
    dims = ((1:cdim)...,) # (mode == 0 ? ((1:cdim)...,) : (cdim,))
    if algo==CUDNN_SOFTMAX_FAST || algo==CUDNN_SOFTMAX_ACCURATE
        y .* (ddx .- sum(y .* ddx, dims=dims))
    elseif algo==CUDNN_SOFTMAX_LOG
        ddx .- sum(ddx .* exp.(y), dims=dims)
    else
        error("Unknown algo: $algo")
    end
end

function TD4(x::DevArray)
    d = ndims(x)
    if d == 4 || d == 5
        TD(x)
    else
        n = size(x,d)
        m = div(length(x),n)
        TD(reshape(x,(1,1,m,n)))
    end
end

#=
mutable structdef enum
{
    CUDNN_SOFTMAX_FAST     = 0,         /* straightforward implementation */
    CUDNN_SOFTMAX_ACCURATE = 1,         /* subtract max from every point to avoid overflow */
    CUDNN_SOFTMAX_LOG      = 2          /* this was introduced at cudnnVersion 3000 */
} cudnnSoftmaxAlgorithm_t;

mutable structdef enum
{
    CUDNN_SOFTMAX_MODE_INSTANCE = 0,   /* compute the softmax over all C, H, W for each N */
    CUDNN_SOFTMAX_MODE_CHANNEL = 1     /* compute the softmax over all C for each H, W, N */
} cudnnSoftmaxMode_t;

=#


text:
import accrue.cryptoerase.runtime.Condition;

class H {
    final C c;
    public int{L /c.c H} y;
    public H(C c) {
	this.c = c;
	int{H} x = 2;
	if (x > 0) {
	    c.c.set();
	}
    }
}

class C {
    static final Condition c = new Condition();
    int{L /c H} x = 0;

    public static void main(String[] args) {
	C c = new C();
	H h = new H(c);
	c.x = 42;
	h.y = c.x;
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_dimension(obj::Euclidean3D)
    return jcall(obj, "getDimension", jint, ())
end

function get_instance(::Type{Euclidean3D})
    return jcall(Euclidean3D, "getInstance", Euclidean3D, ())
end

function get_sub_space(obj::Euclidean3D)
    return jcall(obj, "getSubSpace", Euclidean2D, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
module SossMLJ

import Distributions
import MLJBase
import MLJModelInterface
import MonteCarloMeasurements
import Soss
import Statistics
import Tables

export SossMLJModel
export predict_particles
export rms_distribution
export rms_expected
export rms_median

include("types.jl")

include("loss-functions.jl")
include("machine-operations.jl")
include("models.jl")
include("particles.jl")
include("predictors.jl")

end # end module SossMLJ


text:
 """
```
pseudo_measurement(m::Model1002{T},
    TTT::Matrix{T}, RRR::Matrix{T}, CCC::Vector{T}) where {T<:AbstractFloat}
```

Assign pseudo-measurement equation (a linear combination of states):

```
x_t = ZZ_pseudo*s_t + DD_pseudo
```
"""
function pseudo_measurement(m::Model1002{T},
                            TTT::Matrix{T},
                            RRR::Matrix{T},
                            CCC::Vector{T}) where {T<:AbstractFloat}

    endo      = m.endogenous_states
    endo_addl = m.endogenous_states_augmented
    pseudo    = m.pseudo_observables

    _n_states = n_states_augmented(m)
    _n_pseudo = n_pseudo_observables(m)

    # Compute TTT^10, used for Expected10YearRateGap, Expected10YearRate, and Expected10YearNaturalRate
    TTT10 = (1/40)*((UniformScaling(1.) - TTT)\(TTT - TTT^41))

    # Initialize pseudo ZZ and DD matrices
    ZZ_pseudo = zeros(_n_pseudo, _n_states)
    DD_pseudo = zeros(_n_pseudo)

    ##########################################################
    ## PSEUDO-OBSERVABLE EQUATIONS
    ##########################################################

    ## Output
    ZZ_pseudo[pseudo[:y_t],endo[:y_t]] = 1.

    ## Flexible Output
    ZZ_pseudo[pseudo[:y_f_t],endo[:y_f_t]] = 1.

    ## Natural Rate
    ZZ_pseudo[pseudo[:NaturalRate],endo[:r_f_t]] = 1.
    DD_pseudo[pseudo[:NaturalRate]]              = 100.0*(m[:rstar]-1.0)

    ## π_t
    ZZ_pseudo[pseudo[:π_t],endo[:π_t]] = 1.
    DD_pseudo[pseudo[:π_t]]            = 100*(m[:π_star]-1);

    ## Output Gap
    ZZ_pseudo[pseudo[:OutputGap],endo[:y_t]] = 1;
    ZZ_pseudo[pseudo[:OutputGap],endo[:y_f_t]] = -1;

    ## Ex Ante Real Rate
    ZZ_pseudo[pseudo[:ExAnteRealRate],endo[:R_t]]  = 1;
    ZZ_pseudo[pseudo[:ExAnteRealRate],endo[:Eπ_t]] = -1;
    DD_pseudo[pseudo[:ExAnteRealRate]]             = m[:Rstarn] - 100*(m[:π_star]-1);

    ## Long Run Inflation
    ZZ_pseudo[pseudo[:LongRunInflation],endo[:π_star_t]] = 1.
    DD_pseudo[pseudo[:LongRunInflation]]                 = 100. *(m[:π_star]-1.)

    ## Marginal Cost
    ZZ_pseudo[pseudo[:MarginalCost],endo[:mc_t]] = 1.

    ## Wages
    ZZ_pseudo[pseudo[:Wages],endo[:w_t]] = 1.

    ## Flexible Wages
    ZZ_pseudo[pseudo[:FlexibleWages],endo[:w_f_t]] = 1.

    ## Hours
    ZZ_pseudo[pseudo[:Hours],endo[:L_t]] = 1.

    ## Flexible Hours
    ZZ_pseudo[pseudo[:FlexibleHours],endo[:L_f_t]] = 1.

    ## z_t
    ZZ_pseudo[pseudo[:z_t], endo[:z_t]] = 1.

    ## Expected 10-Year Rate Gap
    ZZ_pseudo[pseudo[:Expected10YearRateGap], :] = TTT10[endo[:R_t], :] - TTT10[endo[:r_f_t], :] - TTT10[endo[:Eπ_t], :]

    ## Nominal FFR
    ZZ_pseudo[pseudo[:NominalFFR], endo[:R_t]] = 1.
    DD_pseudo[pseudo[:NominalFFR]] = m[:Rstarn]

    ## Expected 10-Year Interest Rate
    ZZ_pseudo[pseudo[:Expected10YearRate], :] = TTT10[endo[:R_t], :]
    DD_pseudo[pseudo[:Expected10YearRate]]    = m[:Rstarn]

    ## Expected 10-Year Natural Rate
    ZZ_pseudo[pseudo[:Expected10YearNaturalRate], :] = TTT10[endo[:r_f_t], :] + TTT10[endo[:Eπ_t], :]
    DD_pseudo[pseudo[:Expected10YearNaturalRate]]    = m[:Rstarn]

    ## Expected Nominal Natural Rate
    ZZ_pseudo[pseudo[:ExpectedNominalNaturalRate], endo[:r_f_t]] = 1.
    ZZ_pseudo[pseudo[:ExpectedNominalNaturalRate], endo[:Eπ_t]]  = 1.
    DD_pseudo[pseudo[:ExpectedNominalNaturalRate]]               = m[:Rstarn]

    ## Nominal Rate Gap
    ZZ_pseudo[pseudo[:NominalRateGap], endo[:R_t]]   = 1.
    ZZ_pseudo[pseudo[:NominalRateGap], endo[:r_f_t]] = -1.
    ZZ_pseudo[pseudo[:NominalRateGap], endo[:Eπ_t]]  = -1.

    ## Labor Productivity Growth
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo[:y_t]]           = 1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:y_t1]]     = -1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo[:z_t]]           = 1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:e_gdp_t]]  = 1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:e_gdp_t1]] = -m[:me_level]
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo[:L_t]]           = -1
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:L_t1]]     = 1.
    DD_pseudo[pseudo[:LaborProductivityGrowth]]                       = 100*(exp(m[:z_star]) - 1)

    return PseudoMeasurement(ZZ_pseudo, DD_pseudo)
end


text:
using CubicEoSDatabase

"""
    load(BrusilovskyEoSComponent; name::AbstractString[, component_dbs])

Return `BrusilovskyEoSComponent` named `name` with parameters from `component_dbs`.

# Arguments
- the first argument is `BrusilovskyEoSComponent`

# Keywords
- `name::AbstractString`: the component name
- `component_dbs`: an iterable of `ComponentDatabase` objects (by default,
(Data.martinez(), Data.brusilovsky_comp()))

# Returns
- `BrusilovskyEoSComponent`: the component generated by merging the EoS properties for the
given name from databases
"""
function load(
    ::Type{<:BrusilovskyEoSComponent};
    name::AbstractString,
    component_dbs = (Data.martinez(), Data.brusilovsky_comp())
)
    comp_properties = foldl(
        (dict, db) -> merge!(dict, getentry(db, name)),
        component_dbs;
        init = Dict{Symbol, Any}()
    )

    return __load_brusilovsky_comp__(; comp_properties...)
end

function __load_brusilovsky_comp__(
    ;
    name::AbstractString,
    molecular_mass,
    number_carbons,
    critical_temperature,
    critical_pressure,
    acentric_factor,
    eos_critical_compressibility,
    eos_critical_omega,
    eos_psi,
    extra_kw...  # unneeded keywords to construct component
)
    __load_brusilovsky_comp__(
        name,
        molecular_mass,
        number_carbons,
        critical_temperature,
        critical_pressure,
        acentric_factor,
        eos_critical_compressibility,
        eos_critical_omega,
        eos_psi
    )
end

function __load_brusilovsky_comp__(
    name::AbstractString,
    molecular_mass,
    number_carbons,
    critical_temperature,
    critical_pressure,
    acentric,
    critical_compressibility::AbstractString,
    critical_omega::AbstractString,
    psi::AbstractString
)
    Omegac  = 0.75001                       # Table (4.14)
    Zc      = 0.3357 - 0.0294 * acentric    # 4.203
    if acentric < 0.4489
        Psi = 1.050 + ( 0.105 + 0.482 * acentric ) * acentric   # 4.204
    else
        Psi = 0.429 + ( 1.004 + 1.561 * acentric ) * acentric   # 4.204
    end
    return BrusilovskyEoSComponent(
        name = name,
        critical_pressure = critical_pressure,
        critical_temperature = critical_temperature,
        acentric_factor = acentric,
        Omegac = Omegac,
        Zc = Zc,
        Psi = Psi,
        molar_mass = molecular_mass,
        carbon_number = number_carbons
    )
end

function __load_brusilovsky_comp__(
    name::AbstractString,
    molecular_mass,
    number_carbons,
    critical_temperature,
    critical_pressure,
    acentric,
    critical_compressibility,
    critical_omega,
    psi
)

    return BrusilovskyEoSComponent(
        ;
        name = name,
        critical_pressure = critical_pressure,
        critical_temperature = critical_temperature,
        acentric_factor = acentric,
        Omegac = critical_omega,
        Zc = critical_compressibility,
        Psi = psi,
        molar_mass = molecular_mass,
        carbon_number = number_carbons
    )
end

"""
    load(BrusilovskyEoSMixture; names[, component_dbs])

Return `BrusilovskyEoSMixture` by loading parameters for individual substances named `names`
from `component_dbs` and their interaction parameters from `mix_eos_db`.

# Arguments
- the first argument is `BrusilovskyEoSMixture`

# Keywords
- `names`: an iterable of `AbstractString`s containing the names of components to load
- `component_dbs`: an iterable of `ComponentDatabase` objects (by default,
`(Data.martinez(), Data.brusilovsky_comp())`)
- `mix_eos_db::MixtureDatabase`: the database with the cross-component interaction
parameters (by default, `Data.brusilovsky_mix()`)

# Returns
- `BrusilovskyEoSMixture`: the EoS for the mixture of the specified substances
"""
function load(
    ::Type{<:BrusilovskyEoSMixture};
    names,
    component_dbs = (Data.martinez(), Data.brusilovsky_comp()),
    mix_eos_db::MixtureDatabase = Data.brusilovsky_mix_adjusted(),
)
    components = [
        load(
            BrusilovskyEoSComponent;
            name = name,
            component_dbs = component_dbs
        ) for name in names
    ]
    corrections = getmatrix(mix_eos_db, names)
    return BrusilovskyEoSMixture(;
        components = components,
        corrections...
    )
end

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


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

import LinearAlgebra: checksquare

## sparse matrix multiplication

*(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB} =
    *(sppromote(A, B)...)
*(A::SparseMatrixCSC{TvA,TiA}, transB::Transpose{<:Any,<:SparseMatrixCSC{TvB,TiB}}) where {TvA,TiA,TvB,TiB} =
    (B = transB.parent; (pA, pB) = sppromote(A, B); *(pA, transpose(pB)))
*(A::SparseMatrixCSC{TvA,TiA}, adjB::Adjoint{<:Any,<:SparseMatrixCSC{TvB,TiB}}) where {TvA,TiA,TvB,TiB} =
    (B = adjB.parent; (pA, pB) = sppromote(A, B); *(pA, adjoint(pB)))
*(transA::Transpose{<:Any,<:SparseMatrixCSC{TvA,TiA}}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB} =
    (A = transA.parent; (pA, pB) = sppromote(A, B); *(transpose(pA), pB))
*(adjA::Adjoint{<:Any,<:SparseMatrixCSC{TvA,TiA}}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB} =
    (A = adjA.parent; (pA, pB) = sppromote(A, B); *(adjoint(pA), pB))
*(transA::Transpose{<:Any,<:SparseMatrixCSC{TvA,TiA}}, transB::Transpose{<:Any,<:SparseMatrixCSC{TvB,TiB}}) where {TvA,TiA,TvB,TiB} =
    (A = transA.parent; B = transB.parent; (pA, pB) = sppromote(A, B); *(transpose(pA), transpose(pB)))
*(adjA::Adjoint{<:Any,<:SparseMatrixCSC{TvA,TiA}}, adjB::Adjoint{<:Any,<:SparseMatrixCSC{TvB,TiB}}) where {TvA,TiA,TvB,TiB} =
    (A = adjA.parent; B = adjB.parent; (pA, pB) = sppromote(A, B); *(adjoint(pA), adjoint(pB)))

function sppromote(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB}
    Tv = promote_type(TvA, TvB)
    Ti = promote_type(TiA, TiB)
    A  = convert(SparseMatrixCSC{Tv,Ti}, A)
    B  = convert(SparseMatrixCSC{Tv,Ti}, B)
    A, B
end

# In matrix-vector multiplication, the correct orientation of the vector is assumed.

function mul!(C::StridedVecOrMat, A::SparseMatrixCSC, B::StridedVecOrMat, α::Number, β::Number)
    A.n == size(B, 1) || throw(DimensionMismatch())
    A.m == size(C, 1) || throw(DimensionMismatch())
    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    nzv = A.nzval
    rv = A.rowval
    if β != 1
        β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
    end
    for k = 1:size(C, 2)
        @inbounds for col = 1:A.n
            αxj = α*B[col,k]
            for j = A.colptr[col]:(A.colptr[col + 1] - 1)
                C[rv[j], k] += nzv[j]*αxj
            end
        end
    end
    C
end
*(A::SparseMatrixCSC{TA,S}, x::StridedVector{Tx}) where {TA,S,Tx} =
    (T = promote_type(TA, Tx); mul!(similar(x, T, A.m), A, x, one(T), zero(T)))
*(A::SparseMatrixCSC{TA,S}, B::StridedMatrix{Tx}) where {TA,S,Tx} =
    (T = promote_type(TA, Tx); mul!(similar(B, T, (A.m, size(B, 2))), A, B, one(T), zero(T)))

function mul!(C::StridedVecOrMat, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, B::StridedVecOrMat, α::Number, β::Number)
    A = adjA.parent
    A.n == size(C, 1) || throw(DimensionMismatch())
    A.m == size(B, 1) || throw(DimensionMismatch())
    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    nzv = A.nzval
    rv = A.rowval
    if β != 1
        β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
    end
    for k = 1:size(C, 2)
        @inbounds for col = 1:A.n
            tmp = zero(eltype(C))
            for j = A.colptr[col]:(A.colptr[col + 1] - 1)
                tmp += adjoint(nzv[j])*B[rv[j],k]
            end
            C[col,k] += α*tmp
        end
    end
    C
end
*(adjA::Adjoint{<:Any,<:SparseMatrixCSC{TA,S}}, x::StridedVector{Tx}) where {TA,S,Tx} =
    (A = adjA.parent; T = promote_type(TA, Tx); mul!(similar(x, T, A.n), adjoint(A), x, one(T), zero(T)))
*(adjA::Adjoint{<:Any,<:SparseMatrixCSC{TA,S}}, B::StridedMatrix{Tx}) where {TA,S,Tx} =
    (A = adjA.parent; T = promote_type(TA, Tx); mul!(similar(B, T, (A.n, size(B, 2))), adjoint(A), B, one(T), zero(T)))

function mul!(C::StridedVecOrMat, transA::Transpose{<:Any,<:SparseMatrixCSC}, B::StridedVecOrMat, α::Number, β::Number)
    A = transA.parent
    A.n == size(C, 1) || throw(DimensionMismatch())
    A.m == size(B, 1) || throw(DimensionMismatch())
    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    nzv = A.nzval
    rv = A.rowval
    if β != 1
        β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
    end
    for k = 1:size(C, 2)
        @inbounds for col = 1:A.n
            tmp = zero(eltype(C))
            for j = A.colptr[col]:(A.colptr[col + 1] - 1)
                tmp += transpose(nzv[j])*B[rv[j],k]
            end
            C[col,k] += α*tmp
        end
    end
    C
end
*(transA::Transpose{<:Any,<:SparseMatrixCSC{TA,S}}, x::StridedVector{Tx}) where {TA,S,Tx} =
    (A = transA.parent; T = promote_type(TA, Tx); mul!(similar(x, T, A.n), transpose(A), x, one(T), zero(T)))
*(transA::Transpose{<:Any,<:SparseMatrixCSC{TA,S}}, B::StridedMatrix{Tx}) where {TA,S,Tx} =
    (A = transA.parent; T = promote_type(TA, Tx); mul!(similar(B, T, (A.n, size(B, 2))), transpose(A), B, one(T), zero(T)))

# For compatibility with dense multiplication API. Should be deleted when dense multiplication
# API is updated to follow BLAS API.
mul!(C::StridedVecOrMat, A::SparseMatrixCSC, B::StridedVecOrMat) =
    mul!(C, A, B, one(eltype(B)), zero(eltype(C)))
mul!(C::StridedVecOrMat, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, B::StridedVecOrMat) =
    (A = adjA.parent; mul!(C, adjoint(A), B, one(eltype(B)), zero(eltype(C))))
mul!(C::StridedVecOrMat, transA::Transpose{<:Any,<:SparseMatrixCSC}, B::StridedVecOrMat) =
    (A = transA.parent; mul!(C, transpose(A), B, one(eltype(B)), zero(eltype(C))))

function (*)(X::StridedMatrix{TX}, A::SparseMatrixCSC{TvA,TiA}) where {TX,TvA,TiA}
    mX, nX = size(X)
    nX == A.m || throw(DimensionMismatch())
    Y = zeros(promote_type(TX,TvA), mX, A.n)
    rowval = A.rowval
    nzval = A.nzval
    @inbounds for multivec_row=1:mX, col = 1:A.n, k=A.colptr[col]:(A.colptr[col+1]-1)
        Y[multivec_row, col] += X[multivec_row, rowval[k]] * nzval[k]
    end
    Y
end

function (*)(D::Diagonal, A::SparseMatrixCSC)
    T = Base.promote_op(*, eltype(D), eltype(A))
    mul!(LinearAlgebra.copy_oftype(A, T), D, A)
end
function (*)(A::SparseMatrixCSC, D::Diagonal)
    T = Base.promote_op(*, eltype(D), eltype(A))
    mul!(LinearAlgebra.copy_oftype(A, T), A, D)
end

# Sparse matrix multiplication as described in [Gustavson, 1978]:
# http://dl.acm.org/citation.cfm?id=355796

*(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} = spmatmul(A,B)
*(A::SparseMatrixCSC{Tv,Ti}, B::Adjoint{<:Any,<:SparseMatrixCSC{Tv,Ti}}) where {Tv,Ti} = spmatmul(A, copy(B))
*(A::SparseMatrixCSC{Tv,Ti}, B::Transpose{<:Any,<:SparseMatrixCSC{Tv,Ti}}) where {Tv,Ti} = spmatmul(A, copy(B))
*(A::Transpose{<:Any,<:SparseMatrixCSC{Tv,Ti}}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} = spmatmul(copy(A), B)
*(A::Adjoint{<:Any,<:SparseMatrixCSC{Tv,Ti}}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} = spmatmul(copy(A), B)
*(A::Adjoint{<:Any,<:SparseMatrixCSC{Tv,Ti}}, B::Adjoint{<:Any,<:SparseMatrixCSC{Tv,Ti}}) where {Tv,Ti} = spmatmul(copy(A), copy(B))
*(A::Transpose{<:Any,<:SparseMatrixCSC{Tv,Ti}}, B::Transpose{<:Any,<:SparseMatrixCSC{Tv,Ti}}) where {Tv,Ti} = spmatmul(copy(A), copy(B))

function spmatmul(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti};
                  sortindices::Symbol = :sortcols) where {Tv,Ti}
    mA, nA = size(A)
    mB, nB = size(B)
    nA==mB || throw(DimensionMismatch())

    colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval
    colptrB = B.colptr; rowvalB = B.rowval; nzvalB = B.nzval
    # TODO: Need better estimation of result space
    nnzC = min(mA*nB, length(nzvalA) + length(nzvalB))
    colptrC = Vector{Ti}(undef, nB+1)
    rowvalC = Vector{Ti}(undef, nnzC)
    nzvalC = Vector{Tv}(undef, nnzC)

    @inbounds begin
        ip = 1
        xb = zeros(Ti, mA)
        x  = zeros(Tv, mA)
        for i in 1:nB
            if ip + mA - 1 > nnzC
                resize!(rowvalC, nnzC + max(nnzC,mA))
                resize!(nzvalC, nnzC + max(nnzC,mA))
                nnzC = length(nzvalC)
            end
            colptrC[i] = ip
            for jp in colptrB[i]:(colptrB[i+1] - 1)
                nzB = nzvalB[jp]
                j = rowvalB[jp]
                for kp in colptrA[j]:(colptrA[j+1] - 1)
                    nzC = nzvalA[kp] * nzB
                    k = rowvalA[kp]
                    if xb[k] != i
                        rowvalC[ip] = k
                        ip += 1
                        xb[k] = i
                        x[k] = nzC
                    else
                        x[k] += nzC
                    end
                end
            end
            for vp in colptrC[i]:(ip - 1)
                nzvalC[vp] = x[rowvalC[vp]]
            end
        end
        colptrC[nB+1] = ip
    end

    deleteat!(rowvalC, colptrC[end]:length(rowvalC))
    deleteat!(nzvalC, colptrC[end]:length(nzvalC))

    # The Gustavson algorithm does not guarantee the product to have sorted row indices.
    Cunsorted = SparseMatrixCSC(mA, nB, colptrC, rowvalC, nzvalC)
    C = SparseArrays.sortSparseMatrixCSC!(Cunsorted, sortindices=sortindices)
    return C
end

## solvers
function fwdTriSolve!(A::SparseMatrixCSCUnion, B::AbstractVecOrMat)
# forward substitution for CSC matrices
    nrowB, ncolB  = size(B, 1), size(B, 2)
    ncol = LinearAlgebra.checksquare(A)
    if nrowB != ncol
        throw(DimensionMismatch("A is $(ncol) columns and B has $(nrowB) rows"))
    end

    aa = getnzval(A)
    ja = getrowval(A)
    ia = getcolptr(A)

    joff = 0
    for k = 1:ncolB
        for j = 1:nrowB
            i1 = ia[j]
            i2 = ia[j + 1] - 1

            # loop through the structural zeros
            ii = i1
            jai = ja[ii]
            while ii <= i2 && jai < j
                ii += 1
                jai = ja[ii]
            end

            # check for zero pivot and divide with pivot
            if jai == j
                bj = B[joff + jai]/aa[ii]
                B[joff + jai] = bj
                ii += 1
            else
                throw(LinearAlgebra.SingularException(j))
            end

            # update remaining part
            for i = ii:i2
                B[joff + ja[i]] -= bj*aa[i]
            end
        end
        joff += nrowB
    end
    B
end

function bwdTriSolve!(A::SparseMatrixCSCUnion, B::AbstractVecOrMat)
# backward substitution for CSC matrices
    nrowB, ncolB = size(B, 1), size(B, 2)
    ncol = LinearAlgebra.checksquare(A)
    if nrowB != ncol
        throw(DimensionMismatch("A is $(ncol) columns and B has $(nrowB) rows"))
    end

    aa = getnzval(A)
    ja = getrowval(A)
    ia = getcolptr(A)

    joff = 0
    for k = 1:ncolB
        for j = nrowB:-1:1
            i1 = ia[j]
            i2 = ia[j + 1] - 1

            # loop through the structural zeros
            ii = i2
            jai = ja[ii]
            while ii >= i1 && jai > j
                ii -= 1
                jai = ja[ii]
            end

            # check for zero pivot and divide with pivot
            if jai == j
                bj = B[joff + jai]/aa[ii]
                B[joff + jai] = bj
                ii -= 1
            else
                throw(LinearAlgebra.SingularException(j))
            end

            # update remaining part
            for i = ii:-1:i1
                B[joff + ja[i]] -= bj*aa[i]
            end
        end
        joff += nrowB
    end
    B
end

ldiv!(L::LowerTriangular{T,<:SparseMatrixCSCUnion{T}}, B::StridedVecOrMat) where {T} = fwdTriSolve!(L.data, B)
ldiv!(U::UpperTriangular{T,<:SparseMatrixCSCUnion{T}}, B::StridedVecOrMat) where {T} = bwdTriSolve!(U.data, B)

(\)(L::LowerTriangular{T,<:SparseMatrixCSCUnion{T}}, B::SparseMatrixCSC) where {T} = ldiv!(L, Array(B))
(\)(U::UpperTriangular{T,<:SparseMatrixCSCUnion{T}}, B::SparseMatrixCSC) where {T} = ldiv!(U, Array(B))
\(A::Transpose{<:Real,<:Hermitian{<:Real,<:SparseMatrixCSC}}, B::Vector) = A.parent \ B
\(A::Transpose{<:Complex,<:Hermitian{<:Complex,<:SparseMatrixCSC}}, B::Vector) = copy(A) \ B
\(A::Transpose{<:Number,<:Symmetric{<:Number,<:SparseMatrixCSC}}, B::Vector) = A.parent \ B

function rdiv!(A::SparseMatrixCSC{T}, D::Diagonal{T}) where T
    dd = D.diag
    if (k = length(dd)) ≠ A.n
        throw(DimensionMismatch("size(A, 2)=$(A.n) should be size(D, 1)=$k"))
    end
    nonz = nonzeros(A)
    @inbounds for j in 1:k
        ddj = dd[j]
        if iszero(ddj)
            throw(LinearAlgebra.SingularException(j))
        end
        for i in nzrange(A, j)
            nonz[i] /= ddj
        end
    end
    A
end

rdiv!(A::SparseMatrixCSC{T}, adjD::Adjoint{<:Any,<:Diagonal{T}}) where {T} =
    (D = adjD.parent; rdiv!(A, conj(D)))
rdiv!(A::SparseMatrixCSC{T}, transD::Transpose{<:Any,<:Diagonal{T}}) where {T} =
    (D = transD.parent; rdiv!(A, D))

## triu, tril

function triu(S::SparseMatrixCSC{Tv,Ti}, k::Integer=0) where {Tv,Ti}
    m,n = size(S)
    if !(-m + 1 <= k <= n + 1)
        throw(ArgumentError(string("the requested diagonal, $k, must be at least ",
            "$(-m + 1) and at most $(n + 1) in an $m-by-$n matrix")))
    end
    colptr = Vector{Ti}(undef, n+1)
    nnz = 0
    for col = 1 : min(max(k+1,1), n+1)
        colptr[col] = 1
    end
    for col = max(k+1,1) : n
        for c1 = S.colptr[col] : S.colptr[col+1]-1
            S.rowval[c1] > col - k && break
            nnz += 1
        end
        colptr[col+1] = nnz+1
    end
    rowval = Vector{Ti}(undef, nnz)
    nzval = Vector{Tv}(undef, nnz)
    A = SparseMatrixCSC(m, n, colptr, rowval, nzval)
    for col = max(k+1,1) : n
        c1 = S.colptr[col]
        for c2 = A.colptr[col] : A.colptr[col+1]-1
            A.rowval[c2] = S.rowval[c1]
            A.nzval[c2] = S.nzval[c1]
            c1 += 1
        end
    end
    A
end

function tril(S::SparseMatrixCSC{Tv,Ti}, k::Integer=0) where {Tv,Ti}
    m,n = size(S)
    if !(-m - 1 <= k <= n - 1)
        throw(ArgumentError(string("the requested diagonal, $k, must be at least ",
            "$(-m - 1) and at most $(n - 1) in an $m-by-$n matrix")))
    end
    colptr = Vector{Ti}(undef, n+1)
    nnz = 0
    colptr[1] = 1
    for col = 1 : min(n, m+k)
        l1 = S.colptr[col+1]-1
        for c1 = 0 : (l1 - S.colptr[col])
            S.rowval[l1 - c1] < col - k && break
            nnz += 1
        end
        colptr[col+1] = nnz+1
    end
    for col = max(min(n, m+k)+2,1) : n+1
        colptr[col] = nnz+1
    end
    rowval = Vector{Ti}(undef, nnz)
    nzval = Vector{Tv}(undef, nnz)
    A = SparseMatrixCSC(m, n, colptr, rowval, nzval)
    for col = 1 : min(n, m+k)
        c1 = S.colptr[col+1]-1
        l2 = A.colptr[col+1]-1
        for c2 = 0 : l2 - A.colptr[col]
            A.rowval[l2 - c2] = S.rowval[c1]
            A.nzval[l2 - c2] = S.nzval[c1]
            c1 -= 1
        end
    end
    A
end

## diff

function sparse_diff1(S::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
    m,n = size(S)
    m > 1 || return SparseMatrixCSC(0, n, fill(one(Ti),n+1), Ti[], Tv[])
    colptr = Vector{Ti}(undef, n+1)
    numnz = 2 * nnz(S) # upper bound; will shrink later
    rowval = Vector{Ti}(undef, numnz)
    nzval = Vector{Tv}(undef, numnz)
    numnz = 0
    colptr[1] = 1
    for col = 1 : n
        last_row = 0
        last_val = 0
        for k = S.colptr[col] : S.colptr[col+1]-1
            row = S.rowval[k]
            val = S.nzval[k]
            if row > 1
                if row == last_row + 1
                    nzval[numnz] += val
                    nzval[numnz]==zero(Tv) && (numnz -= 1)
                else
                    numnz += 1
                    rowval[numnz] = row - 1
                    nzval[numnz] = val
                end
            end
            if row < m
                numnz += 1
                rowval[numnz] = row
                nzval[numnz] = -val
            end
            last_row = row
            last_val = val
        end
        colptr[col+1] = numnz+1
    end
    deleteat!(rowval, numnz+1:length(rowval))
    deleteat!(nzval, numnz+1:length(nzval))
    return SparseMatrixCSC(m-1, n, colptr, rowval, nzval)
end

function sparse_diff2(a::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
    m,n = size(a)
    colptr = Vector{Ti}(undef, max(n,1))
    numnz = 2 * nnz(a) # upper bound; will shrink later
    rowval = Vector{Ti}(undef, numnz)
    nzval = Vector{Tv}(undef, numnz)

    z = zero(Tv)

    colptr_a = a.colptr
    rowval_a = a.rowval
    nzval_a = a.nzval

    ptrS = 1
    colptr[1] = 1

    n == 0 && return SparseMatrixCSC(m, n, colptr, rowval, nzval)

    startA = colptr_a[1]
    stopA = colptr_a[2]

    rA = startA : stopA - 1
    rowvalA = rowval_a[rA]
    nzvalA = nzval_a[rA]
    lA = stopA - startA

    for col = 1:n-1
        startB, stopB = startA, stopA
        startA = colptr_a[col+1]
        stopA = colptr_a[col+2]

        rowvalB = rowvalA
        nzvalB = nzvalA
        lB = lA

        rA = startA : stopA - 1
        rowvalA = rowval_a[rA]
        nzvalA = nzval_a[rA]
        lA = stopA - startA

        ptrB = 1
        ptrA = 1

        while ptrA <= lA && ptrB <= lB
            rowA = rowvalA[ptrA]
            rowB = rowvalB[ptrB]
            if rowA < rowB
                rowval[ptrS] = rowA
                nzval[ptrS] = nzvalA[ptrA]
                ptrS += 1
                ptrA += 1
            elseif rowB < rowA
                rowval[ptrS] = rowB
                nzval[ptrS] = -nzvalB[ptrB]
                ptrS += 1
                ptrB += 1
            else
                res = nzvalA[ptrA] - nzvalB[ptrB]
                if res != z
                    rowval[ptrS] = rowA
                    nzval[ptrS] = res
                    ptrS += 1
                end
                ptrA += 1
                ptrB += 1
            end
        end

        while ptrA <= lA
            rowval[ptrS] = rowvalA[ptrA]
            nzval[ptrS] = nzvalA[ptrA]
            ptrS += 1
            ptrA += 1
        end

        while ptrB <= lB
            rowval[ptrS] = rowvalB[ptrB]
            nzval[ptrS] = -nzvalB[ptrB]
            ptrS += 1
            ptrB += 1
        end

        colptr[col+1] = ptrS
    end
    deleteat!(rowval, ptrS:length(rowval))
    deleteat!(nzval, ptrS:length(nzval))
    return SparseMatrixCSC(m, n-1, colptr, rowval, nzval)
end

diff(a::SparseMatrixCSC, dim::Integer)= dim==1 ? sparse_diff1(a) : sparse_diff2(a)

## norm and rank
vecnorm(A::SparseMatrixCSC, p::Real=2) = vecnorm(view(A.nzval, 1:nnz(A)), p)

function norm(A::SparseMatrixCSC,p::Real=2)
    m, n = size(A)
    if m == 0 || n == 0 || isempty(A)
        return float(real(zero(eltype(A))))
    elseif m == 1 || n == 1
        # TODO: compute more efficiently using A.nzval directly
        return norm(Array(A), p)
    else
        Tnorm = typeof(float(real(zero(eltype(A)))))
        Tsum = promote_type(Float64,Tnorm)
        if p==1
            nA::Tsum = 0
            for j=1:n
                colSum::Tsum = 0
                for i = A.colptr[j]:A.colptr[j+1]-1
                    colSum += abs(A.nzval[i])
                end
                nA = max(nA, colSum)
            end
            return convert(Tnorm, nA)
        elseif p==2
            throw(ArgumentError("2-norm not yet implemented for sparse matrices. Try norm(Array(A)) or norm(A, p) where p=1 or Inf."))
        elseif p==Inf
            rowSum = zeros(Tsum,m)
            for i=1:length(A.nzval)
                rowSum[A.rowval[i]] += abs(A.nzval[i])
            end
            return convert(Tnorm, maximum(rowSum))
        end
    end
    throw(ArgumentError("invalid p-norm p=$p. Valid: 1, Inf"))
end

# TODO rank

# cond
function cond(A::SparseMatrixCSC, p::Real=2)
    if p == 1
        normAinv = normestinv(A)
        normA = norm(A, 1)
        return normA * normAinv
    elseif p == Inf
        normAinv = normestinv(copy(A'))
        normA = norm(A, Inf)
        return normA * normAinv
    elseif p == 2
        throw(ArgumentError("2-norm condition number is not implemented for sparse matrices, try cond(Array(A), 2) instead"))
    else
        throw(ArgumentError("second argument must be either 1 or Inf, got $p"))
    end
end

function normestinv(A::SparseMatrixCSC{T}, t::Integer = min(2,maximum(size(A)))) where T
    maxiter = 5
    # Check the input
    n = checksquare(A)
    F = factorize(A)
    if t <= 0
        throw(ArgumentError("number of blocks must be a positive integer"))
    end
    if t > n
        throw(ArgumentError("number of blocks must not be greater than $n"))
    end
    ind = Vector{Int64}(undef, n)
    ind_hist = Vector{Int64}(undef, maxiter * t)

    Ti = typeof(float(zero(T)))

    S = zeros(T <: Real ? Int : Ti, n, t)

    function _rand_pm1!(v)
        for i in eachindex(v)
            v[i] = rand()<0.5 ? 1 : -1
        end
    end

    function _any_abs_eq(v,n::Int)
        for vv in v
            if abs(vv)==n
                return true
            end
        end
        return false
    end

    # Generate the block matrix
    X = Matrix{Ti}(undef, n, t)
    X[1:n,1] = 1
    for j = 2:t
        while true
            _rand_pm1!(view(X,1:n,j))
            yaux = X[1:n,j]' * X[1:n,1:j-1]
            if !_any_abs_eq(yaux,n)
                break
            end
        end
    end
    rmul!(X, inv(n))

    iter = 0
    local est
    local est_old
    est_ind = 0
    while iter < maxiter
        iter += 1
        Y = F \ X
        est = zero(real(eltype(Y)))
        est_ind = 0
        for i = 1:t
            y = norm(Y[1:n,i], 1)
            if y > est
                est = y
                est_ind = i
            end
        end
        if iter == 1
            est_old = est
        end
        if est > est_old || iter == 2
            ind_best = est_ind
        end
        if iter >= 2 && est <= est_old
            est = est_old
            break
        end
        est_old = est
        S_old = copy(S)
        for j = 1:t
            for i = 1:n
                S[i,j] = Y[i,j]==0 ? one(Y[i,j]) : sign(Y[i,j])
            end
        end

        if T <: Real
            # Check whether cols of S are parallel to cols of S or S_old
            for j = 1:t
                while true
                    repeated = false
                    if j > 1
                        saux = S[1:n,j]' * S[1:n,1:j-1]
                        if _any_abs_eq(saux,n)
                            repeated = true
                        end
                    end
                    if !repeated
                        saux2 = S[1:n,j]' * S_old[1:n,1:t]
                        if _any_abs_eq(saux2,n)
                            repeated = true
                        end
                    end
                    if repeated
                        _rand_pm1!(view(S,1:n,j))
                    else
                        break
                    end
                end
            end
        end

        # Use the conjugate transpose
        Z = F' \ S
        h_max = zero(real(eltype(Z)))
        h = zeros(real(eltype(Z)), n)
        h_ind = 0
        for i = 1:n
            h[i] = norm(Z[i,1:t], Inf)
            if h[i] > h_max
                h_max = h[i]
                h_ind = i
            end
            ind[i] = i
        end
        if iter >=2 && ind_best == h_ind
            break
        end
        p = sortperm(h, rev=true)
        h = h[p]
        permute!(ind, p)
        if t > 1
            addcounter = t
            elemcounter = 0
            while addcounter > 0 && elemcounter < n
                elemcounter = elemcounter + 1
                current_element = ind[elemcounter]
                found = false
                for i = 1:t * (iter - 1)
                    if current_element == ind_hist[i]
                        found = true
                        break
                    end
                end
                if !found
                    addcounter = addcounter - 1
                    for i = 1:current_element - 1
                        X[i,t-addcounter] = 0
                    end
                    X[current_element,t-addcounter] = 1
                    for i = current_element + 1:n
                        X[i,t-addcounter] = 0
                    end
                    ind_hist[iter * t - addcounter] = current_element
                else
                    if elemcounter == t && addcounter == t
                        break
                    end
                end
            end
        else
            ind_hist[1:t] = ind[1:t]
            for j = 1:t
                for i = 1:ind[j] - 1
                    X[i,j] = 0
                end
                X[ind[j],j] = 1
                for i = ind[j] + 1:n
                    X[i,j] = 0
                end
            end
        end
    end
    return est
end

# kron

function kron(a::SparseMatrixCSC{Tv,Ti}, b::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
    numnzA = nnz(a)
    numnzB = nnz(b)

    numnz = numnzA * numnzB

    mA,nA = size(a)
    mB,nB = size(b)

    m,n = mA*mB, nA*nB

    colptr = Vector{Ti}(undef, n+1)
    rowval = Vector{Ti}(undef, numnz)
    nzval = Vector{Tv}(undef, numnz)

    colptr[1] = 1

    colptrA = a.colptr
    colptrB = b.colptr
    rowvalA = a.rowval
    rowvalB = b.rowval
    nzvalA = a.nzval
    nzvalB = b.nzval

    col = 1

    @inbounds for j = 1:nA
        startA = colptrA[j]
        stopA = colptrA[j+1]-1
        lA = stopA - startA + 1

        for i = 1:nB
            startB = colptrB[i]
            stopB = colptrB[i+1]-1
            lB = stopB - startB + 1

            ptr_range = (1:lB) .+ (colptr[col]-1)

            colptr[col+1] = colptr[col] + lA * lB
            col += 1

            for ptrA = startA : stopA
                ptrB = startB
                for ptr = ptr_range
                    rowval[ptr] = (rowvalA[ptrA]-1)*mB + rowvalB[ptrB]
                    nzval[ptr] = nzvalA[ptrA] * nzvalB[ptrB]
                    ptrB += 1
                end
                ptr_range = ptr_range .+ lB
            end
        end
    end
    SparseMatrixCSC(m, n, colptr, rowval, nzval)
end

function kron(A::SparseMatrixCSC{Tv1,Ti1}, B::SparseMatrixCSC{Tv2,Ti2}) where {Tv1,Ti1,Tv2,Ti2}
    Tv_res = promote_type(Tv1, Tv2)
    Ti_res = promote_type(Ti1, Ti2)
    A = convert(SparseMatrixCSC{Tv_res,Ti_res}, A)
    B = convert(SparseMatrixCSC{Tv_res,Ti_res}, B)
    return kron(A,B)
end

kron(A::SparseMatrixCSC, B::VecOrMat) = kron(A, sparse(B))
kron(A::VecOrMat, B::SparseMatrixCSC) = kron(sparse(A), B)

function kron(x::SparseVector{Tv,Ti},y::SparseVector{Tv,Ti}) where {Tv,Ti}
    nnzx = nnz(x)
    nnzy = nnz(y)
    nnzz = nnzx*nnzy # number of nonzeros in new vector
    nzind = Vector{Ti}(undef, nnzz) # the indices of nonzeros
    nzval = Vector{Tv}(undef, nnzz) # the values of nonzeros
    @inbounds for i = 1:nnzx, j = 1:nnzy
        this_ind = (i-1)*nnzy+j
        nzind[this_ind] = (x.nzind[i]-1)*y.n + y.nzind[j]
        nzval[this_ind] = x.nzval[i] * y.nzval[j]
    end
    return SparseVector(x.n*y.n,nzind,nzval)
end

function kron(x::SparseVector{Tv1,Ti1}, y::SparseVector{Tv2,Ti2}) where {Tv1,Ti1,Tv2,Ti2}
    Tv_res = promote_type(Tv1, Tv2)
    Ti_res = promote_type(Ti1, Ti2)
    x2 = convert(SparseVector{Tv_res,Ti_res}, x)
    y2 = convert(SparseVector{Tv_res,Ti_res}, y)
    return kron(x2,y2)
end

kron(x::SparseVector{Tv,Ti}, y::AbstractVector) where {Tv,Ti} = kron(x, sparse(y))
kron(x::AbstractVector, y::SparseVector{Tv,Ti}) where {Tv,Ti} = kron(sparse(x), y)

## det, inv, cond

inv(A::SparseMatrixCSC) = error("The inverse of a sparse matrix can often be dense and can cause the computer to run out of memory. If you are sure you have enough memory, please convert your matrix to a dense matrix.")

# TODO

## scale methods

# Copy colptr and rowval from one sparse matrix to another
function copyinds!(C::SparseMatrixCSC, A::SparseMatrixCSC)
    if C.colptr !== A.colptr
        resize!(C.colptr, length(A.colptr))
        copyto!(C.colptr, A.colptr)
    end
    if C.rowval !== A.rowval
        resize!(C.rowval, length(A.rowval))
        copyto!(C.rowval, A.rowval)
    end
end

# multiply by diagonal matrix as vector
function mul!(C::SparseMatrixCSC, A::SparseMatrixCSC, D::Diagonal{<:Vector})
    m, n = size(A)
    b    = D.diag
    (n==length(b) && size(A)==size(C)) || throw(DimensionMismatch())
    copyinds!(C, A)
    Cnzval = C.nzval
    Anzval = A.nzval
    resize!(Cnzval, length(Anzval))
    for col = 1:n, p = A.colptr[col]:(A.colptr[col+1]-1)
        @inbounds Cnzval[p] = Anzval[p] * b[col]
    end
    C
end

function mul!(C::SparseMatrixCSC, D::Diagonal{<:Vector}, A::SparseMatrixCSC)
    m, n = size(A)
    b    = D.diag
    (m==length(b) && size(A)==size(C)) || throw(DimensionMismatch())
    copyinds!(C, A)
    Cnzval = C.nzval
    Anzval = A.nzval
    Arowval = A.rowval
    resize!(Cnzval, length(Anzval))
    for col = 1:n, p = A.colptr[col]:(A.colptr[col+1]-1)
        @inbounds Cnzval[p] = Anzval[p] * b[Arowval[p]]
    end
    C
end

function mul!(C::SparseMatrixCSC, A::SparseMatrixCSC, b::Number)
    size(A)==size(C) || throw(DimensionMismatch())
    copyinds!(C, A)
    resize!(C.nzval, length(A.nzval))
    mul!(C.nzval, A.nzval, b)
    C
end

function mul!(C::SparseMatrixCSC, b::Number, A::SparseMatrixCSC)
    size(A)==size(C) || throw(DimensionMismatch())
    copyinds!(C, A)
    resize!(C.nzval, length(A.nzval))
    mul!(C.nzval, b, A.nzval)
    C
end

function rmul!(A::SparseMatrixCSC, b::Number)
    rmul!(A.nzval, b)
    return A
end
function lmul!(b::Number, A::SparseMatrixCSC)
    lmul!(b, A.nzval)
    return A
end

function \(A::SparseMatrixCSC, B::AbstractVecOrMat)
    m, n = size(A)
    if m == n
        if istril(A)
            if istriu(A)
                return \(Diagonal(Vector(diag(A))), B)
            else
                return \(LowerTriangular(A), B)
            end
        elseif istriu(A)
            return \(UpperTriangular(A), B)
        end
        if ishermitian(A)
            return \(Hermitian(A), B)
        end
        return \(lufact(A), B)
    else
        return \(qrfact(A), B)
    end
end
for (xformtype, xformop) in ((:Adjoint, :adjoint), (:Transpose, :transpose))
    @eval begin
        function \(xformA::($xformtype){<:Any,<:SparseMatrixCSC}, B::AbstractVecOrMat)
            A = xformA.parent
            m, n = size(A)
            if m == n
                if istril(A)
                    if istriu(A)
                        return \($xformop(Diagonal(Vector(diag(A)))), B)
                    else
                        return \($xformop(LowerTriangular(A)), B)
                    end
                elseif istriu(A)
                    return \($xformop(UpperTriangular(A)), B)
                end
                if ishermitian(A)
                    return \($xformop(Hermitian(A)), B)
                end
                return \($xformop(lufact(A)), B)
            else
                return \($xformop(qrfact(A)), B)
            end
        end
    end
end

function factorize(A::SparseMatrixCSC)
    m, n = size(A)
    if m == n
        if istril(A)
            if istriu(A)
                return Diagonal(A)
            else
                return LowerTriangular(A)
            end
        elseif istriu(A)
            return UpperTriangular(A)
        end
        if ishermitian(A)
            return factorize(Hermitian(A))
        end
        return lufact(A)
    else
        return qrfact(A)
    end
end

# function factorize(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) where Ti
#     F = cholfact(A)
#     if LinearAlgebra.issuccess(F)
#         return F
#     else
#         ldltfact!(F, A)
#         return F
#     end
# end
function factorize(A::LinearAlgebra.RealHermSymComplexHerm{Float64,<:SparseMatrixCSC})
    F = cholfact(A)
    if LinearAlgebra.issuccess(F)
        return F
    else
        ldltfact!(F, A)
        return F
    end
end

chol(A::SparseMatrixCSC) = error("Use cholfact() instead of chol() for sparse matrices.")
lu(A::SparseMatrixCSC) = error("Use lufact() instead of lu() for sparse matrices.")
eig(A::SparseMatrixCSC) = error("Use IterativeEigensolvers.eigs() instead of eig() for sparse matrices.")

function Base.cov(X::SparseMatrixCSC; dims::Int=1, corrected::Bool=true)
    vardim = dims
    a, b = size(X)
    n, p = vardim == 1 ? (a, b) : (b, a)

    # The covariance can be decomposed into two terms
    # 1/(n - 1) ∑ (x_i - x̄)*(x_i - x̄)' = 1/(n - 1) (∑ x_i*x_i' - n*x̄*x̄')
    # which can be evaluated via a sparse matrix-matrix product

    # Compute ∑ x_i*x_i' = X'X using sparse matrix-matrix product
    out = Matrix(Base.unscaled_covzm(X, vardim))

    # Compute x̄
    x̄ᵀ = mean(X, dims=vardim)

    # Subtract n*x̄*x̄' from X'X
    @inbounds for j in 1:p, i in 1:p
        out[i,j] -= x̄ᵀ[i] * x̄ᵀ[j]' * n
    end

    # scale with the sample size n or the corrected sample size n - 1
    return rmul!(out, inv(n - corrected))
end

max_stars_repo_name: NymanLauri/julia

max_issues_repo_name: NymanLauri/julia

max_forks_repo_name: NymanLauri/julia


text:
doc"""
    DiscreteUniform(a,b)

A *Discrete uniform distribution* is a uniform distribution over a consecutive sequence of integers between `a` and `b`, inclusive.

$P(X = k) = 1 / (b - a + 1) \quad \text{for } k = a, a+1, \ldots, b.$

```julia
DiscreteUniform(a, b)   # a uniform distribution over {a, a+1, ..., b}

params(d)       # Get the parameters, i.e. (a, b)
span(d)         # Get the span of the support, i.e. (b - a + 1)
probval(d)      # Get the probability value, i.e. 1 / (b - a + 1)
minimum(d)      # Return a
maximum(d)      # Return b
```

External links

* [Discrete uniform distribution on Wikipedia](http://en.wikipedia.org/wiki/Uniform_distribution_(discrete))
"""
immutable DiscreteUniform <: DiscreteUnivariateDistribution
    a::Int
    b::Int
    pv::Float64

    function DiscreteUniform(a::Real, b::Real)
        @check_args(DiscreteUniform, a <= b)
        new(a, b, 1.0 / (b - a + 1))
    end
    DiscreteUniform(b::Real) = DiscreteUniform(0, b)
    DiscreteUniform() = new(0, 1, 0.5)
end

@distr_support DiscreteUniform d.a d.b

### Parameters

span(d::DiscreteUniform) = d.b - d.a + 1
probval(d::DiscreteUniform) = d.pv
params(d::DiscreteUniform) = (d.a, d.b)

### Show

show(io::IO, d::DiscreteUniform) = show(io, d, (:a, :b))

### Statistics

mean(d::DiscreteUniform) = middle(d.a, d.b)

median(d::DiscreteUniform) = middle(d.a, d.b)

var(d::DiscreteUniform) = (span(d)^2 - 1.0) / 12.0

skewness(d::DiscreteUniform) = 0.0

function kurtosis(d::DiscreteUniform)
    n2 = span(d)^2
    -1.2 * (n2 + 1.0) / (n2 - 1.0)
end

entropy(d::DiscreteUniform) = log(span(d))

mode(d::DiscreteUniform) = d.a
modes(d::DiscreteUniform) = [d.a:d.b]

### Evaluation

cdf(d::DiscreteUniform, x::Int) = (x < d.a ? 0.0 :
                                   x > d.b ? 1.0 :
                                   (floor(Int,x) - d.a + 1.0) * d.pv)

pdf(d::DiscreteUniform, x::Int) = insupport(d, x) ? d.pv : 0.0

logpdf(d::DiscreteUniform, x::Int) = insupport(d, x) ? log(d.pv) : -Inf

pdf(d::DiscreteUniform) = fill(probval(d), span(d))

function _pdf!(r::AbstractArray, d::DiscreteUniform, rgn::UnitRange)
    vfirst = round(Int, first(rgn))
    vlast = round(Int, last(rgn))
    vl = max(vfirst, d.a)
    vr = min(vlast, d.b)
    if vl > vfirst
        for i = 1:(vl - vfirst)
            r[i] = 0.0
        end
    end
    fm1 = vfirst - 1
    if vl <= vr
        pv = d.pv
        for v = vl:vr
            r[v - fm1] = pv
        end
    end
    if vr < vlast
        for i = (vr-vfirst+2):length(rgn)
            r[i] = 0.0
        end
    end
    return r
end

function _logpdf!(r::AbstractArray, d::DiscreteUniform, x::AbstractArray)
    lpv = log(probval(d))
    for i = 1:length(x)
        @inbounds r[i] = insupport(d, x[i]) ? lpv : -Inf
    end
    return r
end

quantile(d::DiscreteUniform, p::Float64) = d.a + floor(Int,p * span(d))

function mgf(d::DiscreteUniform, t::Real)
    a, b = d.a, d.b
    u = b - a + 1
    t == 0 ? 1.0 : (exp(t*a) * expm1(t*u)) / (u*expm1(t))
end

function cf(d::DiscreteUniform, t::Real)
    a, b = d.a, d.b
    u = b - a + 1
    t == 0 ? complex(1.0) : (im*cos(t*(a+b)/2) + sin(t*(a-b-1)/2)) / (u*sin(t/2))
end

### Sampling

rand(d::DiscreteUniform) = randi(d.a, d.b)

# Fit model

function fit_mle{T <: Real}(::Type{DiscreteUniform}, x::AbstractArray{T})
    if isempty(x)
        throw(ArgumentError("x cannot be empty."))
    end

    xmin = xmax = x[1]
    for i = 2:length(x)
        @inbounds xi = x[i]
        if xi < xmin
            xmin = xi
        elseif xi > xmax
            xmax = xi
        end
    end

    DiscreteUniform(xmin, xmax)
end


text:
rescaled(x::Array{Float32,1},i::Int) = Float32(i) .+ x.*(1.0f0/(2.0f0*maximum(abs.(x))))
rescaled(x::Array{Float64,1},i::Int) = Float64(i) .+ x.*(1.0/(2.0*maximum(abs.(x))))

# Formatted X labels
function xfmt(xmi::Int64, xma::Int64, fmt::String, N::Int64)
  dt = (xma-xmi)
  tzcorr = Libc.TmStruct(time())._11

  if fmt == "auto"
    # timespan < one minute
    x0 = div(xmi, 1000000)
    x1 = div(xma, 60000000)
    if x0 == x1
      fmt = "%S"
      xstr = string("Seconds from %s:00 [UTC]", Libc.strftime("%Y-%m-%d %H:%M", xmi*μs - tzcorr))
    else
      # timespan < one hour
      x0 = div(x0, 60)
      x1 = div(x1, 60)
      if x0 == x1
        fmt = "%M:%S"
        xstr = string("MM:SS from %s:00:00 [UTC]", Libc.strftime("%Y-%m-%d %H", xmi*μs - tzcorr))
      else
        # timespan < one day
        x0 = div(x0, 24)
        x1 = div(x1, 24)
        if x0 == x1
          fmt = "%H:%M:%S"
          xstr = string("Time [UTC], ", Libc.strftime("%Y-%m-%d", xmi*μs - tzcorr))
        else
          # timespan < one year
          if div(xmi, 31536000000000) == div(xma, 31536000000000)
            fmt = "%j, %T"
            xstr = string("JDN (", Libc.strftime("%Y", xmi*μs - tzcorr), "), Time [UTC]")
          else
            fmt ="%Y-%m-%d %H:%M:%S"
            xstr = string("DateTime [UTC]")
          end
        end
      end
    end
    # plot!(xlabel = xstr)
    PyPlot.xlabel(xstr,
    fontsize=12.0, color="black", fontweight="bold", family="serif")
  else
    # plot!(xlabel = string("Time [", replace(fmt, "%" => ""), "]"))
    PyPlot.xlabel(string("Time [", replace(fmt, "%" => ""), "]"),
    fontsize=12.0, color="black", fontweight="bold", family="serif")
  end

  dt = dt / (N-1)
  xt = Array{Float64,1}(undef, N)
  xl = Array{String,1}(undef, N)
  t = xmi
  for i = 1:N
    setindex!(xt, t, i)
    setindex!(xl, Libc.strftime(fmt, xt[i]*μs), i)
    t = t + dt
  end
  PyPlot.xlim(xmi, xma)
  PyPlot.xticks(xt, xl)
  return nothing
end


text:
FastRational{I1}(numden::Tuple{I2,I2}) where {I1<:Integer, I2<:Integer} = FastRational{I1}(numden[1]//numden[2])
Base.float(x::FastRational) = x.num / x.den

Base.Float64(x::FastRational) = float(x)
Base.Float32(x::FastRational) = Float32(float(x))
Base.Float16(x::FastRational) = Float16(float(x))
Base.BigFloat(x::FastRational) = BigFloat(x.num) / BigFloat(x.den)
Base.BigInt(x::FastRational) = isinteger(x) ? BigInt(x.num) ÷ BigInt(x.den) : throw(InexactError())

function FastRational(x::F; tol::Real=eps(x)) where F<:AbstractFloat
    !isfinite(x) && throw(DomainError("finite values only"))
    return FastRational(rationalize(x), tol=tol)
end
function FastRational{T}(x::F; tol::Real=eps(x)) where {T, F<:AbstractFloat}
    !isfinite(x) && throw(DomainError("finite values only"))
    return FastRational{T}(rationalize(T, x, tol=tol))
end

promote_rule(::Type{FastRational{S}}, ::Type{FastRational{T}}) where {T,S} = FastRational{promote_type(S,T)}
promote_rule(::Type{FastRational{S}}, ::Type{T}) where {T<:SUN,S} = FastRational{promote_type(S,T)}
promote_rule(::Type{FastRational{S}}, ::Type{T}) where {T,S} = promote_type(S,T)
promote_rule(::Type{FastRational{S}}, ::Type{Rational{T}}) where {T,S} = FastRational{promote_type(S,T)}


text:
using Pkg
Pkg.add("Suppressor")
Image = nothing
using Suppressor

function tupper_bitmap(k::BigInt)
    y_dr = divrem(fld(k, 17), 2)
    bin = []
    for i = 1:1802
        push!(bin, y_dr[2] == 1)
        y_dr = divrem(y_dr[1], 2)
    end
    return reshape(bin, (17, 106))
end

function visualize(bitmap)
    str = []
    for y = 1:17
        for x = 1:106
            push!(str, bitmap[y,107-x] ? 'X' : ' ')
        end
        push!(str, '\n')
    end
    return join(str)
end

function create_image(bitmap, col, type, file)
    img = Image.new("RGB", (17,106), "white")
    pixels = img.load()

    for y = 1:17
        for x = 1:106
            if bitmap[y,107-x]
                pixels[y,107-x] = col ? (Int64(trunc((y/17)*255)),Int64(trunc((x/106)*255)),100) : (0,0,0)
            end
        end
    end

    img.transpose(Image.ROTATE_90).save(file*"."*type)
end

function create_bin_file(bitmap, file)
end

params = Dict()
i = 1

while i < length(ARGS)
    global i
    if ARGS[i] == "--to-k"
        params["k"] = "NaN"
    else
        if ARGS[i] == "--from-k"
            params["k"] = ARGS[i+=1]
        elseif ARGS[i] == "--out"
            params["out"] = ARGS[i+=1]
        elseif ARGS[i] == "--in"
            params["in"] = ARGS[i+=1]
        elseif ARGS[i] == "--type"
            params["type"] = ARGS[i+=1]
        else
            exit(1)
        end
    end
    i += 1
end

if length(params) != 3
    exit(2)
end
if !haskey(params, "type")
    exit(3)
end
if !haskey(params, "k")
    exit(4)
end

if params["k"] == "NaN"
    if !haskey(params, "in")
        exit(5)
    end
else
    if !haskey(params, "out")
        exit(6)
    end

    try
        params["k"] = parse(BigInt, params["k"])
    catch cx
        exit(7)
    end
end

if params["type"] == "image" || params["type"] == "gay"
    @suppress begin
        global Image
        Pkg.add("PyCall")
        Pkg.add("Conda")
        using Conda
        Conda.add("Pillow")
        using PyCall
        Image = pyimport("PIL.Image")
    end
end

if typeof(params["k"]) == BigInt
    bitmap = tupper_bitmap(params["k"])
    if params["type"] == "text"
        print(visualize(bitmap))
    elseif params["type"] == "image"
        file_desc = split(params["out"], ".")
        @suppress create_image(bitmap, false, file_desc[2], file_desc[1])
    elseif params["type"] == "gay"
        file_desc = split(params["out"], ".")
        @suppress create_image(bitmap, true, file_desc[2], file_desc[1])
    end
end


text:
# ---
# title: 1288. Remove Covered Intervals
# id: problem1288
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Greedy, Sort, Line Sweep
# link: <https://leetcode.com/problems/remove-covered-intervals/description/>
# hidden: true
# ---
#
# Given a list of `intervals`, remove all intervals that are covered by another
# interval in the list.
#
# Interval `[a,b)` is covered by interval `[c,d)` if and only if `c <= a` and `b
# <= d`.
#
# After doing so, return _the number of remaining intervals_.
#
#
#
# **Example 1:**
#
#
#
#     Input: intervals = [[1,4],[3,6],[2,8]]
#     Output: 2
#     Explanation: Interval [3,6] is covered by [2,8], therefore it is removed.
#
#
# **Example 2:**
#
#
#
#     Input: intervals = [[1,4],[2,3]]
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: intervals = [[0,10],[5,12]]
#     Output: 2
#
#
# **Example 4:**
#
#
#
#     Input: intervals = [[3,10],[4,10],[5,11]]
#     Output: 2
#
#
# **Example 5:**
#
#
#
#     Input: intervals = [[1,2],[1,4],[3,4]]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= intervals.length <= 1000`
#   * `intervals[i].length == 2`
#   * `0 <= intervals[i][0] < intervals[i][1] <= 10^5`
#   * All the intervals are **unique**.
#
#
## @lc code=start
using LeetCode

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


text:
mutable struct VectorCache{T}
    ind_fill::Int64
    ind_max::Int64
    vec::Vector{T}
    function VectorCache{T}() where {T}
        isbitstype(T) || error("This cache requires an isbits type.")
        ind_max = 64
        new{T}(-9999, ind_max, Vector{T}(undef, ind_max))
    end
end

function expand!(vc::VectorCache{T}) where {T}  # TODO: make this function more elegant
    vc.ind_max += vc.ind_max
    resize!(vc.vec, vc.ind_max)
    return nothing
end

function returnNext(vc::VectorCache{T}) where {T}
    (vc.ind_fill == vc.ind_max) && expand!(vc)
    vc.ind_fill += 1
    return vc.vec[vc.ind_fill]
end

function addCacheItem!(vc::VectorCache{T}, item_T::T) where {T}
  vc.ind_fill += 1
  (vc.ind_max < vc.ind_fill) && expand!(vc)
  vc.vec[vc.ind_fill] = item_T
  return nothing
end

function Base.empty!(vc::VectorCache{T}) where {T}
  vc.ind_fill = 0
  return nothing
end

Base.isempty(vc::VectorCache{T}) where {T} = (vc.ind_fill == 0)
Base.length(vc::VectorCache{T}) where {T} = vc.ind_fill
Base.@propagate_inbounds Base.getindex(vc::VectorCache{T}, i::Int) where {T} = vc.vec[i]


text:
module JuFast5

using HDF5, Dates

export f5Reader, checkIntegrity, getAttr, getDataset,
       getASICId, getASICversion, getExpStartTime,
	   getRunId, getVersion, getSequencingKit,
	   getFlowcellType, isLocalBasecalled,
	   getChannelNumber, getSamplingRate, getDigitisation,
	   getFlowcellID, getDeviceID, getDeviceType,
	   getHeatsinkTemp, getSampleID, getProtocolRunID,
	   getRunID, getVersion, getProtocolsVersion,
	   getExperimentDurationSet, getUserFilenameInput,
	   getExperimentType, getFilename, getReadName,
	   getDuration, getMedianB4, getReadID, getReadNumber,
	   getStartMux, getStartTime

const requiredGroup = ["Raw", "UniqueGlobalKey", "PreviousReadInfo"]
const context_tags = "UniqueGlobalKey/context_tags"
const tracking_id = "UniqueGlobalKey/tracking_id"
const channel_id = "UniqueGlobalKey/channel_id"
const raw_reads = "Raw/Reads"

function checkIntegrity(hdf5Obj::HDF5File)
	for g in requiredGroup
		if !exists(hdf5Obj, g)
			error("Incomplete integrity! This Fast5 file lacks the group - ", g)
		end
	end
	nothing
end

function f5Reader(filename::AbstractString)
	f5 = h5open(filename, "r")
	checkIntegrity(f5)
	return f5
end

function getAttr(hdf5Obj::HDF5File, group::AbstractString, attr::AbstractString)
	if ! exists(hdf5Obj, group)
		error("Please check if given group is correct!")
	end

	if ! (attr in names(attrs(hdf5Obj[group])))
		error("Please check if given attribute name is correct!")
	end

	attrInfo = read(attrs(hdf5Obj[group]), attr)
	return attrInfo
end

# Get attributes from Raw/Reads/Read_n

function getReadName(hdf5Obj)
	return string(raw_reads, "/", iterate(keys(read(hdf5Obj[raw_reads])))[1])
end

function getDuration(hdf5Obj::HDF5File)
	attribute = "duration"
	path = getReadName(hdf5Obj)
	duration = getAttr(hdf5Obj, path, attribute)
	return duration
end

function getMedianB4(hdf5Obj::HDF5File)
	attribute = "median_before"
	path = getReadName(hdf5Obj)
	medianb4 = getAttr(jdf5Obj, path, attribute)
	return medianb4
end

function getReadID(hdf5Obj::HDF5File)
	attribute = "read_id"
	path = getReadName(hdf5Obj)
	readid = getAttr(hdf5Obj, path, attribute)
	return readid
end

function getReadNumber(hdf5Obj::HDF5File)
	attribute = "read_number"
	path = getReadName(hdf5Obj)
	readnumber = getAttr(hdf5Obj, path, attribute)
	return readnumber
end

function getStartMux(hdf5Obj::HDF5File)
	attribute = "start_mux"
	path = getReadName(hdf5Obj)
	startmux = getAttr(hdf5Obj, path, attribute)
	return startmux
end

function getStartTime(hdf5Obj::HDF5File)
	attribute = "start_time"
	path = getReadName(hdf5Obj)
	starttime = getAttr(hdf5Obj, path, attribute)
	return starttime
end

# Get attributes from UniqueGlobalKey/channel_id

function getChannelNumber(hdf5Obj::HDF5File)
	attribute = "channel_number"
	channelnumber = getAttr(hdf5Obj, channel_id, attribute)
	return channelnumber
end

function getSamplingRate(hdf5Obj::HDF5File)
	attribute = "sampling_rate"
	samplingrate = getAttr(hdf5Obj, channel_id, attribute)
	return samplingrate
end

function getDigitisation(hdf5Obj::HDF5File)
	attribute = "digitisation"
	digitisation = getAttr(hdf5Obj, channel_id, attribute)
	return digitisation
end

# Get attributes from UniqueGlobalKey/tracking_id

function getFlowcellID(hdf5Obj::HDF5File)
	attribute = "flow_cell_id"
	flowcellid = getAttr(hdf5Obj, tracking_id, attribute)
end

function getDeviceID(hdf5Obj::HDF5File)
	attribute = "device_id"
	deviceid = getAttr(hdf5Obj, tracking_id, attribute)
	return deviceid
end

function getDeviceType(hdf5Obj::HDF5File)
	attribute = "device_type"
	devicetype = getAttr(hdf5Obj, tracking_id, attribute)
	return devicetype
end

function getHeatsinkTemp(hdf5Obj::HDF5File)
	attribute = "heatsink_temp"
	headsinktemp = getAttr(hdf5Obj, tracking_id, attribute)
	return heatsinktemp
end

function getProtocolRunID(hdf5Obj::HDF5File)
	attribute = "protocol_run_id"
	protocolrunid = getAttr(hdf5Obj, tracking_id, attribute)
	return protocolrunid
end

function getProtocolsVersion(hdf5Obj::HDF5File)
	attribute = "protocols_version"
	protocolversion = getAttr(hdf5Obj, tracking_id, attribute)
	return protocolversion
end

function getSampleId(hdf5Obj::HDF5File)
	attribute = "sample_id"
	sampleid = getAttr(hdf5Obj, tracking_id, attribute)
	return sampleid
end

function getASICId(hdf5Obj::HDF5File)
	attribute = "asic_id"
	asicid = getAttr(hdf5Obj, tracking_id, attribute)
	return asicid
end

function getASICversion(hdf5Obj::HDF5File)
	attribute = "asic_version"
	asicversion = getAttr(hdf5Obj, tracking_id, attribute)
	return asicversion
end

function getExpStartTime(hdf5Obj::HDF5File)
	attribute = "exp_start_time"
	exp_start_time = getAttr(hdf5Obj, tracking_id, attribute)
	start_time = replace(exp_start_time, r"Z$" => s"")
	timestamp::DateTime = Dates.DateTime(start_time)
	return timestamp
end

function getRunId(hdf5Obj::HDF5File)
	attribute = "run_id"
	runid = getAttr(hdf5Obj, tracking_id, attribute)
	return runid
end

function getVersion(hdf5Obj::HDF5File)
	attribute = "version"
	version = getAttr(hdf5Obj, tracking_id, attribute)
	return version
end

# Get attributes from UniqueGlobalKey/context_tags

function getSequencingKit(hdf5Obj::HDF5File)
	attribute = "sequencing_kit"
	seqKit = getAttr(hdf5Obj, context_tags, attribute)
	return seqKit
end

function getFlowcellType(hdf5Obj::HDF5File)
	attribute = "flowcell_type"
	flowcell = getAttr(hdf5Obj, context_tags, attribute)
	return flowcell
end

function isLocalBasecalled(hdf5Obj::HDF5File)
	attribute = "local_basecalling"
	islocalbasecalled = getAttr(hdf5Obj, context_tags, attribute)
	if islocalbasecalled == "1"
		return true
	else
		return false
	end
end

function getExperimentDurationSet(hdf5Obj::HDF5File)
	attribute = "experiment_duration_set"
	exp_duration = getAttr(htf5Obj, context_tags, attribute)
	return parse(Int16, exp_duration)
end

function getUserFilenameInput(hdf5Obj::HDF5File)
	attribute = "user_filename_input"
	userinput = getAttr(hdf5Obj, context_tags, attribute)
	return userinput
end

function getExperimentType(hdf5Obj::HDF5File)
	attribute = "experiment_type"
	exptype = getAttr(hdf5Obj, context_tags, attribute)
	return exptype
end

function getFileName(hdf5Obj::HDF5File)
	attribute = "filename"
	filename = getAttr(hdf5Obj, context_tags, attribute)
	return filename
end

# Get datasets from Raw/Reads/Read_x

function getDataset(hdf5Obj::HDF5File)
	dataset = get_datasets(hdf5Obj)
	if length(dataset) > 1
		println("\033[1;33m[Warning] There is more than one read in this fast5 file.\033[0m")
	end

	if length(get_datasets(hdf5Obj)) < 1
		error("No data exists in this fast5 file.")
	end

	return read(dataset[1])
end

end


text:
using Pkg

Pkg.add(PackageSpec(url="https://github.com/JuliaLang/TOML.jl.git"))

using TOML

for k=keys(TOML.parsefile("./Project.toml")["deps"])
    Pkg.add(k)
end


text:
"""
`Linear`: Performs the operation w * x [.+ b]
# Constructors
    `Linear(output::Int, input::Int; kwargs)`

# Fields
    `w`: Weight
    `b`: Bias

# Keywords
      `bias=true` Whether or not to add bias
      `winit=xavier`  Weight initialization `(eltype, dims...)->array`
      `binit=zeros`   Bias initialization `(eltype, dims...)->array`
      `dtype=Float32` Element type of weight and bias

# Forward execution
    `forward(ctx, l::Linear, x)`
    `@mc l(x)`
    `@run l(x)

  If `ndims(x)` > 2, `mat(x)` is called to reshape `x`.
"""
type Linear <: KnetModule
    w::Param
    b::Union{Param, Void}
end

function Linear(output::Int, input::Int;
                bias=true,
                winit=xavier,
                binit=zeros,
                dtype=Float32)
    w = Param(winit(dtype, output, input))
    b = bias ? Param(binit(dtype, output, 1)) : nothing
    return Linear(w, b)
end

#function forward(ctx, l::Linear, x)
function (l::Linear)(ctx, x)
    if ndims(x) > 2;  x = mat(x); end
    o = val(ctx, l.w) * x
    if l.b !== nothing
        o = o .+ val(ctx, l.b)
    end
    return o
end


text:
"Set of small tests for all inference algorithms"
module TestNamespace

using Omega
using UnicodePlots
using Test
using InteractiveUtils: subtypes

function simple(ALG, op, v = 1.0)
  μ = normal(0.0, 1.0)
  x = normal(μ, 1.0)
  samples = rand(μ, op(x, v), 10; alg = ALG)
end

function testall()
  algs = [ALG() for ALG in subtypes(Omega.Inference.SamplingAlgorithm)]
  for ALG in filter(Omega.isapproximate, algs), op in [<ₛ, >ₛ, ==ₛ]
    println("Testing $ALG on $op")
    simple(ALG, op)
  end
end

testall()

end


text:
using Clapeyron, Test, Unitful

@testset "All tests" begin
    include("test_database.jl")
    include("test_solvers.jl")
    include("test_differentials.jl")
    include("test_models.jl")
    include("test_methods.jl")
end


text:
export RK4_naive
export step!

"""
    RK4_naive()

Runge-Kutta fourth order solver.

A naive version of `RK4`, without argument since no pre-allocation is performed.

"""
struct RK4_naive <: TimeSolver end

function step!(s  :: RK4_naive,
               m :: AbstractModel ,
               U  ,
               dt )

    U0 = copy(U)
    m.f!( U0 )
    U1 = copy(U0)

    U0 .= U .+ dt/2 .* U1
    m.f!( U0 )
    U2 = copy(U0)

    U0 .= U .+ dt/2 .* U2
    m.f!( U0 )
    U3 = copy(U0)

    U0 .= U .+ dt .* U3
    m.f!( U0 )
    U4 = copy(U0)

    U .+= dt/6 .* (U1 + 2*U2 + 2*U3 + U4 )

end