Merge pull request #73 from FastNFT/fnft_nsep-improvements

fnft_nsep improvements

Author: wahls

examples/mex_fnft_nsep_example.m

@@ -13,7 +13,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, 2021.
 % Shrinivas Chimmalgi (TU Delft) 2020.
 
 % This examples demonstrates how the nonlinear Fourier transform with
@@ -52,7 +52,9 @@
 % interval [-3j, 3j]. Furthermore, there are degenerate spines
 % of length zero at the degenerate points of the main spectrum.
 
-spines = mex_fnft_nsep(q, T, kappa, 'phase_shift', phase_shift, 'points_per_spine', 100);
+filt_box = [-10 10 -10 10]; % manual filtering provides a speed up
+spines = mex_fnft_nsep(q, T, kappa, 'phase_shift', phase_shift, ...
+    'points_per_spine', 100, 'filt_manual', filt_box, 'loc_max_evals', 100);
 
 % Increase the number of points per spine above to improve the
 % resolution of the spine (at the cost of increased run times).

include/fnft_nsep.h

@@ -14,7 +14,7 @@
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 *
 * Contributors:
-* Sander Wahls (TU Delft) 2017-2018, 2020.
+* Sander Wahls (TU Delft) 2017-2018, 2020-21.
 * Shrinivas Chimmalgi (TU Delft) 2020.
 */
 
@@ -43,16 +43,21 @@
  *  fnft_nsep_opts_loc_SUBSAMPLE_AND_REFINE: Similar approach as for
  *  fnft_nsev_opts_dsloc_SUBSAMPLE_AND_REFINE (see
  *  \link fnft_nsev_opts_t::bound_state_localization \endlink.)\n\n
+ *  fnft_nsep_opts_loc_NEWTON: Refine initial guesses for the main and/or
+ *  auxiliary spectrum by the user using Newton's method.\n\n
+ *  fnft_nsev_opts_dsloc_SUBSAMPLE_AND_REFINE (see
+ *  \link fnft_nsev_opts_t::bound_state_localization \endlink.)\n\n
  *  fnft_nsep_opts_loc_GRIDSEARCH: Uses a grid search to localize roots. Can
  *  only find main and auxiliary spectrum points on the real axis. In the
- *  defocusing case, the main spectrum is always real. The implemented grid 
+ *  defocusing case, the main spectrum is always real. The implemented grid
  *  search gurantees only linear convergence.\n\n
  *  fnft_nsep_opts_loc_MIXED: Uses the SUBSAMPLE_AND_REFINE method to find the
  *  non-real parts of the spectra and the GRIDSEARCH method to find the real
  *  parts.
  */
 typedef enum {
     fnft_nsep_loc_SUBSAMPLE_AND_REFINE,
+    fnft_nsep_loc_NEWTON,
     fnft_nsep_loc_GRIDSEARCH,
     fnft_nsep_loc_MIXED
 } fnft_nsep_loc_t;
@@ -211,6 +216,10 @@ fnft_nsep_opts_t fnft_nsep_default_opts();
  *       - fnft_nse_discretization_2SPLIT8B
  *       - fnft_nse_discretization_4SPLIT4A
  *
+ * For the Newton refinement mode only, one of the following should be used instead:
+ *      - fnft_nse_discretization_BO
+ *      - fnft_nse_discretization_CF4_2
+ *
  * @param[in] D Number of samples. Should be power of two and \f$ D\geq 2\f$.
  * @param[in] q Array of length D, contains samples \f$ q(t_n)=q(x_0, t_n) \f$,
  *  where \f$ t_n = T[0] + nL/D \f$, where \f$L=T[1]-T[0]\f$ is the period and
@@ -222,18 +231,25 @@ fnft_nsep_opts_t fnft_nsep_default_opts();
  * @param[in] phase_shift Real scalar constant. It is the change in the phase
  * over one quasi-period,\f$ \angle(q(t+L)/q(t))\f$. For periodic signals it will be 0.
  * @param[in,out] K_ptr Upon entry, *K_ptr should contain the length of the array
- *  main_spec. Upon return, *K_ptr contains the number of actually detected
+ *  main_spec (except for Newton localization, see below). Upon return,
+ *  *K_ptr contains the number of actually detected
  *  points in the main spectrum. If the length of the array main_spec was not
  *  sufficient to store all of the detected points in the main spectrum, a
  *  warning is printed and as many points as possible are returned instead.
  *  Note that in order to skip the computation of the main spectrum completely,
  *  it is not sufficient to pass *K_ptr==NULL. Instead, pass main_spec==NULL.
+ *  EXCEPTION: If opts->localization==fnft_nsep_loc_NEWTON, then *K_ptr should
+ *  contain the number of initial guess per spine point (see opts->points_per_spine).
  * @param[out] main_spec Array. Upon return, the routine has stored the detected
  *  main specrum points (i.e., the points for which the trace of the monodromy
  *  matrix is either +2 or -2) in the first *K_ptr entries of this array.
  *  If NULL is passed instead, the main spectrum will not be computed.
  *  Has to be preallocated by the user. The user can choose an arbitrary
  *  length. Typically, D is a good choice.
+ *  If opts->localization==fnft_nsep_loc_NEWTON, then main_spec should be of
+ *  the form [g1_1, ..., gK_1, g1_2, ..., gK_2, ..., g1_P, ..., gK_P] when
+ *  the routine is called, where gk_n is the k-th initial guess for the n-th
+ *  spine point, K=*K_ptr, and P=opts->points_per_spine.
  * @param[in,out] M_ptr Upon entry, *M_ptr should contain the length of the array
  *  aux_spec. Upon return, *M_ptr contains the number of actually detected
  *  points in the auxiliary spectrum. If the length of the array aux_spec was

matlab/mex_fnft_nsep.c

@@ -14,7 +14,7 @@
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 *
 * Contributors:
-* Sander Wahls (TU Delft) 2017-2018, 2020.
+* Sander Wahls (TU Delft) 2017-2018, 2020-2021.
 * Shrinivas Chimmalgi (TU Delft) 2020.
 */
 
@@ -30,17 +30,17 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
     double * T;
     double phase_shift = 0.0;
     size_t M;
-    double complex * main_spec;
+    double complex * main_spec = NULL;
     size_t K;
     double complex * aux_spec = NULL;
     int * sheet_indices = NULL;
     int kappa;
     size_t upsampling_factor = 1;
-    size_t i;
+    size_t i, j;
     ptrdiff_t k;
     double *re, *im;
     double *msr, *msi, *asr, *asi;
-    char msg[128]; // buffer for error messages
+    char msg[256]; // buffer for error messages
     fnft_nsep_opts_t opts;
     int ret_code;
 
@@ -145,6 +145,51 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 
             opts.localization = fnft_nsep_loc_SUBSAMPLE_AND_REFINE;
 
+        } else if ( strcmp(str, "loc_newton") == 0 ) {
+
+            opts.localization = fnft_nsep_loc_NEWTON;
+
+            /* Extract initial guesses for the main spectrum */
+            if (k+1 == nrhs || 
+                    (!mxIsEmpty(prhs[k+1]) &&
+                        (!mxIsComplex(prhs[k+1]) || mxGetM(prhs[k+1]) != 1))) {
+                snprintf(msg, sizeof msg, "'loc_newton' should be followed by two complex row vectors of initial guesses for the main and auxiliary spectrum, respectively. Try passing complex(...).");
+                goto on_error;
+            }
+            K = mxGetN(prhs[k+1]);
+            main_spec = mxMalloc(K * sizeof(FNFT_COMPLEX));
+            if (K > 0 && main_spec == NULL) {
+                snprintf(msg, sizeof msg, "Out of memory.");
+                goto on_error;
+            }
+            re = mxGetPr(prhs[k+1]);
+            im = mxGetPi(prhs[k+1]);
+            for (j=0; j<K; j++)
+                main_spec[j] = re[j] + I*im[j];
+
+            /* Extract initial guesses for the aux spectrum */
+            if (k+2 == nrhs || 
+                    (!mxIsEmpty(prhs[k+2]) &&
+                        (!mxIsComplex(prhs[k+2]) || mxGetM(prhs[k+2]) != 1))) {
+                snprintf(msg, sizeof msg, "'loc_newton' should be followed by two complex row vectors of initial guesses for the main and auxiliary spectrum, respectively. Try passing complex(...).");
+                goto on_error;
+            }
+            M = mxGetN(prhs[k+2]);
+            aux_spec = mxMalloc(M * sizeof(FNFT_COMPLEX));
+            if (M > 0 && aux_spec == NULL) {
+                snprintf(msg, sizeof msg, "Out of memory.");
+                goto on_error;
+            }
+            re = mxGetPr(prhs[k+2]);
+            im = mxGetPi(prhs[k+2]);
+            for (j=0; j<M; j++)
+                aux_spec[j] = re[j] + I*im[j];
+
+            /* Increase k to account for the two vectors of initial guesses */
+            k += 2;
+
+            opts.discretization = fnft_nse_discretization_BO;
+
         } else if ( strcmp(str, "loc_gridsearch") == 0 ) {
 
             opts.localization = fnft_nsep_loc_GRIDSEARCH;
@@ -204,13 +249,24 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 
     /* Allocate memory */
     q = mxMalloc(D * sizeof(double complex));
-    K = (opts.points_per_spine)*upsampling_factor*D + 1;
-    M = upsampling_factor*D;
-    main_spec = mxMalloc(K * sizeof(double complex));
-    aux_spec = mxMalloc(M * sizeof(double complex));
-    sheet_indices = mxMalloc(M * sizeof(int));
-    if ( q == NULL || main_spec == NULL || aux_spec == NULL
-    || sheet_indices == NULL ) {
+    if (M > 0)
+        sheet_indices = mxMalloc(M * sizeof(int));
+    if (opts.localization != fnft_nsep_loc_NEWTON) { /* Not Newton */
+        K = (opts.points_per_spine)*upsampling_factor*D + 1;
+        M = upsampling_factor*D;
+        main_spec = mxMalloc(K * sizeof(double complex));
+        aux_spec = mxMalloc(M * sizeof(double complex));
+    } else { /* Newton */
+        if (K%opts.points_per_spine != 0) {
+            snprintf(msg, sizeof msg, "The lengths of the first initial guess vector for Newton localization has to be a multiple of points_per_spine (=%u).",
+                (unsigned int)opts.points_per_spine);
+            goto on_error;
+        }
+        K /= opts.points_per_spine;
+    }
+
+    if ( q == NULL || (K>0 && main_spec == NULL) || (M > 0 && aux_spec == NULL)
+         || (M > 0 && sheet_indices == NULL) ) {
         snprintf(msg, sizeof msg, "Out of memory.");
         goto on_error;
     }

matlab/mex_fnft_nsep.m

@@ -18,12 +18,24 @@
 %   It is possible to provide additional inputs. These come either in the
 %   form of a single string or of a string followed by a value.
 %
-%    'loc_subsample_and_refine' Localize spectra by applying fast eigenmethod
+%    'loc_subsample_and_refine' Localize spectra by applying a fast eigenmethod
 %                   to a subsampled version of the signal, then refine using
-%                   the complete signal.
+%                   the complete signal (default).
 %    'loc_gridsearch' Localize real spectra by using a FFT-based grid search.
 %    'loc_mixed'    Localize real spectra using grid search and non-real using
-%                   subsample and refine (default).
+%                   subsample and refine.
+%    'loc_newton'   Localize spectra by refining initial guesses using Newton's
+%                   method (currently, the BO discretization is always used).
+%                   This argument must be followed by two complex row vectors
+%                   of initial guesses for the main and auxiliary spectrum,
+%                   respectively. (Alternatively, it is also possible to pass
+%                   [] if the corresponding part of the spectrum should be
+%                   skipped.) The length of the first vector must be a
+%                   multiple of points_per_spine (default two, see below).
+%                   It should be of the form
+%                     [g1_1,...,gK_1, g1_2, ..., gK_2, ..., g1_P, ..., gK_P],
+%                   where P=points_per_spine and gk_n denotes the k-th initial
+%                   guess for the n-th spine point.
 %    'loc_max_evals' Maximum number of monodromy matrix evaluations after which
 %                   refinement is stopped. (The acutual number might be slightly
 %                   higher.) This argument must be followed by a real
@@ -49,11 +61,11 @@
 %                   quasi-period. This argument must be followed by a real
 %                   scalar i.e. angle(q(T(2))/q(T(1))). The default value 
 %                   is 0.0, which assumes that q is periodic.
-%   'discr_modal'   Use modified Ablowitz-Ladik discretization.
-%   'discr_2split2A' Use split Boffetta-Osborne discretization (default).
-%   'discr_2split4A' Use fifth order splitting with 4th degree polynomial.
-%   'discr_2split4B' Use fifth order splitting with 2nd degree polynomial.
-%   'discr_4split4B' Fourth-order method. Uses fifth order splitting with
+%    'discr_modal'   Use modified Ablowitz-Ladik discretization.
+%    'discr_2split2A' Use split Boffetta-Osborne discretization (default).
+%    'discr_2split4A' Use fifth order splitting with 4th degree polynomial.
+%    'discr_2split4B' Use fifth order splitting with 2nd degree polynomial.
+%    'discr_4split4B' Fourth-order method. Uses fifth order splitting with
 %                    2nd degree polynomial.
 %    'quiet'        Turns off messages generated by then FNFT C library.
 %                   (To turn off the messages generated by the mex
@@ -79,6 +91,6 @@
 % along with this program. If not, see <http://www.gnu.org/licenses/>.
 %
 % Contributors:
-% Sander Wahls (TU Delft) 2018, 2020.
+% Sander Wahls (TU Delft) 2018, 2020-2021.
 % Shrinivas Chimmalgi (TU Delft) 2020.