| Download
Improved version of the octave interface.
Project: Many Interacting Worlds
Path: octave.py
Views: 349r"""1Interface to GNU Octave23GNU Octave is a free software (GPL) MATLAB-like program with numerical4routines for integrating, solving systems of equations, special5functions, and solving (numerically) differential equations. Please see6http://octave.org/ for more details.78The commands in this section only work if you have the optional9"octave" interpreter installed and available in your PATH. It's not10necessary to install any special Sage packages.1112EXAMPLES::1314sage: octave.eval('2+2') # optional - octave15'ans = 4'1617sage: a = octave(10) # optional - octave18sage: a**10 # optional - octave191e+102021LOG: - creation (William Stein) - ? (David Joyner, 2005-12-18) -22Examples (David Joyner, 2005-01-03)2324Computation of Special Functions25--------------------------------2627Octave implements computation of the following special functions28(see the maxima and gp interfaces for even more special29functions)::3031airy32Airy functions of the first and second kind, and their derivatives.33airy(0,x) = Ai(x), airy(1,x) = Ai'(x), airy(2,x) = Bi(x), airy(3,x) = Bi'(x)34besselj35Bessel functions of the first kind.36bessely37Bessel functions of the second kind.38besseli39Modified Bessel functions of the first kind.40besselk41Modified Bessel functions of the second kind.42besselh43Compute Hankel functions of the first (k = 1) or second (k = 2) kind.44beta45The Beta function,46beta (a, b) = gamma (a) * gamma (b) / gamma (a + b).47betainc48The incomplete Beta function,49erf50The error function,51erfinv52The inverse of the error function.53gamma54The Gamma function,55gammainc56The incomplete gamma function,5758For example,5960::6162sage: octave("airy(3,2)") # optional - octave634.1006864sage: octave("beta(2,2)") # optional - octave650.16666766sage: octave("betainc(0.2,2,2)") # optional - octave670.10468sage: octave("besselh(0,2)") # optional - octave69(0.223891,0.510376)70sage: octave("besselh(0,1)") # optional - octave71(0.765198,0.088257)72sage: octave("besseli(1,2)") # optional - octave731.5906474sage: octave("besselj(1,2)") # optional - octave750.57672576sage: octave("besselk(1,2)") # optional - octave770.13986678sage: octave("erf(0)") # optional - octave79080sage: octave("erf(1)") # optional - octave810.84270182sage: octave("erfinv(0.842)") # optional - octave830.99831584sage: octave("gamma(1.5)") # optional - octave850.88622786sage: octave("gammainc(1.5,1)") # optional - octave870.776878889The Octave interface reads in even very long input (using files) in90a robust manner::9192sage: t = '"%s"'%10^10000 # ten thousand character string.93sage: a = octave.eval(t + ';') # optional - octave, < 1/100th of a second94sage: a = octave(t) # optional - octave9596Note that actually reading ``a`` back out takes forever. This *must*97be fixed as soon as possible, see :trac:`940`.9899Tutorial100--------101102EXAMPLES::103104sage: octave('4+10') # optional - octave10514106sage: octave('date') # optional - octave; random output10718-Oct-2007108sage: octave('5*10 + 6') # optional - octave10956110sage: octave('(6+6)/3') # optional - octave1114112sage: octave('9')^2 # optional - octave11381114sage: a = octave(10); b = octave(20); c = octave(30) # optional - octave115sage: avg = (a+b+c)/3 # optional - octave116sage: avg # optional - octave11720118sage: parent(avg) # optional - octave119Octave120121::122123sage: my_scalar = octave('3.1415') # optional - octave124sage: my_scalar # optional - octave1253.1415126sage: my_vector1 = octave('[1,5,7]') # optional - octave127sage: my_vector1 # optional - octave1281 5 7129sage: my_vector2 = octave('[1;5;7]') # optional - octave130sage: my_vector2 # optional - octave131113251337134sage: my_vector1 * my_vector2 # optional - octave13575136"""137138#*****************************************************************************139# Copyright (C) 2005 William Stein <[email protected]>140#141# Distributed under the terms of the GNU General Public License (GPL)142#143# This code is distributed in the hope that it will be useful,144# but WITHOUT ANY WARRANTY; without even the implied warranty of145# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU146# General Public License for more details.147#148# The full text of the GPL is available at:149#150# http://www.gnu.org/licenses/151#*****************************************************************************152153import os154from expect import Expect, ExpectElement155from sage.misc.misc import verbose156157158class Octave(Expect):159r"""160Interface to the Octave interpreter.161162EXAMPLES::163164sage: octave.eval("a = [ 1, 1, 2; 3, 5, 8; 13, 21, 33 ]") # optional - octave165'a =\n\n 1 1 2\n 3 5 8\n 13 21 33\n\n'166sage: octave.eval("b = [ 1; 3; 13]") # optional - octave167'b =\n\n 1\n 3\n 13\n\n'168sage: octave.eval("c=a \\ b") # solves linear equation: a*c = b # optional - octave; random output169'c =\n\n 1\n 7.21645e-16\n -7.21645e-16\n\n'170sage: octave.eval("c") # optional - octave; random output171'c =\n\n 1\n 7.21645e-16\n -7.21645e-16\n\n'172"""173174def __init__(self, maxread=100, script_subdirectory=None, logfile=None, server=None, server_tmpdir=None,175seed=None):176"""177EXAMPLES::178179sage: octave == loads(dumps(octave))180True181"""182Expect.__init__(self,183name = 'octave',184# We want the prompt sequence to be unique to avoid confusion with syntax error messages containing >>>185prompt = 'octave\:\d+> ',186# We don't want any pagination of output187command = "sage-native-execute octave --no-line-editing --silent --eval 'PS2(PS1());more off' --persist",188maxread = maxread,189server = server,190server_tmpdir = server_tmpdir,191script_subdirectory = script_subdirectory,192restart_on_ctrlc = False,193verbose_start = False,194logfile = logfile,195eval_using_file_cutoff=100)196self._seed = seed197198def set_seed(self, seed=None):199"""200Sets the seed for the random number generator201for this octave interpreter.202203EXAMPLES::204205sage: o = Octave() # optional - octave206sage: o.set_seed(1) # optional - octave2071208sage: [o.rand() for i in range(5)] # optional - octave209[ 0.134364, 0.847434, 0.763775, 0.255069, 0.495435]210"""211if seed is None:212seed = self.rand_seed()213self.eval("rand('state',%d)" % seed)214self._seed = seed215return seed216217def __reduce__(self):218"""219EXAMPLES::220221sage: octave.__reduce__()222(<function reduce_load_Octave at 0x...>, ())223"""224return reduce_load_Octave, tuple([])225226def _read_in_file_command(self, filename):227"""228EXAMPLES::229230sage: filename = tmp_filename()231sage: octave._read_in_file_command(filename)232'source("...");'233"""234return 'source("%s");'%filename235236def _quit_string(self):237"""238EXAMPLES::239240sage: octave._quit_string()241'quit;'242"""243return 'quit;'244245def _install_hints(self):246"""247Returns hints on how to install Octave.248249EXAMPLES::250251sage: print octave._install_hints()252You must get ...253"""254return """255You must get the program "octave" in order to use Octave256from Sage. You can read all about Octave at257http://www.gnu.org/software/octave/258259LINUX / WINDOWS (colinux):260Do apt-get install octave as root on your machine261(or, in Windows, in the colinux console).262263OS X:264* This website has links to binaries for OS X PowerPC265and OS X Intel builds of the latest version of Octave:266http://hpc.sourceforge.net/267Once you get the tarball from there, go to the / directory268and type269tar zxvf octave-intel-bin.tar.gz270to extract it to usr/local/... Make sure /usr/local/bin271is in your PATH. Then type "octave" and verify that272octave starts up.273* Darwin ports and fink have Octave as well.274"""275def _eval_line(self, line, reformat=True, allow_use_file=False,276wait_for_prompt=True, restart_if_needed=False):277"""278EXAMPLES::279280sage: print octave._eval_line('2+2') #optional - octave281ans = 4282"""283if not wait_for_prompt:284return Expect._eval_line(self, line)285if line == '':286return ''287if self._expect is None:288self._start()289if allow_use_file and len(line)>3000:290return self._eval_line_using_file(line)291try:292E = self._expect293# debug294# self._synchronize(cmd='1+%s\n')295verbose("in = '%s'"%line,level=3)296E.sendline(line)297E.expect(self._prompt)298out = E.before299# debug300verbose("out = '%s'"%out,level=3)301except EOF:302if self._quit_string() in line:303return ''304except KeyboardInterrupt:305self._keyboard_interrupt()306if reformat:307if 'syntax error' in out:308raise SyntaxError(out)309out = "\n".join(out.splitlines()[1:])310return out311312def _keyboard_interrupt(self):313print "CntrlC: Interrupting %s..."%self314if self._restart_on_ctrlc:315try:316self._expect.close(force=1)317except pexpect.ExceptionPexpect as msg:318raise pexpect.ExceptionPexpect( "THIS IS A BUG -- PLEASE REPORT. This should never happen.\n" + msg)319self._start()320raise KeyboardInterrupt("Restarting %s (WARNING: all variables defined in previous session are now invalid)"%self)321else:322self._expect.send('\003') # control-c323#self._expect.expect(self._prompt)324#self._expect.expect(self._prompt)325raise KeyboardInterrupt("Ctrl-c pressed while running %s"%self)326327def quit(self, verbose=False):328"""329EXAMPLES::330331sage: o = Octave()332sage: o._start() # optional - octave333sage: o.quit(True) # optional - octave334Exiting spawned Octave process.335"""336# Don't bother, since it just hangs in some cases, and it337# isn't necessary, since octave behaves well with respect338# to signals.339if not self._expect is None:340if verbose:341print "Exiting spawned %s process." % self342return343344def _start(self):345"""346Starts the Octave process.347348EXAMPLES::349350sage: o = Octave() # optional - octave351sage: o.is_running() # optional - octave352False353sage: o._start() # optional - octave354sage: o.is_running() # optional - octave355True356"""357Expect._start(self)358self.eval("page_screen_output=0;")359self.eval("format none;")360# set random seed361self.set_seed(self._seed)362363def _equality_symbol(self):364"""365EXAMPLES::366367sage: octave('0 == 1') # optional - octave3680369sage: octave('1 == 1') # optional - octave3701371"""372return '=='373374def _true_symbol(self):375"""376EXAMPLES::377378sage: octave('1 == 1') # optional - octave3791380"""381return '1'382383def _false_symbol(self):384"""385EXAMPLES::386387sage: octave('0 == 1') # optional - octave3880389"""390return '0'391392def set(self, var, value):393"""394Set the variable ``var`` to the given ``value``.395396EXAMPLES::397398sage: octave.set('x', '2') # optional - octave399sage: octave.get('x') # optional - octave400' 2'401"""402cmd = '%s=%s;'%(var,value)403out = self.eval(cmd)404if out.find("error") != -1 or out.find("Error") != -1:405raise TypeError("Error executing code in Octave\nCODE:\n\t%s\nOctave ERROR:\n\t%s"%(cmd, out))406407def get(self, var):408"""409Get the value of the variable ``var``.410411EXAMPLES::412413sage: octave.set('x', '2') # optional - octave414sage: octave.get('x') # optional - octave415' 2'416"""417s = self.eval('%s'%var)418i = s.find('=')419return s[i+1:]420421def clear(self, var):422"""423Clear the variable named var.424425EXAMPLES::426427sage: octave.set('x', '2') # optional - octave428sage: octave.clear('x') # optional - octave429sage: octave.get('x') # optional - octave430"error: 'x' undefined near line ... column 1"431"""432self.eval('clear %s'%var)433434def console(self):435"""436Spawn a new Octave command-line session.437438This requires that the optional octave program be installed and in439your PATH, but no optional Sage packages need be installed.440441EXAMPLES::442443sage: octave_console() # not tested444GNU Octave, version 2.1.73 (i386-apple-darwin8.5.3).445Copyright (C) 2006 John W. Eaton.446...447octave:1> 2+3448ans = 5449octave:2> [ctl-d]450451Pressing ctrl-d exits the octave console and returns you to Sage.452octave, like Sage, remembers its history from one session to453another.454"""455octave_console()456457def version(self):458"""459Return the version of Octave.460461OUTPUT: string462463EXAMPLES::464465sage: octave.version() # optional - octave; random output depending on version466'2.1.73'467"""468return octave_version()469470def solve_linear_system(self, A, b):471r"""472Use octave to compute a solution x to A\*x = b, as a list.473474INPUT:475476- ``A`` -- mxn matrix A with entries in `\QQ` or `\RR`477478- ``b`` -- m-vector b entries in `\QQ` or `\RR` (resp)479480OUTPUT: A list x (if it exists) which solves M\*x = b481482EXAMPLES::483484sage: M33 = MatrixSpace(QQ,3,3)485sage: A = M33([1,2,3,4,5,6,7,8,0])486sage: V3 = VectorSpace(QQ,3)487sage: b = V3([1,2,3])488sage: octave.solve_linear_system(A,b) # optional - octave (and output is slightly random in low order bits)489[-0.33333299999999999, 0.66666700000000001, -3.5236600000000002e-18]490491AUTHORS:492493- David Joyner and William Stein494"""495m = A.nrows()496if m != len(b):497raise ValueError("dimensions of A and b must be compatible")498from sage.matrix.all import MatrixSpace499from sage.rings.all import QQ500MS = MatrixSpace(QQ,m,1)501b = MS(list(b)) # converted b to a "column vector"502sA = self.sage2octave_matrix_string(A)503sb = self.sage2octave_matrix_string(b)504self.eval("a = " + sA )505self.eval("b = " + sb )506soln = octave.eval("c = a \\ b")507soln = soln.replace("\n\n ","[")508soln = soln.replace("\n\n","]")509soln = soln.replace("\n",",")510sol = soln[3:]511return eval(sol)512513514def sage2octave_matrix_string(self, A):515"""516Return an octave matrix from a Sage matrix.517518INPUT: A Sage matrix with entries in the rationals or reals.519520OUTPUT: A string that evaluates to an Octave matrix.521522EXAMPLES::523524sage: M33 = MatrixSpace(QQ,3,3)525sage: A = M33([1,2,3,4,5,6,7,8,0])526sage: octave.sage2octave_matrix_string(A) # optional - octave527'[1, 2, 3; 4, 5, 6; 7, 8, 0]'528529AUTHORS:530531- David Joyner and William Stein532"""533return str(A.rows()).replace('), (', '; ').replace('(', '').replace(')','')534535def de_system_plot(self, f, ics, trange):536r"""537Plots (using octave's interface to gnuplot) the solution to a538`2\times 2` system of differential equations.539540INPUT:541542543- ``f`` - a pair of strings representing the544differential equations; The independent variable must be called x545and the dependent variable must be called y.546547- ``ics`` - a pair [x0,y0] such that x(t0) = x0, y(t0)548= y0549550- ``trange`` - a pair [t0,t1]551552553OUTPUT: a gnuplot window appears554555EXAMPLES::556557sage: octave.de_system_plot(['x+y','x-y'], [1,-1], [0,2]) # not tested -- does this actually work (on OS X it fails for me -- William Stein, 2007-10)558559This should yield the two plots `(t,x(t)), (t,y(t))` on the560same graph (the `t`-axis is the horizontal axis) of the561system of ODEs562563.. math::564565x' = x+y, x(0) = 1;\qquad y' = x-y, y(0) = -1, \quad\text{for}\quad 0 < t < 2.566"""567eqn1 = f[0].replace('x','x(1)').replace('y','x(2)')568eqn2 = f[1].replace('x','x(1)').replace('y','x(2)')569fcn = "function xdot = f(x,t) xdot(1) = %s; xdot(2) = %s; endfunction"%(eqn1, eqn2)570self.eval(fcn)571x0_eqn = "x0 = [%s; %s]"%(ics[0], ics[1])572self.eval(x0_eqn)573t_eqn = "t = linspace(%s, %s, 200)'"%(trange[0], trange[1])574self.eval(t_eqn)575x_eqn = 'x = lsode("f",x0,t);'576self.eval(x_eqn)577self.eval("plot(t,x)")578579def _object_class(self):580"""581EXAMPLES::582583sage: octave._object_class()584<class 'sage.interfaces.octave.OctaveElement'>585"""586return OctaveElement587588589octave_functions = set()590591def to_complex(octave_string, R):592r"""593Helper function to convert octave complex number594595TESTS::596597sage: from sage.interfaces.octave import to_complex598sage: to_complex('(0,1)', CDF)5991.0*I600sage: to_complex('(1.3231,-0.2)', CDF)6011.3231 - 0.2*I602"""603real, imag = octave_string.strip('() ').split(',')604return R(float(real), float(imag))605606class OctaveElement(ExpectElement):607def _get_sage_ring(self):608r"""609TESTS::610611sage: octave('1')._get_sage_ring() # optional - octave612Real Double Field613sage: octave('I')._get_sage_ring() # optional - octave614Complex Double Field615sage: octave('[]')._get_sage_ring() # optional - octave616Real Double Field617"""618if self.isinteger():619import sage.rings.integer_ring620return sage.rings.integer_ring.ZZ621elif self.isreal():622import sage.rings.real_double623return sage.rings.real_double.RDF624elif self.iscomplex():625import sage.rings.complex_double626return sage.rings.complex_double.CDF627else:628raise TypeError("no Sage ring associated to this element.")629630def __nonzero__(self):631r"""632Test whether this element is nonzero.633634EXAMPLES::635636sage: bool(octave('0')) # optional - octave637False638sage: bool(octave('[]')) # optional - octave639False640sage: bool(octave('[0,0]')) # optional - octave641False642sage: bool(octave('[0,0,0;0,0,0]')) # optional - octave643False644645sage: bool(octave('0.1')) # optional - octave646True647sage: bool(octave('[0,1,0]')) # optional - octave648True649sage: bool(octave('[0,0,-0.1;0,0,0]')) # optional - octave650True651"""652return str(self) != ' [](0x0)' and any(x != '0' for x in str(self).split())653654def _matrix_(self, R=None):655r"""656Return Sage matrix from this octave element.657658EXAMPLES::659660sage: A = octave('[1,2;3,4.5]') # optional - octave661sage: matrix(A) # optional - octave662[1.0 2.0]663[3.0 4.5]664sage: _.base_ring() # optional - octave665Real Double Field666667sage: A = octave('[I,1;-1,0]') # optional - octave668sage: matrix(A) # optional - octave669[1.0*I 1.0]670[ -1.0 0.0]671sage: _.base_ring() # optional - octave672Complex Double Field673674sage: A = octave('[1,2;3,4]') # optional - octave675sage: matrix(ZZ, A) # optional - octave676[1 2]677[3 4]678sage: A = octave('[1,2;3,4.5]') # optional - octave679sage: matrix(RR, A) # optional - octave680[1.00000000000000 2.00000000000000]681[3.00000000000000 4.50000000000000]682"""683oc = self.parent()684if not self.ismatrix():685raise TypeError('not an octave matrix')686if R is None:687R = self._get_sage_ring()688689s = str(self).strip('\n ')690w = [u.strip().split(' ') for u in s.split('\n')]691nrows = len(w)692ncols = len(w[0])693694if self.iscomplex():695w = [[to_complex(x,R) for x in row] for row in w]696697from sage.matrix.all import MatrixSpace698s = str(self).strip()699v = s.split('\n ')700nrows = len(v)701if nrows == 0:702return MatrixSpace(R,0,0)(0)703ncols = len(v[0].split())704M = MatrixSpace(R, nrows, ncols)705v = sum([[x for x in w.split()] for w in v], [])706return M(v)707708def _vector_(self, R=None):709r"""710Return Sage vector from this octave element.711712EXAMPLES::713714sage: A = octave('[1,2,3,4]') # optional - octave715sage: vector(ZZ, A) # optional - octave716(1, 2, 3, 4)717sage: A = octave('[1,2.3,4.5]') # optional - octave718sage: vector(A) # optional - octave719(1.0, 2.3, 4.5)720sage: A = octave('[1,I]') # optional - octave721sage: vector(A) # optional - octave722(1.0, 1.0*I)723"""724oc = self.parent()725if not self.isvector():726raise TypeError('not an octave vector')727if R is None:728R = self._get_sage_ring()729730s = str(self).strip('\n ')731w = s.strip().split(' ')732nrows = len(w)733734if self.iscomplex():735w = [to_complex(x, R) for x in w]736737from sage.modules.free_module import FreeModule738return FreeModule(R, nrows)(w)739740def _scalar_(self):741"""742Return Sage scalar from this octave element.743744EXAMPLES::745746sage: A = octave('2833') # optional - octave747sage: As = A.sage(); As # optional - octave7482833.0749sage: As.parent() # optional - octave750Real Double Field751752sage: B = sqrt(A) # optional - octave753sage: Bs = B.sage(); Bs # optional - octave75453.2259755sage: Bs.parent() # optional - octave756Real Double Field757758sage: C = sqrt(-A) # optional - octave759sage: Cs = C.sage(); Cs # optional - octave76053.2259*I761sage: Cs.parent() # optional - octave762Complex Double Field763"""764if not self.isscalar():765raise TypeError("not an octave scalar")766767R = self._get_sage_ring()768if self.iscomplex():769return to_complex(str(self), R)770else:771return R(str(self))772773def _sage_(self):774"""775Try to parse the octave object and return a sage object.776777EXAMPLES::778779sage: A = octave('2833') # optional - octave780sage: A.sage() # optional - octave7812833.0782sage: B = sqrt(A) # optional - octave783sage: B.sage() # optional - octave78453.2259785sage: C = sqrt(-A) # optional - octave786sage: C.sage() # optional - octave78753.2259*I788sage: A = octave('[1,2,3,4]') # optional - octave789sage: A.sage() # optional - octave790(1.0, 2.0, 3.0, 4.0)791sage: A = octave('[1,2.3,4.5]') # optional - octave792sage: A.sage() # optional - octave793(1.0, 2.3, 4.5)794sage: A = octave('[1,2.3+I,4.5]') # optional - octave795sage: A.sage() # optional - octave796(1.0, 2.3 + 1.0*I, 4.5)797"""798if self.isscalar():799return self._scalar_()800elif self.isvector():801return self._vector_()802elif self.ismatrix():803return self._matrix_()804else:805raise NotImplementedError('octave type is not recognized')806807# An instance808octave = Octave()809810def reduce_load_Octave():811"""812EXAMPLES::813814sage: from sage.interfaces.octave import reduce_load_Octave815sage: reduce_load_Octave()816Octave817"""818return octave819820821def octave_console():822"""823Spawn a new Octave command-line session.824825This requires that the optional octave program be installed and in826your PATH, but no optional Sage packages need be installed.827828EXAMPLES::829830sage: octave_console() # not tested831GNU Octave, version 2.1.73 (i386-apple-darwin8.5.3).832Copyright (C) 2006 John W. Eaton.833...834octave:1> 2+3835ans = 5836octave:2> [ctl-d]837838Pressing ctrl-d exits the octave console and returns you to Sage.839octave, like Sage, remembers its history from one session to840another.841"""842os.system('octave')843844845def octave_version():846"""847Return the version of Octave installed.848849EXAMPLES::850851sage: octave_version() # optional - octave; and output is random852'2.9.12'853"""854return str(octave('version')).strip()855856857