Merge pull request #41 from Dirack/feature/issue/37

Feature/issue/37

.gitignore

@@ -18,3 +18,4 @@ Fig/
 *.pdf
 *.x
 convgraph
+mtconvgraph

Mvfsacrsnh.c

@@ -1,6 +1,6 @@
 /* Zero offset CRS parameter inversion of RN, RNIP, BETA with Very Fast Simulated Annealing (VFSA) Global Optimization. This program uses the Non-Hyperbolic CRS approximation (Fomel and Kazinnik, 2013) to fit seismic data cube and get the best parameters using semblance criteria.
 
-Package version: 2.0.2
+Package version: 2.1.1
 
 Programmer: Rodolfo A. C. Neves (Dirack) 13/08/2021
 
@@ -62,9 +62,10 @@ int main(int argc, char* argv[])
 	bool interval;
 	bool half; // Use half-offset instead of offset
 	bool get_convergence_graph; // Option to generate a convergence graph
+	bool get_mt_convergence_graph; // Option to generate a multi thread convergence graph
+	float **mtgraph=NULL;
 	char strerr[50];
 	bool useprecseed;
-	struct timeval start;
 
 	/* RSF files I/O */  
 	sf_file in; /* Seismic data cube A(m,h,t) */
@@ -73,6 +74,7 @@ int main(int argc, char* argv[])
 	sf_file parameters=NULL;
 	sf_file out; /* RN, RNIP, BETA, Semblance, C0, Temp0, t0, m0 */
 	sf_file outgraph=NULL; /* Convergence graph */
+	sf_file mtoutgraph=NULL; /* Multi thread convergence graph */
 
 	/* RSF files axis */
 	sf_axis ax,ay,az;
@@ -85,6 +87,9 @@ int main(int argc, char* argv[])
 	if(!sf_getbool("getgraph",&get_convergence_graph)) get_convergence_graph=false;
 	/* Generate convergence graph (y/n) */
 
+	if(!sf_getbool("getmtgraph",&get_mt_convergence_graph)) get_mt_convergence_graph=false;
+	/* Generate multi thread convergence graph (y/n) */
+
 	if(!sf_getbool("half",&half)) half=false;
 	/* Use half-offset coordinates (y/n) */
 
@@ -164,6 +169,11 @@ int main(int argc, char* argv[])
 		prepareConvergenceGraphFile(outgraph = sf_output("convgraph"),get_convergence_graph,repeat,itmax);
 	}
 
+	if(get_mt_convergence_graph){
+		prepareMTConvergenceGraphFile(mtoutgraph = sf_output("mtconvgraph"),get_mt_convergence_graph,repeat,itmax);
+		mtgraph = sf_floatalloc2(itmax,repeat);
+	}
+
 	if(varlim)
 		loadParametersFilesVectors(parametersFilesVectors,parametersFilesLabels,nt0,nm0);
 
@@ -222,7 +232,7 @@ int main(int argc, char* argv[])
 
 				#pragma omp parallel for \
 				private(i,q,temp,c,RN,RNIP,BETA,cnew,semb) \
-				shared(semb0,otsemb,otrn,otrnip,otbeta) \
+				shared(semb0,otsemb,otrn,otrnip,otbeta,mtgraph) \
 				schedule(dynamic)
 				for(i=0;i<repeat;i++){
 
@@ -317,6 +327,12 @@ int main(int argc, char* argv[])
 							}	
 
 							if(get_convergence_graph) sf_floatwrite(&otsemb,1,outgraph);
+							if(get_mt_convergence_graph){
+								#pragma omp critical(write_semblance_per_thread)
+								{
+								mtgraph[i][q] = otsemb;
+								}
+							}
 
 					} /* loop over iterations */
 
@@ -352,4 +368,6 @@ int main(int argc, char* argv[])
 	sf_putfloat(out,"C0",c0);
 	sf_putfloat(out,"temp0",temp0);
 	sf_floatwrite(otm[0][0],4*nt0*nm0,out);
+
+	if(get_mt_convergence_graph) sf_floatwrite(mtgraph[0],itmax*repeat,mtoutgraph);
 }

debug/gdb_scripts/gdb_test_mtgraph

@@ -0,0 +1,16 @@
+# gdb_script (GDB Script)
+# 
+# Purpose: Run test in GDB shell.
+# 
+# Site: https://dirack.github.io
+# 
+# Programmer: Rodolfo A C Neves (Dirack) 04/01/2021
+# 
+# Email: rodolfo_profissional@hotmail.com
+# 
+# License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.
+
+cd test_data_cube/
+#set args <dataCube.rsf om0=5 dm0=0.1 nm0=1 ot0=1.1 dt0=0.1 nt0=1 v0=1.5 verb=y repeat=1 > teste.rsf
+set args <dataCube.rsf itmax=1000 getmtgraph=y mtconvgraph=graph.rsf half=y temp0=10 c0=0.5 om0=3 dm0=1 nm0=1 ot0=1.1 dt0=0.1 nt0=1 v0=1.5 verb=y repeat=2 useprecseed=y > teste.rsf
+show args

test/libraries/Test_vfsacrsnh_lib.c

@@ -171,6 +171,7 @@ void convergence_graph_file_preparation(){
 	int n1;
 	float d1;
 	float o1;
+	char label[50];
 
 	prepareConvergenceGraphFile(outgraph = sf_output("convgraph"),get_convergence_graph,repeat,itmax);
 	TEST_ASSERT_TRUE(sf_histint(outgraph,"n1",&n1));
@@ -179,6 +180,11 @@ void convergence_graph_file_preparation(){
 	TEST_ASSERT_FLOAT_WITHIN(0.001,d1,1.0);
 	TEST_ASSERT_TRUE(sf_histfloat(outgraph,"o1",&o1));
 	TEST_ASSERT_FLOAT_WITHIN(0.001,o1,0.0);
+
+	strcpy(label,sf_histstring(outgraph,"label1"));
+	TEST_ASSERT_EQUAL_STRING("Iteration",label);
+	strcpy(label,sf_histstring(outgraph,"label2"));
+	TEST_ASSERT_EQUAL_STRING("Semblance",label);
 }
 
 void check_parameters_file_dimensions(){
@@ -303,6 +309,38 @@ void check_load_parameters_file_data(){
 	}
 }
 
+void mt_convergence_graph_file_preparation(){
+/*< Test multi thread convergence graph file preparation. Test if the parameters set is correct >*/
+	sf_file outgraph=NULL;
+	int itmax = 500, repeat = 5;
+	bool get_mt_convergence_graph=true;
+	int n1, n2;
+	float d1, d2;
+	float o1, o2;
+	char label[50];
+
+	prepareMTConvergenceGraphFile(outgraph = sf_output("mtconvgraph"),get_mt_convergence_graph,repeat,itmax);
+	TEST_ASSERT_TRUE(sf_histint(outgraph,"n1",&n1));
+	TEST_ASSERT_EQUAL(n1,itmax);
+	TEST_ASSERT_TRUE(sf_histfloat(outgraph,"d1",&d1));
+	TEST_ASSERT_FLOAT_WITHIN(0.001,d1,1.0);
+	TEST_ASSERT_TRUE(sf_histfloat(outgraph,"o1",&o1));
+	TEST_ASSERT_FLOAT_WITHIN(0.001,o1,0.0);
+
+	TEST_ASSERT_TRUE(sf_histint(outgraph,"n2",&n2));
+	TEST_ASSERT_EQUAL(n2,repeat);
+	TEST_ASSERT_TRUE(sf_histfloat(outgraph,"d2",&d2));
+	TEST_ASSERT_FLOAT_WITHIN(0.001,d1,1.0);
+	TEST_ASSERT_TRUE(sf_histfloat(outgraph,"o2",&o2));
+	TEST_ASSERT_FLOAT_WITHIN(0.001,o1,0.0);
+
+	strcpy(label,sf_histstring(outgraph,"label1"));
+	TEST_ASSERT_EQUAL_STRING("Iteration",label);
+	strcpy(label,sf_histstring(outgraph,"label2"));
+	TEST_ASSERT_EQUAL_STRING("Semblance",label);
+}
+
+
 int main(int argc, char* argv[]){
 
 	/* Redirect the stdin to datacube file */
@@ -329,6 +367,7 @@ int main(int argc, char* argv[]){
 	RUN_TEST(check_parameters_file_dimensions);
 	RUN_TEST(check_load_data_cube_file_dimensions);
 	RUN_TEST(check_load_parameters_file_data);
+	RUN_TEST(mt_convergence_graph_file_preparation);
 
 	return UNITY_END();
 }

test/test_parameters/mt_convergence_graph/SConstruct

@@ -0,0 +1,106 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#
+# SConstruct  (Madagascar Script)
+#
+# Purpose: Obtain zero offset CRS parameters using VFSA
+# for a set of (t0,m0) pairs.
+#
+# Site: https://www.geofisicando.com
+#
+# Programmer: Rodolfo A. C. Neves (Dirack) 13/08/2021
+#
+# Email: rodolfo_profissional@hotmail.com
+#
+# License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.
+
+# Madagascar library
+from rsf.proj import *
+
+# Python math library
+import math
+import subprocess
+
+__author__="Rodolfo Dirack <rodolfo_profissional@hotmail.com>"
+__version__="1.0"
+__site__="https://www.geofisicando.com"
+
+Help('''
+    Obtain zero offset CRS parameters using VFSA for a set of (t0,m0) pairs
+    =======================================================================
+
+    Obtain RN, RNIP and BETA parameters using Very Fast Simulated Annealing (VFSA)
+    global optimization
+
+    Author: %s
+    Version: %s
+    Site: %s
+    ''' % (__author__,__version__,__site__))
+
+Help('''
+### Kirchhoff modeling ###
+
+Use Kirchhoff modeling in a Gaussian reflector linear velocity model
+velocity increases with depth with a 0.5 velocity gradient
+''')
+
+Flow('gaussianReflector',None,
+     '''
+     math d1=0.01 n1=2001 o1=-5 unit1=km label1=Offset
+     output="4-3*exp(-(x1-5)^2/9)"
+     ''')
+
+# Define gaussian velocity Model
+Flow('velocityModel','gaussianReflector',
+     '''
+     window min1=0 max1=10 |
+     spray axis=1 n=451 d=0.01 o=0 label=Depth unit=km |
+     math output="1.5+0.5*x1+0.0*x2"
+     ''')
+
+Flow('reflectorDip','gaussianReflector','math output="2/3*(x1-5)*input" ')
+
+# Kirchhoff Modeling
+Flow('dataCube','gaussianReflector reflectorDip',
+     '''
+     kirmod cmp=y dip=${SOURCES[1]} 
+     nh=161 dh=0.025 h0=0
+     ns=401 ds=0.025 s0=0
+     freq=10 dt=0.004 nt=1001
+     vel=1.5 gradz=0.5 gradx=0.0 verb=y |
+     put d2=0.0125 label3="CMP" unit3="Km" label2="Offset" unit2="Km" label1=Time unit1=s
+     ''')
+
+om0=3
+dm0=1
+nm0=1
+ot0=1.1
+dt0=0.1
+nt0=1
+ot0=1.1
+v0=1.5
+repeat=1
+
+Help('''
+### Very Fast Simulated Annealing (VFSA) ###
+
+Obtain RN, RNIP and BETA parameters for nm0=%g om0=%g dm0=%g nt0=%g ot0=%g dt0=%g v0=%g
+using non-hyperbolic CRS aproximation and VFSA repeated %d times
+''' % (nm0,om0,dm0,nt0,ot0,dt0,v0,repeat))
+
+ngraphs = int(ARGUMENTS.get("ngraphs",5))
+
+parameter = 'crsParameters'
+graph = 'graph'
+
+# Very Fast Simulated Annealling Global Optimization (VFSA)
+Flow([parameter,graph],'dataCube',
+        '''
+        vfsacrsnh itmax=1000
+        getmtgraph=y mtconvgraph=${TARGETS[1]} half=y temp0=10 c0=0.5 om0=%g dm0=%g nm0=%d ot0=%g dt0=%g nt0=%d v0=%g verb=y repeat=%d
+        useprecseed=y
+        ''' % (om0,dm0,nm0,ot0,dt0,nt0,v0,ngraphs))
+
+Result(graph,'scale axis=1 rscale=-1 | graph title="Semblance - Using %d threads" symbol=*' % ngraphs)
+
+End()

vfsacrsnh_lib.c

@@ -230,7 +230,7 @@ void prepareConvergenceGraphFile(sf_file outgraph, bool get_convergence_graph, i
 	sf_putfloat(outgraph,"d1",1);
 	sf_putfloat(outgraph,"o1",0);
 	sf_putstring(outgraph,"title","Convergence graph");
-	sf_putstring(outgraph,"label1","iteration");
+	sf_putstring(outgraph,"label1","Iteration");
 	sf_putstring(outgraph,"label2","Semblance");
 }
 
@@ -302,3 +302,17 @@ void loadParametersFilesVectors(float** parametersFilesVectors[6],char* labels[6
 		sf_fileclose(parametersFiles);
 	}
 }
+
+void prepareMTConvergenceGraphFile(sf_file outgraph, bool get_mt_convergence_graph, int repeat, int itmax)
+/*< Check parameters (repeat should be equal to 1 when generating a convergence graph)  and Prepare the convergence graph file >*/
+{
+	sf_putint(outgraph,"n1",itmax);
+	sf_putfloat(outgraph,"d1",1);
+	sf_putfloat(outgraph,"o1",0);
+	sf_putint(outgraph,"n2",repeat);
+	sf_putfloat(outgraph,"o2",0.);
+	sf_putfloat(outgraph,"d2",1.);
+	sf_putstring(outgraph,"label1","Iteration");
+	sf_putstring(outgraph,"label2","Semblance");
+}
+

vfsacrsnh_lib.h

@@ -103,4 +103,8 @@ bool checkAndLoadDataCubeDimensionsReturnStrError(sf_file in,
 void loadParametersFilesVectors(float** parametersFilesVectors[6],char* labels[6], int nt0, int nm0);
 /*< Load parameters file vector >*/
 
+
+void prepareMTConvergenceGraphFile(sf_file outgraph, bool get_mt_convergence_graph, int repeat, int itmax);
+/*< Check parameters (repeat should be equal to 1 when generating a convergence graph)  and Prepare the convergence graph file >*/
+
 #endif