In [3]:
higham(2, 30)
Out[3]:
In [5]:
higham(2, 5)
Out[5]:
In [17]:
higham(2, 60)
Out[17]:
In [14]:
sqrt_eps = 1.5e-8
t = 2
L =
In [15]:
higham(t + sqrt_eps, L) - higham(t, L)
Out[15]:
In [ ]:
def delta(f, h, t0):
return fabs(f(t0 - h) - 2*f(t0) + f(t0 + h))
def mu(delta):
tau1 = 100
tau2 = .1
ha = epsf**.25
In [2]:
import numpy as np
In [3]:
fval = np.array([0, 1, 2, 3, 4])
In [5]:
fval.size
Out[5]:
In [6]:
fval[0:2]
Out[6]:
In [7]:
from math import sqrt, fabs
import numpy as np
def higham(t, L):
f = t
for k in range(1, L + 1):
f = sqrt(f)
for k in range(1, L + 1):
f = f**2
f = f**2
return f
#
# Determines the noise of a function from the function values
#
# [fnoise,level,inform] = ECnoise(nf,fval)
#
# The user must provide the function value at nf equally-spaced points.
# For example, if nf = 7, the user could provide
#
# f(x-3h), f(x-2h), f(x-h), f(x), f(x+h), f(x+2h), f(x+3h)
#
# in the array fval. Although nf >= 4 is allowed, the use of at least
# nf = 7 function evaluations is recommended.
#
# Noise will not be detected by this code if the function values differ
# in the first digit.
#
# If noise is not detected, the user should increase or decrease the
# spacing h according to the ouput value of inform. In most cases,
# the subroutine detects noise with the initial value of h.
#
# On exit:
# fnoise is set to an estimate of the function noise;
# fnoise is set to zero if noise is not detected.
#
# level is set to estimates for the noise. The k-th entry is an
# estimate from the k-th difference.
#
# inform is set as follows:
# inform = 1 Noise has been detected.
# inform = 2 Noise has not been detected; h is too small.
# Try 100*h for the next value of h.
# inform = 3 Noise has not been detected; h is too large.
# Try h/100 for the next value of h.
#
# Argonne National Laboratory
# Jorge More' and Stefan Wild. November 2009.
def ECNoise(nf, fval):
level = np.zeros((nf-1))
dsgn = np.zeros((nf-1))
fnoise = 0.0
gamma = 1.0 # = gamma(0)
# Compute the range of function values.
fmin = np.amin(fval)
fmax = np.amax(fval)
if (fmax-fmin)/max(fabs(fmax), fabs(fmin)) > .1:
inform = 3
return fnoise, level, inform
# Construct the difference table.
for j in range(nf-1):
for i in range(nf-(j+1)):
fval[i] = fval[i+1] - fval[i]
# h is too small only when half the function values are equal.
if (j==0 and sum([fval[k] == 0 for k in range(nf - 1)]) >= nf/2):
inform = 2
return fnoise, level, inform
gamma = 0.5*((j+1.)/(2.*(j+1.)-1.))*gamma
# Compute the estimates for the noise level.
level[j] = sqrt(gamma*np.mean(np.square(fval[0:nf-(j+1)])))
# Determine differences in sign.
emin = np.amin(fval[0:nf-(j+1)])
emax = np.amax(fval[0:nf-(j+1)])
if (emin*emax < 0.0):
dsgn[j] = 1
# Determine the noise level.
for k in range(nf-3):
emin = np.amin(level[k:k+3])
emax = np.amax(level[k:k+3])
if (emax<=4*emin and dsgn[k]):
fnoise = level[k]
inform = 1
return fnoise, level, inform
# If noise not detected then h is too large.
inform = 3
return fnoise, level, inform
In [8]:
fval = np.zeros((7))
sqrt_eps = 1.5e-8
initial_h = 100*sqrt_eps
for i,j in zip(range(-3,4), range(0,7)):
fval[j] = higham(2 + i*initial_h, 30)
fnoise, level, inform = ECNoise(7, fval)
In [9]:
inform
Out[9]:
In [10]:
fnoise
Out[10]:
In [27]:
level
Out[27]:
In [ ]: