function [mean_specificity, std_specificity, all_specificities] = find_specificity(model, training_data, n_iterations, ref_points_vec, varargin)
% FIND_SPECIFICITY: Find the specificity of a model.
%
% Code written by Katherine Smith, 2003
%
%    GENERAL
%
%      [mean_specificity, std_specificity, all_specificities] =
%       find_specificity(model, training_data, n_iterations, 
%       ref_points_vec, varargin)
%
%    INPUT/S
%
%      -model:
%           The model for which specificity is evaluated.
%
%      -training_data:
%           The training data used to evalaute the specificity
%           of that model.
%
%      -n_iterations:
%           ?The number of iterations used in the evaluation process?
%
%      -ref_points_vec:
%           ?Reference points?
%
%      -varargin:
%           ?Variance argument input?
%           
%    OUTPUT/S
%
%      -mean_specificity:
%           The mean specificity.
%
%      -std_specificity:
%           The standard deviation of specificities.
%
%      -all_specificities:
%           All the specificity values found.
%
%    PENDING WORK
%
%      -Resolve meaning of inputs and outputs.
%
%    KNOWN BUG/S
%
%      -None.
%
%    COMMENT/S
%
%      -Smith: Construct examples from model and measure 
%       the distance from the nearest example in training set
%      -Full definition given by Rh. Davies.
%
%    RELATED FUNCTION/S
%
%      FIND_GENERALISABILITY
%
%    ABOUT
%
%      -Created:     November 23rd, 2003
%      -Last update: Novermber 25th, 2003
%      -Revision:    0.0.2
%      -Author:      R. S. Schestowitz, University of Manchester
% ==============================================================

args = u_packargs(varargin, 0,'n_modes', size(model.params,2));

best_yet = ones(n_iterations,1)*size(training_data,1); %construct vector of as many ones as iternations and multiply these ones by the size of the training data
stds = sqrt(model.variances); % get square roots of model variances
stds = real(stds); % convert these to real numbers
for i=1:n_iterations % for input number of iterations
  params = zeros(1,size(model.params,2)); 
         % paramters set to 0, their size being the length of second bunch of model parameters
  invalid_params = 1; % set initially to invalid
  while(invalid_params) % and iterate while no validity proven
    params(1:args.n_modes) = randn(1,args.n_modes).*stds(1:args.n_modes);
      % for all modes, set them to a number in the range of their total number * the variance in the position of number of modes
    [example,points] = construct_model_example(params, model, ref_points_vec);
     % create a few examples given the parameters
    if(~is_invalid(points)) % if all our valid
      invalid_params = 0;  % break out of loop
    end
  end

  for j=1:size(training_data,2)
    distance = msd(example, training_data(:,j)); % get mean squred distance
    if(distance < best_yet(i)) % and try to find the closest yet
      best_yet(i) = distance; 
    end
  end
  if(best_yet(i) > 10)
    disp('stop here');
  end
end
mean_specificity = mean(best_yet);
std_specificity = std(best_yet);
all_specificities = best_yet;
