******************************************************
*                    INDUPDAT                        *
*        Selection Index Updating for Maximum        *
*     Economic Return Using the EM Type Algorithm    *
*                       By                           * 
*                   Shizhong Xu                      *
*                Rutgers University                  *
*                  July 4, 1992                      *
******************************************************
;
OPTION LINESIZE=79;
PROC IML;
* P IS AN NXN PHENOTYPIC VARIANCE-COVARIANCE MATRIX;
P={
10816.000  13344.240 148.928  -.593  66.56,
13344.240  19881     252.39   -.932  112.8,
148.928    252.39    12.816    .010   4.296,
-.593      -.932     .010      .003   0,
66.56      112.8     4.296     0      64};
* G IS AN NXK GENETIC COVARIANCE MATRIX;
G={
3244.8  4133.961 50.982 -.459  166.787 -8.825  -27.798,
4133.961 6958.35 89.59  -.726  220.981 -23.261 -58.153,
50.982   89.59   3.204  -.014  -1.622  .555    1.248,
-.459    -.726   -.014  .002  .155     0      0,
-8.825   -23.261  .555    0    -2.16    9.6    14.04};
* W IS A KX1 VECTOR OF ECONOMIC WEIGHTS;
W={0, .0011, .0148, 0, 0, .5, .5};
* TOTAL PROPORTION SELECTED ;
PT=1/15;
* N IS AN 1XM VECTOR OF NUMBER OF TRAITS WITHIN STAGES;
N={2 2 1};
* MAXIMIZE GAIN OR RATIO OR PROFIT?;
MAX='GAIN';
* VCOST IS AN NX1 VECTOR OF COST PER TRAIT;
VCOST={.15, .150 ,.57, 1.15,5.4};
KAPA=50;
BETA=1;

*PROGRAM STARTS;

m=ncol(n);
c=n`;
ii=n[1]; c[1]=sum(vcost[1:ii]);
do jj=2 to m;
ii=ii+n[jj];
c[jj]=sum(vcost[ii-n[jj]+1:ii]);
end;
PI=3.1415926;
X=PROBIT(1-PT); T=EXP(-X*X/2)/SQRT(2*PI);I=T/PT;
PRINT 'VARIANCE COVARIANCE MATRICES AND ECONOMIC 
WEIGHTS';
PRINT P; PRINT G; PRINT W c;
PRINT 'TOTAL PROPORTION AND SELECTION INTENSITY';
PRINT PT I;
B=P**-1*G*W;
DG=G`*B*I/SQRT(B`*P*B); DH=W`*DG;
PRINT 'ONE STAGE INDEX SELECTION';
E=PT*BETA*KAPA;
COST=C[+];
PROFIT=E*DH-COST;
RATIO=DH/COST;
START MOD;
IF MAX='RATIO' THEN GOTO R;
IF MAX='PROFIT' THEN GOTO P;
PRINT B DG DH;GOTO G;
R:PRINT B DG DH COST RATIO;
GOTO G;
P:PRINT B DG DH COST PROFIT;

G:I1=I;
B=P[1:N[1],1:N[1]]**-1*G[1:N[1],]*W; I=N[1];
DO J=2 TO M;
I=I+N[J]  ;
QII=P[1:I,1:I]; QIJ=P[1:I-N[J],1:I]; AI=G[1:I,];
BI=(I(I)-QII**-1*QIJ`*B*(B`*QIJ*QII**-1*QIJ`*B)**-1*B`*QIJ)
   *QII**-1*AI*W;
B=(B//J(N[J],J-1,0))||BI;
END; VI=B`*P*B;
B=B*SQRT(DIAG(VI))**-1;
A=W`*G`*B;
PRINT 'MULTI-STAGE SELECTION';
PRINT B;

START SC(K,Q,C,SC); SC=0; CQ=1;
DO I=1 TO K;
SC=SC+C[I]*CQ;
CQ=CQ*Q[I];
END;
FINISH SC;

start dz(q,dz);
u=probit(j(nrow(q),1,1)-q);
dz=exp(-.5*u##2)/sqrt(2*3.1415926)/q;
finish dz;

q=j(m,1,pt##(1/m));
call sc(m,q,c,cost);
call dz(q,dz);
IF MAX='RATIO' THEN GOTO R1;
IF MAX='PROFIT' THEN GOTO P1;
profit=a*dz;GOTO G1;
R1:profit=a*dz/cost;
GOTO G1;
P1:profit=e*a*dz-cost;

g1:iter=0;
a:iter=iter+1;
do i=1 to m-1;
do j=i+1 to m;
r=q[i]*q[j];
qi=q[i];qj=q[j];
call sc(m,q,c,cost);
call dz(q,dz);
IF MAX='RATIO' THEN GOTO R2;
IF MAX='PROFIT' THEN GOTO P2;
prf0=a*dz;GOTO G2;
R2:prf0=a*dz/cost;
GOTO G2;
P2:prf0=e*a*dz-cost;

g2:vb=.9999;
va=r+.0001;
b:d=(vb-va)/5;
if d<.0001 then goto c;
do f=va to vb by d;
q[i]=f; q[j]=r/q[i];
call sc(m,q,c,cost);
call dz(q,dz);
IF MAX='RATIO' THEN GOTO R3;
IF MAX='PROFIT' THEN GOTO P3;
prf1=a*dz;GOTO G3;
R3:prf1=a*dz/cost;
GOTO G3;
P3:prf1=e*a*dz-cost;

g3:if prf1>prf0 then do;
prf0=prf1; qi=q[i]; qj=q[j];
end;end;
if (qi-d)>va then va=qi-d;
if (qi+d)<vb then vb=qi+d;
goto b;
c:q[i]=qi; q[j]=qj;
end; end;
if ((prf1-profit)<.000001 | iter>50) then goto f;
profit=prf1;
goto a;
f: profit=prf1;
ratio=profit;
dg=g`*b*dz;
dh=w`*dg;
u=probit(j(m,1,1)-q);
print 'Final Results';
IF MAX='RATIO' THEN GOTO R4;
IF MAX='PROFIT' THEN GOTO P4;
print iter u q dz dg dh;goto g4;
R4:print iter u q dz dg dh cost ratio;
GOTO G4;
P4:print iter u q dz dg dh cost profit;
g4:finish mod;run mod;
print 'B   - transformation matrix from y to z';
print 'ITER- number of iterations             ';
print 'U   - truncation points                ';
print 'Q   - proportions selected             ';
print 'DZ  - selection differentials          ';
print 'DG  - genetic gains                    ';
print 'DH  - total genetic gain               ';
print 'COST- cost of measurement              ';
print 'RATIO  - gain to cost ratio            ';
print 'PROFIT - profit                        ';