spline_h.m

Incorporating Prior Knowledge in Cubic Spline Approximation - Application to the Identification of R

function [P1,P2]  = spline_h(cold,tt,cv,reler,ter,knots,cycnum)
%
% cold -> valodi gorbe
% tt ->  idovektor a 12 v. 31 ponthoz
% cc -> 12x3 v. 31x3 matrix adatokhoz
% reler -> rel.zaj 0.05 0.10
% ter -> zaj tipusa: 0=fuggetlen,1=fuggo
% knots -> csomopontok
% cycnum -> ciklusszam
% P1 <- p1 vector
% P2 <- p2 vector
%

seeki = 0; randn('seed',seeki);
P1 = [];
P2 = [];

for II = 1:cycnum,
 pause(0.01); 
 n = length(knots); %number of nodes
 nc = 3; %number of components
 %+noise
 ca = cv(:,1);
 cb = cv(:,2);
 cc = cv(:,3);
 if ter==0,
  ca = ca + randn(size(ca)).*ca*reler;
  cb = cb + randn(size(cb)).*cb*reler;
  cc = cc + randn(size(cc)).*cc*reler;
 else
  rr = randn(size(ca));
  ca = ca + rr.*ca*reler;
  cb = cb + rr.*cb*reler;
  cc = cc + rr.*cc*reler;
 end
 %eq. hard const. /Ca+Cb+Cc = 1, dCa/dt+...+dCc/dt = 0/
 A=repmat(eye(n*2),1,nc);
 b=repmat([1 0],1,n)';
 %/Ca(1) = 1, Cb(1) = 0, Cc(1) = 0/
 A(1,1) = 1;
 A(2*n+1,2*n+1) = 1;
 A(2*n+2,4*n+1) = 1;
 b(1,1) = 1;
 b(2*n+1,1) = 0;
 b(2*n+2,1) = 0;
 %neq. hard const. /dCa/dt < 0, dCb/dt>0 dCc/dt>0 /
 Anv = zeros(n,2*n);
 for i = 1:n
  Anv(i,2*i) = 1;
 end
 An = zeros(n*nc,2*n*nc);
 An(1:n,1:2*n) = Anv;
 An(n+1:2*n,2*n+1:4*n) = -Anv;
 An(2*n+1:3*n,4*n+1:6*n) = -Anv;
 bn = zeros(n*nc,1);
 %no soft. constr.
 s = [];
 lambda1 = 0;
 lambda2 = 0;
 w = [0 0 0];
 sum1 = 0;
 sum2 = 0;
 %Spline
 warning off
 [spa1,spb1,spc1] = fspline3([tt,tt,tt],[ca,cb,cc],knots,s,lambda1,lambda2,w,sum1,sum2,A,b,An,bn);
 warning on
 %P1 and P2
 te = [0:0.02:2]';
 cnew = zeros(length(te),3);
 cnew(:,1) = drawspline(spa1,te);
 cnew(:,2) = drawspline(spb1,te);
 cnew(:,3) = drawspline(spc1,te);
 P2 = [P2, sum(sum((cnew-cold).*(cnew-cold)))];
 t = te;
 c = cnew;
 for j=1:length(te),
  ic(j,1) = intspline(spa1,[0 t(j)]);
  ic(j,2) = intspline(spb1,[0 t(j)]);
  ic(j,3) = intspline(spc1,[0 t(j)]);
 end
 libarycommand = 3;
 rate_libary;
 P1 = [P1, norm((b'-[1 0.5 10 5])./[1 0.5 10 5])];
end