/****************************************************************/
/*	paropt		03.01.1991				*/
/****************************************************************/
/*	Short Description :					*/
/*	Parameter optimization au program.			*/
/****************************************************************/
/*	Keywords :						*/
/*	optimization, parray, paropt				*/
/****************************************************************/
/*	Description/Usage :					*/
/*								*/
/*	The program starts with the current data set (which	*/
/*	must not have PROCNO=999) and asks for the parameter to	*/
/*	be optimized (type for instance D 1, or P 1, etc.), its	*/
/*	initial value, its increment and the number of		*/
/*	experiments. The program performs the first experiment	*/
/*	with the initial value of the parameter, and changes	*/
/*	the parameter according to the entered increment for	*/
/*	the next cycle. Before the AU program is started, the	*/
/*	current plot region must be defined to the spectrum	*/
/*	region that is to be displayed after each experiment.	*/
/*	This region is stored after each experiment under	*/
/*	PROCNO=999 in the same file, so that the results of the	*/
/*	experiments can be viewed simultaneously. This data set	*/
/*	is displayed while the experiments are going on. The	*/
/*	horizontal scaling should be reset (using the mouse) to	*/
/*	see the entire result .					*/
/****************************************************************/
/*	Author(s) :						*/
/*	Name		: Bruno Guigas				*/
/*	Organisation	: Bruker Analytik			*/
/*	Email		: bruno.guigas@bruker.de		*/
/****************************************************************/
/*	Name		Date	Modification:			*/
/*	bg		921001	changed "2 parname" in		*/
/*				"0 parname"			*/
/*	bg		940609	"ft_mod handling inserted	*/
/*	bg		940914	parameter increments were	*/
/*				limited to 10e-6. Fixed.	*/
/*	bg		950302	handling of u/m units inserted:	*/
/*				m=1.0e-3, u=1.0e-6, e.g.	*/
/*				d1=100m=100millisec.		*/
/*	bg		950817	status si parameter of 999 file	*/
/*				set correctly.			*/
/*	bg		950817	removed last change, set DATMOD	*/
/*	bg		960409	added ZGSAFETY handling,	*/
/*				made SPECSIZ larger.		*/
/*	bg		960916	check total file size,		*/
/*				re-calculate offset1.		*/
/*	bg		970910	ZGSAFETY handling corrected	*/
/*	bg		971126	no VIEWDATA if started with	*/
/*				option no			*/
/****************************************************************/
/*
$Id: paropt,v 1.12 1999/03/29 11:42:14 gsc Exp $
*/

#define USE_AULIBP
#define SPECSIZ 1024*1024 /* maxsize of spec */

int ne, nemax, si, offset1, offset2, npoints, totsize, sprocno, mprocno, i;
int ncneu, ncref, specbuf[SPECSIZ], maxval, minval, bytorder, ftmod;
int zgsafety;
int refreshflag = 1;
float f1p, f2p, offset;
double swp, sf, p1, p1start, p1inc;
char parnam[10], parval[20], parinc[20], *envpnt;

USELASTPARS;
GETCURDATA;
if(procno==999)
{
    Proc_err(DEF_ERR_OPT, "paropt must not be started with PROCNO 999.\n");
    return(-1);
}

if ( !strncmp (cmd,"no",2) )
  refreshflag = 0;

/* Turn zg safety off if on. Turn it back on at the end. */
envpnt = getenv("UXNMR_SAFETY");
zgsafety = 0;
if(envpnt != NULL)
{
    if(strcmp(envpnt, "on") == 0)
    {
	zgsafety = 1;
	CPR_exec("env set UXNMR_SAFETY=off", WAIT_TERM);
    }
}

(void)strcpy(parnam, "p1");
GETSTRING("Enter parameter to modify: ", parnam);
parval[0] = '\0';
GETSTRING("Enter initial parameter value: ", parval);
parinc[0] = '\0';
GETSTRING("Enter parameter increment: ", parinc);
ne=16;	/* init. ne */
GETINT("Enter # of experiments: ", ne);
(void)sscanf(parval, "%lf", &p1start); 
(void)sscanf(parinc, "%lf", &p1inc);

i=strlen(parval)-1;
if(i>0)
{
    f1=1.0;
    if( (parval[i]=='m') || (parval[i]=='D') ) f1=1.0e-3;
    if( (parval[i]=='u') || (parval[i]=='U') ) f1=1.0e-6;
    p1start *= f1;
    i1=strlen(parinc)-1;
    if(i1==0) p1inc *= f1;
    if(i1>0)
    {
	if( (parinc[i1]!='m') && (parinc[i1]!='D') ) 
	{
	    if( (parinc[i1]!='u') && (parinc[i1]!='U') )  p1inc *= f1;
	}
    }
}
i=strlen(parinc)-1;
if(i>0)
{
    if( (parinc[i]=='m') || (parinc[i]=='D') )  p1inc *= 1.0e-3;
    if( (parinc[i]=='u') || (parinc[i]=='U') )  p1inc *= 1.0e-6;
}

p1=p1start;   


/* get some required parameters */
FETCHPARS("OFFSET", &offset);
FETCHPARS("SI", &si);
FETCHPARS("SF", &sf);
FETCHPARS("SW_p", &swp);
FETCHPARS("BYTORDP", &bytorder);	/* get BYTORDP of original data set */
FETCHPLPAR("F1P", &f1p);
FETCHPLPAR("F2P", &f2p);

/* calculate file offsets of plot region */
offset1=(int)(((offset-f1p)*si*sf/swp)+0.5);
offset2=(int)(((offset-f2p)*si*sf/swp)+0.5);
npoints=offset2-offset1+1; 	/* number of points in region */
totsize=ne*npoints;		/* total size of region file */
if(totsize>SPECSIZ)
{
    (void)sprintf(text,"spectrum size becomes too large (> %dk).",SPECSIZ/1024);
    STOPMSG(text);
}

mprocno=999;			/* This is procno of region file */
sprocno=procno;			/* save current procno */
RPROCNO(mprocno);		/* switch to region file */
STOREPAR("SI", totsize);	/* set its real size */
FETCHPAR("SI", &i2);		/* get it back */
if(i2<totsize) i2 *= 2;
if(i2>SPECSIZ)
{
    (void)sprintf(text,"spectrum size becomes too large (> %dk).",SPECSIZ/1024);
    STOPMSG(text);
}
STOREPAR("SI", i2); 		/* if rounded too small */
STOREPAR("DATMOD", 0); 		/* set DATMOD to raw */
STOREPARS("BYTORDP", bytorder);	/* set to BYTORDP of original data set */
ZF;				/* zero spectrum */

STOREPARS("SW_p", swp);		/* insert meaningful sw */

/* now perform ne experiments */
USECURPARS;
TIMES(ne)
    RPROCNO(sprocno)		/* back to acqu. data set */
    (void)sprintf(text, "%s %13.10lf",parnam, p1);/* make CPR command */
    SETCURDATA
    CPR_exec(text, 1);		/* set parameter */
    ZG;				/* acquires data */
    EFP;			/* process */
    FETCHPARS("NC_proc", &ncneu);	/* get NC_proc of new spectrum */
    FETCHPARS("FT_mod", &ftmod);
    if(loopcount1==0) ncref=ncneu; /* init. reference NC */

    /* store F1/F2-reg. in region file; re-calculate offset1 in case a*/
    /* parameter was varied which changes the OFFSET of the spectrum. */
    FETCHPARS("OFFSET", &offset);
    offset1=(int)(((offset-f1p)*si*sf/swp)+0.5);
    if(storeregion(loopcount1, mprocno, offset1,
		 npoints, ncneu-ncref, i2, &specbuf[0])==-1) abort;

    if(ncneu>ncref) ncref=ncneu;	/* set new reference NC */
    RPROCNO(mprocno);		/* switch to mulspec */
    if ( refreshflag )
    {
      VIEWDATA
    }
    sleep(1);
    p1+=p1inc;			/* increment parameter value */
END;
nemax=getmaximum(ne, mprocno, npoints, &specbuf[0],
		 &maxval, &minval); /* search max. in region file */
Proc_err(DEF_ERR_OPT, "paropt finished.\nMax. intensity found at experiment %d.\n %s=%13.10lf",
		nemax, parnam, p1start+p1inc*(double)(nemax-1) );
STOREPARS("YMAX_p", maxval);
STOREPARS("YMIN_p", minval);
STOREPARS("NC_proc", ncref);
STOREPARS("FT_mod", ftmod);

if(zgsafety == 1)
	CPR_exec("env set UXNMR_SAFETY=on", WAIT_TERM);
QUIT



/********* store the plot region in separate file **********/
#define REGSIZ  256*1024 /* maxsize of a region */

storeregion(int curexp, 
	    int mprocno, 
	    int offset1, 
	    int npoints, 
	    int ncdiff, 
	    int filsiz, 
	    int* specbuf
    )
{
    int ret, fd, regbuf[REGSIZ], reim, ii1, ii2;
    char fnam[200], spec[5];
/**************** read plot region in buffer **********************/
    for(reim=0; reim<=1; reim++)
    {
	(void)strcpy(spec, "1r");
	if(reim==1) (void)strcpy(spec, "1i");
	(void)sprintf(fnam, "%s/data/%s/nmr/%s/%d/pdata/%d/%s",
		      disk,user,name,expno,procno,spec);
	/* make source file name */
	
	fd=open(fnam, 0);		/* open source file for read */
	if(fd==-1)
	{
	    Proc_err(DEF_ERR_OPT, "File does not exist:\n%s", fnam);
	    return(-1);
	}
	
	lseek(fd, (long)(offset1*4), 0); /* position of plot region */
	if(npoints > REGSIZ)
	{
	    Proc_err(DEF_ERR_OPT, "F1-F2 region size too large (> 16k).");
	    close(fd);
	    return(-1);
	}
	ret=read(fd, regbuf, npoints*4); /* read plot region in buffer */
	if(ret==-1)
	{
	    Proc_err(DEF_ERR_OPT, "Error when reading file:\n%s", fnam);
	    close(fd);
	    return(-1);
	}
	close(fd);
/******************* end read *************************************/

/*********** rescale new or old data depending on NC_proc *********/
	if(curexp==0) goto writebuffer;	/* 1st region, nothing to rescale */
	if(ncdiff==0) goto writebuffer;	/* same NCs, nothing to rescale */
	if(ncdiff<0)	/* then new region must be rescaled */
	{
	    for(ii2=0; ii2<-ncdiff; ii2++)
		for(ii1=0; ii1<npoints; ii1++) regbuf[ii1] /= 2; 
	}
	if(ncdiff>0)	/* then old file must be rescaled */
	{
	    (void)sprintf(fnam, "%s/data/%s/nmr/%s/%d/pdata/%d/%s",
			  disk,user,name,expno,mprocno,spec);  /* make destination file name */
	    
	    fd=open(fnam, 2);		/* open destination file for r+w */
	    if(fd==-1)
	    {
		Proc_err(DEF_ERR_OPT, "File does not exist:\n%s", fnam);
		return(-1);
	    }
	    if(filsiz > SPECSIZ)
	    {
		Proc_err(DEF_ERR_OPT, "spectrum size too large (> 1024k).");
		close(fd);
		return(-1);
	    }
	    ret=read(fd, specbuf, filsiz*4); /* read entire file */
	    if(ret==-1)
	    {
		Proc_err(DEF_ERR_OPT, "Error when reading file:\n%s", fnam);
		close(fd);
		return(-1);
	    }
	    for(ii2=0; ii2<ncdiff; ii2++)
		for(ii1=0; ii1<filsiz; ii1++) specbuf[ii1] /= 2; /* scale file */
	    lseek(fd, 0L, 0);		/* rewind */
	    ret=write(fd, specbuf, filsiz*4); /* write file */
	    if(ret != filsiz*4)
	    {
		Proc_err(DEF_ERR_OPT, "Error when writing file:\n%s.\n%d bytes written", fnam, ret);
		close(fd);
		return(-1);
	    }
	    close(fd);
	} /* end if ncdiff>0 */
/******************* end rescale **********************************/

/*********** write buffer into region *********************/
writebuffer:
	(void)sprintf(fnam, "%s/data/%s/nmr/%s/%d/pdata/%d/%s",
		      disk,user,name,expno,mprocno,spec);
/* make destination file name */

	fd=open(fnam, 1);		/* open destination file for write */
	if(fd==-1)
	{
	    Proc_err(DEF_ERR_OPT, "File does not exist:\n%s", fnam);
	    return(-1);
	}
	
	lseek(fd, (long)(curexp*npoints*4), 0);/* where plot region goes now */
	ret=write(fd, regbuf, npoints*4); /* write buffer to dest. file */
	if(ret != npoints*4)
	{
	    Proc_err(DEF_ERR_OPT, "Error when writing file:\n%s.\n%d bytes written", fnam, ret);
	    close(fd);
	    return(-1);
	}
	close(fd);
    } /* end re-im-loop */
    return(0);
} /* end storeregion */
/*********************************************************************/


/*********************************************************************/
/* search for maximum */
getmaximum(int ne, 
	   int mprocno, 
	   int npoints, 
	   int *specbuf,
	   int *maxval,
	   int *minval
    )
{
    int fd, i, k, max, min, imax, curval, ret;
    char fnam[200];
    
    (void)sprintf(fnam, "%s/data/%s/nmr/%s/%d/pdata/%d/%s",
		  disk,user,name,expno,mprocno,"1r");
    /* make region file name */
    
    fd=open(fnam, 0);		/* open region file for r */
    if(fd==-1)
    {
	Proc_err(DEF_ERR_OPT, "File does not exist:\n%s", fnam);
	return(-1);
    }
    ret=read(fd, specbuf, ne*npoints*4); /* read all regions in buffer */
    if(ret==-1)
    {
	Proc_err(DEF_ERR_OPT, "Error when reading file:\n%s", fnam);
	close(fd);
	return(-1);
    }
    close(fd);
    
/* max. search */
    max=0; imax=0; min=2147483646; /* 2^31 */
    for(i=0; i<ne; i++)
    {
	for(k=0; k<npoints; k++)
	{
	    curval=specbuf[npoints*i+k];
	    if(curval > max) {max=curval; imax=i+1;}
	    if(curval < min) {min=curval;}
	} /* end for k */
    } /* end for i */
    *maxval=max; *minval=min;
    return(imax);
} /* end getmaximum */