function [C C_axes SVD_info] = covar(varargin)
% [C C_axes] = covar(A,B,...,'axes1',A_axes,'axes2',B_axes,...,'power',power_val,'Z',Z_flag)
% [C C_axes covar_info] = covar(A,B,...,'axes1',A_axes,'axes2',B_axes,...,'power',power_val,'Z',Z_flag)
% where the first n elements are the n matrices to covary (the covariance matrix 
% output will be between the first two matrices), the rest of the elements
% are options or 'option',value pairs ... options include
%
%     'axes1', axes_A information from reading an NMR pipe file containing matrix A
%     'axes2', axes_B information from reading an NMR pipe file containing matrix B
%     'axesN', axes_N information from reading an NMR pipe file containing the n'th matrix input
%     'power', p      power by which to raise generalized covariance matrix (default is 0.5)
%     'Z', Z_flag     0 (default) for no transform, 1 for Z-matrix, 2 for 'mrs' 

% some default values
default_axes_value = zeros(27,4); %% will need to replace with more appropriate value?
default_value = -1;
non_default_value = 1;

num_args_in = nargin;

if(num_args_in < 3)
    help covar
    C = 0;
    C_axes = 0;
    SVD_info = 0;
    return;
end

start_opt_idx = num_args_in+1;

% first args (until first non-numeric arg) are input matrices ... count
% these and collect them into a cell array

for(idx=1:num_args_in)
	cur_arg = varargin{idx}; % current arg in is idx'th element of cell array varargin

	if(isnumeric(cur_arg))
		num_mats_in = idx; % number of matrices input by user
	else
		num_option_pairs = (num_args_in - num_mats_in)/2;
		start_opt_idx = idx;
		break;
    end

	input_matrices{idx} = cur_arg; % need to have this as cell array as not all matrices will be same size
end

% defaults for optional vals
axes = permute(repmat(default_axes_value,[1 1 num_mats_in]),[3 1 2]);
        % tile default axes information to have one matrix for each
        % spectrum matrix input -- need to tile in last dim then permute to
        % have desired tiling in first dim
Z_flag = 0;
power_val = 0.5;

sigma_A = 0;
sigma_B = 0;

% parse inputs beyond matrix arguments -- inputs are in the form of option name,
% value pairs

for(idx=start_opt_idx:2:num_args_in)

	opt_name = varargin{idx};
	opt_value = varargin{idx+1};
    
    opt_name_trunc = opt_name(1:min(4,length(opt_name))); % first four chars ID option type
    opt_name_idx = opt_name(5:length(opt_name)); % next chars (if any) may be indexes
                                                 % if length < 5, result is empty string 
	switch opt_name_trunc
		case {'axes'} % axes1, axes2 or ...
            axes(str2num(opt_name_idx),:,:) = opt_value;
		case {'powe'}	% powe rather than power as we only have first four chars 
			power_val = opt_value;       % (having stripped the numbers in case the result was axesN)
		case {'Z'} 
			Z_flag = opt_value;
        % add more cases here as options are expanded
	end

end

for(i=1:num_mats_in)
	cur_dim_info = squeeze(axes(i,27,:)); % vector with 1 for acceptor dims
	cur_size_row = squeeze(axes(i,26,:));
      cur_acceptor_dims = find(cur_dim_info == 1);
	cur_donor_dims = find(cur_dim_info == -1);
%    input_array = reshape(input_matrices{i}, cur_size_row'); % add extra singleton dims as code below assumes 4D tensor
      input_array_permuted = permute(input_matrices{i}, [cur_acceptor_dims' cur_donor_dims']); % place acceptor dims first
	reshaped_mat = reshape(input_array_permuted,[prod(cur_size_row(cur_acceptor_dims)) prod(cur_size_row(cur_donor_dims))]);  
                                                     % reshape array into an acceptors x donors matrix

    [m n] = size(reshaped_mat);
    m_lengths(i) = m; % potentially need to save these for output
    
    if(i > 1) % already have stacked spectrum
        if(n ~= size(stacked_spectrum, 2))
           error_line1 = 'dimensions to be covaried do not match in size:'
           error_line2 = 'spectra are incompatable for covariance'
            return; % can't run covariance in such a case
        end

		if(i == 2)
            m_0 = rows2use(2);
			cols2use = [1+m_0 (m+m_0)];
            
            if(numel(cur_acceptor_dims) == 1)
                new_C_axes = squeeze(axes(i,:,cur_acceptor_dims))';
            else
                C_axes = squeeze(axes(i,:,cur_acceptor_dims));
            end
            new_C_axes(27,:) = -1; % new "donor" dims are acceptor dims of second spectrum input
            C_axes = [C_axes new_C_axes];
        end
    
        sigma_B = median(abs(reshaped_mat(ceil((n*m).*rand(1000,1)))))/norminv(0.75);
        stacked_spectrum = [stacked_spectrum; reshaped_mat];

    else % need to initialize stacked_spectrum -- also initialize rows2use and C_axes
      sigma_A = median(abs(reshaped_mat(ceil((n*m).*rand(1000,1)))))/norminv(0.75);
      stacked_spectrum = reshaped_mat;
	  rows2use = [1 m];
      if(numel(cur_acceptor_dims) == 1)
        C_axes = squeeze(axes(i,:,cur_acceptor_dims))';
      else
        C_axes = squeeze(axes(i,:,cur_acceptor_dims));
      end
    end
    %% should we clean up vars we don't need or let them persist?
end

cur_num_dims = size(C_axes,2);
if(num_mats_in == 1)
    cols2use = rows2use;
    C_axes = [C_axes C_axes];
    C_axes(27,(cur_num_dims+1):4) = -1;
    cur_num_dims = 2*cur_num_dims;
    sigma_B = sigma_A;
end
C_axes(:,(cur_num_dims+1):4) = 0; % pad C_axes to complete axis info array
C_axes(26,(cur_num_dims+1):4) = 1; % but pad dim counts to one
C_axes(27,(cur_num_dims+1):4) = -1; % and pad axes as donor axes

[U S V] = svd(stacked_spectrum, 'econ'); % SVD implimentation of (generalized) covariance NMR
Nd = size(stacked_spectrum, 2);

svd_power = power_val*2; % e.g. power_val = 0.5 covariance corresponds to U*S^1*U'

numel_C = prod(max(C_axes(26,:),1));

if(numel_C <= (2048*2048)) % i.e. less than 32 MB ... we should have enough memory for a 32 MB matrix  
	if(Z_flag == 2)
        stacked_spectrum = max(abs(stacked_spectrum),min(sigma_A,sigma_B)/10000); % make sure spectra are positive
        Amax = max(stacked_spectrum(rows2use(1):rows2use(2),:),[],2);
        Bmax = (max(stacked_spectrum(cols2use(1):cols2use(2),:),[],2)').^(-1);
    end
    clear stacked_spectrum; % we don't need this anymore, save some memory
    
    if(nargout == 3) % save SVD info: should we check for memory here?
		SVD_info.rows2use = rows2use;
		SVD_info.cols2use = cols2use;
		SVD_info.U = U;
		SVD_info.V = V;
		SVD_info.S = S;
		SVD_info.svd_power = svd_power;
        SVD_info.Z_flag = Z_flag;
	end

	clear V; % we don't need it anymore, so free memory
	US = U(rows2use(1):rows2use(2),:)*(S.^power_val);
	USprime = (U(cols2use(1):cols2use(2),:)*(S.^power_val))';
    % Uprime = (U(cols2use(1):cols2use(2),:))';
    
    clear S;
    clear U;

	C = US*USprime;
      %% should we reshape C to have the correct tensor-rank, etc.?
    if(Z_flag == 1)
        var_A = sigma_A^2;
        var_B = sigma_B^2;
        nn = -Nd*var_A*var_B;
        Adiag = diag(US*US');
        Bdiag = diag(USprime'*USprime)';
        clear US;
        clear USprime;
        if(abs(power_val - 1.0) < 0.001)
            C = C./sqrt(nn + var_B*repmat(Adiag,1,size(C,2)) + var_A*repmat(Bdiag,size(C,1),1));
        elseif(abs(power_val - 0.5) < 0.001)
            var_C1 = nn + var_B*repmat(Adiag,1,size(C,2)) + var_A*repmat(Bdiag,size(C,1),1);
            linear_factor = 1./( repmat(Adiag.^(0.5),1,size(C,2))  +  repmat(Bdiag.^(0.5),size(C,1),1) );
            if(max(abs(rows2use-cols2use)) < 1) % C has a diagonal in it, so we must subtract that
                C = C - diag(diag(C));
            end
            sigma_C2 = linear_factor.*sqrt(var_C1 - 2*((linear_factor.*var_C1)*(linear_factor.*(C*C))));
            C = C./sigma_C2;
        else
            error_line3 = 'Z transform currently supported only for \lambda = 1 or \lambda = 0.5'
        end
    elseif(Z_flag == 2)
        clear US;
        clear USprime;
        W = exp(abs(log(Amax*Bmax)));
        C = C./W;
    end
    
else
    clear stacked_spectrum; % we don't need this anymore, save some memory
    
    file_key = rand*100;
	C = ['covar_SVD_' num2str(ceil(file_key)) '.tmp']; % C is name of file with all info 
	temp_file_out = fopen(C,'w');         % instead of covar result

	% format: %
	% header: 7 32 bit floating point numbers containing size of S U, V and C_axes then svd_power then Z_flag, sigma_A and sigma_B %
	%         followed by C_axes information for desired covar spec %
	%         then 4 32 bit floats indicating ranges of rows and cols of U to use for covar recon. %
	% S (note this is first!) %
	% U %
	% V %

	fwrite(temp_file_out, size(S),     'float32');
	fwrite(temp_file_out, size(U),     'float32');
	fwrite(temp_file_out, size(V),     'float32');
	fwrite(temp_file_out, size(C_axes),'float32'); 

	fwrite(temp_file_out, svd_power,   'float32');
    fwrite(temp_file_out, Z_flag,   'float32');
    fwrite(temp_file_out, sigma_A,   'float32');
    fwrite(temp_file_out, sigma_B,   'float32');
    
	fwrite(temp_file_out, C_axes(:),'float32'); 

	fwrite(temp_file_out, rows2use, 'float32');
	fwrite(temp_file_out, cols2use, 'float32');

	fwrite(temp_file_out, S(:), 'float32');
	fwrite(temp_file_out, U(:), 'float32');
	fwrite(temp_file_out, V(:), 'float32');

	fclose(temp_file_out);
    
    SVD_info = C_axes;
end