%Title: "A unimodular demand type which is not a basis change of = substitutes"=0A= %Last Edit:27/8/15=0A= %For further details, please contact:=0A= % Timothy O'Connor, timothy.oconnor@economics.ox.ac.uk=0A= % Elizabeth Baldwin, e.c.baldwin@lse.ac.uk=0A= =0A= %Introduction=0A= =0A= %This program checks whether or not there exists a basis change for a=0A= %SPECIFIC 4x9 matrix whose resulting basis change will have max 1 = positive and=0A= %max 1 negative entry in each column. As this is for a specific 4x9 = matrix,=0A= %this program is NOT immediately generalizable for any 4x9 matrix.=0A= =0A= %We start by initiliazing all of our variables (1) before finding all=0A= %possible combinations of the first four columns of our 4x9 matrix with=0A= %<=3D2 nonzero entries in (2) and (3). We throw out from = consideration any=0A= %matrix who would not have the first four columns being invertible (4). = (5)=0A= %takes the invertible matrices and fills out the rest of the columns.=0A= %(6),(7), and (8) will filter based on if a matrix generates >2 NZ = entries=0A= %in the later part of the matrix. (9) throws out any collection of = columns=0A= %that were originally invertible which are no longer invertible. (10) = will=0A= %then check to see if the last group of candidate matrices can have at = most=0A= %1 positive and at most 1 negative entry in each column.=0A= =0A= %Notation:=0A= %A =3D 4x4 matrix that features the first four columns of the matrix in = question=0A= %Dt =3D full matrix whose columns are the vectors of the demand type=0A= %P =3D basis change on Dt, so that P.Dt is the matrix that would have at = most one positive and at most one negative entry in each column.=0A= =0A= %The unimodular demand type, Dt, that we are investigating is given by:=0A= =0A= a=3D[1,0,0,0];a=3Da';=0A= b=3D[0,1,0,0];b=3Db';=0A= c=3D[0,0,1,0];c=3Dc';=0A= d=3D[1,0,0,1];d=3Dd';=0A= A=3D[a,b,c,d];=0A= Dt=3D[A,d-a+b,d-a+c,d+b,d+c,d-a+b+c];=0A= =0A= %% 1. Initialization of the variables/matrices=0A= =0A= =0A= %We start by creating the possible basis change matrix P by producing = symbolic variables=0A= %that are constrained to the reals.=0A= syms p1_1 real p1_2 real p1_3 real p1_4 real p2_1 real p2_2 real p2_3 = real p2_4 real p3_1 real p3_2 real p3_3 real p3_4 real p4_1 real p4_2 = real p4_3 real p4_4 real=0A= p=3D[p1_1,p1_2,p1_3,p1_4;p2_1,p2_2,p2_3,p2_4;p3_1,p3_2,p3_3,p3_4;p4_1,p4_= 2,p4_3,p4_4];=0A= =0A= %% 2. Initialization of set of matrices with number of nonzero entries = <=3D2=0A= % We gather all possible ways a 4x4 matrix can have at most 2 nonzero = entries in each column=0A= % Once we have found them, we will then assume this matrix has the form = "PA" and set the remaining entries of PA to zero.=0A= =0A= =0A= perms=3Dnchoosek(4,2);%Total number of ways to have 2 Nonzero (NZ) = entries in a column=0A= perms_1=3Dnchoosek(4,1);%Total number of ways to have 1 NZ entry in a = colum=0A= perm_point=3Dnchoosek((1:4),2);%List of combinations for the 2 NZ = selections=0A= perm1_point=3Dnchoosek((1:4),1);%List of combinations for the 1 NZ = selection=0A= max_perms=3Dperms^4+perms_1*perms^3*4+6*perms^2*perms_1^2+4*perms*perms_1= ^3+perms_1^4;%total amount of combinations=0A= max_matrix_collection=3Dzeros(4,4,max_perms);%Preallocating space=0A= =0A= =0A= =0A= %We collect all the possible matrices with 2 nonzero entries in each = column=0A= %such that we cycle through perm_point for each column=0A= n=3D1;=0A= for i=3D1:perms=0A= i_p=3Dperm_point(i,:);=0A= for j=3D1:perms=0A= j_p=3Dperm_point(j,:);=0A= for k=3D1:perms=0A= k_p=3Dperm_point(k,:);=0A= for h=3D1:perms=0A= h_p=3Dperm_point(h,:);=0A= blank_z=3Dzeros(4,4);=0A= blank_z(i_p,1)=3D1;%We input 1 as these will be the = locations that are NZ=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;%Collect the = matrix=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= =0A= %Now we collect all the combinations with 3 columns having 2 nonzero=0A= %entries and 1 column having 1 nonzero entry=0A= =0A= %4th column has 1 nonzero=0A= for i=3D1:perms=0A= i_p=3Dperm_point(i,:);=0A= for j=3D1:perms=0A= j_p=3Dperm_point(j,:);=0A= for k=3D1:perms=0A= k_p=3Dperm_point(k,:);=0A= for h=3D1:perms_1=0A= h_p=3Dperm1_point(h);=0A= blank_z=3Dzeros(4,4);%reset=0A= blank_z(i_p,1)=3D1;=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= %3rd column has 1 nonzero=0A= for i=3D1:perms=0A= i_p=3Dperm_point(i,:);=0A= for j=3D1:perms=0A= j_p=3Dperm_point(j,:);=0A= for k=3D1:perms_1=0A= k_p=3Dperm1_point(k);=0A= for h=3D1:perms=0A= h_p=3Dperm_point(h,:);=0A= blank_z=3Dzeros(4,4);%reset=0A= blank_z(i_p,1)=3D1;=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= %2nd column has 1 nonzero=0A= for i=3D1:perms=0A= i_p=3Dperm_point(i,:);=0A= for j=3D1:perms_1=0A= j_p=3Dperm1_point(j);=0A= for k=3D1:perms=0A= k_p=3Dperm_point(k,:);=0A= for h=3D1:perms=0A= h_p=3Dperm_point(h,:);=0A= blank_z=3Dzeros(4,4);%reset=0A= blank_z(i_p,1)=3D1;=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= %1st column has 1 nonzero=0A= for i=3D1:perms_1=0A= i_p=3Dperm1_point(i);=0A= for j=3D1:perms=0A= j_p=3Dperm_point(j,:);=0A= for k=3D1:perms=0A= k_p=3Dperm_point(k,:);=0A= for h=3D1:perms=0A= h_p=3Dperm_point(h,:);=0A= blank_z=3Dzeros(4,4);%reset=0A= blank_z(i_p,1)=3D1;=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= =0A= %Now we collect all the combinations with 2 columns having 2 nonzero=0A= %entries and 2 columnns having 1 nonzero entry=0A= =0A= %1st and 2nd have 1 nonzero=0A= for i=3D1:perms_1=0A= i_p=3Dperm1_point(i);%1st column has 1=0A= for j=3D1:perms_1=0A= j_p=3Dperm1_point(j);%2nd column has 1=0A= for k=3D1:perms=0A= k_p=3Dperm_point(k,:);=0A= for h=3D1:perms=0A= h_p=3Dperm_point(h,:);=0A= blank_z=3Dzeros(4,4);%reset=0A= blank_z(i_p,1)=3D1;=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= %1st and 3rd have 1 nonzero=0A= for i=3D1:perms_1=0A= i_p=3Dperm1_point(i);=0A= for j=3D1:perms=0A= j_p=3Dperm_point(j,:);=0A= for k=3D1:perms_1=0A= k_p=3Dperm1_point(k);=0A= for h=3D1:perms=0A= h_p=3Dperm_point(h,:);=0A= blank_z=3Dzeros(4,4);%reset=0A= blank_z(i_p,1)=3D1;=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= %1st and 4th have 1 nonzero=0A= for i=3D1:perms_1=0A= i_p=3Dperm1_point(i);=0A= for j=3D1:perms=0A= j_p=3Dperm_point(j,:);=0A= for k=3D1:perms=0A= k_p=3Dperm_point(k,:);=0A= for h=3D1:perms_1=0A= h_p=3Dperm1_point(h);=0A= blank_z=3Dzeros(4,4);%reset=0A= blank_z(i_p,1)=3D1;=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= %2nd and 3rd have 1 nonzero=0A= for i=3D1:perms=0A= i_p=3Dperm_point(i,:);=0A= for j=3D1:perms_1=0A= j_p=3Dperm1_point(j);=0A= for k=3D1:perms_1=0A= k_p=3Dperm1_point(k);=0A= for h=3D1:perms=0A= h_p=3Dperm_point(h,:);=0A= blank_z=3Dzeros(4,4);%reset=0A= blank_z(i_p,1)=3D1;=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= %2nd and 4th have 1 nonzero=0A= for i=3D1:perms=0A= i_p=3Dperm_point(i,:);=0A= for j=3D1:perms_1=0A= j_p=3Dperm1_point(j);=0A= for k=3D1:perms=0A= k_p=3Dperm_point(k,:);=0A= for h=3D1:perms_1=0A= h_p=3Dperm1_point(h);=0A= blank_z=3Dzeros(4,4);%reset=0A= blank_z(i_p,1)=3D1;=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= %3rd and 4th have 1 nonzero=0A= for i=3D1:perms=0A= i_p=3Dperm_point(i,:);=0A= for j=3D1:perms=0A= j_p=3Dperm_point(j,:);=0A= for k=3D1:perms_1=0A= k_p=3Dperm1_point(k);=0A= for h=3D1:perms_1=0A= h_p=3Dperm1_point(h);=0A= blank_z=3Dzeros(4,4);%reset=0A= blank_z(i_p,1)=3D1;=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= =0A= %And now we collect all the combinations with 3 columns having 1 nonzero=0A= %and one column having 2 nonzero entries=0A= =0A= %4th has 2 nonzero=0A= for i=3D1:perms_1=0A= i_p=3Dperm1_point(i);=0A= for j=3D1:perms_1=0A= j_p=3Dperm1_point(j);=0A= for k=3D1:perms_1=0A= k_p=3Dperm1_point(k);=0A= for h=3D1:perms=0A= h_p=3Dperm_point(h,:);=0A= blank_z=3Dzeros(4,4);%reset=0A= blank_z(i_p,1)=3D1;=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= %3rd has 2 nonzero=0A= for i=3D1:perms_1=0A= i_p=3Dperm1_point(i);=0A= for j=3D1:perms=0A= j_p=3Dperm_point(j,:);=0A= for k=3D1:perms_1=0A= k_p=3Dperm1_point(k);=0A= for h=3D1:perms_1=0A= h_p=3Dperm1_point(h);=0A= blank_z=3Dzeros(4,4);%reset=0A= blank_z(i_p,1)=3D1;=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= %2nd has 2 nonzero=0A= for i=3D1:perms_1=0A= i_p=3Dperm1_point(i);=0A= for j=3D1:perms_1=0A= j_p=3Dperm1_point(j);=0A= for k=3D1:perms=0A= k_p=3Dperm_point(k,:);=0A= for h=3D1:perms_1=0A= h_p=3Dperm1_point(h);=0A= blank_z=3Dzeros(4,4);%reset=0A= blank_z(i_p,1)=3D1;=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= %1st has 2 nonzero=0A= for i=3D1:perms=0A= i_p=3Dperm_point(i,:);=0A= for j=3D1:perms_1=0A= j_p=3Dperm1_point(j);=0A= for k=3D1:perms_1=0A= k_p=3Dperm1_point(k);=0A= for h=3D1:perms_1=0A= h_p=3Dperm1_point(h);=0A= blank_z=3Dzeros(4,4);%reset=0A= blank_z(i_p,1)=3D1;=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= =0A= %Now collect with all four having 1 NZ=0A= for i=3D1:perms_1=0A= i_p=3Dperm_point(i);=0A= for j=3D1:perms_1=0A= j_p=3Dperm1_point(j);=0A= for k=3D1:perms_1=0A= k_p=3Dperm1_point(k);=0A= for h=3D1:perms_1=0A= h_p=3Dperm1_point(h);=0A= blank_z=3Dzeros(4,4);%reset=0A= blank_z(i_p,1)=3D1;=0A= blank_z(j_p,2)=3D1;=0A= blank_z(k_p,3)=3D1;=0A= blank_z(h_p,4)=3D1;=0A= max_matrix_collection(:,:,n)=3Dblank_z;=0A= n=3Dn+1;=0A= end=0A= end=0A= end=0A= end=0A= =0A= %And now we have all of the matrices.=0A= %Let us now put it into our actual format (with the pi_j's)=0A= =0A= %% 3. Now we take the set of matrices and express it in "PA" matrix form.=0A= =0A= %Preallocate space for the combinations of "PA" - Warning, this is where = the=0A= %computing time starts to increase as we loop over symbolic variables=0A= A_combos=3D sym(zeros(4,4,max_perms));=0A= =0A= %This gets us our "PA" matrix for all the combinations of columns with = <=3D2=0A= %nonzero entries=0A= for n=3D1:max_perms=0A= for i=3D1:4=0A= for j=3D1:4=0A= if(max_matrix_collection(i,j,n)=3D=3D1)%If this entry is NZ, = then input the correct pi_j value=0A= A_combos(i,j,n)=3Ddot(p(i,:),A(:,j));=0A= else=0A= A_combos(i,j,n)=3D0;=0A= end=0A= end=0A= end=0A= end=0A= =0A= %This makes sure that the appropriate substitutions are in place. If pi_1=0A= %is equal to zero and if pi_1+pi_4 is NZ then we should just have pi_4 = by itself=0A= for n=3D1:max_perms=0A= for i=3D1:4=0A= for j=3D1:3=0A= if A_combos(i,j,n)=3D=3D0=0A= A_combos(:,:,n)=3Dsubs(A_combos(:,:,n),p(i,j),0);=0A= end=0A= end=0A= end=0A= end=0A= =0A= %% 4. Now we filter all the "PA" matrices that would not be invertible =0A= =0A= %Preallocating space for our invertible collection=0A= count=3D0;=0A= for n=3D1:max_perms=0A= if rank(A_combos(:,:,n))=3D=3D4=0A= count=3Dcount+1;=0A= end=0A= end=0A= inv_A_combos=3D sym(zeros(4,4,count));=0A= =0A= count=3D0;=0A= for n=3D1:max_perms=0A= if rank(A_combos(:,:,n))=3D=3D4%If this matrix is invertible, we add = it to our list=0A= count=3Dcount+1;=0A= inv_A_combos(:,:,count)=3DA_combos(:,:,n);=0A= end=0A= end=0A= =0A= %% 5. Now we fill out our collection to include columns 5:9. That is, = we now work with the collection of matrices of the form "P.Dt".=0A= full_combos=3Dsym(zeros(4,9,count));=0A= for n=3D1:count=0A= %From our original matrix, if our first four columns are: a,b,c,d, = then:=0A= = c5=3Dinv_A_combos(:,4,n)-inv_A_combos(:,1,n)+inv_A_combos(:,2,n);%column = five is d-a+b=0A= = c6=3Dinv_A_combos(:,4,n)-inv_A_combos(:,1,n)+inv_A_combos(:,3,n);%column = six is d-a+c=0A= c7=3Dinv_A_combos(:,4,n)+inv_A_combos(:,2,n);%c7 is d+b=0A= c8=3Dinv_A_combos(:,4,n)+inv_A_combos(:,3,n);%c8 is d+c=0A= = c9=3Dinv_A_combos(:,4,n)-inv_A_combos(:,1,n)+inv_A_combos(:,2,n)+inv_A_co= mbos(:,3,n);%c9 is d-a+b+c=0A= filling=3D[c5,c6,c7,c8,c9];=0A= full_combos(:,:,n)=3D[inv_A_combos(:,:,n),filling];=0A= end=0A= =0A= %% 6. Now we check the amount of NZ entries in each column.=0A= =0A= alive=3D0;%new count of matrices that "survive" the filter=0A= =0A= %We take our set and filter out the matrices with greater than=0A= %2 NZ entries within a column. One important thing to note that if we = have=0A= %pi_j, for j=3D1,2,3, by itself in an entry, we know that it is NZ. = However,=0A= %this is not necessarily the case if we have pi_4 by itself. If all that = we=0A= %know is that pi_1+pi_4 is NZ, this does not tell us whether or not pi_4 = is=0A= %zero. Therefore, when we count the number of NZ entries in each column, = we=0A= %can add to our count if we have a pi_j, j~=3D4, by itself, and pi_4 will=0A= %only increase our count if pi_1 is zero.=0A= =0A= % It is much faster to first count how many matrices will remain after = this =0A= % step, preallocate the required space, and then collect those matrices. = Here=0A= % we perform that first count.=0A= for k=3D1:count=0A= flag=3D0;%Will flag on if we have more than two nonzero entries in a = column. We will then not add that matrix to our next list.=0A= for j=3D5:9=0A= nz=3D0;%A counting term for the number of nonzero entries=0A= for i=3D1:4=0A= str=3Dchar(full_combos(i,j,k));%We string the entry=0A= len=3Dlength(str);%Find the length of the string so that we = can see what the last character is=0A= if len<6&&len>1%Provided that the entry is not = just a "0"->length=3D=3D1, or an additive entry->length>5=0A= if str(len)~=3D'4'%If we have a pi_j by itself and = j~=3D4, we know it is a NZ entry=0A= nz=3Dnz+1;=0A= elseif full_combos(i,j,k)=3D=3Dfull_combos(i,4,k)%If = pi_4 is the only entry in the 4th column, then pi_1=3D0 and so pi_4 is NZ=0A= nz=3Dnz+1;=0A= end=0A= elseif len>6%We know look for a double that we know would = be NZ. The only double that we would know is NZ would be the double = found in the fourth column (pi_1+pi_4)=0A= if full_combos(i,j,k)=3D=3Dfull_combos(i,4,k)=0A= nz=3Dnz+1;=0A= end=0A= end=0A= end=0A= if nz>2=0A= flag=3D1;=0A= end=0A= if full_combos(:,j,k)=3D=3D[0;0;0;0]%After step five of filling = out the rest of the columns, there may be a column with all zeros=0A= flag=3D1;=0A= end=0A= end=0A= if flag=3D=3D0%If we do not have an columns with >2 NZ entries = then we increase our count.=0A= alive=3Dalive+1;=0A= end=0A= end=0A= =0A= %Now that we have the count, preallocate, then write in data - this next = part is the same loop as above, we just write in the data now.=0A= first_filter_num=3Dalive;=0A= first_filter=3Dsym(zeros(4,9,first_filter_num));=0A= alive=3D0;=0A= for k=3D1:count=0A= flag=3D0;=0A= for j=3D5:9=0A= nz=3D0;=0A= for i=3D1:4=0A= str=3Dchar(full_combos(i,j,k));=0A= len=3Dlength(str);=0A= if len<6&&len>2=0A= if str(len)~=3D'4'=0A= nz=3Dnz+1;=0A= elseif full_combos(i,j,k)=3D=3Dfull_combos(i,4,k)=0A= nz=3Dnz+1;=0A= end=0A= elseif len>6=0A= if full_combos(i,j,k)=3D=3Dfull_combos(i,4,k)=0A= nz=3Dnz+1;=0A= end=0A= end=0A= end=0A= if nz>2=0A= flag=3D1;=0A= end=0A= if full_combos(:,j,k)=3D=3D[0;0;0;0]=0A= flag=3D1;=0A= end=0A= end=0A= if flag=3D=3D0=0A= alive=3Dalive+1;=0A= first_filter(:,:,alive)=3Dfull_combos(:,:,k);=0A= end=0A= end=0A= =0A= %% 7. Now subtitute for the matrices with >2 entries when 2 are for = sure NZ=0A= =0A= %If a column has for sure 2 nonzero entries and then other entries that = are=0A= %not formally set to zero, we can now formally set to zero the other=0A= %entries. For example, the entries in a column are:=0A= %[p1_2+p1_4,p2_1,p3_3,p4_2+p2_4]. As the second and third entry are for=0A= %sure NZ we can formally set to zero the values of the first and fourth = entry (throughout=0A= %the entire matrix). Our column is then: [0,p2_1,p3_3,0].=0A= =0A= for k=3D1:first_filter_num=0A= for j=3D5:9=0A= nz=3D0;=0A= flag=3D0;=0A= for i=3D1:4=0A= str=3Dchar(first_filter(i,j,k));%As before, we string the = entries and inspect the last character. If we have pi_4, then it will = only be for sure NZ if we know that pi_1 is zero=0A= len=3Dlength(str);%Store the location of the last character=0A= if len<6&&len>1%If the length is greater than = 1 such that it is not "0" and less than 6 such that it is a single = element (pi_j or -pi_j)=0A= if str(len)~=3D'4'=0A= nz=3Dnz+1;=0A= elseif first_filter(i,j,k)=3D=3Dfirst_filter(i,4,k)%Else = if pi_4 is by itself in the fourth column:=0A= nz=3Dnz+1;=0A= end=0A= elseif len>6%We know look for a double that we know would = be NZ. The only double that we would know is NZ would be the double = found in the fourth column (pi_1+pi_4)=0A= if first_filter(i,j,k)=3D=3Dfirst_filter(i,4,k)=0A= nz=3Dnz+1;=0A= end=0A= end=0A= end=0A= if nz=3D=3D2=0A= flag=3D1;=0A= end=0A= if flag=3D=3D1=0A= %If we have two nonzero entries in this column, we will not = go=0A= %back and formally set to zero any entries that we did not=0A= %formally know if they were NZ or not=0A= for i=3D1:4=0A= str=3Dchar(first_filter(i,j,k));=0A= len=3Dlength(str);=0A= if len>6 && = first_filter(i,j,k)~=3Dfirst_filter(i,4,k)%If this entry is a double and = it is NOT the double whose formal value we could know (ie, not pi_1+pi_4)=0A= = first_filter(:,:,k)=3Dsubs(first_filter(:,:,k),first_filter(i,j,k),0);%Th= en we formally set this entry to zero and we make this substitution = throughout the entire matrix=0A= elseif len<6 &&len>1 && = str(len)=3D=3D'4' && = first_filter(i,j,k)~=3Dfirst_filter(i,4,k)%If the entry is pi_4 where we = did not know whether it was NZ prior. This would be the case if we only = know that pi_1+pi_4 is NZ=0A= = first_filter(:,:,k)=3Dsubs(first_filter(:,:,k),first_filter(i,j,k),0);%Th= en we formally set set pi_4 to zero and we make this substitution = throughout the entire matrix=0A= end=0A= end=0A= end=0A= end=0A= end=0A= =0A= %% 8. Now we filter again based on >2 NZ=0A= =0A= %After our substitutions made in 7, we once again filter based on a count=0A= %of the NZ entries in each column. This step is identical to step 6.=0A= alive=3D0;=0A= for k=3D1:first_filter_num=0A= flag=3D0;=0A= for j=3D5:9=0A= nz=3D0;=0A= for i=3D1:4=0A= str=3Dchar(first_filter(i,j,k));=0A= len=3Dlength(str);=0A= if len<6&&len>2=0A= if str(len)~=3D'4'=0A= nz=3Dnz+1;=0A= elseif first_filter(i,j,k)=3D=3Dfirst_filter(i,4,k)=0A= nz=3Dnz+1;=0A= end=0A= elseif len>6=0A= if first_filter(i,j,k)=3D=3Dfirst_filter(i,4,k)=0A= nz=3Dnz+1;=0A= end=0A= end=0A= end=0A= if nz>2=0A= flag=3D1;=0A= end=0A= end=0A= if flag=3D=3D0=0A= alive=3Dalive+1;=0A= end=0A= end=0A= second_filter_num=3Dalive;=0A= second_filter=3Dsym(zeros(4,9,second_filter_num));=0A= =0A= alive=3D0;=0A= for k=3D1:first_filter_num=0A= flag=3D0;=0A= for j=3D5:9=0A= nz=3D0;=0A= for i=3D1:4=0A= str=3Dchar(first_filter(i,j,k));=0A= len=3Dlength(str);=0A= if len<6&&len>2=0A= if str(len)~=3D'4'=0A= nz=3Dnz+1;=0A= elseif first_filter(i,j,k)=3D=3Dfirst_filter(i,4,k)=0A= nz=3Dnz+1;=0A= end=0A= elseif len>6=0A= if first_filter(i,j,k)=3D=3Dfirst_filter(i,4,k)=0A= nz=3Dnz+1;=0A= end=0A= end=0A= end=0A= if nz>2=0A= flag=3D1;=0A= end=0A= end=0A= if flag=3D=3D0=0A= alive=3Dalive+1;=0A= second_filter(:,:,alive)=3Dfirst_filter(:,:,k);=0A= end=0A= end=0A= =0A= %% 9. Now filter on Invertible grounds=0A= =0A= =0A= %From this set, we know that every invertible subset of the original 9=0A= %vectors (the columns of Dt) must also be invertible after Dt has been = acted on by P. This is because the product of two invertible matrices is=0A= %invertible.=0A= possib=3Dnchoosek(9,4);%The amount of ways that one can make a 4x4 = matrix from 9 columns=0A= possib_list=3Dnchoosek(1:9,4);%The combinations of columns to make a 4x4 = matrix=0A= count=3D0;=0A= =0A= for n=3D1:possib=0A= check=3DDt(:,possib_list(n,:));=0A= if rank(check)=3D=3D4%If it is invertible we increase our count=0A= count=3Dcount+1;=0A= end=0A= end=0A= =0A= inv_num=3Dcount;=0A= inv_list=3Dzeros(inv_num,4);%Preallocate the space for the lists of = combinations of columns that are invertible.=0A= count=3D0;=0A= for n=3D1:possib=0A= check=3DDt(:,possib_list(n,:));=0A= if rank(check)=3D=3D4=0A= count=3Dcount+1;=0A= inv_list(count,:)=3Dpossib_list(n,:);%Stores the combination of = columns that are invertible=0A= end=0A= end=0A= =0A= %Now that we know which combinations of columns are invertible in our=0A= %original matrix, we check to see if the SAME combination of columns is=0A= %invertible within our list of candidate matrices. =0A= count=3D0;=0A= for k=3D1:second_filter_num=0A= flag=3D0;=0A= for c=3D1:inv_num=0A= mat=3Dsecond_filter(:,inv_list(c,:),k);%Selects the combination = of columns that should be invertible=0A= if rank(mat)~=3D4%If this is not invertible, we flag.=0A= flag=3D1;=0A= end=0A= end=0A= if flag=3D=3D0%If each combination is invertible, we add it to our = count=0A= count=3Dcount+1;=0A= end=0A= end=0A= =0A= third_filter_num=3Dcount;=0A= third_filter=3Dsym(zeros(4,9,third_filter_num));%Preallocate our space = and run the loop once more to store the information=0A= count=3D0;=0A= for k=3D1:second_filter_num=0A= flag=3D0;=0A= for c=3D1:inv_num=0A= mat=3Dsecond_filter(:,inv_list(c,:),k);=0A= if rank(mat)~=3D4=0A= flag=3D1;=0A= end=0A= end=0A= if flag=3D=3D0=0A= count=3Dcount+1;=0A= third_filter(:,:,count)=3Dsecond_filter(:,:,k);=0A= end=0A= end =0A= =0A= %% 10. And now we filter among the combinations that cannot have at most = 1 positive and at most 1 negative value in each column=0A= =0A= %The general idea is that we will go through each column and make a=0A= %"relationship" between pairs of entries in each column. For instance,=0A= %if one column is [0,p2_2,0,p4_4]', we know that p2_2 and p4_4 must have=0A= %opposite signs (if one is positive, the other is negative). If the = column=0A= %was rather [0,p2_2,0,-p4_4], we would then know that p2_2 and p4_4 have = the=0A= %same sign (both are positive or both are negative).=0A= =0A= %We build these relationships for each column and see if there is a = contradiction.=0A= %For example, imagine that our 4x9 matrix includes the following column = vectors:=0A= %=0A= % [0,p2_2,0,p4_4]' (1)=0A= % [0,0,p3_3,p4_4]' (2)=0A= % [0,p2_2,p3_3,0]' (3)=0A= %=0A= %From (1) we know that p2_2 and p4_4 are of opposite sign. From (2) we = know that=0A= %p3_3 and p4_4 are of opposite sign. We then have a contradiction in (3) = as =0A= %(3) says that p2_2 and p3_3 are of opposite sign, and yet (1) and (2) = combined=0A= %say that p2_2 and p3_3 must be of the same sign. This sort of = contradiction provides=0A= %the foundation for our final filter.=0A= =0A= %We first make sure that there are 2 entries in each column that are not = formally zero. We=0A= %do not in fact need two non-formally zero entries in each column, but a=0A= %"relationship" can only be formed when there are two non-formally zero = entries. It could=0A= %be the case that there are 3 symbolic zeros, but as it so happens, after=0A= %all of the filtering it is the case that we have exactly two entries = that=0A= %are formally zero. We show this here:=0A= check=3D0;=0A= for n=3D1:third_filter_num=0A= for j=3D1:9=0A= nz_num=3Dlength(find(third_filter(:,j,n)));%counts the zeros in = the column=0A= if nz_num~=3D2%If we do not have two entries formally set to = zero:=0A= check=3Dcheck+1;=0A= end=0A= end=0A= end=0A= =0A= %As one can see, check=3D=3D0 such that we have exactly two non-formally = zero entries=0A= %in each column. =0A= =0A= %A second check that we will have to make (to ensure that the next loop = is=0A= %specified correctly) is to make sure that we do not have pi_4's alone = in a=0A= %column where we do not know whether or not it is NZ. As we add=0A= %relationships based on having the same or opposite sign, we want to make=0A= %sure that we do not make a comparison with pi_4 when it may be zero.=0A= %Therefore:=0A= check2=3D0;=0A= for k=3D1:third_filter_num=0A= flag=3D0;=0A= for j=3D5:9=0A= for i=3D1:4=0A= str=3Dchar(third_filter(i,j,n));%As before, we string the = entries and inspect the last character.=0A= len=3Dlength(str);%Store the location of the last character=0A= if len<6&&len>1%If the length is greater than = 1 such that it is not "0" and less than 6 such that it is a single = element (pi_j or -pi_j)=0A= if str(len)=3D=3D'4'=0A= if third_filter(i,j,n)~=3Dthird_filter(i,4,n)%Such = that whether pi_4=3D0 is then unknown=0A= flag=3D1;=0A= end=0A= end=0A= end=0A= end=0A= end=0A= if flag=3D=3D1=0A= check2=3Dcheck2+1;=0A= end=0A= end=0A= %As we can see, check2=3D=3D0. Therefore, we do not have any = indeterminate=0A= %pi_4's. All the pi_4's that are alone in a column will be definitively = NZ.=0A= =0A= =0A= for n=3D1:third_filter_num=0A= %We first build our collection of same/opposite relations between our=0A= %individual elements. "Opps" is a list of element pairs that we know=0A= %are opposite sign from each other. "Same" is a list of element pairs=0A= %that we know have the same sign. Note that we only add single = element=0A= %relations. We shall skip doubles as we will have enough information=0A= %from the single element pairs.=0A= =0A= =0A= flag=3D0;%Will flag on if we have a contradiction=0A= opps_n=3D0;%The count of opposite relations=0A= same_n=3D0;%The count of same relations=0A= for j=3D1:9=0A= indx=3Dfind(third_filter(:,j,n)~=3D0);%We find all of our = nonzero entries=0A= str1=3Dchar(third_filter(indx(1),j,n));%We string each entry=0A= str2=3Dchar(third_filter(indx(2),j,n));=0A= if length(str1)+length(str2)=3D=3D9 %With a length of nine, we = know that only one of the entries has a negative sign in front=0A= same_n=3Dsame_n+1;%As one has a negative sign, we know that = the elements must share the SAME sign=0A= = same(same_n,:)=3D[third_filter(indx(1),j,n),third_filter(indx(2),j,n)];%W= e add this relation to our SAME list=0A= elseif = length(str1)+length(str2)=3D=3D8||length(str1)+length(str2)=3D=3D10=0A= %If the length is 8, neither element has a negative sign. If = it=0A= %is 10, they both have a negative sign. Either way, we know=0A= %that the elements then must be of OPPOSITE sign=0A= opps_n=3Dopps_n+1;=0A= = opps(opps_n,:)=3D[third_filter(indx(1),j,n),third_filter(indx(2),j,n)];=0A= end=0A= end=0A= %Now we clean so that we do not get negative signs in front of = elements=0A= for i=3D1:same_n=0A= for j=3D1:2=0A= negcheck=3Dchar(same(i,j));=0A= if length(negcheck)=3D=3D5=0A= same(i,j)=3D(-1)*same(i,j);=0A= end=0A= end=0A= end=0A= for i=3D1:opps_n=0A= for j=3D1:2=0A= negcheck=3Dchar(opps(i,j));=0A= if length(negcheck)=3D=3D5=0A= opps(i,j)=3D(-1)*opps(i,j);=0A= end=0A= end=0A= end=0A= %Now we do a first check to see whether or not there will be any=0A= %contradictions directly=0A= flip=3D[opps(:,2),opps(:,1)];%also check the flip as p1_1, p1_2 = would be the same as p2_1, p1_1=0A= check=3Dintersect(same,opps,'rows');=0A= check2=3Dintersect(same,flip,'rows');=0A= if ~isempty(check)||~isempty(check2)=0A= %If the intersection of the same and opposite relation is = non-empty=0A= %such that a pair of elements is said to have both the same and=0A= %opposite sign, we know this matrix fails.=0A= flag=3D1;%And then we will skip the next loop and move on to the = next matrix.=0A= end=0A= =0A= %We will now add to our list of same/opposites via transitivity=0A= loop_count=3D1;=0A= %The loop_count will count the number of times that we cycle through=0A= %our same and opposite lists, adding new relationship via = transitivity.=0A= %We have an arbitrary limit on the loop count so that we do not loop=0A= %indefinitely if there was some error prior.=0A= while flag=3D=3D0&&loop_count<5=0A= for i=3D1:length(opps)-1=0A= %As we will look at occurences of a selected element later = in =0A= %the same and opposite lists, we do not look for more =0A= %occurences when we are at the end of the list.=0A= for j=3D1:2=0A= search=3Dopps(i,j);=0A= %'Search' is the first element and we will look for more=0A= %occurences of this element down the list in opposite and=0A= %same. 'Pairing' is the element that we know has an=0A= %opposite relationship with 'search'. =0A= if j=3D=3D1=0A= pairing=3Dopps(i,2);=0A= else=0A= pairing=3Dopps(i,1);=0A= end=0A= newopps=3Dopps((i+1):length(opps),1:2);%Searches down = the list from where we are currently=0A= [r,c]=3Dfind(newopps=3D=3Dsearch);%returns index within = the opposite to build up same, as the opposite of opposites is the same=0A= [r1,c1]=3Dfind(same=3D=3Dsearch);%returns index within = same to build up opps, as the opposite of same is opposite. = =0A= if ~isempty(r)%if we find an entry down the opposite = list such that 'search' is located somewhere down the list=0A= new_rels=3Dsym(zeros(length(r),2));%Then the new = relations that we make will be added to same.=0A= for b=3D1:length(r)%As we may have more than = occurence of the 'search' element=0A= if c(b)=3D=3D2%If 'search' is in the second = column, make the new relationship with the element in the first column=0A= = new_rels(b,1:2)=3D[pairing,newopps(r(b),1)];%stores the new relationship=0A= else%If 'search' is in the first column, make = the new relationship with the element in the second column=0A= new_rels(b,1:2)=3D[pairing,newopps(r(b),2)];=0A= end=0A= end=0A= same=3D[same;new_rels];%As the opposite of opposite = is same, we add to our same list.=0A= end=0A= if ~isempty(r1)%if we find an entry down the same list...=0A= new_rels=3Dsym(zeros(length(r),2));%Then the new = relations that we make will be added to opposite.=0A= for b=3D1:length(r1)=0A= if c1(b)=3D=3D2=0A= = new_rels(b,1:2)=3D[pairing,same(r1(b),1)];%stores the new relationship=0A= else=0A= new_rels(b,1:2)=3D[pairing,same(r1(b),2)];=0A= end=0A= end=0A= opps=3D[opps;new_rels];%As the opposite sign of the = same sign is opposite, we add to our opposite list.=0A= end=0A= end=0A= end=0A= flip=3D[opps(:,2),opps(:,1)];%We flip the list again before = checking intersections=0A= check=3Dintersect(same,opps,'rows');=0A= check2=3Dintersect(same,flip,'rows');=0A= if ~isempty(check)||~isempty(check2)=0A= %If the intersection of the same and opposite relation is = non-empty=0A= %such that a pair of elements is said to have both the same and=0A= %opposite sign, we know this matrix fails.=0A= flag=3D1;=0A= end=0A= loop_count=3Dloop_count+1;=0A= end =0A= =0A= if flag=3D=3D0%if we looped through 5 times and still were not able = to find a contradiction..=0A= candidate_matrix=3Dn;=0A= disp('Error: Matrix #n was not thrown out of consideration');=0A= %Therefore, if there is a printed message, we know that we were=0A= %unable to throw out the n'th matrix in third_filter.=0A= end=0A= end=0A= =0A= =0A= =0A= =0A= =0A= =0A=