Index: /issm/trunk-jpl/src/m/netcdf/README.txt =================================================================== --- /issm/trunk-jpl/src/m/netcdf/README.txt (revision 27891) +++ /issm/trunk-jpl/src/m/netcdf/README.txt (revision 27891) @@ -0,0 +1,87 @@ +The write_netCDF and read_netCDF modules provide a convenient way to save and restore the state of a model class instance +in binary format via NetCDF4. This allows users to store the class state on disk and retrieve it later, facilitating seamless +transitions between Python and MATLAB environments. + +To save a model, call either write_netCDF.py or write_netCDF.m depending on whether your class is in matlab or python. +To read a saved model, call either read_netCDF.py or read_netCDF.m depending on what language you prefer to use the model in. +If you would like to log the names and locations of variables being stored, add the argument verbose = True (verbose = true for matlab). + +Usage Instructions: + + Python: + - Saving a model: + from write_netCDF import write_netCDF + + md = bamg(model(), foo.csv, .01) + + write_netCDF(md, 'adress_to_save/../filename.nc') + + - Reading a model: + from read_netCDF import read_netCDF + + md = read_netCDF('adress_to_file/../filename.nc') + + Verbose examples: + write_netCDF(md, adress_to_save/../filename.nc, verbose = True) + md = read_netCDF(adress_to_file/../filename.nc, verbose = True) + + MATLAB: + - Saving a model: + + write_netCDF(md, adress_to_save/../filename.nc); + + - Reading a model: + + md = read_netCDF(adress_to_file/../filename.nc); + + Verbose examples: + write_netCDF(md, adress_to_save/../filename.nc, verbose = true); + + or: + + write_netCDF(md, adress_to_save/../filename.nc, verbose); + md = read_netCDF(adress_to_file/../filename.nc, verbose = true); + +Dependencies: + Python: + - NumPy + - NetCDF4 / NetCDF4.Dataset + - The model() class + - results.solution / results.solutionstep / results.resultsdakota + - inversion.inversion / inversion.m1qn3inversion / inversion.taoinversion + + MATLAB: + - The model() class + - inversion.inversion / inversion.m1qn3inversion / inversion.taoinversion + + +Additional Information: + +There are currently datatypes that both write_netCDF and read_netCDF modules may not be able to handle. These datatypes might +include lists with multiple datatypes (ie, ['number', 1, 'letter', a, 'color', 'blue']), lists of dicts ect. + +To add functionality for these additional cases, one must simply create a function to handle the case and call it using a +conditional case within the create_var() function. To read the data from the NetCDF4 file, add the case to the +copy_variable_data_to_new_model() function in read_netCDF so that the data can be added to a new model() instance. + +Known issues: + +Unlike Python, MATLAB doesn't utilize subclasses in its model class. This leads to a loss of certain subclass instances. +For instance, the results.solutionstep() class poses a known issue. In MATLAB, there's no direct equivalent. The fields in +'md.results' in MATLAB might correspond to instances of resultsdakota(), solution(), or solutionstep() in Python, but +because those classes don't exist in MATLAB, there is no way for python to know which instance it needs. + +The current workaround, while not theoretically sound, involves searching for the class name string in MATLAB's 'results' +field names. For instance, 'md.results.TransientSolution' is recorded as a solution() class instance. However, problems arise +in cases like 'md.results.StressbalanceSolution', where the code notes a solution() instance, while in Python, it should be a +solutionstep() instance. + +So far, there have been no recorded problems swapping a solutionstep() instance for a solution() instance. + +Potential solutions are: + + - Restructure both Python and MATLAB solve frameworks. In Python, when creating an md.results. instance, + embed 'solutionstep' in the class instance name. + >> This solution is very involved, and would include the tedious modification of >5 files in total + - Create a hash table linking solutions with their corresponding 'md.results.' for reference when saving models to + the netCDF file. Index: /issm/trunk-jpl/src/m/netcdf/read_netCDF.m =================================================================== --- /issm/trunk-jpl/src/m/netcdf/read_netCDF.m (revision 27891) +++ /issm/trunk-jpl/src/m/netcdf/read_netCDF.m (revision 27891) @@ -0,0 +1,326 @@ +%{ +Given a NetCDF4 file, this set of functions will perform the following: + 1. Enter each group of the file. + 2. For each variable in each group, update an empty model with the variable's data + 3. Enter nested groups and repeat + + +If the model you saved has subclass instances that are not in the standard model() class +you can: + 1. Copy lines 30-35, set the "results" string to the name of the subclass instance, + 2. Copy and modify the make_results_subclasses() function to create the new subclass + instances you need. +From there, the rest of this script will automatically create the new subclass +instance in the model you're writing to and store the data from the netcdf file there. +%} + + +function model_copy = read_netCDF(filename, varargin) + if nargin > 1 + verbose = true; + else + verbose = false; + end + + if verbose + fprintf('NetCDF42C v1.1.14\n'); + end + % make a model framework to fill that is in the scope of this file + model_copy = model(); + + % Check if path exists + if exist(filename, 'file') + if verbose + fprintf('Opening %s for reading\n', filename); + end + + % Open the given netCDF4 file + NCData = netcdf.open(filename, 'NOWRITE'); + % Remove masks from netCDF data for easy conversion: NOT WORKING + %netcdf.setMask(NCData, 'NC_NOFILL'); + + % see if results is in there, if it is we have to instantiate some classes + try + results_group_id = netcdf.inqNcid(NCData, "results"); + model_copy = make_results_subclasses(model_copy, NCData, verbose); + catch + end % 'results' group doesn't exist + + % see if inversion is in there, if it is we may have to instantiate some classes + try + inversion_group_id = netcdf.inqNcid(NCData, "inversion"); + model_copy = check_inversion_class(model_copy, NCData, verbose); + catch + end % 'inversion' group doesn't exist + + % loop over first layer of groups in netcdf file + for group = netcdf.inqGrps(NCData) + group_id = netcdf.inqNcid(NCData, netcdf.inqGrpName(group)); + %disp(netcdf.inqGrpNameFull(group_id)) + % hand off first level to recursive search + model_copy = walk_nested_groups(group_id, model_copy, NCData, verbose); + end + + % Close the netCDF file + netcdf.close(NCData); + if verbose + disp('Model Successfully Copied') + end + else + fprintf('File %s does not exist.\n', filename); + end +end + + +function model_copy = make_results_subclasses(model_copy, NCData, verbose) + resultsGroup = netcdf.inqNcid(NCData, "results"); + variables = netcdf.inqVarIDs(resultsGroup); + for name = variables + class_instance = netcdf.inqVar(resultsGroup, name); + class_instance_names_raw = netcdf.getVar(resultsGroup, name, 'char').'; + class_instance_names = cellstr(class_instance_names_raw); + for index = 1:numel(class_instance_names) + class_instance_name = class_instance_names{index}; + model_copy.results = setfield(model_copy.results, class_instance_name, struct()); + end + %model_copy.results = setfield(model_copy.results, class_instance, class_instance_name); + end + model_copy = model_copy; + if verbose + disp('Successfully recreated results structs:') + for fieldname = string(fieldnames(model_copy.results)) + disp(fieldname) + end + end +end + + +function model_copy = check_inversion_class(model_copy, NCData, verbose) + % get the name of the inversion class: either inversion or m1qn3inversion or taoinversion + inversionGroup = netcdf.inqNcid(NCData, "inversion"); + varid = netcdf.inqVarID(inversionGroup, 'inversion_class_name'); + inversion_class = convertCharsToStrings(netcdf.getVar(inversionGroup, varid,'char')); + if strcmp(inversion_class, 'm1qn3inversion') + model_copy.inversion = m1qn3inversion(); + if verbose + disp('Successfully created inversion class instance: m1qn3inversion') + end + elseif strcmp(inversion_class, 'taoinversion') + model_copy.inversion = taoinversion(); + if verbose + disp('Successfully created inversion class instance: taoinversion') + end + else + if verbose + disp('No inversion class was found') + end + end + model_copy = model_copy; +end + + +function model_copy = walk_nested_groups(group_location_in_file, model_copy, NCData, verbose) + % we search the current group level for variables by getting this struct + variables = netcdf.inqVarIDs(group_location_in_file); + + % from the variables struct get the info related to the variables + for variable = variables + [varname, xtype, dimids, numatts] = netcdf.inqVar(group_location_in_file, variable); + + % keep an eye out for nested structs: + if strcmp(varname, 'this_is_a_nested') + is_nested = true; + model_copy = copy_nested_struct(group_location_in_file, model_copy, NCData, verbose); + elseif strcmp(varname, 'solution') + % band-aid pass.. + else + model_copy = copy_variable_data_to_new_model(group_location_in_file, varname, xtype, model_copy, NCData, verbose); + end + end + + % try to find groups in current level, if it doesn't work it's because there is nothing there + %try + % if it's a nested struct the function copy_nested_struct has already been called + if logical(exist('is_nested', 'var')) + % do nothing + else + % search for nested groups in the current level to feed back to this function + groups = netcdf.inqGrps(group_location_in_file); + if not(isempty(groups)) + for group = groups + group_id = netcdf.inqNcid(group_location_in_file, netcdf.inqGrpName(group)); + %disp(netcdf.inqGrpNameFull(group_id)) + model_copy = walk_nested_groups(group, model_copy, NCData, verbose); + end + end + end + %catch % no nested groups here + %end +end + + + +function model_copy = copy_nested_struct(group_location_in_file, model_copy, NCData, verbose) + %{ + A common multidimensional struct array is the 1xn md.results.TransientSolution struct. + The process to recreate is as follows: + 1. Get the name of the struct from group name + 2. Get the fieldnames from the subgroups + 3. Recreate the struct with fieldnames + 4. Populate the fields with their respective values + %} + + % step 1 + name_of_struct = netcdf.inqGrpName(group_location_in_file); + + % step 2 + subgroups = netcdf.inqGrps(group_location_in_file); % numerical cell array with ID's of subgroups + % get single subgroup's data + single_subgroup_ID = subgroups(1); + subgroup_varids = netcdf.inqVarIDs(single_subgroup_ID); + fieldnames = {}; + for variable = subgroup_varids + [varname, xtype, dimids, numatts] = netcdf.inqVar(single_subgroup_ID, variable); + fieldnames{end+1} = varname; + end + + % step 3 + address_in_model_raw = split(netcdf.inqGrpNameFull(group_location_in_file), '/'); + address_in_model = address_in_model_raw{2}; + + % we cannot assign a variable to represent this object as MATLAB treats all variables as copies + % and not pointers to the same memory address + % this means that if address_in_model has more than 1 layer, we need to modify the code. For now, + % we just hope this will do. An example of a no-solution would be model().abc.def.ghi.field + + model_copy.(address_in_model).(name_of_struct) = struct(); + % for every fieldname in the subgroup, create an empty field + for fieldname = string(fieldnames) + model_copy.(address_in_model).(name_of_struct).(fieldname) = {}; + end + + % use repmat to make the struct array multidimensional along the fields axis + number_of_dimensions = numel(subgroups); + model_copy.(address_in_model).(name_of_struct) = repmat(model_copy.(address_in_model).(name_of_struct), 1, number_of_dimensions); + + % step 4 + % for every layer of the multidimensional struct array, populate the fields + for current_layer = 1:number_of_dimensions + % choose subgroup + current_layer_subgroup_ID = subgroups(current_layer); + % get all vars + current_layer_subgroup_varids = netcdf.inqVarIDs(current_layer_subgroup_ID); + % get individual vars and set fields at layer current_layer + for varid = current_layer_subgroup_varids + [varname, xtype, dimids, numatts] = netcdf.inqVar(current_layer_subgroup_ID, varid); + data = netcdf.getVar(current_layer_subgroup_ID, varid); + + % netcdf uses Row Major Order but MATLAB uses Column Major Order so we need to transpose all arrays w/ more than 1 dim + if all(size(data)~=1) || xtype == 2 + data = data.'; + end + + % set the field + model_copy.(address_in_model).(name_of_struct)(current_layer).(varname) = data; + %address_to_struct_in_model = setfield(address_to_struct_in_model(current_layer), varname, data) + end + model_copy.(address_in_model).(name_of_struct)(current_layer); + if verbose + fprintf("Successfully loaded layer %s to multidimension struct array\n", num2str(current_layer)) + end + end + model_copy = model_copy; + if verbose + fprintf('Successfully recreated multidimensional structure array %s in md.%s\n', name_of_struct, address_in_model) + end +end + + + + +%{ +Since there are two types of objects that MATLAB uses (classes and structs), we have to check +which object we're working with before we can set any fields/attributes of it. After this is completed, +we can write the data to that location in the model. +%} + +function model_copy = copy_variable_data_to_new_model(group_location_in_file, varname, xtype, model_copy, NCData, verbose) + %disp(varname) + % this is an inversion band-aid + if strcmp(varname, 'inversion_class_name') || strcmp(varname, 'name_of_struct') || strcmp(varname, 'solution') + % we don't need this + else + % putting try/catch here so that any errors generated while copying data are logged and not lost by the try/catch in walk_nested_groups function + try + %disp(netcdf.inqGrpNameFull(group_location_in_file)) + %disp(class(netcdf.inqGrpNameFull(group_location_in_file))) + address_to_attr = strrep(netcdf.inqGrpNameFull(group_location_in_file), '/', '.'); + varid = netcdf.inqVarID(group_location_in_file, varname); + data = netcdf.getVar(group_location_in_file, varid); + + + % if we have an empty string + if xtype == 2 && isempty(all(data)) + data = cell(char()); + % if we have an empty cell-char array + elseif numel(data) == 1 && xtype == 3 && data == -32767 + data = cell(char()); + elseif isempty(all(data)) + data = [] + end + % band-aid for some cell-char-arrays: + if xtype == 2 && strcmp(data, 'default') + data = {'default'}; + end + + % netcdf uses Row Major Order but MATLAB uses Column Major Order so we need to transpose all arrays w/ more than 1 dim + if all(size(data)~=1) || xtype == 2 + data = data.'; + end + + % if we have a list of strings + if xtype == 2 + try + if strcmp(netcdf.getAtt(group_location_in_file, varid, "type_is"), 'cell_array_of_strings') + data = cellstr(data); + end + catch + % no attr found so we pass + end + end + + % the issm c compiler does not work with int64 datatypes, so we need to convert those to int16 + % reference this (very hard to find) link for netcdf4 datatypes: https://docs.unidata.ucar.edu/netcdf-c/current/netcdf_8h_source.html + %xtype + if xtype == 10 + arg_to_eval = ['model_copy', address_to_attr, '.', varname, ' = ' , 'double(data);']; + eval(arg_to_eval); + %disp('Loaded int64 as int16') + else + arg_to_eval = ['model_copy', address_to_attr, '.', varname, ' = data;']; + eval(arg_to_eval); + end + + if verbose + full_addy = netcdf.inqGrpNameFull(group_location_in_file); + %disp(xtype) + %class(data) + fprintf('Successfully loaded %s to %s\n', varname, full_addy); + end + + catch ME %ME is an MException struct + % Some error occurred if you get here. + fprintf(1,'There was an error with %s! \n', varname) + errorMessage = sprintf('Error in function %s() at line %d.\n\nError Message:\n%s', ME.stack.name, ME.stack.line, ME.message); + fprintf(1, '%s\n', errorMessage); + uiwait(warndlg(errorMessage)); + %line = ME.stack.line + %fprintf(1,'There was an error with %s! \n', varname) + %fprintf('The message was:\n%s\n',ME.message); + %fprintf(1,'The identifier was:\n%s\n',ME.identifier); + + % more error handling... + end + end + model_copy = model_copy; +end Index: /issm/trunk-jpl/src/m/netcdf/read_netCDF.py =================================================================== --- /issm/trunk-jpl/src/m/netcdf/read_netCDF.py (revision 27891) +++ /issm/trunk-jpl/src/m/netcdf/read_netCDF.py (revision 27891) @@ -0,0 +1,301 @@ +# imports +from netCDF4 import Dataset +import numpy as np +import numpy.ma as ma +from os import path, remove +from model import * +import re +from results import * +from m1qn3inversion import m1qn3inversion +from taoinversion import taoinversion +from collections import OrderedDict + + +''' +Given a NetCDF4 file, this set of functions will perform the following: + 1. Enter each group of the file. + 2. For each variable in each group, update an empty model with the variable's data + 3. Enter nested groups and repeat +''' + + +# make a model framework to fill that is in the scope of this file +model_copy = model() + +def read_netCDF(filename, verbose = False): + if verbose: + print('NetCDF42C v1.1.13') + + # check if path exists + if path.exists(filename): + if verbose: + print('Opening {} for reading'.format(filename)) + else: pass + + # open the given netCDF4 file + NCData = Dataset(filename, 'r') + # remove masks from numpy arrays for easy conversion + NCData.set_auto_mask(False) + else: + print('The file you entered does not exist or cannot be found in the current directory') + return print() + + # continuation of band-aid for results class + try: + NCData.groups['results'] + make_results_subclasses(NCData, verbose) + except: + pass + + # similarly, we need to check and see if we have an m1qn3inversion class instance + try: + NCData.groups['inversion'] + check_inversion_class(NCData, verbose) + except: + pass + + # walk through each group looking for subgroups and variables + for group in NCData.groups.keys(): + if 'debris' in group: + pass + else: + # have to send a custom name to this function: filename.groups['group'] + name = "NCData.groups['" + str(group) + "']" + walk_nested_groups(name, NCData, verbose) + + if verbose: + print("Model Successfully Loaded.") + return model_copy + + +def make_results_subclasses(NCData, verbose = False): + ''' + There are 3 possible subclasses: solution, solutionstep, resultsdakota. + In the NetCDF file these are saved as a list of strings. Ie, say there are 2 + instances of solution under results, StressbalanceSolution and TransientSolution. + In the NetCDF file we would see solution = "StressbalanceSolution", "TransientSolution" + To deconstruct this, we need to iteratively assign md.results.StressbalanceSolution = solution() + and md.results.TransientSolution = solution() and whatever else. + ''' + # start with the subclasses + for subclass in NCData.groups['results'].variables.keys(): + class_instance = subclass + '()' + + # now handle the instances + for instance in NCData.groups['results'].variables[subclass][:]: + # this is an ndarray of numpy bytes_ that we have to convert to strings + class_instance_name = instance.tobytes().decode('utf-8').strip() + # from here we can make new subclasses named as they were in the model saved + setattr(model_copy.results, class_instance_name, eval(class_instance)) + if verbose: + print(f'Successfully created results subclass instance {class_instance} named {class_instance_name}.') + + +def check_inversion_class(NCData, verbose = False): + # get the name of the inversion class: either inversion or m1qn3inversion or taoinversion + inversion_class_is = NCData.groups['inversion'].variables['inversion_class_name'][:][...].tobytes().decode() + if inversion_class_is == 'm1qn3inversion': + # if it is m1qn3inversion we need to instantiate that class since it's not native to model() + model_copy.inversion = m1qn3inversion(model_copy.inversion) + if verbose: + print('Conversion successful') + elif inversion_class_is == 'taoinversion': + # if it is taoinversion we need to instantiate that class since it's not native to model() + model_copy.inversion = taoinverion() + if verbose: + print('Conversion successful') + else: pass + + +def walk_nested_groups(group_location_in_file, NCData, verbose = False): + # first, we enter the group by: filename.groups['group_name'] + # second we search the current level for variables: filename.groups['group_name'].variables.keys() + # at this step we check for multidimensional structure arrays and filter them out + # third we get nested group keys by: filename.groups['group_name'].groups.keys() + # if a variables exists, copy the data to the model framework by calling copy function + # if a nested groups exist, repeat all + + for variable in eval(group_location_in_file + '.variables.keys()'): + if variable == 'this_is_a_nested' and 'results' in group_location_in_file: + # have to do some string deconstruction to get the name of the class instance/last group from 'NetCDF.groups['group1'].groups['group1.1']' + pattern = r"\['(.*?)'\]" + matches = re.findall(pattern, group_location_in_file) + name_of_struct = matches[-1] #eval(group_location_in_file + ".variables['solution']") + copy_multidimensional_results_struct(group_location_in_file, name_of_struct, NCData) + istruct = True + + elif variable == 'this_is_a_nested' and 'qmu' in group_location_in_file: + print('encountered qmu structure that is not yet supported.') + # have to do some string deconstruction to get the name of the class instance/last group from 'NetCDF.groups['group1'].groups['group1.1']' + #pattern = r"\['(.*?)'\]" + #matches = re.findall(pattern, group_location_in_file) + #name_of_struct = matches[-1] #eval(group_location_in_file + ".variables['solution']") + #name_of_struct = eval(group_location_in_file + ".variables['']") + #copy_multidimensional_qmu_struct(group_location_in_file, name_of_struct) + isstruct = True + + else: + location_of_variable_in_file = group_location_in_file + ".variables['" + str(variable) + "']" + # group_location_in_file is like filename.groups['group1'].groups['group1.1'].groups['group1.1.1'] + # Define the regex pattern to match the groups within brackets + pattern = r"\['(.*?)'\]" + # Use regex to find all matches and return something like 'group1.group1.1.group1.1.1 ...' where the last value is the name of the variable + matches = re.findall(pattern, location_of_variable_in_file) + variable_name = matches[-1] + location_of_variable_in_model = '.'.join(matches[:-1]) + copy_variable_data_to_new_model(location_of_variable_in_file, location_of_variable_in_model, variable_name, NCData, verbose=verbose) + + if 'istruct' in locals(): + pass + else: + for nested_group in eval(group_location_in_file + '.groups.keys()'): + new_nested_group = group_location_in_file + ".groups['" + str(nested_group) + "']" + walk_nested_groups(new_nested_group, NCData, verbose=verbose) + + + +''' + MATLAB has Multidimensional Structure Arrays in 2 known classes: results and qmu. + The python classes results.py and qmu.py emulate this MATLAB object in their own + unique ways. The functions in this script will assign data to either of these + classes such that the final structure is compatible with its parent class. +''' + +def copy_multidimensional_results_struct(group_location_in_file, name_of_struct, NCData, verbose = False): + ''' + A common multidimensional array is the 1xn md.results.TransientSolution object. + + The way that this object emulates the MATLAB mutli-dim. struct. array is with + the solution().steps attr. which is a list of solutionstep() instances + The process to recreate is as follows: + 1. Get instance of solution() with solution variable (the instance is made in make_results_subclasses) + 2. For each subgroup, create a solutionstep() class instance + 2a. Populate the instance with the key:value pairs + 2b. Append the instance to the solution().steps list + ''' + # step 1 + class_instance_name = name_of_struct + + # for some reason steps is not already a list + setattr(model_copy.results.__dict__[class_instance_name], 'steps', list()) + + steps = model_copy.results.__dict__[class_instance_name].steps + + # step 2 + layer = 1 + for subgroup in eval(group_location_in_file + ".groups.keys()"): + solutionstep_instance = solutionstep() + # step 2a + subgroup_location_in_file = group_location_in_file + ".groups['" + subgroup + "']" + for key in eval(subgroup_location_in_file + ".variables.keys()"): + value = eval(subgroup_location_in_file + ".variables['" + str(key) + "'][:]") + setattr(solutionstep_instance, key, value) + # step 2b + steps.append(solutionstep_instance) + if verbose: + print('Succesfully loaded layer ' + str(layer) + ' to results.' + str(class_instance_name) + ' struct.') + else: pass + layer += 1 + + if verbose: + print('Successfully recreated results structure ' + str(class_instance_name)) + + +def copy_variable_data_to_new_model(location_of_variable_in_file, location_of_variable_in_model, variable_name, NCData, verbose = False): + # as simple as navigating to the location_of_variable_in_model and setting it equal to the location_of_variable_in_file + # NetCDF4 has a property called "_FillValue" that sometimes saves empty lists, so we have to catch those + FillValue = -9223372036854775806 + try: + # results band-aid... + #print(str(location_of_variable_in_model + '.' + variable_name)) + if str(location_of_variable_in_model + '.' + variable_name) in ['results.solutionstep', 'results.solution', 'results.resultsdakota']: + pass + # qmu band-aid + elif 'qmu.statistics.method' in str(location_of_variable_in_model + '.' + variable_name): + pass + # handle any strings: + elif 'char' in eval(location_of_variable_in_file + '.dimensions[0]'): + setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, eval(location_of_variable_in_file + '[:][...].tobytes().decode()')) + # handle ndarrays + lists + elif len(eval(location_of_variable_in_file + '[:]'))>1: + # check for bool + try: # there is only one datatype assigned the attribute 'units' and that is bool, so anything else will go right to except + if eval(location_of_variable_in_file + '.units') == 'bool': + setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, np.array(eval(location_of_variable_in_file + '[:]'), dtype = bool)) + else: + setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, eval(location_of_variable_in_file + '[:]')) + except: + setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, eval(location_of_variable_in_file + '[:]')) + # catch everything else + else: + # check for FillValue. use try/except because try block will only work on datatypes like int64, float, single element lists/arrays ect and not nd-arrays/n-lists etc + try: + # this try block will only work on single ints/floats/doubles and will skip to the except block for all other cases + var_to_save = eval(location_of_variable_in_file + '[:][0]') # note the [0] on the end + if FillValue == var_to_save: + setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, []) + else: + if var_to_save.is_integer(): + setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, int(var_to_save)) + else: + # we have to convert numpy datatypes to native python types with .item() + setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, var_to_save.item()) + except: + setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, eval(location_of_variable_in_file + '[:]')) + except AttributeError: + copy_variable_data_to_new_model_dict(location_of_variable_in_file, location_of_variable_in_model, NCData, verbose=verbose) + + if verbose: + print('Successfully loaded ' + location_of_variable_in_model + '.' + variable_name + ' into model.') + + +def copy_variable_data_to_new_model_dict(location_of_variable_in_file, location_of_variable_in_model, NCData, verbose = False): + # as simple as navigating to the location_of_variable_in_model and setting it equal to the location_of_variable_in_file + + # NetCDF4 has a property called "_FillValue" that sometimes saves empty lists, so we have to catch those + FillValue = -9223372036854775806 + + # the key will be the last item in the location + key = ''.join(location_of_variable_in_model.split('.')[-1]) + + # update the location to point to the dict instead of the dict key + location_of_variable_in_model = '.'.join(location_of_variable_in_model.split('.')[:-1]) + + # verify we're working with a dict: + if isinstance(eval('model_copy.' + location_of_variable_in_model), OrderedDict): + dict_object = eval('model_copy.' + location_of_variable_in_model) + + # handle any strings: + if 'char' in eval(location_of_variable_in_file + '.dimensions[0]'): + data = eval(location_of_variable_in_file + '[:][...].tobytes().decode()') + dict_object.update({key: data}) + + # handle ndarrays + lists + elif len(eval(location_of_variable_in_file + '[:]'))>1: + # check for bool + try: # there is only one datatype assigned the attribute 'units' and that is bool, so anything else will go right to except + if eval(location_of_variable_in_file + '.units') == 'bool': + data = np.array(eval(location_of_variable_in_file + '[:]'), dtype = bool) + dict_object.update({key: data}) + else: + data = eval(location_of_variable_in_file + '[:]') + dict_object.update({key: data}) + except: + data = eval(location_of_variable_in_file + '[:]') + dict_object.update({key: data}) + # catch everything else + else: + # check for FillValue. use try/except because try block will only work on datatypes like int64, float, single element lists/arrays ect and not nd-arrays/n-lists etc + try: + # this try block will only work on single ints/floats/doubles and will skip to the except block for all other cases + if FillValue == eval(location_of_variable_in_file + '[:][0]'): + dict_object.update({key: []}) + else: + # we have to convert numpy datatypes to native python types with .item() + var_to_save = eval(location_of_variable_in_file + '[:][0]') # note the [0] on the end + dict_object.update({key: var_to_save.item()}) + except: + data = eval(location_of_variable_in_file + '[:]') + dict_object.update({key: data}) + else: + print(f"Unrecognized object was saved to NetCDF file and cannot be reconstructed: {location_of_variable_in_model}") Index: /issm/trunk-jpl/src/m/netcdf/write_netCDF.m =================================================================== --- /issm/trunk-jpl/src/m/netcdf/write_netCDF.m (revision 27891) +++ /issm/trunk-jpl/src/m/netcdf/write_netCDF.m (revision 27891) @@ -0,0 +1,685 @@ +%{ +Given a model, this set of functions will perform the following: + 1. Enter each nested class of the model. + 2. View each attribute of each nested class. + 3. Compare state of attribute in the model to an empty model class. + 4. If states are identical, pass. + 5. Otherwise, create nested groups named after class structure. + 6. Create variable named after class attribute and assign value to it. +%} + + +function write_netCDF(model_var, filename, varargin) + if nargin > 2 + verbose = true; + else + verbose = false; + end + if verbose + disp('MATLAB C2NetCDF4 v1.1.14'); + end + + % model_var = class object to be saved + % filename = path and name to save file under + + % Create a NetCDF file to write to + NetCDF = make_NetCDF(filename, verbose); + + % Create an instance of an empty model class to compare model_var against + empty_model = model(); + + % Walk through the model_var class and compare subclass states to empty_model + walk_through_model(model_var, empty_model, NetCDF, verbose); + + % in order to handle some subclasses in the results class, we have to utilize this band-aid + % there will likely be more band-aids added unless a class name library is created with all class names that might be added to a model + try + % if results had meaningful data, save the name of the subclass and class instance + netcdf.inqNcid(NetCDF,'results'); + results_subclasses_bandaid(model_var, NetCDF, verbose); + % otherwise, ignore + catch + end + + netcdf.close(NetCDF); + if verbose + disp('Model successfully saved as NetCDF4'); + end +end + + + +function NetCDF = make_NetCDF(filename, verbose) + % matlab can't handle input in the jupyter interface, so we just yell at the user to rename + % their file if needed + % If file already exists delete / rename it + if exist(filename, 'file') == 2 + fprintf('File %s already exists\n', filename); + disp('Please rename your file.') + return + + % If so, inquire for a new name or to delete the existing file + %newname = input('Give a new name or input "delete" to replace: ', 's'); + + %if strcmpi(newname, 'delete') + %delete filename; + %else + %fprintf('New file name is %s\n', newname); + %filename = newname; + %end + else + % Otherwise create the file and define it globally so other functions can call it + + NetCDF = netcdf.create(filename, 'NETCDF4'); + netcdf.putAtt(NetCDF, netcdf.getConstant('NC_GLOBAL'), 'history', ['Created ', datestr(now)]); + netcdf.defDim(NetCDF, 'Unlim', netcdf.getConstant('NC_UNLIMITED')); % unlimited dimension + netcdf.defDim(NetCDF, 'float', 1); % single integer dimension + netcdf.defDim(NetCDF, 'int', 1); % single float dimension + + if verbose + fprintf('Successfully created %s\n', filename); + end + + return + end +end + + +%{ + Since python uses subclass instances and MATLAB uses fields, we need to guess which subclass instance python will need + given the name of the sub-field in MATLAB. We make this guess based on the name of the MATLAB subfield that will contain + the name of the python subclass instance. For example, md.results.StressbalanceSolution is an subfield in MATLAB, + but a class instance of solution(). Notice that StressbalanceSolution contains the name "Solution" in it. This is what + we will save to the netCDF file for python to pick up. +%} + +function results_subclasses_bandaid(model_var, NetCDF, verbose) + + % The results class may have nested fields within it, so we need to record the name of + % the nested field as it appears in the model that we're trying to save + quality_control = {}; + + % Access the results subclass of model_var + results_var = model_var.results; + + % get the results group id so we can write to it + groupid = netcdf.inqNcid(NetCDF,'results'); + + % Loop through each class instance in results + class_instance_names = fieldnames(results_var); + + % we save lists of instances to the netcdf + solutions = {}; + solutionsteps = {}; + resultsdakotas = {}; + + for i = 1:numel(class_instance_names) + class_instance_name = class_instance_names{i}; + % there are often mutliple instances of the same class/struct so we have to number them + % Check to see if there is a solutionstep class instance + if contains(class_instance_name, 'solutionstep',IgnoreCase=true) + quality_control{end+1} = 1; + solutionsteps{end+1} = class_instance_name; + if verbose + disp('Successfully stored class python subclass instance: solutionstep') + end + end + + % Check to see if there is a solution class instance + if contains(class_instance_name, 'solution',IgnoreCase=true) + quality_control{end+1} = 1; + solutions{end+1} = class_instance_name; + if verbose + disp('Successfully stored class python subclass instance: solution') + end + end + + % Check to see if there is a resultsdakota class instance + if contains(class_instance_name, 'resultsdakota',IgnoreCase=true) + quality_control{end+1} = 1; + resultsdakotas{end+1} = class_instance_name; + if verbose + disp('Successfully stored class python subclass instance: resultsdakota') + end + end + end + if ~isempty(solutionsteps) + write_cell_with_strings('solutionstep', solutionsteps, groupid, NetCDF, verbose) + end + if ~isempty(solutions) + write_cell_with_strings('solution', solutions, groupid, NetCDF, verbose) + end + if ~isempty(resultsdakotas) + write_cell_with_strings('resultsdakota', resultsdakotas, groupid, NetCDF, verbose) + end + + + + % Check if all class instances were processed correctly + if numel(quality_control) ~= numel(class_instance_names) + disp('Error: The class instance within your model.results class is not currently supported by this application'); + disp(class(results_var.(class_instance_name))); + else + if verbose + disp('The results class was successfully stored on disk'); + end + end +end + + + +function walk_through_model(model_var, empty_model, NetCDF, verbose) + % Iterate over first layer of model_var attributes and assume this first layer is only classes fundamental to the model() class + % note that groups are the same as class instances/subfields in this context + groups = fieldnames(model_var); + for group = 1:numel(groups) + % now this new variable takes the form model.mesh , model.damage etc. + model_subclass = model_var.(groups{group}); + empty_model_subclass = empty_model.(groups{group}); + % Now we can recursively walk through the remaining subclasses + list_of_layers = {groups{group}}; + walk_through_subclasses(model_subclass, empty_model_subclass, list_of_layers, empty_model, NetCDF, verbose); + end +end + + +function walk_through_subclasses(model_subclass, empty_model_subclass, given_list_of_layers, empty_model, NetCDF, verbose) + % Recursivley iterate over each subclass' attributes and look for more subclasses and variables with relevant data + % model_subclass is an object (ie, md.mesh.elements) + % list_of_layers is a cell array of subclasses/attributes/fields so that we can copy the structure into netcdf (ie, {'mesh', 'elements'}) + % need to check if inversion or m1qn3inversion or taoinversion class + if numel(given_list_of_layers) == 1 + if strcmp(given_list_of_layers{1}, 'inversion') + create_group(model_subclass, given_list_of_layers, NetCDF, verbose); + check_inversion_class(model_subclass, NetCDF, verbose); + end + end + + % Use try/except since model_subclass is either a subclass/struct w/ props/fields or it's not, no unknown exceptions + try + % look for children - this is where the catch would be called + children = fieldnames(model_subclass); + + % if there are children, loop through them and see if we need to save any data + for child = 1:numel(children) + % record our current location + list_of_layers = given_list_of_layers; + current_child = children{child}; + list_of_layers{end+1} = current_child; + + % this is the value of the current location in the model (ie, md.mesh.elements) + location_of_child = model_subclass.(current_child); + + % if the empty model does not have this attribute, it's because it's new so we save it to netcdf + % there are 2 cases: the location is a struct, the location is a class + if isstruct(model_subclass) + % if the current field is a nested struct assume it has valuable data that needs to be saved + if isstruct(location_of_child) && any(size(location_of_child) > 1) + create_group(location_of_child, list_of_layers, NetCDF, verbose); + + % this would mean that the layer above the layer we're interested in is a struct, so + % we can navigate our empty model as such + elseif isfield(empty_model_subclass, current_child) + % the layer we're interested in does exist, we just need to compare states + location_of_child_in_empty_model = empty_model_subclass.(current_child); + + % if the current attribute is a numerical array assume it has valuable data that needs to be saved + if isnumeric(location_of_child) && logical(numel(location_of_child) > 1) + create_group(location_of_child, list_of_layers, NetCDF, verbose); + % if the attributes are identical we don't need to save anything + elseif (all(isnan(location_of_child)) && all(isnan(location_of_child_in_empty_model))) || isempty(setxor(location_of_child, location_of_child_in_empty_model)) + walk_through_subclasses(location_of_child, location_of_child_in_empty_model, list_of_layers, empty_model, NetCDF, verbose); + % if the attributes are not the same we need to save ours + else + % THE ORDER OF THESE LINES IS CRITICAL + walk_through_subclasses(location_of_child, location_of_child_in_empty_model, list_of_layers, empty_model, NetCDF, verbose); + create_group(location_of_child, list_of_layers, NetCDF, verbose); + end + % this would mean that the layer we're interested in is not fundamental to the model architecture + % and thus needs to be saved to the netcdf + else + walk_through_subclasses(location_of_child, empty_model_subclass, list_of_layers, empty_model, NetCDF, verbose); + create_group(location_of_child, list_of_layers, NetCDF, verbose); + end + % this would mean it's not a struct, and must be a class/subclass + % we now check the state of the class property + else + try + if isprop(empty_model_subclass, current_child) + % the layer we're interested in does exist, we just need to compare states + location_of_child_in_empty_model = empty_model_subclass.(current_child); + % if the current attribute is a numerical array assume it has valuable data that needs to be saved + if isnumeric(location_of_child) && logical(numel(location_of_child) > 1) + create_group(location_of_child, list_of_layers, NetCDF, verbose); + + elseif iscell(location_of_child) + % if the attributes are identical we don't need to save anything + if isempty(setxor(location_of_child, location_of_child_in_empty_model)) + % pass + else + % otherwise we need to save + walk_through_subclasses(location_of_child, empty_model_subclass, list_of_layers, empty_model, NetCDF, verbose); + create_group(location_of_child, list_of_layers, NetCDF, verbose); + end + elseif (all(isnan(location_of_child)) && all(isnan(location_of_child_in_empty_model))) + walk_through_subclasses(location_of_child, location_of_child_in_empty_model, list_of_layers, empty_model, NetCDF, verbose); + % if the attributes are not the same we need to save ours + else + % THE ORDER OF THESE LINES IS CRITICAL + walk_through_subclasses(location_of_child, location_of_child_in_empty_model, list_of_layers, empty_model, NetCDF, verbose); + create_group(location_of_child, list_of_layers, NetCDF, verbose); + end + else + walk_through_subclasses(location_of_child, empty_model_subclass, list_of_layers, empty_model, NetCDF, verbose); + create_group(location_of_child, list_of_layers, NetCDF, verbose); + end + catch + walk_through_subclasses(location_of_child, empty_model_subclass, list_of_layers, empty_model, NetCDF, verbose); + create_group(location_of_child, list_of_layers, NetCDF, verbose); + end + end + end + catch ME + % If the caught error is a fieldname error, it's just saying that a variable has no fields and thus can be ignored + if strcmp(ME.identifier, 'MATLAB:fieldnames:InvalidInput') + % do nothing + % this is if we come accross instances/subfields in our model that are not fundamental to the model class (ie, taoinversion) + elseif strcmp(ME.identifier, 'MATLAB:UndefinedFunction') + walk_through_subclasses(location_of_child, empty_model_subclass, given_list_of_layers, empty_model, NetCDF, verbose); + create_group(location_of_child, list_of_layers, NetCDF, verbose); + % If it's a different error, rethrow it to MATLAB's default error handling + else + disp(ME.identifier) + disp(given_list_of_layers) + rethrow(ME); + end + end +end + + +function create_group(location_of_child, list_of_layers, NetCDF, verbose) + % location_of_child is an object + % list_of_layers is a list like {'inversion', 'StressbalanceSolution','cost_functions_coefficients'} + % first we make the group at the highest level (ie, inversion) + group_name = list_of_layers{1}; + variable_name = list_of_layers{end}; + + % if the group is already made, get it's ID instead of creating it again + try % group hasn't been made + group = netcdf.defGrp(NetCDF, group_name); + catch % group was already made + group = netcdf.inqNcid(NetCDF, group_name); + end + + % if the data is nested, create nested groups to match class structure + if numel(list_of_layers) > 2 + % the string() method is really important here since matlab apparently can't handle the infinite complexity of a string without the string method. + for name = string(list_of_layers(2:end-1)) + % the group levels may have already been made + try % group hasn't been made + group = netcdf.defGrp(group, name); + catch % group was already made + group = netcdf.inqNcid(group, name); + end + end + end + % we may be dealing with an object + try + % if this line works, we're still dealing with a struct and need lower levels + if isempty(fieldnames(location_of_child)) + % do nothing + elseif isstruct(location_of_child) && any(size(location_of_child) > 1) + % we have a nested struct + copy_nested_struct(variable_name, location_of_child, group, NetCDF, verbose) + end + catch + % if that line doesn't work, it means we're dealing with raw data + % not a nested struct, so we can pass + create_var(variable_name, location_of_child, group, NetCDF, verbose); + end +end + + + +function copy_nested_struct(parent_struct_name, address_of_struct, group, NetCDF, verbose) + %{ + This function takes a struct of structs and saves them to netcdf. + + To do this, we get the number of dimensions (substructs) of the parent struct. + Next, we iterate through each substruct and record the data. + For each substruct, we create a subgroup of the main struct. + For each variable, we create dimensions that are assigned to each subgroup uniquely. + %} + + if verbose + disp("Beginning transfer of nested MATLAB struct to the NetCDF") + end + % make a new subgroup to contain all the others: + group = netcdf.defGrp(group, parent_struct_name); + + % make sure other systems can flag the nested struct type + dimID = netcdf.defDim(group, 'struct', 6); + string_var = netcdf.defVar(group, 'this_is_a_nested', "NC_CHAR", dimID); + uint_method=uint8('struct').'; + method_ID = char(uint_method); + netcdf.putVar(group, string_var, method_ID); + + % other systems know the name of the parent struct because it's covered by the results/qmu functions above + + % 'a' will always be 1 and is not useful to us + [a, no_of_dims] = size(address_of_struct); + + for substruct = 1:no_of_dims + % we start by making subgroups with nice names like "TransientSolution_substruct_44" + name_of_subgroup = ['1x', num2str(substruct)]; + subgroup = netcdf.defGrp(group, name_of_subgroup); + + % do some housekeeping to keep track of the current layer + current_substruct = address_of_struct(substruct); + substruct_fields = fieldnames(current_substruct)'; % transpose because matlab only interates over n x 1 arrays + + % now we need to iterate over each variable of the nested struct and save it to this new subgroup + for variable_name = string(substruct_fields) + address_of_child = current_substruct.(variable_name); + create_var(variable_name, address_of_child, subgroup, NetCDF, verbose); + end + end + if verbose + fprintf('Succesfully transferred nested MATLAB struct %s to the NetCDF\n', parent_struct_name) + end +end + + + +% ironically inversion does not have the same problem as results as inversion subfields +% are actually subclasses and not fields +function check_inversion_class(model_var, NetCDF, verbose) + + % Define a persistent variable to ensure this function is only run once + persistent executed; + % Check if the function has already been executed + if isempty(executed) + if verbose + disp('Deconstructing Inversion class instance') + end + % Need to make sure that we have the right inversion class: inversion, m1qn3inversion, taoinversion + groupid = netcdf.inqNcid(NetCDF,'inversion'); + + if isa(model_var, 'm1qn3inversion') + write_string_to_netcdf('inversion_class_name', 'm1qn3inversion', groupid, NetCDF, verbose); + if verbose + disp('Successfully saved inversion class instance m1qn3inversion') + end + elseif isa(model_var, 'taoinversion') + write_string_to_netcdf('inversion_class_name', 'taoinversion', groupid, NetCDF, verbose); + if verbose + disp('Successfully saved inversion class instance taoinversion') + end + else + write_string_to_netcdf('inversion_class_name', 'inversion', groupid, NetCDF, verbose); + if verbose + disp('Successfully saved inversion class instance inversion') + end + end + % Set the persistent variable to indicate that the function has been executed + executed = true; + end +end + + +function create_var(variable_name, address_of_child, group, NetCDF, verbose) + % There are lots of different variable types that we need to handle from the model class + + % get the dimensions we'll need + intdim = netcdf.inqDimID(NetCDF,'int'); + floatdim = netcdf.inqDimID(NetCDF,'float'); + unlimdim = netcdf.inqDimID(NetCDF,'Unlim'); + + % This first conditional statement will catch numeric arrays (matrices) of any dimension and save them + if any(size(address_of_child)>1) && ~iscellstr(address_of_child) && ~ischar(address_of_child) + write_numeric_array_to_netcdf(variable_name, address_of_child, group, NetCDF, verbose); + + % check if it's a string + elseif ischar(address_of_child) + write_string_to_netcdf(variable_name, address_of_child, group, NetCDF, verbose); + + % or an empty variable + elseif isempty(address_of_child) + variable = netcdf.defVar(group, variable_name, "NC_DOUBLE", intdim); + + % or a list of strings + elseif iscellstr(address_of_child) || iscell(address_of_child) && ischar(address_of_child{1}) + write_cell_with_strings(variable_name, address_of_child, group, NetCDF, verbose) + + % or an empty list + elseif iscell(address_of_child) && isempty(address_of_child) || isa(address_of_child, 'double') && isempty(address_of_child) + variable = netcdf.defVar(group, variable_name, "NC_INT", intdim); + netcdf.putVar(group,variable, -32767); + + % or a bool + elseif islogical(address_of_child) + % netcdf4 can't handle bool types like true/false so we convert all to int 1/0 and add an attribute named units with value 'bool' + variable = netcdf.defVar(group, variable_name, 'NC_SHORT', intdim); + netcdf.putVar(group,variable,int8(address_of_child)); + % make sure other systems can flag the bool type + netcdf.putAtt(group,variable,'units','bool'); + + % or a regular list + elseif iscell(address_of_child) + disp('made list w/ unlim dim') + variable = netcdf.defVar(group, variable_name, "NC_DOUBLE", unlimdim); + netcdf.putVar(group,variable,address_of_child); + + % or a float + elseif isfloat(address_of_child) && numel(address_of_child) == 1 + variable = netcdf.defVar(group, variable_name, "NC_DOUBLE", floatdim); + netcdf.putVar(group,variable,address_of_child); + + % or a int + elseif mod(address_of_child,1) == 0 || isinteger(address_of_child) && numel(address_of_child) == 1 + variable = netcdf.defVar(group, variable_name, "NC_SHORT", intdim); + netcdf.putVar(group,variable,address_of_child); + + % anything else... (will likely need to add more cases; ie dict) + else + try + variable = netcdf.defVar(group, variable_name, "NC_DOUBLE", unlimdim); + netcdf.putVar(group,variable,address_of_child); + catch ME + disp(ME.message); + disp(['Datatype given: ', class(address_of_child)]); + end + end + if verbose + fprintf('Successfully transferred data from %s to the NetCDF\n', variable_name); + end +end + + +function write_cell_with_strings(variable_name, address_of_child, group, NetCDF, verbose) + %{ + Write cell array (ie {'one' 'two' 'three'}) to netcdf + %} + + if isempty(address_of_child) + % if the char array is empty, save an empty char + name_of_dimension = ['char', num2str(0)]; + try + dimID = netcdf.defDim(group, name_of_dimension, 0); + catch + dimID = netcdf.inqDimID(group, name_of_dimension); + end + % Now we can make a variable in this dimension: + string_var = netcdf.defVar(group, variable_name, "NC_CHAR", [dimID]); + % we leave empty now + else + % covert data to char array + method_ID = char(address_of_child); + + % make dimensions + [rows, cols] = size(method_ID); + + IDDim1 = netcdf.defDim(group,'cols',cols); + IDDim2 = netcdf.defDim(group,'rows',rows); + + % create the variable slot + IDVarId = netcdf.defVar(group,variable_name,'NC_CHAR', [IDDim1 IDDim2]); + + % save the variable + netcdf.putVar(group, IDVarId, method_ID'); %transpose + + % tell other platforms that this is a cell of strings + netcdf.putAtt(group, IDVarId, 'type_is','cell_array_of_strings'); + end +end + + +function write_string_to_netcdf(variable_name, address_of_child, group, NetCDF, verbose) + % netcdf and strings don't get along.. we have to do it 'custom': + + the_string_to_save = address_of_child; + + if isempty(the_string_to_save) + % if the char array is empty, save an empty char + name_of_dimension = ['char', num2str(0)]; + try + dimID = netcdf.defDim(group, name_of_dimension, 0); + catch + dimID = netcdf.inqDimID(group, name_of_dimension); + end + % Now we can make a variable in this dimension: + string_var = netcdf.defVar(group, variable_name, "NC_CHAR", [dimID]); + % we leave empty now + else + % convert string to + uint_method=uint8(the_string_to_save).'; + method_ID = char(uint_method); + length_of_the_string = numel(method_ID); + + % Convert the string to character data using string array + %str_out = char(the_string_to_save) + + % Determine the length of the string + %length_of_the_string = numel(str_out) + + % Check if the dimension already exists, and if not, create it + name_of_dimension = ['char', num2str(length_of_the_string)]; + try + dimID = netcdf.defDim(group, name_of_dimension, length_of_the_string); + catch + dimID = netcdf.inqDimID(group, name_of_dimension); + end + % Now we can make a variable in this dimension: + string_var = netcdf.defVar(group, variable_name, "NC_CHAR", [dimID]); + % Finally, we can write the variable (always transpose for matlab): + netcdf.putVar(group, string_var, method_ID); + end + + if verbose + disp(['Successfully transferred data from ', variable_name, ' to the NetCDF']); + end +end + + +function write_numeric_array_to_netcdf(variable_name, address_of_child, group, NetCDF, verbose) + % get the dimensions we'll need + intdim = netcdf.inqDimID(NetCDF,'int'); + floatdim = netcdf.inqDimID(NetCDF,'float'); + unlimdim = netcdf.inqDimID(NetCDF,'Unlim'); + + typeis = class(address_of_child); + + if isa(typeis, 'logical') + % because matlab transposes all data into and out of netcdf and because we want cross-platform-compat + % we need to transpose data before it goes into netcdf + data = address_of_child.'; + + % make the dimensions + dimensions = []; + for dimension = size(data) + dim_name = ['dim',int2str(dimension)]; + % if the dimension already exists we can't have a duplicate + try + dimID = netcdf.defDim(group, dim_name, dimension); + catch + dimID = netcdf.inqDimID(group, dim_name); + end + % record the dimension for the variable + dimensions(end+1) = dimID; + end + + % write the variable + netcdf.putVar(group,variable,data); + + % make sure other systems can flag the bool type + netcdf.putAtt(group,variable,'units','bool'); + + % handle all other datatypes here + else + % sometimes an array has just 1 element in it, we account for those cases here: + if numel(address_of_child) == 1 + if isinteger(address_of_child) + variable = netcdf.defVar(group, variable_name, "NC_SHORT", intdim); + netcdf.putVar(group,variable,address_of_child); + elseif isa(address_of_child, 'double') || isa(address_of_child, 'float') + variable = netcdf.defVar(group, variable_name, "NC_DOUBLE", floatdim); + netcdf.putVar(group,variable,address_of_child); + else + disp('Encountered single datatype that was not float64 or int64, saving under unlimited dimension, may cause errors.') + variable = netcdf.defVar(group, variable_name, "NC_DOUBLE", unlimdim); + netcdf.putVar(group,variable,address_of_child); + end + % this is in case of lists so that python doesn't get a (nx1) numpy array and instead gets an n-element list + elseif any(size(address_of_child)==1) + % because matlab transposes all data into and out of netcdf and because we want cross-platform-compat + % we need to transpose data before it goes into netcdf + data = address_of_child.'; + + % make the dimensions + dimensions = []; + for dimension = size(data) + if dimension ~= 1 + dim_name = ['dim',int2str(dimension)]; + % if the dimension already exists we can't have a duplicate + try + dimID = netcdf.defDim(group, dim_name, dimension); + catch + dimID = netcdf.inqDimID(group, dim_name); + end + % record the dimension for the variable + dimensions(end+1) = dimID; + end + end + % create the variable + variable = netcdf.defVar(group, variable_name, "NC_DOUBLE",dimensions); + + % write the variable + netcdf.putVar(group,variable,data); + + % This catches all remaining arrays: + else + % because matlab transposes all data into and out of netcdf and because we want cross-platform-compat + % we need to transpose data before it goes into netcdf + data = address_of_child.'; + + % make the dimensions + dimensions = []; + for dimension = size(data) + dim_name = ['dim',int2str(dimension)]; + % if the dimension already exists we can't have a duplicate + try + dimID = netcdf.defDim(group, dim_name, dimension); + catch + dimID = netcdf.inqDimID(group, dim_name); + end + % record the dimension for the variable + dimensions(end+1) = dimID; + end + % create the variable + variable = netcdf.defVar(group, variable_name, "NC_DOUBLE",dimensions); + + % write the variable + netcdf.putVar(group,variable,data); + end + end +end Index: /issm/trunk-jpl/src/m/netcdf/write_netCDF.py =================================================================== --- /issm/trunk-jpl/src/m/netcdf/write_netCDF.py (revision 27891) +++ /issm/trunk-jpl/src/m/netcdf/write_netCDF.py (revision 27891) @@ -0,0 +1,496 @@ +# imports +import netCDF4 +from netCDF4 import Dataset +import numpy as np +import numpy.ma as ma +import time +import os +from model import * +from results import * +from m1qn3inversion import m1qn3inversion +from taoinversion import taoinversion +#import OrderedStruct + + +''' +Given a md, this set of functions will perform the following: + 1. Enter each nested class of the md. + 2. View each attribute of each nested class. + 3. Compare state of attribute in the model to an empty model class. + 4. If states are identical, pass. + 5. Otherwise, create nested groups named after class structure. + 6. Create variable named after class attribute and assign value to it. +''' + + +def write_netCDF(md, filename: str, verbose = False): + if verbose: + print('Python C2NetCDF4 v1.1.14') + else: pass + ''' + md = model() class instance to be saved + filename = path and name to save file under + verbose = T/F muted or show log statements. Naturally muted + ''' + + # Create a NetCDF file to write to + NetCDF = make_NetCDF(filename, verbose) + + # Create an instance of an empty md class to compare md_var against + empty_model = model() + + # Walk through the md class and compare subclass states to empty_model + walk_through_model(md, empty_model, NetCDF, verbose) + + # in order to handle some subclasses in the results class, we have to utilize this band-aid + # there will likely be more band-aids added unless a class name library is created with all class names that might be added to a md + try: + # if results has meaningful data, save the name of the subclass and class instance + NetCDF.groups['results'] + results_subclasses_bandaid(md, NetCDF, verbose) + # otherwise, ignore + except KeyError: + pass + + NetCDF.close() + if verbose: + print('Model successfully saved as NetCDF4') + else: pass + + +def results_subclasses_bandaid(md, NetCDF, verbose = False): + # since the results class may have nested classes within it, we need to record the name of the + # nested class instance variable as it appears in the md that we're trying to save + quality_control = [] + + # we save lists of instances to the netcdf + solutions = [] + solutionsteps = [] + resultsdakotas = [] + + for class_instance_name in md.results.__dict__.keys(): + if verbose: + print(class_instance_name) + # for each class instance in results, see which class its from and record that info in the netcdf to recreate structure later + # check to see if there is a solutionstep class instance + if isinstance(md.results.__dict__[class_instance_name],solutionstep): + quality_control.append(1) + solutionsteps.append(class_instance_name) + + # check to see if there is a solution class instance + if isinstance(md.results.__dict__[class_instance_name],solution): + quality_control.append(1) + solutions.append(class_instance_name) + + # check to see if there is a resultsdakota class instance + if isinstance(md.results.__dict__[class_instance_name],resultsdakota): + quality_control.append(1) + resultsdakotas.append(class_instance_name) + + if solutionsteps != []: + write_string_to_netcdf(variable_name=str('solutionstep'), address_of_child=solutionsteps, group=NetCDF.groups['results'], list=True, NetCDF=NetCDF, verbose=verbose) + + if solutions != []: + write_string_to_netcdf(variable_name=str('solution'), address_of_child=solutions, group=NetCDF.groups['results'], list=True, NetCDF=NetCDF, verbose=verbose) + + if resultsdakotas != []: + write_string_to_netcdf(variable_name=str('resultsdakota'), address_of_child=resultsdakotas, group=NetCDF.groups['results'], list=True, NetCDF=NetCDF, verbose=verbose) + + + if len(quality_control) != len(md.results.__dict__.keys()): + print('Error: The class instance within your md.results class is not currently supported by this application') + print(type(md.results.__dict__[class_instance_name])) + else: + if verbose: + print('The results class was successfully stored on disk') + else: pass + + +def make_NetCDF(filename: str, verbose = False): + # If file already exists delete / rename it + if os.path.exists(filename): + print('File {} allready exist'.format(filename)) + + # If so, inqure for a new name or to do delete the existing file + newname = input('Give a new name or "delete" to replace: ') + + if newname == 'delete': + os.remove(filename) + else: + print(('New file name is {}'.format(newname))) + filename = newname + else: + # Otherwise create the file and define it globally so other functions can call it + NetCDF = Dataset(filename, 'w', format='NETCDF4') + NetCDF.history = 'Created ' + time.ctime(time.time()) + NetCDF.createDimension('Unlim', None) # unlimited dimension + NetCDF.createDimension('float', 1) # single integer dimension + NetCDF.createDimension('int', 1) # single float dimension + + if verbose: + print('Successfully created ' + filename) + + return NetCDF + + +def walk_through_model(md, empty_model, NetCDF, verbose= False): + # Iterate over first layer of md_var attributes and assume this first layer is only classes + for group in md.__dict__.keys(): + address = md.__dict__[group] + empty_address = empty_model.__dict__[group] + # we need to record the layers of the md so we can save them to the netcdf file + layers = [group] + + # Recursively walk through subclasses + walk_through_subclasses(address, empty_address, layers, NetCDF, empty_model, verbose) + + +def walk_through_subclasses(address, empty_address, layers: list, NetCDF, empty_model, verbose = False): + # See if we have an object with keys or a not + try: + address.__dict__.keys() + is_object = True + except: is_object = False # this is not an object with keys + + if is_object: + # enter the subclass, see if it has nested classes and/or attributes + # then compare attributes between mds and write to netCDF if they differ + # if subclass found, walk through it and repeat + for child in address.__dict__.keys(): + # record the current location + current_layer = layers.copy() + current_layer.append(child) + + # navigate to child in each md + address_of_child = address.__dict__[child] + + # if the current object is a results. object and has the steps attr it needs special treatment + if isinstance(address_of_child, solution) and len(address_of_child.steps) != 0: + create_group(address_of_child, current_layer, is_struct = True, NetCDF=NetCDF, verbose = verbose) + + # if the variable is an array, assume it has relevant data (this is because the next line cannot evaluate "==" with an array) + elif isinstance(address_of_child, np.ndarray): + create_group(address_of_child, current_layer, is_struct = False, NetCDF=NetCDF, verbose = verbose) + + # see if the child exists in the empty md. If not, record it in the netcdf + else: + try: + address_of_child_in_empty_class = empty_address.__dict__[child] + # if that line worked, we can see how the mds' attributes at this layer compare: + + # if the attributes are identical we don't need to save anything + if address_of_child == address_of_child_in_empty_class: + walk_through_subclasses(address_of_child, address_of_child_in_empty_class, current_layer, NetCDF, empty_model, verbose) + + # If it has been modified, record it in the NetCDF file + else: + create_group(address_of_child, current_layer, is_struct = False, NetCDF=NetCDF, verbose = verbose) + walk_through_subclasses(address_of_child, address_of_child_in_empty_class, current_layer, NetCDF, empty_model, verbose) + + except KeyError: # record in netcdf and continue to walk thru md + walk_through_subclasses(address_of_child, empty_address, current_layer, NetCDF, empty_model, verbose) + create_group(address_of_child, current_layer, is_struct = False, NetCDF=NetCDF, verbose = verbose) + else: pass + + +def create_group(address_of_child, layers, is_struct = False, NetCDF=None, verbose = False): + + # Handle the first layer of the group(s) + group_name = layers[0] + try: + group = NetCDF.createGroup(str(group_name)) + except: + group = NetCDF.groups[str(group_name)] + + # need to check if inversion or m1qn3inversion class + if group_name == 'inversion': + check_inversion_class(address_of_child, NetCDF, verbose) + else: pass + + # if the data is nested, create nested groups to match class structure + if len(layers) > 2: + for name in layers[1:-1]: + try: + group = group.createGroup(str(name)) + except: + group = NetCDF.groups[str(name)] + else: pass + + # Lastly, handle the variable(s) + if is_struct: + parent_struct_name = layers[-1] + copy_nested_results_struct(parent_struct_name, address_of_child, group, NetCDF, verbose) + + else: + variable_name = layers[-1] + create_var(variable_name, address_of_child, group, NetCDF, verbose) + + +def singleton(func): + """ + A decorator to ensure a function is only executed once. + """ + def wrapper(*args, **kwargs): + if not wrapper.has_run: + wrapper.result = func(*args, **kwargs) + wrapper.has_run = True + return wrapper.result + wrapper.has_run = False + wrapper.result = None + return wrapper + + +@singleton +def check_inversion_class(address_of_child, NetCDF, verbose = False): + # need to make sure that we have the right inversion class: inversion, m1qn3inversion, taoinversion + if isinstance(address_of_child, m1qn3inversion): + write_string_to_netcdf(variable_name=str('inversion_class_name'), address_of_child=str('m1qn3inversion'), group=NetCDF.groups['inversion'], NetCDF=NetCDF, verbose = verbose) + if verbose: + print('Successfully saved inversion class instance ' + 'm1qn3inversion') + elif isinstance(address_of_child, taoinversion): + write_string_to_netcdf(variable_name=str('inversion_class_name'), address_of_child=str('taoinversion'), group=NetCDF.groups['inversion'], NetCDF=NetCDF, verbose = verbose) + if verbose: + print('Successfully saved inversion class instance ' + 'taoinversion') + else: + write_string_to_netcdf(variable_name=str('inversion_class_name'), address_of_child=str('inversion'), group=NetCDF.groups['inversion'], NetCDF=NetCDF, verbose = verbose) + if verbose: + print('Successfully saved inversion class instance ' + 'inversion') + + +def copy_nested_results_struct(parent_struct_name, address_of_struct, group, NetCDF, verbose = False): + ''' + This function takes a solution class instance and saves the solutionstep instances from .steps to the netcdf. + + To do this, we get the number of dimensions (substructs) of the parent struct. + Next, we iterate through each substruct and record the data. + For each substruct, we create a subgroup of the main struct. + For each variable, we create dimensions that are assigned to each subgroup uniquely. + ''' + if verbose: + print("Beginning transfer of nested MATLAB struct to the NetCDF") + + # make a new subgroup to contain all the others: + group = group.createGroup(str(parent_struct_name)) + + # make sure other systems can flag the nested struct type + write_string_to_netcdf('this_is_a_nested', 'struct', group, list=False, NetCDF=NetCDF, verbose = verbose) + + # other systems know the name of the parent struct because it's covered by the results/qmu functions above + no_of_dims = len(address_of_struct) + for substruct in range(0, no_of_dims): + # we start by making subgroups with nice names like "1x4" + name_of_subgroup = '1x' + str(substruct) + subgroup = group.createGroup(str(name_of_subgroup)) + + # do some housekeeping to keep track of the current layer + current_substruct = address_of_struct[substruct] + substruct_fields = current_substruct.__dict__.keys() + + # now we need to iterate over each variable of the nested struct and save it to this new subgroup + for variable in substruct_fields: + address_of_child = current_substruct.__dict__[variable] + create_var(variable, address_of_child, subgroup, NetCDF, verbose = verbose) + + if verbose: + print(f'Successfully transferred struct {parent_struct_name} to the NetCDF\n') + + +def create_var(variable_name, address_of_child, group, NetCDF, verbose = False): + # There are lots of different variable types that we need to handle from the md class + + # This first conditional statement will catch numpy arrays of any dimension and save them + if isinstance(address_of_child, np.ndarray): + write_numpy_array_to_netcdf(variable_name, address_of_child, group, NetCDF, verbose=verbose) + + # check if it's an int + elif isinstance(address_of_child, int) or isinstance(address_of_child, np.integer): + variable = group.createVariable(variable_name, int, ('int',)) + variable[:] = address_of_child + + # or a float + elif isinstance(address_of_child, float) or isinstance(address_of_child, np.floating): + variable = group.createVariable(variable_name, float, ('float',)) + variable[:] = address_of_child + + # or a string + elif isinstance(address_of_child, str): + write_string_to_netcdf(variable_name, address_of_child, group, NetCDF, verbose=verbose) + + #or a bool + elif isinstance(address_of_child, bool) or isinstance(address_of_child, np.bool_): + # netcdf4 can't handle bool types like True/False so we convert all to int 1/0 and add an attribute named units with value 'bool' + variable = group.createVariable(variable_name, int, ('int',)) + variable[:] = int(address_of_child) + variable.units = "bool" + + # or an empty list + elif isinstance(address_of_child, list) and len(address_of_child)==0: + variable = group.createVariable(variable_name, int, ('int',)) + + # or a list of strings -- this needs work as it can only handle a list of 1 string + elif isinstance(address_of_child,list) and isinstance(address_of_child[0],str): + for string in address_of_child: + write_string_to_netcdf(variable_name, string, group, list=True, NetCDF=NetCDF, verbose=verbose) + + # or a regular list + elif isinstance(address_of_child, list): + variable = group.createVariable(variable_name, type(address_of_child[0]), ('Unlim',)) + variable[:] = address_of_child + + # anything else... (will likely need to add more cases; ie helpers.OrderedStruct) + else: + try: + variable = group.createVariable(variable_name, type(address_of_child), ('Unlim',)) + variable[:] = address_of_child + print(f'Unrecognized variable was saved {variable_name}') + except TypeError: pass # this would mean that we have an object, so we just let this continue to feed thru the recursive function above + except Exception as e: + print(f'There was error with {variable_name} in {group}') + print("The error message is:") + print(e) + print('Datatype given: ' + str(type(address_of_child))) + + if verbose: + print(f'Successfully transferred data from {variable_name} to the NetCDF') + + +def write_string_to_netcdf(variable_name, address_of_child, group, list=False, NetCDF=None, verbose = False): + # netcdf and strings dont get along.. we have to do it 'custom': + # if we hand it an address we need to do it this way: + if list == True: + """ + Save a list of strings to a NetCDF file. + + Convert a list of strings to a numpy.char_array with utf-8 encoded elements + and size rows x cols with each row the same # of cols and save to NetCDF + as char array. + """ + try: + strings = address_of_child + # get dims of array to save + rows = len(strings) + cols = len(max(strings, key = len)) + + # Define dimensions for the strings + rows_name = 'rows' + str(rows) + cols_name = 'cols' + str(cols) + try: + group.createDimension(rows_name, rows) + except: pass + + try: + group.createDimension(cols_name, cols) + except: pass + + # Create a variable to store the strings + string_var = group.createVariable(str(variable_name), 'S1', (rows_name, cols_name)) + + # break the list into a list of lists of words with the same length as the longest word: + # make words same sizes by adding spaces + modded_strings = [word + ' ' * (len(max(strings, key=len)) - len(word)) for word in strings] + # encoded words into list of encoded lists + new_list = [[s.encode('utf-8') for s in word] for word in modded_strings] + + # make numpy char array with dims rows x cols + arr = np.chararray((rows, cols)) + + # fill array with list of encoded lists + for i in range(len(new_list)): + arr[i] = new_list[i] + + # save array to netcdf file + string_var[:] = arr + + if verbose: + print(f'Saved {len(modded_strings)} strings to {variable_name}') + + except Exception as e: + print(f'Error: {e}') + + else: + the_string_to_save = address_of_child + length_of_the_string = len(the_string_to_save) + numpy_datatype = 'S' + str(length_of_the_string) + str_out = netCDF4.stringtochar(np.array([the_string_to_save], dtype=numpy_datatype)) + + # we'll need to make a new dimension for the string if it doesn't already exist + name_of_dimension = 'char' + str(length_of_the_string) + try: + group.createDimension(name_of_dimension, length_of_the_string) + except: pass + # this is another band-aid to the results sub classes... + try: + # now we can make a variable in this dimension: + string = group.createVariable(variable_name, 'S1', (name_of_dimension)) + #finally we can write the variable: + string[:] = str_out + #except RuntimeError: pass + except Exception as e: + print(f'There was an error saving a string from {variable_name}') + print(e) + + +def write_numpy_array_to_netcdf(variable_name, address_of_child, group, NetCDF, verbose = False): + # to make a nested array in netCDF, we have to get the dimensions of the array, + # create corresponding dimensions in the netCDF file, then we can make a variable + # in the netCDF with dimensions identical to those in the original array + + # start by getting the data type at the lowest level in the array: + typeis = address_of_child.dtype + + # catch boolean arrays here + if typeis == bool: + # sometimes an array has just 1 element in it, we account for those cases here: + if len(address_of_child) == 1: + variable = group.createVariable(variable_name, int, ('int',)) + variable[:] = int(address_of_child) + variable.units = "bool" + else: + # make the dimensions + dimensions = [] + for dimension in np.shape(address_of_child): + dimensions.append(str('dim' + str(dimension))) + # if the dimension already exists we can't have a duplicate + try: + group.createDimension(str('dim' + str(dimension)), dimension) + except: pass # this would mean that the dimension already exists + + # create the variable: + variable = group.createVariable(variable_name, int, tuple(dimensions)) + # write the variable: + variable[:] = address_of_child.astype(int) + variable.units = "bool" + + # handle all other datatypes here + else: + # sometimes an array has just 1 element in it, we account for those cases here: + if len(address_of_child) == 1: + if typeis is np.dtype('float64'): + variable = group.createVariable(variable_name, typeis, ('float',)) + variable[:] = address_of_child[0] + elif typeis is np.dtype('int64'): + variable = group.createVariable(variable_name, typeis, ('int',)) + variable[:] = address_of_child[0] + else: + print(f'Encountered single datatype from {variable_name} that was not float64 or int64, saving under unlimited dimension, may cause errors.') + variable = group.createVariable(variable_name, typeis, ('Unlim',)) + variable[:] = address_of_child[0] + + # This catches all arrays/lists: + else: + # make the dimensions + dimensions = [] + for dimension in np.shape(address_of_child): + dimensions.append(str('dim' + str(dimension))) + # if the dimension already exists we can't have a duplicate + try: + group.createDimension(str('dim' + str(dimension)), dimension) + except: pass # this would mean that the dimension already exists + + # create the variable: + variable = group.createVariable(variable_name, typeis, tuple(dimensions)) + + # write the variable: + variable[:] = address_of_child + +