// calibrate.cpp : Defines the entry point for the console application.
//

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

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 ValE
{
	double dVal;
	double dErr;
	ValE(double v, double err): dVal(v), dErr(err) {}
};

class Calibrator
{
public:
	Calibrator() {}
	virtual ValE X(ValE y) = 0;
};
typedef Calibrator* PCalibrator;

class Langmuir : public Calibrator
{
	double m_dYmax;
	double m_dK;
public:
	Langmuir(double dYmax, double dK) : m_dYmax(dYmax), m_dK(dK) {}
	virtual ValE X(ValE y) {return ValE(y.dVal / (m_dYmax - y.dVal) / m_dK, m_dYmax / (m_dYmax - y.dVal) * y.dErr);} 
};

class Freundlich : public Calibrator
{
	double m_da;
	double m_db;
public:
	Freundlich(double a, double b) : m_da(a), m_db(b) {}
	virtual ValE X(ValE y) {return ValE(pow(y.dVal / m_da, 1 / m_db), y.dErr / m_db);} 
};

class Linear : public Calibrator
{
	double m_da;
	double m_db;
public:
	Linear(double a, double b) : m_da(a), m_db(b) {}
	virtual ValE X(ValE y) {return ValE((y.dVal - m_da) / m_db, y.dErr*y.dVal/(y.dVal - m_da));} 
};


int main(int argc, char* argv[])
{
try
{
#if _DEBUG
	__asm{int 3}
#endif
	typedef map<string, CURVEID> NAMEDCURVES;
	typedef map<string, PCalibrator> PROBESENSOR;
	if(argc == 1)
	{
		cerr << "concentration  < calibrated_probes.txt [-freund {F}] [-linear {LN}] [-langmuir {LG}] -intens sample1.txt sample2.txt...  > out.txt\n";
		cerr << "Sample input format:  ExperimentID probeID signalIntens stderror\n";
		cerr << "Calibration input format:  ProbeID a b R2 CurveID (see below)\n";
		cerr << "-freund {F}\tFreundlich calibration y=exp(a)*x^b; input marker, e.g. F, provided by user\n";
		cerr << "-linear {LN}\tLinear calibration y = a+b*x; input marker, e.g. LN, provided by user\n";
		cerr << "-langmuir {LG}\tLangmuir calibration y = Y_max*K*x/(1+K*x); input marker, e.g. LG, provided by user\n";
		cerr << "a = 1/(K*Y_max); b = 1/Y_max\n";
		return -1;
	}
	
	NAMEDCURVES CurveIDs;

	if(SearchParameter(argc, argv, _FREUNDLICH) != -1)
		CurveIDs[argv[SearchParameter(argc, argv, _FREUNDLICH) +1] ] = FREUNDLICH_ID;
	if(SearchParameter(argc, argv, _LANGMUIR) != -1)
		CurveIDs[argv[SearchParameter(argc, argv, _LANGMUIR) +1] ] = LANGMUIR_ID;
	if(SearchParameter(argc, argv, _LINEAR) != -1)
		CurveIDs[argv[SearchParameter(argc, argv, _LINEAR) +1 ] ] = LINEAR_ID;

	PROBESENSOR AllSensors;
	while (!cin.eof())
	{
		string sProbeID, sCurveID;
		double a,b,r2;
		cin >> sProbeID >> a >> b >> r2 >> sCurveID;
		if(sProbeID == "")
			continue;

		switch (CurveIDs[sCurveID])
		{
		case LANGMUIR_ID:
			AllSensors[sProbeID] = new Langmuir(1/b, b / a);
			break;
		case FREUNDLICH_ID:
			AllSensors[sProbeID] = new Freundlich(exp(a), b);
			break;
		case LINEAR_ID:
			AllSensors[sProbeID] = new Linear(a, b);
			break;
		}
	}
	for(int i=SearchParameter(argc, argv, _INTENS)+1; i < argc; i++)
	{
		ifstream ifsRecords(argv[i]);
		if(ifsRecords.bad())
			throw 2;
		while(!ifsRecords.eof())
		{
			string sHybID, sProbeID;
			double dAvIntens, dIntensErr;
			ifsRecords >> sHybID >> sProbeID >> dAvIntens >> dIntensErr;
			if (sHybID == "")
				continue;
			PROBESENSOR::const_iterator PSI = AllSensors.find(sProbeID);
			if(PSI == AllSensors.end())
				continue;
			ValE conc = PSI->second->X(ValE(dAvIntens, dIntensErr));
			cout << sHybID << '\t' << sProbeID << '\t' << conc.dVal << '\t' << conc.dErr << '\n';
		}
	}

}
	catch(int e)
	{
		if(e == 1)
			cerr << "\nInconsistent parameters";
		return 1;
		if(e == 2)
			cerr << "\nFile error";
		return 2;
	}
	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

*/