Commit 93b1cb81 authored by Conor McCoid's avatar Conor McCoid
Browse files

Tetra: algo v1 for higher dim

parent 4961681e
...@@ -29,10 +29,11 @@ S=prod(S); ...@@ -29,10 +29,11 @@ S=prod(S);
Z=U(:,S==1); Z=U(:,S==1);
%--Intersections--% %--Intersections--%
D=zeros(2^n,2^(n+1)); % storage for determinants
H=cell(1); J=H; T=H; H=cell(1); J=H; T=H;
H{1}=X; J{1}=eye(n); T{1}=zeros(n,1); H{1}=X; J{1}=eye(n); T{1}=zeros(n,1);
for i=1:n-1 for i=1:n-1
[H,J,T]=SUB_simplices_Unfold_v1_20211214(H,J,T,X); [H,J,T]=SUB_simplices_Unfold_v1_20211214(H,J,T,X,D);
if isempty(H) % if no intersections are calculated for a given step then there are no more intersections to be found if isempty(H) % if no intersections are calculated for a given step then there are no more intersections to be found
return return
end end
......
function [h,D]=SUB_simplices_IntersectHyperplane_v1_20211217(D,J,T,X)
% 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
%--Initialize outputs--%
n=length(T); m=nnz(T)+1;
h=zeros(n,1);
%--Determine which coordinates of the intersection are needed--%
indT=1:n; indT=indT(T==0);
%--Denominator is the same for all--%
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
end
for i=indT
%--Convert T with i to coordinates of D--%
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));
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
end
\ No newline at end of file
function [H_new,J_new,T_new] = SUB_simplices_Orphan_v2_20211214(H,J,T,y,X) function [H_new,J_new,T_new,D] = SUB_simplices_Orphan_v2_20211214(H,J,T,y,D,X)
% ORPHAN determines adjacent intersections and generates index set % 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,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. % H_new between the n-simplex X and the hyperplanes indexed by T and y.
...@@ -27,7 +27,7 @@ for i=1:m-1 ...@@ -27,7 +27,7 @@ for i=1:m-1
% calculate intersection % calculate intersection
m_new=m_new+1; m_new=m_new+1;
J_new(:,m_new)=J(:,i) | J(:,j); J_new(:,m_new)=J(:,i) | J(:,j);
% H_new(:,m_new)=SUB_simplices_Intersection_v1_?(X,J_new(:,m_new),T_new); [H_new(:,m_new),D]=SUB_simplices_IntersectHyperplane_v1_20211217(D,J_new(:,m_new),T_new,X);
end end
end end
end end
......
function [H_new,J_new,T_new] = SUB_simplices_Unfold_v1_20211214(H_all,J_all,T_all,X) function [H_new,J_new,T_new] = SUB_simplices_Unfold_v1_20211214(H_all,J_all,T_all,X,D)
% UNFOLD marches through combinations of hyperplane intersections % 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 % [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. % intersections H_new with indices of X J_new and indices of Y T_new.
...@@ -22,12 +22,13 @@ function [H_new,J_new,T_new] = SUB_simplices_Unfold_v1_20211214(H_all,J_all,T_al ...@@ -22,12 +22,13 @@ 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); 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; 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) while ~isempty(H_all)
H=H_all{1}; J=J_all{1}; T=T_all{1}; H=H_all{1}; J=J_all{1}; T=T_all{1};
indy=1:n; indy=indy(T==0); indy=1:n; indy=indy(T==0);
for y=indy for y=indy
m_new=m_new+1; % how many new hyperplane combos have we computed? 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}]=SUB_simplices_Orphan_v2_20211214(H,J,T,y,X); [H_new{m_new},J_new{m_new},T_new{m_new},D]=SUB_simplices_Orphan_v2_20211214(H,J,T,y,D,X);
end end
m=length(T_all); indT=zeros(m,1); m=length(T_all); indT=zeros(m,1);
for i=1:m for i=1:m
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment