diff --git a/cmake/FindGFlags.cmake b/cmake/FindGFlags.cmake index 9473cf23..6b2197f9 100644 --- a/cmake/FindGFlags.cmake +++ b/cmake/FindGFlags.cmake @@ -20,4 +20,4 @@ include (FindPackageHandleStandardArgs) find_package_handle_standard_args (Gflags DEFAULT_MSG GFLAGS_LIBRARIES GFLAGS_INCLUDE_DIRS) mark_as_advanced(GFLAGS_LIBRARY_DEBUG GFLAGS_LIBRARY_RELEASE - GFLAGS_LIBRARY GFLAGS_INCLUDE_DIR GFLAGS_ROOT_DIR) \ No newline at end of file + GFLAGS_LIBRARY GFLAGS_INCLUDE_DIR GFLAGS_ROOT_DIR) \ No newline at end of file diff --git a/src/fitting_ops.cxx b/src/fitting_ops.cxx index 399cd1ad..8d09731d 100644 --- a/src/fitting_ops.cxx +++ b/src/fitting_ops.cxx @@ -352,15 +352,21 @@ void _fit_noise_buffer(const bp::object & f, const bp::object & pxx, T* c_data = (T*)c_buf->buf; const int c_stride = c_buf->strides[0] / sizeof(T); - if constexpr (std::is_same::value) { + if constexpr (std::is_same::value) { + // Copy f to double + double f_double[nsamps]; + + std::transform(f_data, f_data + nsamps, f_double, + [](float value) { return static_cast(value); }); + // Get frequency bounds auto [lowf_i, fwhite_i, fwhite_size] = - _get_frequency_limits(f_data, lowf, fwhite_lower, fwhite_upper, nsamps); + _get_frequency_limits(f_double, lowf, fwhite_lower, fwhite_upper, nsamps); // Fit in logspace double log_f[nsamps]; for (int i = 0; i < nsamps; ++i) { - log_f[i] = std::log10(f_data[i]); + log_f[i] = std::log10(f_double[i]); } #pragma omp parallel for @@ -369,30 +375,30 @@ void _fit_noise_buffer(const bp::object & f, const bp::object & pxx, int poff = i * p_stride; int coff = i * c_stride; - double* pxx_det = pxx_data + ioff; + // Copy pxx row to double + double pxx_det[nsamps]; + + std::transform(pxx_data + ioff, pxx_data + ioff + nsamps, pxx_det, + [](float value) { return static_cast(value); }); + + // Cast implicitly on assignment T* p_det = p_data + poff; T* c_det = c_data + coff; - _fit_noise(f_data, log_f, pxx_det, p_det, c_det, - ndets, nsamps, lowf_i, fwhite_i, fwhite_size, - tol, niter, epsilon); + _fit_noise(f_double, log_f, pxx_det, p_det, + c_det, ndets, nsamps, lowf_i, fwhite_i, + fwhite_size, tol, niter, epsilon); } } - else if constexpr (std::is_same::value) { - // Copy f to double - double f_double[nsamps]; - - std::transform(f_data, f_data + nsamps, f_double, - [](float value) { return static_cast(value); }); - + else if constexpr (std::is_same::value) { // Get frequency bounds auto [lowf_i, fwhite_i, fwhite_size] = - _get_frequency_limits(f_double, lowf, fwhite_lower, fwhite_upper, nsamps); + _get_frequency_limits(f_data, lowf, fwhite_lower, fwhite_upper, nsamps); // Fit in logspace double log_f[nsamps]; for (int i = 0; i < nsamps; ++i) { - log_f[i] = std::log10(f_double[i]); + log_f[i] = std::log10(f_data[i]); } #pragma omp parallel for @@ -401,19 +407,13 @@ void _fit_noise_buffer(const bp::object & f, const bp::object & pxx, int poff = i * p_stride; int coff = i * c_stride; - // Copy pxx row to double - double pxx_det[nsamps]; - - std::transform(pxx_data + ioff, pxx_data + ioff + nsamps, pxx_det, - [](float value) { return static_cast(value); }); - - // Cast implicitly on assignment + double* pxx_det = pxx_data + ioff; T* p_det = p_data + poff; T* c_det = c_data + coff; - _fit_noise(f_double, log_f, pxx_det, p_det, - c_det, ndets, nsamps, lowf_i, fwhite_i, - fwhite_size, tol, niter, epsilon); + _fit_noise(f_data, log_f, pxx_det, p_det, c_det, + ndets, nsamps, lowf_i, fwhite_i, fwhite_size, + tol, niter, epsilon); } } } @@ -436,20 +436,21 @@ void fit_noise(const bp::object & f, const bp::object & pxx, bp::object & p, bp: } } + PYBINDINGS("so3g") { bp::def("fit_noise", fit_noise, "fit_noise(f, pxx, p, c, lowf, fwhite_lower, fwhite_upper, tol, niter, epsilon)" "\n" - "Fits noise model with white and 1/f noise to the PSD of signal. Uses a MLE " + "Fits noise model with white and 1/f components to the PSD of signal. Uses a MLE " "method that minimizes a log likelihood. OMP is used to parallelize across dets (rows)." "\n" "Args:\n" - " f: frequency array (float32/64) with dimensions (nsamps,)\n" - " pxx: PSD array (float32/64) with dimensions (nsamps,)\n" - " p: parameter array (float32/64) with dimensions (nparams,). This is modified in place.\n" - " c: uncertaintiy array (float32/64) with dimensions (nsamps,). This is modified in place.\n" - " lowf: Frequency below which estimate of 1/f noise index and knee are estimated for initial " + " f: frequency array (float32/64) with dimensions (nsamps,).\n" + " pxx: PSD array (float32/64) with dimensions (ndets, nsamps).\n" + " p: parameter array (float32/64) with dimensions (ndets, nparams). This is modified in place.\n" + " c: uncertaintiy array (float32/64) with dimensions (ndets, nparams). This is modified in place.\n" + " lowf: Frequency below which the 1/f noise index and fknee are estimated for initial " " guess passed to least_squares fit (float64).\n" " fwhite_lower: Lower frequency used to estimate white noise for initial guess passed to " " least_squares fit (float64).\n" @@ -458,5 +459,5 @@ PYBINDINGS("so3g") " tol: absolute tolerance for minimization (float64).\n" " niter: total number of iterations to run minimization for (int).\n" " epsilon: Value to perturb gradients by when calculating uncertainties with the inverse " - " Hessian matrix (float64).\n"); + " Hessian matrix (float64). Affects minimization only.\n"); } \ No newline at end of file