If both sdpt3 and sedumi fail, should we

• Abort
• Try a new random point in the bisection interval
• Mark the point infeasible and bisect, but return with a warning

Since xUnit is not avalable anymore, we should switch to the standard unit test framework of MATLAB. I can give it a shot in the next few days.

Consider reformatting the LICENSE file so GitHub can automatically identify it as the BSD 3-clause license with licensee.
JOSS Review 2514

I think all the functionality from the low-storage and ssp directories has been incorporated into 'general'. So we should delete the obsolete directories and just move 'general' up into the root directory.

JOSS Review 2514

In addition to the Contributing section in the documentation it would be good to have a top-level CONTRIBUTING file as well so users see a link to it when they open a PR. It would also be good to note if you prefer users create a GitHub issue or use some other means for reporting problems.

To control what information is written to the screen.

Currently, the option solveorderconditions cannot be used for low-storage schemes:

>> rk=rk_opt(5, 2, '3Sstar', 'acc', 'num_starting_points', 2, 'solveorderconditions', 1, 'algorithm', 'interior-point')
Error using  * 
Incorrect dimensions for matrix multiplication. Check that the number of columns in the
first matrix matches the number of rows in the second matrix. To perform elementwise
multiplication, use '.*'.

Error in order_conditions (line 16)
lin_tau = Aeq*x'-beq';

Error in rk_opt>@(x)order_conditions(x,class,s,p,Aeq,beq) (line 113)
        x(i,:) = fsolve(@(x) order_conditions(x,class,s,p,Aeq,beq), x(i,:));

Error in fsolve (line 242)
            fuser = feval(funfcn{3},x,varargin{:});

Error in rk_opt (line 113)
        x(i,:) = fsolve(@(x) order_conditions(x,class,s,p,Aeq,beq), x(i,:));

Caused by:
    Failure in initial objective function evaluation. FSOLVE cannot continue.

While the error message clearly points to the point of failure and suggests the workaround to not use solveorderconditions for low-storage methods, it is still nice to document this problem (done via this issue) and probably also to fix it if needed.

To cite RK-Opt, it might be nice to have a DOI from https://zenodo.org/. If you want, I can set up a .zenodo.json file.

JOSS Review 2514

I think the documentation looks excellent. The only overarching comment I have is that the function documentation style is a little inconsistent. Some functions have a brief description and a bulleted list of inputs/outputs with a description for each item while others just have the brief description which may not cover all the inputs/outputs for the function.

Below is a list of minor comments and typos or style inconsistencies I noticed in reading through the docs. I've linked to the corresponding line in the function comments since there does not seem to be an easy way to do that for the rst files.

• The style for variable names in the input/output lists is some times italic and some times not inconsistent e.g., rk_obj.
• k is missing from the list of inputs in the rst docs.
  

  RK-Opt/RK-coeff-opt/nonlinear_constraints.m

  Line 15 in bbf5649

  % * :math:`k`: Number of steps for multi-step, mlti-stage schemes.

• "Methods" to "methods."

  RK-Opt/RK-coeff-opt/oc_albrecht.m

  Line 4 in bbf5649

  % Order conditions for SSP RK Methods.

• Changing "RKMs" to "RK methods" would be clearer

  RK-Opt/RK-coeff-opt/oc_butcher.m

  Line 4 in bbf5649

  % Order conditions for RKMs.

• The meaning of solveoderconditions=1 is not clear here it would be good to have a note to see the documentation for rk_opt
• "initialize" to "initialized"

  RK-Opt/RK-coeff-opt/rk_opt.m

  Line 25 in bbf5649

  % * startvec: vector of the initial guess ('random' = random approach; 'smart' = smart approach; alternatively, the user can provide the startvec array. By default startvec is initialize with random numbers.

• "controll" to "control"

  RK-Opt/RK-coeff-opt/rk_opt.m

  Line 39 in bbf5649

  % * problem_class: class of problems for which the RK is designed ('linear' or 'nonlinear' problems). This option changes the type of order conditions check, i.e. linear or nonlinear order conditions controll. The default value is 'nonlinear'.

• In the rst docs, "tored" to "stored"

  RK-Opt/RK-coeff-opt/unpack_rk.m

  Line 6 in bbf5649

  % The coefficients are stored in a single vector x as::

• writeField does not have a description

  RK-Opt/RK-coeff-opt/write_field.m

  Line 2 in bbf5649

  % function writeField(writeFid,name,value)

• ". different" to "for different"

  RK-Opt/am_radius-opt/multi_R_opt.m

  Line 6 in bbf5649

  % of optimal values like those that appear in [ketcheson2009]_.

• In the rst docs this line and the next run together

  RK-Opt/am_radius-opt/multi_R_opt.m

  Line 9 in bbf5649

  % k = [k1, k2, ..., kK]^T, K = length(k), ith-element = # of steps

• Not sure if this indentation was intentional

  RK-Opt/am_radius-opt/multi_R_opt.m

  Line 21 in bbf5649

  % integration scheme, one is usually interested in understanding the

• Note the inputs p and q are arrays of polynomial coefficients for the numerator and denominator respectively

  RK-Opt/am_radius-opt/radimpfast.m

  Line 2 in bbf5649

  % function rad=radimpfast(p,q)

• R is missing from outputs in Rkp, Rkp_dw, Rkp_imp, and Rkp_imp_dw e.g.,

  RK-Opt/am_radius-opt/Rkp.m

  Line 10 in bbf5649

  % Outputs:

• I think this was meant to be a bulleted list like Rkp and Rkp_imp

  RK-Opt/am_radius-opt/Rkp_dw.m

  Line 7 in bbf5649

  % Inputs: k = # of steps

• missing ) after first F, missing backslash on second beta, and I believe t is intended to be \tilde{}

  RK-Opt/am_radius-opt/Rkp_dw.m

  Line 13 in bbf5649

  % `u_n = \sum_{j=0}^{k-1} (\alpha[j] + \beta[j] F(u_{n-k+j} + tbeta[j] tF(u_{n-k+j}))`

• tbeta is missing from the list of outputs

  RK-Opt/am_radius-opt/Rkp_imp_dw.m

  Line 2 in bbf5649

  % function [R,alpha,beta]=Rkp_imp_dw(k,p)

• I think this was meant to be a bulleted list like Rkp and Rkp_imp

  RK-Opt/am_radius-opt/Rkp_imp_dw.m

  Line 7 in bbf5649

  % Inputs: k = # of steps

• Fill in the general stability function

  RK-Opt/am_radius-opt/Rskp.m

  Line 19 in bbf5649

  % (Fill in)

• In the rst docs "epends" to "Depends"

  RK-Opt/am_radius-opt/Rskp.m

  Line 21 in bbf5649

  %Depends on MATLAB's optimization toolbox for the LP solver

• I really like the example calls included in the opt_poly_bisect documentation especially the second example that refers to a figure in a publication. Where relevant/applicable, it would be great to have similar examples in the docs for other functions.
• It would be nice if the references in the spectrum documentation linked to the articles (or a full citation).

  RK-Opt/polyopt/spectrum.m

  Line 11 in bbf5649

  % * 'Niegemann-ellipse' and 'Niegemann-circle': See Niegemann 2011

• Add "c should be a column vector of length m"

  RK-Opt/RKtools/am_radius.m

  Line 8 in bbf5649

  % and `b` should be a column vector of length m.

• "10^-{10})" to "10^{-10})."

  RK-Opt/RKtools/am_radius.m

  Line 10 in bbf5649

  % Accuracy can be changed by modifying the value of eps (default `10^-{10}`)

• In the rst docs the equations after "such that" are missing

  RK-Opt/RKtools/am_radius.m

  Line 15 in bbf5649

  % such that

• rk_stabfun uses rk.A and rk.b and plotstabreg accesses rk.Ahat, rk.bhat, and rk.chat. It would be good to document what fields can/should be in an rk variable.
• I noticed the contents of the dwrk_opt and low_storage directories are not documented are these pieces of RK-Opt that are in development?

Matrix dimensions do not agree

Here are some minor suggestions for the documentation.

• In https://rk-opt.readthedocs.io/en/latest/started.html section "Obtaining RK-Opt":
  In the line "add the subdirectory RK-Opt/RKtools", it would be nice to mention how it is done (addpath).

• In https://rk-opt.readthedocs.io/en/latest/RK-coeff-opt.html Line 5:
  Extra tick around the code snippet results_rkopt = runtests('test_rkopt.m'); table(results_rkopt).

• Same comment for test code snippet in https://rk-opt.readthedocs.io/en/latest/polyopt.html

• In https://rk-opt.readthedocs.io/en/latest/RK-coeff-opt.html section "rk_opt" in the fist "Note:" block:
  Capitalize the first letter. Also it would be nice add a reference pointing toward the said numerical experiment.

• In https://rk-opt.readthedocs.io/en/latest/RK-coeff-opt.html
  Please follow a single convention to specify the input and output variables associated with all the functions. For example, the description of the input variables for the functions oc_butcher, unpack_lsrk are scattered across the page. It would be helpful for the reader to find those descriptions near the function description, even if they are redundant. Maybe stick to the style of categorizing "Input" and "Output" variables followed in https://rk-opt.readthedocs.io/en/latest/am_radius-opt.html

• Same comment for various function in https://rk-opt.readthedocs.io/en/latest/RKtools.html and some of the function in https://rk-opt.readthedocs.io/en/latest/am_radius-opt.html

• In https://rk-opt.readthedocs.io/en/latest/RK-coeff-opt.html section "oc_ksrk", the warning could be in a block instead of ...warning.

• In https://rk-opt.readthedocs.io/en/latest/polyopt.html section "opt_poly_bisect":
  In "solvers:" spelling for "licenses"
  In "lam_func:" the line "see Figure 10 of the CAMCoS paper": should add direct link to reference.

• In https://rk-opt.readthedocs.io/en/latest/RKtools.html, last line of math mode:
  A period to complete the sentence.

JOSS Review 2514

Ideally, we could find some way to merge the two unpack_*.m functions.

It is now possible to call fmincon from Python's OpenOpt:

http://openopt.org/fmincon

We can use this to make a Python version of RK-coeff-opt, which will interface much more closely with NodePy.

JOSS Review 2514

In multi_am_radius_2_rk_opt.m there is a missing fgets that causes an error when using the example input file given in the comments (or the empty line before the first separation line in the example input needs to be removed).

--- a/examples/multi_am_radius_2_rk_opt.m
+++ b/examples/multi_am_radius_2_rk_opt.m
@@ -157,6 +156,9 @@ function [s,k,p,coeffs_monomial]= read_data(read_fid,class)
 %
 % Read number of stage and order of the method
 
+    % Read white line
+    tline = fgets(read_fid);
+
     % Read separation line
     tline = fgets(read_fid);    

There also a missing ; earlier in the file.

--- a/examples/multi_am_radius_2_rk_opt.m
+++ b/examples/multi_am_radius_2_rk_opt.m
@@ -64,7 +64,7 @@ for i_stab_poly = 1:number_stab_poly
         % The order of accuracy is set to one to re-use completely the 
         % functions implemented in the RK-coeff-opt to impose just the 
         % linear order conditions.
-        h_p = 1
+        h_p = 1;
         
         % Polynomial coefficients indicies
         poly_coeff_ind = 2:s;

It would be helpful to include the example input in the directory so a user can run the example without having to create the input file first.

For multi_stab_fun_2_rk_opt.m what is parsed by read_data didn't seem to match what is described in the comments and what is used in the multi_stab_fun_2_rk_opt function. I made the following changes

--- a/examples/multi_stab_fun_2_rk_opt.m
+++ b/examples/multi_stab_fun_2_rk_opt.m
@@ -29,7 +29,7 @@ function multi_rk = multi_stab_fun_2_rk_opt(input_file_name,class,objective,vara
 for i_stab_poly = 1:number_stab_poly
 
     % Read number of stage and order of accuracy of the RK scheme
-    [s,p,fp,d]= read_data(read_fid);
+    [s,p,fp,poly_coeff_val]= read_data(read_fid);
         
     % Set tall tree numbers (indices) and tall tree values that will be 
     % used to enforce the stability coefficients. Indeed, when "s" is 
@@ -52,14 +52,15 @@ for i_stab_poly = 1:number_stab_poly
     end
     
     % Polynomial coefficients indices
-    poly_coeff_ind = p+1:s;
-
-    % Polynomial coefficients values
-    poly_coeff_val = d(7+p+1:length(d));
+    if (fp > 0)
+        poly_coeff_ind = p+fp:s;
+    else
+        poly_coeff_ind = [];
+    end
 
     % Call to rk_opt
     rk = rk_opt(s,p,class,objective,'poly_coeff_ind',poly_coeff_ind,...
-                'poly_coeff_val',poly_coeff_val,varargin{:})
+                'poly_coeff_val',poly_coeff_val,varargin{:});
     
 end
 
@@ -105,8 +106,8 @@ end
 
 % =========================================================================
 
-function [s,p,fp,d]= read_data(read_fid)
-%function [s,p,fp]= read_data(read_fid)
+function [s,p,fp,poly_coeff]= read_data(read_fid)
+%function [s,p,fp,poly_coeff]= read_data(read_fid)
 %
 % Read number of stage and order of accuracy of the method
 
@@ -127,6 +128,9 @@ p = floor(d(2));
 % the order of accuracy of the scheme for nonlinear problems.
 fp = floor(d(3));
 
+% Polynomial coefficients values
+poly_coeff = d(4:length(d));
+
 end
 % =========================================================================

and ran the example with with the following input file:

stability poly
1

stages order free params. poly coeff

5      4     1            1 1 0.5  0.166666666666667 0.041666666666667 0.009615384615385

I'm sure it this is the intended usage but it seemed to work as expected. I think it would be helpful in include a sample input file for this example as well.

For both of the above examples it would also be helpful to include some comments about the expected output or maybe have an example output file for a user to compare against.

I also tried out the examples given in the opt_poly_bisect documentation as well as other setups from KA12. I've included the commands I ran below in case anyone else also wants to try reproducing plots from KA12. I'm not sure if it was intentional but I noticed the gap size is fixed at 20 when using the 'gap' option in spectrum.

%% Doc example
lam = spectrum('realaxis',500);
s = 10; p = 2;
[h,poly_coeff] = opt_poly_bisect(lam,s,p,'chebyshev')
plotstabreg_func(poly_coeff,[1])

%% KA12 - Fig 3 top
lam = spectrum('realaxis',6400);
s = 20; p = 4;
[h,poly_coeff] = opt_poly_bisect(lam,s,p,'chebyshev')
plotstabreg_func(poly_coeff,[1])

%% KA12 - Fig 3 bottom
lam = spectrum('realaxis',6400);
s = 20; p = 10;
[h,poly_coeff] = opt_poly_bisect(lam,s,p,'chebyshev')
plotstabreg_func(poly_coeff,[1],[-140 1 -10 10])

%% KA12 - Fig 4 left
lam = spectrum('imagaxis',3200);
s = 7; p = 1;
[h,poly_coeff] = opt_poly_bisect(lam,s,p,'rotated chebyshev')
plotstabreg_func(poly_coeff,[1])

%% KA12 - Fig 4 right
lam = spectrum('imagaxis',3200);
s = 7; p = 4;
[h,poly_coeff] = opt_poly_bisect(lam,s,p,'rotated chebyshev')
plotstabreg_func(poly_coeff,[1])

%% KA12 - Fig 6 left
lam = spectrum('disk',500);
s = 8; p = 3;
[h,poly_coeff] = opt_poly_bisect(lam,s,p,'monomial')
plotstabreg_func(poly_coeff,[1])

%% KA12 - Fig 6 right
lam = spectrum('disk',500);
s = 15; p = 4;
[h,poly_coeff] = opt_poly_bisect(lam,s,p,'monomial')
plotstabreg_func(poly_coeff,[1])

%% KA12 - Fig 8 (size of gap is hard coded to 20)
lam = spectrum('gap',500);
s = 6; p = 1;
[h,poly_coeff] = opt_poly_bisect(lam,s,p,'chebyshev')
plotstabreg_func(poly_coeff,[1],[-50 2 -10 10])

%% KA12 - Fig 10 top
lam_func = @(kappa) spectrum('rectangle',100,kappa,1)
[h,poly_coeff] = opt_poly_bisect(lam,10,1,'chebyshev','lam_func',lam_func)
plotstabreg_func(poly_coeff,[1])

%% KA12 - Fig 10 bottom
lam_func = @(kappa) spectrum('rectangle',100,kappa,10)
[h,poly_coeff] = opt_poly_bisect(lam,20,1,'chebyshev','lam_func',lam_func)

Add test suite in am_radius-opt.
Test:

-Rksp
-Rk_imp
-Rkp_dw
-Rk_imp_dw

1. Currently, we have

   Matlab scripts to search for Runge-Kutta methods that are optimal in terms of SSP coefficient

   as brief description of RK-Opt on GitHub but

   A package for the design of numerical ODE solvers

   in the draft.

2. Currently, the title on GitHub is

   RK-opt

   but we use

   RK-Opt

   in the draft. You can just rename the repository on GitHub - capitalization does not matter in URLs (and GitHub would forward all traffic otherwise).

Example s = 2
p = 1

The size of the optimization vector does not allow the assignment operation at line 103.

If this case and similar ones we should throw an exception and print an error message.

JOSS Review 2514

The Automated design of Runge-Kutta methods section in the documentation includes some nice motivational material about the use of additional stages to improve method properties and the intractability of finding optimal solutions for the larger set of coefficients by hand. I think it would be good to include some of this at the start of the "Statement of need" section of the JOSS software paper.

With Matlab 2020b and CVX Build 1148, when I run test.m, I get a few failed tests. Please see below.

                         Name                          Passed    Failed    Incomplete    Duration      Details   
    _______________________________________________    ______    ______    __________    ________    ____________

    {'test_rkopt/TestSSP32'                       }    false     true        true         0.14553    {1x1 struct}
    {'test_rkopt/TestSSP43'                       }    false     true        true         0.18959    {1x1 struct}
    {'test_rkopt/TestSDIRKSSP22'                  }    false     true        true         0.17585    {1x1 struct}
    {'test_rkopt/TestDIRKSSP22'                   }    false     true        true         0.23504    {1x1 struct}
    {'test_rkopt/TestRK22Acc'                     }    false     true        true         0.19243    {1x1 struct}
    {'test_rkopt/TestSSPTSRK2'                    }    true      false       false        0.54324    {1x1 struct}
    {'test_rkopt/TestSSP32Multistart'             }    false     true        true        0.077666    {1x1 struct}
    {'test_rkopt/TestSSP32MultistartParallel'     }    false     true        true          18.747    {1x1 struct}
    {'test_rkopt/TestSSP43Multistart'             }    false     true        true         0.14351    {1x1 struct}
    {'test_rkopt/TestSSP43MultistartParallel'     }    false     true        true          20.586    {1x1 struct}
    {'test_rkopt/TestSDIRKSSP22Multistart'        }    false     true        true         0.15462    {1x1 struct}
    {'test_rkopt/TestSDIRKSSP22MultistartParallel'}    false     true        true          18.075    {1x1 struct}
    {'test_rkopt/TestDIRKSSP22Multistart'         }    false     true        true         0.17769    {1x1 struct}
    {'test_rkopt/TestDIRKSSP22MultistartParallel' }    false     true        true          17.247    {1x1 struct}
    {'test_rkopt/TestRK22AccMultistart'           }    false     true        true         0.22198    {1x1 struct}
    {'test_rkopt/TestRK22AccMultistartParallel'   }    false     true        true          19.283    {1x1 struct}
    {'test_rkopt/TestRK323SAcc'                   }    false     true        true          3.9953    {1x1 struct}
    {'test_rkopt/TestEmbeddedRK4353Ss'            }    true      false       false        0.00768    {1x1 struct}

rk_opt is designed to run with a number of initial guesses and a number of processes, and the intent is that each parallel process will use an independent (random) set of initial guesses. By default, MATLAB generates the same sequence of random numbers in each new session, so we use rng('shuffle') in order to do searches that are truly independent across different sessions, different machines, etc. However, I just found out (see https://www.mathworks.com/help/parallel-computing/control-random-number-streams-on-workers.html):

Because rng('shuffle') seeds the random number generator based on the current time, do not use this command to set the random number stream on different workers if you want to ensure independent streams. This is especially true when the command is sent to multiple workers simultaneously, such as inside a parfor, spmd, or a communicating job. For independent streams on the workers, use the default behavior; or if that is not sufficient for your needs, consider using a unique substream on each worker using RandStream.

I've read the help on RandStream but haven't understood how to use it for what we want. @Sondar74 @ranocha @abhibsws if you can figure this out, let me know or feel free to submit a PR. Until then, it seems that setting np>1 actually does nothing.