#include "libscl.h"

using namespace std;
using namespace scl;

int main(int argc, char** argp, char** envp)
{
  INTEGER n = 50;
  INTEGER K = 3;

  realmat y(n,1);
  realmat X(n,K);
  realmat b(K,1);

  b[1] = 3.0;
  b[2] = 2.0;
  b[3] = 1.0;

  REAL sd = 0.1;

  INT_32BIT seed = 770116;

  /*
  ofstream dout("quantreg.dat");
  if (!dout) error("Error, tstqreg, cannot open quantreg.dat to write");
  */

  for (INTEGER t=1; t<=n; ++t) {
    X(t,1) = 1.0;
    X(t,2) = (t - REAL(n)/2.0)/REAL(n);
    X(t,3) = pow(X(t,2),2);
    y[t] = X(t,1)*b[1] + X(t,2)*b[2] + X(t,3)*b[3] + sd*unsk(seed);
    /*
    dout << fmt('f',20,8,y[t]) << fmt('f',5,1,X(t,1)) 
         << fmt('f',20,8,X(t,2)) << fmt('f',20,8,X(t,3)) << '\n';
    */
  } 

  /*
  dout.close();

  ifstream din("quantreg.dat");
  if (!din) error("Error, tstqreg, cannot open quantreg.dat to read");

  for (INTEGER t=1; t<=n; ++t) {
    din >> y[t] >> X(t,1) >> X(t,2) >> X(t,3);
  } 

  din.close();
  */

  cout << starbox("/true b//") 
       << b << '\n';

  cout << starbox("/quantile regression b//") << '\n';

  REAL p = 0.05;
  b = quantreg(y,X,p);
  realmat e = y - X*b;
  REAL sum = 0.0;
  for (INTEGER t=1; t<=n; ++t) if (e[t] > 0) ++sum;
  cout << "\n\t" << " p = " << p 
       << "\n\t" << " sum = " << sum 
       << '\n' << b << '\n';

  p = 0.25;
  b = quantreg(y,X,p);
  e = y - X*b;
  sum = 0.0;
  for (INTEGER t=1; t<=n; ++t) if (e[t] > 0) ++sum;
  cout << "\n\t" << " p = " << p 
       << "\n\t" << " sum = " << sum 
       << '\n' << b << '\n';

  p = 0.5;
  b = quantreg(y,X,p);
  e = y - X*b;
  sum = 0.0;
  for (INTEGER t=1; t<=n; ++t) if (e[t] > 0) ++sum;
  cout << "\n\t" << " p = " << p 
       << "\n\t" << " sum = " << sum 
       << '\n' << b << '\n';

  p = 0.75;
  b = quantreg(y,X,p);
  e = y - X*b;
  sum = 0.0;
  for (INTEGER t=1; t<=n; ++t) if (e[t] > 0) ++sum;
  cout << "\n\t" << " p = " << p 
       << "\n\t" << " sum = " << sum 
       << '\n' << b << '\n';

  p = 0.95;
  b = quantreg(y,X,p);
  e = y - X*b;
  sum = 0.0;
  for (INTEGER t=1; t<=n; ++t) if (e[t] > 0) ++sum;
  cout << "\n\t" << " p = " << p 
       << "\n\t" << " sum = " << sum 
       << '\n' << b << '\n';

  return 0;
}
