Fixed additive relationship matrix creation

MATLAB/beefblup.m

@@ -61,13 +61,29 @@ data = sortrows(data,2);
 strings = [1 3 4];
 data(:,strings) = cellfun(@num2str, data(:,strings), 'UniformOutput', false);
 
-% Create a lookup lambda function to find the animal represented by a
+% Define fields to hold id values for animals and their parents
-% certain id
+ids = char(data{:,1});
-animal2row = @(id) find(strcmp(data(:,1), id));
+damids = char(data{:,3});
-row2animal = @(rownum) [data{rownum,1}];
+sireids = char(data{:,4});
-ids = [data{:,1}];
 numanimals = length(data(:,1));
 
+% Define fields to hold the index values for animals and their parents
+dam = zeros(numanimals,1);
+sire = zeros(numanimals,1);
+
+% Find all row numbers where an animal was a parent
+for i=1:numanimals
+   % Find all animals that this animal birthed
+   dammatch = ismember(damids, ids(i,:), 'rows');
+   damindexes = find(dammatch == 1);
+   dam(damindexes) = i;
+   
+   % Find all animals that this animal sired
+   sirematch = ismember(sireids, ids(i,:), 'rows');
+   sireindexes = find(sirematch == 1);
+   sire(sireindexes) = i;
+end
+
 % Store column numbers that need to be deleted
 % Column 6 contains an intermediate Excel calculation and always needs to
 % be deleted
@@ -181,53 +197,34 @@ disp('Creating the additive relationship matrix...')
 % Create an empty matrix for the additive relationship matrix
 A = zeros(numanimals, numanimals);
 
-% Create lambdas to find sire and dam of each animal
-id2dam = @(id) [data{animal2row(num2str(id)), 3}];
-id2sire = @(id) [data{animal2row(num2str(id)), 4}];
-row2dam = @(rownum) [data{rownum, 3}];
-row2sire = @(rownum) [data{rownum, 4}];
-rowofdam = @(id) animal2row(id2dam(id));
-rowofsire = @(id) animal2row(id2sire(id));
-row2rowofdam = @(rownum) rowofdam(row2animal(rownum));
-row2rowofsire = @(rownum) rowofsire(row2animal(rownum));
-
 % Create the additive relationship matrix by the FORTRAN method presented
 % by Henderson
 for i = 1:numanimals
-   id = row2animal(i);
+    if dam(i) ~= 0 && sire(i) ~= 0
-   if isempty(row2rowofsire(i)) && isempty(row2rowofdam(i))
+       for j = 1:(i-1)
+          A(j,i) = 0.5*(A(j,sire(i))+A(j,dam(i)));
+          A(i,j) = A(j,i);
+       end
+       A(i,i) = 1 + 0.5*A(sire(i),dam(i));
+    elseif dam(i) ~= 0 && sire(i) == 0
+        for j = 1:(i-1)
+           A(j,i) = 0.5*A(j,dam(i));
+           A(i,j) = A(j,i);
+        end
+        A(i,i) = 1;
+    elseif dam(i) == 0 && sire(i) ~=0
+        for j = 1:(i-1)
+           A(j,i) = 0.5*A(j,sire(i));
+           A(i,j) = A(j,i);
+        end
+        A(i,i) = 1;            
+    else
         for j = 1:(i-1)
            A(j,i) = 0;
+           A(i,j) = 0;
         end
         A(i,i) = 1;
-       continue
     end
-   
-   if isempty(row2rowofsire(i)) && not(isempty(row2dam(i)))
-       for j = 1:(i-1)
-           A(j,i) = A(j,row2rowofdam(i));
-           A(i,j) = A(j,row2rowofdam(i));
-       end
-       A(i,i) = 1;
-       continue
-   end
-   
-   if isempty(row2rowofdam(i)) && not(isempty(row2rowofsire(i)))
-       for j=1:(i-1)
-           A(j,i) = A(j,row2rowofsire(i));
-           A(i,j) = A(j,row2rowofsire(i));
-       end
-       A(i,i) = 1;
-       continue
-   end
-   
-   for j=1:(i-1)
-       A(j,i) = 0.5.*(A(j,row2rowofsire(i)) + A(j,row2rowofdam(i)));
-       A(i,j) = 0.5.*(A(j,row2rowofsire(i)) + A(j,row2rowofdam(i)));
-   end
-   A(i,i) = 1 + 0.5.*A(row2rowofsire(i), row2rowofdam(i));
-   continue
-     
 end
 
 disp('Creating the additive relationship matrix... Done!')

MATLAB/uniquenan.m

@@ -1,10 +0,0 @@
-% uniquenan
-% Serves the same purpose as UNIQUE, but ensures any NaN fields are not
-% counted
-function y = uniquenan(x)
-y = unique(x);
-if any(isnan(y))
-    y(isnan(y)) = [];
-    y(end + 1) = NaN;
-end
-end