Changeset 6226


Ignore:
Timestamp:
10/11/10 10:05:52 (15 years ago)
Author:
Mathieu Morlighem
Message:

Multivariable optimization now available (still a problem to be fixed in Gradient)

Location:
issm/trunk/src
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp

    r6213 r6226  
    4343                parameters->AddObject(new DoubleParam(MeanVelEnum,iomodel->meanvel));
    4444                parameters->AddObject(new DoubleParam(CmNoiseDmpEnum,iomodel->cm_noisedmp));
    45                 parameters->AddObject(new DoubleParam(CmMinEnum,iomodel->cm_min));
    46                 parameters->AddObject(new DoubleParam(CmMaxEnum,iomodel->cm_max));
     45
    4746                parameters->AddObject(new BoolParam(CmGradientEnum,iomodel->cm_gradient));
    4847
    … …  
    5049                IoModelFetchData(&iomodel->cm_responses,NULL,NULL,iomodel_handle,"cm_responses");
    5150                IoModelFetchData(&iomodel->cm_jump,NULL,NULL,iomodel_handle,"cm_jump");
     51                IoModelFetchData(&iomodel->cm_min,NULL,NULL,iomodel_handle,"cm_min");
     52                IoModelFetchData(&iomodel->cm_max,NULL,NULL,iomodel_handle,"cm_max");
    5253                IoModelFetchData(&iomodel->optscal,NULL,NULL,iomodel_handle,"optscal");
    5354                IoModelFetchData(&iomodel->maxiter,NULL,NULL,iomodel_handle,"maxiter");
    … …  
    5556                parameters->AddObject(new DoubleVecParam(CmResponsesEnum,iomodel->cm_responses,iomodel->nsteps));
    5657                parameters->AddObject(new DoubleVecParam(CmJumpEnum,iomodel->cm_jump,iomodel->nsteps));
    57                 parameters->AddObject(new DoubleVecParam(OptScalEnum,iomodel->optscal,iomodel->nsteps));
     58                parameters->AddObject(new DoubleMatParam(OptScalEnum,iomodel->optscal,iomodel->nsteps,iomodel->num_control_type));
     59                parameters->AddObject(new DoubleVecParam(CmMinEnum,iomodel->cm_min,iomodel->num_control_type));
     60                parameters->AddObject(new DoubleVecParam(CmMaxEnum,iomodel->cm_max,iomodel->num_control_type));
    5861                parameters->AddObject(new DoubleVecParam(MaxIterEnum,iomodel->maxiter,iomodel->nsteps));
    5962
    6063                xfree((void**)&iomodel->cm_responses);
    6164                xfree((void**)&iomodel->cm_jump);
     65                xfree((void**)&iomodel->cm_min);
     66                xfree((void**)&iomodel->cm_max);
    6267                xfree((void**)&iomodel->optscal);
    6368                xfree((void**)&iomodel->maxiter);

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

    r6222 r6226  
    11341134        int*   control_type=NULL;
    11351135        Input* input=NULL;
    1136         double cm_min,cm_max;
     1136        double *cm_min=NULL;
     1137        double *cm_max=NULL;
    11371138
    11381139        /*retrieve some parameters: */
    1139         this->parameters->FindParam(&cm_min,CmMinEnum);
    1140         this->parameters->FindParam(&cm_max,CmMaxEnum);
     1140        this->parameters->FindParam(&cm_min,NULL,CmMinEnum);
     1141        this->parameters->FindParam(&cm_max,NULL,CmMaxEnum);
    11411142        this->parameters->FindParam(&num_controls,NumControlsEnum);
    11421143        this->parameters->FindParam(&control_type,NULL,ControlTypeEnum);
    … …  
    11551156
    11561157                ((ControlInput*)input)->UpdateValue(scalar);
    1157                 input->Constrain(cm_min,cm_max);
     1158                input->Constrain(cm_min[i],cm_max[i]);
    11581159                if (save_parameter) ((ControlInput*)input)->SaveValue();
    11591160

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

    r6222 r6226  
    13421342        int*   control_type=NULL;
    13431343        Input* input=NULL;
    1344         double cm_min,cm_max;
     1344        double *cm_min=NULL;
     1345        double *cm_max=NULL;
    13451346
    13461347        /*retrieve some parameters: */
    1347         this->parameters->FindParam(&cm_min,CmMinEnum);
    1348         this->parameters->FindParam(&cm_max,CmMaxEnum);
     1348        this->parameters->FindParam(&cm_min,NULL,CmMinEnum);
     1349        this->parameters->FindParam(&cm_max,NULL,CmMaxEnum);
    13491350        this->parameters->FindParam(&num_controls,NumControlsEnum);
    13501351        this->parameters->FindParam(&control_type,NULL,ControlTypeEnum);
    … …  
    13641365
    13651366                ((ControlInput*)input)->UpdateValue(scalar);
    1366                 input->Constrain(cm_min,cm_max);
     1367                input->Constrain(cm_min[i],cm_max[i]);
    13671368                if (save_parameter) ((ControlInput*)input)->SaveValue();
    13681369

  • issm/trunk/src/c/objects/IoModel.cpp

    r6213 r6226  
    8888        xfree((void**)&this->weights);
    8989        xfree((void**)&this->cm_jump);
     90        xfree((void**)&this->cm_min);
     91        xfree((void**)&this->cm_max);
    9092        xfree((void**)&this->optscal);
    9193        xfree((void**)&this->maxiter);
    … …  
    166168        IoModelFetchData(&this->tolx,iomodel_handle,"tolx");
    167169        IoModelFetchData(&this->cm_noisedmp,iomodel_handle,"cm_noisedmp");
    168         IoModelFetchData(&this->cm_min,iomodel_handle,"cm_min");
    169         IoModelFetchData(&this->cm_max,iomodel_handle,"cm_max");
    170170        IoModelFetchData(&this->cm_gradient,iomodel_handle,"cm_gradient");
    171171        IoModelFetchData(&this->eps_res,iomodel_handle,"eps_res");
    … …  
    306306        this->maxiter=NULL;
    307307        this->cm_noisedmp=0;
    308         this->cm_min=0;
    309         this->cm_max=0;
     308        this->cm_min=NULL;
     309        this->cm_max=NULL;
    310310        this->cm_gradient=0;
    311311        this->verbose=0;

  • issm/trunk/src/c/objects/IoModel.h

    r6213 r6226  
    114114                double  stokesreconditioning;
    115115                double  cm_noisedmp;
    116                 double  cm_min;
    117                 double  cm_max;
     116                double* cm_min;
     117                double* cm_max;
    118118                int     cm_gradient;;
    119119

  • issm/trunk/src/c/objects/Params/DoubleParam.cpp

    r6213 r6226  
    167167}
    168168/*}}}*/
     169/*FUNCTION DoubleParam::GetParameterValue(double** pdoublearray,int* pM){{{1*/
     170void DoubleParam::GetParameterValue(double** pdoublearray,int* pM){
     171#ifdef _SERIAL_
     172        double* output=NULL;
     173
     174        output=(double*)xmalloc(1*sizeof(double));
     175        *output=(double)value;
     176
     177        /*Assign output podoubleers:*/
     178        if(pM) *pM=1;
     179        *pdoublearray=output;
     180#else
     181        ISSMERROR("Double param of enum %i (%s) cannot return an array of double",enum_type,EnumToString(enum_type));
     182#endif
     183}
     184/*}}}*/
     185/*FUNCTION DoubleParam::GetParameterValue(double** pdoublearray,int* pM,int* pN){{{1*/
     186void DoubleParam::GetParameterValue(double** pdoublearray,int* pM,int* pN){
     187#ifdef _SERIAL_
     188        double* output=NULL;
     189
     190        output=(double*)xmalloc(1*sizeof(double));
     191        *output=(double)value;
     192
     193        /*Assign output podoubleers:*/
     194        if(pM) *pM=1;
     195        if(pN) *pN=1;
     196        *pdoublearray=output;
     197#else
     198        ISSMERROR("Double param of enum %i (%s) cannot return an array of double",enum_type,EnumToString(enum_type));
     199#endif
     200}
     201/*}}}*/
    169202/*FUNCTION DoubleParam::SetMatlabField{{{1*/
    170203#ifdef _SERIAL_

  • issm/trunk/src/c/objects/Params/DoubleParam.h

    r6213 r6226  
    5656                void  GetParameterValue(char** pstring){ISSMERROR("Double param of enum %i (%s) cannot return a string",enum_type,EnumToString(enum_type));}
    5757                void  GetParameterValue(char*** pstringarray,int* pM){ISSMERROR("Double param of enum %i (%s) cannot return a string arrayl",enum_type,EnumToString(enum_type));}
    58                 void  GetParameterValue(double** pdoublearray,int* pM){ISSMERROR("Double param of enum %i (%s) cannot return a double array",enum_type,EnumToString(enum_type));}
    59                 void  GetParameterValue(double** pdoublearray,int* pM, int* pN){ISSMERROR("Double param of enum %i (%s) cannot return a double array",enum_type,EnumToString(enum_type));}
     58                void  GetParameterValue(double** pdoublearray,int* pM);
     59                void  GetParameterValue(double** pdoublearray,int* pM, int* pN);
    6060                void  GetParameterValue(double*** parray, int* pM,int** pmdims, int** pndims){ISSMERROR("Double param of enum %i (%s) cannot return a matrix array",enum_type,EnumToString(enum_type));}
    6161                void  GetParameterValue(Vec* pvec){ISSMERROR("Double param of enum %i (%s) cannot return a Vec",enum_type,EnumToString(enum_type));}

  • issm/trunk/src/c/solutions/control_core.cpp

    r6213 r6226  
    3030        int*    control_type = NULL;
    3131        double* responses=NULL;
    32         double* optscal=NULL;
    3332        double* maxiter=NULL;
    3433        double* cm_jump=NULL;
    … …  
    5251        femmodel->parameters->FindParam(&control_type,NULL,ControlTypeEnum);
    5352        femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum);
    54         femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum);
    5553        femmodel->parameters->FindParam(&maxiter,NULL,MaxIterEnum);
    5654        femmodel->parameters->FindParam(&cm_jump,NULL,CmJumpEnum);
    … …  
    107105
    108106                _printf_("%s\n","      updating parameter using optimized search scalar..."); //true means update save controls
    109                 InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar*optscal[n],true);
     107                InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,true);
    110108               
    111109                if(controlconvergence(J,responses,eps_cm,n)) break;
    … …  
    133131        xfree((void**)&control_type);
    134132        xfree((void**)&responses);
    135         xfree((void**)&optscal);
    136133        xfree((void**)&maxiter);
    137134        xfree((void**)&cm_jump);

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

    r6213 r6226  
    1616       
    1717        /*parameters: */
    18         bool control_steady;
    19         int  num_controls;
    20         int* control_type=NULL;
    21         int  verbose;
     18        bool    control_steady;
     19        int     verbose;
     20        int     num_controls;
     21        int    *control_type   = NULL;
     22        double *optscal_list   = NULL;
    2223
    2324        /*Intermediaries*/
    … …  
    3031        femmodel->parameters->FindParam(&num_controls,NumControlsEnum);
    3132        femmodel->parameters->FindParam(&control_type,NULL,ControlTypeEnum);
     33        femmodel->parameters->FindParam(&optscal_list,NULL,NULL,OptScalEnum);
    3234        femmodel->parameters->FindParam(&verbose,VerboseEnum);
    3335
    … …  
    5153                VecFree(&gradient);
    5254
     55                /*Scale gradient for current step and current parameter*/
     56                VecScale(new_gradient,optscal_list[num_controls*step+i]);
     57
    5358                /*plug back into inputs: */
    5459                ControlInputSetGradientx(femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmodel->parameters,control_type[i],new_gradient);

  • issm/trunk/src/c/solutions/objectivefunctionC.cpp

    r6213 r6226  
    3030        FemModel  *femmodel       = NULL;
    3131        int        n;
    32         double    *optscal        = NULL;
    3332        double    *responses      = NULL;
    3433        int        solution_type;
    … …  
    4342        n=optargs->n;
    4443
    45         femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum);
    4644        femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum);
    4745        femmodel->parameters->FindParam(&isstokes,IsStokesEnum);
    … …  
    6159
    6260        /*update parameter according to scalar: */ //false means: do not save control
    63         InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar*optscal[n],false);
     61        InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,false);
    6462
    6563        /*Run diagnostic with updated inputs: */
    … …  
    8280        /*Free ressources:*/
    8381        xfree((void**)&responses);
    84         xfree((void**)&optscal);
    8582        return J;
    8683}

  • issm/trunk/src/m/model/ismodelselfconsistent.m

    r6176 r6226  
    170170end
    171171%}}}
    172 %SCALAR {{{
    173 if ~isscalar(md.control_type),
    174         error('model not consistent: md.control_type should be a scalar');
    175 end
    176 %}}}
    177172%PARAMETEROUTPUT {{{1
    178173if md.numoutput~=length(md.parameteroutput),
    … …  
    344339
    345340        %CONTROL TYPE
    346         if ~isnumeric(md.control_type),
    347                 error('model not consistent: control_type should be an enum');
    348         end
     341        num_controls=numel(md.control_type);
    349342        checkvalues(md,{'control_type'},[DhDtEnum DragCoefficientEnum RheologyBbarEnum VxEnum VyEnum]);
    350343
    351344        %LENGTH CONTROL FIELDS
    352         fields={'maxiter','optscal','cm_responses','cm_jump'};
     345        fields={'maxiter','cm_responses','cm_jump'};
    353346        checksize(md,fields,[md.nsteps 1]);
     347        fields={'optscal'};
     348        checksize(md,fields,[md.nsteps num_controls]);
     349        fields={'cm_min','cm_max'};
     350        checksize(md,fields,[num_controls 1]);
    354351
    355352        %RESPONSES

  • issm/trunk/src/m/model/marshall.m

    r6214 r6226  
    119119WriteData(fid,md.maxiter,'Mat','maxiter');
    120120WriteData(fid,md.cm_noisedmp,'Scalar','cm_noisedmp');
    121 WriteData(fid,md.cm_min,'Scalar','cm_min');
    122 WriteData(fid,md.cm_max,'Scalar','cm_max');
     121WriteData(fid,md.cm_min,'Mat','cm_min');
     122WriteData(fid,md.cm_max,'Mat','cm_max');
    123123WriteData(fid,md.cm_gradient,'Integer','cm_gradient');
    124124WriteData(fid,md.eps_res,'Scalar','eps_res');

  • issm/trunk/src/m/solutions/control_core.m

    r6214 r6226  
    1313        solution_type=femmodel.parameters.SolutionType;
    1414        responses=femmodel.parameters.CmResponses;
    15         optscal=femmodel.parameters.OptScal;
    1615        maxiter=femmodel.parameters.MaxIter;
    1716        cm_jump=femmodel.parameters.CmJump;
    … …  
    6665
    6766                displaystring('\n%s',['      updating parameter using optimized search scalar:']);
    68                 [femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlUpdate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,search_scalar*optscal(n),1);
     67                [femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlUpdate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,search_scalar,1);
    6968
    7069                disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(J(n))]);

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

    r6214 r6226  
    2525        control_type=femmodel.parameters.ControlType;
    2626        control_steady=femmodel.parameters.ControlSteady;
     27        optscal=femmodel.parameters.OptScal;
    2728
    2829        for i=1:num_controls,
    … …  
    4546                end
    4647
     48                %Scale Gradient
     49                new_gradient=optscal(step,i)*new_gradient;
     50
    4751                %plug back into inputs:
    4852                [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);

  • issm/trunk/src/m/solutions/objectivefunctionC.m

    r6214 r6226  
    44conserve_loads=true;
    55%recover some parameters
    6 optscal=femmodel.parameters.OptScal(n);
    76response=femmodel.parameters.CmResponses(n);
    87analysis_type=femmodel.parameters.AnalysisType;
    … …  
    2019
    2120%Use search scalar to shoot parameter in the gradient direction:
    22 [femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlUpdate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,search_scalar*optscal,0);
     21[femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlUpdate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,search_scalar,0);
    2322
    2423%Run diagnostic with updated inputs:
Note: See TracChangeset for help on using the changeset viewer.