X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Foptim-1.2.0%2Fcdiff.m;fp=octave_packages%2Foptim-1.2.0%2Fcdiff.m;h=7d1cab17a00711b2d8bec4502b385d103430f42e;hp=0000000000000000000000000000000000000000;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hpb=1705066eceaaea976f010f669ce8e972f3734b05 diff --git a/octave_packages/optim-1.2.0/cdiff.m b/octave_packages/optim-1.2.0/cdiff.m new file mode 100644 index 0000000..7d1cab1 --- /dev/null +++ b/octave_packages/optim-1.2.0/cdiff.m @@ -0,0 +1,155 @@ +## Copyright (C) 2002 Etienne Grossmann +## +## 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 . + +## c = cdiff (func,wrt,N,dfunc,stack,dx) - Code for num. differentiation +## = "function df = dfunc (var1,..,dvar,..,varN) .. endfunction +## +## Returns a string of octave code that defines a function 'dfunc' that +## returns the derivative of 'func' with respect to it's 'wrt'th +## argument. +## +## The derivatives are obtained by symmetric finite difference. +## +## dfunc()'s return value is in the same format as that of ndiff() +## +## func : string : name of the function to differentiate +## +## wrt : int : position, in argument list, of the differentiation +## variable. Default:1 +## +## N : int : total number of arguments taken by 'func'. +## If N=inf, dfunc will take variable argument list. +## Default:wrt +## +## dfunc : string : Name of the octave function that returns the +## derivatives. Default:['d',func] +## +## stack : string : Indicates whether 'func' accepts vertically +## (stack="rstack") or horizontally (stack="cstack") +## arguments. Any other string indicates that 'func' +## does not allow stacking. Default:'' +## +## dx : real : Step used in the symmetric difference scheme. +## Default:10*sqrt(eps) +## +## See also : ndiff, eval, todisk + +function c = cdiff (func,wrt,nargs,dfunc,stack,dx) + +if nargin<2, + wrt = 1 ; +end +if nargin<3, + nargs = wrt ; +end +if nargin<4 || strcmp(dfunc,""), + dfunc = ["d",func] ; + if exist(dfunc)>=2, + printf(["cdiff : Warning : name of derivative not specified\n",\ + " and canonic name '%s' is already taken\n"],\ + dfunc); + ## keyboard + end +end +if nargin<5, stack = "" ; end +if nargin<6, dx = 10*sqrt(eps) ; end + +## verbose = 0 ; +## build argstr = "var1,..,dvar,...var_nargs" +if isfinite (nargs) + argstr = sprintf("var%i,",1:nargs); +else + argstr = [sprintf("var%i,",1:wrt),"...,"]; +end + +argstr = strrep(argstr,sprintf("var%i",wrt),"dvar") ; +argstr = argstr(1:length(argstr)-1) ; + +if strcmp("cstack",stack) , # Horizontal stacking ################ + + calstr = "reshape (kron(ones(1,2*ps), dvar(:))+[-dx*eye(ps),dx*eye(ps)], sz.*[1,2*ps])"; + calstr = strrep(argstr,"dvar",calstr) ; + calstr = sprintf("%s(%s)",func,calstr) ; + + calstr = sprintf(strcat(" res = %s;\n", + " pr = prod (size(res)) / (2*ps);\n", + " res = reshape (res,pr,2*ps);\n", + " df = (res(:,ps+1:2*ps)-res(:,1:ps)) / (2*dx);\n"), + calstr) ; + + +elseif strcmp("rstack",stack), # Vertical stacking ################## + + calstr = "kron(ones(2*ps,1),dvar)+dx*[-dv;dv]" ; + calstr = strrep(argstr,"dvar",calstr) ; + calstr = sprintf("%s(%s)",func,calstr) ; + + calstr = sprintf(strcat(" dv = kron (eye(sz(2)), eye(sz(1))(:));\n",\ + " res = %s;\n",\ + " sr = size(res)./[2*ps,1];\n",\ + " pr = prod (sr);\n",\ + " df = (res(sr(1)*ps+1:2*sr(1)*ps,:)-res(1:sr(1)*ps,:))/(2*dx);\n",\ + " scramble = reshape (1:pr,sr(2),sr(1))';\n",\ + " df = reshape (df',pr,ps)(scramble(:),:);\n"),\ + calstr) ; + ## sayif(verbose,"cdiff : calstr='%s'\n",calstr) ; +else # No stacking ######################## + calstr = sprintf("%s (%s)",func,argstr) ; + ## "func(var1,dvar%sdv(:,%d:%d),...,varN)," + ## calstr = strrep(calstr,"dvar","dvar%sdv(:,(i-1)*sz(2)+1:i*sz(2))")(:)'; + + calstr = strrep(calstr,"dvar","dvar%sdv")(:)'; + + ## func(..,dvar+dv(:,1:sz(2)),..) - func(..) + calstr = strcat(calstr,"-",calstr) ; ## strcat(calstr,"-",calstr) ; + calstr = sprintf(calstr,"+","-") ; + tmp = calstr ; + ## sayif(verbose,"cdiff : calstr='%s'\n",calstr) ; + calstr = sprintf(strcat(" dv = zeros (sz); dv(1) = dx;\n",\ + " df0 = %s;\n",\ + " sr = size (df0);\n",\ + " df = zeros(prod (sr),ps); df(:,1) = df0(:);\n",\ + " for i = 2:ps,\n",\ + " dv(i) = dx; dv(i-1) = 0;\n",\ + " df(:,i) = (%s)(:);\n",\ + " end;\n",\ + " df ./= 2*dx;\n" + ), + calstr, tmp) ; + + + ## sayif(verbose,"cdiff : calstr='%s'\n",calstr) ; + + ## "func(var1,reshape(dvar(1:NV,1),SZ1,SZ2),...,varN)," , + ## "func(var1,reshape(dvar(1:NV,2),SZ1,SZ2),...,varN)," , ... + ## "func(var1,reshape(dvar(1:NV,NP),SZ1,SZ2),...,varN)" + ## sayif(verbose,"cdiff : calstr='%s'\n",calstr) ; +end +argstr = strrep (argstr, "...", "varargin"); +calstr = strrep (calstr, "...", "varargin{:}"); + +c = sprintf(strcat("function df = %s (%s)\n",\ + " ## Numerical differentiation of '%s' wrt to it's %d'th argument\n",\ + " ## This function has been written by 'cdiff()'\n",\ + " dx = %e;\n",\ + " sz = size (dvar);\n",\ + " ps = prod (sz);\n",\ + "%s",\ + "endfunction\n"),\ + dfunc,argstr,\ + func,wrt,\ + dx,\ + calstr) ; +