switch to edgesendverts description of connectivity.

fastalgorithms · Mar 29, 2024 · cb5e2e4 · cb5e2e4
1 parent d7432c6
commit cb5e2e4
Show file tree

Hide file tree

Showing 7 changed files with 74 additions and 160 deletions.
diff --git a/chunkie/@chunkgraph/chunkgraph.m b/chunkie/@chunkgraph/chunkgraph.m
@@ -15,10 +15,11 @@
 % 
     properties(SetAccess=public)
         verts
-        edge2verts
+        edgesendverts
         echnks
         regions
         vstruc
+        v2emat
     end
 
     properties(SetAccess=public)
@@ -33,15 +34,31 @@
     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
@@ -78,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);
@@ -98,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');
@@ -119,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 = {};
@@ -208,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'};
@@ -275,7 +289,11 @@
             sourceinfo.d2= d2s;
             sourceinfo.w = ws;
         end
-    end   
+
+        % defined in other files 
+        spmat = build_v2emat(obj)
+    end
+
     methods(Static)
 
     end

diff --git a/chunkie/@chunkgraph/findregions.m b/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

diff --git a/chunkie/@chunkgraph/vertextract.m b/chunkie/@chunkgraph/vertextract.m
@@ -26,9 +26,9 @@
 % author: Jeremy Hoskins
 
 % extract the indices of the edges which terminate at ivert.
-ieplus = find(cgrph.edge2verts(:,ivert) ==  1);
+ieplus = find(cgrph.edgesendverts(2,:) ==  ivert);
 % extract the indices of the edges which begin at ivert.
-ieminus  = find(cgrph.edge2verts(:,ivert) == -1);
+ieminus  = find(cgrph.edgesendverts(1,:) == ivert);
 
 % for each incoming edge, get the tangent vector near the end (at the 
 % last discretization node)
@@ -50,8 +50,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
diff --git a/chunkie/mergeregions.m b/chunkie/mergeregions.m
@@ -3,7 +3,7 @@
     rgnout = {};
 
     e2 = rgn2{1}{1}(1);
-    v2 = find(cgrph.edge2verts(:,abs(e2))==1);
+    v2 = cgrph.edgesendverts(2,abs(e2));
     v2 = cgrph.verts(:,v2);
 
     irgn = 0;
@@ -27,7 +27,7 @@
 
     if (irgn == 0)
         e1 = rgn1{1}{1}(1);
-        v1 = find(cgrph.edge2verts(:,abs(e1))==1);
+        v1 = cgrph.edgesendverts(2,abs(e1));
         v1 = cgrph.verts(:,v1);
 
         irgn = 0;

diff --git a/chunkie/regioninside.m b/chunkie/regioninside.m
@@ -1,7 +1,7 @@
 function [isinside] = regioninside(cgrph,rgn1,rgn2)
 
     e2 = rgn2{1}{1}(1);
-    v2 = find(cgrph.edge2verts(:,abs(e2))==1);
+    v2 = cgrph.edgesendverts(2,abs(e2));
     v2 = cgrph.verts(:,v2);
     irgn = 0;
     for ii=2:numel(rgn1)

diff --git a/devtools/test/chunkgrphconstructTest.m b/devtools/test/chunkgrphconstructTest.m
@@ -45,60 +45,17 @@
 verts = exp(1i*2*pi*(0:4)/5);
 verts = [real(verts);imag(verts)];
 
+% legacy connectivity info
 edge2verts = [-1, 1, 0, 0, 0; ...
                0,-1, 1, 0, 0; ...
                0, 0,-1, 1, 0; ...
                0, 0, 0,-1, 1; ...
                1, 0, 0, 0,-1];
 edge2verts = sparse(edge2verts);
-amp = 0.5;
-frq = 6;
-fchnks    = {};
-for icurve = 1:size(edge2verts,1)
-    fchnks{icurve} = @(t) circulararc(t,cpars{1});
-    fchnks{icurve} = @(t) sinearc(t,amp,frq);
-end
-%fchnks = {};
-
-prefs      = [];
-% prefs.chsmall = 1d-4;
-[cgrph] = chunkgraph(verts,edge2verts,fchnks,prefs);
-
-vstruc = procverts(cgrph);
-rgns = findregions(cgrph);
-cgrph = balance(cgrph);
-
-
-
-ncurve = 1;
-
-% set angles for top and bottom curve
-theta = [pi/2.2 pi/2.4];
-
-% set parametrization start and end points for chunkie objects
-tab = [-theta(1)/2 -theta(2)/2; theta(1)/2 theta(2)/2];
-
-% set parameters for curve parametrization
-vert = [a b; 0 0];
-cpars = cell(1,ncurve);
-cpars{1}.v0 = vert(:,1); cpars{1}.v1 = vert(:,2);
-cpars{1}.theta = theta(1); cpars{1}.ifconvex = 0; cpars{1}.islocal = -1;
-
-cpars{2}.v0 = vert(:,1); cpars{2}.v1 = vert(:,2);
-cpars{2}.theta = theta(2); cpars{2}.ifconvex = 2; cpars{2}.islocal = -1;
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+% connectivity info
+edgesendverts = [1:5; [2:5 1]];
 
-verts = exp(1i*2*pi*(0:4)/5);
-verts = [real(verts);imag(verts)];
-
-edge2verts = [-1, 1, 0, 0, 0; ...
-               0,-1, 1, 0, 0; ...
-               0, 0,-1, 1, 0; ...
-               0, 0, 0,-1, 1; ...
-               1, 0, 0, 0,-1];
-edge2verts = sparse(edge2verts);
 amp = 0.5;
 frq = 6;
 fchnks    = {};
@@ -109,86 +66,15 @@
 %fchnks = {};
 
 prefs      = [];
-% prefs.chsmall = 1d-4;
-[cgrph] = chunkgraph(verts,edge2verts,fchnks,prefs);
-
-vstruc = procverts(cgrph);
-rgns = findregions(cgrph);
-cgrph = balance(cgrph);
-
-
-
-zk = 1.0;
-fkern = @(s,t) -2*chnk.helm2d.kern(zk,s,t,'d');
-
-opts = [];
-opts.nonsmoothonly = false;
-opts.rcip = true;
-[sysmat] = chunkermat(cgrph,fkern,opts);
-sysmat = sysmat + eye(size(sysmat,2));
-
-% generate some targets...
-
-xs = -2:0.01:2;
-ys = -2:0.01:2;
-[X,Y] = meshgrid(xs,ys);
-targs = [X(:).';Y(:).'];
-
-srcinfo = [];
-srcinfo.sources = cgrph.r(:,:);
-w = weights(cgrph);
-n = normals(cgrph);
-
-% a quick hack to find the interior points
-
-srcinfo.dipstr = w(:).';
-srcinfo.dipvec = n(:,:); 
-eps = 1E-8;
-pg  = 0;
-pgt = 1;
-[U] = lfmm2d(eps,srcinfo,pg,targs,pgt);
-U = U.pottarg;
-inds = find(abs(U-2*pi)<pi/10);
-
-%%%%%%%%%%%%%%%%%%
-% generate the right hand side
-
-x0 = 1.3;
-y0 = 0.9;
-
-srcinfo = [];
-srcinfo.sources = cgrph.r(:,:);
-w = weights(cgrph);
-n = normals(cgrph);
-
-s = [];
-s.r = [x0;y0];
-t = [];
-t.r = cgrph.r(:,:);
-rhs = chnk.helm2d.kern(zk,s,t,'s');
-dens = sysmat\rhs;
-dens = -2*dens;
-
-srcinfo.dipstr = (w(:).*dens(:)).';
-srcinfo.dipvec = n(:,:); 
-fints = hfmm2d(eps,zk,srcinfo,pg,targs(:,inds),pgt);
-
-t.r = targs(:,inds);
-true_sol = chnk.helm2d.kern(zk,s,t,'s');
-
-usol = zeros(size(targs,2),1);
-uerr = usol;
-usol(inds) = fints.pottarg;
-uerr(inds) = usol(inds)-true_sol;
-usol = reshape(usol,size(X));
-uerr = reshape(uerr,size(X));
-h = pcolor(X,Y,log10(abs(uerr)));
-set(h,'EdgeColor','None'); hold on;
-plot(cgrph,'w-','LineWidth',2);
-caxis([-16,0])
-colorbar
 
+[cgrph1] = chunkgraph(verts,edge2verts,fchnks,prefs);
+[cgrph2] = chunkgraph(verts,edgesendverts,fchnks,prefs);
+
+cgrph1 = balance(cgrph1);
+cgrph2 = balance(cgrph2);
 
+assert(all(cgrph1.verts(:) == cgrph2.verts(:)))
+assert(all(all(cgrph1.v2emat == cgrph2.v2emat)))
 
 function [r,d,d2] = circulararc(t,cpars)
 %%circulararc

diff --git a/devtools/test/chunkgrphregionTest.m b/devtools/test/chunkgrphregionTest.m
@@ -7,6 +7,7 @@
 
 %clear all
 
+clearvars
 addpaths_loc();
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%