niszetの日記

10年目エンジニアが最近勉強したことを忘れないようにメモをする

coord_polar()+geom_path()で直線を引くには…

注:そのままでは出来ません…

昨日の記事、 niszet.hatenablog.com

について、kohske‏さんから下記のようなアドバイスをいただきました。

ありがたや…

リンク先の記事、

stackoverflow.com

のコードをコピペして動かしてみてください。記事の中で完結しているのでコピペして動かす方が手っ取り早いのです。

coord_polar()改め、coord_radar()を使用してgeom_line()を使ったときにちゃんと直線になることがわかると思います。

昨日の記事に対してやってみる(ただしこちらはデータ順に線を引きたいのでgeom_pathを使用。詳細は前回の記事にて)と、

library(tidyverse)

d=data.frame(height=c(1,2,2,3,4), weight=c(1,3,4,4,2))

dd <- d %>%  mutate(r=sqrt((height^2.0)+(weight^2.0)),theta=atan(height/weight))

p2 = ggplot() +
    geom_path(data=dd, mapping=aes(x=theta, y=r)) +
    geom_point(data=dd, mapping=aes(x=theta, y=r), size=8, fill="white", shape=21) +
    geom_text(data=dd,mapping=aes(x=theta, y=r, label=seq(1,nrow(dd))))

coord_radar <- function (theta = "x", start = 0, direction = 1) {
    theta <- match.arg(theta, c("x", "y"))
    r <- if (theta == "x") "y" else "x"
    ggproto("CordRadar", CoordPolar, theta = theta, r = r, start = start, 
            direction = sign(direction),
            is_linear = function(coord) TRUE)
}

p2 + coord_radar(theta="x",start = -pi/2, direction = -1) +ylim(0,6)+scale_x_continuous(limits=c(0,2*pi),breaks=seq(0,2*pi,by=pi/2))

で、めでたく

f:id:niszet:20170713230612p:plain

を得ることが出来ました。めでたしめでたし。

Coordの中身をちょっと見てみる。

・・・で、終わってしまうのはなんなので、ちょっとcoordの中身を見ていきます。

ggplot2のオブジェクト、関数の中身を見ていく方法については、先日も紹介したkazutanさんのこちらの資料が良いと思います。

ggplotのオブジェクトから眺めてみる

基本はstr()で見ていき、$で要素にアクセスしていきます。
また、関数は()をつけなければ(そしてRで書かれていれば)中身のソースを読むことが出来るので、それらを駆使してみていきます。場合によって、

hoge <- edit(hoge)

のようにして関数を修正して挙動を調べることもできます。今回はやってませんけど…

やってみよう。

さて、上のcoord_radarとcoord_polarの違いは何かといえば、

is_linear = function(coord) TRUE

の有り無しだけです。実際にcoord_polarの中身を見てみると

coord_polar

# function (theta = "x", start = 0, direction = 1) 
# {
#     theta <- match.arg(theta, c("x", "y"))
#     r <- if (theta == "x") 
#         "y"
#     else "x"
#     ggproto(NULL, CoordPolar, theta = theta, r = r, start = start, 
#         direction = sign(direction))
# }
# <environment: namespace:ggplot2>

同じですね。なお、namespaceはとりあえず今は考えなくて良いです。
このis_linearによって、直線になったり、対象の座標軸にそってうねったりするようです。実際にこれを使っている部分はまだ見れていません。課題。

なお、このis_linearはfunctionとして定義しています。実際にやればわかりますがTRUEを代入ではダメです。使う場所でis_linearを関数として呼び出しているためです。
引数のcoordは抜いても動きましたが、副作用がないかはわかっていません。入れておくべきでしょうね。

さて、ggprotoでCoordPolarを渡しているので、次はこれを見ていきます。
こんな感じ。

str(CoordPolar)

# Classes 'CoordPolar', 'Coord', 'ggproto' <ggproto object: Class CoordPolar, Coord>
#     aspect: function
#     distance: function
#     is_linear: function
#     labels: function
#     range: function
#     render_axis_h: function
#     render_axis_v: function
#     render_bg: function
#     render_fg: function
#     train: function
#     transform: function
#     super:  <ggproto object: Class Coord> 

うーむ、よくわからない。わからないが、それぞれの要素に関数が入っていて、親がCoordなのはわかりますね。
このis_linearにアクセスすると、

CoordPolar$is_linear

# <ggproto method>
#   <Wrapper function>
#     function (...) 
# f(...)
# 
#   <Inner function (f)>
#     function () 
# FALSE

で、FALSEが返ってくることがわかります。 後日ちゃんとソース読んで確認しますが、RStudioを使っているとCoordPolar$ まで入力すると補完候補が出ると思います。このときに一覧に出てこないものはこのオブジェクトでは定義されておらず、super(親)で定義されているものが継承されているのだと思っています。

なお、親の親の親の…と先祖をたどっていって、どこかでTRUEが定義されていて、そのあとFALSEで上書きされていなければ、TRUEを返すはず(たぶん)です。
ちょっとこの辺りもちゃんと調べます。課題。

さて、Coord*** 系の、RStudioで入力補完の候補に出てきたものについて、is_linearがT/Fどちらを返すのかを見てみました。

Coord
# Coord系のすべての元になるオブジェクト。FALSEを返す。これがデフォルトの挙動

CoordCartesian$is_linear 
# 唯一、is_linearがTRUEを返す

CoordFixed 
# superがCoordCartesianなので、is_linearも継承。TRUEを返す

CoordFlip
# Coordをそのまま継承。FALSEを返す

CoordMap
# Coordをそのまま継承。FALSEを返す

CoordPolar
# Coordをそのまま継承。FALSEを返す。前回の記事で使用したのはコレ。

CoordQuickmap
# superがCoordCartesianなので、is_linearも継承。TRUEを返す。CoordMapと名前が似ているが中身は違う点に注意が必要かも?

CoordTrans
# Coordをそのまま継承。FALSEを返す

という結果でした。MapとQuickMapで違う、とか、自分でTRUEを返しているのはCartesianだけで、あとはコレを継承しているからっぽい、とかちょっと面白い。 継承~は前記の通り、僕の認識違いの可能性もあるので、ご注意。

まとめ

is_linearについて簡単に調べるだけでもここまでボリュームが…。
…なので、オレオレgeomを作ってみたい!という方は必要になった部分から少しずつ見ていくのが良いと思いました。出来るだけggplot2と真正面から戦わない方が…
今回、サムネイルはちゃんとなっているはず…u_riboさんにアドバイスいただきました。ありがとうございました。

そして、今回はいつもにもまして勢いで書いております。ここ変だよ!とかご指摘大歓迎です。よろしくお願いします。

使っていると段々とわかってきて、わかってくると面白くなってきますね。

Enjoy!!

coord_polarを使ってみる(追記)(修正)

自作geomのその前に

complex型を使う際には何かとcoord_polarの出番が多そう。ということで、ちょっと試してみました。
今回はこちらを参考にしています。

ggplot2 Quick Reference: coord | Software and Programmer Efficiency Research Group

パイチャートというのでしょうか、あの帯状の図は目にするのですが、coord_polarを使用してgeom_lineやgeom_pointを使用している例はあまり見てない気がします。
あたしこのパイ(チャート)嫌いなのよね…

さて、下記コードの前半は上記のサイトからのコピペですが、最後のプロット時にxlim,ylimを与えて、原点が含まれるようにしています。

library(ggplot2)

d=data.frame(height=c(1,2,2,3,4), weight=c(1,3,4,4,2))

p = ggplot() +
    geom_line(data=d, mapping=aes(x=height, y=weight)) +
    geom_point(data=d, mapping=aes(x=height, y=weight), size=8, fill="white", shape=21) +
    geom_text(data=d,mapping=aes(x=height, y=weight, label=seq(1,nrow(d))))

p+coord_equal() +xlim(0,5)+ylim(0,5)

結果はこう。
f:id:niszet:20170712234040p:plain

xlimで(ylim)はなくて、下のようにscale_x**、scale_y**でも良いです。上記のオブジェクトpをstr()など使って見てみるとわかりますが、デフォルトだとlimitsが入ってないのですね。
原点を含みたいのであれば必ずそれが入るようにxlim,ylimを指定してあげればよいということになりますね(max, min の中に0が入っているか?のチェックが必要かな)

さて、coord_polarを使ってみます。

dd <- d %>%  mutate(r=sqrt((height^2.0)+(weight^2.0)),theta=atan(height/weight))

p2 = ggplot() +
    geom_line(data=dd, mapping=aes(x=theta, y=r)) +
    geom_point(data=dd, mapping=aes(x=theta, y=r), size=8, fill="white", shape=21) +
    geom_text(data=dd,mapping=aes(x=theta, y=r, label=seq(1,nrow(dd))))

p2 + coord_polar(theta="x",start = -pi/2, direction = -1) +ylim(0,6)+scale_x_continuous(limits=c(0,2*pi),breaks=seq(0,2*pi,by=pi/2))

f:id:niszet:20170712234747p:plain

geom_lineがなんか変ですが、点については多分意図したものが出来たはず…。最初の図と点の位置はあっている(と思う)ので。
geom_lineはthetaの小さい順につないでいるように見えますね。これは要改良…。

とりあえず思いついたまま書いたのであまり良くないコードになってる気がします。後日清書しよ…

次はこれをcomplex型に使ってみます。今日は時間切れ。日付変わっちゃったけど。

追記(修正)

寝てる間にgeom_lineじゃなくてgeom_pathなら良いのでは?と気づきました。

p2 = ggplot() +
    geom_path(data=dd, mapping=aes(x=theta, y=r)) +
    geom_point(data=dd, mapping=aes(x=theta, y=r), size=8, fill="white", shape=21) +
    geom_text(data=dd,mapping=aes(x=theta, y=r, label=seq(1,nrow(dd))))

p2 + coord_polar(theta="x",start = -pi/2, direction = -1) +ylim(0,6)+scale_x_continuous(limits=c(0,2*pi),breaks=seq(0,2*pi,by=pi/2))

結果はこうです。意図通りっぽい

f:id:niszet:20170713075912p:plain

以下、マニュアルより
geom_path() connects the observations in the order in which they appear in the data.
geom_line() connects them in order of the variable on the x axis.

こちらで全文が見れます。 Connect observations — geom_path • ggplot2

やはり知っているつもりでも、細かいところはこまめにヘルプを見なくてはいけませんね。

(R小ネタ) 前回の続き

小ネタです…

ggplot2のgeom_segmentは、aesとしてx,y,xend,yendをそれぞれ与えることで、各データを(x,y) (xend, yend)を結ぶ線分の形でプロットしてくれるようです。

詳しくはこちらをご参照ください
Line segments and curves — geom_segment • ggplot2

R4DSで使われていたflightsというデータセットがうまく使えそうかな?と思ったのでテストしてみます。

# install.packages("nycflights13")

library(nycflights13)
library(tidyverse)

flights %>%  na.omit() %>% ggplot()+geom_segment(aes(x=arr_time,y=dest,xend=arr_time+air_time,yend=dest))

これで、このような図が作れます。

f:id:niszet:20170711231405p:plain

前回と同じような図が出来ますね。ちょっと棒が長い気がするけど。(ひょっとしたら重なってるかもですね)
もうちょっと良いデータセットを探そうかな…xとxendにdateの方を使えばよい気もしますが…(日付型はちょっと…)

ということで(?)データからMIDIを作る方を今後は作ってみます。関数名はwrite_smf()になりますね。

次回進捗は3か月後かな…。writeはちょっと難しい(smfの仕様の把握が)ので時間がかかりそうです。
tuneRも、midiのreadはあるのですが、writeはないのですよね。なければ作ればいい…。作らなくてもいいけど…

(R進捗)rsound/csound関係進捗(追記)

あまり進捗ないですけど…

csoundと、rsoundと、rcppと…と調べることが広範囲にわたるのであまり進捗がないです。 とりあえず、進んだ分はメモ。

windows版だけかもしれないですが…

先日、rsoundの方でwaveでのファイル出力が出来ることに気づきました。その際、生成したファイルのヘッダにはちゃんとファイルサイズに関する情報が含まれていました。

Csound側で同じコードを実行してもファイルサイズは含まれていなかったで、何かの設定(の抜け)を疑いましたが見つからず。rsoundはcsound APIを呼んでいる形だから、Qtが悪いのでは…と思っています。
RStudio、nvimに続き、CsoundもQtで悩まされるのか…

rsoundで変数を渡す方法

rsoundの方、パラメータをどうやって渡すのかを悩んでいました。csound()は中でas.characterを読んでいるので、その中身を取り出して色々見ていましたが、たとえば周波数440Hzを入力する際に、パラメータとして渡したければ、下記のように`kkk=“=440"とやると良いです。変数名=値の形で与えると、csoundに渡すときに代入文にならないので、=を付けてあげればよいのでした。

orc <- create_orchestra( sr    = 44100,
                         ksmps = 32,
                         nchnls = 2,
                         `0dbfs` = 1,
                         instrument = list(
                             create_instrument(kkk="=440",aout = "vco2 0.5, kkk", outs = "aout, aout"))
)

sco <- create_score(sections = list(
    list(score_i(1, 0, 1))
))

csound_impl(as.character(orc),as.character(sco),paste0("-o","test.wav"))

同様に、これも動きます。

orc <- create_orchestra( sr    = 44100,
                         ksmps = 32,
                         nchnls = 2,
                         `0dbfs` = 1,
                         instrument = list(
                             create_instrument(aout = "vco2 0.5, p4", outs = "aout, aout"))
)

sco <- create_score(sections = list(
    list(score_i(1, 0, 1, 440))
))

csound_impl(as.character(orc),as.character(sco),paste0("-o","test.wav"))

しかしまだ、f系の音が出せていません。うーむ…まだ仕様がわからないのでした。

今後に期待

とりあえず、コードいじらなくてもまだ色々できそうなので、ちょっと試していこうと思っています。csound側の仕様を理解しないといけませんが…。
しかし…日曜にR触ってないで、お料理教室に行った方が…

追記

下記のように記述することで、foscilのサンプルも動きました。力尽きたので、csound()を経由するのはまた後日…?

foscil

eeee <- "0dbfs = 1\nksmps = 32\nnchnls = 2\nsr = 44100\n\ninstr 1\nkcps = 440\nkcar = 1\nkmod = p4\nkndx line 0, p3, 20\naout foscil .5, kcps, kcar, kmod, kndx, 1\nouts aout, aout\nendin"

ssss <- "f 1 0 16384 10 1\n\ni 1 0  9 .01\ni 1 10 .  1\ni 1 20 . 1.414\ni 1 30 5 2.05\ne"

csound_impl(eeee,ssss,paste0("-o","test.wav"))

(R進捗)midiを読んでplotするところまではできた…

残件沢山です

とりあえず、midiを読み、全データを一旦保持しつつ、noteのon/offだけを抜き出してdata.frameとして、ggplot2+geom_segmentで図にするところまでたどり着きました。

結果:

f:id:niszet:20170709122508p:plain

dfはそれぞれの音の高さ、なり始める時間、なり終わる時間などが入っています。geom_segmentはx,y,xend,yendを指定することで、その間に線を引いてくれるようです。欲しかったgeomそのものでした。

ざっくり、下記のようにしてプロットしています。

df %>% ggplot()+geom_segment(aes((x=start_time,y=typexend=end_time,yend=type))

これからまだまだデータ構造含めて修正してしまうので、細かい使い方はまだ書けないです。

(R雑記)パッケージ作成関連をいろいろ試しています。。。

進捗ダメです

ちょっとお仕事の方が忙しく、Rを触っている余裕がなかったのです。
後ほど書きますが、本業はR関係ないので、余暇に頑張らないとR力が上がらないのです。。。

rmusicwork

とりあえず最終的には音楽作りたいのですが、それ以前にR力が足りなさすぎるので色々とお勉強。
今日はroxygen2パッケージによる、パッケージのマニュアル作成関係を色々と触っていました。 github.com

R3.4にしたので色々入ってない

R3.3から3.4にするときに色々とパッケージが消えてしまいました。まぁ必要になったら入れなおせばよいのですが。

www.oreilly.co.jp

を参考に、RStudionのツールバーのBuild-Configure Build Tools… からroxygenの設定をしようとしたのですが、そもそもそんな項目がない。
その後、roxygen2が入ってないことがわかりました。入ってない場合、

devtools::document()

の時点で、入れるようにメッセージが出てきたかと思います。素直に入れましょう。そうすると、先の設定の部分にroxygen関係の設定が現れます。Build & Reload時にドキュメント生成もすることにして、作業抜けがなくなりました。
ま、まだ中身をしっかりしないといけないのですがね…!

今回はこちらも参考にさせていただきました~

qiita.com

devtools: Rパッケージ作成支援 - Heavy Watal

英語でも良ければ、Hadleyのこちらを。オライリー本は大体がweb上に公開されていますよね。でも本で読んだ方が私は好きなのですが。

Welcome · R packages

export

以前は外部に見せたくない関数名の先頭には.を付けていました。今回、exportを付けることで明示的に外部に見せる/見せないを管理できるようになったし、マニュアルぽいものも作れるようになってきて、少しずつR力の高まりを感じます。。。 …が、Rに限らないけど言語習得は時間がかかりますね…むずぃ…

はやくRおじさんになりたい(後ろ半分は達成)

Rでbode線図を描いている事例を見つけました

Rで回路設計について調べるのは、検索性が悪いと思うのですよねー…

電気回路ではRという文字を抵抗を表すシンボルに使うため、うまいこと検索でひっかけられないのでした。

最近はR lang という形で"lang"と入れるようにして検索してます。そうしたら見つけてしまった次の記事:

LC LOW-PASS FILTER ANALYSIS IN R paulorenato.com

これ、Rでやろうと思っていた小ネタの一つだったのですが、ggplot2使っている、complex型使っている、で完全に被ってました。むむむ。 ちょっとだけ解説というか、自分なりの(無駄)アレンジを入れてみます。

ちょっとだけ解説

Rには便利な型が色々とありますね。しかし、仲間はずれがいます。rawとcomplex型です。Rについては結構色々本を読んでいると思うものの、この二つについての解説はほとんど見たことがない気がします。曰く「あまり使わない」…

まぁ、普通の使い方をしていて使うシーンはあまりないとは思うのですが。 自分はraw型はWaveのファイルを読むために既に使っていますが、complex型はまだ使ってなかったです。

rawはバイナリファイルを読み書きするには必須ですが、complex型は複素数。統計分野で虚数なんて出てくるのでしょうか?使われなさそうですね。

電気電子回路や制御の分野だとフィルタの設計などでよく使うのです。最初にあげたサイトもLCフィルタについての記事となっています。
さて、ちょっと使ってみましょうか…。

まずは普通にplotします

plot(exp(-1i*seq(0,1,0.01)*2*pi))

plotはdata frameにしなくても表示できるので、ちゃちゃっと確認する際によく使いますね。
結果は下図のとおり。

f:id:niszet:20170629203826p:plain:w300

これは横軸に実数部、縦軸に虚数部が描かれていますね。 これを明示的に行うと、

data.frame(Re(exp(-1i*seq(0,1,0.01)*2*pi)),Im(exp(-1i*seq(0,1,0.01)*2*pi))) %>% plot()

で、同じ図が得られる。

つづいて、ggplotで同じことをすると、

library(tidyverse)


data.frame(re_x=Re(exp(-1i*seq(0,1,0.01)*2*pi)),im_x=Im(exp(-1i*seq(0,1,0.01)*2*pi))) %>% ggplot(aes(x=re_x,y=im_x))+geom_point()

で表示できます。

f:id:niszet:20170629203846p:plain:w300

さて、ここまでは実部と虚部をReとImを使用してそれぞれプロットしていたので、complex型をそのままy軸にプロットしているわけではない点に注意です。横軸が実数、縦軸が虚数という形。
y軸側をcomplexにしたらどうなるのでしょう?plotで描くなら下記のようになりますか

data_frame(x=seq(0,1,0.01),y=exp(-1i*seq(0,1,0.01)*2*pi)) %>% plot()

一応、下のように図が出ますが、なんか変ですね。 f:id:niszet:20170629203835p:plain:w300

このプロットの実行時には下記のWarningが出ます。

Warning message:
In xy.coords(x, y, xlabel, ylabel, log) :
  imaginary parts discarded in coercion

虚数部無視したよってことでしょうか。exp(ix)=cos(x+)isin(x)ってことで、余弦波となっているのですね。ふむ。 ggplot2の場合、同じようにすると

data.frame(x=seq(0,1,0.01),y=exp(-1i*seq(0,1,0.01)*2*pi)) %>%  ggplot(aes(x=x,y=y))+geom_point()

これはエラーが出てプロットされません。エラーは下記のような内容。複素数のままだと大小関係わからないから確かにプロットできませんよね。エラーで止まるggplot2の方が親切かもしれません。

# Don't know how to automatically pick scale for object of type complex. Defaulting to continuous.
# Error: Discrete value supplied to continuous scale

複素数は実部と虚部の2つをどうにかして表す必要があるので、最初にプロットしたように分けてプロットすることが考えられますが、もう一つ、絶対値と角度の形でも表現できます。極座標形式でしたっけ?それを使用したものがボード線図なのです。続く。

ようやくbode線図について。

ボード線図、あるいはボーデ線図と呼ばれる、周波数に対してのゲインと位相の関係をプロットした図については説明を省略…。
周波数応答系はわかりやすいというのは難しい気がするなぁ。一応参考になるかもなサイトを。 www.kairo-nyumon.com

さて、最初に貼ったサイトからソースを転載しますが、

fseries <- function(min, max, perdec) {
  expmin = floor(log10(min))
  expmax = floor(log10(max))
  step = 1/perdec
  return(10^(seq(expmin,expmax,by= step)))
}
Hf <- function(f, R, L ,C) {
  s = as.complex(2i*pi*f)
  Zser = as.complex( R + s*L)
  Zpar = as.complex(1/(s*C))
  return(Zpar/(Zpar+Zser))
}

ここまでは同じものを。先の関数は対数で見たときに等間隔になるようにする関数、後の方は伝達関数そのものを書いていますね(たぶん)

ここからはちょっと手を加えています。 ゲイン側は同じものを使用していますが、それに加えて位相分をArgによって取り出し、単位をradからdegへと変換しています。

library(ggplot2)
library(scales)
f = fseries(100,1e6,100) #100Hz to 1MHz, 100pts/decade
R = 0.06; L = 1E-6; C = 47E-6
amplitude = 20*log10(Mod(Hf(f,R,L,C))) #amplitude in dB = 2olog(|H(f)|)
df <- data.frame(x = f, y = amplitude)
p1 <- ggplot(df,aes(x,y)) + geom_line(lwd = 1, col = 'blue') +
  scale_y_continuous(breaks=seq(-65,25,by = 5)) + 
  labs(x = "f [Hz]", y = "|H(f)| [dB]") + 
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),labels =     trans_format("log10", math_format(10^.x)))+
  annotation_logticks(sides = "bt")

phase = Arg(Hf(f,R,L,C))/pi*180 #phase in deg = arg(H(f))
df <- data.frame(x = f, y = phase)
p2 <- ggplot(df,aes(x,y)) + geom_line(lwd = 1, col = 'blue') +
    scale_y_continuous(breaks=seq(-180,180,by = 30)) + 
    labs(x = "f [Hz]", y = "Arg(H(f)) [deg]") + 
    scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),labels =     trans_format("log10", math_format(10^.x)))+
    annotation_logticks(sides = "bt")

grid.arrange(p1,p2)

これで、(私)見慣れたbode線図の形となりました。

f:id:niszet:20170629203854p:plain

最後の、プロットを並べて表示についてはこちらを参考にさせていただきました。
id.fnshr.info

tidyverseからマスクされてますよってwarningが出ましたが、今回は使わないので良しとしました。

# library("gridExtra")
# 
#  次のパッケージを付け加えます: ‘gridExtra’ 
# 
#  以下のオブジェクトは ‘package:dplyr’ からマスクされています: 
# 
#      combine 

これ、geom化しようかと考えてはいたのですが、上記の通りある程度処理を加えればなんとかなってしまうので、わざわざgeomを作る必要はないのかな、と思いました。
ggplot2の中身の理解のために、これを題材に作るのは良いかもしれません。課題です。

終わりに

ということで、誰が得するのかわかりませんが、Rでも複素数を扱えて、bode線図も書くことが出来ますよ、という紹介でした。
これをRの得意な統計の分野とつなげられれば、使い道があるんですが…

しかし、リンク先にもあるように、

Well, this started as one of those “I'll do it just because I can" ...

です。Rは意外と色々なことができるので、色々試してみてください。

というわけで…

Enjoy!