Merge pull request #57 from Dirack/develop/2.2.3

Develop/2.2.3

Author: Dirack

Mvfsacrsnh.c

@@ -64,6 +64,7 @@ int main(int argc, char* argv[])
 	float **mtgraph=NULL;
 	char strerr[50];
 	bool useprecseed;
+	int mMAX, hMAX; // Option to define aperture in CMP and Offset
 
 	/* RSF files I/O */  
 	sf_file in; /* Seismic data cube A(m,h,t) */
@@ -158,6 +159,12 @@ int main(int argc, char* argv[])
 	if(!sf_getbool("useprecseed",&useprecseed)) useprecseed=false;
 	/* y: Use a more precise seed for vfsa; n: Use time(NULL)*/
 
+	if(!sf_getint("mmax",&mMAX)) mMAX=50;
+	/* CMP aperture */
+
+	if(!sf_getint("hmax",&hMAX)) hMAX=50;
+	/* Offset aperture */
+
 	if(! sf_getbool("verb",&verb)) verb=false;
 	/* y: active mode; n: quiet mode */
 
@@ -183,6 +190,15 @@ int main(int argc, char* argv[])
 	sf_floatread(t[0][0],nm*nh*nt,in);
 	sf_fileclose(in);
 
+	if(!validBoundariesAndApertureForCMP(nm0, om0, dm0, mMAX, nm, om, dm, strerr))
+		sf_error("%s",strerr);
+
+	if(!validBoundariesAndApertureForOffset(nh,oh,dh,hMAX,strerr))
+		sf_error("%s",strerr);
+
+	if(!validBoundariesAndApertureForT0(nt0, ot0, dt0, nt, ot, dt, strerr))
+		sf_error("%s",strerr);
+
 	if (verb) {
 		sf_warning("Active mode on!!!");
 		sf_warning("Input file parameters: ");
@@ -270,7 +286,7 @@ int main(int argc, char* argv[])
 						otrn=c[0];
 						otrnip=c[1];
 						otbeta=c[2];
-						semb0=semblance(m0,dm,om,oh,dh,dt,nt,t0,v0,c[0],c[1],c[2],t,half);
+						semb0=semblance(m0,dm,om,oh,dh,dt,nt,t0,v0,c[0],c[1],c[2],t,mMAX,hMAX,half);
 						otsemb = semb0;
 					}
 
@@ -289,7 +305,7 @@ int main(int argc, char* argv[])
 
 						semb=0;
 
-						semb=semblance(m0,dm,om,oh,dh,dt,nt,t0,v0,RN,RNIP,BETA,t,half);
+						semb=semblance(m0,dm,om,oh,dh,dt,nt,t0,v0,RN,RNIP,BETA,t,mMAX,hMAX,half);
 
 							/* VFSA parameters convergence condition */		
 							if(fabs(semb) > fabs(semb0) ){

docs/VERSION.md

@@ -1 +1 @@
-2.2.2
+2.2.3

infra/scripts/version_check.sh

@@ -2,7 +2,7 @@
 #
 # Check current version based on active develop branch name
 
-CURRENTVERSION="2.2.2"
+CURRENTVERSION="2.2.3"
 
 VERSIONDOC=$(cat docs/VERSION.md)
 

test/bug/Makefile

@@ -10,7 +10,7 @@
 # 
 # License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.
 
-TESTS=test_issue_32 test_issue_46
+TESTS=test_issue_32 test_issue_46 test_issue_52
 
 all: $(TESTS)
 

test/bug/issue_52/Makefile

@@ -0,0 +1,38 @@
+# Makefile
+#
+# Purpose: Tests for bug related in issue 52.
+# 
+# Site: https://www.geofisicando.com
+# 
+# Programmer: Rodolfo A C Neves (Dirack) 27/02/2024
+# 
+# Email: [email protected]
+# 
+# License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.
+
+# Uncomment the line below to run on a CENTOS system. Or call make as 'make CENTOS=1'
+#CENTOS=1 
+
+DEP = ../../../vfsacrsnh_lib.c
+GDBDEP = ../../../Mvfsacrsnh.c $(DEP)
+UNITY = ../../Unity/unity.c
+
+# Default is to run the test
+all:	Test_bug_issue_52.x
+	@./Test_bug_issue_52.x \
+		; rm Test_bug_issue_52.x
+
+%.x:	%.c
+ifdef CENTOS
+	@gcc -I$$RSFSRC/include -I../.. -I/usr/include/openblas -L$$RSFSRC/lib $< $(DEP) $(UNITY) -o $@ -lm -lrsf -lopenblas -g
+else
+	@gcc -I$$RSFSRC/include -I../.. -L$$RSFSRC/lib $< $(DEP) $(UNITY) -o $@ -lm -lrsf
+endif
+
+clean:
+	rm *.x
+
+help:
+	@echo "Compilation of the tests\nRun 'make' to unit test"
+
+.PHONY:	clean help all

test/bug/issue_52/README.md

@@ -0,0 +1,3 @@
+# Tests for libraries and functions
+
+Those tests run using the program make

test/bug/issue_52/Test_bug_issue_52.c

@@ -0,0 +1,101 @@
+/*
+* Test_bug_issue_52.c (C)
+* 
+* Purpose: Test of bug related in issue 52.
+* 
+* Site: https://www.geofisicando.com
+* 
+* Programmer: Rodolfo A C Neves (Dirack) 27/02/2024
+* 
+* Email: [email protected]
+* 
+* License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.
+*/
+
+#include "Unity/unity.h"
+#include "../vfsacrsnh_lib.h"
+#include <rsf.h>
+#define ITMAX 100
+#define hMAX 50
+#define mMAX 50
+
+void setUp(){};
+
+void tearDown(){};
+
+void testFunctionValueInsideBoundaries(){
+/*< Function to assure that value is inside boundaries >*/
+	TEST_ASSERT_TRUE(checkValueInsideBoundaries(1.0,0.,2.0));	
+	TEST_ASSERT_FALSE(checkValueInsideBoundaries(1.0,1.1,2.0));	
+	TEST_ASSERT_FALSE(checkValueInsideBoundaries(1.0,0.,0.5));	
+}
+
+void checkValidBoundariesAndApertureForCMP(){
+/*< Function to assure that m0 values are inside boundaries to avoid segmentation fault errors >*/
+
+	char strerr[100];
+	int nm=401;
+	float dm=0.025, om=0.;
+	int nm0=1;
+	float dm0=1., om0=3;
+	int aperture=50;
+
+	TEST_ASSERT_TRUE(validBoundariesAndApertureForCMP(nm0, om0, dm0, aperture, nm, om, dm, strerr));
+
+	om0=1.;
+	TEST_ASSERT_FALSE(validBoundariesAndApertureForCMP(nm0, om0, dm0, aperture, nm, om, dm, strerr));
+	TEST_ASSERT_EQUAL_STRING("Invalid boundaries for CMP m0min=-0.250000 m0max=2.250000. Please check aperture!",strerr);
+
+	om0=9.;
+	TEST_ASSERT_FALSE(validBoundariesAndApertureForCMP(nm0, om0, dm0, aperture, nm, om, dm, strerr));
+	TEST_ASSERT_EQUAL_STRING("Invalid boundaries for CMP m0min=7.750000 m0max=10.250000. Please check aperture!",strerr);
+}
+
+void checkValidBoundariesAndApertureForOffset(){
+/*< Function to assure that Offset values are inside boundaries to avoid segmentation fault errors >*/
+
+        char strerr[100];
+        int nh=161;
+        float dh=0.0125, oh=0.;
+        int aperture=50;
+
+        TEST_ASSERT_TRUE(validBoundariesAndApertureForOffset(nh, oh, dh, aperture, strerr));
+
+        oh=0.;
+	aperture=162;
+        TEST_ASSERT_FALSE(validBoundariesAndApertureForOffset(nh, oh, dh, aperture, strerr));
+        TEST_ASSERT_EQUAL_STRING("Invalid boundaries for Offset hmin=0.000000 hmax=2.025000. Please check aperture!",strerr);
+
+}
+
+void checkValidBoundariesAndApertureForT0(){
+/*< Function to assure that t0 values are inside boundaries to avoid segmentation fault errors >*/
+
+        char strerr[100];
+        int nt=1001;
+        float dt=0.004, ot=0.;
+        int nt0=1;
+        float dt0=1., ot0=3;
+
+        TEST_ASSERT_TRUE(validBoundariesAndApertureForT0(nt0, ot0, dt0, nt, ot, dt, strerr));
+
+        ot0=4.004;
+        TEST_ASSERT_FALSE(validBoundariesAndApertureForT0(nt0, ot0, dt0, nt, ot, dt, strerr));
+        TEST_ASSERT_EQUAL_STRING("Invalid boundaries for time t0 t0min=4.004000 t0max=4.004000. Please check aperture!",strerr);
+
+        ot0=-0.004;
+        TEST_ASSERT_FALSE(validBoundariesAndApertureForT0(nt0, ot0, dt0, nt, ot, dt, strerr));
+        TEST_ASSERT_EQUAL_STRING("Invalid boundaries for time t0 t0min=-0.004000 t0max=-0.004000. Please check aperture!",strerr);
+}
+
+
+int main(int argc, char* argv[]){
+
+	UNITY_BEGIN();
+	RUN_TEST(testFunctionValueInsideBoundaries);
+	RUN_TEST(checkValidBoundariesAndApertureForCMP);
+	RUN_TEST(checkValidBoundariesAndApertureForOffset);
+	RUN_TEST(checkValidBoundariesAndApertureForT0);
+
+	return UNITY_END();
+}

test/libraries/Test_vfsacrsnh_lib.c

@@ -16,10 +16,12 @@
 #include "../vfsacrsnh_lib.h"
 #include <rsf.h>
 #define ITMAX 100
+#define hMAX 50
+#define mMAX 50
 
 /* Seimic datacube parameters to be used in tests */
 float*** t; // it stores seismic datacube
-float returnedValues[2*mMAX+1][hMAX]; // nonHyperbolicCRSapp returned values
+float **returnedValues; // nonHyperbolicCRSapp returned values
 int nt=1001, nh=401, nm=161; // datacube dimensions
 sf_file in; // datacube RSF file
 
@@ -124,7 +126,7 @@ void nonHyperbolicCRSapp_function_for_pre_calculated_values(){
 	int h;
 	int half=1;
 
-	nonHyperbolicCRSapp(returnedValues,m0,dm,om,dh,oh,t0,v0,RN,RNIP,BETA,half);
+	nonHyperbolicCRSapp(returnedValues,mMAX,hMAX,m0,dm,om,dh,oh,t0,v0,RN,RNIP,BETA,half);
 
 	/* In the CMP m0 and h0, time is equal to t0 */
 	/* CMP 51 is m0, that is the midle of CMP window */
@@ -148,7 +150,7 @@ void semblance_return_value_between_0_1(){
 	int half=1;
 
 	for(i=0;i<100;i++){
-		semb=semblance(m0,dm,om,oh,dh,dt,nt,t0,v0,RN,RNIP,BETA,t,half);
+		semb=semblance(m0,dm,om,oh,dh,dt,nt,t0,v0,RN,RNIP,BETA,t,mMAX,hMAX,half);
 		TEST_ASSERT(semb>=0.0 && semb <= 1.0);
 	}
 }
@@ -340,7 +342,6 @@ void mt_convergence_graph_file_preparation(){
 	TEST_ASSERT_EQUAL_STRING("Semblance",label);
 }
 
-
 int main(int argc, char* argv[]){
 
 	/* Redirect the stdin to datacube file */
@@ -353,6 +354,8 @@ int main(int argc, char* argv[]){
 	t=sf_floatalloc3(nt,nh,nm);
 	sf_floatread(t[0][0],nm*nh*nt,in);
 
+	returnedValues = sf_floatalloc2(hMAX,2*mMAX+1);
+
 	UNITY_BEGIN();
 	RUN_TEST(signal_function_return_1_for_positive_values_and_zero);
 	RUN_TEST(signal_function_return_minus_1_for_negative_values);

test/test_parameters/rnip_beta_plane_interface/Test_parameter_values.sh

@@ -45,7 +45,7 @@ check_value(){
 }
 
 check_value test_rnip_interface1.txt "1.0" "0.1"
-check_value test_rnip_interface2.txt "1.95" "0.1"
+check_value test_rnip_interface2.txt "1.95" "0.12"
 check_value test_beta_interface1.txt "0.00" "0.1"
 check_value test_beta_interface2.txt "0.00" "0.1"
 

vfsacrsnh_lib.c

@@ -13,8 +13,6 @@
 
 */
 
-#define hMAX 50 // Max of samples to stack in half-offset
-#define mMAX 50 // Max of samples to stack in CMP
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
@@ -109,7 +107,9 @@ VFSA algorithm disturb parameters step.
 
 }
 
-void nonHyperbolicCRSapp(	float t[2*mMAX+1][hMAX] /* non-hyperbolic CRS traveltime surface */,
+void nonHyperbolicCRSapp(	float **t /* non-hyperbolic CRS traveltime surface */,
+				int mMAX /* CMP aperture */,
+				int hMAX /* Offset aperture */,
 				float m0 /* Central CMP of the approximation */,
 				float dm, /* CMP sampling */
 				float om /* CMP axis origin */,
@@ -169,6 +169,8 @@ float semblance(float m0 /* Central CMP of the approximation */,
 		float RNIP /* RNIP, CRS parameter */,
 		float BETA /* BETA, CRS parameter */,
 		float*** t /* reflection data cube A(m,h,t) */,
+		int mMAX, /* CMP aperture */
+		int hMAX, /* Offset aperture */
 		bool half /* Half-offset coordinates */)
 /*< Calculate semblance between the Non Hyperbolic CRS approximation surface and reflection data >*/
 {
@@ -179,13 +181,15 @@ float semblance(float m0 /* Central CMP of the approximation */,
 	float amplitudeSquaredSampleSum=0.;
 	float semblance=0;
 	int tetai;
-	float teta[2*mMAX+1][hMAX];
+	float **teta;
 	int m0_index_init, m0_index_end;
 
+	teta = sf_floatalloc2(hMAX,2*mMAX+1);
+
 	m0_index_init = (int)(m0/dm)-mMAX;
 	m0_index_end = (int)(m0/dm)+mMAX;
 
-	nonHyperbolicCRSapp(teta,m0,dm,om,dh,oh,t0,v0,RN,RNIP,BETA,half);
+	nonHyperbolicCRSapp(teta,mMAX,hMAX,m0,dm,om,dh,oh,t0,v0,RN,RNIP,BETA,half);
 
 	for (im=m0_index_init; im < m0_index_end; im++){
 
@@ -213,6 +217,54 @@ float semblance(float m0 /* Central CMP of the approximation */,
 
 }
 
+bool checkValueInsideBoundaries(float value, float min, float max)
+/*< check if the value is inside the boundaries min < value < max >*/
+{
+	return (min < value && value < max);
+}
+
+bool validBoundariesAndApertureForCMP(int nm0, float om0, float dm0, int mMAX, int nm, float om, float dm, char *strerr)
+/*< Check if CMP coordinates for m0s are inside data cube space to avoid segmentation fault errors >*/
+{
+	float m0min, m0max;
+
+	m0min = om0-(dm*mMAX);
+	m0max = om0+dm0*(nm0-1)+(dm*mMAX);
+
+	if(m0min < om || m0max > (om+dm*(nm-1))){
+		sprintf(strerr,"Invalid boundaries for CMP m0min=%f m0max=%f. Please check aperture!",m0min,m0max);
+		return false;
+	}
+
+	return true;
+}
+
+bool validBoundariesAndApertureForOffset(int nh, float oh, float dh, int hMAX, char *strerr)
+/*< Check if Offset coordinates are inside data cube space to avoid segmentation fault errors >*/
+{
+	if((oh+dh*hMAX) > (oh+dh*(nh-1))){
+		sprintf(strerr,"Invalid boundaries for Offset hmin=%f hmax=%f. Please check aperture!",oh,oh+dh*hMAX);
+		return false;
+	}
+
+	return true;
+}
+
+bool validBoundariesAndApertureForT0(int nt0, float ot0, float dt0, int nt, float ot, float dt, char *strerr)
+/*< Check if t0 coordinates are inside data cube space to avoid segmentation fault errors >*/
+{
+	float t0min = ot0;
+	float t0max = ot0+dt0*(nt0-1);
+
+	if(t0min < ot || t0max > ot+dt*(nt-1)){
+		sprintf(strerr,"Invalid boundaries for time t0 t0min=%f t0max=%f. Please check aperture!",t0min,t0max);
+		return false;
+	}
+
+	return true;
+}
+
+
 bool repeatOptionEqual1ForGetConvergenceGraphTrue(bool get_convergence_graph, int repeat)
 /*< Repeat option should be equal one when get_convergence_graph flag is on to avoid multiple thread >*/
 {