autowrap.py 42 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178
  1. """Module for compiling codegen output, and wrap the binary for use in
  2. python.
  3. .. note:: To use the autowrap module it must first be imported
  4. >>> from sympy.utilities.autowrap import autowrap
  5. This module provides a common interface for different external backends, such
  6. as f2py, fwrap, Cython, SWIG(?) etc. (Currently only f2py and Cython are
  7. implemented) The goal is to provide access to compiled binaries of acceptable
  8. performance with a one-button user interface, e.g.,
  9. >>> from sympy.abc import x,y
  10. >>> expr = (x - y)**25
  11. >>> flat = expr.expand()
  12. >>> binary_callable = autowrap(flat)
  13. >>> binary_callable(2, 3)
  14. -1.0
  15. Although a SymPy user might primarily be interested in working with
  16. mathematical expressions and not in the details of wrapping tools
  17. needed to evaluate such expressions efficiently in numerical form,
  18. the user cannot do so without some understanding of the
  19. limits in the target language. For example, the expanded expression
  20. contains large coefficients which result in loss of precision when
  21. computing the expression:
  22. >>> binary_callable(3, 2)
  23. 0.0
  24. >>> binary_callable(4, 5), binary_callable(5, 4)
  25. (-22925376.0, 25165824.0)
  26. Wrapping the unexpanded expression gives the expected behavior:
  27. >>> e = autowrap(expr)
  28. >>> e(4, 5), e(5, 4)
  29. (-1.0, 1.0)
  30. The callable returned from autowrap() is a binary Python function, not a
  31. SymPy object. If it is desired to use the compiled function in symbolic
  32. expressions, it is better to use binary_function() which returns a SymPy
  33. Function object. The binary callable is attached as the _imp_ attribute and
  34. invoked when a numerical evaluation is requested with evalf(), or with
  35. lambdify().
  36. >>> from sympy.utilities.autowrap import binary_function
  37. >>> f = binary_function('f', expr)
  38. >>> 2*f(x, y) + y
  39. y + 2*f(x, y)
  40. >>> (2*f(x, y) + y).evalf(2, subs={x: 1, y:2})
  41. 0.e-110
  42. When is this useful?
  43. 1) For computations on large arrays, Python iterations may be too slow,
  44. and depending on the mathematical expression, it may be difficult to
  45. exploit the advanced index operations provided by NumPy.
  46. 2) For *really* long expressions that will be called repeatedly, the
  47. compiled binary should be significantly faster than SymPy's .evalf()
  48. 3) If you are generating code with the codegen utility in order to use
  49. it in another project, the automatic Python wrappers let you test the
  50. binaries immediately from within SymPy.
  51. 4) To create customized ufuncs for use with numpy arrays.
  52. See *ufuncify*.
  53. When is this module NOT the best approach?
  54. 1) If you are really concerned about speed or memory optimizations,
  55. you will probably get better results by working directly with the
  56. wrapper tools and the low level code. However, the files generated
  57. by this utility may provide a useful starting point and reference
  58. code. Temporary files will be left intact if you supply the keyword
  59. tempdir="path/to/files/".
  60. 2) If the array computation can be handled easily by numpy, and you
  61. do not need the binaries for another project.
  62. """
  63. import sys
  64. import os
  65. import shutil
  66. import tempfile
  67. from pathlib import Path
  68. from subprocess import STDOUT, CalledProcessError, check_output
  69. from string import Template
  70. from warnings import warn
  71. from sympy.core.cache import cacheit
  72. from sympy.core.function import Lambda
  73. from sympy.core.relational import Eq
  74. from sympy.core.symbol import Dummy, Symbol
  75. from sympy.tensor.indexed import Idx, IndexedBase
  76. from sympy.utilities.codegen import (make_routine, get_code_generator,
  77. OutputArgument, InOutArgument,
  78. InputArgument, CodeGenArgumentListError,
  79. Result, ResultBase, C99CodeGen)
  80. from sympy.utilities.iterables import iterable
  81. from sympy.utilities.lambdify import implemented_function
  82. from sympy.utilities.decorator import doctest_depends_on
  83. _doctest_depends_on = {'exe': ('f2py', 'gfortran', 'gcc'),
  84. 'modules': ('numpy',)}
  85. class CodeWrapError(Exception):
  86. pass
  87. class CodeWrapper:
  88. """Base Class for code wrappers"""
  89. _filename = "wrapped_code"
  90. _module_basename = "wrapper_module"
  91. _module_counter = 0
  92. @property
  93. def filename(self):
  94. return "%s_%s" % (self._filename, CodeWrapper._module_counter)
  95. @property
  96. def module_name(self):
  97. return "%s_%s" % (self._module_basename, CodeWrapper._module_counter)
  98. def __init__(self, generator, filepath=None, flags=[], verbose=False):
  99. """
  100. generator -- the code generator to use
  101. """
  102. self.generator = generator
  103. self.filepath = filepath
  104. self.flags = flags
  105. self.quiet = not verbose
  106. @property
  107. def include_header(self):
  108. return bool(self.filepath)
  109. @property
  110. def include_empty(self):
  111. return bool(self.filepath)
  112. def _generate_code(self, main_routine, routines):
  113. routines.append(main_routine)
  114. self.generator.write(
  115. routines, self.filename, True, self.include_header,
  116. self.include_empty)
  117. def wrap_code(self, routine, helpers=None):
  118. helpers = helpers or []
  119. if self.filepath:
  120. workdir = os.path.abspath(self.filepath)
  121. else:
  122. workdir = tempfile.mkdtemp("_sympy_compile")
  123. if not os.access(workdir, os.F_OK):
  124. os.mkdir(workdir)
  125. oldwork = os.getcwd()
  126. os.chdir(workdir)
  127. try:
  128. sys.path.append(workdir)
  129. self._generate_code(routine, helpers)
  130. self._prepare_files(routine)
  131. self._process_files(routine)
  132. mod = __import__(self.module_name)
  133. finally:
  134. sys.path.remove(workdir)
  135. CodeWrapper._module_counter += 1
  136. os.chdir(oldwork)
  137. if not self.filepath:
  138. try:
  139. shutil.rmtree(workdir)
  140. except OSError:
  141. # Could be some issues on Windows
  142. pass
  143. return self._get_wrapped_function(mod, routine.name)
  144. def _process_files(self, routine):
  145. command = self.command
  146. command.extend(self.flags)
  147. try:
  148. retoutput = check_output(command, stderr=STDOUT)
  149. except CalledProcessError as e:
  150. raise CodeWrapError(
  151. "Error while executing command: %s. Command output is:\n%s" % (
  152. " ".join(command), e.output.decode('utf-8')))
  153. if not self.quiet:
  154. print(retoutput)
  155. class DummyWrapper(CodeWrapper):
  156. """Class used for testing independent of backends """
  157. template = """# dummy module for testing of SymPy
  158. def %(name)s():
  159. return "%(expr)s"
  160. %(name)s.args = "%(args)s"
  161. %(name)s.returns = "%(retvals)s"
  162. """
  163. def _prepare_files(self, routine):
  164. return
  165. def _generate_code(self, routine, helpers):
  166. with open('%s.py' % self.module_name, 'w') as f:
  167. printed = ", ".join(
  168. [str(res.expr) for res in routine.result_variables])
  169. # convert OutputArguments to return value like f2py
  170. args = filter(lambda x: not isinstance(
  171. x, OutputArgument), routine.arguments)
  172. retvals = []
  173. for val in routine.result_variables:
  174. if isinstance(val, Result):
  175. retvals.append('nameless')
  176. else:
  177. retvals.append(val.result_var)
  178. print(DummyWrapper.template % {
  179. 'name': routine.name,
  180. 'expr': printed,
  181. 'args': ", ".join([str(a.name) for a in args]),
  182. 'retvals': ", ".join([str(val) for val in retvals])
  183. }, end="", file=f)
  184. def _process_files(self, routine):
  185. return
  186. @classmethod
  187. def _get_wrapped_function(cls, mod, name):
  188. return getattr(mod, name)
  189. class CythonCodeWrapper(CodeWrapper):
  190. """Wrapper that uses Cython"""
  191. setup_template = """\
  192. from setuptools import setup
  193. from setuptools import Extension
  194. from Cython.Build import cythonize
  195. cy_opts = {cythonize_options}
  196. {np_import}
  197. ext_mods = [Extension(
  198. {ext_args},
  199. include_dirs={include_dirs},
  200. library_dirs={library_dirs},
  201. libraries={libraries},
  202. extra_compile_args={extra_compile_args},
  203. extra_link_args={extra_link_args}
  204. )]
  205. setup(ext_modules=cythonize(ext_mods, **cy_opts))
  206. """
  207. _cythonize_options = {'compiler_directives':{'language_level' : "3"}}
  208. pyx_imports = (
  209. "import numpy as np\n"
  210. "cimport numpy as np\n\n")
  211. pyx_header = (
  212. "cdef extern from '{header_file}.h':\n"
  213. " {prototype}\n\n")
  214. pyx_func = (
  215. "def {name}_c({arg_string}):\n"
  216. "\n"
  217. "{declarations}"
  218. "{body}")
  219. std_compile_flag = '-std=c99'
  220. def __init__(self, *args, **kwargs):
  221. """Instantiates a Cython code wrapper.
  222. The following optional parameters get passed to ``setuptools.Extension``
  223. for building the Python extension module. Read its documentation to
  224. learn more.
  225. Parameters
  226. ==========
  227. include_dirs : [list of strings]
  228. A list of directories to search for C/C++ header files (in Unix
  229. form for portability).
  230. library_dirs : [list of strings]
  231. A list of directories to search for C/C++ libraries at link time.
  232. libraries : [list of strings]
  233. A list of library names (not filenames or paths) to link against.
  234. extra_compile_args : [list of strings]
  235. Any extra platform- and compiler-specific information to use when
  236. compiling the source files in 'sources'. For platforms and
  237. compilers where "command line" makes sense, this is typically a
  238. list of command-line arguments, but for other platforms it could be
  239. anything. Note that the attribute ``std_compile_flag`` will be
  240. appended to this list.
  241. extra_link_args : [list of strings]
  242. Any extra platform- and compiler-specific information to use when
  243. linking object files together to create the extension (or to create
  244. a new static Python interpreter). Similar interpretation as for
  245. 'extra_compile_args'.
  246. cythonize_options : [dictionary]
  247. Keyword arguments passed on to cythonize.
  248. """
  249. self._include_dirs = kwargs.pop('include_dirs', [])
  250. self._library_dirs = kwargs.pop('library_dirs', [])
  251. self._libraries = kwargs.pop('libraries', [])
  252. self._extra_compile_args = kwargs.pop('extra_compile_args', [])
  253. self._extra_compile_args.append(self.std_compile_flag)
  254. self._extra_link_args = kwargs.pop('extra_link_args', [])
  255. self._cythonize_options = kwargs.pop('cythonize_options', self._cythonize_options)
  256. self._need_numpy = False
  257. super().__init__(*args, **kwargs)
  258. @property
  259. def command(self):
  260. command = [sys.executable, "setup.py", "build_ext", "--inplace"]
  261. return command
  262. def _prepare_files(self, routine, build_dir=os.curdir):
  263. # NOTE : build_dir is used for testing purposes.
  264. pyxfilename = self.module_name + '.pyx'
  265. codefilename = "%s.%s" % (self.filename, self.generator.code_extension)
  266. # pyx
  267. with open(os.path.join(build_dir, pyxfilename), 'w') as f:
  268. self.dump_pyx([routine], f, self.filename)
  269. # setup.py
  270. ext_args = [repr(self.module_name), repr([pyxfilename, codefilename])]
  271. if self._need_numpy:
  272. np_import = 'import numpy as np\n'
  273. self._include_dirs.append('np.get_include()')
  274. else:
  275. np_import = ''
  276. includes = str(self._include_dirs).replace("'np.get_include()'",
  277. 'np.get_include()')
  278. code = self.setup_template.format(
  279. ext_args=", ".join(ext_args),
  280. np_import=np_import,
  281. include_dirs=includes,
  282. library_dirs=self._library_dirs,
  283. libraries=self._libraries,
  284. extra_compile_args=self._extra_compile_args,
  285. extra_link_args=self._extra_link_args,
  286. cythonize_options=self._cythonize_options)
  287. Path(os.path.join(build_dir, 'setup.py')).write_text(code)
  288. @classmethod
  289. def _get_wrapped_function(cls, mod, name):
  290. return getattr(mod, name + '_c')
  291. def dump_pyx(self, routines, f, prefix):
  292. """Write a Cython file with Python wrappers
  293. This file contains all the definitions of the routines in c code and
  294. refers to the header file.
  295. Arguments
  296. ---------
  297. routines
  298. List of Routine instances
  299. f
  300. File-like object to write the file to
  301. prefix
  302. The filename prefix, used to refer to the proper header file.
  303. Only the basename of the prefix is used.
  304. """
  305. headers = []
  306. functions = []
  307. for routine in routines:
  308. prototype = self.generator.get_prototype(routine)
  309. # C Function Header Import
  310. headers.append(self.pyx_header.format(header_file=prefix,
  311. prototype=prototype))
  312. # Partition the C function arguments into categories
  313. py_rets, py_args, py_loc, py_inf = self._partition_args(routine.arguments)
  314. # Function prototype
  315. name = routine.name
  316. arg_string = ", ".join(self._prototype_arg(arg) for arg in py_args)
  317. # Local Declarations
  318. local_decs = []
  319. for arg, val in py_inf.items():
  320. proto = self._prototype_arg(arg)
  321. mat, ind = [self._string_var(v) for v in val]
  322. local_decs.append(" cdef {} = {}.shape[{}]".format(proto, mat, ind))
  323. local_decs.extend([" cdef {}".format(self._declare_arg(a)) for a in py_loc])
  324. declarations = "\n".join(local_decs)
  325. if declarations:
  326. declarations = declarations + "\n"
  327. # Function Body
  328. args_c = ", ".join([self._call_arg(a) for a in routine.arguments])
  329. rets = ", ".join([self._string_var(r.name) for r in py_rets])
  330. if routine.results:
  331. body = ' return %s(%s)' % (routine.name, args_c)
  332. if rets:
  333. body = body + ', ' + rets
  334. else:
  335. body = ' %s(%s)\n' % (routine.name, args_c)
  336. body = body + ' return ' + rets
  337. functions.append(self.pyx_func.format(name=name, arg_string=arg_string,
  338. declarations=declarations, body=body))
  339. # Write text to file
  340. if self._need_numpy:
  341. # Only import numpy if required
  342. f.write(self.pyx_imports)
  343. f.write('\n'.join(headers))
  344. f.write('\n'.join(functions))
  345. def _partition_args(self, args):
  346. """Group function arguments into categories."""
  347. py_args = []
  348. py_returns = []
  349. py_locals = []
  350. py_inferred = {}
  351. for arg in args:
  352. if isinstance(arg, OutputArgument):
  353. py_returns.append(arg)
  354. py_locals.append(arg)
  355. elif isinstance(arg, InOutArgument):
  356. py_returns.append(arg)
  357. py_args.append(arg)
  358. else:
  359. py_args.append(arg)
  360. # Find arguments that are array dimensions. These can be inferred
  361. # locally in the Cython code.
  362. if isinstance(arg, (InputArgument, InOutArgument)) and arg.dimensions:
  363. dims = [d[1] + 1 for d in arg.dimensions]
  364. sym_dims = [(i, d) for (i, d) in enumerate(dims) if
  365. isinstance(d, Symbol)]
  366. for (i, d) in sym_dims:
  367. py_inferred[d] = (arg.name, i)
  368. for arg in args:
  369. if arg.name in py_inferred:
  370. py_inferred[arg] = py_inferred.pop(arg.name)
  371. # Filter inferred arguments from py_args
  372. py_args = [a for a in py_args if a not in py_inferred]
  373. return py_returns, py_args, py_locals, py_inferred
  374. def _prototype_arg(self, arg):
  375. mat_dec = "np.ndarray[{mtype}, ndim={ndim}] {name}"
  376. np_types = {'double': 'np.double_t',
  377. 'int': 'np.int_t'}
  378. t = arg.get_datatype('c')
  379. if arg.dimensions:
  380. self._need_numpy = True
  381. ndim = len(arg.dimensions)
  382. mtype = np_types[t]
  383. return mat_dec.format(mtype=mtype, ndim=ndim, name=self._string_var(arg.name))
  384. else:
  385. return "%s %s" % (t, self._string_var(arg.name))
  386. def _declare_arg(self, arg):
  387. proto = self._prototype_arg(arg)
  388. if arg.dimensions:
  389. shape = '(' + ','.join(self._string_var(i[1] + 1) for i in arg.dimensions) + ')'
  390. return proto + " = np.empty({shape})".format(shape=shape)
  391. else:
  392. return proto + " = 0"
  393. def _call_arg(self, arg):
  394. if arg.dimensions:
  395. t = arg.get_datatype('c')
  396. return "<{}*> {}.data".format(t, self._string_var(arg.name))
  397. elif isinstance(arg, ResultBase):
  398. return "&{}".format(self._string_var(arg.name))
  399. else:
  400. return self._string_var(arg.name)
  401. def _string_var(self, var):
  402. printer = self.generator.printer.doprint
  403. return printer(var)
  404. class F2PyCodeWrapper(CodeWrapper):
  405. """Wrapper that uses f2py"""
  406. def __init__(self, *args, **kwargs):
  407. ext_keys = ['include_dirs', 'library_dirs', 'libraries',
  408. 'extra_compile_args', 'extra_link_args']
  409. msg = ('The compilation option kwarg {} is not supported with the f2py '
  410. 'backend.')
  411. for k in ext_keys:
  412. if k in kwargs.keys():
  413. warn(msg.format(k))
  414. kwargs.pop(k, None)
  415. super().__init__(*args, **kwargs)
  416. @property
  417. def command(self):
  418. filename = self.filename + '.' + self.generator.code_extension
  419. args = ['-c', '-m', self.module_name, filename]
  420. command = [sys.executable, "-c", "import numpy.f2py as f2py2e;f2py2e.main()"]+args
  421. return command
  422. def _prepare_files(self, routine):
  423. pass
  424. @classmethod
  425. def _get_wrapped_function(cls, mod, name):
  426. return getattr(mod, name)
  427. # Here we define a lookup of backends -> tuples of languages. For now, each
  428. # tuple is of length 1, but if a backend supports more than one language,
  429. # the most preferable language is listed first.
  430. _lang_lookup = {'CYTHON': ('C99', 'C89', 'C'),
  431. 'F2PY': ('F95',),
  432. 'NUMPY': ('C99', 'C89', 'C'),
  433. 'DUMMY': ('F95',)} # Dummy here just for testing
  434. def _infer_language(backend):
  435. """For a given backend, return the top choice of language"""
  436. langs = _lang_lookup.get(backend.upper(), False)
  437. if not langs:
  438. raise ValueError("Unrecognized backend: " + backend)
  439. return langs[0]
  440. def _validate_backend_language(backend, language):
  441. """Throws error if backend and language are incompatible"""
  442. langs = _lang_lookup.get(backend.upper(), False)
  443. if not langs:
  444. raise ValueError("Unrecognized backend: " + backend)
  445. if language.upper() not in langs:
  446. raise ValueError(("Backend {} and language {} are "
  447. "incompatible").format(backend, language))
  448. @cacheit
  449. @doctest_depends_on(exe=('f2py', 'gfortran'), modules=('numpy',))
  450. def autowrap(expr, language=None, backend='f2py', tempdir=None, args=None,
  451. flags=None, verbose=False, helpers=None, code_gen=None, **kwargs):
  452. """Generates Python callable binaries based on the math expression.
  453. Parameters
  454. ==========
  455. expr
  456. The SymPy expression that should be wrapped as a binary routine.
  457. language : string, optional
  458. If supplied, (options: 'C' or 'F95'), specifies the language of the
  459. generated code. If ``None`` [default], the language is inferred based
  460. upon the specified backend.
  461. backend : string, optional
  462. Backend used to wrap the generated code. Either 'f2py' [default],
  463. or 'cython'.
  464. tempdir : string, optional
  465. Path to directory for temporary files. If this argument is supplied,
  466. the generated code and the wrapper input files are left intact in the
  467. specified path.
  468. args : iterable, optional
  469. An ordered iterable of symbols. Specifies the argument sequence for the
  470. function.
  471. flags : iterable, optional
  472. Additional option flags that will be passed to the backend.
  473. verbose : bool, optional
  474. If True, autowrap will not mute the command line backends. This can be
  475. helpful for debugging.
  476. helpers : 3-tuple or iterable of 3-tuples, optional
  477. Used to define auxiliary functions needed for the main expression.
  478. Each tuple should be of the form (name, expr, args) where:
  479. - name : str, the function name
  480. - expr : sympy expression, the function
  481. - args : iterable, the function arguments (can be any iterable of symbols)
  482. code_gen : CodeGen instance
  483. An instance of a CodeGen subclass. Overrides ``language``.
  484. include_dirs : [string]
  485. A list of directories to search for C/C++ header files (in Unix form
  486. for portability).
  487. library_dirs : [string]
  488. A list of directories to search for C/C++ libraries at link time.
  489. libraries : [string]
  490. A list of library names (not filenames or paths) to link against.
  491. extra_compile_args : [string]
  492. Any extra platform- and compiler-specific information to use when
  493. compiling the source files in 'sources'. For platforms and compilers
  494. where "command line" makes sense, this is typically a list of
  495. command-line arguments, but for other platforms it could be anything.
  496. extra_link_args : [string]
  497. Any extra platform- and compiler-specific information to use when
  498. linking object files together to create the extension (or to create a
  499. new static Python interpreter). Similar interpretation as for
  500. 'extra_compile_args'.
  501. Examples
  502. ========
  503. Basic usage:
  504. >>> from sympy.abc import x, y, z
  505. >>> from sympy.utilities.autowrap import autowrap
  506. >>> expr = ((x - y + z)**(13)).expand()
  507. >>> binary_func = autowrap(expr)
  508. >>> binary_func(1, 4, 2)
  509. -1.0
  510. Using helper functions:
  511. >>> from sympy.abc import x, t
  512. >>> from sympy import Function
  513. >>> helper_func = Function('helper_func') # Define symbolic function
  514. >>> expr = 3*x + helper_func(t) # Main expression using helper function
  515. >>> # Define helper_func(x) = 4*x using f2py backend
  516. >>> binary_func = autowrap(expr, args=[x, t],
  517. ... helpers=('helper_func', 4*x, [x]))
  518. >>> binary_func(2, 5) # 3*2 + helper_func(5) = 6 + 20
  519. 26.0
  520. >>> # Same example using cython backend
  521. >>> binary_func = autowrap(expr, args=[x, t], backend='cython',
  522. ... helpers=[('helper_func', 4*x, [x])])
  523. >>> binary_func(2, 5) # 3*2 + helper_func(5) = 6 + 20
  524. 26.0
  525. Type handling example:
  526. >>> import numpy as np
  527. >>> expr = x + y
  528. >>> f_cython = autowrap(expr, backend='cython')
  529. >>> f_cython(1, 2) # doctest: +ELLIPSIS
  530. Traceback (most recent call last):
  531. ...
  532. TypeError: Argument '_x' has incorrect type (expected numpy.ndarray, got int)
  533. >>> f_cython(np.array([1.0]), np.array([2.0]))
  534. array([ 3.])
  535. """
  536. if language:
  537. if not isinstance(language, type):
  538. _validate_backend_language(backend, language)
  539. else:
  540. language = _infer_language(backend)
  541. # two cases 1) helpers is an iterable of 3-tuples and 2) helpers is a
  542. # 3-tuple
  543. if iterable(helpers) and len(helpers) != 0 and iterable(helpers[0]):
  544. helpers = helpers if helpers else ()
  545. else:
  546. helpers = [helpers] if helpers else ()
  547. args = list(args) if iterable(args, exclude=set) else args
  548. if code_gen is None:
  549. code_gen = get_code_generator(language, "autowrap")
  550. CodeWrapperClass = {
  551. 'F2PY': F2PyCodeWrapper,
  552. 'CYTHON': CythonCodeWrapper,
  553. 'DUMMY': DummyWrapper
  554. }[backend.upper()]
  555. code_wrapper = CodeWrapperClass(code_gen, tempdir, flags if flags else (),
  556. verbose, **kwargs)
  557. helps = []
  558. for name_h, expr_h, args_h in helpers:
  559. helps.append(code_gen.routine(name_h, expr_h, args_h))
  560. for name_h, expr_h, args_h in helpers:
  561. if expr.has(expr_h):
  562. name_h = binary_function(name_h, expr_h, backend='dummy')
  563. expr = expr.subs(expr_h, name_h(*args_h))
  564. try:
  565. routine = code_gen.routine('autofunc', expr, args)
  566. except CodeGenArgumentListError as e:
  567. # if all missing arguments are for pure output, we simply attach them
  568. # at the end and try again, because the wrappers will silently convert
  569. # them to return values anyway.
  570. new_args = []
  571. for missing in e.missing_args:
  572. if not isinstance(missing, OutputArgument):
  573. raise
  574. new_args.append(missing.name)
  575. routine = code_gen.routine('autofunc', expr, args + new_args)
  576. return code_wrapper.wrap_code(routine, helpers=helps)
  577. @doctest_depends_on(exe=('f2py', 'gfortran'), modules=('numpy',))
  578. def binary_function(symfunc, expr, **kwargs):
  579. """Returns a SymPy function with expr as binary implementation
  580. This is a convenience function that automates the steps needed to
  581. autowrap the SymPy expression and attaching it to a Function object
  582. with implemented_function().
  583. Parameters
  584. ==========
  585. symfunc : SymPy Function
  586. The function to bind the callable to.
  587. expr : SymPy Expression
  588. The expression used to generate the function.
  589. kwargs : dict
  590. Any kwargs accepted by autowrap.
  591. Examples
  592. ========
  593. >>> from sympy.abc import x, y
  594. >>> from sympy.utilities.autowrap import binary_function
  595. >>> expr = ((x - y)**(25)).expand()
  596. >>> f = binary_function('f', expr)
  597. >>> type(f)
  598. <class 'sympy.core.function.UndefinedFunction'>
  599. >>> 2*f(x, y)
  600. 2*f(x, y)
  601. >>> f(x, y).evalf(2, subs={x: 1, y: 2})
  602. -1.0
  603. """
  604. binary = autowrap(expr, **kwargs)
  605. return implemented_function(symfunc, binary)
  606. #################################################################
  607. # UFUNCIFY #
  608. #################################################################
  609. _ufunc_top = Template("""\
  610. #include "Python.h"
  611. #include "math.h"
  612. #include "numpy/ndarraytypes.h"
  613. #include "numpy/ufuncobject.h"
  614. #include "numpy/halffloat.h"
  615. #include ${include_file}
  616. static PyMethodDef ${module}Methods[] = {
  617. {NULL, NULL, 0, NULL}
  618. };""")
  619. _ufunc_outcalls = Template("*((double *)out${outnum}) = ${funcname}(${call_args});")
  620. _ufunc_body = Template("""\
  621. #ifdef NPY_1_19_API_VERSION
  622. static void ${funcname}_ufunc(char **args, const npy_intp *dimensions, const npy_intp* steps, void* data)
  623. #else
  624. static void ${funcname}_ufunc(char **args, npy_intp *dimensions, npy_intp* steps, void* data)
  625. #endif
  626. {
  627. npy_intp i;
  628. npy_intp n = dimensions[0];
  629. ${declare_args}
  630. ${declare_steps}
  631. for (i = 0; i < n; i++) {
  632. ${outcalls}
  633. ${step_increments}
  634. }
  635. }
  636. PyUFuncGenericFunction ${funcname}_funcs[1] = {&${funcname}_ufunc};
  637. static char ${funcname}_types[${n_types}] = ${types}
  638. static void *${funcname}_data[1] = {NULL};""")
  639. _ufunc_bottom = Template("""\
  640. #if PY_VERSION_HEX >= 0x03000000
  641. static struct PyModuleDef moduledef = {
  642. PyModuleDef_HEAD_INIT,
  643. "${module}",
  644. NULL,
  645. -1,
  646. ${module}Methods,
  647. NULL,
  648. NULL,
  649. NULL,
  650. NULL
  651. };
  652. PyMODINIT_FUNC PyInit_${module}(void)
  653. {
  654. PyObject *m, *d;
  655. ${function_creation}
  656. m = PyModule_Create(&moduledef);
  657. if (!m) {
  658. return NULL;
  659. }
  660. import_array();
  661. import_umath();
  662. d = PyModule_GetDict(m);
  663. ${ufunc_init}
  664. return m;
  665. }
  666. #else
  667. PyMODINIT_FUNC init${module}(void)
  668. {
  669. PyObject *m, *d;
  670. ${function_creation}
  671. m = Py_InitModule("${module}", ${module}Methods);
  672. if (m == NULL) {
  673. return;
  674. }
  675. import_array();
  676. import_umath();
  677. d = PyModule_GetDict(m);
  678. ${ufunc_init}
  679. }
  680. #endif\
  681. """)
  682. _ufunc_init_form = Template("""\
  683. ufunc${ind} = PyUFunc_FromFuncAndData(${funcname}_funcs, ${funcname}_data, ${funcname}_types, 1, ${n_in}, ${n_out},
  684. PyUFunc_None, "${module}", ${docstring}, 0);
  685. PyDict_SetItemString(d, "${funcname}", ufunc${ind});
  686. Py_DECREF(ufunc${ind});""")
  687. _ufunc_setup = Template("""\
  688. from setuptools.extension import Extension
  689. from setuptools import setup
  690. from numpy import get_include
  691. if __name__ == "__main__":
  692. setup(ext_modules=[
  693. Extension('${module}',
  694. sources=['${module}.c', '${filename}.c'],
  695. include_dirs=[get_include()])])
  696. """)
  697. class UfuncifyCodeWrapper(CodeWrapper):
  698. """Wrapper for Ufuncify"""
  699. def __init__(self, *args, **kwargs):
  700. ext_keys = ['include_dirs', 'library_dirs', 'libraries',
  701. 'extra_compile_args', 'extra_link_args']
  702. msg = ('The compilation option kwarg {} is not supported with the numpy'
  703. ' backend.')
  704. for k in ext_keys:
  705. if k in kwargs.keys():
  706. warn(msg.format(k))
  707. kwargs.pop(k, None)
  708. super().__init__(*args, **kwargs)
  709. @property
  710. def command(self):
  711. command = [sys.executable, "setup.py", "build_ext", "--inplace"]
  712. return command
  713. def wrap_code(self, routines, helpers=None):
  714. # This routine overrides CodeWrapper because we can't assume funcname == routines[0].name
  715. # Therefore we have to break the CodeWrapper private API.
  716. # There isn't an obvious way to extend multi-expr support to
  717. # the other autowrap backends, so we limit this change to ufuncify.
  718. helpers = helpers if helpers is not None else []
  719. # We just need a consistent name
  720. funcname = 'wrapped_' + str(id(routines) + id(helpers))
  721. workdir = self.filepath or tempfile.mkdtemp("_sympy_compile")
  722. if not os.access(workdir, os.F_OK):
  723. os.mkdir(workdir)
  724. oldwork = os.getcwd()
  725. os.chdir(workdir)
  726. try:
  727. sys.path.append(workdir)
  728. self._generate_code(routines, helpers)
  729. self._prepare_files(routines, funcname)
  730. self._process_files(routines)
  731. mod = __import__(self.module_name)
  732. finally:
  733. sys.path.remove(workdir)
  734. CodeWrapper._module_counter += 1
  735. os.chdir(oldwork)
  736. if not self.filepath:
  737. try:
  738. shutil.rmtree(workdir)
  739. except OSError:
  740. # Could be some issues on Windows
  741. pass
  742. return self._get_wrapped_function(mod, funcname)
  743. def _generate_code(self, main_routines, helper_routines):
  744. all_routines = main_routines + helper_routines
  745. self.generator.write(
  746. all_routines, self.filename, True, self.include_header,
  747. self.include_empty)
  748. def _prepare_files(self, routines, funcname):
  749. # C
  750. codefilename = self.module_name + '.c'
  751. with open(codefilename, 'w') as f:
  752. self.dump_c(routines, f, self.filename, funcname=funcname)
  753. # setup.py
  754. with open('setup.py', 'w') as f:
  755. self.dump_setup(f)
  756. @classmethod
  757. def _get_wrapped_function(cls, mod, name):
  758. return getattr(mod, name)
  759. def dump_setup(self, f):
  760. setup = _ufunc_setup.substitute(module=self.module_name,
  761. filename=self.filename)
  762. f.write(setup)
  763. def dump_c(self, routines, f, prefix, funcname=None):
  764. """Write a C file with Python wrappers
  765. This file contains all the definitions of the routines in c code.
  766. Arguments
  767. ---------
  768. routines
  769. List of Routine instances
  770. f
  771. File-like object to write the file to
  772. prefix
  773. The filename prefix, used to name the imported module.
  774. funcname
  775. Name of the main function to be returned.
  776. """
  777. if funcname is None:
  778. if len(routines) == 1:
  779. funcname = routines[0].name
  780. else:
  781. msg = 'funcname must be specified for multiple output routines'
  782. raise ValueError(msg)
  783. functions = []
  784. function_creation = []
  785. ufunc_init = []
  786. module = self.module_name
  787. include_file = "\"{}.h\"".format(prefix)
  788. top = _ufunc_top.substitute(include_file=include_file, module=module)
  789. name = funcname
  790. # Partition the C function arguments into categories
  791. # Here we assume all routines accept the same arguments
  792. r_index = 0
  793. py_in, _ = self._partition_args(routines[0].arguments)
  794. n_in = len(py_in)
  795. n_out = len(routines)
  796. # Declare Args
  797. form = "char *{0}{1} = args[{2}];"
  798. arg_decs = [form.format('in', i, i) for i in range(n_in)]
  799. arg_decs.extend([form.format('out', i, i+n_in) for i in range(n_out)])
  800. declare_args = '\n '.join(arg_decs)
  801. # Declare Steps
  802. form = "npy_intp {0}{1}_step = steps[{2}];"
  803. step_decs = [form.format('in', i, i) for i in range(n_in)]
  804. step_decs.extend([form.format('out', i, i+n_in) for i in range(n_out)])
  805. declare_steps = '\n '.join(step_decs)
  806. # Call Args
  807. form = "*(double *)in{0}"
  808. call_args = ', '.join([form.format(a) for a in range(n_in)])
  809. # Step Increments
  810. form = "{0}{1} += {0}{1}_step;"
  811. step_incs = [form.format('in', i) for i in range(n_in)]
  812. step_incs.extend([form.format('out', i, i) for i in range(n_out)])
  813. step_increments = '\n '.join(step_incs)
  814. # Types
  815. n_types = n_in + n_out
  816. types = "{" + ', '.join(["NPY_DOUBLE"]*n_types) + "};"
  817. # Docstring
  818. docstring = '"Created in SymPy with Ufuncify"'
  819. # Function Creation
  820. function_creation.append("PyObject *ufunc{};".format(r_index))
  821. # Ufunc initialization
  822. init_form = _ufunc_init_form.substitute(module=module,
  823. funcname=name,
  824. docstring=docstring,
  825. n_in=n_in, n_out=n_out,
  826. ind=r_index)
  827. ufunc_init.append(init_form)
  828. outcalls = [_ufunc_outcalls.substitute(
  829. outnum=i, call_args=call_args, funcname=routines[i].name) for i in
  830. range(n_out)]
  831. body = _ufunc_body.substitute(module=module, funcname=name,
  832. declare_args=declare_args,
  833. declare_steps=declare_steps,
  834. call_args=call_args,
  835. step_increments=step_increments,
  836. n_types=n_types, types=types,
  837. outcalls='\n '.join(outcalls))
  838. functions.append(body)
  839. body = '\n\n'.join(functions)
  840. ufunc_init = '\n '.join(ufunc_init)
  841. function_creation = '\n '.join(function_creation)
  842. bottom = _ufunc_bottom.substitute(module=module,
  843. ufunc_init=ufunc_init,
  844. function_creation=function_creation)
  845. text = [top, body, bottom]
  846. f.write('\n\n'.join(text))
  847. def _partition_args(self, args):
  848. """Group function arguments into categories."""
  849. py_in = []
  850. py_out = []
  851. for arg in args:
  852. if isinstance(arg, OutputArgument):
  853. py_out.append(arg)
  854. elif isinstance(arg, InOutArgument):
  855. raise ValueError("Ufuncify doesn't support InOutArguments")
  856. else:
  857. py_in.append(arg)
  858. return py_in, py_out
  859. @cacheit
  860. @doctest_depends_on(exe=('f2py', 'gfortran', 'gcc'), modules=('numpy',))
  861. def ufuncify(args, expr, language=None, backend='numpy', tempdir=None,
  862. flags=None, verbose=False, helpers=None, **kwargs):
  863. """Generates a binary function that supports broadcasting on numpy arrays.
  864. Parameters
  865. ==========
  866. args : iterable
  867. Either a Symbol or an iterable of symbols. Specifies the argument
  868. sequence for the function.
  869. expr
  870. A SymPy expression that defines the element wise operation.
  871. language : string, optional
  872. If supplied, (options: 'C' or 'F95'), specifies the language of the
  873. generated code. If ``None`` [default], the language is inferred based
  874. upon the specified backend.
  875. backend : string, optional
  876. Backend used to wrap the generated code. Either 'numpy' [default],
  877. 'cython', or 'f2py'.
  878. tempdir : string, optional
  879. Path to directory for temporary files. If this argument is supplied,
  880. the generated code and the wrapper input files are left intact in
  881. the specified path.
  882. flags : iterable, optional
  883. Additional option flags that will be passed to the backend.
  884. verbose : bool, optional
  885. If True, autowrap will not mute the command line backends. This can
  886. be helpful for debugging.
  887. helpers : 3-tuple or iterable of 3-tuples, optional
  888. Used to define auxiliary functions needed for the main expression.
  889. Each tuple should be of the form (name, expr, args) where:
  890. - name : str, the function name
  891. - expr : sympy expression, the function
  892. - args : iterable, the function arguments (can be any iterable of symbols)
  893. kwargs : dict
  894. These kwargs will be passed to autowrap if the `f2py` or `cython`
  895. backend is used and ignored if the `numpy` backend is used.
  896. Notes
  897. =====
  898. The default backend ('numpy') will create actual instances of
  899. ``numpy.ufunc``. These support ndimensional broadcasting, and implicit type
  900. conversion. Use of the other backends will result in a "ufunc-like"
  901. function, which requires equal length 1-dimensional arrays for all
  902. arguments, and will not perform any type conversions.
  903. References
  904. ==========
  905. .. [1] https://numpy.org/doc/stable/reference/ufuncs.html
  906. Examples
  907. ========
  908. Basic usage:
  909. >>> from sympy.utilities.autowrap import ufuncify
  910. >>> from sympy.abc import x, y
  911. >>> import numpy as np
  912. >>> f = ufuncify((x, y), y + x**2)
  913. >>> type(f)
  914. <class 'numpy.ufunc'>
  915. >>> f([1, 2, 3], 2)
  916. array([ 3., 6., 11.])
  917. >>> f(np.arange(5), 3)
  918. array([ 3., 4., 7., 12., 19.])
  919. Using helper functions:
  920. >>> from sympy import Function
  921. >>> helper_func = Function('helper_func') # Define symbolic function
  922. >>> expr = x**2 + y*helper_func(x) # Main expression using helper function
  923. >>> # Define helper_func(x) = x**3
  924. >>> f = ufuncify((x, y), expr, helpers=[('helper_func', x**3, [x])])
  925. >>> f([1, 2], [3, 4])
  926. array([ 4., 36.])
  927. Type handling with different backends:
  928. For the 'f2py' and 'cython' backends, inputs are required to be equal length
  929. 1-dimensional arrays. The 'f2py' backend will perform type conversion, but
  930. the Cython backend will error if the inputs are not of the expected type.
  931. >>> f_fortran = ufuncify((x, y), y + x**2, backend='f2py')
  932. >>> f_fortran(1, 2)
  933. array([ 3.])
  934. >>> f_fortran(np.array([1, 2, 3]), np.array([1.0, 2.0, 3.0]))
  935. array([ 2., 6., 12.])
  936. >>> f_cython = ufuncify((x, y), y + x**2, backend='Cython')
  937. >>> f_cython(1, 2) # doctest: +ELLIPSIS
  938. Traceback (most recent call last):
  939. ...
  940. TypeError: Argument '_x' has incorrect type (expected numpy.ndarray, got int)
  941. >>> f_cython(np.array([1.0]), np.array([2.0]))
  942. array([ 3.])
  943. """
  944. if isinstance(args, Symbol):
  945. args = (args,)
  946. else:
  947. args = tuple(args)
  948. if language:
  949. _validate_backend_language(backend, language)
  950. else:
  951. language = _infer_language(backend)
  952. helpers = helpers if helpers else ()
  953. flags = flags if flags else ()
  954. if backend.upper() == 'NUMPY':
  955. # maxargs is set by numpy compile-time constant NPY_MAXARGS
  956. # If a future version of numpy modifies or removes this restriction
  957. # this variable should be changed or removed
  958. maxargs = 32
  959. helps = []
  960. for name, expr, args in helpers:
  961. helps.append(make_routine(name, expr, args))
  962. code_wrapper = UfuncifyCodeWrapper(C99CodeGen("ufuncify"), tempdir,
  963. flags, verbose)
  964. if not isinstance(expr, (list, tuple)):
  965. expr = [expr]
  966. if len(expr) == 0:
  967. raise ValueError('Expression iterable has zero length')
  968. if len(expr) + len(args) > maxargs:
  969. msg = ('Cannot create ufunc with more than {0} total arguments: '
  970. 'got {1} in, {2} out')
  971. raise ValueError(msg.format(maxargs, len(args), len(expr)))
  972. routines = [make_routine('autofunc{}'.format(idx), exprx, args) for
  973. idx, exprx in enumerate(expr)]
  974. return code_wrapper.wrap_code(routines, helpers=helps)
  975. else:
  976. # Dummies are used for all added expressions to prevent name clashes
  977. # within the original expression.
  978. y = IndexedBase(Dummy('y'))
  979. m = Dummy('m', integer=True)
  980. i = Idx(Dummy('i', integer=True), m)
  981. f_dummy = Dummy('f')
  982. f = implemented_function('%s_%d' % (f_dummy.name, f_dummy.dummy_index), Lambda(args, expr))
  983. # For each of the args create an indexed version.
  984. indexed_args = [IndexedBase(Dummy(str(a))) for a in args]
  985. # Order the arguments (out, args, dim)
  986. args = [y] + indexed_args + [m]
  987. args_with_indices = [a[i] for a in indexed_args]
  988. return autowrap(Eq(y[i], f(*args_with_indices)), language, backend,
  989. tempdir, args, flags, verbose, helpers, **kwargs)