# Source code for mpi4py_fft.io.nc_file

    import os
    import numpy as np
    from mpi4py import MPI
    from .file_base import FileBase

    # https://github.com/Unidata/netcdf4-python/blob/master/examples/mpi_example.py
    # Note. Not using groups because Visit does not understand it

    __all__ = ('NCFile',)

    comm = MPI.COMM_WORLD

    class NCFile(FileBase):
        """Class for writing data to NetCDF4 format

        Parameters
        ----------
        ncname : str
            Name of netcdf file to be created
        domain : Sequence, optional
            An optional spatial mesh or domain to go with the data.
            Sequence of either

                - 2-tuples, where each 2-tuple contains the (origin, length)
                  of each dimension, e.g., (0, 2*pi).
                - Arrays of coordinates, e.g., np.linspace(0, 2*pi, N). One
                  array per dimension.
        mode : str
            ``r``, ``w`` or ``a`` for read, write or append. Default is ``a``.
        clobber : bool, optional
            If True (default), opening a file with mode='w' will clobber an
            existing file with the same name. If False, an exception will be
            raised if a file with the same name already exists.
        kw : dict, optional
            Optional additional keyword arguments used when creating the file
            used to store data.

        Note
        ----
        Each class instance creates one unique NetCDF4-file, with one step-counter.
        It is possible to store multiple fields in each file, but all snapshots of
        the fields must be taken at the same time. If you want one field stored
        every 10th timestep and another every 20th timestep, then use two different
        class instances with two different filenames ``ncname``.
        """
        def __init__(self, ncname, domain=None, mode='a', clobber=True, **kw):
            FileBase.__init__(self, ncname, domain=domain)
            from netCDF4 import Dataset
            # netCDF4 does not seem to handle 'a' if the file does not already exist
            if mode == 'a' and not os.path.exists(ncname):
                mode = 'w'
            self.f = Dataset(ncname, mode=mode, clobber=clobber, parallel=True,
                             comm=comm, **kw)
            self.dims = None
            if 'time' not in self.f.variables:
                self.f.createDimension('time', None)
                self.f.createVariable('time', float, ('time'))
            self.close()

        def _check_domain(self, group, field):
            N = field.global_shape[field.rank:]
            if self.domain is None:
                self.domain = []
                for i in range(field.dimensions):
                    self.domain.append(np.linspace(0, 2*np.pi, N[i]))

            assert len(self.domain) == field.dimensions
            if len(self.domain[0]) == 2:
                d = self.domain
                self.domain = []
                for i in range(field.dimensions):
                    self.domain.append(np.linspace(d[i][0], d[i][1], N[i]))

            self.dims = ['time']
            for i in range(field.rank):
                ind = 'ijk'[i]
                self.dims.append(ind)
                if not ind in self.f.variables:
                    self.f.createDimension(ind, field.dimensions)
                    n = self.f.createVariable(ind, float, (ind))
                    n[:] = np.arange(field.dimensions)

            for i in range(field.dimensions):
                xyz = 'xyzrst'[i]
                self.dims.append(xyz)
                if not xyz in self.f.variables:
                    self.f.createDimension(xyz, N[i])
                    nc_xyz = self.f.createVariable(xyz, float, (xyz))
                    nc_xyz[:] = self.domain[i]

            self.f.sync()

        @staticmethod
        def backend():
            return 'netcdf4'

        def open(self, mode='r+'):
            from netCDF4 import Dataset
            self.f = Dataset(self.filename, mode=mode, parallel=True, comm=comm)

        def write(self, step, fields, **kw):
            """Write snapshot ``step`` of ``fields`` to NetCDF4 file

            Parameters
            ----------
            step : int
                Index of snapshot.
            fields : dict
                The fields to be dumped to file. (key, value) pairs are group name
                and either arrays or 2-tuples, respectively. The arrays are complete
                arrays to be stored, whereas 2-tuples are arrays with associated
                *global* slices.
            as_scalar : boolean, optional
                Whether to store rank > 0 arrays as scalars. Default is False.

            Example
            -------
            >>> from mpi4py import MPI
            >>> from mpi4py_fft import PFFT, NCFile, newDistArray
            >>> comm = MPI.COMM_WORLD
            >>> T = PFFT(comm, (15, 16, 17))
            >>> u = newDistArray(T, forward_output=False, val=1)
            >>> v = newDistArray(T, forward_output=False, val=2)
            >>> f = NCFile('ncfilename.nc', mode='w')
            >>> f.write(0, {'u': [u, (u, [slice(None), 4, slice(None)])],
            ...             'v': [v, (v, [slice(None), 5, 5])]})
            >>> f.write(1, {'u': [u, (u, [slice(None), 4, slice(None)])],
            ...             'v': [v, (v, [slice(None), 5, 5])]})

            This stores the following datasets to the file ``ncfilename.nc``.
            Using in a terminal 'ncdump -h ncfilename.nc', one gets::

                netcdf ncfilename {
                dimensions:
                        time = UNLIMITED ; // (2 currently)
                        x = 15 ;
                        y = 16 ;
                        z = 17 ;
                variables:
                        double time(time) ;
                        double x(x) ;
                        double y(y) ;
                        double z(z) ;
                        double u(time, x, y, z) ;
                        double u_slice_4_slice(time, x, z) ;
                        double v(time, x, y, z) ;
                        double v_slice_5_5(time, x) ;
                }

            """
            self.open()
            nc_t = self.f.variables.get('time')
            nc_t.set_collective(True)
            it = nc_t.size
            if step in nc_t.__array__(): # If already stored at this step previously
                it = np.argwhere(nc_t.__array__() == step)[0][0]
            else:
                nc_t[it] = step
            FileBase.write(self, it, fields, **kw)
            self.close()

        def read(self, u, name, **kw):
            step = kw.get('step', 0)
            self.open()
            s = u.local_slice()
            s = (step,) + s
            u[:] = self.f[name][s]
            self.close()

        def _write_slice_step(self, name, step, slices, field, **kw):
            assert name not in self.dims # Crashes if user tries to name fields x, y, z, .
            rank = field.rank
            slices = list((slice(None),)*rank + tuple(slices))
            slname = self._get_slice_name(slices[rank:])
            s = field.local_slice()
            slices, inside = self._get_local_slices(slices, s)
            sp = np.nonzero([isinstance(x, slice) for x in slices])[0]
            sf = np.take(s, sp)
            sdims = ['time'] + list(np.take(self.dims, np.array(sp)+1))
            fname = "_".join((name, slname))
            if fname not in self.f.variables:
                h = self.f.createVariable(fname, field.dtype, sdims)
            else:
                h = self.f.variables[fname]
            h.set_collective(True)

            h[step] = 0 # collectively create dataset
            h.set_collective(False)
            sf = tuple([step] + list(sf))
            sl = tuple(slices)
            if inside:
                h[sf] = field[sl]
            h.set_collective(True)
            self.f.sync()

        def _write_group(self, name, u, step, **kw):
            assert name not in self.dims # Crashes if user tries to name fields x, y, z, .
            s = u.local_slice()
            if name not in self.f.variables:
                h = self.f.createVariable(name, u.dtype, self.dims)
            else:
                h = self.f.variables[name]
            h.set_collective(True)
            s = (step,) + s
            h[s] = u
            self.f.sync()