added integration

src/grid/composite.jl

@@ -2,11 +2,9 @@ module CompositeGrid
 
 export LogDensedGrid, Composite
 
-using StaticArrays, FastGaussQuadrature
-include("simple.jl")
-using .SimpleGrid
+using StaticArrays, FastGaussQuadrature, CompositeGrids
 
-struct Composite{T<:AbstractFloat,PG,SG} <: ClosedGrid
+struct Composite{T<:AbstractFloat,PG,SG} <: SimpleGrid.ClosedGrid
     bound::SVector{2,T}
     size::Int
     grid::Vector{T}
@@ -46,7 +44,7 @@ struct Composite{T<:AbstractFloat,PG,SG} <: ClosedGrid
 end
 
 function Base.floor(grid::Composite{T,PG,SG}, x) where {T,PG,SG}
-    if SG<:ClosedGrid
+    if SG<:SimpleGrid.ClosedGrid
         i = floor(grid.panel, x)
         return grid.inits[i]-1+floor(grid.subgrids[i],x)
     end
@@ -67,16 +65,16 @@ Base.lastindex(grid::Composite) = grid.size
 
 function CompositeLogGrid(type, bound, N, minterval, d2s, order, T=Float64)
     if type == :cheb
-        SubGridType = BaryCheb{T}
+        SubGridType = SimpleGrid.BaryCheb{T}
     elseif type == :gauss
-        SubGridType = GaussLegendre{T}
+        SubGridType = SimpleGrid.GaussLegendre{T}
     elseif type == :uniform
-        SubGridType = Uniform{T}
+        SubGridType = SimpleGrid.Uniform{T}
     else
         error("$type not implemented!")
     end
 
-    panel = Log{T}(bound, N, minterval, d2s)
+    panel = SimpleGrid.Log{T}(bound, N, minterval, d2s)
     #println("logpanel:",panel.grid)
     subgrids = Vector{SubGridType}([])
 
@@ -85,17 +83,17 @@ function CompositeLogGrid(type, bound, N, minterval, d2s, order, T=Float64)
         push!(subgrids, SubGridType(_bound,order))
     end
 
-    return Composite{T, Log{T},SubGridType}(panel,subgrids)
+    return Composite{T, SimpleGrid.Log{T},SubGridType}(panel,subgrids)
 
 end
 
 function LogDensedGrid(type, bound, dense_on, N, minterval, order, T=Float64)
     if type == :cheb
-        SubGridType = BaryCheb{T}
+        SubGridType = SimpleGrid.BaryCheb{T}
     elseif type == :gauss
-        SubGridType = GaussLegendre{T}
+        SubGridType = SimpleGrid.GaussLegendre{T}
     elseif type == :uniform
-        SubGridType = Uniform{T}
+        SubGridType = SimpleGrid.Uniform{T}
     else
         error("$type not implemented!")
     end
@@ -161,15 +159,15 @@ function LogDensedGrid(type, bound, dense_on, N, minterval, order, T=Float64)
         push!(d2slist, true)
     end
 
-    panel = Arbitrary{T}(panelgrid)
+    panel = SimpleGrid.Arbitrary{T}(panelgrid)
     #println("panel:",panel.grid)
-    subgrids = Vector{Composite{T, Log{T}, SubGridType}}([])
+    subgrids = Vector{Composite{T, SimpleGrid.Log{T}, SubGridType}}([])
 
     for i in 1:length(panel.grid)-1
         push!(subgrids, CompositeLogGrid(type, [panel[i],panel[i+1]], N, minterval, d2slist[i], order))
     end
 
-    return Composite{T, Arbitrary{T},Composite{T, Log{T}, SubGridType}}(panel,subgrids)
+    return Composite{T, SimpleGrid.Arbitrary{T},Composite{T, SimpleGrid.Log{T}, SubGridType}}(panel,subgrids)
 
 end
 

src/grid/interpolate.jl

@@ -1,13 +1,22 @@
 module Interp
 
-using StaticArrays, FastGaussQuadrature
+using StaticArrays, FastGaussQuadrature, CompositeGrids
 
 include("chebyshev.jl")
 
-include("simple.jl")
-using .SimpleGrid
-include("composite.jl")
-using .CompositeGrid
+# include("simple.jl")
+# using .SimpleGrid
+# include("composite.jl")
+# using .CompositeGrid
+
+abstract type InterpStyle end
+struct FloorInterp <: InterpStyle end
+struct ChebInterp <: InterpStyle end
+struct CompositeInterp <: InterpStyle end
+
+InterpStyle(::Type) = FloorInterp()
+InterpStyle(::Type{<:SimpleGrid.BaryCheb}) = ChebInterp()
+InterpStyle(::Type{<:CompositeGrid.Composite}) = CompositeInterp()
 
 """
    linear1D(data,xgrid, x) 
@@ -46,44 +55,114 @@ linear interpolation of data(x)
     return gx
 end
 
-function interp1D(data, xgrid, x)
+function interp1D(data, xgrid::T, x) where {T}
+    interp1D(InterpStyle(T), data, xgrid, x)
+end
+
+function interp1D(::FloorInterp,data, xgrid, x)
     return linear1D(data, xgrid, x)
 end
 
-function interp1D(data, xgrid::BaryCheb, x)
+function interp1D(::ChebInterp, data, xgrid, x)
     return barycheb(xgrid.size, x, data, xgrid.weight, xgrid.grid)
 end
 
-function interp1D(data, xgrid::Composite, x)
+function interp1D(::CompositeInterp,data, xgrid, x)
     i = floor(xgrid.panel, x)
     head, tail = xgrid.inits[i], xgrid.inits[i]+xgrid.subgrids[i].size-1
     return interp1D(data[head:tail], xgrid.subgrids[i], x)
 end
 
-function interpGrid(data, xgrid, grid)
+function interpGrid(data, xgrid::T, grid) where {T}
+    interpGrid(InterpStyle(T), data, xgrid, grid)
+end
+
+function interpGrid(::Union{FloorInterp,ChebInterp}, data, xgrid, grid)
     ff = zeros(eltype(data), length(grid))
     for (xi, x) in enumerate(grid)
         ff[xi] = interp1D(data, xgrid, x)
     end
     return ff
 end
 
-
-function interpGrid(data, xgrid::Composite, grid)
+function interpGrid(::CompositeInterp, data, xgrid, grid)
     ff = zeros(eltype(data), length(grid))
 
     init, curr = 1, 1
-    for pi in xgrid.panel.size-1
+    for pi in 1:xgrid.panel.size-1
         if grid[init]< xgrid.panel[pi+1]
             head, tail = xgrid.inits[pi], xgrid.inits[pi]+xgrid.subgrids[pi].size-1
-            while grid[curr]<xgrid.panel[pi+1]
+            while grid[curr]<xgrid.panel[pi+1] && curr<length(grid)
                 curr += 1
             end
-            ff[init:curr-1] = interpGrid(data[head:tail], xgrid.subgrids[pi], grid[init:curr-1])
+            if grid[curr]<xgrid.panel[pi+1] && curr==length(grid)
+                ff[init:curr] = interpGrid(data[head:tail], xgrid.subgrids[pi], grid[init:curr])
+            else
+                ff[init:curr-1] = interpGrid(data[head:tail], xgrid.subgrids[pi], grid[init:curr-1])
+            end
+            # println(data[head:tail])
+            # println(xgrid.subgrids[pi].grid)
+            # println(grid[init:curr-1])
+            # println(ff[init:curr-1])
             init = curr
         end
     end
     return ff
 end
 
+abstract type IntegrateStyle end
+struct WeightIntegrate <: IntegrateStyle end
+struct NoIntegrate <: IntegrateStyle end
+struct CompositeIntegrate <: IntegrateStyle end
+
+IntegrateStyle(::Type) = NoIntegrate()
+IntegrateStyle(::Type{<:SimpleGrid.GaussLegendre}) = WeightIntegrate()
+IntegrateStyle(::Type{<:SimpleGrid.Uniform}) = WeightIntegrate()
+IntegrateStyle(::Type{<:SimpleGrid.Arbitrary}) = WeightIntegrate()
+IntegrateStyle(::Type{<:SimpleGrid.Log}) = WeightIntegrate()
+IntegrateStyle(::Type{<:CompositeGrid.Composite}) = CompositeIntegrate()
+
+function integrate1D(data, xgrid::T) where {T}
+    return integrate1D(IntegrateStyle(T), data, xgrid)
+end
+
+function integrate1D(::NoIntegrate, data, xgrid)
+    return 0.0
+    result = eltype(data)(0.0)
+
+    grid = xgrid.grid
+    for i in 1:xgrid.size
+        if i==1
+            weight = 0.5*(grid[2]-grid[1])
+        elseif i==xgrid.size
+            weight = 0.5*(grid[end]-grid[end-1])
+        else
+            weight = 0.5*(grid[i+1]-grid[i-1])
+        end
+        result += data[i]*weight
+    end
+    return result
+end
+
+function integrate1D(::WeightIntegrate, data, xgrid)
+    result = eltype(data)(0.0)
+
+    for i in 1:xgrid.size
+        result += data[i]*xgrid.weight[i]
+    end
+    return result
+end
+
+function integrate1D(::CompositeIntegrate, data, xgrid)
+    result = eltype(data)(0.0)
+
+    for pi in 1:xgrid.panel.size-1
+        head, tail = xgrid.inits[pi], xgrid.inits[pi]+xgrid.subgrids[pi].size-1
+        result += integrate1D( data[head:tail],xgrid.subgrids[pi])
+        currgrid = xgrid.subgrids[pi]
+    end
+    return result
+
+end
+
 end

src/grid/simple.jl

@@ -14,11 +14,22 @@ struct Arbitrary{T<:AbstractFloat} <: ClosedGrid
     bound::SVector{2,T}
     size::Int
     grid::Vector{T}
+    weight::Vector{T}
 
     function Arbitrary{T}(grid) where {T<:AbstractFloat}
         bound = [grid[1],grid[end]]
         size = length(grid)
-        return new{T}(bound, size, grid)
+        weight = similar(grid)
+        for i in 1:size
+            if i==1
+                weight[1] = 0.5*(grid[2]-grid[1])
+            elseif i==size
+                weight[end] = 0.5*(grid[end]-grid[end-1])
+            else
+                weight[i] = 0.5*(grid[i+1]-grid[i-1])
+            end
+        end
+        return new{T}(bound, size, grid, weight)
     end
 end
 
@@ -41,13 +52,23 @@ struct Uniform{T<:AbstractFloat} <: ClosedGrid
     bound::SVector{2,T}
     size::Int
     grid::Vector{T}
+    weight::Vector{T}
 
     function Uniform{T}(bound, N) where {T<:AbstractFloat}
         Ntot = N - 1
         interval = (bound[2]-bound[1])/Ntot
         grid = bound[1] .+ Vector(1:N) .* interval .- ( interval )
-
-        return new{T}(bound, N, grid)
+        weight = similar(grid)
+        for i in 1:N
+            if i==1
+                weight[1] = 0.5*(grid[2]-grid[1])
+            elseif i==N
+                weight[end] = 0.5*(grid[end]-grid[end-1])
+            else
+                weight[i] = 0.5*(grid[i+1]-grid[i-1])
+            end
+        end
+        return new{T}(bound, N, grid, weight)
     end
 end
 
@@ -115,6 +136,7 @@ struct Log{T<:AbstractFloat} <: ClosedGrid
     bound::SVector{2,T}
     size::Int
     grid::Vector{T}
+    weight::Vector{T}
 
     λ::T
     d2s::Bool
@@ -135,7 +157,18 @@ struct Log{T<:AbstractFloat} <: ClosedGrid
         end
         grid[1] = bound[1]
         grid[end] = bound[2]
-        return new{T}(bound, N, grid, λ, d2s)
+        weight = similar(grid)
+        for i in 1:N
+            if i==1
+                weight[1] = 0.5*(grid[2]-grid[1])
+            elseif i==N
+                weight[end] = 0.5*(grid[end]-grid[end-1])
+            else
+                weight[i] = 0.5*(grid[i+1]-grid[i-1])
+            end
+        end
+
+        return new{T}(bound, N, grid,weight, λ, d2s)
     end
 end
 

test/interpolate.jl

@@ -90,7 +90,7 @@
     end
 
     @testset "DensedLog" begin
-        β = π
+        β = 4
         tgrid = CompositeGrid.LogDensedGrid(:cheb, [0.0, β], [0.0, 0.5β, β], 4, 0.001, 4)
         # tugrid = Grid.Uniform{Float64,33}(0.0, β, (true, true))
         # kugrid = Grid.Uniform{Float64,33}(0.0, maxK, (true, true))
@@ -125,8 +125,10 @@
         @test abs(f(t) - fbar) < 3.e-6 # linear interpolation, so error is δK+δt
 
         tlist = rand(10) * β
+        tlist = sort(tlist)
+        println(tlist)
         ff = Interp.interpGrid(data, tgrid, tlist)
-        # println(tlist)
+
 
         for (ti, t) in enumerate(tlist)
             fbar = Interp.interp1D(data, tgrid, t)
@@ -136,5 +138,35 @@
         end
     end
 
+    @testset "Integrate" begin
+        β = 1.0
+        tgrid = SimpleGrid.GaussLegendre{Float64}([0.0, β], 4)
+        # tugrid = Grid.Uniform{Float64,33}(0.0, β, (true, true))
+        # kugrid = Grid.Uniform{Float64,33}(0.0, maxK, (true, true))
+        f(t) = t
+        data = zeros(tgrid.size)
+        for (ti, t) in enumerate(tgrid.grid)
+            data[ti] = f(t)
+        end
+        println(tgrid.grid)
+        println(data)
+        println(tgrid.weight)
+        println(sum(data.*tgrid.weight))
+        int_result = Interp.integrate1D(data, tgrid)
+        @test abs(int_result - 0.5) < 3.e-6
+
+        β = 1.0
+        tgrid = CompositeGrid.LogDensedGrid(:gauss, [0.0, β], [0.0, 0.5β, β], 2, 0.001, 3)
+        # tugrid = Grid.Uniform{Float64,33}(0.0, β, (true, true))
+        # kugrid = Grid.Uniform{Float64,33}(0.0, maxK, (true, true))
+        data = zeros(tgrid.size)
+        for (ti, t) in enumerate(tgrid.grid)
+            data[ti] = f(t)
+        end
+        println(tgrid.grid)
+        println(data)
+        int_result = Interp.integrate1D(data, tgrid)
+        @test abs(int_result - 0.5) < 3.e-6
+    end
 end