diff --git a/README.md b/README.md index 545b5ab..ef77c27 100644 --- a/README.md +++ b/README.md @@ -22,12 +22,12 @@ and optimized for integration. The description is attached in the comments in th # 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, # minimum interval length of log grid - 5 # N of bottom layer + type=:gauss,# The top layer grid is :gauss, optimized for integration. For interpolation use :cheb + bound=[0.0, β],# The grid is defined on [0.0, β] + dense_at=[0.0, β],# and is densed at 0.0 and β, as given by 2nd and 3rd parameter. + N=5,# N of log grid + minterval=0.005, # minimum interval length of log grid + order=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: @@ -110,11 +110,11 @@ For example: 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 +but the grid points are distributed according to Chebyshev nodes. 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 +but the grid points are distributed according to Gauss Legendre quadrature. 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 @@ -156,5 +156,38 @@ given, but only linear interpolation is implemented, even when some of the grids 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. +Differenciation is provided for 1D grid, and a high precision algorithm is implemented for BaryCheb grids. + A detailed manual can be found [here](https://numericaleft.github.io/CompositeGrids.jl/dev/lib/interpolate/). +An example of interpolation and differenciation is shown below: +```julia +using CompositeGrids +β = π + +# Generating a log densed composite grid with LogDensedGrid() +tgrid = CompositeGrid.LogDensedGrid( + type=:cheb,# The top layer grid is :cheb + bound=[0.0, β],# The grid is defined on [0.0, β] + dense_at=[0.0, β],# and is densed at 0.0 and β, as given by 2nd and 3rd parameter. + N=5,# N of log grid + minterval=0.005, # minimum interval length of log grid + order=5 # N of bottom layer +) + +# function to be represented: +f(t) = sin(t) +# numerical value on grid points: +data = [f(t) for t in tgrid] + +# integrate with integrate1D(): +sin1 = Interp.interp1D(data, tgrid, 1.0) +dsin1 = Interp.differentiate1D(data, tgrid, 1.0) + +println("result=", (sin1, dsin1)) +println("comparing to:", (sin(1.0), cos(1.0))) +``` +``` +result=(0.8414425112056995, 0.5400742649805592) +comparing to:(0.8414709848078965, 0.5403023058681398) +``` \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index 5489321..0874c3c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -11,20 +11,19 @@ makedocs(; format=Documenter.HTML(; prettyurls=get(ENV, "CI", "false") == "true", canonical="https://numericaleft.github.io/CompositeGrids.jl", - assets=String[], + assets=String[] ), pages=[ - "Home" => "index.md", - "Manual" => Any[ - ], + "Home" => "README.md", + "API reference" => "index.md", "Library" => Any[ "lib/simple.md", "lib/composite.md", "lib/interpolate.md", - # map(s -> "lib/$(s)", sort(readdir(joinpath(@__DIR__, "src/lib")))) - # "Internals" => map(s -> "lib/$(s)", sort(readdir(joinpath(@__DIR__, "src/lib")))) + # map(s -> "lib/$(s)", sort(readdir(joinpath(@__DIR__, "src/lib")))) + # "Internals" => map(s -> "lib/$(s)", sort(readdir(joinpath(@__DIR__, "src/lib")))) ] - ], + ] ) deploydocs(; diff --git a/docs/src/README.md b/docs/src/README.md new file mode 100644 index 0000000..ef77c27 --- /dev/null +++ b/docs/src/README.md @@ -0,0 +1,193 @@ +[![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( + type=:gauss,# The top layer grid is :gauss, optimized for integration. For interpolation use :cheb + bound=[0.0, β],# The grid is defined on [0.0, β] + dense_at=[0.0, β],# and is densed at 0.0 and β, as given by 2nd and 3rd parameter. + N=5,# N of log grid + minterval=0.005, # minimum interval length of log grid + order=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]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 distributed according to Chebyshev nodes. 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 distributed according to Gauss Legendre quadrature. 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, τ grids need to be densed around +0 and β, 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. + +Differenciation is provided for 1D grid, and a high precision algorithm is implemented for BaryCheb grids. + +A detailed manual can be found [here](https://numericaleft.github.io/CompositeGrids.jl/dev/lib/interpolate/). + +An example of interpolation and differenciation is shown below: +```julia +using CompositeGrids +β = π + +# Generating a log densed composite grid with LogDensedGrid() +tgrid = CompositeGrid.LogDensedGrid( + type=:cheb,# The top layer grid is :cheb + bound=[0.0, β],# The grid is defined on [0.0, β] + dense_at=[0.0, β],# and is densed at 0.0 and β, as given by 2nd and 3rd parameter. + N=5,# N of log grid + minterval=0.005, # minimum interval length of log grid + order=5 # N of bottom layer +) + +# function to be represented: +f(t) = sin(t) +# numerical value on grid points: +data = [f(t) for t in tgrid] + +# integrate with integrate1D(): +sin1 = Interp.interp1D(data, tgrid, 1.0) +dsin1 = Interp.differentiate1D(data, tgrid, 1.0) + +println("result=", (sin1, dsin1)) +println("comparing to:", (sin(1.0), cos(1.0))) +``` +``` +result=(0.8414425112056995, 0.5400742649805592) +comparing to:(0.8414709848078965, 0.5403023058681398) +``` \ No newline at end of file diff --git a/src/CompositeGrids.jl b/src/CompositeGrids.jl index 5d0c876..14efab1 100644 --- a/src/CompositeGrids.jl +++ b/src/CompositeGrids.jl @@ -1,8 +1,8 @@ module CompositeGrids using StaticArrays -include("old/grid.jl") -export Grid +#include("old/grid.jl") +#export Grid include("grid/chebyshev.jl") export BaryChebTools diff --git a/src/grid/chebyshev.jl b/src/grid/chebyshev.jl index 354afc8..b3cdb66 100644 --- a/src/grid/chebyshev.jl +++ b/src/grid/chebyshev.jl @@ -27,6 +27,16 @@ function barychebinit(n) end return x, w end +function barychebinit(::Type{T}, n) where {T} + x = zeros(T, n) + w = similar(x) + for i in 1:n + c = (2i - 1)π / (2n) + x[n - i + 1] = cos(c) + w[n - i + 1] = (-1)^(i - 1) * sin(c) + end + return x, w +end function vandermonde(x) # generate vandermonde matrix for x @@ -40,6 +50,18 @@ function vandermonde(x) return vmat end +function vandermonde(::Type{T}, x) where {T} + # generate vandermonde matrix for x + n = length(x) + vmat = zeros(T, (n, n)) + for i in 1:n + for j in 1:n + vmat[i,j] = x[i]^(j-1) + end + end + return vmat +end + function invvandermonde(order) # generate inverse of vandermonde matrix for cheb x of given order n = order @@ -54,6 +76,21 @@ function invvandermonde(order) return inv(transpose(vmat)) end +function invvandermonde(::Type{T}, order) where {T} + # generate inverse of vandermonde matrix for cheb x of given order + n = order + vmat = zeros(T, (n, n)) + for i in 1:n + for j in 1:n + c = (2n-2i+1)π/(2n) + x = cos(c) + vmat[i,j] = x^(j-1) + end + end + return inv(transpose(vmat)) +end + + @inline function divfactorial(j, i) if i == 0 return 1 @@ -83,6 +120,25 @@ function weightcoef(a, i::Int, n) return b end + +function weightcoef(::Type{T}, a, i::Int, n) where {T} + # integrate when i=1; differentiate when i=-1 + b = zeros(T, 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 + else + b[j] = 0 + end + end + return b +end + + + @inline function calcweight(invmat, b) return invmat*b end @@ -101,6 +157,8 @@ Reference: Berrut, J.P. and Trefethen, L.N., 2004. Barycentric lagrange interpol - Interpolation result """ function barycheb(n, x, f, wc, xc) + #print("type of data in cheb $(typeof.([n,x,f,wc,xc]))\n") + #print("interp: $(f[1])\n") for j in 1:n if x == xc[j] return f[j] @@ -186,39 +244,48 @@ function chebdiff(n, x, f, invmat) return sum(intw .* f) end -struct BaryCheb1D{N} +struct BaryCheb1D{N,T} # wrapped barycheb 1d grid, x in [-1, 1] - x::SVector{N, Float64} - w::SVector{N, Float64} - invmat::SMatrix{N, N, Float64} + x::SVector{N, T} + w::SVector{N, T} + invmat::SMatrix{N, N, T} +end - function BaryCheb1D(N::Int) - x, w = barychebinit(N) - invmat = invvandermonde(N) +function BaryCheb1D(N::Int,T) + x, w = barychebinit(T,N) + invmat = invvandermonde(T,N) - return new{N}(x, w, invmat) - end + return BaryCheb1D{N,T}(x, w, invmat) end +function BaryCheb1D(N::Int) + x, w = barychebinit(Float64,N) + invmat = invvandermonde(Float64,N) + + return BaryCheb1D{N,Float64}(x, w, invmat) +end + +BaryCheb1D{N}(x, w, invmat) where{N} = BaryCheb1D{N,Float64}(x, w, invmat) + Base.getindex(bc::BaryCheb1D, i) = bc.x[i] -function interp1D(data, xgrid::BaryCheb1D{N}, x) where {N} +function interp1D(data, xgrid::BaryCheb1D{N,T}, x) where {N,T} return barycheb(N, x, data, xgrid.w, xgrid.x) end -function integrate1D(data, xgrid::BaryCheb1D{N}; x1=-1, x2=1) where {N} +function integrate1D(data, xgrid::BaryCheb1D{N,T}; x1=-1, x2=1) where {N,T} return chebint(N, x1, x2, data, xgrid.invmat) end -function interpND(data, xgrid::BaryCheb1D{N}, xs) where {N} +function interpND(data, xgrid::BaryCheb1D{N,T}, xs) where {N,T} return barychebND(N, xs, data, xgrid.w, xgrid.x, length(xs)) end -function integrateND(data, xgrid::BaryCheb1D{N}, x1s, x2s) where {N} +function integrateND(data, xgrid::BaryCheb1D{N,T}, x1s, x2s) where {N,T} DIM = length(x1s) @assert DIM == length(x2s) - intws = zeros(Float64, (DIM, N)) + intws = zeros(T, (DIM, N)) for i in 1:DIM wc = weightcoef(x2s[i], 1, N) .- weightcoef(x1s[i], 1, N) intws[i, :] = calcweight(xgrid.invmat, wc) @@ -237,10 +304,10 @@ function integrateND(data, xgrid::BaryCheb1D{N}, x1s, x2s) where {N} return result end -function integrateND(data, xgrid::BaryCheb1D{N}, DIM) where {N} +function integrateND(data, xgrid::BaryCheb1D{N,T}, DIM) where {N,T} @assert N ^ DIM == length(data) - intws = zeros(Float64, (DIM, N)) + intws = zeros(T, (DIM, N)) wc = weightcoef(1.0, 1, N) .- weightcoef(-1.0, 1, N) intw = calcweight(xgrid.invmat, wc) diff --git a/src/grid/composite.jl b/src/grid/composite.jl index 33840fe..25b3f8a 100644 --- a/src/grid/composite.jl +++ b/src/grid/composite.jl @@ -211,6 +211,8 @@ function CompositeLogGrid(type, bound, N, minterval, d2s, order, T=Float64, invV end +CompositeLogGrid(; type, bound, N, minterval, d2s, order, T=Float64, invVandermonde=SimpleG.invvandermonde(order)) = CompositeLogGrid(type, bound, N, minterval, d2s, order, T, invVandermonde) + """ function LogDensedGrid(type, bound, dense_at, N, minterval, order, T=Float64) @@ -242,17 +244,21 @@ function LogDensedGrid(type, bound, dense_at, N, minterval, order, T=Float64; is @assert bound[1] <= dense_at[1] <= dense_at[end] <= bound[2] dp = Vector{T}([]) + # construct panel points + # combine two panel points when they are closer than 2minterval + # this is to ensure 1. subgrids are constructed on a interval that is not too small + # and 2. the minimal interval is taken around both points for i in 1:length(dense_at) if i == 1 - if abs(dense_at[i] - bound[1]) < minterval + if abs(dense_at[i] - bound[1]) <= 2minterval push!(dp, bound[1]) - elseif abs(dense_at[i] - bound[2]) < minterval + elseif abs(dense_at[i] - bound[2]) <= 2minterval push!(dp, bound[2]) else push!(dp, dense_at[i]) end elseif i != length(dense_at) - if abs(dense_at[i] - dp[end]) < minterval + if abs(dense_at[i] - dp[end]) <= 2minterval if dp[end] != bound[1] dp[end] = (dense_at[i] + dense_at[i-1]) / 2.0 end @@ -260,13 +266,13 @@ function LogDensedGrid(type, bound, dense_at, N, minterval, order, T=Float64; is push!(dp, dense_at[i]) end else - if abs(dense_at[i] - bound[2]) < minterval - if abs(dp[end] - bound[2]) < minterval + if abs(dense_at[i] - bound[2]) <= 2minterval + if abs(dp[end] - bound[2]) <= 2minterval dp[end] = bound[2] else push!(dp, bound[2]) end - elseif abs(dense_at[i] - dp[end]) < minterval + elseif abs(dense_at[i] - dp[end]) <= 2minterval if dp[end] != bound[1] dp[end] = (dense_at[i] + dense_at[i-1]) / 2.0 end @@ -326,4 +332,6 @@ function LogDensedGrid(type, bound, dense_at, N, minterval, order, T=Float64; is end +LogDensedGrid(; type, bound, dense_at, N, minterval, order, T=Float64, isperiodic=false) = LogDensedGrid(type, bound, dense_at, N, minterval, order, T; isperiodic=isperiodic) + end diff --git a/src/grid/interpolate.jl b/src/grid/interpolate.jl index 5ae060e..2595d37 100644 --- a/src/grid/interpolate.jl +++ b/src/grid/interpolate.jl @@ -452,7 +452,7 @@ function interp1D(data::AbstractMatrix, xgrid::T, x; axis=1, method=InterpStyle( end end -function interp1D(data::AbstractVector, xgrid::T, x; axis=1, method=InterpStyle(T)) where {DT,T} +function interp1D(data::AbstractVector, xgrid::T, x; axis=1, method=InterpStyle(T)) where {T} return interp1D(method, data, xgrid, x) end diff --git a/src/grid/simple.jl b/src/grid/simple.jl index f3857be..c46a19b 100644 --- a/src/grid/simple.jl +++ b/src/grid/simple.jl @@ -79,6 +79,11 @@ struct Arbitrary{T<:Real,BT} <: AbstractGrid{T} # allow customized bound that's different from [grid[1], grid[end]] @assert bound[1] <= grid[1] @assert bound[2] >= grid[end] + @assert allunique(grid) "Grid points should be unique!" + if !(issorted(grid)) + @warn "Input grid is not sorted, using sorted instead." + grid = sort(grid) + end BT = typeof(BoundType(BTIN, bound, [grid[1], grid[end]])) if BT == PeriodicBound && (bound[1] == grid[1] && bound[2] == grid[end]) size = length(grid) - 1 @@ -168,7 +173,7 @@ Base.length(grid::AbstractGrid) = grid.size Base.size(grid::AbstractGrid) = (grid.size,) Base.size(grid::AbstractGrid, I::Int) = grid.size -Base.view(grid::AbstractGrid, inds...) where {N} = Base.view(grid.grid, inds...) +Base.view(grid::AbstractGrid, inds...) = Base.view(grid.grid, inds...) # set is not allowed for grids Base.getindex(grid::AbstractGrid, i) = grid.grid[i] Base.firstindex(grid::AbstractGrid) = 1 @@ -228,8 +233,10 @@ create Uniform grid. """ function Uniform{T,BTIN}(bound, N; gpbound=bound) where {T<:AbstractFloat,BTIN} + @assert bound[1] < bound[2] @assert bound[1] <= gpbound[1] @assert bound[2] >= gpbound[2] + @assert N > 0 BT = typeof(BoundType(BTIN, bound, gpbound)) if BT == PeriodicBound && (bound[1] == gpbound[1] && bound[2] == gpbound[2]) Ntot = N @@ -324,6 +331,8 @@ struct BaryCheb{T<:AbstractFloat} <: AbstractGrid{T} create BaryCheb grid. """ function BaryCheb{T}(bound, N) where {T<:AbstractFloat} + @assert bound[1] < bound[2] + @assert N > 1 order = N x, w = barychebinit(order) grid = zeros(T, N) @@ -372,6 +381,8 @@ struct GaussLegendre{T<:AbstractFloat} <: AbstractGrid{T} create GaussLegendre grid. """ function GaussLegendre{T}(bound, N) where {T<:AbstractFloat} + @assert bound[1] < bound[2] + @assert N > 1 order = N x, w = gausslegendre(order) grid = zeros(T, N) @@ -419,6 +430,9 @@ struct Log{T<:AbstractFloat} <: AbstractGrid{T} create Log grid. """ function Log{T}(bound, N, minterval, d2s) where {T<:AbstractFloat} + @assert bound[1] < bound[2] + @assert N > 1 + @assert 0 < minterval < bound[2] - bound[1] grid = zeros(T, N) M = N - 2 λ = (minterval / (bound[2] - bound[1]))^(1.0 / M) @@ -449,6 +463,8 @@ struct Log{T<:AbstractFloat} <: AbstractGrid{T} end end +Log{T}(; bound, N, minterval, d2s) where {T<:AbstractFloat} = Log{T}(bound, N, minterval, d2s) + function Base.show(io::IO, grid::Log; isSimplified=false) if isSimplified print(io, diff --git a/src/old/grid.jl b/src/old/grid.jl deleted file mode 100644 index 7705dbd..0000000 --- a/src/old/grid.jl +++ /dev/null @@ -1,530 +0,0 @@ -module Grid - -# export Log, UniformGrid, tauGrid, fermiKGrid, boseKGrid - -using StaticArrays - -@enum GridType LOG UNIFORM UNILOG - -struct Coeff{T <: AbstractFloat} - bound::SVector{2,T} - idx::SVector{2,T} - λ::T - a::T - b::T - - function Coeff{T}(bound, idx, λ, dense2sparse::Bool) where {T} - λ /= abs(idx[2] - idx[1]) - if dense2sparse == false - bound = (bound[2], bound[1]) - idx = (idx[2], idx[1]) - λ = -λ - end - _l1, _l2 = T(1.0), exp(λ * (idx[2] - idx[1])) - b = (bound[2] - bound[1]) / (_l2 - _l1) - a = (bound[1] * _l2 - bound[2] * _l1) / (_l2 - _l1) - return new{T}(bound, idx, λ, a, b) - end -end - -struct UniLog{T <: AbstractFloat} - bound::SVector{2,T} - idx::SVector{2,T} - isopen::SVector{2,Bool} - M::Int - N::Int - λ::T - d2s::Bool - - function UniLog{T}(bound, init, minterval::T,M::Int,N::Int, dense2sparse::Bool,isopen = @SVector[false,true]) where {T} - @assert N*minterval=1 - pos = l.d2s ? l.idx[2] : l.idx[1] - return Base.floor(Int,pos) - end - - i_m= Base.floor(log(norm)/log(l.λ)) - if i_m>=l.M - if l.d2s - result = Base.floor(head+norm/l.λ^(l.M)*l.N) - return (l.isopen[1]&&result==head) ? (result+1) : result - else - result = Base.floor(tail-norm/l.λ^(l.M)*l.N) - return (l.isopen[2]&&result==tail) ? (result-1) : result - end - end - i_n= Base.floor((norm-l.λ^(i_m+1))/(l.λ^(i_m)-l.λ^(i_m+1))*l.N) - if l.d2s - result = Base.floor(head+(l.M-i_m)*l.N+i_n) - return (l.isopen[2]&&result==tail) ? result-1 : result - else - result = Base.floor((tail-(l.M-i_m)*l.N-i_n)-1) - return (l.isopen[1]&&result==head) ? result+1 : result - end -end - -function _grid(l::Coeff{T}, pos) where {T} - return l.a + l.b * exp(l.λ * (T(pos) - l.idx[1])) -end - -function _floor(l::Coeff{T}, x::T) where {T} - pos = l.idx[1] + log((x - l.a) / l.b) / l.λ - return Base.floor(Int, pos) -end - -function checkOrder(grid) - for idx = 2:length(grid) - @assert grid[idx - 1] < grid[idx] "The grid at $idx is not in the increase order: \n$grid" - end -end - -struct Log{T,SIZE,SEG} # create a log grid of the type T with SIZE grids and SEG of segments - grid::MVector{SIZE,T} - size::Int - head::T - tail::T - segment::SVector{SEG,Int} # ends of each segments - coeff::SVector{SEG,Coeff{T}} - isopen::SVector{2,Bool} - - function Log{T,SIZE,SEG}(coeff, range, isopen) where {T <: AbstractFloat,SIZE,SEG} - @assert SIZE > 1 "Size must be large than 1" - for ri = 2:SEG - @assert range[ri - 1][end] + 1 == range[ri][1] "ranges must be connected to each other" - end - @assert range[1][1] == 1 "ranges should start with the idx 1" - @assert range[end][end] == SIZE "ranges ends with $(range[end][end]), expected $SIZE" - - grid, segment = [], [] - for s = 1:SEG - push!(segment, range[s][end]) - for idx in range[s] - push!(grid, _grid(coeff[s], idx)) - end - end - head, tail = grid[1], grid[end] - isopen[1] && (grid[1] += eps(T) * 1e4) - isopen[2] && (grid[end] -= eps(T) * 1e4) - checkOrder(grid) - return new{T,SIZE,SEG}(grid, SIZE, head, tail, segment, coeff, isopen) - end -end - -function Base.floor(grid::Log{T,SIZE,2}, x) where {T,SIZE} - (grid.head <= x <= grid.tail) || error("$x is out of the uniform grid range!") - segment = grid.segment - if x < grid[2] - return 1 - elseif x < grid[segment[1]] - return _floor(grid.coeff[1], x) - elseif x < grid[segment[1] + 1] - return segment[1] - elseif x < grid[end - 1] - return _floor(grid.coeff[2], x) - else - return SIZE - 1 - end -end - -function Base.floor(grid::Log{T,SIZE,3}, x) where {T,SIZE} - (grid.head <= x <= grid.tail) || error("$x is out of the uniform grid range!") - segment = grid.segment - if x < grid[2] - return 1 - elseif x < grid[segment[1]] - return _floor(grid.coeff[1], x) - elseif x < grid[segment[1] + 1] - return segment[1] - elseif x < grid[segment[2]] - return _floor(grid.coeff[2], x) - elseif x < grid[segment[2] + 1] - return segment[2] - elseif x < grid[end - 1] - return _floor(grid.coeff[3], x) - else - return SIZE - 1 - end -end - -Base.getindex(grid::Log, i) = grid.grid[i] -Base.firstindex(grid::Log) = 1 -Base.lastindex(grid::Log) = grid.size - -""" - Uniform{Type,SIZE} - -Create a uniform Grid with a given type and size - -# Member: -- `β`: inverse temperature -- halfLife: the grid is densest in the range (0, halfLife) and (β-halfLife, β) -- size: the Grid size -- grid: vector stores the grid -- size: size of the grid vector -- head: grid head -- tail: grid tail -- δ: distance between two grid elements -- isopen: if isopen[1]==true, then grid[1] will be slightly larger than the grid head. Same for the tail. -""" -struct Uniform{T,SIZE} - grid::SVector{SIZE,T} - size::Int - head::T - tail::T - δ::T - isopen::SVector{2,Bool} - - """ - Uniform{Type,SIZE}(head, tail, isopen) - - Create a uniform Grid with a given type and size - - # Arguments: - - head: the starting point of the grid - - tail: the end of the grid - - isopen: if isopen[1]==true, then grid[1]=head+eps; If isopen[2]==true, then grid[2]=tail-eps. Otherwise, grid[1]==head / grid[2]==tail - """ - function Uniform{T,SIZE}(head, tail, isopen) where {T <: AbstractFloat,SIZE} - @assert SIZE > 1 "Size must be large than 1" - grid = Array(LinRange(T(head), T(tail), SIZE)) - isopen[1] && (grid[1] += eps(T)) - isopen[2] && (grid[end] -= eps(T)) - return new{T,SIZE}(grid, SIZE, head, tail, (tail - head) / (SIZE - 1), isopen) -end -end - -function Base.floor(grid::Uniform, x) - - (grid.head <= x <= grid.tail) || error("$x is out of the uniform grid range!") - - if grid[2] <= x < grid[end - 1] - return floor(Int, (x - grid.head) / grid.δ + 1) - elseif x < grid[2] - return 1 - else - return grid.size - 1 -end -end - -Base.getindex(grid::Uniform, i) = grid.grid[i] -Base.firstindex(grid::Uniform) = 1 -Base.lastindex(grid::Uniform) = grid.size - - -struct UniLogs{T<:AbstractFloat,SIZE,SEG} - grid::MVector{SIZE,T} - size::Int - head::T - tail::T - unilogs::SVector{SEG,UniLog{T}} - segment::SVector{SEG,T} # ends of each segments - segindex::SVector{SEG,Int} - isopen::SVector{2,Bool} - - function UniLogs{T,SIZE,SEG}(bounds, minterval::T,M::Int,N::Int, Isopen = @SVector[true,true], issparse = @SVector[false,false]) where {T<:AbstractFloat,SIZE,SEG} - @assert SEG > 0 - size = (M+1)*N*SEG + 1 - @assert SIZE == size - - grid, segment,segindex = [], [],[] - unilogs = [] - if issparse[1]==false - for s = 1:SEG - if s%2==1 - bound = @SVector[ bounds[(s+1)÷2], (bounds[(s+1)÷2+1]+bounds[(s+1)÷2])/2] - if s == SEG && issparse[2] - bound = @SVector[ bounds[(s+1)÷2],bounds[(s+1)÷2+1]] - end - init = 1 + (M+1)*N*(s-1) - push!(segment,bound[2]) - push!(segindex, 1 + (M+1)*N*(s)) - isopen = @SVector[false, true] - if s == SEG && issparse[2] - isopen = @SVector[false, false] - end - g = UniLog{T}(bound,init,minterval,M,N,true,isopen) - push!(unilogs,g) - for idx=g.idx[1]:g.idx[2] - push!(grid, _grid(g, idx)) - end - else - bound = @SVector[ (bounds[s÷2+1]+bounds[s÷2])/2,bounds[s÷2+1]] - init = 1 + (M+1)*N*(s-1) - push!(segment, bound[2]) - push!(segindex, 1 + (M+1)*N*(s)) - isopen = @SVector[false,s==SEG ? false : true] - g = UniLog{T}(bound,init,minterval,M,N,false,isopen) - push!(unilogs,g) - for idx=g.idx[1]:g.idx[2] - push!(grid, _grid(g, idx)) - end - end - end - else - for s = 1:SEG - if s%2==1 - bound = @SVector[(bounds[(s+1)÷2+1]+bounds[(s+1)÷2])/2,bounds[(s+1)÷2+1]] - if s == 1 - bound = @SVector[ bounds[1], bounds[2]] - end - init = 1 + (M+1)*N*(s-1) - push!(segment,bound[2]) - push!(segindex, 1 + (M+1)*N*(s)) - isopen = @SVector[false, s==SEG ? false : true] - g = UniLog{T}(bound,init,minterval,M,N,false,isopen) - push!(unilogs,g) - for idx=g.idx[1]:g.idx[2] - push!(grid, _grid(g, idx)) - end - else - bound = @SVector[bounds[s÷2+1], (bounds[s÷2+1]+bounds[s÷2+2])/2] - if s == SEG && issparse[2] - bound = @SVector[bounds[s÷2+1], bounds[(s+2)÷2+1]] - end - init = 1 + (M+1)*N*(s-1) - push!(segment, bound[2]) - push!(segindex, 1 + (M+1)*N*(s)) - isopen = @SVector[false,s==SEG ? false : true] - g = UniLog{T}(bound,init,minterval,M,N,true,isopen) - push!(unilogs,g) - for idx=g.idx[1]:g.idx[2] - push!(grid, _grid(g, idx)) - end - end - end - end - - head, tail = grid[1], grid[end] - Isopen[1] && (grid[1] += eps(T) * 1e4) - Isopen[2] && (grid[end] -= eps(T) * 1e4) - checkOrder(grid) - return new{T,size,SEG}(grid, size, head, tail,unilogs, segment,segindex, Isopen) - end -end - -function Base.floor(grid::UniLogs{T,SIZE,SEG}, x) where {T,SIZE,SEG} -# (grid.head <= x <= grid.tail) || error("$x is out of the uniform grid range!") - segment = grid.segment - seg = 0 - for i=1:SEG - if x=0 -- N: each interval of the 'log' grid is further divided uniformly into N part, N>=1 -""" -@inline function tauUL(β, minterval, M::Int,N::Int, type=Float64) - seg = 2 - size = (M+1)*N*seg+1 - bounds = @SVector[0.0,β] - return UniLogs{Float64,size,seg}(bounds,minterval,M,N) -end - -@inline function tauUL(para, M::Int,N::Int, type=Float64) - return tauUL(para.β, para.τ_min, M, N, type) -end -""" - fermiK(Kf, maxK, halfLife, size::Int, kFi = floor(Int, 0.5size), type = Float64) - -Create a logarithmic fermionic K Grid, which is densest near the Fermi momentum ``k_F`` - -#Arguments: -- Kf: Fermi momentum -- maxK: the upper bound of the grid -- halfLife: the grid is densest in the range (Kf-halfLife, Kf+halfLife) -- size: the Grid size -- kFi: index of Kf -""" -@inline function fermiK(Kf, maxK, halfLife, size::Int, kFi=floor(Int, 0.5size), type=Float64) - size = Int(size) - c1 = Grid.Coeff{type}([0.0, Kf], [1.0, kFi], 1.0 / halfLife, false) - r1 = 1:kFi - - c2 = Grid.Coeff{type}([Kf, maxK], [kFi - 1.0, size], 1.0 / halfLife, true) - r2 = kFi + 1:size - return Log{type,size,2}([c1, c2], [r1, r2], [true, false]) -end - -""" - fermiKUL(Kf, maxK, minterval, M, N, type = Float64) - -Create a unilog fermionic K Grid, which is densest near the Fermi momentum `k_F` -the grid has 2 segments, 2*(M+1)*N+1 points. - -#Arguments: -- Kf: Fermi momentum -- maxK: the upper bound of the grid -- minterval: the minimal interval of the grid -- M: (M+1) is the number of the 'log' grid of each segment of the grids, M>=0 -- N: each interval of the 'log' grid is further divided uniformly into N part, N>=1 -""" - -@inline function fermiKUL(Kf, maxK, minterval, M::Int,N::Int, type=Float64) - seg = 2 - size = (M+1)*N*seg+1 - bounds = @SVector[0.0,Kf,maxK] - return UniLogs{Float64,size,seg}(bounds,minterval,M,N,[true,false],[true,true]) -end - -@inline function fermiKUL(para, M::Int,N::Int, type=Float64) - return fermiKUL(para.kF, para.k_max, para.k_min, M, N, type) -end - - - -""" - boseK(Kf, maxK, halfLife, size::Int, kFi = floor(Int, 0.5size), twokFi = floor(Int, 2size / 3), type = Float64) - -Create a logarithmic bosonic K Grid, which is densest near the momentum `0` and `2k_F` - -#Arguments: -- Kf: Fermi momentum -- maxK: the upper bound of the grid -- halfLife: the grid is densest in the range (0, Kf+halfLife) and (2Kf-halfLife, 2Kf+halfLife) -- size: the Grid size -- kFi: index of Kf -- twokFi: index of 2Kf -""" -@inline function boseK( - Kf, - maxK, - halfLife, - size::Int, - kFi=floor(Int, size / 3), - twokFi=floor(Int, 2size / 3), - type=Float64, -) - size = Int(size) - λ = 1.0 / halfLife - c1 = Grid.Coeff{type}([0.0, Kf], [1.0, kFi], λ, true) - r1 = 1:kFi - - c2 = Grid.Coeff{type}([Kf, 2.0 * Kf], [kFi - 1.0, twokFi], λ, false) - r2 = kFi + 1:twokFi - - c3 = Grid.Coeff{type}([2.0 * Kf, maxK], [twokFi - 1.0, size], λ, true) - r3 = twokFi + 1:size - - K = Grid.Log{type,size,3}([c1, c2, c3], [r1, r2, r3], [true, false]) -return K -end - -""" - boseKUL(Kf, maxK, minterval, M, N, type = Float64) - -Create a unilog bosonic K Grid, which is densest near the momentum `0` and `2k_F` -the grid has 3 segments, 3*(M+1)*N+1 points. - -#Arguments: -- Kf: Fermi momentum -- maxK: the upper bound of the grid -- minterval: the minimal interval of the grid -- M: (M+1) is the number of the 'log' grid of each segment of the grids, M>=0 -- N: each interval of the 'log' grid is further divided uniformly into N part, N>=1 -""" -@inline function boseKUL( - Kf, - maxK, - minterval, - M, - N, - type=Float64 -) - seg = 3 - size = (M+1)*N*seg+1 - bounds = @SVector[0.0,2Kf,maxK] - return UniLogs{Float64,size,seg}(bounds,minterval,M,N,[true,false],[false,true]) -end - -@inline function boseKUL(para, M::Int,N::Int, type=Float64) - return boseKUL(para.kF, para.k_max, para.k_min, M, N, type) -end - - -include("interpolate.jl") - -end diff --git a/test/grid.jl b/test/grid.jl index 8a9ffdf..faa409c 100644 --- a/test/grid.jl +++ b/test/grid.jl @@ -12,6 +12,12 @@ end @testset "ArbitraryGrid" begin + + # resulted grid should be sorted + println("Testing Arbitrary constructor, a warning should appear.") + arbitrary = SimpleGrid.Arbitrary{Float64}([0.0, 0.2, 0.6, 0.8, 0.4, 1.0]) + @test arbitrary.grid == [0.0, 0.2, 0.4, 0.6, 0.8, 1.0] + uniform = SimpleGrid.Uniform{Float64}([0.0, 1.0], 10) arbitrary = SimpleGrid.Arbitrary{Float64}(uniform.grid) println(arbitrary) @@ -121,7 +127,8 @@ end @testset "LogGrid" begin - loggrid = SimpleGrid.Log{Float64}([0.0, 1.0], 6, 0.001, true) + # loggrid = SimpleGrid.Log{Float64}([0.0, 1.0], 6, 0.001, true) + loggrid = SimpleGrid.Log{Float64}(bound=[0.0, 1.0], N=6, minterval=0.001, d2s=true) println(loggrid) # println(SimpleGrid.denseindex(loggrid)) @test floor(loggrid, 0.0) == 1 @@ -204,7 +211,13 @@ end @test floor(comp, comp[end]) == comp.size - 1 @test floor(comp, 1.0) == comp.size - 1 - comp = CompositeGrid.LogDensedGrid(:uniform, [0.0, 10.0], [0.0, 1.0, 1.0, 2.0, 2.000001], 4, 0.001, 4) + # comp = CompositeGrid.LogDensedGrid(:uniform, [0.0, 10.0], [0.0, 1.0, 1.0, 2.0, 2.000001], 4, 0.001, 4) + comp = CompositeGrid.LogDensedGrid(type=:uniform, + bound=[0.0, 10.0], + dense_at=[0.0, 1.0, 1.0, 2.0, 2.000001], + N=4, + minterval=0.001, + order=4) # println(comp.grid) # println(comp.inits) # println(CompositeGrid.denseindex(comp)) @@ -229,6 +242,13 @@ end # println(comp.grid) comp = CompositeGrid.LogDensedGrid(:uniform, [0.0, 10.0], [0.5, 10.0], 4, 0.001, 4) # println(comp.grid) + + # test edge case when two panel points are exactly minterval away + # raw panel will be [0.0, 5.0,5.001,5.002,10.0], + # after construction it will be [0.0, 5.001, 10.0] + comp = CompositeGrid.LogDensedGrid(:uniform, [0.0, 10.0], [5.0, 5.002], 4, 0.001, 4) + @test length(comp.panel) == 3 + @test isapprox(comp.panel[2], 5.001, atol=0.001) end end