octave.py 25 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711
  1. """
  2. Octave (and Matlab) code printer
  3. The `OctaveCodePrinter` converts SymPy expressions into Octave expressions.
  4. It uses a subset of the Octave language for Matlab compatibility.
  5. A complete code generator, which uses `octave_code` extensively, can be found
  6. in `sympy.utilities.codegen`. The `codegen` module can be used to generate
  7. complete source code files.
  8. """
  9. from __future__ import annotations
  10. from typing import Any
  11. from sympy.core import Mul, Pow, S, Rational
  12. from sympy.core.mul import _keep_coeff
  13. from sympy.core.numbers import equal_valued
  14. from sympy.printing.codeprinter import CodePrinter
  15. from sympy.printing.precedence import precedence, PRECEDENCE
  16. from re import search
  17. # List of known functions. First, those that have the same name in
  18. # SymPy and Octave. This is almost certainly incomplete!
  19. known_fcns_src1 = ["sin", "cos", "tan", "cot", "sec", "csc",
  20. "asin", "acos", "acot", "atan", "atan2", "asec", "acsc",
  21. "sinh", "cosh", "tanh", "coth", "csch", "sech",
  22. "asinh", "acosh", "atanh", "acoth", "asech", "acsch",
  23. "erfc", "erfi", "erf", "erfinv", "erfcinv",
  24. "besseli", "besselj", "besselk", "bessely",
  25. "bernoulli", "beta", "euler", "exp", "factorial", "floor",
  26. "fresnelc", "fresnels", "gamma", "harmonic", "log",
  27. "polylog", "sign", "zeta", "legendre"]
  28. # These functions have different names ("SymPy": "Octave"), more
  29. # generally a mapping to (argument_conditions, octave_function).
  30. known_fcns_src2 = {
  31. "Abs": "abs",
  32. "arg": "angle", # arg/angle ok in Octave but only angle in Matlab
  33. "binomial": "bincoeff",
  34. "ceiling": "ceil",
  35. "chebyshevu": "chebyshevU",
  36. "chebyshevt": "chebyshevT",
  37. "Chi": "coshint",
  38. "Ci": "cosint",
  39. "conjugate": "conj",
  40. "DiracDelta": "dirac",
  41. "Heaviside": "heaviside",
  42. "im": "imag",
  43. "laguerre": "laguerreL",
  44. "LambertW": "lambertw",
  45. "li": "logint",
  46. "loggamma": "gammaln",
  47. "Max": "max",
  48. "Min": "min",
  49. "Mod": "mod",
  50. "polygamma": "psi",
  51. "re": "real",
  52. "RisingFactorial": "pochhammer",
  53. "Shi": "sinhint",
  54. "Si": "sinint",
  55. }
  56. class OctaveCodePrinter(CodePrinter):
  57. """
  58. A printer to convert expressions to strings of Octave/Matlab code.
  59. """
  60. printmethod = "_octave"
  61. language = "Octave"
  62. _operators = {
  63. 'and': '&',
  64. 'or': '|',
  65. 'not': '~',
  66. }
  67. _default_settings: dict[str, Any] = dict(CodePrinter._default_settings, **{
  68. 'precision': 17,
  69. 'user_functions': {},
  70. 'contract': True,
  71. 'inline': True,
  72. })
  73. # Note: contract is for expressing tensors as loops (if True), or just
  74. # assignment (if False). FIXME: this should be looked a more carefully
  75. # for Octave.
  76. def __init__(self, settings={}):
  77. super().__init__(settings)
  78. self.known_functions = dict(zip(known_fcns_src1, known_fcns_src1))
  79. self.known_functions.update(dict(known_fcns_src2))
  80. userfuncs = settings.get('user_functions', {})
  81. self.known_functions.update(userfuncs)
  82. def _rate_index_position(self, p):
  83. return p*5
  84. def _get_statement(self, codestring):
  85. return "%s;" % codestring
  86. def _get_comment(self, text):
  87. return "% {}".format(text)
  88. def _declare_number_const(self, name, value):
  89. return "{} = {};".format(name, value)
  90. def _format_code(self, lines):
  91. return self.indent_code(lines)
  92. def _traverse_matrix_indices(self, mat):
  93. # Octave uses Fortran order (column-major)
  94. rows, cols = mat.shape
  95. return ((i, j) for j in range(cols) for i in range(rows))
  96. def _get_loop_opening_ending(self, indices):
  97. open_lines = []
  98. close_lines = []
  99. for i in indices:
  100. # Octave arrays start at 1 and end at dimension
  101. var, start, stop = map(self._print,
  102. [i.label, i.lower + 1, i.upper + 1])
  103. open_lines.append("for %s = %s:%s" % (var, start, stop))
  104. close_lines.append("end")
  105. return open_lines, close_lines
  106. def _print_Mul(self, expr):
  107. # print complex numbers nicely in Octave
  108. if (expr.is_number and expr.is_imaginary and
  109. (S.ImaginaryUnit*expr).is_Integer):
  110. return "%si" % self._print(-S.ImaginaryUnit*expr)
  111. # cribbed from str.py
  112. prec = precedence(expr)
  113. c, e = expr.as_coeff_Mul()
  114. if c < 0:
  115. expr = _keep_coeff(-c, e)
  116. sign = "-"
  117. else:
  118. sign = ""
  119. a = [] # items in the numerator
  120. b = [] # items that are in the denominator (if any)
  121. pow_paren = [] # Will collect all pow with more than one base element and exp = -1
  122. if self.order not in ('old', 'none'):
  123. args = expr.as_ordered_factors()
  124. else:
  125. # use make_args in case expr was something like -x -> x
  126. args = Mul.make_args(expr)
  127. # Gather args for numerator/denominator
  128. for item in args:
  129. if (item.is_commutative and item.is_Pow and item.exp.is_Rational
  130. and item.exp.is_negative):
  131. if item.exp != -1:
  132. b.append(Pow(item.base, -item.exp, evaluate=False))
  133. else:
  134. if len(item.args[0].args) != 1 and isinstance(item.base, Mul): # To avoid situations like #14160
  135. pow_paren.append(item)
  136. b.append(Pow(item.base, -item.exp))
  137. elif item.is_Rational and item is not S.Infinity:
  138. if item.p != 1:
  139. a.append(Rational(item.p))
  140. if item.q != 1:
  141. b.append(Rational(item.q))
  142. else:
  143. a.append(item)
  144. a = a or [S.One]
  145. a_str = [self.parenthesize(x, prec) for x in a]
  146. b_str = [self.parenthesize(x, prec) for x in b]
  147. # To parenthesize Pow with exp = -1 and having more than one Symbol
  148. for item in pow_paren:
  149. if item.base in b:
  150. b_str[b.index(item.base)] = "(%s)" % b_str[b.index(item.base)]
  151. # from here it differs from str.py to deal with "*" and ".*"
  152. def multjoin(a, a_str):
  153. # here we probably are assuming the constants will come first
  154. r = a_str[0]
  155. for i in range(1, len(a)):
  156. mulsym = '*' if a[i-1].is_number else '.*'
  157. r = r + mulsym + a_str[i]
  158. return r
  159. if not b:
  160. return sign + multjoin(a, a_str)
  161. elif len(b) == 1:
  162. divsym = '/' if b[0].is_number else './'
  163. return sign + multjoin(a, a_str) + divsym + b_str[0]
  164. else:
  165. divsym = '/' if all(bi.is_number for bi in b) else './'
  166. return (sign + multjoin(a, a_str) +
  167. divsym + "(%s)" % multjoin(b, b_str))
  168. def _print_Relational(self, expr):
  169. lhs_code = self._print(expr.lhs)
  170. rhs_code = self._print(expr.rhs)
  171. op = expr.rel_op
  172. return "{} {} {}".format(lhs_code, op, rhs_code)
  173. def _print_Pow(self, expr):
  174. powsymbol = '^' if all(x.is_number for x in expr.args) else '.^'
  175. PREC = precedence(expr)
  176. if equal_valued(expr.exp, 0.5):
  177. return "sqrt(%s)" % self._print(expr.base)
  178. if expr.is_commutative:
  179. if equal_valued(expr.exp, -0.5):
  180. sym = '/' if expr.base.is_number else './'
  181. return "1" + sym + "sqrt(%s)" % self._print(expr.base)
  182. if equal_valued(expr.exp, -1):
  183. sym = '/' if expr.base.is_number else './'
  184. return "1" + sym + "%s" % self.parenthesize(expr.base, PREC)
  185. return '%s%s%s' % (self.parenthesize(expr.base, PREC), powsymbol,
  186. self.parenthesize(expr.exp, PREC))
  187. def _print_MatPow(self, expr):
  188. PREC = precedence(expr)
  189. return '%s^%s' % (self.parenthesize(expr.base, PREC),
  190. self.parenthesize(expr.exp, PREC))
  191. def _print_MatrixSolve(self, expr):
  192. PREC = precedence(expr)
  193. return "%s \\ %s" % (self.parenthesize(expr.matrix, PREC),
  194. self.parenthesize(expr.vector, PREC))
  195. def _print_Pi(self, expr):
  196. return 'pi'
  197. def _print_ImaginaryUnit(self, expr):
  198. return "1i"
  199. def _print_Exp1(self, expr):
  200. return "exp(1)"
  201. def _print_GoldenRatio(self, expr):
  202. # FIXME: how to do better, e.g., for octave_code(2*GoldenRatio)?
  203. #return self._print((1+sqrt(S(5)))/2)
  204. return "(1+sqrt(5))/2"
  205. def _print_Assignment(self, expr):
  206. from sympy.codegen.ast import Assignment
  207. from sympy.functions.elementary.piecewise import Piecewise
  208. from sympy.tensor.indexed import IndexedBase
  209. # Copied from codeprinter, but remove special MatrixSymbol treatment
  210. lhs = expr.lhs
  211. rhs = expr.rhs
  212. # We special case assignments that take multiple lines
  213. if not self._settings["inline"] and isinstance(expr.rhs, Piecewise):
  214. # Here we modify Piecewise so each expression is now
  215. # an Assignment, and then continue on the print.
  216. expressions = []
  217. conditions = []
  218. for (e, c) in rhs.args:
  219. expressions.append(Assignment(lhs, e))
  220. conditions.append(c)
  221. temp = Piecewise(*zip(expressions, conditions))
  222. return self._print(temp)
  223. if self._settings["contract"] and (lhs.has(IndexedBase) or
  224. rhs.has(IndexedBase)):
  225. # Here we check if there is looping to be done, and if so
  226. # print the required loops.
  227. return self._doprint_loops(rhs, lhs)
  228. else:
  229. lhs_code = self._print(lhs)
  230. rhs_code = self._print(rhs)
  231. return self._get_statement("%s = %s" % (lhs_code, rhs_code))
  232. def _print_Infinity(self, expr):
  233. return 'inf'
  234. def _print_NegativeInfinity(self, expr):
  235. return '-inf'
  236. def _print_NaN(self, expr):
  237. return 'NaN'
  238. def _print_list(self, expr):
  239. return '{' + ', '.join(self._print(a) for a in expr) + '}'
  240. _print_tuple = _print_list
  241. _print_Tuple = _print_list
  242. _print_List = _print_list
  243. def _print_BooleanTrue(self, expr):
  244. return "true"
  245. def _print_BooleanFalse(self, expr):
  246. return "false"
  247. def _print_bool(self, expr):
  248. return str(expr).lower()
  249. # Could generate quadrature code for definite Integrals?
  250. #_print_Integral = _print_not_supported
  251. def _print_MatrixBase(self, A):
  252. # Handle zero dimensions:
  253. if (A.rows, A.cols) == (0, 0):
  254. return '[]'
  255. elif S.Zero in A.shape:
  256. return 'zeros(%s, %s)' % (A.rows, A.cols)
  257. elif (A.rows, A.cols) == (1, 1):
  258. # Octave does not distinguish between scalars and 1x1 matrices
  259. return self._print(A[0, 0])
  260. return "[%s]" % "; ".join(" ".join([self._print(a) for a in A[r, :]])
  261. for r in range(A.rows))
  262. def _print_SparseRepMatrix(self, A):
  263. from sympy.matrices import Matrix
  264. L = A.col_list()
  265. # make row vectors of the indices and entries
  266. I = Matrix([[k[0] + 1 for k in L]])
  267. J = Matrix([[k[1] + 1 for k in L]])
  268. AIJ = Matrix([[k[2] for k in L]])
  269. return "sparse(%s, %s, %s, %s, %s)" % (self._print(I), self._print(J),
  270. self._print(AIJ), A.rows, A.cols)
  271. def _print_MatrixElement(self, expr):
  272. return self.parenthesize(expr.parent, PRECEDENCE["Atom"], strict=True) \
  273. + '(%s, %s)' % (expr.i + 1, expr.j + 1)
  274. def _print_MatrixSlice(self, expr):
  275. def strslice(x, lim):
  276. l = x[0] + 1
  277. h = x[1]
  278. step = x[2]
  279. lstr = self._print(l)
  280. hstr = 'end' if h == lim else self._print(h)
  281. if step == 1:
  282. if l == 1 and h == lim:
  283. return ':'
  284. if l == h:
  285. return lstr
  286. else:
  287. return lstr + ':' + hstr
  288. else:
  289. return ':'.join((lstr, self._print(step), hstr))
  290. return (self._print(expr.parent) + '(' +
  291. strslice(expr.rowslice, expr.parent.shape[0]) + ', ' +
  292. strslice(expr.colslice, expr.parent.shape[1]) + ')')
  293. def _print_Indexed(self, expr):
  294. inds = [ self._print(i) for i in expr.indices ]
  295. return "%s(%s)" % (self._print(expr.base.label), ", ".join(inds))
  296. def _print_KroneckerDelta(self, expr):
  297. prec = PRECEDENCE["Pow"]
  298. return "double(%s == %s)" % tuple(self.parenthesize(x, prec)
  299. for x in expr.args)
  300. def _print_HadamardProduct(self, expr):
  301. return '.*'.join([self.parenthesize(arg, precedence(expr))
  302. for arg in expr.args])
  303. def _print_HadamardPower(self, expr):
  304. PREC = precedence(expr)
  305. return '.**'.join([
  306. self.parenthesize(expr.base, PREC),
  307. self.parenthesize(expr.exp, PREC)
  308. ])
  309. def _print_Identity(self, expr):
  310. shape = expr.shape
  311. if len(shape) == 2 and shape[0] == shape[1]:
  312. shape = [shape[0]]
  313. s = ", ".join(self._print(n) for n in shape)
  314. return "eye(" + s + ")"
  315. def _print_lowergamma(self, expr):
  316. # Octave implements regularized incomplete gamma function
  317. return "(gammainc({1}, {0}).*gamma({0}))".format(
  318. self._print(expr.args[0]), self._print(expr.args[1]))
  319. def _print_uppergamma(self, expr):
  320. return "(gammainc({1}, {0}, 'upper').*gamma({0}))".format(
  321. self._print(expr.args[0]), self._print(expr.args[1]))
  322. def _print_sinc(self, expr):
  323. #Note: Divide by pi because Octave implements normalized sinc function.
  324. return "sinc(%s)" % self._print(expr.args[0]/S.Pi)
  325. def _print_hankel1(self, expr):
  326. return "besselh(%s, 1, %s)" % (self._print(expr.order),
  327. self._print(expr.argument))
  328. def _print_hankel2(self, expr):
  329. return "besselh(%s, 2, %s)" % (self._print(expr.order),
  330. self._print(expr.argument))
  331. # Note: as of 2015, Octave doesn't have spherical Bessel functions
  332. def _print_jn(self, expr):
  333. from sympy.functions import sqrt, besselj
  334. x = expr.argument
  335. expr2 = sqrt(S.Pi/(2*x))*besselj(expr.order + S.Half, x)
  336. return self._print(expr2)
  337. def _print_yn(self, expr):
  338. from sympy.functions import sqrt, bessely
  339. x = expr.argument
  340. expr2 = sqrt(S.Pi/(2*x))*bessely(expr.order + S.Half, x)
  341. return self._print(expr2)
  342. def _print_airyai(self, expr):
  343. return "airy(0, %s)" % self._print(expr.args[0])
  344. def _print_airyaiprime(self, expr):
  345. return "airy(1, %s)" % self._print(expr.args[0])
  346. def _print_airybi(self, expr):
  347. return "airy(2, %s)" % self._print(expr.args[0])
  348. def _print_airybiprime(self, expr):
  349. return "airy(3, %s)" % self._print(expr.args[0])
  350. def _print_expint(self, expr):
  351. mu, x = expr.args
  352. if mu != 1:
  353. return self._print_not_supported(expr)
  354. return "expint(%s)" % self._print(x)
  355. def _one_or_two_reversed_args(self, expr):
  356. assert len(expr.args) <= 2
  357. return '{name}({args})'.format(
  358. name=self.known_functions[expr.__class__.__name__],
  359. args=", ".join([self._print(x) for x in reversed(expr.args)])
  360. )
  361. _print_DiracDelta = _print_LambertW = _one_or_two_reversed_args
  362. def _nested_binary_math_func(self, expr):
  363. return '{name}({arg1}, {arg2})'.format(
  364. name=self.known_functions[expr.__class__.__name__],
  365. arg1=self._print(expr.args[0]),
  366. arg2=self._print(expr.func(*expr.args[1:]))
  367. )
  368. _print_Max = _print_Min = _nested_binary_math_func
  369. def _print_Piecewise(self, expr):
  370. if expr.args[-1].cond != True:
  371. # We need the last conditional to be a True, otherwise the resulting
  372. # function may not return a result.
  373. raise ValueError("All Piecewise expressions must contain an "
  374. "(expr, True) statement to be used as a default "
  375. "condition. Without one, the generated "
  376. "expression may not evaluate to anything under "
  377. "some condition.")
  378. lines = []
  379. if self._settings["inline"]:
  380. # Express each (cond, expr) pair in a nested Horner form:
  381. # (condition) .* (expr) + (not cond) .* (<others>)
  382. # Expressions that result in multiple statements won't work here.
  383. ecpairs = ["({0}).*({1}) + (~({0})).*(".format
  384. (self._print(c), self._print(e))
  385. for e, c in expr.args[:-1]]
  386. elast = "%s" % self._print(expr.args[-1].expr)
  387. pw = " ...\n".join(ecpairs) + elast + ")"*len(ecpairs)
  388. # Note: current need these outer brackets for 2*pw. Would be
  389. # nicer to teach parenthesize() to do this for us when needed!
  390. return "(" + pw + ")"
  391. else:
  392. for i, (e, c) in enumerate(expr.args):
  393. if i == 0:
  394. lines.append("if (%s)" % self._print(c))
  395. elif i == len(expr.args) - 1 and c == True:
  396. lines.append("else")
  397. else:
  398. lines.append("elseif (%s)" % self._print(c))
  399. code0 = self._print(e)
  400. lines.append(code0)
  401. if i == len(expr.args) - 1:
  402. lines.append("end")
  403. return "\n".join(lines)
  404. def _print_zeta(self, expr):
  405. if len(expr.args) == 1:
  406. return "zeta(%s)" % self._print(expr.args[0])
  407. else:
  408. # Matlab two argument zeta is not equivalent to SymPy's
  409. return self._print_not_supported(expr)
  410. def indent_code(self, code):
  411. """Accepts a string of code or a list of code lines"""
  412. # code mostly copied from ccode
  413. if isinstance(code, str):
  414. code_lines = self.indent_code(code.splitlines(True))
  415. return ''.join(code_lines)
  416. tab = " "
  417. inc_regex = ('^function ', '^if ', '^elseif ', '^else$', '^for ')
  418. dec_regex = ('^end$', '^elseif ', '^else$')
  419. # pre-strip left-space from the code
  420. code = [ line.lstrip(' \t') for line in code ]
  421. increase = [ int(any(search(re, line) for re in inc_regex))
  422. for line in code ]
  423. decrease = [ int(any(search(re, line) for re in dec_regex))
  424. for line in code ]
  425. pretty = []
  426. level = 0
  427. for n, line in enumerate(code):
  428. if line in ('', '\n'):
  429. pretty.append(line)
  430. continue
  431. level -= decrease[n]
  432. pretty.append("%s%s" % (tab*level, line))
  433. level += increase[n]
  434. return pretty
  435. def octave_code(expr, assign_to=None, **settings):
  436. r"""Converts `expr` to a string of Octave (or Matlab) code.
  437. The string uses a subset of the Octave language for Matlab compatibility.
  438. Parameters
  439. ==========
  440. expr : Expr
  441. A SymPy expression to be converted.
  442. assign_to : optional
  443. When given, the argument is used as the name of the variable to which
  444. the expression is assigned. Can be a string, ``Symbol``,
  445. ``MatrixSymbol``, or ``Indexed`` type. This can be helpful for
  446. expressions that generate multi-line statements.
  447. precision : integer, optional
  448. The precision for numbers such as pi [default=16].
  449. user_functions : dict, optional
  450. A dictionary where keys are ``FunctionClass`` instances and values are
  451. their string representations. Alternatively, the dictionary value can
  452. be a list of tuples i.e. [(argument_test, cfunction_string)]. See
  453. below for examples.
  454. human : bool, optional
  455. If True, the result is a single string that may contain some constant
  456. declarations for the number symbols. If False, the same information is
  457. returned in a tuple of (symbols_to_declare, not_supported_functions,
  458. code_text). [default=True].
  459. contract: bool, optional
  460. If True, ``Indexed`` instances are assumed to obey tensor contraction
  461. rules and the corresponding nested loops over indices are generated.
  462. Setting contract=False will not generate loops, instead the user is
  463. responsible to provide values for the indices in the code.
  464. [default=True].
  465. inline: bool, optional
  466. If True, we try to create single-statement code instead of multiple
  467. statements. [default=True].
  468. Examples
  469. ========
  470. >>> from sympy import octave_code, symbols, sin, pi
  471. >>> x = symbols('x')
  472. >>> octave_code(sin(x).series(x).removeO())
  473. 'x.^5/120 - x.^3/6 + x'
  474. >>> from sympy import Rational, ceiling
  475. >>> x, y, tau = symbols("x, y, tau")
  476. >>> octave_code((2*tau)**Rational(7, 2))
  477. '8*sqrt(2)*tau.^(7/2)'
  478. Note that element-wise (Hadamard) operations are used by default between
  479. symbols. This is because its very common in Octave to write "vectorized"
  480. code. It is harmless if the values are scalars.
  481. >>> octave_code(sin(pi*x*y), assign_to="s")
  482. 's = sin(pi*x.*y);'
  483. If you need a matrix product "*" or matrix power "^", you can specify the
  484. symbol as a ``MatrixSymbol``.
  485. >>> from sympy import Symbol, MatrixSymbol
  486. >>> n = Symbol('n', integer=True, positive=True)
  487. >>> A = MatrixSymbol('A', n, n)
  488. >>> octave_code(3*pi*A**3)
  489. '(3*pi)*A^3'
  490. This class uses several rules to decide which symbol to use a product.
  491. Pure numbers use "*", Symbols use ".*" and MatrixSymbols use "*".
  492. A HadamardProduct can be used to specify componentwise multiplication ".*"
  493. of two MatrixSymbols. There is currently there is no easy way to specify
  494. scalar symbols, so sometimes the code might have some minor cosmetic
  495. issues. For example, suppose x and y are scalars and A is a Matrix, then
  496. while a human programmer might write "(x^2*y)*A^3", we generate:
  497. >>> octave_code(x**2*y*A**3)
  498. '(x.^2.*y)*A^3'
  499. Matrices are supported using Octave inline notation. When using
  500. ``assign_to`` with matrices, the name can be specified either as a string
  501. or as a ``MatrixSymbol``. The dimensions must align in the latter case.
  502. >>> from sympy import Matrix, MatrixSymbol
  503. >>> mat = Matrix([[x**2, sin(x), ceiling(x)]])
  504. >>> octave_code(mat, assign_to='A')
  505. 'A = [x.^2 sin(x) ceil(x)];'
  506. ``Piecewise`` expressions are implemented with logical masking by default.
  507. Alternatively, you can pass "inline=False" to use if-else conditionals.
  508. Note that if the ``Piecewise`` lacks a default term, represented by
  509. ``(expr, True)`` then an error will be thrown. This is to prevent
  510. generating an expression that may not evaluate to anything.
  511. >>> from sympy import Piecewise
  512. >>> pw = Piecewise((x + 1, x > 0), (x, True))
  513. >>> octave_code(pw, assign_to=tau)
  514. 'tau = ((x > 0).*(x + 1) + (~(x > 0)).*(x));'
  515. Note that any expression that can be generated normally can also exist
  516. inside a Matrix:
  517. >>> mat = Matrix([[x**2, pw, sin(x)]])
  518. >>> octave_code(mat, assign_to='A')
  519. 'A = [x.^2 ((x > 0).*(x + 1) + (~(x > 0)).*(x)) sin(x)];'
  520. Custom printing can be defined for certain types by passing a dictionary of
  521. "type" : "function" to the ``user_functions`` kwarg. Alternatively, the
  522. dictionary value can be a list of tuples i.e., [(argument_test,
  523. cfunction_string)]. This can be used to call a custom Octave function.
  524. >>> from sympy import Function
  525. >>> f = Function('f')
  526. >>> g = Function('g')
  527. >>> custom_functions = {
  528. ... "f": "existing_octave_fcn",
  529. ... "g": [(lambda x: x.is_Matrix, "my_mat_fcn"),
  530. ... (lambda x: not x.is_Matrix, "my_fcn")]
  531. ... }
  532. >>> mat = Matrix([[1, x]])
  533. >>> octave_code(f(x) + g(x) + g(mat), user_functions=custom_functions)
  534. 'existing_octave_fcn(x) + my_fcn(x) + my_mat_fcn([1 x])'
  535. Support for loops is provided through ``Indexed`` types. With
  536. ``contract=True`` these expressions will be turned into loops, whereas
  537. ``contract=False`` will just print the assignment expression that should be
  538. looped over:
  539. >>> from sympy import Eq, IndexedBase, Idx
  540. >>> len_y = 5
  541. >>> y = IndexedBase('y', shape=(len_y,))
  542. >>> t = IndexedBase('t', shape=(len_y,))
  543. >>> Dy = IndexedBase('Dy', shape=(len_y-1,))
  544. >>> i = Idx('i', len_y-1)
  545. >>> e = Eq(Dy[i], (y[i+1]-y[i])/(t[i+1]-t[i]))
  546. >>> octave_code(e.rhs, assign_to=e.lhs, contract=False)
  547. 'Dy(i) = (y(i + 1) - y(i))./(t(i + 1) - t(i));'
  548. """
  549. return OctaveCodePrinter(settings).doprint(expr, assign_to)
  550. def print_octave_code(expr, **settings):
  551. """Prints the Octave (or Matlab) representation of the given expression.
  552. See `octave_code` for the meaning of the optional arguments.
  553. """
  554. print(octave_code(expr, **settings))