Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Initial review and refactoring of C++ code #2104

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .clang-tidy
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,7 @@ Checks: |
-readability-function-cognitive-complexity,
-cppcoreguidelines-avoid-non-const-global-variables,
-modernize-use-override,
-bugprone-easily-swappable-parameters,

WarningsAsErrors: '*'
HeaderFilterRegex: '*'
Expand Down
78 changes: 50 additions & 28 deletions include/CPnumerics.h
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
#include <vector>
#include <set>
#include <cfloat>
#include <stdlib.h> // For abs
#include <cstdlib> // For abs
#include <algorithm> // For max
#include <numeric>
#include <cmath>
Expand Down Expand Up @@ -53,8 +53,9 @@ class Spline
{
public:
/** An empty, invalid spline */
Spline() {}
Spline() = default;

// TODO: consider moving implementations to the .cpp file instead
/** A spline with x and y values */
Spline(const std::vector<X>& x, const std::vector<Y>& y) {
if (x.size() != y.size()) {
Expand All @@ -67,11 +68,17 @@ class Spline
return;
}

typedef typename std::vector<X>::difference_type size_type;
using size_type = typename std::vector<X>::difference_type;

size_type n = y.size() - 1;

std::vector<Y> b(n), d(n), a(n), c(n + 1), l(n + 1), u(n + 1), z(n + 1);
std::vector<Y> b(n);
std::vector<Y> d(n);
std::vector<Y> a(n);
std::vector<Y> c(n + 1);
std::vector<Y> l(n + 1);
std::vector<Y> u(n + 1);
std::vector<Y> z(n + 1);
std::vector<X> h(n + 1);

l[0] = Y(1);
Expand Down Expand Up @@ -100,14 +107,16 @@ class Spline
mElements.push_back(Element(x[i], y[i], b[i], c[i], d[i]));
}
}
virtual ~Spline() {}
virtual ~Spline() = default; // TODO why declare a virtual Dtor here?

Y operator[](const X& x) const {
return interpolate(x);
}

Y interpolate(const X& x) const {
if (mElements.size() == 0) return Y();
if (mElements.empty()) {
return Y();
}

typename std::vector<element_type>::const_iterator it;
it = std::lower_bound(mElements.begin(), mElements.end(), element_type(x));
Expand All @@ -120,7 +129,9 @@ class Spline

/* Evaluate at multiple locations, assuming xx is sorted ascending */
std::vector<Y> interpolate_vec(const std::vector<X>& xx) const {
if (mElements.size() == 0) return std::vector<Y>(xx.size());
if (mElements.empty()) {
return std::vector<Y>(xx.size());
}

typename std::vector<X>::const_iterator it;
typename std::vector<element_type>::const_iterator it2;
Expand Down Expand Up @@ -161,7 +172,7 @@ class Spline
Y a, b, c, d;
};

typedef Element element_type;
using element_type = Element;
std::vector<element_type> mElements;
};

Expand Down Expand Up @@ -192,7 +203,8 @@ std::vector<T> linspace(T xmin, T xmax, std::size_t n) {
template <typename T>
std::vector<T> log10space(T xmin, T xmax, std::size_t n) {
std::vector<T> x(n, 0.0);
T logxmin = log10(xmin), logxmax = log10(xmax);
T logxmin = log10(xmin);
T logxmax = log10(xmax);

for (std::size_t i = 0; i < n; ++i) {
x[i] = exp((logxmax - logxmin) / (n - 1) * i + logxmin);
Expand All @@ -203,7 +215,8 @@ std::vector<T> log10space(T xmin, T xmax, std::size_t n) {
template <typename T>
std::vector<T> logspace(T xmin, T xmax, std::size_t n) {
std::vector<T> x(n, 0.0);
T logxmin = log(xmin), logxmax = log(xmax);
T logxmin = log(xmin);
T logxmax = log(xmax);

for (std::size_t i = 0; i < n; ++i) {
x[i] = exp((logxmax - logxmin) / (n - 1) * i + logxmin);
Expand All @@ -221,8 +234,10 @@ std::vector<T> logspace(T xmin, T xmax, std::size_t n) {
*/
template <typename T>
void bisect_vector(const std::vector<T>& vec, T val, std::size_t& i) {
T rL, rM, rR;
std::size_t N = vec.size(), L = 0, R = N - 1, M = (L + R) / 2;
std::size_t N = vec.size();
std::size_t L = 0;
std::size_t R = N - 1;
std::size_t M = (L + R) / 2;
// Move the right limits in until they are good
while (!ValidNumber(vec[R])) {
if (R == 1) {
Expand All @@ -237,11 +252,13 @@ void bisect_vector(const std::vector<T>& vec, T val, std::size_t& i) {
}
L++;
}
rL = vec[L] - val;
rR = vec[R] - val;
T rL = vec[L] - val;
T rR = vec[R] - val;
T rM;
while (R - L > 1) {
if (!ValidNumber(vec[M])) {
std::size_t MR = M, ML = M;
std::size_t MR = M;
std::size_t ML = M;
// Move middle-right to the right until it is ok
while (!ValidNumber(vec[MR])) {
if (MR == vec.size() - 1) {
Expand Down Expand Up @@ -301,27 +318,31 @@ void bisect_vector(const std::vector<T>& vec, T val, std::size_t& i) {
*/
template <typename T>
void bisect_segmented_vector_slice(const std::vector<std::vector<T>>& mat, std::size_t j, T val, std::size_t& i) {
T rL, rM, rR;
std::size_t N = mat[j].size(), L = 0, R = N - 1, M = (L + R) / 2;
std::size_t N = mat[j].size();
std::size_t L = 0;
std::size_t R = N - 1;
std::size_t M = (L + R) / 2;
// Move the right limits in until they are good
while (!ValidNumber(mat[R][j])) {
if (R == 1) {
throw CoolProp::ValueError("All the values in bisection vector are invalid");
}
R--;
}
rR = mat[R][j] - val;
T rR = mat[R][j] - val;
// Move the left limits in until they are good
while (!ValidNumber(mat[L][j])) {
if (L == mat.size() - 1) {
throw CoolProp::ValueError("All the values in bisection vector are invalid");
}
L++;
}
rL = mat[L][j] - val;
T rL = mat[L][j] - val;
T rM;
while (R - L > 1) {
if (!ValidNumber(mat[M][j])) {
std::size_t MR = M, ML = M;
std::size_t MR = M;
std::size_t ML = M;
// Move middle-right to the right until it is ok
while (!ValidNumber(mat[MR][j])) {
if (MR == mat.size() - 1) {
Expand Down Expand Up @@ -443,17 +464,19 @@ class SplineClass
public:
double a, b, c, d;
SplineClass() : Nconstraints(0), A(4, std::vector<double>(4, 0)), B(4, 0), a(_HUGE), b(_HUGE), c(_HUGE), d(_HUGE) {}
bool build(void);
bool build();
bool add_value_constraint(double x, double y);
void add_4value_constraints(double x1, double x2, double x3, double x4, double y1, double y2, double y3, double y4);
bool add_derivative_constraint(double x, double dydx);
double evaluate(double x);
double evaluate(double x) const;
};

/// from http://stackoverflow.com/a/5721830/1360263
template <class T>
T factorial(T n) {
if (n == 0) return 1;
if (n == 0) {
return 1;
}
return n * factorial(n - 1);
}
/// see https://proofwiki.org/wiki/Nth_Derivative_of_Mth_Power
Expand Down Expand Up @@ -521,11 +544,10 @@ T CubicInterp(T x0, T x1, T x2, T x3, T f0, T f1, T f2, T f3, T x) {
Lagrange cubic interpolation as from
http://nd.edu/~jjwteach/441/PdfNotes/lecture6.pdf
*/
T L0, L1, L2, L3;
L0 = ((x - x1) * (x - x2) * (x - x3)) / ((x0 - x1) * (x0 - x2) * (x0 - x3));
L1 = ((x - x0) * (x - x2) * (x - x3)) / ((x1 - x0) * (x1 - x2) * (x1 - x3));
L2 = ((x - x0) * (x - x1) * (x - x3)) / ((x2 - x0) * (x2 - x1) * (x2 - x3));
L3 = ((x - x0) * (x - x1) * (x - x2)) / ((x3 - x0) * (x3 - x1) * (x3 - x2));
T L0 = ((x - x1) * (x - x2) * (x - x3)) / ((x0 - x1) * (x0 - x2) * (x0 - x3));
T L1 = ((x - x0) * (x - x2) * (x - x3)) / ((x1 - x0) * (x1 - x2) * (x1 - x3));
T L2 = ((x - x0) * (x - x1) * (x - x3)) / ((x2 - x0) * (x2 - x1) * (x2 - x3));
T L3 = ((x - x0) * (x - x1) * (x - x2)) / ((x3 - x0) * (x3 - x1) * (x3 - x2));
return L0 * f0 + L1 * f1 + L2 * f2 + L3 * f3;
};
/** /brief Calculate the first derivative of the function using a cubic interpolation form
Expand Down
8 changes: 5 additions & 3 deletions include/CoolPropLib.h
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,7 @@ EXPORT_CODE void CONVENTION Props1SImulti(const char* Outputs, char* backend, co
* \note If there is an error, a huge value will be returned, you can get the error message by doing something like get_global_param_string("errstring",output)
*/
EXPORT_CODE double CONVENTION PropsSI(const char* Output, const char* Name1, double Prop1, const char* Name2, double Prop2, const char* Ref);

/**
*\overload
*\sa \ref CoolProp::PropsSImulti(const std::vector<std::string>& Outputs, const std::string& Name1, const std::vector<double>& Prop1,
Expand All @@ -123,13 +124,14 @@ EXPORT_CODE double CONVENTION PropsSI(const char* Output, const char* Name1, dou
EXPORT_CODE void CONVENTION PropsSImulti(const char* Outputs, const char* Name1, double* Prop1, const long size_Prop1, const char* Name2,
double* Prop2, const long size_Prop2, char* backend, const char* FluidNames, const double* fractions,
const long length_fractions, double* result, long* resdim1, long* resdim2);

/**
*\overload
*\sa \ref CoolProp::PhaseSI(const std::string &, double, const std::string &, double, const std::string&)
*
* \note This function returns the phase string in pre-allocated phase variable. If buffer is not large enough, no copy is made
*/
EXPORT_CODE long CONVENTION PhaseSI(const char* Name1, double Prop1, const char* Name2, double Prop2, const char* Ref, char* phase, int n);
EXPORT_CODE long CONVENTION PhaseSI(const char* Name1, double Prop1, const char* Name2, double Prop2, const char* FluidName, char* phase, int n);

/**
*\overload
Expand Down Expand Up @@ -204,7 +206,7 @@ EXPORT_CODE int CONVENTION set_reference_stateD(const char* Ref, double T, doubl
* \note If there is an error, a huge value will be returned, you can get the error message by doing something like get_global_param_string("errstring",output)
*/
EXPORT_CODE void CONVENTION propssi_(const char* Output, const char* Name1, const double* Prop1, const char* Name2, const double* Prop2,
const char* Ref, double* output);
const char* FluidName, double* output);

/// Convert from degrees Fahrenheit to Kelvin (useful primarily for testing)
EXPORT_CODE double CONVENTION F2K(double T_F);
Expand All @@ -219,7 +221,7 @@ EXPORT_CODE long CONVENTION get_param_index(const char* param);
*
* @returns index The index as a long. If input is invalid, returns -1
*/
EXPORT_CODE long CONVENTION get_input_pair_index(const char* param);
EXPORT_CODE long CONVENTION get_input_pair_index(const char* pair);
/** \brief Redirect all output that would go to console (stdout) to a file
*/
EXPORT_CODE long CONVENTION redirect_stdout(const char* file);
Expand Down
9 changes: 5 additions & 4 deletions include/DataStructures.h
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ namespace CoolProp {

struct SimpleState
{
double rhomolar, T, p, hmolar, smolar, umolar, Q;
double rhomolar{}, T{}, p{}, hmolar{}, smolar{}, umolar{}, Q{};
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There seems to be some problem with the examples. I have not looked at this in detail, yet.
https://github.com/CoolProp/CoolProp/actions/runs/3735565861/jobs/6592680936#step:8:145

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this just forces default initialization, no?

SimpleState() {
fill(_HUGE);
}
Expand All @@ -28,7 +28,7 @@ struct SimpleState
umolar = v;
Q = v;
}
bool is_valid() {
bool is_valid() const {
return ValidNumber(rhomolar) && ValidNumber(T) && ValidNumber(hmolar) && ValidNumber(p);
}
};
Expand Down Expand Up @@ -61,6 +61,7 @@ struct SsatSimpleState : public SimpleState
/// The structure is taken directly from the AbstractState class.
//
// !! If you add a parameter, update the map in the corresponding CPP file !!
// TODO: make that an **enum class**
enum parameters
{
INVALID_PARAMETER = 0,
Expand Down Expand Up @@ -469,8 +470,8 @@ enum backends
};

/// Convert a string into the enum values
void extract_backend_families(std::string backend_string, backend_families& f1, backend_families& f2);
void extract_backend_families_string(std::string backend_string, backend_families& f1, std::string& f2);
void extract_backend_families(const std::string& backend_string, backend_families& f1, backend_families& f2);
void extract_backend_families_string(const std::string& backend_string, backend_families& f1, std::string& f2);
std::string get_backend_string(backends backend);

#if !defined(NO_FMTLIB) && FMT_VERSION >= 90000
Expand Down
Loading