From 5334b8d34b3a385e0fd1421c852556c15f1da1c0 Mon Sep 17 00:00:00 2001 From: Jared Males Date: Mon, 6 Jul 2026 00:58:19 +0200 Subject: [PATCH] codex hunting segfaults --- include/ao/analysis/fourierTemporalPSD.hpp | 254 +++++++++++++++------ include/ao/analysis/speckleAmpPSD.hpp | 69 +++++- include/math/ft/fftT.hpp | 7 + include/sigproc/averagePeriodogram.hpp | 5 + 4 files changed, 256 insertions(+), 79 deletions(-) diff --git a/include/ao/analysis/fourierTemporalPSD.hpp b/include/ao/analysis/fourierTemporalPSD.hpp index de0aea1fd..63a587fbb 100644 --- a/include/ao/analysis/fourierTemporalPSD.hpp +++ b/include/ao/analysis/fourierTemporalPSD.hpp @@ -770,14 +770,61 @@ int fourierTemporalPSD::analyzePSDGrid( const std::string &subDir bool writeXfer ) { + if( m_aosys == nullptr ) + { + std::cerr << "fourierTemporalPSD::analyzePSDGrid: m_aosys is null.\n"; + return -1; + } + + if( mnMax <= 0 ) + { + std::cerr << "fourierTemporalPSD::analyzePSDGrid: mnMax must be > 0.\n"; + return -1; + } + + if( mnCon < 0 ) + { + std::cerr << "fourierTemporalPSD::analyzePSDGrid: mnCon must be >= 0.\n"; + return -1; + } + + if( lpNc < 0 ) + { + std::cerr << "fourierTemporalPSD::analyzePSDGrid: lpNc must be >= 0.\n"; + return -1; + } + + if( lifetimeTrials < 0 ) + { + std::cerr << "fourierTemporalPSD::analyzePSDGrid: lifetimeTrials must be >= 0.\n"; + return -1; + } + + if( mags.empty() ) + { + std::cerr << "fourierTemporalPSD::analyzePSDGrid: at least one guide star magnitude is required.\n"; + return -1; + } + + if( m_aosys->tauWFS() <= 0 ) + { + std::cerr << "fourierTemporalPSD::analyzePSDGrid: tauWFS must be > 0.\n"; + return -1; + } + std::string dir = psdDir + "/" + subDir; /*** Dump Params to file ***/ - mkdir( dir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH ); + ioutils::createDirectories( dir ); std::ofstream fout; std::string fn = dir + '/' + "params.txt"; fout.open( fn ); + if( !fout.good() ) + { + std::cerr << "fourierTemporalPSD::analyzePSDGrid: could not open " << fn << " for writing.\n"; + return -1; + } fout << "#---------------------------\n"; m_aosys->dumpAOSystem( fout ); @@ -787,9 +834,15 @@ int fourierTemporalPSD::analyzePSDGrid( const std::string &subDir fout << "# mnCon = " << mnCon << '\n'; fout << "# lpNc = " << lpNc << '\n'; fout << "# mags = "; - for( size_t i = 0; i < mags.size() - 1; ++i ) - fout << mags[i] << ","; - fout << mags[mags.size() - 1] << '\n'; + for( size_t i = 0; i < mags.size(); ++i ) + { + if( i > 0 ) + { + fout << ","; + } + fout << mags[i]; + } + fout << '\n'; fout << "# lifetimeTrials = " << lifetimeTrials << '\n'; fout << "# uncontrolledLifetimes = " << uncontrolledLifetimes << '\n'; fout << "# writePSDs = " << std::boolalpha << writePSDs << '\n'; @@ -803,6 +856,35 @@ int fourierTemporalPSD::analyzePSDGrid( const std::string &subDir realT tauWFS = m_aosys->tauWFS(); realT deltaTau = m_aosys->deltaTau(); + std::vector gridFreq; + if( getGridFreq( gridFreq, psdDir ) < 0 ) + { + return -1; + } + + if( gridFreq.size() < 2 ) + { + std::cerr << "fourierTemporalPSD::analyzePSDGrid: PSD frequency grid must have at least two samples.\n"; + return -1; + } + + size_t imax = 0; + while( imax < gridFreq.size() && gridFreq[imax] <= 0.5 * fs ) + { + ++imax; + } + + if( imax < gridFreq.size() - 1 && gridFreq[imax] <= 0.5 * fs * ( 1.0 + 1e-7 ) ) + { + ++imax; + } + + if( imax < 2 ) + { + std::cerr << "fourierTemporalPSD::analyzePSDGrid: PSD frequency grid has fewer than two samples below the loop Nyquist frequency.\n"; + return -1; + } + std::vector fms; sigproc::makeFourierModeFreqs_Rect( fms, 2 * mnMax ); @@ -813,18 +895,16 @@ int fourierTemporalPSD::analyzePSDGrid( const std::string &subDir gains.resize( 2 * mnMax + 1, 2 * mnMax + 1 ); vars.resize( 2 * mnMax + 1, 2 * mnMax + 1 ); speckleLifetimes.resize( 2 * mnMax + 1, 2 * mnMax + 1 ); - - gains( mnMax, mnMax ) = 0; - vars( mnMax, mnMax ) = 0; - speckleLifetimes( mnMax, mnMax ) = 0; + gains.setZero(); + vars.setZero(); + speckleLifetimes.setZero(); gains_lp.resize( 2 * mnMax + 1, 2 * mnMax + 1 ); vars_lp.resize( 2 * mnMax + 1, 2 * mnMax + 1 ); speckleLifetimes_lp.resize( 2 * mnMax + 1, 2 * mnMax + 1 ); - - gains_lp( mnMax, mnMax ) = 0; - vars_lp( mnMax, mnMax ) = 0; - speckleLifetimes_lp( mnMax, mnMax ) = 0; + gains_lp.setZero(); + vars_lp.setZero(); + speckleLifetimes_lp.setZero(); bool doLP = false; if( lpNc > 1 ) @@ -843,15 +923,15 @@ int fourierTemporalPSD::analyzePSDGrid( const std::string &subDir { for( size_t s = 0; s < mags.size(); ++s ) { - std::string psdOutDir = std::format("{}/outputPSDS_{}_si",dir, mags[s]); + std::string psdOutDir = std::format("{}/outputPSDs_{}_si",dir, mags[s]); //dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_si"; - mkdir( psdOutDir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH ); + ioutils::createDirectories( psdOutDir ); if( doLP ) { - std::string psdOutDir = std::format("{}/outputPSDS_{}_lp",dir, mags[s]); + std::string psdOutDir = std::format("{}/outputPSDs_{}_lp",dir, mags[s]); //dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_lp"; - mkdir( psdOutDir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH ); + ioutils::createDirectories( psdOutDir ); } } } @@ -898,6 +978,8 @@ int fourierTemporalPSD::analyzePSDGrid( const std::string &subDir m_aosys->optd( m_aosys->optd() ); // just trigger a re-calc realT strehl = m_aosys->strehl(); + bool gridReadFailed = false; + #pragma omp parallel { realT localMag = mags[s]; @@ -918,28 +1000,13 @@ int fourierTemporalPSD::analyzePSDGrid( const std::string &subDir // PSDs. //**< Get the frequency grid, and nyquist limit it to f_s/2 - getGridPSD( tfreq, tPSDp, psdDir, 0, 1 ); // To get the freq grid - - size_t imax = 0; - while( tfreq[imax] <= 0.5 * fs ) - { - ++imax; - if( imax > tfreq.size() - 1 ) - break; - } - - if( imax < tfreq.size() - 1 && tfreq[imax] <= 0.5 * fs * ( 1.0 + 1e-7 ) ) - { - ++imax; - } + tfreq.assign( gridFreq.begin(), gridFreq.begin() + imax ); if( writePSDs ) { - tfreqHF.assign( tfreq.begin(), tfreq.end() ); + tfreqHF.assign( gridFreq.begin(), gridFreq.end() ); } - tfreq.erase( tfreq.begin() + imax, tfreq.end() ); - tPSDpPOL.resize( tfreq.size() ); // pre=allocate //**> @@ -969,26 +1036,23 @@ int fourierTemporalPSD::analyzePSDGrid( const std::string &subDir std::vector> ETFxn; std::vector> NTFxn; - if( lifetimeTrials > 0 ) + if( lifetimeTrials > 0 || writeXfer ) { ETFxn.resize( tfreq.size() ); NTFxn.resize( tfreq.size() ); + } - if( writeXfer ) - { - std::string tfOutFile = std::format("{}/outputTF_{}_si", dir, mags[s]); - //dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_si/"; - ioutils::createDirectories( tfOutFile ); - } + if( writeXfer ) + { + std::string tfOutFile = std::format("{}/outputTF_{}_si", dir, mags[s]); + //dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_si/"; + ioutils::createDirectories( tfOutFile ); if( doLP ) { - if( writeXfer ) - { - std::string tfOutFile = std::format("{}/outputTF_{}_lp", dir, mags[s]); - //dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_lp/"; - ioutils::createDirectories( tfOutFile ); - } + tfOutFile = std::format("{}/outputTF_{}_lp", dir, mags[s]); + //dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_lp/"; + ioutils::createDirectories( tfOutFile ); } } @@ -1028,10 +1092,29 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; else { - realT k = sqrt( m * m + n * n ) / m_aosys->D(); - //**< Get the open-loop turb. PSD - getGridPSD( tPSDp, psdDir, m, n ); + if( getGridPSD( tPSDp, psdDir, m, n ) < 0 ) + { +#pragma omp critical + { + gridReadFailed = true; + } + watcher.incrementAndOutputStatus(); + continue; + } + + if( tPSDp.size() != gridFreq.size() ) + { +#pragma omp critical + { + std::cerr << "fourierTemporalPSD::analyzePSDGrid: PSD size mismatch for mode (" + << m << ", " << n << "). Expected " << gridFreq.size() + << " samples, got " << tPSDp.size() << ".\n"; + gridReadFailed = true; + } + watcher.incrementAndOutputStatus(); + continue; + } // Get integral of entire open-loop PSD var0 = sigproc::psdVar( tfreq, tPSDp ); @@ -1174,6 +1257,7 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; else { gopt_lp = 0; + var_lp = var; } } else @@ -1260,23 +1344,36 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; if( i == 0 ) // Write freq on the first one { - std::string fOutFile = tfOutFile + "freq.binv"; + std::string fOutFile = tfOutFile + "/freq.binv"; ioutils::writeBinVector( fOutFile, tfreq ); } } if( lifetimeTrials > 0 ) { - speckleAmpPSD( spfreq, sppsd, tfreq, tPSDp, ETFxn, tPSDn, NTFxn, lifetimeTrials ); - realT spvar = sigproc::psdVar( spfreq, sppsd ); + if( speckleAmpPSD( spfreq, sppsd, tfreq, tPSDp, ETFxn, tPSDn, NTFxn, lifetimeTrials ) < 0 ) + { +#pragma omp critical + { + gridReadFailed = true; + } + } + else + { + realT spvar = sigproc::psdVar( spfreq, sppsd ); - realT splifeT = 100.0; - realT error; + realT splifeT = 100.0; + realT error; - realT tau = pvm( error, spfreq, sppsd, splifeT ) * ( splifeT ) / spvar; + realT tau = 0; + if( spvar > 0 ) + { + tau = pvm( error, spfreq, sppsd, splifeT ) * ( splifeT ) / spvar; + } - speckleLifetimes( mnMax + m, mnMax + n ) = tau; - speckleLifetimes( mnMax - m, mnMax - n ) = tau; + speckleLifetimes( mnMax + m, mnMax + n ) = tau; + speckleLifetimes( mnMax - m, mnMax - n ) = tau; + } } if( doLP ) @@ -1315,23 +1412,36 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; if( i == 0 ) // Write freq on the first one { - std::string fOutFile = tfOutFile + "freq.binv"; + std::string fOutFile = tfOutFile + "/freq.binv"; ioutils::writeBinVector( fOutFile, tfreq ); } } if( lifetimeTrials > 0 ) { - speckleAmpPSD( spfreq, sppsd, tfreq, tPSDp, ETFxn, tPSDn, NTFxn, lifetimeTrials ); - realT spvar = sigproc::psdVar( spfreq, sppsd ); + if( speckleAmpPSD( spfreq, sppsd, tfreq, tPSDp, ETFxn, tPSDn, NTFxn, lifetimeTrials ) < 0 ) + { +#pragma omp critical + { + gridReadFailed = true; + } + } + else + { + realT spvar = sigproc::psdVar( spfreq, sppsd ); - realT splifeT = 100.0; - realT error; + realT splifeT = 100.0; + realT error; - realT tau = pvm( error, spfreq, sppsd, splifeT ) * ( splifeT ) / spvar; + realT tau = 0; + if( spvar > 0 ) + { + tau = pvm( error, spfreq, sppsd, splifeT ) * ( splifeT ) / spvar; + } - speckleLifetimes_lp( mnMax + m, mnMax + n ) = tau; - speckleLifetimes_lp( mnMax - m, mnMax - n ) = tau; + speckleLifetimes_lp( mnMax + m, mnMax + n ) = tau; + speckleLifetimes_lp( mnMax - m, mnMax - n ) = tau; + } } } // if(doLP) } // if( (lifetimeTrials > 0 || writeXfer) && ( uncontrolledLifetimes || inside )) @@ -1439,6 +1549,11 @@ std::cerr << __FILE__ << " " << __LINE__ << "\n"; } // omp for i..nModes } // omp Parallel + if( gridReadFailed ) + { + return -1; + } + Eigen::Array cim; fits::fitsFile ff; @@ -1780,7 +1895,7 @@ int fourierTemporalPSD::intensityPSD( } else { - for( int q = 0; q < ETFsi.size(); ++q ) + for( size_t q = 0; q < ETFsi[i].size(); ++q ) { ETFsi[i][q] = 1; NTFsi[i][q] = 0; @@ -1791,7 +1906,8 @@ int fourierTemporalPSD::intensityPSD( } sigproc::averagePeriodogram tavgPgram( - sz2Sided / 1., + sz2Sided, + 0, 1 / fs ); // this is just to get the size right, per-thread instances below std::vector> spPSDs; spPSDs.resize( nModes ); @@ -1862,7 +1978,7 @@ int fourierTemporalPSD::intensityPSD( N_nm.resize( sz2Sided ); // Periodogram averager - sigproc::averagePeriodogram avgPgram( sz2Sided / 1., 1 / fs ); //, 0, 1); + sigproc::averagePeriodogram avgPgram( sz2Sided, 0, 1 / fs ); //, 0, 1); avgPgram.win( sigproc::window::hann ); // The periodogram output @@ -1987,7 +2103,7 @@ int fourierTemporalPSD::intensityPSD( // Filter it // Set sqrt(PSD), just a pointer switch - pfilt.psdSqrt( &sqrtNPSD[pp], tfreq[1] - tfreq[0] ); + nfilt.psdSqrt( &sqrtNPSD[pp], tfreq[1] - tfreq[0] ); nfilt.filter( N_n ); nfilt.filter( N_nm ); diff --git a/include/ao/analysis/speckleAmpPSD.hpp b/include/ao/analysis/speckleAmpPSD.hpp index 74d681164..04ece2191 100644 --- a/include/ao/analysis/speckleAmpPSD.hpp +++ b/include/ao/analysis/speckleAmpPSD.hpp @@ -8,8 +8,14 @@ #ifndef speckleAmpPSD_hpp #define speckleAmpPSD_hpp +#include +#include #include +#ifdef _OPENMP +#include +#endif + #include #include "../../math/constants.hpp" @@ -64,6 +70,25 @@ int speckleAmpPSD( ) { + if( N <= 0 ) + { + std::cerr << "speckleAmpPSD: N must be > 0.\n"; + return -1; + } + + if( freq.size() < 2 ) + { + std::cerr << "speckleAmpPSD: frequency grid must have at least two samples.\n"; + return -1; + } + + if( fmPSD.size() != freq.size() || nPSD.size() != freq.size() || fmXferFxn.size() != freq.size() || + nXferFxn.size() != freq.size() ) + { + std::cerr << "speckleAmpPSD: input vectors must all match the frequency grid size.\n"; + return -1; + } + std::vector psd2, npsd2; std::vector> xfer2, nxfer2; @@ -89,7 +114,9 @@ int speckleAmpPSD( int Nsamp = 1.0 * psd2.size(); int NsampStart = 0.5 * Nwd - 0.5 * Nsamp; - sigproc::averagePeriodogram globalAvgPgram( Nsamp * 0.1, dt ); + size_t avgLen = std::max( 1, static_cast( Nsamp * 0.1 ) ); + + sigproc::averagePeriodogram globalAvgPgram( avgLen, 0, dt ); spPSD.resize( globalAvgPgram.size() ); for( size_t n = 0; n < spPSD.size(); ++n ) spPSD[n] = 0; @@ -103,7 +130,11 @@ int speckleAmpPSD( // and the noise variance realT nVar = sigproc::psdVar( freq, nPSD ); +#ifdef _OPENMP +#pragma omp parallel if( !omp_in_parallel() ) +#else #pragma omp parallel +#endif { // Filters for imposing the PSDs mx::sigproc::psdFilter filt; @@ -140,7 +171,7 @@ int speckleAmpPSD( std::vector vn( Nsamp ); // Periodogram averager - sigproc::averagePeriodogram avgPgram( Nsamp * 0.1, dt ); //, 0, 1); + sigproc::averagePeriodogram avgPgram( avgLen, 0, dt ); //, 0, 1); avgPgram.win( sigproc::window::hann ); // The temporary periodogram @@ -161,9 +192,16 @@ int speckleAmpPSD( filt( fm_n ); math::vectorMeanSub( fm_n ); realT actvar = math::vectorVariance( fm_n ); - realT norm = sqrt( fmVar / actvar ); - for( size_t q = 0; q < fm_n.size(); ++q ) - fm_n[q] *= norm; + if( actvar > 0 && fmVar > 0 ) + { + realT norm = sqrt( fmVar / actvar ); + for( size_t q = 0; q < fm_n.size(); ++q ) + fm_n[q] *= norm; + } + else + { + std::fill( fm_n.begin(), fm_n.end(), 0 ); + } // And move it to the Fourier domain fft( tform1.data(), fm_n.data() ); @@ -173,11 +211,19 @@ int speckleAmpPSD( nfilt.filter( N_nm ); realT Nactvar = 0.5 * ( math::vectorVariance( N_n ) + math::vectorVariance( N_nm ) ); - norm = sqrt( nVar / Nactvar ); - for( size_t q = 0; q < fm_n.size(); ++q ) - N_n[q] *= norm; - for( size_t q = 0; q < fm_n.size(); ++q ) - N_nm[q] *= norm; + if( Nactvar > 0 && nVar > 0 ) + { + realT norm = sqrt( nVar / Nactvar ); + for( size_t q = 0; q < fm_n.size(); ++q ) + N_n[q] *= norm; + for( size_t q = 0; q < fm_n.size(); ++q ) + N_nm[q] *= norm; + } + else + { + std::fill( N_n.begin(), N_n.end(), 0 ); + std::fill( N_nm.begin(), N_nm.end(), 0 ); + } // And move them to the Fourier domain fft( Ntform1.data(), N_n.data() ); @@ -225,7 +271,10 @@ int speckleAmpPSD( // Calculate PSD of the speckle amplitude if( !noPSD ) + { + std::fill( tpgram.begin(), tpgram.end(), 0 ); avgPgram( tpgram, vn ); + } // Accumulate #pragma omp critical diff --git a/include/math/ft/fftT.hpp b/include/math/ft/fftT.hpp index ab8c0ba98..7bb7a7f3c 100644 --- a/include/math/ft/fftT.hpp +++ b/include/math/ft/fftT.hpp @@ -240,7 +240,14 @@ void fftT::destroyPlan() { if( m_plan ) { +#ifndef MX_FFTW_NOOMP +#pragma omp critical + { +#endif fftw_destroy_plan( m_plan ); +#ifndef MX_FFTW_NOOMP + } +#endif } m_plan = 0; diff --git a/include/sigproc/averagePeriodogram.hpp b/include/sigproc/averagePeriodogram.hpp index 674e6d025..3f7f149a3 100644 --- a/include/sigproc/averagePeriodogram.hpp +++ b/include/sigproc/averagePeriodogram.hpp @@ -405,6 +405,11 @@ void averagePeriodogram::operator()( realT *pgram, const realT *ts, size_ std::cerr << "averagePeriodogram: Window size not correct.\n"; } + for( size_t j = 0; j < m_size; ++j ) + { + pgram[j] = 0; + } + int Navg = sz / m_nOver; while( Navg * m_nOver + m_avgLen > sz )