stratinv.c.bac

su 的源代码库

	 /* linear interpolating */
	 a = (t1Aux[i + 1] - t1Aux[i]) / (offAux[i + 1] - offAux[i]);
	 b = t1Aux[i] - a * offAux[i];
	 j = k;
	 while (j <= offAux[i + 1])
	 {
	    t1Mute[j - 1] = a * j + b;
		/* DD
		fprintf(stdout, "%d %f\n", j, t1Mute[j - 1]);*/
	    j++;
	 }
      }
      /* freeing */
      free1float(t1Aux);
      free1float(offAux);
   }

	/* DD 
	exit(-1);*/

   /* number of misfit samples */ 
   nDM = NINT((t2 - t1) / dt) + 1;
   nTotalSamples = nDM * info->nR; 

   if (info->verbose)
      fprintf(stderr, "Number of data points in misfit computation: %d\n",
	      nTotalSamples);
   
   /* # of iterations */
   if (!getparint("maxiter", &maxIter)) maxIter = 4;
   if (!getparfloat("stdnoise", &noiseVar)) 
      err("Specify noise standard deviation\n");
   
   /* some hard-coded parameters */
   fR = 1; info->wR = 2 * PI * fR;   /* reference frequency */
   
   /* how many layers */
   fp = fopen(modelFile,"r");
   if (fp == NULL)
      err("No model file!\n");
   
   info->nL = 0;
   depth = 0;
   while (fscanf(fp, "%f %f %f %f %f %f\n", 
		 &aux, &aux, &aux, &aux, &aux, &aux) != EOF)
      info->nL++;
   info->nL--;
   fclose(fp);
   
   /* considering the unknown layers */
   limRange = NINT((limZ[1] - limZ[0]) / dZ);
   /* DD */
   fprintf(stderr,"Number of layers - 1 (info->nL): %d\n", info->nL);
   fprintf(stderr,"Target width: %d\n", limRange);

   /* basic time-frequency stuff */
   info->nSamples = NINT(tMod / dt) + 1;
   info->nSamples = npfar(info->nSamples);

   tMod = dt * (info->nSamples - 1);
   info->dF = 1. / (tMod);

   /* adjusting f1 and f2 */
   aux = info->dF;
   while (aux < info->f1) aux += info->dF;
   info->f1 = aux;
   while (aux < info->f2) aux += info->dF;
   info->f2 = aux;
   
   info->nF = NINT(info->f2 / info->dF) - NINT(info->f1 / info->dF) + 1; 
   if (info->nF%2 == 0) 
   {
      info->f2 += info->dF;
      info->nF++;
   }
   initF = NINT(info->f1 / info->dF);
   lastF = NINT(info->f2 / info->dF);
   if (info->nFreqProc == 0)
   {
      info->nFreqProc = NINT((float) info->nF / (float) nProc + .5);
      if (info->nFreqProc > info->nF) info->nFreqProc = info->nF;
   }
   else
      while (info->nFreqProc > info->nF) info->nFreqProc /= 2;
   nFreqPart = NINT((float) info->nF / (float) info->nFreqProc + .5);

   /* DD 
   fprintf(stderr, "nF -> %d f1 %f f2 %f df %f\n", 
	   info->nF, info->f1, info->f2, info->dF);*/

   /* memory allocation */
   alpha = alloc1float(info->nL + 1);
   beta = alloc1float(info->nL + 1);
   rho = alloc1float(info->nL + 1);
   if (vpFrechet) alphaMean = alloc1float(info->nL + 1);
   if (vsFrechet) betaMean = alloc1float(info->nL + 1);
   if (rhoFrechet) rhoMean = alloc1float(info->nL + 1);
   if (vpFrechet) alpha0 = alloc1float(info->nL + 1);
   if (vsFrechet) beta0 = alloc1float(info->nL + 1);
   if (rhoFrechet) rho0 = alloc1float(info->nL + 1);
   qP = alloc1float(info->nL + 1);
   qS = alloc1float(info->nL + 1);
   thick = alloc1float(info->nL + 1);

   /* DD 
   fprintf(stderr, "nDM : %d\n", nDM);*/
   
   dataObs = alloc2float(nDM, info->nR);
   resCD = alloc2complex(info->nF, info->nR);
   CD = alloc1float(nDM * (nDM + 1) / 2);

   if (PRIOR)
   {
      if(vpFrechet)
         CMvP = alloc1float(limRange * (limRange + 1) / 2);
      if(vsFrechet)
         CMvS = alloc1float(limRange * (limRange + 1) / 2);
      if(rhoFrechet)
         CMrho = alloc1float(limRange * (limRange + 1) / 2);
   }

   grad0 = alloc1float(numberPar * limRange);
   grad1 = alloc1float(numberPar * limRange);
   search = alloc1float(numberPar * limRange);

   /* distributed stuff */
   processes = alloc1int(nProc);
   procInfo = alloc2int(2, nProc);
   statusFreq = alloc2int(3, nFreqPart);

   /* defining frequency partitions */
   for (k = initF, i = 0; i < nFreqPart; i++, k += info->nFreqProc)
   {
      statusFreq[i][0] = k;
      statusFreq[i][1] = MIN(k + info->nFreqProc - 1, lastF);
      statusFreq[i][2] = 0;
   }
          
   /* reading the model file */
   fp = fopen(modelFile,"r");
   if (info->verbose)
      fprintf(stderr,"  Thickness     rho     vP     qP    vS     qS\n");
   
   for (k = 0; k < info->nL + 1; k++)
   {
      fscanf(fp, "%f %f %f %f %f %f\n", &thick[k], &rho[k], &alpha[k], 
	     &qP[k], &beta[k], &qS[k]);
      if (info->verbose)
	 fprintf(stderr,"   %7.4f      %4.3f   %3.2f  %5.1f  %3.2f  %5.1f\n",
		 thick[k], rho[k], alpha[k], qP[k], beta[k], qS[k]);
      if (vpFrechet) alphaMean[k] = alpha[k];
      if (vsFrechet) betaMean[k] = beta[k];
      if (rhoFrechet) rhoMean[k] = rho[k];
   }
   fclose(fp);

   /* setting lim[0] and lim[1] */
   for (depth = thick[0], i = 1; i <= info->nL; depth += thick[i], i++)
   {
      if (NINT(depth / dZ) <= NINT(limZ[0] / dZ)) lim[0] = i;
      if (NINT(depth / dZ) < NINT(limZ[1] / dZ)) lim[1] = i;
   }
   lim[1]++;

   /* DD 
   fprintf(stderr, "lim[0] %d lim[1] %d\n", lim[0], lim[1]);*/
   
   /* some modeling parameters */
   /* slowness increment */
   info->dU = (info->u2 - info->u1) / (float) info->nU;

   /* imaginary part of frequency for damping wrap-around */
   info->tau = log(info->tau) / tMod;
   if (info->tau > TAUMAX) info->tau = TAUMAX;
   /* DD 
   fprintf(stderr, "TAU %f\n", info->tau);*/

   /* reading data and model covariance matrixes */
   inputCovar(corrDataFile, corrModelFile);
   
   /* starting inverse procedure */
   /* opening report file */
   fp = fopen("report", "w");
   if (vpFrechet) fprintf(fp, "Inversion with respect to p-wave velocity\n");
   if (vsFrechet) fprintf(fp, "Inversion with respect to s-wave velocity\n");
   if (rhoFrechet) fprintf(fp, "Inversion with respect to density\n");
   fclose(fp);

   /* some initializations */
   modCount = 0;
   gradCount = 0;
   notDone = 1;
   dataIsFit = 0;
   oFNorm = 1;

   /* reading "observed" data */
   inputData(dataFile);

   /* more terms in INFO structure */
   info->lim[0] = lim[0];
   info->lim[1] = lim[1];
   info->limRange = limRange;
   info->numberPar = numberPar;
   info->vpFrechet = vpFrechet;
   info->vsFrechet = vsFrechet;
   info->rhoFrechet = rhoFrechet;
   
   /* modeling data at input model and computing its gradient */
   sprintf(info->id, "Initial guess");
   oF0 = modeling(); 
   sprintf(info->id, "Gradient at initial guess");
   gradient(grad0); 
   modCount++; gradCount++;

   /* computing the search direction */
   for (i = 0; i < numberPar * limRange; i++)
   {
      grad1[i] = grad0[i];
   }
   slope = newSearch(grad0, grad1, search);
      
   /* reporting */
   fp = fopen("report", "a");
   fprintf(fp,"-----------------------\n");
   fprintf(fp,"gradient at iteration [OF:%d][GR:%d]:\n", 
	   modCount, gradCount);
   for (i = 0; i < numberPar * limRange; i++)
      fprintf(fp, "grad[%d]: %f\n", i, grad0[i]);
   fprintf(fp,"search at iteration [OF:%d][GR:%d]:\n", 
	   modCount, gradCount);
   for (i = 0; i < numberPar * limRange; i++)
      fprintf(fp, "search[%d]: %f\n", i, search[i]);
   fprintf(fp,"-----------------------\n");
   fclose(fp);
  
   while (notDone)
   {
      /* computing step length via a line search */
      /* just the real part of the search direction is used */
      oF1 = lineSearch(search, oF0, slope);  
	 
      /* computing the gradient at the new model */
      sprintf(info->id, "Gradient at update model");
      gradient(grad1);
      gradCount++;

      /* and the new search direction */
      slope = newSearch(grad0, grad1, search);
      
      /* reporting */
      fp = fopen("report", "a");
      fprintf(fp,"-----------------------\n");
      fprintf(fp,"gradient at iteration [OF:%d][GR:%d]:\n", 
	      modCount, gradCount);
      for (i = 0; i < numberPar * limRange; i++)
	 fprintf(fp, "grad[%d]: %f\n", i, grad1[i]);
      fprintf(fp,"search at iteration [OF:%d][GR:%d]:\n", 
	      modCount, gradCount);
      for (i = 0; i < numberPar * limRange; i++)
	 fprintf(fp, "search[%d]: %f\n", i, search[i]);
      fprintf(fp,"-----------------------\n");
      fclose(fp);

      /* checking stopping criterion (just the gradient and */
      /* objective function for now ) */
      notDone = stop(grad0, grad1, oF0, oF1, maxIter);

      /* switching gradients */
      auxP = grad0;
      grad0 = grad1;
      grad1 = auxP;

      aux = oF0;
      oF0 = oF1;  
      oF1 = aux;   
   }
}