octave_packages/optim-1.2.0/cdiff.m

   1 ## Copyright (C) 2002 Etienne Grossmann <etienne@egdn.net>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## c = cdiff (func,wrt,N,dfunc,stack,dx) - Code for num. differentiation
  17 ##   = "function df = dfunc (var1,..,dvar,..,varN) .. endfunction
  18 ##
  19 ## Returns a string of octave code that defines a function 'dfunc' that
  20 ## returns the derivative of 'func' with respect to it's 'wrt'th
  21 ## argument.
  22 ##
  23 ## The derivatives are obtained by symmetric finite difference.
  24 ##
  25 ## dfunc()'s return value is in the same format as that of  ndiff()
  26 ##
  27 ## func  : string : name of the function to differentiate
  28 ##
  29 ## wrt   : int    : position, in argument list, of the differentiation
  30 ##                  variable.                                Default:1
  31 ##
  32 ## N     : int    : total number of arguments taken by 'func'.
  33 ##                  If N=inf, dfunc will take variable argument list.
  34 ##                                                         Default:wrt
  35 ##
  36 ## dfunc : string : Name of the octave function that returns the
  37 ##                   derivatives.                   Default:['d',func]
  38 ##
  39 ## stack : string : Indicates whether 'func' accepts vertically
  40 ##                  (stack="rstack") or horizontally (stack="cstack")
  41 ##                  arguments. Any other string indicates that 'func'
  42 ##                  does not allow stacking.                Default:''
  43 ##
  44 ## dx    : real   : Step used in the symmetric difference scheme.
  45 ##                                                  Default:10*sqrt(eps)
  46 ##
  47 ## See also : ndiff, eval, todisk
  48
  49 function c = cdiff (func,wrt,nargs,dfunc,stack,dx)
  50
  51 if nargin<2,
  52   wrt = 1 ;
  53 end
  54 if nargin<3,
  55   nargs = wrt ;
  56 end
  57 if nargin<4 || strcmp(dfunc,""),
  58   dfunc = ["d",func] ;
  59   if exist(dfunc)>=2,
  60     printf(["cdiff : Warning : name of derivative not specified\n",\
  61             "        and canonic name '%s' is already taken\n"],\
  62            dfunc);
  63     ## keyboard
  64   end
  65 end
  66 if nargin<5, stack = "" ; end
  67 if nargin<6, dx = 10*sqrt(eps)  ; end
  68
  69 ## verbose = 0 ;
  70 ## build argstr = "var1,..,dvar,...var_nargs"
  71 if isfinite (nargs)
  72   argstr = sprintf("var%i,",1:nargs);
  73 else
  74   argstr = [sprintf("var%i,",1:wrt),"...,"];
  75 end
  76
  77 argstr = strrep(argstr,sprintf("var%i",wrt),"dvar") ;
  78 argstr = argstr(1:length(argstr)-1) ;
  79
  80 if strcmp("cstack",stack) ,     # Horizontal stacking ################
  81
  82   calstr = "reshape (kron(ones(1,2*ps), dvar(:))+[-dx*eye(ps),dx*eye(ps)], sz.*[1,2*ps])";
  83   calstr = strrep(argstr,"dvar",calstr) ;
  84   calstr = sprintf("%s(%s)",func,calstr) ;
  85
  86   calstr = sprintf(strcat("  res = %s;\n",
  87                           "  pr = prod (size(res)) / (2*ps);\n",
  88                           "  res = reshape (res,pr,2*ps);\n",
  89                           "  df = (res(:,ps+1:2*ps)-res(:,1:ps)) / (2*dx);\n"),
  90                    calstr) ;
  91
  92
  93 elseif strcmp("rstack",stack),  # Vertical stacking ##################
  94
  95   calstr = "kron(ones(2*ps,1),dvar)+dx*[-dv;dv]" ;
  96   calstr = strrep(argstr,"dvar",calstr) ;
  97   calstr = sprintf("%s(%s)",func,calstr) ;
  98
  99   calstr = sprintf(strcat("  dv = kron (eye(sz(2)), eye(sz(1))(:));\n",\
 100                           "  res = %s;\n",\
 101                           "  sr = size(res)./[2*ps,1];\n",\
 102                           "  pr = prod (sr);\n",\
 103                           "  df = (res(sr(1)*ps+1:2*sr(1)*ps,:)-res(1:sr(1)*ps,:))/(2*dx);\n",\
 104                           "  scramble = reshape (1:pr,sr(2),sr(1))';\n",\
 105                           "  df = reshape (df',pr,ps)(scramble(:),:);\n"),\
 106                    calstr) ;
 107   ## sayif(verbose,"cdiff : calstr='%s'\n",calstr) ;
 108 else                            # No stacking ########################
 109   calstr = sprintf("%s (%s)",func,argstr) ;
 110   ## "func(var1,dvar%sdv(:,%d:%d),...,varN),"
 111   ## calstr = strrep(calstr,"dvar","dvar%sdv(:,(i-1)*sz(2)+1:i*sz(2))")(:)';
 112
 113   calstr = strrep(calstr,"dvar","dvar%sdv")(:)';
 114
 115   ## func(..,dvar+dv(:,1:sz(2)),..) - func(..)
 116   calstr = strcat(calstr,"-",calstr) ; ## strcat(calstr,"-",calstr) ;
 117   calstr = sprintf(calstr,"+","-") ;
 118   tmp = calstr ;
 119   ## sayif(verbose,"cdiff : calstr='%s'\n",calstr) ;
 120   calstr = sprintf(strcat("  dv = zeros (sz); dv(1) = dx;\n",\
 121                           "  df0 = %s;\n",\
 122                           "  sr = size (df0);\n",\
 123                           "  df = zeros(prod (sr),ps); df(:,1) = df0(:);\n",\
 124                           "  for i = 2:ps,\n",\
 125                           "     dv(i) = dx; dv(i-1) = 0;\n",\
 126                           "     df(:,i) = (%s)(:);\n",\
 127                           "  end;\n",\
 128                           "  df ./= 2*dx;\n"
 129                           ),
 130                    calstr, tmp) ;
 131
 132
 133   ## sayif(verbose,"cdiff : calstr='%s'\n",calstr) ;
 134
 135   ## "func(var1,reshape(dvar(1:NV,1),SZ1,SZ2),...,varN)," ,
 136   ## "func(var1,reshape(dvar(1:NV,2),SZ1,SZ2),...,varN)," , ...
 137   ## "func(var1,reshape(dvar(1:NV,NP),SZ1,SZ2),...,varN)"
 138   ## sayif(verbose,"cdiff : calstr='%s'\n",calstr) ;
 139 end
 140 argstr = strrep (argstr, "...", "varargin");
 141 calstr = strrep (calstr, "...", "varargin{:}");
 142
 143 c = sprintf(strcat("function df = %s (%s)\n",\
 144                    "  ## Numerical differentiation of '%s' wrt to it's %d'th argument\n",\
 145                    "  ## This function has been written by 'cdiff()'\n",\
 146                    "  dx = %e;\n",\
 147                    "  sz = size (dvar);\n",\
 148                    "  ps = prod (sz);\n",\
 149                    "%s",\
 150                    "endfunction\n"),\
 151             dfunc,argstr,\
 152             func,wrt,\
 153             dx,\
 154             calstr) ;
 155