$$\phantom{******** If you see this on the webpage then the browser could not locate *********}$$
$$\phantom{******** jsMath/easy/load.js or the variable root is set wrong in this file *********}$$
$$\newcommand{\vector}[1]{\left[\begin{array}{c} #1 \end{array}\right]}$$ $$\newenvironment{matrix}{\left[\begin{array}{cccccccccc}} {\end{array}\right]}$$ $$\newcommand{\A}{{\cal A}}$$ $$\newcommand{\W}{{\cal W}}$$

/* Copyright 2008-2014 Research Foundation State University of New York   */

/* This file is part of QUB Express.                                      */

/* QUB Express is free software; you can redistribute it and/or modify    */
/* it under the terms of the GNU General Public License as published by   */
/* the Free Software Foundation, either version 3 of the License, or      */
/* (at your option) any later version.                                    */

/* QUB Express is distributed in the hope that it will be useful,         */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of         */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          */
/* GNU General Public License for more details.                           */

/* You should have received a copy of the GNU General Public License,     */
/* named LICENSE.txt, in the QUB Express program directory.  If not, see  */
/* <http://www.gnu.org/licenses/>.                                        */

/*

*/

#ifdef _WIN32
  #include <windows.h>
  #include <time.h>
  #define sleep Sleep
#else
  #include <stdlib.h>
  #define BOOL int
  #define TRUE 1
  #define FALSE 0
#endif

#include <iostream>
using std::cout;
using std::cerr;
using std::endl;
#include <math.h>
#include <string.h>
#include <vector>
#include <map>
using std::vector;
using std::map;
using std::pair;

#include "qubfast.h"
#include "qub_filter.h"
#include "running_mean.h"



// for each stim level we'll average at most the final Ntail responses.
// user must feed chunks of data in the reverse order of the file
// (though within each chunk the samples are in usual forward order)

class qub_dr_equil
{
public:
  int Ntail;
  map<float, RunningMean> stats;
};


extern "C" QUBFAST_API void* qub_doseresponse_equil_init(int Ntail)
{
  qub_dr_equil *dr = new qub_dr_equil;
  dr->Ntail = Ntail;
  return dr;
}

extern "C" QUBFAST_API void qub_doseresponse_equil_free(void *dr)
{
  if ( dr )
    delete (qub_dr_equil *) dr;
}

extern "C" QUBFAST_API void qub_doseresponse_equil_add_data(void *drx, int Ndata, float *stim, float *resp)
{
  qub_dr_equil *dr = (qub_dr_equil*) drx;
  int Ntail = dr->Ntail;
  
  for (int i=Ndata-1; i>=0; --i) {
    RunningMean& stats = dr->stats[stim[i]];
    if ( stats.n < Ntail )
      stats.add(resp[i]);
  }
}

extern "C" QUBFAST_API int qub_doseresponse_equil_result(void *drx, int max_results, float *dose, float *resp, float *resp_std)
{
  qub_dr_equil *dr = (qub_dr_equil*) drx;
  int N = 0;

  vector<float> tail_dose;
  vector<RunningMean*> tailed;
  for (map<float, RunningMean>::iterator ii = dr->stats.begin(); ii != dr->stats.end(); ++ii) {
    if ( ii->second.n == dr->Ntail ) {
      tail_dose.push_back(ii->first);
      tailed.push_back(&(ii->second));
    }
  }
  int skip = ((int) ceil(((float)tailed.size()) * 1.0 / max_results + 0.5)) - 1;
  if ( tailed.size() > 1 ) {
    for (size_t i=0; i<tailed.size(); i += skip + 1) {
      if ( N == max_results )
	break;
      dose[N] = tail_dose[i];
      resp[N] = (float) tailed[i]->mean;
      resp_std[N] = (float) tailed[i]->std();
      ++N;
    }
  } else {
    skip = ((int) ceil(((float)dr->stats.size()) * 1.0 / max_results + 0.5)) - 1;
    if ( skip >= 0 ) {
      for (map<float, RunningMean>::iterator ii = dr->stats.begin(); ii != dr->stats.end(); ++ii) {
	if ( N == max_results )
	  break;
	dose[N] = ii->first;
	resp[N] = (float) ii->second.mean;
	resp_std[N] = (float) ii->second.std();
	++N;
	for (int i=0; (i<skip) && (ii != dr->stats.end()); ++i)
	  ++ii;
      }
    }
  }
  return N;
}



// filter, then for each stim level take the peak
// also find the global minimum
// chunks can be in any order, and should be complete segments

class qub_dr_peak
{
public:
  float sampling_sec;
  float filter_kHz;
  float sign;
  float minval;
  map<float, float> maxvals;
};


extern "C" QUBFAST_API void* qub_doseresponse_peak_init(double sampling_sec, double filter_kHz, double sign)
{
  qub_dr_peak *dr = new qub_dr_peak;
  dr->sampling_sec = (float) sampling_sec;
  dr->filter_kHz = (float) filter_kHz;
  dr->sign = (sign >= 0.0) ? 1.0 : -1.0;
  //cerr << "init: sign = " << sign << " -> " << dr->sign << endl;
  dr->minval = 1e20f;
  return dr;
}

extern "C" QUBFAST_API void qub_doseresponse_peak_free(void *dr)
{
  if ( dr )
    delete (qub_dr_peak *) dr;
}

extern "C" QUBFAST_API void qub_doseresponse_peak_add_data(void *drx, int Ndata, float *stim, float *resp)
{
  qub_dr_peak *dr = (qub_dr_peak*) drx;
  vector<float> filtered;
  float *ff = resp;
  
  if ( dr->filter_kHz ) {
    filtered.resize(Ndata);
    ff = &(filtered[0]);
    qub_filterdataf(resp, ff, Ndata, dr->sampling_sec, (float) (dr->filter_kHz*1e3));
  }
  
  for (int i=0; i<Ndata; ++i) {
    float r = dr->sign * ff[i];
    if ( r < dr->minval )
      dr->minval = r;

    float d = stim[i];
    map<float,float>::iterator maxi = dr->maxvals.find(d);
    if ( maxi == dr->maxvals.end() ) {
      dr->maxvals[d] = r;
    } else if ( maxi->second < r ) {
      maxi->second = r;
    }
  }
}

extern "C" QUBFAST_API int qub_doseresponse_peak_result(void *drx, int max_results, float *dose, float *resp)
{
  qub_dr_peak *dr = (qub_dr_peak*) drx;
  int N = 0;

  int skip = ((int) ceil(((float)dr->maxvals.size()) * 1.0 / max_results + 0.5)) - 1;
  if ( skip >= 0 ) {
    for (map<float, float>::iterator ii = dr->maxvals.begin(); ii != dr->maxvals.end(); ++ii) {
      if ( N == max_results )
	break;
      dose[N] = ii->first;
      resp[N] = dr->sign * (ii->second - dr->minval);
      //cerr << dose[N] << ": " << resp[N] << " = " << dr->sign << " * (" << ii->second << " - " << dr->minval << ")" << endl;
      ++N;
      for (int i=0; (i<skip) && (ii != dr->maxvals.end()); ++i)
	++ii;
    }
  }
  return N;
}