text:
import AbstractAlgebra.Generic: FreeModuleElem

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

# Private functions (will eventually change these functions into internal/private functions)

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

@doc Markdown.doc"""
    projective_quadratic_triple(Q::MatrixElem{T}, L::MatrixElem{T}, c::T, k) -> Tuple

For 1 <= k <= length(L), this function returns a tuple, where the first item in the tuple is
a kxk (positive definite) symmetric matrix, the next entry is a column vector of length k
and the last entry is a rational number.

THEORY:
R, (x1,x2,x3) = PolynomialRing(QQ, ["x1", "x2", "x3"])
Given a Polynomial ring of n variables we proceed as follows:
Generate a list of projective quadratic triples of n-1, ..., 1 variables,
given an n-variabled quadratic triple QT=[Q,L,c], where Q is a nxn symmetric
matrix, L is a column vector of length n and c is a rational number.
NOTE: if k is an int < n then it returns the k-variabled
projective quadratic triple.
"""
function projective_quadratic_triple(Q::MatrixElem{T}, L::MatrixElem{T}, c::T, k) where T <: RingElem
  n = nrows(Q)
  r1 = Vector{T}(undef,n-1); m1 = Vector{T}(undef,n-1);
  Q1 = Vector{typeof(Q)}(undef,n-1); p1 = Vector{typeof(L)}(undef,n-1);
  L1 = Vector{typeof(L)}(undef,n-1);
  for i in 1:n-1
    r1[i] = Q[n+1-i, n+1-i]
    m1[i] = L[n+1-i,1:1][1]
    Q1[i] = Q[1:n-i, 1:n-i] #n-i symmetric matrix extracted from the matrix Q
    p1[i] = Q[1:n-i, n+1-i:n+1-i] #Q[ri:rj, ci:cj] returns enrties from a matrix in the respective rows r_ and columns c_
    L1[i] = L[1:n-i,1:1] #column vector of length n-i
  end
  #------------------------------------------------
  c1 = c; pQT = Vector{typeof(Q)}(undef,n-1); pLL = Vector{typeof(L)}(undef,n-1); cc1 = Vector{T}();
  for j in 1:length(Q1)
    pt1 = transpose(p1[j])
    pQT[j] = Q1[j] - r1[j]^-1 * (p1[j] * pt1)
    pLL[j] = L1[j] - ((m1[j] * r1[j]^-1) * p1[j])
    c1 = c1 - (m1[j]^2 * r1[j]^-1) #rational number
    push!(cc1, c1)
  end
  #------------------------------------------------
  if k < n
    return pQT[n-k], pLL[n-k], cc1[n-k]
  else
    return Q,L,c
  end
end

@doc Markdown.doc"""
    range_ellipsoid_dim1(Q::MatrixElem{T}, L::MatrixElem{T}, c::T) ->  StepRange{fmpz, fmpz}

Returns a finite set of integers for which the inhomogeneous quadratic function,
given a one variabled quadratic triple, is less than or equal to zero.
"""
function range_ellipsoid_dim1(Q::MatrixElem{T}, L::MatrixElem{T}, c::T) where T <: RingElem
  @assert nrows(Q) == 1
  @assert L[1]^2 - Q[1]*c >= 0
  #x1, x2 are roots of q s.t x1 <= x2
  #x1 = (-L[1] - sqrt(L[1]^2 - Q[1]*c)) / Q[1]
  #x2 = (-L[1] + sqrt(L[1]^2 - Q[1]*c)) / Q[1]
  # return the integers in the interval [x1,x2]
  sqrt_floor(a::fmpz) = isqrt(a)
  a = -L[1] // Q[1]; b = (L[1]^2 - Q[1]*c) // Q[1]^2;
  cc = lcm(denominator(a), denominator(b))
  aa = FlintZZ(a * cc); bb = FlintZZ(b * cc^2);
  l = round(fmpz, (aa - sqrt_floor(bb))//cc, RoundUp)
  r = round(fmpz, (aa + sqrt_floor(bb))//cc, RoundDown)
  return l:r
end

@doc Markdown.doc"""
    positive_quadratic_triple(a::fmpz, Q::MatrixElem{T}, L::MatrixElem{T}, c::T) -> Tuple

This function computes a n-1 variabled quadraric triple i*(a, QT), given an n>1
variabled projective quadratic triple QT and an integer a.
Return value: an n-1 variabled positive quadratic triple i*(a, QT)
"""
function positive_quadratic_triple(a::fmpz, Q::MatrixElem{T}, L::MatrixElem{T}, c::T) where T <: RingElem
  n = nrows(Q)
  @assert n > 1
  #------------------------------------------------
  Q1 = Q[2:n, 2:n] #n-1 symmetric matrix extracted from the matrix Q
  L1 = a * Q[2:n, 1:1] + L[2:n,1:1] #Q[ri:rj, ci:cj] -> enrties from a matrix in the respective rows r_ and columns c_, L[i:n,1:1] column vector of length n-i
  c1 = a^2 * Q[1, 1] + 2 * a * L[1] + c #rational number
  #------------------------------------------------
  return Q1, L1, c1 #n-1 variabled quadratic triple
end

@doc Markdown.doc"""
    positive_quadratic_triple2(aa::Vector{fmpz}, Q::MatrixElem{T},
                               L::MatrixElem{T}, c::T) -> Tuple

NOTE: aa is an array of length k.

This function computes an n-k variabled quadratic triple i*(aa, QT), given an n
variabled quadratic triple QT and a tuple aa = {a1, ..., ak}, where k < n.
Return value: an n-k variabled positive quadratic triple i*(aa, QT).
"""
function positive_quadratic_triple2(aa::Vector{fmpz}, Q::MatrixElem{T}, L::MatrixElem{T}, c::T) where T <: RingElem
  QT = Q, L, c
  k = length(aa)
  @req length(L) >= 2 "The input quadratic triple should be at least 2 variabled."
  @req k < length(L) "The first input parameter should have length one less than the third input parameter"
  for i in 1:k
    R = positive_quadratic_triple(aa[i], QT[1], QT[2], QT[3])
    QT = R[1], R[2], R[3]
  end
  return QT
end

@doc Markdown.doc"""
    list_positive_quadratic_triple(aa::Vector{fmpz}, Q::MatrixElem{T},
                                   L::MatrixElem{T}, c::T) -> Array

Computes a list (for v = 1) QT_m^1 == i*(aa, QT_m^0), for aa in E(QT_1^0)
"""
function list_positive_quadratic_triple(aa::Vector{fmpz}, Q::MatrixElem{T}, L::MatrixElem{T}, c::T) where T <: RingElem
  n = nrows(Q)
  v = length(aa)+1
  t = n-length(aa)
  ListI = Vector(undef, t)
  for m in v:n
    P = projective_quadratic_triple(Q, L, c, m) #QT_m^0
    ListI[m-v+1] = positive_quadratic_triple2(aa, P[1], P[2], P[3])
  end
  return ListI
end

@doc Markdown.doc"""
    list_positive_quadratic_triple2(b::fmpz, ListIv::Vector{Any}) -> Array

For a fixed v, this function computes a list of QT_m^v := i*(b, QT_m^{v-1}),
where b is in E(QT_v^{v-1}).
"""
function list_positive_quadratic_triple2(b::fmpz, ListIv::Vector{Any}) #ListIv is the list of quadratic triples QT_m^{v-1}
  L = Vector(undef, length(ListIv)-1);
  for i in 1:length(ListIv)-1
    L[i] = positive_quadratic_triple(b, ListIv[i+1][1], ListIv[i+1][2], ListIv[i+1][3])
  end
  return L
end

@doc Markdown.doc"""
    posQuadTriple_intVectors(QT::Vector{Any}, E::Vector{Any}) -> Tuple

INPUT:
  A list of quadratic triples QT a list of tuples E of the form {a_1, ..., a_k}
OUTPUT:
  A new list of quadratic triples in the one dimension higher, and an updated
  list E of the form {a_1, ..., a_k, a_{k+1}}
"""
function posQuadTriple_intVectors(QT::Vector{Any}, E::Vector{Any})
  # this function could be rewritten for better performance at some point
  QTT1 = [];
  QTT = Array{Array}(undef, length(QT))
  bbb = [];
  for j in 1:length(QT) #for a fixed v
    b = range_ellipsoid_dim1(QT[j][1][1], QT[j][1][2], QT[j][1][3]) #E(QT_v^{v-1}) = {b_1, ..., b_N}
    b1 = collect(b)
    QT1 = Array{Array}(undef, length(b1))
    for k in 1:length(b1)
      QT1[k] = list_positive_quadratic_triple2(b1[k], QT[j]) #QT_m^v for m = v+1, ..., n
      for kk in E
        vb = vcat(kk, b1[k])
        if !(vb in bbb)
          push!(bbb,vb) #bb is now a list of tuples [x_1, ..., x_{v-1},b[k]] that satisfiy a inhomogeneous quad function for a v variabled quad triple
        end
      end
    end
    QTT[j] = QT1
  end
  for j in 1:length(QT)
    append!(QTT1, QTT[j])
  end
  return QTT1, bbb
end

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

# Public functions

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

@doc Markdown.doc"""
    convert_type(G::MatrixElem{T}, K::MatrixElem{T}, d::T) -> Tuple{ZLat, MatrixElem{T}, T}
Where T is a concrete type, e.g. fmpz, fmpq, etc.
Converts a quadratic triple QT = [Q, K, d] to the input values required for closest vector problem (CVP).
"""
function convert_type(G::MatrixElem{T}, K::MatrixElem{T}, d::T) where T <: RingElem
  if G[1,1] > 0
    Q = G
  else
    Q = -G
  end
  vector = -inv(Q) * K
  upperbound = (-transpose(K) * inv(Q) * -K)[1] - d
  Lattice = Zlattice(gram = Q)
  return Lattice, vector, upperbound
end

@doc Markdown.doc"""
    convert_type(L::ZLat, v::MatrixElem{T}, c::T) -> Tuple{fmpq_mat, fmpq_mat, fmpq}

Where T is a concrete type, e.g. fmpz, fmpq, etc.
Converts the input values from closest vector enumeration (CVE) to the corresponding quadratic triple QT = [Q, K, d].
"""
function convert_type(L::ZLat, v::MatrixElem{T}, c::T) where T <: RingElem
  V = ambient_space(L)
  Q = gram_matrix(V)
  K = -Q*v
  v = vec([i for i in v])
  d = inner_product(V,v,v)-c
  return Q, K, d
end

@doc Markdown.doc"""
    closest_vectors(Q::MatrixElem{T}, L::MatrixElem{T}, c::T; equal::Bool=false, sorting::Bool=false)
                                                    -> Vector{Vector{fmpz}}

Return all the integer vectors `x` of length n such that the inhomogeneous
quadratic function `q_{QT}(x) := xQx + 2xL + c <= 0` corresponding to an n variabled
quadratic triple. If the optional argument ``equal = true``, it return
all vectors `x` such that `q_{QT}(x) = 0`. By default ``equal = false``.
If the argument ``sorting = true``, then we get a a list of sorted vectors.
The Default value for ``sorting`` is set to ``false``.
"""
function closest_vectors(G::MatrixElem{T}, L::MatrixElem{T}, c::T; equal::Bool=false, sorting::Bool=false) where T <: RingElem
  #1 < v <= n+1, a = [a_1, ..., a_{v-1}] int tuple &
  @req det(G) != 0 "The symmetric matrix is not definite."
  if G[1,1] > 0
    Q = G
  else
    Q = -G
  end
  n = nrows(Q)
  P = projective_quadratic_triple(Q, L, c, 1)
  a = range_ellipsoid_dim1(P[1], P[2], P[3])#E(QT_1^0)
  aa = collect(a)
  QT = Vector(undef, length(aa))
  for i in 1:length(aa)
    QT[i] = list_positive_quadratic_triple(fmpz[aa[i]], Q, L, c) #List  QT_m^1 i*(aa[i], QT_m^0) for m = 2,...,n. Here v = 1
  end
  QTT = []; bbb = []; EE = Array{Array{fmpz,1},1}();
  list_b1 = []
  for j in 1:length(QT) #v=2
    b = range_ellipsoid_dim1(QT[j][1][1], QT[j][1][2], QT[j][1][3])#E(QT_2^1) = {b_1, ..., b_N}
    b1 = collect(b)
    QT1 = Vector(undef, length(b1))
    for kk in aa
      for k in 1:length(b1)
        QT1[k] = list_positive_quadratic_triple2(b1[k], QT[j]) #QT_m^2 for m = 3, ..., n
        bb = [kk, b1[k]]
        if !(bb in bbb)
        push!(bbb,bb) #bbb is a list of tuples [a,b] that satisfiy a inhomogeneous quad function for a 2 variabled quad triple
        end
      end
    end
    QT[j] = QT1
  end
  for j in 1:length(QT)
    append!(QTT, QT[j])
  end
  QTT1 = QTT;
  E = bbb
  for v in 3:n
    QTT1, E = posQuadTriple_intVectors(QTT1, E)
  end
  if !equal
    for k in E
      t = length(k)
      if (transpose(matrix(QQ,t,1,k))*Q*matrix(QQ,t,1,k))[1] + (2*transpose(matrix(QQ,t,1,k))*L)[1] + c <= 0
        push!(EE, k)
      end
    end
    if !sorting
      return EE
    else
      return sort!(EE)
    end
  else
    for k in E
      t = length(k)
      if (transpose(matrix(QQ,t,1,k))*Q*matrix(QQ,t,1,k))[1] + (2*transpose(matrix(QQ,t,1,k))*L)[1] + c == 0
        push!(EE, k)
      end
    end
    if !sorting
      return EE
    else
      return sort!(EE)
    end
  end
end

@doc Markdown.doc"""
    closest_vectors(L:ZLat, v:MatrixElem{T}, upperbound::T; equal::Bool=false, sorting::Bool=false)
                                                    -> Vector{Vector{fmpz}}

Return all vectors `x` in `L` such that `b(v-x,v-x) <= c`, where `b` is the bilinear form on `L`.
If the optional argument ``equal = true`` then it return all vectors `x` in `L` such that `b(v-x,v-x) = c`.
By default ``equal = false``. If the argument ``sorting = true``, then we get a a list of sorted vectors.
The Default value for ``sorting`` is set to ``false``.
"""
function closest_vectors(L::ZLat, v::MatrixElem{T} , upperbound::T; equal::Bool=false, sorting::Bool=false) where T <: RingElem
  epsilon = QQ(1//10)   # some number > 0, not sure how it influences performance
  d = size(v)[1]
  if isdefinite(L) == false
    throw("Zlattice is indefinite.")
  end
  if rank(L) != d
    throw("Zlattice must have the same rank as the length of the vector in the second argument.")
  end
  g1 = gram_matrix(L)
  if g1[1,1] > 0
    G1 = g1
  else
    G1 = -g1
  end
  e = matrix(QQ, 1, 1, [upperbound//3+epsilon])
  G = diagonal_matrix(G1, e)
  B = diagonal_matrix(basis_matrix(L),matrix(QQ,1,1,[1]))
  for i in 1:d
    B[end,i] = -v[i]
  end
  N = Zlattice(B,gram=G)

  delta = QQ(4//3)*upperbound + epsilon
  sv = Hecke.short_vectors(N, delta)
  cv = Array{Array{fmpz,1},1}()
  for a in sv
    a = a[1]
    if a[end] == 0
    continue
    end
    if a[end] == -1
      a = -a
    end
    x = a[1:end-1]
    push!(cv, x)
  end
  V = ambient_space(L)
  v1 = Vector{T}()
  for i in v
    push!(v1, i)
  end
  cv2 = Array{Array{fmpz,1},1}()
  if !equal
    for x in cv
      t = x - v1
      dist = inner_product(V,t,t)
      if dist <= upperbound
          push!(cv2,x)
      end
    end
    if !sorting
      return cv2
    else
      return sort!(cv2)
    end
  else
    for x in cv
      t = x - v1
      dist = inner_product(V,t,t)
      if dist == upperbound
        push!(cv2,x)
      end
    end
    if !sorting
      return cv2
    else
      return sort!(cv2)
    end
  end
end


text:
# ---
# title: Resize
# cover: resize.gif
# ---

# Set the static size of the image

using Augmentor
using ImageShow, ImageCore

img_in = testpattern(RGB, ratio=0.5)

mosaicview(
    img_in,
    augment(img_in, Resize(240, 320));
    fillvalue=colorant"white", nrow=1, npad=10
)

# ## References

#md # ```@docs
#md # Resize
#md # ```

## save covers #src
using ImageMagick #src
include(joinpath("..", "assets", "utilities.jl")) #src
cover = make_gif(testpattern(RGB, ratio=0.5), Resize(240, 320), 2) #src
ImageMagick.save("resize.gif", cover; fps=1) #src

ext: jl


text:
############################################################################################
include("unconstrained.jl")
include("constrained.jl")
include("fixed.jl")

include("distribution.jl")
include("bijector.jl")

include("simplex.jl")
include("corrmatrix.jl")
include("covmatrix.jl")

############################################################################################
# Export


text:
module CFF

using Dates
using BibInternal
using BibParser
using JSONSchema: Schema, validate
using YAML

export parse_file

const CFF_VERSIONS = Set(["1.2.0"])
const HELP_URL     = "https://github.com/citation-file-format/citation-file-format/blob/main/schema-guide.md"
const PACKAGE_ROOT = pkgdir(BibParser)
const schemas      = Dict{String,Schema}()
const current_id   = 0

"""
Simple error struct used to abort parsing if the CFF version is not supported or invalid.
"""
struct UnsupportedCFFVersion <: Any
    version::String
end

"""
    parse_file(path::String; id::String = "") ->
        Tuple{
            Union{
                BibInternal.Entry,
                UnsupportedCFFVersion,
                YAML.ParserError,
                Any
              },
            Bool
        }

Parse a CFF file located at `path`. Return a tuple with result and boolean status.
On success, a BibInternal.Entry is returned. On failure, the error object is returned.
"""
function parse_file(path; id = "")
    try
        content = YAML.load_file(path; dicttype=Dict{String, Any})
        cff_version = content["cff-version"]
        if !(cff_version in CFF_VERSIONS)
            throw(UnsupportedCFFVersion(cff_version))
        end

        if !haskey(schemas, cff_version)
            schemas[cff_version] = Schema(read(joinpath(PACKAGE_ROOT, "src", "cff", "schema-$(cff_version).json"), String))
        end

        errors = validate(content, schemas[cff_version])
        if !isnothing(errors)
            @warn "Invalid CFF file; see schema guide for details: $(HELP_URL)"
            @debug errors

            return errors, false
        end

        if haskey(content, "preferred-citation")
            content = content["preferred-citation"]
        end

        names     = add_names(content)
        access    = add_access(content)
        # no booktitle in CFF
        booktitle = ""
        date      = add_date(content)
        # editors appear only in references in CFF
        editors   = []
        eprint    = add_eprint(content)
        title     = content["title"] # always exists
        id        = isempty(id) ? generate_id(names, title, date.year, access.doi) : id
        in_       = add_in(content)
        fields    = Dict()
        type_     = add_type(content)

        BibInternal.Entry(
            access, names, booktitle, date, editors, eprint, id, in_, fields, title, type_
        ), true
    catch err
        if isa(err, YAML.ParserError)
            print("Parse error invalid YAML: ")
        end
        @warn err

        err, false
    end
end

"""
    add_names(content::Dict{String, Any}) -> Vector{BibInternal.Name}
"""
function add_names(content)
    [parse_name(author) for author in content["authors"]]
end

"""
    add_access(content::Dict{String, Any}) -> BibInternal.Access
"""
function add_access(content)
    doi = get(content, "doi", "")
    if haskey(content, "identifiers")
        # if there is no top level doi, take the one from the first identifier
        # the previous value is overwritten since this is the behavior of other implementations (cff-converter-python)
        for identifier in content["identifiers"]
            if get(identifier, "type", "") == "doi"
                doi = identifier["value"]
                break
            end
        end
    end
    url = get(content, "repository-code", "")
    if isempty(url) && haskey(content, "url")
        url = content["url"]
    end

    BibInternal.Access(doi, "", url)
end

"""
    add_date(content::Dict{String, Any}) -> BibInternal.Date
"""
function add_date(content)
    date = get(content, "date-released", "")
    if isempty(date)
        return BibInternal.Date("", "", "")
    end

    date = Date(date, dateformat"yyyy-mm-dd")

    BibInternal.Date(string(Dates.day(date)), string(Dates.month(date)), string(Dates.year(date)))
end

"""
    add_eprint(content::Dict{String, Any}) -> BibInternal.Eprint
"""
function add_eprint(content)
    # need to determine what conventions are widely used in CFFs for this field
    # this is not implemented by other converters yet
    BibInternal.Eprint("", "", "")
end

"""
    generate_id(names::Vector{BibInternal.Name}, title::String, year::String, doi::String) -> String
"""
function generate_id(names, title, year, doi)
    separator = "-"
    replace_dict = reduce(merge, map(x -> Dict([x => ""]), collect("-\$£%&(){}+!?/\\:;'\"~#")))
    replace_dict[' '] = separator
    replace_func = str -> join(map(c -> get(replace_dict, c, "$c"), collect(str)))

    dash_title = replace_func(title)[1:5]
    names_prefix = join(
        map(
            x -> replace_func(x),
            map(BibParser.name_to_string, first(names, 2))
        ),
        separator
    )

    prefix = join([names_prefix, dash_title, year], separator)
    if isempty(doi)
        current_id += 1
        "$(prefix)$(separator)$(current_id)"
    else
        "$(prefix)$(separator)$(doi)"
    end
end

"""
    add_in(content::Dict{String, Any}) -> BibInternal.In
"""
function add_in(content)
    start     = get(content, "start", "")
    finish    = get(content, "end", "")
    pages     = start == finish || isempty(finish) ? start : "{$(start)--$(finish)}"
    publisher = get(content, "publisher", Dict())
    publisher = get(publisher, "name", "")

    BibInternal.In(
        "",                                # address
        "",                                # chapter
        "",                                # edition
        "",                                # institution
        get(content, "journal", ""),
        string(get(content, "issue", "")), # number
        "",                                # organization
        pages,
        publisher,
        "",                                # school
        "",                                # series
        string(get(content, "volume", ""))
    )
end

const CFF_TO_BIBTEX_TYPES = Dict{String, String}(
    [
        "article"           => "article",
        "book"              => "book",
        "manual"            => "manual",
        "unpublished"       => "unpublished",
        "conference"        => "proceedings",
        "proceedings"       => "proceedings",
        "conference-paper"  => "proceedings",
        "magazine-article"  => "article",
        "newspaper-article" => "article",
        "pamphlet"          => "booklet"
    ]
)

"""
    add_type(content::Dict{String, Any}) -> String
"""
function add_type(content)
    get(CFF_TO_BIBTEX_TYPES, content["type"], "misc")
end

"""
    parse_name(Dict{String, Any} -> BibInternal.Name
"""
function parse_name(author)
    if haskey(author, "name")
        BibInternal.Name("", author["name"], "", "", "")
    else
        BibInternal.Name(
            get(author, "name-particle", ""),
            get(author, "family-names", ""),
            get(author, "name-suffix", ""),
            get(author, "given-names", ""),
            ""  # there is no middle name in CFF
        )
    end
end

end # module


text:
# These methods to create a Cholesky directly from the factorisation will be in Julia 1.7
# https://github.com/JuliaLang/julia/pull/39352
if VERSION < v"1.7"
    function LinearAlgebra.Cholesky(L::LowerTriangular{T}) where {T}
        return Cholesky{T,typeof(L.data)}(L.data, 'L', 0)
    end
    function LinearAlgebra.Cholesky(U::UpperTriangular{T}) where {T}
        return Cholesky{T,typeof(U.data)}(U.data, 'U', 0)
    end
end

_chol_cov(q::AbstractMvNormal) = cholesky(Symmetric(cov(q)))
_chol_cov(q::MvNormal) = cholesky(q.Σ)


text:
# --------------------------------------------------------------------------
# ACE1.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

module ACE1

using Reexport
@reexport using JuLIP

# external imports that are useful for all submodules
include("extimports.jl")

include("auxiliary.jl")
include("prototypes.jl")

# basic polynomial building blocks
include("polynomials/sphericalharmonics.jl")
include("polynomials/transforms.jl"); @reexport using ACE1.Transforms
include("polynomials/orthpolys.jl"); @reexport using ACE1.OrthPolys

# The One-particle basis is the first proper building block
include("oneparticlebasis.jl")

include("grapheval.jl")

# the permutation-invariant basis: this is a key building block
# for other bases but can also be a useful export itself
include("pibasis.jl")
include("pipot.jl")

# rotation-invariant site potentials (incl the ACE model)
include("rpi/rpi.jl")
@reexport using ACE1.RPI

# pair potentials + repulsion
include("pairpots/pair.jl");
@reexport using ACE1.PairPotentials

# lots of stuff related to random samples:
#  - random configurations
#  - random potentials
#  ...
include("random.jl")
@reexport using ACE1.Random

include("utils.jl")
@reexport using ACE1.Utils

include("utils/pure.jl")

include("compat/compat.jl")

include("export/export.jl")
include("export/export_multispecies.jl")

include("testing/testing.jl")

# - bond model
# - pure basis
# - real basis
# - regularisers
# - descriptors
# - random potentials

end # module


text:
using MultivariateStats
using Interpolations

function interp_one_obs_order_for_pca_1d(data::EOST, log_lambda_out::AbstractRange{T1}; order::Integer, obs::Integer,
      output::AbstractArray{T2,1} = zeros(eltype(data.flux),length(log_lambda_out)) ,
      z_bc::Real = zero(T1),
      blaze::Poly{T4} = Poly([one(T1)])
      ) where { EOST <: AbstractEchelleObservationSet, T1<:Real, T2<:Real, T4<:Real}
    @assert length(log_lambda_out) == length(output)
    @assert 1<=order<=num_orders(data)
    @assert 1<=obs<=num_obs(data)
    pixel_range = 1:size(data.flux,1)
    interpolator = extrapolate(
            interpolate( (log.(lambda_ssb_from_obs(data.lambda[:,order,obs],z_bc)), ),
            data.flux[:,order,obs]./polyval(blaze,pixel_range), Gridded(Linear()) ),
                  Flat() )
    output .= interpolator(log_lambda_out)
    return output
end

function interp_obs_order_for_pca_1d(data::EOST, log_lambda_out::AbstractRange{T1}; order::Integer,
      output::AbstractArray{T2,2} = zeros(eltype(data.flux),(length(log_lambda_out),num_obs(data))),
      z_bc::AbstractVector{T3} = zeros(T1,num_obs(data)),
      blaze::Poly{T4} = Poly([one(T1)])
      ) where { EOST <: AbstractEchelleObservationSet, T1<:Real, T2<:Real, T3<:Real, T4<:Real}
   obs_range = 1:num_obs(data)
   map(i->interp_one_obs_order_for_pca_1d(data,log_lambda_out,order=order,obs=i,
      output=view(output,:,i), z_bc=z_bc[i], blaze=blaze), obs_range)
   return output
end

function fit_pca_obs_order_interpolated(data::EOST, log_lambda_out::AbstractRange{T1}; order::Integer,
         z_bc::AbstractVector{T3} = zeros(T1,num_obs(data)),
         maxoutdim::Integer=10, pratio::Real=1-1e-6
         ) where { EOST <: AbstractEchelleObservationSet, T1<:Real, T3<:Real }
   obs_range = 1:num_obs(data)
   obs_interp = interp_obs_order_for_pca_1d(data, log_lambda_out, order=order, z_bc=z_bc)
   pca = fit(PCA,obs_interp, maxoutdim=10,pratio=1-1e-6)
   return pca
end

#=
data_for_pca = interp_obs_order_for_pca_1d(data,res2d[1].log_lambda_range, order=1)
pca1 = fit(PCA,data_for_pca, maxoutdim=10,pratio=1-1e-6)
principalvars(pca1)./tvar(pca1)
=#

#=
function interp_obs_for_pca_1d(data::EOST, log_lambda_out::AbstractRange{T1}; order::Integer, obs::Integer,
    output::AbstractArray{T2,1} = zeros(eltype(data.flux),(length(log_lambda_out),num_obs(data) ) ) ,
    z_bc::AbstractVector{T3} = zeros(T1,num_obs(data)),
    pixel_range::AbstractRange{T4} = 1:size(data.lambda,1),
    blaze::Poly{T5} = Poly([one(T2)])
    ) where { EOST <: AbstractEchelleObservationSet, T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real}
=#


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Viewer for the GBMCTS tree.
"""
module GBMCTSView

export TreeView, logstep!

using POMDPs
using ExprSearch, GBMCTS
using RLESUtils, TreeToJSON
using Plots
using Iterators

type TreeView
  film::Animation
  steps::Vector{JDict}
end

TreeView() = TreeView(Animation(), JDict[])

function get_name(s::StateNode)
  n_text = "n=[" * join(s.n, ",") * "]"
  Q_text = "Q=[" * join(map(x->signif(x, 4), s.Q), ",") * "]"
  na_text = "n_actions=$(length(s.action_map))"
  sp_text = "sp_visited=[" * join(map(x->!isnull(x), s.nextstates), ",")  * "]"
  expl_text = "fully_explored=$(s.explored)"
  return join([n_text, Q_text, na_text, sp_text, expl_text], "\\\\")
end

function logstep!(view::TreeView, tree::Dict{UInt64,StateNode}, startstate::State, curstate::State)

  function get_children(s::StateNode)
    valids = filter(x -> !isnull(x), s.nextstates)
    return imap(x -> ("", tree[x]), valids)
  end

  function get_color(s::StateNode)
    color = "steelblue" #default
    if all(s.n .== 0) color = "red" end #completely unvisited
    if s == tree[startstate] color = "grey" end #start state
    if s.explored color = "purple" end #fully explored
    if length(s.n) == 0 color = "black" end #terminal
    if s == tree[curstate] color = "orange" end #current state
    return color
  end

  vc = VisCalls(get_name, get_children)
  d = to_jdict(tree[startstate], vc, color=get_color)
  push!(view.steps, d)

  #animation frame
  #p = plottree(d)
  #push!(view.film, p)

  return d
end

end # module

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
# Test that input validation works where possible

@test TruncateStd(1) isa TruncateStd{<:Int}
@test TruncateStd(1.0) isa TruncateStd{<:Real}
@test_throws DomainError TruncateStd(0) isa TruncateStd{<:Int}
@test_throws DomainError TruncateStd(0.0) isa TruncateStd{<:Real}
@test_throws DomainError TruncateStd(-1) isa TruncateStd{<:Int}
@test_throws DomainError TruncateStd(-1.2) isa TruncateStd{<:Real}

@test TruncateQuantiles(0.1, 0.9) isa TruncateQuantiles{<:Real, <:Real}
@test TruncateQuantiles(0.0, 1.0) isa TruncateQuantiles{<:Real, <:Real}
@test_throws DomainError TruncateQuantiles(-0.1, 0.9)
@test_throws DomainError TruncateQuantiles(0.2, 1.9)
@test_throws DomainError TruncateQuantiles(0.3, 0.1)

@test TruncateRange(-10, 10) isa TruncateRange{<:Int, <:Int}
@test TruncateRange(-10.0, 10) isa TruncateRange{<:Real, <:Int}
@test_throws DomainError TruncateRange(5, 3)


text:
## EXPERIMENTAL
_b_ints = (511, 261632, 133955584, 68585259008, 35115652612096, 17979214137393152, 9205357638345293824) # for the last one the first bit is 0 (signbit)

@inline FLIPSIGN(x, y) = flipsign(x,y)
@inline FLIPSIGN(::Missing, ::Missing) = missing
@inline SHIFT(x,n) = x >> n
@inline SHIFT(::Missing, n) = missing

Base.@propagate_inbounds function _convert_for_countsort!(conv_x, x, bid, _b_ints, lo, hi)
    shift_val = (bid-1)*9
    for i in lo:hi
        conv_x[i] = SHIFT(abs(x[i]) & _b_ints[bid] , shift_val)
    end
    for i in lo:hi
        conv_x[i] = FLIPSIGN(conv_x[i], sign(x[i]))
    end
end

function _radixlike_sort_int_missatright!(x, conv_x, idx, cpy_idx, _b_ints, lo, hi)
    for j in 1:7
        cpy_idx .= idx
        _convert_for_countsort!(conv_x, x, j, _b_ints, lo, hi)
        _ds_sort_int_missatright_nopermx!(conv_x[idx], idx, cpy_idx, lo, hi, 1023, -512, Val(Int32))
    end
end


text:
module SparseGrids

import FastGaussQuadrature: gausshermite
import IterTools: collect, product

include("grids.jl")
include("nested.jl")

# package code goes here
export
	sparsegrid,
	tensorgrid,
	listNdq,
	combvec,
	kpn

end # module


text:
"""
    Hm... Power method should work well for symmatric tridiagonal matrix that is
    always changing its size... In fact it does, let's try inverse iterations
    instead.
"""

"""
    Left and right boundary can actually be infinite.
    Because of left and right bound, we have to restrict it to type Abstract Floats so that it's
    comparable, which is going to be used a lot for search the eigen values of the matrix.
"""
function InversePowerIterSearch(
        A::AbstractMatrix{T},
        left_bound::T,
        right_bound::T
    ) where {T<: AbstractFloat}
    # Get the mid point.
    @assert size(A, 1) == size(A, 2) "It has to be square matrix."
    @assert left_bound != right_bound "The interval must be non-degenerate."
    if left_bound > right_bound
        left_bound, right_bound = right_bound, left_bound
    end
    if left_bound == -Inf && right_bound == Inf
        error("Left and right bound cannot both be Inf, then it's unbounded. ")
    end
    v = randn(T, size(A, 1))    # previous iterations
    u = similar(v)                 # current iterations
    Au = similar(v)
    t = similar(v)              # temp
    MaxItr = 100
    while MaxItr > 0
        # Choose the right point.
        if left_bound == -Inf || right_bound == Inf
            θ = left_bound == -Inf ? 2*left_bound : 2*right_bound
            if θ == -Inf || θ == Inf
                error("Inverse Iteration Igen Search failed, bound limit reached for half opened interval.  ")
            end
        else
            θ = (left_bound + right_bound)/2
        end
        println("θ: $θ ∈ [$left_bound, $right_bound]")
        # Inverse Power iterations.
        λ = θ
        while λ == θ || norm(Au - λ*u, 1) > 1e-10
            println(λ)
            t .= u
            u .= (A - λ*I)\v
            u ./= norm(u)
            if norm(u) == NaN
                u .= v
                println("coverted to $λ")
                break
            else
                v .= t
                Au = A*u
                λ = dot(u, Au)
            end
        end
        # Interval shrink.
        if λ > left_bound && λ < right_bound

            return λ, u
        else
            if λ <= left_bound
                left_bound = θ

            else
                right_bound = θ
            end
        end
        MaxItr -= 1
    end
    if MaxItr == 0
       error("Inverse Power iterations method reached maximal iterations without convergence. ")
    end
end

# Basic Tests and shit
using LinearAlgebra
A = Diagonal(sort(randn(10)))
display(A)
λ, v = InversePowerIterSearch(A, A[1, 1], A[3, 3])
print("λ = $λ")
# λ, v = InversePowerIterSearch(A, -Inf, A[2, 2])
# print("λ = $λ")


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

runGenConfig("data/cities/small/4/gen_config.xml", overwriteOutputPath = true, doPrint = false)

@testset "route update" begin
	# Test updating of route with updateRouteToTime! function.
	# In particular, tests updating a route to a time at which it will be very close to a node.
	# Passes if assertions in routing code all pass.

	sim = initSim("data/cities/small/4/sim_config.xml")

	# shorthand
	net = sim.net
	fNodes = net.fGraph.nodes
	numFNodes = length(fNodes)

	@assert(numFNodes > length(net.rGraph.nodes)) # check that some nodes are only in fGraph

	updateRoute!(route::Route, time::Float) = if (route.recentUpdateTime <= time && route.startTime <= time) JEMSS.updateRouteToTime!(net, route, time) end

	# test route between each pair of nodes, with and without off-road travel
	startTime = 0.0
	minArcTime = minimum(fNetTravel -> minimum(fNetTravel.arcTimes), net.fNetTravels)
	offRoadSpeed = 60*24 # 60 km/hr, though any positive value will do here
	offRoadDists = [0.0, offRoadSpeed * minArcTime]
	for i = 1:numFNodes, j = 1:numFNodes, offRoadDist in offRoadDists
		route = Route()
		offRoadTime = offRoadDist / offRoadSpeed
		JEMSS.initRoute!(sim, route; startLoc = fNodes[i].location, startFNode = i, startFNodeDist = offRoadDist)
		JEMSS.changeRoute!(sim, route, highPriority, startTime + offRoadTime, fNodes[j].location, j)
		updateRoute!(route, startTime)
		updateRoute!(route, (startTime + route.startFNodeTime) / 2)
		while route.status != routeAfterEndNode
			time = route.nextFNodeTime
			dt = (time + 1) * eps(Float)
			updateRoute!(route, time - dt) # just before next node
			updateRoute!(route, time) # just at next node
			updateRoute!(route, time + dt) # just after next node
			updateRoute!(route, time + minArcTime / 2) # part way between next node and future next node
		end
		updateRoute!(route, route.endTime)
	end
	@test true
end

@testset "route distance" begin
	# Test calcRouteDistance! function.
	# Compare against a calculation that relies on each arc distance being equal to the calculated distance between the arc's nodes.
	# Assumes that shortestPathDistance function is correct (this is tested in test_network.jl).

	sim = initSim("data/cities/small/1/sim_config.xml")

	# shorthand
	net = sim.net
	fNodes = net.fGraph.nodes

	# check that arc distance is equal to distance between two nodes, measured by normDist function
	JEMSS.normDist(arc::Arc) = normDist(sim.map, fNodes[arc.fromNodeIndex].location, fNodes[arc.toNodeIndex].location)
	@assert(all(arc -> isapprox(normDist(arc), arc.distance, rtol = eps()), net.fGraph.arcs))

	# set absolute error tolerance between the two distance calculations
	atol = eps() * (sim.map.xRange * sim.map.xScale + sim.map.yRange * sim.map.yScale) # eps * upper bound on route distance

	# returns distance of route travelled up to given time,
	# assuming that arc distance is equal to calculated distance between arc's nodes
	function calcRouteDistance2(net::Network, map::Map, route::Route, time::Float)
		fNodes = net.fGraph.nodes # shorthand
		dist = 0.0
		currentLoc = JEMSS.getRouteCurrentLocation!(net, route, time)
		if route.recentFNode == nullIndex
			dist += normDist(map, route.startLoc, currentLoc) # off-road distance before startFNode
		else
			dist += normDist(map, route.startLoc, fNodes[route.startFNode].location) # off-road distance before startFNode
			dist += shortestPathDistance(net, route.travelModeIndex, route.startFNode, route.recentFNode)
			dist += normDist(map, fNodes[route.recentFNode].location, currentLoc) # off-road distance after recentFNode
		end
		return dist
	end

	dt = 0.0001 # time step
	endTime = min(1.0, sim.calls[end].arrivalTime)
	allPass = true
	for t = sim.startTime : dt : endTime
		simulate!(sim, time = t)
		for amb in sim.ambulances
			route = amb.route # shorthand
			if route.travelModeIndex != nullIndex
				dist1 = JEMSS.calcRouteDistance!(sim, route, t)
				dist2 = calcRouteDistance2(net, sim.map, route, t)
				allPass &= isapprox(dist1, dist2, atol = atol) # compare distances
			end
		end
	end
	@test allPass
end


text:
module ZeroInflation

using Distributions
using LogExpFunctions
using Random
using Requires
using UnPack

export zeroinflated

struct ZeroInflated{D1<:UnivariateDistribution,D2<:UnivariateDistribution,S<:ValueSupport} <: UnivariateDistribution{S}
    original::D1
    mix::D2
    function ZeroInflated(d::UnivariateDistribution, m::UnivariateDistribution)
        new{typeof(d),typeof(m),Distributions.value_support(typeof(d))}(d, m)
    end
end

Distributions.params(d::ZeroInflated) = tuple(Distributions.params(d.original)..., Distributions.params(d.mix)...)
Distributions.minimum(d::ZeroInflated) = min(minimum(d.original), 0)
Distributions.maximum(d::ZeroInflated) = max(maximum(d.original), 0)

function Distributions.logpdf(d::ZeroInflated, x::Real)
    @unpack original, mix = d
    if x == 0
        loglik = [
            logpdf(mix, 1),
            logpdf(mix, 0) + logpdf(original, 0)
        ]
        return logsumexp(loglik)
    else
        return logpdf(mix, 0) + logpdf(original, x)
    end
end

function Distributions.logpdf(d::ZeroInflated, x::Vector{T}) where {T<:Real}
    @unpack original, mix = d
    N_zero = sum(x .== 0)
    N_nonzero = length(x) - N_zero

    loglik = N_zero * logsumexp([
        logpdf(mix, 1),
        logpdf(mix, 0) + logpdf(original, 0)
    ])
    loglik += N_nonzero * logpdf(mix, 0)
    loglik += sum(logpdf.(original, filter(!iszero, x)))
    return loglik
end

Distributions.pdf(d::ZeroInflated, x::Real) = exp(logpdf(d, x))

function Distributions.rand(rng::Random.AbstractRNG, d::ZeroInflated)
    @unpack original, mix = d
    c = rand(rng, mix)
    return c == 1 ? zero(eltype(original)) : rand(rng, original)
end

"""
    zeroinflated(d::UnivariateDistribution, p::Real)
    zeroinflated(d::UnivariateDistribution, m::UnivariateDistribution)

Construct a zero-inflated distribution. The mixing distribution can be specified using
a mixing probability *p* or a custom mixing distribution *m* (e.g. `Bernoulli(p)``).
"""

zeroinflated(d::UnivariateDistribution, p::Real) = ZeroInflated(d, Bernoulli(p))
zeroinflated(d::UnivariateDistribution, m::UnivariateDistribution) = ZeroInflated(d, m)

end


text:
"""
This script builds the Pollen.jl documentation so that it can be loaded
by the frontend. It accepts one argument: the path where the generated
files should be stored.

    > julia docs/make.jl DIR

Use `./serve.jl` for interactive development.
"""

# Create target folder
isempty(ARGS) && error("Please pass a file path to make.jl:\n\t> julia docs/make.jl DIR ")
DIR = abspath(mkpath(ARGS[1]))

# Create Project
project = include("project.jl")

@info "Rewriting documents..."
Pollen.rewritesources!(project)

@info "Writing to disk at \"$DIR\"..."
Pollen.build(
    FileBuilder(
        JSONFormat(),
        DIR,
    ),
    project,
)


text:
function elastancefn!(system::CVSystem,n::Int64)
    # time within heart cycle
    if system.solverparams.numbeats > 0
        tp = system.t[n+1] - sum(system.heart.activation.th[1:end-1]);
    else
        tp = system.t[n+1];
    end

    # heart activation function
    # if tp < system.heart.activation.tce[system.solverparams.numbeats+1]
    #     phi = (system.heart.activation.a*sin(π*tp/
    #         system.heart.activation.tce[system.solverparams.numbeats+1])-
    #         system.heart.activation.b*sin(2*π*tp/
    #         system.heart.activation.tce[system.solverparams.numbeats+1]));
    # else
    #     phi = 0;
    # end
    g1 = (tp/system.heart.activation.tau1)^system.heart.activation.m1;
    g2 = (tp/system.heart.activation.tau2)^system.heart.activation.m2;

    # ventricular elastance
    # system.heart.lv.E[n+1] = (system.heart.lv.Emin*(1-phi) +
    #     system.heart.lv.Emax[system.solverparams.numbeats+1]*phi);
    # system.heart.rv.E[n+1] = (system.heart.rv.Emin*(1-phi) +
    #     system.heart.rv.Emax[system.solverparams.numbeats+1]*phi);
    system.heart.lv.E[n+1] = (system.heart.lv.Emax[system.solverparams.numbeats+1]-
        system.heart.lv.Emin)*system.heart.activation.k[system.solverparams.numbeats+1]*
        (g1/(1+g1))*(1/(1+g2))+system.heart.lv.Emin;
    system.heart.rv.E[n+1] = (system.heart.rv.Emax[system.solverparams.numbeats+1]-
        system.heart.rv.Emin)*system.heart.activation.k[system.solverparams.numbeats+1]*
        (g1/(1+g1))*(1/(1+g2))+system.heart.rv.Emin;
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
using Test, Colors, ColorSchemes, ColorTypes, FixedPointNumbers

monalisa = ColorScheme([
    RGB(0.05482025926320272, 0.016508952654741622, 0.019315160361063788),
    RGB(0.07508160782698388, 0.034110215845969745, 0.039708343938094984),
    RGB(0.10884977211887092, 0.033667530751245296, 0.026120424375656533),
    RGB(0.10025110094110237, 0.05342427394738222, 0.04975936729231899),
    RGB(0.11004568002009293, 0.06764950003139521, 0.07202128202310687),
    RGB(0.1520114897984492, 0.06721701384356317, 0.04758612657624729),
    RGB(0.16121466572057147, 0.10737190368841328, 0.07491505937992286),
    RGB(0.2272468746270438, 0.09450818887496519, 0.053122482545649836),
    RGB(0.24275776450376843, 0.14465569383748178, 0.09254885719488251),
    RGB(0.19832488479851235, 0.16827798680930195, 0.08146721610879516),
    RGB(0.29030547394827216, 0.1566704731433784, 0.06955958896758961),
    RGB(0.3486958875330028, 0.14413808439049522, 0.06517845643634491),
    RGB(0.2631529920611145, 0.22896210929698424, 0.1119250237167965),
    RGB(0.35775151767110114, 0.23955578484799914, 0.08566681526152695),
    RGB(0.42895506355552904, 0.19814294026377038, 0.07315576139822164),
    RGB(0.3359280058835734, 0.30177882691623686, 0.14764230985832),
    RGB(0.5168174153887967, 0.2588008525490645, 0.07751817567374263),
    RGB(0.44056726473192726, 0.3387984774995975, 0.10490250831857457),
    RGB(0.4048595970607235, 0.40823989479512734, 0.2096109034699151),
    RGB(0.619694338941659, 0.33787470822764315, 0.0871136546089913),
    RGB(0.5108290351302369, 0.41506713362977327, 0.13590312315603137),
    RGB(0.5272516131642648, 0.4706039514608196, 0.21392546020040532),
    RGB(0.5942622209175139, 0.47822315473126586, 0.14678522310513448),
    RGB(0.735266714513005, 0.4318652289706696, 0.1049661472744881),
    RGB(0.6201870982552801, 0.5227924127640037, 0.2167074150596878),
    RGB(0.6929049533440698, 0.5663098519207086, 0.18551505068207655),
    RGB(0.6814114992549445, 0.5814898147520997, 0.27039081549715527),
    RGB(0.8500397772474145, 0.5401215248181611, 0.1362117676724628),
    RGB(0.7575520588269891, 0.6334254649343621, 0.25145144950124687),
    RGB(0.8164723313500291, 0.6970150665478066, 0.32242062463720045),
    RGB(0.9330273170314637, 0.6651641943114455, 0.19865164906805746),
    RGB(0.9724409077178674, 0.7907008712807734, 0.2851364857083522)],
    "testing",
    "colors from Leonardo da Vinci's Mona Lisa")

@testset "basic tests" begin
    @test length(monalisa) == 32
    @test length(monalisa.colors) == 32
    # test that sampling schemes yield different values
    @inferred get(monalisa, 0.0)
    @test get(monalisa, 0.0) != get(monalisa, 0.5)
    @test monalisa[end] ≈ monalisa[1.0]
end

# getinverse() tests are now in ColorSchemes

@testset "conversion tests" begin
    # convert an Array{T,2} to an RGB image
    tmp = @inferred get(monalisa, rand(10, 10))
    @test typeof(tmp) == Array{ColorTypes.RGB{Float64}, 2}

    # test conversion with default clamp
    x = [0.0 1.0 ; -1.0 2.0]
    y = @inferred get(monalisa, x)
    @test y[1,1] ≈ y[2,1]
    @test y[1,2] ≈ y[2,2]

    # test conversion with symbol clamp
    y2 = @inferred get(monalisa, x, :clamp)

    @test y2[1] ≈ y[1]
    @test y2[2] ≈ y[2]
    @test y2[end] ≈ y[end]

    # test conversion with symbol extrema
    y2 = @inferred get(monalisa, x, :extrema)
    @test y2[2,1] ≈ y[1, 1]   # Minimum now becomes one edge of ColorScheme
    @test y2[2,2] ≈ y[1, 2]   # Maximum now becomes other edge of ColorScheme
    @test y2[1,1] !== y2[2, 1] # Inbetween values or now different

    # test conversion with manually supplied range
    y3 = @inferred get(monalisa, x, (-1.0, 2.0))
    @test y3[1] ≈ y2[1]
    @test y3[2] ≈ y2[2]
    @test y3[end] ≈ y2[end]

    # test empty range (#43)
    y4 = @inferred get(monalisa, 0.4, (0.0, 0.0))
    @test 0.3 < y4.r < 0.4
    @test 0.2 < y4.g < 0.25
    @test 0.1 < y4.b < 0.15

    # test gray value (#23)
    c = @inferred get(monalisa, Gray(N0f16(1.0)))
    @test typeof(c) == RGB{Float64}
    @test c.r > 0.95
    @test c.g > 0.75
    @test c.b < 0.3
    c = @inferred get(monalisa, Gray24(N0f16(1.0)))
    @test typeof(c) == RGB{Float64}
    @test c.r > 0.95
    # Booleans
    @inferred get(monalisa, false)

    @test get(monalisa, 0.0) ≈  get(monalisa, false)
    @test get(monalisa, 1.0) ≈ get(monalisa, true)
end

@testset "misc tests" begin
    # test with steplen (#17)
    r  = range(0, stop=5, length=10)
    y  = @inferred get(monalisa, r)
    y2 = @inferred get(monalisa, collect(r))

    @test isapprox(y[1], y2[1]) # fixes #61

    # test for specific value
    val = 0.2
    y   = @inferred get(monalisa, [val])
    y2  = @inferred get(monalisa, val)
    @test isapprox(y[1], y2)

    col = @inferred get(reverse(monalisa), 0.0)
    @test col.r > 0.9
    @test col.g > 0.7
    @test col.b > 0.2

    # test iteration of a colorscheme
    counter = 1
    for (n, i) in enumerate(monalisa)
        @test counter == n
        counter += 1
    end
    @test counter == 33

    # test findcolorschemes()
    @test length(findcolorscheme("rainbow")) > 0
end

# these won't error but they don't yet work correctly either
# until `get()` can handle nono-normalized color values
@testset "tests with 255-based scheme" begin
    monalisa1 = deepcopy(monalisa)
    # set all colors to 255
    for c in eachindex(monalisa1.colors)
        monalisa1.colors[c] = RGB(monalisa1.colors[c].r * 255, monalisa1.colors[c].g * 255, monalisa1.colors[c].b * 255)
    end
    @test length(monalisa1) ≈ 32
    @test length(monalisa1.colors) ≈ 32

    @test 12 < monalisa1[1].r < 14
    @test 3  < monalisa1[1].g < 5
    @test 3  < monalisa1[1].b < 6

    # test that sampling schemes yield different values
    @inferred get(monalisa1, 0.0)
    @test get(monalisa1, 0.0) != get(monalisa1, 0.5)
    @test monalisa1[end] ≈ monalisa1[1.0]
    tmp = @inferred get(monalisa1, rand(10, 10))
    @test typeof(tmp) == Array{ColorTypes.RGB{Float64}, 2}

    # test conversion with default clamp
    x = [0.0 1.0 ; -1.0 2.0]
    y = @inferred get(monalisa1, x)
    @test y[1,1] ≈ y[2,1]
    @test y[1,2] ≈ y[2,2]

    # test conversion with symbol clamp
    y2 = @inferred get(monalisa1, x, :clamp)

    @test isapprox(y2[1], y[1])
    @test isapprox(y2[2], y[2])
    @test isapprox(y2[end], y[end])

    # test conversion with symbol centered
    y2 = @inferred get(monalisa1, x, :centered)
    @test y2[2,2] ≈ y[1, 2]   # Maximum now becomes one edge of ColorScheme
    @test y2[1,1] ≈ get(monalisa1, 0.5) # Zero value is centered
    @test y2[1,2] ≈ get(monalisa1, 0.75)
    @test y2[2,1] ≈ get(monalisa1, 0.25)

    # test conversion with symbol extrema
    y2 = @inferred get(monalisa1, x, :extrema)
    @test y2[2,1] ≈ y[1, 1]   # Minimum now becomes one edge of ColorScheme
    @test y2[2,2] ≈ y[1, 2]   # Maximum now becomes other edge of ColorScheme
    @test y2[1,1] !== y2[2, 1] # Inbetween values now different

    # test conversion with manually supplied range
    y3 = @inferred get(monalisa1, x, (-1.0, 2.0))
    @test isapprox(y3[1], y2[1])
    @test isapprox(y3[2], y2[2])
    @test isapprox(y3[end], y2[end])

    # test error on unknown rangescale
    @test_throws ArgumentError get(monalisa1, x, :foo)

    # tests concatenation
    l1 = length(monalisa)
    l2 = length(monalisa * monalisa)
    @test l2 == 2l1
end

max_stars_repo_name: cormullion/ColorSchemes

max_issues_repo_name: cormullion/ColorSchemes

max_forks_repo_name: cormullion/ColorSchemes


text:
using ArgParse
using JSON3
using Sockets

function parse_commandline(args)
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--ne"
            help = "entity name"
        "--var"
            help = "variable name"
            required = true
        "--val"
            help = "variable value"
            arg_type = Float64
        "--ts"
            help = "timestamp"
            arg_type = Int
            default = Int(floor(time()))
    end

    return parse_args(args, s)
end

function script_main()
    parsed_args = parse_commandline(ARGS)
    sock = connect(3333)
    write(sock, JSON3.write(Dict(
            "ne" => parsed_args["ne"],
            "ts" => parsed_args["ts"],
            parsed_args["var"] => parsed_args["val"],
        )) * "\n")

    response = readline(sock)
    while response !== "DONE"
        msg = JSON3.read(response)
        println(msg)
        response = readline(sock)
    end
    close(sock)
end

script_main()


text:
using CUDA
using BenchmarkTools
using GemmKernels
using Statistics
using Printf
using LinearAlgebra

function summarise_bench(b)
    mu = mean(b.times)
    sigma = std(b.times)

    if mu < 1e3
        mu, sigma, unit = mu, sigma, "ns"
    elseif mu < 1e6
        mu, sigma, unit = mu / 1e3, sigma / 1e3, "μs"
    elseif mu < 1e9
        mu, sigma, unit = mu / 1e6, sigma / 1e6, "ms"
    else
        mu, sigma, unit = mu / 1e9, sigma / 1e9, "s"
    end

    println(string(@sprintf("(%.2f ± %.2f)", mu, sigma), " ", unit))
end

M = 4096
N = 4096
K = 4096

transpose_a = false
transpose_b = false

a_h = rand(Float16, M);
b_h = rand(Float16, (K, N)) / sqrt(Float16(K))
c_h = rand(Float32, (M, N))

# Transpose input if necessary
a_h = transpose_a ? transpose(a_h) : a_h
b_h = transpose_b ? transpose(b_h) : b_h

a_gk     = CuArray(a_h)
a_cublas = CuArray(Array(Diagonal(a_h)))
b        = CuArray(b_h)
c        = CuArray(c_h)

function bench_cublas(a, b, c, M, N, K, transpose_a, transpose_b, num_iterations = 100)
    CUBLAS.cublasSetMathMode(CUBLAS.handle(), CUBLAS.CUBLAS_TENSOR_OP_MATH)

    CUDA.@sync begin
        for i = 1 : num_iterations
            CUDA.CUBLAS.gemmEx!(!transpose_a ? 'N' : 'T', !transpose_b ? 'N' : 'T', Float32(1), a, b, Float32(1), c)
        end
    end
end

function bench_gemmkernels(a, b, c, M, N, K, transpose_a, transpose_b, num_iterations = 100)
    conf = GemmKernels.get_config(
                                  gemm_shape = (M = M, N = N, K = K),
                                  operator = Operator.WMMAOp{16, 16, 16},
                                  global_a_layout = Layout.Diagonal{Float16},
                                  global_b_layout = transpose_b ? Layout.AlignedRowMajor{Float16} : Layout.AlignedColMajor{Float16},

                                  global_c_layout = Layout.AlignedColMajor{Float32},
                                  global_d_layout = Layout.AlignedColMajor{Float32},

                                  shared_a_layout = Layout.Padded{Layout.AlignedColMajor{Float16}, 8},

                                  is_a_col_major = !transpose_a,
                                  is_b_col_major = !transpose_b,
                                 )
    CUDA.@sync begin
        for i = 1 : num_iterations
            GemmKernels.matmul(a, b, c, c, conf;)
        end
    end
end


text:
using Requests
# import Requests: get

const BOSS_TOKEN = ENV["INTERN_TOKEN"]

headers = Dict("Authorization"          => "Token $(BOSS_TOKEN)",
                        "content-type"  => "application/blosc",
                        "Accept"        => "application/blosc")

resp = Requests.get(URI("https://api.theboss.io/v0.8/cutout/team2_waypoint"*
                "/pinky10/em/0/10000:10005/10000:10006/0:1/");
                headers = headers)

fieldnames(resp)
data = resp.data

using Blosc
Blosc.set_num_threads(div(Sys.CPU_CORES,2))
data = Blosc.decompress(UInt8, data)

# data = reshape(data,(5,6,4))

length(data)

@show resp


text:
function htmlheader(data)
  header = "
\n"
  if typeof(data) <: DataFrames.DataFrame
    for field in names(data)
      header = string(header,"",field,"	\n")
    end
  else
    for j in 1:size(data)[2]
      header = string(header,"",j,"	\n")
    end
  end
  header = string(header,"
\n")
  return header
end

function htmlbody(data)
  body = string()
  for i in 1:size(data)[1]
    body = string(body,"
\n")
    for j in 1:size(data)[2]
      body = string(body,"",data[i,j],"	\n")
    end
    body = string(body,"
\n")
  end
  return body
end

function writehtml(data)

preamble = "
<!DOCTYPE html>\n
\n
\n
\n
  <meta charset=\"UTF-8\">\n
  <title> mytable
  </title>\n
<\head>\n
\n
\n
\n

\n"

ending = "
\n
  \n
  \n
  \n
  \n
  "

html = string(preamble,htmlheader(data),htmlbody(data),ending);

return html

end


text:
### single layer initWB

#TODO create initializers dependent for each method

"""
    initialize W and B for layer with inputs of size of (nl_1) and layer size
        of (nl)

    returns:
        W: of size of (nl, nl_1)
"""
function initWB!(
    cLayer::FCLayer,
    p::Type{T} = Float64::Type{Float64};
    He = true,
    coef = 0.01,
    zro = false,
) where {T}

    s = (cLayer.channels, cLayer.prevLayer.channels)
    if He
        coef = sqrt(2 / s[end])
    end
    if zro
        W = zeros(T, s...)
    else
        W = T.(randn(s...) .* coef)
    end
    B = zeros(T, (s[1], 1))

    cLayer.W, cLayer.B = W, B
    cLayer.dW = zeros(T, s...)
    cLayer.dB = deepcopy(B)
    return nothing
end #initWB

function initWB!(
    cLayer::Activation,
    p::Type{T} = Float64::Type{Float64};
    He = true,
    coef = 0.01,
    zro = false,
) where {T}

    return nothing
end

function initWB!(
    cLayer::Flatten,
    p::Type{T} = Float64::Type{Float64};
    He = true,
    coef = 0.01,
    zro = false,
) where {T}

    return nothing
end

function initWB!(
    cLayer::ConcatLayer,
    p::Type{T} = Float64::Type{Float64};
    He = true,
    coef = 0.01,
    zro = false,
) where {T}

    return nothing
end

function initWB!(
    cLayer::Input,
    p::Type{T} = Float64::Type{Float64};
    He = true,
    coef = 0.01,
    zro = false,
) where {T}

    # cLayer.inputS = cLayer.outputS = size(cLayer.A)
    return nothing
end

function initWB!(
    cLayer::BatchNorm,
    p::Type{T} = Float64::Type{Float64};
    He = true,
    coef = 0.01,
    zro = false,
) where {T}

    cn = cLayer.prevLayer.channels

    if He
        coef = sqrt(2 / cn)
    end
    N = length(cLayer.prevLayer.inputS)
    normDim = cLayer.dim
    #bring the batch size into front
    S = (cLayer.prevLayer.outputS[end], cLayer.prevLayer.outputS[1:end-1]...)
    paramS = Array{Integer,1}([1])
    for i=2:N
        if S[i] < 1 || (i-1) <= normDim
            push!(paramS, 1)
        else
            push!(paramS, S[i])
        end #if S[i] < 1
    end #for

    if !zro
        W = T.(randn(paramS...) .* coef)
    else
        W = zeros(T, paramS...)
    end

    if !(size(cLayer.W) == size(W))
        cLayer.W = W
        cLayer.dW = zeros(T, paramS...)
        cLayer.B = zeros(T, paramS...)
        cLayer.dB = zeros(T, paramS...)

    end #if !(size(cLayer.W) == size(W))
    return nothing
end #function initWB!(cLayer::BatchNorm

export initWB!

### deepInitWB!

"""
    initialize W's and B's using

    inputs:
        X := is the input of the neural Network
        outLayer := is the output Layer or the current layer
                    of initialization
        cnt := is a counter to determinde the current step
                and its an internal variable

        kwargs:
            He := is a true/false array, whether to use the He **et al.** initialization
                    or not

            coef := when not using He **et al.** initialization use this coef
                    to multiply with the random numbers initialization
            zro := true/false variable whether to initialize W with zeros or not

"""
function deepInitWB!(
    outLayer::Layer,
    cnt = -1;
    He = true,
    coef = 0.01,
    zro = false,
    dtype = Float64,
)
    if cnt < 0
        cnt = outLayer.forwCount + 1
    end

    if dtype == nothing
        T = eltype(X)
    else
        T = dtype
    end

    prevLayer = outLayer.prevLayer
    forwCount = outLayer.forwCount
    if outLayer isa Input
        if forwCount < cnt
            outLayer.forwCount += 1
            initWB!(outLayer, T; He = He, coef = coef, zro = zro)
        end #if forwCount < cnt
    elseif isa(outLayer, MILayer)
        if forwCount < cnt
            outLayer.forwCount += 1
            for prevLayer in outLayer.prevLayer
                deepInitWB!(
                    prevLayer,
                    cnt;
                    He = He,
                    coef = coef,
                    zro = zro,
                    dtype = dtype,
                )
            end #for prevLayer in outLayer.prevLayer

        end #if forwCount < cnt

    else #if prevLayer == nothing
        if forwCount < cnt
            outLayer.forwCount += 1
            deepInitWB!(
                prevLayer,
                cnt;
                He = He,
                coef = coef,
                zro = zro,
                dtype = dtype,
            )

            initWB!(outLayer, T; He = He, coef = coef, zro = zro)
        end #if forwCount < cnt

    end #if prevLayer == nothing

    return nothing
end #deepInitWB

export deepInitWB!

### single layer initVS

function initVS!(
    cLayer::FCLayer,
    optimizer::Symbol
    )

    cLayer.V[:dw] = deepcopy(cLayer.dW)
    cLayer.V[:db] = deepcopy(cLayer.dB)
    if optimizer == :adam
        cLayer.S = deepcopy(cLayer.V)
    end #if optimizer == :adam

    return nothing
end #function initVS!

function initVS!(
    cLayer::IoA,
    optimizer::Symbol
    ) where {IoA <: Union{Input, Activation, Flatten, ConcatLayer}}

    return nothing
end #

function initVS!(
    cLayer::BatchNorm,
    optimizer::Symbol,
    )

    T = typeof(cLayer.W)
    if (! haskey(cLayer.V, :dw)) || (size(cLayer.V[:dw]) != size(cLayer.dW))
        cLayer.V = Dict(:dw=>deepcopy(cLayer.dW), :db=>deepcopy(cLayer.dB))
        cLayer.S = deepcopy(cLayer.V)
    end

    return nothing
end

export initVS!

### deepInitVS!

function deepInitVS!(outLayer::Layer, optimizer::Symbol, cnt::Integer = -1)

    if cnt < 0
        cnt = outLayer.forwCount + 1
    end

    prevLayer = outLayer.prevLayer

    if optimizer == :adam || optimizer == :momentum
        if outLayer isa Input
            if outLayer.forwCount < cnt
                outLayer.forwCount += 1
                initVS!(outLayer, optimizer)
            end #if outLayer.forwCount < cnt
        elseif isa(outLayer, MILayer)
            if outLayer.forwCount < cnt
                outLayer.forwCount += 1
                for prevLayer in outLayer.prevLayer
                    deepInitVS!(prevLayer, optimizer, cnt)
                end #for
            end #if outLayer.forwCount < cnt
        else #if prevLayer == nothing
            if outLayer.forwCount < cnt
                outLayer.forwCount += 1
                deepInitVS!(prevLayer, optimizer, cnt)
                initVS!(outLayer, optimizer)
            end #if outLayer.forwCount < cnt
        end #if prevLayer == nothing

    end #if optimizer==:adam || optimizer==:momentum

    return nothing
end #function deepInitVS!

export deepInitVS!


text:
#Tensors in Tucker format + functions

export ttensor, randttensor
export coresize, cp_als, display, ewprod, full, hadcten, hosvd, hosvd1, hosvd2, hosvd3, hosvd4, innerprod, isequal, lanczos, lanczos_tridiag, mhadtv, minus, mrank
export msvdvals, mtimes, mttkrp, ndims, nrank, nvecs, permutedims, plus, randrange, randsvd, reorth, reorth!, size, tenmat, ttm, ttv, uminus, norm

"""
    ttensor(cten,fmat)

Tensor in Tucker format defined by its core tensor and factor matrices.
For ttensor X, X.isorth=true if factor matrices are othonormal.
"""
mutable struct ttensor
	cten::Array{<:Number}
	fmat::MatrixCell
	isorth::Bool
	function ttensor(cten::Array{<:Number},fmat::MatrixCell,isorth::Bool)
		for A in fmat
			if norm(A'*A-eye(size(A,2)))>(size(A,1)^2)*eps()
				isorth=false
			end
		end
		new(cten,fmat,isorth)
	end
end
#ttensor(cten::Array{<:Number},fmat::MatrixCell,isorth::Bool)=ttensor(cten,fmat,isorth)
ttensor(cten::Array{<:Number},fmat::MatrixCell)=ttensor(cten,fmat,true)
ttensor(cten::Array{<:Number},mat::Matrix{<:Number},isorth::Bool)=ttensor(cten,collect(mat),isorth)
ttensor(cten::Array{<:Number},mat::Matrix{<:Number})=ttensor(cten,collect(mat),true)
#If array of matrices isn't defined as MatrixCell, but as M=[M1,M2,...,Mn]:
ttensor(cten::Array{<:Number},fmat::Array{Matrix{T},1},isorth::Bool) where T<:Number=ttensor(cten,MatrixCell(fmat),isorth)
ttensor(cten::Array{<:Number},fmat::Array{Matrix{T},1}) where T<:Number=ttensor(cten,MatrixCell(fmat),true)

"""
    randttensor(I::Vector,R::Vector)
    randttensor(I::Integer,R::Integer,N::Integer)

Create random ttensor of size I and multilinear rank R, or of order N and size I × ⋯ × I and mulilinear rank (R,...,R).
"""
function randttensor(sz::AbstractVector{<:Integer},R::AbstractVector{<:Integer})
  @assert(size(sz)==size(R),"Size and rank should be of same length.")
  cten=randn(tuple(R...)) #create radnom core tensor
  fmat=[randn(sz[n],R[n]) for n=1:length(sz)] #create random factor matrices
  ttensor(cten,fmat)
end
randttensor(sz::Integer,R::Integer,N::Integer)=randttensor(repeat([sz],N),repeat([R],N));
#For input defined as tuples or nx1 matrices - randttensor(([I,I,I],[R,R,R]))
function randttensor(arg...)
  randttensor([arg[1]...],[arg[2]...])
end

"""
    coresize(X)

Size of core tensor of a ttensor.
"""
function coresize(X::ttensor)
  size(X.cten)
end

#Compute a CP decomposition with R components of a tensor X. **Documentation in tensor.jl.
function cp_als(X::ttensor,R::Integer;init="rand",tol=1e-4,maxit=1000,dimorder=[])
    N=ndims(X)
    nr=norm(X)
    K=ktensor
    if length(dimorder) == 0
        dimorder=collect(1:N)
    end
    fmat=MatrixCell(undef,N)
    if isa(init,Vector) || isa(init,MatrixCell)
        @assert(length(init)==N,"Wrong number of initial matrices.")
        for n in dimorder[2:end]
            @assert(size(init[n])==(size(X,n),R),"$(n)-th initial matrix is of wrong size.")
            fmat[n]=init[n]
        end
    elseif init=="rand"
        [fmat[n]=rand(size(X,n),R) for n in dimorder[2:end]]
    elseif init=="eigs" || init=="nvecs"
        [fmat[n]=nvecs(X,n,R) for n in dimorder[2:end]]
    else
        error("Initialization method wrong.")
    end
    G = zeros(R,R,N); #initalize gramians
    [G[:,:,n]=fmat[n]'*fmat[n] for n in dimorder[2:end]]
    fit=0
    for k=1:maxit
        fitold=fit
        lambda=[]
        for n in dimorder
            fmat[n]=mttkrp(X,fmat,n)
            W=reshape(prod(G[:,:,setdiff(collect(1:N),n)],dims=3),Val(2))
            fmat[n]=fmat[n]/W
            if k == 1
                lambda = sqrt.(sum(fmat[n].^2,dims=1))[:] #2-norm
            else
                lambda = maximum(maximum(abs.(fmat[n]),dims=1),dims=1)[:] #max-norm
            end
            fmat[n] = fmat[n]./lambda'
            G[:,:,n] = fmat[n]'*fmat[n]
        end
        K=ktensor(lambda,fmat)
        if nr==0
            fit=norm(K)^2-2*innerprod(X,K)
        else
            nr_res=sqrt.(abs.(nr^2+norm(K)^2 .-2*innerprod(X,K)))
            fir=1 .-nr_res/nr
        end
        fitchange=abs.(fitold-fit)
        if k>1 && fitchange<tol
            break
        end
    end
    arrange!(K)
    fixsigns!(K)
    K
end

"""
---
TensorToolbox:

    display(X::ttensor[,name])
    display(X::ktensor[,name])
    display(X::htensor[,name])

Displays a tensor X of a given name.
"""
function display(X::ttensor,name="ttensor")
  print("Tucker tensor of size ",size(X)," with core tensor of size ",coresize(X))
  if X.isorth == true
    print(" with orthonormal factor matrices")
  end
  print(":\n")
  println("$name.cten: ")
  show(stdout, "text/plain", X.cten)
  for n=1:ndims(X)
      println("\n\n$name.fmat[$n]:")
      show(stdout, "text/plain", X.fmat[n])
  end
end

"""
---
TensorToolbox:

    full(X::ttensor)
    full(X::ktensor)
    full(X::htensor)

Make full tensor out of a decomposed tensor.
"""
function full(X::ttensor)
  ttm(X.cten,X.fmat)
end

"""
   ewprod(X::ttensor,Y::ttensor)
   ewprod(X::TTtensor,Y::TTtensor)

Element-wise product of two ttensors/TTtensors.
"""
function ewprod(X1::ttensor,X2::ttensor)
  @assert(size(X1) == size(X2))
  fmat=MatrixCell(undef,ndims(X1)) #initilize factor matrix
  n=1
  for (A1,A2) in zip(X1.fmat,X2.fmat)
      fmat[n]=khatrirao(A1,A2,'t')
      n+=1
	 end
  cten=tkron(X1.cten,X2.cten) #Kronecker product of core tensors
  ttensor(cten,fmat)
end

"""
    hadcten(X,Y,fmat)

Core tensor of Hadamard product of two ttensors with given factor matrices.
"""
function hadcten(X1::ttensor,X2::ttensor,fmat::MatrixCell)
  N=ndims(X1);
  C=MatrixCell(undef,N)
  for n=1:N
    C[n]=fmat[n]'*khatrirao(X1.fmat[n],X2.fmat[n],'t');
  end
  cten=krontm(X1.cten,X2.cten,C)
end
hadcten(X1::ttensor,X2::ttensor,fmat::Array{AbstractMatrix{T}}) where {T<:Number}=hadcten(X1,X2,MatrixCell(fmat))
#HOSVD for a ttensor. **Documentation in tensor.jl
function hosvd(X::ttensor;method="svd",reqrank=[],eps_abs=[],eps_rel=[])
  F=hosvd(X.cten,method=method,reqrank=reqrank,eps_abs=eps_abs,eps_rel=eps_rel)
  fmat=MatrixCell(undef,ndims(X))  #fmat=Array{Matrix,ndims(X)}
  [fmat[n]=X.fmat[n]*F.fmat[n] for n=1:ndims(X)]
  reorth(ttensor(F.cten,fmat))
end

"""
   hosvd1(X,Y; <keyword arguments>)

Hadamard product of ttensors X and Y as ttensor. Creates the product and calls hosvd.
See also: hosvd2, hosvd3, hosvd4.

## Arguments:
- `method` ∈ {"svd","lanczos","randsvd"}. Method for SVD. Default: "randsvd".
- `reqrank::Vector`: Requested mutlilinear rank. Optional.
- `eps_abs::Number/Vector`: Drop singular values (of mode-n matricization) below eps_abs. Optional.
- `eps_rel::Number/Vector`: Drop singular values (of mode-n matricization) below eps_rel*sigma_1. Optional.
- `p::Integer`: Oversampling parameter. Defaul p=10.
"""
function hosvd1(X1::ttensor,X2::ttensor;method="randsvd",reqrank=[],eps_abs=[],eps_rel=[],p=10)
  Xprod=full(X1).*full(X2);
  hosvd(Xprod,method=method,reqrank=reqrank,eps_abs=eps_abs,eps_rel=eps_rel,p=p)
end

"""
   hosvd2(X,Y; <keyword arguments>)

Hadamard product of ttensors X and Y as ttensor. Orthogonalizes factor matrices from structure and calls hosvd on updated core tensor.
See also: hosvd1, hosvd3, hosvd4.

## Arguments:
- `method` ∈ {"svd","lanczos","randsvd"} Method for SVD. Default: "randsvd".
- `reqrank::Vector`: Requested mutlilinear rank. Optional.
- `eps_abs::Number/Vector`: Drop singular values (of mode-n matricization) below eps_abs. Optional.
- `eps_rel::Number/Vector`: Drop singular values (of mode-n matricization) below eps_rel*sigma_1. Optional.
- `p::Integer`: Oversampling parameter. Defaul p=10.
"""
function hosvd2(X1::ttensor,X2::ttensor;method="randsvd",reqrank=[],eps_abs=[],eps_rel=[],p=10)
  @assert(size(X1) == size(X2))
  N=ndims(X1)
  Q=MatrixCell(undef,N)
  R=MatrixCell(undef,N)
  n=1
  for (A1,A2) in zip(X1.fmat,X2.fmat)
    Ahad=khatrirao(A1,A2,'t')
    Qm,R[n]=qr(Ahad)
	Q[n]=Matrix(Qm)
    n+=1
	 end
  X=hosvd(krontm(X1.cten,X2.cten,R),method=method,reqrank=reqrank,eps_abs=eps_abs,eps_rel=eps_rel,p=p);
  cten=X.cten;
  fmat=MatrixCell(undef,N)
  [fmat[n]=Q[n]*X.fmat[n] for n=1:N];
  ttensor(cten,fmat)
end

"""
   hosvd3(X,Y; <keyword arguments>)

Hadamard product of ttensors X and Y as ttensor. Structure exploiting, works with (X ∗ Y)ₙ(X ∗ Y)ₙᵀ matrices.
See also: hosvd1, hosvd2, hosvd4.

## Arguments:
- `method` ∈ {"lanczos","randsvd"} Structure exploiting method for SVD. Default: "lanczos".
- `reqrank::Vector`: Requested mutlilinear rank. Optional.
- `variant` ∈ {'A','B'} Variant of multiplication (X ∗ Y)ₙ(X ∗ Y)ₙᵀ. Default: 'B'.
- `eps_abs::Number/Vector`: Drop singular values (of mode-n matricization) below eps_abs. Optional.
- `eps_rel::Number/Vector`: Drop singular values (of mode-n matricization) below eps_rel*sigma_1. Optional.
- `p::Integer`: Oversampling parameter. Defaul p=10.
"""
function hosvd3(X1::ttensor,X2::ttensor;method="lanczos",reqrank=[],variant='B',eps_abs=[],eps_rel=[],p=10)
  @assert(size(X1) == size(X2))
  N=ndims(X1)
  Ahad=MatrixCell(undef,N) #initilize factor matrices
  if method != "lanczos" && method != "randsvd"
    error("Incorect method name.")
  end
  reqrank=check_vector_input(reqrank,N,0);
  eps_abs=check_vector_input(eps_abs,N,1e-8);
  eps_rel=check_vector_input(eps_rel,N,0);
  @assert(N==length(reqrank),"Dimensions mismatch.")
  for n=1:N
    if method=="lanczos"
      Ahad[n],S=lanczos(X1,X2,n,variant=variant,reqrank=reqrank[n],tol=eps_abs[n],p=p)
    elseif method=="randsvd"
      Ahad[n],S=randsvd(X1,X2,n,variant=variant,reqrank=reqrank[n],tol=eps_abs[n],p=p)
    end
    if reqrank[n] == 0
      eps_rel[n] != 0 ?  tol=eps_rel[n]*S[1] : tol=eps_abs[n];
      K=findall(x-> x>tol ? true : false,S)
      Ahad[n]=Ahad[n][:,K];
    end
  end
  core=hadcten(X1,X2,Ahad)
  ttensor(core,Ahad)
end

"""
   hosvd4(X,Y; <keyword arguments>)

Hadamard product of ttensors X and Y as ttensor. Uses rank-1 randomized algorithm for finding range of (X ∗ Y)ₙ.
If reqrank defined, calls additonal hosvd on updated core tensor.
See also: hosvd1, hosvd2, hosvd3.

## Arguments:
- `method` ∈ {"svd","lanczos","randsvd"} Method for SVD. Default: "svd".
- `reqrank::Vector`: Requested mutlilinear rank. Optional.
- `eps_abs::Number/Vector`: Drop singular values (of mode-n matricization) below eps_abs. Optional.
- `eps_rel::Number/Vector`: Drop singular values (of mode-n matricization) below eps_rel*sigma_1. Optional.
- `p::Integer`: Oversampling parameter. Defaul p=10.
"""
function hosvd4(X1::ttensor,X2::ttensor;method="svd",reqrank=[],eps_abs=[],eps_rel=[],p=10)
  @assert(size(X1) == size(X2))
  N=ndims(X1)
  reqrank=check_vector_input(reqrank,N,0)
  eps_abs=check_vector_input(eps_abs,N,1e-8)
  eps_rel=check_vector_input(eps_rel,N,0)
  Q=MatrixCell(undef,N) #range approximation of tenmat(X1.*X2,n)
  #KR=MatrixCell(undef,N); #transpose Khatri-Rao product of X1.fmat and X2.fmat
  fmat=MatrixCell(undef,N)
  #[KR[n]=khatrirao(X1.fmat[n],X2.fmat[n],'t') for n=1:N]
  KR=khatrirao(X1.fmat,X2.fmat,'t')
  for n=1:N
    Q[n]=randrange(X1.cten,X2.cten,KR,n,reqrank=reqrank[n],tol=eps_abs[n],p=p)
  end
  [fmat[n]=Q[n]'*KR[n] for n=1:N]
  H=krontm(X1.cten,X2.cten,fmat)
  if length(reqrank) != 0 #fixed-rank problem
    Htucker=hosvd(H,reqrank=reqrank,method=method,eps_abs=eps_abs,eps_rel=eps_rel)
    [fmat[n]=Q[n]*Htucker.fmat[n] for n=1:N]
    return ttensor(Htucker.cten,fmat)
  else
    return ttensor(H,Q)
  end
end

#Inner product of two ttensors. **Documentation in tensor.jl.
function innerprod(X1::ttensor,X2::ttensor)
	@assert size(X1) == size(X2)
	if prod(size(X1.cten)) > prod(size(X2.cten))
		innerprod(X2,X1)
	else
    N=ndims(X1)
    fmat=MatrixCell(undef,N)
    [fmat[n]=X1.fmat[n]'*X2.fmat[n] for n=1:N]
		innerprod(X1.cten,ttm(X2.cten,fmat))
	end
end

"""
    isequal(X::ttensor,Y::ttensor)
    isequal(X::ktensor,Y::ktensor)
    isequal(X::ktensor,Y::htensor)

Two tensors in decomposed format are equal if they have equal components. Same as: X==Y.
"""
function isequal(X1::ttensor,X2::ttensor)
  if (X1.cten == X2.cten) && (X1.fmat == X2.fmat)
    true
  else
    false
  end
end
==(X1::ttensor,X2::ttensor) =isequal(X1,X2)

"""
    lanczos(X,Y,n; <keyword arguments>)

Structure exploiting Lanczos based SVD - computes left singular vectors and singular values of n-mode matricization (X ∗ Y)ₙ.
Works with matrix (X ∗ Y)ₙ(X ∗ Y)ₙᵀ.

## Arguments:
- `reqrank::Integer`: Requested rank. Optional.
- `variant` ∈ {'A','B'} Variant of multiplication (X ∗ Y)ₙ(X ∗ Y)ₙᵀ. Default: 'B'.
- `tol::Number/Vector`: Tolerance. Default: 1e-8.
- `maxit`: Maximal number of iterations. Default: 1000.
- `p::Integer`: Oversampling parameter. Defaul p=10.
"""
function lanczos(X1::ttensor,X2::ttensor,mode::Integer;reqrank=0,variant='B',tol=1e-8,maxit=1000,p=10)
  @assert(size(X1)==size(X2),"Dimensions mismatch")
  Q,T=lanczos_tridiag(X1,X2,mode,reqrank=reqrank,variant=variant,tol=tol,maxit=maxit,p=p)
  E=eigen(T,tol,Inf);
  U=E.vectors[:,end:-1:1];
  S=sqrt.(abs.(E.values[end:-1:1]));
  if reqrank!=0
    if reqrank > size(U,2)
      @warn "Required rank for mode $mode exceeds actual rank, the resulting rank will be ",size(U,2),". Try changing tolerance.";
    else
      U=U[:,1:reqrank];
      S=S[1:reqrank];
    end
  end
  U=Q*U;
  U,S
end

"""
    lanczos_tridiag(X,Y,n; <keyword arguments>)

Structure exploiting Lanczos tridiagonalization algorithm for n-mode matricization A=(X ∗ Y)ₙ(X ∗ Y)ₙᵀ.
Returns orthonormal Q and symmetric tridiagonal T such that A≈Q*T*Q'.

## Arguments:
- `reqrank::Integer`: Requested rank. Optional.
- `variant` ∈ {'A','B'} Variant of multiplication (X ∗ Y)ₙ(X ∗ Y)ₙᵀ. Default: 'B'.
- `tol::Number/Vector`: Tolerance. Default: 1e-8.
- `maxit`: Maximal number of iterations. Default: 1000.
- `p::Integer`: Oversampling parameter. Defaul p=10.
"""
function lanczos_tridiag(X1::ttensor,X2::ttensor,mode::Integer;reqrank=0,variant='B',tol=1e-8,maxit=1000,p=10)
  sz=size(X1)
  m=sz[mode]
  n=prod(deleteat!(copy([sz...]),mode))
  K=min(m,maxit);
  if reqrank!=0
    K=min(K,reqrank+p);
  end
  α=zeros(K)
  β=zeros(K)
  v=randn(m);
  q=v/norm(v);
  Q=zeros(m,1)
  Q[:,1]=q;
  k=0; #needed if stopping criterion is met
  for k=1:K
    r=mhadtv(X1,X2,Q[:,k],mode,variant=variant)
    α[k]=dot(r,Q[:,k])
    r=r-α[k]*Q[:,k]
    [r=r-Q*(Q'*r) for i=1:3]
    β[k]=norm(r)
    if β[k] < tol
        K=k
      break
    end
    if k!=K
      Q=[Q r/β[k]]
    end
  end
  T=SymTridiagonal(α[1:K], β[1:K-1])
  Q,T
end

"""
    mhadtv(X,Y,v,n,t='b';variant='B')

Mode-n matricized Hadamard product of ttensors X and Y times vector v.
- `t='b'`:  (X ∗ Y)ₙ(X ∗ Y)ₙᵀv. Variant of multiplication (X ∗ Y)ₙ(X ∗ Y)ₙᵀ can be specified to 'A' or 'B'.
- `t='n'`:  (X ∗ Y)ₙᵀv.
- `t='t'`:  (X ∗ Y)ₙᵀv.
If v is a matrix, multiply column by column.
"""
function mhadtv(X1::ttensor,X2::ttensor,v::AbstractVector{<:Number},n::Integer,t='b';variant='B')
  @assert(size(X1)==size(X2),"Dimensions mismatch")
  sz=size(X1)
  r1=coresize(X1)
  r2=coresize(X2)
  R=[r1...].*[r2...]
  N=setdiff(1:ndims(X1),n) #all indices but n

  #X1=G₁ ×₁ A₁ ×₂ ... ×ₗ Aₗ
  #X2=G₂ ×₁ B₁ ×₂ ... ×ₗ Bₗ

  if t=='t'
    @assert(length(v) == sz[n],"Vector v is of inappropriate size.")
    w1=krtv(copy(X1.fmat[n]'),copy(X2.fmat[n]'),v); #w1=(Aₖ' ⨀ Bₖ')*v
    W1=mkrontv(X1.cten,X2.cten,w1,n,'t') #W1=tenmat(G₁ ⨂ G₂,n)'*w1
    for k in N
      #W1=copy(reshape(W1,R[k],round.(Int,prod(size(W1))/R[k])))
      W1=copy(reshape(W1,R[k],:))
      W2=tkrtv(X1.fmat[k],X2.fmat[k],W1) #vec(W2)=(Aₖ ⨀' Bₖ)*vec(W1)
      W1=copy(W2')
    end
    vec(W1)
  elseif t=='n'
    @assert(length(v) == prod(deleteat!(copy(collect(sz)),n)),"Vector v is of inappropriate size.")
    W1=v
    for k in N
      W1=copy(reshape(W1,sz[k],:))
      W2=krtv(copy(X1.fmat[k]'),copy(X2.fmat[k]'),W1) #W2=(Aₖ' ⨀ Bₖ')*W1
      W1=copy(W2')
    end
    W2=mkrontv(X1.cten,X2.cten,vec(W1),n) #W1=tenmat(G₁ ⨂ G₂),n)*vec(W2)
    tkrtv(X1.fmat[n],X2.fmat[n],W2) #(Aₖ ⨀' Bₖ)*W2
  elseif t=='b'
    @assert(length(v) == sz[n],"Vector v is of inappropriate size.")
    if variant == 'A'    #use when prod(I[N])-prod(R[N]) < 0
      mhadtv(X1,X2,mhadtv(X1,X2,v,n,'t'),n,'n')
    elseif variant == 'B'
      w1=krtv(copy(X1.fmat[n]'),copy(X2.fmat[n]'),v); #w1=(Aₖ' ⨀ Bₖ')*v
      W1=mkrontv(X1.cten,X2.cten,w1,n,'t') #W1=tenmat(G₁ ⨂ G₂,n)'*w1
      for k in N
        #W1=copy(reshape(W1,R[k],round.(Int,prod(size(W1))/R[k])))
        W1=copy(reshape(W1,R[k],:))
        W2=tkrtv(X1.fmat[k],X2.fmat[k],W1)
        W1=copy(krtv(copy(X1.fmat[k]'),copy(X2.fmat[k]'),W2)')
      end
      W2=mkrontv(X1.cten,X2.cten,vec(W1),n) #W2=tenmat(G₁ ⨂ G₂),n)*vec(W1)
      tkrtv(X1.fmat[n],X2.fmat[n],W2) #(Aₖ' ⨀ Bₖ')'*W2
    else
      error("Variant should be either 'A' or 'B'.")
    end
  end
end
function mhadtv(X1::ttensor,X2::ttensor,M::AbstractMatrix{<:Number},n::Integer,t='b';variant='B')
  @assert(size(X1)==size(X2),"Dimensions mismatch")
    if sort(collect(size(vec(M))))[1]==1
        return mhadtv(X1,X2,vec(M),n,t,variant=variant);
  end
  sz=size(X1)
  In=sz[n]
  Im=prod(deleteat!(copy([sz...]),n))
  if t=='n'
    @assert(size(M,1) == Im, "Dimensions mismatch")
    Mprod=zeros(In,size(M,2))
  elseif t=='t'
    @assert(size(M,1) == In, "Dimensions mismatch")
    Mprod=zeros(Im,size(M,2))
  elseif t=='b'
    @assert(size(M,1) == In, "Dimensions mismatch")
    Mprod=zeros(In,size(M,2))
  end
  [Mprod[:,j]=  mhadtv(X1,X2,M[:,j],n,t,variant=variant) for j=1:size(M,2)]
  Mprod
end

function mhadtm(X1::ttensor,X2::ttensor,M::AbstractMatrix{<:Number},n::Integer,t='b';variant='B')
  @warn "Function mhadtm is depricated. Use mhadtv."
  mhadtm(X1,X2,M,n,t,variant=variant)
end

"""
    minus(X::ttensor,Y::ttensor)
    minus(X::ktensor,Y::ktensor)
    minus(X::htensor,Y::htensor)

Subtraction of two tensors. Same as: X-Y.
"""
function minus(X1::ttensor,X2::ttensor)
  1*X1+(-1)*X2
end
-(X1::ttensor,X2::ttensor) =minus(X1,X2)

#Multilinear rank of a ttensor. **Documentation in tensor.jl.
function mrank(X::ttensor)
  ntuple(n->nrank(X,n),ndims(X))
end
function mrank(X::ttensor,tol::Number)
   ntuple(n->nrank(X,n,tol),ndims(X))
end

"""
   msvdvals(X,n)

Singular values of mode-n matricization of a ttensor calculated directly. See also: nvecs.
"""
function msvdvals(X::ttensor,n::Integer)
  if X.isorth != true
    Y=reorth(X)
  end
  Gn=tenmat(Y.cten,n)
  svdvals(Gn)
end

"""
    mtimes(a,X)

Scalar times ttensor. Same as: a*X.
"""
function mtimes(α::Number,X::ttensor)
	ttensor(α*X.cten,X.fmat);
end
*(α::Number,X::ttensor)=mtimes(α,X)
*(X::ttensor,α::Number)=*(α,X)

#Matricized ttensor times Khatri-Rao product. **Documentation in tensor.jl.
function mttkrp(X::ttensor,M::MatrixCell,n::Integer)
  N=ndims(X)
  @assert(length(M) == N,"Wrong number of matrices.")
  modes=setdiff(1:N,n)
  sz=[size(X)...]
  K=size(M[modes[1]],2)
  @assert(!any(map(Bool,[size(M[m],2)-K for m in modes])),"Matrices must have the same number of columns")
  @assert(!any(map(Bool,[size(M[m],1)-sz[m] for m in modes])),"Matrices are of wrong size")
  fmat=MatrixCell(undef,N-1)
  i=1
  for m in modes
      fmat[i]=X.fmat[m]'*M[m]
      i+=1
  end
  Y=mttkrp(X.cten,fmat,n)
  X.fmat[n]*Y
end
mttkrp(X::ttensor,M::Array{Matrix{T}},n::Integer) where {T<:Number}=mttkrp(X,MatrixCell(M),n)

#Number of modes of a ttensor. **Documentation in Base.
function ndims(X::ttensor)
	ndims(X.cten)
end

#n-rank of a ttensor. **Documentation in tensor.jl.
function nrank(X::ttensor,n::Integer)
  rank(X.fmat[n])
end
function nrank(X::ttensor,n::Integer,tol::Number)
  rank(X.fmat[n],tol)
end

#Computes the r leading left singular vectors of mode-n matricization of a tensor X. **Documentation in tensor.jl.
function nvecs(X::ttensor,n::Integer,r=0;flipsign=false)
  if r==0
    r=size(X,n)
  end
  N=ndims(X)
  V=MatrixCell(undef,N)
  V[n]=X.fmat[n]
  for m in setdiff(1:N,n)
    V[m]=X.fmat[m]'*X.fmat[m]
  end
  H=ttm(X.cten,V)
  Hn=tenmat(H,n)
  Gn=tenmat(X.cten,n)
  #V=eigs(Symmetric(Hn*Gn'*X.fmat[n]'),nev=r,which=:LM)[2] #has bugs
  V=eigen(Symmetric(Hn*Gn'*X.fmat[n]')).vectors[:,end:-1:end-r+1]
  if flipsign
      maxind = findmax(abs.(V),1)[2]
      for i = 1:r
          ind=ind2sub(size(V),maxind[i])
          if V[ind...] < 0
             V[:,ind[2]] = V[:,ind[2]] * -1
          end
      end
  end
  V
end

#Permute dimensions of a ttensor. **Documentation in Base.
function permutedims(X::ttensor,perm::AbstractVector{<:Integer})
  @assert(collect(1:ndims(X))==sort(perm),"Invalid permutation")
  cten=permutedims(X.cten,perm)
  fmat=X.fmat[perm]
  ttensor(cten,fmat)
end

"""
    plus(X::ttensor,Y::ttensor)
    plus(X::ktensor,Y::ktensor)
    plus(X::htensor,Y::htensor)

Addition of two tensors. Same as: X+Y.
"""
function plus(X1::ttensor,X2::ttensor)
  @assert(size(X1) == size(X2),"Dimension mismatch.")
  fmat=[[X1.fmat[n] X2.fmat[n]] for n=1:ndims(X1)] #concatenate factor matrices
  cten=blockdiag(X1.cten,X2.cten)
  ttensor(cten,fmat)
end
+(X1::ttensor,X2::ttensor)=plus(X1,X2)

"""
    randrange(X,Y,n; <keyword arguments>)

Structure exploiting randomized range approximation of n-mode matricization of Hadamard product (X ∗ Y)ₙ, where X and Y are ttensors.

### Arguments:
- `reqrank::Integer`: Requested rank. Optional.
- `tol::Number/Vector`: Tolerance. Default: 1e-8.
- `maxit`: Maximal number of iterations. Default: 1000.
- `r`: Number of samples for stopping criterion. Default: r=10.
- `p::Integer`: Oversampling parameter. Defaul p=10.
"""
function randrange(X1::ttensor,X2::ttensor,mode::Integer;tol=1e-8,maxit=1000,reqrank=0,p=10,r=10)
  @assert(size(X1) == size(X2))
  N=ndims(X1)
  KR=MatrixCell(undef,N); #transpose Khatri-Rao product of X1.fmat and X2.fmat
  [KR[n]=khatrirao(X1.fmat[n],X2.fmat[n],'t') for n=1:N]
  randrange(X1.cten,X2.cten,KR,mode,tol=tol,maxit=maxit,reqrank=reqrank,p=p,r=r)
end

function randrange(C1::AbstractArray{<:Number,N},C2::AbstractArray{<:Number,N},KR::MatrixCell,mode::Integer;tol=1e-8,maxit=1000,reqrank=0,p=10,r=10) where N
  sz=zeros(Int,N)
  [sz[n]= size(KR[n],1) for n=1:N]
  m=sz[mode]
  n=prod(deleteat!(copy([sz...]),mode))
  remmodes=setdiff(1:N,mode);
  y=VectorCell(undef,N-1);
  if reqrank!=0
    Y=zeros(m,reqrank+p);
    for i=1:reqrank+p
      [y[k]=randn(size(KR[remmodes[k]],1)) for k=1:N-1]
      w=krontkron(reverse(KR[remmodes]),reverse(y),'t')
      Y[:,i]=KR[mode]*mkrontv(C1,C2,w,mode)
    end
    #Q=full(qrfact(Y)[:Q]);
    Q=Matrix(qr(Y).Q)
  else
    maxit=min(m,n,maxit);
    rangetol=tol*sqrt.(pi/2)/10;
    Y=zeros(m,r);
    for i=1:r
      [y[k]=randn(size(KR[remmodes[k]],1)) for k=1:N-1]
      w=krontkron(reverse(KR[remmodes]),reverse(y),'t')
      Y[:,i]=KR[mode]*mkrontv(C1,C2,w,mode);
    end
    j=0;
    Q=zeros(m,0);
    maxcolnorm=maximum([norm(Y[:,i]) for i=1:r])
    while maxcolnorm > rangetol && j<maxit
      j+=1;
      p=Q'*Y[:,j];
      Y[:,j]-=Q*p;
      q=Y[:,j]/norm(Y[:,j]);
      Q=[Q q];
      [y[k]=randn(size(KR[remmodes[k]],1)) for k=1:N-1]
      w=krontkron(reverse(KR[remmodes]),reverse(y),'t')
      w=KR[mode]*mkrontv(C1,C2,w,mode);
      p=Q'*w;
      Y=[Y w-Q*p]; #Y[:,j+r]=w-Q*p;
      p=q'*Y[:,j+1:j+r-1]
      Y[:,j+1:j+r-1]-=q*p;
      maxcolnorm=maximum([norm(Y[:,i]) for i=j+1:j+r])
    end
  end
  Q
end
randrange(C1::AbstractArray{<:Number,N},C2::AbstractArray{<:Number,N},KR::Array{Matrix{T}},mode::Integer;tol=1e-8,maxit=1000,reqrank=0,p=10,r=10) where {T<:Number,N}=randrange(C1,C2,MatrixCell(KR),mode,tol,maxit,reqrank,p,r)

"""
    randsvd(X,Y,n; <keyword arguments>)

Structure exploiting randomized SVD - computes left singular vectors and singular values of n-mode matricization (X ∗ Y)ₙ.
Works with matrix (X ∗ Y)ₙ(X ∗ Y)ₙᵀ.

## Arguments:
- `reqrank::Integer`: Requested rank. Optional.
- `variant` ∈ {'A','B'} Variant of multiplication (X ∗ Y)ₙ(X ∗ Y)ₙᵀ. Default: 'B'.
- `tol::Number/Vector`: Tolerance. Default: 1e-8.
- `maxit`: Maximal number of iterations. Default: 1000.
- `r`: Number of samples for stopping criterion. Default: r=10.
- `p::Integer`: Oversampling parameter. Defaul p=10.
"""
function randsvd(X1::ttensor,X2::ttensor,mode::Integer;variant='B',tol=1e-8,maxit=1000,reqrank=0,p=10,r=10)
  @assert(size(X1)==size(X2),"Dimensions mismatch")
  sz=size(X1)
  m=sz[mode]
  n=prod(deleteat!(copy([sz...]),mode))
  if reqrank!=0
     Y=mhadtv(X1,X2,randn(m,reqrank+p),mode,variant=variant);  #Y=A*(A'*randn(m,reqrank+p));
     #Q=full(qrfact(Y)[:Q]);
     Q=Matrix(qr(Y).Q)
  else
    maxit=min(m,n,maxit);
    rangetol=tol*sqrt.(pi/2)/10;
    Y=mhadtv(X1,X2,randn(m,r),mode,variant=variant);  #Y=A*(A'*randn(m,r));
    j=0;
    Q=zeros(m,0);
    maxcolnorm=maximum([norm(Y[:,i]) for i=1:r])
    while maxcolnorm > rangetol && j<maxit
      j+=1;
      p=Q'*Y[:,j];
      Y[:,j]-=Q*p;
      q=Y[:,j]/norm(Y[:,j]);
      Q=[Q q];
      w=mhadtv(X1,X2,randn(m),mode,variant=variant); #w=A*(A'*randn(m));
      p=Q'*w;
      Y=[Y w-Q*p]; #Y[:,j+r]=w-Q*p;
      p=q'*Y[:,j+1:j+r-1]
      Y[:,j+1:j+r-1]-=q*p;
      maxcolnorm=maximum([norm(Y[:,i]) for i=j+1:j+r])
    end
  end
  B=mhadtv(X1,X2,Q,mode,'t');#B=A'*Q;
  B=Symmetric(B'*B);
  #or (faster for small rank):
  #B=mhadtv(X1,X2,Q,mode,'n');
  #B=Symmetric(Q'*B);
  E=eigen(B,tol,Inf);
  U=E.vectors[:,end:-1:1];
  S=sqrt.(abs.(E.values[end:-1:1]));
  if reqrank != 0
    if reqrank > size(U,2)
      @warn "Required rank for mode $mode exceeds actual rank, the resulting rank will be smaller."
    else
      U=U[:,1:reqrank];
      S=S[1:reqrank];
    end
  end
  U=Q*U;
  U,S
end

"""
      reorth(X::ttensor)
      reorth(X::htensor)
      reorth(X::TTtensor[,direction,full])

Orthogonalize factor matrices of a tensor.
"""
function reorth(X::ttensor)
  N=ndims(X)
	if X.isorth
		X
	else
		Q=MatrixCell(undef,N)
		R=MatrixCell(undef,N)
    n=1;
		for A in X.fmat
			Qt,Rt=qr(A)
			Q[n]=Matrix(Qt)
			R[n]=Rt
      n+=1
		end
		ttensor(ttm(X.cten,R),Q)
	end
end

"""
    reorth!(X::ttensor)
    reorth!(X::htensor)
    reorth(X::TTtensor[,direction,full])

Orthogonalize factor matrices of a tensor. Rewrite ttensor.
"""
function reorth!(X::ttensor)
	if X.isorth != true
		for n=1:ndims(X)
			Q,R=qr(X.fmat[n])
			X.fmat[n]=Matrix(Q)
			X.cten=ttm(X.cten,R,n)
		end
    X.isorth=true;
  end
  X
end

#Size of a ttensor. **Documentation in Base.
function size(X::ttensor)
	tuple([size(X.fmat[n],1) for n=1:ndims(X)]...)
end
#Size of n-th mode of a ttensor.
function size(X::ttensor,n::Integer)
  size(X.fmat[n],1)
end

function Base.show(io::IO,X::ttensor)
    display(X)
end

#Mode-n matricization of a ttensor. **Documentation in tensor.jl.
tenmat(X::ttensor,n::Integer)=tenmat(full(X),n)

#ttensor times matrix (n-mode product). **Documentation in tensor.jl.
#t='t' transposes matrices
function ttm(X::ttensor,M::MatrixCell,modes::AbstractVector{<:Integer},t='n')
  if t=='t'
	  [M[n]=M[n]' for n=1:length(M)]
  end
  @assert(length(modes)<=length(M),"Too few matrices")
  @assert(length(M)<=ndims(X),"Too many matrices")
  if length(modes)<length(M)
    M=M[modes]
  end
  fmat=copy(X.fmat)
  for n=1:length(modes)
    fmat[modes[n]]=M[n]*X.fmat[modes[n]]
  end
    ttensor(X.cten,fmat)
end
ttm(X::ttensor,M::MatrixCell,modes::AbstractRange{<:Integer},t='n')=ttm(X,M,collect(modes),t)
ttm(X::ttensor,M::AbstractMatrix{<:Number},n::Integer,t='n') where{T<:Number}=ttm(X,MatrixCell([M]),[n],t)
ttm(X::ttensor,M::MatrixCell,t::Char)=ttm(X,M,1:length(M),t)
ttm(X::ttensor,M::MatrixCell)=ttm(X,M,1:length(M))
function ttm(X::ttensor,M::MatrixCell,n::Integer,t='n')
	if n>0
		ttm(X,M[n],n,t)
	else
		modes=setdiff(1:ndims(X),-n)
		ttm(X,M,modes,t)
	end
end
#If array of matrices isn't defined as MatrixCell, but as M=[M1,M2,...,Mn]:
ttm(X::ttensor,M::Array{Matrix{T}},modes::AbstractVector{<:Integer},t='n') where {T<:Number}=ttm(X,MatrixCell(M),modes,t)
ttm(X::ttensor,M::Array{Matrix{T}},modes::AbstractRange{<:Integer},t='n') where {T<:Number}=ttm(X,MatrixCell(M),modes,t)
ttm(X::ttensor,M::Array{Matrix{T}},t::Char) where {T<:Number}=ttm(X,MatrixCell(M),t)
ttm(X::ttensor,M::Array{Matrix{T}}) where {T<:Number}=ttm(X,MatrixCell(M))
ttm(X::ttensor,M::Array{Matrix{T}},n::Integer,t='n') where {T<:Number}=ttm(X,MatrixCell(M),n,t)

#ttensor times vector (n-mode product). **Documentation in tensor.jl.
#t='t' transposes matrices
function ttv(X::ttensor,V::VectorCell,modes::AbstractVector{<:Integer})
  N=ndims(X)
  remmodes=setdiff(1:N,modes)
  fmat=VectorCell(undef,N)
  if length(modes) < length(V)
    V=V[modes]
  end
  for n=1:length(modes)
    fmat[modes[n]]=X.fmat[modes[n]]'*V[n]
  end
  cten=ttv(X.cten,fmat,modes)
  if isempty(remmodes)
    cten
  else
    ttensor(cten,X.fmat[remmodes])
  end
end
ttv(X::ttensor,v::AbstractVector{<:Number},n::Integer)=ttv(X,VectorCell([v]),[n])
ttv(X::ttensor,V::VectorCell,modes::AbstractRange{<:Integer})=ttv(X,V,collect(modes))
ttv(X::ttensor,V::VectorCell)=ttv(X,V,1:length(V))
function ttv(X::ttensor,V::VectorCell,n::Integer)
	if n>0
		ttm(X,V[n],n)
	else
		modes=setdiff(1:ndims(X),-n)
		ttm(X,A,modes)
	end
end
#If array of vectors isn't defined as VectorCell, but as V=[v1,v2,...,vn]:
ttv(X::ttensor,V::Array{Vector{T}},modes::AbstractVector{<:Integer}) where {T<:Number}=ttv(X,VectorCell(V),modes)
ttv(X::ttensor,V::Array{Vector{T}},modes::AbstractRange{<:Integer}) where {T<:Number}=ttv(X,VectorCell(V),modes)
ttv(X::ttensor,V::Array{Vector{T}}) where {T<:Number}=ttv(X,VectorCell(V))
ttv(X::ttensor,V::Array{Vector{T}},n::Integer) where {T<:Number}=ttv(X,VectorCell(V),n)

"""
   uminus(X::ttensor)
   uminus(X::ktensor)
   uminus(X::htensor)

Unary minus. Same as: (-1)*X.
"""
uminus(X::ttensor)=mtimes(-1,X)
-(X::ttensor)=uminus(X)

#Frobenius norm of a ttensor. **Documentation in Base.
function norm(X::ttensor)
	if prod(size(X)) > prod(size(X.cten))
		if X.isorth
			norm(X.cten)
		else
			R=MatrixCell(undef,ndims(X))
			for n=1:ndims(X)
				#R[n]=qrfact(X.fmat[n])[:R]
				R[n]=qr(X.fmat[n]).R
			end
			norm(ttm(X.cten,R))
		end
	else
		norm(full(X))
	end
end


text:
module SfAtmosphere

using Unitful
using ..SfPlanetary, ..SfThermo, ..SfMatter

import PhysicalConstants.CODATA2018: k_B, G, R, g_n
import ..SfUnits: MolarMass

export scale_height, jeans_escape_timescale, jeans_parameter

scale_height(t, m::Unitful.Mass, g = g_n) = (k_B * t) / (m * g) |> u"km"

scale_height(t, m::MolarMass, g = g_n) = (R * t) / (m * g) |> u"km"

"""
    jeans_escape_timescale(T::Unitful.Temperature, exobase::Unitful.Length, M::Unitful.Mass, R::Unitful.Length, m::Unitful.Mass)

Approximate Jeans escape timescale for a gas of atomic mass `m` on a planet of mass `M` and radius `R`, with an exobase at altitude `xh` and temperature `T`.

Working from http://cococubed.asu.edu/code_pages/jeans_escape.shtm
"""
function jeans_escape_timescale(T::Unitful.Temperature, xh::Unitful.Length, M::Unitful.Mass, R::Unitful.Length, m::Unitful.Mass)
   g = (G*M)/(R+xh)^2 # gravity at exobase
   H = (k_B*T)/(m*g) # scale height for gas
   v_peak=maxwell_boltzmann_peak_speed(m, T) # peak of maxwell-boltzmann distribution
   v_esc=sqrt((2G*M)/(R+xh)) # escape velocity of planet at exosphere altitude
   λ =(v_esc/v_peak)^2
   v_jeans = v_peak * ((1 + λ)*exp(-λ))/sqrt(4π)
   return H/v_jeans |> u"yr"
end

"""
    jeans_escape_timescale(T::Unitful.Temperature, xh::Unitful.Length, body::AbstractBody, m::Unitful.Mass)

Approximate Jeans escape timescale for a gas of atomic mass `m` on planet `body`, with an exobase at altitude `xh` and temperature `T`.

Working from http://cococubed.asu.edu/code_pages/jeans_escape.shtm
"""
jeans_escape_timescale(T::Unitful.Temperature, xh::Unitful.Length, body::AbstractBody, m::Unitful.Mass) =
	jeans_escape_timescale(T, xh, mass(body), radius(body), m)

# https://arxiv.org/ftp/arxiv/papers/1009/1009.5110.pdf
# https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2008GL036513
jeans_parameter(m_planet::Unitful.Mass, m_molecule::Unitful.Mass, t_exosphere::Unitful.Temperature, r_exosphere::Unitful.Length) =
	(G * m_planet * m_molecule) / (k_B * t_exosphere * r_exosphere) |> u"m/m"

jeans_parameter(body::AbstractBody, m_molecule::Unitful.Mass, t_exosphere::Unitful.Temperature, r_exosphere::Unitful.Length) =
	jeans_parameter(mass(body), m_molecule, t_exosphere, r_exosphere)

end

max_stars_repo_name: icantbelieveitsnotphysics/SfPhysics

max_issues_repo_name: icantbelieveitsnotphysics/SfPhysics

max_forks_repo_name: icantbelieveitsnotphysics/SfPhysics


text:
module Gomagorithm

include("SOneTo.jl")

end # module


text:
using LazyArtifacts

struct DictionaryEntry
    trad::String
    simp::String
    pinyin::String
    senses::Vector{String}
end

function Base.print(io::IO, entry::DictionaryEntry)
    char_str = entry.trad == entry.simp ? entry.trad : "$(entry.trad) ($(entry.simp))"
    print(io, "$char_str: [$(entry.pinyin)]\n")
    print(io, join(map(w -> "\t" * w, entry.senses), "\n"))
    return nothing
end

struct ChineseDictionary
    entries::Dict{String, Vector{DictionaryEntry}}

    """
        ChineseDictionary([filename])

    Load a text-based dictionary either from the default dictionary file or from the provided
    filename. The format of the text file must be the same as
    [that used by the CC-CEDICT project](https://cc-cedict.org/wiki/format:syntax) for
    compatibility reasons.

    For general use, it's the easiest to just use the default dictionary (from the CC-CEDICT project).
    This is loaded if you don't specify a filename.
    """
    function ChineseDictionary(filename=joinpath(artifact"cedict", "cedict_ts.u8"))
        dict = Dict{String, Vector{DictionaryEntry}}()
        pattern = r"^([^#\s]+) ([^\s]+) \[(.*)\] /(.+)/$"

        for line in eachline(filename)
            m = match(pattern, line)
            m === nothing && continue

            trad, simp, pinyin, defns = String.(m.captures)
            entry = DictionaryEntry(trad, simp, pinyin, split(defns, "/"))

            dict[trad] = push!(get(dict, trad, []), entry)
            simp != trad && (dict[simp] = push!(get(dict, simp, []), entry))
        end

        new(dict)
    end
end

# iteration
Base.iterate(dict::ChineseDictionary) = iterate(dict.entries)
Base.IteratorSize(::Type{ChineseDictionary}) = HasLength()
Base.IteratorEltype(::Type{ChineseDictionary}) = HasEltype()
Base.length(dict::ChineseDictionary) = length(dict.entries)
Base.eltype(::Type{ChineseDictionary}) = Pair{String, Vector{DictionaryEntry}}

# indexing
Base.getindex(dict::ChineseDictionary, i) = getindex(dict.entries, i)
Base.setindex!(dict::ChineseDictionary, v, i) = setindex!(dict.entries, v, i)

# dictionaries
Base.keys(dict::ChineseDictionary) = keys(dict.entries)
Base.values(dict::ChineseDictionary) = values(dict.entries)
Base.haskey(dict::ChineseDictionary, key) = haskey(dict.entries, key)


text:
using ITensorBenchmarks

@testset "runbenchmarks" begin
  runbenchmarks(stdout; benchmarks = ["trg", "ctmrg", "dmrg_1d"], test = 1:2, blas_num_threads = 1:2)
end