Skip to content

Commit

Permalink
refactoring
Browse files Browse the repository at this point in the history
  • Loading branch information
@ebertolazzi
ebertolazzi committed Feb 28, 2024
1 parent 5fa9220 commit a2f4ed4
Show file tree
Hide file tree
Showing 35 changed files with 622 additions and 634 deletions.
34 changes: 17 additions & 17 deletions src/SplineAkima.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,20 +74,20 @@ namespace Splines {

void
Akima_build(
real_type const * X,
real_type const * Y,
real_type * Yp,
integer npts
real_type const X[],
real_type const Y[],
real_type Yp[],
integer npts
) {

if ( npts == 2 ) { // solo 2 punti, niente da fare
Yp[0] = Yp[1] = (Y[1]-Y[0])/(X[1]-X[0]);
} else {
Malloc_real mem("Akima_build");
real_type * m = mem.malloc( size_t(npts+3) );
real_type * m{ mem.malloc( size_t(npts+3) ) };

// calcolo slopes (npts-1) intervals + 4
for ( size_t i = 1; i < size_t(npts); ++i )
for ( size_t i{1}; i < size_t(npts); ++i )
m[i+1] = (Y[i]-Y[i-1])/(X[i]-X[i-1]);

// extra slope at the boundary
Expand All @@ -97,16 +97,16 @@ namespace Splines {
m[size_t(npts+2)] = 2*m[size_t(npts+1)]-m[size_t(npts)];

// minimum delta slope
real_type epsi = 0;
for ( size_t i = 0; i < size_t(npts+2); ++i ) {
real_type epsi{0};
for ( size_t i{0}; i < size_t(npts+2); ++i ) {
real_type dm = std::abs(m[i+1]-m[i]);
if ( dm > epsi ) epsi = dm;
}
epsi *= 1E-8;

// 0 1 2 3 4---- n-1 n n+1 n+2
// + + + + +
for ( size_t i = 0; i < size_t(npts); ++i )
for ( size_t i{0}; i < size_t(npts); ++i )
Yp[i] = akima_one( epsi, m[i], m[i+1], m[i+2], m[i+3] );
}
}
Expand All @@ -117,14 +117,14 @@ namespace Splines {

void
AkimaSpline::build() {
string msg = fmt::format("AkimaSpline[{}]::build():", m_name );
string msg{ fmt::format("AkimaSpline[{}]::build():", m_name ) };
UTILS_ASSERT(
m_npts > 1, "{} npts = {} not enought points\n", msg, m_npts
);
Utils::check_NaN( m_X, (msg+" X ").c_str(), m_npts, __LINE__, __FILE__ );
Utils::check_NaN( m_Y, (msg+" Y ").c_str(), m_npts, __LINE__, __FILE__ );
integer ibegin = 0;
integer iend = 0;
integer ibegin{0};
integer iend{0};
do {
// cerca intervallo monotono strettamente crescente
while ( ++iend < m_npts && m_X[iend-1] < m_X[iend] ) {}
Expand All @@ -149,17 +149,17 @@ namespace Splines {
// gc["ydata"]
//
*/
string where = fmt::format("AkimaSpline[{}]::setup():", m_name );
GenericContainer const & gc_x = gc("xdata",where.c_str());
GenericContainer const & gc_y = gc("ydata",where.c_str());
string where{ fmt::format("AkimaSpline[{}]::setup():", m_name ) };
GenericContainer const & gc_x{ gc("xdata",where.c_str()) };
GenericContainer const & gc_y{ gc("ydata",where.c_str()) };

vec_real_type x, y;
{
std::string ff = fmt::format( "{}, field `xdata'", where );
std::string ff{ fmt::format( "{}, field `xdata'", where ) };
gc_x.copyto_vec_real( x, ff.c_str() );
}
{
std::string ff = fmt::format( "{}, field `ydata'", where );
std::string ff{ fmt::format( "{}, field `ydata'", where ) };
gc_y.copyto_vec_real( y, ff.c_str() );
}
this->build( x, y );
Expand Down
92 changes: 46 additions & 46 deletions src/SplineAkima2D.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,10 +113,10 @@ namespace Splines {
real_type volatility_factor = dz0*dz0 + dz1*dz1 + dz2*dz2 + dz3*dz3;
// epsi value used to decide whether or not
// the volatility factor is essentially zero.
real_type epsi = (z0*z0+z1*z1+z2*z2+z3*z3)*1.0E-12;
real_type epsi{ (z0*z0+z1*z1+z2*z2+z3*z3)*1.0E-12 };
// Accumulates the weighted primary estimates and their weights.
if ( volatility_factor > epsi ) { // finite weight.
real_type WT = 1.0/ (volatility_factor*SXX);
real_type WT{ 1.0/ (volatility_factor*SXX) };
RF[1] += WT*primary_estimate;
RF[0] += WT;
} else { // infinite weight.
Expand Down Expand Up @@ -146,22 +146,22 @@ namespace Splines {

real_type by0[4], by1[4], CY1A[4], CY2A[4], CY3A[4], SYA[4], SYYA[4];

real_type X0 = X[4];
real_type Y0 = Y[4];
real_type Z00 = Z[4][4];
real_type X0{ X[4] };
real_type Y0{ Y[4] };
real_type Z00{ Z[4][4] };

// Initial setting
real_type DXF[2] = {0,0};
real_type DXI[2] = {0,0};
real_type DYF[2] = {0,0};
real_type DYI[2] = {0,0};
real_type DXYF[2] = {0,0};
real_type DXYI[2] = {0,0};

for ( integer k = 0; k < 4; ++k ) {
int j1 = 4 + stencil[0][k];
int j2 = 4 + stencil[1][k];
int j3 = 4 + stencil[2][k];
real_type DXF[2]{0,0};
real_type DXI[2]{0,0};
real_type DYF[2]{0,0};
real_type DYI[2]{0,0};
real_type DXYF[2]{0,0};
real_type DXYI[2]{0,0};

for ( integer k{0}; k < 4; ++k ) {
integer j1{ 4 + stencil[0][k] };
integer j2{ 4 + stencil[1][k] };
integer j3{ 4 + stencil[2][k] };
if ( j1 >= int(jmin) && j3 <= int(jmax) )
estimate(
Z00, Z[4][j1], Z[4][j2], Z[4][j3],
Expand All @@ -171,19 +171,19 @@ namespace Splines {
);
}

for ( integer kx = 0; kx < 4; ++kx ) {
for ( integer kx{0}; kx < 4; ++kx ) {
int i1 = 4 + stencil[0][kx];
int i2 = 4 + stencil[1][kx];
int i3 = 4 + stencil[2][kx];

if ( i1 < int(imin) || i3 > int(imax) ) continue;

real_type X1 = X[i1] - X0;
real_type X2 = X[i2] - X0;
real_type X3 = X[i3] - X0;
real_type Z10 = Z[i1][4];
real_type Z20 = Z[i2][4];
real_type Z30 = Z[i3][4];
real_type X1{ X[i1] - X0 };
real_type X2{ X[i2] - X0 };
real_type X3{ X[i3] - X0 };
real_type Z10{ Z[i1][4] };
real_type Z20{ Z[i2][4] };
real_type Z30{ Z[i3][4] };
real_type CX1, CX2, CX3, SX, SXX, B00X, B10;
estimate(
Z00, Z10, Z20, Z30,
Expand All @@ -193,22 +193,22 @@ namespace Splines {
B00X, B10, DXF, DXI
);

for ( int ky = 0; ky < 4; ++ky ) {
int j1 = 4 + stencil[0][ky];
int j2 = 4 + stencil[1][ky];
int j3 = 4 + stencil[2][ky];
for ( int ky{0}; ky < 4; ++ky ) {
integer j1{ 4 + stencil[0][ky] };
integer j2{ 4 + stencil[1][ky] };
integer j3{ 4 + stencil[2][ky] };
if ( j1 < int(jmin) || j3 > int(jmax) ) continue;

real_type Y1 = Y[j1] - Y0;
real_type Y2 = Y[j2] - Y0;
real_type Y3 = Y[j3] - Y0;
real_type CY1 = CY1A[ky];
real_type CY2 = CY2A[ky];
real_type CY3 = CY3A[ky];
real_type SY = SYA[ky];
real_type SYY = SYYA[ky];
real_type B00Y = by0[ky];
real_type B01 = by1[ky];
real_type Y1 { Y[j1] - Y0 };
real_type Y2 { Y[j2] - Y0 };
real_type Y3 { Y[j3] - Y0 };
real_type CY1 { CY1A[ky] };
real_type CY2 { CY2A[ky] };
real_type CY3 { CY3A[ky] };
real_type SY { SYA[ky] };
real_type SYY { SYYA[ky] };
real_type B00Y{ by0[ky] };
real_type B01 { by1[ky] };

real_type Z01 = Z[4][j1], Z02 = Z[4][j2], Z03 = Z[4][j3];
real_type Z11 = Z[i1][j1], Z12 = Z[i1][j2], Z13 = Z[i1][j3];
Expand Down Expand Up @@ -295,7 +295,7 @@ namespace Splines {
void
Akima2Dspline::make_spline() {

size_t nn = size_t( m_nx*m_ny );
size_t nn{ size_t( m_nx*m_ny ) };
m_mem_bicubic.reallocate( 3*nn );
m_DX = m_mem_bicubic( nn );
m_DY = m_mem_bicubic( nn );
Expand All @@ -307,26 +307,26 @@ namespace Splines {

real_type x_loc[9], y_loc[9], z_loc[9][9];

for ( size_t i0 = 0; i0 < size_t(m_nx); ++i0 ) {
for ( size_t i0{0}; i0 < size_t(m_nx); ++i0 ) {
size_t imin = 4 > i0 ? 4-i0 : 0;
size_t imax = size_t(m_nx) < 5+i0 ? 3+(size_t(m_nx)-i0) : 8;

for ( size_t i = imin; i <= imax; ++i ) x_loc[i] = m_X[i+i0-4]-m_X[i0];
for ( size_t i{imin}; i <= imax; ++i ) x_loc[i] = m_X[i+i0-4]-m_X[i0];

for ( size_t j0 = 0; j0 < size_t(m_ny); ++j0 ) {
for ( size_t j0{0}; j0 < size_t(m_ny); ++j0 ) {
size_t jmin = 4 > j0 ? 4-j0 : 0;
size_t jmax = size_t(m_ny) < 5+j0 ? 3+(size_t(m_ny)-j0) : 8;

for ( size_t j = jmin; j <= jmax; ++j ) y_loc[j] = m_Y[j+j0-4]-m_Y[j0];
for ( size_t j{jmin}; j <= jmax; ++j ) y_loc[j] = m_Y[j+j0-4]-m_Y[j0];

for ( size_t i = imin; i <= imax; ++i )
for ( size_t j = jmin; j <= jmax; ++j )
for ( size_t i{imin}; i <= imax; ++i )
for ( size_t j{jmin}; j <= jmax; ++j )
z_loc[i][j] = m_Z[size_t(ipos_C(integer(i+i0-4),
integer(j+j0-4),
m_ny))];

// if not enough points, extrapolate
size_t iadd = 0, jadd = 0;
size_t iadd{0}, jadd{0};
if ( imax < 3+imin ) {
x_loc[imin-1] = 2*x_loc[imin] - x_loc[imax];
x_loc[imax+1] = 2*x_loc[imax] - x_loc[imin];
Expand All @@ -344,7 +344,7 @@ namespace Splines {
real_type x0 = x_loc[imin];
real_type x1 = x_loc[imin+1];
real_type x2 = x_loc[imax];
for ( size_t j = jmin; j <= jmax; ++j ) {
for ( size_t j{jmin}; j <= jmax; ++j ) {
real_type z0 = z_loc[imin][j];
real_type z1 = z_loc[imin+1][j];
real_type z2 = z_loc[imax][j];
Expand Down
32 changes: 16 additions & 16 deletions src/SplineBessel.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,19 +48,19 @@ namespace Splines {

void
Bessel_build(
real_type const * X,
real_type const * Y,
real_type * Yp,
integer npts
real_type const X[],
real_type const Y[],
real_type Yp[],
integer npts
) {

size_t n = size_t(npts > 0 ? npts-1 : 0);
size_t n{ size_t(npts > 0 ? npts-1 : 0) };

Malloc_real mem("Bessel_build");
real_type * m = mem.malloc( size_t(n+1) );
real_type * m{ mem.malloc( size_t(n+1) ) };

// calcolo slopes
for ( size_t i = 0; i < n; ++i )
for ( size_t i{0}; i < n; ++i )
m[i] = (Y[i+1]-Y[i])/(X[i+1]-X[i]);

if ( npts == 2 ) { // caso speciale 2 soli punti
Expand All @@ -69,7 +69,7 @@ namespace Splines {

} else {

for ( size_t i = 1; i < n; ++i ) {
for ( size_t i{1}; i < n; ++i ) {
real_type DL = X[i] - X[i-1];
real_type DR = X[i+1] - X[i];
Yp[i] = (DR*m[i-1]+DL*m[i])/((DL+DR));
Expand All @@ -86,16 +86,16 @@ namespace Splines {

void
BesselSpline::build() {
string msg = fmt::format("BesselSpline[{}]::build():", m_name );
string msg{ fmt::format("BesselSpline[{}]::build():", m_name ) };
UTILS_ASSERT(
m_npts > 1,
"{} npts = {} not enought points\n",
msg, m_npts
);
Utils::check_NaN( m_X, (msg+" X").c_str(), m_npts, __LINE__, __FILE__ );
Utils::check_NaN( m_Y, (msg+" Y").c_str(), m_npts, __LINE__, __FILE__ );
integer ibegin = 0;
integer iend = 0;
integer ibegin{0};
integer iend{0};
do {
// cerca intervallo monotono strettamente crescente
while ( ++iend < m_npts && m_X[iend-1] < m_X[iend] ) {}
Expand All @@ -118,17 +118,17 @@ namespace Splines {
// gc["ydata"]
//
*/
string where = fmt::format("BesselSpline[{}]::setup( gc ):", m_name );
GenericContainer const & gc_x = gc("xdata",where.c_str());
GenericContainer const & gc_y = gc("ydata",where.c_str());
string where{ fmt::format("BesselSpline[{}]::setup( gc ):", m_name ) };
GenericContainer const & gc_x{ gc("xdata",where.c_str()) };
GenericContainer const & gc_y{ gc("ydata",where.c_str()) };

vec_real_type x, y;
{
std::string ff = fmt::format( "{}, field `xdata'", where );
std::string ff{ fmt::format( "{}, field `xdata'", where ) };
gc_x.copyto_vec_real ( x, ff.c_str() );
}
{
std::string ff = fmt::format( "{}, field `ydata'", where );
std::string ff{ fmt::format( "{}, field `ydata'", where ) };
gc_y.copyto_vec_real ( y, ff.c_str() );
}
this->build( x, y );
Expand Down
14 changes: 7 additions & 7 deletions src/SplineBiCubic.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,22 +39,22 @@ namespace Splines {

void
BiCubicSpline::make_spline() {
size_t nn = size_t(m_nx*m_ny);
size_t nn{ size_t(m_nx*m_ny) };
m_mem_bicubic.reallocate( 3*nn );
m_DX = m_mem_bicubic( nn );
m_DY = m_mem_bicubic( nn );
m_DXY = m_mem_bicubic( nn );

// calcolo derivate
PchipSpline sp;
for ( integer j = 0; j < m_ny; ++j ) {
for ( integer j{0}; j < m_ny; ++j ) {
sp.build( m_X, 1, &m_Z[size_t(this->ipos_C(0,j))], m_ny, m_nx );
for ( integer i = 0; i < m_nx; ++i )
for ( integer i{0}; i < m_nx; ++i )
m_DX[size_t(this->ipos_C(i,j))] = sp.yp_node(i);
}
for ( integer i = 0; i < m_nx; ++i ) {
for ( integer i{0}; i < m_nx; ++i ) {
sp.build( m_Y, 1, &m_Z[size_t(this->ipos_C(i,0))], 1, m_ny );
for ( integer j = 0; j < m_ny; ++j )
for ( integer j{0}; j < m_ny; ++j )
m_DY[size_t(this->ipos_C(i,j))] = sp.yp_node(j);
}
std::fill_n( m_DXY, nn, 0 );
Expand All @@ -65,8 +65,8 @@ namespace Splines {
void
BiCubicSpline::write_to_stream( ostream_type & s ) const {
fmt::print( "Nx = {} Ny = {}\n", m_nx, m_ny );
for ( integer i = 1; i < m_nx; ++i ) {
for ( integer j = 1; j < m_ny; ++j ) {
for ( integer i{1}; i < m_nx; ++i ) {
for ( integer j{1}; j < m_ny; ++j ) {
size_t i00 = size_t(this->ipos_C(i-1,j-1));
size_t i10 = size_t(this->ipos_C(i,j-1));
size_t i01 = size_t(this->ipos_C(i-1,j));
Expand Down
Loading

0 comments on commit a2f4ed4

Please sign in to comment.