From 6ae7566c9b92032fc20a49f2d2f5c75e6e088289 Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Thu, 8 Apr 2021 08:54:47 +0200
Subject: [PATCH 01/30] requirements: astropy version 4.1 instead of newest one
 (erfa bug)

---
 requirements.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/requirements.txt b/requirements.txt
index 7e45fe779..eb0c804b4 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,4 @@
-astropy>=3.1.1
+astropy=4.1
 joblib>=0.13.1
 matplotlib>=3.0.2
 modopt>=1.2.0

From 8c8b8d3861f8592b8a299928c8a1259625a73f6e Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Thu, 8 Apr 2021 09:34:47 +0200
Subject: [PATCH 02/30] requirements astropy version fix

---
 requirements.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/requirements.txt b/requirements.txt
index eb0c804b4..bf39cccbe 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,4 @@
-astropy=4.1
+astropy==4.1
 joblib>=0.13.1
 matplotlib>=3.0.2
 modopt>=1.2.0

From b96b8e0285f4782baa41b0d305d748a5c3968d97 Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Tue, 13 Apr 2021 12:48:18 +0000
Subject: [PATCH 03/30] small changes in job subm scripts

---
 scripts/python/canfar_run_analyse.py  | 6 +++---
 scripts/sh/canfar_submit_selection.sh | 2 +-
 scripts/sh/job_sp.bash                | 2 +-
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/scripts/python/canfar_run_analyse.py b/scripts/python/canfar_run_analyse.py
index 5cb613735..820cc7be7 100755
--- a/scripts/python/canfar_run_analyse.py
+++ b/scripts/python/canfar_run_analyse.py
@@ -171,13 +171,13 @@ def get_status(tile_num):
             final_cat_found = False
             with open(out_name) as out_file:
                 for line in out_file:
-                    m = re.match('Upload final_cat to.*(\d+) files', line)
+                    m = re.match('Upload.*final_cat.*return value = (\.*)', line)
                     if m:
                         final_cat_found = True
                         if m[1] == '0':
-                            status = res_noout, '0 final_cat output files'
+                            status = res_ok, 'successful upload of final_cat'
                         else:
-                            status = res_ok, 'success, {} final_cat output files'.format(m[1])
+                            status = res_noout, 'failed upload of final_cat'
                         break
 
             if final_cat_found == False:
diff --git a/scripts/sh/canfar_submit_selection.sh b/scripts/sh/canfar_submit_selection.sh
index 2e7fa280f..a2fe5b7d3 100644
--- a/scripts/sh/canfar_submit_selection.sh
+++ b/scripts/sh/canfar_submit_selection.sh
@@ -11,7 +11,7 @@
 
 # Job file for one tile
 SP_ROOT=$HOME/shapepipe
-sp_job="$SP_ROOT/scripts/sh/canfar_sp.bash"
+sp_job="$SP_ROOT/scripts/sh/job_sp.bash"
 
 
 # Functions
diff --git a/scripts/sh/job_sp.bash b/scripts/sh/job_sp.bash
index 2740a7363..be17542fc 100755
--- a/scripts/sh/job_sp.bash
+++ b/scripts/sh/job_sp.bash
@@ -16,7 +16,7 @@
 ## Default values
 do_env=0
 job=255
-psf='mccd'
+psf='psfex'
 retrieve='vos'
 nsh_step=3500
 nsh_max=-1

From 6fa5d2fca32954dc85d1073bfc4fb21f440d9ca4 Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Wed, 21 Apr 2021 18:17:55 +0200
Subject: [PATCH 04/30] libs: astropy 4.1, no mccd, installation works

---
 environment.yml  | 5 ++---
 requirements.txt | 3 +--
 2 files changed, 3 insertions(+), 5 deletions(-)

diff --git a/environment.yml b/environment.yml
index a8749e078..081095abc 100644
--- a/environment.yml
+++ b/environment.yml
@@ -1,10 +1,10 @@
-name: shapepipe
+name: sptest
 channels:
   - conda-forge
   - defaults
 dependencies:
   - python=3.7
-  - astropy>=3.0
+  - astropy=4.1
   - galsim>=2.1.4
   - joblib>=0.13
   - matplotlib>=3.0.2
@@ -24,6 +24,5 @@ dependencies:
     - PyQt5>=5.15.1
     - pyqtgraph>=0.11.0
     - python-pysap>=0.0.4
-    - mccd>=0.0.3
     - treecorr>=4.1.5
     - git+https://github.com/tobias-liaudat/Stile
diff --git a/requirements.txt b/requirements.txt
index 7e45fe779..843774856 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,4 @@
-astropy>=3.1.1
+astropy==4.1
 joblib>=0.13.1
 matplotlib>=3.0.2
 modopt>=1.2.0
@@ -14,5 +14,4 @@ galsim>=2.1.4
 PyQt5>=5.15.1
 pyqtgraph>=0.11.0
 python-pysap>=0.0.4
-mccd>=0.0.3
 treecorr>=4.1.5

From c2f395dcc6222f649678a0556522cf57697fd414 Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Mon, 26 Apr 2021 13:56:59 +0200
Subject: [PATCH 05/30] job submission script: replace 'find' to avoid stderr

---
 scripts/sh/job_sp.bash | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/scripts/sh/job_sp.bash b/scripts/sh/job_sp.bash
index 2740a7363..b77419f88 100755
--- a/scripts/sh/job_sp.bash
+++ b/scripts/sh/job_sp.bash
@@ -361,17 +361,20 @@ if [[ $do_job != 0 ]]; then
 
   ### The following are very a bad hacks to get additional input file paths
   if [ "$psf" == "psfex" ]; then
-    input_psfex=`find . -name star_split_ratio_80-*.psf | head -n 1`
+    #input_psfex=`find . -name star_split_ratio_80-*.psf | head -n 1`
+    input_psfex=``ls -1 ./output/*/psfex_runner/output/star_split_ratio_80*.psf | head -n 1`
     command_sp "ln -s `dirname $input_psfex` input_psfex" "Link psfex output"
   else
     input_psf_mccd=`find . -name "fitted_model*.npy" | head -n 1`
     command_sp "ln -s `dirname $input_psf_mccd` input_psf_mccd" "Link MCCD output"
   fi
 
-  input_split_exp=`find output -name flag-*.fits | head -n 1`
+  #input_split_exp=`find output -name flag-*.fits | head -n 1`
+  input_split_exp=`ls -1 ./output/*/split_exp_runner/output/flag*.fits | head -n 1`
   command_sp "ln -s `dirname $input_split_exp` input_split_exp" "Link split_exp output"
 
-  input_sextractor=`find . -name sexcat_sexcat-*.fits | head -n 1`
+  #input_sextractor=`find . -name sexcat_sexcat-*.fits | head -n 1`
+  input_sextractor=`ls -1 ./output/*/sextractor_runner/output/sexcat_sexcat*.fits | head -n 1`
   command_sp "ln -s `dirname $input_sextractor` input_sextractor" "Link sextractor output"
 
 fi

From 63fca070d3849de4fcc830feec76a04aadd5eecc Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Sun, 2 May 2021 09:17:52 +0200
Subject: [PATCH 06/30] restored main version of req/env

---
 environment.yml  | 5 +++--
 requirements.txt | 3 ++-
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/environment.yml b/environment.yml
index 081095abc..a8749e078 100644
--- a/environment.yml
+++ b/environment.yml
@@ -1,10 +1,10 @@
-name: sptest
+name: shapepipe
 channels:
   - conda-forge
   - defaults
 dependencies:
   - python=3.7
-  - astropy=4.1
+  - astropy>=3.0
   - galsim>=2.1.4
   - joblib>=0.13
   - matplotlib>=3.0.2
@@ -24,5 +24,6 @@ dependencies:
     - PyQt5>=5.15.1
     - pyqtgraph>=0.11.0
     - python-pysap>=0.0.4
+    - mccd>=0.0.3
     - treecorr>=4.1.5
     - git+https://github.com/tobias-liaudat/Stile
diff --git a/requirements.txt b/requirements.txt
index 843774856..7e45fe779 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,4 @@
-astropy==4.1
+astropy>=3.1.1
 joblib>=0.13.1
 matplotlib>=3.0.2
 modopt>=1.2.0
@@ -14,4 +14,5 @@ galsim>=2.1.4
 PyQt5>=5.15.1
 pyqtgraph>=0.11.0
 python-pysap>=0.0.4
+mccd>=0.0.3
 treecorr>=4.1.5

From 417c4fccd0ada7e68db9f74a1523aef3a24c9219 Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Sun, 2 May 2021 09:32:50 +0200
Subject: [PATCH 07/30] job_sp: find -> ls

---
 scripts/sh/job_sp.bash | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/scripts/sh/job_sp.bash b/scripts/sh/job_sp.bash
index b77419f88..82cd40b3e 100755
--- a/scripts/sh/job_sp.bash
+++ b/scripts/sh/job_sp.bash
@@ -362,7 +362,7 @@ if [[ $do_job != 0 ]]; then
   ### The following are very a bad hacks to get additional input file paths
   if [ "$psf" == "psfex" ]; then
     #input_psfex=`find . -name star_split_ratio_80-*.psf | head -n 1`
-    input_psfex=``ls -1 ./output/*/psfex_runner/output/star_split_ratio_80*.psf | head -n 1`
+    input_psfex=`ls -1 ./output/*/psfex_runner/output/star_split_ratio_80*.psf | head -n 1`
     command_sp "ln -s `dirname $input_psfex` input_psfex" "Link psfex output"
   else
     input_psf_mccd=`find . -name "fitted_model*.npy" | head -n 1`

From 5a7da0b401069f1c3be887035fbb217fc8685286 Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Fri, 7 May 2021 09:20:50 +0200
Subject: [PATCH 08/30] post_proc mccd to improve

---
 scripts/python/merge_final_cat.py | 2 +-
 scripts/sh/psf_residuals.bash     | 1 +
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/scripts/python/merge_final_cat.py b/scripts/python/merge_final_cat.py
index 42ba95f99..be89824bb 100755
--- a/scripts/python/merge_final_cat.py
+++ b/scripts/python/merge_final_cat.py
@@ -305,7 +305,7 @@ def main(argv=None):
     print('Saving final catalogue')
     np.save('final_cat.npy', d)
 
-    mgs = '{} catalog files merged with success'.format(count)
+    msg = '{} catalog files merged with success'.format(count)
     if param.verbose:
         print(msg)
     print(msg, file=f_log)
diff --git a/scripts/sh/psf_residuals.bash b/scripts/sh/psf_residuals.bash
index 4929350b6..1c6d5d488 100755
--- a/scripts/sh/psf_residuals.bash
+++ b/scripts/sh/psf_residuals.bash
@@ -114,6 +114,7 @@ fi
 n_skipped=0
 n_created=0
 FILES=output/*/${runner}/output/${psfval_file_base}*
+echo ${FILES[@]}
 for val in ${FILES[@]}; do
     base=`basename $val`
     link_s "$pwd/$val" "$dir_individual/$base"

From 8bd9b11bfb9cd58db5c8fffe5fd4da36dd6cc1ff Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Fri, 7 May 2021 11:04:03 +0200
Subject: [PATCH 09/30] merge MCCD star cat script: added required X, Y

---
 scripts/python/merge_star_cat_mccd.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/scripts/python/merge_star_cat_mccd.py b/scripts/python/merge_star_cat_mccd.py
index 936526a74..9957d35f3 100755
--- a/scripts/python/merge_star_cat_mccd.py
+++ b/scripts/python/merge_star_cat_mccd.py
@@ -175,6 +175,7 @@ def create_full_cat(input_dir, input_file_pattern, output_path, verbose=False):
         print('Found {} files'.format(len(starcat_names)))
 
     # Initialise columns to be saved
+    x, y = [], []
     ra, dec = [], []
     g1_psf, g2_psf, size_psf = [], [], []
     g1, g2, size = [], [], []
@@ -187,6 +188,9 @@ def create_full_cat(input_dir, input_file_pattern, output_path, verbose=False):
         starcat_j = fits.open(name)
 
         # positions
+        x += list(starcat_j[1].data['GLOB_POSITION_IMG_LIST'][:, 0])
+        y += list(starcat_j[1].data['GLOB_POSITION_IMG_LIST'][:, 1])
+
         ra += list(starcat_j[1].data['RA_LIST'])
         dec += list(starcat_j[1].data['DEC_LIST'])
 
@@ -211,7 +215,7 @@ def create_full_cat(input_dir, input_file_pattern, output_path, verbose=False):
 
     # Collect columns
     # convert back to sigma for consistency
-    data = {'RA': ra, 'DEC': dec,
+    data = {'X': x, 'Y': y, 'RA': ra, 'DEC': dec,
             'E1_PSF_HSM': g1_psf, 'E2_PSF_HSM': g2_psf, 'SIGMA_PSF_HSM': np.sqrt(size_psf),
             'E1_STAR_HSM': g1, 'E2_STAR_HSM': g2, 'SIGMA_STAR_HSM': np.sqrt(size),
             'FLAG_PSF_HSM': flag_psf, 'FLAG_STAR_HSM': flag_star,

From e9caaebcc1b790141c99205a5cc0ef4e6dc2df37 Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Fri, 4 Jun 2021 08:20:50 +0200
Subject: [PATCH 10/30] cfis -> utilities

---
 scripts/python/create_sample_results.py         | 2 +-
 shapepipe/utilities/__init__.py                 | 2 +-
 {scripts/python => shapepipe/utilities}/cfis.py | 0
 3 files changed, 2 insertions(+), 2 deletions(-)
 rename {scripts/python => shapepipe/utilities}/cfis.py (100%)

diff --git a/scripts/python/create_sample_results.py b/scripts/python/create_sample_results.py
index 6c03145e7..640468237 100755
--- a/scripts/python/create_sample_results.py
+++ b/scripts/python/create_sample_results.py
@@ -244,7 +244,7 @@ def main(argv=None):
 
     input_IDs = read_ID_list(param.input_IDs, verbose=param.verbose)
 
-    result_base_names = ['psfex', 'psfex_interp_exp', 'setools_mask', 'setools_stat', 'setools_plot',
+    result_base_names = ['psfex_interp_exp', 'setools_mask', 'setools_stat', 'setools_plot',
                          'pipeline_flag', 'final_cat', 'logs']
 
     if os.path.isdir(param.output_dir):
diff --git a/shapepipe/utilities/__init__.py b/shapepipe/utilities/__init__.py
index 7e488988a..828282ecd 100644
--- a/shapepipe/utilities/__init__.py
+++ b/shapepipe/utilities/__init__.py
@@ -9,4 +9,4 @@
 
 """
 
-__all__ = ['canfar', 'file_system']
+__all__ = ['canfar', 'file_system', 'cfis']
diff --git a/scripts/python/cfis.py b/shapepipe/utilities/cfis.py
similarity index 100%
rename from scripts/python/cfis.py
rename to shapepipe/utilities/cfis.py

From 9dc57d3d133704b29c2c289d215501437df310a8 Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Fri, 4 Jun 2021 09:30:32 +0200
Subject: [PATCH 11/30] scripts updated vis-a-vis cfis import

---
 scripts/python/create_sample_results.py |  2 +-
 scripts/python/merge_final_cat.py       | 23 +++++++----------------
 scripts/python/merge_star_cat_mccd.py   | 20 ++++++++++----------
 scripts/python/merge_star_cat_psfex.py  |  2 +-
 4 files changed, 19 insertions(+), 28 deletions(-)

diff --git a/scripts/python/create_sample_results.py b/scripts/python/create_sample_results.py
index 640468237..122937be3 100755
--- a/scripts/python/create_sample_results.py
+++ b/scripts/python/create_sample_results.py
@@ -17,7 +17,7 @@
 import io
 from contextlib import redirect_stdout
 from optparse import OptionParser
-import cfis
+from shapepipe.utilities import cfis
 from shapepipe.utilities.file_system import mkdir
 
 
diff --git a/scripts/python/merge_final_cat.py b/scripts/python/merge_final_cat.py
index be89824bb..d0929d4f5 100755
--- a/scripts/python/merge_final_cat.py
+++ b/scripts/python/merge_final_cat.py
@@ -25,7 +25,7 @@
 
 from tqdm import tqdm
 
-import cfis
+from shapepipe.utilities import cfis
 
 
 class param:
@@ -270,16 +270,8 @@ def main(argv=None):
     for key in d_tmp.dtype.names:
         d[key] = d_tmp[key]
     count = count + 1
-    print('File \'{}\' copied ({}/{})'.format(lpath[0], count, len(lpath)))
-
-    #new_dt = np.dtype(d_tmp.dtype.descr + [('TILE_ID', '>i4')])
-    #d = np.zeros(d_tmp.shape, dtype=new_dt)
-
-    #if 'TILE_ID' in d_tmp.dtype.names:
-        #d['TILE_ID'].fill(int(''.join(re.findall('\d+', l[0]))))
-
-    # Read all final catalogues and merge
-    #for i in tqdm(lpath[1:], total=len(lpath)-1):
+    if param.verbose:
+        print('File \'{}\' copied ({}/{})'.format(lpath[0], count, len(lpath)))
 
     for i in lpath[1:]:
         if ('final_cat' not in i) | ('.npy' in i):
@@ -292,17 +284,16 @@ def main(argv=None):
             for key in d_tmp.dtype.names:
                 dd[key] = d_tmp[key]
             count = count + 1
-            print('File \'{}\' copied ({}/{})'.format(i, count, len(lpath)))
-
-            #if 'TILE_ID' in d_tmp.dtype.names:
-                #dd['TILE_ID'].fill(int(''.join(re.findall('\d+', i))))
+            if param.verbose:
+                print('File \'{}\' copied ({}/{})'.format(i, count, len(lpath)))
 
             d = np.concatenate((d, dd))
         except:
             print('Error while copying file \'{}\''.format(i))
 
     # Save merged catalogue as numpy binary file
-    print('Saving final catalogue')
+    if param.verbose:
+        print('Saving final catalogue')
     np.save('final_cat.npy', d)
 
     msg = '{} catalog files merged with success'.format(count)
diff --git a/scripts/python/merge_star_cat_mccd.py b/scripts/python/merge_star_cat_mccd.py
index 9957d35f3..03179cfb9 100755
--- a/scripts/python/merge_star_cat_mccd.py
+++ b/scripts/python/merge_star_cat_mccd.py
@@ -191,20 +191,20 @@ def create_full_cat(input_dir, input_file_pattern, output_path, verbose=False):
         x += list(starcat_j[1].data['GLOB_POSITION_IMG_LIST'][:, 0])
         y += list(starcat_j[1].data['GLOB_POSITION_IMG_LIST'][:, 1])
 
-        ra += list(starcat_j[1].data['RA_LIST'])
-        dec += list(starcat_j[1].data['DEC_LIST'])
+        ra += list(starcat_j[1].data['RA_LIST'][:])
+        dec += list(starcat_j[1].data['DEC_LIST'][:])
 
         # shapes (convert sigmas to R^2)
-        g1_psf += list(starcat_j[1].data['PSF_MOM_LIST'][0])
-        g2_psf += list(starcat_j[1].data['PSF_MOM_LIST'][1])
-        size_psf += list(starcat_j[1].data['PSF_MOM_LIST'][2]**2)
-        g1 += list(starcat_j[1].data['STAR_MOM_LIST'][0])
-        g2 += list(starcat_j[1].data['STAR_MOM_LIST'][1])
-        size += list(starcat_j[1].data['STAR_MOM_LIST'][2]**2)
+        g1_psf += list(starcat_j[1].data['PSF_MOM_LIST'][:,0])
+        g2_psf += list(starcat_j[1].data['PSF_MOM_LIST'][:,1])
+        size_psf += list(starcat_j[1].data['PSF_MOM_LIST'][:,2]**2)
+        g1 += list(starcat_j[1].data['STAR_MOM_LIST'][:,0])
+        g2 += list(starcat_j[1].data['STAR_MOM_LIST'][:,1])
+        size += list(starcat_j[1].data['STAR_MOM_LIST'][:,2]**2)
 
         # flags
-        flag_psf += list(starcat_j[1].data['PSF_MOM_LIST'][3])
-        flag_star += list(starcat_j[1].data['STAR_MOM_LIST'][3])
+        flag_psf += list(starcat_j[1].data['PSF_MOM_LIST'][:,3])
+        flag_star += list(starcat_j[1].data['STAR_MOM_LIST'][:,3])
 
         # CCD number
         ccd_nb = list(starcat_j[1].data['CCD_ID_LIST'])
diff --git a/scripts/python/merge_star_cat_psfex.py b/scripts/python/merge_star_cat_psfex.py
index c4d636594..6408ee82d 100755
--- a/scripts/python/merge_star_cat_psfex.py
+++ b/scripts/python/merge_star_cat_psfex.py
@@ -24,7 +24,7 @@
 from optparse import OptionParser
 from astropy.io import fits
 
-import cfis
+from shapepipe.utilities import cfis
 from shapepipe.pipeline import file_io as sc
 
 

From 6bd9c66020c6706013ed5bf617cefda422865db1 Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Thu, 10 Jun 2021 21:34:55 +0200
Subject: [PATCH 12/30] added utilities function for 'vls'

---
 scripts/jupyter/plot_spectro_areas.ipynb |  2 +-
 scripts/python/canfar_avail_results.py   | 26 +-----------
 shapepipe/utilities/canfar.py            | 53 ++++++++++++++++++++++--
 3 files changed, 53 insertions(+), 28 deletions(-)

diff --git a/scripts/jupyter/plot_spectro_areas.ipynb b/scripts/jupyter/plot_spectro_areas.ipynb
index 43e14f6c7..16cf32cdb 100644
--- a/scripts/jupyter/plot_spectro_areas.ipynb
+++ b/scripts/jupyter/plot_spectro_areas.ipynb
@@ -1352,7 +1352,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.9"
+   "version": "3.8.8"
   }
  },
  "nbformat": 4,
diff --git a/scripts/python/canfar_avail_results.py b/scripts/python/canfar_avail_results.py
index 17cece671..27b940b64 100755
--- a/scripts/python/canfar_avail_results.py
+++ b/scripts/python/canfar_avail_results.py
@@ -19,15 +19,11 @@
 
 from optparse import OptionParser
 
-from shapepipe.utilities.canfar import vosHandler
+from shapepipe.utilities.canfar import dir_list
 
 import cfis
 
 
-# Global VLS definition
-vls = vosHandler('vls')
-
-
 def params_default():
     """Set default parameter values.
 
@@ -205,25 +201,7 @@ def check_results(ID_files, input_vos, result_base_names, n_complete, verbose=Fa
     cmd = 'vls'
     vos_dir = 'vos:cfis/{}'.format(input_vos)
 
-    sys.argv = []
-    sys.argv.append(cmd)
-    sys.argv.append(vos_dir)
-    #sys.argv.append('> .vls.tmp')
-    #if certfile:
-        #sys.argv.append('--certfile={}'.format(certfile))
-
-    if verbose:
-        print('Getting vos directory content from vls...')
-    f = io.StringIO()
-
-    try:
-        with redirect_stdout(f):
-            vls()
-    except:
-        print('Error during vls command')
-        raise
- 
-    vls_out = f.getvalue()
+    vls_out = dir_list(vos_dir)
 
     n_found = {}
     n_IDs = {}
diff --git a/shapepipe/utilities/canfar.py b/shapepipe/utilities/canfar.py
index e3dd92cc0..3588a7c49 100644
--- a/shapepipe/utilities/canfar.py
+++ b/shapepipe/utilities/canfar.py
@@ -10,6 +10,8 @@
 
 import os
 import sys  
+from io import StringIO
+from contextlib import redirect_stdout
 
 try:
     import vos.commands as vosc
@@ -110,14 +112,59 @@ def download(source, target, verbose=False):
         status, True/False or success/failure
     """
 
+    cmd = 'vcp'
+
     if not os.path.exists(target):
-        sys.argv = ['vcp', source, target]
+        sys.argv = [cmd, source, target]
         if verbose:
             print('Downloading file {} to {}...'.format(source, target))
-        vcp = vosHandler('vcp')
+        vcp = vosHandler(cmd)
         vcp()
         if verbose:
             print('Download finished.')
     else:
         if verbose:
-            print('Target file {} exists, skipping download.'.format(target))
\ No newline at end of file
+            print('Target file {} exists, skipping download.'.format(target))
+
+
+def dir_list(path, verbose=False):
+    """list
+
+    List content of path on vos
+
+    Parameters
+    ----------
+    path : string
+        path on vos, starts with 'vos:cfis/...'
+    verbose : bool, optional, default=False
+        verbose output if True
+
+    Raises
+    ------
+    HTTPError, KeyError
+
+    Returns
+    -------
+    vls_out : array of string
+        file or directory at path
+    """
+
+    cmd = 'vls'
+    sys.argv = [cmd, path]
+    vls = vosHandler(cmd)
+
+    if verbose:
+        print('Getting vos directory content from vls...')
+
+    f = StringIO()
+
+    try:
+        with redirect_stdout(f):
+            vls()
+    except:
+        print('Error during vls command')
+        raise
+
+    vls_out = f.getvalue()
+    print(vls_out)
+    return vls_out.split('\n')

From 977a6cc11c5a2b707cd52626dbecf3290bf622cb Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Thu, 10 Jun 2021 21:40:51 +0200
Subject: [PATCH 13/30] removed comments

---
 scripts/sh/job_sp.bash | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/scripts/sh/job_sp.bash b/scripts/sh/job_sp.bash
index 54a33b56c..00bb78ff7 100755
--- a/scripts/sh/job_sp.bash
+++ b/scripts/sh/job_sp.bash
@@ -361,7 +361,6 @@ if [[ $do_job != 0 ]]; then
 
   ### The following are very a bad hacks to get additional input file paths
   if [ "$psf" == "psfex" ]; then
-    #input_psfex=`find . -name star_split_ratio_80-*.psf | head -n 1`
     input_psfex=`ls -1 ./output/*/psfex_runner/output/star_split_ratio_80*.psf | head -n 1`
     command_sp "ln -s `dirname $input_psfex` input_psfex" "Link psfex output"
   else
@@ -369,11 +368,9 @@ if [[ $do_job != 0 ]]; then
     command_sp "ln -s `dirname $input_psf_mccd` input_psf_mccd" "Link MCCD output"
   fi
 
-  #input_split_exp=`find output -name flag-*.fits | head -n 1`
   input_split_exp=`ls -1 ./output/*/split_exp_runner/output/flag*.fits | head -n 1`
   command_sp "ln -s `dirname $input_split_exp` input_split_exp" "Link split_exp output"
 
-  #input_sextractor=`find . -name sexcat_sexcat-*.fits | head -n 1`
   input_sextractor=`ls -1 ./output/*/sextractor_runner/output/sexcat_sexcat*.fits | head -n 1`
   command_sp "ln -s `dirname $input_sextractor` input_sextractor" "Link sextractor output"
 

From aef4fb2a9a0f295bde482a9fd510a22f9a8b040d Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Thu, 10 Jun 2021 21:42:12 +0200
Subject: [PATCH 14/30] removed print line

---
 scripts/sh/psf_residuals.bash | 1 -
 1 file changed, 1 deletion(-)

diff --git a/scripts/sh/psf_residuals.bash b/scripts/sh/psf_residuals.bash
index 1c6d5d488..4929350b6 100755
--- a/scripts/sh/psf_residuals.bash
+++ b/scripts/sh/psf_residuals.bash
@@ -114,7 +114,6 @@ fi
 n_skipped=0
 n_created=0
 FILES=output/*/${runner}/output/${psfval_file_base}*
-echo ${FILES[@]}
 for val in ${FILES[@]}; do
     base=`basename $val`
     link_s "$pwd/$val" "$dir_individual/$base"

From 1ea15a6cac33ad56902ea714cf2e3a87d1885b86 Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Thu, 10 Jun 2021 21:43:31 +0200
Subject: [PATCH 15/30] conflict solved in canfar.py

---
 shapepipe/utilities/canfar.py | 15 ---------------
 1 file changed, 15 deletions(-)

diff --git a/shapepipe/utilities/canfar.py b/shapepipe/utilities/canfar.py
index 6f12d4b7a..3588a7c49 100644
--- a/shapepipe/utilities/canfar.py
+++ b/shapepipe/utilities/canfar.py
@@ -94,10 +94,6 @@ def __call__(self, *args, **kwargs):
             raise vosError('Error in VOs command: {}'
                            ''.format(self._command.__name__))
 
-<<<<<<< HEAD
-=======
-
->>>>>>> upstream/master
 def download(source, target, verbose=False):
     """Download file from vos.
 
@@ -116,7 +112,6 @@ def download(source, target, verbose=False):
         status, True/False or success/failure
     """
 
-<<<<<<< HEAD
     cmd = 'vcp'
 
     if not os.path.exists(target):
@@ -124,20 +119,12 @@ def download(source, target, verbose=False):
         if verbose:
             print('Downloading file {} to {}...'.format(source, target))
         vcp = vosHandler(cmd)
-=======
-    if not os.path.exists(target):
-        sys.argv = ['vcp', source, target]
-        if verbose:
-            print('Downloading file {} to {}...'.format(source, target))
-        vcp = vosHandler('vcp')
->>>>>>> upstream/master
         vcp()
         if verbose:
             print('Download finished.')
     else:
         if verbose:
             print('Target file {} exists, skipping download.'.format(target))
-<<<<<<< HEAD
 
 
 def dir_list(path, verbose=False):
@@ -181,5 +168,3 @@ def dir_list(path, verbose=False):
     vls_out = f.getvalue()
     print(vls_out)
     return vls_out.split('\n')
-=======
->>>>>>> upstream/master

From c82d584d1ac4fec337fad85a3158854d2006587f Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Thu, 10 Jun 2021 21:46:59 +0200
Subject: [PATCH 16/30] quick cleanup of cfis.py

---
 shapepipe/utilities/cfis.py | 25 +++----------------------
 1 file changed, 3 insertions(+), 22 deletions(-)

diff --git a/shapepipe/utilities/cfis.py b/shapepipe/utilities/cfis.py
index 4813562b6..20242362c 100644
--- a/shapepipe/utilities/cfis.py
+++ b/shapepipe/utilities/cfis.py
@@ -8,11 +8,9 @@
 :Authors: Martin Kilbinger
 
 :Date: 19/01/2018
-"""
-
-# Compability with python2.x for x>6
-from __future__ import print_function
 
+:Package: ShapePipe
+"""
 
 import re
 import sys
@@ -32,7 +30,6 @@
 
 from shapepipe.utilities.file_system import mkdir
 
-
 unitdef = 'degree'
 
 # Maybe define class for these constants?
@@ -60,7 +57,6 @@ def var_list(self, **kwds):
         return vars(self)
 
 
-
 class CfisError(Exception):
    """ MyError
 
@@ -107,7 +103,6 @@ def __init__(self, name, ra, dec, exp_time=-1, valid='Unknown'):
         else:
             self.valid = valid
 
-
     def cut(self, no_cuts=False):
         """Return True (False) if image does (not) need to be cut from selection.
 
@@ -156,8 +151,6 @@ def get_ID(self):
         else:
             return '{}.{}'.format(m[1], m[2])
 
-
-
     def print(self, file=sys.stdout, base_name=False, name_only=True, ID_only=False):
         """Print image information as ascii Table column
 
@@ -198,7 +191,6 @@ def print(self, file=sys.stdout, base_name=False, name_only=True, ID_only=False)
             print(' {:5d} {:8s}'.format(self.exp_time, self.valid), end='', file=file)
         print(file=file)
 
-
     def print_header(self, file=sys.stdout):
         """Print header for ascii Table output
 
@@ -388,7 +380,6 @@ def check_error_stop(ex_list, verbose=True, stop=False):
     else:
         s = sum([abs(i) for i in ex_list])
 
-
     # Evaluate exit codes
     if s > 0:
         n_ex = sum([1 for i in ex_list if i != 0])
@@ -475,7 +466,6 @@ def symlink(src, dst, verbose=False):
     os.symlink(src, dst)
 
 
-
 def print_color(color, txt, file=sys.stdout, end='\n'):
     """Print text with color. If not supported, print standard text.
 
@@ -712,7 +702,6 @@ def get_tile_coord_from_nixy(nix, niy):
     return ra, dec
 
 
-
 def get_tile_name(nix, niy, band, image_type='tile', input_format='full'):
     """Return tile name for given tile numbers.
 
@@ -789,7 +778,6 @@ def get_tile_number(tile_name):
     return nix, niy
 
 
-
 def get_log_file(path, verbose=False):
     """Return log file content
 
@@ -818,7 +806,6 @@ def get_log_file(path, verbose=False):
     return log
 
 
-
 def check_ra(ra):
     """Range check of right ascension.
 
@@ -841,7 +828,6 @@ def check_ra(ra):
     return 0
 
 
-
 def check_dec(dec):
     """Range check of declination.
 
@@ -863,7 +849,6 @@ def check_dec(dec):
     return 0
 
 
-
 def get_Angle(str_coord):
     """Return Angles ra, dec from coordinate string
 
@@ -888,7 +873,6 @@ def get_Angle(str_coord):
     return a_ra, a_dec
 
 
-
 def get_Angle_arr(str_coord, num=-1, wrap=True, verbose=False):
     """Return array of Angles from coordinate string
 
@@ -926,7 +910,6 @@ def get_Angle_arr(str_coord, num=-1, wrap=True, verbose=False):
     return angles
 
 
-
 def read_list(fname, col=None):
     """Read list from ascii file.
 
@@ -955,6 +938,7 @@ def read_list(fname, col=None):
         file_list = dat[col]
 
     file_list.sort()
+
     return file_list
 
 
@@ -1141,7 +1125,6 @@ def get_image_list(inp, band, image_type, col=None, input_format='full', verbose
     return img_list
 
 
-
 def exclude(f, exclude_list):
     """Return True if f is on exclude_list
 
@@ -1248,7 +1231,6 @@ def log_get_tile_nums(log):
     return set(tile_nums)
 
 
-
 def log_get_exp_num(log, exp_name, k_img, k_weight, k_flag):
     """Return exposure number from log file for given exposure name and HDU numbers.
 
@@ -1271,7 +1253,6 @@ def log_get_exp_num(log, exp_name, k_img, k_weight, k_flag):
         exposure number, None if not found
     """
 
-
     for line in log:
         this_exp_name = log_line_get_entry(line, 'exp_name')
         this_k_img    = log_line_get_entry(line, 'k_img')

From 571fb3ddc5daf17bd9fd27a063ae05cdca3af6ae Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Thu, 10 Jun 2021 22:13:58 +0200
Subject: [PATCH 17/30] passing tests

---
 shapepipe/utilities/canfar.py |   3 +-
 shapepipe/utilities/cfis.py   | 212 +++++++++++++++++-----------------
 2 files changed, 106 insertions(+), 109 deletions(-)

diff --git a/shapepipe/utilities/canfar.py b/shapepipe/utilities/canfar.py
index 3588a7c49..375697cc9 100644
--- a/shapepipe/utilities/canfar.py
+++ b/shapepipe/utilities/canfar.py
@@ -9,7 +9,7 @@
 """
 
 import os
-import sys  
+import sys
 from io import StringIO
 from contextlib import redirect_stdout
 
@@ -94,6 +94,7 @@ def __call__(self, *args, **kwargs):
             raise vosError('Error in VOs command: {}'
                            ''.format(self._command.__name__))
 
+
 def download(source, target, verbose=False):
     """Download file from vos.
 
diff --git a/shapepipe/utilities/cfis.py b/shapepipe/utilities/cfis.py
index 20242362c..fb58b4372 100644
--- a/shapepipe/utilities/cfis.py
+++ b/shapepipe/utilities/cfis.py
@@ -34,13 +34,13 @@
 
 # Maybe define class for these constants?
 size = {}
-size['tile']     = 0.5
-size['weight']   = 0.5
+size['tile'] = 0.5
+size['weight'] = 0.5
 size['exposure'] = 1.0
 
 # Cut criteria for exposures
-exp_time_min     = 95
-flag_valid       = 'V'
+exp_time_min = 95
+flag_valid = 'V'
 
 
 class param:
@@ -91,14 +91,14 @@ def __init__(self, name, ra, dec, exp_time=-1, valid='Unknown'):
             image information
         """
 
-        self.name     = name
-        self.ra       = ra
-        self.dec      = dec
-        if exp_time == None:
+        self.name = name
+        self.ra = ra
+        self.dec = dec
+        if exp_time is None:
             self.exp_time = -1
         else:
             self.exp_time = exp_time
-        if valid == None:
+        if valid is None:
             self.valid = 'Unknown'
         else:
             self.valid = valid
@@ -118,7 +118,7 @@ def cut(self, no_cuts=False):
         """
 
         # Do not cut if no_cuts flag is set
-        if no_cuts == True:
+        if no_cuts:
             return False
 
         # Cut if exposure time smaller than minimum (and not flagged as unknown or n/a)
@@ -145,7 +145,7 @@ def get_ID(self):
         NoneType : if name does not match to ID pattern
         """
 
-        m = re.search('(\d{3}).{1}(\d{3})', self.name)
+        m = re.search(r'(\d{3}).{1}(\d{3})', self.name)
         if m is None:
             raise NoneType('No ID match in file name {}'.format(name))
         else:
@@ -176,7 +176,7 @@ def print(self, file=sys.stdout, base_name=False, name_only=True, ID_only=False)
             name = self.name
 
         if ID_only:
-            m = re.search('\d{3}.\d{3}', name)
+            m = re.search(r'\d{3}.\d{3}', name)
             if m is None:
                 raise NoneType('No ID match in file name {}'.format(name))
             else:
@@ -262,11 +262,10 @@ def run_cmd_old(cmd_list, run=True, verbose=True, stop=False, parallel=True, fil
     if verbose is True and len(cmd_list) > 1:
         print('Running {} commands, parallel = {}'.format(len(cmd_list), parallel))
 
-
-    ex_list   = []
+    ex_list = []
     pipe_list = []
-    out_list  = []
-    err_list  = []
+    out_list = []
+    err_list = []
 
     if env is None:
         env = os.environ.copy()
@@ -284,7 +283,7 @@ def run_cmd_old(cmd_list, run=True, verbose=True, stop=False, parallel=True, fil
                     print_color('blue', 'Skipping command \'{}\', file \'{}\' exists'.format(cmd, file_list[i]))
             else:
                 if verbose is True:
-                        print_color('green', 'Running command \'{0}\''.format(cmd))
+                    print_color('green', 'Running command \'{0}\''.format(cmd))
 
                 # Run command
                 try:
@@ -322,7 +321,6 @@ def run_cmd_old(cmd_list, run=True, verbose=True, stop=False, parallel=True, fil
         out_list.append(out)
         err_list.append(err)
 
-
     if parallel is True:
         for i, pipe in enumerate(pipe_list):
             pipe.wait()
@@ -396,7 +394,6 @@ def check_error_stop(ex_list, verbose=True, stop=False):
         if stop is True:
             sys.exit(s)
 
-
     return s
 
 
@@ -441,7 +438,7 @@ def log_command(argv, name=None, close_no_return=True):
 
     print('', file=f)
 
-    if close_no_return == False:
+    if not close_no_return:
         return f
 
     if name != 'sys.stdout' and name != 'sys.stderr':
@@ -487,12 +484,11 @@ def print_color(color, txt, file=sys.stdout, end='\n'):
 
     try:
         import colorama
-        colors = {'red'    : colorama.Fore.RED,
-                  'green'  : colorama.Fore.GREEN,
-                  'blue'   : colorama.Fore.BLUE,
-                  'yellow' : colorama.Fore.YELLOW,
-                  'black'  : colorama.Fore.BLACK,
-                 }
+        colors = {'red': colorama.Fore.RED,
+                  'green': colorama.Fore.GREEN,
+                  'blue': colorama.Fore.BLUE,
+                  'yellow': colorama.Fore.YELLOW,
+                  'black': colorama.Fore.BLACK}
 
         if colors[color] is None:
             col = colorama.Fore.BLACK
@@ -538,9 +534,9 @@ def my_string_split(string, num=-1, verbose=False, stop=False, sep=None):
         return None
 
     if sep is None:
-        has_space      = string.find(' ')
+        has_space = string.find(' ')
         has_underscore = string.find('_')
-        has_dot        = string.find('.')
+        has_dot = string.find('.')
 
         if has_space != -1:
             my_sep = ' '
@@ -561,7 +557,7 @@ def my_string_split(string, num=-1, verbose=False, stop=False, sep=None):
 
     res = string.split(my_sep)
 
-    if num != -1 and num != len(res) and stop==True:
+    if num != -1 and num != len(res) and stop:
         raise CfisError('String \'{}\' has length {}, required is {}'.format(string, len(res), num))
 
     return res
@@ -592,40 +588,39 @@ def get_file_pattern(pattern, band, image_type, want_re=True, ext=True):
     """
 
     if pattern == '':
-        if image_type in ('exposure', 'exposure_flag', 'exposure_flag.fz', \
-		'exposure_weight', 'exposure_weight.fz'):
-            pattern_base = '\d{7}p'
+        if image_type in ('exposure', 'exposure_flag', 'exposure_flag.fz',
+                          'exposure_weight', 'exposure_weight.fz'):
+            pattern_base = r'\d{7}p'
         else:
-            pattern_base  = 'CFIS.*\.{}'.format(band)
+            pattern_base = r'CFIS.*\.{}'.format(band)
     else:
         pattern_base = pattern
 
-
     if ext:
         if image_type == 'exposure':
-            pattern_out  = '{}\.fits\.fz'.format(pattern_base)
+            pattern_out = r'{}\.fits\.fz'.format(pattern_base)
         elif image_type == 'exposure_flag':
-            pattern_out  = '{}\.flag\.fits'.format(pattern_base)
+            pattern_out = r'{}\.flag\.fits'.format(pattern_base)
         elif image_type == 'exposure_flag.fz':
-            pattern_out  = '{}\.flag\.fits\.fz'.format(pattern_base)
+            pattern_out = r'{}\.flag\.fits\.fz'.format(pattern_base)
         elif image_type == 'exposure_weight':
-            pattern_out  = '{}\.weight\.fits'.format(pattern_base)
+            pattern_out = r'{}\.weight\.fits'.format(pattern_base)
         elif image_type == 'exposure_weight.fz':
-            pattern_out  = '{}\.weight\.fits\.fz'.format(pattern_base)
+            pattern_out = r'{}\.weight\.fits\.fz'.format(pattern_base)
         elif image_type == 'tile':
-            pattern_out = '{}\.fits'.format(pattern_base)
+            pattern_out = r'{}\.fits'.format(pattern_base)
         elif image_type == 'cat':
-            pattern_out = '{}\.cat'.format(pattern_base)
+            pattern_out = r'{}\.cat'.format(pattern_base)
         elif image_type == 'weight':
-            pattern_out = '{}\.weight\.fits'.format(pattern_base)
+            pattern_out = r'{}\.weight\.fits'.format(pattern_base)
         elif image_type == 'weight.fz':
-            pattern_out = '{}\.weight\.fits\.fz'.format(pattern_base)
+            pattern_out = r'{}\.weight\.fits\.fz'.format(pattern_base)
         else:
             raise CfisError('Invalid type \'{}\''.format(image_type))
     else:
         pattern_out = pattern_base
 
-    if want_re == False:
+    if not want_re:
         pattern_out = pattern_out.replace('\\', '')
 
     return pattern_out
@@ -642,7 +637,7 @@ def get_tile_number_from_coord(ra, dec, return_type=str):
     dec: Angle
         declination
     return type: <type 'type'>
-	return type, int or str
+        return type, int or str
 
     Returns
     -------
@@ -656,7 +651,6 @@ def get_tile_number_from_coord(ra, dec, return_type=str):
     yi = int(np.rint(y))
 
     x = ra.degree * np.cos(dec.radian) * 2.0
-    #x = ra.degree * 2 * np.cos(y/2 / 180 * np.pi - np.pi/2)
     xi = int(np.rint(x))
     if xi == 720:
         xi = 0
@@ -728,12 +722,12 @@ def get_tile_name(nix, niy, band, image_type='tile', input_format='full'):
         if input_format == 'ID_only':
             tile_base = '{:03d}.{:03d}'.format(nix, niy)
         else:
-    	    tile_base = 'CFIS.{:03d}.{:03d}.{}'.format(nix, niy, band)
+            tile_base = 'CFIS.{:03d}.{:03d}.{}'.format(nix, niy, band)
     elif type(nix) is str and type(niy) is str:
         if input_format == 'ID_only':
             tile_base = '{}.{}'.format(nix, niy)
         else:
-    	    tile_base = 'CFIS.{}.{}.{}'.format(nix, niy, band)
+            tile_base = 'CFIS.{}.{}.{}'.format(nix, niy, band)
     else:
         raise CfisError('Invalid type for input tile numbers {}, {}'.format(nix, niy))
 
@@ -768,8 +762,8 @@ def get_tile_number(tile_name):
         tile number for y
     """
 
-    m = re.search('(\d{3})[\.-](\d{3})', tile_name)
-    if m == None or len(m.groups()) != 2:
+    m = re.search(r'(\d{3})[\.-](\d{3})', tile_name)
+    if m is None or len(m.groups()) != 2:
         raise CfisError('Image name \'{}\' does not match tile name syntax'.format(tile_name))
 
     nix = m.groups()[0]
@@ -798,7 +792,7 @@ def get_log_file(path, verbose=False):
         raise CfisError('Log file \'{}\' not found'.format(path))
 
     f_log = open(path, 'r')
-    log   = f_log.readlines()
+    log = f_log.readlines()
     if verbose:
         print('Reading log file, {} lines found'.format(len(log)))
     f_log.close()
@@ -867,7 +861,7 @@ def get_Angle(str_coord):
 
     ra, dec = my_string_split(str_coord, num=2, stop=True)
 
-    a_ra  = Angle(ra)
+    a_ra = Angle(ra)
     a_dec = Angle(dec)
 
     return a_ra, a_dec
@@ -904,7 +898,7 @@ def get_Angle_arr(str_coord, num=-1, wrap=True, verbose=False):
         if wrap:
             c = SkyCoord(ra, dec)
         else:
-            c = param(ra = Angle(ra), dec = Angle(dec))
+            c = param(ra=Angle(ra), dec=Angle(dec))
         angles.append(c)
 
     return angles
@@ -932,7 +926,7 @@ def read_list(fname, col=None):
         f.close()
     else:
         import pandas as pd
-        dat = pd.read_csv(fname, sep='\s+', dtype='string', header=None)
+        dat = pd.read_csv(fname, sep=r'\s+', dtype='string', header=None)
         if col not in dat:
             col = int(col)
         file_list = dat[col]
@@ -1014,8 +1008,8 @@ def get_exposure_info(logfile_name, verbose=False):
     for line in f:
         dat = re.split(' |', line)
         name = dat[0]
-        ra   = Angle(' hours'.format(dat[8]))
-        dec  = Angle(' degree'.format(dat[9]))
+        ra = Angle(' hours'.format(dat[8]))
+        dec = Angle(' degree'.format(dat[9]))
         valid = dat[21]
 
         img = image(name, ra, dec, valid=valid)
@@ -1048,47 +1042,47 @@ def get_image_list(inp, band, image_type, col=None, input_format='full', verbose
         image list
     """
 
-    file_list     = []
-    ra_list       = []
-    dec_list      = []
+    file_list = []
+    ra_list = []
+    dec_list = []
     exp_time_list = []
-    valid_list    = []
+    valid_list = []
 
     if os.path.isdir(inp):
         if col is not None:
             raise CfisError('Column name (-c option) only valid if input is file')
 
         # Read file names from directory listing
-        inp_type  = 'dir'
+        inp_type = 'dir'
         file_list = glob.glob('{}/*'.format(os.path.abspath(inp)))
 
     elif os.path.isfile(inp):
         if image_type in ('tile', 'weight', 'weight.fz'):
             # File names in single-column ascii file
-            inp_type  = 'file'
+            inp_type = 'file'
             file_list = read_list(inp, col=col)
         elif image_type == 'exposure':
             # File names and coordinates in ascii file
-            inp_type  = 'file'
+            inp_type = 'file'
             dat = ascii.read(inp)
 
             if len(dat.keys()) == 3:
                 # File is exposure + coord list (obtained from get_coord_CFIS_pointings.py)
                 file_list = dat['Pointing']
-                ra_list   = dat['R.A.[degree]']
-                dec_list  = dat['Declination[degree]']
+                ra_list = dat['R.A.[degree]']
+                dec_list = dat['Declination[degree]']
             elif len(dat.keys()) == 12:
                 # File is log file, e.g. from http://www.cfht.hawaii.edu/Science/CFIS-DATA/logs/MCLOG-CFIS.r.qso-elx.log
                 # Default file separator is '|'
                 for d in dat:
                     file_list.append('d{}p.fits.fz'.format(d['col1']))
-                    ra  = re.split('\s*', d['col4'])[0]
-                    dec = re.split('\s*', d['col4'])[1]
+                    ra = re.split(r'\s*', d['col4'])[0]
+                    dec = re.split(r'\s*', d['col4'])[1]
                     ra_list.append(Angle('{} hours'.format(ra)).degree)
                     dec_list.append(dec)
                     exp_time = int(d['col5'])
                     exp_time_list.append(exp_time)
-                    valid = re.split('\s*', d['col11'])[2]
+                    valid = re.split(r'\s*', d['col11'])[2]
                     valid_list.append(valid)
             else:
                 raise CfisError('Wrong file format, #columns={}, has to be 3 or 12'.format(len(dat.keys())))
@@ -1101,9 +1095,9 @@ def get_image_list(inp, band, image_type, col=None, input_format='full', verbose
     # Filter file list to match CFIS image pattern
     img_list = []
     if input_format == 'ID_only':
-        pattern = get_file_pattern('\d{3}.\d{3}', band, image_type, ext=False)
+        pattern = get_file_pattern(r'\d{3}.\d{3}', band, image_type, ext=False)
     else:
-        pattern = get_file_pattern('CFIS.\d{{3}}.\d{{3}}\.{}'.format(band), band, image_type)
+        pattern = get_file_pattern(r'CFIS.\d{{3}}.\d{{3}}\.{}'.format(band), band, image_type)
 
     for img in image_list:
 
@@ -1119,7 +1113,7 @@ def get_image_list(inp, band, image_type, col=None, input_format='full', verbose
         if len(m) != 0:
             img_list.append(img)
 
-    if verbose == True and len(img_list) > 0:
+    if verbose and len(img_list) > 0:
         print('{} image files found in input {} \'{}\''.format(len(img_list), inp_type, inp))
 
     return img_list
@@ -1255,14 +1249,14 @@ def log_get_exp_num(log, exp_name, k_img, k_weight, k_flag):
 
     for line in log:
         this_exp_name = log_line_get_entry(line, 'exp_name')
-        this_k_img    = log_line_get_entry(line, 'k_img')
+        this_k_img = log_line_get_entry(line, 'k_img')
         this_k_weight = log_line_get_entry(line, 'k_weight')
-        this_k_flag   = log_line_get_entry(line, 'k_flag')
+        this_k_flag = log_line_get_entry(line, 'k_flag')
 
-        if this_exp_name == exp_name  and \
-            int(this_k_img) == k_img and \
-            int(this_k_weight) == k_weight and \
-            int(this_k_flag) == k_flag:
+        if this_exp_name == exp_name \
+                and int(this_k_img) == k_img \
+                and int(this_k_weight) == k_weight \
+                and int(this_k_flag) == k_flag:
             return log_line_get_entry(line, 'exp_num')
 
     # No matching entry found
@@ -1298,11 +1292,11 @@ def find_image_at_coord(images, coord, band, image_type, no_cuts=False, input_fo
 
     ra, dec = get_Angle(coord)
 
-    if verbose == True:
+    if verbose:
         print('Looking for image at coordinates {}, {}'.format(ra, dec))
 
     if image_type in ('tile', 'weight', 'weight.fz'):
-        nix, niy  = get_tile_number_from_coord(ra, dec, return_type=int)
+        nix, niy = get_tile_number_from_coord(ra, dec, return_type=int)
         tile_name = get_tile_name(nix, niy, band, image_type, input_format=input_format)
 
         img_found = []
@@ -1312,16 +1306,16 @@ def find_image_at_coord(images, coord, band, image_type, no_cuts=False, input_fo
                 ra_c, dec_c = get_tile_coord_from_nixy(nix, niy)
                 if img.ra is not None or img.dec is not None:
                     raise CfisError('Coordinates in image are already '
-                                         'set to {}, {}, cannot update to {}, {}'.\
-                                         format(img.ra, img.dec, ra_c, dec_c))
+                                    'set to {}, {}, cannot update to {}, {}'
+                                    ''.format(img.ra, img.dec, ra_c, dec_c))
                 img.ra = ra_c
                 img.dec = dec_c
                 img_found.append(img)
 
         if len(img_found) != 0:
-                pass
+            pass
         else:
-            if verbose == True:
+            if verbose:
                 print('Tile with numbers ({}, {}) not found'.format(nix, niy))
 
         if len(img_found) > 1:
@@ -1333,18 +1327,18 @@ def find_image_at_coord(images, coord, band, image_type, no_cuts=False, input_fo
         img_found = []
         for img in images:
             # Check distance along ra and dec from image center
-            sc_img_same_ra  = SkyCoord(ra, img.dec)
+            sc_img_same_ra = SkyCoord(ra, img.dec)
             sc_img_same_dec = SkyCoord(img.ra, dec)
-            distance_ra  = sc_input.separation(sc_img_same_dec)
+            distance_ra = sc_input.separation(sc_img_same_dec)
             distance_dec = sc_input.separation(sc_img_same_ra)
             if distance_ra.degree < size[image_type]/2 and distance_dec.degree < size[image_type]/2:
-                if img.cut(no_cuts=no_cuts) == False:
+                if not img.cut(no_cuts=no_cuts):
                     img_found.append(img)
 
         if len(img_found) != 0:
-                pass
+            pass
         else:
-            if verbose == True:
+            if verbose:
                 print('No exposure image found')
 
     else:
@@ -1390,7 +1384,7 @@ def find_images_in_area(images, angles, band, image_type, no_cuts=False, verbose
     else:
         ra_bounds = [[angles[0].ra, angles[1].ra]]
 
-    if verbose == True:
+    if verbose:
         print('Looking for all images within rectangle, '
               'dec=({}, {}), '
               ''.format(angles[0].dec, angles[1].dec), end='')
@@ -1401,9 +1395,9 @@ def find_images_in_area(images, angles, band, image_type, no_cuts=False, verbose
     if image_type in ('tile', 'weight', 'weight.fz'):
         for img in images:
             nix, niy = get_tile_number(img.name)
-            ra, dec  = get_tile_coord_from_nixy(nix, niy)
+            ra, dec = get_tile_coord_from_nixy(nix, niy)
 
-            if  dec.is_within_bounds(angles[0].dec, angles[1].dec):
+            if dec.is_within_bounds(angles[0].dec, angles[1].dec):
                 within = False
 
                 # Check whether image is in any of the ra bound pairs
@@ -1413,14 +1407,14 @@ def find_images_in_area(images, angles, band, image_type, no_cuts=False, verbose
                         break
                 if within:
                     if img.ra is None or img.dec is None:
-                        img.ra  = ra
+                        img.ra = ra
                         img.dec = dec
 
                     found.append(img)
 
     elif image_type == 'exposure':
         for img in images:
-            if  img.dec.is_within_bounds(angles[0].dec, angles[1].dec):
+            if img.dec.is_within_bounds(angles[0].dec, angles[1].dec):
                 within = False
 
                 # Check whether image is in any of the ra bound pairs
@@ -1429,13 +1423,13 @@ def find_images_in_area(images, angles, band, image_type, no_cuts=False, verbose
                         within = True
                         break
                 if within:
-                    if img.cut(no_cuts=no_cuts) == False:
+                    if not img.cut(no_cuts=no_cuts):
                         found.append(img)
 
     else:
         raise CfisError('Image type \'{}\' not implemented yet'.format(image_type))
 
-    if verbose == True:
+    if verbose:
         print('{} images found in area'.format(len(found)))
 
     return found
@@ -1458,7 +1452,7 @@ def plot_init():
 
 
 def plot_area(images, angles, image_type, outbase, interactive, col=None, show_numbers=False,
-              show_circle=True, ax=None, lw=None, save=True, dxy=0):
+              show_circle=True, show_area_border=False, ax=None, lw=None, save=True, dxy=0):
     """Plot images within area.
 
     Parameters
@@ -1477,6 +1471,8 @@ def plot_area(images, angles, image_type, outbase, interactive, col=None, show_n
         color
     show_circle : bool, optional, default True
         plot circle center and circumference around area if True
+    show_area_border : bool, optional, default False
+        plot rectangle around area
     ax : axes, optional, default None
         Init axes if None
     lw : float, optional, default None
@@ -1504,7 +1500,7 @@ def plot_area(images, angles, image_type, outbase, interactive, col=None, show_n
     # Field center
     n_ima = len(images)
     if n_ima > 0:
-        ra_c  = sum([img.ra for img in images])/float(n_ima)
+        ra_c = sum([img.ra for img in images])/float(n_ima)
         dec_c = sum([img.dec for img in images])/float(n_ima)
         cos_dec_c = np.cos(dec_c)
     else:
@@ -1526,18 +1522,18 @@ def plot_area(images, angles, image_type, outbase, interactive, col=None, show_n
         radius = 0
 
     if col:
-        c  = col
+        c = col
     else:
         c = color[image_type]
 
     for img in images:
-        x  = img.ra.degree
-        y  = img.dec.degree
+        x = img.ra.degree
+        y = img.dec.degree
         nix, niy = get_tile_number(img.name)
         if show_numbers:
             plt.text(x, y, '{}.{}'.format(nix, niy), fontsize=3,
-                    horizontalalignment='center',
-                    verticalalignment='center')
+                     horizontalalignment='center',
+                     verticalalignment='center')
 
         # Image boundary
         dx = size[image_type] / 2 / cos_dec_c
@@ -1545,9 +1541,9 @@ def plot_area(images, angles, image_type, outbase, interactive, col=None, show_n
         cx, cy = square_from_centre(x, y, dx, dy, dxy=dxy)
         ax.plot(cx, cy, '-', color=c, linewidth=my_lw)
 
-    # Area border
-    #cx, cy = square_from_corners(angles[0], angles[1])
-    #ax.plot(cx, cy, 'r-.', linewidth=my_lw)
+    if show_area_border:
+        cx, cy = square_from_corners(angles[0], angles[1])
+        ax.plot(cx, cy, 'r-.', linewidth=my_lw)
 
     plt.xlabel('R.A. [degree]')
     plt.ylabel('Declination [degree]')
@@ -1560,8 +1556,8 @@ def plot_area(images, angles, image_type, outbase, interactive, col=None, show_n
         xm = (angles[1].ra.degree + angles[0].ra.degree) / 2
         dx = angles[1].ra.degree - angles[0].ra.degree
     else:
-        dx = max(360 - angles[0].ra.deg,  angles[1].ra.deg) * 2
-        xm = ( (angles[0].ra.deg - 360) + angles[1].ra.deg ) / 2
+        dx = max(360 - angles[0].ra.deg, angles[1].ra.deg) * 2
+        xm = ((angles[0].ra.deg - 360) + angles[1].ra.deg) / 2
     ym = (angles[1].dec.degree + angles[0].dec.degree) / 2
     dy = angles[1].dec.degree - angles[0].dec.degree
     lim = max(dx, dy)
@@ -1572,7 +1568,7 @@ def plot_area(images, angles, image_type, outbase, interactive, col=None, show_n
         print('Saving plot to {}'.format(outname))
         plt.savefig(outname)
 
-    if interactive == True:
+    if interactive:
         plt.show()
 
     return ra_c, dec_c, radius

From f1c414db084c4ddcc5448ef7eaeae6a2cbed5d60 Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Tue, 15 Jun 2021 15:16:03 +0200
Subject: [PATCH 18/30] removed comment

---
 shapepipe/utilities/canfar.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/shapepipe/utilities/canfar.py b/shapepipe/utilities/canfar.py
index 375697cc9..15ea051c9 100644
--- a/shapepipe/utilities/canfar.py
+++ b/shapepipe/utilities/canfar.py
@@ -167,5 +167,5 @@ def dir_list(path, verbose=False):
         raise
 
     vls_out = f.getvalue()
-    print(vls_out)
+
     return vls_out.split('\n')

From c9af7be7cd982641471b5aecd5aa16963a1380ab Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Tue, 15 Jun 2021 16:34:33 +0200
Subject: [PATCH 19/30] first revision spread model module

---
 .../SpreadModel_package/SpreadModel_script.py | 260 ++++++++++++++++++
 shapepipe/modules/spread_model_runner.py      | 244 ++--------------
 shapepipe/utilities/__init__.py               |   2 +-
 shapepipe/utilities/galaxy.py                 |  36 +++
 4 files changed, 315 insertions(+), 227 deletions(-)
 create mode 100644 shapepipe/modules/SpreadModel_package/SpreadModel_script.py
 create mode 100644 shapepipe/utilities/galaxy.py

diff --git a/shapepipe/modules/SpreadModel_package/SpreadModel_script.py b/shapepipe/modules/SpreadModel_package/SpreadModel_script.py
new file mode 100644
index 000000000..65fc0b256
--- /dev/null
+++ b/shapepipe/modules/SpreadModel_package/SpreadModel_script.py
@@ -0,0 +1,260 @@
+# -*- coding: utf-8 -*-
+
+"""SPREAD MODEL SCRIPT
+
+Class to compute the spread model, criterion to select galaxies
+
+:Author: Axel Guinot
+
+:Date: 2019, 2020
+
+:Package: ShapePipe
+
+"""
+
+
+import numpy as np
+import galsim
+from sqlitedict import SqliteDict
+
+from shapepipe.pipeline import file_io as io
+from shapepipe.utilities import galaxy
+
+
+def get_sm(obj_vign, psf_vign, model_vign, weight_vign):
+    """ Get Spread model
+
+    This method compute the spread moel for an object.
+
+    Parameters
+    ----------
+    obj_vign : numpy.ndarray
+        Vignet of the object.
+    psf_vign : numpy.ndarray
+        Vignet of the gaussian model of the PSF.
+    model_vign : numpy.ndarray
+        Vignet of the galaxy model.
+    weight_vign : numpy.ndarray
+        Vignet of the weight at the object position.
+
+    Returns
+    -------
+    sm : float
+        Spread model value.
+    sm_err : float
+        Spread model error value.
+    """
+
+    m = (obj_vign > -1e29) & (weight_vign > 0)
+    w = m.astype(float)
+
+    t_v = model_vign.ravel()
+    g_v = obj_vign.ravel()
+    psf_v = psf_vign.ravel()
+
+    noise_v = (1 / weight_vign).ravel()
+    noise_v[np.isinf(noise_v)] = 0
+    w_v = w.ravel()
+
+    tg = np.sum(t_v*w_v*g_v)
+    pg = np.sum(psf_v*w_v*g_v)
+    tp = np.sum(t_v*w_v*psf_v)
+    pp = np.sum(psf_v*w_v*psf_v)
+
+    tnt = np.sum(t_v*noise_v*t_v*w_v)
+    pnp = np.sum(psf_v*noise_v*psf_v*w_v)
+    tnp = np.sum(t_v*noise_v*psf_v*w_v)
+    err = tnt * pg**2 + pnp * tg**2 - 2*tnp * pg * tg
+
+    if pg > 0:
+        sm = (tg / pg) - (tp / pp)
+    else:
+        sm = 1
+
+    if (pg > 0) & (err > 0):
+        sm_err = np.sqrt(err) / pg**2
+    else:
+        sm_err = 1
+
+    return sm, sm_err
+
+
+def get_model(sigma, flux, img_shape, pixel_scale=0.186):
+    """ Get model
+
+    This method computes :
+     - an exponential galaxy model with scale radius = 1/16 FWHM
+     - a gaussian model for the PSF
+
+    Parameters
+    ----------
+    sigma : float
+        Sigma of the PSF (in pixel units).
+    flux : float
+        Flux of the galaxy for the model.
+    img_shape : list
+        Size of the output vignet [xsize, ysize].
+    pixel_scale : float, optional, default=0.186
+        Pixel scale to use for the model (in arcsec)
+
+    Returns
+    -------
+    gal_vign : numpy.ndarray
+        Vignet of the galaxy model.
+    psf_vign : numpy.ndarray
+        Vignet of the PSF model.
+    """
+
+    scale_radius = 1/16 * galaxy.sigma_to_fhwm(sigma, pixel_scale=pixel_scale)
+    gal_obj = galsim.Exponential(scale_radius=scale_radius, flux=flux)
+
+    psf_obj = galsim.Gaussian(sigma=sigma*pixel_scale)
+
+    gal_obj = galsim.Convolve(gal_obj, psf_obj)
+
+    gal_vign = gal_obj.drawImage(nx=img_shape[0], ny=img_shape[1],
+                                 scale=pixel_scale).array
+    psf_vign = psf_obj.drawImage(nx=img_shape[0], ny=img_shape[1],
+                                 scale=pixel_scale).array
+
+    return gal_vign, psf_vign
+
+
+class SpreadModel(object):
+    """SpreadModel initialisation.
+
+    Parameters
+    ----------
+    sex_cat_path : string
+        path to SExtractor catalogue
+    psf_cat_path : string
+        path to PSF catalogue
+    weight_cat_path : string
+        path to weight catalogue
+    output_path : string
+        output file path of pasted catalog
+    w_log :
+        log file
+    pixel_scale : float
+        pixel scale in arcsec
+    output_mode : str
+        must be in ['new', 'add'].
+         'new' will create a new catalog with : [number, mag, sm, sm_err]
+         'add' will output a copy of the input SExtractor with the column sm
+          and sm_err.
+    """
+
+    def __init__(self, sex_cat_path, psf_cat_path, weight_cat_path,
+                 output_path, w_log, pixel_scale, output_mode):
+
+        self._sex_cat_path = sex_cat_path
+        self._psf_cat_path = psf_cat_path
+        self._weight_cat_path = weight_cat_path
+        self._output_path = output_path
+        self._w_log = w_log
+        self._pixel_scale = pixel_scale
+        self._output_mode = output_mode
+
+    def process(self):
+        """Process
+        Process the spread model computation
+        """
+
+        # Get data
+        sex_cat = io.FITSCatalog(sex_cat_path, SEx_catalog=True)
+        sex_cat.open()
+        obj_id = np.copy(sex_cat.get_data()['NUMBER'])
+        obj_vign = np.copy(sex_cat.get_data()['VIGNET'])
+        obj_mag = None
+        if output_mode == 'new':
+            obj_mag = np.copy(sex_cat.get_data()['MAG_AUTO'])
+        sex_cat.close()
+
+        psf_cat = SqliteDict(psf_cat_path)
+
+        weight_cat = io.FITSCatalog(weight_cat_path, SEx_catalog=True)
+        weight_cat.open()
+        weigh_vign = weight_cat.get_data()['VIGNET']
+        weight_cat.close()
+
+        # Get spread model
+        skip_obj = False
+        spread_model_final = []
+        spread_model_err_final = []
+        for i, id_tmp in enumerate(obj_id):
+            sigma_list = []
+
+            if psf_cat[str(id_tmp)] == 'empty':
+                spread_model_final.append(-1)
+                spread_model_err_final.append(1)
+                continue
+
+            psf_expccd_name = list(psf_cat[str(id_tmp)].keys())
+
+            for expccd_name_tmp in psf_expccd_name:
+                sigma_list.append(psf_cat[str(id_tmp)][expccd_name_tmp]['SHAPES']['SIGMA_PSF_HSM'])
+
+            obj_vign_tmp = obj_vign[i]
+            obj_flux_tmp = 1.
+            obj_sigma_tmp = np.mean(sigma_list)
+            obj_weight_tmp = weigh_vign[i]
+            obj_model_tmp, obj_psf_tmp = get_model(obj_sigma_tmp,
+                                                   obj_flux_tmp,
+                                                   obj_vign_tmp.shape,
+                                                   pixel_scale)
+
+            obj_sm, obj_sm_err = get_sm(obj_vign_tmp,
+                                        obj_psf_tmp,
+                                        obj_model_tmp,
+                                        obj_weight_tmp)
+
+            spread_model_final.append(obj_sm)
+            spread_model_err_final.append(obj_sm_err)
+
+        spread_model_final = np.array(spread_model_final)
+        spread_model_err_final = np.array(spread_model_err_final)
+
+        psf_cat.close()
+
+        self.save_results(spread_model_final,
+                          spread_model_err_final,
+                          obj_mag,
+                          obj_id)
+
+    def save_results(self, sm, sm_err, mag, number):
+        """Save Results
+
+        Save output catalogue with spread model and errors.
+
+        Parameters
+        ----------
+        sm : numpy.ndarray
+            Value of the spread model for all objects.
+        sm_err : numpy.ndarray
+            Value of the spread model error for all objects.
+        mag : numpy.ndarray
+            Magnitude of all objects (only for new catalog).
+        number : numpy.ndarray
+            Id of all objects (only for new catalog).
+        """
+
+        if self.output_mode == 'new':
+            new_cat = io.FITSCatalog(self.output_path, SEx_catalog=True,
+                                     open_mode=io.BaseCatalog.OpenMode.ReadWrite)
+            dict_data = {'NUMBER': number,
+                         'MAG': mag,
+                         'SPREAD_MODEL': sm,
+                         'SPREADERR_MODEL': sm_err}
+            new_cat.save_as_fits(data=dict_data, sex_cat_path=self.sex_cat_path)
+        elif self.output_mode == 'add':
+            ori_cat = io.FITSCatalog(self.sex_cat_path, SEx_catalog=True)
+            ori_cat.open()
+            new_cat = io.FITSCatalog(self.output_path, SEx_catalog=True,
+                                     open_mode=io.BaseCatalog.OpenMode.ReadWrite)
+            ori_cat.add_col('SPREAD_MODEL', sm, new_cat=True, new_cat_inst=new_cat)
+            ori_cat.close()
+            new_cat.open()
+            new_cat.add_col('SPREADERR_MODEL', sm_err)
+            new_cat.close()
+        else:
+            ValueError('Mode must be in [new, add].')
diff --git a/shapepipe/modules/spread_model_runner.py b/shapepipe/modules/spread_model_runner.py
index 58db61149..98a63fab7 100644
--- a/shapepipe/modules/spread_model_runner.py
+++ b/shapepipe/modules/spread_model_runner.py
@@ -1,177 +1,27 @@
 # -*- coding: utf-8 -*-
 
-"""SPREAD MODEL RU+NNER
+"""SPREAD MODEL RUNNER
 
-This module compute the spread model.
+This module computes the spread model.
 
 :Author: Axel Guinot
 
-"""
+:Date: 2019, 2020
 
-import numpy as np
-import galsim
+"""
 
 from shapepipe.modules.module_decorator import module_runner
-from shapepipe.pipeline import file_io as io
-from sqlitedict import SqliteDict
-
-
-def get_sm(obj_vign, psf_vign, model_vign, weight_vign):
-    """ Get Spread model
-
-    This method compute the spread moel for an object.
-
-    Parameters
-    ----------
-    obj_vign : numpy.ndarray
-        Vignet of the object.
-    psf_vign : numpy.ndarray
-        Vignet of the gaussian model of the PSF.
-    model_vign : numpy.ndarray
-        Vignet of the galaxy model.
-    weight_vign : numpy.ndarray
-        Vignet of the weight at the object position.
-
-    Returns
-    -------
-    sm : float
-        Spread model value.
-    sm_err : float
-        Spread model error value.
-
-    """
-
-    m = (obj_vign > -1e29) & (weight_vign > 0.)
-    w = m.astype(float)
-
-    t_v = model_vign.ravel()
-    g_v = obj_vign.ravel()
-    psf_v = psf_vign.ravel()
-
-    noise_v = (1. / weight_vign).ravel()
-    noise_v[np.isinf(noise_v)] = 0.
-    w_v = w.ravel()
-
-    tg = np.sum(t_v*w_v*g_v)
-    pg = np.sum(psf_v*w_v*g_v)
-    tp = np.sum(t_v*w_v*psf_v)
-    pp = np.sum(psf_v*w_v*psf_v)
-
-    tnt = np.sum(t_v*noise_v*t_v*w_v)
-    pnp = np.sum(psf_v*noise_v*psf_v*w_v)
-    tnp = np.sum(t_v*noise_v*psf_v*w_v)
-    err = tnt * pg**2. + pnp * tg**2. - 2. * tnp * pg * tg
-
-    if pg > 0.:
-        sm = (tg / pg) - (tp / pp)
-    else:
-        sm = 1.
-
-    if (pg > 0.) & (err > 0.):
-        sm_err = np.sqrt(err) / pg**2.
-    else:
-        sm_err = 1.
-
-    return sm, sm_err
-
-
-def get_model(sigma, flux, img_shape, pixel_scale=0.186):
-    """ Get model
-
-    This method compute :
-     - an exponential galaxy model with scale radius = 1/16 FWHM
-     - a gaussian model for the PSF
-
-    Parameters
-    ----------
-    sigma : float
-        Sigma of the PSF (in pixel units).
-    flux : float
-        Flux of the galaxy for the model.
-    img_shape : list
-        Size of the output vignet [xsize, ysize].
-    pixel_scale : float
-        Pixel scale to use for the model (in arcsec/pixel)
-
-    Returns
-    -------
-    gal_vign : numpy.ndarray
-        Vignet of the galaxy model.
-    psf_vign : numpy.ndarray
-        Vignet of the PSF model.
-
-    """
 
-    gal_obj = galsim.Exponential(scale_radius=1./16.*sigma*2.355*pixel_scale,
-                                 flux=flux)
 
-    psf_obj = galsim.Gaussian(sigma=sigma*pixel_scale)
-
-    gal_obj = galsim.Convolve(gal_obj, psf_obj)
-
-    gal_vign = gal_obj.drawImage(nx=img_shape[0], ny=img_shape[1],
-                                 scale=pixel_scale).array
-    psf_vign = psf_obj.drawImage(nx=img_shape[0], ny=img_shape[1],
-                                 scale=pixel_scale).array
-
-    return gal_vign, psf_vign
-
-
-def save_results(sex_cat_path, output_path, sm, sm_err, mag, number,
-                 mode='new'):
-    """ Save results
-
-    Save the spread model results.
-
-    Parameters
-    ----------
-    sex_cat_path : str
-        Path of the original SExtractor catalog.
-    output_path : str
-        Path of the output catalog.
-    sm : numpy.ndarray
-        Value of the spread model for all objects.
-    sm_err : numpy.ndarray
-        Value of the spread model error for all objects.
-    mag : numpy.ndarray
-        Magnitude of all objects (only for new catalog).
-    number : numpy.ndarray
-        Id of all objects (only for new catalog).
-    mode : str
-        Must be in ['new', 'add'].
-        'new' will create a new catalog with : [number, mag, sm, sm_err]
-        'add' will output a copy of the input SExtractor with the column sm
-        and sm_err.
-
-    """
-
-    if mode == 'new':
-        new_cat = io.FITSCatalog(output_path, SEx_catalog=True,
-                                 open_mode=io.BaseCatalog.OpenMode.ReadWrite)
-        dict_data = {'NUMBER': number,
-                     'MAG': mag,
-                     'SPREAD_MODEL': sm,
-                     'SPREADERR_MODEL': sm_err}
-        new_cat.save_as_fits(data=dict_data, sex_cat_path=sex_cat_path)
-    elif mode == 'add':
-        ori_cat = io.FITSCatalog(sex_cat_path, SEx_catalog=True)
-        ori_cat.open()
-        new_cat = io.FITSCatalog(output_path, SEx_catalog=True,
-                                 open_mode=io.BaseCatalog.OpenMode.ReadWrite)
-        ori_cat.add_col('SPREAD_MODEL', sm, new_cat=True, new_cat_inst=new_cat)
-        ori_cat.close()
-        new_cat.open()
-        new_cat.add_col('SPREADERR_MODEL', sm_err)
-        new_cat.close()
-    else:
-        ValueError('Mode must be in [new, add].')
-
-
-@module_runner(input_module=['sextractor_runner', 'psfex_interp_runner_me',
-                             'vignetmaker_runner'], version='1.0',
-               file_pattern=['sexcat', 'galaxy_psf', 'weight_vign'],
-               file_ext=['.fits', '.sqlite', '.fits'],
-               depends=['numpy', 'galsim'])
+@module_runner(
+    version='1.0',
+    input_module=['sextractor_runner', 'psfex_interp_runner_me',
+                  'vignetmaker_runner'],
+    file_pattern=['sexcat', 'galaxy_psf', 'weight_vign'],
+    file_ext=['.fits', '.sqlite', '.fits'],
+    depends=['numpy', 'galsim'],
+    run_method='parallel'
+)
 def spread_model_runner(input_file_list, run_dirs, file_number_string,
                         config, w_log):
 
@@ -189,70 +39,12 @@ def spread_model_runner(input_file_list, run_dirs, file_number_string,
     pixel_scale = config.getfloat('SPREAD_MODEL_RUNNER', 'PIXEL_SCALE')
     output_mode = config.get('SPREAD_MODEL_RUNNER', 'OUTPUT_MODE')
 
-    file_name = suffix + 'sexcat_sm' + file_number_string + '.fits'
-    output_path = run_dirs['output'] + '/' + file_name
-
-    # Get data
-    sex_cat = io.FITSCatalog(sex_cat_path, SEx_catalog=True)
-    sex_cat.open()
-    obj_id = np.copy(sex_cat.get_data()['NUMBER'])
-    obj_vign = np.copy(sex_cat.get_data()['VIGNET'])
-    obj_mag = None
-    if output_mode == 'new':
-        obj_mag = np.copy(sex_cat.get_data()['MAG_AUTO'])
-    sex_cat.close()
-
-    # psf_cat = np.load(psf_cat_path, allow_pickle=True).item()
-    psf_cat = SqliteDict(psf_cat_path)
-
-    weight_cat = io.FITSCatalog(weight_cat_path, SEx_catalog=True)
-    weight_cat.open()
-    weigh_vign = weight_cat.get_data()['VIGNET']
-    weight_cat.close()
-
-    # Get spread model
-    skip_obj = False
-    spread_model_final = []
-    spread_model_err_final = []
-    for i, id_tmp in enumerate(obj_id):
-        sigma_list = []
-
-        if psf_cat[str(id_tmp)] == 'empty':
-            spread_model_final.append(-1)
-            spread_model_err_final.append(1)
-            continue
-
-        psf_expccd_name = list(psf_cat[str(id_tmp)].keys())
-
-        for expccd_name_tmp in psf_expccd_name:
-            sigma_list.append(psf_cat[str(id_tmp)][expccd_name_tmp]['SHAPES']['SIGMA_PSF_HSM'])
-
-        obj_vign_tmp = obj_vign[i]
-        obj_flux_tmp = 1.
-        obj_sigma_tmp = np.mean(sigma_list)
-        obj_weight_tmp = weigh_vign[i]
-        obj_model_tmp, obj_psf_tmp = get_model(obj_sigma_tmp,
-                                               obj_flux_tmp,
-                                               obj_vign_tmp.shape,
-                                               pixel_scale)
-
-        obj_sm, obj_sm_err = get_sm(obj_vign_tmp,
-                                    obj_psf_tmp,
-                                    obj_model_tmp,
-                                    obj_weight_tmp)
-
-        spread_model_final.append(obj_sm)
-        spread_model_err_final.append(obj_sm_err)
-
-    spread_model_final = np.array(spread_model_final)
-    spread_model_err_final = np.array(spread_model_err_final)
+    file_name = f'{suffix}sexcat_sm{file_number_string}.fits'
+    output_path = f'{run_dirs["output"]}/{file_name}'
 
-    psf_cat.close()
+    inst = sm.SpreadModel(sex_cat_path, psf_cat_path, weight_cat_path,
+                          output_path, w_log, pixel_scale, output_mode)
 
-    save_results(sex_cat_path, output_path, spread_model_final,
-                 spread_model_err_final,
-                 obj_mag,
-                 obj_id,
-                 output_mode)
+    inst.process()
 
     return None, None
diff --git a/shapepipe/utilities/__init__.py b/shapepipe/utilities/__init__.py
index 828282ecd..0c25bb8b1 100644
--- a/shapepipe/utilities/__init__.py
+++ b/shapepipe/utilities/__init__.py
@@ -9,4 +9,4 @@
 
 """
 
-__all__ = ['canfar', 'file_system', 'cfis']
+__all__ = ['canfar', 'file_system', 'cfis', 'galaxy']
diff --git a/shapepipe/utilities/galaxy.py b/shapepipe/utilities/galaxy.py
new file mode 100644
index 000000000..43ffd9ee5
--- /dev/null
+++ b/shapepipe/utilities/galaxy.py
@@ -0,0 +1,36 @@
+# -*- coding: utf-8 -*-
+
+"""GALAXY TOOLS
+
+This module defines methods to deal with galaxy images.
+
+:Author: Martin Kilbinger <martin.kilbinger@cea.fr>
+
+:Date: 06/2021
+
+:Package: ShapePipe
+
+"""
+
+
+def sigma_to_fwhm(sigma, pixel_size=1):
+    """sigma to fwhm
+
+    Transform from size sigma to FWHM.
+    The conversion factor corresponds to a 1D Gaussian profile.
+
+    Parameters
+    ----------
+    sigma : (array of) float
+        input size(s)
+    pixel_size : float, optional, default=1
+        pixel size in arcsec, set to 1 if no scaling
+        required
+
+    Returns
+    -------
+    fwhm : (array of) float
+        output fwhm(s)
+    """
+
+    return sigma * 2.355 * pixel_size

From c30caddd4542c420c7cd2f4e55c70c568d66c64e Mon Sep 17 00:00:00 2001
From: martinkilbinger <martin.kilbinger@cea.fr>
Date: Thu, 17 Jun 2021 10:21:33 +0200
Subject: [PATCH 20/30] job_sp script: added vos upload dir to cmd line options

---
 scripts/sh/job_sp.bash | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/scripts/sh/job_sp.bash b/scripts/sh/job_sp.bash
index 00bb78ff7..a7b8d9ebc 100755
--- a/scripts/sh/job_sp.bash
+++ b/scripts/sh/job_sp.bash
@@ -18,6 +18,7 @@ do_env=0
 job=255
 psf='psfex'
 retrieve='vos'
+RESULTS='cosmostat/kilbinger/results'
 nsh_step=3500
 nsh_max=-1
 nsh_jobs=8
@@ -40,6 +41,8 @@ usage="Usage: $(basename "$0") [OPTIONS] TILE_ID_1 [TILE_ID_2 [...]]
     \tPSF model, one in ['psfex'|'mccd'], default='$psf'\n
    -r, --retrieve METHOD\n
    \tmethod to retrieve images, one in ['vos'|'symlink]', default='$retrieve'\n
+   -o, --output_dir\n
+   \toutput (upload) directory on vos:cfis, default='$RESULTS'\n
    --nsh_jobs NJOB\n
    \tnumber of shape measurement parallel jobs, default=$nsh_jobs\n
    --nsh_step NSTEP\n
@@ -129,9 +132,6 @@ fi
 # SExtractor library bug work-around
 export PATH="$PATH:$VM_HOME/bin"
 
-# Results upload subdirectory on vos
-RESULTS=results_$psf
-
 ## Path variables used in shapepipe config files
 
 # Run path and location of input image directories
@@ -242,7 +242,7 @@ function upload() {
       fi
    fi
    tar czf ${base}_${ID}.tgz ${upl[@]}
-   command "$VCP ${base}_${ID}.tgz vos:cfis/cosmostat/kilbinger/$RESULTS" "Upload log tar ball"
+   command "$VCP ${base}_${ID}.tgz vos:cfis/$RESULTS" "Upload log tar ball"
 }
 
 # Upload log files

From 71c94739b1e6d4ff2b035ecbc33669755fd2bde2 Mon Sep 17 00:00:00 2001
From: martinkilbinger <martin.kilbinger@cea.fr>
Date: Thu, 17 Jun 2021 11:20:58 +0200
Subject: [PATCH 21/30] sm revised module, minor bug fixed

---
 .../SpreadModel_package/SpreadModel_script.py | 27 +++++++++----------
 shapepipe/modules/spread_model_runner.py      |  4 ++-
 2 files changed, 15 insertions(+), 16 deletions(-)

diff --git a/shapepipe/modules/SpreadModel_package/SpreadModel_script.py b/shapepipe/modules/SpreadModel_package/SpreadModel_script.py
index 65fc0b256..a695b664d 100644
--- a/shapepipe/modules/SpreadModel_package/SpreadModel_script.py
+++ b/shapepipe/modules/SpreadModel_package/SpreadModel_script.py
@@ -133,8 +133,6 @@ class SpreadModel(object):
         path to weight catalogue
     output_path : string
         output file path of pasted catalog
-    w_log :
-        log file
     pixel_scale : float
         pixel scale in arcsec
     output_mode : str
@@ -145,13 +143,12 @@ class SpreadModel(object):
     """
 
     def __init__(self, sex_cat_path, psf_cat_path, weight_cat_path,
-                 output_path, w_log, pixel_scale, output_mode):
+                 output_path, pixel_scale, output_mode):
 
         self._sex_cat_path = sex_cat_path
         self._psf_cat_path = psf_cat_path
         self._weight_cat_path = weight_cat_path
         self._output_path = output_path
-        self._w_log = w_log
         self._pixel_scale = pixel_scale
         self._output_mode = output_mode
 
@@ -161,18 +158,18 @@ def process(self):
         """
 
         # Get data
-        sex_cat = io.FITSCatalog(sex_cat_path, SEx_catalog=True)
+        sex_cat = io.FITSCatalog(self._sex_cat_path, SEx_catalog=True)
         sex_cat.open()
         obj_id = np.copy(sex_cat.get_data()['NUMBER'])
         obj_vign = np.copy(sex_cat.get_data()['VIGNET'])
         obj_mag = None
-        if output_mode == 'new':
+        if self._output_mode == 'new':
             obj_mag = np.copy(sex_cat.get_data()['MAG_AUTO'])
         sex_cat.close()
 
-        psf_cat = SqliteDict(psf_cat_path)
+        psf_cat = SqliteDict(self._psf_cat_path)
 
-        weight_cat = io.FITSCatalog(weight_cat_path, SEx_catalog=True)
+        weight_cat = io.FITSCatalog(self._weight_cat_path, SEx_catalog=True)
         weight_cat.open()
         weigh_vign = weight_cat.get_data()['VIGNET']
         weight_cat.close()
@@ -201,7 +198,7 @@ def process(self):
             obj_model_tmp, obj_psf_tmp = get_model(obj_sigma_tmp,
                                                    obj_flux_tmp,
                                                    obj_vign_tmp.shape,
-                                                   pixel_scale)
+                                                   self._pixel_scale)
 
             obj_sm, obj_sm_err = get_sm(obj_vign_tmp,
                                         obj_psf_tmp,
@@ -238,18 +235,18 @@ def save_results(self, sm, sm_err, mag, number):
             Id of all objects (only for new catalog).
         """
 
-        if self.output_mode == 'new':
-            new_cat = io.FITSCatalog(self.output_path, SEx_catalog=True,
+        if self._output_mode == 'new':
+            new_cat = io.FITSCatalog(self._output_path, SEx_catalog=True,
                                      open_mode=io.BaseCatalog.OpenMode.ReadWrite)
             dict_data = {'NUMBER': number,
                          'MAG': mag,
                          'SPREAD_MODEL': sm,
                          'SPREADERR_MODEL': sm_err}
-            new_cat.save_as_fits(data=dict_data, sex_cat_path=self.sex_cat_path)
-        elif self.output_mode == 'add':
-            ori_cat = io.FITSCatalog(self.sex_cat_path, SEx_catalog=True)
+            new_cat.save_as_fits(data=dict_data, sex_cat_path=self._sex_cat_path)
+        elif self._output_mode == 'add':
+            ori_cat = io.FITSCatalog(self._sex_cat_path, SEx_catalog=True)
             ori_cat.open()
-            new_cat = io.FITSCatalog(self.output_path, SEx_catalog=True,
+            new_cat = io.FITSCatalog(self._output_path, SEx_catalog=True,
                                      open_mode=io.BaseCatalog.OpenMode.ReadWrite)
             ori_cat.add_col('SPREAD_MODEL', sm, new_cat=True, new_cat_inst=new_cat)
             ori_cat.close()
diff --git a/shapepipe/modules/spread_model_runner.py b/shapepipe/modules/spread_model_runner.py
index 98a63fab7..ecaa10401 100644
--- a/shapepipe/modules/spread_model_runner.py
+++ b/shapepipe/modules/spread_model_runner.py
@@ -11,6 +11,8 @@
 """
 
 from shapepipe.modules.module_decorator import module_runner
+from shapepipe.modules.SpreadModel_package import SpreadModel_script as sm
+
 
 
 @module_runner(
@@ -43,7 +45,7 @@ def spread_model_runner(input_file_list, run_dirs, file_number_string,
     output_path = f'{run_dirs["output"]}/{file_name}'
 
     inst = sm.SpreadModel(sex_cat_path, psf_cat_path, weight_cat_path,
-                          output_path, w_log, pixel_scale, output_mode)
+                          output_path, pixel_scale, output_mode)
 
     inst.process()
 

From ed0d12c2d3cbf666178dafd5d1b734b8e19ca6b0 Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Tue, 29 Jun 2021 13:38:57 +0200
Subject: [PATCH 22/30] FWHM calculation: constant higher ornst higher prec

---
 shapepipe/utilities/galaxy.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/shapepipe/utilities/galaxy.py b/shapepipe/utilities/galaxy.py
index 43ffd9ee5..212cd7f4c 100644
--- a/shapepipe/utilities/galaxy.py
+++ b/shapepipe/utilities/galaxy.py
@@ -33,4 +33,7 @@ def sigma_to_fwhm(sigma, pixel_size=1):
         output fwhm(s)
     """
 
-    return sigma * 2.355 * pixel_size
+    # cst = 2 * sqrt(2 * ln(2))
+    cst = 2.35482004503
+
+    return sigma * cst * pixel_size

From 192b69ca32a9c45d3983a64f3cf8261a646d21cd Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Wed, 21 Jul 2021 18:01:38 +0200
Subject: [PATCH 23/30] sm develop new syntax

---
 .../SpreadModel_package/SpreadModel_script.py | 257 ------------------
 shapepipe/modules/spread_model_runner.py      |  41 ++-
 2 files changed, 30 insertions(+), 268 deletions(-)
 delete mode 100644 shapepipe/modules/SpreadModel_package/SpreadModel_script.py

diff --git a/shapepipe/modules/SpreadModel_package/SpreadModel_script.py b/shapepipe/modules/SpreadModel_package/SpreadModel_script.py
deleted file mode 100644
index a695b664d..000000000
--- a/shapepipe/modules/SpreadModel_package/SpreadModel_script.py
+++ /dev/null
@@ -1,257 +0,0 @@
-# -*- coding: utf-8 -*-
-
-"""SPREAD MODEL SCRIPT
-
-Class to compute the spread model, criterion to select galaxies
-
-:Author: Axel Guinot
-
-:Date: 2019, 2020
-
-:Package: ShapePipe
-
-"""
-
-
-import numpy as np
-import galsim
-from sqlitedict import SqliteDict
-
-from shapepipe.pipeline import file_io as io
-from shapepipe.utilities import galaxy
-
-
-def get_sm(obj_vign, psf_vign, model_vign, weight_vign):
-    """ Get Spread model
-
-    This method compute the spread moel for an object.
-
-    Parameters
-    ----------
-    obj_vign : numpy.ndarray
-        Vignet of the object.
-    psf_vign : numpy.ndarray
-        Vignet of the gaussian model of the PSF.
-    model_vign : numpy.ndarray
-        Vignet of the galaxy model.
-    weight_vign : numpy.ndarray
-        Vignet of the weight at the object position.
-
-    Returns
-    -------
-    sm : float
-        Spread model value.
-    sm_err : float
-        Spread model error value.
-    """
-
-    m = (obj_vign > -1e29) & (weight_vign > 0)
-    w = m.astype(float)
-
-    t_v = model_vign.ravel()
-    g_v = obj_vign.ravel()
-    psf_v = psf_vign.ravel()
-
-    noise_v = (1 / weight_vign).ravel()
-    noise_v[np.isinf(noise_v)] = 0
-    w_v = w.ravel()
-
-    tg = np.sum(t_v*w_v*g_v)
-    pg = np.sum(psf_v*w_v*g_v)
-    tp = np.sum(t_v*w_v*psf_v)
-    pp = np.sum(psf_v*w_v*psf_v)
-
-    tnt = np.sum(t_v*noise_v*t_v*w_v)
-    pnp = np.sum(psf_v*noise_v*psf_v*w_v)
-    tnp = np.sum(t_v*noise_v*psf_v*w_v)
-    err = tnt * pg**2 + pnp * tg**2 - 2*tnp * pg * tg
-
-    if pg > 0:
-        sm = (tg / pg) - (tp / pp)
-    else:
-        sm = 1
-
-    if (pg > 0) & (err > 0):
-        sm_err = np.sqrt(err) / pg**2
-    else:
-        sm_err = 1
-
-    return sm, sm_err
-
-
-def get_model(sigma, flux, img_shape, pixel_scale=0.186):
-    """ Get model
-
-    This method computes :
-     - an exponential galaxy model with scale radius = 1/16 FWHM
-     - a gaussian model for the PSF
-
-    Parameters
-    ----------
-    sigma : float
-        Sigma of the PSF (in pixel units).
-    flux : float
-        Flux of the galaxy for the model.
-    img_shape : list
-        Size of the output vignet [xsize, ysize].
-    pixel_scale : float, optional, default=0.186
-        Pixel scale to use for the model (in arcsec)
-
-    Returns
-    -------
-    gal_vign : numpy.ndarray
-        Vignet of the galaxy model.
-    psf_vign : numpy.ndarray
-        Vignet of the PSF model.
-    """
-
-    scale_radius = 1/16 * galaxy.sigma_to_fhwm(sigma, pixel_scale=pixel_scale)
-    gal_obj = galsim.Exponential(scale_radius=scale_radius, flux=flux)
-
-    psf_obj = galsim.Gaussian(sigma=sigma*pixel_scale)
-
-    gal_obj = galsim.Convolve(gal_obj, psf_obj)
-
-    gal_vign = gal_obj.drawImage(nx=img_shape[0], ny=img_shape[1],
-                                 scale=pixel_scale).array
-    psf_vign = psf_obj.drawImage(nx=img_shape[0], ny=img_shape[1],
-                                 scale=pixel_scale).array
-
-    return gal_vign, psf_vign
-
-
-class SpreadModel(object):
-    """SpreadModel initialisation.
-
-    Parameters
-    ----------
-    sex_cat_path : string
-        path to SExtractor catalogue
-    psf_cat_path : string
-        path to PSF catalogue
-    weight_cat_path : string
-        path to weight catalogue
-    output_path : string
-        output file path of pasted catalog
-    pixel_scale : float
-        pixel scale in arcsec
-    output_mode : str
-        must be in ['new', 'add'].
-         'new' will create a new catalog with : [number, mag, sm, sm_err]
-         'add' will output a copy of the input SExtractor with the column sm
-          and sm_err.
-    """
-
-    def __init__(self, sex_cat_path, psf_cat_path, weight_cat_path,
-                 output_path, pixel_scale, output_mode):
-
-        self._sex_cat_path = sex_cat_path
-        self._psf_cat_path = psf_cat_path
-        self._weight_cat_path = weight_cat_path
-        self._output_path = output_path
-        self._pixel_scale = pixel_scale
-        self._output_mode = output_mode
-
-    def process(self):
-        """Process
-        Process the spread model computation
-        """
-
-        # Get data
-        sex_cat = io.FITSCatalog(self._sex_cat_path, SEx_catalog=True)
-        sex_cat.open()
-        obj_id = np.copy(sex_cat.get_data()['NUMBER'])
-        obj_vign = np.copy(sex_cat.get_data()['VIGNET'])
-        obj_mag = None
-        if self._output_mode == 'new':
-            obj_mag = np.copy(sex_cat.get_data()['MAG_AUTO'])
-        sex_cat.close()
-
-        psf_cat = SqliteDict(self._psf_cat_path)
-
-        weight_cat = io.FITSCatalog(self._weight_cat_path, SEx_catalog=True)
-        weight_cat.open()
-        weigh_vign = weight_cat.get_data()['VIGNET']
-        weight_cat.close()
-
-        # Get spread model
-        skip_obj = False
-        spread_model_final = []
-        spread_model_err_final = []
-        for i, id_tmp in enumerate(obj_id):
-            sigma_list = []
-
-            if psf_cat[str(id_tmp)] == 'empty':
-                spread_model_final.append(-1)
-                spread_model_err_final.append(1)
-                continue
-
-            psf_expccd_name = list(psf_cat[str(id_tmp)].keys())
-
-            for expccd_name_tmp in psf_expccd_name:
-                sigma_list.append(psf_cat[str(id_tmp)][expccd_name_tmp]['SHAPES']['SIGMA_PSF_HSM'])
-
-            obj_vign_tmp = obj_vign[i]
-            obj_flux_tmp = 1.
-            obj_sigma_tmp = np.mean(sigma_list)
-            obj_weight_tmp = weigh_vign[i]
-            obj_model_tmp, obj_psf_tmp = get_model(obj_sigma_tmp,
-                                                   obj_flux_tmp,
-                                                   obj_vign_tmp.shape,
-                                                   self._pixel_scale)
-
-            obj_sm, obj_sm_err = get_sm(obj_vign_tmp,
-                                        obj_psf_tmp,
-                                        obj_model_tmp,
-                                        obj_weight_tmp)
-
-            spread_model_final.append(obj_sm)
-            spread_model_err_final.append(obj_sm_err)
-
-        spread_model_final = np.array(spread_model_final)
-        spread_model_err_final = np.array(spread_model_err_final)
-
-        psf_cat.close()
-
-        self.save_results(spread_model_final,
-                          spread_model_err_final,
-                          obj_mag,
-                          obj_id)
-
-    def save_results(self, sm, sm_err, mag, number):
-        """Save Results
-
-        Save output catalogue with spread model and errors.
-
-        Parameters
-        ----------
-        sm : numpy.ndarray
-            Value of the spread model for all objects.
-        sm_err : numpy.ndarray
-            Value of the spread model error for all objects.
-        mag : numpy.ndarray
-            Magnitude of all objects (only for new catalog).
-        number : numpy.ndarray
-            Id of all objects (only for new catalog).
-        """
-
-        if self._output_mode == 'new':
-            new_cat = io.FITSCatalog(self._output_path, SEx_catalog=True,
-                                     open_mode=io.BaseCatalog.OpenMode.ReadWrite)
-            dict_data = {'NUMBER': number,
-                         'MAG': mag,
-                         'SPREAD_MODEL': sm,
-                         'SPREADERR_MODEL': sm_err}
-            new_cat.save_as_fits(data=dict_data, sex_cat_path=self._sex_cat_path)
-        elif self._output_mode == 'add':
-            ori_cat = io.FITSCatalog(self._sex_cat_path, SEx_catalog=True)
-            ori_cat.open()
-            new_cat = io.FITSCatalog(self._output_path, SEx_catalog=True,
-                                     open_mode=io.BaseCatalog.OpenMode.ReadWrite)
-            ori_cat.add_col('SPREAD_MODEL', sm, new_cat=True, new_cat_inst=new_cat)
-            ori_cat.close()
-            new_cat.open()
-            new_cat.add_col('SPREADERR_MODEL', sm_err)
-            new_cat.close()
-        else:
-            ValueError('Mode must be in [new, add].')
diff --git a/shapepipe/modules/spread_model_runner.py b/shapepipe/modules/spread_model_runner.py
index ecaa10401..29f84103e 100644
--- a/shapepipe/modules/spread_model_runner.py
+++ b/shapepipe/modules/spread_model_runner.py
@@ -2,7 +2,7 @@
 
 """SPREAD MODEL RUNNER
 
-This module computes the spread model.
+Module runner for ``spread_model``
 
 :Author: Axel Guinot
 
@@ -11,24 +11,33 @@
 """
 
 from shapepipe.modules.module_decorator import module_runner
-from shapepipe.modules.SpreadModel_package import SpreadModel_script as sm
-
+from shapepipe.modules.spread_mdel_package import spread_model as sm
 
 
 @module_runner(
     version='1.0',
-    input_module=['sextractor_runner', 'psfex_interp_runner_me',
-                  'vignetmaker_runner'],
+    input_module=[
+        'sextractor_runner',
+        'psfex_interp_runner_me',
+        'vignetmaker_runner'
+    ],
     file_pattern=['sexcat', 'galaxy_psf', 'weight_vign'],
     file_ext=['.fits', '.sqlite', '.fits'],
     depends=['numpy', 'galsim'],
     run_method='parallel'
 )
-def spread_model_runner(input_file_list, run_dirs, file_number_string,
-                        config, w_log):
-
+def spread_model_runner(
+    input_file_list,
+    run_dirs,
+    file_number_string,
+    config,
+    w_log
+):
+
+    # Get input files
     sex_cat_path, psf_cat_path, weight_cat_path = input_file_list
 
+    # Get file suffix (optional)
     if config.has_option('SPREAD_MODEL_RUNNER', 'SUFFIX'):
         suffix = config.get('SPREAD_MODEL_RUNNER', 'SUFFIX')
         if (suffix.lower() != 'none') & (suffix != ''):
@@ -38,15 +47,25 @@ def spread_model_runner(input_file_list, run_dirs, file_number_string,
     else:
         suffix = ''
 
+    # Get pixel scale and output mode
     pixel_scale = config.getfloat('SPREAD_MODEL_RUNNER', 'PIXEL_SCALE')
     output_mode = config.get('SPREAD_MODEL_RUNNER', 'OUTPUT_MODE')
 
+    # Set output file path
     file_name = f'{suffix}sexcat_sm{file_number_string}.fits'
     output_path = f'{run_dirs["output"]}/{file_name}'
 
-    inst = sm.SpreadModel(sex_cat_path, psf_cat_path, weight_cat_path,
-                          output_path, pixel_scale, output_mode)
-
+    # Create spread model class instance
+    inst = sm.SpreadModel(
+        sex_cat_path,
+        psf_cat_path,
+        weight_cat_path,
+        output_path,
+        pixel_scale,
+        output_mode
+    )
+
+    # Process spread model computation
     inst.process()
 
     return None, None

From a479594936762858969b0fc2b2fe6b6d0de755d1 Mon Sep 17 00:00:00 2001
From: martinkilbinger <martin.kilbinger@cea.fr>
Date: Thu, 22 Jul 2021 17:41:37 +0200
Subject: [PATCH 24/30] typo fixed

---
 shapepipe/modules/spread_model_runner.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/shapepipe/modules/spread_model_runner.py b/shapepipe/modules/spread_model_runner.py
index 29f84103e..2630ac5bd 100644
--- a/shapepipe/modules/spread_model_runner.py
+++ b/shapepipe/modules/spread_model_runner.py
@@ -11,7 +11,7 @@
 """
 
 from shapepipe.modules.module_decorator import module_runner
-from shapepipe.modules.spread_mdel_package import spread_model as sm
+from shapepipe.modules.spread_model_package import spread_model as sm
 
 
 @module_runner(

From e1e265f824aef62a7069c515ece201cc04081a29 Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Thu, 22 Jul 2021 17:44:32 +0200
Subject: [PATCH 25/30]  added sm script

---
 .../spread_model_package/spread_model.py      | 309 ++++++++++++++++++
 1 file changed, 309 insertions(+)
 create mode 100644 shapepipe/modules/spread_model_package/spread_model.py

diff --git a/shapepipe/modules/spread_model_package/spread_model.py b/shapepipe/modules/spread_model_package/spread_model.py
new file mode 100644
index 000000000..877cd9083
--- /dev/null
+++ b/shapepipe/modules/spread_model_package/spread_model.py
@@ -0,0 +1,309 @@
+# -*- coding: utf-8 -*-
+
+"""SPREAD MODEL SCRIPT
+
+Class to compute the spread model, criterion to select galaxies
+
+:Author: Axel Guinot
+
+:Date: 2019, 2020
+
+:Package: ShapePipe
+
+"""
+
+
+import numpy as np
+import galsim
+from sqlitedict import SqliteDict
+
+from shapepipe.pipeline import file_io as io
+from shapepipe.utilities import galaxy
+
+
+def get_sm(obj_vign, psf_vign, model_vign, weight_vign):
+    """ Get Spread model
+
+    This method compute the spread moel for an object.
+
+    Parameters
+    ----------
+    obj_vign : numpy.ndarray
+        Vignet of the object.
+    psf_vign : numpy.ndarray
+        Vignet of the gaussian model of the PSF.
+    model_vign : numpy.ndarray
+        Vignet of the galaxy model.
+    weight_vign : numpy.ndarray
+        Vignet of the weight at the object position.
+
+    Returns
+    -------
+    sm : float
+        Spread model value.
+    sm_err : float
+        Spread model error value.
+    """
+
+    # Mask invalid pixels
+    m = (obj_vign > -1e29) & (weight_vign > 0)
+    w = m.astype(float)
+
+    # Set noise as inverse weight
+    noise_v = (1 / weight_vign).ravel()
+
+    # Remove infinite noise pixels
+    noise_v[np.isinf(noise_v)] = 0
+
+    # Transform 2D vignets to 1D vectors
+    t_v = model_vign.ravel()
+    g_v = obj_vign.ravel()
+    psf_v = psf_vign.ravel()
+    w_v = w.ravel()
+
+    # Compute scalar products used in spread model
+    tg = np.sum(t_v*w_v*g_v)
+    pg = np.sum(psf_v*w_v*g_v)
+    tp = np.sum(t_v*w_v*psf_v)
+    pp = np.sum(psf_v*w_v*psf_v)
+
+    tnt = np.sum(t_v*noise_v*t_v*w_v)
+    pnp = np.sum(psf_v*noise_v*psf_v*w_v)
+    tnp = np.sum(t_v*noise_v*psf_v*w_v)
+    err = tnt * pg**2 + pnp * tg**2 - 2*tnp * pg * tg
+
+    # Compute spread model
+    if pg > 0:
+        sm = (tg / pg) - (tp / pp)
+    else:
+        sm = 1
+
+    if (pg > 0) & (err > 0):
+        sm_err = np.sqrt(err) / pg**2
+    else:
+        sm_err = 1
+
+    return sm, sm_err
+
+
+def get_model(sigma, flux, img_shape, pixel_scale=0.186):
+    """ Get model
+
+    This method computes
+     - an exponential galaxy model with scale radius = 1/16 FWHM
+     - a gaussian model for the PSF
+
+    Parameters
+    ----------
+    sigma : float
+        Sigma of the PSF (in pixel units).
+    flux : float
+        Flux of the galaxy for the model.
+    img_shape : list
+        Size of the output vignet [xsize, ysize].
+    pixel_scale : float, optional, default=0.186
+        Pixel scale to use for the model (in arcsec)
+
+    Returns
+    -------
+    gal_vign : numpy.ndarray
+        Vignet of the galaxy model.
+    psf_vign : numpy.ndarray
+        Vignet of the PSF model.
+    """
+
+    # Get scale radius
+    scale_radius = 1/16 * galaxy.sigma_to_fhwm(sigma, pixel_scale=pixel_scale)
+
+    # Get galaxy model
+    gal_obj = galsim.Exponential(scale_radius=scale_radius, flux=flux)
+
+    # Get PSF
+    psf_obj = galsim.Gaussian(sigma=sigma*pixel_scale)
+
+    # Convolve both
+    gal_obj = galsim.Convolve(gal_obj, psf_obj)
+
+    # Draw galaxy and PSF on vignets
+    gal_vign = gal_obj.drawImage(
+        nx=img_shape[0],
+        ny=img_shape[1],
+        scale=pixel_scale
+    ).array
+
+    psf_vign = psf_obj.drawImage(
+        nx=img_shape[0],
+        ny=img_shape[1],
+        scale=pixel_scale
+    ).array
+
+    return gal_vign, psf_vign
+
+
+class SpreadModel(object):
+    """SpreadModel class.
+
+    Parameters
+    ----------
+    sex_cat_path : string
+        path to SExtractor catalogue
+    psf_cat_path : string
+        path to PSF catalogue
+    weight_cat_path : string
+        path to weight catalogue
+    output_path : string
+        output file path of pasted catalog
+    pixel_scale : float
+        pixel scale in arcsec
+    output_mode : str
+        must be in ['new', 'add'].
+         'new' will create a new catalog with : [number, mag, sm, sm_err]
+         'add' will output a copy of the input SExtractor with the column sm
+          and sm_err.
+    """
+
+    def __init__(
+        self,
+        sex_cat_path,
+        psf_cat_path,
+        weight_cat_path,
+        output_path,
+        pixel_scale,
+        output_mode
+    ):
+
+        self._sex_cat_path = sex_cat_path
+        self._psf_cat_path = psf_cat_path
+        self._weight_cat_path = weight_cat_path
+        self._output_path = output_path
+        self._pixel_scale = pixel_scale
+        self._output_mode = output_mode
+
+    def process(self):
+        """Process
+        Process the spread model computation
+        """
+
+        # Get data
+        sex_cat = io.FITSCatalog(self._sex_cat_path, SEx_catalog=True)
+        sex_cat.open()
+        obj_id = np.copy(sex_cat.get_data()['NUMBER'])
+        obj_vign = np.copy(sex_cat.get_data()['VIGNET'])
+        obj_mag = None
+        if self._output_mode == 'new':
+            obj_mag = np.copy(sex_cat.get_data()['MAG_AUTO'])
+        sex_cat.close()
+
+        psf_cat = SqliteDict(self._psf_cat_path)
+
+        weight_cat = io.FITSCatalog(self._weight_cat_path, SEx_catalog=True)
+        weight_cat.open()
+        weigh_vign = weight_cat.get_data()['VIGNET']
+        weight_cat.close()
+
+        # Get spread model
+        skip_obj = False
+        spread_model_final = []
+        spread_model_err_final = []
+        for i, id_tmp in enumerate(obj_id):
+            sigma_list = []
+
+            if psf_cat[str(id_tmp)] == 'empty':
+                spread_model_final.append(-1)
+                spread_model_err_final.append(1)
+                continue
+
+            psf_expccd_name = list(psf_cat[str(id_tmp)].keys())
+
+            for expccd_name_tmp in psf_expccd_name:
+                psf_cat_id_ccd = psf_cat[str(id_tmp)][expccd_name_tmp]
+                sigma_list.append(
+                    psf_cat_id_ccd['SHAPES']['SIGMA_PSF_HSM']
+                )
+
+            obj_vign_tmp = obj_vign[i]
+            obj_flux_tmp = 1.
+            obj_sigma_tmp = np.mean(sigma_list)
+            obj_weight_tmp = weigh_vign[i]
+            obj_model_tmp, obj_psf_tmp = get_model(
+                obj_sigma_tmp,
+                obj_flux_tmp,
+                obj_vign_tmp.shape,
+                self._pixel_scale
+            )
+
+            obj_sm, obj_sm_err = get_sm(
+                obj_vign_tmp,
+                obj_psf_tmp,
+                obj_model_tmp,
+                obj_weight_tmp
+            )
+
+            spread_model_final.append(obj_sm)
+            spread_model_err_final.append(obj_sm_err)
+
+        spread_model_final = np.array(spread_model_final)
+        spread_model_err_final = np.array(spread_model_err_final)
+
+        psf_cat.close()
+
+        self.save_results(
+            spread_model_final,
+            spread_model_err_final,
+            obj_mag,
+            obj_id
+        )
+
+    def save_results(self, sm, sm_err, mag, number):
+        """Save Results
+
+        Save output catalogue with spread model and errors.
+
+        Parameters
+        ----------
+        sm : numpy.ndarray
+            Value of the spread model for all objects.
+        sm_err : numpy.ndarray
+            Value of the spread model error for all objects.
+        mag : numpy.ndarray
+            Magnitude of all objects (only for new catalog).
+        number : numpy.ndarray
+            Id of all objects (only for new catalog).
+        """
+
+        if self._output_mode == 'new':
+            new_cat = io.FITSCatalog(
+                self._output_path,
+                SEx_catalog=True,
+                open_mode=io.BaseCatalog.OpenMode.ReadWrite
+            )
+            dict_data = {
+                'NUMBER': number,
+                'MAG': mag,
+                'SPREAD_MODEL': sm,
+                'SPREADERR_MODEL': sm_err
+            }
+            new_cat.save_as_fits(
+                data=dict_data,
+                sex_cat_path=self._sex_cat_path
+            )
+        elif self._output_mode == 'add':
+            ori_cat = io.FITSCatalog(self._sex_cat_path, SEx_catalog=True)
+            ori_cat.open()
+            new_cat = io.FITSCatalog(
+                self._output_path,
+                SEx_catalog=True,
+                open_mode=io.BaseCatalog.OpenMode.ReadWrite
+            )
+            ori_cat.add_col(
+                'SPREAD_MODEL',
+                sm,
+                new_cat=True,
+                new_cat_inst=new_cat
+            )
+            ori_cat.close()
+            new_cat.open()
+            new_cat.add_col('SPREADERR_MODEL', sm_err)
+            new_cat.close()
+        else:
+            raise ValueError('Mode must be in [new, add].')

From bb7777766ae656c8d263a3f13a91a93bb10cc3ed Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Thu, 22 Jul 2021 17:45:26 +0200
Subject: [PATCH 26/30]  added sm script

---
 shapepipe/modules/spread_model_package/__init__.py | 11 +++++++++++
 1 file changed, 11 insertions(+)
 create mode 100644 shapepipe/modules/spread_model_package/__init__.py

diff --git a/shapepipe/modules/spread_model_package/__init__.py b/shapepipe/modules/spread_model_package/__init__.py
new file mode 100644
index 000000000..0d7f2f8d8
--- /dev/null
+++ b/shapepipe/modules/spread_model_package/__init__.py
@@ -0,0 +1,11 @@
+# -*- coding: utf-8 -*-
+
+"""VIGNET MAKER PACKAGE
+
+This package contains the module(s) for ``spread_model``.
+
+:Author: Axel Guinot
+
+"""
+
+__all__ = ['spread_model']

From 2c55174b9fb97d034f23113023cd73762c93fdcc Mon Sep 17 00:00:00 2001
From: martinkilbinger <martin.kilbinger@cea.fr>
Date: Thu, 22 Jul 2021 21:28:21 +0200
Subject: [PATCH 27/30] removed vignetmaker_runner2 from config file, to be
 tested with other PR

---
 example/cfis/config_tile_SxMiViSmVi.ini | 2 +-
 example/cfis/config_tile_SxPiViSmVi.ini | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/example/cfis/config_tile_SxMiViSmVi.ini b/example/cfis/config_tile_SxMiViSmVi.ini
index ddf3c04ac..2664d35c2 100644
--- a/example/cfis/config_tile_SxMiViSmVi.ini
+++ b/example/cfis/config_tile_SxMiViSmVi.ini
@@ -21,7 +21,7 @@ RUN_NAME = run_sp_SxMiViSmVi
 # Module name, single string or comma-separated list of valid module runner names
 MODULE = sextractor_runner, mccd_interp_runner,
          vignetmaker_runner, spread_model_runner,
-         vignetmaker_runner2
+         vignetmaker_runner
 
 # Parallel processing mode, SMP or MPI
 MODE = SMP
diff --git a/example/cfis/config_tile_SxPiViSmVi.ini b/example/cfis/config_tile_SxPiViSmVi.ini
index 3f6f2ccaa..d45428e41 100644
--- a/example/cfis/config_tile_SxPiViSmVi.ini
+++ b/example/cfis/config_tile_SxPiViSmVi.ini
@@ -21,7 +21,7 @@ RUN_NAME = run_sp_SxPsViSmVi
 # Module name, single string or comma-separated list of valid module runner names
 MODULE = sextractor_runner, psfex_interp_runner,
          vignetmaker_runner, spread_model_runner,
-         vignetmaker_runner2
+         vignetmaker_runner
 
 # Parallel processing mode, SMP or MPI
 MODE = SMP

From 194523b6e0c70aae6469e25e12e2977cc669cc75 Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Fri, 30 Jul 2021 11:57:50 +0200
Subject: [PATCH 28/30] PR review changes

---
 scripts/python/merge_final_cat.py             | 35 ++++++++++++-------
 .../spread_model_package/spread_model.py      |  5 +++
 2 files changed, 28 insertions(+), 12 deletions(-)

diff --git a/scripts/python/merge_final_cat.py b/scripts/python/merge_final_cat.py
index 356696159..d58988313 100755
--- a/scripts/python/merge_final_cat.py
+++ b/scripts/python/merge_final_cat.py
@@ -81,14 +81,25 @@ def parse_options(p_def):
     usage  = "%prog [OPTIONS]"
     parser = OptionParser(usage=usage)
 
-    parser.add_option('-i', '--input_path', dest='input_path', type='string',
-                      default=p_def.input_path,
-                      help='input path, default=\'{}\''.format(p_def.input_path))
+    parser.add_option(
+        '-i',
+        '--input_path',
+        dest='input_path',
+        type='string',
+        default=p_def.input_path,
+        help=f'input path, default=\'{p_def.input_path}\''
+    )
     parser.add_option('-p', '--param_path', dest='param_path', type='string',
                       default=None,
                       help='parameter file path, default=None')
 
-    parser.add_option('-v', '--verbose', dest='verbose', action='store_true', help='verbose output')
+    parser.add_option(
+        '-v',
+        '--verbose',
+        dest='verbose',
+        action='store_true',
+        help='verbose output'
+    )
 
     options, args = parser.parse_args()
 
@@ -176,7 +187,7 @@ def read_param_file(path, verbose=False):
 
     if verbose:
         if len(param_list) > 0: 
-            print('Copying {} columns'.format(len(param_list)), end='')
+            print(f'Copying {len(param_list)} columns', end='')
         else:
             print('Copying all columns', end='')
         print(' into final catalogue')
@@ -259,7 +270,7 @@ def main(argv=None):
         lpath.append(os.path.join(path, this_l))
 
     if param.verbose:
-        print('{} final catalog files found'.format(len(lpath)))
+        print(f'{len(lpath)} final catalog files found')
 
     count = 0
 
@@ -271,25 +282,25 @@ def main(argv=None):
         d[key] = d_tmp[key]
     count = count + 1
     if param.verbose:
-        print('File \'{}\' copied ({}/{})'.format(lpath[0], count, len(lpath)))
+        print(f'File \'lpath[0]{}\' copied ({count}/{len(lpath)})')
 
-    for i in lpath[1:]:
-        if ('final_cat' not in i) | ('.npy' in i):
+    for idx in lpath[1:]:
+        if ('final_cat' not in idx) | ('.npy' in idx):
             continue
 
         try:
-            d_tmp = get_data(i, 1, param.param_list)
+            d_tmp = get_data(idx, 1, param.param_list)
 
             dd = np.zeros(d_tmp.shape, dtype=d_tmp.dtype)
             for key in d_tmp.dtype.names:
                 dd[key] = d_tmp[key]
             count = count + 1
             if param.verbose:
-                print('File \'{}\' copied ({}/{})'.format(i, count, len(lpath)))
+                print(f'File \'{idx}\' copied ({count}/{len(lpath)})')
 
             d = np.concatenate((d, dd))
         except:
-            print('Error while copying file \'{}\''.format(i))
+            print(f'Error while copying file \'{idx}\'')
 
     # Save merged catalogue as numpy binary file
     if param.verbose:
diff --git a/shapepipe/modules/spread_model_package/spread_model.py b/shapepipe/modules/spread_model_package/spread_model.py
index 877cd9083..66214d764 100644
--- a/shapepipe/modules/spread_model_package/spread_model.py
+++ b/shapepipe/modules/spread_model_package/spread_model.py
@@ -269,6 +269,11 @@ def save_results(self, sm, sm_err, mag, number):
             Magnitude of all objects (only for new catalog).
         number : numpy.ndarray
             Id of all objects (only for new catalog).
+
+        Raises
+        ------
+        ValueError
+            For incorrect output mod
         """
 
         if self._output_mode == 'new':

From d8336452ba61706297cdbe667dd9628e6ab6cfc7 Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Fri, 30 Jul 2021 11:59:16 +0200
Subject: [PATCH 29/30] PR review change requests

---
 .../modules/spread_model_package/spread_model.py     | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/shapepipe/modules/spread_model_package/spread_model.py b/shapepipe/modules/spread_model_package/spread_model.py
index 66214d764..074d43e5f 100644
--- a/shapepipe/modules/spread_model_package/spread_model.py
+++ b/shapepipe/modules/spread_model_package/spread_model.py
@@ -145,13 +145,13 @@ class SpreadModel(object):
 
     Parameters
     ----------
-    sex_cat_path : string
+    sex_cat_path : str
         path to SExtractor catalogue
-    psf_cat_path : string
+    psf_cat_path : str
         path to PSF catalogue
-    weight_cat_path : string
+    weight_cat_path : str
         path to weight catalogue
-    output_path : string
+    output_path : str
         output file path of pasted catalog
     pixel_scale : float
         pixel scale in arcsec
@@ -205,7 +205,7 @@ def process(self):
         skip_obj = False
         spread_model_final = []
         spread_model_err_final = []
-        for i, id_tmp in enumerate(obj_id):
+        for idx id_tmp in enumerate(obj_id):
             sigma_list = []
 
             if psf_cat[str(id_tmp)] == 'empty':
@@ -224,7 +224,7 @@ def process(self):
             obj_vign_tmp = obj_vign[i]
             obj_flux_tmp = 1.
             obj_sigma_tmp = np.mean(sigma_list)
-            obj_weight_tmp = weigh_vign[i]
+            obj_weight_tmp = weigh_vign[idx]
             obj_model_tmp, obj_psf_tmp = get_model(
                 obj_sigma_tmp,
                 obj_flux_tmp,

From 1078b675651879de76f9cd2f8a0f54b51d712bcd Mon Sep 17 00:00:00 2001
From: Martin Kilbinger <martin.kilbinger@cea.fr>
Date: Fri, 30 Jul 2021 13:12:38 +0200
Subject: [PATCH 30/30] added module_config_sec

---
 shapepipe/modules/spread_model_runner.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/shapepipe/modules/spread_model_runner.py b/shapepipe/modules/spread_model_runner.py
index 2630ac5bd..09858f7df 100644
--- a/shapepipe/modules/spread_model_runner.py
+++ b/shapepipe/modules/spread_model_runner.py
@@ -31,6 +31,7 @@ def spread_model_runner(
     run_dirs,
     file_number_string,
     config,
+    module_config_sec,
     w_log
 ):
 
@@ -38,8 +39,8 @@ def spread_model_runner(
     sex_cat_path, psf_cat_path, weight_cat_path = input_file_list
 
     # Get file suffix (optional)
-    if config.has_option('SPREAD_MODEL_RUNNER', 'SUFFIX'):
-        suffix = config.get('SPREAD_MODEL_RUNNER', 'SUFFIX')
+    if config.has_option(module_config_sec, 'SUFFIX'):
+        suffix = config.get(module_config_sec, 'SUFFIX')
         if (suffix.lower() != 'none') & (suffix != ''):
             suffix = suffix + '_'
         else:
@@ -48,8 +49,8 @@ def spread_model_runner(
         suffix = ''
 
     # Get pixel scale and output mode
-    pixel_scale = config.getfloat('SPREAD_MODEL_RUNNER', 'PIXEL_SCALE')
-    output_mode = config.get('SPREAD_MODEL_RUNNER', 'OUTPUT_MODE')
+    pixel_scale = config.getfloat(module_config_sec, 'PIXEL_SCALE')
+    output_mode = config.get(module_config_sec, 'OUTPUT_MODE')
 
     # Set output file path
     file_name = f'{suffix}sexcat_sm{file_number_string}.fits'