/*!\file:  gradient_core.cpp
 * \brief compute inverse method gradient direction.
 */ 

#include "../toolkits/toolkits.h"
#include "../classes/objects/objects.h"
#include "../shared/shared.h"
#include "../io/io.h"
#include "../EnumDefinitions/EnumDefinitions.h"
#include "./solutions.h"
#include "../modules/modules.h"
#include "../include/include.h"
#include "../solvers/solvers.h"

void gradient_core(FemModel* femmodel,int step,bool orthogonalize){ 

	/*Intermediaries*/
	IssmDouble  norm_inf;
	IssmDouble *norm_list    = NULL;
	Vector*     new_gradient = NULL;
	Vector*     gradient     = NULL;
	Vector*     old_gradient = NULL;

	/*Compute gradient*/
	if(VerboseControl()) _pprintLine_("   compute cost function gradient");
	Gradjx(&gradient,&norm_list,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters);

	if (orthogonalize){
		if(VerboseControl()) _pprintLine_("   orthogonalization");
		ControlInputGetGradientx(&old_gradient,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters);
		Orthx(&new_gradient,gradient,old_gradient); xdelete(&old_gradient); xdelete(&gradient);
	}
	else{ 
		new_gradient=gradient;
	}

	/*Check that gradient is clean*/
	norm_inf=new_gradient->Norm(NORM_INF);
	if(norm_inf<=0)    _error2_("||∂J/∂α||∞ = 0    gradient norm is zero");
	if(xIsNan<IssmDouble>(norm_inf))_error2_("||∂J/∂α||∞ = NaN  gradient norm is NaN");

	/*plug back into inputs: */
	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,new_gradient);
	xdelete(&new_gradient);

	/*Scale Gradients*/
	ControlInputScaleGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,norm_list,step);

	/*Clean up and return*/
	xDelete<IssmDouble>(norm_list);
}
