/*--------------------------------------------------------------------------
 * g2ox.ox - functions mapped from Gauss to Ox
 *
 *       (C) Jurgen Doornik 1996-2001
 *
 * Note, this code requires array_base 1:
 */
#pragma array_base(1)
/*--------------------------------------------------------------------------*/

#ifndef G2OXOX_INCLUDED
#define G2OXOX_INCLUDED

#include <oxstd.h>
#include <oxfloat.h>
#include <oxprob.h>
#include <g2ox.h>
#import <maximize>				   /* required for hessp, gradp and qnewton */
#import <solvenle>						 			/* required for eqsolve */

enum
{   RS_BINARY_SKIP, RS_ASCII_OBS, RS_ASCII_VAR, RS_BINARY_VAR,
    RS_PRINT_OBS,   RS_FMT_VAR,   RS_XLS_OBS,   RS_WKS_OBS,
    RS_UNKNOWN,     RS_FAILED,    RS_DHT_DAT,   RS_SIZEMAT,
    RS_MAT,         RS_BSIZEMAT,  RS_BMAT,      RS_FMT,
 	RS_DTA,   		RS_BIN_DAT,   RS_CSV_OBS
};
loaddata(const file, const type, const cvar, const cobs, const fskip,
    const adata, const anames, const afrequency,
    const ayear1, const aperiod1, const ayear2, const aperiod2);

static decl s_detl_value = 0.0;

#ifdef OX_GAUSS
namespace gauss {
#endif

decl _svderr = 1, _svdtol = 1e-13;

__loadd__(sName, const iMode);

static __isscalar__(const x)
{
    return ::sizerc(x) == 1;
}

static __trap1error__(sMsg, const code)
{
	if (trapchk(1))
		return error(code);

	if (!isstring(sMsg))
	{
		if (sMsg <= 40) sMsg = "Matrix not positive definite";
		else if (sMsg <= 50) sMsg = "Matrix singular";
		else if (sMsg <= 60) sMsg = "Maximum number of iterations reached";
		else if (sMsg <= 70) sMsg = "Singular values not all computed";
		else sMsg = "??";
	}
	::oxrunerror("Error: " + sMsg, 1);
	exit(50);
}
static __solpd__(const b, const a, const errorcode)
{
    decl ml, md, iwarn = oxwarning(0), iret;
    iret = ::decldl(a, &ml, &md);
	oxwarning(iwarn);

	if (iret != 1)
		return __trap1error__(40, errorcode);

	::s_detl_value = 1 / double(::prodr(md));
    return ::solveldl(ml, md, b);
}
#ifdef OX_GAUSS
__strcmp__(const s1, const s2);

__opdivide__(const y, const x)
{
	if (::sizerc(x) <= 1 || ::sizerc(y) <= 1)
	{
		return y ./ x;
	}

	decl r = ::sizer(x), c = ::sizec(x);
	if (r == c)											/* solve xw=y for w */
	{
		decl ml, mu, mp, dg, sign, iwarn = oxwarning(0), iret;
		iret = ::declu(x, &ml, &mu, &mp);
		oxwarning(iwarn);
		if (iret == 0)
		{
			return __trap1error__(10, 10);
		}

		dg = ::diagonal(mu);
	
		sign = ::sumr(mp[2][] .!= ::range(0, r-1));
		sign = ::imod(sign, 2) == 1 ? -1 : 1;
		::s_detl_value = double(::prodr(dg)) * sign;
		return ::solvelu(ml, mu, mp, y);
	}
	if (r > c)									  /* solve (x'x)w=x'y for w */
	{
		if (trapchk(2))
		{	decl w = ::isdotnan(w) .? 0 .: w;
		    __solpd__(w'y, w'w, 40);
		}
		return __solpd__(x'y, x'x, 40);
	}

	return y/x;									   /* will produce an error */
}
__opstrarrayhcat__(const x, const y)
{
	/* x and y have already been converted to arrays */
	decl i, r = ::sizeof(x), s = new array[r];
	if (r != ::sizeof(y))
		::println("Runtime error: rows don't match for $~");
	for (i = 1; i <= r; ++i)
		s[i] = array(x[i]) ~ y[i];
	return s;
}
__opstrplus__(const x, const y)
{
	if (::ismatrix(x) || ::isdouble(x) || ::isint(x))
	{
		decl xx = matrix(x), i, rxx = rows(xx), ry = rows(y);
		if (rxx == 1 && ry > 1)
			for (i = ry, xx = matrix(y); i >= 1; --i)
				xx[i] = __opstrplus__(string(x), y[i]);
		else if (ry == 1)
			for (i = rxx; i >= 1; --i)
				xx[i] = __opstrplus__(string(xx[i]), y);
		else
			for (i = rxx; i >= 1; --i)
				xx[i] = __opstrplus__(string(xx[i]), y[i]);
		return xx;
	}

	if (::isstring(x))
	{	if (::isstring(y))
			return x + y;
		if (::isarray(y))
		{
			::println("Runtime error: string $+ string-array");
			return x + y;
		}	
		if (::sizerc(y) == 1)
			return x + string(y);
		decl r = ::sizer(y), c = ::sizec(y), i, j, as = new array[r];
		if (c == 1)
		{	for (i = 1; i <= r; ++i)
				as[i] = x + string(y[i][1]);
		}
		else
		{	for (i = 1; i <= r; ++i)
				for (j = 1, as[i] = new array[c]; j <= c; ++j)
					as[i][j] = x + string(y[i][j]);
		}
		return as;
	}
	return x + y;
}
static __cvt2array(const m)
{
	decl r = ::sizer(m), c = ::sizec(m), i, j, as;
	if (::isstring(m))
		as = {array(m)};
	else if (::isarray(m))
	{
		if (r && ::isarray(m[1]))
			return m;
		for (i = 1, as = new array[c][1]; i <= c; ++i)
			as[i][1] = m[i];
	}
	else
	{	for (i = 1, as = new array[r][c]; i <= r; ++i)
			for (j = 1; j <= c; ++j)
				as[i][j] = string(m[i][j]);
	}
	return as;
}
static __opdollardotbase__(const x, const y)
{
	decl ax = __cvt2array(x), ay = __cvt2array(y);
	decl rx = rows(ax), cx = cols(ax), ry = rows(ay), cy = cols(ay), i, j, m;

//println("ax", ax);	
//println("ay", ay);
	if (rx == 1 && cx == 1)
	{
		for (i = 1, m = new matrix[ry][cy]; i <= ry; ++i)
			for (j = 1; j <= cy; ++j)
				m[i][j] = __strcmp__(ax[1][1], ay[i][j]);
	}
	else if (ry == 1 && cy == 1)
	{
		for (i = 1, m = new matrix[rx][cx]; i <= rx; ++i)
			for (j = 1; j <= cx; ++j)
				m[i][j] = __strcmp__(ax[i][j], ay[1][1]);
	}
	else if (cx == 1 && ry == 1)
	{
		for (i = 1, m = new matrix[rx][cy]; i <= rx; ++i)
			for (j = 1; j <= cy; ++j)
				m[i][j] = __strcmp__(ax[i][1], ay[1][j]);
	}
	else if (cx == 1)
	{
		for (i = 1, m = new matrix[ry][cy]; i <= ry; ++i)
			for (j = 1; j <= cy; ++j)
				m[i][j] = __strcmp__(ax[i][1], ay[i][j]);
	}
	else if (ry == 1)
	{
		for (i = 1, m = new matrix[rx][cx]; i <= rx; ++i)
			for (j = 1; j <= cx; ++j)
				m[i][j] = __strcmp__(ax[i][j], ay[i][1]);
	}
	else
	{
		for (i = 1, m = new matrix[rx][cx]; i <= rx; ++i)
			for (j = 1; j <= cx; ++j)
				m[i][j] = __strcmp__(ax[i][j], ay[i][j]);
	}
	return m;
}
__opdollardot__(const x, const y, const op)
{
	decl m = __opdollardotbase__(x, y);
	if (op == '<')		 return m .< 0;
	else if (op == 'l')	 return m .<= 0;
	else if (op == '=')	 return m .== 0;
	else if (op == 'g')	 return m .>= 0;
	else if (op == '>')	 return m .> 0;
	else if (op == '!')	 return m .!= 0;
}
__opdollar__(const x, const y, const op)
{
	decl m = __opdollardotbase__(x, y);
	if (op == '<')		 return m < 0;
	else if (op == 'l')	 return m <= 0;
	else if (op == '=')	 return m == 0;
	else if (op == 'g')	 return m >= 0;
	else if (op == '>')	 return m > 0;
	else if (op == '!')	 return m != 0;
}
#else
closeall()
{
}
#endif

__printunsup__(const s)
{
    println(s, "() unsupported");
	return {};
}

#ifdef OX_GAUSS
__balance__(const ma, const adg);
balance(const mx)
{
    decl a, s;
    a = __balance__(mx, &s);
    return { a, ::diag(s) };
}
#else
balance(const mx)
{
    __printunsup__("balance");
	return <>;
}
#endif
band(const ma, const n)
{
	return ::diagonal(ma', n)';	 /* ma' to be compatible for non-sym matrix */
}
bandchol(const mb)
{
	decl ml, md, cb = ::columns(mb), cn = ::rows(mb), i, j;

	::decldlband(mb', &ml, &md);

	md = 1 ./ ::sqrt(md);
    for (i = cb - 1, ml[cb][] .*= md, j = 1; i >= 1; --i, ++j)
	{
		/* make Choleski factor; remember: arraybase is 1 here */
		ml[i][] .*= (::zeros(1,j) ~ md[ : cn - j]);
	}
	return ml';
}
bandcholsol(const ma, const mb)
{
	decl ml, md, cb = ::columns(mb), cn = ::rows(mb), i, j;

	ml = mb';
	md = 1 ./ ml[cb][];
    for (i = cb - 1, ml[cb][] = 1, j = 1; i >= 1; --i, ++j)
	{
		/* make LD factors; remember: arraybase is 1 here */
		ml[i][] .*= (::zeros(1,j) ~ md[ : cn - j]);
	}
	return ::solveldlband(ml, sqr(md), ma);
}
bandltsol(const ma, const mb)
{
	return bandsolpd(ma, mb);
}
bandrv(const mb)
{
	decl ma, cb = ::columns(mb), cn = ::rows(mb), i, j, k, m;

	ma = zeros(cn, cn);
    for (i = 1, m = -cb + 1; i <= cn; ++i, m++)
        for (j = 1 + max(0,m), k = j - m; j < i; ++j, ++k)
            ma[i][j] = mb[i][k];
	ma = ::setupper(ma, ma', mb[][cb]);
	
	return ma;
}
bandsolpd(const ma, const mb)
{
	decl ml, md;

	::decldlband(mb', &ml, &md);
	return ::solveldlband(ml, md, ma);
}
base10(const x)
{
    decl mantissa, pow;

	if (x .== 0)
		return {0, 1};
    pow = ::floor(::log10(::fabs(x)));
    mantissa = x ./ 10^pow;

return { mantissa, pow };
}
besselj(const n, const x)
{
    return ::bessel(x, 'J', ::trunc(n));
}
bessely(const n, const x)
{
    return ::bessel(x, 'Y', ::trunc(n));
}

cdfbeta(const x, const a, const b)
{
	return ::probbeta(x, a, b);
}
cdfbvn(const h, const k, const rho)
{
	return ::probbvn(h, k, rho);
}
cdfbvn2(const h, const dh, const k, const dk, const r)
{
    __printunsup__("cdfbvn2");
	return .NaN;
}
cdfbvn2e(const h, const dh, const k, const dk, const r)
{
    __printunsup__("cdfbvn2e");
	return .NaN;
}
cdfchinc(const x, const nu, const d)
{
	return ::probchi(x, nu, ::sqr(d));
}
cdffnc(const a, const m, const n, const d)
{
	return ::probf(a, m, n, ::sqr(d));
}
cdfgam(const x, const r)
{
	return ::probgamma(r, x, ::ones(x));
}
cdfmvn(const x, const r)
{
    __printunsup__("cdfmvn");
	return .NaN;
}
cdfn2(const x, const dx)
{
	return probn(x + dx) - probn(x);
}
cdftci(const p, const n)
{
	return ::quant(1 - p, n);
}
cdftnc(const x, const n, const d)
{
	return ::probt(x, n, d);
}
cdftvn(const x1, const x2, const x3, const rho12, const rho23, const rho31)
{
	decl sigma = ::unit(3);
	sigma[2][1] = sigma[1][2] = rho12;
	sigma[3][1] = sigma[1][3] = rho31;
	sigma[3][2] = sigma[2][3] = rho23;
	return ::probmvn(x1~x2~x3, sigma);
}
#ifndef OX_GAUSS
cdir(const s)
{
    return ::getcwd();
}
#endif
changedir(const s)
{
    ::chdir(s);
    return ::getcwd();
}
chol(const mx)
{
	decl iwarn = oxwarning(0);
	decl p = ::choleski(mx);
	oxwarning(iwarn);
	if (p == 0)
		return __trap1error__(10, 10);
	::s_detl_value = double(::sqr(::prodr(::diagonal(p))));
	return p';                    /* upper triangular in Gauss */
}
choldn(const mp, const mx)
{
	if (::columns(mx) == 1)
		return ::choleski ( (mp'mp - mx*mx') ');
    return ::choleski( (mp'mp - mx'mx)' );
}
cholsol(const mb, mp)
{
    decl vdi;

    mp = mp';
    vdi = 1 ./ ::diagonal(mp);
    return ::solveldl(mp .* vdi, sqr(vdi), mb);
}
cholup(const mp, const mx)
{
	if (::columns(mx) == 1)
		return ::choleski ( (mp'mp + mx*mx') ');
    return ::choleski((mp'mp + mx'mx)');
}
chrs(const mx)
{
    decl s, i, j, k, r = ::rows(mx), c = ::columns(mx);

    if (__isscalar__(mx))
    {
        s = new string[2];
        s[1] = mx;
        s[2] = 0;
        return s;
    }
    s = new string[r * c + 1];

    for (i = k = 1; i <= r; ++i)
        for (j = 1; j <= c; ++j, ++k)
            s[k] = int(mx[i][j]);
    s[k] = 0;

    return s;
}
#ifndef OX_COMPLEX
cmadd(const xr, const xi, const yr, const yi)
{
	return {xr + yr, xi + yi};
}
cmcplx(const x)
{
	return {x, ::zeros(x)};
}
cmcplx2(const x1, const x2)
{
	return {x1, ::zeros(x1), x2, ::zeros(x2)};
}
cmdiv(const xr, const xi, const yr, const yi)
{
	decl r = ::sizer(xr), c = ::sizec(xr), z;
	z = ::cdiv(xr[]' | xi[]', yr[]' | yi[]');
	return {::reshape(z, r, c), ::reshape(z[r * c + 1 :], r, c)};
}
cmemult(const xr, const xi, const yr, const yi)
{
	return {xr .* yr - xi .* yi, xr .* yi + xi .* yr};
}
cmimag(const xr, const xi)
{
	return xi;
}
cminv(const xr, const xi)
{
    decl inv, zr, zi;
    inv = invert(xr) * xi;
    zr = invert(xr + xi * inv);
    zi = -inv * zr;
    return  {zr, zi};
}
cmmult(const xr, const xi, const yr, const yi)
{
	return {xr * yr - xi * yi, xr * yi + xi * yr};
}
cmreal(const xr, const xi)
{
	return xr;
}
cmsoln(const br, const bi, const ar, const ai)
{
    decl r = ::sizer(br), z, ml, mu, mp;

	::declu((ar ~ -ai) | (ai ~ ar), &ml, &mu, &mp);
	z = ::solvelu(ml, mu, mp, br | bi);
    return {z[:r][], z[r + 1:][]};
}
cmsub(const xr, const xi, const yr, const yi)
{
	return {xr - yr, xi - yi};
}
cmtrans(const xr, const xi)
{
	return {xr', -xi'};
}
#endif

code(const me, const v)
{
    decl i, j, k, vy, rv = ::sizer(v), r = ::sizer(me), c = ::sizec(me);

    vy = ::zeros(r, 1);
    for (i = 1; i <= r; ++i)
    {
        for (j = 1, k = rv; j <= c; ++j)
            if (me[i][j]) k = j;
        vy[i][1] = v[k][1];
    }
    return vy;
}
color(const c)
{
	return c;
}
complex(const xr, const xi)
{
    __printunsup__("complex");
    return xr | xi;
}
con(const r, const c)
{
    decl x;
    scan("?", "%#m", r, c, &x);
    return x;
}
cond(const mx)
{
	decl mw = svds(mx), vmin;
    vmin = min(mw);
    return vmin == 0 || ::isnan(mw) ? 1e300 : mw[1];
}
conj(const mx)
{
    __printunsup__("conj");
	return <>;
}
cons()
{
    decl s;
    scan("?", "%s", &s);
    return s;
}
conv(const b, const x, f, el)
{
    decl cv = ::polymul(b', x');

    if (el == 0) el = ::columns(cv);
    if (f == 0) f = 1;

    return cv[][f:el]';
}
coreleft()
{
	return 2^15;
}
corrvc(const mvc)
{
    decl sdi = 1 ./ ::sqrt(::diagonal(mvc));

    return sdi .* mvc .* sdi';
}
corrm(const mm)
{
    decl ct = mm[1][1];
    decl means = mm[2:][1]/ct;

    return corrvc( mm[2:][2:]/ct - means * means');
}
corrx(const mx)
{
	decl z = mx - ::meanc(mx);
    return corrvc(z'z);
}
counts(const vx, const v)
{
    decl x = ::countc(::vec(vx), ::vec(v)');
    return x[1 : ::rows(x) - 1][];      // drop count above upper bound
}
countwts(const x, const v, const w)
{
	decl vx = vec(x), vv = vec(v), vw = vec(w);
	if (::sizec(vx) != 1 || ::sizec(vv) != 1 || ::sizec(vw) != 1 ||
		::sizer(vx) != ::sizer(vw))
	{
		println("countwts dimension error");
		exit(1);
	}
	decl rx = ::sizer(vx), rv = ::sizer(vv), i, count;
	count = ::zeros(rv, 1);
	for (i = 1; i <= rx; ++i)
		count += (vx[i] .<= vv) * vw[i];
	return ::diff0(count, 1);
}
crossprd(const x, const y)
{
    decl v1, v2, v3;
    v1 = x[2][] .* y[3][] - x[3][] .* y[2][];
    v2 = x[3][] .* y[1][] - x[1][] .* y[3][];
    v3 = x[1][] .* y[2][] - x[2][] .* y[1][];
    return v1 | v2 | v3;
}
crout(const mx)
{
	decl ml, mu, dg;
	::declu(mx, &ml, &mu, 0);
	dg = ::diagonal(mu);
	::s_detl_value = double(::prodr(dg));
	// NB: Gauss has diagonal of ones in U, Ox has it in L
	return ::setupper(dg .* ml, mu ./ dg');
}
croutp(const mx)
{
	decl ml, mu, mp, dg, p, q, r = ::sizer(mx), sign;
	::declu(mx, &ml, &mu, &mp);
	dg = ::diagonal(mu);

	sign = ::sumr(mp[2][] .!= ::range(0, r-1));
	sign = ::imod(sign, 2) == 1 ? -1 : 1;
	::s_detl_value = double(::prodr(dg)) * sign;
	
	q = mp[1][] + 1;
	p = ::sortbyc(q' ~ ::range(1, r)', 1)[][2];
	// NB: Gauss has diagonal of ones in U, Ox has it in L
	return (::setupper(dg .* ml, mu ./ dg')[q][]) | p';
}
csrcol()
{
    __printunsup__("csrcol");
	return 1;
}
csrlin()
{
    __printunsup__("csrlin");
	return 1;
}
csrtype(const mx)
{
    return 1;
}
cumprodc(const mx)
{
    return ::cumprod(mx);
}
cumsumc(const mx)
{
    return ::cumulate(mx);
}
cvtos(const mas)
{
	decl s, i;
	for (i = 1, s = ""; i <= ::sizerc(mas); ++i)
		s ~= string(mas[i]);
	return s;
}
#ifdef OX_GAUSS
date()
#else
_g_date()
#endif
{
	decl t = ::timing(::today(), 1);
    return t[ : 3]' | ((t[4] * 60 + t[5]) * 60 + t[6]) * 100;
}
datestr(vt)
{
	decl d, m, y;
	if (vt == 0)
		vt = ::timing(::today(), 1);
	y = vt[1], m = vt[2], d = vt[3];
	y = ::imod(y, 100);
	return sprint("%2d", m, "/", "%02d", d, "/", "%02d", y);
}
datestring(vt)
{
	decl d, m, y;
	if (vt == 0)
		vt = ::timing(::today(), 1);
	y = vt[1], m = vt[2], d = vt[3];
	return sprint("%2d", m, "/", "%02d", d, "/", "%04d", y);
}
datestrymd(vt)
{
	decl d, m, y;
	if (vt == 0)
		vt = ::timing(::today(), 1);
	y = vt[1], m = vt[2], d = vt[3];
	return sprint("%04d", y, "%02d", m, "%02d", d);
}
dayinyr(vt)
{
	if (vt == 0)
		vt = ::timing(::today(), 1);
	return ::dayofcalender(vt[1], vt[2], vt[3]) - ::dayofcalender(vt[1], 1, 1) + 1;
}
delif(const mx, vif)
{
    decl del = ::vecindex(!vif);
    if (::sizerc(del) == 0)
        return M_NAN;

    return mx[del][];
}
design(const mx)
{
    decl i, my, r = ::sizer(mx);

    my = ::zeros(r, ::round(::max(mx[][1])));

    for (i = 1; i <= r; ++i)
        my[i][::round(mx[i][1])] = 1;

    return my;
}
#ifdef OX_GAUSS
det(const mx)
{
	::s_detl_value = determinant(mx);
	return ::s_detl_value;
}
#endif
detl()
{
	return ::s_detl_value;
}
#ifdef OX_GAUSS
#ifdef OX_COMPLEX
dfft(const mx)
{
    return (::dfft(mx[]') ') / sizerc(mx);
}
#else
dfft(const xr, const xi)
{
	decl f;
	if (xi == 0)
		f = (::dfft(xr[]') ') / sizerc(xr);
	else
		f = (::dfft(xr[]' | xi[]') ') / sizerc(xr);
    return {f[][1], f[][2]};
}
#endif
#else
_g_dfft(const mx)
{
    return (::dfft(mx') ') / sizerc(mx);
}
#endif
#ifdef OX_GAUSS
#ifdef OX_COMPLEX
dffti(const mx)
{
    return (::dfft(mx', 1) ') * sizer(mx);
}
#else
dffti(const xr, const xi)
{
	decl f = (::dfft(xr' | xi', 1) ') * sizer(xr);
    return {f[][1], f[][2]};
}
#endif
#else
dffti(const mx)
{
    return (::dfft(mx', 1) ') * sizer(mx);
}
#endif
dfree(const drive)
{
	return 2^31;
}
#ifdef OX_GAUSS
diag(const ma)
{
    return ::diagonal(ma)';
}
#else
diagonalincol(const ma)
{
    return ::diagonal(ma)';
}
#endif
dotfeq(const ma, const mb)
{
    return isdotfeq(ma, mb);
}
dotfge(const ma, const mb)
{
    return isdotfeq(ma, mb) || ma .> mb;
}
dotfgt(const ma, const mb)
{
    return isdotfeq(ma + ::fuzziness(0), mb) || ma .> mb;
}
dotfle(const ma, const mb)
{
    return isdotfeq(ma, mb) || ma .< mb;
}
dotflt(const ma, const mb)
{
    return isdotfeq(ma - ::fuzziness(0), mb) || ma .< mb;
}
dotfne(const ma, const mb)
{
    return !isdotfeq(ma, mb);
}
#ifdef OX_GAUSS
__dos__(const s)
#else
dos(const s)
#endif
{
	::systemcall(s);
}
dstat(sName, vars)
{
	decl x, asx, asx_dbl, r, i;
	if (double(sName) == 0)
	{
		x = vars;
		r = ::sizec(x);
		for (i = 1, asx = new array[r];	i <= r; ++i)
			asx[i] = ::sprint("X", i);
	}
	else
	{	[x, asx] = __loadd__(sName, 2);
		if (!(__isscalar__(vars) && vars == 0))
		{
			decl sel = ::deleter(indcv(vars, asx));
			x = x[][sel];
			asx = array(asx[sel]);
		}
	}
	decl meanx = ::meanc(x)', varx = ::varc(x)', minx = ::minc(x)',
		maxx = ::maxc(x)', misx = ::isdotnan(x);

	if (__output == 1)
		::println("%r", asx, "%c", {"mean", "var", "min", "max", "valid obs"},
			meanx ~ varx ~ minx ~ maxx ~ ::sumc(!misx)');

	for (i = 1, r = sizeof(asx), asx_dbl = new matrix[r]; i <= r; ++i)
		asx_dbl[i] = double(asx[i]);

	return {asx_dbl, meanx, varx, ::sqrt(varx), minx, maxx,
		::sumc(!misx)', ::sumc(misx)'};
}
__dummy__(mx, const v, const bOpen)
{
	decl r, c, k, d, i;
	mx = ::deleter(mx);
    if (::sizerc(mx) == 0 || ::sizec(mx) > 1)
		return <>;

	r = ::sizer(mx);
    c = ::sizec(mx);
	k = ::sizerc(v);
    d = ::zeros(r, k + bOpen);
	
    d[][1] = mx .<= v[1];
    for (i = 2;	i <= k; ++i)
        d[][i] = v[i-1] .< mx .&& mx .<= v[i];
	if (bOpen)
	    d[][k+1] = v[k] .< mx;
    return d;
}
dummy(const mx, const v)
{
	return __dummy__(mx, v, 1);
}
dummybr(const mx, const v)
{
	return __dummy__(mx, v, 0);
}
dummydn(mx, const v, const p)
{
	mx = __dummy__(mx, v, 1);
	mx[1][p] = M_NAN;
	return ::deletec(mx);
}
editm(const mx)
{
    __printunsup__("editm");
	return mx;
}
eig(const mx)
{
    decl ev, iret, iwarn = oxwarning(0);
    iret = ::eigen(mx, &ev);
	oxwarning(iwarn);
	if (iret != 0)
		return __trap1error__(60, 60);
    return ev';
}
eigcg(const mr, const mi)
{
    __printunsup__("eigcg");
	return <>;
}
eigcg2(const mr, const mi)
{
    __printunsup__("eigcg2");
	return <>;
}
eigch(const mr, const mi)
{
    __printunsup__("eigch");
}
eigch2(const mr, const mi)
{
    __printunsup__("eigch2");
	return <>;
}
eigh(const mx)
{
    decl ev, iret, iwarn = oxwarning(0);
    iret = ::eigensym(mx, &ev);
	oxwarning(iwarn);
	if (iret != 0)
		return __trap1error__(60, 60);
    return ::reverser(ev)';
}
eighv(const mx)
{
    decl ev, em, iret, iwarn = oxwarning(0);
    iret = ::eigensym(mx, &ev, &em);
	oxwarning(iwarn);
	if (iret != 0)
		return __trap1error__(60, 60);
    return {::reverser(ev)', ::reverser(em)};
}
eigv(const mx)
{
    decl ev, em, r = ::sizer(mx), iret, iwarn = oxwarning(0);
    iret = ::eigen(mx, &ev, &em);
	oxwarning(iwarn);
	if (iret != 0)
		return __trap1error__(60, 60);
    return {ev', em};
}
eigrg(const mx)
{
    decl ev, iret, iwarn = oxwarning(0);
    iret = ::eigen(mx, &ev);
	oxwarning(iwarn);
	if (iret != 0)
		return __trap1error__(60, 60);
	if (::rows(ev) == 1)
		return {ev', ::zeros(ev)'};
    return {ev[1][]', ev[2][]'};
}
eigrg2(const mx)
{
    decl ev, em, r = ::sizer(mx), iret, iwarn = oxwarning(0);
    iret = ::eigen(mx, &ev, &em);
	oxwarning(iwarn);
	if (iret != 0)
		return __trap1error__(60, 60);
	if (::rows(ev) == 1)
		return {ev', ::zeros(ev)', em, ::zeros(em)};
    return {ev[1][]', ev[2][]', em[ : r][], em[r+1 : ][]};
}
eigrs(const mx)
{
    return eigh(mx);
}
eigrs2(const mx)
{
    return eighv(mx);
}
end()
{
	closeall();
    ::exit(0);
}
envget(const s)
{
    return ::getenv(s);
}
decl __macheps = DBL_EPSILON;
decl _eqs_typicalx = 0;
decl _eqs_typicalf = 0;
decl _eqs_iterinfo = 0;
decl _eqs_jacobianproc = 0;
decl _eqs_maxiters = 100;
decl _eqs_steptol  = DBL_EPSILON^(2/3);
static decl _eqs_solvenlefunc_;

_eqs_func_(const avF, const vX)
{
	avF[1] = _eqs_solvenlefunc_(vX);
	return 1;
}
_eqs_funcjac_(const amJac, const vX)
{
	amJac[1] = _eqs_jacobianproc(vX);
	return 1;
}

eqsolve(const func, const start)
{
	decl x = start, iret;
	
	_eqs_solvenlefunc_ = func;
	iret = SolveNLE(_eqs_func_, &x, ::isfunction(_eqs_jacobianproc) ? 0 : -1,
		_eqs_funcjac_, 1e3 * __tol, 1e4 * __tol, _eqs_maxiters, _eqs_iterinfo);

	if (iret == MAX_CONV || iret == MAX_WEAK_CONV)
		iret = 1;
	else if (iret == MAX_MAXIT)
		iret = 4;
	else if (iret == MAX_LINE_FAIL)
		iret = 5;
	else if (iret == MAX_FUNC_FAIL)
		iret = 7;

	if (__output)
	{
		decl fx = _eqs_solvenlefunc_(x);
		::println("\n========================================================");
		::println("eqsolve: result=", int(iret), " sqrt(f'f/2)=", sqrt(double(fx'fx)/2), " __tol=",__tol);
		::println("%c", {"start","roots","f(roots)"}, start ~ x ~ fx);
		::println("========================================================");
	}
	return {x, iret};
}
eqsolveset()
{
	_eqs_typicalx = 0;
	_eqs_typicalf = 0;
	_eqs_iterinfo = 0;
	_eqs_jacobianproc = 0;
	_eqs_maxiters = 100;
	_eqs_steptol  = DBL_EPSILON^(2/3);
    return {};
}
erfc(const x)
{
    return 1 - erf(x);
}
error(i)
{
	decl s = string(.NaN);
	i = int(i);
	if (i >= 0 && i <= 65535)
	{
#ifdef OX_BIG_ENDIAN
		s[7] = ::idiv(i, 256);
		s[8] = ::imod(i, 256);
#else
		s[1] = ::idiv(i, 256);
		s[2] = ::imod(i, 256);
#endif
	}
	return double(s);
}
etdays(const vt1, const vt2)
{
    return ::dayofcalender(vt2[1], vt2[2], vt2[3]) - ::dayofcalender(vt1[1], vt1[2], vt1[3]);
}
ethsec(const vt1, const vt2)
{
	decl hsecs = etdays(vt1, vt2) * 8640000.0;
	return hsecs + vt2[4] - vt1[4];
}
etstr(hsecs)
{
	decl days, hrs, mins, secs, s;
	days = int(hsecs / 8640000);  hsecs -= days * 8640000.0;
	hrs	 = int(hsecs / 360000);   hsecs -= hrs * 360000;
	mins = int(hsecs / 6000);     hsecs -= mins * 6000;
	
	s = "";
	if (days)
		s ~= ::sprint(days, " days  ");
	if (hrs)
		s ~= ::sprint(hrs, " hours  ");
	if (mins)
		s ~= ::sprint(mins, " minutes  ");
	return s ~ ::sprint("%.2f", hsecs / 100, " seconds");
}
exctsmpl(const infile, const outfile, const percent)
{
    __printunsup__("exctsmpl");
	return <>;
}
exec(const program, const cmdline)
{
	if (sizeof(cmdline))
		::systemcall(program ~ " " ~ cmdline);
	else
		::systemcall(program);
}
export(x, fname, namelist)
{
    return __printunsup__("export");
}
exportf(dataset, fname, namelist)
{
    return __printunsup__("exportf");
}
fge(const ma, const mb)
{
    return isfeq(ma, mb) || ma > mb;
}
fgt(const ma, const mb)
{
    return isfeq(ma + fuzziness(0), mb) || ma > mb;
}
fle(const ma, const mb)
{
    return isfeq(ma, mb) || ma < mb;
}
flt(const ma, const mb)
{
    return isfeq(ma - fuzziness(0), mb) || ma < mb;
}
fne(const ma, const mb)
{
    return !isfeq(ma, mb);
}
#ifdef OX_GAUSS
fft(const mx)
#else
_g_fft(const mx)
#endif
{
	decl f = ::fft(mx[]');
	return f' / ::sizec(f);
}
ffti(const mx)
{
    return ::fft(mx', 1)' * ::sizer(mx);
}
fftm(const mx, const dim)
{
    __printunsup__("fftm");
	return <>;
}
fftmi(const mx, const dim)
{
    __printunsup__("fftmi");
	return <>;
}
fftn(const mx)
{
	return fft(mx);
}
fileinfo(const fspec)
{
    return __printunsup__("fileinfo");
}
files(const s)
{
	decl as = filesa(s), x, sf, idot;

	x = zeros(::sizeof(as), 2);
	for (decl i = 1; i <= ::sizeof(as); ++i)
	{
		sf = as[i];
		idot = strfind(sf, '.');
		if (idot == -1)
			x[i][1] = string(sf);
		else
		{
			x[i][1] = string(sf[ : idot - 1]);
			x[i][2] = string(s[idot : ]);
		}
	}
	return x;
}
#ifndef OX_GAUSS
filesa(const s)
{
    return __printunsup__("filesa");
}
#endif
static __getfilename__(const sPath, const sName, const sExt, const bCheck)
{
	decl ch = sName[1], sfile, dotpos;
	if (ch == '\\' || ch == '/' || ch == '~' || sName[2] == ':')
		sfile = sName;
	else
		sfile = sPath + sName;
	dotpos = ::strfindr(sfile, '.');
	if (dotpos == -1 || ::strfind(sfile[dotpos : ], <'\\','/'>) != -1)
	{
		if (bCheck)		// try if exists without extension
		{
			decl f = ::fopen(sfile);
			if (::isfile(f))
			{
				::fclose(f);
				return sfile;
			}
		}
		sfile += sExt;
	}
return sfile;		
}
static __createdat__(sName, const vnames, const cColumns)
{
	decl file, as = {}, i, cx, cn, sbase = "X";

	if (!::isstring(sName) || sName == "")
		sName = "temp.dat";
	else
		sName = __getfilename__("", sName, ".dat", FALSE);
	file = ::fopen(sName, "wbF");
	if (!::isfile(file))
		return 0;

	if (::ismatrix(vnames) || ::isdouble(vnames))
	{
		cn = ::sizer(vnames);
		for (i = 1, as = new array[cn]; i <= cn; ++i)
			as[i] = string(vnames[i]);
	}
	else if (::isarray(vnames))
		as = vec(vnames);
	else if (::isstring(vnames))
	{
		if (cColumns == 1)
			as = array(vnames);
		else
			sbase = vnames;
	}
	cn = ::sizec(as);
	cx = cColumns;
	decl sfmt = sprint("%", "0", cx == 0 ? 1 : 1 + int(::log10(cx)), "d");
	for (i = cn + 1; i <= cx; ++i)
		as |= sprint(sbase, sfmt, i);

	fwrite(file, as);
	return file;
}	

#ifdef OX_GAUSS
__fclearerror__(const file);
__ferror__(const file);
__fgets__(const file, const n, ...);
__strerror__();

static decl s_FileList, s_cOpenFiles = -1;
const decl MAX_FILES = 12;
close(const ifileno)
{
	if (s_cOpenFiles == -1)
		return;
	if (ifileno <= 0 || ifileno > MAX_FILES)
	{	::println("invalid file handle");
		return;
	}
	if (::isfile(s_FileList[ifileno]))
	{
		::fclose(s_FileList[ifileno]);
		s_FileList[ifileno] = 0;
		--s_cOpenFiles;
	}
	return 0;
}
closeall()
{
	decl i;
	if (s_cOpenFiles == -1)
		return;
	for (i = 1; i <= MAX_FILES; ++i)
	{
		if (::isfile(s_FileList[i]))
		{
			::fclose(s_FileList[i]);
			s_FileList[i] = 0;
		}
	}
	s_cOpenFiles = 0;
}
output(const sfilename, const imode);	// forward declaration
__new__()
{
	closeall();
	output(0, 'f');
}
static __finsertfile__(const file)
{
	decl i;

	if (s_cOpenFiles == -1)					   /* allocate array with files */
	{	s_FileList = new array[MAX_FILES];
		for (i = 1; i <= MAX_FILES; ++i)
			s_FileList[i] = 0;
		s_cOpenFiles = 0;
	}
	
	if (::isfile(file))						     /* successful: find a slot */
	{
		for (i = 1; i <= MAX_FILES; ++i)
		{
			if (!::isfile(s_FileList[i]))					 /* unused slot */
			{
				s_FileList[i] = file;
				++s_cOpenFiles;
				return i;
			}
		}
		::fclose(file);
	}
	return 0;
}
create(const sOptions, const iComplex, const sName, const iWith, const vnames,
	const col, const dtype, const vtype)
{
	return __finsertfile__(__createdat__(sName, vnames, col));
}
open(const sName, const iMode, const iVarindxi)
{
	return __finsertfile__( ::fopen(__getfilename__("", sName, ".dat", FALSE),
		iMode == 0 ? "rbF" : iMode == 1 ? "abF" : "r+bF"));
}
static __finvalidhandle__(const ifileno)
{
	if (ifileno <= 0 || ifileno > MAX_FILES || !::isfile(s_FileList[ifileno]))
	{	::println("invalid file handle");
		return 1;
	}
	return 0;
}
eof(const ifileno)
{
	if (__finvalidhandle__(ifileno))
		return 1;
	return ::feof(s_FileList[ifileno]);
}
fcheckerr(const ifileno)
{
	if (__finvalidhandle__(ifileno))
		return 0;
	return __ferror__(s_FileList[ifileno]) != 0;
}
fclearerr(const ifileno)
{
	if (__finvalidhandle__(ifileno))
		return;
	__fclearerror__(s_FileList[ifileno]);
}
fstrerror()
{
	return __strerror__();	/* NB: doesn't clear error status, unlike Gauss */
}
fopen(const sfilename, const mode)
{
	return __finsertfile__(::fopen(sfilename, mode));
}
fflush(const ifileno)
{
	return ::fflush(s_FileList[ifileno]);
}
fgets(const ifileno, const n)
{
	return __fgets__(s_FileList[ifileno], n);
}
fgetsa(const ifileno, const n)
{
	decl as = {}, i, s;
	for (i = 1; !::feof(s_FileList[ifileno]) && i <= n; ++i)
	{
		::fscan(s_FileList[ifileno], "%z", &s);
		as |= s ~ "\n";
	}
	return as;
}
fgetsat(const ifileno, const n)
{
	decl as = {}, i, s;
	for (i = 1; !::feof(s_FileList[ifileno]) && i <= n; ++i)
	{
		::fscan(s_FileList[ifileno], "%z", &s);
		as |= s;
	}
	return as;
}
fgetst(const ifileno, const n)
{
	return __fgets__(s_FileList[ifileno], n, 0);
}
fputs(const ifileno, const sa)
{
	if (::isarray(sa))
	{
		decl i, c;
		for (i = 1, c = 0; i <= ::sizeof(sa); ++i)
			c += fputs(ifileno, sa[i]);
		return c;
	}
	else if (::isstring(sa))
	{	fprint(s_FileList[ifileno], sa);
		return 1;
	}
	return 0;
}
fputst(const ifileno, const sa)
{
	if (::isarray(sa))
	{
		decl i, c;
		for (i = 1, c = 0; i <= ::sizeof(sa); ++i)
			c += fputst(ifileno, sa[i]);
		return c;
	}
	else if (::isstring(sa))
	{	fprintln(s_FileList[ifileno], sa);
		return 1;
	}
	return 0;
}
fseek(const ifileno, const offset, const base)
{
	return ::fseek(s_FileList[ifileno], 'c', offset + base);
}
ftell(const ifileno)
{
	return ::fseek(s_FileList[ifileno]);
}
colsf(const ifileno)
{
    return ::columns(s_FileList[ifileno]);
}
rowsf(const ifileno)
{
    return ::rows(s_FileList[ifileno]);
}
readr(const ifileno, r)
{
	decl x, cr = rowsf(ifileno), cread;
	if (r > cr) r = cr;
    cread = ::fread(s_FileList[ifileno], &x, 'f', r);
	if (cread <= 0)
		x = <>;
	else if (cread < r)
		x = x[:cread][];
	return x;
}
seekr(const ifileno, const r)
{
    ::fseek(s_FileList[ifileno], 'f', r - 1);
}
writer(const ifileno, const x)
{
    fwrite(s_FileList[ifileno], x);
}
#else
fcheckerr(const ifileno)
{
	return 0;
}
fclearerr(const ifileno)
{
}
fstrerror()
{
	return "";
}
fgets(const file, const n)
{
	decl s;
	fscan(file, "%z", &s);
	return s;
}
fgetst(const file, const n)
{
	decl s;
	fscan(file, "%z", &s);
	if (::sizeof(s) > 0 && s[::sizeof(s)] == "\n")
		return s[: ::sizeof(s) - 1];
	return s;
}
_g_fseek(const file, const offset, const base)
{
	return fseek(file, 'c', offset + base);
}
colsf(const fh)
{
    return columns(fh);
}
rowsf(const f)
{
    return rows(f);
}
readr(const f1, const r)
{
	decl x;
    ::fread(f1, &x, 'f', r);
	return x;
}
seekr(const fh, const r)
{
    fseek(fh, 'f', r - 1);
}
writer(const fh, const x)
{
    fwrite(fh, x);
}
#endif

ftocv(const mx, wid, prec)
{
	decl r = ::sizer(mx), c = ::sizec(mx), y = ::zeros(r, c), i, j, sfmt;

	if (wid > 8)
	{	::println("ftocv warning: width must be <= 8");
		wid = 8;
	}
	if (prec > 8)
	{	::println("ftocv warning: precision must be <= 8");
		prec = 8;
	}
	sfmt = ::sprint("%", "0", wid, ".", prec, "f");
	for (i = 1; i <= r; ++i)
		for (j = 1; j <= c; ++j)
			y[i][j] = double(sprint(sfmt, mx[i][j]));
	return y;
}
ftos(const mx, const fmt, const wid, const prec)
{
	decl sfmt, ipos = strfind(fmt, "%"), s1 = "", s2 = "", ch;
	if (ipos < 0)
	{	return ftos(mx, "%"~fmt, wid, prec);
	}
	if (ipos > 1)
		s1 = fmt[ : ipos - 1];

	sfmt = "%";
	++ipos; ch = fmt[ipos];
	if (ch == '-' || ch == '+')
	{	sfmt ~= string(ch);
		++ipos;  ch = fmt[ipos];
	}
	if (ch == ' ' || ch == '#')
	{	sfmt ~= string(ch);
		++ipos;  ch = fmt[ipos];
	}
	if (ch == '*')
	{	sfmt ~= sprint(int(wid));
		++ipos;  ch = fmt[ipos];
	}
	if (ch == '.')
	{	++ipos;  ch = fmt[ipos];
		if (ch == '*')
		{	sfmt ~= sprint(".", int(prec));
			++ipos;  ch = fmt[ipos];
		}
	}
	if (ch == 'l')
	{	++ipos;  ch = fmt[ipos];
	}
	sfmt ~= ch;
	++ipos;
	if (ipos < sizeof(fmt))
		s2 = fmt[ipos : ];
		
	return sprint(s1, sfmt, double(mx), s2);
}
gammaii(const r, const p)
{
	return ::quangamma(p, r, ::ones(r));
}
decl __output = 1;
decl __tol = 1e-5;
gausset()
{
	__output = 1;
	__tol = 1e-5;
	return {};
}
getf(const filename, const mode)
{
	// strictly: always do binary to keep \r\n, mode 0 is only about Ctrl+Z
	decl file = ::fopen(filename, mode == 0 ? "r" : "rb"), s;
	if (!::isfile(file))
	{
		::println("getf() error: ", filename, " not found");
		return "";
	}
	::fread(file, &s, 's');
	::fclose(file);
	return s;
}
getname(const sName)
{
	return __loadd__(sName, 1);
}
getnr(const nset, const ncols)
{
    __printunsup__("getnr");
	return <>;
}
getpath(const pfname)
{
    __printunsup__("getpath");
	return "";
}

static decl g_jacder1func;
g_jacder1map(const avF, const vP)
{
    avF[1] = g_jacder1func(vP);
    return 1;
}
gradp(const f, const vp)
{
	decl grad;
    g_jacder1func = f;
	NumJacobianEx(g_jacder1map, vp, <>, TRUE, &grad);
	return grad;
}
hasimag(const x)
{
	return FALSE;
}
header(const procname, const dataset, const ver)
{
    __printunsup__("header");
	return "";
}
hess(const x)
{
    __printunsup__("hess");
	return <>;
}

static decl g_numder2func;
g_numder2map(const vP, const adFunc, const avScore, const amHessian)
{
    adFunc[1] = g_numder2func(vP);
    return 1;
}
hessp(const f, const vp)
{
    decl hess;

    g_numder2func = f;
    Num2Derivative(g_numder2map, vp, &hess);
    return hess;
}
imag(const x)
{
	return <>;
}
import(const fname, const range, const sheet)
{
	if (strfind(fname, ".xls") == -1)
	{						   
		println("Warning: this version of import() only supports .xls files");
		return 0;
	}
	decl names,x,i,ifreq,iyear1,iperiod1,iyear2,iperiod2; 
	::loaddata(fname, RS_XLS_OBS, 0, 0, 0, &x, &names, &ifreq, &iyear1, &iperiod1,
		&iyear2, &iperiod2);
	if (sheet>1) println("Warning: read only the first sheet");
	if (isstring(range))
	{
		decl r = ::sizeof(range);
		decl double_point, point_point, before_double_point=0, first_col=0, last_col=0, first_row=0, last_row=0;
		decl range_before, range_after;
	
		double_point=strfind(range, ":");
		point_point=strfind(range, "..");
	
		if (double_point!=-1)
		{
			range_before=range[1:double_point-1];
			range_after=range[double_point+1:];
		}
		if (point_point!=-1)
		{
			range_before=strupr(range[1:point_point-1]);
			range_after=strupr(range[point_point+2:]);
		}
		decl r_a = new matrix[sizeof(range_after)][1], r_b = new matrix[sizeof(range_before)][1];
		for (i = 1; i <= sizeof(range_before); ++i)
		{
			if ((range_before[i]>64)&&(range_before[i]<91))	 // A = 65, Z=90
				first_col=range_before[i]-64+(i-1)*26;
			r_b[i]=range_before[i];
		}
		for (i = 1; i <= sizeof(range_after); ++i)
		{
			if ((range_after[i]>64)&&(range_after[i]<91))	 // A = 65, Z=90
				last_col=range_after[i]-64+(i-1)*26;
			r_a[i]=range_after[i];
		}
		decl first_r=selectifr(r_b, (r_b.<64))-48;
		decl last_r=selectifr(r_a, (r_a.<64))-48;
		for (i=1; i<=rows(first_r);++i)
			first_row+=first_r[i]*10^(rows(first_r)-i);
		for (i=1; i<=rows(last_r);++i)
			last_row+=last_r[i]*10^(rows(last_r)-i);
		if ((first_row>last_row)||(first_col>last_col))
		{
			oxrunerror("Warning: range of import() misspecified", 0);
		}
		if (first_row>1) 
			return { x[first_row-1:last_row-1][first_col:last_col], names };
		else
		{
			println("Warning: names are not recovered in this version of import()");			
			return { x[first_row:last_row-1][first_col:last_col], names };
		}
	}
	else
		return { x, names };
}
importf(const fname, const dataset, const range, const sheet)
{
    return __printunsup__("importf");
}
indcv(const what, const where)
{
#ifdef OX_GAUSS
	decl aswhat, aswhere, r, i, w;
	if (::isstring(what))
		aswhat = array(what);
	else
	{	r = ::isarray(what) ? ::sizeof(what) : ::sizer(what);
		for (i = 1, aswhat = new array[r]; i <= r; ++i)
			aswhat[i] = string(what[i]);
	}
	if (::isstring(where))
		aswhere = array(where);
	else
	{	r = ::isarray(where) ? ::sizeof(where) : ::sizer(where);
		for (i = 1, aswhere = new array[r]; i <= r; ++i)
			aswhere[i] = string(where[i]);
	}
	w = ::strifind(aswhere, aswhat);
	w = w .== -1 .? M_NAN .: w;
	return w';
#else	
	decl w = strifind(array(where), array(what));
	w = w .== -1 .? M_NAN .: w;
	return w';
#endif
}
indexcat(const x, const v)
{
    if (__isscalar__(v))
        return ::vecindex(x .== v);
    else
       return ::vecindex(x .> v[1] .&& x .<= v[2]);
}
indices(const dataset, const vars)
{
	decl all_x, sel_x, asx, i, var = <>, k, vars_double, asx_double;
	[all_x, asx] = __loadd__(dataset, 2);
	asx = array(asx);
	k = ::sizec(all_x);			
	if (k == 0)
		return __trap1error__("indices() no data", 1);
	
	decl rr = ::sizeof(asx);
	decl r = ::sizer(vars);

	for (i = 1, asx_double = new matrix[rr][1]; i <= rr; ++i)
		asx_double[i][] = double(asx[i]);	

	if (isstring(vars))
		vars_double = double(vars);
	else
	{
		for (i = 1, vars_double = new matrix[r][1]; i <= r; ++i)
			vars_double[i][] = double(vars[i]);	
	}

	sel_x = indcv(vars_double, asx_double);
	
	if (::isnan(sel_x) && sumc(::isdotnan(sel_x)) == r)
	{
		if (::maxc(vars_double) > k || ::minc(vars_double) < 1)
		{
			return __trap1error__("indices() out of range", 2);
		}
		var = asx_double[vars_double];
		sel_x = indcv(var, asx_double);
	}						 
	else
	{
		decl dot = ::isdotnan(sel_x);
		if (sumc(dot) < k && sumc(dot) > 0)
		{
//			decl unrecov = selectifr(vars_double, isdotnan(sel_x));
//			for (i = 1; i <= sizerc(unrecov); ++i)
//				println(string(unrecov[i]));	
			return __trap1error__("indices() some variables are not in the dataset", 2);
		}
		var = vars_double;
//		for (i = 1, var = new array[r]; i <= r; ++i)
//			var[i] = string(vars_double[i][]);	
	}
	return { var, sel_x };
}
indices2(const dataset, const var1, const var2)
{
	decl vars1, vars2, ind1, ind2, ret;
	[ vars1, ind1 ] = indices(dataset, var1);
	[ vars2, ind2 ] = indices(dataset, var2);
	return { vars1, ind1, vars2, ind2 };
}
indnv(const what, const where)
{
    decl i, j, f, d, r1 = rows(what), r2 = rows(where);

    for (i = 1, f = constant(.NaN, r1, 1); i <= r1; ++i)
    {
        d = what[i];

        for (j = 1; j <= r2; ++j)
            if (d == where[j])
            {
				f[i] = j;
				break;
			}
    }
    return f;
}
intgrat2(const f, const xl, const gl)
{
    __printunsup__("intgrat2");
	return .NaN;
}
intgrat3(const f, const xl, const gl, const hl)
{
    __printunsup__("intgrat3");
	return .NaN;
}
decl _intord = 12;
intquad1(const f, const xl)
{
    __printunsup__("intquad() unsupported, use quadpack\n");
	return .NaN;
}
intquad2(const f, const xl, const yl)
{
    __printunsup__("intquad2() unsupported, use quadpack\n");
	return .NaN;
}
intquad3(const f, const xl, const yl, const zl)
{
    __printunsup__("intquad3() unsupported, use quadpack\n");
	return .NaN;
}
intrleav(const infile1, const infile2, const outfile, const keyvar, const keytyp)
{
    return __printunsup__("intrleav");
}
__sortchmat__(const x)
{
	decl i, c = ::sizerc(x), asx = new array[c], idx;
	for (i = 1; i <= c; ++i)
		asx[i] = string(x[i]);
	idx = ::sortcindex(asx);
	return x[idx];
}
intrsect(const v1, const v2, const flag)
{
    decl x = ::intersection(v1, v2);
	if (::sizerc(x) == 0)
		return M_NAN;
	if (flag == 0)		// need to sort them alphabetically rtaher the numerically
		x = __sortchmat__(x);
	if (::isnan(v1) && ::isnan(v2))
		return M_NAN | x';
	return x';		
}
intsimp(const f, const xl, const tol)
{
    __printunsup__("intsimp() unsupported, use quadpack");
	return .NaN;
}
#ifdef OX_GAUSS
inv(const mx)
{
	decl logdet, sign, minv;
	if (__isscalar__(mx))
		return 1 / mx;
	decl iwarn = oxwarning(0);
	minv = ::invert(mx, &logdet, &sign);
	oxwarning(iwarn);
	
	if (sign == 0)
		return __trap1error__(50, 50);

	::s_detl_value = ::exp(logdet) * sign;
	return minv;
}
invpd(const mx)
{
	decl logdet, minv, m = norm(mx);
	decl iwarn = oxwarning(0);
	minv = ::invertsym(mx, &logdet);
	oxwarning(iwarn);
	if (minv == 0)
		return __trap1error__(20, 20);

	::s_detl_value = ::exp(logdet);
	return minv;
}
#endif
invswp(const x)
{
    return 1 / x;
}
iscplx(const x)
{
    return 0;
}
iscplxf(const x)
{
    return 0;
}
keyw()
{
	decl c;
	scan("%c", &c);
	return c;
}
key()
{	/* should return 0 if no key available, so not really supported */
	return keyw();
}
lag1(const x)
{
    return ::lag0(x, 1, M_NAN);
}
lagn(const x, const n)
{
    return ::lag0(x, n, M_NAN);
}
lncdfbvn(const h, const k, const rho)
{
	return ::log(::probbvn(h, k, rho));
}
lncdfbvn2(const h, const dh, const k, const dk, const r)
{
    __printunsup__("lncdfbvn2");
	return .NaN;
}
lncdfmvn(const x, const r)
{
    __printunsup__("lncdfmvn");
	return .NaN;
}
lncdfn(const x)
{
	return ::log(::probn(x));
}
lncdfn2(const x, const dx)
{
	return ::log(cdfn2(x, dx));
}
lncdfnc(const x)
{
	return ::log(::tailn(x));
}
lnfact(const mx)
{
    return loggamma(mx + 1);
}
lnpdfmvn(const x, const s)
{
	decl sinv, logdet;
	sinv = ::invertsym(s, &logdet);
    return -0.5*(rows(s) * ::log(M_2PI) + logdet + x'sinv*x);
}
lnpdfn(const x)
{
    return -::log(M_SQRT2PI) - ::sqr(x) / 2;
}
#ifdef OX_GAUSS
static decl s_sLoadsPath__ = "", s_sLoadmPath__ = "", s_sLoadfPath__ = "";
static __setpath__(const asOld, const sPath)
{
	if (::isstring(sPath) && sPath != "")
	{
		decl s = sPath, ch;
		ch = s[sizeof(s)]; 
		if (ch != '\\' && ch != '/')
			s += "/";
		asOld[1] = s;
	}
}
__loadpath__(const iType, const sPath)
{
	if (iType == 'm')
		__setpath__(&s_sLoadmPath__, sPath);
	else if (iType == 's')
		__setpath__(&s_sLoadsPath__, sPath);
}
loadm(const r, const c, const sName)
{
	if (!::isstring(sName) || sName == "")
		return;
	decl sfile = __getfilename__(s_sLoadmPath__, sName, ".fmt", TRUE);
	decl len, x = <>;

	len = ::sizeof(sfile);
	if (len > 3 && ::strlwr(sfile[len - 3 : len]) == ".fmt")
		x = ::loadmat(sfile);
	else if (r == 0 || c == 0)
		return <>;
	else if (r == -2)	// file created with save, no [] used 
		x = loadd(sfile);
	else				// yes: [] or [r,c] follow variable name, is asci
	{
		decl file = ::fopen(sfile, "r");
		if (!::isfile(file))
		{
			::println("Could not open file for loading: ", sfile);
			return <>;
		}
		::fscan(file, "%#m", -1, -1, &x);
		::fclose(file);
		if (r > 0 && c > 0)
			x = ::reshape(x, r, c);
	}
	return x;
}
loads(const r, const c, const sName)
{
	if (!::isstring(sName) || sName == "")
		return;
	decl sfile = __getfilename__(s_sLoadsPath__, sName, ".fst", TRUE);
	decl file = ::fopen(sfile, "rb"), type, dummy, s;
	if (!::isfile(file))
	{
		::println("Could not open file for loading: ", sfile);
		return "";
	}
	::fread(file, &s, 'c', 4);	   /* read first 4 bytes and check for type */
	type = -1;
	if (s[1] == 0xdf && s[2] == 0xdf)					/* small string v89 */
		type = 0;
	else if (s[1] == 0xdf && s[2] == 0xee)			 /* extended string v89 */
		type = 1;
	else if (s[1] == 0 && s[2] == 0)				/* check for string v92 */
	{
		::fread(file, &s, 'c', 4);
		if (s[3] == 0xcf && s[4] == 0xee)					/* small string v89 */
			type = 2;
	}
	if (type == -1)
	{
		::println("Invalid .fst format: ", sfile);
		return "";
	}
	if (type == 2)
		::fread(file, &dummy, 'c', 120);
	else
		::fread(file, &dummy, 'c', 12);
	
	::fread(file, &s, 's');
	::fclose(file);
	return s;
}
#endif
__loadd__(sName, const iMode)
{
	if (!::isstring(sName) || sName == "")
		sName = "temp.dat";
	else
		sName = __getfilename__("", sName, ".dat", TRUE);
	decl file = ::fopen(sName, "rbf"), x, as;
	if (!::isfile(file))
	{
		::oxrunerror(::sprint("Could not open file for loading: ", sName), 1);
		return iMode == 2 ? {<>, 0} : 0;
	}
	::fread(file, &as, 's');
	::fread(file, &x, 'f', ::rows(file));
	::fclose(file);
	if (iMode == 2)
		return {x, as};
	return iMode == 1 ? as : x;
}
loadd(sName)
{
	return __loadd__(sName, 0);
}
#ifdef OX_GAUSS
lower(s)
{
	if (::isstring(s))
		return ::strlwr(s);
	decl r = ::sizer(s), c = ::sizec(s), i, j;
	for (i = 1; i <= r; ++i)
		for (j = 1; j <= c; ++j)
			s[i][j] = double(::strlwr(string(s[i][j])));
	return s;
}
#endif
lowmat1(const mx)
{
    return ::setdiagonal(::lower(mx), 1);
}
ltrisol(const mb, const ml)
{
	// unlike lusol this doesn't check for row-pivoted L
	return ::solvelu(ml, 0, 0, mb);
}
lu(const mx)
{
	decl ml, mu, mp;
	::declu(mx, &ml, &mu, &mp);
	return {ml[mp[1][] + 1][], mu};
}
lusol(const mb, const ml, const mu)
{
	return ::solvelu(ml, mu, -1, mb);
}
makevars(const x, const vnames, const xnames)
{
    return __printunsup__("makevars");
}
#ifdef OX_GAUSS
maxc(const x)
{
    return ::maxc(x)';
}
maxindc(const x)
{
    return ::maxcindex(x)';
}
#endif
maxvec()
{
	return 2^31;
}
mbesselei(const x, const n, const alpha)
{
	decl mx, i;
	for (mx = <>, i = 0; i < n; ++i)
		mx ~= ::bessel(x, "ie", i + alpha);
	return mx;
}
mbesselei0(const x)
{
	return ::bessel(x, "ie", 0);
}
mbesselei1(const x)
{
	return ::bessel(x, "ie", 1);
}
mbesseli(const x, const n, const alpha)
{
	decl mx, i;
	for (mx = <>, i = 0; i < n; ++i)
		mx ~= ::bessel(x, 'i', i + alpha);
	return mx;
}
mbesseli0(const x)
{
	return ::bessel(x, 'i', 0);
}
mbesseli1(const x)
{
	return ::bessel(x, 'i', 1);
}

#ifdef OX_GAUSS
meanc(const x)
{
    return ::meanc(x)';
}
#endif
#ifdef OX_GAUSS
median(const x)
{
    return ::quantilec(x)';
}
#endif
medit(const x, const xv,const xfmt)
{
    return __printunsup__("medit");
}
mergeby(const infile1, const infile2, const outfile, const keytyp)
{
    return __printunsup__("mergeby");
}
mergevar(const vnames)
{
    return __printunsup__("mergevar");
}
#ifdef OX_GAUSS
minc(const x)
{
    return ::minc(x)';
}
minindc(const x)
{
    return ::mincindex(x)';
}
#endif
miss(const x, const v)
{
    return x .== v .? M_NAN .: x;
}
missex(const x, const e)
{
    return e .? M_NAN .: x;
}
missrv(const x, const v)
{
    return ::isdotnan(x) .? v .: x;
}
moment(const a, const b)
{
	if (b == 1)
	{
		decl x = ::deleter(a);
		if (::sizerc(x) == 0)
			return 0;
	    return x'x;
	}
	else if (b == 2)
	{
	    decl x = ::isdotnan(a) .? 0 .: a;
	    return x'x;
	}
    return a'a;
}
momentd(const dataset, const vars)
{
    __printunsup__("momentd");
	return <>;
}
nametype(const vname, const vtype)
{
    return __printunsup__("nametype");
}
nextn(const n0)
{
    __printunsup__("nextn");
	return .NaN;
}
nextnevn(const n0)
{
    __printunsup__("nextnevn");
	return .NaN;
}
ndpchk(const x)
{
    return 0;
}
ndpclex()
{
    return {};
}
ndpcntrl(const n, const m)
{
    return 0;
}
null(const x)
{
    return ::nullspace(x');
}
null1(const x, const dataset)
{
    __printunsup__("null1");
	return <>;
}
static olsOutput(const par, const parnames, const nbpar, const T) 
{ 
	decl i, mx, tval, abstval, tprob, astval, tp;
	print("%29s", "Coefficient", "%11s", "Std.Error"); 
	println("%9s", "t-value", "%8s", "t-prob"); 
	for (i = 1; i < sizer(par)+1; ++i)								 
   	{
       	print("%-17s", string(parnames[i]), "%#12.6f", par[i][1]); 
   		print("%#11.6f", par[i][2]); 
    	tval = par[i][3]; 
		abstval = fabs(tval); 
    	tprob = 2 * tailt(abstval, T-sizer(par));  
   		println(abstval <= 1e-4 ? "%9.2f" : abstval >= 1e3 ? "%#9.4g" 
	   		: "%#9.4g", tval, "%8.4f", tprob); 
	} 
}

decl __con = 1, __miss = 0, __altnam, __row, __olsres;
ols(const dataset, const depvar, const indvars)	  // TO BE CONTINUED
{
	decl Y,X,beta,covmat,XX,e,Omega,Q,k,T,i,SE,sigma2,V,stb,cor,dwstat,
	     Tstat,mpar,names,R2,found_missing=0,parnames={},Yname;
	if (!isstring(dataset))
	{
		Y=depvar; X=indvars;
		if (__con==1)
			X=1~indvars;
		else
			X=indvars;
		for (i=1; i < columns(X)+1; ++i)		parnames |= sprint("X",i);
	}
	else	 
	{
		decl data, asx; 
		[data, asx] = __loadd__(dataset, 2);
		decl col_Y=indcv(depvar, asx);
		decl col_X=indcv(indvars, asx);
		Y=data[][col_Y];
		if (__con==1)
		{
			parnames={};
			parnames|=sprint("CONSTANT");
			X=1~data[][col_X];
		}
		else
		{
			parnames={};
			X=data[][col_X];
		}
		Yname=asx[col_Y];
		for (i=1;i<columns(X);++i)
			parnames|=sprint(asx[col_X[i]]);	
	}
	decl miss=::sumr(isdotnan(Y~X));
	if (::maxc(miss)>0)
		++__miss;
	Y=::selectifr(Y,miss .== 0);
	X=::selectifr(X,miss .== 0);

	k=::columns(X);T=::rows(Y);
	::olsc(Y, X,&beta, &covmat, &XX);		 // covmat = 1/(X'X) and XX = X'X 
	e = Y - X * beta;
	sigma2=(e'*e)/(T-k);
	V=sigma2*covmat;
	SE=sqrt(diagonal(V));
	Tstat= beta'./ SE; 
	mpar = beta ~ SE' ~ Tstat';
	if (__output)
		olsOutput(mpar, parnames, k, T);
	R2=1-(e'*e)/((Y-meanc(Y))'*(Y-meanc(Y)));
	stb=999*ones(k,1);
	dwstat=cor=999;														 
	return { parnames , XX, beta , stb , V , SE' , sqrt(sigma2) , cor , R2 , e , dwstat };
}																						  
olsqr(const y, const x)
{
    decl vb;
    olsc(y, x, &vb);
    return vb;
}
olsqr2(const y, const x)
{
    decl vb;
    olsc(y, x, &vb);
    return {vb, y - x*vb, x*vb};
}
optn(const n0)
{
    __printunsup__("optn");
	return .NaN;
}
optnevn(const n0)
{
    __printunsup__("optnevn");
	return .NaN;
}
orth(x)
{
	decl mht, mr, mp, c;

	if (::sizer(x) < ::sizec(x))
		x = x[][ : ::sizer(x)];  // strange

	::decqr(x, &mht, &mr, &mp);
	mht = ::decqrmul(mht, ::unit(::sizer(x)));
	c = ::sumr(::diagonal(mr) .!= 0);
	return mht[ : c][]';
}
#ifdef OX_GAUSS
static decl s_sLogFilename = "output.out";
static decl s_sLogMode = "";
output(const sfilename, const imode)
{
	if (::isstring(sfilename))
		s_sLogFilename = sfilename;
	else if (imode == 0)
		::eprint("Log file: ", s_sLogFilename, " status: ", s_sLogMode, "\n");
		
	if (imode == 'o')		 /* on */
		::fopen(s_sLogFilename, "la"), s_sLogMode = "on";
    else if (imode == 'f')	 /* off */
		::fclose("l"), s_sLogMode = "off";
    else if (imode == 'r')	 /* reset */
		::fopen(s_sLogFilename, "l"), s_sLogMode = "reset";
}
#endif
packr(const x)
{
    return deleter(x);
}

static __isletter__(const ch)
{
	if ( (ch >= 'a' && ch <= 'z') || (ch >= 'A' && ch <= 'Z') || ch == '_'
		|| (ch >= '0' && ch <= '9'))
		return TRUE;
	return FALSE;		
}
static __isalpha__(const ch)
{
	if ( (ch >= 'a' && ch <= 'z') || (ch >= 'A' && ch <= 'Z') || ch == '_')
		return TRUE;
	return FALSE;		
}
static __isdigit__(const ch)
{
	return (ch >= '0' && ch <= '9') ? TRUE : FALSE;
}
static __NextTokLen__(const s, const sDelim)
{
	decl i = 1, j = 1, len = ::sizeof(s);
		
	while (i <= len)
	{
		while (s[i] == ' ') ++i;
		j = i + 1;

		if ((s[i] >= '0' && s[i] <= '9') || (s[i] == '.' && __isdigit__(s[i + 1])))
		{
			decl c = 0;
			::sscan(s[i:] , "%*g%n", &c);
			j = i + c;
		}
		else if (__isalpha__(s[i]))
			while (__isletter__(s[j])) ++j;
		else
			while (!__isletter__(s[j]) && s[j] != ' ') ++j;
	
		if (j - i > 1 || ::strfind(sDelim, s[i]) == -1)
			break;
		++i;
	}
	return i > len ? 0 : i~j;
}
static __strtok__(s, const sdelim)
{
	decl ipos, as = {};
	ipos = __NextTokLen__(s, sdelim);
	while (ipos[2] > ipos[1])
	{
		as ~= s[ipos[1] : ipos[2] - 1];
		if (ipos[2] >= sizeof(s))
			break;
		s = s[ipos[2] : ];
		ipos = __NextTokLen__(s, sdelim);
	}
	return as;
}
parse(const str, const delim)
{
	decl sdelim, mfound = <>, i, c = ::sizerc(delim), s, astok;

	for (i = 1, sdelim = ""; i <= c; )
	{	s = string(delim[i]);	   /* only use first character of delimiter */
		if (::sizeof(s))
			sdelim ~= s[1], ++i;
		else
			--c;
	}
	astok = __strtok__(str, sdelim ~ " ");
	
	c = ::sizeof(astok);
	if (c == 0)
		return M_NAN;

	for (i = 1, mfound = new matrix[c]; i <= c; ++i)
		mfound[i] = double(astok[i]);
	return mfound;
}
pause(const isec)
{
    decl start = timer();

    do {}
    while ( (timer() - start) / 100 < isec);
	return {};
}
pi()
{
    return M_PI;
}
pinv(const x)
{
	if (__isscalar__(x))
		return 1 / x;
	decl mu, mw, mv, mwi;
	[mu, mw, mv] = svd2(x);
	mwi = ::diagonal(mw);
	mwi = mwi .< _svdtol .? 0 .: 1 ./ mwi;
	if (::sumr(mwi) == 0)
		return __trap1error__("Singular values all zero", 33);
	return (mv .* mwi) * mu';
}
polychar(const x)
{
    decl ev;
    ::eigen(x, &ev);
    return ::polymake(ev)';
}
polyeval(const x, const c)
{
    decl rx = ::sizer(x), cx = ::sizec(x),rc = ::sizer(c), cc = ::sizec(c);
	decl i, p, c1, c1not0, y, yc;

	if (rx == 1 && cx == 1)
	{
		if (rc == 1)
			return c;
		c1not0 = c[1][] .== 0;
		c1 = c1not0 .? 1 .: c[1][];
		y = recsercp(x, (c ./ c1)[2:][]);
        p = rc - 1;
        return (c1 .* y[p][] - x^p .* c1not0)';
	}
    else if (rx > 1 && rx == cx && cc == 1)
	{
		for (i = 2, y = ::setdiagonal(::zeros(rc,rc), c[1]); i <= rc; ++i)
            y = y * x + ::setdiagonal(::zeros(rc,rc), c[i]);
        return y;
	}
    else
        ::println("Error in polyeval()");
}

static __polint__(const xa, const ya, const x)
{
	decl i, m, ns, den, dif, dift, ho, hp, w, c, d, n = ::sizerc(xa), y, dy;

	c = d = ya;
	dif = ::fabs(x - xa);
    ns = ::mincindex(dif[]);
	dif = dif[ns];
	
	for (m = 1, y = ya[ns--]; m < n; ++m, y += dy)
	{
		for (i = 1; i <= n - m; ++i)
		{
			ho = xa[i] - x;
			hp = xa[i + m] - x;
			w = c[i + 1] - d[i];
			den = ho - hp;
			if (den == 0.0)
			{	::println("Error in polint: two xa's are identical");
				return <>;
			}
			den = w / den;
			d[i] = hp * den;
			c[i] = ho * den;
		}
		dy = 2 * ns < (n - m) ? c[ns + 1] : d[ns--];
	}
	return y ~ dy;
}
__polyint__(const xa, const ya, const x, const poldeg, const aPolerr)
{
    decl y, yh, yl, n = ::sizer(xa);

	if (n - 1 <= poldeg)
        yl = 1, yh = n;
    else
	{
	    decl nfor = ::idiv(poldeg, 2);
	    decl nbak = n - nfor;
	    decl ns = ::mincindex(::fabs(x - xa));
		
		if (ns - nbak <= 1)
	        yl = 1, yh = yl + poldeg;
    	else if (ns + nfor >= n)
        	yl = n - poldeg, yh = n;
	    else
    	    yl = ns - nbak, yh = ns + nfor;
	}
	y = __polint__(xa[yl : yh][1], ya[yl : yh][1], x);
	if (y == <>)
		return <-1;-1;-1>;
    aPolerr[1] = y[2];
    return y[1];
}
decl _poldeg = 6;
decl _polerr;
polyint(const xa, const ya, const x)
{
	return __polyint__(xa, ya, x, _poldeg, &_polerr);
}
#ifdef OX_GAUSS
polymake(const x)
{
    return ::polymake(x') ';
}
#endif
polymat(const x, const p)
{
	decl m, i;
	for (i = 2, m = x; i <= p; ++i)
		m ~= x .^ p;
	return m;		
}
polymult(const c1, const c2)
{
    return ::polymul(c1', c2') ';
}
#ifdef OX_GAUSS
polyroot(const x)
{
	decl r, i = 1;
	while (x[i] == 0) ++i;
	::polyroots(x[i : ]', &r);
	if (r[2][] == 0)
		return r[1][]';
    return r';
}
#endif
#ifndef OX_GAUSS
printfm(const x, const mask, const fmt)
{
    ::print(x);
}
printfmt(const x, const mask)
{
    ::print(x);
}
#endif
#ifdef OX_GAUSS
prodc(const x)
{
    return ::prodc(x)';
}
#endif
putf(const f, const str, const start, slen, const mode, const append)
{
	if ((mode != 0 && mode != 1) || str == "")
		return 5;
	if (append != 0 && append != 1)
		return 6;
	if (f == "")
		return 1;
	decl smode = append ? "a" : "w";
	decl file = ::fopen(f, mode == 0 ? smode : smode ~ "b");
	if (!::isfile(file))
		return 2;
	decl len = ::sizeof(str);
	if (slen > len - start) slen = len - start + 1;

	::fwrite(file, str[start : start + slen - 1]);
	::fclose(file);
	return 0;
}

decl _qn_relgradtol = 1e-5;
decl _qn_gradproc;
decl _qn_maxiters = 10000;
decl _qn_printiters;
decl _qn_parnames;
decl _qn_printresults;
static decl _qn_maxbfgsfunc_;

static _qn_func_(const vP, const adFunc, const avScore, const amHessian)
{
	adFunc[1] = -_qn_maxbfgsfunc_(vP);
	if (isarray(avScore))
		avScore[1] = -_qn_gradproc(vP);
	return 1;
}
qnewton(const func, vP)
{
	decl iret, f, vg;
	_qn_maxbfgsfunc_ = func;

	iret = ::MaxBFGS(_qn_func_, &vP, &f, 0, !isfunction(_qn_gradproc),
		_qn_relgradtol, _qn_relgradtol * 100, _qn_maxiters, _qn_printiters);

	if (iret == MAX_CONV || iret == MAX_WEAK_CONV)
		iret = 0;
	else if (iret == MAX_MAXIT)
		iret = 2;
	else if (iret == MAX_LINE_FAIL)
		iret = 5;
	else if (iret == MAX_FUNC_FAIL)
		iret = 3;

	if (isfunction(_qn_gradproc))
		vg = _qn_gradproc(vP);
	else
		Num1Derivative(_qn_func_, vP, &vg);
		
	return {vP, f, vg, iret};
}

qprog(const start, const q, const r, const a, const b, const c, const d, const bnds)
{
    return __printunsup__("qprog");
}
qqr(const x)
{
	decl mht, mr, r = ::sizer(x), c = ::sizec(x);
	::decqr(x, &mht, &mr, 0);
	mht = ::decqrmul(mht, ::unit(r));
	return { mht[ : min(r,c)][]', mr};
}
qqre(const x)
{
	decl mht, mr, mp, r = ::sizer(x), c = ::sizec(x);
	::decqr(x, &mht, &mr, &mp);
	mht = ::decqrmul(mht, ::unit(r));
	return { mht[ : min(r,c)][]', mr, mp[1][]' + 1};
}
qqrep(const x, const pvt)
{
	return qqre(x);
}
qr(const x)
{
	decl mht, mr;
	::decqr(x, &mht, &mr, 0);
	return mr;
}
qre(const x)
{
	decl mht, mr, mp;
	::decqr(x, &mht, &mr, &mp);
	return { mr, mp[1][]' + 1};
}
qrep(const x,const pvt)
{
	return qre(x);
}
qrsol(const b, const r)
{
	return solvelu(0, r, 0, b);
}
qrtsol(const b, const r)
{
	return solvelu(r, 0, 0, b);
}
qtyr(const y, const x)
{
	decl mht, mr, mp;
	::decqr(x, &mht, &mr, 0);
	mht = ::decqrmul(mht, y);
	return {mht, mr};
}
qtyre(const y, const x)
{
	decl mht, mr, mp;
	::decqr(x, &mht, &mr, &mp);
	mht = ::decqrmul(mht, y);
	return {mht, mr, mp[1][]' + 1};
}
qtyrep(const y, const x, const pvt)
{
    return qtyre(y, x);
}
quantile(const x, const e)
{
	decl rx = ::sizer(x), cx = ::sizec(x), re = ::sizer(e), m, iw, z, w, i;

	if (min(e) < 0 || max(e) > 1)
		::println("Error in quantile: need 0 <= level <= 1");
	z = rx * e;
    iw = ::floor(z);
    w = z - iw;

    for (i = 1, m = ::zeros(re, cx); i <= cx; ++i)
	{
        z = ::sortc(x[][i]);
        m[][i] = z[iw] + w .* (z[iw + 1] - z[iw]);
	}
	return m;
//	return ::quantilec(x, e');	better, because it interpolates
}
quantiled(const dataset, const x, const e)
{
    __printunsup__("quantiled");
	return <>;
}
qyr(const y, const x)
{
	decl mht, mr;
	::decqr(x, &mht, &mr, 0);
	mht = ::decqrmul(::reversec(mht), y);
	return {mht, mr};
}
qyre(const y, const x)
{
	decl mht, mr, mp;
	::decqr(x, &mht, &mr, &mp);
	mht = ::decqrmul(::reversec(mht), y);
	return {mht, mr, mp[1][]' + 1};
}
qyrep(const y, const x, const pvt)
{
    return qyre(y, x);
}
#ifdef OX_GAUSS
rank(const x)
{
    decl mw = svd(x);
    return ::isnan(mw) ? 0 : ::sumc(::fabs(mw) .> _svdtol);
}
#endif
rankindx(x,flag)
{
	// NB seems that flag is ignored,
	// otherwise use sortindc instead of sortind when flag=0
	// also note: using sortind to handle missing values
	return ::sortbyc(::range(1,::sizer(x))' ~ sortind(x), 2)[][1];
}
real(const x)
{
	return x;
}
recode(const x, const e, const v)
{
	return x .* !(::maxc(e') ') + e * v;
}
recserar(const x, const x0, const a)
{
    return ::cumsum(x, a', x0);
}
recsercp(x, z)
{
    z = matrix(z);
    if (__isscalar__(x))
        x *= ::ones(::rows(z), ::columns(z));

    return ::cumprod(x, 1, z);
}
recserrc(const x, const z)
{
	decl i, y, xv = x[]', r = ::sizer(z), t0;

	y = ::zeros(z);
	y[1][] = ::fmod(int(xv), int(z[1][]));
	t0 = ::trunc(xv ./ z[1][]);
	
	for (i = 2; i <= r; ++i)
	{	y[i][] = ::fmod(t0, int(z[i][]));
		t0 = ::trunc(t0 ./ z[i][]);
	}		
	return y;
}
#ifdef OX_GAUSS
reshape(const ma, const r, const c)
{
	if (::isarray(ma))
	{
		decl i, j, k, vma, cma, x;
		vma = vec(ma);	// don't use ::vec here
		cma = ::sizeof(vma);
		for (i = k = 1, x = new array[r][c]; i <= r; ++i)
		{
			for (j = 1; j <= c; ++j, ++k)
			{
				if (k > cma)
					k = 1;
				x[i][j] = vma[k];
			}
		}
		return x;
	}
	else if (::isstring(ma))
	{
		decl i, j, x;
		for (i = 1, x = new array[r][c]; i <= r; ++i)
			for (j = 1; j <= c; ++j)
				x[i][j] = ma;
		return x;
	}
	else
		return ::reshape(ma, r, c);
}
#endif
rfft(const mx)
{
	decl f = ::fft(mx[]');
	return f' / ::sizec(f);
}
rffti(const mx)
{
    return ::fft(mx', 1)' * ::sizer(mx);
}
rfftip(const mx)
{
    return rffti(mx);
}
rfftn(const mx)
{
    return rfft(mx);
}
rfftnp(const mx)
{
    return rfftn(mx);
}
rfftp(const mx)
{
    return rfft(mx);
}
#ifdef OX_GAUSS
rndbeta(const r, const c, const a, const b)
{
	return ::ranbeta(r, c, a, b);
}
#endif
rndgam(const r, const c, const alpha)
{
	return ::rangamma(r, c, alpha, ::ones(alpha));
}
#ifdef OX_GAUSS
rndnb(const r, const c, const k, const p)
{
	return ::rannegbin(r, c, k, p);
}
#endif
rndns(const r, const c, const s)
{
    ::ranseed(max(s,1));
    return ::rann(r, c);
}
rndp(const r, const c, const lambda)
{
	return ::ranpoisson(r, c, lambda);
}
rndseed(const s)
{
    ::ranseed(max(s,1));
}
rndus(const r, const c, const s)
{
    ::ranseed(max(s,1));
    return ::ranu(r, c);
}
rndvm(const r, const c, const m, const k)
{
	return ::fmod(::ranmises(r, c, k) + m, M_2PI);
}
rotater(const x, const c)
{
	if (::sizer(c) == 1)
	{
		decl cx = ::sizec(x);
		decl cc = imod(c[1][1], cx);
		if (cc == 0)
			return x;
		else if (cc > 0)
			return (x[][cx - cc + 1 : ] ~ x[][ : cx - cc]);
		else
			return (x[][1 - cc : ] ~ x[][ : -cc]);
	}
	else
	{	decl i, y;
		for (i = ::sizer(x), y = x; i >= 1; --i)
			y[i][] = rotater(y[i][], c[i]);
		return y;
	}
}
rref(x)
{
    decl i, j, k, pr, aj, m = ::sizer(x), n = ::sizec(x);
	decl tol = max(m,n) * 2.24e-16 * norm(x);
	
	for (i = j = 1; i <= m && j <= n; ++j)
	{
        aj = ::fabs(x[i : m][j]);
        k = ::maxcindex(aj) + i - 1;
        if (max(aj) <= tol)
            x[i : m][j] = ::zeros(m - i + 1, 1);
        else
		{
            x[i~k][j : ] = x[k~i][j : ];
            pr = x[i][j : ] / x[i][j];
            x[][j : ] -= x[][j] .* pr;
            x[i][j : ] = pr;
            ++i;
        }
    }
    return x;
}
#ifdef OX_GAUSS
static decl s_sSavePath__ = "";
save(const sOption, const sPath, const sName, const data)
{
	decl sfile, ch, type = 2;

	if (sOption == "/v89")
		type = 0;

	__setpath__(&s_sSavePath__, sPath);

	if (!::isstring(sName) || sName == "")
		return;

	if (::isstring(data))				  /* save as extended v89 or as v96 */
	{
		sfile = __getfilename__(s_sSavePath__, sName, ".fst", FALSE);
		decl file = ::fopen(sfile, "wb"), len, lendiv;
		if (!::isfile(file))
		{
			::println("Could not open file for saving: ", sfile);
			return;
		}
		len = ::sizeof(data);
		lendiv = ::idiv(len + 1, 8);
		::fprint(file, "%c", 0xdf, "%c", 0xee, "%c", 0x1d, "%c", 0x00);
		::fwrite(file, 1);
		::fwrite(file, len + 1);
		::fwrite(file, (lendiv + 1) * 8);		
	    ::fwrite(file, data);
		lendiv = (lendiv + 1) * 8 - len;
		for (; lendiv > 0; --lendiv)
			::fprint(file, "%c", 0);		
	    ::fclose(file);
	}
	else
	{
		sfile = __getfilename__(s_sSavePath__, sName, ".fmt", FALSE);
		savemat(sfile, data, type == 2 ? 1 : 0);
	}		
	//::println("Saved to: ", sfile);
}
#endif
saved(const x, sName, const vnames)
{
	decl file = __createdat__(sName, vnames, ::sizec(x));
	if (!::isfile(file))
		return 0;

	fwrite(file, x);
	fclose(file);
return 1;
}
scalerr(d)
{
	if (!__isscalar__(d))
		return 0;
	if (::ismatrix(d) && d == <>)
		return 65535;
	decl s = string(d), i0, i1;
#ifdef OX_BIG_ENDIAN
	i0 = s[7], i1 = s[8];
#else
	i0 = s[1], i1 = s[2];
#endif
	if (i0 < 0) i0 = 256 + i0;
	if (i1 < 0) i1 = 256 + i1;
	return i0 * 256 + i1;
}
scalmiss(const x)
{
    return __isscalar__(x) ? ::isnan(x) : 0;
}
schtoc(const sch, const trans)
{
    __printunsup__("schtoc");
	return <>;
}
schur(const x)
{
    __printunsup__("schur");
	return <>;
}
#ifdef OX_GAUSS
__screen__(s)
{
	s = ::strlwr(s);
	if (s == "off")
		::oxprintlevel(0);
	else if (s == "on")
		::oxprintlevel(1);
}
#endif
selif(const x,e)
{
    return delif(x, !e);
}
seqa(const start, const inc, const m)
{
	decl st = start, in = inc, mm = m;
	if (::ismatrix(st)) st = double(st);
	if (::ismatrix(in)) in = double(in);
	if (::ismatrix(mm)) mm = int(mm);
	if (in == 0)
		return ::constant(st, mm, 1); 
	else
		return st + ::range(0, in * (mm - 1), in)';
}

seqm(start, const inc, const m)
{
    decl i, res = ::zeros(m,1);

    for (i = 1; i <= m; ++i, start *= inc)
        res[i][1] = start;
    return res;
}
setcnvrt(const type, const ans)
{
}
setdif(const v1, const v2, flag)
{
	if (flag == 1)
		return ::exclusion(v1, v2)';
	else if (flag == 0)
		return sortcc(::exclusion(v1, v2), 1)';
	else if (flag == 2)	// case insensitive
		return setdif(upper(v1), upper(v2), 0);
}
setvars(const dataset)
{
    return __printunsup__("setvars");
}
setvmode(const x)
{
    return 0;
}
shiftr(const x, const c, const d)
{
	if (::sizer(c) == 1 && ::sizer(d) == 1)
	{
		decl cx = ::sizec(x), rx = ::sizer(x), cc = c[1][1];
		if (cc == 0)
			return x;
		else if (::fabs(cc) >= cx)
			return ::constant(d, rx, cx);
		else if (cc > 0)
			return (::constant(d, rx, cc) ~ x[][ : cx - cc]);
		else
			return (x[][1 - cc : ] ~ ::constant(d, rx, cx + cc - 1));
	}
	else
	{	decl i, y, vc = c, vd = d, rx = ::sizer(x), rc = ::sizer(vc), rd = ::sizer(vd);
		if (rx == 1 && (rc > rx || rd > rx) )
		{
			rx = max(rc, rd);
			y = reshape(x, rx, ::sizec(x)); 
		}
		else
			y = x;
		if (rc == 1) vc = ::constant(vc, rx, 1);
		if (rd == 1) vd = ::constant(vd, rx, 1);
		for (i = rx; i >= 1; --i)
			y[i][] = shiftr(y[i][], vc[i], vd[i]);
		return y;
	}
}
sleep(const isecs)
{
	return pause(isecs);
}
solpd(const b, const a)
{
    return __solpd__(b, a, 30);
}
#ifdef OX_GAUSS
sortc(const x, const c)
{
	if (c == 1 && ::sizec(x) == 1)
		return ::sortc(x);
	return ::sortbyc(x, c);
}
#endif
sortcc(const x, const c)
{
	if (::ismatrix(x))
		return x[sortindc(x[][c])][];
	else if (::isarray(x))
	{
		if (::sizeof(x) && ::isarray(x[1]))
			return x[::sortcindex(x[][c])][];
		else
			return ::sortc(x);
	}
}
sortd(const infile, const outfile, const keyvar, const keytyp)
{
    return __printunsup__("sortd");
}
sorthc(const x, const c)
{
	if (c == 1 && ::sizec(x) == 1)
		return ::sortc(x);
	return ::sortbyc(x, c);
}
sorthcc(const x, const c)
{
    return sortcc(x, c);
}
sortind(const x)
{
	return ::sortcindex(x);
}
sortindc(const x)
{
	if (::isarray(x))
	{
		if (::sizeof(x) && ::isarray(x[1]))
			return ::sortcindex(x[][1]);
		else
			return ::sortcindex(x);
	}
	decl i, as, c = ::sizer(x);
	for (i = 1, as = new array[c]; i <= c; ++i)
		as[i] = string(x[i]);
    return ::sortcindex(as);
}
#ifdef OX_GAUSS
sortmc(const x, const vc)
{
	// don't understand the Gauss implementation, which seems wrong
	return ::sortbyc(x, vc);
}
#endif
spline1d(const x, const y, const d, const s, const sigma, const g)
{
    __printunsup__("spline1d");
	return <>;
}
spline2d(const x, const y, const z, const sigma, const g)
{
    __printunsup__("spline2d");
	return <>;
}
sqpsolve(const func, const start)
{
    __printunsup__("sqpsolve");
	return <>;
}
stdc(const x)
{
    return ::sqrt(::varc(x) * ::sizer(x) / (::sizer(x) - 1))';
}
stocv(const s)
{
    decl len = sizeof(s), mc = <>, i;
	for (i = 1; i <= len; i += 8)
		mc |= double(s[i : min(i + 7, len)]);
	return mc;
}
stof(s)			   			    /* currently no support for complex numbers */
{
    if (::isstring(s))
	{
		decl x = <>, d;
		while (::sscan(&s, "%g", &d) == 1)
			x |= d;
		return x;
	}
	decl r = ::sizer(s), c = ::sizec(s), i, j, d;
	for (i = 1; i <= r; ++i)
		for (j = 1; j <= c; ++j)
		{
			d = 0;
			if (::sscan(string(s[i][j]), "%g", &d) == 1)
				s[i][j] = d;
		}
	return s;
}
stop()
{
    ::exit(0);
}
strindx(where, what, const start)
{
	if (!::isstring(what))
		what = string(what[1][1]);
	if (!::isstring(where))
		where = string(where[1][1]);
	decl ipos = ::strfind(where[start : ], what);
	return ipos == -1 ? 0 : ipos + start - 1;
}
strlen(const s)
{
	if (::isstring(s))
		return ::sizeof(s);
	decl r = ::sizer(s), c = ::sizec(s), i, j, len = ::zeros(r,c);
	for (i = 1; i <= r; ++i)
		for (j = 1; j <= c; ++j)
			len[i][j] = ::sizeof(string(s[i][j]));

	return len;
}
strput(substr, str, const offs)
{
	if (!::isstring(substr))
		substr = string(substr[1][1]);
	if (!::isstring(str))
		str = string(str[1][1]);
	if (::sizeof(str) - offs < ::sizeof(substr))
		return str[ : offs - 1] ~ substr;
	str[offs:] = substr;
    return str;
}
strrindx(where, what, const start)
{
	if (!::isstring(what))
		what = string(what[1][1]);
	if (!::isstring(where))
		where = string(where[1][1]);
	decl ipos = ::strfindr(where[start : ], what);
	return ipos == -1 ? 0 : ipos + start - 1;
}
strsect(str, const start, const len)
{
	if (!::isstring(str))
		str = string(str[1][1]);
		
	if (len <= 0 || ::sizeof(str) - start < 0)
		return "";
	if (::sizeof(str) - start < len)
		return str[start : ];
    return str[start : start + len - 1];
}
#ifdef OX_GAUSS
submat(const x, const r, const c)
#else
_g_submat(const x, const r, const c)
#endif
{
    if (r == 0 && c == 0)
        return x;
    if (r == 0)
        return x[][c];
    if (c == 0)
        return x[r][];
    return x[r][c];
}
subscat(const x, v, s)
{
	if (__isscalar__(v))
		return x .== v .? s .: x;
    decl idx = sumr(x .> v'), w, bnan = isnan(v[1]), dmis;
	if (bnan) v = v[2:], dmis = s[1], s = s[2:];
	w = selectrc(s' .* ::ones(x) ~ x, ::range(1,::sizer(x))-1, idx)';
	if (bnan) w = ::isdotnan(x) .? dmis .: w;
	return w;
}
substute(const x, const v, s)
{
	if (::isstring(s)) s = double(s);
    return v .== 0 .? x .: s;
}
#ifdef OX_GAUSS
sumc(const x)
{
    return ::sumc(x)';
}
#endif
static decl s_svdtrap__ = 1;
svd(const mx)
{
    decl mu, mw, iwarn = oxwarning(0), iret, svdtrap = s_svdtrap__;
	s_svdtrap__ = 1;
    _svderr = iret = ::decsvd(mx, &mu, &mw);
	oxwarning(iwarn);
	if (iret != 0)
	{	_svderr = error(_svderr);
		return svdtrap ? __trap1error__(70, iret) : _svderr;
	}
    return mw';
}
svd2(const mx)
{
    decl mu, mw, mv, iwarn = oxwarning(0), iret, svdtrap;
	s_svdtrap__ = 1;
    _svderr = iret = ::decsvd(mx, &mu, &mw, &mv);
	oxwarning(iwarn);
	if (iret != 0)
	{	_svderr = error(_svderr);
		return svdtrap ? __trap1error__(70, iret) : _svderr;
	}
    return {mu, ::diag(mw), mv};
}
svd1(const mx)
{
	decl rx = ::sizer(mx), cx = ::sizec(mx);
	if (rx > cx)
	{	decl mw = mx ~ ::zeros(rx, rx - cx), mu, mv;
		[mu, mw, mv] = svd2(mw);
    	return {mu, mw[:cx][:cx], mv[:cx][:cx]};
	}
	return svd2(mx);
}
svdcusv(const mx)
{
	s_svdtrap__ = 0;
	return svd2(mx);
}
svds(const mx)
{
	s_svdtrap__ = 0;
	return svd(mx);
}
svdusv(const mx)
{
	s_svdtrap__ = 0;
	return svd1(mx);
}
sysstate(const vcase, const y)
{
#ifdef OX_GAUSS
	if (vcase == 4)
	{
		decl s = s_sSavePath__;
		if (::isstring(y))
			__setpath__(&s_sSavePath__, y);
		return s;
	}
	else if (vcase == 5)
	{
		decl s = s_sLoadmPath__;
		if (::isstring(y))
			__setpath__(&s_sLoadmPath__, y);
		return s;
	}
	else if (vcase == 6)
	{
		decl s = s_sLoadfPath__;
		if (::isstring(y))
			__setpath__(&s_sLoadfPath__, y);
		return s;
	}
	else if (vcase == 7)
	{
		decl s = s_sLoadsPath__;
		if (::isstring(y))
			__setpath__(&s_sLoadsPath__, y);
		return s;
	}
#endif	
}
system(const n)
{
    ::exit(n);
}
tab(const col)
{
    return {};
}
tempname(path, pre, ext)
{
	if (!::isstring(ext))
		ext = "";
	else if (ext[1] == '.')
		ext = ext[2 : ];
	if (!::isstring(pre) || pre == "") pre = "TMP";
	if (!::isstring(path))
		path = ::getcwd();
	
	decl file, filename, i;
	i = ::sizeof(path) ? path[::sizeof(path)] : 0;
	if (i && i != '\\' && i != '/')
		path += "/";
	for (i = 0; i < 1000; ++i)
	{
		filename = ::sprint(path, pre, "%05d", i, ".", ext);
		file = ::fopen(filename, "r");
		if (!::isfile(file))
		{
			file = ::fopen(filename, "a");
			if (!::isfile(file))
				return "";
			::fclose(file);
			::fremove(filename);
			return filename;
		}
		::fclose(file);
	}
	return "";
}
#ifdef OX_GAUSS
time()
#else
_g_time()
#endif
{
	decl t = ::timing(::today(), 1);
    return t[4 : ]' | 0;	 /* doesn;t do milliseconds! */
}
timestr(const x)
{
	decl s = x;
#ifdef OX_GAUSS
	if (s == 0) s = time();
#else
	if (s == 0) s = _g_time();
#endif
    return sprint("%2d", s[1], ":", "%02d", s[2], ":", "%02d", s[3]);
}
tocart(const xy)
{
    __printunsup__("tocart");
	return <>;
}
token(const s)
{
	decl i = 1, i0, len = ::sizeof(s);
	while (s[i] == ',' || s[i] == ' ') ++i;
	i0 = i;
	while (i <= len && s[i] != ',' && s[i] != ' ') ++i;
	return {i > i0 ? s[i0 : i - 1] : "", i <= len ? s[i : ] : ""};
}
topolar(const xy)
{
    __printunsup__("topolar");
	return <>;
}
static decl __trap__ = 0;
trap(const trpval, ...)
{
	decl va = va_arglist();
	if (::sizeof(va) == 0)
		__trap__ = trpval;
	else
		__trap__ = ::binor(::binand(__trap__, ::bincomp(int(va[1]))), trpval);
}
trapchk(const m)
{
	return ::binand(__trap__, int(m));
}
trim(const x, const top, const bot)
{
    return trimr(x, top, bot);
}
trimr(const x, const top, const bot)
{
	return x[top + 1 : ::sizer(x) - bot][];
}
type(const x)
{
    if (::isstring(x))
        return 13;
    else if (::ismatrix(x))
        return 6;
    else if (::isfunction(x))
        return 8;
    else if (::isarray(x))
        return 15;
return 0;
}
#ifdef OX_GAUSS
extern __typecv__(const s);
typecv(const x)
{
    if (::isstring(x))
		return __typecv__(x);
	decl i, ret;
	for (i = 1, ret = <>; i <= ::sizerc(x); ++i)
		ret |= __typecv__(string(x[i]));
	return ret;
}
#else
typecv(const x)
{
    __printunsup__("typecv");
	return <>;
}
#endif
typef(const x)
{
    __printunsup__("typef");
	return .NaN;
}
union(const v1, const v2, const flag)
{
    decl x = ::unique(v1 | v2);
	if (flag == 0)		// need to sort them alphabetically
		x = __sortchmat__(x);	
	if (::isnan(v1) || ::isnan(v2))
		return M_NAN | x';
	return x';
}
uniqindx(const v1, const flag)
{
	decl x = union(v1, <>, flag);
	return ::vecrindex(v1, x, 1);
}
#ifdef OX_GAUSS
unique(const v1, const flag)
{
   	return union(v1, <>, flag);
}
#else
g_unique(const v1, const flag)
{
    return ::unique(v1);
}
#endif
upmat1(const mx)
{
    return ::setdiagonal(::upper(mx), 1);
}
#ifdef OX_GAUSS
upper(const s)
{
	if (::isstring(s))
		return ::strupr(s);
	decl r = ::sizer(s), c = ::sizec(s), i, j, w = s;
	for (i = 1; i <= r; ++i)
		for (j = 1; j <= c; ++j)
			w[i][j] = double(::strupr(string(s[i][j])));
	return w;
}
#endif
utrisol(const b, const u)
{
	return solvelu(0, u, 0, b);
}
vals(s)
{
	if (!::isstring(s))
		s = string(s[1][1]);
    decl r = ::sizeof(s), x = new matrix[r][1], i;
	for (i = 1; i <= r; ++i)
		x[i] = s[i];
	return x;	
}
#ifndef OX_GAUSS
varget(const s)
{
    __printunsup__("varget");
	return <>;
}
varput(const x, const n)
{
    return __printunsup__("varput");
}
#endif
vargetl(const s)
{
    return __printunsup__("vargetl");
}
varputl(const x, const n)
{
    return __printunsup__("varputl");
}
vartypef(const names)
{
    return __printunsup__("vartype");
}
vartype(const names)
{
    return vartypef(names);
}
vcm(const m)
{
    decl ct = m[1][1], mu = m[2:][1];
    return (m[2:][2:] - mu*mu' / ct) / (ct - 1);
}
vcx(const m)
{
    decl ct = ::sizer(m);
    return ::variance(m) * ct / (ct - 1);
}
#ifdef OX_GAUSS
vec(const ma)
{
	if (::isarray(ma))
	{
		decl r = ::sizeof(ma), c, i, j, k, x;
		if (!r || !::isarray(ma[1]))
			return ma;
		c = ::sizeof(ma[1]);
		for (i = k = 1, x = new array[r * c]; i <= c; ++i)
			for (j = 1; j <= r; ++j, ++k)
				x[k] = ma[j][i];
		return x;
	}
	else if (::isstring(ma))
		return ma;
	else
		return ::vec(ma);
}
vecr(const ma)
{
	if (::isarray(ma))
	{
		decl r = ::sizeof(ma), c, i, j, k, x;
		if (!r || !::isarray(ma[1]))
			return ma;
		c = ::sizeof(ma[1]);
		for (i = k = 1, x = new array[r * c]; i <= r; ++i)
			for (j = 1; j <= c; ++j, ++k)
				x[k] = ma[i][j];
		return x;
	}
	else if (::isstring(ma))
		return ma;
	else
		return ::vecr(ma);
}
#endif
vget(const dbuf, const name)
{
    return __printunsup__("vget");
}
vlist(const dbuf)
{
    return __printunsup__("vlist");
}
vnamecv(const dbuf)
{
    return __printunsup__("vnamecv");
}
vput(const dbuf, const x, const name)
{
    return __printunsup__("vput");
}
vread(const dbuf, const xname)
{
    return __printunsup__("vread");
}
vtypecv(const dbuf)
{
    return __printunsup__("vtypecv");
}
wait()
{
	decl i;
    scan("%c", &i);
	return "";
}
waitc()
{
	return wait();		// not quite the same: waitc clears keyboard first
}
xpnd(const v)
{
	return ::reflect(::unvech(::reversec(::vec(v))));
//	return ::unvech(v);
}


#ifdef OX_GAUSS
} /* end of gauss namespace */
#endif

#endif /* G2OXOX_INCLUDED */
