text:
import Base: size, show, @propagate_inbounds

# Define the concept of basis functions, where each function is indexed by an integer
# For instance, (1, x, ψ0, ψ1,..., ψn) defines a basis where the index:
# 0 corresponds to the constant function
# 1 corresponds to the linear function
# n+2 corresponds to the n-th order physicist Hermite function

export Basis,
       vander!,
       vander,
       ProHermiteBasis,
       PhyHermiteBasis,
       CstProHermiteBasis,
       CstPhyHermiteBasis,
       CstLinProHermiteBasis,
       CstLinPhyHermiteBasis,
       iszerofeatureactive

"""
$(TYPEDEF)

An abstract type to define basis of functions

## Constructors
- `Basis(m)`

## Construct new basis of functions

A specific feature basis `MyBasis` must be a subtype of `Basis`, and the following routines must be implemented:
* `Base.show(io::IO, B::MyBasis)`` (optional, but desired)
* `@propagate_inbounds Base.getindex(B::MyBasis, i::Int64)`
* `vander!(dV, B::MyBasis, maxi::Int64, k::Int64, x)`
* `iszerofeatureactive(B::MyBasis) = Bool`
"""
abstract type Basis end

Base.size(B::Basis) = B.m

"""
$(TYPEDEF)

The basis composed of the probabilistic Hermite functions

## Fields
$(TYPEDFIELDS)
"""
struct ProHermiteBasis <: Basis
    m::Int64
end

@propagate_inbounds Base.getindex(B::ProHermiteBasis, i::Int64) = FamilyScaledProHermite[i]

function Base.show(io::IO, B::ProHermiteBasis)
    println(io,"Basis of "*string(B.m)*" functions: 0th -> "*string(B.m-1)*"th degree Probabilistic Hermite function")
end

"""
$(TYPEDEF)

The basis composed of the physicist Hermite functions

## Fields
$(TYPEDFIELDS)
"""
struct PhyHermiteBasis <: Basis
    m::Int64
end

@propagate_inbounds Base.getindex(B::PhyHermiteBasis, i::Int64) = FamilyScaledPhyHermite[i]

function Base.show(io::IO, B::PhyHermiteBasis)
    println(io,"Basis of "*string(B.m)*" functions: 0th -> "*string(B.m-1)*"th degree Physicist Hermite function")
end

"""
$(TYPEDEF)

The basis composed of the constant function and probabilistic Hermite functions

## Fields
$(TYPEDFIELDS)
"""
struct CstProHermiteBasis <: Basis
    m::Int64
end

@propagate_inbounds Base.getindex(B::CstProHermiteBasis, i::Int64) = i==1 ? FamilyProPolyHermite[1] : FamilyScaledProHermite[i-1]

function Base.show(io::IO, B::CstProHermiteBasis)
    println(io,"Basis of "*string(B.m)*" functions: Constant, 0th -> "*string(B.m-2)*"th degree Probabilistic Hermite function")
end

"""
$(TYPEDEF)

The basis composed of the constant function and physicist Hermite functions

## Fields
$(TYPEDFIELDS)
"""
struct CstPhyHermiteBasis <: Basis
    m::Int64
end

@propagate_inbounds Base.getindex(B::CstPhyHermiteBasis, i::Int64) = i==1 ? FamilyProPolyHermite[1] : FamilyScaledPhyHermite[i-1]

function Base.show(io::IO, B::CstPhyHermiteBasis)
    println(io,"Basis of "*string(B.m)*" functions: Constant, 0th -> "*string(B.m-2)*"th degree Physicist Hermite function")
end

"""
$(TYPEDEF)

The basis composed of the constant function, the identity, and probabilistic Hermite functions

## Fields
$(TYPEDFIELDS)
"""
struct CstLinProHermiteBasis <: Basis
    m::Int64
end

@propagate_inbounds Base.getindex(B::CstLinProHermiteBasis, i::Int64) = i==1 ? FamilyProPolyHermite[1] : i==2 ? FamilyProPolyHermite[2] : FamilyScaledProHermite[i-2]

function Base.show(io::IO, B::CstLinProHermiteBasis)
    println(io,"Basis of "*string(B.m)*" functions: Constant, Linear, 0th -> "*string(B.m-3)*"th degree Probabilistic Hermite function")
end

"""
$(TYPEDEF)

The basis composed of the constant function, the identity, and physicist Hermite functions

## Fields
$(TYPEDFIELDS)
"""
struct CstLinPhyHermiteBasis <: Basis
    m::Int64
end

@propagate_inbounds Base.getindex(B::CstLinPhyHermiteBasis, i::Int64) = i==1 ? FamilyProPolyHermite[1] : i==2 ? FamilyProPolyHermite[2] : FamilyScaledPhyHermite[i-2]

function Base.show(io::IO, B::CstLinPhyHermiteBasis)
    println(io,"Basis of "*string(B.m)*" functions: Constant, Linear 0th -> "*string(B.m-3)*"th degree Physicist Hermite function")
end

# Implementation of vander! for the different bases
"""
    vander!(dV, B::T, m::Int64, k::Int64, x) where T <: Basis

Compute the Vandermonde matrix of the basis `B` for the vector `x`up to the m-th feature of the basis
"""

"""
    vander!(dV, B::ProHermiteBasis, maxi::Int64, k::Int64, x)

Compute the Vandermonde matrix for the vector `x`
"""
function vander!(dV, B::ProHermiteBasis, m::Int64, k::Int64, x)
    N = size(x,1)
    @assert size(dV) == (N, m+1) "Wrong dimension of the Vander matrix"
    dVshift = view(dV,:,1:m+1)
    vander!(dVshift, FamilyScaledProHermite[m+1], k, x)
    return dV
end

"""
    vander!(dV, B::PhyHermiteBasis, maxi::Int64, k::Int64, x)

Compute the Vandermonde matrix for the vector `x`
"""
function vander!(dV, B::PhyHermiteBasis, m::Int64, k::Int64, x)
    N = size(x,1)
    @assert size(dV) == (N, m+1) "Wrong dimension of the Vander matrix"
    dVshift = view(dV,:,1:m+1)
    vander!(dVshift, FamilyScaledPhyHermite[m+1], k, x)
    return dV
end

"""
    vander!(dV, B::T, m::Int64, k::Int64, x) where T <: Basis

Compute the Vandermonde matrix of the basis `B` for the vector `x`up to the m-th feature of the basis
"""
function vander!(dV, B::CstProHermiteBasis, m::Int64, k::Int64, x)
    N = size(x,1)
    @assert size(dV) == (N, m+1) "Wrong dimension of the Vander matrix"

    col0 = view(dV,:,1)
    if k==0
         fill!(col0, 1.0)
    else
         fill!(col0, 0.0)
    end
    if m == 0
        return dV
    end
    dVshift = view(dV,:,2:m+1)
    vander!(dVshift, FamilyScaledProHermite[m], k, x)
     # .= vander(B.f[maxi+1], k, x)
    return dV
end

"""
    vander!(dV, B::CstPhyHermiteBasis, maxi::Int64, k::Int64, x)

Compute the Vandermonde matrix for the vector `x`
"""
function vander!(dV, B::CstPhyHermiteBasis, m::Int64, k::Int64, x)
    N = size(x,1)
    @assert size(dV) == (N, m+1) "Wrong dimension of the Vander matrix"

    col0 = view(dV,:,1)
    if k==0
         fill!(col0, 1.0)
    else
         fill!(col0, 0.0)
    end
    if m == 0
        return dV
    end
    dVshift = view(dV,:,2:m+1)
    vander!(dVshift, FamilyScaledPhyHermite[m], k, x)
     # .= vander(B.f[maxi+1], k, x)
    return dV
end

"""
    vander!(dV, B::CstLinProHermiteBasis, maxi::Int64, k::Int64, x)

Compute the Vandermonde matrix for the vector `x`
"""
function vander!(dV, B::CstLinProHermiteBasis, m::Int64, k::Int64, x)
    N = size(x,1)
    @assert size(dV) == (N, m+1) "Wrong dimension of the Vander matrix"

    if m >= 0
        # Constant feature
        col0 = view(dV,:,1)
        if k==0
             fill!(col0, 1.0)
        else
             fill!(col0, 0.0)
        end
    end

    if m >= 1
        # Linear eature
        col1 = view(dV,:,2)
        if k==0
            copy!(col1, x)
        elseif k==1
            fill!(col1, 1.0)
        else
            fill!(col1, 0.0)
        end
    end

    if m < 2
        return dV
    end
    dVshift = view(dV,:,3:m+1)
    vander!(dVshift, FamilyScaledProHermite[m-1], k, x)
     # .= vander(B.f[maxi+1], k, x)
    return dV
end

"""
    vander!(dV, B::CstLinPhyHermiteBasis, maxi::Int64, k::Int64, x)

Compute the Vandermonde matrix for the vector `x`
"""
function vander!(dV, B::CstLinPhyHermiteBasis, m::Int64, k::Int64, x)
    N = size(x,1)
    @assert size(dV) == (N, m+1) "Wrong dimension of the Vander matrix"

    if m >= 0
        # Constant feature
        col0 = view(dV,:,1)
        if k==0
             fill!(col0, 1.0)
        else
             fill!(col0, 0.0)
        end
    end

    if m >= 1
        # Linear eature
        col1 = view(dV,:,2)
        if k==0
            copy!(col1, x)
        elseif k==1
            fill!(col1, 1.0)
        else
            fill!(col1, 0.0)
        end
    end

    if m < 2
        return dV
    end
    dVshift = view(dV,:,3:m+1)
    vander!(dVshift, FamilyScaledPhyHermite[m-1], k, x)
     # .= vander(B.f[maxi+1], k, x)
    return dV
end

vander!(dV, B::Union{PhyHermiteBasis, ProHermiteBasis}, k::Int64, x) = vander!(dV, B, B.m-1, k, x)
vander!(dV, B::Union{CstPhyHermiteBasis, CstProHermiteBasis}, k::Int64, x) = vander!(dV, B, B.m-1, k, x)
vander!(dV, B::Union{CstLinPhyHermiteBasis, CstLinProHermiteBasis}, k::Int64, x) = vander!(dV, B, B.m-1, k, x)

vander(B::Union{PhyHermiteBasis, ProHermiteBasis}, m::Int64, k::Int64, x) = vander!(zeros(size(x,1),m+1), B, m, k, x)
vander(B::Union{CstPhyHermiteBasis, CstProHermiteBasis}, m::Int64, k::Int64, x) = vander!(zeros(size(x,1),m+1), B, m, k, x)
vander(B::Union{CstLinPhyHermiteBasis, CstLinProHermiteBasis}, m::Int64, k::Int64, x) = vander!(zeros(size(x,1),m+1), B, m, k, x)

# Return the k-th derivative of the Vandermonde matrix of the basis B evaluated at the samples x
vander(B::Basis, k::Int64, x) = vander!(zeros(size(x,1),B.m), B, k, x)

"""
iszerofeatureactive(B::MyBasis) = Bool

Determine if the zeroth order feature of the basis `B` is a constant (return `true`).
For instance, the first dimension is active in a feature of index [0 1] if `B = ProHermiteBasis`, but not active if `B = CstProHermiteBasis`.
"""
iszerofeatureactive(B::Union{PhyHermiteBasis, ProHermiteBasis}) = true
iszerofeatureactive(B::Union{CstPhyHermiteBasis, CstProHermiteBasis, CstLinPhyHermiteBasis, CstLinProHermiteBasis}) = false

# """
#     vander!(dV, B, maxi, k, x)
#
# Compute the Vandermonde matrix for the vector `x`
# """
# function vander!(dV, B::Basis, maxi::Int64, k::Int64, x)
#     N = size(x,1)
#     @assert size(dV) == (N, maxi+1) "Wrong dimension of the Vander matrix"
#
#     col0 = view(dV,:,1)
#     if k==0
#          fill!(col0, 1.0)
#     else
#          fill!(col0, 0.0)
#     end
#     if maxi == 0
#         return dV
#     end
#     dVshift = view(dV,:,2:maxi+1)
#     vander!(dVshift, FamilyScaledProHermite[maxi], k, x)
#      # .= vander(B.f[maxi+1], k, x)
#     return dV
# end
# @propagate_inbounds Base.getindex(B::Basis, i::Int64) = B.f[i]

#
# function CstPhyHermiteBasis(m::Int64; scaled::Bool = true)
#     if scaled == true
#         return Basis(ntuple(i -> i==1 ? FamilyProPolyHermite[1] : FamilyScaledPhyHermite[i-1], m+2))
#     else
#         return Basis(ntuple(i -> i==1 ? FamilyProPolyHermite[1] : FamilyPhyHermite[i-1], m+2))
#     end
# end
#
# function CstLinPhyHermiteBasis(m::Int64; scaled::Bool = true)
#     if scaled == true
#         return Basis(ntuple(i -> i<=2 ? FamilyProPolyHermite[i] : FamilyScaledPhyHermite[i-2], m+3))
#     else
#         return Basis(ntuple(i -> i==1 ? FamilyProPolyHermite[i] : FamilyPhyHermite[i-2], m+3))
#     end
# end
#
# function CstLinProHermiteBasis(m::Int64; scaled::Bool = true)
#     if scaled == true
#         return Basis(ntuple(i -> i<=2 ? FamilyProPolyHermite[i] : FamilyScaledProHermite[i-2], m+3))
#     else
#         return Basis(ntuple(i -> i==1 ? FamilyProPolyHermite[i] : FamilyProHermite[i-2], m+3))
#     end
# end

# function CstProHermiteBasis(m::Int64; scaled::Bool = false)
#     f = zeros(ParamFcn, m+2)
#     # f[1] = 1.0
#     f[1] = FamilyProPolyHermite[1]
#     for i=0:m
#         f[2+i] = ProHermite(i; scaled = scaled)
#     end
#     return Basis(f)
# end
#
# function CstLinPhyHermiteBasis(m::Int64; scaled::Bool = false)
#     f = zeros(ParamFcn, m+3)
#     # f[1] = 1.0
#     f[1] = FamilyProPolyHermite[1]
#     # f[1] = x
#     f[2] = FamilyProPolyHermite[2]
#     for i=0:m
#         f[3+i] = PhyHermite(i; scaled = scaled)
#     end
#     return Basis(f)
# end
#
# function CstLinProHermiteBasis(m::Int64; scaled::Bool = false)
#     f = zeros(ParamFcn, m+3)
#     # f[1] = 1.0
#     f[1] = FamilyProPolyHermite[1]
#     # f[1] = x
#     f[2] = FamilyProPolyHermite[2]
#     for i=0:m
#         f[3+i] = ProHermite(i; scaled = scaled)
#     end
#     return Basis(f)
# end

# (F::Array{ParamFcn,1})(x::T) where {T <: Real} = map!(fi->fi(x), zeros(T, size(F,1)), F)
# (B::Basis)(x::T) where {T<:Real} = B.f(x)

# (B::Basis{m})(x::T) where {m, T<:Real} = map!(fi->fi(x), zeros(T, m), B.f)

# @propagate_inbounds Base.getindex(F::Array{T,1}, i::Int) where {T<:ParamFcn} = getindex(F,i)
# @propagate_inbounds Base.setindex!(F::Array{T,1}, v::ParamFcn, i::Int) where {T<:ParamFcn} = setindex!(F,v,i)

# @propagate_inbounds Base.getindex(B::Basis, i::Int) = getindex(B.f,i)
# @propagate_inbounds Base.setindex!(B::Basis, v::ParamFcn, i::Int) = setindex!(B.f,v,i)
#

# function vander!(dV, B::Basis{m}, maxi::Int64, k::Int64, x) where {m}
#     N = size(x,1)
#     @assert size(dV) == (N, maxi+1) "Wrong dimension of the Vander matrix"
#
#     col0 = view(dV,:,1)
#      if k==0
#          fill!(col0, 0.0)
#      else
#          fill!(col0, )
#     @inbounds for i=1:maxi+1
#         col = view(dV,:,i)
#
#         if i==1
#             if k==0
#                 fill!(col, 1.0)
#             else
#                 fill!(col , 0.0)
#             end
#         else
#             if typeof(B.f[i]) <: Union{PhyHermite, ProHermite}
#                 # Store the k-th derivative of the i-th order Hermite polynomial
#                 derivative!(col, B.f[i], k, x)
#             elseif typeof(B[i]) <: Union{PhyPolyHermite, ProPolyHermite} && degree(B[i])>0
#                     Pik = derivative(B.f[i], k)
#                 #     # In practice the first component of the basis will be constant,
#                 #     # so this is very cheap to construct the derivative vector
#                     @. col = Pik(x)
#             end
#         end
#     end
#     return dV
# end


text:
# Autogenerated wrapper script for MLIR_jll for aarch64-apple-darwin
export libMLIR

using Zlib_jll
using libLLVM_jll
JLLWrappers.@generate_wrapper_header("MLIR")
JLLWrappers.@declare_library_product(libMLIR, "@rpath/libMLIR.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, libLLVM_jll)
    JLLWrappers.@init_library_product(
        libMLIR,
        "lib/libMLIR.dylib",
        nothing,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
function gaussradau( n::Int32 )
#RADAUPTS   Gauss-Legendre-Radau Quadrature Nodes and Weights

    if ( n == 1 )
        x = ([-1.0], [2.0])
    elseif ( n == 2 )
        x = ([-1.0, 1/3], [.5, 1.5])
    else
    # Compute via GaussJacobi:
        x = gaussjacobi( n - 1, 0.0, 1.0 )
        w = x[2]; x = x[1];
        x = [-1.0, x]
        w = [2.0/n^2, w./(1.0 + x[2:end])]
        x = (x, w)
    end
    return x
end


text:
function plot(this_plot::SimplePlot, that_plot::SimplePlot, varargs...; kwargs...)
  cur_plots = [
    this_plot, that_plot, varargs...
  ]

  composite_plot = CompositePlot(cur_plots; kwargs...)

  composite_plot
end


text:
## As defined in
## Mendo, L., 2016. An asymptotically optimal Bernoulli factory for certain
## functions that can be expressed as power series. arXiv:1612.08923
function _mendoPower(f::F, a::Float64, rng::RNG) where {F<:Function,
  RNG<:AbstractRNG}
  k::Int64 = 1
  X::Bool = false
  V::Bool = false
  flips::Int64 = 0
  while true
    X = f()
    flips += 1
    X && return true, flips
    V = rand(rng) < a/k
    V && return false, flips
    k += 1
  end
end


text:
function _helper_read_edges(file,last_vertex = -1,start_vertex = 0)
  edges = Int[0,3]
  try
    edges = readdlm(file,Int)
  catch e
    print(e)
    return Union{}
  end

  if (last_vertex != -1)
    v = collect(start_vertex:last_vertex)
  else
    v = sort(unique([edges[1:end,1];edges[1:end,2]]))
  end
  vertxs = Dict()
  for vert in v
    vertxs[vert] = 1
  end
  return v,vertxs,edges
end

function read_edges_from_file_1(file;last_vertex= -1,directed = false,start_vertex = 0,key_of_edge_attr="cost")

  v,vertxs,edges = _helper_read_edges(file,last_vertex,start_vertex)
  e = ExEdge[]
  for i = 1 : size(edges)[1]
    if (haskey(vertxs,edges[i,1])) && (haskey(vertxs,edges[i,2]))
      edge = ExEdge(i,edges[i,1],edges[i,2])
      edge.attributes[key_of_edge_attr] = edges[i,3]
      push!(e,edge)
    end
  end
  new_graph = graph(Int64[],ExEdge{Int64}[],is_directed = directed)
  for vertex in v
    add_vertex!(new_graph,vertex)
  end
  for edge in e
      add_edge!(new_graph,edge)
  end
  return new_graph
end

function read_edges_from_file_no_selfloops(file;last_vertex= -1,directed = false,start_vertex = 0,key_of_edge_attr="cost")

  v,vertxs,edges = _helper_read_edges(file,last_vertex,start_vertex)
  e = ExEdge[]
  for i = 1 : size(edges)[1]
    if (haskey(vertxs,edges[i,1])) && (haskey(vertxs,edges[i,2]))
      if (edges[i,1] != edges[i,2])
        edge = ExEdge(i,edges[i,1],edges[i,2])
        edge.attributes[key_of_edge_attr] = edges[i,3]
        push!(e,edge)
      end
    end
  end
  new_graph = graph(Int64[],ExEdge{Int64}[],is_directed = directed)
  for vertex in v
    add_vertex!(new_graph,vertex)
  end
  for edge in e
      add_edge!(new_graph,edge)
  end
  return new_graph
end

function read_edges_from_file_no_parallel_e(file;last_vertex= -1,directed = false,start_vertex = 0,key_of_edge_attr="cost")

  v,vertxs,edges = _helper_read_edges(file,last_vertex,start_vertex)
  e = ExEdge[]
  already_in = Dict()
  for i = 1 : size(edges)[1]
    if (haskey(vertxs,edges[i,1])) && (haskey(vertxs,edges[i,2]))
      if (haskey(already_in,sort([edges[i,1],edges[i,2]])) == false)
        already_in[sort([edges[i,1],edges[i,2]])] = 1
        edge = ExEdge(i,edges[i,1],edges[i,2])
        edge.attributes[key_of_edge_attr] = edges[i,3]
        push!(e,edge)
      end
    end
  end
  new_graph = graph(Int64[],ExEdge{Int64}[],is_directed = directed)
  for vertex in v
    add_vertex!(new_graph,vertex)
  end
  for edge in e
      add_edge!(new_graph,edge)
  end
  return new_graph
end

function read_edges_from_file_no_selfl_or_parallel_e(file;last_vertex= -1,directed = false,start_vertex = 0,key_of_edge_attr="cost")

  v,vertxs,edges = _helper_read_edges(file,last_vertex,start_vertex)
  e = ExEdge[]
  already_in = Dict()
  for i = 1 : size(edges)[1]
    if (haskey(vertxs,edges[i,1])) && (haskey(vertxs,edges[i,2]))
      if (edges[i,1] != edges[i,2])
        if (haskey(already_in,sort([edges[i,1],edges[i,2]])) == false)
          already_in[sort([edges[i,1],edges[i,2]])] = 1
          edge = ExEdge(i,edges[i,1],edges[i,2])
          edge.attributes[key_of_edge_attr] = edges[i,3]
          push!(e,edge)
        end
      end
    end
  end
  new_graph = graph(Int64[],ExEdge{Int64}[],is_directed = directed)
  for vertex in v
    add_vertex!(new_graph,vertex)
  end
  for edge in e
      add_edge!(new_graph,edge)
  end
  return new_graph
end

function read_edges_from_file(file::AbstractString;last_vertex::Int64 = -1,directed::Bool = false,start_vertex::Int64 = 0,key_of_edge_attr::AbstractString="cost",selfloops::Bool = true,parallel_edges::Bool =true )
  if (selfloops) && (parallel_edges)
    read_edges_from_file_1(file,last_vertex = last_vertex,directed = directed,key_of_edge_attr = key_of_edge_attr)
  elseif (selfloops == false) && (parallel_edges)
    read_edges_from_file_no_selfloops(file,last_vertex = last_vertex,directed = directed,key_of_edge_attr = key_of_edge_attr)
  elseif (selfloops) && (parallel_edges == false)
    read_edges_from_file_no_parallel_e(file,last_vertex = last_vertex,directed = directed,key_of_edge_attr = key_of_edge_attr)
  else
    read_edges_from_file_no_selfl_or_parallel_e(file,last_vertex = last_vertex,directed = directed,key_of_edge_attr = key_of_edge_attr)
  end
end

function save_graph_to_file(g::GenericGraph;file::AbstractString = "graph.txt",edge_key_att::AbstractString = "")

  edgs = edges(g)
  arr = Array(Int64,length(edgs),3)
  if edge_key_att == ""
    for i=1:length(edgs)
    arr[i,1] = edgs[i].source
    arr[i,2] = edgs[i].target
    arr[i,3] = 1
    end
  else
  for i=1:length(edgs)
    arr[i,1] = edgs[i].source
    arr[i,2] = edgs[i].target
    arr[i,3] = edgs[i].attributes[edge_key_att]
  end
  end
  writedlm(file,arr," ")

end

max_stars_repo_name: johnybx/GraphsExt

max_issues_repo_name: johnybx/GraphsExt

max_forks_repo_name: johnybx/GraphsExt


text:
#####################################
#State###############################
#####################################

immutable Ket{b,T} <: State{Ket{b,T}}
	label::T
end

immutable Bra{b,T} <: State{Bra{b,T}}
	label::T
end

Ket{T}(b,label::T) = Ket{b,T}(label)
ket(b,label) = Ket(b,label)

Bra{T}(b,label::T) = Bra{b,T}(label)
bra(b,label) = Bra(b,label)

typealias Single Union(Bra, Ket)

immutable Tensor{S<:Single} <: State{S}
	states::Vector{S}
	Tensor{K<:Ket}(v::Vector{K}) = new(v)
	Tensor{B<:Bra}(v::Vector{B}) = new(v)
end

tensor() = error("no method tensor()")
tensor(s::State) = s
tensor(s::State...) = tensor(collect(s))

tensor{S<:Single}(v::Vector{S}) = Tensor{S}(v)
tensor(s::Array) = reduce(tensor, s)

for t=(:Ket,:Bra)
	@eval begin
	tensor(a::($t), b::($t)) = Tensor{($t)}([a,b])
	tensor{S<:($t)}(a::($t), b::Tensor{S}) = Tensor{($t)}([a,separate(b)])
	tensor{S<:($t)}(a::Tensor{S}, b::($t)) = Tensor{($t)}([separate(a),b])
	tensor{A<:($t), B<:($t)}(a::Tensor{A}, b::Tensor{B}) = Tensor{($t)}([separate(a),separate(b)])

	tensor{S<:($t)}(a::S, b::S) = Tensor{S}([a,b])
	tensor{S<:($t)}(a::S, b::Tensor{S}) = Tensor{S}([a,separate(b)])
	tensor{S<:($t)}(a::Tensor{S}, b::S) = Tensor{S}([separate(a),b])
	tensor{S<:($t)}(a::Tensor{S}, b::Tensor{S}) = Tensor{S}([separate(a),separate(b)])
	end
end

tensor{K<:Ket, B<:Bra}(a::State{K},b::State{B}) = error("KindMismatch: cannot perform tensor($a, $b)")
tensor{K<:Ket, B<:Bra}(b::State{B},a::State{K}) = error("KindMismatch: cannot perform tensor($b, $a)")

Tensor(args...) = tensor(args...)

#####################################
#Misc Functions######################
#####################################

in(s::Single, t::Tensor) = false
for kb=(:Ket,:Bra)
	@eval begin
	in{S1<:($kb),S2<:($kb)}(s::S1, t::Tensor{S2}) = in(s, t.states)

	bsym{b,T}(s::Type{($kb){b,T}}) = b
	bsym{b,T}(s::($kb){b,T}) = b
	bsym{b}(s::Type{($kb){b}}) = b
	bsym{b}(s::($kb){b}) = b
	bsym(s::Type{($kb)}) = Any
	bsym(s::($kb)) = Any

	labeltype{b,T}(t::Type{($kb){b,T}}) = T
	labeltype{b}(t::Type{($kb){b}}) = Any
	labeltype(t::Type{($kb)}) = Any

	kind{b,T}(s::Type{State{($kb){b,T}}}) = ($kb){b}
	kind{b}(s::Type{State{($kb){b}}}) = ($kb){b}
	kind(s::Type{State{($kb)}}) = ($kb)
	kind{b,T}(s::State{($kb){b,T}}) = ($kb){b}
	kind{b}(s::State{($kb){b}}) = ($kb){b}
	kind(s::State{($kb)}) = ($kb)
	end
end

labeltype{S}(s::Type{Tensor{S}}) = Vector{labeltype(S)}
labeltype(s::State) = labeltype(typeof(s))

label(s::Single) = s.label
bsym(s::Tensor) = map(bsym, s.states)
label(s::Tensor) = [label(i) for i in s.states]

copy{S<:Single}(s::S) = S(copy(s.label))
copy{S<:Tensor}(s::S) = S(copy(s.states))

labeltype(s::State) = labeltype(typeof(s))
in(s::State, t::State) = in(s, t.states)

eltype{S<:Single}(t::Type{S}) = S
eltype{S}(t::Type{Tensor{S}}) = S
eltype(s::State) = eltype(typeof(s))

dual(t::Type{Ket}) = Bra
dual(t::Type{Bra}) = Ket
dual{b,T}(t::Type{Ket{b,T}}) = Bra{b,T}
dual{b,T}(t::Type{Bra{b,T}}) = Ket{b,T}
dual{S}(t::Type{Tensor{S}}) = Tensor{dual(S)}

isequal{S<:Single}(a::S, b::S) = isequal(label(a), label(b))
=={S<:Single}(a::S,b::S) = label(a)==label(b)

isequal{S<:Tensor}(a::S, b::S) = isequal(a.states, b.states)
=={S<:Tensor}(a::S, b::S) = a.states==b.states

hash(s::Tensor) = sum(map(hash, s.states))

ctranspose{S<:Single}(s::S) = dual(S)(label(s))
ctranspose{S<:Tensor}(s::S) = dual(S)([ctranspose(i) for i in s.states])

getindex(s::Tensor, x) = s.states[x]
separate(s::Single) = s
separate(s::Tensor) = s.states

labeldelta(a::State, b::State) = 0 #default to 0
labeldelta(a::Single, b::Single) = label(a)==label(b) ? 1 : 0
labeldelta(a::Tensor, b::Tensor) = label(a)==label(b) ? 1 : 0

isdual(a::State, b::State) = false #default to false
isdual{b,T}(x::Bra{b,T},y::Ket{b,T}) = label(x)==label(y)
isdual{b,T}(x::Ket{b,T},y::Bra{b,T}) = isdual(y,x)
isdual{b,T}(x::Tensor{Bra{b,T}},y::Tensor{Ket{b,T}}) = label(x)==label(y)
isdual{b,T}(x::Tensor{Ket{b,T}},y::Tensor{Bra{b,T}}) = isdual(y,x)
isdual{b}(x::Tensor{Bra{b}},y::Tensor{Ket{b}}) = label(x)==label(y)
isdual{b}(x::Tensor{Ket{b}},y::Tensor{Bra{b}}) = isdual(y,x)
isdual(x::Tensor{Bra},y::Tensor{Ket}) = label(x)==label(y) && samebasis(x,y)
isdual(x::Tensor{Ket},y::Tensor{Bra}) = isdual(y,x)

for op=(:length, :endof)
	@eval ($op)(s::Tensor) = $(op)(s.states)
end

svec(bsym, arr::Array{Any}) = [ket(bsym, i) for i in arr]
svec{T}(bsym, arr::Array{T}) = map(Ket{bsym, T}, arr)

#####################################
#Show Functions######################
#####################################

reprlabel(s::Single) = "$(repr(label(s))):$(bsym(s))"
function reprlabel(s::Tensor)
	str = "$(reprlabel(s.states[1]))"
	for i=2:length(s.states)
		str = "$(str), $(reprlabel(s.states[i]))"
	end
	return str
end

show{K<:Ket}(io::IO, s::State{K}) = print(io, "| $(reprlabel(s)) $rang")
show{B<:Bra}(io::IO, s::State{B}) = print(io, "$lang $(reprlabel(s)) |")

#####################################
#Arithmetic Operations###############
#####################################

inner(x::Bra, y::Ket) = InnerProduct(x,y)
inner{b,T1,T2}(x::Bra{b,T1}, y::Ket{b,T2}) = labeldelta(x, y)
inner{b,T1,T2}(x::Tensor{Bra{b,T1}}, y::Tensor{Ket{b,T2}}) = labeldelta(x,y)
inner{b,T}(x::Bra{b}, y::Tensor{Ket{b,T}}) = labeldelta(x,y[1]) * tensor(y.states[2:end])
inner{b,T}(x::Tensor{Bra{b,T}}, y::Ket{b}) = tensor(y.states[2:end]) * labeldelta(x[1],y)
inner{K<:Ket}(a::Bra, b::Tensor{K}) = inner(a,b[1])*tensor(b[2:end]...)
inner{B<:Bra}(a::Tensor{B}, b::Ket) = inner(a[1],b)*tensor(a[2:end]...)

inner{K<:Ket}(a::Bra, b::Tensor{K}, i::Int) = inner(a, b[i])*tensor(vcat(b[1:i-1], b[i+1:end])...)
inner{B<:Bra}(a::Tensor{B}, b::Ket, i::Int) = tensor(vcat(a[1:i-1], a[i+1:end])...)*inner(a[i], b)

function inner{B<:Bra,K<:Ket,S<:Single}(a::Tensor{B}, b::Tensor{K}, i::Int, target::Type{S}=Ket)
	if target<:Ket
		return inner(a, b[i])*tensor(vcat(b[1:i-1], b[i+1:end])...)
	else
		return tensor(vcat(a[1:i-1], a[i+1:end])...)*inner(a[i], b)
	end
end

function inner{B<:Bra,K<:Ket}(a::Tensor{B}, b::Tensor{K})
	if samebasis(a,b)
		return labeldelta(a,b)
	else
		for s in a.states
			b = *(s,b)
		end
		return b
	end
end

for t=(:Ket, :Bra)
	@eval begin
	kron{A<:($t),B<:($t)}(a::State{A}, b::State{B}) = tensor(a,b)
	end
end
kron{K<:Ket,B<:Bra}(a::State{K}, b::State{B}) = OuterProduct(a,b)
kron{B<:Bra,K<:Ket}(a::State{B}, b::State{K}) = OuterProduct(b,a)

*{B<:Bra,K<:Ket}(a::State{B}, b::State{K}) = inner(a,b)
*{K<:Ket,B<:Bra}(a::State{K}, b::State{B}) = kron(a,b)


text:
#######################################################################################
# DATA Configuration: these are used to provide default values in case options in the #
# data configuration files are missing.                                               #
#######################################################################################
const DEFAULT_BUILD_SUMMARY = false  # whether to summarize text before indexing
const DEFAULT_SUMMARY_NS = 1  # Default number of sentences for a summary
const DEFAULT_STEM_WORDS = false  # whether to stem words or not

# Search
const DEFAULT_VECTORS = :bm25  # can be :count, :tf, :tfidf, :bm25, :word2vec, :glove or :conceptnet
const DEFAULT_VECTORS_TRANSFORM = :none  # can be :none, :lsa or :rp
const DEFAULT_VECTORS_DIMENSION = 1  # can be any positive Int
const DEFAULT_VECTORS_ELTYPE = Float32
const DEFAULT_SEARCH_INDEX = :hnsw  # can be :naive, :brutetree, :kdtree or :hnsw
const DEFAULT_EMBEDDINGS_KIND = :binary  # can be :binary or :text
const DEFAULT_DOC2VEC_METHOD = :boe  # can be :boe, :sif, :borep or :cpmean
const DEFAULT_SEARCH_METHOD = :exact  #can be :exact or :regex
const DEFAULT_HEURISTIC  = nothing #  i.e. :hamming, :levenshtein (nothing for no suggestions)
const DEFAULT_BM25_KAPPA = 2  # default value for BM25 κ parameter
const DEFAULT_BM25_BETA = 0.75  # default value for BM25 β parameter
const DEFAULT_SIF_ALPHA = 0.01  # default value for α parameter of the SIF doc2vec method
const DEFAULT_BOREP_DIMENSION = 1024  # default BOREP embedder output dimensionality
const DEFAULT_BOREP_POOLING_FUNCTION = :sum  # pooling function for BOREP embeddera i.e. :sum, :max
const DEFAULT_DISC_NGRAM = 2  # DisC embedder n-gram parameter
const DEFAULT_SCORE_ALPHA = 0.5  # default value of the α parameter of the score transformation
const DEFAULT_OOV_POLICY = :none  # values are :none, :large_vector
const DEFAULT_OOV_VAL = 1000  # default value for OOV vectors
# Results
const DEFAULT_RESULT_AGGREGATION_STRATEGY = :mean  # can be :minimum, :maximum, :mean, :median, :product

# Text stripping flags
const DEFAULT_TEXT_STRIP_FLAGS = strip_case | strip_punctuation | strip_articles |
                                 strip_prepositions | strip_whitespace |
                                 strip_corrupt_utf8 | strip_accents
const DEFAULT_QUERY_STRIP_FLAGS = strip_case | strip_punctuation | strip_articles |
                                  strip_prepositions | strip_whitespace |
                                  strip_corrupt_utf8 | strip_accents

# Caching options
const DEFAULT_CACHE_DIRECTORY = nothing
const DEFAULT_CACHE_COMPRESSION = "none"

# Query parsing
const DEFAULT_QUERY_PARSING_SEPARATOR=":"
const DEFAULT_QUERY_GENERATION_FIELDS=Symbol[]

# SearchEnv-related options
const DEFAULT_DB_ID_KEY = :garamond_linear_id
const DEFAULT_DATA_LOADER_NAME = :noop_loader
const DEFAULT_DATA_SAMPLER_NAME = :noop_sampler
const DEFAULT_RANKER_NAME = :noop_ranker
const DEFAULT_RECOMMENDER_NAME = :noop_recommender
const DEFAULT_INPUT_PARSER_NAME = :noop_input_parser

#################
# SEARCH ENGINE #
#################
# TODO(Corneliu): Check here if any of them can also be data configuration options
#                 i.e. max matches, edit distance, document type
const DEFAULT_MAX_MATCHES = 1  # maximum number of matches that can be retrned
const DEFAULT_RESPONSE_SIZE = DEFAULT_MAX_MATCHES  # maximum number of matches that are returned for a request
const DEFAULT_RESPONSE_PAGE = 1  # which response page i.e. batch of response_size results to display
const DEFAULT_SORT_FIELDS = Symbol[]
const DEFAULT_SORT_REVERSE = false
const DEFAULT_MAX_SUGGESTIONS = 0  # maximum number of overall suggestions
const DEFAULT_CUSTOM_WEIGHTS = Dict{Symbol,Any}()  # default custom searcher weights
const DEFAULT_PDFTOTEXT_PROGRAM = "/usr/bin/pdftotext"  # program to convert PDFs to text
const DEFAULT_MAX_EDIT_DISTANCE = 2  # maximum edit distance for which to return suggestions

# DYNAMIC Configuration: These defaults can be through run-time options of the
#                        Garamond CLI client/server utilities
const DEFAULT_LOG_LEVEL = Logging.Info
const DEFAULT_LOGGER = ConsoleLogger

# STATIC Configuration: These constants cannot be overwritten
const DEFAULT_GARAMONDRC_FILE = expanduser("~/.garamondrc.jl")
const DEFAULT_ENVOP_SLEEP_INTERVAL = 10  # 10 seconds to wait between environment operation commands
const HEURISTIC_TO_DISTANCE = Dict(  # heuristic to distance object mapping
    :levenshtein => StringDistances.Levenshtein(),
    :dameraulevenshtein => StringDistances.DamerauLevenshtein(),
    :hamming => StringDistances.Hamming(),
    :jaro => StringDistances.Jaro())
const DEFAULT_DISTANCE = HEURISTIC_TO_DISTANCE[:jaro]  # default distance
const DEFAULT_PARSER_CONFIG = nothing
const DEFAULT_INDEXABLE_FIELDS = nothing
const DEFAULT_EMBEDDABLE_FIELDS = nothing
# Dictionaries for String <=>Languages.Language / Languages.Languages <=> String
# conversion
const STR_TO_LANG = Dict("english"=>Languages.English,
                         "french"=>Languages.French,
                         "german"=>Languages.German,
                         "italian"=>Languages.Italian,
                         "finnish"=>Languages.Finnish,
                         "dutch"=>Languages.Dutch,
                         "portuguese"=>Languages.Portuguese,
                         "spanish"=>Languages.Spanish,
                         "russian"=>Languages.Russian,
                         "serbian"=>Languages.Serbian,# and Languages.Croatian
                         "swedish"=>Languages.Swedish,
                         "czech"=>Languages.Czech,
                         "polish"=>Languages.Polish,
                         "bulgarian"=>Languages.Bulgarian,
                         "esperanto"=>Languages.Esperanto,
                         "hungarian"=>Languages.Hungarian,
                         "greek"=>Languages.Greek,
                         "norwegian"=>Languages.Nynorsk,
                         "slovene"=>Languages.Slovene,
                         "romanian"=>Languages.Romanian,
                         "vietnamese"=>Languages.Vietnamese,
                         "latvian"=>Languages.Latvian,
                         "turkish"=>Languages.Turkish,
                         "danish"=>Languages.Danish,
                         "arabic"=>Languages.Arabic,
                         "persian"=>Languages.Persian,
                         "korean"=>Languages.Korean,
                         "thai"=>Languages.Thai,
                         "georgian"=>Languages.Georgian,
                         "hebrew"=>Languages.Hebrew,
                         "telugu"=>Languages.Telugu,
                         "estonian"=>Languages.Estonian,
                         "hindi"=>Languages.Hindi,
                         "lithuanian"=>Languages.Lithuanian,
                         "ukrainian"=>Languages.Ukrainian,
                         "belarusian"=>Languages.Belarusian,
                         "swahili"=>Languages.Swahili,
                         "urdu"=>Languages.Urdu,
                         "kurdish"=>Languages.Kurdish,
                         "azerbaijani"=>Languages.Azerbaijani,
                         "tamil"=>Languages.Tamil
                        )

const LANG_TO_STR = Dict((v=>k) for (k,v) in STR_TO_LANG)
const SUPPORTED_LANGUAGES = [Languages.English,
                             Languages.German,
                             Languages.Romanian,
                             Languages.French,
                             Languages.Italian,
                             Languages.Dutch]
const DEFAULT_LANGUAGE = Languages.English
const DEFAULT_LANGUAGE_STR=LANG_TO_STR[DEFAULT_LANGUAGE]
const DEFAULT_TOKENIZER = :stringanalysis
const DEFAULT_NGRAM_COMPLEXITY = 1

##########################
# OTHER USEFUL CONSTANTS #
##########################
const DEFAULT_VERSION = "0.2.0-dev"
const DEFAULT_VERSION_DATE = "2020-02-07"
const DEFAULT_VERSION_COMMIT = "b8fa82f*"


text:
# TODO: INNER AND GEOMETRIC PRODUCT DEDICATED TO BLADES
# TODO: NEW STRUTUCT VERSOR
##############################################################################
# Multivector convertion                                                     #
##############################################################################
"""
```
multivector(b::BasisBlade) (Liga function)
multivector(a::Number) (Liga function)
multivector(v::Vector{Float64}) (Liga function)
multivector(A::Blade) (Liga function)
```
This function is used to create a MultiVector type
from numbers, vectors, basis blade or Blades.

## Example (in 𝔾₃₁)

```julia-repl
julia> multivector(e1)
```
returns the multivector `(1.0e1)`.

## Example (in 𝔾₃₁)

```julia-repl
julia> multivector([1.0,2.0,3.0])
```
returns the multivector `(1.0 e1)+(2.0 e2)+(3.0 e3)`.

## Example (in 𝔾₃₁)

```julia-repl
julia> multivector(Blade([[1.0,2.0,3.0]],1,true))
```
returns the multivector `(1.0 e1)+(2.0 e2)+(3.0 e3)`.
"""
function multivector(a::BasisBlade)
	ms=spzeros(GIdLayout[4])
	if a.index<0
		ms[-a.index]=-1.0
	else
		ms[a.index]=1.0
	end
	return MultiVector(ms)
end
function multivector(a::Number)
	ms=spzeros(GIdLayout[4])
	ms[1]=Float64(a)
	return MultiVector(ms)
end

# rever estas funções
function multivector(v::Vector{Float64})
	ms=spzeros(GIdLayout[4])
	for i=1:length(v)
		ms[i+1]=v[i]
	end
	return MultiVector(ms)
end
function multivector(A::Blade)
	ms=multivector(A.vectors[1])
	for i=2:A.grade
		ms=outer(ms,multivector(A.vectors[i]))
	end
	if A.scalar≈1.0
		return ms
	else
		return A.scalar*ms
	end
end

##############################################################################
#  Blade convertion and constructor                                          #
##############################################################################
"""
```
blade(scalar,args::Vector{Float64}...) (Liga Function)
```
This function is useful for two purposes:
 1. Is a simple to create a Blade type element;
 2. Convert a Float64 vector into Blade type element.

 ## Example (in 𝔾₃₁)

 ```julia-repl
julia> blade(1.0,[1.0,2.0],[1.0,-1.0])
```
returns the Blade type element `1.0*[1.0,2.0]∧[1.0,-1.0]`.
"""
function blade(scalar,args::Vector{Float64}...)
	nargs=length(args)
	v=Vector{Vector{Float64}}(undef,nargs)
	for i=1:nargs
		v[i]=args[i]
	end
	return Blade(scalar,v,nargs,true)
end

##############################################################################
# Grade like functions                                                       #
##############################################################################
"""
```
grade(b::BasisBlade) (Liga function)
```
This function is used to determine the grade of a basis blade element

## Example (in 𝔾₃₁)

```julia-repl
julia> grade(e123)
```
return the value (Int) `3`.
"""
function grade(b::BasisBlade)
	if b.index==0
		return 0
	else
		return length(findall(LogicalBasisBlade[abs(b.index)]))
	end
end

"""
```
gradeplus(b::BasisBlade) (Liga function)
```
This function is used to determine the grade plus a Basis Blade element.
It returns an integer.
"""
function gradeplus(b::BasisBlade)
	if b.index==0
		return 0
	else
		return length(findall(LogicalBasisBlade[abs(b.index)][1:GIdLayout[1]]))
	end
end

"""
```
grademinus(b::BasisBlade) (Liga function)
```
This function is used to determine the grade minus of a basis blade element
It returns an integer.
"""
function grademinus(b::BasisBlade)
	if b.index==0
		return 0
	else
		return length(findall(LogicalBasisBlade[abs(b.index)][GIdLayout[1]+1:GIdLayout[1]+GIdLayout[2]]))
	end
end

"""
```
gradeprojection(b::BasisBlade,k::Int) (Liga function)
```
This function is used to determine the grade projection of a basis blade.
This function returns a Basis Blade element.

## Example (in 𝔾₃₁)

```julia-repl
julia> gradeprojection(-e123,3)
```
returns the Basis Blade `-e123`
"""
function gradeprojection(b::BasisBlade,k::Int)
	if grade(b)==k
		return b
	else
		return BasisBlade(0)
	end
end

"""
```
gradeprojection(b::MultiVector,k::Int) (Liga function)
```
This function is used to determine the grade projection of a multivector.
This function returns a MultiVector element.

## Example (in 𝔾₃₁)

```julia-repl
julia> gradeprojection(1.0*id+e2-2.0*e12p,3)
```
returns the multivector  `- ( 2.0 e12+ )`.
"""
function gradeprojection(b::MultiVector,k::Int)
	ind=findnz(b.comp)
	nnv=nnz(b.comp)
	ms=spzeros(GIdLayout[4])
	for i=1:nnv
		bb=gradeprojection(BasisBlade(ind[1][i]),k)
		if bb.index!=0
			ms[bb.index]=ind[2][i]
		end
	end
	return MultiVector(ms)
end

############################################################################
# Geometric product functions                                              #
############################################################################

"""
```
geometric(a::GAAbstractType,b::GAAbstractType) (Liga function)
```
This function is used to calculate the geometric product of
GAAbstractType elements (Basis Blade, MultiVector or Blades).
It returns a MultiVector type.

## Example (in 𝔾₃)

```julia-repl
julia> geometric(e12,e12)
```
returns a multivector  `-1.0id`.

Note that we can use other structs in arguments.

## Example (in 𝔾₃)

```julia-repl
julia> geometric(1.0+3.0*e12,e1)
```
returns `( 1.0 e1 ) - ( 3.0 e2)`.

## Example (in 𝔾₄₁ conformal)
```julia-repl
julia> u=2.0*id + 1.5*e1 -3.0*e2 +4.0*e∞ -3.0*e0

julia> v=3.0ep

julia> geometric(u,v)
```
returns `( 16.5 id ) + ( 6.0 e+ ) + ( 4.5 e1+ ) - ( 9.0 e2+ ) - ( 7.5 e+- )`.

You can use the infix operator ∘ instead of geometric command.

## Example (in 𝔾₃)

```julia-repl
julia> u = 1.0+3.0*e12

julia> v = e1

julia> u ∘ v
```
"""
function geometric(a::GAAbstractType,b::GAAbstractType)
	if !isa(a,MultiVector)
		a=multivector(a)
	end
	if !isa(b,MultiVector)
		b=multivector(b)
	end
	ind1=findnz(a.comp)
	ind2=findnz(b.comp)
	nnv1=nnz(a.comp)
	nnv2=nnz(b.comp)
	ms=spzeros(GIdLayout[4])
	for i=1:nnv1
		vl=ind1[2][i]
		for j=1:nnv2
			ind=BBgeoprodTable[ind1[1][i],ind2[1][j]].index
			if ind>0
				ms[ind]=ms[ind]+vl*ind2[2][j]
			end
			if ind<0
				ms[-ind]=ms[-ind]-vl*ind2[2][j]
			end
		end
	end
	return MultiVector(ms)
end

"""
```
∘(a::GAAbstractType,b::GAAbstractType) (Liga function)
```
It is an infix operator for geometric product. This function is
used to calculate the geometric product of GAAbstractType elements
(Basis Blade, MultiVector or Blades). It returns a MultiVector type.

## Example (in 𝔾₃)

```julia-repl
julia> u = 1.0+3.0*e12

julia> v = e1

julia> u ∘ v
```
"""
function Base.:∘(a::GAAbstractType,b::GAAbstractType)
	return geometric(a,b)
end

################################################################################
# inner product                                                                #
################################################################################

"""
```
inner(a::GAAbstractType,b::GAAbstractType) (Liga function)
```
This function is used to determine the inner product
of GAAbstractType elements (Basis Blade, MultiVector or Blades).

## Example (in 𝔾₃)

```julia-repl
julia> inner(1.0+3.0*e12,e1)
```
returns `- ( 3.0 e2 )`.

## Example (in 𝔾₄₁ conformal)

```julia-repl
julia> u=2.0*id + 1.5*e1 -3.0*e2 +4.0*e∞ -3.0*e0

julia> v=3.0*ep

julia> inner(u,v)
```
returns `( 16.5 id )`.

An alternative way to run this functions is using the infix
operator (⋅).
"""
function inner(a::GAAbstractType,b::GAAbstractType)
	if !isa(a,MultiVector)
		a=multivector(a)
	end
	if !isa(b,MultiVector)
		b=multivector(b)
	end
	ind1=findnz(a.comp)
	ind2=findnz(b.comp)
	nnv1=nnz(a.comp)
	nnv2=nnz(b.comp)
    ms=spzeros(GIdLayout[4])
	for i=1:nnv1
		vl=ind1[2][i]
		for j=1:nnv2
			ind=BBinnerprodTable[ind1[1][i],ind2[1][j]].index
			if ind>0
				ms[ind]=ms[ind]+vl*ind2[2][j]
			end
			if ind<0
				ms[-ind]=ms[-ind]-vl*ind2[2][j]
			end
		end
    end
    return MultiVector(ms)
end

"""
```
⋅(a::GAAbstractType,b::GAAbstractType) (Liga function)
```
This function is used to determine the inner product
of GAAbstractType elements (Basis Blade or MultiVector).

## Example (in 𝔾₄₁ conformal)

```julia-repl
julia> u=2.0*id + 1.5*e1 -3.0*e2 +4.0*e∞ -3.0*e0

julia> v=3.0*ep

julia> u⋅v
```
returns `( 16.5 id )`.
"""
function ⋅(a::GAAbstractType,b::GAAbstractType)
	return inner(a,b)
end

##############################################################################
# outer product                                                              #
##############################################################################

"""
```
outer(a::GAAbstractType,b::GAAbstractType) (Liga function)
```
This function is used to determine the outer product
of GAAbstractType elements (Basis Blade, MultiVector or Blades).

## Example (in 𝔾₃)

```julia-repl
julia> outer(1.0+3.0*e12,e1)
```
returns `( 1.0 e1 )`.

## Example (in 𝔾₄₁ conformal)

```julia-repl
julia> u=2.0*id + 1.5*e1 -3.0*e2 +4.0*e∞ -3.0*e0

julia> v=3.0*ep

julia> outer(u,v)
```
returns `( 6.0 e+ ) + ( 4.5 e1+ ) - ( 9.0 e2+ ) - ( 7.5 e+- )`.

An alternative way to run this function is using the infix
operator (^).

```julia-repl
julia> u=2.0*id + 1.5*e1 -3.0*e2 +4.0*e∞ -3.0*e0

julia> v=3.0*ep

julia> u^v
```
"""
function outer(a::GAAbstractType,b::GAAbstractType)
	if isa(a,Blade) && isa(b,Blade)
		len = length(a.vectors[1])
		if a.grade+b.grade>GIdLayout[1]+GIdLayout[2]
			return Blade(0.0,[zeros(len)],0,true)
		else
			A=zeros(len,a.grade+b.grade)
			for i=1:a.grade
				A[:,i]=a.vectors[i]
			end
			for i=a.grade+1:b.grade
				A[:,i]=b.vectors[i]
			end
			if rank(A)<a.grade+b.grade
				return Blade(0.0,[zeros(len)],0,true)
			else
				return Blade(a.scalar*b.scalar,[a.vectors;b.vectors],a.grade+b.grade,true)
			end
		end
	else
		if !isa(a,MultiVector)
			a=multivector(a)
		end
		if !isa(b,MultiVector)
			b=multivector(b)
		end
		ind1=findnz(a.comp)
		ind2=findnz(b.comp)
		nnv1=nnz(a.comp)
		nnv2=nnz(b.comp)
    	ms=spzeros(GIdLayout[4])
		for i=1:nnv1
			vl=ind1[2][i]
			for j=1:nnv2
				ind=BBouterprodTable[ind1[1][i],ind2[1][j]].index
				if ind>0
					ms[ind]=ms[ind]+vl*ind2[2][j]
				end
				if ind<0
					ms[-ind]=ms[-ind]-vl*ind2[2][j]
				end
			end
    	end
		return MultiVector(ms)
	end
end

"""
```
^(a::GAAbstractType,b::GAAbstractType) (Liga function)
```
In an infix operator to outer product.
This function is used to determine the outer product
of GAAbstractType elements (Basis Blade, MultiVector or Blades).

## Example (in 𝔾₃)

```julia-repl

julia> 1.0+3.0*e12^e1

```
returns `( 1.0 e1 )`.

## Example (in 𝔾₄₁ conformal)

```julia-repl
julia> u=2.0*id + 1.5*e1 -3.0*e2 +4.0*e∞ -3.0*e0

julia> v=3.0*ep

julia> u^v
```
returns `( 6.0 e+ ) + ( 4.5 e1+ ) - ( 9.0 e2+ ) - ( 7.5 e+- )`.

"""
function Base.:^(a::GAAbstractType,b::GAAbstractType)
	return outer(a,b)
end

#############################################################################
# sum and difference operations                                             #
#############################################################################
function Base.:+(a::Number,b::GAAbstractType)
	a=multivector(a)
	if !isa(b,MultiVector)
		b=multivector(b)
	end
	return MultiVector(a.comp+b.comp)
end
function Base.:+(b::GAAbstractType,a::Number)
	a=multivector(a)
	if !isa(b,MultiVector)
		b=multivector(b)
	end
	return MultiVector(a.comp+b.comp)
end
function Base.:+(a::GAAbstractType,b::GAAbstractType)
	if !isa(a,MultiVector)
		a=multivector(a)
	end
	if !isa(b,MultiVector)
		b=multivector(b)
	end
	return MultiVector(a.comp+b.comp)
end
function Base.:-(a::Number,b::GAAbstractType)
	a=multivector(a)
	if !isa(b,MultiVector)
		b=multivector(b)
	end
	return MultiVector(a.comp-b.comp)
end
function Base.:-(b::GAAbstractType,a::Number)
	a=multivector(a)
	if !isa(b,MultiVector)
		b=multivector(b)
	end
	return MultiVector(b.comp-a.comp)
end
function Base.:-(a::GAAbstractType,b::GAAbstractType)
	if !isa(a,MultiVector)
		a=multivector(a)
	end
	if !isa(b,MultiVector)
		b=multivector(b)
	end
	return MultiVector(a.comp-b.comp)
end
#Basis blade signal
function Base.:-(a::BasisBlade)
	if a.index>0
		return BasisBlade(-a.index)
	else
		return BasisBlade(a.index)
	end
end

##############################################################################
# scalar product                                                             #
##############################################################################
"""
```
scalar(a::Number,b::GAAbstractType) (Liga function)
```
This function is used to define the most simple multivector structure.

## Example (in 𝔾₃₁)
```julia-repl
julia> u = scalar(0.13,e1)
```
returns a multivector defined by [0.0,0.13,0.0,...,0.0] (see more about
MultiVectors). Consequently it is a constructor if b is a blade basis. If b is
a MultiVector, it return the scalar product. For example,
```julia-repl
julia> u = scalar(0.13,e2-e1+2.0)
```
returns the multivector ( 0.26 id ) - (0.13 e1 ) + ( 0.13 e2 ).
An alternative way to use this command is by the simbol ``*``. For example,
```julia-repl
julia> u = 0.13*(e2-e1+2.0)
```
"""
function scalar(a::Number,b::GAAbstractType)
	if !isa(b,MultiVector)
		b=multivector(b)
	end
	return MultiVector(a.*(b.comp))
end
function scalar(b::GAAbstractType,a::Number)
	return scalar(a,b)
end
"""
```
scalar(A::Blade,B::Blade) (Liga function)
```
This function returns a scalar product of
two blades
"""
function scalar(A::Blade,B::Blade)
	if A.grade==B.grade
		return inner(A,B)
	else
		return 0.0
	end
end
"""
```
*(a::Number,b::GAAbstractType) (Liga function)
```
This function is used to define the most simple multivector structure.

## Example (in 𝔾₃₁)

```julia-repl
julia> u = scalar(0.13,e1)
```
returns a multivector defined by [0.0,0.13,0.0,...,0.0] (see more about
MultiVectors). Consequently it is a constructor if b is a blade basis. If b is
a MultiVector, it return the scalar product. For example,
```julia-repl
julia> u = scalar(0.13,e2-e1+2.0)
```
returns the multivector ( 0.26 id ) - (0.13 e1 ) + ( 0.13 e2 ).
An alternative way to use this command is by the simbol ``*``. For example,
```julia-repl
julia> u = 0.13*(e2-e1+2.0)
```
"""
function Base.:*(a::Number,b::GAAbstractType)
	return scalar(a,b)
end
function Base.:*(b::GAAbstractType,a::Number)
	return scalar(a,b)
end
function Base.:*(A::Blade,B::Blade)
	return scalar(A,B)
end
##############################################################################
# Involutions                                                                #
##############################################################################

"""
```
reverse(b::BasisBlade) (Liga function)
```
This function calculates the reverse of a Basis Blade.
For example,
```julia-repl
julia> u = reverse(e12)
```
returns the BasisBlade `-e12`.
"""
function reverse(b::BasisBlade)
	k=grade(b)
	if (-1.0)^(k*(k-1)/2.0)>0.0
		return b
	else
		return BasisBlade(-b.index)
	end
end

"""
```
reverse(V::MultiVector) (Liga function)
```
This function calculates the reverse of a MultiVector.
It returns a new MultiVector.

## Example (in 𝔾₃₁)

```julia-repl
julia> u=1.0*id -2.0*e1 -5.0*e∞ +3.0*e∞0
julia> reverse(u)
```
returns the multivector `(1.0 id) - ( 2.0 e1) - ( 5.0 e∞ ) - ( 3.0 e∞0 )`.
"""
function reverse(b::MultiVector)
	ind=findnz(b.comp)
	nnv=nnz(b.comp)
	ms=spzeros(GIdLayout[4])
	for i=1:nnv
		a=reverse(BasisBlade(ind[1][i]))
		if a.index>0
			ms[a.index]=ind[2][i]
		else
			ms[-a.index]=-ind[2][i]
		end
	end
	return MultiVector(ms)
end

"""
```
reverse(B::Blade) (Liga function)
```
This function calculates the reverse of a Blade.
It returns a new Blade.

## Example (in 𝔾₃₁)
```julia-repl
julia> B=Blade([[1.0,0.0,0.0],[0.0,1.0,0.0]],2,true)
julia> reverse(B)
```
returns `-1.0*[1.0, 0.0, 0.0]∧[0.0, 1.0, 0.0]`.
"""
function reverse(A::Blade)
	#tmp=Vector{Float64}(undef,GIdLayout[])
	A.scalar=(-1.0)^((A.grade*(A.grade-1))/2.0)
	return A
end

"""
```
conjugate(b::BasisBlade) (Liga function)
```
This function calculates the conjugate of a Basis Blade.

## Example (in 𝔾₃₁)

```julia-repl
julia> u = conjugate(e12)
```
returns the BasisBlade `-e12`.
"""
function conjugate(b::BasisBlade)
	#pergunta: como fica no conformal?
	r=grademinus(b)
	if (-1.0)^r>0.0
		return reverse(b)
	else
		return -reverse(b)
	end
end

"""
```
conjugate(V::MultiVector) (Liga function)
```
This function calculates the conjugate of a MultiVector
and return a multivector.

## Example (in 𝔾₃₁)

```julia-repl
julia> u = -1.0+3.0*e1-1.0*e12

julia> conjugate(u)
```
returns the MultiVector `- ( 1.0 id ) + ( 3.0 e1 ) + ( 1.0 e12 )`.
"""
function conjugate(b::MultiVector)
	ind=findnz(b.comp)
	nnv=nnz(b.comp)
	ms=spzeros(GIdLayout[4])
	for i=1:nnv
		a=conjugate(BasisBlade(ind[1][i]))
		if a.index>0
			ms[a.index]=ind[2][i]
		else
			ms[-a.index]=-ind[2][i]
		end
	end
	return MultiVector(ms)
end

function conjugatevector(v::Vector{Float64})
	for i=2:length(v)+1
		if grademinus(BasisBlade(i))==1
			v[i-1]=-v[i-1]
		end
	end
	return v
end

"""
```
conjugate(B::Blade) (Liga function)
```
This function calculates the conjugate of a Blade
and return a new Blade.

## Example (in 𝔾₃₂)

```julia-repl
julia> u = blade([1.0,2.0,3.0,4.0,5.0],[-1.0,-1.0,-1.0,-1.0,3.0])

julia> conjugate(u)
```
returns the Blade `[-1.0, -1.0, -1.0, 1.0, -3.0]∧[1.0, 2.0, 3.0, -4.0, -5.0]`.
"""
function conjugate(A::Blade)
	v=Vector{Vector{Float64}}(undef,A.grade)
	for i=1:A.grade
		v[i]=conjugatevector(A.vectors[A.grade+1-i])
	end
	return Blade(v,A.grade,true)
end

###############################################################################
#               Duality                                                       #
###############################################################################
"""
```
dual(b::GAAbstractType) (Liga function)
```
This function calculates the dual of a GAAbstractType
and returns a multivector element.

## Example (in 𝔾₃₁)

```julia-repl
julia> u = dual(e12)
```
returns the multivector ( 1.0 e+- ).

## Example (in 𝔾₃₂)

```julia-repl
julia> b=blade([1.0,0.0,3.0,0.0,5.0],[0.0,-1.0,-1.0,-1.0,0.0])
```

returns the multivector  `- ( 5.0 e123 ) - ( 5.0 e124 ) + ( 5.0 e134 ) - ( 3.0 e125 ) - ( 1.0 e235 ) - ( 3.0 e145 ) - ( 1.0 e245 ) + ( 1.0 e345 )`.
"""
function dual(b::GAAbstractType)
	if !isa(b,MultiVector)
		b=multivector(b)
	end
	#b=multivector(a)
	ind=findnz(b.comp)
	nnv=nnz(b.comp)
	ms=ind[2][1]*bbdual(BasisBlade(ind[1][1]))
	for i=2:nnv
		ms=ms+ind[2][i]*bbdual(BasisBlade(ind[1][i]))
	end
	return ms
end

"""
```
bbdual(b::BasisBlade) (Liga function)
```
This function calculates the dual of a BasisBlade type
and returns a multivector element. It is an auxiliary function.
"""
function bbdual(a::BasisBlade)
	#cj=conjugate(I)
	#if cj.index<0
	#	return -bbgeometric(a,BasisBlade(-cj.index))
	#else
	#	return  bbgeometric(a,cj)
	#end
	return geometric(a,conjugate(I))
end
#function dual(b::MultiVector)
#	ind=findnz(b.comp)
#	nnv=nnz(b.comp)
#	ms=ind[2][1]*bbdual(BasisBlade(ind[1][1]))
#	for i=2:nnv
#		ms=ms+ind[2][i]*bbdual(BasisBlade(ind[1][i]))
#	end
#	return ms
#end
#function dual(A::Blade)
#	return dual(multivector(A))
#end

###############################################################################
#               Magnitude                                                     #
###############################################################################
"""
```
magnitude(A::GAAbstractType) (Liga function)
```
This function calculates the magnitude of a general GAAbstractType
and returns a scalar.

## Example in 𝔾₄₁ (layout(4,1,"GA"))

```julia-repl
julia> magnitude(2.0*id -2.0e3)
```
returns `2.8284271247461903`.
"""
function magnitude(A::GAAbstractType)
	if !isa(A,MultiVector)
		A=multivector(A)
	end
	return norm(A.comp,2)
end

###############################################################################
#              Inverse                                                        #
###############################################################################
"""
```
inverse(b::BasisBlade) (Liga function)
```
This function calculates the inverse of a Basis Blade.
For example, in layout(3,1,"GA")
```julia-repl
julia> u = inverse(e12)
```
returns the BasisBlade -e12. Note that conjugate (for basis blade) is equal to
inverse command.
"""
function inverse(a::BasisBlade)
	return conjugate(a)
end
"""
```
inverse(b::MultiVector) (Liga function)
```
This function calculates the inverse of a multivector element.
For example, in layout(3,1,"GA")
```julia-repl
julia>  u = inverse(e1+e12-(1.5*en))
```
returns the multivector ( 0.5714 e1 ) - ( 0.0952 e- ) - ( 0.5714 e12 ) - ( 0.7619 e2- ).

"""
function inverse(a::MultiVector)
	len=GIdLayout[4]
	A=spzeros(len,len)
	for i=1:len
		A[i,:]=a.comp'*Γ[i]
	end
	b=zeros(len)
	b[1]=1.0
	#display(det(A))
	return MultiVector(sparse(A\b))
end

#need to check!
function inverse(A::Blade)
	B=reverse(A)
	B.vectors[1]=(1.0/geometric(A,B).comp[1])*B.vectors[1]
	return Blade(B.vectors,A.grade,true)
end
###############################################################################
# embedding                                                                   #
###############################################################################
"""
```
embedding(v::Vector{Float64}) (Liga function)
```
This function makes sense just in Conformal or Projective layout.
It was define to embedding a real vector in corresponding enviromment.
For exemplo, let us suppose that layout(4,1,"Conformal") was defined.
So, it is possible embedding a vector (like [1.0,2.0,3.0]) 	in conformal
space as the following.

```julia-repl
julia>  embedding([1.0,2.0,3.0])
```
returns the multivector ( 1.0 e1 ) + ( 2.0 e2 ) + ( 3.0 e3 ) + ( 6.5 e+ ) + ( 7.5 e- ).
If projective space was defined (like layout(4,0,"Projective")) the same vector
becomes ( 1.0 e1 ) + ( 2.0 e2 ) + ( 3.0 e3 ) + ( 1.0 e4 ).
"""
function embedding(x::Vector{Float64})
	if GIdLayout[3]=="GA"
		error("function defined but not enable for this layout.")
	else
		ms=spzeros(GIdLayout[4])
		n=length(x)
		if n!=GIdLayout[1]-1
			error("something wrong with layout dimension and vector dimension")
		else
			if GIdLayout[3]=="Projective"
				ms[2:n+1]=x
				ms[n+2]=1.0
				return MultiVector(ms)
			end
			if GIdLayout[3]=="Conformal"
				ms[2:n+1]=x
				a=0.5*(x'*x)
				ms[n+2]=a-0.5
				ms[n+3]=a +0.5
				return MultiVector(ms)
			end
		end
	end
end
###############################################################################
# convertions                                                                 #
###############################################################################
function Base.convert(::Type{Float64}, x::MultiVector)
	return x.comp[1]
end
function Base.Float64(x::MultiVector)
	return convert(Float64,x)
end

##############################################################################
# copy function                                                              #
# ############################################################################
function Base.copy(A::MultiVector)
	return MultiVector(copy(A.comp))
end
function Base.copy(A::Blade)
	return Blade(copy(A.vectors),copy(A.grade),copy(A.checked))
end

##############################################################################
# equal function                                                             #
# ############################################################################
function Base.isequal(A::MultiVector,B::MultiVector)
	return isequal(A.comp,B.comp)
end

function Base.:(==)(A::MultiVector,B::MultiVector)
	return A.comp==B.comp
end
##############################################################################
# approx function                                                            #
##############################################################################
function Base.:≈(A::MultiVector,B::MultiVector)
	return isapprox(A.comp,B.comp,atol=LigaPrecision)
end


text:
using FastTransforms
using Base.Test

c = collect(linspace(-big(1.0),1,16))
@test norm(fft(c) - fft(map(Float64,c))) < 3Float64(norm(c))*eps()
@test norm(ifft(c) - ifft(map(Float64,c))) < 3Float64(norm(c))*eps()

c = collect(linspace(-big(1.0),1.0,201))
@test norm(ifft(fft(c))-c) < 200norm(c)eps(BigFloat)

p = plan_dct(c)
@test norm(dct(c) - p*c) == 0

pi = plan_idct!(c)
@test norm(pi*dct(c) - c) < 400norm(c)*eps(BigFloat)

@test norm(dct(c)-dct(map(Float64,c)),Inf) < 10eps()

cc = cis(c)
@test norm(dct(cc)-dct(map(Complex{Float64},cc)),Inf) < 10eps()


text:
@using LinearAlgebraicRepresentation as LAR
@using ViewerGL as GL
using SparseArrays
include("../../src/Lar.jl")
LarA = Lar.Arrangement

model = Lar.Model([
    0.0 1.0 0.0 0.5 0.1
    0.0 0.0 1.0 0.1 0.5
])
Lar.addModelCells!(model, 1, SparseArrays.sparse(Array{Int8, 2}([
    -1  1  0  0  0
    -1  0  1  0  0
     0 -1  1  0  0
    -1  0  0  1  0
    -1  0  0  0  1
     0  0  0 -1  1
])))
Lar.addModelCells!(model, 2, SparseArrays.sparse(Array{Int8, 2}([
     1 -1  1  0  0  0
     0  0  0  1 -1  1
])))
Lar.view2DModel(model)
arrangedModel = LarA.planar_arrangement(model)
Lar.view2DModel(arrangedModel, 1)


text:
using GenericLinearAlgebra
using LinearAlgebra
using SpecialFunctions
using Dates
using MAT
using Roots
using FFTW
using PyPlot
using PolynomialRoots
using Printf

function Remez(targetf, parity, degree, xapp, lef = -1, rig = 1, maxiter = 10, sample_size = 20, eps = 1e-12)
    #--------------------------------------------------------------------------------------------------------------
    # Input:
    #     targetf: target funciton
    #     parity: the parity of approximation function, 0 -- even, 1 -- odd, 2 -- no constraint
    #     degree: the degree of approximation polynomial
    #     xapp: initial points to approximate, if provide a single zero, will choose roots of Chebyshev polynomials
    #     lef, rig: left and right endpoints of approximation interval
    #     maxiter: Max iteration
    #     sample_size: sample size used in Remez algorithm
    #     eps: accuracy of approximation
    #
    # Output:
    #     Expansion of approximation polynomial under Chebyshev basis
    #
    # We note that our implementation does not strictly follow the reference. The algorithm may crush when
    # the degree is very large or the problem is ill-conditioned.
    #
    #------------------------------------------------------------------------------------------------------------
    #
    # Reference:
    #     E. W. Cheney.
    #     Introduction to approximation theory
    #
    # Author: X. Meng
    # Version 1.0 .... 02/2020
    #
    #------------------------------------------------------------------------------------------------------------

    if(xapp==0)
        xapp = big.(cos.(collect(range(pi/2,stop=0,length=degree+1))))
    end
    iter = 0
    eps2 = 1e-18
    comp = ones(degree+1,1)
    comp[2:2:end] .= -1
    while(true)
        iter += 1

        # Step1: Find the best approximation on n+1 given points by solving a linear equation

        A = zeros(BigFloat,degree+1,degree+1)
        b = zeros(BigFloat,degree+1,1)
        for i=1:degree
            if(parity==0)
                deg = 2*(i-1)
            elseif(parity==1)
                deg = 2*i-1
            else
                deg = i-1
            end
            for j=1:degree+1
                A[j,i] = chebyshev(xapp[j],deg)
            end
        end
        A[:,end] = comp
        for i=1:degree+1
            b[i] = targetf(xapp[i])
        end
        sol = A\b

        # Step2: Find roots of residual function

        eps3 = findmin([eps2,1e-4*abs(sol[end])])[1]
        xroot = zeros(BigFloat,degree,1)
        for i=1:degree
            rootl = xapp[i]
            rootr = xapp[i+1]
            pm = (-1)^(i)*sign(sol[end])
            froot(x) = chebyshevfunc(x,sol[1:end-1],parity,-targetf(x))
            rootiter = 0
            xroot[i] =  brent(rootl,rootr,froot,eps3,eps3)
            if(xroot[i]==Inf)
                return
            end
        end

        # Step3: In each pair of adjacent roots, find a point x such that the absolute value of
        #        residual function is maximized. In addition, values of residual function alternate
        #        in sign.

        xappnew = copy(xapp)
        maxtot = abs(sol[end])
        r_max = -Inf
        y_max = -Inf

        for i=1:degree+1
            maxapp = abs(sol[end])
            pm = (-1)^(i)*sign(sol[end])
            if(i==1)
                lend = lef
            else
                lend = xroot[i-1]
            end
            if(i==degree+1)
                rend = rig
            else
                rend = xroot[i]
            end
            exh = collect(range(lend,stop=rend,length=sample_size))
            for j=1:length(exh)
                fval = chebyshevfunc(exh[j],sol[1:end-1],parity,-targetf(exh[j]))
                if(abs(fval)>r_max)
                    r_max = abs(fval)
                end
                if(fval*pm>y_max)
                    y_max = fval*pm
                end
                if(fval*pm>maxapp)
                    maxapp = fval*pm
                    maxtot = maximum([maxtot,maxapp])
                    xappnew[i] = exh[j]
                end
            end
        end
        xapp = copy(xappnew)

        if(abs(r_max-y_max)>1e-12)
            @printf("Warning: the interpolation points maybe incorrect\n")
        end
        @printf("The %3d-th itertion: previous error is %5.4e, L_inf approximation error is %5.4e\n",iter,sol[end],maxtot)
        if(maxtot<eps||iter>=maxiter||abs((sol[end]-maxtot)/sol[end])<1e-4)
            return sol[1:end-1]
        end
    end
    end

function chebyshev(x,n) # T_n(x)
    if(abs(x)<=1)
        return cos(n*acos(x))
    elseif(x>1)
        return cosh(n*acosh(x))
    else
        return (-1)^n*cosh(n*acosh(-x))
    end
end

function chebyshevfunc(x,sol,parity,init)
    # Compute the value of a summation of Chebyshev polynomials at x
    y = init
    for i=1:length(sol)
        if(parity==0)
            deg = 2*(i-1)
        elseif(parity==1)
            deg = 2*i-1
        else
            deg = i-1
        end
        y += sol[i]*chebyshev(x,deg)
    end
    return y
end

function brent(a,b,f,tol1,tol2)
    # brent method for finding roots on a given interval [a,b]
    # f(a)f(b) should be less than 0, tol1, tol2 are stopping criteria
    fa = f(a)
    fb = f(b)
    if(f(a)*f(b)>0)
        println("Error, f(a)f(b)>0")
        return Inf
    end
    if(abs(f(a))<abs(f(b)))
        tmp = a
        tmp2 = fa
        a = b
        fa = fb
        b = tmp
        fb = tmp2
    end
    c = a
    fc = fa
    s = b
    d = 1e-10
    iter = 0
    mflag = true
    while(true)
        iter += 1
        if(iter>1000)
            println("Brent method: reaches max iteration.")
            return Inf
        end
        if(abs(b-a)<tol2||abs(f(s))<tol1)
            return s
        end
        if(fa!=fc&&fb!=fc)
            s = a*fb*fc/((fa-fb)*(fa-fc))+b*fa*fc/((fb-fa)*(fb-fc))+c*fa*fb/((fc-fa)*(fc-fb))
        else
            s = b-fb*(b-a)/(fb-fa)
        end
        if((s>=b)||(s<=(3*a+b)/4)||((abs(s-b)*2>abs(b-c))&&mflag)||((abs(s-b)*2>abs(d-c))&&(!mflag))||(tol1>abs(b-c)&&mflag)||((tol1>abs(d-c))&&(!mflag)))
            s = (a+b)/2
            mflag = true
        else
            mflag = false
        end
        fs = f(s)
        d = c
        c = b
        fc = fb
        if(fa*fs<0)
            b = s
            fb = fs
        else
            a = s
            fa = fs
        end
        if(abs(f(a))<abs(f(b)))
            tmp = a
            tmp2 = fa
            a = b
            fa = fb
            b = tmp
            fb = tmp2
        end
    end
end

d = Remez(x->1/x, 1, 10, 0, 1, 10)

using DelimitedFiles

# open file in append mode
file = open("example_file.txt", "w")

# writing to a file using write() method
writedlm(file, d)

# We need to close the file in order to write the content from the disk to file
close(file)

max_stars_repo_name: Aaron-Gregory/RE-QSVT-HHL

max_issues_repo_name: Aaron-Gregory/RE-QSVT-HHL

max_forks_repo_name: Aaron-Gregory/RE-QSVT-HHL


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

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

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

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

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

check_likelihoodtype(ltype::Symbol) = Symbol(ltype) == :Regression || Symbol(ltype) == :Classification || Symbol(ltype) == :Event

function treat_likelihood(likelihood::Expr)

end

function treat_params(params)

end

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

const AGP=AugmentedGaussianProcesses

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

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

end

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

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

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

        function AGP.pdf(l::$(lname),y::Real,f::Real)
            C(l)*exp(g(l,y)*f)*φ(l,α(l,y)-β(l,y)*f+γ(l,y)*f^2)
        end

        function AGP.logpdf(l::$(lname),y::Real,f::Real)
            log(C(l))+g(l,y)*f+log(φ(l,α(l,y)-β(l,y)*f+γ(l,y)*f^2))
        end

        function C(l::$(lname){T}) where {T}
            C()
        end

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

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

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

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

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

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

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

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

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

        function AGP.compute_proba(l::$(lname){T},y::AbstractVector,μ::AbstractVector{T},σ²::AbstractVector{T}) where {T<:Real}
            if typeof(l) <: RegressionLikelihood
                return μ,max.(σ²,zero(T)).+var(l)
            elseif typeof(l) <: ClassificationLikelihood
                N = length(μ)
                pred = zeros(T,N)
                sig_pred = zeros(T,N)
                for i in 1:N
                    x = AGP.pred_nodes.*sqrt(max(σ²[i],zero(T))).+μ[i]
                    pred[i] = dot(AGP.pred_weights,AGP.pdf.(l,1,x))
                    sig_pred[i] = dot(AGP.pred_weights,AGP.pdf.(l,1,x).^2) - pred[i]^2
                end
                return pred,sig_pred
            else
                @error "Prediction not implemented yet, please file an issue"
            end
        end

        ### Local Updates Section ###

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

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

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

        ### Natural Gradient Section ###

        @inline AGP.∇E_μ(l::$(lname),::AOptimizer,y::AbstractVector) where {T} = (g(l,y)+l.θ.*β(l,y),)
        @inline AGP.∇E_Σ(l::$(lname),::AOptimizer,y::AbstractVector) where {T} = (l.θ.*γ(l,y),)

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

        function AugmentedKL(model::AbstractGP{T,<:$(lname),<:AnalyticVI}) where {T}
            -dot(l.c²,l.θ)-sum(log,φ.(l,c²))
        end

        ### Gradient Section ###

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

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

    end
end
end


text:
"""
    cnse_boundary_state!(::NumericalFluxFirstOrder, ::Impenetrable{FreeSlip}, ::CNSE3D)
apply free slip boundary condition for velocity
sets reflective ghost point
"""
@inline function cnse_boundary_state!(
    ::NumericalFluxFirstOrder,
    ::Impenetrable{FreeSlip},
    ::CNSE3D,
    ::TurbulenceClosure,
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
    args...,
)
    state⁺.ρ = state⁻.ρ

    ρu⁻ = state⁻.ρu
    state⁺.ρu = ρu⁻ - 2 * n⁻ ⋅ ρu⁻ .* SVector(n⁻)

    return nothing
end

"""
    cnse_boundary_state!(::NumericalFluxGradient, ::Impenetrable{FreeSlip}, ::CNSE3D)
apply free slip boundary condition for velocity
sets non-reflective ghost point
"""
function cnse_boundary_state!(
    ::NumericalFluxGradient,
    ::Impenetrable{FreeSlip},
    ::CNSE3D,
    ::ConstantViscosity,
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
    args...,
)
    state⁺.ρ = state⁻.ρ

    ρu⁻ = state⁻.ρu
    state⁺.ρu = ρu⁻ - n⁻ ⋅ ρu⁻ .* SVector(n⁻)

    return nothing
end

"""
    shallow_normal_boundary_flux_second_order!(::NumericalFluxSecondOrder, ::Impenetrable{FreeSlip}, ::CNSE3D)
apply free slip boundary condition for velocity
apply zero numerical flux in the normal direction
"""
function cnse_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{FreeSlip},
    ::CNSE3D,
    ::ConstantViscosity,
    state⁺,
    gradflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    aux⁻,
    t,
    args...,
)
    state⁺.ρu = state⁻.ρu
    gradflux⁺.ν∇u = n⁻ * (@SVector [-0, -0, -0])'

    return nothing
end

"""
    cnse_boundary_state!(::NumericalFluxFirstOrder, ::Impenetrable{NoSlip}, ::CNSE3D)
apply no slip boundary condition for velocity
sets reflective ghost point
"""
@inline function cnse_boundary_state!(
    ::NumericalFluxFirstOrder,
    ::Impenetrable{NoSlip},
    ::CNSE3D,
    ::TurbulenceClosure,
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
    args...,
)
    state⁺.ρ = state⁻.ρ
    state⁺.ρu = -state⁻.ρu

    return nothing
end

"""
    cnse_boundary_state!(::NumericalFluxGradient, ::Impenetrable{NoSlip}, ::CNSE3D)
apply no slip boundary condition for velocity
set numerical flux to zero for U
"""
@inline function cnse_boundary_state!(
    ::NumericalFluxGradient,
    ::Impenetrable{NoSlip},
    ::CNSE3D,
    ::ConstantViscosity,
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
    args...,
)
    FT = eltype(state⁺)
    state⁺.ρu = @SVector zeros(FT, 3)

    return nothing
end

"""
    cnse_boundary_state!(::NumericalFluxSecondOrder, ::Impenetrable{NoSlip}, ::CNSE3D)
apply no slip boundary condition for velocity
sets ghost point to have no numerical flux on the boundary for U
"""
@inline function cnse_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{NoSlip},
    ::CNSE3D,
    ::ConstantViscosity,
    state⁺,
    gradflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    aux⁻,
    t,
    args...,
)
    state⁺.ρu = -state⁻.ρu
    gradflux⁺.ν∇u = gradflux⁻.ν∇u

    return nothing
end

"""
    cnse_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::Penetrable{FreeSlip}, ::CNSE3D)
no mass boundary condition for penetrable
"""
cnse_boundary_state!(
    ::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    ::Penetrable{FreeSlip},
    ::CNSE3D,
    ::ConstantViscosity,
    _...,
) = nothing

"""
    cnse_boundary_state!(::NumericalFluxSecondOrder, ::Penetrable{FreeSlip}, ::CNSE3D)
apply free slip boundary condition for velocity
apply zero numerical flux in the normal direction
"""
function cnse_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Penetrable{FreeSlip},
    ::CNSE3D,
    ::ConstantViscosity,
    state⁺,
    gradflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    aux⁻,
    t,
    args...,
)
    state⁺.ρu = state⁻.ρu
    gradflux⁺.ν∇u = n⁻ * (@SVector [-0, -0, -0])'

    return nothing
end

"""
    cnse_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::Impenetrable{MomentumFlux}, ::HBModel)
apply kinematic stress boundary condition for velocity
applies free slip conditions for first-order and gradient fluxes
"""
function cnse_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    ::Impenetrable{<:MomentumFlux},
    model::CNSE3D,
    turb::TurbulenceClosure,
    args...,
)
    return cnse_boundary_state!(
        nf,
        Impenetrable(FreeSlip()),
        model,
        turb,
        args...,
    )
end

"""
    cnse_boundary_state!(::NumericalFluxSecondOrder, ::Impenetrable{MomentumFlux}, ::HBModel)
apply kinematic stress boundary condition for velocity
sets ghost point to have specified flux on the boundary for ν∇u
"""
@inline function cnse_boundary_state!(
    ::NumericalFluxSecondOrder,
    bc::Impenetrable{<:MomentumFlux},
    model::CNSE3D,
    ::ConstantViscosity,
    state⁺,
    gradflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    aux⁻,
    t,
    args...,
)
    state⁺.ρu = state⁻.ρu
    gradflux⁺.ν∇u = n⁻ * bc.drag(state⁻, aux⁻, t)'

    return nothing
end

"""
    cnse_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::Penetrable{MomentumFlux}, ::HBModel)
apply kinematic stress boundary condition for velocity
applies free slip conditions for first-order and gradient fluxes
"""
function cnse_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    ::Penetrable{<:MomentumFlux},
    model::CNSE3D,
    turb::TurbulenceClosure,
    args...,
)
    return cnse_boundary_state!(
        nf,
        Penetrable(FreeSlip()),
        model,
        turb,
        args...,
    )
end

"""
    cnse_boundary_state!(::NumericalFluxSecondOrder, ::Penetrable{MomentumFlux}, ::HBModel)
apply kinematic stress boundary condition for velocity
sets ghost point to have specified flux on the boundary for ν∇u
"""
@inline function cnse_boundary_state!(
    ::NumericalFluxSecondOrder,
    bc::Penetrable{<:MomentumFlux},
    shallow::CNSE3D,
    ::ConstantViscosity,
    state⁺,
    gradflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    aux⁻,
    t,
    args...,
)
    state⁺.ρu = state⁻.ρu
    gradflux⁺.ν∇u = n⁻ * bc.drag(state⁻, aux⁻, t)'

    return nothing
end


text:
# maintenance.jl

println("Updating all packages\n")
@time Pkg.update()

println("Done!\n")

max_stars_repo_name: anthony-wang/julia-scripts

max_issues_repo_name: anthony-wang/julia-scripts

max_forks_repo_name: anthony-wang/julia-scripts