$$\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-2011 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 "qub_adaptive_resample.h"
#include "qub_mean.h"
#include "qub_eventmaker.h"

typedef struct {
  int Cls, N;
} NCls;

typedef struct {
  int First, Last;
  double Mean, Std;
} ResIval;

class ResTree {
public:
  int First, Last;
  double Mean, Std;
  ResTree *Left, *Right;
  
  ResTree(int aFirst, int aLast)
    : First(aFirst), Last(aLast), Left(NULL), Right(NULL)
  {}
  ~ResTree() {
    if ( Left )   delete Left;
    if ( Right )  delete Right;
  }
  
  void split() {
    // Assert((Left = nil) and (Right = nil), 'TAdResTree.Split: not a leaf node');
    int count = Last - First + 1;
    splitBefore(First + (count / 2));
  }
  void splitBefore(int RightFirst) {
    Left  = new ResTree(First, RightFirst-1);
    Right = new ResTree(RightFirst, Last);
  }
  void calcMeanStd(float* data) {
    int n = Last - First + 1;
    qub_mean_rss(data+First, n, &Mean, &Std);
    if ( n <= 1 )
      Std = 0.0;
    else
      Std = sqrt(Std / (n-1));
  }
  void splitStd(float* data, double MaxStd) {
    if ( ! Left ) {
      calcMeanStd(data);
      if ( Std > MaxStd )
	split();
    }
    if ( Left ) {
      Left->splitStd(data, MaxStd);
      Right->splitStd(data, MaxStd);
    }
  }
};

class ResTreeDFS {
private:
  bool LeavesOnly;
  vector<ResTree*> ToVisit;
  
  void push(ResTree* node) { ToVisit.push_back(node); }
  ResTree* pop() {
    ResTree* result = ToVisit.back();
    ToVisit.pop_back();
    return result;
  }

public:
  ResTree* Here;
  
  ResTreeDFS(ResTree* Root, bool aLeavesOnly)
    : LeavesOnly(aLeavesOnly), Here(NULL)
  {
    push(Root);
  }
  bool next() { // returns false if no more
    bool result = ToVisit.size();
    if ( result ) {
      Here = pop();
      if ( Here->Left ) {
	push(Here->Right);
	push(Here->Left);
	if ( LeavesOnly )
	  result = next();
      }
    }
    return result;
  }
};

//******************************************************************************

extern "C" QUBFAST_API int qub_adaptive_resample(float *data, int ndata, double MaxStd,
                                                 float *out_mean, float *out_std,
						 int *out_first, int *out_last, int *out_closest)
  // Fills out_* with a reduced form of data, as intervals with all out_std[i] <= MaxStd (binary decomposition).
  // @param data: sampled input array
  // @param ndata: length of data array
  // @param MaxStd: desired ceiling for out_std[i]
  // @param out_mean: array of each interval's mean value; must have room for up to ndata intervals
  // @param out_std: array of each interval's standard deviation
  // @param out_first: array of each interval's first sample index, relative to start of data
  // @param out_last: array of each interval's final sample index
  // @param out_closest: array of closest sample index to mean in each interval

  // @returns: number of intervals
{
  ResTree *tree = new ResTree(0, ndata-1);
  tree->splitStd(data, MaxStd);
  
  int out_i=0;
  ResTreeDFS *trav = new ResTreeDFS(tree, true);
  while ( trav->next() ) {
    out_first[out_i] = trav->Here->First;
    out_last[out_i]  = trav->Here->Last;
    out_mean[out_i]  = (float) trav->Here->Mean;
    if ( out_std )
      out_std[out_i]   = (float) trav->Here->Std;
    if ( out_closest) {
      out_closest[out_i] = out_first[out_i];
      double closest_diff = fabs(data[out_first[out_i]] - out_mean[out_i]);
      for (int j=out_first[out_i]; j<=out_last[out_i]; ++j) {
	double diff = fabs(data[j] - out_mean[out_i]);
	if ( diff < closest_diff ) {
	  out_closest[out_i] = j;
	  closest_diff = diff;
	}
      }
    }
    ++out_i;
  }
  delete trav;
  delete tree;
  return out_i;
}


//******************************************************************************


class qub_adres
{
public:
  double amp_delta;
  double min_dur_ms;
  double sampling_ms;
  QUB_EventMaker_ByGauss eventMaker;

  qub_adres(double _amp_delta, double _min_dur_ms, double _sampling_ms)
    : amp_delta(_amp_delta), min_dur_ms(_min_dur_ms), sampling_ms(_sampling_ms)
  {}
  
  ~qub_adres()
  {}
  
  void add(float *samples, int sample_count)
  {
    vector<float> means(sample_count), stds(sample_count);
    vector<int> ff(sample_count), ll(sample_count);
    int Nevent = qub_adaptive_resample(samples, sample_count, amp_delta, & means[0], & stds[0], & ff[0], & ll[0], NULL);
    for (int i=0; i<Nevent; ++i) {
      eventMaker.add_event(means[i], ((float)(sample_count-1))*squared(stds[i]), ll[i]-ff[i]+1);
    }
  }
  
  void done_add()
  {
    eventMaker.done_add();
    eventMaker.combine_all(amp_delta);
  }
  
  int get_class_count()
  {
    return eventMaker.class_count();
  }
  
  int get_dist(double *amps, double *stds)
  {
    eventMaker.get_dist(& amps[0], & stds[0]);
    return eventMaker.class_count();
  }
  
  int get_next_dwells(int offset_in_data, int sample_count, int *firsts, int *lasts, int *classes, float *durations)
  {
    int ndwt = eventMaker.get_next_dwells(offset_in_data, sample_count, firsts, lasts, classes);
    if ( min_dur_ms > 0.0 ) {
      ndwt = delete_short_events(ndwt, firsts, lasts, classes, int(ceil(min_dur_ms / sampling_ms)));
    }
    if ( durations ) {
      for (int i=0; i<ndwt; ++i)
	durations[i] = (float) ((lasts[i] - firsts[i] + 1) * sampling_ms);
    }
    return ndwt;
  }
};


QUBFAST_API void*  qub_adres_create(double amp_delta, double min_dur_ms, double sampling_ms)
{
  return new qub_adres(amp_delta, min_dur_ms, sampling_ms);
}

QUBFAST_API void   qub_adres_destroy(void *adres)
{
  delete (qub_adres*) adres;
}

QUBFAST_API void   qub_adres_add(void *adres, float *samples, int count)
{
  ((qub_adres*) adres)->add(samples, count);
}

QUBFAST_API void   qub_adres_done_add(void *adres)
{
  ((qub_adres*) adres)->done_add();
}

QUBFAST_API int    qub_adres_get_class_count(void *adres)
{
  return ((qub_adres*) adres)->get_class_count();
}

QUBFAST_API int    qub_adres_get_dist(void *adres, double *amps, double *stds)
{
  return ((qub_adres*) adres)->get_dist(amps, stds);
}

QUBFAST_API int    qub_adres_get_next_dwells(void *adres, int offset_in_data, int sample_count,
					       int *firsts, int *lasts, int *classes, float *durations)
{
  return ((qub_adres*) adres)->get_next_dwells(offset_in_data, sample_count, firsts, lasts, classes, durations);
}