% FACTORS: Estimates loadings of a set of secondary factors of a covariance
%          matrix, given a primary factor and a corresponding boolean matrix
%          of specifying variables included in the secondary factor(s).
%          Secondary factors by default are orthogonal (with the primary factor
%          and with one another), but this constraint can be relaxed.
%          Also returns the residual covariance matrix.
%          Error handling is done by function WRIGHT.
%
%     Syntax: [fact,rc] = factors(c,pfactor,secnd,orthog,suppress)
%
%         c =       [n x n] residual covariance or correlation matrix.
%         pfactor = n-length vector of loadings of primary factor.
%         secnd =     [n x s] matrix of boolean values (T/F = 1/0)
%                     indicating by 1's the submatrices of variables
%                     to be included in the secondary factors (s of them).
%         orthog =  flag indicating whether secondary factors are to be:
%                     0 - oblique to primary factor and one another;
%                     1 - orthogonal to primary factor but oblique to other
%                           secondary factors;
%                     2 - orthogonal to primary factor and all other secondary
%                           factors;
%                     [default = 2].
%         suppress = boolean flag indicating whether warning messages
%                     are to be suppressed (invoked for bootstrapping and
%                     randomization).
%         ----------------------------------------------------------------
%         fact =    [n x s+1] matrix of factor loadings (primary + secondary).
%         rc =      [n x n] matrix of residual covariances.
%

function [fact,rc] = factors(c,pfactor,secnd,orthog,suppress)
  TRUE = 1; FALSE = 0;

  [nvar,p] = size(c);
  [n,nsec] = size(secnd);

  % Estimate loadings from covariance matrix.  Initial guess is that
  % loadings are equal to the root-mean-square covariance
  % (including variance), by column, orthogonal to other factors.

  % First pass: optimize secondary factors separately with respect to
  % primary factor only.

  fact = zeros(nvar,nsec);
  init_load = [];

  for s = 1:nsec                      % Initial loadings for secondary factors
    secvect = find(secnd(:,s));           % Transform current factor into var list

    if (length(secvect)==2)             % 2 variables: exact solution
      cc = c(secvect(1),secvect(2));      % Covar between the 2 variables

      if (orthog==0)                      % If factor can be oblique to primary,
        f = ones(2,1) * sqrt(abs(cc));      % Each loading is sqrt of covar
        if (cc<0)                           % If covar is negative, one loading
          f(2) = -f(2);                     %   is negative
        end;
      else                                % If should be orthogonal to primary,
        a1 = pfactor(secvect(1));           % Exact solution possible only if:
        a2 = pfactor(secvect(2));           %   primary factor is general,
        caa = -cc*a2/a1;                    %   secondary factor is bipolar.
        if (caa > 0)
          f = ones(2,1);
          f(1) = sqrt(caa);
          f(2) = cc / f(1);
        else                                % But not possible if covar is pos,
          f = ones(2,1) * sqrt(abs(cc));      % so each loading is sqrt of covar
          if (~suppress)
            disp('    Warning: 2-variable secondary factor must be oblique');
            disp('             to primary factor when residual covariance');
            disp('             is positive.');
            suppress = TRUE;
          end;
        end;
      end;

    else                                % >2 variables: estimated solution
      redc = c(secvect,secvect);          % Reduced covar matrix for this factor
      g = min(sqrt(abs(mean(redc)))',sqrt(diag(redc)) );         % Initial loadings
      f = fmins('factorsf',g,[],[],c,pfactor,secnd(:,s),orthog); % Optimize this factor
    end;

    fact(secvect,s) = f;                % Insert into factor-loadings matrix
    init_load = [init_load;f];          % Append onto initial-loadings vector
  end;

  % Optional second pass: simultaneously optimize 2+ secondary factors if they
  % are intercorrelated.  Secondary factors of two variables having loadings
  % of the same sign are omitted from second-pass optimization because they
  % must be oblique to any primary factor.

  if (orthog)
    ns = sum(secnd);                      % Omit sec factors with 2 vars of same sign

    omit = [];
    keep = [];
    for s = 1:nsec
      if (ns(s) > 2)
        keep = [keep s];
      else
        if (abs(sum(sign(fact(:,s)))))
          omit = [omit s];
        else
          keep = [keep s];
        end;
      end;
    end;

    if (nsec-length(omit) > 1)          % Simultaneously optimize >1 factor
      if (length(omit)>0)                 % If 2-var factors to be omitted,
        nsec_copy = nsec;                 %   Make copies of full matrices
        sec_copy = secnd;
        fact_copy = fact;

        nsec = nsec - length(omit);       %   Delete oblique 2-var factors
        secnd(:,omit) = [];
        fact(:,omit) = [];

        factors_omitted = TRUE;
      else
        factors_omitted = FALSE;
      end;

      vc = abs(vectcorr(fact));           % Vector corrs among secondary factors
      sumcorr = sum(trilow(vc));

      if (sumcorr > eps)                  % Optimize if any corrs > 0
        f = fmins('factorsf',init_load,[],[],c,pfactor,secnd,orthog);
      end;

      low = 1;                            % Unpack loadings vector into matrix
      for s = 1:nsec                      %   of secondary-factor loadings
        secvect = find(secnd(:,s));
        high = low+length(secvect)-1;
        fact(secvect,s) = f(low:high);
        low = high+1;
      end;
    end;

    if (factors_omitted)                % If 2-var factors have been omitted,
      nsec = nsec_copy;                 %   restore them
      secnd = sec_copy;
  
      fact_copy(:,keep) = fact;
      fact = fact_copy;
    end;
  end;  % if orthog

%  r2 = vectcorr(pfactor,fact).^2;     % Independent variances of secondary factors
%  rc = c;                             % Residual covariances
  for s = 1:nsec
    f = fact(:,s);
    ff = f*f';                          % Reproduced covariances
%    for i = 1:nvar
%      if (abs(f(i))>eps)
%       rc(i,i) = rc(i,i) - ff(i,i);    % Fully reproduce diagonal elements
%        ff(i,i) = 0;                    % Conditionally reproduce off-diag elements
%      end;
%    end;
%    rc = rc - ff*(1-r2(s));
    rc = c - ff;
  end;

  fact = [pfactor fact];              % Return primary + secondary factor loadings

  return;
