out aax;
%off nat;
on rounded;

procedure delta(ii,jj);
begin 
if ii eq jj then return 1 else return 0; 
end;

procedure addone(x); rrr:=x+1;

% random numbers with uniform distribution generated in Mathematica 
%no=Table[Random[Real],{i,1,100 }]

noise:=
{0.7611189572266488, 0.809932949360867, 
  0.32754272053822, 0.3711809230728502, 
  0.2198155852281288, 0.2810910296348176, 
  0.5994259353572776, 0.987082750593675, 
  0.6503306993502929, 0.6922800675720912, 
  0.7549523349478609, 0.445705988068098}$

n:=2;
m:=2;
%--------------------------------------------------
A(1):=10$ A(2):=30$
AA:=TP(MAT( (A(1), A(2) )));
%--------------------------------------------------

matrix KM(n,n);
for tt:=1:n do for pp:=1:m do km(tt,pp):=k(tt,pp);

k(1,1):=0.6$ k(1,2):=0.4$
k(2,1):=0.4$ k(2,2):=0.6$

%-----------------
RR:=TP(MAT( (R(1), R(2) )));
NNO:=TP(MAT( (NO(1), NO(2) )));
NO(1):=100; NO(2):=20;

RR:= KM*AA + NNO;


%--------------- matrix formulation ---------------------------------
matrix RM(n,n);
for ii:=1:n do rm(ii,ii):=r(ii);


%---- SYMMETRY -------------------------------------
matrix ONE(n,n);
for ii:=1:n do ONE(ii,ii):=1;


matrix SS(n,n);
for tt:=1:n do for pp:=1:m do SS(tt,pp):=S(tt,pp);

S(1,2):= - RR(1,1)*(S(1,1) - 1)/RR(2,1);
S(2,1):= - RR(2,1)*(S(2,2) - 1)/RR(1,1);

S(1,1):=2; S(2,2):=2;

SS*RR - RR;

%---- SYMMETRYX -------------------------------------



FACTOR A;

SS;
SS**-1;
                             

SS*KM*AA - KM*AA;
SS*NNO - NNO;

%----------------
SSX:=SS-ONE;
SSX**-1;


SSX*KM*AA;
SSX*NNO;

KM*AA;
NNO;


YY:=TP(MAT( (Y(1), Y(2) )));

SSX*YY;

shut aax;
bye;
end;