www.gusucode.com > mbcexpr 工具箱 matlab 源码程序 > mbcexpr/@cgnormaliser/BP_opt.m
function [LT, cost, OK, msg, varargout] = BP_opt(LT,om, sf) %BP_OPT % Copyright 2000-2013 The MathWorks, Inc. and Ford Global Technologies, Inc. % called by bpopt. % If nargout>3, then varargout contains the spline data that we develop. % % The method proceeds as follows: firstly evaluate the model over the chosen grid and then use this grid and graph % to generate a spline approximation to the graph. To determine whether the current choice of breakpoints is any % good we need to create a lookup table based on the new breakpoints that approximates the model. To do this evaluate % the spline at the new breakpoint positions and use the resulting matrix as the values matrix in the lookup table. % The optimising function then evaluates the lookup table over the chosen grid and subtracts this from the model values % at these values seeking to minimise the difference. mod = get(sf, 'model'); % are all the variables in the table also in the equation? [var, ~, otherVariables]= getvariables(LT,mod); [saveothervar, OK, msg] = setVariables(LT, otherVariables,om); if ~OK if nargout > 4 varargout{1} = []; varargout{2} = []; cost = Inf; end resetVariables(LT, otherVariables,saveothervar); return end % (1) initialise stuff BP = LT.Breakpoints; Xexpr = LT.Xexpr; endpointflag = get(om, 'FixEndPoints'); infoflag = get(om, 'UpdateBPHistory'); % (2) Determine values of x and y that correspond to the chosen grid [savetablevar, OK, msg] = setVariables(LT, var,om); if ~OK if nargout > 4 varargout{1} = []; varargout{2} = []; end cost = Inf; resetVariables(LT, [otherVariables var],[saveothervar savetablevar]); return end abnormalflag = get(om, 'abnormalflag'); if abnormalflag % get variables to use in the table if there was a choice xstr = get(om, 'xVariable'); for i = 1:length(var) if strcmp(var(i).getname, xstr) xvar = var(i); else % if it is not to be used as x, use the set point var(i).info = var(i).set('value', var(i).get('setpoint')); % I don't set this back again.... end end mainvar = xvar; else x = var(1); mainvar = x; end % (3) sets up the comparison matrix M = evaluateaverage(mod.info,mainvar); % (4) Compute approximation spline Xk = Xexpr.eval; % values being fed into normaliser if nargout > 4 varargout{1} = M(:); varargout{2} = Xk; end Vx = LT.Values; % Values array of normaliser Rx = Xexpr.eval; % X inputs to normaliser BPLx = LT.BPLocks; % locked breakpoints if endpointflag == 1 BPLx([1;end]) = 1; % Lock endpoints if needed end Z = BP(BPLx == 0); % values of variables breakpoints Bex = BP(end); % largest breakpoint Bix = BP(1); % minimum BP Zt1 = (BP-Bix)/(Bex-Bix); % Transform BP's to [0,1] indx = (1:length(Zt1))'; Z1L = Zt1(BPLx == 1); % transformed values of locked BP's Idx = indx(BPLx == 1);% index of locked BP's if BPLx(1)==0 % if bottom is not locked, then add it to the locked transformed BP's Z1L = [Zt1(1);Z1L]; Idx = [0;Idx]; end if BPLx(end) == 0 % ditto top Z1L = [Z1L;Zt1(end)]; Idx = [Idx;length(Zt1)]; end if BPLx(1)==0 % if bottom is not locked, then add it to the locked transformed BP's Z1L = [Zt1(end);Z1L]; Idx = [1;Idx]; end if BPLx(end) == 0 % ditto top Z1L = [Z1L;Zt1(end)]; Idx = [Idx;length(Zt1)]; end Zjx = []; for i = 2:length(Z1L) if Idx(i-1)+1<=Idx(i)-1 Zjx = [Zjx;eval(cgmathsobject,'jupp',Zt1(Idx(i-1)+1:Idx(i)-1),Z1L(i-1),Z1L(i))]; % perform Jupp t'form on BP values % between successive locked BP's end end linear1fH = gethandle(cgmathsobject,'linear1'); invjuppfH = gethandle(cgmathsobject,'invjupp'); juppfilterfH = gethandle(cgmathsobject,'juppfilter'); fopts = optimset(optimset('lsqnonlin'),... 'Display','none',... 'PrecondBandWidth',0); if numel(M)<numel(Zjx) % use Levenberg-Marquardt if there is not enough data fopts= optimset(fopts,'LargeScale','off'); %ignore bounds LB=[]; else LB = sqrt(eps)*ones(length(Z),1); end % normalize cost function based on starting values M = M(:); c0 = 1; Zjx(Zjx<=0)= 0.5; z0=i_cost(Zjx,Idx); c0=sqrt(sum(z0.^2)); if ~isfinite(c0) || c0<sqrt(eps) % if the initial cost is too small then we won't try any any scaling c0 = 1; end Zout = lsqnonlin(@i_cost,Zjx,LB,[],fopts,Idx); [~,BP] = undoJupp(Zout,Idx); % clean up if infoflag == 1 LT = set(LT,'breakpoints',{BP(:),'Optimized'}); else LT = setBPunofficial(LT, BP(:)); end resetVariables(LT, [otherVariables var],[saveothervar savetablevar]); cost = Inf; OK = 1; msg = ''; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % cost function function Zout = i_cost(Z,Idx) %adapted from BreakpointAnalysis/Optimisation/BPValopt1 % OUT = BPVALOPT1(Z,Vx,RX,XkS,M,Idx,ct,BPcx,ZLxBex,Bix) % Breakpoint optimisation function for one normaliser to be used with lsqnonlin. Arguments are as follows: % % Z - current Jupped breakpoint values, % Vx - values field of x normaliser, % Rx - X values being used in optimisation, % V - % Xk - % M - cell of comparison matrices, % Idx - Indices of locked x breakpoints, % ZLx - Scaled value of locked x breakpoints, % Bex - value of x endpoint, % Bix - value of x starting point, % % Undo Jupp transform penx = 1; [Znx,BP1] = undoJupp(Z,Idx); % Problem we now face is that although what follows is designed to prevent coalescing BP's, the user % may have forced and locked some BP's to be repeated, in which case the following lines computing penx and peny % will be a source of eternal infinities the like of which will confuse and confound our poor optimiser and lead to % an ever present hour glass marking time until a reboot. So to prevent this undesirable state of affairs and return % things to arrowlike sanity we will attempt to strip out the naughty BP's before computing penx and peny. Why not smack % it with a unique and be done you say, it is left as an exercise to the reader to work out why this frankly stinks - but % as a hint, what if the optimiser puts one BP on top of another, a unique would not allow us to penalise this. if Idx(1) == 0 Idx(1) = []; Idx(end) = []; Znx = feval(juppfilterfH,Znx,Idx); penx = penx*(1-0.1*sum(log((length(Znx)+1)/2*diff([-1;Znx;2])))); else Znx = feval(juppfilterfH,Znx,Idx); penx = penx*(1-0.1*sum(log((Idx(end)+1)/2*diff([-1;Znx;2])))); end n = max(Vx); BPI1 = feval(linear1fH,Vx,BP1,0:n); % find interpolated BP's VAL = feval(linear1fH,Xk,M,BPI1); % get values field corresponding to these new BP's E = feval(linear1fH,BPI1,VAL,Rx); % approximation to model based on these new BP's opt = M-E(:); % error estimate Zout = opt*penx/c0; end function [Znx,BP1] = undoJupp(Z,Idx) % Undo Jupp transform Znx = Z1L(1); count = 0; for j = 2:length(Idx) ztemp = []; if Idx(j-1)+1<=Idx(j)-1 ztemp = feval(invjuppfH,Z(count+1:count+Idx(j)-Idx(j-1)-1),Z1L(j-1),Z1L(j)); count = count+Idx(j)-Idx(j-1)-1; end Znx = [Znx;ztemp;Z1L(j)]; end BP1 = Bix+(Bex-Bix)*(Znx); % undo t'form to [0,1] end end