Merge branch 'master' into scicade

.travis/setup.sh

@@ -3,7 +3,3 @@ set -ev
 if [[ $TRAVIS_BRANCH == 'development' ]]; then
   julia -e 'Pkg.checkout("BasisFunctions","development")'
 fi
-
-if [[ $TRAVIS_BRANCH == 'residual' ]]; then
-  julia -e 'Pkg.checkout("BasisFunctions","residual")'
-fi

src/FrameFun.jl

@@ -84,14 +84,15 @@ import BasisFunctions: AbstractSubGrid, IndexSubGrid, is_subindex, supergrid,
     similar_subgrid
 
 import BasisFunctions: Gram, DualGram, MixedGram, DiscreteGram, DiscreteDualGram, DiscreteMixedGram
+import BasisFunctions: dual
 
 import BasisFunctions: discrete_approximation_operator, continuous_approximation_operator
 
 ###############################
 ## Exhaustive list of exports
 ###############################
 # from funs.jl
-export ExpFun, ChebyFun, Fun, SetFun, sampling_grid, domain, abserror
+export ExpFun, ChebyFun, Fun, SetFun, sampling_grid, domain, abserror, maxerror, L2error
 
 # from subgrid.jl
 export MaskedGrid
@@ -146,6 +147,7 @@ include("fun/funs.jl")
 
 
 include("approximation/fe_solvers.jl")
+include("approximation/continuous_solver.jl")
 include("approximation/lowranksolver.jl")
 include("approximation/fastsolver.jl")
 include("approximation/smoothsolver.jl")

src/approximation/continuous_solver.jl

@@ -0,0 +1,36 @@
+
+function continuous_approximation_operator(dest::ExtensionFrame; sampling_factor=1, solver=FrameFun.FE_DirectSolver, options...)
+    # since the other one is not very efficient (this one isnt either), concider this not as a general case
+    (sampling_factor ≈ 1) &&
+        (return ContinuousSolverPlan(solver(MixedGram(dest; options...); options...), continuous_normalization(dest; options...)))
+    src = resize(dest, round(Int, sampling_factor*length(dest)))
+    println(sampling_factor)
+    ContinuousSolverPlan(solver(MixedGram(dest, src; options...); options...), continuous_normalization(src; options...))
+end
+
+continuous_normalization(set::FunctionSet; options...) = DualGram(set; options...)
+continuous_normalization(frame::ExtensionFrame; options...) = DualGram(basis(frame); options...)
+
+immutable ContinuousSolverPlan{T} <: AbstractOperator{T}
+    src                     :: FunctionSet
+    dest                    :: FunctionSet
+    mixedgramsolver         :: FE_Solver
+    normalizationofb        :: AbstractOperator
+
+    scratch                 :: Vector{T}
+    mixedgram               :: AbstractOperator
+    ContinuousSolverPlan{T}(src::FunctionSet, dest::FunctionSet, mixedgramsolver::FE_Solver, normalizationofb::AbstractOperator) where {T} =
+        new(src, dest, mixedgramsolver, normalizationofb, zeros(T, length(src)), mixedgramsolver.op)
+end
+
+ContinuousSolverPlan(solver::FE_Solver{T}, normalization::AbstractOperator) where {T} =
+    ContinuousSolverPlan(src(solver), dest(solver), solver, normalization)
+
+ContinuousSolverPlan(src::FunctionSet, dest::FunctionSet, solver::FE_Solver{T}, normalization::AbstractOperator) where {T} =
+    ContinuousSolverPlan{T}(src, dest, solver, normalization)
+
+function apply!(s::ContinuousSolverPlan, coef_dest, coef_src)
+    apply!(s.normalizationofb, s.scratch, coef_src)
+    coef_dest[:] = s.mixedgramsolver*s.scratch
+    # println(norm(s.mixedgram*coef_dest-s.scratch))
+end

src/approximation/fastsolver.jl

@@ -27,7 +27,7 @@ struct FE_ProjectionSolver{ELT} <: FE_Solver{ELT}
     end
 end
 
-FE_ProjectionSolver(op::AbstractOperator, scaling; options...) =
+FE_ProjectionSolver(op::AbstractOperator; scaling=nothing, options...) =
     FE_ProjectionSolver{eltype(op)}(op, scaling; options...)
 
 function plunge_operator(op, scaling)

src/approximation/fe_approx.jl

@@ -55,12 +55,16 @@ function oversampled_evaluation_operator(S::FunctionSet, D::Domain; sampling_fac
         BG = boundary(grid(lB), D)
         op = [op; grid_evaluation_operator(S,gridspace(B,BG),BG)]
     end
-    (op,length(lB))
+    (op,scaling_factor(lB))
 end
 
+scaling_factor(S::FunctionSet) = length(S)
+scaling_factor(S::DerivedSet) = scaling_factor(superset(S))
+scaling_factor(S::ChebyshevBasis) = length(S)/2
+
 function discrete_approximation_operator(set::ExtensionFrame; solver = default_frame_solver(domain(set), basis(set)), options...)
     (op, scaling) = oversampled_evaluation_operator(basis(set),domain(set);options...)
-    solver(op, scaling; options...)
+    solver(op; scaling=scaling, options...)
 end
 
 primarybasis(set::FunctionSet) = set
@@ -71,15 +75,15 @@ end
 struct FE_BestSolver
 end
 
-function FE_BestSolver(op::AbstractOperator, scaling; verbose= false, options...)
+function FE_BestSolver(op::AbstractOperator; scaling=NaN, verbose= false, options...)
     if has_transform(src(op))
         R = estimate_plunge_rank(op)
         if R < size(op, 2)/2
             verbose && println("Estimated plunge rank $R smaller than $(size(op,2))/2 -> projection solver ")
-            FE_ProjectionSolver(op, scaling; verbose=verbose,options...)
+            FE_ProjectionSolver(op; scaling=scaling, verbose=verbose,options...)
         else
             verbose && println("Estimated plunge rank $R greater than $(size(op,2))/2 -> direct solver ")
-            FE_DirectSolver(op, scaling; verbose=verbose,options...)
+            FE_DirectSolver(op; verbose=verbose,options...)
         end
     else
         # Don't bother with a fast algorithm if there is no fast transform

src/approximation/fe_solvers.jl

@@ -17,59 +17,22 @@ src(s::FE_Solver) = dest(op(s))
 
 dest(s::FE_Solver) = src(op(s))
 
-
-
 struct FE_DirectSolver{ELT} <: FE_Solver{ELT}
     op      ::  AbstractOperator
     QR      ::  Factorization
 
-    function FE_DirectSolver{ELT}(op::AbstractOperator,scaling) where ELT
+    function FE_DirectSolver{ELT}(op::AbstractOperator) where ELT
         new(op, qrfact(matrix(op),Val{true}))
     end
 end
 
-FE_DirectSolver{ELT}(op::AbstractOperator{ELT}, scaling; options...) =
-    FE_DirectSolver{eltype(op)}(op,scaling)
+FE_DirectSolver{ELT}(op::AbstractOperator{ELT}; options...) =
+    FE_DirectSolver{eltype(op)}(op)
 
 function apply!(s::FE_DirectSolver, coef_dest, coef_src)
     coef_dest[:] = s.QR \ coef_src
 end
 
-immutable ContinuousDirectSolver{T} <: AbstractOperator{T}
-  src                     :: FunctionSet
-  mixedgramfactorization  :: Factorization
-  normalizationofb        :: AbstractOperator
-  scratch                 :: Array{T,1}
-end
-
-dest(s::ContinuousDirectSolver) = s.src
-
-ContinuousDirectSolver(frame::ExtensionFrame; options...) =
-    ContinuousDirectSolver{eltype(frame)}(frame, qrfact(matrix(MixedGram(frame; options...)),Val{true}), DualGram(basis(frame); options...), zeros(eltype(frame),length(frame)))
-
-function apply!(s::ContinuousDirectSolver, coef_dest, coef_src)
-  apply!(s.normalizationofb, s.scratch, coef_src)
-  coef_dest[:] = s.mixedgramfactorization \ s.scratch
-end
-
-immutable ContinuousTruncatedSolver{T} <: AbstractOperator{T}
-  src                     :: FunctionSet
-  mixedgramsvd            :: AbstractOperator
-  normalizationofb        :: AbstractOperator
-  scratch                 :: Array{T,1}
-end
-
-dest(s::ContinuousTruncatedSolver) = s.src
-
-function ContinuousTruncatedSolver(frame::ExtensionFrame; cutoff=1e-5, fullsvd=false, options...)
-    fullsvd? solver = FrameFun.ExactTruncatedSvdSolver: solver = FrameFun.TruncatedSvdSolver
-    ContinuousTruncatedSolver{eltype(frame)}(frame, solver(MixedGram(frame; options...); cutoff=cutoff, options...), DualGram(basis(frame); options...), zeros(eltype(frame),length(frame)))
-end
-
-function apply!(s::ContinuousTruncatedSolver, coef_dest, coef_src)
-  apply!(s.normalizationofb, s.scratch, coef_src)
-  coef_dest[:] = s.mixedgramsvd*s.scratch
-end
 ## abstract FE_IterativeSolver <: FE_Solver
 
 

src/approximation/lowranksolver.jl

@@ -6,7 +6,7 @@ If A is MxN, and of rank R, the cost of constructing this solver is MR^2.
 
 For tensor product operators it returns a decomposition of the linearized system
 """
-struct TruncatedSvdSolver{ELT} <: AbstractOperator{ELT}
+struct TruncatedSvdSolver{ELT} <: FE_Solver{ELT}
     # Keep the original operator
     op          ::  AbstractOperator
     # Random matrix
@@ -21,37 +21,43 @@ struct TruncatedSvdSolver{ELT} <: AbstractOperator{ELT}
     scratch_src ::  Array{ELT,1}
     scratch_dest::  Array{ELT,1}
     smallcoefficients :: Bool
+    smalltol :: Float64
 
-    function TruncatedSvdSolver{ELT}(op::AbstractOperator; cutoff = default_cutoff(problem), R = 5, growth_factor = sqrt(2), verbose = false, smallcoefficients=false,options...) where ELT
+
+    function TruncatedSvdSolver{ELT}(op::AbstractOperator; cutoff = default_cutoff(op), R = 5, growth_factor = 2, verbose = false, smallcoefficients=false,smalltol=10,options...) where ELT
         finished=false
         USV = ()
         R = min(R, size(op,2))
         random_matrix = map(ELT, rand(size(op,2), R))
         C = apply_multiple(op, random_matrix)
         c = cond(C)
-        # TODO change things with cold back
-        # cold = 1
+# <<<<<<< HEAD
+#         # TODO change things with cold back
+#         # cold = 1
+#         m = maximum(abs.(C))
+#         while (c < 1/cutoff) && (R<size(op,2)) #&& (c/cold > 10)
+# =======
+        cold = cutoff
         m = maximum(abs.(C))
-        while (c < 1/cutoff) && (R<size(op,2)) #&& (c/cold > 10)
+        while (c < 1/cutoff) && (R<size(op,2)) && (c>cold*10)
+            verbose && println("c : $c\t cold : $cold\t cutoff : $cutoff")
             verbose && println("Solver truncated at R = ", R, " dof out of ",size(op,2))
             R0 = R
             R = min(round(Int,growth_factor*R),size(op,2))
-            verbose && println("Solver truncated at R = ", R, " dof out of ",size(op,2))
-            verbose && println(c," ",m," ",cutoff)
             extra_random_matrix = map(ELT, rand(size(op,2), R-R0))
             Cextra = apply_multiple(op, extra_random_matrix)
             random_matrix = [random_matrix extra_random_matrix]
             # Extra condition: condition number has to increase by a significant amount each step, otherwise, possibly well conditioned.
             # cold = c
             C = [C Cextra]
             c = cond(C)
-            m = maximum(abs.(C))
+
         end
+        verbose && println("c : $c\t cold : $cold\t cutoff : $cutoff")
         verbose && println("Solver truncated at R = ", R, " dof out of ",size(op,2))
-            verbose && println(c," ",m," ",cutoff)
         USV = LAPACK.gesdd!('S',C)
         S = USV[2]
-        maxind = findlast(S.>cutoff*m)
+        maxind = findlast(S.>(maximum(S)*cutoff))
         Sinv = 1./S[1:maxind]
         y = zeros(ELT, size(USV[3],1))
         sy = zeros(ELT, maxind)
@@ -61,7 +67,7 @@ struct TruncatedSvdSolver{ELT} <: AbstractOperator{ELT}
 
         scratch_src = zeros(ELT, length(dest(op)))
         scratch_dest = zeros(ELT, length(src(op)))
-        new(op, W, USV[1][:,1:maxind]',Sinv,USV[3][1:maxind,:]',y,sy, scratch_src, scratch_dest,smallcoefficients)
+        new(op, W, USV[1][:,1:maxind]',Sinv,USV[3][1:maxind,:]',y,sy, scratch_src, scratch_dest,smallcoefficients,smalltol)
     end
 end
 
@@ -80,7 +86,7 @@ function apply!(s::TruncatedSvdSolver, coef_dest::Vector, coef_src::Vector)
         s.sy[i]=s.sy[i]*s.Sinv[i]
     end
     if s.smallcoefficients
-        s.sy[abs(s.sy).>100*maximum(abs(coef_src))]=0
+        s.sy[abs.(s.sy).>s.smalltol*maximum(abs.(coef_src))]=0
     end
     A_mul_B!(s.y, s.V, s.sy)
 
@@ -94,7 +100,7 @@ function apply!(s::TruncatedSvdSolver, coef_dest, coef_src)
         s.sy[i]=s.sy[i]*s.Sinv[i]
     end
     if s.smallcoefficients
-        s.sy[abs(s.sy).>100*maximum(abs(coef_src))]=0
+        s.sy[abs.(s.sy).>s.smalltol*maximum(abs.(coef_src))]=0
     end
     A_mul_B!(s.y, s.V, s.sy)
     apply!(s.W, s.scratch_dest, s.y)
@@ -108,7 +114,7 @@ The cost is equal to the computation of the SVD of A.
 
 For tensor product operators it returns a decomposition of the linearized system
 """
-struct ExactTruncatedSvdSolver{ELT} <: AbstractOperator{ELT}
+struct ExactTruncatedSvdSolver{ELT} <: FE_Solver{ELT}
     # Keep the original operator
     op          ::  AbstractOperator
     # Decomposition
@@ -120,7 +126,7 @@ struct ExactTruncatedSvdSolver{ELT} <: AbstractOperator{ELT}
     sy          ::  Array{ELT,1}
     scratch_src ::  Array{ELT,1}
 
-    function ExactTruncatedSvdSolver{ELT}(op::AbstractOperator; cutoff = default_cutoff(problem), verbose = false,options...) where ELT
+    function ExactTruncatedSvdSolver{ELT}(op::AbstractOperator; cutoff = default_cutoff(op), verbose = false,options...) where ELT
         C = matrix(op)
         m = maximum(abs.(C))
         USV = LAPACK.gesdd!('S',C)

src/approximation/smoothsolver.jl

@@ -22,8 +22,8 @@ struct FE_SmoothProjectionSolver{ELT} <: FE_Solver{ELT}
     function FE_SmoothProjectionSolver{ELT}(op::AbstractOperator, scaling; cutoff = default_cutoff(op), cutoffv=sqrt(cutoff), R = estimate_plunge_rank(op), verbose=false,  options...) where ELT
         plunge_op = plunge_operator(op, scaling)
         # Create Random matrices
-        TS1 = TruncatedSvdSolver(plunge_op*op; cutoff = cutoff, options...)
-        TS2 = TruncatedSvdSolver(op'*plunge_op; cutoff = cutoffv, options...)
+        TS1 = TruncatedSvdSolver(plunge_op*op; cutoff = cutoff, verbose=verbose,R=R,options...)
+        TS2 = TruncatedSvdSolver(op'*plunge_op; cutoff = cutoffv, verbose=verbose, R=R, options...)
         # D = Sobolev operator
         D = IdxnScalingOperator(src(op); options...)
         AD = inv(D)
@@ -42,7 +42,7 @@ struct FE_SmoothProjectionSolver{ELT} <: FE_Solver{ELT}
 end
 
 
-FE_SmoothProjectionSolver(op::AbstractOperator, scaling; options...) =
+FE_SmoothProjectionSolver(op::AbstractOperator; scaling=nothing, options...) =
         FE_SmoothProjectionSolver{eltype(op)}(op, scaling; options...)
 
 function apply!(s::FE_SmoothProjectionSolver, destarg, src, coef_dest, coef_src)

src/diffequation.jl

@@ -20,8 +20,9 @@ end
 struct DirichletBC
     dRhs   :: Function
     D      :: Domain
-    function DirichletBC(dRhs=default_boundary_condition :: Function, D=FullSpace())
-        new(dRhs,D)
+    factor :: Number
+    function DirichletBC(dRhs=default_boundary_condition :: Function, D=FullSpace(), factor=1.0)
+        new(dRhs,D,factor)
     end
 end
 
@@ -42,7 +43,7 @@ default_boundary_condition(x,y,z) = 0
 
 function operator(BC :: DirichletBC, S::FunctionSet, G::AbstractGrid, D::Domain)
     G = subgrid(G,BC.D)
-    grid_evaluation_operator(S,DiscreteGridSpace(G,eltype(S)),G)
+    BC.factor*grid_evaluation_operator(S,DiscreteGridSpace(G,eltype(S)),G)
 end
 
 function operator(BC :: NeumannBC, S::FunctionSet{2}, G::AbstractGrid{2}, D::Domain2d)
@@ -124,8 +125,8 @@ function rhs(D::DiffEquation; incboundary = false, options...)
     if incboundary
         push!(rhs,sample(BG,D.DRhs, eltype(src(op))))
     end
-    for BC in D.BCs
-        push!(rhs,sample(subgrid(BG,BC.D),BC.dRhs, eltype(src(op))))
+    for i = 1:length(D.BCs)
+        push!(rhs,sample(subgrid(BG,D.BCs[i].D),D.BCs[i].dRhs, eltype(src(op))))
     end
     MultiArray(rhs)
 end
@@ -135,8 +136,7 @@ function solve(D::DiffEquation, solver=FE_ProjectionSolver; options...)
     Adiff= inv_diagonal(D.Diff)
     b = rhs(D; options...)
     OP = operator(D; options...)
-    A = solver(OP, length(lB); options...)
+    A = solver(OP; scaling=length(lB), options...)
     coef  = A * b
     SetFun(D.D, dest(A), Adiff*coef)
 end
-

src/domains/extensions.jl

@@ -101,8 +101,8 @@ function boundingbox(d::UnionDomain)
 end
 
 function boundingbox(d::IntersectionDomain)
-    left = SVector(maximum(hcat(map(leftendpoint,elements(d))...),2))
-    right = SVector(minimum(hcat(map(rightendpoint,elements(d))...),2))
+    left = SVector(maximum(hcat(map(leftendpoint,elements(d))...),2)...)
+    right = SVector(minimum(hcat(map(rightendpoint,elements(d))...),2)...)
     boundingbox(left,right)
 end