📄 ga_real_imp.cpp
字号:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <conio.h>
#include <ctype.h>
#define MSAMP 2000
#define MNANG 10
#define MNSAM MNANG*MSAMP
#define MNWLT 200
#define MAXRAND 32768
void ang1_ARAT(double* ar,double zpp,double vpp,double zss,double vss,int m,int nctrl,int kps);
void reflect(double ref[MSAMP],double Vp[MSAMP],double Vs[MSAMP],double Den[MSAMP],double angle,double Vp0,double Vs0,double Den0,int nt);
void readdata(double seis[MNSAM],double Time[MSAMP],double Vp[MSAMP],double Vs[MSAMP],double Den[MSAMP],double wavlet[MNWLT],double ang[MNANG],int nt[1],int nang[1],int nwlt[1], double Vp0[1],double Vs0[1],double Den0[1],double search_range[1],int inter_number[1]);
double cc(int n,int m);
void readdata(double seis[MNSAM],double Time[MSAMP],double Vp[MSAMP],double Vs[MSAMP],double Den[MSAMP],double wavlet[MNWLT],double ang[MNANG],int nt[1],int nang[1],int nwlt[1], double Vp0[1],double Vs0[1],double Den0[1],double search_range[1],int inter_number[1])
{
FILE *fp1,*fp2,*fp3;
// char Log_Title[120];
double temp1;
int i;
fp1=fopen("ga_para21470.dat","r");
fp2=fopen("ga_log21470.dat","r");
fp3=fopen("ga_seis21470_norm.dat","r");
fscanf(fp1,"%d%d%d%lf%lf%lf%lf%d",&nt[0],&nang[0],&nwlt[0], &Vp0[0],&Vs0[0], &Den0[0],&search_range[0],&inter_number[0]);
for (i=0;i<nang[0];i++)
{
fscanf(fp1,"%lf",&ang[i]);
}
for (i=0;i<nwlt[0];i++)
{
fscanf(fp1,"%lf",&wavlet[i]);
wavlet[i]=wavlet[i];
}
for (i=0;i<nt[0];i++)
{
fscanf(fp2,"%lf %lf %lf %lf",&Time[i],&Vp[i],&Vs[i],&Den[i]);
}
// for (i=0;i<nt[0];i++)
// {
// fscanf(fp3,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",&temp1,&temp2,&temp3,&seis[i],&seis[i+nt[0]],&seis[i+2*nt[0]],&seis[i+3*nt[0]],&seis[i+4*nt[0]],&seis[i+5*nt[0]],&seis[i+6*nt[0]]);
// }
for (i=0;i<nt[0]*nang[0];i++)
fscanf(fp3,"%lf %lf", &temp1,&seis[i]);
}
void reflect(double ref[MSAMP],double Vp[MSAMP],double Vs[MSAMP],double Den[MSAMP],double angle,double Vp0,double Vs0,double Den0,int nt)
{
double Vp1,Vp2,Vs1,Vs2,b1,b2;
double zpp,vpp,zss,vss;
double aa,a2,a4,a6;
double ar[11];
int i ;
//
for(i=0;i<MSAMP;i++)
ref[i]=0.0;
//
for(i=0;i<nt;i++){
if(i>0){
Vp1=Vp[i-1];
Vs1=Vs[i-1];
b1=Den[i-1];
}
else{
Vp1=Vp0;
Vs1=Vs0;
b1=Den0;
}
Vp2=Vp[i];
Vs2=Vs[i];
b2=Den[i];
//
zpp=b1*Vp1/(b2*Vp2);
vpp=Vp1/Vp2;
zss=b1*Vs1/(b2*Vs2);
vss=Vs2/Vp2;
//Zoeppritz
ang1_ARAT(ar,zpp,vpp,zss,vss,3,0,0);
aa=sin(angle);
a2=aa*aa;
a4=a2*a2;
a6=a2*a4;
ref[i]=ar[1]+ar[3]*a2+ar[5]*a4+ar[7]*a6;
}
}
//*****************************************
//Ο撞fold悼訁
//*****************************************
void fold(double ref[MSAMP],double wavlet[MNWLT],double seis0[MSAMP],int nt, int nwlt)
{
int i,nn,k,j;
nn=nt+nwlt-1;
for(i=0;i<nn;i++)
seis0[i]=0.0;
for(i=0;i<nt;i++)
for(j=0;j<nwlt;j++){
k=i+j;
seis0[k]=ref[i]*wavlet[j]+seis0[k];
}
for(i=0;i<nt;i++)
seis0[i]=seis0[i+nwlt/2];
}
//*****************************************
//夭覝Ο撞
//*****************************************
void forward(double syn[MNSAM],double Vp[MSAMP],double Vs[MSAMP],double Den[MSAMP],double ang[MNANG],double wavlet[MNWLT],int nt,int nang,int nwlt,double Vp0,double Vs0,double Den0)
{
double angle, seis0[MSAMP], ref[MSAMP];
int i,j;
for(i=0;i<nang;i++){
angle=ang[i]*3.1415926/180.0;
reflect(ref,Vp,Vs,Den,angle,Vp0,Vs0,Den0,nt);
fold(ref,wavlet,seis0,nt,nwlt);
for(j=0;j<nt;j++){
int j1;
j1=j+i*nt;
syn[j1]=seis0[j];//seis0 from fold()
}
}
}
//*****************************************
//Ο撞GA_waveinvs悼訁
//*****************************************
#define Samples 31
void GA_waveinvs(double seis[MNSAM],double Vp[MSAMP],double Vs[MSAMP],double Den[MSAMP],double wavlet[MNWLT],double ang[MNANG],int nt,int nang,int nwlt, double Vp0,double Vs0,double Den0,double search_range,int inter_number)
{
double Vpr[MSAMP],Vsr[MSAMP],Denr[MSAMP]; //store bind referance
double syn[MNSAM]; //forward results
double e1=0.005;
// unsigned long n;
// unsigned long i,j,k,kkk;
int i,j,kkk,n;
//The Range of Parameters Vp, Vs, Den.
double Vp_Range[MSAMP],Vs_Range[MSAMP],Den_Range[MSAMP];
double D_Vp[MSAMP],D_Vs[MSAMP],D_Den[MSAMP];
//Control possiblility 佬磿負
// double Repr_p=1.0,Cross_p=0.45,Muta_p=0.45;
double mut_p;
//Coding Parameters
// unsigned long CodeVp1,CodeVs1,CodeDen1,CodeVp2,CodeVs2,CodeDen2;
unsigned long MaxCodeVp=256,MaxCodeVs=256,MaxCodeDen=256;
unsigned long M_Bite=8;
//Initial Samples
unsigned long Min_Sample_Number=-2,Max_Sample_Number=-1;
//Object Parameters
double Vp_Object[MSAMP],Vs_Object[MSAMP],Den_Object[MSAMP];
// unsigned long Code1,Code2,Code3;
// unsigned long Loca1;
// double pp,pp1,pp2,pp3;
// Objects Functions
double Objects[200];
double Min_Objects0,Min_Objects=1.0e20,sum11=0;
double Max_Objects=1.0e-20;
double Objects_var=0,Objects_mean=0;
// Fitness Functions
// double Fitness[201],Fitness_p[201];
//
double Ge_Vp[MSAMP][Samples],Ge_Vs[MSAMP][Samples],Ge_Den[MSAMP][Samples];
double Ge_Vp0[MSAMP][Samples],Ge_Vs0[MSAMP][Samples],Ge_Den0[MSAMP][Samples];
double Ge_Vp1[MSAMP][Samples/2],Ge_Vs1[MSAMP][Samples/2],Ge_Den1[MSAMP][Samples/2];
double Ge_VpM[MSAMP],Ge_VsM[MSAMP],Ge_DenM[MSAMP];
double irand;//irand1;
// int irand2;
int irand1,p;
n=nt*nang;
//search_range=0.125;
printf(" \n come here 3.1");
printf("n=%d,nt=%d,nang=%d",n,nt,nang);
getch();
// prepare for the ga ,
for(i=0;i<nt;i++)
{
// vpr vsr denr : reference parameters;
Vpr[i]=Vp[i];
Vsr[i]=Vs[i];
Denr[i]=Den[i];
// the search range, possible change or permitalbe range;
Vp_Range[i]=search_range*Vp[i];
Vs_Range[i]=search_range*Vs[i];
Den_Range[i]=search_range*Den[i];
// steps that ga used
D_Vp[i]=Vp_Range[i]/MaxCodeVp;
D_Vs[i]=Vs_Range[i]/MaxCodeVs;
D_Den[i]=Den_Range[i]/MaxCodeDen;
}
// the above: Calculate the interval of parameters
//the following : produce the initial population(16 samples, 16 persons)
// samples: individual number of the group and nt the number of the elements of one individual;
/// in fact, there are three groups in this GA algorithm;
// one: Ge_Vp, 2: Ge_Vs, 3: Ge_Den;
// the following program is right,
// randomly produce all the individuals
for(i=0;i<Samples;i++)
{
for(j=0;j<nt;j++)
{
Ge_Vp[j][i]=Vpr[j]-Vp_Range[j]/2.0+rand()%MaxCodeVp*D_Vp[j];
Ge_Vs[j][i]=Vsr[j]-Vs_Range[j]/2.0+rand()%MaxCodeVs*D_Vs[j];
Ge_Den[j][i]=Denr[j]-Den_Range[j]/2.0+rand()%MaxCodeDen*D_Den[j];
}
}
// Object Functions
for(i=0;i<n;i++)
sum11=sum11+fabs(seis[i])/n;
// individual circle
for(i=0;i<Samples;i++)
{
for(j=0;j<nt;j++)
{
Vp[j]=Ge_Vp[j][i];
Vs[j]=Ge_Vs[j][i];
Den[j]=Ge_Den[j][i];
}
// for every individual, calculate it's fitness
forward(syn,Vp,Vs,Den,ang,wavlet,nt,nang,nwlt,Vp0,Vs0,Den0);
// object is related to fitness calculation
Objects[i]=0;
//
for(j=0;j<n;j++)
{
Objects[i]=Objects[i]+fabs(seis[j]-syn[j])/n;
}
// 每一个个体偏离进归一化；
Objects[i]=Objects[i]/sum11;
// 所有偏差之和
Objects_mean=Objects_mean+Objects[i];
// 偏差的最小值计算，并记录下号码；
// search for the best individual
if(Objects[i]<Min_Objects)
{
Min_Objects=Objects[i];
Min_Sample_Number=i;
}
// printf("Objects=%f\n",Objects[i]);
}
printf(" \n come here 3.3\n\n");
//getch();
// 全部个体偏差的平均值
Objects_mean=Objects_mean/Samples;
printf("Min_Objects=%f\n",Min_Objects);
// 如果最优的个体满足，则输出结果，
if(Min_Objects<e1)
{// Output results
for(j=0;j<nt;j++)
{
Vp[j]=Ge_Vp[j][Min_Sample_Number];
Vs[j]=Ge_Vs[j][Min_Sample_Number];
Den[j]=Ge_Den[j][Min_Sample_Number];
}
// output the result
return;
}
// Min_Object s 0 上下代之间的关系；
Min_Objects0=Min_Objects;
// Ge_*记录上一代； Ge_*0 下一代
// Begin to the GA
//Dai xunhuan
for(kkk=0;kkk<inter_number;kkk++)
{
//***********最优个体位于最后面Ge_*[][Sample-1]及复制操作********
for(j=0;j<nt;j++)
{
//复制
Ge_Vp0[j][0]=Ge_Vp[j][Min_Sample_Number];
Ge_Vs0[j][0]=Ge_Vs[j][Min_Sample_Number];
Ge_Den0[j][0]=Ge_Den[j][Min_Sample_Number];
/*
Vp_Object[j]= Ge_Vp0[j][0];
Vs_Object[j]=Ge_Vs0[j][0];
Den_Object[j]=Ge_Den0[j][0];
*/
//交换
Ge_VpM[j]=Ge_Vp[j][Min_Sample_Number];
Ge_VsM[j]=Ge_Vs[j][Min_Sample_Number];
Ge_DenM[j]=Ge_Den[j][Min_Sample_Number];
Ge_Vp[j][Min_Sample_Number]=Ge_Vp[j][Samples-1];
Ge_Vs[j][Min_Sample_Number]=Ge_Vs[j][Samples-1];
Ge_Den[j][Min_Sample_Number]=Ge_Den[j][Samples-1];
Ge_Vp[j][Samples-1]=Ge_VpM[j];
Ge_Vs[j][Samples-1]=Ge_VsM[j];
Ge_Den[j][Samples-1]=Ge_DenM[j];
}
//*****************************************
//******************个体之间相互竞争选择优秀个体*******
p=0;
for(i=1;i<=Samples-2;i=i+2)
{
if(Objects[i-1]<=Objects[i])
{ for(j=0;j<nt;j++)
{
Ge_Vp1[j][p]= Ge_Vp[j][i-1];
Ge_Vs1[j][p]= Ge_Vs[j][i-1];
Ge_Den1[j][p]= Ge_Den[j][i-1];
}
p=p+1;
}
else if(Objects[i-1]>Objects[i])
{ for(j=0;j<nt;j++)
{
Ge_Vp1[j][p]= Ge_Vp[j][i];
Ge_Vs1[j][p]= Ge_Vs[j][i];
Ge_Den1[j][p]= Ge_Den[j][i];
}
p=p+1;
}
}
//****************Crossover****************************
for(i=0;i<=Samples/2-1;i++)
{
for(j=0;j<nt;j++)
{
//parameter Vp
irand=rand()/(32768.0);
irand1=(int)rand()/(32768.)*Samples/2;
Ge_Vp0[j][i+1]=irand*Ge_Vp1[j][irand1]+(1-irand)*Ge_Vp1[j][i];
//parameter Vs
irand=rand()/(32768.0);
irand1=(int)rand()/(32768.)*Samples/2;
Ge_Vs0[j][i+1]=irand*Ge_Vs1[j][irand1]+(1-irand)*Ge_Vs1[j][i];
//parameter Density
irand=rand()/(32768.0);
irand1=(int)rand()/(32768.)*Samples/2;
Ge_Den0[j][i+1]=irand*Ge_Den1[j][irand1]+(1-irand)*Ge_Den1[j][i];
}
}
//*****************Mutation****************************
for(i=Samples/2+1;i<Samples;i++)
{
for(j=0;j<nt;j++)
{
if(Min_Objects>0.9)
mut_p=12;
else if(Min_Objects<=0.9&&Min_Objects>0.7)
mut_p=30*Min_Objects-15;
else if(Min_Objects<=0.7&&Min_Objects>0.6)
mut_p=40*Min_Objects-22;
//parameter Vp
Ge_Vp0[j][i]=(Ge_Vp0[j][0])+(rand()/(32768.0)*mut_p-mut_p/2)*D_Vp[j];
//parameter Vs
Ge_Vs0[j][i]=(Ge_Vs0[j][0])+(rand()/(32768.0)*mut_p-mut_p/2)*D_Vs[j];
//parameter Density
Ge_Den0[j][i]=(Ge_Den0[j][0])+(rand()/(32768.0)*mut_p-mut_p/2)*D_Den[j];
}
}
Min_Objects=1.0e20;
Objects_mean=0.0;
for(i=0;i<Samples;i++)
{
for(j=0;j<nt;j++)
{
Vp[j]=Ge_Vp0[j][i];
Vs[j]=Ge_Vs0[j][i];
Den[j]=Ge_Den0[j][i];
}
forward(syn,Vp,Vs,Den,ang,wavlet,nt,nang,nwlt,Vp0,Vs0,Den0);
Objects[i]=0;
for(j=0;j<n;j++){
Objects[i]=Objects[i]+fabs(seis[j]-syn[j])/n;
}
Objects[i]=Objects[i]/sum11;
Objects_mean=Objects_mean+Objects[i];
if(Objects[i]<Min_Objects)
{
Min_Objects=Objects[i];
Min_Sample_Number=i;
}
if(Objects[i]>Max_Objects)
{
Max_Objects=Objects[i];
Max_Sample_Number=i;
}
// printf("Objects=%f\n",Objects[i]);
}
Objects_mean=Objects_mean/Samples;
printf("kkk=%d,Min_Objects=%f,Min_Objects0=%f\n",kkk,Min_Objects,Min_Objects0);
if(Min_Objects<e1)
{// Output results
for(j=0;j<nt;j++)
{
Vp[j]=Ge_Vp[j][Min_Sample_Number];
Vs[j]=Ge_Vs[j][Min_Sample_Number];
Den[j]=Ge_Den[j][Min_Sample_Number];
}
return;
}
if(Min_Objects<Min_Objects0)
{
Min_Objects0=Min_Objects;
for(j=0;j<nt;j++)
{
Vp_Object[j]=Ge_Vp[j][Min_Sample_Number];
Vs_Object[j]=Ge_Vs[j][Min_Sample_Number];
Den_Object[j]=Ge_Den[j][Min_Sample_Number];
}
}
/*
else
{
Objects[Max_Sample_Number]=Min_Objects0;
for(j=0;j<nt;j++)
{
Ge_Vp[j][Max_Sample_Number]=Vp_Object[j];
Ge_Vs[j][Max_Sample_Number]=Vs_Object[j];
Ge_Den[j][Max_Sample_Number]=Den_Object[j];
}
}
*/
for(i=0;i<Samples;i++)
{
for(j=0;j<nt;j++)
{
Ge_Vp[j][i]=Ge_Vp0[j][i];
Ge_Vs[j][i]=Ge_Vs0[j][i];
Ge_Den[j][i]=Ge_Den0[j][i];
}
}
}//end of kkk
// in two cases, output vp,vs,den,
for(j=0;j<nt;j++){
Vp[j]=Vp_Object[j];
Vs[j]=Vs_Object[j];
Den[j]=Den_Object[j];
}
}
//************************************
//main function
//************************************
int main(int argc, char* argv[])
{
int j;
double Vp[MSAMP],Vs[MSAMP],Den[MSAMP],Time[MSAMP];//common vpsdn
double seis[MNSAM],syn[MNSAM]; //forward results
double ang[MNANG],wavlet[MNWLT];
double Vp0,Vs0,Den0,search_range;
int nt,nang,nwlt,inter_number;//nt:samples;nang:angles;nwlt: wavelets;
double VVp0[1],VVs0[1],DDen0[1],Ssearch_range[1];
int Nnt[1],Nnang[1],Nnwlt[1],Iinter_number[1];
FILE *fp1,*fp2;
printf(" \n come here 1");
// getch();
readdata(seis,Time,Vp,Vs,Den,wavlet,ang,Nnt,Nnang,Nnwlt,VVp0,VVs0,DDen0,Ssearch_range,Iinter_number);
printf(" \n come here 2");
//getch();
nt=Nnt[0];
nang=Nnang[0];
nwlt=Nnwlt[0];
Vp0=VVp0[0];
Vs0=VVs0[0];
Den0=DDen0[0];
search_range=Ssearch_range[0];
inter_number=Iinter_number[0];
printf(" \n come here 3");
getch();
GA_waveinvs(seis,Vp,Vs,Den,wavlet,ang,nt,nang,nwlt,Vp0,Vs0,Den0,search_range,inter_number);
printf(" \n come here 4");
//getch();
fp1=fopen("W21470_res.dat","w");//讋⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -