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 #only import if it does not conflict
15
16useGraphics = False #without test
17#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
18#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
19try: #only if called from test suite
20 from modelUnitTests import exudynTestGlobals #for globally storing test results
21 useGraphics = exudynTestGlobals.useGraphics
22except:
23 class ExudynTestGlobals:
24 pass
25 exudynTestGlobals = ExudynTestGlobals()
26#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
27
28esym = exu.symbolic
29import numpy as np
30import math
31
32#reset counters:
33esym.Real.__newCount = 0
34esym.Real.__deleteCount = 0
35esym.Vector.__newCount = 0
36esym.Vector.__deleteCount = 0
37esym.Matrix.__newCount = 0
38esym.Matrix.__deleteCount = 0
39
40SymReal = esym.Real
41
42def Not(value):
43 return not value
44
45def PyRealName(name, value):
46 return value
47
48def PyNamedVector(name, value):
49 return np.array(value)
50
51def Evaluate(value):
52 #print(str(value))
53 if (type(value) == esym.Real
54 or type(value) == esym.Vector
55 or type(value) == esym.Matrix
56 ):
57 return value.Evaluate()
58 if type(value) == np.ndarray:
59 value = value.copy() #requires copy, to store results!
60 return value
61
62def Min(a,b):
63 return a if a<b else b
64def Max(a,b):
65 return a if a>b else b
66
67math.sign = np.sign #to have function in math
68math.Not = Not
69math.abs = np.abs
70math.mod = math.fmod
71math.min = Min
72math.max = Max
73math.round = np.round
74math.ceil = np.ceil
75math.floor = np.floor
76math.acosh=np.arccosh #numpy agrees to 16 digits with std::acosh, math.acosh not
77
78
79scalarTests = True
80vectorTests = True
81debug = False
82caseSymbolic = 0 #exudyn.symbolic
83casePython = 1 #Python
84
85listFunctions = ['abs', 'sign', 'Not',
86 'round', 'ceil', 'floor', 'exp',
87 'sqrt', 'log',
88 'sin', 'cos', 'tan', 'asin', 'acos', 'atan',
89 'sinh', 'cosh', 'tanh', 'asinh', 'acosh', 'atanh',]
90
91listBinFunctions = ['pow', 'atan2', 'mod', 'min', 'max']
92
93
94# some integers, zeros, Reals, etc.
95listNumbers = [-3,0,2,-0.13157932734543,0.,0.2345473645342,math.pi]
96# listNumbers = [-0.13157932734543,0.,math.pi]
97listResults = []
98currentResult = [] #0=symbolic, 1=math/numpy
99cntWrong = 0
100cntTests = 0
101sumResults = 0.
102if scalarTests:
103 #for recording in []:
104 for recording in [True,False]:
105 # for recording in [True]:
106 if debug: exu.Print('***********')
107 if recording:
108 if debug: exu.Print('WITH RECORDING')
109 else:
110 if debug: exu.Print('WITHOUT RECORDING')
111 for number in listNumbers:
112 if debug: exu.Print('NUMBER=',number)
113
114 esym.SetRecording(recording)
115
116 result = [[],[]]
117 for case in [0,1]:
118 if case == 0:
119 Real = SymReal
120 NamedReal = SymReal
121 lib = esym
122 else:
123 Real = float
124 NamedReal = PyRealName
125 lib = math
126
127 a = Real(number)
128 b = NamedReal("b",math.pi)
129 c = NamedReal("c",-0.1)
130 result[case] += [Evaluate(a)]
131 result[case] += [Evaluate(b)]
132 result[case] += [Evaluate(c)]
133
134 d = a+lib.sin(a/b)**2+lib.tan(a)*lib.cos(a)*c*3
135 result[case] += [Evaluate(d)]
136 result[case] += [Evaluate(a+b)]
137 result[case] += [Evaluate(b+a)]
138 result[case] += [Evaluate(b-a)]
139 result[case] += [Evaluate(a-b)]
140 result[case] += [Evaluate(b*a)]
141 result[case] += [Evaluate(a*b)]
142 result[case] += [Evaluate(a/b)]
143 result[case] += [Evaluate(a/3)]
144 result[case] += [Evaluate(-a)]
145 result[case] += [Evaluate(+a)]
146 result[case] += [Evaluate(b**a)]
147 # result[case] += [Evaluate(1**a)] #does not WORK!
148
149 # +=,-=
150 d = Real(0)
151 result[case] += [Evaluate(d)]
152 d += a
153 result[case] += [Evaluate(d)]
154 #Real delete count wrong from here:
155 d -= a*b
156 result[case] += [Evaluate(d)]
157 d *= a
158 result[case] += [Evaluate(d)]
159 if (float(a) != 0.):
160 d /= a
161 result[case] += [Evaluate(d)]
162
163 f = a < b
164 result[case] += [Evaluate(f)]
165 f = a <= b
166 result[case] += [Evaluate(f)]
167 f = a > b
168 result[case] += [Evaluate(f)]
169 f = a >= b
170 result[case] += [Evaluate(f)]
171 f = a == b
172 result[case] += [Evaluate(f)]
173 f = a != b
174 result[case] += [Evaluate(f)]
175
176 #functions
177 for funcStr in listFunctions:
178 if funcStr=='log':
179 if number <= 0: continue
180 if funcStr=='sqrt':
181 if number < 0: continue
182 if funcStr=='acosh':
183 if number < 1 : continue
184 if funcStr=='atanh':
185 if abs(number) >= 1 : continue
186 if funcStr=='acos' or funcStr=='asin':
187 if abs(number) > 1 : continue
188
189 func = getattr(lib, funcStr)
190 # if debug: exu.Print(funcStr)
191 d = func(a)
192 if debug: exu.Print(' '+funcStr+'('+str(Evaluate(a))+')=',Evaluate(d))
193 result[case] += [Evaluate(d)]
194
195 #binary functions
196 for funcStr in listBinFunctions:
197 func = getattr(lib, funcStr)
198 #if debug: exu.Print(funcStr)
199 d = func(a,2)
200 if debug: exu.Print(' '+funcStr+'('+str(Evaluate(a))+','+str(Evaluate(b))+')=',Evaluate(d))
201 result[case] += [Evaluate(d)]
202
203 #ifthenelse
204 if case==0:
205 f = lib.IfThenElse(a, a+1, b+a)
206 else:
207 f = a+1 if a else b+a
208 result[case] += [Evaluate(f)]
209
210
211 result = np.array(result).T.tolist()
212 cntTests += len(result)
213
214 for res in result:
215 s = 'correct'
216 sumResults += res[0]
217 if res[0] != res[1]:
218 s = 'WRONG:'
219 s += str(res[0]-res[1])
220 cntWrong+=1
221 if debug: exu.Print('. res sym:',res[0], ', res Py:', res[1], s)
222
223#%%++++++++++++++++++++++++++++++++++++++++++++++++++++
224#Vector-Matrix
225SymVector = esym.Vector
226SymMatrix = esym.Matrix
227
228if vectorTests:
229 for recording in [False,True]:
230 # for recording in [True]:
231 if debug: exu.Print('***********')
232 if recording:
233 if debug: exu.Print('WITH RECORDING')
234 else:
235 if debug: exu.Print('WITHOUT RECORDING')
236
237 esym.SetRecording(recording)
238
239 result = [[],[]]
240 for case in [1,0]:
241 if debug: exu.Print('*********\ncase:',case)
242 if case == 0:
243 Real = SymReal
244 NamedReal = SymReal
245 lib = esym
246 Vector = esym.Vector
247 NamedVector = esym.Vector
248 Matrix = esym.Matrix
249 else:
250 Real = float
251 NamedReal = PyRealName
252 lib = math
253 Vector = np.array
254 NamedVector = PyNamedVector
255 Matrix = np.array
256
257 a = NamedReal("a",0.5)
258 b = NamedReal("b",math.pi)
259
260 v = Vector([1.3])
261 w = Vector([4.2])
262 result[case] += [Evaluate(v+w)]
263
264 v = NamedVector('vv',[1,2])
265 w = Vector([4.2-1,-1.2-1])
266 x = Vector([0.,0.])
267
268 if case == 0:
269 w.SetVector([4.2,-1.2])
270 d = 3.3*v+a*w+(v*w)*v
271 x[0] += d.NormL2()
272 x[1] += Evaluate(v == v )
273 x[1] += Evaluate(v == Vector([1,2.1]) )
274 x[1] += Evaluate(v != Vector([1,2.1]) )
275 x[1] += Evaluate(v == Vector([1,2.]) )
276 x[1] += Evaluate(v != Vector([1,2.]) )
277 else:
278 w = np.array([4.2,-1.2])
279 d = 3.3*v+a*w+(v@w)*v
280 x[0] += np.linalg.norm(d)
281 x[1] += (v == v ).all()
282 x[1] += (v == Vector([1,2.1]) ).all()
283 x[1] += (v != Vector([1,2.1]) ).any()
284 x[1] += (v == Vector([1,2.]) ).all()
285 x[1] += (v != Vector([1,2.]) ).any()
286
287 result[case] += [Evaluate(x)]
288
289 result[case] += [Evaluate(d)]
290 result[case] += [Evaluate(v+w)]
291 result[case] += [Evaluate(v-w)]
292 if case==0:
293 n = w.NumberOfItems()
294 result[case] += [Evaluate(v*w)]
295 # print('v0:',x)
296 if case==1:
297 n = len(w)
298 result[case] += [Evaluate(v@w)]
299 # print('v1:',x)
300 x = Vector([n,1.1])
301
302 result[case] += [Evaluate(x)]
303 result[case] += [Evaluate(v)]
304 result[case] += [Evaluate(-v)]
305 result[case] += [Evaluate(a*v)]
306 result[case] += [Evaluate(v*a)]
307 result[case] += [Evaluate(3.3*v)]
308 result[case] += [Evaluate(v*3.3)]
309
310 v = Vector([-0.33,0.347,1.5])
311 w = Vector([4.2,-1.2+10,7.7])
312 w[1] = -1.2
313 result[case] += [Evaluate(w)]
314
315 result[case] += [Evaluate(d)]
316 result[case] += [Evaluate(v+w)]
317 result[case] += [Evaluate(v-w)]
318
319 if case==0:
320 result[case] += [Evaluate(v*w)]
321 result[case] += [Evaluate(v.MultComponents(w))]
322
323 if case==1:
324 result[case] += [Evaluate(v@w)]
325 result[case] += [Evaluate(v*w)]
326
327 result[case] += [Evaluate(v)]
328 result[case] += [Evaluate(-v)]
329 result[case] += [Evaluate(a*v)]
330 result[case] += [Evaluate(v*a)]
331 result[case] += [Evaluate(3.3*v)]
332 result[case] += [Evaluate(v*3.3)]
333
334 u = Vector([0.,0.,0.])
335 u += 2*v
336 # print('case'+str(case)+': 2*v=',Evaluate(2*v),', v=',v,', u=',u)
337 u = u + 2*v
338 result[case] += [Evaluate(u)]
339 u -= v
340 result[case] += [Evaluate(u)]
341 result[case] += [Evaluate(u==v)]
342 result[case] += [Evaluate(u!=v)]
343 u *= a
344 result[case] += [Evaluate(u)]
345 u *= 1/3
346 result[case] += [Evaluate(u)]
347 result[case] += [Evaluate(u==v)]
348 result[case] += [Evaluate(u!=v)]
349
350 #MATRIX MATRIX MATRIX
351 M = Matrix(np.array([[2.,0.1,0.33],
352 [-0.1,2.3,0.7],
353 [0,0.34,1.8]]))
354 N = Matrix(np.array([[1.,0.3,-0.33],
355 [-0.9,1.3,-0.7],
356 [0,0.64,-1.8]]))
357 result[case] += [Evaluate(N)]
358
359 if case==0:
360 M02 = M[0,2].Evaluate()
361 M20 = M[2,0].Evaluate()
362
363 if (M02 != M.Get(0,2).Evaluate()):
364 M02 = -1000
365 #print('*****************ERROR')
366 nr = M.NumberOfRows()
367 nc = M.NumberOfColumns()
368
369 N.SetMatrix(np.array([[(nr+nc)/6,0.3,-M02],
370 [-0.9,1.3,-0.7],
371 [M20,0.64+100,-1.8]]))
372 N[2,1] = 0.64
373
374 result[case] += [Evaluate(N)]
375
376 result[case] += [Evaluate(M)]
377 result[case] += [Evaluate(a*M)]
378 result[case] += [Evaluate(0.5*M)]
379 result[case] += [Evaluate(M*a)]
380 result[case] += [Evaluate(M*0.5)]
381 result[case] += [Evaluate(M+N)]
382 result[case] += [Evaluate(M-N)]
383 if case==0:
384 result[case] += [Evaluate(M*N)]
385 result[case] += [Evaluate(M*v)]
386 result[case] += [Evaluate(v*M)]
387 if case==1:
388 result[case] += [Evaluate(M@N)]
389 result[case] += [Evaluate(M@v)]
390 result[case] += [Evaluate(v@M)]
391
392 P = 0.*M
393 result[case] += [Evaluate(P)]
394
395 #matrix delete counts wrong starting here:
396 P += N
397 result[case] += [Evaluate(P)]
398 P -= M
399 result[case] += [Evaluate(P)]
400 P *= 1.3
401 result[case] += [Evaluate(P)]
402
403
404 # result = np.array(result).T.tolist() #doesnt work for mixed components
405 cntTests += len(result[0])
406 for i in range(len(result[0])):
407 res = [result[0][i],result[1][i]]
408 s = 'correct'
409 sumResults += np.linalg.norm(res[0])
410 #if (res[0] != res[1]).any(): #problem with 1e-16 errors
411 wrong = False
412 if np.linalg.norm(res[0] - res[1]) > 1e-15:
413 s = '\\diff:\n'
414 s += str(res[0]-res[1])
415 cntWrong+=1
416 wrong=True
417 if debug or wrong:
418 exu.Print('. res sym:\n',res[0], ',\n res Py:\n', res[1], s)
419
420#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
421#cleanup and check new/delete
422if True:
423 del a,b,c,d,f
424 del u,v,w,x
425 del M,N,P
426
427 #check sum of new/delete
428 newDelSum = 0
429 newDelSum += esym.Real.__newCount-esym.Real.__deleteCount
430 newDelSum += esym.Vector.__newCount-esym.Vector.__deleteCount
431 newDelSum += esym.Matrix.__newCount-esym.Matrix.__deleteCount
432 if newDelSum != 0:
433 exu.Print('*******************')
434 exu.Print('* WARNING *')
435 exu.Print('New-Del=',newDelSum)
436 exu.Print('*******************')
437 sumResults += newDelSum
438
439exu.Print('Real.newCount=',esym.Real.__newCount)
440exu.Print('Real.deleteCount=',esym.Real.__deleteCount)
441
442exu.Print('Vector.newCount=',esym.Vector.__newCount)
443exu.Print('Vector.deleteCount=',esym.Vector.__deleteCount)
444
445exu.Print('Matrix.newCount=',esym.Matrix.__newCount)
446exu.Print('Matrix.deleteCount=',esym.Matrix.__deleteCount)
447
448#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
449
450exu.Print('\nfinished',cntTests,'tests')
451exu.Print('WRONG results (after vectorTests):',cntWrong,'\n')
452
453#
454u = sumResults/1000
455exu.Print('u=',u)
456exu.Print('solution of symbolicModuleTest=',u)
457
458# result for 10000 steps; identical for both UF cases
459exudynTestGlobals.testError = u - (0.9480053738744615)
460exudynTestGlobals.testResult = u