Index: /issm/trunk-jpl/src/m/classes/issmsettings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/issmsettings.py (revision 23771) +++ /issm/trunk-jpl/src/m/classes/issmsettings.py (revision 23772) @@ -3,92 +3,95 @@ from WriteData import WriteData + class issmsettings(object): - """ - ISSMSETTINGS class definition + """ + ISSMSETTINGS class definition - Usage: - issmsettings=issmsettings(); - """ + Usage: + issmsettings=issmsettings(); + """ - def __init__(self): # {{{ - self.results_on_nodes = [] - self.io_gather = 0 - self.lowmem = 0 - self.output_frequency = 0 - self.coupling_frequency = 0 - self.recording_frequency = 0 - self.waitonlock = 0 - self.solver_residue_threshold = 0 + def __init__(self): # {{{ + self.results_on_nodes = [] + self.io_gather = 0 + self.lowmem = 0 + self.output_frequency = 0 + self.coupling_frequency = 0 + self.recording_frequency = 0 + self.waitonlock = 0 + self.solver_residue_threshold = 0 - #set defaults - self.setdefaultparameters() + #set defaults + self.setdefaultparameters() - #}}} - def __repr__(self): # {{{ - string=" general issmsettings parameters:" + #}}} + def __repr__(self): # {{{ + string = " general issmsettings parameters:" - string="%s\n%s"%(string,fielddisplay(self,"results_on_nodes","list of output for which results will be output for all the nodes of each element, Use 'all' for all output on nodes.")) - string="%s\n%s"%(string,fielddisplay(self,"io_gather","I/O gathering strategy for result outputs (default 1)")) - string="%s\n%s"%(string,fielddisplay(self,"lowmem","is the memory limited ? (0 or 1)")) - string="%s\n%s"%(string,fielddisplay(self,"output_frequency","frequency at which results are saved in all solutions with multiple time_steps")) - string="%s\n%s"%(string,fielddisplay(self,"sb_coupling_frequency","frequency at which StressBalance solver is coupled (default 1)")) - string="%s\n%s"%(string,fielddisplay(self,"recording_frequency","frequency at which the runs are being recorded, allowing for a restart")) - string="%s\n%s"%(string,fielddisplay(self,"waitonlock","maximum number of minutes to wait for batch results, or return 0")) - string="%s\n%s"%(string,fielddisplay(self,"solver_residue_threshold","throw an error if solver residue exceeds this value (NaN to deactivate)")) - return string - #}}} - def setdefaultparameters(self): # {{{ - - #are we short in memory ? (0 faster but requires more memory) - self.lowmem=0 + string = "%s\n%s" % (string, fielddisplay(self, "results_on_nodes", "list of output for which results will be output for all the nodes of each element, Use 'all' for all output on nodes.")) + string = "%s\n%s" % (string, fielddisplay(self, "io_gather", "I/O gathering strategy for result outputs (default 1)")) + string = "%s\n%s" % (string, fielddisplay(self, "lowmem", "is the memory limited ? (0 or 1)")) + string = "%s\n%s" % (string, fielddisplay(self, "output_frequency", "frequency at which results are saved in all solutions with multiple time_steps")) + string = "%s\n%s" % (string, fielddisplay(self, "sb_coupling_frequency", "frequency at which StressBalance solver is coupled (default 1)")) + string = "%s\n%s" % (string, fielddisplay(self, "recording_frequency", "frequency at which the runs are being recorded, allowing for a restart")) + string = "%s\n%s" % (string, fielddisplay(self, "waitonlock", "maximum number of minutes to wait for batch results, or return 0")) + string = "%s\n%s" % (string, fielddisplay(self, "solver_residue_threshold", "throw an error if solver residue exceeds this value (NaN to deactivate)")) + return string + #}}} - #i/o: - self.io_gather=1 + def setdefaultparameters(self): # {{{ - #results frequency by default every step - self.output_frequency=1 + #are we short in memory ? (0 faster but requires more memory) + self.lowmem = 0 - #coupling frequency of the stress balance solver by default every step - self.sb_coupling_frequency=1 - - #checkpoints frequency, by default never: - self.recording_frequency=0 + #i/o: + self.io_gather = 1 + #results frequency by default every step + self.output_frequency = 1 - #this option can be activated to load automatically the results - #onto the model after a parallel run by waiting for the lock file - #N minutes that is generated once the solution has converged - #0 to deactivate - self.waitonlock=2**31-1 + #coupling frequency of the stress balance solver by default every step + self.sb_coupling_frequency = 1 - #throw an error if solver residue exceeds this value - self.solver_residue_threshold=1e-6; + #checkpoints frequency, by default never: + self.recording_frequency = 0 - return self - #}}} - def checkconsistency(self,md,solution,analyses): # {{{ - md = checkfield(md,'fieldname','settings.results_on_nodes','stringrow',1) - md = checkfield(md,'fieldname','settings.io_gather','numel',[1],'values',[0,1]) - md = checkfield(md,'fieldname','settings.lowmem','numel',[1],'values',[0,1]) - md = checkfield(md,'fieldname','settings.output_frequency','numel',[1],'>=',1) - md = checkfield(md,'fieldname','settings.sb_coupling_frequency','numel',[1],'>=',1) - md = checkfield(md,'fieldname','settings.recording_frequency','numel',[1],'>=',0) - md = checkfield(md,'fieldname','settings.waitonlock','numel',[1]) - md = checkfield(md,'fieldname','settings.solver_residue_threshold','numel',[1],'>',0) + #this option can be activated to load automatically the results + #onto the model after a parallel run by waiting for the lock file + #N minutes that is generated once the solution has converged + #0 to deactivate + self.waitonlock = 2 ** 31 - 1 - return md - # }}} - def marshall(self,prefix,md,fid): # {{{ - WriteData(fid,prefix,'data',self.results_on_nodes,'name','md.settings.results_on_nodes','format','StringArray') - WriteData(fid,prefix,'object',self,'class','settings','fieldname','io_gather','format','Boolean') - WriteData(fid,prefix,'object',self,'class','settings','fieldname','lowmem','format','Boolean') - WriteData(fid,prefix,'object',self,'class','settings','fieldname','output_frequency','format','Integer') - WriteData(fid,prefix,'object',self,'class','settings','fieldname','sb_coupling_frequency','format','Integer') - WriteData(fid,prefix,'object',self,'class','settings','fieldname','recording_frequency','format','Integer') - - if self.waitonlock>0: - WriteData(fid,prefix,'name','md.settings.waitonlock','data',True,'format','Boolean') - else: - WriteData(fid,prefix,'name','md.settings.waitonlock','data',False,'format','Boolean') - WriteData(fid,prefix,'object',self,'fieldname','solver_residue_threshold','format','Double') - # }}} + #throw an error if solver residue exceeds this value + self.solver_residue_threshold = 1e-6 + + return self + #}}} + + def checkconsistency(self, md, solution, analyses): # {{{ + md = checkfield(md, 'fieldname', 'settings.results_on_nodes', 'stringrow', 1) + md = checkfield(md, 'fieldname', 'settings.io_gather', 'numel', [1], 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'settings.lowmem', 'numel', [1], 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'settings.output_frequency', 'numel', [1], '>=', 1) + md = checkfield(md, 'fieldname', 'settings.sb_coupling_frequency', 'numel', [1], '>=', 1) + md = checkfield(md, 'fieldname', 'settings.recording_frequency', 'numel', [1], '>=', 0) + md = checkfield(md, 'fieldname', 'settings.waitonlock', 'numel', [1]) + md = checkfield(md, 'fieldname', 'settings.solver_residue_threshold', 'numel', [1], '>', 0) + + return md + # }}} + + def marshall(self, prefix, md, fid): # {{{ + WriteData(fid, prefix, 'data', self.results_on_nodes, 'name', 'md.settings.results_on_nodes', 'format', 'StringArray') + WriteData(fid, prefix, 'object', self, 'class', 'settings', 'fieldname', 'io_gather', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'class', 'settings', 'fieldname', 'lowmem', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'class', 'settings', 'fieldname', 'output_frequency', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'class', 'settings', 'fieldname', 'sb_coupling_frequency', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'class', 'settings', 'fieldname', 'recording_frequency', 'format', 'Integer') + + if self.waitonlock > 0: + WriteData(fid, prefix, 'name', 'md.settings.waitonlock', 'data', True, 'format', 'Boolean') + else: + WriteData(fid, prefix, 'name', 'md.settings.waitonlock', 'data', False, 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'fieldname', 'solver_residue_threshold', 'format', 'Double') + # }}} Index: /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py (revision 23771) +++ /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py (revision 23772) @@ -1,212 +1,202 @@ -from netCDF4 import Dataset, stringtochar +from netCDF4 import Dataset import numpy as np import time import collections -from mesh2d import * -from mesh3dprisms import * -from results import * from os import path, remove -def export_netCDF(md,filename): - #Now going on Real treatment - if path.exists(filename): - print(('File {} allready exist'.format(filename))) - newname=eval(input('Give a new name or "delete" to replace: ')) - if newname=='delete': - remove(filename) - else: - print(('New file name is {}'.format(newname))) - filename=newname - - NCData=Dataset(filename, 'w', format='NETCDF4') - NCData.description = 'Results for run' + md.miscellaneous.name - NCData.history = 'Created ' + time.ctime(time.time()) - - #define netCDF dimensions - try: - StepNum=np.shape(dict.values(md.results.__dict__))[1] - except IndexError: - StepNum=1 - Dimension1=NCData.createDimension('DimNum1',StepNum)#time is first - DimDict={len(Dimension1):'DimNum1'} - dimindex=1 - - dimlist=[2,md.mesh.numberofelements,md.mesh.numberofvertices,np.shape(md.mesh.elements)[1]] - for i in range(0,4): - if dimlist[i] not in list(DimDict.keys()): - dimindex+=1 - NewDim=NCData.createDimension('DimNum'+str(dimindex),dimlist[i]) - DimDict[len(NewDim)]='DimNum'+str(dimindex) - - typelist=[bool,str,str,int,float,complex, - collections.OrderedDict, - np.int64,np.ndarray,np.float64] - groups=dict.keys(md.__dict__) - #get all model classes and create respective groups - for group in groups: - NCgroup=NCData.createGroup(str(group)) - #In each group gather the fields of the class - fields=dict.keys(md.__dict__[group].__dict__) - - #looping on fields - for field in fields: - #Special treatment for list fields - if type(md.__dict__[group].__dict__[field])==list: - StdList=False - if len(md.__dict__[group].__dict__[field])==0: - StdList=True - else: - StdList=type(md.__dict__[group].__dict__[field][0]) in typelist - NCgroup.__setattr__('classtype', md.__dict__[group].__class__.__name__) - if StdList: #this is a standard or empty list just proceed - Var=md.__dict__[group].__dict__[field] - DimDict=CreateVar(NCData,Var,field,NCgroup,DimDict) - else: #this is a list of fields, specific treatment needed - Listsize=len(md.__dict__[group].__dict__[field]) - Subgroup=NCgroup.createGroup(str(field)) - Subgroup.__setattr__('classtype',md.__dict__[group].__dict__[field].__class__.__name__) - for listindex in range(0,Listsize): - try: - Listgroup=Subgroup.createGroup(str(md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__['name'])) - except KeyError: - for naming in ['step']: - Listgroup=Subgroup.createGroup(str(md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__[naming])) - except AttributeError: - Listgroup=Subgroup.createGroup(str(md.__dict__[group].__dict__[field].__class__.__name__)+str(listindex)) - Listgroup.__setattr__('classtype',md.__dict__[group].__dict__[field].__getitem__(listindex).__class__.__name__) - try: - subfields=dict.keys(md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__) - except AttributeError: - subfields=dict.keys(md.__dict__[group].__dict__[field].__getitem__(listindex)) - for subfield in subfields: - if subfield!='outlog': - try: - Var=md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__[subfield] - except AttributeError: - Var=md.__dict__[group].__dict__[field].__getitem__(listindex)[subfield] - DimDict=CreateVar(NCData,Var,subfield,Listgroup,DimDict,md.__dict__[group],field,listindex) - - #No subgroup, we directly treat the variable - elif type(md.__dict__[group].__dict__[field]) in typelist or field=='bamg': - NCgroup.__setattr__('classtype', md.__dict__[group].__class__.__name__) - Var=md.__dict__[group].__dict__[field] - DimDict=CreateVar(NCData,Var,field,NCgroup,DimDict) - elif md.__dict__[group].__dict__[field] is None: - print(( 'field md.{}.{} is None'.format(group,field))) - #do nothing - #if it is a masked array - elif type(md.__dict__[group].__dict__[field]) is np.ma.core.MaskedArray: - NCgroup.__setattr__('classtype', md.__dict__[group].__class__.__name__) - Var=md.__dict__[group].__dict__[field].data - DimDict=CreateVar(NCData,Var,field,NCgroup,DimDict) - else: - NCgroup.__setattr__('classtype', str(group)) - Subgroup=NCgroup.createGroup(str(field)) - Subgroup.__setattr__('classtype',md.__dict__[group].__class__.__name__) - subfields=dict.keys(md.__dict__[group].__dict__[field].__dict__) - - for subfield in subfields: - if str(subfield)!='outlog': - Var=md.__dict__[group].__dict__[field].__dict__[subfield] - DimDict=CreateVar(NCData,Var,subfield,Subgroup,DimDict) - - NCData.close() - -#============================================================================ -#Define the variables -def CreateVar(NCData,var,field,Group,DimDict,*step_args): - #grab type - try: - val_type=str(var.dtype) - except AttributeError: - val_type=type(var) - #grab dimension - try: - val_shape=dict.keys(var) - except TypeError: - val_shape=np.shape(var) - - TypeDict = {float:'f8', - 'float64':'f8', - np.float64:'f8', - int:'i8', - 'int64':'i8', - np.int64:'i8', - str:str, - dict:str} - - val_dim=np.shape(val_shape)[0] - - #Now define and fill up variable - #treating scalar string or bool as atribute - if val_type in [str,str,bool]: - Group.__setattr__(str(field).swapcase(), str(var)) - #treating list as string table - elif val_type==list: - dimensions,DimDict=GetDim(NCData,var,val_shape,DimDict,val_dim) - #try to get the type from the first element - try: - nctype=TypeDict[type(var[0])] - except IndexError: - nctype=str #most probably an empty list take str for that - ncvar = Group.createVariable(str(field),nctype,dimensions,zlib=True) - if val_shape==0: - ncvar= [] - else: - for elt in range(0,val_shape[0]): - ncvar[elt] = var[elt] - #treating bool tables as string tables - elif val_type=='bool': - dimensions,DimDict=GetDim(NCData,var,val_shape,DimDict,val_dim) - ncvar = Group.createVariable(str(field),str,dimensions,zlib=True) - for elt in range(0,val_shape[0]): - ncvar[elt] = str(var[elt]) - #treating dictionaries as tables of strings - elif val_type==collections.OrderedDict or val_type==dict: - dimensions,DimDict=GetDim(NCData,var,val_shape,DimDict,val_dim) - ncvar = Group.createVariable(str(field),str,dimensions,zlib=True) - for elt in range(0,val_dim): - ncvar[elt,0]=dict.keys(var)[elt] - ncvar[elt,1]=str(dict.values(var)[elt]) #converting to str to avoid potential problems - #Now dealing with numeric variables - elif val_type in [float,'float64',np.float64,int,'int64']: - dimensions,DimDict=GetDim(NCData,var,val_shape,DimDict,val_dim) - ncvar = Group.createVariable(str(field),TypeDict[val_type],dimensions,zlib=True) - try: - nan_val=np.isnan(var) - if nan_val.all(): - ncvar [:] = 'NaN' - else: - ncvar[:] = var - except TypeError: #type does not accept nan, get vallue of the variable - ncvar[:] = var - else: - print(('WARNING type "{}" is unknown for "{}.{}"'.format(val_type,Group.name,field))) - return DimDict - -#============================================================================ -#retriev the dimension tuple from a dictionnary -def GetDim(NCData,var,shape,DimDict,i): - output=[] - #grab dimension - for dim in range(0,i): #loop on the dimensions - if type(shape[0])==int: - try: - output=output+[str(DimDict[shape[dim]])] #test if the dimension allready exist - except KeyError: #if not create it - if (shape[dim])>0: - index=len(DimDict)+1 - NewDim=NCData.createDimension('DimNum'+str(index),(shape[dim])) - DimDict[len(NewDim)]='DimNum'+str(index) - output=output+[str(DimDict[shape[dim]])] - elif type(shape[0])==str or type(shape[0])==str:#dealling with a dictionnary - try: - #dimension5 is 2 to treat with dict - output=[str(DimDict[np.shape(shape)[0]])]+[DimDict[2]] - except KeyError: - index=len(DimDict)+1 - NewDim=NCData.createDimension('DimNum'+str(index),np.shape(shape)[0]) - DimDict[len(NewDim)]='DimNum'+str(index) - output=[str(DimDict[np.shape(dict.keys(var))[0]])]+[DimDict[2]] - break - return tuple(output), DimDict + +def export_netCDF(md, filename): + if path.exists(filename): + print('File {} allready exist'.format(filename)) + newname = eval(input('Give a new name or "delete" to replace: ')) + if newname == 'delete': + remove(filename) + else: + print(('New file name is {}'.format(newname))) + filename = newname + NCData = Dataset(filename, 'w', format='NETCDF4') + NCData.description = 'Results for run' + md.miscellaneous.name + NCData.history = 'Created ' + time.ctime(time.time()) + # define netCDF dimensions + try: + StepNum = np.shape(dict.values(md.results.__dict__))[1] + except IndexError: + StepNum = 1 + Dimension1 = NCData.createDimension('DimNum1', StepNum) # time is first + DimDict = {len(Dimension1): 'DimNum1'} + dimindex = 1 + dimlist = [2, md.mesh.numberofelements, md.mesh.numberofvertices, np.shape(md.mesh.elements)[1]] + for i in range(0, 4): + if dimlist[i] not in list(DimDict.keys()): + dimindex += 1 + NewDim = NCData.createDimension('DimNum'+str(dimindex), dimlist[i]) + DimDict[len(NewDim)] = 'DimNum' + str(dimindex) + typelist = [bool, str, str, int, float, complex, + collections.OrderedDict, + np.int64, np.ndarray, np.float64] + groups = dict.keys(md.__dict__) + # get all model classes and create respective groups + for group in groups: + NCgroup = NCData.createGroup(str(group)) + # In each group gather the fields of the class + fields = dict.keys(md.__dict__[group].__dict__) + # looping on fields + for field in fields: + # Special treatment for list fields + if type(md.__dict__[group].__dict__[field]) == list: + StdList = False + if len(md.__dict__[group].__dict__[field]) == 0: + StdList = True + else: + StdList = type(md.__dict__[group].__dict__[field][0]) in typelist + NCgroup.__setattr__('classtype', md.__dict__[group].__class__.__name__) + if StdList: # this is a standard or empty list just proceed + Var = md.__dict__[group].__dict__[field] + DimDict = CreateVar(NCData, Var, field, NCgroup, DimDict) + else: # this is a list of fields, specific treatment needed + Listsize = len(md.__dict__[group].__dict__[field]) + Subgroup = NCgroup.createGroup(str(field)) + Subgroup.__setattr__('classtype', md.__dict__[group].__dict__[field].__class__.__name__) + for listindex in range(0, Listsize): + try: + Listgroup = Subgroup.createGroup(str(md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__['name'])) + except KeyError: + for naming in ['step']: + Listgroup = Subgroup.createGroup(str(md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__[naming])) + except AttributeError: + Listgroup=Subgroup.createGroup(str(md.__dict__[group].__dict__[field].__class__.__name__)+str(listindex)) + Listgroup.__setattr__('classtype',md.__dict__[group].__dict__[field].__getitem__(listindex).__class__.__name__) + try: + subfields=dict.keys(md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__) + except AttributeError: + subfields=dict.keys(md.__dict__[group].__dict__[field].__getitem__(listindex)) + for subfield in subfields: + if subfield!='outlog': + try: + Var=md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__[subfield] + except AttributeError: + Var = md.__dict__[group].__dict__[field].__getitem__(listindex)[subfield] + DimDict = CreateVar(NCData,Var,subfield,Listgroup,DimDict,md.__dict__[group],field,listindex) + # No subgroup, we directly treat the variable + elif type(md.__dict__[group].__dict__[field]) in typelist or field == 'bamg': + NCgroup.__setattr__('classtype', md.__dict__[group].__class__.__name__) + Var = md.__dict__[group].__dict__[field] + DimDict = CreateVar(NCData, Var, field, NCgroup, DimDict) + elif md.__dict__[group].__dict__[field] is None: + print('field md.{}.{} is None'.format(group, field)) + # do nothing + # if it is a masked array + elif type(md.__dict__[group].__dict__[field]) is np.ma.core.MaskedArray: + NCgroup.__setattr__('classtype', md.__dict__[group].__class__.__name__) + Var = md.__dict__[group].__dict__[field].data + DimDict = CreateVar(NCData,Var,field,NCgroup,DimDict) + else: + NCgroup.__setattr__('classtype', str(group)) + Subgroup = NCgroup.createGroup(str(field)) + Subgroup.__setattr__('classtype', md.__dict__[group].__class__.__name__) + subfields = dict.keys(md.__dict__[group].__dict__[field].__dict__) + for subfield in subfields: + if str(subfield) != 'outlog': + Var = md.__dict__[group].__dict__[field].__dict__[subfield] + DimDict = CreateVar(NCData, Var, subfield, Subgroup, DimDict) + NCData.close() + +# ============================================================================ +# Define the variables + + +def CreateVar(NCData, var, field, Group, DimDict, *step_args): + # grab type + try: + val_type = str(var.dtype) + except AttributeError: + val_type = type(var) + # grab dimension + if val_type in [collections.OrderedDict, dict]: + val_shape = len(var) + val_dim = 2 + else: + val_shape = np.shape(var) + val_dim = np.shape(val_shape)[0] + TypeDict = {float: 'f8', + 'float64': 'f8', + np.float64: 'f8', + int: 'i8', + 'int64': 'i8', + np.int64: 'i8', + str: str, + dict: str} + # Now define and fill up variable + # treating scalar string or bool as atribute + if val_type in [str, bool]: + Group.__setattr__(str(field).swapcase(), str(var)) + # treating list as string table + elif val_type == list: + dimensions, DimDict = GetDim(NCData, val_shape, val_type, DimDict, val_dim) + # try to get the type from the first element + try: + nctype = TypeDict[type(var[0])] + except IndexError: + nctype = str # most probably an empty list take str for that + ncvar = Group.createVariable(str(field), nctype, dimensions, zlib=True) + if val_shape == 0: + ncvar = [] + else: + for elt in range(0, val_shape[0]): + ncvar[elt] = var[elt] + # treating bool tables as string tables + elif val_type == 'bool': + dimensions, DimDict = GetDim(NCData, val_shape, val_type, DimDict, val_dim) + ncvar = Group.createVariable(str(field), str, dimensions, zlib=True) + for elt in range(0, val_shape[0]): + ncvar[elt] = str(var[elt]) + # treating dictionaries as tables of strings + elif val_type in [collections.OrderedDict, dict]: + dimensions, DimDict = GetDim(NCData, val_shape, val_type, DimDict, val_dim) + ncvar = Group.createVariable(str(field), str, dimensions, zlib=True) + for elt, key in enumerate(dict.keys(var)): + ncvar[elt, 0] = key + ncvar[elt, 1] = str(var[key]) # converting to str to avoid potential problems + # Now dealing with numeric variables + elif val_type in [float, 'float64', np.float64, int, 'int64']: + dimensions, DimDict = GetDim(NCData, val_shape, val_type, DimDict, val_dim) + ncvar = Group.createVariable(str(field), TypeDict[val_type], dimensions, zlib=True) + try: + nan_val = np.isnan(var) + if nan_val.all(): + ncvar[:] = 'NaN' + else: + ncvar[:] = var + except TypeError: # type does not accept nan, get vallue of the variable + ncvar[:] = var + else: + print(('WARNING type "{}" is unknown for "{}.{}"'.format(val_type, Group.name, field))) + return DimDict + +# ============================================================================ +# retriev the dimension tuple from a dictionnary + + +def GetDim(NCData, val_shape, val_type, DimDict, val_dim): + output = [] + if val_type in [collections.OrderedDict, dict]: # dealling with a dictionnary + try: + output = [str(DimDict[val_shape])] # first try to get the coresponding dimension if ti exists + output = output + [DimDict[2]] # dimension5 is 2 to treat with dict + except KeyError: + index = len(DimDict) + 1 # if the dimension does not exist, increment naming + NewDim = NCData.createDimension('DimNum'+str(index), val_shape) # create dimension + DimDict[len(NewDim)] = 'DimNum' + str(index) # and update the dimension dictionary + output = [str(DimDict[val_shape])] + [DimDict[2]] # now proceed with the shape of the value + else: + # loop on dimensions + for dim in range(0, val_dim): # loop on the dimensions + try: + output = output + [str(DimDict[val_shape[dim]])] # test if the dimension allready exist + except KeyError: # if not create it + if (val_shape[dim]) > 0: + index = len(DimDict) + 1 + NewDim = NCData.createDimension('DimNum' + str(index), (val_shape[dim])) + DimDict[len(NewDim)] = 'DimNum' + str(index) + output = output + [str(DimDict[val_shape[dim]])] + return tuple(output), DimDict Index: /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/read_netCDF.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/read_netCDF.py (revision 23771) +++ /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/read_netCDF.py (revision 23772) @@ -1,25 +1,20 @@ from netCDF4 import Dataset -import time -import collections -from os import path, remove def netCDFRead(filename): - - def walktree(data): - keys = list(data.groups.keys()) - yield keys - for key in keys: - for children in walktree(data.groups[str(key)]): - yield children - - if path.exists(filename): - print(('Opening {} for reading '.format(filename))) - NCData=Dataset(filename, 'r') - class_dict={} - - for children in walktree(NCData): - for child in children: - class_dict[str(child)]=str(getattr(NCData.groups[str(child)],'classtype')+'()') + def walktree(data): + keys = list(data.groups.keys()) + yield keys + for key in keys: + for children in walktree(data.groups[str(key)]): + yield children - print(class_dict) - + if path.exists(filename): + print(('Opening {} for reading '.format(filename))) + NCData=Dataset(filename, 'r') + class_dict={} + + for children in walktree(NCData): + for child in children: + class_dict[str(child)]=str(getattr(NCData.groups[str(child)],'classtype')+'()') + + print(class_dict) Index: /issm/trunk-jpl/src/m/contrib/defleurian/paraview/enveloppeVTK.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/paraview/enveloppeVTK.py (revision 23771) +++ /issm/trunk-jpl/src/m/contrib/defleurian/paraview/enveloppeVTK.py (revision 23772) @@ -28,5 +28,5 @@ if os.path.exists(filename): print(('File {} allready exist'.format(filename))) - newname=eval(input('Give a new name or "delete" to replace: ')) + newname=input('Give a new name or "delete" to replace: ') if newname=='delete': filelist = glob.glob(filename+'/*') Index: /issm/trunk-jpl/src/m/io/loadmodel.py =================================================================== --- /issm/trunk-jpl/src/m/io/loadmodel.py (revision 23771) +++ /issm/trunk-jpl/src/m/io/loadmodel.py (revision 23772) @@ -2,39 +2,40 @@ #hack to keep python 2 compatipility try: - #py3 import - from dbm.ndbm import whichdb + #py3 import + from dbm import whichdb except ImportError: - #py2 import - from whichdb import whichdb + #py2 import + from whichdb import whichdb from netCDF4 import Dataset + def loadmodel(path): - """ - LOADMODEL - load a model using built-in load module + """ + LOADMODEL - load a model using built-in load module - check that model prototype has not changed. if so, adapt to new model prototype. + check that model prototype has not changed. if so, adapt to new model prototype. - Usage: - md=loadmodel(path) - """ + Usage: + md=loadmodel(path) + """ - #check existence of database (independent of file extension!) - if whichdb(path): - #do nothing - pass - else: - try: - NCFile=Dataset(path,mode='r') - NCFile.close() - pass - except RuntimeError: - raise IOError("loadmodel error message: file '%s' does not exist" % path) - # try: - #recover model on file and name it md - struc=loadvars(path) - name=[key for key in list(struc.keys())] - if len(name)>1: - raise IOError("loadmodel error message: file '%s' contains several variables. Only one model should be present." % path) + #check existence of database (independent of file extension!) + if whichdb(path): + #do nothing + pass + else: + try: + NCFile = Dataset(path, mode='r') + NCFile.close() + pass + except RuntimeError: + raise IOError("loadmodel error message: file '%s' does not exist" % path) + # try: + #recover model on file and name it md + struc = loadvars(path) + name = [key for key in list(struc.keys())] + if len(name) > 1: + raise IOError("loadmodel error message: file '%s' contains several variables. Only one model should be present." % path) - md=struc[name[0]] - return md + md = struc[name[0]] + return md Index: /issm/trunk-jpl/src/m/io/loadvars.py =================================================================== --- /issm/trunk-jpl/src/m/io/loadvars.py (revision 23771) +++ /issm/trunk-jpl/src/m/io/loadvars.py (revision 23772) @@ -1,254 +1,258 @@ import shelve -import os.path -import numpy as np +import numpy as np from netCDF4 import Dataset -from netCDF4 import chartostring from re import findall -from os import path from collections import OrderedDict from model import * -#hack to keep python 2 compatipility +#hack to keep python 2 compatibility try: - #py3 import - from dbm.ndbm import whichdb + #py3 import + from dbm import whichdb except ImportError: - #py2 import - from whichdb import whichdb + #py2 import + from whichdb import whichdb + def loadvars(*args): - """ - LOADVARS - function to load variables to a file. - - This function loads one or more variables from a file. The names of the variables - must be supplied. If more than one variable is specified, it may be done with - a list of names or a dictionary of name as keys. The output type will correspond - to the input type. All the variables in the file may be loaded by specifying only - the file name. - - Usage: - a=loadvars('shelve.dat','a') - [a,b]=loadvars('shelve.dat',['a','b']) - nvdict=loadvars('shelve.dat',{'a':None,'b':None}) - nvdict=loadvars('shelve.dat') - - """ - - filename='' - nvdict={} - - if len(args) >= 1 and isinstance(args[0],str): - filename=args[0] - if not filename: - filename='/tmp/shelve.dat' - - else: - raise TypeError("Missing file name.") - - if len(args) >= 2 and isinstance(args[1],str): # (filename,name) - for name in args[1:]: - nvdict[name]=None - - elif len(args) == 2 and isinstance(args[1],list): # (filename,[names]) - for name in args[1]: - nvdict[name]=None - - elif len(args) == 2 and isinstance(args[1],dict): # (filename,{names:values}) - nvdict=args[1] - - elif len(args) == 1: # (filename) - pass - - else: - raise TypeError("Unrecognized input arguments.") - - if whichdb(filename): - print(("Loading variables from file '%s'." % filename)) - - my_shelf = shelve.open(filename,'r') # 'r' for read-only - if nvdict: - for name in list(nvdict.keys()): - try: - nvdict[name] = my_shelf[name] - print(("Variable '%s' loaded." % name)) - except KeyError: - value = None - print(("Variable '%s' not found." % name)) - - else: - for name in list(my_shelf.keys()): - nvdict[name] = my_shelf[name] - print(("Variable '%s' loaded." % name)) - - my_shelf.close() - - else: - try: - NCFile=Dataset(filename,mode='r') - NCFile.close() - except RuntimeError: - raise IOError("File '%s' not found." % filename) - - classtype,classtree=netCDFread(filename) - nvdict['md']=model() - NCFile=Dataset(filename,mode='r') - for mod in dict.keys(classtype): - if np.size(classtree[mod])>1: - curclass=NCFile.groups[classtree[mod][0]].groups[classtree[mod][1]] - if classtype[mod][0]=='list': - keylist=[key for key in curclass.groups] - try: - steplist=[int(key) for key in curclass.groups] - except ValueError: - steplist=[int(findall(r'\d+',key)[0]) for key in keylist] - indexlist=[index*(len(curclass.groups)-1)/max(1,max(steplist)) for index in steplist] - listtype=curclass.groups[keylist[0]].classtype - if listtype=='dict': - nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]] = [OrderedDict() for i in range(max(1,len(curclass.groups)))] - else: - nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]] = [getattr(__import__(listtype),listtype)() for i in range(max(1,len(curclass.groups)))] - Tree=nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]][:] - else: - nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]] = getattr(classtype[mod][1],classtype[mod][0])() - Tree=nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]] - else: - curclass=NCFile.groups[classtree[mod][0]] - nvdict['md'].__dict__[mod] = getattr(classtype[mod][1],classtype[mod][0])() - Tree=nvdict['md'].__dict__[classtree[mod][0]] - #treating groups that are lists - for i in range(0,max(1,len(curclass.groups))): - if len(curclass.groups)>0: - listclass=curclass.groups[keylist[i]] - else: - listclass=curclass - for var in listclass.variables: - if var not in ['errlog','outlog']: - varval=listclass.variables[str(var)] - vardim=varval.ndim - try: - val_type=str(varval.dtype) - except AttributeError: - val_type=type(varval) - if vardim==0: - if type(Tree)==list: - t=indexlist[i] - if listtype=='dict': - Tree[t][str(var)]=varval[0] - else: - Tree[t].__dict__[str(var)]=varval[0] - else: - if str(varval[0])=='': - Tree.__dict__[str(var)]=[] - elif varval[0]=='True': - Tree.__dict__[str(var)]=True - elif varval[0]=='False': - Tree.__dict__[str(var)]=False - else: - Tree.__dict__[str(var)]=varval[0] - - elif vardim==1: - if varval.dtype==str: - if varval.shape[0]==1: - Tree.__dict__[str(var)]=[str(varval[0]),] - elif 'True' in varval[:] or 'False' in varval[:]: - Tree.__dict__[str(var)]=np.asarray([V=='True' for V in varval[:]],dtype=bool) - else: - Tree.__dict__[str(var)]=[str(vallue) for vallue in varval[:]] - else: - if type(Tree)==list: - t=indexlist[i] - if listtype=='dict': - Tree[t][str(var)]=varval[:] - else: - Tree[t].__dict__[str(var)]=varval[:] - else: - try: - #some thing specifically require a list - mdtype=type(Tree.__dict__[str(var)]) - except KeyError: - mdtype=float - if mdtype==list: - Tree.__dict__[str(var)]=[mdval for mdval in varval[:]] - else: - Tree.__dict__[str(var)]=varval[:] - elif vardim==2: - #dealling with dict - if varval.dtype==str: #that is for toolkits wich needs to be ordered - if any(varval[:,0]=='toolkit'): #toolkit definition have to be first - Tree.__dict__[str(var)]=OrderedDict([('toolkit', str(varval[np.where(varval[:,0]=='toolkit')[0][0],1]))]) - - strings1=[str(arg[0]) for arg in varval if arg[0]!='toolkits'] - strings2=[str(arg[1]) for arg in varval if arg[0]!='toolkits'] - Tree.__dict__[str(var)].update(list(zip(strings1, strings2))) - else: - if type(Tree)==list: - #t=int(keylist[i][-1])-1 - t=indexlist[i] - if listtype=='dict': - Tree[t][str(var)]=varval[:,:] - else: - Tree[t].__dict__[str(var)]=varval[:,:] - else: - Tree.__dict__[str(var)]=varval[:,:] - elif vardim==3: - if type(Tree)==list: - t=indexlist[i] - if listtype=='dict': - Tree[t][str(var)]=varval[:,:,:] - else: - Tree[t].__dict__[str(var)]=varval[:,:,:] - else: - Tree.__dict__[str(var)]=varval[:,:,:] - else: - print('table dimension greater than 3 not implemented yet') - for attr in listclass.ncattrs(): - if attr!='classtype': #classtype is for treatment, don't get it back - if type(Tree)==list: - t=indexlist[i] - if listtype=='dict': - Tree[t][str(attr).swapcase()]=str(listclass.getncattr(attr)) - else: - Tree[t].__dict__[str(attr).swapcase()]=str(listclass.getncattr(attr)) - else: - Tree.__dict__[str(attr).swapcase()]=str(listclass.getncattr(attr)) - if listclass.getncattr(attr)=='True': - Tree.__dict__[str(attr).swapcase()]=True - elif listclass.getncattr(attr)=='False': - Tree.__dict__[str(attr).swapcase()]=False - NCFile.close() - if len(args) >= 2 and isinstance(args[1],str): # (value) - value=[nvdict[name] for name in args[1:]] - return value - - elif len(args) == 2 and isinstance(args[1],list): # ([values]) - value=[nvdict[name] for name in args[1]] - return value - - elif (len(args) == 2 and isinstance(args[1],dict)) or (len(args) == 1): # ({names:values}) - return nvdict + """ + LOADVARS - function to load variables to a file. + + This function loads one or more variables from a file. The names of the variables + must be supplied. If more than one variable is specified, it may be done with + a list of names or a dictionary of name as keys. The output type will correspond + to the input type. All the variables in the file may be loaded by specifying only + the file name. + + Usage: + a=loadvars('shelve.dat','a') + [a,b]=loadvars('shelve.dat',['a','b']) + nvdict=loadvars('shelve.dat',{'a':None,'b':None}) + nvdict=loadvars('shelve.dat') + + """ + + filename = '' + nvdict = {} + + if len(args) >= 1 and isinstance(args[0], str): + filename = args[0] + if not filename: + filename = '/tmp/shelve.dat' + + else: + raise TypeError("Missing file name.") + + if len(args) >= 2 and isinstance(args[1], str): # (filename,name) + for name in args[1:]: + nvdict[name] = None + + elif len(args) == 2 and isinstance(args[1], list): # (filename,[names]) + for name in args[1]: + nvdict[name] = None + + elif len(args) == 2 and isinstance(args[1], dict): # (filename,{names:values}) + nvdict = args[1] + + elif len(args) == 1: # (filename) + pass + + else: + raise TypeError("Unrecognized input arguments.") + + if whichdb(filename): + print("Loading variables from file {}.".format(filename)) + my_shelf = shelve.open(filename, 'r') # 'r' for read-only + if nvdict: + for name in list(nvdict.keys()): + try: + nvdict[name] = my_shelf[name] + print(("Variable '%s' loaded." % name)) + except KeyError: + value = None + print("Variable '{}' not found.".format(name)) + + else: + for name in list(my_shelf.keys()): + nvdict[name] = my_shelf[name] + print(("Variable '%s' loaded." % name)) + my_shelf.close() + + else: + try: + NCFile = Dataset(filename, mode='r') + NCFile.close() + except RuntimeError: + raise IOError("File '{}' not found.".format(filename)) + + classtype, classtree = netCDFread(filename) + nvdict['md'] = model() + NCFile = Dataset(filename, mode='r') + for mod in dict.keys(classtype): + #print('-Now treating classtype {}'.format(mod)) + if np.size(classtree[mod]) > 1: + curclass = NCFile.groups[classtree[mod][0]].groups[classtree[mod][1]] + if classtype[mod][0] == 'list': + keylist = [key for key in curclass.groups] + try: + steplist = [int(key) for key in curclass.groups] + except ValueError: + steplist = [int(findall(r'\d+', key)[0]) for key in keylist] + indexlist = [index * (len(curclass.groups) - 1) / max(1, max(steplist)) for index in steplist] + listtype = curclass.groups[keylist[0]].classtype + if listtype == 'dict': + nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]] = [OrderedDict() for i in range(max(1, len(curclass.groups)))] + else: + nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]] = [getattr(__import__(listtype), listtype)() for i in range(max(1, len(curclass.groups)))] + Tree = nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]][:] + else: + nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]] = getattr(classtype[mod][1], classtype[mod][0])() + Tree = nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]] + else: + curclass = NCFile.groups[classtree[mod][0]] + nvdict['md'].__dict__[mod] = getattr(classtype[mod][1], classtype[mod][0])() + Tree = nvdict['md'].__dict__[classtree[mod][0]] + #treating groups that are lists + for i in range(0, max(1, len(curclass.groups))): + if len(curclass.groups) > 0: + listclass = curclass.groups[keylist[i]] + else: + listclass = curclass + for var in listclass.variables: + #print("treating var {}".format(var)) + if var not in ['errlog', 'outlog']: + varval = listclass.variables[str(var)] + vardim = varval.ndim + if vardim == 0: + if type(Tree) == list: + t = int(indexlist[i]) + if listtype == 'dict': + Tree[t][str(var)] = varval[0] + + else: + Tree[t].__dict__[str(var)] = varval[0] + + else: + if str(varval[0]) == '': #no value + Tree.__dict__[str(var)] = [] + + elif varval[0] == 'True': #treatin bool + Tree.__dict__[str(var)] = True + + elif varval[0] == 'False': #treatin bool + Tree.__dict__[str(var)] = False + + else: + Tree.__dict__[str(var)] = varval[0].item() + + elif vardim == 1: + if varval.dtype == str: + if varval.shape[0] == 1: + Tree.__dict__[str(var)] = [str(varval[0]), ] + + elif 'True' in varval[:] or 'False' in varval[:]: + Tree.__dict__[str(var)] = np.asarray([V == 'True' for V in varval[:]], dtype=bool) + + else: + Tree.__dict__[str(var)] = [str(vallue) for vallue in varval[:]] + + else: + if type(Tree) == list: + t = int(indexlist[i]) + if listtype == 'dict': + Tree[t][str(var)] = varval[:] + + else: + Tree[t].__dict__[str(var)] = varval[:] + + else: + try: + #some thing specifically require a list + mdtype = type(Tree.__dict__[str(var)]) + except KeyError: + mdtype = float + if mdtype == list: + Tree.__dict__[str(var)] = [mdval for mdval in varval[:]] + + else: + Tree.__dict__[str(var)] = varval[:].data + + elif vardim == 2: + #dealling with dict + if varval.dtype == str: #that is for toolkits wich needs to be ordered + if any(varval[:, 0] == 'toolkit'): #toolkit definition have to be first + Tree.__dict__[str(var)] = OrderedDict([('toolkit', str(varval[np.where(varval[:, 0] == 'toolkit')[0][0], 1]))]) + strings1 = [str(arg[0]) for arg in varval if arg[0] != 'toolkits'] + strings2 = [str(arg[1]) for arg in varval if arg[0] != 'toolkits'] + Tree.__dict__[str(var)].update(list(zip(strings1, strings2))) + else: + if type(Tree) == list: + t = int(indexlist[i]) + if listtype == 'dict': + Tree[t][str(var)] = varval[:, :] + else: + Tree[t].__dict__[str(var)] = varval[:, :] + else: + Tree.__dict__[str(var)] = varval[:, :].data + elif vardim == 3: + if type(Tree) == list: + t = int(indexlist[i]) + if listtype == 'dict': + Tree[t][str(var)] = varval[:, :, :] + else: + Tree[t].__dict__[str(var)] = varval[:, :, :] + else: + Tree.__dict__[str(var)] = varval[:, :, :].data + else: + print('table dimension greater than 3 not implemented yet') + for attr in listclass.ncattrs(): + if attr != 'classtype': #classtype is for treatment, don't get it back + if type(Tree) == list: + t = int(indexlist[i]) + if listtype == 'dict': + Tree[t][str(attr).swapcase()] = str(listclass.getncattr(attr)) + else: + Tree[t].__dict__[str(attr).swapcase()] = str(listclass.getncattr(attr)) + else: + Tree.__dict__[str(attr).swapcase()] = str(listclass.getncattr(attr)) + if listclass.getncattr(attr) == 'True': + Tree.__dict__[str(attr).swapcase()] = True + elif listclass.getncattr(attr) == 'False': + Tree.__dict__[str(attr).swapcase()] = False + NCFile.close() + if len(args) >= 2 and isinstance(args[1], str): # (value) + value = [nvdict[name] for name in args[1:]] + return value + + elif len(args) == 2 and isinstance(args[1], list): # ([values]) + value = [nvdict[name] for name in args[1]] + return value + + elif (len(args) == 2 and isinstance(args[1], dict)) or (len(args) == 1): # ({names:values}) + return nvdict def netCDFread(filename): - print(('Opening {} for reading '.format(filename))) - NCData=Dataset(filename, 'r') - class_dict={} - class_tree={} - - for group in NCData.groups: - if len(NCData.groups[group].groups)>0: - for subgroup in NCData.groups[group].groups: - classe=str(group)+'.'+str(subgroup) - class_dict[classe]=[str(getattr(NCData.groups[group].groups[subgroup],'classtype')),] - if class_dict[classe][0] not in ['dict','list','cell']: - class_dict[classe].append(__import__(class_dict[classe][0])) - class_tree[classe]=[group,subgroup] - else: - classe=str(group) - try: - class_dict[classe]=[str(getattr(NCData.groups[group],'classtype')),] - if class_dict[classe][0] not in ['dict','list','cell']: - class_dict[classe].append(__import__(class_dict[classe][0])) - class_tree[classe]=[group,] - except AttributeError: - print(('group {} is empty'.format(group))) - NCData.close() - return class_dict,class_tree + print(('Opening {} for reading '.format(filename))) + NCData = Dataset(filename, 'r') + class_dict = {} + class_tree = {} + + for group in NCData.groups: + if len(NCData.groups[group].groups) > 0: + for subgroup in NCData.groups[group].groups: + classe = str(group) + '.' + str(subgroup) + class_dict[classe] = [str(getattr(NCData.groups[group].groups[subgroup], 'classtype')), ] + if class_dict[classe][0] not in ['dict', 'list', 'cell']: + class_dict[classe].append(__import__(class_dict[classe][0])) + class_tree[classe] = [group, subgroup] + else: + classe = str(group) + try: + class_dict[classe] = [str(getattr(NCData.groups[group], 'classtype')), ] + if class_dict[classe][0] not in ['dict', 'list', 'cell']: + class_dict[classe].append(__import__(class_dict[classe][0])) + class_tree[classe] = [group, ] + except AttributeError: + print(('group {} is empty'.format(group))) + NCData.close() + return class_dict, class_tree Index: /issm/trunk-jpl/src/m/os/issmscpin.py =================================================================== --- /issm/trunk-jpl/src/m/os/issmscpin.py (revision 23771) +++ /issm/trunk-jpl/src/m/os/issmscpin.py (revision 23772) @@ -54,9 +54,8 @@ #string to copy multiple files using scp: string="'{"+','.join([str(x) for x in packages])+"}'" - if port: - subprocess.call('scp -P %d %s@localhost:%s %s/. ' % (port,login,os.path.join(path,string),os.getcwd()),shell=True) + subprocess.call('scp -P {} {}@localhost:{} {}/. '.format(port,login,os.path.join(path,string),os.getcwd()),shell=True) else: - subprocess.call('scp -T %s@%s:%s %s/.' % (login,host,os.path.join(path,string),os.getcwd()),shell=True) + subprocess.call('scp -T {}@{}:{} {}/.'.format(login,host,os.path.join(path,string),os.getcwd()),shell=True) #check scp worked for package in packages: Index: /issm/trunk-jpl/src/m/plot/plot_overlay.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_overlay.py (revision 23771) +++ /issm/trunk-jpl/src/m/plot/plot_overlay.py (revision 23772) @@ -67,5 +67,5 @@ # normalize array arr=arr/np.float(np.max(arr.ravel())) - arr=1.-arr # somehow the values got flipped + arr=1.-arr # somehow the values got flipped if options.getfieldvalue('overlayhist',0)==1: Index: /issm/trunk-jpl/src/m/solve/WriteData.py =================================================================== --- /issm/trunk-jpl/src/m/solve/WriteData.py (revision 23771) +++ /issm/trunk-jpl/src/m/solve/WriteData.py (revision 23772) @@ -1,375 +1,368 @@ +from struct import pack, error import numpy as np -from struct import pack,unpack import pairoptions -def WriteData(fid,prefix,*args): - """ - WRITEDATA - write model field in binary file - - Usage: - WriteData(fid,varargin) - """ - - #process options - options=pairoptions.pairoptions(*args) - - #Get data properties - if options.exist('object'): - #This is an object field, construct enum and data - obj = options.getfieldvalue('object') - fieldname = options.getfieldvalue('fieldname') - classname = options.getfieldvalue('class',str(type(obj)).rsplit('.')[-1].split("'")[0]) - name = options.getfieldvalue('name',prefix+'.'+fieldname); - if options.exist('data'): - data = options.getfieldvalue('data') - else: - data = getattr(obj,fieldname) - else: - #No processing required - data = options.getfieldvalue('data') - name = options.getfieldvalue('name') - - datatype = options.getfieldvalue('format') - mattype = options.getfieldvalue('mattype',0) #only required for matrices - timeserieslength = options.getfieldvalue('timeserieslength',-1) - - #Process sparse matrices -# if issparse(data), -# data=full(data); -# end - - #Scale data if necesarry - if options.exist('scale'): - data=np.array(data) - scale = options.getfieldvalue('scale') - if np.size(data) > 1 and np.ndim(data) > 1 and np.size(data,0)==timeserieslength: - data[0:-1,:] = scale*data[0:-1,:] - else: - data = scale*data - if np.size(data) > 1 and np.size(data,0)==timeserieslength: - yts = options.getfieldvalue('yts') - if np.ndim(data) > 1: - data[-1,:] = yts*data[-1,:] - else: - data[-1] = yts*data[-1] - - #Step 1: write the enum to identify this record uniquely - fid.write(pack('i',len(name))) - fid.write(pack('{}s'.format(len(name)),name.encode())) - - - #Step 2: write the data itself. - if datatype=='Boolean': # {{{ -# if len(data) !=1: -# raise ValueError('fi eld %s cannot be marshalled as it has more than one element!' % name[0]) - - #first write length of record - fid.write(pack('q',4+4)) #1 bool (disguised as an int)+code - - #write data code: - fid.write(pack('i',FormatToCode(datatype))) - - #now write integer - fid.write(pack('i',int(data))) #send an int, not easy to send a bool - # }}} - - elif datatype=='Integer': # {{{ -# if len(data) !=1: -# raise ValueError('field %s cannot be marshalled as it has more than one element!' % name[0]) - - #first write length of record - fid.write(pack('q',4+4)) #1 integer + code - - #write data code: - fid.write(pack('i',FormatToCode(datatype))) - - #now write integer - fid.write(pack('i',int(data))) #force an int, - # }}} - - elif datatype=='Double': # {{{ -# if len(data) !=1: -# raise ValueError('field %s cannot be marshalled as it has more than one element!' % name[0]) - - #first write length of record - fid.write(pack('q',8+4)) #1 double+code - - #write data code: - fid.write(pack('i',FormatToCode(datatype))) - - #now write double - fid.write(pack('d',data)) - # }}} - - elif datatype=='String': # {{{ - #first write length of record - fid.write(pack('q',len(data)+4+4)) #string + string size + code - - #write data code: - fid.write(pack('i',FormatToCode(datatype))) - - #now write string - fid.write(pack('i',len(data))) - fid.write(pack('{}s'.format(len(data)),data.encode())) - # }}} - - elif datatype in ['IntMat','BooleanMat']: # {{{ - - if isinstance(data,(int,bool)): - data=np.array([data]) - elif isinstance(data,(list,tuple)): - data=np.array(data).reshape(-1,) - if np.ndim(data) == 1: - if np.size(data): - data=data.reshape(np.size(data),) - else: - data=data.reshape(0,0) - - #Get size - s=data.shape - #if matrix = NaN, then do not write anything - if np.ndim(data)==2 and np.product(s)==1 and np.all(np.isnan(data)): - s=(0,0) - - #first write length of record - fid.write(pack('q',4+4+8*np.product(s)+4+4)) #2 integers (32 bits) + the double matrix + code + matrix type - - #write data code and matrix type: - fid.write(pack('i',FormatToCode(datatype))) - fid.write(pack('i',mattype)) - - #now write matrix - if np.ndim(data) == 1: - fid.write(pack('i',s[0])) - fid.write(pack('i',1)) - for i in range(s[0]): - fid.write(pack('d',float(data[i]))) #get to the "c" convention, hence the transpose - else: - fid.write(pack('i',s[0])) - fid.write(pack('i',s[1])) - for i in range(s[0]): - for j in range(s[1]): - fid.write(pack('d',float(data[i][j]))) #get to the "c" convention, hence the transpose - # }}} - - elif datatype=='DoubleMat': # {{{ - - if isinstance(data,(bool,int,float)): - data=np.array([data]) - elif isinstance(data,(list,tuple)): - data=np.array(data).reshape(-1,) - if np.ndim(data) == 1: - if np.size(data): - data=data.reshape(np.size(data),) - else: - data=data.reshape(0,0) - - #Get size - s=data.shape - #if matrix = NaN, then do not write anything - if np.ndim(data)==1 and np.product(s)==1 and np.all(np.isnan(data)): - s=(0,0) - - #first write length of record - recordlength=4+4+8*np.product(s)+4+4; #2 integers (32 bits) + the double matrix + code + matrix type - if recordlength > 4**31 : - raise ValueError('field {} cannot be marshalled because it is larger than 4^31 bytes!'.format(enum)) - - fid.write(pack('q',recordlength)) #2 integers (32 bits) + the double matrix + code + matrix type - - #write data code and matrix type: - fid.write(pack('i',FormatToCode(datatype))) - fid.write(pack('i',mattype)) - - #now write matrix - if np.ndim(data) == 1: - fid.write(pack('i',s[0])) - fid.write(pack('i',1)) - for i in range(s[0]): - fid.write(pack('d',float(data[i]))) #get to the "c" convention, hence the transpose - else: - fid.write(pack('i',s[0])) - fid.write(pack('i',s[1])) - for i in range(s[0]): - for j in range(s[1]): - fid.write(pack('d',float(data[i][j]))) #get to the "c" convention, hence the transpose - # }}} - - elif datatype=='CompressedMat': # {{{ - - if isinstance(data,(bool,int,float)): - data=np.array([data]) - elif isinstance(data,(list,tuple)): - data=np.array(data).reshape(-1,) - if np.ndim(data) == 1: - if np.size(data): - data=data.reshape(np.size(data),) - else: - data=data.reshape(0,0) - - #Get size - s=data.shape - if np.ndim(data) == 1: - n2=1 - else: - n2=s[1] - - #if matrix = NaN, then do not write anything - if np.ndim(data)==1 and np.product(s)==1 and np.all(np.isnan(data)): - s=(0,0) - n2=0 - - #first write length of record - recordlength=4+4+8+8+1*(s[0]-1)*n2+8*n2+4+4 #2 integers (32 bits) + the matrix + code + matrix type - if recordlength > 4**31 : - raise ValueError('field %s cannot be marshalled because it is larger than 4^31 bytes!' % enum) - - fid.write(pack('q',recordlength)) #2 integers (32 bits) + the matrix + code + matrix type - - #write data code and matrix type: - fid.write(pack('i',FormatToCode(datatype))) - fid.write(pack('i',mattype)) - - #Write offset and range - A = data[0:s[0]-1] - offsetA = A.min() - rangeA = A.max() - offsetA - - if rangeA == 0: - A = A*0 - else: - A = (A-offsetA)/rangeA*255. - - #now write matrix - if np.ndim(data) == 1: - fid.write(pack('i',s[0])) - fid.write(pack('i',1)) - fid.write(pack('d',float(offsetA))) - fid.write(pack('d',float(rangeA))) - for i in range(s[0]-1): - fid.write(pack('B',int(A[i]))) - - fid.write(pack('d',float(data[s[0]-1]))) #get to the "c" convention, hence the transpose - - elif np.product(s) > 0: - fid.write(pack('i',s[0])) - fid.write(pack('i',s[1])) - fid.write(pack('d',float(offsetA))) - fid.write(pack('d',float(rangeA))) - for i in range(s[0]-1): - for j in range(s[1]): - fid.write(pack('B',int(A[i][j]))) #get to the "c" convention, hence the transpose - - for j in range(s[1]): - fid.write(pack('d',float(data[s[0]-1][j]))) - - # }}} - - elif datatype=='MatArray': # {{{ - - #first get length of record - recordlength=4+4 #number of records + code - for matrix in data: - if isinstance(matrix,(bool,int,float)): - matrix=np.array([matrix]) - elif isinstance(matrix,(list,tuple)): - matrix=np.array(matrix).reshape(-1,) - if np.ndim(matrix) == 1: - if np.size(matrix): - matrix=matrix.reshape(np.size(matrix),) - else: - matrix=matrix.reshape(0,0) - - s=matrix.shape - recordlength+=4*2+np.product(s)*8 #row and col of matrix + matrix of doubles - - #write length of record - fid.write(pack('q',recordlength)) - - #write data code: - fid.write(pack('i',FormatToCode(datatype))) - - #write data, first number of records - fid.write(pack('i',len(data))) - - for matrix in data: - if isinstance(matrix,(bool,int,float)): - matrix=np.array([matrix]) - elif isinstance(matrix,(list,tuple)): - matrix=np.array(matrix).reshape(-1,) - if np.ndim(matrix) == 1: - matrix=matrix.reshape(np.size(matrix),) - - s=matrix.shape - - if np.ndim(matrix) == 1: - fid.write(pack('i',s[0])) - fid.write(pack('i',1)) - for i in range(s[0]): - fid.write(pack('d',float(matrix[i]))) #get to the "c" convention, hence the transpose - else: - fid.write(pack('i',s[0])) - fid.write(pack('i',s[1])) - for i in range(s[0]): - for j in range(s[1]): - fid.write(pack('d',float(matrix[i][j]))) - # }}} - - elif datatype=='StringArray': # {{{ - - #first get length of record - recordlength=4+4 #for length of array + code - for string in data: - recordlength+=4+len(string) #for each string - - #write length of record - fid.write(pack('q',recordlength)) - - #write data code: - fid.write(pack('i',FormatToCode(datatype))) - - #now write length of string array - fid.write(pack('i',len(data))) - - #now write the strings - for string in data: - fid.write(pack('i',len(string))) - fid.write(pack('{}s'.format(len(string)),string.encode())) - # }}} - - else: # {{{ - raise TypeError('WriteData error message: data type: {} not supported yet! ({})'.format(datatype,enum)) - # }}} - -def FormatToCode(datatype): # {{{ - """ - This routine takes the datatype string, and hardcodes it into an integer, which - is passed along the record, in order to identify the nature of the dataset being - sent. - """ - - if datatype=='Boolean': - code=1 - elif datatype=='Integer': - code=2 - elif datatype=='Double': - code=3 - elif datatype=='String': - code=4 - elif datatype=='BooleanMat': - code=5 - elif datatype=='IntMat': - code=6 - elif datatype=='DoubleMat': - code=7 - elif datatype=='MatArray': - code=8 - elif datatype=='StringArray': - code=9 - elif datatype=='CompressedMat': - code=10 - else: - raise InputError('FormatToCode error message: data type not supported yet!') - - return code + +def WriteData(fid, prefix, *args): + """ + WRITEDATA - write model field in binary file + + Usage: + WriteData(fid, varargin) + """ + + #process options + options = pairoptions.pairoptions(*args) + + #Get data properties + if options.exist('object'): + #This is an object field, construct enum and data + obj = options.getfieldvalue('object') + fieldname = options.getfieldvalue('fieldname') + name = options.getfieldvalue('name', prefix + '.' + fieldname) + if options.exist('data'): + data = options.getfieldvalue('data') + else: + data = getattr(obj, fieldname) + else: + #No processing required + data = options.getfieldvalue('data') + name = options.getfieldvalue('name') + + datatype = options.getfieldvalue('format') + mattype = options.getfieldvalue('mattype', 0) #only required for matrices + timeserieslength = options.getfieldvalue('timeserieslength', -1) + + #Process sparse matrices +# if issparse(data), +# data = full(data); +# end + + #Scale data if necesarry + if options.exist('scale'): + data = np.array(data) + scale = options.getfieldvalue('scale') + if np.size(data) > 1 and np.ndim(data) > 1 and np.size(data, 0) == timeserieslength: + data[0:-1, :] = scale * data[0:-1, :] + else: + data = scale * data + if np.size(data) > 1 and np.size(data, 0) == timeserieslength: + yts = options.getfieldvalue('yts') + if np.ndim(data) > 1: + data[-1, :] = yts * data[-1, :] + else: + data[-1] = yts * data[-1] + + #Step 1: write the enum to identify this record uniquely + fid.write(pack('i', len(name))) + fid.write(pack('{}s'.format(len(name)), name.encode())) + + #Step 2: write the data itself. + if datatype == 'Boolean': # {{{ + #first write length of record + fid.write(pack('q', 4 + 4)) #1 bool (disguised as an int) + code + #write data code: + fid.write(pack('i', FormatToCode(datatype))) + + #now write bool as an integer + try: + fid.write(pack('i', int(data))) #send an int, not easy to send a bool + except error as Err: + raise ValueError('field {} cannot be marshaled, {}'.format(name, Err)) + # }}} + + elif datatype == 'Integer': # {{{ + #first write length of record + fid.write(pack('q', 4 + 4)) #1 integer + code + + #write data code: + fid.write(pack('i', FormatToCode(datatype))) + + #now write integer + try: + fid.write(pack('i', int(data))) #force an int, + except error as Err: + raise ValueError('field {} cannot be marshaled, {}'.format(name, Err)) + # }}} + + elif datatype == 'Double': # {{{ + #first write length of record + fid.write(pack('q', 8 + 4)) #1 double + code + + #write data code: + fid.write(pack('i', FormatToCode(datatype))) + + #now write double + try: + fid.write(pack('d', data)) + except error as Err: + raise ValueError('field {} cannot be marshaled, {}'.format(name, Err)) + # }}} + + elif datatype == 'String': # {{{ + #first write length of record + fid.write(pack('q', len(data) + 4 + 4)) #string + string size + code + + #write data code: + fid.write(pack('i', FormatToCode(datatype))) + + #now write string + fid.write(pack('i', len(data))) + fid.write(pack('{}s'.format(len(data)), data.encode())) + # }}} + + elif datatype in ['IntMat', 'BooleanMat']: # {{{ + if isinstance(data, (int, bool)): + data = np.array([data]) + elif isinstance(data, (list, tuple)): + data = np.array(data).reshape(-1,) + if np.ndim(data) == 1: + if np.size(data): + data = data.reshape(np.size(data),) + else: + data = data.reshape(0, 0) + + #Get size + s = data.shape + #if matrix = NaN, then do not write anything + if np.ndim(data) == 2 and np.product(s) == 1 and np.all(np.isnan(data)): + s = (0, 0) + + #first write length of record + recordlength = 4 + 4 + 8 * np.product(s) + 4 + 4 #2 integers (32 bits) + the double matrix + code + matrix type + fid.write(pack('q', recordlength)) + + #write data code and matrix type: + fid.write(pack('i', FormatToCode(datatype))) + fid.write(pack('i', mattype)) + + #now write matrix + fid.write(pack('i', s[0])) + try: + fid.write(pack('i', s[1])) + except IndexError: + fid.write(pack('i', 1)) + for i in range(s[0]): + if np.ndim(data) == 1: + fid.write(pack('d', float(data[i]))) #get to the "c" convention, hence the transpose + else: + for j in range(s[1]): + fid.write(pack('d', float(data[i][j]))) #get to the "c" convention, hence the transpose + # }}} + + elif datatype == 'DoubleMat': # {{{ + + if isinstance(data, (bool, int, float)): + data = np.array([data]) + elif isinstance(data, (list, tuple)): + data = np.array(data).reshape(-1,) + if np.ndim(data) == 1: + if np.size(data): + data = data.reshape(np.size(data),) + else: + data = data.reshape(0, 0) + + #Get size + s = data.shape + #if matrix = NaN, then do not write anything + if np.ndim(data) == 1 and np.product(s) == 1 and np.all(np.isnan(data)): + s = (0, 0) + + #first write length of record + recordlength = 4 + 4 + 8 * np.product(s) + 4 + 4 #2 integers (32 bits) + the double matrix + code + matrix type + + try: + fid.write(pack('q', recordlength)) + except error as Err: + raise ValueError('Field {} can not be marshaled, {}, with "number" the lenght of the record.'.format(name, Err)) + + #write data code and matrix type: + fid.write(pack('i', FormatToCode(datatype))) + fid.write(pack('i', mattype)) + #now write matrix + fid.write(pack('i', s[0])) + try: + fid.write(pack('i', s[1])) + except IndexError: + fid.write(pack('i', 1)) + for i in range(s[0]): + if np.ndim(data) == 1: + fid.write(pack('d', float(data[i]))) #get to the "c" convention, hence the transpose + else: + for j in range(s[1]): + fid.write(pack('d', float(data[i][j]))) #get to the "c" convention, hence the transpose + # }}} + + elif datatype == 'CompressedMat': # {{{ + if isinstance(data, (bool, int, float)): + data = np.array([data]) + elif isinstance(data, (list, tuple)): + data = np.array(data).reshape(-1,) + if np.ndim(data) == 1: + if np.size(data): + data = data.reshape(np.size(data),) + else: + data = data.reshape(0, 0) + + #Get size + s = data.shape + if np.ndim(data) == 1: + n2 = 1 + else: + n2 = s[1] + + #if matrix = NaN, then do not write anything + if np.ndim(data) == 1 and np.product(s) == 1 and np.all(np.isnan(data)): + s = (0, 0) + n2 = 0 + + #first write length of record + recordlength = 4 + 4 + 8 + 8 + 1 * (s[0] - 1) * n2 + 8 * n2 + 4 + 4 #2 integers (32 bits) + the matrix + code + matrix type + try: + fid.write(pack('q', recordlength)) + except error as Err: + raise ValueError('Field {} can not be marshaled, {}, with "number" the lenght of the record.'.format(name, Err)) + + #write data code and matrix type: + fid.write(pack('i', FormatToCode(datatype))) + fid.write(pack('i', mattype)) + + #Write offset and range + A = data[0:s[0] - 1] + offsetA = A.min() + rangeA = A.max() - offsetA + + if rangeA == 0: + A = A * 0 + else: + A = (A - offsetA) / rangeA * 255. + + #now write matrix + fid.write(pack('i', s[0])) + try: + fid.write(pack('i', s[1])) + except IndexError: + fid.write(pack('i', 1)) + fid.write(pack('d', float(offsetA))) + fid.write(pack('d', float(rangeA))) + if np.ndim(data) == 1: + for i in range(s[0] - 1): + fid.write(pack('B', int(A[i]))) + fid.write(pack('d', float(data[s[0] - 1]))) #get to the "c" convention, hence the transpose + + elif np.product(s) > 0: + for i in range(s[0] - 1): + for j in range(s[1]): + fid.write(pack('B', int(A[i][j]))) #get to the "c" convention, hence the transpose + + for j in range(s[1]): + fid.write(pack('d', float(data[s[0] - 1][j]))) + + # }}} + + elif datatype == 'MatArray': # {{{ + + #first get length of record + recordlength = 4 + 4 #number of records + code + for matrix in data: + if isinstance(matrix, (bool, int, float)): + matrix = np.array([matrix]) + elif isinstance(matrix, (list, tuple)): + matrix = np.array(matrix).reshape(-1,) + if np.ndim(matrix) == 1: + if np.size(matrix): + matrix = matrix.reshape(np.size(matrix),) + else: + matrix = matrix.reshape(0, 0) + + s = matrix.shape + recordlength += 4 * 2 + np.product(s) * 8 #row and col of matrix + matrix of doubles + + #write length of record + fid.write(pack('q', recordlength)) + + #write data code: + fid.write(pack('i', FormatToCode(datatype))) + + #write data, first number of records + fid.write(pack('i', len(data))) + + for matrix in data: + if isinstance(matrix, (bool, int, float)): + matrix = np.array([matrix]) + elif isinstance(matrix, (list, tuple)): + matrix = np.array(matrix).reshape(-1,) + if np.ndim(matrix) == 1: + matrix = matrix.reshape(np.size(matrix),) + + s = matrix.shape + + fid.write(pack('i', s[0])) + try: + fid.write(pack('i', s[1])) + except IndexError: + fid.write(pack('i', 1)) + for i in range(s[0]): + if np.ndim(matrix) == 1: + fid.write(pack('d', float(matrix[i]))) #get to the "c" convention, hence the transpose + else: + for j in range(s[1]): + fid.write(pack('d', float(matrix[i][j]))) + # }}} + + elif datatype == 'StringArray': # {{{ + #first get length of record + recordlength = 4 + 4 #for length of array + code + for string in data: + recordlength += 4 + len(string) #for each string + + #write length of record + fid.write(pack('q', recordlength)) + + #write data code: + fid.write(pack('i', FormatToCode(datatype))) + + #now write length of string array + fid.write(pack('i', len(data))) + + #now write the strings + for string in data: + fid.write(pack('i', len(string))) + fid.write(pack('{}s'.format(len(string)), string.encode())) + # }}} + + else: # {{{ + raise TypeError('WriteData error message: data type: {} not supported yet! ({})'.format(datatype, name)) + # }}} + + +def FormatToCode(datatype): # {{{ + """ + This routine takes the datatype string, and hardcodes it into an integer, which + is passed along the record, in order to identify the nature of the dataset being + sent. + """ + if datatype == 'Boolean': + code = 1 + elif datatype == 'Integer': + code = 2 + elif datatype == 'Double': + code = 3 + elif datatype == 'String': + code = 4 + elif datatype == 'BooleanMat': + code = 5 + elif datatype == 'IntMat': + code = 6 + elif datatype == 'DoubleMat': + code = 7 + elif datatype == 'MatArray': + code = 8 + elif datatype == 'StringArray': + code = 9 + elif datatype == 'CompressedMat': + code = 10 + else: + raise IOError('FormatToCode error message: data type not supported yet!') + return code # }}} Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py =================================================================== --- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 23771) +++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 23772) @@ -4,294 +4,301 @@ import results as resultsclass -def parseresultsfromdisk(md,filename,iosplit): - if iosplit: - saveres=parseresultsfromdiskiosplit(md,filename) - else: - saveres=parseresultsfromdiskioserial(md,filename) - return saveres - -def parseresultsfromdiskioserial(md,filename): # {{{ - #Open file - try: - fid=open(filename,'rb') - except IOError as e: - raise IOError("loadresultsfromdisk error message: could not open '{}' for binary reading.".format(filename)) - - #initialize results: - saveres=[] - - #Read fields until the end of the file. - loadres=ReadData(fid,md) - - counter=0 - check_nomoresteps=0 - step=loadres['step'] - - while loadres: - #check that the new result does not add a step, which would be an error: - if check_nomoresteps: - if loadres['step']>=1: - raise TypeError("parsing results for a steady-state core, which incorporates transient results!") - - #Check step, increase counter if this is a new step - if(step!=loadres['step'] and loadres['step']>1): - counter = counter + 1 - step = loadres['step'] - - #Add result - if loadres['step']==0: - #if we have a step = 0, this is a steady state solution, don't expect more steps. - index = 0; - check_nomoresteps=1 - elif loadres['step']==1: - index = 0 - else: - index = counter; - - if index > len(saveres)-1: - for i in range(len(saveres)-1,index-1): - saveres.append(None) - saveres.append(resultsclass.results()) - elif saveres[index] is None: - saveres[index]=resultsclass.results() - - #Get time and step - if loadres['step'] != -9999.: - saveres[index].__dict__['step']=loadres['step'] - if loadres['time'] != -9999.: - saveres[index].__dict__['time']=loadres['time'] - - #Add result - saveres[index].__dict__[loadres['fieldname']]=loadres['field'] - - #read next result - loadres=ReadData(fid,md) - - fid.close() - - return saveres - # }}} -def parseresultsfromdiskiosplit(md,filename): # {{{ - - #Open file - try: - fid=open(filename,'rb') - except IOError as e: - raise IOError("loadresultsfromdisk error message: could not open '{}' for binary reading.".format(filename)) - - saveres=[] - - #if we have done split I/O, ie, we have results that are fragmented across patches, - #do a first pass, and figure out the structure of results - loadres=ReadDataDimensions(fid) - while loadres: - - #Get time and step - if loadres['step'] > len(saveres): - for i in range(len(saveres),loadres['step']-1): - saveres.append(None) - saveres.append(resultsclass.results()) - setattr(saveres[loadres['step']-1],'step',loadres['step']) - setattr(saveres[loadres['step']-1],'time',loadres['time']) - - #Add result - setattr(saveres[loadres['step']-1],loadres['fieldname'],float('NaN')) - - #read next result - loadres=ReadDataDimensions(fid) - - #do a second pass, and figure out the size of the patches - fid.seek(0) #rewind - loadres=ReadDataDimensions(fid) - while loadres: - - #read next result - loadres=ReadDataDimensions(fid) - - #third pass, this time to read the real information - fid.seek(0) #rewind - loadres=ReadData(fid,md) - while loadres: - - #Get time and step - if loadres['step']> len(saveres): - for i in range(len(saveres),loadres['step']-1): - saveres.append(None) - saveres.append(saveresclass.saveres()) - setattr(saveres[loadres['step']-1],'step',loadres['step']) - setattr(saveres[loadres['step']-1],'time',loadres['time']) - - #Add result - setattr(saveres[loadres['step']-1],loadres['fieldname'],loadres['field']) - - #read next result - loadres=ReadData(fid,md) - - #close file - fid.close() - - return saveres - # }}} - -def ReadData(fid,md): # {{{ - """ - READDATA - ... - - Usage: - field=ReadData(fid,md) - """ - - #read field - try: - length=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - fieldname=struct.unpack('{}s'.format(length),fid.read(length))[0][:-1] - fieldname=fieldname.decode() #strings are booleans when stored so need to be converted back - time=struct.unpack('d',fid.read(struct.calcsize('d')))[0] - step=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - datatype=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - M=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - if datatype==1: - field=np.array(struct.unpack('{}d'.format(M),fid.read(M*struct.calcsize('d'))),dtype=float) - - elif datatype==2: - field=struct.unpack('{}s'.format(M),fid.read(M))[0][:-1] - field=field.decode() - - elif datatype==3: - N=struct.unpack('i',fid.read(struct.calcsize('i')))[0] -# field=transpose(fread(fid,[N M],'double')); - field=np.zeros(shape=(M,N),dtype=float) - for i in range(M): - field[i,:]=struct.unpack('{}d'.format(N),fid.read(N*struct.calcsize('d'))) - - elif datatype==4: - N=struct.unpack('i',fid.read(struct.calcsize('i')))[0] -# field=transpose(fread(fid,[N M],'int')); - field=np.zeros(shape=(M,N),dtype=int) - for i in range(M): - field[i,:]=struct.unpack('{}i'.format(N),fid.read(N*struct.calcsize('i'))) - - else: - raise TypeError("cannot read data of datatype {}".format(datatype)) - - #Process units here FIXME: this should not be done here! - yts=md.constants.yts - if fieldname=='BalancethicknessThickeningRate': - field = field*yts - elif fieldname=='HydrologyWaterVx': - field = field*yts - elif fieldname=='HydrologyWaterVy': - field = field*yts - elif fieldname=='Vx': - field = field*yts - elif fieldname=='Vy': - field = field*yts - elif fieldname=='Vz': - field = field*yts - elif fieldname=='Vel': - field = field*yts - elif fieldname=='BasalforcingsGroundediceMeltingRate': - field = field*yts - elif fieldname=='BasalforcingsFloatingiceMeltingRate': - field = field*yts - elif fieldname=='TotalFloatingBmb': - field = field/10.**12*yts #(GigaTon/year) - elif fieldname=='TotalFloatingBmbScaled': - field = field/10.**12*yts #(GigaTon/year) - elif fieldname=='TotalGroundedBmb': - field = field/10.**12*yts #(GigaTon/year) - elif fieldname=='TotalGroundedBmbScaled': - field = field/10.**12*yts #(GigaTon/year) - elif fieldname=='TotalSmb': - field = field/10.**12*yts #(GigaTon/year) - elif fieldname=='TotalSmbScaled': - field = field/10.**12*yts #(GigaTon/year) - elif fieldname=='SmbMassBalance': - field = field*yts - elif fieldname=='SmbPrecipitation': - field = field*yts - elif fieldname=='SmbRunoff': - field = field*yts - elif fieldname=='SmbEvaporation': - field = field*yts; - elif fieldname=='SmbRefreeze': - field = field*yts; - elif fieldname=='SmbEC': - field = field*yts - elif fieldname=='SmbAccumulation': - field = field*yts - elif fieldname=='SmbMelt': - field = field*yts - elif fieldname=='SmbMAdd': - field = field*yts - elif fieldname=='CalvingCalvingrate': - field = field*yts - elif fieldname=='LoveKernelsReal' or fieldname=='LoveKernelsImag': - nlayer = md.materials.numlayers - degmax = md.love.sh_nmax - nfreq = md.love.nfreq - #for numpy 1.8+ only - #temp_field = np.full((degmax+1,nfreq,nlayer+1,6),0.0) - temp_field = np.empty((degmax+1,nfreq,nlayer+1,6)) - temp_field.fill(0.0) - for ii in range(degmax+1): - for jj in range(nfreq): - for kk in range(nlayer+1): - ll = ii*(nlayer+1)*6 + (kk*6+1) - for mm in range(6): - temp_field[ii,jj,kk,mm] = field[ll+mm-1,jj] - field=temp_field - - if time !=-9999: - time = time/yts - - saveres=OrderedDict() - saveres['fieldname']=fieldname - saveres['time']=time - saveres['step']=step - saveres['field']=field - - except struct.error as e: - saveres=None - - return saveres - # }}} + +def parseresultsfromdisk(md, filename, iosplit): + if iosplit: + saveres = parseresultsfromdiskiosplit(md, filename) + else: + saveres = parseresultsfromdiskioserial(md, filename) + return saveres + + +def parseresultsfromdiskioserial(md, filename): # {{{ + #Open file + try: + fid = open(filename, 'rb') + except IOError as e: + raise IOError("loadresultsfromdisk error message: could not open '{}' for binary reading.".format(filename)) + + #initialize results: + saveres = [] + + #Read fields until the end of the file. + loadres = ReadData(fid, md) + + counter = 0 + check_nomoresteps = 0 + step = loadres['step'] + + while loadres: + #check that the new result does not add a step, which would be an error: + if check_nomoresteps: + if loadres['step'] >= 1: + raise TypeError("parsing results for a steady-state core, which incorporates transient results!") + + #Check step, increase counter if this is a new step + if(step != loadres['step'] and loadres['step'] > 1): + counter = counter + 1 + step = loadres['step'] + + #Add result + if loadres['step'] == 0: + #if we have a step = 0, this is a steady state solution, don't expect more steps. + index = 0 + check_nomoresteps = 1 + elif loadres['step'] == 1: + index = 0 + else: + index = counter + + if index > len(saveres) - 1: + for i in range(len(saveres) - 1, index - 1): + saveres.append(None) + saveres.append(resultsclass.results()) + elif saveres[index] is None: + saveres[index] = resultsclass.results() + + #Get time and step + if loadres['step'] != -9999.: + saveres[index].__dict__['step'] = loadres['step'] + if loadres['time'] != -9999.: + saveres[index].__dict__['time'] = loadres['time'] + + #Add result + saveres[index].__dict__[loadres['fieldname']] = loadres['field'] + + #read next result + loadres = ReadData(fid, md) + + fid.close() + + return saveres + # }}} + + +def parseresultsfromdiskiosplit(md, filename): # {{{ + #Open file + try: + fid = open(filename, 'rb') + except IOError as e: + raise IOError("loadresultsfromdisk error message: could not open '{}' for binary reading.".format(filename)) + + saveres = [] + + #if we have done split I/O, ie, we have results that are fragmented across patches, + #do a first pass, and figure out the structure of results + loadres = ReadDataDimensions(fid) + while loadres: + + #Get time and step + if loadres['step'] > len(saveres): + for i in range(len(saveres), loadres['step'] - 1): + saveres.append(None) + saveres.append(resultsclass.results()) + setattr(saveres[loadres['step'] - 1], 'step', loadres['step']) + setattr(saveres[loadres['step'] - 1], 'time', loadres['time']) + + #Add result + setattr(saveres[loadres['step'] - 1], loadres['fieldname'], float('NaN')) + + #read next result + loadres = ReadDataDimensions(fid) + + #do a second pass, and figure out the size of the patches + fid.seek(0) #rewind + loadres = ReadDataDimensions(fid) + while loadres: + + #read next result + loadres = ReadDataDimensions(fid) + + #third pass, this time to read the real information + fid.seek(0) #rewind + loadres = ReadData(fid, md) + while loadres: + + #Get time and step + if loadres['step'] > len(saveres): + for i in range(len(saveres), loadres['step'] - 1): + saveres.append(None) + saveres.append(saveresclass.saveres()) + setattr(saveres[loadres['step'] - 1], 'step', loadres['step']) + setattr(saveres[loadres['step'] - 1], 'time', loadres['time']) + + #Add result + setattr(saveres[loadres['step'] - 1], loadres['fieldname'], loadres['field']) + + #read next result + loadres = ReadData(fid, md) + + #close file + fid.close() + + return saveres + # }}} + + +def ReadData(fid, md): # {{{ + """ + READDATA - + + Usage: + field = ReadData(fid, md) + """ + + #read field + try: + length = struct.unpack('i', fid.read(struct.calcsize('i')))[0] + fieldname = struct.unpack('{}s'.format(length), fid.read(length))[0][:-1] + fieldname = fieldname.decode() #strings are binaries when stored so need to be converted back + time = struct.unpack('d', fid.read(struct.calcsize('d')))[0] + step = struct.unpack('i', fid.read(struct.calcsize('i')))[0] + datatype = struct.unpack('i', fid.read(struct.calcsize('i')))[0] + M = struct.unpack('i', fid.read(struct.calcsize('i')))[0] + if datatype == 1: + field = np.array(struct.unpack('{}d'.format(M), fid.read(M * struct.calcsize('d'))), dtype=float) + + elif datatype == 2: + field = struct.unpack('{}s'.format(M), fid.read(M))[0][:-1] + field = field.decode() + + elif datatype == 3: + N = struct.unpack('i', fid.read(struct.calcsize('i')))[0] + #field = transpose(fread(fid, [N M], 'double')); + field = np.zeros(shape=(M, N), dtype=float) + for i in range(M): + field[i, :] = struct.unpack('{}d'.format(N), fid.read(N * struct.calcsize('d'))) + + elif datatype == 4: + N = struct.unpack('i', fid.read(struct.calcsize('i')))[0] + # field = transpose(fread(fid, [N M], 'int')); + field = np.zeros(shape=(M, N), dtype=int) + for i in range(M): + field[i, :] = struct.unpack('{}i'.format(N), fid.read(N * struct.calcsize('i'))) + + else: + raise TypeError("cannot read data of datatype {}".format(datatype)) + + #Process units here FIXME: this should not be done here! + yts = md.constants.yts + if fieldname == 'BalancethicknessThickeningRate': + field = field * yts + elif fieldname == 'HydrologyWaterVx': + field = field * yts + elif fieldname == 'HydrologyWaterVy': + field = field * yts + elif fieldname == 'Vx': + field = field * yts + elif fieldname == 'Vy': + field = field * yts + elif fieldname == 'Vz': + field = field * yts + elif fieldname == 'Vel': + field = field * yts + elif fieldname == 'BasalforcingsGroundediceMeltingRate': + field = field * yts + elif fieldname == 'BasalforcingsFloatingiceMeltingRate': + field = field * yts + elif fieldname == 'TotalFloatingBmb': + field = field / 10.**12 * yts #(GigaTon/year) + elif fieldname == 'TotalFloatingBmbScaled': + field = field / 10.**12 * yts #(GigaTon/year) + elif fieldname == 'TotalGroundedBmb': + field = field / 10.**12 * yts #(GigaTon/year) + elif fieldname == 'TotalGroundedBmbScaled': + field = field / 10.**12 * yts #(GigaTon/year) + elif fieldname == 'TotalSmb': + field = field / 10.**12 * yts #(GigaTon/year) + elif fieldname == 'TotalSmbScaled': + field = field / 10.**12 * yts #(GigaTon/year) + elif fieldname == 'SmbMassBalance': + field = field * yts + elif fieldname == 'SmbPrecipitation': + field = field * yts + elif fieldname == 'SmbRunoff': + field = field * yts + elif fieldname == 'SmbEvaporation': + field = field * yts + elif fieldname == 'SmbRefreeze': + field = field * yts + elif fieldname == 'SmbEC': + field = field * yts + elif fieldname == 'SmbAccumulation': + field = field * yts + elif fieldname == 'SmbMelt': + field = field * yts + elif fieldname == 'SmbMAdd': + field = field * yts + elif fieldname == 'CalvingCalvingrate': + field = field * yts + elif fieldname == 'LoveKernelsReal' or fieldname == 'LoveKernelsImag': + nlayer = md.materials.numlayers + degmax = md.love.sh_nmax + nfreq = md.love.nfreq + #for numpy 1.8+ only + #temp_field = np.full((degmax+1, nfreq, nlayer+1, 6), 0.0) + temp_field = np.empty((degmax + 1, nfreq, nlayer + 1, 6)) + temp_field.fill(0.0) + for ii in range(degmax + 1): + for jj in range(nfreq): + for kk in range(nlayer + 1): + ll = ii * (nlayer + 1) * 6 + (kk * 6 + 1) + for mm in range(6): + temp_field[ii, jj, kk, mm] = field[ll + mm - 1, jj] + field = temp_field + + if time != -9999: + time = time / yts + + saveres = OrderedDict() + saveres['fieldname'] = fieldname + saveres['time'] = time + saveres['step'] = step + saveres['field'] = field + + except struct.error as e: + saveres = None + + return saveres + # }}} + + def ReadDataDimensions(fid): # {{{ - """ - READDATADIMENSIONS - read data dimensions, step and time, but not the data itself. - - Usage: - field=ReadDataDimensions(fid) - """ - - #read field - try: - length=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - fieldname=struct.unpack('{}s'.format(length),fid.read(length))[0][:-1] - time=struct.unpack('d',fid.read(struct.calcsize('d')))[0] - step=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - dtattype=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - M=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - N=1 #default - if datatype==1: - fid.seek(M*8,1) - elif datatype==2: - fid.seek(M,1) - elif datatype==3: - N=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - fid.seek(N*M*8,1) - else: - raise TypeError("cannot read data of datatype {}".format(datatype)) - - saveres=OrderedDict() - saveres['fieldname']=fieldname - saveres['time']=time - saveres['step']=step - saveres['M']=M - saveres['N']=N - - except struct.error as e: - saveres=None - - return saveres - # }}} + """ + READDATADIMENSIONS - read data dimensions, step and time, but not the data itself. + + Usage: + field = ReadDataDimensions(fid) + """ + + #read field + try: + length = struct.unpack('i', fid.read(struct.calcsize('i')))[0] + fieldname = struct.unpack('{}s'.format(length), fid.read(length))[0][:-1] + time = struct.unpack('d', fid.read(struct.calcsize('d')))[0] + step = struct.unpack('i', fid.read(struct.calcsize('i')))[0] + datatype = struct.unpack('i', fid.read(struct.calcsize('i')))[0] + M = struct.unpack('i', fid.read(struct.calcsize('i')))[0] + N = 1 #default + if datatype == 1: + fid.seek(M * 8, 1) + elif datatype == 2: + fid.seek(M, 1) + elif datatype == 3: + N = struct.unpack('i', fid.read(struct.calcsize('i')))[0] + fid.seek(N * M * 8, 1) + else: + raise TypeError("cannot read data of datatype {}".format(datatype)) + + saveres = OrderedDict() + saveres['fieldname'] = fieldname + saveres['time'] = time + saveres['step'] = step + saveres['M'] = M + saveres['N'] = N + + except struct.error as Err: + print(Err) + saveres = None + + return saveres + # }}}