Navigation

Source code for sympy.printing.octave

"""Octave (and Matlab) code printerThe `OctaveCodePrinter` converts SymPy expressions into Octave expressions.It uses a subset of the Octave language for Matlab compatibility.A complete code generator, which uses `octave_code` extensively, can be foundin `sympy.utilities.codegen`. The `codegen` module can be used to generatecomplete source code files."""from__future__importprint_function,divisionfromsympy.coreimportMul,Pow,S,Rationalfromsympy.core.compatibilityimportstring_types,rangefromsympy.core.mulimport_keep_coefffromsympy.printing.codeprinterimportCodePrinter,Assignmentfromsympy.printing.precedenceimportprecedencefromreimportsearch# List of known functions. First, those that have the same name in# SymPy and Octave. This is almost certainly incomplete!known_fcns_src1=["sin","cos","tan","asin","acos","atan","atan2","sinh","cosh","tanh","asinh","acosh","atanh","log","exp","erf","gamma","sign","floor","csc","sec","cot","coth","acot","acoth","erfc","besselj","bessely","besseli","besselk","erfinv","erfcinv","factorial"]# These functions have different names ("Sympy": "Octave"), more# generally a mapping to (argument_conditions, octave_function).known_fcns_src2={"Abs":"abs","ceiling":"ceil","conjugate":"conj","DiracDelta":"dirac","Heaviside":"heaviside",}

[docs]classOctaveCodePrinter(CodePrinter):""" A printer to convert expressions to strings of Octave/Matlab code. """printmethod="_octave"language="Octave"_operators={'and':'&','or':'|','not':'~',}_default_settings={'order':None,'full_prec':'auto','precision':16,'user_functions':{},'human':True,'contract':True,'inline':True,}# Note: contract is for expressing tensors as loops (if True), or just# assignment (if False). FIXME: this should be looked a more carefully# for Octave.def__init__(self,settings={}):super(OctaveCodePrinter,self).__init__(settings)self.known_functions=dict(zip(known_fcns_src1,known_fcns_src1))self.known_functions.update(dict(known_fcns_src2))userfuncs=settings.get('user_functions',{})self.known_functions.update(userfuncs)def_rate_index_position(self,p):returnp*5def_get_statement(self,codestring):return"%s;"%codestringdef_get_comment(self,text):return"% {0}".format(text)def_declare_number_const(self,name,value):return"{0} = {1};".format(name,value)def_format_code(self,lines):returnself.indent_code(lines)def_traverse_matrix_indices(self,mat):# Octave uses Fortran order (column-major)rows,cols=mat.shapereturn((i,j)forjinrange(cols)foriinrange(rows))def_get_loop_opening_ending(self,indices):open_lines=[]close_lines=[]foriinindices:# Octave arrays start at 1 and end at dimensionvar,start,stop=map(self._print,[i.label,i.lower+1,i.upper+1])open_lines.append("for %s = %s:%s"%(var,start,stop))close_lines.append("end")returnopen_lines,close_linesdef_print_Mul(self,expr):# print complex numbers nicely in Octaveif(expr.is_numberandexpr.is_imaginaryandexpr.as_coeff_Mul()[0].is_integer):return"%si"%self._print(-S.ImaginaryUnit*expr)# cribbed from str.pyprec=precedence(expr)c,e=expr.as_coeff_Mul()ifc<0:expr=_keep_coeff(-c,e)sign="-"else:sign=""a=[]# items in the numeratorb=[]# items that are in the denominator (if any)ifself.ordernotin('old','none'):args=expr.as_ordered_factors()else:# use make_args in case expr was something like -x -> xargs=Mul.make_args(expr)# Gather args for numerator/denominatorforiteminargs:if(item.is_commutativeanditem.is_Powanditem.exp.is_Rationalanditem.exp.is_negative):ifitem.exp!=-1:b.append(Pow(item.base,-item.exp,evaluate=False))else:b.append(Pow(item.base,-item.exp))elifitem.is_RationalanditemisnotS.Infinity:ifitem.p!=1:a.append(Rational(item.p))ifitem.q!=1:b.append(Rational(item.q))else:a.append(item)a=aor[S.One]a_str=[self.parenthesize(x,prec)forxina]b_str=[self.parenthesize(x,prec)forxinb]# from here it differs from str.py to deal with "*" and ".*"defmultjoin(a,a_str):# here we probably are assuming the constants will come firstr=a_str[0]foriinrange(1,len(a)):mulsym='*'ifa[i-1].is_numberelse'.*'r=r+mulsym+a_str[i]returnriflen(b)==0:returnsign+multjoin(a,a_str)eliflen(b)==1:divsym='/'ifb[0].is_numberelse'./'returnsign+multjoin(a,a_str)+divsym+b_str[0]else:divsym='/'ifall([bi.is_numberforbiinb])else'./'return(sign+multjoin(a,a_str)+divsym+"(%s)"%multjoin(b,b_str))def_print_Pow(self,expr):powsymbol='^'ifall([x.is_numberforxinexpr.args])else'.^'PREC=precedence(expr)ifexpr.exp==S.Half:return"sqrt(%s)"%self._print(expr.base)ifexpr.is_commutative:ifexpr.exp==-S.Half:sym='/'ifexpr.base.is_numberelse'./'return"1"+sym+"sqrt(%s)"%self._print(expr.base)ifexpr.exp==-S.One:sym='/'ifexpr.base.is_numberelse'./'return"1"+sym+"%s"%self.parenthesize(expr.base,PREC)return'%s%s%s'%(self.parenthesize(expr.base,PREC),powsymbol,self.parenthesize(expr.exp,PREC))def_print_MatPow(self,expr):PREC=precedence(expr)return'%s^%s'%(self.parenthesize(expr.base,PREC),self.parenthesize(expr.exp,PREC))def_print_Pi(self,expr):return'pi'def_print_ImaginaryUnit(self,expr):return"1i"def_print_Exp1(self,expr):return"exp(1)"def_print_GoldenRatio(self,expr):# FIXME: how to do better, e.g., for octave_code(2*GoldenRatio)?#return self._print((1+sqrt(S(5)))/2)return"(1+sqrt(5))/2"def_print_NumberSymbol(self,expr):ifself._settings["inline"]:returnself._print(expr.evalf(self._settings["precision"]))else:# assign to a variable, perhaps more readable for longer programreturnsuper(OctaveCodePrinter,self)._print_NumberSymbol(expr)def_print_Assignment(self,expr):fromsympy.functions.elementary.piecewiseimportPiecewisefromsympy.tensor.indexedimportIndexedBase# Copied from codeprinter, but remove special MatrixSymbol treatmentlhs=expr.lhsrhs=expr.rhs# We special case assignments that take multiple linesifnotself._settings["inline"]andisinstance(expr.rhs,Piecewise):# Here we modify Piecewise so each expression is now# an Assignment, and then continue on the print.expressions=[]conditions=[]for(e,c)inrhs.args:expressions.append(Assignment(lhs,e))conditions.append(c)temp=Piecewise(*zip(expressions,conditions))returnself._print(temp)ifself._settings["contract"]and(lhs.has(IndexedBase)orrhs.has(IndexedBase)):# Here we check if there is looping to be done, and if so# print the required loops.returnself._doprint_loops(rhs,lhs)else:lhs_code=self._print(lhs)rhs_code=self._print(rhs)returnself._get_statement("%s = %s"%(lhs_code,rhs_code))def_print_Infinity(self,expr):return'inf'def_print_NegativeInfinity(self,expr):return'-inf'def_print_NaN(self,expr):return'NaN'def_print_list(self,expr):return'{'+', '.join(self._print(a)forainexpr)+'}'_print_tuple=_print_list_print_Tuple=_print_listdef_print_BooleanTrue(self,expr):return"true"def_print_BooleanFalse(self,expr):return"false"def_print_bool(self,expr):returnstr(expr).lower()# Could generate quadrature code for definite Integrals?#_print_Integral = _print_not_supporteddef_print_MatrixBase(self,A):# Handle zero dimensions:if(A.rows,A.cols)==(0,0):return'[]'elifA.rows==0orA.cols==0:return'zeros(%s, %s)'%(A.rows,A.cols)elif(A.rows,A.cols)==(1,1):# Octave does not distinguish between scalars and 1x1 matricesreturnself._print(A[0,0])elifA.rows==1:return"[%s]"%A.table(self,rowstart='',rowend='',colsep=' ')elifA.cols==1:# note .table would unnecessarily equispace the rowsreturn"[%s]"%"; ".join([self._print(a)forainA])return"[%s]"%A.table(self,rowstart='',rowend='',rowsep=';\n',colsep=' ')def_print_SparseMatrix(self,A):fromsympy.matricesimportMatrixL=A.col_list();# make row vectors of the indices and entriesI=Matrix([[k[0]+1forkinL]])J=Matrix([[k[1]+1forkinL]])AIJ=Matrix([[k[2]forkinL]])return"sparse(%s, %s, %s, %s, %s)"%(self._print(I),self._print(J),self._print(AIJ),A.rows,A.cols)# FIXME: Str/CodePrinter could define each of these to call the _print# method from higher up the class hierarchy (see _print_NumberSymbol).# Then subclasses like us would not need to repeat all this._print_Matrix= \
_print_DenseMatrix= \
_print_MutableDenseMatrix= \
_print_ImmutableMatrix= \
_print_ImmutableDenseMatrix= \
_print_MatrixBase_print_MutableSparseMatrix= \
_print_ImmutableSparseMatrix= \
_print_SparseMatrixdef_print_MatrixElement(self,expr):returnself._print(expr.parent)+'(%s, %s)'%(expr.i+1,expr.j+1)def_print_MatrixSlice(self,expr):defstrslice(x,lim):l=x[0]+1h=x[1]step=x[2]lstr=self._print(l)hstr='end'ifh==limelseself._print(h)ifstep==1:ifl==1andh==lim:return':'ifl==h:returnlstrelse:returnlstr+':'+hstrelse:return':'.join((lstr,self._print(step),hstr))return(self._print(expr.parent)+'('+strslice(expr.rowslice,expr.parent.shape[0])+', '+strslice(expr.colslice,expr.parent.shape[1])+')')def_print_Indexed(self,expr):inds=[self._print(i)foriinexpr.indices]return"%s(%s)"%(self._print(expr.base.label),", ".join(inds))def_print_Idx(self,expr):returnself._print(expr.label)def_print_Identity(self,expr):return"eye(%s)"%self._print(expr.shape[0])def_print_hankel1(self,expr):return"besselh(%s, 1, %s)"%(self._print(expr.order),self._print(expr.argument))def_print_hankel2(self,expr):return"besselh(%s, 2, %s)"%(self._print(expr.order),self._print(expr.argument))# Note: as of 2015, Octave doesn't have spherical Bessel functionsdef_print_jn(self,expr):fromsympy.functionsimportsqrt,besseljx=expr.argumentexpr2=sqrt(S.Pi/(2*x))*besselj(expr.order+S.Half,x)returnself._print(expr2)def_print_yn(self,expr):fromsympy.functionsimportsqrt,besselyx=expr.argumentexpr2=sqrt(S.Pi/(2*x))*bessely(expr.order+S.Half,x)returnself._print(expr2)def_print_airyai(self,expr):return"airy(0, %s)"%self._print(expr.args[0])def_print_airyaiprime(self,expr):return"airy(1, %s)"%self._print(expr.args[0])def_print_airybi(self,expr):return"airy(2, %s)"%self._print(expr.args[0])def_print_airybiprime(self,expr):return"airy(3, %s)"%self._print(expr.args[0])def_print_Piecewise(self,expr):ifexpr.args[-1].cond!=True:# We need the last conditional to be a True, otherwise the resulting# function may not return a result.raiseValueError("All Piecewise expressions must contain an ""(expr, True) statement to be used as a default ""condition. Without one, the generated ""expression may not evaluate to anything under ""some condition.")lines=[]ifself._settings["inline"]:# Express each (cond, expr) pair in a nested Horner form:# (condition) .* (expr) + (not cond) .* (<others>)# Expressions that result in multiple statements won't work here.ecpairs=["({0}).*({1}) + (~({0})).*(".format(self._print(c),self._print(e))fore,cinexpr.args[:-1]]elast="%s"%self._print(expr.args[-1].expr)pw=" ...\n".join(ecpairs)+elast+")"*len(ecpairs)# Note: current need these outer brackets for 2*pw. Would be# nicer to teach parenthesize() to do this for us when needed!return"("+pw+")"else:fori,(e,c)inenumerate(expr.args):ifi==0:lines.append("if (%s)"%self._print(c))elifi==len(expr.args)-1andc==True:lines.append("else")else:lines.append("elseif (%s)"%self._print(c))code0=self._print(e)lines.append(code0)ifi==len(expr.args)-1:lines.append("end")return"\n".join(lines)

[docs]defoctave_code(expr,assign_to=None,**settings):r"""Converts `expr` to a string of Octave (or Matlab) code. The string uses a subset of the Octave language for Matlab compatibility. Parameters ========== expr : Expr A sympy expression to be converted. assign_to : optional When given, the argument is used as the name of the variable to which the expression is assigned. Can be a string, ``Symbol``, ``MatrixSymbol``, or ``Indexed`` type. This can be helpful for expressions that generate multi-line statements. precision : integer, optional The precision for numbers such as pi [default=16]. user_functions : dict, optional A dictionary where keys are ``FunctionClass`` instances and values are their string representations. Alternatively, the dictionary value can be a list of tuples i.e. [(argument_test, cfunction_string)]. See below for examples. human : bool, optional If True, the result is a single string that may contain some constant declarations for the number symbols. If False, the same information is returned in a tuple of (symbols_to_declare, not_supported_functions, code_text). [default=True]. contract: bool, optional If True, ``Indexed`` instances are assumed to obey tensor contraction rules and the corresponding nested loops over indices are generated. Setting contract=False will not generate loops, instead the user is responsible to provide values for the indices in the code. [default=True]. inline: bool, optional If True, we try to create single-statement code instead of multiple statements. [default=True]. Examples ======== >>> from sympy import octave_code, symbols, sin, pi >>> x = symbols('x') >>> octave_code(sin(x).series(x).removeO()) 'x.^5/120 - x.^3/6 + x' >>> from sympy import Rational, ceiling, Abs >>> x, y, tau = symbols("x, y, tau") >>> octave_code((2*tau)**Rational(7, 2)) '8*sqrt(2)*tau.^(7/2)' Note that element-wise (Hadamard) operations are used by default between symbols. This is because its very common in Octave to write "vectorized" code. It is harmless if the values are scalars. >>> octave_code(sin(pi*x*y), assign_to="s") 's = sin(pi*x.*y);' If you need a matrix product "*" or matrix power "^", you can specify the symbol as a ``MatrixSymbol``. >>> from sympy import Symbol, MatrixSymbol >>> n = Symbol('n', integer=True, positive=True) >>> A = MatrixSymbol('A', n, n) >>> octave_code(3*pi*A**3) '(3*pi)*A^3' This class uses several rules to decide which symbol to use a product. Pure numbers use "*", Symbols use ".*" and MatrixSymbols use "*". A HadamardProduct can be used to specify componentwise multiplication ".*" of two MatrixSymbols. There is currently there is no easy way to specify scalar symbols, so sometimes the code might have some minor cosmetic issues. For example, suppose x and y are scalars and A is a Matrix, then while a human programmer might write "(x^2*y)*A^3", we generate: >>> octave_code(x**2*y*A**3) '(x.^2.*y)*A^3' Matrices are supported using Octave inline notation. When using ``assign_to`` with matrices, the name can be specified either as a string or as a ``MatrixSymbol``. The dimenions must align in the latter case. >>> from sympy import Matrix, MatrixSymbol >>> mat = Matrix([[x**2, sin(x), ceiling(x)]]) >>> octave_code(mat, assign_to='A') 'A = [x.^2 sin(x) ceil(x)];' ``Piecewise`` expressions are implemented with logical masking by default. Alternatively, you can pass "inline=False" to use if-else conditionals. Note that if the ``Piecewise`` lacks a default term, represented by ``(expr, True)`` then an error will be thrown. This is to prevent generating an expression that may not evaluate to anything. >>> from sympy import Piecewise >>> pw = Piecewise((x + 1, x > 0), (x, True)) >>> octave_code(pw, assign_to=tau) 'tau = ((x > 0).*(x + 1) + (~(x > 0)).*(x));' Note that any expression that can be generated normally can also exist inside a Matrix: >>> mat = Matrix([[x**2, pw, sin(x)]]) >>> octave_code(mat, assign_to='A') 'A = [x.^2 ((x > 0).*(x + 1) + (~(x > 0)).*(x)) sin(x)];' Custom printing can be defined for certain types by passing a dictionary of "type" : "function" to the ``user_functions`` kwarg. Alternatively, the dictionary value can be a list of tuples i.e., [(argument_test, cfunction_string)]. This can be used to call a custom Octave function. >>> from sympy import Function >>> f = Function('f') >>> g = Function('g') >>> custom_functions = { ... "f": "existing_octave_fcn", ... "g": [(lambda x: x.is_Matrix, "my_mat_fcn"), ... (lambda x: not x.is_Matrix, "my_fcn")] ... } >>> mat = Matrix([[1, x]]) >>> octave_code(f(x) + g(x) + g(mat), user_functions=custom_functions) 'existing_octave_fcn(x) + my_fcn(x) + my_mat_fcn([1 x])' Support for loops is provided through ``Indexed`` types. With ``contract=True`` these expressions will be turned into loops, whereas ``contract=False`` will just print the assignment expression that should be looped over: >>> from sympy import Eq, IndexedBase, Idx, ccode >>> len_y = 5 >>> y = IndexedBase('y', shape=(len_y,)) >>> t = IndexedBase('t', shape=(len_y,)) >>> Dy = IndexedBase('Dy', shape=(len_y-1,)) >>> i = Idx('i', len_y-1) >>> e = Eq(Dy[i], (y[i+1]-y[i])/(t[i+1]-t[i])) >>> octave_code(e.rhs, assign_to=e.lhs, contract=False) 'Dy(i) = (y(i + 1) - y(i))./(t(i + 1) - t(i));' """returnOctaveCodePrinter(settings).doprint(expr,assign_to)

defprint_octave_code(expr,**settings):"""Prints the Octave (or Matlab) representation of the given expression. See `octave_code` for the meaning of the optional arguments. """print(octave_code(expr,**settings))