Updating miniapp (efc4d348) · Commits · Berrill, Mark A / RayTrace-miniapp

src/AtomicModel/interp.cpp

+202 −73

Original line number	Diff line number	Diff line
		@@ -21,7 +21,7 @@ double interp::interp_linear( size_t N, const double xi, const double yi, doub
		double dx = ( x - xi[i - 1] ) / ( xi[i] - xi[i - 1] );
		double dx2 = 1.0 - dx;
		double y = dx2 * yi[i - 1] + dx * yi[i];
		return ( y );
		return y;
		}


		@@ -206,98 +206,227 @@ double interp::calc_width( size_t n, const double x, const double y )
		}


		// Calculate coefficients for the bisection method
		double interp::bisection_coeff(
		size_t N, const double x_in, const double r_in, std::array<double, 2> &range )
		// Get the 2 points on either side of a sign change
		static inline void getPoints( size_t N, const double x_in, const double r_in,
		std::array<double, 4> &x, std::array<double, 4> &f, std::array<double, 2> &range )
		{
		if ( N < 2 )
		throw std::logic_error( "Error: N<2" );
		if ( N == 2 ) {
		range[0] = x_in[0];
		range[1] = x_in[1];
		if ( range[0] > range[1] )
		std::swap( range[0], range[1] );
		return 0.5 * ( x_in[0] + x_in[1] );
		}
		// Get the sign of the function such that f(-inf) is negitive
		double sign = 1;
		// Get the sign of f(-inf)
		double x_min = 1e200;
		bool sign = false;
		for ( size_t i = 0; i < N; i++ ) {
		if ( x_in[i] < x_min ) {
		x_min = x_in[i];
		sign = r_in[i] < 0 ? 1 : -1;
		sign = r_in[i] >= 0;
		}
		}
		// Copy the closest two points on each side
		double x[4] = { -1e200, -1e200, 1e200, 1e200 };
		double f[4] = { -1e200, -1e200, 1e200, 1e200 };
		// Get the closest two points on either side (if they exist)
		x = { -1e210, -1e210, 1e210, 1e210 };
		f = { 1e210, 1e210, 1e210, 1e210 };
		for ( size_t i = 0; i < N; i++ ) {
		double x0 = x_in[i];
		double f0 = r_in[i] * sign;
		if ( f0 < 0.0 ) {
		if ( x0 > x[0] ) {
		x[0] = x0;
		f[0] = f0;
		}
		if ( ( r_in[i] >= 0 ) == sign ) {
		if ( x_in[i] > x[0] ) {
		x[0] = x_in[i];
		f[0] = r_in[i];
		if ( x[0] > x[1] ) {
		std::swap( x[0], x[1] );
		std::swap( f[0], f[1] );
		}
		} else {
		if ( x0 < x[3] ) {
		x[3] = x0;
		f[3] = f0;
		}
		} else {
		if ( x_in[i] < x[3] ) {
		x[3] = x_in[i];
		f[3] = r_in[i];
		if ( x[2] > x[3] ) {
		std::swap( x[2], x[3] );
		std::swap( f[2], f[3] );
		}
		}
		}
		if ( x[1] > x[2] ) {
		// We have multiple zeros, try to keep one
		if ( x[0] < x[2] ) {
		range[0] = x[0];
		range[1] = x[2];
		return 0.5 * ( x[0] + x[2] );
		} else if ( x[1] < x[3] ) {
		range[0] = x[1];
		range[1] = x[3];
		return 0.5 * ( x[1] + x[3] );
		} else {
		throw std::logic_error( "unable to find valid range" );
		}
		// Fix the sign and order of x such that f(-inf) is negitive
		if ( f[1] > 0.0 ) {
		f[0] = -f[0];
		f[1] = -f[1];
		f[2] = -f[2];
		f[3] = -f[3];
		}
		range[0] = x[1];
		range[1] = x[2];
		}


		// Calculate coefficients for the bisection method
		double interp::bisection_coeff(
		size_t N, const double x_in, const double r_in, std::array<double, 2> &range )
		{
		if ( N < 2 )
		throw std::logic_error( "Error: N<2" );
		if ( N == 2 ) {
		range[0] = x_in[0];
		range[1] = x_in[1];
		if ( range[0] > range[1] )
		std::swap( range[0], range[1] );
		return 0.5 * ( x_in[0] + x_in[1] );
		}
		// Get the points on each side of the sign change
		std::array<double, 4> x, f;
		getPoints( N, x_in, r_in, x, f, range );
		// If we have taken too many iterations, default to basic bisection
		if ( N > 30 ) {
		// We are having difficulties, switch back to simple bisection
		return 0.5 * ( x[1] + x[2] );
		}
		if ( fabs( f[0] ) > 1e100 ) {
		// Use an uneven bisection to try an ensure we have 2 points on each side
		double r = exp( -( N - 1. ) );
		return ( 1.0 - r ) * x[1] + r * x[2];
		} else if ( fabs( f[3] ) > 1e100 ) {
		// Use an uneven bisection to try an ensure we have 2 points on each side
		double r = exp( -( N - 1. ) );
		return ( 1.0 - r ) * x[2] + r * x[1];
		}
		// Use the 4 closest points to extimate the zero-crossing
		if ( f[0] > f[1] ) {
		std::swap( x[0], x[1] );
		std::swap( f[0], f[1] );
		} else if ( f[2] > f[3] ) {
		std::swap( x[2], x[3] );
		std::swap( f[2], f[3] );
		return 0.5 * ( range[0] + range[1] );
		}
		// Try to get 2 points on each side
		if ( fabs( f[0] ) > 1e200 \|\| fabs( f[3] ) > 1e200 ) {
		if ( fabs( f[0] ) > 1e200 )
		return 0.65 * x[1] + 0.35 * x[2];
		else
		return 0.35 * x[1] + 0.65 * x[2];
		}
		// Use the 4 closest points to estimate the zero-crossing
		auto xi = x;
		auto fi = f;
		if ( fi[0] > fi[1] ) {
		std::swap( xi[0], xi[1] );
		std::swap( fi[0], fi[1] );
		}
		if ( fi[2] > fi[3] ) {
		std::swap( xi[2], xi[3] );
		std::swap( fi[2], fi[3] );
		}
		double y0 = pchip( 4, 2, fi.data(), xi.data(), 0.0 );
		// Use the points on either side to get additional estimates
		double y1 = x[0] - f[0] * ( x[1] - x[0] ) / ( f[1] - f[0] );
		double y2 = x[2] - f[2] * ( x[3] - x[2] ) / ( f[3] - f[2] );
		// Compute the standard deviation of the guesses
		double ym = 0.333333333333 * ( y0 + y1 + y2 );
		double s2 = ( y0 - ym ) * ( y0 - ym ) + ( y1 - ym ) * ( y1 - ym ) + ( y1 - ym ) * ( y1 - ym );
		double tol = 0.1 * ( range[1] - range[0] ) * ( range[1] - range[0] );
		// Choose the optimal value for the interpolation
		double y, tol2 = 0.1;
		if ( s2 < tol ) {
		// We are well approximated by a polynomial
		y = y0;
		tol2 = 0.02;
		} else if ( y1 > x[1] && y2 > x[1] && y1 < x[2] && y2 < x[2] ) {
		// Use the average of the linear approximations
		y = 0.5 * ( y1 + y2 );
		} else {
		// Default to regular bisection
		y = 0.5 * ( x[1] + x[2] );
		}
		double y = pchip( 4, 2, f, x, 0.0 );
		// Check that y is within the range
		y = std::min( y, x[2] );
		y = std::max( y, x[1] );
		// Adjust y to move it toward the center
		y = 0.95 * y + 0.05 * ( 0.5 * ( x[1] + x[2] ) );
		// Make sure y is in the range, adjusting if necessary
		y = std::max( y, ( 1 - tol2 ) * range[0] + tol2 * range[1] );
		y = std::min( y, tol2 * range[0] + ( 1 - tol2 ) * range[1] );
		if ( y != y )
		y = 0.5 * ( x[1] + x[2] );
		y = 0.5 * ( range[0] + range[1] );
		return y;
		}


		// Modified bisection root finding
		double interp::bisection(
		std::function<double( double )> fun, double lb, double ub, double tol1, double tol2, int *N )
		{
		if ( ub <= lb )
		throw std::logic_error( "Error: lb <= ub" );
		std::array<double, 2> range = { lb, ub };
		double x[101] = { 0 }, f[101] = { 0 };
		x[0] = lb;
		x[1] = ub;
		f[0] = fun( x[0] );
		f[1] = fun( x[1] );
		if ( fabs( f[0] ) < tol1 \|\| fabs( f[1] ) < tol1 )
		return fabs( f[0] ) < tol1 ? x[0] : x[1];
		if ( ( f[0] < 0 && f[1] < 0 ) \|\| ( f[0] > 0 && f[1] > 0 ) )
		throw std::logic_error( "Error: sign(f(lb)) == sign(f(ub))" );
		size_t i = 2;
		while ( true ) {
		// Get the next guess
		x[i] = bisection_coeff( i, x, f, range );
		if ( ( range[1] - range[0] ) < tol2 )
		break;
		// Evaluate f
		f[i] = fun( x[i] );
		if ( fabs( f[i] ) < tol1 )
		break;
		if ( i >= 100 )
		throw std::logic_error( "Excessive number of iterations" );
		i++;
		}
		if ( N )
		*N = i + 1;
		return x[i];
		}


		// Modified fixed_point solve
		double interp::fixed_point( std::function<double( double )> fun, double x0, double lb, double ub,
		double tol1, double tol2, int *N )
		{
		if ( ub <= lb )
		throw std::logic_error( "Error: lb <= ub" );
		std::array<double, 2> range = { lb, ub };
		double x[256] = { 0 }, f[256] = { 0 }, r[256] = { 0 };
		x[0] = x0;
		f[0] = fun( x[0] );
		r[0] = f[0] - x[0];
		if ( fabs( r[0] ) < tol1 ) {
		if ( N )
		*N = 1;
		return x[0];
		}
		int i = 1;
		bool set_lb = false;
		int N1 = 1;
		int N2 = 0;
		while ( true ) {
		// Get the next guess
		if ( ( std::min( N1, N2 ) > 0 && i >= 5 ) \|\| std::min( N1, N2 ) >= 2 ) {
		// Perform bisection
		x[i] = bisection_coeff( i, x, r, range );
		} else if ( std::min( N1, N2 ) > 0 && i >= 5 ) {
		// Try to bound the solution
		if ( !set_lb ) {
		set_lb = true;
		x[i] = lb;
		} else {
		x[i] = ub;
		}
		} else {
		// Perform fixed point iteration
		x[i] = f[i - 1];
		if ( x[i] <= range[0] \|\| x[i] >= range[1] ) {
		if ( std::min( N1, N2 ) > 0 ) {
		x[i] = bisection_coeff( i, x, r, range );
		} else if ( !set_lb ) {
		x[i] = lb;
		set_lb = true;
		} else {
		x[i] = ub;
		}
		}
		}
		if ( ( range[1] - range[0] ) < tol2 )
		break;
		// Evaluate f
		f[i] = fun( x[i] );
		r[i] = f[i] - x[i];
		if ( fabs( r[i] ) < tol1 )
		break;
		// Check if we have bounded the solution
		if ( ( r[0] >= 0 ) == ( r[i] >= 0 ) )
		N1++;
		else
		N2++;
		if ( std::min( N1, N2 ) == 1 ) {
		std::array<double, 4> xi, fi;
		getPoints( i + 1, x, r, xi, fi, range );
		}
		if ( i >= 250 )
		throw std::logic_error( "Excessive number of iterations" );
		i++;
		}
		if ( N )
		*N = i + 1;
		return x[i];
		}

src/AtomicModel/interp.h

+34 −2

Original line number	Diff line number	Diff line
		@@ -368,6 +368,21 @@ inline void unique( size_t nx, const double X, size_t ny, double Y, size_t I
		double calc_width( size_t n, const double x, const double y );


		/*!
		* @brief Function to solve for the zero of f(x)
		* @details This function solves for the solution of f(x)=0 using a modified bisection method.
		* @param fun The function to solve of the form f(x,options), where
		* options is a pointer to any additional data needed by the function
		* @param lb The lower bound to search
		* @param ub The upper bound to search
		* @param tol1 The tolerance of the function evaluation to stop: abs(f(x)) < tol1
		* @param tol2 The tolerance of x to stop: abs(x1-x2) < tol2
		* @param[out] N Optional argument to return the number function evaluations used
		*/
		double bisection( std::function<double( double )> fun, double lb, double ub, double tol1,
		double tol2, int *N = nullptr );


		/*!
		* @brief Function to solve for the zero of f(x)
		* @details This function solves for the solution of f(x)=0 using a modified bisection method.
		@@ -380,8 +395,25 @@ double calc_width( size_t n, const double x, const double y );
		* @param options Any additional options that are needed by the function
		*/
		template<class... Args>
		double bisection( std::function<double( double, Args... )> fun, double lb, double ub, double tol1,
		double tol2, Args... options );
		double bisection( const std::function<double( double, Args... )> &fun, double lb, double ub,
		double tol1, double tol2, Args... options );


		/*!
		* @brief Function to solve x = f(x)
		* @details This function solves x = f(x) using a bounded fixed point method or bisection
		* depending on the convergence of f(x)
		* @param fun The function to solve of the form f(x,options), where
		* options is a pointer to any additional data needed by the function
		* @param[in] x0 The initial guess
		* @param[in] lb The lower bound to search
		* @param[in] ub The upper bound to search
		* @param[in] tol1 The tolerance of the function evaluation to stop: abs(f(x)) < tol1
		* @param[in] tol2 The tolerance of x to stop: abs(x1-x2) < tol2
		* @param[out] N Optional argument to return the number function evaluations used
		*/
		double fixed_point( std::function<double( double )> fun, double x0, double lb, double ub,
		double tol1, double tol2, int *N = nullptr );


		/*!

src/AtomicModel/interp.hpp

+620 −710

Original line number	Diff line number	Diff line
		@@ -165,79 +165,55 @@ HOST_DEVICE inline size_t interp::findfirstsingle( const TYPE *X, size_t size_X,

		// Subroutine to perform a quicksort
		template<class T>
		void interp::quicksort( size_t a_size, T *arr )
		void interp::quicksort( size_t n, T *arr )
		{
		long int n = static_cast<long int>( a_size );
		if ( n <= 1 )
		return;
		bool test;
		long int i, ir, j, jstack, k, l, istack[100];
		T a, tmp_a;
		jstack = 0;
		l = 0;
		ir = n - 1;
		int64_t jstack = 0;
		int64_t l = 0;
		int64_t ir = n - 1;
		int64_t istack[100];
		while ( 1 ) {
		if ( ir - l < 7 ) { // Insertion sort when subarray small enough.
		for ( j = l + 1; j <= ir; j++ ) {
		a = arr[j];
		test = true;
		for ( i = j - 1; i >= 0; i-- ) {
		if ( arr[i] < a ) {
		arr[i + 1] = a;
		test = false;
		if ( ir - l < 7 ) { // Insertion sort when subarray small enough
		for ( int64_t j = l + 1; j <= ir; j++ ) {
		for ( int64_t i = j - 1; i >= l; i-- ) {
		if ( arr[i] < arr[i + 1] )
		break;
		}
		arr[i + 1] = arr[i];
		}
		if ( test ) {
		i = l - 1;
		arr[i + 1] = a;
		std::swap( arr[i], arr[i + 1] );
		}
		}
		if ( jstack == 0 )
		return;
		ir = istack[jstack]; // Pop stack and begin a new round of partitioning.
		ir = istack[jstack]; // Pop stack and begin a new round of partitioning
		l = istack[jstack - 1];
		jstack -= 2;
		} else {
		k = ( l + ir ) / 2; // Choose median of left, center and right elements as partitioning
		// element a. Also rearrange so that a(l) < a(l+1) < a(ir).
		tmp_a = arr[k];
		arr[k] = arr[l + 1];
		arr[l + 1] = tmp_a;
		if ( arr[l] > arr[ir] ) {
		tmp_a = arr[l];
		arr[l] = arr[ir];
		arr[ir] = tmp_a;
		}
		if ( arr[l + 1] > arr[ir] ) {
		tmp_a = arr[l + 1];
		arr[l + 1] = arr[ir];
		arr[ir] = tmp_a;
		}
		if ( arr[l] > arr[l + 1] ) {
		tmp_a = arr[l];
		arr[l] = arr[l + 1];
		arr[l + 1] = tmp_a;
		}
		// Scan up to find element > a
		j = ir;
		a = arr[l + 1]; // Partitioning element.
		// Choose median of left, center and right elements as partitioning element a
		// Also rearrange so that a(l) ? a(l+1) ? a(ir)
		int64_t k = ( l + ir ) / 2;
		std::swap( arr[k], arr[l + 1] );
		if ( arr[l] > arr[ir] )
		std::swap( arr[l], arr[ir] );
		if ( arr[l + 1] > arr[ir] )
		std::swap( arr[l + 1], arr[ir] );
		if ( arr[l] > arr[l + 1] )
		std::swap( arr[l], arr[l + 1] );
		// Partitioning element
		int64_t j = ir;
		int64_t i;
		for ( i = l + 2; i <= ir; i++ ) {
		if ( arr[i] < a )
		if ( arr[i] < arr[l + 1] )
		continue;
		while ( arr[j] > a ) // Scan down to find element < a.
		while ( arr[j] > arr[l + 1] ) // Scan down to find element < a
		j--;
		if ( j < i )
		break; // Pointers crossed. Exit with partitioning complete.
		tmp_a = arr[i]; // Exchange elements of both arrays.
		arr[i] = arr[j];
		arr[j] = tmp_a;
		break; // Pointers crossed. Exit with partitioning complete
		std::swap( arr[i], arr[j] ); // Exchange elements of both arrays
		}
		arr[l + 1] = arr[j]; // Insert partitioning element in both arrays.
		arr[j] = a;
		// Insert partitioning element in both arrays
		std::swap( arr[l + 1], arr[j] );
		jstack += 2;
		// Push pointers to larger subarray on stack, process smaller subarray immediately.
		// Push pointers to larger subarray on stack, process smaller subarray immediately
		if ( ir - i + 1 >= j - l ) {
		istack[jstack] = ir;
		istack[jstack - 1] = i;
		@@ -250,106 +226,74 @@ void interp::quicksort( size_t a_size, T *arr )
		}
		}
		}
		template<class type>
		inline void interp::quicksort( std::vector<type> &X )
		{
		quicksort( X.size(), X.data() );
		}


		// Subroutine to perform a quicksort
		template<class T1, class T2>
		void interp::quicksort( size_t size, T1 arr, T2 brr )
		void interp::quicksort( size_t n, T1 arr, T2 brr )
		{
		long int n = static_cast<long int>( size );
		if ( n <= 1 )
		return;
		bool test;
		long int i, ir, j, jstack, k, l, istack[100];
		T1 a, tmp_a;
		T2 b, tmp_b;
		jstack = 0;
		l = 0;
		ir = n - 1;
		int64_t jstack = 0;
		int64_t l = 0;
		int64_t ir = n - 1;
		int64_t istack[100];
		while ( 1 ) {
		if ( ir - l < 7 ) { // Insertion sort when subarray small enough.
		for ( j = l + 1; j <= ir; j++ ) {
		a = arr[j];
		b = brr[j];
		test = true;
		for ( i = j - 1; i >= 0; i-- ) {
		if ( arr[i] < a ) {
		arr[i + 1] = a;
		brr[i + 1] = b;
		test = false;
		if ( ir - l < 7 ) { // Insertion sort when subarray small enough
		for ( int64_t j = l + 1; j <= ir; j++ ) {
		for ( int64_t i = j - 1; i >= l; i-- ) {
		if ( arr[i] < arr[i + 1] )
		break;
		}
		arr[i + 1] = arr[i];
		brr[i + 1] = brr[i];
		}
		if ( test ) {
		i = l - 1;
		arr[i + 1] = a;
		brr[i + 1] = b;
		std::swap( arr[i], arr[i + 1] );
		std::swap( brr[i], brr[i + 1] );
		}
		}
		if ( jstack == 0 )
		return;
		ir = istack[jstack]; // Pop stack and begin a new round of partitioning.
		ir = istack[jstack]; // Pop stack and begin a new round of partitioning
		l = istack[jstack - 1];
		jstack -= 2;
		} else {
		k = ( l + ir ) / 2; // Choose median of left, center and right elements as partitioning
		// element a. Also rearrange so that a(l) ? a(l+1) ? a(ir).
		tmp_a = arr[k];
		arr[k] = arr[l + 1];
		arr[l + 1] = tmp_a;
		tmp_b = brr[k];
		brr[k] = brr[l + 1];
		brr[l + 1] = tmp_b;
		// Choose median of left, center and right elements as partitioning element a
		// Also rearrange so that a(l) ? a(l+1) ? a(ir)
		int64_t k = ( l + ir ) / 2;
		std::swap( arr[k], arr[l + 1] );
		std::swap( brr[k], brr[l + 1] );
		if ( arr[l] > arr[ir] ) {
		tmp_a = arr[l];
		arr[l] = arr[ir];
		arr[ir] = tmp_a;
		tmp_b = brr[l];
		brr[l] = brr[ir];
		brr[ir] = tmp_b;
		std::swap( arr[l], arr[ir] );
		std::swap( brr[l], brr[ir] );
		}
		if ( arr[l + 1] > arr[ir] ) {
		tmp_a = arr[l + 1];
		arr[l + 1] = arr[ir];
		arr[ir] = tmp_a;
		tmp_b = brr[l + 1];
		brr[l + 1] = brr[ir];
		brr[ir] = tmp_b;
		std::swap( arr[l + 1], arr[ir] );
		std::swap( brr[l + 1], brr[ir] );
		}
		if ( arr[l] > arr[l + 1] ) {
		tmp_a = arr[l];
		arr[l] = arr[l + 1];
		arr[l + 1] = tmp_a;
		tmp_b = brr[l];
		brr[l] = brr[l + 1];
		brr[l + 1] = tmp_b;
		}
		// Scan up to find element > a
		j = ir;
		a = arr[l + 1]; // Partitioning element.
		b = brr[l + 1];
		std::swap( arr[l], arr[l + 1] );
		std::swap( brr[l], brr[l + 1] );
		}
		// Partitioning element
		int64_t j = ir;
		int64_t i;
		for ( i = l + 2; i <= ir; i++ ) {
		if ( arr[i] < a )
		if ( arr[i] < arr[l + 1] )
		continue;
		while ( arr[j] > a ) // Scan down to find element < a.
		while ( arr[j] > arr[l + 1] ) // Scan down to find element < a
		j--;
		if ( j < i )
		break; // Pointers crossed. Exit with partitioning complete.
		tmp_a = arr[i]; // Exchange elements of both arrays.
		arr[i] = arr[j];
		arr[j] = tmp_a;
		tmp_b = brr[i];
		brr[i] = brr[j];
		brr[j] = tmp_b;
		}
		arr[l + 1] = arr[j]; // Insert partitioning element in both arrays.
		arr[j] = a;
		brr[l + 1] = brr[j];
		brr[j] = b;
		break; // Pointers crossed. Exit with partitioning complete
		std::swap( arr[i], arr[j] ); // Exchange elements of both arrays
		std::swap( brr[i], brr[j] );
		}
		// Insert partitioning element in both arrays
		std::swap( arr[l + 1], arr[j] );
		std::swap( brr[l + 1], brr[j] );
		jstack += 2;
		// Push pointers to larger subarray on stack, process smaller subarray immediately.
		// Push pointers to larger subarray on stack, process smaller subarray immediately
		if ( ir - i + 1 >= j - l ) {
		istack[jstack] = ir;
		istack[jstack - 1] = i;
		@@ -362,17 +306,8 @@ void interp::quicksort( size_t size, T1 arr, T2 brr )
		}
		}
		}


		template<class type>
		inline void quicksort( std::vector<type> &X )
		{
		quicksort( X.size(), X.data() );
		}


		template<class type1, class type2>
		inline void quicksort( std::vector<type1> &X, std::vector<type2> &Y )
		inline void interp::quicksort( std::vector<type1> &X, std::vector<type2> &Y )
		{
		if ( X.size() != Y.size() )
		throw std::logic_error( "Error: X.size() != Y>size()" );
		@@ -460,36 +395,11 @@ inline void interp::unique(

		// Modified bisection root finding
		template<class... Args>
		double interp::bisection( std::function<double( double, Args... )> fun, double lb, double ub,
		double interp::bisection( const std::function<double( double, Args... )> &fun, double lb, double ub,
		double tol1, double tol2, Args... options )
		{
		if ( ub <= lb )
		throw std::logic_error( "Error: lb <= ub" );
		std::array<double, 2> range = { lb, ub };
		double x[101] = { 0 }, f[101] = { 0 };
		x[0] = lb;
		x[1] = ub;
		f[0] = fun( x[0], options... );
		f[1] = fun( x[1], options... );
		if ( fabs( f[0] ) < tol1 \|\| fabs( f[1] ) < tol1 )
		return fabs( f[0] ) < tol1 ? x[0] : x[1];
		if ( ( f[0] < 0 && f[1] < 0 ) \|\| ( f[0] > 0 && f[1] > 0 ) )
		throw std::logic_error( "Error: sign(f(lb)) == sign(f(ub))" );
		size_t i = 2;
		while ( true ) {
		// Get the next guess
		x[i] = bisection_coeff( i, x, f, range );
		if ( ( range[1] - range[0] ) < tol2 )
		break;
		// Evaluate f
		f[i] = fun( x[i], options... );
		if ( fabs( f[i] ) < tol1 )
		break;
		if ( i >= 100 )
		throw std::logic_error( "Excessive number of iterations" );
		i++;
		}
		return x[i];
		auto fun2 = [fun, options...]( double x ) { return fun( x, options... ); };
		return interp::bisection( fun2, lb, ub, tol1, tol2 );
		}

src/CMakeLists.txt

+50 −64

File changed.

Preview size limit exceeded, changes collapsed.

src/CreateImage.cpp

+14 −11

File changed.

Preview size limit exceeded, changes collapsed.