// calibrate.cpp : Defines the entry point for the console application.
//
// Developed by Alex E Pozhitkov, 2012, University of Washington

//#include "stdafx.h"
#include <stdio.h>
#include <iostream>
#include <fstream>
#include <map>
#include <valarray>
#include <vector>
#include <string>
#include <algorithm>
#define _TABOUTPUT string("-T")
#define _FREUNDLICH string("-freund")
#define _LINEAR string("-linear")
#define _LANGMUIR string("-langmuir")
#define _R2FILTER string("-R2")

using namespace std;

enum CURVEID
{
	FREUNDLICH_ID,
	LINEAR_ID,
	LANGMUIR_ID
};

int SearchParameter(int argc, char* argv[], const string& sSwitch)
{
	int iParamIndex = 1;
	for (iParamIndex = 1; iParamIndex < argc; iParamIndex++)
	{
		if(sSwitch == argv[iParamIndex]) 
			break;
	}
	return (iParamIndex == argc) ? -1 : iParamIndex;
}

struct POINT
{
	double m_dX;
	double m_dY;
	POINT(double x, double y) : m_dX(x), m_dY(y) {}
};

inline double getX(const POINT& pt)
{
	return pt.m_dX;
}

inline double getY(const POINT& pt)
{
	return pt.m_dY;
}

inline double power2(double d)
{
	return d * d;
}

struct LINFITPARAMS
{
	double a;
	double b;
	double R2;
};

struct PROBECALIBRATION
{
	CURVEID CurveTypeID;
	LINFITPARAMS lfParameters;
	PROBECALIBRATION(CURVEID type_id, LINFITPARAMS lfp) : CurveTypeID(type_id), lfParameters(lfp) {}
	bool operator< (const PROBECALIBRATION &s) const {return lfParameters.R2 < s.lfParameters.R2;}
};
 

bool LinearFit(const valarray<double>& vx, const valarray<double>& vy, LINFITPARAMS* pLinFit)
{
	int S = vx.size(); //number of points
	if(S != vy.size())
		return false;
	double Sx = vx.sum();
	double Sxx = (vx * vx).sum();
	double Sy = vy.sum();
	double Sxy = (vx * vy).sum();
	double X = Sx / S;  // mean x
	double Y = Sy / S;  // mean y
	double sdX = sqrt(((vx - X).apply(power2)).sum() / (S - 1)) ;
	double sdY = sqrt(((vy - Y).apply(power2)).sum() / (S - 1)) ;
	pLinFit->a = (Sxx * Sy - Sx * Sxy) / (S * Sxx - Sx * Sx);
	pLinFit->b = (S * Sxy - Sx * Sy) / (S * Sxx - Sx * Sx);
	pLinFit->R2 = pow((((Sxy - Y * Sx) - X * Sy) + S * X * Y) / ((sdX * sdY) * (S - 1)), 2);
	return true;
}

int main(int argc, char* argv[])
{
try
{
#if _DEBUG
	__asm{int 3}
#endif
	string sLine;
	double dAvSi;
	double dConc;
	double dErr;
	float fR2cutOff = 0;
	string sProbeID;
	typedef vector<POINT> POINTS;
	typedef map<string, POINTS> PROBEMEASUREMENTS;
	typedef map<string, vector<PROBECALIBRATION> > CALIBRATEDPROBES;

	PROBEMEASUREMENTS Calibrations;

	if(argc == 1)
	{
		cerr << "calibrate  < calibration_series.txt [-T] [-R2 {0.9}] [-freund {F}] [-linear {LN}] [-langmuir {LG}]  > out.txt\n";
		cerr << "\nInput format:  concentration probeID signalIntens stderror\n";
		cerr << "-T\ttabulated output for inspection elsewhere\n";
		cerr << "-R2 {0.9}\tR2 cut off, e.g 0.9, provided by user\n";
		cerr << "-freund {F}\tFreundlich calibration; output marker, e.g. F, provided by user\n";
		cerr << "-linear {LN}\tLinear calibration; output marker, e.g. LN, provided by user\n";
		cerr << "-langmuir {LG}\tLangmuir calibration; output marker, e.g. LG, provided by user\n";
		cerr << "When more than 1 calibration curve suggested, the one with highest R2 is saved in the output.";
		return -1;
	}

	if(SearchParameter(argc, argv, _R2FILTER) != -1)
		sscanf(argv[SearchParameter(argc, argv, _R2FILTER) + 1], "%f", &fR2cutOff);

	while (!cin.eof())
	{
		cin >> dConc >> sProbeID >> dAvSi >> dErr;
		Calibrations[sProbeID].push_back(POINT(dConc, dAvSi));
	}

	if(SearchParameter(argc, argv, _TABOUTPUT) != -1)
	{
		for(PROBEMEASUREMENTS::const_iterator CI = Calibrations.begin(); CI != Calibrations.end(); CI++)
		{
			cout << CI->first <<"\tX\t";
			for(POINTS::const_iterator PI = CI->second.begin(); PI != CI->second.end(); PI++)
				cout << PI->m_dX <<'\t';
			cout << '\n' << CI->first <<"\tY\t";
			for(POINTS::const_iterator PI = CI->second.begin(); PI != CI->second.end(); PI++)
				cout << PI->m_dY <<'\t';
			cout << '\n';
		}
		return 0;
	}
	CALIBRATEDPROBES calibProbes;
	map<CURVEID, string> CurveIDsOutput;
	if(SearchParameter(argc, argv, _FREUNDLICH) != -1)
	{
		CurveIDsOutput[FREUNDLICH_ID] = argv[SearchParameter(argc, argv, _FREUNDLICH) +1 ];
		for(PROBEMEASUREMENTS::const_iterator CI = Calibrations.begin(); CI != Calibrations.end(); CI++)
		{
			valarray<double> vdX(CI->second.size());
			valarray<double> vdY(CI->second.size());
			transform(CI->second.begin(), CI->second.end(), &vdX[0], getX);
			transform(CI->second.begin(), CI->second.end(), &vdY[0], getY);
			LINFITPARAMS Lfp;
			LinearFit(log(vdX), log(vdY), &Lfp);
			if(Lfp.R2 > fR2cutOff)
				calibProbes[CI->first].push_back(PROBECALIBRATION(FREUNDLICH_ID, Lfp));
		}
	}
	if(SearchParameter(argc, argv, _LINEAR) != -1)
	{
		CurveIDsOutput[LINEAR_ID] = argv[SearchParameter(argc, argv, _LINEAR) +1 ];
		for(PROBEMEASUREMENTS::const_iterator CI = Calibrations.begin(); CI != Calibrations.end(); CI++)
		{
			valarray<double> vdX(CI->second.size());
			valarray<double> vdY(CI->second.size());
			transform(CI->second.begin(), CI->second.end(), &vdX[0], getX);
			transform(CI->second.begin(), CI->second.end(), &vdY[0], getY);
			LINFITPARAMS Lfp;
			LinearFit(vdX, vdY, &Lfp);
			if(Lfp.R2 > fR2cutOff)
				calibProbes[CI->first].push_back(PROBECALIBRATION(LINEAR_ID, Lfp));
		}
	}
	if(SearchParameter(argc, argv, _LANGMUIR) != -1)
	{
		CurveIDsOutput[LANGMUIR_ID] = argv[SearchParameter(argc, argv, _LANGMUIR) +1 ];
		for(PROBEMEASUREMENTS::const_iterator CI = Calibrations.begin(); CI != Calibrations.end(); CI++)
		{
			valarray<double> vdX(CI->second.size());
			valarray<double> vdY(CI->second.size());
			transform(CI->second.begin(), CI->second.end(), &vdX[0], getX);
			transform(CI->second.begin(), CI->second.end(), &vdY[0], getY);
			LINFITPARAMS Lfp;
			LinearFit(vdX, vdX/vdY, &Lfp);
			if(Lfp.R2 > fR2cutOff)
				calibProbes[CI->first].push_back(PROBECALIBRATION(LANGMUIR_ID, Lfp));
		}
	}
	for(CALIBRATEDPROBES::const_iterator cpi = calibProbes.begin(); cpi != calibProbes.end(); cpi++)
	{
		vector<PROBECALIBRATION>::const_iterator p = max_element(cpi->second.begin(), cpi->second.end());
		cout << cpi->first << '\t' << p->lfParameters.a << '\t' << p->lfParameters.b << '\t' << p->lfParameters.R2 << '\t' << CurveIDsOutput[p->CurveTypeID]; 
		cout << '\n';
	}
}
	catch(int e)
	{
		if(e == 1)
			cerr << "\nInconsistent parameters";
		return 1;
	}
	return 0;
}

/*
--model: y = a+b*x
SELECT	ID, 
		(Sxx * Sy - Sx * Sxy) / (S * Sxx - Sx * Sx) AS a, 
		(S * Sxy - Sx * Sy) / (S * Sxx - Sx * Sx) AS b, 
		POWER((((Sxy - Y * Sx) - X * Sy) + S * X * Y) / ((sdX * sdY) * (S - 1)), 2) AS R2
FROM    (
		SELECT	ID, 
				COUNT(ID) AS S, 
				SUM(X) AS Sx, 
				SUM(X * X) AS Sxx, 
				SUM(Y) AS Sy, 
				SUM(X * Y) AS Sxy, 
				AVG(X) AS X, 
				AVG(Y) AS Y, 
				STDEV(X) AS sdX, 
                STDEV(Y) AS sdY
		FROM          dbo.fitXY
		GROUP BY ID) AS FitXYParams_1

*/