function eval_1d_app_model_obj_fn
% EVAL_1D_APP_MODEL_OBJ_FN: Evaluate 1-D appearance model objective function. In
%                           practice it is a procedure that carries out many
%                           experiments that evaluate the objective function.
%
% Code written by Katherine Smith, 2003
%
%    GENERAL
%
%      eval_1d_app_model_obj_fn
%
%    INPUT/S
%
%      -
%           
%    OUTPUT/S
%
%      -
%
%    PENDING WORK
%
%      -
%
%    KNOWN BUG/S
%
%      -
%
%    COMMENT/S
%
%      
%
%    RELATED FUNCTION/S
%
%      There are some near duplicates in the same directory.
%
%    ABOUT
%
%      -Created:     November 23rd, 2003
%      -Last update: Novermber 30th, 2003
%      -Revision:    0.0.1
%      -Author:      R. S. Schestowitz, University of Manchester
% ==============================================================

       % Smith: make images and points
n_images = 100;
       % Number of images to use for evaluation
n_sets = 5;
       % Number of sets of images (for averaging results); deafult was 50
image_width = 50;
       % Width of the images
bump_width = 0.2;
       % Width of the bump in the images
ref_shift = 0.2;
max_shift = 0.2;
step = 0.1;
       % Warping parameters??
los = zeros(n_images,floor((max_shift-ref_shift)/step));
his = zeros(n_images,floor((max_shift-ref_shift)/step));
       % Will hold the number of positions where bump is high and low
min_error = 0;
max_error = 1;
error_step = 0.1;
       % See section on <Warping images badly, 
       % then modelling poorly warped images> where these are used

n_error_steps = (max_error-min_error)/error_step;
range_score = zeros(n_sets, n_error_steps, 11);
factor_step = 0.1;
factor_min = 0;
factor_max = 0.4;
blurred = 0;
       % blurred images?? (Boolean)
spec_iters = 10;
       % iterations for specificity
gen_iters = 10;
       % iterations for generalisability
build_unwarped_model = 0;
build_perfect_warp_model = 0;
find_gen_and_spec = 1;
       % Booleans that are self-explanatory. They contol of the flow of this
       % routine as it performs many unrelated tests

        
       % The following is similar to build_1d_model. REORGANISE?
for this_set = 1:n_sets
  tic;
  ['calculating set ', num2str(this_set)]
  % for each set of n_sets images, build a model
  % calculate obj fn values from these models, plot the mean with error
  % bars given by 1 stddev away
  [imagelist images points his_no_norm los_no_norm] = make_1d_images(n_images, image_width, bump_width);
  % smooth images
  window = 9;
%    images = average_smooth(images, window); blurred = 1;
%    images = gaussian_smooth(images, window); blurred = 1;

  ref_hi_no_norm = his_no_norm(1);
                                      % Set reference to be image 1??
  ref_lo_no_norm = los_no_norm(1);
  ref_points_vec_no_norm = points(:,1);

  % Smith: normalise points from -1 to 1
  points = -1 + 2*(points-1)/(image_width-1);
  los = -1 + 2*(los_no_norm-1)/(image_width-1);
  his = -1 + 2*(his_no_norm-1)/(image_width-1);
         % Is this some normalisation that is required for reasonable warps


  keep_pts = 0.999999;
  keep_int = 0.8;
  ref_points_vec = points(:,1);
  ref_images_vec = images(:,1);
  ref_hi = his(1);
  ref_lo = los(1);
          % again use image 1 as reference?? Or is it just initialisation?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Smith: Modelling unwarped images

  [intensity_model.pcs,intensity_model.variances,intensity_model.params, intensity_model.mean, intensity_model.var_r, intensity_model.total_var, intensity_model.impdata] = st_pca(images', keep_int);
  [shape_model.pcs,shape_model.variances,shape_model.params, shape_model.mean, shape_model.var_r, shape_model.total_var, shape_model.impdata] = st_pca(points', keep_pts);
          % perform PCA on the images
  intensity_model.sd = std(intensity_model.params);
  shape_model.sd = std(shape_model.params);
          % get standard deviation from the models (of shape and intensity)
  intensity_total_vars = intensity_model.total_var;
  shape_total_vars = shape_model.total_var;
          % get total variance for each model

  if (build_unwarped_model),
    % Smith: combine params - initially normalise to variance described in model
    total_shape_variance = sum(shape_model.variances);
    total_intensity_variance = sum(intensity_model.variances);
          % get total variances
    combined_model.shape_weight = ones(1,size(shape_model.params,2));
          % set weight Ws to initially be 1??
    [combined_model.pcs, combined_model.variances, combined_model.params, combined_model.mean, combined_model.total_var, combined_model.impdata] = make_combined_model(shape_model, intensity_model, combined_model.shape_weight);
          % get combined model
    combined_model.sd = std(combined_model.params);
    combined_model.intensity_model = intensity_model;
    combined_model.shape_model = shape_model;
    combined_model.n_shape_modes = size(shape_model.pcs,2);
    combined_model.label = 'No warp';
             % set some parameters of the model
    [combined_model.mean_specificity, combined_model.std_specificity] = find_specificity(combined_model, images, spec_iters, ref_points_vec);
             % find the specificity of the combined model
    fig_title = 'No warp: ';

%  show_intensity_model(intensity_model, 2, fig_title)

      %This isn't very interesting so don't show it!
%  show_shape_model(shape_model, ref_points_vec, ref_images_vec, 2, white_width)

    show_combined_model(combined_model, ref_points_vec, 2, 5, fig_title)
    disp(['Score of unwarped model: ', num2str(measure_model(combined_model.variances, 50))]);
  end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Smith: Warping images, then modelling warped images
% now warp them using the known correspondences and build another model.

  for i=1:n_images
    % now bung in the right deformations to get them back
    warped_points(:,i) = -1 + 2*(linear_warp(ref_points_vec_no_norm, los_no_norm(i), his_no_norm(i), ref_lo_no_norm, ref_hi_no_norm)-1)/(image_width-1);
            % warping images? Why -1 + 2*warp()???
    warped_points_vec = warped_points(:,i);
            % make a simple copy (duplicate)
    if any(warped_points_vec > 1)
      disp('stop!');
        % RSS: Why should warp have this effect anyway??
    end
    image_vec = images(:,i);
    % resample at the warped points - point on reg. grid has value of point at same index in warpy grid
    % need to interpolate
    warped_images(:,i) = interp1(ref_points_vec, image_vec, warped_points_vec, 'linear',0);
         % apply the warp to the images.
  end
  
  if(build_perfect_warp_model)
	  % This flaf appears at the initialisation section of this function
	  % This is similar to the section above? --> MERGER??
    [w_intensity_model.pcs,w_intensity_model.variances,w_intensity_model.params, w_intensity_model.mean, w_intensity_model.var_r, w_intensity_model.total_var, w_intensity_model.impdata] = st_pca(warped_images', keep_int);
    [w_shape_model.pcs,w_shape_model.variances,w_shape_model.params, w_shape_model.mean, w_shape_model.var_r, w_shape_model.total_var, w_shape_model.impdata] = st_pca(warped_points', keep_pts);
    w_intensity_model.sd = std(w_intensity_model.params);
    w_shape_model.sd = std(w_shape_model.params);
         % As above
    
    % Smith: combine params - initially normalise to variance described in model
    total_shape_variance = sum(w_shape_model.variances);
    total_intensity_variance = sum(w_intensity_model.variances);
    w_c_model.shape_weight = ones(1,size(w_shape_model.pcs,2));
    if(total_shape_variance ~= 0)
      w_c_model.shape_weight(:) = sqrt(total_intensity_variance/total_shape_variance);
    end
           % RSS: What is this about??
    [w_c_model.pcs, w_c_model.variances, w_c_model.params, w_c_model.mean, w_c_model.total_var, w_c_model.impdata] = make_combined_model(w_shape_model, w_intensity_model, w_c_model.shape_weight);
      % combined model
    w_c_model.sd = sqrt(w_c_model.variances);
    w_c_model.intensity_model = w_intensity_model;
    w_c_model.shape_model = w_shape_model;
    w_c_model.n_shape_modes = size(w_shape_model.pcs,2);
    w_c_model.label = 'Perfect warp';
%    [w_c_model.mean_specificity, w_c_model.std_specificity] = find_specificity(w_c_model, images, spec_iters, ref_points_vec);

    fig_title = 'Aligned: ';
    show_combined_model(w_c_model, ref_points_vec, 2, 5, fig_title)
    disp(['Score of perfectly warped model: ', num2str(measure_model(w_c_model.variances, 50))]);
  end
           % RSS: What is it that makes it 'perfect'?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Warping images badly, then modelling poorly warped images

% Again this repeats some of the above in an inefficient way

  error_ctr = 1;
  for error = min_error:error_step:max_error
    ['calculating error ', num2str(error)]
%    for error = -1*(image_width*0.2):image_width*0.2
    for im=1:n_images
      errors(error_ctr) = error;
      % Smith: now bung in deformations to get them back(ish)
      badly_warped_points(:,im) = perturb_points_cps(warped_points(:,im), error);
      badly_warped_points_vec = badly_warped_points(:,im);
      image_vec = images(:,im);
      % Smith: resample at the warped points - point on reg. grid has value of point at same index in warpy grid
      % need to interpolate
      badly_warped_images(:,im) = interp1(ref_points_vec, image_vec, badly_warped_points_vec, 'linear', 0);
      % RSS: apply these bad warps by interpolating.
    end
    [b_w_intensity_model.pcs,b_w_intensity_model.variances,b_w_intensity_model.params, b_w_intensity_model.mean, b_w_intensity_model.var_r, b_w_intensity_model.total_var, b_w_intensity_model.impdata] = st_pca(badly_warped_images', keep_int);
    [b_w_shape_model.pcs,b_w_shape_model.variances,b_w_shape_model.params, b_w_shape_model.mean, b_w_shape_model.var_r, b_w_shape_model.total_var, b_w_shape_model.impdata] = st_pca(badly_warped_points', keep_pts);
    b_w_intensity_model.sd = std(b_w_intensity_model.params);
    b_w_shape_model.sd = std(b_w_shape_model.params);
           % Create the models and get the standard deviation of these models
    % Smith: try using some different normalisation factors
    factor_ctr = 1;
          % This <factor> appears to be the coefficient Ws
    for i=5:10
      factor = 10^i * 1/(10^10); 
          % Silly numbers... 'i' could range from 1 to 5 with same effect 
      models(this_set, error_ctr, factor_ctr).model = build_model_from_models(b_w_shape_model, b_w_intensity_model, ['constant: ', num2str(factor)],'constant', factor);
        % create a model from the models available
      range_score(this_set, error_ctr, factor_ctr) = measure_model(models(this_set, error_ctr, factor_ctr).model.variances, 50);
        % and score it
      if(find_gen_and_spec)
        [models(this_set, error_ctr, factor_ctr).mean_specificity, models(this_set, error_ctr, factor_ctr).std_specificity] = find_specificity(models(this_set, error_ctr, factor_ctr).model, images, spec_iters, ref_points_vec);
        [models(this_set, error_ctr, factor_ctr).mean_generalisability, models(this_set, error_ctr, factor_ctr).std_generalisability] = find_generalisability(models(this_set, error_ctr, factor_ctr).model, gen_iters);
      end
      factor_ctr = factor_ctr+1;
    end

    % Smith: variance
    % RSS: Isn't the following the same as above?
    
    total_shape_variance = sum(b_w_shape_model.variances);
    total_intensity_variance = sum(b_w_intensity_model.variances);
    if(total_shape_variance ~= 0)
      factor = sqrt(total_intensity_variance/total_shape_variance);
    end
    models(this_set, error_ctr, factor_ctr).model = build_model_from_models(b_w_shape_model, b_w_intensity_model, ['variance: ', num2str(factor)],'constant', factor);
    if(find_gen_and_spec)
      [models(this_set, error_ctr, factor_ctr).mean_specificity, models(this_set, error_ctr, factor_ctr).std_specificity] = find_specificity(models(this_set, error_ctr, factor_ctr).model, images, spec_iters, ref_points_vec);
      [models(this_set, error_ctr, factor_ctr).mean_generalisability, models(this_set, error_ctr, factor_ctr).std_generalisability] = find_generalisability(models(this_set, error_ctr, factor_ctr).model, gen_iters);
    end
    range_score(this_set, error_ctr, factor_ctr) = measure_model(models(this_set, error_ctr, factor_ctr).model.variances, 50);
    factor_ctr = factor_ctr+1;

    % Smith: edges
    % RSS: Not sure what this is yet...
    edges = abs(diff(badly_warped_images,1,2));
    factor = sum(edges(:))/size(edges(:),1);
    models(this_set, error_ctr, factor_ctr).model = build_model_from_models(b_w_shape_model, b_w_intensity_model, ['edges: ', num2str(factor)],'constant', factor);
    range_score(this_set, error_ctr, factor_ctr) = measure_model(models(this_set, error_ctr, factor_ctr).model.variances, 50);
    if(find_gen_and_spec)
      [models(this_set, error_ctr, factor_ctr).mean_specificity, models(this_set, error_ctr, factor_ctr).std_specificity] = find_specificity(models(this_set, error_ctr, factor_ctr).model, images, spec_iters, ref_points_vec);
      [models(this_set, error_ctr, factor_ctr).mean_generalisability, models(this_set, error_ctr, factor_ctr).std_generalisability] = find_generalisability(models(this_set, error_ctr, factor_ctr).model, gen_iters);
    end
    factor_ctr = factor_ctr+1;
    
    shape_modes_ctr(error_ctr,this_set) = size(b_w_shape_model.pcs,2);
    intensity_modes_ctr(error_ctr,this_set) = size(b_w_intensity_model.pcs,2);
    error_ctr = error_ctr+1;
  end % error
  t = toc;
  disp(['Time for set ',num2str(this_set),': ',num2str(t)]);
end % set

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot graphs showing mean generalisability and specificity 
% against degree of misregistration

if(find_gen_and_spec)
   for j=1:error_ctr-1
      error_labels(j).label = num2str(errors(j));
      for k=1:factor_ctr-1
        spec_labels(k).label = models(1,1,k).model.label;
        for stupid = 1:n_sets
          temp_spec_means(stupid,k) = models(stupid,j,k).mean_specificity;
          temp_spec_stds(stupid,k) = models(stupid,j,k).std_specificity;
          temp_gen_means(stupid,k) = models(stupid,j,k).mean_generalisability;
          temp_gen_stds(stupid,k) = models(stupid,j,k).std_generalisability;
        end
        spec_means(k,j) = mean(mean(temp_spec_means));
        spec_stds(k,j) = mean(mean(temp_spec_stds));
        gen_means(k,j) = mean(mean(temp_gen_means));
        gen_stds(k,j) = mean(mean(temp_gen_stds));
      end
      error_spec_means(j,:) = mean(temp_spec_means);
      error_gen_means(j,:) = mean(temp_gen_means);
   end
   spec_v_gen_fig = figure;, figure(spec_v_gen_fig), plot(error_spec_means', error_gen_means'), title('Specificity vs generalisability'),...
   xlabel('Specificity'),ylabel('Generalisability');
   cmap = colormap(jet(factor_ctr-1));
   for i=1:factor_ctr-1
     colour = cmap(i,:);
     hold on, plot(error_spec_means(:,i),error_gen_means(:,i),'Marker','o','MarkerFaceColor',colour,'LineStyle','none');
   end
	
   legend(error_labels(:).label,-1);
   cb = colorbar;
   for i=1:factor_ctr-1
      cb_ticks(i) = i+0.5;
      cb_labels{i} = models(1,1,i).model.label;
   end
   set(cb,'YTick',cb_ticks);
   set(cb,'YTickLabel',cb_labels);
	%spec_v_gen_fig = figure;, figure(spec_v_gen_fig), plot(mean(temp_spec_means), mean(temp_gen_means)), title(['Specificity vs generalisability, error ',num2str(j)]),...
	%  xlabel('Specificity'),ylabel('Generalisability');
	%legend(spec_labels(:).label,-1);
   spec_fig = figure;,figure(spec_fig),errorbar(log(spec_means)', spec_stds'./spec_means'), ...
   title('Specificity'),xlabel('Error'),ylabel('Log Specificity (msd of rand example from training example)');
   legend(spec_labels(:).label,-1);
   gen_fig = figure;,figure(gen_fig),errorbar(log(gen_means)', gen_stds'./gen_means'), ...
   title('Generalisability'),xlabel('Error'),ylabel('Log Generalisability (msd of example constructed from model to example');
   legend(spec_labels(:).label,-1);
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot graphs showing log of determinant of c params
% against degree of misregistration

if(blurred == 1)
  fig = figure, cmap = jet(256);, title('Blurred');
else
  fig = figure, cmap = jet(256);, title('Unblurred');
end
line_ctr = 1;
mode_fig = figure;,figure(mode_fig),errorbar([mean(shape_modes_ctr,2) mean(intensity_modes_ctr,2)], [std(shape_modes_ctr,1,2) std(intensity_modes_ctr,1,2)]),title('number of modes'),xlabel('Max error'),ylabel('number of modes');
legend('Shape', 'Intensity');

%for(factor = 0:factor_step:factor_max) 
mean_score = zeros(size(models,2), size(models,3));
for(i = 1:size(models,3)) 
  mean_score(:,i) = mean(range_score(:,:,i))';
  stddev(:,i) = std(range_score(:,:,i))';
  legend_labels(line_ctr).name = models(1,1,line_ctr).model.label;
  line_ctr = line_ctr+1;
end

figure(fig),hold on,lines(line_ctr).line = errorbar(log(mean_score), stddev);
for i=1:size(models,3)
  shape_weights(i) = models(1,1,i).model.shape_weight(1);
end
legend(legend_labels(:).name,-1);
xlabel('Max error');
ylabel('log of determinant of covariance matrix');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot graphs showing generalisability and specificity
% against number of modes retained for each degree of misregistration

if(find_gen_and_spec)
  for i=1:error_ctr-1
      for(j=1:factor_ctr-1)
         for(k=1:n_sets)
           num_modes(j,k) = size(models(k,i,j).model.pcs,2);
         end
      end
      most_modes = max(max(num_modes));
      mean_specs = zeros(factor_ctr-1,most_modes);
      std_specs = zeros(factor_ctr-1,most_modes);
      for j=1:factor_ctr-1
         [a,b] = show_spec_v_modes(models(:,i,j), images, spec_iters, ref_points_vec);
         mean_specs(j,1:size(a,2)) = a;
         std_specs(j,1:size(b,2)) = b;
      %  show_gen_v_modes(models, gen_iters);
      end
      figure,errorbar(mean(mean_specs),mean(std_specs)),xlabel('Number of modes'),ylabel('Specificity'), title(['Error ',num2str(i)]);
      clear mean_specs;
      clear std_specs;
  end
end