function [B]=misscorrmatrix(A,mode);
%B=misscorrmatrix(A,mode);
%
% Claus A. Andersson, claus@andersson.dk, Copyrighted (C) 98
%
% This file derives a correlation matrix of a matrix that has
% missing values. Note that all columns must have at least two 
% non-missing observations in common with all other columns.
% Otherwise the returned correlation mapping has nans
% on the corresponding positions.
%
% A      :  Matrix (Samples x Variables) with missing values.
% mode   :  Defines the mode of operation
%           0: Correlation (std&mean), [-1:1] (default)
%           1: Covariances (mean), [-inf:inf]
%
% Note: If any variables are constant over the samples, the results
%       will have nans in them.  
%
% Two steps to make a mapping:
%   B=misscorrmatrix(A); %Step 1
%   imagesc(B);colorbar;grid on; %Step 2


if ~exist('mode'),
    mode=0;
end;

i_A=size(A,1);
j_A=size(A,2);
isnanA=isnan(A);
sumisnanA=sum(isnanA);
B=nan*zeros(j_A);

nan_list=find(sumisnanA);
true_list=find(~sumisnanA);

for i=1:j_A,
    true_ivec=~isnanA(:,i);
    for j=i+1:j_A,
        true_jvec=~isnanA(:,j);
        list=true_ivec.*true_jvec;
        divis=1/sum(list);
        idx=find(list);
              
        ivec=A(idx,i);
        jvec=A(idx,j);
        
        ivec=ivec-sum(ivec)*divis;
        jvec=jvec-sum(jvec)*divis;
        
        if mode==0,
            ivecsum=sum(ivec.^2);
            if ivecsum>eps,
                ivec=ivec/sqrt(ivecsum);
            else
                ivec=ivec*0;
            end;
            
            jvecsum=sum(jvec.^2);
            if jvecsum>eps,
                jvec=jvec/sqrt(jvecsum);
            else
                jvec=jvec*0;
            end;
        end;
        
        B(i,j)=ivec'*jvec;
    end;
end;

