問
- 名前のリストが入ったテキストファイルを読み込み、ソートする。
- アルファベットを1から26までの整数に対応付けてその和を計算し、さらに名前の順位を乗ずることで「名前のスコア」を算出する。
- 例えばCOLINなら3 + 15 + 12 + 9 + 14 = 53で、順位は938なので、53*938=49714がCOLINのスコアになる。
- 全ての名前のスコアの和を求めよ。
解答
データの読み込みと基本的な文字列操作と演算ができれば特に問題はない…。
続きを読むデータの読み込みと基本的な文字列操作と演算ができれば特に問題はない…。
続きを読む昔作った資料をもとに話す機会があったので、スペインの施設園芸について手持ちの文献を再度調べていた。
スペインの南部のアルメリア県にはハウスが集中している地域がいくつかあって、特にアルメリア市の西部に位置するアドラ市とエルエヒド市を中心とするエリアは明確にハウスが密集している。衛星写真に切り替えて確認してみてほしい。
このエリアの施設面積だけで2万haあると言われており、日本全国の施設面積の1/3程度に相当する。
この地域でこのように施設園芸が普及するまでの経緯が面白い。もともとこの地域はブドウ(ワイン用ではなく生食用だったそうだ)の栽培が盛んで、棚仕立てでブドウを栽培していたのだが、ある時生産者がこの棚の上にビニールを載せることでその下で色々な野菜を作れる、ということに気付いたそうだ。
そこから施設園芸が広まったので、ハウスの形状も独特であった。今もこの地域で最も普及しているハウスはRaspa y Amagado(ラスパ・イ・アマガド)と呼ばれるもので、ベースはビニールの載ったブドウの棚なのだが、雨水を排水するために屋根がゆるくギザギザと波打っている。Raspaは魚の骨、Amagadoは屈曲を意味し、天井部の高い部分がRaspaで谷になっている部分がAmagadoなのだそうだ。
この話を「施設と園芸」の2012年の夏号(No.158)で読んで、一度見てみたいものだ(写真は載っていたがモノクロで分かりにくかった)と思っていたのだが、昨日改めて航空写真などを確認していてそう言えばストリートビューというものがあったと気付き、適当な場所を確認してみたら記事の説明そのもののRaspa y Amagadoを確認することができた。天井部の高い部分はまさに背骨という感じだ。
上に示した場所の場合は最新の写真が2011年7月のものだが、ハウスの中は空っぽだ。この地域では夏は生産物単価が下がるのと、ハウス内部が高温になって作業が大変という理由から基本的に9月から翌年5月までの作型で栽培をするそうだ。同じ場所の2009年1月の写真を遡って確認できるので見てみると、ハウスの中にトマトが植わっているのが分かる。
夏の間は前作の片付けや次作の準備を進めつつビーチで優雅に過ごすそうだ。
ハウスはなんだか安っぽく見える(実際建造費は安いのかもしれない)が、開口部はかなり目合の細かいネットで完全に覆われているのが伺える。出入り口は基本的に二重扉となっており、外部からの害虫の侵入を防止している。中には送風装置を備えて体に付着した害虫を払い落とすような仕組みを備えた施設もあるそうだ。日本の先進地域並の装備である。
加えてスペインの農業で特徴的なのは天敵の利用が一般的であるという点だ。天敵利用は高度な技術を必要とするが、ヨーロッパでは有機農業への需要が高いこともあり、天敵利用技術が発達しており、天敵などの生物農薬を生産販売するメーカーも複数ある。販売する苗にあらかじめ天敵を定着しておくということも行われているそうだが、これは生産者の天敵利用に対する深い理解と技術力がなければ到底不可能なことだ。天敵利用に関しては明らかに日本の先を行っている。
かつては気候に恵まれ安く野菜が生産できるという点だけが利点で、オランダのような施設園芸の先進国は天敵利用技術に先んじていたため、多少値段が高くても有意性を維持できていた。というようなことが古い教科書には書いてある。
そのスペインに天敵利用が急速に広まり、HACCPやグローバルGAPを取得する生産者も続々出てきて施設園芸の先進国も苦しい立場に…というのも少し昔の話で、今ではモロッコやイスラエルが農業新興国として安くて高品質な野菜を輸出するようになっており、今度はスペインの立場も危うくなっているそうだ。
などということが書いてあった記事に関心していたのももう5年も前の話だ。5年もたてば状況も変わっている可能性が高い。また時間をみつけて新しい情報を確認したいところだ。
技術的に難しい部分とか特に無いのでこちらで。
Wikipediaで確認しよう。ちなみに英語版のほうが詳しいぞ。
あまりきちんと説明するとボロがでるのでいい加減に説明すると、暖かいものほど強い光を出して、暖かいものほど光の波長のピークが短波長側にずれるというやつだ。
こういうのよく見かけるだろう。
この手のプロットだと英語版Wikipedia(Planck’s law - Wikipedia)にあるやつが割と好きだ。太陽の表面温度くらいの放射スペクトルがちょうど可視光になってるのがよく分かる。
温度も連続変数のはずだ。だとすれば、なんかこう、グラデーションっぽい感じでプロットしてほしい。
今回はggplot2
だけです。
library(ggplot2)
k = 1.38e-23 # Boltzmann constant h = 6.62e-34 # Planck constant c = 3e8 # speed of light p_rad <- function(lambda, t){ 2*h*c^2/(lambda^5) * 1/(exp((h*c)/(lambda*k*t)) - 1) } lim = 5000 t = lim - (1:(lim^(1/1.3)))^1.3 lambda = 10^-9 * 1:2000 param = expand.grid(t=t, lambda=lambda) u = cbind(param, rad = p_rad(param$lambda, param$t))
温度と波長について少しずつ値を変えながら色々な組合せでの放射強度を計算していくわけだけれども、温度については先程のプロットからも分かるように、高温側では少し値が変わっただけで放射強度も大きく変わり、低温側ではそれほど変わらないという関係がある。
従って、高温側の隙間を小さく、低温側の隙間を大きくするような感じで調整してやると多分ポイントのムラが減る。
また、「色々な組み合わせ」はexpand.grid
でやればパッと作れる。最初for
とrbind
で適当にやったら大変時間がかかった。
jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000")) ggplot(u, aes(x = lambda, y = rad, color = t)) + geom_point(size = 0.5) + scale_color_gradientn(name = "Temperature(K)", colors = jet.colors(7)) + theme_classic() + labs(x = "Wavelength (m)", y = "Intensity (W/sr/m^2)") + theme(legend.position = "right")
jetというのは少し前までMATLABのデフォルトだったカラーマップで、色々なところでよく見かける。作り方はUsing jet colormap in R and ggplot2 – Hi!!を参考にした。
そしてgeom_point
でプロットする。とても、とても時間がかかる。
ほかに上手いやり方があるような気がして仕方がないが、思いつかなかった。
ほ場に有機質資材、例えば牛ふん堆肥であったり稲わらであったりというものを施用した場合、これらに含まれている窒素というのは必ずしもすぐには使いやすい形にならない。
有機質資材中の窒素は基本的には微生物の分解を受けることで植物が利用しやすい形態となる。有機質資材はそれぞれ特有の割合で炭素と窒素を含んでいて、これをC/N比と言う。そして、有機物を分解する微生物も炭素と窒素を体に含んでいて、微生物のC/N比はだいたい20くらいである。従って、C/N比が20より小さいときは、微生物が有機物を体内に取り込むと、窒素のほうが多いので窒素が土壌中に放出される。アバウトな説明だがこれが有機物の肥料効果が現れる流れである。従って、C/N比が低い有機物ほど肥料効果が高く、早い。
では20より大きいと炭素が放出されるかというとそうではなくて、微生物は土壌中にもともとあった窒素を利用して増殖をはじめる。従って土壌中の窒素分が吸われてしまい、負の肥料効果とでも言うべきものが現れる。C/N比が大きな有機物を投入することで土壌中の窒素がなくなってしまう現象を窒素飢餓と呼ぶ。
しかし、微生物に取り込まれた窒素も、微生物が死ねば土壌中に放出される。炭素は呼吸に伴って炭酸ガスとして減少していくので、基本的には有機物中の窒素はいつかは全て放出される(ガス化する場合もあるが)。
したがって、毎年有機物を入れ続ければ、土壌中に放出される窒素の量は年々増えていき、放出量は投入量に漸近していく。しかし、平衡状態に達するには通常長い年月がかかる。何年か連用したときの窒素放出量を知るすべはないだろうか。
土壌中の有機物の分解モデルはいくつか考えられているものがあって、有機物が半減期の異なるいくつかの種類のコンパートメントの複合したものであるという前提でモデル化したもので比較的良好な結果がえられているものが多い。
今回はhttp://agriknowledge.affrc.go.jp/RN/2010322673のモデルを使う*1。
Rでやっていこう。
まず、paramというデータテーブルが次のように用意してあるものとする。
資材名 | a | c | f |
---|---|---|---|
乾燥牛ふん | 0.10 | 0.58 | 0.32 |
牛ふん堆肥 | 0.04 | 0.15 | 0.81 |
牛ふんオガクズ堆肥 | 0.39 | 0.11 | 0.85 |
豚ぷん堆肥 | 0.02 | 0.39 | 0.59 |
豚ぷんオガクズ堆肥 | 0.02 | 0.34 | 0.64 |
バーク堆肥 | 0.14 | -0.08 | 0.96 |
未熟稲ワラ堆肥 | -0.06 | 0.26 | 0.80 |
中熟稲ワラ堆肥 | 0.03 | 0.18 | 0.79 |
完熟稲ワラ堆肥 | 0.15 | 0.11 | 0.74 |
稲ワラ | -0.40 | 0.95 | 0.45 |
小麦ワラ | -1.83 | 1.63 | 1.20 |
モミガラ | 0.19 | -0.11 | 0.92 |
水稲根 | -0.35 | 0.63 | 0.72 |
オガクズ | -0.03 | -2.77 | 3.80 |
製紙カス | -1.15 | 0.49 | 1.66 |
発酵製紙カス | 0.43 | -0.37 | 0.94 |
余剰汚泥 | 0.65 | 0.28 | 0.07 |
こんな感じで。data.tableは気分。
library(dplyr) library(data.table) library(ggplot2) library(ggrepel)
毎年の有機物投入量を1として、t年後の残存率、累積率、放出率は、
## functions --------------- residual_rate <- function(a, c, f, t){ # 残存率 a * 0.01^t + c * 0.63^t + f * 0.955^t } accum_rate <- function(a, c, f, t){ # 累積率 a * 0.01 * (1-0.01^t)/(1-0.01) + c * 0.63 * (1-0.63^t)/(1-0.63) + f * 0.955 * (1-0.955^t)/(1-0.955) } release_rate <- function(a, c, f, t){ #放出率 1 - (a * 0.01^t + c * 0.63^t + f * 0.955^t) }
そしてtをベクトルで与えられるように関数を整備して、プロットに備える。
sim_fun_gen <- function(t, name, fun){ tmp_p <- param[資材名==name] a <- tmp_p$a c <- tmp_p$c f <- tmp_p$f sapply(t, function(t) fun(a, c, f, t)) } resid_rate_sim <- function(t, name){ sim_fun_gen(t, name, residual_rate) } accum_rate_sim <- function(t, name){ sim_fun_gen(t, name, accum_rate) } rel_rate_sim <- function(t, name){ sim_fun_gen(t, name, release_rate) }
実際に一番知りたいのは窒素の放出量だ。連用何年目に何kg窒素が出てくるのかということだ。
放出量は、放出率に実際に資材が持っている窒素の量を乗ずれば出てくる。
t = 0:50 shizai <- c("豚ぷん堆肥", "豚ぷんオガクズ堆肥", "牛ふん堆肥", "牛ふんオガクズ堆肥", "稲ワラ", "中熟稲ワラ堆肥", "小麦ワラ", "オガクズ") n_rate <- c(2.7, 1.4, 1.1, 0.8, 0.6, 0.5, 0.4, 0.1) names(n_rate) <- shizai result <- data.table(資材名=character(0), 連用年数=numeric(0), 窒素放出量=numeric(0)) for(i in 1:length(n_rate)){ result <- rbind(result, data.table(資材名 = names(n_rate)[i], 連用年数 = t, 窒素放出量 = rel_rate_sim(t, names(n_rate)[i])* n_rate[i] * 10) ) }
result %>% ggplot(aes(x = 連用年数, y = 窒素放出量, color = 資材名)) + geom_line() + theme_bw(base_family = "Osaka") + geom_text_repel( data = subset(result, 連用年数==max(連用年数)), aes(label = 資材名), family = "Osaka", nudge_x = 10, size = 3, segment.size = 0.25, segment.alpha = 0.5 ) + lims(x = c(min(result$連用年数), max(result$連用年数)*1.3)) + theme(legend.position = "none") + geom_hline(aes(yintercept = 0), lty = 2, cex = 0.1) + labs(title = "有機質資材を1t/10aずつ連用した水田での無機体窒素放出経過", x = "連用年数(年)", y = "無機体窒素放出量(kg/10a)")
どの資材も平衡状態とみなせるような状態になるまでに50年以上はかかる。
そして、C/N比の低い資材は最初窒素を吸収する性質があることもわかる。オガクズなどは30年近く連用してやっと窒素が出てくる。このような資材を使う場合は窒素を補ってやらないと窒素飢餓が起きる危険性がある
一方、豚ぷん堆肥などはもともと窒素含量が多いというのもあるが、5年も連用しないうちに窒素の放出量は10kg/10aを超え、無視できない量になっている。このような資材の場合は投入量を制限するか、連用を控えるなどの対策をしなければ窒素過剰*2となる危険性があることが分かるだろう。
*1:実は元の論文はあまり確認してなくて、パラメータ等は堆肥・有機質肥料の基礎知識を参考にした
gmpで解いた。
どの数のときにどれだけの文字数がどこに追加されるかを考えて、1000個のベクトルにまとめて加算していった。
最初綴りを間違えていて解けなかった。恥ずかしい。
文字列の読み込みは直接書くよりも一旦ファイルに保存してread.table()とかでやるほうが楽なようだ。
あとは動的計画法でボトムアップ方式で解く。
月初のみチェックすれば良いので、各月の日数を加算していけばよい。
gmpで解いた。
問題自体は難しくないし計算時間もかからない。文字列の読み込みと添字操作の練習。
gmpを使って無理矢理やっても2〜3秒だが、overviewにより効率的な方法が書いてある。約数の数同士の関係をちゃんと考えてやれば1秒以内に終わる。
ただの和なのでgmpで普通に解ける。積とかなら工夫が居るかもしれない。
試した範囲では普通にメモしながらrepeat回していくのが結局一番早かった。それでも10秒くらいかかる。
ただの組み合わせ。
ゴンペルツ曲線では相対成長速度と重量の関係はどうなってるのかと思って適当に回帰してみたら
が良く当てはまったので、この仮定(相対成長速度は重量の対数に比例して減少する)のもとに微分方程式
を解いてみると、ちゃんとゴンペルツ関数が出て来る。
ここで
であり、に時間に関するパラメータが分離されているので、始点の問題を棚上げできる。
あとは適当に。
gompertz <- function(t, c1, a, b){ # RGR vs Time c1 * exp(-(exp(a-b*t)/b)) } t <- seq(0, 20, 0.5) a <- 0.3 b <- 0.25 c1 <- 10 g <- gompertz(t, c1, a, b) plot(t, g) r_t <- function(w, t) sum(coef(lm(w ~ t + I(t^2)))[2:3] * c(1, 2*t[2])) / w[2] r <- numeric(length(g)) r <- NA for(i in 1:length(g)){ r[i] <- r_t(g[c(i-1, i, i+1)], t[c(i-1, i, i+1)]) } gompertz2 <- function(t, c1, a, b){ # RGR vs Weight exp(a/b - exp(b*c1 - b*t)/b) } result <- lm(r ~ log(g)) result_coef <- coef(result) f <- function(c1) (sum(gompertz2(t, c1, result_coef[1], -result_coef[2]) - g))^2 c1 <- optimize(f, interval = c(-50, 100))$minimum points(t, gompertz2(t, c1, result_coef[1], -result_coef[2]), col = "red", type = "l")
植物の重量に上限があったとして、相対成長速度が上限と今の重量の差に比例するとしてみよう。
これは適当な定数をもってきて書き換えると、要するに相対成長速度が重量の一次式で表現できるということになる。
これを相対成長速度の定義に代入する。
そしてについて解くとこうなる。
これは実はロジスティック方程式である(定数を整理してみると良い)。
ととに適当な値を入れてプロットしてみよう。
a <- 0.4 b <- 0.2 c1 <- -15 t <- seq(0, 50, 1) d <- logistic(t, c1, a, b) plot(t, d)
そして、重量と相対成長速度の関係は当然だが一次式になっている。
r_t <- function(w, t) sum(coef(lm(w ~ t + I(t^2)))[2:3] * c(1, 2*t[2])) / w[2] r <- numeric(length(d)) r <- NA for(i in 1:length(d)){ r[i] <- r_t(d[c(i-1, i, i+1)], t[c(i-1, i, i+1)]) } plot(d, r)
点が不等間隔に並んでいるのも分かるだろう。重量の上限値でもって相対成長速度は0になるので、負の値になることはない。
ここからパラメータとを求めるのは簡単で、単に回帰すれば良い。またはoptimizeを使えば求められる。ただし、は曲線の立ち上がりを左右に移動させるだけなので、データのどこをと定めるかによって左右される。本質的なのはとだ。
plot(t, d) result <- lm(r ~ d) result_coef <- coef(result) result_coef # > 0.4074815 -0.2040751 f <- function(c1) (sum(logistic(t, c1, result_coef[1], -result_coef[2]) - d))^2 c1 <- optimize(f, interval = c(-100, 100))$minimum points(t, logistic(t, c1, result_coef[1], -result_coef[2]), col = "red", type = "l")
今回の例では相対成長速度と時間を切り離している。そのため、例えば測定データのセットが複数あった場合、相対成長速度vs重量というデータに変換してしまうことで複数の結果を束ねて解析することが容易になる。
複数の葉や果実の生長、肥大を測定した場合に、どこをとして測定を開始するのかというのは難しい問題である。似たような大きさのものを選んだとしても、特にサイズが小さいうちは指数関数的に生育が進むので、初期の誤差が生育曲線の位相の大きなずれとなって現れてくる。この問題を棚上げしたまま共通のパラメータを推定できるのは大きなメリットだろう。
ちなみにゴンペルツ曲線とロジスティック曲線はそもそも形が違う(ゴンペルツ曲線は対称ではない)ので、その点は注意が必要である。
前回はRGRを重量の対数の変化速度として求めたが、そもそもRGRの定義は
であったので、複数のデータポイントがあれば何らかの方法で平滑化して微分すれば、より精度良くRGRが推定できるだろう。例えば前後2点を含めた3点をラグランジュ補間して得た二次関数を微分するとか。
Rでラグランジュ補間をするパッケージは多分あるが面倒なので二次関数を当てはめてしまう。
r_t <- function(w, t) sum(coef(lm(w ~ t + I(t^2)))[2:3] * c(1, 2*t[2])) / w[2] r <- numeric(length(g)) r <- NA for(i in 1:length(g)){ r[i] <- r_t(g[c(i-1, i, i+1)], t[c(i-1, i, i+1)]) }
そしてあとは前回のヤツを多少いじって実行する。
result <- glm(r ~ t, family = gaussian(link="log"), start = c(1, 0)) result_coef <- coef(result) plot(t, g) f <- function(c1) (sum(gompertz(t, c1, result_coef[1], -result_coef[2]) - g))^2 c1 <- optimize(f, interval = c(0, 100))$minimum points(t, gompertz(t, c1, result_coef[1], -result_coef[2]), col = "red", type = "l")
すると少しはマシな推定値が出る。
ただ、繰り返しになるがゴンペルツ曲線を求めるなら最初からゴンペルツ曲線をデータに当てはめたほうが良い。例えばgに適当な誤差を与えて推定を繰り返してみると、マズさが分かるだろう。
特定のタイミングにおけるRGRの推定値がほしければ、ゴンペルツ関数を微分したものに推定したパラメータを渡せば良い。微分した結果を関数として使うにはderiv()を使って引数にfunc=TRUEを指定する。
dg <- deriv(~ c1 * exp(-(exp(a-b*t)/b)), "t", func = TRUE) attr(dg(t), "gradient") r <- attr(dg(t), "gradient")[,1] / gompertz(t, c1, a, b) plot(r)