Changeset 23231

Timestamp:

09/08/18 13:03:10 (7 years ago)

Author:

kruegern

Message:

NEW: added Dakota/QMU Chaco tests (234,235,413,414,417,444), minor bug fixes, left out 418 due to unresolved memory issues

Location:

issm/trunk-jpl

Files:

• src/m/modules/Chaco.py (modified) (1 diff)
• src/m/modules/Scotch.py (modified) (1 diff)
• src/m/partition/partitioner.py (modified) (2 diffs)
• src/m/qmu/dakota_in_write.py (modified) (1 diff)
• src/m/qmu/dakota_out_parse.py (modified) (1 diff)
• src/m/qmu/process_qmu_response_data.m (modified) (1 diff)
• src/wrappers/Chaco/Chaco.cpp (modified) (1 diff)
• src/wrappers/python/io/FetchPythonData.cpp (modified) (12 diffs)
• src/wrappers/python/io/pythonio.h (modified) (1 diff)
• test/NightlyRun/test234.py (added)
• test/NightlyRun/test235.py (added)
• test/NightlyRun/test413.py (added)
• test/NightlyRun/test414.py (added)
• test/NightlyRun/test417.py (added)
• test/NightlyRun/test444.py (added)

issm/trunk-jpl/src/m/modules/Chaco.py

@@ -5 +5 @@
 
    Usage:
-      [assgn] = Chaco(A,vwgts,ewgts,x,y,z,options,nparts,goal);
+      assgn = Chaco(A,vwgts,ewgts,x,y,z,options,nparts,goal);
 
    A:                   Input adjacency matrix

issm/trunk-jpl/src/m/modules/Scotch.py

@@ -2 +2 @@
 
 def Scotch(*varargin):
-'''SCOTCH - Scotch partitioner
+        '''SCOTCH - Scotch partitioner
 
    Usage:

issm/trunk-jpl/src/m/partition/partitioner.py

@@ -3 +3 @@
 import MatlabFuncs as m
 from adjacency import *
-#from Chaco import *
+from Chaco import *
 #from Scotch import *
-#from MeshPartition import *
+from MeshPartition import *
 from project3d import *
+from mesh2d import *
 
 def partitioner(md,*varargin):
@@ -60 +61 @@
 
         if m.strcmpi(package,'chaco'):
-                raise RuntimeError('Chaco is not currently supported for this function')
+                #raise RuntimeError('Chaco is not currently supported for this function')
 
                 #  default method (from chaco.m)
-                #method=np.array([1,1,0,0,1,1,50,0,.001,7654321]).reshape(-1,1)
-                #method[0]=3    #  global method (3=inertial (geometric))
-                #method[2]=0    #  vertex weights (0=off, 1=on)
+                method=np.array([1,1,0,0,1,1,50,0,.001,7654321])
+                method[0]=3    #  global method (3=inertial (geometric))
+                method[2]=0    #  vertex weights (0=off, 1=on)
 
                 #specify bisection
-                #method[5]=options.getfieldvalue('section')#  ndims (1=bisection, 2=quadrisection, 3=octasection)
+                method[5]=options.getfieldvalue('section')#  ndims (1=bisection, 2=quadrisection, 3=octasection)
 
                 #are we using weights? 
-                #if m.strcmpi(options.getfieldvalue('weighting'),'on'):
-                        #weights=np.floor(md.qmu.vertex_weight/min(md.qmu.vertex_weight))
-                        #method[2]=1
-                #else:
-                        #weights=[]
+                if m.strcmpi(options.getfieldvalue('weighting'),'on'):
+                        weights=np.floor(md.qmu.vertex_weight/min(md.qmu.vertex_weight))
+                        method[2]=1
+                else:
+                        weights=[]
         
+                method = method.reshape(-1,1)   # transpose to 1x10 instead of 10
+
                 #  partition into nparts
-                #if isinstance(md.mesh,mesh2d):
-                        #part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x, md.mesh.y,np.zeros((md.mesh.numberofvertices,)),method,npart,[]).T+1 #index partitions from 1 up. like metis.
-                #else:
-                        #part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x, md.mesh.y,md.mesh.z,method,npart,[]).T+1 #index partitions from 1 up. like metis.
-        
+                if isinstance(md.mesh,mesh2d):
+                        part=np.array(Chaco(md.qmu.adjacency,weights,np.array([]),md.mesh.x, md.mesh.y,np.zeros((md.mesh.numberofvertices,)),method,npart,np.array([]))).T+1 #index partitions from 1 up. like metis.
+                else:
+                        part=np.array(Chaco(md.qmu.adjacency,weights,np.array([]),md.mesh.x, md.mesh.y,md.mesh.z,method,npart,np.array([]))).T+1 #index partitions from 1 up. like metis.
         
         elif m.strcmpi(package,'scotch'):

issm/trunk-jpl/src/m/qmu/dakota_in_write.py

@@ -62 +62 @@
         filei2=fullfile(pathstr,name+ext)
 
-        print 'Opening Dakota input file \''+filei2 + '\''
+        print 'Opening Dakota input file \''+filei2 + '\'.'
         try:
                 with open(filei2,'w+') as fidi:

issm/trunk-jpl/src/m/qmu/dakota_out_parse.py

@@ -241 +241 @@
         dmax95 =prctile_issm(data,95,0)
 
+        # Note: the following line may cause the following warning
+        #       (should not crash or invalidate results) when one of
+        #       the inputs does not change with respect to the
+        #       other/s causing an internal divide-by-zero error
+
+        #/usr/local/lib/python2.7/dist-packages/numpy/lib/function_base.py:3163:
+        #       RuntimeWarning: invalid value encountered in true_divide
+        #       c /= stddev[:,None]
+
+        #       (and/or the same but with "c /= stddev[None, :]")
+
         dcorrel=np.corrcoef(data.T)
 

issm/trunk-jpl/src/m/qmu/process_qmu_response_data.m

@@ -43 +43 @@
         %ok, process the domains named in qmu_mass_flux_profiles,  to build a list of segments (MatArray)
         md.qmu.mass_flux_segments=cell(num_mass_flux,1);
-        md.qmu.mass_flux_segments
         for i=1:num_mass_flux,
                 md.qmu.mass_flux_segments{i}=MeshProfileIntersection(md.mesh.elements,md.mesh.x,md.mesh.y,[md.qmu.mass_flux_profile_directory '/' md.qmu.mass_flux_profiles{i}]);
         end
-        md.qmu.mass_flux_segments
 end

issm/trunk-jpl/src/wrappers/Chaco/Chaco.cpp

@@ -47 +47 @@
         /*Fetch Data*/
         FetchChacoData(&nvtxs,&adjacency,&start,&ewgts,A_IN,EWGTS_IN);
-        FetchData(&vwgts,&nterms,VWGTS_IN); 
-        FetchData(&x,&nterms,X_IN); 
-        FetchData(&y,&nterms,Y_IN); 
-        FetchData(&z,&nterms,Z_IN); 
-        FetchData(&in_options,&nterms,OPTNS_IN); 
+        FetchData(&vwgts,&nterms,VWGTS_IN);
+        FetchData(&x,&nterms,X_IN);
+        FetchData(&y,&nterms,Y_IN);
+        FetchData(&z,&nterms,Z_IN);
+        FetchData(&in_options,&nterms,OPTNS_IN);
         for (i=0;i<(nterms<10?nterms:10);i++) options[i]=in_options[i]; //copy in_options into default options
-        FetchData(&npart,NPARTS_IN); 
+        FetchData(&npart,NPARTS_IN);
         //int * nparts=xNew<int>(1); nparts[0]=npart; //weird Chacox interface ain't it?
-        FetchData(&goal,&nterms,GOAL_IN); 
+        FetchData(&goal,&nterms,GOAL_IN);
 
         /*Allocate output: */

issm/trunk-jpl/src/wrappers/python/io/FetchPythonData.cpp

@@ -32 +32 @@
         }
         else if (PyInt_Check(py_float))
-         dscalar=(double)PyInt_AsLong(py_float);
+                dscalar=(double)PyInt_AsLong(py_float);
         else if (PyBool_Check(py_float))
                 dscalar=(double)PyLong_AsLong(py_float);
@@ -132 +132 @@
         double* dmatrix=NULL;
         double* matrix=NULL;
-        int M,N;
+        int M=0;
+        int N=0;
         int ndim;
         npy_intp*  dims=NULL;
@@ -236 +237 @@
         /*output: */
         int* matrix=NULL;
-        int M,N;
+        int M=0;
+        int N=0;
         int ndim;
         npy_intp*  dims=NULL;
@@ -325 +327 @@
         bool* bmatrix=NULL;
         bool* matrix=NULL;
-        int M,N;
+        int M=0;
+        int N=0;
         int ndim;
         npy_intp*  dims=NULL;
@@ -411 +414 @@
         double* dvector=NULL;
         double* vector=NULL;
-        int M;
+        int M=0;
         int ndim;
         npy_intp*  dims=NULL;
@@ -499 +502 @@
         /*output: */
         float* vector=NULL;
-        int M;
+        int M=0;
         int ndim;
         npy_intp*  dims=NULL;
@@ -539 +542 @@
                                 /*transform into int vector: */
                                 vector=xNew<float>(M);
-                                for(i=0;i<M;i++)vector[i]=(float)lvector[i];
+                                for(i=0;i<M;i++)vector[i]=(float)dvector[i];
                         }
 
@@ -566 +569 @@
                 }
                 else
-                 vector=NULL;
+                        vector=NULL;
         }
         else{
@@ -584 +587 @@
         /*output: */
         int* vector=NULL;
-        int M;
+        int M=0;
         int ndim;
         npy_intp*  dims=NULL;
@@ -624 +627 @@
                                 /*transform into int vector: */
                                 vector=xNew<int>(M);
-                                for(i=0;i<M;i++)vector[i]=(int)lvector[i];
+                                for(i=0;i<M;i++)vector[i]=(int)dvector[i];
                         }
 
@@ -672 +675 @@
         bool* bvector=NULL;
         bool* vector=NULL;
-        int M;
+        int M=0;
         int ndim;
         npy_intp*  dims=NULL;
@@ -923 +926 @@
         *pcontours=contours;
 }
-/*}}}*/
-void FetchChacoData(int* pnvtxs,int** padjacency,int** pstart,float** pewgts,PyObject* A_IN,PyObject* EWGTS_IN){
-        _error_("Nathan... I need your help here");
-}
+
+void FetchChacoData(int* pnvtxs,int** padjacency,int** pstart,float** pewgts,PyObject* A_IN, PyObject* EWGTS_IN){/*{{{*/
+
+        double* a = ((double*)PyArray_DATA((PyArrayObject*)A_IN));
+        int ndim  = PyArray_NDIM((PyArrayObject*)A_IN);
+        int nzmax = PyArray_CountNonzero((PyArrayObject*)A_IN);
+
+        /*Fetch adjacency matrix:*/
+        int      nvtxs       = PyArray_DIMS((PyArrayObject*)A_IN)[1];
+
+        int* mwstart = xNewZeroInit<int>(nvtxs+1);
+        pyGetJc(a,nvtxs,mwstart);
+
+        int* mwadjacency = xNewZeroInit<int>(nzmax);
+        pyGetIr(a,nvtxs,nzmax,mwadjacency);
+        
+        int* start = xNew<int>(nvtxs+1);
+        for(int i=0;i<nvtxs+1;i++) start[i]=(int)mwstart[i];
+        int* adjacency = xNew<int>(nzmax);
+        for(int i=0; i<nzmax; i++) adjacency[i]=(int)mwadjacency[i];
+
+        /*Get edges weights*/
+        int nedges = start[nvtxs];
+        float* ewgts = NULL;
+        int size = PyArray_SIZE((PyArrayObject*)EWGTS_IN);
+        //size may be 1 and still be empty;
+        //presumably size of edge_weights input will never be exactly 1
+        if(size != 1 && size != 0){
+                ewgts = xNewZeroInit<float>(nedges);
+                for(int i = 0; i<nedges;i++){
+                        ewgts[i] = (float)a[i];
+                }
+        }
+
+        /*Assign output pointers*/
+        *pnvtxs     = nvtxs;
+        *padjacency = adjacency;
+        *pstart     = start;
+        *pewgts     = ewgts;
+}
+/*}}}*/
+
+void pyGetJc(double* a, int nvtxs, int* Jc){
+/*{{{*/
+        //get the number of non-zero elements in each row, adding to the prior number;
+        //always starts with 0 and has length: number_of_rows_in(a)+1 == nvtxs+1
+        // eg: 0, 1, 3, 6, 8, 13, ...
+        // for: 1 in the first row, 2 in the second, 3 in the third, etc.
+        int c = 0;
+        Jc[c++] = 0;
+
+        for(int i=0;i<nvtxs;i++){
+                for(int j=0;j<nvtxs;j++){
+                        if ((int)a[i+(j*nvtxs)] != 0){
+                                Jc[c] += 1;
+                        }
+                }
+                c++;
+                Jc[c] = Jc[c-1];
+        }
+        return;
+}
+/*}}}*/
+
+void pyGetIr(double* a, int nvtxs, int nzmax, int* Ir){
+/*{{{*/
+        //get the row-wise position of each non-zero element;
+        //has length: number_of_non_zero_elements_in(a) == nzmax
+        // eg: 10, 24, 25, 4, ...
+        // the 10th, 24th, and 25th elements in the first row, the 4th in the second row
+        // using pyGetJc to determine which indices correspond to which row
+        int r = 0;
+
+        for(int i=0;i<nvtxs;i++){
+                for(int j=0;j<nvtxs;j++){
+                        if ((int)a[i+(j*nvtxs)] != 0){
+                                Ir[r] = j;
+                                r++;
+                        }
+                }
+        }
+        return;
+}
+/*}}}*/
 
 /*Python version dependent: */

issm/trunk-jpl/src/wrappers/python/io/pythonio.h

@@ -51 +51 @@
 void FetchChacoData(int* pnvtxs,int** padjacency,int** pstart,float** pewgts,PyObject* A_IN,PyObject* EWGTS_IN);
 
+void pyGetJc(double* a, int nvtxs, int* Jc);
+void pyGetIr(double* a, int nvtxs, int nzmax, int* Ir);
+
 int CheckNumPythonArguments(PyObject* inputs,int NRHS, void (*function)( void ));