function slide=circustent
% This demonstrates that the shape of a circus tent is 
% determined by a constrained optimization problem
% that we can solve using the large-scale optimization
% functionality in the Optimization Toolbox.

%   Copyright 1990-2000 The MathWorks, Inc.
%   $Revision: 1.7 $  $Date: 2000/06/16 22:23:29 $

if nargout<1,
   playshow circustent
else
   %========== Slide 1 ==========
   
   slide(1).code={
      'watchon;',
      'load tentdata;',
      'set(gcf,''name'',''circustent'');',
      'set(gcf,''doublebuffer'',''on'',''render'',''zbuffer'');',
      'set(gcf,''defaultaxesposition'',get(gca,''position''));',
      'cla;',
      's=surfl(aproxtent);',
      '%Draw sketch of tent',
      'set(s,''edgecolor'',''none'');',
      'shading interp;',
      'surface(zeros(size(aproxtent)));',
      'patch([15 18 15],[15 18 15],[.35 .35 .4],''r'');',
      'view([-20,30]);',
      'colormap(gray);',
      'axis tight;',
      'axis off;',
      'rotate3d on;',
      'drawnow;',
      'watchoff;',
      ''};
   slide(1).text={
      ['The shape of a circus tent is determined by a constrained ',...
            'optimization problem.  We can solve this problem with ',...
            'the large-scale optimization functionality in the ',...
            'Optimization Toolbox'],
      ''};
   
   %========== Slide 2 ==========
   
   slide(2).code={
      'watchon;',
      'delete(gca);',
      'gca;',
      'axis([1 30 1 30 0 .5]);',
      'axis off;',
      'rotate3d on;',
      'largeL = zeros(36);',
      'mask = [6 7 30 31];',
      'largeL(mask,mask) = .3*ones(4);',
      'largeL(18:19,18:19) = .5*ones(2);',
      'xx = [1:5,5:6,6:15,15:16,16:25,25:26,26:30];',
      '[XX,YY] = meshgrid(xx) ;',
      'surface(XX,YY,largeL,''facecolor'',[.5 .5 .5],''edgecolor'',''none'');',
      'light;',
      'colormap(gray);',
      'view([-20,30]);',
      'title(''The set of tent poles'')',
      'drawnow;',
      'watchoff;',
      ''};
   slide(2).text={
      ['Imagine building a circus tent to cover a square lot. ',...
            'The tent has five poles that will be covered with an ',...
            'elastic material. From this structure, we want to find ',...
            'the natural shape of the tent. This natural shape corresponds to ',...
            'the minimum ',...
            'of a certain energy function computed ',...
            'from the surface position and ',...
            'squared norm of its gradient.'],
      ''};
   
   %========== Slide 3 ==========
   
   slide(3).code={
      'watchon;',
      'L = zeros(30); E = ones(2);',
      'L(15:16,15:16) = .5*E;',
      'L(5:6,5:6) = .3*E; L(25:26,5:6) = .3*E;',
      'L(5:6,25:26) = .3*E; L(25:26,25:26) = .3*E;',
      '',
      'surface(L,''facecolor'',''none'',''edgecolor'',''m'');',
      'title(''Lower Bound Constraint Surface'')',
      'drawnow;',
      'watchoff;',
      '' };
   slide(3).text={
      'The supporting poles determine a lower bound for the tent:',
      '>> L = zeros(30);  E = ones(2);',
      '>> L(15:16,15:16) = .5*E;',
      '>> L(5:6,5:6) = .3*E; L(25:26,5:6) = .3*E;',
      '>> L(5:6,25:26) = .3*E; L(25:26,25:26) = .3*E;',
      'We can visualize the constraint by plotting it as a magenta mesh:',
      '>> mesh(L,''edgecolor'',''m'')',
      ''};
   
   %========== Slide 4 ==========
   
   slide(4).code={
      'watchon;',
      'sstart=.5*ones(30,30);',
      'cla; axis([1 30 1 30 0 .5]);view([-20,30]);',
      'surface(L,''facecolor'',''none'',''edgecolor'',''m'');',
      ['surface(sstart,''facecolor'',''none'',''linestyle'','...
            '''none'',''marker'',''.'',''markeredgecolor'',''b'')'],
      'title(''Initial Value (blue) and Lower Bound (magenta)'');',
      'drawnow;',
      'watchoff;',
      '' };
   slide(4).text={
      ['To solve this problem, we will find the height of the ',...
            'optimized surface at a finite number of points. Our ',...
            'initial guess, sstart, is shown in blue.'],
      '>> sstart=.5*ones(30,30);',
      '>> s= surf(sstart)',
      '>> set(s,''FaceColor'',''none'',''LineStyle'',''none'', ...',
      '       ''Marker'', ''.'', ''MarkerEdgeColor'',''b'');',
      ''};
   
   %========== Slide 5 ==========
   
   slide(5).code={
      'watchon;',
      'delete(findobj(gcf,''tag'',''deleteme''));',
      'low = reshape(L,900,1);',
      'xstart=reshape(sstart,900,1);',
      '%Illustrate the reordering',
      '%Draw grid points',
      'xx = 0:4;',
      '[X Y] = meshgrid(xx,xx);',
      'gpts = plot(X(:),Y(:),''b.'');',
      'set(gpts,''markersize'',10);',
      'axis off; axis([-2 12 -1.5 5.5]);',
      'hold on',
      '%Draw arrow',
      'l(1)=line([7.5 6.5],[2 2.5]);',
      'l(2)=line([7.5 6.5],[2 1.5]);',
      'l(3)=line([7.5 5.5],[2 2]);',
      'set(l,''color'',''b'');',
      '%Draw vector',
      'yy = 0.2*xx;',
      'zz = [-1.5+yy,yy,1.5+yy,3+yy,4.5+yy];',
      'vect = plot(9*ones(25,1),zz,''b.'');',
      'set(vect,''markersize'',9);',
      'axis off;',
      'hold off;',
      'drawnow;',
      'watchoff;',
      '' };
   slide(5).text={
      ['In order to formulate the problem as a standard optimization ',...
            'problem, we resize both the matrices into vectors.'],
      'L representing the initial values:',
      '>> low = reshape(L,900,1);',
      'sstart representing the lower boundary constraint',
      '>> xstart=reshape(sstart,900,1);',
      ''};
   
   %========== Slide 6 ==========
   
   slide(6).code={
      'watchon;',
      'str=''min { c''''*x+0.5*x''''*H*x : x > low }'';',
      ['h=uicontrol(''style'',''text'',''string'',str,',...
            '''units'',''normalized'',',...
            '''position'',get(gcf,''defaultaxesposition''),',...
            '''fontsize'',20,',...
            '''backgroundcolor'',get(gcf,''color''),',...
            '''tag'',''deleteme'');'],   
      'H = delsq(numgrid(''S'',30+2));',
      'h = 1/(30-1);',
      'c = -h^2*ones(30^2,1);',
      'delete(gca)',
      'drawnow;',
      'watchoff;',
      '' };
   slide(6).text={
      ['The surface formed by the elastic membrane is determined ',...
            'by the linearly constrained optimization problem'],
      '            min{ c''*x+0.5*x''*H*x : low <= x }',
      ['where c''*x + 0.5*x''*H*x is the discrete approximation ',...
            'of the energy function. H and c are as follows:'],
      '>> H = delsq(numgrid(''S'',30+2));',
      '>> h = 1/30-1; >> c = -h^2*ones(30,1);',
      ''};
   
   %========== Slide 7 ==========
   
   slide(7).code={
      'watchon;',
      'delete(findobj(gcf,''tag'',''deleteme''));',
      'spy(H);',
      'title(''Structure of Hessian Matrix'');',
      'drawnow;',
      'watchoff;',
      '' };
   slide(7).text={
      ['Each point of the energy function is only affected by ',...
            'its immediate neighbors.  Consequently the Hessian ',...
            'matrix H is sparse and has a special structure.'],
      '>> spy(H)',
      ['Because H is sparse we can use a large-scale ',...
         'algorithm to solve the optimization problem.'],
   ''};


%========== Slide 8 ==========

slide(8).code={
   'watchon;',
   '%Draw baby versions of the X-G plots shown in the LSOT progress info',
   'plot(xXX3,XX3,''b.'',xXX1,XX1,''r.'',xXX2,XX2,''r.'');',
   'set(gca,''YTick'',[-1 1]);',
   'set(gca,''YTickLabel'',{''lower'';''upper''});',
   'axis([1 900 -1 1]);',
   'title(''Relative position of x(i) to upper and lower bounds (log-scale)'');',
   'drawnow;',
   'watchoff;',
   '' };
slide(8).text={
   ['It will take about 14 iterations to solve this ',...
         'problem. At each iteration a progress information ',...
         'window will be updated. Before starting the ',...
         'optimization, we will look at the two plots ',...
         'found in the window.'],
      ['An example of the first plot is shown above.  Each component ',...
         'is plotted against its position relative to the lower ',...
         'and upper bounds. If a component is at any of its ',...
         'bounds it is plotted in red.  Otherwise, it is strictly ',...
         'between its bounds and plotted in blue.'],
   ''};

%========== Slide 9 ==========

slide(9).code={
   'watchon;',
   '%Now show the G plot',
   'currplot=plot(xGG2,GG2,''b.'',xGG1,GG1,''r.'');',
   'title(''Relative gradient scaled to the range -1 to 1'');',
   'legend(''abs(grad) > tol'',''abs(grad) <= tol'',4);',
   'axis([1 900 -1 1]);',
   'set(gca,''YTick'',[-1 0 1]);',
   'ylabel(''gradient'')',
   '%Delete LSOT windows (in case we''re going backwards)',
   'delete(findobj(''type'',''figure'',''Name'',''Algorithm Performance Statistics''));',
   'delete(findobj(''type'',''figure'',''Name'',''Progress Information''));',
   'drawnow;',
   'watchoff;',
   '' };
slide(9).text={
   ['Here is an example of the second plot found on the ',... 
         'status window. ',...  
         'The second plot shows the normalized components of the ',...
         'gradient g.  The components that are close to ',...
         'zero (within a certain tolerance) are plotted in ',...
         'red.  The rest are plotted in blue.'],
   ''};

%========== Slide 10 ==========

slide(10).code={
   'watchon;',
   'delete(gca)',
   'options = optimset(''LargeScale'',''on'',''display'',''off'',''ShowStatusWindow'',''iterplus'');',
   'x = quadprog(H, c, [], [], [], [], low, [], xstart, options);',
   %   'load x;',
   %'figure(findobj(''type'',''figure'',''name'',''circustent''));',   
   'drawnow;',
   'watchoff;', 
   ''};
slide(10).text={
   'First we set the options for QUADPROG:',
   '>> options = optimset(''LargeScale'',''on'', ''Display'', ''off'', ...',
   '         ''ShowStatusWindow'', ''iterplus'');',
   '',
   'Now we solve the problem by using the routine QUADPROG:',
   '>> x = quadprog(H, c, [], [], [], [], low, [], xstart, options);',
   ''};

%========== Slide 11 ==========

slide(11).code={
   'watchon;',
   % Delete LSOT windows'
   'delete(findobj(''type'',''figure'',''Name'',''Algorithm Performance Statistics''));',
   'delete(findobj(''type'',''figure'',''Name'',''Progress Information''));',
   % reshape solution
   'S = reshape(x,30,30);',
   % lefthand plot (starting)
   'subplot(1,2,1);',
   'surf(L,''facecolor'',[.5 .5 .5]);',
   'surface(sstart,''edgecolor'',''b'',''facecolor'',''none'');',
   'title(''Starting Surface'')',
   'axis off',
   'axis tight;',
   'view([-20,30]);',
   % righthand plot (solution)
   'subplot(1,2,2);',
   'surf(L,''facecolor'',[.5 .5 .5]);',
   'surface(S,''edgecolor'',''b'',''facecolor'',''none'');',
   'title(''Solution Surface'')',
   'axis off',
   'axis tight;',
   'view([-20,30]);',
   'drawnow;',
   'watchoff;',
   ''};
slide(11).text={
   ['Now obtain the surface solution by going back to the ',...
         'original ordering'],
   '>> S = reshape(x,30,30);',
   'Plot this solution in blue mesh:',
   '>> mesh(S,''edgecolor'',''b'',''facecolor'',''none'');',
   ''};

%========== Slide 12 ==========

slide(12).code={
   'watchon;',
   '%Draw tent',
   'subplot(1,1,1)',
   'surf(L,''facecolor'',[0 0 0]);',
   'hold on;',
   'surfl(S);',
   'hold off;',
   'axis tight;',
   'axis off;',
   'view([-20,30]);',
   'drawnow;',
   'watchoff;',
   ''};
slide(12).text={
   'Let''s see the circus tent:',
   '>> surfl(S)',
   ''};
end