//***********************************************************
//
//  NXV4 "Math-on-Paper" software
//
//  package : MTX class _____________________________________
//            
//            supports: move semantix
//                      operators + - * / L,R-val context
//                                += -= *= /=
//                      fabs()
//                      ~  -      unary operators L,R-val
//                      mtx<arithmetic> i.e.     int
//                                               float 
//                                               double
//                                               real
//                                           cpx<int> 
//                                           cpx<float> 
//                                           cpx<double> 
//                                           cpx<real> 
//                      type_casts between the types 
//                      cout << overload 
//
//                      TOTAL = 5084 instantiated op/f's
//                      _____________________________________
//
//
//  authors : Mihai-Tiberiu Dima (1)
//            Maria Dima         (2)
//
//             Meshcheryakov Laboratory of Information Techn.
//             Joint Institute for Nuclear Research
//             Dubna-Russia
//
//            on leave from:
//
//             Hyperion Univ. Bucharest-Romania
//
//
//  date    : MTX4/v1.0.alfa   / Mon Mar 07 12:45:54 CET 2022
//
//***********************************************************
#include <cmath>
#include <iostream>

#include "cpx.hh"
#include "vec.hh"
#include "mtx.hh"

using namespace std                                         ;
using namespace cpx4                                        ;
using namespace vec4                                        ;
using namespace mtx4                                        ;


// application for linear regression
//
// M |u> = M_uu |u>
//
// where:
//  - M_ij = < dx_i dx_j >
//  - |u>  = line direction
//  - M_uu = long axis of nebula
//
// to use this example do the following:
//   rm test.cc
//   ln -s appl.cc test.cc
//   make test
//   make run

  using fltx = float                                        ;
  using dblx = double                                       ;
  using real = long double                                  ;

int main()
{

 // vec4::utf8x = true                                         ;
 // set of points, average = 0

 auto r1 = vec<dblx>( 0.1,  1.2,  3.0)                      ;
 auto r2 = vec<dblx>( 0.0, -1.3, -3.1)                      ;
 auto r3 = vec<dblx>(-0.2, -2.4, -6.0)                      ;
 auto r4 = vec<dblx>( 0.1,  2.5,  6.1)                      ;
 auto r5 = vec<dblx>( 0.4,  4.8, 12.0)                      ;
 auto r6 = vec<dblx>(-0.2, -4.6,-12.4)                      ;
 auto r7 = vec<dblx>( 0.5,  6.0, 15.0)                      ;
 auto r8 = vec<dblx>(-0.7, -6.2,-15.4)                      ;

 auto M  = (r1^r1) + (r2^r2) + (r3^r3) + (r4^r4) +
           (r5^r5) + (r6^r6) + (r7^r7) + (r8^r8)            ;

 M /= 8                                                     ;

 auto E  = eigen(M)                                         ;

 vec<dblx> u1(E.x11, E.x21, E.x31)                          ;
 vec<dblx> u2(E.x12, E.x22, E.x32)                          ;
 vec<dblx> u3(E.x13, E.x23, E.x33)                          ;

 auto L1 = ( u1 | M | u1 )                                  ;
 auto L2 = ( u2 | M | u2 )                                  ;
 auto L3 = ( u3 | M | u3 )                                  ;

 auto         LL = L1                                       ;
 auto         LA = L2                                       ;
 auto         LB = L3                                       ;
 auto         uu = u1                                       ;
 if (L2>LL)  {LA = LL                                       ;
              LL = L2                                       ;
              uu = u2                                       ;}
 if (L3>LL)  {LB = LL                                       ;
              LL = L3                                       ;
              uu = u3                                       ;}
 
  cout                                             << endl  ;
  cout << noboolalpha                              << endl  ;

  cout                                             << endl  ;
  cout << " _____________________________________" << endl  ;
  cout                                             << endl  ;
  cout << " The nebula has long axis " << sqrt(LL) << endl  ; 
  cout << " along direction "          << uu       << endl  ;
  cout                                             << endl  ;
  cout << " The nebula waist is " << sqrt(LA) << " x "
                                  << sqrt(LB)      << endl  ; 
  cout << " _____________________________________" << endl  ;
  cout                                             << endl  ;
  cout                                             << endl  ;
  cout                                             << endl  ;
  cout                                             << endl  ;
  cout << boolalpha                                << endl  ;
  cout << " The nebula mtx is M = " << M << endl   << endl  ; 
  cout                                             << endl  ;
  cout                                             << endl  ;
  cout                                             << endl  ;
  cout << noboolalpha                              << endl  ;
  cout << "("
       << cpx<int >(1,1)   << ") * (" 
       <<  vec<dblx>(1,2,0) << ") = "
       << (cpx<int >(1,1)) * vec<dblx>(1,2,0)      << endl  ;
  cout                                             << endl  ;
  cout                                             << endl  ;
  cout << boolalpha                                << endl  ;
  cout <<  cpx<int >(1,1)   << " * " 
       <<  vec<dblx>(1,2,0) << " = "
       << (cpx<int >(1,1)) * vec<dblx>(1,2,0)      << endl  ;
  cout                                             << endl  ;
  cout                                             << endl  ;
  cout                                             << endl  ;
 
}