add chunkerflam utility, passes some limited tests

chunkie/+chnk/+flam/nproxy_square.m

@@ -9,23 +9,26 @@
 %   width - box/cube width (of center box, proxy surface at
 %                1.5*width)
 %   opts - options structure
-%       opts.eps - tolerance (default 1e-13)
+%       opts.rank_or_tol - tolerance or maximum rank (default 1e-13)
+%       opts.eps - alias for rank_or_tol for backward compat (default 1e-13)
 %       opts.nsrc - number of sources to use in test (default 200)  
 %
 % Output:
-%   npxy - number of points on side of box to be passed
-%       to proxy surface routine
+%   npxy - number of points on perimeter of box to be sent to proxy routine
 %
 
   nsrc = 200;
-  eps = 1e-13;
+  rank_or_tol = 1e-13;
 
   if nargin < 3
     opts = [];
   end
 
   if isfield(opts,'eps')
-    eps = opts.eps;
+    rank_or_tol = opts.eps;
+  end
+  if isfield(opts,'rank_or_tol')
+    rank_or_tol = opts.rank_or_tol;
   end
   if isfield(opts,'nsrc')
     nsrc = opts.nsrc;
@@ -50,7 +53,7 @@
 
     mat = kern(srcinfo,targinfo);
 
-    [sk,~] = id(mat,eps);
+    [sk,~] = id(mat,rank_or_tol);
 
     if length(sk) < min(nsrc,npxy)
       npxy = (floor((length(sk)-1)/4+0.1)+1)*4;

chunkie/chunkerflam.m

@@ -0,0 +1,196 @@
+function [F] = chunkerflam(chnkr,kern,dval,opts)
+%CHUNKERFLAM build the requested FLAM compressed representation 
+% (e.g. a recursive skeletonization factorization) of the system matrix 
+% for given kernel and chunker description of boundary. This routine
+% has the same quadrature options as CHUNKERMAT.
+%
+% Syntax: sysmat = chunkerflam(chnkr,kern,opts)
+%
+% Input:
+%   chnkr - chunker object describing boundary
+%   kern  - kernel function. By default, this should be a function handle
+%           accepting input of the form kern(srcinfo,targinfo), where 
+%           srcinfo and targinfo are in the ptinfo struct format, i.e.
+%                ptinfo.r - positions (2,:) array
+%                ptinfo.d - first derivative in underlying
+%                     parameterization (2,:)
+%                ptinfo.d2 - second derivative in underlying
+%                     parameterization (2,:)
+%   dval - (default 0.0) float or float array. Let A be the matrix 
+%           corresponding to on-curve convolution with the provided kernel. 
+%           If a scalar is provided, the system matrix is 
+%                   A + dval*eye(size(A))  
+%           If a vector is provided, it should be length size(A,1). The
+%           system matrix is then
+%                   A + diag(dval)
+%
+% Optional input:
+%   opts  - options structure. available options (default settings)
+%           opts.flamtype = string ('rskelf'), type of compressed 
+%                           representation to compute. Available:
+%                               
+%                               The recursive skeletonization routines
+%                               should be sufficient curves which are not
+%                               approximately space filling.
+%
+%                               - 'rskelf', recursive skeletonization
+%                               factorization. Can be immediately used 
+%                               for system solves (via RSKELF_SV) and 
+%                               determinants (via RSKELF_LOGDET)
+%                               - 'rskel', recursive skeletonization
+%                               compression. Can be immediately used for
+%                               matvecs or embedded in a sparse matrix.
+%                               Not recommended unless you have a
+%                               compelling reason.
+%           opts.useproxy = boolean (true), use a proxy function to speed
+%                           up the FLAM compression. It may be desirable to
+%                           turn this off if there appear to be precision
+%                           issues.
+%
+%           opts.quad = string ('ggqlog'), specify quadrature routine to 
+%                       use. 
+%
+%                       - 'ggqlog' uses a generalized Gaussian quadrature 
+%                       designed for logarithmically singular kernels and 
+%                       smooth kernels with removable singularities
+%                       - 'native' selects standard scaled Gauss-Legendre 
+%                       quadrature for native functions
+%           opts.l2scale = boolean (false), if true scale rows by 
+%                           sqrt(whts) and columns by 1/sqrt(whts)
+%           opts.occ = integer "occupancy" parameter (200) determines
+%                   how many sources are in any leaf of the tree used to
+%                   sort points. Determines the smaller dimensions of the
+%                   first set of compressions. Sent to FLAM
+%           opts.rank_or_tol = integer or float "rank_or_tol" 
+%                   parameter (1e-14). Lower precision increases speed.
+%                   Sent to FLAM
+%
+% Output:
+%   F - the requested FLAM compressed representation of the 
+%           system matrix
+%
+% Examples:
+%   F = chunkermat(chnkr,kern,dval); % standard options
+%   sysmat = chunkermat(chnkr,kern,dval,opts);
+%
+
+F = [];
+
+if length(chnkr) > 1
+    chnkr = merge(chnkr);
+end
+
+if nargin < 3
+    dval = 0.0;
+end
+
+if nargin < 4
+    opts = [];
+end
+
+quad = 'ggqlog';
+flamtype = 'rskelf';
+useproxy = true;
+occ = 200;
+rank_or_tol = 1e-14;
+
+
+if or(chnkr.nch < 1,chnkr.k < 1)
+    warning('empty chunker, doing nothing')
+    sp = [];
+    return
+end
+
+% determine operator dimensions using first two points
+
+srcinfo = []; targinfo = [];
+srcinfo.r = chnkr.r(:,1); srcinfo.d = chnkr.d(:,1); 
+srcinfo.d2 = chnkr.d2(:,1);
+i2 = min(2,chnkr.npt);
+targinfo.r = chnkr.r(:,i2); targinfo.d = chnkr.d(:,i2); 
+targinfo.d2 = chnkr.d2(:,i2);
+
+ftemp = kern(srcinfo,targinfo);
+opdims = size(ftemp);
+
+if (length(dval) == 1)
+    dval = dval*ones(opdims(1)*chnkr.npt,1);
+end
+
+if (length(dval) ~= opdims(1)*chnkr.npt)
+    warning('provided dval array is length %d. must be scalar or length %d',...
+        length(dval),opdims(1)*chnkr.npt);
+    return
+end
+
+% get opts from struct if available
+
+if isfield(opts,'quad')
+    quad = opts.quad;
+end
+if isfield(opts,'l2scale')
+    l2scale = opts.l2scale;
+end
+
+if isfield(opts,'occ')
+    occ = opts.occ;
+end
+if isfield(opts,'rank_or_tol')
+    rank_or_tol = opts.rank_or_tol;
+end
+
+% check if chosen FLAM routine implemented before doing real work
+
+if ~ (strcmpi(flamtype,'rskelf') || strcmpi(flamtype,'rskel'))
+    warning('selected flamtype %s not available, doing nothing',flamtype)
+    return
+end
+
+% get nonsmooth quadrature
+
+if strcmpi(quad,'ggqlog')
+
+    type = 'log';
+    sp = chnk.quadggq.buildmattd(chnkr,kern,opdims,type);
+
+elseif strcmpi(quad,'native')
+
+    sp = sparse(chnkr.npt,chnkr.npt);
+
+else
+    warning('specified quadrature method not available');
+    return;
+end
+
+sp = sp + spdiags(dval,0,chnkr.npt,chnkr.npt);
+
+% prep and call flam
+
+wts = weights(chnkr);
+xflam = chnkr.r(:,:);
+
+width = max(max(chnkr)-min(chnkr));
+optsnpxy = []; optsnpxy.rank_or_tol = rank_or_tol;
+optsnpxy.nsrc = occ;
+
+npxy = chnk.flam.nproxy_square(kern,width,optsnpxy);
+
+matfun = @(i,j) chnk.flam.kernbyindex(i,j,chnkr,wts,kern,opdims,sp);
+[pr,ptau,pw,pin] = chnk.flam.proxy_square_pts(npxy);
+
+if strcmpi(flamtype,'rskelf')
+    ifaddtrans = true;
+    pxyfun = @(x,slf,nbr,l,ctr) chnk.flam.proxyfun(slf,nbr,l,ctr,chnkr,wts, ...
+        kern,opdims,pr,ptau,pw,pin,ifaddtrans);
+    F = rskelf(matfun,xflam,occ,rank_or_tol,pxyfun);
+end
+
+if strcmpi(flamtype,'rskel') 
+    pxyfunr = @(rc,rx,cx,slf,nbr,l,ctr) chnk.flam.proxyfunr(rc,rx,slf,nbr,l, ...
+        ctr,chnkr,wts,kern,opdims,pr,ptau,pw,pin);
+    F = rskel(matfun,xflam,xflam,occ,rank_or_tol,pxyfunr);
+end
+
+
+
+end

chunkie/chunkermat.m

@@ -28,7 +28,6 @@
 %                       - 'native' selects standard scaled Gauss-Legendre 
 %                       quadrature for native functions
 %
-%           opts.quadorder = integer (chnkr.k), desired quadrature order.
 %           opts.nonsmoothonly = boolean (false), if true, only compute the
 %                         entries for which a special quadrature is used
 %                         (e.g. self and neighbor interactoins) and return
@@ -53,7 +52,6 @@
     opts = [];
 end
 
-quadorder = chnkr.k;
 quad = 'ggqlog';
 nonsmoothonly = false;
 l2scale = false;
@@ -76,9 +74,6 @@
 
 % get opts from struct if available
 
-if isfield(opts,'quadorder')
-    quadorder = opts.quadorder;
-end
 if isfield(opts,'quad')
     quad = opts.quad;
 end

devtools/test/flamutilitiesTest.m

@@ -103,11 +103,24 @@
 ifaddtrans = true;
 pxyfun = @(x,slf,nbr,l,ctr) chnk.flam.proxyfun(slf,nbr,l,ctr,chnkr,wts, ...
     fkern,opdims,pr,ptau,pw,pin,ifaddtrans);
+
+
 start = tic; F = rskelf(matfun,xflam,200,1e-14,pxyfun); t1 = toc(start);
+F = chunkerflam(chnkr,fkern,-0.5);
+
+fprintf('%5.2e s : time for flam rskelf compress\n',t1)
+
+pxyfunr = @(rc,rx,cx,slf,nbr,l,ctr) chnk.flam.proxyfunr(rc,rx,slf,nbr,l, ...
+        ctr,chnkr,wts,fkern,opdims,pr,ptau,pw,pin);
+
+opts = [];
+start = tic; F2 = rskel(matfun,xflam,xflam,200,1e-14,pxyfunr,opts); t1 = toc(start);
+
+fprintf('%5.2e s : time for flam rskel compress\n',t1)
 
 afun = @(x) rskelf_mv(F,x);
 
-fprintf('%5.2e s : time for flam rskelf compress\n',t1)
+
 
 err2 = norm(sys2-sys,'fro')/norm(sys,'fro');
 fprintf('%5.2e   : fro error of build \n',err2)