invtools.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

/*  Function inputData()                                        */
/*                                                              */
/*  Seismic data input and transformation to the frequency      */
/*  domain                                                      */
/*                                                              */
/*  Input parameters:                                           */
/*  dataFile...............file name of input data              */
/*  nSamples...............number of samples per trace          */
/*                         global                               */
/*  nR.....................number of offsets                    */
/*                         global                               */
/*                                                              */
/*  Output parameters:                                          */
/*  dataObs................observed data in time domain         */
/*                         global                               */
/*                                                              */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
#include "stratInv.h"
segy tr;		        /* reading data */
void inputData(char* dataFile)
{
   /* declaration of variables */
   int iS, iR, iF, iF1, iF2;    /* generic counters */
   int ns;			/* # of samples */
   int wL;                      /* window length */
   float *buffer = NULL;	/* to input data */
   float window;                /* windowing purposes */
   complex *bufferC = NULL;	/* to Fourier transform the input data */
   FILE *fp;			/* input file */

   /* memory for bufferC */
   bufferC = alloc1complex(info->nSamples / 2 + 1);
   
   fp = fopen(dataFile,"r");
   if (fp == NULL)
      err("Can't open input data file!\n");

   for (iR = 0; iR < info->nR; iR++)
   {
      fgettr(fp, &tr);
      ns = tr.ns;
      /* DD 
      fprintf(stderr, "ns %d\n", ns);*/

      /* allocating memory */
      if (iR == 0) buffer = alloc1float(MAX(ns, info->nSamples));

      /* reseting */
      for (iS = 0; iS < MAX(ns, info->nSamples); iS++) buffer[iS] = 0;
      memcpy(buffer, tr.data, ns * FSIZE);
      
      /* buffer -> dataObs and compensating for complex frequency */
      for (iS = 0; iS < info->nSamples; iS++)
      {
	 buffer[iS] *= exp(-info->tau * iS * dt);
	 /* DD 
	 fprintf(stderr, "buffer[%d] : %f\n", iS, buffer[iS]);*/
      }

      /* going to the Fourier domain */
      pfarc(-1, info->nSamples, buffer, bufferC);
      
      /* windowing (PERC_WINDOW) spectrum */
      iF1 = NINT(info->f1 / info->dF);
      iF2 = NINT(info->f2 / info->dF);
      wL = info->nF * PERC_WINDOW / 2;
      wL = 2 * wL + 1;
      for (iS = 0, iF = 0; iF < info->nSamples / 2 + 1; iF++)
      {
	 window = 0;
	 if (iF < iF1 || iF >= iF2)
	 {
	    bufferC[iF] = cmplx(0, 0);
	 }
	 else if (iF - iF1 < (wL - 1) / 2)
	 {
	    window =
	       .42 - .5 * cos(2 * PI * (float) iS / ((float) (wL - 1))) +
		  .08 * cos(4 * PI * (float) iS / ((float) (wL - 1)));
	    bufferC[iF].r *= window; bufferC[iF].i *= window;
	    iS++;
	 }
	 else if (iF - iF1 >= info->nF - (wL - 1) / 2)
	 {
	    iS++;
	    window =
	       .42 - .5 * cos(2 * PI * (float) iS / ((float) (wL - 1))) +
		  .08 * cos(4 * PI * (float) iS / ((float) (wL - 1)));
	    bufferC[iF].r *= window; bufferC[iF].i *= window;
	 }
      }

      /* going back to time domain */
      pfacr(1, info->nSamples, bufferC, buffer);

      /* copying to dataObs within target window and scaling */
      for (iF = 0, iS = NINT(t1 / dt); iS <= NINT(t2 / dt); iS++, iF++)
      {
	 dataObs[iR][iF] = (scaleData * buffer[iS]) / (float) info->nSamples;
	 /* DD 
	 fprintf(stderr, "%d %d %f %f %f %f\n", iR, iF, dataObs[iR][iF], 
		 info->f1, info->f2, scaleData);*/
      }
   }
   /* DD 
   fprintf(stderr, "energy %f\n", auxm1 / (nDM * info->nR));
   fwrite(&dataObs[0][0], sizeof(float), nDM * info->nR, stdout);
   exit(-1);*/
   
   /* freeing memory */
   free1float(buffer);
   free1complex(bufferC);

   fclose(fp); 
}
/*  Function inputCovar()                                       */
/*                                                              */
/*  Input of data and model covariance matrixes                 */
/*                                                              */
/*  Input parameters:                                           */
/*  corrDataFile...............file name with data covariance   */
/*  corrModelFile[0 1 2].......file name with model covariance  */
/*                             [0] : (P-wave velocities)        */
/*                             [1] : (S-wave velocities)        */
/*                             [2] : (Density)                  */
/*                             used in case PRIOR = 1           */
/*                                                              */
/*  Output parameters:                                          */
/*  CD.....................data covariance matrix               */
/*  CM.....................model covariance matrix              */
/*                         outputted in case PRIOR = 1          */ 
/*                         (CMvP, CMvS, CMrho)                  */
/*                                                              */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
void inputCovar(char* corrDataFile, char* corrModelFile[3])
{
   /* declaration of variables */
   int iS, iS1, i;              /* generic counters */
   float *buffer = NULL;	/* to input covariances */
   FILE *fp;			/* input file */

   noiseVar = 0;
   /* inputting data covariances */
   if (DATACOV)
   {
      fp = fopen(corrDataFile, "r");
      if (fp == NULL)
	 err("Can't read data covariance file %s\n", corrDataFile);
   }
   
   /* allocating memory for auxiliary buffer */
   if (DATACOV) buffer = alloc1float(nDM);
   
   for (i = 0, iS = 0; iS < nDM; iS++)
   {
      if (DATACOV)
      {
	 fread(buffer, sizeof(float), nDM, fp);
         if (noiseVar < 1. / buffer[iS]) 
            noiseVar = 1. / buffer[iS];
      }
      for (iS1 = iS; iS1 < nDM; iS1++, i++)
      {
	 if (DATACOV)
	    CD[i] = buffer[iS1];
	 else
	    if (iS == iS1) CD[i] = 1;
	    else CD[i] = 0;
      }
   }
   
   /* freeing memory */
   if (DATACOV) free1float(buffer);

   if (PRIOR)
   {
      /* allocating memory for auxiliary buffer */
      buffer = alloc1float(limRange);
      
      if (vpFrechet)
      {
	 fp = fopen(corrModelFile[0], "r");
	 if (fp == NULL)
	    err("Can't read model covariance file %s\n", corrModelFile[0]);
	 
	 for (i = 0, iS = 0; iS < limRange; iS++)
	 {
	    fread(buffer, sizeof(float), limRange, fp);
	    for (iS1 = iS; iS1 < limRange; iS1++, i++)
	    {
	       CMvP[i] = buffer[iS1];
	    }
	 }
	 fclose(fp);
      }

      if (vsFrechet)
      {
	 fp = fopen(corrModelFile[1], "r");
	 if (fp == NULL)
	    err("Can't read model covariance file %s\n", corrModelFile[1]);
	 
	 for (i = 0, iS = 0; iS < limRange; iS++)
	 {
	    fread(buffer, sizeof(float), limRange, fp);
	    for (iS1 = iS; iS1 < limRange; iS1++, i++)
	    {
	       CMvS[i] = buffer[iS1];
	    }
	 }
	 fclose(fp);
      }

      if (rhoFrechet)
      {
	 fp = fopen(corrModelFile[2], "r");
	 if (fp == NULL)
	    err("Can't read model covariance file %s\n", corrModelFile[2]);
	 
	 for (i = 0, iS = 0; iS < limRange; iS++)
	 {
	    fread(buffer, sizeof(float), limRange, fp);
	    for (iS1 = iS; iS1 < limRange; iS1++, i++)
	    {
	       CMrho[i] = buffer[iS1];
	    }
	 }
	 fclose(fp);
      }
      /* freeing memory */
      free1float(buffer);
   }
}
/*  Function lineSearch()                                       */
/*                                                              */
/*  Perform a cubic line search looking for a "reasonable"      */
/*  step length                                                 */
/*                                                              */
/*  Input parameters:                                           */
/*  search.................search direction                     */
/*  fxc....................objective function value at          */
/*                         initial model                        */
/*  slope..................gradient . search                    */
/*                                                              */
/*  Output parameters:                                          */
/*  alpha, beta or rho....updated model                         */
/*                        global                                */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
/*  Adapted from Doug Hart's implementation                     */
float lineSearch(float *search, float fxc, float slope)
{
   int cnt=0, i;                         /* counters */
   int flag;                             /* ==1, minimum relative change */
                                         /* quitting LS */
   static float lambda = INITIAL_LAMBDA; 	
					 /* initial step length */
   float fOld;       
   float lambdaOld;                      /* used to restore the model */
   float lambda2, lambdaprev, lambdaprev2, lambdatemp;
   float fnew, fprev;
   float f2, f1;
   float c, cm11, cm12, cm21, cm22;
   float a, b, disc;                     /* used in interpolation */
 
   /*if (lambda < LAMBDA_MIN) lambda = LAMBDA_MIN;*/

   /* copying model to auxiliary buffers */
   for (i = 0; i <= info->nL; i++)
   {
      if (vpFrechet) alpha0[i] = alpha[i];
      if (vsFrechet) beta0[i] = beta[i];
      if (rhoFrechet) rho0[i] = rho[i];
   }
   
   /* new trial point a full step away */
   flag = update(lambda, search);
   if (flag) return(fxc);    /* marginal change */
   
   /* evaluating new model */
   sprintf(info->id, "Modeling at line search");
   fnew = modeling(); fOld = fnew;
   cnt++; modCount++;
   
   /* Goldstein's test for to prevent small step sizes */
   /* See Fletcher, Practical Methods of Optimization, page 26ff. */
   while (fnew < (fxc + (1 - ALPHA) * lambda * slope)) 
   {
      fOld = fnew;
      lambdaOld = lambda;
      lambda *= 1.5;   /* increase step size by factor of 1.5 */
      flag = update(lambda, search);
      if (flag) return(fnew);              /* marginal change */
      sprintf(info->id, "Modeling at line search");
      fnew = modeling();      /* evaluating new model */
      cnt++; modCount++;
      if (cnt >= MAX_ITER_LS) 
      {
	 if (fOld < fnew)
	 {
	    restore(lambdaOld, search);
	    fnew = fOld;
	 }
	 return(fnew);
      }
   }

   if (fOld < fnew)
   {
      restore(lambdaOld, search);
      lambda = lambdaOld;
      fnew = fOld;
   }

   /* Armijo's test for lambda too large (the backtracking algorithm) */
   /* See Dennis and Schnabel, Numerical Methods for Unconstrained   */
   /* Optimization and Nonlinear Equations, page 126ff.              */
   if(fnew > (fxc + ALPHA * lambda * slope)) 
   {
      fOld = fnew;
      lambdaOld = lambda;
      /* first try a quadratic fit */   
      lambda2 = lambda * lambda;
      f1 = fnew - fxc - slope * lambda;
      lambdatemp = -slope * lambda2 / (2.0 * f1);
      lambdaprev = lambda;
      fprev      = fnew;
      if(lambdatemp < (0.1 * lambda))
         lambda *= 0.1;
      else
         lambda  = lambdatemp;

      /* evaluating new model */
      flag = update(lambda, search);
      if (flag) return(fnew);   /* marginal change */
      sprintf(info->id, "Modeling at line search");
      fnew = modeling();
      cnt++; modCount++;

      if (cnt >= MAX_ITER_LS) 
      {
	 if (fOld < fnew)
	 {
	    restore(lambdaOld, search);
	    fnew = fOld;
	 }
	 return(fnew);
      }
      while(fnew > (fxc + ALPHA * lambda * slope)) 
      { 
	 fOld = fnew;
	 lambdaOld = lambda;
         /* try cubic models if quadratic failed */
         lambda2     = lambda * lambda;
         f1          = fnew - fxc - slope * lambda;
         lambdaprev2 = lambdaprev * lambdaprev;
         f2          = fprev - fxc - lambdaprev * slope;
         c           = 1.0 / (lambda - lambdaprev);
         cm11        = 1.0 / lambda2;
         cm12        = -1.0 / lambdaprev2;
         cm21        = -lambdaprev / lambda2;
         cm22        = lambda/lambdaprev2;
         a           = c * (cm11 * f1 + cm12 * f2);
         b           = c * (cm21 * f1 + cm22 * f2);
         disc        = b * b - 3.0 * a * slope;
         if((fabs(a) > MINFLOAT) && (disc > MINFLOAT)) 
	                                     /* legitimate cubic fit */
            lambdatemp = (-b + sqrt(disc)) / (3.0 * a);
         else                                /* degenerate quadratic fit */
            lambdatemp = - slope * lambda2 / (2.0 * f1);
         if(lambdatemp >= 0.5 * lambda)
            lambdatemp = 0.5 * lambda;
         if(lambdatemp < (0.1 * lambda))
            lambda *= 0.1;
         else
            lambda  = lambdatemp;

	 /* evaluating new model */
	 flag = update(lambda, search);
	 if (flag) return(fnew);    /* marginal change */
	 sprintf(info->id, "Modeling at line search");
	 fnew = modeling();   
	 cnt++; modCount++;
	 if (cnt >= MAX_ITER_LS) 
	 {
	    if (fOld < fnew)
	    {
	       restore(lambdaOld, search);
	       fnew = fOld;