自然は意外と単純なのか？
このワークシートはMath by Codeの一部です。
アート、背景、実装、バリエーションの順に見ていきましょう。
前回は花火のような粒子の動きを
運動法則のようなコードで
作ることができた。
今回は、生物の成長や動きのように、
一見複雑なアートが、
単純なルールに偶然がからむだけで疑似的に作られることを体験しよう。
今回のテーマは、
単純なルール ＋ 偶然
この発想や現象は、
さまざまなところで使われている。
ゲーム、音楽、賭け事、
いろいろあるね。
背景にあるのは、
単純なオブジェクト生成の「ルール」と
偶然を生み出す「乱数」だ。
１つずつ見てみよう。
１．背景
ランダムウォーク、千鳥足というのがある。
酔っ払いがまっすぐ前に歩けず、ふらふらと右へ左へ向きが定まらずに進むようすをさす。
このルールは単純で、進行方向は右前か左前かのどちらか。
どっちに行くかは予測できない。
ここに乱数が使える。
＜乱数の数学＞
乱数は数列だ。
一次合同式を作ってみる。
初項（乱数の種）をr0、
rn=(a rn-1+b) mod M
としよう。
数列が整数としたら、法Mによる剰余は0からM-1のM通りだから、
いつかは剰余が同じになる。
だから、デタラメな出現に見えるようにするには、Mを巨大な数にすればよいね。
この発想のように、一次合同式を使った疑似乱数を一様乱数uniform randomという。
なぜ、一様というかの理由は簡単で、０からM-1までの数が計算で出現するので、
基本的には、ゾーンに区切れば、出現の確率はどれも同じくらいになるはずだから。
質問：初項seedから一様乱数を生成するコードはどう作ればよいでしょうか。
一次合同式をf(x)=ax+b(mod M) として、初項ｘ=seedとするとき、
geogebraなら、理論上はList=Iteration( f, seed, M)とします。
現実的には、f(x)=ax + bとおき、g(x) =f- floor( f /M) Mとします。
なぜか、Mod(f, M)とすると文法エラーがおきるバグがあるみたいです。
l1=iteration(g, seed, M)とすると、初項seedの数列がM+1番目までのリストになります。
そこで、l3=take(l1, 2, M+1)とすることで、seedのあとのM個の数列ができますね。
l2=unique(l3)とすると、重複のないリストになるね。
このl2のサイズがMであるなら、M個の乱数が重複なく出そろうことになるね。
たとえば、a=2, b=3, M=7, seed=3とすると、集合｛０，２，３}で、かたよる。
たとえば、a=5, b=3, M=8, seed=3とすると、集合｛0,1,2,3,4,5,6,7}で、出そろう。
「aを８で割って５余る数、
bを奇数、
Mを2のべき乗の数
にすると、一様にM個の剰余が出そろう。」
と言われています。
下のアプレットで、数値を変えて実験してみよう。
一様乱数（実験用）
＜乱数と言語＞
乱数は偶然性を生み出すのにとても便利だ。
まず、おもな用途は
A. 基本は０と１の間の小数を生み出す。
B. それを拡大して、０とaの間の小数を生み出す。
C. aからb未満の小数を一様に出す。
D. さいころをふるように、整数のaからbまでが出る。b-a+1通りの整数を出す。
実際の言語では、どうなっているのだろうか？
・processingでは、pythonでもp5.jsでも同じ書き方。モジュールのインポートはなし。
A. ０以上１未満の小数を出す。random()
B. 0以上a未満の小数を出す。 random(a)
B. aからｂ未満の小数をだす。random(a,b)
C. aからｂまでの整数をだす。int(random(a,b+1))
・ふつうのpythonではimport randomをして
A. ０以上１未満の小数を出す。random.random()
B. ０以上a未満の小数を出す。random.random(0, a)
C. aからｂ未満の小数をだす。random.uniform(a, b)
D. aからｂまでの整数をだす。random.randint(a, b)
使う関数が一種類ならば、from random import randintのようにかけば、random.はつけなくてよい。
・ふつうのjavascriptでは、モジュールなしだが、Math.をつける。
A. ０以上１未満の小数を出す。Math.random()
B. ０以上a未満の小数を出す。Math.random()*a
C. aからｂ未満の小数をだす。Math.random()*(b-a)+a
D. aからｂまでの整数をだす。Math.floor(Math.random()*(b-a+1))+a
それぞれの言語によって特徴があります。
processing関係の乱数はrandom関数1本で済みます。便利すぎですが、一般的ではありません。
ふつうのpythonは利用者も多いので、この関数を覚えるのがいいと思います。
一方でふつうのjavascriptも利用者が多いのですが、いちいち関数を加工する必要があります。
まあ、B,Cタイプなら、その場対応でなんとかなるでしょう。
しかし、Dタイプを使うときは、
アロー関数などで、次のように定義しておくと便利でしょう。
let randint = (a,b) => Math.floor(Math.random()*(b-a+1))+a;
こうするとrandint(1,6)と書くだけで、さいころの目関数として使えますね。
＜正規乱数＞
乱数はゲームなどでは一様乱数ですむことが多いかもしれません。
でも、現実の世界では、一様よりも、中心にかたまる、正規分布の方が自然だ。
だからシミュレーションなどでは、正規分布する乱数が使われることが多い。
質問：0以上１未満を返す一様乱数random()から平均m,標準偏差sdの正規乱数を作るには、どうしたらよいでしょうか。
一様乱数はｘが０以上１未満のとき、確率密度関数f(x)=1/(1-0)=1となるね。
連続変数の分散と平均の出し方は、離散変数のΣを∫にしただけだから、Σで思い出して、積分に直そう。
平均はE(x)=Σｘp(x)=m,分散、つまり、偏差2乗の平均 V(x)=1/nΣ(x-m)2=.......=E(x2)-m2
だから、
ｘの平均はE(X)=
ｘ2乗の平均はE(X2)=
だから、ｘの分散は
ということは、
「中心極限定理」から、どんな分布でも、平均値の確率分布は正規分布に収束するので、
一様乱数の平均値を振り続けると、正規分布に収束するでしょう。
一様乱数randomを12回ふった値の合計resは、
平均が、分散が に近づき、その標準偏差がsd＝１の平方根＝１となる。
見かけは積分記号で一見すごそうにみえるけど、
ただのΣ計算、ただのたし算を連続変数化しただけで、もともと線形な処理だ。
だから、resに線形な処理をして、分散を変形できるはずだ。
たとえば、resから6を引いて平均を０にする。
さらに、標準偏差sd、平均mを指定した正規分布になるようにするには、
res*sd+mをすると、正規乱数として返してくれるね。
たとえば、pythonでやってみよう。
#[IN]Python
#====================
# norm from uniform
import matplotlib.pyplot as plt
import pandas as pd
from random import random
counter = 1000
def norm(ave, sd ):
res = 0
ans=0 ;
for i in range(12):
ans = random()
res +=ans
res -=6
return res *sd + ave
def docalc():
ave = 50
sd = 10
data=[]
for i in range(counter):
datai= norm(ave, sd)
data.append(datai)
#print(data)
df = pd.Series(data)
df.hist(bins=8, range=[20, 80])
plt.show()
docalc()
のような処理を、言語の文法に合わせて実行すればよいでしょう。
一様分布関数から、正規分布に近い出力ができるね。
[OUT]
＜ルール＞
絵の部品になる図形に単純なルールを決めよう。
（ルール０）
・1個の点を画面の下（南）におく。
・点は1回で北東か北西に進む。どちらにするかはランダム。
・進んだ跡を画面にかく。
このルールでかくと、ランダムウォークに見えるアートができるでしょう。
（ルール１）
・円の1個目は中央におく。
・2個目からの円は平面にデタラメにばらまく。
・追加の円は、一番近い円に移動してから円を画面にかく。
このルールでかくと、冒頭の「植物の枝の成長」のように見えるアートができます。
（ルール２）
・円の1個目は中央におく。
・2個目からの円は平面にデタラメにばらまく。
・追加の円は、他の円に重ならない大きさの円を画面にかく。
このルールでかくと、円の隙間がどんどん埋められていくアートができるでしょう。
（ルール３）
・円は８方位に毎回ランダムに移動できる。
・移動したら、移動したあとを画面にかく。
このルールでかくと、ランダムウォークの8方位バージョンのアートができるでしょう。
あとは、図形や線にランダムまたは、位置に応じた色をつけることで、
さらにアートに変化がでるでしょう。
２．実装
ルール１を実装することで、
冒頭の図をかいてみよう。
＜pyhon in processing＞
'use strict';
maxCount = 5000 # max count of the cirlces
currentCount = 1
x = []
y = []
r = []
def setup():
size(800, 800)
strokeWeight(0.5)
# first circle
x.append( width/ 2)
y.append( height/2)
r.append(10)
def draw():
global currentCount, x,y,r
background(255)
newR = random(1, 7)
newX = random(newR, width - newR)
newY = random(newR, height - newR)
closestDist = 99999999999999999999 #max of distance
closestIndex = 0
for i in range(currentCount):
newDist = dist(newX, newY, x[i], y[i])
if newDist < closestDist:
closestDist = newDist
closestIndex = i
angle = atan2(newY - y[closestIndex],newX - x[closestIndex])
x.append( x[closestIndex] + cos(angle)*(r[closestIndex] + newR))
y.append( y[closestIndex] + sin(angle)*(r[closestIndex] + newR))
r.append( newR )
currentCount +=1
for i in range(currentCount):
fill(255*x[i]/width,255*y[i]/height,255*r[i]/7) #rule for coloring
ellipse(x[i], y[i], r[i]*2, r[i]*2) #draw circle
if currentCount >= maxCount:
noLoop()
単純で行数も少なく冒頭のキレイなアートができる。
しかし、意味が分かりにくいので、javascriptでクラスを使ってかき、
geogebraで動くようにしよう。
前回と同様に、javascriptで円のクラスを定義する。クラス利用（生成と描画）を分離してかく。
（円の定義）
コンストラクター（生成関数）で決めよう。
半径Rは１から７、X座標はー(画面横幅-円の半径)から(画面横幅-円の半径)とするこで、
円が画面から切れないようにする。Y座標も同様にしよう。
円の色、R,G,Bは255を最大に対して、200,100,50にして、茶色っぽくする。
円の連番は０番にかりに入れる。
・setId関数で円の連番をあとで設定できるようにする。
・setRXY関数で円のサイズR,座標X,Yをあとで変更できるようにする。
・neaness(other)関数で、他の円otherとの距離の2乗を求める。一番近い円を探すときに使おう。
・approach(other)関数で、他の円otherとの中心を結ぶ線にそって、中心どうしの距離が、
半径の和（this.R+other.R)になるまでくっつくときに、自分のX,Y座標を計算して更新する。つまり、くっつける。このスライドをするための角度angleを中心を結ぶ線を斜辺とする直角三角形のｘ方向の長さとｙ方向の長さから、そのatanから求めます。Math.atan2関数が使えます。
下の図のように、
this.Xは、半径の和のcos(angle)だけother.Xたし、
this.Yは、半径和のsin(angle)だけother.Yにたします。
・draw(app)関数で、1個の円をかきます。
画面のappletのオブジェクトを受け取り、アプレット画面に対する描画関係のAPI関数を使う。
円を1個かき、fillを１にして中もそめる。ラベルなし、軸なし、格子線なしにする。
（円の利用）
リセットボタンを押す(OnClicked)と、setup関数を実行します。
画面には隠れてますがスラーダーtがアニメーションで動き、最新情報(OnUpdate)でdraw関数を実行します。このdrawは円の１こずつのdrawではなく、画面全体のdrawにかかわります。
・setup関数で画面サイズをwidth,heigtともに600になるようにしつつ、原点が中央にくるようにします。
1個目の円を描く前に、書いてある円をぜんぶ削除します。1個目の円は半径10で原点を中心にかく。
・draw関数で、追加の円をかきます。追加円と既存円との近さを求めて、一番近い円の番号が
closestIIndexです。この番号の円に追加した円をapproachします。
それから、追加円のdrawをしましょう。
・円の集まりリストcirclesに円のインスタンスを追加したり、円の番号をセットするsetIDなどの
タイミングがずれると、ちがう絵になったりするので、このあたりの検証は
前回紹介した、VSコードでgeogebra のwebへの埋め込みなどでデバックするといいね。
冗長にはなるけれども、手続き処理をする書き方よりも、改造がしやすい利点があるね。
approach法
＜geogebra by javascript ＞
let windowggbApplet = ggbApplet;
'use strict';
let circles = []; // 円たち
const maxCount = 5000; // max count of the cirlces
let currentCount = 1;
let width = 600;
let height = 600;
let randint = (a,b) => Math.floor(Math.random()*(b-a+1))+a;
function ggbOnInit() {
setup();
}
class Circle {
constructor() {
this.R = randint(1,7);
this.X = randint(this.R - width , width - this.R);
this.Y = randint(this.R - height, height - this.R);
this.colR= 200;
this.colG= 100;
this.colB= 50;
this.num = 0;
}
setId(num){
this.num =num;
}
setRXY(r,x,y) {
this.R = r;
this.X = x;
this.Y = y;
}
nearness(other){
return Math.pow((this.X- other.X),2) + Math.pow((this.Y- other.Y),2);
}
aproach(other){
let angle = Math.atan2(this.Y - other.Y, this.X - other.X);
this.X= other.X + Math.cos(angle) * (other.R + this.R);
this.Y= other.Y + Math.sin(angle) * (other.R + this.R);
this.colR= this.X *255/width+122;
this.colG= this.Y *255/height+122;
this.colB= this.R *255/7+122;
}
draw(app){
let posinfo = `c${this.num} = Circle(Point({${this.X},${this.Y}}),${this.R})`;
app.evalCommand(posinfo);
app.setColor("c"+this.num, this.colR, this.colG, this.colB);
app.setFilling("c"+this.num,1);//fillingout
app.setLabelVisible("c"+this.num,false);//NO label
app.setAxesVisible(false,false);//No Axes
app.setGridVisible(false);//No Grid
}
}
function setup() {
// 既存の点cオブジェクトあれば削除
if (circles.length > 0) {
for (let i = 0; i < circles.length; i++) {
windowggbApplet.deleteObject("c" + i);
}
}
currentCount = 1;
circles = [ ];
windowggbApplet.setCoordSystem(- width/2, width/2, - height/2, height/2);
//1個目の円をかく。
let firstC = new Circle();
firstC.setId(0);
firstC.setRXY(10,0,0);
circles[0] = firstC;
firstC.draw(windowggbApplet);
}
function draw() {
//円の追加
let closestDist = Number.MAX_VALUE;
let closestIndex = 0;
let newC = new Circle();
// 追加円に、もっとも近い円を探す。
for (let i = 0; i < currentCount; i++) {
let target = circles[i];
let newDist = newC.nearness(target);
if (newDist < closestDist) {
closestDist = newDist;
closestIndex = i;
}
}
//追加円をかく。
newC.setId(currentCount);
newC.aproach(circles[closestIndex]);
circles[currentCount] = newC;
newC.draw(windowggbApplet);
currentCount++;
}
３．バリエーション
質問：ルール２のように、円が重ならないように増えていくようにするにはどのようなコードにしたらよいでしょうか。
円をばらまくところはルール１と同じです。
ただ、円に重ならないようにするしくみを作る必要がありますね。
まず、円のコンストラクタで半径を50にしますが、最低半径は２にします。
そうすると、追加したい円thisと既存の円otherとの隙間がないと追加したものを無効にしましょう。
centerDist(other){ return Math.sqrt(Math.pow((this.X- other.X),2) + Math.pow((this.Y- other.Y),2));}で中心距離を求めます。
circumDist(other){return this.centerDist(other) - other.R -2;}周囲と周囲との最大のすきまを求めます。
noRest(other){return this.circumDist(other)<0 ? true: false;}は隙間がないときにtrueです。
もし、既存の円全部を調べて、追加するすきまがあれば、
一番隙間の狭い円に対して、標準半径50と最低半径２の間でくっつけらたらくっつける調整をします。
fit(other){
let maxR = this.centerDist(other) - other.R;
let res= (minR<=maxR)? true: false;
if(res){
if (maxR < this.R) {
this.R = maxR;
}
}
return res;
}
可能な最大半径をmaxRとして、minR=2以下ならfitは成功とします。さらに、maxRが50を切っているなら
半径を50からmaxRに更新して、隙間にぴったり詰め込めばよいですね。
理論上はこれを、drawで実行すると、毎回円を追加するたびに、隙間に入ります。
実際、VSコードとｊｓでwebにgeogebraアプレットを埋め込みができます。
しかし、これを素のgeogebraアプレットのスクリプト記述に貼り付けても
必ず動く保証はありません。
ここでは、上のアプレットをコピーして、グローバルjavaスクリプトに貼り付けただけですが、
最初は動きませんでした。しかし、スライダーの最新情報のjavascriptにalert("draw ");を入れたら動きました。
意味がわかりませんね。
同一シートにjavascriptのアプレットが２つ以上あると内部動作が不安定になるとか、裏で動いているgeogebraの謎の機構によるものと思われます。
参考までに、下に貼り付けたコードをのせておきますね。
重ならない円が増殖する
let Applet = ggbApplet;
'use strict';
let circles = []; //
const maxCount = 5000; // max count of the cirlces
let currentCount = 1;
let width = 600;
let height = 600;
let randint = (a,b) => Math.floor(Math.random()*(b-a+1))+a;
let minR = 3;
function ggbOnInit() {
setup();
}
class Circle {
constructor() {
this.R = 50;
this.X = randint(this.R - width , width - this.R);
this.Y = randint(this.R - height, height - this.R);
this.colR= 200;
this.colG= 100;
this.colB= 50;
this.num = 0;
}
setId(num){
this.num =num;
}
setRXY(r,x,y) {
this.R = r;
this.X = x;
this.Y = y;
}
centerDist(other){
//中心間の距離
return Math.sqrt(Math.pow((this.X- other.X),2) + Math.pow((this.Y- other.Y),2));
}
circumDist(other){
//周囲と周囲との最大のすきま
return this.centerDist(other) - other.R -2;
}
noRest(other){
//重なるときtrue
return this.circumDist(other)<0 ? true: false;
}
fit(other){
//一番周囲が近い他円に対して、半径の最適化が成功するとtrue
let maxR = this.centerDist(other) - other.R;
//console.log("this.R,maxR=",this.R,maxR);
let res= (minR<=maxR)? true: false;
if(res){
if (maxR < this.R) {
this.R = maxR;
}
console.log("this.R=",this.R);
}
this.colR= this.X *255/width+122;
this.colG= this.Y *255/height+122;
this.colB= this.R *255/7+122;
return res;
}
draw(app){
let posinfo = `c${this.num} = Circle(Point({${this.X},${this.Y}}),${this.R})`;
app.evalCommand(posinfo);
app.setColor("c"+this.num, this.colR, this.colG, this.colB);
app.setFilling("c"+this.num, 1);//filling
app.setLabelVisible("c"+this.num,false);//NO label
app.setAxesVisible(false,false);//No Axes
app.setGridVisible(false);//No Grid
}
}
function setup() {
//
if (circles.length > 0) {
for (let i = 0; i < circles.length; i++) {
Applet.deleteObject("c" + i);
}
}
currentCount = 1;
circles = []
Applet.setCoordSystem(- width/2, width/2, - height/2, height/2);
//1個目の円をかく。
let firstC = new Circle();
firstC.setId(0);
firstC.setRXY(50,0,0);
circles[0] = firstC;
firstC.draw(Applet);
}
function draw() {
//円の追加
let closestRest = Number.MAX_VALUE;
let closestIndex = 0;
let newC = new Circle();
let IntersectOthers = false;
// 他円との重なりチェック
for (let i = 0; i < currentCount; i++) {
let target = circles[i];
IntersectOthers ||= newC.noRest(target);
}
//console.log("intersect is ",IntersectOthers)
if (IntersectOthers==false){
// 追加円に、もっとも近い円を探す。
for (let i = 0; i < currentCount; i++) {
let target = circles[i];
let newRest = newC.circumDist(target);
if (newRest < closestRest && newRest >= 0) {
closestRest = newRest;
closestIndex = i;
}
}
//半径が範囲内であれば、追加円をかく。
if (newC.fit(circles[closestIndex])){
newC.setId(currentCount);
//console.log("(this , near) =",newC.num, closestIndex);
circles[currentCount] = newC;
newC.draw(Applet);
currentCount++;
}
}
}