Merge branch 'master' of https://github.com/fastalgorithms/chunkie in…

…to docs

chunkie/@chunker/max.m

@@ -17,4 +17,4 @@
 
 % author: Travis Askham ([email protected])
 
-rmax = max(reshape(obj.r,obj.dim,obj.k*obj.nch),[],2);
+rmax = max(real(reshape(obj.r,obj.dim,obj.k*obj.nch)),[],2);

chunkie/@chunker/min.m

@@ -17,4 +17,4 @@
 
 % author: Travis Askham ([email protected])
 
-rmin = min(reshape(obj.r,obj.dim,obj.k*obj.nch),[],2);
+rmin = min(real(reshape(obj.r,obj.dim,obj.k*obj.nch)),[],2);

chunkie/@chunkgraph/build_v2emat.m

@@ -0,0 +1,42 @@
+function A = build_v2emat(obj)
+%V2EMAT get a "vertex to edge" matrix for the chunkgraph object
+% 
+% input: obj - a chunkgraph object 
+% 
+% output: A - a nedges x nverts sparse matrix 
+%  A(i,j) = 1 if edge i starts at vertex j
+%  A(i,j) = -1 if edge i ends at vertex j
+%  A(i,j) = 2 if edge i starts and ends at vertex j
+%
+
+nverts = size(obj.verts,2);
+nedges = size(obj.edgesendverts,2);
+
+
+ii = zeros(2*nedges,1);
+jj = zeros(2*nedges,1);
+vv = zeros(2*nedges,1);
+nf = 0;
+
+for i = 1:nedges
+    j1 = obj.edgesendverts(1,i);
+    j2 = obj.edgesendverts(2,i);
+    if j1 == j2
+        nf = nf + 1;
+        ii(nf) = i;
+        jj(nf) = j1;
+        vv(nf) = 2;
+    else
+        nf = nf + 1;
+        ii(nf) = i;
+        jj(nf) = j1;
+        vv(nf) = -1;
+        nf = nf + 1;
+        ii(nf) = i;
+        jj(nf) = j2;
+        vv(nf) = 1;
+    end
+end
+
+A = sparse(ii(1:nf),jj(1:nf),vv(1:nf),nedges,nverts);
+

chunkie/@chunkgraph/chunkgraph.m

@@ -15,10 +15,11 @@
 % 
     properties(SetAccess=public)
         verts
-        edge2verts
+        edgesendverts
         echnks
         regions
         vstruc
+        v2emat
     end
 
     properties(SetAccess=public)
@@ -33,16 +34,36 @@
     end
 
     methods
-        function obj = chunkgraph(verts,edge2verts,fchnks,cparams)
+        function obj = chunkgraph(verts,edgesendverts,fchnks,cparams)
             if (nargin == 0)
                 return
             end
             if (numel(verts)==0)
                 return
             end
-            obj.verts      = verts;
-            obj.edge2verts = edge2verts;
+            obj.verts = verts;
+
+            nverts = size(verts(:,:),2);
+            assert(nverts == size(edgesendverts,2),'edge specification not compatible with number of vertices');
+            if (size(edgesendverts,1) ~= 2)
+                nedge = size(edgesendverts,1);
+                edgevertends_new = zeros(2,nedge);
+                nedge = size(edgesendverts,1);
+                for i = 1:nedge
+                    edgevertends_new(1,i) = find(edgesendverts(i,:) == -1);
+                    edgevertends_new(2,i) = find(edgesendverts(i,:) == 1);
+                end
+                edgesendverts = edgevertends_new;
+                assert(all(edgesendverts(:) ~= 0),'edge specification had an error');
+            end 
+
+            obj.edgesendverts = edgesendverts;
+            obj.v2emat = build_v2emat(obj);
             obj.echnks     = chunker.empty;
+
+            if nargin < 3
+                fchnks = [];
+            end
 
             if (nargin < 4)
                 cploc = [];
@@ -74,15 +95,11 @@
             pref.nchmax = 10000;
             pref.k = 16;
 
-            if (size(verts,2) ~= size(edge2verts,2))
-                error('Incompatible vertex and edge sizes');
-            end
-
             echnks = chunker.empty();
-            for i=1:size(edge2verts,1)
+            for i=1:size(edgesendverts,2)
                 if (numel(fchnks)<i || isempty(fchnks{i}))
-                    i1 = find(edge2verts(i,:)==-1);
-                    i2 = find(edge2verts(i,:)==1);
+                    i1 = edgesendverts(1,i);
+                    i2 = edgesendverts(2,i);
                     v1 = verts(:,i1);
                     v2 = verts(:,i2);
                     fcurve = @(t) chnk.curves.linefunc(t,v1,v2);
@@ -94,8 +111,10 @@
                     [vs,~,~] =fchnks{i}([0,1]);
                     chnkr = chunkerfunc(fchnks{i},cploc,pref);
                     chnkr = sort(chnkr);
-                    vfin0 = verts(:,find(edge2verts(i,:)==-1));
-                    vfin1 = verts(:,find(edge2verts(i,:)== 1));
+                    i1 = edgesendverts(1,i);
+                    i2 = edgesendverts(2,i);
+                    vfin0 = verts(:,i1);
+                    vfin1 = verts(:,i2);
                     r0 = vs(:,1);
                     r1 = vfin0;
                     scale = norm(vfin1-vfin0,'fro')/norm(vs(:,2)-vs(:,1),'fro');
@@ -115,9 +134,8 @@
             obj.wts = weights(obj);
             %[regions] = findregions(obj);
             %obj.regions = regions;
-
-            adjmat = edge2verts'*edge2verts;
-            g = graph(adjmat);
+
+            g = graph(edgesendverts(1,:),edgesendverts(2,:));
             ccomp = conncomp(g);
 
             chnkcomp = {};
@@ -204,10 +222,10 @@
             validateattributes(val,classes,{})
             obj.verts = val;
         end
-        function obj = set.edge2verts(obj,val)
+        function obj = set.edgesendverts(obj,val)
             classes = {'numeric'};
             validateattributes(val,classes,{})
-            obj.edge2verts = val;
+            obj.edgesendverts = val;
         end    
         function obj = set.echnks(obj,val)
             classes = {'chunker'};
@@ -271,7 +289,11 @@
             sourceinfo.d2= d2s;
             sourceinfo.w = ws;
         end
-    end   
+
+        % defined in other files 
+        spmat = build_v2emat(obj)
+    end
+
     methods(Static)
 
     end

chunkie/@chunkgraph/findregions.m

@@ -27,18 +27,17 @@
 else
     vstruc = procverts(obj);
 end    
-nedge  = size(obj.edge2verts,1);
-e2v    = obj.edge2verts;
-
-% each edge belongs to two regions (going in opposite directions)
-edges = [1:nedge,-(1:nedge)];
+nedge  = size(obj.edgesendverts,2);
 
 % if the user has provided iverts, get the reduced edge list
 if (nargin >1)
-    e2vtmp = e2v(:,iverts);
-    [iinds,jinds] = find(e2vtmp ~= 0);
+    v2etmp = obj.v2emat(:,iverts);
+    [iinds,jinds] = find(v2etmp ~= 0);
     iinds = unique(iinds);
     edges = [iinds,-iinds];
+else
+    % each edge belongs to two regions (going in opposite directions)
+    edges = [1:nedge,-(1:nedge)];
 end    
 
 regions = {};
@@ -56,8 +55,14 @@
     estart = enum;
     ecycle = [enum];
 
-    iv0 = find(e2v(abs(enum),:)==-sign(enum));
-    ivc = find(e2v(abs(enum),:)== sign(enum));
+    if enum > 0
+        iv0 = obj.edgesendverts(1,abs(enum));
+        ivc = obj.edgesendverts(2,abs(enum));
+    else
+        iv0 = obj.edgesendverts(2,abs(enum));
+        ivc = obj.edgesendverts(1,abs(enum));
+    end
+
     ifdone = false;
 
     while (~ifdone)
@@ -72,7 +77,11 @@
             esign = vstruc{ivc}{2}(1);
         end
         enum = -esign*enext;
-        ivc = find(e2v(abs(enum),:)== sign(enum));
+        if enum > 0
+            ivc = obj.edgesendverts(2,abs(enum));
+        else
+            ivc = obj.edgesendverts(1,abs(enum));
+        end
         if (enum == estart)
             ifdone = true;
         else

chunkie/@chunkgraph/vertextract.m

@@ -25,10 +25,14 @@
 
 % author: Jeremy Hoskins
 
+irel = find(cgrph.v2emat(:,ivert) ~= 0);
 % extract the indices of the edges which terminate at ivert.
-ieplus = find(cgrph.edge2verts(:,ivert) ==  1);
+ieplus = find(cgrph.edgesendverts(2,irel) ==  ivert);
+ieplus = irel(ieplus); ieplus = ieplus(:).';
 % extract the indices of the edges which begin at ivert.
-ieminus  = find(cgrph.edge2verts(:,ivert) == -1);
+ieminus  = find(cgrph.edgesendverts(1,irel) == ivert);
+ieminus = irel(ieminus); ieminus = ieminus(:).';
+
 
 % for each incoming edge, get the tangent vector near the end (at the 
 % last discretization node)
@@ -50,8 +54,8 @@
 [angs,isrtededges] = sort(angs);
 
 % compute inds and isgn using isrtedges
-inds = [ieplus',ieminus'];
-isgn = [ones(size(ieplus')),-ones(size(ieminus'))];
+inds = [ieplus,ieminus];
+isgn = [ones(size(ieplus)),-ones(size(ieminus))];
 inds = inds(isrtededges);
 isgn = isgn(isrtededges);
 end

chunkie/chunkerflam.m

@@ -244,7 +244,7 @@
     mmax = max([mmax,max(chnkr)],[],2);
     mmin = min([mmin,min(chnkr)],[],2);
 end
-
+xflam = real(xflam);
 width = max(mmax-mmin);
 
 chnkrtotal = merge(chnkrs);

chunkie/chunkerfunc.m

@@ -19,6 +19,8 @@
 %	cparams - curve parameters structure (defaults)
 %       cparams.ta = left end of t interval (0)
 %       cparams.tb = right end of t interval (2*pi)
+%       cparams.tsplits = set of initial break points for discretization in
+%                 parameter space (should be in [ta,tb])
 %       cparams.ifclosed = flag determining if the curve
 %           is to be interpreted as a closed curve (true)
 %       cparams.chsmall = max size of end intervals if
@@ -139,19 +141,42 @@
 
 ab = zeros(2,nchmax);
 adjs = zeros(2,nchmax);
-ab(1,1)=ta;
-ab(2,1)=tb;
-nch=1;
+
+if (isfield(cparams,'tsplits'))
+    tsplits = cparams.tsplits;
+    tsplits = [tsplits(:); ta; tb];
+else
+    tsplits = [ta;tb];
+end
+
+tsplits = sort(unique(tsplits),'ascend');
+lab = length(tsplits);
+if (lab-1 > nchmax)
+    error(['nchmax exceeded in chunkerfunc on initial splits.\n ',...
+        'try increasing nchmax in preference struct']);
+end
+if (any(tsplits > tb) || any(tsplits < ta))
+    error(['tsplits outside of interval of definition.\n', ...
+          'check definition of splits, ta and tb']);
+end
+
+ab(1,1:(lab-1)) = tsplits(1:end-1);
+ab(2,1:(lab-1)) = tsplits(2:end);
+
+nch=lab-1;
+adjs(1,1:nch) = 0:(nch-1);
+adjs(2,1:nch) = 2:(nch+1);
+
 if ifclosed
-    adjs(1,1)=1;
-    adjs(2,1)=1;
+    adjs(1,1)=nch;
+    adjs(2,nch)=1;
 else
     adjs(1,1)=-1;
-    adjs(2,1)=-1;
+    adjs(2,nch)=-1;
 end
 nchnew=nch;
 
-maxiter_res=nchmax;
+maxiter_res=nchmax-nch;
 
 rad_curr = 0;
 for ijk = 1:maxiter_res

chunkie/chunkermat.m

@@ -485,11 +485,12 @@
 
         [Pbc,PWbc,starL,circL,starS,circS,ilist,starL1,circL1] = ...
             chnk.rcip.setup(ngl,ndim,nedge,isstart);
-
-        % this might need to be fixed in triple junction case
+        optsrcip = opts;
+        optsrcip.nonsmoothonly = false;
+
         R = chnk.rcip.Rcompchunk(chnkrs,iedgechunks,kern,ndim, ...
             Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,starL1,circL1,... 
-            sbclmat,sbcrmat,lvmat,rvmat,u,opts);
+            sbclmat,sbcrmat,lvmat,rvmat,u,optsrcip);
 
         sysmat_tmp = inv(R) - eye(2*ngl*nedge*ndim);
         if (~nonsmoothonly)