//  Created by Joe Becker on 3/7/05.
//  Copyright 2005 __MyCompanyName__. All rights reserved.

#import "GmtShipTrack.h"

@implementation GmtShipTrack

- (void)dealloc
{
	[gmtShipTrackRecordList release];
	[min release];
	[max release];
	[super dealloc];	
}

- (id)init
{
	if (self = [super init]) {
		gmtShipTrackRecordList = [[NSMutableArray alloc] init];
		carterCount = 0;
	}
	return self;
}

-(id)resetStats {
	
	numGood			= [[GmtShipTrackRecord alloc] initTo: 0];
	sumOfX			= [[GmtShipTrackRecord alloc] initTo: NAN];
	sumOfXSqr		= [[GmtShipTrackRecord alloc] initTo: NAN];
	sumOfY			= [[GmtShipTrackRecord alloc] initTo: NAN];
	sumOfYSqr		= [[GmtShipTrackRecord alloc] initTo: NAN];
	sumOfXY			= [[GmtShipTrackRecord alloc] initTo: NAN];
	sumOfDeltaY		= [[GmtShipTrackRecord alloc] initTo: NAN];
	sumOfXDeltaY	= [[GmtShipTrackRecord alloc] initTo: NAN];
	sumOfDeltaYSqr	= [[GmtShipTrackRecord alloc] initTo: NAN];
	
	mean			= [[GmtShipTrackRecord alloc] initTo: NAN];
	stdDev			= [[GmtShipTrackRecord alloc] initTo: NAN];
	slope			= [[GmtShipTrackRecord alloc] initTo: NAN];
	offset			= [[GmtShipTrackRecord alloc] initTo: NAN];
	
	[mean   retain];
	[stdDev retain];
	[slope  retain];
	[offset retain];
	
	[numGood	retain];
	[sumOfX		retain];
	[sumOfXSqr	retain];
	[sumOfY		retain];
	[sumOfYSqr	retain];
	[sumOfXY	retain];
	[sumOfDeltaY retain];
	[sumOfXDeltaY retain];
	[sumOfDeltaYSqr retain];
	
	return self;
}

-(id)readFromPath:(NSString *)filePath
{	
	[self resetStats];
	
	// unlike other stats: we want the max, min of all non-NaN data so graph to full extent "looks right"
	max    = [[GmtShipTrackRecord alloc] initTo: NAN];
	min    = [[GmtShipTrackRecord alloc] initTo: NAN];
	
	NSString *trackData = [NSString stringWithContentsOfFile:filePath];
	NSScanner *scanner  = [NSScanner scannerWithString: trackData];
	
	[scanner scanUpToString:@"\n" intoString:&firstLineOfFile];
	[scanner scanUpToString:@"\n" intoString:&secondLineOfFile];
	
	[firstLineOfFile retain];
	[secondLineOfFile retain];
	
	while ([scanner isAtEnd] == NO) 
	{
		GmtShipTrackRecord *rec;
		NSString *stringVal;
		
		if (![scanner scanUpToString:@"\n" intoString:&stringVal]) {
			NSRunAlertPanel(@"GmtShipTrack.m",
							@"String could not be read",nil,nil,nil);
			return nil;
		}
		
		rec = [[GmtShipTrackRecord alloc] initFromString:stringVal];
		[rec setTag:GMT_TAG];
		[gmtShipTrackRecordList addObject:rec];
		
		max = [self maxRecordOf:max:rec];
		min = [self minRecordOf:min:rec];
	}
	
	// Calculate mean, standard deviation, slope and offset
	
	[self calculateSumsForStats];
	
	// this retain is crucial since we are adding the stat record to another object, they would be freed.
	
	[max    retain];
	[min    retain];
	return self;
}

- (id)writeToPath:(NSString *)filePath
{	
	id obj;
	
	NSMutableString	*myOut = [NSMutableString stringWithCapacity:1000];
	
	[myOut appendString:firstLineOfFile];
	[myOut appendString:@"\n"];
	[myOut appendString:secondLineOfFile];	
	[myOut appendString:@"\n"];
	
	NSEnumerator *en = [gmtShipTrackRecordList objectEnumerator];
	while (obj = [en nextObject]) {
		[myOut appendString:[obj stringRep]];
	}
	
	if (![myOut writeToFile:filePath atomically:YES]) {
		NSRunAlertPanel(@"GmtShipTrack.m",
						@"file and/or string could not be written",nil,nil,nil);
		return nil;
	}
	
	return self;	
}

- (id)gmtShipTrackRecordList {
	return gmtShipTrackRecordList;
}

- (BOOL) isNaN:(double)x 
{	
	if (x <= 0 || x >= 0)
	{
		// x is not a NaN
		return NO;
	}
	else
	{
		return YES;
	}
}

- (id)setDoubleWithNaN:(id)cell:(double)data
{
	if ([self isNaN:data]){
		[cell setStringValue:@"NaN"];
	} else {
		[cell setDoubleValue:data];
	}
	return cell;
}

- (double) addAllowingForNAN:(double)a:(double)b
{
	if ([self isNaN:a] && [self isNaN:b])
	{
		// both inputs are NAN
		return NAN;			
	}
	else if ([self isNaN:a])
	{
		// a is NAN
		return b;			
	}
	else if ([self isNaN:b])
	{
		// b is NAN
		return a;			
	}
	else 			
		// Neither a or b is NaN
	{
		return a+b;
	}
}


- (double) maxOrMinAllowingForNAN:(BOOL)wantMax:(double)a:(double)b
{
	if ([self isNaN:a] && [self isNaN:b])
	{
		// both inputs are NAN
		return NAN;			
	}
	else if ([self isNaN:a])
	{
		// a is NAN
		return b;			
	}
	else if ([self isNaN:b])
	{
		// b is NAN
		return a;			
	}
	else 			
		// Neither a or b is NaN
	{
		if (wantMax)
		{
			if(a > b)
				return a;
			else
				return b;
		}
		else
		{
			if(a < b)
				return a;
			else
				return b;
		}
	}
}

- (GmtShipTrackRecord *)maxMinRecordOf:(BOOL)wantMax:(GmtShipTrackRecord *)a:(GmtShipTrackRecord *)b
{
	GmtShipTrackRecord *tmp = [[GmtShipTrackRecord alloc] initTo:NAN];
	int i;
	
	for (i=TIME_FIELD; i<=CC_FIELD; i++)
	{
		[tmp setData:i:[self maxOrMinAllowingForNAN:wantMax:[a data:i]:[b data:i]]];
	}
	
	return tmp;
}

- (GmtShipTrackRecord *)maxRecordOf:(GmtShipTrackRecord *)a:(GmtShipTrackRecord *)b
{
	return [self maxMinRecordOf:YES:a:b];
}

- (GmtShipTrackRecord *)minRecordOf:(GmtShipTrackRecord *)a:(GmtShipTrackRecord *)b
{
	return [self maxMinRecordOf:NO:a:b];
}

- (void)calculateSumsForStats
{
	GmtShipTrackRecord	*obj= [gmtShipTrackRecordList objectAtIndex:0];
	NSEnumerator		*en = [gmtShipTrackRecordList objectEnumerator];
	int					i;
		
	[self resetStats];
	
	while (obj = [en nextObject])
	{
		for (i=TIME_FIELD; i<=CC_FIELD; i++)
		{
			double actual	= [obj data:i];
			double predicted= [obj predictedData:i];
			
			//FIXME: this is a hack Dave wanted so that we don't use blatently bad predicted when calculating stats
			switch (i) {
				case TC_FIELD : {
					predicted = [obj pt];
					if (predicted > -4.0) {
						predicted = NAN;						
					}
					break;
				}				
				default:
					break;
			}
			
			double x = predicted;
			double y = actual;
			double deltaY = actual-predicted;

			int maskedFlag = ([obj flag] & ~1);
			int mask = (1<<i);
			
			maskedFlag &= mask;
			
			if (maskedFlag == 0 && ![self isNaN:x] && ![self isNaN:y]) { // do not include flagged data
				[numGood		setData:i:[self addAllowingForNAN:[numGood			data:i]	: 1]];
				[sumOfX			setData:i:[self addAllowingForNAN:[sumOfX			data:i]	: x]];
				[sumOfY			setData:i:[self addAllowingForNAN:[sumOfY			data:i]	: y]];
				[sumOfXY		setData:i:[self addAllowingForNAN:[sumOfXY			data:i]	: x*y]];
				[sumOfXSqr		setData:i:[self addAllowingForNAN:[sumOfXSqr		data:i] : x*x]];
				[sumOfYSqr		setData:i:[self addAllowingForNAN:[sumOfYSqr		data:i] : y*y]];

				[sumOfDeltaY	setData:i:[self addAllowingForNAN:[sumOfDeltaY		data:i]	: deltaY]];
				[sumOfXDeltaY	setData:i:[self addAllowingForNAN:[sumOfXDeltaY		data:i]	: x*deltaY]];
				[sumOfDeltaYSqr	setData:i:[self addAllowingForNAN:[sumOfDeltaYSqr	data:i] : deltaY*deltaY]];
			}
		}
	}
	
	for (i=TIME_FIELD; i<=CC_FIELD; i++)
	{
		// code looks like hell with all the [sumOfX data:i] stuff...
		double cnt		= [numGood data:i];
		double X		= [sumOfX data:i];
		double Y		= [sumOfY data:i];
		double DeltaY	= [sumOfDeltaY data:i];
		double XY		= [sumOfXY data:i];
		double XDeltaY	= [sumOfXDeltaY data:i];
		double XSqr		= [sumOfXSqr data:i];
		double YSqr		= [sumOfYSqr data:i];
		double DeltaYSqr= [sumOfDeltaYSqr data:i];
		
		if (cnt == 0)
		{
			[mean	setData:i:	NAN];
			[stdDev	setData:i:	NAN];
			[slope	setData:i:	NAN];
			[offset	setData:i:	NAN];
		} else {
			[mean	setData:i:	Y/cnt];
			[stdDev	setData:i:	sqrt(YSqr/cnt - Y*Y/cnt/cnt)];
			[slope	setData:i:	(cnt*XY - X*Y )/(cnt*XSqr - X*X)];
			[offset	setData:i:	(Y*XSqr - X*XY)/(cnt*XSqr - X*X)];

			[mean	setData:i:	DeltaY/cnt];	// per Dave Sandwell
			[stdDev	setData:i:	sqrt(DeltaYSqr/cnt - DeltaY*DeltaY/cnt/cnt)];	// per Dave
			[slope	setData:i:	Y/X];			// per Walter Smith
			[offset	setData:i:	(DeltaY*XSqr - X*XDeltaY)/(cnt*XSqr - X*X)];	// common sense?
			[offset	setData:i:	0];				// per Dave
		}
	}
}

- (GmtShipTrackRecord *)max {return max;}
- (GmtShipTrackRecord *)min {return min;}
- (GmtShipTrackRecord *)mean {return mean;}
- (GmtShipTrackRecord *)stdDev {return stdDev;}
- (GmtShipTrackRecord *)slope {return slope;}
- (GmtShipTrackRecord *)offset {return offset;}
- (GmtShipTrackRecord *)numGood {return numGood;}

- (GmtShipTrackRecord *)objectAt:(int)dataField:(float)x:(float)y
{
	id obj;
	GmtShipTrackRecord *rec;
	NSEnumerator *en = [gmtShipTrackRecordList objectEnumerator];
	
	// Only time needs to have an offset from first record added
	if (dataField == TIME_FIELD) {
		x +=[[gmtShipTrackRecordList objectAtIndex:0] time];
	}
	
	while (obj = [en nextObject]) {
		rec = obj;
		BOOL hit;
		const float epsilion = 0.001;
		switch (dataField) {
			case TIME_FIELD:hit =									fabs([rec time]-x) <  epsilion; break;
			case CC_FIELD : hit =									fabs([rec cc] - y) <  epsilion; break;
			case TT_FIELD : hit =									fabs([rec tt] - y) <  epsilion; break;
			case LAT_FIELD: hit = fabs([rec lon]- x) <  epsilion && fabs([rec lat]- y) <  epsilion; break;
			case TC_FIELD : hit = fabs([rec pt] - x) <  epsilion && fabs([rec tc] - y) <  epsilion; break;
			case GF_FIELD : hit = fabs([rec gs] - x) <  epsilion && fabs([rec gf] - y) <  epsilion; break;
			case GT_FIELD : hit = fabs([rec gs] - x) <  epsilion && fabs([rec gt] - y) <  epsilion; break;
			case MA_FIELD : hit = fabs([rec pm] - x) <  epsilion && fabs([rec ma] - y) <  epsilion; break;
			default:
				hit = NO;
				NSRunAlertPanel(nil,@"GmtShipTrack.m objectAt: User selected a unexpected data tag.",nil,nil,nil);
		}
		if (hit)
			return rec;
	}
	return nil;
}

- (void)carterCorrectionAddNotSub:(BOOL)addNotSub
{
	GmtShipTrackRecord *obj;
	NSEnumerator *en = [gmtShipTrackRecordList objectEnumerator];

	while (obj = [en nextObject]) {
	
		if ([self isNaN:[obj data:TC_FIELD]])
			continue;
		
		if (addNotSub) {
			[obj setData:TC_FIELD :
				[self addAllowingForNAN:[obj tc] :  [obj cc]]];
		} else {
			[obj setData:TC_FIELD :
				[self addAllowingForNAN:[obj tc] : -[obj cc]]];
		}
	}
	[self calculateSumsForStats];
}

- (void)addCarterCorrection {
	carterCount++;
	[self carterCorrectionAddNotSub:YES];
}

- (void)subCarterCorrection {
	carterCount--;
	[self carterCorrectionAddNotSub:NO];
}

- (double)carterCount {
	return carterCount;
}

- (double)setCarterCount:(double)cnt {
	return carterCount = cnt;
}

- (NSRect)extent:(int)field
{
	float minTime = [[self min] time];
	float maxTime = [[self max] time];
	
	float minData = [[self min] data:field];
	float maxData = [[self max] data:field];
	
	return 	NSMakeRect(0, minData, maxTime-minTime, maxData-minData);
}

@end
