source: tol/tol_pkg/MulVarFunAprox/test/test_0001/test.tol @ 345a4d7

/* -*- mode: C++ -*- */
//////////////////////////////////////////////////////////////////////////////
// FILE   : test.tol
// PURPOSE: Package MulVarFunAprox test
//////////////////////////////////////////////////////////////////////////////

Text email_ = "vdebuen@tol-project.org"; 
Text link_ = ""; 
Text summary_ = "Test for mutivariant function aproximation";

//Write here the initial test TOL code if needed
Real numErr0 = Copy(NError);
Real numWar0 = Copy(NWarning);

//Random seed setting
Real PutRandomSeed(0);
Real rndSeed = GetRandomSeed(0);
WriteLn("Current Random Seed = "<<rndSeed);

//#Require MulVarFunAprox;
#Embed "../../../_embed_all.tol";

Real n = 2;


Set Base.VFun = {SetOfCode(
  VMatrix p00 (VMatrix x) { Rand(VRows(x),1,1,1) },
  VMatrix p10 (VMatrix x) { SubCol(x,[[1]]) },
  VMatrix p01 (VMatrix x) { SubCol(x,[[2]]) },
  VMatrix p20 (VMatrix x) { RPow(SubCol(x,[[1]]),2) },
  VMatrix p02 (VMatrix x) { RPow(SubCol(x,[[2]]),2) },
  VMatrix p11 (VMatrix x) { SubCol(x,[[1]])$*SubCol(x,[[2]]) }
)};

Set Base.VDerFst = {[[
{SetOfCode(
  VMatrix p00 (VMatrix x) { Rand(VRows(x),1,0,0) },
  VMatrix p10 (VMatrix x) { Rand(VRows(x),1,1,1) },
  VMatrix p01 (VMatrix x) { Rand(VRows(x),1,0,0) },
  VMatrix p20 (VMatrix x) { SubCol(x,[[1]])*2 },
  VMatrix p02 (VMatrix x) { Rand(VRows(x),1,0,0) },
  VMatrix p11 (VMatrix x) { SubCol(x,[[2]]) }
)};
{SetOfCode(
  VMatrix p00 (VMatrix x) { Rand(VRows(x),1,0,0) },
  VMatrix p10 (VMatrix x) { Rand(VRows(x),1,0,0) },
  VMatrix p01 (VMatrix x) { Rand(VRows(x),1,1,1) },
  VMatrix p20 (VMatrix x) { Rand(VRows(x),1,0,0) },
  VMatrix p02 (VMatrix x) { SubCol(x,[[2]])*2 },
  VMatrix p11 (VMatrix x) { SubCol(x,[[1]]) }
)}
]]};

Set Base.VDerSnd = {[[
{[[
  {SetOfCode(
    VMatrix p00 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p10 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p01 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p20 (VMatrix x) { Rand(VRows(x),1,2,2) },
    VMatrix p02 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p11 (VMatrix x) { Rand(VRows(x),1,0,0) }
  )};
  {SetOfCode(
    VMatrix p00 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p10 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p01 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p20 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p02 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p11 (VMatrix x) { Rand(VRows(x),1,1,1) }
  )}
]]},
{[[
  {SetOfCode(
    VMatrix p00 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p10 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p01 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p20 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p02 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p11 (VMatrix x) { Rand(VRows(x),1,1,1) }
  )};
  {SetOfCode(
    VMatrix p00 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p10 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p01 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p20 (VMatrix x) { Rand(VRows(x),1,0,0) },
    VMatrix p02 (VMatrix x) { Rand(VRows(x),1,2,2) },
    VMatrix p11 (VMatrix x) { Rand(VRows(x),1,0,0) }
  )}
]]} ]]};

Real K = Card(Base.VFun);

/* */

Real J = K*10000;

VMatrix X = Group("ConcatColumns",For(1,n,VMatrix(Real i)
{
  Real min = Rand(0,0);
  Real len = Rand(5,10)*(i^(1/n));
  Rand(J,1,min,min+len)
}));


VMatrix fun.eval.intra.base(VMatrix X) 
{
  Set x = For(1,n,VMatrix(Real i) { SubCol(X,[[i]]) });
  (x[1]*0+1)*Rand(-10,10)+
  x[1]*Rand(-10,10)+
  x[2]*Rand(-10,10)+
  x[1]^2*Rand(-10,10)+
  x[2]^2*Rand(-10,10)+
  x[1]$*x[2]*Rand(-10,10)
};

VMatrix fun.eval.extra.base(VMatrix X) 
{
  Set x = For(1,n,VMatrix(Real i) { SubCol(X,[[i]]) });
  Sin(x[1])+Cos(x[1])
};

VMatrix Fx = fun.eval.extra.base(X);

Real lengthRatio = 1*(J/(K*n))^(1/n);

Set X.col.info = For(1,n,Set(Real i)
{[[
  VMatrix col = SubCol(X,[[i]]);
  Real min = VMatMin(col);
  Real max = VMatMax(col);
  Real length = max-min;
  Real grid.size = 1*length/lengthRatio
]]});

Real grid.size = SetAvr(Traspose(Extract(X.col.info,5))[1]);
Real radius = grid.size / Sqrt(2);

Set range = For(1,n,Set(Real i)
{
  Set info = X.col.info[i];
  Real points = 1+Ceil(info::length/grid.size);
  Real exceed = (points-1)*grid.size-info::length;
  Matrix rnk = Tra(DifEq(1/(1-B),Rand(points,1,1,1))-Rand(points,1,1,1));   
  Set MatSet((rnk*grid.size)+info::min-exceed/2)[1]
});

VMatrix Y = Mat2VMat(SetMat(Group("CartProd",range)));
Real H = VRows(Y);
Real N = H*K;

VMatrix Y.min = Group("ConcatColumns", For(1,n,VMatrix(Real i)
{
  Real aux = range[i][1];
  Rand(H,1,aux,aux)
}));

VMatrix Y.max = Group("ConcatColumns", For(1,n,VMatrix(Real i)
{
  Real aux = range[i][Card(range[i])];
  Rand(H,1,aux,aux)
}));


VMatrix Y.index = Mat2VMat(SetMat(Group("CartProd",For(1,n,Set(Real i)
{
  Range(1,Card(range[i]),1)
}))));

VMatrix show.X.Y = 
  ((X*0) <<  Y   ) |
  ( X    << (Y*0)); 

VMatrix closest.Y.index(VMatrix A)
{
  Real m = VRows(A);
  VMatrix cumulated = Rand(m,1,1,1);
  VMatrix A.position = Rand(m,1,1,1)+SetSum(For(1,n,VMatrix(Real i_)
  {
    Real i = n-i_+1;
    Real cRi = Card(range[i]);
    Set info = X.col.info[i];
    VMatrix Ai = SubCol(A,[[i]]);
    VMatrix Yi.min = Rand(m,1,range[i][1],range[i][1]);
    VMatrix Ri = Round((Ai-Yi.min)/grid.size);
    VMatrix minRi = Rand(m,1,0,0);
    VMatrix maxRi = Rand(m,1,cRi-1,cRi-1);
    VMatrix aux = Min(maxRi,Max(minRi,Ri));
    VMatrix pos = aux $* cumulated;
    VMatrix cumulated := cumulated*cRi;
    pos
  }))
};

VMatrix Y.closest.Y.index = closest.Y.index(Y);
VMatrix Y.closest.Y = SubRow(Y,VCol2Set(Y.closest.Y.index));
Real ok.1 = VMatMax(Abs(Y-Y.closest.Y))==0;

VMatrix Y.neighbour.pos = Min(Y.max, Y + Rand(H,n,grid.size,grid.size));
VMatrix Y.neighbour.neg = Max(Y.min, Y - Rand(H,n,grid.size,grid.size));

VMatrix Y.neighbourhood = DifEq(1/(1-B),Rand(H,1,1,1)) |
{
  Set aux.1 = BinGroup("<<", For(1,n,Set(Real i)
  {
    VMatrix pre  = SubCol(Y,Range(1,  i-1,1));
    VMatrix post = SubCol(Y,Range(i+1,n,  1));
    SetOfVMatrix(pre | SubCol(Y.neighbour.neg,[[i]]) | post,
                 pre | SubCol(Y.neighbour.pos,[[i]]) | post)
  }));
  Group("ConcatColumns",EvalSet(aux.1,closest.Y.index))
};

Set Y.neighbourhood.extrictly = For(1,H,Set(Real h) 
{ 
  VCol2Set(SubRow(Y.neighbourhood,[[h]]))  - [[h]]
});

Real X.maxNeighbour = 1+2*n;

VMatrix X.closest.Y.index = closest.Y.index(X);
VMatrix X.closest.Y = SubRow(Y,VCol2Set(X.closest.Y.index));

VMatrix X.j=DifEq(1/(1-B),Rand(J,1,1,1));

VMatrix X.h=SubRow(Y.neighbourhood, VCol2Set(X.closest.Y.index));

Set X.dif = For(1,X.maxNeighbour,VMatrix(Real i)
{
  X - SubRow(Y,VCol2Set(SubCol(X.h,[[i]])))
});

Set X.tao = For(1,X.maxNeighbour,VMatrix(Real i)
{
  VMatrix X.norm2 = RPow(X.dif[i],2)*Rand(n,1,1,1);
  RPow(X.norm2,0.5)/radius
});

Set X.tao2 = For(1,X.maxNeighbour,VMatrix(Real i)
{
  RPow(X.tao[i],2)
});

Set X.match = For(1,X.maxNeighbour,VMatrix(Real i)
{
  LE(X.tao[i],Rand(J,1,1,1))
});

Set X.lambda = For(1,X.maxNeighbour,VMatrix(Real i)
{
  RPow(X.match[i] $* (Rand(J,1,1,1)-X.tao2[i]),4)
});

VMatrix X.lambda.sum = SetSum(X.lambda);

Set X.weight = For(1,X.maxNeighbour,VMatrix(Real i)
{
  X.lambda[i] $/ X.lambda.sum
});

Matrix QX.triplet = Group("ConcatRows", For(1,X.maxNeighbour,Matrix(Real i)
{
  VMatrix X.hiK = (SubCol(X.h,[[i]])-1)*K;
  Group("ConcatRows", For(1,K,Matrix(Real k)
  {
    Code bvf = Base.VFun[k];
    VMatrix pk = bvf(X.dif[i]);
    VMatrix qk = X.weight[i]$*pk; 
    Set sel.rows = MatQuery::SelectMatch(NE(Drop(Abs(qk),1.E-15),qk*0));
    VMat2Mat(
      SubRow(X.j,sel.rows) |
     (SubRow(X.hiK,sel.rows)+k) |
      SubRow(qk,sel.rows))
  }))
}));

VMatrix QX = Pack(Convert(Triplet(QX.triplet, J, N),"Cholmod.R.Sparse"));

VMatrix displ.Z.1 = 
  Rand(H,1,-0.1,0.1)*grid.size/2 |
  Rand(H,1, 0.9,1.1)*grid.size/2;
VMatrix displ.Z.2 = 
  Rand(H,1, 0.9,1.1)*grid.size/2 |
  Rand(H,1,-0.1,0.1)*grid.size/2;
VMatrix displ.Z.3 = 
  Rand(H,1,-0.9,-1.1)*grid.size/2 |
  Rand(H,1,-0.9,-1.1)*grid.size/2;

VMatrix Z.all = (Y+displ.Z.1) << (Y+displ.Z.2) << (Y+displ.Z.3);
VMatrix Z.Y.1.all = Y << Y << Y;
VMatrix Z.Y.2.all = (Y+displ.Z.1*2) << (Y+displ.Z.2*2) << (Y+displ.Z.3*2);

VMatrix Z = Z.all;
VMatrix Z.Y.1 = Z.Y.1.all;
VMatrix Z.Y.2 = Z.Y.2.all;
VMatrix Z.h.1 = closest.Y.index(Z.Y.1);
VMatrix Z.h.2 = closest.Y.index(Z.Y.2);


VMatrix show.X.Y.Z = 
  ( X    << (Y*0) << (Z*0)) |
  ((X*0) <<  Y    << (Z*0)) | 
  ((X*0) << (Y*0) <<  Z   ); 

VMatrix show.Y.Z = 
  ( Y    << (Z*0)) | 
  ((Y*0) <<  Z   ); 


Real T = VRows(Z);
Real prior.Z.order = 2;
Real Mz = Case(
  prior.Z.order< 0, 0,
  prior.Z.order==0, T,
  prior.Z.order==1, T*(1+n*K),
  prior.Z.order==2, T*(1+(n^2)*K));

Real Z.maxNeighbour = If(Mz,2,0);
VMatrix Z.j=DifEq(1/(1-B),Rand(T,1,1,1));

VMatrix Z.h=Z.h.1 | Z.h.2;

Set Z.dif = For(1,Z.maxNeighbour,VMatrix(Real i)
{
  Z - SubRow(Y,VCol2Set(SubCol(Z.h,[[i]])))
});

Set Z.tao = For(1,Z.maxNeighbour,VMatrix(Real i)
{
  VMatrix Z.norm2 = RPow(Z.dif[i],2)*Rand(n,1,1,1);
  RPow(Z.norm2,0.5)/radius
});

Set Z.match = For(1,Z.maxNeighbour,VMatrix(Real i)
{
  LE(Z.tao[i],Rand(T,1,1,1))
});

Matrix PZ.triplet = Group("ConcatRows", For(1,Z.maxNeighbour,Matrix(Real i)
{
  VMatrix Z.hiK = (SubCol(Z.h,[[i]])-1)*K;
  Group("ConcatRows", For(1,K,Matrix(Real k)
  {
    Real cumRows = 0;
    Code bvf = Base.VFun[k];
    VMatrix pk = bvf(Z.dif[i])*If(i==1,+1,-1);
    Set sel.rows = MatQuery::SelectMatch(NE(Drop(Abs(pk),1.E-15),pk*0));
    VMatrix col.idx = (SubRow(Z.hiK,sel.rows)+k);
    {
      VMatrix row.idx = SubRow(Z.j,sel.rows)+cumRows;
      Real cumRows := cumRows + T;
      VMat2Mat(row.idx | col.idx | SubRow(pk,sel.rows))
    }  << If(prior.Z.order<1,Rand(0,3,0,0),
    Group("ConcatRows",For(1,n,Matrix(Real d)
    {
      Code bvfd = Base.VDerFst[d][k];
      VMatrix pkd = bvfd(Z.dif[i])*If(i==1,+1/2,-1/2);
      VMatrix row.idx = SubRow(Z.j,sel.rows)+cumRows;
      Real cumRows := cumRows + T;
      VMat2Mat(row.idx | col.idx | SubRow(pkd,sel.rows))
    })))<< If(prior.Z.order<2,Rand(0,3,0,0),
    Group("ConcatRows",For(1,n,Matrix(Real d1)
    {
      Group("ConcatRows",For(d1,n,Matrix(Real d2)
      {
        Code bvfd2 = Base.VDerSnd[d1][d2][k];
        VMatrix pkd2 = bvfd2(Z.dif[i])*If(i==1,+1/4,-1/4);
        VMatrix row.idx = SubRow(Z.j,sel.rows)+cumRows;
        Real cumRows := cumRows + T;
        VMat2Mat(row.idx | col.idx | SubRow(pkd2,sel.rows))
      }))
    })))
  }))
}));


Real M = J+Mz;

VMatrix PZ = Pack(Convert(Triplet(PZ.triplet, Mz, N),"Cholmod.R.Sparse"));

VMatrix Fz = Rand(Mz,1,0,0);
VMatrix F = Fx<<Fz;
VMatrix Q.all = QX << PZ;
Set Q.active = MatQuery::SelectMatch(
  NE(Tra(Rand(1,M,1,1)*Abs(Q.all)),Rand(N,1,0,0)));

VMatrix Q = SubCol(Q.all,Q.active);

VMatrix LQ = CholeskiFactor(Q,"XtX");

VMatrix alpha = CholeskiSolve(LQ,Tra(Tra(F)*Q),"PtLLtP");
VMatrix F.frc = SubRow(Q*alpha,Range(1,J,1));
VMatrix F.cmp = Fx | F.frc;
VMatrix E = Fx - F.frc;
Real E.avr = VMatAvr(E);
Real E.ste = VMatStDs(E);

Set View([[
Real n, K, J, H, 
X.maxNeighbour, Z.maxNeighbour,
T, N, Mz,M, E.avr, E.ste, "\n"]],"");

/*
Matrix Q_ = VMat2Mat(Q);

Set svd = SVD(Q_);
Matrix svd.D = Tra(SubDiag(svd[2],0));

/*
Real numErr1 = Copy(NError);
Real numWar1 = Copy(NWarning);

Set partialResults_ = [[numErr0, numErr1, 
                        ok.1, ok.2]];

//This is a messure of the success of the test 
Real quality_ = And(numErr1 == numErr0, 
                    numWar1 == numWar0, 
                    ok.1, ok.2);

//Return the results 
Set resultStr_ = @strTestStatus(summary_, link_, quality_,
                  "Partial results = "<<partialResults_,
                  "NO DBServerType", "NO DBServerHost", "NO DBAlias",
                  email_);
WriteLn(""<<resultStr_);
resultStr_;

/* */

/*
VMatrix Z.min = Group("ConcatColumns", For(1,n,VMatrix(Real i)
{
  Real aux = range[i][1];
  Rand(2*H,1,aux,aux)
}));

VMatrix Z.max = Group("ConcatColumns", For(1,n,VMatrix(Real i)
{
  Real aux = range[i][Card(range[i])];
  Rand(2*H,1,aux,aux)
}));

VMatrix Z.interior = EQ(
  And(GE(Z.all,Z.min), LE(Z.all,Z.max))*Rand(n,1,1,1),
  Rand(2*H,1,n,n));

VMatrix Z = SubRow(Z.all,MatQuery::SelectMatch(Z.interior));
VMatrix Z.Y.1 = SubRow(Z.Y.1.all,MatQuery::SelectMatch(Z.interior));
VMatrix Z.Y.2 = SubRow(Z.Y.2.all,MatQuery::SelectMatch(Z.interior));
VMatrix Z.h.1 = closest.Y.index(Z.Y.1);
VMatrix Z.h.2 = closest.Y.index(Z.Y.2);
*/

/*
Real deg = 3;

Real dim = 2;
Real maxDeg = 3;
VMatrix candidates = Mat2VMat(SetMat(Range(0,maxDeg,1)^dim));
VMatrix match = LE(candidates*Rand(dim,1,1,1),Rand(VRows(candidates),1,maxDeg,maxDeg));
VMatrix degrees = SubRow(candidates, MatQuery::SelectMatch(match));
Real K = VRows(degrees);
Class @PolynBasic
{
  Real k;
  VMatrix fun (VMatrix x) 
  {
    SetProd(For(1,dim,VMatrix(Real i)
    {
      
    })
  };
  
};
Set base = For(1,K,NameBlock(Real k)
{
});
/* */