#include "sgm_esl.h"

int SGM_ESL::Routine(char* parafile){
	double object_func, delta_grav;
	string temp_str;
	//初始化char数组
	char run_type[1024], data_type[1024], msh_name[1024], msh_out_name[1024], obs_name[1024], obsout_name[1024];
	strcpy(run_type,"NULL");
	strcpy(data_type,"NULL");
	strcpy(msh_name,"NULL");
	strcpy(msh_out_name,"NULL");
	strcpy(obs_name,"NULL");
	strcpy(obsout_name,"NULL");
	strcpy(init_model_name_,"NULL");
	//初始化参数
	int regional_tag = -1;
	double layer_depth = 50, layer_Depth = 100;
	alpha_ = 1.0;
	beta_ = 1.0;
	iter_times_ = 300;
	I0_ = loc_I_ = reg_I_ = 90;
	D0_ = loc_D_ = reg_D_ = 0;

	ifstream infile;
	if (open_infile(infile,parafile)) return -1;
	while(getline(infile,temp_str)){
		if (*(temp_str.begin()) == '#') continue;
		else{
			if (1 == sscanf(temp_str.c_str(),"run-type=%[^\n]",run_type))
				clog << "run type:\t" << run_type << endl;
			else if (1 == sscanf(temp_str.c_str(),"data-type=%[^\n]",data_type))
				clog << "data type:\t" << data_type << endl;
			else if (1==sscanf(temp_str.c_str(),"mshin-file=%[^\n]",msh_name))
				clog << "mshin file:\t" << msh_name << endl;
			else if (1==sscanf(temp_str.c_str(),"mshout-file=%[^\n]",msh_out_name))
				clog << "mshout file:\t" << msh_out_name << endl;
			else if (1==sscanf(temp_str.c_str(),"obsin-file=%[^\n]",obs_name))
				clog << "obsin file:\t" << obs_name << endl;
			else if (1==sscanf(temp_str.c_str(),"obsout-file=%[^\n]",obsout_name))
				clog << "obsout file:\t" << obsout_name << endl;
			else if (1==sscanf(temp_str.c_str(),"initial-model=%[^\n]",init_model_name_))
				clog << "initial model:\t" << init_model_name_ << endl;
			else if (1==sscanf(temp_str.c_str(),"depth-weighting=%lf",&beta_))
				clog << "depth weighting:\t" << beta_ << endl;
			else if (1==sscanf(temp_str.c_str(),"volume-weighting=%lf",&alpha_))
				clog << "volume weighting:\t" << alpha_ << endl;
			else if (1==sscanf(temp_str.c_str(),"regional-tag=%d",&regional_tag))
				clog << "regional tag:\t" << regional_tag << endl;
			else if (2==sscanf(temp_str.c_str(),"layer-depth=%lf/%lf",&layer_depth,&layer_Depth))
				clog << "layer depth:\t" << layer_depth << "-->" << layer_Depth << endl;
			else if (6==sscanf(temp_str.c_str(),"magnetic-para=%lf/%lf/%lf/%lf/%lf/%lf",&I0_,&D0_,&loc_I_,&loc_D_,&reg_I_,&reg_D_))
				clog << "incident/decline angle of geomagnetic field:\t" << I0_ << "/" << D0_ << endl 
				<< "incident/decline angle of magnetization (local):\t" << loc_I_ << "/" << loc_D_ << endl 
				<< "incident/decline angle of magnetization (regional):\t" << reg_I_ << "/" << reg_D_ << endl;
			else if (1==sscanf(temp_str.c_str(),"iteration-para=%d",&iter_times_))
				clog << "iteration times: " << iter_times_ << endl;
			else
			{
				cerr << BOLDRED << "error ==> " << RESET << "undefined parameter: " << temp_str << endl;
				return -1;
			}
		}
	}
	infile.close();
	clog << "================================" << endl;

	if (InitModBlocks(msh_name,regional_tag,layer_depth,layer_Depth)) return -1;
	if (InitObsPoints(obs_name)) return -1;

	if (!strcmp(run_type,"inversion")){
		//计算核函数
		if (!strcmp(data_type,"gravity")) CalGKernel();
		else if (!strcmp(data_type,"deltaT")) CalDeltaTKernel();
		else{
			cerr << BOLDRED << "error ==> " << RESET << "unknown data type: " << data_type << endl;
			return -1;
		} 
		//计算加权矩阵
		CalWvWz();
		CalWp(1.0);
		//首先按1.0计算数据拟合差矩阵
		CalWd(1.0);
		//计算total_feild_对应的PartB
		CalPartB(total_field_);
		//反演total_field_
		OptimizeCG_Precondition(total_field_);
		OutModBlocks(msh_out_name);
		OutModEleData(msh_out_name,"\"equivalent source\"",invert_model_);
		//输出分离的局部与区域场
		CalPartialModel(regional_field_);
		if (OutObsPoints(obsout_name, regional_field_)) return -1;
	}
	else if (!strcmp(run_type,"separate")){
		//计算核函数
		if (!strcmp(data_type,"gravity")) CalGKernel();
		else if (!strcmp(data_type,"deltaT")) CalDeltaTKernel();
		else{
			cerr << BOLDRED << "error ==> " << RESET << "unknown data type: " << data_type << endl;
			return -1;
		} 
		//计算加权矩阵
		CalWvWz();
		CalWp(1.0);
		//首先按1.0计算数据拟合差矩阵
		CalWd(1.0);
		//计算total_feild_对应的PartB
		CalPartB(total_field_);
		//反演total_field_
		OptimizeCG_Precondition(total_field_);
		//正演区域和局部异常
		CalPartialModel(regional_field_,regional_model_list_);
		CalPartialModel(object_field_,object_model_list_);

		//按0.5计算数据拟合差矩阵
		CalWd(0.5);
		//计算regional_feild_对应的PartB
		CalPartB(regional_field_);
		//反演regional_feild_
		OptimizeCG_Precondition(regional_field_);
		//正演regional_feild_区域和局部异常
		CalPartialModel(regional_regional_field_,regional_model_list_);
		CalPartialModel(regional_object_field_,object_model_list_);

		//计算object_feild_对应的PartB
		CalPartB(object_field_);
		//反演object_feild_
		OptimizeCG_Precondition(object_field_);
		//正演object_feild_区域和局部异常
		CalPartialModel(object_regional_field_,regional_model_list_);
		CalPartialModel(object_object_field_,object_model_list_);

		// 临时变量 用于收集每次迭代中区域与目标异常的变化
		/*
		int it_times = 0;
		string tmp_str, tmp_str2 = "inter_sep_ret_";
		stringstream tmp_ss;
		string inter_ret_str;

		// 临时代码 用于输出每次迭代中区域与目标异常
		tmp_ss.clear(); tmp_ss << it_times; tmp_ss >> tmp_str;
		inter_ret_str = tmp_str2 + tmp_str + ".txt";
		if (OutObsPoints(inter_ret_str.c_str(),regional_field_,object_field_,regional_object_field_)) return -1;
		it_times++;
		*/
		//计算目标函数
		object_func = 0.0;
		for (int i = 0; i < obs_points_num_; i++){
			object_func += pow((regional_object_field_[i] - object_regional_field_[i])/(sqrt(2)*obs_points_[i].dev),2);
		}
		object_func /= obs_points_num_;

		while(object_func > 1.0){
			clog << "Objective function = " << object_func << endl;
			//合并新的局部与区域场
			for (int i = 0; i < obs_points_num_; i++){
				regional_field_[i] = regional_regional_field_[i] + object_regional_field_[i];
				object_field_[i] = regional_object_field_[i] + object_object_field_[i];
			}

			//计算regional_feild_对应的PartB
			CalPartB(regional_field_);
			//反演regional_feild_
			OptimizeCG_Precondition(regional_field_);
			//正演regional_feild_区域和局部异常
			CalPartialModel(regional_regional_field_,regional_model_list_);
			CalPartialModel(regional_object_field_,object_model_list_);

			//计算object_feild_对应的PartB
			CalPartB(object_field_);
			//反演object_feild_
			OptimizeCG_Precondition(object_field_);
			//正演object_feild_区域和局部异常
			CalPartialModel(object_regional_field_,regional_model_list_);
			CalPartialModel(object_object_field_,object_model_list_);

			// 临时代码 用于输出每次迭代中区域与目标异常
			/*
			tmp_ss.clear(); tmp_ss << it_times; tmp_ss >> tmp_str;
			inter_ret_str = tmp_str2 + tmp_str + ".txt";
			if (OutObsPoints(inter_ret_str.c_str(),regional_field_,object_field_,regional_object_field_)) return -1;
			it_times++;
			*/

			//计算目标函数
			object_func = 0.0;
			for (int i = 0; i < obs_points_num_; i++){
				object_func += pow((regional_object_field_[i] - object_regional_field_[i])/(sqrt(2)*obs_points_[i].dev),2);
			}
			object_func /= obs_points_num_;
		}
		clog << "Objective function = " << object_func << endl;

		double total_engery = 0.0;
		double mean_engery = 0.0;
		for (int i = 0; i < obs_points_num_; i++){
			total_engery += pow(regional_object_field_[i],2);
			mean_engery += regional_object_field_[i];
		}
		clog << "Total energy of the remaining field = " << sqrt(total_engery) << endl;
		clog << "RMS of misfit (mGal|nT) = " << sqrt(total_engery/obs_points_num_) << endl;

		total_engery = 0.0;
		mean_engery /= obs_points_num_;
		for (int i = 0; i < obs_points_num_; i++){
			total_engery += pow(regional_object_field_[i] - mean_engery,2);
		}
		clog << "SD of misfit (mGal|nT) = " << sqrt(total_engery/obs_points_num_) << endl;

		//输出分离的局部与区域场
		if (OutObsPoints(obsout_name,regional_field_,object_field_,regional_object_field_)) return -1;
	}
	else if (!strcmp(run_type,"RTP")){
		//计算核函数
		if (!strcmp(data_type,"deltaT")) CalDeltaTKernel();
		else{
			cerr << BOLDRED << "error ==> " << RESET << "unsupported data type: " << data_type << endl;
			return -1;
		} 
		//计算加权矩阵
		CalWvWz();
		CalWp(1.0);
		//首先按1.0计算数据拟合差矩阵
		CalWd(1.0);
		//计算total_feild_对应的PartB
		CalPartB(total_field_);
		//反演total_field_
		OptimizeCG_Precondition(total_field_);
		//重新计算正演核矩阵 垂直极化的情况
		I0_ = loc_I_ = reg_I_ = 90.0;
		D0_ = loc_D_ = reg_D_ = 0.0;
		CalDeltaTKernel();
		//输出分离的局部与区域场
		CalPartialModel(regional_field_);
		if (OutObsPoints(obsout_name, regional_field_)) return -1;
	}
	else cerr << BOLDRED << "error ==> " << RESET << "unknown running type: " << run_type << endl;

	return 0;
}