Changes to interpolate GB as a whole

.gitignore

@@ -11,7 +11,7 @@ InputDataError.txt
 *.pdf
 *.pyc
 *.HOLD
-fort.8
+fort.*
 mod/
 obj/
 Data/

Configuration/Interpolation/SpatialFcns.cfg

@@ -20,5 +20,3 @@ Function, dim=x, shape=Logistic, precon=1
 Function, dim=y, shape=Logistic, precon=1
 Function, dim=x, shape=Logistic, precon=2
 Function, dim=y, shape=Logistic, precon=2
-Function, dim=x, shape=Logistic, precon=3
-Function, dim=y, shape=Logistic, precon=3

PythonScripts/GUI/GeoSAM/GeoSams.py

@@ -252,7 +252,10 @@ def Run_Sim(self):
                 if returnCode == 0:
                     messagebox.showinfo("GeoSAMS/UK/Plotting", f'ALL DONE\n{args}')
                 else:
-                    messagebox.showerror("UK", f'Failed\n{args}\nReturn Code = {returnCode}')
+                    if returnCode == 99: 
+                        messagebox.showerror("UK", f'Failed\n{args}\nReturn Code = {returnCode}\nMath Failure. Try changing Spatial Functions')
+                    else:
+                        messagebox.showerror("UK", f'Failed\n{args}\nReturn Code = {returnCode}')
             else:
                 messagebox.showerror("GeoSAM Sim", f'Failed\n{result.args}\nReturn Code = {result.returncode}')
 
@@ -304,14 +307,18 @@ def InterpAndPlotResults(self):
         # Determine which if any file suffixes are needed
         if self.domainName=='MA': self.savedByStratum = False
         if self.savedByStratum:
-            # Only GB and AL[L] use savedByStratum
+            # Only GB and AL use savedByStratum
             if self.domainName=='GB':
-                rgn = ['_SW', '_N', '_S', '_W']
+                # rgn = ['_SW', '_N', '_S', '_W']
+                rgn = ['_GB']
             else:
                 # This would be AL
-                rgn = ['_SW', '_N', '_S', '_W', '_MA']
+                #rgn = ['_SW', '_N', '_S', '_W', '_MA']
+                rgn = ['_GB', '_MA']
         else:
-            # This would be just MA
+            # This would be just MA and, since MA forces savedByStratum to false, it does NOT append the domain suffix
+            # i.e. Data/X_Y_LPUE_MA2015_0_MA.csv
+            #                             ^^
             rgn = ['']
 
         # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -349,10 +356,13 @@ def InterpAndPlotResults(self):
                         gridFile = 'MAxyzLatLon.csv'
                         ukCfgFile = 'UK_MA.cfg'
                     else:
-                        gridFile = 'GBxyzLatLon' + r + '.csv'
+                        #gridFile = 'GBxyzLatLon' + r + '.csv'
+                        gridFile = 'GBxyzLatLon.csv'
                         ukCfgFile = 'UK_GB.cfg'
                 else:
-                    gridFile = self.domainName+'xyzLatLon' + r + '.csv'
+                    # DEPRECATE: if we no longer need to separate GB into sub regions
+                    #gridFile = self.domainName+'xyzLatLon' + r + '.csv'
+                    gridFile = self.domainName+'xyzLatLon.csv'
                 cmd = [ex, ukCfgFile, self.domainName, obsFile, gridFile, zArg]
                 result = subprocess.run(cmd)
                 if (result.returncode != 0):

SRC/Globals.f90

@@ -20,7 +20,6 @@ module globals
 integer, parameter :: max_num_years  = 50
 integer, parameter :: max_num_areas = 25
 integer, parameter :: max_sides = 8
-integer, parameter :: num_regions = 5
 integer, parameter :: region_none=0
 integer, parameter :: region_N=1
 integer, parameter :: region_S=2
@@ -43,7 +42,7 @@ module globals
 integer, parameter :: write_dev = 63
 
 ! arbritrary value chosen as smallest with 2 digit exponent
-real(dp), parameter :: zero_thresh = 1.0D-99 
+real(dp), parameter :: zero_threshold = 1.0D-99 
 
 ! ASIN incorrectly produces error
 ! Error: Fortran 2003: Elemental function as initialization expression with non-integer/non-character arguments
@@ -103,25 +102,6 @@ logical function Is_Leap_Year (yr)
     endif
 endfunction
 
-!---------------------------------------------------------------------------------------------------
-!> Purpose: Computes SQRT( SUM( X(1:N)^2 / N) )
-!> @param[in] x a real array 
-!> @param[in] n the number of elements in x
-!---------------------------------------------------------------------------------------------------
-real(dp) function Compute_RMS(x, n)
-    implicit none
-    integer, intent(in) :: n
-    real(dp), intent(in) :: x(:)
-    real(dp) result 
-
-    result = sqrt(sum(x(1:n)**2) / float(n))
-    if (is_nan(result)) then
-        write(*,*) term_red, 'Compute_RMS FAILED; NaN', term_blk
-        STOP 99
-    endif
-    Compute_RMS = result
-endfunction Compute_RMS
-
 !---------------------------------------------------------------------------------------------------
 !> isnan does not work nor does (x/=x)
 !>
@@ -134,18 +114,6 @@ logical function is_nan(x)
     is_nan = (buf(8:10) .EQ. 'NaN')
 endfunction is_nan
 
-!---------------------------------------------------------------------------------------------------
-!> Purpose: Computes Artithmetic Mean
-!> @param[in] x a real array 
-!> @param[in] n the number of elements in x
-!---------------------------------------------------------------------------------------------------
-real(dp) function Compute_MEAN(x, n)
-    implicit none
-    integer, intent(in) :: n
-    real(dp), intent(in) :: x(:)
-    Compute_MEAN = sum(x(1:n) / float(n))
-endfunction Compute_MEAN
-
 !==============================
 ! Changes a string to upper case
 !==============================
@@ -240,6 +208,11 @@ function matrixinv(x,n)
             enddo
         enddo
     enddo   
-endfunction matrixinv
+
+    if (matrixinv(1,1)/=matrixinv(1,1)) then
+        write(*,*) term_red, 'matrixinv FAILED', term_blk
+        STOP 99
+    endif 
+    endfunction matrixinv
 
 end module globals

SRC/IORoutines.f90

@@ -264,7 +264,13 @@ subroutine Write_Column_CSV(n,f,header,file_name,append)
     write(write_dev, '(A)'//NEW_LINE(cr)) trim(output_str)
     do k=1,n
         read(read_dev,'(A)',iostat=io) input_str
-        write(output_str,'(A,A,(ES14.7 : ))') trim(input_str),',',f(k)
+        ! sometimes f(k) takes on the value of E-311, which can not be read back in
+        ! seems the minimum value is E-300, even with using format ES15.7E3 
+        if (((f(k) > 0.D0) .AND. (f(k) < zero_threshold)) .OR. (f(k) .GT. 1.D0/zero_threshold)) then
+            write(output_str,'(A,A,(ES14.7 : ))') trim(input_str),',',0.D0
+        else
+            write(output_str,'(A,A,(ES14.7 : ))') trim(input_str),',',f(k)
+        endif
         write(write_dev, '(A)'//NEW_LINE(cr)) trim(output_str)
     enddo
     close(write_dev)
@@ -284,7 +290,13 @@ subroutine Write_Column_CSV(n,f,header,file_name,append)
     open(write_dev,file=file_name)
     write(write_dev, '(A)') header
     do k=1,n
-        write(write_dev,fmtstr) f(k)
+        ! sometimes f(k) takes on the value of E-311, which can not be read back in
+        ! seems the minimum value is E-300, even with using format ES15.7E3 
+        if (((f(k) > 0.D0) .AND. (f(k) < zero_threshold)) .OR. (f(k) .GT. 1.D0/zero_threshold)) then
+            write(write_dev,fmtstr) 0.D0
+        else
+            write(write_dev,fmtstr) f(k)
+        endif
     enddo
     close(write_dev)
 endif

SRC/ScallopPopDensity.f90

@@ -610,8 +610,11 @@ subroutine Write_Column_CSV_By_Stratum(n,f, lat, lon, stratum, header,file_name,
 character(5000) output_str
 integer, parameter :: temp_dev = 70
 integer, parameter :: appd_dev = 80
-integer offset
-character(3) :: rgn(num_regions) = (/ '_N ', '_S ', '_SW', '_W ', '_MA'/)
+integer, parameter :: num_regions = 2
+!character(3) :: rgn(num_regions) = (/ '_N ', '_S ', '_SW', '_W ', '_MA'/)
+character(3) :: rgn(num_regions) = (/ '_GB', '_MA'/)
+integer offset, region
+
 if (append) then
     ! read existing files by region0
     ! create a temp file by region to write output then move temp to existing file_name
@@ -624,16 +627,22 @@ subroutine Write_Column_CSV_By_Stratum(n,f, lat, lon, stratum, header,file_name,
         write(temp_dev+offset, '(A)'//NEW_LINE(cr)) trim(output_str)
     enddo
     do k=1,n
-        offset = Get_Region(lat(k), lon(k), stratum(k))
-        if (offset > 0) then
+        ! offset = Get_Region(lat(k), lon(k), stratum(k))
+        region = Get_Region(lat(k), lon(k), stratum(k))
+        if (region > 0) then
+            if (region < region_MA) then
+                offset = 1
+            else
+                offset = 2
+            endif
             read(appd_dev+offset,'(A)',iostat=io) input_str
             if (f(k) < 0.0) then
                 write(*,'(A,A,A,A,A,A,A,I5,A)') term_yel, 'WARNING: Negative value in: ', term_blk, &
                 &  file_name//trim(rgn(offset))//'.csv', term_yel, ' line: ', term_blk, k
             endif
             ! sometimes f(k) takes on the value of E-311, which can not be read back in
             ! seems the minimum value is E-300, even with using format ES15.7E3 
-            if ((f(k) > 0.D0) .AND. (f(k) < zero_thresh)) then
+            if ((f(k) > 0.D0) .AND. (f(k) < zero_threshold)) then
                 write(output_str,'(A,A,(ES14.7 : ))') trim(input_str),',',0.D0
             else
                 write(output_str,'(A,A,(ES14.7 : ))') trim(input_str),',',f(k)

UKsrc/NonLinearSpatialFcn.f90

@@ -325,7 +325,7 @@ subroutine NLSF_Least_Sq_Fit(obs, nlsf, save_data)
 
     ! initialize residual 0
     residual_vect(:) = obs%field(1:num_obs_points)
-    write(*,'(A,F10.6)') 'residual  0: ', Compute_RMS(residual_vect(:), num_obs_points)
+    write(*,'(A,F10.6)') 'residual  0: ', sqrt(sum(residual_vect(1:num_obs_points)**2) / float(num_obs_points))
 
     do j = 1, nsf_local
         if(use_orig_data) residual_vect(:) = obs%field(1:num_obs_points)
@@ -341,8 +341,7 @@ subroutine NLSF_Least_Sq_Fit(obs, nlsf, save_data)
             field_precond(:) = NLSF_Eval_Semivariance(obs, nlsf(k))
         endif
         residual_vect = NLSF_Fit_Function(obs, nlsf(j), residual_vect, field_precond)
-
-        nlsf(j)%rms = Compute_RMS(residual_vect(:), num_obs_points)
+        nlsf(j)%rms = sqrt(sum(residual_vect(1:num_obs_points)**2) / float(num_obs_points))
         rms(j) = nlsf(j)%rms
         residuals(:, j) = residual_vect(:)
         write(*,'(A,I2,A,F10.6, A, A3, A, I2, A, I2)')'residual ', j, ': ', nlsf(j)%rms, &
@@ -375,8 +374,9 @@ subroutine NLSF_Least_Sq_Fit(obs, nlsf, save_data)
         nsf_local = min(nsf_local,nsflim_local)
 
         write(*,'(A,'//trim(buf)//'F10.6)') 'function rms:', nlsf(1:nsf_local)%rms
-        mu = Compute_MEAN(obs%field(1:num_obs_points), num_obs_points)
-        rms0 = Compute_RMS(obs%field(:) - mu, num_obs_points)
+        mu = sum(obs%field(1:num_obs_points)) / float(num_obs_points)
+        rms0   = sqrt(sum((obs%field(1:num_obs_points) - mu  )**2) / float(num_obs_points))
+
         j = 1
         if (save_data) call write_csv(num_obs_points,nsf_local,residuals,'residuals0.csv',num_obs_points, .false.)
         do while( ( nlsf(j)%rms .lt. 0.9_dp * rms0 + 0.1_dp * nlsf(1)%rms ) .and. (j+1 .lt. nsf_local) )
@@ -549,7 +549,8 @@ function NLSF_Fit_Function(obs, nlsf, y, f)
 ! spf: penalizes integral [d2 f /d x2 f(x)]^2
 ! set spf=0 for no roughness penalty
 !----------------------------------------------
-rms = Compute_RMS(y(1:num_points), num_points)
+rms = sqrt( sum(y(1:num_points)**2) / float(num_points) )
+
 nlsf%lambda = nlsf%lambda_min
 smoothness_penalty = NLSF_Smooth_Penalty(nlsf)
 spf = rms/smoothness_penalty
@@ -568,7 +569,7 @@ function NLSF_Fit_Function(obs, nlsf, y, f)
         nlsf%lambda = lambda(k)
         s(:) = NLSF_Eval_Semivariance(obs, nlsf)
         lpr = LSF_Simple_Linear_Regression(y(1:num_points), s(1:num_points) * f(1:num_points), num_points)
-        rms = Compute_RMS(y(1:num_points) - lpr(1:num_points), num_points)
+        rms = sqrt( sum((y(1:num_points) - lpr(1:num_points))**2) / float(num_points) )
         smoothness_penalty = NLSF_Smooth_Penalty(nlsf) ! penalize roughness analytic approximation
         err = rms + spf * smoothness_penalty
         if (err.lt.err_min)then

UKsrc/UniversalKriging.f90

@@ -148,9 +148,8 @@ program UKsimulation
 
 
 !SPECIAL CASE, ALL OBSERVED DATA POINTS ARE 0.0
-
-if (sum(obs%field(1:num_obs_points)) .GE. zero_thresh) then
-    SF = Compute_MEAN(obs%field(1:num_obs_points), num_obs_points) / 5.D0
+if (sum(obs%field(1:num_obs_points)) .GE. zero_threshold) then
+    SF = (sum(obs%field(1:num_obs_points)) / float(num_obs_points)) / 5.D0
     obs%field(1:num_obs_points) = log((one_scallop_per_tow + obs%field(1:num_obs_points)) / SF)
 
     !-----------------------------------------------------------------------------------------------------
@@ -182,7 +181,7 @@ program UKsimulation
     spatial_fcn = Krig_Eval_Spatial_Function(obs, num_spat_fcns, num_obs_points, nlsf, save_data)
     residual = LSF_Generalized_Least_Squares&
     &       (obs%field, spatial_fcn, residual_cov, num_obs_points, num_spat_fcns, beta, Cbeta, save_data)
-    write(*,'(A,F10.6)')'Ordinary Least Sq Residual:', Compute_RMS( residual(:), num_obs_points)
+    write(*,'(A,F10.6)')'Ordinary Least Sq Residual:', sqrt(sum(residual(1:num_obs_points)**2) / float(num_obs_points))
     if (save_data) then
         call Write_Vector_Scalar_Field(num_obs_points, residual, 'OLSresidual.txt')
         call Write_Vector_Scalar_Field(num_obs_points, obs%field, 'data.txt')
@@ -211,7 +210,7 @@ program UKsimulation
     ! COMPUTED BUT NOT USED except for GLSres
     trndOBS =  matmul(spatial_fcn, beta)
     resOBS(1:num_obs_points) = obs%field(1:num_obs_points) - trndOBS(1:num_obs_points)
-    write(*,'(A,F10.6)')'GLSres:', Compute_RMS(resOBS(:), num_obs_points)
+    write(*,'(A,F10.6)')'GLSres:', sqrt(sum(resOBS(1:num_obs_points)**2) / float(num_obs_points))
 
     deallocate(spatial_fcn)
     deallocate( beta, Cbeta)

mfiles/PlotLatLonGridSurvey.m

@@ -337,24 +337,26 @@ function SetColorbar(isOctave)
 end % function
 
 function PlotGrid(domain, thisTitle, isOctave, surveyLon, surveyLat, surveyData, gridLon, gridLat, gridData)
+gridPtSize = 3;
+survPtSize = 12;
 f = figure('Name', thisTitle);
 if isOctave
     pSurvey = scatter(surveyLon, surveyLat, surveyData, surveyData, "filled");
-    set(pSurvey, 'sizedata', 10); % size of dots
+    set(pSurvey, 'sizedata', survPtSize); % size of dots
     hold on
 
     pGrid = scatter(gridLon, gridLat, gridData, gridData, "filled");
-    set(pGrid, 'sizedata', 3); % size of dots
+    set(pGrid, 'sizedata', gridPtSize); % size of dots
 
     % enlarge figure
     figure(f,"position",get(0,"screensize"))
 else
     pSurvey = geoscatter(surveyLat, surveyLon, surveyData, surveyData, "filled");
-    pSurvey.SizeData = 10; % size of dots
+    pSurvey.SizeData = survPtSize; % size of dots
     hold on
 
     pGrid = geoscatter(gridLat, gridLon, gridData, gridData, "filled");
-    pGrid.SizeData = 3; % size of dots
+    pGrid.SizeData = gridPtSize; % size of dots
 	% enlarge figure
     if strcmp(domain, 'GB')
     	f.OuterPosition = [1963.4 -221.4 1500 1087.2];