source: issm/oecreview/Archive/18296-19100/ISSM-18615-18616.diff@ 19102

Last change on this file since 19102 was 19102, checked in by Mathieu Morlighem, 10 years ago

NEW: added 18296-19100
File size: 5.8 KB

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

     
    1212        bool        control_analysis;
    1313        int         inversiontype;
    1414        int         nsteps;
    15         int         num_control_type;
    16         int         num_cm_responses;
     15        int         num_controls;
     16        int         num_costfunc;
    1717        int        *control_type     = NULL;
    1818        int        *maxiter          = NULL;
    1919        int        *cm_responses     = NULL;
    2020        IssmDouble *cm_jump          = NULL;
    2121        IssmDouble *optscal          = NULL;
     22        IssmDouble *control_scaling_factors = NULL;
    2223
    2324        /*retrieve some parameters: */
    2425        iomodel->Constant(&control_analysis,InversionIscontrolEnum);
    … …  
    4041                }
    4142
    4243                /*Now, recover fit, optscal and maxiter as vectors: */
    43                 iomodel->FetchData(&control_type,NULL,&num_control_type,InversionControlParametersEnum);
    44                 iomodel->FetchData(&cm_responses,NULL,&num_cm_responses,InversionCostFunctionsEnum);
    45                 parameters->AddObject(new IntVecParam(InversionControlParametersEnum,control_type,num_control_type));
    46                 parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,cm_responses,num_cm_responses));
     44                iomodel->FetchData(&control_type,NULL,&num_controls,InversionControlParametersEnum);
     45                iomodel->FetchData(&cm_responses,NULL,&num_costfunc,InversionCostFunctionsEnum);
     46                parameters->AddObject(new IntVecParam(InversionControlParametersEnum,control_type,num_controls));
     47                parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,cm_responses,num_costfunc));
    4748
    4849                /*Inversion type specifics*/
    4950                switch(inversiontype){
    … …  
    5354                                iomodel->FetchData(&cm_jump,&nsteps,NULL,InversionStepThresholdEnum);
    5455                                iomodel->FetchData(&optscal,NULL,NULL,InversionGradientScalingEnum);
    5556                                iomodel->FetchData(&maxiter,NULL,NULL,InversionMaxiterPerStepEnum);
    56                                 parameters->AddObject(new DoubleMatParam(InversionGradientScalingEnum,optscal,nsteps,num_control_type));
     57                                parameters->AddObject(new DoubleMatParam(InversionGradientScalingEnum,optscal,nsteps,num_controls));
    5758                                parameters->AddObject(new DoubleVecParam(InversionStepThresholdEnum,cm_jump,nsteps));
    5859                                parameters->AddObject(new IntVecParam(InversionMaxiterPerStepEnum,maxiter,nsteps));
    5960                                break;
    … …  
    7273                                parameters->AddObject(iomodel->CopyConstantObject(InversionGttolEnum));
    7374                                parameters->AddObject(iomodel->CopyConstantObject(InversionMaxstepsEnum));
    7475                                parameters->AddObject(iomodel->CopyConstantObject(InversionMaxiterEnum));
     76                                iomodel->FetchData(&control_scaling_factors,NULL,NULL,InversionControlScalingFactorsEnum);
     77                                parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_controls));
    7578                                break;
    7679                        case 3:/*Validation*/
    7780                                break;

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

     
    3232        double       f,dxmin,gttol;
    3333        int          maxsteps,maxiter;
    3434        int          intn,num_controls,solution_type;
     35        IssmDouble  *scaling_factors = NULL;
    3536        IssmDouble  *X  = NULL;
    3637        IssmDouble  *G  = NULL;
    3738
    … …  
    4243        femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum);
    4344        femmodel->parameters->FindParam(&dxmin,InversionDxminEnum);
    4445        femmodel->parameters->FindParam(&gttol,InversionGttolEnum);
     46        femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
    4547        femmodel->parameters->SetParam(false,SaveResultsEnum);
    4648
    4749        /*Initialize M1QN3 parameters*/
    … …  
    7375        long n = long(intn);
    7476        G = xNew<double>(n);
    7577
     78        /*Scale control for M1QN3*/
     79        if(num_controls!=1) _error_("not supported yet...");
     80        for(long i=0;i<n;i++){
     81                X[i] = X[i]/scaling_factors[0];
     82        }
     83
    7684        /*Allocate m1qn3 working arrays (see doc)*/
    7785        long      m   = 100;
    7886        long      ndz = 4*n+m*(2*n+1);
    … …  
    131139        FemModel   *femmodel  = (FemModel*)dzs;
    132140
    133141        /*Recover number of cost functions responses*/
    134         int num_responses;
     142        int         num_responses;
     143        int         num_controls;
     144        IssmDouble* scaling_factors = NULL;
    135145        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
     146        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     147        femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
    136148
    137149        /*Constrain input vector*/
    138150        IssmDouble  *XL = NULL;
    139151        IssmDouble  *XU = NULL;
    140152        GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
    141153        GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
     154        if(num_controls!=1) _error_("not supported yet");
    142155        for(long i=0;i<*n;i++){
     156                X[i] = X[i]*scaling_factors[0];
    143157                if(X[i]>XU[i]) X[i]=XU[i];
    144158                if(X[i]<XL[i]) X[i]=XL[i];
    145159        }
    … …  
    161175        femmodel->CostFunctionx(pf,&Jlist,NULL);
    162176        _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
    163177
    164 
    165 
    166178        if(indic==0){
    167179                /*dry run, no gradient required*/
    168180
    … …  
    189201
    190202        /*Constrain Gradient*/
    191203        IssmDouble  Gnorm = 0.;
     204        if(num_controls!=1) _error_("not supported yet");
    192205        for(long i=0;i<*n;i++){
    193206                if(X[i]>=XU[i]) G[i]=0.;
    194207                if(X[i]<=XL[i]) G[i]=0.;
     208                G[i] = G[i]*scaling_factors[0];
     209                X[i] = X[i]/scaling_factors[0];
    195210                Gnorm += G[i]*G[i];
    196211        }
    197212        Gnorm = sqrt(Gnorm);
    … …  
    202217        _printf0_("\n");
    203218
    204219        /*Clean-up and return*/
    205                 xDelete<IssmDouble>(Jlist);
     220        xDelete<IssmDouble>(Jlist);
    206221        xDelete<IssmDouble>(XU);
    207222        xDelete<IssmDouble>(XL);
    208223}
Note: See TracBrowser for help on using the repository browser.