symbolicModuleTest.py
You can view and download this file on Github: symbolicModuleTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Tests for symbolic scalar, vector and matrix
5#
6# Author: Johannes Gerstmayr
7# Date: 2023-12-14
8#
9# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
10#
11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12import exudyn as exu
13from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
14import exudyn.graphics as graphics
15
16testIsActive = exu.sys.get('testIsActive', True)
17
18esym = exu.symbolic
19import numpy as np
20import math
21
22#reset counters:
23esym.Real.__newCount = 0
24esym.Real.__deleteCount = 0
25esym.Vector.__newCount = 0
26esym.Vector.__deleteCount = 0
27esym.Matrix.__newCount = 0
28esym.Matrix.__deleteCount = 0
29
30SymReal = esym.Real
31
32def Not(value):
33 return not value
34
35def PyRealName(name, value):
36 return value
37
38def PyNamedVector(name, value):
39 return np.array(value)
40
41def Evaluate(value):
42 #print(str(value))
43 if (type(value) == esym.Real
44 or type(value) == esym.Vector
45 or type(value) == esym.Matrix
46 ):
47 return value.Evaluate()
48 if type(value) == np.ndarray:
49 value = value.copy() #requires copy, to store results!
50 return value
51
52def Min(a,b):
53 return a if a<b else b
54def Max(a,b):
55 return a if a>b else b
56
57math.sign = np.sign #to have function in math
58math.Not = Not
59math.abs = np.abs
60math.mod = math.fmod
61math.min = Min
62math.max = Max
63math.round = np.round
64math.ceil = np.ceil
65math.floor = np.floor
66math.acosh=np.arccosh #numpy agrees to 16 digits with std::acosh, math.acosh not
67
68
69scalarTests = True
70vectorTests = True
71debug = False
72caseSymbolic = 0 #exudyn.symbolic
73casePython = 1 #Python
74
75listFunctions = ['abs', 'sign', 'Not',
76 'round', 'ceil', 'floor', 'exp',
77 'sqrt', 'log',
78 'sin', 'cos', 'tan', 'asin', 'acos', 'atan',
79 'sinh', 'cosh', 'tanh', 'asinh', 'acosh', 'atanh',]
80
81listBinFunctions = ['pow', 'atan2', 'mod', 'min', 'max']
82
83
84# some integers, zeros, Reals, etc.
85listNumbers = [-3,0,2,-0.13157932734543,0.,0.2345473645342,math.pi]
86# listNumbers = [-0.13157932734543,0.,math.pi]
87listResults = []
88currentResult = [] #0=symbolic, 1=math/numpy
89cntWrong = 0
90cntTests = 0
91sumResults = 0.
92if scalarTests:
93 #for recording in []:
94 for recording in [True,False]:
95 # for recording in [True]:
96 if debug: exu.Print('***********')
97 if recording:
98 if debug: exu.Print('WITH RECORDING')
99 else:
100 if debug: exu.Print('WITHOUT RECORDING')
101 for number in listNumbers:
102 if debug: exu.Print('NUMBER=',number)
103
104 esym.SetRecording(recording)
105
106 result = [[],[]]
107 for case in [0,1]:
108 if case == 0:
109 Real = SymReal
110 NamedReal = SymReal
111 lib = esym
112 else:
113 Real = float
114 NamedReal = PyRealName
115 lib = math
116
117 a = Real(number)
118 b = NamedReal("b",math.pi)
119 c = NamedReal("c",-0.1)
120 result[case] += [Evaluate(a)]
121 result[case] += [Evaluate(b)]
122 result[case] += [Evaluate(c)]
123
124 d = a+lib.sin(a/b)**2+lib.tan(a)*lib.cos(a)*c*3
125 result[case] += [Evaluate(d)]
126 result[case] += [Evaluate(a+b)]
127 result[case] += [Evaluate(b+a)]
128 result[case] += [Evaluate(b-a)]
129 result[case] += [Evaluate(a-b)]
130 result[case] += [Evaluate(b*a)]
131 result[case] += [Evaluate(a*b)]
132 result[case] += [Evaluate(a/b)]
133 result[case] += [Evaluate(a/3)]
134 result[case] += [Evaluate(-a)]
135 result[case] += [Evaluate(+a)]
136 result[case] += [Evaluate(b**a)]
137 # result[case] += [Evaluate(1**a)] #does not WORK!
138
139 # +=,-=
140 d = Real(0)
141 result[case] += [Evaluate(d)]
142 d += a
143 result[case] += [Evaluate(d)]
144 #Real delete count wrong from here:
145 d -= a*b
146 result[case] += [Evaluate(d)]
147 d *= a
148 result[case] += [Evaluate(d)]
149 if (float(a) != 0.):
150 d /= a
151 result[case] += [Evaluate(d)]
152
153 f = a < b
154 result[case] += [Evaluate(f)]
155 f = a <= b
156 result[case] += [Evaluate(f)]
157 f = a > b
158 result[case] += [Evaluate(f)]
159 f = a >= b
160 result[case] += [Evaluate(f)]
161 f = a == b
162 result[case] += [Evaluate(f)]
163 f = a != b
164 result[case] += [Evaluate(f)]
165
166 #functions
167 for funcStr in listFunctions:
168 if funcStr=='log':
169 if number <= 0: continue
170 if funcStr=='sqrt':
171 if number < 0: continue
172 if funcStr=='acosh':
173 if number < 1 : continue
174 if funcStr=='atanh':
175 if abs(number) >= 1 : continue
176 if funcStr=='acos' or funcStr=='asin':
177 if abs(number) > 1 : continue
178
179 func = getattr(lib, funcStr)
180 # if debug: exu.Print(funcStr)
181 d = func(a)
182 if debug: exu.Print(' '+funcStr+'('+str(Evaluate(a))+')=',Evaluate(d))
183 result[case] += [Evaluate(d)]
184
185 #binary functions
186 for funcStr in listBinFunctions:
187 func = getattr(lib, funcStr)
188 #if debug: exu.Print(funcStr)
189 d = func(a,2)
190 if debug: exu.Print(' '+funcStr+'('+str(Evaluate(a))+','+str(Evaluate(b))+')=',Evaluate(d))
191 result[case] += [Evaluate(d)]
192
193 #ifthenelse
194 if case==0:
195 f = lib.IfThenElse(a, a+1, b+a)
196 else:
197 f = a+1 if a else b+a
198 result[case] += [Evaluate(f)]
199
200
201 result = np.array(result).T.tolist()
202 cntTests += len(result)
203
204 for res in result:
205 s = 'correct'
206 sumResults += res[0]
207 #a RELATIVE tolerance, for the same reason as in the vector section below (#2509,
208 #the successor of #2501): this compares two independent implementations of a
209 #transcendental function - the symbolic module's and Python's - and exact equality
210 #between them holds only by luck. It held until numpy 2.5.3, where acosh(2) comes
211 #out 1 ulp apart (1.3169578969248166 against ...68) and the model reported two
212 #wrong results. What this test is about is that the symbolic module computes the
213 #right function, not that two libraries round identically.
214 if abs(res[0] - res[1]) > 1e-14*max(abs(res[0]), abs(res[1]), 1.):
215 s = 'WRONG:'
216 s += str(res[0]-res[1])
217 cntWrong+=1
218 if debug: exu.Print('. res sym:',res[0], ', res Py:', res[1], s)
219
220#%%++++++++++++++++++++++++++++++++++++++++++++++++++++
221#Vector-Matrix
222SymVector = esym.Vector
223SymMatrix = esym.Matrix
224
225if vectorTests:
226 for recording in [False,True]:
227 # for recording in [True]:
228 if debug: exu.Print('***********')
229 if recording:
230 if debug: exu.Print('WITH RECORDING')
231 else:
232 if debug: exu.Print('WITHOUT RECORDING')
233
234 esym.SetRecording(recording)
235
236 result = [[],[]]
237 for case in [1,0]:
238 if debug: exu.Print('*********\ncase:',case)
239 if case == 0:
240 Real = SymReal
241 NamedReal = SymReal
242 lib = esym
243 Vector = esym.Vector
244 NamedVector = esym.Vector
245 Matrix = esym.Matrix
246 else:
247 Real = float
248 NamedReal = PyRealName
249 lib = math
250 Vector = np.array
251 NamedVector = PyNamedVector
252 Matrix = np.array
253
254 a = NamedReal("a",0.5)
255 b = NamedReal("b",math.pi)
256
257 v = Vector([1.3])
258 w = Vector([4.2])
259 result[case] += [Evaluate(v+w)]
260
261 v = NamedVector('vv',[1,2])
262 w = Vector([4.2-1,-1.2-1])
263 x = Vector([0.,0.])
264
265 if case == 0:
266 w.SetVector([4.2,-1.2])
267 d = 3.3*v+a*w+(v*w)*v
268 x[0] += d.NormL2()
269 x[1] += Evaluate(v == v )
270 x[1] += Evaluate(v == Vector([1,2.1]) )
271 x[1] += Evaluate(v != Vector([1,2.1]) )
272 x[1] += Evaluate(v == Vector([1,2.]) )
273 x[1] += Evaluate(v != Vector([1,2.]) )
274 else:
275 w = np.array([4.2,-1.2])
276 d = 3.3*v+a*w+(v@w)*v
277 x[0] += np.linalg.norm(d)
278 x[1] += (v == v ).all()
279 x[1] += (v == Vector([1,2.1]) ).all()
280 x[1] += (v != Vector([1,2.1]) ).any()
281 x[1] += (v == Vector([1,2.]) ).all()
282 x[1] += (v != Vector([1,2.]) ).any()
283
284 result[case] += [Evaluate(x)]
285
286 result[case] += [Evaluate(d)]
287 result[case] += [Evaluate(v+w)]
288 result[case] += [Evaluate(v-w)]
289 if case==0:
290 n = w.NumberOfItems()
291 result[case] += [Evaluate(v*w)]
292 # print('v0:',x)
293 if case==1:
294 n = len(w)
295 result[case] += [Evaluate(v@w)]
296 # print('v1:',x)
297 x = Vector([n,1.1])
298
299 result[case] += [Evaluate(x)]
300 result[case] += [Evaluate(v)]
301 result[case] += [Evaluate(-v)]
302 result[case] += [Evaluate(a*v)]
303 result[case] += [Evaluate(v*a)]
304 result[case] += [Evaluate(3.3*v)]
305 result[case] += [Evaluate(v*3.3)]
306
307 v = Vector([-0.33,0.347,1.5])
308 w = Vector([4.2,-1.2+10,7.7])
309 w[1] = -1.2
310 result[case] += [Evaluate(w)]
311
312 result[case] += [Evaluate(d)]
313 result[case] += [Evaluate(v+w)]
314 result[case] += [Evaluate(v-w)]
315
316 if case==0:
317 result[case] += [Evaluate(v*w)]
318 result[case] += [Evaluate(v.MultComponents(w))]
319
320 if case==1:
321 result[case] += [Evaluate(v@w)]
322 result[case] += [Evaluate(v*w)]
323
324 result[case] += [Evaluate(v)]
325 result[case] += [Evaluate(-v)]
326 result[case] += [Evaluate(a*v)]
327 result[case] += [Evaluate(v*a)]
328 result[case] += [Evaluate(3.3*v)]
329 result[case] += [Evaluate(v*3.3)]
330
331 u = Vector([0.,0.,0.])
332 u += 2*v
333 # print('case'+str(case)+': 2*v=',Evaluate(2*v),', v=',v,', u=',u)
334 u = u + 2*v
335 result[case] += [Evaluate(u)]
336 u -= v
337 result[case] += [Evaluate(u)]
338 result[case] += [Evaluate(u==v)]
339 result[case] += [Evaluate(u!=v)]
340 u *= a
341 result[case] += [Evaluate(u)]
342 u *= 1/3
343 result[case] += [Evaluate(u)]
344 result[case] += [Evaluate(u==v)]
345 result[case] += [Evaluate(u!=v)]
346
347 #MATRIX MATRIX MATRIX
348 M = Matrix(np.array([[2.,0.1,0.33],
349 [-0.1,2.3,0.7],
350 [0,0.34,1.8]]))
351 N = Matrix(np.array([[1.,0.3,-0.33],
352 [-0.9,1.3,-0.7],
353 [0,0.64,-1.8]]))
354 result[case] += [Evaluate(N)]
355
356 if case==0:
357 M02 = M[0,2].Evaluate()
358 M20 = M[2,0].Evaluate()
359
360 if (M02 != M.Get(0,2).Evaluate()):
361 M02 = -1000
362 #print('*****************ERROR')
363 nr = M.NumberOfRows()
364 nc = M.NumberOfColumns()
365
366 N.SetMatrix(np.array([[(nr+nc)/6,0.3,-M02],
367 [-0.9,1.3,-0.7],
368 [M20,0.64+100,-1.8]]))
369 N[2,1] = 0.64
370
371 result[case] += [Evaluate(N)]
372
373 result[case] += [Evaluate(M)]
374 result[case] += [Evaluate(a*M)]
375 result[case] += [Evaluate(0.5*M)]
376 result[case] += [Evaluate(M*a)]
377 result[case] += [Evaluate(M*0.5)]
378 result[case] += [Evaluate(M+N)]
379 result[case] += [Evaluate(M-N)]
380 if case==0:
381 result[case] += [Evaluate(M*N)]
382 result[case] += [Evaluate(M*v)]
383 result[case] += [Evaluate(v*M)]
384 if case==1:
385 result[case] += [Evaluate(M@N)]
386 result[case] += [Evaluate(M@v)]
387 result[case] += [Evaluate(v@M)]
388
389 P = 0.*M
390 result[case] += [Evaluate(P)]
391
392 #matrix delete counts wrong starting here:
393 P += N
394 result[case] += [Evaluate(P)]
395 P -= M
396 result[case] += [Evaluate(P)]
397 P *= 1.3
398 result[case] += [Evaluate(P)]
399
400
401 # result = np.array(result).T.tolist() #doesnt work for mixed components
402 cntTests += len(result[0])
403 for i in range(len(result[0])):
404 res = [result[0][i],result[1][i]]
405 s = 'correct'
406 sumResults += np.linalg.norm(res[0])
407 #if (res[0] != res[1]).any(): #problem with 1e-16 errors
408 wrong = False
409 #the tolerance is RELATIVE to the magnitude being compared; with the absolute 1e-15 that
410 #stood here, a value of magnitude 9.7476 was compared to within 1e-15 - less than one
411 #ulp (2*eps*9.75 = 4.3e-15), so whether the test passed depended on the last bit of a
412 #numpy sum. It did pass under numpy 2.4.6 and failed under 2.2.4 with a difference of
413 #1.78e-15, on a byte-identical exudyn binary (#2501). A real symbolic/numeric mismatch
414 #is O(1); the few ulp of a different summation order are not what this test is about.
415 scale = max(np.linalg.norm(res[0]), np.linalg.norm(res[1]), 1.)
416 if np.linalg.norm(res[0] - res[1]) > 1e-14*scale:
417 s = '\\diff:\n'
418 s += str(res[0]-res[1])
419 cntWrong+=1
420 wrong=True
421 if debug or wrong:
422 exu.Print('. res sym:\n',res[0], ',\n res Py:\n', res[1], s)
423
424#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
425#checks (#2479): everything above compares NUMBERS
426#against Python's math module. These check what nothing checked at all - Diff, the VariableSet,
427#named variables without recording, __str__ and the error paths. They are written as a FUNCTION so
428#that every symbolic object they create is released on return, which keeps the new/delete balance
429#below intact; and they only count into cntWrong, so that the test result stays what it was.
430def ExtendedSymbolicChecks():
431 wrong = 0
432
433 def Check(condition, info):
434 if not condition:
435 exu.Print('*** wrong: ' + info)
436 return 0 if condition else 1
437
438 #--- Diff: compared against a central difference of the same expression
439 esym.SetRecording(True)
440 for value in [0.3, 1.7, -0.8]:
441 x = SymReal('x', value)
442 cases = [(x*x*x, 3*value**2),
443 (esym.sin(x)*esym.cos(x), math.cos(2*value)),
444 (esym.exp(x)/x, math.exp(value)/value - math.exp(value)/value**2),
445 (esym.sqrt(x*x+1), value/math.sqrt(value**2+1)),
446 ]
447 for (f, analytic) in cases:
448 wrong += Check(abs(f.Diff(x) - analytic) < 1e-12, 'Diff value ' + str(value))
449
450 #pow is differentiated as exp(exponent*log(base)), so its derivative is NaN for base<=0 -
451 #a deliberate answer, not a failure, and nothing said so before
452 f = esym.pow(x, 4)
453 if value > 0:
454 wrong += Check(abs(f.Diff(x) - 4*value**3) < 1e-12, 'Diff pow value ' + str(value))
455 else:
456 wrong += Check(math.isnan(f.Diff(x)), 'Diff pow of a negative base is NaN')
457
458 #the derivative follows the variable, just as the value does
459 f = esym.sin(x*x)
460 x.SetValue(0.5)
461 wrong += Check(abs(f.Diff(x) - math.cos(0.25)*1.) < 1e-12, 'Diff after SetValue')
462
463 #a variable the expression does not contain contributes nothing
464 y = SymReal('y', value)
465 wrong += Check(f.Diff(y) == 0., 'Diff of an unrelated variable')
466
467 #--- the VariableSet, which no test touched so far
468 variables = esym.VariableSet()
469 wrong += Check(variables.NumberOfItems() == 0, 'VariableSet starts empty')
470 variables.Add('alpha', 1.25)
471 variables.Add('beta', -2.)
472 wrong += Check(variables.NumberOfItems() == 2, 'VariableSet counts its variables')
473 wrong += Check(variables.Exists('alpha'), 'VariableSet knows alpha')
474 wrong += Check(not variables.Exists('gamma'), 'VariableSet does not invent gamma')
475 wrong += Check(variables.Get('alpha').Evaluate() == 1.25, 'VariableSet returns the value')
476 variables.Set('alpha', 3.5)
477 wrong += Check(variables.Get('alpha').Evaluate() == 3.5, 'VariableSet overwrites a value')
478 wrong += Check(sorted(variables.GetNames()) == ['alpha', 'beta'], 'VariableSet lists its names')
479 variables.Reset()
480 wrong += Check(variables.NumberOfItems() == 0, 'VariableSet is empty after Reset')
481
482 #--- with recording OFF a named variable is a plain number: no tree, and no link
483 esym.SetRecording(False)
484 x = SymReal('x', 2.)
485 f = x*x
486 wrong += Check(f.Evaluate() == 4., 'unrecorded expression is evaluated immediately')
487 x.SetValue(5.)
488 wrong += Check(f.Evaluate() == 4., 'an unrecorded expression does NOT follow its variable')
489 wrong += Check(float(str(f)) == 4., 'an unrecorded expression prints its value')
490 esym.SetRecording(True)
491
492 #--- __str__ prints the expression tree
493 x = SymReal('x', 2.)
494 wrong += Check(str(x) == 'x', 'a named variable prints its name')
495 wrong += Check(str(x+x) == '(x + x)', 'a sum prints both operands')
496 wrong += Check(str(esym.sqrt(x)) == 'sqrt(x)', 'a function prints its name')
497
498 #--- the error paths. A failed operation used to leak its expression nodes, which showed up in
499 #the new/delete balance below; since #2481 it frees them, so nothing has to be compensated here.
500 #ONLY the regular module: exudynCPPfast is compiled with __FAST_EXUDYN_LINALG, which removes
501 #exactly these range checks by design, so there the two cases below do NOT raise - that is the
502 #trade the fast module makes, not a defect (found when the suite first ran against it, #2495)
503 hasRangeChecks = '[FAST]' not in exu.config.Version(True)
504 if hasRangeChecks:
505 try:
506 esym.Vector([1., 2.]) * esym.Vector([1., 2., 3.])
507 wrong += Check(False, 'a vector product of different sizes must raise')
508 except Exception:
509 pass
510
511 try:
512 v = esym.Vector([1., 2.])
513 v[5]
514 wrong += Check(False, 'an index beyond the vector must raise')
515 except Exception:
516 pass
517
518 return wrong
519
520cntWrong += ExtendedSymbolicChecks()
521
522#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
523#cleanup and check new/delete
524if True:
525 del a,b,c,d,f
526 del u,v,w,x
527 del M,N,P
528
529 #check sum of new/delete
530 newDelSum = 0
531 newDelSum += esym.Real.__newCount-esym.Real.__deleteCount
532 newDelSum += esym.Vector.__newCount-esym.Vector.__deleteCount
533 newDelSum += esym.Matrix.__newCount-esym.Matrix.__deleteCount
534 if newDelSum != 0:
535 exu.Print('*******************')
536 exu.Print('* WARNING *')
537 exu.Print('New-Del=',newDelSum)
538 exu.Print('*******************')
539 sumResults += newDelSum
540
541exu.Print('Real.newCount=',esym.Real.__newCount)
542exu.Print('Real.deleteCount=',esym.Real.__deleteCount)
543
544exu.Print('Vector.newCount=',esym.Vector.__newCount)
545exu.Print('Vector.deleteCount=',esym.Vector.__deleteCount)
546
547exu.Print('Matrix.newCount=',esym.Matrix.__newCount)
548exu.Print('Matrix.deleteCount=',esym.Matrix.__deleteCount)
549
550#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
551
552exu.Print('\nfinished',cntTests,'tests')
553exu.Print('WRONG results (after vectorTests):',cntWrong,'\n')
554
555#a wrong comparison used to move the 14th digit and pass; counting it into the result makes
556#the test fail instead. cntWrong is 0 in a passing run, so the reference is unchanged (#2479)
557u = sumResults/1000 + cntWrong
558exu.Print('u=',u)
559exu.Print('solution of symbolicModuleTest=',u)
560
561# result for 10000 steps; identical for both UF cases
562exu.sys['testResult'] = u