text:
#######################################################################
# Contains code from  Catlab.jl:
# https://raw.githubusercontent.com/AlgebraicJulia/Catlab.jl/master/src/graphics/Graphviz.jl
#
# That license for that code is:
#
# The MIT License
#
# Copyright (c) 2017-2020: Evan Patterson.
#
# 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.
#######################################################################
""" AST and pretty printer for Graphviz's DOT language.

References:

- DOT grammar: http://www.graphviz.org/doc/info/lang.html
- DOT language guide: http://www.graphviz.org/pdf/dotguide.pdf
"""

# AST
#####

abstract type Expression end
abstract type Statement <: Expression end

""" AST type for Graphviz's "HTML-like" node labels.

For now, the HTML content is just a string.
"""
struct Html
  content::String
end
Base.print(io::IO, html::Html) = print(io, html.content)

const AttributeValue = Union{String,Html}
const Attributes = OrderedDict{Symbol,AttributeValue}

as_attributes(attrs::Attributes) = attrs
as_attributes(d::OrderedDict) = Attributes(Symbol(k) => d[k] for k in keys(d))
as_attributes(d::AbstractDict) =
  Attributes(Symbol(k) => d[k] for k in sort!(collect(keys(d))))

@with_kw_noshow struct Graph <: Expression
  name::String
  directed::Bool
  prog::String = "dot"
  stmts::Vector{Statement} = Statement[]
  graph_attrs::Attributes = Attributes()
  node_attrs::Attributes = Attributes()
  edge_attrs::Attributes = Attributes()
end

Graph(name::String, stmts::Vector{Statement}; kw...) =
  Graph(; name=name, directed=false, stmts=stmts, kw...)
Graph(name::String, stmts::Vararg{Statement}; kw...) =
  Graph(; name=name, directed=false, stmts=collect(stmts), kw...)
Digraph(name::String, stmts::Vector{Statement}; kw...) =
  Graph(; name=name, directed=true, stmts=stmts, kw...)
Digraph(name::String, stmts::Vararg{Statement}; kw...) =
  Graph(; name=name, directed=true, stmts=collect(stmts), kw...)

@with_kw_noshow struct Subgraph <: Statement
  name::String = "" # Subgraphs can be anonymous
  stmts::Vector{Statement} = Statement[]
  graph_attrs::Attributes = Attributes()
  node_attrs::Attributes = Attributes()
  edge_attrs::Attributes = Attributes()
end

Subgraph(stmts::Vector{Statement}; kw...) = Subgraph(; stmts=stmts, kw...)
Subgraph(stmts::Vararg{Statement}; kw...) = Subgraph(; stmts=collect(stmts), kw...)
Subgraph(name::String, stmts::Vector{Statement}; kw...) =
  Subgraph(; name=name, stmts=stmts, kw...)
Subgraph(name::String, stmts::Vararg{Statement}; kw...) =
  Subgraph(; name=name, stmts=collect(stmts), kw...)

struct Node <: Statement
  name::String
  attrs::Attributes
end
Node(name::String, attrs::AbstractDict) = Node(name, as_attributes(attrs))
Node(name::String; attrs...) = Node(name, attrs)

struct NodeID <: Expression
  name::String
  port::String
  anchor::String
  NodeID(name::String, port::String="", anchor::String="") = new(name, port, anchor)
end

struct Edge <: Statement
  path::Vector{NodeID}
  attrs::Attributes
end
Edge(path::Vector{NodeID}, attrs::AbstractDict) = Edge(path, as_attributes(attrs))
Edge(path::Vector{NodeID}; attrs...) = Edge(path, attrs)
Edge(path::Vararg{NodeID}; attrs...) = Edge(collect(path), attrs)
Edge(path::Vector{String}, attrs::AbstractDict) = Edge(map(NodeID, path), attrs)
Edge(path::Vector{String}; attrs...) = Edge(map(NodeID, path), attrs)
Edge(path::Vararg{String}; attrs...) = Edge(map(NodeID, collect(path)), attrs)

# Bindings
##########

""" Run a Graphviz program.

Invokes Graphviz through its command-line interface. If the `Graphviz_jll`
package is installed and loaded, it is used; otherwise, Graphviz must be
installed on the local system.

For bindings to the Graphviz C API, see the the package
[GraphViz.jl](https://github.com/Keno/GraphViz.jl). At the time of this writing,
GraphViz.jl is unmaintained.
"""
function run_graphviz(io::IO, graph::Graph; prog::Union{String,Nothing}=nothing,
  format::String="json0")
  if isnothing(prog)
    prog = graph.prog
  end
  @assert prog in ("dot", "neato", "fdp", "sfdp", "twopi", "circo")
  if USE_GV_JLL[]
    fun = getfield(Graphviz_jll, Symbol(prog))
    fun() do exe
      open(`$exe -T$format`, io, write=true) do gv
        pprint(gv, graph)
      end
    end
  else
    open(`$prog -T$format`, io, write=true) do gv
      pprint(gv, graph)
    end
  end
end
function run_graphviz(graph::Graph; kw...)
  io = IOBuffer()
  run_graphviz(io, graph; kw...)
  seekstart(io)
end

function Base.show(io::IO, ::MIME"image/svg+xml", graph::Graph)
  run_graphviz(io, graph, format="svg")
end

# Pretty-print
##############

""" Pretty-print the Graphviz expression.
"""
pprint(expr::Expression) = pprint(stdout, expr)
pprint(io::IO, expr::Expression) = pprint(io, expr, 0)

function pprint(io::IO, graph::Graph, n::Int)
  indent(io, n)
  print(io, graph.directed ? "digraph " : "graph ")
  print(io, graph.name)
  println(io, " {")
  pprint_attrs(io, graph.graph_attrs, n + 2; pre="graph", post=";\n")
  pprint_attrs(io, graph.node_attrs, n + 2; pre="node", post=";\n")
  pprint_attrs(io, graph.edge_attrs, n + 2; pre="edge", post=";\n")
  for stmt in graph.stmts
    pprint(io, stmt, n + 2, directed=graph.directed)
    println(io)
  end
  indent(io, n)
  println(io, "}")
end

function pprint(io::IO, subgraph::Subgraph, n::Int; directed::Bool=false)
  indent(io, n)
  if isempty(subgraph.name)
    println(io, "{")
        else
    print(io, "subgraph ")
    print(io, subgraph.name)
    println(io, " {")
  end
  pprint_attrs(io, subgraph.graph_attrs, n + 2; pre="graph", post=";\n")
  pprint_attrs(io, subgraph.node_attrs, n + 2; pre="node", post=";\n")
  pprint_attrs(io, subgraph.edge_attrs, n + 2; pre="edge", post=";\n")
  for stmt in subgraph.stmts
    pprint(io, stmt, n + 2, directed=directed)
    println(io)
  end
  indent(io, n)
  print(io, "}")
end

function pprint(io::IO, node::Node, n::Int; directed::Bool=false)
  indent(io, n)
  print(io, node.name)
  pprint_attrs(io, node.attrs)
  print(io, ";")
end

function pprint(io::IO, node::NodeID, n::Int)
  print(io, node.name)
  if !isempty(node.port)
    print(io, ":")
    print(io, node.port)
  end
  if !isempty(node.anchor)
    print(io, ":")
    print(io, node.anchor)
end
end

function pprint(io::IO, edge::Edge, n::Int; directed::Bool=false)
  indent(io, n)
  for (i, node) in enumerate(edge.path)
    if i > 1
      print(io, directed ? " -> " : " -- ")
    end
    pprint(io, node, n)
  end
  pprint_attrs(io, edge.attrs)
  print(io, ";")
end

function pprint_attrs(io::IO, attrs::Attributes, n::Int=0;
                      pre::String="", post::String="")
  if !isempty(attrs)
    indent(io, n)
        print(io, pre)
    print(io, " [")
    for (i, (key, value)) in enumerate(attrs)
      if (i > 1) print(io, ",") end
      print(io, key)
      print(io, "=")
      print(io, value isa Html ? "<" : "\"")
      print(io, value)
      print(io, value isa Html ? ">" : "\"")
    end
    print(io, "]")
    print(io, post)
  end
end

indent(io::IO, n::Int) = print(io, " "^n)

#######################################################################
# following is adapted from Petri.jl
# https://github.com/mehalter/Petri.jl
#######################################################################

const graph_attrs = Attributes(:rankdir => "LR")
const node_attrs  = Attributes(:shape => "plain", :style => "filled", :color => "white")
const edge_attrs  = Attributes(:splines => "splines")

function edgify(δ, i, reverse::Bool)
    attr = Attributes()
    return map(δ) do p
        val = String(p[1].f.name)
      weight = "$(p[2])"
      attr = Attributes(:label => weight, :labelfontsize => "6")
      return Edge(reverse ? ["rx_$i", "$val"] :
                            ["$val", "rx_$i"], attr)
    end
end

# make distinguished edge based on rate constant
function edgifyrates(rxs, specs)
    es = Edge[]
    for (i, rx) in enumerate(rxs)
        deps = rx.rate isa Number ? Any[] : get_variables(rx.rate, specs)
        for dep in deps
            val = String(dep.f.name)
            attr = Attributes(:color => "#d91111", :style => "dashed")
            e = Edge(["$val", "rx_$i"], attr)
            push!(es, e)
        end
    end
    es
end

# create edges from one reaction complex to another reaction complex
function edgifycomplex(δ,attr)
    return map(δ) do p
        return Edge([p[1], p[2]] , attr)
    end
end

# modify vector of string of complexes into graphviz compatible strings
function modifystrcomp(strcomp::Vector{String})
  for i in 1:length(strcomp)
    #   cannot allow (t) as per graphviz
      if occursin("(t)",strcomp[i])
          strcomp[i] = replace(strcomp[i], "(t)" => "")
      end
      #   prettify NUll with ∅
      if occursin("0",strcomp[i])
          strcomp[i] = replace(strcomp[i], "0" => "∅")
      end
  end
  strcomp = "<".*strcomp.*">"
end

"""
    complexgraph(rn::ReactionSystem; complexdata=reactioncomplexes(rn))

Creates a Graphviz graph of the [`ReactionComplex`](@ref)s in `rn`. Reactions
correspond to arrows and reaction complexes to blue circles.

Notes:
- Black arrows from complexes to complexes indicate reactions whose rate is a
  parameter or a `Number`. i.e. `k, A --> B`.
- Red dashed arrows from complexes to complexes indicate reactions whose rate
  depends on species. i.e. `k*C, A --> B` for `C` a species.
- Requires the Graphviz jll to be installed, or Graphviz to be installed and
  commandline accessible.
"""
function complexgraph(rn::ReactionSystem; complexdata=reactioncomplexes(rn))
    rxs   = reactions(rn);
    specs = species(rn);
    complexes,B = complexdata;
    fun = rcel -> specs[rcel.speciesid]*rcel.speciesstoich;
    compfun(rc) = rc == Catalyst.ReactionComplex{Int64}[] ? 0 : sum(fun, rc);

    strcomp = [string(compfun(rc)) for rc in complexes];

    newstrcomp = modifystrcomp(strcomp)
    compnodes = [Node(str, Attributes(:shape => "circle",:color => "#6C9AC3")) for str in newstrcomp]

    edges = []
    for (i,r) in enumerate(rxs)
        subcomp = newstrcomp[argmin(@view B[:,i])]
        prodcomp = newstrcomp[argmax(@view B[:,i])]
        deps = get_variables(r.rate, specs)
        if deps != Any[]
            attr = Attributes(:color => "#d91111", :style => "dashed")
            push!(edges, edgifycomplex(zip([subcomp],[prodcomp]),attr))
        else
            attr = Attributes()
            push!(edges,edgifycomplex(zip([subcomp],[prodcomp]),attr))
        end
    end
    stmts2 = Vector{Statement}()
    append!(stmts2, compnodes)
    append!(stmts2, collect(Iterators.flatten(edges)))
    g = Digraph("G", stmts2; graph_attrs=graph_attrs, node_attrs=node_attrs,edge_attrs=edge_attrs)
    return g
end

"""
    Graph(rn::ReactionSystem)

Converts a [`ReactionSystem`](@ref) into a Graphviz graph.
Reactions correspond to small green circles, and species to blue circles.

Notes:
- Black arrows from species to reactions indicate reactants, and are labelled
  with their input stoichiometry.
- Black arrows from reactions to species indicate products, and are labelled
  with their output stoichiometry.
- Red arrows from species to reactions indicate that species is used within the
  rate expression. For example, in the reaction `k*A, B --> C`, there would be a
  red arrow from `A` to the reaction node. In `k*A, A+B --> C`, there would be
  red and black arrows from `A` to the reaction node.
- Requires the Graphviz jll to be installed, or Graphviz to be installed and
  commandline accessible.
"""
function Graph(rn::ReactionSystem)
    rxs   = reactions(rn)
    specs = species(rn)
    statenodes = [Node(string(s.f.name), Attributes(:shape => "circle", :color => "#6C9AC3")) for s in specs]
    transnodes = [Node(string("rx_$i"), Attributes(:shape => "point", :color => "#E28F41", :width => ".1")) for (i, r) in enumerate(rxs)]

    stmts = vcat(statenodes, transnodes)
    edges = map(enumerate(rxs)) do (i, r)
      vcat(edgify(zip(r.substrates, r.substoich), i, false),
           edgify(zip(r.products, r.prodstoich), i, true))
    end
    es = edgifyrates(rxs, specs)
    (!isempty(es)) && push!(edges, es)

    stmts2 = Vector{Statement}()
    append!(stmts2, stmts)
    append!(stmts2, collect(Iterators.flatten(edges)))
    g = Digraph("G", stmts2; graph_attrs=graph_attrs, node_attrs=node_attrs, edge_attrs=edge_attrs)
    return g
end

"""
    savegraph(g::Graph, fname, fmt="png")

Given a `Graph` generated by [`Graph`](@ref), save the graph to the file with
name `fname` and extension `fmt`.

Notes:
- `fmt="png"` is the default output format.
- Requires the Graphviz jll to be installed, or Graphviz to be installed and commandline accessible.
"""
function savegraph(g::Graph, fname, fmt="png")
    open(fname, "w") do io
        run_graphviz(io, g, format=fmt)
    end
    nothing
end


text:
for func in [:clblasSasum, :clblasDasum, :clblasScasum, :clblasDzsum]
    @eval @blasfun $func(N::Csize_t, asum::cl.CL_mem, offAsum::Csize_t,
                              X::cl.CL_mem, offx::Csize_t, incx::Cint,
                              scratch_buff::cl.CL_mem,
                              n_queues::cl.CL_uint,
                              queues::Ptr{cl.CL_command_queue},
                              n_events_in_wait_list::cl.CL_uint,
                              event_wait_list::Ptr{cl.CL_event},
                              events::Ptr{cl.CL_event})

    @eval @blasfun2 $func(N::Csize_t, asum::cl.CL_mem, offAsum::Csize_t,
                               X::cl.CL_mem, offx::Csize_t, incx::Cint,
                               scratch_buff::cl.CL_mem)
end

ext: jl


text:
export batched_mul, ⊠,  batched_vec
export batched_transpose, batched_adjoint

include("./batchedadjtrans.jl")

using LinearAlgebra: BlasFloat, Transpose, Adjoint, AdjOrTransAbsMat

_unbatch(A) = A
_unbatch(A::BatchedAdjOrTrans) = parent(A)

"""
    batched_mul(A, B) -> C
    A ⊠ B  # \\boxtimes

Batched matrix multiplication. Result has `C[:,:,k] == A[:,:,k] * B[:,:,k]` for all `k`.
If `size(B,3) == 1` then instead `C[:,:,k] == A[:,:,k] * B[:,:,1]`, and similarly for `A`.

To transpose each matrix, apply `batched_transpose` to the array,
or `batched_adjoint` for conjugate-transpose:

```jldoctest
julia> A, B = randn(2,5,17), randn(5,9,17);

julia> A ⊠ B |> size
(2, 9, 17)

julia> batched_adjoint(A) |> size
(5, 2, 17)

julia> batched_mul(A, batched_adjoint(randn(9,5,17))) |> size
(2, 9, 17)

julia> A ⊠ randn(5,9,1) |> size
(2, 9, 17)

julia> batched_transpose(A) == PermutedDimsArray(A, (2,1,3))
true
```

The equivalent `PermutedDimsArray` may be used in place of `batched_transpose`.
Other permutations are also handled by BLAS,
provided that the batch index `k` is not the first dimension of the underlying array.
Thus `PermutedDimsArray(::Array, (1,3,2))` and `PermutedDimsArray(::Array, (3,1,2))` are fine.

However, `A = PermutedDimsArray(::Array, (3,2,1))` is not acceptable to BLAS,
since the batch dimension is the contiguous one: `stride(A,3) == 1`.
This will be copied, as doing so is faster than `batched_mul_generic!`.

Both this `copy` and `batched_mul_generic!` produce `@debug` messages,
and setting for instance `ENV["JULIA_DEBUG"] = NNlib` will display them.
"""
function batched_mul(A::AbstractArray{T1, 3}, B::AbstractArray{T2, 3}) where {T1, T2}
    size(A, 3) == size(B, 3) || size(A, 3) == 1 || size(B, 3) == 1 ||
        throw(DimensionMismatch("batch size mismatch: A != B"))
    _batched_mul(storage_typejoin(A, B), A, B)
end

const ⊠ = batched_mul

function _batched_mul(::Type, A, B)
    T = promote_type(eltype(A), eltype(B))
    C = similar(A, T, (size(A, 1), size(B, 2), max(size(A, 3), size(B, 3))))
    batched_mul!(C, A, B)
    C
end
function _batched_mul(::Type{<:DenseArray{T}}, A, B) where {T<:BlasFloat}
    C = similar(A, T, (size(A, 1), size(B, 2), max(size(A, 3), size(B, 3))))
    batched_mul!(C, _copy_if_faster(A), _copy_if_faster(B))
    C
end

function _copy_if_faster(X::AbstractArray{<:Number, 3})
    is_strided(X) || return X
    if Base.stride(X, 3) == 1 && Base.stride(X, 1) != 1
        @debug "copying to avoid batched_mul_generic!" typeof(X) size(X) strides(X)
        return copy(X)
    end
    X
end
function _copy_if_faster(X::BatchedAdjoint{<:Complex})
    Xbase = _unbatch(X)
    is_strided(Xbase) || return X
    if Base.stride(Xbase, 1) != 1
        @debug "copying to avoid batched_mul_generic!" typeof(X) size(X) strides(_unbatch(X))
        return copy(X) # or batched_adjoint(copy(Xbase)), may be better on GPU?
    end
    X
end

# Gradient, allowing that size(A,3)==1 means it's "broadcasted" out to size(B,3)

function rrule(::typeof(batched_mul), A::AbstractArray{<:Any,3}, B::AbstractArray{<:Any,3})
    function batched_mul_pullback(Δ)
        Athunk = @thunk begin
            tmp = batched_mul(Δ, batched_adjoint(B))
            size(A,3) == 1 ? sum(tmp, dims=3) : tmp
        end
        Bthunk = @thunk begin
            tmp = batched_mul(batched_adjoint(A), Δ)
            size(B,3) == 1 ? sum(tmp, dims=3) : tmp
        end
        return (NoTangent(), Athunk, Bthunk)
    end
    batched_mul(A, B), batched_mul_pullback
end

"""
    batched_mul(A::Array{T,3}, B::Matrix)
    batched_mul(A::Matrix, B::Array{T,3})
    A ⊠ B

This is always matrix-matrix multiplication, but
either `A` or `B` may lack a batch index.

* When `B` is a matrix, result has `C[:,:,k] == A[:,:,k] * B[:,:]` for all `k`.

* When `A` is a matrix, then `C[:,:,k] == A[:,:] * B[:,:,k]`.
  This can also be done by reshaping and calling `*`,
  for instance `A ⊡ B` using TensorCore.jl, but is implemented here using
  `batched_gemm` instead of `gemm`.

```jldoctest
julia> randn(16,8,32) ⊠ randn(8,4) |> size
(16, 4, 32)

julia> randn(16,8,32) ⊠ randn(8,4,1) |> size  # equivalent
(16, 4, 32)

julia> randn(16,8) ⊠ randn(8,4,32) |> size
(16, 4, 32)
```

See also `batched_vec` to regard `B` as a batch of vectors, `A[:,:,k] * B[:,k]`.
"""
batched_mul(A::AbstractArray{T,3} where T, B::AbstractMatrix) = _semi_batched_mul(A,B)

# Simplify signature of batched_mul by hiding dispatch on Adjoint etc:

_semi_batched_mul(A::AbstractArray{<:Any,3}, B::AbstractMatrix) =
    batched_mul(A, reshape(B, size(B)..., 1))

_semi_batched_mul(A::AbstractArray{<:Any,3}, B::Adjoint{<:Number,<:AbstractMatrix}) =
    batched_mul(A, batched_adjoint(reshape(parent(B), size(parent(B))..., 1)))

_semi_batched_mul(A::AbstractArray{<:Any,3}, B::Transpose{<:Number,<:AbstractMatrix}) =
    batched_mul(A, batched_transpose(reshape(parent(B), size(parent(B))..., 1)))

batched_mul(A::AbstractMatrix, B::AbstractArray{T,3} where T) = _semi_batched_mul(A,B)

_semi_batched_mul(A::AbstractMatrix, B::AbstractArray{<:Any,3}) =
    batched_mul(reshape(A, size(A)..., 1), B)

_semi_batched_mul(A::Adjoint{<:Number,<:AbstractMatrix}, B::AbstractArray{<:Any,3}) =
    batched_mul(batched_adjoint(reshape(parent(A), size(parent(A))..., 1)), B)

_semi_batched_mul(A::Transpose{<:Number,<:AbstractMatrix}, B::AbstractArray{<:Any,3}) =
    batched_mul(batched_transpose(reshape(parent(A), size(parent(A))..., 1)), B)

"""
    batched_vec(A::Array{T,3}, B::Matrix)
    batched_vec(A::Array{T,3}, b::Vector)

Batched matrix-vector multiplication:
the result has `C[:,:,k] == A[:,:,k] * B[:,k]` for all `k`,
or else `C[:,:,k] == A[:,:,k] * b` for `b::Vector`.

With the same argument types, `batched_mul(A, B)` would regard `B` as
a fixed matrix, not a batch of vectors. Both reshape and then
call `batched_mul(::Array{T,3}, ::Array{T,3})`.

```jldoctest
julia> A, B, b = randn(16,8,32), randn(8,32), randn(8);

julia> batched_vec(A,B) |> size
(16, 32)

julia> batched_vec(A,b) |> size
(16, 32)
```
"""
function batched_vec(A::AbstractArray{T,3} where T, B::AbstractMatrix)
    # If B is transposed, then stride=1 is the batch dim, so we will end up copying anyway:
    if B isa AdjOrTransAbsMat{<:BlasFloat, <:StridedMatrix}
        return batched_vec(A, copy(B))
    end
    reshape(batched_mul(A, reshape(B, size(B,1), 1, size(B,2))), size(A,1), size(A,3))
end

batched_vec(A::AbstractArray{T,3} where T, b::AbstractVector) =
    reshape(batched_mul(A, reshape(b, length(b), 1, 1)), size(A,1), size(A,3))

"""
    batched_mul!(C, A, B) -> C
    batched_mul!(C, A, B, α=1, β=0)

In-place batched matrix multiplication, equivalent to
`mul!(C[:,:,k], A[:,:,k], B[:,:,k], α, β)` for all `k`.
If `size(B,3) == 1` then every batch uses `B[:,:,1]` instead.

This will call `batched_gemm!` whenever possible. For real arrays this means that,
for `X ∈ [A,B,C]`, either `strides(X,1)==1` or `strides(X,2)==1`, the latter may
be caused by `batched_transpose` or by for instance `PermutedDimsArray(::Array, (3,1,2))`.
Unlike `batched_mul` this will never make a copy.

For complex arrays, the wrapper made by `batched_adjoint` must be outermost to be seen.
In this case the strided accepted by BLAS are more restricted, if `stride(C,1)==1` then
only `stride(AorB::BatchedAdjoint,2) == 1` is accepted.
"""
function batched_mul!(C::AbstractArray{T,3}, A::AbstractArray{<:Any,3}, B::AbstractArray{<:Any,3},
        α::Number=one(T), β::Number=zero(T)) where {T}
    _batched_mul!(storage_typejoin(C,A,B), C, A, B, α, β)
    C
end

_batched_mul!(::Type, C, A, B, α::Number, β::Number) = batched_mul_generic!(C, A, B, α, β)

_batched_mul!(::Type{DT}, C, A, B, α::Number, β::Number) where {DT<:DenseArray{T}} where {T<:BlasFloat} =
    _batched_try_gemm!(DT, C, A, B, α, β)

function _batched_try_gemm!(::Type{DT}, C, A, B, α::Number, β::Number) where {DT<:DenseArray{T}} where {T<:BlasFloat}

    alpha, beta = promote(α, β, zero(T))
    alpha isa T && beta isa T || return batched_mul_generic!(C, A, B, α, β)

    are_strided(_unbatch(A), _unbatch(B)) || return batched_mul_generic!(C, A, B, α, β)
    C isa StridedArray || return batched_mul_generic!(C, A, B, α, β)

    blasA, transA = if A isa BatchedAdjoint && T <: Complex
        Base.stride(parent(A),1) == 1 || return batched_mul_generic!(C, A, B, α, β)
        parent(A), 'C'
    elseif Base.stride(A,2) == 1 && size(A,1) > 1
        batched_transpose(A), 'T'
    elseif Base.stride(A,1) == 1
        A, 'N'
    elseif Base.stride(A,2) == 1  # This is awful, but exhaustively tested. Issues 268, 282.
        batched_transpose(A), 'T'
    else
        return batched_mul_generic!(C, A, B, α, β)
    end

    blasB, transB = if B isa BatchedAdjoint && T <: Complex
        Base.stride(parent(B),1) == 1 || return batched_mul_generic!(C, A, B, α, β)
        parent(B), 'C'
    elseif Base.stride(B,2) == 1 && size(B,1) > 1
        batched_transpose(B), 'T'
    elseif Base.stride(B,1) == 1
        B, 'N'
    elseif Base.stride(B,2) == 1
        batched_transpose(B), 'T'
    else
        return batched_mul_generic!(C, A, B, α, β)
    end

    _batched_gemm!(DT, transA, transB, alpha, blasA, blasB, beta, C)
    C
end

_batched_gemm!(::Type{<:Array}, transA::Char, transB::Char, α::Number, A, B, β::Number, C) =
    batched_gemm!(transA, transB, α, A, B, β, C)

_BATCHED_LIST = [
    (:(AbstractArray{<:Any, 3}), :identity),
    (:BatchedTranspose,          :transpose),
    (:BatchedAdjoint,            :adjoint),
]
for (TA, fA) in _BATCHED_LIST, (TB, fB) in _BATCHED_LIST

    @eval function batched_mul_generic!(C::AbstractArray{T, 3}, A::$TA, B::$TB,
            α::Number=one(T), β::Number=zero(T)) where {T}

        size(A, 3) == size(C, 3) || size(A, 3) == 1 || throw(DimensionMismatch("batch size mismatch: A != C"))
        size(B, 3) == size(C, 3) || size(B, 3) == 1 || throw(DimensionMismatch("batch size mismatch: B != C"))
        @debug "calling fallback method for batched_mul!" typeof(A) size(A) typeof(B) size(B) typeof(C)

        Abase, Bbase = _unbatch(A), _unbatch(B)
        sA, oA = size(A,3) == 1 ? (0,1) : (1,0)
        sB, oB = size(B,3) == 1 ? (0,1) : (1,0)

        @inbounds for k in 1:size(C,3)
            @views mul!(C[:,:,k], $fA(Abase[:,:,k*sA+oA]), $fB(Bbase[:,:,k*sB+oB]), α, β)
        end
        C
    end

end

"""
    storage_type(A) -> Type

Removes all wrappers to return the `Array` or `CuArray` (or whatever) type within.
```
julia> view(reshape(ones(10)',2,5),:, 3:4) |> storage_type
Array{Float64,1}

julia> reshape(sparse(rand(10)), 5,2) |> storage_type
SparseVector{Float64,Int64}
```
"""
function storage_type(A::AbstractArray)
    P = parent(A)
    typeof(A) === typeof(P) ? typeof(A) : storage_type(P)
end
storage_type(A) = typeof(A)

"""
    storage_typejoin(A, B, C, ...) -> Type

Reduces with `Base.promote_typejoin`, in order that this conveys useful information
for dispatching to BLAS. It does not tell you what container to allocate:
```
julia> storage_typejoin(rand(2), rand(Float32, 2))
Array{T,1} where T

julia> eltype(ans) <: LinearAlgebra.BlasFloat
false

julia> storage_typejoin(rand(2), rand(2,3), rand(2,3,4))
Array{Float64,N} where N
```
"""
storage_typejoin(A, Bs...) = Base.promote_typejoin(storage_type(A), storage_typejoin(Bs...))
storage_typejoin(A) = storage_type(A)

"""
    is_strided(A::AbstractArray) -> Bool

This generalises `A isa StridedArray` to treat wrappers like `A::PermutedDimsArray`,
for which it returns `is_strided(parent(A))`.

It returns `true` for `CuArray`s, and `PermutedDimsArray`s of those.

Other wrappers (defined outside Base, LinearAlgebra) are assumed not to break
strided-ness, and hence also return `is_strided(parent(A))`.
This correctly handles things like `NamedDimsArray` wihch don't alter indexing.
However, it's a little pessimistic in that e.g. a `view` of such a container will return
`false`, even in cases where the same `view` of `parent(A)` would be a `StridedArray`.
"""
is_strided(A::StridedArray) = true
is_strided(A) = false
function is_strided(A::AbstractArray)
    M = parentmodule(typeof(A))
    if parent(A) === A # SparseMatrix, StaticArray, etc
        false
    elseif M === Base || M === Core || M ===LinearAlgebra
        # bad reshapes, etc, plus Diagonal, UpperTriangular, etc.
        false
    else
        is_strided(parent(A)) # PermutedDimsArray, NamedDimsArray
    end
end

is_strided(A::BatchedAdjoint) = eltype(A) <: Real && is_strided(parent(A))
is_strided(A::BatchedTranspose) = is_strided(parent(A))

is_strided(A::LinearAlgebra.Transpose) = is_strided(parent(A))
is_strided(A::LinearAlgebra.Adjoint) = eltype(A) <: Real && is_strided(parent(A))
# This needs Compat 3.14, for any Julia < 1.6

are_strided(As...) = mapfoldl(is_strided, &, As; init=true)


text:
# a composite type to test for more complex element types
@eval struct RGB{T}
    r::T
    g::T
    b::T
end

n = 256

@testset "constructors" begin
    # static
    @on_device @cuStaticSharedMem(Float32, 1)
    @on_device @cuStaticSharedMem(Float32, (1,2))
    @on_device @cuStaticSharedMem(Tuple{Float32, Float32}, 1)
    @on_device @cuStaticSharedMem(Tuple{Float32, Float32}, (1,2))
    @on_device @cuStaticSharedMem(Tuple{RGB{Float32}, UInt32}, 1)
    @on_device @cuStaticSharedMem(Tuple{RGB{Float32}, UInt32}, (1,2))

    # dynamic
    @on_device @cuDynamicSharedMem(Float32, 1)
    @on_device @cuDynamicSharedMem(Float32, (1, 2))
    @on_device @cuDynamicSharedMem(Tuple{Float32, Float32}, 1)
    @on_device @cuDynamicSharedMem(Tuple{Float32, Float32}, (1,2))
    @on_device @cuDynamicSharedMem(Tuple{RGB{Float32}, UInt32}, 1)
    @on_device @cuDynamicSharedMem(Tuple{RGB{Float32}, UInt32}, (1,2))

    # dynamic with offset
    @on_device @cuDynamicSharedMem(Float32, 1, 8)
    @on_device @cuDynamicSharedMem(Float32, (1,2), 8)
    @on_device @cuDynamicSharedMem(Tuple{Float32, Float32}, 1, 8)
    @on_device @cuDynamicSharedMem(Tuple{Float32, Float32}, (1,2), 8)
    @on_device @cuDynamicSharedMem(Tuple{RGB{Float32}, UInt32}, 1, 8)
    @on_device @cuDynamicSharedMem(Tuple{RGB{Float32}, UInt32}, (1,2), 8)
end

@testset "dynamic shmem" begin

@testset "statically typed" begin
    function kernel(d, n)
        t = threadIdx().x
        tr = n-t+1

        s = @cuDynamicSharedMem(Float32, n)
        s[t] = d[t]
        sync_threads()
        d[t] = s[tr]

        return
    end

    a = rand(Float32, n)
    d_a = CuArray(a)

    @cuda threads=n shmem=n*sizeof(Float32) kernel(d_a, n)
    @test reverse(a) == Array(d_a)
end

@testset "parametrically typed" begin
    @testset for T in [Int32, Int64, Float16, Float32, Float64]
        function kernel(d::CuDeviceArray{T}, n) where {T}
            t = threadIdx().x
            tr = n-t+1

            s = @cuDynamicSharedMem(T, n)
            s[t] = d[t]
            sync_threads()
            d[t] = s[tr]

            return
        end

        a = rand(T, n)
        d_a = CuArray(a)

        @cuda threads=n shmem=n*sizeof(T) kernel(d_a, n)
        @test reverse(a) == Array(d_a)
    end
end

@testset "alignment" begin
    # bug: used to generate align=12, which is invalid (non pow2)
    function kernel(v0::T, n) where {T}
        shared = @cuDynamicSharedMem(T, n)
        @inbounds shared[Cuint(1)] = v0
        return
    end

    n = 32
    typ = typeof((0f0, 0f0, 0f0))
    @cuda shmem=n*sizeof(typ) kernel((0f0, 0f0, 0f0), n)
end

end

@testset "static shmem" begin

@testset "statically typed" begin
    function kernel(d, n)
        t = threadIdx().x
        tr = n-t+1

        s = @cuStaticSharedMem(Float32, 1024)
        s2 = @cuStaticSharedMem(Float32, 1024)  # catch aliasing

        s[t] = d[t]
        s2[t] = 2*d[t]
        sync_threads()
        d[t] = s[tr]

        return
    end

    a = rand(Float32, n)
    d_a = CuArray(a)

    @cuda threads=n kernel(d_a, n)
    @test reverse(a) == Array(d_a)
end

@testset "parametrically typed" begin
    @testset for typ in [Int32, Int64, Float16, Float32, Float64]
        function kernel(d::CuDeviceArray{T}, n) where {T}
            t = threadIdx().x
            tr = n-t+1

            s = @cuStaticSharedMem(T, 1024)
            s2 = @cuStaticSharedMem(T, 1024)  # catch aliasing

            s[t] = d[t]
            s2[t] = d[t]
            sync_threads()
            d[t] = s[tr]

            return
        end

        a = rand(typ, n)
        d_a = CuArray(a)

        @cuda threads=n kernel(d_a, n)
        @test reverse(a) == Array(d_a)
    end
end

@testset "alignment" begin
    # bug: used to generate align=12, which is invalid (non pow2)
    function kernel(v0::T) where {T}
        shared = CUDA.@cuStaticSharedMem(T, 32)
        @inbounds shared[Cuint(1)] = v0
        return
    end

    @cuda kernel((0f0, 0f0, 0f0))
end

end

@testset "dynamic shmem consisting of multiple arrays" begin

# common use case 1: dynamic shmem consists of multiple homogeneous arrays
#                    -> split using `view`
@testset "homogeneous" begin
    function kernel(a, b, n)
        t = threadIdx().x
        tr = n-t+1

        s = @cuDynamicSharedMem(eltype(a), 2*n)

        sa = view(s, 1:n)
        sa[t] = a[t]
        sync_threads()
        a[t] = sa[tr]

        sb = view(s, n+1:2*n)
        sb[t] = b[t]
        sync_threads()
        b[t] = sb[tr]

        return
    end

    a = rand(Float32, n)
    d_a = CuArray(a)

    b = rand(Float32, n)
    d_b = CuArray(b)

    @cuda threads=n shmem=2*n*sizeof(Float32) kernel(d_a, d_b, n)
    @test reverse(a) == Array(d_a)
    @test reverse(b) == Array(d_b)
end

# common use case 2: dynamic shmem consists of multiple heterogeneous arrays
#                    -> construct using pointer offset
@testset "heterogeneous" begin
    function kernel(a, b, n)
        t = threadIdx().x
        tr = n-t+1

        sa = @cuDynamicSharedMem(eltype(a), n)
        sa[t] = a[t]
        sync_threads()
        a[t] = sa[tr]

        sb = @cuDynamicSharedMem(eltype(b), n, n*sizeof(eltype(a)))
        sb[t] = b[t]
        sync_threads()
        b[t] = sb[tr]

        return
    end

    a = rand(Float32, n)
    d_a = CuArray(a)

    b = rand(Int64, n)
    d_b = CuArray(b)

    @cuda threads=n shmem=(n*sizeof(Float32) + n*sizeof(Int64)) kernel(d_a, d_b, n)
    @test reverse(a) == Array(d_a)
    @test reverse(b) == Array(d_b)
end

end


text:
@testset "html blocks" begin
    b = """
        {{ A }} <!-- B --> <script> C </script>
        """ |> FP.default_html_tokenizer |> FP.default_html_blockifier
    @test b[1].name == :DBB
    @test b[2].name == :COMMENT
    @test b[3].name == :SCRIPT

    @test FP.content(b[1]) == " A "
end


text:
using Immortal
using Test

@testset "Immortal.jl" begin
    # Write your tests here.
    #TODO: написать тесты
end


text:
"""
    Hamming()

Creates the Hamming distance

The Hamming distance is defined as the number of characters that do not match
"""
struct Hamming <: StringMetric end

function (dist::Hamming)(s1, s2; max_dist::Union{Integer, Nothing} = nothing)
    (s1 === missing) | (s2 === missing) && return missing
    out = abs(length(s2) - length(s1))
    for (ch1, ch2) in zip(s1, s2)
        out += ch1 != ch2
        if max_dist !== nothing
            out > max_dist && return Int(max_dist + 1)
        end
    end
    return out
end

"""
    Jaro()

Creates the Jaro distance

The Jaro distance is defined as

``1 - (m / |s1| + m / |s2| + (m - t) / m) / 3``

where ``m`` is the number of matching characters and
``t`` is half the number of transpositions.
"""
struct Jaro <: StringSemiMetric end

## http://alias-i.com/lingpipe/docs/api/com/aliasi/spell/JaroWinklerDistance.html
function (dist::Jaro)(s1, s2)
    (s1 === missing) | (s2 === missing) && return missing
    len1, len2 = length(s1), length(s2)
    if len1 > len2
        s1, s2 = s2, s1
        len1, len2 = len2, len1
    end
    # If both iterators empty, formula in Wikipedia gives 1, but it makes more sense to set it to s1 == s2
    len2 > 0 || return Float64(s1 == s2)
    d = max(0, div(len2, 2) - 1)
    flag = fill(false, len2)
    ch1_match = Vector{eltype(s1)}()
    for (i1, ch1) in enumerate(s1)
        for (i2, ch2) in enumerate(s2)
            # for each character in s1, greedy search of matching character in s2 within a distance d
            i2 >= i1 - d || continue
            i2 <= i1 + d || break
            if ch1 == ch2 && !flag[i2]
                flag[i2] = true
                push!(ch1_match, ch1)
                break
            end
        end
    end
    if isempty(ch1_match)
        return 1.0
    else
        #  m counts number matching characters
        m = length(ch1_match)
        # t/2 counts number transpositions
        t = 0
        i1 = 0
        for (i2, ch2) in enumerate(s2)
            if flag[i2]
                i1 += 1
                @inbounds t += ch2 != ch1_match[i1]
            end
        end
        return 1 - (m / len1 + m / len2 + (m - 0.5 * t) / m) / 3
    end
end

"""
    JaroWinkler(;p = 0.1, threshold = 0.3, maxlength = 4)

Creates the JaroWinkler distance

The JaroWinkler distance is defined as the Jaro distance, which is multiplied by
``(1-min(l,  maxlength) * p)`` as long as it is lower than `threshold`, and where `l` denotes the length of the common prefix.
"""
struct JaroWinkler <: StringSemiMetric
    p::Float64          # scaling factor. Default to 0.1
    threshold::Float64  # boost limit. Default to 0.3
    maxlength::Integer  # max length of common prefix. Default to 4
end

JaroWinkler(; p = 0.1, threshold = 0.3, maxlength = 4) = JaroWinkler(p, threshold, maxlength)

## http://alias-i.com/lingpipe/docs/api/com/aliasi/spell/JaroWinklerDistance.html
function (dist::JaroWinkler)(s1, s2)
    (s1 === missing) | (s2 === missing) && return missing
    out = Jaro()(s1, s2)
    if out <= dist.threshold
        l = common_prefix(s1, s2)[1]
        out = (1 - min(l, dist.maxlength) * dist.p) * out
    end
    return out
end

"""
    Levenshtein()

Creates the Levenshtein distance

The Levenshtein distance is the minimum number of operations (consisting of insertions, deletions,
substitutions of a single character) required to change one string into the other.
"""
struct Levenshtein <: StringMetric end

## Source: http://blog.softwx.net/2014/12/optimizing-levenshtein-algorithm-in-c.html
# Return max_dist + 1 if distance higher than max_dist
# to differentiate distance equal to max_dist or not, which is important for find fctions.
function (dist::Levenshtein)(s1, s2; max_dist::Union{Integer, Nothing} = nothing)
    (s1 === missing) | (s2 === missing) && return missing
    len1, len2 = length(s1), length(s2)
    if len1 > len2
        s1, s2 = s2, s1
        len1, len2 = len2, len1
    end
    if max_dist !== nothing
        len2 - len1 > max_dist && return Int(max_dist + 1)
    end
    # prefix common to both strings can be ignored
    k = common_prefix(s1, s2)
    k == len1 && return len2 - k
    # first row of matrix set to distance between "" and s2[1:i]
    v = collect(1:(len2-k))
    current = 0
    for (i1, ch1) in enumerate(s1)
        i1 > k || continue
        left = current = i1 - k - 1
        if max_dist !== nothing
            value_lb = left - 1
        end
        for (i2, ch2) in enumerate(s2)
            i2 > k || continue
            above = current
            # cost on diagonal (substitution)
            current = left
            @inbounds left = v[i2 - k]
            if ch1 != ch2
                # minimum between substitution, deletion and insertion
                current = min(current + 1, above + 1, left + 1)
            end
            if max_dist !== nothing
                value_lb = min(value_lb, left)
            end
            @inbounds v[i2 - k] = current
        end
        if max_dist !== nothing
            value_lb > max_dist && return Int(max_dist + 1)
        end
    end
    if max_dist !== nothing
        current > max_dist && return Int(max_dist + 1 )
    end
    return current
end

"""
        OptimalStringAlignement()

    Creates the OptimalStringAlignement distance (also known ad the unrestricted DamerauLevenshtein distance).

    It is the minimum number of operations (consisting of insertions,
    deletions or substitutions of a single character, or transposition of two adjacent characters)
    required to change one string into the other.

    The distance differs slightly from the Damerau-Levenshtein algorithm by imposing
    the restriction that no substring is edited more than once. So for example, "CA" to "ABC" has an edit
    distance of 2 by a complete application of Damerau-Levenshtein, but a distance of 3 by this method that
    uses the optimal string alignment algorithm. In particular, the restricted distance does not satisfy
    the triangle inequality.
"""
struct OptimalStringAlignement <: StringSemiMetric end

## http://blog.softwx.net/2015/01/optimizing-damerau-levenshtein_15.html
# Return max_dist + 1 if distance higher than max_dist
function (dist::OptimalStringAlignement)(s1, s2; max_dist::Union{Integer, Nothing} = nothing)
    (s1 === missing) | (s2 === missing) && return missing
    len1, len2 = length(s1), length(s2)
    if len1 > len2
        s1, s2 = s2, s1
        len1, len2 = len2, len1
    end
    if max_dist !== nothing
        len2 - len1 > max_dist && return Int(max_dist + 1)
    end
    k = common_prefix(s1, s2)
    k == len1 && return len2 - k
    v = collect(1:(len2 - k))
    w = similar(v)
    prevch1, prevch2 = first(s1), first(s2)
    if max_dist !== nothing
        i2_start = 0
        i2_end = max_dist
    end
    current = 0
    for (i1, ch1) in enumerate(s1)
        i1 > k || (prevch1 = ch1 ; continue)
        left = i1 - k - 1
        current = i1 - k
        nextTransCost = 0
        if max_dist !== nothing
            i2_start += i1 - k - 1 + len2 - len1 > max_dist
            i2_end += i2_end < len2
        end
        for (i2, ch2) in enumerate(s2)
            i2 > k || (prevch2 = ch2 ; continue)
            # no need to look beyond window of lower right diagonal - max distance cells
            # lower right diag is i1 - (len2 - len1)) and the upper left diagonal + max_dist cells (upper left is i1)
            if max_dist !== nothing
                (k + i2_start < i2 < 1 + k + i2_end) || (prevch2 = ch2 ; continue)
            end
            above = current
            thisTransCost = nextTransCost
            @inbounds nextTransCost = w[i2 - k]
            @inbounds w[i2 - k] = current = left
            @inbounds left = v[i2 - k]
            if ch1 != ch2
                # minimum between substitution, deletion and insertion
                current = min(current + 1, above + 1, left + 1)
                if i1 > k + 1 && i2 > k + 1 && ch1 == prevch2 && prevch1 == ch2
                    thisTransCost += 1
                    current = min(current, thisTransCost)
                end
            end
            @inbounds v[i2 - k] = current
            prevch2 = ch2
        end
        if max_dist !== nothing
            v[i1 - k + len2 - len1] > max_dist && return Int(max_dist + 1)
        end
        prevch1 = ch1
    end
    if max_dist !== nothing
        current > max_dist && return Int(max_dist + 1)
    end
    return Int(current)
end

"""
    DamerauLevenshtein()

Creates the DamerauLevenshtein distance

The DamerauLevenshtein distance is the minimum number of operations (consisting of insertions,
deletions or substitutions of a single character, or transposition of two adjacent characters)
required to change one string into the other.
"""
struct DamerauLevenshtein <: StringMetric end

# https://en.wikipedia.org/wiki/Damerau–Levenshtein_distance
# https://www.lemoda.net/text-fuzzy/damerau-levenshtein/
# Compared to Levenshtein & Restricted distance, cannot get by with only two vectors since transposition can be global
function (dist::DamerauLevenshtein)(s1, s2)
    (s1 === missing) | (s2 === missing) && return missing
    len1, len2 = length(s1), length(s2)
    if len1 > len2
        s1, s2 = s2, s1
        len1, len2 = len2, len1
    end
    k = common_prefix(s1, s2)
    k == len1 && return len2 - k
    # da[ch1] will store last spotted position of ch1 in s1
    da = Dict{eltype(s1), UInt32}()
    sizehint!(da, len1 - k)
    # distm[i1+1, i2+1] will store the distance between first i1 elements of s1 and first i2 elements of s2
    distm = zeros(UInt32, len1 + 1 - k, len2 + 1 - k)
    distm[:, 1] = 0:(len1-k)
    distm[1, :] = 0:(len2-k)
    for (i1, ch1) in enumerate(s1)
        i1 > k || continue
        # j2 is last spotted position of ch1 in s2
        # j1 will be last spotted position of ch2 in s1
        j2 = 0
        for (i2, ch2) in enumerate(s2)
            i2 > k || continue
            if ch1 == ch2
                @inbounds distm[i1 + 1 - k, i2 + 1 - k] = distm[i1 - k, i2 - k]
                j2 = i2
            else
                # minimum between substitution, deletion and insertion
                @inbounds pre = min(distm[i1 - k, i2 - k] + one(UInt32),
                                    distm[i1 + 1 - k, i2 - k] + one(UInt32),
                                    distm[i1 - k, i2 + 1 - k] + one(UInt32))
                # minimum wrt transposition --- avoid lookup if we know transposition won't be chosen
                # either because we're treating first character of s1 or because ch1 has not been spotted in s2 yet
                j1 = (i1 == k + 1 || j2 == 0) ? 0 : get(da, ch2, 0)
                if j1 > 0
                    @inbounds pre = min(pre, distm[j1 - k, j2 - k] + (i1 - j1 - 1) + 1 + (i2 - j2 - 1))
                end
                @inbounds distm[i1 + 1 - k, i2 + 1 - k] = pre
            end
        end
        da[ch1] = i1
    end
    return Int(distm[end, end])
end

"""
    RatcliffObershelp()

Creates the RatcliffObershelp distance

The distance between two strings is defined as one minus  the number of matching characters
divided by the total number of characters in the two strings. Matching characters are those
in the longest common subsequence plus, recursively, matching characters in the unmatched
region on either side of the longest common subsequence.
"""
struct RatcliffObershelp <: StringSemiMetric end

function (dist::RatcliffObershelp)(s1, s2)
    (s1 === missing) | (s2 === missing) && return missing
    len1, len2 = length(s1), length(s2)
    n_matched = length_matching_blocks(s1, s2, 1, 1, len1, len2)
    len1 + len2 == 0 ? 0.0 : 1 - 2 * n_matched / (len1 + len2)
end

function length_matching_blocks(s1, s2, start1::Integer, start2::Integer, end1::Integer, end2::Integer)
    # p is just a storage vector which will be reused
    p = zeros(Int, max(end1 - start1, end2 - start2) + 1)
    length_matching_blocks!(p, s1, s2, start1, start2, end1, end2)
end

function length_matching_blocks!(p::Vector{Int}, s1, s2, start1::Integer, start2::Integer, end1::Integer, end2::Integer)
    end1 >= start1 || return 0
    end2 >= start2 || return 0
    j1, j2, len = longest_common_pattern!(p, s1, s2, start1, start2, end1, end2)
    # exit if there is no common substring
    len == 0 && return 0
    return len +
    length_matching_blocks!(p, s1, s2, start1, start2, j1 - 1, j2 - 1) +
    length_matching_blocks!(p, s1, s2, j1 + len, j2 + len, end1, end2)
end

function longest_common_pattern!(p, s1, s2, start1, start2, end1, end2)
    if end1 - start1 > end2 - start2
        j2, j1, len = longest_common_pattern!(p, s2, s1, start2, start1, end2, end1)
    else
        j1, j2, len = 0, 0, 0
        fill!(p, 0)
        # p[i2-start2+1] stores the startingindex of the longest
        # common pattern up to i2 with prevch1 as last matching character
        for (i1, ch1) in enumerate(s1)
            i1 >= start1 || continue
            i1 <= end1 || break
            oldj2 = 0
            for (i2, ch2) in enumerate(s2)
                i2 >= start2 || continue
                i2 <= end2 || break
                if ch1 != ch2
                    newj2 = 0
                else
                    newj2 = oldj2 > 0 ? oldj2 : i2
                    newlen = i2 - newj2 + 1
                    if newlen > len
                        j1, j2, len = i1 - newlen + 1, newj2, newlen
                    end
                end
                p[i2 - start2 + 1], oldj2 = newj2, p[i2 - start2 + 1]
            end
        end
    end
    return j1, j2, len
end


text:
include("../examples/lorenz63.jl")
include("../src/lss.jl")
using LinearAlgebra
using JLD
function refine_parameter(n, n_gd_steps)
	#create epsilon orbit.
	n_runup = 5000
    s = [10., 28., 8/3]
    u0 = rand(3,1)
	u0 = lorenz63(u0, s, n_runup)[end,:,:]
	u_trj = lorenz63(u0, s, n)[:,:,1]
	n, d = size(u_trj)
	eps = 0.0025
	n_noise = 2
	noise = eps*randn(d)
	z_obs = u_trj[:,3]
	u_trj[1,:] .+= noise
	u1 = lorenz63(reshape(u_trj[1,:], d, 1), s, n-1)[:,:,1]
	u_trj = u1
	# run gradient descent to minimize observation error.
	n_gd_steps = 5000
	gamma = 1.0
	dJds = zeros(n_gd_steps)
	d_u = 2
	u = zeros(n, d, n_gd_steps)
	# compute_sens assumes many objective functions
	dJ = zeros(1,d,n)
	error = zeros(n_gd_steps)
	for i = 1:n_gd_steps
		# Set up LSS
		# J is now sum of squared observation error.
		dJ[1,3,:] = -2/n*(z_obs .- u_trj[:,3])
    	du_trj = dlorenz63(u_trj, s)
        X = perturbation(u_trj,s) #ith col in T_{u_{i+1}} M
    	#f = vectorField(u_trj,s)
    	f = zeros(d,n)
    	y, xi = lss(du_trj, X, f, s, d_u)
		dJds[i] = compute_sens(y, xi, dJ, f)[1]
    	#println(dJds[i])
		# Make u_trj the computed shadowing trajectory
		ds = gamma*dJds[i]
		s[2] = s[2] - ds
		u_trj = u_trj - ds*y'
		error[i] =  sum((u_trj[:,3] .- z_obs).^2)/n
		@show error[i]
		u[:,:,i] = u_trj
	end
	return z_obs, u[:,3,:], error
end
function refine_parameter_and_trajectory(n, n_gd_steps)
	#create epsilon orbit.
	n_runup = 5000
    s = [10., 28., 8/3]
    u0 = rand(3,1)
	u0 = lorenz63(u0, s, n_runup)[end,:,:]
	u_trj = lorenz63(u0, s, n)[:,:,1]
	n, d = size(u_trj)
	eps = 0.001
	n_noise = 2
	noise = eps*randn(d)
	z_obs = u_trj[:,3]
	u_trj[1,:] .+= noise
	u1 = lorenz63(reshape(u_trj[1,:], d, 1), s, n-1)[:,:,1]
	u_trj = u1
	# run gradient descent to minimize observation error.
	#n_gd_steps = 5000
	gamma = 0.1
	dJds = zeros(n_gd_steps)
	d_u = 2
	u = zeros(n, d, n_gd_steps)
	# compute_sens assumes many objective functions
	dJ = zeros(1,d,n)
	error = zeros(n_gd_steps)
	n_pert = 100
	flag = 1
	for i = 1:n_gd_steps
		if flag == 0
			break
		end
		# Set up LSS
		# J is now sum of squared observation error.
		dJ[1,3,:] = -2/n*(z_obs .- u_trj[:,3])
    	du_trj = dlorenz63(u_trj, s)
        X = perturbation(u_trj,s) #ith col in T_{u_{i+1}} M
    	#f = vectorField(u_trj,s)
    	f = zeros(d,n)
    	y, xi = lss(du_trj, X, f, s, d_u)
		dJds[i] = compute_sens(y, xi, dJ, f)[1]
    	#println(dJds[i])
		# Make u_trj the computed shadowing trajectory
		ds = gamma*dJds[i]
		s[2] = s[2] - ds
		#u_trj = u_trj - ds*y'
		u2 = u_trj[n_pert,:] .- ds*y[:,n_pert]
		u_trj[n_pert:end,:] = lorenz63(reshape(u2,d,1), s, n-n_pert)[:,:,1]
		error[i] =  sum((u_trj[n_pert:end,3] .- z_obs[n_pert:end]).^2)/(n-n_pert)
		@show error[i]
		if i==2
			if error[i] > error[i-1]
				flag = 0
				break
			end
		end
		u[:,:,i] = u_trj
	end
	return z_obs, u[:,3,:], error, flag
end

function assimilate()
	n_repeat = 1
	n = 2000
	n_gd_steps = 5000
	z_obs = zeros(n, n_repeat)
	z_trj = zeros(n, n_gd_steps, n_repeat)
	errors = zeros(n_gd_steps, n_repeat)
	for i = 1:n_repeat
		res1, res2, res3 = refine_parameter(n-1,
											n_gd_steps)
		z_obs[:,i] = res1
		z_trj[:,:,i] = res2
		errors[:,i] = res3
	end
	println("Errors are ")
	@show errors[:,1]
	save("../data/l63_asmln_thru_param_estn.jld",
		 "z_obs", z_obs, "z_prd", z_trj, "msq_err",
		 errors)

	return z_obs, z_trj, errors
end
function assimilate_parameter_and_trajectory()
	n_repeat = 100
	n = 2000
	n_gd_steps = 100
	z_obs = zeros(n, n_repeat)
	z_trj = zeros(n, n_gd_steps, n_repeat)
	errors = zeros(n_gd_steps, n_repeat)
	i = 1
	while i <= n_repeat
		@show i
		res1, res2, res3, flag = refine_parameter_and_trajectory(n-1,
											n_gd_steps)
		if flag==1
			z_obs[:,i] = res1
			z_trj[:,:,i] = res2
			errors[:,i] = res3
			i = i + 1
		end
	end
	println("Errors are ")
	@show errors[:,1]
	save("../data/l63_asmln_ngd200_n2000.jld",
		 "z_obs", z_obs, "z_prd", z_trj, "msq_err",
		 errors)

	return z_obs, z_trj, errors
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
module TestShowPattern
using PatternDispatch.Recode

@show (@qpat 1)
@show (@qpat 1::Int)
@show (@qpat 1::Float64)

println()
@show (@qpat x::Any)
@show (@qpat ::Any)
@show (@qpat ::Union{})
@show (@qpat (x::Union{},y::Int))

println()
@show (@qpat x)
@show (@qpat ::Int)
@show (@qpat x::Int)
@show (@qpat x~y)
@show (@qpat x~y::Int)

println()
@show (@qpat (x,y))
@show (@qpat (x,y)::Tuple)
@show (@qpat z~(x::Int,1))

end


text:
"builds the hamiltonian of BLG with a random model implemented
on the edges of the system see `rectangle` we also add two possible SC configurations:
A and B in the manuscript.
The amount of disorder is parametrized by a value `η` which determines the average number
of atom vacancies in the edges.
    `η = 0` -> crystalographic edges
we can control the rotation of the lattices w.r.t the SC leads using the keyword `θ`"

function rectangle_randombounds_sc(p, θ = 0, η = 0.; selfy = false, ϕ0 = 0)
    (; Ls, Δ, Ws) = p
    (; model0, field!, modelinter) = modelS(p)
    lat_top, lat_bot = latBLG(p, θ)

    # LOCAL SELF-ENERGY MODEL
    xmin, xmax = extrema(r->r[1], sitepositions(Quantica.combine(lat_top, lat_bot)))
    ymin, ymax = extrema(r->r[2], sitepositions(Quantica.combine(lat_top, lat_bot)))
    self_region(r) = ifelse(selfy == false, !(xmin+Ls <= r[1] <= xmax - Ls),
        !(xmin+Ls <= r[1] <= xmax - Ls) || !(ymin+Ws <= r[2] <= ymax - Ws))
    diagphi(ϕ) = Diagonal(SA[cis(ϕ), cis(ϕ), cis(-ϕ), cis(-ϕ)])
    sCself! = @onsite!((o, r; ϕ) -> o +
        self_region(r) * Δ * diagphi(sign(r[1])* ϕ/4) * σyτy * diagphi(sign(r[1])*ϕ/4)')

    # RANDOM MODEL (we randomly replace sites within self_energy by vacancies)
    vacancies(r) = sample([false, true], StatsBase.Weights([1-η, η])) * self_region(r)

    # HAMILTONIAN BUILD
    h_top = lat_top |> hamiltonian(model0; orbitals = Val(4)) |>
        unitcell(mincoordination = 5, region = !vacancies)
    h_bot = lat_bot |> hamiltonian(model0; orbitals = Val(4)) |>
        unitcell(mincoordination = 5, region = !vacancies)
     ph = Quantica.combine(h_top, h_bot; coupling = modelinter) |>
        parametric(field!, sCself!)
    ## checking vacancy probability
    if η !=0
        sitesinhdis = length(collect(siteindices(ph(ϕ=0))))
        println(sitesinhdis)
        normalsites, allsites =
            sitesinhord(lat_top, lat_bot, self_region, model0, modelinter)
        println("vacancy probability: ", 1-(sitesinhdis-normalsites)/(allsites-normalsites))
    else nothing
    end
    return ph(ϕ = ϕ0)
end

"""
    sitesinhord(latt, latb, self_region, model, modelinter)
returns the number of normal sites and the total number of sites in the non-disordered lattice
    see: rectangle_randombounds_sc()
"""
function sitesinhord(latt, latb, self_region, model, modelinter)
    h_top = latt |> hamiltonian(model; orbitals = Val(4))
    h_bot = latb |> hamiltonian(model; orbitals = Val(4))
    normalsites = length(collect(siteindices(Quantica.combine(h_top |>
        unitcell(mincoordination = 5, region = !self_region),
            h_bot |> unitcell(mincoordination = 5, region = !self_region);
                coupling = modelinter))))
    allsites = length(collect(siteindices(Quantica.combine(h_top |>
        unitcell(mincoordination = 5), h_bot |>
            unitcell(mincoordination = 5); coupling = modelinter))))
    return normalsites, allsites
end

# Deprecated terms in H (explicit construction of the leads)
# ONSITE SC MODEL (transparencies... disabled by default)
# SCo! = @onsite!((o, r) -> o + smooth_method(r) * Δ * σyτy)
# SCh! = @hopping!((t, r, dr) -> t * (1-smooth_method(r));
#     sublats = (:At, :Bt ,:Ab, :Bb, :A, :B) .=> (:At, :Bt ,:Ab, :Bb, :A, :B))
# SCτ! = @hopping!((t, r, dr) -> t * τ; region = (r, dr) -> xor(regionS(r - dr/2),
#     regionS(r + dr/2)))
    # contact model (Josephson vs rectangular mask architectures)
    # smooth_sides(r) = iszero(d) ? ifelse(abs(r[1])<Ln/2 && abs(r[2])<W/2, 0.0, 1.0) :
    #         1 - (1 + tanh((Ln/2-abs(r[1]))/d))/2
    # smooth_rectangle(r) = iszero(d) ? ifelse(abs(r[1])<Ln/2 && abs(r[2])<W/2, 0.0, 1.0) :
    #             1 - (1 + tanh((Ln/2-abs(r[1]))/d))*(1 + tanh((W/2-abs(r[2]))/d))/4
    # smooth_method(r) = ifelse(sidecontacts == true , smooth_sides(r), smooth_rectangle(r))


text:
using InvertibleNetworks, LinearAlgebra, Test
using PyPlot, Flux, Random, JLD
import Flux.Optimise.update!

# Random seed
Random.seed!(100)

# Load MNIST
X = load("../mnist_training.jld")["X"]

# Data
nx = 32
ny = 32
n_in = 1
batchsize = 10
k = 3   # kernel size
s = 1   # stride
p = 1   # padding

# Network
nscales = 2
steps_per_scale = 10
H = NetworkHyperbolic(nx, ny, n_in, batchsize, nscales, steps_per_scale; α=1f-1, hidden_factor=8, ncenter=10)

#############################################################################################################
# Training

# Loss
function loss(H, X)
    Y, logdet = H.forward(X)
    f = -log_likelihood(Y) - logdet
    ΔY = -∇log_likelihood(Y)
    H.backward(ΔY, Y)
    return f
end

# Training
ntrain = 60000
maxiter = 1000
opt = Flux.ADAM(1f-2)
fval = zeros(Float32, maxiter)
Params = get_params(H)

for j=1:maxiter

    # Evaluate objective and gradients
    idx = randperm(ntrain)[1:batchsize]
    fval[j] = loss(H, X[:,:,:,idx])
    mod(j, 10) == 0 && (print("Iteration: ", j, "; f = ", fval[j], "\n"))

    # Update params
    for p in Params
        update!(opt, p.data, p.grad)
    end
    clear_grad!(H)
end

# Test
testsize = 1
idx = randperm(ntrain)[1]
Xi = X[:,:,:,idx:idx]
Yi = H.forward(Xi)[1]

Yi_ = randn(Float32, nx, ny, n_in, testsize)
Xi_ = H.inverse(Yi_)

figure(figsize=[8,8])
i=1
subplot(2,2,1); imshow(Xi[:,:,1,1]); title(L"Data space: $x \sim \hat{p}_X$")
subplot(2,2,2); imshow(Yi[:,:,1,1]); title(L"Latent space: $z = f(x)$")
subplot(2,2,3); imshow(Xi_[:,:,1,i]); title(L"Data space: $x = f^{-1}(z)$")
subplot(2,2,4); imshow(Yi_[:,:,1,i]); title(L"Latent space: $z \sim \hat{p}_Z$")

max_stars_repo_name: slimgroup/LoopUnrolledSeismicImaging

max_issues_repo_name: slimgroup/LoopUnrolledSeismicImaging

max_forks_repo_name: slimgroup/LoopUnrolledSeismicImaging


text:
struct FINITO_adaptive_iterable{R<:Real,C<:RealOrComplex{R},Tx<:AbstractArray{C},Tf,Tg}
    F::Array{Tf}            # smooth term
    g::Tg                   # nonsmooth term
    x0::Tx                  # initial point
    N::Int                  # of data points in the finite sum problem
    L::Maybe{Union{Array{R},R}} # Lipschitz moduli of nabla f_i
    tol::R                  # coordinate-wise tolerance
    tol_b::R                # tolerance for the adaptive part
    sweeping::Int8          # 1, 2, 3 for rand, cyclic, shuffled
    α::R                    # in (0, 1), e.g.: 0.99
end

mutable struct FINITO_adaptive_state{R<:Real,Tx}
    s::Array{Tx}            # table of x_j stacked as array of arrays
    ∇f::Array{Tx}           # table of gradients
    γ::Array{R}             # stepsize parameter
    hat_γ::R                # average γ
    indr::Array{Int}        # coordinate selection index set
    fi_x::Array{R}          # value of smooth term
    av::Tx                  # the running average
    z::Tx                   # zbar
    # some extra placeholders
    τ::R                    # stepsize for linesearch
    res::Tx                 # residual (for termination)
    γ_b::R
    ind::Array{Int}         # remaining coordinates (for termination)
    idx::Int                # idx to be updated
    idxr::Int               # running idx
end

function FINITO_adaptive_state(
    s::Array{Tx},
    ∇f::Array{Tx},
    γ::Array{R},
    hat_γ::R,
    indr,
    fi_x,
    av,
    z,
) where {R,Tx}
    return FINITO_adaptive_state{R,Tx}(
        s,
        ∇f,
        γ,
        hat_γ,
        indr,
        fi_x,
        av,
        z,
        R(1.0),
        zero(av),
        R(0.0),
        copy(indr),
        Int(0),
        Int(0),
    )
end

function Base.iterate(iter::FINITO_adaptive_iterable{R,C,Tx}) where {R,C,Tx}
    N = iter.N
    ind = collect(1:N) # full index set
    # computing the gradients and updating the table s
    s = Vector{Tx}(undef, 0)
    ∇f = fill(iter.x0, (N,))
    fi_x = zeros(R, N)    # compute the cost for the case of lineasearch
    for i = 1:N
        ∇f[i], fi_x[i] = gradient(iter.F[i], iter.x0)
        push!(s, copy(iter.x0)) # table of x_i
    end
    # updating the stepsize
    γ = zeros(R, N)
    for i = 1:N
        L_int = zeros(N)
        xeps = iter.x0 .+ one(R)
        grad_f_xeps, f_xeps = gradient(iter.F[i], xeps)
        nmg = norm(grad_f_xeps - ∇f[i])
        t = 1
        while nmg < eps(R)  # in case xeps has the same gradient
            println("initial upper bound for L too small")
            xeps = iter.x0 .+ rand(t * [-1, 1], size(iter.x0))
            grad_f_xeps, f_xeps = gradient(iter.F[i], xeps)
            nmg = norm(grad_f_xeps - ∇f[i])
            t *= 2
        end
        # decide how to initialize!
        L_int[i] = nmg / (t * sqrt(length(iter.x0)))
        L_int[i] /= iter.N
        γ[i] = iter.α / (L_int[i])
    end
    #initializing the vectors
    hat_γ = 1 / sum(1 ./ γ)
    av = hat_γ * (sum(s ./ γ) - sum(∇f) / (length(ind))) # the running average
    z, ~ = prox(iter.g, av, hat_γ)

    state = FINITO_adaptive_state(s, ∇f, γ, hat_γ, ind, fi_x, av, z)

    return state, state
end

function Base.iterate(
    iter::FINITO_adaptive_iterable{R},
    state::FINITO_adaptive_state{R},
) where {R}

    # select an index
    if iter.sweeping == 1 # uniformly random
        state.idxr = rand(1:iter.N)
    elseif iter.sweeping == 2  # cyclic
        state.idxr = mod(state.idxr, iter.N) + 1
    elseif iter.sweeping == 3  # shuffled cyclic
        if state.idx == iter.N
            state.ind = randperm(iter.N)
            state.idx = 1
        else
            state.idx += 1
        end
        state.idxr = state.ind[state.idx]
    end

    @. state.res = state.z - state.s[state.idxr]
    # backtrack γ (warn if γ gets too small)
    while true
        if state.γ[state.idxr] < iter.tol_b / iter.N
            @warn "parameter `γ` became too small ($(state.γ))"
            return nothing
        end
        ~, fi_z = gradient(iter.F[state.idxr], state.z)
        fi_model =
            state.fi_x[state.idxr] +
            real(dot(state.∇f[state.idxr], state.res)) +
            (0.5 * iter.N * iter.α / state.γ[state.idxr]) * (norm(state.res)^2)
        tol = 10 * eps(R) * (1 + abs(fi_z))
        R(fi_z) <= fi_model + tol && break

        state.γ_b = state.γ[state.idxr]
        state.γ[state.idxr] *= 0.8
        # update hat_γ, av, z
        state.av ./= state.hat_γ
        @. state.av += state.s[state.idxr] / state.γ[state.idxr]
        @. state.av -= state.s[state.idxr] / state.γ_b
        state.hat_γ = 1 / (1 / state.hat_γ + 1 / state.γ[state.idxr] - 1 / state.γ_b)
        state.av .*= state.hat_γ
        prox!(state.z, iter.g, state.av, state.hat_γ) # compute prox(barz)
        @. state.res = state.z - state.s[state.idxr]
    end
    # perform the main steps
    @. state.av += (state.hat_γ / state.γ[state.idxr]) * (state.z .- state.s[state.idxr])
    state.s[state.idxr] .= state.z  #update x_i
    @. state.av += (state.hat_γ / iter.N) * state.∇f[state.idxr]
    state.fi_x[state.idxr] = gradient!(state.∇f[state.idxr], iter.F[state.idxr], state.z)
    @. state.av -= (state.hat_γ / iter.N) * state.∇f[state.idxr]
    prox!(state.z, iter.g, state.av, state.hat_γ)

    return state, state
end

solution(state::FINITO_adaptive_state) = state.z

#TODO list
## the initial guess for L may be modified
## ensuring envelope is lower bounded
## minibatch not supported

max_issues_repo_name: kul-forbes/CIAOAlgorithms

max_forks_repo_name: kul-forbes/CIAOAlgorithms


text:
################################################################
# INFERENCE METHODS FOR SPN OBJECTS
################################################################

"""
Marginal inference with incomplete evidence. Queries of type P(E=e),
where E ⊆ X.
"""
function marginalInference!(spn::SumProductNetwork, x::AbstractVector, e::BitVector)
	setInput!(spn, x, e)
	return eval!(spn)
end

"""
Marginal inference with complete evidence. Queries of type P(X=x).
"""
function marginalInference!(spn::SumProductNetwork, x::AbstractVector)
	d = length(x)
	e::BitVector = trues(d)
	return marginalInference!(spn, x, e)
end

"""
Conditional inference. Queries of type P(Q=q | E=e).

We have P(q|e) = P(q,e)/P(e), so it can be computed with two upward passes.
Because we compute probabilities in logspace, we have
log(P(q|e)) = log(P(q,e)/P(e)) = log(P(q,e)) - log(P(e)).
"""
function conditionalInference!(spn::SumProductNetwork, x::AbstractVector, q::BitVector, e::BitVector)
	if !(length(x) == length(q) == length(e))
		throw(ArgumentError("the input x, the query bitvector, and the evidence bitvector must have the same length"))
	elseif sum(q .& e) > 0
		throw(ArgumentError("a query variable is also part of the evidence"))
	end
	qe = q .| e
	log_joint_prob = marginalInference!(spn, x, qe)
	log_evidence_prob = marginalInference!(spn, x, e)
	return log_joint_prob - log_evidence_prob
end

"""
Given values for some set of variables E ⊆ X, find the most
likely assignment to all remaining variables Q = X\E.
"""
function mpeInference!(spn::SumProductNetwork, x::AbstractVector, e::BitVector)

end


text:
function raise(errcode, message_buffer)
    if errcode[] != 0
        if errcode[] == 1
            error("CoolProp: ", unsafe_string(convert(Ptr{UInt8}, pointer(message_buffer))))
        elseif errcode[] == 2
            error("CoolProp: message buffer too small")
        else # == 3
            error("CoolProp: unknown error")
        end
    end
end
# ---------------------------------
#        Low-level access
# ---------------------------------

"""
    AbstractState_factory(backend::AbstractString, fluids::AbstractString)

Generate an AbstractState instance return an integer handle to the state class generated to be used in the other low-level accessor functions.

# Arguments
* `backend`: The backend you will use could be: `["HEOS", "REFPROP", "INCOMP", "IF97", "TREND", "HEOS&TTSE", "HEOS&BICUBIC", "SRK", "PR", "VTPR"]` etc.
* `fluids`: '&' delimited list of fluids. To get a list of possible values call `get_global_param_string(key)` with `key` one of the following: `["FluidsList", "incompressible_list_pure", "incompressible_list_solution", "mixture_binary_pairs_list", "predefined_mixtures"]`, also there is a list in CoolProp online documentation [List of Fluids](http://www.coolprop.org/fluid_properties/PurePseudoPure.html#list-of-fluids), or simply type `?CoolProp_fluids`

# Example
```julia
julia> HEOS = AbstractState_factory("HEOS", "R245fa");
julia> TTSE = AbstractState_factory("HEOS&TTSE", "R245fa");
julia> BICU = AbstractState_factory("HEOS&BICUBIC", "R245fa");
julia> SRK = AbstractState_factory("SRK", "R245fa");
julia> PR = AbstractState_factory("PR", "R245fa");
```
"""
function AbstractState_factory(backend::AbstractString, fluids::AbstractString)
    AbstractState = ccall( (:AbstractState_factory, "CoolProp"), Clong, (Cstring, Cstring, Ref{Clong}, Ptr{UInt8}, Clong), backend, fluids, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return AbstractState
end

"""
    AbstractState_free(handle::Clong)

Release a state class generated by the low-level interface wrapper.

# Arguments
* `handle`: The integer handle for the state class stored in memory
"""
function AbstractState_free(handle::Clong)
    ccall( (:AbstractState_free, "CoolProp"), Void, (Clong, Ref{Clong}, Ptr{UInt8}, Clong), handle, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return nothing
end

"""
    AbstractState_set_fractions(handle::Clong, fractions::Array{Float64})

Set the fractions (mole, mass, volume) for the AbstractState.

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `fractions`: The array of fractions

# Example
```julia
julia> handle = AbstractState_factory("HEOS", "Water&Ethanol");
julia> pq_inputs = get_input_pair_index("PQ_INPUTS");
julia> t = get_param_index("T");
julia> AbstractState_set_fractions(handle, [0.4, 0.6]);
julia> AbstractState_update(handle, pq_inputs, 101325, 0);
julia> AbstractState_keyed_output(handle, t)
352.3522212991724
julia> AbstractState_free(handle);
```
"""
function AbstractState_set_fractions(handle::Clong, fractions::Array{Float64})
    ccall( (:AbstractState_set_fractions, "CoolProp"), Void, (Clong, Ptr{Cdouble}, Clong, Ref{Clong}, Ptr{UInt8}, Clong), handle, fractions, length(fractions), errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return nothing
end

"""
    AbstractState_update(handle::Clong, input_pair::Clong, value1::Real, value2::Real)
    AbstractState_update(handle::Clong, input_pair::AbstractString, value1::Real, value2::Real)

Update the state of the AbstractState.

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `input_pair::Clong`: The integer value for the input pair obtained from get_input_pair_index(param::AbstractString)
* `input_pair::AbstractString`: The name of an input pair
* `value1`: The first input value
* `value2`: The second input value

# Example
```julia
julia> handle = AbstractState_factory("HEOS", "Water&Ethanol");
julia> pq_inputs = get_input_pair_index("PQ_INPUTS");
julia> t = get_param_index("T");
julia> AbstractState_set_fractions(handle, [0.4, 0.6]);
julia> AbstractState_update(handle, pq_inputs, 101325, 0);
julia> AbstractState_keyed_output(handle, t)
352.3522212991724
julia> AbstractState_free(handle);
```
"""
function AbstractState_update(handle::Clong, input_pair::Clong, value1::Real, value2::Real)
    ccall( (:AbstractState_update, "CoolProp"), Void, (Clong, Clong, Cdouble, Cdouble, Ref{Clong}, Ptr{UInt8}, Clong), handle, input_pair, value1, value2, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return nothing
end

function AbstractState_update(handle::Clong, input_pair::AbstractString, value1::Real, value2::Real)
    AbstractState_update(handle::Clong, get_input_pair_index(input_pair), value1::Real, value2::Real)
    return nothing
end

"""
    AbstractState_keyed_output(handle::Clong, param::Clong)

Get an output value from the `AbstractState` using an integer value for the desired output value.

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `param::Clong`: param The integer value for the parameter you want

# Note
See `AbstractState_output`
"""
function AbstractState_keyed_output(handle::Clong, param::Clong)
    output = ccall( (:AbstractState_keyed_output, "CoolProp"), Cdouble, (Clong, Clong, Ref{Clong}, Ptr{UInt8}, Clong), handle, param, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    if output == -Inf
        error("CoolProp: no correct state has been set with AbstractState_update")
    end
    return output
end

"""
    AbstractState_output(handle::Clong, param::AbstractString)

Get an output value from the `AbstractState` using an integer value for the desired output value. It is a convenience function that call `AbstractState_keyed_output`

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `param::AbstractString`: The name for the parameter you want
"""
function AbstractState_output(handle::Clong, param::AbstractString)
    return AbstractState_keyed_output(handle, get_param_index(param))
end

"""
    AbstractState_specify_phase(handle::Clong, phase::AbstractString)

Specify the phase to be used for all further calculations.

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `phase`: The string with the phase to use. Possible candidates are listed bellow:

Phase name                 |Condition
:--------------------------|-----------------------
phase_liquid               |
phase_gas                  |
phase_twophase             |
phase_supercritical        |
phase_supercritical_gas    |p < pc, T > Tc
phase_supercritical_liquid |p > pc, T < Tc
phase_critical_point       |p = pc, T = Tc
phase_unknown              |
phase_not_imposed          |

# Example
```julia
julia> heos = AbstractState_factory("HEOS", "Water");
# Do a flash call that is a very low density state point, definitely vapor
julia> @time AbstractState_update(heos, "DmolarT_INPUTS", 1e-6, 300);
  0.025233 seconds (5.23 k allocations: 142.283 KB)
# Specify the phase - for some inputs (especially density-temperature), this will result in a
# more direct evaluation of the equation of state without checking the saturation boundary
julia> AbstractState_specify_phase(heos, "phase_gas");
# We try it again - a bit faster
julia> @time AbstractState_update(heos, "DmolarT_INPUTS", 1e-6, 300);
  0.000050 seconds (5 allocations: 156 bytes)
julia> AbstractState_free(heos);
```
"""
function AbstractState_specify_phase(handle::Clong, phase::AbstractString)
    ccall( (:AbstractState_specify_phase, "CoolProp"), Void, (Clong, Cstring, Ref{Clong}, Ptr{UInt8}, Clong), handle, phase, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return nothing
end

"""
    AbstractState_unspecify_phase(handle::Clong)

Unspecify the phase to be used for all further calculations.

# Arguments
* `handle`: The integer handle for the state class stored in memory
"""
function AbstractState_unspecify_phase(handle::Clong)
    ccall( (:AbstractState_unspecify_phase, "CoolProp"), Void, (Clong, Ref{Clong}, Ptr{UInt8}, Clong), handle, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return nothing
end

"""
    AbstractState_update_and_common_out(handle::Clong, input_pair::Clong, value1::Array{Float64}, value2::Array{Float64}, length::Integer, T::Array{Float64}, p::Array{Float64}, rhomolar::Array{Float64}, hmolar::Array{Float64}, smolar::Array{Float64})
    AbstractState_update_and_common_out(handle::Clong, input_pair::AbstractString, value1::Array{Float64}, value2::Array{Float64}, length::Integer, T::Array{Float64}, p::Array{Float64}, rhomolar::Array{Float64}, hmolar::Array{Float64}, smolar::Array{Float64})

Update the state of the AbstractState and get an output value five common outputs (temperature, pressure, molar density, molar enthalpy and molar entropy) from the AbstractState using pointers as inputs and output to allow array computation.

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `input_pair::Clong`: The integer value for the input pair obtained from get_input_pair_index
* `input_pair::AbstractString`:
* `value1`: The pointer to the array of the first input parameters
* `value2`: The pointer to the array of the second input parameters
* `length`: The number of elements stored in the arrays (both inputs and outputs MUST be the same length)
* `T`: The pointer to the array of temperature
* `p`: The pointer to the array of pressure
* `rhomolar`: Array of molar density
* `hmolar`: The array of molar enthalpy
* `smolar`: Trray of molar entropy

# Example
```julia
julia> handle = AbstractState_factory("HEOS", "Water&Ethanol");
julia> pq_inputs = get_input_pair_index("PQ_INPUTS");
julia> AbstractState_set_fractions(handle, [0.4, 0.6]);
julia> T = [0.0]; p = [0.0]; rhomolar = [0.0]; hmolar = [0.0]; smolar = [0.0];
julia> AbstractState_update_and_common_out(handle, pq_inputs, [101325.0], [0.0], 1, T, p, rhomolar, hmolar, smolar);
julia> AbstractState_free(handle);
```
"""
function AbstractState_update_and_common_out(handle::Clong, input_pair::Clong, value1::Array{Float64}, value2::Array{Float64}, length::Integer, T::Array{Float64}, p::Array{Float64}, rhomolar::Array{Float64}, hmolar::Array{Float64}, smolar::Array{Float64})
    ccall( (:AbstractState_update_and_common_out, "CoolProp"), Void, (Clong, Clong, Ref{Cdouble}, Ref{Cdouble}, Clong, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Clong}, Ptr{UInt8}, Clong), handle, input_pair, value1, value2, length, T, p, rhomolar, hmolar, smolar, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return T, p, rhomolar, hmolar, smolar
end

function AbstractState_update_and_common_out(handle::Clong, input_pair::AbstractString, value1::Array{Float64}, value2::Array{Float64}, length::Integer, T::Array{Float64}, p::Array{Float64}, rhomolar::Array{Float64}, hmolar::Array{Float64}, smolar::Array{Float64})
    return AbstractState_update_and_common_out(handle, get_input_pair_index(input_pair), value1, value2, length, T, p, rhomolar, hmolar, smolar)
end

function AbstractState_update_and_common_out(handle::Clong, input_pair::Clong, value1::Array{Float64}, value2::Array{Float64}, length::Integer)
    T, p, rhomolar, hmolar, smolar = [fill(NaN,length) for i=1:5]
    return AbstractState_update_and_common_out(handle, input_pair, value1, value2, length, T, p, rhomolar, hmolar, smolar)
end

function AbstractState_update_and_common_out(handle::Clong, input_pair::AbstractString, value1::Array{Float64}, value2::Array{Float64}, length::Integer)
    return AbstractState_update_and_common_out(handle, get_input_pair_index(input_pair), value1, value2, length)
end

"""
    AbstractState_update_and_1_out(handle::Clong, input_pair::Clong, value1::Array{Float64}, value2::Array{Float64}, length::Integer, output::Clong, out::Array{Float64})
    AbstractState_update_and_1_out(handle::Clong, input_pair::AbstractString, value1::Array{Float64}, value2::Array{Float64}, length::Integer, output::AbstractString, out::Array{Float64})

Update the state of the AbstractState and get one output value (temperature, pressure, molar density, molar enthalpy and molar entropy) from the AbstractState using pointers as inputs and output to allow array computation.

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `input_pair::Clong`: The integer value for the input pair obtained from get_input_pair_index
* `input_pair::AbstractString`:
* `value1`: The pointer to the array of the first input parameters
* `value2`: The pointer to the array of the second input parameters
* `length`: The number of elements stored in the arrays (both inputs and outputs MUST be the same length)
* `output`: The indice for the output desired
* `out`: The array for output
"""
function AbstractState_update_and_1_out(handle::Clong, input_pair::Clong, value1::Array{Float64}, value2::Array{Float64}, length::Integer, output::Clong, out::Array{Float64})
    ccall( (:AbstractState_update_and_1_out, "CoolProp"), Void, (Clong, Clong, Ref{Cdouble}, Ref{Cdouble}, Clong, Clong, Ref{Cdouble}, Ref{Clong}, Ptr{UInt8}, Clong), handle, input_pair, value1, value2, length, output, out, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return out
end

function AbstractState_update_and_1_out(handle::Clong, input_pair::AbstractString, value1::Array{Float64}, value2::Array{Float64}, length::Integer, output::AbstractString, out::Array{Float64})
    return AbstractState_update_and_1_out(handle, get_input_pair_index(input_pair), value1, value2, length, get_param_index(output), out)
end

function AbstractState_update_and_1_out(handle::Clong, input_pair::Clong, value1::Array{Float64}, value2::Array{Float64}, length::Integer, output::Clong)
    out = fill(NaN,length)
    return AbstractState_update_and_1_out(handle, input_pair, value1, value2, length, output, out)
end

function AbstractState_update_and_1_out(handle::Clong, input_pair::AbstractString, value1::Array{Float64}, value2::Array{Float64}, length::Integer, output::AbstractString)
    return AbstractState_update_and_1_out(handle, get_input_pair_index(input_pair), value1, value2, length, get_param_index(output))
end

"""
    AbstractState_update_and_5_out(handle::Clong, input_pair::Clong, value1::Array{Float64}, value2::Array{Float64}, length::Integer, outputs::Array{Clong}, out1::Array{Float64}, out2::Array{Float64}, out3::Array{Float64}, out4::Array{Float64}, out5::Array{Float64})
    AbstractState_update_and_5_out{S<:AbstractString}(handle::Clong, input_pair::AbstractString, value1::Array{Float64}, value2::Array{Float64}, length::Integer, outputs::Array{S}, out1::Array{Float64}, out2::Array{Float64}, out3::Array{Float64}, out4::Array{Float64}, out5::Array{Float64})

Update the state of the AbstractState and get an output value five common outputs (temperature, pressure, molar density, molar enthalpy and molar entropy) from the AbstractState using pointers as inputs and output to allow array computation.

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `input_pair::Clong`: The integer value for the input pair obtained from get_input_pair_index
* `input_pair::AbstractString`:
* `value1`: The pointer to the array of the first input parameters
* `value2`: The pointer to the array of the second input parameters
* `length`: The number of elements stored in the arrays (both inputs and outputs MUST be the same length)
* `outputs`: The 5-element vector of indices for the outputs desired
* `out1`: The array for the first output
* `out2`: The array for the second output
* `out3`: The array for the third output
* `out4`: The array for the fourth output
* `out5`: The array for the fifth output
"""
function AbstractState_update_and_5_out(handle::Clong, input_pair::Clong, value1::Array{Float64}, value2::Array{Float64}, length::Integer, outputs::Array{Clong}, out1::Array{Float64}, out2::Array{Float64}, out3::Array{Float64}, out4::Array{Float64}, out5::Array{Float64})
    ccall( (:AbstractState_update_and_5_out, "CoolProp"), Void, (Clong, Clong, Ref{Cdouble}, Ref{Cdouble}, Clong, Ref{Clong}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Clong}, Ptr{UInt8}, Clong), handle, input_pair, value1, value2, length, outputs, out1, out2, out3, out4, out5, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return out1, out2, out3, out4, out5
end

function AbstractState_update_and_5_out{S<:AbstractString}(handle::Clong, input_pair::AbstractString, value1::Array{Float64}, value2::Array{Float64}, length::Integer, outputs::Array{S}, out1::Array{Float64}, out2::Array{Float64}, out3::Array{Float64}, out4::Array{Float64}, out5::Array{Float64})
    outputs_key = [get_param_index(outputs[k]) for k = 1:5]
    return AbstractState_update_and_5_out(handle, get_input_pair_index(input_pair), value1, value2, length, outputs_key, out1, out2, out3, out4, out5)
end

function AbstractState_update_and_5_out(handle::Clong, input_pair::Clong, value1::Array{Float64}, value2::Array{Float64}, length::Integer, outputs::Array{Clong})
    out1, out2, out3, out4, out5 = [fill(NaN,length) for i=1:5]
    return AbstractState_update_and_5_out(handle, input_pair, value1, value2, length, outputs, out1, out2, out3, out4, out5)
end

function AbstractState_update_and_5_out{S<:AbstractString}(handle::Clong, input_pair::AbstractString, value1::Array{Float64}, value2::Array{Float64}, length::Integer, outputs::Array{S})
    outputs_key = [get_param_index(outputs[k]) for k = 1:5]
    return AbstractState_update_and_5_out(handle, get_input_pair_index(input_pair), value1, value2, length, outputs_key)
end

"""
    AbstractState_set_binary_interaction_double(handle::Clong, i::Int, j::Int, parameter::AbstractString, value::Float64)

Set binary interraction parrameter for diffrent mixtures model e.g.: "linear", "Lorentz-Berthelot"

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `i`: indice of the first fluid of the binary pair
* `j`: indice of the second fluid of the binary pair
* `parameter`: string wit the name of the parameter, e.g.: "betaT", "gammaT", "betaV", "gammaV"
* `value`: the value of the binary interaction parameter

# Example
```julia
julia> handle = AbstractState_factory("HEOS", "Water&Ethanol");
julia> AbstractState_set_binary_interaction_double(handle, 0, 1, "betaT", 0.987);
julia> pq_inputs = get_input_pair_index("PQ_INPUTS");
julia> t = get_param_index("T");
julia> AbstractState_set_fractions(handle, [0.4, 0.6]);
julia> AbstractState_update(handle, pq_inputs, 101325, 0);
julia> AbstractState_keyed_output(handle, t)
349.32634425309755
julia> AbstractState_free(handle);
```
"""
function AbstractState_set_binary_interaction_double(handle::Clong, i::Integer, j::Integer, parameter::AbstractString, value::Real)
    ccall( (:AbstractState_set_binary_interaction_double, "CoolProp"), Void, (Clong, Clong, Clong, Cstring, Cdouble, Ref{Clong}, Ptr{UInt8}, Clong), handle, i, j, parameter, value, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return nothing
end

"""
    AbstractState_set_cubic_alpha_C(handle::Clong, i::Integer, parameter::AbstractString, c1::Real, c2::Real, c3::Real)

Set cubic's alpha function parameters.

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `i`: indice of the fluid the parramter should be applied too (for mixtures)
* `parameter`: the string specifying the alpha function to use, e.g. "TWU" for the Twu or "MC" for Mathias-Copeman alpha function.
* `c1`: the first parameter for the alpha function
* `c2`: the second parameter for the alpha function
* `c3`: the third parameter for the alpha function
"""
function AbstractState_set_cubic_alpha_C(handle::Clong, i::Integer, parameter::AbstractString, c1::Real, c2::Real, c3::Real)
    ccall( (:AbstractState_set_cubic_alpha_C, "CoolProp"), Void, (Clong, Clong, Cstring, Cdouble, Cdouble, Cdouble, Ref{Clong}, Ptr{UInt8}, Clong), handle, i, parameter, c1, c2, c3, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return nothing
end

"""
    AbstractState_set_fluid_parameter_double(handle::Clong, i::Integer, parameter::AbstractString, value::Real)

Set some fluid parameter (ie volume translation for cubic). Currently applied to the whole fluid not to components.

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `i`: indice of the fluid the parramter should be applied to (for mixtures)
* `parameter`: string wit the name of the parameter, e.g. "c", "cm", "c_m" for volume translation parrameter.
* `value`: the value of the parameter
"""
function AbstractState_set_fluid_parameter_double(handle::Clong, i::Integer, parameter::AbstractString, value::Real)
    ccall( (:AbstractState_set_fluid_parameter_double, "CoolProp"), Void, (Clong, Clong, Cstring, Cdouble, Ref{Clong}, Ptr{UInt8}, Clong), handle, i, parameter, value, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return nothing
end

"""
    AbstractState_first_saturation_deriv(handle::Clong, of::Clong, wrt::Clong)

Calculate a saturation derivative from the AbstractState using integer values for the desired parameters.

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `of`: The parameter of which the derivative is being taken
* `wrt`: The derivative with with respect to this parameter

# Example
```julia
julia> as = AbstractState_factory("HEOS", "Water");
julia> AbstractState_update(as, "PQ_INPUTS", 15e5, 0);
julia> AbstractState_first_saturation_deriv(as, get_param_index("Hmolar"), get_param_index("P"))
0.0025636362140578207
```

# Ref
double CoolProp::AbstractState_first_saturation_deriv(const long handle, const long Of, const long Wrt, long* errcode, char* message_buffer, const long buffer_length);
"""
function AbstractState_first_saturation_deriv(handle::Clong, of::Clong, wrt::Clong)
    output = ccall( (:AbstractState_first_saturation_deriv, "CoolProp"), Cdouble, (Clong, Clong, Clong, Ref{Clong}, Ptr{UInt8}, Clong), handle, of, wrt, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    if output == -Inf
        error("CoolProp: no correct state has been set with AbstractState_update")
    end
    return output
end

"""
    AbstractState_first_partial_deriv(handle::Clong, of::Clong, wrt::Clong, constant::Clong)

Calculate the first partial derivative in homogeneous phases from the AbstractState using integer values for the desired parameters.

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `of`: The parameter of which the derivative is being taken
* `Wrt`: The derivative with with respect to this parameter
* `Constant`: The parameter that is not affected by the derivative

# Example
```julia
julia> as = AbstractState_factory("HEOS", "Water");
julia> AbstractState_update(as, "PQ_INPUTS", 15e5, 0);
julia> AbstractState_first_partial_deriv(as, get_param_index("Hmolar"), get_param_index("P"), get_param_index("S"))
2.07872526058326e-5
julia> AbstractState_first_partial_deriv(as, get_param_index("Hmolar"), get_param_index("P"), get_param_index("D"))
5.900781297636475e-5
```

# Ref
double CoolProp::AbstractState_first_partial_deriv(const long handle, const long Of, const long Wrt, const long Constant, long* errcode, char* message_buffer, const long buffer_length);
"""
function AbstractState_first_partial_deriv(handle::Clong, of::Clong, wrt::Clong, constant::Clong)
    output = ccall( (:AbstractState_first_partial_deriv, "CoolProp"), Cdouble, (Clong, Clong, Clong, Clong, Ref{Clong}, Ptr{UInt8}, Clong), handle, of, wrt, constant, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    if output == -Inf
        error("CoolProp: no correct state has been set with AbstractState_update")
    end
    return output
end

"""
    AbstractState_build_phase_envelope(handle::Clong, level::AbstractString)

Build the phase envelope.

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `level`: How much refining of the phase envelope ("none" to skip refining (recommended) or "veryfine")

# Note
If there is an error in an update call for one of the inputs, no change in the output array will be made

# Ref
CoolPRop::AbstractState_build_phase_envelope(const long handle, const char* level, long* errcode, char* message_buffer, const long buffer_length);
"""
function AbstractState_build_phase_envelope(handle::Clong, level::AbstractString)
    ccall( (:AbstractState_build_phase_envelope, "CoolProp"), Void, (Clong, Cstring, Ref{Clong}, Ptr{UInt8}, Clong), handle, level, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return nothing
end

"""
    AbstractState_get_phase_envelope_data(handle::Clong, length::Integer, T::Array{Float64}, p::Array{Float64}, rhomolar_vap::Array{Float64}, rhomolar_liq::Array{Float64}, x::Array{Float64}, y::Array{Float64})

Get data from the phase envelope for the given mixture composition.

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `length`: The number of elements stored in the arrays (both inputs and outputs MUST be the same length)
* `T`: The pointer to the array of temperature (K)
* `p`: The pointer to the array of pressure (Pa)
* `rhomolar_vap`: The pointer to the array of molar density for vapor phase (m^3/mol)
* `rhomolar_liq`: The pointer to the array of molar density for liquid phase (m^3/mol)
* `x`: The compositions of the "liquid" phase (WARNING: buffer should be Ncomp*Npoints in length, at a minimum, but there is no way to check buffer length at runtime)
* `y`: The compositions of the "vapor" phase (WARNING: buffer should be Ncomp*Npoints in length, at a minimum, but there is no way to check buffer length at runtime)

# Example
```julia
julia> HEOS=AbstractState_factory("HEOS","Methane&Ethane");
julia> length=200;
julia> t=zeros(length);p=zeros(length);x=zeros(2*length);y=zeros(2*length);rhomolar_vap=zeros(length);rhomolar_liq=zeros(length);
julia> AbstractState_set_fractions(HEOS, [0.2, 1 - 0.2])
julia> AbstractState_build_phase_envelope(HEOS, "none")
julia> AbstractState_get_phase_envelope_data(HEOS, length, t, p, rhomolar_vap, rhomolar_liq, x, y)
julia> AbstractState_free(HEOS)
```

# Note
If there is an error in an update call for one of the inputs, no change in the output array will be made

# Ref
CoolProp::AbstractState_get_phase_envelope_data(const long handle, const long length, double* T, double* p, double* rhomolar_vap, double* rhomolar_liq, double* x, double* y, long* errcode, char* message_buffer, const long buffer_length);
"""
function AbstractState_get_phase_envelope_data(handle::Clong, length::Integer, T::Array{Float64}, p::Array{Float64}, rhomolar_vap::Array{Float64}, rhomolar_liq::Array{Float64}, x::Array{Float64}, y::Array{Float64})
    ccall( (:AbstractState_get_phase_envelope_data, "CoolProp"), Void, (Clong, Clong, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Clong}, Ptr{UInt8}, Clong), handle, length, T, p, rhomolar_vap, rhomolar_liq, x, y, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return T, p, rhomolar_vap, rhomolar_liq, x, y
end

function AbstractState_get_phase_envelope_data(handle::Clong, length::Integer, ncomp::Integer)
    T, p, rhomolar_vap, rhomolar_liq = [fill(NaN,length) for i=1:5]
    x, y = [fill(NaN,length*ncomp) for i=1:2]
    return AbstractState_get_phase_envelope_data(handle, length, T, p, rhomolar_vap, rhomolar_liq, x, y)
end

"""
    AbstractState_build_spinodal(handle::Clong)

Build the spinodal.

# Arguments
* `handle`: The integer handle for the state class stored in memory

# Ref
CoolProp::AbstractState_build_spinodal(const long handle, long* errcode, char* message_buffer, const long buffer_length);
"""
function AbstractState_build_spinodal(handle::Clong)
    ccall( (:AbstractState_build_spinodal, "CoolProp"), Void, (Clong, Ref{Clong}, Ptr{UInt8}, Clong), handle, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return nothing
end

"""
    AbstractState_get_spinodal_data(handle::Clong, length::Integer, tau::Array{Float64}, dalta::Array{Float64}, m1::Array{Float64})

Get data for the spinodal curve.

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `length`: The number of elements stored in the arrays (all outputs MUST be the same length)
* `tau`: The pointer to the array of reciprocal reduced temperature
* `delta`: The pointer to the array of reduced density
* `m1`: The pointer to the array of M1 values (when L1=M1=0, critical point)

# Note
If there is an error, no change in the output arrays will be made

# Example
julia> HEOS=AbstractState_factory("HEOS","Methane&Ethane");
julia> AbstractState_set_fractions(HEOS, [0.1, 0.9]);
julia> AbstractState_build_spinodal(HEOS);
julia> tau, delta, m1 = AbstractState_get_spinodal_data(HEOS, 127);

# Ref
CoolProp::AbstractState_get_spinodal_data(const long handle, const long length, double* tau, double* delta, double* M1, long* errcode, char* message_buffer, const long buffer_length);
"""
function AbstractState_get_spinodal_data(handle::Clong, length::Integer, tau::Array{Float64}, delta::Array{Float64}, m1::Array{Float64})
    ccall( (:AbstractState_get_spinodal_data, "CoolProp"), Void, (Clong, Clong, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Clong}, Ptr{UInt8}, Clong), handle, length, tau, delta, m1, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return tau, delta, m1;
end

function AbstractState_get_spinodal_data(handle::Clong, length::Integer)
    tau, delta, m1 = [fill(NaN,length) for i=1:3]
    return AbstractState_get_spinodal_data(handle, length, tau, delta, m1)
end

"""
    abstractState_all_critical_points(handle::Clong, length::Integer, T::Array{Float64}, p::Array{Float64}, rhomolar::Array{Float64}, stable::Array{Clong})

Calculate all the critical points for a given composition.

# Arguments
* `handle`: The integer handle for the state class stored in memory
* `length`: The length of the buffers passed to this function
* `T`: The pointer to the array of temperature (K)
* `p`: The pointer to the array of pressure (Pa)
* `rhomolar`: The pointer to the array of molar density (m^3/mol)
* `stable`: The pointer to the array of boolean flags for whether the critical point is stable (1) or unstable (0)

# Note
If there is an error in an update call for one of the inputs, no change in the output array will be made

# Ref
CoolProp::AbstractState_all_critical_points(const long handle, const long length, double* T, double* p, double* rhomolar, long* stable, long* errcode, char* message_buffer, const long buffer_length);
"""
function AbstractState_all_critical_points(handle::Clong, length::Integer, T::Array{Float64}, p::Array{Float64}, rhomolar::Array{Float64}, stable::Array{Clong})
    ccall( (:AbstractState_all_critical_points, "CoolProp"), Void, (Clong, Clong, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Clong}, Ref{Clong}, Ptr{UInt8}, Clong), handle, length, T, p, rhomolar, stable, errcode, message_buffer::Array{UInt8, 1}, buffer_length)
    raise(errcode, message_buffer)
    return T, p, rhomolar, stable
end

function AbstractState_all_critical_points(handle::Clong, length::Integer)
    T, p, rhomolar = [fill(NaN,length) for i=1:3]
    stable = zeros(Clong, length)
    return  AbstractState_all_critical_points(handle, length, T, p, rhomolar, stable)
end


text:
include("ops/op_fp_dd.jl")
include("ops/op_fpfp_dd.jl")
include("ops/op_fpdd_dd.jl")
include("ops/op_dd_dd.jl")
include("ops/op_ddfp_dd.jl")
include("ops/op_ddsi_dd.jl")
include("ops/op_dddd_dd.jl")
include("ops/op_fp_db.jl")
include("ops/op_db_db.jl")
include("ops/op_fpfp_db.jl")
include("ops/op_fpdb_db.jl")
include("ops/op_dbfp_db.jl")
include("ops/op_dbsi_db.jl")
include("ops/op_dbdb_db.jl")
include("ops/arith.jl")


text:
# A simple quasi-static problem, based on an L-BFGS minimization
using DiscreteElasticRods
DER = DiscreteElasticRods
using BenchmarkTools
using DelimitedFiles

# Parameters
E = 5e6
ν = 0.48
L = 0.30
R = 0.003

n_array = [10, 20, 30, 40, 50, 75, 100, 200, 500]
time_array = zeros(2, length(n_array))
mem_array = zeros(2, length(n_array))

for k = 1:length(n_array)
    ns = n_array[k]

    # Create a straight rod with NS seg
    p1 = [0., 0., 0.]
    p2 = [L, 0., 0.]
    d0 = [0., 0., 1.]
    ref_rod = DER.straight_rod(p1, p2, d0, ns)

    # Reference kinematics and properties
    props = DER.elliptical_stiffness(E, ν, R, R)
    ref_strain = DER.full_kinematics(ref_rod)

    # q ordering: [θ2, x3, θ3, x4, ..., x, θ]
    ncore = 4*ns - 7
    nq = ncore + 7

    # non - allocating
    function param_2_rod!(Δx, Δθ, q)

        # Rearrange core DOFs
        Δx[1, 3:end-2] = q[2:4:ncore]
        Δx[2, 3:end-2] = q[3:4:ncore]
        Δx[3, 3:end-2] = q[4:4:ncore]
        Δθ[2:end-1] = q[1:4:ncore]

        # Cantilever BC: do nothing
        # Free BC: fill in the end of the coords
        Δx[:,end-1] = q[ncore+1:ncore+3]
        Δθ[end] = q[ncore+4]
        Δx[:,end] = q[ncore+5:ncore+7]

    end

    T = Float64
    cache1 = zeros(T, 6, ns+1)
    cache2 = zeros(T, 24, ns)
    Δx = zeros(T, 3,ns+2)
    Δθ = zeros(T, 1,ns+1)
    rod = DER.allocate_rod(T, ns)
    strain = DER.allocate_strain(T,ns)
    caches = (zeros(6, ns+1), zeros(9, ns))
    p_cache = (Δx, Δθ, rod, strain, caches)

    function total_energy_nonalloc(q::Vector{T}) where T

        Δx, Δθ, rod, strain, caches = p_cache
        cache_6ne, cache_9ns = caches

        param_2_rod!(Δx, Δθ, q)
        delta = DER.rod_delta(Δx, Δθ)
        DER.copy!(rod, ref_rod)
        DER.rod_update!(rod, delta, cache)
        DER.full_kinematics!(strain, rod, (cache_6ne, cache_9ns))
        DER.elastic_energy(ref_strain, strain, props, cache_9ns)
    end

    # Allocating
    function param_2_rod(q)
        T = eltype(q)
        Δx = zeros(T, 3,ns+2)
        Δθ = zeros(T, 1,ns+1)
        param_2_rod!(Δx, Δθ, q)
        Δx, Δθ
    end

    function total_energy(q)
        Δx, Δθ = param_2_rod(q)
        rod_coords = DER.rod_delta(Δx, Δθ)

        new_rod = DER.rod_update(ref_rod, rod_coords)
        strains = DER.full_kinematics(new_rod)
        DER.elastic_energy(ref_strain, strains, props)
    end

    # Benchmark Both
    q = rand(nq)
    total_energy(q)
    total_energy_nonalloc(q)
    yes_alloc = @benchmark total_energy($q)
    no_alloc = @benchmark total_energy_nonalloc($q)

    time_array[1,k] = minimum(yes_alloc).time/1e9
    mem_array[1,k] = memory(yes_alloc)
    time_array[2,k] = minimum(no_alloc).time/1e9
    mem_array[2,k] = memory(no_alloc)
end

writedlm("evaluation_timing.txt", time_array, ',')
writedlm("evaluation_memory.txt", mem_array, ',')