Merge branch 'develop'

.gitignore

@@ -47,6 +47,11 @@ pkgconfig
 # Temp files
 *~
 
+# Index files
+TAGS
+.DS_Store
+Thumbs.db
+
 # Python cache files
 *.pyc
 *__pycache__
@@ -55,6 +60,7 @@ pyamazed/redshift_wrap.cpp
 pyamazed/redshift_wrap.h
 pyamazed/redshift.py
 *.egg*
+fixed-wheel
 
 swig_ext.py
 

.gitlab-ci.yml

@@ -14,10 +14,11 @@ job:
   script:
     - mkdir build
     - cd build
-    - cmake .. -DBUILD_SHARED_LIBS=ON -DCMAKE_BUILD_TYPE=Release
+    - cmake .. -DBUILD_SHARED_LIBS=ON -DCMAKE_BUILD_TYPE=Release -DBUILD_TESTS=ON
     - make -j4
+    - make install
+    - make test
     - cd ..
     - python3.6 -m venv --system-site-packages venv
     - source venv/bin/activate
     - pip3.6 install .
-    - make -C build tests

CHANGES.md

@@ -0,0 +1,16 @@
+# 0.6.0 (2019-04-11)
+
+## New Features
+
+* Added the new option `tplfitauto` to `continuumcomponent` parameter. With this option, the code selects the best method (`linemodel` or `fullmodel`) for redshift determination depending on the signal to noise ratio (SNR) of all template fitting. If the best SNR is greater than 50 then the `fullmodel` is selected. Otherwise the `linemodel` is selected
+* Added the new option `svdlcp2` to `secondpasslcfittingmethod`. Enable second degree polynomial on line support on second pass (correct continuum removal residue)
+* Implemented new fit options for LyAlpha line: `lyafit.asymfitmin`, `lyafit.asymfitmax`, `lyafit.asymfitstep`, `lyafit.widthfitmin`, `lyafit.widthfitmax`, `lyafit.widthfitstep`, `lyafit.deltafitmin`, `lyafit.deltafitmax`, `lyafit.deltafitstep`
+
+## API Changes
+None
+
+## Bug Fixes
+None
+
+## Other Changes and Additions
+None

CMakeLists.txt

@@ -19,12 +19,28 @@ ELSE()
   SET(Boost_USE_STATIC_LIBS ON)
 ENDIF()
 
-FIND_PACKAGE(Boost COMPONENTS filesystem system thread timer chrono program_options unit_test_framework)
-FIND_PACKAGE(Cfitsio)
-FIND_PACKAGE(GSL 2.1)
-FIND_PACKAGE(FFTW)
+LIST(FIND MISSING_THIRDPARTIES "boost" find_boost)
+LIST(FIND MISSING_THIRDPARTIES "cfitsio" find_cfitsio)
+LIST(FIND MISSING_THIRDPARTIES "gsl" find_gsl)
+LIST(FIND MISSING_THIRDPARTIES "fftw" find_fftw)
 
-IF(NOT GSL_FOUND)
+IF(NOT find_boost GREATER -1)
+  FIND_PACKAGE(Boost COMPONENTS filesystem system thread timer chrono program_options unit_test_framework)
+ENDIF()
+
+IF(NOT find_cfitsio GREATER -1)
+  FIND_PACKAGE(Cfitsio)
+ENDIF()
+
+IF(NOT find_gsl GREATER -1)
+  FIND_PACKAGE(GSL 2.1)
+ENDIF()
+
+IF(NOT find_fftw GREATER -1)
+  FIND_PACKAGE(FFTW)
+ENDIF()
+
+IF(NOT GSL_FOUND AND NOT find_gsl GREATER -1)
   # older versions of cmake don't have FindGSL
   FIND_PACKAGE(PkgConfig)
   pkg_check_modules(GSL gsl>=2.1)
@@ -99,13 +115,17 @@ if(CCACHE_FOUND)
   set(ENV{CXX} "ccache ${CMAKE_CXX_COMPILER}")
 endif(CCACHE_FOUND)
 
+IF(MISSING_THIRDPARTIES)
+  LIST(REMOVE_DUPLICATES MISSING_THIRDPARTIES)
+ENDIF()
+
 #--------------------------------------------------------
 # Parse version, this will be used later to configure
 # cpf-core-redshift-version.cmake with the appropriate version number
 #--------------------------------------------------------
 
-FILE (STRINGS ${ROOT_DIR}/VERSION VERSION)
-MESSAGE( STATUS "Version used: ${VERSION}" )
+FILE (STRINGS ${ROOT_DIR}/VERSION PROJECT_VERSION)
+MESSAGE( STATUS "Version used: ${PROJECT_VERSION}" )
 
 #SET(CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} -Wall -Wno-sign-compare -Wno-parentheses")
 IF(BUILD_SHARED_LIBS)
@@ -136,8 +156,8 @@ LINK_DIRECTORIES( ${cpf-redshift_LINK_DIR} )
 # Exported configuration variables
 #--------------------------------------------------------
 
-IF( NOT (CFITSIO_FOUND AND FFTW_FOUND AND GSL_FOUND AND Boost_FOUND) )
-  MESSAGE( STATUS "Some library not found on system. Using thirdparty/" )
+IF( MISSING_THIRDPARTIES )
+  MESSAGE( STATUS "Some library not found on system. Using thirdparty/ for ${MISSING_THIRDPARTIES}" )
   set(PARALLEL_BUILD_FLAG "-j4")
   ADD_SUBDIRECTORY( tools )
 ELSE()
@@ -152,11 +172,6 @@ if (OPENMP_FOUND)
     set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
 endif()
 
-#-------------------------  -------------------------------
-# Extract version number
-#--------------------------------------------------------
-FILE (STRINGS ${ROOT_DIR}/VERSION VERSION_NUMBER)
-
 #--------------------------------------------------------
 # Compile library
 #--------------------------------------------------------
@@ -248,11 +263,11 @@ CONFIGURE_PACKAGE_CONFIG_FILE(
 )
 
 # Configure version file (this is used by find_package to determine a matching version)
-CONFIGURE_FILE(
-       ${CMAKE_CURRENT_SOURCE_DIR}/cmake/cpf-redshift-config-version.cmake.in
-       ${CMAKE_CURRENT_BINARY_DIR}/cmake/cpf-redshift-config-version.cmake
-       @ONLY
-)
+WRITE_BASIC_PACKAGE_VERSION_FILE(
+  ${CMAKE_CURRENT_BINARY_DIR}/cmake/cpf-redshift-config-version.cmake
+  VERSION ${VERSION}
+  COMPATIBILITY SameMajorVersion
+  )
 
 CONFIGURE_FILE(
 	${CMAKE_CURRENT_SOURCE_DIR}/cmake/version.h.in
@@ -262,6 +277,7 @@ CONFIGURE_FILE(
 # Generate TAGS file if etags is installed
 find_program(ETAGS NAMES etags ctags)
 if(ETAGS)
-  execute_process(COMMAND "find ${ROOT_DIR}/RedshiftLibrary/ ( -name *.cpp -o -name *.h -o -name *.hpp -o -name *.py -o -name *.c ) -exec ${ETAGS} -a ${ROOT_DIR}/TAGS {} +"
-    WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
+  message("Running " ${ETAGS})
+  execute_process(COMMAND find RedshiftLibrary/ ( -name *.cpp -o -name *.h -o -name *.hpp -o -name *.py -o -name *.c ) -not -name redshift.py -exec ${ETAGS} -a TAGS {} +
+    WORKING_DIRECTORY ${ROOT_DIR})
 endif(ETAGS)

README.md

@@ -46,6 +46,7 @@ Use `brew` as packet manager on MacOS:
 
 Use [Anaconda](https://www.anaconda.com/) as python3 provider and then install python dependencies with `pip`:
 
+    pip3 install numpy
     pip3 install astropy
 
 

RedshiftLibrary/RedshiftLibrary/linemodel/continuummodelsolution.h

@@ -25,6 +25,7 @@ struct CContinuumModelSolution
     Float64             tplMerit;
     Float64             tplDtm;
     Float64             tplMtm;
+    Float64             tplLogPrior;
 
     //polynom
     TFloat64List    pCoeffs;

RedshiftLibrary/RedshiftLibrary/linemodel/elementlist.h

@@ -16,6 +16,8 @@
 #include <RedshiftLibrary/linemodel/modelspectrumresult.h>
 #include <RedshiftLibrary/linemodel/element.h>
 
+#include <RedshiftLibrary/statistics/pdfz.h>
+
 #include <RedshiftLibrary/spectrum/template/catalog.h>
 #include <RedshiftLibrary/linemodel/templatesfitstore.h>
 
@@ -62,7 +64,7 @@ class CLineModelElementList
     void EstimateSpectrumContinuum(Float64 opt_enhance_lines, const TFloat64Range &lambdaRange);
 
     void LoadFitContinuumOneTemplate(const TFloat64Range& lambdaRange, const CTemplate& tpl);
-    void LoadFitContinuum(const TFloat64Range& lambdaRange, Int32 icontinuum);
+    void LoadFitContinuum(const TFloat64Range& lambdaRange, Int32 icontinuum, Int32 autoSelect);
     void setRedshift(Float64 redshift, bool reinterpolatedContinuum);
     Int32 ApplyContinuumOnGrid(const CTemplate& tpl, Float64 zcontinuum);
     Bool SolveContinuum(const CSpectrum& spectrum,
@@ -79,16 +81,20 @@ class CLineModelElementList
                         Float64 &fitDustCoeff,
                         Int32 &fitMeiksinIdx,
                         Float64& fitDtM,
-                        Float64& fitMtM);
+                        Float64& fitMtM,
+                        Float64 &fitLogprior);
     std::string getFitContinuum_tplName();
     Float64 getFitContinuum_tplAmplitude();
+    Float64 getFitContinuum_snr();
     Float64 getFitContinuum_tplMerit();
     void setFitContinuum_tplAmplitude(Float64 tplAmp, std::vector<Float64> polyCoeffs);
     Float64 getFitContinuum_tplIsmDustCoeff();
     Float64 getFitContinuum_tplIgmMeiksinIdx();
     Float64* getPrecomputedGridContinuumFlux();
     void SetContinuumComponent(std::string component);
     Int32 SetFitContinuum_FitStore(CTemplatesFitStore* fitStore);
+    Int32 SetFitContinuum_PriorHelper(CPriorHelperContinuum* priorhelper);
+    void SetFitContinuum_SNRMax(Float64 snr_max);
     void SetFitContinuum_Option(Int32 opt);
     Int32 GetFitContinuum_Option();
     void SetFitContinuum_FitValues(std::string tplfit_name,
@@ -99,6 +105,7 @@ class CLineModelElementList
                                    Float64 tplfit_continuumredshift,
                                    Float64 tplfit_dtm,
                                    Float64 tplfit_mtm,
+                                   Float64 tplfit_logprior,
                                    std::vector<Float64> polyCoeffs);
 
     Int32 LoadFitContaminantTemplate(const TFloat64Range& lambdaRange, const CTemplate& tpl);
@@ -117,9 +124,11 @@ class CLineModelElementList
     Float64 getTplshape_bestAmplitude();
     Int32 getTplshape_count();
     std::vector<Float64> getTplshape_priors();
+    std::vector<CPdfz::SPriorZ> getTplshape_priorsPz();
     std::vector<Float64> GetChisquareTplshape();
     std::vector<Float64> GetScaleMargTplshape();
     std::vector<bool> GetStrongELPresentTplshape();
+    std::vector<Int32> GetNLinesAboveSNRTplshape();
 
     Int32 GetNElements();
     Int32 GetModelValidElementsNDdl();
@@ -168,9 +177,10 @@ class CLineModelElementList
 
     void reinitModel();
     void refreshModel(Int32 lineTypeFilter=-1);
-    void reinitModelUnderElements(std::vector<UInt32> filterEltsIdx, Int32 lineIdx=-1 );
+    void getModel(CSpectrumFluxAxis& modelfluxAxis, Int32 lineTypeFilter=-1);
+    void reinitModelUnderElements(std::vector<UInt32> filterEltsIdx, Int32 lineIdx=-1);
     void refreshModelInitAllGrid();
-    void refreshModelUnderElements(std::vector<UInt32> filterEltsIdx, Int32 lineIdx=-1 );
+    void refreshModelUnderElements(std::vector<UInt32> filterEltsIdx, Int32 lineIdx=-1);
     void refreshModelDerivVelUnderElements(std::vector<UInt32> filterEltsIdx);
     void refreshModelDerivVelAbsorptionUnderElements(std::vector<UInt32> filterEltsIdx);
     void refreshModelDerivVelEmissionUnderElements(std::vector<UInt32> filterEltsIdx);
@@ -193,9 +203,11 @@ class CLineModelElementList
     Float64 getScaleMargCorrection(Int32 idxLine=-1);
     Float64 getContinuumScaleMargCorrection();
     Float64 getStrongerMultipleELAmpCoeff();
+    std::vector<std::string> getLinesAboveSNR(Float64 snrcut=3.5);
     Float64 getCumulSNRStrongEL();
     Float64 getCumulSNROnRange( TInt32Range idxRange );
     bool GetModelStrongEmissionLinePresent();
+    bool GetModelHaStrongest();
     Float64 getModelErrorUnderElement(UInt32 eltId);
     Float64 getContinuumMeanUnderElement(UInt32 eltId);
     Int32 LoadModelSolution(const CLineModelSolution&  modelSolution);
@@ -207,6 +219,7 @@ class CLineModelElementList
     Float64 GetContinuumError(Int32 eIdx, Int32 subeIdx);
     Int32 GetFluxDirectIntegration(TInt32List eIdx_list,
                                    TInt32List subeIdx_list,
+                                   Int32 opt_cont_substract_abslinesmodel,
                                    Float64& fluxdi,
                                    Float64& snrdi);
     const CSpectrumFluxAxis&    GetModelContinuum() const;
@@ -243,21 +256,34 @@ class CLineModelElementList
     TStringList GetModelRulesLog();
 
     Int32 setPassMode(Int32 iPass);
-
+    void SetForcedisableTplratioISMfit(bool opt);
 
     CRayCatalogsTplShape* m_CatalogTplShape;
     std::vector<Float64> m_ChisquareTplshape;
     std::vector<std::vector<Float64>> m_FittedAmpTplshape;
     std::vector<std::vector<Float64>> m_FittedErrorTplshape;
     std::vector<std::vector<Float64>> m_MtmTplshape;
     std::vector<std::vector<Float64>> m_DtmTplshape;
+    std::vector<std::vector<Float64>> m_LyaAsymCoeffTplshape;
+    std::vector<std::vector<Float64>> m_LyaWidthCoeffTplshape;
 
     bool m_enableAmplitudeOffsets;
     Float64 m_LambdaOffsetMin = -400.0;
     Float64 m_LambdaOffsetMax = 400.0;
     Float64 m_LambdaOffsetStep = 25.0;
     bool m_enableLambdaOffsetsFit;
 
+    bool m_opt_lya_forcefit=false;
+    bool m_opt_lya_forcedisablefit=false;
+    Float64 m_opt_lya_fit_asym_min=0.0;
+    Float64 m_opt_lya_fit_asym_max=4.0;
+    Float64 m_opt_lya_fit_asym_step=1.0;
+    Float64 m_opt_lya_fit_width_min=1.;
+    Float64 m_opt_lya_fit_width_max=4.;
+    Float64 m_opt_lya_fit_width_step=1.;
+    Float64 m_opt_lya_fit_delta_min=0.;
+    Float64 m_opt_lya_fit_delta_max=0.;
+    Float64 m_opt_lya_fit_delta_step=1.;
 
     Int32 m_opt_fitcontinuum_maxCount = 2;
     bool m_opt_firstpass_forcedisableMultipleContinuumfit=true;
@@ -283,6 +309,7 @@ class CLineModelElementList
                                                  Int32 polyOrder=-1);
 
     bool m_forceDisableLyaFitting;
+    bool m_forceLyaFitting=false;
     Int32 setLyaProfile( Float64 redshift, const CSpectrumSpectralAxis& spectralAxis );
 
     std::vector<UInt32> getSupportIndexes(std::vector<UInt32> EltsIdx);
@@ -304,6 +331,7 @@ class CLineModelElementList
 
     std::vector<Float64> m_ScaleMargCorrTplshape;
     std::vector<bool> m_StrongELPresentTplshape;
+    std::vector<Int32> m_NLinesAboveSNRTplshape;
 
     Float64 m_Redshift;
 
@@ -364,8 +392,12 @@ class CLineModelElementList
     Float64 m_fitContinuum_tplFitRedshift; // only used with m_fitContinuum_option==2 for now
     Float64 m_fitContinuum_tplFitDtM;
     Float64 m_fitContinuum_tplFitMtM;
+    Float64 m_fitContinuum_tplFitLogprior;
+    Float64 m_fitContinuum_tplFitSNRMax=0.0;
     std::vector<Float64> m_fitContinuum_tplFitPolyCoeffs;   // only used with m_fitContinuum_option==2 for now
     bool m_forcedisableMultipleContinuumfit=false;
+    Float64 m_fitContinuum_tplFitAlpha=0.;
+    CPriorHelperContinuum* m_fitContinuum_priorhelper;
 
     bool m_lmfit_noContinuumTemplate;
     bool m_lmfit_bestTemplate;

RedshiftLibrary/RedshiftLibrary/linemodel/linemodelextremaresult.h

@@ -54,6 +54,7 @@ class CLineModelExtremaResult : public COperatorResult
 
     //
     TFloat64List            StrongELSNR;
+    std::vector<std::vector<std::string>>            StrongELSNRAboveCut;
     TFloat64List            bic;    // bayesian information criterion for each extrema
     std::vector<CContinuumIndexes::TContinuumIndexList> ContinuumIndexes; //continuum indexes for each extrema
     std::vector<CMask>      OutsideLinesMask;   //Mask with 0 under the lines and 1 anywhere else
@@ -73,6 +74,9 @@ class CLineModelExtremaResult : public COperatorResult
     TFloat64List            FittedTplDustCoeff;     //Calzetti dustcoeff for the best template fitted for continuum
     std::vector<Int32>      FittedTplMeiksinIdx;    //Meiksin igm index for the best template fitted for continuum
     TFloat64List      FittedTplRedshift;    //Redshift for the best template fitted for continuum
+    TFloat64List      FittedTplDtm;    //DTM for the best template fitted for continuum
+    TFloat64List      FittedTplMtm;    //MTM for the best template fitted for continuum
+    TFloat64List      FittedTplLogPrior;    //log prior for the best template fitted for continuum
     std::vector<TFloat64List>      FittedTplpCoeffs;    //poly coeffs for the best template fitted for continuum
 
     //template ratio

RedshiftLibrary/RedshiftLibrary/linemodel/linemodelsolution.h

@@ -31,6 +31,7 @@ struct CLineModelSolution
     Float64 lfHa;
     Float64 snrOII;
     Float64 lfOII;
+    Int32 NLinesAboveSnrCut;
 
     std::vector<Float64> LambdaObs;  //observed position in Angstrom
     std::vector<Float64> Velocity;  //dispersion velocity in km/s

RedshiftLibrary/RedshiftLibrary/linemodel/modelcontinuumfittingresult.h

@@ -21,7 +21,13 @@ class CModelContinuumFittingResult : public COperatorResult
 
 public:
 
-    CModelContinuumFittingResult(Float64 _redshift, std::string _name, Float64 _merit, Float64 _amp, Float64 _ismCoeff, Int32 _igmIndex);
+    CModelContinuumFittingResult(Float64 _redshift,
+                                 std::string _name,
+                                 Float64 _merit,
+                                 Float64 _amp,
+                                 Float64 _ismCoeff,
+                                 Int32 _igmIndex,
+                                 Float64 _fitting_snr);
     CModelContinuumFittingResult();
     virtual ~CModelContinuumFittingResult();
 
@@ -41,6 +47,9 @@ class CModelContinuumFittingResult : public COperatorResult
     Float64 Amp;
     Float64 IsmCoeff;
     Int32   IgmIndex;
+
+    //fitting info
+    Float64 Fitting_snr;
 };
 
 }

RedshiftLibrary/RedshiftLibrary/linemodel/templatesfitstore.h

@@ -26,6 +26,7 @@ class CTemplatesFitStore
         Float64 fitAmplitude;
         Float64 fitDtM;
         Float64 fitMtM;
+        Float64 logprior;
     };
     typedef SValues TemplateFitValues;
 
@@ -39,7 +40,8 @@ class CTemplatesFitStore
              Float64 merit,
              Float64 fitAmplitude,
              Float64 fitDtM,
-             Float64 fitMtM);
+             Float64 fitMtM,
+             Float64 logprior);
 
     void prepareRedshiftList();
     void initFitValues();