dm output 'clear';
dm log 'clear';	

/****************************************************************************/
/*	PREVOD KORELACNI MATICE NA POZITIVNE DEFINITNI							*/
/*														Ing. Zdenka Vesela	*/
/*														14.11.2006			*/
/****************************************************************************/

/*  Upraveno pro potřeby výpočtu plemenných hodnot pro přírůstky býků v odchovnách

Hana Krejčová
8. ledna 2007   */
/******************************* ZADANI MAKER *******************************/
%LET a=5;	  /*... zadani rozmeru matice (pocet korelovanych vlastnosti ...*/
%LET max=1000000;  /*... maximalni pocet iteraci						 ...*/
/************** ZADANI ADRESARE A NAZVU UKLADANE KORELACNI MATICE ***********/
filename mat_gen 'e:/hana/matice/mat_gen';

proc iml;
/*********************** ZADANI KORELACNI MATICE genetické****************************/
r={
0.1100E-02  0.1390E-03 -0.3168E-03  0.3383E-03  0.3496E-03,
  0.1390E-03  0.4587E-03  0.1078E-03 -0.5268E-04  0.2988E-03,
 -0.3168E-03  0.1078E-03  0.6784E-03  0.2050E-03  0.3439E-03,
  0.3383E-03 -0.5268E-04  0.2050E-03  0.3530E-03  0.2792E-03,
  0.3496E-03  0.2988E-03  0.3439E-03  0.2792E-03  0.5411E-03
};
/*********************** ZADANI MATICE S POCTY ZVIRAT************************/
p={
1 1 1 1 1,
1 1 1 1 1,
1 1 1 1 1,
1 1 1 1 1,
1 1 1 1 1
};
/****************************************************************************/
/*... 1. Vypocet matice vlastnich vektoru (eigenvectors) U
	a diagonalni matice vlastnich hodnot  (eigenvalues) D
	korelacni matice R										...*/
do i=1 to &max by 1;
call eigen(b,u,r);
d=diag(b);
d1=d;
if d[&a,&a]>0 then do;
	create r from r;
	append from r;		/*... Vytvoreni souboru s korelacni matici ...*/
	print r;
	print i;
	r=j(1,1,0);
	d=j(1,1,0);
	p=j(1,1,0);
	u=j(1,1,0);
	a=root(r); /*... Kontrola správnosti ...*/
end;
do j=1 to &a by 1;
	if d1[j,j]<0.0001 then do;
	d1[j,j]=0.0002;
	end;
end;
trd=trace(d);
trd1=trace(d1);
trdd1=trd/trd1;
dd=j(&a,&a,0);
do j=1 to &a by 1;
	dd[j,j]=d1[j,j]*trdd1;
end;
w=j(&a);
do j=1 to &a by 1;
	do k=1 to &a by 1;
		w[j,k]=1/p[j,k];
	end;
end;
/*... Vypocet nove korelacni matice ...*/
r1=r-(r-(u*dd*u`))#w;
r2=j(&a);
do j=1 to &a by 1;
	do k=1 to &a by 1;
		r2[j,k]=r1[j,k]/sqrt(r1[j,j]*r1[k,k]);
	end;
end;
r=r2;
end;

data mat_gen; set r; file mat_gen; 
	put col1 col2 col3 col4 col5 ;
	 /*... Polde potreby doplnit nebo umazat sloupce ...*/
run;
quit;


filename mat_tp 'e:/hana/matice/mat_tp';

proc iml;
/*********************** ZADANI KORELACNI MATICE prostřeďové****************************/
r={
0.8092E-03  0.1297E-02  0.1078E-02  0.6046E-03 -0.3391E-04,
  0.1297E-02  0.2119E-02  0.1724E-02  0.9623E-03 -0.9404E-04,
  0.1078E-02  0.1724E-02  0.1441E-02  0.8088E-03 -0.3626E-04,
  0.6046E-03  0.9623E-03  0.8088E-03  0.4597E-03 -0.1271E-04,
 -0.3391E-04 -0.9404E-04 -0.3626E-04 -0.1271E-04  0.5279E-04

};
/*********************** ZADANI MATICE S POCTY ZVIRAT************************/
p={
1 1 1 1 1,
1 1 1 1 1,
1 1 1 1 1,
1 1 1 1 1,
1 1 1 1 1
};
/****************************************************************************/
/*... 1. Vypocet matice vlastnich vektoru (eigenvectors) U
	a diagonalni matice vlastnich hodnot  (eigenvalues) D
	korelacni matice R										...*/
do i=1 to &max by 1;
call eigen(b,u,r);
d=diag(b);
d1=d;
if d[&a,&a]>0 then do;
	create r from r;
	append from r;		/*... Vytvoreni souboru s korelacni matici ...*/
	print r;
	print i;
	r=j(1,1,0);
	d=j(1,1,0);
	p=j(1,1,0);
	u=j(1,1,0);
	a=root(r); /*... Kontrola správnosti ...*/
end;
do j=1 to &a by 1;
	if d1[j,j]<0.0001 then do;
	d1[j,j]=0.0002;
	end;
end;
trd=trace(d);
trd1=trace(d1);
trdd1=trd/trd1;
dd=j(&a,&a,0);
do j=1 to &a by 1;
	dd[j,j]=d1[j,j]*trdd1;
end;
w=j(&a);
do j=1 to &a by 1;
	do k=1 to &a by 1;
		w[j,k]=1/p[j,k];
	end;
end;
/*... Vypocet nove korelacni matice ...*/
r1=r-(r-(u*dd*u`))#w;
r2=j(&a);
do j=1 to &a by 1;
	do k=1 to &a by 1;
		r2[j,k]=r1[j,k]/sqrt(r1[j,j]*r1[k,k]);
	end;
end;
r=r2;
end;

data mat_tp; set r; file mat_tp; 
	put col1 col2 col3 col4 col5 ;
	 /*... Polde potreby doplnit nebo umazat sloupce ...*/
run;
quit;