Solving the signedness problem with integers and releasing version 2.0.0

* Fixes #268; Solving the signedness problem with integers.

* Small clean-up.

* Switch to use size_t everywhere.

* Fixing numpy seed in the example script for easier debugging.

* CHANGELOG updated for python3 release.

* Mentioning JOSS paper update and v1.2.2 mentioned in the CHANGELOG.

CHANGELOG.rst

@@ -33,13 +33,13 @@ and this project adheres to
 .. Attribution
 .. ^^^^^^^^^^^
 
-[v2.0.0] - 2023-02-03
+[v2.0.0] - 2023-02-16
 ~~~~~~~~~~~~~~~~~~~~~
 
 Summary
 ^^^^^^^
 
-* This major release candidate migrates X-PSI from Python2 (X-PSI v1.2.1 or lower) to Python3 (X-PSI v2.0 and higher), with corresponding updates and improvements to all documentation and tutorials.  
+* This major release migrates X-PSI from Python2 (X-PSI v1.2.1 or lower) to Python3 (X-PSI v2.0 and higher), with corresponding updates and improvements to all documentation and tutorials.
 
 Fixed
 ^^^^^
@@ -53,32 +53,36 @@ Added
 * Post-processing - adding names of parameters across diagonal in corner plots
 * Extra yticks options for plotting functions in the tutorials
 * `--noopenmp` install option for Mac Users
-* Added option to fix the random seed for the synthetic data generation in Python3 version
+* Added option to fix the random seed for the synthetic data generation in Python3 version.
+* Added option to plot y-axis in the residuals in a user selected scale (e.g., either log or lin).
 
 Changed
 ^^^^^^^
 
-* Modified all X-PSI routines to work in Python3
+* Modified all X-PSI routines to work in Python3.
 * General Documentation (Applications, Team and Acknowledgements, Citation, Future pages) updated for both Python2 and Python3 documentation branches.
-* Installation and tutorial pages modified for Python3
-* Module generator updated for Python3 and documentation added
-* Projection tool updated for Python3 and documentation added
-* Github actions modified to work in Python3
-* Github actions modified to use mamba with install commands on one line to improve speed
-* Updated references in the documentation and tutorial notebooks
-* CustomInstrument channel_edges argument now changed to mandatory in tutorial notebooks and examples
-* X-PSI Postprocessing now supports up-to-date versions of Nestcheck, Getdist. 
+* Installation and tutorial pages modified for Python3.
+* Module generator updated for Python3 and documentation added.
+* Projection tool updated for Python3 and documentation added.
+* Github actions modified to work in Python3.
+* Github actions modified to use mamba with install commands on one line to improve speed.
+* Updated references in the documentation and tutorial notebooks.
+* CustomInstrument channel_edges argument now changed to mandatory in tutorial notebooks and examples.
+* X-PSI Postprocessing now supports up-to-date versions of NestCheck and GetDist.
+* Specified the integer types to be always size_t in Cython files in those integer comparisons that raised warnings for different signedness of integers.
+* The JOSS paper has been updated to link to published version.
+* A final Python2 release of X-PSI (v1.2.2) was created in the Python2 branch to match the JOSS publication.
 
 Deprecated
 ^^^^^^^^^^
 
-* The Python2 version of X-PSI (v1.2.1) is now considered deprecated, although documentation and tutorials are still available.
+* The Python2 version of X-PSI (v1.2.2) is now considered deprecated, although documentation and tutorials are still available.
 
 Removed
 ^^^^^^^
 
-* Removed requirement of FFMPEG for Animations in tutorials
-* Suppressed printf() statements from c code in tutorial notebooks
+* Removed requirement of FFMPEG for Animations in tutorials.
+* Suppressed printf() statements from c code in tutorial notebooks.
 
 Attribution
 ^^^^^^^^^^^

examples/examples_fast/Modules/main.py

@@ -6,6 +6,8 @@
 import xpsi
 from xpsi.global_imports import gravradius
 
+np.random.seed(xpsi._rank+10)
+
 import time
 
 from CustomInstrument import CustomInstrument

xpsi/cellmesh/integrator.pyx

@@ -207,16 +207,16 @@ def integrate(size_t numThreads,
         cos_deflection = np.zeros((deflection.shape[0],
                                    deflection.shape[1]),
                                    dtype = np.double)
-        for i in range(deflection.shape[0]):
-            for j in range(deflection.shape[1]):
+        for i in range(<size_t>deflection.shape[0]):
+            for j in range(<size_t>deflection.shape[1]):
                 cos_deflection[i,j] = cos(deflection[i, N_R - j - 1])
 
         cos_alpha_alt = np.zeros((cos_alpha.shape[0],
                                   cos_alpha.shape[1]),
                                   dtype = np.double)
 
-        for i in range(cos_alpha.shape[0]):
-            for j in range(cos_alpha.shape[1]):
+        for i in range(<size_t>cos_alpha.shape[0]):
+            for j in range(<size_t>cos_alpha.shape[1]):
                 cos_alpha_alt[i,j] = cos_alpha[i, N_R - j - 1]
 
     if image_order_limit is not None:
@@ -274,12 +274,12 @@ def integrate(size_t numThreads,
         j = 0
         # use this to decide whether or not to compute parallel:
         # does the local vicinity of the parallel contain radiating material?
-        while j < cellArea.shape[1]:
+        while j < <size_t>cellArea.shape[1]:
             if CELL_RADIATES[i,j] == 1:
                 break
             j = j + 1
 
-        if j == cellArea.shape[1]:
+        if j == <size_t>cellArea.shape[1]:
             continue
 
         gsl_interp_accel_reset(accel_alpha[T])
@@ -561,7 +561,7 @@ def integrate(size_t numThreads,
                     gsl_interp_init(interp_GEOM[T], phase_ptr, GEOM_ptr, N_L)
 
                     j = 0
-                    while j < cellArea.shape[1] and terminate[T] == 0:
+                    while j < <size_t>cellArea.shape[1] and terminate[T] == 0:
                         if CELL_RADIATES[i,j] == 1:
                             phi_shift = phi[i,j]
                             for k in range(N_P):

xpsi/cellmesh/integrator_for_azimuthal_invariance.pyx

@@ -212,16 +212,16 @@ def integrate(size_t numThreads,
         cos_deflection = np.zeros((deflection.shape[0],
                                    deflection.shape[1]),
                                    dtype = np.double)
-        for i in range(deflection.shape[0]):
-            for j in range(deflection.shape[1]):
+        for i in range(<size_t>deflection.shape[0]):
+            for j in range(<size_t>deflection.shape[1]):
                 cos_deflection[i,j] = cos(deflection[i, N_R - j - 1])
 
         cos_alpha_alt = np.zeros((cos_alpha.shape[0],
                                   cos_alpha.shape[1]),
                                   dtype = np.double)
 
-        for i in range(cos_alpha.shape[0]):
-            for j in range(cos_alpha.shape[1]):
+        for i in range(<size_t>cos_alpha.shape[0]):
+            for j in range(<size_t>cos_alpha.shape[1]):
                 cos_alpha_alt[i,j] = cos_alpha[i, N_R - j - 1]
 
     if image_order_limit is not None:
@@ -279,13 +279,13 @@ def integrate(size_t numThreads,
         j = 0
         # use this to decide whether or not to compute parallel:
         # Does the local vicinity of the parallel contain radiating material?
-        while j < cellArea.shape[1]:
+        while j < <size_t>cellArea.shape[1]:
             if CELL_RADIATES[i,j] == 1:
                 J = j
                 break
             j = j + 1
 
-        if j == cellArea.shape[1]:
+        if j == <size_t>cellArea.shape[1]:
             continue
 
         gsl_interp_accel_reset(accel_alpha[T])
@@ -557,7 +557,7 @@ def integrate(size_t numThreads,
                         gsl_interp_init(interp_PROFILE[T], phase_ptr, profile_ptr, N_L)
 
                         j = 0
-                        while j < cellArea.shape[1] and terminate[T] == 0:
+                        while j < <size_t>cellArea.shape[1] and terminate[T] == 0:
                             if CELL_RADIATES[i,j] == 1:
                                 phi_shift = phi[i,j]
                                 for k in range(N_P):

xpsi/cellmesh/integrator_for_time_invariance.pyx

@@ -164,16 +164,16 @@ def integrate(size_t numThreads,
         cos_deflection = np.zeros((deflection.shape[0],
                                    deflection.shape[1]),
                                    dtype = np.double)
-        for i in range(deflection.shape[0]):
-            for j in range(deflection.shape[1]):
+        for i in range(<size_t>deflection.shape[0]):
+            for j in range(<size_t>deflection.shape[1]):
                 cos_deflection[i,j] = cos(deflection[i, N_R - j - 1])
 
         cos_alpha_alt = np.zeros((cos_alpha.shape[0],
                                   cos_alpha.shape[1]),
                                   dtype = np.double)
 
-        for i in range(cos_alpha.shape[0]):
-            for j in range(cos_alpha.shape[1]):
+        for i in range(<size_t>cos_alpha.shape[0]):
+            for j in range(<size_t>cos_alpha.shape[1]):
                 cos_alpha_alt[i,j] = cos_alpha[i, N_R - j - 1]
 
     if image_order_limit is not None:

xpsi/likelihoods/_poisson_likelihood_given_background.pyx

@@ -104,7 +104,7 @@ def poisson_likelihood_given_background(double exposure_time,
             raise TypeError('An iterable is required to specify component-by-'
                             'component positivity.')
         else:
-            if len(allow_negative) != num_components:
+            if <size_t>len(allow_negative) != num_components:
                 raise ValueError('Number of allow_negative declarations does '
                                  'not match the number of components..')
 
@@ -130,7 +130,7 @@ def poisson_likelihood_given_background(double exposure_time,
     cdef gsl_interp *inter_ptr = NULL
     cdef accel *acc_ptr = NULL
 
-    for i in range(STAR.shape[0]):
+    for i in range(<size_t>STAR.shape[0]):
         for p in range(num_components):
             signal = components[p]
             signal_phase_set = component_phases[p]
@@ -144,7 +144,7 @@ def poisson_likelihood_given_background(double exposure_time,
             gsl_interp_init(interp_ptr, phases_ptr, signal_ptr,
                             signal_phase_set.shape[0])
 
-            for j in range(phases.shape[0] - 1):
+            for j in range(<size_t>phases.shape[0] - 1):
                 a = phases[j] + phase_shift
                 b = phases[j+1] + phase_shift
 
@@ -180,7 +180,7 @@ def poisson_likelihood_given_background(double exposure_time,
                     if _val > 0.0 or _allow_negative[p] == 1:
                         STAR[i,j] += _val
 
-        for j in range(phases.shape[0] - 1): # interpolant safety procedure
+        for j in range(<size_t>phases.shape[0] - 1): # interpolant safety procedure
             if STAR[i,j] < 0.0:
                 STAR[i,j] = 0.0
 
@@ -195,8 +195,8 @@ def poisson_likelihood_given_background(double exposure_time,
         double LOGLIKE = 0.0, EXPEC = 0.0
         double n = <double>(phases.shape[0] - 1)
 
-    for i in range(STAR.shape[0]):
-        for j in range(STAR.shape[1]):
+    for i in range(<size_t>STAR.shape[0]):
+        for j in range(<size_t>STAR.shape[1]):
             EXPEC = (STAR[i,j] + background[i,j]/n) * exposure_time
             LOGLIKE -= EXPEC
             LOGLIKE += counts[i,j] * log(EXPEC)

xpsi/likelihoods/default_background_marginalisation.pyx

@@ -56,8 +56,8 @@ def precomputation(int[:,::1] data):
         size_t i, j
         double[::1] precomp = np.zeros(data.shape[0], dtype = np.double)
 
-    for i in range(data.shape[0]):
-        for j in range(data.shape[1]):
+    for i in range(<size_t>data.shape[0]):
+        for j in range(<size_t>data.shape[1]):
             precomp[i] += gsl_sf_lnfact(<unsigned int>(data[i,j]))
 
         precomp[i] *= -1.0
@@ -79,7 +79,7 @@ ctypedef struct args:
 
 cdef double marginal_integrand(double B, void *params) nogil:
 
-    cdef int j
+    cdef size_t j
     cdef double c, x = 0.0
     cdef args *a  = <args*> params
 
@@ -99,7 +99,7 @@ cdef double marginal_integrand(double B, void *params) nogil:
 
 cdef double delta(double B, void *params) nogil:
 
-    cdef int j
+    cdef size_t j
     cdef double x = 0.0, y = 0.0
     cdef args *a  = <args*> params
 
@@ -304,7 +304,7 @@ def eval_marginal_likelihood(double exposure_time,
             raise TypeError('An iterable is required to specify component-by-'
                             'component positivity.')
         else:
-            if len(allow_negative) != num_components:
+            if <size_t> len(allow_negative) != num_components:
                 raise ValueError('Number of allow_negative declarations does '
                                  'not match the number of components..')
 
@@ -325,7 +325,7 @@ def eval_marginal_likelihood(double exposure_time,
     cdef gsl_interp *inter_ptr = NULL
     cdef accel *acc_ptr = NULL
 
-    for i in range(STAR.shape[0]):
+    for i in range(<size_t> STAR.shape[0]):
         for p in range(num_components):
             pulse = components[p]
             pulse_phase_set = component_phases[p]
@@ -339,7 +339,7 @@ def eval_marginal_likelihood(double exposure_time,
             gsl_interp_init(interp_ptr, phases_ptr, pulse_ptr,
                             pulse_phase_set.shape[0])
 
-            for j in range(phases.shape[0] - 1):
+            for j in range(<size_t> (phases.shape[0] - 1)):
                 pa = phases[j] + phase_shift
                 pb = phases[j+1] + phase_shift
 
@@ -375,13 +375,13 @@ def eval_marginal_likelihood(double exposure_time,
                     if _val > 0.0 or _allow_negative[p] == 1:
                         STAR[i,j] += _val
 
-        for j in range(phases.shape[0] - 1): # interpolant safety procedure
+        for j in range(<size_t> (phases.shape[0] - 1)): # interpolant safety procedure
             if STAR[i,j] < 0.0:
                 STAR[i,j] = 0.0
 
         av_DATA = 0.0; av_STAR = 0.0
 
-        for j in range(phases.shape[0] - 1):
+        for j in range(<size_t> (phases.shape[0] - 1)):
             STAR[i,j] *= n
             if background is not None:
                 STAR[i,j] += _background[i,j]

xpsi/pixelmesh/RK_IP2S_tracer.pyx

@@ -160,7 +160,7 @@ cdef int RK(_RAY *const RAY,
         if status == GSL_SUCCESS:
             RAY.NUMSTEPS += 1
 
-            if RAY.NUMSTEPS >= RAY.MAXSTEPS:
+            if <size_t> RAY.NUMSTEPS >= RAY.MAXSTEPS:
                 printf("\n\nSteps exceeded... %d", RAY.NUMSTEPS)
                 printf("\nx: %.8e; y: %.8e", X_IP, Y_IP)
                 printf("\ny[]: %.8e, %.8e, %.8e, %.8e, %.8e, %.8e",
@@ -296,7 +296,7 @@ cdef int RK(_RAY *const RAY,
                 RAY.NUM_SINGULARITY_STEPS += 1
                 #printf("\nSingularity steps: %d", RAY.NUM_SINGULARITY_STEPS)
 
-                if RAY.NUM_SINGULARITY_STEPS == RAY.MAXSTEPS:
+                if <size_t>RAY.NUM_SINGULARITY_STEPS == RAY.MAXSTEPS:
                     RAY.EVOLVE = 0
                     break
 

xpsi/pixelmesh/integrator.pyx

@@ -365,8 +365,8 @@ def integrate(size_t numThreads,
 
                             if cache_intensities == 1:
                                 for search in range(N_CP):
-                                    if k == cache_phase_indices[search]:
-                                        if cache_counter_E < N_CE and p == cache_energy_indices[cache_counter_E]:
+                                    if k == <size_t> cache_phase_indices[search]:
+                                        if cache_counter_E < N_CE and p == <size_t> cache_energy_indices[cache_counter_E]:
                                             if single_precision_cache == 0:
                                                 IMAGE_DUB[cache_phase_indices[search], cache_counter_E, INDEX + 1] = I_E / energies[p]
                                             else:
@@ -397,17 +397,17 @@ def integrate(size_t numThreads,
 
                     if cache_intensities == 1:
                         for search in range(N_CP):
-                            if k == cache_phase_indices[search]:
-                                if cache_counter_E < N_CE and p == cache_energy_indices[cache_counter_E]:
+                            if k == <size_t> cache_phase_indices[search]:
+                                if cache_counter_E < N_CE and p == <size_t> cache_energy_indices[cache_counter_E]:
                                     if single_precision_cache == 0:
                                         IMAGE_DUB[cache_phase_indices[search], cache_counter_E, 0] = I_E / energies[p]
                                     else:
                                         IMAGE[cache_phase_indices[search], cache_counter_E, 0] = I_E / energies[p]
 
                 if cache_intensities == 1:
                     for search in range(N_CP):
-                        if k == cache_phase_indices[search]:
-                            if cache_counter_E < N_CE and p == cache_energy_indices[cache_counter_E]:
+                        if k == <size_t> cache_phase_indices[search]:
+                            if cache_counter_E < N_CE and p == <size_t> cache_energy_indices[cache_counter_E]:
                                 cache_counter_E = cache_counter_E + 1
 
                 integrated_flux[k,p] *= SEMI_MINOR * SEMI_MAJOR

xpsi/surface_radiation_field/__init__.pyx

@@ -452,9 +452,9 @@ def effective_gravity(double[::1] cos_colatitude,
     cdef double[::1] gravity = np.zeros(cos_colatitude.shape[0],
                                         dtype=np.double)
 
-    cdef size_t i
+    cdef unsigned int i
 
-    for i in range(gravity.shape[0]):
+    for i in range(<size_t>gravity.shape[0]):
         gravity[i] = effectiveGravity(cos_colatitude[i],
                                       R_eq[i],
                                       zeta[i],

xpsi/surface_radiation_field/archive/hot/numerical_fbeam.pyx

@@ -399,7 +399,7 @@ cdef double eval_hot(size_t THREAD,
         I_hot = anorm*(1.0+abb*((E)**cbb)*mu+bbb*((E)**dbb)*mu**2)*I_nsx   
     if emodel==3:          
         #integral using trapezoidal rule (consider changing to Gauss quadrature)          
-        for imu in range(int(nimu)):
+        for imu in range(<size_t> nimu):
             mu_imu = mu_imu+(1.0/nimu)
             if imu==0 or imu==nimu-1:
                     dmu = (0.5/nimu)

xpsi/tools/energy_integrator.pyx

@@ -82,13 +82,13 @@ def energy_integrator(size_t N_Ts,
         T = threadid()
         cpy = _signal[T]
 
-        for j in range(energies.shape[0]):
+        for j in range(<size_t>energies.shape[0]):
             cpy[j] = pow(10.0, energies[j]) * signal[j,i] * log(10.0)
 
         gsl_interp_accel_reset(acc[T])
         gsl_interp_init(interp[T], &(energies[0]), cpy, energies.shape[0])
 
-        for j in range(energy_edges.shape[0] - 1):
+        for j in range(<size_t>energy_edges.shape[0] - 1):
             if energy_edges[j + 1] > max_energy:
                 upper_energy = max_energy
             else: