# Approximates Accuracies for multiple traits given single-trait accuracies;
# works for any 3 trait combination

# Programmed in Octave by Ignacy Misztal, 1998
# Example from paper by Strabel et al. in J. Animal Science 79:833-839.

# Pro

g0=[5 2 1
2 3 1
1 1 6];

r0=[8 3 2
3 5 4
2 4 9];

#r0=[8 0 0
#0 5 0
#0 0 9];
#g0=[5 0 0
#0 3 0
#0 0 6];

# EON from records only
eonr=[.8 .5 0];

#EON from pedigrees only
eonp=[2.3 3.4 1.2];


# Calculate single-trait variance ratios for conversion from eonp to epnp
alpha_a=(diag(r0) ./ diag(g0))';
alpha_s=((diag(r0)+.75*diag(g0))  ./ (.25*diag(g0)))';

# convert eonp to epnp
epnp=((alpha_a+eonp) .* alpha_s) ./ alpha_a - alpha_s ;


# MT contribution due to records
o12=sqrt(diag(eonr));
p=o12 !=0;
w1=o12*pinv(p*r0*p)*o12

# MT contribution to the animal itself
w2=pinv(g0)


# MT contribution due to pedigrees

nt=size(epnp,2);
w3=zeros(nt);

for i=1:nt
  p=zeros(nt);min=1e10;
  for j=1:nt
    if (epnp(j) !=0) 
       p(j,j)=1;
       if (epnp(j) <min) 
         min=epnp(j);
       endif
    endif
  endfor
  r00=pinv(p*r0*p);
  #w3=w3+min*pinv(g0)/3*(eye(nt)-pinv(eye(nt)+.75*r00*g0))
  #w3=w3+min*(eye(nt)-2/3*pinv(4.0/3*pinv(g0)+r00)*pinv(g0))
  w3=w3+min*(eye(nt)-pinv(pinv(g0)+0.75*r00)*pinv(g0))
  pinv(pinv(g0)+0.75*r00)
  epnp=max(0,epnp-min);
endfor 

w3=pinv(g0)/3*w3


w=w1+w2+w3 
acc=diag(eye(3)-diag(diag(inv(w)))*inv(diag(diag(g0))))