Continuing migration of updated pygeosx tools

Author: cssherman

pygeosx_tools_package/pygeosx_tools/geophysics/fiber.py

@@ -7,7 +7,8 @@
 
 
 class Fiber():
-    def __init__(self):
+
+    def __init__( self ):
         self.time = []
         self.channel_position = []
         self.gage_length = -1.0
@@ -19,9 +20,9 @@ def __init__(self):
 
 
 class FiberAnalysis():
-    def __init__(self):
+
+    def __init__( self ):
         """
         InSAR Analysis class
         """
         self.set_names = []
-

pygeosx_tools_package/pygeosx_tools/geophysics/insar.py

@@ -0,0 +1,86 @@
+from mpi4py import MPI
+import numpy as np
+
+# Get the MPI rank
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+
+
+def get_global_array_range( local_values ):
+    # 1D arrays will return a scalar, ND arrays an array
+    N = np.shape( local_values )
+    local_min = 1e100
+    local_max = -1e100
+    if ( len( N ) > 1 ):
+        local_min = np.zeros( N[ 1 ] ) + 1e100
+        local_max = np.zeros( N[ 1 ] ) - 1e100
+
+    # For >1D arrays, keep the last dimension
+    query_axis = 0
+    if ( len( N ) > 2 ):
+        query_axis = tuple( [ ii for ii in range( 0, len( N ) - 1 ) ] )
+
+    # Ignore zero-length results
+    if len( local_values ):
+        local_min = np.amin( local_values, axis=query_axis )
+        local_max = np.amax( local_values, axis=query_axis )
+
+    # Gather the results onto rank 0
+    all_min = comm.gather( local_min, root=0 )
+    all_max = comm.gather( local_max, root=0 )
+    global_min = 1e100
+    global_max = -1e100
+    if ( rank == 0 ):
+        global_min = np.amin( np.array( all_min ), axis=0 )
+        global_max = np.amax( np.array( all_max ), axis=0 )
+
+    # Broadcast
+    global_min = comm.bcast( global_min, root=0 )
+    global_max = comm.bcast( global_max, root=0 )
+
+    return global_min, global_max
+
+
+def gather_array( local_values, allgather=False, concatenate=True ):
+    # Find buffer size
+    N = np.shape( local_values )
+    M = np.prod( N )
+    all_M = []
+    max_M = 0
+    if allgather:
+        all_M = comm.allgather( M )
+        max_M = np.amax( all_M )
+    else:
+        all_M = comm.gather( M, root=0 )
+        if ( rank == 0 ):
+            max_M = np.amax( all_M )
+        max_M = comm.bcast( max_M, root=0 )
+
+    # Pack the array into a buffer
+    send_buff = np.zeros( max_M )
+    send_buff[ :M ] = np.reshape( local_values, ( -1 ) )
+    receive_buff = np.zeros( ( comm.size, max_M ) )
+
+    # Gather the buffers
+    if allgather:
+        comm.Allgather( [ send_buff, MPI.DOUBLE ], [ receive_buff, MPI.DOUBLE ] )
+    else:
+        comm.Gather( [ send_buff, MPI.DOUBLE ], [ receive_buff, MPI.DOUBLE ], root=0 )
+
+    # Unpack the buffers
+    all_values = []
+    R = list( N )
+    R[ 0 ] = -1
+    if ( ( rank == 0 ) | allgather ):
+        # Reshape each rank's contribution
+        for ii in range( comm.size ):
+            if ( all_M[ ii ] > 0 ):
+                tmp = np.reshape( receive_buff[ ii, :all_M[ ii ] ], R )
+                all_values.append( tmp )
+
+        # Concatenate into a single array
+        if concatenate:
+            if ( len( all_values ) ):
+                all_values = np.concatenate( all_values, axis=0 )
+
+    return all_values

pygeosx_tools_package/pygeosx_tools/geophysics/microseismic.py

@@ -0,0 +1,32 @@
+def patch_numpy_where():
+    """
+    Some packages (numpy, etc.) use np.where in a way that isn't
+    compliant with the API.  In most cases this is fine, but when the
+    geosx environment is initialized this can cause critical errors with
+    nan inputs.  This patch checks and updates the inputs to where
+    """
+    import numpy as np
+    np.where_original_fn = np.where
+    print( 'patching np.where', flush=True )
+
+    def variable_as_array_type( x, target_shape ):
+        if x is None:
+            return
+
+        if isinstance( x, ( np.ndarray, list ) ):
+            return x
+
+        y = np.empty( target_shape )
+        y[...] = x
+        return y
+
+    def flexible_where( condition, x=None, y=None ):
+        s = np.shape( condition )
+        x = variable_as_array_type( x, s )
+        y = variable_as_array_type( y, s )
+        return np.where_original_fn( condition, x, y )
+
+    np.where = flexible_where
+
+
+patch_numpy_where()

pygeosx_tools_package/pygeosx_tools/parallel_io.py

@@ -0,0 +1,115 @@
+import os
+import warnings
+import matplotlib.pyplot as plt
+from matplotlib import cm
+import matplotlib.colors as mcolors
+
+
+class HighResPlot():
+
+    def __init__( self ):
+        self.handle = plt.figure()
+        self.isAdjusted = False
+        self.colorspace = 'RGB'
+        self.compress = False
+        self.output_format = 'png'
+        self.font_weight = 'normal'
+        self.axis_font_size = 8
+        self.legend_font_size = 8
+        self.size = ( 3.5, 2.625 )
+        self.dpi = 200
+        self.apply_tight_layout = True
+        self.disable_antialiasing = False
+        self.use_imagemagick = False
+        self.apply_font_size()
+
+    def __del__( self ):
+        try:
+            self.handle.clf()
+            plt.close( self.handle.number )
+        except:
+            pass
+
+    def set_active( self ):
+        plt.figure( self.handle.number )
+
+    def reset( self ):
+        tmp = self.handle.number
+        self.set_active()
+        self.handle.clf()
+        del self.handle
+        self.handle = plt.figure( tmp )
+        self.apply_font_size()
+        self.isAdjusted = False
+
+    def apply_font_size( self ):
+        self.font = { 'weight': self.font_weight, 'size': self.axis_font_size }
+        plt.rc( 'font', **self.font )
+
+    def apply_format( self ):
+        self.set_active()
+        self.handle.set_size_inches( self.size )
+        self.apply_font_size()
+
+        for ax in self.handle.get_axes():
+            tmp = ax.get_legend()
+            if tmp:
+                ax.legend( loc=tmp._loc, fontsize=self.legend_font_size )
+
+        if self.apply_tight_layout:
+            with warnings.catch_warnings():
+                warnings.simplefilter( 'ignore' )
+                self.handle.tight_layout( pad=1.5 )
+                self.apply_tight_layout = False
+
+    def save( self, fname ):
+        self.apply_format()
+        if plt.isinteractive():
+            plt.draw()
+            plt.pause( 0.1 )
+
+        if self.use_imagemagick:
+            # Matplotlib's non-pdf layout is sometimes odd,
+            # so use the imagemagick convert as a back-up option
+            plt.savefig( '%s.pdf' % ( fname ), format='pdf', dpi=self.dpi )
+
+            if ( self.output_format != 'pdf' ):
+                cmd = 'convert -density %i %s.pdf' % ( self.dpi, fname )
+
+                if ( self.colorspace == 'CMYK' ):
+                    cmd += ' -colorspace CMYK'
+                if ( self.colorspace == 'gray' ):
+                    cmd += ' -colorspace gray'
+                if ( self.compress ):
+                    cmd += ' -compress lzw'
+                if ( self.disable_antialiasing ):
+                    cmd += ' +antialias'
+
+                cmd += ' %s.%s' % ( fname, self.output_format )
+                os.system( cmd )
+                os.system( 'rm %s.pdf' % fname )
+        else:
+            plt.savefig( '%s.%s' % ( fname, self.output_format ), format=self.output_format, dpi=self.dpi )
+
+
+def get_modified_jet_colormap():
+    # Add an updated colormap
+    cdict = cm.get_cmap( 'jet' ).__dict__[ '_segmentdata' ].copy()
+    for k in cdict.keys():
+        tmp_seq = list( cdict[ k ] )
+        tmp_final = list( tmp_seq[ 0 ] )
+        tmp_final[ 1 ] = 1.0
+        tmp_final[ 2 ] = 1.0
+        tmp_seq[ 0 ] = tuple( tmp_final )
+        cdict[ k ] = tuple( tmp_seq )
+    return mcolors.LinearSegmentedColormap( 'jet_mod', cdict )
+
+
+def get_periodic_line_style( ii ):
+    col = [ 'k', 'b', 'c', 'g', 'y', 'r', 'm' ]
+    mark = [ '-', '--', ':', '-.', 'o', '*', '^', '+' ]
+
+    jj = ii % ( len( col ) * len( mark ) )
+    style = col[ jj % len( col ) ] + mark[ int( jj / len( col ) ) ]
+
+    return style

pygeosx_tools_package/pygeosx_tools/patches.py

@@ -1,13 +1,9 @@
 import sys
 import numpy as np
-from mpi4py import MPI
 import matplotlib.pyplot as plt
 import pylvarray
 import pygeosx
-
-# Get the MPI rank
-comm = MPI.COMM_WORLD
-rank = comm.Get_rank()
+from pygeosx_tools import parallel_io
 
 
 def get_wrapper( problem, target_key, write_flag=False ):
@@ -52,7 +48,7 @@ def get_wrapper_par( problem, target_key, allgather=False, ghost_key='' ):
     Returns:
         np.ndarray: The wrapper as a numpy ndarray
     """
-    if ( comm.size == 1 ):
+    if ( parallel_io.comm.size == 1 ):
         # This is a serial problem
         return get_wrapper( problem, target_key )
 
@@ -66,45 +62,7 @@ def get_wrapper_par( problem, target_key, allgather=False, ghost_key='' ):
             ghost_values = get_wrapper( problem, ghost_key )
             local_values = local_values[ ghost_values < -0.5 ]
 
-        # Find buffer size
-        N = np.shape( local_values )
-        M = np.prod( N )
-        all_M = []
-        max_M = 0
-        if allgather:
-            all_M = comm.allgather( M )
-            max_M = np.amax( all_M )
-        else:
-            all_M = comm.gather( M, root=0 )
-            if ( rank == 0 ):
-                max_M = np.amax( all_M )
-            max_M = comm.bcast( max_M, root=0 )
-
-        # Pack the array into a buffer
-        send_buff = np.zeros( max_M )
-        send_buff[ :M ] = np.reshape( local_values, ( -1 ) )
-        receive_buff = np.zeros( ( comm.size, max_M ) )
-
-        # Gather the buffers
-        if allgather:
-            comm.Allgather( [ send_buff, MPI.DOUBLE ], [ receive_buff, MPI.DOUBLE ] )
-        else:
-            comm.Gather( [ send_buff, MPI.DOUBLE ], [ receive_buff, MPI.DOUBLE ], root=0 )
-
-        # Unpack the buffers
-        all_values = []
-        R = list( N )
-        R[ 0 ] = -1
-        if ( ( rank == 0 ) | allgather ):
-            # Reshape each rank's contribution
-            for ii in range( comm.size ):
-                if ( all_M[ ii ] > 0 ):
-                    tmp = np.reshape( receive_buff[ ii, :all_M[ ii ] ], R )
-                    all_values.append( tmp )
-
-            # Concatenate into a single array
-            all_values = np.concatenate( all_values, axis=0 )
-        return all_values
+        return parallel_io.gather_array( local_values, allgather=allgather )
 
 
 def gather_wrapper( problem, key, ghost_key='' ):
@@ -147,34 +105,7 @@ def get_global_value_range( problem, key ):
         tuple: The global min/max of the target
     """
     local_values = get_wrapper( problem, key )
-
-    # 1D arrays will return a scalar, ND arrays an array
-    N = np.shape( local_values )
-    local_min = 1e100
-    local_max = -1e100
-    if ( len( N ) > 1 ):
-        local_min = np.zeros( N[ 1 ] ) + 1e100
-        local_max = np.zeros( N[ 1 ] ) - 1e100
-
-    # For >1D arrays, keep the last dimension
-    query_axis = 0
-    if ( len( N ) > 2 ):
-        query_axis = tuple( [ ii for ii in range( 0, len( N ) - 1 ) ] )
-
-    # Ignore zero-length results
-    if len( local_values ):
-        local_min = np.amin( local_values, axis=query_axis )
-        local_max = np.amax( local_values, axis=query_axis )
-
-    # Gather the results onto rank 0
-    all_min = comm.gather( local_min, root=0 )
-    all_max = comm.gather( local_max, root=0 )
-    global_min = 1e100
-    global_max = -1e100
-    if ( rank == 0 ):
-        global_min = np.amin( np.array( all_min ), axis=0 )
-        global_max = np.amax( np.array( all_max ), axis=0 )
-    return global_min, global_max
+    return parallel_io.get_global_array_range( local_values )
 
 
 def print_global_value_range( problem, key, header, scale=1.0, precision='%1.4f' ):
@@ -195,7 +126,7 @@ def print_global_value_range( problem, key, header, scale=1.0, precision='%1.4f'
     global_min *= scale
     global_max *= scale
 
-    if ( rank == 0 ):
+    if ( parallel_io.rank == 0 ):
         if isinstance( global_min, np.ndarray ):
             min_str = ', '.join( [ precision % ( x ) for x in global_min ] )
             max_str = ', '.join( [ precision % ( x ) for x in global_max ] )
@@ -313,12 +244,12 @@ def get_matching_wrapper_path( problem, filters ):
     search_datastructure_wrappers_recursive( problem, filters, matching_paths )
 
     if ( len( matching_paths ) == 1 ):
-        if ( rank == 0 ):
+        if ( parallel_io.rank == 0 ):
             print( 'Found matching wrapper: %s' % ( matching_paths[ 0 ] ) )
         return matching_paths[ 0 ]
 
     else:
-        if ( rank == 0 ):
+        if ( parallel_io.rank == 0 ):
             print( 'Error occured while looking for wrappers:' )
             print( 'Filters: [%s]' % ( ', '.join( filters ) ) )
             print( 'Matching wrappers: [%s]' % ( ', '.join( matching_paths ) ) )
@@ -360,7 +291,7 @@ def plot_history( records, output_root='.', save_figures=True, show_figures=True
         save_figures (bool): Flag to indicate whether figures should be saved (default = True)
         show_figures (bool): Flag to indicate whether figures should be drawn (default = False)
     """
-    if ( rank == 0 ):
+    if ( parallel_io.rank == 0 ):
         for k in records.keys():
             if ( k != 'time' ):
                 # Set the active figure