From d48a6bae3751c957e855edd967bf303be898dda3 Mon Sep 17 00:00:00 2001 From: Xiansheng Cai Date: Wed, 26 Jul 2023 15:44:59 -0400 Subject: [PATCH] rm old --- README.md | 5 +- docs/src/README.md | 7 +- src/old/chebyshev.jl | 72 --------------- src/old/interpolate.jl | 194 ---------------------------------------- src/old/old_com_grid.jl | 111 ----------------------- src/old/pgrid.jl | 115 ------------------------ 6 files changed, 9 insertions(+), 495 deletions(-) delete mode 100644 src/old/chebyshev.jl delete mode 100644 src/old/interpolate.jl delete mode 100644 src/old/old_com_grid.jl delete mode 100644 src/old/pgrid.jl diff --git a/README.md b/README.md index bc1b88f..ab91dca 100644 --- a/README.md +++ b/README.md @@ -91,7 +91,10 @@ The `SimpleGrid` module in `CompositeGrids.jl` offers various basic grids that s - **Log Grid:** Defined by the boundary, number of grid points, minimum interval, and direction. It generates a log-dense grid based on the provided parameters. An O(1) floor function is provided. For example: ```julia using CompositeGrids - loggrid = SimpleGrid.Log{Float64}([0.0,1.0], 6, 0.0001, true) + loggrid = SimpleGrid.Log{Float64}( + bound=[0.0,1.0], N=6, + minterval=0.0001, + d2s = true) # dense to sparse println(loggrid.grid) ``` ``` diff --git a/docs/src/README.md b/docs/src/README.md index 1a16d06..ab91dca 100644 --- a/docs/src/README.md +++ b/docs/src/README.md @@ -1,6 +1,6 @@ [![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) +[![img](https://codecov.io/gh/numericalEFT/CompositeGrids.jl/branch/master/graph/badge.svg?token=WN6HO1XASY)](https://codecov.io/gh/numericaleft/CompositeGrids.jl) # Introduction @@ -91,7 +91,10 @@ The `SimpleGrid` module in `CompositeGrids.jl` offers various basic grids that s - **Log Grid:** Defined by the boundary, number of grid points, minimum interval, and direction. It generates a log-dense grid based on the provided parameters. An O(1) floor function is provided. For example: ```julia using CompositeGrids - loggrid = SimpleGrid.Log{Float64}([0.0,1.0], 6, 0.0001, true) + loggrid = SimpleGrid.Log{Float64}( + bound=[0.0,1.0], N=6, + minterval=0.0001, + d2s = true) # dense to sparse println(loggrid.grid) ``` ``` diff --git a/src/old/chebyshev.jl b/src/old/chebyshev.jl deleted file mode 100644 index 9f7502c..0000000 --- a/src/old/chebyshev.jl +++ /dev/null @@ -1,72 +0,0 @@ -""" -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 -""" -function barychebinit(n) - x = zeros(Float64, 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 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 -""" -function barycheb(n, x, f, wc, xc) - for j in 1:n - if x == xc[j] - return f[j] - end - end - - num, den = 0.0, 0.0 - for j in 1:n - q = wc[j] / (x - xc[j]) - num += q * f[j] - den += q - end - return num / den -end - -function barycheb2(n, x, f, wc, xc) - for j in 1:n - if x == xc[j] - return f[j, :] - end - end - - den = 0.0 - num = zeros(eltype(f), size(f)[2]) - for j in 1:n - q = wc[j] / (x - xc[j]) - num += q .* f[j, :] - den += q - end - return num ./ den -end \ No newline at end of file diff --git a/src/old/interpolate.jl b/src/old/interpolate.jl deleted file mode 100644 index b5565fc..0000000 --- a/src/old/interpolate.jl +++ /dev/null @@ -1,194 +0,0 @@ -using Statistics -""" - linear1D(data,xgrid, x) - -linear interpolation of data(x) - -#Arguments: -- xgrid: one-dimensional grid of x -- data: one-dimensional array of data -- x: x -""" - -@inline function linear1D(data, xgrid, x) - - xarray = xgrid.grid - - xi0,xi1 = 0,0 - if(x<=xarray[firstindex(xgrid)]) - xi0=1 - xi1=2 - elseif(x>=xarray[lastindex(xgrid)]) - xi0=lastindex(xgrid)-1 - xi1=xi0+1 - else - xi0=floor(xgrid,x) - xi1=xi0+1 - end - - dx0, dx1 = x - xarray[xi0], xarray[xi1] - x - - d0, d1 = data[xi0], data[xi1] - - g = d0 * dx1 + d1 * dx0 - - gx = g / (dx0 + dx1) - return gx -end - - -""" - linear2D(data, xgrid, ygrid, x, y) - -linear interpolation of data(x, y) - -#Arguments: -- xgrid: one-dimensional grid of x -- ygrid: one-dimensional grid of y -- data: two-dimensional array of data -- x: x -- y: y -""" -@inline function linear2D(data, xgrid, ygrid, x, y) - - xarray, yarray = xgrid.grid, ygrid.grid - - # if ( - # x <= xarray[firstindex(xgrid)] || - # x >= xarray[lastindex(xgrid)] || - # y <= yarray[firstindex(ygrid)] || - # y >= yarray[lastindex(ygrid)] - # ) - # return 0.0 - # end - # xi0, yi0 = floor(xgrid, x), floor(ygrid, y) - # xi1, yi1 = xi0 + 1, yi0 + 1 - - xi0,xi1,yi0,yi1 = 0,0,0,0 - if(x<=xarray[firstindex(xgrid)]) - xi0=1 - xi1=2 - elseif(x>=xarray[lastindex(xgrid)]) - xi0=lastindex(xgrid)-1 - xi1=xi0+1 - else - xi0=floor(xgrid,x) - xi1=xi0+1 - end - - if(y<=yarray[firstindex(ygrid)]) - yi0=1 - yi1=2 - elseif(y>=yarray[lastindex(ygrid)]) - yi0=lastindex(ygrid)-1 - yi1=yi0+1 - else - yi0=floor(ygrid,y) - yi1=yi0+1 - end - - dx0, dx1 = x - xarray[xi0], xarray[xi1] - x - dy0, dy1 = y - yarray[yi0], yarray[yi1] - y - - d00, d01 = data[xi0, yi0], data[xi0, yi1] - d10, d11 = data[xi1, yi0], data[xi1, yi1] - - g0 = data[xi0, yi0] * dx1 + data[xi1, yi0] * dx0 - g1 = data[xi0, yi1] * dx1 + data[xi1, yi1] * dx0 - - gx = (g0 * dy1 + g1 * dy0) / (dx0 + dx1) / (dy0 + dy1) - return gx -end - - -""" -testInterpolation1D(func, grid1, grid2) - -Test interpolation on grid1 with function func, compare results with func(grid2). -Return max deviation and std of deviation. - -#Arguments: -- func: function needs to be interpolated -- grid1: grid to be tested -- grid2: grid to be interpolated -""" -function testInterpolation1D(func, grid1, grid2, isRelativeError=false, find=false) - data1 = [func(grid1[i]) for i in 1:(grid1.size)] - - data2 = [func(grid2[i]) for i in 1:(grid2.size)] - data2_int = [linear1D(data1, grid1, grid2[i]) for i in 1:(grid2.size)] - - d = data2-data2_int - if isRelativeError - d = abs.(d ./ maximum( abs.(data2) )) - else - d = abs.( d ) - end - # d_max = maximum( d ) - d_max, d_max_index = findmax(d) - d_std = std(d) - if find - return d_max, d_std, d_max_index - else - return d_max, d_std - end -end - -# function testInterpolation1D_rel(func, grid1, grid2) -# data1 = [func(grid1[i]) for i in 1:(grid1.size)] - -# data2 = [func(grid2[i]) for i in 1:(grid2.size)] -# data2_int = [linear1D(data1, grid1, grid2[i]) for i in 1:(grid2.size)] - - # d_rel = abs.((data2-data2_int) ./ data2) -# d_max = maximum(d_rel) -# d_std = std(d_rel) -# return d_max, d_std -# end - -""" -optimizeUniLog(generator, para, MN, func) - -Return a pair of parameters (M, N) which give optimized unilog grid -generated by generator that minimize absolute error when interpolate func. - -#Arguments: -- generator: function generator(para, M, N) that generates unilog grid. -- para: parameters needed by generator -- MN: the result is constrained by (M+1)*N 1 - if floor(Int, MN/(Nb+1))==Na - # maximize Nb for current Na - Nb = Nb + 1 - continue - end - - # test (Na-1, Nb) and (Nb-1, Na) - d1 = f(Na-1, Nb) - d2 = f(Nb-1, Na) - - if d1<=d2 && d1 k.panel[kpidx + 1] - # if q is too large, move k panel to the next - # println("before $q, $kpidx, $(k.panel[kpidx]) -> $(k.panel[kpidx + 1])") - while q > kgrid.panel[kpidx + 1] - kpidx += 1 - end - # println("after $q, $kpidx, $(k.panel[kpidx]) -> $(k.panel[kpidx + 1])") - head, tail = idx(kpidx, 1, order), idx(kpidx, order, order) - fx = @view f[head:tail] # all F in the same kpidx-th K panel - x = @view k.grid[head:tail] - w = @view k.wgrid[head:tail] - @assert kpidx <= k.Np - end - ff[qi] = barycheb(order, q, fx, w, x) # the interpolation is independent with the panel length - # @assert k.panel[kpidx - 1][order] <= q <= k.panel[pidx + 1][1] "$q for kpidx=$kpidx with $x" - end - return ff -end - -if abspath(PROGRAM_FILE) == @__FILE__ - grid=CompositeGrid([0.0, 1.0, 2.0], 4, :cheb) - println(grid.grid) -end diff --git a/src/old/pgrid.jl b/src/old/pgrid.jl deleted file mode 100644 index 91d5911..0000000 --- a/src/old/pgrid.jl +++ /dev/null @@ -1,115 +0,0 @@ -using CompositeGrids: Grid -using StaticArrays - -# -# generate a p grid for integration. dense on k and kf. -# - -struct PGrid{T<:AbstractFloat, SIZE} - k::T - kF::T - iscombined::Int - grid::MVector{SIZE,T} - size::Int - head::T - tail::T - unilogs::SVector{4,Grid.UniLog{T}} - - function PGrid{T,SIZE}(k, kF, maxk, mink, MS) where {T<:AbstractFloat, SIZE} - size = sum(MS)+1 - - grid, segment,segindex = [],[],[] - unilogs = [] - iscombined = 0 - - if abs(k-kF)<2mink - k1 = (k+kF)/2.0 - iscombined = 2 - M1 = MS[1]+MS[2]-1 - M2 = MS[3]+MS[4]-1 - - g1 = Grid.UniLog{T}([0.0, k1], 1, mink, M1, 1, false, [false, true]) - push!(unilogs,g1) - push!(unilogs,g1) - for idx=g1.idx[1]:g1.idx[2] - push!(grid, Grid._grid(g1, idx)) - end - - g2 = Grid.UniLog{T}([k1, maxk], M1+1, mink, M2, 1, true, [false, false]) - push!(unilogs,g2) - push!(unilogs,g2) - for idx=g2.idx[1]:g2.idx[2] - push!(grid, Grid._grid(g2, idx)) - end - - elseif abs(k)<2mink - iscombined = 1 - - g0 = Grid.UniLog{T}([0.0, kF/2.0], 1, mink, MS[1]+MS[2]-1, 1, true, [false, true]) - push!(unilogs,g0) - push!(unilogs,g0) - for idx=g0.idx[1]:g0.idx[2] - push!(grid, Grid._grid(g0, idx)) - end - - g1 = Grid.UniLog{T}([kF/2.0, kF], MS[1]+MS[2]+1, mink, MS[3]-1, 1, false, [false, true]) - push!(unilogs,g1) - for idx=g1.idx[1]:g1.idx[2] - push!(grid, Grid._grid(g1, idx)) - end - - g2 = Grid.UniLog{T}([kF, maxk], MS[1]+MS[2]+MS[3]+1, mink, MS[4]-1, 1, true, [false, false]) - push!(unilogs,g2) - for idx=g2.idx[1]:g2.idx[2] - push!(grid, Grid._grid(g2, idx)) - end - - else - k1, k2 = minimum([k,kF]), maximum([k,kF]) - bounds = [[0.0, k1],[k1,(k1+k2)/2.0], [(k1+k2)/2.0, k2],[k2, maxk]] - d2s = [false, true, false, true] - isopen = [[false, true],[false, true],[false, true],[false, false]] - init = 1 - - for i in 1:4 - g = Grid.UniLog{T}(bounds[i], init, mink, MS[i]-1, 1, d2s[i], isopen[i]) - push!(unilogs, g) - for idx=g.idx[1]:g.idx[2] - push!(grid, Grid._grid(g, idx)) - end - init += MS[i] - end - - end - - head, tail = grid[1], grid[end] - return new{T,size}(k, kF, iscombined, grid, size, head, tail, unilogs) - end - -end - - -""" - function pGrid(k, kF, maxk, mink, MS) - -Create a logarithmic p Grid, which is densest near the Fermi momentum ``k_F`` and a given momentum ``k`` - -#Arguments: -- k: a given momentum -- kF: Fermi momentum -- maxK: the upper bound of the grid -- mink: the minimum interval of the grid -- MS: a list of # of points in each segments, the total size will be sum(MS)+1 -""" -@inline function pGrid(k, kF, maxk, mink, MS) - return PGrid{Float64, sum(MS)+1}(k, kF, maxk, mink, MS) -end - -if abspath(PROGRAM_FILE) == @__FILE__ - pgrid = pGrid(0.001, 1.0, 10.0, 0.001, [4,4,4,4]) - println(pgrid.grid) - pgrid = pGrid(0.999, 1.0, 10.0, 0.001, [4,4,4,4]) - println(pgrid.grid) - pgrid = pGrid(2.0, 1.0, 10.0, 0.001, [4,4,4,4]) - println(pgrid.grid) -end