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

module Doctrine
export Monocl, MonoclCategory, MonoclError, Ob, Hom, SubOb, SubHom,
  dom, codom, subob_dom, subob_codom, subob_id, subhom_id,
  compose, compose2, id, otimes, munit, opow, braid, mcopy, delete, pair,
  hom, ev, curry, coerce, construct, MonoclElem, to_wiring_diagram

using AutoHashEquals

using Catlab
import Catlab.Theories: CategoryExpr, ObExpr, HomExpr,
  SymmetricMonoidalCategory, Ob, Hom, dom, codom, compose, id, otimes, munit,
  braid, mcopy, delete, pair, hom, ev, curry

using Catlab.WiringDiagrams
import Catlab.WiringDiagrams: Box, WiringDiagram, to_wiring_diagram

# Cartesian (closed) category
#############################

""" Doctrine of *cartesian category*

This signature differs from the official Catlab doctrine by allowing `mcopy`
terms of size greater than 2.
"""
@signature CartesianCategory{Ob,Hom} <: SymmetricMonoidalCategory{Ob,Hom} begin
  opow(A::Ob, n::Int)::Ob

  mcopy(A::Ob, n::Int)::Hom(A,opow(A,n))
  delete(A::Ob)::Hom(A,munit())
end

""" Doctrine of *cartesian closed category*

This signature is identical to the official Catlab doctrine, except for
inheriting from a different doctrine.
"""
@signature CartesianClosedCategory{Ob,Hom} <: CartesianCategory{Ob,Hom} begin
  hom(A::Ob, B::Ob)::Ob
  ev(A::Ob, B::Ob)::Hom(otimes(hom(A,B),A),B)
  curry(A::Ob, B::Ob, f::Hom(otimes(A,B),C))::Hom(A,hom(B,C)) ⊣ (C::Ob)
end

# Monocl category
#################

struct MonoclError <: Exception
  message::String
end

""" Doctrine for Monocl: MONoidal Ontology and Computer Language

A doctrine of monoidal categories derived from the doctrine of cartesian closed
categories with implicit conversion of types.
"""
@signature MonoclCategory{Ob,Hom,SubOb,SubHom} <: CartesianClosedCategory{Ob,Hom} begin
  """ Subobject relation.

  The domain object is a subobject of the codomain object. Alternatively,
  in PLT jargon, the domain type is a subtype of the codomain type.
  """
  SubOb(dom::Ob, codom::Ob)::TYPE

  """ Submorphism relation.

  The domain morphism is a submorphism of the codomain morphism. This means
  there is a natural transformation, whose components are subobject morphisms,
  having as (co)domain the functor on the interval category determined by the
  (co)domain morphism.

  Alternatively, in PLT jargon, the domain function is a specialization of the
  codomain function, a form of ad-hoc polymorphism.
  """
  SubHom(dom::Hom(A0,B0), codom::Hom(A,B),
         subob_dom::SubOb(A0,A), subob_codom::SubOb(B0,B))::TYPE ⊣
    (A0::Ob, B0::Ob, A::Ob, B::Ob)

  """ Coercion morphism of type A to type B.
  """
  coerce(sub::SubOb(A,B))::Hom(A,B) ⊣ (A::Ob, B::Ob)

  """ Constructor for instances of type A with data f: A -> B.

  The semantics of this term are:
     compose(construct(f), f) = id(B)
  """
  construct(f::Hom(A,B))::Hom(B,A) ⊣ (A::Ob, B::Ob)

  # Category of subobject morphisms.
  # TODO: Support internal homs.
  # XXX: Cannot reuse `id` for subobjects because cannot dispatch on return type.
  subob_id(A::Ob)::SubOb(A,A)
  compose(f::SubOb(A,B), g::SubOb(B,C))::SubOb(A,C) ⊣ (A::Ob, B::Ob, C::Ob)
  otimes(f::SubOb(A,B), g::SubOb(C,D))::SubOb(otimes(A,C),otimes(B,D)) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob)

  # 2-category of submorphism natural transformations.
  subhom_id(f::Hom(A,B))::SubHom(f,f,subob_id(dom(f)),subob_id(codom(f))) ⊣
    (A::Ob, B::Ob)
  compose(α::SubHom(f,g,α0,α1), β::SubHom(g,h,β0,β1))::SubHom(f,h,compose(α0,β0),compose(α1,β1)) ⊣
    (A0::Ob, B0::Ob, A::Ob, B::Ob, A1::Ob, B1::Ob,
     f::Hom(A0,B0), g::Hom(A,B), h::Hom(A1,B1),
     α0::SubOb(A0,A), α1::SubOb(B0,B), β0::SubOb(A,A1), β1::SubOb(B,B1))
  compose2(α::SubHom(f,g,α0,α1), β::SubHom(h,k,β0,β1))::SubHom(compose(f,h),compose(g,k),α0,β1) ⊣
    (A0::Ob, B0::Ob, C0::Ob, A::Ob, B::Ob, C::Ob,
     f::Hom(A0,B0), g::Hom(A,B), h::Hom(B0,C0), k::Hom(B,C),
     α0::SubOb(A0,A), α1::SubOb(B0,B), β0::SubOb(B0,B), β1::SubOb(C0,C))
end

""" Syntax system for Monocl: MONoidal Ontology and Computer Language
"""
@syntax Monocl{ObExpr,HomExpr,CategoryExpr,CategoryExpr} MonoclCategory begin
  # TODO: Implicit conversion is not yet implemented, so we have disabled
  # domain checks when composing morphisms!
  # TODO: Domain checks when composing submorphisms need only check domain
  # and codomain of subobjects because subobjects form a pre-order.
  compose(f::Hom, g::Hom) = associate_unit(new(f,g; strict=false), id)
  compose(f::SubOb, g::SubOb) = associate_unit(new(f,g; strict=true), subob_id)
  compose(f::SubHom, g::SubHom) = associate_unit(new(f,g; strict=false), subhom_id)
  compose2(f::SubHom, g::SubHom) = associate_unit(new(f,g; strict=false), subhom_id)

  otimes(A::Ob, B::Ob) = associate_unit(new(A,B), munit)
  otimes(f::Hom, g::Hom) = associate(new(f,g))
  otimes(f::SubOb, g::SubOb) = associate(new(f,g))
  opow(A::Ob, n::Int) = otimes([A for i=1:n])

  # TODO: Enforce pre-order, not just reflexivity.
  coerce(sub::SubOb) = dom(sub) == codom(sub) ? id(dom(sub)) : new(sub)
end

SubOb(dom::Monocl.Ob, codom::Monocl.Ob) = SubOb(nothing, dom, codom)
SubHom(dom::Monocl.Hom, codom::Monocl.Hom,
       subob_dom::Monocl.SubOb, subob_codom::Monocl.SubOb) =
  SubHom(nothing, dom, codom, subob_dom, subob_codom)

""" Pairing of two (or more) morphisms.

Pairing is possible in any cartesian category. This method differs from the
standard Catlab definition by allowing coercion of the common domain object.
"""
function pair(A::Monocl.Ob, fs::Vector{Monocl.Hom})
  compose(mcopy(A,length(fs)), otimes(fs))
end
function pair(fs::Vector{Monocl.Hom})
  A = dom(first(fs))
  @assert all(dom(f) == A for f in fs)
  pair(A, fs)
end
pair(A::Monocl.Ob, fs::Vararg{Monocl.Hom}) = pair(A, collect(Monocl.Hom,fs))
pair(fs::Vararg{Monocl.Hom}) = pair(collect(Monocl.Hom,fs))

mcopy(A::Monocl.Ob) = mcopy(A,2)
construct(A::Monocl.Ob) = construct(delete(A))

# Monocl category of elements
#############################

""" Object in Monocl's category of elements.
"""
@auto_hash_equals struct MonoclElem
  ob::Union{Monocl.Ob,Nothing}
  value::Any
  MonoclElem(ob, value=nothing) = new(ob, value)
end

# Monocl wiring diagrams
########################

function Box(f::Monocl.Hom)
  Box(f, collect(dom(f)), collect(codom(f)))
end
function WiringDiagram(dom::Monocl.Ob, codom::Monocl.Ob)
  WiringDiagram(collect(dom), collect(codom))
end

function to_wiring_diagram(expr::Monocl.Hom)
  functor((Ports, WiringDiagram, Monocl.SubOb), expr;
    terms = Dict(
      :Ob => (expr) -> Ports([expr]),
      :Hom => (expr) -> singleton_diagram(Box(expr)),
      :SubOb => identity,
      :coerce => (expr) -> to_wiring_diagram(first(expr)),
      :construct => (expr) -> singleton_diagram(Box(expr)),
    )
  )
end

function to_wiring_diagram(sub::Monocl.SubOb)
  A, B = collect(dom(sub)), collect(codom(sub))
  @assert length(A) == length(B)
  f = WiringDiagram(A, B)
  add_wires!(f, ((input_id(f),i) => (output_id(f),i) for i in eachindex(A)))
  return f
end

end

max_stars_repo_name: IBM/semanticflowgraph

max_issues_repo_name: epatters/semanticflowgraph

max_forks_repo_name: IBM/semanticflowgraph


text:
module ExogenousStates

using Parameters
import QuantEcon: MarkovChain, stationary_distributions

abstract type StateSpace end

struct EndogenousStateSpace{T0,T1,T2,T3,T4} <: StateSpace
  grids::T0
  grid::T1 # grid
  indices::T2
  linear_indices::T3
  size::T4
end

struct ExogenousStateSpace{T0,T1,T2,T3,T4,T5,T6} <: StateSpace
  grids::T0
  grid::T1 # grid
  indices::T2
  linear_indices::T3
  size::T4
  mc::T5
  dist::T6 # stationary distribution
end

Base.size(ss::StateSpace) = ss.size
Base.length(ss::StateSpace) = prod(size(ss))
Base.keys(ss::StateSpace) = keys(ss.grid[1])
MarkovChain(exo::ExogenousStateSpace) = exo.mc
dimension(ss::StateSpace) = length(size(ss))
linear_indices(ss::StateSpace) = ss.linear_indices

function EndogenousStateSpace(grids_nt::NamedTuple)
  keys_ = keys(grids_nt)

  size_ = Tuple(length.(collect(grids_nt)))

  grid0 = Iterators.product(grids_nt...)
  gridNT = NamedTuple{keys_}.(grid0)

  indices0 = Iterators.product([1:n for n in size_]...)
  indicesNT = NamedTuple{keys_}.(indices0)

  EndogenousStateSpace(grids_nt, gridNT, indicesNT, LinearIndices(size_), size_)

end

using Test

@testset "Endogenous state space" begin
  a_min = 0.0
  h_min = eps()

  a = LinRange(a_min, 7.5, 10)
  h = LinRange(h_min, 5.0, 15)

  endo = EndogenousStateSpace((a=a, h=h))

  @test size(endo) == (10, 15)
  @test length(endo) == 10 * 15
  @test keys(endo) == (:a, :h)
  @test endo.grid[1,1] == (a=a_min, h=h_min)
  @test endo.grid[5,7] == (a=a[5], h=h[7])
  @test endo.indices[3,14] == (a=3, h=14)
end

function ExogenousStateSpace(vec_mc)
  grids = getproperty.(vec_mc, :state_values)
  size = Tuple(length.(grids))

  mc = product(vec_mc...)

  grid = mc.state_values

  indices0 = collect(Iterators.product([1:n for n in size]...))
  indicesNT = NamedTuple{keys(grid[1])}.(indices0)

  ExogenousStateSpace(grids, grid, indicesNT, LinearIndices(size), size, mc, stationary_distributions(mc)[1])
end

named_grid(grid, name) = NamedTuple{(name,)}.(Tuple.(Ref.(grid)))

function add_name(mc::MarkovChain, name::Symbol)
  MarkovChain(mc.p, named_grid(mc.state_values, name))
end

function MarkovChain(p, state_values, name::Symbol)
  MarkovChain(p, named_grid(state_values, name))
end

## Vector of independent Markov chains x1, x2

function product(mc1::MarkovChain, mc2::MarkovChain)

  combined_grid = [merge(s1, s2) for s2 in mc2.state_values for s1 in mc1.state_values]

  MarkovChain(kron(mc2.p, mc1.p), combined_grid)
end

function product(mc_vec...)
  #@show eltype(mc_vec)
  @assert eltype(mc_vec) <: MarkovChain
  if length(mc_vec) == 1
     return mc_vec[1]
  else
    mc12 = product(mc_vec[1], mc_vec[2])
  end
  if length(mc_vec) == 2
    return mc12
  else
    return product(mc12, mc_vec[3], mc_vec[4:end]...)
  end
end

export EndogenousStateSpace, ExogenousStateSpace, dimension, linear_indices

end

ext: jl


text:
"""
    MRAxis(f0, dt, numpoints; ppm0=4.7)
    MRAxis(experiment, numpoints; ppm0=4.7)

A combined time and frequency axis object for an MR experiment.

`f0` is the spectrometer frequency in MHz.  `dt` is the time between ADC
samples in seconds (somewhat bizarrely called "dwell time" in Siemens
terminology).

The number of points in FID acquisition is given by `numpoints` - for now we
resist trying to guess this from `expt` until more is known about the
structure there.

`ppm0` is the nominal PPM of the spectrometer reference frequency which is set
to 4.7 by default.  That is, we assume the spectrometer is tuned to the water
resonance, and 4.7 is a chosen nominal chemical shift of water relative to TMS.
Note that this is only approximate because the true shift can vary based on pH,
temperature, etc and you should calibrate your spectral axes against known
metabolites.
"""
struct MRAxis
    f0::Float64    # MHz
    dt::Float64    # s
    n::Int         #
    ppm0::Float64  # ppm
end

MRAxis(f0, dt, numpoints; ppm0=4.7) = MRAxis(f0, dt, numpoints, ppm0)

function MRAxis(expt::MRExperiment, numpoints; ppm0=4.7)
    # The 1e-6 conversion here is so that we have f0 in MHz. (As a consequence,
    # the factor of 1e6 in unit conversion to ppm happens automatically.)
    # TODO: Argubly ugly, and perhaps should just be in Hz.
    f0 = 1e-6 * expt.metadata["sTXSPEC.asNucleusInfo[0].lFrequency"]
    # The Siemens name "DwellTime" is probably a stretched analogy with radar
    # signal processing terminology.
    dt = 1e-9 * expt.metadata["sRXSPEC.alDwellTime[0]"]
    MRAxis(f0, dt, numpoints, ppm0)
end

# TODO
# Base.show(io::IO, a::MRAxis)

hertz_to_ppm(a::MRAxis, f) = a.ppm0 - f/a.f0
ppm_to_hertz(a::MRAxis, ppm) = (a.ppm0 - ppm)*a.f0

time_axis(a::MRAxis) = a.dt*(0:a.n-1)

function frequency_axis(a::MRAxis)
    # Assume frequency axis is for the data `fftshift(fft(fid))`
    # Offset of zero frequency sample after fftshift is:
    zero_offset = floor(Int, a.n/2)
    (-zero_offset:a.n-1-zero_offset) / (a.n*a.dt)
end

function frequency_axis_ppm(a::MRAxis)
    @warn "Is it right that this flips the axis direction?"
    hertz_to_ppm.(Ref(a), frequency_axis(a))
end

function frequency_axis(a::AxisArray, axis=Axis{:time})
    t = AxisArrays.axes(a, axis).val
    n = length(t)
    zero_offset = floor(Int, n/2)
    (-zero_offset:n-1-zero_offset) * uconvert.(u"Hz", one(eltype(t))/(n*step(t)))
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
###############################################################################
#
#   Fraction.jl : generic fraction fields
#
###############################################################################

export FractionField, num, den

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

parent_type(::Type{Frac{T}}) where T <: RingElem = FracField{T}

elem_type(::Type{FracField{T}}) where {T <: RingElem} = Frac{T}

doc"""
    base_ring{T}(S::Nemo.FracField{T})
> Return the base ring $R$ of the given fraction field.
"""
base_ring(a::Nemo.FracField{T}) where T <: RingElem = a.base_ring::parent_type(T)

doc"""
    base_ring{T}(r::Nemo.FracElem)
> Return the base ring $R$ of the fraction field that the supplied
> element $a$ belongs to.
"""
base_ring(a::Nemo.FracElem) = base_ring(parent(a))

doc"""
    parent(a::Nemo.FracElem)
> Return the parent object of the given fraction element.
"""
parent(a::Nemo.FracElem) = a.parent

function isdomain_type(::Type{T}) where {S <: RingElement, T <: Nemo.FracElem{S}}
   return isdomain_type(S)
end

function isexact_type(a::Type{T}) where {S <: RingElement, T <: Nemo.FracElem{S}}
   return isexact_type(S)
end

doc"""
    characteristic{T <: RingElem}(R::Nemo.FracField{T})
> Return the characteristic of the given field.
"""
function characteristic(R::Nemo.FracField{T}) where T <: RingElem
   return characteristic(base_ring(R))
end

function check_parent(a::Nemo.FracElem, b::Nemo.FracElem)
   parent(a) != parent(b) && error("Incompatible rings in fraction field operation")
end

###############################################################################
#
#   Constructors
#
###############################################################################

function //(x::T, y::T) where {T <: RingElem}
   iszero(y) && throw(DivideError())
   g = gcd(x, y)
   z = Frac{T}(divexact(x, g), divexact(y, g))
   try
      z.parent = FracDict[R]
   catch
      z.parent = FractionField(parent(x))
   end
   return z
end

//(x::T, y::Nemo.FracElem{T}) where {T <: RingElem} = parent(y)(x)//y

//(x::Nemo.FracElem{T}, y::T) where {T <: RingElem} = x//parent(x)(y)

###############################################################################
#
#   Basic manipulation
#
###############################################################################

function Base.hash(a::Nemo.FracElem, h::UInt)
   b = 0x8a30b0d963237dd5%UInt
   return xor(b, hash(numerator(a), h), hash(denominator(a), h), h)
end

function numerator(a::Nemo.FracElem)
   u = canonical_unit(a.den)
   return divexact(a.num, u)
end

function denominator(a::Nemo.FracElem)
   u = canonical_unit(a.den)
   return divexact(a.den, u)
end

doc"""
    zero(R::Nemo.FracField)
> Return $0/1$ in the given fraction field.
"""
zero(R::Nemo.FracField) = R(0)

doc"""
    one(R::Nemo.FracField)
> Return $1/1$ in the given fraction field.
"""
one(R::Nemo.FracField) = R(1)

doc"""
    iszero(a::Nemo.FracElem)
> Return `true` if the supplied element $a$ is zero in the fraction field it
> belongs to, otherwise return `false`.
"""
iszero(a::Nemo.FracElem) = iszero(numerator(a))

doc"""
    isone(a::Nemo.FracElem)
> Return `true` if the supplied element $a$ is one in the fraction field it
> belongs to, otherwise return `false`.
"""
isone(a::Nemo.FracElem) = numerator(a) == denominator(a)

doc"""
    isunit(a::Nemo.FracElem)
> Return `true` if the supplied element $a$ is invertible in the fraction field
> it belongs to, i.e. the numerator is nonzero, otherwise return `false`.
"""
isunit(a::Nemo.FracElem) = !iszero(numerator(a))

function deepcopy_internal(a::Frac{T}, dict::ObjectIdDict) where {T <: RingElem}
   v = Frac{T}(deepcopy(numerator(a)), deepcopy(denominator(a)))
   v.parent = parent(a)
   return v
end

###############################################################################
#
#   Canonicalisation
#
###############################################################################

canonical_unit(a::Nemo.FracElem) = a

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, x::Nemo.FracElem)
   u = canonical_unit(denominator(x))
   n = divexact(numerator(x), u)
   d = divexact(denominator(x), u);
   if d != 1 && needs_parentheses(n)
      print(io, "(")
   end
   print(io, n)
   if d != 1
      if needs_parentheses(n)
         print(io, ")")
      end
      print(io, "//")
      if needs_parentheses(d)
         print(io, "(")
      end
      print(io, d)
      if needs_parentheses(d)
         print(io, ")")
      end
   end
end

function show(io::IO, a::Nemo.FracField)
   print(io, "Fraction field of ", base_ring(a))
end

needs_parentheses(x::Nemo.FracElem) = isone(denominator(x)) && needs_parentheses(numerator(x))

isnegative(x::Nemo.FracElem) = !needs_parentheses(numerator(x)) && isnegative(numerator(x))

show_minus_one(::Type{Nemo.FracElem{T}}) where {T <: RingElem} = show_minus_one(T)

###############################################################################
#
#   Unary operators
#
###############################################################################

doc"""
    -(a::Nemo.FracElem)
> Return $-a$.
"""
function -(a::Nemo.FracElem)
   return parent(a)(-numerator(a), denominator(a))
end

###############################################################################
#
#   Binary operators
#
###############################################################################

doc"""
    +{T <: RingElem}(a::Nemo.FracElem{T}, b::Nemo.FracElem{T})
> Return $a + b$.
"""
function +(a::Nemo.FracElem{T}, b::Nemo.FracElem{T}) where {T <: RingElem}
   check_parent(a, b)
   n = numerator(a)*denominator(b) + numerator(b)*denominator(a)
   d = denominator(a)*denominator(b)
   g = gcd(n, d)
   return parent(a)(divexact(n, g), divexact(d, g))
end

doc"""
    -{T <: RingElem}(a::Nemo.FracElem{T}, b::Nemo.FracElem{T})
> Return $a - b$.
"""
function -(a::Nemo.FracElem{T}, b::Nemo.FracElem{T}) where {T <: RingElem}
   check_parent(a, b)
   n = numerator(a)*denominator(b) - numerator(b)*denominator(a)
   d = denominator(a)*denominator(b)
   g = gcd(n, d)
   return parent(a)(divexact(n, g), divexact(d, g))
end

doc"""
    *{T <: RingElem}(a::Nemo.FracElem{T}, b::Nemo.FracElem{T})
> Return $a\times b$.
"""
function *(a::Nemo.FracElem{T}, b::Nemo.FracElem{T}) where {T <: RingElem}
   check_parent(a, b)
   g1 = gcd(numerator(a), denominator(b))
   g2 = gcd(numerator(b), denominator(a))
   n = divexact(numerator(a), g1)*divexact(numerator(b), g2)
   d = divexact(denominator(a), g2)*divexact(denominator(b), g1)
   return parent(a)(n, d)
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

doc"""
    *(a::Nemo.FracElem, b::Union{Integer, Rational, AbstractFloat})
> Return $a\times b$.
"""
function *(a::Nemo.FracElem, b::Union{Integer, Rational, AbstractFloat})
   c = base_ring(a)(b)
   g = gcd(denominator(a), c)
   n = numerator(a)*divexact(c, g)
   d = divexact(denominator(a), g)
   return parent(a)(n, d)
end

doc"""
    *(a::Union{Integer, Rational, AbstractFloat}, b::Nemo.FracElem)
> Return $a\times b$.
"""
function *(a::Union{Integer, Rational, AbstractFloat}, b::Nemo.FracElem)
   c = base_ring(b)(a)
   g = gcd(denominator(b), c)
   n = numerator(b)*divexact(c, g)
   d = divexact(denominator(b), g)
   return parent(b)(n, d)
end

doc"""
    *{T <: RingElem}(a::Nemo.FracElem{T}, b::T)
> Return $a\times b$.
"""
function *(a::Nemo.FracElem{T}, b::T) where {T <: RingElem}
   g = gcd(denominator(a), b)
   n = numerator(a)*divexact(b, g)
   d = divexact(denominator(a), g)
   return parent(a)(n, d)
end

doc"""
    *{T <: RingElem}(a::T, b::Nemo.FracElem{T})
> Return $a\times b$.
"""
function *(a::T, b::Nemo.FracElem{T}) where {T <: RingElem}
   g = gcd(denominator(b), a)
   n = numerator(b)*divexact(a, g)
   d = divexact(denominator(b), g)
   return parent(b)(n, d)
end

doc"""
    +(a::Nemo.FracElem, b::Union{Integer, Rational, AbstractFloat})
> Return $a + b$.
"""
function +(a::Nemo.FracElem, b::Union{Integer, Rational, AbstractFloat})
   n = numerator(a) + denominator(a)*b
   d = denominator(a)
   g = gcd(n, d)
   return parent(a)(divexact(n, g), divexact(d, g))
end

doc"""
    -(a::Nemo.FracElem, b::Union{Integer, Rational, AbstractFloat})
> Return $a - b$.
"""
function -(a::Nemo.FracElem, b::Union{Integer, Rational, AbstractFloat})
   n = numerator(a) - denominator(a)*b
   d = denominator(a)
   g = gcd(n, d)
   return parent(a)(divexact(n, g), divexact(d, g))
end

doc"""
    +(a::Union{Integer, Rational, AbstractFloat}, b::Nemo.FracElem)
> Return $a + b$.
"""
+(a::Union{Integer, Rational, AbstractFloat}, b::Nemo.FracElem) = b + a

doc"""
    -(a::Union{Integer, Rational, AbstractFloat}, b::Nemo.FracElem)
> Return $a - b$.
"""
function -(a::Union{Integer, Rational, AbstractFloat}, b::Nemo.FracElem)
   n = a*denominator(b) - numerator(b)
   d = denominator(b)
   g = gcd(n, d)
   return parent(b)(divexact(n, g), divexact(d, g))
end

doc"""
    +{T <: RingElem}(a::Nemo.FracElem{T}, b::T)
> Return $a + b$.
"""
function +(a::Nemo.FracElem{T}, b::T) where {T <: RingElem}
   n = numerator(a) + denominator(a)*b
   d = denominator(a)
   g = gcd(n, d)
   return parent(a)(divexact(n, g), divexact(d, g))
end

doc"""
    -{T <: RingElem}(a::Nemo.FracElem{T}, b::T)
> Return $a - b$.
"""
function -(a::Nemo.FracElem{T}, b::T) where {T <: RingElem}
   n = numerator(a) - denominator(a)*b
   d = denominator(a)
   g = gcd(n, d)
   return parent(a)(divexact(n, g), divexact(d, g))
end

doc"""
    +{T <: RingElem}(a::T, b::Nemo.FracElem{T})
> Return $a + b$.
"""
+(a::T, b::Nemo.FracElem{T}) where {T <: RingElem} = b + a

doc"""
    -{T <: RingElem}(a::T, b::Nemo.FracElem{T})
> Return $a - b$.
"""
function -(a::T, b::Nemo.FracElem{T}) where {T <: RingElem}
   n = a*denominator(b) - numerator(b)
   d = denominator(b)
   g = gcd(n, d)
   return parent(b)(divexact(n, g), divexact(d, g))
end

###############################################################################
#
#   Comparisons
#
###############################################################################

doc"""
    =={T <: RingElem}(x::Nemo.FracElem{T}, y::Nemo.FracElem{T})
> Return `true` if $x == y$ arithmetically, otherwise return `false`. Recall
> that power series to different precisions may still be arithmetically
> equal to the minimum of the two precisions.
"""
function ==(x::Nemo.FracElem{T}, y::Nemo.FracElem{T}) where {T <: RingElem}
   check_parent(x, y)
   return (denominator(x) == denominator(y) && numerator(x) == numerator(y)) || (numerator(x)*denominator(y) == denominator(x)*numerator(y))
end

doc"""
    isequal{T <: RingElem}(x::Nemo.FracElem{T}, y::Nemo.FracElem{T})
> Return `true` if $x == y$ exactly, otherwise return `false`. This function is
> useful in cases where the numerators and denominators of the fractions are
> inexact, e.g. power series. Only if the power series are precisely the same,
> to the same precision, are they declared equal by this function.
"""
function isequal(x::Nemo.FracElem{T}, y::Nemo.FracElem{T}) where {T <: RingElem}
   if parent(x) != parent(y)
      return false
   end
   return isequal(numerator(x)*denominator(y), denominator(x)*numerator(y))
end

###############################################################################
#
#   Ad hoc comparisons
#
###############################################################################

doc"""
    ==(x::Nemo.FracElem, y::Union{Integer, Rational, AbstractFloat})
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
function ==(x::Nemo.FracElem, y::Union{Integer, Rational, AbstractFloat})
   return (isone(denominator(x)) && numerator(x) == y) || (numerator(x) == denominator(x)*y)
end

doc"""
    ==(x::Union{Integer, Rational, AbstractFloat}, y::Nemo.FracElem)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::Union{Integer, Rational, AbstractFloat}, y::Nemo.FracElem) = y == x

doc"""
    =={T <: RingElem}(x::Nemo.FracElem{T}, y::T)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
function ==(x::Nemo.FracElem{T}, y::T) where {T <: RingElem}
   return (isone(denominator(x)) && numerator(x) == y) || (numerator(x) == denominator(x)*y)
end

doc"""
    =={T <: RingElem}(x::T, y::Nemo.FracElem{T})
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::T, y::Nemo.FracElem{T}) where {T <: RingElem} = y == x

###############################################################################
#
#   Inversion
#
###############################################################################

doc"""
    inv(a::Nemo.FracElem)
> Return the inverse of the fraction $a$.
"""
function inv(a::Nemo.FracElem)
   iszero(numerator(a)) && throw(DivideError())
   return parent(a)(denominator(a), numerator(a))
end

###############################################################################
#
#   Exact division
#
###############################################################################

doc"""
    divexact{T <: RingElem}(a::Nemo.FracElem{T}, b::Nemo.FracElem{T})
> Return $a/b$.
"""
function divexact(a::Nemo.FracElem{T}, b::Nemo.FracElem{T}) where {T <: RingElem}
   check_parent(a, b)
   g1 = gcd(numerator(a), numerator(b))
   g2 = gcd(denominator(b), denominator(a))
   n = divexact(numerator(a), g1)*divexact(denominator(b), g2)
   d = divexact(denominator(a), g2)*divexact(numerator(b), g1)
   return parent(a)(n, d)
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

doc"""
    divexact(a::Nemo.FracElem, b::Union{Integer, Rational, AbstractFloat})
> Return $a/b$.
"""
function divexact(a::Nemo.FracElem, b::Union{Integer, Rational, AbstractFloat})
   b == 0 && throw(DivideError())
   c = base_ring(a)(b)
   g = gcd(numerator(a), c)
   n = divexact(numerator(a), g)
   d = denominator(a)*divexact(c, g)
   return parent(a)(n, d)
end

doc"""
    divexact(a::Union{Integer, Rational, AbstractFloat}, b::Nemo.FracElem)
> Return $a/b$.
"""
function divexact(a::Union{Integer, Rational, AbstractFloat}, b::Nemo.FracElem)
   iszero(b) && throw(DivideError())
   c = base_ring(b)(a)
   g = gcd(numerator(b), c)
   n = denominator(b)*divexact(c, g)
   d = divexact(numerator(b), g)
   return parent(b)(n, d)
end

doc"""
    divexact{T <: RingElem}(a::Nemo.FracElem{T}, b::T)
> Return $a/b$.
"""
function divexact(a::Nemo.FracElem{T}, b::T) where {T <: RingElem}
   iszero(b) && throw(DivideError())
   g = gcd(numerator(a), b)
   n = divexact(numerator(a), g)
   d = denominator(a)*divexact(b, g)
   return parent(a)(n, d)
end

doc"""
    divexact{T <: RingElem}(a::T, b::Nemo.FracElem{T})
> Return $a/b$.
"""
function divexact(a::T, b::Nemo.FracElem{T}) where {T <: RingElem}
   iszero(b) && throw(DivideError())
   g = gcd(numerator(b), a)
   n = denominator(b)*divexact(a, g)
   d = divexact(numerator(b), g)
   return parent(b)(n, d)
end

function divides(a::Nemo.FracElem{T}, b::Nemo.FracElem{T}) where {T <: RingElem}
   iszero(b) && error("Division by zero in divides")
   return true, divexact(a, b)
end

###############################################################################
#
#   Powering
#
###############################################################################

doc"""
    ^(a::Nemo.FracElem, b::Int)
> Return $a^b$.
"""
function ^(a::Nemo.FracElem{T}, b::Int) where {T <: RingElem}
   if b < 0
      a = inv(a)
      b = -b
   end
   return parent(a)(numerator(a)^b, denominator(a)^b)
end

###############################################################################
#
#   GCD
#
###############################################################################

doc"""
    gcd{T <: RingElem}(a::Nemo.FracElem{T}, b::Nemo.FracElem{T})
> Return a greatest common divisor of $a$ and $b$ if one exists. N.B: we define
> the GCD of $a/b$ and $c/d$ to be gcd$(ad, bc)/bd$, reduced to lowest terms.
> This requires the existence of a greatest common divisor function for the
> base ring.
"""
function gcd(a::Nemo.FracElem{T}, b::Nemo.FracElem{T}) where {T <: RingElem}
   check_parent(a, b)
   n = gcd(numerator(a)*denominator(b), denominator(a)*numerator(b))
   d = denominator(a)*denominator(b)
   g = gcd(n, d)
   return parent(a)(divexact(n, g), divexact(d, g))
end

################################################################################
#
#   Remove and valuation
#
################################################################################

doc"""
    remove{T <: RingElem}(z::Nemo.FracElem{T}, p::T)
> Return the tuple $n, x$ such that $z = p^nx$ where $x$ has valuation $0$ at
> $p$.
"""
function remove(z::Nemo.FracElem{T}, p::T) where {T <: RingElem}
   iszero(z) && error("Not yet implemented")
   v, d = remove(denominator(z), p)
   w, n = remove(numerator(z), p)
   return w-v, n//d
end

doc"""
    valuation{T <: RingElem}(z::Nemo.FracElem{T}, p::T)
> Return the valuation of $z$ at $p$.
"""
function valuation(z::Nemo.FracElem{T}, p::T) where {T <: RingElem}
   v, _ = remove(z, p)
   return v
end

###############################################################################
#
#   Unsafe operators and functions
#
###############################################################################

function zero!(c::Nemo.FracElem)
   c.num = zero!(c.num)
   if !isone(c.den)
      c.den = one(parent(c))
   end
   return c
end

function mul!(c::Nemo.FracElem{T}, a::Nemo.FracElem{T}, b::Nemo.FracElem{T}) where {T <: RingElem}
   g1 = gcd(numerator(a), denominator(b))
   g2 = gcd(numerator(b), denominator(a))
   c.num = divexact(numerator(a), g1)*divexact(numerator(b), g2)
   c.den = divexact(denominator(a), g2)*divexact(denominator(b), g1)
   return c
end

function addeq!(c::Nemo.FracElem{T}, a::Nemo.FracElem{T}) where {T <: RingElem}
   n = c.num*denominator(a) + numerator(a)*c.den
   c.den = mul!(c.den, c.den, denominator(a))
   g = gcd(n, c.den)
   c.num = divexact(n, g)
   c.den = divexact(c.den, g)
   return c
end

function add!(c::Nemo.FracElem{T}, a::Nemo.FracElem{T}, b::Nemo.FracElem{T}) where {T <: RingElem}
   n = c.num*denominator(a) + numerator(a)*c.den
   d = c.den*denominator(a)
   g = gcd(n, d)
   c.num = divexact(n, g)
   c.den = divexact(d, g)
   return c
end

function addeq!(c::Nemo.FracElem{T}, a::Nemo.FracElem{T}, b::Nemo.FracElem{T}) where {T <: RingElem}
   n = numerator(b)*denominator(a) + numerator(a)*denominator(b)
   c.den = mul!(c.den, denominator(b), denominator(a))
   g = gcd(n, d)
   c.num = divexact(n, g)
   c.den = divexact(c.den, g)
   return c
end

###############################################################################
#
#   Random functions
#
###############################################################################

function rand(S::Nemo.FracField{T}, v...) where {T <: RingElem}
   R = base_ring(S)
   n = rand(R, v...)
   d = R()
   while d == 0
      d = rand(R, v...)
   end
   return S(n, d)
end

###############################################################################
#
#   Promotion rules
#
###############################################################################

promote_rule(::Type{Frac{T}}, ::Type{Frac{T}}) where T <: RingElement = Frac{T}

function promote_rule(::Type{Frac{T}}, ::Type{U}) where {T <: RingElem, U <: RingElem}
   promote_rule(T, U) == T ? Frac{T} : Union{}
end

###############################################################################
#
#   Parent object call overloading
#
###############################################################################

function (a::FracField{T})(b::RingElement) where {T <: RingElement}
   return a(base_ring(a)(b))
end

function (a::FracField{T})() where {T <: RingElement}
   z = Frac{T}(zero(base_ring(a)), one(base_ring(a)))
   z.parent = a
   return z
end

function (a::FracField{T})(b::T) where {T <: RingElement}
   parent(b) != base_ring(a) && error("Could not coerce to fraction")
   z = Frac{T}(b, one(base_ring(a)))
   z.parent = a
   return z
end

function (a::FracField{T})(b::T, c::T) where {T <: RingElement}
   parent(b) != base_ring(a) && error("Could not coerce to fraction")
   parent(c) != base_ring(a) && error("Could not coerce to fraction")
   z = Frac{T}(b, c)
   z.parent = a
   return z
end

function (a::FracField{T})(b::T, c::Union{Integer, Rational, AbstractFloat}) where {T <: RingElement}
   parent(b) != base_ring(a) && error("Could not coerce to fraction")
   z = Frac{T}(b, base_ring(a)(c))
   z.parent = a
   return z
end

function (a::FracField{T})(b::Union{Integer, Rational, AbstractFloat}, c::T) where {T <: RingElement}
   parent(c) != base_ring(a) && error("Could not coerce to fraction")
   z = Frac{T}(base_ring(a)(b), c)
   z.parent = a
   return z
end

function (a::FracField{T})(b::Union{Integer, Rational, AbstractFloat}) where {T <: RingElement}
   z = Frac{T}(base_ring(a)(b), one(base_ring(a)))
   z.parent = a
   return z
end

function (a::FracField{T})(b::Integer, c::Integer) where {T <: RingElement}
   z = Frac{T}(base_ring(a)(b), base_ring(a)(c))
   z.parent = a
   return z
end

function (a::FracField{T})(b::Frac{T}) where {T <: RingElement}
   a != parent(b) && error("Could not coerce to fraction")
   return b
end

###############################################################################
#
#   FractionField constructor
#
###############################################################################

doc"""
    FractionField(R::Nemo.Ring; cached=true)
> Return the parent object of the fraction field over the given base ring $R$.
> If `cached == true` (the default), the returned parent object is cached so
> that it will always be returned by a call to the constructor when the same
> base ring $R$ is supplied.
"""
function FractionField(R::Nemo.Ring; cached=true)
   R2 = R
   T = elem_type(R)

   return FracField{T}(R, cached)
end


text:
# Metida
# Copyright © 2019-2020 Vladimir Arnautov aka PharmCat <mail@pharmcat.net>

using Metida

include("test.jl")


text:
# write your own tests here
using SimilaritySearch
using Test

include("testkvp.jl")
include("testknr.jl")


text:
#unum-macro.jl
#=
#creates a macro for inputting unums in a convenient fashion.
macro unum(prefix, parts...)
  environment_known::Bool = false
  #for now, set {4,6} as a default environment.
  ess::Integer = 4
  fss::Integer = 6
  uexpr::Expr
  #figure out the environment
  if prefix.head == :cell1d
    if (length(prefix.args) == 1)
      (prefix.args[1] != :auto) && throw(ParseError("invalid environment prefix"))
    else
      (length(prefix.args) > 2) && throw(ParseError("invalid environment prefix"))
      (length(prefix.args) == 2) && ((!isinteger(prefix.args[1])) || (!isinteger(prefix.args[2]))) && throw(ParseError("invalid environment prefix"))
      ess = prefix.args[1]
      fss = prefix.args[2]
      (ess > 6) && throw(ParseError("esizesize > 6 not currently supported"))
      (fss > 16) && throw(ParseError("fsizesize > 16 not currently supported"))
      environment_known = true
    end
    isempty(parts) && throw(ParseError("unum value missing"))
    uexpr = parts[1]
  else
    uexpr = prefix
  end
  isempty(uexpr) && return
  Meta.show_sexpr(parts)
end
=#

#dummy functions which will be used to set unum properties
↓() = nothing
⋯() = nothing
exact() = nothing
ulp() = nothing
export exact, ulp, ↓, ⋯

function __f2ubit(c::Function)
  ((c == ↓) || (c == exact)) && return z16
  ((c == ⋯) || (c == ulp)) && return UNUM_UBIT_MASK
  throw(ParseError("invalid ubit identifier supplied in unum definition"))
end

#alternative fallback, not using a macro
type __unum_shim
  ESS::Integer
  FSS::Integer
  val::Float64
  ubit::UInt16
end

colon(a::Array{Any,1}, b::Float64, c::Function) = __unum_shim(a[1], a[2], b, __f2ubit(c))
function colon(a::__unum_shim, b::Integer, c::Integer)
  x = convert(Unum{a.ESS, a.FSS}, a.val)
  unum(Unum{a.ESS, a.FSS}, UInt16(c), UInt16(b), (x.flags & UNUM_SIGN_MASK) | a.ubit, x.fraction, x.exponent)
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using RvSpectML
using DataFrames, CSV

function ccf_total(order_list_timeseries::AbstractChunkListTimeseries, line_list_df::DataFrame, pipeline::PipelinePlan; recalc::Bool = false,
                  output_fn_suffix::String = "", range_no_mask_change::Real=RvSpectMLBase.max_bc, ccf_mid_velocity::Real=0.0, v_step::Real=250.0,
                  mask_scale_factor::Real=1, mask_type::Symbol = :tophat, use_old::Bool = false, use_pixel_vars::Bool = false )
    if need_to(pipeline,:ccf_total) || recalc
      if verbose println("# Computing CCF.")  end
      @assert !need_to(pipeline,:extract_orders)
      @assert !need_to(pipeline,:clean_line_list_tellurics)
      if mask_type == :tophat
          mask_shape = CCF.TopHatCCFMask(order_list_timeseries.inst, scale_factor=mask_scale_factor)
      elseif mask_type == :gaussian
          mask_shape = CCF.GaussianCCFMask(order_list_timeseries.inst, σ_scale_factor=mask_scale_factor)
      elseif mask_type == :supergaussian
            mask_shape = CCF.SuperGaussianCCFMask(order_list_timeseries.inst, σ_scale_factor=mask_scale_factor)
      elseif mask_type == :halfcos
          mask_shape = CCF.CosCCFMask(order_list_timeseries.inst, scale_factor=mask_scale_factor)
      else
        @error("Requested mask shape (" * string(mask_type) * " not avaliable.")
      end

      line_list = CCF.BasicLineList(line_list_df.lambda, line_list_df.weight)
      ccf_plan = CCF.BasicCCFPlan(mask_shape = mask_shape, line_list=line_list, midpoint=ccf_mid_velocity, range_no_mask_change=range_no_mask_change, step=v_step)
      v_grid = CCF.calc_ccf_v_grid(ccf_plan)
      if use_old
        @time ccfs = CCF.calc_ccf_chunklist_timeseries_old(order_list_timeseries, ccf_plan)
      else
        @time ccfs = CCF.calc_ccf_chunklist_timeseries(order_list_timeseries, ccf_plan, use_pixel_vars=use_pixel_vars)
      end
      if save_data(pipeline, :ccf_total)
         CSV.write(joinpath(output_dir,target_subdir * "_ccfs" * output_fn_suffix * ".csv"),Tables.table(ccfs',header=Symbol.(v_grid)))
         #CSV.write(joinpath(output_dir,target_subdir * "_ccfs_expr.csv"),Tables.table(ccfs_expr',header=Symbol.(v_grid)))
      end
      set_cache!(pipeline, :ccf_total, (ccfs=ccfs, v_grid=v_grid) )
      dont_need_to!(pipeline,:ccf_total)
    end

    if has_cache(pipeline,:ccf_total) return read_cache(pipeline,:ccf_total)
    else   @error("Invalid pipeline state.")          end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1Handler <: SwaggerModel
    exec::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ExecAction } # spec name: exec
    httpGet::Any # spec type: Union{ Nothing, IoK8sApiCoreV1HTTPGetAction } # spec name: httpGet
    tcpSocket::Any # spec type: Union{ Nothing, IoK8sApiCoreV1TCPSocketAction } # spec name: tcpSocket

    function IoK8sApiCoreV1Handler(;exec=nothing, httpGet=nothing, tcpSocket=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1Handler, Symbol("exec"), exec)
        setfield!(o, Symbol("exec"), exec)
        validate_property(IoK8sApiCoreV1Handler, Symbol("httpGet"), httpGet)
        setfield!(o, Symbol("httpGet"), httpGet)
        validate_property(IoK8sApiCoreV1Handler, Symbol("tcpSocket"), tcpSocket)
        setfield!(o, Symbol("tcpSocket"), tcpSocket)
        o
    end
end # type IoK8sApiCoreV1Handler

const _property_map_IoK8sApiCoreV1Handler = Dict{Symbol,Symbol}(Symbol("exec")=>Symbol("exec"), Symbol("httpGet")=>Symbol("httpGet"), Symbol("tcpSocket")=>Symbol("tcpSocket"))
const _property_types_IoK8sApiCoreV1Handler = Dict{Symbol,String}(Symbol("exec")=>"IoK8sApiCoreV1ExecAction", Symbol("httpGet")=>"IoK8sApiCoreV1HTTPGetAction", Symbol("tcpSocket")=>"IoK8sApiCoreV1TCPSocketAction")
Base.propertynames(::Type{ IoK8sApiCoreV1Handler }) = collect(keys(_property_map_IoK8sApiCoreV1Handler))
Swagger.property_type(::Type{ IoK8sApiCoreV1Handler }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1Handler[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1Handler }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1Handler[property_name]

function check_required(o::IoK8sApiCoreV1Handler)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1Handler }, name::Symbol, val)
end


text:
# Facto.jl
# Like a Herbrand Universe and Quereis on it

global flog = stdout

"""
is t a var?
"""
isvar(N, t)= t in N

"""
index of var in N
"""
function vix(N,t)
    findall(map(x->x==t,N))
end

"""
apply a substitution (N,V) to a atomic term t
"""
function subst(N,t,V)
    if !isvar(N,t); return t end
    ix = vix(N, t)
    if V[ix[1]] == N[ix[1]]; return t
    else; return V[ix[1]]
    end
end

"""
apply a substitution to q
"""
function substitution(N,q,V)
    map(t->subst(N,t,V),q)
end

"""
apply a substitution to Q
"""
function substitutionQ(N,Q,V)
    map(q->substitution(N,q,V),Q)
end

"""
apply a unify t1 with t2
"""
function unifytt(N,t1,t2)
 if t1 == t2; t2
 elseif isvar(N,t1); t2
 elseif isvar(N,t2); t1
 else false
 end
end

"""
unify of q and f
"""
function unifyqf(N,q,f)
 R=deepcopy(N)
 if length(q)!=length(f); return false end
 for ix in 1:length(q)
  t1=q[ix]
  t2=f[ix]
  nt=subst(N,t1,R)
  if nt == t2; continue end
  if !isvar(N,nt)
   if !isvar(N,t2); return false end
  else # nt!=t2
   os = R[vix(N,t1)][1]
   vs = N[vix(N,t1)][1]
   if (os != t2) && (os != vs)
    return false
   end
   for i in vix(N,nt)
     R[i] = t2
   end
#   R[vix(N,nt)[1]]=t2
  end
 end
 return R
end

"""
apply a unify m1 to m2
"""
function unifymm(N,m1,m2)
println(flog,"unifymm<m1,m2=$m1,$m2")
 if m2 == N; return m1 end
 nm = []
 for i in 1:length(m1)
  t1,t2 = m1[i],m2[i]
  utt=unifytt(N,t1,t2)
  if utt == false; return false end
  push!(nm,utt)
 end
 return nm
end

"""
unify m and M
"""
function unifymM(N,m,M)
 if N==m; return M end
 if isempty(M); return [m] end
 rM=[]
 for m2 in M
  rm = unifymm(N,m,m2)
  if rm == false
   continue
  elseif isempty(rm)
   push!(rM,m)
  elseif rm==N
   continue
  else
   push!(rM,rm)
  end
 end
 rM
end

"""
unify M and M'
"""
function unifyMM(N,M1,M2)
 if isempty(M1); return M2 end
 if isempty(M2); return M1 end
 rM =[]
 for m1 in M1
  M12=unifymM(N,m1,M2)
  append!(rM,M12)
 end
 return rM
end

#=-
"""
unify members of M
"""
function unifyM(N, M)
 if isempty(M); return [N] end
 newM = [M[1]]
 if length(M)==1; return newM end
 for m in M[2:end]
  newM = unifymM(N,m,newM)
 end
 newM
end
==#

"""
find solution in K for q
"""
function solveq(N,q,K)
 M=[]
 for k in K
   am = unifyqf(N,q,k)
   if am==false; continue
   else push!(M,am) end
 end
 uM=unique(M)
 if isempty(uM); return false end
 return uM
end

"""
solve q with M
"""
function solveqM(N,q,M,K)
 rM = []
 if isempty(M)
  rM = solveq(N,q,K)
  return rM
 end
 for m in M
  bm=deepcopy(m)
  qq = substitution(N,q,bm)

  nM = solveq(N,qq,K)
  if nM==false; continue end

  rrM = unifymM(N,bm,nM)  ## replace bm with rrM
  if rrM==false; continue end

   append!(rM, rrM)
 end
 urM=unique(rM)
 if isempty(urM); return false end
 return urM
end

"""
make all models of qxf
"""
function makemodels(N,Q,K)
 MQ = Array{Any,2}(undef,length(K),length(Q))
 for ki in 1:length(K)
 for qi in 1:length(Q)
  MQ[ki,qi] = unifyqf(N,Q[qi],K[ki])
 end
 end
 return MQ
end

"""
reduction from MQ to models
"""
function reduction(N,mp,MFQ,qi)
println(flog,"reduction<mp,qi=$mp,$qi")
 if qi > size(MFQ,2); return [mp] end
 rM=[]
 for fi in 1:length(MFQ[:,qi])
println(flog,">fi=$fi")
  if MFQ[fi,qi] == false
   continue
  end
  mi = unifymm(N,mp, MFQ[fi,qi])
println(flog,"unifymm>mi=$mi")
  if mi == false
   continue
  end

  Mi = reduction(N,mi,MFQ,qi+1)
println(flog,"reduction>Mi=$Mi")
  if Mi==false
    return false
  else
    append!(rM, Mi)
  end
 end
 return rM
end

"""
new solveQ solve Q
"""
function solveQ(N,Q,K)
println(flog,"solveQ<Q=$Q")
 MFQ=makemodels(N,Q,K)
#println(flog,"makemodels>MFQ=$MFQ")
 rr=reduction(N,N,MFQ,1)
println(flog,"solveQ>rr=$rr")
 return rr
end

### END of main funcions
"""
statements of query true
"""
function findFact(N,q,K)
 M=solveq(N,q,K)
 if isempty(M) || M == false
   print(statement(N,[q],N))
   println(" IS FALSE")
 else
   return(statements(N,[q],M))
 end
 return(:negative)
end

### UI
"""
convert q to statement
"""
statement(N,q,m) = substitution(N,q,m)

"""
make statements from Q and M
"""
function statements(N,Q,M)
 rM = []
 for am in M
 for aq in Q
  push!(rM, statement(N,aq,am))
 end
 end
 unique(rM)
end
"""
string statement
"""
function stringT(F)
 r = ""
 for T in F
  for t in T
   r = r * " " * string(t)
  end
 r = r * ".\n"
 end
 r
end

"""
print F in string mode
"""
function printF(F)
 if length(F) == 0
  println("[]")
 else
  print(stringT(F))
 end
end

"""
print queries  on models
"""
function printQM(N,Q,M)
 ss =statements(N,Q,M)
 printF(ss)
 println()
end

"""
print model
"""
function printm(m)
 if length(m) == 0
  println("[]")
  return
 end

 for t in m[1:end-1]
  print("$t,")
 end
 println("$(m[end])")
end

"""
print models
"""
function printM(M)
 for m in M
  printm(m)
 end
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
__precompile__()

module PortAudio

using SampledSignals
using Devectorize
using RingBuffers
using Compat
import Compat: UTF8String, view

# Get binary dependencies loaded from BinDeps
include( "../deps/deps.jl")
include("libportaudio.jl")

export PortAudioStream

# These sizes are all in frames
# the block size is what we request from portaudio if no blocksize is given.
# The ringbuffer and pre-fill will be twice the blocksize
const DEFAULT_BLOCKSIZE=4096
# data is passed to and from the ringbuffer in chunks with this many frames
# it should be at most the ringbuffer size, and must evenly divide into the
# the underlying portaudio buffer size. E.g. if PortAudio is running with a
# 2048-frame buffer period, the chunk size can be 2048, 1024, 512, 256, etc.
const CHUNKSIZE=128

function __init__()
    # initialize PortAudio on module load
    swallow_stderr() do
        Pa_Initialize()
    end

    # the portaudio callbacks are parametric on the sample type
    global const pa_callbacks = Dict{Type, Ptr{Void}}()

    for T in (Float32, Int32, Int16, Int8, UInt8)
        pa_callbacks[T] = cfunction(portaudio_callback, Cint,
            (Ptr{T}, Ptr{T}, Culong, Ptr{Void}, Culong,
            Ptr{CallbackInfo{T}}))
    end
end

function versioninfo(io::IO=STDOUT)
    println(io, Pa_GetVersionText())
    println(io, "Version Number: ", Pa_GetVersion())
end

type PortAudioDevice
    name::UTF8String
    hostapi::UTF8String
    maxinchans::Int
    maxoutchans::Int
    defaultsamplerate::Float64
    idx::PaDeviceIndex
end

PortAudioDevice(info::PaDeviceInfo, idx) = PortAudioDevice(
        unsafe_string(info.name),
        unsafe_string(Pa_GetHostApiInfo(info.host_api).name),
        info.max_input_channels,
        info.max_output_channels,
        info.default_sample_rate,
        idx)

function devices()
    ndevices = Pa_GetDeviceCount()
    infos = PaDeviceInfo[Pa_GetDeviceInfo(i) for i in 0:(ndevices - 1)]
    PortAudioDevice[PortAudioDevice(info, idx-1) for (idx, info) in enumerate(infos)]
end

# not for external use, used in error message printing
devnames() = join(["\"$(dev.name)\"" for dev in devices()], "\n")

"""Give a pointer to the given field within a Julia object"""
function fieldptr{T}(obj::T, field::Symbol)
    fieldnum = findfirst(fieldnames(T), field)
    offset = fieldoffset(T, fieldnum)
    FT = fieldtype(T, field)

    Ptr{FT}(pointer_from_objref(obj) + offset)
end

# we want this to be immutable so we can stack allocate it
immutable CallbackInfo{T}
    inchannels::Int
    inbuf::LockFreeRingBuffer{T}
    outchannels::Int
    outbuf::LockFreeRingBuffer{T}
    synced::Bool
end

# paramaterized on the sample type and sampling rate type
type PortAudioStream{T, U}
    samplerate::U
    blocksize::Int
    stream::PaStream
    sink # untyped because of circular type definition
    source # untyped because of circular type definition
    bufinfo::CallbackInfo{T} # immutable data used in the portaudio callback

    # this inner constructor is generally called via the top-level outer
    # constructor below
    function PortAudioStream(indev::PortAudioDevice, outdev::PortAudioDevice,
            inchans, outchans, sr, blocksize, synced)
        inchans = inchans == -1 ? indev.maxinchans : inchans
        outchans = outchans == -1 ? outdev.maxoutchans : outchans
        inparams = (inchans == 0) ?
            Ptr{Pa_StreamParameters}(0) :
            Ref(Pa_StreamParameters(indev.idx, inchans, type_to_fmt[T], 0.0, C_NULL))
        outparams = (outchans == 0) ?
            Ptr{Pa_StreamParameters}(0) :
            Ref(Pa_StreamParameters(outdev.idx, outchans, type_to_fmt[T], 0.0, C_NULL))
        this = new(sr, blocksize, C_NULL)
        finalizer(this, close)
        this.sink = PortAudioSink{T, U}(outdev.name, this, outchans, blocksize*2)
        this.source = PortAudioSource{T, U}(indev.name, this, inchans, blocksize*2)
        if synced && inchans > 0 && outchans > 0
            # we've got a synchronized duplex stream. initialize with the output buffer full
            write(this.sink, SampleBuf(zeros(T, blocksize*2, outchans), sr))
        end
        this.bufinfo = CallbackInfo(inchans, this.source.ringbuf,
                                    outchans, this.sink.ringbuf, synced)
        this.stream = swallow_stderr() do
            Pa_OpenStream(inparams, outparams, float(sr), blocksize,
                paNoFlag, pa_callbacks[T], fieldptr(this, :bufinfo))
        end

        Pa_StartStream(this.stream)

        this
    end
end

# this is the top-level outer constructor that all the other outer constructors
# end up calling
function PortAudioStream(indev::PortAudioDevice, outdev::PortAudioDevice,
        inchans=2, outchans=2; eltype=Float32, samplerate=-1, blocksize=DEFAULT_BLOCKSIZE, synced=false)
    if samplerate == -1
        sampleratein = rationalize(indev.defaultsamplerate) * Hz;
        samplerateout = rationalize(outdev.defaultsamplerate) * Hz;
        if inchans > 0 && outchans > 0 && sampleratein != samplerateout
            error("""
            Can't open duplex stream with mismatched samplerates (in: $sampleratein, out: $samplerateout).
                   Try changing your sample rate in your driver settings or open separate input and output
                   streams""")
        elseif inchans > 0
            samplerate = sampleratein
        else
            samplerate = samplerateout
        end
    end
    PortAudioStream{eltype, typeof(samplerate)}(indev, outdev, inchans, outchans, samplerate, blocksize, synced)
end

# handle device names given as streams
function PortAudioStream(indevname::AbstractString, outdevname::AbstractString, args...; kwargs...)
    indev = nothing
    outdev = nothing
    for d in devices()
        if d.name == indevname
            indev = d
        end
        if d.name == outdevname
            outdev = d
        end
    end
    if indev == nothing
        error("No device matching \"$indevname\" found.\nAvailable Devices:\n$(devnames())")
    end
    if outdev == nothing
        error("No device matching \"$outdevname\" found.\nAvailable Devices:\n$(devnames())")
    end

    PortAudioStream(indev, outdev, args...; kwargs...)
end

# if one device is given, use it for input and output, but set inchans=0 so we
# end up with an output-only stream
function PortAudioStream(device::PortAudioDevice, inchans=2, outchans=2; kwargs...)
    PortAudioStream(device, device, inchans, outchans; kwargs...)
end
function PortAudioStream(device::AbstractString, inchans=2, outchans=2; kwargs...)
    PortAudioStream(device, device, inchans, outchans; kwargs...)
end

# use the default input and output devices
function PortAudioStream(inchans=2, outchans=2; kwargs...)
    inidx = Pa_GetDefaultInputDevice()
    indevice = PortAudioDevice(Pa_GetDeviceInfo(inidx), inidx)
    outidx = Pa_GetDefaultOutputDevice()
    outdevice = PortAudioDevice(Pa_GetDeviceInfo(outidx), outidx)
    PortAudioStream(indevice, outdevice, inchans, outchans; kwargs...)
end

function Base.close(stream::PortAudioStream)
    if stream.stream != C_NULL
        Pa_StopStream(stream.stream)
        Pa_CloseStream(stream.stream)
        close(stream.source)
        close(stream.sink)
        stream.stream = C_NULL
    end

    nothing
end

Base.isopen(stream::PortAudioStream) = stream.stream != C_NULL

SampledSignals.samplerate(stream::PortAudioStream) = stream.samplerate
Base.eltype{T, U}(stream::PortAudioStream{T, U}) = T

Base.read(stream::PortAudioStream, args...) = read(stream.source, args...)
Base.read!(stream::PortAudioStream, args...) = read!(stream.source, args...)
Base.write(stream::PortAudioStream, args...) = write(stream.sink, args...)
Base.write(sink::PortAudioStream, source::PortAudioStream, args...) = write(sink.sink, source.source, args...)
Base.flush(stream::PortAudioStream) = flush(stream.sink)

function Base.show(io::IO, stream::PortAudioStream)
    println(io, typeof(stream))
    println(io, "  Samplerate: ", samplerate(stream))
    print(io, "  Buffer Size: ", stream.blocksize, " frames")
    if nchannels(stream.sink) > 0
        print(io, "\n  ", nchannels(stream.sink), " channel sink: \"", stream.sink.name, "\"")
    end
    if nchannels(stream.source) > 0
        print(io, "\n  ", nchannels(stream.source), " channel source: \"", stream.source.name, "\"")
    end
end

# Define our source and sink types
for (TypeName, Super) in ((:PortAudioSink, :SampleSink),
                          (:PortAudioSource, :SampleSource))
    @eval type $TypeName{T, U} <: $Super
        name::UTF8String
        stream::PortAudioStream{T, U}
        chunkbuf::Array{T, 2}
        ringbuf::LockFreeRingBuffer{T}
        nchannels::Int

        function $TypeName(name, stream, channels, ringbufsize)
            # portaudio data comes in interleaved, so we'll end up transposing
            # it back and forth to julia column-major
            chunkbuf = zeros(T, channels, CHUNKSIZE)
            ringbuf = LockFreeRingBuffer(T, ringbufsize * channels)
            new(name, stream, chunkbuf, ringbuf, channels)
        end
    end
end

SampledSignals.nchannels(s::Union{PortAudioSink, PortAudioSource}) = s.nchannels
SampledSignals.samplerate(s::Union{PortAudioSink, PortAudioSource}) = samplerate(s.stream)
SampledSignals.blocksize(s::Union{PortAudioSink, PortAudioSource}) = s.stream.blocksize
Base.eltype{T, U}(::Union{PortAudioSink{T, U}, PortAudioSource{T, U}}) = T
Base.close(s::Union{PortAudioSink, PortAudioSource}) = close(s.ringbuf)

function Base.show{T <: Union{PortAudioSink, PortAudioSource}}(io::IO, stream::T)
    println(io, T, "(\"", stream.name, "\")")
    print(io, nchannels(stream), " channels")
end

function Base.flush(sink::PortAudioSink)
    while nwritable(sink.ringbuf) < length(sink.ringbuf)
        wait(sink.ringbuf)
    end
end

function SampledSignals.unsafe_write(sink::PortAudioSink, buf::Array, frameoffset, framecount)
    nwritten = 0
    while nwritten < framecount
        while nwritable(sink.ringbuf) == 0
            wait(sink.ringbuf)
        end
        # in 0.4 transpose! throws an error if the range is a UInt
        writable = div(nwritable(sink.ringbuf), nchannels(sink))
        towrite = Int(min(writable, CHUNKSIZE, framecount-nwritten))
        # make a buffer of interleaved samples
        transpose!(view(sink.chunkbuf, :, 1:towrite),
                   view(buf, (1:towrite)+nwritten+frameoffset, :))
        write(sink.ringbuf, sink.chunkbuf, towrite*nchannels(sink))

        nwritten += towrite
    end

    nwritten
end

function SampledSignals.unsafe_read!(source::PortAudioSource, buf::Array, frameoffset, framecount)
    nread = 0
    while nread < framecount
        while nreadable(source.ringbuf) == 0
            wait(source.ringbuf)
        end
        # in 0.4 transpose! throws an error if the range is a UInt
        readable = div(nreadable(source.ringbuf), nchannels(source))
        toread = Int(min(readable, CHUNKSIZE, framecount-nread))
        read!(source.ringbuf, source.chunkbuf, toread*nchannels(source))
        # de-interleave the samples
        transpose!(view(buf, (1:toread)+nread+frameoffset, :),
                   view(source.chunkbuf, :, 1:toread))

        nread += toread
    end

    nread
end

# This is the callback function that gets called directly in the PortAudio
# audio thread, so it's critical that it not interact with the Julia GC
function portaudio_callback{T}(inptr::Ptr{T}, outptr::Ptr{T},
        nframes, timeinfo, flags, userdata::Ptr{CallbackInfo{T}})
    info = unsafe_load(userdata)
    # if there are no channels, treat it as if we can write as many 0-frame channels as we want
    framesreadable = info.outchannels > 0 ? div(nreadable(info.outbuf), info.outchannels) : nframes
    frameswritable = info.inchannels > 0 ? div(nwritable(info.inbuf), info.inchannels) : nframes
    if info.synced
        framesreadable = min(framesreadable, frameswritable)
        frameswritable = framesreadable
    end
    towrite = min(frameswritable, nframes) * info.inchannels
    toread = min(framesreadable, nframes) * info.outchannels

    read!(info.outbuf, outptr, toread)
    write(info.inbuf, inptr, towrite)

    if framesreadable < nframes
        outsamples = nframes * info.outchannels
        # xrun, copy zeros to outbuffer
        # TODO: send a notification to an error msg ringbuf
        memset(outptr+sizeof(T)*toread, 0, sizeof(T)*(outsamples-toread))
    end

    paContinue
end

"""Call the given function and discard stdout and stderr"""
function swallow_stderr(f)
    origerr = STDERR
    (errread, errwrite) = redirect_stderr()
    result = f()
    redirect_stderr(origerr)
    close(errwrite)
    close(errread)

    result
end

memset(buf, val, count) = ccall(:memset, Ptr{Void},
    (Ptr{Void}, Cint, Csize_t),
    buf, val, count)

end # module PortAudio


text:
using Phylo
using Phylo: Rootedness, Rooted, TreeType, TraversalOrder
using Phylo: AbstractNode, AbstractBranch, AbstractTree
using SimpleTraits
using Unitful

@traitdef PreferNodeObjects{X}
@traitimpl PreferNodeObjects{X} <- _prefernodeobjects(X)

@traitdef PreferBranchObjects{X}
@traitimpl PreferBranchObjects{X} <- _preferbranchobjects(X)

@traitdef MatchNodeType{T, N}
@traitimpl MatchNodeType{T, N} <- _matchnodetype(T, N)

@traitdef MatchNodeTypes{T, N1, N2}
@traitimpl MatchNodeTypes{T, N1, N2} <- _matchnodetypes(T, N1, N2)

@traitdef HoldsNodeData{T, ND}
@traitimpl HoldsNodeData{T, ND} <- (_nodedatatype(T) ≡ ND)

@traitdef MatchBranchType{T, B}
@traitimpl MatchBranchType{T, B} <- _matchbranchtype(T, B)

@traitdef MatchBranchNodeType{T, B, N}
@traitimpl MatchBranchNodeType{T, B, N} <- _matchbranchnodetype(T, N)

@traitdef MatchTreeNameType{T, TN}
@traitimpl MatchTreeNameType{T, TN} <- _matchtreenametype(T, TN)

"""
    _prefernodeobjects(::Type{<:AbstractTree})
    _prefernodeobjects(::Type{<:AbstractNode})

Does this tree or node type prefer nodes to be objects or names? Must be
implemented for every node type.
"""
function _prefernodeobjects end
_prefernodeobjects(::Type{<:AbstractTree{TT, RT, NL, N, B}}) where
{TT, RT, NL, N, B} = _prefernodeobjects(N)

"""
    _preferbranchobjects(::Type{<:AbstractTree})
    _preferbranchobjects(::Type{<:AbstractBranch})

Does this tree or branch type prefer branches to be objects or names? Must be
implemented for every branch type.
"""
function _preferbranchobjects end
_preferbranchobjects(::Type{<:AbstractTree{TT, RT, NL, N, B}}) where
{TT, RT, NL, N, B} = _preferbranchobjects(B)

"""
    _matchnodetype(::Type{<:AbstractTree{TT, RT, NL, N, B}}, ::Type{N})
    _matchnodetype(::Type{<:AbstractTree{TT, RT, NL, N, B}}, ::Type{NL})

Does this tree type prefer the node or node label type provided?
"""
function _matchnodetype end
_matchnodetype(::Type{<:AbstractTree{TT, RT, NL, N, B}},
               ::Type{N}) where {TT, RT, NL, N, B} = _prefernodeobjects(N)
_matchnodetype(::Type{<:AbstractTree{TT, RT, NL, N, B}},
               ::Type{NL}) where {TT, RT, NL, N, B} = !_prefernodeobjects(N)
_matchnodetypes(::Type{<:AbstractTree}, ::Type{<:Any}, ::Type{<:Any}) = false
_matchnodetypes(::Type{T}, ::Type{N}, ::Type{N}) where {T, N} =
    _prefernodeobjects(T, N)

    """
    _matchbranchtype(::Type{<:AbstractTree}, ::Type{<:AbstractBranch})

Does this tree type prefer the branch or branch label type provided?
"""
function _matchbranchtype end
_matchbranchtype(::Type{<:AbstractTree{TT, RT, NL, N, B}},
                 ::Type{B}) where {TT, RT, NL, N, B} =
                     _preferbranchobjects(B)
_matchbranchtype(::Type{<:AbstractTree{TT, RT, NL, N, B}},
                 ::Type{Int}) where {TT, RT, NL, N, B} =
                     !_preferbranchobjects(B)

"""
    _matchbranchnodetype(::Type{<:AbstractTree},
                         ::Type{<:AbstractBranch},
                         ::Type{<:AbstractNode})

Does this tree type prefer the branch and node types provided?
"""
function _matchbranchnodetype end
_matchbranchtype(::Type{T}, ::Type{B}, ::Type{N}) where
    {TT, RT, NL, N, B, T <: AbstractTree{TT, RT, NL, N, B}} =
                    _matchbranchtype(T, B) && _matchnodetype(T, N)

"""
    _matchtreenametype(::Type{<:AbstractTree}, ::Type{X})

Does this tree type prefer the node or node label type provided?
"""
_matchtreenametype(::Type{T}, ::Type{TN}) where {T <: AbstractTree, TN} =
    (TN ≡ _treenametype(T))

# Need to be able to create new node and branch labels
function _newlabel(ids::Vector{Label}) where Label <: Integer
    return isempty(ids) ? 1 : maximum(ids) + 1
end
function _newlabel(ids::Vector{Vector{Label}}) where Label <: Integer
    return mapreduce(_newlabel, max, 1, ids)
end
_newlabel(names::Vector{String}, prefix::String) =
    prefix * "$(_newlabelfree(names, prefix))"
_newlabel(names::Vector{Vector{String}}, prefix::String) =
    prefix * "$(mapreduce(_newlabelfree, max, 1, names))"

function _newlabelfree(names::Vector{String}, prefix)
    num = length(names) + 1
    if "$prefix$num" ∉ names
        return num
    end
    len = length(prefix)
    goodnames = [name[(len - 1):end] for name in names
                 if length(name) > len && name[1:len] == prefix]
    num = length(goodnames) + 1
    while ("$num" ∈ goodnames)
        num += 1
    end
    return num
end

"""
    _newbranchlabel(tree::AbstractTree)

Returns a new unique branch name for a tree.
"""
function _newbranchlabel end
_newbranchlabel(tree::AbstractTree) = _newlabel(collect(getbranchnames(tree)))

"""
    _newnodelabel(tree::AbstractTree)

Returns a new unique node name for a tree.
"""
function _newnodelabel end
_newnodelabel(tree::AbstractTree{TT, RT, String, N, B}) where
{TT, RT, N, B} = _newlabel(collect(getnodenames(tree)), "Node ")
_newnodelabel(tree::AbstractTree{TT, RT, I, N, B}) where
    {TT, RT, I <: Integer, N, B} = _newlabel(collect(getnodenames(tree)))

# AbstractTree type methods
"""
    _treenametype(::Type{AbstractTree})

Returns the label type for a tree type. Must be implemented for any tree type.
"""
function _treenametype end
_treenametype(::Type{<: AbstractTree{OneTree}}) = Int

# AbstractTree methods
"""
    _ntrees(::AbstractTree)

Returns the number of trees in an object. Must be implemented for any ManyTrees
type.
"""
function _ntrees end
_ntrees(tree::AbstractTree{OneTree}) = 1
_ntrees(tree::AbstractTree{ManyTrees}) = length(_gettreenames(tree))

"""
    _gettrees(::AbstractTree)

Returns the trees in an object. Must be implemented for any ManyTrees type.
"""
function _gettrees end
_gettrees(tree::AbstractTree{OneTree}) = [tree]

"""
    _gettreenames(::AbstractTree)

Returns the names for the trees. Can be implemented for any ManyTrees type.
"""
function _gettreenames end
_gettreenames(tree::AbstractTree) = _gettreename.(_gettrees(tree))

"""
    _gettreename(::AbstractTree)

Returns the name for a single tree. Should be implemented for any OneTree
type where they have names.
"""
function _gettreename end
_gettreename(tree::AbstractTree{OneTree}) = 1

"""
    _gettree(::Pair{Label, AbstractTree})
    _gettree(::AbstractTree, id)

Returns a tree - either itself if it is a single tree, or the single tree
in a set with label id. Must be implemented for any ManyTrees type.
"""
function _gettree end
_gettree(treepair::Pair{Label, <: AbstractTree{OneTree}}) where Label =
    treepair[2]
@traitfn function _gettree(tree::T,
                           id::TN) where {T <: AbstractTree{OneTree}, TN;
                                          MatchTreeNameType{T, TN}}
    id == _gettreename(tree) || throw(BoundsError(tree, id))
    return tree
end

"""
    _getnodes(tree::AbstractTree{OneTree})

Returns a vector of nodes for a OneTree tree. Either _getnodes() must be
implemented for any OneTree tree type.
"""
function _getnodes end
_getnodes(tree::AbstractTree{OneTree}, order::TraversalOrder) =
    _traversal(tree, order)

"""
    _getnodenames(tree::AbstractTree{OneTree})

Returns a vector of node names for a OneTree tree. Can
be implemented for any OneTree tree type, especially PreferNodeObjects trees.
"""
function _getnodenames end
@traitfn _getnodenames(tree::T, order::TraversalOrder) where
{T <: AbstractTree{OneTree}; !PreferNodeObjects{T}} = _getnodes(tree)
@traitfn _getnodenames(tree::T, order::TraversalOrder) where
{T <: AbstractTree{OneTree}; PreferNodeObjects{T}} =
    _getnodename.(tree, _getnodes(tree))

"""
    _nnodes(::AbstractTree)

Returns the number of nodes (internal nodes and leaves) in a single tree. May
be implemented for any OneTree tree type (otherwise infers from _getnodes()).
"""
_nnodes(tree::AbstractTree{OneTree}) = length(_getnodes(tree, anyorder))

"""
    _getleafnames(::AbstractTree)

Returns the leaf names of a tree. May be implemented for any
tree type (otherwise determined from _getnodenames() and _isleaf() functions).
"""
function _getleafnames end
_getleafnames(tree::AbstractTree{OneTree}, order::TraversalOrder) =
    [_getnodename(tree, node)
     for node in _traversal(tree, order) if _isleaf(tree, node)]
_getleafnames(tree::AbstractTree{ManyTrees}, order::TraversalOrder) =
    _getleafnames(_gettree(tree, first(_gettreenames(tree))), order)

"""
    _getleaves(::AbstractTree)

Returns the leaves (tips) of a single tree. May be implemented for any
OneTree type (otherwise determined from _getnodes() and _isleaf() functions).
"""
function _getleaves end
_getleaves(tree::AbstractTree{OneTree}, order::TraversalOrder) =
    [node for node in _traversal(tree, order) if _isleaf(tree, node)]

"""
    _nleaves(::AbstractTree)

Returns the number of leaves (tips) in a tree. May be implemented for any
tree type (otherwise determined from the _getleafnames() function).
"""
_nleaves(tree::AbstractTree) = length(_getleafnames(tree, anyorder))

"""
    _getroots(::AbstractTree)

Returns the root(s) of a tree. May be implemented for any
OneTree type (otherwise determined from _getnodes() and _isroot() functions).
"""
_getroots(tree::AbstractTree{OneTree, <: Rooted}) =
    [node for node in _getnodes(tree) if _isroot(tree, node)]

"""
    _getroot(::AbstractTree)

Returns the unique root of a rooted tree. May be implemented for any
OneTree type (otherwise determined from _getroots()).
"""
function _getroot end
function _getroot(tree::AbstractTree{OneTree, <: Rooted})
    @assert _nroots(tree) == 1 "More than one root for tree " *
        "(found $(_nroots(tree)))"
    return first(_getroots(tree))
end

"""
    _nroots(::AbstractTree)

Returns the number of roots (subtrees) in a OneTree tree. May be implemented
for any ManyRoots type (otherwise infers from _getroots()).
"""
function _nroots end
_nroots(tree::AbstractTree{OneTree, <: Rooted}) = length(_getroots(tree))
_nroots(tree::AbstractTree{OneTree, Unrooted}) = 0

"""
    _getnode(::AbstractTree, id)

Returns the node or name associated with id (which could be a name,
a node or a pair) from a tree. Must be implemented for any PreferNodeObjects
tree and node label type.
"""
function _getnode end
@traitfn _getnode(tree::T,
                  node::N) where {RT, NL, N,
                                  T <: AbstractTree{OneTree, RT, NL, N};
                                  PreferNodeObjects{T}} = node
@traitfn _getnode(tree::T,
                  pair::Pair{NL, N}) where {RT, NL, N,
                                            T <: AbstractTree{OneTree,
                                                              RT, NL, N};
                                            PreferNodeObjects{T}} = pair[2]
@traitfn _getnode(tree::T,
                  nodename::NL) where {RT, NL,
                                       T <: AbstractTree{OneTree, RT, NL};
                                       !PreferNodeObjects{T}} = nodename
@traitfn _getnode(tree::T,
                  pair::Pair{NL, N}) where {RT, NL, N,
                                            T <: AbstractTree{OneTree,
                                                              RT, NL, N};
                                            !PreferNodeObjects{T}} = pair[1]

"""
    _getnodename(::AbstractTree, id)

Returns the name of a node associated with id (which could be a name, a node
or a pair) from a tree. Must be implemented for PreferNodeObjects tree types.
"""
function _getnodename end
_getnodename(::AbstractTree{OneTree, RT, NL, N},
             pair::Pair{NL, N}) where {RT, NL, N} = pair[1]
_getnodename(::AbstractTree{OneTree, RT, NL, N},
             nodename::NL) where {RT, NL, N} = nodename

"""
    _hasnode(tree::AbstractTree, node[name])

Does the tree contain this node? Must be implemented for any PreferNodeObjects
tree type with a node label.
"""
function _hasnode end
_hasnode(tree::AbstractTree{OneTree, RT, NL},
         node::NL) where {RT, NL} = node ∈ _getnodenames(tree)
_hasnode(tree::AbstractTree{OneTree, RT, NL, N},
         node::N) where {RT, NL, N} = node ∈ _getnodes(tree)

"""
    _createnode!(tree::AbstractTree, nodename)

Must be implemented for any AbstractTree subtype.
"""
function _createnode! end

"""
    _deletenode!(tree::AbstractTree, nodename)

Must be implemented for any AbstractTree subtype.
"""
function _deletenode! end

"""
    _getbranch(::AbstractTree, id)

Returns the branch or name associated with id (which could be a name,
a branch or a pair) from a tree. Must be implemented for any PreferBranchObjects
tree and branch label type.
"""
function _getbranch end
@traitfn _getbranch(tree::T,
                    branch::B) where {RT, NL, N, B,
                                      T <: AbstractTree{OneTree, RT, NL, N, B};
                                      PreferBranchObjects{T}} = branch
@traitfn _getbranch(tree::T,
                    pair::Pair{Int, B}) where {RT, NL, N, B,
                                               T <: AbstractTree{OneTree, RT,
                                                                 NL, N, B};
                                               PreferBranchObjects{T}} = pair[2]
@traitfn _getbranch(tree::T,
                    branchname::Int) where {T <: AbstractTree{OneTree};
                                            !PreferBranchObjects{T}} =
                                                branchname
@traitfn _getbranch(tree::T,
                    pair::Pair{Int, B}) where {RT, NL, N, B,
                                               T <: AbstractTree{OneTree,
                                                                 RT, NL, N, B};
                                               !PreferBranchObjects{T}} =
                                                   pair[1]

"""
    _getbranchname(::AbstractTree, id)

Returns the name of a branch associated with id (which could be a name, a branch
or a pair) from a tree. Must be implemented for PreferBranchObjects tree types.
"""
function _getbranchname end
_getbranchname(::AbstractTree{OneTree, RT, NL, N},
               pair::Pair{NL, N}) where {RT, NL, N} = pair[1]
_getbranchname(::AbstractTree{OneTree, RT},
               branchname::Int) where RT = branchname

"""
    _getbranches(tree::AbstractTree)

Returns a vector of branches for a OneTree tree. Either _getbranches() or
_getbranchnames() must be implemented for any OneTree tree type.
"""
function _getbranches end

"""
    _nbranches(::AbstractTree)

Returns the number of branches in a single tree. May be implemented for any
OneTree tree type (otherwise infers from _getbranches()).
"""
_nbranches(tree::AbstractTree{OneTree}) = length(_getbranches(tree))

"""
    _getbranchnames(tree::AbstractTree)

Returns a vector of branch names for a OneTree tree. Either _getbranches() or
_getbranchnames() must be implemented for any OneTree tree type.
"""
function _getbranchnames end
@traitfn _getbranchnames(tree::T) where {T <: AbstractTree{OneTree};
                                         !PreferBranchObjects{T}} =
                                             _getbranches(tree)
@traitfn _getbranchnames(tree::T) where {T <: AbstractTree{OneTree};
                                         PreferBranchObjects{T}} =
                                             [_getbranchname(tree, branch)
                                              for branch in _getbranches(tree)]

"""
    _hasbranch(tree::AbstractTree, node[name])

Does the tree contain this branch? Must be implemented for any
PreferBranchObjects tree type with a branch label.
"""
function _hasbranch end
_hasbranch(tree::AbstractTree{OneTree},
           branch::Int) = branch ∈ _getbranchnames(tree)
@traitfn _hasbranch(tree::T,
                    branch::B) where {RT, NL, N, B,
                                      T <: AbstractTree{OneTree, RT, NL, N, B};
                                      PreferBranchObjects{T}} =
                                          branch ∈ _getbranches(tree)

"""
    _createbranch!(tree::AbstractTree, source, destination,
                   length, data)

Create a new branch and add it to a tree. Must be implemented for any
AbstractTree subtype.
"""
function _createbranch! end

"""
    _deletebranch!(tree::AbstractTree, branch)

Delete a branch, reoving it from a tree. Must be implemented for any
AbstractTree subtype.
"""
function _deletebranch! end

"""
    _hasrootheight(::AbstractTree)

"""
function _hasrootheight(::AbstractTree)
    return false
end

"""
    _getrootheight(::AbstractTree)

"""
function _getrootheight end

"""
    _setrootheight!(::AbstractTree, value)

"""
function _setrootheight! end

"""
    _clearrootheight!(::AbstractTree)

"""
function _clearrootheight! end

"""
    _validate!(::AbstractTree)

"""
function _validate! end
_validate!(::AbstractTree) = true

"""
    _traversal(tree::AbstractTree, order::TraversalOrder, todo, sofar)

Return an iterable object containing nodes in given order - preorder, inorder,
postorder or breadthfirst
"""
function _traversal end
_traversal(tree::AbstractTree{OneTree, <: Rooted}, order::TraversalOrder) =
    (order == anyorder ? _getnodes(tree) :
     _traversal(tree, order, collect(_getroots(tree))))
function _traversal(tree::AbstractTree{OneTree, <: Rooted},
                    order::TraversalOrder, todo,
                    sofar = eltype(todo)[])
    while !isempty(todo)
        if order == Phylo.breadthfirst
            append!(sofar, todo)
            children = eltype(todo)[]
            for node in todo
                append!(children, getchildren(tree, node))
            end
            todo = children
        else
            now = popfirst!(todo)
            children = getchildren(tree, now)
            if isempty(children)
                push!(sofar, now)
            else
                if order == Phylo.preorder
                    push!(sofar, now)
                    sofar = _traversal(tree, order, children, sofar)
                elseif order == Phylo.postorder
                    sofar = _traversal(tree, order, children, sofar)
                    push!(sofar, now)
                elseif order == Phylo.inorder
                    child = popfirst!(children)
                    sofar = _traversal(tree, order, [child], sofar)
                    push!(sofar, now)
                    sofar = _traversal(tree, order, children, sofar)
                end
            end
        end
    end
    return sofar
end

# AbstractNode methods
"""
    _isleaf(tree::AbstractTree, node)

Is the node a leaf? Does not need to be implemented for any node type –
inferred from _outdegree or _degree - unless tree knows which nodes are
leaves and not nodes.
"""
function _isleaf end
_isleaf(tree::AbstractTree{OneTree, <: Rooted}, node) =
    _outdegree(tree, _getnode(tree, node)) == 0
function _isleaf(tree::AbstractTree{OneTree, Unrooted}, node)
    _degree(tree, node) == 0 && return true
    _degree(tree, node) > 1 && return false
    return missing
end

"""
    _isroot(tree::AbstractTree, node)

Is the node a root node of the tree?
"""
function _isroot end
_isroot(tree::AbstractTree{OneTree, <: Rooted}, node) =
    !_hasinbound(tree, node)
_isroot(tree::AbstractTree{OneTree, Unrooted}, node) = false

"""
    _isinternal(tree::AbstractTree, node)

Is the node internal to the tree?
"""
function _isinternal end
_isinternal(tree::AbstractTree{OneTree, <: Rooted}, node) =
    _outdegree(tree, node) > 0 && _hasinbound(tree, node)
_isinternal(tree::AbstractTree{OneTree, Unrooted}, node) =
    _degree(tree, node) < 2 ? false : missing

"""
    _isunattached(tree::AbstractTree, node)

Does the node currently form its own (sub)tree?
"""
_isunattached(tree::AbstractTree, node) = _degree(tree, node) == 0

"""
    _indegree(tree::AbstractTree, node)

In degree of node. Can be implemented for rooted nodes, otherwise inferred
from _hasinbound.
"""
function _indegree end
_indegree(tree::AbstractTree{OneTree, <: Rooted}, node) =
    Int(_hasinbound(tree, node))
_indegree(tree::AbstractTree{OneTree, Unrooted}, node) =
    _degree(tree, node) == 0 ? 0 : missing

"""
    _hasinboundspace(tree::AbstractTree, node::AbstractNode)

Is there space for a new inbound connection on a node?
"""
function _hasinboundspace end
_hasinboundspace(tree::AbstractTree{OneTree}, node) =
    !_hasinbound(tree, node)

"""
    _outdegree(tree::AbstractTree, node::AbstractNode)

Out degree of node.
"""
function _outdegree end
_outdegree(tree::AbstractTree{OneTree, <: Rooted}, node) =
    length(_getoutbounds(tree, node))
_outdegree(tree::AbstractTree{OneTree, Unrooted},
           node::AbstractNode{Unrooted}) =
               _degree(tree, node) == 0 ? 0 : missing

"""
    _hasoutboundspace(tree::AbstractTree, node::AbstractNode)

Is there space for a new outbound connection on a node? Must be implemented if
a node has a limit on the number of outbound connections (eg for a binary tree)
"""
function _hasoutboundspace end
_hasoutboundspace(::AbstractTree{OneTree, <: Rooted}, _) = true

"""
    _hasspace(tree::AbstractTree, node::AbstractNode)

Is there space for a new connection on a node? Must be implemented if
a node has a limit on the number of connections (eg for a binary tree)
"""
function _hasspace end
_hasspace(tree::AbstractTree{OneTree, <: Rooted}, node) =
    _hasinboundspace(tree, node) | _hasoutboundspace(tree, node)
_hasspace(::AbstractTree{OneTree, Unrooted}, _) = true

"""
    _degree(tree::AbstractTree, node::AbstractNode)

Degree of node. Must be implemented for Unrooted nodes, otherwise can be
inferred from indegree and outdegree.
"""
function _degree end
_degree(tree::AbstractTree{OneTree, <: Rooted}, node) =
    _indegree(tree, node) + _outdegree(tree, node)

"""
    _hasinbound(tree::AbstractTree, node::AbstractNode)

Must be implemented for any AbstractNode subtype.
"""
function _hasinbound end

"""
    _getinbound(tree::AbstractTree, node::AbstractNode)

Get the inbound connection. Must be implemented for any rooted AbstractNode subtype.
"""
function _getinbound end

"""
    _getparent(tree::AbstractTree, node)

Return the parent node for this node. Can be implemented for Rooted node types.
"""
function _getparent end
_getparent(tree::AbstractTree{OneTree}, node) =
    _getnode(tree, _src(tree, _getinbound(tree, node)))

"""
    _addinbound!(tree::AbstractTree, node::AbstractNode, inbound)

Adds a branch to the input of a rooted node. Must be implemented for any rooted
AbstractNode subtype unless this happens when a branch is created.
"""
function _addinbound! end

"""
    _removeinbound!(tree::AbstractTree, node::AbstractNode, inbound)

Removes a branch from the input of a rooted node. Must be implemented for any
rooted AbstractNode subtype unless this happens when a branch is deleted.
"""
function _removeinbound! end

"""
    _getoutbounds(tree::AbstractTree, node::AbstractNode)

Returns the outbound connections of a rooted node. Must be implemented for any
rooted AbstractNode subtype.
"""
function _getoutbounds end

"""
    _addoutbound!(tree::AbstractTree, node::AbstractNode, branch)

Add an outbound branch to a rooted node. Must be implemented for any Rooted
AbstractNode subtype unless this happens when a branch is created.
"""
function _addoutbound! end

"""
    _removeoutbound!(tree::AbstractTree, node::AbstractNode, branch)

Remove an outbound branch from a rooted node. Must be implemented for any
AbstractNode subtype unless this happens when a branch is deleted.
"""
function _removeoutbound! end

"""
    _getchildren(tree::AbstractTree, node)
    _getchildren(tree::AbstractTree, nodename)

Return the child node(s) for this node. May be implemented for any rooted
AbstractNode subtype.
"""
function _getchildren end
_getchildren(tree::AbstractTree{OneTree, <: Rooted}, node) =
    [_dst(tree, branch) for branch in _getoutbounds(tree, node)]

"""
    _getconnections(tree::AbstractTree, node::AbstractNode)

Returns all of the connections of a node. Must be implemented for any unrooted
AbstractNode subtype, can be inferred from _getinbound and _getoutbounds for
a rooted node.
"""
function _getconnections end
_getconnections(tree::AbstractTree{OneTree}, node) =
    _hasinbound(tree, node) ?
    append!([_getinbound(tree, node)], _getoutbounds(tree, node)) :
    _getoutbounds(tree, node)

"""
    _getsiblings(tree::AbstractTree, node::AbstractNode)

Returns all of the siblings of a node. May be implemented for any
AbstractNode subtype, can be inferred from _getparent and _getchildren for
a rooted node or _getconnections for an unrooted node.
"""
function _getsiblings end
_getsiblings(tree::AbstractTree{OneTree, <: Rooted}, node) =
    _hasinbound(tree, node) ?
    append!([_getparent(tree, node)], _getchildren(tree, node)) :
    _getchildren(tree, node)
_getsiblings(tree::AbstractTree{OneTree, Unrooted}, node) =
    [_conn(tree, branch, node) for branch in _getconnections(tree, node)]

"""
    _addconnection!(tree::AbstractTree, node::AbstractNode, branch)

Add a connection to an unrooted node. Must be implemented for any unrooted
AbstractNode subtype unless this happens when a branch is added.
"""
function _addconnection! end
_addconnection!(::AbstractTree{OneTree, <: Rooted}, _, _) =
    error("Direction of connection must be specified for a rooted node")

"""
    _removeconnection!(tree::AbstractTree, node::AbstractNode, branch)

Remove a connection from an unrooted node. Must be implemented for any Unrooted
AbstractNode subtype unless this happens when a branch is deleted.
"""
function _removeconnection! end
_removeconnection!(::AbstractTree{OneTree, <: Rooted}, _, _) =
    error("Direction of connection must be specified for a rooted node")

"""
    _hasheight(tree::AbstractTree, nodename)

"""
function _hasheight(::AbstractTree, _)
    return false
end

"""
    _getheight(tree::AbstractTree, nodename)

"""
function _getheight(::AbstractTree, _)
    return missing
end

"""
    _setheight!(::AbstractTree, nodename, value)

"""
function _setheight!(::AbstractTree, _, value)
    throw(NullException())
    return value
end

# AbstractBranch methods
"""
    _src(branch::AbstractBranch)

Return source node for a branch. Must be implemented for any rooted
AbstractBranch subtype.
"""
function _src end

"""
    _dst(branch::AbstractBranch)

Return destination node for a branch. Must be implemented for any rooted
AbstractBranch subtype.
"""
function _dst end

"""
    _conns(tree::AbstractTree, branch::AbstractBranch)

Return a vector of connections for a branch. Must be implemented for any
Unrooted AbstractBranch subtype, otherwise can combine _src and _dst.
"""
function _conns end
_conns(tree::AbstractTree{OneTree, <: Rooted}, branch) =
    [_src(tree, branch), _dst(tree, branch)]

"""
    _conn(branch::AbstractBranch, exclude::AbstractNode)

Return the connection for a branch that isn't the `exclude` node. May be
implemented for any Unrooted AbstractBranch subtype, otherwise will use
_conns.
"""
function _conn end
function _conn(tree::AbstractTree{OneTree, <: Rooted}, branch, exclude)
    other = [node for node in _conns(tree, branch) if node !== exclude]
    @assert length(other) == 1 _getnodename(tree, exclude) *
        " is not connected to branch $(_getbranchname(tree, branch))"
    return first(other)
end

"""
    _getlength

Return length of a branch. May be implemented for any AbstractBranch
subtype.
"""
function _getlength end
_getlength(::AbstractTree, _) = missing

"""
    _leafinfotype(::Type{<:AbstractTree})

Returns the type of the leaf info data.
"""
function _leafinfotype end
_leafinfotype(::Type{<:AbstractTree}) = Nothing

"""
    _nodedatatype(::Type{<:AbstractTree})

Returns the type of the node info data.
"""
function _nodedatatype end
_nodedatatype(::Type{<:AbstractTree}) = Nothing

"""
    _branchdims(::Type{<:AbstractTree})

Returns the dimensions of the branch lengths for the tree.
"""
function _branchdims end
_branchdims(::Type{<:AbstractTree}) = NoDims

"""
    _branchdatatype(::Type{<:AbstractTree})

Returns the type of the branch info data.
"""
function _branchdatatype end
_branchdatatype(::Type{<:AbstractTree}) = Nothing

function _getnodedata end
_getnodedata(tree::AbstractTree{OneTree}, node, label) =
    _getnodedata(tree, node)[label]
function _setnodedata! end
_setnodedata!(tree::AbstractTree{OneTree}, node, label, value) =
    (_getnodedata(tree, node)[label] = value)

function _getbranchdata end
_getbranchdata(tree::AbstractTree, branch, label) =
    _getbranchdata(tree, branch)[label]
function _setbranchdata! end
_setbranchdata!(tree::AbstractTree, branch, label, value) =
    (_getbranchdata(tree, branch)[label] = value)

function _getleafinfo end
function _setleafinfo! end

"""
    _gettreeinfo(tree::AbstractTree)
    _gettreeinfo(tree::AbstractTree, treename)

Returns the info data associated with the tree(s).
"""
function _gettreeinfo end
_gettreeinfo(tree::AbstractTree{OneTree}) =
    Dict(_gettreename(tree) => Dict{String, Any}())
function _gettreeinfo(tree::AbstractTree{OneTree}, treename)
    @assert _gettreename(tree) == treename "No tree called $treename"
    return Dict{String, Any}()
end

"""
    _resetleaves!(::AbstractTree)

Fixes leaf naming after creation or deletion of nodes or branches. Must be
implemented by tree types where this is handled separately.
"""
function _resetleaves! end
_resetleaves!(::AbstractTree) = true


text:
using LinearMaps
using IterativeSolvers
using Statistics
import Base: @kwdef

# Mixing rules: (ρin, ρout) => ρnext, where ρout is produced by diagonalizing the
# Hamiltonian at ρin These define the basic fix-point iteration, that are then combined with
# acceleration methods (eg anderson). For the mixing interface we use `δF = ρout - ρin` and
# `δρ = ρnext - ρin`, such that the mixing interface is `mix(mixing, basis, δF; kwargs...) -> δρ`
# with the user being assumed to add this to ρin to get ρnext.
# All these methods attempt to approximate the inverse Jacobian of the SCF step,
# ``J^-1 = (1 - χ0 (vc + K_{xc}))^-1``, where vc is the Coulomb and ``K_{xc}`` the
# exchange-correlation kernel. Note that "mixing" is sometimes used to refer to the combined
# process of formulating the fixed-point and solving it; we call "mixing" only the first part
# The notation in this file follows Herbst, Levitt arXiv:2009.01665

@doc raw"""
Simple mixing: ``J^{-1} ≈ α``
"""
@kwdef struct SimpleMixing
    α::Real = 0.8
end
@timing "SimpleMixing" function mix(mixing::SimpleMixing, ::PlaneWaveBasis, δFs...; kwargs...)
    T = eltype(δFs[1])
    map(δFs) do δF  # Apply to both δF{total} and δF{spin} in the same way
        isnothing(δF) && return nothing
        T(mixing.α) * δF
    end
end

@doc raw"""
Kerker mixing: ``J^{-1} ≈ \frac{α |G|^2}{k_{TF}^2 + |G|^2}``
where ``k_{TF}`` is the Thomas-Fermi wave vector. For spin-polarized calculations
by default the spin density is not preconditioned. Unless a non-default value
for ``ΔDOS`` is specified. This value should roughly be the expected difference in density
of states (per unit volume) between spin-up and spin-down.

Notes:
- Abinit calls ``1/k_{TF}`` the dielectric screening length (parameter *dielng*)
"""
@kwdef struct KerkerMixing
    # Default parameters suggested by Kresse, Furthmüller 1996 (α=0.8, kTF=1.5Å⁻¹)
    # DOI 10.1103/PhysRevB.54.11169
    α::Real    = 0.8
    kTF::Real  = 0.8  # == sqrt(4π (DOS_α + DOS_β) / Ω)
    ΔDOS::Real = 0.0  # == (DOS_α - DOS_β) / Ω
end

@timing "KerkerMixing" function mix(mixing::KerkerMixing, basis::PlaneWaveBasis,
                                    δF::RealFourierArray, δFspin=nothing; kwargs...)
    T    = eltype(δF)
    G²   = [sum(abs2, G) for G in G_vectors_cart(basis)]
    kTF  = T.(mixing.kTF)
    ΔDOS = T.(mixing.ΔDOS)

    # For Kerker the model dielectric written as a 2×2 matrix in spin components is
    #     1 - [-DOSα      0] * [1 1]
    #         [    0  -DOSβ]   [1 1] * (4π/G²)
    # which maps (δρα, δρβ)ᵀ to (δFα, δFβ)ᵀ and where DOSα and DOSβ is the density
    # of states per unit volume in the spin-up and spin-down channels. After basis
    # transformation to a mapping (δρtot, δρspin)ᵀ to (δFtot, δFspin)ᵀ this becomes
    #     [(G² + kTF²)    0]
    #     [ 4π * ΔDOS    G²] / G²
    # where we defined kTF² = 4π * (DOSα + DOSβ) and ΔDOS = DOSα - DOSβ.
    # Gaussian elimination on this matrix yields for the linear system ε δρ = δF
    #     δρtot  = G² δFtot / (G² + kTF²)
    #     δρspin = δFspin - 4π * ΔDOS / (G² + kTF²) δFtot

    δρtot    = @. δF.fourier * G² / (kTF^2 + G²)
    δρtot[1] = δF.fourier[1]  # Copy DC component, otherwise it never gets updated
    if basis.model.n_spin_components == 1
        from_fourier(basis, T(mixing.α) * δρtot), nothing
    else
        δρspin = @. δFspin.fourier - δF.fourier * (4π * ΔDOS) / (kTF^2 + G²)
        from_fourier(basis, T(mixing.α) * δρtot), from_fourier(basis, T(mixing.α) * δρspin)
    end
end

@doc raw"""
The same as [`KerkerMixing`](@ref), but the Thomas-Fermi wavevector is computed
from the current density of states at the Fermi level.
"""
@kwdef struct KerkerDosMixing
    α::Real = 0.8
end
@timing "KerkerDosMixing" function mix(mixing::KerkerDosMixing, basis::PlaneWaveBasis,
                                       δF, δFspin=nothing; εF, ψ, eigenvalues, kwargs...)
    if basis.model.temperature == 0
        return mix(SimpleMixing(α=mixing.α), basis, δF, δFspin)
    else
        n_spin = basis.model.n_spin_components
        Ω = basis.model.unit_cell_volume
        dos  = [compute_dos(εF, basis, eigenvalues, spins=[σ]) / Ω for σ in 1:n_spin]
        kTF  = sqrt(4π * sum(dos))
        ΔDOS = n_spin == 2 ? dos[1] - dos[2] : 0.0
        mix(KerkerMixing(α=mixing.α, kTF=kTF, ΔDOS=ΔDOS), basis, δF, δFspin)
    end
end

@doc raw"""
We use a simplification of the Resta model DOI 10.1103/physrevb.16.2717 and set
``χ_0(q) = \frac{C_0 G^2}{4π (1 - C_0 G^2 / k_{TF}^2)}``
where ``C_0 = 1 - ε_r`` with ``ε_r`` being the macroscopic relative permittivity.
We neglect ``K_\text{xc}``, such that
``J^{-1} ≈ α \frac{k_{TF}^2 - C_0 G^2}{ε_r k_{TF}^2 - C_0 G^2}``

By default it assumes a relative permittivity of 10 (similar to Silicon).
`εr == 1` is equal to `SimpleMixing` and `εr == Inf` to `KerkerMixing`.
The mixing is applied to ``ρ`` and ``ρ_\text{spin}`` in the same way.
"""
@kwdef struct DielectricMixing
    α::Real   = 0.8
    kTF::Real = 0.8
    εr::Real  = 10
end
@timing "DielectricMixing" function mix(mixing::DielectricMixing, basis::PlaneWaveBasis,
                                        δFs...; kwargs...)
    T = eltype(δFs[1])
    εr = T(mixing.εr)
    kTF = T(mixing.kTF)
    εr == 1               && return mix(SimpleMixing(α=mixing.α), basis, δFs...)
    εr > 1 / sqrt(eps(T)) && return mix(KerkerMixing(α=mixing.α, kTF=kTF), basis, δFs...)

    C0 = 1 - εr
    Gsq = [sum(abs2, G) for G in G_vectors_cart(basis)]
    map(δFs) do δF  # Apply to both δF{total} and δF{spin} in the same way
        isnothing(δF) && return nothing
        δρ    = @. T(mixing.α) * δF.fourier * (kTF^2 - C0 * Gsq) / (εr * kTF^2 - C0 * Gsq)
        δρ[1] = δF.fourier[1]  # Copy DC component, otherwise it never gets updated
        from_fourier(basis, δρ)
    end
end

@doc raw"""
The model for the susceptibility is
```math
\begin{aligned}
    χ_0(r, r') &= (-D_\text{loc}(r) δ(r, r') + D_\text{loc}(r) D_\text{loc}(r') / D) \\
    &+ \sqrt{L(x)} \text{IFFT} \frac{C_0 G^2}{4π (1 - C_0 G^2 / k_{TF}^2)} \text{FFT} \sqrt{L(x)}
\end{aligned}
```
where ``C_0 = 1 - ε_r``, ``D_\text{loc}`` is the local density of states,
``D`` is the density of states and
the same convention for parameters are used as in [`DielectricMixing`](@ref).
Additionally there is the real-space localisation function `L(r)`.
For details see Herbst, Levitt 2020 arXiv:2009.01665

Important `kwargs` passed on to [`χ0Mixing`](@ref)
- `α`: Damping parameter
- `RPA`: Is the random-phase approximation used for the kernel (i.e. only Hartree kernel is
  used and not XC kernel)
- `verbose`: Run the GMRES in verbose mode.
"""
function HybridMixing(;εr=1.0, kTF=0.8, localisation=identity, kwargs...)
    χ0terms = [DielectricModel(εr=εr, kTF=kTF, localisation=localisation), LdosModel()]
    χ0Mixing(; χ0terms=χ0terms, kwargs...)
end

@doc raw"""
Generic mixing function using a model for the susceptibility composed of the sum of the `χ0terms`.
For valid `χ0terms` See the subtypes of `χ0Model`. The dielectric model is solved in
real space using a GMRES. Either the full kernel (`RPA=false`) or only the Hartree kernel
(`RPA=true`) are employed. `verbose=true` lets the GMRES run in verbose mode
(useful for debugging).
"""
@kwdef struct χ0Mixing
    α::Real   = 0.8
    RPA::Bool = true       # Use RPA, i.e. only apply the Hartree and not the XC Kernel
    χ0terms   = χ0Model[Applyχ0Model()]  # The terms to use as the model for χ0
    verbose::Bool = false  # Run the GMRES verbosely
end

@timing "χ0Mixing" function mix(mixing::χ0Mixing, basis, δF_tot::RealFourierArray,
                                δF_spin=nothing; ρin, ρ_spin_in, kwargs...)
    T = eltype(δF_tot)
    n_spin = basis.model.n_spin_components
    @assert basis.model.spin_polarization in (:none, :spinless, :collinear)

    # Initialise χ0terms and remove nothings (terms that don't yield a contribution)
    χ0applies = [χ0(basis; ρin=ρin, ρ_spin_in=ρ_spin_in, kwargs...) for χ0 in mixing.χ0terms]
    χ0applies = [apply for apply in χ0applies if !isnothing(apply)]

    # If no applies left, do not bother running GMRES and directly do simple mixing
    isempty(χ0applies) && return mix(SimpleMixing(α=mixing.α), basis, δF_tot, δF_spin)

    # Solve J δρ = δF with J = (1 - χ0 vc) and χ_0 given as the sum of the χ0terms
    devec(x) = reshape(x, size(δF_tot)..., n_spin)
    function Jop(x)
        δF = devec(x)  # [:, :, :, 1] is spin-up (or total), [:, :, :, 2] is spin-down

        # Apply Kernel (just vc for RPA and (vc + K_{xc}) if not RPA)
        @views if n_spin == 1
            x_tot  = from_real(basis, δF[:, :, :, 1])
            x_spin = nothing
        else
            x_tot  = from_real(basis, δF[:, :, :, 1] .+ δF[:, :, :, 2])
            x_spin = from_real(basis, δF[:, :, :, 1] .- δF[:, :, :, 2])
        end
        δV = apply_kernel(basis, x_tot, x_spin; ρ=ρin, ρspin=ρ_spin_in, RPA=mixing.RPA)

        # set DC of δV to zero (δV[1] is spin-up or total, δV[2] is spin-down)
        δV_DC = mean(mean(δV[σ].real) for σ in 1:n_spin)
        δV[1].real .-= δV_DC
        n_spin == 2 && (δV[2].real .-= δV_DC)

        JδF = copy(δF)
        for apply_term! in χ0applies
            apply_term!(JδF, δV, -1)  # JδF .-= χ0 * δV
        end
        vec(JδF .-= mean(JδF))  # Zero DC component in total density response
    end

    if n_spin == 1
        δF_updown = δF_tot.real
    else
        δF_updown = cat((δF_tot.real .+ δF_spin.real) ./ 2,  # spin-up
                        (δF_tot.real .- δF_spin.real) ./ 2,  # spin-down
                        dims=4)
    end
    δF_updown .-= mean(δF_updown)  # Zero DC of δF_updown
    J = LinearMap(Jop, length(δF_updown))
    x = gmres(J, vec(δF_updown), verbose=mixing.verbose)
    # TODO Further improvement: Adapt tolerance of gmres to norm(ρ_out - ρ_in)

    δρ = T(mixing.α) .* devec(x)  # Apply damping
    δρ .+= (sum(δF_tot.real) - sum(δρ)) / length(δF_tot)  # Set DC from δF

    @views if n_spin == 1
        from_real(basis, δρ[:, :, :, 1]), nothing
    else
        (from_real(basis, δρ[:, :, :, 1] .+ δρ[:, :, :, 2]),
         from_real(basis, δρ[:, :, :, 1] .- δρ[:, :, :, 2]))
    end
end