text:
# Copyright 2020 Digital Domain 3.0
#
# Licensed under the Apache License, Version 2.0 (the "Apache License")
# with the following modification; you may not use this file except in
# compliance with the Apache License and the following modification to it:
# Section 6. Trademarks. is deleted and replaced with:
#
# 6. Trademarks. This License does not grant permission to use the trade
#    names, trademarks, service marks, or product names of the Licensor
#    and its affiliates, except as required to comply with Section 4(c) of
#    the License and to reproduce the content of the NOTICE file.
#
# You may obtain a copy of the Apache License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the Apache License with the above modification is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the Apache License for the specific
# language governing permissions and limitations under the Apache License.

using Test
using DiscreteDifferentialGeometry
using Multivectors
using HalfEdges
using LinearAlgebra
using StaticArrays
using Base.Iterators

using DiscreteDifferentialGeometry: cotan

zerofun(x,y...) = zero(x)

include("KForms_runtests.jl")

module G2
  using Multivectors
  # need last param to be true for differential forms
  @generate_basis("++",false,true,true)
end
using .G2

module G3
  using Multivectors
  @generate_basis("+++",false,true,true)
end

using .G3

function cube()
  P = [0 0 0
       1 0 0
       0 0 1
       1 0 1
       1 1 1
       0 1 1
       1 1 0
       0 1 0
      ]
  P = SVector{3}.(eachrow(Float64.(P)))
  (Topology([1,2,3, 2,4,3, 3,4,6, 4,5,6, 2,5,4, 2,7,5, 7,8,5, 5,8,6, 1,3,6, 1,6,8, 1,8,2, 2,8,7]), P)
end

# not actually a D8, two tetrahedron glued together
function D8()
  w = sqrt(2.0); peak = 0.8164965809277261; base = -0.4082482904638631;
  P = [ -w/2 base 0.0
       w/2 base 0.0
       0.0 peak 0.0
       0.0 0.0 peak-base
       0.0 0.0 base-peak ]

  P = SVector{3}.(eachrow(Float64.(P)))
  (Topology([1,2,4, 3,4,2, 1,4,3, 1,5,2, 1,3,5, 5,3,2]), P)
end

@testset begin "interpolation 3D"
  e₁, e₂, e₃  = alle( G3, 3)[1:3]

  # an equlateral triangle
  l = -0.5e₁ + 0.0e₂; r = 0.5e₁+0.0e₂; t = 0.0e₁ + e₂(sqrt(3.0)/2.0);

  # integrate whitney interpolating form over the edge opposite the bottom left corner
  # should end up with the same value we set on the edge ( 1.0 )
  ddg_ = DiscreteDifferentialGeometry
  dx = 0.005;  rtval = 42.0;
  edge_rt = sum(map(x->ddg_.interpolate((l,r,t), [rtval,0.0,0.0], r + (t-r)*x)⋅(t-r)*dx, 0.0:dx:(1.0-dx)))
  @test edge_rt ≈ rtval
  @test norm(ddg_.interpolate((l,r,t), [1.0,1.0,1.0], (r+l+t)/3.0)) ≈ 0.0
  # interpolate gradient of function with value of 1.0 at l and 0.0 at opposite edge
  gradl = ddg_.interpolate((l,r,t), [0.0,1.0,1.0], (r+l+t)/3.0)
  @test gradl⋅(1.0e₁+1.0e₂) < -0.1
  @test norm(gradl⋅⋆(t-r)) ≈ 0.0

end

function interpolate_ga(pa::KVector,pb,pc,va,vb,vc)
  tri𝐼 = (pb-pa)∧(pc-pa)
  tri𝐼 = tri𝐼/Multivectors.norm_sqr(tri𝐼)
  vb*((pa-pc)⋅tri𝐼) + va*((pc-pb)⋅tri𝐼) + vc*((pb-pa)⋅tri𝐼)
end

interpolate_ga(pa::SVector{Float64},pb,pc,ea,eb,ec) = interpolate_ga(KVector(pa), KVector(pb), KVector(pc),ea,eb,ec)

function interpolate_vec(pa,pb,pc,va,vb,vc)
  N = (pb-pa)×(pc-pa)
  twoA = norm(N)
  N = N/twoA
  1.0/twoA*(N×(va*(pc-pb)+vb*(pa-pc)+vc*(pb-pa)))
end

@testset begin "operators"
  e₁, e₂, e₃  = alle( G3, 3)[1:3]
  e¹, e², e³ = alld( G3, 3)[1:3]
  𝐼 = pseudoscalar(e₁)

  #== hex grid

  1.732    1---2---3
          / \ / \ / \
  0.866  4---5---6---7
        / \ / \ / \ / \
  0    8---9---A---B---C
        \ / \ / \ / \ /
         D---E---F---16
          \ / \ / \ /
  -1.732   17--18--19
  ==#

  A = 10; B = 11; C = 12; D = 13; E = 14; F = 15;
  h = 1.732

  topo_hex = Topology([1,4,5, 2,1,5, 2,5,6, 3,2,6, 3,6,7,
                       4,8,9, 5,4,9, 5,9,A, 6,5,A, 6,A,B, 7,6,B, 7,B,C,
                       9,8,D, 9,D,E, A,9,E, A,E,F, B,A,F, B,F,16, C,B,16,
                       E,D,17, E,17,18, F,E,18, F,18,19, 16,F,19])

  P_hex =
    SVector{3, Float64}.([[-1,h,0],[0,h,0],[1,h,0],
                          [-1.5,.5h,0],[-.5,.5h,0],[.5,.5h,0],[1.5,.5h,0],
                          [-2,0,0],[-1,0,0],[0,0,0],[1,0,0],[2,0,0],
                          [-1.5,-.5h,0],[-.5,-.5h,0],[.5,-.5h,0],[1.5,-.5h,0],
                          [-1,-h,0],[0,-h,0],[1,-h,0]])
  P = P_hex; topo = topo_hex
  V = nvertices(topo); E = nedges(topo); F = nfaces(topo);

  @test iszero(𝑑(topo, 𝑑(topo))) == true
  @test size(array(𝑑(topo))) == (E, V)
  @test size(array(𝑑(topo,𝑑(topo)))) == (F, V)
  @test ⋆(topo,P,⋆(topo,P)).D |> tr ≈ nvertices(topo)
  @test ⋆(topo,P,⋆(topo,P)).D |> det ≈ 1.0

  @test size(array(⋆(topo, P))) == (V, V)
  @test size(array(⋆(topo, P, ⋆(topo, P)))) == (V, V)

  @test size(array(⋆(topo, P, 𝑑(topo)))) == (E, V)
  @test size(array(⋆(topo, P, ⋆(topo, P, 𝑑(topo))))) == (E, V)

  𝑑mesh = op->𝑑(topo, op)
  smesh = op->⋆(topo, P, op)

  # navigate this structure :
  #==
   Ω₀---𝑑₀-→Ω₁----𝑑₁-→Ω₂
   ↑        ↑         ↑
 ⋆₀|⋆₀⁻   ⋆₁|⋆₁⁻    ⋆₂|⋆₂⁻
   ↓        ↓         ↓
   Ω⋆₂←-𝑑₁--Ω⋆₁←-𝑑₀---Ω⋆₀
  ==#
  # around the outside
  @test 𝑑(topo) |> 𝑑mesh |> smesh |> 𝑑mesh |> 𝑑mesh |> smesh |> size∘array == (V,V)
  # accross top dipping down twice
  @test 𝑑(topo) |> smesh |> smesh |> 𝑑mesh |> smesh |> smesh |> size∘array == (F,V)
  # left loop
  @test 𝑑(topo) |> smesh |> 𝑑mesh |> size∘array == (V,V)
  # right loop
  @test 𝑑(topo) |> 𝑑mesh |> smesh |> 𝑑mesh |> smesh |> size∘array == (E,V)

  # solve a simple heat diffusion step
  dt = (1.0/24.0)
  L = array(Δ(topo, P))*dt + I
  b = zeros(nvertices(topo)); b[1] = 100.0; b[19] = -75.0;
  x = L\b
  @test x[10] > 0.0 # middle. more heat than cold
  @test x[5] > 0.0  # close to hotspot
  @test x[15] < 0.0 # close to coldspot
  @test x[8] > 0.0 # closer to hotspot
  @test poisson(L, b, zeros(length(b)), []) ≈ x
  @test poisson(L, b, ones(length(b)), [5])[5] == 1.0
  @test poisson(L, b, ones(length(b)), [5])[6] != 1.0

  ## geodesic distance.
  ## 3 steps
  ## 1) integrate heat flow ( heat diffusion ).

  #    Going to use symmetric laplacian, could use L from before if we wanted to
  Lc = array(Δcell(topo, P))

  # can get vertex area matrix M from the hodge star taking values of 1.0 on a vertex to the vertex cell
  M = ⋆(topo, P, PrimalKForm{0}(ones(nvertices(topo)))) |> array

  # or the hodge star method that implies
  @test ⋆(topo,P) == ⋆(topo, P, PrimalKForm{0}(ones(nvertices(topo))))

  L = M + dt*Lc
  @test issymmetric(Lc) == true
  @test issymmetric(L) == true

  source = 5;
  b = zeros(nvertices(topo)); b[source] = 1.0
  # can leave this ⋆ out if desired, since entries of b are integrated dirac delta functions
  # leaving it in to show general way of using Δcell rather than Δ
  b = ⋆(topo,P)*b

  u = L\b

  # u holds the temperatures at each vertex

  # with vector algebra
  vec∇face = map( HalfEdges.polygons(topo) ) do abc
    a,b,c = abc
    pa,pb,pc = P[abc]
    N = (pb-pa)×(pc-pa)
    twoA = norm(N)
    N = N/twoA
    1.0/twoA*(N×(u[a]*(pc-pb)+u[b]*(pa-pc)+u[c]*(pb-pa)))
  end

  alledge = HalfEdges.edges(topo);
  allface = HalfEdges.polygons(topo);

  # a 1-form taking values from vertices to edges.
  dϕ = 𝑑(topo, PrimalKForm{0}(ones(nvertices(topo))))
  # apply to u (values on vertices) and get gradient of u on edges
  ∇u = dϕ*u

  # interpolate with whitney 1-forms to triangle barycentres
  ∇face = map(enumerate(IncidentEdges(topo, FaceHandle))) do (iface, edgeindices)
    abc = (p->KVector(p, 1𝐼)).(P[allface[iface]])
    # triangles vertices and edges are ordered so first edge points at first vertex
    # so need to loop almost all the way around
    nedge = length(edgeindices)
    oei = collect(take(drop(cycle(edgeindices), nedge-1), nedge))
    # order so value on edge is the one opposite each vertex
    edgeval = ∇u[oei]
    faceabc = allface[iface]
    orient = map(ie->HalfEdges.orientation(faceabc, alledge[ie]), oei)
    edgeval = orient .* edgeval
    centrep = sum(abc)/3.0
    interpolate( abc, edgeval, centrep )
  end

  # similar to closed-form vector algebra but using Geometric Algebra
  ga∇face = map( HalfEdges.polygons(topo) ) do abc
    a,b,c = abc
    pa,pb,pc = map(p->KVector(p, 1.0𝐼), P[abc])
    tri𝐼 = (pb-pa)∧(pc-pa); tri𝐼 *= 1.0/Multivectors.norm_sqr(tri𝐼);
    (u[b]*((pa-pc)⋅tri𝐼) + u[a]*((pc-pb)⋅tri𝐼) + u[c]*((pb-pa)⋅tri𝐼))
  end

  @test abs(sum(sum(vec∇face - coords.(∇face)))) < eps(1.0)*nvertices(topo)
  @test abs(sum(sum(vec∇face - coords.(ga∇face)))) < eps(1.0)*nvertices(topo)

  ## 3) find divergence of our vector field at vertices

  # go from dual triangle ( vertex in centre ) -> dual edge -> primal edge
  he_ = HalfEdges

  X = -map(normalize_safe, ∇face)

  vecX = map(coords, X)
  # vector algebra calculation
  vecdivX = map( map(h->OneRing(h, topo), he_.vertices(topo)) ) do vring
    vringⱼ = Iterators.filter(h->!he_.isboundary(topo,h), vring)

    cote(h) = ( fh = he_.face(topo, h);
                     Xⱼ = vecX[fh];
                     ea = P[he_.tail(topo, h)] - P[he_.head(topo, h)];
                     h = he_.next(topo, h);
                     eb = P[he_.head(topo, h)] - P[he_.tail(topo, h)];
                     ei = eb-ea;
                     cotθa = cotan(-eb, -ei);
                     cotθb = cotan(ei, -ea);
                     cotθa*(ea⋅Xⱼ) + cotθb*(eb⋅Xⱼ))

    0.5*mapreduce(cote, + ,vringⱼ)
  end

  # the distance function from our source point
  #vecϕ = Lc\vecdivX
  vecϕ = poisson(Lc, vecdivX, zeros(length(P)), [source])

  # with differential operators
  #
  # ∬𝑑ω = ⨚ω
  # each time you encounter an (evaluation of) exterior derivative of a form,
  # replace any evaluation over a simplex σ by a direct evaluation of the form itself over the boundary of σ
  # div F = ∇⋅F = ⋆𝑑⋆(F♭)

  #!me not quit working
  #==
  measuringform = ⋆(topo, P, 𝑑(topo))*array(P)
  valueform = 𝑑(topo, CellKForm{0}(ones(nfaces(topo))))*array(vecX)
  hecirculate = zeros(nhalfedges(topo));
  for (i,hh) in enumerate(UniqueHalfEdges(topo))
    hecirculate[hh] = measuringform[i,:]⋅valueform[i,:]
    hecirculate[he_.opposite(topo,hh)] = -hecirculate[hh]
  end
  decdivX = map(map(vh->OneRing(vh, topo), he_.vertices(topo))) do oner
    mapreduce(hh->hecirculate[hh], +, oner)
  end

  decϕ = poisson(Lc, decdivX, zeros(length(P)), [source])
  ==#

  #==
  sdf = ⋆(topo, P, 𝑑(topo, CellKForm{0}(ones(nfaces(topo)))))
  sdfX = array(sdf)*array(vecX)
  XX = SVector{3}[]
  resize!(XX, nhalfedges(topo))
  for (i,hh) in enumerate(UniqueHalfEdges(topo))
    XX[hh] = SVector{3}(sdfX[i,:])
    XX[he_.opposite(topo, hh)] = -XX[hh]
  end

  decX = map(map(vh->OneRing(vh, topo), he_.vertices(topo))) do oner
    mapreduce(hh->XX[hh], +, oner)
  end
  decdivX = sum.(decX)
  #decϕ = Lc\decdivX
  decϕ = poisson(Lc, decdivX, 10.0*ones(length(P)), [source])
  ==#

  # start with primal 1-form.
  # since it's a 1-form we will ultimately need ωᵢ⋅Xᵢ
  # this could be accomplished

  #== WIP
  gc = array(vecX)

  divX = ⋆(topo, P, 𝑑(topo, CellKForm{0}(gc[:,1])))
  divY = ⋆(topo, P, 𝑑(topo, CellKForm{0}(gc[:,2])))
  divZ = ⋆(topo, P, 𝑑(topo, CellKForm{0}(gc[:,3])))

  divG = array(divX+divY+divZ)*ones(nfaces(topo))

  divVert = zeros(nvertices(topo))
  for (iedge, uhe) in enumerate(UniqueHalfEdges(topo))
    b = head(topo, uhe)
    a = tail(topo, uhe)

    if a < b
      divVert[a] += divG[iedge]
      divVert[b] -= divG[iedge]
    else
      divVert[b] += divG[iedge]
      divVert[a] -= divG[iedge]
    end

  end

  b = -divVert
  ϕ = Lc\b
  ==#

    #!me need barycentric dual simplex for this to work ( otherwise ⋆₁⁻¹ has divide by zero )
  #== a 3,2,3,2,3 grid

  1   1---2---3
      |\ / \ /|
      | 4---5 |
      |/|\ /|\|
  0   6 | 7 | 8
      |\|/ \|/|
      | 9---A |
      |/ \ / \|
  -1  B---C---D
  y
  ↑→x -1  0   1
  ==#
  #==
  topo = Topology( [4,2,1,4,5,2,5,3,2,  6,4,1,7,5,4,5,8,3, 6,9,4, 9,7,4,7,10,5,10,8,5,
             11,9,6,9,10,7,10,13,8,   11,12,9,9,12,10,12,13,10] )

  @test iszero(𝑑(topo, 𝑑(topo))) == true

  P_flat_centred =
    SVector{3, Float64}.([[-1,1,0],[0,1,0],[1,1,0],
                         [-0.5,0.5,0],[0.5,0.5,0],
                         [-1,0,0],[0,0,0],[1,0,0],
                         [-0.5,-0.5,0],[0.5,-0.5,0],
                         [-1,-1,0],[0,-1,0],[1,-1,0]])
  P_flat_centred_2D =
    SVector{2, Float64}.([[-1,1],[0,1],[1,1],
                         [-0.5,0.5],[0.5,0.5],
                         [-1,0],[0,0],[1,0],
                         [-0.5,-0.5],[0.5,-0.5],
                         [-1,-1],[0,-1],[1,-1]])

  P = P_flat_centred_2D
  ==#

end

#!me cool example. need to work out some issues still
# kinda works but there appears to be some discretization artifacts visible on sphere.
function Δbarycentre(topo, P)
  he_ = HalfEdges
  Δrows = map(1:nvertices(topo)) do i
    oringᵢ = OneRing(topo, VertexHandle(i))
    map( oringᵢ ) do hⱼ
      vᵢ = P[he_.tail(topo, hⱼ)]
      a = P[he_.head(topo, hⱼ)]
      b = P[he_.head(topo, he_.next(topo, hⱼ))]

      vab = (vᵢ+a+b)/3.0

      ohⱼ = he_.opposite(topo, hⱼ)
      if !he_.isboundary(topo, ohⱼ)
        c = P[he_.head(topo, he_.next(topo, ohⱼ))]
        vca = (vᵢ+c+a)/3.0
      else
        vca = (vᵢ+a)/2.0
      end

      # a rotor taking objects in the tangent space of vᵢ from barycentre of ∠vab to barycentre of ∠vca
      #!me normalizing or any scaling needed here?
      rotor = (vab-vᵢ)/(vca-vᵢ)
      (he_.head(topo, hⱼ), rotor)
    end
  end

end

vector(k::KA) where {K<:KVector, KA<:AbstractVector{K}} = SVector{3}.(coords.(k))
function curvature_rotor(topo, Pga)
  Lbi = Δbarycentre(topo, Pga)
  Lbidot(i,P) = mapreduce( ((vj, rot),)-> rot*(Pga[vj]-Pga[i]) - (Pga[vj]-Pga[i]), +, Lbi[i])
  M = ⋆(topo, vector(Pga), PrimalKForm{0}(ones(nvertices(topo)))) |> array
  bicurve = inv(M)*map(i->Lbidot(i, Pga), 1:length(Pga))
  su = norm.(bicurve)
end

#==

  vector(k::KA) where {K<:KVector, KA<:AbstractVector{K}} = SVector{3}.(coords.(k))
  @generate_basis("+++")
  𝐼 = pseudoscalar(e₁)
  Pga = map(p->KVector(p, 1𝐼), (P))

  Lbi = Δbarycentre(topo, Pga)
  Lbidot(i,P) = mapreduce( ((vj, rot),)-> normalize(rot)*(Pga[i]-Pga[vj]), +, Lbi[i])
  M = ⋆(topo, vector(Pga), PrimalKForm{0}(ones(nvertices(topo)))) |> array
  bicurve = inv(M)*map(i->Lbidot(i, Pga), 1:length(Pga))
  su = norm.(bicurve)
  drawvertexx(topo, vector(Pga), su, [])

==#
#==
## Example: Exploring Discrete Laplace-Beltrami Operators

To illustrate the utillity of Multivectors as a tool for prototyping and exploring Geometric Algebra in CGI applications, we develop a Laplace-Beltrami operator using Multivector Types.

The core idea is that we can use the action of rotors on gradient vectors to compute their divergence in the neighbourhood of a vertex.  Mapping between the smooth setting and the discrete, each rotor on a face is conceptually operating in the tangent space generated by the differential of scalar function at the vertex.
From this argument we can see that, in the limit, all the rotors are acting in the same plane on gradient vectors in that same plane.  The important observation is that we don't need half angles or the sandwich operator to apply these tangent rotors.  The half-angle and sandwich operator are so the vector elements perpendicular to the plane of rotation are correctly handled.  We simply apply the rotors using the geometric product.  This makes the application rotors linear and simple to compute.

To compute the rotors for each incident edge ( triangle pair ) we just need the vertex, a start point and an end point.  The ratio of the two offset points is the desired rotor.  Note that the cotan formula used in the standard formulation can be interpreted as the aspect _ratio_ of integrated area contribution.

Now it is somewhat trivial to explore different discretizations.  All we need to do is compute the actual geometric points involved.  Circumcentre, orthocentre, barycentre, are all just points to our algorithm.

    function ΔGA(topo, P, cellvert)
      he_ = HalfEdges
      Δrows = map(1:nvertices(topo)) do i
        oringᵢ = OneRing(topo, VertexHandle(i))
        map( oringᵢ ) do hⱼ
          # points for triangle 1
          vᵢ = P[he_.tail(topo, hⱼ)]
          a  = P[he_.head(topo, hⱼ)]
          b_left = P[he_.head(topo, he_.next(topo, hⱼ))]

          # points for triangle 2
          ohⱼ = he_.opposite(topo, hⱼ)
          b_right = P[he_.head(topo, he_.next(topo, ohⱼ))]

          # point on triangle 1 (left of edge)
          v1 = cellvert(vᵢ, a, b_left)

          # point on triangle 2 (right of edge)
          v2 = cellvert(vᵢ, a, b_right)

          # a rotor taking objects in the tangent space of vᵢ from right triangle to left
          rotor = (v1-vᵢ)/(v2-vᵢ)

          # store the index of the vertex at end of edge and the rotor across the edge
          (he_.head(topo, hⱼ), rotor)
        end
      end
    end

    barycentre(v, a, b) = (v+a+b)/3.0
==#


text:
"""
joint_degree_seq(N::UnipartiteNetwork)
    N: Unipartite directed simple network
Returns the number of preys and predators for all species in N (i.e. the joint degree sequence)
"""
function joint_degree_seq(N::UnipartiteNetwork)
    # get the joint degree sequence
    kin = vec(sum(N.edges, dims=1)) # number of predators for all species (column sums)
    kout = vec(sum(N.edges, dims=2)) # number of preys for all species (row sums)
    k = (kin = kin, kout = kout) # combine in and out-degree sequences
    return k
end

ext: jl


text:
"""
`readtable(filepath)` read matlab mat or csv files, and returns a DataFrame
"""
function readtable(filepath)
    fname = basename(filepath);
	if occursin(r"(\.txt|\.csv)",fname)
		table_i = CSV.read(filepath, DataFrame);
	else
		dict_i = matopen(filepath) do io
			read(io)
		end # load xxx.mat and return a Dict
		table_i = DataFrame();
		for (k, v) in dict_i
			v2 = view(v, :); # a non-copied way of [v...]
			insertcols!(table_i, k => v2);
		end
	end

    return table_i
end


text:
function jl_CMGR_CreateCMgr()
    Dim = 100
    cmp_ref = Ref{CnstrMgrPointer}()
    ccall(
        (:CMGR_CreateCMgr, LIB_CVRPSEP),
        Cvoid,
        (Ptr{CnstrMgrPointer}, Cint),
        cmp_ref, Dim
    )
    return cmp_ref[]
end
init_CMP = jl_CMGR_CreateCMgr # alias

function jl_CMGR_MoveCnstr!(MyCutsCMP::CnstrMgrPointer, MyOldCutsCMP::CnstrMgrPointer)
    size = unsafe_load(MyCutsCMP).Size
    for i in 0:size-1
        ccall(
            (:CMGR_MoveCnstr, LIB_CVRPSEP),
            Cvoid,
            (CnstrMgrPointer, CnstrMgrPointer, Cint, Cint),
            MyCutsCMP, MyOldCutsCMP, i, 0
        )
    end
    unsafe_load(MyCutsCMP).Size = 0
end


text:
const VecTuple{N,T}= Union{NTuple{N,T}, Vector{NTuple{N,T}}}

# This is a bit inefficient. The more verbose implementation below is more efficient and
#   allows use of the internal expand!
#expand(I::Vector{NTuple{N,T}}) where {N,T}= [map(a->a[i], I) for i in 1:N]
expand(I::Vector{NTuple{N,T}}) where {N,T}= begin
    r= Vector{Vector{T}}(undef,N)
    expand!(r,I)
    return r
end
expand(I::NTuple{N,T}) where {N,T}= [I[i] for i in 1:N]
expand!(r,I::Vector{NTuple{N,T}}) where {N,T}= foreach(i-> (r[i]= map(a->a[i],I)), 1:N)
expand!(r,I::NTuple{N,T}) where {N,T}= foreach(i-> (r[i]= I[i]), 1:N)

# NOTE this could probably be a macro
# OPTIMIZE Less runtime cost?
joinTuples(x,y,z...)= (x..., joinTuples(y,z...)...)
joinTuples(x,y)= (x..., y...)
joinTuples(x)= x

cartesianIdx(iPre::VecTuple{N1,T},iPost::VecTuple{N2,T}) where {N1,N2,T<:Integer}= (
    vec( map( x->CartesianIndex(joinTuples(x...)), Iterators.product(iPre,iPost) ) )
        ::Vector{CartesianIndex{N1+N2}})

#using Lazy
#cartesianIdx(iPre::VecTuple{N1,T},iPost::VecTuple{N2,T}) where {N1,N2,T<:Integer}=
#    (@>> Iterators.product(iPre,iPost) begin
#      map( x->CartesianIndex(joinTuples(x...)) )
#      vec
#    end )::Vector{CartesianIndex{N1+N2}}

vecTuple_2_tupleVec(i::Vector{NTuple{N,Int}}) where N=
    NTuple{N,Vector{Int}}(expand(i))
vecTuple_2_tupleVec(i)= i


text:
@doc raw"""
    Line3

An object `l` of the data type [`Line3`](@ref) is a directed straight line in
the three-dimensional Euclidean space ``Ε^3``.
"""
Line3

"""
    Line3(p::Point3, q::Point3)

Introduces a line `l` passing through the points `p` and `q`.

Line `l` is directed from `p` to `q`.
"""
Line3(p::Point3, q::Point3)

"""
    Line3(p::Point3, d::Direction3)

Introduces a line `l` passing through point `p` with direction `d`.
"""
Line3(p::Point3, d::Direction3)

"""
    Line3(p::Point3, d::Vector3)

Introduces a line `l` passing through point `p` and oriented by `v`.
"""
Line3(p::Point3, v::Vector3)

"""
    Line3(s::Segment3)

Introduces a line `l` supporting the segment `s`, oriented from source to
target.
"""
Line3(s::Segment3)

"""
    Line3(r::Ray3)

Introduces a line `l` supporting the ray `r`, with same orientation.
"""
Line3(r::Ray3)

"""
    ==(l::Line3, h::Line3)

Test for equality: two lines are equal, iff they have a non empty intersection
and the same direction.
"""
==(l::Line3, h::Line3)

"""
    point(l::Line3, i::Real)

Returns an arbitrary point on `l`.

It holds `point(l, i) == point(l, j)`, iff `i == j`.
"""
point(l::Line3, i::Real)

@cxxdereference point(l::Line3, i::Real) = point(l, convert(FT, i))

"""
    projection(l::Line3, p::Point3)

Returns the orthogonal projection of `p` onto `l`.
"""
projection(l::Line3, p::Point3)

"""
    is_degenerate(l::Line3)

Returns `true` iff line `l` is degenerated to a point.
"""
is_degenerate(l::Line3)

"""
    has_on(l::Line3, p::Point3)

Returns `true`, iff `p` lies properly on `l`.
"""
has_on(l::Line3, p::Point3)

"""
    perpendicular_plane(l::Line3, p::Point3)

Returns the plane perpendicular to `l` passing through `p`.
"""
perpendicular_plane(l::Line3, p::Point3)

"""
    to_vector(l::Line3)

Returns a vector having the same direction as `l`.
"""
to_vector(l::Line3)

"""
    direction(l::Line3)

Returns the direction of `l`.
"""
direction(l::Line3)

"""
    opposite(l::Line3)

Returns the line with opposite direction.
"""
opposite(l::Line3)

"""
    transform(l::Line3, t::AffTransformation2)

Returns the line obtained by applying `t` on a point on `l` and the direction of
`l`.
"""
transform(l::Line3, t::AffTransformation2)


text:
exact_Yᵏs = Dict(((1,0), (1,0), 0) => r -> -1.0*r*exp(-2*r) + 1.0 - 1.0*exp(-2*r),
                 ((2,0), (1,0), 0) => r -> 0.314269680527354*r^2*exp(-3*r/2) + 0.20951312035157*r*exp(-3*r/2),
                 ((3,0), (1,0), 0) => r -> -0.0320750149549792*r^3*exp(-4*r/3) + 0.144337567297406*r^2*exp(-4*r/3) + 0.108253175473055*r*exp(-4*r/3),
                 ((4,0), (1,0), 0) => r -> 0.00166666666666667*r^4*exp(-5*r/4) - 0.0293333333333333*r^3*exp(-5*r/4) + 0.0864*r^2*exp(-5*r/4) + 0.06912*r*exp(-5*r/4),
                 ((1,0), (2,0), 0) => r -> 0.314269680527354*r^2*exp(-3*r/2) + 0.20951312035157*r*exp(-3*r/2),
                 ((2,0), (2,0), 0) => r -> -0.125*r^3*exp(-r) - 0.25*r^2*exp(-r) - 0.75*r*exp(-r) + 1.0 - 1.0*exp(-r),
                 ((3,0), (2,0), 0) => r -> 0.0145154947720485*r^4*exp(-5*r/6) - 0.0203216926808679*r^3*exp(-5*r/6) + 0.0600941483562806*r^2*exp(-5*r/6) + 0.0721129780275368*r*exp(-5*r/6),
                 ((4,0), (2,0), 0) => r -> -0.000818410626373319*r^5*exp(-3*r/4) + 0.0103665346007287*r^4*exp(-3*r/4) - 0.0174594266959641*r^3*exp(-3*r/4) + 0.0320089489426009*r^2*exp(-3*r/4) + 0.0426785985901346*r*exp(-3*r/4),
                 ((2,1), (2,1), 0) => r -> -0.0416666666666667*r^3*exp(-r) - 0.25*r^2*exp(-r) - 0.75*r*exp(-r) + 1.0 - 1.0*exp(-r),
                 ((3,1), (2,1), 0) => r -> 0.00592592592592593*r^4*exp(-5*r/6) + 0.0213333333333333*r^3*exp(-5*r/6) + 0.0512*r^2*exp(-5*r/6) + 0.06144*r*exp(-5*r/6),
                 ((4,1), (2,1), 0) => r -> -0.000366004358815785*r^5*exp(-3*r/4) + 0.0024400290587719*r^4*exp(-3*r/4) + 0.00976011623508759*r^3*exp(-3*r/4) + 0.0260269766269002*r^2*exp(-3*r/4) + 0.0347026355025336*r*exp(-3*r/4),
                 ((1,0), (3,0), 0) => r -> -0.0320750149549792*r^3*exp(-4*r/3) + 0.144337567297406*r^2*exp(-4*r/3) + 0.108253175473055*r*exp(-4*r/3),
                 ((2,0), (3,0), 0) => r -> 0.0145154947720485*r^4*exp(-5*r/6) - 0.0203216926808679*r^3*exp(-5*r/6) + 0.0600941483562806*r^2*exp(-5*r/6) + 0.0721129780275368*r*exp(-5*r/6),
                 ((3,0), (3,0), 0) => r -> -0.00182898948331047*r^5*exp(-2*r/3) + 0.00548696844993141*r^4*exp(-2*r/3) - 0.0493827160493827*r^3*exp(-2*r/3) - 0.148148148148148*r^2*exp(-2*r/3) - 0.555555555555556*r*exp(-2*r/3) + 1.0 - 1.0*exp(-2*r/3),
                 ((4,0), (3,0), 0) => r -> 0.000109098690323059*r^6*exp(-7*r/12) - 0.00135594086544373*r^5*exp(-7*r/12) + 0.00788516350569616*r^4*exp(-7*r/12) - 0.00985132547138708*r^3*exp(-7*r/12) + 0.0212097668780658*r^2*exp(-7*r/12) + 0.0363596003623985*r*exp(-7*r/12),
                 ((2,1), (3,1), 0) => r -> 0.00592592592592593*r^4*exp(-5*r/6) + 0.0213333333333333*r^3*exp(-5*r/6) + 0.0512*r^2*exp(-5*r/6) + 0.06144*r*exp(-5*r/6),
                 ((3,1), (3,1), 0) => r -> -0.000914494741655236*r^5*exp(-2*r/3) - 0.00274348422496571*r^4*exp(-2*r/3) - 0.0246913580246914*r^3*exp(-2*r/3) - 0.148148148148148*r^2*exp(-2*r/3) - 0.555555555555556*r*exp(-2*r/3) + 1.0 - 1.0*exp(-2*r/3),
                 ((4,1), (3,1), 0) => r -> 5.97558136842097e-5*r^6*exp(-7*r/12) - 0.000324388702857139*r^5*exp(-7*r/12) + 0.00112194588958108*r^4*exp(-7*r/12) + 0.00577000743213127*r^3*exp(-7*r/12) + 0.0197828826244501*r^2*exp(-7*r/12) + 0.0339135130704858*r*exp(-7*r/12),
                 ((3,2), (3,2), 0) => r -> -0.000182898948331047*r^5*exp(-2*r/3) - 0.00274348422496571*r^4*exp(-2*r/3) - 0.0246913580246914*r^3*exp(-2*r/3) - 0.148148148148148*r^2*exp(-2*r/3) - 0.555555555555555*r*exp(-2*r/3) + 1.0 - 1.0*exp(-2*r/3),
                 ((4,2), (3,2), 0) => r -> 1.54288847492012e-5*r^6*exp(-7*r/12) + 0.000132247583564581*r^5*exp(-7*r/12) + 0.000906840573014273*r^4*exp(-7*r/12) + 0.00466375151835912*r^3*exp(-7*r/12) + 0.0159900052058027*r^2*exp(-7*r/12) + 0.0274114374956618*r*exp(-7*r/12),
                 ((1,0), (4,0), 0) => r -> 0.00166666666666667*r^4*exp(-5*r/4) - 0.0293333333333333*r^3*exp(-5*r/4) + 0.0864*r^2*exp(-5*r/4) + 0.06912*r*exp(-5*r/4),
                 ((2,0), (4,0), 0) => r -> -0.000818410626373319*r^5*exp(-3*r/4) + 0.0103665346007287*r^4*exp(-3*r/4) - 0.0174594266959641*r^3*exp(-3*r/4) + 0.0320089489426009*r^2*exp(-3*r/4) + 0.0426785985901346*r*exp(-3*r/4),
                 ((3,0), (4,0), 0) => r -> 0.000109098690323059*r^6*exp(-7*r/12) - 0.00135594086544373*r^5*exp(-7*r/12) + 0.00788516350569616*r^4*exp(-7*r/12) - 0.00985132547138708*r^3*exp(-7*r/12) + 0.0212097668780658*r^2*exp(-7*r/12) + 0.0363596003623985*r*exp(-7*r/12),
                 ((4,0), (4,0), 0) => r -> -6.78168402777778e-6*r^7*exp(-r/2) + 0.000135633680555556*r^6*exp(-r/2) - 0.00146484375*r^5*exp(-r/2) + 0.00390625*r^4*exp(-r/2) - 0.0234375*r^3*exp(-r/2) - 0.09375*r^2*exp(-r/2) - 0.4375*r*exp(-r/2) + 1.0 - 1.0*exp(-r/2),
                 ((2,1), (4,1), 0) => r -> -0.000366004358815785*r^5*exp(-3*r/4) + 0.0024400290587719*r^4*exp(-3*r/4) + 0.00976011623508759*r^3*exp(-3*r/4) + 0.0260269766269002*r^2*exp(-3*r/4) + 0.0347026355025336*r*exp(-3*r/4),
                 ((3,1), (4,1), 0) => r -> 5.97558136842097e-5*r^6*exp(-7*r/12) - 0.000324388702857139*r^5*exp(-7*r/12) + 0.00112194588958108*r^4*exp(-7*r/12) + 0.00577000743213127*r^3*exp(-7*r/12) + 0.0197828826244501*r^2*exp(-7*r/12) + 0.0339135130704858*r*exp(-7*r/12),
                 ((4,1), (4,1), 0) => r -> -4.06901041666667e-6*r^7*exp(-r/2) + 4.8828125e-5*r^6*exp(-r/2) - 0.000423177083333333*r^5*exp(-r/2) - 0.00130208333333333*r^4*exp(-r/2) - 0.0130208333333333*r^3*exp(-r/2) - 0.09375*r^2*exp(-r/2) - 0.4375*r*exp(-r/2) + 1.0 - 1.0*exp(-r/2),
                 ((3,2), (4,2), 0) => r -> 1.54288847492012e-5*r^6*exp(-7*r/12) + 0.000132247583564581*r^5*exp(-7*r/12) + 0.000906840573014273*r^4*exp(-7*r/12) + 0.00466375151835912*r^3*exp(-7*r/12) + 0.0159900052058027*r^2*exp(-7*r/12) + 0.0274114374956618*r*exp(-7*r/12),
                 ((4,2), (4,2), 0) => r -> -1.35633680555556e-6*r^7*exp(-r/2) - 5.42534722222222e-6*r^6*exp(-r/2) - 9.765625e-5*r^5*exp(-r/2) - 0.00130208333333333*r^4*exp(-r/2) - 0.0130208333333333*r^3*exp(-r/2) - 0.09375*r^2*exp(-r/2) - 0.4375*r*exp(-r/2) + 1.0 - 1.0*exp(-r/2),
                 ((4,3), (4,3), 0) => r -> -1.93762400793651e-7*r^7*exp(-r/2) - 5.42534722222222e-6*r^6*exp(-r/2) - 9.765625e-5*r^5*exp(-r/2) - 0.00130208333333333*r^4*exp(-r/2) - 0.0130208333333333*r^3*exp(-r/2) - 0.09375*r^2*exp(-r/2) - 0.4375*r*exp(-r/2) + 1.0 - 1.0*exp(-r/2),
                 ((2,1), (1,0), 1) => r -> -0.544331053951817*r^2*exp(-3*r/2) - 1.45154947720485*r*exp(-3*r/2) - 1.93539930293979*exp(-3*r/2) + 1.29026620195986/r - 1.29026620195986*exp(-3*r/2)/r,
                 ((3,1), (1,0), 1) => r -> 0.0680413817439772*r^3*exp(-4*r/3) - 0.102062072615966*r^2*exp(-4*r/3) - 0.459279326771846*r*exp(-4*r/3) - 0.688918990157769*exp(-4*r/3) + 0.516689242618327/r - 0.516689242618327*exp(-4*r/3)/r,
                 ((4,1), (1,0), 1) => r -> -0.00387298334620742*r^4*exp(-5*r/4) + 0.0526725735084209*r^3*exp(-5*r/4) - 0.0322232214404457*r^2*exp(-5*r/4) - 0.237956096790984*r*exp(-5*r/4) - 0.380729754865574*exp(-5*r/4) + 0.304583803892459/r - 0.304583803892459*exp(-5*r/4)/r,
                 ((2,1), (2,0), 1) => r -> 0.21650635094611*r^3*exp(-r) + 0.866025403784439*r^2*exp(-r) + 2.59807621135332*r*exp(-r) + 5.19615242270663*exp(-r) - 5.19615242270663/r + 5.19615242270663*exp(-r)/r,
                 ((3,1), (2,0), 1) => r -> -0.03079201435678*r^4*exp(-5*r/6) - 0.0492672229708481*r^3*exp(-5*r/6) - 0.280823170933834*r^2*exp(-5*r/6) - 1.06417201617032*r*exp(-5*r/6) - 2.55401283880876*exp(-5*r/6) + 3.06481540657052/r - 3.06481540657052*exp(-5*r/6)/r,
                 ((4,1), (2,0), 1) => r -> 0.00190181443578183*r^5*exp(-3*r/4) - 0.0164823917767758*r^4*exp(-3*r/4) - 0.00845250860347478*r^3*exp(-3*r/4) - 0.078890080299098*r^2*exp(-3*r/4) - 0.360640367081591*r*exp(-3*r/4) - 0.961707645550909*exp(-3*r/4) + 1.28227686073455/r - 1.28227686073455*exp(-3*r/4)/r,
                 ((1,0), (2,1), 1) => r -> -0.544331053951817*r^2*exp(-3*r/2) - 1.45154947720485*r*exp(-3*r/2) - 1.93539930293979*exp(-3*r/2) + 1.29026620195986/r - 1.29026620195986*exp(-3*r/2)/r,
                 ((2,0), (2,1), 1) => r -> 0.21650635094611*r^3*exp(-r) + 0.866025403784439*r^2*exp(-r) + 2.59807621135332*r*exp(-r) + 5.19615242270663*exp(-r) - 5.19615242270663/r + 5.19615242270663*exp(-r)/r,
                 ((3,0), (2,1), 1) => r -> -0.0251415744421884*r^4*exp(-5*r/6) - 0.015084944665313*r^3*exp(-5*r/6) - 0.0859841845922842*r^2*exp(-5*r/6) - 0.325834804770761*r*exp(-5*r/6) - 0.782003531449827*exp(-5*r/6) + 0.938404237739792/r - 0.938404237739792*exp(-5*r/6)/r,
                 ((3,2), (2,1), 1) => r -> -0.00795046391999925*r^4*exp(-5*r/6) - 0.0763244536319928*r^3*exp(-5*r/6) - 0.435049385702359*r^2*exp(-5*r/6) - 1.64860819845104*r*exp(-5*r/6) - 3.95665967628251*exp(-5*r/6) + 4.74799161153901/r - 4.74799161153901*exp(-5*r/6)/r,
                 ((4,0), (2,1), 1) => r -> 0.00141752878633286*r^5*exp(-3*r/4) - 0.0151203070542171*r^4*exp(-3*r/4) - 0.00252005117570286*r^3*exp(-3*r/4) - 0.0235204776398933*r^2*exp(-3*r/4) - 0.107522183496655*r*exp(-3*r/4) - 0.286725822657747*exp(-3*r/4) + 0.382301096876997/r - 0.382301096876997*exp(-3*r/4)/r,
                 ((4,2), (2,1), 1) => r -> 0.000633938145260609*r^5*exp(-3*r/4) + 0.000845250860347479*r^4*exp(-3*r/4) - 0.0112700114712997*r^3*exp(-3*r/4) - 0.105186773732131*r^2*exp(-3*r/4) - 0.480853822775454*r*exp(-3*r/4) - 1.28227686073455*exp(-3*r/4) + 1.70970248097939/r - 1.70970248097939*exp(-3*r/4)/r,
                 ((2,1), (3,0), 1) => r -> -0.0251415744421884*r^4*exp(-5*r/6) - 0.015084944665313*r^3*exp(-5*r/6) - 0.0859841845922842*r^2*exp(-5*r/6) - 0.325834804770761*r*exp(-5*r/6) - 0.782003531449827*exp(-5*r/6) + 0.938404237739792/r - 0.938404237739792*exp(-5*r/6)/r,
                 ((3,1), (3,0), 1) => r -> 0.00387987259910314*r^5*exp(-2*r/3) + 0.104756560175785*r^3*exp(-2*r/3) + 0.628539361054709*r^2*exp(-2*r/3) + 2.82842712474619*r*exp(-2*r/3) + 8.48528137423857*exp(-2*r/3) - 12.7279220613579/r + 12.7279220613579*exp(-2*r/3)/r,
                 ((4,1), (3,0), 1) => r -> -0.000253522446428548*r^6*exp(-7*r/12) + 0.00213683204846919*r^5*exp(-7*r/12) - 0.0118457069611869*r^4*exp(-7*r/12) - 0.0245022422714471*r^3*exp(-7*r/12) - 0.177708570320386*r^2*exp(-7*r/12) - 0.930546695495837*r*exp(-7*r/12) - 3.19044581312858*exp(-7*r/12) + 5.469335679649/r - 5.469335679649*exp(-7*r/12)/r,
                 ((1,0), (3,1), 1) => r -> 0.0680413817439772*r^3*exp(-4*r/3) - 0.102062072615966*r^2*exp(-4*r/3) - 0.459279326771846*r*exp(-4*r/3) - 0.688918990157769*exp(-4*r/3) + 0.516689242618327/r - 0.516689242618327*exp(-4*r/3)/r,
                 ((2,0), (3,1), 1) => r -> -0.03079201435678*r^4*exp(-5*r/6) - 0.0492672229708481*r^3*exp(-5*r/6) - 0.280823170933834*r^2*exp(-5*r/6) - 1.06417201617032*r*exp(-5*r/6) - 2.55401283880876*exp(-5*r/6) + 3.06481540657052/r - 3.06481540657052*exp(-5*r/6)/r,
                 ((3,0), (3,1), 1) => r -> 0.00387987259910314*r^5*exp(-2*r/3) + 0.104756560175785*r^3*exp(-2*r/3) + 0.628539361054709*r^2*exp(-2*r/3) + 2.82842712474619*r*exp(-2*r/3) + 8.48528137423857*exp(-2*r/3) - 12.7279220613579/r + 12.7279220613579*exp(-2*r/3)/r,
                 ((3,2), (3,1), 1) => r -> 0.00122692344444433*r^5*exp(-2*r/3) + 0.011042310999999*r^4*exp(-2*r/3) + 0.0828173324999922*r^3*exp(-2*r/3) + 0.496903994999953*r^2*exp(-2*r/3) + 2.23606797749979*r*exp(-2*r/3) + 6.70820393249937*exp(-2*r/3) - 10.0623058987491/r + 10.0623058987491*exp(-2*r/3)/r,
                 ((4,0), (3,1), 1) => r -> -0.000231433271238018*r^6*exp(-7*r/12) + 0.00218208512881559*r^5*exp(-7*r/12) - 0.0102019564464106*r^4*exp(-7*r/12) - 0.0109468612028152*r^3*exp(-7*r/12) - 0.0793948175149231*r^2*exp(-7*r/12) - 0.41574013535087*r*exp(-7*r/12) - 1.42539474977441*exp(-7*r/12) + 2.44353385675613/r - 2.44353385675613*exp(-7*r/12)/r,
                 ((4,2), (3,1), 1) => r -> -0.000103500105348671*r^6*exp(-7*r/12) - 0.000266143128039439*r^5*exp(-7*r/12) - 0.00349788111137549*r^4*exp(-7*r/12) - 0.0338925126319613*r^3*exp(-7*r/12) - 0.245813827880159*r^2*exp(-7*r/12) - 1.28717058962701*r*exp(-7*r/12) - 4.41315630729261*exp(-7*r/12) + 7.56541081250162/r - 7.56541081250162*exp(-7*r/12)/r,
                 ((2,1), (3,2), 1) => r -> -0.00795046391999925*r^4*exp(-5*r/6) - 0.0763244536319928*r^3*exp(-5*r/6) - 0.435049385702359*r^2*exp(-5*r/6) - 1.64860819845104*r*exp(-5*r/6) - 3.95665967628251*exp(-5*r/6) + 4.74799161153901/r - 4.74799161153901*exp(-5*r/6)/r,
                 ((3,1), (3,2), 1) => r -> 0.00122692344444433*r^5*exp(-2*r/3) + 0.011042310999999*r^4*exp(-2*r/3) + 0.0828173324999922*r^3*exp(-2*r/3) + 0.496903994999953*r^2*exp(-2*r/3) + 2.23606797749979*r*exp(-2*r/3) + 6.70820393249937*exp(-2*r/3) - 10.0623058987491/r + 10.0623058987491*exp(-2*r/3)/r,
                 ((4,1), (3,2), 1) => r -> -8.01708368692231e-5*r^6*exp(-7*r/12) - 4.58119067824132e-5*r^5*exp(-7*r/12) - 0.000602099346283145*r^4*exp(-7*r/12) - 0.00583400608820935*r^3*exp(-7*r/12) - 0.0423125716287711*r^2*exp(-7*r/12) - 0.221564011437929*r*exp(-7*r/12) - 0.759648039215756*exp(-7*r/12) + 1.30225378151272/r - 1.30225378151272*exp(-7*r/12)/r,
                 ((4,3), (3,2), 1) => r -> -1.74947108800558e-5*r^6*exp(-7*r/12) - 0.000359891195246862*r^5*exp(-7*r/12) - 0.00472999856610161*r^4*exp(-7*r/12) - 0.0458310420069473*r^3*exp(-7*r/12) - 0.332400964006431*r^2*exp(-7*r/12) - 1.74057232061549*r*exp(-7*r/12) - 5.96767652782454*exp(-7*r/12) + 10.2303026191278/r - 10.2303026191278*exp(-7*r/12)/r,
                 ((2,1), (4,0), 1) => r -> 0.00141752878633286*r^5*exp(-3*r/4) - 0.0151203070542171*r^4*exp(-3*r/4) - 0.00252005117570286*r^3*exp(-3*r/4) - 0.0235204776398933*r^2*exp(-3*r/4) - 0.107522183496655*r*exp(-3*r/4) - 0.286725822657747*exp(-3*r/4) + 0.382301096876997/r - 0.382301096876997*exp(-3*r/4)/r,
                 ((3,1), (4,0), 1) => r -> -0.000231433271238018*r^6*exp(-7*r/12) + 0.00218208512881559*r^5*exp(-7*r/12) - 0.0102019564464106*r^4*exp(-7*r/12) - 0.0109468612028152*r^3*exp(-7*r/12) - 0.0793948175149231*r^2*exp(-7*r/12) - 0.41574013535087*r*exp(-7*r/12) - 1.42539474977441*exp(-7*r/12) + 2.44353385675613/r - 2.44353385675613*exp(-7*r/12)/r,
                 ((4,1), (4,0), 1) => r -> 1.57592095792945e-5*r^7*exp(-r/2) - 0.000252147353268712*r^6*exp(-r/2) + 0.00252147353268712*r^5*exp(-r/2) + 0.0605153647844909*r^3*exp(-r/2) + 0.484122918275927*r^2*exp(-r/2) + 2.90473750965556*r*exp(-r/2) + 11.6189500386223*exp(-r/2) - 23.2379000772445/r + 23.2379000772445*exp(-r/2)/r,
                 ((1,0), (4,1), 1) => r -> -0.00387298334620742*r^4*exp(-5*r/4) + 0.0526725735084209*r^3*exp(-5*r/4) - 0.0322232214404457*r^2*exp(-5*r/4) - 0.237956096790984*r*exp(-5*r/4) - 0.380729754865574*exp(-5*r/4) + 0.304583803892459/r - 0.304583803892459*exp(-5*r/4)/r,
                 ((2,0), (4,1), 1) => r -> 0.00190181443578183*r^5*exp(-3*r/4) - 0.0164823917767758*r^4*exp(-3*r/4) - 0.00845250860347478*r^3*exp(-3*r/4) - 0.078890080299098*r^2*exp(-3*r/4) - 0.360640367081591*r*exp(-3*r/4) - 0.961707645550909*exp(-3*r/4) + 1.28227686073455/r - 1.28227686073455*exp(-3*r/4)/r,
                 ((3,0), (4,1), 1) => r -> -0.000253522446428548*r^6*exp(-7*r/12) + 0.00213683204846919*r^5*exp(-7*r/12) - 0.0118457069611869*r^4*exp(-7*r/12) - 0.0245022422714471*r^3*exp(-7*r/12) - 0.177708570320386*r^2*exp(-7*r/12) - 0.930546695495837*r*exp(-7*r/12) - 3.19044581312858*exp(-7*r/12) + 5.469335679649/r - 5.469335679649*exp(-7*r/12)/r,
                 ((3,2), (4,1), 1) => r -> -8.01708368692231e-5*r^6*exp(-7*r/12) - 4.58119067824132e-5*r^5*exp(-7*r/12) - 0.000602099346283145*r^4*exp(-7*r/12) - 0.00583400608820935*r^3*exp(-7*r/12) - 0.0423125716287711*r^2*exp(-7*r/12) - 0.221564011437929*r*exp(-7*r/12) - 0.759648039215756*exp(-7*r/12) + 1.30225378151272/r - 1.30225378151272*exp(-7*r/12)/r,
                 ((4,0), (4,1), 1) => r -> 1.57592095792945e-5*r^7*exp(-r/2) - 0.000252147353268712*r^6*exp(-r/2) + 0.00252147353268712*r^5*exp(-r/2) + 0.0605153647844909*r^3*exp(-r/2) + 0.484122918275927*r^2*exp(-r/2) + 2.90473750965556*r*exp(-r/2) + 11.6189500386223*exp(-r/2) - 23.2379000772445/r + 23.2379000772445*exp(-r/2)/r,
                 ((4,2), (4,1), 1) => r -> 7.04773277819367e-6*r^7*exp(-r/2) - 2.81909311127747e-5*r^6*exp(-r/2) + 0.000451054897804395*r^5*exp(-r/2) + 0.00541265877365274*r^4*exp(-r/2) + 0.0541265877365274*r^3*exp(-r/2) + 0.433012701892219*r^2*exp(-r/2) + 2.59807621135332*r*exp(-r/2) + 10.3923048454133*exp(-r/2) - 20.7846096908265/r + 20.7846096908265*exp(-r/2)/r,
                 ((2,1), (4,2), 1) => r -> 0.000633938145260609*r^5*exp(-3*r/4) + 0.000845250860347479*r^4*exp(-3*r/4) - 0.0112700114712997*r^3*exp(-3*r/4) - 0.105186773732131*r^2*exp(-3*r/4) - 0.480853822775454*r*exp(-3*r/4) - 1.28227686073455*exp(-3*r/4) + 1.70970248097939/r - 1.70970248097939*exp(-3*r/4)/r,
                 ((3,1), (4,2), 1) => r -> -0.000103500105348671*r^6*exp(-7*r/12) - 0.000266143128039439*r^5*exp(-7*r/12) - 0.00349788111137549*r^4*exp(-7*r/12) - 0.0338925126319613*r^3*exp(-7*r/12) - 0.245813827880159*r^2*exp(-7*r/12) - 1.28717058962701*r*exp(-7*r/12) - 4.41315630729261*exp(-7*r/12) + 7.56541081250162/r - 7.56541081250162*exp(-7*r/12)/r,
                 ((4,1), (4,2), 1) => r -> 7.04773277819367e-6*r^7*exp(-r/2) - 2.81909311127747e-5*r^6*exp(-r/2) + 0.000451054897804395*r^5*exp(-r/2) + 0.00541265877365274*r^4*exp(-r/2) + 0.0541265877365274*r^3*exp(-r/2) + 0.433012701892219*r^2*exp(-r/2) + 2.59807621135332*r*exp(-r/2) + 10.3923048454133*exp(-r/2) - 20.7846096908265/r + 20.7846096908265*exp(-r/2)/r,
                 ((4,3), (4,2), 1) => r -> 1.53794137780447e-6*r^7*exp(-r/2) + 2.46070620448716e-5*r^6*exp(-r/2) + 0.000344498868628202*r^5*exp(-r/2) + 0.00413398642353842*r^4*exp(-r/2) + 0.0413398642353842*r^3*exp(-r/2) + 0.330718913883074*r^2*exp(-r/2) + 1.98431348329844*r*exp(-r/2) + 7.93725393319377*exp(-r/2) - 15.8745078663875/r + 15.8745078663875*exp(-r/2)/r,
                 ((3,2), (4,3), 1) => r -> -1.74947108800558e-5*r^6*exp(-7*r/12) - 0.000359891195246862*r^5*exp(-7*r/12) - 0.00472999856610161*r^4*exp(-7*r/12) - 0.0458310420069473*r^3*exp(-7*r/12) - 0.332400964006431*r^2*exp(-7*r/12) - 1.74057232061549*r*exp(-7*r/12) - 5.96767652782454*exp(-7*r/12) + 10.2303026191278/r - 10.2303026191278*exp(-7*r/12)/r,
                 ((4,2), (4,3), 1) => r -> 1.53794137780447e-6*r^7*exp(-r/2) + 2.46070620448716e-5*r^6*exp(-r/2) + 0.000344498868628202*r^5*exp(-r/2) + 0.00413398642353842*r^4*exp(-r/2) + 0.0413398642353842*r^3*exp(-r/2) + 0.330718913883074*r^2*exp(-r/2) + 1.98431348329844*r*exp(-r/2) + 7.93725393319377*exp(-r/2) - 15.8745078663875/r + 15.8745078663875*exp(-r/2)/r,
                 ((3,2), (1,0), 2) => r -> -0.0507150516208487*r^3*exp(-4*r/3) - 0.228217732293819*r^2*exp(-4*r/3) - 0.684653196881458*r*exp(-4*r/3) - 1.54046969298328*exp(-4*r/3) - 2.31070453947492*exp(-4*r/3)/r + 1.73302840460619/r^2 - 1.73302840460619*exp(-4*r/3)/r^2,
                 ((4,2), (1,0), 2) => r -> 0.00372677996249965*r^4*exp(-5*r/4) - 0.020869967789998*r^3*exp(-5*r/4) - 0.114486680447989*r^2*exp(-5*r/4) - 0.366357377433566*r*exp(-5*r/4) - 0.879257705840557*exp(-5*r/4) - 1.40681232934489*exp(-5*r/4)/r + 1.12544986347591/r^2 - 1.12544986347591*exp(-5*r/4)/r^2,
                 ((3,2), (2,0), 2) => r -> 0.0229510124219699*r^4*exp(-5*r/6) + 0.174427694406971*r^3*exp(-5*r/6) + 1.05758265240437*r^2*exp(-5*r/6) + 5.07639673154099*r*exp(-5*r/6) + 18.2750282335475*exp(-5*r/6) + 43.8600677605141*exp(-5*r/6)/r - 52.6320813126169/r^2 + 52.6320813126169*exp(-5*r/6)/r^2,
                 ((4,2), (2,0), 2) => r -> -0.00183002179407892*r^5*exp(-3*r/4) + 0.00122001452938595*r^4*exp(-3*r/4) + 0.0146401743526314*r^3*exp(-3*r/4) + 0.104107906507601*r^2*exp(-3*r/4) + 0.555242168040538*r*exp(-3*r/4) + 2.22096867216215*exp(-3*r/4) + 5.92258312576574*exp(-3*r/4)/r - 7.89677750102099/r^2 + 7.89677750102099*exp(-3*r/4)/r^2,
                 ((2,1), (2,1), 2) => r -> -0.208333333333333*r^3*exp(-r) - 1.25*r^2*exp(-r) - 5.0*r*exp(-r) - 15.0*exp(-r) - 30.0*exp(-r)/r + 30.0/r^2 - 30.0*exp(-r)/r^2,
                 ((3,1), (2,1), 2) => r -> 0.0296296296296296*r^4*exp(-5*r/6) + 0.106666666666667*r^3*exp(-5*r/6) + 0.512*r^2*exp(-5*r/6) + 2.4576*r*exp(-5*r/6) + 8.84736*exp(-5*r/6) + 21.233664*exp(-5*r/6)/r - 25.4803968/r^2 + 25.4803968*exp(-5*r/6)/r^2,
                 ((4,1), (2,1), 2) => r -> -0.00183002179407892*r^5*exp(-3*r/4) + 0.0122001452938595*r^4*exp(-3*r/4) + 0.0292803487052628*r^3*exp(-3*r/4) + 0.104107906507601*r^2*exp(-3*r/4) + 0.555242168040538*r*exp(-3*r/4) + 2.22096867216215*exp(-3*r/4) + 5.92258312576574*exp(-3*r/4)/r - 7.89677750102099/r^2 + 7.89677750102099*exp(-3*r/4)/r^2,
                 ((4,3), (2,1), 2) => r -> -0.000399343494989788*r^5*exp(-3*r/4) - 0.00532457993319718*r^4*exp(-3*r/4) - 0.0468563034121352*r^3*exp(-3*r/4) - 0.318054908009645*r^2*exp(-3*r/4) - 1.69629284271811*r*exp(-3*r/4) - 6.78517137087242*exp(-3*r/4) - 18.0937903223265*exp(-3*r/4)/r + 24.1250537631019/r^2 - 24.1250537631019*exp(-3*r/4)/r^2,
                 ((3,2), (3,0), 2) => r -> -0.0028918862918778*r^5*exp(-2*r/3) - 0.0173513177512668*r^4*exp(-2*r/3) - 0.156161859761401*r^3*exp(-2*r/3) - 1.17121394821051*r^2*exp(-2*r/3) - 7.02728368926307*r*exp(-2*r/3) - 31.6227766016838*exp(-2*r/3) - 94.8683298050514*exp(-2*r/3)/r + 142.302494707577/r^2 - 142.302494707577*exp(-2*r/3)/r^2,
                 ((4,2), (3,0), 2) => r -> 0.000243952087818557*r^6*exp(-7*r/12) - 0.000104550894779382*r^5*exp(-7*r/12) + 0.0106417875043299*r^4*exp(-7*r/12) + 0.110725798646556*r^3*exp(-7*r/12) + 0.951250981920879*r^2*exp(-7*r/12) + 6.52286387602889*r*exp(-7*r/12) + 33.54615707672*exp(-7*r/12) + 115.015395691611*exp(-7*r/12)/r - 197.169249757048/r^2 + 197.169249757048*exp(-7*r/12)/r^2,
                 ((2,1), (3,1), 2) => r -> 0.0296296296296296*r^4*exp(-5*r/6) + 0.106666666666667*r^3*exp(-5*r/6) + 0.512*r^2*exp(-5*r/6) + 2.4576*r*exp(-5*r/6) + 8.84736*exp(-5*r/6) + 21.233664*exp(-5*r/6)/r - 25.4803968/r^2 + 25.4803968*exp(-5*r/6)/r^2,
                 ((3,1), (3,1), 2) => r -> -0.00457247370827618*r^5*exp(-2*r/3) - 0.0137174211248285*r^4*exp(-2*r/3) - 0.185185185185185*r^3*exp(-2*r/3) - 1.48148148148148*r^2*exp(-2*r/3) - 8.88888888888889*r*exp(-2*r/3) - 40.0*exp(-2*r/3) - 120.0*exp(-2*r/3)/r + 180.0/r^2 - 180.0*exp(-2*r/3)/r^2,
                 ((4,1), (3,1), 2) => r -> 0.000298779068421049*r^6*exp(-7*r/12) - 0.00162194351428569*r^5*exp(-7*r/12) + 0.0108779971033296*r^4*exp(-7*r/12) + 0.0725013977341712*r^3*exp(-7*r/12) + 0.553920713484602*r^2*exp(-7*r/12) + 3.79831346389441*r*exp(-7*r/12) + 19.5341835285998*exp(-7*r/12) + 66.974343526628*exp(-7*r/12)/r - 114.813160331362/r^2 + 114.813160331362*exp(-7*r/12)/r^2,
                 ((4,3), (3,1), 2) => r -> 6.51989379575165e-5*r^6*exp(-7*r/12) + 0.000950041667380954*r^5*exp(-7*r/12) + 0.0116879075718127*r^4*exp(-7*r/12) + 0.121203875241187*r^3*exp(-7*r/12) + 1.04057961182678*r^2*exp(-7*r/12) + 7.13540305252647*r*exp(-7*r/12) + 36.6963585558504*exp(-7*r/12) + 125.816086477201*exp(-7*r/12)/r - 215.684719675202/r^2 + 215.684719675202*exp(-7*r/12)/r^2,
                 ((1,0), (3,2), 2) => r -> -0.0507150516208487*r^3*exp(-4*r/3) - 0.228217732293819*r^2*exp(-4*r/3) - 0.684653196881458*r*exp(-4*r/3) - 1.54046969298328*exp(-4*r/3) - 2.31070453947492*exp(-4*r/3)/r + 1.73302840460619/r^2 - 1.73302840460619*exp(-4*r/3)/r^2,
                 ((2,0), (3,2), 2) => r -> 0.0229510124219699*r^4*exp(-5*r/6) + 0.174427694406971*r^3*exp(-5*r/6) + 1.05758265240437*r^2*exp(-5*r/6) + 5.07639673154099*r*exp(-5*r/6) + 18.2750282335475*exp(-5*r/6) + 43.8600677605141*exp(-5*r/6)/r - 52.6320813126169/r^2 + 52.6320813126169*exp(-5*r/6)/r^2,
                 ((3,0), (3,2), 2) => r -> -0.0028918862918778*r^5*exp(-2*r/3) - 0.0173513177512668*r^4*exp(-2*r/3) - 0.156161859761401*r^3*exp(-2*r/3) - 1.17121394821051*r^2*exp(-2*r/3) - 7.02728368926307*r*exp(-2*r/3) - 31.6227766016838*exp(-2*r/3) - 94.8683298050514*exp(-2*r/3)/r + 142.302494707577/r^2 - 142.302494707577*exp(-2*r/3)/r^2,
                 ((3,2), (3,2), 2) => r -> -0.000914494741655235*r^5*exp(-2*r/3) - 0.0137174211248285*r^4*exp(-2*r/3) - 0.135802469135802*r^3*exp(-2*r/3) - 1.03703703703704*r^2*exp(-2*r/3) - 6.22222222222222*r*exp(-2*r/3) - 28.0*exp(-2*r/3) - 84.0*exp(-2*r/3)/r + 126.0/r^2 - 126.0*exp(-2*r/3)/r^2,
                 ((4,0), (3,2), 2) => r -> 0.000172500175581118*r^6*exp(-7*r/12) - 0.000591429173420976*r^5*exp(-7*r/12) + 0.00253469645751847*r^4*exp(-7*r/12) + 0.0259746036979988*r^3*exp(-7*r/12) + 0.222473929071335*r^2*exp(-7*r/12) + 1.52553551363201*r*exp(-7*r/12) + 7.84561121296464*exp(-7*r/12) + 26.8992384444502*exp(-7*r/12)/r - 46.1129801904861/r^2 + 46.1129801904861*exp(-7*r/12)/r^2,
                 ((4,2), (3,2), 2) => r -> 7.71444237460059e-5*r^6*exp(-7*r/12) + 0.000661237917822907*r^5*exp(-7*r/12) + 0.00589446372459277*r^4*exp(-7*r/12) + 0.0536331424611299*r^3*exp(-7*r/12) + 0.447720145762475*r^2*exp(-7*r/12) + 3.07008099951412*r*exp(-7*r/12) + 15.7889879975012*exp(-7*r/12) + 54.1336731342897*exp(-7*r/12)/r - 92.8005825159253/r^2 + 92.8005825159253*exp(-7*r/12)/r^2,
                 ((3,2), (4,0), 2) => r -> 0.000172500175581118*r^6*exp(-7*r/12) - 0.000591429173420976*r^5*exp(-7*r/12) + 0.00253469645751847*r^4*exp(-7*r/12) + 0.0259746036979988*r^3*exp(-7*r/12) + 0.222473929071335*r^2*exp(-7*r/12) + 1.52553551363201*r*exp(-7*r/12) + 7.84561121296464*exp(-7*r/12) + 26.8992384444502*exp(-7*r/12)/r - 46.1129801904861/r^2 + 46.1129801904861*exp(-7*r/12)/r^2,
                 ((4,2), (4,0), 2) => r -> -1.51643064880357e-5*r^7*exp(-r/2) + 0.000121314451904285*r^6*exp(-r/2) - 0.00145577342285143*r^5*exp(-r/2) - 0.0116461873828114*r^4*exp(-r/2) - 0.139754248593737*r^3*exp(-r/2) - 1.39754248593737*r^2*exp(-r/2) - 11.1803398874989*r*exp(-r/2) - 67.0820393249937*exp(-r/2) - 268.328157299975*exp(-r/2)/r + 536.656314599949/r^2 - 536.656314599949*exp(-r/2)/r^2,
                 ((2,1), (4,1), 2) => r -> -0.00183002179407892*r^5*exp(-3*r/4) + 0.0122001452938595*r^4*exp(-3*r/4) + 0.0292803487052628*r^3*exp(-3*r/4) + 0.104107906507601*r^2*exp(-3*r/4) + 0.555242168040538*r*exp(-3*r/4) + 2.22096867216215*exp(-3*r/4) + 5.92258312576574*exp(-3*r/4)/r - 7.89677750102099/r^2 + 7.89677750102099*exp(-3*r/4)/r^2,
                 ((3,1), (4,1), 2) => r -> 0.000298779068421049*r^6*exp(-7*r/12) - 0.00162194351428569*r^5*exp(-7*r/12) + 0.0108779971033296*r^4*exp(-7*r/12) + 0.0725013977341712*r^3*exp(-7*r/12) + 0.553920713484602*r^2*exp(-7*r/12) + 3.79831346389441*r*exp(-7*r/12) + 19.5341835285998*exp(-7*r/12) + 66.974343526628*exp(-7*r/12)/r - 114.813160331362/r^2 + 114.813160331362*exp(-7*r/12)/r^2,
                 ((4,1), (4,1), 2) => r -> -2.03450520833333e-5*r^7*exp(-r/2) + 0.000244140625*r^6*exp(-r/2) - 0.00260416666666667*r^5*exp(-r/2) - 0.0104166666666667*r^4*exp(-r/2) - 0.151041666666667*r^3*exp(-r/2) - 1.5625*r^2*exp(-r/2) - 12.5*r*exp(-r/2) - 75.0*exp(-r/2) - 300.0*exp(-r/2)/r + 600.0/r^2 - 600.0*exp(-r/2)/r^2,
                 ((4,3), (4,1), 2) => r -> -4.43965434236638e-6*r^7*exp(-r/2) - 3.55172347389311e-5*r^6*exp(-r/2) - 0.000568275755822897*r^5*exp(-r/2) - 0.00795586058152055*r^4*exp(-r/2) - 0.0954703269782467*r^3*exp(-r/2) - 0.954703269782467*r^2*exp(-r/2) - 7.63762615825973*r*exp(-r/2) - 45.8257569495584*exp(-r/2) - 183.303027798234*exp(-r/2)/r + 366.606055596467/r^2 - 366.606055596467*exp(-r/2)/r^2,
                 ((1,0), (4,2), 2) => r -> 0.00372677996249965*r^4*exp(-5*r/4) - 0.020869967789998*r^3*exp(-5*r/4) - 0.114486680447989*r^2*exp(-5*r/4) - 0.366357377433566*r*exp(-5*r/4) - 0.879257705840557*exp(-5*r/4) - 1.40681232934489*exp(-5*r/4)/r + 1.12544986347591/r^2 - 1.12544986347591*exp(-5*r/4)/r^2,
                 ((2,0), (4,2), 2) => r -> -0.00183002179407892*r^5*exp(-3*r/4) + 0.00122001452938595*r^4*exp(-3*r/4) + 0.0146401743526314*r^3*exp(-3*r/4) + 0.104107906507601*r^2*exp(-3*r/4) + 0.555242168040538*r*exp(-3*r/4) + 2.22096867216215*exp(-3*r/4) + 5.92258312576574*exp(-3*r/4)/r - 7.89677750102099/r^2 + 7.89677750102099*exp(-3*r/4)/r^2,
                 ((3,0), (4,2), 2) => r -> 0.000243952087818557*r^6*exp(-7*r/12) - 0.000104550894779382*r^5*exp(-7*r/12) + 0.0106417875043299*r^4*exp(-7*r/12) + 0.110725798646556*r^3*exp(-7*r/12) + 0.951250981920879*r^2*exp(-7*r/12) + 6.52286387602889*r*exp(-7*r/12) + 33.54615707672*exp(-7*r/12) + 115.015395691611*exp(-7*r/12)/r - 197.169249757048/r^2 + 197.169249757048*exp(-7*r/12)/r^2,
                 ((3,2), (4,2), 2) => r -> 7.71444237460059e-5*r^6*exp(-7*r/12) + 0.000661237917822907*r^5*exp(-7*r/12) + 0.00589446372459277*r^4*exp(-7*r/12) + 0.0536331424611299*r^3*exp(-7*r/12) + 0.447720145762475*r^2*exp(-7*r/12) + 3.07008099951412*r*exp(-7*r/12) + 15.7889879975012*exp(-7*r/12) + 54.1336731342897*exp(-7*r/12)/r - 92.8005825159253/r^2 + 92.8005825159253*exp(-7*r/12)/r^2,
                 ((4,0), (4,2), 2) => r -> -1.51643064880357e-5*r^7*exp(-r/2) + 0.000121314451904285*r^6*exp(-r/2) - 0.00145577342285143*r^5*exp(-r/2) - 0.0116461873828114*r^4*exp(-r/2) - 0.139754248593737*r^3*exp(-r/2) - 1.39754248593737*r^2*exp(-r/2) - 11.1803398874989*r*exp(-r/2) - 67.0820393249937*exp(-r/2) - 268.328157299975*exp(-r/2)/r + 536.656314599949/r^2 - 536.656314599949*exp(-r/2)/r^2,
                 ((4,2), (4,2), 2) => r -> -6.78168402777778e-6*r^7*exp(-r/2) - 2.71267361111111e-5*r^6*exp(-r/2) - 0.000651041666666667*r^5*exp(-r/2) - 0.0104166666666667*r^4*exp(-r/2) - 0.130208333333333*r^3*exp(-r/2) - 1.3125*r^2*exp(-r/2) - 10.5*r*exp(-r/2) - 63.0*exp(-r/2) - 252.0*exp(-r/2)/r + 504.0/r^2 - 504.0*exp(-r/2)/r^2,
                 ((2,1), (4,3), 2) => r -> -0.000399343494989788*r^5*exp(-3*r/4) - 0.00532457993319718*r^4*exp(-3*r/4) - 0.0468563034121352*r^3*exp(-3*r/4) - 0.318054908009645*r^2*exp(-3*r/4) - 1.69629284271811*r*exp(-3*r/4) - 6.78517137087242*exp(-3*r/4) - 18.0937903223265*exp(-3*r/4)/r + 24.1250537631019/r^2 - 24.1250537631019*exp(-3*r/4)/r^2,
                 ((3,1), (4,3), 2) => r -> 6.51989379575165e-5*r^6*exp(-7*r/12) + 0.000950041667380954*r^5*exp(-7*r/12) + 0.0116879075718127*r^4*exp(-7*r/12) + 0.121203875241187*r^3*exp(-7*r/12) + 1.04057961182678*r^2*exp(-7*r/12) + 7.13540305252647*r*exp(-7*r/12) + 36.6963585558504*exp(-7*r/12) + 125.816086477201*exp(-7*r/12)/r - 215.684719675202/r^2 + 215.684719675202*exp(-7*r/12)/r^2,
                 ((4,1), (4,3), 2) => r -> -4.43965434236638e-6*r^7*exp(-r/2) - 3.55172347389311e-5*r^6*exp(-r/2) - 0.000568275755822897*r^5*exp(-r/2) - 0.00795586058152055*r^4*exp(-r/2) - 0.0954703269782467*r^3*exp(-r/2) - 0.954703269782467*r^2*exp(-r/2) - 7.63762615825973*r*exp(-r/2) - 45.8257569495584*exp(-r/2) - 183.303027798234*exp(-r/2)/r + 366.606055596467/r^2 - 366.606055596467*exp(-r/2)/r^2,
                 ((4,3), (4,3), 2) => r -> -9.68812003968254e-7*r^7*exp(-r/2) - 2.71267361111111e-5*r^6*exp(-r/2) - 0.000511532738095238*r^5*exp(-r/2) - 0.0076264880952381*r^4*exp(-r/2) - 0.0933779761904762*r^3*exp(-r/2) - 0.9375*r^2*exp(-r/2) - 7.5*r*exp(-r/2) - 45.0*exp(-r/2) - 180.0*exp(-r/2)/r + 360.0/r^2 - 360.0*exp(-r/2)/r^2,
                 ((4,3), (1,0), 3) => r -> -0.00197202659436654*r^4*exp(-5*r/4) - 0.0126209702039458*r^3*exp(-5*r/4) - 0.0605806569789401*r^2*exp(-5*r/4) - 0.24232262791576*r*exp(-5*r/4) - 0.775432409330433*exp(-5*r/4) - 1.86103778239304*exp(-5*r/4)/r - 2.97766045182886*exp(-5*r/4)/r^2 + 2.38212836146309/r^3 - 2.38212836146309*exp(-5*r/4)/r^3,
                 ((4,3), (2,0), 3) => r -> 0.000968356512192217*r^5*exp(-3*r/4) + 0.0109747071381785*r^4*exp(-3*r/4) + 0.103291361300503*r^3*exp(-3*r/4) + 0.826330890404025*r^2*exp(-3*r/4) + 5.5088726026935*r*exp(-3*r/4) + 29.380653881032*exp(-3*r/4) + 117.522615524128*exp(-3*r/4)/r + 313.393641397675*exp(-3*r/4)/r^2 - 417.858188530233/r^3 + 417.858188530233*exp(-3*r/4)/r^3,
                 ((3,2), (2,1), 3) => r -> -0.0185510824799983*r^4*exp(-5*r/6) - 0.178090391807983*r^3*exp(-5*r/6) - 1.28225082101748*r^2*exp(-5*r/6) - 7.69350492610488*r*exp(-5*r/6) - 36.9288236453034*exp(-5*r/6) - 132.943765123092*exp(-5*r/6)/r - 319.065036295421*exp(-5*r/6)/r^2 + 382.878043554506/r^3 - 382.878043554506*exp(-5*r/6)/r^3,
                 ((4,2), (2,1), 3) => r -> 0.00147918900560809*r^5*exp(-3*r/4) + 0.00197225200747745*r^4*exp(-3*r/4),
                 ((4,3), (3,0), 3) => r -> -0.000129087311236578*r^6*exp(-7*r/12) - 0.00149372460145184*r^5*exp(-7*r/12) - 0.0205209308342312*r^4*exp(-7*r/12) - 0.246312138361854*r^3*exp(-7*r/12) - 2.53349628029335*r^2*exp(-7*r/12) - 21.7156824025145*r*exp(-7*r/12) - 148.907536474385*exp(-7*r/12) - 765.810187582551*exp(-7*r/12)/r - 2625.63492885446*exp(-7*r/12)/r^2 + 4501.08844946479/r^3 - 4501.08844946479*exp(-7*r/12)/r^3,
                 ((3,2), (3,1), 3) => r -> 0.0028628213703701*r^5*exp(-2*r/3) + 0.0257653923333309*r^4*exp(-2*r/3) + 0.257653923333309*r^3*exp(-2*r/3) + 2.31888530999978*r^2*exp(-2*r/3) + 17.3916398249984*r*exp(-2*r/3) + 104.34983894999*exp(-2*r/3) + 469.574275274956*exp(-2*r/3)/r + 1408.72282582487*exp(-2*r/3)/r^2 - 2113.0842387373/r^3 + 2113.0842387373*exp(-2*r/3)/r^3,
                 ((4,2), (3,1), 3) => r -> -0.000241500245813565*r^6*exp(-7*r/12) - 0.000621000632092025*r^5*exp(-7*r/12) - 0.0152588726742612*r^4*exp(-7*r/12) - 0.194664687937418*r^3*exp(-7*r/12) - 2.00226536164202*r^2*exp(-7*r/12) - 17.1622745283602*r*exp(-7*r/12) - 117.68416819447*exp(-7*r/12) - 605.23286500013*exp(-7*r/12)/r - 2075.08410857187*exp(-7*r/12)/r^2 + 3557.28704326607/r^3 - 3557.28704326607*exp(-7*r/12)/r^3,
                 ((2,1), (3,2), 3) => r -> -0.0185510824799983*r^4*exp(-5*r/6) - 0.178090391807983*r^3*exp(-5*r/6) - 1.28225082101748*r^2*exp(-5*r/6) - 7.69350492610488*r*exp(-5*r/6) - 36.9288236453034*exp(-5*r/6) - 132.943765123092*exp(-5*r/6)/r - 319.065036295421*exp(-5*r/6)/r^2 + 382.878043554506/r^3 - 382.878043554506*exp(-5*r/6)/r^3,
                 ((3,1), (3,2), 3) => r -> 0.0028628213703701*r^5*exp(-2*r/3) + 0.0257653923333309*r^4*exp(-2*r/3) + 0.257653923333309*r^3*exp(-2*r/3) + 2.31888530999978*r^2*exp(-2*r/3) + 17.3916398249984*r*exp(-2*r/3) + 104.34983894999*exp(-2*r/3) + 469.574275274956*exp(-2*r/3)/r + 1408.72282582487*exp(-2*r/3)/r^2 - 2113.0842387373/r^3 + 2113.0842387373*exp(-2*r/3)/r^3,
                 ((4,1), (3,2), 3) => r -> -0.000187065286028187*r^6*exp(-7*r/12) - 0.000106894449158964*r^5*exp(-7*r/12) - 0.00690232728855025*r^4*exp(-7*r/12) - 0.0921473782137682*r^3*exp(-7*r/12) - 0.947801604484473*r^2*exp(-7*r/12) - 8.12401375272406*r*exp(-7*r/12) - 55.7075228758221*exp(-7*r/12) - 286.495831932799*exp(-7*r/12)/r - 982.271423769598*exp(-7*r/12)/r^2 + 1683.89386931931/r^3 - 1683.89386931931*exp(-7*r/12)/r^3,
                 ((4,3), (3,2), 3) => r -> -4.08209920534635e-5*r^6*exp(-7*r/12) - 0.000839746122242677*r^5*exp(-7*r/12) - 0.0122363006383933*r^4*exp(-7*r/12) - 0.148069520330137*r^3*exp(-7*r/12) - 1.52300078053856*r^2*exp(-7*r/12) - 13.0542924046162*r*exp(-7*r/12) - 89.5151479173682*exp(-7*r/12) - 460.363617860751*exp(-7*r/12)/r - 1578.38954695114*exp(-7*r/12)/r^2 + 2705.81065191625/r^3 - 2705.81065191625*exp(-7*r/12)/r^3,
                 ((4,3), (4,0), 3) => r -> 8.02419675442114e-6*r^7*exp(-r/2) + 3.20967870176846e-5*r^6*exp(-r/2) + 0.000962903610530536*r^5*exp(-r/2) + 0.0154064577684886*r^4*exp(-r/2) + 0.21569040875884*r^3*exp(-r/2) + 2.58828490510608*r^2*exp(-r/2) + 25.8828490510608*r*exp(-r/2) + 207.062792408487*exp(-r/2) + 1242.37675445092*exp(-r/2)/r + 4969.50701780368*exp(-r/2)/r^2 - 9939.01403560736/r^3 + 9939.01403560736*exp(-r/2)/r^3,
                 ((3,2), (4,1), 3) => r -> -0.000187065286028187*r^6*exp(-7*r/12) - 0.000106894449158964*r^5*exp(-7*r/12) - 0.00690232728855025*r^4*exp(-7*r/12) - 0.0921473782137682*r^3*exp(-7*r/12) - 0.947801604484473*r^2*exp(-7*r/12) - 8.12401375272406*r*exp(-7*r/12) - 55.7075228758221*exp(-7*r/12) - 286.495831932799*exp(-7*r/12)/r - 982.271423769598*exp(-7*r/12)/r^2 + 1683.89386931931/r^3 - 1683.89386931931*exp(-7*r/12)/r^3,
                 ((4,2), (4,1), 3) => r -> 1.64447098157852e-5*r^7*exp(-r/2) - 6.5778839263141e-5*r^6*exp(-r/2) + 0.00171024982084166*r^5*exp(-r/2) + 0.0231541514206256*r^4*exp(-r/2) + 0.315738428463077*r^3*exp(-r/2) + 3.78886114155692*r^2*exp(-r/2) + 37.8886114155692*r*exp(-r/2) + 303.108891324553*exp(-r/2) + 1818.65334794732*exp(-r/2)/r + 7274.61339178928*exp(-r/2)/r^2 - 14549.2267835786/r^3 + 14549.2267835786*exp(-r/2)/r^3,
                 ((2,1), (4,2), 3) => r -> 0.00147918900560809*r^5*exp(-3*r/4) + 0.00197225200747745*r^4*exp(-3*r/4),
                 ((3,1), (4,2), 3) => r -> -0.000241500245813565*r^6*exp(-7*r/12) - 0.000621000632092025*r^5*exp(-7*r/12) - 0.0152588726742612*r^4*exp(-7*r/12) - 0.194664687937418*r^3*exp(-7*r/12) - 2.00226536164202*r^2*exp(-7*r/12) - 17.1622745283602*r*exp(-7*r/12) - 117.68416819447*exp(-7*r/12) - 605.23286500013*exp(-7*r/12)/r - 2075.08410857187*exp(-7*r/12)/r^2 + 3557.28704326607/r^3 - 3557.28704326607*exp(-7*r/12)/r^3,
                 ((4,1), (4,2), 3) => r -> 1.64447098157852e-5*r^7*exp(-r/2) - 6.5778839263141e-5*r^6*exp(-r/2) + 0.00171024982084166*r^5*exp(-r/2) + 0.0231541514206256*r^4*exp(-r/2) + 0.315738428463077*r^3*exp(-r/2) + 3.78886114155692*r^2*exp(-r/2) + 37.8886114155692*r*exp(-r/2) + 303.108891324553*exp(-r/2) + 1818.65334794732*exp(-r/2)/r + 7274.61339178928*exp(-r/2)/r^2 - 14549.2267835786/r^3 + 14549.2267835786*exp(-r/2)/r^3,
                 ((4,3), (4,2), 3) => r -> 3.58852988154377e-6*r^7*exp(-r/2) + 5.74164781047003e-5*r^6*exp(-r/2) + 0.000947371888727555*r^5*exp(-r/2) + 0.0148134513510127*r^4*exp(-r/2) + 0.206699321176921*r^3*exp(-r/2) + 2.48039185412305*r^2*exp(-r/2) + 24.8039185412305*r*exp(-r/2) + 198.431348329844*exp(-r/2) + 1190.58808997907*exp(-r/2)/r + 4762.35235991626*exp(-r/2)/r^2 - 9524.70471983253/r^3 + 9524.70471983253*exp(-r/2)/r^3,
                 ((1,0), (4,3), 3) => r -> -0.00197202659436654*r^4*exp(-5*r/4) - 0.0126209702039458*r^3*exp(-5*r/4) - 0.0605806569789401*r^2*exp(-5*r/4) - 0.24232262791576*r*exp(-5*r/4) - 0.775432409330433*exp(-5*r/4) - 1.86103778239304*exp(-5*r/4)/r - 2.97766045182886*exp(-5*r/4)/r^2 + 2.38212836146309/r^3 - 2.38212836146309*exp(-5*r/4)/r^3,
                 ((2,0), (4,3), 3) => r -> 0.000968356512192217*r^5*exp(-3*r/4) + 0.0109747071381785*r^4*exp(-3*r/4) + 0.103291361300503*r^3*exp(-3*r/4) + 0.826330890404025*r^2*exp(-3*r/4) + 5.5088726026935*r*exp(-3*r/4) + 29.380653881032*exp(-3*r/4) + 117.522615524128*exp(-3*r/4)/r + 313.393641397675*exp(-3*r/4)/r^2 - 417.858188530233/r^3 + 417.858188530233*exp(-3*r/4)/r^3,
                 ((3,0), (4,3), 3) => r -> -0.000129087311236578*r^6*exp(-7*r/12) - 0.00149372460145184*r^5*exp(-7*r/12) - 0.0205209308342312*r^4*exp(-7*r/12) - 0.246312138361854*r^3*exp(-7*r/12) - 2.53349628029335*r^2*exp(-7*r/12) - 21.7156824025145*r*exp(-7*r/12) - 148.907536474385*exp(-7*r/12) - 765.810187582551*exp(-7*r/12)/r - 2625.63492885446*exp(-7*r/12)/r^2 + 4501.08844946479/r^3 - 4501.08844946479*exp(-7*r/12)/r^3,
                 ((3,2), (4,3), 3) => r -> -4.08209920534635e-5*r^6*exp(-7*r/12) - 0.000839746122242677*r^5*exp(-7*r/12) - 0.0122363006383933*r^4*exp(-7*r/12) - 0.148069520330137*r^3*exp(-7*r/12) - 1.52300078053856*r^2*exp(-7*r/12) - 13.0542924046162*r*exp(-7*r/12) - 89.5151479173682*exp(-7*r/12) - 460.363617860751*exp(-7*r/12)/r - 1578.38954695114*exp(-7*r/12)/r^2 + 2705.81065191625/r^3 - 2705.81065191625*exp(-7*r/12)/r^3,
                 ((4,0), (4,3), 3) => r -> 8.02419675442114e-6*r^7*exp(-r/2) + 3.20967870176846e-5*r^6*exp(-r/2) + 0.000962903610530536*r^5*exp(-r/2) + 0.0154064577684886*r^4*exp(-r/2) + 0.21569040875884*r^3*exp(-r/2) + 2.58828490510608*r^2*exp(-r/2) + 25.8828490510608*r*exp(-r/2) + 207.062792408487*exp(-r/2) + 1242.37675445092*exp(-r/2)/r + 4969.50701780368*exp(-r/2)/r^2 - 9939.01403560736/r^3 + 9939.01403560736*exp(-r/2)/r^3,
                 ((4,2), (4,3), 3) => r -> 3.58852988154377e-6*r^7*exp(-r/2) + 5.74164781047003e-5*r^6*exp(-r/2) + 0.000947371888727555*r^5*exp(-r/2) + 0.0148134513510127*r^4*exp(-r/2) + 0.206699321176921*r^3*exp(-r/2) + 2.48039185412305*r^2*exp(-r/2) + 24.8039185412305*r*exp(-r/2) + 198.431348329844*exp(-r/2) + 1190.58808997907*exp(-r/2)/r + 4762.35235991626*exp(-r/2)/r^2 - 9524.70471983253/r^3 + 9524.70471983253*exp(-r/2)/r^3,
                 )


text:
function Base.show(io::IO, obj::DoseTime)
    print(io, "Dose - $(obj.dose); Time - $(obj.time); Tau - $(obj.tau)")

end
function Base.show(io::IO, obj::ElimRange)
    print(io, "Elimination range: $(obj.kelstart) - $(obj.kelend) ")
    if length(obj.kelexcl) > 0
        print(io, "Exclusions: $(obj.kelexcl[1])")
        if length(obj.kelexcl) > 1 for i = 2:length(obj.kelexcl) print(io, ", $(obj.kelexcl[i])") end end
        print(io, ".")
    else
        print(io, "No exclusion.")
    end
end
function Base.show(io::IO, obj::KelData)
    println(io, "Elimination table:")
    header = ["Strat time", "End time", "a", "b", "r²", "Adjusted r²"]
    mt = metida_table(obj.s, obj.e, obj.a, obj.b, obj.r, obj.ar; names =  Tuple(Symbol.(header)))
    PrettyTables.pretty_table(io, mt; tf = PrettyTables.tf_compact, header = header)
end

# PK Subject

function Base.show(io::IO, obj::PKSubject)
    println(io, "  Pharmacokinetic subject")
    println(io, "Observations: $(length(obj)); ", obj.dosetime)
    println(io,  obj.kelrange)
    PrettyTables.pretty_table(io, metida_table(obj.time, obj.obs; names = (:Time, :Concentration)); tf = PrettyTables.tf_compact)

end
function Base.show(io::IO, obj::UPKSubject)
    println(io, "  Pharmacokinetic subject (urine)")
    println(io, "Observations: $(length(obj)); ", obj.dosetime)
    println(io,  obj.kelrange)
    PrettyTables.pretty_table(io, metida_table(getindex.(obj.time, 1), getindex.(obj.time, 2), obj.obs, obj.vol); tf = PrettyTables.tf_compact, header  = ["Start time", "End time", "Concentration", "Volume"])

end
function Base.show(io::IO, obj::DataSet{Subj}) where Subj <: PKSubject
    println(io, "DataSet: Pharmacokinetic subject")
    println(io, "Length: $(length(obj))")
    for i = 1:length(obj)
        println(io, "Subject $(i): ", obj[i].id)
    end
end

function Base.show(io::IO, obj::T) where T <: NCAResult
    println(io, "  Pharmacokinetic subject NCA result")
    PrettyTables.pretty_table(io, obj.result; header = ["Parameter", "Value"], tf = PrettyTables.tf_compact)
end

function Base.show(io::IO, obj::DataSet{Res}) where Res <: NCAResult
    println(io, "DataSet: Pharmacokinetic subject NCA result")
    println(io, "Length: $(length(obj))")
    for i = 1:length(obj)
        println(io, "Subject $(i): ", obj[i].data.id)
    end
end


text:
@snn_kw struct SpikingSynapseParameter{FT=Float32}
    τpre::FT = 20ms
    τpost::FT = 20ms
    Wmax::FT = 0.01
    ΔApre::FT = 0.01 * Wmax
    ΔApost::FT = -ΔApre * τpre / τpost * 1.05
end

@snn_kw mutable struct SpikingSynapse{VIT=Vector{Int32},VFT=Vector{Float32},VBT=Vector{Bool}}
    param::SpikingSynapseParameter = SpikingSynapseParameter()
    rowptr::VIT # row pointer of sparse W
    colptr::VIT # column pointer of sparse W
    I::VIT      # postsynaptic index of W
    J::VIT      # presynaptic index of W
    index::VIT  # index mapping: W[index[i]] = Wt[i], Wt = sparse(dense(W)')
    W::VFT  # synaptic weight
    tpre::VFT = zero(W) # presynaptic spiking time
    tpost::VFT = zero(W) # postsynaptic spiking time
    Apre::VFT = zero(W) # presynaptic trace
    Apost::VFT = zero(W) # postsynaptic trace
    fireI::VBT # postsynaptic firing
    fireJ::VBT # presynaptic firing
    g::VFT # postsynaptic conductance
    records::Dict = Dict()
end

"""
[Spking Synapse](https://brian2.readthedocs.io/en/2.0b4/resources/tutorials/2-intro-to-brian-synapses.html)
"""
SpikingSynapse

function SpikingSynapse(pre, post, sym; σ = 0.0, p = 0.0, kwargs...)
    w = σ * sprand(post.N, pre.N, p)
    rowptr, colptr, I, J, index, W = dsparse(w)
    fireI, fireJ = post.fire, pre.fire
    g = getfield(post, sym)
    SpikingSynapse(;@symdict(rowptr, colptr, I, J, index, W, fireI, fireJ, g)..., kwargs...)
end

function forward!(c::SpikingSynapse, param::SpikingSynapseParameter)
    @unpack colptr, I, W, fireJ, g = c
    @inbounds for j in 1:(length(colptr) - 1)
        if fireJ[j]
            for s in colptr[j]:(colptr[j+1] - 1)
                g[I[s]] += W[s]
            end
        end
    end
end

function plasticity!(c::SpikingSynapse, param::SpikingSynapseParameter, dt::Float32, t::Float32)
    @unpack rowptr, colptr, I, J, index, W, tpre, tpost, Apre, Apost, fireI, fireJ, g = c
    @unpack τpre, τpost, Wmax, ΔApre, ΔApost = param
    @inbounds for j in 1:(length(colptr) - 1)
        if fireJ[j]
            for s in colptr[j]:(colptr[j+1] - 1)
                Apre[s] *= exp32(- (t - tpre[s]) / τpre)
                Apost[s] *= exp32(- (t - tpost[s]) / τpost)
                Apre[s] += ΔApre
                tpre[s] = t
                W[s] = clamp(W[s] + Apost[s], 0f0, Wmax)
            end
        end
    end
    @inbounds for i in 1:(length(rowptr) - 1)
        if fireI[i]
            for st in rowptr[i]:(rowptr[i+1] - 1)
                s = index[st]
                Apre[s] *= exp32(- (t - tpre[s]) / τpre)
                Apost[s] *= exp32(- (t - tpost[s]) / τpost)
                Apost[s] += ΔApost
                tpost[s] = t
                W[s] = clamp(W[s] + Apre[s], 0f0, Wmax)
            end
        end
    end
end


text:
@testset "Functions to extend" begin

    @testset "predicates" begin
        @test GroupsCore.isabelian isa Function
        @test GroupsCore.issolvable isa Function
        @test GroupsCore.isnilpotent isa Function
        @test GroupsCore.isperfect isa Function
    end

    @testset "generic constructions" begin
        @test GroupsCore.derivedsubgroup isa Function
        @test GroupsCore.center isa Function
        @test GroupsCore.socle isa Function
        @test GroupsCore.sylowsubgroup isa Function
    end

    @testset "general subgroups" begin
        @test GroupsCore.centralizer isa Function
        @test GroupsCore.normalizer isa Function
        @test GroupsCore.stabilizer isa Function
    end
end


text:
using PyCall

let dependencies = ["catboost==0.24.8"]
    pip = pyimport("pip")
    pip.main(["install"; split(get(ENV, "PIPFLAGS", "")); dependencies])
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule rocPRIM_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("rocPRIM")
JLLWrappers.@generate_main_file("rocPRIM", UUID("52935e6f-76c5-5ebb-b227-36676f75be9c"))
end  # module rocPRIM_jll


text:
"""
	countapp(v,u)

Counts how many times each value in `u` appears in `v`. The Function returns a `Vector{Float64}` containing the element count.

# Examples
```
julia> countapp([1,2,3,2,1,2],[2,3,1])
3-element Array{Float64,1}:
 3.0
 1.0
 2.0

julia> countapp([1,2,3,2,1,2],[1,2,4,5,1])
5-element Array{Float64,1}:
 2.0
 3.0
 0.0
 0.0
 2.0
```
"""
function countapp(v::T where T<:AbstractVector{V}, u::S where S<:AbstractVector{V}) where V
	countapp!(zeros(Float64, length(u)), v, u)
end

function countapp!(out::Vector{Float64}, v::T where T<:AbstractVector{V}, u::S where S<:AbstractVector{V}) where V
	m = length(u)
	@assert length(out) == m "[countapp!] output vector size must be $m."

	@inbounds for i in 1:m				# for after 'u' first: reduce allocations significantly
		out[i] = 0.0
		@simd for vi in v
			out[i] +=(vi==u[i])
		end
	end
	return out
end

"""
	countappw(v,u [,w,val])

Counts how many times each value in `u` appears in `v`. Each value in `v` can have an associated weight in the vector `w`. The counts
corresponding to a value in `u` not appearing in `v` can receive a cusom value `val` (default `0.0`).
The function returns a `Vector{Float64}` containing the weighted normalized element count.

# Examples
```
julia> j4pr.countappw([1,2,3,2,1,2],[1,2,4,5,1])
5-element Array{Float64,1}:
 0.333333
 0.5
 0.0
 0.0
 0.333333

julia> j4pr.countappw([1,2,3,2,1,2],[1,2,4,5,1],ones(6),-1.0)
5-element Array{Float64,1}:
 2.0
 3.0
-1.0
-1.0
 2.0
```
"""
function countappw(v::T where T<:AbstractVector{V}, u::S where S<:AbstractVector{V},
			w::Vector{Float64}=fill(1/length(v),length(v)), val::Float64=0.0) where V
	countappw!(zeros(Float64, length(u)), v, u, w, val)
end

function countappw!(out::Vector{Float64}, v::T where T<:AbstractVector{V}, u::S where S<:AbstractVector{V},
	   		w::Vector{Float64}=fill(1/length(v),length(v)), val::Float64=0.0) where V
	m = length(u)
	n = length(v)
	@assert length(out) == m "[countappw!] output vector size must be $m."

	mask = falses(m)			# vector that keeps trak of the modified positions in `out`; non-modified elements get val at the end
	@inbounds for i in 1:m			# for after 'u' first: reduce allocations significantly
		out[i] = 0.0
		@simd for j in 1:n
			b=ifelse(v[j]==u[i],true,false)
			mask[i] |= b
			out[i] += b*w[j]
		end
	end
	out[.!mask] = val
	return out
end

"""
	gini(p)

Compute the gini impurity of an array `p`.
"""
function gini(p::AbstractArray{T}) where T<:Real
	s = zero(T)
	@inbounds @simd for i in 1:length(p)
		s += p[i]^2
	end
	return 1-s
end

"""
	misclassification(p)

Compute  the misclassification impurity of an array `p`.
"""
function misclassification(p::AbstractArray{T}) where T<:Real
	m = zero(T)
	@inbounds for i in 1:length(p)
		if p[i] > m m=p[i] end
	end
	return 1-m
end

"""
	linearsplit(v, n [;count=0, prop=0])

Generate a sorted vector of `n` linearly spaced values starting from the values of an
input vector `v`. The input vector may be trimmed using the parameters `prop` and `count`
as in `StatsBase.trim`.
"""
function linearsplit(v::T where T<:AbstractVector, n::Int; prop::Float64=0.0, count::Int=0)
	vt = trim(v, prop=prop, count=count)		# trim
	n = min(n, length(unique(vt)))			# number of points
	L = linspace(minimum(vt),maximum(vt),n+1); 	# generate space
	return collect( (L.+L.step.hi/2.0)[1:end-1] )	# take midpoints
end

"""
	densitysplitv, n [;count=0, prop=0])

Generate a sorted vector of `n` values, spaced according to the value density of an
input vector `v`. The input vector may be trimmed using the parameters `prop` and `count`
as in `StatsBase.trim`.
"""
function densitysplit(v::T where T<:AbstractVector, n::Int; prop::Float64=0.0, count::Int=0)
	if length(v)/length(unique(v)) < n
		v .+= 10*rand(length(v))*eps()
	end
	vt = trim(v, prop=prop, count=count)		# trim
	n = min(n, length(unique(vt)))			# number of points
	vt = sort(sample(vt, n, replace=false))		# sort, sample
	vo = Vector{Float64}(n)				# preallocate output
	@inbounds @simd for i in 1:n-1
		vo[i] = vt[i] + (vt[i+1]-vt[i])/2	# threshold position
	end
	vo[n] = vt[n]-eps()				# the last point
	return vo
end

"""
	confusionmatrix(predictions, references [;kwargs])

Creates and returns the confusion matrix corresponding to the `predictions`
and `references` input label vectors. The two vectors must have the same element
type and be of the same size. The returned confusion matrix is a `Matrix{Float64}`.

# Keyword arguments
  * `showmatrix::Bool` specifies whether the matrix should be printed (default `false`)
  * `normalize::Bool` specifies whether the matrix should normalized with respect to the number of labels (default `false`)
  * `positive` is the value of the positive or target class; if its value is `nothing`, the
  confusion matrix is calculating considering all classes while is a value of the same type
  as the element type of `predictions` and `references` is specified (the value must be present
  in the `references` vector, it will consider that class having the label `true` and
  all other classes as having the label `false`

# Examples
```
julia> using j4pr

julia> references = ["a","a","c","c","b","b"];

julia> predictions = ["a","b","c","c","a","a"];

julia> confusionmatrix(predictions,references)
3×3 Array{Float64,2}:
 1.0  2.0  0.0
 1.0  0.0  0.0
 0.0  0.0  2.0

julia> confusionmatrix(predictions,references;normalize=true)
3×3 Array{Float64,2}:
 0.5  1.0  0.0
 0.5  0.0  0.0
 0.0  0.0  1.0

julia> confusionmatrix(predictions,references;normalize=true,positive="a")
2×2 Array{Float64,2}:
 0.5  0.5
 0.5  0.5

julia> confusionmatrix(predictions,references;normalize=true,positive="a",showmatrix=true);

reference labels (columns), "a" is "true":
 "true"  "false"
------------
0.5   0.5
0.5   0.5
------------
```
"""
function confusionmatrix(predictions::AbstractArray{T}, references::AbstractArray{T};
			showmatrix::Bool=false, normalize::Bool=false, positive=nothing) where T
	@assert length(predictions) == length(references) "[confusionmatrix] The predicted and reference labels should have the same length."

	# If positive class is specified, binarize labels
	_binarize_(predictions,references,::Void) =
		predictions, references, sort(unique(predictions)), sort(unique(references))

	_binarize_(predictions::AbstractArray{T}, references::AbstractArray{T}, positive::T) = begin
		@assert positive in references "[confusionmatrix] $(positive) is not in the reference label vector."
		yb = falses(length(predictions))
		yrb = falses(length(references))
		yb[predictions.==positive] = true
		yrb[references.==positive] = true
		return yb, yrb, sort(unique(yb),rev=true), sort(unique(yrb),rev=true)
	end
	y, yr, yu, yru = _binarize_(predictions, references, positive)

	# Construct confusion matrix
	C = length(yru)
	@assert issubset(yu,yru) "[confusionmatrix] The predicted labels should be a subset of the reference labels."
	cm::Matrix{Float64} = zeros(C,C)
	@inbounds for (j,cr) in enumerate(yru)
		for (i,cp) in enumerate(yru)
			cm[i,j] = sum((yr .== cr) .& (y .== cp))
		end
	end

	# Check for normalization
	if normalize
		# Loop through the classes and normalize the columns
		# of the confusion matrix with respect to their sum
		# i.e. the sum of each column should be 1.0
		for j in 1:C
			cm[:,j]/=sum(yr .==yru[j])
		end
	end

	# Check if the matrix should be nicely printed or not
	if showmatrix
		println()
		if !(positive isa Void)
			println("reference labels (columns), \"$(positive)\" is \"true\":")
		else
			println("reference labels (columns):")
		end
		lsize=ceil(Int, log10(length(y)))+2
		println(sprint((io::IO,v)->map(x->print(io,lpad(" \"$x\" ",lsize)),v), yru))

		println(repeat("-", (lsize+3)*C))
		for i in 1:size(cm,1)
			for j in 1:size(cm,2)
				print(lpad("$(cm[i,j])   ",lsize))
			end
			println()
		end
		println(repeat("-", (lsize+3)*C))
	end
	return cm
end

# Function that generates integers based on  priors, that sum up to a given number N;
# N random numbers are generated and the cumulative distribution is thresholded so that
# a fraction of the total sum (given by the priors) falls closely to a cardinal in the
# cumulative distribution
#
# Input:
# 	N 	- sum of the numbers
# 	priors 	- prior probabilities (must sum up to one)
# Output:
#	out 	- vector of integers that sum up to N and has the same size as the priors vector
function pintgen(n::Int, priors::Vector{T} where T<:Real)
	if (checkpriors(priors) != true)
		error("[pintgen] Priors must summate to one and be >=0, <=1.")
	end
	m = length(priors)									# Length of the priors vector (and also of output vector)
	out = zeros(m) 										# Initialize output
	v = rand(n) 										# Generate random variable vector
	csv = cumsum(v)
	sv = csv[n]
	csp = cumsum(priors)
	for i = 1 : m
		pos = find(csv .<= csp[i] * sv);						# Find positions where the cumulative sum is smaller
		if isempty(pos)  								# than the total sum multiplied by the prior
			out[i] = 0
		else
			out[i] = maximum(pos)							# Return the last position or 0, depending wether any
		end 										# position has been found or not
		if (i > 1) out[i] = out[i] - sum(out[1:i-1]) end				# Compensate positions for using the cumulative priors
	end
	return round.(Int,out)
end

# Version of the function that ignores priors and returns whatever the first argument is;
# Basically a wrapper so that both specifying class numbers and total sample size works transparently
function pintgen(n::Vector{Int}, priors::Vector{T} where T<:Real)
	@assert checkpriors(priors) && (size(n) == size(priors))
		"[pintgen] Priors are not correct or size mismatch between class sizes and priors."
	return n
end


text:
# test errors
# prepare_dataframe_for_dtm

@testset verbose = true "errorTesting" begin

    thisfile = joinpath(datadir, "GermanMotorPremiums", "data1Small.csv")
    @test isfile(thisfile)
    @time df_tmp = CSV.read(thisfile, DataFrame,strict=true, types=eltypesData1,  copycols=true);

    selected_explanatory_vars = ["PLZ_WOHNORT","ART_DES_WOHNEIGENTUM","GEBURTSDATUM","FAMILIENSTAND","NATIONALITAET","GESCHLECHT","FINANZIERUNGSART","STADT","KENNZEICHEN"]

# keep 10 largest PLZ only
    vorig = deepcopy(df_tmp[!,:PLZ_WOHNORT])
    counts, freqs, vals, keep = getCounts(vorig, threshold=10)
    vnew = df_tmp[!,:PLZ_WOHNORT]
    mapToOther!(vnew, keep, 9999999)

##################################################
# Prepare the data
##################################################

# examples of test_throws
    @test_throws ErrorException error("some text")
    @test_throws ErrorException("some text") error("some text") # this is more exact and also tests the text thrown
    @test_throws AssertionError @assert false

    dtmtable, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars);

# @test_throws AssertionError prepare_dataframe_for_dtm!(df_tmp,treat_as_categorical_variable=["PLZ_WOHNORT"],weightcol="EXPOSURE",numocl="LOSS20HALF",trnvalcol="does_not_exist",denomcol="PREMIUM66",independent_vars=selected_explanatory_vars);
    @test_throws MethodError prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numocl="LOSS20HALF", trnvalcol="does_not_exist", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars);

    @test_throws AssertionError prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars);
    @test_throws AssertionError prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="does_not_exist", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars);

    @test_throws AssertionError prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="does_not_exist", independent_vars=selected_explanatory_vars);

    @test_throws AssertionError prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="does_not_exist", numcol="LOSS20HALF", independent_vars=selected_explanatory_vars);

    @test_throws AssertionError prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], keycol="does_not_exist", numcol="LOSS20HALF", independent_vars=selected_explanatory_vars);

    @test_throws ErrorException prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["does_not_exist","PLZ_WOHNORT"], numcol="LOSS20HALF", independent_vars=selected_explanatory_vars);

# @test_throws Base.UVError prepare_dataframe_for_dtm!(df_tmp,directory="K:\\does_not_existIUDFLKJSDFEAA\\aaEFIEF_EFDXxT\\",treat_as_categorical_variable=["PLZ_WOHNORT"],numcol="LOSS20HALF",independent_vars=selected_explanatory_vars);
# @test_throws Base.IOError prepare_dataframe_for_dtm!(df_tmp,directory="K:\\does_not_existIUDFLKJSDFEAA\\aaEFIEF_EFDXxT\\",treat_as_categorical_variable=["PLZ_WOHNORT"],numcol="LOSS20HALF",independent_vars=selected_explanatory_vars);
# @test_throws AssertionError prepare_dataframe_for_dtm!(df_tmp,directory="K:\\does_not_existIUDFLKJSDFEAA\\aaEFIEF_EFDXxT\\",treat_as_categorical_variable=["PLZ_WOHNORT"],numcol="LOSS20HALF",independent_vars=selected_explanatory_vars);

    @test_throws AssertionError prepare_dataframe_for_dtm!(df_tmp, valpct=-22.1, treat_as_categorical_variable=["PLZ_WOHNORT"], numcol="LOSS20HALF", independent_vars=selected_explanatory_vars);

end # testset


text:
with_repo_access("test isconnected") do test_repo, path
    skip, dlfile = false, tempname()
    try
       download("https://github.com", dlfile)
    catch
        skip = true
    finally
        if isfile(dlfile)
            rm(dlfile)
        end
    end
    if !skip
        remote = GitRemote(test_repo, "git://github.com/libgit2/libgit2.git")
        connect(remote, :fetch) do r
            @test connected(r)
        end
        @test disconnected(remote)
    end
end

with_repo_access("test remote names") do test_repo, path
    ns = remote_names(test_repo)
    @test sort(ns) == sort(["test_remote", "libgit2"])
end

with_repo_access("test remotes") do test_repo, path
    rs = remotes(test_repo)
    @test sort(map(r -> name(r), rs)) == sort(["test_remote", "libgit2"])
end

with_repo_access("test remote new name") do test_repo, path
    remote = GitRemote(test_repo, "git://github.com/libgit2/libgit2.git")
    @test name(remote) == nothing
    @test "git://github.com/libgit2/libgit2.git" == url(remote)
end

with_repo_access("test remote new invalid url") do test_repo, path
    #TODO: throw custom LibGitError ?
    @test_throws ArgumentError GitRemote(test_repo, "libgit2")
end

with_repo_access("test remote url set") do test_repo, path
    new_url = "git://github.com/libgit2/TestGitRepository.git"
    remote = GitRemote(test_repo, "git://github.com/libgit2/libgit2.git")
    set_url!(remote, new_url)
    @test url(remote) == new_url
end

with_repo_access("test set remote invalid url") do test_repo, path
    new_url = "invalid"
    remote = GitRemote(test_repo, "git://github.com/libgit2/libgit2.git")
    @test_throws ArgumentError set_url!(remote, new_url)
end

with_repo_access("test remote push url") do test_repo, path
    @test push_url(lookup_remote(test_repo, "test_remote")) == "git://github.com/libgit2/TestEmptyRepository.git"
    @test push_url(lookup_remote(test_repo, "libgit2")) == nothing
end

with_repo_access("test remote push url set") do test_repo, path
    new_url = "git://github.com/libgit2/TestGitRepository.git"
    remote = GitRemote(test_repo, "git://github.com/libgit2/libgit2.git")

    @test push_url(remote) == nothing
    set_push_url!(remote, new_url)
    @test push_url(remote) == new_url
end

with_repo_access("test remote push url set invalid") do test_repo, path
    new_url = "upstream"
    remote = GitRemote(test_repo, "git://github.com/libgit2/libgit2.git")
    @test_throws ArgumentError set_push_url!(remote, new_url)
end

with_repo_access("test remote fetch refspecs") do test_repo, path
    remote = lookup_remote(test_repo, "test_remote")
    @test ["+refs/heads/*:refs/remotes/test_remote/*"] == fetch_refspecs(remote)
    @test isempty(push_refspecs(lookup_remote(test_repo, "libgit2")))
end

with_repo_access("test remote push refspecs") do test_repo, path
    remote = lookup_remote(test_repo, "test_remote")
    @test ["refs/heads/*:refs/heads/testing/*"] == push_refspecs(remote)
    @test isempty(push_refspecs(lookup_remote(test_repo, "libgit2")))
end

with_repo_access("test remote add fetch") do test_repo, path
    remote = GitRemote(test_repo, "git://github.com/libgit2/libgit2.git")
    add_fetch!(remote, "+refs/heads/*:refs/remotes/test/*")
    @test ["+refs/heads/*:refs/remotes/test/*"] == fetch_refspecs(remote)
end

with_repo_access("test remote add push") do test_repo, path
    remote = GitRemote(test_repo, "git://github.com/libgit2/libgit2.git")
    add_push!(remote, "refs/heads/*:refs/heads/test/*")
    @test ["refs/heads/*:refs/heads/test/*"] == push_refspecs(remote)
end

with_repo_access("test remote clear refspecs") do test_repo, path
    remote = lookup_remote(test_repo, "test_remote")
    clear_refspecs!(remote)

    @test isempty(push_refspecs(remote))
    @test isempty(fetch_refspecs(remote))
end

with_repo_access("test remote lookup") do test_repo, path
    remote = lookup_remote(test_repo, "libgit2")
    @test url(remote)  == "git://github.com/libgit2/libgit2.git"
    @test name(remote) == "libgit2"
end

with_repo_access("test remote lookup missing") do test_repo, path
    @test lookup_remote(test_repo, "missing_remote") == nothing
end

with_repo_access("test remote lookup invalid") do test_repo, path
    @test_throws LibGitError{:Config,:InvalidSpec} lookup_remote(test_repo, "*\?")
end

#TODO: make this work on windows
@unix_only begin
    with_tmp_repo_access("test remote add") do test_repo, path
        remote_add!(test_repo, "upstream", "git://github.com/libgit2/libgit2.git")
        remote = lookup_remote(test_repo, "upstream")
        @test name(remote) == "upstream"
        @test url(remote)  == "git://github.com/libgit2/libgit2.git"
    end
end

with_tmp_repo_access("test remote add with invalid url") do test_repo, path
    @test_throws ArgumentError remote_add!(test_repo, "upstream", "libgit2")
end

with_tmp_repo_access("test remote url set") do test_repo, path
    new_url = "git://github.com/l?#!@#ibgit2/TestGitRepository.git"
    remote = lookup_remote(test_repo, "origin")
    set_url!(remote, new_url)
    save!(remote)
    @test url(lookup_remote(test_repo, "origin")) == new_url
end

with_tmp_repo_access("test remote rename") do test_repo, path
    remote = lookup_remote(test_repo, "origin")
    rename!(remote, "new_remote_name")
    @test lookup_remote(test_repo, "new_remote_name") != nothing
end

with_tmp_repo_access("test rename invalid name") do test_repo, path
    remote = lookup_remote(test_repo, "origin")
    @test_throws LibGitError{:Config,:InvalidSpec} rename!(remote, "/?")
end

with_tmp_repo_access("test remote rename exists") do test_repo, path
    remote = lookup_remote(test_repo, "origin")
    @test_throws LibGitError{:Config,:Exists} rename!(remote, "origin")
end

with_tmp_repo_access("test remote rename error callback") do test_repo, path
    GitConfig(test_repo)["remote.origin.fetch"]  = "+refs/*:refs/*"
    remote = lookup_remote(test_repo, "origin")
    @test ["+refs/*:refs/*"] == rename!(remote, "test_remote")
end

sandboxed_clone_test("testrepo.git", "test remote push single ref") do test_repo, test_clone, path
    create_ref(test_repo, "refs/heads/unit_test", Oid("8496071c1b46c854b31185ea97743be6a8774479"))
    remote = lookup_remote(test_repo, "origin")
    result = push!(test_clone, remote,
                   ["refs/heads/master",
                    "refs/heads/master:refs/heads/foobar",
                    "refs/heads/unit_test"])
    @test isempty(result)
    @test target(lookup_ref(test_clone, "refs/heads/foobar")) == Oid("a65fedf39aefe402d3bb6e24df4d4f5fe4547750")
    @test target(lookup_ref(test_clone, "refs/heads/unit_test")) == Oid("8496071c1b46c854b31185ea97743be6a8774479")
end

sandboxed_clone_test("testrepo.git", "test push non bare raise error") do test_repo, test_clone, path
    GitConfig(test_clone)["core.bare"] = false
    create_ref(test_repo, "refs/heads/unit_test",
               Oid("8496071c1b46c854b31185ea97743be6a8774479"))
    remote = lookup_remote(test_repo, "origin")
    #TODO: better errors
    @test_throws ErrorException push!(test_clone, remote, ["refs/heads/master"])
end

sandboxed_clone_test("testrepo.git", "test remote push non ff raise error") do test_repo, test_clone, path
    create_ref(test_repo, "refs/heads/unit_test",
               Oid("8496071c1b46c854b31185ea97743be6a8774479"))
    remote = lookup_remote(test_repo, "origin")
    #TODO: better errors
    @test_throws ErrorException push!(test_clone, remote,
                                      ["refs/heads/unit_test:refs/heads/master"])
    @test target(lookup_ref(test_clone, "refs/heads/master")) == Oid("a65fedf39aefe402d3bb6e24df4d4f5fe4547750")
end

sandboxed_clone_test("testrepo.git", "test remote push non ff raise no error") do test_repo, test_clone, path
    create_ref(test_repo, "refs/heads/unit_test",
               Oid("8496071c1b46c854b31185ea97743be6a8774479"))
    remote = lookup_remote(test_repo, "origin")
    result = push!(test_clone, remote, ["+refs/heads/unit_test:refs/heads/master"])
    @test isempty(result)
    @test target(lookup_ref(test_clone, "refs/heads/master")) == Oid("8496071c1b46c854b31185ea97743be6a8774479")
end

remote_transport_test("test remote disconnect") do test_repo, test_remote
    connect(test_remote, :fetch)
    @test connected(test_remote)
    disconnect(test_remote)
    @test disconnected(test_remote)
end

remote_transport_test("test remote ls") do test_repo, test_remote
    connect(test_remote, :fetch) do r
        rheads = LibGit2.ls(r)
        @test length(rheads) == 7
        rhead = rheads[1]
        @test rhead.islocal  == false
        @test rhead.id  != nothing
        @test rhead.lid  == nothing
    end
end

remote_transport_test("test remote git fetch") do test_repo, test_remote
    connect(test_remote, :fetch) do r
        download(r) |> update_tips!
    end
    @test Oid(revparse(test_repo, "origin/master")) == Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    @test test_repo["36060c5"] != nothing
end