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

#include <modelbase.h>


/*-------------------------- Modelbase : Database -----------------------------*/
Modelbase::Modelbase()
{
    if (oxversion() < 300)
        println("Requires Ox 3.00 or newer!");

    this->Database();
    ClearModel();
    SetMethod(0);
    SetPrint(TRUE);
    m_iResult = -1;

    m_cY = m_cX = m_cYlag = m_cPar = m_cTinit = m_cTforc = 0;
	m_sDbName = "";
	m_dForcCrit = 2.0;

    SetParCount(0);

    println(GetPackageName(), " package version ", GetPackageVersion(),
        ", object created on ", date());
}
Modelbase::GetDbName()
{
	return m_sDbName;
}
Modelbase::SetDbName(const sName)
{
	m_sDbName = sName;
}
///////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////
// database overrides to set model status
Modelbase::Select(const iGroup, const aSel)
{
    ClearModel();
    Database::Select(iGroup, aSel);
}
Modelbase::DeSelect()
{
    ClearModel();
    Database::DeSelect();
}
Modelbase::SetSample(const iYear1, const iPeriod1, const iYear2, const iPeriod2)
{
    SetSelSample(iYear1, iPeriod1, iYear2, iPeriod2);
}
Modelbase::SetSelSample(const iYear1, const iPeriod1, const iYear2, const iPeriod2)
{
    ClearModel();
    Database::SetSample(iYear1, iPeriod1, iYear2, iPeriod2);
}
Modelbase::Grow(const cTadd)
{
    decl i, j, ct;

	if (cTadd <= 0)
		return;

	Database::Grow(cTadd);
	
    ct = rows(m_mData);
	j = strfind(m_asNames, "Constant");
	if (j >= 0)
		m_mData[][j] = 1;

	j = strfind(m_asNames, "Trend");
	if (j >= 0)
	{
		decl t0 = m_mData[0][j];
		if (isnan(m_mData[0][j])) t0 = 1;
	
	    m_mData[][j] = range(t0, t0 + ct)';
	}
	
    if (m_iFreq == 1)
		return;

    decl season = floor(
        (fmod(range(1,ct)'+GetPeriod1()+1+m_iFreq, m_iFreq) + 1)/m_iFreq);

	j = strfind(m_asNames, "Seasonal");
	if (j >= 0)
        m_mData[][j] = season;
	j = strfind(m_asNames, "CSeasonal");
	if (j >= 0)
        m_mData[][j] = season - 1 / m_iFreq;
}
///////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////
// OxPack functions
Modelbase::SendVarStatus()
{
    return
        {{ "&Y variable", 'Y', STATUS_GROUP + STATUS_ENDOGENOUS, Y_VAR},
         { "&X variable", 'X', STATUS_GROUP, X_VAR}};
}
Modelbase::SendSpecials()
{
    return {"Constant", "Trend", "Seasonal", "CSeasonal" };
}   
Modelbase::SendFunctions()
{
    return 0;
}   
Modelbase::SendMethods()
{
    return
        {{ "Maximum Likelihood", 0, FALSE, 0}};
}
Modelbase::SendMenu(const sMenu)
{
    return 0;
}
Modelbase::SendResults(const sType)
{
    return 0;
}
Modelbase::ReceiveData()
{
    // get database sample info
    decl freq, year1, period1, year2, period2, idet, it1, cxt;
    [freq, year1, period1, year2, period2] = "OxPackGetData"("Sample");

    // create the database, allowing for lagged seasonals
    Database();
	idet = "OxPackGetData"("Deterministic");
    Create(freq, year1 - (idet > 1 ? 1 : 0), period1, year2, period2);
	it1 = idet > 1 ? freq : 0;
    Append("OxPackGetData"("Matrix"), "OxPackGetData"("Names"), it1);
	SetDbName("OxPackGetData"("DbName"));

	cxt = columns(m_mData) + 1;	 // index of trend created by Deterministic
    Deterministic(idet);
	if (it1 && columns(m_mData) > cxt && m_asNames[cxt] == "Trend")
	{
		m_mData[][cxt] -= it1;	 // ensure trend is always 1 for first obs
	}
	m_iDeterministic = idet;
	
    DeSelect();
}
Modelbase::ForceYlag(const iYgroup)
{
	decl yset = selectifc(m_mVarsel, m_mSelgroup .== iYgroup .&& m_mLagsel .== 0);
	m_mSelgroup[vecindex(sumc(m_mVarsel .== yset'))] = iYgroup;
}
Modelbase::ReceiveModel()
{
    // get selection of database variables
    Select(Y_VAR, "OxPackGetData"("SelGroup", Y_VAR));
    Select(X_VAR, "OxPackGetData"("SelGroup", X_VAR));
	ForceYlag(Y_VAR);

	// get selected sample
    decl freq, year1, period1, year2, period2;
    [year1, period1, year2, period2] = "OxPackGetData"("SelSample");
    ForceSelSample(year1, period1, year2, period2);

    // get method
    decl imethod;
    [imethod, m_cTforc, m_bRecursive] = "OxPackGetData"("Method");
    SetMethod(imethod);
    m_cTinit = "OxPackGetData"("Init");
}
Modelbase::SendDialog(const sDialog)
{
    if (sDialog == "OP_OPTIONS")
    {
        decl deps1, deps2, iprint, iitmax, bcompact;
        [iitmax, iprint, bcompact] = GetMaxControl();       
        [deps1, deps2] = GetMaxControlEps();        
        return
        {
            { "Maximization Settings" },
            { "Maximum number of iterations", CTL_INT,    iitmax,  "itmax" },
            { "Write results for every",      CTL_INT,    iprint,  "iprint" },
            { "Write in compact form",        CTL_CHECK,  bcompact,"compact"},
            { "Convergence" },
            { "Strong convergence tolerance", CTL_DOUBLE, deps1,   "eps1"},
            { "Weak convergence tolerance",   CTL_DOUBLE, deps2,   "eps2"},
            { "Default" },
            { "Reset default",                CTL_CHECK,  0,       "default"}
        };
    }
    return 0;
}   
Modelbase::ReceiveDialog(const sDialog, const asOptions, const aValues)
{
    if (sDialog == "OP_OPTIONS")
    {
        if (aValues[5] == 0)
        {   MaxControl(aValues[0], aValues[1], aValues[2]);     
            MaxControlEps(aValues[3], aValues[4]);
        }
        else
        {   MaxControl(1000, 0, 0);     
            MaxControlEps(1e-4, 5e-3);
        }
    }
	return 1;
}   
///////////////////////////////////////////////////////////////////////
// general functions
Modelbase::GetPackageName()
{
    return "Modelbase";
}
Modelbase::GetPackageVersion()
{
    return "1.0";
}
Modelbase::IsUnivariate()
{
    return TRUE;
}
Modelbase::GetModelStatus()
{
    return m_iModelStatus;
}
Modelbase::SetModelStatus(const iModelStatus)
{
    m_iModelStatus = iModelStatus;
}
Modelbase::ClearModel()
{
    m_iModelStatus = MS_NONE;
    ClearEstimation();
}
Modelbase::ClearEstimation()
{
    m_mCovar = m_mCovarRobust = -1;           /* indicates not yet computed */
    m_dLogLik = M_NAN;
    m_iResult = MAX_NOCONV;
}
Modelbase::SetPrint(fPrint)
{
    m_fPrint = fPrint;
}
Modelbase::SetResult(const iResult)
{
    m_iResult = iResult;
}
Modelbase::GetResult()
{
    return m_iResult;
}
Modelbase::SetMethod(const iMethod)
{
    m_iMethod = iMethod;
}
Modelbase::GetMethod()
{
    return m_iMethod;
}
Modelbase::GetLogLik()
{
    return m_dLogLik;
}
Modelbase::GetResiduals()
{
    return <>;
}
Modelbase::GetResVar()
{
    return <>;
}
Modelbase::GetcT()
{
    return m_cT;
}
Modelbase::GetcDfLoss()
{
    return sizerc(GetFreePar());
}
///////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////
// parameter related
Modelbase::SetParCount(const cPar, ...)
{
    decl args = va_arglist();
    
    if (m_cPar > 0 && sizeof(args) > 0 && args[0])  // add
    {
        if (cPar > 0)                               // adding cPar
        {   m_vPar = m_vPar[ : m_cPar - 1] | zeros(cPar, 1);
            m_vFixedPar = m_vFixedPar[ : m_cPar - 1] | zeros(cPar, 1);
            m_vIsFreePar = m_vIsFreePar[ : m_cPar - 1] | ones(cPar, 1);
            m_cPar += cPar;
        }
        else                                        // removing cPar
        {   m_cPar += cPar;
            m_vPar = m_vPar[ : m_cPar - 1];
            m_vFixedPar = m_vFixedPar[ : m_cPar - 1];
            m_vIsFreePar = m_vIsFreePar[ : m_cPar - 1];
        }
    }
    else                                            // set
    {   m_cPar = cPar;
        m_vPar = m_vFixedPar = zeros(cPar, 1);
        m_vIsFreePar = ones(cPar, 1);
    }
    m_vIdxFreePar = vecindex(m_vIsFreePar);
}
Modelbase::GetParCount()
{
    return m_cPar;
}
Modelbase::GetFreeParCount()
{
	return sizerc(m_vIdxFreePar);
}
Modelbase::FixPar(const iP, const dFix)
{
    m_vIsFreePar[iP] = 0;
    m_vPar[iP] = m_vFixedPar[iP] = dFix;
    m_vIdxFreePar = vecindex(m_vIsFreePar);
}
Modelbase::FreePar(const iP)
{
    if (iP < 0)
	{
		if (m_cPar)
			m_vIsFreePar[] = 1;
	}
    else
        m_vIsFreePar[iP] = 1;
    m_vIdxFreePar = vecindex(m_vIsFreePar);
}
Modelbase::GetParStatus()
{
    return { m_cPar, m_vIsFreePar, m_vFixedPar };
}
Modelbase::GetFreePar()
{
    return matrix(m_vPar[m_vIdxFreePar]);
}
Modelbase::ResetFixedPar()
{
    m_vPar = m_vIsFreePar .? m_vPar .: m_vFixedPar;
}
Modelbase::MapParToFree(const vPar)
{
    return matrix( vec(vPar)[m_vIdxFreePar] );
}
Modelbase::SetFreePar(const vParFree)
{
    m_vPar[m_vIdxFreePar] = vec(vParFree);
}
Modelbase::GetPar()
{
    return m_vPar;
}
Modelbase::SetPar(const vPar)
{
    if (m_cPar == 0)
        SetParCount(sizerc(vPar));
    if (m_cPar)
	    m_vPar[] = vec(vPar);       // total size must match
}
Modelbase::GetParNames()
{
    decl i, asy = {}, asx = {};
    GetGroupLagNames(Y_VAR, 1, 100, &asy);
    if (m_cX)
        GetGroupNames(X_VAR, &asx);
    for (i = sizeof(asx); i < m_cPar; ++i)
        asx ~= sprint("Par ", "%-2d", i + 1);
    return asy ~ asx;
}
Modelbase::GetParTypes()
{
    return 0;
}
Modelbase::GetFreeParNames()
{
	decl as = GetParNames();
	if (isarray(as))
		return array(as[m_vIdxFreePar]);
	else
		return 0;
}
Modelbase::GetLongRunNames()
{
    return 0;
}
Modelbase::GetLongRunInfo()
{
    return 0;
}
///////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////
// covariance/standard errors
Modelbase::GetStdErr()
{
    if (m_mCovar == -1)
        Covar();
    decl mstderr = zeros(1, m_cPar);
    mstderr[][m_vIdxFreePar] = sqrt(diagonal(m_mCovar));

    return mstderr';
}
Modelbase::GetStdErrRobust()
{
    if (m_mCovar == -1)
        Covar();
    if (m_mCovarRobust == -1 || sizerc(m_mCovarRobust) == 0)
        return <>;
    decl mstderr = zeros(1, m_cPar);
    mstderr[][m_vIdxFreePar] = sqrt(diagonal(m_mCovarRobust));

    return mstderr';
}
Modelbase::GetCovar()
{
    if (m_mCovar == -1)
        Covar();
    return m_mCovar;
}
Modelbase::GetCovarRobust()
{
    if (m_mCovar == -1)
        Covar();
    if (m_mCovarRobust == -1)
        return <>;
    return m_mCovarRobust;
}
///////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////
// functions to move up the model status chain
Modelbase::InitData()
{
    ClearModel();

    if (m_iT1sel < 0)
        SetSample(-1, 1, -1, 1);

    m_iT1est = m_iT1sel;  m_iT2est = m_iT2sel;

    m_mY = GetGroupLag(Y_VAR, 0, 0);
    m_cY = columns(m_mY);

    if (m_cY > 1 && IsUnivariate())
    {   println("Only univariate models allowed");
        return FALSE;
    }
    else if (m_cY == 0)
    {   println("No Y variable selected");
        return FALSE;
    }

    m_mX = GetGroupLag(Y_VAR, 1, 1000);
	m_cYlag = columns(m_mX);
    m_mX ~= GetGroup(X_VAR);
    m_cX = columns(m_mX);
        
    m_cT = m_iT2est - m_iT1est + 1;
    if (m_cT <= 2)
    {   println("Only ", m_cT, " observations, this is not enough.");
        return FALSE;
    }
    m_iModelStatus = MS_DATA;

    return TRUE;
}
Modelbase::InitPar()
{
    if (m_iModelStatus < MS_DATA)
        InitData();
    if (m_iModelStatus < MS_DATA)                              /* need data */
        return FALSE;

    ClearEstimation();
    ResetFixedPar();
    m_iModelStatus = MS_PARAMS;

    return TRUE;
}

Modelbase::SetStartPar(const vParFree)
{
    if (!Modelbase::InitPar())
        return FALSE;
    SetFreePar(vParFree);
    
    return TRUE;
}

Modelbase::DoEstimation(vPar)
{
    return vPar;
}

Modelbase::Estimate()
{
    decl  vpstart, vpfree, estout;

    ClearEstimation();
    
    if (m_iModelStatus < MS_DATA && !InitData())
    {   eprint("Failed to get data.\n");
        return FALSE;
    }
    if (m_iModelStatus < MS_PARAMS && !InitPar())
    {   eprint("Failed to get starting values.\n");
        return FALSE;
    }   
    m_iModelStatus = MS_EST_FAILED;
    vpstart = GetFreePar();         // map parameters to format used in iterative procedure
    estout = DoEstimation(vpstart); // do the estimation
    vpfree = isarray(estout) ? estout[0] : estout;
    
    SetFreePar(vpfree);             // map estimated parameters to normal format

    if (m_iResult >= MAX_CONV && m_iResult < MAX_MAXIT)
        m_iModelStatus = MS_ESTIMATED;
    else
        m_iModelStatus = MS_EST_FAILED;

    if (m_fPrint)
    {   Output();
        if (isarray(estout))
            OutputMax(estout[1], m_iResult, vpstart, estout[2]);
    }
    return m_iModelStatus == MS_ESTIMATED;
}

Modelbase::Covar()
{
}

Modelbase::Output()
{
    if (!OutputHeader(GetPackageName()))     // returns FALSE if no estimation
        return;

    OutputPar();
    OutputLogLik();
}
Modelbase::OutputHeader(const sTitle)
{
    println("\n---- ", sTitle, " ----");

    if (m_iT1sel >= 0)
	{
		if (m_iFreq > 1)
			println("The estimation sample is:  ",
    	        ObsYear(m_iT1est), " (", ObsPeriod(m_iT1est), ") - ",
        	    ObsYear(m_iT2est), " (", ObsPeriod(m_iT2est), ")");
		else
			println("The estimation sample is:  ",
    	        ObsYear(m_iT1est), " - ", ObsYear(m_iT2est));
	}
    else
        println("No sample selected yet.");

    if (m_iModelStatus != MS_ESTIMATED)
    {
        if (m_iModelStatus == MS_EST_FAILED)
        {   println("Model estimation failed: ", MaxConvergenceMsg(m_iResult));
            print("Free coefficients:", GetFreePar()');
        }
        else
            println("Model not yet estimated.");
        return FALSE;       
    }
    decl asy;
    GetGroupLagNames(Y_VAR, 0, 0, &asy);
	if (sizeof(asy) == 1)
	{   print("The dependent variable is: ", asy[0]);
		if (m_sDbName != "")
			println("  (", m_sDbName, ")");
		else
			print("\n");
	
	}
    return TRUE;
}
Modelbase::OutputPar()
{
    decl i, tval, abstval, tprob;
    decl aspar = GetParNames();
    decl aspartypes = GetParTypes();
	decl cpartypes = isarray(aspartypes) ? sizeof(aspartypes) : 0;
    decl vstderr = GetStdErr(), vrobstderr = GetStdErrRobust(), bcovarrobust = FALSE;
    decl ct = GetcT();
    decl mpar, vp = GetPar();
    decl cp = rows(vp);
    decl cdfloss = GetcDfLoss();

    mpar = vp ~ vstderr ~ vp ./ vstderr;
    if (rows(vrobstderr) > 1)
    {   mpar ~= vrobstderr ~ vp ./ vrobstderr;
        bcovarrobust = TRUE;
    }
    print("%29s", "Coefficient", "%11s", "Std.Error");
    if (bcovarrobust)
        print("%11s", "robust-SE");
    println("%9s", "t-value", "%8s", "t-prob");
    
    for (i = 0; i < cp; ++i)
    {
        if (!m_vIsFreePar[i] && mpar[i][0] == 0)
            continue;
        tval = mpar[i][2];
		if (i < cpartypes)
	        print("%-13s", aspar[i], " ", "%-2s", aspartypes[i], "%#13.6g", mpar[i][0]);
		else
	        print("%-16s", aspar[i], "%#13.6g", mpar[i][0]);
        if (!m_vIsFreePar[i])
            println("    (fixed)");
        else if (mpar[i][1] > 0)
        {
            print("%#11.4g", mpar[i][1]);
            if (bcovarrobust)
            {
                tval = mpar[i][4];
                print("%#11.4g", mpar[i][3]);
            }
			abstval = fabs(tval);
            tprob = 2 * tailt(abstval, ct - cdfloss);
            println(abstval <= 1e-4 ? "%9.2f" : abstval >= 1e3 ? "%#9.4g"
				: "%#9.3g", tval, "%8.3f", tprob);
        }
        else
            print("\n");
    }
}
Modelbase::OutputLogLik()
{
    decl ct = GetcT();
    decl cdfloss = GetcDfLoss();

    println("\nlog-likelihood", "%15.9g", m_dLogLik);
    println("no. of observations", "%10d", ct,
          "  no. of parameters  ", "%10d", cdfloss);
    println("AIC.T         ", "%15.9g", -2 * m_dLogLik + 2 * cdfloss,
          "  AIC           ", "%15.9g", (-2 * m_dLogLik + 2 * cdfloss) / ct);
    if (m_cY == 1)
    {
        decl asy;
        GetGroupLagNames(Y_VAR, 0, 0, &asy);
        println("%-16s", sprint("mean(", asy[0],")"), "%13.6g", double(meanc(m_mY)),
                "%-18s", sprint("  var(", asy[0],")"),  "%13.6g", double(varc(m_mY)));
    }
}
Modelbase::OutputMax(const sMethod, const iResult, const vPstart, const bNumerical)
{
    if (iResult >= 0)
    {
        decl deps1, deps2;
        [deps1, deps2] = GetMaxControlEps();
        println(sMethod, " using ", bNumerical ? "numerical" : "analytical",
            " derivatives (eps1=", deps1, "; eps2=", deps2, "):\n",
            MaxConvergenceMsg(iResult));
    }
    if (sizerc(vPstart))
        print("Used starting values:", vec(vPstart)', "\n");
}

///////////////////////////////////////////////////////////////////////
// Model Evaluation

Modelbase::GetForecastData(const iGroup, const mnLag, const mxLag,
	const cTforc, ...)
{
	decl mx, t1x, t1, t2, cx;
	decl t2selsave = m_iT2sel;					   /* save selection sample */
	decl itfrom = m_iT2est + 1;

	decl va = va_arglist();
	if (sizeof(va) > 0)
		itfrom = int(va[0]);
	
    m_iT2sel = rows(m_mData) - 1;		 /* full sample to get whole series */

	if (mnLag <= mxLag)
    	mx = GetGroupLag(iGroup, mnLag, mxLag);
	else
    	mx = GetGroup(iGroup);		

	m_iT2sel = t2selsave; 						  /* reset selection sample */
	
	t1x = t1 = itfrom - m_iT1sel;
	t2 = min(rows(mx) + m_iT1sel - m_iT1est - 1, t1 + cTforc - 1);
	cx = columns(mx);
	
	if (cx == 0 || t2 < t1)						 /* no post-estimation data */
		return <>;

    return mx[t1 : t2][];
}

Modelbase::PrintTestVal(const dTest, const cR, const cTdf, const sLabel)
{
    decl pval;

    if (cR <= 0)
        return;
    if (cTdf > 0)
    {   
        pval = tailf(dTest, cR, cTdf);
        print(sLabel, " F(", cR, ",", cTdf, ") = ", dTest, " [", "%.4f", pval, "]");
    }
    else
    {   pval = tailchi(dTest, cR);
        print(sLabel, " Chi^2(", cR, ") = ", dTest, " [", "%.4f", pval, "]");
    }
    println(pval < 0.01 ? " **" : pval < 0.05 ? " *" : "");
}
Modelbase::TestRestrictions(mR, ...)
{
    decl cr, vr, msel, test, ctdf = 0, as, bsubset, blongrun, bgenres, btestrob,
		testrob, mrtr, vrbr, vp, pval, stype;
    decl args = va_arglist();

	blongrun = bgenres = btestrob = FALSE;
	
    if (sizeof(args) == 0 || sizerc(args[0]) == 0)
    {
        msel = mR;
        mR = unit(sizerc(mR))[vecindex(msel)][]; // mR selects vars for exclusion
        bsubset = TRUE;                          // translate to selection matrix
        cr = sizer(mR);
        vr = zeros(cr, 1);
    }
    else
    {   bsubset = FALSE;
        cr = sizer(mR);
        vr = vec(args[0]);
        if (sizerc(vr) != cr)
        {
            println("Need: rows(mR) == sizerc(vR)");
            return;
        }
    }
    if (sizeof(args) >= 3)
	{
		if (!isstring(args[2]))
			return;
		bgenres = TRUE;
	}
    if (sizeof(args) > 1 && args[1] == 1)
    {
        blongrun = TRUE;
        [as, vp, mrtr] = GetLongRunInfo();
    }
    else
    {	as = GetFreeParNames();
		vp = GetFreePar();
        mrtr = GetCovar();
    }
	if (sizerc(mrtr) == 0)
	{
		println("*** Error: cannot test restrictions");
		return;
	}
    mrtr = invertgen(mR * mrtr * mR', 30);        /* mrtr = R * Var(b) * R' */
	if (bgenres)
		vrbr = vr;
	else
	    vrbr = mR * vp - vr;                               /* vrbr = Rb - r */
    test = double(vrbr' * mrtr * vrbr);     /* (Rb-r)' * inv(mrtr) * (Rb-r) */
        
    if (!blongrun && rows(m_mCovarRobust) && columns(mR) == rows(m_mCovarRobust))
	{
	    mrtr = invertgen(mR * m_mCovarRobust * mR', 30);
		vrbr = bgenres ? vr : mR * vp - vr;
	    testrob = double(vrbr' * mrtr * vrbr);
		btestrob = TRUE;
	}

	decl bdo_f = FALSE;

	if (m_fPrint)
	{
		if (bsubset)
	    {   println("\nTest for excluding", blongrun ? " from the long run" : "",
	            ": ", as[vecindex(msel)]);
			stype = "Subset";
	    }
	    else if (!bgenres)
	    {   println("\nTest for linear restrictions (Rb=r)", blongrun ? " on the long run" : "",
	            ":");
	        println("R matrix", "%c", as, mR, "r vector", vr');
			stype = "LinRes";
	    }
		else
	    {   println("\nTest for general restrictions:");
			println(args[2]);
			stype = "GenRes";
		}
		if (bdo_f)
			test /= cr, ctdf = GetcT() - GetcDfLoss();
	    PrintTestVal(test, cr, ctdf, stype);
		if (btestrob)
		{	println("Using robust standard errors:");
		    PrintTestVal(testrob, cr, ctdf, stype);
		}
	}
	
	if (ctdf > 0)
        pval = tailf(test, cr, ctdf);
    else
    	pval = tailchi(test, cr);

	decl retval = test | pval;
	if (btestrob)
	{
		if (ctdf > 0)
	        pval = tailf(testrob, cr, ctdf);
	    else
	    	pval = tailchi(testrob, cr);
		retval ~= testrob | pval;
	}		
	return retval;
}
/*------------------------ END Modelbase : Database ---------------------------*/

Modelbase::Buffering(const bBufferOn)
{
	// default does nothing
}
Modelbase::LoadOptions()
{
    decl deps1, deps2, iprint, iitmax, bcompact;
    [iitmax, iprint, bcompact] = GetMaxControl();       
    [deps1, deps2] = GetMaxControlEps();        

	iitmax  = "OxPackReadProfileInt"("ModelBase", "itmax", iitmax);
	iprint  = "OxPackReadProfileInt"("ModelBase", "iprint", iprint);
	bcompact= "OxPackReadProfileInt"("ModelBase", "bcompact", bcompact);
	deps1   = "OxPackReadProfileDouble"("ModelBase", "deps1", deps1);
	deps2   = "OxPackReadProfileDouble"("ModelBase", "deps2", deps2);

	MaxControl(iitmax, iprint, bcompact);     
    MaxControlEps(deps1, deps2);
}
Modelbase::SaveOptions()
{
    decl deps1, deps2, iprint, iitmax, bcompact;
    [iitmax, iprint, bcompact] = GetMaxControl();       
    [deps1, deps2] = GetMaxControlEps();        
	"OxPackWriteProfileInt"("ModelBase", "itmax", iitmax);
	"OxPackWriteProfileInt"("ModelBase", "iprint", iprint);
	"OxPackWriteProfileInt"("ModelBase", "bcompact", bcompact);
	"OxPackWriteProfileDouble"("ModelBase", "deps1", deps1);
	"OxPackWriteProfileDouble"("ModelBase", "deps2", deps2);
}

Modelbase::BatchVarStatus(const sTypes, const vcTypes)
{
	return 0;
}
Modelbase::BatchMethod(const sMethod)
{
	return 0;
}
Modelbase::Batch(const sBatch, ...)
{
	if (sBatch == "option")
	{
		decl args = va_arglist(), cargs;

		if (isarray(args))
			args = args[0];
		cargs = sizeof(args);

		if (cargs == 2)
		{	decl sopt = strlwr(args[0]);
			
			if (sopt == "maxit")
		        MaxControl(args[1], -1, -1);
			else if (sopt == "print")
		        MaxControl(-1, args[1], -1);
			else if (sopt == "compact")
		        MaxControl(-1, -1, args[1]);
			else if (sopt == "strong")
		        MaxControlEps(args[1], 0);
			else if (sopt == "weak")
		        MaxControlEps(0, args[1]);
			return 1;
		}
		if (cargs)	
			println("Batch error: unkown argument for option: ", args[0]);
	}
	else
		println("Batch error: unkown command: ", sBatch);

return 0;	
}
Modelbase::GetBatchModelSettings()
{
	return "";
}
