#include <oxstd.h>
#include <pcfiml.h>

#pragma  link("maximize.oxo")
#pragma  link("database.oxo")
#pragma  link("pcfiml.oxo")
//#pragma  link("nortest.oxo")


olsnorm(my, mx, amb, amxtxi, amxtx)
{
    decl ml, md;

    decldl(amxtx[0] = mx'mx, &ml, &md);
    amb[0] = solveldl(ml, md, mx'my);
    amxtxi[0] = solveldl(ml, md, unit(columns(mx)));

//  print(amxtx[0], amxtxi[0], amb[0]);
}
olssvd(my, mx, amb, amxtxi, amxtx)
{
    decl mu, mw, mv;

    decsvd(mx, &mu, &mw, &mv);
    amxtx[0] = mv*diag(mw.^2)*mv';
    amxtxi[0] = mv*diag(mw.^-2)*mv';
    amb[0] = mv*diag(mw.^-1)*(mu'my);
//  amxtx[0] = (mv.*(mw.^2))*mv';
//  amxtxi[0] = (mv./(mw.^2))*mv';
//  amb[0] = (mv./mw)*(mu'my);
//  print(amxtx[0], amxtxi[0], amb[0]);
}

decl mX, mY, cT, cY, cX, fTrans;
generate(const cy, const ck, const ct)
{
    mX = ranu(ct,ck);
    mY = rann(ct,cy) / 10 + mX * ones(ck,1);
    cT = ct;  cY = cy;  cX = ck;
    fTrans = FALSE;
}
transpose()
{
    mX = mX';
    mY = mY';
    fTrans = TRUE;
}
olsnortime(const crep)
{
    decl mb, mxtx, mxtxi, i, time;

    time = timer();
    for (i = 0; i < crep; i++)
        olsnorm(mY, mX, &mb, &mxtxi, &mxtx);

    print("norm: k=", cX, " T=", cT, " n=", cY,
          " rep=", crep, " time lapsed: ", timespan(time), "\n");
}
olsqrtime(const crep)
{
    decl mb, mxtx, mxtxi, i, time;

    time = timer();
    if (fTrans)
        for (i = 0; i < crep; i++)
            olsr(mY, mX, &mb, &mxtxi, &mxtx);
    else
        for (i = 0; i < crep; i++)
            olsc(mY, mX, &mb, &mxtxi, &mxtx);
    print("qr:   k=", cX, " T=", cT, " n=", cY,
          " rep=", crep, " time lapsed: ", timespan(time), "\n");
}
olsinvtime(const crep)
{
    decl mb, mxtx, mxtxi, i, time;

    time = timer();
    if (fTrans)
        for (i = 0; i < crep; i++)
            ols2r(mY, mX, &mb, &mxtxi, &mxtx);
    else
        for (i = 0; i < crep; i++)
            ols2c(mY, mX, &mb, &mxtxi, &mxtx);

    print("inv:  k=", cX, " T=", cT, " n=", cY,
          " rep=", crep, " time lapsed: ", timespan(time), "\n");
}
olssvdtime(const crep)
{
    decl mb, mxtx, mxtxi, i, time;

    time = timer();
    for (i = 0; i < crep; i++)
        olssvd(mY, mX, &mb, &mxtxi, &mxtx);

    print("acc:  k=", cX, " T=", cT, " n=", cY,
          " rep=", crep, " time lapsed: ", timespan(time), "\n");
}
main()
{
	print("data in columns\n");
    generate(1,5,50);
    olsinvtime(2000);  olsqrtime(2000);  //olssvdtime(2000);
    generate(1,5,500);
    olsinvtime(200);   olsqrtime(200);   //olssvdtime(200);
    generate(1,5,5000);
    olsinvtime(20);    olsqrtime(20);    //olssvdtime(20);

	print("data in rows\n");
    generate(1,5,50);   transpose();
    olsinvtime(2000);  olsqrtime(2000);  //olssvdtime(2000);
    generate(1,5,500);  transpose();
    olsinvtime(200);   olsqrtime(200);   //olssvdtime(200);
    generate(1,5,5000); transpose();
    olsinvtime(20);    olsqrtime(20);    //olssvdtime(20);

	print("data in columns\n");
    generate(4,5,50);
    olsinvtime(2000);  olsqrtime(2000);  //olssvdtime(2000);
    generate(4,5,500);
    olsinvtime(200);   olsqrtime(200);   //olssvdtime(200);
    generate(4,5,5000);
    olsinvtime(20);    olsqrtime(20);    //olssvdtime(20);

	print("data in rows\n");
    generate(4,5,50);   transpose();
    olsinvtime(2000);  olsqrtime(2000);  //olssvdtime(2000);
    generate(4,5,500);  transpose();
    olsinvtime(200);   olsqrtime(200);   //olssvdtime(200);
    generate(4,5,5000); transpose();
    olsinvtime(20);    olsqrtime(20);    //olssvdtime(20);

	print("data in columns\n");
    generate(4,5,50);
    decl pcfiml = new PcFiml();
    pcfiml->Create(1, 1, 1, 50, 1);                    // create the database
    pcfiml->Append(mY, {"y1", "y2", "y3", "y4"}, 0);
    pcfiml->Append(mX, {"x1", "x2", "x3", "x4", "x5"}, 0);
    pcfiml->Select(Y_VAR, { "y1", 0, 0, "y2", 0, 0, "y3", 0, 0, "y4", 0, 0} );
    pcfiml->Select(X_VAR, { "x1", 0, 0, "x2", 0, 0, "x3", 0, 0, "x4", 0, 0, "x5", 0, 0} );
    pcfiml->SetSample(1, 1, 1, 50);
    pcfiml->SetPrint(FALSE);

    decl i, crep = 2000, time = timer();
    for (i = 0; i < crep; i++)
        pcfiml->Estimate();

    print("pcf:  k=", cX, " T=", cT, " n=", cY,
          " rep=", crep, " time lapsed: ", timespan(time), "\n");
}

