Changeset 20653

Timestamp:

05/26/16 00:34:45 (9 years ago)

Author:

Mathieu Morlighem

Message:

CHG: testing new IoModel using jekins (no AD capability on my local machine)

Location:

issm/trunk-jpl/src/c

Files:

• Makefile.am (modified) (1 diff)
• classes/IoModel.cpp (modified) (44 diffs)
• classes/IoModel.h (modified) (4 diffs)
• main/kriging.cpp (modified) (1 diff)
• modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/Makefile.am

@@ -26 +26 @@
                                         ./classes/gauss/GaussTetra.cpp\
                                         ./classes/gauss/GaussPenta.cpp\
+                                        ./classes/IoModel.cpp\
                                         ./classes/FemModel.cpp\
                                         ./classes/Loads/Friction.cpp\
                                         ./classes/Inputs/TransientInput.cpp\
                                         ./classes/Constraints/SpcTransient.cpp\
-                                        ./classes/IndependentObject.cpp\
                                         ./classes/DependentObject.cpp\
                                         ./classes/DofIndexing.cpp\
-                                        ./classes/IoModel.cpp\
                                         ./classes/Contours.cpp\
                                         ./classes/Nodes.cpp\

issm/trunk-jpl/src/c/classes/IoModel.cpp

@@ -3 +3 @@
  * into ISSM, from Matlab, or through a binary file opened for reading.
  */
+
+/*CODES:
+ * 1: boolean constant
+ * 2: integer constant
+ * 3: IssmDouble constant
+ * 5: boolean vector
+ * 6: int vector
+ * 7: IssmDouble vector*/
 
 #ifdef HAVE_CONFIG_H
@@ -19 +27 @@
 #include "../shared/shared.h"
 
+/*NEW NEW NEW NEW NEW NEW NEW NEW NEW NEW*/
+IoConstant::IoConstant(){/*{{{*/
+        this->isindependent = false;
+        this->data_enum     = -1;
+        this->name          = NULL;
+        this->constant      = NULL;
+}
+/*}}}*/
+IoConstant::~IoConstant(){/*{{{*/
+        xDelete<char>(this->name);
+        delete this->constant;
+}
+/*}}}*/
+IoConstant::IoConstant(bool value,int enum_in){/*{{{*/
+        this->isindependent = false;
+        this->data_enum     = enum_in;
+        this->name          = NULL;
+        this->constant      = new BoolParam(enum_in,value);
+}
+/*}}}*/
+IoConstant::IoConstant(int value,int enum_in){/*{{{*/
+        this->isindependent = false;
+        this->data_enum     = enum_in;
+        this->name          = NULL;
+        this->constant      = new IntParam(enum_in,value);
+}
+/*}}}*/
+IoConstant::IoConstant(IssmDouble value,int enum_in){/*{{{*/
+        this->isindependent = false;
+        this->data_enum     = enum_in;
+        this->name          = NULL;
+        this->constant      = new DoubleParam(enum_in,value);
+}
+/*}}}*/
+IoConstant::IoConstant(char* value,int enum_in){/*{{{*/
+        this->isindependent = false;
+        this->data_enum     = enum_in;
+        this->name          = NULL;
+        this->constant      = new StringParam(enum_in,value);
+}
+/*}}}*/
+
+IoData::IoData(){/*{{{*/
+        this->isindependent = false;
+        this->data_enum     = -1;
+        this->name          = NULL;
+        this->code          = -1;
+        this->layout        = -1;
+        this->M             = 0;
+        this->N             = 0;
+        this->data          = NULL;
+}
+/*}}}*/
+IoData::~IoData(){/*{{{*/
+        xDelete<char>(this->name);
+        xDelete<IssmDouble>(this->data);
+}
+/*}}}*/
+IoData::IoData(IssmDouble* matrix,int code_in,int layout_in,int M_in,int N_in,int enum_in){/*{{{*/
+        this->isindependent = false;
+        this->data_enum     = enum_in;
+        this->name          = NULL;
+        this->code          = code_in;
+        this->layout        = layout_in;
+        this->M             = M_in;
+        this->N             = N_in;
+        this->data          = matrix; /*do not copy*/
+        _assert_(code_in==5 ||  code_in==6 || code_in==7);
+}
+/*}}}*/
+
+void  IoModel::AddConstant(IoConstant* in_constant){/*{{{*/
+
+        _assert_(in_constant);
+
+        /*Go through dataset of constant and check whether it already exists */
+        vector<IoConstant*>::iterator iter;
+
+        for(iter=constants.begin();iter<constants.end();iter++){
+                if((*iter)->data_enum==in_constant->data_enum){
+                        delete in_constant;
+                        return;
+                }
+        }
+
+        this->constants.push_back(in_constant);
+}
+/*}}}*/
+void  IoModel::AddConstantIndependent(IoConstant* in_constant){/*{{{*/
+
+        _assert_(in_constant);
+
+        /*Set constant as independent*/
+        in_constant->isindependent = true;
+
+        /*Add to constnats*/
+        this->AddConstant(in_constant);
+}
+/*}}}*/
+void  IoModel::AddData(IoData* in_data){/*{{{*/
+
+        _assert_(in_data);
+
+        /*Go through dataset of data and check whether it already exists */
+        vector<IoData*>::iterator iter;
+
+        for(iter=data.begin();iter<data.end();iter++){
+                if((*iter)->data_enum==in_data->data_enum){
+                        delete in_data;
+                        return;
+                }
+        }
+
+        this->data.push_back(in_data);
+}
+/*}}}*/
+void  IoModel::AddDataIndependent(IoData* in_data){/*{{{*/
+
+        _assert_(in_data);
+
+        /*Set data as independent*/
+        in_data->isindependent = true;
+
+        /*Add to constnats*/
+        this->AddData(in_data);
+}
+/*}}}*/
+void  IoModel::FetchIndependentConstant(int* pXcount,IssmPDouble* X,int name){/*{{{*/
+
+        /*recover my_rank:*/
+        int my_rank=IssmComm::GetRank();
+
+        /*recover Xcount if X is not NULL:*/
+        int Xcount = 0;
+        if(X) Xcount=*pXcount;
+
+        #ifdef _HAVE_ADOLC_ //cannot come here unless you are running AD mode, from DeclaredIndependents:
+
+        /*output: */
+        IssmPDouble  pscalar;
+        IssmDouble   scalar; //same as pscalar, except it's an ADOLC independent variable
+        int          code;
+
+        /*Set file pointer to beginning of the data: */
+        fid=this->SetFilePointerToData(&code,NULL,name);
+        if(code!=3) _error_("expecting a IssmDouble for enum " << EnumToStringx(name));
+
+        /*We have to read a scalar from disk. First read the dimensions of the scalar, then the scalar: */
+        if(my_rank==0){
+                if(fread(&pscalar,sizeof(IssmPDouble),1,fid)!=1)_error_("could not read scalar ");
+
+                /*Now, before we even broadcast this to other nodes, declare the scalar  as an independent variable!. If we 
+                 *have been supplied an X vector, use it instead of what we just read: */
+                if(X){
+                        scalar<<=X[Xcount];
+                }
+                else{
+                        scalar<<=pscalar;
+                }
+        }
+
+        ISSM_MPI_Bcast(&scalar,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
+        this->AddConstantIndependent(new IoConstant(scalar,name));
+
+        /*increment offset into X vector, now that we have read 1 value:*/
+        Xcount++; *pXcount=Xcount;
+        #endif
+}
+/*}}}*/
+void  IoModel::FetchIndependentData(int* pXcount,IssmPDouble* X,int name){/*{{{*/
+
+        /*recover my_rank:*/
+        int my_rank=IssmComm::GetRank();
+
+        /*recover Xcount if X is not NULL:*/
+        int Xcount = 0;
+        if(X) Xcount=*pXcount;
+
+        #ifdef _HAVE_ADOLC_ //cannot come here unless you are running AD mode, from DeclaredIndependents:
+
+        /*Intermediaries*/
+        int M,N;
+        IssmPDouble* buffer=NULL; //a buffer to read the data from disk
+        IssmDouble* matrix=NULL; //our independent variable
+        int code,layout;
+
+        /*Set file pointer to beginning of the data: */
+        fid=this->SetFilePointerToData(&code,&layout,name);
+        if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(name));
+
+        /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */
+        /*numberofelements: */
+        if(my_rank==0){  
+                if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix ");
+        }
+        ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
+
+        if(my_rank==0){  
+                if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix ");
+        }
+        ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
+
+        /*Now allocate matrix: */
+        if(M*N){
+                buffer=xNew<IssmPDouble>(M*N);
+                matrix=xNew<IssmDouble>(M*N);
+
+                /*Read matrix on node 0, then broadcast: */
+                if(my_rank==0){  
+                        if(fread(buffer,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
+
+                        /*Now, before we even broadcast this to other nodes, declare the whole matrix as a independent variable!
+                          If we have been supplied an X vector, use it instead of what we just read: */
+                        if(X){
+                                for(int i=0;i<M*N;i++) matrix[i]<<=X[Xcount+i];  /*<<= ADOLC overloaded operator to declare independent*/
+                        }
+                        else{
+                                for(int i=0;i<M*N;i++) matrix[i]<<=buffer[i];
+                        }
+                }
+                ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
+
+                xDelete<IssmPDouble>(buffer);
+        }
+        else _error_("cannot declare the independent variable " << EnumToStringx(name) <<  "if it's empty!");
+
+        /*Add to data as independent*/
+        this->AddDataIndependent(new IoData(matrix,code,layout,M,N,name));
+
+        /*increment offset into X vector, now that we have read M*N values:*/
+        Xcount+=M*N; *pXcount=Xcount;
+        #endif
+}
+/*}}}*/
+void  IoModel::FindConstant(bool* pvalue,int constant_enum){/*{{{*/
+
+        /*Intermediary*/
+        vector<IoConstant*>::iterator iter;
+
+        for(iter=constants.begin();iter<constants.end();iter++){
+                IoConstant* ioconstant=*iter;
+
+                if(ioconstant->data_enum==constant_enum){
+                        ioconstant->constant->GetParameterValue(pvalue);
+                        return;
+                }
+        }
+
+        for(vector<IoConstant*>::iterator iter=constants.begin();iter<constants.end();iter++) (*iter)->constant->Echo();
+        _error_("Could not find constant \""<<EnumToStringx(constant_enum) <<"\"");
+}
+/*}}}*/
+void  IoModel::FindConstant(int* pvalue,int constant_enum){/*{{{*/
+
+        /*Intermediary*/
+        vector<IoConstant*>::iterator iter;
+
+        for(iter=constants.begin();iter<constants.end();iter++){
+                IoConstant* ioconstant=*iter;
+
+                if(ioconstant->data_enum==constant_enum){
+                        ioconstant->constant->GetParameterValue(pvalue);
+                        return;
+                }
+        }
+
+        _error_("Could not find constant \""<<EnumToStringx(constant_enum) <<"\"");
+}
+/*}}}*/
+void  IoModel::FindConstant(IssmDouble* pvalue,int constant_enum){/*{{{*/
+
+        /*Intermediary*/
+        vector<IoConstant*>::iterator iter;
+
+        for(iter=constants.begin();iter<constants.end();iter++){
+                IoConstant* ioconstant=*iter;
+
+                if(ioconstant->data_enum==constant_enum){
+                        ioconstant->constant->GetParameterValue(pvalue);
+                        return;
+                }
+        }
+
+        _error_("Could not find constant \""<<EnumToStringx(constant_enum) <<"\"");
+}
+/*}}}*/
+void  IoModel::FindConstant(char** pvalue,int constant_enum){/*{{{*/
+
+        /*Intermediary*/
+        vector<IoConstant*>::iterator iter;
+
+        for(iter=constants.begin();iter<constants.end();iter++){
+                IoConstant* ioconstant=*iter;
+
+                if(ioconstant->data_enum==constant_enum){
+                        ioconstant->constant->GetParameterValue(pvalue);
+                        return;
+                }
+        }
+
+        _error_("Could not find constant \""<<EnumToStringx(constant_enum) <<"\"");
+}
+/*}}}*/
+int   IoModel::NumIndependents(void){/*{{{*/
+
+        /*Initialize output*/
+        int num_independents = 0;
+
+        /*Process constants*/
+        for(vector<IoConstant*>::iterator iter=constants.begin();iter<constants.end();iter++){
+                if((*iter)->isindependent){
+                        num_independents+= 1;
+                }
+        }
+
+        /*Process data*/
+        for(vector<IoData*>::iterator iter=data.begin();iter<data.end();iter++){
+                if((*iter)->isindependent){
+                        num_independents+= (*iter)->M*(*iter)->N;
+                }
+        }
+
+        /*return*/
+        return num_independents;
+}
+/*}}}*/
+void  IoModel::FillIndependents(IssmDouble* xp){/*{{{*/
+
+        _assert_(xp);
+
+        /*Initialize local num ind*/
+        int local_num_ind = 0;
+
+        /*Process constants*/
+        for(vector<IoConstant*>::iterator iter=constants.begin();iter<constants.end();iter++){
+                if((*iter)->isindependent){
+                        (*iter)->constant->GetParameterValue(&xp[local_num_ind]);
+                        local_num_ind += 1;
+                }
+        }
+
+        /*Process data*/
+        for(vector<IoData*>::iterator iter=data.begin();iter<data.end();iter++){
+                if((*iter)->isindependent){
+                        for(int i=0;i<(*iter)->M*(*iter)->N;i++){
+                                xp[local_num_ind+i] = (*iter)->data[i];
+                        }
+                        local_num_ind += (*iter)->M*(*iter)->N;
+                }
+        }
+
+        _assert_(local_num_ind == this->NumIndependents());
+}
+/*}}}*/
+/*NEW NEW NEW NEW NEW NEW NEW NEW NEW NEW*/
+
 IoModel::IoModel(){/*{{{*/
+
         this->fid=NULL;
         this->solution_enum=-1;
-        this->data=NULL;
-        this->independents=NULL;
-        this->independent_objects=NULL;
-        this->constants=NULL;
 
         this->my_elements=NULL;
@@ -65 +426 @@
         this->solution_enum = solution_enum_in;
 
-        /*Initialize data: */
-        this->data=xNew<IssmDouble*>(MaximumNumberOfDefinitionsEnum);
-        for(int i=0;i<MaximumNumberOfDefinitionsEnum;i++) this->data[i]=NULL;
-
         /*If we are running in AD mode, we need to start the trace and declare our independent variables now, 
          *and prevent them from being erased during successive calls to iomodel->FetchConstants, iomodel->FetchData and 
@@ -76 +433 @@
 
         /*Initialize and read constants:*/
-        this->constants=new Parameters();
-        this->FetchConstants(); /*this routine goes through the input file, and fetches bools, ints, IssmDoubles and strings only, nothing memory intensive*/
-
-        /*Now some very specific piece of logic. We want to know whether we are going to carry out an autodiff run. There are 
-         *several cases: either trace is true (which bypasses the isautodiff choice), or trace is false but isautodiff is on, in which 
-         *case two cases are possible: are we running a control method or not? To make things simpler, we are going to pin down everyting 
-         *on isautodiff. By the way, we could not do this before we had the constants dataset initialized! {{{*/
-        this->constants->FindParam(&autodiff,AutodiffIsautodiffEnum);
-        this->constants->FindParam(&iscontrol,InversionIscontrolEnum);
-
-        if(trace)autodiff=true;
+        this->FetchConstants(); /*this routine goes through the input file, and fetches bool, int, IssmDouble and string only, nothing memory intensive*/
+
+        /*Is this an autodiff run?*/
+        this->FindConstant(&autodiff,AutodiffIsautodiffEnum);
+        this->FindConstant(&iscontrol,InversionIscontrolEnum);
+        if(trace){
+                autodiff=true;
+        }
         else{
-                if(autodiff && !iscontrol) autodiff=true;
-                else autodiff=false;
-        }
-        this->constants->SetParam(autodiff,AutodiffIsautodiffEnum);
-        /*}}}*/
+                if(autodiff && !iscontrol)
+                 autodiff=true;
+                else
+                 autodiff=false;
+        }
+        this->AddConstant(new IoConstant(autodiff,AutodiffIsautodiffEnum));
 
         /*Initialize permanent data: */
@@ -119 +474 @@
 IoModel::~IoModel(){/*{{{*/
 
-        /*Delete independents*/
-        if(this->independents){
-                for(int i=0;i<MaximumNumberOfDefinitionsEnum;i++){
-                        if(this->independents[i]){
-                                IssmDouble* array=this->data[i];
-                                xDelete<IssmDouble>(array);
-                        }
-                }
-        }
-
-        /*checks in debugging mode*/
-        #if defined(_ISSM_DEBUG_) && !defined(_HAVE_ADOLC_)
-        if(this->data){
-                for(int i=0;i<MaximumNumberOfDefinitionsEnum;i++){
-                        if(this->data[i]){
-                                _printf0_("WARNING: previous pointer of " << EnumToStringx(i) << " has not been freed (DeleteData has not been called)\n");
-                        }
-                }
-        }
-        #endif
-
-        if(this->constants) delete this->constants;
-
-        xDelete<IssmDouble*>(this->data);
-        xDelete<bool>(this->independents);
-        if(this->independent_objects)delete this->independent_objects;
+        /*Delete constants*/
+        vector<IoConstant*>::iterator iter1;
+        for(iter1=constants.begin();iter1<constants.end();iter1++){
+                delete *iter1;
+        }
+        this->constants.clear();
+
+        /*Delete data*/
+        vector<IoData*>::iterator iter2;
+        for(iter2=data.begin();iter2<data.end();iter2++){
+                #if defined(_ISSM_DEBUG_)
+                if(!(*iter2)->isindependent){
+                        _printf0_("WARNING: IoData " << EnumToStringx((*iter2)->data_enum) << " has not been freed (DeleteData has not been called)\n");
+                }
+                #endif
+                delete *iter2;
+        }
+        this->data.clear();
 
         xDelete<bool>(this->my_elements);
@@ -244 +591 @@
 }
 /*}}}*/
-void  IoModel::FindConstant(bool* poutput,int constant_enum){/*{{{*/
-
-        _assert_(constant_enum>=0);
-        _assert_(this->constants);
-
-        this->constants->FindParam(poutput,constant_enum);
-}
-/*}}}*/
-void  IoModel::FindConstant(int* poutput,int constant_enum){/*{{{*/
-
-        _assert_(constant_enum>=0);
-        _assert_(this->constants);
-
-        this->constants->FindParam(poutput,constant_enum);
-}
-/*}}}*/
-void  IoModel::FindConstant(IssmDouble* poutput,int constant_enum){/*{{{*/
-
-        _assert_(constant_enum>=0);
-        _assert_(this->constants);
-
-        this->constants->FindParam(poutput,constant_enum);
-}
-/*}}}*/
-void  IoModel::FindConstant(char** poutput,int constant_enum){/*{{{*/
-
-        _assert_(constant_enum>=0);
-        _assert_(this->constants);
-
-        this->constants->FindParam(poutput,constant_enum);
-}
-/*}}}*/
 Param* IoModel::CopyConstantObject(int constant_enum){/*{{{*/
 
-        _assert_(this->constants);
-
-        /*Find constant*/
-        //FIXME _error_("Needs to be FIXED");
-        Param* param=this->constants->FindParamObject(constant_enum);
-        if(!param) _error_("Constant " << EnumToStringx(constant_enum) << " not found in iomodel");
-
-        return xDynamicCast<Param*>(param->copy());
+        /*Intermediary*/
+        vector<IoConstant*>::iterator iter;
+
+        for(iter=constants.begin();iter<constants.end();iter++){
+                IoConstant* ioconstant=*iter;
+
+                if(ioconstant->data_enum==constant_enum){
+                        return ioconstant->constant->copy();
+                }
+        }
+
+        _error_("Constant " << EnumToStringx(constant_enum) << " not found in iomodel");
+        return NULL;
 }
 /*}}}*/
 IssmDouble* IoModel::Data(int data_enum){/*{{{*/
 
-        _assert_(data_enum<MaximumNumberOfDefinitionsEnum);
-        _assert_(data_enum>=0);
-
-        return this->data[data_enum];
+        /*Intermediary*/
+        vector<IoData*>::iterator iter;
+
+        for(iter=data.begin();iter<data.end();iter++){
+                IoData* iodata=*iter;
+                if(iodata->data_enum==data_enum) return iodata->data;
+        }
+
+        return NULL;
 }
 /*}}}*/
@@ -338 +663 @@
 void  IoModel::DeclareIndependents(bool trace,IssmPDouble* X){/*{{{*/
 
-        int  i;
-        bool autodiff = false;
-        bool iscontrol = false;
+        bool autodiff,iscontrol;
         int  num_independent_objects;
         int  Xcount=0;
@@ -347 +670 @@
         int *types = NULL;
 
-        int  dummy;
-
         /*Initialize array detecting whether data[i] is an independent AD mode variable: */
-        this->independents=xNew<bool>(MaximumNumberOfDefinitionsEnum);
-        for(i=0;i<MaximumNumberOfDefinitionsEnum;i++) this->independents[i]=false;
-
         this->FetchData(&autodiff,AutodiffIsautodiffEnum);
         this->FetchData(&iscontrol,InversionIscontrolEnum);
@@ -359 +677 @@
 
                 #ifdef _HAVE_ADOLC_
-                /*build dataset made of independent objects:*/
-                this->independent_objects=new DataSet();
                 this->FetchData(&num_independent_objects,AutodiffNumIndependentObjectsEnum);
                 if(num_independent_objects){
-                        this->FetchData(&names,&dummy,&dummy,AutodiffIndependentObjectNamesEnum);
-                        this->FetchData(&types,&dummy,&dummy,AutodiffIndependentObjectTypesEnum);
+                        this->FetchData(&names,NULL,NULL,AutodiffIndependentObjectNamesEnum);
+                        this->FetchData(&types,NULL,NULL,AutodiffIndependentObjectTypesEnum);
 
                         /*create independent objects, and at the same time, fetch the corresponding independent variables, 
                          *and declare them as such in ADOLC: */
-                        for(i=0;i<num_independent_objects;i++){
-
-                                IndependentObject* independent_object=NULL;
-                                independent_object=new IndependentObject(names[i],types[i]);
-
-                                /*add to independent_objects dataset:*/
-                                this->independent_objects->AddObject(independent_object);
-
-                                /*now go fetch the independent variable: */
-                                independent_object->FetchIndependent(this,&Xcount,X); //supply the pointer to iomodel.
+                        for(int i=0;i<num_independent_objects;i++){
+
+                                if(types[i]==0){
+                                        /*Scalar*/
+                                        this->FetchIndependentConstant(&Xcount,X,names[i]);
+                                }
+                                else if(types[i]==1){
+                                        /* vector:*/
+                                        this->FetchIndependentData(&Xcount,X,names[i]);
+                                }
+                                else{
+                                        _error_("Independent cannot be of size " << types[i]);
+                                }
                         }
                         xDelete<int>(names);
@@ -387 +706 @@
                 #endif
         }
-        else this->independent_objects=NULL;
-
 }
 /*}}}*/
 void  IoModel::DeleteData(int num,...){/*{{{*/
 
+        /*Intermediaries*/
         va_list ap;
         int     dataenum;
-
-        /*Go through the entire list of enums and delete the corresponding data from the iomodel-data dataset: */
+        vector<IoData*>::iterator iter;
+
+        /*Go through the entire list of data and delete the corresponding data from the iomodel-data dataset: */
         va_start(ap,num);
         for(int i=0;i<num;i++){
                 dataenum=va_arg(ap,int);
-                _assert_(dataenum<MaximumNumberOfDefinitionsEnum);
-
-                /*do not erase independent variables for the AD mode computations!: */
-                if (!this->independents[dataenum]) xDelete<IssmDouble>(this->data[dataenum]);
+
+                for(iter=data.begin();iter<data.end();iter++){
+                        IoData* iodata=*iter;
+                        if(iodata->data_enum==dataenum && !iodata->isindependent){
+                                delete *iter;
+                                this->data.erase(iter);
+                                break;
+                        }
+                }
         }
         va_end(ap);
 } /*}}}*/
-void  IoModel::DeleteData(IssmDouble* vector, int dataenum){/*{{{*/
-
-        /*do not erase independent variables for the AD mode computations!: */
-        if(vector)if (!this->independents[dataenum]) xDelete<IssmDouble>(vector);
-
+void  IoModel::DeleteData(IssmDouble* vector_in, int dataenum){/*{{{*/
+
+        vector<IoData*>::iterator iter;
+
+        /*do not do anything if pointer is NULL*/
+        if(!vector_in) return;
+
+        /*do not delete if this is an independent variable*/
+        for(iter=data.begin();iter<data.end();iter++){
+                IoData* iodata=*iter;
+                if(iodata->data_enum==dataenum && iodata->isindependent){
+                        return;
+                }
+        }
+
+        /*Go ahead and delete*/
+        xDelete<IssmDouble>(vector_in);
 } /*}}}*/
 void  IoModel::DeleteData(char*** pstringarray, int numstrings, int dataenum){/*{{{*/
@@ -429 +765 @@
 void  IoModel::FetchConstants(void){/*{{{*/
 
-        int my_rank;
-
         /*record descriptions; */
         int record_enum;
@@ -437 +771 @@
 
         /*records: */
-        int  booleanint=0;
-        int  integer=0;
-        IssmPDouble pscalar=0;
-        IssmDouble scalar=0;
-        char* string=NULL;
-        int   string_size;
+        int          booleanint  = 0;
+        int          integer     = 0;
+        IssmPDouble  pscalar     = 0;
+        IssmDouble   scalar      = 0;
+        char        *string      = NULL;
+        int          string_size;
 
         /*recover my_rank:*/
-        my_rank=IssmComm::GetRank();
+        int my_rank=IssmComm::GetRank();
 
         /*Check that some fields have been allocated*/
         _assert_(this->fid || my_rank);
-        _assert_(this->constants);
 
         /*Go find in the binary file, the position of the data we want to fetch: */
@@ -487 +820 @@
                                                 /*create BoolParam: */
                                                 if(record_enum!=MaximumNumberOfDefinitionsEnum && record_enum!=MaximumNumberOfDefinitionsEnum+1)
-                                                 this->constants->AddObject(new BoolParam(record_enum,(bool)booleanint)); //cast to boolean
+                                                 this->AddConstant(new IoConstant((bool)booleanint,record_enum)); //cast to boolean
 
                                                 break;
@@ -496 +829 @@
 
                                                 /*create IntParam: */
-                                                this->constants->AddObject(new IntParam(record_enum,integer));
+                                                this->AddConstant(new IoConstant(integer,record_enum));
 
                                                 break;
                                         case 3:
-                                                /*Read the scalar and broadcast it to other cpus. However, if this record has already 
-                                                 * been read in DeclareIndependents, we should not re-read it, but grab it from the iomodel->data 
-                                                 * slots: */
-                                                if(this->independents[record_enum]){
-                                                        scalar=*(this->data[record_enum]);
-                                                }
-                                                else{
-                                                        if(fread(&pscalar,sizeof(IssmPDouble),1,this->fid)!=1) _error_("could not read scalar ");
-                                                        ISSM_MPI_Bcast(&pscalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
-                                                        scalar=pscalar;
-                                                }
-
-                                                /*create DoubleParam: */
-                                                this->constants->AddObject(new DoubleParam(record_enum,scalar));
-
+                                                  {
+                                                        /*IssmDouble, check whether it is already there (from "declare independents")*/
+                                                        bool exists = false;
+                                                        vector<IoConstant*>::iterator iter;
+                                                        for(iter=constants.begin();iter<constants.end();iter++){
+                                                                IoConstant* ioconstant=*iter;
+                                                                if(ioconstant->data_enum==record_enum){
+                                                                        exists = true;
+                                                                        break;
+                                                                }
+                                                        }
+                                                        if(!exists){
+                                                                if(fread(&pscalar,sizeof(IssmPDouble),1,this->fid)!=1) _error_("could not read scalar ");
+                                                                ISSM_MPI_Bcast(&pscalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
+                                                                scalar=pscalar;
+
+                                                                /*create DoubleParam: */
+                                                                this->AddConstant(new IoConstant(scalar,record_enum));
+                                                        }
+                                                  }
                                                 break;
                                         case 4: 
@@ -535 +873 @@
 
                                                 /*Add string to parameters: */
-                                                this->constants->AddObject(new StringParam(record_enum,string));
+                                                this->AddConstant(new IoConstant(string,record_enum));
 
                                                 /*Free string*/
                                                 xDelete<char>(string);
-
                                                 break;
                                         case 5: 
-                                                        /*We are not interested in this record, too memory intensive. Skip it: */
-                                                        /*skip: */
-                                                        fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer
-                                                        fseek(fid,record_length,SEEK_CUR);
-                                                        break;
                                         case 6: 
-                                                        /*We are not interested in this record, too memory intensive. Skip it: */
-                                                        /*skip: */
-                                                        fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer
-                                                        fseek(fid,record_length,SEEK_CUR);
-                                                        break;
                                         case 7: 
-                                                        /*We are not interested in this record, too memory intensive. Skip it: */
-                                                        /*skip: */
-                                                        fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer
-                                                        fseek(fid,record_length,SEEK_CUR);
-                                                        break;
-
                                         case 8: 
-                                                        /*We are not interested in this record, too memory intensive. Skip it: */
-                                                        /*skip: */
-                                                        fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer
-                                                        fseek(fid,record_length,SEEK_CUR);
-                                                        break;
-
                                         case 9: 
                                                         /*We are not interested in this record, too memory intensive. Skip it: */
@@ -573 +888 @@
                                                         fseek(fid,record_length,SEEK_CUR);
                                                         break;
-
                                         default: 
                                                 _error_("unknown record type:" << record_code); 
-                                                break;;
+                                                break;
                                 }
                         }
@@ -597 +911 @@
                                         /*create BoolParam: */
                                         if(record_enum!=MaximumNumberOfDefinitionsEnum && record_enum!=MaximumNumberOfDefinitionsEnum+1)
-                                         this->constants->AddObject(new BoolParam(record_enum,(bool)booleanint)); //cast to a boolean
+                                         this->AddConstant(new IoConstant((bool)booleanint,record_enum)); //cast to a boolean
                                         break;
 
@@ -605 +919 @@
 
                                         /*create IntParam: */
-                                        this->constants->AddObject(new IntParam(record_enum,integer));
+                                        this->AddConstant(new IoConstant(integer,record_enum));
 
                                         break;
                                 case 3:
                                         /*scalar. get it from cpu 0 */
-                                        ISSM_MPI_Bcast(&pscalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
-                                        scalar=pscalar;
-                                        /*create DoubleParam: */
-                                        this->constants->AddObject(new DoubleParam(record_enum,scalar));
-
+                                          {
+                                                /*IssmDouble, check whether it is already there (from "declare independents")*/
+                                                bool exists = false;
+                                                vector<IoConstant*>::iterator iter;
+                                                for(iter=constants.begin();iter<constants.end();iter++){
+                                                        IoConstant* ioconstant=*iter;
+                                                        if(ioconstant->data_enum==record_enum){
+                                                                exists = true;
+                                                                break;
+                                                        }
+                                                }
+                                                if(!exists){
+                                                        ISSM_MPI_Bcast(&pscalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+                                                        scalar=pscalar;
+                                                        /*create DoubleParam: */
+                                                        this->AddConstant(new IoConstant(scalar,record_enum));
+                                                }
+                                          }
                                         break;
                                 case 4: 
@@ -630 +957 @@
                                         }
                                         /*Add string to parameters: */
-                                        this->constants->AddObject(new StringParam(record_enum,string));
+                                        this->AddConstant(new IoConstant(string,record_enum));
 
                                         /*Free string*/
@@ -644 +971 @@
                                 default: 
                                         _error_("unknown record type:" << record_code); 
-                                        break;;
+                                        break;
                                 }
 
@@ -653 +980 @@
 /*}}}*/
 void  IoModel::FetchData(bool* pboolean,int data_enum){/*{{{*/
-
-        int my_rank;
 
         /*output: */
@@ -661 +986 @@
 
         /*recover my_rank:*/
-        my_rank=IssmComm::GetRank();
+        int my_rank=IssmComm::GetRank();
 
         /*Set file pointer to beginning of the data: */
@@ -682 +1007 @@
 void  IoModel::FetchData(int* pinteger,int data_enum){/*{{{*/
 
-        int my_rank;
-
         /*output: */
         int   integer;
@@ -689 +1012 @@
 
         /*recover my_rank:*/
-        my_rank=IssmComm::GetRank();
+        int my_rank=IssmComm::GetRank();
 
         /*Set file pointer to beginning of the data: */
@@ -709 +1032 @@
 void  IoModel::FetchData(IssmDouble* pscalar,int data_enum){/*{{{*/
 
-        int my_rank;
-
         /*output: */
         IssmPDouble   scalar;
@@ -716 +1037 @@
 
         /*recover my_rank:*/
-        my_rank=IssmComm::GetRank();
+        int my_rank=IssmComm::GetRank();
 
         /*Set file pointer to beginning of the data: */
@@ -735 +1056 @@
 /*}}}*/
 void  IoModel::FetchData(char** pstring,int data_enum){/*{{{*/
-
-        int my_rank;
 
         /*output: */
@@ -744 +1063 @@
 
         /*recover my_rank:*/
-        my_rank=IssmComm::GetRank();
+        int my_rank=IssmComm::GetRank();
 
         /*Set file pointer to beginning of the data: */
@@ -781 +1100 @@
 /*}}}*/
 void  IoModel::FetchData(int** pmatrix,int* pM,int* pN,int data_enum){/*{{{*/
-
-        int my_rank;
         int i,j;
 
@@ -792 +1109 @@
 
         /*recover my_rank:*/
-        my_rank=IssmComm::GetRank();
+        int my_rank=IssmComm::GetRank();
 
         /*Set file pointer to beginning of the data: */
@@ -850 +1167 @@
 void  IoModel::FetchData(IssmDouble** pmatrix,int* pM,int* pN,int data_enum){/*{{{*/
 
-        int my_rank;
-
+        /*First, look if has already been loaded (might be an independent variable)*/
+        vector<IoData*>::iterator iter;
+        for(iter=data.begin();iter<data.end();iter++){
+                IoData* iodata=*iter;
+                if(iodata->data_enum==data_enum){
+                        *pmatrix=iodata->data;
+                        if(pM) *pM=iodata->M;
+                        if(pN) *pN=iodata->N;
+                        return;
+                }
+        }
+         
         /*output: */
-        int M,N;
-        IssmPDouble* matrix=NULL;
-        int code=0;
+        int          M,N;
+        IssmPDouble *matrix = NULL;
+        int          code   = 0;
 
         /*recover my_rank:*/
-        my_rank=IssmComm::GetRank();
+        int my_rank=IssmComm::GetRank();
 
         /*Set file pointer to beginning of the data: */
@@ -868 +1195 @@
         /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */
         /*numberofelements: */
-        if(my_rank==0){  
+        if(my_rank==0){
                 if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix ");
         }
@@ -888 +1215 @@
                 ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
 
-                _assert_(this->independents);
-                if(this->independents[data_enum]){
-                        /*this data has already been checked out! So cancel all that we've done here, and return 
-                         * the data[data_enum] directly: */
-                        *pmatrix=this->data[data_enum];
-                }
-                else{
-                        *pmatrix=xNew<IssmDouble>(M*N);
-                        for(int i=0;i<M*N;++i) (*pmatrix)[i]=matrix[i];
-                }
+                *pmatrix=xNew<IssmDouble>(M*N);
+                for(int i=0;i<M*N;++i) (*pmatrix)[i]=matrix[i];
                 xDelete<IssmPDouble>(matrix);
         }
@@ -903 +1222 @@
           *pmatrix=NULL;
         /*Assign output pointers: */
-        if (pM)*pM=M;
-        if (pN)*pN=N;
+        if(pM) *pM=M;
+        if(pN) *pN=N;
 }
 /*}}}*/
 void  IoModel::FetchData(char*** pstrings,int* pnumstrings,int data_enum){/*{{{*/
-
-        int my_rank;
 
         int i;
@@ -923 +1240 @@
 
         /*recover my_rank:*/
-        my_rank=IssmComm::GetRank();
+        int my_rank=IssmComm::GetRank();
 
         /*Set file pointer to beginning of the data: */
@@ -975 +1292 @@
 
         int i;
-
-        int my_rank;
-
         /*output: */
         IssmDouble** matrices=NULL;
@@ -990 +1304 @@
 
         /*recover my_rank:*/
-        my_rank=IssmComm::GetRank();
+        int my_rank=IssmComm::GetRank();
 
         /*Set file pointer to beginning of the data: */
@@ -1105 +1419 @@
 
         va_list     ap;
-        int         dataenum;
+        int         dataenum,code,layout;
         IssmDouble *matrix   = NULL;
         int         M,N;
-        int         i;
+        bool        exists;
+        vector<IoData*>::iterator iter;
 
         /*Go through the entire list of enums and fetch the corresponding data. Add it to the iomodel->data dataset. Everything
@@ -1114 +1429 @@
 
         va_start(ap,num);
-        for(i=0; i<num; i++){
+        for(int i=0; i<num; i++){
 
                 dataenum=va_arg(ap, int);
-                _assert_(dataenum<MaximumNumberOfDefinitionsEnum); 
-                _assert_(dataenum>=0);
-
-                if (this->independents[dataenum]){
+                exists = false;
+
+                for(iter=data.begin();iter<data.end();iter++){
+                        IoData* iodata=*iter;
+                        if(iodata->data_enum==dataenum){
+                                /*Already there, no need to fetch it*/
+                                _assert_(iodata->isindependent);
+                                exists = true;
+                                break;
+                        }
+                }
+
+                if(exists){
                         /*this data has already been checked out! Continue: */
                         continue;
                 }
-
-                /*Some checks in debugging mode*/
-                /*{{{*/
-                #ifdef _ISSM_DEBUG_
-                _assert_(dataenum<MaximumNumberOfDefinitionsEnum);
-                if(this->data[dataenum]){
-                        _error_("Info: trying to fetch " << EnumToStringx(dataenum) << " but previous pointer has not been freed (DeleteData has not been called)");
-                }
-                #endif
-                /*}}}*/
-
-                /*Add to this->data: */
-                this->FetchData(&matrix,&M,&N,dataenum);
-                this->data[dataenum]=matrix;
+                else{
+                        /*Add to this->data: */
+                        this->SetFilePointerToData(&code,&layout,dataenum);
+                        this->FetchData(&matrix,&M,&N,dataenum);
+                        this->AddData(new IoData(matrix,code,layout,M,N,dataenum));
+                }
         }
         va_end(ap);
-
 }
 /*}}}*/
@@ -1373 +1688 @@
                                 ISSM_MPI_Bcast(pmatrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
 
-                                _assert_(this->independents);
-                                if(this->independents[data_enum]){
-                                        /*this data has already been checked out! So cancel all that we've done here, and return 
-                                         * the data[data_enum] directly: */
-                                        matrix=this->data[data_enum];
-                                }
-                                else{
+                                //if(this->independents[data_enum]){ FIXME
+                                //      /*this data has already been checked out! So cancel all that we've done here, and return 
+                                //       * the data[data_enum] directly: */
+                                //      matrix=this->data[data_enum];
+                                //}
+                                //else{
                                         matrix=xNew<IssmDouble>(M*N);
                                         for (int i=0;i<M*N;++i) matrix[i]=pmatrix[i];
-                                }
+                                //}
                                 xDelete<IssmPDouble>(pmatrix);
                         }
@@ -1479 +1793 @@
                                 ISSM_MPI_Bcast(pmatrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
 
-                                _assert_(this->independents);
-                                if(this->independents[data_enum]){
-                                        /*this data has already been checked out! So cancel all that we've done here, and return 
-                                         * the data[data_enum] directly: */
-                                        matrix=this->data[data_enum];
-                                        for (int i=0;i<M*N;++i) integer_matrix[i]=reCast<int>(matrix[i]);
-                                }
-                                else{
+                                //if(this->independents[data_enum]){ FIXME
+                                //      /*this data has already been checked out! So cancel all that we've done here, and return 
+                                //       * the data[data_enum] directly: */
+                                //      matrix=this->data[data_enum];
+                                //      for (int i=0;i<M*N;++i) integer_matrix[i]=reCast<int>(matrix[i]);
+                                //}
+                                //else{
                                         for (int i=0;i<M*N;++i) integer_matrix[i]=pmatrix[i];
-                                }
+                                //}
                                 xDelete<IssmPDouble>(pmatrix);
                         }
@@ -1525 +1838 @@
 void  IoModel::FetchDataToInput(Elements* elements,int vector_enum,IssmDouble default_value){/*{{{*/
 
+        /*First, look whether it is not already loaded in this->data*/
+        vector<IoData*>::iterator iter;
+        for(iter=data.begin();iter<data.end();iter++){
+                IoData* iodata=*iter;
+                if(iodata->data_enum==vector_enum){
+                        _assert_(iodata->code==7);
+                        for(int i=0;i<elements->Size();i++){
+                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                                element->InputCreate(iodata->data,this,iodata->M,iodata->N,iodata->layout,vector_enum,iodata->code);//we need i to index into elements.
+                        }
+                        return;
+                }
+        }
+
         /*intermediary: */
         int         code,vector_layout;
@@ -1544 +1871 @@
         }
 
-        /*Free ressources. Pay attention to not freeing an AD mode independent variable though!:*/
-        if(!this->independents[vector_enum]) xDelete<IssmDouble>(doublearray);
+        /*Free ressources*/
+        xDelete<IssmDouble>(doublearray);
 }
 /*}}}*/
 void  IoModel::FetchDataToInput(Elements* elements,int vector_enum){/*{{{*/
+
+        /*First, look whether it is not already loaded in this->data*/
+        vector<IoData*>::iterator iter;
+        for(iter=data.begin();iter<data.end();iter++){
+                IoData* iodata=*iter;
+                if(iodata->data_enum==vector_enum){
+                        for(int i=0;i<elements->Size();i++){
+                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                                element->InputCreate(iodata->data,this,iodata->M,iodata->N,iodata->layout,vector_enum,iodata->code);//we need i to index into elements.
+                        }
+                        return;
+                }
+        }
 
         /*intermediary: */
@@ -1615 +1955 @@
                         break;
         }
-        /*Free ressources. Pay attention to not freeing an AD mode independent variable though!:*/
-        if (!this->independents[vector_enum]) xDelete<IssmDouble>(doublearray);
+        /*Free ressources*/
+        xDelete<IssmDouble>(doublearray);
         xDelete<char>(string);
 }

issm/trunk-jpl/src/c/classes/IoModel.h

@@ -9 +9 @@
 
 #include "../shared/Enum/Enum.h"
+#include <vector>
 
 class Parameters;
@@ -15 +16 @@
 class Option;
 
+class IoConstant { /*holds single IssmDouble, int, bool and char from input*/
+        public:
+                bool   isindependent;
+                int    data_enum;
+                char*  name;
+                Param* constant; 
+
+                ~IoConstant();
+                IoConstant();
+                IoConstant(bool value,int enum_in);
+                IoConstant(int value,int enum_in);
+                IoConstant(IssmDouble value,int enum_in);
+                IoConstant(char* value,int enum_in);
+};
+
+class IoData { /*holds temporary data (array), memory intensive*/
+        public:
+                bool        isindependent;
+                int         data_enum;
+                char*       name;
+                int         M,N;
+                int         code;
+                int         layout;
+                IssmDouble* data;
+
+                ~IoData();
+                IoData();
+                IoData(IssmDouble* matrix,int code,int layout_in,int M,int N,int enum_in);
+};
+
 class IoModel {
 
         private: 
-                Parameters *constants;   //this dataset holds all IssmDouble, int, bool and char from input
+                std::vector<IoConstant*> constants; //this dataset holds all IssmDouble, int, bool and char from input
+                std::vector<IoData*>     data;      //this dataset holds temporary data, memory intensive
+
+                /*for AD mode: to keep track of our independent variables we fetch:*/
+                //bool    *independents;
+                //DataSet *independent_objects;
 
         public:
-                IssmDouble **data;   //this dataset holds temporary data, memory intensive.
-
                 /*pointer to input file*/
                 FILE *fid;
@@ -55 +89 @@
                 int constraintcounter;   //keep track of how many constraints are being created in each analysis
 
-                /*for AD mode: to keep track of our independent variables we fetch:*/
-                bool    *independents;
-                DataSet *independent_objects;
-
                 /*Methods*/
                 ~IoModel();
@@ -64 +94 @@
                 IoModel(FILE* iomodel_handle,int solution_enum_in,bool trace,IssmPDouble* X);
 
+                /*NEW*/
+                void        AddConstant(IoConstant* constant_in);
+                void        AddConstantIndependent(IoConstant* constant_in);
+                void        AddData(IoData* data_in);
+                void        AddDataIndependent(IoData* data_in);
+                void        FindConstant(bool* pvalue,int constant_enum);
+                void        FindConstant(int* pvalue,int constant_enum);
+                void        FindConstant(IssmDouble* pvalue,int constant_enum);
+                void        FindConstant(char **pvalue,int constant_enum);
+                void        FetchIndependentConstant(int* pXcount,IssmPDouble* X,int name);
+                void        FetchIndependentData(int* pXcount,IssmPDouble* X,int name);
+                int         NumIndependents();
+                void        FillIndependents(IssmDouble* xp);
+
                 /*Input/Output*/
                 void        CheckEnumSync(void);
-                void        FindConstant(bool *poutput,int constant_enum);
-                void        FindConstant(int *poutput,int constant_enum);
-                void        FindConstant(IssmDouble *poutput,int constant_enum);
-                void        FindConstant(char **poutput,int constant_enum);
                 Param      *CopyConstantObject(int constant_enum);
                 IssmDouble *Data(int dataenum);

issm/trunk-jpl/src/c/main/kriging.cpp

@@ -145 +145 @@
         iomodel->fid=fid;
         iomodel->CheckEnumSync();
-        iomodel->independents=xNew<bool>(MaximumNumberOfDefinitionsEnum); for(int i=0;i<MaximumNumberOfDefinitionsEnum;i++) iomodel->independents[i]=false;
         iomodel->FetchData(&x,&M,&N,0);        nobs=M*N;
         iomodel->FetchData(&y,&M,&N,1);        _assert_(M*N==nobs);

issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp

@@ -133 +133 @@
                 /*Independents have already been recovered in iomodel->DeclareIndependents. Just do some more processing. 
                  *In particular, figure out num_independents, and create the state vector xp, or size num_independents x 1 :*/
-                num_ind=0;
-                for(i=0;i<iomodel->independent_objects->Size();i++){
-                        IndependentObject* ind=(IndependentObject*)iomodel->independent_objects->GetObjectByOffset(i);
-                        num_ind+=ind->NumIndependents();
-                }
+                num_ind=iomodel->NumIndependents();
+                parameters->AddObject(new IntParam(AutodiffNumIndependentsEnum,num_ind));
+
                 if(num_ind){
                         xp=xNew<IssmDouble>(num_ind);
-                        xp_backup=xp;
-                        for(i=0;i<iomodel->independent_objects->Size();i++){
-                                IndependentObject* ind=(IndependentObject*)iomodel->independent_objects->GetObjectByOffset(i);
-                                ind->FillIndependents(iomodel->data,xp);
-                                local_num_ind=ind->NumIndependents(); xp=xp+local_num_ind;
-                        }
-                        xp=xp_backup; parameters->AddObject(new DoubleVecParam(AutodiffXpEnum,xp,num_ind));
+                        iomodel->FillIndependents(xp);
+                        parameters->AddObject(new DoubleVecParam(AutodiffXpEnum,xp,num_ind));
                         xDelete<IssmDouble>(xp);
                 }
-                parameters->AddObject(new IntParam(AutodiffNumIndependentsEnum,num_ind));
-
-                /*Don't forget to copy  iomodel->independent_objects to parameters: */
-                parameters->AddObject(new DataSetParam(AutodiffIndependentObjectsEnum,iomodel->independent_objects));
                 /*}}}*/
         }