/*****************************************************************************
File:           fast_v3.c   November 9
Authors:        Ioan Tabus   (tabus@cs.tut.fi)
		Bogdan Dumitrescu
                Doina Petrescu
                Alin Alecu


Purpose:        Fast stack filter design


Copyright 1998   Ioan Tabus, Bogdan Dumitrescu, Doina Petrescu, Alin Alecu. 
		All Rights Reserved.

Reference:
@ARTICLE{TFI,
   author = { I. T\u{a}bu\c{s} and D. Petrescu and M. Gabbouj  },
   title = {A Training Framework for  Stack and {Boolean} Filtering-{F}ast optimal design procedures and robustness case study},
   journal = TRIP # " Special Issue on Nonlinear Image Processing",
   volume = "IP-5",
   year = "1996",
   month = "June",
   pages = "809--826",
}


These programs are supplied free of charge for research purposes only,
and may not be sold or incorporated into any commercial product.  There is
ABSOLUTELY NO WARRANTY of any sort, nor any undertaking that they are
fit for ANY PURPOSE WHATSOEVER.  Use them at your own risk.  If you do
happen to find a bug, or have modifications to suggest, please report
it to  tabus@cs.tut.fi.  The copyright notice above and this statement of conditions 
must remain an integral part of each and every copy made of this file.

/****************************************************************************/
/*                                                                          */
/*                  Fast (sub)optimal stack filter design                   */ 
/*                                                                          */
/*  Inputs:                                                                 */
/*    - name of coefficients file  (e.g. as written by coeff1.c)            */
/*    - name of output filter file                                          */
/*  Output:                                                                 */ 
/*    - filter file  (it will be written in boolean tabular form)           */
/*                                                                          */
/*  Allows insertion of: - desired window size                              */
/*                                                                          */
/*                                                                          */
/****************************************************************************/

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>

#define REC_COMP_COST		/* gains are computed recursively in Hasse */
				/* keep it active !!!			   */

#define COMBINE_SONS		/* make an extra step trying to combine sons
				   of a node, in order to improve MAE.     */
				/* Small improvement of MAE...		   */

/*#define RUN_LP*/		/* try to run LP                           */
				/* disactivate for large window size !     */


#define N  21   /* Max. size of window.
                      Eg:  5   9   10     17     19       21  */
 
#define LG 2097152  /* Max. no. of coefficients.
                      Eg: 32 512 1024 131072 524288  2097152 */

#define NUM 258064

#define LN 300   /* dimension of cost matrix                                */
#define LNN 1000   


typedef unsigned char byte;

byte f[LG], fp[LG], fpocket[LG], fnon[LG];

int costs[LG]; /* Cost coefficients read from disk                          */

int C0;
double J0,J2,Jmin,J,Jbool[N+1],Jstack[N+1];

#ifdef RUN_LP
double  Rcost[LN]; /* Costs matrix, used in LP                              */
#endif

int NR,   /* Current window size                                            */
    LG1;

clock_t cl1;

double evalJc(); /* Evaluates Jstack                                        */
void comp_boool_stack(); /* Computes the Boolean and stack filter           */ 
byte comp_LP();  /* Procedure which implements LP                           */
void calculate(); /* Computes lg, ln and lnn                                */

int winsize, /* Actual window size                                          */
    lg,      /* 2^winsize                                                   */
    ln,      /* lg+2                                                        */
    lnn;     /* winsize*lg/2 + 1                                            */

int  hasse[LG];	/* Level-ordered undecided nodes in Hasse diagram           */

#ifdef REC_COMP_COST
byte visited[LG]; /* used for proper recursion in the Hasse diagram         */
#endif

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


void calculate(int size)
{
  lg=pow(2,size);
  ln=lg+2;
  lnn=(size*lg/2) + 1;
}


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


int	main( argc, argv)
int argc;
char * argv[];

{

   FILE *fpt;
   char namec[40], namef[40],str[10];


   int nn,i,j,j1,j2,jocker,depl,p,it,imin,sw,t,st,best,v[N],cn,tmp;
   int jmax,k;
   double s,pivot;


/* Input window size                                                       */
        

   printf("\n Input window size: ");
   scanf("%d",&winsize);
   calculate(winsize);

/* Reading cost coefficients from disk file                                */

    if(argc==3) strcpy(namec, argv[1]);
    else { printf("\n Name of coefficients file: ");
           scanf("%s",namec);
         }
    fpt=fopen(namec,"r");
    if(fpt==NULL) {printf("\n Can't open %s\n", namec); return 0;}

    if(argc==3) strcpy(namef,argv[2]);
    else { printf("\n Name of filter file: ");
           scanf("%s", namef);
         }
    for(nn=0;nn<lg;nn++)
     { if(fscanf(fpt,"%d",&costs[nn])!=1)
          { fclose(fpt); printf("\n Error reading %s\n",namec); 
            exit(-1);
          }
      }

    fscanf(fpt,"%d",&C0); 

    fclose(fpt);



          /***************************
         *     Main Recursions     * 
         ***************************/


    cl1=clock();


    LG1=lg;

    for(NR=winsize;NR<=winsize;NR++)
      {

       printf("\n *****  NR=  %d  *****\n",NR);

/* Compute boolean and stack filter                                         */

       comp_boool_stack();

/* Evaluate Jstack                                                          */
       Jstack[NR]=evalJc(fp);
       printf(" Jpos = %.5lf \n ",Jstack[NR]);
       printf(" time  %.3lf \n",((double)clock()-cl1)/CLOCKS_PER_SEC);  
       LG1=LG1*2;

      }

      LG1=lg;


/* Now try to solve LP problem                                              */
    #ifdef RUN_LP
     if(comp_LP()){printf("\n LP solution could not be computed \n");}
    #endif


  printf("\n  time  %.3lf\n",((double)clock()-cl1)/CLOCKS_PER_SEC);  


/* Write filter to disk file                                                */

     fpt=fopen(namef,"w");
     if(fpt==NULL) {printf("\n Error writing filter\n"); return 0;}
     for(i=0;i<lg;i++)
       { fprintf(fpt,"%d \n",fp[i]);
       }                              
     fclose(fpt); 

}


/****************************************************************************/
#ifdef REC_COMP_COST

int	comp_cost( int nod )
{
  int c, v, l, fiu;

  c = costs[nod];
  visited[nod] = 1;

  v = 1;
  for ( l=0 ; l < winsize ; l++ )
  {
    fiu = nod | v;
    if ( (fiu > nod) && ! fpocket[fiu] && ! visited[fiu] )
      c += comp_cost(fiu);
    v <<= 1;
  }

  return c;
}

void	reset_visited( int nod )
{
  int v, l, fiu;

  visited[nod] = 0;

  v = 1;
  for ( l=0 ; l < winsize ; l++ )
  {
    fiu = nod | v;
    if ( (fiu > nod) && ! fpocket[fiu] )
      reset_visited(fiu);
    v <<= 1;
  }

}

void	set_sons_decided(int nod, int no_fp)
{
  int v, l, fiu;

  fpocket[nod] = 1;
  if ( ! no_fp )
    fp[nod] = 1;

  v = 1;
  for ( l=0 ; l < winsize ; l++ )
  {
    fiu = nod | v;
    if ( (fiu > nod) && ! fpocket[fiu] )
      set_sons_decided(fiu, no_fp);
    v <<= 1;
  }
}

#endif /* #ifdef REC_COMP_COST */

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

 /* Compute the boolean and stack filter */                          
void comp_boool_stack()

{
 int i,j,j1,j2,depl,it,imin,sw,t,st,best,ifrez,Improv;
 double Jdelta,J1;
 
  int 		nniv[N+1], cniv[N+1], nod, fiu, v, nnivc, l, k, ii, kk, ll, vv;
  int		castig1;
 #ifdef RUN_LP
  char		fpocket1[LG];
  int		costs1[LG];
 #endif




 /************
* Stage II  *
*************/

/* f will store decided nodes, namely f[i]=1 for cost[i]<0 and f[i]=0 for   */
/* cost[i]>0                                                                */
/* fnon will store undecided nodes                                          */
                          
 for(i=0;i<LG1;i++)
   {  fnon[i]=0;
      if(costs[i]>0) {f[i]=0;}
          else {f[i]=1; if (costs[i]==0)  fnon[i]=1; }
    }

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

/* Evaluate J                                                               */

 Jbool[NR]=evalJc(f);
   printf(" Jbool %.5lf \n ",Jbool[NR]);   


/* Test if the Boolean filter is stack                                      */


/***********************
* Stage III-1 + III-2  *
***********************/
 best=1;

 for(i=LG1-1;i>=0;i--)
   { fpocket[i]=0;fp[i]=0;

/* if f[node]=1 (decided on 1), verify all nodes that stack under. If we   */
/*  find one that has f[node']=0, then stack condition not respected       */

     if(f[i])
     { fp[i]=1;fpocket[i]=1;
       depl=1;
       for(j=0;j<NR;j++)
         { j2=(depl | i); 
           depl<<=1;
           if(j2>i) {if(!fp[j2]){fp[i]=0;fpocket[i]=0;best=0;break;} } 
                         /*j2 lays under i*/
         }
     }
   }
  

/* exit with fp=set I1                                                     */

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


/* Evaluate J                                                              */

 J0=evalJc(fp);
 printf("  J(f^+_max) %.5lf \n ",J0);

 
 if(best==0) 
   {
    printf(" Boolean isn't stack!\n");
                         

/*******************
* Stage III-3       *
*******************/

/* for all nodes node_i                                                    */


  for(i=0;i<LG1;i++)     /* Look for compulsory fp=0 */
   {    


/* if it is decided on 0 or it is undecided, if we find a node that stacks */
/* over and is decided on 1, then clearly node_i does not belong to the    */
/* set I0                                                                  */
     if((!f[i])||(fnon[i]))
     { depl=1;fpocket[i]=1;
       for(j=0;j<NR;j++)
         { 
           j2=(~(depl) & i); 
           depl<<=1;
           if(j2<i) { if((f[j2]&&(!fnon[j2]) ) || (!fpocket[j2])) 
             { fpocket[i]=0; }} /* j2 is over j1  */
         }
     }
   }

  /* ifrez=0; */

  /* Exit with fpocket = I0 */
  /*******************/

  /*ifrez will give us how many undecided nodes we have                    */
  /*for(i=0;i<LG1;i++) { if(fpocket[i]==0) ifrez++;}
  printf(" ifrez %d \n ",ifrez);*/
    it=0; 


/*******************
* Stage III-5       *
*******************/

   /* for (k=0 ; k < LG1 ; k++ )
     if ( ! fpocket[k] && fnon[k] )
       printf("%d ", fp[k]); */

  #ifdef RUN_LP
   for (k=0 ; k < LG1 ; k++ )
     costs1[k] = costs[k];
  #endif

  #ifdef RUN_LP
   for ( kk=0 ; kk < 2 ; kk++ )			/* try to decide 1, then to 0 */
  #else
   for ( kk=0 ; kk < 1 ; kk++ )			/* try to decide only to 1 */
  #endif
   {
   #ifdef RUN_LP
    if ( kk )
    {

      for ( i=0 ; i < LG1 ; i++ )
        fpocket1[i] = fpocket[i];

      for ( i=0 ; i < LG1 ; i++ )
      {
	j = i ^ (LG1-1);
        costs[j] = - costs1[i];
        fpocket[j] = fpocket1[i];
      }
    }
   #endif

    for ( k=0 ; k <= winsize ; k++ )
      nniv[k] = 0;

    /* first pass: compute number of undecided nodes per level */
    for ( i=0 ; i < LG1 ; i++ )
    {
      if ( ! fpocket[i] )
      {
        v = 0;
        ii=i;
        for ( j=0 ; j < winsize ; j++ )
        {
          v += (ii & 1);
          ii >>= 1;
        }
        ++nniv[v];
      }
    }

    /* start address in hasse for levels */
    cniv[winsize] = 0;
    for ( k=winsize-1 ; k >=0 ; k-- )
    {
      cniv[k] = cniv[k+1] + nniv[k+1];
      /*printf("cniv(k): %d\n", cniv[k] );*/
    }

    /* second pass: put nodes in appropriate positions in hasse */
    for ( i=0, l=0 ; i < LG1 ; i++ )
    {
      if ( ! fpocket[i] )
      {
        v = 0;
        ii=i;
        for ( j=0 ; j < winsize ; j++ )
        {
          v += (ii & 1);
          ii >>= 1;
        }
	hasse[ cniv[v]++ ] = i;
 	++l;
      }
    }

    printf("kk: %d    undecided: %d\n", kk, l);


    while (1)
      {
       it++;
       J0=evalJc(fp);
       printf("  it= %d    J(improv) %.5lf ",it,J0);
       printf(" time  %.3lf  \n",((double)clock()-cl1)/CLOCKS_PER_SEC);  


      
      Improv=0;


      nnivc = 0;
      for ( k=winsize ; k >=0 ; nnivc += nniv[k], k-- )
      {
	if ( nniv[k] )
	{
	  for ( j=nnivc ; j < nniv[k] + nnivc ; j++ )
	  {
	    nod = hasse[j];
	    if ( ! fpocket[nod] )
	    {
	      castig1 = costs[nod];

	      /* if current node is a candidate for improvment */
	      if ( castig1 <= 0 )
	      {
	       #ifdef REC_COMP_COST
		reset_visited(nod);
		castig1 = comp_cost( nod );
	       #else
                for ( l=0 ; l < j ; l++ )
                {
                  fiu = hasse[l];
                  if ( (fiu | nod) == fiu && ! fpocket[fiu] )
		    castig1 +=  costs[fiu];
                }
	       #endif
              }

	      /* if there is an actual improvment */
	      if ( castig1 < 0 )
	      {
		Improv = 1;

	       #ifdef REC_COMP_COST
		set_sons_decided(nod, kk);
	       #else
                for ( l=0 ; l <= j ; l++ )
	  	{
		  fiu = hasse[l];
                  if ( (fiu | nod) == fiu )
		  {
		    fpocket[fiu] = 1;
		    if ( !kk )
		      fp[fiu] = 1;
		  }
		}
	       #endif
	      }
	    }
	  }
	}
      }


      /* if(  (!Improv) | ( (J0-evalJc(fp))<0.0000001*J0 )  ) break; */
      if(  !Improv ) break;
              
    } /* end while (1) */
   }


  /* try sons combinations to improve MAE */

  #ifdef COMBINE_SONS
   nnivc = 0;
   for ( k=winsize ; k >=0 ; nnivc += nniv[k], k-- )
   {
     if ( nniv[k] )
     {
       for ( j=nnivc ; j < nniv[k] + nnivc ; j++ )
       {
         nod = hasse[j];
         reset_visited(nod);
 
         castig1 = 0;
         for ( l=0, v=1 ; l < winsize ; l++, v<<=1 )
         {
           fiu = nod | v;
           if ( fiu > nod && !fpocket[fiu] )
           {
             castig1 += comp_cost(fiu);
             if ( castig1 < 0 )         /* the good combination of sons */
             {
               /* set decided all sons and their descendents */
               for ( ll=0, vv=1 ; ll <= l ; ll++, vv<<=1 )
               {
                 fiu = nod | vv;
                 if ( fiu > nod && !fpocket[fiu] )
                   set_sons_decided(fiu, 0);
               }
 
               reset_visited(nod);
               castig1 = 0;
             }
           }
         }
       }
     }
   }
  #endif /* #ifdef COMBINE_SONS */
 


  #ifdef RUN_LP
   for (k=0 ; k < LG1 ; k++ )
   {
     fpocket1[k] = fpocket[k^(LG1-1)];
     costs[k] = costs1[k];
   }

   for (k=0 ; k < LG1 ; k++ )
     fpocket[k] = fpocket1[k];
  #endif

 /* Test if the  filter is stack                                           */
 /*best=1;
 for(i=lg-1;i>=0;i--)
   { 
     if(fp[i])
     {
       depl=1;
       for(j=0;j<winsize;j++)
         { j2=(depl | i); 
           depl<<=1;
           if(j2>i) {if(!fp[j2]){ best=0;break;} } 
         }
     }
   }
  
  if(best==0) 
    printf(" Stack is not a stack!\n");
  else
    printf(" Stack is stack!\n"); */


  }

 else {
         for(i=0;i<LG1;i++) fp[i]=f[i];
      }
}


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

/* Evaluate J(c,cost) starting from  canonical terms                       */
double evalJc( cof)
unsigned char cof[];


{
 unsigned i1;
 double Jloc;

 Jloc=C0;
 for(i1=0;i1<LG1;i1++)
   if(cof[i1]) Jloc+=costs[i1];
 return(Jloc/NUM);
}

 

 
/***************************************************************************/
#ifdef RUN_LP

/* Compute LP */
byte comp_LP()

{
   
byte  fpo[LG],cannot,best;
double costlp[LG];
char R[LNN][LN]; 
int LGG1,LNN1,LN1,
    coresp[LG],/* contains undecided nodes                                 */
    cor[LG];   /* cor[i] contains the position j for coresp[j]=node i      */
unsigned zerov[LG], /* store nonbasic variables for simplex                */
         basis[LNN]; /* store basic variables for simplex                  */
unsigned f1[LG],f2[LG];
 int nn,i,j,j1,j2,jocker,depl,p,it,tmp;
 int jmax,k;
 double s,pivot;


 
 /* Enter the simplex with fpocket=0 - undecided node                      */
 /*                        fpocket=1 - decided node                        */

 /* construct coresp and cor                                               */
 LGG1=0;
 for(j=0;j<lg;j++)
    /* if(!fpocket[j] && fp[j]==0) */
    if (!fpocket[j] ) 
    {
      if ( LGG1 < LN )
        coresp[LGG1]=j;
      cor[j]=LGG1;
      LGG1++;
      /* if(LGG1>=LN)   {cannot=1;return(cannot);} */
    }

 printf("LGG1: %d\n", LGG1);

 if(LGG1>=LN)   {cannot=1;return(cannot);}

/* Construct R */

/* for every ni belonging to coresp, search for every node nj that stacks  */
/* under. If fpocket[nj]=1 (decided node) then R[k,i]=1 and R[k,LGG1]=1 ,  */
/* after which we select another ni                                        */
/*        If fpocket[nj]=0 (undecided node) then R[k,i]=1, R[k,j]=-1,      */
/*        R[k,LGG1]=0 after which go to next nj                            */
 LNN1=0;
 for(i=0;i<LGG1;i++)
   { j1=coresp[i];   /* j1 and i are the principal  characters  */
     jocker=1;
     depl=1;
     for(j=0;j<winsize;j++)
       { j2=depl|j1; depl<<=1;
         if(j2!=j1)   /* j2 lays under j1  */
           {  
             if(fpocket[j2])
                     {if(jocker)
                         { for(p=0;p<LGG1;p++) R[LNN1][p]=0;
                           R[LNN1][LGG1]=1;
                           R[LNN1][i]=1;
                           LNN1++;
                           if(LNN1>=LNN)   {cannot=1;return(cannot);}
                           jocker=0;
                         }
                      }
             else
               { for(p=0;p<=LGG1;p++) R[LNN1][p]=0;
                 R[LNN1][cor[j2]]=-1;
                 R[LNN1][i]=1;
                 if(LNN1>=LNN)   {cannot=1;return(cannot);}
                 LNN1++;
               }
           }
       }
   }
LN1=LGG1+2;



/* Construct zerov                                                         */
 for(it=0;it<LGG1;it++)
   zerov[it]=it;

 /* Construct basis                                                        */
 for(j=0;j<LNN1;j++)
   basis[j]=LGG1+j;


 /* Reduce fp so as to contain only decisions for ni belonging to Iu       */
 /* (undecided class). Modify Rcost in the same manner                     */
for(k=0;k<lg;k++) fpo[k]=fp[k];
for(k=0;k<LGG1;k++)  
    {
    fp[k]=fpo[coresp[k]];
    Rcost[k]=-costs[coresp[k]];
    }
 Rcost[LGG1]=0.;
 Rcost[LGG1+1]=0.;




 /* initialize simplex */
 /* for all ni belonging to undecided class, if fp[ni]=1 then find in R a   */
 /* row on which we have on column k a 1 and on the last column a 0 (the    */
 /* form fk-fj <=0 ). If we don't find any, search for a row on which we    */
 /* have on column k a 1 and on the last column a 1 (the form fk<=1)        */
 /* when we find one of these two types of rows, permute node k from zerov  */
 /* with the node from basis corresponding to the row we found.             */
 /* Add a column to R, and write 1 in this column on the position           */
 /* corresponding to that row. For all rows except the selected one, if on  */
 /* position k we have a 1, subtract from that row the initial row we found.*/
 /* If, on the other hand, we have a -1 on position k, add to  that row     */
 /* the initial row.                                                        */ 

for(k=LGG1-1;k>=0;k--)   
   { if(fp[k]==1)
       { jmax=-1;
         for(j=0;j<LNN1;j++)
           if((R[j][k]==1)&&(R[j][LGG1]==0)) {jmax=j;break;} 
         if(jmax<0)
          {  for(j=0;j<LNN1;j++)
               if((R[j][k]==1)&&(R[j][LGG1]==1)) {jmax=j;break;}
          }
        if(jmax<0) {printf("\n You are degenerate!!!\n"); return;}
        tmp=zerov[k];
        zerov[k]=basis[jmax];
        basis[jmax]=tmp;
        for(j=0;j<LNN1;j++)
            R[j][LGG1+1]=0;
        R[jmax][LGG1+1]=1;
        Rcost[LGG1+1]=0.;
        for(j=0;j<LNN1;j++)
         if(j!=jmax)
           { if(R[j][k]==1) for(i=0;i<LN1;i++) R[j][i]-=R[jmax][i];
             else if(R[j][k]==-1) for(i=0;i<LN1;i++) R[j][i]+=R[jmax][i]; 
           }
        pivot=Rcost[k];
        for(i=0;i<LN1;i++) 
                           Rcost[i]-=pivot*R[jmax][i];
        for(i=0;i<LNN1;i++) R[i][k]=R[i][LGG1+1];
        Rcost[k]=Rcost[LGG1+1];
       }
   }

/* Reconstruct the function                                                 */
  for(i=0;i<lg;i++) fp[i]=fpo[i];  

  for(i=0;i<LGG1;i++)
     if (zerov[i]<LGG1) fp[coresp[zerov[i]]]=0;
  for(i=0;i<LNN1;i++)
     if(basis[i]<LGG1) fp[coresp[basis[i]]]=R[i][LGG1];
 

  /* Evaluate J                                                             */

 printf("\n LNN1=%d  LGG1=%d  J=%.5lf \n",LNN1,LGG1,evalJc(fp));

 while(1)   /*  Optimization */
   { 
     k=0;
     s=Rcost[0];
     for(i=1;i<LGG1;i++)
       if(Rcost[i]>s) {k=i;s=Rcost[i];}
     
     if(Rcost[k]<=0) break;   /*  All done!  */

     /* identical algorithm as that of initialization */
     else 
       { jmax=-1;
         for(j=0;j<LNN1;j++)
           if((R[j][k]==1)&&(R[j][LGG1]==0)) {jmax=j;break;} 
         if(jmax<0)
          {  for(j=0;j<LNN1;j++)
               if((R[j][k]==1)&&(R[j][LGG1]==1)) {jmax=j;break;}
          }
        if(jmax<0) {printf("\n You are degenerate!!!\n"); return;}
        tmp=zerov[k];
        zerov[k]=basis[jmax];
        basis[jmax]=tmp;
        for(j=0;j<LNN1;j++)
            R[j][LGG1+1]=0;
        R[jmax][LGG1+1]=1;
        Rcost[LGG1+1]=0.;
        for(j=0;j<LNN1;j++)
         if(j!=jmax)
           { if(R[j][k]==1) for(i=0;i<LN1;i++) R[j][i]-=R[jmax][i];
             else if(R[j][k]==-1) for(i=0;i<LN1;i++) R[j][i]+=R[jmax][i]; 
           }
       pivot=Rcost[k];
       for(i=0;i<LN1;i++) Rcost[i]-=pivot*R[jmax][i];
        for(i=0;i<LNN1;i++) R[i][k]=R[i][LGG1+1];
        Rcost[k]=Rcost[LGG1+1];
       
      }
   }
   

 /* Reconstruct the function                                               */
  for(i=0;i<lg;i++) fp[i]=fpo[i];   

  for(i=0;i<LGG1;i++)
     if (zerov[i]<LGG1) fp[coresp[zerov[i]]]=0;
  for(i=0;i<LNN1;i++)
     if(basis[i]<LGG1) fp[coresp[basis[i]]]=R[i][LGG1];

  /* Evaluate J                                                            */
  printf("\n LP: LNN1=%d  LGG1=%d  J=%.5lf \n",LNN1,LGG1,evalJc(fp));


 /* Test if the  filter is stack                                           */
  best=1;
 for(i=lg-1;i>=0;i--)
   { fpocket[i]=0; 
     if(fp[i])
     {  fpocket[i]=1;
       depl=1;
       for(j=0;j<winsize;j++)
         { j2=(depl | i); 
           depl<<=1;
           if(j2>i) {if(!fp[j2]){ fpocket[i]=0;best=0;break;} } 
                 /*j2 lays under i*/
         }
     }
   }
  
if(best==0) 
    printf(" Stack is not a stack!\n");
  else
    printf(" Stack is stack!\n");


  

  printf("\n LNN1=%d  LGG1=%d  J=%.5lf \n",LNN1,LGG1,evalJc(fp));

  
  printf("  %.5lf",evalJc(fp));
  printf("   %.3lf",((double)clock()-cl1)/CLOCKS_PER_SEC);   

return(0);

}

#endif /* #ifdef RUN_LP */