Commit 07104e31 authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Implementation de lecture partielle de tableaux depuis un dataset de HDF5,...

Implementation de lecture partielle de tableaux depuis un dataset de HDF5, methode HDF5ArrayHandler<T>::ReadSubset() et operateur >> defini pour Range_Array_Grp<T> , Reza
parent 97b716c5
......@@ -21,6 +21,32 @@ namespace SOPHYA {
/*!
\ingroup HDF5IOServer
\brief HDF5 I/O handler for array objects
\warning : see the note below concerning SOPHYA array dimensions and HDF5 convention concerning the arrays.
\sa SOPHYA::TArray<T>
\verbatim
The SOPHYA TArray<T> has the following indexing and memory layout convention.
Sizes and array element indexes are ordered as (sizeX, sizeY, sizeZ, sizeT, sizeU)
and (indX, indY, indZ, ...) , where elements along X axis are contiguous in memory,
for a packed array.
Consider the following declations:
TArray<int> IA(sX, sY, sZ); // SOPHYA c++ array object declaration
int itab[sZ][sY][sX]; // standard c/c++ native array declaration
These two arrays have the same memory layout, and the element
IA(ix, jy, kz) and itab[kz][jy][ix]
would have the same offset relative to the first element.
One sees that the order of sizes or array index elements are reversed between
the list of arguments for SOPHYA arrays using the bracket operator (sX,sY,sZ)
and native c/c++ arrays with square bracket operator [sZ][sY][sX].
In HDF5, array size are labeled N1, N2, N3 ... matching the order of writing
(from left to right) of the c/c++ array declaration [N1][N2][N3].
So an HDF data set with size N1 x N2 x N3 is mapped to a SOPHYA
N1 x N2 x N3 -> TArray<T>(N3, N2, N1)
\endverbatim
*/
template <class T>
......@@ -28,6 +54,7 @@ class HDF5ArrayHandler : public HDF5HandlerInterface {
public :
HDF5ArrayHandler() { dobj=NULL; ownobj=true; }
HDF5ArrayHandler(TArray< T > & obj) { dobj = &obj; ownobj=false; }
HDF5ArrayHandler(TArray< T > * pobj) { dobj = pobj; ownobj=false; }
virtual ~HDF5ArrayHandler() { if (ownobj && dobj) delete dobj; }
virtual AnyDataObj* DataObj() { return(dobj); }
......@@ -101,14 +128,22 @@ public :
inline operator T&() { return(*dobj); }
//----- Ecriture
//! Writes the complete array as an IMAGE HDU
/*! \brief Writes the complete array as an IMAGE HDU
\warning: a packed array is created for writing if the original array is not packed
*/
virtual void Write(HDF5InOutFile& ohs, const char* name) const
{
if (( dobj == NULL) || (dobj->Size() < 1))
throw NullPtrError("HDF5ArrayHandler<T>::Write() dobj=NULL or dobj->Size()=0");
if (!(dobj->IsPacked()))
throw NotAvailableOperation("HDF5ArrayHandler<T>::Write() HDF5 I/O NOT supported for unpacked arrays");
TArray<T> * pdobj = dobj;
TArray<T> parr;
if (!(dobj->IsPacked())) {
// cout<<"SOPHYA::HDF5ArrayHandler<T>::Write() creating a packed array before writing it to file"<<endl;
// throw NotAvailableOperation("HDF5ArrayHandler<T>::Write() HDF5 I/O NOT supported for unpacked arrays");
parr=dobj->PackElements();
pdobj=&parr;
}
// -- creating the dataset
hsize_t sizes[BASEARRAY_MAXNDIMS] = {0,0,0,0,0};
/* Note that in HDF5 convention, the LAST dimension is the FASTEST changing, contiguous in memory.
......@@ -136,7 +171,7 @@ public :
}
ohs.WriteAttribute(dataset,"SOPHYA_ClassName",mtv);
// writing array data
ohs.FullDataSetWrite(dataset, HDF5Types::NativeDataType(x), dobj->Data());
ohs.FullDataSetWrite(dataset, HDF5Types::NativeDataType(x), pdobj->Data());
// writing the DVList Info() object as attributes
ohs.WriteAttributes(dataset, dobj->Info());
return;
......@@ -202,6 +237,92 @@ public :
HDF5IdWrapper::closeDataTypeForCompound(fcdtyp);
return;
}
//----- Lecture
/*! \brief Resize the array and reads a subset of the dataset to the array.
The dataset subset (or slab) and hence the array size are defined by the \b subs argument.
the \b subs element order matches the one of the HDF5 dataset.
If the last argument (fgcompact) is true, CompactAllDimensions() is called on the array after read.
\sa SOPHYA::Range_Array_Grp<T>
\sa SOPHYA::RangeArrayGroup
*/
virtual void ReadSubset(HDF5InOutFile& ihs, const char* objpath, std::vector<Range>& subs, bool fgcompact=true)
{
H5O_type_t otype=ihs.GetObjectType(objpath);
if (otype != H5O_TYPE_DATASET) throw HDF5IOException("HDF5ArrayHandler<T>::Read() HDF5 object not a DataSet");
HDF5IdWrapper datatype;
HDF5IdWrapper dataspace;
HDF5IdWrapper dataset = ihs.OpenDataSet(objpath, datatype, dataspace);
T x;
hid_t rdmemtyp=HDF5Types::NativeDataType(x);
hid_t dtm = HDF5Types::DataType(x);
bool fg1 = datatype.IsSameDataType(dtm);
bool fg2 = (fg1?true:datatype.IsCompatibleDataType(dtm));
if (!fg1 && !fg2)
throw TypeMismatchExc("HDF5ArrayHandler<T>::ReadSubset() / Incompatible datatypes between file and object");
hid_t fcdtyp=H5I_INVALID_HID;
if (datatype.IsCompound()) {
if (!fg1)
throw TypeMismatchExc("HDF5ArrayHandler<T>::ReadSubset()/Compound type in file, NOT SAME as DataType(T)");
fcdtyp=HDF5IdWrapper::getCompatibleDataTypeForCompound(rdmemtyp, datatype.get_hid());
rdmemtyp=fcdtyp;
}
if (!fg1 && fg2) {
cerr << "HDF5ArrayHandler<T>::ReadSubset()/Warning: reading data from file with conversion"<<endl;
cerr << " dataset name="<<objpath<<" object data type T="<<DataTypeInfo<T>::getTypeName()<<endl;
}
// Setting array dimensions
int rank = dataspace.getDSRank();
if (rank > BASEARRAY_MAXNDIMS)
throw HDF5IOException("HDF5ArrayHandler<T>::ReadSubset() dataset rank > max SOPHYA BASEARRAY_MAXNDIMS=5");
if (rank != (int)subs.size())
throw HDF5IOException("HDF5ArrayHandler<T>::ReadSubset() rank != subs.size() (incompatible specification of the dataset subset)");
hsize_t elcnt[BASEARRAY_MAXNDIMS]={0,0,0,0,0};
hsize_t sz[BASEARRAY_MAXNDIMS]={0,0,0,0,0};
dataspace.getDSSizes(sz);
for(int_4 ir=0; ir<rank; ir++) {
if (sz[ir]<1) throw HDF5IOException("HDF5ArrayHandler<T>::ReadSubset() dataset full size <= 0");
}
hsize_t offset[BASEARRAY_MAXNDIMS]={0,0,0,0,0};
sa_size_t asz[BASEARRAY_MAXNDIMS]={0,0,0,0,0};
/* Note that in HDF5 convention, the LAST dimension is the FASTEST changing, contiguous in memory.
int A[N1][N2][N3] --> HDF5 data set with size[0]=N1, size[1]=N2, size[2]=N3
correspond to a Sophya array TArray<int>(N3,N2,N1) */
size_t fullsz=1;
int_4 irs=rank;
for(int_4 ir=0; ir<rank; ir++) {
irs--;
if (subs[ir]==Range::all()) { asz[irs] = (sa_size_t)sz[ir]; elcnt[ir]=sz[ir]; offset[ir]=0; }
else if (subs[ir]==Range::first()) { asz[irs] = (sa_size_t)1; elcnt[ir]=1; offset[ir]=0; }
else if (subs[ir]==Range::last()) { asz[irs] = (sa_size_t)1; elcnt[ir]=1; offset[ir]=sz[ir]-1; }
else {
asz[irs] = subs[ir].Size(); elcnt[ir]=(hsize_t)subs[ir].Size(); offset[ir]=(hsize_t)subs[ir].Start();
if (offset[ir]+elcnt[ir]>sz[ir])
throw HDF5IOException("HDF5ArrayHandler<T>::ReadSubset() - Range ERROR (offset+elcnt)>dataset_size");
}
}
if ( dobj == NULL ) {
dobj = new TArray<T> ;
ownobj = true;
}
//DEL cout<<"*DBG*HDF5ArrayHandler<T>::ReadSubset() - setting array size, rank="<<rank<<" sizes=";
//DEL for(int ij=0; ij<rank; ij++) cout<<asz[ij]<<" x "; cout<<endl;
dobj->SetSize(rank, asz);
// Reading array data
ihs.DataSetRead(dataset, offset, elcnt, rdmemtyp, dobj->Data());
// Reading attributes
ihs.ReadAttributes(dataset, dobj->Info());
dobj->Info().DeleteKey("SOPHYA_HDF5Handler_ClassName");
dobj->Info().DeleteKey("SOPHYA_ClassName");
HDF5IdWrapper::closeDataTypeForCompound(fcdtyp);
if (fgcompact) dobj->CompactAllDimensions();
return;
}
protected :
TArray<T> * dobj;
......@@ -224,8 +345,9 @@ inline HDF5InOutFile& operator << (HDF5InOutFile& os, TArray<T> const & arr)
return os;
}
/*!
\brief operator >> overload to read a TArray<T> from a HDF5InOutFile file wrapper object
\brief operator >> overload to read a TArray<T> from a HDF5InOutFile using a file wrapper object
Use the HDF5NameTag object (pull it from the HDF5InOutFile stream) to define the dataset name to be read
*/
......@@ -239,6 +361,21 @@ inline HDF5InOutFile& operator >> (HDF5InOutFile& is, TArray<T> & arr)
}
/*!
\brief operator >> overload to read a subset of a dataset of an HDF5InOutFile int an array
Use the HDF5NameTag object (pull it from the HDF5InOutFile stream) to define the dataset name to be read
*/
template <class T>
inline HDF5InOutFile& operator >> (HDF5InOutFile& is, Range_Array_Grp<T> sub_arr)
{
if (!sub_arr.arrp_) throw NullPtrError("operator >>(HDF5InOutFile, Range_Array_Grp<T>> ) null array pointer");
HDF5ArrayHandler<T> fio(sub_arr.arrp_);
string san=is.GetNextObjectOpenName();
fio.ReadSubset(is, san.c_str(), sub_arr.vr_);
return(is);
}
} // Fin du namespace SOPHYA
#endif /* HDF5ARRHAND_H_SEEN */
......
......@@ -15,6 +15,8 @@
int th5_debug(string fname, int sx, int sy, int prtlev=0, HDF5Types::ByteOrder bo=HDF5Types::BO_Auto);
//---- Testing I/O for arrays
int th5_arrays(string fname, int sx, int sy, int prtlev=0, HDF5Types::ByteOrder bo=HDF5Types::BO_Auto);
//---- Testing partial I/O for arrays
int th5_subarrays(string fname, int prtlev=0, HDF5Types::ByteOrder bo=HDF5Types::BO_Auto);
//---- Testing I/O for NTuples
int th5_ntuple(string fname, int IMAX, int JMAX, bool fgdouble, int blk, bool fgautoname,
int prtlev=0, HDF5Types::ByteOrder bo=HDF5Types::BO_Auto);
......@@ -60,20 +62,24 @@ int main(int narg, char* arg[])
if (fname=="debug") {
rc = th5_debug(fname, sx, sy, prt, h5bo);
rc += th5_subarrays(fname, prt, h5bo);
}
else {
rc = th5_arrays(fname, sx, sy, prt, h5bo);
PrtTim("Check arrays done ");
int rcsa = th5_subarrays(fname, prt, h5bo);
PrtTim("Check sub-arrays done ");
if (rcsa) rc+=8;
int rcnt = th5_ntuple(fname, sx, sy, false, blk, false, prt, h5bo);
if (rcnt) rc+=8;
rcnt = th5_ntuple(fname, sx, sy, true, blk, true, prt, h5bo);
if (rcnt) rc+=16;
rcnt = th5_ntuple(fname, sx, sy, true, blk, true, prt, h5bo);
if (rcnt) rc+=32;
PrtTim("Check NTuple done ");
rcnt = th5_segdb(fname, sx, sy, prt);
if (rcnt) rc+=32;
if (rcnt) rc+=64;
PrtTim("Check SegDataBlock<T> done ");
rcnt = th5_dtable(fname, sx, sy, prt, h5bo);
if (rcnt) rc+=64;
if (rcnt) rc+=128;
PrtTim("Check DataTable done ");
}
......@@ -291,7 +297,128 @@ int th5_arrays(string fname, int sx, int sy, int prtlev, HDF5Types::ByteOrder b
return rc;
}
inline int ck_diff_ia(TArray<int_4>& a, TArray<int_4>& b, int_4& imin, int_4& imax)
{
TArray<int_4> dIA=a-b;
dIA.MinMax(imin,imax);
if ((imin!=0)||(imax!=0)) return 1;
return 0;
}
/* -- Fonction -- */
int th5_subarrays(string fname, int prtlev, HDF5Types::ByteOrder bo)
{
int rc=0;
string sbo="BO_Auto";
if (bo==HDF5Types::BO_LittleEndian) sbo="BO_LittleEndian";
if (bo==HDF5Types::BO_BigEndian) sbo="BO_BigEndian";
cout << "======= th5_subarrays(): testing array partial I/O with HDF5IOServer (int_4)"<<endl;
sa_size_t sz=12, sy=25, sx=30;
// We create a integer array SizeX x SizeY
TArray<int_4> IA(sx,sy,sy);
for(sa_size_t k=0; k<sz; k++) {
for(sa_size_t j=0; j<sy; j++)
for(sa_size_t i=0; i<sx; i++)
IA(i,j,k) = 1000*(k+1)+200*(j+1)+3*(i+1);
}
TArray<int_4> sIAz3 = IA(Range::all(), Range::all(), Range(3));
TArray<int_4> sIAz7 = IA(Range::all(), Range::all(), Range(7));
TArray<int_4> sIAxya = IA(Range(12,25), Range(6,20), Range::last());
TArray<int_4> sIAxyb = IA(Range(12,25), Range(6,20), Range(5,7));
fname = fname + "_suba.h5";
{
if (prtlev>0) cout << "1/ Writing five arrays (IA,sIAz3,sIAz7,sIAxya,sIAxyb) to HDF5 file: "<<fname<<endl;
HDF5InOutFile h5os(fname, HDF5InOutFile::CreateOverwrite);
h5os.setDefaultByteOrdering(bo);
h5os << HDF5NameTag("IA") << IA;
h5os << HDF5NameTag("sIAz3") << sIAz3;
h5os << HDF5NameTag("sIAz7") << sIAz7;
h5os << HDF5NameTag("sIAxya") << sIAxya;
h5os << HDF5NameTag("sIAxyb") << sIAxyb;
}
{
if (prtlev>0) cout << "2/ Reading five arrays from HDF5 file: "<<fname<<endl;
TArray<int_4> IAH, sIAz3H, sIAz7H, sIAxyaH, sIAxybH;
HDF5InOutFile h5is(fname, HDF5InOutFile::RO);
h5is >> HDF5NameTag("IA") >> IAH;
h5is >> HDF5NameTag("sIAz3") >> sIAz3H;
h5is >> HDF5NameTag("sIAz7") >> sIAz7H;
h5is >> HDF5NameTag("sIAxya") >> sIAxyaH;
h5is >> HDF5NameTag("sIAxyb") >> sIAxybH;
if (prtlev > 1) {
cout << " ----- Array IAH From HDF5 = \n " << IAH << endl;
cout << " ----- Array sIAz3H From HDF5 = \n " << sIAz3H << endl;
cout << " ----- Array sIAz7H From HDF5 = \n " << sIAz7H << endl;
cout << " ----- Array sIAxyaH From HDF5 = \n " << sIAxyaH << endl;
cout << " ----- Array sIAxybH From HDF5 = \n " << sIAxybH << endl;
}
int_4 imin,imax;
int rcck=0;
if (ck_diff_ia(IA,IAH,imin,imax)) {
if (prtlev>0) cout << "3.a/ Error in reading IA , Diff: min="<<imin<<" max="<<imax<<endl;
rcck += 1;
}
if (ck_diff_ia(sIAz3,sIAz3H,imin,imax)) {
if (prtlev>0) cout << "3.b/ Error in reading sIAz3 , Diff: min="<<imin<<" max="<<imax<<endl;
rcck += 2;
}
if (ck_diff_ia(sIAz7,sIAz7H,imin,imax)) {
if (prtlev>0) cout << "3.c/ Error in reading sIAz7 , Diff: min="<<imin<<" max="<<imax<<endl;
rcck += 2;
}
if (ck_diff_ia(sIAxya,sIAxyaH,imin,imax)) {
if (prtlev>0) cout << "3.d/ Error in reading sIAxya , Diff: min="<<imin<<" max="<<imax<<endl;
rcck += 2;
}
if (ck_diff_ia(sIAxyb,sIAxybH,imin,imax)) {
if (prtlev>0) cout << "3.e/ Error in reading sIAxyb , Diff: min="<<imin<<" max="<<imax<<endl;
rcck += 2;
}
if (rcck==0) cout << "3/ Check 5 array read done -> OK"<<endl;
else cout << "3/ Check 5 array read done -> ERROR rcck="<<rcck<<endl;
if (prtlev>0) cout << "4/ Reading sub-arrays from HDF5 file: "<<fname<<endl;
TArray<int_4> srIAz3H, srIAz7H, srIAxyaH, srIAxybH;
h5is >> HDF5NameTag("IA") >> RangeArrayGroup(Range::all(),Range::all(),Range(3),srIAz3H);
h5is >> HDF5NameTag("IA") >> RangeArrayGroup(Range::all(),Range::all(),Range(7),srIAz7H);
h5is >> HDF5NameTag("IA") >> RangeArrayGroup(Range(12,25),Range(6,20),Range::last(),srIAxyaH);
h5is >> HDF5NameTag("IA") >> RangeArrayGroup(Range(12,25),Range(6,20),Range(5,7),srIAxybH);
if (prtlev > 1) {
cout << " ----- Array srIAz3H From HDF5 (partial read) = \n " << srIAz3H << endl;
cout << " ----- Array srIAz7H From HDF5 (partial read) = \n " << srIAz7H << endl;
cout << " ----- Array srIAxyaH From HDF5 (partial read) = \n " << srIAxyaH << endl;
cout << " ----- Array srIAxybH From HDF5 (partial read) = \n " << srIAxybH << endl;
}
rcck=0;
if (ck_diff_ia(sIAz3,srIAz3H,imin,imax)) {
if (prtlev>0) cout << "4.a/ Error in partial reading sIAz3 , Diff: min="<<imin<<" max="<<imax<<endl;
rcck += 16;
}
if (ck_diff_ia(sIAz7,srIAz7H,imin,imax)) {
if (prtlev>0) cout << "4.b/ Error in partial reading sIAz7 , Diff: min="<<imin<<" max="<<imax<<endl;
rcck += 16;
}
if (ck_diff_ia(sIAxya,srIAxyaH,imin,imax)) {
if (prtlev>0) cout << "4.c/ Error in partial reading sIAxya , Diff: min="<<imin<<" max="<<imax<<endl;
rcck += 16;
}
if (ck_diff_ia(sIAxyb,srIAxybH,imin,imax)) {
if (prtlev>0) cout << "4.d/ Error in partial reading sIAxyb , Diff: min="<<imin<<" max="<<imax<<endl;
rcck += 16;
}
if (rcck==0) cout << "5/ Check array partial read done (4 arrays) -> OK"<<endl;
else cout << "5/ Check array partial read done (4 arrays) -> ERROR rcck="<<rcck<<endl;
}
return rc;
}
/* -- Fonction -- */
int th5_ntuple(string fname, int IMAX, int JMAX, bool fgdouble, int blk, bool fgautoname, int prtlev, HDF5Types::ByteOrder bo)
{
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment