text:
abstract type AbstractPoissonProblem end

mutable struct PoissonProblem{C,InvLaplacian,Diag,Vec,Yt,Metric} <: AbstractPoissonProblem
    k::C
    T⁻¹::InvLaplacian # Factorized Laplacian
    r⁻¹::Diag
    rhs::Vec # Right-hand side of Poisson problem
    Y::Yt
    rᵏ::Vec
    rᵏ⁺¹::Vec
    rₘₐₓ⁻²ᵏ⁺¹::C
    S::Metric
    ∫ρ::C # Integrated mutual density
end

Base.eltype(::PoissonProblem{C}) where C = C

function Base.show(io::IO, p::PoissonProblem{C}) where C
    m,n = size(p.T⁻¹)
    write(io, "$(m)×$(n) PoissonProblem $C, multipole = $(p.k)")
end

isrealtype(::Type{T}) where {T<:Real} = true
isrealtype(::Type{T}) where {T<:Complex} = false

"""
    get_double_laplacian(R, k)

Return Laplacian (multiplied by `2`) for partial wave `k` of the
basis `R`, `-∂ᵣ² + k(k+1)/r²`.
"""
function get_double_laplacian(R,k,::Type{T}) where T
    D = Derivative(axes(R,1))
    Tᵏ = R'D'D*R
    r = axes(R,1)
    if Tᵏ isa BlockSkylineMatrix
        Tᵏ = Matrix(Tᵏ)
    end
    V = R'*QuasiDiagonal(k*(k+1)./r.^2)*R
    Tᵏ -= V
    isrealtype(T) ? Tᵏ : complex(Tᵏ)
end

struct PoissonCache{M,F,Diag,V₁,V₂}
    T::M
    T⁻¹::F
    r⁻¹::Diag
    rᵏ::V₁
    rᵏ⁺¹::V₂
end

function get_or_create_poisson_cache!(poisson_cache, k, R, ::Type{C}) where C
    k in keys(poisson_cache) && return poisson_cache[k]
    T = get_double_laplacian(R,k,C)
    T⁻¹ = factorize(T isa BlockSkylineMatrix ? Matrix(T) : T)
    r = axes(R,1)

    f = isrealtype(C) ? identity : complex
    r⁻¹ = f(R'*QuasiDiagonal(1 ./ r)*R)
    rᵏ   = f(R \ r.^k)
    rᵏ⁺¹ = f(R \ r.^(k+1))

    pc = PoissonCache(T, T⁻¹, r⁻¹, rᵏ, rᵏ⁺¹)
    poisson_cache[k] = pc
    pc
end

"""
    PoissonProblem(R, k[, ::Type{C}=eltype(R); poisson_cache])

Create the Poisson problem of order `k` for densities expanded over
`R`.
"""
function PoissonProblem(R, k, ::Type{C}=eltype(R);
                        Y=zeros(C, size(R,2)),
                        poisson_cache::Dict{Int,<:PoissonCache} = Dict{Int,PoissonCache}()) where C
    pc = get_or_create_poisson_cache!(poisson_cache, k, R, C)

    r = axes(R,1)
    rₘₐₓ⁻²ᵏ⁺¹ = inv(rightendpoint(r.domain)^(2k+1))

    PoissonProblem(C(k), pc.T⁻¹, pc.r⁻¹, similar(Y), Y,
                   pc.rᵏ, pc.rᵏ⁺¹, C(rₘₐₓ⁻²ᵏ⁺¹), R'R, zero(C))
end

#=
This method computes the integral \(Y^k(n\ell,n'\ell';r)/r\) where
\[\begin{aligned}
Y^k(n\ell,n'\ell';r)&\equiv
r\int_0^\infty
U^k(r,s)
P^*(n\ell;s)
P(n'\ell';s)
\mathrm{d}s\\
&=
\int_0^r
\left(\frac{s}{r}\right)^k
P^*(n\ell;s)
P(n'\ell';s)
\mathrm{d}s+
\int_r^\infty
\left(\frac{r}{s}\right)^{k+1}
P^*(n\ell;s)
P(n'\ell';s)
\mathrm{d}s
\end{aligned}\]
through the solution of Poisson's problem
\[\left[\frac{\mathrm{d}^2}{\mathrm{d}r^2} -
 \frac{k(k+1)}{r^2}\right]
Y^k(n\ell,n'\ell';r)
=-\frac{2k+1}{r}
P^*(n\ell;r)P(n'\ell';r)\]
subject to the boundary conditions
\[\begin{cases}
Y^k(n\ell,n'\ell';0) &= 0,\\
\frac{\mathrm{d}}{\mathrm{d}r}
Y^k(n\ell,n'\ell';r) &=
-\frac{k}{r}
Y^k(n\ell,n'\ell';r),\quad
r \to \infty.
\end{cases}\]

References:

- Fischer, C. F., & Guo, W. (1990). Spline Algorithms for the
  Hartree-Fock Equation for the Helium Ground State. Journal of
  Computational Physics, 90(2),
  486–496. http://dx.doi.org/10.1016/0021-9991(90)90176-2

- McCurdy, C. W., Baertschy, M., & Rescigno, T. N. (2004). Solving the
  Three-Body Coulomb Breakup Problem Using Exterior Complex
  Scaling. Journal of Physics B: Atomic, Molecular and Optical
  Physics, 37(17),
  137–187. http://dx.doi.org/10.1088/0953-4075/37/17/r01

=#

function solve!(poisson::PoissonProblem, ρ::AbstractVector)
    mul!(poisson.rhs, poisson.r⁻¹, ρ, -(2poisson.k+1), false)
    ldiv!(poisson.Y, poisson.T⁻¹, poisson.rhs)

    # Add in homogeneous contribution
    poisson.∫ρ = dot(poisson.rᵏ, poisson.S, ρ)
    s = poisson.∫ρ*poisson.rₘₐₓ⁻²ᵏ⁺¹
    poisson.Y .+= s .* poisson.rᵏ⁺¹
end

solve!(poisson::PoissonProblem, ρ::Density) =
    solve!(poisson, ρ.ρ)

mutable struct AsymptoticPoissonProblem{PP<:PoissonProblem, B,
                                        I<:AbstractRange,
                                        VO₁, VO₂, VO₃} <: AbstractPoissonProblem
    pp::PP # Poisson problem of the inner region
    R̃::B # Basis of the inner region
    inner::I # Range of inner region
    Y::VO₁ # Solution
    Ytail::VO₂ # View of the asymptotic part of Y′
    Ỹ::VO₃ # Asymptotic part of the solution
end

"""
    AsymptoticPoissonProblem(R, k, u, v, R̃[; w′=similar(U)])

Create the Poisson problem of order `k` for the mutual density `u†(r)
.* v(r)`; the Poisson problem is solved numerically within the domain
of `R̃` and an asymptotic solution is used outside. For this to be
valid, `u` or `v` have to vanish before the end of `R̃`. `w′` is a
`MulQuasiVector` of the same kind as `u` and `v`, and may optionally
be provided as a pre-allocated vector, in case it is e.g. the diagonal
of a potential matrix.
"""
function AsymptoticPoissonProblem(R, k::Int,
                                  R̃::AbstractQuasiMatrix,
                                  ::Type{C}=eltype(R);
                                  kwargs...) where C
    r = axes(R, 1)
    # It is assumed that grid spacings &c agree.
    axes(R̃,1).domain ⊆ r.domain &&
        axes(R̃,2) ⊆ axes(R,2) ||
        throw(ArgumentError("$(R̃) not a subset of $(R)"))
    inner = 1:size(R̃,2)
    tail = inner[end]+1:size(R,2)

    Y = zeros(C, size(R,2))
    pp = PoissonProblem(R̃, k, C; Y=view(Y, inner), kwargs...)

    Ỹ = (R \ r .^ (-k))[tail]

    AsymptoticPoissonProblem(pp, R̃, inner, Y, view(Y, tail), Ỹ)
end

Base.eltype(poisson::AsymptoticPoissonProblem) = eltype(poisson.pp)

function solve!(poisson::AsymptoticPoissonProblem, ρ::Density)
    solve!(poisson.pp, view(ρ.ρ, poisson.inner))

    # This is useful for the case we the asymptotic tail is not
    # actually needed, e.g. when applied to an orbital of compact
    # support.
    if !isempty(poisson.Ỹ)
        # Copy over asymptotic solution and weight it by the charge
        # density within the inner region.
        poisson.Ytail .= poisson.pp.∫ρ .* poisson.Ỹ
    end

    poisson.Y
end

export AbstractPoissonProblem, PoissonProblem, solve!, AsymptoticPoissonProblem


text:
"""
Core genie configuration / settings functionality.
"""
module Configuration

using Genie, YAML

export is_dev, is_prod, is_test, env, cache_enabled, Settings, DEV, PROD, TEST, IN_REPL
export LOG_LEVEL_VERBOSITY_VERBOSE, LOG_LEVEL_VERBOSITY_MINIMAL

# app environments
const DEV   = "dev"
const PROD  = "prod"
const TEST  = "test"

# log levels
const LOG_LEVEL_VERBOSITY_VERBOSE = :verbose
const LOG_LEVEL_VERBOSITY_MINIMAL = :minimal

# defaults
const IN_REPL = false
const GENIE_VERSION = v"0.7"

"""
    is_dev()  :: Bool
    is_prod() :: Bool
    is_test() :: Bool

Set of utility functions that return whether or not the current environment is development, production or testing.

# Examples
```julia
julia> Configuration.is_dev()
true

julia> Configuration.is_prod()
false
```
"""
is_dev()  :: Bool = (Genie.config.app_env == DEV)
is_prod() :: Bool = (Genie.config.app_env == PROD)
is_test() :: Bool = (Genie.config.app_env == TEST)

"""
    env() :: String

Returns the current Genie environment.

# Examples
```julia
julia> Configuration.env()
"dev"
```
"""
env() :: String = Genie.config.app_env

"""
    cache_enabled() :: Bool

Indicates whether or not the app has caching enabled (`cache_duration > 0`).
"""
cache_enabled() :: Bool = (Genie.config.cache_duration > 0)

"""
    mutable struct Settings

App configuration - sets up the app's defaults. Individual options are overwritten in the corresponding environment file.
"""
mutable struct Settings
  server_port::Int
  server_host::String
  server_workers_count::Int
  server_document_root::String
  server_handle_static_files::Bool
  server_signature::String

  app_env::String

  cors_headers::Dict{String,String}
  cors_allowed_origins::Vector{String}

  suppress_output::Bool
  output_length::Int

  tasks_folder::String
  test_folder::String

  log_folder::String

  cache_folder::String
  cache_adapter::Symbol
  cache_duration::Int
  cache_table::String

  log_router::Bool
  log_level::Symbol
  log_verbosity::Symbol
  log_formatted::Bool
  log_cache::Bool
  log_views::Bool

  assets_path::String
  assets_serve::Bool
  assets_fingerprinted::Bool

  tests_force_test_env::Bool

  session_auto_start::Bool
  session_key_name::String
  session_storage::Symbol
  session_folder::String
  session_table::String

  inflector_irregulars::Vector{Tuple{String,String}}

  html_template_engine::Symbol
  json_template_engine::Symbol

  flax_compile_templates::Bool

  lookup_ip::Bool

  run_as_server::Bool

  websocket_server::Bool

  renderer_default_layout_file::Symbol

  Settings(;
            server_port                 = 8000, # default port for binding the web server
            server_host                 = "0.0.0.0",
            server_workers_count        = 1,
            server_document_root        = "public",
            server_handle_static_files  = true,
            server_signature            = "Genie/$GENIE_VERSION/Julia/$VERSION",

            app_env       = ENV["GENIE_ENV"],

            cors_headers  = Dict{String,String}(
              "Access-Control-Allow-Origin"       => "", # ex: "*" or "http://mozilla.org"
              "Access-Control-Expose-Headers"     => "", # ex: "X-My-Custom-Header, X-Another-Custom-Header"
              "Access-Control-Max-Age"            => "86400", # 24 hours
              "Access-Control-Allow-Credentials"  => "", # "true" or "false"
              "Access-Control-Allow-Methods"      => "", # ex: "POST, GET"
              "Access-Control-Allow-Headers"      => "", # ex: "X-PINGOTHER, Content-Type"
            ),
            cors_allowed_origins = String[],

            suppress_output = false,
            output_length   = 10_000, # where to truncate strings in console

            task_folder       = joinpath("task"),
            test_folder       = joinpath("test"),

            log_folder        = joinpath("log"),

            cache_folder      = joinpath("cache"),
            cache_adapter     = :FileCacheAdapter,
            cache_duration    = 0,
            cache_table       = "storage_caches",

            log_router    = false,
            log_level     = :debug,
            log_verbosity = LOG_LEVEL_VERBOSITY_VERBOSE,
            log_formatted = true,
            log_cache     = true,
            log_views     = true,

            assets_path           = "/",
            assets_serve          =  true,
            assets_fingerprinted  = false,

            tests_force_test_env = true,

            session_auto_start  = true,
            session_key_name    = "__GENIESID",
            session_storage     = :File,
            session_folder      = joinpath("session"),
            session_table       = "storage_sessions",

            inflector_irregulars = Tuple{String,String}[],

            html_template_engine = :Flax,
            json_template_engine = :Flax,

            flax_compile_templates = false,

            lookup_ip = true,

            run_as_server = false,

            websocket_server = false,

            renderer_default_layout_file = :app,
        ) =
              new(
                  server_port, server_host,
                  server_workers_count, server_document_root, server_handle_static_files, server_signature,
                  app_env,
                  cors_headers, cors_allowed_origins,
                  suppress_output, output_length,
                  task_folder, test_folder,
                  log_folder,
                  cache_folder, cache_adapter, cache_duration, cache_table,
                  log_router,
                  log_level, log_verbosity, log_formatted, log_cache, log_views,
                  assets_path, assets_serve, assets_fingerprinted,
                  tests_force_test_env,
                  session_auto_start, session_key_name, session_storage, session_folder, session_table,
                  inflector_irregulars,
                  html_template_engine, json_template_engine,
                  flax_compile_templates,
                  lookup_ip,
                  run_as_server,
                  websocket_server,
                  renderer_default_layout_file
                )
end

end

ext: jl


text:
function NeQuickModel(arg0::Vector{jdouble})
    return NeQuickModel((Vector{jdouble},), arg0)
end

function NeQuickModel(arg0::Vector{jdouble}, arg1::TimeScale)
    return NeQuickModel((Vector{jdouble}, TimeScale), arg0, arg1)
end

function get_parameters_drivers(obj::NeQuickModel)
    return jcall(obj, "getParametersDrivers", List, ())
end

function path_delay(obj::NeQuickModel, arg0::FieldSpacecraftState, arg1::TopocentricFrame, arg2::jdouble, arg3::Vector{RealFieldElement})
    return jcall(obj, "pathDelay", RealFieldElement, (FieldSpacecraftState, TopocentricFrame, jdouble, Vector{RealFieldElement}), arg0, arg1, arg2, arg3)
end

function path_delay(obj::NeQuickModel, arg0::SpacecraftState, arg1::TopocentricFrame, arg2::jdouble, arg3::Vector{jdouble})
    return jcall(obj, "pathDelay", jdouble, (SpacecraftState, TopocentricFrame, jdouble, Vector{jdouble}), arg0, arg1, arg2, arg3)
end

function stec(obj::NeQuickModel, arg0::AbsoluteDate, arg1::GeodeticPoint, arg2::GeodeticPoint)
    return jcall(obj, "stec", jdouble, (AbsoluteDate, GeodeticPoint, GeodeticPoint), arg0, arg1, arg2)
end

function stec(obj::NeQuickModel, arg0::FieldAbsoluteDate, arg1::FieldGeodeticPoint, arg2::FieldGeodeticPoint)
    return jcall(obj, "stec", RealFieldElement, (FieldAbsoluteDate, FieldGeodeticPoint, FieldGeodeticPoint), arg0, arg1, arg2)
end


text:
using AcFDTD
using Base.Test

schm = IISO()
show(schm)
schm = SLF()
show(schm)
schm = IWB()
show(schm)

env = FDTDEnv(0.1,IISO())
show(env)
env = FDTDEnv(SLF(),1e4)
show(env)
geo = CuboidRoom(1.,2.,3.,50*ones(6), env)
show(geo)
geo = CuboidRoom(1.,2.,3.,50*ones(6), FDTDEnv(0.1,IISO()))
show(geo)
#geo = LShapedRoom(1.,2.,3.,0.5,1.5,50*ones(6), FDTDEnv(0.1,IISO()))
#show(geo)
#geo = LShapedRoom(10,20,14,5,18,50*ones(6), FDTDEnv(0.1,IISO()))
#show(geo)

size(geo)

max_stars_repo_name: nantonel/JuFdtd

max_issues_repo_name: nantonel/JuFdtd

max_forks_repo_name: nantonel/JuFdtd


text:
#This will capture broadcasts. Inelegant but reasonable.
@capture(ex, f_(a_,b_)) && startswith(string(f), '.')

#This will capture indexing:
@capture(ex, A_[ijk__])

#This will capture indexed assignment
@capture(ex, A_[ijk__] = body_)

# Correct formatting/help from Michael Abbott on Slack


text:
#normalization factor
abstract type YLMNorm{T} end

struct Schmidt{T} <: YLMNorm{T}; end
struct Laplace{T} <: YLMNorm{T}; end
struct Nonorm{T} <: YLMNorm{T}; end
struct Full{T} <: YLMNorm{T}; end

function ylmcoeff(N::Schmidt{T},l::Int64, m::Int64) where T

	k = one(T)
	for i in (l-m+1):(l+m)
	  k *= i
	end
	return sqrt(1/k)
end

function ylmcoeff(N::Laplace{T}, l::Int64, m::Int64) where T

  k = one(T)
  for i in (l-m+1):(l+m)
    k *= i
  end

  return sqrt((2*l+1) / (4*T(pi)*k))
end

function ylmcoeff(N::Full{T}, l::Int64, m::Int64) where T
	fac = m==0 ? one(T) : sqrt(2*one(T))
	return ylmcoeff(Laplace{T}(),l,m)*sqrt(4T(pi))*fac
end

function ylmcoeff(N::Nonorm{T}, l::Int64, m::Int64) where T
 return one(T)
end

function cossinpoly(m::Int64, x::Variable, y::Variable)

  sum = zero(x*y)
  for j in 0:div(m,2)
	  k = 2j
    sum += ((-1)^j)*binomial(m, k)*y^k*x^(m-k)
  end
  return sum
end

function sinsinpoly(m::Int64, x::Variable, y::Variable)

  sum = zero(x*y)
  for j in 0:div((m-1),2)
	  k = 2j+1
    sum += ((-1)^j)*binomial(m, k)*y^k*x^(m-k)
  end
  return sum
end

# #legendre polynomial:
# P(l::BigInt,x) = differentiate((x^2 - 1)^l,x,l)
#
# #associated legendre polynomial (without sin(θ)^m)
# P(l::BigInt,m::BigInt,x) = 1//(big(2)^l*factorial(l))*differentiate(P(l,x),x,m)

"""
    ylm(l::Int64, m::Int64, x::Variable, y::Variable, z::Variable)
*Description:*  Calculation of the spherical harmonic for a given order (l,m) in Cartesian coordinates\\

*Input:*  `l`       - Degree of the spherical harmonic\\
          `m`       - Order of the spherical harmonic\\
          `x, y, z` - Cartesian coordinates\\

*Output:*  Spherical harmonic polynomial
"""
function ylm(l::Int64, m::Int64, x::Variable, y::Variable, z::Variable;
			 norm::YLMNorm{T}=Laplace{Float64}(),real=true) where T

	if abs(m) > l
	throw(DomainError(m,"-l <= m <= l expected, but m = $m and l = $l."))
	end

	p = (z^2 - 1)^l

	for i = 1:l+abs(m)
	  c = i <= l ? 1/(2one(T)*i) : one(T)
	  p = c*differentiate(p, z)
	end

	p *= (-1)^m

	if real
		if m > 0
			out = ylmcoeff(norm, l, m)*cossinpoly(m,x,y)*p
			if norm != Nonorm{T}()
				out *= sqrt(2one(T))
			end
			return out
		elseif m < 0
			out = ylmcoeff(norm, l, abs(m))*sinsinpoly(abs(m),x,y)*p
			if norm != Nonorm{T}()
				out *= sqrt(2one(T))
			end
			return out
		else
			return ylmcoeff(norm, l, 0)*p
		end
	else
		return ylmcoeff(norm, l, m)*(cossinpoly(m,x,y)+sign(m)*im*sinsinpoly(m,x,y))*p
	end
end

# multiplying r^l*ylm(x,y,z)
function rlylm(l::Int64, m::Int, x::Variable, y::Variable, z::Variable;
			   norm::YLMNorm{T}=Laplace{Float64}(),real=true) where T
	p = ylm(l,m,x,y,z;norm=norm,real=real)
	tout = []
	# Zerlegung des Polynoms in Terme:
	for t in terms(p)
		deg = degree(monomial(t)) # Gibt den gesamten Grad des Monoms an
		degR = l-deg # durch das Kürzen ergibt sich ein Grad von l-deg fuer r
		push!(tout,(x^2+y^2+z^2)^div(degR,2)*t) # r² wird durch x²+y²+z² ersetzt
	end

	return polynomial(tout)
end

# solid harmonics
function rlm(l::Int64, m::Int64, x::Variable, y::Variable, z::Variable;
			 norm::YLMNorm{T}=Laplace{Float64}(),real=true) where T
	rlm = rlylm(l,m,x,y,z; norm=norm,real=real)
	if norm != Nonorm{T}()
		rlm = sqrt(4one(T)*T(pi)/(2*l+1))*rlm
	end
	return rlm
end


text:
include("mdp_tests.jl")
include("task_generator_tests.jl")


text:
##############
# Issue type #
##############

type Issue <: GitHubType
    id::Nullable{Int}
    number::Nullable{Int}
    comments::Nullable{Int}
    title::Nullable{GitHubString}
    state::Nullable{GitHubString}
    body::Nullable{GitHubString}
    user::Nullable{Owner}
    assignee::Nullable{Owner}
    closed_by::Nullable{Owner}
    created_at::Nullable{Dates.DateTime}
    updated_at::Nullable{Dates.DateTime}
    closed_at::Nullable{Dates.DateTime}
    labels::Nullable{Vector{Dict}}
    milestone::Nullable{Dict}
    pull_request::Nullable{PullRequest}
    url::Nullable{HttpCommon.URI}
    html_url::Nullable{HttpCommon.URI}
    labels_url::Nullable{HttpCommon.URI}
    comments_url::Nullable{HttpCommon.URI}
    events_url::Nullable{HttpCommon.URI}
    locked::Nullable{Bool}
end

Issue(data::Dict) = json2github(Issue, data)
Issue(number::Real) = Issue(Dict("number" => number))

namefield(issue::Issue) = issue.number

###############
# API Methods #
###############

function issue(repo, issue_obj; options...)
    result = gh_get_json("/repos/$(name(repo))/issues/$(name(issue_obj))"; options...)
    return Issue(result)
end

function issues(repo; options...)
    results, page_data = gh_get_paged_json("/repos/$(name(repo))/issues"; options...)
    return map(Issue, results), page_data
end

function create_issue(repo; options...)
    result = gh_post_json("/repos/$(name(repo))/issues"; options...)
    return Issue(result)
end

function edit_issue(repo, issue; options...)
    result = gh_patch_json("/repos/$(name(repo))/issues/$(name(issue))"; options...)
    return Issue(result)
end


text:
# Bio.Align
# =========
#
# Module for sequence alignments.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module Align

export
    Operation,
    AlignmentAnchor,
    Alignment,
    AlignedSequence,
    seq2ref,
    ref2seq,
    ismatchop,
    isinsertop,
    isdeleteop,
    cigar,
    OP_MATCH,
    OP_INSERT,
    OP_DELETE,
    OP_SKIP,
    OP_SOFT_CLIP,
    OP_HARD_CLIP,
    OP_PAD,
    OP_SEQ_MATCH,
    OP_SEQ_MISMATCH,
    OP_BACK,
    OP_START,
    # alignment types
    GlobalAlignment,
    SemiGlobalAlignment,
    OverlapAlignment,
    LocalAlignment,
    EditDistance,
    HammingDistance,
    LevenshteinDistance,
    # substitution matrices
    AbstractSubstitutionMatrix,
    SubstitutionMatrix,
    DichotomousSubstitutionMatrix,
    EDNAFULL,
    PAM30,
    PAM70,
    PAM250,
    BLOSUM45,
    BLOSUM50,
    BLOSUM62,
    BLOSUM80,
    BLOSUM90,
    # alignment models
    AbstractScoreModel,
    AffineGapScoreModel,
    AbstractCostModel,
    CostModel,
    # pairwise alignment
    PairwiseAlignment,
    count_matches,
    count_mismatches,
    count_insertions,
    count_deletions,
    count_aligned,
    PairwiseAlignmentResult,
    pairalign,
    score,
    distance,
    alignment,
    hasalignment,
    # HTS
    SAMHeader,
    SAMRecord,
    SAMReader,
    SAMWriter,
    BAMRecord,
    BAMReader,
    BAMWriter,
    BAI,
    AuxDataDict,
    ismapped,
    refname,
    nextrefname,
    refindex,
    nextrefindex,
    leftposition,
    rightposition,
    nextleftposition,
    mappingquality,
    flag,
    templatelength,
    seqname,
    cigar,
    sequence,
    seqlength,
    qualities,
    header,
    isoverlapping,  # re-export from Bio.Intervals
    # SAM flags
    SAM_FLAG_PAIRED,
    SAM_FLAG_PROPER_PAIR,
    SAM_FLAG_UNMAP,
    SAM_FLAG_MUNMAP,
    SAM_FLAG_REVERSE,
    SAM_FLAG_MREVERSE,
    SAM_FLAG_READ1,
    SAM_FLAG_READ2,
    SAM_FLAG_SECONDARY,
    SAM_FLAG_QCFAIL,
    SAM_FLAG_DUP,
    SAM_FLAG_SUPPLEMENTARY

import Bio
using Bio.Seq
using Bio.Intervals
import Bio.StringFields: StringField
import IntervalTrees

include("operations.jl")
include("anchors.jl")
include("alignment.jl")
include("alignedseq.jl")

include("types.jl")
include("submat.jl")
include("models.jl")
include("pairwise/pairalign.jl")
include("hts/hts.jl")

include("deprecated.jl")

end


text:
using IterTools, LinearAlgebra, SparseArrays

using MultivariatePolynomials: AbstractVariable, AbstractMonomialLike,
      exponents, variables, powers, monomials, coefficient

using LazySets: Hyperrectangle, dim, low, high
using IntervalArithmetic: interval


text:
module SignatureGB

const Info = Dict{String, Int}
new_info() = Dict([("arit_ops_groebner", 0),
                   ("arit_ops_module", 0),
                   ("arit_ops_interred", 0),
                   ("interred_mat_size_rows", 0),
                   ("interred_mat_size_cols", 0),
                   ("num_zero_red", 0),
                   ("max_deg_reached", 0),
                   ("gb_size_bef_interred", 0),
                   ("gb_size_aft_interred", 0)])

mutable struct timings
    time_core_loop::Float64
    time_intermed_clean::Float64
    time_final_clean::Float64
end

include("./useful.jl")
include("./context.jl")
include("./termorder.jl")
include("./polynomials.jl")
include("./coefficients.jl")
include("./monomialtable.jl")
include("./sigtable.jl")
include("./f5data.jl")
include("./pairs.jl")
include("./symbolicpp.jl")
include("./reduction.jl")
include("./interreduction.jl")
include("./gen_example_file.jl")

export f5, decompose

# converting a vector of singular polynomials into our own data structures
function f5setup(I::Vector{P};
                 mod_order=:POT,
                 mon_order=:GREVLEX,
                 index_type=UInt32,
                 mask_type=UInt32,
                 pos_type=UInt32,
                 trace_sig_tail_tags = Symbol[],
                 max_remasks=3,
                 kwargs...) where {P <: AA.MPolyElem}

    R = parent(first(I))
    if mon_order == :GREVLEX
        order = Grevlex(Singular.nvars(R))
    else
        error("only grevlex order currently supported")
    end
    if mod_order != :POT
        error("only position over term order currently supported")
    end
    dat = f5data(I, mod_order = mod_order, trace_sig_tail_tags = trace_sig_tail_tags,
                 index_type = index_type, mask_type = mask_type,
                 pos_type = pos_type, order = order,
                 max_remasks = max_remasks)
end

# build initial pairset, basis and syzygies
function pairs_and_basis(dat::F5Data,
                         basis_length;
                         start_gen = 1)

    R = dat.R
    ctx = dat.ctx
    G = new_basis(ctx, basis_length)
    pos_t = pos_type(ctx)
    for i in 1:(start_gen - 1)
        lm = leadingmonomial(ctx, (pos_t(i), ctx.po.mo(R(1))))
        push!(G[pos_t(i)], (ctx.po.mo(R(1)), lm))
    end
    H = new_syz(ctx, basis_length)
    pairs = pairset(ctx)
    [pair!(ctx, pairs, ctx(i, R(1))) for i in start_gen:basis_length]
    G, H, pairs
end

function f5core!(dat::F5Data{I, SΓ},
                 G::Basis{I, M},
                 H::Syz{I, M},
                 pairs::PairSet{I, M, SΓ};
                 select = :deg_and_pos,
                 new_elems = new_elems_f5!,
                 interreduction = true,
                 select_both = true,
                 continous_cleanup = true,
                 verbose = 0) where {I, M, SΓ <: SigPolynomialΓ{I, M}}

    ctx = dat.ctx
    R = dat.R
    use_max_sig = (select == :deg_and_pos || select == :pos)
    orig_length = length(ctx.ord_indices)

    info_hashmap = Dict([(i, new_info()) for i in keys(ctx.ord_indices)])
    mat_cnt = 0
    global_zero_red_count = 0
    times = timings(0.0, 0.0, 0.0)
    verbose_stats = verbose > 0
    verbose_mat = verbose > 1

    curr_pos = pos(ctx, first(pairs)[1])
    curr_pos_key = first(pairs)[1][2][1]
    curr_tag = gettag(ctx, first(pairs)[1])

    non_zero_cond = I[]

    while !(isempty(pairs))
        #- INTERREDUCTION -#
        timed = @elapsed begin
            indx = pos(ctx, first(pairs)[1])
            if indx != curr_pos && mod_order(dat.ctx) == :POT

                # here we do several steps when the index jumps:
                # possibly remasking
                if dat.remasks_left > 0
                    remask!(dat.ctx.po.mo.table)
                    dat.remasks_left -= 1
                end

                # if we are in a decomposition computation we add non-zero conditions
                if curr_tag == :g_gen
                    non_zero_cond_local_sigs = H[curr_pos_key]
                    non_zero_cond_local = eltype(ctx.po)[]
                    for m in non_zero_cond_local_sigs
                        sig = (curr_pos_key, m)
                        try
                            push!(non_zero_cond_local, project(ctx, (curr_pos_key, m)))
                        catch KeyError
                            continue
                        end
                    end
                    if !(isempty(non_zero_cond_local))
                        h = random_lin_comb(ctx.po, non_zero_cond_local)
                        new_posit_key = register!(ctx, h, info_hashmap)
                        push!(non_zero_cond, new_posit_key)
                    end
                end

                if curr_tag == :f
                    for (j, i) in enumerate(non_zero_cond)
                        new_pos!(ctx, i, curr_pos + I(j), curr_pos_key, :h)
                        G[i] = Tuple{M, M}[]
                        H[i] = M[]
                        pair!(ctx, pairs, (i, one(ctx.po.mo)))
                    end
                end

                if curr_tag == :h
                    G[curr_pos_key] = Tuple{M, M}[]
                end

                # interreducing
                info_hashmap[curr_pos_key]["gb_size_bef_interred"] = gb_size(ctx, G)
                if interreduction && indx > 2
                    verbose_mat && println("-----")
                    verbose_mat && println("INTERREDUCING")
                    G, arit_ops, mat_size = interreduce(ctx, G, H, verbose = verbose_mat)
                    info_hashmap[curr_pos_key]["arit_ops_interred"] += arit_ops
                    info_hashmap[curr_pos_key]["interred_mat_size_rows"] = mat_size[1]
                    info_hashmap[curr_pos_key]["interred_mat_size_cols"] = mat_size[2]
                    verbose_mat && println("-----")
                end
                info_hashmap[curr_pos_key]["gb_size_aft_interred"] = gb_size(ctx, G)
            end
            # update current index, current tag and current index key
            curr_pos = indx
            curr_tag = gettag(ctx, first(pairs)[1])
            curr_pos_key = first(pairs)[1][2][1]

            #- PAIR SELECTION -#
            total_num_pairs = length(pairs)
            to_reduce, are_pairs, nselected, sig_degree = select!(ctx, pairs, Val(select), select_both = select_both)

            #- SYMBOLIC PP -#
            symbolic_pp_timed  = @timed symbolic_pp!(ctx, to_reduce, G, H,
                                                     use_max_sig = use_max_sig,
                                                     sig_degree = sig_degree,
                                                     max_sig_pos = curr_pos,
                                                     are_pairs = are_pairs,
                                                     enable_lower_pos_rewrite = !(interreduction))
            done = symbolic_pp_timed.value
            symbolic_pp_time = symbolic_pp_timed.time
            mat = f5matrix(ctx, done, collect(to_reduce), curr_pos, curr_pos_key, curr_tag,
                           trace_sig_tail_tags = dat.trace_sig_tail_tags,
                           enable_lower_pos_rewrite = !(interreduction))
            mat_size = (length(rows(mat)), length(mat.tbl))
            mat_dens = sum([length(rw) for rw in rows(mat)]) / (mat_size[1] * mat_size[2])

            #- REDUCTION -#
            reduction_dat = @timed reduction!(mat)

            # storing computational information
            zero_red_count = length(filter(sig_rw -> isempty(sig_rw[2]), mat.sigs_rows))
            info_hashmap[curr_pos_key]["arit_ops_groebner"] += reduction_dat.value[1]
            info_hashmap[curr_pos_key]["arit_ops_module"] += reduction_dat.value[2]
            info_hashmap[curr_pos_key]["num_zero_red"] += zero_red_count
            if curr_tag == :h
                global_zero_red_count += zero_red_count
            end

            #- PAIR GENERATION AND ADDING ELEMENTS TO BASIS/SYZYGIES -#
            pair_gen_dat = @timed new_elems(ctx, mat, pairs, G, H, info_hashmap, non_zero_cond,
                                            enable_lower_pos_rewrite = !(interreduction))
            pairs = pair_gen_dat.value
            pair_gen_time = pair_gen_dat.time

            #- PRINTING INFORMATION -#
            if verbose_mat
                mat_cnt += 1
                println("selected $(nselected) / $(total_num_pairs) pairs, in position $(curr_pos), tagged $(curr_tag), sig-degree of sel. pairs: $(sig_degree)")
                println("symbolic pp took $(symbolic_pp_time) secs.")
                println("Matrix $(mat_cnt): $(reduction_dat.time) secs reduction / size = $(mat_size) / density = $(mat_dens)")
                if !(iszero(zero_red_count))
                    println("$(zero_red_count) zero reductions at sig-degree $(sig_degree) / position $(curr_pos)")
                end
                verbose_mat && println("Pair generation took $(pair_gen_time) seconds.")
            end

            if isempty(pairs) && curr_tag == :f && !(isempty(non_zero_cond))
                for (j, i) in enumerate(non_zero_cond)
                    new_pos!(ctx, i, curr_pos + I(j), curr_pos_key, :h)
                    pair!(ctx, pairs, (i, one(ctx.po.mo)))
                    G[i] = Tuple{M, M}[]
                    H[i] = M[]
                end
                if !(isempty(pairs))
                    if interreduction
                        verbose_mat && println("interreducing")
                        info_hashmap[curr_pos_key]["gb_size_bef_interred"] = gb_size(ctx, G)
                        G, arit_ops, mat_size = interreduce(ctx, G, H, verbose = verbose_mat)
                        info_hashmap[curr_pos_key]["interred_mat_size_rows"] = mat_size[1]
                        info_hashmap[curr_pos_key]["interred_mat_size_cols"] = mat_size[2]
                        info_hashmap[curr_pos_key]["arit_ops_interred"] += arit_ops
                        info_hashmap[curr_pos_key]["gb_size_aft_interred"] = gb_size(ctx, G)
                    end
                    curr_pos = pos(ctx, first(pairs)[1])
                    curr_tag = gettag(ctx, first(pairs)[1])
                    curr_pos_key = first(pairs)[1][2][1]
                end
            end
        end
        # Split up timings between core loop and cleanup
        if curr_tag in [:f, :g, :g_gen]
            times.time_core_loop += timed
        elseif iszero(f_left(ctx, curr_pos))
            times.time_final_clean += timed
        else
            times.time_intermed_clean += timed
        end
    end
    verbose_mat && println("-----")

    if ctx.ord_indices[curr_pos_key][:tag] == :h
        G[curr_pos_key] = Tuple{M, M}[]
    end

    # interreducing
    timed = @elapsed begin
        if interreduction
            verbose_mat && println("FINAL INTERREDUCTION STEP")
            info_hashmap[curr_pos_key]["gb_size_bef_interred"] = gb_size(ctx, G)
            G, arit_ops, mat_size = interreduce(ctx, G, H, verbose = verbose_mat)
            info_hashmap[curr_pos_key]["interred_mat_size_rows"] = mat_size[1]
            info_hashmap[curr_pos_key]["interred_mat_size_cols"] = mat_size[2]
            info_hashmap[curr_pos_key]["arit_ops_interred"] += arit_ops
            info_hashmap[curr_pos_key]["gb_size_aft_interred"] = gb_size(ctx, G)
        end
    end
    times.time_core_loop += timed

    #- PRINTING INFORMATION -#
    verbose_stats && println("-----")
    num_arit_operations_groebner = 0
    num_arit_operations_module_overhead = 0
    num_arit_operations_interreduction = 0
    num_arit_operations_cleanup = 0
    kys = sort(collect(keys(info_hashmap)), by = k -> ctx.ord_indices[k][:position])
    for k in kys
        info = info_hashmap[k]
        arit_ops_groebner = info["arit_ops_groebner"]
        arit_ops_module = info["arit_ops_module"]
        arit_ops_interred = info["arit_ops_interred"]
        interred_mat_size_rows = info["interred_mat_size_rows"]
        interred_mat_size_cols = info["interred_mat_size_cols"]
        num_zero_red = info["num_zero_red"]
        gb_size_bef_interred = info["gb_size_bef_interred"]
        gb_size_aft_interred = info["gb_size_aft_interred"]
        max_degree = info["max_deg_reached"]
        generator_poly = ctx((k, one(ctx.po.mo)))[:poly]
        isempty(generator_poly.mo) ? deg = -1 : deg = degree(ctx.po, generator_poly)
        num_arit_operations_groebner += arit_ops_groebner
        num_arit_operations_module_overhead += arit_ops_module
        num_arit_operations_interreduction += arit_ops_interred
        if ctx.ord_indices[k][:tag] in [:h, :p, :p_gen]
            num_arit_operations_cleanup += (arit_ops_groebner + arit_ops_interred + arit_ops_module)
        end
        if verbose_stats
            if k <= orig_length
                println("INFO for index $(ctx.ord_indices[k][:position]), original index $(k), tagged $(ctx.ord_indices[k][:tag]), degree $(deg)")
            else
                println("INFO for index $(ctx.ord_indices[k][:position]), tagged $(ctx.ord_indices[k][:tag]), degree $(deg)")
            end
            println("Arithmetic operations in GB computation:              $(arit_ops_groebner)")
            println("Arithmetic operations in module overhead:             $(arit_ops_module)")
            println("Arithmetic operations in interreduction:              $(arit_ops_interred)")
            println("Size of interreduction matrix:                        $((interred_mat_size_rows, interred_mat_size_cols))")
            println("GB size before interreduction:                        $(gb_size_bef_interred)")
            println("GB size after interreduction:                         $(gb_size_aft_interred)")
            println("Maximal degree reached:                               $(max_degree)")
            println("Number of zero reductions:                            $(num_zero_red)")
            println("-----")
        end
    end
    if verbose_stats
        total_num_arit_ops = num_arit_operations_groebner + num_arit_operations_module_overhead + num_arit_operations_interreduction
        println("total arithmetic operations (groebner):               $(num_arit_operations_groebner)")
        println("total arithmetic operations (module overhead):        $(num_arit_operations_module_overhead)")
        println("total arithmetic operations (interreduction):         $(num_arit_operations_interreduction)")
        println("final number of arithmetic operations (core loop):    $(total_num_arit_ops - num_arit_operations_cleanup)")
        println("final number of arithmetic operations (cleanup step): $(num_arit_operations_cleanup)")
        println("final number of arithmetic operations (total):        $(total_num_arit_ops)")
        println("final size of GB:                                     $(gb_size(dat.ctx, G))")
        println("-----")
        println("time spent in core loop:                              $(times.time_core_loop)")
        println("time spent in intermediate cleanup:                   $(times.time_intermed_clean)")
        println("time spent in final cleanup:                          $(times.time_final_clean)")
    end

    verbose_stats && println("saturation added $(global_zero_red_count) elements not in the original ideal.")
    total_num_arit_ops = num_arit_operations_groebner + num_arit_operations_module_overhead + num_arit_operations_interreduction

    G, total_num_arit_ops
end

function f5(I::Vector{P};
            start_gen = 1,
            trace_sig_tail_tags = Symbol[],
            mod_order=:POT,
            mon_order=:GREVLEX,
            index_type=UInt32,
            mask_type=UInt32,
            pos_type=UInt32,
            select = :deg_and_pos,
            verbose = 0,
            interreduction = true,
            max_remasks = 3,
            kwargs...) where {P <: AA.MPolyElem}

    R = parent(first(I))
    dat = f5setup(I, mod_order = mod_order,
                  mon_order = mon_order, index_type = index_type,
                  mask_type = mask_type, pos_type = pos_type,
                  trace_sig_tail_tags = trace_sig_tail_tags,
                  max_remasks = max_remasks, kwargs...)
    G, H, pairs = pairs_and_basis(dat, length(I), start_gen = start_gen)
    G, total_num_arit_ops = f5core!(dat, G, H, pairs, select = select, interreduction = interreduction, verbose = verbose)
    [R(dat.ctx, (i, g[1])) for i in keys(G) for g in G[i]]
end

function decompose(I::Vector{P};
                   start_gen = 1,
                   mod_order=:POT,
                   mon_order=:GREVLEX,
                   index_type=UInt32,
                   mask_type=UInt32,
                   pos_type=UInt32,
                   select = :deg_and_pos,
                   verbose = 0,
                   interreduction = true,
                   max_remasks = 3,
                   kwargs...) where {P <: AA.MPolyElem}

    R = parent(first(I))
    dat = f5setup(I, mod_order = mod_order,
                  mon_order = mon_order, index_type = index_type,
                  mask_type = mask_type, pos_type = pos_type,
                  trace_sig_tail_tags = [:f, :g_gen, :h],
                  max_remasks = max_remasks, kwargs...)
    G, H, pairs = pairs_and_basis(dat, length(I), start_gen = start_gen)
    G, total_num_arit_ops = f5core!(dat, G, H, pairs, select = select, verbose = verbose,
                                    new_elems = new_elems_decomp!, select_both = false, interreduction = interreduction)
    [R(dat.ctx, (i, g[1])) for i in keys(G) for g in G[i]]
end
end


text:
commands = ["F10", "N3", "F7", "R90", "F11"]
commands = split(strip(read("input.txt", String)), "\n")

mutable struct Position
    x::Int64
    y::Int64
    bearing::Int64
end

pos = Position(0, 0, 90)
commandmap = Dict(
    'N' => v -> pos.y += v,
    'S' => v -> pos.y -= v,
    'W' => v -> pos.x -= v,
    'E' => v -> pos.x += v,
    'R' => v -> pos.bearing = (pos.bearing + v) % 360,
    'L' => v -> pos.bearing = (360 + pos.bearing - v) % 360,
    'F' => v -> commandmap["NESW"[(pos.bearing ÷ 90) + 1]](v))

for command in commands
    commandmap[command[1]](parse(Int, command[2:end]))
    println(pos)
end

print(abs(pos.x) + abs(pos.y)) # 521

ship = Position(0, 0, 0)
waypoint = Position(10, 1, 0)

function rotate(n)
    prev = Position(waypoint.x, waypoint.y, 0)
    if n == 90
        waypoint.x = prev.y
        waypoint.y = -prev.x
    elseif n == 180
        waypoint.x = -waypoint.x
        waypoint.y = -waypoint.y
    elseif n == 270
        waypoint.x = -prev.y
        waypoint.y = prev.x
    end
end

commanddict = Dict(
    'N' => v -> waypoint.y += v,
    'S' => v -> waypoint.y -= v,
    'W' => v -> waypoint.x -= v,
    'E' => v -> waypoint.x += v,
    'R' => rotate,
    'L' => v -> rotate(360 - v),
    'F' => function (v)
    diff = (x = waypoint.x - ship.x,)
    ship.x += waypoint.x * v
    ship.y += waypoint.y * v
end)

for command in commands
    commanddict[command[1]](parse(Int, command[2:end]))
end

print(abs(ship.x) + abs(ship.y))

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
abstract type FluidFlowProblem end
abstract type SteadyStateProblem <: FluidFlowProblem end
abstract type TimeDependantProblem <: FluidFlowProblem end

abstract type InitialValueProblem <: FluidFlowProblem end
abstract type DoubleDistributionProblem <: FluidFlowProblem end

has_external_force(::FluidFlowProblem) = false
function velocity_gradient(
    ::FluidFlowProblem,
    x::T,
    y::T,
    timestep::Real = 0.0,
) where {T <: Real}
    return zeros(T, 2, 2)
end

function range(problem::FluidFlowProblem)
    Δx = problem.domain_size[1] / problem.NX
    Δy = problem.domain_size[2] / problem.NY

    x_range = range(Δx / 2, stop = problem.domain_size[1] - Δx / 2, length = problem.NX)
    y_range = range(Δy / 2, stop = problem.domain_size[1] - Δy / 2, length = problem.NY)

    return x_range, y_range
end

boundary_conditions(problem::FluidFlowProblem) = BoundaryCondition[]

function deviatoric_tensor(
    q::Quadrature,
    problem::FluidFlowProblem,
    x::T,
    y::T,
    time::Real = 0.0,
) where {T <: Real}
    a = velocity_gradient(problem, x, y, time)
    ν = viscosity(problem)

    σ = -ν * [
        2 * a[1, 1] a[1, 2] + a[2, 1]
        a[1, 2] + a[2, 1] 2 * a[2, 2]
    ]

    return σ
end

function pressure_tensor(
    q::Quadrature,
    problem::FluidFlowProblem,
    x::T,
    y::T,
    time::Real = 0.0,
) where {T <: Real}
    A = problem.A
    B = problem.B

    p = pressure(q, problem, x, y, time) * I
    return p - deviatoric_tensor(q, problem, x, y, time)
end

function force(
    problem::FluidFlowProblem,
    x_idx::Int,
    y_idx::Int,
    time::Real = 0.0,
) where {Int <: Integer}
    x_range, y_range = range(problem)

    x = x_range[x_idx]
    y = y_range[y_idx]

    return force(problem, x, y, time)
end

is_steady_state(problem::FluidFlowProblem) = problem.static
is_time_dependant(problem::FluidFlowProblem) = !problem.static

# Dimensionless
viscosity(problem::FluidFlowProblem) = problem.ν * delta_x(problem)^2 / delta_t(problem)
heat_diffusion(problem::FluidFlowProblem) =
    problem.κ * delta_x(problem)^2 / delta_t(problem)
reynolds(problem::FluidFlowProblem) = problem.NY * problem.u_max / problem.ν

function delta_t(problem::FluidFlowProblem)
    return delta_x(problem) * problem.u_max
end

function delta_x(problem::FluidFlowProblem)
    if (problem.NX > problem.NY)
        return problem.domain_size[1] * (1 / problem.NX)
    end

    return problem.domain_size[2] * (1 / problem.NY)
end

lattice_viscosity(problem) = problem.ν #::FluidFlowProblem)
lattice_density(q, problem::FluidFlowProblem, x, y, t = 0.0) = density(q, problem, x, y, t)
lattice_velocity(q, problem::FluidFlowProblem, x, y, t = 0.0) =
    problem.u_max * velocity(problem, x, y, t)
lattice_pressure(q, problem::FluidFlowProblem, x, y, t = 0.0) =
    problem.u_max^2 * pressure(q, problem, x, y, t)
lattice_force(problem::FluidFlowProblem, x, y, t = 0.0) =
    problem.u_max * delta_t(problem) * force(problem, x, y, t)
lattice_temperature(q, problem::FluidFlowProblem, x, y, t = 0.0) =
    pressure(q, problem, x, y) / density(q, problem, x, y)

dimensionless_viscosity(problem) = problem.ν * delta_x(problem)^2 / delta_t(problem)
dimensionless_density(problem::FluidFlowProblem, ρ) = ρ
dimensionless_velocity(problem::FluidFlowProblem, u) = u / problem.u_max
dimensionless_pressure(q, problem::FluidFlowProblem, p) = p# (p - 1.0) / (q.speed_of_sound_squared * problem.u_max^2 )
dimensionless_temperature(q, problem::FluidFlowProblem, T) = T #* q.speed_of_sound_squared
dimensionless_force(problem::FluidFlowProblem, F) = F / (problem.u_max * delta_t(problem))
dimensionless_stress(problem::FluidFlowProblem, σ) = begin
    # Rescale to dimensionless number (TODO check why problem.u_max)
    factor = 1 / problem.u_max^2

    return σ * factor
end

function equilibrium(q::Quadrature, problem::FluidFlowProblem, x::T, y::T) where {T <: Real}
    equilibrium(
        q,
        lattice_density(q, problem, x, y),
        lattice_velocity(q, problem, x, y),
        lattice_temperature(q, problem, x, y),
    )
end

include("taylor_green_vortex.jl")
include("decaying_shear_flow.jl")
include("poiseuille.jl")
include("couette_flow.jl")
include("lid_driven_cavity.jl")
include("linear_hydrodynamics_modes.jl")
include("second_order_convergence.jl")

# error(::Val{:density}, node, solution) = density(node) - density(solution)


text:
abstract type AbstractObservable end
abstract type AbstractSimpleObservable <: AbstractObservable end
abstract type AbstractCompositeObservable <: AbstractObservable end

const axis_dict = Dict{Symbol,Int64}(:x=>1, :y=>2, :z=>3)

struct QuantumObservable{TO<:AbstractSparseMatrix} <: AbstractSimpleObservable
    O::TO
    temp_vector::Vector{Complex{Float64}}
end
function QuantumObservable(Dh::Int64, q_spins::Vector{Int64}, sigma::Int64)
    assert(sigma in 1:3)
    for i in eachindex(q_spins)
        assert(2^q_spins[i]<=Dh)
    end
    QuantumObservable(build_Spin_Operator(Dh, q_spins, sigma), Vector{Complex{Float64}}(Dh))
end

struct ClassicalObservable <: AbstractSimpleObservable
    sigma::Int64
    cl_spins::Vector{Int64}
end

struct EmptyObservable <: AbstractSimpleObservable end

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

function (QO::QuantumObservable)(v::Vector{Complex{Float64}})
    A_mul_B!(QO.temp_vector, QO.O, v)
    return real(dot(v, QO.temp_vector))
end
function (QO::QuantumObservable)(v1::Vector{Complex{Float64}}, v2::Vector{Complex{Float64}})
    A_mul_B!(QO.temp_vector, QO.O, v2)
    return real(dot(v1, QO.temp_vector))
end
function (CO::ClassicalObservable)(v::VectorOfArray{Float64, 2, Vector{Vector{Float64}}})
    s = zero(Float64)
    sigma = CO.sigma
    for cl_spin in CO.cl_spins
        s += v[cl_spin, sigma]
    end
    return s
end
function (CO::ClassicalObservable)(v::VectorOfArray{T, 2, Vector{Vector{T}}}) where T<:Complex{Float64}
    s = zero(Float64)
    sigma = CO.sigma
    for cl_spin in CO.cl_spins
        s += real(v[cl_spin, sigma])
    end
    return s
end

@inline (EO::EmptyObservable)(v) = 0.0

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

convert_axis(axis) = error("no method for such input type of axis")
convert_axis(axis::Symbol) = axis_dict[axis]

struct ExactObservable{T1,T2,TO} <: AbstractCompositeObservable
    axis::T1
    name::T2
    QO::QuantumObservable{TO}
    id::Int64
    function ExactObservable(A::ExactApprox, spins::Vector{NTuple{2,Int64}}, axis::T1, name::T2, id::Int64) where{T1,T2}
        QO = QuantumObservable(get_Dh(A), get_positions(A, spins), convert_axis(axis))
                               new{T1,T2,typeof(QO.O)}(axis, name, QO, id)
    end
end
@inline ExactObservable(A, spins, axis, name) = ExactObservable(A, spins, axis, name, axis_dict[axis]+1)

struct ClusteredObservable{T1,T2,TO} <: AbstractCompositeObservable
    axis::T1
    name::T2
    QO::QuantumObservable{TO}
    function ClusteredObservable(A::ClusteredApprox, spins::Vector{NTuple{2,Int64}}, axis::T1, name::T2) where{T1,T2}
        spins_c = Set{Int64}()
        for pos in get_positions(A, spins)
            push!(spins_c, pos[1])
        end
        spins_cc = Int64[spins_c...]
        sort!(spins_cc)
        QO = QuantumObservable(get_Dh(A), spins_cc, convert_axis(axis))
        new{T1,T2,typeof(QO.O)}(axis, name, QO)
    end
end

struct PureClassicalObservable{T1,T2} <: AbstractCompositeObservable
    axis::T1
    name::T2
    CO::ClassicalObservable
    function PureClassicalObservable(A::PureClassicalApprox, spins::Vector{NTuple{2,Int64}}, axis::T1, name::T2) where{T1,T2}
        new{T1,T2}(axis, name, ClassicalObservable(convert_axis(axis), get_positions(A, spins)))
    end
end

struct HybridObservable{T1,T2, O1<:Union{QuantumObservable, EmptyObservable},O2<:Union{ClassicalObservable, EmptyObservable}} <: AbstractCompositeObservable
    axis::T1
    name::T2
    c::Float64
    QO::O1
    CO::O2
    function HybridObservable(A::HybridApprox, spins::Vector{NTuple{2,Int64}}, axis::T1, name::T2) where{T1,T2}
        q_poss, cl_poss = get_positions(A, spins)
        if !isempty(q_poss)
            QO = QuantumObservable(get_Dh(A), q_poss, convert_axis(axis))
        else
            QO = EmptyObservable()
        end
        if !isempty(cl_poss)
            CO = ClassicalObservable(convert_axis(axis), cl_poss)
        else
            CO = EmptyObservable()
        end
        c = sqrt(get_Dh(A)+1)/2
        new{T1,T2, typeof(QO), typeof(CO)}(axis, name, c, QO, CO)
    end
end

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

@inline (EO::ExactObservable)(v::Vector{Complex{Float64}}) = EO.QO(v)
@inline (EO::ExactObservable)(voa::VectorOfArray{Complex{Float64},2,Array{Array{Complex{Float64},1},1}}) = EO.QO(voa[1],voa[EO.id])

@inline (PCO::PureClassicalObservable)(v::VectorOfArray{Float64, 2, Vector{Vector{Float64}}}) = PCO.CO(v)
function (CLO::ClusteredObservable)(v::VectorOfArray{Complex{Float64}, 2, Vector{Vector{Complex{Float64}}}})
    s = zero(Float64)
    for i in eachindex(v)
        @inbounds s = s + CLO.QO(v[i])
    end
    return s
end
@inline function (HO::HybridObservable)(v::ArrayPartition{Complex{Float64}, Tuple{Vector{Complex{Float64}}, VectorOfArray{T, 2, Vector{Vector{T}}}}}) where T<:Union{Float64,Complex{Float64}}
    return HO.c*HO.QO(v.x[1]) + HO.CO(v.x[2])
end

@inline function (HO::HybridObservable)(v::ArrayPartition{Complex{Float64},Tuple{Array{Complex{Float64},1},VectorOfArray{Float64,2,Array{Array{Float64,1},1}},VectorOfArray{Float64,2,Array{Array{Float64,1},1}}}})
    return HO.c*HO.QO(v.x[1]) + HO.CO(v.x[2])
end

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

@inline build_Observable(A::ExactApprox, spins::Vector{NTuple{2,Int64}}, axis, name) = ExactObservable(A, spins, axis, name)
@inline build_Observable(A::ExactApprox, spins::Vector{NTuple{2,Int64}}, axis, name, id::Int64) = ExactObservable(A, spins, axis, name, id)
@inline build_Observable(A::ClusteredApprox, spins::Vector{NTuple{2,Int64}}, axis, name) = ClusteredObservable(A, spins, axis, name)
@inline build_Observable(A::PureClassicalApprox, spins::Vector{NTuple{2,Int64}}, axis, name) = PureClassicalObservable(A, spins, axis, name)
@inline build_Observable(A::HybridApprox, spins::Vector{NTuple{2,Int64}}, axis, name) = HybridObservable(A, spins, axis, name)

@inline build_Observable(A::AbstractApproximation, name::Symbol, axis::Symbol) = build_Observable(A, get_spins_by_name(A, name), axis, name)
@inline build_Observable(A::AbstractApproximation, s_array::SpinArray, axis::Symbol) = build_Observable(A, s_array.spins, axis, s_array.name)

@inline build_Observable(A::ExactApprox, name::Symbol, axis::Symbol, id::Int64) = build_Observable(A, get_spins_by_name(A, name), axis, name, id)
@inline build_Observable(A::ExactApprox, s_array::SpinArray, axis::Symbol, id::Int64) = build_Observable(A, s_array.spins, axis, s_array.name, id)

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

struct ObservablesSet{T<:Tuple} <: AbstractObservable
    Observables::T
    len::Int64
    ObservablesSet(O_tuple::T) where {T} = new{T}(O_tuple, length(O_tuple))
end
ObservablesSet(O::AbstractCompositeObservable) = ObservablesSet((O,))

@inline ObservablesSet(obs...) = ObservablesSet(obs)
@inline function (OSET::ObservablesSet)(t, v, integrator)
    data = Vector{Float64}(OSET.len)
    for i in eachindex(OSET.Observables)
        data[i] = OSET.Observables[i](v)
    end
    return data
end

@inline return_type(O::AbstractCompositeObservable) = Float64
@inline return_type(O::ObservablesSet) = Vector{Float64}

#########################################
#Extra machinery to define couplings for calculating correlation functions

struct ConvRules
    str_vec::Vector{String}
    rules::Vector{NTuple{2,Int64}}
end

@inline get_rules(cfrules::ConvRules) = cfrules.rules

#For Exact Approximation I should make a special type of CFVals and dispatch on ExactApprox to choose this particular type.

struct CFVals
    str_vec::Vector{String}
    ts::Vector{Float64}
    data::VectorOfArray{Float64, 2, Vector{Vector{Float64}}}
    errors::VectorOfArray{Float64, 2, Vector{Vector{Float64}}}
end

struct CFVals0
    str_vec::Vector{String}
    ts::Vector{Float64}
    data::VectorOfArray{Float64, 2, Vector{Vector{Float64}}}
end

struct CFArray
    str_vec::Vector{String}
    ts::Vector{Float64}
    datas::Dict{String, Matrix{Float64}}
end

import Base.similar
import Base.zeros

@inline similar(cf::CFVals) = CFVals(copy(cf.str_vec), copy(cf.ts), similar(cf.data), similar(cf.errors))
@inline similar(cf0::CFVals0) = CFVals0(copy(cf0.str_vec), copy(cf0.ts), similar(cf0.data))

function zeros(cf::CFVals)
    temp = similar(cf)
    fill!(temp.data, 0.0)
    fill!(temp.errors, 0.0)
    return temp
end
function zeros(cf0::CFVals0)
    temp = similar(cf0)
    fill!(temp.data, 0.0)
    return temp
end

function CFVals(rules::ConvRules, ts::Vector{Float64})
    data = VectorOfArray([zeros(length(ts)) for rule in get_rules(rules)])
    errors = VectorOfArray([zeros(length(ts)) for rule in get_rules(rules)])
    CFVals(rules.str_vec, ts, data, errors)
end
function CFVals0(rules::ConvRules, ts::Vector{Float64})
    data = VectorOfArray([zeros(length(ts)) for rule in get_rules(rules)])
    CFVals0(rules.str_vec, ts, data)
end

function make_oset(A::AbstractApproximation, links::NTuple{D, NTuple{2,Union{Symbol, SpinArray}}}, axes::Vector{Symbol}) where D
    os = Set{Tuple{Union{Symbol, SpinArray},Symbol}}()
    for axis in axes
        for link in links
            push!(os, (link[1], axis))
            push!(os, (link[2], axis))
        end
    end
    os_a = [os...]
    sort!(os_a)
    return os_a
end

function make_oset(A::ExactApprox, linkss::NTuple{D, NTuple{2,Union{Symbol, SpinArray}}}, axes::Vector{Symbol}) where D
    links = ((:all,:all),)
    os = Set{Tuple{Symbol,Symbol,Int64}}()
    for axis_i in eachindex(axes)
        for link in links
            push!(os, (link[1], axes[axis_i], axis_i+1))
        end
    end
    os_a = [os...]
    sort!(os_a)
    return os_a
end

function set_CorrFuncs(A::AbstractApproximation, links::NTuple{D, NTuple{2,Union{Symbol, SpinArray}}}, axes::Vector{Symbol}) where D
    os_a = make_oset(A, links, axes)
    print("$os_a\n")
    OSET = ObservablesSet([build_Observable(A, O...) for O in os_a]...)
    rules = Vector{NTuple{2,Int64}}()
    links_str = Vector{String}()
    for axis in axes
        for i in eachindex(links)
            push!(links_str, get_string((links[i], axis)))
            push!(rules, ( findfirst(os_a, (links[i][1], axis)) , findfirst(os_a, (links[i][2], axis)) ) )
        end
    end
    return ConvRules(links_str, rules), OSET
end
set_CorrFuncs(A::CompositeApproximation, links::NTuple{D, NTuple{2,Union{Symbol, SpinArray}}}, axes::Vector{Symbol}) where D = set_CorrFuncs(A.A1, links, axes)

@inline set_CorrFuncs(A::abstractApproximation, links::NTuple{2,Union{Symbol, SpinArray}} = (:all,:all), axes::Vector{Symbol} = [:x, :y,:z]) = set_CorrFuncs(A, tuple(links), axes)

function calculateCorrFunc!(A::ExactApprox, saved_values::SavedValues, rules::ConvRules, vals::CFVals)
    savevals = saved_values.saveval
    rule = get_rules(rules)
    for i in eachindex(rule)
        r1 = rule[i][1]
        @inbounds for j = 1:length(vals.data[1])
            vals.data[j,i] += savevals[j][i]
        end
    end
end

function calculateCorrFunc!(A::AbstractApproximation, saved_values::SavedValues, rules::ConvRules, vals::CFVals, iter)
     meanvals = Vector{Vector{Float64}}()
     savevals = saved_values.saveval
     l = length(savevals)
     nmeans = length(savevals[1])
     for m in Base.OneTo(nmeans)
         v = Vector{Float64}(l)
         for i in Base.OneTo(l)
             v[i] = savevals[i][m]
         end
         push!(meanvals, v)
     end
     #####
     data = vals.data
     errors = vals.errors
     rule = get_rules(rules)
     for i in eachindex(rule)
        r1 = rule[i][1]
        r2 = rule[i][2]
        if r1 == r2
            convolute_s!(meanvals[r1], data[i], errors[i], iter)
        else
            convolute_d!(meanvals[r1], meanvals[r2], data[i], errors[i], iter)
        end
     end
end

function calculateCorrFunc!(A::AbstractApproximation, saved_values::SavedValues, rules::ConvRules, vals::CFVals0, iter)
     meanvals = Vector{Vector{Float64}}()
     savevals = saved_values.saveval
     l = length(savevals)
     nmeans = length(savevals[1])
     for m in Base.OneTo(nmeans)
         v = Vector{Float64}(l)
         for i in Base.OneTo(l)
             v[i] = savevals[i][m]
         end
         push!(meanvals, v)
     end
     #####
     data = vals.data
     rule = get_rules(rules)
     for i in eachindex(rule)
        r1 = rule[i][1]
        r2 = rule[i][2]
        if r1 == r2
            convolute_s!(meanvals[r1], data[i], iter)
        else
            convolute_d!(meanvals[r1], meanvals[r2], data[i], iter)
        end
     end
end

@inline calculateCorrFunc!(A::CompositeApproximation, saved_values::SavedValues, rules::ConvRules, vals::T, iter) where {T<:Union{CFVals, CFVals0}} = calculateCorrFunc!(A.A1, saved_values, rules, vals, iter)

function convolute_d!(mean1::Vector{Float64}, mean2::Vector{Float64}, data::Vector{Float64})
    delta = 0
    stop = length(mean1)
    maxDelta = length(data) - 1
    @inbounds while delta <= maxDelta
        stmd = stop-delta
        val = dot(mean1, 1:stmd, mean2, (1+delta):stop)
        val += dot(mean2, 1:stmd, mean1, (1+delta):stop)
        data[delta+1] += val*0.5/stmd
        delta += 1
    end
end

function convolute_s!(mean::Vector{Float64}, data::Vector{Float64})
    delta = 0
    stop = length(mean)
    maxDelta = length(data) - 1
    @inbounds while delta <= maxDelta
        stmd = stop-delta
        data[delta+1] += dot(mean, 1:stmd, mean, (1+delta):stop)/stmd
        delta += 1
    end
end

function convolute_d!(mean1::Vector{Float64}, mean2::Vector{Float64}, data::Vector{Float64}, errors::Vector{Float64}, iter)
    delta = 1
    stop = length(mean1)
    maxDelta = length(data)
    @inbounds while delta <= maxDelta
        stmd = stop-delta+1
        val = dot(mean1, 1:stmd, mean2, delta:stop)
        val += dot(mean2, 1:stmd, mean1, delta:stop)
        val *= 0.5/stmd
        d1 = (val - data[delta])
        data[delta] += d1/iter
        d2 = (val - data[delta])
        errors[delta] += d1*d2
        delta += 1
    end
end

function convolute_s!(mean::Vector{Float64}, data::Vector{Float64}, errors::Vector{Float64}, iter)
    delta = 1
    stop = length(mean)
    maxDelta = length(data)
    @inbounds while delta <= maxDelta
        stmd = stop-delta+1
        val = dot(mean, 1:stmd, mean, delta:stop)/stmd
        d1 = (val-data[delta])
        data[delta] += d1/iter
        d2 = (val-data[delta])
        errors[delta] += d1*d2
        delta += 1
    end
end

function convolute_d!(mean1::Vector{Float64}, mean2::Vector{Float64}, data::Vector{Float64}, iter)
    delta = 1
    stop = length(mean1)
    maxDelta = length(data)
    @inbounds while delta <= maxDelta
        stmd = stop-delta+1
        val = dot(mean1, 1:stmd, mean2, delta:stop)
        val += dot(mean2, 1:stmd, mean1, delta:stop)
        val *= 0.5/stmd
        data[delta] = val
        delta += 1
    end
end

function convolute_s!(mean::Vector{Float64}, data::Vector{Float64}, iter)
    delta = 1
    stop = length(mean)
    maxDelta = length(data)
    @inbounds while delta <= maxDelta
        stmd = stop-delta+1
        data[delta] = dot(mean, 1:stmd, mean, delta:stop)/stmd
        delta += 1
    end
end


text:
function to_string(msg)
    if isa(msg, AbstractString)
        msg # pass-through

    elseif isa(msg, Expr) || isa(msg, Symbol) # && !isempty(msg)
        # message is an expression needing evaluating
        :(Main.Base.string($(esc(msg))))

    else
        Main.Base.string(msg)

    end #if
end #function

"""
custom abstract `Exception` that is the parent of all Exception in the `GeoEfficiency` package
"""
abstract type GeoException <: Exception end

import Base: showerror

showerror(io::IO, err::GeoException) = print(io, nameof(typeof(err)), ": ", err.msg)

"""
custom `Exception` indicating invalid radiation detector dimensions
"""
struct  InValidDetectorDim <: GeoException
   	msg::AbstractString
end

"""

    @validateDetector cond [text]

throw an [`InValidDetectorDim`](@ref) if `cond` is `false`.
Message `text` is optionally displayed upon validation failure.

# Examples
```jldoctest
julia> @validateDetector iseven(3) "3 is an odd number!"
ERROR: InValidDetectorDim: 3 is an odd number!

julia> @validateDetector isodd(3) "What even are numbers?"
```
"""
macro validateDetector(ex, msgs...)
    msg = isempty(msgs) ? "'$ex' is not satisfied" : to_string(msgs[1])
    return :($(esc(ex)) ? $(nothing) : throw(InValidDetectorDim($msg)))
end

"""
custom `Exception` indicating a source-to-detector geometry which may be valid but not implemented yet
"""
struct  NotImplementedError <: GeoException
   	msg::AbstractString
end

"""

    @notImplementedError [msg]

custom macro to throw [`NotImplementedError`](@ref) `Exception`.
"""
macro notImplementedError(msgs...)
    msg = isempty(msgs) ? "" : to_string(msgs[1])
    return :(throw(NotImplementedError($msg)))
end

"""
custom `Exception` indicating a not valid source to detector geometry
"""
struct  InValidGeometry <: GeoException
   	msg::AbstractString
end

"""

    @inValidGeometry [msg]

custom macro to throw [`NotImplementedError`](@ref) `Exception`.
"""
macro inValidGeometry(msgs...)
    msg = isempty(msgs) ? "" : to_string(msgs[1])
    return :(throw(InValidGeometry($msg)))
end


text:
export NormalOp, normalOperator

struct NormalOp{S,D}
  parent::S
  weights::D
end

function Base.copy(S::NormalOp)
  return NormalOp(copy(S.parent), S.weights)
end

function normalOperator(S, W=opEye())
  return NormalOp(S,W)
end

function Base.size(S::NormalOp)
  return (S.parent.ncol, S.parent.ncol)
end

function Base.size(S::NormalOp, dim)
  if dim == 1 || dim == 2
    return S.parent.ncol
  else
    error()
  end
end

function LinearAlgebra.mul!(x, S::NormalOp, b)
  x .= S * b
  return x
end

# Generic fallback -> TODO avoid allocations
function Base.:*(N::NormalOp, x::AbstractVector)

  #@info size(x)  size(N.parent)   size(N.weights)

  d = adjoint(N.parent)*(N.weights*(N.parent*x))
  return d
end


text:
""" 2-category of finitely presented categories.

This module is for the 2-category **Cat** what the module [`FinSets](@ref) is
for the category **Set**: a finitary, combinatorial setting where explicit
calculations can be carried out. We emphasize that the prefix `Fin` means
"finitely presented," not "finite," as finite categories are too restrictive a
notion for many purposes. For example, the free category on a graph is finite if
and only if the graph is DAG, which is a fairly special condition. This usage of
`Fin` is also consistent with `FinSet` because for sets, being finite and being
finitely presented are equivalent.
"""
module FinCats
export FinCat, FinCatGraph, Path, ob_generator, hom_generator,
  ob_generator_name, hom_generator_name, ob_generators, hom_generators,
  equations, is_discrete, is_free, graph, edges, src, tgt, presentation,
  FinFunctor, FinDomFunctor, is_functorial, collect_ob, collect_hom, force,
  FinTransformation, components, is_natural, is_initial

using AutoHashEquals
using Reexport
using StaticArrays: SVector
using DataStructures: IntDisjointSets, in_same_set, num_groups

@reexport using ..Categories
using ...GAT, ...Present, ...Syntax
import ...Present: equations
using ...Theories: Category, Schema, ObExpr, HomExpr, AttrExpr, AttrTypeExpr
import ...Theories: dom, codom, id, compose, ⋅, ∘
using ...CSetDataStructures, ...Graphs
import ...Graphs: edges, src, tgt, enumerate_paths
import ..Categories: CatSize, ob, hom, ob_map, hom_map, component, op

# Categories
############

""" Size of a finitely presented category.
"""
struct FinCatSize <: CatSize end

""" A finitely presented (but not necessarily finite!) category.
"""
const FinCat{Ob,Hom} = Cat{Ob,Hom,FinCatSize}

FinCat(g::HasGraph, args...; kw...) = FinCatGraph(g, args...; kw...)
FinCat(pres::Presentation, args...; kw...) =
  FinCatPresentation(pres, args...; kw...)

""" Object generators of finitely presented category.

The object generators of finite presented category are almost always the same as
the objects. In principle, however, it is possible to have equations between
objects, so that there are fewer objects than object generators.
"""
function ob_generators end

""" Morphism generators of finitely presented category.
"""
function hom_generators end

""" Coerce or look up object generator in a finitely presented category.

Because object generators usually coincide with objects, the default method for
[`ob`](@ref) in finitely presented categories simply calls this function.
"""
function ob_generator end

ob(C::FinCat, x) = ob_generator(C, x)

""" Coerce or look up morphism generator in a finitely presented category.

Since morphism generators often have a different data type than morphisms (e.g.,
in a free category on a graph, the morphism generators are edges and the
morphisms are paths), the return type of this function is generally different
than that of [`hom`](@ref).
"""
function hom_generator end

""" Name of object generator, if any.

When object generators have names, this function is a one-sided inverse to
[`ob_generator`](@ref) in that `ob_generator(C, ob_generator_name(C, x)) == x`.
"""
function ob_generator_name end

""" Name of morphism generator, if any.

When morphism generators have names, this function is a one-sided inverse to
[`hom_generator`](@ref). See also: [`ob_generator_name`](@ref).
"""
function hom_generator_name end

""" Is the category discrete?

A category is *discrete* if it is has no non-identity morphisms.
"""
is_discrete(C::FinCat) = isempty(hom_generators(C))

""" Is the category freely generated?
"""
is_free(C::FinCat) = isempty(equations(C))

# Opposite FinCats
#-----------------

const OppositeFinCat{Ob,Hom} = OppositeCat{Ob,Hom,FinCatSize}

ob_generators(C::OppositeFinCat) = ob_generators(C.cat)
hom_generators(C::OppositeFinCat) = hom_generators(C.cat)

ob_generator(C::OppositeCat, x) = ob_generator(C.cat, x)
hom_generator(C::OppositeCat, f) = hom_generator(C.cat, f)

# Categories on graphs
######################

""" Abstract type for category backed by finite generating graph.
"""
abstract type FinCatGraph{G,Ob,Hom} <: FinCat{Ob,Hom} end

""" Generating graph for a finitely presented category.
"""
graph(C::FinCatGraph) = C.graph

ob_generators(C::FinCatGraph) = vertices(graph(C))
hom_generators(C::FinCatGraph) = edges(graph(C))

ob_generator(C::FinCatGraph, x) = all(has_vertex(graph(C), x)) ? x :
  error("Vertex $x not contained in graph $(graph(C))")
hom_generator(C::FinCatGraph, f) = all(has_edge(graph(C), f)) ? f :
  error("Edge $f not contained in graph $(graph(C))")

ob_generator(C::FinCatGraph, x::Union{AbstractString,Symbol}) =
  vertex_named(graph(C), x)
hom_generator(C::FinCatGraph, f::Union{AbstractString,Symbol}) =
  edge_named(graph(C), f)
ob_generator_name(C::FinCatGraph, x) = vertex_name(graph(C), x)
hom_generator_name(C::FinCatGraph, f) = edge_name(graph(C), f)

# FIXME: These functions should go somewhere else, maybe in `Graphs`.
# Better yet, we should have a notion of "primary key" to avoid this.
vertex_name(G::HasGraph, v) = v
edge_name(G::HasGraph, e) = e
function vertex_named end
function edge_named end

function Base.show(io::IO, C::FinCatGraph)
  print(io, "FinCat(")
  show(io, graph(C))
  print(io, ", [")
  join(io, equations(C), ", ")
  print(io, "])")
end

# Paths in graphs
#----------------

""" Path in a graph.

The path is allowed to be empty but always has definite start and end points
(source and target vertices).
"""
@auto_hash_equals struct Path{V,E,Edges<:AbstractVector{E}}
  edges::Edges
  src::V
  tgt::V
end
edges(path::Path) = path.edges
src(path::Path) = path.src
tgt(path::Path) = path.tgt

function Path(g::HasGraph, es::AbstractVector)
  !isempty(es) || error("Nonempty edge list needed for nontrivial path")
  all(e -> has_edge(g, e), es) || error("Path has edges not contained in graph")
  Path(es, src(g, first(es)), tgt(g, last(es)))
end

function Path(g::HasGraph, e)
  has_edge(g, e) || error("Edge $e not contained in graph $g")
  Path(SVector(e), src(g,e), tgt(g,e))
end

function Base.empty(::Type{Path}, g::HasGraph, v::T) where T
  has_vertex(g, v) || error("Vertex $v not contained in graph $g")
  Path(SVector{0,T}(), v, v)
end

Base.reverse(p::Path) = Path(reverse(edges(p)), tgt(p), src(p))

function Base.vcat(p1::Path, p2::Path)
  tgt(p1) == src(p2) ||
    error("Path start/end points do not match: $(tgt(p1)) != $(src(p2))")
  Path(vcat(edges(p1), edges(p2)), src(p1), tgt(p2))
end

function Base.show(io::IO, path::Path)
  print(io, "Path(")
  show(io, edges(path))
  print(io, ": $(src(path)) → $(tgt(path)))")
end

""" Abstract type for category whose morphisms are paths in a graph.

(Or equivalence classes of paths in a graph, but we compute with
"""
const FinCatPathGraph{G,V,E} = FinCatGraph{G,V,Path{V,E}}

dom(C::FinCatPathGraph, e) = src(graph(C), e)
dom(C::FinCatPathGraph, path::Path) = src(path)
codom(C::FinCatPathGraph, e) = tgt(graph(C), e)
codom(C::FinCatPathGraph, path::Path) = tgt(path)

id(C::FinCatPathGraph, x) = empty(Path, graph(C), x)
compose(C::FinCatPathGraph, fs...) =
  reduce(vcat, coerce_path(graph(C), f) for f in fs)

hom(C::FinCatPathGraph, f) = coerce_path(graph(C), f)

coerce_path(g::HasGraph, path::Path) = path
coerce_path(g::HasGraph, x) = Path(g, x)

# Free category on graph
#-----------------------

""" Free category generated by a finite graph.

The objects of the free category are vertices in the graph and the morphisms are
(possibly empty) paths.
"""
@auto_hash_equals struct FreeCatGraph{G<:HasGraph} <: FinCatPathGraph{G,Int,Int}
  graph::G
end

FinCatGraph(g::HasGraph) = FreeCatGraph(g)

is_free(::FreeCatGraph) = true

function Base.show(io::IO, C::FreeCatGraph)
  print(io, "FinCat(")
  show(io, graph(C))
  print(io, ")")
end

# Category on graph with equations
#---------------------------------

""" Category presented by a finite graph together with path equations.

The objects of the category are vertices in the graph and the morphisms are
paths, quotiented by the congruence relation generated by the path equations.
See (Spivak, 2014, *Category theory for the sciences*, §4.5).
"""
@auto_hash_equals struct FinCatGraphEq{G<:HasGraph,Eqs<:AbstractVector{<:Pair}} <:
    FinCatPathGraph{G,Int,Int}
  graph::G
  equations::Eqs
end

equations(C::FinCatGraphEq) = C.equations

function FinCatGraph(g::HasGraph, eqs::AbstractVector)
  eqs = map(eqs) do eq
    lhs, rhs = coerce_path(g, first(eq)), coerce_path(g, last(eq))
    (src(lhs) == src(rhs) && tgt(lhs) == tgt(rhs)) ||
      error("Paths $lhs and $rhs in equation do not have equal (co)domains")
    lhs => rhs
  end
  FinCatGraphEq(g, eqs)
end

# Symbolic categories
#####################

""" Category defined by a `Presentation` object.

The presentation type can, of course, be a category (`Theories.Category`). It
can also be a schema (`Theories.Schema`). In this case, the schema's objects and
attribute types are regarded as the category's objects and the schema's
morphisms, attributes, and attribute types as the category's morphisms (where
the attribute types are identity morphisms). When the schema is formalized as a
profunctor whose codomain category is discrete, this amounts to taking the
collage of the profunctor.
"""
@auto_hash_equals struct FinCatPresentation{T,Ob,Hom} <: FinCat{Ob,Hom}
  presentation::Presentation{T}
end

function FinCatPresentation(pres::Presentation{T}) where T
  S = pres.syntax
  FinCatPresentation{T,S.Ob,S.Hom}(pres)
end
function FinCatPresentation(pres::Presentation{Schema})
  S = pres.syntax
  Ob = Union{S.Ob, S.AttrType}
  Hom = Union{S.Hom, S.Attr, S.AttrType}
  FinCatPresentation{Schema,Ob,Hom}(pres)
end

presentation(C::FinCatPresentation) = C.presentation

ob_generators(C::FinCatPresentation) = generators(presentation(C), :Ob)
ob_generators(C::FinCatPresentation{Schema}) = let P = presentation(C)
  vcat(generators(P, :Ob), generators(P, :AttrType))
end

hom_generators(C::FinCatPresentation) = generators(presentation(C), :Hom)
hom_generators(C::FinCatPresentation{Schema}) = let P = presentation(C)
  vcat(generators(P, :Hom), generators(P, :Attr))
end

equations(C::FinCatPresentation) = equations(presentation(C))

ob_generator(C::FinCatPresentation, x) = ob(C, presentation(C)[x])
ob_generator(C::FinCatPresentation, x::GATExpr{:generator}) = ob(C, x)
ob_generator_name(C::FinCatPresentation, x::GATExpr{:generator}) = first(x)

hom_generator(C::FinCatPresentation, f) = hom(C, presentation(C)[f])
hom_generator(C::FinCatPresentation, f::GATExpr{:generator}) = hom(C, f)
hom_generator_name(C::FinCatPresentation, f::GATExpr{:generator}) = first(f)

ob(C::FinCatPresentation, x::GATExpr) =
  gat_typeof(x) == :Ob ? x : error("Expression $x is not an object")
ob(C::FinCatPresentation{Schema}, x::GATExpr) =
  gat_typeof(x) ∈ (:Ob, :AttrType) ? x :
    error("Expression $x is not an object or attribute type")

hom(C::FinCatPresentation, f) = hom_generator(C, f)
hom(C::FinCatPresentation, fs::AbstractVector) =
  mapreduce(f -> hom(C, f), compose, fs)
hom(C::FinCatPresentation, f::GATExpr) =
  gat_typeof(f) == :Hom ? f : error("Expression $f is not a morphism")
hom(C::FinCatPresentation{Schema}, f::GATExpr) =
  gat_typeof(f) ∈ (:Hom, :Attr, :AttrType) ? f :
    error("Expression $f is not a morphism or attribute")

id(C::FinCatPresentation{Schema}, x::AttrTypeExpr) = x
compose(C::FinCatPresentation{Schema}, f::AttrTypeExpr, g::AttrTypeExpr) =
  (f == g) ? f : error("Invalid composite of attribute type identities: $f != $g")

function Base.show(io::IO, C::FinCatPresentation)
  print(io, "FinCat(")
  show(io, presentation(C))
  print(io, ")")
end

# Functors
##########

""" A functor out of a finitely presented category.
"""
const FinDomFunctor{Dom<:FinCat,Codom<:Cat} = Functor{Dom,Codom}

FinDomFunctor(maps::NamedTuple{(:V,:E)}, dom::FinCatGraph, codom::Cat) =
  FinDomFunctor(maps.V, maps.E, dom, codom)

function FinDomFunctor(ob_map, dom::FinCat, codom::Cat{Ob,Hom}) where {Ob,Hom}
  is_discrete(dom) ||
    error("Morphism map omitted by domain category is not discrete: $dom")
  FinDomFunctor(ob_map, empty(ob_map, Hom), dom, codom)
end
FinDomFunctor(ob_map, ::Nothing, dom::FinCat, codom::Cat) =
  FinDomFunctor(ob_map, dom, codom)

function hom_map(F::FinDomFunctor, path::Path)
  C, D = dom(F), codom(F)
  path = decompose(C, path)
  mapreduce(e -> hom_map(F, e), (gs...) -> compose(D, gs...),
            edges(path), init=id(D, ob_map(F, src(path))))
end
decompose(C::FinCatGraph, path::Path) = path
decompose(C::OppositeCat, path::Path) = reverse(path)

function hom_map(F::FinDomFunctor, f::GATExpr{:compose})
  C, D = dom(F), codom(F)
  mapreduce(f -> hom_map(F, f), (gs...) -> compose(D, gs...), decompose(C, f))
end
decompose(C::FinCatPresentation, f::GATExpr{:compose}) = args(f)
decompose(C::OppositeCat, f::GATExpr{:compose}) = reverse(decompose(C.cat, f))

function hom_map(F::FinDomFunctor, f::GATExpr{:id})
  id(codom(F), ob_map(F, dom(f)))
end

(F::FinDomFunctor)(expr::ObExpr) = ob_map(F, expr)
(F::FinDomFunctor)(expr::HomExpr) = hom_map(F, expr)

Categories.show_type_constructor(io::IO, ::Type{<:FinDomFunctor}) =
  print(io, "FinDomFunctor")

""" Collect assignments of functor's object map as a vector.
"""
collect_ob(F::FinDomFunctor) = map(x -> ob_map(F, x), ob_generators(dom(F)))

""" Collect assignments of functor's morphism map as a vector.
"""
collect_hom(F::FinDomFunctor) = map(f -> hom_map(F, f), hom_generators(dom(F)))

""" Is the purported functor on a presented category functorial?

This function checks that functor preserves domains and codomains. When
`check_equations` is `true` (the default), it also checks that the functor
preserves all equations in the domain category. Note that in some cases this may
not be possible.

See also: [`is_natural`](@ref).
"""
function is_functorial(F::FinDomFunctor; check_equations::Bool=false)
  C, D = dom(F), codom(F)
  all(hom_generators(C)) do f
    g = hom_map(F, f)
    dom(D,g) == ob_map(F, dom(C,f)) && codom(D,g) == ob_map(F, codom(C,f))
  end || return false

  if check_equations
    all(equations(C)) do (lhs, rhs)
      is_hom_equal(D, hom_map(F, lhs), hom_map(F, rhs))
    end || return false
  end

  true
end

function Base.map(F::Functor{<:FinCat,<:TypeCat}, f_ob, f_hom)
  C = dom(F)
  FinDomFunctor(map(x -> f_ob(ob_map(F, x)), ob_generators(C)),
                map(f -> f_hom(hom_map(F, f)), hom_generators(C)), C)
end

""" A functor between finitely presented categories.
"""
const FinFunctor{Dom<:FinCat,Codom<:FinCat} = FinDomFunctor{Dom,Codom}

FinFunctor(maps, dom::FinCat, codom::FinCat) = FinDomFunctor(maps, dom, codom)
FinFunctor(ob_map, hom_map, dom::FinCat, codom::FinCat) =
  FinDomFunctor(ob_map, hom_map, dom, codom)
FinFunctor(ob_map, hom_map, dom::Presentation, codom::Presentation) =
  FinDomFunctor(ob_map, hom_map, FinCat(dom), FinCat(codom))

Categories.show_type_constructor(io::IO, ::Type{<:FinFunctor}) =
  print(io, "FinFunctor")

# Mapping-based functors
#-----------------------

""" Functor out of a finitely presented category given by explicit mappings.
"""
@auto_hash_equals struct FinDomFunctorMap{Dom<:FinCat, Codom<:Cat,
    ObMap, HomMap} <: FinDomFunctor{Dom,Codom}
  ob_map::ObMap
  hom_map::HomMap
  dom::Dom
  codom::Codom
end

function FinDomFunctor(ob_map::Union{AbstractVector,AbstractDict},
                       hom_map::Union{AbstractVector,AbstractDict},
                       dom::FinCat, codom::Cat)
  length(ob_map) == length(ob_generators(dom)) ||
    error("Length of object map $ob_map does not match domain $dom")
  length(hom_map) == length(hom_generators(dom)) ||
    error("Length of morphism map $hom_map does not match domain $dom")

  ob_map = mappairs(x -> ob_key(dom, x), y -> ob(codom, y), ob_map)
  hom_map = mappairs(f -> hom_key(dom, f), g -> hom(codom, g), hom_map)
  FinDomFunctorMap(ob_map, hom_map, dom, codom)
end
FinDomFunctor(ob_map::Union{AbstractVector{Ob},AbstractDict{<:Any,Ob}},
              hom_map::Union{AbstractVector{Hom},AbstractDict{<:Any,Hom}},
              dom::FinCat) where {Ob,Hom} =
  FinDomFunctor(ob_map, hom_map, dom, TypeCat(Ob, Hom))

ob_key(C::FinCat, x) = ob_generator(C, x)
hom_key(C::FinCat, f) = hom_generator(C, f)

# Use generator names, rather than generators themselves, for Dict keys.
ob_key(C::FinCatPresentation, x) = presentation_key(x)
hom_key(C::FinCatPresentation, f) = presentation_key(f)
presentation_key(name::Union{AbstractString,Symbol}) = name
presentation_key(expr::GATExpr{:generator}) = first(expr)

Categories.do_ob_map(F::FinDomFunctorMap, x) = F.ob_map[ob_key(F.dom, x)]
Categories.do_hom_map(F::FinDomFunctorMap, f) = F.hom_map[hom_key(F.dom, f)]

collect_ob(F::FinDomFunctorMap) = values(F.ob_map)
collect_hom(F::FinDomFunctorMap) = values(F.hom_map)

op(F::FinDomFunctorMap) = FinDomFunctorMap(F.ob_map, F.hom_map,
                                           op(dom(F)), op(codom(F)))

""" Force evaluation of lazily defined function or functor.
"""
function force(F::FinDomFunctor)
  C = dom(F)
  FinDomFunctor(make_map(x -> ob_map(F, x), ob_generators(C)),
                make_map(f -> hom_map(F, f), hom_generators(C)), C)
end
force(F::FinDomFunctorMap) = F

function Categories.do_compose(F::FinDomFunctorMap, G::FinDomFunctorMap)
  FinDomFunctorMap(mapvals(x -> ob_map(G, x), F.ob_map),
                   mapvals(f -> hom_map(G, f), F.hom_map), dom(F), codom(G))
end

function Base.show(io::IO, F::T) where T <: FinDomFunctorMap
  Categories.show_type_constructor(io, T); print(io, "(")
  show(io, F.ob_map)
  print(io, ", ")
  show(io, F.hom_map)
  print(io, ", ")
  Categories.show_domains(io, F)
  print(io, ")")
end

# Natural transformations
#########################

""" A natural transformation whose domain category is finitely generated.

This type is for natural transformations ``α: F ⇒ G: C → D`` such that the
domain category ``C`` is finitely generated. Such a natural transformation is
given by a finite amount of data (one morphism in ``D`` for each generating
object of ``C``) and its naturality is verified by finitely many equations (one
equation for each generating morphism of ``C``).
"""
const FinTransformation{C<:FinCat,D<:Cat,Dom<:FinDomFunctor,Codom<:FinDomFunctor} =
  Transformation{C,D,Dom,Codom}

FinTransformation(F, G; components...) = FinTransformation(components, F, G)

""" Components of a natural transformation.
"""
components(α::FinTransformation) =
  make_map(x -> component(α, x), ob_generators(dom_ob(α)))

""" Is the transformation between `FinDomFunctors` a natural transformation?

This function uses the fact that to check whether a transformation is natural,
it suffices to check the naturality equations on a generating set of morphisms
of the domain category. In some cases, checking the equations may be expensive
or impossible. When the keyword argument `check_equations` is `false`, only the
domains and codomains of the components are checked.

See also: [`is_functorial`](@ref).
"""
function is_natural(α::FinTransformation; check_equations::Bool=true)
  F, G = dom(α), codom(α)
  C, D = dom(F), codom(F) # == dom(G), codom(G)
  all(ob_generators(C)) do c
    α_c = α[c]
    dom(D, α_c) == ob_map(F,c) && codom(D, α_c) == ob_map(G,c)
  end || return false

  if check_equations
    all(hom_generators(C)) do f
      Ff, Gf = hom_map(F,f), hom_map(G,f)
      α_c, α_d = α[dom(C,f)], α[codom(C,f)]
      is_hom_equal(D, compose(D, α_c, Gf), compose(D, Ff, α_d))
    end || return false
  end

  true
end

function check_transformation_domains(F::Functor, G::Functor)
  # XXX: Equality of TypeCats is too strict, so for now we are punting on
  # (co)domain checks in that case.
  C, C′ = dom(F), dom(G)
  (C isa TypeCat && C′ isa TypeCat) || C == C′ ||
    error("Mismatched domains in functors $F and $G")
  D, D′ = codom(F), codom(G)
  (D isa TypeCat && D′ isa TypeCat) || D == D′ ||
    error("Mismatched codomains in functors $F and $G")
  (C, D)
end

# Mapping-based transformations
#------------------------------

""" Natural transformation with components given by explicit mapping.
"""
@auto_hash_equals struct FinTransformationMap{C<:FinCat,D<:Cat,
    Dom<:FinDomFunctor{C,D},Codom<:FinDomFunctor,Comp} <: FinTransformation{C,D,Dom,Codom}
  components::Comp
  dom::Dom
  codom::Codom
end

function FinTransformation(components::Union{AbstractVector,AbstractDict},
                           F::FinDomFunctor, G::FinDomFunctor)
  C, D = check_transformation_domains(F, G)
  length(components) == length(ob_generators(C)) ||
    error("Incorrect number of components in $components for domain category $C")
  components = mappairs(x -> ob_key(C,x), f -> hom(D,f), components)
  FinTransformationMap(components, F, G)
end

component(α::FinTransformationMap, x) = α.components[ob_key(dom_ob(α), x)]
components(α::FinTransformationMap) = α.components

op(α::FinTransformationMap) = FinTransformationMap(components(α),
                                                   op(codom(α)), op(dom(α)))

function Categories.do_compose(α::FinTransformationMap, β::FinTransformation)
  F = dom(α)
  D = codom(F)
  FinTransformationMap(mapvals(α.components, keys=true) do c, f
                         compose(D, f, component(β, c))
                       end, F, codom(β))
end

function Categories.do_composeH(F::FinDomFunctorMap, β::Transformation)
  G, H = dom(β), codom(β)
  FinTransformationMap(mapvals(c -> component(β, c), F.ob_map),
                       compose(F, G, strict=false), compose(F, H, strict=false))
end

function Categories.do_composeH(α::FinTransformationMap, H::Functor)
  F, G = dom(α), codom(α)
  FinTransformationMap(mapvals(f -> hom_map(H, f), α.components),
                       compose(F, H, strict=false), compose(G, H, strict=false))
end

function Base.show(io::IO, α::FinTransformationMap)
  print(io, "FinTransformation(")
  show(io, components(α))
  print(io, ", ")
  Categories.show_domains(io, α, recurse=false)
  print(io, ", ")
  Categories.show_domains(io, dom(α))
  print(io, ")")
end

# Initial functors
##################

"""
Dual to a [final functor](https://ncatlab.org/nlab/show/final+functor), an
*initial functor* is one for which pulling back diagrams along it does not
change the limits of these diagrams.

This amounts to checking, for a functor C->D, that, for every object d in
Ob(D), the comma category (F/d) is connected.
"""
function is_initial(F::FinFunctor)::Bool
  Gₛ, Gₜ = graph(dom(F)), graph(codom(F))
  pathₛ, pathₜ = enumerate_paths.([Gₛ, Gₜ])

  function connected_nonempty_slice(t::Int)::Bool
    paths_into_t = incident(pathₜ, t, :tgt)
    # Generate slice objects
    ob_slice = Pair{Int,Vector{Int}}[] # s ∈ Ob(S) and a path ∈ T(F(s), t)
    for s in vertices(Gₛ)
      paths_s_to_t = incident(pathₜ, ob_map(F,s), :src) ∩ paths_into_t
      append!(ob_slice, [s => pathₜ[p, :eprops] for p in paths_s_to_t])
    end

    # Empty case
    isempty(ob_slice) && return false

    """
    For two slice objects (m,pₘ) and (n,pₙ) check for a morphism f ∈ S(M,N) such
    that there is a commutative triangle pₘ = f;pₙ
    """
    function check_pair(i::Int, j::Int)::Bool
      (m,pₘ), (n,pₙ) = ob_slice[i], ob_slice[j]
      es = incident(pathₛ, m, :src) ∩ incident(pathₛ, n, :tgt)
      paths = pathₛ[es, :eprops]
      return any(f -> pₘ == vcat(edges.(hom_map(F,f))..., pₙ), paths)
    end

    # Use check_pair to determine pairwise connectivity
    connected = IntDisjointSets(length(ob_slice)) # sym/trans/refl closure
    obs = 1:length(ob_slice)
    for (i,j) in Base.Iterators.product(obs, obs)
      if !in_same_set(connected, i, j) && check_pair(i,j)
        union!(connected, i, j)
      end
    end
    return num_groups(connected) == 1
  end

  # Check for each t ∈ T whether F/t is connected
  return all(connected_nonempty_slice, 1:nv(Gₜ))
end

# Dict utilities
################

# Something like this should be built into Julia, but unfortunately isn't.

function mappairs(kmap, vmap, pairs::T) where T<:AbstractDict
  (dicttype(T))(kmap(k) => vmap(v) for (k,v) in pairs)
end
mappairs(kmap, vmap, vec::AbstractVector) = map(vmap, vec)

function mapvals(f, pairs::T; keys::Bool=false, iter::Bool=false) where T<:AbstractDict
  (iter ? identity : dicttype(T))(if keys
    (k => f(k,v) for (k,v) in pairs)
  else
    (k => f(v) for (k,v) in pairs)
  end)
end
function mapvals(f, collection; keys::Bool=false, iter::Bool=false)
  do_map = iter ? Iterators.map : map
  if keys
    do_map(f, eachindex(collection), values(collection))
  else
    do_map(f, values(collection))
  end
end

dicttype(::Type{T}) where T <: AbstractDict = T.name.wrapper
dicttype(::Type{<:Iterators.Pairs}) = Dict

make_map(f, xs::UnitRange{Int}) = map(f, xs)
make_map(f, xs) = Dict(x => f(x) for x in xs)

end

max_forks_repo_name: epatters/CompCat


text:
using ApproxFun, MultivariateOrthogonalPolynomials, Plots, GLVisualize
    import MultivariateOrthogonalPolynomials: plan_evaluate

window = glscreen()
    @async GLWindow.waiting_renderloop(window)

S = WeightedTriangle(1,1,1)

Δ = Laplacian(S)
Dx = Derivative(S,[1,0])
Dy = Derivative(S,[0,1])

h=0.0001
ε=0.01

A=I-h*(ε*Δ+Dx)

QR=qrfact(A)
C=cache(eye(A))  # conversion operator

# Plotting grid
PT = ProductTriangle(1,1,1)
    X, Y = points(PT,30,30)
    X = Float32.(X)
    Y = Float32.(Y)

X = linspace(0f0,1f0,30)*ones(Float32,30)'
Y = X'
Y = (1-X).*Y

# Initial condition

u=u0=Fun(S,4randn(10));
    P = plan_evaluate(u)
    Z = Matrix{Float32}(P.(X,Y))
    vis = visualize((X,Y,Z), :surface)
    _view(vis)

u=u0=Fun(S,4randn(10));
    P = plan_evaluate(u)
    Z = Matrix{Float32}(P.(X,Y))
    vis = visualize((X,Y,Z), :surface)
    GLAbstraction.set_arg!(vis, :position_z, Z)
    yield()

Z

for k=1:500
    u=\(QR,C*u;tolerance=1E-5)
    chop!(u,1E-7)
    if mod(k,2) == 0
        P = plan_evaluate(u)
        Z = Matrix{Float32}(P.(X,Y))
        GLAbstraction.set_arg!(vis, :position_z, Z)
    end
    yield()
end


text:
using Configs
using Test

@testset "Unit tests" begin
    include("./unit.jl")
end
@testset "Functional tests" begin
    include("./functional.jl")
end

max_stars_repo_name: tisztamo/Configs

max_issues_repo_name: tisztamo/Configs

max_forks_repo_name: tisztamo/Configs


text:
"""
An abstraction that knows how to update a model and compute a search
direction (gradient estimate).
"""
@compat abstract type SearchDirection <: LearningStrategy end

type GradientAverager <: SearchDirection
    ∇avg::Vector{Float64}
    GradientAverager() = new()
end

function init(ga::GradientAverager, model)
    ga.∇avg = zeros(length(grad(model)))
end

# for a single observation, just return ∇
function search_direction(model, ga::GradientAverager, obs)
    update!(model, obs)
    grad(model)
end

# for a minibatch, compute the average gradient
function search_direction(model, ga::GradientAverager, batch::AbstractObsIterator)
    fill!(ga.∇avg, 0.0)
    scalar = 1 / nobs(batch)
    ∇ = grad(model)
    for obs in batch
        update!(model, obs)

        # add to the total param change for this gl/gradient
        @simd for i in 1:length(∇)
            @inbounds ga.∇avg[i] += ∇[i] * scalar
        end
    end
    ga.∇avg
end

# -------------------------------------------------------------

"""
A sub-learner which can update model parameters using a search direction (which might be an estimate
    of the gradient), with LearningRate lr and ParamUpdater pu (SGD, Adam, etc).

Note: we might update the model's internal state while computing the SearchDirection.
"""
immutable GradientLearner{LR <: LearningRate, PU <: ParamUpdater, SD <: SearchDirection} <: LearningStrategy
    lr::LR
    pu::PU
    sd::SD
end

function GradientLearner(lr::LearningRate = FixedLR(1e-1),
                         pu::ParamUpdater = RMSProp())
    GradientLearner(lr, pu, GradientAverager())
end
# function GradientLearner(pu::ParamUpdater,
#                          lr::LearningRate = FixedLR(1e-3),
#                          sd::SearchDirection = GradientAverager())
#     GradientLearner(lr, pu, sd)
# end
function GradientLearner(lr::Number,
                         pu::ParamUpdater = RMSProp(),
                         sd::SearchDirection = GradientAverager())
    GradientLearner(FixedLR(lr), pu, sd)
end

function pre_hook(gl::GradientLearner, model)
    init(gl.pu, model)
    init(gl.sd, model)
end

# one iteration update
function learn!(model, gl::GradientLearner, data)
    θ = params(model)
    ∇ = grad(model)

    # # optional setup before iteration update
    # before_grad_calc(θ, gl.pu, ∇)

    # get this iterations search direction (gradient estimate)
    sd = search_direction(model, gl.sd, data)

    # update the params using the search direction
    update!(θ, gl.pu, sd, value(gl.lr))
end

max_stars_repo_name: JuliaML/StochasticOptimization

max_issues_repo_name: JuliaML/StochasticOptimization

max_forks_repo_name: JuliaML/StochasticOptimization


text:
# NOTE:  backend should implement `html_body` and `html_head`

# CREDIT: parts of this implementation were inspired by @joshday's PlotlyLocal.jl

function standalone_html(plt::PlottingObject; title::AbstractString = get(plt.plotargs, :window_title, "Plots.jl"))
    """
    <!DOCTYPE html>
    <html>
        <head>
            <title>$title</title>
            $(html_head(plt))
        </head>
        <body>
            $(html_body(plt))
        </body>
    </html>
    """
end

function open_browser_window(filename::AbstractString)
    @osx_only   return run(`open $(filename)`)
    @linux_only return run(`xdg-open $(filename)`)
    @windows_only return run(`$(ENV["COMSPEC"]) /c start $(filename)`)
    warn("Unknown OS... cannot open browser window.")
end

function standalone_html_window(plt::PlottingObject; kw...)
    html = standalone_html(plt; kw...)
    # println(html)
    filename = string(tempname(), ".html")
    output = open(filename, "w")
    write(output, html)
    close(output)
    open_browser_window(filename)
end


text:
module CartesianDiscreteModelsTests

using Gridap
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Geometry: CartesianMap
using Gridap.Io
using Test

domain = (0,1,0,1,0,1)
partition = (3,3,3)
model = CartesianDiscreteModel(domain,partition)
test_discrete_model(model)
@test is_oriented(get_grid(model)) == true
@test is_oriented(get_grid_topology(model)) == true

grid_topology = get_grid_topology(model)
cell_to_pindices = get_cell_permutations(grid_topology)

labels = get_face_labeling(model)
@test num_tags(labels) == 28
@test num_entities(labels) == 27

map = get_cell_map(get_grid(model))
@test isa(map,CartesianMap)

desc = get_cartesian_descriptor(model)

function polar(q)
  r, t, z = q
  x = r*cos(t)
  y = r*sin(t)
  Point(x,y,z)
end

domain = (1,2,0,pi,0,0.5)
partition = (10,30,4)
model = CartesianDiscreteModel(domain,partition;map=polar)
test_discrete_model(model)
@test is_oriented(get_grid(model)) == true

model2 = from_dict(DiscreteModel,to_dict(model))
test_discrete_model(model2)

model3 = CartesianDiscreteModel(desc,CartesianIndex(2,2,2),CartesianIndex(3,3,3))
test_discrete_model(model3)

#using Gridap.Visualization
#writevtk(model2,"model2")

end # module


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCellDG
include("../useful_routines.jl")

polyorder = 1
numqp = 2

basis = LagrangeTensorProductBasis(2, polyorder)
x0 = [0.0, 0.0]
widths = [3.0, 1.0]
nelements = [3, 1]
dim, nodesperelement = size(interpolation_points(basis))
mesh = CutCellDG.CGMesh([0.0, 0.0], [3.0, 1.0], [3, 1], nodesperelement)
normal = [1.0, 0.0]
x0 = [0.5, 0.0]
levelset =
    CutCellDG.LevelSet(x -> plane_distance_function(x, normal, x0), mesh, basis)
nodalcoordinates = CutCellDG.nodal_coordinates(mesh)

tol = 1e-3
perturbation = 0.0

cellsign = CutCellDG.cell_sign!(levelset, tol, perturbation)
@test allequal(cellsign, [0, 1, 1])

normal = [1.0, 0.0]
x0 = [1.0, 0.0]
tol = 1e-3
perturbation = 1e-2
levelset =
    CutCellDG.LevelSet(x -> plane_distance_function(x, normal, x0), mesh, basis)
cellsign = CutCellDG.cell_sign!(levelset, tol, perturbation)

@test allapprox(
    CutCellDG.coefficients(levelset, 1),
    [-1.0 + perturbation, -1.0 + perturbation, perturbation, perturbation],
)
@test allapprox(
    CutCellDG.coefficients(levelset, 2),
    [perturbation, perturbation, 1.0, 1.0],
)
@test allapprox(CutCellDG.coefficients(levelset, 3), [1.0, 1.0, 2.0, 2.0])

basis = HermiteTensorProductBasis(2)
normal = [1.0, 0.0]
x0 = [1.0, 0.0]
tol = 1e-3
perturbation = 1e-2
levelset = CutCellDG.LevelSet(
    x -> plane_distance_function(x, normal, x0),
    mesh,
    basis,
)
cellsign = CutCellDG.cell_sign!(levelset,tol,perturbation)
testcellsign = [0,1,1]
@test allequal(testcellsign,cellsign)
coeffs1 = reshape(CutCellDG.coefficients(levelset,1),4,:)
jac = CutCellDG.jacobian(mesh)
testcoeffs1 = [-1+perturbation   -1+perturbation   perturbation   perturbation
                0.0               0.0              0.0            0.0
                jac[1]            jac[1]           jac[1]         jac[1]
                0.0               0.0              0.0            0.0]
@test allapprox(coeffs1,testcoeffs1,1e4eps())
coeffs2 = reshape(CutCellDG.coefficients(levelset,2),4,:)
testcoeffs2 = [perturbation    perturbation   1.0            1.0
               0.0             0.0            0.0            0.0
               jac[1]          jac[1]         jac[1]         jac[1]
               0.0             0.0            0.0            0.0]
@test allapprox(coeffs2,testcoeffs2,1e4eps())
coeffs3 = reshape(CutCellDG.coefficients(levelset,3),4,:)
testcoeffs3 = [1.0     1.0    2.0    2.0
               0.0     0.0    0.0    0.0
               jac[1]  jac[1] jac[1] jac[1]
               0.0     0.0    0.0    0.0]
@test allapprox(coeffs3,testcoeffs3,1e4eps())

basis = LagrangeTensorProductBasis(2,1)
refpoints = interpolation_points(basis)
dgmesh = CutCellDG.DGMesh(x0,widths,nelements,refpoints)
posactivenodeids = CutCellDG.active_node_ids(dgmesh,+1,cellsign)
@test allequal(posactivenodeids, 1:12)
negactivenodeids = CutCellDG.active_node_ids(dgmesh,-1, cellsign)
@test allequal(negactivenodeids, 1:4)

totalnumnodes = CutCellDG.number_of_nodes(dgmesh)
cutmeshnodeids =
    CutCellDG.cut_mesh_nodeids(posactivenodeids, negactivenodeids, totalnumnodes)
testcutmeshnodeids = [
    1 2 3 4 5 6 7 8 9 10 11 12
    13 14 15 16 0 0 0 0 0 0 0 0
]
@test allequal(cutmeshnodeids, testcutmeshnodeids)

cutmesh = CutCellDG.CutMesh(dgmesh, levelset)
@test allequal(CutCellDG.nodal_connectivity(cutmesh, +1, 1), [1, 2, 3, 4])
@test allequal(CutCellDG.nodal_connectivity(cutmesh, -1, 1), [13, 14, 15, 16])


text:
using Base.Test
using SymEngine
CDenseMatrix = SymEngine.CDenseMatrix

@vars x

# constructors
A = [x 1 2; 3 x 4; 5 6 x]
M = convert(CDenseMatrix, A)
M = CDenseMatrix(A)

# abstract interface
@test M[1,1] == A[1,1]
for i in eachindex(M)
    @test M[i] == A[i]
end
M[1,1] = x^2
@test M[1,1] == x^2
M[1,1] = x

# inherits from DenseArray
@test all(M + 3 .== A + 3)
@test all(M * 3 .== A * 3)
@test all(M + M .== A + A)
@test all(M * M .== A * A)
@test all(M' .== A')

# generic det
@test prod([subs(det(M) - det(A), x, i) == 0 for i in 2:10]) == true
@test inv(M) - inv(A) == zeros(Basic, 3,3)

# factorizations
@test lu(M) == lu(A)

A = [x 1 2; 0 x 4; 0 0 x]
b = [1, 2, 3]
M = convert(CDenseMatrix, A)
out = M \ b
@test M * out - b == zeros(Basic, 3, 1)


text:
using Test
using ExtensibleScheduler

@testset "Action" begin

    function return_noparam()
        "Returned from return_noparam"
    end

    function return_args(x)
        "Returned from return_args with $x"
    end

    function return_kwargs(; a="default")
        "Returned from return_kwargs with a=$a"
    end

    @testset "noparam" begin
        action = Action(return_noparam)
        expected = "Returned from return_noparam"
        @test run(action) == expected

        action = Action(return_noparam, ()...; Dict()...)
        @test run(action) == expected
    end

    @testset "args" begin
        action = Action(return_args, 3)
        expected = "Returned from return_args with 3"
        @test run(action) == expected
    end

    @testset "kwargs" begin
        action = Action(return_kwargs; Dict(:a=>5)...)
        expected = "Returned from return_kwargs with a=5"
        @test run(action) == expected
    end

end