@*
HDF5Format.tmpl

Created by Graham Dennis on 2009-01-24.

Copyright (c) 2009-2012, Graham Dennis

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.

*@
@extends xpdeint.Features._HDF5Format
@from xpdeint.Geometry.SplitUniformDimensionRepresentation import SplitUniformDimensionRepresentation
@from xpdeint.CallOnceGuards import callOnceGuard, callOncePerInstanceGuard

@def description: HDF5 output format

@attr $name = 'hdf5'
@attr $mpiSafe = True

@@callOnceGuard
@def includes
  @#
  @super
  @#
#if defined(HAVE_HDF5_HL)
  #include <hdf5_hl.h>
#endif
@end def

@def writeOutFunctionImplementationBody($dict)
  @#
${writeOutFunctionImplementationBegin(dict)}@slurp

  @set $featureOrdering = ['Driver']
  @set $featureDict = dict.copy()
  @set $featureDict['extraIndent'] = 0
${insertCodeForFeatures('binaryWriteOutBegin', featureOrdering, featureDict)}@slurp
  @set $extraIndent = featureDict['extraIndent']

${writeOutFunctionContents(dict), extraIndent=extraIndent}@slurp

${insertCodeForFeaturesInReverseOrder('binaryWriteOutEnd', featureOrdering, featureDict)}@slurp

${writeOutFunctionImplementationEnd(dict)}@slurp
  @#
@end def

@def truncateOutputFiles($baseFilename)
char _dataFilename[200];
snprintf(_dataFilename, 200, "%s.h5", ${baseFilename});

H5Fclose(H5Fcreate(_dataFilename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT));
@end def


@def writeOutFunctionContents($dict)
  @#
  @set $fp = dict['fp']
  @set $baseFilename = dict['baseFilename']
  @set $groupID = dict['groupID']
  @set $field = dict['field']
  @set $basis = dict['basis']
  @set $dependentVariables = dict['dependentVariables']
  @set $componentCount = 0
  @for $variable in $dependentVariables
    @set $componentCount += len($variable.vector.components)
    @if $variable.vector.type == 'complex'
      @set $componentCount += len($variable.vector.components)
    @end if
  @end for
  @set dict['componentCount'] = componentCount
  @#
char _h5Filename[200];
snprintf(_h5Filename, 200, "%s.h5", ${baseFilename});

/* Open the file */
hid_t hdf5_file = H5Fopen(_h5Filename, H5F_ACC_RDWR, H5P_DEFAULT);
if (hdf5_file < 0) {
  _LOG(_WARNING_LOG_LEVEL, "Failed to open HDF5 file '%s', will try to create it.", _h5Filename);
  hdf5_file = H5Fcreate(_h5Filename, H5F_ACC_EXCL, H5P_DEFAULT, H5P_DEFAULT);
  if (hdf5_file < 0) {
    _LOG(_ERROR_LOG_LEVEL, "Failed to create HDF5 file '%s'. Bailing.", _h5Filename);
  }
}

/* Create the group for this data */
hid_t group;
if (!H5Lexists(hdf5_file, "/${groupID}", H5P_DEFAULT))
  group = H5Gcreate(hdf5_file, "/${groupID}", H5P_DEFAULT);
else
  group = H5Gopen(hdf5_file, "/${groupID}");

if ($fp) {
  fprintf($fp, "    <Stream><Metalink Format=\"HDF5\" Type=\"Remote\" Group=\"/${groupID}\"/>\n");
  fprintf($fp, "%s.h5\n", ${baseFilename});
  fprintf($fp, "    </Stream>\n");
}

/* Create the coordinate data sets */
hsize_t coordinate_length;
hid_t coordinate_dataspace;
  @for dim in field.dimensions
    @set $dimRep = dim.inBasis(basis)
coordinate_length = ${dimRep.globalLattice};
coordinate_dataspace = H5Screate_simple(1, &coordinate_length, NULL);
    @set $dataType = {'real': 'H5T_NATIVE_REAL', 'long': 'H5T_NATIVE_LONG'}[$dimRep.type]
hid_t dataset_${dimRep.name};
if (!H5Lexists(hdf5_file, "/${groupID}/${dimRep.name}", H5P_DEFAULT))
  dataset_${dimRep.name} = H5Dcreate(hdf5_file, "/${groupID}/${dimRep.name}", ${dataType}, coordinate_dataspace, H5P_DEFAULT);
else
  dataset_${dimRep.name} = H5Dopen(hdf5_file, "/${groupID}/${dimRep.name}");
    @if isinstance(dimRep, SplitUniformDimensionRepresentation)
      @set $dimArrayName = c'${dimRep.name}_data'
${dimRep.type}* ${dimArrayName} = (${dimRep.type}*)xmds_malloc(${dimRep.globalLattice} * sizeof(${dimRep.type}));
for (long _i0 = 0; _i0 < ${dimRep.globalLattice}; _i0++) {
  ${dimArrayName}[_i0] = ${dimRep.arrayName}[(_i0 + (${dimRep.globalLattice}+1)/2) % ${dimRep.globalLattice}];
}
    @else
      @set $dimArrayName = dimRep.arrayName
    @end if
H5Dwrite(dataset_${dimRep.name}, $dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, ${dimArrayName});
#if defined(HAVE_HDF5_HL)
  H5DSset_scale(dataset_${dimRep.name}, "${dimRep.name}");
#endif

    @if isinstance(dimRep, SplitUniformDimensionRepresentation)
xmds_free(${dimArrayName});
    @end if
H5Sclose(coordinate_dataspace);
  @end for

hsize_t file_dims[] = {${', '.join(dim.inBasis(basis).globalLattice for dim in field.dimensions)}};
hid_t file_dataspace = H5Screate_simple(${len(field.dimensions)}, file_dims, NULL);

  @for variable in dependentVariables
    @if $variable.vector.type == 'real'
      @set $variable['separatedComponents'] = list(enumerate($variable.components))
    @else
      @set $components = []
      @set variable['separatedComponents'] = components
      @for offset, componentName in enumerate($variable.components)
        @silent components.extend([(2*offset, componentName + 'R'), (2*offset+1, componentName + 'I')])
      @end for
    @end if
    @for offset, componentName in $variable.separatedComponents
hid_t dataset_${componentName};
if (!H5Lexists(hdf5_file, "/${groupID}/${componentName}", H5P_DEFAULT))
  dataset_${componentName} = H5Dcreate(hdf5_file, "/${groupID}/${componentName}", H5T_NATIVE_REAL, file_dataspace, H5P_DEFAULT);
else
  dataset_${componentName} = H5Dopen(hdf5_file, "/${groupID}/${componentName}");
#if defined(HAVE_HDF5_HL)
      @for dimNum, dim in enumerate(field.dimensions)
  H5DSattach_scale(dataset_${componentName}, dataset_${dim.inBasis(basis).name}, ${dimNum});
      @end for
#endif
    @end for
  @end for
  @for dim in field.dimensions
    @set $dimRep = dim.inBasis(basis)
H5Dclose(dataset_${dimRep.name});
  @end for

  @# This is where all of the magic MPI code goes
  @set $featureOrdering = ['Driver']
  @set $featureDict = dict.copy()
  @set $featureDict['extraIndent'] = 0
${insertCodeForFeatures('binaryWriteOutWriteDataBegin', $featureOrdering, featureDict)}@slurp
  @set $extraIndent = featureDict['extraIndent']
  @silent dict['operation'] = 'write'
  @silent dict['variables'] = dict['dependentVariables']

if (${field.sizeInBasis(basis)}) {
  ${processData(dict), autoIndent=True, extraIndent=extraIndent}@slurp
}

  @# This is where the rest of the magic MPI code goes
${insertCodeForFeaturesInReverseOrder('binaryWriteOutWriteDataEnd', $featureOrdering, featureDict)}@slurp

  @for variable in dependentVariables
    @for offset, componentName in $variable.separatedComponents
H5Dclose(dataset_${componentName});
    @end for
  @end for

H5Sclose(file_dataspace);
H5Gclose(group);
H5Fclose(hdf5_file);
  @#
@end def