The following are the current version of the grading functions used for grading the quality of the antiderivative with reference to the optimal antiderivative included in the test suite.
There is a version for Maple and for Mathematica/Rubi. There is a version for grading Sympy and version for use with Sagemath.
The following are links to the current source code.
The following are the listings of source code of the grading functions.
(* Original version thanks to Albert Rich emailed on 03/21/2017 *) (* ::Package:: *) (* Nasser: April 7,2022. add second output which gives reason for the grade *) (* Small rewrite of logic in main function to make it*) (* match Maple's logic. No change in functionality otherwise*) (* ::Subsection:: *) (*GradeAntiderivative[result,optimal]*) (* ::Text:: *) (*If result and optimal are mathematical expressions, *) (* GradeAntiderivative[result,optimal] returns*) (* "F" if the result fails to integrate an expression that*) (* is integrable*) (* "C" if result involves higher level functions than necessary*) (* "B" if result is more than twice the size of the optimal*) (* antiderivative*) (* "A" if result can be considered optimal*) GradeAntiderivative[result_,optimal_] := Module[{expnResult,expnOptimal,leafCountResult,leafCountOptimal,finalresult}, expnResult = ExpnType[result]; expnOptimal = ExpnType[optimal]; leafCountResult = LeafCount[result]; leafCountOptimal = LeafCount[optimal]; (*Print["expnResult=",expnResult," expnOptimal=",expnOptimal];*) If[expnResult<=expnOptimal, If[Not[FreeQ[result,Complex]], (*result contains complex*) If[Not[FreeQ[optimal,Complex]], (*optimal contains complex*) If[leafCountResult<=2*leafCountOptimal, finalresult={"A","none"} ,(*ELSE*) finalresult={"B","Both result and optimal contain complex but leaf count is larger than twice the leaf count of optimal. $"<>ToString@leafCountResult<>"$ vs. $2("<>ToString@leafCountOptimal<>")="<>ToString[2*leafCountOptimal]<>"$."} ] ,(*ELSE*) finalresult={"C","Result contains complex when optimal does not."} ] ,(*ELSE*)(*result does not contains complex*) If[leafCountResult<=2*leafCountOptimal, finalresult={"A","none"} ,(*ELSE*) finalresult={"B","Leaf count is larger than twice the leaf count of optimal. $"<>ToString@leafCountResult<>"$ vs. $2("<>ToString@leafCountOptimal<>")="<>ToString[2*leafCountOptimal]<>"$."} ] ] ,(*ELSE*) (*expnResult>expnOptimal*) If[FreeQ[result,Integrate] && FreeQ[result,Int], finalresult={"C","Result contains higher order function than in optimal. Order "<>ToString[expnResult]<>" vs. order "<>ToString[expnOptimal]<>" in optimal."} , finalresult={"F","Contains unresolved integral."} ] ]; finalresult ] (* ::Text:: *) (*The following summarizes the type number assigned an *) (*expression based on the functions it involves*) (*1 = rational function*) (*2 = algebraic function*) (*3 = elementary function*) (*4 = special function*) (*5 = hyperpergeometric function*) (*6 = appell function*) (*7 = rootsum function*) (*8 = integrate function*) (*9 = unknown function*) ExpnType[expn_] := If[AtomQ[expn], 1, If[ListQ[expn], Max[Map[ExpnType,expn]], If[Head[expn]===Power, If[IntegerQ[expn[[2]]], ExpnType[expn[[1]]], If[Head[expn[[2]]]===Rational, If[IntegerQ[expn[[1]]] || Head[expn[[1]]]===Rational, 1, Max[ExpnType[expn[[1]]],2]], Max[ExpnType[expn[[1]]],ExpnType[expn[[2]]],3]]], If[Head[expn]===Plus || Head[expn]===Times, Max[ExpnType[First[expn]],ExpnType[Rest[expn]]], If[ElementaryFunctionQ[Head[expn]], Max[3,ExpnType[expn[[1]]]], If[SpecialFunctionQ[Head[expn]], Apply[Max,Append[Map[ExpnType,Apply[List,expn]],4]], If[HypergeometricFunctionQ[Head[expn]], Apply[Max,Append[Map[ExpnType,Apply[List,expn]],5]], If[AppellFunctionQ[Head[expn]], Apply[Max,Append[Map[ExpnType,Apply[List,expn]],6]], If[Head[expn]===RootSum, Apply[Max,Append[Map[ExpnType,Apply[List,expn]],7]], If[Head[expn]===Integrate || Head[expn]===Int, Apply[Max,Append[Map[ExpnType,Apply[List,expn]],8]], 9]]]]]]]]]] ElementaryFunctionQ[func_] := MemberQ[{ Exp,Log, Sin,Cos,Tan,Cot,Sec,Csc, ArcSin,ArcCos,ArcTan,ArcCot,ArcSec,ArcCsc, Sinh,Cosh,Tanh,Coth,Sech,Csch, ArcSinh,ArcCosh,ArcTanh,ArcCoth,ArcSech,ArcCsch },func] SpecialFunctionQ[func_] := MemberQ[{ Erf, Erfc, Erfi, FresnelS, FresnelC, ExpIntegralE, ExpIntegralEi, LogIntegral, SinIntegral, CosIntegral, SinhIntegral, CoshIntegral, Gamma, LogGamma, PolyGamma, Zeta, PolyLog, ProductLog, EllipticF, EllipticE, EllipticPi },func] HypergeometricFunctionQ[func_] := MemberQ[{Hypergeometric1F1,Hypergeometric2F1,HypergeometricPFQ},func] AppellFunctionQ[func_] := MemberQ[{AppellF1},func]
# File: GradeAntiderivative.mpl # Original version thanks to Albert Rich emailed on 03/21/2017 #Nasser 03/22/2017 Use Maple leaf count instead since buildin #Nasser 03/23/2017 missing 'ln' for ElementaryFunctionQ added #Nasser 03/24/2017 corrected the check for complex result #Nasser 10/27/2017 check for leafsize and do not call ExpnType() # if leaf size is "too large". Set at 500,000 #Nasser 12/22/2019 Added debug flag, added 'dilog' to special functions # see problem 156, file Apostol_Problems #Nasser 4/07/2022 add second output which gives reason for the grade GradeAntiderivative := proc(result,optimal) local leaf_count_result, leaf_count_optimal, ExpnType_result, ExpnType_optimal, debug:=false; leaf_count_result:=leafcount(result); #do NOT call ExpnType() if leaf size is too large. Recursion problem if leaf_count_result > 500000 then return "B","result has leaf size over 500,000. Avoiding possible recursion issues."; fi; leaf_count_optimal := leafcount(optimal); ExpnType_result := ExpnType(result); ExpnType_optimal := ExpnType(optimal); if debug then print("ExpnType_result",ExpnType_result," ExpnType_optimal=",ExpnType_optimal); fi; # If result and optimal are mathematical expressions, # GradeAntiderivative[result,optimal] returns # "F" if the result fails to integrate an expression that # is integrable # "C" if result involves higher level functions than necessary # "B" if result is more than twice the size of the optimal # antiderivative # "A" if result can be considered optimal #This check below actually is not needed, since I only #call this grading only for passed integrals. i.e. I check #for "F" before calling this. But no harm of keeping it here. #just in case. if not type(result,freeof('int')) then return "F","Result contains unresolved integral"; fi; if ExpnType_result<=ExpnType_optimal then if debug then print("ExpnType_result<=ExpnType_optimal"); fi; if is_contains_complex(result) then if is_contains_complex(optimal) then if debug then print("both result and optimal complex"); fi; if leaf_count_result<=2*leaf_count_optimal then return "A",""; else return "B",cat("Both result and optimal contain complex but leaf count of result is larger than twice the leaf count of optimal. ", convert(leaf_count_result,string)," vs. $2 (", convert(leaf_count_optimal,string)," ) = ",convert(2*leaf_count_optimal,string),"$."); end if else #result contains complex but optimal is not if debug then print("result contains complex but optimal is not"); fi; return "C","Result contains complex when optimal does not."; fi; else # result do not contain complex # this assumes optimal do not as well. No check is needed here. if debug then print("result do not contain complex, this assumes optimal do not as well"); fi; if leaf_count_result<=2*leaf_count_optimal then if debug then print("leaf_count_result<=2*leaf_count_optimal"); fi; return "A",""; else if debug then print("leaf_count_result>2*leaf_count_optimal"); fi; return "B",cat("Leaf count of result is larger than twice the leaf count of optimal. $", convert(leaf_count_result,string),"$ vs. $2(", convert(leaf_count_optimal,string),")=",convert(2*leaf_count_optimal,string),"$."); fi; fi; else #ExpnType(result) > ExpnType(optimal) if debug then print("ExpnType(result) > ExpnType(optimal)"); fi; return "C",cat("Result contains higher order function than in optimal. Order ", convert(ExpnType_result,string)," vs. order ", convert(ExpnType_optimal,string),"."); fi; end proc: # # is_contains_complex(result) # takes expressions and returns true if it contains "I" else false # #Nasser 032417 is_contains_complex:= proc(expression) return (has(expression,I)); end proc: # The following summarizes the type number assigned an expression # based on the functions it involves # 1 = rational function # 2 = algebraic function # 3 = elementary function # 4 = special function # 5 = hyperpergeometric function # 6 = appell function # 7 = rootsum function # 8 = integrate function # 9 = unknown function ExpnType := proc(expn) if type(expn,'atomic') then 1 elif type(expn,'list') then apply(max,map(ExpnType,expn)) elif type(expn,'sqrt') then if type(op(1,expn),'rational') then 1 else max(2,ExpnType(op(1,expn))) end if elif type(expn,'`^`') then if type(op(2,expn),'integer') then ExpnType(op(1,expn)) elif type(op(2,expn),'rational') then if type(op(1,expn),'rational') then 1 else max(2,ExpnType(op(1,expn))) end if else max(3,ExpnType(op(1,expn)),ExpnType(op(2,expn))) end if elif type(expn,'`+`') or type(expn,'`*`') then max(ExpnType(op(1,expn)),max(ExpnType(rest(expn)))) elif ElementaryFunctionQ(op(0,expn)) then max(3,ExpnType(op(1,expn))) elif SpecialFunctionQ(op(0,expn)) then max(4,apply(max,map(ExpnType,[op(expn)]))) elif HypergeometricFunctionQ(op(0,expn)) then max(5,apply(max,map(ExpnType,[op(expn)]))) elif AppellFunctionQ(op(0,expn)) then max(6,apply(max,map(ExpnType,[op(expn)]))) elif op(0,expn)='int' then max(8,apply(max,map(ExpnType,[op(expn)]))) else 9 end if end proc: ElementaryFunctionQ := proc(func) member(func,[ exp,log,ln, sin,cos,tan,cot,sec,csc, arcsin,arccos,arctan,arccot,arcsec,arccsc, sinh,cosh,tanh,coth,sech,csch, arcsinh,arccosh,arctanh,arccoth,arcsech,arccsch]) end proc: SpecialFunctionQ := proc(func) member(func,[ erf,erfc,erfi, FresnelS,FresnelC, Ei,Ei,Li,Si,Ci,Shi,Chi, GAMMA,lnGAMMA,Psi,Zeta,polylog,dilog,LambertW, EllipticF,EllipticE,EllipticPi]) end proc: HypergeometricFunctionQ := proc(func) member(func,[Hypergeometric1F1,hypergeom,HypergeometricPFQ]) end proc: AppellFunctionQ := proc(func) member(func,[AppellF1]) end proc: # u is a sum or product. rest(u) returns all but the # first term or factor of u. rest := proc(u) local v; if nops(u)=2 then op(2,u) else apply(op(0,u),op(2..nops(u),u)) end if end proc: #leafcount(u) returns the number of nodes in u. #Nasser 3/23/17 Replaced by build-in leafCount from package in Maple leafcount := proc(u) MmaTranslator[Mma][LeafCount](u); end proc:
#Dec 24, 2019. Nasser M. Abbasi: # Port of original Maple grading function by # Albert Rich to use with Sympy/Python #Dec 27, 2019 Nasser. Added `RootSum`. See problem 177, Timofeev file # added 'exp_polar' from sympy import * def leaf_count(expr): #sympy do not have leaf count function. This is approximation return round(1.7*count_ops(expr)) def is_sqrt(expr): if isinstance(expr,Pow): if expr.args[1] == Rational(1,2): return True else: return False else: return False def is_elementary_function(func): return func in [exp,log,ln,sin,cos,tan,cot,sec,csc, asin,acos,atan,acot,asec,acsc,sinh,cosh,tanh,coth,sech,csch, asinh,acosh,atanh,acoth,asech,acsch ] def is_special_function(func): return func in [ erf,erfc,erfi, fresnels,fresnelc,Ei,Ei,Li,Si,Ci,Shi,Chi, gamma,loggamma,digamma,zeta,polylog,LambertW, elliptic_f,elliptic_e,elliptic_pi,exp_polar ] def is_hypergeometric_function(func): return func in [hyper] def is_appell_function(func): return func in [appellf1] def is_atom(expn): try: if expn.isAtom or isinstance(expn,int) or isinstance(expn,float): return True else: return False except AttributeError as error: return False def expnType(expn): debug=False if debug: print("expn=",expn,"type(expn)=",type(expn)) if is_atom(expn): return 1 elif isinstance(expn,list): return max(map(expnType, expn)) #apply(max,map(ExpnType,expn)) elif is_sqrt(expn): if isinstance(expn.args[0],Rational): #type(op(1,expn),'rational') return 1 else: return max(2,expnType(expn.args[0])) #max(2,ExpnType(op(1,expn))) elif isinstance(expn,Pow): #type(expn,'`^`') if isinstance(expn.args[1],Integer): #type(op(2,expn),'integer') return expnType(expn.args[0]) #ExpnType(op(1,expn)) elif isinstance(expn.args[1],Rational): #type(op(2,expn),'rational') if isinstance(expn.args[0],Rational): #type(op(1,expn),'rational') return 1 else: return max(2,expnType(expn.args[0])) #max(2,ExpnType(op(1,expn))) else: return max(3,expnType(expn.args[0]),expnType(expn.args[1])) #max(3,ExpnType(op(1,expn)),ExpnType(op(2,expn))) elif isinstance(expn,Add) or isinstance(expn,Mul): #type(expn,'`+`') or type(expn,'`*`') m1 = expnType(expn.args[0]) m2 = expnType(list(expn.args[1:])) return max(m1,m2) #max(ExpnType(op(1,expn)),max(ExpnType(rest(expn)))) elif is_elementary_function(expn.func): #ElementaryFunctionQ(op(0,expn)) return max(3,expnType(expn.args[0])) #max(3,ExpnType(op(1,expn))) elif is_special_function(expn.func): #SpecialFunctionQ(op(0,expn)) m1 = max(map(expnType, list(expn.args))) return max(4,m1) #max(4,apply(max,map(ExpnType,[op(expn)]))) elif is_hypergeometric_function(expn.func): #HypergeometricFunctionQ(op(0,expn)) m1 = max(map(expnType, list(expn.args))) return max(5,m1) #max(5,apply(max,map(ExpnType,[op(expn)]))) elif is_appell_function(expn.func): m1 = max(map(expnType, list(expn.args))) return max(6,m1) #max(5,apply(max,map(ExpnType,[op(expn)]))) elif isinstance(expn,RootSum): m1 = max(map(expnType, list(expn.args))) #Apply[Max,Append[Map[ExpnType,Apply[List,expn]],7]], return max(7,m1) elif str(expn).find("Integral") != -1: m1 = max(map(expnType, list(expn.args))) return max(8,m1) #max(5,apply(max,map(ExpnType,[op(expn)]))) else: return 9 #main function def grade_antiderivative(result,optimal): #print ("Enter grade_antiderivative for sagemath") #print("Enter grade_antiderivative, result=",result," optimal=",optimal) leaf_count_result = leaf_count(result) leaf_count_optimal = leaf_count(optimal) #print("leaf_count_result=",leaf_count_result) #print("leaf_count_optimal=",leaf_count_optimal) expnType_result = expnType(result) expnType_optimal = expnType(optimal) if str(result).find("Integral") != -1: grade = "F" grade_annotation ="" else: if expnType_result <= expnType_optimal: if result.has(I): if optimal.has(I): #both result and optimal complex if leaf_count_result <= 2*leaf_count_optimal: grade = "A" grade_annotation ="" else: grade = "B" grade_annotation ="Both result and optimal contain complex but leaf count of result is larger than twice the leaf count of optimal. "+str(leaf_count_result)+" vs. $2 ("+str(leaf_count_optimal)+") = "+ str(2*leaf_count_optimal)+"$." else: #result contains complex but optimal is not grade = "C" grade_annotation ="Result contains complex when optimal does not." else: # result do not contain complex, this assumes optimal do not as well if leaf_count_result <= 2*leaf_count_optimal: grade = "A" grade_annotation ="" else: grade = "B" grade_annotation ="Leaf count of result is larger than twice the leaf count of optimal. "+str(leaf_count_result)+" vs. $2 ("+str(leaf_count_optimal)+") = "+ str(2*leaf_count_optimal)+"$." else: grade = "C" grade_annotation ="Result contains higher order function than in optimal. Order "+str(ExpnType_result)+" vs. order "+str(ExpnType_optimal)+"." #print("Before returning. grade=",grade, " grade_annotation=",grade_annotation) return grade, grade_annotation
#Dec 24, 2019. Nasser: Ported original Maple grading function by # Albert Rich to use with Sagemath. This is used to # grade Fricas, Giac and Maxima results. #Dec 24, 2019. Nasser: Added 'exp_integral_e' and 'sng', 'sin_integral' # 'arctan2','floor','abs','log_integral' #June 4, 2022 Made default grade_annotation "none" instead of "" due # issue later when reading the file. #July 14, 2022. Added ellipticF. This is until they fix sagemath, then remove it. from sage.all import * from sage.symbolic.operators import add_vararg, mul_vararg debug=False; def tree_size(expr): r""" Return the tree size of this expression. """ #print("Enter tree_size, expr is ",expr) if expr not in SR: # deal with lists, tuples, vectors return 1 + sum(tree_size(a) for a in expr) expr = SR(expr) x, aa = expr.operator(), expr.operands() if x is None: return 1 else: return 1 + sum(tree_size(a) for a in aa) def is_sqrt(expr): if expr.operator() == operator.pow: #isinstance(expr,Pow): if expr.operands()[1]==1/2: #expr.args[1] == Rational(1,2): if debug: print ("expr is sqrt") return True else: return False else: return False def is_elementary_function(func): #debug=False m = func.name() in ['exp','log','ln', 'sin','cos','tan','cot','sec','csc', 'arcsin','arccos','arctan','arccot','arcsec','arccsc', 'sinh','cosh','tanh','coth','sech','csch', 'arcsinh','arccosh','arctanh','arccoth','arcsech','arccsch','sgn', 'arctan2','floor','abs' ] if debug: if m: print ("func ", func , " is elementary_function") else: print ("func ", func , " is NOT elementary_function") return m def is_special_function(func): #debug=False if debug: print ("type(func)=", type(func)) m= func.name() in ['erf','erfc','erfi','fresnel_sin','fresnel_cos','Ei', 'Ei','Li','Si','sin_integral','Ci','cos_integral','Shi','sinh_integral' 'Chi','cosh_integral','gamma','log_gamma','psi,zeta', 'polylog','lambert_w','elliptic_f','elliptic_e','ellipticF', 'elliptic_pi','exp_integral_e','log_integral'] if debug: print ("m=",m) if m: print ("func ", func ," is special_function") else: print ("func ", func ," is NOT special_function") return m def is_hypergeometric_function(func): return func.name() in ['hypergeometric','hypergeometric_M','hypergeometric_U'] def is_appell_function(func): return func.name() in ['hypergeometric'] #[appellf1] can't find this in sagemath def is_atom(expn): #debug=False if debug: print ("Enter is_atom, expn=",expn) if not hasattr(expn, 'parent'): return False #thanks to answer at https://ask.sagemath.org/question/49179/what-is-sagemath-equivalent-to-atomic-type-in-maple/ try: if expn.parent() is SR: return expn.operator() is None if expn.parent() in (ZZ, QQ, AA, QQbar): return expn in expn.parent() # Should always return True if hasattr(expn.parent(),"base_ring") and hasattr(expn.parent(),"gens"): return expn in expn.parent().base_ring() or expn in expn.parent().gens() return False except AttributeError as error: print("Exception,AttributeError in is_atom") print ("cought exception" , type(error).__name__ ) return False def expnType(expn): if debug: print (">>>>>Enter expnType, expn=", expn) print (">>>>>is_atom(expn)=", is_atom(expn)) if is_atom(expn): return 1 elif type(expn)==list: #isinstance(expn,list): return max(map(expnType, expn)) #apply(max,map(ExpnType,expn)) elif is_sqrt(expn): if type(expn.operands()[0])==Rational: #type(isinstance(expn.args[0],Rational): return 1 else: return max(2,expnType(expn.operands()[0])) #max(2,expnType(expn.args[0])) elif expn.operator() == operator.pow: #isinstance(expn,Pow) if type(expn.operands()[1])==Integer: #isinstance(expn.args[1],Integer) return expnType(expn.operands()[0]) #expnType(expn.args[0]) elif type(expn.operands()[1])==Rational: #isinstance(expn.args[1],Rational) if type(expn.operands()[0])==Rational: #isinstance(expn.args[0],Rational) return 1 else: return max(2,expnType(expn.operands()[0])) #max(2,expnType(expn.args[0])) else: return max(3,expnType(expn.operands()[0]),expnType(expn.operands()[1])) #max(3,expnType(expn.operands()[0]),expnType(expn.operands()[1])) elif expn.operator() == add_vararg or expn.operator() == mul_vararg: #isinstance(expn,Add) or isinstance(expn,Mul) m1 = expnType(expn.operands()[0]) #expnType(expn.args[0]) m2 = expnType(expn.operands()[1:]) #expnType(list(expn.args[1:])) return max(m1,m2) #max(ExpnType(op(1,expn)),max(ExpnType(rest(expn)))) elif is_elementary_function(expn.operator()): #is_elementary_function(expn.func) return max(3,expnType(expn.operands()[0])) elif is_special_function(expn.operator()): #is_special_function(expn.func) m1 = max(map(expnType, expn.operands())) #max(map(expnType, list(expn.args))) return max(4,m1) #max(4,m1) elif is_hypergeometric_function(expn.operator()): #is_hypergeometric_function(expn.func) m1 = max(map(expnType, expn.operands())) #max(map(expnType, list(expn.args))) return max(5,m1) #max(5,m1) elif is_appell_function(expn.operator()): m1 = max(map(expnType, expn.operands())) #max(map(expnType, list(expn.args))) return max(6,m1) #max(6,m1) elif str(expn).find("Integral") != -1: #this will never happen, since it #is checked before calling the grading function that is passed. #but kept it here. m1 = max(map(expnType, expn.operands())) #max(map(expnType, list(expn.args))) return max(8,m1) #max(5,apply(max,map(ExpnType,[op(expn)]))) else: return 9 #main function def grade_antiderivative(result,optimal): if debug: print ("Enter grade_antiderivative for sagemath") print("Enter grade_antiderivative, result=",result) print("Enter grade_antiderivative, optimal=",optimal) print("type(anti)=",type(result)) print("type(optimal)=",type(optimal)) leaf_count_result = tree_size(result) #leaf_count(result) leaf_count_optimal = tree_size(optimal) #leaf_count(optimal) #if debug: print ("leaf_count_result=", leaf_count_result, "leaf_count_optimal=",leaf_count_optimal) expnType_result = expnType(result) expnType_optimal = expnType(optimal) if debug: print ("expnType_result=", expnType_result, "expnType_optimal=",expnType_optimal) if expnType_result <= expnType_optimal: if result.has(I): if optimal.has(I): #both result and optimal complex if leaf_count_result <= 2*leaf_count_optimal: grade = "A" grade_annotation ="none" else: grade = "B" grade_annotation ="Both result and optimal contain complex but leaf count of result is larger than twice the leaf count of optimal. "+str(leaf_count_result)+" vs. $2 ("+str(leaf_count_optimal)+") = "+ str(2*leaf_count_optimal)+"$." else: #result contains complex but optimal is not grade = "C" grade_annotation ="Result contains complex when optimal does not." else: # result do not contain complex, this assumes optimal do not as well if leaf_count_result <= 2*leaf_count_optimal: grade = "A" grade_annotation ="none" else: grade = "B" grade_annotation ="Leaf count of result is larger than twice the leaf count of optimal. "+str(leaf_count_result)+" vs. $2 ("+str(leaf_count_optimal)+") = "+ str(2*leaf_count_optimal)+"$." else: grade = "C" grade_annotation ="Result contains higher order function than in optimal. Order "+str(expnType_result)+" vs. order "+str(expnType_optimal)+"." print("Before returning. grade=",grade, " grade_annotation=",grade_annotation) return grade, grade_annotation