/*--------------------------------------------------------------------------
 * simula.ox - code for Simulation class
 *
 *       (C) Jurgen Doornik 1994-2001
 *
 * 12- 2- 01: added support for two-sided tests
 *--------------------------------------------------------------------------*/

#include <oxstd.h>
#include <oxfloat.h>
#include <simula.h>
#import <database>

Simulation::Simulation(const mT, const mxT, const cRep, const fCommon,
    const dSeed, const mPvalue, const mTrueParam)
{
    m_mT =        vecr(matrix(mT))';
    m_mxT =       mxT;
    m_cRep =      cRep;
    m_fCommon =   fCommon;
	
	if (dSeed == -1 || isarray(dSeed))
	    m_dSeed = dSeed;
	else
	{	m_dSeed = ranseed(0);
		if (isarray(m_dSeed))
	    	m_dSeed[0] = dSeed;
		else
	    	m_dSeed = dSeed;
	}

	setPvalues(mPvalue);
    m_mTrueParam = vecr(mTrueParam);				 /* force column vector */

	m_bRecursive = m_bStore = FALSE;
	m_mCoefMeanRec = m_mCoefMeanSqRec = m_mRejectRec =
		m_mCritValRec = m_mTestMomRec = 0;
	m_iPlotRep =  -1;

    m_asTest = new array[0];
    m_asCoef = new array[0];
}
Simulation::SetRecursive(const bRecursive)
{
	m_bRecursive = bRecursive;
}
Simulation::SetStore(const bStore)
{
	m_bStore = bStore;
}
Simulation::SetPlotRep(const iPlotRep)
{
	m_iPlotRep = iPlotRep;
}
Simulation::SetTestNames(const asNames)
{
    decl i;

    m_asTest = new array[sizeof(asNames)];
    for (i = 0; i < sizeof(asNames); ++i)
        m_asTest[i] = asNames[i];
}
Simulation::SetCoefNames(const asNames)
{
    decl i;

    m_asCoef = new array[sizeof(asNames)];
    for (i = 0; i < sizeof(asNames); ++i)
        m_asCoef[i] = asNames[i];
}

Simulation::Simulate()
{
    decl ct, i, t, itno, creject, isto, coefs, pvalues, tests;
    decl fcoefs, fpvalues, ftests;

    ct = columns(m_mT);
    ranseed(m_dSeed);                                       /* set the seed */
    for (i = m_cRejected = 0; i < ct; ++i)
    {
        t = m_mT[0][i];
        if (m_fCommon)                         /* if common random numbers: */
            ranseed(m_dSeed);                               /* set the seed */

        for (itno = 0; itno < m_cRep; ++itno)  /* iterate over replications */
        {
            creject = 0;
            while (!Generate(itno, t, m_mxT))  /* generate this replication */
            {   ++m_cRejected;  ++creject;
                if (creject > 10 && imod(creject, 10) == 0)
                    eprint("Warning: ", creject, " experiments rejected for replication ", itno, "\n");
            }

            if (itno == 0)        /* do initialization in first replication */
            {
				doFirstRep(&ftests, &fpvalues, &fcoefs);
            }

			isto = m_bStore ? itno : 0;
            if (fcoefs)
            {   coefs = GetCoefficients();
                m_mCoefVal[][isto] = coefs;
                m_mCoefMean += coefs;
                m_mCoefMeanSq += coefs .^ 2;
            }
            if (fpvalues)
                m_mReject += GetPvalues() .< m_mPvalueExp;
            if (ftests)
                m_mTestVal[][itno] = GetTestStatistics();

			if (m_iPlotRep > 0 && imod(itno + 1, m_iPlotRep) == 0)
				Plot(itno, t);
        }

        if (ftests)
        {   
			if (sizerc(m_mIdxTwoSided) == 0)
				 m_mCritVal	= quantiler(m_mTestVal, 1 - m_mPvalue);
			else
			{	decl q = quantiler(m_mTestVal, 1 - m_mPvalue ~ 1 - m_mPvalue / 2 ~ m_mPvalue / 2);
				m_mCritVal = q[][ : sizec(m_mPvalue) - 1];
				m_mCritVal[m_mIdxTwoSided][] = q[m_mIdxTwoSided][sizec(m_mPvalue) : ];
				m_mCritVal1 = q[][2 * sizec(m_mPvalue) :];
			}
        }

        m_mCoefMean /= m_cRep;
        m_mCoefMeanSq /= m_cRep;
        m_mReject /= m_cRep;

		if (m_bRecursive)
			stoRec(i, ftests, fpvalues, fcoefs);
		if (!m_bRecursive || i == ct - 1)
			report(t, ftests, fpvalues, fcoefs);

		if (m_iPlotRep >= 0)
			Plot(m_cRep, t);
    }
}

Simulation::setPvalues(const mPvalue)
{
	decl i;

    m_mPvalue = vecr(mPvalue)';						    /* force row vector */
	m_mPvalueExp = 0;							 /* used for dot comparison */

	m_asPvalue = m_asPvalue1 = new array[columns(mPvalue)];
    for (i = 0; i < columns(mPvalue); ++i)
	{
        m_asPvalue[i] = sprint(mPvalue[0][i] * 100, "%");
        m_asPvalue1[i] = sprint(mPvalue[0][i] * 50, "%");
	}
    m_mIdxTwoSided = m_mIdxOneSided = <>;
}

Simulation::IsTwoSided()
{
	return <>;
}
Simulation::doFirstRep(const abTests, const abPvalues, const abCoefs)
{
    decl fcoefs, fpvalues, ftests, r, c, mval;

	fcoefs = fpvalues = ftests = FALSE;
	m_mCoefMean = m_mCoefMeanSq = m_mReject = m_mCritVal = m_mCritVal1 = 0;
	m_mTestVal = m_mCoefVal = 0;

	mval = GetCoefficients();
	if (isdouble(mval)) mval = matrix(mval);
    r = rows(mval);
    c = columns(mval);
    if (r > 0 && c > 0)
    {   fcoefs = TRUE;
        if (r != rows(m_mTrueParam))	/* wrong #rows: truncate/pad with 0 */
			m_mTrueParam = shape(m_mTrueParam, r, 1);
        m_mCoefVal = zeros(r, m_bStore ? m_cRep : 1);
    }
	
	mval = GetPvalues();
	if (isdouble(mval)) mval = matrix(mval);
    r = rows(mval);
    c = columns(mval);
    if (r > 0 && c > 0)
    {   fpvalues = TRUE;
        m_mPvalueExp = ones(r, 1) ** m_mPvalue;
    }
	
	mval = GetTestStatistics();
	if (isdouble(mval)) mval = matrix(mval);
    r = rows(mval);
    c = columns(mval);
    if (r > 0 && c > 0)
    {
		ftests = TRUE;
        m_mTestVal = zeros(r, m_cRep);
		m_mIdxOneSided = range(0, r - 1);
		
		decl istwosided = IsTwoSided();
		if (sizerc(istwosided) == r)
		{
			m_mIdxTwoSided = vecindex(istwosided);
			m_mIdxOneSided = vecindex(!istwosided);
		}
    }
	abTests[0] = ftests;
	abPvalues[0] = fpvalues;
	abCoefs[0] = fcoefs;
}

Simulation::report(const t, const bTests, const bPvalues, const bCoefs)
{
	decl m;
	
    print("\nT=", t, ", M=", m_cRep, ", seed=", m_dSeed,
        (m_fCommon) ? " (common)" : "", "\n");
    if (m_cRejected)
        print(m_cRejected, " experiments rejected!\n");
    if (bTests)
    {
		m = standardize(m_mTestVal') ';
		decl n = columns(m);
		print("\nmoments of test statistics",
              "%r", m_asTest, "%c", {"mean", "std.dev", "skewness", "ex.kurtosis"},
			  	meanr(m_mTestVal) ~ sqrt(varr(m_mTestVal)) ~
				sumr(m .^ 3) / n ~ sumr(m .^ 4) / n - 3);

		if (sizerc(m_mIdxOneSided))
		{	print("\ncritical values (tail quantiles)",
    	          "%r", array(m_asTest[m_mIdxOneSided]), "%c", m_asPvalue,
				 m_mCritVal[m_mIdxOneSided][]);
		}
		if (sizerc(m_mIdxTwoSided))
		{	print("\ncritical values (two sided: left tail quantiles)",
	              "%r", array(m_asTest[m_mIdxTwoSided]), "%c", m_asPvalue1,
				  m_mCritVal1[m_mIdxTwoSided][]);
			print("\ncritical values (two sided: right tail quantiles)",
    	          "%r", array(m_asTest[m_mIdxTwoSided]), "%c", m_asPvalue1,
				  m_mCritVal[m_mIdxTwoSided][]);
		}
	}
    if (bPvalues)
        print("\nrejection frequencies",
              "%r", m_asTest ~ {"[ASE]"}, "%c", m_asPvalue, m_mReject |
			  	sqrt(m_mPvalue .* (1 - m_mPvalue) / m_cRep));
    if (bCoefs)
    {
		m = m_mCoefMean ~ sqrt(m_mCoefMeanSq - m_mCoefMean .^ 2);
        print("\nmoments of estimates",
              "%r", m_asCoef,
              "%c", { "mean", "MCSD"},
              m);
		m = m_mCoefMean - m_mTrueParam
            ~ sqrt( (m_mCoefMeanSq - m_mCoefMean .^ 2) / m_cRep)
            ~ sqrt(m_mCoefMeanSq - 2 * m_mCoefMean .* m_mTrueParam + m_mTrueParam .^ 2)
			~ m_mTrueParam;
        print("\nbiases of estimates",
              "%r", m_asCoef,
              "%c", { "mean bias", "MCSE", "RMSE", "true value"},
              m);
    }
}
Simulation::stoRec(const i, const bTests, const bPvalues, const bCoefs)
{
	if (i == 0)
	{
		decl ct = columns(m_mT);

		if (bCoefs)
		{	m_mCoefMeanRec = constant(M_NAN, ct, rows(m_mCoefMean));
			m_mCoefMeanSqRec = constant(M_NAN, ct, rows(m_mCoefMeanSq));
		}
        if (bPvalues)
			m_mRejectRec = constant(M_NAN, ct, rows(m_mReject) * columns(m_mReject));
		if (bTests)
		{	m_mCritValRec = constant(M_NAN, ct, rows(m_mCritVal));
			m_mTestMomRec = constant(M_NAN, ct, 2 * rows(m_mTestVal));
		}
	}
    if (bCoefs)
	{	m_mCoefMeanRec[i][] = m_mCoefMean';
		m_mCoefMeanSqRec[i][] = m_mCoefMeanSq';
	}
    if (bPvalues)
		m_mRejectRec[i][] = vec(m_mReject') ';
	if (bTests)
	{
		m_mCritValRec[i][] = m_mCritVal';
		m_mTestMomRec[i][] = meanr(m_mTestVal)' ~ sqrt(varr(m_mTestVal))';
	}
}

Simulation::GetCoefficients()
{
    return 0;
}
Simulation::GetPvalues()
{
    return 0;
}
Simulation::GetTestStatistics()
{
    return 0;
}
Simulation::Plot(const iRep, const iT)
{
}

Simulation::SaveIn7(const sFilename)
{
    decl fstore = FALSE, i;
    decl db = new Database();

	db.Create(1, 1, 1, m_cRep, 1);

    if (columns(m_mTestVal) == m_cRep)
    {
		db.Append(m_mTestVal', m_asTest, 0);
        fstore = TRUE;
    }
    if (columns(m_mCoefVal) == m_cRep)
    {
		db.Append(m_mCoefVal', m_asCoef, 0);
        fstore = TRUE;
    }
	if (fstore)
    	db.SaveIn7(sFilename);
    delete db;
	
	return fstore;
}
Simulation::SaveRecIn7(const sFilename)
{
    decl fstore = FALSE, ct = columns(m_mT), i, j;
    decl db = new Database();

	db.Create(1, 1, 1, ct, 1);

    db.Append(m_mT', {"T"}, 0);

    if (rows(m_mCoefMeanRec) == ct && rows(m_mCoefMeanSqRec) == ct)
    {
		for (i = 0; i < columns(m_mCoefMeanRec); ++i)
		{	
			if (!isnan(m_mTrueParam[i]))
			{	decl mcsd = sqrt(m_mCoefMeanSqRec[][i] - m_mCoefMeanRec[][i] .^ 2);
				decl bias = m_mCoefMeanRec[][i] - m_mTrueParam[i];
				decl rmse = sqrt(sqr(mcsd) + sqr(bias));
				
        		db.Append(m_mCoefMeanRec[][i], m_asCoef[i]);
        		db.Append(mcsd, "MCSD_" ~ m_asCoef[i]);
				db.Append(bias, "Bias_" ~ m_asCoef[i]);
				db.Append(rmse, "RMSE_" ~ m_asCoef[i]);
				fstore = TRUE;
			}
		}
    }
	
    if (rows(m_mCritValRec) == ct)
    {   db.Append(m_mCritValRec, m_asTest, 0);
        fstore = TRUE;
    }
    if (rows(m_mRejectRec) == ct)
    {
		decl c = sizeof(m_asTest), as = new array[c], sproc;
		
		for (j = 0; j < columns(m_mPvalue); ++j)
		{	sproc = sprint("%.2g", m_mPvalue[0][j]);
			if (sproc[0] == '0') sproc = sproc[1:];
			for (i = 0; i < c; ++i)
				as[i] = sprint(m_asTest[i], ":", sproc);
			db.Append(m_mRejectRec[][j*c : (j+1)*c-1], as, 0);
		}
        fstore = TRUE;
    }
    if (rows(m_mTestMomRec) == ct)
    {
		decl c = columns(m_mTestMomRec) / 2;
		for (i = 0; i < c; ++i)
		{	
			db.Append(m_mTestMomRec[][i], "mean_" ~ m_asTest[i]);
			db.Append(m_mTestMomRec[][c + i], "sd_" ~ m_asTest[i]);
			fstore = TRUE;
		}
	}
	if (fstore)
    	db.SaveIn7(sFilename);
    delete db;
	
	return fstore;
}
