(testmodels-symbolicmoduletest)= # symbolicModuleTest.py You can view and download this file on Github: [symbolicModuleTest.py](https://github.com/jgerstmayr/EXUDYN/tree/master/python/TestModels/symbolicModuleTest.py) ```{code-block} python :linenos: #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # This is an EXUDYN example # # Details: Tests for symbolic scalar, vector and matrix # # Author: Johannes Gerstmayr # Date: 2023-12-14 # # 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. # #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ import exudyn as exu from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities import exudyn.graphics as graphics testIsActive = exu.sys.get('testIsActive', True) esym = exu.symbolic import numpy as np import math #reset counters: esym.Real.__newCount = 0 esym.Real.__deleteCount = 0 esym.Vector.__newCount = 0 esym.Vector.__deleteCount = 0 esym.Matrix.__newCount = 0 esym.Matrix.__deleteCount = 0 SymReal = esym.Real def Not(value): return not value def PyRealName(name, value): return value def PyNamedVector(name, value): return np.array(value) def Evaluate(value): #print(str(value)) if (type(value) == esym.Real or type(value) == esym.Vector or type(value) == esym.Matrix ): return value.Evaluate() if type(value) == np.ndarray: value = value.copy() #requires copy, to store results! return value def Min(a,b): return a if ab else b math.sign = np.sign #to have function in math math.Not = Not math.abs = np.abs math.mod = math.fmod math.min = Min math.max = Max math.round = np.round math.ceil = np.ceil math.floor = np.floor math.acosh=np.arccosh #numpy agrees to 16 digits with std::acosh, math.acosh not scalarTests = True vectorTests = True debug = False caseSymbolic = 0 #exudyn.symbolic casePython = 1 #Python listFunctions = ['abs', 'sign', 'Not', 'round', 'ceil', 'floor', 'exp', 'sqrt', 'log', 'sin', 'cos', 'tan', 'asin', 'acos', 'atan', 'sinh', 'cosh', 'tanh', 'asinh', 'acosh', 'atanh',] listBinFunctions = ['pow', 'atan2', 'mod', 'min', 'max'] # some integers, zeros, Reals, etc. listNumbers = [-3,0,2,-0.13157932734543,0.,0.2345473645342,math.pi] # listNumbers = [-0.13157932734543,0.,math.pi] listResults = [] currentResult = [] #0=symbolic, 1=math/numpy cntWrong = 0 cntTests = 0 sumResults = 0. if scalarTests: #for recording in []: for recording in [True,False]: # for recording in [True]: if debug: exu.Print('***********') if recording: if debug: exu.Print('WITH RECORDING') else: if debug: exu.Print('WITHOUT RECORDING') for number in listNumbers: if debug: exu.Print('NUMBER=',number) esym.SetRecording(recording) result = [[],[]] for case in [0,1]: if case == 0: Real = SymReal NamedReal = SymReal lib = esym else: Real = float NamedReal = PyRealName lib = math a = Real(number) b = NamedReal("b",math.pi) c = NamedReal("c",-0.1) result[case] += [Evaluate(a)] result[case] += [Evaluate(b)] result[case] += [Evaluate(c)] d = a+lib.sin(a/b)**2+lib.tan(a)*lib.cos(a)*c*3 result[case] += [Evaluate(d)] result[case] += [Evaluate(a+b)] result[case] += [Evaluate(b+a)] result[case] += [Evaluate(b-a)] result[case] += [Evaluate(a-b)] result[case] += [Evaluate(b*a)] result[case] += [Evaluate(a*b)] result[case] += [Evaluate(a/b)] result[case] += [Evaluate(a/3)] result[case] += [Evaluate(-a)] result[case] += [Evaluate(+a)] result[case] += [Evaluate(b**a)] # result[case] += [Evaluate(1**a)] #does not WORK! # +=,-= d = Real(0) result[case] += [Evaluate(d)] d += a result[case] += [Evaluate(d)] #Real delete count wrong from here: d -= a*b result[case] += [Evaluate(d)] d *= a result[case] += [Evaluate(d)] if (float(a) != 0.): d /= a result[case] += [Evaluate(d)] f = a < b result[case] += [Evaluate(f)] f = a <= b result[case] += [Evaluate(f)] f = a > b result[case] += [Evaluate(f)] f = a >= b result[case] += [Evaluate(f)] f = a == b result[case] += [Evaluate(f)] f = a != b result[case] += [Evaluate(f)] #functions for funcStr in listFunctions: if funcStr=='log': if number <= 0: continue if funcStr=='sqrt': if number < 0: continue if funcStr=='acosh': if number < 1 : continue if funcStr=='atanh': if abs(number) >= 1 : continue if funcStr=='acos' or funcStr=='asin': if abs(number) > 1 : continue func = getattr(lib, funcStr) # if debug: exu.Print(funcStr) d = func(a) if debug: exu.Print(' '+funcStr+'('+str(Evaluate(a))+')=',Evaluate(d)) result[case] += [Evaluate(d)] #binary functions for funcStr in listBinFunctions: func = getattr(lib, funcStr) #if debug: exu.Print(funcStr) d = func(a,2) if debug: exu.Print(' '+funcStr+'('+str(Evaluate(a))+','+str(Evaluate(b))+')=',Evaluate(d)) result[case] += [Evaluate(d)] #ifthenelse if case==0: f = lib.IfThenElse(a, a+1, b+a) else: f = a+1 if a else b+a result[case] += [Evaluate(f)] result = np.array(result).T.tolist() cntTests += len(result) for res in result: s = 'correct' sumResults += res[0] #a RELATIVE tolerance, for the same reason as in the vector section below (#2509, #the successor of #2501): this compares two independent implementations of a #transcendental function - the symbolic module's and Python's - and exact equality #between them holds only by luck. It held until numpy 2.5.3, where acosh(2) comes #out 1 ulp apart (1.3169578969248166 against ...68) and the model reported two #wrong results. What this test is about is that the symbolic module computes the #right function, not that two libraries round identically. if abs(res[0] - res[1]) > 1e-14*max(abs(res[0]), abs(res[1]), 1.): s = 'WRONG:' s += str(res[0]-res[1]) cntWrong+=1 if debug: exu.Print('. res sym:',res[0], ', res Py:', res[1], s) #%%++++++++++++++++++++++++++++++++++++++++++++++++++++ #Vector-Matrix SymVector = esym.Vector SymMatrix = esym.Matrix if vectorTests: for recording in [False,True]: # for recording in [True]: if debug: exu.Print('***********') if recording: if debug: exu.Print('WITH RECORDING') else: if debug: exu.Print('WITHOUT RECORDING') esym.SetRecording(recording) result = [[],[]] for case in [1,0]: if debug: exu.Print('*********\ncase:',case) if case == 0: Real = SymReal NamedReal = SymReal lib = esym Vector = esym.Vector NamedVector = esym.Vector Matrix = esym.Matrix else: Real = float NamedReal = PyRealName lib = math Vector = np.array NamedVector = PyNamedVector Matrix = np.array a = NamedReal("a",0.5) b = NamedReal("b",math.pi) v = Vector([1.3]) w = Vector([4.2]) result[case] += [Evaluate(v+w)] v = NamedVector('vv',[1,2]) w = Vector([4.2-1,-1.2-1]) x = Vector([0.,0.]) if case == 0: w.SetVector([4.2,-1.2]) d = 3.3*v+a*w+(v*w)*v x[0] += d.NormL2() x[1] += Evaluate(v == v ) x[1] += Evaluate(v == Vector([1,2.1]) ) x[1] += Evaluate(v != Vector([1,2.1]) ) x[1] += Evaluate(v == Vector([1,2.]) ) x[1] += Evaluate(v != Vector([1,2.]) ) else: w = np.array([4.2,-1.2]) d = 3.3*v+a*w+(v@w)*v x[0] += np.linalg.norm(d) x[1] += (v == v ).all() x[1] += (v == Vector([1,2.1]) ).all() x[1] += (v != Vector([1,2.1]) ).any() x[1] += (v == Vector([1,2.]) ).all() x[1] += (v != Vector([1,2.]) ).any() result[case] += [Evaluate(x)] result[case] += [Evaluate(d)] result[case] += [Evaluate(v+w)] result[case] += [Evaluate(v-w)] if case==0: n = w.NumberOfItems() result[case] += [Evaluate(v*w)] # print('v0:',x) if case==1: n = len(w) result[case] += [Evaluate(v@w)] # print('v1:',x) x = Vector([n,1.1]) result[case] += [Evaluate(x)] result[case] += [Evaluate(v)] result[case] += [Evaluate(-v)] result[case] += [Evaluate(a*v)] result[case] += [Evaluate(v*a)] result[case] += [Evaluate(3.3*v)] result[case] += [Evaluate(v*3.3)] v = Vector([-0.33,0.347,1.5]) w = Vector([4.2,-1.2+10,7.7]) w[1] = -1.2 result[case] += [Evaluate(w)] result[case] += [Evaluate(d)] result[case] += [Evaluate(v+w)] result[case] += [Evaluate(v-w)] if case==0: result[case] += [Evaluate(v*w)] result[case] += [Evaluate(v.MultComponents(w))] if case==1: result[case] += [Evaluate(v@w)] result[case] += [Evaluate(v*w)] result[case] += [Evaluate(v)] result[case] += [Evaluate(-v)] result[case] += [Evaluate(a*v)] result[case] += [Evaluate(v*a)] result[case] += [Evaluate(3.3*v)] result[case] += [Evaluate(v*3.3)] u = Vector([0.,0.,0.]) u += 2*v # print('case'+str(case)+': 2*v=',Evaluate(2*v),', v=',v,', u=',u) u = u + 2*v result[case] += [Evaluate(u)] u -= v result[case] += [Evaluate(u)] result[case] += [Evaluate(u==v)] result[case] += [Evaluate(u!=v)] u *= a result[case] += [Evaluate(u)] u *= 1/3 result[case] += [Evaluate(u)] result[case] += [Evaluate(u==v)] result[case] += [Evaluate(u!=v)] #MATRIX MATRIX MATRIX M = Matrix(np.array([[2.,0.1,0.33], [-0.1,2.3,0.7], [0,0.34,1.8]])) N = Matrix(np.array([[1.,0.3,-0.33], [-0.9,1.3,-0.7], [0,0.64,-1.8]])) result[case] += [Evaluate(N)] if case==0: M02 = M[0,2].Evaluate() M20 = M[2,0].Evaluate() if (M02 != M.Get(0,2).Evaluate()): M02 = -1000 #print('*****************ERROR') nr = M.NumberOfRows() nc = M.NumberOfColumns() N.SetMatrix(np.array([[(nr+nc)/6,0.3,-M02], [-0.9,1.3,-0.7], [M20,0.64+100,-1.8]])) N[2,1] = 0.64 result[case] += [Evaluate(N)] result[case] += [Evaluate(M)] result[case] += [Evaluate(a*M)] result[case] += [Evaluate(0.5*M)] result[case] += [Evaluate(M*a)] result[case] += [Evaluate(M*0.5)] result[case] += [Evaluate(M+N)] result[case] += [Evaluate(M-N)] if case==0: result[case] += [Evaluate(M*N)] result[case] += [Evaluate(M*v)] result[case] += [Evaluate(v*M)] if case==1: result[case] += [Evaluate(M@N)] result[case] += [Evaluate(M@v)] result[case] += [Evaluate(v@M)] P = 0.*M result[case] += [Evaluate(P)] #matrix delete counts wrong starting here: P += N result[case] += [Evaluate(P)] P -= M result[case] += [Evaluate(P)] P *= 1.3 result[case] += [Evaluate(P)] # result = np.array(result).T.tolist() #doesnt work for mixed components cntTests += len(result[0]) for i in range(len(result[0])): res = [result[0][i],result[1][i]] s = 'correct' sumResults += np.linalg.norm(res[0]) #if (res[0] != res[1]).any(): #problem with 1e-16 errors wrong = False #the tolerance is RELATIVE to the magnitude being compared; with the absolute 1e-15 that #stood here, a value of magnitude 9.7476 was compared to within 1e-15 - less than one #ulp (2*eps*9.75 = 4.3e-15), so whether the test passed depended on the last bit of a #numpy sum. It did pass under numpy 2.4.6 and failed under 2.2.4 with a difference of #1.78e-15, on a byte-identical exudyn binary (#2501). A real symbolic/numeric mismatch #is O(1); the few ulp of a different summation order are not what this test is about. scale = max(np.linalg.norm(res[0]), np.linalg.norm(res[1]), 1.) if np.linalg.norm(res[0] - res[1]) > 1e-14*scale: s = '\\diff:\n' s += str(res[0]-res[1]) cntWrong+=1 wrong=True if debug or wrong: exu.Print('. res sym:\n',res[0], ',\n res Py:\n', res[1], s) #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #checks (#2479): everything above compares NUMBERS #against Python's math module. These check what nothing checked at all - Diff, the VariableSet, #named variables without recording, __str__ and the error paths. They are written as a FUNCTION so #that every symbolic object they create is released on return, which keeps the new/delete balance #below intact; and they only count into cntWrong, so that the test result stays what it was. def ExtendedSymbolicChecks(): wrong = 0 def Check(condition, info): if not condition: exu.Print('*** wrong: ' + info) return 0 if condition else 1 #--- Diff: compared against a central difference of the same expression esym.SetRecording(True) for value in [0.3, 1.7, -0.8]: x = SymReal('x', value) cases = [(x*x*x, 3*value**2), (esym.sin(x)*esym.cos(x), math.cos(2*value)), (esym.exp(x)/x, math.exp(value)/value - math.exp(value)/value**2), (esym.sqrt(x*x+1), value/math.sqrt(value**2+1)), ] for (f, analytic) in cases: wrong += Check(abs(f.Diff(x) - analytic) < 1e-12, 'Diff value ' + str(value)) #pow is differentiated as exp(exponent*log(base)), so its derivative is NaN for base<=0 - #a deliberate answer, not a failure, and nothing said so before f = esym.pow(x, 4) if value > 0: wrong += Check(abs(f.Diff(x) - 4*value**3) < 1e-12, 'Diff pow value ' + str(value)) else: wrong += Check(math.isnan(f.Diff(x)), 'Diff pow of a negative base is NaN') #the derivative follows the variable, just as the value does f = esym.sin(x*x) x.SetValue(0.5) wrong += Check(abs(f.Diff(x) - math.cos(0.25)*1.) < 1e-12, 'Diff after SetValue') #a variable the expression does not contain contributes nothing y = SymReal('y', value) wrong += Check(f.Diff(y) == 0., 'Diff of an unrelated variable') #--- the VariableSet, which no test touched so far variables = esym.VariableSet() wrong += Check(variables.NumberOfItems() == 0, 'VariableSet starts empty') variables.Add('alpha', 1.25) variables.Add('beta', -2.) wrong += Check(variables.NumberOfItems() == 2, 'VariableSet counts its variables') wrong += Check(variables.Exists('alpha'), 'VariableSet knows alpha') wrong += Check(not variables.Exists('gamma'), 'VariableSet does not invent gamma') wrong += Check(variables.Get('alpha').Evaluate() == 1.25, 'VariableSet returns the value') variables.Set('alpha', 3.5) wrong += Check(variables.Get('alpha').Evaluate() == 3.5, 'VariableSet overwrites a value') wrong += Check(sorted(variables.GetNames()) == ['alpha', 'beta'], 'VariableSet lists its names') variables.Reset() wrong += Check(variables.NumberOfItems() == 0, 'VariableSet is empty after Reset') #--- with recording OFF a named variable is a plain number: no tree, and no link esym.SetRecording(False) x = SymReal('x', 2.) f = x*x wrong += Check(f.Evaluate() == 4., 'unrecorded expression is evaluated immediately') x.SetValue(5.) wrong += Check(f.Evaluate() == 4., 'an unrecorded expression does NOT follow its variable') wrong += Check(float(str(f)) == 4., 'an unrecorded expression prints its value') esym.SetRecording(True) #--- __str__ prints the expression tree x = SymReal('x', 2.) wrong += Check(str(x) == 'x', 'a named variable prints its name') wrong += Check(str(x+x) == '(x + x)', 'a sum prints both operands') wrong += Check(str(esym.sqrt(x)) == 'sqrt(x)', 'a function prints its name') #--- the error paths. A failed operation used to leak its expression nodes, which showed up in #the new/delete balance below; since #2481 it frees them, so nothing has to be compensated here. #ONLY the regular module: exudynCPPfast is compiled with __FAST_EXUDYN_LINALG, which removes #exactly these range checks by design, so there the two cases below do NOT raise - that is the #trade the fast module makes, not a defect (found when the suite first ran against it, #2495) hasRangeChecks = '[FAST]' not in exu.config.Version(True) if hasRangeChecks: try: esym.Vector([1., 2.]) * esym.Vector([1., 2., 3.]) wrong += Check(False, 'a vector product of different sizes must raise') except Exception: pass try: v = esym.Vector([1., 2.]) v[5] wrong += Check(False, 'an index beyond the vector must raise') except Exception: pass return wrong cntWrong += ExtendedSymbolicChecks() #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #cleanup and check new/delete if True: del a,b,c,d,f del u,v,w,x del M,N,P #check sum of new/delete newDelSum = 0 newDelSum += esym.Real.__newCount-esym.Real.__deleteCount newDelSum += esym.Vector.__newCount-esym.Vector.__deleteCount newDelSum += esym.Matrix.__newCount-esym.Matrix.__deleteCount if newDelSum != 0: exu.Print('*******************') exu.Print('* WARNING *') exu.Print('New-Del=',newDelSum) exu.Print('*******************') sumResults += newDelSum exu.Print('Real.newCount=',esym.Real.__newCount) exu.Print('Real.deleteCount=',esym.Real.__deleteCount) exu.Print('Vector.newCount=',esym.Vector.__newCount) exu.Print('Vector.deleteCount=',esym.Vector.__deleteCount) exu.Print('Matrix.newCount=',esym.Matrix.__newCount) exu.Print('Matrix.deleteCount=',esym.Matrix.__deleteCount) #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ exu.Print('\nfinished',cntTests,'tests') exu.Print('WRONG results (after vectorTests):',cntWrong,'\n') #a wrong comparison used to move the 14th digit and pass; counting it into the result makes #the test fail instead. cntWrong is 0 in a passing run, so the reference is unchanged (#2479) u = sumResults/1000 + cntWrong exu.Print('u=',u) exu.Print('solution of symbolicModuleTest=',u) # result for 10000 steps; identical for both UF cases exu.sys['testResult'] = u ```