#include <stdio.h>

#include <exception>

#include "haar.h"
#include "line.h"
#include "daub.h"
#include "yahooTS.h"
#include "dbl_sort.h"
#include "histogram.h"
#include "wavethresh.h"

/** \file

  Test wavelet histogramming


  The documentation in this file is formatted for doxygen
  (see www.doxygen.org).

<h4>
   Copyright and Use
</h4>

   You may use this source code without limitation and without
   fee as long as you include:

<blockquote>
     This software was written and is copyrighted by Ian Kaplan, Bear
     Products International, www.bearcave.com, 2002.
</blockquote>

   This software is provided "as is", without any warranty or
   claim as to its usefulness.  Anyone who uses this source code
   uses it at their own risk.  Nor is any support provided by
   Ian Kaplan and Bear Products International.

   Please send any bug fixes or suggested source changes to:

<pre>
     iank@bearcave.com
</pre>

  @author Ian Kaplan

 */

/**
  Print out a vector of doubles, one per line
 */
void pr_vec( const double *v, const size_t N )
{
  for (size_t i = 0; i < N; i++) {
    // printf("%7.4f, %3d\n", v[i], i );
    printf("%7.6f\n", v[i] );
  }
}



/**
  Get a 1-day return time series calculated from a stock price
  (for example, the open or close price).  Here return is

    ret<sub>t</sub> = (p<sub>t</sub> - p<sub>t-1</sub>) / p<sub>t</sub>

  Where p is the stock price and t is in trading days.
  This formula gives the fractional return on a daily basis.

  To calculate N return values, there must be N+1 stock price
  values.

  The function will return true if the return time series is
  calculated without problem, false otherwise.

 */
bool get1DayRetTimeSeries(const char *ticker, 
			  double *ret, 
			  const size_t N )
{
  bool rslt = true;
  yahooTS ts( "..\\data\\equities\\" );

  size_t n = N+1;
  double *price = new double[ N+1 ];

  // fetch N+1 close price values
  ts.getTS( ticker, price, n, yahooTS::Close );
  if (n == N+1) {
    for (size_t i = 1; i < N+1; i++) {
      ret[i-1] = (price[i] - price[i-1])/price[i-1];
    }
  }
  else {
    rslt = false;
  }

  delete [] price;

  return rslt;
} // get1DayRetTimeSeries




/**
  Return the absolute value of "val"
 */
double abs( const double val )
{
  double new_val = val;

  if (val < 0)
    new_val = -val;

  return val;
} // abs



int
main()
{
  const size_t N = 512;
  double ret[ N ];

  if (get1DayRetTimeSeries( "ibm", ret, N )) {

    const size_t num_binz = N;
    histogram::bin_vec binz( num_binz );
    histogram histo;
    size_t i;

    histo.calculate( ret, N, binz );

    // Normalized haar wavelet
    // haar<histogram::bin_vec> w;

    // Daubechies D4 wavelet
    // Daubechies<histogram::bin_vec> w;

    // Unnormalized Linear interpolation wavelet
    line<histogram::bin_vec> w;

    w.forwardTrans( binz, N );

    const size_t start = 16;
    const size_t size = N - start;
    wavethresh<histogram::bin_vec>::thresh( binz, start, size );

    w.inverseTrans( binz, N );

    // convert frequency to probability
    for (i = 0; i < num_binz; i++) {
      binz[i] = binz[i] / static_cast<double>( N );
    }


    for (i = 0; i < num_binz; i++) {
	printf("%7.4f %7.4f\n", binz.start(i), binz[i] );
    }

  }

  return 0;
}