program contract(input,output);
{ read nhood.v, apply contraction M, write contracted.v }
const
  hdim=13;
  edim=27;
  lmax=29;
  pdeg=100;
  vdim=1001;

%include 'types.h';

var
  cov,beta2,sigma,rho,tau: scalar;
  vkappa,v,vin,vout,v0,v1,v2,v3,v4: vector;
  m: matrix;

%include 'reps.i';
%include 'scalar.i';
%include 'vector.i';
%include 'linear.i';
%include 'hiemaps.i';

begin { contract }
  sinit;
  writeln('contract: reading parameters and matrix m');
  sread('rho.s',rho);
  sread('beta2.s',beta2);
  sread('cov.s',cov);
  sread('sigma.s',sigma);
  sread('tau.s',tau);
  mread('m.m',m);
  writeln;

  writeln('contract: reading vectors');
  vread('kappa.v',rho,vkappa);
  vread('approx.v',rho,v);
  vread('difference.v',rho,v0);
  vread('nhood.v',rho,vin);
  writeln;

  writeln('contract: compute terms linear in nhood');
  dmm1(v,vin,v1);
  writeln;
  dmm2(v,vin,v2);
  writeln;
  dmm3(v,vin,v3);
  writeln;

  writeln('contract: compute term quadratic in nhood');
  nn(vdim,vin,vout);
  vmprod(m,vout,v4);
  writeln;

  writeln('contract: add all to constant term');
  vout:=v0;
  vadd(v4,vout);
  vadd(v3,vout);
  vadd(v2,vout);
  vadd(v1,vout);
  vwrite('contracted.v',vout);
  writeln;

  writeln('contract: checking image of approx+nhood');
  if vsub(vout,vin) then writeln('contract: ok')
  else writeln('contract: error - not mapped into itself');

  sdone('contract');
end { contract }.