ngen_wsta.c

su 的源代码库

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

#include "extern.h"
#include "pvm.h"
#include "stdio.h"
#include "su.h"
#include "segy.h"
#include "header.h"
#include "values.h"

/************************************************  file ngen_wsta.c  *******/
/* 111494 */

segy tr;		/* SEGY DATA */

/*    
    Static external variables for routine eval and Go_uphill
*/

static float **CROSS;		/* crosscorrelation	*/
static int ***POINTER;		/* set of pointers	*/
static int *FOLDCMP;		/* fold of each CMP	*/
static int **DERIVATIVE;	/* to compute the derivatives */
static int ***POINT_CROSS;      /* pointers             */
static float *taper;		/* taper for trend detection */
static int **CONSIST;           /* for the penalty function */
static int n_of_consist;        /* """ """ """"""" """""""" */
static char msg[80];            /* used when printing error message     */

eval(how_many, genes)
int how_many;	/* How many members will be evaluated		*/
int genes;      /* number of elements in vect                   */
{
	int *SOURCE, *RECEIVER;		/* POINTERS		*/
	static int firstflag = 1;
	double static_j, static_i; 	/* definig total statics*/
	register int i, j, k;
	int NTRACES;			/* # of traces read	*/
	int shift;
	int ipoint;
	int isample, isource, itrace, ireceiver, icmp, cmp_current, cmp_last;
				/* counters and similar stuff		*/
	int i_to_calc;		/* defines member to be evaluated	*/
        char suff[10];          /* defines data file                    */
	FILE *fp;		/* Data file				*/

	double aux;		/* auxiliar quantity			*/
/*
    Setting the true solution
*/
	if (firstflag) 
	{
                fp = fopen(datafile,"r");

                if (fp == NULL)
                {
                        sprintf(msg, "Can't open %s", datafile);
                        Error(msg);
                }
		if (verbose)
                	fprintf(stderr,"Subpopulation %d read data file %s\n",instance, datafile);

/*
    Memory allocation for vector:
	POINTER[0...NCMP-1][0...FOLD-1][SOURCE_ID, RECEIVER_ID]
*/
		POINTER = alloc3int(2,MAXFOLD,NCMP);
		if (POINTER == NULL) Error("Allocation failed for POINTER");
/*
    Source vector used to index source statics
*/
		SOURCE = alloc1int(NSOURCES);
		if (SOURCE == NULL) Error("Allocation failed for SOURCE");
/*
    Receiver vector used to index receiver statics
*/
		RECEIVER = alloc1int(NRECEIVERS);
		if (RECEIVER == NULL) Error("Allocation failed for RECEIVER");
/*
    Foldcmp vector used for fold of each CMP 
*/
		FOLDCMP = alloc1int(NCMP);
		if (FOLDCMP == NULL) Error("Allocation failed for FOLDCMP");
/*
    Allocating memory for correlation and for consistency vector
*/
                CROSS = alloc2float(TOTAL_LAG+1, MAXFOLD*(MAXFOLD-1)/2);
                if (CROSS == NULL) Error("Allocation failed for CROSS");

                CONSIST = alloc2int(2, NSOURCES);
                if (CONSIST == NULL) Error("Allocation failed for CONSIST");

                to_filter = alloc2int(3, NRECEIVERS);
                if (to_filter == NULL) Error("Allocation failed for to_filter");

		taper = alloc1float(NSOURCES + NRECEIVERS);
                if (taper == NULL) Error("Allocation failed for taper");
/*
    Reading headers. That's just what I need here...
    And copying traces and headers to temporary files 
    NOTE THAT CMPs MUST BE IN ASCENDING ORDER 
*/
		NTRACES = 0;
		isource = 0;
		ireceiver = 0;
		icmp = -1;
		cmp_current = 0;
		cmp_last = 0;
		while (fgettr(fp,&tr)) { 
		 	cmp_current = tr.cdp;
			if (cmp_current != cmp_last)
			{
				if (icmp >= 0) FOLDCMP[icmp] = itrace;
				icmp++;
				itrace = 0;
				cmp_last = cmp_current;
			}	
		 	POINTER[icmp][itrace][0] = tr.sx;	
	                POINTER[icmp][itrace][1] = tr.gx;

			if (icmp == 0 && itrace == 0)
				SOURCE[isource] = POINTER[icmp][itrace][0];

			else if (!(icmp == 0 && itrace == 0) && SOURCE[isource] != POINTER[icmp][itrace][0])
			{
				i = isource - 1;
				while (i >= 0 && SOURCE[i] != POINTER[icmp][itrace][0])
					i--;
				if (i < 0)
				{
					isource++;
					SOURCE[isource] = POINTER[icmp][itrace][0];
				}
			}

                        if (icmp == 0 && itrace == 0)
                                RECEIVER[ireceiver] = POINTER[icmp][itrace][1];

			else if (!(icmp == 0 && itrace == 0) && RECEIVER[ireceiver] != POINTER[icmp][itrace][1])
			{
				i = ireceiver - 1;
				while (i >= 0 && RECEIVER[i] != POINTER[icmp][itrace][1])
					i--;
				if (i < 0)
				{
					ireceiver++;
					RECEIVER[ireceiver] = POINTER[icmp][itrace][1];
				}
			}

			itrace++;
			NTRACES++;
		};
		FOLDCMP[icmp] = itrace;	/* last CMP */
		fclose(fp);
		if (verbose)
		{
			fprintf(stderr,"Subpopulation %d has %d traces read in %d CMPs\n",instance,NTRACES,icmp+1);
			fprintf(stderr,"Subpopulation %d has %d source statics and %d receiver statics\n",instance,isource+1,ireceiver+1);
		}
/*
    Enhancing the pointers
*/
		for (icmp = 0; icmp < NCMP; icmp++)
		{
			for (itrace = 0; itrace < FOLDCMP[icmp]; itrace++) 
			{
				isource = 0;
				while (SOURCE[isource] != POINTER[icmp][itrace][0])
                                        isource++;

				POINTER[icmp][itrace][0] = isource;

				ireceiver = 0;
			 	while (RECEIVER[ireceiver] != POINTER[icmp][itrace][1])
					ireceiver++;	
				POINTER[icmp][itrace][1] = ireceiver;
			}
		}
/*
    Detecting where source and receiver positions are the same.
    In this case I will constrain them to be similar by a 
    Lagrange multiplier
*/
		for (n_of_consist = 0, isource = 0; isource < NSOURCES; isource++) 
		{
			ireceiver = 0;
			while (ireceiver < NRECEIVERS && SOURCE[isource] != 
							 RECEIVER[ireceiver])
				ireceiver++;

			if (ireceiver != NRECEIVERS)
			{
				CONSIST[n_of_consist][0] = isource;
				CONSIST[n_of_consist][1] = ireceiver;
/*
    Source statics pointed by isource and receiver statics pointed
    by ireceivers should be very similar 
*/
				n_of_consist++;
			}
		}
/*
    Now computing the taper for trend detection
*/
		for (i = 0; i < NSOURCES + NRECEIVERS; i++)
			taper[i] = sin(PI * i / ((float) NSOURCES + NRECEIVERS - 1));
/*
    Computing the right indexes for interpolation
*/
		for (ireceiver = 0; ireceiver < NRECEIVERS; ireceiver++)
		{
			i = 0;
			while (i < NRECEIVERS && 
			       (RECEIVER[ireceiver] - dx) != RECEIVER[i])
				i++;

                        j = 0;
                        while (j < NRECEIVERS && 
                               (RECEIVER[ireceiver] + dx) != RECEIVER[j])
                        	j++;

			if (i != NRECEIVERS)
				to_filter[ireceiver][0] = NSOURCES + i;
			else
				to_filter[ireceiver][0] = NSOURCES + ireceiver;

			if (j != NRECEIVERS)
				to_filter[ireceiver][1] = NSOURCES + j;
			else
				to_filter[ireceiver][1] = NSOURCES + ireceiver;
			
			if (i == NRECEIVERS || j == NRECEIVERS)
				to_filter[ireceiver][2] = 2; 
			else
				to_filter[ireceiver][2] = 3;
		}
/*
    Freeing memory
*/
		free1int(SOURCE);
		free1int(RECEIVER);	/* Neither for this pointers	   */

		if (genes != NSOURCES + NRECEIVERS)
			Error("Improper bit string at objective function evaluation");
		firstflag = 0;		/* end of initialization	   */
	}
/*
    Resetting...
*/

	for (i=0; i < how_many; i++)
		eval_returned[i] = 0.;
	
	rewind(Xfp);
	for (icmp = 0; icmp < NCMP; icmp++)
	{
	   fread(CROSS[0],sizeof(float),(TOTAL_LAG+1)*(MAXFOLD*(MAXFOLD-1)/2),Xfp);
	   for (i_to_calc=0; i_to_calc < how_many; i_to_calc++)
           {
		for (ipoint = 0, i = 0; i < FOLDCMP[icmp] - 1; i++)
		{
			for (j = i + 1; j < FOLDCMP[icmp]; j++, ipoint++)
			{
/*  
    Defining static shift for j-th and i-th traces
*/
			   static_j = to_be_calculated[i_to_calc][POINTER[icmp][j][0]] + to_be_calculated[i_to_calc][NSOURCES+POINTER[icmp][j][1]];
			   static_i = to_be_calculated[i_to_calc][POINTER[icmp][i][0]] + to_be_calculated[i_to_calc][NSOURCES+POINTER[icmp][i][1]];

			   shift = TOTAL_LAG/2 + NINT(static_j - static_i);
/*
    That's the lag. note the proper shift.
*/
			   eval_returned[i_to_calc] -= (double) CROSS[ipoint][shift];
			}
		}
           }
	}
	return;
}

/**** this routine will provide a conj. grad. search	       */
/**** for each subpopulation         		   Denver Basin*/

int Go_uphill(best, OBJ_PAST, genes)
/*
     element best will be a initial guess for a conjugate gradient search
*/
double **best; 		/* floating point representation                */
double *OBJ_PAST;	/* initial performances				*/
int genes;		/* # of elements in vect			*/
{
	static int firstflag = 1;
	double **GRADIENT_BEF;			/* the previous iteration */
	double **SEARCH;			/* search direction	*/
	double **GRADIENT; 			/* obvious...		*/
	double **NEW;				/* NEW: current solution*/
	double *OBJ_CUR;			/* current performance  */
	double aux;				/* auxiliar		*/
	double *temp;				/* swapping		*/
	double *descent;			/* to reset SEARCH 	*/
	double *BETA;				/* steps		*/
	double change;				/* relat change in fitness */
	double aux1, aux2;			/* auxiliary quantities */
	register int i, j, jj;			/* counters		*/
	int how_many, i_to_calc;		/* for tape management	*/
	int *DONE;				/* stopping flag per member */
	int *ITERATIONS;			/* iterations / member */
	int DONE_TOTAL=0;			/* overall stopping flag */
	int number_of_iterations=0;		/* # of iter. performed */
	int first_iter=1;			/* first iterat. flag	*/
	int ipoint_base, iparam;
	int itrace, icmp;        		/* pointers             */
/*
    Setting the pointers for the gradient calculation 
*/
	if (firstflag) 
	{
/*
    Allocating memory for DERIVATIVE, and POINT_CROSS
*/
                DERIVATIVE = alloc2int(NSOURCES + NRECEIVERS, NCMP);
                if (DERIVATIVE == NULL) Error("Allocation failed for DERIVATIVE");
		POINT_CROSS = alloc3int(MAXFOLD, NSOURCES + NRECEIVERS, NCMP);
                if (POINT_CROSS == NULL) Error("Allocation failed for POINT_CROSS");
		for (icmp = 0; icmp < NCMP; icmp++)
		{
			for (iparam = 0; iparam < NSOURCES + NRECEIVERS; iparam++)
			{
				itrace = 0;
				if (iparam < NSOURCES)	/* source statics */
				{
					while(itrace < FOLDCMP[icmp] && POINTER[icmp][itrace][0] != iparam)
						itrace++;
					if (itrace == FOLDCMP[icmp])
						DERIVATIVE[icmp][iparam] = -1; 
/*
    In the above situation the parameter IPARAM does not show up in any one
    of the relative shifts in this CMP
*/
					else
						DERIVATIVE[icmp][iparam] = itrace;
				}
				else			/* receiver statics */
				{
					while(itrace < FOLDCMP[icmp] && (POINTER[icmp][itrace][1] + NSOURCES) != iparam)
						itrace++;
					if (itrace == FOLDCMP[icmp])
						DERIVATIVE[icmp][iparam] = -1; 
/*
    In the above situation the parameter IPARAM does not show up in any one
    of the relative shifts in this CMP
*/
					else 
						DERIVATIVE[icmp][iparam] = itrace;
				}

				if (DERIVATIVE[icmp][iparam] != -1)
				{
					for (j = 0, ipoint_base = 0; j < itrace; j++)
	/* this will be used later */		ipoint_base += FOLDCMP[icmp] - 1 - j;

/*
    Computing pointers to use in the derivatives of the Xcorrelation. These
    pointers will able the search of the CORRECT Xcorrelation to be 
    differentiated
*/
					for (i = 0, POINT_CROSS[icmp][iparam][i] = 0; i < FOLDCMP[icmp]; i++, POINT_CROSS[icmp][iparam][i] = 0)
					{
						if (i < itrace)
						{
						/* computing POINT_CROSS */
							for (j = 0; j < i; j++)
								POINT_CROSS[icmp][iparam][i] += FOLDCMP[icmp] - 1 - j;
							POINT_CROSS[icmp][iparam][i] += itrace - i - 1;

					/* POINT_CROSS is the pointer for */
					/* computing the derivative of the */
   					/* cross_correlation 		   */
						}
						else if (i > itrace)
							POINT_CROSS[icmp][iparam][i] = ipoint_base + i - itrace - 1;
					}	/* derivative exists */
				}	/* all traces of CMP icmp */
			}	/* all parameters */
		}	/* for all cmps */
/*
    with the information computed above one should be able to
    calculate the derivatives 
*/
		firstflag = 0;		/* end of initialization	   */
	}
/*
    Necessary memory allocation
*/
        NEW = alloc2double(NSOURCES + NRECEIVERS, Popsize);
        if (NEW == NULL) Error("Allocation failed for NEW");

        GRADIENT = alloc2double(NSOURCES + NRECEIVERS, Popsize);
        if (GRADIENT == NULL) Error("Allocation failed for GRADIENT");

        GRADIENT_BEF = alloc2double(NSOURCES + NRECEIVERS, Popsize);
        if (GRADIENT_BEF == NULL) Error("Allocation failed for GRADIENT_BEF");

        SEARCH = alloc2double(NSOURCES + NRECEIVERS, Popsize);
        if (SEARCH == NULL) Error("Allocation failed for SEARCH");

        OBJ_CUR = alloc1double(Popsize);
        if (OBJ_CUR == NULL) Error("Allocation failed for OBJ_CUR");

        BETA = alloc1double(Popsize);
        if (BETA == NULL) Error("Allocation failed for BETA");

        descent = alloc1double(Popsize);
        if (descent == NULL) Error("Allocation failed for descent");

        DONE = alloc1int(Popsize);
        if (DONE == NULL) Error("Allocation failed for DONE");

        ITERATIONS = alloc1int(Popsize);
        if (ITERATIONS == NULL) Error("Allocation failed for ITERATIONS");
/*
    necessary resetting
*/
	for (i_to_calc = 0; i_to_calc < Popsize; i_to_calc++)
	{
		DONE[i_to_calc] = 0;
		BETA[i_to_calc] = 0.;
		ITERATIONS[i_to_calc] = 0;

        	for (iparam = 0; iparam < NSOURCES + NRECEIVERS; iparam++)
        	{       /* that's reseting  */
                	GRADIENT[i_to_calc][iparam] = 0.;
                	GRADIENT_BEF[i_to_calc][iparam] = 0.;
                	SEARCH[i_to_calc][iparam] = 0.;
        	}
	}
/*
    the gradient will be calculated now for the 1st iteration
*/
	Derivatives(DONE, GRADIENT, best);
/*
    BETA, the step to calculate the CONJUGATE GRADIENT directions
*/
	/* DONE and DONE_TOTAL control the iterations */
	while (!DONE_TOTAL)
	{
		/* search directions */
		for (i_to_calc = 0; i_to_calc < Popsize; i_to_calc++)
		{
			if (!DONE[i_to_calc])
			{
				for (aux1=0., aux2=0., iparam = 0; iparam < NSOURCES + NRECEIVERS; iparam++)
				{	
					aux1 += GRADIENT[i_to_calc][iparam]*(GRADIENT[i_to_calc][iparam]-GRADIENT_BEF[i_to_calc][iparam]);
					aux2 += GRADIENT_BEF[i_to_calc][iparam]*GRADIENT_BEF[i_to_calc][iparam];
				}
				BETA[i_to_calc] = aux1 / aux2;
/* 
    Computing the new search direction 
*/
				for (descent[i_to_calc] = 0., iparam = 0; iparam < NSOURCES + NRECEIVERS; iparam++)
				{
					SEARCH[i_to_calc][iparam] = -GRADIENT[i_to_calc][iparam] + BETA[i_to_calc]*SEARCH[i_to_calc][iparam];
					descent[i_to_calc] += GRADIENT[i_to_calc][iparam] * SEARCH[i_to_calc][iparam];
				}
				if (descent[i_to_calc] > 0. || first_iter)
				{	/* SEARCH = -GRADIENT */
					/* always true if 1st iter */
					for (descent[i_to_calc] = 0., iparam = 0; iparam < NSOURCES + NRECEIVERS; iparam++)
					{
						SEARCH[i_to_calc][iparam] = -GRADIENT[i_to_calc][iparam];
						descent[i_to_calc] += GRADIENT[i_to_calc][iparam] * SEARCH[i_to_calc][iparam];
					}
				}
                        }       /* closing the if for done */
                }       /* closing the i_to_calc loop */
/* 
	Doing line search
*/