File: numgradient.cc

package info (click to toggle)
octave-optim 1.5.3-2
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 2,368 kB
  • sloc: cpp: 1,437; makefile: 185; perl: 169; xml: 29; sh: 3
file content (163 lines) | stat: -rw-r--r-- 4,828 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
// Copyright (C) 2004, 2006 Michael Creel <michael.creel@uab.es>
//
// This program 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.
//
// This program 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
// this program; if not, see <http://www.gnu.org/licenses/>.

// numgradient: numeric central difference gradient

#include <oct.h>
#include <octave/parse.h>
#include <octave/lo-mappers.h>
#include <octave/Cell.h>
#include <float.h>

#include "error-helpers.h"

// argument checks
static bool
any_bad_argument(const octave_value_list& args)
{
  if (!args(0).is_string())
    {
      _p_error("numgradient: first argument must be string holding objective function name");
      return true;
    }

  if (!args(1).OV_ISCELL ())
    {
      _p_error("numgradient: second argument must cell array of function arguments");
      return true;
    }

  // minarg, if provided
  if (args.length() == 3)
    {
      int tmp;
      bool err;
      SET_ERR (tmp = args(2).int_value(), err);
      if (err)
        {
          _p_error ("numgradient: 3rd argument, if supplied,  must an integer\n\
that specifies the argument wrt which differentiation is done");
          return true;
        }
          if ((tmp > args(1).length ()) || (tmp < 1))
            {
              _p_error("numgradient: 3rd argument must be a positive integer that indicates \n\
which of the elements of the second argument is the\n\
one to differentiate with respect to");
              return true;
            }
    }
  return false;
}


DEFUN_DLD(numgradient, args, , "numgradient(f, {args}, minarg)\n\
\n\
Numeric central difference gradient of f with respect\n\
to argument \"minarg\".\n\
* first argument: function name (string)\n\
* second argument: all arguments of the function (cell array)\n\
* third argument: (optional) the argument to differentiate w.r.t.\n\
        (scalar, default=1)\n\
\n\
\"f\" may be vector-valued. If \"f\" returns\n\
an n-vector, and the argument is a k-vector, the gradient\n\
will be an nxk matrix\n\
\n\
Example:\n\
function a = f(x);\n\
        a = [x'*x; 2*x];\n\
endfunction\n\
numgradient(\"f\", {ones(2,1)})\n\
ans =\n\
\n\
  2.00000  2.00000\n\
  2.00000  0.00000\n\
  0.00000  2.00000\n\
")
{
  int nargin = args.length();
  if (!((nargin == 2)|| (nargin == 3)))
    {
      error("numgradient: you must supply 2 or 3 arguments");
      return octave_value_list();
    }

  // check the arguments
  if (any_bad_argument( args))
    {
      error ("error in numgradient");
      return octave_value_list();
    }

  std::string f (args(0).string_value());
  Cell f_args_cell (args(1).cell_value());
  octave_value_list f_args, f_return;
  Matrix obj_value, obj_left, obj_right;
  double SQRT_EPS, p, delta, diff;
  int i, j, k, n, minarg, test;

  // Default values for controls
  minarg = 1; // by default, first arg is one over which we minimize

  // copy cell contents over to octave_value_list to use feval()
  k = f_args_cell.numel ();
  f_args(k - 1); // resize only once
  for (i = 0; i<k; i++) f_args(i) = f_args_cell(i);

  // check which arg w.r.t which we need to differentiate
  if (args.length() == 3) minarg = args(2).int_value();
  Matrix parameter = f_args(minarg - 1).matrix_value();

  // initial function value
  f_return = OCTAVE__FEVAL (f, f_args);
  obj_value = f_return(0).matrix_value();

  n = obj_value.rows(); // find out dimension
  k = parameter.rows();
  Matrix derivative(n, k);
  Matrix columnj;

  for (j=0; j<k; j++)
    {
      // get 1st derivative by central difference
      p = parameter(j);

      // determine delta for finite differencing
      SQRT_EPS = sqrt(DBL_EPSILON);
      diff = exp(log(DBL_EPSILON)/3);
      test = (fabs(p) + SQRT_EPS) * SQRT_EPS > diff;
      if (test) delta = (fabs(p) + SQRT_EPS) * SQRT_EPS;
      else delta = diff;

      // right side
      parameter(j) = p + delta;
      f_args(minarg - 1) = parameter;
      f_return = OCTAVE__FEVAL (f, f_args);
      obj_right = f_return(0).matrix_value();

      // left size
      parameter(j) = p - delta;
      f_args(minarg - 1) = parameter;
      f_return = OCTAVE__FEVAL (f, f_args);
      obj_left = f_return(0).matrix_value();

      parameter(j) = p;  // restore original parameter
      columnj = (obj_right - obj_left) / (2*delta);
      for (i=0; i<n; i++) derivative(i, j) = columnj(i);
    }

  return octave_value(derivative);
}