Tetra: unfolding now appears to be working correctly after several failed attempts

Research/Edge Intersection/Tetrahedra/ALGO_ProjectionMatrix3D_v2_20211222.m

+function M = ALGO_ProjectionMatrix3D_v2_20211222(Na,Ta,Nb,Tb)
+% PROJECTIONMATRIX3D calculates the projection matrix from tesselation 1 to
+% tesselation 2
+%   M=ProjectionMatrix3D(N1,T1,N2,T2) constructs the projection matrix M
+%   from the nodes N1 to the nodes N2, with the tesselations T1 and T2,
+%   respectively.
+
+%===Advancing front algorithm (WIP)===%
+na=size(Ta,1); nb=size(Tb,1);
+cand=[1;1];
+listb=1;
+flagb=zeros(nb+1,1); flagb([1,end])=1;
+while ~isempty(listb)
+    Tetb=listb(1); listb=listb(2:end);
+    lista=cand(1,cand(2,:)==Tetb);
+    lista=unique(lista);
+    flaga=zeros(na+1,1);
+    flaga(end)=1; flaga(lista)=1;
+    while ~isempty(lista)
+        Teta=lista(1); lista=lista(2:end);
+        [P,n]=ALGO_simplices_IntersectSimplices_v1_20211216(Nb(:,Tb(Tetb,1:4)),Na(:,Ta(Teta,1:4)));
+        if ~isempty(P)
+            nhbrsa=Ta(Teta,5:8);
+            lista=[lista, nhbrsa(flaga(nhbrsa)==0)];
+            flaga(nhbrsa)=1;
+            nhbrsb=Tb(Tetb,5:8);
+            nhbrsb=nhbrsb(n==1);  % nhbrs of Tetb which intersect Teta
+            listb=[listb, nhbrsb(flagb(nhbrsb)==0)];
+            flagb(nhbrsb)=1;
+            nhbrsb=[Teta*ones(size(nhbrsb));nhbrsb];
+            cand=[cand, nhbrsb];
+        end
+    end
+end
+end
\ No newline at end of file

Research/Edge Intersection/Tetrahedra/EX_3DHecht_v1_20201103.m

@@ -2,12 +2,13 @@
 th=pi/10; disp=0.1*ones(3,1);
 [N,T,N2]=MESH_d20_20201103(th,th,th,disp);
 
-figure(1)
+figure(2)
 clf
 PLOT_TetraMesh_20201103(N,T,N2,T);
-ALGO_ProjectionMatrix3D_v1_20201104(N,T,N2,T);
+InterfaceMatrix3d(N,T,N2,T);
 
-figure(2)
+figure(1)
 clf
 PLOT_TetraMesh_20201103(N,T,N2,T);
-InterfaceMatrix3d(N,T,N2,T);
\ No newline at end of file
+% ALGO_ProjectionMatrix3D_v1_20201104(N,T,N2,T);
+ALGO_ProjectionMatrix3D_v2_20211222(N,T,N2,T);
\ No newline at end of file

Research/Edge Intersection/Tetrahedra/EX_3DHecht_v2_20211222.m

+%===EX_3DHecht_v1_20201103===%
+th=pi/10; disp=0.1*ones(3,1);
+[N,T,N2]=MESH_d20_20201103(th,th,th,disp);
+U=N2(:,T(1,1:4)); V=N(:,T(1,1:4));
+
+figure(1)
+clf
+PLOT_TetraMesh_20201103(N,T(1,:),N2,T(1,:));
+ALGO_TetraIntersect_v1_20201103(V,U);
+
+figure(2)
+clf
+PLOT_TetraMesh_20201103(N,T(1,:),N2,T(1,:));
+% ALGO_ProjectionMatrix3D_v1_20201104(N,T,N2,T);
+ALGO_simplices_IntersectSimplices_v1_20211216(V,U);
\ No newline at end of file

Research/Edge Intersection/Tetrahedra/Higher dimensions/ALGO_simplices_IntersectSimplices_v1_20211216.m

-function [Z]=ALGO_simplices_IntersectSimplices_v1_20211216(U,V)
+function [Z,nhbrs]=ALGO_simplices_IntersectSimplices_v1_20211216(V,U)
 % INTERSECTSIMPLICES calculates the intersection between two simplices of
 % arbitrary dimension
 
@@ -10,17 +10,28 @@ function [Z]=ALGO_simplices_IntersectSimplices_v1_20211216(U,V)
 % m-face of X does not have a maximum number of intersections at a given
 % step.
 
+% Also missing is the compatibility with an advancing front algorithm.
+% Without it, this cannot be considered part of the PANG family. We need to
+% add a variable n that indicates which neighbours of X (or Y?) also
+% intersect Y (or X?).
+% On this matter, the order so far prescribed for the 2D and 3D simplices
+% is to have the e_0 plane last and the others in order.
+
 n=length(U); % there are n vertices of U (and V) which means the dimension
              %  of the space is n-1
              
 %--Change of Coordinates--%
-v0=V(:,n); A=V(:,1:n-1)-V(:,n);
+v0=V(:,1); A=V(:,2:end)-v0;
+global X
 X=A \ (U-v0);
-X=[X; 1-ones(1,n-1)*X]; % tack on the x*e_0 'coordinates'
+X=[1-ones(1,n-1)*X; X]; % tack on the x*e_0 'coordinates'
+                        % the row order in X must match the column order
+                        % for binV, ie. for the vertices of V
 
 %--Vertices of X in Y--%
+nhbrs=zeros(1,n);     % which neighbours of V also intersect U
 S=X>=0;
-if prod(prod(S,2))==0   % if there is a given coordinate where all X values are negative
+if prod(sum(S,2))==0   % if there is a given coordinate where all X values are negative
   disp('Simplices do not intersect') % then there is a hyperplane that entirely divides X and Y
   Z=[];
   return
@@ -29,19 +40,26 @@ S=prod(S);
 Z=U(:,S==1);
 
 %--Intersections--%
+global D
 D=zeros(2^n,2^(n+1)); % storage for determinants
 H=cell(1); J=H; T=H;
 H{1}=X; J{1}=eye(n); T{1}=zeros(n,1);
-for i=1:n-1
-  [H,J,T]=SUB_simplices_Unfold_v1_20211214(H,J,T,X,D);
+for i=1:n-2
+  [H,J,T]=SUB_simplices_Unfold_v1_20211214(H,J,T);
   if isempty(H)         % if no intersections are calculated for a given step then there are no more intersections to be found
     return
   end
-  S=prod(H>=0); Z=[Z,A*H(1:n,S==1)+v0];
+  for j=1:length(T)     % given this unpacking it is probably more efficient to do this step within the subfunctions
+    Hj=H{j}; S=prod(Hj>=0); Z=[Z,A*Hj(2:end,S==1)+v0];
+    if sum(S)~=0             % check which hyperplanes of V intersect with edges of U
+      ind=1:n; ind=ind(T{j}==1)-1; ind(ind==0)=n;
+      nhbrs(ind)=1;
+    end
+  end
 end
 
 %--Vertices of Y in X--%
-binV=zeros(n,1); indV=binV; % another probably terrible way of doing things... see notes on v1 of Unfold and Orphan for thoughts
+binV=zeros(n,1); indV=binV; % another probably terrible way of doing things... see notes on v1 of Unfold and v2 of Orphan for thoughts
 while ~isempty(H)
   H1=H{1}; T1=T{1};
   indT=1:n; indT=indT(T1==0); % at most two values in ind
@@ -60,6 +78,12 @@ while ~isempty(H)
       indT(i)=1;
     end
   end
-  H=H{indT==0}; T=T{indT==0};
+  H=H(indT==0); T=T(indT==0);
 end
-Z=[Z,V(:,binV)];
\ No newline at end of file
+Z=[Z,V(:,binV==1)];
+if sum(binV)>=n-1
+  nhbrs=ones(1,n);
+end
+
+PLOT_PatchPolyhedron_v1_20201105(Z)
+pause(0)
\ No newline at end of file

Research/Edge Intersection/Tetrahedra/Higher dimensions/SUB_simplices_IntersectHyperplane_v1_20211217.m

-function [h,D]=SUB_simplices_IntersectHyperplane_v1_20211217(D,J,T,X)
+function h=SUB_simplices_IntersectHyperplane_v1_20211217(J,T)
 % INTERSECTHYPERPLANE computes the intersection of an m-face of X and a
 % collection of hyperplanes of Y
 
 % Uses Cramer's rule with a look-up table to find determinants indexed by J
 % and T
 
+global X D % recall X and D from IntersectSimplices
+
 %--Initialize outputs--%
 n=length(T); m=nnz(T)+1;
 h=zeros(n,1);
@@ -16,15 +18,28 @@ indT=1:n; indT=indT(T==0);
 j=bit2int(J,n);
 k1=bit2int([1;T],n+1);
 if D(j,k1)==0
-    D(j,k1)=det([ones(m,1), X(J==1,T==1)]); % could also be calculated with other entries of D now that they're indexed with binary
+    D(j,k1)=det([ones(m,1), X(T==1,J==1)']); % could also be calculated with other entries of D now that they're indexed with binary
 end
 
 for i=indT
     %--Convert T with i to coordinates of D--%
+    count=sum(T(1:i));
     T_new=T; T_new(i)=1;
     k=bit2int([0;T_new],n+1);
     if D(j,k)==0
-        D(j,k)=det(X(J==1,T_new==1));
+        D(j,k)=det(X(T_new==1,J==1));
     end
-    h(i)=(-1)^(i-1) * D(j,k)/D(j,k1); % need to multiply by a sign to indicate column swaps between num and denom
+    h(i)=(-1)^count * D(j,k)/D(j,k1); % need to multiply by a sign to indicate column swaps between num and denom
+end
+end
+
+function k=bit2int(v,n)
+% BIT2INT converts binary string to integer
+
+k=0;
+for i=1:n
+    if v(i)==1
+        k=k + 2^(n-i);
+    end
+end
 end
\ No newline at end of file

Research/Edge Intersection/Tetrahedra/Higher dimensions/SUB_simplices_Orphan_v2_20211214.m

-function [H_new,J_new,T_new,D] = SUB_simplices_Orphan_v2_20211214(H,J,T,y,D,X)
+function [H_new,J_new,T_new] = SUB_simplices_Orphan_v2_20211214(H,J,T,y)
 % ORPHAN determines adjacent intersections and generates index set
 %   [H_new,J_new,T_new] = Orphan(H,J,T,A,y,X) finds all intersections
 %   H_new between the n-simplex X and the hyperplanes indexed by T and y.
@@ -9,7 +9,6 @@ function [H_new,J_new,T_new,D] = SUB_simplices_Orphan_v2_20211214(H,J,T,y,D,X)
 %           of H lie;
 %       T - hyperplanes of Y where the intersections H lie;
 %       y - next hyperplane to intersect the convex hull of H with;
-%       X - position of the vertices of X.
 %
 %   Outputs:
 %       H_new - intersections between X and Y indexed by J_new and T_new;
@@ -27,7 +26,7 @@ for i=1:m-1
             % calculate intersection
                 m_new=m_new+1;
                 J_new(:,m_new)=J(:,i) | J(:,j);
-                [H_new(:,m_new),D]=SUB_simplices_IntersectHyperplane_v1_20211217(D,J_new(:,m_new),T_new,X);
+                H_new(:,m_new)=SUB_simplices_IntersectHyperplane_v1_20211217(J_new(:,m_new),T_new);
             end
         end
     end

Research/Edge Intersection/Tetrahedra/Higher dimensions/SUB_simplices_Unfold_v1_20211214.m

-function [H_new,J_new,T_new] = SUB_simplices_Unfold_v1_20211214(H_all,J_all,T_all,X,D)
+function [H_new,J_new,T_new,m_new] = SUB_simplices_Unfold_v1_20211214(H_all,J_all,T_all)
 % UNFOLD marches through combinations of hyperplane intersections
 %   [H_new,J_new,T_new]=Unfold(H_all,J_all,T_all) computes the new set of
 %   intersections H_new with indices of X J_new and indices of Y T_new.
@@ -23,19 +23,23 @@ function [H_new,J_new,T_new] = SUB_simplices_Unfold_v1_20211214(H_all,J_all,T_al
 T=T_all{1}; n=length(T); d=nnz(T)+1; m_max=nchoosek(n,d);
 H_new=cell(m_max,1); J_new=H_new; T_new=H_new; m_new=0;
 % D=zeros(m_max,nchoosek(n+1,d)); % possibly a faster way to get second input from first
-while ~isempty(H_all)
-    H=H_all{1}; J=J_all{1}; T=T_all{1};
-    indy=1:n; indy=indy(T==0);
+m=length(T_all); indTy=zeros(n,m);
+for i=1:m
+    H=H_all{i}; J=J_all{i}; T=T_all{i};
+    indTy(T==0,i)=indTy(T==0,i)+1;
+    indy=1:n; indy=indy(indTy(:,i)==1);
     for y=indy
-        m_new=m_new+1;      % how many new hyperplane combos have we computed?
-        [H_new{m_new},J_new{m_new},T_new{m_new},D]=SUB_simplices_Orphan_v2_20211214(H,J,T,y,D,X);
+        [H_temp,J_temp,T_temp]=SUB_simplices_Orphan_v2_20211214(H,J,T,y);
+        if ~isempty(H_temp)
+            m_new=m_new+1;      % how many new hyperplane combos have we computed?
+            H_new{m_new}=H_temp; J_new{m_new}=J_temp; T_new{m_new}=T_temp;
+        end
     end
-    m=length(T_all); indT=zeros(m,1);
-    for i=1:m
-        if nnz(T_all{i}-T)<=2
-            indT(i)=1;      % which hyperplane combos are nearby?
+    for j=i+1:m
+        temp=T_all{j}-T;
+        if nnz(temp)==2
+            indTy(temp~=0,j)=-1; % remove hyperplane combinations already completed in this iteration
         end
     end
-    H_all=H_all(indT==0); J_all=J_all(indT==0); T_all=T_all(indT==0); % remove hyperplane combinations already completed in this iteration
 end
-H_new=H_new{1:m_new}; J_new=J_new{1:m_new}; T_new=T_new{1:m_new};
\ No newline at end of file
+H_new=H_new(1:m_new); J_new=J_new(1:m_new); T_new=T_new(1:m_new);
\ No newline at end of file

Research/Edge Intersection/Tetrahedra/mccoid2020tetrahedra.tex

@@ -1337,4 +1337,19 @@ where $I_\Gamma$ has columns $\vec{e}_\gamma$ for $\gamma \in \Gamma$.
 
 We could also calculate numerators and denominators separately and combine them as needed.
 
+\subsubsection{Notes on current implementation}
+
+Current implementation (IntersectSimplices v1, Unfold v1, Orphan v2 and IntersectHyperplane v1) is in working order and compatible with ProjectionMatrix3D.
+It is not clear if the unfolding process has any overlaps / is in any way efficient.
+Certainly the orphaning could be done more smoothly but we're stuck on certain graph problems.
+
+Also missing is the check of signs from previous generation of intersections indicating signs of next generation.
+If two intersections in a given generation have the same sign in a given coordinate then the sign of any intersection in the next generation with the same vertices has the same sign in that coordinate.
+As explained somewhere above if there are fewer than $m$ intersections between $m$ hyperplanes of $Y$ and an $m$--face of $X$ then the orientation of the face is predetermined by other intersections.
+
+The size of the matrix $D$ which holds the values of the various determinants is much larger than it needs to be.
+Its current size is determined by the easy way of translating between the binary strings for $J$ and $\Gamma$ and integer coordinates (literally convert binary string to integer).
+This leaves a lot of blank spaces in both directions.
+It is almost certainly preferable to a) come up with a new bijection between the binary strings and integers and b) use sparse matrix for $D$ (but still global?).
+
 \end{document}
\ No newline at end of file