Merge pull request #76 from FastNFT/kdvv_improvements

Kdvv improvements

Author: wahls

CHANGELOG.md

@@ -1,14 +1,36 @@
 # Changelog
 
+## [0.5.0]
+
+### Added
+
+- A NFT routine fnft_manakovv for the Manakov equation with vanishing boundaries was added (continuous spectrum only).
+- The periodic NFT routine fnft_nsep now also supports pure Newton refinement.
+- The routine fnft_kdvv now can now also compute the discrete spectrum. To locate the bound states, either Newton's method or a grid search with additional Newton refinements are available.
+- The slow scattering methods for AKNS-type systems now include a normalization procedure to deal with numerical overflow (enabled by default).
+- The routine fnft__kdv_finvscatter has been added. The plan is to later use it for a fast inverse KdV NFT.
+
+### Changed
+
+- New criteria for stopping Newton iteration in fnft_nsep.
+- Simplified some tests to reduce run times.
+- The code for AKNS scattering has been overhauled.
+
+### Fixed
+
+- Several memory leaks have been fixed.
+
 ## [0.4.1] -- 2020-07-13
 
 ### Changed
 
 - Number of samples for fnft_nsep again has to be a power of two.
+- misc_resample no longer issues a warning when the signal appears to be undersampled. This gave the wrong impression that CFx_y discrizations suffer more in such scenarios than the other ones, which do not use this routine.
 
 ### Fixed
 
-- misc_downsample could return incorrect values for first_last_index[1]
+- misc_downsample could return incorrect values for first_last_index[1].
+- Some errors in fnft_nsev become meaningless when bound_states==NULL and should not be risen in that case.
 
 ## [0.4.0] -- 2020-07-08
 

CMakeLists.txt

@@ -20,9 +20,9 @@ cmake_minimum_required(VERSION 3.17.3)
 project(fnft C)
 
 set(FNFT_VERSION_MAJOR 0)
-set(FNFT_VERSION_MINOR 4)
-set(FNFT_VERSION_PATCH 1)
-set(FNFT_VERSION_SUFFIX "") # should not be longer than FNFT_SUFFIX_MAXLEN
+set(FNFT_VERSION_MINOR 5)
+set(FNFT_VERSION_PATCH 0)
+set(FNFT_VERSION_SUFFIX "dev") # should not be longer than FNFT_SUFFIX_MAXLEN
 set(FNFT_VERSION ${FNFT_VERSION_MAJOR}.${FNFT_VERSION_MINOR}.${FNFT_VERSION_PATCH}${FNFT_VERSION_SUFFIX})
 
 include(CheckIncludeFiles)
@@ -223,4 +223,4 @@ if (WITH_MATLAB)
       # -DMEX_DOUBLE_HANDLE avoids the new complex interleaved API
     endforeach()
   endif()
-endif()
+endif()

INSTALL.md

@@ -78,7 +78,7 @@ in the PowerShell to install Scoop. Then, run the commands
 
     scoop install git
     scoop install cmake
-    scoop install gcc
+    scoop install mingw
 
 in the PowerShell to install the required tools. Change (for example) to your home directory,
 download FNFT, and change into the new dir:

README.md

@@ -118,3 +118,4 @@ The algorithms in FNFT utilize ideas from the following references. More informa
 - P. J. Prins and S. Wahls, ["Soliton Phase Shift Calculation for the Korteweg–De Vries Equation"](https://doi.org/10.1109/ACCESS.2019.2932256), IEEE Access, vol. 7, pp. 122914--122930, July 2019.
 - S. Medvedev, I. Vaseva, I. Chekhovskoy and M. Fedoruk, ["Exponential fourth order schemes for direct Zakharov-Shabat problem"](https://doi.org/10.1364/OE.377140), Optics Express, vol. 28, pp. 20--39, 2020.
 - J. Mertsching, ["Quasiperiodie Solutions of the Nonlinear Schroedinger Equation"](https://doi.org/10.1002/prop.2190350704), Fortschritte der Physik, vol. 35, pp. 519--536, 1987.
+- L. de Vries, ["Fast Numerical Nonlinear Fourier Transform Algorithms for the Manakov Equation"](http://resolver.tudelft.nl/uuid:0276e693-3408-4472-9749-b754c2114183"), Master thesis, TU Delft, 2021.

examples/mex_fnft_kdvv_example.m

@@ -39,7 +39,10 @@
 
 %%% Compute the nonlinear Fourier transform %%%
 
-[contspec, bound_states, norming_constants] = mex_fnft_kdvv(q, T, XI);
+dxi = sqrt(max(q)) / 1000; % use approximately 1000 grid points for the
+% bound state localization step
+[contspec, bound_states, norming_constants] = mex_fnft_kdvv(q, T, XI, ...
+    'grid_spacing', dxi);
 % mex_fnft_kdvv has many options => run "help mex_fnft_kdvv" to learn more
 
 %%% Plot the results %%%

examples/mex_fnft_kdvv_example_4.m

@@ -46,8 +46,11 @@
 
 t = linspace(T(1),T(2),D);
 q = double(q_fun(t)); % signal samples
+dxi = sqrt(max(q)) / 1000; % use approximately 1000 grid points for the
+% bound state localization step
 [~,bound_states_computed,normconsts_computed]=mex_fnft_kdvv(q, T, XI,... 
-    'discr_CF4_2','bsloc_gridsearch_refine','skip_cs', 'bsloc_niter',20);
+    'discr_CF4_2','bsloc_gridsearch_refine', 'grid_spacing', dxi, ...
+    'skip_cs', 'bsloc_niter',20);
 
 %%% Plot results %%%
 

include/fnft_kdvv.h

@@ -14,7 +14,7 @@
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 *
 * Contributors:
-* Sander Wahls (TU Delft) 2017-2018.
+* Sander Wahls (TU Delft) 2017-2018, 2023.
 * Shrinivas Chimmalgi (TU Delft) 2019-2020.
 * Peter J Prins (TU Delft) 2020-2021.
 */
@@ -54,9 +54,13 @@
  *  \f$ O(niter (*K\_ptr) D) \f$. \n \n
  *  fnft_kdvv_bsloc_GRIDSEARCH_AND_REFINE: The algorithm evaluates \f$ a(\xi) \f$ on the grid
  *  \f$ \xi = \{j nexttoward(0,1), jh,2jh,3jh,\ldots,(M-3)jh\,(M-2)jh,j nexttoward((M-1)h,0)\\}\f$,
- *  where \f$ h:= \sqrt{c \max_t q(t)} / (M-1)\f$, where \f$ M=1000\f$ and
- *  \f$ c=1\f$ for all second order discretizations, `fnft_kdv_discretization_4SPLIT4A'/'B'('_VANILLA'), and `fnft_kdv_discretization_CF4_2`(`_VANILLA`); \f$ c \f$ is approximately 2 for
- *  other discretizations. The sign changes of \f$ a(\xi) \f$ on this grid are used as initial
+ *  where \f$ h:= \sqrt{c \max_t q(t)} / (M-1)\f$, where \f$ c=1\f$ for all second order discretizations,
+ *  `fnft_kdv_discretization_4SPLIT4A'/'B'('_VANILLA'), and `fnft_kdv_discretization_CF4_2`(`_VANILLA`);
+ *  \f$ c \f$ is approximately 2 for other discretizations. The constant \f$ M \f$ is chosen
+ *  such that the distance between the neighboring grid points is not larger than the
+ *  parameter \link fnft_kdvv_opts_t::grid_spacing \endlink. This parameter therefore must
+ *  be set if this algorithm is used.
+ *  The sign changes of \f$ a(\xi) \f$ on this grid are used as initial
  *  guesses for the bound states, which are then refined as in `fnft_kdvv_bsloc_NEWTON`.
  */
 typedef enum {
@@ -119,7 +123,7 @@ typedef enum {
  * @var fnft_kdvv_opts_t::niter
  *  Number of Newton iterations to be carried out when either the
  *  fnft_kdvv_bsloc_NEWTON, or the fnft_kdvv_bsloc_GRIDSEARCH_AND_REFINE method
- *  is used.
+ *  is used. Can be zero or positive.
  *
  * @var fnft_kdvv_opts_t::discspec_type
  *  Controls how \link fnft_kdvv \endlink fills the array
@@ -141,7 +145,7 @@ typedef enum {
  *  Controls which discretization is applied to the continuous-time Zakharov-
  *  Shabat scattering problem. See \link fnft_kdv_discretization_t \endlink.\n\n
  *
- *  @var fnft_kdvv_opts_t::richardson_extrapolation_flag
+ * @var fnft_kdvv_opts_t::richardson_extrapolation_flag
  *  Controls whether Richardson extrapolation is applied to try and improve
  *  the accuracy of the computed spectrum. First approximation is computed
  *  as usual using all the supplied samples. A second approximation is computed
@@ -154,7 +158,12 @@ typedef enum {
  *  worse accuracy compared to the first approximation.
  *  By default, Richardson extrapolation is disabled (i.e., the
  *  flag is zero). To enable, set the flag to one.
- */
+ *
+ * @var fnft_kdvv_opts_t::grid_spacing
+ *   Grid spacing parameter for the fnft_kdvv_bsloc_GRIDSEARCH_AND_REFINE method,
+ *   where the number of grid points is chosen such that the distance between two
+ *   consequtive grid points is not larger than grid_spacing, which has to be
+ *   positive. This option has to be set if GRIDSEARCH_AND_REFINE is used. */
 typedef struct {
     fnft_kdvv_bsloc_t bound_state_localization;
     FNFT_UINT niter;
@@ -163,6 +172,7 @@ typedef struct {
     FNFT_INT normalization_flag;
     fnft_kdv_discretization_t discretization;
     FNFT_UINT richardson_extrapolation_flag;
+    FNFT_REAL grid_spacing;
 } fnft_kdvv_opts_t;
 
 /**

include/private/fnft__akns_scatter.h

@@ -14,9 +14,10 @@
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 *
 * Contributors:
-* Sander Wahls (TU Delft) 2017-2018.
+* Sander Wahls (TU Delft) 2017-2018, 2023.
 * Shrinivas Chimmalgi (TU Delft) 2017-2020.
 * Peter J. Prins (TU Delft) 2021.
+* Sander Wahls (KIT), 2023.
 */
 
 /**
@@ -59,6 +60,10 @@
  * If derivative_flag=1 returns [S11 S12 S21 S22 S11' S12' S21' S22'] in 
  * result where S11' is the derivative of S11 w.r.t to \f$\lambda\f$.
  * Should be preallocated with size 4*K or 8*K accordingly.
+ * @param[in,out] W Pass an array of size K. Upon exit, it contains scaling factors that
+ * arise due to an internal normalization of the scattering matrix (to deal with potential
+ * overflow issues). The result has to be scaled by the power of 2 to obtain the final
+ * values. No normalization is carried out if NULL is passed.
  * @param[in] discretization The type of discretization to be used. Should be of type 
  * \link fnft__akns_discretization_t \endlink. Not all akns_discretization_t 
  * discretizations are supported. Check \link fnft__akns_discretization_t \endlink 
@@ -79,6 +84,7 @@ FNFT_INT fnft__akns_scatter_matrix(FNFT_UINT const D,
                                    FNFT_UINT const K,
                                    FNFT_COMPLEX const * const lambda,
                                    FNFT_COMPLEX * const result,
+                                   FNFT_INT * const W,
                                    fnft__akns_discretization_t discretization,
                                    fnft__akns_pde_t const PDE,
                                    FNFT_UINT const vanilla_flag,
@@ -115,12 +121,20 @@ FNFT_INT fnft__akns_scatter_matrix(FNFT_UINT const D,
  * The \f$b(\lambda)\f$ are calculated using the criterion from
  * Prins and Wahls, <a href="https://doi.org/10.1109/ACCESS.2019.2932256">&quot;
  * Soliton Phase Shift Calculation for the Korteweg–De Vries Equation,&quot;</a>.
+ * @param[in,out] Ws Pass an array of size K. Upon exit, it contains scaling factors that
+ * arise due to an internal normalization of the scattering process (to deal with potential
+ * overflow issues). The returned values for a and a_prime still have to be multiplied with
+ * corresponding power of two (i.e. POW(2, Ws[i])) to obtain the final values. Note
+ * that this is not required for the values of b. No normalization is carried out if NULL
+ * is passed.
  * @param[in] discretization The type of discretization to be used. Should be of type
  * \link fnft__akns_discretization_t \endlink. Not all akns_discretization_t discretizations are supported.
  * Check \link fnft_nse_discretization_t \endlink for list of supported types.
  * @param[in] PDE The partial differential equation for which the calculation
  * has to be done. Should be of type \link fnft__akns_pde_t \endlink.
- * @param[in] vanilla_flag For calculations for the KdV equation, pass 1 for the original mapping to the AKNS framework with r=-1. Pass 0 for the alternative mapping with q=-1. Unused for NSE.
+ * @param[in] vanilla_flag For calculations for the KdV equation, pass 1 for the
+ * original mapping to the AKNS framework with r=-1. Pass 0 for the alternative
+ * mapping with q=-1. Unused for NSE.
  * @param[in] skip_b_flag If set to 1 the routine will not compute \f$b(\lambda)\f$.
  * @return \link FNFT_SUCCESS \endlink or one of the FNFT_EC_... error codes
  *  defined in \link fnft_errwarn.h \endlink.
@@ -136,6 +150,7 @@ FNFT_INT akns_scatter_bound_states(FNFT_UINT const D,
                                    FNFT_COMPLEX * const a_vals,
                                    FNFT_COMPLEX * const aprime_vals,
                                    FNFT_COMPLEX * const b,
+                                   FNFT_INT * Ws,
                                    fnft__akns_discretization_t const discretization,
                                    fnft__akns_pde_t const PDE,
                                    FNFT_UINT const vanilla_flag,

include/private/fnft__errwarn.h

@@ -14,7 +14,7 @@
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 *
 * Contributors:
-* Sander Wahls (TU Delft) 2017-2018.
+* Sander Wahls (TU Delft) 2017-2018, 2023.
 */
 
 /**
@@ -100,6 +100,13 @@
  */
 #define FNFT__CHECK_RETCODE(ret_code, label) {if (ret_code!=FNFT_SUCCESS) {ret_code=FNFT__E_SUBROUTINE(ret_code);goto label;}}
 
+/**
+ * Macro for checking if memory has been allocated. If var!=NULL, a
+ * FNFT__E_NOMEM error is assigned to ret_code and the macro uses goto to jump to label.
+ * @ingroup private_errwarn
+ */
+#define FNFT__CHECK_NOMEM(var,ret_code,label) {if (var==NULL) {ret_code=FNFT__E_NOMEM;goto label;}}
+
 #ifdef FNFT_ENABLE_SHORT_NAMES
 #define WARN(msg)           FNFT__WARN(msg)
 #define E_NOMEM             FNFT__E_NOMEM
@@ -112,6 +119,7 @@
 #define E_SANITY_CHECK_FAILED(msg) FNFT__E_SANITY_CHECK_FAILED(msg)
 #define E_ASSERTION_FAILED  FNFT__E_ASSERTION_FAILED
 #define CHECK_RETCODE(ret_code,label) FNFT__CHECK_RETCODE(ret_code,label)
+#define CHECK_NOMEM(var,ret_code,label) FNFT__CHECK_NOMEM(var,ret_code,label)
 #endif
 
 /* --- Auxiliary macros and functions. Do not use directly. --- */

include/private/fnft__kdv_scatter.h

@@ -14,9 +14,10 @@
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 *
 * Contributors:
-* Sander Wahls (TU Delft) 2017-2018.
+* Sander Wahls (TU Delft) 2017-2018, 2023.
 * Shrinivas Chimmalgi (TU Delft) 2017.
 * Peter J Prins (TU Delft) 2020.
+* Sander Wahls (KIT) 2023.
 */
 
 /**
@@ -66,6 +67,12 @@
  * The \f$b(\lambda)\f$ are calculated using the criterion from
  * Prins and Wahls, <a href="https://doi.org/10.1109/ACCESS.2019.2932256">&quot;
  * Soliton Phase Shift Calculation for the Korteweg–De Vries Equation,&quot;</a>.
+ * @param[in,out] Ws Pass an array of size K. Upon exit, it contains scaling factors that
+ * arise due to an internal normalization of the scattering process (to deal with potential
+ * overflow issues). The returned values for a and a_prime still have to be multiplied with
+ * corresponding power of two (i.e. POW(2, Ws[i])) to obtain the final values. Note
+ * that this is not required for the values of b. No normalization is carried out if NULL
+ * is passed.
  * @param[in] discretization The type of discretization to be used. Should be of type
  * \link fnft_kdv_discretization_t \endlink. Not all kdv_discretization_t discretizations are supported.
  * Check \link fnft_kdv_discretization_t \endlink for list of supported types.
@@ -77,7 +84,7 @@
 FNFT_INT fnft__kdv_scatter_bound_states(const FNFT_UINT D, FNFT_COMPLEX const * const q,
     FNFT_COMPLEX const * const r, FNFT_REAL const * const T, FNFT_UINT const K,
     FNFT_COMPLEX * const bound_states, FNFT_COMPLEX * const a_vals,
-    FNFT_COMPLEX * const aprime_vals, FNFT_COMPLEX * const b,
+    FNFT_COMPLEX * const aprime_vals, FNFT_COMPLEX * const b, FNFT_INT * const Ws,
     fnft_kdv_discretization_t const discretization, FNFT_UINT const skip_b_flag);
 
 /**
@@ -106,6 +113,10 @@ FNFT_INT fnft__kdv_scatter_bound_states(const FNFT_UINT D, FNFT_COMPLEX const *
  * If derivative_flag=1 returns [S11 S12 S21 S22 S11' S12' S21' S22'] in
  * result where S11' is the derivative of S11 w.r.t to lambda.
  * Should be preallocated with size 4*K or 8*K accordingly.
+ * @param[in,out] W Pass an array of size K. Upon exit, it contains scaling factors that
+ * arise due to an internal normalization of the scattering matrix (to deal with potential
+ * overflow issues). The result has to be scaled by the power of 2 to obtain the final
+ * values. No normalization is carried out if NULL is passed.
  * @param[in] discretization The type of discretization to be used. Should be of type
  * \link fnft_kdv_discretization_t \endlink. Not all kdv_discretization_t discretizations are supported.
  * Check \link fnft_kdv_discretization_t \endlink for list of supported types.
@@ -116,10 +127,11 @@ FNFT_INT fnft__kdv_scatter_bound_states(const FNFT_UINT D, FNFT_COMPLEX const *
  *  defined in \link fnft_errwarn.h \endlink.
  * @ingroup kdv
  */
-FNFT_INT fnft__kdv_scatter_matrix(const FNFT_UINT D, FNFT_COMPLEX const * const q, FNFT_COMPLEX const * const r,
+FNFT_INT fnft__kdv_scatter_matrix(const FNFT_UINT D, 
+    FNFT_COMPLEX const * const q, FNFT_COMPLEX const * const r,
     const FNFT_REAL eps_t, const FNFT_INT kappa, const FNFT_UINT K,
-    FNFT_COMPLEX const * const lambda,
-    FNFT_COMPLEX * const result, fnft_kdv_discretization_t const discretization,
+    FNFT_COMPLEX const * const lambda, FNFT_COMPLEX * const result, 
+    FNFT_INT * const W, fnft_kdv_discretization_t const discretization,
     const FNFT_UINT derivative_flag);
 
 #ifdef FNFT_ENABLE_SHORT_NAMES