/*
 * objectivefunctionC.c:
 */
#include "../../../config.h"

#if defined(_PARALLEL_) && defined(_HAVE_PETSC_)

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

#undef __FUNCT__ 
#define __FUNCT__ "objectivefunctionC"

double objectivefunctionC(double* search_vector,double fit,double optscal,FemModel*  femmodel,ParameterInputs* inputs){
	
		
	/*Error management: */
	int noerr=1;
	int i;
	int j;
	double*          parameter=NULL;
	double*          newparameter=NULL;
	Vec*             vec_gradient=NULL;
	double*          gradient=NULL;
	Vec*             u_g=NULL;
	double*          u_g_double=NULL;

	/*output: */
	double           J;

	/*inputs recovered from the FemModel: */
	WorkspaceParams* workspaceparams=NULL;
	BatchParams*        batchparams=NULL;
	DataSet*            bgpdt=NULL;
	DataSet*            bgpdtb=NULL;
	DataSet*            est=NULL;
	DataSet*            lst=NULL;
	DataSet*            ept=NULL;
	DataSet*            mpt=NULL;
	DataSet*            geom3=NULL;
	Mat*                G_mn=NULL;
	Vec*                y_s=NULL;
	int                 flag_y_s=-1;
	int                 analysis_type=-1;

	/*Recover model inputs: */
	workspaceparams=femmodel->workspaceparams;
	batchparams=femmodel->batchparams;
	bgpdt=femmodel->bgpdt;
	bgpdtb=femmodel->bgpdtb;
	est=femmodel->est;
	lst=femmodel->lst;
	ept=femmodel->ept;
	mpt=femmodel->mpt;
	geom3=femmodel->geom3;
	G_mn=femmodel->G_mn;
	y_s=femmodel->y_s;
	flag_y_s=femmodel->flag_y_s;
	analysis_type=femmodel->batchparams->analysis_type_enum;

	//Go through parameters, and update along gradients, multiplying by search_vector.
	for(i=0;i<workspaceparams->num_control_parameters;i++){

		
		char* control_type=workspaceparams->control_types[i];

		/*Get parameter, and gradient for the parameter: */
		parameter=WorkspaceParamsGetParameter(workspaceparams,control_type);
		vec_gradient=WorkspaceParamsGetParameterGradient(workspaceparams,control_type);
	
	
		/*serialize gradient: */
		VecToMPISerial(&gradient,vec_gradient);


		/*Ok, for this parameter, we have a direction. Update the parameter along 
		 * the direction, using the search_vector value. Now, because parameter comes
		 * directly from a workspaceparams pointer, we do not want to modify it. Make a copy
		 * before updating it.*/
		newparameter=xmalloc(workspaceparams->gsize*sizeof(double));
		memcpy(newparameter,parameter,workspaceparams->gsize*sizeof(double));


		for(j=0;j<(int)(workspaceparams->gsize/6);j++){
			newparameter[6*j+0]=parameter[6*j+0]-search_vector[i]*optscal*gradient[6*j+0];
		}
		

		/*Add newparameter to inputs: */
		ParameterInputsAddFromMat(inputs,newparameter,workspaceparams->gsize,control_type);

		/*Now, call on  each parameter's self constraining routine: */
		WorkspaceParamsConstrain(workspaceparams,control_type);
	}

	//Run diagnostic with updated parameters.
	cielodiagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,femmodel);
	VecToMPISerial(&u_g_double,u_g);

	//Compute misfit for this velocity field. 
	ParameterInputsAddFromDouble(inputs,fit,"fit");
	Misfitx( &J, femmodel->bgpdt, femmodel->bgpdtb, femmodel->est, femmodel->ept, femmodel->mpt,femmodel->geom3, u_g_double, workspaceparams->u_g_obs,inputs, femmodel->batchparams->analysis_type_enum);

	/*Free ressources: */
	xfree((void**)&gradient);
	xfree((void**)&newparameter);
	xfree((void**)&u_g_double);

	return J;
}
#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_)

