From c87b8fa4653cc28e016b29d296fbf50388f3e7c7 Mon Sep 17 00:00:00 2001
From: Houtan Bastani <houtan.bastani@ens.fr>
Date: Tue, 31 Jan 2012 10:19:49 +0100
Subject: [PATCH] remove files of the form: *.m[012]
---
MatlabFiles/CSMINWEL.M0 | 271 ----------------------------------------
MatlabFiles/ERRORS.m1 | 91 --------------
MatlabFiles/HESSCD.M1 | 55 --------
MatlabFiles/IMPULSE.m1 | 58 ---------
MatlabFiles/bfgsi.m0 | 24 ----
MatlabFiles/bfgsi.m1 | 24 ----
MatlabFiles/csminit.m0 | 163 ------------------------
MatlabFiles/csminit.m1 | 147 ----------------------
MatlabFiles/csminit.m2 | 158 -----------------------
MatlabFiles/csminwel.m1 | 206 ------------------------------
MatlabFiles/csminwel.m2 | 271 ----------------------------------------
MatlabFiles/gradcd.m1 | 53 --------
12 files changed, 1521 deletions(-)
delete mode 100755 MatlabFiles/CSMINWEL.M0
delete mode 100755 MatlabFiles/ERRORS.m1
delete mode 100755 MatlabFiles/HESSCD.M1
delete mode 100755 MatlabFiles/IMPULSE.m1
delete mode 100755 MatlabFiles/bfgsi.m0
delete mode 100755 MatlabFiles/bfgsi.m1
delete mode 100755 MatlabFiles/csminit.m0
delete mode 100755 MatlabFiles/csminit.m1
delete mode 100755 MatlabFiles/csminit.m2
delete mode 100755 MatlabFiles/csminwel.m1
delete mode 100755 MatlabFiles/csminwel.m2
delete mode 100755 MatlabFiles/gradcd.m1
diff --git a/MatlabFiles/CSMINWEL.M0 b/MatlabFiles/CSMINWEL.M0
deleted file mode 100755
index d3e8005..0000000
--- a/MatlabFiles/CSMINWEL.M0
+++ /dev/null
@@ -1,271 +0,0 @@
-function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel(fcn,x0,H0,grad,crit,nit, ...
- P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15,P16,P17,...
- P18,P19,P20)
-% [fhat,xhat,ghat,Hhat,itct,fcount,retcodehat] = csminwel(fcn,x0,H0,grad,crit,nit, ...
-% P1,P2,P3,P4)
-%
-% x0(estimated parameter) SHOULD BE A ROW VECTOR.
-%
-% fcn and grad are strings naming functions. If grad is null, numerical
-% derivatives are used. The P parameters need not be present. They are passed
-% to fcn as additional arguments.
-%
-% Reindented and modified to allow compilation (comment out use of tailstr, eval) by CAS,
-% 7/22/96
-%
-[nx,no]=size(x0);
-nx=max(nx,no);
-Verbose=1;
-NumGrad= ( ~isstr(grad) | length(grad)==0);
-done=0;
-itct=0;
-fcount=0;
-snit=100;
-tailstr = ')';
-stailstr = [];
-% Lines below make the number of Pi's optional. This is inefficient, though, and precludes
-% use of the matlab compiler. Without them, we use feval and the number of Pi's must be
-% changed with the editor for each application. Places where this is required are marked
-% with ARGLIST comments
-%for i=nargin-6:-1:1
-% tailstr=[ ',P' num2str(i) tailstr];
-% stailstr=[' P' num2str(i) stailstr];
-%end
-%f0 = eval([fcn '(x0' tailstr]);
-%ARGLIST
-f0 = feval(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
-% disp('first fcn in csminwel.m ----------------') % Jinill on 9/5/95
-if f0 > 1e50, disp('Bad initial parameter.'), return, end
-if NumGrad
- if length(grad)==0
- %[g badg] = eval(['numgrad(fcn, x0' tailstr]);
- %ARGLIST
- [g badg] = numgrad(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
- else
- badg=any(find(grad==0));
- g=grad;
- end
- %numgrad(fcn,x0,P1,P2,P3,P4);
-else
- %[g badg] = eval([grad '(x0' tailstr]);
- %ARGLIST
- [g badg] = feval(grad,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
-end
-x=x0;
-f=f0;
-H=H0;
-cliff=0;
-while ~done
- g1=[]; g2=[]; g3=[];
- %addition fj. 7/6/94 for control
- disp('-----------------')
- disp('-----------------')
- %disp('f and x at the beginning of new iteration')
- disp('f at the beginning of new iteration')
- f
- %if size(x,1)>size(x,2)
- % x'
- %else
- % x
- %end
- %-------------------------
- itct=itct+1;
- disp(' going to start csminit.m to produce f1')
- %[f1 x1 fc retcode1] = eval(['csminit(fcn,x,f,g,badg,H' tailstr]);
- %ARGLIST
- [f1 x1 fc retcode1] = csminit(fcn,x,f,g,badg,H,P1,P2,P3,P4,P5,P6,P7,...
- P8,P9,P10,P11,P12,P13);
- % itct=itct+1;
- fcount = fcount+fc;
- % erased on 8/4/94
- % if (retcode == 1) | (abs(f1-f) < crit)
- % done=1;
- % end
- % if itct > nit
- % done = 1;
- % retcode = -retcode;
- % end
- if retcode1 ~= 1
- if retcode1==2 | retcode1==4
- wall1=1; badg1=1;
- else
- if NumGrad
- %[g1 badg1] = eval(['numgrad(fcn, x1' tailstr]);
- %ARGLIST
- [g1 badg1] = numgrad(fcn, x1,P1,P2,P3,P4,P5,P6,P7,P8,P9,...
- P10,P11,P12,P13);
- else
- %[g1 badg1] = eval([grad '(x1' tailstr]);
- %ARGLIST
- [g1 badg1] = feval(grad, x1,P1,P2,P3,P4,P5,P6,P7,P8,P9,...
- P10,P11,P12,P13);
- end
- wall1=badg1;
- % g1
- badg1;
- %eval(['save g1 g1 x1 f1 ' stailstr]);
- %ARGLIST
- save g1 g1 x1 f1 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
- end
- if wall1 % & (~done) by Jinill
- % Bad gradient or back and forth on step length. Possibly at
- % cliff edge. Try perturbing search direction.
- %
- %fcliff=fh;xcliff=xh;
- Hcliff=H+diag(diag(H).*rand(nx,1));
- disp('Cliff. Perturbing search direction.')
- %[f2 x2 fc retcode2] = eval(['csminit(fcn,x,f,g,badg,Hcliff' tailstr]);
- %ARGLIST
- [f2 x2 fc retcode2] = csminit(fcn,x,f,g,badg,Hcliff,P1,P2,P3,P4,...
- P5,P6,P7,P8,P9,P10,P11,P12,P13);
- fcount = fcount+fc; % put by Jinill
- if f2 < f
- if retcode2==2 | retcode2==4
- wall2=1; badg2=1;
- else
- if NumGrad
- %[g2 badg2] = eval(['numgrad(fcn, x2' tailstr]);
- %ARGLIST
- [g2 badg2] = numgrad(fcn, x2,P1,P2,P3,P4,P5,P6,P7,P8,...
- P9,P10,P11,P12,P13);
- else
- %[g2 badg2] = eval([grad '(x2' tailstr]);
- %ARGLIST
- [g2 badg2] = feval(grad,x2,P1,P2,P3,P4,P5,P6,P7,P8,...
- P9,P10,P11,P12,P13);
- end
- wall2=badg2;
- % g2
- badg2
- %eval(['save g2 g2 x2 f2 ' stailstr]);
- %ARGLIST
- save g2 g2 x2 f2 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
- end
- if wall2
- disp('Cliff again. Try traversing')
- if norm(x2-x1) < 1e-13
- f3=f; x3=x; badg3=1;
- else
- gcliff=((f2-f1)/((norm(x2-x1))^2))*(x2-x1)';
- %[f3 x3 fc retcode3] = eval(['csminit(fcn,x,f,gcliff,0,eye(nx)' tailstr]);
- %ARGLIST
- [f3 x3 fc retcode3] = csminit(fcn,x,f,gcliff,0,eye(nx),P1,P2,P3,...
- P4,P5,P6,P7,P8,...
- P9,P10,P11,P12,P13);
- fcount = fcount+fc; % put by Jinill
- if retcode3==2 | retcode3==4
- wall3=1; badg3=1;
- else
- if NumGrad
- %[g3 badg3] = eval(['numgrad(fcn, x3' tailstr]);
- %ARGLIST
- [g3 badg3] = numgrad(fcn, x3,P1,P2,P3,P4,P5,P6,P7,P8,...
- P9,P10,P11,P12,P13);
- else
- %[g3 badg3] = eval([grad '(x3' tailstr]);
- %ARGLIST
- [g3 badg3] = feval(grad,x3,P1,P2,P3,P4,P5,P6,P7,P8,...
- P9,P10,P11,P12,P13);
- end
- wall3=badg3;
- % g3
- badg3
- %eval(['save g3 g3 x3 f3 ' stailstr]);
- %ARGLIST
- save g3 g3 x3 f3 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
- end
- end
- else
- f3=f; x3=x; badg3=1;
- end
- else
- f3=f; x3=x; badg3=1;
- end
- else
- % normal iteration, no walls, or else we're finished here.
- f2=f; f3=f; badg2=1; badg3=1;
- end
- end
- %how to pick gh and xh
- if f3<f & badg3==0
- ih=3
- fh=f3;xh=x3;gh=g3;badgh=badg3;retcodeh=retcode3;
- elseif f2<f & badg2==0
- ih=2
- fh=f2;xh=x2;gh=g2;badgh=badg2;retcodeh=retcode2;
- elseif f1<f & badg1==0
- ih=1
- fh=f1;xh=x1;gh=g1;badgh=badg1;retcodeh=retcode1;
- else
- [fh,ih] = min([f1,f2,f3]);
- ih
- %eval(['xh=x' num2str(ih) ';'])
- if size(x1,2)>1
- xi=[x1;x2;x3];
- xh=xi(ih,:);
- else
- xi=[x1 x2 x3];
- xh=xi(:,ih);
- end
- %eval(['gh=g' num2str(ih) ';'])
- %eval(['retcodeh=retcode' num2str(ih) ';'])
- retcodei=[retcode1,retcode2,rectode3];
- retcodeh=retcodei(ih);
- if gh==[]
- if NumGrad
- %[gh badgh] = eval(['numgrad(fcn, xh' tailstr]);
- %ARGLIST
- [gh badgh] = numgrad(fcn,xh,P1,P2,P3,P4,P5,P6,P7,P8,...
- P9,P10,P11,P12,P13);
- else
- %[gh badgh] = eval([grad '(xh' tailstr]);
- %ARGLIST
- [gh badgh] = feval(grad,xh,P1,P2,P3,P4,P5,P6,P7,P8,...
- P9,P10,P11,P12,P13);
- end
- end
- badgh=1;
- end
- %end of picking
- %ih
- %fh
- %xh
- %gh
- %badgh
- stuck = (abs(fh-f) < crit);
- if (~badg)&(~badgh)&(~stuck)
- H = bfgsi(H,gh-g,xh-x);
- end
- if Verbose
- disp('----')
- disp(['Improvement on iteration ' int2str(itct) ' = ' num2str(f-fh)])
- end
- if Verbose
- if itct > nit
- disp('iteration count termination')
- done = 1;
- elseif stuck
- disp('improvement < crit termination')
- done = 1;
- end
- rc=retcodeh;
- if rc == 1
- disp('zero gradient')
- elseif rc == 6
- disp('smallest step still improving too slow, reversed gradient')
- elseif rc == 5
- disp('largest step still improving too fast')
- elseif (rc == 4) | (rc==2)
- disp('back and forth on step length never finished')
- elseif rc == 3
- disp('smallest step still improving too slow')
- end
- end
- f=fh;
- x=xh;
- g=gh;
- badg=badgh;
-end
-% what about making an m-file of 10 lines including numgrad.m
-% since it appears three times in csminwel.m
-
\ No newline at end of file
diff --git a/MatlabFiles/ERRORS.m1 b/MatlabFiles/ERRORS.m1
deleted file mode 100755
index 19e06be..0000000
--- a/MatlabFiles/ERRORS.m1
+++ /dev/null
@@ -1,91 +0,0 @@
-function [vd,str,imf] = errors(Bh,swish,nn)
-% Computing variance decompositions and impulse functions with
-% [vd,str,imf] = errors(Bh,swish,nn)
-% where imf and vd is of the same format as in RATS, that is to say:
-% Column: nvar responses to 1st shock,
-% nvar responses to 2nd shock, and so on.
-% Row: steps of impulse responses.
-% vd: variance decompositions
-% str: standard errors of each variable, steps-by-nvar
-% imf: impulse response functions
-% Bh is the estimated reduced form coefficient in the form
-% Y(T*nvar) = XB + U, X: T*k, B: k*nvar. The matrix
-% form or dimension is the same as "Bh" from the function "sye";
-% swish is the inv(A0) in the structural model A(L)y(t) = e(t).
-% nn is the numbers of inputs [nvar,lags,# of impulse responses].
-
-nvar = nn(1);
-lags = nn(2);
-imstep = nn(3); % number of steps for impulse responses
-
-Ah = Bh';
-% Row: nvar equations
-% Column: 1st lag (with nvar variables) to lags (with nvar variables) + const = k.
-
-imf = zeros(imstep,nvar*nvar);
-vd = imf;
-% Column: nvar responses to 1st shock, nvar responses to 2nd shock, and so on.
-% Row: steps of impulse responses.
-str = zeros(imstep,nvar); % initializing standard errors of each equation
-M = zeros(nvar*(lags+1),nvar);
-% Stack M0;M1;M2;...;Mlags
-M(1:nvar,:) = swish;
-Mtem = M(1:nvar,:); % temporary M -- impulse responses.
-%
-Mvd = Mtem.^2; % saved for the cumulative sum later
-Mvds = (sum(Mvd'))';
-str(1,:) = sqrt(Mvds'); % standard errors of each equation
-Mvds = Mvds(:,ones(size(Mvds,1),1));
-Mvdtem = (100.0*Mvd) ./ Mvds; % tempoary Mvd -- variance decompositions
-% first or initial responses to
-% one standard deviation shock (or forecast errors).
-% Row: responses; Column: shocks
-%
-% * put in the form of "imf"
-imf(1,:) = Mtem(:)';
-vd(1,:) = Mvdtem(:)';
-
-t = 1;
-ims1 = min([imstep-1 lags]);
-while t <= ims1
- Mtem = zeros(nvar,nvar);
- for k = 1:t
- Mtem = Ah(:,nvar*(k-1)+1:nvar*k)*M(nvar*(t-k)+1:nvar*(t-k+1),:) + Mtem;
- % Row: nvar equations, each for the nvar variables at tth lag
- end
- % ** impulse response functions
- M(nvar*t+1:nvar*(t+1),:) = Mtem;
- imf(t+1,:) = Mtem(:)';
- % stack imf with each step, Row: 6 var to 1st shock, 6 var to 2nd shock, etc.
- % ** variance decompositions
- Mvd = Mvd + Mtem.^2; % saved for the cumulative sum later
- Mvds = (sum(Mvd'))';
- str(t+1,:) = sqrt(Mvds'); % standard errors of each equation
- Mvds = Mvds(:,ones(size(Mvds,1),1));
- Mvdtem = (100.0*Mvd) ./ Mvds; % tempoary Mvd -- variance decompositions
- vd(t+1,:) = Mvdtem(:)';
- % stack vd with each step, Row: 6 var to 1st shock, 6 var to 2nd shock, etc.
- t= t+1;
-end
-
-for t = lags+1:imstep-1
- M(1:nvar*lags,:) = M(nvar+1:nvar*(lags+1),:);
- Mtem = zeros(nvar,nvar);
- for k = 1:lags
- Mtem = Ah(:,nvar*(k-1)+1:nvar*k)*M(nvar*(lags-k)+1:nvar*(lags-k+1),:) + Mtem;
- % Row: nvar equations, each for the nvar variables at tth lag
- end
- % ** impulse response functions
- M(nvar*lags+1:nvar*(lags+1),:) = Mtem;
- imf(t+1,:) = Mtem(:)';
- % stack imf with each step, Row: 6 var to 1st shock, 6 var to 2nd shock, etc.
- % ** variance decompositions
- Mvd = Mvd + Mtem.^2; % saved for the cumulative sum later
- Mvds = (sum(Mvd'))';
- str(t+1,:) = sqrt(Mvds'); % standard errors of each equation
- Mvds = Mvds(:,ones(size(Mvds,1),1));
- Mvdtem = (100.0*Mvd) ./ Mvds; % tempoary Mvd -- variance decompositions
- vd(t+1,:) = Mvdtem(:)';
- % stack vd with each step, Row: 6 var to 1st shock, 6 var to 2nd shock, etc.
-end
-
\ No newline at end of file
diff --git a/MatlabFiles/HESSCD.M1 b/MatlabFiles/HESSCD.M1
deleted file mode 100755
index 017db98..0000000
--- a/MatlabFiles/HESSCD.M1
+++ /dev/null
@@ -1,55 +0,0 @@
-function grdd = hesscd(fcn,x0,grdh)
-% computing numerical hessian using a central difference with
-% function grdd = hesscd(fcn,x0,grdh)
-%
-% fcn: a string naming the objective function.
-% x0: a column vector n*1, at which point the hessian is evaluated.
-% grdh: step size.
-% grdd: hessian matrix (second derivative), n*n.
-
-stps = eps^(1/3);
-% eps: floating point relative accuracy or machine precision: 2.22e-16
-% stps: step size recommended by Dennis and Schnabel: 6.006e-6
-
-% ** initializations
-k = size(x0,1);
-grdd = zeros(k);
-
-% ** Computation of stepsize (dh)
-if all(grdh)
- dh = grdh;
-else
- ax0 = abs(x0);
- if all(x0)
- dax0 = x0 ./ ax0;
- else
- dax0 = 1;
- end
- dh = stps * (max([ax0 (1e-2)*ones(k,1)]'))' .* dax0;
-end
-
-xdh = x0 + dh;
-dh = xdh - x0; % This increases precision slightly
-dhm = dh(:,ones(k,1));
-ee = eye(k) .* dhm;
-
-i = 1;
-while i <= k
- j = i;
- while j <= k
-
- grdd(i,j) = (feval(fcn,x0 + ee(:,i) + ee(:,j)) - ...
- feval(fcn,x0 - ee(:,i) + ee(:,j)) - ...
- feval(fcn,x0 + ee(:,i) - ee(:,j)) + ...
- feval(fcn,x0 - ee(:,i) - ee(:,j))) ...
- / (4 * dh(i) * dh(j));
-
- if i ~= j
- grdd(j,i) = grdd(i,j);
- end
-
- j = j+1;
- end
- i = i+1;
-end
-
diff --git a/MatlabFiles/IMPULSE.m1 b/MatlabFiles/IMPULSE.m1
deleted file mode 100755
index bc8490a..0000000
--- a/MatlabFiles/IMPULSE.m1
+++ /dev/null
@@ -1,58 +0,0 @@
-function imf = impulse(Bh,swish,nn)
-% Computing impulse functions with
-% imf = impulse(Bh,swish,nn)
-% where imf is in a format that is the SAME as in RATS.
-% Column: nvar responses to 1st shock,
-% nvar responses to 2nd shock, and so on.
-% Row: steps of impulse responses.
-% Bh is the estimated reduced form coefficient in the form
-% Y(T*nvar) = XB + U, X: T*k, B: k*nvar. The matrix
-% form or dimension is the same as "Bh" from the function "sye";
-% swish is the inv(A0) in the structural model A(L)y(t) = e(t).
-% nn is the numbers of inputs [nvar,lags,# of impulse responses].
-
-nvar = nn(1);
-lags = nn(2);
-imstep = nn(3); % number of steps for impulse responses
-
-Ah = Bh';
-% Row: nvar equations
-% Column: 1st lag (with nvar variables) to lags (with nvar variables) + const = k.
-
-imf = zeros(imstep,nvar*nvar);
-% Column: nvar responses to 1st shock, nvar responses to 2nd shock, and so on.
-% Row: steps of impulse responses.
-M = zeros(nvar*(lags+1),nvar);
-% Stack M0;M1;M2;...;Mlags
-M(1:nvar,:) = swish;
-Mtem = M(1:nvar,:); % temporary M.
-% first (initial) responses to 1 standard deviation shock. Row: responses; Column: shocks
-% * put in the form of "imf"
-imf(1,:) = Mtem(:)';
-
-t = 1;
-ims1 = min([imstep-1 lags]);
-while t <= ims1
- Mtem = zeros(nvar,nvar);
- for k = 1:t
- Mtem = Ah(:,nvar*(k-1)+1:nvar*k)*M(nvar*(t-k)+1:nvar*(t-k+1),:) + Mtem;
- % Row: nvar equations, each for the nvar variables at tth lag
- end
- M(nvar*t+1:nvar*(t+1),:) = Mtem;
- imf(t+1,:) = Mtem(:)';
- % stack imf with each step, Row: 6 var to 1st shock, 6 var to 2nd shock, etc.
- t= t+1;
-end
-
-for t = lags+1:imstep-1
- M(1:nvar*lags,:) = M(nvar+1:nvar*(lags+1),:);
- Mtem = zeros(nvar,nvar);
- for k = 1:lags
- Mtem = Ah(:,nvar*(k-1)+1:nvar*k)*M(nvar*(lags-k)+1:nvar*(lags-k+1),:) + Mtem;
- % Row: nvar equations, each for the nvar variables at tth lag
- end
- M(nvar*lags+1:nvar*(lags+1),:) = Mtem;
- imf(t+1,:) = Mtem(:)';
- % stack imf with each step, Row: 6 var to 1st shock, 6 var to 2nd shock, etc.
-end
-
\ No newline at end of file
diff --git a/MatlabFiles/bfgsi.m0 b/MatlabFiles/bfgsi.m0
deleted file mode 100755
index 800fefa..0000000
--- a/MatlabFiles/bfgsi.m0
+++ /dev/null
@@ -1,24 +0,0 @@
-function H = bfgsi(H0,dg,dx)
-% H = bfgsi(H0,dg,dx)
-% dg is previous change in gradient; dx is previous change in x;
-% 6/8/93 version that updates inverse hessian instead of hessian
-% itself.
-if size(dg,2)>1
- dg=dg';
-end
-if size(dx,2)>1
- dx=dx';
-end
-Hdg = H0*dg;
-dgdx = dg'*dx;
-if (abs(dgdx) >1e-12)
- H = H0 + (1+(dg'*Hdg)/dgdx)*(dx*dx')/dgdx - (dx*Hdg'+Hdg*dx')/dgdx;
-else
- disp('bfgs update failed.')
- disp(['|dg| = ' num2str(sqrt(dg'*dg)) '|dx| = ' num2str(sqrt(dx'*dx))]);
- disp(['dg''*dx = ' num2str(dgdx)])
- disp(['|H*dg| = ' num2str(Hdg'*Hdg)])
- H=H0;
-end
-save H.dat H;
-
\ No newline at end of file
diff --git a/MatlabFiles/bfgsi.m1 b/MatlabFiles/bfgsi.m1
deleted file mode 100755
index 8090b32..0000000
--- a/MatlabFiles/bfgsi.m1
+++ /dev/null
@@ -1,24 +0,0 @@
-function H = bfgsi(H0,dg,dx)
-% H = bfgsi(H0,dg,dx)
-% dg is previous change in gradient; dx is previous change in x;
-% 6/8/93 version that updates inverse hessian instead of hessian
-% itself.
-if size(dg,2)>1
- dg=dg';
-end
-if size(dx,2)>1
- dx=dx';
-end
-Hdg = H0*dg;
-dgdx = dg'*dx;
-if (abs(dgdx) >1e-12)
- H = H0 + (1+(dg'*Hdg)/dgdx)*(dx*dx')/dgdx - (dx*Hdg'+Hdg*dx')/dgdx;
-else
- disp('bfgs update failed.')
- disp(['|dg| = ' num2str(sqrt(dg'*dg)) '|dx| = ' num2str(sqrt(dx'*dx))]);
- disp(['dg''*dx = ' num2str(dgdx)])
- disp(['|H*dg| = ' num2str(Hdg'*Hdg)])
- H=H0;
-end
-save H.dat H
-
diff --git a/MatlabFiles/csminit.m0 b/MatlabFiles/csminit.m0
deleted file mode 100755
index 1af99f8..0000000
--- a/MatlabFiles/csminit.m0
+++ /dev/null
@@ -1,163 +0,0 @@
-function [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,...
- P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15,P16,P17,...
- P18,P19,P20)
-% [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,...
-% P1,P2,P3,P4,P5,P6,P7,P8)
-% retcodes: 0, normal step. 5, largest step still improves too fast.
-% 4,2 back and forth adjustment of stepsize didn't finish. 3, smallest
-% stepsize still improves too slow. 6, no improvement found. 1, zero
-% gradient.
-%---------------------
-% Modified 7/22/96 to omit variable-length P list, for efficiency and compilation.
-% Places where the number of P's need to be altered or the code could be returned to
-% its old form are marked with ARGLIST comments.
-%
-% Fixed 7/17/93 to use inverse-hessian instead of hessian itself in bfgs
-% update.
-%
-% Fixed 7/19/93 to flip eigenvalues of H to get better performance when
-% it's not psd.
-%
-tailstr = ')';
-for i=nargin-6:-1:1
- tailstr=[ ',P' num2str(i) tailstr];
-end
-%ANGLE = .03;
-ANGLE = .005;
-%THETA = .03;
-THETA = .1;
-FCHANGE = 1000;
-MINLAMB = 1e-9;
-% fixed 7/15/94
-% MINDX = .0001;
-% MINDX = 1e-6;
-MINDFAC = .01;
-fcount=0;
-lambda=1;
-lambdahat=1;
-xhat=x0;
-f=f0;
-fhat=f0;
-g = g0;
-gnorm = norm(g);
-%
-if (gnorm < 1.e-12) & ~badg % put ~badg 8/4/94
- retcode =1;
- % gradient convergence
-else
- % with badg true, we don't try to match rate of improvement to directional
- % derivative. We're satisfied just to get some improvement in f.
- %
- %if(badg)
- % dx = -g*FCHANGE/(gnorm*gnorm);
- % dxnorm = norm(dx);
- % if dxnorm > 1e12
- % disp('Bad, small gradient problem.')
- % dx = dx*FCHANGE/dxnorm;
- % end
- %else
- % Gauss-Newton step;
- %---------- Start of 7/19/93 mod ---------------
- %[v d] = eig(H0);
- %toc
- %d=max(1e-10,abs(diag(d)));
- %d=abs(diag(d));
- %dx = -(v.*(ones(size(v,1),1)*d'))*(v'*g);
-% toc
- dx = -H0*g;
-% toc
- dxnorm = norm(dx);
- if dxnorm > 1e12
- disp('Near-singular H problem.')
- dx = dx*FCHANGE/dxnorm;
- end
- dfhat = dx'*g0;
- %end
- %
- %
- if ~badg
- % test for alignment of dx with gradient and fix if necessary
- a = -dfhat/(gnorm*dxnorm);
- if a<ANGLE
- dx = dx - (ANGLE*dxnorm/gnorm+dfhat/(gnorm*gnorm))*g;
- dfhat = dx'*g;
- dxnorm = norm(dx);
- disp(sprintf('Correct for low angle: %g',a))
- end
- end
- end
- disp(sprintf('Predicted improvement: %18.9f',-dfhat/2))
- %
- % Have OK dx, now adjust length of step (lambda) until min and
- % max improvement rate criteria are met.
- done=0;
- factor=3;
- shrink=1;
- while ~done
- if size(x0,2)>1
- dxtest=x0+dx'*lambda;
- else
- dxtest=x0+dx*lambda;
- end
- % home
- f = eval([fcn '(dxtest' tailstr]);
- %ARGLIST
- %f = feval(fcn,dxtest,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
- % f = feval(fcn,x0+dx*lambda,P1,P2,P3,P4,P5,P6,P7,P8);
- disp(sprintf('lambda = %10.5g; f = %20.7f',lambda,f ))
- %debug
- %disp(sprintf('Improvement too great? f0-f: %g, criterion: %g',f0-f,-(1-THETA)*dfhat*lambda))
- if f<fhat
- fhat=f;
- xhat=dxtest;
- lambdahat = lambda;
- end
- fcount=fcount+1;
- if (~badg & (f0-f < max([-THETA*dfhat*lambda 0]))) | ...
- (badg & (f0-f) < 0)
- if ~shrink
- shrink=1;
- factor=factor^.6;
- %if (abs(lambda)*(factor-1)*dxnorm < MINDX) | (abs(lambda)*(factor-1) < MINLAMB)
- if abs(factor-1)<MINDFAC
- retcode=2;
- done=1;
- end
- end
- lambda=lambda/factor;
- if abs(lambda) < MINLAMB
- if (lambda > 0) & (f0 <= fhat)
- % try going against gradient, which may be inaccurate
- lambda = -lambda*factor^6
- else
- if lambda < 0
- retcode = 6;
- else
- retcode = 3;
- end
- done = 1;
- end
- end
- elseif ~badg & ( (lambda > 0) & (f0-f > -(1-THETA)*dfhat*lambda) )
- if shrink
- shrink=0;
- factor = factor^.6;
- %if ( abs(lambda)*(factor-1)*dxnorm< MINDX ) | ( abs(lambda)*(factor-1)< MINLAMB)
- if abs(factor-1) < MINDFAC
- retcode = 4;
- done=1;
- end
- end
- lambda=lambda*factor;
- if abs(lambda) > 1e20;
- retcode = 5;
- done =1;
- end
- else
- done=1;
- retcode=0;
- end
- end
-end
-disp(sprintf('Norm of dx %10.5g', dxnorm))
-
\ No newline at end of file
diff --git a/MatlabFiles/csminit.m1 b/MatlabFiles/csminit.m1
deleted file mode 100755
index 6990dc0..0000000
--- a/MatlabFiles/csminit.m1
+++ /dev/null
@@ -1,147 +0,0 @@
-function [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,...
- P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15,P16,P17,...
- P18,P19,P20)
-% [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,...
-% P1,P2,P3,P4,P5,P6,P7,P8)
-% retcodes: 0, normal step. 5, largest step still improves too fast.
-% 4,2 back and forth adjustment of stepsize didn't finish. 3, smallest
-% stepsize still improves too slow. 6, no improvement found. 1, zero
-% gradient.
-%---------------------
-% Fixed 7/17/93 to use inverse-hessian instead of hessian itself in bfgs
-% update.
-%
-% Fixed 7/19/93 to flip eigenvalues of H to get better performance when
-% it's not psd.
-%
-tailstr = ')';
-for i=nargin-6:-1:1
- tailstr=[ ',P' num2str(i) tailstr];
-end
-ANGLE = .03;
-%%THETA = .03;
-THETA = 0.1; % TZ, mod 8/11/96
-FCHANGE = 1000;
-MINLAMB = 1e-9;
-% fixed 7/15/94
-% MINDX = .0001;
-MINDX = 1e-6;
-fcount=0;
-lambda=1;
-lambdahat=1;
-xhat=x0;
-f=f0;
-fhat=f0;
-g = g0;
-gnorm = norm(g);
-%
-if (gnorm < 1.e-12) & ~badg % put ~badg 8/4/94
- retcode =1;
- % gradient convergence
-else
- % with badg true, we don't try to match rate of improvement to directional
- % derivative. We're satisfied just to get some improvement in f.
- %
- %if(badg)
- % dx = -g*FCHANGE/(gnorm*gnorm);
- % dxnorm = norm(dx);
- % if dxnorm > 1e12
- % disp('Bad, small gradient problem.')
- % dx = dx*FCHANGE/dxnorm;
- % end
- %else
- % Gauss-Newton step;
- %---------- Start of 7/19/93 mod ---------------
- [v d] = eig(H0);
- %% 1e5*diag(d)'; % Tao Zha
- d=max(1e-10,abs(diag(d)));
- dx = -v*diag(d)*v'*g;
- %dx = -H0*g;
- dxnorm = norm(dx);
- if dxnorm > 1e12
- disp('Near-singular H problem.')
- dx = dx*FCHANGE/dxnorm;
- end
- dfhat = dx'*g0;
- %end
- %
- %
- if ~badg
- % test for alignment of dx with gradient and fix if necessary
- a = -dfhat/(gnorm*dxnorm);
- if a<ANGLE
- dx = dx - (ANGLE*dxnorm/gnorm+dfhat/(gnorm*gnorm))*g;
- dfhat = dx'*g;
- dxnorm = norm(dx);
- end
- end
- %
- % Have OK dx, now adjust length of step (lambda) until min and
- % max improvement rate criteria are met.
- done=0;
- factor=3;
- shrink=1;
- while ~done
- if size(x0,2)>1
- dxtest=x0+dx'*lambda;
- else
- dxtest=x0+dx*lambda;
- end
- % home
- f = eval([fcn '(dxtest' tailstr]);
- % f = feval(fcn,x0+dx*lambda,P1,P2,P3,P4,P5,P6,P7,P8);
- fprintf('lambda = %12g; f = %12g\n',lambda,f )
- if f<fhat
- fhat=f;
- xhat=dxtest;
- lambdahat = lambda;
- end
- fcount=fcount+1;
- if (~badg & (f0-f < max([-THETA*dfhat*lambda 0]))) | ...
- (badg & (f0-f) < 0)
- if ~shrink
- shrink=1;
- factor=factor^.6;
- %if abs(lambda)*(factor-1)*dxnorm < MINDX
- if abs(lambda)*(factor-1) < MINLAMB
- retcode=2;
- done=1;
- end
- end
- lambda=lambda/factor;
- %if abs(lambda) < MINLAMB%
- if abs(lambda)*dxnorm < MINDX
- if (lambda > 0) & (f0 <= fhat)
- % try going against gradient, which may be inaccurate
- lambda = -lambda*factor^6
- else
- if lambda < 0
- retcode = 6;
- else
- retcode = 3;
- end
- done = 1;
- end
- end
- elseif ~badg & ( (lambda > 0) & (f0-f > -(1-THETA)*dfhat*lambda) )
- if shrink
- shrink=0;
- factor = factor^.6;
- %if abs(lambda)*(factor-1)*dxnorm< MINDX
- if abs(lambda)*(factor-1)< MINLAMB
- retcode = 4;
- done=1;
- end
- end
- lambda=lambda*factor;
- if abs(lambda) > 1e20;
- retcode = 5;
- done =1;
- end
- else
- done=1;
- retcode=0;
- end
- end
-end
-
\ No newline at end of file
diff --git a/MatlabFiles/csminit.m2 b/MatlabFiles/csminit.m2
deleted file mode 100755
index eabf9af..0000000
--- a/MatlabFiles/csminit.m2
+++ /dev/null
@@ -1,158 +0,0 @@
-function [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,...
- P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15,P16,P17,...
- P18,P19,P20)
-% [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,...
-% P1,P2,P3,P4,P5,P6,P7,P8)
-% retcodes: 0, normal step. 5, largest step still improves too fast.
-% 4,2 back and forth adjustment of stepsize didn't finish. 3, smallest
-% stepsize still improves too slow. 6, no improvement found. 1, zero
-% gradient.
-%---------------------
-% Modified 7/22/96 to omit variable-length P list, for efficiency and compilation.
-% Places where the number of P's need to be altered or the code could be returned to
-% its old form are marked with ARGLIST comments.
-%
-% Fixed 7/17/93 to use inverse-hessian instead of hessian itself in bfgs
-% update.
-%
-% Fixed 7/19/93 to flip eigenvalues of H to get better performance when
-% it's not psd.
-%
-tailstr = ')';
-for i=nargin-6:-1:1
- tailstr=[ ',P' num2str(i) tailstr];
-end
-ANGLE = 0.005; % mod 8/12/96
-%ANGLE = .03;
-THETA = 0.1; % mod 8/12/96
-%THETA = .03;
-FCHANGE = 1000;
-MINLAMB = 1e-9;
-% fixed 7/15/94
-% MINDX = .0001;
-MINDX = 1e-6;
-fcount=0;
-lambda=1;
-lambdahat=1;
-xhat=x0;
-f=f0;
-fhat=f0;
-g = g0;
-gnorm = norm(g);
-%
-if (gnorm < 1.e-12) & ~badg % put ~badg 8/4/94
- retcode =1;
- % gradient convergence
-else
- % with badg true, we don't try to match rate of improvement to directional
- % derivative. We're satisfied just to get some improvement in f.
- %
- %if(badg)
- % dx = -g*FCHANGE/(gnorm*gnorm);
- % dxnorm = norm(dx);
- % if dxnorm > 1e12
- % disp('Bad, small gradient problem.')
- % dx = dx*FCHANGE/dxnorm;
- % end
- %else
- % Gauss-Newton step;
- %---------- Start of 7/19/93 mod ---------------
- %[v d] = eig(H0);
- %toc
- %d=max(1e-10,abs(diag(d)));
- %d=abs(diag(d));
- %dx = -(v.*(ones(size(v,1),1)*d'))*(v'*g);
-% toc
- dx = -H0*g;
-% toc
- dxnorm = norm(dx);
- if dxnorm > 1e12
- disp('Near-singular H problem.')
- dx = dx*FCHANGE/dxnorm;
- end
- dfhat = dx'*g0;
- %end
- %
- %
- if ~badg
- % test for alignment of dx with gradient and fix if necessary
- a = -dfhat/(gnorm*dxnorm);
- if a<ANGLE
- dx = dx - (ANGLE*dxnorm/gnorm+dfhat/(gnorm*gnorm))*g;
- dfhat = dx'*g;
- dxnorm = norm(dx);
- end
- end
- %
- % Have OK dx, now adjust length of step (lambda) until min and
- % max improvement rate criteria are met.
- done=0;
- factor=3;
- shrink=1;
- while ~done
- if size(x0,2)>1
- dxtest=x0+dx'*lambda;
- else
- dxtest=x0+dx*lambda;
- end
- % home
- f = eval([fcn '(dxtest' tailstr]);
- %ARGLIST
- %f = feval(fcn,dxtest,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
- % f = feval(fcn,x0+dx*lambda,P1,P2,P3,P4,P5,P6,P7,P8);
- fprintf('lambda = %12g; f = %12g\n',lambda,f )
- if f<fhat
- fhat=f;
- xhat=dxtest;
- lambdahat = lambda;
- end
- fcount=fcount+1;
- if (~badg & (f0-f < max([-THETA*dfhat*lambda 0]))) | ...
- (badg & (f0-f) < 0)
- if ~shrink
- shrink=1;
- factor=factor^.6;
- %if abs(lambda)*(factor-1)*dxnorm < MINDX
- if abs(lambda)*(factor-1) < MINLAMB
- retcode=2;
- done=1;
- end
- end
- lambda=lambda/factor;
- %if abs(lambda) < MINLAMB%
- if abs(lambda)*dxnorm < MINDX
- if (lambda > 0) & (f0 <= fhat)
- % try going against gradient, which may be inaccurate
- lambda = -lambda*factor^6
- else
- if lambda < 0
- retcode = 6;
- else
- retcode = 3;
- end
- done = 1;
- end
- end
- elseif ~badg & ( (lambda > 0) & (f0-f > -(1-THETA)*dfhat*lambda) )
- if shrink
- shrink=0;
- factor = factor^.6;
- %if abs(lambda)*(factor-1)*dxnorm< MINDX
- if abs(lambda)*(factor-1)< MINLAMB
- retcode = 4;
- done=1;
- end
- end
- lambda=lambda*factor;
- if abs(lambda) > 1e20;
- retcode = 5;
- done =1;
- end
- else
- done=1;
- retcode=0;
- end
- end
-end
-disp(sprintf('Norm of dx %d', dxnorm))
-
\ No newline at end of file
diff --git a/MatlabFiles/csminwel.m1 b/MatlabFiles/csminwel.m1
deleted file mode 100755
index a4c2315..0000000
--- a/MatlabFiles/csminwel.m1
+++ /dev/null
@@ -1,206 +0,0 @@
-function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel(fcn,x0,H0,grad,crit,nit, ...
- P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15,P16,P17,...
- P18,P19,P20)
-% [fhat,xhat,ghat,Hhat,itct,fcount,retcodehat] = csminwel(fcn,x0,H0,grad,crit,nit, ...
-% P1,P2,P3,P4,P5,P6,P7,P8)
-%
-% x0(estimated parameter) SHOULD BE A ROW VECTOR.
-%
-% fcn and grad are strings naming functions. If grad is null, numerical
-% derivatives are used. The P parameters need not be present. They are passed
-% to fcn as additional arguments.
-%
-%
-[nx,no]=size(x0);
-nx=max(nx,no);
-Verbose=1;
-NumGrad= ( ~isstr(grad) | length(grad)==0);
-done=0;
-itct=0;
-fcount=0;
-snit=100;
-tailstr = ')';
-stailstr = [];
-for i=nargin-6:-1:1
- tailstr=[ ',P' num2str(i) tailstr];
- stailstr=[' P' num2str(i) stailstr];
-end
-f0 = eval([fcn '(x0' tailstr]);
-% disp('first fcn in csminwel.m ----------------') % Jinill on 9/5/95
-% f0 = feval(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8);
-if f0 > 1e50, disp('Bad initial parameter.'), return, end
-if NumGrad
- if length(grad)==0
- [g badg] = eval(['numgrad(fcn, x0' tailstr]);
- else
- badg=any(find(grad==0));
- g=grad;
- end
- %numgrad(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8);
-else
- [g badg] = eval([grad '(x0' tailstr]);
- % feval(grad,x0,P1,P2,P3,P4,P5,P6,P7,P8);
-end
-x=x0;
-f=f0;
-H=H0;
-cliff=0;
-while ~done
-g1=[]; g2=[]; g3=[];
-%addition fj. 7/6/94 for control
-disp('-----------------')
-disp('-----------------')
-disp('f and x at the beginning of new iteration')
-f
-if size(x,1)>size(x,2)
- x'
-else
- x
-end
-%-------------------------
- itct=itct+1;
- disp(' going to start csminit.m to produce f1')
- [f1 x1 fc retcode1] = eval(['csminit(fcn,x,f,g,badg,H' tailstr]);
- % itct=itct+1;
- fcount = fcount+fc;
- % erased on 8/4/94
- % if (retcode == 1) | (abs(f1-f) < crit)
- % done=1;
- % end
- % if itct > nit
- % done = 1;
- % retcode = -retcode;
- % end
- if retcode1 ~= 1
- if retcode1==2 | retcode1==4
- wall1=1; badg1=1;
- else
- if NumGrad, [g1 badg1] = eval(['numgrad(fcn, x1' tailstr]);
- else, [g1 badg1] = eval([grad '(x1' tailstr]);
- end
- wall1=badg1;
- % g1
- badg1
- eval(['save hxg.dat g1 x1 f1 ' stailstr]);
- end
- if wall1 % & (~done) by Jinill
- % Bad gradient or back and forth on step length. Possibly at
- % cliff edge. Try perturbing search direction.
- %
- %fcliff=fh;xcliff=xh;
- Hcliff=H+diag(diag(H).*rand(nx,1));
- disp('Cliff. Perturbing search direction.')
- [f2 x2 fc retcode2] = eval(['csminit(fcn,x,f,g,badg,Hcliff' tailstr]);
- fcount = fcount+fc; % put by Jinill
- if f2 < f
- if retcode2==2 | retcode2==4
- wall2=1; badg2=1;
- else
- if NumGrad, [g2 badg2] = eval(['numgrad(fcn, x2' tailstr]);
- else, [g2 badg2] = eval([grad '(x2' tailstr]);
- end
- wall2=badg2;
- % g2
- badg2
- eval(['save g2 g2 x2 f2 ' stailstr]);
- end
- if wall2
- disp('Cliff again. Try traversing')
- if norm(x2-x1) < 1e-13
- f3=f; x3=x; badg3=1;
- else
- gcliff=((f2-f1)/((norm(x2-x1))^2))*(x2-x1)';
- [f3 x3 fc retcode3] = eval(['csminit(fcn,x,f,gcliff,0,eye(nx)' tailstr]);
- fcount = fcount+fc; % put by Jinill
- if retcode3==2 | retcode3==4
- wall3=1; badg3=1;
- else
- if NumGrad, [g3 badg3] = eval(['numgrad(fcn, x3' tailstr]);
- else, [g3 badg3] = eval([grad '(x3' tailstr]);
- end
- wall3=badg3;
- % g3
- badg3
- eval(['save g3 g3 x3 f3 ' stailstr]);
- end
- end
- else
- f3=f; x3=x; badg3=1;
- end
- else
- f3=f; x3=x; badg3=1;
- end
-
- else
- % normal iteration, no walls, or else we're finished here.
- f2=f; f3=f; badg2=1; badg3=1;
- end
- end
-
- %how to pick gh and xh
- if f3<f & badg3==0
- ih=3
- fh=f3;xh=x3;gh=g3;badgh=badg3;retcodeh=retcode3;
- elseif f2<f & badg2==0
- ih=2
- fh=f2;xh=x2;gh=g2;badgh=badg2;retcodeh=retcode2;
- elseif f1<f & badg1==0
- ih=1
- fh=f1;xh=x1;gh=g1;badgh=badg1;retcodeh=retcode1;
- else
- [fh,ih] = min([f1,f2,f3]);
- ih
- eval(['xh=x' num2str(ih) ';'])
- eval(['gh=g' num2str(ih) ';'])
- eval(['retcodeh=retcode' num2str(ih) ';'])
- if gh==[]
- if NumGrad, [gh badgh] = eval(['numgrad(fcn, xh' tailstr]);
- else, [gh badgh] = eval([grad '(xh' tailstr]);
- end
- end
- badgh=1;
- end
- %end of picking
- %ih
- %fh
- %xh
- %gh
- %badgh
-
- stuck = (abs(fh-f) < crit);
- if (~badg)&(~badgh)&(~stuck), H = bfgsi(H,gh-g,xh-x); end
- if Verbose
- disp('----')
- disp(['Improvement on iteration ' int2str(itct) ' = ' num2str(f-fh)])
- end
- if Verbose
- if itct > nit
- disp('iteration count termination')
- done = 1;
- elseif stuck
- disp('improvement < crit termination')
- done = 1;
- end
- rc=retcodeh;
- if rc == 1
- disp('zero gradient')
- elseif rc == 6
- disp('smallest step still improving too slow, reversed gradient')
- elseif rc == 5
- disp('largest step still improving too fast')
- elseif (rc == 4) | (rc==2)
- disp('back and forth on step length never finished')
- elseif rc == 3
- disp('smallest step still improving too slow')
- end
- end
-
-f=fh;
-x=xh;
-g=gh;
-badg=badgh;
-end
-
-% what about making an m-file of 10 lines including numgrad.m
-% since it appears three times in csminwel.m
-
\ No newline at end of file
diff --git a/MatlabFiles/csminwel.m2 b/MatlabFiles/csminwel.m2
deleted file mode 100755
index 0397abb..0000000
--- a/MatlabFiles/csminwel.m2
+++ /dev/null
@@ -1,271 +0,0 @@
-function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel(fcn,x0,H0,grad,crit,nit, ...
- P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15,P16,P17,...
- P18,P19,P20)
-% [fhat,xhat,ghat,Hhat,itct,fcount,retcodehat] = csminwel(fcn,x0,H0,grad,crit,nit, ...
-% P1,P2,P3,P4)
-%
-% x0(estimated parameter) SHOULD BE A ROW VECTOR.
-%
-% fcn and grad are strings naming functions. If grad is null, numerical
-% derivatives are used. The P parameters need not be present. They are passed
-% to fcn as additional arguments.
-%
-% Reindented and modified to allow compilation (comment out use of tailstr, eval) by CAS,
-% 7/22/96
-%
-[nx,no]=size(x0);
-nx=max(nx,no);
-Verbose=1;
-NumGrad= ( ~isstr(grad) | length(grad)==0);
-done=0;
-itct=0;
-fcount=0;
-snit=100;
-tailstr = ')';
-stailstr = [];
-% Lines below make the number of Pi's optional. This is inefficient, though, and precludes
-% use of the matlab compiler. Without them, we use feval and the number of Pi's must be
-% changed with the editor for each application. Places where this is required are marked
-% with ARGLIST comments
-for i=nargin-6:-1:1
- tailstr=[ ',P' num2str(i) tailstr];
- stailstr=[' P' num2str(i) stailstr];
-end
-f0 = eval([fcn '(x0' tailstr]);
-%ARGLIST
-%f0 = feval(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
-% disp('first fcn in csminwel.m ----------------') % Jinill on 9/5/95
-if f0 > 1e50, disp('Bad initial parameter.'), return, end
-if NumGrad
- if length(grad)==0
- [g badg] = eval(['numgrad(fcn, x0' tailstr]);
- %ARGLIST
- %[g badg] = numgrad(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
- else
- badg=any(find(grad==0));
- g=grad;
- end
- %numgrad(fcn,x0,P1,P2,P3,P4);
-else
- [g badg] = eval([grad '(x0' tailstr]);
- %ARGLIST
- %[g badg] = feval(grad,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
-end
-x=x0;
-f=f0;
-H=H0;
-cliff=0;
-while ~done
- g1=[]; g2=[]; g3=[];
- %addition fj. 7/6/94 for control
- disp('-----------------')
- disp('-----------------')
- %disp('f and x at the beginning of new iteration')
- disp('f at the beginning of new iteration')
- f
- %if size(x,1)>size(x,2)
- % x'
- %else
- % x
- %end
- %-------------------------
- itct=itct+1;
- disp(' going to start csminit.m to produce f1')
- [f1 x1 fc retcode1] = eval(['csminit(fcn,x,f,g,badg,H' tailstr]);
- %ARGLIST
- %[f1 x1 fc retcode1] = csminit(fcn,x,f,g,badg,H,P1,P2,P3,P4,P5,P6,P7,...
- % P8,P9,P10,P11,P12,P13);
- % itct=itct+1;
- fcount = fcount+fc;
- % erased on 8/4/94
- % if (retcode == 1) | (abs(f1-f) < crit)
- % done=1;
- % end
- % if itct > nit
- % done = 1;
- % retcode = -retcode;
- % end
- if retcode1 ~= 1
- if retcode1==2 | retcode1==4
- wall1=1; badg1=1;
- else
- if NumGrad
- [g1 badg1] = eval(['numgrad(fcn, x1' tailstr]);
- %ARGLIST
- %[g1 badg1] = numgrad(fcn, x1,P1,P2,P3,P4,P5,P6,P7,P8,P9,...
- % P10,P11,P12,P13);
- else
- [g1 badg1] = eval([grad '(x1' tailstr]);
- %ARGLIST
- %[g1 badg1] = feval(grad, x1,P1,P2,P3,P4,P5,P6,P7,P8,P9,...
- % P10,P11,P12,P13);
- end
- wall1=badg1;
- % g1
- badg1;
- eval(['save g1 g1 x1 f1 ' stailstr]);
- %ARGLIST
- %save g1 g1 x1 f1 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
- end
- if wall1 % & (~done) by Jinill
- % Bad gradient or back and forth on step length. Possibly at
- % cliff edge. Try perturbing search direction.
- %
- %fcliff=fh;xcliff=xh;
- Hcliff=H+diag(diag(H).*rand(nx,1));
- disp('Cliff. Perturbing search direction.')
- [f2 x2 fc retcode2] = eval(['csminit(fcn,x,f,g,badg,Hcliff' tailstr]);
- %ARGLIST
- %[f2 x2 fc retcode2] = csminit(fcn,x,f,g,badg,Hcliff,P1,P2,P3,P4,...
- % P5,P6,P7,P8,P9,P10,P11,P12,P13);
- fcount = fcount+fc; % put by Jinill
- if f2 < f
- if retcode2==2 | retcode2==4
- wall2=1; badg2=1;
- else
- if NumGrad
- [g2 badg2] = eval(['numgrad(fcn, x2' tailstr]);
- %ARGLIST
- %[g2 badg2] = numgrad(fcn, x2,P1,P2,P3,P4,P5,P6,P7,P8,...
- % P9,P10,P11,P12,P13);
- else
- [g2 badg2] = eval([grad '(x2' tailstr]);
- %ARGLIST
- %[g2 badg2] = feval(grad,x2,P1,P2,P3,P4,P5,P6,P7,P8,...
- % P9,P10,P11,P12,P13);
- end
- wall2=badg2;
- % g2
- badg2
- eval(['save g2 g2 x2 f2 ' stailstr]);
- %ARGLIST
- %save g2 g2 x2 f2 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
- end
- if wall2
- disp('Cliff again. Try traversing')
- if norm(x2-x1) < 1e-13
- f3=f; x3=x; badg3=1;
- else
- gcliff=((f2-f1)/((norm(x2-x1))^2))*(x2-x1)';
- [f3 x3 fc retcode3] = eval(['csminit(fcn,x,f,gcliff,0,eye(nx)' tailstr]);
- %ARGLIST
- %[f3 x3 fc retcode3] = csminit(fcn,x,f,gcliff,0,eye(nx),P1,P2,P3,...
- % P4,P5,P6,P7,P8,...
- % P9,P10,P11,P12,P13);
- fcount = fcount+fc; % put by Jinill
- if retcode3==2 | retcode3==4
- wall3=1; badg3=1;
- else
- if NumGrad
- [g3 badg3] = eval(['numgrad(fcn, x3' tailstr]);
- %ARGLIST
- %[g3 badg3] = numgrad(fcn, x3,P1,P2,P3,P4,P5,P6,P7,P8,...
- % P9,P10,P11,P12,P13);
- else
- [g3 badg3] = eval([grad '(x3' tailstr]);
- %ARGLIST
- %[g3 badg3] = feval(grad,x3,P1,P2,P3,P4,P5,P6,P7,P8,...
- % P9,P10,P11,P12,P13);
- end
- wall3=badg3;
- % g3
- badg3
- eval(['save g3 g3 x3 f3 ' stailstr]);
- %ARGLIST
- %save g3 g3 x3 f3 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
- end
- end
- else
- f3=f; x3=x; badg3=1;
- end
- else
- f3=f; x3=x; badg3=1;
- end
- else
- % normal iteration, no walls, or else we're finished here.
- f2=f; f3=f; badg2=1; badg3=1;
- end
- end
- %how to pick gh and xh
- if f3<f & badg3==0
- ih=3
- fh=f3;xh=x3;gh=g3;badgh=badg3;retcodeh=retcode3;
- elseif f2<f & badg2==0
- ih=2
- fh=f2;xh=x2;gh=g2;badgh=badg2;retcodeh=retcode2;
- elseif f1<f & badg1==0
- ih=1
- fh=f1;xh=x1;gh=g1;badgh=badg1;retcodeh=retcode1;
- else
- [fh,ih] = min([f1,f2,f3]);
- ih
- %eval(['xh=x' num2str(ih) ';'])
- if size(x1,2)>1
- xi=[x1;x2;x3];
- xh=xi(ih,:);
- else
- xi=[x1 x2 x3];
- xh=xi(:,ih);
- end
- %eval(['gh=g' num2str(ih) ';'])
- %eval(['retcodeh=retcode' num2str(ih) ';'])
- retcodei=[retcode1,retcode2,retcode3];
- retcodeh=retcodei(ih);
- if gh==[]
- if NumGrad
- [gh badgh] = eval(['numgrad(fcn, xh' tailstr]);
- %ARGLIST
- %[gh badgh] = numgrad(fcn,xh,P1,P2,P3,P4,P5,P6,P7,P8,...
- % P9,P10,P11,P12,P13);
- else
- [gh badgh] = eval([grad '(xh' tailstr]);
- %ARGLIST
- %[gh badgh] = feval(grad,xh,P1,P2,P3,P4,P5,P6,P7,P8,...
- % P9,P10,P11,P12,P13);
- end
- end
- badgh=1;
- end
- %end of picking
- %ih
- %fh
- %xh
- %gh
- %badgh
- stuck = (abs(fh-f) < crit);
- if (~badg)&(~badgh)&(~stuck)
- H = bfgsi(H,gh-g,xh-x);
- end
- if Verbose
- disp('----')
- disp(['Improvement on iteration ' int2str(itct) ' = ' num2str(f-fh)])
- end
- if Verbose
- if itct > nit
- disp('iteration count termination')
- done = 1;
- elseif stuck
- disp('improvement < crit termination')
- done = 1;
- end
- rc=retcodeh;
- if rc == 1
- disp('zero gradient')
- elseif rc == 6
- disp('smallest step still improving too slow, reversed gradient')
- elseif rc == 5
- disp('largest step still improving too fast')
- elseif (rc == 4) | (rc==2)
- disp('back and forth on step length never finished')
- elseif rc == 3
- disp('smallest step still improving too slow')
- end
- end
- f=fh;
- x=xh;
- g=gh;
- badg=badgh;
-end
-% what about making an m-file of 10 lines including numgrad.m
-% since it appears three times in csminwel.m
-
\ No newline at end of file
diff --git a/MatlabFiles/gradcd.m1 b/MatlabFiles/gradcd.m1
deleted file mode 100755
index dc07b4a..0000000
--- a/MatlabFiles/gradcd.m1
+++ /dev/null
@@ -1,53 +0,0 @@
-function grdd = gradcd(fcn,x0,grdh)
-% computes numerical gradient of a single-valued function or Jacobian
-% matrix of a vector-valued function using a central difference with
-% function grdd = gradcd(fcn,x0,grdh)
-%
-% fcn: a string naming a vector-valued function (f:n*1 -> k*1).
-% x0: a column vector n*1, at which point the hessian is evaluated.
-% grdh: step size, n*1.
-% grdd: Jacobian matrix, k*n.
-
-stps = eps^(1/3);
-% eps: floating point relative accuracy or machine precision: 2.22e-16
-% stps: step size recommended by Dennis and Schnabel: 6.006e-6
-
-% ** initializations
-n = size(x0,1);
-k = size(feval(fcn,x0),1); % dimension of "fcn"
-
-% ** Computation of stepsize (dh)
-if all(grdh)
- dh = grdh;
-else
- ax0 = abs(x0);
- if all(x0)
- dax0 = x0 ./ ax0;
- else
- dax0 = 1;
- end
- dh = stps * (max([ax0 ones(n,1)]'))' .* dax0;
-end
-
-xdh = x0 + dh;
-dh = xdh - x0; % This increases precision slightly
-%
-argplus = x0(:,ones(n,1));
-dnum = 1:n+1:n^2; % positions of diagonals in vec(argplus).
-argplus(dnum) = xdh; % replace the diagonals of "argplus" by "xdh".
-%
-argminus = x0(:,ones(n,1));
-argminus(dnum) = x0-dh; % replace the diagonals of "argplus" by "xdh".
-
-grdd = zeros(k,n); % preallocate to speed the loop.
-i = 0;
-while i ~= n
- i = i+1;
- grdd(:,i) = feval(fcn,argplus(:,i)) - feval(fcn,argminus(:,i));
-end
-dhm = dh(:,ones(k,1));
-dhm = dhm'; % k*n
-grdd = grdd ./ (2*dhm);
-
-
-
\ No newline at end of file
--
GitLab