Merge pull request #77 from FastNFT/nsev_improvements

Nsev improvements

Author: wahls

CHANGELOG.md

@@ -4,16 +4,19 @@
 
 ### 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.
+- A NFT routine fnft_manakovv for the Manakov equation with vanishing boundaries was added (continuous spectrum only).
 - The slow scattering methods for AKNS-type systems now include a normalization procedure to deal with numerical overflow (enabled by default).
+- The periodic NFT routine fnft_nsep now also supports pure Newton refinement.
+- The vanishing NFT routine fnft_nsev now also supports manual filtering.
 - 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.
+- New criteria for stopping Newton iterations in fnft_nsep and fnft_nsev.
+- The tolerance for the Newton refinements in fnft_nsev can now be set by the user.
+- The default number of iterations for the Newton refinements in fnft_nsev is now 100, and the user is warned if the number was too small.
+- Reduced some tests to reduce run times.
 - The code for AKNS scattering has been overhauled.
 
 ### Fixed

Doxyfile

@@ -77,7 +77,7 @@ Please use the [issue tracker](https://github.com/FastNFT/FNFT/issues) to report
 
 ## Contributors
 
-* Sander Wahls, TU Delft
+* Sander Wahls, KIT (since July 2023) and TU Delft (before)
 * Shrinivas Chimmalgi, TU Delft
 * Peter J. Prins, TU Delft
 * Marius Brehler, TU Dortmund

README.md

@@ -53,7 +53,7 @@
  *  initial guesses for the bound states are available. The complexity is
  *  \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$,
+ *  \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$ 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

include/fnft_kdvv.h

@@ -17,6 +17,7 @@
 * Sander Wahls (TU Delft) 2017-2018.
 * Shrinivas Chimmalgi (TU Delft) 2019-2020.
 * Peter J Prins (TU Delft) 2020-2021.
+* Sander Wahls (KIT) 2023.
 */
 
 /**
@@ -47,12 +48,15 @@
  *  to each other are merged. \n \n
  *  fnft_nsev_bsfilt_FULL: Bound states in physically implausible regions and
  *  outside the region based on the step-size of the supplied samples are
- *  rejected.
+ *  rejected. \n \n
+ *  fnft_nsev_opts_filt_MANUAL: Only points within the specified
+ *  \link fnft_nsev_opts_t::bounding_box \endlink are kept. \n\n
  */
 typedef enum {
     fnft_nsev_bsfilt_NONE,
     fnft_nsev_bsfilt_BASIC,
-    fnft_nsev_bsfilt_FULL
+    fnft_nsev_bsfilt_FULL,
+    fnft_nsev_bsfilt_MANUAL
 } fnft_nsev_bsfilt_t;
 
 /**
@@ -63,8 +67,7 @@ typedef enum {
  *  https://arxiv.org/abs/1611.02435 and https://github.com/eiscor/eiscor)
  *  with \f$ O(D^2) \f$ complexity is used to detect the roots of
  *  \f$ a(\lambda) \f$. (Note: FNFT incorporates a development version of this
- *  routine as no release was available yet.) This method is relatively slow,
- *  but very reliable. \n \n
+ *  routine as no release was available yet.) This method is relatively slow. \n \n
  *  fnft_nsev_bsloc_NEWTON: Newton's method is used to refine a given set of initial guesses.
  *  The discretization used for the the refinement is one of the base methods \link fnft_nse_discretization_t.h \endlink.
  *  The number of iterations is specified through the field \link fnft_nsev_opts_t::niter
@@ -76,7 +79,7 @@ typedef enum {
  *  initial guesses for the bound states are available. The complexity is
  *  \f$ O(niter (*K\_ptr) D) \f$. \n \n
  *  fnft_nsev_bsloc_SUBSAMPLE_AND_REFINE: This method offers a good compromise
- *  between the other two. The method automatically finds initial guesses for
+ *  between the previous two. The method automatically finds initial guesses for
  *  the NEWTON method by first applying the FAST_EIGENVALUE method to a
  *  subsampled version of the signal. Second these initial guesses are refined
  *  using the NEWTON method. The number of samples of the subsampled signal can
@@ -158,9 +161,13 @@ typedef enum {
  *   \link fnft_nsev_bsloc_t \endlink for details.
  *
  * @var fnft_nsev_opts_t::niter
- *  Number of Newton iterations to be carried out when either the
- *  fnft_nsev_bsloc_NEWTON or the fnft_nsev_bsloc_SUBSAMPLE_AND_REFINE method
- *  is used.
+ *  For fnft_nsev_bsloc_NEWTON and fnft_nsev_bsloc_SUBSAMPLE_AND_REFINE: Maximum
+ *  number of Newton iterations to be carried out.
+ * 
+ *  @var fnft_nsev_opts_t::tol
+ *  Some bound state localization methods such as Newton have tolerance
+ *  parameters. If this value is negative, these parameters will be chosen
+ *  automatically. Set to a non-negative value to chose the tolerance manually.
  *
  * @var fnft_nsev_opts_t::discspec_type
  *  Controls how \link fnft_nsev \endlink fills the array
@@ -194,18 +201,26 @@ typedef enum {
  *  such as discontinuous signals, applying Richardson extrapolation may result in
  *  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.
+ *  flag is zero). To enable, set the flag to one.\n\n
+ *
+ * @var fnft_nsev_opts_t::bounding_box
+ *  Array of four reals. Defines a box in the complex plane that is used for
+ *  manual filtering: \n
+ *  bounding_box[0] <= real(lambda) <= bounding_box[1] \n
+ *  bounding_box[2] <= imag(lambda) <= bounding_box[3] \n
  */
 typedef struct {
     fnft_nsev_bsfilt_t bound_state_filtering;
     fnft_nsev_bsloc_t bound_state_localization;
     FNFT_UINT niter;
+    FNFT_REAL tol;
     FNFT_UINT Dsub;
     fnft_nsev_dstype_t discspec_type;
     fnft_nsev_cstype_t contspec_type;
     FNFT_INT normalization_flag;
     fnft_nse_discretization_t discretization;
     FNFT_UINT richardson_extrapolation_flag;
+    FNFT_REAL bounding_box[4];
 } fnft_nsev_opts_t;
 
 /**
@@ -215,14 +230,16 @@ typedef struct {
  * @returns A \link fnft_nsev_opts_t \endlink object with the following options.\n
  *  bound_state_filtering = fnft_nsev_bsfilt_FULL\n
  *  bound_state_localization = fnft_nsev_bsloc_SUBSAMPLE_AND_REFINE\n
- *  niter = 10\n
+ *  niter = 100\n
+ *  tol = -1.0\n
  *  discspec_type = fnft_nsev_dstype_NORMING_CONSTANTS\n
  *  contspec_type = fnft_nsev_cstype_REFLECTION_COEFFICIENT\n
  *  normalization_flag = 1\n
  *  discretization = fnft_nse_discretization_2SPLIT4B\n
  *  richardson_extrapolation_flag = 0\n
+ *  bounding_box = {NAN, NAN, NAN, NAN}\n
  *
-  * @ingroup fnft
+ * @ingroup fnft
  */
 fnft_nsev_opts_t fnft_nsev_default_opts();
 
@@ -377,6 +394,7 @@ FNFT_INT fnft_nsev(const FNFT_UINT D, FNFT_COMPLEX const * const q,
 #define nsev_bsfilt_NONE fnft_nsev_bsfilt_NONE
 #define nsev_bsfilt_BASIC fnft_nsev_bsfilt_BASIC
 #define nsev_bsfilt_FULL fnft_nsev_bsfilt_FULL
+#define nsev_bsfilt_MANUAL fnft_nsev_bsfilt_MANUAL
 #define nsev_bsloc_FAST_EIGENVALUE fnft_nsev_bsloc_FAST_EIGENVALUE
 #define nsev_bsloc_NEWTON fnft_nsev_bsloc_NEWTON
 #define nsev_bsloc_SUBSAMPLE_AND_REFINE fnft_nsev_bsloc_SUBSAMPLE_AND_REFINE

include/fnft_nsev.h

@@ -117,7 +117,7 @@ FNFT_INT fnft__akns_scatter_matrix(FNFT_UINT const D,
  * @param[out] a_vals Array of length K, contains the values of \f$a(\lambda)\f$.
  * @param[out] aprime_vals Array of length K, contains the values of
  * \f$ a'(\lambda) = \frac{\partial a(\lambda)}{\partial \lambda}\f$.
- * @param[out] b Array of length K, contains the values of \f$b(\lambda)\f$.
+ * @param[out] b_vals Array of length K, contains the values of \f$b(\lambda)\f$.
  * 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>.
@@ -149,7 +149,7 @@ FNFT_INT akns_scatter_bound_states(FNFT_UINT const D,
                                    FNFT_COMPLEX const * const bound_states,
                                    FNFT_COMPLEX * const a_vals,
                                    FNFT_COMPLEX * const aprime_vals,
-                                   FNFT_COMPLEX * const b,
+                                   FNFT_COMPLEX * const b_vals,
                                    FNFT_INT * Ws,
                                    fnft__akns_discretization_t const discretization,
                                    fnft__akns_pde_t const PDE,

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, 2022.
 * Shrinivas Chimmalgi (TU Delft) 2019-2020.
 * Peter J. Prins (2021).
+* Sander Wahls (KIT) 2023.
 */
 
 #include <string.h>
@@ -164,7 +165,21 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 
             /* Increase k to account for vector of initial guesses */
     	    k++;
-            
+
+        } else if ( strcmp(str, "bsloc_tol") == 0 ) {
+
+            /* Extract desired number of iterations */
+            if ( k+1 == nrhs || !mxIsDouble(prhs[k+1])
+            || mxGetNumberOfElements(prhs[k+1]) != 1
+                    || mxGetScalar(prhs[k+1]) < 0.0 ) {
+                snprintf(msg, sizeof msg, "'bsloc_tol' should be followed by a non-negative real scalar.");
+                goto on_error;
+            }
+            opts.tol = (FNFT_REAL)mxGetScalar(prhs[k+1]);
+
+            /* Increase k to account for vector of initial guesses */
+    	    k++;
+
         } else if ( strcmp(str, "bsloc_subsamp_refine") == 0 ) {
 
             opts.bound_state_localization = fnft_nsev_bsloc_SUBSAMPLE_AND_REFINE;
@@ -195,6 +210,26 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
         } else if ( strcmp(str, "bsfilt_full") == 0 ) {
 
             opts.bound_state_filtering = fnft_nsev_bsfilt_FULL;
+
+        } else if ( strcmp(str, "bsfilt_manual") == 0 ) {
+
+            opts.bound_state_filtering = fnft_nsev_bsfilt_MANUAL;
+
+            /* Extract bounding box for manual filltering */
+            if ( k+1 == nrhs || !mxIsDouble(prhs[k+1])
+            || mxGetM(prhs[k+1]) != 1 || mxGetN(prhs[k+1]) != 4) {
+                snprintf(msg, sizeof msg, "'filt_manual' should be followed by a real row vector of length four. See the help.");
+                goto on_error;
+            }
+            double const * const tmp = mxGetPr(prhs[k+1]);
+            opts.bounding_box[0] = (FNFT_REAL)tmp[0];
+            opts.bounding_box[1] = (FNFT_REAL)tmp[1];
+            opts.bounding_box[2] = (FNFT_REAL)tmp[2];
+            opts.bounding_box[3] = (FNFT_REAL)tmp[3];
+
+            /* Increase k to account for bounding box vector */
+            k++;
+
         } else if ( strcmp(str, "RE") == 0 ) {
 
             opts.richardson_extrapolation_flag  = 1;
@@ -314,7 +349,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 
             opts.discretization = fnft_nse_discretization_4SPLIT4B;
 
-
         // Slow discretizations
         } else if ( strcmp(str, "discr_BO") == 0 ) {
 
@@ -374,8 +408,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
         if (bound_states == NULL) {
             K = fnft_nsev_max_K(D, &opts);
             if (K == 0) {
-                snprintf(msg, sizeof msg, "Discretization does not support the chosen bound state localization method.");
-                goto on_error;
+                //snprintf(msg, sizeof msg, "Discretization does not support the chosen bound state localization method.");
+                //goto on_error;
+                K = D;
             }
             bound_states = mxMalloc(K * sizeof(FNFT_COMPLEX));
         }

matlab/mex_fnft_nsev.c

@@ -36,7 +36,9 @@
 %                   below) is set by the algorithm. Not followed by a
 %                   value.
 %   'bsloc_niter'   Number of iterations to be carried by Newton's method.
-%                   Followed by a positive integer.
+%                   Followed by a non-negative integer.
+%   'bsloc_tol'     Tolerance used for stopping Newton's method. Followed
+%                   by a non-negative real number.
 %   'bsloc_Dsub'    The desired number of samples for the subsampled signal
 %                   in the 'subsamp_refine' method. Less samples in the
 %                   subsampled stage result in faster execution time, but
@@ -48,6 +50,10 @@
 %                   bound states in the lower half plane.
 %   'bsfilt_full'   Full bound state filtering. Additionally removes
 %		            physically implausible bound states.
+%   'bsfilt_manual' Manual filtering. Remove bound states outside of a box
+%                   provided by the user. Must be followed by a real vector
+%                   of the form [min_real, max_real, min_imag, max_imag]
+%                   that specifies the box.
 %   'RE'            Use Richardson extrpolation to improve accuracy. The
 %                   approximations of the nonlinear Fourier spectrum are
 %                   calcuated using all given samples and again with half of
@@ -133,6 +139,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, 2022.
 % Shrinivas Chimmalgi (TU Delft) 2019-2020.
 % Peter J. Prins (TU Delft) 2021.
+% Sander Wahls (KIT) 2023.