/*
 * cielocontrol.c:
 */

#include "../include/cielo.h"
#include "../modules.h"
#include "./parallel.h"

#undef __FUNCT__ 
#define __FUNCT__ "cielocontrol"

int main(int argc,char* *argv){
	
	int i,n;
	
	/*I/O: */
	FILE* fid=NULL;
	char* inputfilename=NULL;
	char* outputfilename=NULL;
	char* lockname=NULL;
	char* analysis_type="control";

	/*Finite element model: */
	#if defined(_PARALLEL_) && defined(_HAVE_PETSC_)
	FemModel femmodel;
	#endif

	/*control : */
	#if defined(_PARALLEL_) && defined(_HAVE_PETSC_)
	Vec*             u_g=NULL;
	#endif
	double*          search_vector=NULL;
	int              status;
	ParameterInputs* inputs=NULL;
	int              reloop;
	
	#if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_))
		_printf_("%s%s\n",__FUNCT__," error message: parallel executable was compiled without support of parallel libraries!");
		return 1;
	#else
	
		/*Initialize MPI environment: */
		PetscInitialize(&argc,&argv,(char *)0,"");  

		/*Size and rank: */
		MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);  
		MPI_Comm_size(MPI_COMM_WORLD,&num_procs); 

		/*Some checks on size of cluster*/
		if (num_procs<=1){
			_printf_("\nSize of MPI COMM WORLD is 1, needs to be at least 2. Include more nodes\n"); 
			PetscFinalize(); 
			return 0;
		}

		/*Recover dbdir, input file name and output file name: */
		dbdir=argv[1];
		inputfilename=argv[2];
		outputfilename=argv[3];
		lockname=argv[4];

		/*Open handle to data on disk: */
		fid=fopen(inputfilename,"rb");
		if(fid==NULL){
			_printf_("%s%s\n",__FUNCT__," error message: could not open file ",inputfilename," for binary reading"); 
			return 0;
		}

		/*Read and create finite element model: */
		if(!CreateFemModel(&femmodel,fid,analysis_type)){
			_printf_("%s\n",__FUNCT__," error message: could not read finite element model!\n");
			return 0;
		}

		/*Initialize inputs: */
		inputs=NewParameterInputs();

		/*For now, only one parameter is allowed: */
		if (femmodel.workspaceparams->num_control_parameters>1){
			_printf_("%s%s\n",__FUNCT__," error message: multiple control parameters not implemented yet!");
			PetscFinalize();
			return 1;
		}

		/*Initialize search vector: */
		search_vector=xmalloc(femmodel.workspaceparams->num_control_parameters*sizeof(double));

		/*Start looping: */
		for(n=0;n<femmodel.workspaceparams->nsteps;n++){
			
			_printf_("\n%s%i%s%i\n","   control method step ",(n+1),"/",femmodel.workspaceparams->nsteps);

			//initialize inputs, ie parameters on which we invert.
			DeleteParameterInputs(&inputs); inputs=NewParameterInputs();
			
			for(i=0;i<femmodel.workspaceparams->num_control_parameters;i++){
				char* control_type=femmodel.workspaceparams->control_types[i];
				ParameterInputsAddFromMat(inputs,WorkspaceParamsGetParameter(femmodel.workspaceparams,control_type),femmodel.workspaceparams->gsize,control_type);
			}
			ParameterInputsAddFromDouble(inputs,femmodel.workspaceparams->fit[n],"fit");

			//update velocity:
			_printf_("%s\n","      updating velocity...");
			VecFree(&u_g); cielodiagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,&femmodel);
			_printf_("%s\n","      done.");

			//call gradJ module 
			_printf_("%s\n","      computing gradJ...");
			GradJCompute(inputs,&femmodel);
			_printf_("%s\n","      done.");
			
			//normalize directions 
			_printf_("%s\n","      normalizing directions...");
			GradJOrth(femmodel.workspaceparams);
			_printf_("%s\n","      done.");


			//search along the direction for the parameter vector that minimizes the misfit J
			_printf_("%s\n","      optimizing...");
			status=GradJSearch(search_vector,&femmodel,n);
			_printf_("%s\n","      done.");
			
			_printf_("%s%g\n","      misfit: ",femmodel.workspaceparams->J[n]);

			//Check the GradJSearch actually resulting in improving the misfit. Otherwise, reloop.
			reloop=GradJCheck(femmodel.workspaceparams,n,status);
			if (reloop){
				printf("reloop\n",reloop);
				n=n-1;
				continue;
			}

			//update parameters with new optimized values
			_printf_("%s\n","      updating parameters...");
			ParameterUpdate(search_vector,n,femmodel.workspaceparams,femmodel.batchparams);
			_printf_("%s\n","      done.");

			//temporary saving 
			if ((n%1)==0){
				_printf_("%s\n","      saving temporary results...");
				OutputControl(femmodel.workspaceparams,femmodel.batchparams,u_g,femmodel.tpartition,outputfilename,"control");
				_printf_("%s\n","      done.");
			}

		}
		
		/*Write results to disk: */
		_printf_("%s\n","      saving final results...");
		OutputControl(femmodel.workspaceparams,femmodel.batchparams,u_g,femmodel.tpartition,outputfilename,"control");
		_printf_("%s\n","      done.");

		/*Write lock file if requested: */
		if (femmodel.batchparams->waitonlock){
			WriteLockFile(lockname);
		}
			
		cleanup_and_return:
		_printf_("exiting.\n");
		
		/*Synchronize everyone before exiting: */
		MPI_Barrier(MPI_COMM_WORLD);

		/*Close MPI libraries: */
		PetscFinalize(); 
		
		return 0; //unix success return;
	#endif
}
