changeset 22625:081a201b77c7 stable

Clean up ode options implementation to follow Octave coding standards. * known_option_names.m: Delete file * scripts/ode/module.mk: Remove known_option_names from build system. * ode23.m, ode45.m: Fix typo in docstring. Correct indentation. Remove trailing whitespace. * AbsRel_Norm.m: Use default for input argument to simplify function. Remove input validation for private, internal function. * odedefaults.m: Add docstring. Use persistent variables for performance. * odemergeopts.m: Fix indentation. * starting_stepsize.m: Show input func as '@func' in docstring.
author Carlo de Falco <carlo.defalco@polimi.it>
date Sat, 15 Oct 2016 10:30:48 +0200
parents 37b7b86f62f2
children 869c02fde46c
files scripts/ode/module.mk scripts/ode/ode23.m scripts/ode/ode45.m scripts/ode/private/AbsRel_Norm.m scripts/ode/private/known_option_names.m scripts/ode/private/odedefaults.m scripts/ode/private/odemergeopts.m scripts/ode/private/starting_stepsize.m
diffstat 8 files changed, 108 insertions(+), 150 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ode/module.mk	Sat Oct 15 02:51:30 2016 -0500
+++ b/scripts/ode/module.mk	Sat Oct 15 10:30:48 2016 +0200
@@ -6,7 +6,6 @@
   scripts/ode/private/AbsRel_Norm.m \
   scripts/ode/private/integrate_adaptive.m \
   scripts/ode/private/kahan.m \
-  scripts/ode/private/known_option_names.m \
   scripts/ode/private/odedefaults.m \
   scripts/ode/private/odemergeopts.m \
   scripts/ode/private/ode_event_handler.m \
--- a/scripts/ode/ode23.m	Sat Oct 15 02:51:30 2016 -0500
+++ b/scripts/ode/ode23.m	Sat Oct 15 10:30:48 2016 +0200
@@ -1,3 +1,4 @@
+## Copyright (C) 2016, Carlo de Falco
 ## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it>
 ## Copyright (C) 2014-2016 Jacopo Corno <jacopo.corno@gmail.com>
 ## Copyright (C) 2013-2016 Roberto Porcu' <roberto.porcu@polimi.it>
@@ -43,8 +44,8 @@
 ##
 ## By default, @code{ode23} uses an adaptive timestep with the
 ## @code{integrate_adaptive} algorithm.  The tolerance for the timestep
-## computation may be changed by using the options @qcode{"RelTol"},
-## and @qcode{"AbsTol"},. 
+## computation may be changed by using the options @qcode{"RelTol"}
+## and @qcode{"AbsTol"}.
 ##
 ## @var{init} contains the initial value for the unknowns.  If it is a row
 ## vector then the solution @var{y} will be a matrix in which each column is
@@ -211,7 +212,7 @@
                                              init, odeopts.AbsTol,
                                              odeopts.RelTol,
                                              strcmp (odeopts.NormControl,
-                                             "on"), odeopts.funarguments);
+                                                     "on"), odeopts.funarguments);
   endif
 
 
@@ -231,18 +232,18 @@
     if (! strcmp (odeopts.MStateDependence, "none")) # constant mass matrices have already
       mass = @(t,x) odeopts.Mass (t, x, odeopts.funarguments{:});
       fun = @(t,x) mass (t, x, odeopts.funarguments{:}) ...
-        \ fun (t, x, odeopts.funarguments{:});
+            \ fun (t, x, odeopts.funarguments{:});
     else                 # if ((! strcmp (odeopts.MStateDependence, "none")) == false)
       mass = @(t) odeopts.Mass (t, odeopts.funarguments{:});
       fun = @(t,x) mass (t, odeopts.funarguments{:}) ...
-        \ fun (t, x, odeopts.funarguments{:});
+            \ fun (t, x, odeopts.funarguments{:});
     endif
   endif
 
-  
+
   solution = integrate_adaptive (@runge_kutta_23, ...
                                  order, fun, trange, init, odeopts);
-    
+
 
   ## Postprocessing, do whatever when terminating integration algorithm
   if (odeopts.haveoutputfunction)  # Cleanup plotter
--- a/scripts/ode/ode45.m	Sat Oct 15 02:51:30 2016 -0500
+++ b/scripts/ode/ode45.m	Sat Oct 15 10:30:48 2016 +0200
@@ -1,3 +1,4 @@
+## Copyright (C) 2016, Carlo de Falco
 ## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it>
 ## Copyright (C) 2014-2016 Jacopo Corno <jacopo.corno@gmail.com>
 ## Copyright (C) 2013-2016 Roberto Porcu' <roberto.porcu@polimi.it>
@@ -41,7 +42,7 @@
 ##
 ## By default, @code{ode45} uses an adaptive timestep with the
 ## @code{integrate_adaptive} algorithm.  The tolerance for the timestep
-## computation may be changed by using the options @qcode{"RelTol"},
+## computation may be changed by using the options @qcode{"RelTol"}
 ## and @qcode{"AbsTol"}.
 ##
 ## @var{init} contains the initial value for the unknowns.  If it is a row
@@ -162,7 +163,7 @@
                                      "InitialSlope", "MaxOrder", "BDF"});
   attributes = rmfield (attributes, {"Jacobian", "JPattern", "Vectorized", ...
                                      "MvPattern", "MassSingular", ...
-                                     "InitialSlope", "MaxOrder", "BDF"}); 
+                                     "InitialSlope", "MaxOrder", "BDF"});
 
   odeopts = odemergeopts (odeopts, defaults, classes, attributes, 'ode45');
 
@@ -195,8 +196,8 @@
                                              init, odeopts.AbsTol,
                                              odeopts.RelTol,
                                              strcmp (odeopts.NormControl,
-                                             "on"), odeopts.funarguments);
-  endif 
+                                                     "on"), odeopts.funarguments);
+  endif
 
 
   if (! isempty (odeopts.Mass) && isnumeric (odeopts.Mass))
@@ -215,18 +216,18 @@
     if (! strcmp (odeopts.MStateDependence, "none")) # constant mass matrices have already
       mass = @(t,x) odeopts.Mass (t, x, odeopts.funarguments{:});
       fun = @(t,x) mass (t, x, odeopts.funarguments{:}) ...
-             \ fun (t, x, odeopts.funarguments{:});
+            \ fun (t, x, odeopts.funarguments{:});
     else                 # if ((! strcmp (odeopts.MStateDependence, "none")) == false)
       mass = @(t) odeopts.Mass (t, odeopts.funarguments{:});
       fun = @(t,x) mass (t, odeopts.funarguments{:}) ...
-             \ fun (t, x, odeopts.funarguments{:});
+            \ fun (t, x, odeopts.funarguments{:});
     endif
   endif
 
- 
+
   solution = integrate_adaptive (@runge_kutta_45_dorpri,
                                  order, fun, trange, init, odeopts);
-  
+
 
   ## Postprocessing, do whatever when terminating integration algorithm
   if (odeopts.haveoutputfunction)  # Cleanup plotter
--- a/scripts/ode/private/AbsRel_Norm.m	Sat Oct 15 02:51:30 2016 -0500
+++ b/scripts/ode/private/AbsRel_Norm.m	Sat Oct 15 10:30:48 2016 +0200
@@ -22,23 +22,10 @@
 ## Undocumented internal function.
 ## @end deftypefn
 
-function retval = AbsRel_Norm (x, x_old, AbsTol, RelTol, normcontrol, y)
-
-  n = length (x);
-
-  if (nargin == 5)
-    y = zeros (size (x));
-  endif
+function retval = AbsRel_Norm (x, x_old, AbsTol, RelTol, normcontrol, y = zeros (size (x)))
 
-  if (length (x_old) != n || length (y) != n)
-    error ("Octave:invalid-input-arg", "invalid dimensions of input arguments");
-  endif
-
-  if ((length (AbsTol) != 1 && length (AbsTol) != n)
-      || (length (RelTol) != 1 && length (RelTol) != n))
-    error ("Octave:invalid-input-arg", "invalid dimensions of input arguments");
-  endif
-
+  n = numel (x);
+  
   sc = AbsTol + max (abs (x), abs (x_old)) .* RelTol;
   if (normcontrol)
     retval = max (abs (x - y) ./ sc);
--- a/scripts/ode/private/known_option_names.m	Sat Oct 15 02:51:30 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,35 +0,0 @@
-## Copyright (C) 2016 Carlo de Falco
-##
-## This file is part of Octave.
-##
-## Octave is free software; you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 3 of the License, or (at
-## your option) any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-## General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {} {@var{val} =} known_option_names ()
-## Return a list of known names for ode options.
-## @seealso{odeset, odeget}
-## @end deftypefn
-
-function ret = known_option_names ()
-
-ret = {"AbsTol"; "BDF"; "Events"; "InitialSlope";
-       "InitialStep"; "Jacobian"; "JConstant"; "JPattern";
-       "Mass"; "MassConstant"; "MassSingular"; "MaxOrder";
-       "MaxStep"; "MStateDependence"; "MvPattern";
-       "NonNegative"; "NormControl"; "OutputFcn"; "OutputSel";
-       "Refine"; "RelTol"; "Stats"; "Vectorized";
-       "TimeStepSize"; "TimeStepNumber"};
-
-endfunction
--- a/scripts/ode/private/odedefaults.m	Sat Oct 15 02:51:30 2016 -0500
+++ b/scripts/ode/private/odedefaults.m	Sat Oct 15 10:30:48 2016 +0200
@@ -1,3 +1,4 @@
+## Copyright (C) 2016, Carlo de Falco
 ## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it>
 ##
 ## This file is part of Octave.
@@ -16,81 +17,85 @@
 ## along with Octave; see the file COPYING.  If not, see
 ## <http://www.gnu.org/licenses/>.
 
+## -*- texinfo -*-
+## @deftypefn {} {[@var{d}, @var{c}, @var{a}] =} odedefaults (@var{n}, @var{t0}, @var{tf})
+## Undocumented internal function.
+## @end deftypefn
+
 function [defaults, classes, attributes] = odedefaults (n, t0, tf)
 
-defaults = odeset ('AbsTol', 1e-6,
-                   'BDF', 'off',
-                   'Events', [],
-                   'InitialSlope', zeros(n,1),
-                   'InitialStep', [],
-                   'Jacobian', [],
-                   'JConstant', 'off',
-                   'JPattern', [],
-                   'Mass', [],
-                   'MassConstant', 'off',
-                   'MassSingular', 'maybe',
-                   'MaxOrder', 5,
-                   'MaxStep', 0.1*abs(t0-tf),
-                   'MStateDependence', 'weak',
-                   'MvPattern', [],
-                   'NonNegative', [],
-                   'NormControl', 'off',
-                   'OutputFcn', [],
-                   'OutputSel', [],
-                   'Refine', 1,  
-                   'RelTol', 1e-3,
-                   'Stats', 'off',
-                   'Vectorized', 'off');
- 
-classes = odeset ('Abstol', {"float"},
-                  'BDF', "char",
-                  'Events', {"function_handle"},
-                  'InitialSlope', {"float"},
-                  'InitialStep', {"float"},
-                  'Jacobian', {"float", "function_handle", "cell"},
-                  'JConstant', "char",
-                  'JPattern', {"float"},
-                  'Mass', {"float", "function_handle"},
-                  'MassConstant', "char",
-                  'MassSingular', "char",
-                  'MaxOrder', {"float"},
-                  'MaxStep', {"float"},
-                  'MStateDependence', "char",
-                  'MvPattern', {"float"},
-                  'NonNegative', {"float"},
-                  'NormControl', "char",
-                  'OutputFcn', {"function_handle"},
-                  'OutputSel', {"float"},
-                  'Refine', {"float"},
-                  'RelTol', {"float"},
-                  'Stats', "char",
-                  'Vectorized', "char");
-
-
-##FIXME: How can I check Jacobian where it's a cell????? Maybe it's better to check it inside the solver
-##FIXME: Vectorized can be a cell of stings
-attributes = odeset ('AbsTol', {"real", "vector", "positive"},
-                     'BDF', {"on", "off"},
-                     'Events', {},
-                     'InitialSlope', {"real", "vector", "numel", n},
-                     'InitialStep', {"positive", "scalar"},
-                     'Jacobian', {},
-                     'JConstant', {"on", "off"},
-                     'JPattern', {"vector"},
-                     'Mass', {},
-                     'MassConstant', {"on", "off"},
-                     'MassSingular', {"no", "maybe", "yes"},
-                     'MaxOrder', {">=", 0, "<=", 5, "integer"},
-                     'MaxStep', {"positive", "scalar", "real"},
-                     'MStateDependence', {"weak", "strong", "none"},
-                     'MvPattern', {"vector"},
-                     'NonNegative', {"vector", "integer", "positive"},
-                     'NormControl', {"on", "off"},
-                     'OutputFcn', {},
-                     'OutputSel', {"vector", "integer", "positive",...
-                                   ">", 0, "<=", n},
-                     'Refine', {"scalar", ">", 0, "integer"},
-                     'RelTol', {"scalar", "positive", "real"},
-                     'Stats', {"on", "off"},
-                     'Vectorized', {"on", "off"});
+  persistent defaults = struct ("AbsTol", 1e-6,
+                                "BDF", "off",
+                                "Events", [],
+                                "InitialSlope", zeros(n,1),
+                                "InitialStep", [],
+                                "Jacobian", [],
+                                "JConstant", "off",
+                                "JPattern", [],
+                                "Mass", [],
+                                "MassConstant", "off",
+                                "MassSingular", "maybe",
+                                "MaxOrder", 5,
+                                "MaxStep", 0.1 * abs (t0-tf),
+                                "MStateDependence", "weak",
+                                "MvPattern", [],
+                                "NonNegative", [],
+                                "NormControl", "off",
+                                "OutputFcn", [],
+                                "OutputSel", [],
+                                "Refine", 1,  
+                                "RelTol", 1e-3,
+                                "Stats", "off",
+                                "Vectorized", "off");
+  
+  defaults.MaxStep = (0.1 * abs (t0-tf));
+  
+  persistent classes = struct ("AbsTol", {{"float"}},
+                               "BDF", "char",
+                               "Events", {{"function_handle"}},
+                               "InitialSlope", {{"float"}},
+                               "InitialStep", {{"float"}},
+                               "Jacobian", {{"float", "function_handle", "cell"}},
+                               "JConstant", "char",
+                               "JPattern", {{"float"}},
+                               "Mass", {{"float", "function_handle"}},
+                               "MassConstant", "char",
+                               "MassSingular", "char",
+                               "MaxOrder", {{"float"}},
+                               "MaxStep", {{"float"}},
+                               "MStateDependence", "char",
+                               "MvPattern", {{"float"}},
+                               "NonNegative", {{"float"}},
+                               "NormControl", "char",
+                               "OutputFcn", {{"function_handle"}},
+                               "OutputSel", {{"float"}},
+                               "Refine", {{"float"}},
+                               "RelTol", {{"float"}},
+                               "Stats", "char",
+                               "Vectorized", "char");
+  
+  persistent attributes = struct ("AbsTol", {{"real", "vector", "positive"}},
+                                  "BDF", {{"on", "off"}},
+                                  "Events", {{}},
+                                  "InitialSlope", {{"real", "vector", "numel", n}},
+                                  "InitialStep", {{"positive", "scalar"}},
+                                  "Jacobian", {{}},
+                                  "JConstant", {{"on", "off"}},
+                                  "JPattern", {{"vector"}},
+                                  "Mass", {{}},
+                                  "MassConstant", {{"on", "off"}},
+                                  "MassSingular", {{"no", "maybe", "yes"}},
+                                  "MaxOrder", {{">=", 0, "<=", 5, "integer"}},
+                                  "MaxStep", {{"positive", "scalar", "real"}},
+                                  "MStateDependence", {{"weak", "strong", "none"}},
+                                  "MvPattern", {{"vector"}},
+                                  "NonNegative", {{"vector", "integer", "positive"}},
+                                  "NormControl", {{"on", "off"}},
+                                  "OutputFcn", {{}},
+                                  "OutputSel", {{"vector", "integer", "positive",...
+                                                 ">", 0, "<=", n}},
+                                  "Refine", {{"scalar", ">", 0, "integer"}},
+                                  "RelTol", {{"scalar", "positive", "real"}},
+                                  "Stats", {{"on", "off"}},
+                                  "Vectorized", {{"on", "off"}});
 endfunction
--- a/scripts/ode/private/odemergeopts.m	Sat Oct 15 02:51:30 2016 -0500
+++ b/scripts/ode/private/odemergeopts.m	Sat Oct 15 10:30:48 2016 +0200
@@ -32,11 +32,11 @@
 
       else
         error ("Octave:invalid-input-arg",
-                [fun_name ": invalid value assigned to field '%s'"], key);
+               [fun_name ": invalid value assigned to field '%s'"], key);
       endif
-      
-    options.(key) = useroptions.(key);
-    
+
+      options.(key) = useroptions.(key);
+
     endif
   endfor
 endfunction
--- a/scripts/ode/private/starting_stepsize.m	Sat Oct 15 02:51:30 2016 -0500
+++ b/scripts/ode/private/starting_stepsize.m	Sat Oct 15 10:30:48 2016 +0200
@@ -17,12 +17,12 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn {} {@var{h} =} starting_stepsize (@var{order}, @var{@@func}, @var{t0}, @var{x0}, @var{AbsTol}, @var{RelTol}, @var{normcontrol})
+## @deftypefn {} {@var{h} =} starting_stepsize (@var{order}, @var{func}, @var{t0}, @var{x0}, @var{AbsTol}, @var{RelTol}, @var{normcontrol})
 ##
 ## Determine a good initial timestep for an ODE solver of order @var{order}
 ## using the algorithm described in reference [1].
 ##
-## The input argument @var{@@func}, is the function describing the differential
+## The input argument @var{func}, is the function describing the differential
 ## equations, @var{t0} is the initial time, and @var{x0} is the initial
 ## condition.  @var{AbsTol} and @var{RelTol} are the absolute and relative
 ## tolerance on the ODE integration taken from an ode options structure.
@@ -67,12 +67,12 @@
        AbsRel_Norm (yh - y, yh - y, AbsTol, RelTol, normcontrol);
 
   if (max (d1, d2) <= 1e-15)
-    h1 = max (1e-6, h0*1e-3);
+    h1 = max (1e-6, h0 * 1e-3);
   else
     h1 = (1e-2 / max (d1, d2)) ^(1 / (order+1));
   endif
 
-  h = min (100*h0, h1);
+  h = min (100 * h0, h1);
 
 endfunction