Changeset 8808


Ignore:
Timestamp:
Apr 2, 2013, 4:42:04 PM (12 years ago)
Author:
wilsonp
Message:

Removed netcdf dependence from fit code and compile. Also fixed fitsmooth to use OpenMP when building AtA directly from point coordinates. Changed the default blocking size for reading files in config.py

Location:
trunk/anuga_core/source/anuga
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/config.py

    r8780 r8808  
    201201optimised_gradient_limiter = True # Use hardwired gradient limiter
    202202
    203 points_file_block_line_size = 5000 # Number of lines read in from a points file
     203points_file_block_line_size = 5e7 # Number of lines read in from a points file
    204204                                  # when blocking
    205205

  • trunk/anuga_core/source/anuga/fit_interpolate/fit.py

    r8763 r8808  
    263263            self.Atz = Atz
    264264        else:
    265             fitsmooth.combine_partial_AtA_Atz(self.AtA, AtA,\
     265            fitsmooth.combine_partial_AtA_Atz(self.AtA, AtA, \
    266266                    self.Atz, Atz, zdim, self.mesh.number_of_nodes)
    267 
    268     def _build_matrix_AtA_Atz_file(self, filename, attribute_name=None, max_read_lines=500,\
    269                                    verbose=False):
    270         """Build:
    271         AtA  m x m  interpolation matrix, and,
    272         Atz  m x a  interpolation matrix where,
    273         m is the number of basis functions phi_k (one per vertex)
    274         a is the number of data attributes
    275 
    276         This algorithm uses a quad tree data structure for fast binning of
    277         data points.
    278 
    279         If Ata is None, the matrices AtA and Atz are created.
    280 
    281         This function can be called again and again, with sub-sets of
    282         the point coordinates.  Call fit to get the results.
    283 
    284         Preconditions
    285         z and points are numeric
    286         Point_coordindates and mesh vertices have the same origin.
    287 
    288         The number of attributes of the data points does not change
    289         """
    290         # PADARN NOTE: THe following block of code is translated from
    291         # things that were being done in the Geospatial_data object
    292         # which is no longer required if data from file in C.
    293         #---------------------------------------------------------
    294         # Default attribute name here seems to only have possibility
    295         # of being 'None'.
    296         G_data = Geospatial_data(filename,
    297                                          max_read_lines=1,
    298                                          load_file_now=True,
    299                                          verbose=verbose)
    300 
    301 
    302         if attribute_name is None:
    303             if G_data.default_attribute_name is not None:
    304                 attribute_name = G_data.default_attribute_name
    305             else:
    306                 attribute_name = G_data.attributes.keys()[0]
    307                 # above line takes the first one from keys
    308 
    309         if verbose is True:
    310             log.critical('Using attribute %s' % attribute_name)
    311             log.critical('Available attributes: %s' % (G_data.attributes.keys()))
    312 
    313         msg = 'Attribute name %s does not exist in data set' % attribute_name
    314         assert G_data.attributes.has_key(attribute_name), msg
    315         #---------------------------------------------------------
    316 
    317         # MAKE SURE READING ABSOLUTE POINT COORDINATES
    318         [AtA, Atz, npts] = fitsmooth.build_matrix_AtA_Atz_fileread(self.root.root, \
    319                self.mesh.number_of_nodes, \
    320                self.mesh.triangles, \
    321                filename, \
    322                attribute_name, \
    323                max_read_lines)
    324         self.point_count = npts
    325         self.AtA = AtA
    326         self.Atz = Atz
    327267
    328268    def fit(self, point_coordinates_or_filename=None, z=None,
    … …  
    330270            point_origin=None,
    331271            attribute_name=None,
    332             max_read_lines=500,
    333             use_blocking_option2=True):
     272            max_read_lines=1e7):
    334273        """Fit a smooth surface to given 1d array of data points z.
    335274
    … …  
    350289                use_blocking_option2 = False
    351290
    352         # NOTE PADARN 12/12/12: Currently trying to get all blocking to be done
    353         # in C, but blocking option 1, which does everything using the python
    354         # interface can be used in the meantime (much slower).
    355 
    356         #-----------------------BLOCKING OPTION 1----------------------------
    357         if use_blocking_option2 is False:
     291        # NOTE PADARN 29/03/13: File reading from C has been removed. Now
     292        # the input is either a set of points, or a filename which is then
     293        # handled by the Geospatial_data object
     294
     295        if verbose:
     296            print 'Fit.fit: Initializing'
     297
     298        # Use blocking to load in the point info
     299        if isinstance(point_coordinates_or_filename, basestring):
     300            msg = "Don't set a point origin when reading from a file"
     301            assert point_origin is None, msg
     302            filename = point_coordinates_or_filename
     303
     304            G_data = Geospatial_data(filename,
     305                                     max_read_lines=max_read_lines,
     306                                     load_file_now=False,
     307                                     verbose=verbose)
     308
     309            for i, geo_block in enumerate(G_data):
     310
     311               # Build the array
     312                points = geo_block.get_data_points(absolute=True)
     313                z = geo_block.get_attributes(attribute_name=attribute_name)
     314
     315                self._build_matrix_AtA_Atz(points, z, attribute_name, verbose)
     316
     317            point_coordinates = None
     318
    358319            if verbose:
    359                 print 'Fit.fit: Initializing'
    360 
    361             # Use blocking to load in the point info
    362             if isinstance(point_coordinates_or_filename, basestring):
    363                 msg = "Don't set a point origin when reading from a file"
    364                 assert point_origin is None, msg
    365                 filename = point_coordinates_or_filename
    366 
    367                 G_data = Geospatial_data(filename,
    368                                          max_read_lines=max_read_lines,
    369                                          load_file_now=False,
    370                                          verbose=verbose)
    371 
    372                 for i, geo_block in enumerate(G_data):
    373 
    374                    # Build the array
    375                     points = geo_block.get_data_points(absolute=True)
    376                     z = geo_block.get_attributes(attribute_name=attribute_name)
    377 
    378                     self._build_matrix_AtA_Atz(points, z, attribute_name, verbose)
    379 
    380                 point_coordinates = None
    381 
    382                 if verbose:
    383                     print ''
    384             else:
    385                 point_coordinates = point_coordinates_or_filename
    386 
    387         #-----------------------BLOCKING OPTION 2----------------------------
     320                print ''
    388321        else:
    389             if verbose:
    390                 print 'Fit.fit: Initializing'
    391 
    392             if isinstance(point_coordinates_or_filename, basestring):
    393                 msg = "Don't set a point origin when reading from a file"
    394                 assert point_origin is None, msg
    395                 filename = point_coordinates_or_filename
    396 
    397                 self._build_matrix_AtA_Atz_file(filename, attribute_name=attribute_name,\
    398                                                 verbose=verbose)
    399 
    400                 point_coordinates = None
    401             else:
    402                 point_coordinates = point_coordinates_or_filename
    403         #--------------------------------------------------------------------
    404 
     322            point_coordinates = point_coordinates_or_filename
     323
     324
     325        # This condition either means a filename was read or the function
     326        # recieved a None as input
    405327        if point_coordinates is None:
    406328            if verbose:
    … …  
    410332            assert self.Atz is not None
    411333
    412             # FIXME (DSG) - do  a message
     334
    413335        else:
    414336            point_coordinates = ensure_absolute(point_coordinates,

  • trunk/anuga_core/source/anuga/fit_interpolate/fitsmooth.c

    r8763 r8808  
    1111#include "sparse_dok.h" /*in utilities */
    1212#include "quad_tree.h"
    13 #include <netcdf.h>
    1413#include "omp.h"
    1514
    … …  
    2423    #define PYVERSION273
    2524#endif
    26 
    27 //------------------------ NET CDF READING ----------------------------
    28 
    29 // Note Padarn 06/12/12: I am removing the use of these functions from
    30 // the workflow of the new code. I think they should be reintroduced, but
    31 // to make sure all the unit tests work properly more checking of filetypes
    32 // etc needs to be done.
    33 // Note Padarn 12/12/12: These have been moved back into the workflow,
    34 // but still there is no way of handling .csv or .txt files efficiently.
    35 
    36 
    37 // Reads a block of a netcdf file into the arrays 'point_coordiantes' and 'point_values'. The
    38 // block is specified by a start point (np_start) and and an end point (np_end)
    39 // NOTE PADARN: FIX - only reads elevation at the moment
    40 int _read_net_cdf_points_block(char * filename, char * attribute_name, int np_start,int np_end,double *point_coordinates,double *point_values){
    41 
    42     //point_coordinates and _values are intitialised outside
    43     int np=np_end-np_start+1;
    44     int status;                       /* error status */
    45     int ncid;                          /* netCDF ID */
    46     int rh_id;                         /* variable ID */
    47     size_t start_points[] = {np_start, 0}; /* start at first value */
    48     size_t count_points[] = {np, 2};
    49     size_t start[] = {np_start};
    50     size_t count[] = {np};
    51    
    52     if ((status = nc_open(filename, NC_NOWRITE, &ncid))) ERR(status);
    53 
    54     if ((status = nc_inq_varid(ncid, "points", &rh_id))) ERR(status);
    55 
    56     /* read coordinates from netCDF variable */
    57     if ((status = nc_get_vara_double(ncid, rh_id, start_points, count_points, point_coordinates))) ERR(status);
    58 
    59     if ((status = nc_inq_varid(ncid, attribute_name, &rh_id))) ERR(status);
    60 
    61     /* read values from netCDF variable */
    62     if ((status = nc_get_vara_double(ncid, rh_id, start, count, point_values))) ERR(status);
    63 
    64     if ((status = nc_close(ncid))) ERR(status);
    65 
    66     return 0;
    67 
    68 }
    69 
    70 // Takes a netcdf file and gets the number of points stored in it (also the number of values for
    71 // attributes specified)
    72 // NOTE PADARN: FIX - uses 'points' to get points, this may be a problem.
    73 int _read_net_cdf_entries(char * filename,int * x){
    74 
    75     int status;                       /* error status */
    76     int ncid;
    77     int p_id;
    78     int ndims;
    79     size_t length = {0};
    80 
    81     if ((status = nc_open(filename, NC_NOWRITE, &ncid))) ERR(status);
    82 
    83     if ((status = nc_inq_varid(ncid, "points", &p_id))) ERR(status);
    84    
    85     if((status = nc_inq_varndims(ncid, p_id, &ndims))) ERR(status);
    86    
    87     int * ndimsize = malloc(sizeof(int)*ndims);
    88 
    89     if((status = nc_inq_vardimid(ncid, p_id, ndimsize))) ERR(status);
    90 
    91     if((status = nc_inq_dimlen(ncid, ndimsize[0], &length))) ERR(status);
    92 
    93     *x = (int)length;
    94 
    95     if ((status = nc_close(ncid))) ERR(status);
    96 
    97     return 0;
    98 
    99 }
    100 
    101 // Takes a netcdf file and gets the x and y corners
    102 int _read_net_cdf_corners(char * filename,double *xllcorner, double *yllcorner){
    103 
    104     int status;                       /* error status */
    105     int ncid;
    106 
    107     if ((status = nc_open(filename, NC_NOWRITE, &ncid))) ERR(status);
    108 
    109     if ((status = nc_get_att_double(ncid, NC_GLOBAL, "xllcorner", xllcorner))) ERR(status);
    110     if ((status = nc_get_att_double(ncid, NC_GLOBAL, "yllcorner", yllcorner))) ERR(status);
    111 
    112     if ((status = nc_close(ncid))) ERR(status);
    113 
    114     return 0;
    115 
    116 }
    117 
    118 //--------------------------------------------------------------------------
    11925
    12026//-------------------------- QUANTITY FITTING ------------------------------
    … …  
    272178// Builds the AtA and Atz interpolation matrix
    273179// and residual. Uses a quad_tree for fast access to the triangles of the mesh.
    274 // This function takes a filename and an attribute name, in netcdf format .pts,
    275 // and reads the attribute value and point location from this file.
    276 int _build_matrix_AtA_Atz_fileread(int N, long * triangles,
    277                       char * filename,
    278                       char * attribute_name,
    279                       int blocksize,
    280                       sparse_dok * AtA,
    281                       double * Atz,quad_tree * quadtree)
    282               {
    283 
    284 
    285     int np;
    286     _read_net_cdf_entries(filename,&np);
    287     //n=10;
    288 
    289     double xllcorner, yllcorner;
    290     _read_net_cdf_corners(filename,&xllcorner,&yllcorner);
    291 
    292     //printf(" xllcorner %g yllcorner %g\n",xllcorner,yllcorner);
    293 
    294 
    295     double * point_coordinates = malloc(2*sizeof(double)*blocksize);
    296     double * point_values = malloc(sizeof(double)*blocksize);
    297     int k;
    298     int left;
    299 
    300     int i;
    301     for(i=0;i<N;i++){
    302         Atz[i]=0;
    303     }
    304 
    305     edge_key_t key;
    306     int w;
    307     int start = 0;
    308     left = np;
    309 
    310    
    311     while(left > 0){
    312         if(left-blocksize<=0){
    313             blocksize = left;
    314         }
    315 
    316         // Read data into arrays
    317         _read_net_cdf_points_block(filename,attribute_name,start,
    318         start+blocksize-1,point_coordinates,point_values);
    319        
    320         start = start + blocksize;
    321 
    322             #pragma omp parallel for private(k,i,key,w)
    323             for(k=0;k<blocksize;k++){
    324 
    325 
    326                 double x = point_coordinates[2*k]+xllcorner;
    327                 double y = point_coordinates[2*k+1]+yllcorner;
    328                 triangle * T = search(quadtree,x,y);
    329 
    330                 if(T!=NULL){
    331                     double * sigma = calculate_sigma(T,x,y);
    332                     int js[3];
    333                     for(i=0;i<3;i++){
    334                         js[i]=triangles[3*(T->index)+i];
    335                     }
    336                    
    337                     #pragma omp critical
    338                     {
    339                     for(i=0;i<3;i++){
    340 
    341                         Atz[js[i]] += sigma[i]*point_values[k];
    342 
    343                         for(w=0;w<3;w++){
    344                            
    345                             key.i=js[i];
    346                             key.j=js[w];
    347 
    348                             add_dok_entry(AtA,key,sigma[i]*sigma[w]);
    349                         }                       
    350                     }
    351                     }
    352                     free(sigma);
    353                     sigma=NULL;
    354                 }
    355 
    356             }
    357        
    358         left=left-blocksize;
    359     }
    360    
    361     return 0;
    362 }
    363 
    364 // Builds the AtA and Atz interpolation matrix
    365 // and residual. Uses a quad_tree for fast access to the triangles of the mesh.
    366180// This function takes a list of point coordinates, and associated point values
    367181// (for any number of attributes).
    … …  
    391205
    392206
    393     //#pragma omp parallel for private(k,i,key,w)
     207    #pragma omp parallel for private(k,i,key,w)
    394208    for(k=0;k<npts;k++){
    395209
    … …  
    651465   
    652466    #endif
    653 
    654 }
    655 
    656 // Read a data file (netcdf format) to build
    657 // AtA and Atz. Returns a pointer to the sparse dok matrix AtA wrapped in a capsule
    658 // object and a python list for the array Atz.
    659 PyObject *build_matrix_AtA_Atz_fileread(PyObject *self, PyObject *args) {
    660 
    661 
    662     PyArrayObject *triangles;
    663     char * file_name;
    664     char * attribute_name;
    665     int block_size;
    666     int N; // Number of triangles
    667     int err;
    668     PyObject *tree;
    669 
    670     // Convert Python arguments to C
    671     if (!PyArg_ParseTuple(args, "OiOssi",&tree, &N,
    672                                             &triangles,
    673                                             &file_name,
    674                                             &attribute_name,
    675                                             &block_size
    676                                             )) {
    677       PyErr_SetString(PyExc_RuntimeError,
    678               "fitsmooth.c: could not parse input");
    679       return NULL;
    680     }
    681 
    682    
    683 
    684     CHECK_C_CONTIG(triangles);
    685 
    686     #ifdef PYVERSION273
    687     quad_tree * quadtree = (quad_tree*) PyCapsule_GetPointer(tree,"quad tree");
    688     #else
    689     quad_tree * quadtree = (quad_tree*) PyCObject_AsVoidPtr(tree);
    690     #endif
    691 
    692     sparse_dok * dok_AtA; // Should be an input argument?
    693     dok_AtA = make_dok();
    694     double * Atz = malloc(sizeof(double)*N);
    695 
    696     err = _build_matrix_AtA_Atz_fileread(N,(long*) triangles->data,
    697                       file_name,
    698                       attribute_name,
    699                       block_size,
    700                       dok_AtA,
    701                       Atz,
    702                       quadtree);
    703 
    704     if (err != 0) {
    705       PyErr_SetString(PyExc_RuntimeError,
    706               "Unknown Error");
    707       return NULL;
    708     }
    709 
    710     #ifdef PYVERSION273
    711     PyObject * AtA_cap =  PyCapsule_New((void*) dok_AtA,
    712                   "sparse dok",
    713                   &delete_dok_cap);
    714     #else
    715     PyObject * AtA_cap =  PyCObject_FromVoidPtr((void*) dok_AtA,
    716                   &delete_dok_cobj);
    717     #endif
    718 
    719     PyObject *Atz_ret = c_double_array_to_list(Atz,N);
    720 
    721     int npts;
    722     err = _read_net_cdf_entries(file_name, &npts);
    723 
    724     PyObject *lst = PyList_New(3);
    725     PyList_SET_ITEM(lst, 0, AtA_cap);
    726     PyList_SET_ITEM(lst, 1, Atz_ret);
    727     PyList_SET_ITEM(lst, 2, PyInt_FromLong((long)npts));
    728     return lst;
    729467
    730468}
    … …  
    1121859    {"build_quad_tree",build_quad_tree, METH_VARARGS, "Print out"},
    1122860    {"build_smoothing_matrix",build_smoothing_matrix, METH_VARARGS, "Print out"},
    1123     {"build_matrix_AtA_Atz_fileread",build_matrix_AtA_Atz_fileread, METH_VARARGS, "Print out"},
    1124861    {"build_matrix_AtA_Atz_points",build_matrix_AtA_Atz_points, METH_VARARGS, "Print out"},
    1125862    {"combine_partial_AtA_Atz",combine_partial_AtA_Atz, METH_VARARGS, "Print out"},

  • trunk/anuga_core/source/anuga/utilities/compile.py

    r8787 r8808  
    1515import os, string, sys
    1616
    17 
    18 NETCDF_LIB_DIR = os.getenv('NETCDF_LIB_DIR', '')
    19 NETCDF_INCLUDE_DIR = os.getenv('NETCDF_INCLUDE_DIR', '')
    20 
    21 #print 'NETCDF_LIB_DIR: ',NETCDF_LIB_DIR
    22 #print 'NETCDF_INCLUDE_DIR: ',NETCDF_INCLUDE_DIR
    23 
    2417I_dirs = ''
    25 if NETCDF_INCLUDE_DIR != '' :
    26     I_dirs = ' -I"%s" ' % NETCDF_INCLUDE_DIR
    27 
    28 netcdf_lib_dirs = ''
    29 if NETCDF_LIB_DIR != '' :
    30     netcdf_lib_dirs = ' -L"%s" ' % NETCDF_LIB_DIR
    31 
    32 #print 'netcdf_lib_dirs: ',netcdf_lib_dirs
    3318
    3419#NumPy ------------------------------------
    … …  
    333318  # Make shared library (*.so or *.dll)
    334319  if FN=="fitsmooth.c":
    335     fitlibs = libs + netcdf_lib_dirs
     320    fitlibs = libs
    336321    if sys.platform == 'win32':   
    337322      if fitlibs is "":
    338         s = '%s -%s %s ../utilities/quad_tree.o ../utilities/sparse_dok.o ../utilities/sparse_csr.o -o %s.%s -lm  -fopenmp netcdf.dll' %(loader, sharedflag, object_files, root1, libext)
     323        s = '%s -%s %s ../utilities/quad_tree.o ../utilities/sparse_dok.o ../utilities/sparse_csr.o -o %s.%s -lm  -fopenmp' %(loader, sharedflag, object_files, root1, libext)
    339324      else:
    340         s = '%s -%s %s ../utilities/quad_tree.o ../utilities/sparse_dok.o ../utilities/sparse_csr.o -o %s.%s %s -lm  -fopenmp netcdf.dll' %(loader, sharedflag, object_files, root1, libext, fitlibs)
     325        s = '%s -%s %s ../utilities/quad_tree.o ../utilities/sparse_dok.o ../utilities/sparse_csr.o -o %s.%s %s -lm  -fopenmp' %(loader, sharedflag, object_files, root1, libext, fitlibs)
    341326    else:   
    342327      if fitlibs is "":
    343         s = '%s -%s %s ../utilities/quad_tree.o ../utilities/sparse_dok.o ../utilities/sparse_csr.o -o %s.%s -lm  -fopenmp -lnetcdf' %(loader, sharedflag, object_files, root1, libext)
     328        s = '%s -%s %s ../utilities/quad_tree.o ../utilities/sparse_dok.o ../utilities/sparse_csr.o -o %s.%s -lm  -fopenmp' %(loader, sharedflag, object_files, root1, libext)
    344329      else:
    345         s = '%s -%s %s ../utilities/quad_tree.o ../utilities/sparse_dok.o ../utilities/sparse_csr.o -o %s.%s %s -lm  -fopenmp -lnetcdf' %(loader, sharedflag, object_files, root1, libext, fitlibs)
     330        s = '%s -%s %s ../utilities/quad_tree.o ../utilities/sparse_dok.o ../utilities/sparse_csr.o -o %s.%s %s -lm  -fopenmp' %(loader, sharedflag, object_files, root1, libext, fitlibs)
    346331  elif FN=="quad_tree_ext.c":
    347332    if libs is "":
Note: See TracChangeset for help on using the changeset viewer.