// wave1.C : example for advanced use of grafix
// using of one-dim wavelets : function wt1 of numerical recipes
// (c) wolf 11/96
// 

#define NRANSI
extern "C" {
#include "nr.h" 
#include "nrutil.h"
}

#include "toolbox.h"

#define NMAX 1024
#define NCEN 333
#define NWID 33

char *help1[] = { 
  "this program shows the approximation of a peak function",
  "with 1-dimensional wavelets (wt1), ","",
  "the button 'itest' uses 'daub4' instead of 'pwt'",
  "the button 'wdisp' shows the wavelets theirselves",
  "button 'spawn' creates a new window of the same shape",
  "the pulldown 'k' allows to select the order of wavelets","",
  "moving the slider of the scrollbar sets different thresholds",
  "the buttons '>' '<' '>>' allow stepwise and continuous changing",
  "for the wavelet coeffs for reconstruction","",
  "the bottom ticks display the indeces of used wavelet coeffs",
  "'dmax' is the maximum approximation error",
  0};


class wave_main : public poll_main {
  scrollbar *scr;
  float *u,*v,*w,*cv;
  int k; // order of wavelets 
  Bool itest; // use daub4 instead of pwt
  Bool wdisp; // show wavelets only 
  
  float frac; // modified with scrollbar
  coord_window *cw;
  info_window *infw;
  window *tickw; // window for ticks  

  void spawn() { new wave_main("spawn",width,height); }  

public:
  wave_main(char *name, int ww = 600, int wh = 300) :
  poll_main(name,ww,wh) {

    new toggle_redraw_button(*mb,"itest",&itest,this);  
    new toggle_redraw_button(*mb,"wdisp",&wdisp,this);
    new instance_button <wave_main> (*mb,"spawn",spawn,this);
   
    u = vector(1,NMAX); // used for determine threshold
    v = vector(1,NMAX); // transformed coeffs
    w = vector(1,NMAX); // original series
    cv = vector(1,NMAX); // back transform
    
    disp_vector **dv = new disp_vector*[3]; // never local ! (spawn fails)
    dv[0] = new disp_vector(w, "w : original vector");
    dv[1] = new disp_vector(cv,"cv : backtransform ");
    dv[2] = 0;
 
    int hc = 60; // leave this room at top for scrollbar
    cw = new nr_coord(*this,NMAX,dv,ww,wh-hc-bh-8,0,hc); 
    cw->define_coord(0,-0.5,NMAX,1.0);  
    tickw = new window(*this,ww,10,1,wh-bh-8,0);
    infw = new info_window(*this,ww,bh,0,wh-bh);    

    char *klist[] = { "4", "12", "20", 0};
    make_radio_menu(*mb,"k",klist,this);
    new help_button(*mb,"help",help1);
    wdisp = False; frac = 20; k = 4; itest = 0; // initials
    for (int i = 1; i <= NMAX; i++)
      w[i] = (i > NCEN-NWID && i < NCEN+NWID ?
	      ((float)(i-NCEN+NWID)*(float)(NCEN+NWID-i))/(NWID*NWID) : 0.0);
    frac = 20;
    scr = new nr_scrollbar(*this,ww,0,bh,9,&frac,this,1,200);

    if (top) main_loop(); else RealizeChildren();
  }

  // action called from radio menu 
  virtual void action(char*, char *value) {
    k = atoi(value); // printf("k %d \n",k);
    solve_it();
  }
  ~wave_main() { 
    free_vector(w,1,NMAX);
    free_vector(v,1,NMAX);
    free_vector(u,1,NMAX);
    free_vector(cv,1,NMAX);
  }

  virtual void solve_it() { // the main callback funct
    // printf("cb %f %d\n",frac,k); 
    // compute coeffs of w -> v, u = |v|
    // could be done only once at init
    if (!itest) pwtset(k); 
    for (int i = 1; i <= NMAX; i++) v[i] = w[i];
    wt1(v, NMAX, 1,itest ? daub4 : pwt); // forward trafo
    for (int i=1; i <= NMAX; i++) u[i] = fabs(v[i]);

    int nfrac = irint(frac);   
    // threshold : use f_select because of name conflict for select
    float thresh = f_select(NMAX-nfrac,NMAX,u); 
    tickw->clear();
    int i,nused = 0;
    for (i = 1; i <= NMAX; i++) {
      if (wdisp) // show wavelets only : inverse transformation of unit vector
	cv[i] = (i == nfrac) ? 1.0 : 0.0; 
      else {
	float av = fabs(v[i]); 
	if (av <= thresh) cv[i] = 0.0; 
	else { cv[i] = v[i]; nused++; tickw->line(i,0,i,5); } 
      }
    }
   
    wt1(cv, NMAX, -1, itest ? daub4 : pwt); // back transform
    if (wdisp) {
      sprintf(infw->info,"display wavelet of order %d index %d",k,nfrac);
    } else {
      float tmp, dmax = 0.0; // get max difference
      for (i = 1; i <= NMAX; i++) if ((tmp=fabs(cv[i]-w[i])) > dmax) dmax=tmp;
      char tmpstr[20]; 
      if (itest) sprintf(tmpstr,"itest"); else sprintf(tmpstr,"k %d",k);
      sprintf(infw->info,"%s NMAX %d nfrac %d nused %d thresh %g dmax %.6f",
	      tmpstr, NMAX,nfrac,nused,thresh,dmax);
    }
    cw->redraw();
    infw->redraw();
  }
  virtual void redraw() { // called from toggle_redraw_button 
    solve_it();
  } 
};

main(int argc, char *argv[]) { 
  new wave_main(argv[0]); 
}