Merge pull request #43 from Dirack/develop/2.1.1

Develop/2.1.1

.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()

test/test_parameters/mt_fixed_temp0_c0/Makefile

@@ -0,0 +1,10 @@
+all: generate_data generate_plot view
+
+generate_data: SConstruct
+	scons
+
+generate_plot: plot_test.py
+	python plot_test.py > surface.asc
+
+view: SConscript
+	scons -f $< view

test/test_parameters/mt_fixed_temp0_c0/SConscript

@@ -0,0 +1,26 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# SConstruct  (Madagascar Script)
+#
+# Purpose: Parameters test of VFSA process.
+#
+# Site: https://dirack.github.io
+# 
+# Version 1.0
+#
+# Programer: Rodolfo A. C. Neves (Dirack) 10/11/2020
+#
+# Email: rodolfo_profissional@hotmail.com
+#
+# License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.
+
+# Madagascar library
+from rsf.proj import *
+
+s="surface"
+Flow(s,s+'.asc','dd type=float form=native | smooth rect1=5 rect2=5 repeat=2')
+
+Result(s,'grey color=j scalebar=y title=Semblance')
+
+End()

test/test_parameters/mt_fixed_temp0_c0/SConstruct

@@ -0,0 +1,85 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# SConstruct  (Madagascar Script)
+#
+# Purpose: Parameters test of VFSA process.
+#
+# Site: https://dirack.github.io
+# 
+# Version 1.0
+#
+# Programer: Rodolfo A. C. Neves (Dirack) 10/11/2020
+#
+# 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
+
+# Modeling: Gaussian reflector in a velocity linear model
+# velocity increases with depth and 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)"
+     ''')
+
+# 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" ')
+
+# Kirchoff 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
+     ''')
+
+
+orep=0
+drep=1
+nrep=10
+
+ont=0.
+dnt=1
+nnt=200
+
+om0=5
+dm0=0.1
+nm0=1
+ot0=1.1
+dt0=0.1
+nt0=1
+v0=1.5
+
+for nti in range(nnt):
+    for nrepi in range(nrep):
+
+        repeat = orep+drep*nrepi
+        nt = ont+dnt*nti
+
+        parameters='crsParameters_rep_%d_nt_%d' % (nrepi,nti)
+        # Very Fast Simulated Aneelling Global Optimization (VFSA)
+        Flow(parameters,'dataCube',
+                '''
+                vfsacrsnh om0=%g dm0=%g nm0=%d ot0=%g dt0=%g nt0=%d v0=%g verb=y repeat=%d
+                half=y
+                c0=0.4 temp0=2.5 itmax=%d
+                ''' % (om0,dm0,nm0,ot0,dt0,nt0,v0,repeat,math.floor(nt)))
+
+End()