<commit_before>#include <vector>
#include "Solvers.h"
#include "math.h"
#include "MatrixMath.h"
#include <iostream>
#include "CoolPropTools.h"

namespace CoolProp{

/**
In this formulation of the Multi-Dimensional Newton-Raphson solver the Jacobian matrix is known.
Therefore, the dx vector can be obtained from 

J(x)dx=-f(x)

for a given value of x.  The pointer to the class FuncWrapperND that is passed in must implement the call() and Jacobian() 
functions, each of which take the vector x. The data is managed using std::vector<double> vectors

@param f A pointer to an subclass of the FuncWrapperND class that implements the call() and Jacobian() functions
@param x0 The initial guess value for the solution
@param tol The root-sum-square of the errors from each of the components
@param maxiter The maximum number of iterations
@param errstring  A string with the returned error.  If the length of errstring is zero, no errors were found
@returns If no errors are found, the solution.  Otherwise, _HUGE, the value for infinity
*/
std::vector<double> NDNewtonRaphson_Jacobian(FuncWrapperND *f, std::vector<double> x0, double tol, int maxiter, std::string *errstring)
{
	int iter=0;
	*errstring=std::string("");
	std::vector<double> f0,v,negative_f0;
	std::vector<std::vector<double> > J;
	double error = 999;
	while (iter==0 || fabs(error)>tol){
		f0 = f->call(x0);
		J = f->Jacobian(x0);
		
		// Negate f0
		negative_f0 = f0;
		for (unsigned int i = 0; i<f0.size(); i++){ negative_f0[i] *= -1;}
		// find v from J*v = -f
		v = linsolve(J, negative_f0);
		// Update the guess
		x0[0] += v[0];
		x0[1] += v[1];
		error = root_sum_square(f0);
		if (iter>maxiter){
			*errstring=std::string("reached maximum number of iterations");
			x0[0]=_HUGE;
		}
		iter++;
	}
	return x0;
}

/**
In the newton function, a 1-D Newton-Raphson solver is implemented using exact solutions.  An initial guess for the solution is provided.

@param f A pointer to an instance of the FuncWrapper1D class that implements the call() function 
@param x0 The inital guess for the solution
@param ftol The absolute value of the tolerance accepted for the objective function
@param maxiter Maximum number of iterations
@param errstring A pointer to the std::string that returns the error from Secant.  Length is zero if no errors are found
@returns If no errors are found, the solution, otherwise the value _HUGE, the value for infinity
*/
double Newton(FuncWrapper1D &f, double x0, double ftol, int maxiter, std::string &errstring)
{
	double x, dx, fval=999;
    int iter=1;
    errstring.clear();
	x = x0;
    while (iter < 2 || fabs(fval) > ftol)
    {
		fval = f.call(x);
		dx = -fval/f.deriv(x);
		x += dx;
		
		if (fabs(dx/x) < 10*DBL_EPSILON)
		{
			return x;
		}

		if (iter>maxiter)
		{
			errstring= "reached maximum number of iterations";
			throw SolutionError(format("Newton reached maximum number of iterations"));
		}
        iter=iter+1;
    }
    return x;
}

/**
In the secant function, a 1-D Newton-Raphson solver is implemented.  An initial guess for the solution is provided.

@param f A pointer to an instance of the FuncWrapper1D class that implements the call() function 
@param x0 The inital guess for the solutionh
@param dx The initial amount that is added to x in order to build the numerical derivative
@param tol The absolute value of the tolerance accepted for the objective function
@param maxiter Maximum number of iterations
@param errstring A pointer to the std::string that returns the error from Secant.  Length is zero if no errors are found
@returns If no errors are found, the solution, otherwise the value _HUGE, the value for infinity
*/
double Secant(FuncWrapper1D &f, double x0, double dx, double tol, int maxiter, std::string &errstring)
{
	double x1=0,x2=0,x3=0,y1=0,y2=0,x,fval=999;
    int iter=1;
	errstring = "";
	
	if (fabs(dx)==0){ errstring="dx cannot be zero"; return _HUGE;}
    while (iter<=2 || fabs(fval)>tol)
    {
        if (iter==1){x1=x0; x=x1;}
        if (iter==2){x2=x0+dx; x=x2;}
        if (iter>2) {x=x2;}
			fval=f.call(x);
            if (!ValidNumber(fval)){throw ValueError("Residual function in secant returned invalid number");};
        if (iter==1){y1=fval;}<commit_msg>Additional error checking and debugging in Solvers<commit_after>#include <vector>
#include "Solvers.h"
#include "math.h"
#include "MatrixMath.h"
#include <iostream>
#include "CoolPropTools.h"

namespace CoolProp{

/**
In this formulation of the Multi-Dimensional Newton-Raphson solver the Jacobian matrix is known.
Therefore, the dx vector can be obtained from 

J(x)dx=-f(x)

for a given value of x.  The pointer to the class FuncWrapperND that is passed in must implement the call() and Jacobian() 
functions, each of which take the vector x. The data is managed using std::vector<double> vectors

@param f A pointer to an subclass of the FuncWrapperND class that implements the call() and Jacobian() functions
@param x0 The initial guess value for the solution
@param tol The root-sum-square of the errors from each of the components
@param maxiter The maximum number of iterations
@param errstring  A string with the returned error.  If the length of errstring is zero, no errors were found
@returns If no errors are found, the solution.  Otherwise, _HUGE, the value for infinity
*/
std::vector<double> NDNewtonRaphson_Jacobian(FuncWrapperND *f, std::vector<double> x0, double tol, int maxiter, std::string *errstring)
{
	int iter=0;
	*errstring=std::string("");
	std::vector<double> f0,v,negative_f0;
	std::vector<std::vector<double> > J;
	double error = 999;
	while (iter==0 || fabs(error)>tol){
		f0 = f->call(x0);
		J = f->Jacobian(x0);
		
		// Negate f0
		negative_f0 = f0;
		for (unsigned int i = 0; i<f0.size(); i++){ negative_f0[i] *= -1;}
		// find v from J*v = -f
		v = linsolve(J, negative_f0);
		// Update the guess
		x0[0] += v[0];
		x0[1] += v[1];
		error = root_sum_square(f0);
		if (iter>maxiter){
			*errstring=std::string("reached maximum number of iterations");
			x0[0]=_HUGE;
		}
		iter++;
	}
	return x0;
}

/**
In the newton function, a 1-D Newton-Raphson solver is implemented using exact solutions.  An initial guess for the solution is provided.

@param f A pointer to an instance of the FuncWrapper1D class that implements the call() function 
@param x0 The inital guess for the solution
@param ftol The absolute value of the tolerance accepted for the objective function
@param maxiter Maximum number of iterations
@param errstring A pointer to the std::string that returns the error from Secant.  Length is zero if no errors are found
@returns If no errors are found, the solution, otherwise the value _HUGE, the value for infinity
*/
double Newton(FuncWrapper1D &f, double x0, double ftol, int maxiter, std::string &errstring)
{
	double x, dx, fval=999;
    int iter=1;
    errstring.clear();
	x = x0;
    while (iter < 2 || fabs(fval) > ftol)
    {
		fval = f.call(x);
		dx = -fval/f.deriv(x);

		if (!ValidNumber(fval)){
            throw ValueError("Residual function in newton returned invalid number");
        };

        x += dx;

		if (fabs(dx/x) < 10*DBL_EPSILON)
		{
			return x;
		}

		if (iter>maxiter)
		{
			errstring= "reached maximum number of iterations";
			throw SolutionError(format("Newton reached maximum number of iterations"));
		}
        iter=iter+1;
    }
    return x;
}

/**
In the secant function, a 1-D Newton-Raphson solver is implemented.  An initial guess for the solution is provided.

@param f A pointer to an instance of the FuncWrapper1D class that implements the call() function 
@param x0 The inital guess for the solutionh
@param dx The initial amount that is added to x in order to build the numerical derivative
@param tol The absolute value of the tolerance accepted for the objective function
@param maxiter Maximum number of iterations
@param errstring A pointer to the std::string that returns the error from Secant.  Length is zero if no errors are found
@returns If no errors are found, the solution, otherwise the value _HUGE, the value for infinity
*/
double Secant(FuncWrapper1D &f, double x0, double dx, double tol, int maxiter, std::string &errstring)
{
    #if defined(COOLPROP_DEEP_DEBUG)
    static std::vector<double> xlog, flog;
    xlog.clear(); flog.clear();
    #endif

	double x1=0,x2=0,x3=0,y1=0,y2=0,x,fval=999;
    int iter=1;
	errstring = "";
	
	if (fabs(dx)==0){ errstring="dx cannot be zero"; return _HUGE;}
    while (iter<=2 || fabs(fval)>tol)
    {
        if (iter==1){x1=x0; x=x1;}
        if (iter==2){x2=x0+dx; x=x2;}
        if (iter>2) {x=x2;}

            fval = f.call(x);

            #if defined(COOLPROP_DEEP_DEBUG)
                xlog.push_back(x);
                flog.push_back(fval);
            #endif

            if (!ValidNumber(fval)){
                throw ValueError("Residual function in secant returned invalid number");
            };
        if (iter==1){y1=fval;}
{'fraction_non_alphanumeric': 0.1063782087209954, 'fraction_numerical': 0.01726263871763255, 'mean_word_length': 3.846279196089082, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 36, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4961947', 'n_tokens_mistral': 3048, 'n_tokens_neox': 2869, 'n_words': 1115}