Rlab 简介:用于科学和工程应用的高级语言

作者:Ian Searle

当我意识到自己花费了太多时间编写 Fortran 和 C 语言程序来进行工程分析、测试数据缩减和继续教育时,我开始使用高级语言。我寻找更好的工程编程工具;除了 Matlab 和一些类 Matlab 程序(均为商业软件)外,我一无所获(那是在 1989 年)。我认为 Matlab 的语言设计不够强大,但如果它的成本更低并且可以在更多平台上使用,我会使用它。我将“业余时间”的研究转向了解释器构建,并开始原型化 Rlab。大约一年后,我将 0.10 版本发布到互联网上征求意见。现在,将近五年过去了,Rlab 经历了重大的变化和改进,这主要归功于世界各地用户的出色帮助。

Rlab 并不试图成为 Matlab 的克隆版本。相反,它借鉴了我认为 Matlab 语言的最佳特性,并提供了改进的语言语法和语义。语法已经得到改进,以允许用户更自由地表达,并减少歧义。变量作用域规则已得到改进,以方便创建更大的程序和程序库。添加了一个异构关联数组,以允许用户创建和操作任意数据结构。

概述

术语“高级语言”有许多不同的含义。在 Rlab 的上下文中,它意味着变量没有显式类型或维度。变量的类型和维度是从用法中推断出来的,并且可以在程序执行期间任意更改。这种自动类型化和动态分配使科学程序员从许多通常与 C 或 Fortran 等语言相关的耗时编程任务中解放出来。

在 Rlab 的情况下,高级也意味着交互式。由于目的之一是为数据操作和可视化提供便捷的方法,因此该程序可以在批处理或交互模式下运行。在 shell 提示符下键入 rlab 即可启动程序。Rlab 为用户提供命令行提示符 >,程序语句由此发出并执行。

您能得到什么?它花费多少?大多数人发现,他们开发一个可用的程序比使用 C 或 Fortran 快得多。此外,内置的图形功能和交互式操作模式促进了数据和方法的实验和探索。您付出的代价是在 某些 情况下程序执行速度较慢。编译语言和解释语言之间的性能测试通常表明,编译语言在所有方面都更快。如果您的 Rlab 程序专门使用标量运算,并且不使用任何内置的线性代数函数,那么它将比 Fortran 或 C 版本慢得多。但是,如果您利用 Rlab 优化的数组运算、内置的线性代数和实用函数,您可能会发现 Rlab 的性能非常好。

虽然我无法在一篇文章中提供完整的介绍,但我可以向您展示足够的语言内容,以帮助您确定它是否是您想要尝试的东西。“Rlab 入门”中提供了更完整的介绍,该入门包含在源代码和二进制发行版中。

数据

数值数组或矩阵是最常用的数据类型。数值数组可以是一维的(向量)或二维的。大多数情况下,矩阵是通过从另一个程序或文件中读取数据来创建的。对于更简单的示例,我们将以交互方式创建矩阵。数组/矩阵语法使用方括号 [ ](如 C 语言)来分隔矩阵,并使用分号将每一行与下一行分隔开来。因此,要在 Rlab 中创建矩阵,用户需要键入

> a = [1, 2; 3, 4]

这将提供以下响应

  a =
        1          2
        3          4

矩阵引用和赋值使用以下语法

> a[1;] = 2 * a[2;2,1]
 a =
        8          6
        3          4

分号分隔行和列规范。缺少行或列规范将选择所有行或列。此外,用户可以按任何顺序指定特定的索引或索引组合。在前面的示例中,a[2;2,1] 选择 a 的第二行,以及该行的第二列和第一列(按该顺序)。因此,a[2;2,1] 的计算结果为 [4,3]。然后将此量乘以 2,并将结果分配给 a 的第一行。

在本文中,所有示例都将使用实数。但是,Rlab 也可以很好地处理复数。要输入复数,可以使用虚数常数 ij(不用担心,您仍然可以使用 ij 作为变量)。例如:a = 1 + 2j 创建实部为 1,虚部为 2 的复数。所有数值运算符和函数(在有意义的情况下)都适用于复数值以及实数值。

传统的数值运算符被重载。当处理标量操作数时,结果是大多数人所期望的。但是,当处理数组/矩阵操作数时,数值运算符的行为有所不同。加法和减法函数在两个操作数之间以逐元素的方式运行。乘法运算符 * 执行矩阵内积。为了说明

> a = [ 1 , 2 , 3 ];
> b = [ 4 , 5 , 6 ];
> a + b
        5          7          9
> a' * b
        4          5          6
        8         10         12
       12         15         18
> a * b'
       32

请注意(在前两行中),语句末尾的 ; 会抑制打印结果。除法运算符 / 对标量操作数执行传统除法,并使用矩阵/向量操作数求解联立线性方程组。存在一组逐元素运算符,用于在数组之间执行逐元素运算。逐元素运算符是双符号运算符;第一个符号始终是 .,因此逐元素除法运算符是 ./,逐元素乘法运算符是 .*

> a = [1,2,3;4,5,6];
> b = 1./a
 b =
        1        0.5      0.333
     0.25        0.2      0.167
> a.*b
        1          1          1
        1          1          1

除了数值数据类之外,还有一个字符串类。字符串数组/矩阵的处理方式与数值数组完全相同

> s = ["how", "to create" ;
>      "a string", "within rlab"]
 s =
how          to create
a string     within rlab

如您所见,字符串数组由任意长度的字符串组成。对数组中字符串的相对大小没有限制。字符串可以使用 + 运算符连接。其他(传统的数值)运算符不适用于字符串数组。

列表

C 语言的优点之一是能够创建数据结构以适应特定的编程任务。Rlab 提供了一个类似的编程结构,称为列表。在其最简单的形式中,列表是单维关联数组。但是,由于列表元素可以是字符串、数字、数组、函数和其他列表,因此它们提供了一个强大的多维工具。列表的创建方式与数值数组类似

> l = << 3*pi ; rand(3,4) ;
>        type = "my-list" >>
 l =
   1            2            type

<< >> 分隔符包含列表元素,; 分隔元素。= 用于为元素分配索引名称(不能是数字),如上面列表的第三个元素所示。如果未为列表元素分配索引,则默认情况下会分配一个数值,如前两个元素所示。

列表的元素通过其索引引用:l.["type"] 返回字符串 "my-list"l.[1] 返回 3*pi 的值。对于显式索引名称,可以使用简写符号:l.type 也将返回字符串 "my-list"

列表索引可以替换为计算结果为字符串或数值标量的表达式,从而允许用户以自动方式访问选定的元素。如果 index 是一个包含字符串 "type" 的变量,则表达式 l.[index] 将访问 type 元素并返回字符串 "my-list"。如果 index 包含数字 1,则 l.[index] 将返回 9.42。列表将在本文后面再次讨论。

开始使用

Rlab 为用户提供用于条件执行的 conditional 和 if 语句,用于循环功能的 for 和 while 语句,以及函数或子例程。当我们继续进行一些示例时,将介绍这些功能。

Rlab 函数有点不寻常,值得关注。函数虽然不是一等公民,但却是对象,用变量引用。因此,要创建和保存函数,必须将其分配给变量。函数参数传递和变量作用域规则旨在方便创建更大的程序和程序库。我们的第一个示例,积分常微分方程将说明其中的一些功能。

一个流行的例子是 Van der Pol 方程,因为它简单、非线性,并且演示了极限环振荡。Rlab 有一个内置的积分引擎和几个“rfile”积分器。rfile 是一个包含 Rlab 程序或语句的文件。积分函数需要以下输入:返回导数值的函数、定义积分区间的开始和结束的值以及初始条件。

整个问题可以交互式执行。但是,我选择将程序放在一个文件 (int.r) 中。这允许我编辑和重新执行,而无需进行大量重复的键入。对于本文,我使用了 shell 转义功能(第一列中的 \)通过 pr 对文件进行 cat 操作,以生成 列表 1

由于 ode 函数积分一阶微分方程,我们必须将 Van der Pol 方程写成

An Introduction to Rlab: A High Level Language for Scientific and Engineering Applications

作为两个一阶方程

An Introduction to Rlab: A High Level Language for Scientific and Engineering Applications

计算并返回

An Introduction to Rlab: A High Level Language for Scientific and Engineering Applications

的函数被编写并分配给变量 vdpol。定义函数后,变量 t0(积分开始时间);tf(积分结束时间);和 x0(初始条件)被初始化(第 9-11 行)。内置计时函数 tictoc 在调用 ode 时使用,以告诉我们积分需要多长时间。

完成此文件后,以下命令将执行它,就像从提示符中输入每一行一样。

"> rfile int
ODE time:      0.620

这个简单的问题运行速度相当快。ode 的输出,即与时间相关的值的矩阵,

An Introduction to Rlab: A High Level Language for Scientific and Engineering Applications

存储在变量 out 中(第 14 行)。在研究微分方程的行为时,便捷的数据可视化是一个真正的优势。在本例中,我们想要查看

An Introduction to Rlab: A High Level Language for Scientific and Engineering Applications

(x) 和

An Introduction to Rlab: A High Level Language for Scientific and Engineering Applications

(xd) 随时间的变化。我们还想查看此问题的相平面,即

An Introduction to Rlab: A High Level Language for Scientific and Engineering Applications

与 的关系图。

An Introduction to Rlab: A High Level Language for Scientific and Engineering Applications

我们可以使用 plot 函数来完成此操作,该函数绘制矩阵列。如果 plot 的输入是单列矩阵,则矩阵行索引将针对横坐标轴绘制,矩阵列元素将针对纵坐标轴绘制。如果 plot 的输入是 N 列矩阵,则第一列中的值将针对横坐标轴绘制,其余列中的值将分别针对横坐标轴绘制。

可以使用命令 plot(out); 创建类似于 图 1 的图。由于 out 的第一列是时间,第二列和第三列是 xxd,我们所要做的就是将未更改的矩阵提供给 plot。如果我们想要绘制相平面,如 图 2 所示,我们需要指定我们想要绘制第三列与第二列的关系图。为此,只需从 out 中提取第二列和第三列,如下所示:plot(out[;2,3]);

Plplot 图形库为 Rlab 提供内置的 2D 和 3D 图形功能。提供了用于大多数常见图形功能的内置函数,例如:2D、3D、网格图、等高线图和直方图。Plplot 支持大多数常见的图形输出设备,最值得注意的是 X-Windows 和 PostScript。

Plplot 图形可能不足以满足需求。在这些情况下,将正确格式化的数据传输到另一个程序的能力至关重要。有几种将 Rlab 与其他程序连接的方法。与大多数优秀的 Unix 应用程序一样,Rlab 将从 stdin 读取并写入 stdout。有一些函数用于以 ASCII 格式写入矩阵,以及一个类似 C 的 [f]printf 函数。还有一个 system 函数,允许 Rlab 程序执行可以使用 Unix shell 完成的任何操作。但是,使与其他程序连接最容易的是写入和读取进程(管道)的功能。

输入/输出

I/O 系统的用户端旨在尽可能简化,同时又不限制需要它的用户的能力。文件和进程通过字符串标识。每个能够读取或写入的函数都会在必要时打开和关闭文件。但是,为特殊情况提供了 openclose 函数。有几个版本的读取和写入函数可以理解 Rlab 的所有内部数据结构,并提供一定程度的 Matlab 文件兼容性。用于特殊 ASCII 输入需求的 getlinestrsplt 函数、用于二进制输入的 fread 以及用于特殊格式化输出的 fprintf 完善了 I/O 函数。

正如您可能期望的那样,字符串 "stdin""stdout""stderr" 指向其 Unix 系统对应项。任何执行文件 I/O 的函数也可以通过管道执行进程 I/O,只需将文件名字符串替换为进程名称字符串即可。进程名称字符串是以 | 作为第一个字符的字符串。字符串的其余部分是任何 shell 命令。Rlab 将创建一个管道,fork 并调用 shell。管道的方向是从用法中推断出来的。此功能使与其他程序之间的数据传输相当简单。Rlab-Gnuplot 接口就是一个很好的例子,它完全用 Rlab 语言编写,使用进程 I/O 功能将数据和命令传输到 Gnuplot。

作为演示,我们将通过一个简单的 X-Geomview 程序接口来探索进程 I/O。X-Geomview 是一个功能强大的 3 维图形引擎,带有交互式 GUI。X-Geomview 可以从文件读取数据,也可以从 stdin 读取数据/命令。X-Geomview 使用 Motif,但静态链接的 Linux 二进制文件(以及源代码)可从 www.geom.umn.edu/software/geomview/docs/geomview.html 获取。

在此示例中,我将生成数据并绘制经典的墨西哥帽。代码列在 列表 2 中。

示例的数据在第 14 行完成;从那里开始,我们只是将数据发送到 X-Geomview 进程。变量 GM 保存一个字符串,其第一个字符是 |,表示应该打开一个进程到字符串的其余部分。以下语句(第 16-21 行)将对象定义发送到 X-Geomview,第 23-30 行包括嵌套的 for 循环,这些循环将多边形顶点坐标发送到 X-Geomview。包含结果的 X-Geomview 窗口的快照显示在 图 3 中。当然,创建此类图的更好方法是创建一个自动化 X-Geomview 接口的函数(这将包含在 Rlab 的下一个版本中)。

操作工作区

高级语言非常适合原型化过程,但在实际“生产”情况下使用程序时,通常会略有不足。在此示例中,我们将假装我们刚刚开发了一种计算时频分布的方法。实际上,我们将使用 Rene van der Heiden 的 tfd 函数,该函数源自 Choi 和 Williams 的论文“使用指数核改进多组分信号的时频表示”。

现在我们要使用 tfd 处理大量数据。由于 tfd 执行速度相当快,我们不希望仅仅为了能够处理大量数据而不得不使用其他语言重新编写它。假设您有很多时间历史数据文件,您希望“推送”通过 tfd。有些文件包含单个事件数据矩阵,而另一些文件包含多个数据矩阵。您希望能够编写一个程序,该程序可以以最少的用户干预处理所有此类文件。对于某些语言来说,困难在于无法操作变量名和隔离数据。

Rlab 通过列表解决此问题。列表允许创建和隔离任意数据结构,并提供一种系统地操作工作区变量和变量名的机制。我之前展示了如何使用字符串访问列表元素。列表也可以使用字符串变量或任何计算结果为字符串的表达式进行索引。

我尚未透露的有趣之处在于,整个工作区都可以被视为列表!通过特殊符号 $$ 授予对工作区的访问权限。您可以将 $$ 用作列表变量的名称。例如,您可以像这样使用余弦函数:$$.["cos"](pi),或:$$.cos(pi)。第一种方法提供了最大的灵活性。现在我们了解了这一特殊功能,我们可以相对轻松地处理我们的问题。该程序将一次读取每个包含数据的文件(它们与模式 *.dat 匹配),计算时频分布,并将每个数据集的分布保存在工作区中。处理完成后,新的工作区将保存在单个文件中以供以后使用。

刚刚描述的程序包含在 列表 3 中。我应该指出几件事

  • 第 1 行:require 语句检查工作区中是否有名为的函数。如果找到该函数,则什么也不会发生。如果未找到该函数,则从磁盘加载它。

  • 第 10 行:getline 函数读取 ASCII 文本文件,并自动将每一行解析为字段(有点像 awk)。字段(字符串或数字)在列表中返回。当 getline 返回一个空列表时(由 length 函数检测到),while 循环终止。

  • 第 12 行:每个文件名都存储在字符串数组中。

  • 第 13 行:readb 函数从每个文件中读取所有数据,并将其存储在列表 $$.[filenm[i]] 中。这是工作区中的一个变量,其名称与文件名相同。例如,如果第一个文件是“x1.dat”,则将创建一个名为“x1.dat”的列表变量。

  • 第 24 行:现在我们将对我们读取的数据进行操作。该程序将循环遍历数组 filenm 中的字符串。

  • 第 26 行:对于每个文件 (i),程序将循环遍历每个文件中包含的所有数据。members 函数返回列表元素名称的字符串数组。

  • 第 28 行:就是这样!$$.[i] 是工作区中的一个列表(每个数据文件一个)。$$.[i].[j] 是第 i 个列表中的第 j 个变量。因此,我们正在计算我们读取的每个文件中的每个矩阵的时频分布。$$.[i].[j+"y"] 为每个执行的时频分布在每个列表中创建一个新变量(与矩阵名称相同,但在末尾附加了“y”)。

总结

还有许多我没有讨论过的功能、函数和用户贡献的程序。特别值得注意的是大量的线性代数函数集。Rlab 为 Netlib 和 Statlib 提供的 LAPACK、FFTPACK 和 RANLIB 库提供了非常方便的接口。

到现在为止,您应该已经看到了足够的语言内容来决定是否值得尝试一下。如果您想了解更多信息,可以查看:www.eskimo.com/~ians/rlab.html。Rlab 的 Linux 二进制文件现在以 RPM 格式分发,并且可以在 ftp://ftp.eskimo.com/u/i/ians/rlab/linux 获得。现在 ELF 发行版已可用,并且已脱离 beta 测试,并且考虑到 ELF 的动态链接有多好,Rlab 二进制文件是使用 ELF 对象文件格式构建的。

致谢

Plplot 图形库由德克萨斯大学的 Maurice LeBrun 博士提供。底层线性代数子例程(LAPACK 和 BLAS)来自 Netlib 存储库。当然,如果没有 GNU 工具,这一切都不可能实现。多年来,许多个人为 Rlab 做出了其他贡献。源代码发行版中的 ACKNOWLEDGMENT 文件试图提及所有人。

在哪里获取 Rlab

Rlab 随附一百多个各种类型的 rfile,其中许多由用户贡献。纽约州波茨坦市克拉克森大学的 J. Layton 教授正在完成 Rlab 控制工具箱的最后工作。曼彻斯特大学/UMIST 的 Higham 教授的测试矩阵工具箱的端口也已可用。有关更多信息,请访问 www.eskimo.com/~ians/rlab.html

Ian Searle 目前在华盛顿州西雅图市的航空航天研究领域工作,并在业余时间从事 Rlab 的工作。

加载 Disqus 评论