経験的Mathematica入門

松田裕幸 matsuda@symbolics.jp (2011.5.6)

code: https://colab.research.google.com/drive/1bEEOaJm5rhlLyAs97Rm5cdT6tXdo-q5z

Mathematicaはプログラミング言語

Mathematicaは超高級電卓だという話を聞いていた。しかし、面白そうなので
『Mathematicaクックブック』を眺めるとなんか雰囲気が違う。
Mathematicaっていったい何者？

Mathematicaはプログラミング言語の１つです。JavaやC++などと比べるとあり得ない機能をたくさん持ち、表記も一風変わっています。しかし、プログラミング言語です。

プログラミング言語というからにはそれを使ってプログラムを書きます。では、Mathematicaで書くプログラムとはどんなものでしょう？

はじめに確認しておかなければならないことは、Mathematicaではプログラムはすべて式の集まりであり、Mathematicaにおいては式は関数と等価です。通常、プログラムで式というと、こんな感じのものを思い浮かべるかもしれません。

$3 - \sqrt{4}$

Mathematicaでも同じように書けますが、この式を関数を使って書き直すと次のようになります。

Subtract[3, Sqrt[4]]

これを分解してみましょう。

Subtract   関数名Subtract(引き算)
[          普通関数というと( )が一般的だがMathematicaでは[  ]を使用
 3,        関数Subtractの第1引数
 Sqrt      関数Subtractの第2引数、同時に関数Sqrt(平方根)
  [
  4        関数Sqrtの引数
  ]
 ]
 注：関数の中にまた関数がある場合、内側から計算していきます。

一般的にプログラムは「実行する」と言われます。しかし、Mathematicaのプログラムは式あるいは関数のため、実行するとは言わず、評価する(evaluate)といいます。ではさっそくこの式を評価したいのですが、その前にMathematicaにおけるプログラミング環境について確認しておく必要があります。拡張子 *.nbを持つこのファイルはノートブックと呼ばれます。そして画面右端にみえる ] マークで囲まれた単位をセルと呼びます。（注：本資料はもともとMathematicaノートブックで作られたものをPythonノートブック*.ipynb上に再現しているためセルに階層は存在しません。）

これで準備が整いました。Mathematicaの式があるセルをクリックし、Shift+Enterキーを押してください（Shitfキーを押したままEnterキーを押します）。

In [ ]:
Subtract[3, Sqrt[4]]
Out[ ]:
1

式の左側に現れたラベルが入力と出力を表します。In[]が評価前の式を表し、Out[]が評価結果を表します。当然、最初に示した式も同じように評価できます。（注： Mathematica ノートブックと異なり、Jupyter ノートブックでは $\sqrt{4}$ を直接入力することはできません。）

In [ ]:
3 - Sqrt[4]
Out[ ]:
1

これをプログラムと言われても違和感を感じるかもしれません。どうみても電卓と同じじゃないか、どこがプログラミング言語なのか？何万行ものMathematicaプログラムなんてあるのか？当然の疑問です。でも答えはYESです。

通常は自分で書いた関数（方法は後述）が数行から数十行。そうした関数を使って、また別の関数を作り、などしていくと、あっというまに数百行に達します。MathematicaはOSや他の言語とのインタフェースも充実していますので、Mathematicaでシステム寄りのプログラムを書こうとすると自然、規模が大きくなっていきます。

Mathematicaの組み込み関数はすべて大文字から始まっているのですぐにわかります。逆に、ユーザが書く関数はこれと区別するために通常小文字で始めるのが基本です。また、すべての関数名はそれ自身が関数の意味を表すよう長い名前が推奨されています。

MathematicaはC++などやJavaなどと異なりインタプリタ言語です。したがって、プログラミング第1歩に登場する"Hello,World"を作る時もmain関数などの必要ありません。

In [ ]:
Print["Hello,World"]

Hello,World

記号計算

ちょっと他の言語では表現できないものの1つが記号計算です。たとえばC言語で次の計算は可能でしょうか？変数に値が入ってないと計算できず、エラーとなるはずです。

In [ ]:
b + 3 +2 a
Out[ ]:
Output

Mathematicaでは変数（正確にはシンボル）は変数のまま式の中に書けます。その一方、そのシンボルに対する値が与えられたときは式全体を評価し値を持つことができます。"/."は置換の意味。"a->9"は/.左辺にあるaすべてを9で置き換えます（後述）。

In [ ]:
b + 3 +2 a /. {a -> 4, b->2} (* 2+3+2*4 *)
Out[ ]:
13

Mathematicaはこの能力を使うことで初等数学の問題を解くことができます。（関数の説明は省略します。大まかにいうと、２つの方程式を使い、グラフを描き、方程式を解くことで交点を求め、今度はその交点をグラフ上に加えている、が分かれば十分です。）

In [ ]:
exp1 = x^2 + y^2 == 1; 
exp2 = y - 2 x^2 + 3/2 == 0;
ContourPlot[
    Evaluate@{exp1, exp2}, {x, -1.5, 1.5}, {y, -1.5, 1.5}, 
 ImageSize -> 200]
Out[ ]:
Output
In [ ]:
pts = Solve[exp1 && exp2, {x, y}](*方程式を解く *)
Out[ ]:
Output

解いた答えを座標として赤点で描きます。円(exp1)と放物線(exp2)が交差する場所が連立方程式の解になっていることがわかります。

In [ ]:
Show[{
  ContourPlot[
       Evaluate@{exp1, exp2}, {x, -1.5, 1.5}, {y, -1.5, 1.5}],
  Graphics[{Red, PointSize[Medium], Point[{x, y} /. pts]}]}, 
  ImageSize -> 200]
Out[ ]:
Output

オブジェクトとリスト

Mathematicaにとってすべての計算対象はオブジェクト呼ばれます。オブジェクトはオブジェクト指向のオブジェクトではなく、関数引数の要素、あるいは関数の計算結果になるという意味でオブジェクトです。たとえば、上記の例では、計算し可視化された画像自体がオブジェクトして変数 img に束縛されています。

そして以降、この img に束縛された画象がオブジェクトとして操作されるのを見ていきます。

In [1]:
img = ExampleData[{"TestImage", "House"}]
Out[1]:
Output

フロントエンドを持たないWolfram Languageであっても、計算し求めた画像を Export 関数によって外部ファイルに保存することができます。

In [2]:
Export[ "house.svg", img]
Out[2]:
house.svg

Mathematicaではデータの基本はオブジェクトですが、複数のオブジェクトを構成するにはリストを使います。もちろん、リストもオブジェクトですので、リストのリストも当然ありえます。リストは集合と異なり、異なるデータタイプが混在できます。(ここで、Graphics、Circleは関数です。Red,Thick, Dashedは関数Circleに対する属性、ImageSizeはグラフィックスの大きさを指定するオプションです。）

In [3]:
List[
  3.14,
  img,
 Graphics[{Red, Thick, Dashed, Circle[]}, ImageSize -> 100]
 ]
Out[3]:
Output

あるいは{ }を用い、次のように記述します。ここに限らず、Mathematicaは同じことを表現するのに複数の表現を持ちます。どれを使うかは使う時の状況、使用者の好みに依存します。

In [4]:
{3.14,  img,  ImageMultiply[img, 3.14]}
Out[4]:
Output
In [5]:
ImageDimensions[img] (* 画像の大きさを縦横のピクセル数で与える *)
Out[5]:
{256, 256}

Mathematicaではオブジェクトを関数によって直接操作し、結果も直接見ることができます。

In [6]:
ImageTake[img, 50]
Out[6]:
Output
In [7]:
ImageRotate[img, Right]
Out[7]:
Output

リストデータを作るにはすでに見たように要素を直接指定する以外に、計算によって作り出すこともできます。次の例は、関数Tableを使い、画像を15°ずつずらしたもののリストを作り出しています。Table内でImageRotate関数を複数回実行し、それぞれの結果からなる画像のリストを用意します。ImageRotateの中のTableパラメータは$\theta$ で、これを０からPiまでPi/6刻みに増やしています。一方、ImageRotateは元の画像 img を$\theta$分回転しますが、その際、回転後の画像の大きさを{100,100}に指定し、かつ、標準で残る黒い背景をTransparentで消しています。

In [8]:
m = Table[
  ImageRotate[img, theta, {100, 100}, 
    Background -> Transparent], {theta, 0.,  Pi, Pi/6}]
Out[8]:
Output

実行可能なハイパーリンク構造を持ったヘルプ

すでに書きましたようにMathematicaには4,000近い関数があります。当然、そのすべてを覚えることはできないし、かつ、覚えていても細かい使い方までは覚えていません。そこでノートブックから興味ある関数を見つける方法をお教えします。? でパターンにマッチする関数名を探します。次の例では関数名に Rotate を含むものをすべて列挙しています。概略を知りたい関数をクリックします。するとオレンジ線の下に関数引数と簡単な説明がでてきます。そして一番最後の >>をクリックするとHelpのドキュメントに移動します。そこから関連する項目に移動することも可能です。

?\*Rotate\*

他のプログラミング言語のHelpと最も異なる点はHelp画面もノートブックで出来ているので、そこで各関数を評価し様々なテストが可能になっていることです。これは極めて凄いことです。Help画面の中で、プログラムも実行できるなんて、他の言語ではありえるでしょうか？実は本当のことをいえば、そういう言語もなくはありませんが、ノートブックという形でプログラム開発環境と透過な形で実現しているのはMathematicaだけだと思います。

注： Wolfram EngineにはHelpが提供されていません。オンライン上のヘルプ https://reference.wolfram.com/language/ を利用します。

リスト操作

この回転画象のリストmの要素数は

In [ ]:
Length[m]
Out[ ]:
7

あるいは

In [ ]:
m // Length
Out[ ]:
7

です。さらには

In [ ]:
Length@m
Out[ ]:
7

なんて芸当も可能です。これはどれがベストな表記ということではなく、プログラム中の位置、あるいはプログラミングの状況に応じて自然と使い分けられていきます。

リスト中の特定の要素を取り出すにはPart関数を使います。ちなみに、多くのプログラミング言語と違い、リストを配列を見なすとindexは0からではなく1から始めます。本来多くのプログラミング言語がindexを0から始めたのは計算効率の理由が1番でそれ以上の理由は無く、0から始めることによって得られるメリットより危険性の方が大きいので、Mathematicaの姿勢は正しいと思います。

In [ ]:
SeedRandom[1]; m=RandomInteger[100, 10] (* 0から100までの整数乱数10個からなるリストを生成 *)
Out[ ]:
{80, 14, 0, 67, 3, 65, 100, 23, 97, 68}
In [ ]:
Part[m, 4]
Out[ ]:
67

では最初の３つ、あるいは最後の２つとかを取り出すには2つの方法があります。

In [ ]:
{Take[m, 3], Take[m, -2]}
Out[ ]:
{{80, 14, 0}, {97, 68}}
In [ ]:
{m[[1 ;; 3]],  m[[-2 ;; -1]]}
Out[ ]:
{{80, 14, 0}, {97, 68}}

Mathematicaのプログラミングに慣れるには、基本的なリスト構成関数、リスト操作関数をしっかり学ぶことをお勧めします。『Mathematicaクックブック』でよく取り上げられているものとして次のものがあります。

In [ ]:
(* 1から10000までの整数乱数20個のリスト *)
SeedRandom[1];  L = RandomInteger[{1, 10000}, {20}]
Out[ ]:
{1793, 8580, 1945, 613, 8962, 3960, 3245, 7216, 3653, 3397, 1401, 9772, 4269, 2343, 3788, 3858, 6917, 5767, 5507, 4396}
In [ ]:
Partition[L, 4]
Out[ ]:
{{1793, 8580, 1945, 613}, {8962, 3960, 3245, 7216}, {3653, 3397, 1401, 9772}, {4269, 2343, 3788, 3858}, {6917, 5767, 5507, 4396}}
In [ ]:
ListConvolve[{2, 3}, {a, b, c, d, e, f}]
Out[ ]:
Output
In [ ]:
Select[L, PrimeQ] (* 素数だけを取り出す *)
Out[ ]:
{613, 6917, 5507}
In [ ]:
Cases[L, a_ /; a > 8000]  (* Lの中の要素で8000より大きな数のリスト。文法の意味は後述 *)
Out[ ]:
{8580, 8962, 9772}
In [ ]:
Count[L, a_ /; a > 8000]
Out[ ]:
3

関数型制御構造

Mathematicaは一般的なプログラミング言語が持つ制御構造をより高度に、より見やすくした表記をいくつも持っています。そしてこうした表記を使うと、Do、For, Whileなどより計算速度が高速になるというメリットもあります。

Nest, NestList

In [ ]:
a = x;
Do[a = f[a], {5}];
a
Out[ ]:
f[f[f[f[f[x]]]]]
In [ ]:
Nest[f, x, 5]
Out[ ]:
f[f[f[f[f[x]]]]]
In [ ]:
NestList[f, x, 5]
Out[ ]:
{x, f[x], f[f[x]], f[f[f[x]]], f[f[f[f[x]]]], f[f[f[f[f[x]]]]]}
In [ ]:
Nest[Sqrt, 2, 4]
Out[ ]:
Output

FixedPoint, FixedPointList

以下の例は具体的な価無しに計算させるのが難しいので、評価させていません。

a = x;
While[a != f[a], a = f[a]];
a
FixedPoint[f, x]
In [ ]:
NestList[1 + Floor[#/2] &, 1000, 20]
Out[ ]:
{1000, 501, 251, 126, 64, 33, 17, 9, 5, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2}
In [ ]:
FixedPoint[1 + Floor[#/2] &, 1000]
Out[ ]:
2

Fold, FoldList

In [ ]:
a = x;
Do[
   a = f[a, i],   {i, 1, 3}];
a
Out[ ]:
f[f[f[x, 1], 2], 3]
In [ ]:
Fold[f, x, {1, 2, 3}]
Out[ ]:
f[f[f[x, 1], 2], 3]

Map

In [ ]:
l = {};
data = {10, 20, 30};
Do[
   AppendTo[l,  f[ data[[i]] ]], {i, 1, 3}];
l
Out[ ]:
{f[10], f[20], f[30]}
In [ ]:
Map[f, data]
Out[ ]:
{f[10], f[20], f[30]}
In [ ]:
Map[f, data] /. {f -> Sqrt}
Out[ ]:
Output

ルール型プログラミング

実はMathematicaでは関数型プログラミングよりルール型プログラミングの方が難しいかもしれません。関数型はあくまで関数型の制御構造のみを理解すればいいのに、ルール型プログラミングに現れるパターンと実際の働きとの連関を掴むのが難しいからです。

しかし、すでに見た例から１つだけ取り上げ、説明してみます。

In [ ]:
L
Count[L, a_ /; a > 8000]
Out[ ]:
Out[76]: 
{1793, 8580, 1945, 613, 8962, 3960, 3245, 7216, 3653, 3397, 1401, 9772, 4269, 2343, 3788, 3858, 6917, 5767, 5507, 4396}
Out[77]: 
3
In [ ]:
count = 0;
Do[If[L[[i]] > 8000, count++], {i, Length[L]}];
count
Out[ ]:
3

_ はパターンを表します。この場合、Lの任意の要素を対象とし、a_ とすることでそれに名前（シンボル） a を付与します。そして\/\; で条件を追加します。これを一般的なプログラミング言語風に書いてみると上のようになります。

MathematicaのHelpドキュメントには膨大な解説 guide/RulesAndPatterns が掲載されています。

ここまで読まれてきたにぜひ、確認していただきたいことがあります。あえてFortran風表記とMathematica風表記を併用させてきましたが、Mathematica風表記が「恰好いい」から、そちらを使いなさい、ではないということです。もちろん、記述が簡潔な方が見やすいし、間違いも少なくなります。しかし、本当の理由は性能です。

Mathematicaのようなインタプリタ言語の場合、コンパイラ言語と異なりシンボル（変数）への参照は常にメモリ参照を引き起こし、速度に対し致命的です。

つまりMathematicaで少しでも速いプログラムを書くための鉄則の１つとして、不必要にシンボルを使わないことです。上記Fortran風コードを見てお分かりのように、多くの変数が使われています。これがコンパイラ言語ならば賢いコンパイラがこうした変数をすべてレジスタ参照に置き換えてくれるので、まったく気にする必要はないのですが、Mathematicaはインタプリタ言語です。ですので、繰り返しの中に不必要に変数を介在させることは望ましくありません。

関数定義

さて最後に自前の関数を作ってみます。

In [ ]:
Clear[f, a, b];

関数 f を定義します。

In [ ]:
f[a_, b_] := Sqrt[a^2 + b^2]

定義した関数 f に引数3, 4を与えます。

In [ ]:
f [3, 4]
Out[ ]:
5

何か変な感じしませんか？特に従来のプログラミング言語しか知らないとなおさらかと思います。上記定義を、一般的な表記に変えてみます。

$f(a,b) = \sqrt{a^{2}+b^{2}}$

そう、こんな感じですよね。では、これに合わせて、上記定義を書き換えてみましょう。

In [ ]:
Clear[f];
f[a, b] = Sqrt[a^2 + b^2]
Out[ ]:
Output

ありゃ、何か変です。でも構わず先に行ってみます。

In [ ]:
f[3, 4]
Out[ ]:
Output

やっぱり変でした（笑）。前者と後者には変更点が２箇所（正確には３箇所）あります。１つは :\= についてです。Mathematicaではシンボルの評価は「置き換え」だと書きました。関数の参照はまさに置き換えです。そのため、通常、関数を定義する際には、関数名とその引数に本体を束縛させます。したがって、最低限、次のような形に戻します。

In [ ]:
Clear[f];
f[a, b] := Sqrt[a^2 + b^2]
f[3, 4]
Out[ ]:
Output

それでも結果は変わりません。

どこがまずいのでしょう。Mathematicaは数値のみならずシンボルも扱えることを思い出してください。そして、関数引数の場所では、代入という作業ではなく、パターンマッチングという作業が起きることを覚えてください。つまり、引数にa, bとあるということは関数 f を呼び出す側にも a と b が必要だということです。

In [ ]:
f[a, b] (* symbol 'a' matches to 'a' and 'b' matches to 'b' *)
Out[ ]:
Output

確かにうまく行きました。しかし、本来は任意の数をマッチさせたかったのに。どうすればいいんでしょうか？

Mathematicaには便利な記号があります。つまり、任意のものをマッチさせる記号があります。それが _です。

In [ ]:
Clear[f];
f[_, _] := Sqrt[a^2 + b^2];
f[3, 4]  (* 3 matches to _ and 4 matched to _ *)
Out[ ]:
Output

今度は確かに引数3と4はマッチしたようです。が、答えは望むものではありません。マッチさせるだけでなく、マッチさせたものをなんらかの方法で保持させなければなりません。そのため、マッチさせた後、保持させるために _ の前にシンボル名を置きます。そして関数本体ではこのシンボル名を使い計算を行うようにすればいいわけです。

In [ ]:
Clear[f];
f[a_, b_] := Sqrt[a^2 + b^2];
f[3, 4] (* a->3, b->4  and calculates Sqrt[a^2 + b^2] *)
Out[ ]:
5

局所変数

関数本体が大きくなってくると、一時的に使いたい変数が出てきます。そうした変数を気楽に使っているとよそに悪さを及ぼします。

In [ ]:
c = 3;
f[a_, b_] := (c = Sqrt[a^2 + b^2]; c*c);
f[3, 4]
Out[ ]:
25
In [ ]:
c
Out[ ]:
5

さて何が問題でしょう？ 関数 f の外にある変数 c の値は元々3でした。それが関数ｆを使ったからといって突然5に変わってしまうのは嬉しくないはずです。なぜこんなことが起きたのでしょう？関数 f の中で一時変数 c を無造作に使ったためです。そうか、では、絶対他では使いそうにない名前を付ければいいのではないだろうか？たとえば、ccccccとか。しかし、プログラム規模が大きくなればなるほど、未使用の確認は大変になってきます。また、後の修正でたまたま同じ名前の変数を使ってしまう可能性も出てきます。つまり、特殊な名前を使うアイデアは神経を使うだけ使って結局は問題を回避出来ません。

究極の解は、一時的に使う、つまり局所的に使う変数を指定する方法を使うことです。MathematicaにはそのためにModuleという機能が用意されています。

In [ ]:
c = 3;
f[a_, b_] := Module[{c = Sqrt[a^2 + b^2]},  c*c];
f[3, 4]
Out[ ]:
25
In [ ]:
c
Out[ ]:
3

複数定義？

一般的なプログラミング言語では同じ関数名が複数の定義を持つことは許されていません。しかし、Mathematicaはそれが可能です。これもパターンマッチングに関係します。Mathematicaではすべてのシンボルは評価において置き換えの対象になります。その置き換えが可能かどうか判断するためにパターンマッチングを使います。ただし、それが関数の場合は引数部分も含めてパターンと考えます。

すると次のケースはどうなるでしょう？　関数名は同じ。そして引数パターンも同じ。ということは後から定義されたものが残り、前に定義されたものは消えます。

In [ ]:
Clear[f]
f[a_, b_] := Sqrt[a^2 + b^2]
f[a_, b_] := Module[{c = Sqrt[a^2 + b^2]}, c*c]
In [ ]:
?f
Out[ ]:
Output

次に同じ関数名だけど引数パターンが異なるものを追加したらどうなるでしょう？ 今度は引数パターンが異なるので、シンボル f には２つの定義が存在することになります。

In [ ]:
f[x_] := x^3
In [ ]:
?f
Out[ ]:
Output

ルール型プログラミングの例

In [ ]:
Clear[f]
f[0] := 1;
f[n_] := n*f[n - 1] (* 再帰関数 *)
In [ ]:
f[3]
Out[ ]:
6

次のフィボナッチ数列を計算する例も再帰関数ですが、ただし再帰の途中で一度計算済みの値を何度も再計算しています。

In [ ]:
fibonacci[0] = 0;
fibonacci[1] = 1;
fibonacci[n_] := fibonacci[n-1] + fibonacci[n-1]
In [ ]:
Table[fibonacci[i], {i, 0, 10}]
Out[ ]:
{0, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512}
In [ ]:
Timing[fibonacci[22]]
Out[ ]:
{2.6605, 2097152}

ほぼ同じ定義ですが、一度計算した値を保持することで、計算時間は n に比例した時間となります。

In [ ]:
Clear[fibonacci]
fibonacci[0] = 0;
fibonacci[1] = 1;
fibonacci[n_] := fibonacci[n]= fibonacci[n-1] + fibonacci[n-1]
In [ ]:
Timing[fibonacci[22]]
Out[ ]:
{0.001231, 2097152}

次の例では関数シンボル名も引数パターンもまったく同じです。ただし、引数に対する条件（　/;　以降に記載）が異なります。こうした条件もパターンマッチングの際、考慮されます。よってこれら３つの定義はまターンマッチングの観点から正確に区別されます。

In [ ]:
Clear[f];
f[a_] := a*2 /; a > 0
f[a_] := a/2 /; a < 0
f[a_] := 0     /; a == 0
In [ ]:
f[-9]
Out[ ]:
Output

なお、この例は次のように書くこともできます。

In [ ]:
Clear[f];
f[a_/; a>0] := a*2
f[a_/; a<0] := a/2
f[a_/; a==0] := 0
In [ ]:
f[-9]
Out[ ]:
Output

おまけ

Mathematica ノートブック上ではノートブック固有の表記　tutorial/TwoDimensionalExpressionInputOverview が使えますが、Jupyterノートブック上のWolfram Engineでは、出力表示のみにしか登場しません。入力はあくまで、キーボードでタイプ可能な文字種に限られます。

In [ ]: