LORENE
poly_regression.C
1 /*
2  * Functions to compute polynomial regression of 1 data.
3  *
4  */
5 
6 /*
7  * Copyright (c) 2022 Jerome Novak
8  *
9  * This file is part of LORENE.
10  *
11  * LORENE is free software; you can redistribute it and/or modify
12  * it under the terms of the GNU General Public License as published by
13  * the Free Software Foundation; either version 2 of the License, or
14  * (at your option) any later version.
15  *
16  * LORENE is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  * GNU General Public License for more details.
20  *
21  * You should have received a copy of the GNU General Public License
22  * along with LORENE; if not, write to the Free Software
23  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24  *
25  */
26 
27 
28 
29 
30 /*
31  * $Id: poly_regression.C,v 1.1 2022/04/15 13:39:25 j_novak Exp $
32  * $Log: poly_regression.C,v $
33  * Revision 1.1 2022/04/15 13:39:25 j_novak
34  * New class Eos_compose_fit to generate fitted EoSs from CompOSE tables.
35  *
36  *
37  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Utilities/poly_regression.C,v 1.1 2022/04/15 13:39:25 j_novak Exp $
38  *
39  */
40 
41 // Lorene headers
42 #include "matrice.h"
43 
44 namespace Lorene {
45 
46  // Function to initialize the transformation matrix Chebyshev -> Taylor expansion
47  //-------------------------------------------------------------------------------
48  Tbl initialize_dd(int ncoef) {
49  Tbl dd(ncoef, ncoef) ;
50  dd.annule_hard() ;
51  dd.set(0,0) = 1. ; dd.set(0, 1) = 0 ;
52  dd.set(1, 0) = 0 ; dd.set(1, 1) = 1. ;
53  for (int i=2; i<ncoef; i++) {
54  dd.set(i,0) = 0 ;
55  dd.set(i, 1) = 0 ;
56  }
57  for (int j=2; j<ncoef; j++) {
58  dd.set(0, j) = -dd(0, j-2) ;
59  for (int i=1; i<ncoef; i++) {
60  dd.set(i, j) = 2.*dd(i-1, j-1) - dd(i, j-2) ; // Recursion formula
61  // of Chebyshev polynomials
62  }
63  }
64 
65  return dd ;
66  }
67 
68 
69  // Computes the polynomial regression of degree n_poly to data (xx, yy)
70  // Result is an array of Chebyshev coefficients.
71  Tbl poly_regression(const Tbl& xx, const Tbl& yy, int n_poly) {
72 
73  assert((xx.get_ndim() == 1) && (yy.get_ndim() == 1)) ;
74  assert( xx.get_taille() == yy.get_taille() ) ;
75 
76  int n_data = xx.get_dim(0) ;
77 
78  Tbl powxj(n_data, 2*n_poly+1) ; powxj.set_etat_qcq() ;
79  for (int i=0; i<n_data; i++)
80  powxj.set(i, 0) = 1. ;
81  for (int j=1; j<2*n_poly+1; j++) {
82  for (int i=0; i<n_data; i++) {
83  powxj.set(i,j) = powxj(i, j-1)*xx(i) ;
84  }
85  }
86 
87  Tbl xpow(2*n_poly+1) ; xpow.annule_hard() ;
88  for (int j=0; j<2*n_poly+1; j++) {
89  for (int i=0; i<n_data; i++) {
90  xpow.set(j) += powxj(i, j) ;
91  }
92  }
93 
94  Tbl rhs(n_poly+1) ; rhs.annule_hard() ;
95  for (int j=0; j<n_poly+1; j++) {
96  for (int i=0; i<n_data; i++) {
97  rhs.set(j) += powxj(i, j)*yy(i) ;
98  }
99  }
100 
101  Matrice mat(n_poly+1, n_poly+1) ; mat.set_etat_qcq() ;
102  mat.set(0,0) = n_data ;
103  for (int i=0; i<=n_poly; i++) {
104  for (int j=0; j<=n_poly; j++) {
105  mat.set(i,j) = xpow(i+j) ;
106  }
107  }
108  mat.set_lu() ;
109 
110  Tbl sol_Taylor = mat.inverse(rhs) ; // Solution expressed as Taylor coefficients
111 
112  Matrice cheb = initialize_dd(n_poly+1) ;
113  cheb.set_lu() ;
114 
115  return cheb.inverse(sol_Taylor) ;
116  }
117 
118 }
Tbl poly_regression(const Tbl &, const Tbl &, int)
Polynomial regression, giving Chebyshev coefficients.
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition: matrice.C:395
Lorene prototypes.
Definition: app_hor.h:67
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:427
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
Definition: tbl.h:420
Matrix handling.
Definition: matrice.h:152
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:423
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:417
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:178
Basic array class.
Definition: tbl.h:164
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:375