Commit cd90a62c authored by myron's avatar myron
Browse files

first test with HighFive

parent e34114e4
......@@ -2753,6 +2753,9 @@ def superR_getVolume_Esynt(mydata):
DD.tofile(DDname)
SS.tofile(SSname)
SS.tofile(SSname)
h5py.File(DSname+".h5","w")["data"] = DS.astype("f")
def superR_getVolume(mydata):
"""
......
......@@ -10,6 +10,11 @@ project (FRSV LANGUAGES CXX C)
include(CheckLanguage)
find_package(HighFive REQUIRED)
find_package(MPI )
if(MPI_C_FOUND)
......@@ -89,6 +94,6 @@ set(CMAKE_CXX_FLAGS_PIPPO "-g -Wextra -Wall ${CMAKE_CXX_FLAGS}" )
add_executable (frsv frsv.cc)
target_link_libraries( frsv PRIVATE yaml-cpp Boost::filesystem Boost::system )
target_link_libraries( frsv PRIVATE yaml-cpp Boost::filesystem Boost::system HighFive )
target_include_directories ( frsv PRIVATE "${PROJECT_SOURCE_DIR}")
......@@ -5,8 +5,16 @@
#include<stdio.h>
#include <boost/filesystem.hpp>
#include<math.h>
#include <highfive/H5File.hpp>
#include <highfive/H5DataSet.hpp>
#include <highfive/H5DataSpace.hpp>
#include <highfive/H5DataType.hpp>
#define assertm(exp, msg) assert(((void)msg, exp))
namespace h5 = HighFive;
struct DimPars_struct {
int NE ;
int NV ;
......@@ -30,7 +38,7 @@ typedef InputFileNames_struct InputFileNames;
float * read_volume(std::string fn, size_t nbytes) {
boost::filesystem::path file(fn);
if(!( boost::filesystem::file_size(file) == nbytes)){
std::cout<< std::string("File ")+fn+std::string(" has not the required size ")<< nbytes<< " it is " <<boost::filesystem::file_size(file)<<std::endl;
std::cerr<< std::string("File ")+fn+std::string(" has not the required size ")<< nbytes<< " it is " <<boost::filesystem::file_size(file)<<std::endl;
}
assert(( boost::filesystem::file_size(file) == nbytes));
float *res = (float*) malloc(nbytes);
......@@ -44,8 +52,8 @@ float * read_volume(std::string fn, size_t nbytes) {
void save_volume(std::string fn, float * X , size_t nbytes) {
FILE * stream = fopen( fn.c_str() ,"w");
size_t bread = fwrite( X, nbytes, 1, stream);
assert(bread==1);
size_t bcount = fwrite( X, nbytes, 1, stream);
assert(bcount==1);
fclose(stream);
}
......@@ -80,16 +88,45 @@ public:
this->if_names = if_names_a;
this->dimpars = dimpars_a;
{
std::string testh5_name(if_names.DSname+".h5" ) ;
h5::File test_file( testh5_name , h5::File::ReadOnly);
h5::DataSet test_dataset = test_file.getDataSet("data");
std::vector<size_t> test_dims = test_dataset.getDimensions();
printf( " %d %d %d %d %d\n", test_dims[0], test_dims[1], test_dims[2], test_dims[3], test_dims[4] );
this->dimpars.NV = test_dims[0] ;
this->dimpars.NROI = test_dims[0] ;
this->dimpars.DIMZ = test_dims[0] ;
this->dimpars.DIMY = test_dims[0] ;
this->dimpars.DIMX = test_dims[0] ;
this->DS = (float*) malloc( size_t(dimpars.NV)*
size_t(dimpars.NROI)*
size_t(dimpars.DIMZ)*
size_t(dimpars.DIMY)*
size_t(dimpars.DIMX) * sizeof(float));
test_dataset.read( this->DS , h5::AtomicType<float>() ) ;
}
this->DS = read_volume( if_names.DSname,
size_t(dimpars.NV)*size_t(dimpars.NROI)*size_t(dimpars.DIMZ)* size_t(dimpars.DIMY)* size_t(dimpars.DIMX) * sizeof(float));
size_t(dimpars.NV)*
size_t(dimpars.NROI)*
size_t(dimpars.DIMZ)*
size_t(dimpars.DIMY)*
size_t(dimpars.DIMX) * sizeof(float));
if(dimpars.ZSTART>=0) {
int dimz = 1+dimpars.ZEND - dimpars.ZSTART ;
float *ds = (float*) malloc(size_t(dimpars.NV)*size_t(dimpars.NROI)*size_t(dimz)* size_t(dimpars.DIMY)* size_t(dimpars.DIMX) * sizeof(float) ) ;
size_t BLOCCO = size_t(dimpars.DIMZ)* size_t(dimpars.DIMY)* size_t(dimpars.DIMX ) ;
size_t blocco = size_t( dimz)* size_t(dimpars.DIMY)* size_t(dimpars.DIMX ) ;
size_t BLOCK = size_t(dimpars.DIMZ)* size_t(dimpars.DIMY)* size_t(dimpars.DIMX ) ;
size_t block = size_t( dimz)* size_t(dimpars.DIMY)* size_t(dimpars.DIMX ) ;
size_t offset = size_t(dimpars.ZSTART)* size_t(dimpars.DIMY)* size_t(dimpars.DIMX) ;
......@@ -97,9 +134,9 @@ public:
for(int ir=0; ir< dimpars.NROI; ir++) {
memcpy(
ds + (iv*dimpars.NROI+ir)*blocco ,
DS + (iv*dimpars.NROI+ir)*BLOCCO+offset,
blocco*sizeof(float)
ds + (iv*dimpars.NROI+ir)*block ,
DS + (iv*dimpars.NROI+ir)*BLOCK+offset,
block*sizeof(float)
);
}
}
......@@ -123,14 +160,28 @@ public:
size_t(dimpars.NROI)*
sizeof(float));
/*{
int cnt=0;
for(int iE=0; iE< dimpars.NE; iE++) {
for(int iV=0; iV< dimpars.NV; iV++) {
for(int iR=0; iR< dimpars.NROI; iR++) {
this->coefficients[cnt] = ( iE==iV );
if(iR!=1)
this->coefficients[cnt] = 0;
cnt++;
}
}
}
}
*/
DimPars *d = & dimpars;
this->Matff = (float*) malloc( d->NE * d->NE * d->NROI * sizeof(float) ) ;
this->MatSffF2 = (float*) malloc( d->NE*d->NE* d->DIMX*d->DIMX * sizeof(float) ) ;
this->settaMatff();
this->setMatff();
this->VectA = (float*) malloc( d->NE * d->DIMZ * d->DIMY * d->DIMX * sizeof(float) ) ;
......@@ -144,12 +195,12 @@ public:
#define SS_addr( iroi, ix1, ix2 ) (( (iroi)*DIMX+ (ix1))*DIMX+ (ix2))
#define PRO_addr( iV, iroi, iz, iy, ix ) ((( ( (iV)*NROI +iroi )*DIMZ + (iz))*DIMY+ (iy))*DIMX+ (ix))
void settaVectA_and_Mat(SolVol target, FreeFactsVol F ) {
settaMatSffF2(F) ;
settaVectA( target, F ) ;
void setVectA_and_Mat(SolVol target, FreeFactsVol F ) {
setMatSffF2(F) ;
setVectA( target, F ) ;
};
void settaVectA(SolVol target, FreeFactsVol F ) {
void setVectA(SolVol target, FreeFactsVol F ) {
int NE = dimpars.NE;
int NV = dimpars.NV;
int NROI = dimpars.NROI;
......@@ -172,8 +223,23 @@ public:
}
}
}
void settaMatff() {
void project_solution(SolVol target){
int NE = dimpars.NE;
int DIMZ = dimpars.DIMZ;
int DIMY = dimpars.DIMY;
int DIMX = dimpars.DIMX;
for(int ie=NE-1 ; ie<NE; ie++) {
for( int iz=0; iz<DIMZ; iz++) {
for( int iy=0; iy<DIMY; iy++) {
for( int ix=0; ix<DIMX; ix++) {
target.ptr[ SOL_addr( ie, iz, iy , ix ) ]= 0.6 ;
}
}
}
}
}
void setMatff() {
int NE = dimpars.NE;
int NV = dimpars.NV;
int NROI = dimpars.NROI;
......@@ -191,7 +257,7 @@ public:
}
void settaMatSffF2(FreeFactsVol F) {
void setMatSffF2(FreeFactsVol F) {
int NE = dimpars.NE;
int NROI = dimpars.NROI;
int DIMX = dimpars.DIMX;
......@@ -211,7 +277,7 @@ public:
}
}
void settaFreeFacts(FreeFactsVol V , SolVol C) {
void setFreeFacts(FreeFactsVol V , SolVol C) {
float buff[100];
{
int NE = dimpars.NE;
......@@ -638,7 +704,8 @@ int main(int argc, char ** argv) {
boost::filesystem::path p(argv[index_input]);
std::string dirname = p.parent_path().string()+"/";
InputFileNames if_names = {
dirname + mockup_config["DSname"].as<std::string>() ,
dirname + mockup_config["DDname"].as<std::string>() ,
......@@ -668,25 +735,32 @@ int main(int argc, char ** argv) {
pb.settoVal(ffacts ,1.0f);
/*
ffacts.ptr[0]=0;
ffacts.ptr[2]=0;
ffacts.ptr[3]=0;
*/
{
double norma = sqrt(pb.scalar(ffacts, ffacts));
pb.Scal( ffacts,1.0/norma ) ;
}
pb.settoVal(X ,0.0f);
for(int iter_c = 0; iter_c<3; iter_c++) {
for(int iter_c = 0; iter_c<5; iter_c++) {
pb.settaVectA_and_Mat( XvectA, ffacts ) ;
pb.setVectA_and_Mat( XvectA, ffacts ) ;
float Lip=1.0;
pb.settoVal(Xtmp ,1.0f);
for(int i=0; i< 30; i++) {
for(int i=0; i< 60; i++) {
if(0) {
pb.copyScal(fftmp, ffacts, 1.0);
pb.settaFreeFacts( fftmp , X) ;
pb.setFreeFacts( fftmp , X) ;
}
pb.applyMatA(grad , Xtmp);
......@@ -698,27 +772,33 @@ int main(int argc, char ** argv) {
}
for(int iter = 0; iter<40; iter++) {
pb.applyMatA(grad , X);
for(int iter = 0; iter<1000; iter++) {
pb.settozero(grad);
pb.applyMatA(grad , X);
pb.gradReg_2_Add(grad, X , betaE, betaV ) ;
pb.AXPBYCR( -1.0/Lip, grad, +1.0/Lip, XvectA, 1.0, X);
// pb.project_solution(X);
double error = pb.scalar(X, grad)/2.0 - pb.scalar(X, XvectA) ;
double merit = pb.scalar(grad, grad) + pb.scalar( XvectA, XvectA) - 2*pb.scalar(grad, XvectA);
printf(" iter %d %e\n", iter, merit);
printf("iter %d %e\n", iter, error);
}
pb.save(X, "solution.raw");
pb.settaFreeFacts( ffacts , X) ;
pb.setFreeFacts( ffacts , X) ;
{
double norma = sqrt(pb.scalar(ffacts, ffacts));
pb.Scal( ffacts,1.0/norma ) ;
}
printf(" setta ok \n");
printf(" set ok \n");
}
pb.save(X, "solution.raw");
......
......@@ -296,7 +296,8 @@ monitor_path_template = datadir + filename_template +"_monitor.hd5:/monitor"
energy_custom_grid = np.array([6746.0,6754.0,6755.5,6756.2,6757.5,6759.3,6762.5,6770.0,6790.5])
# energy_custom_grid = np.array([6746.0,6754.0,6755.5,6756.2,6757.5,6759.3,6762.5,6770.0,6790.5])
energy_custom_grid = np.array([6745.0,6754.0,6755.5,6756.2,6757.5,6759.3,6762.5,6770.0,6791])
energy_exp_grid = h5py.File( datadir +(filename_template%1)+".nxs" ,"r")["/root.spyc.config1d_RIXS_0024/scan_data/actuator_1_1"][()]
......@@ -327,7 +328,7 @@ reference_file = reference_target_file
if(1): # SAMPLE extraction
if(0): # SAMPLE extraction
extract_sample_givenrois(
roi_path = roi_path ,
data_path_template = data_path_template ,
......@@ -351,7 +352,7 @@ if(1): # INTERPOLATION ESYNTH
if(1):
if(0):
os.system("rm %s"%scalarprods_target_file)
for iE in range(Edim) :
get_scalars( iE = iE,
......@@ -365,7 +366,7 @@ scalarprods_file = scalarprods_target_file
# ### ESYNTH
if(0):
if(1):
get_volume_Esynt( scalarprods_file = scalarprods_file,
interpolation_file = interpolation_infos_file)
......
Supports Markdown
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