source: issm/oecreview/Archive/17984-18295/ISSM-18127-18128.diff@ 18296

Last change on this file since 18296 was 18296, checked in by Mathieu Morlighem, 11 years ago

Added 17984-18295
File size: 54.6 KB

  • ../trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp

     
    1919#include "./types.h"
    2020#include "./isnan.h"
    2121
    22 void BrentSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,void*),void* optargs){
     22void BrentSearch(IssmDouble** pJ,OptPars optpars,IssmDouble* X0,IssmDouble (*f)(IssmDouble*,void*),IssmDouble (*g)(IssmDouble**,IssmDouble*,void*),void* usr){
    2323
    2424        /* This routine is optimizing a given function using Brent's method
    2525         * (Golden or parabolic procedure)*/
    2626
    2727        /*Intermediary*/
     28        int        iter;
    2829        IssmDouble si,gold,intervalgold,oldintervalgold;
    29         IssmDouble parab_num,parab_den;
    30         IssmDouble distance,cm_jump;
     30        IssmDouble parab_num,parab_den,distance;
    3131        IssmDouble fxmax,fxmin,fxbest;
    3232        IssmDouble fx,fx1,fx2;
    33         IssmDouble xmax,xmin,xbest;
    34         IssmDouble x,x1,x2,xm;
     33        IssmDouble x,x1,x2,xm,xbest;
    3534        IssmDouble tol1,tol2,seps;
    3635        IssmDouble tolerance = 1.e-4;
    37         int        maxiter ,iter;
    38         bool       loop= true,goldenflag;
    3936
    4037        /*Recover parameters:*/
    41         xmin=optpars->xmin;
    42         xmax=optpars->xmax;
    43         maxiter=optpars->maxiter;
    44         cm_jump=optpars->cm_jump;
     38        int         nsteps  = optpars.nsteps;
     39        int         nsize   = optpars.nsize;
     40        IssmDouble  xmin    = optpars.xmin;
     41        IssmDouble  xmax    = optpars.xmax;
     42        int        *maxiter = optpars.maxiter;
     43        IssmDouble *cm_jump = optpars.cm_jump;
    4544
    46         /*initialize counter and get response at the boundaries*/
    47         iter=0;
    48         fxmin = (*f)(xmin,optargs);
    49         if (xIsNan<IssmDouble>(fxmin)) _error_("Function evaluation returned NaN");
    50         cout<<setprecision(5);
    51         if(VerboseControl()) _printf0_("\n");
    52         if(VerboseControl()) _printf0_("       Iteration         x           f(x)       Tolerance         Procedure\n");
    53         if(VerboseControl()) _printf0_("\n");
    54         if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmin<<"  "<<setw(12)<<fxmin<<"           N/A         boundary\n");
    55         fxmax = (*f)(xmax,optargs);
    56         if (xIsNan<IssmDouble>(fxmax)) _error_("Function evaluation returned NaN");
    57         if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmax<<"  "<<setw(12)<<fxmax<<"           N/A         boundary\n");
     45        /*Initialize gradient and controls*/
     46        IssmDouble* G = NULL;
     47        IssmDouble* J = xNew<IssmDouble>(nsteps);
     48        IssmDouble* X = xNew<IssmDouble>(nsize);
    5849
    59         /*test if jump option activated and xmin==0*/
    60         if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && (fxmax/fxmin)<cm_jump){
    61                 *psearch_scalar=xmax;
    62                 *pJ=fxmax;
    63                 return;
    64         }
     50        /*start iterations*/
     51        for(int n=0;n<nsteps;n++){
    6552
    66         /*initialize optimization variables*/
    67         seps=sqrt(DBL_EPSILON);    //precision of a IssmDouble
    68         distance=0.0;              //new_x=old_x + distance
    69         gold=0.5*(3.0-sqrt(5.0));  //gold = 1 - golden ratio
    70         intervalgold=0.0;          //distance used by Golden procedure
     53                /*Reset some variables*/
     54                iter = 0;
     55                xmin = 0.;
     56                xmax = 1.;
     57                bool loop = true;
     58                cout<<setprecision(5);
    7159
    72         /*1: initialize the values of the 4 x needed (x1,x2,x,xbest)*/
    73         x1=xmin+gold*(xmax-xmin);
    74         x2=x1;
    75         xbest=x1;
    76         x=xbest;
     60                /*Get current Gradient at xmin=0*/
     61                if(VerboseControl()) _printf0_("\n" << "   step " << n+1 << "/" << nsteps << "\n");
     62                fxmin = (*g)(&G,X0,usr); if(xIsNan<IssmDouble>(fxmin)) _error_("Function evaluation returned NaN");
     63                if(VerboseControl()) _printf0_("\n");
     64                if(VerboseControl()) _printf0_("       Iteration         x           f(x)       Tolerance         Procedure\n");
     65                if(VerboseControl()) _printf0_("\n");
     66                if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmin<<"  "<<setw(12)<<fxmin<<"           N/A         boundary\n");
     67               
     68                /*Get f(xmax)*/
     69                for(int i=0;i<nsize;i++) X[i]=X0[i]+xmax*G[i];
     70                fxmax = (*f)(X,usr); if (xIsNan<IssmDouble>(fxmax)) _error_("Function evaluation returned NaN");
     71                if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmax<<"  "<<setw(12)<<fxmax<<"           N/A         boundary\n");
    7772
    78         /*2: call the function to be evaluated*/
    79         fxbest = (*f)(x,optargs);
    80         if(xIsNan<IssmDouble>(fxbest)) _error_("Function evaluation returned NaN");
    81         iter=iter+1;
     73                /*test if jump option activated and xmin==0*/
     74                if(!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && (fxmax/fxmin)<cm_jump[n]){
     75                        for(int i=0;i<nsize;i++) X0[i]=X0[i]+xmax*G[i];
     76                        xDelete<IssmDouble>(G);
     77                        J[n]=fxmax;
     78                        continue;
     79                }
    8280
    83         /*3: update the other variables*/
    84         fx1=fxbest;
    85         fx2=fxbest;
    86         /*xm is always in the middle of a and b*/
    87         xm=0.5*(xmin+xmax);                           
    88         /*update tolerances*/
    89         tol1=seps*sqrt(pow(xbest,2))+tolerance/3.0;
    90         tol2=2.0*tol1;
     81                /*initialize optimization variables*/
     82                seps=sqrt(DBL_EPSILON);    //precision of a IssmDouble
     83                distance=0.0;              //new_x=old_x + distance
     84                gold=0.5*(3.0-sqrt(5.0));  //gold = 1 - golden ratio
     85                intervalgold=0.0;          //distance used by Golden procedure
    9186
    92         /*4: print result*/
    93         if(VerboseControl())
    94          _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<xbest<<"  "<<setw(12)<<fxbest<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<"         initial\n");
    95         if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && ((fxbest/fxmin)<cm_jump)){
    96                 if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump << "\n");
    97                 loop=false;
    98         }
     87                /*1: initialize the values of the 4 x needed (x1,x2,x,xbest)*/
     88                x1=xmin+gold*(xmax-xmin);
     89                x2=x1;
     90                xbest=x1;
     91                x=xbest;
    9992
    100         while(loop){
     93                /*2: call the function to be evaluated*/
     94                for(int i=0;i<nsize;i++) X[i]=X0[i]+x*G[i];
     95                fxbest = (*f)(X,usr); if(xIsNan<IssmDouble>(fxbest)) _error_("Function evaluation returned NaN");
     96                iter++;
    10197
    102                 goldenflag=true;
     98                /*3: update the other variables*/
     99                fx1=fxbest;
     100                fx2=fxbest;
     101                /*xm is always in the middle of a and b*/
     102                xm=0.5*(xmin+xmax);                           
     103                /*update tolerances*/
     104                tol1=seps*sqrt(pow(xbest,2))+tolerance/3.0;
     105                tol2=2.0*tol1;
    103106
    104                 // Is a parabolic fit possible ?
    105                 if (sqrt(pow(intervalgold,2))>tol1){
     107                /*4: print result*/
     108                if(VerboseControl())
     109                 _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<xbest<<"  "<<setw(12)<<fxbest<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<"         initial\n");
     110                if (!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && ((fxbest/fxmin)<cm_jump[n])){
     111                        if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump[n] << "\n");
     112                        loop=false;
     113                }
    106114
    107                         // Yes, so fit parabola
    108                         goldenflag=false;
    109                         parab_num=(xbest-x1)*(xbest-x1)*(fxbest-fx2)-(xbest-x2)*(xbest-x2)*(fxbest-fx1);;
    110                         parab_den=2.0*(xbest-x1)*(fxbest-fx2)-2.0*(xbest-x2)*(fxbest-fx1);
     115                while(loop){
    111116
    112                         //reverse p if necessary
    113                         if(parab_den>0.0){
    114                                 parab_num=-parab_num;
    115                         }
    116                         parab_den=sqrt(pow(parab_den,2));
    117                         oldintervalgold=intervalgold;
    118                         intervalgold=distance;
     117                        bool goldenflag=true;
    119118
    120                         // Is the parabola acceptable (we use seps here because in some configuration parab_num==parab_den*(xmax-xbest)
    121                         // and the result is not repeatable anymore
    122                         if (( sqrt(pow(parab_num,2)) < sqrt(pow(0.5*parab_den*oldintervalgold,2))) &&
    123                                                 (parab_num>parab_den*(xmin-xbest)+seps) &&
    124                                                 (parab_num<parab_den*(xmax-xbest)-seps)){
     119                        // Is a parabolic fit possible ?
     120                        if (sqrt(pow(intervalgold,2))>tol1){
    125121
    126                                 // Yes, parabolic interpolation step
    127                                 distance=parab_num/parab_den;
    128                                 x=xbest+distance;
     122                                // Yes, so fit parabola
     123                                goldenflag=false;
     124                                parab_num=(xbest-x1)*(xbest-x1)*(fxbest-fx2)-(xbest-x2)*(xbest-x2)*(fxbest-fx1);;
     125                                parab_den=2.0*(xbest-x1)*(fxbest-fx2)-2.0*(xbest-x2)*(fxbest-fx1);
    129126
    130                                 // f must not be evaluated too close to min_x or max_x
    131                                 if (((x-xmin)<tol2) || ((xmax-x)<tol2)){
    132                                         if ((xm-xbest)<0.0) si=-1;
    133                                         else                si=1;
    134                                         //compute new distance
    135                                         distance=tol1*si;
     127                                //reverse p if necessary
     128                                if(parab_den>0.0){
     129                                        parab_num=-parab_num;
    136130                                }
     131                                parab_den=sqrt(pow(parab_den,2));
     132                                oldintervalgold=intervalgold;
     133                                intervalgold=distance;
     134
     135                                // Is the parabola acceptable (we use seps here because in some configuration parab_num==parab_den*(xmax-xbest)
     136                                // and the result is not repeatable anymore
     137                                if (( sqrt(pow(parab_num,2)) < sqrt(pow(0.5*parab_den*oldintervalgold,2))) &&
     138                                                        (parab_num>parab_den*(xmin-xbest)+seps) &&
     139                                                        (parab_num<parab_den*(xmax-xbest)-seps)){
     140
     141                                        // Yes, parabolic interpolation step
     142                                        distance=parab_num/parab_den;
     143                                        x=xbest+distance;
     144
     145                                        // f must not be evaluated too close to min_x or max_x
     146                                        if (((x-xmin)<tol2) || ((xmax-x)<tol2)){
     147                                                if ((xm-xbest)<0.0) si=-1;
     148                                                else                si=1;
     149                                                //compute new distance
     150                                                distance=tol1*si;
     151                                        }
     152                                }
     153                                else{
     154                                        // Not acceptable, must do a golden section step
     155                                        goldenflag=true;
     156                                }
    137157                        }
    138                         else{
    139                                 // Not acceptable, must do a golden section step
    140                                 goldenflag=true;
    141                         }
    142                 }
    143158
    144                 //Golden procedure
    145                 if(goldenflag){
    146                         // compute the new distance d
    147                         if(xbest>=xm){
    148                                 intervalgold=xmin-xbest;   
     159                        //Golden procedure
     160                        if(goldenflag){
     161                                // compute the new distance d
     162                                if(xbest>=xm){
     163                                        intervalgold=xmin-xbest;   
     164                                }
     165                                else{
     166                                        intervalgold=xmax-xbest; 
     167                                }
     168                                distance=gold*intervalgold;
    149169                        }
    150                         else{
    151                                 intervalgold=xmax-xbest; 
    152                         }
    153                         distance=gold*intervalgold;
    154                 }
    155170
    156                 // The function must not be evaluated too close to xbest
    157                 if(distance<0) si=-1;
    158                 else           si=1;
    159                 if(sqrt(pow(distance,2))>tol1) x=xbest+si*sqrt(pow(distance,2));
    160                 else                           x=xbest+si*tol1;
     171                        // The function must not be evaluated too close to xbest
     172                        if(distance<0) si=-1;
     173                        else           si=1;
     174                        if(sqrt(pow(distance,2))>tol1) x=xbest+si*sqrt(pow(distance,2));
     175                        else                           x=xbest+si*tol1;
    161176
    162                 //evaluate function on x
    163                 fx = (*f)(x,optargs);
    164                 if(xIsNan<IssmDouble>(fx)) _error_("Function evaluation returned NaN");
    165                 iter=iter+1;
     177                        //evaluate function on x
     178                        for(int i=0;i<nsize;i++) X[i]=X0[i]+x*G[i];
     179                        fx = (*f)(X,usr); if(xIsNan<IssmDouble>(fx)) _error_("Function evaluation returned NaN");
     180                        iter++;
    166181
    167                 // Update a, b, xm, x1, x2, tol1, tol2
    168                 if (fx<=fxbest){
    169                         if (x>=xbest) xmin=xbest;
    170                         else          xmax=xbest;
    171                         x1=x2;    fx1=fx2;
    172                         x2=xbest; fx2=fxbest;
    173                         xbest=x;  fxbest=fx;
    174                 }
    175                 else{ // fx > fxbest
    176                         if (x<xbest) xmin=x;
    177                         else         xmax=x;
    178                         if ((fx<=fx2) || (x2==xbest)){
    179                                 x1=x2; fx1=fx2;
    180                                 x2=x;  fx2=fx;
     182                        // Update a, b, xm, x1, x2, tol1, tol2
     183                        if (fx<=fxbest){
     184                                if (x>=xbest) xmin=xbest;
     185                                else          xmax=xbest;
     186                                x1=x2;    fx1=fx2;
     187                                x2=xbest; fx2=fxbest;
     188                                xbest=x;  fxbest=fx;
    181189                        }
    182                         else if ( (fx <= fx1) || (x1 == xbest) || (x1 == x2) ){
    183                                 x1=x;  fx1=fx;
     190                        else{ // fx > fxbest
     191                                if (x<xbest) xmin=x;
     192                                else         xmax=x;
     193                                if ((fx<=fx2) || (x2==xbest)){
     194                                        x1=x2; fx1=fx2;
     195                                        x2=x;  fx2=fx;
     196                                }
     197                                else if ( (fx <= fx1) || (x1 == xbest) || (x1 == x2) ){
     198                                        x1=x;  fx1=fx;
     199                                }
    184200                        }
    185                 }
    186                 xm = 0.5*(xmin+xmax);
    187                 tol1=seps*pow(pow(xbest,2),0.5)+tolerance/3.0;
    188                 tol2=2.0*tol1;
    189                 if(VerboseControl())
    190                  _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<x<<"  "<<setw(12)<<fx<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<
    191                                          "         "<<(goldenflag?"golden":"parabolic")<<"\n");
     201                        xm = 0.5*(xmin+xmax);
     202                        tol1=seps*pow(pow(xbest,2),0.5)+tolerance/3.0;
     203                        tol2=2.0*tol1;
     204                        if(VerboseControl())
     205                         _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<x<<"  "<<setw(12)<<fx<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<
     206                                                 "         "<<(goldenflag?"golden":"parabolic")<<"\n");
    192207
    193                 /*Stop the optimization?*/
    194                 if (sqrt(pow(xbest-xm,2)) < (tol2-0.5*(xmax-xmin))){
    195                         if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'tolx'=" << tolerance << "\n");
    196                         loop=false;
     208                        /*Stop the optimization?*/
     209                        if (sqrt(pow(xbest-xm,2)) < (tol2-0.5*(xmax-xmin))){
     210                                if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'tolx'=" << tolerance << "\n");
     211                                loop=false;
     212                        }
     213                        else if (iter>=maxiter[n]){
     214                                if(VerboseControl()) _printf0_("      exiting: Maximum number of iterations has been exceeded  ('maxiter'=" << maxiter[n] << ")\n");
     215                                loop=false;
     216                        }
     217                        else if (!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && ((fxbest/fxmin)<cm_jump[n])){
     218                                if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump[n] << "\n");
     219                                loop=false;
     220                        }
     221                        else{
     222                                //continue
     223                                loop=true;
     224                        }
     225                }//end while
     226
     227                //Now, check that the value on the boundaries are not better than current fxbest
     228                if (fxbest>fxmin){
     229                        xbest=optpars.xmin; fxbest=fxmin;
    197230                }
    198                 else if (iter>=maxiter){
    199                         if(VerboseControl()) _printf0_("      exiting: Maximum number of iterations has been exceeded  ('maxiter'=" << maxiter << ")\n");
    200                         loop=false;
     231                if (fxbest>fxmax){
     232                        xbest=optpars.xmax; fxbest=fxmax;
    201233                }
    202                 else if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && ((fxbest/fxmin)<cm_jump)){
    203                         if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump << "\n");
    204                         loop=false;
    205                 }
    206                 else{
    207                         //continue
    208                         loop=true;
    209                 }
    210         }//end while
    211234
    212         //Now, check that the value on the boundaries are not better than current fxbest
    213         if (fxbest>fxmin){
    214                 xbest=optpars->xmin; fxbest=fxmin;
     235                /*Assign output pointers: */
     236                for(int i=0;i<nsize;i++) X0[i]=X0[i]+xbest*G[i];
     237                xDelete<IssmDouble>(G);
     238                J[n]=fxbest;
    215239        }
    216         if (fxbest>fxmax){
    217                 xbest=optpars->xmax; fxbest=fxmax;
    218         }
    219 
    220         /*Assign output pointers: */
    221         *psearch_scalar=xbest;
    222         *pJ=fxbest;
     240       
     241        /*return*/
     242        xDelete<IssmDouble>(X);
     243        *pJ=J;
    223244}

  • ../trunk-jpl/src/c/shared/Numerics/numerics.h

     
    2929
    3030int         min(int a,int b);
    3131int         max(int a,int b);
    32 void        BrentSearch(IssmDouble *psearch_scalar,IssmDouble*pJ,OptPars*optpars,IssmDouble (*f)(IssmDouble,void*),void* optargs);
     32void        BrentSearch(IssmDouble** pJ,OptPars optpars,IssmDouble* X0,IssmDouble (*f)(IssmDouble*,void*),IssmDouble (*g)(IssmDouble**,IssmDouble*,void*),void* usr);
    3333void        cross(IssmDouble *result,IssmDouble*vector1,IssmDouble*vector2);
    3434void        XZvectorsToCoordinateSystem(IssmDouble *T,IssmDouble*xzvectors);
    3535int         cubic(IssmDouble a, IssmDouble b, IssmDouble c, IssmDouble d, double X[3], int *num);

  • ../trunk-jpl/src/c/shared/Numerics/OptPars.h

     
    99
    1010struct OptPars{
    1111
    12         IssmDouble xmin;
    13         IssmDouble xmax;
    14         IssmDouble cm_jump;
    15         int maxiter;
     12        IssmDouble  xmin;
     13        IssmDouble  xmax;
     14        IssmDouble *cm_jump;
     15        int* maxiter;
     16        int  nsteps;
     17        int  nsize;
    1618
    1719};
    1820

  • ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp

     
    1515        int         num_control_type;
    1616        int         num_cm_responses;
    1717        int        *control_type     = NULL;
     18        int        *maxiter          = NULL;
    1819        int        *cm_responses     = NULL;
    1920        IssmDouble *cm_jump          = NULL;
    2021        IssmDouble *optscal          = NULL;
    21         IssmDouble *maxiter          = NULL;
    2222
    2323        /*retrieve some parameters: */
    2424        iomodel->Constant(&control_analysis,InversionIscontrolEnum);
    … …  
    5555                                iomodel->FetchData(&maxiter,NULL,NULL,InversionMaxiterPerStepEnum);
    5656                                parameters->AddObject(new DoubleMatParam(InversionGradientScalingEnum,optscal,nsteps,num_control_type));
    5757                                parameters->AddObject(new DoubleVecParam(InversionStepThresholdEnum,cm_jump,nsteps));
    58                                 parameters->AddObject(new DoubleVecParam(InversionMaxiterPerStepEnum,maxiter,nsteps));
     58                                parameters->AddObject(new IntVecParam(InversionMaxiterPerStepEnum,maxiter,nsteps));
    5959                                break;
    6060                        case 1:/*TAO*/
    6161                                parameters->AddObject(iomodel->CopyConstantObject(InversionFatolEnum));
    … …  
    8282
    8383                xDelete<int>(control_type);
    8484                xDelete<int>(cm_responses);
     85                xDelete<int>(maxiter);
    8586                iomodel->DeleteData(cm_jump,InversionStepThresholdEnum);
    8687                iomodel->DeleteData(optscal,InversionGradientScalingEnum);
    87                 iomodel->DeleteData(maxiter,InversionMaxiterPerStepEnum);
    8888        }
    8989}

  • ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp

     
    66#include "../../shared/shared.h"
    77#include "../../toolkits/toolkits.h"
    88
    9 void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,IssmDouble* vector){
     9void SetControlInputsFromVectorx(FemModel* femmodel,IssmDouble* vector){
    1010
    1111        int  num_controls;
    1212        int *control_type = NULL;
    1313
    1414        /*Retrieve some parameters*/
    15         parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    16         parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     15        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     16        femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
    1717
    1818        for(int i=0;i<num_controls;i++){
    19                 for(int j=0;j<elements->Size();j++){
    20                         Element* element=(Element*)elements->GetObjectByOffset(j);
     19                for(int j=0;j<femmodel->elements->Size();j++){
     20                        Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
    2121                        element->SetControlInputsFromVector(vector,control_type[i],i);
    2222                }
    2323        }
    … …  
    2525        xDelete<int>(control_type);
    2626}
    2727
    28 void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,Vector<IssmDouble>* vector){
     28void SetControlInputsFromVectorx(FemModel* femmodel,Vector<IssmDouble>* vector){
    2929
    30         IssmDouble* serial_vector=NULL;
    31 
    32         serial_vector=vector->ToMPISerial();
    33 
    34         SetControlInputsFromVectorx(elements,nodes, vertices, loads, materials, parameters,serial_vector);
    35 
    36         /*Free ressources:*/
     30        IssmDouble* serial_vector=vector->ToMPISerial();
     31        SetControlInputsFromVectorx(femmodel,serial_vector);
    3732        xDelete<IssmDouble>(serial_vector);
    3833}

  • ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h

     
    77#include "../../classes/classes.h"
    88
    99/* local prototypes: */
    10 void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,Vector<IssmDouble>* vector);
    11 void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,IssmDouble* vector);
     10void SetControlInputsFromVectorx(FemModel* femmodel,Vector<IssmDouble>* vector);
     11void SetControlInputsFromVectorx(FemModel* femmodel,IssmDouble* vector);
    1212
    1313#endif

  • ../trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp

     
    4949                gradient->AXPY(gradient_list[i],1.0);
    5050                delete gradient_list[i];
    5151        }
     52        //gradient->Echo();
     53        //_error_("S");
    5254
    5355        /*Check that gradient is clean*/
    5456        norm_inf=gradient->Norm(NORM_INF);

  • ../trunk-jpl/src/c/cores/controlvalidation_core.cpp

     
    7272                for(int i=0;i<n;i++) X[i] = X0[i] + alpha;
    7373
    7474                /*Calculate j(k+alpha delta k) */
    75                 SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
     75                SetControlInputsFromVectorx(femmodel,X);
    7676                solutioncore(femmodel);
    7777                femmodel->CostFunctionx(&j,NULL,NULL);
    7878

  • ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp

     
    105105        }
    106106
    107107        /*Get solution*/
    108         SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
     108        SetControlInputsFromVectorx(femmodel,X);
    109109        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
    110110        femmodel->OutputControlsx(&femmodel->results);
    111111        femmodel->results->AddObject(new GenericExternalResult<double>(JEnum,f,1,0));
    … …  
    130130        int         solution_type;
    131131        FemModel   *femmodel  = (FemModel*)dzs;
    132132
    133         /*Recover responses*/
    134         int         num_responses;
    135         int        *responses = NULL;
    136         femmodel->parameters->FindParam(&responses,&num_responses,InversionCostFunctionsEnum);
     133        /*Recover number of cost functions responses*/
     134        int num_responses;
     135        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    137136
    138137        /*Constrain input vector*/
    139138        IssmDouble  *XL = NULL;
    … …  
    146145        }
    147146
    148147        /*Update control input*/
    149         SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
     148        SetControlInputsFromVectorx(femmodel,X);
    150149
    151150        /*Recover some parameters*/
    152151        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    … …  
    169168
    170169        if(indic==0){
    171170                /*dry run, no gradient required*/
    172                 xDelete<int>(responses);
    173171                xDelete<IssmDouble>(XU);
    174172                xDelete<IssmDouble>(XL);
    175173                return;
    … …  
    192190        }
    193191
    194192        /*Clean-up and return*/
    195         xDelete<int>(responses);
    196193        xDelete<IssmDouble>(XU);
    197194        xDelete<IssmDouble>(XL);
    198195}

  • ../trunk-jpl/src/c/cores/control_core.cpp

     
    1010#include "../solutionsequences/solutionsequences.h"
    1111
    1212/*Local prototypes*/
    13 bool controlconvergence(IssmDouble J, IssmDouble tol_cm);
    14 IssmDouble objectivefunction(IssmDouble search_scalar,void* optargs);
     13IssmDouble FormFunction(IssmDouble* X,void* usr);
     14IssmDouble FormFunctionGradient(IssmDouble** pG,IssmDouble* X,void* usr);
     15typedef struct {
     16        FemModel* femmodel;
     17        int       nsize;
     18} AppCtx;
    1519
    1620void control_core(FemModel* femmodel){
    1721
    1822        int     i;
    1923
    2024        /*parameters: */
    21         int        num_controls;
     25        int        num_controls,nsize;
    2226        int        nsteps;
    2327        IssmDouble tol_cm;
    2428        int        solution_type;
    … …  
    2630        bool       dakota_analysis = false;
    2731
    2832        int        *control_type   = NULL;
    29         IssmDouble *maxiter        = NULL;
     33        int*        maxiter        = NULL;
    3034        IssmDouble *cm_jump        = NULL;
    3135
    3236        /*intermediary: */
    33         IssmDouble search_scalar = 1.;
    3437        OptPars    optpars;
    3538
    3639        /*Solution and Adjoint core pointer*/
    … …  
    6063        if(VerboseControl()) _printf0_("   preparing initial solution\n");
    6164        if(isFS) solutioncore(femmodel);
    6265
    63         /*Initialize cost function: */
    64         J=xNew<IssmDouble>(nsteps);
     66        /*Get initial guess*/
     67        Vector<IssmDouble> *Xpetsc = NULL;
     68        GetVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
     69        IssmDouble* X0 = Xpetsc->ToMPISerial();
     70        Xpetsc->GetSize(&nsize);
     71        delete Xpetsc;
    6572
    6673        /*Initialize some of the BrentSearch arguments: */
    67         optpars.xmin=0; optpars.xmax=1;
     74        optpars.xmin    = 0;
     75        optpars.xmax    = 1;
     76        optpars.nsteps  = nsteps;
     77        optpars.nsize   = nsize;
     78        optpars.maxiter = maxiter;
     79        optpars.cm_jump = cm_jump;
    6880
    69         /*Start looping: */
    70         for(int n=0;n<nsteps;n++){
     81        /*Initialize function argument*/
     82        AppCtx usr;
     83        usr.femmodel = femmodel;
     84        usr.nsize    = nsize;
    7185
    72                 /*Display info*/
    73                 if(VerboseControl()) _printf0_("\n" << "   control method step " << n+1 << "/" << nsteps << "\n");
     86        /*Call Brent optimization*/
     87        BrentSearch(&J,optpars,X0,&FormFunction,&FormFunctionGradient,(void*)&usr);
    7488
    75 
    76                 /*In steady state inversion, compute new temperature field now*/
    77                 if(solution_type==SteadystateSolutionEnum) solutioncore(femmodel);
    78 
    79                 if(VerboseControl()) _printf0_("   compute adjoint state:\n");
    80                 adjointcore(femmodel);
    81                 gradient_core(femmodel,n,search_scalar==0.);
    82 
    83                 if(VerboseControl()) _printf0_("   optimizing along gradient direction\n");
    84                 optpars.maxiter=reCast<int,IssmDouble>(maxiter[n]); optpars.cm_jump=cm_jump[n];
    85                 BrentSearch(&search_scalar,J+n,&optpars,&objectivefunction,(void*)femmodel);
    86 
    87                 if(VerboseControl()) _printf0_("   updating parameter using optimized search scalar\n"); //true means update save controls
    88                 InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,true);
    89 
    90                 if(controlconvergence(J[n],tol_cm)) break;
     89        if(VerboseControl()) _printf0_("   preparing final solution\n");
     90        IssmDouble  *XL = NULL;
     91        IssmDouble  *XU = NULL;
     92        GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
     93        GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
     94        for(long i=0;i<nsize;i++){
     95                if(X0[i]>XU[i]) X0[i]=XU[i];
     96                if(X0[i]<XL[i]) X0[i]=XL[i];
    9197        }
    92 
    93         if(VerboseControl()) _printf0_("   preparing final solution\n");
     98        xDelete<IssmDouble>(XU);
     99        xDelete<IssmDouble>(XL);
     100        SetControlInputsFromVectorx(femmodel,X0);
    94101        femmodel->parameters->SetParam(true,SaveResultsEnum);
    95102        solutioncore(femmodel);
    96103
    … …  
    110117
    111118        /*Free ressources: */
    112119        xDelete<int>(control_type);
    113         xDelete<IssmDouble>(maxiter);
     120        xDelete<int>(maxiter);
    114121        xDelete<IssmDouble>(cm_jump);
    115122        xDelete<IssmDouble>(J);
     123        xDelete<IssmDouble>(X0);
    116124}
    117 bool controlconvergence(IssmDouble J, IssmDouble tol_cm){
    118125
    119         bool converged=false;
     126IssmDouble FormFunction(IssmDouble* X,void* usrvoid){
    120127
    121         /*Has convergence been reached?*/
    122         if (!xIsNan<IssmDouble>(tol_cm) && J<tol_cm){
    123                 converged=true;
    124                 if(VerboseConvergence()) _printf0_("      Convergence criterion reached: J = " << J << " < " << tol_cm);
    125         }
    126 
    127         return converged;
    128 }
    129 
    130 IssmDouble objectivefunction(IssmDouble search_scalar,void* optargs){
    131 
    132128        /*output: */
    133129        IssmDouble J;
    134130
    135131        /*parameters: */
    136         int        solution_type,analysis_type;
    137         bool       isFS       = false;
     132        int        solution_type,analysis_type,num_responses;
     133        bool       isFS           = false;
    138134        bool       conserve_loads = true;
    139         FemModel  *femmodel       = (FemModel*)optargs;
     135        AppCtx*    usr = (AppCtx*)usrvoid;
     136        FemModel  *femmodel  = usr->femmodel;
     137        int        nsize     = usr->nsize;
    140138
    141139        /*Recover parameters: */
    142140        femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
    143141        femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    144142        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
     143        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    145144
    146         /*set analysis type to compute velocity: */
     145        /*Constrain input vector*/
     146        IssmDouble  *XL = NULL;
     147        IssmDouble  *XU = NULL;
     148        GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
     149        GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
     150        for(long i=0;i<nsize;i++){
     151                if(X[i]>XU[i]) X[i]=XU[i];
     152                if(X[i]<XL[i]) X[i]=XL[i];
     153        }
     154
     155        /*Update control input*/
     156        SetControlInputsFromVectorx(femmodel,X);
     157
     158        /*solve forward: */
    147159        switch(solution_type){
    148160                case SteadystateSolutionEnum:
     161                        femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
     162                        stressbalance_core(femmodel);   //We need a 3D velocity!! (vz is required for the next thermal run)
     163                        break;
    149164                case StressbalanceSolutionEnum:
    150165                        femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
     166                        solutionsequence_nonlinear(femmodel,conserve_loads);
    151167                        break;
    152168                case BalancethicknessSolutionEnum:
    153169                        femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
     170                        solutionsequence_linear(femmodel);
    154171                        break;
    155172                case BalancethicknessSoftSolutionEnum:
    156                         femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
     173                        /*NOTHING*/
    157174                        break;
    158175                case Balancethickness2SolutionEnum:
    159176                        femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
     177                        solutionsequence_linear(femmodel);
    160178                        break;
    161179                default:
    162180                        _error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
    163181        }
    164182
    165         /*update parameter according to scalar: */ //false means: do not save control
    166         InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,false);
     183        /*Compute misfit for this velocity field.*/
     184        IssmDouble* Jlist = NULL;
     185        femmodel->CostFunctionx(&J,&Jlist,NULL);
     186        //_printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
    167187
    168         /*Run stressbalance with updated inputs: */
    169         if (solution_type==SteadystateSolutionEnum){
    170                 stressbalance_core(femmodel);   //We need a 3D velocity!! (vz is required for the next thermal run)
     188        /*Retrieve objective functions independently*/
     189        //for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
     190        //_printf0_("\n");
     191
     192        /*Free ressources:*/
     193        xDelete<IssmDouble>(XU);
     194        xDelete<IssmDouble>(XL);
     195        xDelete<IssmDouble>(Jlist);
     196        return J;
     197}
     198IssmDouble FormFunctionGradient(IssmDouble** pG,IssmDouble* X,void* usrvoid){
     199
     200        /*output: */
     201        IssmDouble J;
     202
     203        /*parameters: */
     204        void (*adjointcore)(FemModel*)=NULL;
     205        int         solution_type,analysis_type,num_responses,num_controls,numvertices;
     206        bool        isFS           = false;
     207        bool        conserve_loads = true;
     208        IssmDouble *scalar_list    = NULL;
     209        IssmDouble *Jlist          = NULL;
     210        IssmDouble *G              = NULL;
     211        IssmDouble *norm_list      = NULL;
     212        AppCtx     *usr            = (AppCtx*)usrvoid;
     213        FemModel   *femmodel       = usr->femmodel;
     214        int         nsize          = usr->nsize;
     215
     216        /*Recover parameters: */
     217        femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
     218        femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     219        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
     220        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
     221        femmodel->parameters->FindParam(&scalar_list,NULL,NULL,InversionGradientScalingEnum);
     222        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);   _assert_(num_controls);
     223        numvertices=femmodel->vertices->NumberOfVertices();
     224
     225        /*Constrain input vector*/
     226        IssmDouble  *XL = NULL;
     227        IssmDouble  *XU = NULL;
     228        GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
     229        GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
     230        for(long i=0;i<nsize;i++){
     231                if(X[i]>XU[i]) X[i]=XU[i];
     232                if(X[i]<XL[i]) X[i]=XL[i];
    171233        }
    172         else if (solution_type==StressbalanceSolutionEnum){
    173                 solutionsequence_nonlinear(femmodel,conserve_loads);
     234
     235        /*Update control input*/
     236        SetControlInputsFromVectorx(femmodel,X);
     237
     238        /*Compute Adjoint*/
     239        AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
     240        adjointcore(femmodel);
     241
     242        /*Compute gradient*/
     243        Gradjx(&G,&norm_list,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     244
     245        /*Compute scaling factor*/
     246        IssmDouble scalar = scalar_list[0]/norm_list[0];
     247        for(int i=1;i<num_controls;i++) scalar=min(scalar,scalar_list[i]/norm_list[i]);
     248
     249        /*Constrain Gradient*/
     250        for(int i=0;i<num_controls;i++){
     251                for(int j=0;j<numvertices;j++){
     252                        G[i*numvertices+j] = scalar*G[i*numvertices+j];
     253                }
    174254        }
    175         else if (solution_type==BalancethicknessSolutionEnum){
    176                 solutionsequence_linear(femmodel);
     255
     256        /*Needed for output results (FIXME: should be placed 6 lines below)*/
     257        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
     258
     259        for(long i=0;i<nsize;i++){
     260                if(X[i]>=XU[i]) G[i]=0.;
     261                if(X[i]<=XL[i]) G[i]=0.;
    177262        }
    178         else if (solution_type==Balancethickness2SolutionEnum){
    179                 solutionsequence_linear(femmodel);
     263
     264        /*solve forward: (FIXME: not needed actually...)*/
     265        switch(solution_type){
     266                case SteadystateSolutionEnum:
     267                        femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
     268                        stressbalance_core(femmodel);   //We need a 3D velocity!! (vz is required for the next thermal run)
     269                        break;
     270                case StressbalanceSolutionEnum:
     271                        femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
     272                        solutionsequence_nonlinear(femmodel,conserve_loads);
     273                        break;
     274                case BalancethicknessSolutionEnum:
     275                        femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
     276                        solutionsequence_linear(femmodel);
     277                        break;
     278                case BalancethicknessSoftSolutionEnum:
     279                        /*NOTHING*/
     280                        break;
     281                case Balancethickness2SolutionEnum:
     282                        femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
     283                        solutionsequence_linear(femmodel);
     284                        break;
     285                default:
     286                        _error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
    180287        }
    181         else if (solution_type==BalancethicknessSoftSolutionEnum){
    182                 /*Don't do anything*/
    183         }
    184         else{
    185                 _error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
    186         }
    187288
    188289        /*Compute misfit for this velocity field.*/
    189         femmodel->CostFunctionx(&J,NULL,NULL);
     290        femmodel->CostFunctionx(&J,&Jlist,NULL);
     291        //_printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
     292        //for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
     293        //_printf0_("\n");
    190294
    191         /*Free ressources:*/
     295        /*Clean-up and return*/
     296        xDelete<IssmDouble>(XU);
     297        xDelete<IssmDouble>(XL);
     298        xDelete<IssmDouble>(norm_list);
     299        xDelete<IssmDouble>(scalar_list);
     300        xDelete<IssmDouble>(Jlist);
     301        *pG = G;
    192302        return J;
    193303}

  • ../trunk-jpl/src/c/cores/controltao_core.cpp

     
    9292        TaoGetSolutionVector(tao,&X->pvector->vector);
    9393        G=new Vector<IssmDouble>(0); VecFree(&G->pvector->vector);
    9494        TaoGetGradientVector(tao,&G->pvector->vector);
    95         SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
     95        SetControlInputsFromVectorx(femmodel,X);
    9696        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
    9797        femmodel->OutputControlsx(&femmodel->results);
    9898        femmodel->results->AddObject(new GenericExternalResult<double*>(femmodel->results->Size()+1,JEnum,user.J,maxiter+3,1,1,0));
    … …  
    112112        TaoDestroy(&tao);
    113113        TaoFinalize();
    114114}
    115 int FormFunctionGradient(TaoSolver tao, Vec Xpetsc, IssmDouble *fcn,Vec G,void *userCtx){
     115int FormFunctionGradient(TaoSolver tao, Vec Xpetsc, IssmDouble *fcn,Vec G,void *uservoid){
    116116
    117117        /*Retreive arguments*/
    118118        int                  solution_type;
    119         AppCtx              *user            = (AppCtx *)userCtx;
     119        AppCtx              *user            = (AppCtx *)uservoid;
    120120        FemModel            *femmodel        = user->femmodel;
    121121        Vector<IssmDouble>  *gradient        = NULL;
    122122        Vector<IssmDouble>  *X               = NULL;
    … …  
    126126
    127127        /*Set new variable*/
    128128        //VecView(X,PETSC_VIEWER_STDOUT_WORLD);
    129         SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
     129        SetControlInputsFromVectorx(femmodel,X);
    130130        delete X;
    131131
    132132        /*Recover some parameters*/

  • ../trunk-jpl/src/c/classes/Elements/Tria.cpp

     
    983983                this->inputs->AddInput(new IntInput(ApproximationEnum,reCast<int>(iomodel->Data(FlowequationElementEquationEnum)[index])));
    984984        }
    985985
    986 
    987986        /*Control Inputs*/
    988987        if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
    989988                for(i=0;i<num_control_type;i++){

  • ../trunk-jpl/src/c/classes/Elements/Penta.cpp

     
    28502850        GradientIndexing(&vertexpidlist[0],control_index);
    28512851
    28522852        /*Get input (either in element or material)*/
     2853        if(control_enum==MaterialsRheologyBbarEnum) control_enum=MaterialsRheologyBEnum;
    28532854        Input* input=inputs->GetInput(control_enum);
    28542855        if(!input) _error_("Input " << EnumToStringx(control_enum) << " not found in element");
    28552856
    … …  
    28642865void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
    28652866
    28662867        IssmDouble  values[NUMVERTICES];
    2867         int     vertexpidlist[NUMVERTICES];
    2868         Input  *input     = NULL;
    2869         Input  *new_input = NULL;
     2868        int         vertexpidlist[NUMVERTICES],control_init;
     2869        Input      *input     = NULL;
     2870        Input      *new_input = NULL;
    28702871
     2872        /*Specific case for depth averaged quantities*/
     2873        control_init=control_enum;
     2874        if(control_enum==MaterialsRheologyBbarEnum){
     2875                control_enum=MaterialsRheologyBEnum;
     2876                if(!IsOnBase()) return;
     2877        }
     2878        if(control_enum==DamageDbarEnum){
     2879                control_enum=DamageDEnum;
     2880                if(!IsOnBase()) return;
     2881        }
     2882
    28712883        /*Get out if this is not an element input*/
    28722884        if(!IsInput(control_enum)) return;
    28732885
    … …  
    28752887        GradientIndexing(&vertexpidlist[0],control_index);
    28762888
    28772889        /*Get values on vertices*/
    2878         for (int i=0;i<NUMVERTICES;i++){
     2890        for(int i=0;i<NUMVERTICES;i++){
    28792891                values[i]=vector[vertexpidlist[i]];
    28802892        }
    28812893        new_input = new PentaInput(control_enum,values,P1Enum);
    … …  
    28862898        }
    28872899
    28882900        ((ControlInput*)input)->SetInput(new_input);
     2901
     2902        if(control_init==MaterialsRheologyBbarEnum){
     2903                this->InputExtrude(control_enum);
     2904        }
     2905        if(control_init==DamageDbarEnum){
     2906                this->InputExtrude(control_enum);
     2907        }
    28892908}
    28902909/*}}}*/
    28912910

  • ../trunk-jpl/src/m/classes/inversion.m

     
    2424                thickness_obs               = NaN
    2525        end
    2626        methods
    27          function createxml(obj,fid) % {{{
    28             fprintf(fid, '<!-- inversion -->\n');           
    29                    
    30             % inversion parameters
    31             fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="inversion parameters">','<section name="inversion" />');                   
    32             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="iscontrol" type="',class(obj.iscontrol),'" default="',convert2str(obj.iscontrol),'">','     <section name="inversion" />','     <help> is inversion activated? </help>','  </parameter>');
    33            
    34             % incompleteadjoing drop-down (0 or 1)
    35             fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="incomplete_adjoint" type="alternative" optional="false">','     <section name="inversion" />','     <help> 1: linear viscosity, 0: non-linear viscosity </help>');
    36             fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
    37             fprintf(fid,'%s\n%s\n','       <option value="1" type="string" default="false"> </option>','</parameter>');
    38            
    39             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="control_parameters" type="',class(obj.control_parameters),'" default="',convert2str(obj.control_parameters),'">','     <section name="inversion" />','     <help> ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''} </help>','  </parameter>');
    40                
    41             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="nsteps" type="',class(obj.nsteps),'" default="',convert2str(obj.nsteps),'">','     <section name="inversion" />','     <help> number of optimization searches </help>','  </parameter>');
    42             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions" type="',class(obj.cost_functions),'" default="',convert2str(obj.cost_functions),'">','     <section name="inversion" />','     <help> indicate the type of response for each optimization step  </help>','  </parameter>');
    43             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions_coefficients" type="',class(obj.cost_functions_coefficients),'" default="',convert2str(obj.cost_functions_coefficients),'">','     <section name="inversion" />','     <help> cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter </help>','  </parameter>');
    44                
    45             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_function_threshold" type="',class(obj.cost_function_threshold),'" default="',convert2str(obj.cost_function_threshold),'">','     <section name="inversion" />','     <help> misfit convergence criterion. Default is 1%, NaN if not applied </help>','  </parameter>');
    46             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="maxiter_per_step" type="',class(obj.maxiter_per_step),'" default="',convert2str(obj.maxiter_per_step),'">','     <section name="inversion" />','     <help> maximum iterations during each optimization step  </help>','  </parameter>');
    47             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="gradient_scaling" type="',class(obj.gradient_scaling),'" default="',convert2str(obj.gradient_scaling),'">','     <section name="inversion" />','     <help> scaling factor on gradient direction during optimization, for each optimization step </help>','  </parameter>');
    48                
    49             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="step_threshold" type="',class(obj.step_threshold),'" default="',convert2str(obj.step_threshold),'">','     <section name="inversion" />','     <help> decrease threshold for misfit, default is 30% </help>','  </parameter>');
    50             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="min_parameters" type="',class(obj.min_parameters),'" default="',convert2str(obj.min_parameters),'">','     <section name="inversion" />','     <help> absolute minimum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
    51             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="max_parameters" type="',class(obj.max_parameters),'" default="',convert2str(obj.max_parameters),'">','     <section name="inversion" />','     <help> absolute maximum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
    52                
    53             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vx_obs" type="',class(obj.vx_obs),'" default="',convert2str(obj.vx_obs),'">','     <section name="inversion" />','     <help> observed velocity x component [m/yr] </help>','  </parameter>');
    54             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vy_obs" type="',class(obj.vy_obs),'" default="',convert2str(obj.vy_obs),'">','     <section name="inversion" />','     <help> observed velocity y component [m/yr]  </help>','  </parameter>');
    55             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vel_obs" type="',class(obj.vel_obs),'" default="',convert2str(obj.vel_obs),'">','     <section name="inversion" />','     <help> observed velocity magnitude [m/yr] </help>','  </parameter>');
    56             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="thickness_obs" type="',class(obj.thickness_obs),'" default="',convert2str(obj.thickness_obs),'">','     <section name="inversion" />','     <help> observed thickness [m]) </help>','  </parameter>');
    57                
    58             fprintf(fid,'%s\n%s\n','</frame>');   
    59            
    60             fprintf(fid,'%s\n%s\n%s\n','<frame key="2" label="Available cost functions">','<section name="inversion" />');                   
    61             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAbsVelMisfit" type="','string','" default="','101','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
    62             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceRelVelMisfit" type="','string','" default="','102','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
    63             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVelMisfit" type="','string','" default="','103','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
    64                
    65             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVxVyMisfit" type="','string','" default="','104','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
    66             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAverageVelMisfit" type="','string','" default="','105','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
    67             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsMisfit" type="','string','" default="','106','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
    68                
    69             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="DragCoefficientAbsGradient" type="','string','" default="','107','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
    70             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="RheologyBbarAbsGradient" type="','string','" default="','108','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
    71             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsGradient" type="','string','" default="','109','">','     <section name="inversion" />','     <help> </help>','  </parameter>');
    72                
    73             fprintf(fid,'%s\n%s\n','</frame>');   
    74        
    75         end % }}}       
     27                function createxml(obj,fid) % {{{
     28                        fprintf(fid, '<!-- inversion -->\n');           
     29
     30                        % inversion parameters
     31                        fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="inversion parameters">','<section name="inversion" />');                   
     32                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="iscontrol" type="',class(obj.iscontrol),'" default="',convert2str(obj.iscontrol),'">','     <section name="inversion" />','     <help> is inversion activated? </help>','  </parameter>');
     33
     34                        % incompleteadjoing drop-down (0 or 1)
     35                        fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="incomplete_adjoint" type="alternative" optional="false">','     <section name="inversion" />','     <help> 1: linear viscosity, 0: non-linear viscosity </help>');
     36                        fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
     37                        fprintf(fid,'%s\n%s\n','       <option value="1" type="string" default="false"> </option>','</parameter>');
     38
     39                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="control_parameters" type="',class(obj.control_parameters),'" default="',convert2str(obj.control_parameters),'">','     <section name="inversion" />','     <help> ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''} </help>','  </parameter>');
     40
     41                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="nsteps" type="',class(obj.nsteps),'" default="',convert2str(obj.nsteps),'">','     <section name="inversion" />','     <help> number of optimization searches </help>','  </parameter>');
     42                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions" type="',class(obj.cost_functions),'" default="',convert2str(obj.cost_functions),'">','     <section name="inversion" />','     <help> indicate the type of response for each optimization step  </help>','  </parameter>');
     43                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions_coefficients" type="',class(obj.cost_functions_coefficients),'" default="',convert2str(obj.cost_functions_coefficients),'">','     <section name="inversion" />','     <help> cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter </help>','  </parameter>');
     44
     45                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_function_threshold" type="',class(obj.cost_function_threshold),'" default="',convert2str(obj.cost_function_threshold),'">','     <section name="inversion" />','     <help> misfit convergence criterion. Default is 1%, NaN if not applied </help>','  </parameter>');
     46                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="maxiter_per_step" type="',class(obj.maxiter_per_step),'" default="',convert2str(obj.maxiter_per_step),'">','     <section name="inversion" />','     <help> maximum iterations during each optimization step  </help>','  </parameter>');
     47                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="gradient_scaling" type="',class(obj.gradient_scaling),'" default="',convert2str(obj.gradient_scaling),'">','     <section name="inversion" />','     <help> scaling factor on gradient direction during optimization, for each optimization step </help>','  </parameter>');
     48
     49                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="step_threshold" type="',class(obj.step_threshold),'" default="',convert2str(obj.step_threshold),'">','     <section name="inversion" />','     <help> decrease threshold for misfit, default is 30% </help>','  </parameter>');
     50                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="min_parameters" type="',class(obj.min_parameters),'" default="',convert2str(obj.min_parameters),'">','     <section name="inversion" />','     <help> absolute minimum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
     51                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="max_parameters" type="',class(obj.max_parameters),'" default="',convert2str(obj.max_parameters),'">','     <section name="inversion" />','     <help> absolute maximum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
     52
     53                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vx_obs" type="',class(obj.vx_obs),'" default="',convert2str(obj.vx_obs),'">','     <section name="inversion" />','     <help> observed velocity x component [m/yr] </help>','  </parameter>');
     54                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vy_obs" type="',class(obj.vy_obs),'" default="',convert2str(obj.vy_obs),'">','     <section name="inversion" />','     <help> observed velocity y component [m/yr]  </help>','  </parameter>');
     55                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vel_obs" type="',class(obj.vel_obs),'" default="',convert2str(obj.vel_obs),'">','     <section name="inversion" />','     <help> observed velocity magnitude [m/yr] </help>','  </parameter>');
     56                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="thickness_obs" type="',class(obj.thickness_obs),'" default="',convert2str(obj.thickness_obs),'">','     <section name="inversion" />','     <help> observed thickness [m]) </help>','  </parameter>');
     57
     58                        fprintf(fid,'%s\n%s\n','</frame>');   
     59
     60                        fprintf(fid,'%s\n%s\n%s\n','<frame key="2" label="Available cost functions">','<section name="inversion" />');                   
     61                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAbsVelMisfit" type="','string','" default="','101','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
     62                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceRelVelMisfit" type="','string','" default="','102','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
     63                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVelMisfit" type="','string','" default="','103','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
     64
     65                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVxVyMisfit" type="','string','" default="','104','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
     66                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAverageVelMisfit" type="','string','" default="','105','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
     67                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsMisfit" type="','string','" default="','106','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
     68
     69                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="DragCoefficientAbsGradient" type="','string','" default="','107','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
     70                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="RheologyBbarAbsGradient" type="','string','" default="','108','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
     71                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsGradient" type="','string','" default="','109','">','     <section name="inversion" />','     <help> </help>','  </parameter>');
     72
     73                        fprintf(fid,'%s\n%s\n','</frame>');   
     74
     75                end % }}}       
    7676                function obj = inversion(varargin) % {{{
    7777                        switch nargin
    7878                                case 0
    … …  
    195195                        WriteData(fid,'object',obj,'fieldname','incomplete_adjoint','format','Boolean');
    196196                        if ~obj.iscontrol, return; end
    197197                        WriteData(fid,'object',obj,'fieldname','nsteps','format','Integer');
    198                         WriteData(fid,'object',obj,'fieldname','maxiter_per_step','format','DoubleMat','mattype',3);
     198                        WriteData(fid,'object',obj,'fieldname','maxiter_per_step','format','IntMat','mattype',3);
    199199                        WriteData(fid,'object',obj,'fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1);
    200200                        WriteData(fid,'object',obj,'fieldname','gradient_scaling','format','DoubleMat','mattype',3);
    201201                        WriteData(fid,'object',obj,'fieldname','cost_function_threshold','format','Double');
Note: See TracBrowser for help on using the repository browser.