/*           功能：模擬退火反演算法                      */
/*           文獻：地球物理反演                          */
/*           作者：王家映教授                            */
/*           出版：高等教育出版社(第二版)                */
/*           編程：蔣龍聰                                */
/*           單位：中國地質大學(武漢)地球物理與空間信息  */
/*           專業：地球探測與信息技術                    */

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>

/*           功能：模擬退火反演算法                      */
#define DS     0.99
#define ERRO   0.50
#define Lambda 6
#define LayerN 4
#define H      5
#define MarkT  10
#define MinT   0.001
#define POINTS 1001
#define SF     0.5

void   SA(double,double);
double ranD(void);
double obfun(double[]);
double OBS[POINTS];
double Temperature;
FILE * SATXT;

double ranD(void)
{
	double ranD=0.0;
	ranD=(rand()%5000)/5000.0;
//	ranD=0.1+(0.6-0.2)*ranD;
	return(ranD);
	
}

double obfun(double x[])
{
	double g[POINTS],b[201],r[POINTS];//定義數組
	double dt=0.5,f0=60;  //數組賦值   dt=0.5ms
	double t[3];//時間
	double v[4];
	double SUM=0.0;
    int mx=1;  
    int i=0,j=0,nw=100;
    double h[3]={50,20,50};   


	v[0]=x[0];
	
	v[1]=x[1];
	
	v[2]=x[2];

	v[3]=x[3];
   
	
	//初始化反射系數數組		  
	for(i=0;i<POINTS;i++)
		r[i]=0.0;

	t[0]=2*h[0]/v[0];
	j=(int)(1000*t[0]/dt);
	r[j]=(v[1]-v[0])/(v[1]+v[0]); //1層 

	
	for(i=1;i<MarkT;i++)
	{
		t[i]=t[i-1]+2*h[i]/v[i];
		j=(int)(1000*t[i]/dt);
		r[j]=(v[i+1]-v[i])/(v[i+1]+v[i]);//2-4層
	}	

	
	for(i=-nw;i<nw+1;i++)
	{
		double a=(0.001*3.1415926*f0*i*dt);
        b[i+nw]=100.0*(1.0-2.0*a*a)*exp(-a*a);//求取子波，0相位，
        
	}

	
	//Convolution   子波和參數做卷積
	for(i=0;i<POINTS;i++)//從第一道循環
	{
		double sum=0.0;
	    for(j=-nw;j<nw+1;j++)
		{
			if(i-j>=0&&i-j<POINTS)
		    sum=sum+b[j+nw]*r[i-j];
		}
		g[i]=sum;
		
	}

	for(i=0;i<POINTS;i++)         //目標函數
		SUM+=(OBS[i]-g[i])*(OBS[i]-g[i]);
	    SUM=sqrt(SUM);
	return(SUM);

}

/*采用保留局部最優解方式*/
void SA(double lbounds[],double ubounds[])
{
	int MaxGen=0;
	int Gen=0;
	int Counter=0;
	int i=0,j=0,k=0;
	double R=0.0;
	double E0=0.0;
	double E1=0.0;
	double DeltaE=0.0;
	double P=0.0;
	double T=0.0;
	double Val=0.0;
	double BestVal=0.0;
	double Local[LayerN];
	double Global[LayerN];
	double CurSol[LayerN];
	double NexSol[LayerN];
	
	T=Temperature;
	MaxGen=(int)(log(MinT/Temperature)/log(DS)+0.5);
	MaxGen=MaxGen;

	//第一步:隨機選擇初始值

	for(i=0;i<(LayerN);i++)
	{
		CurSol[i]=(ubounds[i])*ranD();
//		cout<<CurSol[i];
	}
	
/*
  	for(i=0;i<LayerN;i++)
	{
		CurSol[i]=500+i*200;
	}
*/
	//初始化局部變量和全局變量值
	for(i=0;i<(LayerN);i++)
	{
		Local[i]=CurSol[i];
		Global[i]=CurSol[i];
	}

	//第二步:在某一溫度下,對當前模型進行擾動,得到新的模型值NextSol
	while(T>MinT)
	{
		
		Gen++;
		for(i=0;i<MarkT;i++)
		{
			
			Counter++;
			R=Gen/(2*MaxGen);

			if( (int) ( ranD()+0.5 )==0 )
			{
				for(i=0;i<LayerN;i++)
					NexSol[i]=CurSol[i]+(ubounds[i]-CurSol[i])*( pow( ranD()*(1.0-R),Lambda) );
			}
			else
			{
				for(i=0;i<LayerN;i++)
					NexSol[i]=CurSol[i]-(CurSol[i]-lbounds[i])*( pow( ranD()*(1.0-R),Lambda) );

			}

/*
			if( (int) ( ranD()+0.5 )==0 )
			{
				for(i=0;i<LayerN;i++)
					NexSol[i]=CurSol[i]+(ubounds[i]-lbounds[i])*( pow( ranD()*(1.0-R),Lambda) );
			}
			else
			{
				for(i=0;i<LayerN;i++)
					NexSol[i]=CurSol[i]-(ubounds[i]-lbounds[i])*( pow( ranD()*(1.0-R),Lambda) );
			}
*/		
			//判斷是否為全局最優解
			if(obfun(NexSol)<obfun(Global))
			{
				for(i=0;i<(LayerN);i++)
				{
					Local[i]=Global[i];
					Global[i]=NexSol[i];
				}
			}

			E1=obfun(NexSol);
			E0=obfun(CurSol);
			DeltaE=E1-E0;
			P=exp(-DeltaE/T);
			P=pow((1-H*P*E1/(E0*E0)),1/H);
		//	cout<<"DeltaE"<<DeltaE<<endl;
			if(DeltaE<0)
			{
				for(i=0;i<(LayerN);i++)
					CurSol[i]=NexSol[i];
			}
			//else if//( exp(-DeltaE/T)<ranD() )
			else if(P<ranD())
			{
				for(i=0;i<(LayerN);i++)
					CurSol[i]=NexSol[i];
			}

			Val=obfun(NexSol);
			BestVal=obfun(Global);
			fprintf(SATXT,"%10.6f   %10.6f",Val,BestVal);
			

			/*               把最優解輸出到文件中去           */
			if(Val<=ERRO||BestVal<=ERRO)
			{
				printf("尋找解得次數為:%5d\n",Counter);
/*				
				if ((salog = fopen("galog.txt","w"))==NULL)
				{
					exit(1);
				}
*/
				fprintf(SATXT,"\n\n");
				for(i=0;i<(LayerN);i++)
					fprintf(SATXT," %10.7f   ",Global[i]);
//				fclose(salog);
				printf("Sucess\n");
				exit(0);
			}
			//第三步:降低溫度
			T=T*DS;
		//	simulated(BestVal);
//			cout<<" T= "<<T<<"  Val= "<<Val<<"  BestVal= "<<BestVal<<" Counter= "<<Counter <<endl;
		}
//		simulated(BestVal);
	}


}

void main(void)
{	
	int i=0,j=0,k=0;
	double tim=0.0,E1=0,E0=0,DeltaE=0;
	double amptitude=0.0;
	double Time[POINTS];
	double CurSol[LayerN],NexSol[LayerN];
	double lbound[LayerN]={1000,800,1000,1000};
	double ubound[LayerN]={3000,1500,3000,3000};

	srand((unsigned int)time(NULL));
/*           從文件中讀取數據                      */

	FILE *INFILE;

	if ((INFILE = fopen("result.txt","r"))==NULL)
	{
		printf("打開的文件不存在,請調試!\n");
		exit(1);
	}
	for (i=0; i<POINTS; i++)
	{
		fscanf(INFILE, "%lf",&amptitude);
        fscanf(INFILE, "%lf",&tim);
		Time[i]=tim;
		OBS[i]=amptitude;	
	}
	for (i=0; i<POINTS; i++)
		printf("F=%f  Data=%f\n",Time[i],OBS[i]);
	
	fclose(INFILE);

	SATXT=fopen("simulate.txt","w");
/*           確定初始溫度值                      */	
	for(i=0;i<(LayerN);i++)
	{
		CurSol[i]=(ubound[i])*ranD();
		fprintf(SATXT,"  %10.7f   ",CurSol[i]);
	}

//	fprintf(salog," \n   ");

	for(i=0;i<LayerN;i++)
	{
		NexSol[i]=CurSol[i]+(2*ranD()-1)*SF;
	}

	E1=obfun(NexSol);
	E0=obfun(CurSol);
	DeltaE=E1-E0;
	Temperature=-DeltaE/log(0.995);
	printf("初始溫度為：%f\n",Temperature);

	SA(lbound,ubound);
   
	fclose(SATXT);
	


}