Merge pull request #27 from numericalEFT/modulize-chebyshev

Modulize chebyshev

README.md

@@ -0,0 +1,160 @@
+[![img](https://img.shields.io/badge/docs-dev-blue.svg)](https://numericaleft.github.io/CompositeGrids.jl/dev/)
+[![img](https://github.com/numericaleft/CompositeGrids.jl/workflows/CI/badge.svg)](https://github.com/numericaleft/CompositeGrids.jl/actions)
+[![img](https://codecov.io/gh/numericalEFT/CompositeGrids.jl/branch/main/graph/badge.svg?token=WN6HO1XASY)](https://codecov.io/gh/numericaleft/CompositeGrids.jl)
+
+
+# Introduction
+
+CompositeGrids gives a unified interface to generate various common 1D grids
+and also the composite grids that is a combination of basic grids,
+together with the floor function, interpolation function and also integration function
+that is optimized for some of the grids.
+
+
+# Quick Start
+
+In the following example we show how to generate a τ grid from 0 to β, log-densed at 0 and β,
+and optimized for integration. The description is attached in the comments in the code.
+
+```julia
+    using CompositeGrids
+    β = 10
+
+    # Generating a log densed composite grid with LogDensedGrid()
+    tgrid = CompositeGrid.LogDensedGrid(
+        :gauss,# The top layer grid is :gauss, optimized for integration. For interpolation use :cheb
+        [0.0, β],# The grid is defined on [0.0, β]
+        [0.0, β],# and is densed at 0.0 and β, as given by 2nd and 3rd parameter.
+        5,# N of log grid
+        0.005, # niminum interval length of log grid
+        5 # N of bottom layer
+    )
+    # The grid has 3 layers.
+    # The top layer is defined by the boundary and densed points. In this case its:
+    println("Top layer:",tgrid.panel.grid)
+    # The middle layer is a log grid with 4 points and minimum interval length 0.001:
+    println("First subgrid of middle layer:",tgrid.subgrids[1].panel.grid)
+    # The bottom layer is a Gauss-Legendre grid with 5 points:
+    println("First subgrid of bottom layer:",tgrid.subgrids[1].subgrids[1].grid)
+
+    # function to be integrated:
+    f(t) = exp(t)+exp(β-t)
+    # numerical value on grid points:
+    data = [f(t) for (ti, t) in enumerate(tgrid.grid)]
+
+    # integrate with integrate1D():
+    int_result = Interp.integrate1D(data, tgrid)
+
+    println("result=",int_result)
+    println("comparing to:",2*(exp(β)-1))
+```
+
+```
+    Top layer:[0.0, 5.0, 10.0]
+    First subgrid of middle layer:[0.0, 0.005000000000000001, 0.05000000000000001, 0.5, 5.0]
+    First subgrid of bottom layer:[0.00023455038515334025, 0.0011538267247357924, 0.0025000000000000005, 0.0038461732752642086, 0.004765449614846661]
+    result=44050.91248775534
+    comparing to:44050.931589613436
+```
+
+# Installation
+
+Static version could be installed via standard package manager with Pkg.add("CompositeGrids").
+
+For developing version, git clone this repo and add with Pkg.develop("directory/of/the/repo").
+
+
+# Manual
+
+
+## Basics
+
+The grids are provided in two modules, SimpleGrid and CompositeGrid. SimpleGrid consists of several
+common 1D grids that is defined straightforward and has simple structure. CompositeGrid defines a
+general type of grids composed by a panel grid and a set of subgrids. The common interface of grids
+are the following:
+
+-   g.bound gives the boundary of the interval of the grid.
+-   g.size gives the total number of grid points.
+-   g.grid gives the array of grid points.
+-   g[i] returns the i-th grid point, same as g.grid[i].
+-   floor(g, x) returns the largest index of grid point where g[i]<x. Return 1 for x<g[1] and (grid.size-1) for x>g[end], so that both floor() and (floor()+1) are valid grid indices.
+
+Interpolation and integration are also provided, with different implemented functions for different grids.
+
+
+## Simple Grids
+
+Various basic grids are designed for use and also as components of composite grids, including:
+Arbitrary, Uniform, Log, BaryCheb, and GaussLegendre.
+
+Arbitrary grid is the most general basic grid, which takes an array and turn it into a grid.
+An O(ln(N)) floor function based on searchsortedfirst() is provided.
+
+Uniform grid is defined by the boundary and number of grid points.
+An O(1) floor function is provided.
+
+Log grid is defined by the boundary, number of grid points, minimum interval, and also the direction.
+A log densed grid is generated according to the parameters provided.
+For example:
+
+```julia
+    using CompositeGrids
+    loggrid = SimpleGrid.Log{Float64}([0.0,1.0], 6, 0.0001, true)
+    println(loggrid.grid)
+```
+```
+    [0.0, 0.00010000000000000005, 0.0010000000000000002, 0.010000000000000002, 0.1, 1.0]
+```
+
+An O(1) floor function is provided.
+
+BaryCheb grid is designed for interpolation. It's defined by the boundary and number of grid points,
+but the grid points are not distributed uniformly. The floor function is not optimized
+so the O(ln(N)) function will be used, but the interpolation is based on an optimized algorithm.
+
+GaussLegendre grid is designed for integration. It's defined by the boundary and number of grid points,
+but the grid points are not distributed uniformly. The floor function is not optimized
+so the O(ln(N)) function will be used. The 1D integration is optimized.
+
+Also notice that there's open grids and closed grids. Closed grids means that the boundary points are
+also grid points, while open grids means the opposite. Only BaryCheb and GaussLegendre are open.
+
+A detailed manual can be found [here](https://numericaleft.github.io/CompositeGrids.jl/dev/lib/simple/).
+
+
+## Composite Grids
+
+Composite grid is a general type of grids where the whole interval is first divided by a panel grid,
+then each interval of a panel grid is divided by a smaller grid in subgrids. Subgrid could also be
+composite grid.
+
+LogDensedGrid is a useful generator of CompositeGrid which gives a general solution when an 1D grid on an
+interval is needed to be log-densed around several points. For example, &tau; grids need to be densed around
+0 and &beta;, and momentum grids need to be densed around Fermi momentum.
+The grid is defined as a three-layer composite grid with the top layer being an Arbitrary grid defined by
+the boundary and densed points, the middle layer a Log grid which is densed at the points required, and the
+bottom layer a grid of three options. Three types are :cheb, :gauss, and :uniform, which corresponds to
+BaryCheb grid for interpolation, GaussLegendre grid for integration, and Uniform grid for general use.
+The floor function is defined recursively, i.e. the floor function of the panel grid is called to find the
+corresponding subgrid, and then the floor function of the subgrid is called to find the result. Since the
+subgrids could also be CompositeGrid, this process continues until the lowest level of the subgrids is reached.
+
+A detailed manual can be found [here](https://numericaleft.github.io/CompositeGrids.jl/dev/lib/composite/).
+
+
+## Interpolation and Integration
+
+Interpolation gives an estimate of the function value at x with given grid and function value on the grid.
+For most of the simple grids the interpolation is given by linear interpolation with the floor function to find
+the corresponding grid points. BaryCheb uses an optimized algorithm for interpolation which makes use of the information
+of all grid points, and thus gives a more precise interpolation with the same number of grid points, given the condition that
+the function itself is smooth enough. For composite grids, the interpolation is done recursively, so that the final result
+depends on the type of lowest level grid. Interpolation for higher dimension where the data is defined on a list of grids is also
+given, but only linear interpolation is implemented, even when some of the grids are BaryCheb.
+
+Integration over 1D grid is also provided. For most of simple grids it's given by linear integral, while for GaussLegendre grid it's
+optimized. For composite grids it's again recursively done so that the method depends on the type of lowest level grids.
+
+A detailed manual can be found [here](https://numericaleft.github.io/CompositeGrids.jl/dev/lib/interpolate/).
+

src/CompositeGrids.jl

@@ -4,6 +4,9 @@ using StaticArrays
 include("old/grid.jl")
 export Grid
 
+include("grid/chebyshev.jl")
+export BaryChebTools
+
 include("grid/simple.jl")
 const SimpleGrid = SimpleG # alias for older convention
 const AbstractGrid, OpenGrid, ClosedGrid = SimpleGrid.AbstractGrid, SimpleGrid.OpenGrid, SimpleGrid.ClosedGrid

src/grid/chebyshev.jl

@@ -1,12 +1,18 @@
+module BaryChebTools
+
+using ..StaticArrays
+
+export BaryCheb1D, interp1D, interpND, integrate1D, integrateND
+######################################################
+#---------------- 1D barycheb ------------------------
+######################################################
+
 """
 barychebinit(n)
-
 Get Chebyshev nodes of first kind and corresponding barycentric Lagrange interpolation weights. 
 Reference: Berrut, J.P. and Trefethen, L.N., 2004. Barycentric lagrange interpolation. SIAM review, 46(3), pp.501-517.
-
 # Arguments
 - `n`: order of the Chebyshev interpolation
-
 # Returns
 - Chebyshev nodes
 - Barycentric Lagrange interpolation weights
@@ -23,6 +29,7 @@ function barychebinit(n)
 end
 
 function vandermonde(x)
+    # generate vandermonde matrix for x
     n = length(x)
     vmat = zeros(Float64, (n, n))
     for i in 1:n
@@ -34,6 +41,7 @@ function vandermonde(x)
 end
 
 function invvandermonde(order)
+    # generate inverse of vandermonde matrix for cheb x of given order
     n = order
     vmat = zeros(Float64, (n, n))
     for i in 1:n
@@ -58,11 +66,13 @@ end
     end
 end
 
+
 function weightcoef(a, i::Int, n)
     # integrate when i=1; differentiate when i=-1
     b = zeros(Float64, n)
     for j in 1:n
         if j+i-1 > 0
+            # b[j] = a^(j+i-1)/(j+i-1)
             b[j] = a^(j+i-1) * divfactorial(j, i)
         elseif j+i-1 == 0
             b[j] = 1
@@ -79,17 +89,14 @@ end
 
 """
 function barycheb(n, x, f, wc, xc)
-
 Barycentric Lagrange interpolation at Chebyshev nodes
 Reference: Berrut, J.P. and Trefethen, L.N., 2004. Barycentric lagrange interpolation. SIAM review, 46(3), pp.501-517.
-
 # Arguments
 - `n`: order of the Chebyshev interpolation
 - `x`: coordinate to interpolate
 - `f`: array of size n, function at the Chebyshev nodes
 - `wc`: array of size n, Barycentric Lagrange interpolation weights
 - `xc`: array of size n, coordinates of Chebyshev nodes
-
 # Returns
 - Interpolation result
 """
@@ -109,6 +116,47 @@ function barycheb(n, x, f, wc, xc)
     return num / den
 end
 
+function barychebND(n, xs, f, wc, xc, DIM)
+    haseq = false
+    eqinds = zeros(Int, DIM)
+    for i in 1:DIM
+        for j in 1:n
+            if xs[i] == xc[j]
+                eqinds[i] = j
+                haseq = true
+            end
+        end
+    end
+
+    if haseq
+        newxs = [xs[i] for i in 1:DIM if eqinds[i]==0]
+        newDIM = length(newxs)
+        if newDIM == 0
+            return f[CartesianIndex(eqinds...)]
+        else
+            newf = view(f, [(i==0) ? (1:n) : (i) for i in eqinds]...)
+            return _barychebND_noneq(n, newxs, newf, wc, xc, newDIM)
+        end
+    else
+        return _barychebND_noneq(n, xs, f, wc, xc, DIM)
+    end
+end
+
+function _barychebND_noneq(n, xs, f, wc, xc, DIM)
+    # deal with the case when there's no xs[i] = xc[j]
+    inds = CartesianIndices(NTuple{DIM, Int}(ones(Int, DIM) .* n))
+    num, den = 0.0, 0.0
+    for (indi, ind) in enumerate(inds)
+        q = 1.0
+        for i in 1:DIM
+            q *= wc[ind[i]] / (xs[i] - xc[ind[i]])
+        end
+        num += q * f[indi]
+        den += q
+    end
+    return num / den
+end
+
 function barycheb2(n, x, f, wc, xc)
     for j in 1:n
         if x == xc[j]
@@ -138,3 +186,75 @@ function chebdiff(n, x, f, invmat)
     return sum(intw .* f)
 end
 
+struct BaryCheb1D{N}
+    # wrapped barycheb 1d grid, x in [-1, 1]
+    x::SVector{N, Float64}
+    w::SVector{N, Float64}
+    invmat::SMatrix{N, N, Float64}
+
+    function BaryCheb1D(N::Int)
+        x, w = barychebinit(N)
+        invmat = invvandermonde(N)
+
+        return new{N}(x, w, invmat)
+    end
+end
+
+Base.getindex(bc::BaryCheb1D, i) = bc.x[i]
+
+function interp1D(data, xgrid::BaryCheb1D{N}, x) where {N}
+    return barycheb(N, x, data, xgrid.w, xgrid.x)
+end
+
+function integrate1D(data, xgrid::BaryCheb1D{N}; x1=-1, x2=1) where {N}
+    return chebint(N, x1, x2, data, xgrid.invmat)
+end
+
+function interpND(data, xgrid::BaryCheb1D{N}, xs) where {N}
+    return barychebND(N, xs, data, xgrid.w, xgrid.x, length(xs))
+end
+
+function integrateND(data, xgrid::BaryCheb1D{N}, x1s, x2s) where {N}
+    DIM = length(x1s)
+    @assert DIM == length(x2s)
+
+    intws = zeros(Float64, (DIM, N))
+    for i in 1:DIM
+        wc = weightcoef(x2s[i], 1, N) .- weightcoef(x1s[i], 1, N)
+        intws[i, :] = calcweight(xgrid.invmat, wc)
+    end
+
+    result = 0.0
+    inds = CartesianIndices(NTuple{DIM, Int}(ones(Int, DIM) .* N))
+    for (indi, ind) in enumerate(inds)
+        w = 1.0
+        for i in 1:DIM
+            w *= intws[i, ind[i]]
+        end
+        result += data[indi] * w
+    end
+
+    return result
+end
+
+function integrateND(data, xgrid::BaryCheb1D{N}, DIM) where {N}
+    @assert N ^ DIM == length(data)
+
+    intws = zeros(Float64, (DIM, N))
+    wc = weightcoef(1.0, 1, N) .- weightcoef(-1.0, 1, N)
+    intw = calcweight(xgrid.invmat, wc)
+
+    result = 0.0
+    inds = CartesianIndices(NTuple{DIM, Int}(ones(Int, DIM) .* N))
+    for (indi, ind) in enumerate(inds)
+        w = 1.0
+        for i in 1:DIM
+            w *= intw[ind[i]]
+        end
+        result += data[indi] * w
+    end
+
+    return result
+end
+
+end