Changeset 6238

Timestamp:

10/11/10 12:28:01 (14 years ago)

Author:

Mathieu Morlighem

Message:

Fixed Gradient scaling (Added new module)

Location:

issm/trunk/src

Files:

• c/Makefile.am (modified) (2 diffs)
• c/modules/ControlInputScaleGradientx (added)
• c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp (added)
• c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.h (added)
• c/modules/Orthx/Orthx.cpp (modified) (1 diff)
• c/modules/modules.h (modified) (1 diff)
• c/objects/Elements/Element.h (modified) (1 diff)
• c/objects/Elements/Penta.cpp (modified) (1 diff)
• c/objects/Elements/Penta.h (modified) (1 diff)
• c/objects/Elements/Tria.cpp (modified) (1 diff)
• c/objects/Elements/Tria.h (modified) (1 diff)
• c/objects/Inputs/ControlInput.cpp (modified) (1 diff)
• c/objects/Inputs/ControlInput.h (modified) (1 diff)
• c/solutions/gradient_core.cpp (modified) (3 diffs)
• m/solutions/gradient_core.m (modified) (3 diffs)
• mex/ControlInputScaleGradient (added)
• mex/ControlInputScaleGradient/ControlInputScaleGradient.cpp (added)
• mex/ControlInputScaleGradient/ControlInputScaleGradient.h (added)
• mex/Makefile.am (modified) (2 diffs)

issm/trunk/src/c/Makefile.am ¶

@@ -551 +551 @@
                                         ./modules/ControlInputSetGradientx/ControlInputSetGradientx.cpp\
                                         ./modules/ControlInputSetGradientx/ControlInputSetGradientx.h\
+                                        ./modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp\
+                                        ./modules/ControlInputScaleGradientx/ControlInputScaleGradientx.h\
                                         ./modules/SystemMatricesx/SystemMatricesx.cpp\
                                         ./modules/SystemMatricesx/SystemMatricesx.h\
@@ -1115 +1117 @@
                                         ./modules/ControlInputSetGradientx/ControlInputSetGradientx.cpp\
                                         ./modules/ControlInputSetGradientx/ControlInputSetGradientx.h\
+                                        ./modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp\
+                                        ./modules/ControlInputScaleGradientx/ControlInputScaleGradientx.h\
                                         ./modules/SystemMatricesx/SystemMatricesx.cpp\
                                         ./modules/SystemMatricesx/SystemMatricesx.h\

issm/trunk/src/c/modules/Orthx/Orthx.cpp ¶

@@ -25 +25 @@
         }
 
-        /*scale to 1: gradient=gradient/max(abs(gradient))*/
-        VecNorm(newgradj,NORM_INFINITY,&norm_new);
-        if (norm_new<=0){
-                ISSMERROR("||∂J/∂α||∞ = 0  gradient is zero");
-        }
-        if(isnan(norm_new)){
-                ISSMERROR("input vector norm is NaN\n");
-        }
-        VecScale(newgradj,1.0/norm_new);
-
         /*Assign correct pointer*/
         *pnewgradj=newgradj;

issm/trunk/src/c/modules/modules.h ¶

@@ -21 +21 @@
 #include "./ControlInputGetGradientx/ControlInputGetGradientx.h"
 #include "./ControlInputSetGradientx/ControlInputSetGradientx.h"
+#include "./ControlInputScaleGradientx/ControlInputScaleGradientx.h"
 #include "./CostFunctionx/CostFunctionx.h"
 #include "./CreateNodalConstraintsx/CreateNodalConstraintsx.h"

issm/trunk/src/c/objects/Elements/Element.h ¶

@@ -60 +60 @@
                 virtual void   ControlInputGetGradient(Vec gradient,int enum_type)=0;
                 virtual void   ControlInputSetGradient(double* gradient,int enum_type)=0;
+                virtual void   ControlInputScaleGradient(int enum_type, double scale)=0;
                 virtual void   ProcessResultsUnits(void)=0;
                 virtual void   MinVel(double* pminvel, bool process_units)=0;

issm/trunk/src/c/objects/Elements/Penta.cpp ¶

@@ -816 +816 @@
         ((ControlInput*)input)->GetGradient(gradient,&doflist1[0]);
 
+}/*}}}*/
+/*FUNCTION Penta::ControlInputScaleGradient{{{1*/
+void Penta::ControlInputScaleGradient(int enum_type,double scale){
+
+        Input* input=NULL;
+
+        if(enum_type==RheologyBbarEnum){
+                input=(Input*)matice->inputs->GetInput(enum_type);
+        }
+        else{
+                input=inputs->GetInput(enum_type);
+        }
+        if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type));
+        if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type));
+
+        ((ControlInput*)input)->ScaleGradient(scale);
 }/*}}}*/
 /*FUNCTION Penta::ControlInputSetGradient{{{1*/

issm/trunk/src/c/objects/Elements/Penta.h ¶

@@ -95 +95 @@
                 void   InputScale(int enum_type,double scale_factor);
                 void   ControlInputGetGradient(Vec gradient,int enum_type);
+                void   ControlInputScaleGradient(int enum_type,double scale);
                 void   ControlInputSetGradient(double* gradient,int enum_type);
                 void   InputToResult(int enum_type,int step,double time);

issm/trunk/src/c/objects/Elements/Tria.cpp ¶

@@ -793 +793 @@
 
 }/*}}}*/
+/*FUNCTION Tria::ControlInputScaleGradient{{{1*/
+void Tria::ControlInputScaleGradient(int enum_type,double scale){
+
+        Input* input=NULL;
+
+        if(enum_type==RheologyBbarEnum){
+                input=(Input*)matice->inputs->GetInput(enum_type);
+        }
+        else{
+                input=inputs->GetInput(enum_type);
+        }
+        if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type));
+        if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type));
+
+        ((ControlInput*)input)->ScaleGradient(scale);
+}/*}}}*/
 /*FUNCTION Tria::ControlInputSetGradient{{{1*/
 void Tria::ControlInputSetGradient(double* gradient,int enum_type){

issm/trunk/src/c/objects/Elements/Tria.h ¶

@@ -96 +96 @@
                 void   InputScale(int enum_type,double scale_factor);
                 void   ControlInputGetGradient(Vec gradient,int enum_type);
+                void   ControlInputScaleGradient(int enum_type,double scale);
                 void   ControlInputSetGradient(double* gradient,int enum_type);
                 void   InputToResult(int enum_type,int step,double time);

issm/trunk/src/c/objects/Inputs/ControlInput.cpp ¶

@@ -266 +266 @@
         if(gradient) gradient->GetVectorFromInputs(gradient_vec,doflist);
 }/*}}}*/
+/*FUNCTION ControlInput::ScaleGradient{{{1*/
+void ControlInput::ScaleGradient(double scaling_factor){
+        if(!gradient) ISSMERROR("Gradient of ControlInput %s not found",EnumToString(enum_type));
+        gradient->Scale(scaling_factor);
+}/*}}}*/
 /*FUNCTION ControlInput::SetGradient{{{1*/
 void ControlInput::SetGradient(Input* gradient_in){

issm/trunk/src/c/objects/Inputs/ControlInput.h ¶

@@ -78 +78 @@
                 ElementResult* SpawnGradient(int step, double time);
                 void GetGradient(Vec gradient_vec,int* doflist);
+                void ScaleGradient(double scale);
                 void SetGradient(Input* gradient_in);
                 void UpdateValue(double scalar);

issm/trunk/src/c/solutions/gradient_core.cpp ¶

@@ -21 +21 @@
         int    *control_type   = NULL;
         double *optscal_list   = NULL;
+        double  optscal,norm_grad;
 
         /*Intermediaries*/
@@ -34 +35 @@
         femmodel->parameters->FindParam(&verbose,VerboseEnum);
 
+        /*Compute and norm gradient of all controls*/
         for (int i=0;i<num_controls;i++){
 
@@ -44 +46 @@
                         _printf_("%s","      orthogonalization...\n");
                         ControlInputGetGradientx(&old_gradient,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,control_type[i]);
-                        Orthx(&new_gradient,gradient,old_gradient);
-                        VecFree(&old_gradient);
+                        Orthx(&new_gradient,gradient,old_gradient); VecFree(&old_gradient); VecFree(&gradient);
                 }
                 else{ 
                         _printf_("%s","      normalizing directions...\n");
-                        Orthx(&new_gradient,gradient,NULL);
+                        Orthx(&new_gradient,gradient,NULL); VecFree(&gradient);
                 }
-                VecFree(&gradient);
 
-                /*Scale gradient for current step and current parameter*/
-                VecScale(new_gradient,optscal_list[num_controls*step+i]);
+                /*Get scaling factor of current control:*/
+                VecNorm(new_gradient,NORM_INFINITY,&norm_grad);
+                if(norm_grad<=0)    ISSMERROR("||∂J/∂α||∞ = 0    gradient norm of J with respect to %s is zero",EnumToString(control_type[i]));
+                if(isnan(norm_grad))ISSMERROR("||∂J/∂α||∞ = NaN  gradient norm of J with respect to %s is NaN" ,EnumToString(control_type[i]));
+                if(i==0 || (optscal_list[num_controls*step+i]/norm_grad)<optscal) optscal=optscal_list[num_controls*step+i]/norm_grad;
 
                 /*plug back into inputs: */
                 ControlInputSetGradientx(femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmodel->parameters,control_type[i],new_gradient);
-
-                /*Free ressources and return:*/
                 VecFree(&new_gradient);
         }
+        printf("norm grad = %g\n",norm_grad);
+        printf("optscal = %g\n",optscal);
+
+        /*Scale Gradients*/
+        for (int i=0;i<num_controls;i++) ControlInputScaleGradientx(femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmodel->parameters,control_type[i],optscal);
 
         /*Clean up and return*/

issm/trunk/src/m/solutions/gradient_core.m ¶

@@ -25 +25 @@
         control_type=femmodel.parameters.ControlType;
         control_steady=femmodel.parameters.ControlSteady;
-        optscal=femmodel.parameters.OptScal;
+        optscal_list=femmodel.parameters.OptScal;
 
         for i=1:num_controls,
@@ -32 +32 @@
                 grad=Gradj(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 
-
-                if control_steady;
+                if control_steady,
                         femmodel=diagnostic_core(femmodel);
                 end
@@ -46 +45 @@
                 end
 
-                %Scale Gradient
-                new_gradient=optscal(step,i)*new_gradient;
+                 %Get scaling factor of current control:
+                 norm_grad=norm(new_gradient,inf);
+                 if(norm_grad<=0),     error(['||∂J/∂α||∞ = 0   gradient norm of J with respect to ' EnumToString(control_type(i))  ' is zero']); end
+                 if(isnan(norm_grad)), error(['||∂J/∂α||∞ = NaN gradient norm of J with respect to ' EnumToString(control_type(i))  ' is NaN' ]); end
+                 if(i==1 | (optscal_list(step,i)/norm_grad)<optscal) optscal=optscal_list(step,i)/norm_grad; end
 
                 %plug back into inputs:
                 [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=ControlInputSetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,  femmodel.parameters,control_type(i),new_gradient);
         end
+
+        %Scale all gradients
+        for i=1:num_controls,
+                [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=ControlInputScaleGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,  femmodel.parameters,control_type(i),optscal);
+        end

issm/trunk/src/mex/Makefile.am ¶

@@ -15 +15 @@
                                 ControlOptimization\
                                 ControlInputGetGradient\
+                                ControlInputScaleGradient\
                                 ControlInputSetGradient\
                                 ContourToMesh \
@@ -133 +134 @@
                                                                                 ControlInputGetGradient/ControlInputGetGradient.h
 
+ControlInputScaleGradient_SOURCES = ControlInputScaleGradient/ControlInputScaleGradient.cpp\
+                                                                                         ControlInputScaleGradient/ControlInputScaleGradient.h
+
 ControlInputSetGradient_SOURCES = ControlInputSetGradient/ControlInputSetGradient.cpp\
                           ControlInputSetGradient/ControlInputSetGradient.h