www.gusucode.com > mbcexpr 工具箱 matlab 源码程序 > mbcexpr/@cgnormfunction/BPinit_curv.m
function [LT,cost,OK, msg] = BPinit_curv(LT,om, sf) %BPINIT_CURV % % LT = BPinit_curv(LT,om, sf) % LT - cgnormfunction % Copyright 2000-2005 The MathWorks, Inc. and Ford Global Technologies, Inc. % see pg 46 of Carl de Boor's 'A Practical guide to splines', % which gives a method of choosing knots for 1D linear splines Norm = LT.Xexpr;% get the normaliser BPx = Norm.get('breakpoints');%get breakpoints BPxval = Norm.get('values'); valx =min(BPxval):max(BPxval); ngBPx = length(valx)-length(intersect(valx, BPxval));% number of green breakpoints in x Xexpr = Norm.get('x');%pointer to x-coordinates of first normalizer eq = get(sf, 'model'); % get pointers to the variables feeding into the LT [var, ~, otherVariables] = getvariables(LT,eq); x = var(1); if length(var) > 1 OK = 0; msg = 'There is more than one variable feeding into the normalizer. Try optimizing instead. '; cost = Inf; return end xval = x.eval;%store the original values in 'x' and 'y' so we can reset at the end saveothervar = setVariables(LT, otherVariables,om); savetablevar = cell(1, length(var)); range = zeros(length(var), 2); for i = 1:length(var) savetablevar{i} = var(i).eval; range(i,1:2) = get(om, ['Range_' var(i).getname]); if range(i,1) >= range(i,2) OK = 0; msg = ['Range_' var(i).getname ' must be strictly increasing.']; cost = Inf; resetVariables(LT, [otherVariables var(1:i)],[saveothervar savetablevar]); return end end BPLx = Norm.get('bplocks');%BPL locked breakpoints nx = length(BPLx);%number of red breakpoints in x % if the locks are empty, assume there are no locks if isempty(nx) nx = length(BPx); BPLx = false(nx,0); end nBPx = nx + ngBPx;%total number of breakpoints % set the values of the variables to be the endpoints of the range x.info = x.set('value',range(1,:)); % get the end breakpoints Xk = Xexpr.eval; % set the end breakpoints if BPLx(1) >0 || BPLx(end) > 0 % if the end points are locked %leave the end breakpoints alone else BPx(1) = min(Xk); BPx(nx) = max(Xk); end lockindx = find(BPLx >0); % lock the end breakpoints lockindx = unique([1; lockindx; nx]); % share the average curvature between the locked breakpoints for i = 1:length(lockindx) -1 % number of breakpoints between two adjacent locked ones nbetweenx = lockindx(i+1) - lockindx(i) - 1; if nbetweenx > 0 % if there is space between locked BPs % Set up grid of test knots, from which we will select a good choice of initial knots % test 10 x (the number of desired knots) possible positions in both x and y numxtestknots = 10*(nbetweenx+2); %10 times the number of bps between the locked ones xtestBP = linspace(BPx(lockindx(i)),BPx(lockindx(i+1)),numxtestknots); % set the xexpr to equal xtestBP, by playing with the variable x [x, OK, msg] = setvalue(Xexpr.info,xtestBP,x); if ~OK OK = 0; msg = 'Problem inverting the normalizer'; resetVariables(LT, [otherVariables var],[saveothervar savetablevar]); x.info = x.set('value',xval); return end % we now take into account the value of eq, at a wider range of values of the other variables M = evaluateaverage(eq.info,x); % check for infs and NaNs anynan = find(isnan(M), 1); anyinfin = find(isinf(M), 1); if ~isempty(anynan) cost = Inf; OK = 0; msg = 'Models returns NaN''s.'; % reset values x.info = x.set('value',xval); resetVariables(LT, otherVariables,saveothervar); return elseif ~isempty(anyinfin) cost = Inf; OK = 0; msg = 'Model gets infinite.'; % reset values x.info = x.set('value',xval); resetVariables(LT, otherVariables,saveothervar); return end xcurvs = diff(M,2); %measure of the total curvature in x between the locked breakpoints totalxcurv= sum(sqrt(abs(xcurvs)))/(numxtestknots-3); %measure of the total curvature in x if totalxcurv == 0 % autospace instead % reset values x.info = x.set('value',xval); resetVariables(LT, otherVariables,saveothervar); [om, OK, msg] =bpinit(LT, BPx(1), BPx(end), valx(end), nBPx); if OK [~, cost, OK, msg] = run(om, LT, []); end if ~OK return end end xcurvsofar = zeros(1, numxtestknots-2); for j=1:numxtestknots-2 xcurvsofar(j) = sum(sqrt(abs(xcurvs(1:j))))/(numxtestknots-3); % a measure of the model curvature from the first breakpoint to the ith breakpoint end Xk = Xexpr.eval; for j=1: nbetweenx %find at what x-value the xcurvsofar exceeds the fraction of the curvature %that should have happened by the jth breakpoint index = find(xcurvsofar > ((BPxval(lockindx(i)+j)-BPxval(lockindx(i)))/(BPxval(lockindx(i+1))-BPxval(lockindx(i))))*totalxcurv, 1); if ~isempty(index) BPx(lockindx(i)+j) = Xk(index); %set this to be the new breakpoint else BPx(lockindx(i)+j) = max(Xk); end end end end BPx = sort(BPx); Norm.info = Norm.set('breakpoints',{BPx(:),'Filled - share ave curvature algorithm.'}); % reset values x.info = x.set('value',xval); resetVariables(LT, otherVariables,saveothervar); cost = Inf; OK =1; msg = '';