basis = 3
dim = 4
index = 0

x = ones(1,36)
x(1:3) = 12
x(4:6) = 13
x(7:9) = 14
x(10:12) = 21
x(13:15) = 23
x(16:18) = 24
x(19:21) = 31
x(22:24) = 32
x(25:27) = 34
x(28:30) = 41
x(31:33) = 42
x(34:36) = 43

absum = zeros(1,6)
leap = index+1;
for j = 1:dim-1
  for k = j+1:dim
      tmp = 0;
      for i = 1: basis
         x(index + (j-1)*basis*(dim-1) + basis*(k-2)+i )
         tmp = tmp + x(index + (j-1)*basis*(dim-1) + basis*(k-2)+i ); 
      end
      for i = 1:basis
         x(index + (k-1)*basis*(dim-1) + basis*(j-1)+i )
         tmp = tmp + x(index + (k-1)*basis*(dim-1) + basis*(j-1)+i );
      end
      absum(leap) = tmp;
      leap = leap + 1
 end
end

% part(II) of absum(var_jk)
index = 0;
absum = zeros(1,6)
x = ones(1,18);
x(1:3) = 12;
x(4:6) = 13;
x(7:9) = 14;
x(10:12) = 23;
x(13:15) = 24;
x(16:18) = 34;

leap = 1;
for j = 1:dim-1
  for k = j+1:dim
      tmp = 0;
      for i = 1:basis
        x(index + basis*(2*dim-j)*(j-1)/2 + basis*(k-j-1)+i)
        tmp = tmp + abs(x(index + basis*(2*dim-j)*(j-1)/2 + basis*(k-j-1)+i));
      end
      absum(leap) = tmp;
      leap = leap + 1
  end
end


