//  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 [is_equal] = CL_isAlmostEqual(A, B, rtol, atol, crit)
// Value comparison with tolerance
//
// Calling Sequence
// [is_equal] = CL_isAlmostEqual(A, B [,rtol, atol, crit])
//
// Description
// <itemizedlist><listitem>
// <p>This function returns %t if a value (A) is equal to a reference value (B) to within a tolerance 
// tol defined by:</p>
// <p>- tol = rtol*max_val + atol,</p>
// <p>where <b>rtol</b> and <b>atol</b> are the relative and absolute tolerance, respectively, 
// and max_val depends on the comparison criterion <b>crit</b> defined by:</p>
// <p>- crit = "element" => the comparison will be done element by element: 
// max_val = max(abs(A),abs(B)) </p>
// <p>- crit = "CLnorm" => the comparison will be done on the norms of the columns vectors: 
// max_val = max(CL_norm(A),CL_norm(B))</p>
// <p></p>
// <p><b>Notes:</b></p>
// <p>- A and B can be matrices (type="constant"), hypermatrices (type="hypermat"), structures or lists 
// (type="struct", "list", "tlist" or "mlist"). </p>
// <p>- This function only works for "double" values (booleans or strings are not handled). </p>
// <p>- If A and B are of different type, different size, or have different field names, %f is returned. </p>
// <p>- Hypermatrices can only be compared using the "element" comparison criterion. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// A: First value used in comparison (the result)
// B: Second value used in comparison (the reference)
// rtol: (optional) Relative precision. Default is 1.e-14 if atol is 0, 0 otherwise.  (1x1)
// atol: (optional) Absolute precision. Default is 0. (1x1)
// crit: (optional, string) Comparison criterion: "element" or "CLnorm". Default is "element".
// is_equal: (boolean) %t if A and B are equal, %f otherwise (1x1)
// 
// Authors
// CNES - DCT/SB
//
// Examples
// A = 1.E-6; 
// B = A + 1.e-12; 
// CL_isAlmostEqual(A, B, rtol=2.e-6) // => %t
// CL_isAlmostEqual(A, B, rtol=5.e-7) // => %f



// Internal function to compare two matrix (type "constant")
// or two hypermatrix (type "hypermat")
// Hypermatrix are only compared with "element" comparaison criterion
function [is_equal] = CL__isAlmostEqual_m_hm(A, B, rtol, atol, crit)

  sA = size(A);
  sB = size(B);
  // If A and B do not have the same dimensions:
  if (~isequal(sA,sB)); is_equal = %f; return; end;  // <-- RETURN !
  
  // If %nan do not match in A and B:
  if (~isequal(find(isnan(A)), find(isnan(B)))) 
    is_equal = %f; 
    return; // <-- RETURN !
  end
  
  // If A and B are empty, they are equal
  if (A == []); return; end; // <-- RETURN !
  
  // Function for comparison
  if (crit == "element")
    // Nan values are not considered
    I = find(~isnan(A));
    err = matrix(abs(A(I)-B(I)),1,-1);
    max_val = matrix(max(abs(A(I)),abs(B(I))),1,-1);
    I1 = find((A(I)==0 & B(I)<>0) | (A(I)<>0 & B(I)==0));
    if (I1 <> [] & atol == 0)
        max_val(I1) = 1;
        CL__warning("Comparison with 0 => Absolute precision used instead of relative");
    end
    
  elseif (crit == "CLnorm")
    // Nan values are not considered
    I = find(~isnan(CL_norm(A)));
    err = matrix(CL_norm(A(:,I)-B(:,I)),1,-1);
    max_val = matrix(max(CL_norm(A(:,I)),CL_norm(B(:,I))),1,-1);
    I1=find((CL_norm(A(:,I))==0 & CL_norm(B(:,I))<>0) | (CL_norm(A(:,I))<>0 & CL_norm(B(:,I))==0));
    if (I1 <> [] & atol == 0)
        max_val(I1) = 1;
        CL__warning("Comparison with 0 => Absolute precision used instead of relative");
    end
    
  end

  tol = rtol .* max_val + atol;
  if (find(err > tol) <> []); is_equal = %f; end;

endfunction


// Internal function for comparing:  
// - 2 structure (type=17)
// - 2 mlist (type=17)
// - 2 tlist(type=16)
// - 2 list(type=15)
// - NB: typeof returns the first string in the first list entry for tlist and mlist
//       -> we use type function instead
function [is_equal] = CL__isAlmostEqual_list(A, B, rtol, atol, crit)
  // Structure, mlist or tlist
  if (type(A) == 17 | type(A) == 16)
    fA = fieldnames(A); // column vector
    fB = fieldnames(B); // column vector
    
    // If field names do not match
    if (~isequal(fA,fB))
      is_equal = %f;  
      return; // <-- RETURN
    end;
    
    for (k = 1 : size(fA,1))
      is_equal = CL_isAlmostEqual(A(fA(k)), B(fA(k)), rtol, atol, crit); // recursive call
      if (~is_equal); return; end
    end

  // List
  elseif (typeof(A) == "list")
    sA = size(A);
    sB = size(B);
    if (sA <> sB) 
      is_equal = %f; 
      return; // <-- RETURN
    end
    
    for (k = 1 : sA)
      is_equal = CL_isAlmostEqual(A(k), B(k), rtol, atol, crit); // recursive call
      if (~is_equal); return; end
    end
  end
endfunction



// ---------------------
// MAIN FUNCTION
// ---------------------

if (argn(2) < 2 | argn(2) > 5)
  CL__error("Wrong number of input arguments");
end

if (~exists("rtol","local")); rtol = 0; end;
if (~exists("atol","local")); atol = 0; end;
if (~exists("crit","local")); crit = "element"; end;

if (rtol == 0 & atol == 0); rtol = 1.e-14; end;

if (crit <> "element" & crit <> "CLnorm")
  CL__error("Invalid value for crit");
end

if (typeof(A) <> typeof(B)); 
  is_equal = %f;  
  CL__warning("Comparison of inconsistent types");
  return;  // < -- RETURN!
end

is_equal = %t;

// ---------------------
// Matrix or hypermatrix
// ---------------------
if (typeof(A) == "constant" | typeof(A) == "hypermat")
   if ((typeof(A) == "hypermat") & (crit == "CLnorm"))  then
       CL__error("Comparison criterion not useable with hypermatrices");
   end
  [is_equal] = CL__isAlmostEqual_m_hm(A,B,rtol,atol,crit);
  
// -----------------------------------------------------------------------
// Structure(type=17) or tlist(type=16) or mlist(type=17) or list(type=15)
// -----------------------------------------------------------------------
// Note : hypermatrix actually is a mlist (type=17) and could also be handled here
//        but it is dealt with in the first "if"
elseif (type(A) == 17 | type(A) == 16 | typeof(A) == "list")
  [is_equal] = CL__isAlmostEqual_list(A,B,rtol,atol,crit);
  
else
  CL__error("Types not handled in comparison");
end

endfunction