From fafbca97584a785fa6b9cdc788ad082fd82d2866 Mon Sep 17 00:00:00 2001
From: Brown <2mx@ornl.gov>
Date: Thu, 7 Apr 2022 14:55:04 -0400
Subject: [PATCH] connect to fitapi and update h5 writing
add better descriptions for indexing
into covariance matrix, and matrix setup
remove TODO in macs5.f90
---
sammy/samtry/tr073/ctest_n | 1 +
sammy/samtry/tr073/ctest_o | 1 +
sammy/samtry/tr073/ctest_p | 1 +
sammy/samtry/tr073/ctest_q | 1 +
sammy/samtry/tr073/ctest_r | 1 +
sammy/src/acs/macs5.f90 | 37 ++++++++++++++++++-------------
sammy/src/endf/CovarianceData.cpp | 2 +-
sammy/src/endf/CovarianceData.h | 24 ++++++++++----------
sammy/src/fff/mfff6.f90 | 2 +-
sammy/src/io/CoverxIO.cpp | 1 -
sammy/src/io/Hdf5IO.cpp | 36 +++++++++++++++++++++++++-----
11 files changed, 71 insertions(+), 36 deletions(-)
diff --git a/sammy/samtry/tr073/ctest_n b/sammy/samtry/tr073/ctest_n
index 952ff414f..ad08fbd58 100755
--- a/sammy/samtry/tr073/ctest_n
+++ b/sammy/samtry/tr073/ctest_n
@@ -1,3 +1,4 @@
t73nnn.inp
t73nnn.acs
t73nnn.ndf
+
diff --git a/sammy/samtry/tr073/ctest_o b/sammy/samtry/tr073/ctest_o
index a8a6539f9..2b4e55862 100755
--- a/sammy/samtry/tr073/ctest_o
+++ b/sammy/samtry/tr073/ctest_o
@@ -3,3 +3,4 @@ t73ooo.acs
t73ooo.grd
t73nnn.ndf
+
diff --git a/sammy/samtry/tr073/ctest_p b/sammy/samtry/tr073/ctest_p
index 67d38fcc5..7b20fc5af 100755
--- a/sammy/samtry/tr073/ctest_p
+++ b/sammy/samtry/tr073/ctest_p
@@ -1,3 +1,4 @@
t73nnn.inp
t73ppp.acs
t73nnn.ndf
+
diff --git a/sammy/samtry/tr073/ctest_q b/sammy/samtry/tr073/ctest_q
index b563181be..784f5b317 100755
--- a/sammy/samtry/tr073/ctest_q
+++ b/sammy/samtry/tr073/ctest_q
@@ -1,3 +1,4 @@
t73nnn.inp
t73qqq.acs
t73nnn.ndf
+
diff --git a/sammy/samtry/tr073/ctest_r b/sammy/samtry/tr073/ctest_r
index d24948175..3218508f2 100755
--- a/sammy/samtry/tr073/ctest_r
+++ b/sammy/samtry/tr073/ctest_r
@@ -1,3 +1,4 @@
t73nnn.inp
t73rrr.acs
t73nnn.ndf
+
diff --git a/sammy/src/acs/macs5.f90 b/sammy/src/acs/macs5.f90
index 3d8f5ea94..290ed3584 100644
--- a/sammy/src/acs/macs5.f90
+++ b/sammy/src/acs/macs5.f90
@@ -235,7 +235,7 @@ module acs5_m
real(8), parameter :: Zero = 0.0d0
integer, parameter :: L0 = 0, L1 = 1, Iu = 59
character(len=80) :: h5filename,coverxfilename
- logical :: mt51 = .true.
+ logical :: mt51 = .true., wantf33 = .false.
! --- class types ---------------
type(GridData) :: grid
type(CoverxIO) :: coverx
@@ -252,16 +252,21 @@ module acs5_m
write(6,*) "If Coverx for File 33 is wanted, type name of H5 parameter"
write(6,*) "covariance file: (else press enter)"
read(5,'(A80)') h5filename
- write(6,*) "Reading file: ", trim(h5filename)
- write(6,*) "Please enter name of ouput Coverx file"
- read(5,'(A80)') coverxfilename
- write(6,*) "Writing file: ", trim(coverxfilename)
+ if( trim(h5filename).ne."" ) then
+ wantf33 = .true.
+ write(6,*) "Reading file: ", trim(h5filename)
+ write(6,*) "Please enter name of ouput Coverx file"
+ read(5,'(A80)') coverxfilename
+ write(6,*) "Writing file: ", trim(coverxfilename)
+ end if
write(6,*) "*******************************************************"
- ! if( filename has values ) else( dont make file 33 )
- ! read covariance from h5 file into covData object
- call h5%initialize()
- call coverx%initialize()
- call h5%readRPCM(covData,trim(h5filename))
+
+ if( wantf33 ) then
+ ! read covariance from h5 file into covData object
+ call h5%initialize()
+ call coverx%initialize()
+ call h5%readRPCM(covData,trim(h5filename))
+ end if
!
CALL Newopn (59, Sam57x, 0)
@@ -279,7 +284,6 @@ module acs5_m
IF (Kktype(I).EQ.5) Ktype(2) = I
! --------------------
! Set up C++ API memory
- ! TODO: should use GridData for all grid memory in acs/fff
! --------------------
call grid%initialize()
if(I.eq.1) then
@@ -340,7 +344,7 @@ module acs5_m
! Set MAT and MT ID numbers for the coverx
! file calculations
! ------------------------------------
- call coverx%setMatMt(Mat,Mt)
+ if( wantf33 ) call coverx%setMatMt(Mat,Mt)
! *** first entry in ENDF/B-VI file -- HEAD record
WRITE (59,10600)
@@ -367,10 +371,11 @@ module acs5_m
! Option to dump Coverx file with pointwise xs covariance
!-----------------------------------------------------
!
-
- call coverx%setCoverxFile(expData,covData)
- call coverx%writeCoverxFile(trim(coverxfilename))
- call coverx%writeFile33(Za,Awr)
+ if( wantf33 ) then
+ call coverx%setCoverxFile(expData,covData)
+ call coverx%writeCoverxFile(trim(coverxfilename))
+ call coverx%writeFile33(Za,Awr)
+ end if
RETURN
!
diff --git a/sammy/src/endf/CovarianceData.cpp b/sammy/src/endf/CovarianceData.cpp
index 7073aece1..ad90435c5 100644
--- a/sammy/src/endf/CovarianceData.cpp
+++ b/sammy/src/endf/CovarianceData.cpp
@@ -11,7 +11,7 @@ namespace sammy{
uParamValues(orig.uParamValues),
uOrigParamValues(orig.uOrigParamValues),
fitStep(orig.fitStep),
- paramIndex(orig.paramIndex){
+ covarIndex(orig.covarIndex){
covariance = NULL;
if (orig.covariance != NULL) {
covariance = std::make_shared<endf::ResonanceCovariance>(*orig.covariance);
diff --git a/sammy/src/endf/CovarianceData.h b/sammy/src/endf/CovarianceData.h
index 0fed6c12a..700a98da9 100644
--- a/sammy/src/endf/CovarianceData.h
+++ b/sammy/src/endf/CovarianceData.h
@@ -176,15 +176,15 @@ namespace sammy{
* Get the number of total parameters: The pup'd and non-pup'd
*/
int getNumTotalParam() const{
- int np = (int)paramIndex.size();
+ int np = (int)covarIndex.size();
if( np != 0 ) return np;
else{
- return numParameters;
+ return numParameters; // covarIndex has not always been set
}
}
void resetTotalParam(int size){
- paramIndex.resize(size);
+ covarIndex.resize(size);
}
/**
@@ -216,10 +216,10 @@ namespace sammy{
* @param parIndex the parameter index
*/
void setCovIndex(int covIndex, int parIndex){
- if ((int)paramIndex.size() < (parIndex+1)){
- paramIndex.resize(parIndex+1);
+ if ((int)covarIndex.size() < (parIndex+1)){
+ covarIndex.resize(parIndex+1);
}
- paramIndex[parIndex] = covIndex;
+ covarIndex[parIndex] = covIndex;
}
/**
@@ -241,20 +241,20 @@ namespace sammy{
if (parIndex < 0 ){
return -2;
}
- if (parIndex >= (int)paramIndex.size()) return parIndex;
+ if (parIndex >= (int)covarIndex.size()) return parIndex;
- return paramIndex[parIndex];
+ return covarIndex[parIndex];
}
/**
* Is the indicated parameter a puped parameter
- * @param parIndex the index to check
+ * @param parIndex the index to user input varied pars
* @return true if a puped parameters, false otherwise
*/
bool isPupedParameter(int parIndex) const{
- if (parIndex < 0 || parIndex >= (int)paramIndex.size()) return false;
- if (paramIndex[parIndex] >= numParameters) return true;
+ if (parIndex < 0 || parIndex >= (int)covarIndex.size()) return false;
+ if (covarIndex[parIndex] >= numParameters) return true;
return false;
}
@@ -353,7 +353,7 @@ namespace sammy{
std::vector<double> fitStep;
/** The index of varied/pup'd parameters to resonance objects */
- std::vector<int> paramIndex;
+ std::vector<int> covarIndex;
/** Indices of parameters that are varied but don't matter, so they can be ignored */
std::vector<int> irrelevant;
diff --git a/sammy/src/fff/mfff6.f90 b/sammy/src/fff/mfff6.f90
index 62c23f9ff..522fe1212 100644
--- a/sammy/src/fff/mfff6.f90
+++ b/sammy/src/fff/mfff6.f90
@@ -46,7 +46,7 @@ module fff6_m
READ (32) (Vrpr(N),N=1,Nnn1)
call covData%getUCovariance(uCov)
- call covData%resetTotalParam(Nvpall1) ! no PUPs in URR
+ call covData%resetTotalParam(Nvpall1) ! TODO: add PUPs in URR (Nvpall1 doesn't inlcude PUPs)
call setCovMatrix(Nvpall1, uCov, Vrpr)
!
READ (32) (A,N=1,Numexc), (Dum(N),N=1,Numexc), (A,N=1,Numexc)
diff --git a/sammy/src/io/CoverxIO.cpp b/sammy/src/io/CoverxIO.cpp
index 633872c8a..763b85865 100644
--- a/sammy/src/io/CoverxIO.cpp
+++ b/sammy/src/io/CoverxIO.cpp
@@ -10,7 +10,6 @@ namespace sammy{
// variable setup
int numGrids = gdList.getLength();
int numPar = cov.getNumTotalParam();
- std::cerr << "Num par = " << numPar << std::endl;
endf::ResonanceCovariance * parCov = cov.getUCovariance();
diff --git a/sammy/src/io/Hdf5IO.cpp b/sammy/src/io/Hdf5IO.cpp
index d3da2fbbe..791f59934 100644
--- a/sammy/src/io/Hdf5IO.cpp
+++ b/sammy/src/io/Hdf5IO.cpp
@@ -115,9 +115,13 @@ namespace sammy {
// ---- set up the parameter with u-param to start ----------
std::vector<double> param;
+ std::vector<int> covind;
int numParam = cov.getNumTotalParam();
endf::ResonanceCovariance * parCov;
- for( int i=0; i<numParam; ++i ) param.push_back(cov.getUParamValue(i));
+ for( int i=0; i<numParam; ++i ){
+ param.push_back(cov.getUParamValue(i));
+ covind.push_back(cov.getCovIndex(i));
+ }
// ---- get appropriate (p||u) covariance and parameters ----
if( rpcmIsReduced ){
parCov = cov.getUCovariance();
@@ -130,9 +134,14 @@ namespace sammy {
// ---- check if length of params match covariance ----------
int numRows = parCov->getNumRows();
if( numParam != numRows || numRows != param.size() ) {
- std::string msg("covariance and parameters don't match in length!");
- throw std::invalid_argument(msg);
+ std::string msg("\ncovariance and params don't match in length!\n");
+ msg += "cov : " + std::to_string(numRows) + "\n";
+ msg += "covpar: " + std::to_string(numParam) + "\n";
+ msg += "param : " + std::to_string(param.size()) + "\n";
+ msg += "indices must be used to match params and cov\n";
+ std::cerr << msg << std::endl;
}
+ std::cout << numRows << "x" << numRows << " matrix..." << std::endl;
double data[numRows][numRows];
for( int i=0; i<numRows; ++i ){
@@ -158,15 +167,25 @@ namespace sammy {
int prank = sizeof(pdim) / sizeof(hsize_t);
std::string paramname("/param/");
+ // --- covind space ----
+ hsize_t idim[1];
+ idim[0] = numParam;
+ int irank = sizeof(idim) / sizeof(hsize_t);
+ std::string covindname("/covind/");
+
H5::DataSpace cspace(crank,cdim);
H5::DataSpace pspace(prank,pdim);
+ H5::DataSpace ispace(irank,idim);
H5::H5File file = H5::H5File(hdf5Filename,H5F_ACC_TRUNC);
- H5::DataSet covdataset = H5::DataSet(file.createDataSet(covarname,H5::PredType::NATIVE_DOUBLE,cspace));
- H5::DataSet pardataset = H5::DataSet(file.createDataSet(paramname,H5::PredType::NATIVE_DOUBLE,pspace));
+ H5::DataSet covdataset = H5::DataSet(file.createDataSet(covarname, H5::PredType::NATIVE_DOUBLE,cspace));
+ H5::DataSet pardataset = H5::DataSet(file.createDataSet(paramname, H5::PredType::NATIVE_DOUBLE,pspace));
+ H5::DataSet inddataset = H5::DataSet(file.createDataSet(covindname,H5::PredType::NATIVE_INT, ispace));
// write param covariance to HDF5 file
covdataset.write( data, H5::PredType::NATIVE_DOUBLE );
// write parameters to HDF5 file
pardataset.write( param.data(), H5::PredType::NATIVE_DOUBLE );
+ // write indices of cov-matrix for params to HDF5 file
+ inddataset.write( covind.data(), H5::PredType::NATIVE_INT );
}
@@ -179,24 +198,31 @@ namespace sammy {
H5::H5File file(hdf5Filename, H5F_ACC_RDONLY);
H5::DataSet h5covSet = file.openDataSet( "covar" );
H5::DataSet h5parSet = file.openDataSet( "param" );
+ H5::DataSet h5indSet = file.openDataSet( "covind" );
H5::DataSpace h5covSpace = h5covSet.getSpace();
H5::DataSpace h5parSpace = h5parSet.getSpace();
+ H5::DataSpace h5indSpace = h5indSet.getSpace();
// --- get dimensions of data to size local arrays
hsize_t cov_dims[2];
hsize_t par_dims[1];
+ hsize_t ind_dims[1];
int ndims = h5covSpace.getSimpleExtentDims( cov_dims, NULL);
ndims = h5parSpace.getSimpleExtentDims( par_dims, NULL);
+ ndims = h5parSpace.getSimpleExtentDims( ind_dims, NULL);
double h5cov[cov_dims[0]][cov_dims[1]];
double h5par[par_dims[0]];
+ double h5ind[ind_dims[0]];
// --- read data into local arrays
h5covSet.read(h5cov,H5::PredType::NATIVE_DOUBLE);
h5parSet.read(h5par,H5::PredType::NATIVE_DOUBLE);
+ h5indSet.read(h5ind,H5::PredType::NATIVE_INT);
// --- load into CovarianceData object
cov.makeNewUCovariance();
endf::ResonanceCovariance * parCov = cov.getUCovariance();
cov.clearUParameterValues();
for( int i=0; i<cov_dims[0]; ++i ){
cov.setUParamValue( i,h5par[i] );
+ cov.setCovIndex ( i,h5ind[i] );
for( int j=0; j<cov_dims[1]; ++j ){
parCov->setCovariance( i,j,h5cov[i][j] );
}
--
GitLab