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