//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

function [corm, sd] = CL_cov2cor(covm)
// Covariance to correlation matrix
//
// Calling Sequence
// [corm, sd] = CL_cov2cor(covm)
//
// Description
// <itemizedlist><listitem>
// <p> Computes the correlation matrix and the standard deviation vector from the covariance matrix.</p> 
// <p> The correlation matrix and standard deviation vector are defined as follows: </p>
// <p> - <b>corm</b>(i,i) = 1 </p>
// <p> - <b>corm</b>(i,j) = <b>covm</b>(i,j) /  (<b>sd</b>(i) * <b>sd</b>(j)) </p>
// <p> - <b>sd</b>(i) = covm(i,i)^(1/2) </p>
// <p></p>
// <p><b>Notes</b>:</p>
// <p>- The covariance matrix should be symetrical (not checked). </p>
// <p>- An error is raised if one diagonal term is not positive. </p>
// <p>- If one diagonal term of the covariance matrix is 0, the corresponding diagonal term of the 
//      correlation matrix is 1; other terms of the correlation matrix on the same row or column
//      are set to %nan. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// covm: Covariance matrix (NxNxK)
// corm: Correlation matrix (NxNxK)
// sd: Standard deviation vector (NxK)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_cor2cov
//
// Examples
// covm = [1, 0, 1.5; 0, 4, -4.8; 1.5, -4.8, 9];  
// [corm, sd] = CL_cov2cor(covm);
// 
// // Consistency check: 
// covm2 = CL_cor2cov(corm, sd) // => covm
// 

// Covariance to correlation matrix (2d case)

function [corm,sd] = cov2cor(covm)
    
  diagcovm = diag(covm);
  
  I = find(diagcovm < 0);
  if (I <> []); CL__error("Invalid term in covariance matrix"); end

  Izero = find(diagcovm == 0);
  // insert %nan to avoid division by 0
  diagcovm(Izero) = %nan;
    
  w = diag(1 ./ sqrt(diagcovm));
  // Correlation matrix
  corm = w * covm * w;
    
  // force 1 on the diagonal 
  // as in R (see https://svn.r-project.org/R/trunk/src/library/stats/R/cor.R)
  N = size(covm,1); 
  corm(sub2ind([N,N], 1:N, 1:N)) = 1;
    
  // standard deviation always defined
  sd = sqrt(diagcovm);
  sd(Izero) = 0;
  
endfunction

    
// Argument checking

// special case: 
// output = [] if input = []
if (covm == [])
  corm = []; 
  sd = []; 
  return; 
end
 
// add 1 => also works for matrices
s = [size(covm), 1];
 
if (s(1) <> s(2))
  CL__error("Input matrix not square");
end
  

if (s(3) == 1)
  // Matrix    
  [corm,sd] = cov2cor(covm);   
    
else
  // Hypermatrix
  corm = zeros(covm);
  N = size(covm, 3);
  sd = zeros(s(1), N);
    
  // loop on 3rd dimension
  for (i = 1 : N)    
    [cormi, sdi] = cov2cor(covm(:,:,i));
    corm(:,:,i) = cormi;
    sd(:,i) = sdi;      
  end    
end
  
endfunction