source: issm/trunk/src/m/consistency/checkfield.py@ 25836

import numpy as np
import os
from re import findall, split
from pairoptions import pairoptions
from operator import attrgetter
import MatlabFuncs as m


def checkfield(md, *args):
    """CHECKFIELD - check field consistency

    Used to check model consistency
        
    Requires:
        'field' or 'fieldname' option. If 'fieldname' is provided, it will 
        retrieve it from the model md. (md.(fieldname))
        If 'field' is provided, it will assume the argument following 'field' 
        is a numeric array.

    Available options:
        - NaN: 1 if check that there is no NaN
        - size: [lines cols], NaN for non checked dimensions, or 'universal' 
        for any input type (nodal, element, time series, etc)
        - > :  greater than provided value
        - >= : greater or equal to provided value
        - < :  smallerthan provided value
        - <=: smaller or equal to provided value
        - < vec:  smallerthan provided values on each vertex
        - timeseries: 1 if check time series consistency (size and time)
        - values: list of acceptable values
        - numel: list of acceptable number of elements
        - cell: 1 if check that is cell
        - empty: 1 if check that non empty
        - message: overloaded error message

    Usage:
        md = checkfield(md, fieldname, options)
    """

    #get options
    options = pairoptions(*args)

    #get field from model
    if options.exist('field'):
        field = options.getfieldvalue('field')
        fieldname = options.getfieldvalue('fieldname', 'no fieldname')
    else:
        fieldname = options.getfieldvalue('fieldname')
        fieldprefix = split(r'\[(.*?)\]', fieldname)[0]
        fieldindexes = findall(r'\[(.*?)\]', fieldname)
        field = attrgetter(fieldprefix)(md)
        for index in fieldindexes:
            try:
                field = field[index.strip("\'")]
            except TypeError:
                field = field[int(index)]  #looking for an index and not a key

    # that works for py2
    #        exec("field = md.{}".format(fieldname))
    #        exec("field = md.{}".format(fieldname), namespace)

    if isinstance(field, (bool, int, float)):
        field = np.array([field])

    #check empty
    if options.exist('empty'):
        if not field:
            md = md.checkmessage(options.getfieldvalue('message', "field '{}' is empty".format(fieldname)))

    #Check size
    if options.exist('size'):
        fieldsize = options.getfieldvalue('size')
        if type(fieldsize) == str:
            if m.strcmp(fieldsize, 'universal'):

                #Check that vector size will not be confusing for ModelProcessorx
                if (md.mesh.numberofvertices == md.mesh.numberofelements):
                    raise Exception('number of vertices is the same as number of elements')
                elif (md.mesh.numberofvertices + 1 == md.mesh.numberofelements):
                    raise Exception('number of vertices + 1 is the same as number of elements')
                elif (md.mesh.numberofvertices == md.mesh.numberofelements + 1):
                    raise Exception('number of vertices is the same as number of elements + 1')

                #Uniform field
                if (field.shape[0] == 1):
                    if (np.ndim(field) > 1 and np.shape(field)[1] != 1):
                        md = md.checkmessage(options.getfieldvalue('message', "field '{}' is not supported".format(fieldname)))

                #vertex oriented input, only one column allowed
                elif (np.shape(field)[0] == md.mesh.numberofvertices):
                    if (np.ndim(field) > 1 and np.shape(field)[1] != 1):
                        md = md.checkmessage(options.getfieldvalue('message', "field '{}' is not supported".format(fieldname)))

                #element oriented input, one or more column (patch) is ok
                elif (np.shape(field)[0] == md.mesh.numberofelements):
                    #nothing to do here (either constant per element, or defined on nodes)
                    pass

                #vertex time series
                elif (np.shape(field)[0] == md.mesh.numberofvertices + 1):
                    if (np.ndim(field) > 1 and np.shape(field)[1] <= 1):
                        md = md.checkmessage(options.getfieldvalue('message', "field '{}' is not supported".format(fieldname)))

                #element time series
                elif (np.shape(field)[0] == md.mesh.numberofelements + 1):
                    if (np.ndim(field) > 1 and np.shape(field)[1] <= 1):
                        md = md.checkmessage(options.getfieldvalue('message', "field '{}' is not supported".format(fieldname)))

                #else not supported
                else:
                    md = md.checkmessage(options.getfieldvalue('message', "field '{}' is not supported".format(fieldname)))

            else:
                raise Exception("fieldsize '{}' not supported yet".format(fieldsize))

        else:
            if len(np.shape(field)) < len(fieldsize):
                md = md.checkmessage(options.getfieldvalue('message', "field {} has size {} but should be size {}".format(fieldname, np.shape(field), fieldsize)))
            else:
                for i in range(np.size(fieldsize)):
                    if (not np.isnan(fieldsize[i])) and (np.shape(field)[i] != fieldsize[i]):
                        md = md.checkmessage(options.getfieldvalue('message', "field {} dimension  # {} should be of size {}".format(fieldname, i, fieldsize[i])))

    #Check numel
    if options.exist('numel'):
        fieldnumel = options.getfieldvalue('numel')
        if (type(fieldnumel) == int and np.size(field) != fieldnumel) or (type(fieldnumel) == list and np.size(field) not in fieldnumel):
            if len(fieldnumel) == 1:
                md = md.checkmessage(options.getfieldvalue('message', "field '{}' size should be {}".format(fieldname, fieldnumel)))
            elif len(fieldnumel) == 2:
                md = md.checkmessage(options.getfieldvalue('message', "field '{}' size should be {} or {}".format(fieldname, fieldnumel[0], fieldnumel[1])))
            else:
                md = md.checkmessage(options.getfieldvalue('message', "field '{}' size should be {}".format(fieldname, fieldnumel)))

    #check NaN
    if options.getfieldvalue('NaN', 0):
        if np.any(np.isnan(field)):
            md = md.checkmessage(options.getfieldvalue('message', "NaN values found in field '{}'".format(fieldname)))

    #check Inf
    if options.getfieldvalue('Inf', 0):
        if np.any(np.isinf(field)):
            md = md.checkmessage(options.getfieldvalue('message', "Inf values found in field '{}'".format(fieldname)))

    #check cell
    if options.getfieldvalue('cell', 0):
        if not isinstance(field, (tuple, list, dict)):
            md = md.checkmessage(options.getfieldvalue('message', "field '{}' should be a tuple, list, or dict".format(fieldname)))

    #check values
    if options.exist('values'):
        fieldvalues = options.getfieldvalue('values')
        if False in m.ismember(field, fieldvalues):
            if len(fieldvalues) == 1:
                md = md.checkmessage(options.getfieldvalue('message', "field '{}' value should be '{}'".format(fieldname, fieldvalues[0])))
            elif len(fieldvalues) == 2:
                md = md.checkmessage(options.getfieldvalue('message', "field '{}' values should be '{}' '{}' or '{}'".format(fieldname, fieldvalues[0], fieldvalues[1])))
            else:
                md = md.checkmessage(options.getfieldvalue('message', "field '{}' should have values in {}".format(fieldname, fieldvalues)))

    #check greater
    if options.exist('>='):
        lowerbound = options.getfieldvalue('>=')
        if type(lowerbound) is str:
            lowerbound = attrgetter(lowerbound)(md)
        if np.size(lowerbound) > 1:  #checking elementwise
            if any(field < lowerbound):
                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))
        else:
            minval = np.nanmin(field)
            if options.getfieldvalue('timeseries', 0):
                minval = np.nanmin(field[:-1])
            elif options.getfieldvalue('singletimeseries', 0):
                if np.size(field) == 1:  #some singletimeseries are just one value
                    minval = field
                else:
                    minval = np.nanmin(field[0])

            if minval < lowerbound:
                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))

    if options.exist('>'):
        lowerbound = options.getfieldvalue('>')
        if type(lowerbound) is str:
            lowerbound = attrgetter(lowerbound)(md)
        if np.size(lowerbound) > 1:  #checking elementwise
            if any(field <= lowerbound):
                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))
        else:
            minval = np.nanmin(field)
            if options.getfieldvalue('timeseries', 0):
                minval = np.nanmin(field[:-1])
            elif options.getfieldvalue('singletimeseries', 0):
                if np.size(field) == 1:  #some singletimeseries are just one value
                    minval = field
                else:
                    minval = np.nanmin(field[0])

            if minval <= lowerbound:
                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))

    #check smaller
    if options.exist('<='):
        upperbound = options.getfieldvalue('<=')
        if type(upperbound) is str:
            upperbound = attrgetter(upperbound)(md)
        if np.size(upperbound) > 1:  #checking elementwise
            if any(field > upperbound):
                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))
        else:
            maxval = np.nanmax(field)
            if options.getfieldvalue('timeseries', 0):
                maxval = np.nanmax(field[:-1])
            elif options.getfieldvalue('singletimeseries', 0):
                if np.size(field) == 1:  #some singletimeseries are just one value
                    maxval = field
                else:
                    maxval = np.nanmax(field[0])
            elif hasattr(field, 'fov_forward_indices'):
                maxval = field.fov_forward_indices[0]
            if maxval > upperbound:
                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))

    if options.exist('<'):
        upperbound = options.getfieldvalue('<')
        if type(upperbound) is str:
            upperbound = attrgetter(upperbound)(md)
        if np.size(upperbound) > 1:  #checking elementwise
            if any(field >= upperbound):
                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))

        else:
            maxval = np.nanmax(field)
            if options.getfieldvalue('timeseries', 0):
                maxval = np.nanmax(field[:-1])
            elif options.getfieldvalue('singletimeseries', 0):
                if np.size(field) == 1:  #some singletimeseries are just one value
                    maxval = field.copy()
                else:
                    maxval = np.nanmax(field[0])

                if maxval >= upperbound:
                    md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))

    #check file
    if options.getfieldvalue('file', 0):
        if not os.path.exists(field):
            md = md.checkmessage("file provided in {}: {} does not exist".format(fieldname, field))

    #Check row of strings
    if options.exist('stringrow'):
        if not isinstance(field, list):
            md = md.checkmessage(options.getfieldvalue('message', "field {} should be a list".format(fieldname)))

    #Check forcings (size and times)
    if options.getfieldvalue('timeseries', 0):
        if field.shape[0] == md.mesh.numberofvertices or field.shape[0] == md.mesh.numberofelements:
            if np.ndim(field) > 1 and not np.size(field, 1) == 1:
                md = md.checkmessage(options.getfieldvalue('message', "field {} should have only one column as there are md.mesh.numberofvertices lines".format(fieldname)))
        elif field.shape[0] == md.mesh.numberofvertices + 1 or field.shape[0] == md.mesh.numberofelements + 1:
            if np.ndim(field) > 1 and not all(field[-1, :] == np.sort(field[-1, :])):
                md = md.checkmessage(options.getfieldvalue('message', "field {} columns should be sorted chronologically".format(fieldname)))
            if np.ndim(field) > 1 and any(field[-1, 0:-1] == field[-1, 1:]):
                md = md.checkmessage(options.getfieldvalue('message', "field {} columns must not contain duplicate timesteps".format(fieldname)))
        else:
            md = md.checkmessage(options.getfieldvalue('message', "field {} should have md.mesh.numberofvertices or md.mesh.numberofvertices + 1 lines".format(fieldname)))

    #Check single value forcings (size and times)
    if options.getfieldvalue('singletimeseries', 0):
        if field.shape[0] == 2:
            if not all(field[-1, :] == np.sort(field[-1, :])):
                md = md.checkmessage(options.getfieldvalue('message', "field {} columns should be sorted chronologically".format(fieldname)))
            if any(field[-1, 0:-1] == field[-1, 1:]):
                md = md.checkmessage(options.getfieldvalue('message', "field {} columns must not contain duplicate timesteps".format(fieldname)))
        elif field.shape[0] == 1:
            if np.ndim(field) > 1 and not np.size(field, 1) == 1:
                md = md.checkmessage(options.getfieldvalue('message', "field {} should be either a scalar or have 2 lines".format(fieldname)))
        else:
            md = md.checkmessage(options.getfieldvalue('message', "field {} should have 2 lines or be a scalar".format(fieldname)))

    return md