// Copyright (C) Jaravine VA, Orekhov VYu (2006), Swedish NMR Centre, Gothenburg.
// 

#define _MAIN_
#define VERS "2006-08-08"



/* Compiler defs below:  these or nothing */
// #define INDEX3ONLY  
/*  #define SGI   */



/* history 
    2006-05-23 spheader.c
   	      + flow of indices of 4 types; 
	      + if the indices are read from the file, additional points will be selected if requested number is more
	      + add minimal indices ;
	      
    2005-09-06 spheader.c
   	      + change flow of indices from N-points package random into a continous 1 point random flow 
	        continuous point selection 1 by 1 (no tries)
		but table is still filled, and table outputs indices with the option 'regular'
	      + add flag to override checking for point repetition - usefull for high sparse match	
   	      + checking/repairing holes before point selection
	      + option for shifting by any number in any dimension

   2005-05-20 spheader4D.c 
              + option for projection reconstruction for Eriks Kupce  PR= ;
   	      + option for reading indices from file with indices and numbers in any format - the flag read. 

   2005-05-05 spheader4D.c  
              + 4D case ;+ modulation due to homonuclear coupling (for any number of indirect dims) ;
	      prepared for usage in vnmr ; 

   2005-04-04 spheader3D.c 
              + 3D case; + mirror image for 0 and 1/2 point cases (for any indirect dim) ;
	    
   2004-12-12 putsparseC2_1D.c - generate sparse (matched acq) header to input dense mdd data file ; 
              + 1D case; several data and functions are changed ; reads mdd file for size pars . 
	      no reading of fids, pars or other data - just printing header files ;
	      shape processing is separated into file 'pros.c' 

   2004-09-09 resdatC2_1D.c - collapse ND shapes; input data are mdd format shapes with amplitudes.
   
*/

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


#define MAXMODE       5
#define MAXNDSIZE    16
#define MAXFNAME    256
#define MAXLINE     512

#define YES           1 
#define NO         (-1)
#define OK            1
#define SPARSEMINIM    0
#define MAXFILLHOLES  (MAXNDSIZE/2)
#define PI            M_PI
#define MAXLISTWORDS  32
#define MAXLISTWORD   128
#define MAXCYCLESCOUNT 100000000

 char comment[MAXLINE];
 char buff[MAXLINE];
 char buff0[MAXLINE];
 char buff1[MAXLINE];
 char buff2[MAXLINE];
 char allargv[MAXLINE];
 char sparseType[MAXLINE];
 char **words;

typedef struct SPARSE { 
    float  inputfraction, resfraction;
    unsigned char *flg, *flgread;     // flg=1 if the point is ON in sparse; flgread=1 the point is read from file
    float  *r,*i;                     // optional data read for 1Ds
    unsigned long **countproj;        // square table NDxdsize[j] of counters of all points projections onto all axes
    unsigned long *divisor,*leftover; // used for translating indices to 0,1 code  
    unsigned int *B,index;            // for current point: Binary array of translated index  (instead of *D[i] ) e.g. 110001=49
    int *indices,*indicesrandom;      // all indices for the dim in random order as they arrived  

    int NPMIN; 		  // min   incremental [evolution] indices span , e.g. from 0 to 64
    int NPMAX; 		  // max  
    
    int NSP,NSP_;         // number of sparse points initially requested (with _ determined after sparsing)
    int ndsize;		  // number of split dimensions
    int dsize[MAXNDSIZE]; // split dimensions for all axes
    int NPMULT;       // product of the split dimensions 

    float expDe, SW, T2, split,J; int ishift;    // exponential decay, spectral width, T2, splitting constant; index shift;
    char  sOption[MAXLISTWORD];		         // string of sparse options
    unsigned char regular,shuffle,over,subset;   // regular,shuffle,over:  index flow 0 or 1
} S_SPARSE;


int       prepareSparseHeader (int argc, char **argv);
int 	  makeSparseMatchedArray(struct SPARSE **sp, int m, int mode);
int       checkSparseHoles      (struct SPARSE **sp, int m, int mode);
int       repairSparseHoles     (struct SPARSE **sp, int m, int mode);
int       printSparseHeader     (struct SPARSE **sp,  int mode, FILE *ofp);
int       printSparsePointIndex (struct SPARSE **sp,  int mode, FILE *ofp, FILE *ofpi);
int       printSparseExp        (struct SPARSE **sp,  int mode, FILE *ofp);
int       printSparseData       (struct SPARSE **sp,  int mode, FILE *ofp);
int       allocArraysinSparse  (struct SPARSE *sp );
unsigned int splitindices2index(struct SPARSE *sp);
int       index2splitindices(struct SPARSE *sp, int i);
int       printsplitindex(int p, int mirep, int m, struct SPARSE **sp, int im, FILE *ofp);
int       SortIntArray(int size, int *a, int *b, unsigned char flag);

char      *read_N_strings_from_string (int N, char *string, char **array);
char      *read_N_inumbers_from_string(int N, char *string, int   *array);
char      *read_N_fnumbers_from_string(int N, char *string, float *array);
char      *get_time_string();
long      get_time_sec(char *strtimevalue, char *format);
char      *strstrcase(const char *s1, const char *s2);
FILE      *Fopen(char *fname, char *fm);



#ifdef _MAIN_
main (int argc, char **argv)
{  
   if(0!= prepareSparseHeader (argc, argv) ) exit(-1);
   exit(0);
}
#endif


int prepareSparseHeader (int argc, char **argv)
{ int    m, i, j, ni, kk, d, seed, toss, mode;    
  char *s,*s1;
  struct SPARSE *sp[MAXMODE];
  char   fout_0[MAXFNAME],fout_1[MAXFNAME],fout_2[MAXFNAME],fout_3[MAXFNAME],fout_4[MAXFNAME],fout_5[MAXFNAME]; 		
  FILE   *f1,*f2,*f3,*f4,*f5;

  char help[]={" option order: \n\
     FILEOUT  MODE  NPMIN  NPMAX NPSPARSE  NDSIZES  DSIZES  SW  T2  SPARSETYPE  SPARSEOPTS \n all arguments should be enclosed in '' \n\n\
   1  FILEOUT  'nameCNH' 				# output root name        \n\
   2  MODE     '3 [reg/shuffle/over/subset] '		# n dimensions, optional option for index random selection (default: shuffle)  \n\n\
                reg	- regular =uniform sampling, must set T2 to high values e.g. T2=99 and J=0  \n\
		shuffle - incremental matched sampling, points are checked for no repetition   \n\
		over	- overriding the repetition check; so sampling can be repeated  \n\
		subset  - set this option if only reading from file; if this option is not set reading from file is done before other options are invoked \n\n\
   3  NPMIN    '0 0 0'  		 	        # min indices [for evolution increments] \n\
   4  NPMAX    '30 60 100'			        # max indices \n\
   5  NSP      '8  16 100' 	  			# n of points in sparse data	\n\
   6  NDSIZES  '5  6  1'   	  			# ndsize	\n\
   7  DSIZES   '2 2 2 2 2   2 2 2 2 2 2  100'      	# dsize[ndsize]  (multiples can be higher than NPMAX) \n\
   8  SW       '6770.087 2500 12001.2' 			# spectral widths [Hz] \n\
   9  T2       '0.050 1 999' 		 		# T2 [sec]  use e.g. 1 or 999 for CT (=almost constant)	\n\
   10 SPARSETYPE 'seed=12345' or 'seed=time' or 'PR=0:1,2:3,...'  or/and 'file=name [format]'   \n\n\
                 # point selection : random seed, from file, PR (projections); \n\
		 # for initialization can read indices from file (if 1st field is not 'int' must define format) 2nd field is Re of data point, 3rd field Im if complex ; \n\n\
   11 SPARSEOPTS 'J=35,ishift=1  CT,f180=n  0'	# options for nonrandom matching to coherence evolution \n\n\
  "};


  char examples[]={" \
spheader tmp '1 reg' '96' '64' '7' '2 2 2  2 2 2 2' '2000' '0.020' 'seed=12305' '0' \n\
spheader tmp '1 reg' '60' '23' '6' '2 2 2 2 2 2' '2000' '0.20' 'seed=12305 file=tmp   '0' \n\
spheader tmp 2 '96 58' '64 58' '7 1' '2 2 2 2 2 2 2  58' '2000 0' '0.020 0' 'seed=1235' '0' \n\
spheader tmp 2 '96 58' '16 58' '7 1' '2 2 2 2 2 2 2  58' '2000 0' '0.020 0' 'seed=1234' '0' \n\
spheader tmpCNH  3 '30  60 32' '15  30 32' '5 6 1' '2 2 2 2 2	2 2 2 2 2 2  32' '3770.087 2500 12001.2' '0.029 1 999' 'seed=12305' '0 CT,f180=n 0'   \n\
spheader tmpCNH  3 '400 60 32' '20 20 32' '9 6 1' '2 2 2 2 2 2 2 2 2  2 2 2 2 2 2  32' '6770.087 2500 12001.2' '0.040 1 999' 'seed=12305' 'J=35.0 CT,f180=n 0'	\n\
spheader tmpHCCH 4 '120 60 60 4' '30 20 20 4' '7 6 6 1' '2 2 2 2 2 2 2   2 2 2 2 2 2  2 2 2 2 2 2  4' '12001 2500 2500 12001' '0.08 0.04 0.04 1' 'seed=12305' '0 f180=n CT,f180=n 0'  \n\
spheader tmpHC   2 '624 1' '624 1' '10 1' '2 2 2 2 2 2 2 2 2 2 1' '12001 2500' '0.08 1' 'file=name' '0 0'  \n\
spheader tmpHCCH 4 '8 4 4 4' '2 3 3 4' '3 2 2 1' '2 2 2   2 2  2 2   4' '12001 2500 2500 12001' '0.08 0.04 0.04 1' 'seed=12305' '0 J=35.0,f180=n  J=35.0,CT,f180=n 0' \n\
spheader tmpCNH  3 '60 30 4' '50 25 4' '6 5 1' '2 2 2 2 2 2   2 2 2 2 2  4' '6770.087 2500 12001.2' '0.020 1 999' 'seed=12305' 'J=35,ishift=1  CT,f180=n  0'  \n\
spheader tmpCNH '3 over' '64 32 4' '50 20 4' '6 5 1' '2 2 2 2 2 2   2 2 2 2 2  4' '6770.087 2500 12001.2' '0.010 1 999' 'seed=12305' 'J=35 CT,f180=n  0' \n\
spheader tmpCNH '3 shuffle' '0 0 0' '64 32 1' '10 32 1' '6 5 1' '2 2 2 2 2 2   2 2 2 2 2  1' '1499.99 1500 12001.2' '1 1 1' 'seed=54321' 'CT,f180=n CT,f180=n 0' \n\
   "};
  

/* printf("\nsizes: %d %d %d  %d %d %d \n", sizeof(char),  sizeof(short), sizeof(int),sizeof(long),sizeof(float),sizeof(double)); 
  */
 
   

 if(argc<11) 
 { 
    printf("\n%s -- make sparse headers for mdd 1-4D sparse data; \n",argv[0]);
    printf("\nUsage:\n%s <fileout> <option2> <option3> ... <optionN>  \n",argv[0]); 
    printf("\n%s\n",help /*,examples */);
    printf("output, the 5 headers are:\n");
    printf(".hdr_1,2 - headers for mdd data files; \n");
    printf(".hdr_3   - decimal indices for vnmr;\n");
    printf(".hdr_4   -  illustration of selection table in ascii text format and dimensions exponent profiles;\n");
    printf("optional .hdr_5 - is printed only if measurements are read into this program (for future developments) \n");
   
    return(-1);
 }

  allargv[0]=0;
  for(i=0;i<argc;i++) { sprintf(buff," '%s' ",argv[i]);    strncat(allargv,buff,MAXLINE); } 

  ni=sscanf(argv[2],"%d",&mode); if(ni<1) return(-1);  // read mode = 1D - 5D

  // allocation
  for(m=0;m<=mode; m++)
  {   if(NULL==(sp[m]=calloc(1,sizeof(struct SPARSE)) )  )  {fprintf(stderr,"Cannot allocate memory\n"); return(-1);}
      strcpy(sp[m]->sOption, "0");
      sp[m]->expDe=999.0;
      sp[m]->split=0.0;
      sp[m]->regular=0; sp[m]->shuffle=1; sp[m]->over=0; sp[m]->subset=0;
      
  }
  words=(char**)calloc(MAXLISTWORDS,sizeof(char*));  if(NULL==words) {fprintf(stderr,"Error: main: Cannot allocate memory\n"); exit(1);}

  if( strstrcase( argv[2], "reg"     )!=NULL) sp[0]->regular=1,sp[0]->shuffle=0,sp[0]->over=0;
  if( strstrcase( argv[2], "shuffle" )!=NULL) sp[0]->regular=0,sp[0]->shuffle=1,sp[0]->over=0;
  if( strstrcase( argv[2], "over"    )!=NULL) sp[0]->regular=0,sp[0]->shuffle=1,sp[0]->over=1;
  if( strstrcase( argv[2], "subset"  )!=NULL) sp[0]->regular=0,sp[0]->shuffle=1,sp[0]->over=0,sp[0]->subset=1;
  
  fprintf(stdout,"regular=%d shuffle=%d over=%d subset=%d \n", sp[0]->regular,sp[0]->shuffle,sp[0]->over,sp[0]->subset);
  
  for(s=argv[3],m=1;m<=mode; m++)
  {   s= read_N_inumbers_from_string(1, s, &(sp[m]->NPMIN) );    // read number of points	 for each mode dimension
  }
  for(s=argv[4],m=1;m<=mode; m++)
  {   s= read_N_inumbers_from_string(1, s, &(sp[m]->NPMAX) );    // read number of points	 for each mode dimension
  }
  for(s=argv[5],m=1;m<=mode; m++)
  {   s= read_N_inumbers_from_string(1, s, &(sp[m]->NSP) );  // read number of desired sparse points	 for each mode dimension
  }
  for(s=argv[6],m=1;m<=mode; m++)
  {   s= read_N_inumbers_from_string(1, s, &(sp[m]->ndsize) );   // read number of split dimensions for each mode dimension
  }
  for(s=argv[7],m=1;m<=mode; m++)
  {   s= read_N_inumbers_from_string(sp[m]->ndsize, s, sp[m]->dsize );   // read split dimensions sizes for each mode dimension
      for(sp[m]->NPMULT=1, d=0; d < sp[m]->ndsize;d++) sp[m]->NPMULT*= sp[m]->dsize[d];   //calc product of split dimensions 
      sp[m]->inputfraction=(float)(sp[m]->NSP)/sp[m]->NPMULT;
  }
  for(s=argv[8],m=1;m<=mode; m++)
  {   s= read_N_fnumbers_from_string(1, s, &(sp[m]->SW) );             
  } 
  for(s=argv[9],m=1;m<=mode; m++)
  {   s= read_N_fnumbers_from_string(1, s, &(sp[m]->T2) );
      sp[m]->expDe= 1.0/(sp[m]->T2*sp[m]->SW);
  }
  
  if(argc>=11)
  {   strncpy(sparseType,argv[10],MAXLINE);
      seed=5;
      if( (s1=strstrcase( sparseType, "seed=" ))!=NULL) 
      {    if(strncmp("time",s1+5,4)==0) {seed=(unsigned int)time(NULL);} else { sscanf(s1+5,"%d", &seed); }
      } 
      if(isdigit(sparseType[0])) sscanf(sparseType,"%d", &seed);   printf("SEED=%d\n", seed);  
      srand(seed);
  }   

  if(argc>=12)
  {  s=argv[11];
     strncpy(sp[0]->sOption, s, MAXLISTWORD);
     read_N_strings_from_string(mode, s, words);
      
     for(m=1;m<=mode; m++)
     {  
        strncpy(sp[m]->sOption, words[m-1], MAXLISTWORD);

        sp[m]->J=0;sp[m]->split=0;
        if( (s1=strstrcase( sp[m]->sOption, "J=" )) != NULL)
        {    read_N_fnumbers_from_string(1, s1+2, &sp[m]->J );   	   
             sp[m]->split = (PI*sp[m]->J)/sp[m]->SW; 
        }
	sp[m]->ishift=0;
        if( (s1=strstrcase( sp[m]->sOption, "ishift=" )) != NULL)
        {    read_N_inumbers_from_string(1, s1+7, &sp[m]->ishift );
        }
     }	
  } 
    
 if(mode==3)
  {    sp[0]->NPMAX= sp[1]->NPMAX* sp[2]->NPMAX;
       sp[0]->NSP = sp[1]->NSP * sp[2]->NSP;
       sp[0]->ndsize=2;
       sp[0]->dsize[0]= sp[1]->NPMULT; sp[0]->dsize[1]= sp[2]->NPMULT;
       sp[0]->NPMULT =  sp[1]->NPMULT * sp[2]->NPMULT;
       sp[0]->inputfraction=(float)(sp[0]->NSP)/sp[0]->NPMAX;
       sp[0]->expDe=  0.01;sp[0]->split= 0.01;
  }
 if(mode==4)
  {    sp[0]->NPMAX= sp[1]->NPMAX* sp[2]->NPMAX* sp[3]->NPMAX;
       sp[0]->NSP = sp[1]->NSP * sp[2]->NSP * sp[3]->NSP ;
       sp[0]->ndsize=2;
       sp[0]->dsize[0]= sp[1]->NPMULT; sp[0]->dsize[1]= sp[2]->NPMULT; sp[0]->dsize[2]= sp[3]->NPMULT;
       sp[0]->NPMULT = sp[1]->NPMULT * sp[2]->NPMULT * sp[3]->NPMULT;
       sp[0]->inputfraction=(float)(sp[0]->NSP)/sp[0]->NPMAX;
       sp[0]->expDe=  0.01;sp[0]->split= 0.01;
  }
  
  
  for(m=0;m<=mode; m++) if(0!=allocArraysinSparse(sp[m])) return(-1);

  switch(mode)
  { case 1: makeSparseMatchedArray(sp,1,mode);
    break;
    case 2: makeSparseMatchedArray(sp,1,mode); makeSparseMatchedArray(sp,2,mode);
    break;
    case 3: makeSparseMatchedArray(sp,0,mode); makeSparseMatchedArray(sp,3,mode); 
    break;
    case 4: makeSparseMatchedArray(sp,0,mode); makeSparseMatchedArray(sp,4,mode); 
    break; 
     
  }


  // open output files
  strcpy(fout_0,argv[1]); 
#ifdef INDEX3ONLY
  sprintf(fout_2,"%s.hdr_2",fout_0);  f2=Fopen(fout_2,"w");
  sprintf(fout_3,"%s.hdr_3",fout_0);  f3=Fopen(fout_3,"w");
#else
  sprintf(fout_1,"%s.hdr_1",fout_0);  f1=Fopen(fout_1,"w");
  sprintf(fout_2,"%s.hdr_2",fout_0);  f2=Fopen(fout_2,"w");
  sprintf(fout_3,"%s.hdr_3",fout_0);  f3=Fopen(fout_3,"w");
  sprintf(fout_4,"%s.hdr_4",fout_0);  f4=Fopen(fout_4,"w");
  sprintf(fout_5,"%s.hdr_5",fout_0);  f5=Fopen(fout_5,"w");
#endif


    
#ifdef INDEX3ONLY
   printSparseIndices(sp,mode,f2,f3);
#else
   printSparseHeader(sp,mode,f1);
   printSparseIndices(sp,mode,f2,f3);
   printSparseExp(sp,mode,f4);
   printSparseData(sp,mode,f5);
#endif
   
 
 fclose(f3); 
 
 return(0);
}


// alloc & calc divisor
int allocArraysinSparse(struct SPARSE *spm)
{  int dim,d,i;  long mult;
   
  // allocate memory
  spm->countproj=(unsigned long**)calloc(spm->ndsize,sizeof(unsigned long));
  for(dim=0; dim<spm->ndsize; dim++) spm->countproj[dim]=(unsigned long*)calloc(spm->dsize[dim],sizeof(unsigned long));

  spm->B=(unsigned int *)calloc(spm->ndsize,sizeof(unsigned int));   
  spm->divisor =(unsigned long *)calloc(spm->ndsize,sizeof(unsigned long)); 
  spm->leftover=(unsigned long *)calloc(spm->ndsize,sizeof(unsigned long));   
  
  //calc divisor   -- fast index first, so start filling last 
  for(dim=spm->ndsize-1; dim>=0; dim--)
  {   for(mult=1, d=0; d < dim;d++) mult*= spm->dsize[d]; 
      spm->divisor[dim]=mult;
  }    

   spm->flg=(unsigned char *)calloc(spm->NPMULT,sizeof(unsigned char));     if(NULL==spm->flg) { fprintf(stderr,"Cannot allocate memory \n"); return(-1); }
   spm->flgread=(unsigned char *)calloc(spm->NPMULT,sizeof(unsigned char)); if(NULL==spm->flgread) { fprintf(stderr,"Cannot allocate memory \n"); return(-1); }
   
  // init sparse points to 0
  for(i=0; i<spm->NPMULT; i++) spm->flg[i]=spm->flgread[i]=0; 
  spm->NSP_=0;
  spm->indices=(int *)calloc(spm->NPMULT,sizeof(int));   	for(i=0;i<spm->NPMULT;i++) spm->indices[i]=0;
  spm->indicesrandom=(int *)calloc(spm->NPMULT,sizeof(int));    for(i=0;i<spm->NPMULT;i++) spm->indicesrandom[i]=0;

  spm->r=(float*)calloc(spm->NPMULT,sizeof(float));   for(i=0;i<spm->NPMULT;i++) spm->r[i]=0;
  spm->i=(float*)calloc(spm->NPMULT,sizeof(float));   for(i=0;i<spm->NPMULT;i++) spm->i[i]=0;
 
  return(0);
}



int makeSparseMatchedArray(struct SPARSE **sp, int m, int mode)
{  int i,ic,  i1, i2, i3, ileft, ip, ncheck; 
   float r1,e1, e1max;
   struct SPARSE *spm;
   char *s1,*s2,*s,field[64];
   int ni,ns,prs[64],hit;
   
   spm=sp[m]; 
   printf("\n\n makeSparseMatchedArray m=%d mode=%d\n",m,mode);
   
   
   printf("\nsparse options:  '%s' \n", sparseType );
   for(i=0;i<=mode;i++)  printf(" m=%d option='%s' decay=%.6f split=%.3f(J=%.1f) \n", i, sp[i]->sOption, sp[i]->expDe, sp[i]->split, sp[i]->J);  
   
   
   e1max=0.0;

   //
   // case of indices supplied in the file=name [format=]'     if file indices are not int must define input format , output will be int 
   // 
   if( (s1=strstr( sparseType, "file=" ))!=NULL && ((m==0 && mode!=1) || (m==1 && mode==1) )  ) 
   {    char  fsparsename[MAXFNAME];  FILE *fpsparse;  
        static char indexformat[128]={"%d %d"};
        float re,im; 
        sscanf(s1+5,"%s", fsparsename);   fpsparse=Fopen(fsparsename,"r");

 	if( (s2=strstr( sparseType,"format="))!=NULL) sscanf(s2+7,"%s", indexformat);
        printf("filename: '%s' input format: '%s'  mode=%d  fpsparse=%ld \n",fsparsename, indexformat, mode , (long)fpsparse);
        
	switch(mode)
	{ case 1: 
	  case 2:  
	        while( fgets(buff,MAXLINE,fpsparse) != NULL ) 
		{    for(i=0;i<strlen(buff);i++) if(buff[i]==',') buff[i]=' ';
 		     ni=sscanf(buff,"%s %s %s",buff0,buff1,buff2);
 		     
		     if(strlen(buff)<=1 || isalpha(buff[0]) || buff[0]=='#' ) continue;
		     ni=3; // how many fields to read
		     switch(ni) 
		     { case 1:  sscanf(buff,indexformat,&i );  
		          		           
		        break;
		       case 2:  sscanf(buff0,indexformat,&i);
		                sscanf(buff1,"%g",&re); 
				spm->r[i]=re; spm->i[i]=0;
		                spm->flg[i]=spm->flgread[i]=1; spm->indicesrandom[spm->NSP_]=i; spm->NSP_++;
			
		        break;
		       case 3:  sscanf(buff0,indexformat,&i);
		                sscanf(buff1,"%g",&re); 
				sscanf(buff2,"%g",&im); 
		                spm->r[i]=re; spm->i[i]=im;
		                spm->flg[i]=spm->flgread[i]=1; spm->indicesrandom[spm->NSP_]=i; spm->NSP_++;
		        break;
		     }
		     
		}
          break;
	  case 3:  
	        while( fgets(buff,MAXLINE,fpsparse) != NULL ) 
		{    
//		     printf("#  '%s'  \n", buff);
		     
		     for(i=0;i<strlen(buff);i++) if(buff[i]==',') buff[i]=' ';
 		     ns=sscanf(buff,indexformat,&i1, &i2); 
		     if( ns!=2 ) { printf("wrong data in line '%s' \n", buff); continue; } 
		     i= i2* sp[1]->NPMAX + i1;
                     spm->flg[i]=spm->flgread[i]=1; spm->indicesrandom[spm->NSP_]=i; spm->NSP_++;
		}
	
          break;
	  
        }

	fclose(fpsparse); 
        spm->resfraction=(float)(spm->NSP_)/spm->NPMAX;
        printf("\nFILE: *** sparse=%.5f (has read %d points from %d ) \n", spm->resfraction, spm->NSP_, spm->NPMAX);
   } 

     
   //
   // special case of projection reconstruction: no holes or size checks, no random - just diagonals as requested in the list after "PR=" 0:1,1:0,1:1,1:2,2:1 
   //

   if( (s1=strstrcase( sparseType, "PR=" )) != NULL && m==0)
   {    // presently any diagonals denoted with integers : 
        // scan pr option 
	for(i=0,ic=0,s=s1+3;   i<64 && ic<strlen(sparseType);   i++, ic+=strlen(field) ) 
  	{   while( ( !isdigit((int)s[ic])  ) && (ic<strlen(sparseType)) ) ic++;
            ni=sscanf(&s[ic],"%d",&prs[i]);  if(ni!=1) break;
	    sprintf(field,"%d",prs[i]); 
  	}
	ni=i;         
        for(i=0; i<spm->NPMAX && ni>0; i++)
        {    spm->flg[i]=0;
	     switch(mode)
	     { case 1: 
	       case 2: 
	       break;
	       case 3:i2= i/sp[1]->NPMAX; i1= i%sp[1]->NPMAX; 
	              for(hit=0,ic=0;ic<ni;ic+=2) 
		      {   if(i2*prs[ic]==i1*prs[ic+1]) {hit=1;break;}
		      }
 	  	      if(hit==1) 
		      { spm->flg[i]=1; spm->indicesrandom[spm->NSP_]=i; spm->NSP_++; 
		      }  
	  	      sp[2]->NSP=1;
	       break;
	       case 4:i3= i/(sp[1]->NPMAX*sp[2]->NPMAX); ileft=i%(sp[1]->NPMAX*sp[2]->NPMAX); i2= ileft/sp[1]->NPMAX; i1= ileft%sp[1]->NPMAX;  
	  	      if(i3==0 || i2==0 || i1==0 || i2==i1 || i3==i2 || i1==i3) 
		      { spm->flg[i]=1; spm->indicesrandom[spm->NSP_]=i; spm->NSP_++; 
		      }    
	  	      sp[2]->NSP=sp[3]->NSP=1;
	       break;
	     }  
        }
        spm->resfraction=(float)(spm->NSP_)/spm->NPMAX;
        printf("%s\nPR: *** sparse=%.5f (%d points from %d ) \n", sparseType, spm->resfraction, spm->NSP_, spm->NPMAX);
	sp[1]->NSP=spm->NSP_;
   }

   
   
   //
   // general case of matched acquisition: with decays, J-homonuclear modulations       
   //

// #  if( (s1=strstrcase( sparseType, "seed=" ))!=NULL )
   {    int i_;
        long count;
        printf("\n\ndim=%d \n ***SEED: Requested  sparse=%.5f (%d points from %d mult=%d) \n", m, spm->inputfraction, spm->NSP, spm->NPMAX, spm->NPMULT);
   
        // here file selection is initially ignored
	if( spm->subset==1 )
	{     for(i=0; i<spm->NPMULT;i++) spm->flg[i]=0;
              spm->NSP_=0; 
	}
	
     	if( spm->ndsize >=2 )	// check & correct holes
      	{    for(i=0; i<MAXFILLHOLES; i++) 
	     {    ncheck=checkSparseHoles(sp,m,mode);  if(ncheck==0) break;
	          if(NO==repairSparseHoles(sp,m,mode) ) break;
	     }
      	}

        // calc exp matched acquisition 
	i_=0; count=0;
        while( spm->NSP_<spm->NSP && count<MAXCYCLESCOUNT)
        {  
#ifdef SGI
             i=((spm->NPMAX-1)*(float)rand()/RAND_MAX); 
#else
             i=((spm->NPMAX  )*(float)rand()/RAND_MAX); 
#endif
             count++;
// printf("SEED: newi=%d  spm->NSP_=%d  spm->NSP=%d spm->NPMAX=%d \n",i, spm->NSP_, spm->NSP, spm->NPMAX);	
  	     if( spm->flg[i]!=0 && sp[0]->over==0 ) { continue;}
	     
  	     /* if( spm->flgread[i]==0 && ((s1=strstr( sparseType, "file=" ))!=NULL) ) {continue;} */ 
	     if( spm->flgread[i]==0 && ((s1=strstr( sparseType, "file=" ))!=NULL) && spm->subset==1 ) {continue;}
	     /* only a subset of read schedule is allowed */
	     
	     switch(mode)
	     { case 1: 
	       case 2:e1=fabs(cos(sp[m]->split*i))*exp( -sp[m]->expDe*i); i1=i; i2=0;
	       break;
	       case 3:i2= i/sp[1]->NPMAX; i1= i%sp[1]->NPMAX;    // index spans input NPMAX ranges, fast dim is i1 or sp[1]
	  	      if(m==0)  e1= fabs(cos(sp[1]->split*i1))*exp(-sp[1]->expDe*i1) * fabs(cos(sp[2]->split*i2))*exp(-sp[2]->expDe*i2); 
	  	      else e1=fabs(cos(sp[m]->split*i))*exp(-sp[m]->expDe*i); 
	       break;
	       case 4:i3= i/(sp[1]->NPMAX*sp[2]->NPMAX); ileft=i%(sp[1]->NPMAX*sp[2]->NPMAX); i2= ileft/sp[1]->NPMAX; i1= ileft%sp[1]->NPMAX;   
	  	      if(m==0)  e1= fabs(cos(sp[1]->split*i1))*exp(-sp[1]->expDe*i1) * fabs(cos(sp[2]->split*i2))*exp(-sp[2]->expDe*i2) * fabs(cos(sp[3]->split*i3))*exp(-sp[3]->expDe*i3); 
	  	      e1=fabs(cos(sp[m]->split*i))*exp( -sp[m]->expDe*i);
	       break;
	       
	     }

	     /* check min of the indices */
	     // printf("- i=%d  i1=%d i2=%d    npmin1=%d npmin2=%d   [%d]\n", i, i1, i2, sp[1]->NPMIN, sp[2]->NPMIN, (int)count  );
	     if( m==0 && ( (mode<=2 && i1<spm->NPMIN) || (mode==3 && (i1<sp[1]->NPMIN || i2<sp[2]->NPMIN)) || (mode==4 && (i1<sp[1]->NPMIN || i2<sp[2]->NPMIN || i3<sp[3]->NPMIN))  ) )
             continue;

             r1=(float)rand()/RAND_MAX;	     if(e1>=e1max) e1max=e1;
	     if(i_>10000 && e1max<0.001) e1=e1*1000;  
             if(e1>=r1 || i<=SPARSEMINIM) 
	     {   spm->flg[i]++; spm->indicesrandom[spm->NSP_]=i; spm->NSP_++;
	     
	         spm->resfraction=(float)(spm->NSP_)/spm->NPMAX;
/*		 
                 printf("\n  %d   seed_sparse= %.5f  ", spm->NSP_,  spm->resfraction);
                 printf(" index= %d (%d,%d)   %.4f   ", i, i1,i2, e1);
*/
             } else {}
        } 		
        spm->resfraction=(float)(spm->NSP_)/spm->NPMAX;
        printf(" ***SEED: sparse request=%.5f  final=%.5f  %d points  selection from 1 to %d  mult=%d \n", spm->inputfraction, spm->resfraction,  spm->NSP_, spm->NPMAX, spm->NPMULT );

    }  
   

   if( spm->ndsize >=2 )  ncheck=checkSparseHoles(sp,m,mode);
   SortIntArray( spm->NSP_, spm->indicesrandom, spm->indices, sp[0]->regular);

 return(0);
}




int checkSparseHoles(struct SPARSE **sp, int m, int mode)
{   int  i,d,im, k,i1,i2,i3, ileft, foundholes;
    struct SPARSE *spm;
    
    spm=sp[m];
    foundholes=0;
    for(im=0;im<=mode;im++) for(d=0; d<sp[im]->ndsize;d++) for(k=0; k<sp[im]->dsize[d];k++) sp[im]->countproj[d][k]=0;  // init countproj

    for(i=0; i<spm->NPMAX;i++)
    {   if( spm->flg[i]==0 ) continue;
	switch(mode)
	{ case 1: 
	  case 2:index2splitindices(spm,i);  for(d=0; d < spm->ndsize;d++) (spm->countproj[d][ spm->B[d] ])++;
	  break;
	  case 3:i2= i/sp[1]->NPMAX; i1= i%sp[1]->NPMAX;
		 index2splitindices(sp[1],i1);    for(d=0; d < sp[1]->ndsize;d++) (sp[1]->countproj[d][ sp[1]->B[d] ])++;
		 index2splitindices(sp[2],i2);    for(d=0; d < sp[2]->ndsize;d++) (sp[2]->countproj[d][ sp[2]->B[d] ])++;
	  break;
	  case 4:i3= i/(sp[1]->NPMAX*sp[2]->NPMAX); ileft=i%(sp[1]->NPMAX*sp[2]->NPMAX); i2= ileft/sp[1]->NPMAX; i1= ileft%sp[1]->NPMAX;   
		 index2splitindices(sp[1],i1);    for(d=0; d < sp[1]->ndsize;d++) (sp[1]->countproj[d][ sp[1]->B[d] ])++;
		 index2splitindices(sp[2],i2);    for(d=0; d < sp[2]->ndsize;d++) (sp[2]->countproj[d][ sp[2]->B[d] ])++;
		 index2splitindices(sp[3],i3);    for(d=0; d < sp[3]->ndsize;d++) (sp[3]->countproj[d][ sp[3]->B[d] ])++;
	  break;
	}  
    }

    printf("   Counting projections from %d points: \n", spm->NSP_);
    switch(mode)
    { case 1: 
      case 2: printf(" for %d dims %d : ",m,spm->ndsize);
              for(d=0; d < spm->ndsize;d++) 
	      {   for(k=0; k<spm->dsize[d];k++) { if(spm->countproj[d][k]==0) foundholes++;  printf(" %ld",spm->countproj[d][k]); };  printf(", ");
	      }
	      printf("\n");
      break;
      case 3: printf(" for 1 dims %d : ",sp[1]->ndsize);
              for(d=0; d < sp[1]->ndsize;d++) 
    	      {   for(k=0; k<sp[1]->dsize[d];k++) { if(sp[1]->countproj[d][k]==0) foundholes++;  printf(" %ld",sp[1]->countproj[d][k]); };  printf(", ");
	      }
	      printf("\n");
    	      printf(" for 2 dims %d : ",sp[2]->ndsize);
	      for(d=0; d < sp[2]->ndsize;d++) 
	      {   for(k=0; k<sp[2]->dsize[d];k++) { if(sp[2]->countproj[d][k]==0) foundholes++;  printf(" %ld",sp[2]->countproj[d][k]); };  printf(", ");
	      }
	      printf("\n");
      break;
      case 4: printf(" for 1 dims %d : ",sp[1]->ndsize);
              for(d=0; d < sp[1]->ndsize;d++) 
    	      {   for(k=0; k<sp[1]->dsize[d];k++) { if(sp[1]->countproj[d][k]==0) foundholes++;  printf(" %ld",sp[1]->countproj[d][k]); };  printf(", ");
	      }
	      printf("\n");
    	      printf(" for 2 dims %d : ",sp[2]->ndsize);
	      for(d=0; d < sp[2]->ndsize;d++) 
	      {   for(k=0; k<sp[2]->dsize[d];k++) { if(sp[2]->countproj[d][k]==0) foundholes++;  printf(" %ld",sp[2]->countproj[d][k]); };  printf(", ");
	      }
	      printf("\n");
    	      printf(" for 3 dims %d : ",sp[3]->ndsize);
	      for(d=0; d < sp[3]->ndsize;d++) 
	      {   for(k=0; k<sp[3]->dsize[d];k++) { if(sp[3]->countproj[d][k]==0) foundholes++;  printf(" %ld",sp[3]->countproj[d][k]); };  printf(", ");
	      }
	      printf("\n");
      break;
    }  

    printf(" Foundholes: %d\n",foundholes);
    return(foundholes);
}  
  


int repairSparseHoles(struct SPARSE **sp, int m, int mode)
{   int  i,d,k1,k2,k3,im, i1,i2,i3;
    struct SPARSE *spm;
    
    spm=sp[m];
    for(im=0;im<=mode;im++) for(d=0; d < sp[im]->ndsize;d++) sp[im]->B[d]=0;  // init all Bs
    k1=k2=i1=i2=0;
    switch(mode)
    { case 1: 
      case 2: for(d=0; d < spm->ndsize;d++) 
	      {   for(k1=0; k1<spm->dsize[d];k1++)  
	          {   if(spm->countproj[d][k1]==0) {spm->B[d]=k1;break;}  // repair d,k as the first found hole
	          }
	      }
	      i=splitindices2index(spm);
      break;
      case 3: for(d=0; d < sp[1]->ndsize;d++) 
    	      {   for(k1=0; k1<sp[1]->dsize[d];k1++)  
    	          {   if(sp[1]->countproj[d][k1]==0) {sp[1]->B[d]=k1;break;}  // repair d,k as the first found hole
	          }
	      }
	      for(d=0; d < sp[2]->ndsize;d++) 
    	      {   for(k2=0; k2<sp[2]->dsize[d];k2++)  
    	          {   if(sp[2]->countproj[d][k2]==0) {sp[2]->B[d]=k2;break;}  // repair d,k as the first found hole
	          }
	      }
	      i1= splitindices2index(sp[1]) ;i2= splitindices2index(sp[2]) ;
	      i= i2* sp[1]->NPMAX + i1;				              // index spans input NPMAX ranges, fast dim is i1
      break;
      case 4: for(d=0; d < sp[1]->ndsize;d++) 
    	      {   for(k1=0; k1<sp[1]->dsize[d];k1++)  
    	          {   if(sp[1]->countproj[d][k1]==0) {sp[1]->B[d]=k1;break;}  // repair d,k as the first found hole
	          }
	      }
	      for(d=0; d < sp[2]->ndsize;d++) 
    	      {   for(k2=0; k2<sp[2]->dsize[d];k2++)  
    	          {   if(sp[2]->countproj[d][k2]==0) {sp[2]->B[d]=k2;break;}  // repair d,k as the first found hole
	          }
	      }
	      for(d=0; d < sp[3]->ndsize;d++) 
    	      {   for(k3=0; k3<sp[3]->dsize[d];k3++)  
    	          {   if(sp[3]->countproj[d][k3]==0) {sp[3]->B[d]=k3;break;}  // repair d,k as the first found hole
	          }
	      }
	      
	      i3= splitindices2index(sp[3]); i2= splitindices2index(sp[2]); i1= splitindices2index(sp[1]);
	      i= i3*sp[1]->NPMAX*sp[2]->NPMAX + i2*sp[1]->NPMAX + i1;			// index spans input NPMAX ranges, fast dim is i1
      break;
      
    }  
    if(i>=spm->NPMAX) 
    {   printf("repairSparseHoles : point %d is outside the range  \n",i); 
        i=(spm->NPMAX-1)*(float)rand()/RAND_MAX; // take random point
    } 

    // take i 
    if( spm->flg[i]==0 && (spm->flgread[i]!=0 || strstr(sparseType,"file=")==NULL)  )
    {   spm->flg[i]++;   spm->indicesrandom[spm->NSP_]=i;  spm->NSP_++; 
        spm->resfraction=(float)(spm->NSP_)/spm->NPMAX;
        printf("  %d   repseed_sparse= %.5f   index= %d    \n", spm->NSP_, spm->resfraction, i);    
    } else 
    {   switch(mode)
        { case 3: printf("Warning: void call to repairSparseHoles :  given point i=%d  hole at dims=%d %d , points i1=%d i2=%d\n",i,k1,k2, i1,i2); break;
          case 4: printf("Warning: void call to repairSparseHoles :  given point i=%d  hole at dims=%d %d %d, points i1=%d i2=%d i3=%d\n",i,k1,k2,k3, i1,i2,i3); break;	
	}
    }

    return(OK);
}



// calc indices from i
int index2splitindices(struct SPARSE *sp, int i)
{   int  d;
  
   sp->index=i;   // fast index first
   for(d=sp->ndsize-1; d>=0; d--)
   {    // razmen i  (checked)
    	if(d==sp->ndsize-1)
   	{  sp->B[d]        = (long)i/sp->divisor[d]; 
   	   sp->leftover[d] = (long)i%sp->divisor[d]; 
    	} else
    	{  sp->B[d]        = (long)sp->leftover[d+1]/sp->divisor[d]; 
    	   sp->leftover[d] = (long)sp->leftover[d+1]%sp->divisor[d]; 
    	}							
  
 
    	if( sp->B[d]<0 || sp->B[d] >= sp->dsize[d] ) 
	{  fprintf(stderr,"Warning:dimout for %d ",i); 
  	   fprintf(stderr,"  d=%d  sp->dsize[d]=%d sp->divisor[d]=%ld   sp->B[d]=%ld \n",  d, sp->dsize[d], (long)sp->divisor[d], (long)sp->B[d] );
	}
   }
   return(0);
}

// calc i  from indices
unsigned int splitindices2index(struct SPARSE *sp)
{   int d;
    sp->index=0;
    for(d=sp->ndsize-1;d>=0;d--) sp->index+= sp->divisor[d]*sp->B[d] ; 
    return(sp->index);
}



// read  N strings from given char* 'string' and store them into provided char** array; 
// return the char pointing to immediately after the last number
char *read_N_strings_from_string(int N, char *string, char **array)
{  int i,ic,ni; char *s,field[MAXLISTWORD];

 	for(i=0,ic=0,s=&string[0];   i<N;   i++,ic += strlen(field) ) 
  	{   while(isspace(s[ic]) ) ic++;
            ni=sscanf(&s[ic],"%s",field);       if(ni!=1) break;
	    array[i]= (char*)calloc(strlen(field)+1,sizeof(char));  if(NULL==array[i]) {fprintf(stderr,"Error: read_N_strings_from_string: Cannot allocate memory\n"); exit(1);}
	    strncpy(array[i], field, MAXLISTWORD); 
	    // printf(" [ %s ] ", array[i]);
  	} 
  	if(i<N) { fprintf(stderr,"number of items '%d' in the string '%s' is less than required '%d' \n",i,string,N); abort();}
	return(&s[ic]);
}

// read  N int_numbers from given char string and store them into provided int array; 
// return the char pointing to immediately after the last number
char *read_N_inumbers_from_string(int N, char *string, int *array)
{  int i,ic,ni; char *s,field[32];

 	for(i=0,ic=0,s=&string[0];   i<N;   i++,ic += strlen(field) ) 
  	{   while(!isdigit(s[ic]) ) ic++;
	    ni=sscanf(&s[ic],"%s",field);       if(ni!=1) break;
            ni=sscanf(field,"%d",&(array[i]));  if(ni!=1) break;
  	} 
  	if(i<N) { fprintf(stderr,"number of items '%d' in the string '%s' is less than required '%d' \n",i,string,N); abort(); }	
	return(&s[ic]);
}

// read  N float_numbers from given char string and store them into provided int array; 
// return the char pointing to immediately after the last number
char *read_N_fnumbers_from_string(int N, char *string, float *array)
{  int i,ic,ni; char *s,field[32];

	for(i=0,ic=0,s=&string[0];   i<N;   i++,ic += strlen(field) ) 
  	{   while(isspace(s[ic]) ) ic++;
	    ni=sscanf(&s[ic],"%s",field);       if(ni!=1) break;
            ni=sscanf(field,"%f",&(array[i]));  if(ni!=1) break;
  	} 
  	if(i<N) { fprintf(stderr,"number of items '%d' in the string '%s' is less than required '%d' \n",i,string,N); abort();}	
	return(&s[ic]);
}




char *get_time_string()
{  static char string[64], *stime;
   static time_t rawtime, t;  
   static struct tm *p;
   static int shift=1; // time shift  Berlin CET 1 ; LA PST -8 ; HK  HKG 8  ;

   time ( &rawtime );
   p = gmtime ( &rawtime );
   t = mktime(p);   
   sprintf(string,"%d/%d/%d %02d:%02d.%02d", p->tm_mday, p->tm_mon+1, p->tm_year+1900, p->tm_hour+shift, p->tm_min, p->tm_sec);
   return(string);
}


long get_time_sec(char *strtimevalue, char *format)
{  static char string[64], *stime;
   static time_t rawtime, t;  
   struct tm *p;
      //   Convert a Data-Plus-Time String to Broken-Down Time and Then into Seconds
      //   The following example demonstrates the use of strptime() to convert a string into broken-down time.
      //   The broken-down time is then converted into seconds since the Epoch using mktime().
      //   time secs since  00:00 0/0/1970
      //   strptime("6 Dec 2001 12:33:45", "%d %b %Y %H:%M:%S", &tm)
      //   strptime("6/12/2001 12:33:45", "%d/%m/%Y %H:%M:%S", &tm)
   
   time ( &rawtime );
   p = gmtime ( &rawtime );

   stime =  (char*)strptime( (char*)strtimevalue, (char*)format, (struct tm *)p ); 

   t = mktime(p);
   
   if ( t== -1 ) { printf(" format does not match input file \n"); }
   else 
   {   /* printf(" input=%s format=%s   asctime=%s  secs=%ld\n", strtimevalue, format, asctime(p), (long)t); */
   }
   return((long)t);
}




int printSparseExp(struct SPARSE **sp, int mode, FILE *ofp)
{  int  i,m; int i1,i2,i3,ileft;
   struct SPARSE *spm;

   switch(mode)
   { case 1:
     case 2:for(m=1;m<=mode; m++)
   	    {	 spm=sp[m]; spm->resfraction=(float)(spm->NSP_)/spm->NPMAX;
   	    	 printf("\n");
		 for(i=0; i<spm->NPMAX;i++) { if(spm->flg[i]!=0) printf("%d",spm->flg[i]); else printf(".");}  
  	    	 printf("\n#Result matched sparse=%.5f (%d points from %d ) \n", spm->resfraction, spm->NSP_, spm->NPMAX);

		 fprintf(ofp,"#sparseExp: m=%d   %12.5f %d %d  expDe=%.5f  mode=%d\n", m, sp[m]->resfraction, sp[m]->NSP_, sp[m]->NPMULT, sp[m]->expDe, mode);
   	    	 for(i=0; i<sp[m]->NPMAX;i++) fprintf(ofp,"%12.5f %12.5f \n",(float)(sp[m]->expDe*i),(float)sp[m]->flg[i]*exp(-sp[m]->expDe*i) );
   	    }
     break;
     case 3:spm=sp[0];
   	    spm->resfraction=(float)(spm->NSP_)/spm->NPMAX;
	    for(i=0; i<spm->NPMAX;i++) 
	    {	
	         if(i%sp[1]->NPMAX==0) printf("\n");
	         if(spm->flg[i]!=0) printf("%d",spm->flg[i]); else printf("."); 
	    
	         if(i%sp[1]->NPMAX==0) fprintf(ofp,"\n");
	         if(spm->flg[i]!=0) fprintf(ofp,"x"); else fprintf(ofp,"."); 
	    };
    	    printf("\n#Result matched sparse=%.5f (%d points from %d ) \n", spm->resfraction, spm->NSP_, spm->NPMAX, spm->expDe);
	    
   	    for(m=0;m<=mode; m++)
  	    {	  fprintf(ofp,"#sparseExp: m=%d   %12.5f %d %d  expDe=%.5f  mode=%d\n", m, sp[m]->resfraction, sp[m]->NSP_, sp[m]->NPMULT, sp[m]->expDe, mode);
   	    	  for(i=0; i<sp[m]->NPMAX;i++) fprintf(ofp,"%12.5f %12.5f \n", (float)(sp[m]->expDe*i), (float)sp[m]->flg[i]*exp(-sp[m]->expDe*i) );
	    }
     break;
     case 4:spm=sp[0];
   	    spm->resfraction=(float)(spm->NSP_)/spm->NPMAX;
	    for(i=0; i<spm->NPMAX;i++) 
	    {	 i3= i/(sp[1]->NPMAX*sp[2]->NPMAX); ileft=i%(sp[1]->NPMAX*sp[2]->NPMAX); i2= ileft/sp[1]->NPMAX; i1= ileft%sp[1]->NPMAX;
	    	    
	         if(i1==0) fprintf(ofp,"\n");
	         if(spm->flg[i]!=0) fprintf(ofp,"x"); else fprintf(ofp,"."); 
	    };
    	    printf("\n#Result matched sparse=%.5f (%d points from %d ) \n", spm->resfraction, spm->NSP_, spm->NPMAX, spm->expDe);
	    
   	    for(m=0;m<=mode; m++)
  	    {	  fprintf(ofp,"#sparseExp:  %12.5f %d %d  expDe=%.5f\n", sp[m]->resfraction, sp[m]->NSP_, sp[m]->NPMULT, sp[m]->expDe);
   	    	  for(i=0; i<sp[m]->NPMAX;i++) fprintf(ofp,"%12.5f %12.5f \n", (float)(sp[m]->expDe*i), (float)sp[m]->flg[i]*exp(-sp[m]->expDe*i) );
	    }
     break;
     
     default: break;
   } 
   return(0);
}


int printSparseData(struct SPARSE **sp, int mode, FILE *ofp)
{  int  i_,i1;

   switch(mode)
   { case 1: case 2:case 3:case 4:
	     for(i_=0; i_<sp[1]->NSP_; i_++)
	     {   i1= sp[1]->indices[i_];     //  fprintf(ofp,"i_=%d i1=%d  i1flg=%d ",i_, i1, sp[1]->flg[i1]);
	         if( sp[1]->flg[i1]==0 ) continue;
		 if( sp[1]->flgread[i1] != 0 ) fprintf(ofp,"%g %g\n", sp[1]->r[i1],sp[1]->i[i1]); else fprintf(ofp,"%d \t 0 \t 0 \n", i1);  
	     };
     break;
 
     default:  printf("\nSorry, this mode (%d core dimensions) is not implemented yet \n",mode);
   } 
          
   return(0);
}


int printSparseHeader(struct SPARSE **sp, int mode, FILE *ofp)
 { int i,m;
   int sumndsize, prodNSP; 
   
/*   
mdd sparse asc complex
spheader3D2005-15-03 (c) 2/4/2005 16:15.36
$
10 1 350
2 2 2 2 2 2 2 2 2 14 
*/   
   fprintf(ofp,"mdd sparse asc complex\n");

   sprintf(comment,"%s (c) %s ",  VERS, allargv); 
   fprintf(ofp,"%s\n$\n",comment);
 
   for(m=1,sumndsize=0, prodNSP=1; m<=mode; m++)
   {      sumndsize += sp[m]->ndsize; prodNSP *= sp[m]->NSP; 
   }
   
   for(m=1; m<=mode; m++)
   {    if(strstrcase(sp[m]->sOption,"CT")!=NULL) 	{  sumndsize+=1;  prodNSP*=2;} 
   	if(strstrcase(sp[m]->sOption,"ishift=")!=NULL) 	{  sumndsize+=1;  prodNSP*=2;}
   }	
   
   fprintf(ofp,"%d %d %d\n",sumndsize,1,prodNSP);
  
   fprintf(ofp," ");
   for(m=1; m<=mode; m++)
   {   
        for(i=0;i<sp[m]->ndsize;i++) { fprintf(ofp,"%d ", sp[m]->dsize[i]);} 
	if(strstrcase(sp[m]->sOption,"CT")!=NULL) 	fprintf(ofp,"2 ");
        if(strstrcase(sp[m]->sOption,"ishift=")!=NULL) 	fprintf(ofp,"  2   ");
   }
   
   fprintf(ofp,"\n");
   printf("Printed header for sOption0='%s', sumndsize=%d and total sparse points prodNSP=%d\n",sp[0]->sOption, sumndsize,prodNSP );
   return(0);
 }


int printSparseIndices(struct SPARSE **sp, int mode, FILE *ofp, FILE *ofpi)
{  int  i,i_,i1,i2,i3,i4, ileft, m; 
   int  p[5], mrep[5];
 
   for(m=1; m<=mode; m++)
   {   mrep[m]=1; if(strstrcase(sp[m]->sOption,"CT")!=NULL) mrep[m]=2; 	 // repeat for CT
   }
 
   /* 
   int shalf,i1shift,mirep;
    
   i1shift=0;shalf=0;
   for(mirep=1, m=1; m<=mode; m++)
   {    if(strstrcase(sp[m]->sOption,"CT")!=NULL) mirep*=2;		  
   }
   if(strstrcase(sp[0]->sOption,"ishift=")!=NULL) mirep=(mirep)*2 ;      
			    i1shift=0,shalf=0;
			    if(strstrcase(sp[0]->sOption,"ishift=")!=NULL){if( (p==2 && mirep==2) || (p>2 && mirep==4) ) i1shift=sp[1]->ishift, shalf=1; if( i1shift+i1>=sp[1]->NPMULT) i1shift=0, shalf=0;}
		            printsplitindex(p,mirep,1,sp,i1+i1shift,ofp); 
			    if(strstrcase(sp[0]->sOption,"ishift=")!=NULL) fprintf(ofp," %1d ",shalf);  
   
   */
   switch(mode)
   { case 1: 
	     for(i_=0; i_<sp[1]->NSP_; i_++)
	     {   i1= sp[1]->indices[i_];
	         if( sp[1]->flg[i1]==0 ) continue;
		 	
		 	for(p[1]=1; p[1]<=mrep[1]; p[1]++, fprintf(ofp,"  ") )  //  CT: first A then E, else just E
		 	{     printsplitindex(p[1],mrep[1],1,sp,i1,ofp); 
			    
		 	};
		 	fprintf(ofp,"\n");
		 	fprintf(ofpi,"%d \n", i1); // simple indices
	     };
     break;
     case 2: for(i2=0; i2<sp[2]->NPMAX;i2++)
   	     {   if( sp[2]->flg[i2]==0 ) continue;

	     	 for(i_=0; i_<sp[1]->NSP_; i_++)
		 {   i1= sp[1]->indices[i_];
		     if( sp[1]->flg[i1]==0 ) continue;
		 	 							  //    'E ES'   or  'A E'  p='0 1' ;  
		 	for(p[2]=1; p[2]<=mrep[2]; p[2]++, fprintf(ofp,"  ") )
			for(p[1]=1; p[1]<=mrep[1]; p[1]++, fprintf(ofp,"  ") )
			{     printsplitindex(p[1],mrep[1],1,sp,i1,ofp); 
 			      printsplitindex(p[2],mrep[2],2,sp,i2,ofp);
			      
		 	};
			fprintf(ofp,"\n");	    
		        if(i2==0)  fprintf(ofpi,"%d ", i1); // simple indices
	     	 };
   	     };
     break;
     case 3:              
             for(i3=0; i3<sp[3]->NPMAX;i3++)
	     {   if( sp[3]->flg[i3]==0 ) continue;

	     	 for(i_=0; i_<sp[0]->NSP_; i_++)
		 {   	i= sp[0]->indices[i_];
		        if( sp[0]->flg[i]==0 ) continue;
			
	     	 	i1= i%sp[1]->NPMAX; i2= i/sp[1]->NPMAX; 
  					 	
		 	for(p[3]=1; p[3]<=mrep[3]; p[3]++ )  //    'E ES'   or  'A E'  p='0 1' ;  'A1A2 A1E2 E1A2 E1E2'    or p='0 1 2 3'
		 	for(p[2]=1; p[2]<=mrep[2]; p[2]++ )
			for(p[1]=1; p[1]<=mrep[1]; p[1]++ )
			{   								   fprintf(ofp," ");
			   /* char c; */
			   /*   c='E';if(p[1]!=mrep[1])c='A'; fprintf(ofp,"C%c->",c); */   printsplitindex(p[1],mrep[1],1,sp,i1,ofp);  
 			   /*   c='E';if(p[2]!=mrep[2])c='A'; fprintf(ofp,"N%c->",c); */   printsplitindex(p[2],mrep[2],2,sp,i2,ofp); 
			   /*   c='E';if(p[3]!=mrep[3])c='A'; fprintf(ofp,"H%c->",c); */   printsplitindex(p[3],mrep[3],3,sp,i3,ofp); 
			   
		 	};
			fprintf(ofp,"\n");
			if(i3==0) fprintf(ofpi,"%d %d\n", i1,i2); // simple indices
		 };	
	     };
     break;
     case 4: for(i4=0; i4<sp[4]->NPMAX;i4++)
	     {   if( sp[4]->flg[i4]==0 ) continue;

	     	 for(i_=0; i_<sp[0]->NSP_; i_++)
		 {   	i= sp[0]->indices[i_];
		        if( sp[0]->flg[i]==0 ) continue;
			
		        i3= i/(sp[1]->NPMAX*sp[2]->NPMAX); ileft=i%(sp[1]->NPMAX*sp[2]->NPMAX); i2= ileft/sp[1]->NPMAX; i1= ileft%sp[1]->NPMAX;
  					 	
		 	for(p[4]=1; p[4]<=mrep[4]; p[4]++, fprintf(ofp,"  ") )  //    'E ES'   or  'A E'  p='0 1' ;  'A1A2 A1E2 E1A2 E1E2'    or p='0 1 2 3'
		 	for(p[3]=1; p[3]<=mrep[3]; p[3]++, fprintf(ofp,"  ") )
			for(p[2]=1; p[2]<=mrep[2]; p[2]++, fprintf(ofp,"  ") )
			for(p[1]=1; p[1]<=mrep[1]; p[1]++, fprintf(ofp,"  ") )
			{     printsplitindex(p[1],mrep[1],1,sp,i1,ofp); 
 			      printsplitindex(p[2],mrep[2],2,sp,i2,ofp);
			      printsplitindex(p[3],mrep[3],3,sp,i3,ofp); 
			      printsplitindex(p[4],mrep[4],4,sp,i4,ofp); 
		 	};
			fprintf(ofp,"\n");
			if(i4==0)  fprintf(ofpi,"%d %d %d\n", i1,i2,i3); // simple indices
		 };	
	     };
     break;
     
     default:  printf("\nSorry, this mode (%d core dimensions) is not implemented yet \n",mode);
   } 
          
   return(0);
}



//
//  print split index, any 'p' times
//                     for ct:   first comes A, second comes E  (A index is reflected to negative side of the axis, E normal ie on positive side)  done for each dim
//                     for ishift the printing is repeated (p*2) with shifted 'i' 
// ct in N2  = E1A2 E1E2      the same result for  ct in C1  = A1E2 E1E2
//    C  N mi0 H       C  N mi1 H    
// p=           0               1
// half         0  (A)          1    (E)


// shifted + ct
// Ca N mi0 s0 H     C N mi1 s0 H     C+Is N mi0 s1 H    C+Is N mi1 s1 H
// p=          

// ctC1  +  ctN2   =  'A1A2 E1A2 A1E2 E1E2'  
//       C 0 N 0 H    C 1 N 0 H       C 0 N 1 H      C 1 N 1 H


int printsplitindex(int p, int mrep, int m, struct SPARSE **sp, int i, FILE *ofp)
{   int half,d,ip;

    ip=i;  half=1;     // E  to the right half, keep input indexing

    if(strstrcase(sp[m]->sOption,"CT")!=NULL)	// mirror image   f180=y t=1/2, else f180=n t=0 (shifted 1 pos right) 
    {   
        if( p!=mrep )
        {   half=0;      // A  to the left half, reflect and correct for f180 (modify indexing)
    	    if(strstrcase(sp[m]->sOption,"f180=y")!=NULL)    
    	    {  		 ip=(sp[m]->NPMULT-1)-i+0;
    	    } else 
    	    {  if(i!=0)  ip=(sp[m]->NPMULT-1)-i+1; else {half=1;ip=i;}      // exception for 0 point= repeat point as is
    	    }
        }  else 
	{  half=1;       // E to the right half
	}
    }	

    index2splitindices(sp[m], ip);
    for(d=0; d<sp[m]->ndsize;d++)  fprintf(ofp,"%d ",(int)sp[m]->B[d]);
       
    // now extra bit CT
    if(strstrcase(sp[m]->sOption,"CT")!=NULL)  	   fprintf(ofp,"%1d ",half); 

    return(0);
}


int SortIntArray(int size, int *a, int *b, unsigned char flag)
{ int k,j;  int tmp;

   // first make a copy to b from a and then sort b ascending 
   for(k=0;k<size;k++) b[k]=a[k];
      
   if(flag!=0)
   {
      for(k=0;k<size;k++)
      for(j=1;j<size;j++)
      { 
        if( b[j-1] > b[j] )
        {   tmp=b[j-1];
	    b[j-1]=b[j];
	    b[j]=tmp;
        };
      };

   }   
   return(0);   
}      


FILE *Fopen(char *fname, char *mode)
{ static FILE *fp;
	
	if(NULL==(fp=fopen(fname,mode)))
	{ fprintf(stderr,"Error: can't open file %s %s\n",fname,mode); exit(1);
	}
	
	return(fp);
}


char      *strstrcase(const char *s1, const char *s2)
{    int i; char string1[1024],string2[1024]; 
     for(i=0;i<strlen(s1) && i<1024;i++) string1[i]=(char)toupper((int)s1[i]); string1[i]=0;
     for(i=0;i<strlen(s2) && i<1024;i++) string2[i]=(char)toupper((int)s2[i]); string2[i]=0;

  return( strstr(string1,string2) );
}