function result = mf(X,model,mode,options);
%MF Matrix factorization for simplifying PCA and Tucker models.
%  Takes a PCA or Tucker model and rotates one mode to optimal simplicity.
%  The input is the data used to fit the original model, the model 
%  structure and the mode to be rotated. The remainder of the model is 
%  counterrotated so that the fit of the model is the same. In PCA 
%  the scores or loadings are counterrotated whereas in Tucker, the core 
%  is counterrotated.
%  The rotation is done by orthomax where default is to set the parameter 
%  gamma to one (corresponds to varimax) and where setting it to
%  zero would correspond to quartimax.
%  Optional structure input (options) can contain any of the fields:
%
%    display   [ {'off'} | 'on']   Governs screen display.
%    gamma     between 0 and 1. One is default.
%
%I/O: result = mf(X,model,mode,options);
%
%See also: PCA, TUCKER

% Copyright © Eigenvector Research, Inc. 1991-2006
% Licensee shall not re-compile, translate or convert "M-files" contained
%  in PLS_Toolbox for use with any software other than MATLAB®, without
%  written permission from Eigenvector Research, Inc.

if nargin == 0; X = 'io'; end
if ischar(X);
  options = [];

  options.plots = 'on';
  options.display = 'on';
  options.preprocessing = {[]};
  options.confidencelimit = 0.95;
  options.gamma = 1;
  
  if nargout==0; evriio(mfilename,X,options); else; result = evriio(mfilename,X,options); end
  return; 
end

if nargin<3
  error(' You must define the mode in which you want to rotate the components in MF')
end

if nargin < 4 | isempty(options);
  options = mf('options');
else
  options = reconopts(options,'mf');
end

if options.gamma>1
  warning('Gamma can at most be one')
  options.gamma = 1;
elseif options.gamma<0
  warning('Gamma must be at least zero')
  options.gamma = 0;
end



if isa(X,'double')      %convert X to DataSet
  X        = dataset(X);
  X.name   = inputname(1);
  X.author = 'MCR';
elseif ~isa(X,'dataset')
  error(['Input X must be class ''double'' or ''dataset''.'])
end
if ndims(X.data)>2
  error(['Input X must contain a 2-way array. Input has ',int2str(ndims(varargin{1}.data)),' modes.'])
end


%Handle Preprocessing
if isempty(options.preprocessing);
  options.preprocessing = {[]};  %reinterpet as empty cell
end
if ~isa(options.preprocessing,'cell');
  options.preprocessing = {options.preprocessing};  %insert into cell
end


if mdcheck(X.data(X.include{1},X.include{2}));
  if strcmp(options.display,'on'); 
    warning('Missing Data Found - Replacing with "best guess". Results may be affected by this action.'); 
  end
  [flag,missmap,replaced] = mdcheck(X.data(X.include{1},X.include{2}));
  X.data(X.include{1},X.include{2}) = replaced;
end

%preprocessing
if ~isempty(options.preprocessing{1});
  [X,options.preprocessing{1}] = preprocess('calibrate',options.preprocessing{1},X);
end


% GO AHEAD
if strcmp(model.modeltype,'PCA')
  [Arot,T,f]=orthomax(model.loads{mode},options.gamma);
  if mode==1
    othermode = 2;
  else
    othermode = 1;
  end
  B = model.loads{othermode}*pinv(T');
  loads{mode,1} = Arot;
  loads{othermode,1} = B;
  
%   % Fit a two-way PARAFAC to obtain a complete structure
%   o = parafac('options');
%   o.constraints = o.constraints(1:2);
%   for j = 1:2
%     o.constraints{j}.fixed.position = ones(size(result.loads{j}));
%     o.constraints{j}.fixed.value = result.loads{j};
%     o.constraints{j}.fixed.weight = -1;
%   end
%   o.display = 'off';
%   o.waitbar = 'off';
%   o.plots = 'off';
%   result = parafac(X,size(result.loads{1},2),o);
%   result.modeltype = 'matrix factorization';
%   det = result.detail;
%   det = rmfield(det,{'stopcrit','critfinal','options','coreconsistency','algo','initialization'});
%   result.detail = det;


  result      = modelstruct('mcr');
  result.date = date;
  result.time = clock;  
  result = copydsfields(X,result);
  result.loads = loads;
  result.rotationtype = ['ORTHOMAX gamma = ',num2str(options.gamma)];
  result.rotationmatrix = T;
  result.detail.data{1}    = X;
  result.pred{1}           = result.loads{1,1}*result.loads{2,1}';
  result.detail.res{1}     = X.data(:,result.detail.includ{2,1}) - result.pred{1};
  result.ssqresiduals{1,1} = result.detail.res{1}.^2;
  result.ssqresiduals{2,1} = NaN*ones(1,size(X.data,2));
  result.ssqresiduals{2,1}(1,result.detail.includ{2,1}) = ...
  sum(result.ssqresiduals{1,1}(result.detail.includ{1,1},:),1); %var SSQs based on cal samples only
  result.ssqresiduals{1,1} = sum(result.ssqresiduals{1,1},2);    %sample SSQs for ALL samples (excluded vars already gone)

  %calculate SSQ table
  %calculate % signal by factor
  for j=1:size(result.loads{1},2);
    sig(j) = sum(sum((result.loads{1}(result.detail.includ{1,1},j)*result.loads{2}(:,j)').^2));
  end
  sig = normaliz(sig,[],1);  %normalize to 100%
  uncap = sum(sum(result.detail.res{1}(result.detail.includ{1,1},:).^2))./sum(sum(mncn(X.data(X.include{1},X.include{2})).^2));

  if uncap>=1;  %happens with really bad fits (e.g. poor constraints or badly processed data)
    uncap = sum(sum(result.detail.res{1}(result.detail.includ{1,1},:).^2))./sum(sum(X.data(X.include{1},X.include{2}).^2));
  end
  if uncap>=1;
    uncap = 0;
  end
  result.detail.ssq = [[1:size(result.loads{1},2)]' sig'*100 (1-uncap)*sig'*100 cumsum((1-uncap)*sig')*100];
  switch options.display
    case 'on'
      ssqtable(result)
  end

  if length(result.detail.includ{1,1})>size(B,2)
    result.detail.tsqlim{1,1} = tsqlim(length(result.detail.includ{1,1}),size(B,2),options.confidencelimit*100);
  else
    result.detail.tsqlim{1,1} = 0;
  end

  try
    switch lower(options.plots)
      case 'final'
        plotloads(result);
        plotscores(result);
      case 'on'
        plotloads(result);
        plotscores(result);
    end
  catch
    warning(lasterr)
  end
  
elseif strcmp(lower(model.modeltype),'tucker')
    error('Not implemented yet')
end


function [Arot,T,f]=orthomax(A,gamma);
%ORHTOMAX 
% Orthomax rotation of A
% based on Kiers (1997), psychometrika, 62, 579-598
% according to Clarkson % Jennrich (1988), psychometrika, 53(2), 251-259
% Remedy of ten Berge (1995), psychometrika, 60, 437-446 not built in
%
% [Arot,T,f]=orthomax(A,gamma);
%
% input
% A      loading-matrix to rotate
% gamma  orthomax parameter
%
% output
% Arot   rotated A
% T      rotation matrix

% Copyright © Eigenvector Research, Inc. 1991-2006
% Licensee shall not re-compile, translate or convert "M-files" contained
%  in PLS_Toolbox for use with any software other than MATLAB®, without
%  written permission from Eigenvector Research, Inc.

ConvCrit=1e-6;
[I,F]=size(A);
Arot = A;

f= sum(Arot(:).^4) - gamma/I*sum((sum(Arot.^2)).^2);
fold=2*f;
it=0;
k3 = -gamma/I;
k4 = 1;

while abs(f-fold)/abs(f)>ConvCrit
  fold=f;
  it=it+1;
  for f1=1:F-1
    for f2=f1+1:F
      x=Arot(:,f1);
      y=Arot(:,f2);

      a3 = .25*(x'*x -y'*y)^2-(x'*y)^2;
      a4 = .25*(x.^2'*x.^2 +y.^2'*y.^2)-1.5*(x.^2'*y.^2);
      b3 = x'*y*(x'*x-y'*y);
      b4 = x.^3'*y-x'*y.^3;

      a = k3*a3+k4*a4;
      b = k3*b3+k4*b4;
      theta = .25*atan2(b,a);

      xnew = x*cos(theta)+y*sin(theta);
      ynew = -x*sin(theta)+y*cos(theta);

      u = x.*x-y.*y;
      v = 2*x.*y;
      u = u-mean(u);
      v = v-mean(v);
      a = 2*I*u'*y;
      b = I*u'*u-I*y'*y;
      c = (a^2+b^2)^(.5);
      v = -sign(a)*( (b+c)/2*c  )^(.5);

      Arot(:,f1) = xnew;
      Arot(:,f2) = ynew;
    end
  end

  f= sum(Arot(:).^4) - gamma/I*sum((sum(Arot.^2)).^2);
end;

T = pinv(A)*Arot;