/*------------------------------------------------------------------------------*
 * File Name: OCN_g12.h															*
 * Creation: SDB 5/22/2001														*
 * Purpose: Origin C Header for NAG functions									*
 * Copyright (c) OriginLab Corp.	2001										*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 *------------------------------------------------------------------------------*/

#ifndef _O_NAG_G12_H
#define _O_NAG_G12_H

//#importdll "ONAG" // NAG DLL prepared by OriginLab
#pragma dll(ONAG)	
#include <NAG\nag_types.h>


/* begin proto */

/**	g12aac
			computes the Kaplan-Meier,(or product-limit),estimates of
			survival probabilities for a sample of failure times.

Example 1:
	Assume we have "Data1" Worksheet which has 4 columns and 10 rows.
	The first column contain 10 data, and we will put the result,
	the time in the second column, survival probability in the third 
	column, and standard deviation in fourth column.  We also have
	"Data2" Worksheet which has 2 columns and 10 rows, these two columns
	contain 10 data in each column.
	
	int n = 10, nd;
	//Attach six Datasets to these 6 columns
	Dataset dt("data1",0);
	Dataset dtp("data1",1);
	Dataset dp("data1",2);
	Dataset dpsig("data1",3);
	Dataset difreq("data2",0);
	Dataset dic("data2",1);
	
	dt.SetSize(10);
	dtp.SetSize(10);
	dp.SetSize(10);
	dpsig.SetSize(10);
	difreq.SetSize(10);
	dic.SetSize(10);
	
	//Dataset cannot be the parameter of function, but vector can be.
	vector t = dt, tp = dtp, p = dp, psig = dpsig;
	vector<int> ifreq = difreq, ic = dic;
	nag_prod_limit_surviv_fn(n, t, ic, ifreq, &nd, tp, p, psig);
	//Put the result to worksheet "Data1" the second, third and fourth column.
	dtp = tp;
	dp = p;
	dpsig = psig;
	
	 
Example 2:
	The remission times for a set of 21 leukemia patients at 18 distinct 
	time points are read in and the Kaplan-Meier estimate computed and 
	printed. For further details see Gross and Clark (1975), page 242.
	
void test_nag_prod_limit_surviv_fn()
{
	double psig[18];
	double p[18];
	double t[18] = {6.0, 6.0, 7.0, 9.0, 10.0, 10.0, 11.0, 13.0, 16.0, 17.0, 
	19.0, 20.0, 22.0, 23.0, 25.0, 32.0, 34.0, 35.0};
	double tp[18];
	int i, n = 18, nd;
	int ifreq[18] = {1, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1}; 
	int ic[18] = {1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1};
	
	nag_prod_limit_surviv_fn(n, t, ic, ifreq, &nd, tp, p, psig);
	
	printf("nd = %d \n",nd);
	printf("\n   Time     Survival 	 Standard\n");
	printf(" 	      probability    deviation\n\n");
	
	for (i = 0; i < nd; ++i)
		printf(" %6.1f%10.3f %10.3f\n", tp[i], p[i], psig[i]);
}

	The output is following:
	
nd = 7 

   Time     Survival 	 Standard
 	      probability    deviation

    6.0       0.857        0.076
    7.0       0.807        0.087
   10.0       0.753        0.096
   13.0       0.690        0.107
   16.0       0.627        0.114
   22.0       0.538        0.128
   23.0       0.448        0.135


Return:
	This function returns NAG error code, 0 if no error.
	
	11:  On entry, n must not be less than 2: n=_value_.
	604: On entry, ic[_value_]=_value_. The censor code for an observation must be either 0 or 1.
	605: On entry, freq[_value_]=_value_. The value of frequency for an observation must be = 0.
	73:  Memory allocation failed.
	74:  An internal error has occurred in this function. Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.
		
	successfully call of the nag_prod_limit_surviv_fn function.		
*/

	int nag_prod_limit_surviv_fn(
		int n,	//the number of failure and censored times give in t.
		const double t[],	//the failure and censored times;these need not be ordered.
		const int ic[],	//contains the censoring code of the ith observation.
		const int freq[],	//indicates whether frequencies are provided for each failure and censored time point.
		int *nd,	//the number of distinct failure times.
		double tp[],	//contains the ith ordered distinct failure time.
		double p[],	//contains the Kaplan-Meier estimate of the survival probability.
		double psig[]	//contains an estimate of the standard deviation.
	);	// General block design or completely randomized design.
	
	
/**	g12bac
		returns parameter estimates and other statistics that are associated with the Cox proportional
		hazards model for fixed covariates

Example1:
	Assume we have "Data1" Worksheet which is double type and has 4 columns 
	and 42 rows.  The first two columns contain 42 data in each column, and we 
	will put the result, the time in the third column, survival function in 
	the fourth column.  We have "Data2" Worksheet which is integer type and has 
	2 columns and 42 rows.  The first column contains 42 data.
	
	double b[2], cov[3], dev, df;
	double omega[42], res[42], sc[2], se[2]; 
	double tol = 5e-5,  v[336], y[42], offset[42];
	int i, i_1, ip = 1, ip1, iprint = 0, irank, isi[42], sz[1] = {1};
	int j, m = 1, maxit = 20, n = 42, nd, ndmax = 42, ns = 0, tdsur = 1, tdv = 8;
	double c_b42  = 0;
	int c_0 = 0;
	string filename = "c:\\test.txt";
	ip1 = ip +1;
	i_1 = n;
	//Attach five Datasets to these five columns
	Dataset dt("data1",0);	
	Dataset dz("data1",1); 
	Dataset dtp("data1",2);
	Dataset dsur("data1",3);	
	Dataset dic("data2",0);
	
	dt.SetSize(42);
	dz.SetSize(42);
	dtp.SetSize(42);
	dsur.SetSize(42);
	//Dataset cannot be the parameter of function, but vector can be.
	vector t = dt, z = dz, tp = dtp, sur = dsur;
	vector<int> ic = dic;

	for(i =0; i < n; i++)
	{
		y[i] = 1.0 -1.0*ic[i];
		v[i+250] = log(t[i]);
		offset[i] = log(t[i]);
	}	
    nag_glm_poisson(Nag_Log, Nag_MeanInclude, n, z, m, m, sz, ip1, y, NULL, offset,
    				c_b42, &dev, &df, b, &irank, se, cov, v, tdv, tol, maxit, c_0, 
    				filename, c_b42); 
	for(i = 0; i < ip; i++)
		b[i] = b[i+1];
	if(irank != ip +1)
		printf("WARNING: covariate not of full rank\n");		
	nag_surviv_cox_model(n, m, ns, z, m, sz, ip, t, ic, NULL, isi, &dev, b, 
						se, sc, cov, res, &nd, tp, sur, tdsur, ndmax, tol,
						maxit, iprint, filename); 
	//put the result to Worksheet "data1", the third time is time and fourth column
	//is suvival function.
	dtp = tp;
	dsur = sur;
      	
		
Example2:
	The data are the remission times for two groups of leukemia patients(see Gross
	and Clark(1975)p242). A dummy variable indicates which group they come from. An
	initial estimate is computed using the exponential model and then the Cox proportional
	hazard model is fitted and parameter estimates and  the survival function are printed	

void test_nag_surviv_cox_model()
{
	double b[2], cov[3], dev, df;
	double omega[42], res[42], sc[2], se[2]; 
	double sur[42];
	double t[42] = {1, 1, 2, 2, 3, 4, 4, 5, 5, 8, 8, 8, 8, 11, 11, 12, 12, 
	15, 17, 22, 23, 6, 6, 6, 7, 10, 13, 16, 22, 23, 6, 9, 10, 11, 17, 19,
	20, 25, 32, 32, 34, 35};
	double z[42] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
	0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
	double tol, tp[42], v[336], y[42], offset[42];
	int i, i_1, ip, ip1, iprint = 0, irank, isi[42], sz[1];
	int ic[42] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
	int j, m = 1, maxit = 20, n = 42, nd, ndmax = 42, ns = 0, tdsur = 1, tdv;
	double c_b42  = 0;
	int c_0 = 0;
	
	//Get the origin folder path.
	char sztemp[100];
	LT_get_str("%Y", sztemp, 100);
	
	//Define a filename and working path.
	string str1, filename;
	str1 = sztemp;
	filename = str1 + "text.txt";
	printf("%s\n", filename);
	
	sz[0] =1;
	ip =1;
	ip1 = ip +1;
	tol = 5e-5;
	i_1 = n;
	tdv = 8;
	for(i =0; i < n; i++)
	{
		y[i] = 1.0 -1.0*ic[i];
		v[i+250] = log(t[i]);
		offset[i] = log(t[i]);
	}
	
    nag_glm_poisson(Nag_Log, Nag_MeanInclude, n, z, m, m, sz, ip1, y, NULL, offset,
    				c_b42, &dev, &df, b, &irank, se, cov, v, tdv, tol, maxit, c_0, 
    				filename, c_b42);
    
	for(i = 0; i < ip; i++)
		b[i] = b[i+1];
	if(irank != ip +1)
		printf("WARNING: covariate not of full rank\n");
		
	nag_surviv_cox_model(n, m, ns, z, m, sz, ip, t, ic, NULL, isi, &dev, b, 
						se, sc, cov, res, &nd, tp, sur, tdsur, ndmax, tol,
						maxit, iprint, filename); 
	printf("\n");
	printf("Parameter          Estimate");
	printf("       Standard Error");
	printf("\n");
	for(i = 0; i < ip; i++)
		printf("%6ld          %8.4f         %8.4f\n",i + 1, b[i], se[i]);
	printf("\n");
	printf("Deviance = %13.4f\n",dev);
	printf("\n");
	printf("    Time     Survivor Function");
	printf("\n");
	ns = 1;
	for(i = 0; i < nd; i++)
	{
		printf("%10.0f",tp[i]);
		for(j = 0; j < ns; j++)
		{
			printf("      %8.4f",sur[i * tdsur +j]);
			if( (j + 1) % 3 == 0)
				printf("\n");
		}
		printf("\n");
	}      	
}

	The output is following:
	
C:\sftest\Gr394\text.txt

Parameter          Estimate       Standard Error
     1           -1.5091           0.4096

Deviance =      172.7592

    Time     Survivor Function
         1        0.9640
         2        0.9264
         3        0.9065
         4        0.8661
         5        0.8235
         6        0.7566
         7        0.7343
         8        0.6506
        10        0.6241
        11        0.5724
        12        0.5135
        13        0.4784
        15        0.4447
        16        0.4078
        17        0.3727
        22        0.2859
        23        0.1908
	
Return:
	This function returns NAG error code, 0 if no error.
	
	11:  On entry, m must not be less than 1: m = <value>.  On entry, n must not be less than 2: n = <value>.  On entry, ns must not be less than 0: ns = <value>.  On entry, tdsur must not be less than 1: tdsur = <value>.  On entry, max_iter must not be less than 0: max_iter = <value>.
	17:  On entry, tdsur = <value> while ns = <value>.  These parameters must satisfy tdsur >= ns.  On entry, tdz = <value> while m = <value>. These parameters must satisfy tdz >= m.
	117: On entry, tol = <value>, machine precision (X02AJC) = <value>.  Constraint: tol >= 10 * machine precision.
	78:  Cannot open file outfile for appending.
	79:  Cannot close file outfile.
	103: On entry, sz[<value>] =<value>.  Constrant: sz[<value>] >= 0.  On entry, ic[<value>] = <value>.  Constraint: ic[<value>] = 0 or 1.  On entry, isi[<value>] = <value>.  Constraint: 0 <= isi[<value>] <= ns.
	607: On entry, ip = <value> and the number of non zero values of sz = <value>.  Constraint: ip = the number of non-zero values of sz.
	608: On entry, the number of values of sz[i] > 0 is <value>, n = <value> and excluded observations with isi[i] = 0 is <value>.  Constraint: the number of values of non-zero sz must be less than n - excluded observations.
	609: On entry, ndmax is = <value> while the output value of nd = <value>.  Constraint: ndmax >= nd.
	610: The matrix of second partial derivatives is singular.  Try different starting values or include fewer covariates.
	611: Overflow has been detected.  Try different starting values.
	612: Convergence has not been achieved in max_iter iterations.  The progress towards convergence can be examined by using by setting iprint to >= 1.  Any non _convergence decreasing the deviance maybe due to a linear combination of covariates being monotonic with time.  Full results are returned.
	613: In the current iteration 10 step halvings have been performed withour decreasing the deviance from the previous iteration.  Convergence is assumed.
	73:  Memory allocation failed.
	74:  An internal error has occurred in this function. Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.
		
	successfully call of the nag_surviv_cox_model function.
*/

int nag_surviv_cox_model(
	int n,	//the number of data points
	int m,	//the number of covariates
	int ns,	//the number of strata
	const double z[],	//the covariates
	int tdz,	//the second dimension of array z;
	const int sz[],	//indicates which subset of covariates is to be included
	int ip,	//the number of covariates included as indicated by sz
	const double t[],	//the vector of n failure censoring times
	const int ic[],	//the status of the individual at time t given in t
	const double omega[],	//if an offset is set, it contain value \omega[i], otherwise it is NULL;
	const int isi[],	//the stratum indicator
	double *dev,	//the deviance
	double b[],	//initial estimates of the covariate coefficient parameters
	double se[],	//the asymptotic standard error of the estimate
	double sc[],	//the value of the score function
	double cov[],	//the variance-covariance matirx of the parameter estimates
	double res[],	//the residuals
	int *nd,	//the number of distinct failure times
	double tp[], //the distinct failure time
	double sur[],	//the estimated survival function
	int tdsur,	//the second dimension of the array
	int ndmax,	//the first dimension of the array
	double tol,	//the accuracy required for the estimation
	int max_iter,	//the maximun number of iterations
	int iprint,	//indicates if the printing information on the iteration is required
	const char *outfile	//the name of the file into which information is to be output.
);









#endif //!_O_NAG_G12_H