% Section 6.5.3, Figure 6.20
% Boyd & Vandenberghe "Convex Optimization"
% Original by Lieven Vandenberghe
% Adapted for CVX by Joelle Skaf - 10/03/05
% (a figure is generated)
%
% Given data u_1,...,u_m and v_1,...,v_m in R, the goal is to fit to the
% data piecewise polynomials with maximum degree 3 (with continuous first
% and second derivatives).
% The [0,1] interval is divided into 3 equal intervals: [-1, -1/3],
% [-1/3,1/3], [1/3,1] with the following polynomials defined on each
% interval respectively:
% p1(t) = x11 + x12*t + x13*t^2 + x14*t^3
% p2(t) = x21 + x22*t + x23*t^2 + x24*t^3
% p3(t) = x31 + x32*t + x33*t^2 + x34*t^3
% L2-norm and Linfty-norm cases are considered

% Input Data
n=4;  % variables per segment
m=40;
randn('state',0);
% generate 50 points ui, vi
u = linspace(-1,1,m);
v = 1./(5+40*u.^2) + 0.1*u.^3 + 0.01*randn(1,m);

a = -1/3;  b = 1/3;  % boundary points
u1 = u(find(u<a)); m1 = length(u1);
u2 = u(find((u >= a) & (u<b)));  m2 = length(u2);
u3 = u(find((u >= b)));  m3 = length(u3);

A1 = vander(u1');   A1 = fliplr(A1(:,m1-n+[1:n]));
A2 = vander(u2');   A2 = fliplr(A2(:,m2-n+[1:n]));
A3 = vander(u3');   A3 = fliplr(A3(:,m3-n+[1:n]));

%L-2 fit
fprintf(1,'Computing splines in the case of L2-norm...');

cvx_begin
    variables x1(n) x2(n) x3(n)
    minimize ( norm( [A1*x1;A2*x2;A3*x3] - v') )
    %continuity conditions at point a
    [1 a a^2   a^3]*x1 == [1 a a^2   a^3]*x2;
    [0 1 2*a 3*a^2]*x1 == [0 1 2*a 3*a^2]*x2;
    [0 0   2 6*a  ]*x1 == [0 0   2 6*a  ]*x2;
    %continuity conditions at point b
    [1 b b^2   b^3]*x2 == [1 b b^2   b^3]*x3;
    [0 1 2*b 3*b^2]*x2 == [0 1 2*b 3*b^2]*x3;
    [0 0   2 6*b  ]*x2 == [0 0   2 6*b  ]*x3;
cvx_end

fprintf(1,'Done! \n');

% L-infty fit
fprintf(1,'Computing splines in the case of Linfty-norm...');

cvx_begin
    variables xl1(n) xl2(n) xl3(n)
    minimize ( norm( [A1*xl1;A2*xl2;A3*xl3] - v', inf) )
    %continuity conditions at point a
    [1 a a^2   a^3]*xl1 == [1 a a^2   a^3]*xl2;
    [0 1 2*a 3*a^2]*xl1 == [0 1 2*a 3*a^2]*xl2;
    [0 0   2 6*a  ]*xl1 == [0 0   2 6*a  ]*xl2;
    %continuity conditions at point b
    [1 b b^2   b^3]*xl2 == [1 b b^2   b^3]*xl3;
    [0 1 2*b 3*b^2]*xl2 == [0 1 2*b 3*b^2]*xl3;
    [0 0   2 6*b  ]*xl2 == [0 0   2 6*b  ]*xl3;
cvx_end

fprintf(1,'Done! \n');

% evaluate the interpolating polynomials using Horner's method
u1s = linspace(-1.0,a,1000)';
p1 = x1(1) + x1(2)*u1s + x1(3)*u1s.^2 + x1(4).*u1s.^3;
p1l1 = xl1(1) + xl1(2)*u1s + xl1(3)*u1s.^2 + xl1(4).*u1s.^3;

u2s = linspace(a,b,1000)';
p2 = x2(1) + x2(2)*u2s + x2(3)*u2s.^2 + x2(4).*u2s.^3;
p2l1 = xl2(1) + xl2(2)*u2s + xl2(3)*u2s.^2 + xl2(4).*u2s.^3;

u3s = linspace(b,1.0,1000)';
p3 = x3(1) + x3(2)*u3s + x3(3)*u3s.^2 + x3(4).*u3s.^3;
p3l1 = xl3(1) + xl3(2)*u3s + xl3(3)*u3s.^2 + xl3(4).*u3s.^3;

us = [u1s;u2s;u3s];
p = [p1;p2;p3];
pl = [p1l1;p2l1;p3l1];
% plot function and cubic splines
d = plot(us,p,'b-',u,v,'go', us,pl,'r--',...
         [-1 -1], [-0.1 0.25], 'k--', [1 1], [-0.1 0.25], 'k--', ...
         [a a], [-0.1 0.25], 'k--', [b b], [-0.1 0.25], 'k--');

title('Approximation using 2 cubic splines');
xlabel('u');
ylabel('f(u)');
legend('L_2 norm','Data points','L_{\infty} norm', 'Location','Best');
% print -deps splineapprox.eps
Computing splines in the case of L2-norm... 
Calling sedumi: 47 variables, 13 equality constraints
   For improved efficiency, sedumi is solving the dual problem.
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
Put 6 free variables in a quadratic cone
eqs m = 13, order n = 5, dim = 49, blocks = 3
nnz(A) = 197 + 0, nnz(ADA) = 169, nnz(L) = 91
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            4.15E+00 0.000
  1 :  -1.74E-01 8.17E-01 0.000 0.1969 0.9000 0.9000   1.71  1  1  1.2E+00
  2 :  -8.96E-02 1.72E-01 0.000 0.2111 0.9000 0.9000   1.25  1  1  3.0E-01
  3 :  -1.16E-01 5.70E-03 0.000 0.0331 0.9900 0.9900   0.91  1  1  4.3E-03
  4 :  -1.17E-01 4.68E-05 0.302 0.0082 0.9990 0.9990   1.00  1  1  4.1E-05
  5 :  -1.17E-01 3.23E-06 0.000 0.0691 0.9900 0.9900   1.00  1  1  4.0E-06
  6 :  -1.17E-01 1.30E-08 0.152 0.0040 0.9990 0.9990   1.00  1  1  1.6E-08
  7 :  -1.17E-01 4.75E-10 0.000 0.0366 0.9900 0.9902   1.00  1  1  1.1E-09

iter seconds digits       c*x               b*y
  7      0.0   8.9 -1.1660335735e-01 -1.1660335750e-01
|Ax-b| =   1.5e-10, [Ay-c]_+ =   2.3E-11, |x|=  2.6e+00, |y|=  1.9e+00

Detailed timing (sec)
   Pre          IPM          Post
0.000E+00    4.000E-02    1.000E-02    
Max-norms: ||b||=1, ||c|| = 2.018994e-01,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 6.24939.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0.116603
Done! 
Computing splines in the case of Linfty-norm... 
Calling sedumi: 126 variables, 53 equality constraints
   For improved efficiency, sedumi is solving the dual problem.
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
Split 6 free variables
eqs m = 53, order n = 133, dim = 133, blocks = 1
nnz(A) = 592 + 0, nnz(ADA) = 577, nnz(L) = 315
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            8.39E+01 0.000
  1 :  -4.68E-01 4.44E+01 0.000 0.5294 0.9000 0.9000   2.90  1  1  3.8E+01
  2 :  -3.12E-02 2.14E+01 0.000 0.4813 0.9000 0.9000  12.63  1  1  1.7E+00
  3 :  -2.65E-02 1.07E+01 0.000 0.4985 0.9000 0.9000   1.83  1  1  7.1E-01
  4 :  -3.14E-02 4.14E+00 0.000 0.3886 0.9000 0.9000   1.17  1  1  2.7E-01
  5 :  -3.16E-02 1.48E+00 0.000 0.3575 0.9000 0.9000   1.10  1  1  9.6E-02
  6 :  -3.20E-02 1.16E-01 0.000 0.0785 0.9900 0.9900   1.02  1  1  7.5E-03
  7 :  -3.20E-02 1.99E-05 0.000 0.0002 0.9999 0.9999   1.00  1  1  
iter seconds digits       c*x               b*y
  7      0.0  15.1 -3.2038334523e-02 -3.2038334523e-02
|Ax-b| =   1.4e-16, [Ay-c]_+ =   1.2E-16, |x|=  7.3e-01, |y|=  2.4e+00

Detailed timing (sec)
   Pre          IPM          Post
0.000E+00    3.000E-02    0.000E+00    
Max-norms: ||b||=1, ||c|| = 2.018994e-01,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 6.03999.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0.0320383
Done!