#include "hdf5.h"
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define FILENAME  "bychunk.h5"
#define TRUE    1
#define FALSE   0

int main()
{
    hid_t propsID, spaceID, typeID, fileID, dsetID;
    hsize_t dsetDims[2] = { 4066763, 973 }; /* Dataset dimension size */
    hsize_t chunk_dims[2] = { 1000, 8 };    /* Chunk dimension size */
    size_t chunksize = (chunk_dims[0] * chunk_dims[1]); /* Size of one chunk */
    int numx;           /* # of chunks in the x dimension */
    int numy;           /* # of chunks in the y dimension */
    double *data;       /* Buffer for writing data */
    herr_t status;      /* Return status */
    int i, j, k;        /* Local index variables */
    hid_t mem_sid;      /* Memory dataspace ID */
    hsize_t start[2];   /* Start of hyperslab */
    hsize_t stride[2];  /* Stride of hyperslab */
    hsize_t count[2];   /* Block count */
    hsize_t block[2];   /* Block sizes */

    time_t tt,*tsec;    /* Time measurement */
    struct tm *tm;
    char msg[50];


    /* Allocate buffer for one chunk */
    data = (double *)malloc(chunksize *sizeof(double));

    /* Initialize the buffer */
    for (i = 0; i< chunksize; i++)
       data[i] = (double)NAN;

    /* Dataset creation property: set chunk layout and filter */
    propsID = H5Pcreate(H5P_DATASET_CREATE);
    status = H5Pset_deflate(propsID, 3);
    status = H5Pset_chunk(propsID, 2, chunk_dims);

    /* Create memory space for one chunk */
    mem_sid = H5Screate_simple(2, chunk_dims, NULL);

    /* Calculate the # of chunks in the x and y dimension */
    /* Round up to the next integer # of chunks, to accomodate partial chunks */
    numx = ((dsetDims[0] + chunk_dims[0]) - 1) /chunk_dims[0];
    numy = ((dsetDims[1] + chunk_dims[1]) - 1) /chunk_dims[1];

    { /* Read the dataset by chunks */

        /* Create dataspace */
        spaceID = H5Screate_simple(2, dsetDims, dsetDims);
    
        /* Create datatype */
        typeID = H5Tcopy(H5T_NATIVE_DOUBLE);

        /* Create the test file */
        fileID = H5Fcreate(FILENAME, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);

        /* Create the dataset */
        dsetID = H5Dcreate(fileID, "dataset", typeID, spaceID, H5P_DEFAULT, propsID, H5P_DEFAULT);

        /* Settings for hyperslab selection */
        start[0] = start[1] = 0;
        stride[0] = 1; stride[1] = 1;
        count[0] = 1; count[1] = 1;

        /* Measure time */
        tt=time((long*)0);
        tsec=&tt;
        tm=gmtime(tsec);
        sprintf(msg,"%02i:%02i:%02i",tm->tm_hour,tm->tm_min,tm->tm_sec);
        printf("begin time %s\n",msg);

        /* # of chunks in the x dimension */
        for(i = 0; i < numx; i++) {
            hbool_t new = FALSE;

            /* Starting position in the y dimension */
            start[1] = 0;

            /* Partial edge chunk or full chunk */
            if(((i+1)*chunk_dims[0]) > dsetDims[0]) {
                block[0] = dsetDims[0] - (i*chunk_dims[0]);
                new = TRUE;
            } else
                block[0] = chunk_dims[0];

            /* # of chunks in the y dimension */
            for(j = 0; j < numy; j++) {
                hid_t new_mem_sid = -1;

                /* Partial edge chunk or full chunk */
                if(((j+1)*chunk_dims[1]) > dsetDims[1]) {
                    block[1] = dsetDims[1] - (j*chunk_dims[1]);
                    new = TRUE;
                } else
                    block[1] = chunk_dims[1];

                status = H5Sselect_hyperslab(spaceID, H5S_SELECT_SET, start, stride, count, block);

                /* Write to the partial edge chunk or the full chunk */
                if(new) {
                    double *newdata;    /* Buffer for the partial chunk */
                    int m;              /* Local index variable */

                    /* Create memory dataspace for the partial edge chunk */
                    new_mem_sid = H5Screate_simple(2, block, NULL);

                    /* Allocate and initialize buffer for the partial edge chunk */
                    newdata = (double *)malloc(block[0] * block[1] *sizeof(double));
                    for (m = 0; m < (block[0]*block[1]); m++)
                        newdata[m] = (double)NAN;

                    /* Write to the partial edge chunk */
                    status = H5Dwrite(dsetID, typeID, new_mem_sid, spaceID, H5P_DEFAULT, newdata);

                    /* Close the memory dataspace */
                    H5Sclose(new_mem_sid);

                    /* Free the buffer */
                    free(newdata);
                } else
                    /* Write to the full chunk  */
                    status = H5Dwrite(dsetID, typeID, mem_sid, spaceID, H5P_DEFAULT, data);

                /* Move to the next chunk in the y dimension */
                start[1] += chunk_dims[1];
            }

            /* Move to the next chunk in the x dimension */
            start[0] += chunk_dims[0];
        }


        /* print out end time in hh:mm:ss */
        tt=time((long*)0);
        tsec=&tt;
        tm=gmtime(tsec);
        sprintf(msg,"%02i:%02i:%02i",tm->tm_hour,tm->tm_min,tm->tm_sec);
        printf("loop %d end %s\n\n",k,msg);

        status = H5Sclose(spaceID);
        status = H5Dclose(dsetID);
        status = H5Fclose(fileID);

        //remove ("myfile.h5");
    }

    H5Sclose(mem_sid);
    free(data);
 }