function slide=molecule
% MOLECULE
%   This demonstration uses the large-scale functionality in 
%   the Optimization Toolbox to solve a two-dimensional Molecule 
%   Conformation Problem.  This type of problem also arises in
%   other applications including satellite ranging and surveying.

% Javier Pena, Matthew J. Simoneau, Mary Ann Branch
% Copyright (c) 1990-98 by The MathWorks, Inc.  
% $Revision: 1.4 $ $Date: 1998/09/09 02:13:25 $

if nargout<1,
  playshow molecule
else
  %========== Slide 1 ==========

slide(1).code={
% The defaultaxesposition needs to be set so subplots behave   
   'if ~isappdata(gcf,''defaultaxesset'');',
   '   set(gcf,''defaultaxesposition'',get(gca,''position''));',
   '   setappdata(gcf,''defaultaxesset'',1);',
   'end;'
   % Draw three spheres connected by cylinders (just for fun)
   'subplot(1,1,1); delete(gca);'
   'hold on;'
   '[cx,cy,cz]=cylinder(.1);'
   'surf(cx,cy,cz+1,''facecolor'',''b'',''linestyle'',''none'');'
   'surf(cz+1,cx,cy,''facecolor'',''b'',''linestyle'',''none'');'
   '[x,y,z]=sphere(40);'
   'surf(x+3,y,z,''facecolor'',''r'',''linestyle'',''none'')'
   'surf(x,y,z+3,''facecolor'',''r'',''linestyle'',''none'')'
   'surf(x,y,z,''facecolor'',''r'',''linestyle'',''none'');'
   'axis equal; axis off;'
   'set(gca,''cameraposition'', [14 24 -23]);'
   'set(gca,''cameraviewangle'', 8);'         
   'set(gca,''cameraupvector'',  [-0.3 -0.6 -.7]);'
   'light;'
   'hold off' };
  
slide(1).text={
     ['This demonstration uses the large-scale functionality in ', ...
           'the Optimization Toolbox to solve a two-dimensional Molecule ',... 
     'Conformation Problem.  This type of problem also arises in ',...
     'other applications including satellite ranging and surveying.']};

%========== Slide 2 ==========

slide(2).code={
   'load molecule;'
   '[smallI,smallJ,smalldist] = find(smallS);'
   
   % Plot the measured distances on the left
   'subplot(1,2,1);'
   'for i = 1:length(smalldist)'
   '   line([0,smalldist(i)],[i i]);'
   '   text(-0.07,i,sprintf(''%d'',smallJ(i)));'
   '   text(smalldist(i)+.02,i,sprintf(''%d'',smallI(i)));'
   'end;'
   ''
   'axis([-.1 1 0 10]); axis off;'
   'title(''Measured distances between Pairs'');'
   
   % Plot the configuration on the right
   'subplot(1,2,2);'
   'smallpts = plot(smallsoln(:,1),smallsoln(:,2),''g.'');'
   'set(smallpts,''markersize'',15);'
   'for i = 1:length(smalldist),'
   ['   line([smallsoln(smallI(i),1) smallsoln(smallJ(i),1)],'...
         '[smallsoln(smallI(i),2) smallsoln(smallJ(i),2)]);']
   'end;' 
   'textpos = smallsoln + [-.09 0;-.08 .08;.04 .01;.04 0;.03 -.07;.03 .07];'
   'for i=1:6,'
   '   text(textpos(i,1),textpos(i,2),sprintf(''%d'',i));'
   'end;'
   'axis([-0.1 1 0 1.1]); axis square; axis off;'
   'title(''Atom Locations'');'};

slide(2).text={
   'MOLECULE CONFORMATION PROBLEM:',
   ['Arrange the N atoms of a molecule in a way that the distances ',...
         'between specified pairs of atoms match experimental data. '],
   '',
   ['A simple example is shown above.  The data is on the left and ',... 
         'the solution on the right.'],
   ''};

%========== Slide 3 ==========

slide(3).code={
   % plot the initial configuration
   'subplot(1,1,1);'
   'N = 25;'
   '[I,J,dist] = find(S);'
   'Xstart = [ [.25 .25]; [.25 .75]; [.75 .25]; reshape(xstart,N-3,2) ];'
   'pts = plot(Xstart(:,1),Xstart(:,2),''g.''); set(pts,''Markersize'',18);'
   'for k = 1:length(dist)',
   ['   h = line([Xstart(I(k),1) Xstart(J(k),1)],'...
         '[Xstart(I(k),2) Xstart(J(k),2)]);'],
   'end;'
   'title(''Initial Locations'');'
   'axis square;'
   '' };

slide(3).text={
   ['An optimization approach to this problem is start off with a ',...
         'random configuration (shown above) and move the atoms until the ',...
         'configuration ',...
         'with the smallest error is found.  Specifically, given '...
         'locations for the atoms:'],
   'X = [x(1:N) y(1:N)]'';',
   ['we associate an error (norm(X(i,:)-X(j,:))^2 - S(i,j)^2)^2 to each ',...
         'pair (i,j) in the set of measured data S. The total error, given by ',...
         'MMOLE(X,S), is the sum of the errors of all pairs in S.  We will ',...
         'find a local minimizer for the error function MMOLE.'],
   ''};

  %========== Slide 4 ==========

slide(4).code={
   % Illustrate the invariances
   'subplot(2,2,1);'
   'pts = plot(smallsoln(:,1),smallsoln(:,2),''g.'');'
   'set(pts,''markersize'',15);'
   'for k = 1:length(smalldist)'
   ['   line([smallsoln(smallI(k),1) smallsoln(smallJ(k),1)],'...
      '[smallsoln(smallI(k),2) smallsoln(smallJ(k),2)]);']
   'end;'
   'axis([0 1 0 1]);'
   'axis square; title(''original'');'
  ''
   'subplot(2,2,2);'
   'pts = plot(1+smallsoln(:,1),2+smallsoln(:,2),''g.'');'
   'set(pts,''markersize'',15);'
   'for k = 1:length(smalldist)'
   ['   line([1+smallsoln(smallI(k),1) 1+smallsoln(smallJ(k),1)],'...
         '[2+smallsoln(smallI(k),2) 2+smallsoln(smallJ(k),2)]);']
   'end;'
   'axis([1 2 2 3]);'
   'axis square; title(''translated'');'
   ''
   'subplot(2,2,3);'
   'pts = plot(smallsoln(:,2),-smallsoln(:,1),''g.'');' 
   'set(pts,''markersize'',15);'
   'for k = 1:length(smalldist)'
   ['   line([smallsoln(smallI(k),2) smallsoln(smallJ(k),2)],'...
         '[-smallsoln(smallI(k),1) -smallsoln(smallJ(k),1)]);']
   'end;'
   'axis([0 1 -1 0]);'
   'axis square; title(''rotated'');'
   ''
   'subplot(2,2,4);'
   'pts = plot(1-smallsoln(:,1),smallsoln(:,2),''g.'');'
   'set(pts,''markersize'',15);'
   'for k = 1:length(smalldist)'
   ['   line([1-smallsoln(smallI(k),1) 1-smallsoln(smallJ(k),1)],'...
      '[smallsoln(smallI(k),2) smallsoln(smallJ(k),2)]);']
   'end;'
   'axis([0 1 0 1]);'
   'axis square; title(''reflected'');'
   'set(findobj(gcf,''type'',''axes''),''Fontsize'',6);'
   '' };

slide(4).text={
   ['Distances are invariant under rigid motion; translations, rotations '... 
         'and reflections do not effect the error function MMOLE.',]
   '',
   ['In order to avoid such motions, we assume that the locations of '...
         'three of the atoms are known. Since these locations are fixed, '...
         'we will actually work with only 2N-6 unknowns.'],
   ''};

  %========== Slide 5 ==========

slide(5).code={
   '',
   '[I,J,dist] = find(S);',
   'cla reset;',
   'subplot(1,1,1);',
   'spy(S+S'');',
   % restoring the renderermode works around a quirk in spy.(g48349)
   'set(gcf,''renderermode'',''auto'');',
   'title(''Measured distances'');' };
slide(5).text={
   'Let''s load an example with 25 atoms:',
   '>> load moledata;',
   '',
   ['The sparse matrix S is a 25 x 25 table.  Each nonzero entry in S '...
         'corresponds to a known distance between a pair of atoms.'],
   '>> spy(S+S'')',
   '',
   ''};

%========== Slide 6 ==========

slide(6).code={
   '[val,g,H] = mmole(xstart,S);',
   'spy(H);',
   % restoring the renderermode works around a quirk in spy.(g48349)
   'set(gcf,''renderermode'',''auto'');',
   'title(''Structure of Hessian Matrix'');',
   '',
   '',
   '' };

slide(6).text={
   ['MMOLE takes the current positions of the atoms and the sparse '...
         'distance matrix (the goal), and returns the distance '...
         'function''s current value (the error), its gradient, and '...
         'its sparse Hessian.'],
   '>> [val,g,H] = mmole(x,S);',
   '>> spy(H);',
   ['Because H is sparse we can use a large-scale optimization algorithm '...
         'to solve this optimization problem.'],
   ''};

%========== Slide 7 ==========

slide(7).code={
   'subplot(1,1,1);',
   'xstart = [ [.25 .25]; [.25 .75]; [.75 .25]; reshape(xstart,N-3,2)];',
   'pts = plot(xstart(:,1),xstart(:,2),''g.''); set(pts,''Markersize'',15);',
   'for k = 1:length(dist)',
   '   energy(k) = (dist(k)^2 - norm(xstart(I(k),:) - xstart(J(k),:))^2)^2;',
   'end;',
   '',
   % Set our own scale of colors'
   '',
   'colormap(energscale);',
   'numcolors = size(energscale,1);',
   'for k = 1:length(dist)',
   ['   h = line([xstart(I(k),1) xstart(J(k),1)],'...
         '[xstart(I(k),2) xstart(J(k),2)]);'],
   ['   set(h,''Color'',energscale(min(numcolors,'...
         'max(1,ceil(numcolors*energy(k)))),:));'],
   'end;',
   'title('' Initial locations and error distribution '');',
   'axis square;',
   'colorbar;',
   
   ['delete(findobj(''type'',''figure'',''Name'', '...
         '''Algorithm Performance Statistics''));'],
   'delete(findobj(''type'',''figure'',''Name'', ''Progress Information''));',
   '' };
slide(7).text={
   ['This is a plot of the atoms in their random starting locations. '...
         'The error between the starting distances and the target '...
         'distances is indicated by the color.'],
   '',
   ['Observe that several edges have a yellow or red color.  This '...
         'indicates there is a big discrepancy between the current '...
         'locations and the experimental data. Our goal is to minimize '...
         'these errors.'],
   ''};

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

slide(8).code={
   'if ~isempty(findobj(gcf,''String'',''AutoPlay'')),'
   '  but(1)=findobj(gcf,''String'',''Next >>'');',
   '  but(2)=findobj(gcf,''String'',''Prev <<'');',
   '  but(3)=findobj(gcf,''String'',''Reset'');',
   '  but(4)=findobj(gcf,''String'',''AutoPlay'');',
   'else'
   '  but=[];'
   'end;'
   'state=get(but,''enable'');',
   'set(but,''enable'',''off'');', 
   'xstart = reshape(xstart(4:N,:),2*N-6,1);',
   ['options = optimset(''largescale'',''on'',''display'',''none'','...
         '''showstatus'',''iterplus'', ''gradobj'','...
         '''on'', ''hessian'',''on'');'],
   'thewindow=gcf;',
   'x = fminunc(''mmole'',xstart,options,S);',
   'figure(thewindow)',
   'for pp=1:length(but),set(but(pp),''enable'',state{pp});end;'};
slide(8).text={
   ['First we set the options for the optimization function FMINUNC. '...
         'We tell it to use the gradient and Hessian calculated '...
         'in mmole.m and show us a final statistics report and '...
         'iteration progress information:'],
   '>> options= optimset(''largescale'',''on'',''display'',',
   ['       ''none'',''showstatus'',''iterplus'',''gradobj'',''on'','... 
         '''hessian'',''on'');'],
   'Now we start the optimization routine FMINUNC',
   '>> x = fminunc(''mmole'',xstart,options,S);',
   'FMINUNC will take about 40 iterations to solve this problem.',
   ''};

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

slide(9).code={
   'subplot(1,1,1)',
   'xsol = [ [.25 .25]; [.25 .75]; [.75 .25]; reshape(x,N-3,2)];',
   'pts = plot(xsol(:,1),xsol(:,2),''g.''); set(pts,''markersize'',15);',
   'for k = 1:length(dist)',
   '   energy(k) = (dist(k)^2 - norm(xsol(I(k),:) - xsol(J(k),:))^2)^2;',
   'end;',
   'maxenergy = max(energy);'
   
   % Set our own scale of colors'
   'colormap(energscale);',
   'numcolors = size(energscale,1);',
   'for k = 1:length(dist)',
   'h = line([xsol(I(k),1) xsol(J(k),1)],[xsol(I(k),2) xsol(J(k),2)]);',
   ['set(h,''Color'','...
         'energscale(min(numcolors,max(1,ceil(numcolors*energy(k)))),:));'],
   'end;',
   'title('' Final error distribution for pairs in S '');',
   'axis([min(xsol(:,1)),max(xsol(:,1)),min(xsol(:,2)),max(xsol(:,2))]);',
   'axis square;',
   'colorbar;',
   
   % Delete LSOT windows
   ['delete(findobj(''type'',''figure'',''Name'', '...
         '''Algorithm Performance Statistics''));'],
   'delete(findobj(''type'',''figure'',''Name'', ''Progress Information''));',
   '',
   '' };
slide(9).text={
   ['This plot shows the atoms in their final positions.  Observe that '...
         'all the line colors are gray.  This means they are all closely '...
         'matching the target distances given by S.'],
   ''};


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

slide(10).code={
   'subplot(1,2,1);',
   'xstart = [ [.25 .25]; [.25 .75]; [.75 .25]; reshape(xstart,N-3,2)];',
   'pts = plot(xstart(:,1),xstart(:,2),''g.''); set(pts,''Markersize'',15);',
   'for k = 1:length(dist)',
   '  energy(k) = (dist(k)^2 - norm(xstart(I(k),:) - xstart(J(k),:))^2)^2;', 
   'end;'
   % Set our own scale of colors
   'colormap(energscale); ',
   'numcolors = size(energscale,1);',
   'for k = 1:length(dist)',
   'h = line([xstart(I(k),1) xstart(J(k),1)],[xstart(I(k),2) xstart(J(k),2)]);',
   ['set(h,''Color'','...
         'energscale(min(numcolors,max(1,ceil(numcolors*energy(k)))),:));'],
   'end;',
   'title('' Initial '');',
   'axis square;',
   'colorbar;',
   'subplot(1,2,2);',
   'pts = plot(xsol(:,1),xsol(:,2),''g.''); set(pts,''markersize'',15);',
   'for k = 1:length(dist)',
   '  energy(k) = (dist(k)^2 - norm(xsol(I(k),:) - xsol(J(k),:))^2)^2; ',
   'end;',
   'maxenergy = max(energy); ',
   % Set our own scale of colors
   'colormap(energscale);',
   'numcolors = size(energscale,1);',
   'for k = 1:length(dist)',
   '  h = line([xsol(I(k),1) xsol(J(k),1)],[xsol(I(k),2) xsol(J(k),2)]);',
   ['  set(h,''Color'','...
         'energscale(min(numcolors,max(1,ceil(numcolors*energy(k)))),:));'],
   'end;',
   'title('' Final '');',
   'axis([min(xsol(:,1)),max(xsol(:,1)),min(xsol(:,2)),max(xsol(:,2))])',
   'axis square;',
   'colorbar;',
};

slide(10).text={
   'Let''s compare the initial and final configuration.'
};


end
