~ubuntu-branches/debian/stretch/uncertainties/stretch : revision 1

1

'''

2

Mathematical operations that generalize many operations from the

3

standard math module so that they also work on numbers with

4

uncertainties.

5

6

Examples:

7

8

from umath import sin

9

10

# Manipulation of numbers with uncertainty:

11

x = uncertainties.num_with_uncert((3, 0.1))

12

print sin(x) # prints 0.141120008...+/-0.098999...

13

14

# The umath functions also work on regular Python floats:

15

print sin(3) # prints 0.141120008... This is a Python float.

16

17

Importing all the functions from this module into the global namespace

18

is possible. This is encouraged when using a Python shell as a

19

calculator. Example:

20

21

import uncertainties

22

from uncertainties.umath import * # Imports tan(), etc.

23

24

x = uncertainties.num_with_uncert((3, 0.1))

25

print tan(x) # tan() is the uncertainties.umath.tan function

26

27

The numbers with uncertainty handled by this module are objects from

28

the uncertainties module, from either the Variable or the

29

AffineScalarFunc class.

30

31

32

Please send feature requests, bug reports, or feedback to this address.

33

34

This software is released under a dual license. (1) The GNU General

35

Public License version 2. (2) Any other license, as long as it is

36

obtained from the original author.'''

37

38

from __future__ import division # Many analytical derivatives depend on this

39

40

# Standard modules

41

import math

42

import sys

43

import itertools

44

import functools

45

import inspect

46

47

# Non-standard modules

48

import uncertainties

49

50

###############################################################################

51

52

# We wrap the functions from the math module so that they keep track of

53

# uncertainties by returning a AffineScalarFunc object.

54

55

# Some functions from the math module cannot be simply adapted to work

56

# with AffineScalarFunc objects (either as their result or as their

57

# arguments):

58

59

# (1) Some functions return a result of a type whose value and

60

# variations (uncertainties) cannot be represented by AffineScalarFunc

61

# (e.g., math.frexp, which returns a tuple). The exception raised

62

# when not wrapping them with wrap() is more obvious than the

63

# one obtained when wrapping them (in fact, the wrapped functions

64

# attempts operations that are not supported, such as calculation a

65

# subtraction on a result of type tuple).

66

67

# (2) Some functions don't take scalar arguments (which can be varied

68

# during differentiation): math.fsum... Such functions can either be:

69

70

# - wrapped in a special way in wrap_math_functions()

71

72

# - excluded from wrapping by adding their name to no_std_wrapping

73

74

# - wrapped in a general way in wrap_math_functions(); in this case,

75

# the function does not have to be mentioned in any location in this

76

# code. The function should function normally, except possibly when

77

# given AffineScalarFunc arguments.

78

79

no_std_wrapping = ['frexp', 'modf', 'fsum'] # Exclude from standard wrapping

80

81

std_wrapped_math_funcs = [] # Math functions wrapped in the standard way

82

83

non_std_wrapped_math_funcs = [] # Math functions wrapped in a non-standard way

84

85

# Function that copies the relevant attributes from generalized

86

# functions from the math module:

87

wraps = functools.partial(functools.update_wrapper,

88

assigned=('__doc__', '__name__'))

89

90

########################################

91

# Special cases:

92

93

# fsum takes a single argument, which cannot be differentiated.

94

# However, each of the arguments inside this single list can

95

# be a variable. We handle this in a specific way:

96

97

try:

98

original_func = math.fsum # Shortcut

99

except AttributeError:

100

pass # math.fsum introduced in Python 2.6

101

else:

102

103

# We wrap math.fsum

104

105

def wrapped_fsum():

106

"""

107

Return an uncertainty aware version of math.fsum, which must

108

be contained in _original_func.

109

"""

110

111

# The fsum function is flattened, in order to use the

112

# wrap() wrapper:

113

114

flat_fsum = lambda *args: original_func(args)

115

116

flat_fsum_wrap = uncertainties.wrap(

117

flat_fsum, itertools.repeat(lambda *args: 1))

118

119

return wraps(lambda arg_list: flat_fsum_wrap(*arg_list),

120

original_func)

121

122

fsum = wrapped_fsum()

123

124

# Wrapping already done:

125

non_std_wrapped_math_funcs.append('fsum')

126

127

########################################

128

# Wrapping of built-in math functions not in no_std_wrapping:

129

130

# Fixed formulas for the derivatives of some functions from the math

131

# module (some functions might not be present in all version of

132

# Python). Singular points are not taken into account. The user

133

# should never give "large" uncertainties: problems could only appear

134

# if this assumption does not hold.

135

136

# Functions not mentioned in _fixed_derivatives have their derivatives

137

# calculated numerically.

138

139

# Functions that have singularities (possibly at infinity) benefit

140

# from analytical calculations (instead of the default numerical

141

# calculation). Even slowly varying functions (e.g., abs()) yield

142

# more precise results when differentiated analytically, because of

143

# the loss of precision in numerical calculations.

144

145

#def log_1arg_der(x):

146

# """

147

# Derivative of log(x) (1-argument form).

148

# """

149

# return 1/x

150

151

def log_der0(*args):

152

"""

153

Derivative of math.log() with respect to its first argument.

154

155

Works whether 1 or 2 arguments are given.

156

"""

157

if len(args) == 1:

158

return 1/args[0]

159

else:

160

return 1/args[0]/math.log(args[1]) # 2-argument form

161

162

# The following version goes about as fast:

163

164

## A 'try' is used for the most common case because it is fast when no

165

## exception is raised:

166

#try:

167

# return log_1arg_der(*args) # Argument number check

168

#except TypeError:

169

# return 1/args[0]/math.log(args[1]) # 2-argument form

170

171

fixed_derivatives = {

172

# In alphabetical order, here:

173

'acos': [lambda x: -1/math.sqrt(1-x**2)],

174

'acosh': [lambda x: 1/math.sqrt(x**2-1)],

175

'asin': [lambda x: 1/math.sqrt(1-x**2)],

176

'asinh': [lambda x: 1/math.sqrt(1+x**2)],

177

'atan': [lambda x: 1/(1+x**2)],

178

'atan2': [lambda y, x: x/(x**2+y**2), # Correct for x == 0

179

lambda y, x: -y/(x**2+y**2)], # Correct for x == 0

180

'atanh': [lambda x: 1/(1-x**2)],

181

'ceil': [lambda x: 0],

182

'copysign': [lambda x, y: (1 if x >= 0 else -1) * math.copysign(1, y),

183

lambda x, y: 0],

184

'cos': [lambda x: -math.sin(x)],

185

'cosh': [math.sinh],

186

'degrees': [lambda x: math.degrees(1)],

187

'exp': [math.exp],

188

'fabs': [lambda x: 1 if x >= 0 else -1],

189

'floor': [lambda x: 0],

190

'hypot': [lambda x, y: x/math.hypot(x, y),

191

lambda x, y: y/math.hypot(x, y)],

192

'ldexp': [lambda x, y: 2**y,

193

# math.ldexp only accepts an integer as its second

194

# argument:

195

None],

196

'log': [log_der0,

197

lambda x, y: -math.log(x, y)/y/math.log(y)],

198

'log10': [lambda x: 1/x/math.log(10)],

199

'log1p': [lambda x: 1/(1+x)],

200

'pow': [lambda x, y: y*math.pow(x, y-1),

201

lambda x, y: math.log(x) * math.pow(x, y)],

202

'radians': [lambda x: math.radians(1)],

203

'sin': [math.cos],

204

'sinh': [math.cosh],

205

'sqrt': [lambda x: 0.5/math.sqrt(x)],

206

'tan': [lambda x: 1+math.tan(x)**2],

207

'tanh': [lambda x: 1-math.tanh(x)**2]

208

}

209

210

# Many built-in functions in the math module are wrapped with a

211

# version which is uncertainty aware:

212

213

this_module = sys.modules[__name__]

214

215

# We do not want to wrap module attributes such as __doc__, etc.:

216

for (name, func) in inspect.getmembers(math, inspect.isbuiltin):

217

218

if name in no_std_wrapping:

219

continue

220

221

if name in fixed_derivatives:

222

derivatives = fixed_derivatives[name]

223

else:

224

derivatives = None # Means: numerical calculation required

225

setattr(this_module, name,

226

wraps(uncertainties.wrap(func, derivatives), func))

227

std_wrapped_math_funcs.append(name)

228

229

###############################################################################

230

# Unit tests

231

232

def test_fixed_derivatives_math_funcs():

233

"""

234

Check of wrapped functions from the math module.

235

"""

236

237

for name in std_wrapped_math_funcs:

238

# print "Checking %s..." % name

239

func = getattr(this_module, name)

240

# Numerical derivatives of func: the nominal value of func() results

241

# is used as the underlying function:

242

numerical_derivatives = uncertainties.NumericalDerivatives(

243

lambda *args: func(*args).nominal_value)

244

uncertainties._compare_derivatives(func, numerical_derivatives)

245

246

def test_compound_expression():

247

"""

248

Test equality between different formulas.

249

"""

250

251

x = uncertainties.num_with_uncert((3, 0.1))

252

253

# Prone to numerical errors (but not much more than floats):

254

assert tan(x) == sin(x)/cos(x)

255

256

257

def test_numerical_example():

258

"Test specific numerical examples"

259

260

x = uncertainties.num_with_uncert((3.14, 0.01))

261

result = sin(x)

262

# In order to prevent big errors such as a wrong, constant value

263

# for all analytical and numerical derivatives, which would make

264

# test_fixed_derivatives_math_funcs() succeed despite incorrect

265

# calculations:

266

assert ("%.6f +/- %.6f" % (result.nominal_value, result.std_dev())

267

== "0.001593 +/- 0.010000")

268

269

# Regular calculations should still work:

270

assert("%.11f" % math.sin(3) == "0.14112000806")

271

272

def test_monte_carlo_comparison():

273

"""

274

Full comparison to a Monte-Carlo calculation.

275

276

Both the nominal values and the covariances are compared between

277

the direct calculation performed in this module and a Monte-Carlo

278

simulation.

279

"""

280

281

try:

282

import numpy

283

import numpy.random

284

except ImportError:

285

raise Exception("Test not performed because Numpy is not available")

286

287

# Works on numpy.arrays of Variable objects (whereas numpy.sin()

288

# does not):

289

sin_array_uncert = numpy.vectorize(sin)

290

291

# Example expression (with correlations, and multiple variables combined

292

# in a non-linear way):

293

def function(x, y):

294

"""

295

Function that takes two NumPy arrays of the same size.

296

"""

297

# The uncertainty due to x is about equal to the uncertainty

298

# due to y:

299

return 10 * x**2 - x * sin_array_uncert(y**3)

300

301

x = uncertainties.num_with_uncert((0.2, 0.01))

302

y = uncertainties.num_with_uncert((10, 0.001))

303

function_result_this_module = function(x, y)

304

nominal_value_this_module = function_result_this_module.nominal_value

305

306

# Covariances "f*f", "f*x", "f*y":

307

covariances_this_module = numpy.array(uncertainties.covariance_matrix(

308

(x, y, function_result_this_module)))

309

310

def monte_carlo_calc(n_samples, ):

311

"""

312

Calculate function(x, y) on n_samples and returns the median,

313

and the covariances between (x, y, function(x, y)).

314

"""

315

# Result of a Monte-Carlo simulation:

316

x_samples = numpy.random.normal(x.nominal_value, x.std_dev(), n_samples)

317

y_samples = numpy.random.normal(y.nominal_value, y.std_dev(), n_samples)

318

function_samples = function(x_samples, y_samples)

319

320

cov_mat = numpy.cov([x_samples, y_samples], function_samples)

321

322

return (numpy.median(function_samples), cov_mat)

323

324

(nominal_value_samples, covariances_samples) = monte_carlo_calc(100000)

325

326

327

## Comparison between both results:

328

329

# The covariance matrices must be close:

330

331

# We rely on the fact that covariances_samples very rarely has

332

# null elements:

333

334

assert numpy.vectorize(uncertainties._numbers_close)(

335

covariances_this_module,

336

covariances_samples,

337

0.05).all(), (

338

"The covariance matrices do not coincide between"

339

" the Monte-Carlo simulation and the direct calculation:\n"

340

"* Monte-Carlo:\n%s\n* Direct calculation:\n%s" %

341

(covariances_samples, covariances_this_module)

342

)

343

344

# The nominal values must be close:

345

assert uncertainties._numbers_close(

346

nominal_value_this_module,

347

nominal_value_samples,

348

# The scale of the comparison depends on the standard

349

# deviation: the nominal values can differ by a fraction of

350

# the standard deviation:

351

sqrt(covariances_samples[2, 2])

352

/ abs(nominal_value_samples) * 0.5), (

353

"The nominal value (%f) does not coincide with that of"

354

" the Monte-Carlo simulation (%f), for a standard deviation of %f." %

355

(nominal_value_this_module,

356

nominal_value_samples,

357

sqrt(covariances_samples[2, 2]))

358

)

359

360

361

def test_math_module():

362

"Operations with the math module"

363

364

x = uncertainties.num_with_uncert((-1.5, 0.1))

365

366

# The exponent must not be differentiated, when calculating the

367

# following (the partial derivative with respect to the exponent

368

# is not defined):

369

assert (x**2).nominal_value == 2.25

370

371

# Regular operations are chosen to be unchanged:

372

assert isinstance(math.sin(3), float)

373

374

# Python >=2.6 functions:

375

376

if sys.version_info[:2] >= (2, 6):

377

378

# factorial must not be "damaged" by this module; in particular,

379

# as of Python 2.6, it does not accept non integral values, and

380

# must therefore not be differentiated:

381

assert math.factorial(4) == 24

382

383

# Boolean functions:

384

assert not isinf(x)

385

386

# Comparison, possibly between an AffineScalarFunc object and a

387

# boolean, which makes things more difficult for this code:

388

assert isinf(x) == False

389

390

# fsum is special because it does not take a fixed number of

391

# variables:

392

assert fsum([x, x]).nominal_value == -3

393

394

###############################################################################

395

# Exported functions:

396

397

__all__ = std_wrapped_math_funcs + non_std_wrapped_math_funcs

398