function [similar] = PerSim(A,B)

% function [similar] = PerSim(A,B)
%
% (1) A and B are 16x16 matrices
% (2) A and B must satisfy A^2 = J and B^2 = J
% (3) A and B have previously known row multiplicities
%       (program is dependent on known row multiplicities)
% (4) function determines if A and B are permutationally
%       similar by first permuted the matrices into
%       decreasing row multiplicity (DRM) form.
% (5) similar equals:
%       1 iff matrices are permutationally similar
%       0 iff matrices are not similar
%       -1 iff matrices fail to conform to above restrictions
%
% Definition: A matrix in DRM form has the property that
%   if row i and row j are equivalent, then either i = j+1,
%   i = j-1, or i = j +/- c where all rows between i and j
%   are equal.  Furthermore, sets of equivalent rows appear
%   in decreasing size from the top to the bottom of the
%   matrix.  In other words, if a set of 3 equivalent rows
%   and a set of 4 equivalent rows appear in a matrix, then
%   the set of 4 rows (rows a, b, c, and d) and 3 rows (rows
%   x, y, and z) have the property that a,b,c,d < x,y,z.
%
% Example of 4x4 matrices in DRM form:
%
%   A = [1 0 1 0;    B = [0 1 0 1;
%        1 0 1 0;         0 1 0 1;
%        0 1 0 1;         1 0 1 0;
%        0 1 0 1;]        1 0 1 0;]
%
% Proposition: permuting A and B into DRM form dramatically
%   decreases the number of permutations necessary to assure
%   A is not permutationally similar to B
%
% NOTE: Program has not yet been completed for all valid
%   row multiplicities!
%

n = 16;
k = sqrt(n);

% check if matrices are correct size

[mA nA] = size(A);
[mB nB] = size(B);

if (mA ~= n)
    reason_for_failure = 'A must be nxn';
    similar = -1;
    return
else
    if (nA ~= n)
        reason_for_failure = 'A must be nxn';
        similar = -1;
        return
    else
        if (mB ~= n)
            reason_for_failure = 'B must be nxn';
            similar = -1;
            return
        else
            if (nB ~= n)
                reason_for_failure = 'B must be nxn';
                similar = -1;
                return
            end
        end
    end
end

% check if given matrix completions yield matrices with same rank
% Note: different rank implies that matrices are not similar

rA = rank(A);
rB = rank(B);

if (rA ~= rB)
    reason_for_failure = 'Matrices have different rank';
    similar = 0;
    return
end

% check if given matrix completions yield matrices that satisfy A^2 = J

squareA = A*A;
squareB = B*B;
squareAsameJ = 1;
squareBsameJ = 1;

for i = 1:n
    for j = 1:n
        if (squareA(i,j) ~= 1)
            squareAsameJ = 0;
        end
        if (squareB(i,j) ~= 1)
            squareBsameJ = 0;
        end
    end
    if (squareAsameJ == 0)
        if (squareBsameJ == 0)
            break
        end
    end
end

if (squareAsameJ == 0)
    if (squareBsameJ == 0)
        reason_for_failure = 'A^2 ~= J and B^2 ~= J';
        similar = -1;
        return
    else
        reason_for_failure = 'A^2 ~= J';
        similar = -1;
        return
    end
else
    if (squareBsameJ == 0)
        reason_for_failure = 'B^2 ~= J';
        similar = -1;
        return
    end
end

% We consider row and column multiplicities and permute
%    A and B into DRM form

uniqueRowMatrixA = zeros(n,n); % stores a copy for each unique row in A
multVectorA = zeros(1,n);      % stores number of rows equal to corresponding row in uniqueRowMatrixA
multMatrixA = zeros(n,k);      % stores row numbers of rows that are equal to the corresponding
                               %   row in uniqueRowMatrixA
counterA = 2;                  % stores current row in uniqueRowMatrixA to where a unique row should
                               %   be copied next
equalA = 1;                    % temp bool variable to show if two rows are equal up to current point
matchedRowA = 0;               % temp variable to hold row # of row in uniqueRowMatrixA that
                               %   matches with current row of A
                               
% first row in A is first unique row
for i = 1:n
    uniqueRowMatrixA(1,i) = A(1,i);
end
multVectorA(1) = 1;
multMatrixA(1,1) = 1;

for i = 2:n % for each remaining row in A

    for r = 1:(counterA-1) % compare to each row in uniqueRowMatrixA
        for s = 1:n        % compare to element in each column
            if (A(i,s) ~= uniqueRowMatrixA(r,s))
                equalA = 0;
                break
            end
        end
        if (equalA == 1)
            matchedRowA = r;
            break
        end
        equalA = 1;
    end
    if (matchedRowA ~= 0)
        multVectorA(r) = multVectorA(r) + 1;
        multMatrixA(r,multVectorA(r)) = i;
    else
        for j = 1:n
            uniqueRowMatrixA(counterA,j) = A(i,j);
        end
        multVectorA(counterA) = multVectorA(counterA) + 1;
        multMatrixA(counterA,1) = i;
        counterA = counterA + 1;
    end
    equalA = 1;
    matchedRowA = 0;
end

uniqueRowMatrixB = zeros(n,n); % stores a copy for each unique row in B
multVectorB = zeros(1,n);      % stores number of rows equal to corresponding row in uniqueRowMatrixB
multMatrixB = zeros(n,k);      % stores row numbers of rows that are equal to the corresponding
                               %   row in uniqueRowMatrixB
counterB = 2;                  % stores current row in uniqueRowMatrixB to where a unique row should
                               %   be copied next
equalB = 1;                    % temp bool variable to show if two rows are equal up to current point
matchedRowB = 0;               % temp variable to hold row # of row in uniqueRowMatrixB that
                               %   matches with current row of B
                               
% first row in B is first unique row
for i = 1:n
    uniqueRowMatrixB(1,i) = B(1,i);
end
multVectorB(1) = 1;
multMatrixB(1,1) = 1;

for i = 2:n % for each remaining row in B

    for r = 1:(counterB-1) % compare to each row in uniqueRowMatrixB
        for s = 1:n        % compare to element in each column
            if (B(i,s) ~= uniqueRowMatrixB(r,s))
                equalB = 0;
                break
            end
        end
        if (equalB == 1)
            matchedRowB = r;
            break
        end
        equalB = 1;
    end
    if (matchedRowB ~= 0)
        multVectorB(r) = multVectorB(r) + 1;
        multMatrixB(r,multVectorB(r)) = i;
    else
        for j = 1:n
            uniqueRowMatrixB(counterB,j) = B(i,j);
        end
        multVectorB(counterB) = multVectorB(counterB) + 1;
        multMatrixB(counterB,1) = i;
        counterB = counterB + 1;
    end
    equalB = 1;
    matchedRowB = 0;
end

Aprime = A';
Bprime = B';

uniqueRowMatrixAprime = zeros(n,n); % stores a copy for each unique row in Aprime
multVectorAprime = zeros(1,n);      % stores number of rows equal to corresponding row in uniqueRowMatrixAprime
multMatrixAprime = zeros(n,k);      % stores row numbers of rows that are equal to the corresponding
                                    %   row in uniqueRowMatrixAprime
counterAprime = 2;                  % stores current row in uniqueRowMatrixAprime to where a unique row should
                                    %   be copied next
equalAprime = 1;                    % temp bool variable to show if two rows are equal up to current point
matchedRowAprime = 0;               % temp variable to hold row # of row in uniqueRowMatrixAprime that
                                    %   matches with current row of Aprime
                               
% first row in Aprime is first unique row
for i = 1:n
    uniqueRowMatrixAprime(1,i) = Aprime(1,i);
end
multVectorAprime(1) = 1;
multMatrixAprime(1,1) = 1;

for i = 2:n % for each remaining row in Aprime

    for r = 1:(counterAprime-1) % compare to each row in uniqueRowMatrixAprime
        for s = 1:n             % compare to element in each column
            if (Aprime(i,s) ~= uniqueRowMatrixAprime(r,s))
                equalAprime = 0;
                break
            end
        end
        if (equalAprime == 1)
            matchedRowAprime = r;
            break
        end
        equalAprime = 1;
    end
    if (matchedRowAprime ~= 0)
        multVectorAprime(r) = multVectorAprime(r) + 1;
        multMatrixAprime(r,multVectorAprime(r)) = i;
    else
        for j = 1:n
            uniqueRowMatrixAprime(counterAprime,j) = Aprime(i,j);
        end
        multVectorAprime(counterAprime) = multVectorAprime(counterAprime) + 1;
        multMatrixAprime(counterAprime,1) = i;
        counterAprime = counterAprime + 1;
    end
    equalAprime = 1;
    matchedRowAprime = 0;
end

uniqueRowMatrixBprime = zeros(n,n); % stores a copy for each unique row in Bprime
multVectorBprime = zeros(1,n);      % stores number of rows equal to corresponding row in uniqueRowMatrixBprime
multMatrixBprime = zeros(n,k);      % stores row numbers of rows that are equal to the corresponding
                                    %   row in uniqueRowMatrixBprime
counterBprime = 2;                  % stores current row in uniqueRowMatrixBprime to where a unique row should
                                    %   be copied next
equalBprime = 1;                    % temp bool variable to show if two rows are equal up to current point
matchedRowBprime = 0;               % temp variable to hold row # of row in uniqueRowMatrixBprime that
                                    %   matches with current row of Bprime
                               
% first row in Bprime is first unique row
for i = 1:n
    uniqueRowMatrixBprime(1,i) = Bprime(1,i);
end
multVectorBprime(1) = 1;
multMatrixBprime(1,1) = 1;

for i = 2:n % for each remaining row in Bprime

    for r = 1:(counterBprime-1) % compare to each row in uniqueRowMatrixB
        for s = 1:n             % compare to element in each column
            if (Bprime(i,s) ~= uniqueRowMatrixBprime(r,s))
                equalBprime = 0;
                break
            end
        end
        if (equalBprime == 1)
            matchedRowBprime = r;
            break
        end
        equalBprime = 1;
    end
    if (matchedRowBprime ~= 0)
        multVectorBprime(r) = multVectorBprime(r) + 1;
        multMatrixBprime(r,multVectorBprime(r)) = i;
    else
        for j = 1:n
            uniqueRowMatrixBprime(counterBprime,j) = Bprime(i,j);
        end
        multVectorBprime(counterBprime) = multVectorBprime(counterBprime) + 1;
        multMatrixBprime(counterBprime,1) = i;
        counterBprime = counterBprime + 1;
    end
    equalBprime = 1;
    matchedRowBprime = 0;
end

sortMultiplicitiesA = sort(multVectorA);
sortMultiplicitiesB = sort(multVectorB);
sortMultiplicitiesAprime = sort(multVectorAprime);
sortMultiplicitiesBprime = sort(multVectorBprime);

A_Row_Mult_Equal_B_Row_Mult = 1;
Aprime_Row_Mult_Equal_Bprime_Row_Mult = 1;

for i = 1:n

    if (sortMultiplicitiesA(i) ~= sortMultiplicitiesB(i))
        A_Row_Mult_Equal_B_Row_Mult = 0;
    end
    
    if (sortMultiplicitiesAprime(i) ~= sortMultiplicitiesBprime(i))
        Aprime_Row_Mult_Equal_Bprime_Row_Mult = 0;
    end
    
end

cont = 0;

if (A_Row_Mult_Equal_B_Row_Mult == 1)
    % A row multiplicity equals B row multiplicity
    if (Aprime_Row_Mult_Equal_Bprime_Row_Mult == 1)
        % A transpose row multiplicity equals B transpose row multiplicity
        consideration = 'Considering if A is similar to B';
        cont = 1;
    end
end

if (cont == 0)
    reason_for_failure = 'Similarity impossible due to column and row multiplicities of A and B';
    similar = 0;
    return
end

% permute A and B into DRM form

newA = zeros(n,n);
newB = zeros(n,n);
count = 1;
count2 = 1;

for i = 0:(k-1)
    for j = 1:n
        if (multVectorA(j) == k - i)
            for f = 1:k
                if (multMatrixA(j,f) ~= 0)
                    listA(count) = multMatrixA(j,f);
                    count = count+1;
                end
            end
            newMultVectorA(count2) = k - i;
            newMultVectorA(count2+1) = 0;
            count2 = count2 + 1;
        end 
    end
end

count = 1;
count2 = 1;

for i = 0:(k-1)
    for j = 1:n
        if (multVectorB(j) == k - i)
            for f = 1:k
                if (multMatrixB(j,f) ~= 0)
                    listB(count) = multMatrixB(j,f);
                    count = count+1;
                end
            end
            newMultVectorB(count2) = k - i;
            newMultVectorB(count2+1) = 0;
            count2 = count2 + 1;
        end 
    end
end

I = eye(n,n);
P_A = I(listA,:);
P_B = I(listB,:);

newA = P_A * A * P_A';
newB = P_B * B * P_B';

if (newA == newB)
    % A and B are similar if initial permutation yields newA == newB
    permutation_matrix = P_A' * P_B;
    similar = 1;
    return
end

% finally, we consider all permutations within row multiplicities, and output
%   any permutation matrix that satisfies A = P*B*P'

numGroupsA4 = 0;
numGroupsA3 = 0;
numGroupsA2 = 0;
numGroupsA1 = 0;

i = 1;
f = 1;
count = 1;

while (newMultVectorA(i) ~= 0)
    if (newMultVectorA(i) == 4)
        for j = 1:newMultVectorA(i)
            groupA4(f,j) = listA(count);
            count = count+1;
        end
        f = f+1;
        numGroupsA4 = numGroupsA4 + 1;
    end
    i = i+1;
end

i = 1;
f = 1;

while (newMultVectorA(i) ~= 0)
    if (newMultVectorA(i) == 3)
        for j = 1:newMultVectorA(i)
            groupA3(f,j) = listA(count);
            count = count+1;
        end
        f = f+1;
        numGroupsA3 = numGroupsA3 + 1;
    end
    i = i+1;
end

i = 1;
f = 1;

while (newMultVectorA(i) ~= 0)
    if (newMultVectorA(i) == 2)
        for j = 1:newMultVectorA(i)
            groupA2(f,j) = listA(count);
            count = count+1;
        end
        f = f+1;
        numGroupsA2 = numGroupsA2 + 1;
    end
    i = i+1;
end

i = 1;
f = 1;

while (newMultVectorA(i) ~= 0)
    if (newMultVectorA(i) == 1)
        for j = 1:newMultVectorA(i)
            groupA1(f,j) = listA(count);
            count = count+1;
        end
        f = f+1;
        numGroupsA1 = numGroupsA1 + 1;
    end
    i = i+1;
end

totalGroupsA = numGroupsA1 + numGroupsA2 + numGroupsA3 + numGroupsA4;

groupA4a = zeros(1,4);
groupA4b = zeros(1,4);
groupA4c = zeros(1,4);
groupA4d = zeros(1,4);
groupA3a = zeros(1,3);
groupA3b = zeros(1,3);
groupA3c = zeros(1,3);
groupA3d = zeros(1,3);
groupA3e = zeros(1,3);
groupA2a = zeros(1,2);
groupA2b = zeros(1,2);
groupA2c = zeros(1,2);
groupA2d = zeros(1,2);
groupA2e = zeros(1,2);
groupA2f = zeros(1,2);
groupA2g = zeros(1,2);
groupA2h = zeros(1,2);
groupA1a = zeros(1,numGroupsA1);

if (numGroupsA4 == 4)
    for i = 1:4
        groupA4a(i) = groupA4(1,i);
        groupA4b(i) = groupA4(2,i);
        groupA4c(i) = groupA4(3,i);
        groupA4d(i) = groupA4(4,i);
    end
    permsA4(:,:,1) = perms( [ groupA4a ] );
    permsA4(:,:,2) = perms( [ groupA4b ] );
    permsA4(:,:,3) = perms( [ groupA4c ] );
    permsA4(:,:,4) = perms( [ groupA4d ] );
end
if (numGroupsA4 == 3)
    for i = 1:4
        groupA4a(i) = groupA4(1,i);
        groupA4b(i) = groupA4(2,i);
        groupA4c(i) = groupA4(3,i);
    end
    permsA4(:,:,1) = perms( [ groupA4a ] );
    permsA4(:,:,2) = perms( [ groupA4b ] );
    permsA4(:,:,3) = perms( [ groupA4c ] );
end
if (numGroupsA4 == 2)
    for i = 1:4
        groupA4a(i) = groupA4(1,i);
        groupA4b(i) = groupA4(2,i);
    end
    permsA4(:,:,1) = perms( [ groupA4a ] );
    permsA4(:,:,2) = perms( [ groupA4b ] );
end
if (numGroupsA4 == 1)
    permsA4 = perms( [ groupA4 ] );
end
if (numGroupsA3 == 5)
    for i = 1:3
        groupA3a(i) = groupA3(1,i);
        groupA3b(i) = groupA3(2,i);
        groupA3c(i) = groupA3(3,i);
        groupA3d(i) = groupA3(4,i);
        groupA3e(i) = groupA3(5,i);
    end
    permsA3(:,:,1) = perms( [ groupA3a ] );
    permsA3(:,:,2) = perms( [ groupA3b ] );
    permsA3(:,:,3) = perms( [ groupA3c ] );
    permsA3(:,:,4) = perms( [ groupA3d ] );
    permsA3(:,:,5) = perms( [ groupA3e ] );
end
if (numGroupsA3 == 4)
    for i = 1:3
        groupA3a(i) = groupA3(1,i);
        groupA3b(i) = groupA3(2,i);
        groupA3c(i) = groupA3(3,i);
        groupA3d(i) = groupA3(4,i);
    end
    permsA3(:,:,1) = perms( [ groupA3a ] );
    permsA3(:,:,2) = perms( [ groupA3b ] );
    permsA3(:,:,3) = perms( [ groupA3c ] );
    permsA3(:,:,4) = perms( [ groupA3d ] );
end
if (numGroupsA3 == 3)
    for i = 1:3
        groupA3a(i) = groupA3(1,i);
        groupA3b(i) = groupA3(2,i);
        groupA3c(i) = groupA3(3,i);
    end
    permsA3(:,:,1) = perms( [ groupA3a ] );
    permsA3(:,:,2) = perms( [ groupA3b ] );
    permsA3(:,:,3) = perms( [ groupA3c ] );
end
if (numGroupsA3 == 2)
    for i = 1:3
        groupA3a(i) = groupA3(1,i);
        groupA3b(i) = groupA3(2,i);
    end
    permsA3(:,:,1) = perms( [ groupA3a ] );
    permsA3(:,:,2) = perms( [ groupA3b ] );
end
if (numGroupsA3 == 1)
    permsA3 = perms( [ groupA3 ] );
end
if (numGroupsA2 == 8)
    for i = 1:2
        groupA2a(i) = groupA2(1,i);
        groupA2b(i) = groupA2(2,i);
        groupA2c(i) = groupA2(3,i);
        groupA2d(i) = groupA2(4,i);
        groupA2e(i) = groupA2(5,i);
        groupA2f(i) = groupA2(6,i);
        groupA2g(i) = groupA2(7,i);
        groupA2h(i) = groupA2(8,i);
    end
    permsA2(:,:,1) = perms( [ groupA2a ] );
    permsA2(:,:,2) = perms( [ groupA2b ] );
    permsA2(:,:,3) = perms( [ groupA2c ] );
    permsA2(:,:,4) = perms( [ groupA2d ] );
    permsA2(:,:,5) = perms( [ groupA2e ] );
    permsA2(:,:,6) = perms( [ groupA2f ] );
    permsA2(:,:,7) = perms( [ groupA2g ] );
    permsA2(:,:,8) = perms( [ groupA2h ] );
end
if (numGroupsA2 == 7)
    for i = 1:2
        groupA2a(i) = groupA2(1,i);
        groupA2b(i) = groupA2(2,i);
        groupA2c(i) = groupA2(3,i);
        groupA2d(i) = groupA2(4,i);
        groupA2e(i) = groupA2(5,i);
        groupA2f(i) = groupA2(6,i);
        groupA2g(i) = groupA2(7,i);
    end
    permsA2(:,:,1) = perms( [ groupA2a ] );
    permsA2(:,:,2) = perms( [ groupA2b ] );
    permsA2(:,:,3) = perms( [ groupA2c ] );
    permsA2(:,:,4) = perms( [ groupA2d ] );
    permsA2(:,:,5) = perms( [ groupA2e ] );
    permsA2(:,:,6) = perms( [ groupA2f ] );
    permsA2(:,:,7) = perms( [ groupA2g ] );
end
if (numGroupsA2 == 6)
    for i = 1:2
        groupA2a(i) = groupA2(1,i);
        groupA2b(i) = groupA2(2,i);
        groupA2c(i) = groupA2(3,i);
        groupA2d(i) = groupA2(4,i);
        groupA2e(i) = groupA2(5,i);
        groupA2f(i) = groupA2(6,i);
    end
    permsA2(:,:,1) = perms( [ groupA2a ] );
    permsA2(:,:,2) = perms( [ groupA2b ] );
    permsA2(:,:,3) = perms( [ groupA2c ] );
    permsA2(:,:,4) = perms( [ groupA2d ] );
    permsA2(:,:,5) = perms( [ groupA2e ] );
    permsA2(:,:,6) = perms( [ groupA2f ] );
end
if (numGroupsA2 == 5)
    for i = 1:2
        groupA2a(i) = groupA2(1,i);
        groupA2b(i) = groupA2(2,i);
        groupA2c(i) = groupA2(3,i);
        groupA2d(i) = groupA2(4,i);
        groupA2e(i) = groupA2(5,i);
    end
    permsA2(:,:,1) = perms( [ groupA2a ] );
    permsA2(:,:,2) = perms( [ groupA2b ] );
    permsA2(:,:,3) = perms( [ groupA2c ] );
    permsA2(:,:,4) = perms( [ groupA2d ] );
    permsA2(:,:,5) = perms( [ groupA2e ] );
end
if (numGroupsA2 == 4)
    for i = 1:2
        groupA2a(i) = groupA2(1,i);
        groupA2b(i) = groupA2(2,i);
        groupA2c(i) = groupA2(3,i);
        groupA2d(i) = groupA2(4,i);
    end
    permsA2(:,:,1) = perms( [ groupA2a ] );
    permsA2(:,:,2) = perms( [ groupA2b ] );
    permsA2(:,:,3) = perms( [ groupA2c ] );
    permsA2(:,:,4) = perms( [ groupA2d ] );
end
if (numGroupsA2 == 3)
    for i = 1:2
        groupA2a(i) = groupA2(1,i);
        groupA2b(i) = groupA2(2,i);
        groupA2c(i) = groupA2(3,i);
    end
    permsA2(:,:,1) = perms( [ groupA2a ] );
    permsA2(:,:,2) = perms( [ groupA2b ] );
    permsA2(:,:,3) = perms( [ groupA2c ] );
end
if (numGroupsA2 == 2)
    for i = 1:2
        groupA2a(i) = groupA2(1,i);
        groupA2b(i) = groupA2(2,i);
    end
    permsA2(:,:,1) = perms( [ groupA2a ] );
    permsA2(:,:,2) = perms( [ groupA2b ] );
end
if (numGroupsA2 == 1)
    permsA2 = perms( [ groupA2 ] );
end
for i = 1:numGroupsA1
    groupA1a(i) = groupA1(i,1);
end
permsA1 = perms( [ groupA1a ] );

% consider row multiplicity vector and test for final permutation matrix

newMultVectorA_temp = zeros(1,n);
for i = 1:totalGroupsA
    newMultVectorA_temp(i) = newMultVectorA(i);
end
newMultVectorA = newMultVectorA_temp;

count = 0;

if (newMultVectorA == [ 4 4 3 3 1 1 0 0 0 0 0 0 0 0 0 0 ])
    for a = 1:factorial(4)
        for b = 1:factorial(4)
            for c = 1:factorial(3)
                for d = 1:factorial(3)
                    for e = 1:factorial(2)
                        for f = 1:factorial(2)
                            for g = 1:factorial(2)
                                temp_list = [ permsA4(a,:,e) permsA4(b,:,3-e) permsA3(c,:,f) permsA3(d,:,3-f) permsA1(g,:) ];
                                temp_perm = I(temp_list,:);
                                temp_newA = temp_perm * A * temp_perm';
                                if (temp_newA == newB)
                                    similar = 1;
                                    permutation_matrix = P_B' * temp_perm;
                                    return
                                end
                                count = count + 1;
                            end
                        end
                    end
                end
            end
        end
    end
end

if (newMultVectorA == [ 4 4 2 2 1 1 1 1 0 0 0 0 0 0 0 0 ])
    for a = 1:factorial(4)
        for b = 1:factorial(4)
            for c = 1:factorial(2)
                for d = 1:factorial(2)
                    for e = 1:factorial(2)
                        for f = 1:factorial(2)
                            for g = 1:factorial(4)
                                temp_list = [ permsA4(a,:,c) permsA4(b,:,3-c) permsA2(d,:,f) permsA2(e,:,3-f) permsA1(g,:) ];
                                temp_perm = I(temp_list,:);
                                temp_newA = temp_perm * A * temp_perm';
                                if (temp_newA == newB)
                                    similar = 1;
                                    permutation_matrix = P_B' * temp_perm;
                                    return
                                end
                                count = count + 1;
                            end
                        end
                    end
                end
            end
        end
    end
end

if (newMultVectorA == [ 4 3 3 3 1 1 1 0 0 0 0 0 0 0 0 0 ])
    orderA3 = perms( [ 1 2 3 ] );
    for a = 1:factorial(4)
        for b = 1:factorial(3)
            for c = 1:factorial(3)
                for d = 1:factorial(3)
                    for e = 1:factorial(3)
                        for f = 1:factorial(3)
                            temp_order = orderA3(e,:);
                            one = temp_order(1);
                            two = temp_order(2);
                            three = temp_order(3);
                            temp_list = [ permsA4(a,:) permsA3(b,:,one) permsA3(c,:,two) permsA3(d,:,three) permsA1(f,:) ];
                            temp_perm = I(temp_list,:);
                            temp_newA = temp_perm * A * temp_perm';
                            if (temp_newA == newB)
                                similar = 1;
                                permutation_matrix = P_B' * temp_perm;
                                return
                            end
                            count = count + 1;
                        end
                    end
                end
            end
        end
    end
end

if (newMultVectorA == [ 4 3 3 2 1 1 1 1 0 0 0 0 0 0 0 0 ])
    for a = 1:factorial(4)
        for b = 1:factorial(3)
            for c = 1:factorial(3)
                for d = 1:factorial(2)
                    for e = 1:factorial(2)
                        for f = 1:factorial(4)
                            temp_list = [ permsA4(a,:) permsA3(b,:,d) permsA3(c,:,3-d) permsA2(e,:) permsA1(f,:) ];
                            temp_perm = I(temp_list,:);
                            temp_newA = temp_perm * A * temp_perm';
                            if (temp_newA == newB)
                                similar = 1;
                                permutation_matrix = P_B' * temp_perm;
                                return
                            end
                            count = count + 1;
                        end
                    end
                end
            end
        end
    end
end

if (newMultVectorA == [ 4 4 2 2 2 2 0 0 0 0 0 0 0 0 0 0 ])
    orderA2 = perms( [ 1 2 3 4 ] );
    for a = 1:factorial(4)
        for b = 1:factorial(4)
            for c = 1:factorial(2)
                for d = 1:factorial(2)
                    for e = 1:factorial(2)
                        for f = 1:factorial(2)
                            for g = 1:factorial(2)
                                for h = 1:factorial(4)
                                    temp_order = orderA2(h,:);
                                    one = temp_order(1);
                                    two = temp_order(2);
                                    three = temp_order(3);
                                    four = temp_order(4);
                                    temp_list = [ permsA4(a,:,c) permsA4(b,:,3-c) permsA2(d,:,one) permsA2(e,:,two) permsA2(f,:,three) permsA2(g,:,four) ];
                                    temp_perm = I(temp_list,:);
                                    temp_newA = temp_perm * A * temp_perm';
                                    if (temp_newA == newB)
                                        similar = 1;
                                        permutation_matrix = P_B' * temp_perm;
                                        return
                                    end
                                    count = count + 1;
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

if (newMultVectorA == [ 4 3 2 2 2 1 1 1 0 0 0 0 0 0 0 0 ])
    orderA2 = perms( [ 1 2 3 ] );
    for a = 1:factorial(4)
        for b = 1:factorial(3)
            for c = 1:factorial(2)
                for d = 1:factorial(2)
                    for e = 1:factorial(2)
                        for f = 1:factorial(3)
                            for g = 1:factorial(3)
                                temp_order = orderA2(f,:);
                                one = temp_order(1);
                                two = temp_order(2);
                                three = temp_order(3);
                                temp_list = [ permsA4(a,:) permsA3(b,:) permsA2(c,:,one) permsA2(d,:,two) permsA2(e,:,three) permsA1(g,:) ];
                                temp_perm = I(temp_list,:);
                                temp_newA = temp_perm * A * temp_perm';
                                if (temp_newA == newB)
                                    similar = 1;
                                    permutation_matrix = P_B' * temp_perm;
                                    return
                                end
                                count = count + 1;
                            end
                        end
                    end
                end
            end
        end
    end
end

if (newMultVectorA == [ 4 2 2 2 1 1 1 1 1 1 0 0 0 0 0 0 ])
    orderA2 = perms( [ 1 2 3 ] );
    for a = 1:factorial(4)
        for b = 1:factorial(2)
            for c = 1:factorial(2)
                for d = 1:factorial(2)
                    for e = 1:factorial(3)
                        for f = 1:factorial(6)
                            temp_order = orderA2(e,:);
                            one = temp_order(1);
                            two = temp_order(2);
                            three = temp_order(3);
                            temp_list = [ permsA4(a,:) permsA2(b,:,one) permsA2(c,:,two) permsA2(d,:,three) permsA1(f,:) ];
                            temp_perm = I(temp_list,:);
                            temp_newA = temp_perm * A * temp_perm';
                            if (temp_newA == newB)
                                similar = 1;
                                permutation_matrix = P_B' * temp_perm;
                                return
                            end
                            count = count + 1;
                        end
                    end
                end
            end
        end
    end
end

if (newMultVectorA == [ 3 3 3 3 1 1 1 1 0 0 0 0 0 0 0 0 ])
    orderA3 = perms( [ 1 2 3 4 ] );
    for a = 1:factorial(3)
        for b = 1:factorial(3)
            for c = 1:factorial(3)
                for d = 1:factorial(3)
                    for e = 1:factorial(4)
                        for f = 1:factorial(4)
                            temp_order = orderA3(e,:);
                            one = temp_order(1);
                            two = temp_order(2);
                            three = temp_order(3);
                            four = temp_order(4);
                            temp_list = [ permsA3(a,:,one) permsA3(b,:,two) permsA3(c,:,three)  permsA3(d,:,four) permsA1(f,:) ];
                            temp_perm = I(temp_list,:);
                            temp_newA = temp_perm * A * temp_perm';
                            if (temp_newA == newB)
                                similar = 1;
                                permutation_matrix = P_B' * temp_perm;
                                return
                            end
                            count = count + 1;
                        end
                    end
                end
            end
        end
    end
end

if (newMultVectorA == [ 3 3 3 2 1 1 1 1 1 0 0 0 0 0 0 0 ])
    orderA3 = perms( [ 1 2 3 ] );
    for a = 1:factorial(3)
        for b = 1:factorial(3)
            for c = 1:factorial(3)
                for d = 1:factorial(3)
                    for e = 1:factorial(2)
                        for f = 1:factorial(5)
                            temp_order = orderA3(d,:);
                            one = temp_order(1);
                            two = temp_order(2);
                            three = temp_order(3);
                            temp_list = [ permsA3(a,:,one) permsA3(b,:,two) permsA3(c,:,three) permsA2(e,:) permsA1(f,:) ];
                            temp_perm = I(temp_list,:);
                            temp_newA = temp_perm * A * temp_perm';
                            if (temp_newA == newB)
                                similar = 1;
                                permutation_matrix = P_B' * temp_perm;
                                return
                            end
                            count = count + 1;
                        end
                    end
                end
            end
        end
    end
end

if (newMultVectorA == [ 4 3 2 2 1 1 1 1 1 0 0 0 0 0 0 0 ])
    for a = 1:factorial(4)
        for b = 1:factorial(3)
            for c = 1:factorial(2)
                for d = 1:factorial(2)
                    for e = 1:factorial(2)
                        for f = 1:factorial(5)
                            temp_list = [ permsA4(a,:) permsA3(b,:) permsA2(c,:,e) permsA2(d,:,3-e) permsA1(f,:) ];
                            temp_perm = I(temp_list,:);
                            temp_newA = temp_perm * A * temp_perm';
                            if (temp_newA == newB)
                                similar = 1;
                                permutation_matrix = P_B' * temp_perm;
                                return
                            end
                            count = count + 1;
                        end
                    end
                end
            end
        end
    end
end

if (newMultVectorA == [ 4 3 2 1 1 1 1 1 1 1 0 0 0 0 0 0 ])
    for a = 1:factorial(4)
        for b = 1:factorial(3)
            for c = 1:factorial(2)
                for d = 1:factorial(7)
                    temp_list = [ permsA4(a,:) permsA3(b,:) permsA2(c,:) permsA1(d,:) ];
                    temp_perm = I(temp_list,:);
                    temp_newA = temp_perm * A * temp_perm';
                    if (temp_newA == newB)
                        similar = 1;
                        permutation_matrix = P_B' * temp_perm;
                        return
                    end
                    count = count + 1;
                end
            end
        end
    end
end

if (newMultVectorA == [ 3 3 2 2 2 2 1 1 0 0 0 0 0 0 0 0 ])
    orderA2 = perms( [ 1 2 3 4 ] );
    for a = 1:factorial(3)
        for b = 1:factorial(3)
            for c = 1:factorial(2)
                for d = 1:factorial(2)
                    for e = 1:factorial(2)
                        for f = 1:factorial(2)
                            for g = 1:factorial(2)
                                for h = 1:factorial(4)
                                    for i = 1:factorial(2)
                                        temp_order = orderA2(h,:);
                                        one = temp_order(1);
                                        two = temp_order(2);
                                        three = temp_order(3);
                                        four = temp_order(4);
                                        temp_list = [ permsA3(a,:,c) permsA3(b,:,3-c) permsA2(d,:,one) permsA2(e,:,two) permsA2(f,:,three) permsA2(g,:,four) permsA1(i,:) ];
                                        temp_perm = I(temp_list,:);
                                        temp_newA = temp_perm * A * temp_perm';
                                        if (temp_newA == newB)
                                            similar = 1;
                                            permutation_matrix = P_B' * temp_perm;
                                            return
                                        end
                                        count = count + 1;
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

if (newMultVectorA == [ 3 2 2 2 1 1 1 1 1 1 1 0 0 0 0 0 ])
    orderA2 = perms( [ 1 2 3 ] );
    for a = 1:factorial(3)
        for b = 1:factorial(2)
            for c = 1:factorial(2)
                for d = 1:factorial(2)
                    for e = 1:factorial(3)
                        for f = 1:factorial(7)
                            temp_order = orderA2(e,:);
                            one = temp_order(1);
                            two = temp_order(2);
                            three = temp_order(3);
                            temp_list = [ permsA3(a,:) permsA2(b,:,one) permsA2(c,:,two) permsA2(d,:,three) permsA1(f,:) ];
                            temp_perm = I(temp_list,:);
                            temp_newA = temp_perm * A * temp_perm';
                            if (temp_newA == newB)
                                similar = 1;
                                permutation_matrix = P_B' * temp_perm;
                                return
                            end
                            count = count + 1;
                        end
                    end
                end
            end
        end
    end
end

if (newMultVectorA == [ 2 2 2 2 2 1 1 1 1 1 1 0 0 0 0 0 ])
    orderA2 = perms( [ 1 2 3 4 5 ] );
    for a = 1:factorial(2)
        for b = 1:factorial(2)
            for c = 1:factorial(2)
                for d = 1:factorial(2)
                    for e = 1:factorial(2)
                        for f = 1:factorial(5)
                            for g = 1:factorial(6)
                                temp_order = orderA2(f,:);
                                one = temp_order(1);
                                two = temp_order(2);
                                three = temp_order(3);
                                four = temp_order(4);
                                five = temp_order(5);
                                temp_list = [ permsA2(a,:,one) permsA2(b,:,two) permsA2(c,:,three) permsA2(d,:,four) permsA2(e,:,five) permsA1(g,:) ];
                                temp_perm = I(temp_list,:);
                                temp_newA = temp_perm * A * temp_perm';
                                if (temp_newA == newB)
                                    similar = 1;
                                    permutation_matrix = P_B' * temp_perm;
                                    return
                                end
                                count = count + 1;
                            end
                        end
                    end
                end
            end
        end
    end
end

if (newMultVectorA == [ 3 3 3 1 1 1 1 1 1 1 0 0 0 0 0 0 ])
    orderA3 = perms( [ 1 2 3 ] );
    for a = 1:factorial(3)
        for b = 1:factorial(3)
            for c = 1:factorial(3)
                for d = 1:factorial(3)
                    for e = 1:factorial(7)
                        temp_order = orderA3(d,:);
                        one = temp_order(1);
                        two = temp_order(2);
                        three = temp_order(3);
                        temp_list = [ permsA3(a,:,one) permsA3(b,:,two) permsA3(c,:,three) permsA1(e,:) ];
                        temp_perm = I(temp_list,:);
                        temp_newA = temp_perm * A * temp_perm';
                        if (temp_newA == newB)
                            similar = 1;
                            permutation_matrix = P_B' * temp_perm;
                            return
                        end
                        count = count + 1;
                    end
                end
            end
        end
    end
end

if (count == 0)
    reason_for_failure = 'Unknown row multiplicity vector'
    newMultVectorA
    similar = -1;
    return
end

reason_for_failure = 'Matrices are not similar';

similar = 0;