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