Index: /issm/trunk/src/m/contrib/musselman/write_netCDF.m
===================================================================
--- /issm/trunk/src/m/contrib/musselman/write_netCDF.m	(revision 27868)
+++ /issm/trunk/src/m/contrib/musselman/write_netCDF.m	(revision 27868)
@@ -0,0 +1,465 @@
+%{
+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)
+    disp('MATLAB C2NetCDF4 v1.1.12');
+    % model_var = class object to be saved
+    % filename = path and name to save file under
+    
+    % Create a NetCDF file to write to
+    make_NetCDF(filename);
+    
+    % Create an instance of an empty model class to compare model_var against
+    global empty_model;
+    empty_model = model();
+
+    % Walk through the model_var class and compare subclass states to empty_model
+    walk_through_model(model_var);
+
+    % 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
+    global NetCDF;
+    try
+        % if results had meaningful data, save the name of the subclass and class instance
+        NetCDF.groups('results');
+        results_subclasses_bandaid(model_var);
+        % otherwise, ignore
+    catch
+    end
+    
+    netcdf.close(NetCDF);
+    disp('Model successfully saved as NetCDF4');
+end
+
+
+function make_NetCDF(filename)
+    % If file already exists delete / rename it
+    if exist(filename, 'file') == 2
+        fprintf('File %s already exists\n', filename);
+    
+        % 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
+    end
+    % Otherwise create the file and define it globally so other functions can call it
+    global NetCDF;
+    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
+    
+
+    fprintf('Successfully created %s\n', filename);
+end
+
+
+function results_subclasses_bandaid(model_var)
+    % 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 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);
+    
+    for i = 1:numel(class_instance_names)
+        class_instance_name = class_instance_names{i};
+        
+        % Check to 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 isa(results_var.(class_instance_name), 'solutionstep')
+            quality_control{end+1} = 1;
+            write_string_to_netcdf('variable_name', 'solutionstep', 'address_of_child', class_instance_name, 'group', groupid);
+        end
+        
+        % Check to see if there is a solution class instance
+        if isa(results_var.(class_instance_name), 'solution')
+            quality_control{end+1} = 1;
+            write_string_to_netcdf('variable_name', 'solution', 'address_of_child', class_instance_name, 'group', groupid);
+        end
+        
+        % Check to see if there is a resultsdakota class instance
+        if isa(results_var.(class_instance_name), 'resultsdakota')
+            quality_control{end+1} = 1;
+            write_string_to_netcdf('variable_name', 'resultsdakota', 'address_of_child', class_instance_name, 'group', groupid);
+        end
+    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
+        disp('The results class was successfully stored on disk');
+    end
+end
+
+
+
+function walk_through_model(model_var)
+    global empty_model;
+    % 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);
+    end
+end
+        
+
+function walk_through_subclasses(model_subclass, empty_model_subclass, given_list_of_layers)
+    % 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'})
+    % Use try/except since model_subclass is either a subclass w/ fields or it's not, no unknown exceptions
+    try 
+        % look for children - this is where the catch would be called
+        children = fieldnames(model_subclass);
+        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
+            try
+                location_of_child_in_empty_model = empty_model_subclass.(current_child);
+            catch
+                % empty model didn't have the attr, so we save it to netcdf
+                create_group(location_of_child, list_of_layers);
+            end
+            
+            % 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);
+                % i don't think I need this line but it's in my python code...
+                % walk_through_subclasses(location_of_child)
+            % 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);
+            % if the attributes are not the same we need to save ours
+            else
+                create_group(location_of_child, list_of_layers);
+                walk_through_subclasses(location_of_child, location_of_child_in_empty_model, list_of_layers);
+            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')
+        % 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')
+        % If it's a different error, rethrow it to MATLAB's default error handling
+        else
+            rethrow(ME);
+        end
+    end
+end 
+        
+
+function create_group(location_of_child, list_of_layers)
+    global NetCDF;
+    % 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};
+    
+    % 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
+    % need to check if inversion or m1qn3inversion or taoinversion class
+    if strcmp(group_name, 'inversion')
+        check_inversion_class(location_of_child);
+    end
+    
+    % if the data is nested, create nested groups to match class structure
+    if numel(list_of_layers) > 2
+        for name = list_of_layers{2:end-1}
+            % again, 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
+
+    % Lastly, handle the variable(s)
+    variable_name = list_of_layers{end};
+    create_var(variable_name, location_of_child, group);
+    
+end
+
+
+
+function check_inversion_class(model_var)
+    global NetCDF;
+    % 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)
+        % 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);
+            disp('Successfully saved inversion class instance m1qn3inversion');
+        elseif isa(model_var, 'taoinversion')
+            write_string_to_netcdf('inversion_class_name', 'taoinversion', groupid);
+            disp('Successfully saved inversion class instance taoinversion');
+        else
+            write_string_to_netcdf('inversion_class_name', 'inversion', groupid);
+            disp('Successfully saved inversion class instance inversion');
+        end
+        % Set the persistent variable to indicate that the function has been executed
+        executed = true;
+    else
+        % Code to skip execution if the function has already run
+    end
+
+
+end
+
+
+
+function create_var(variable_name, address_of_child, group)
+    global NetCDF;
+    % 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) && ~ischar(address_of_child)
+        write_numeric_array_to_netcdf(variable_name, address_of_child, group);
+    
+    % check if it's an int
+    elseif mod(address_of_child,1) == 0 || isinteger(address_of_child) && numel(address_of_child) == 1
+        variable = netcdf.defVar(group, variable_name, "NC_INT", intdim);
+        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 string
+    elseif ischar(address_of_child)
+        write_string_to_netcdf(variable_name, address_of_child, group);
+
+    % 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', 'int');
+        netcdf.putVar(group,variable,logical(address_of_child));
+        % make sure other systems can flag the bool type
+        netcdf.putAtt(group,variable,'units','bool');
+        
+    % 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_SHORT", intdim);
+
+    % or a list of strings -- this needs work as it can only handle a list of 1 string
+    elseif iscell(address_of_child) && ischar(address_of_child{1})
+        for i = 1:numel(address_of_child)
+            write_string_to_netcdf(variable_name, address_of_child{i}, group, list = true);
+        end
+
+    % 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);
+
+    % 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);
+            disp('Used Unlim Dim');
+        catch ME
+            disp(ME.message);
+            disp(['Datatype given: ', class(address_of_child)]);
+        end
+    end
+
+    disp(['Successfully transferred data from ', variable_name, ' to the NetCDF']);
+end
+
+
+
+
+function write_string_to_netcdf(variable_name, address_of_child, group, list)
+    % netcdf and strings don't get along.. we have to do it 'custom':
+    global NetCDF;
+    % If 'list' is not provided, assume it is false
+    if nargin < 4
+        list = false;
+    end
+    the_string_to_save = address_of_child;
+
+    % 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:
+    netcdf.putVar(group, string_var, method_ID);
+
+    disp(['Successfully transferred data from ', variable_name, ' to the NetCDF']);
+end
+
+
+
+function write_numeric_array_to_netcdf(variable_name, address_of_child, group)
+    global NetCDF;
+    % 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
+