Merge pull request #75 from fastalgorithms/feature/issue-remove-chnkr…

…-dot-h

Feature/issue remove chnkr dot h

chunkie/+chnk/+quadba/buildmattd.m

@@ -11,7 +11,6 @@
 adj = chnkr.adj;
 d = chnkr.d;
 d2 = chnkr.d2;
-h = chnkr.h;
 
 [~,~,u] = lege.exps(k);
 

chunkie/+chnk/adapgausswts.m

@@ -8,15 +8,14 @@
 % WARNING: this routine currently assumes that the kernel function is
 %   translation invariant to recenter (improves stability). 
 %
-% Syntax: [mat,maxrecs,numints,iers] = adapgausswts(r,d,d2,h,ct,bw,j, ...
+% Syntax: [mat,maxrecs,numints,iers] = adapgausswts(r,d,d2,ct,bw,j, ...
 %   rt,dt,d2t,kern,opdims,t,w,opts)
 %
 % Input:
 %   r - chnkr nodes
 %   d - chnkr derivatives at nodes
 %   n - chnkr normals at nodes
 %   d2 - chnkr 2nd derivatives at nodes
-%   h - lengths of chunks in parameter space
 %   ct - Legendre nodes at order of chunker
 %   bw - barycentric interpolation weights for Legendre nodes at order of
 %   chunker

chunkie/+chnk/chunk_nearparam.m

@@ -24,6 +24,7 @@
 
 maxnewt = 15;
 thresh0 = 1.0d-14;
+iextra = 3;
 
 if nargin < 4 || isempty(t) || nargin < 5 || isempty(u)
     [t,~,u] = lege.exps(chnkr.k);
@@ -40,12 +41,6 @@
     thresh0 = opts.thresh;
 end
 
-dist = Inf;
-
-if nargin < 3
-    ich = 1:nch;
-end
-
 dim = size(pts,1);
 npts = numel(pts)/dim;
 pts = reshape(pts,dim,npts);
@@ -92,24 +87,28 @@
 
     rdiff0 = r0(:) - ref;
     dprime = rdiff0.'*d0(:);
-    dprime2 = d0(:).'*d0(:) + d0(:).'*d20(:);
+    dprime2 = (d0(:).'*d0(:) + rdiff0(:).'*d20(:));
 
     newtsuccess = false;
+    ifend = 0;
 
     for iii = 1:maxnewt
 
         % Newton step
 
         deltat = - dprime/dprime2;
 
-        t0 = t0+deltat;
+        t1 = t0 + deltat;
 
-        if (t0 > 1.0)
-            t0 = 1;
+        if (t1 > 1.0)
+            t1 = 1;
         end
-        if (t0 < -1.0)
-            t0 = -1;
+        if (t1 < -1.0)
+            t1 = -1;
         end
+
+        deltat = t1-t0;
+        t0 = t1;
 
         all0 = (lege.exev(t0,cfs));
         r0 = all0(1:dim);
@@ -118,9 +117,13 @@
 
         rdiff0 = r0(:)-ref(:);
         dprime = rdiff0.'*d0(:);
-        dprime2 = d0(:).'*d0(:) + d0(:).'*d20(:);
+        dprime2 = d0(:).'*d0(:) + rdiff0(:).'*d20(:);
+
+        if min(abs(dprime),abs(deltat)) < thresh
+            ifend = ifend+1;
+        end
 
-        if abs(dprime) < thresh
+        if (ifend >= iextra)
             newtsuccess = true;
             break;
         end
@@ -196,7 +199,7 @@
                 d20 = d21;
                 rdiff0 = rdiff1;
                 dprime = rdiff0.'*d0(:);
-                dprime2 = d0(:).'*d0(:) + d0(:).'*d20(:);
+                dprime2 = d0(:).'*d0(:);
                 dist0 = dist1;
                 lam = lam/3;
                 if abs(dprime) < thresh

chunkie/+chnk/funcuni.m

@@ -104,16 +104,18 @@
 for i = 1:nch
     a=ab(1,i);
     b=ab(2,i);
+    h = (b-a)/2;
 
     ts = a + (b-a)*(xs+1)/2;
     [out{:}] = fcurve(ts);
     chnkr.r(:,:,i) = reshape(out{1},dim,k);
-    chnkr.d(:,:,i) = reshape(out{2},dim,k);
-    chnkr.d2(:,:,i) = reshape(out{3},dim,k);
-    chnkr.h(i) = (b-a)/2;
+    chnkr.d(:,:,i) = reshape(out{2},dim,k)*h;
+    chnkr.d2(:,:,i) = reshape(out{3},dim,k)*h*h;
 end
 
 chnkr.adj = adjs(:,1:nch);
+chnkr.wtsstor(:,1:nch) = weights(chnkr);
+chnkr.nstor(:,:,1:nch) = normals(chnkr);
 
 end
 

chunkie/@chunker/arclengthdens.m

@@ -1,10 +1,9 @@
 function ds = arclengthdens(chnkr)
 %ARCLENGTHDENS arc length density on chunks
 %
-% warning: this takes the length of the ith chunk in parameter space to 
-%   be 2*chnkr.h(i) as opposed to 2. Thus the smooth integration rule
+%   The smooth integration rule
 %   on the ith chunk is 
-%       ds(:,i).*w*chnkr.h(i)
+%       ds(:,i).*w
 %   where w are the standard Legendre weights of appropriate order and 
 %   ds is the output of this routine
 %

chunkie/@chunker/arclengthfun.m

@@ -13,10 +13,9 @@
 istart = 1;
 
 A = lege.intmat(chnkr.k);
-[~,w] = lege.exps(chnkr.k);
 ds = arclengthdens(chnkr);
-chunklens = sum((w(:)*(chnkr.h(:).')).*ds,1);
-s = (A*ds).*(chnkr.h(:).');
+chunklens = chunklen(chnkr);
+s = A*ds;
 
 for i = 1:ncomp
     nch = nchs(i);

chunkie/@chunker/move.m

@@ -19,4 +19,4 @@
 
 obj.wts = weights(obj);
 
-end
+end

chunkie/@chunker/nearest.m

@@ -40,6 +40,7 @@
 maxnewt = 15;
 thresh0 = 1.0d-9;
 iextra = 3;
+nch = chnkr.nch;
 
 if nargin < 5 || isempty(u)
     [~,~,u] = lege.exps(chnkr.k);
@@ -69,124 +70,15 @@
     % grab chunk data
     ii = ich(i);
     r = chnkr.rstor(:,:,ii);
-    d = chnkr.dstor(:,:,ii);
-    d2 = chnkr.d2stor(:,:,ii);
-    h = chnkr.hstor(ii);
-    rc = u*(r.');
-    dc = u*(d.');
-    d2c = u*(d2.');
 
-    % starting guess on over sampled grid
-    rover = (ainterpover*(r.')).';
-    rdist0 = sqrt(sum((rover-ref(:)).^2,1));
-    [disti,ind] = min(rdist0);
-    tt = xover(ind);
-
-    % check the endpoints
-    ts = [-1;1];
-    rs = (lege.exev(ts,rc)).';
-    rsdist = sqrt(sum((rs-ref).^2,1));
-    [dists,inds] = min(rsdist);
-    if dists < disti
-        % continue if endpoints closer
-        ri = rs(:,inds);
-        if dists < dist
-            dist = dists;
-            ds = (lege.exev(ts,dc)).';
-            d2s = (lege.exev(ts,d2c)).';
-            di = ds(:,inds);
-            d2i = d2s(:,inds);
-            rn = ri; dn = di; d2n = d2i;
-            tn = ts(inds);
-            ichn = ii;
-        end
-        continue;
-    end
-
-    % run levenberg
-
-    t0 = tt;
-
-    ifend = 0;
-
-    r0 = (lege.exev(t0,rc)).';
-
-    rdiff0 = ref(:) - r0(:);
-
-    drc = lege.derpol(rc);
-    %d2rc = lege.derpol(drc);
-    thresh = thresh0;
-
-    for iii = 1:maxnewt
-        %d0 = (lege.exev(t0,dc)).'*h;
-        d0 = (lege.exev(t0,drc)).'; % actual derivative of chunk
-%        d20 = (lege.exev(t0,d2rc)).'; % actual 2nd derivative of chunk
-
-        % newton info
-%         
-%         dprime = rdiff0.'*d0(:);
-%         dprime2 = (d0(:).'*d0(:) + d0(:).'*d20(:));
-
-        % levenberg instead
-
-        deltat = d0(:)\(rdiff0(:));
-        t1 = t0+deltat;
-
-        if (t1 > 1.0)
-            t1 = (1.0+t0)/2;
-        end
-        if (t1 < -1.0)
-            t1 = (-1.0+t0)/2;
-        end
-
-        r1 = (lege.exev(t1,rc)).';
-        rdiff1 = ref(:)-r1(:);
-
-        err = min(abs(deltat),abs(sum(d0(:).*rdiff0(:))));
-        if err < thresh
-            ifend = ifend+1;
-        end
-
-        %sum(rdiff0(:).*d0(:))
-        %abs(deltat)
-
-        t0 = t1;
-        rdiff0 = rdiff1;
-
-%         if (err < thresh)
-%             ifend = ifend + 1;
-%         end
-% 
-        %t0 = t1;
-
-        if (ifend >= iextra)
-            break;
-        end
-    end
-
-    % done levenberg
-
-    ri = (lege.exev(t0,rc)).';
-    disti = sqrt(sum((ri(:)-ref(:)).^2,1));
-    if (ifend < iextra)
-        warning('newton failed in nearest');
-%          deltat
-%          disti
-%          jac
-%          t0
-%          rhs
-%          ifend
-%          sum(jac.*rhs)
-    else
-%        fprintf('success\n')
-    end
-
+    [ti,ri,di,d2i,disti] = chnk.chunk_nearparam(r,ref,opts,chnkr.tstor,u);
+    disti = sqrt(disti);
     if disti < dist
         dist = disti;
         rn = ri;
-        dn = (lege.exev(t1,dc)).'*h;
-        d2n = (lege.exev(t1,d2c)).'*h*h;
-        tn = t1;
+        dn = di;
+        d2n = d2i;
+        tn = ti;
         ichn = ii;
     end
 end

chunkie/@chunker/whts.m

@@ -1,12 +1,8 @@
-function wchnk = whts(chnkr)
+function wts = whts(chnkr)
 
 
 warning('whts is deprecated and will be removed, use weights instead');
 
-  k = chnkr.k;
-  nch = chnkr.nch;
-  [~,w] = lege.exps(k);
-  wchnk = reshape(sqrt(sum((chnkr.d).^2,1)),k,nch);
-  wchnk = wchnk.*bsxfun(@times,w(:),(chnkr.h(:)).');
+    wts = weights(chnkr);
 
 end

chunkie/chunkerinterior.m

@@ -99,7 +99,6 @@
     chnkr.d2(:,:,istart:iend) = fliplr(chnkr.d2(:,:,1:nch));
     chnkr.d2(1,:,istart:iend) = chnkr.d2(1,:,istart:iend);
 
-    chnkr.h(istart:iend) = chnkr.h(1:nch);
     chnkr.wts = weights(chnkr);
     chnkr.n = normals(chnkr);
     chnkr.adj(1,:) = 0:chnkr.nch-1;

chunkie/chunkerkernevalmat.m

@@ -273,28 +273,13 @@
     d = chnkr.d;
     n = chnkr.n;
     d2 = chnkr.d2;
-    h = chnkr.h;
 
     for i = 1:nch
         jmat = 1 + (i-1)*k*opdims(2);
         jmatend = i*k*opdims(2);
 
-        mat(:,jmat:jmatend) =  chnk.adapgausswts(r,d,n,d2,h,ct,bw,i,targs, ...
-                    targd,targd2,kern,opdims,t,w,opts);
-
-        js1 = jmat:jmatend;
-        js1 = repmat( (js1(:)).',1,opdims(1)*numel(ji));
-
-        indji = (ji-1)*opdims(1);
-        indji = repmat( (indji(:)).', opdims(1),1) + ( (1:opdims(1)).');
-        indji = indji(:);
-        indji = repmat(indji,1,opdims(2)*k);
-
-        iend = istart+numel(mat1)-1;
-        is(istart:iend) = indji(:);
-        js(istart:iend) = js1(:);
-        vs(istart:iend) = mat1(:);
-        istart = iend+1;
+        mat(:,jmat:jmatend) =  chnk.adapgausswts(r,d,n,d2,ct,bw,i,targs, ...
+                    targd,targn,targd2,kern,opdims,t,w,opts);
     end
 
 else

chunkie/chunkermat.m

@@ -556,14 +556,13 @@
 d = chnkr.d;
 n = chnkr.n;
 d2 = chnkr.d2;
-h = chnkr.h;
 
 for i = 1:nch
     jmat = 1 + (i-1)*k*opdims(2);
     jmatend = i*k*opdims(2);
 
     [ji] = find(flag(:,i));
-    mat1 =  chnk.adapgausswts(r,d,n,d2,h,ct,bw,i,targs(:,ji), ...
+    mat1 =  chnk.adapgausswts(r,d,n,d2,ct,bw,i,targs(:,ji), ...
                 targd(:,ji),targn(:,ji),targd2(:,ji),kernev,opdims,t,w,opts);
 
     js1 = jmat:jmatend;