Page MenuHome
Differential D6807 Diff 21646 extern/opennurbs/opennurbs_optimize.cpp

Changeset View

Standalone View

View Options

extern/opennurbs/opennurbs_optimize.cpp

  • This file was added.
/* $NoKeywords: $ */
/*
//
// Copyright (c) 1993-2012 Robert McNeel & Associates. All rights reserved.
// OpenNURBS, Rhinoceros, and Rhino3D are registered trademarks of Robert
// McNeel & Associates.
//
// THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED WARRANTY.
// ALL IMPLIED WARRANTIES OF FITNESS FOR ANY PARTICULAR PURPOSE AND OF
// MERCHANTABILITY ARE HEREBY DISCLAIMED.
//
// For complete openNURBS copyright information see <http://www.opennurbs.org>.
//
////////////////////////////////////////////////////////////////
*/
#include "opennurbs.h"
#if !defined(ON_COMPILING_OPENNURBS)
// This check is included in all opennurbs source .c and .cpp files to insure
// ON_COMPILING_OPENNURBS is defined when opennurbs source is compiled.
// When opennurbs source is being compiled, ON_COMPILING_OPENNURBS is defined
// and the opennurbs .h files alter what is declared and how it is declared.
#error ON_COMPILING_OPENNURBS must be defined when compiling opennurbs
#endif
int ON_FindLocalMinimum(
int (*f)(void*,double,double*,double*), void* farg,
double ax, double bx, double cx,
double rel_stepsize_tol, double abs_stepsize_tol, int max_it,
double *t_addr
)
/* Use Brent's algorithm (with derivative) to Find a (local) minimum of a function
*
* INPUT:
* ax, bx, cx a bracketed minimum satisfying conditions 1 and 2.
* 1) either ax < bx < cx or cx < bx < ax.
* 2) f(bx) < f(ax) and f(bx) < f(ax).
* farg
* pointer passed to function f()
* f
* evaluation function with prototype
* int f(void* farg,double t,double* ft,double* dft)
* f(farg,t,&ft,&dft) should compute ft = value of function at t
* and dft = value of derivative at t.
* -1: failure
* 0: success
* 1: |f(x)| is small enough - TL_NRdbrent() will return *t_addr = x
* and the return code 1.
* rel_stepsize_tol, abs_stepsize_tol (0 < rel_stepsize_tol < 1 and 0 < abs_stepsize_tol)
* rel_stepsize_tol is a fractional tolerance and abs_stepsize_tol is an absolute tolerance
* that determine the minimum step size for a given iteration.
* minimum delta t = rel_stepsize_tol*|t| + abs_stepsize_tol.
* When in doubt, use
* rel_stepsize_tol = ON_EPSILON
* abs_stepsize_tol = 1/2*(desired absolute precision for *t_addr).
* max_it ( >= 2)
* maximum number of iterations to permit (when in doubt use 100)
* Closest Point to bezier minimizations typically take < 30
* iterations.
*
* OUTPUT:
* *t_addr abcissa of a local minimum between ax and cx.
* 0: failure
* 1: success
* 2: After max_iteration_cnt iterations the tolerance restrictions
* where not satisfied. Try increasing max_it, rel_stepsize_tol and/or abs_stepsize_tol
* or use the value of (*t_addr) with extreme caution.
*/
{
// See Numerical Recipes in C's dbrent() for a description of the basic algorithm
int rc,ok1,ok2;
double a,b,d,d1,d2,du,dv,dw,dx,e,fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,v,w,x,xm;
d=e=0.0;
if ( 0 == t_addr )
{
ON_ERROR("t_addr is nullptr");
return 0;
}
*t_addr = bx;
if ( max_it < 2 )
{
ON_ERROR("max_it must be >= 2");
return 0;
}
if ( !ON_IsValid(rel_stepsize_tol) || rel_stepsize_tol <= 0.0 || rel_stepsize_tol >= 1.0 )
{
ON_ERROR("rel_stepsize_tol must be strictly between 0.0 and 1.0");
return 0;
}
if ( !ON_IsValid(abs_stepsize_tol) || abs_stepsize_tol <= 0.0 )
{
ON_ERROR("abs_stepsize_tol must be > 0");
return 0;
}
a=(ax < cx ? ax : cx);
b=(ax > cx ? ax : cx);
x=w=v=bx;
rc = f(farg,x,&fx,&dx);
if (rc) {
// f() returned nonzero return code which means we need to bailout
if ( rc < 0 ) {
ON_ERROR("ON_FindLocalMinimum() f() failed to evaluate.");
}
*t_addr = x;
return rc>0 ? 1 : 0; // return 1 means f() said result is good enough, return = 0 means f() failed
}
fw=fv=fx;
dw=dv=dx;
while(max_it--) {
xm=0.5*(a+b);
tol1=rel_stepsize_tol*fabs(x)+abs_stepsize_tol;
tol2=2.0*tol1;
if (fabs(x-xm) <= (tol2-0.5*(b-a))) {
// further adjustments to x are smaller than stepsize tolerance
*t_addr=x;
return 1;
}
if (fabs(e) > tol1) {
d1=2.0*(b-a);
d2=d1;
if (dw != dx) d1=(w-x)*dx/(dx-dw);
if (dv != dx) d2=(v-x)*dx/(dx-dv);
u1=x+d1;
u2=x+d2;
ok1 = (a-u1)*(u1-b) > 0.0 && dx*d1 <= 0.0;
ok2 = (a-u2)*(u2-b) > 0.0 && dx*d2 <= 0.0;
olde=e;
e=d;
if (ok1 || ok2) {
if (ok1 && ok2)
d=(fabs(d1) < fabs(d2) ? d1 : d2);
else if (ok1)
d=d1;
else
d=d2;
if (fabs(d) <= fabs(0.5*olde)) {
u=x+d;
if (u-a < tol2 || b-u < tol2)
{d = (xm >= x) ? tol1 : -tol1;}
} else {
d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
}
} else {
d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
}
} else {
d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
}
if (fabs(d) >= tol1) {
u=x+d;
rc = f(farg,u,&fu,&du);
}
else {
u = (d >= 0.0) ? x+tol1 : x-tol1;
rc = f(farg,u,&fu,&du);
if (rc >= 0 && fu > fx) {
// tweaking x any more increases function value - x is a numerical minimum
*t_addr=x;
return 1;
}
}
if (rc) {
// f() returned nonzero return code which means we need to bailout
if ( rc < 0 ) {
ON_ERROR("ON_FindLocalMinimum() f() failed to evaluate.");
}
else {
*t_addr = (fu < fx) ? u : x;
}
return rc>0 ? 1 : 0;
}
if (fu <= fx) {
if (u >= x) a=x; else b=x;
v=w;fv=fw;dv=dw;
w=x;fw=fx;dw=dx;
x=u;fx=fu;dx=du;
} else {
if (u < x) a=u; else b=u;
if (fu <= fw || w == x) {
v=w;fv=fw;dv=dw;
w=u;fw=fu;dw=du;
} else if (fu < fv || v == x || v == w) {
v=u;fv=fu;dv=du;
}
}
}
*t_addr = x; // best known answer
ON_ERROR("ON_FindLocalMinimum() failed to converge");
return 2; // 2 means we failed to converge
}
ON_LocalZero1::ON_LocalZero1()
: m_t0(ON_UNSET_VALUE)
, m_t1(ON_UNSET_VALUE)
, m_f_tolerance(0.0)
, m_t_tolerance(0.0)
, m_k(nullptr)
, m_k_count(0)
, m_s0(ON_UNSET_VALUE)
, m_f0(ON_UNSET_VALUE)
, m_s1(ON_UNSET_VALUE)
, m_f1(ON_UNSET_VALUE)
{}
ON_LocalZero1::~ON_LocalZero1()
{}
bool
ON_LocalZero1::BracketZero( double s0, double f0,
double s1, double f1,
int level )
{
double s, f, d;
// private helper for FindSearchDomain()
if ( (f0 <= 0.0 && f1 >= 0.0) || (f0 >= 0.0 && f1 <= 0.0)
|| fabs(f0) <= m_f_tolerance || fabs(f1) <= m_f_tolerance ) {
m_t0 = s0;
m_t1 = s1;
return true;
}
if ( level++ <= 8 ) {
s = 0.5*s0+s1;
if ( s0 < s && s < s1 && Evaluate(s,&f,&d,0) ) {
if ( f*d >= 0.0 ) {
// search left side first
if ( BracketZero(s0,f0,s,f,level ) ) {
m_s0 = s0;
m_f0 = f0;
m_s1 = s;
m_f1 = f;
return true;
}
if ( BracketZero(s,f,s1,f1,level ) ) {
m_s0 = s;
m_f0 = f;
m_s1 = s1;
m_f1 = f1;
return true;
}
}
else {
// search right side first
if ( BracketZero(s,f,s1,f1,level ) ) {
m_s0 = s;
m_f0 = f;
m_s1 = s1;
m_f1 = f1;
return true;
}
if ( BracketZero(s0,f0,s,f,level ) ) {
m_s0 = s0;
m_f0 = f0;
m_s1 = s;
m_f1 = f;
return true;
}
}
}
}
return false;
}
bool
ON_LocalZero1::BracketSpan( double s0, double f0, double s1, double f1 )
{
int i0, i1, i;
double fm, fp;
bool rc = true;
if ( m_k && m_k_count >= 3 ) {
i0 = ON_SearchMonotoneArray(m_k,m_k_count,s0);
if ( i0 < 0 )
i0 = 0;
i1 = ON_SearchMonotoneArray(m_k,m_k_count,s1);
if ( i1 >= m_k_count )
i1 = m_k_count-1;
while ( i1 >= 0 && s1 == m_k[i1] ) {
i1--;
}
i0++;
while ( i0 < m_k_count-1 && m_k[i0] == m_k[i0+1] )
i0++;
if ( i0 <= i1 ) {
// we have s0 < m_k[i0] <= ... <= m_k[i1] < s1
Evaluate( m_k[i0], &fm, nullptr,-1 ); // gaurd against C0 discontinuities
Evaluate( m_k[i0], &fp, nullptr, 1 );
if ( (f0 <= 0.0 && fm >= 0.0) || (f0 >= 0.0 && fm <= 0.0) ) {
m_s1 = m_k[i0];
m_f1 = fm;
}
else if ( (f1 <= 0.0 && fp >= 0.0) || (f1 >= 0.0 && fp <= 0.0) ) {
m_s0 = m_k[i0];
m_f0 = fp;
if ( i0 < i1 ) {
Evaluate( m_k[i1], &fm, nullptr, -1 );
Evaluate( m_k[i1], &fp, nullptr, 1 );
if ( (f1 <= 0.0 && fp >= 0.0) || (f1 >= 0.0 && fp <= 0.0) ) {
m_s0 = m_k[i1];
m_f0 = fp;
}
else if ( (f0 <= 0.0 && fm >= 0.0) || (f0 >= 0.0 && fm <= 0.0) ) {
m_s1 = m_k[i1];
m_f1 = fm;
// we have s0 = m_k[i0] < m_k[i1] = s1 and a bonafide sign change
while (i0+1 < i1) {
// bisect m_k[i0],...,m_k[i1] until we get a sign change between
// m_k[i],m_k[i+1]. We need to do this in order to make sure
// we are passing a C2 function to the repeated zero finders.
i = (i0+i1)>>1;
Evaluate( m_k[i], &fm, nullptr, -1 );
Evaluate( m_k[i], &fp, nullptr, 1 );
if ( (f0 <= 0.0 && fm >= 0.0) || (f0 >= 0.0 && fm <= 0.0) ) {
m_s1 = m_k[i];
m_f1 = fm;
i1 = i;
while ( i1>0 && m_k[i1-1]==m_k[i1])
i1--;
}
else if ( (f1 <= 0.0 && fp >= 0.0) || (f1 >= 0.0 && fp <= 0.0) ) {
m_s0 = m_k[i];
m_f0 = fp;
i0 = i;
while ( i0 < m_k_count-2 && m_k[i0] == m_k[i0+1] )
i0++;
}
else {
// discontinuous sign change across m_k[i]
rc = false;
break;
}
}
}
else {
// discontinuous sign change across m_k[i1]
rc = false;
}
}
}
else {
// discontinuous sign change across m_k[i0]
rc = false;
}
}
}
return rc;
}
bool ON_LocalZero1::FindZero( double* t )
{
// Find values of m_s0 and m_s1 between m_t0 and m_t1 such that
// f(m_t0) and f(m_t1) have different signs
if ( !ON_IsValid(m_t0) )
{
if ( !ON_IsValid(m_t1) )
{
ON_ERROR("Illegal input - m_t0 and m_t1 are not valid.");
return false;
}
m_s0 = m_s1 = m_t1;
}
else if ( !ON_IsValid(m_t1) )
{
m_s0 = m_s1 = m_t0;
}
else if ( m_t0 <= m_t1 )
{
m_s0 = m_t0;
m_s1 = m_t1;
}
else if ( m_t1 < m_t0 )
{
m_s0 = m_t1;
m_s1 = m_t0;
}
else
{
ON_ERROR("Illegal input - m_t0 and m_t1 are not valid.");
return false;
}
if ( m_s0 == m_s1 )
{
if ( Evaluate( m_s0, &m_f0, nullptr, 1 ) )
{
m_f1 = m_f0;
if ( fabs(m_f0) <= m_f_tolerance )
{
*t = m_s0;
return true;
}
ON_ERROR("Illegal input - m_t0 = m_t1 and the function value is not zero at m_t0.");
return false;
}
ON_ERROR("Evaluation failed.");
return false;
}
if ( !Evaluate( m_s0, &m_f0, nullptr, 1 ) )
{
ON_ERROR("Evaluation failed at m_s0.");
return false;
}
if ( !Evaluate( m_s1, &m_f1, nullptr, -1 ) )
{
ON_ERROR("Evaluation failed at m_s1.");
return false;
}
if ( !BracketZero( m_s0, m_f0, m_s1, m_f1 ) )
{
ON_ERROR("Unable to bracket a zero of the function.");
return false;
}
if ( fabs(m_f0) <= m_f_tolerance && fabs(m_f0) <= fabs(m_f1) )
{
// |f(s0)| <= user specified stopping tolerance
*t = m_s0;
return true;
}
if ( fabs(m_f1) <= m_f_tolerance )
{
// |f(s1)| <= user specified stopping tolerance
*t = m_s1;
return true;
}
if ( !BracketSpan( m_s0, m_f0, m_s1, m_f1 ) )
{
ON_ERROR("Unable to bracket the function in a span of m_k[]. m_k[] may be invalid.");
return false;
}
if ( !NewtonRaphson( m_s0, m_f0, m_s1, m_f1, 128, t ) )
{
ON_ERROR("Newton-Raphson failed to converge. Is your function C2?");
return false;
}
return true;
}
bool ON_LocalZero1::NewtonRaphson( double s0, double f0,
double s1, double f1,
int maxit, double* t )
{
// private function - input must satisfy
//
// 1) t is not nullptr
//
// 2) maxit >= 2
//
// 3) s0 != s1,
//
// 4) f0 = f(s0), f1 = f(s1),
//
// 5) either f0 < 0.0 and f1 > 0.0, or f0 > 0.0 and f1 < 0.0
//
// 6) f() is C0 on [s0,s1] and C2 on (s0,s1)
//
// 7) ds_tolerance >= 0.0 - the search for a zero will stop
// if the zero in bracketed in an interval that is no
// larger than ds_tolerance.
//
// When in doubt, ds_tolerance = (fabs(s0) + fabs(s1))*ON_EPSILON
// works well.
double s, f, d, x, ds, prevds;
if ( fabs(f0) <= m_f_tolerance && fabs(f0) <= fabs(f1) ) {
// |f(s0)| <= user specified stopping tolerance
*t = s0;
return true;
}
if ( fabs(f1) <= m_f_tolerance ) {
// |f(s1)| <= user specified stopping tolerance
*t = s1;
return true;
}
if ( f0 > 0.0 ) {
x = s0; s0 = s1; s1 = x;
x = f0; f0 = f1; f1 = x;
}
s = 0.5*(s0+s1);
if ( !Evaluate( s, &f, &d, 0 ) ) {
*t = (fabs(f0) <= fabs(f1)) ? s0 : s1;
return false;
}
if ( fabs(f) <= m_f_tolerance ) {
// |f(s)| <= user specified stopping tolerance
*t = s;
return true;
}
if ( f1 <= 0.0 ) {
*t = (fabs(f0) <= fabs(f1)) ? s0 : s1;
return false;
}
ds = fabs(s1-s0);
prevds = 0.0;
while (maxit--) {
if ( (f+(s0-s)*d)*(f+(s1-s)*d) > 0.0 // true if NR's line segment doesn't cross zero inside interval
|| fabs(2.0*f) > fabs(prevds*d) // true if expected decrease in function value from previous step didn't happen
) {
// bisect
prevds = ds;
ds = 0.5*(s1-s0);
s = s0+ds;
if ( s == s0 ) {
// interval is too small to be divided using double division
if ( fabs(f1) < fabs(f0) ) {
s = s1;
}
*t = s;
return true;
}
}
else {
// Newton iterate
prevds = ds;
ds = -f/d;
x = s;
s += ds;
if ( s == x ) {
// Newton step size < smallest double than can be added to s
if ( fabs(f0) < fabs(f) ) {
f = f0;
s = s0;
}
if ( fabs(f1) < fabs(f) ) {
s = s1;
}
*t = s;
return true;
}
}
if ( !Evaluate( s, &f, &d, 0 ) ) {
*t = (fabs(f0) <= fabs(f1)) ? s0 : s1; // emergency bailout
return false;
}
if ( fabs(f) <= m_f_tolerance ) {
// |f(s)| <= user specified stopping tolerance
if ( fabs(f0) < fabs(f) ) {
f = f0;
*t = s0;
}
if ( fabs(f1) < fabs(f) ) {
*t = s1;
}
return true;
}
if ( f < 0.0 ) {
f0 = f; // f0 needed for emergency bailout
s0 = s;
}
else { // f > 0.0
f1 = f; // f1 needed for emergency bailout
s1 = s;
}
if ( fabs(s1-s0) <= m_t_tolerance ) {
// a root has been bracketed to an interval that is small enough
// to satisify user.
*t = (fabs(f0) <= fabs(f1)) ? s0 : s1;
return true;
}
}
*t = (fabs(f0) <= fabs(f1)) ? s0 : s1; // emergency bailout
return false;
}