survivalanalysis.c

图像处理的压缩算法

/*------------------------------------------------------------------------------*
 * File Name:	SurvivalAnalysis.c: Contains OriginC source code for all		*
 * 	Survival Analysis features.													*
 * Creation: GJL 1/28/2002														*
 * Purpose: Origin C file														*
 * Copyright ( c ) OriginLab Corp.	2000-2002									*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 *	GJL 1/28/2002: Split off source code from OStat.c							*
 *	CP/GJL 7/21/02 v7.0348 OSTAT_COMPILE_PROBLEM_WITH_TYPE_ERROR_MSG			*
 *		I have added Type_ErrorMsg1 function to replace all those that need		*
 *		second string arguments													*
 *	GJL 11/14/02 v7.0434 USE_OC_WKS_TYPE_NOT_LT									*
 *------------------------------------------------------------------------------*/

////////////////////////////////////////////////////////////////////////////////////
//
// System includes
#include <Origin.h>
//#include <OC_const.h>
//#include <string.h>
//#include <data.h>
//#include <math.h>
//#include <utilities.h>
//#include <stdio.h>
//#include <page.h>

// Includes definitions of SurvivalAnalysis constants, structures, and non-localized strings.
//  Also includes NAG definitions and prototypes of NAG functions used by SurvivalAnalysis.
#include "SurvivalAnalysis.h"

// Includes prototypes of Application Utilities.
#include "App_Utils.h"

// Includes definitions of all SurvivalAnalysis.c localized strings. $ causes CodeBuilder to look for 
// correct version of Local.h depending on Windows Regional settings
#include "$Local.h"
//
////////////////////////////////////////////////////////////////////////////////////

/**
		Implements the Kaplan-Meier (Product-Limit) Estimator. The computational
		engine is the NAG function nag_prod_limit_surviv_fn (g12aac). Additional
		reference is Chapter 2 of Applied Survival Analysis: Regression Modeling
		of Time to Event Data (1998), Hosmer, D.W. and Lemeshow, S. (1998) John
		Wiley and Sons Inc., New York, NY.
	Example:
		See [KaplanMeier] section of KaplanMeier.OGS for sample call.
	Parameters:
		nPts=Number of data points in Time and Censor datasets
		nResults=Bitwise integer specifying output options
			1=Output Survivorship Function in Resuls Log
			2=Send output to temporary worksheet in addition to Results Log
			4=Create plot of Survivorship Function (implemented in KaplanMeier.OGS)
			8=Plot Errors or Confidence Limits along with Survivorship Function  
		strSurvivalWKS=Name of temporary worksheet holding sorted copies of all
			required input (Time and Censor) and output datasets   
	Return:
		Depending on nResults, outputs Summary of Event and Censored Times, the 
		Survivorship Function, and Quartile Estimates to the Results Log and to
		a worksheet. Returns SA_NO_ERROR on successful exit or an SA_ error code.
*/
int saKaplanMeier( int nPts, int nResults, double dConLevel, string strSurvivalWKS )
{
	int iSize, iErr, nd;
	string strDatasetName;
	BOOL bErr;

	// Read Time dataset (including failure and censored times) into vector to pass to NAG
	strDatasetName = strSurvivalWKS + SA_TIME_DATASET_NAME;	
	Dataset dsTIME( strDatasetName );
	iSize = dsTIME.GetSize();
	vector<double> vTIME;
	vTIME.SetSize( iSize );
	vTIME = dsTIME;	
	dsTIME.Detach(); // No longer need input dataset

	// Read Censor dataset into vector to pass to NAG	
	strDatasetName = strSurvivalWKS + SA_CENSOR_DATASET_NAME;
	Dataset dsCENSOR( strDatasetName );
	iSize = dsCENSOR.GetSize();
	vector<int> vCENSOR;
	vCENSOR.SetSize( iSize );
	vCENSOR = dsCENSOR;
	
	// Perform basic stats on Censor dataset (all 0's and 1's) for summary output
	Dataset *pdsCENSOR = &dsCENSOR;
	BasicStats bsStat;
	BasicStats *pbsStat = &bsStat;
	bErr = Data_sum( pdsCENSOR, pbsStat );
	if( bErr == FALSE )
		return SA_DATA_SUM_ERROR1;
	
	dsCENSOR.Detach(); // No longer need input dataset

	// Contains distinct ordered failure times on output from NAG (does not include duplicates or censored times)
	vector<double> vTP;
	vTP.SetSize( iSize );

	// Contains Kaplan-Meier estimate of Survival Probability for each time in vTP on output from NAG 
	vector<double> vP;
	vP.SetSize( iSize );

	// Contains Kaplan-Meier estimate of error (standard deviation) of each probability in vP on output from NAG
	vector<double> vPSIG;
	vPSIG.SetSize( iSize );

	// From <NAG\OCN_g12.h>: g12aac nag_prod_limit_surviv_fn: Compute Kaplan-Meier estimate of
	// Survivorship Function. Frequency is NULL (count each point once; nd is the number of distinct
	// failure times output in vTP.
	iErr = nag_prod_limit_surviv_fn( nPts, vTIME, vCENSOR, NULL, &nd, vTP, vP, vPSIG );
	if( iErr != NE_NOERROR )
		return SA_PROD_LIMIT_ERROR;

	// Compute Quartile Estimates (Percentiles)
	double ardPercentile[3];
	iErr = saKaplanMeier_Compute_Percentiles( nd, vTP, vP, ardPercentile );
	if( iErr != SA_NO_ERROR )
		return iErr;
	
	// Compute Errors or Confidence Limits (upper and lower) for Percentiles
	double ardPerLLim[3];
	double ardPerULim[3];
	iErr = saKaplanMeier_Compute_Con_Limits( dConLevel, nd, vTP, vP, vPSIG, ardPerLLim, ardPerULim );
	if( iErr != SA_NO_ERROR )
		return iErr;

	// Output results in Results Log
	iErr = saKaplanMeier_Output_to_ResultsLog( nResults, nPts, vTIME, vCENSOR, pbsStat->total, vP, vPSIG, dConLevel,
		 ardPercentile, ardPerLLim, ardPerULim );
	if( iErr != SA_NO_ERROR )
		return iErr;

	// Output results in worksheet? If user wants worksheet or plot write everything out to worksheet
	if( nResults > 1 )
	{
		// Function to output results in worksheet
		iErr = saKaplanMeier_Output_to_Worksheet( strSurvivalWKS, nPts, vTIME, vCENSOR, nd , vP, vPSIG, dConLevel,
			 ardPercentile, ardPerLLim, ardPerULim );
		if( iErr != SA_NO_ERROR )
			return iErr;
	}

	return SA_NO_ERROR;
}

/**
		Computes Quartile Estimates (the 25th, 50th, and 75th Percentiles) for the Kaplan-Meier (Product-Limit)
		Estimator function.
	Example:
		See function saKaplanMeier above for sample call.
	Parameters:
		nd=Number of distinct failure times in vTP (nd==size of vTP)
		vTP=Vector containing ordered failure times output by NAG (does not include duplicate or censored times)
		vP=Vector containing Kaplan-Meier estimate of Survival Probability for each time in vTP output by NAG
		ardPercentile=Three element array containing the 25th, 50th and 75th Percentiles
	Return:
		Array of doubles holding the 25th, 50th, and 75th percentiles of the Kaplan-Meier estimate
		of the Survivorship Function. Returns SA_NO_ERROR on successful exit.
*/
int saKaplanMeier_Compute_Percentiles( int nd, vector<double> &vTP, vector<double> &vP, double *ardPercentile )
{
	int ii, jj, iPercInterp, ivPofii, iP;
	double dPercInterp;

	// Get LabTalk system.graph.percinterp setting (from graph tab in Tools:Options DB).  ==1 when
	//  percentiles falling between two elements should be averaged. 
	LT_evaluate( "system.graph.percinterp", &dPercInterp );
	iPercInterp = dPercInterp;

	// Initialize values to negative time as flag for missing value
	ardPercentile[0] = -1;
	ardPercentile[1] = -1;
	ardPercentile[2] = -1;	
	jj = 0;
	iP = 75000; // Integerized percentile threshold for 25th percentile
	
	// Loop on three percentiles (25th, 50th, and 75th) as long as all non-missing value
	//  surviviorship function values have not been exausted
	for( ii = 0; ( ii < nd ) && ( jj < 3 ); ii++ )
	{
		// Allow for equality within tolerance by rounding to 10,000'sth place and integerizing
		ivPofii = vP[ii] * 100000 + 0.5;

		// If survivorship function value is less than or equal to percentile threshold...
		if( ivPofii <= iP )
		{
			// If equal to percentile threshold and interpolate option is enabled then try to interpolate...
			if( ivPofii == iP && iPercInterp )
			{
				// If next element exists to interpolate with then interpolate...
				if( ii + 1 < nd )  
					ardPercentile[jj] = ( vTP[ii] + vTP[ii + 1] ) / 2; // Interpolate
				// Else leave as missing value by not setting
			}
			else
				ardPercentile[jj] = vTP[ii]; // Don't interpolate
			ii--; // Decrement loop index to start looking for next percentile with current time value
			jj++; // Got current percentile...now get next one
			iP = 75000 - jj * 25000; // Compute integerized threshold for next percentile
		}
	}

	return SA_NO_ERROR;
}

/**
		Computes Upper and Lower Confidence Limits for the Percentiles of the Kaplan-Meier (Product-Limit)
		Estimator function.
	Example:
		See function saKaplanMeier above for sample call.
	Parameters:
		dConLevel=Confidence level
		nd=Number of distinct failure times in vTP (nd==size of vTP)
		vTP=Vector containing distinct ordered failure times output by NAG (does not include duplicate or 
			censored times)
		vP=Vector containing Kaplan-Meier estimate of Survival Probability for each time in vTP output by NAG
		vPSIG=Vector containing standard deviation of each probability in vP output by NAG	
		ardPerLLim=Three element array containing a lower confidence limit for the 25th, 50th, and 75th
			percentiles
		ardPerULim=Three element array containing an upper confidence limit for the 25th, 50th, and 75th
			percentiles
	Return:
		Two arrays of doubles holding the upper and lower confidence limits for the 25th, 50th, and 75th
		percentiles of the Kaplan-Meier estimate of the Survivorship Function. Returns SA_NO_ERROR on
		successful exit.
*/
int saKaplanMeier_Compute_Con_Limits( double dConLevel, int nd, vector<double> &vTP, vector<double> &vP,
	 vector<double> &vPSIG, double *ardPerLLim, double *ardPerULim )
{
	int ii, jj;
	double dP, dChi;
	NagError neErr;
	NagError *pneFail = &neErr;		// NAG error structure

	// Initialize all to -1 for later output as NANUM if not successfully computed
	ardPerLLim[0] = -1;
	ardPerLLim[1] = -1;
	ardPerLLim[2] = -1;	
	ardPerULim[0] = -1;
	ardPerULim[1] = -1;
	ardPerULim[2] = -1;

	// From <NAG\OCN_g01.h>: g01fcc nag_deviates_chi_sq: Compute Chi-Square deviate 
	dChi = nag_deviates_chi_sq( dConLevel, 1, pneFail );
	if( pneFail->code != NE_NOERROR )
		return SA_CHI_SQ_ERROR1;
	
	// Loop on three percentiles 25th, 50th, and 75th
	for( jj = 0; jj < 3; jj++ )
	{
		// Compute non-integerized threshold for current percentile
		dP = 0.75 - jj * 0.25;
		
		// Loop through all distict failure times
		for( ii = 0; ii < nd; ii++ )
		{
			// If condition is TRUE then time value may be upper or lower limit 
			if( fabs( vP[ii] - dP ) <= sqrt( dChi * pow( vPSIG[ii],2 ) ) )
			{
				if( ardPerLLim[jj] < 0 ) // Lower confidence limit is the first time the condition is TRUE
					ardPerLLim[jj] = vTP[ii];
				if( ii + 1 < nd ) // Upper confidence limit is the next time value the last time the condition is TRUE 
					ardPerULim[jj] = vTP[ii + 1];				
			} // End if condition TRUE
		} // End loop through all failure times
	} // End loop on three percentiles
	
	return SA_NO_ERROR;
}

/**
		Outputs the results of the Kaplan-Meier (Product-Limit) Estimator to the Results Log.
	Example:
		See the function saKaplanMeier above for sample call.
	Parameters:
		nResults=Bitwise flag (bit 1) indicating whether to output the Survivorship Function
			to the Results Log
		nPts=Number of data points in Time and Censor datasets
		vTIME=Vector of input Time values
		vCENSOR=Vector of input Censor values
		dCensoredPts=Number of censored data points
		vP=Vector containing Survivorship Function values
		vPSIG=Vector containing Survivorship Function Error (standard deviation) values
		dConLevel=Confidence Level
		ardPercentile=Array containing 25th, 50th, and 75th Percentile values
		ardPerLLim=Array containing Lower Limits for Percentile values
		ardPerULim=Array containing Upper Limits for Percentile values
	Return:
		Outputs the results of the Kaplan-Meier Estimator (including a Main Header and Summary
		table, a Survivorship Function with Errors table, and a Quartile Estimates table) to the
		Results Log and returns SA_NO_ERROR on successful exit or an SA_ error code.
*/
int saKaplanMeier_Output_to_ResultsLog( int nResults, int nPts, vector<double> &vTIME, vector<int> &vCENSOR, double dCensoredPts,
	 vector<double> &vP, vector<double> &vPSIG, double dConLevel, double *ardPercentile, double *ardPerLLim,
	 double *ardPerULim )
{
	int ii, jj, iErr;
	string strDB;
	string strTimeData;
	string strCensorData;
	string strCensorValue;
	char szTemp[40];
	double dPrevUncensoredTime;

	string strOut;	
	int iRow, iCol;