Home > atmlab > gridcreation > annealing > find_best_freq_set_anneal.m

find_best_freq_set_anneal

PURPOSE ^

FIND_BEST_FREQ_SET_ANNEAL

SYNOPSIS ^

function [sb, Wb, eb, h] = find_best_freq_set_anneal(H, y_mono, n, C)

DESCRIPTION ^

 FIND_BEST_FREQ_SET_ANNEAL

 Find best combination of n frequencies. (The number of
 frequencies to use is fix and given.)

 The idea of this function is from three sources mostly:
 a. http://en.wikipedia.org/wiki/Simulated_annealing. 
 b. Numerical recipes.
 c. Frank Evans description of his pseudo-K approach.

 The algorithm, very roughly, is this:
 1. Pick random selections of frequency
 2. For each frequency set, calculate a weight matrix by
    multilinear regression.
 3. Optimize for the best frequency set by simulated annealing

 The optional input structure C can be used to set some control
 parameters for the simulated annealing algorithm.

 C.Nblock: The calculation is done in blocks. After each block,
           the statistics are evaluated, and the temperature is
           adjusted. This gives the maximum number of moves
           (iterations) in one block. 

 C.Nsucc:  Desired number of successfull moves in one block. If
           Nsucc is reached, the block is terminated, even if
           Nblock is not reached. This makes the algorithm faster
           at high temperatures, where most moves are accepted.

 C.Tfact:  If the error did not decrease in the last block
           (compared to the block before), then the new
           temperature T will be given by T*Tfact.

 C.verb:   If 0, not output.
           If 1, only text output.
           If 2, text and plot output.
 
 C.use_rel_error: Flag, if true then use relative error instead
                  of absolute error.


 FORMAT [sb, Wb, eb, h] = find_best_freq_set_anneal(H,y_mono,n[,C])

 OUT sb = The solution as a logical array.
     Wb = The associated weight matrix.
     eb = The associated error.
     h  = A structure containing the history of the annealing
          run. This is useful for making plots of the convergence
          behaviour.

 IN  H        Sensor response matrix
     y_mono   A batch of monochromatic Tbs (dimension must match
              y_ref and H)
     n        Size of frequency set to test
     C        An optional structure with control parameters.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

find_best_freq_set_anneal.m

SOURCE CODE ^

0001 % FIND_BEST_FREQ_SET_ANNEAL
0002 %
0003 % Find best combination of n frequencies. (The number of
0004 % frequencies to use is fix and given.)
0005 %
0006 % The idea of this function is from three sources mostly:
0007 % a. http://en.wikipedia.org/wiki/Simulated_annealing.
0008 % b. Numerical recipes.
0009 % c. Frank Evans description of his pseudo-K approach.
0010 %
0011 % The algorithm, very roughly, is this:
0012 % 1. Pick random selections of frequency
0013 % 2. For each frequency set, calculate a weight matrix by
0014 %    multilinear regression.
0015 % 3. Optimize for the best frequency set by simulated annealing
0016 %
0017 % The optional input structure C can be used to set some control
0018 % parameters for the simulated annealing algorithm.
0019 %
0020 % C.Nblock: The calculation is done in blocks. After each block,
0021 %           the statistics are evaluated, and the temperature is
0022 %           adjusted. This gives the maximum number of moves
0023 %           (iterations) in one block.
0024 %
0025 % C.Nsucc:  Desired number of successfull moves in one block. If
0026 %           Nsucc is reached, the block is terminated, even if
0027 %           Nblock is not reached. This makes the algorithm faster
0028 %           at high temperatures, where most moves are accepted.
0029 %
0030 % C.Tfact:  If the error did not decrease in the last block
0031 %           (compared to the block before), then the new
0032 %           temperature T will be given by T*Tfact.
0033 %
0034 % C.verb:   If 0, not output.
0035 %           If 1, only text output.
0036 %           If 2, text and plot output.
0037 %
0038 % C.use_rel_error: Flag, if true then use relative error instead
0039 %                  of absolute error.
0040 %
0041 %
0042 % FORMAT [sb, Wb, eb, h] = find_best_freq_set_anneal(H,y_mono,n[,C])
0043 %
0044 % OUT sb = The solution as a logical array.
0045 %     Wb = The associated weight matrix.
0046 %     eb = The associated error.
0047 %     h  = A structure containing the history of the annealing
0048 %          run. This is useful for making plots of the convergence
0049 %          behaviour.
0050 %
0051 % IN  H        Sensor response matrix
0052 %     y_mono   A batch of monochromatic Tbs (dimension must match
0053 %              y_ref and H)
0054 %     n        Size of frequency set to test
0055 %     C        An optional structure with control parameters.
0056 
0057 % 2008-09-25 Created by Stefan Buehler.
0058 
0059 function [sb, Wb, eb, h] = find_best_freq_set_anneal(H, y_mono, n, C)
0060 
0061 
0062 %------------- Set overall control parameters. ----------
0063 
0064 % Default options
0065 def = struct(...
0066     'Nblock', 10000, ...       % Maximum number of cases in one block.
0067     'Nsucc',  100, ...         % Desired number of successfull cases in one block.
0068     'Tfact',  0.9, ...         % Temperature reduction factor.
0069     'verb',   1, ...           % Verbosity.
0070     'use_rel_error', false ... % Use absolute error by default.
0071     );        
0072 
0073 % Check input
0074 if nargin < 4
0075   % Control parameter structure missing, use default.
0076   C = def;
0077 else
0078   % User gave C structure, check which elements are present.
0079   if ~isstruct(C)
0080     error('Input argument ''C'' is not a structure.')
0081   end
0082   fs = {'Nblock', 'Nsucc', 'Tfact', 'verb', 'use_rel_error'};
0083   for nm=1:length(fs)
0084     if ~isfield(C,fs{nm})
0085       C.(fs{nm}) = def.(fs{nm}); 
0086     end
0087   end
0088 end
0089 
0090 % Control structure is now set for sure. Assign to convenient
0091 % internal names.
0092 
0093 % Maximum number of cases in one block:
0094 Nblock = C.Nblock;
0095 
0096 % Alternatively required number of successfull cases in one block:
0097 Nsucc = C.Nsucc;
0098 
0099 % Temperature reduction factor:
0100 Tfact  = C.Tfact;
0101 
0102 % Verbosity:
0103 verb = C.verb;
0104 
0105 % Relative error or absolute error minimization:
0106 use_rel_error = C.use_rel_error;
0107 
0108 %--------------------------------------------------------
0109 
0110 
0111 % Important dimensions:
0112 
0113 x = size(H);
0114 n_channels = x(1);
0115 n_fmono    = x(2);
0116 
0117 % The first dimension of y_mono is the number of monochromatic
0118 % frequencies times the number of viewing angles, since spectra for
0119 % different viewing angles are appended in the output vector y.
0120 
0121 x = size(y_mono);
0122 n_views = x(1)/n_fmono;
0123 n_cases = x(2);
0124 
0125 if verb>0
0126     disp(sprintf('==================================================='));
0127     disp(sprintf('Dimensions of the input fields \n'));
0128     disp(sprintf('==================================================='));
0129     disp(sprintf('H:      n_channels: %d, n_fmono %d\n',n_channels,n_fmono));...
0130     disp(sprintf('y_mono: n_cases: %d, n_y_mono: %d n_views: %d\n',n_cases,x(1),n_views));
0131     disp(sprintf('==================================================='));
0132 end
0133     
0134 
0135 if(mod(x(1),n_fmono) ~= 0)
0136     error('check dimensions/size of ''y_mono''!');
0137 end
0138 % FIXME: check orientation of y_mono
0139 clear x;
0140 
0141 % Carrying the view dimension through the entire calculation would
0142 % make everything very complicated. Especially the calculation of
0143 % regression weights would no longer be straightforward.
0144 % So instead, we choose a different approach here: We treat the
0145 % measurement at different views as if they were different
0146 % atmospheric cases. That means we have to do some re-shaping of
0147 % the input variable y_mono here.
0148 % Original y_mono[n_fmono*n_views, n_cases]
0149 % New      y_mono[n_fmono, n_cases*n_views]
0150 
0151 y_mono_resh = reshape(y_mono,[n_fmono, n_cases*n_views]);
0152 
0153 
0154 % Reference y:
0155 y_ref = H * y_mono_resh;
0156 
0157 % Initial guess.
0158 s = logical(zeros(1,n_fmono));
0159 for i=1:n
0160   % Pick random frequency among those that are not yet active.
0161   r = pick_random_freq(~s);
0162 
0163   % Activate.
0164   s(r) = 1;
0165 end
0166 
0167 
0168 % History inside one block:
0169 e_shorthist = zeros(1,Nblock);
0170 
0171 
0172 % Find out initial temperature. We do this by doing Nblock random
0173 % state changes.
0174 for i=1:Nblock
0175   sn = neighbour(s);  
0176   W  = weights(sn, y_ref, H, y_mono_resh);
0177   e  = test_freq_set(sn, y_ref, W, y_mono_resh, use_rel_error);
0178   e_shorthist(i) = e;
0179   s = sn;
0180 end
0181 % From Vicente et al.
0182 % We want to sete Tstart such that at the beginning almost all
0183 % moves are accepted. (Divide by ln(acceptance probability).)
0184 T = -mean(abs(diff(e_shorthist)))/log(0.99);
0185 
0186 % Save mean energy and other statistical parameters of this initial
0187 % block as start of the history series
0188 
0189 h.t_hist     = NaN(1,1000);
0190 h.e_hist     = NaN(1,1000);
0191 h.e_std_hist = NaN(1,1000);
0192 h.e_min_hist = NaN(1,1000);
0193 h.e_max_hist = NaN(1,1000);
0194 
0195 h.t_hist(1)     = T;
0196 h.e_hist(1)     = mean(e_shorthist);
0197 h.e_std_hist(1) = std(e_shorthist);
0198 h.e_min_hist(1) = min(e_shorthist);
0199 h.e_max_hist(1) = max(e_shorthist);
0200 
0201 if verb>0
0202   disp(sprintf('mean(abs(delta_e)) = %g', mean(abs(diff(e_shorthist)))));
0203   disp(sprintf('T(start) = %g', T));
0204 end
0205 
0206 % Initialize best state energy:
0207 eb = h.e_max_hist(1);
0208 
0209 
0210 k = 1;
0211 
0212 
0213 go_on = true;
0214 while go_on
0215   
0216   % Do a block of steps at a time
0217   n_succ = 0;                           % Count successful steps
0218   for i=1:Nblock
0219 
0220     % Select neighbour.
0221     sn = neighbour(s);  
0222     
0223     % Calculate weights.
0224     W  = weights(sn, y_ref, H, y_mono_resh);
0225 
0226     % Calculate error (= energy).
0227     en = test_freq_set(sn, y_ref, W, y_mono_resh, use_rel_error);
0228     
0229     % Energy difference:
0230     de = en - e;
0231     
0232     % Update best state, if this one is better.
0233     if en < eb
0234       sb = sn;
0235       Wb = W;
0236       eb = en;
0237     end
0238 
0239     % Probability of this state:
0240     Pk = exp(-de/T);
0241     
0242     if rand(1,1) < Pk
0243       s = sn;
0244       e = en;      
0245       
0246       n_succ = n_succ + 1;   
0247     end
0248     
0249     e_shorthist(i) = e;    
0250 
0251     if n_succ >= Nsucc
0252       break
0253     end
0254   end
0255     
0256   k = k+1;
0257   
0258   h.t_hist(k)     = T;
0259   h.e_hist(k)     = mean(e_shorthist(1:i));
0260   h.e_std_hist(k) = std(e_shorthist(1:i));
0261   h.e_min_hist(k) = min(e_shorthist(1:i));
0262   h.e_max_hist(k) = max(e_shorthist(1:i));
0263 
0264   if verb>0
0265     if (use_rel_error)
0266       disp(sprintf('T = %g, e_min = %g, e_mean = %g, e_max = %g [fractional errors]',...
0267                    h.t_hist(k),...    
0268                    h.e_min_hist(k),...    
0269                    h.e_hist(k),...    
0270                    h.e_max_hist(k) ))
0271     else
0272       disp(sprintf('T = %g K, e_min = %g K, e_mean = %g K, e_max = %g K',...
0273                    h.t_hist(k),...    
0274                    h.e_min_hist(k),...    
0275                    h.e_hist(k),...    
0276                    h.e_max_hist(k) ))
0277     end
0278   end
0279     
0280   % Display results
0281   if verb>1
0282     figure(1)
0283     semilogy([h.e_min_hist(1:k); ...
0284               h.e_hist(1:k); ...
0285               h.e_max_hist(1:k)]');
0286     xlabel('Iteration')
0287     if (use_rel_error)
0288       ylabel('RMS error [fractional]')
0289     else
0290       ylabel('RMS error [K]')
0291     end
0292   
0293     figure(2)
0294     semilogy(h.t_hist(1:k));
0295     xlabel('Iteration')
0296     if (use_rel_error)
0297       ylabel('Temperature [fractional]')
0298     else
0299       ylabel('Temperature [K]')
0300     end
0301       
0302     figure(3)
0303     loglog(h.t_hist(1:k), h.e_hist(1:k), '.');
0304     xlabel('Temperature [K]')
0305     if (use_rel_error)
0306       ylabel('RMS error [fractional]')
0307     else
0308       ylabel('RMS error [K]')
0309     end
0310   end
0311 
0312   % Should we decrease T? - Not if the mean error is still
0313   % decreasing.
0314   if h.e_hist(k) >= h.e_hist(k-1)
0315     T = T * Tfact;
0316   end
0317   
0318   % Should we continue? - We use as stop criterion that there were
0319   % no more successful moves
0320   if n_succ == 0
0321     go_on = false;
0322   end
0323 end
0324   
0325 
0326 if verb>0
0327   disp('Best combination:')
0328   if (use_rel_error)
0329     disp(sprintf('RMS error:   %g [fractional]', eb))
0330   else
0331     disp(sprintf('RMS error:   %g [K]', eb))
0332   end
0333 end
0334

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005