%% options copyright owner = Dirk Krause copyright year = 2019-xxxx SPDX-License-Identifier: BSD-3-Clause %% header /** @file dk4itmin.h Iteration algorithms for minimum search. Iteration can be used on functions unimodal on interval [a;b] only. This means the function must have exactly one local minimum on interval [a;b], no local maximum. To search for a maximum of f(x) simply search for the minimum of the 0-f(x) function. */ #ifndef DK4ITER_H_INCLUDED #if DK4_BUILDING_DKTOOLS4 #include "dk4iter.h" #else #include #endif #endif #ifdef __cplusplus extern "C" { #endif /** Run iteration for minimum search. @param d Address of result array, 2 elements: x-value and y-value. @param pp Address of variable to store number of passes, may be NULL. @param fct Function to find minimum for. The function must be unimodal on [a;b], the dk4itmin_interval() function can not check this. @param ps Parameter set for function, may be NULL. @param a Left interval border. @param b Right interval border. Required: f(c) #if DK4_HAVE_ASSERT_H #ifndef ASSERT_H_INCLUDED #include #define ASSERT_H_INCLUDED 1 #endif #endif #if DK4_HAVE_STDLIB_H #ifndef STDLIB_H_INCLUDED #include #define STDLIB_H_INCLUDED 1 #endif #endif #if DK4_HAVE_LIMITS_H #ifndef LIMITS_H_INCLUDED #include #define LIMITS_H_INCLUDED 1 #endif #endif #if DK4_HAVE_STDINT_H #ifndef STDINT_H_INCLUDED #include #define STDINT_H_INCLUDED 1 #endif #endif #ifndef DK4MEM_H_INCLUDED #include "dk4mem.h" #endif #ifndef DK4ITER_H_INCLUDED #include "dk4iter.h" #endif #ifndef DK4MIN_H_INCLUDED #include "dk4itmin.h" #endif #ifndef DK4MATH_H_INCLUDED #include "dk4math.h" #endif static int dk4itmin_check_tolerances( double eps_y, double y, double yold, double eps_x, double x, double xold ) { int back = 0; if ((!(isgreater(eps_y, 0.0))) || (isgreater(eps_y, fabs(y - yold)))) { if ((!(isgreater(eps_x, 0.0))) || (isgreater(eps_x, fabs(x - xold)))) { back = 1; } } return back; } static int dk4itmin_golden_section( double *pd, unsigned long *pp, dk4_iter_fct_t *fct, void const *ps, double a, double fa, double b, double fb, dk4_iter_ctx_t const *ctx ) { double c = 0.0; /* Probe position c */ double fc = 0.0; /* Function value at c */ double d = 0.0; /* Probe position d */ double fd = 0.0; /* Function value at d */ double f1 = 0.0; /* Factor for probe calculation */ double f2 = 0.0; /* Factor for probe calculation */ double x = 0.0; /* Current probe point */ double xold = 0.0; /* Previous probe point */ double y = 0.0; /* Current probe value */ double yold = 0.0; /* Previous probe value */ double max = 0.0; /* Maximum of fa and fb */ unsigned long passno = 0UL; /* Current pass number */ int res = 0; /* Operations result */ int cc = 1; /* 1=go on, 0=finished, -1=abort */ int back = DK4_ITER_RESULT_E_ARGS; $? "+ dk4itmin_golden_section" /* Initialize allowed maximum function value */ /* Initialize factors for probe calculation */ f1 = sqrt(5.0) - 1.0; $? ". f1 = %g", f1 f2 = 3.0 - sqrt(5.0); $? ". f2 = %g", f2 max = ((fa > fb) ? (fa) : (fb)); /* Initial probe points */ c = (f1 * a + f2 * b) / 2.0; d = (f2 * a + f1 * b) / 2.0; res = (*fct)(&fc, c, ps); if ((0 == res) || (!(dk4ma_is_finite(fc)))) { $? "! fc calculation" back = DK4_ITER_RESULT_E_FCT; goto finished; } res = (*fct)(&fd, d, ps); if ((0 == res) || (!(dk4ma_is_finite(fd)))) { $? "! fd calculation" back = DK4_ITER_RESULT_E_FCT; goto finished; } $? ". A ( %g , %g )", a, fa $? ". B ( %g , %g )", b, fb $? ". C ( %g , %g )", c, fc $? ". D ( %g , %g )", d, fd if (fc > max) { $? "! fc not unimodal" back = DK4_ITER_RESULT_E_CONV; goto finished; } if (fd > max) { $? "! fd not unimodal" back = DK4_ITER_RESULT_E_CONV; goto finished; } /* Iteration loop */ do { /* Increase pass number */ if (ULONG_MAX > passno) { passno++; } $? ". pass %lu", passno $? ". A ( %g , %g )", a, fa $? ". B ( %g , %g )", b, fb $? ". C ( %g , %g )", c, fc $? ". D ( %g , %g )", d, fd /* Save probe point and value */ xold = x; yold = y; /* Process the probe points */ if (fc < fd) { $? ". switch to [a;d]" b = d; d = c; c = (f1 * a + f2 * b) / 2.0; fb = fd; fd = fc; res = (*fct)(&fc, c, ps); if ((0 == res) || (!(dk4ma_is_finite(fc)))) { back = DK4_ITER_RESULT_E_FCT; $? "! fc calculation" cc = -1; } else { max = ((fa > fb) ? (fa) : (fb)); if (fc > max) { $? "! fc not unimodal" back = DK4_ITER_RESULT_E_CONV; cc = -1; } else { x = c; $? ". x = %g", x y = fc; $? ". y = %g", y } } } else { $? ". switch to [c;b]" a = c; c = d; d = (f2 * a + f1 * b) / 2.0; fa = fc; fc = fd; res = (*fct)(&fd, d, ps); if ((0 == res) || (!(dk4ma_is_finite(fd)))) { back = DK4_ITER_RESULT_E_FCT; $? "! fd calculation" cc = -1; } else { max = ((fa > fb) ? (fa) : (fb)); if (fd > max) { $? "! fd not unimodal" back = DK4_ITER_RESULT_E_CONV; cc = -1; } else { x = d; $? ". x = %g", x y = fd; $? ". y = %g", y } } } /* Check for reached tolerance */ if (1 == cc) { if (1UL < passno) { if (x == xold) { cc = 0; $? ". exact result found" } else { if (0 == ctx->exact) { res = dk4itmin_check_tolerances( ctx->eps_y, y, yold, ctx->eps_x, x, xold ); if (0 != res) { cc = 0; $? ". tolerances reached" } #if TRACE_DEBUG else { $? ". must continue" } #endif } } } #if TRACE_DEBUG else { $? ". first pass, cannot check yet" } #endif } /* Abort if number of passes exceeded */ if (1 == cc) { if ((0UL < ctx->maxpass) && (passno >= ctx->maxpass)) { back = DK4_ITER_RESULT_E_PASSES; cc = -1; $? "! too many passes" } #if TRACE_DEBUG else { $? ". can continue" } #endif } } while (1 == cc); /* Transfer results on success */ if (0 == cc) { $? ". success" pd[0] = x; pd[1] = y; back = DK4_ITER_RESULT_SUCCESS; if (NULL != pp) { *pp = passno; } } finished: $? "- dk4itmin_golden_section %d", back return back; } int dk4itmin_interval( double *d, unsigned long *pp, dk4_iter_fct_t *fct, void const *ps, double a, double b, dk4_iter_ctx_t const *ctx ) { dk4_iter_ctx_t mctx; /* Iteration context copy */ double fa = 0.0; /* Value at left border */ double fb = 0.0; /* Value at right border */ int res = 0; /* Operation result */ int back = DK4_ITER_RESULT_E_ARGS; $? "+ dk4itmin_interval" /* Check function cal arguments */ #if DK4_USE_ASSERT assert(NULL != d) assert(NULL != fct) assert(a < b) #endif if ((NULL == d) || (NULL == fct)) { $? "! d or fct" goto finished; } if (a >= b) { $? "! zero length interval" goto finished; } /* Calculate values for the three given points */ res = (*fct)(&fa, a, ps); if ((0 == res) || (!(dk4ma_is_finite(fa)))) { back = DK4_ITER_RESULT_E_FCT; $? "! function calculation for a" goto finished; } res = (*fct)(&fb, b, ps); if ((0 == res) || (!(dk4ma_is_finite(fb)))) { back = DK4_ITER_RESULT_E_FCT; $? "! function calculation for b" goto finished; } /* Copy or initialize context */ if (NULL != ctx) { $? ". use provided context" DK4_MEMCPY(&mctx,ctx,sizeof(dk4_iter_ctx_t)); } else { $? ". no context, use defaults" dk4iter_ctx_init(&mctx); } /* Choose algorithm depending on context At this time only golden section available */ back = dk4itmin_golden_section(d, pp, fct, ps, a, fa, b, fb, &mctx); finished: $? "- dk4itmin_interval %d", back return back; } /* vim: set ai sw=4 ts=4 : */