salvageon

結城浩さんのcodeIQの問題 https://codeiq.jp/ace/yuki_hiroshi/q1215

フィードバックをいただきました。両方正解でした。
Rを使っていた人は2名だけのようでしたが，ビッグデータの取り扱いは統計ツールを日常的に使う人にも関係は深いと思うので，もっとR使いがいてもいいのではないでしょうか。私はgmpというパッケージで超長整数を扱って計算しました。
DB1は単調増加とわかった段階で，微分方程式の数値解析のオイラー法みたいな計算をしましたが，解説ではバイナリサーチをしていました。この場合はオイラー法のほうがAPIをたたく回数は減ると思いますが…

解答集
http://togetter.com/li/758302

以下にanswer.txtをさらします。（画面からはみ出る部分に改行を追加）

V406435859539156181269150751031
V1101943557675920722238136981003
ENV: 最終的に「R」。途中経過で電卓，Excel, 少しだけpython
POINT: DB1はキーが単調増加。DB2はキーを100bitの2進数で表し
いちばん下の1の桁を上下反転させたものに残りの桁を加えたもの
の順に並んでいる。
（キーQを割り切れる最大の2のべき乗を2^k(k>=0)とした時，
2^(99-k)+(Q-2^k)/2^(k+1)がindexになる）

Rのソースと実行結果


    

  • -
    • 

#DB=1
library(gmp)
db0 <- as.bigz(0) #小さいほう
db1 <- pow.bigz(2,100) - as.bigz(1) #大きいほう
target <- as.bigz("208050656559285601386927895421059705239114932023754")
 #クォーテーションマークでくくることが重要
VV <- ""
k <- 0

urltext <- paste("http://salvageon.textfile.org/?db=1&index=", db0, sep="")
dd <- read.table(file=url(urltext), header=F, sep=" ", colClasses = "character")
 #多倍長整数も文字列で取り込む
ddl <- data.frame(dd)
k <- k + 1
ky0 <- as.bigz(substr(dd[,3],2,100))
Sys.sleep(1.0)
urltext <- paste("http://salvageon.textfile.org/?db=1&index=", db1, sep="")
dd <- read.table(file=url(urltext), header=F, sep=" ", colClasses = "character")
ddl <- rbind(ddl, dd)
k <- k + 1
ky1 <- as.bigz(substr(dd[,3],2,100))
Sys.sleep(1.0)

while (VV == "") {
    dbx <- as.bigz(sub.bigz(target,ky0)/sub.bigz(ky1,ky0)*as.bigq(sub.bigz(db1,db0))+_
    as.bigq(db0))
 #均一の増加率だと仮定して探し出すキーのインデックスをあたりをつける
    urltext <- paste("http://salvageon.textfile.org/?db=1&index=", dbx, sep="")
    dd <- read.table(file=url(urltext), header=F, sep=" ", colClasses = "character")
    ddl <- rbind(ddl, dd)
    k <- k + 1
    kyx <- as.bigz(substr(dd[,3],2,100))
    dbx
    Sys.sleep(1.0)

    if (sign(sub.bigz(kyx,target)) == 1) {
        db1 <- dbx
        ky1 <- kyx
    } else if (sign(sub.bigz(kyx,target)) == -1) {
        db0 <- dbx
        ky0 <- kyx
    } else {
        VV <- dd[,4]
    }
}

    

  • -
    • 

> VV
[1] "V406435859539156181269150751031"
> ddl
  V1                              V2
1  1                               0
2  1 1267650600228229401496703205375
3  1       2565070352901388243030490
4  1       2565070352901388243030491
                                                          V3
1                             K19584500653053662232211568267
2 K102818053067020370282264996429296647301172110315939147734
3       K208050656559285601386927806221480835207699089545741
4       K208050656559285601386927895421059705239114932023754
                               V4                              V5
1 V830542486163201924733591581247 1267650600228229401496703205376
2 V764581954396429898637988092086 1267650600228229401496703205376
3 V400774553816906363998664980394 1267650600228229401496703205376
4 V406435859539156181269150751031 1267650600228229401496703205376

    

  • -
    • 

#最初の数個を見た限りでは単調増加とはいえkeyの階差数列はばらつきが
 大きいようでしたが，はさみうちをするつもりで書いたwhile文ではわず
 か2回目で解にたどり着きました。
#多倍長整数の計算でうっかりnumericにしてしまうと32bitで丸められて
 しまって正しい数字にならないというバグを探すのにてこずりました。
 targetをクォーテーションマークでくくっていなかったために変な数(bigz)
 になっていました。

    

  • -
    • 

> target <- as.bigz("208050656559285601386927895421059705239114932023754")
> target
Big Integer ('bigz') :
[1] 208050656559285601386927895421059705239114932023754
> target <- as.bigz(208050656559285601386927895421059705239114932023754)
> target
Big Integer ('bigz') :
[1] 208050656559285588701770946804922019477298484346880
> 

    

  • -
    • 



    

  • -
    • 

#DB=2
db0 <- as.bigz(0)
db1 <- pow.bigz(2,100) - as.bigz(1)
target <- as.bigz("2023636070998557444542586045")
VV <- ""

urltext <- paste("http://salvageon.textfile.org/?db=2&index=", db0, sep="")
dd<-read.table(file=url(urltext), header=F, sep=" ", colClasses = "character")
ddl <- data.frame(dd)

#規則性を発見する過程は略
for (i in 1:20) {
    dbx <- as.bigz(i)
    urltext <- paste("http://salvageon.textfile.org/?db=2&index=", dbx, sep="")
    dd<-read.table(file=url(urltext), header=F, sep=" ", colClasses = "character")
    ddl <- rbind(ddl,dd)
    Sys.sleep(1.0)
}

dbx <- add.bigz(pow.bigz(2,99),as.bigz(sub.bigz(target,as.bigz(1))/as.bigz(2)))
urltext <- paste("http://salvageon.textfile.org/?db=2&index=", dbx, sep="")
dd<-read.table(file=url(urltext), header=F, sep=" ", colClasses = "character")
ddl <- rbind(ddl,dd)
kyx <- as.bigz(substr(dd[,3],2,100))
VV <- dd[,4]

    

  • -
    • 

> VV
[1] "V1101943557675920722238136981003"
> ddl
   V1                             V2                               V3
1   2                              0                               K0
2   2                              1  K633825300114114700748351602688
3   2                              2  K316912650057057350374175801344
4   2                              3  K950737950171172051122527404032
5   2                              4  K158456325028528675187087900672
6   2                              5  K475368975085586025561263702016
7   2                              6  K792281625142643375935439503360
8   2                              7 K1109194275199700726309615304704
9   2                              8   K79228162514264337593543950336
10  2                              9  K237684487542793012780631851008
11  2                             10  K396140812571321687967719751680
12  2                             11  K554597137599850363154807652352
13  2                             12  K713053462628379038341895553024
14  2                             13  K871509787656907713528983453696
15  2                             14 K1029966112685436388716071354368
16  2                             15 K1188422437713965063903159255040
17  2                             16   K39614081257132168796771975168
18  2                             17  K118842243771396506390315925504
19  2                             18  K198070406285660843983859875840
20  2                             19  K277298568799925181577403826176
21  2                             20  K356526731314189519170947776512
22  2 634837118149613979470622895710    K2023636070998557444542586045
                                 V4                              V5
1   V866608106806207094188269706010 1267650600228229401496703205376
2   V179347330550907655201591936271 1267650600228229401496703205376
3   V967874566490056905763590096976 1267650600228229401496703205376
4    V21007428639505136878697302990 1267650600228229401496703205376
5   V499903880406975167914626368467 1267650600228229401496703205376
6   V137612525995897340648000428116 1267650600228229401496703205376
7  V1013924898829306772380170269614 1267650600228229401496703205376
8  V1252404645056427291182533085987 1267650600228229401496703205376
9   V520595999339570965993883078828 1267650600228229401496703205376
10 V1256891284967277336714134496424 1267650600228229401496703205376
11  V710302931554296640526866843240 1267650600228229401496703205376
12  V545375329928202174295039839589 1267650600228229401496703205376
13  V848644394970802704508212815012 1267650600228229401496703205376
14  V871991163218223896080356075737 1267650600228229401496703205376
15  V958747346069852079160957398348 1267650600228229401496703205376
16  V805754670819605302288978300850 1267650600228229401496703205376
17  V115917505363134018513900007273 1267650600228229401496703205376
18  V154823019080372587041596024019 1267650600228229401496703205376
19  V108729423996271630976080975567 1267650600228229401496703205376
20  V886213339588600673938994101775 1267650600228229401496703205376
21  V334205435694911673982823977556 1267650600228229401496703205376
22 V1101943557675920722238136981003 1267650600228229401496703205376


    

  • -
    • 

#最初の20個と最後の1個を見て，データの個数が2の100乗でindexとkeyが
 重複しないことから2進数の100桁の数字の桁の入れ替えで作られている
 ようだとわかりました。

実はRで解く前にwindowsの電卓（超長整数を扱えるようだ）を使っていろいろ試していたのですが，DB2では2進数と10進数の変換に以下のサイトが役立ちました。
http://hogehoge.tk/tool/number.html
ありがとうございます。

カレー問題

えー、たいへんひさしぶりの更新です。
年度末繁忙期の現実逃避に、結城浩さんのCodeIQのカレー問題をエクセルだけで解いてみました。

解答に書いたメモを転載。本当は実際のエクセルファイルの絵もあるといいのですが、とりあえず。

      • -

これならExcel（表計算ソフト）だけでも出来そうだなと取りかかりましたが，30分では終わらず1時間ぐらいかかってしまいました。

  • blendlist.txtをExcelで開く。スペース区切りで2列に分け，上に1行挿入して列名（spice1,spice2）を入れる。
  • （このままpivot集計をしてみたらspice1とspice2が交換可能であることに気づいた）
  • spice1のリスト[A2:A2183]の下のセル[A2184]にspice2[B2:B2183]をコピーして貼り付け，spice2の下[B2184]にspice1[A2:A2183]を追加。
  • Excel2003の「データ→ピボットテーブルとピボットグラフ　レポート」で新しいワークシート(Sheet1)に行タイトルspice1と列タイトルspice2，項目spice2の個数でピボット集計。対角線に対して対称で，blendlistに存在する組み合わせに1，ない場合は空白の表ができる。spiceの種類は128種であることがわかる。
  • ピボットテーブル[A3:DZ133](Sheet1)をコピーして新しいワークシート(Sheet2)の左上セル[A1]に形式を選択して貼り付け（数値）。あとの操作をやりやすくするため。
  • 出現頻度が一番高いspiceを探すため，「総計」列の最大値を探すと44が2つある。とりあえず"Minnlabic"を最初の「指標種」とする。
  • Sheet1の"Minnlabic"の行（58番目なのでピボットテーブルの行62）を表の下[Sheet1!行135]にコピーして貼り付け。
  • spiceの名前[A5:A132]をコピーして[A137]に値貼り付け。
  • [B137]に"=B5*B$135"を入力，コピーして[B137:DY264]に貼り付け。[DZ137]に"=sum(B137:DY137)"と入力。コピーして[DZ138:DZ264]に貼り付け。Minnlabicとの組み合わせでポイントになるspiceを共通に持っている数を数える。（28以上と3以下の両極端に分かれるぞ。しめしめ）
  • Sheet1の[DZ137:DZ264]をコピーしてSheet2の[EA3]のセルに数値で貼り付け。また[A2]のセルに形式を選択して貼り付け（数値，「行列を入れ替える」にチェック）。
  • Sheet2の表[A3:EA130]を選択，EA列の降順で並べ替え。続けて[A1:DY130]を選択。並べ替え→オプション「列単位」で行1の降順で並べ替え。
  • この時点で表の対角線に対する"1"の並び方でblendlistが大きく3つのグループになることがわかる。
  • "Minnlabic"と相性のよいスパイスを含むブレンドが28以上となるスパイスのグループをグループ1とする。グループ1に含まれないスパイスについて，例えば"Anewatry"を次の「指標種」に選び，同様の処理を繰り返すことが可能。
  • しかし，並べ替えた表を見ると，"Minnlabic"と相性のよいスパイスを含むブレンドが0のスパイスの中でLoveaniic, Masema, Peregrama, Rinawatry, Xuomaの5種が，「"Minnlabic"と相性のよいスパイスを含むブレンドが1〜3」のスパイスと相性がよく，残りは別のグループを形成することがわかるので，手作業で5種を"Minnlabic"と相性のよいスパイスを含むブレンドが0のスパイスの中でで最初に来るように並べ替えを行った。（行選択，カットして貼り付け。列選択，カットして貼り付け）
  • Group1（"Minnlabic"と相性のよいスパイスを含むブレンドが28以上），Group2（"Minnlabic"と相性のよいスパイスを含むブレンドが1〜3＋上記5種），Group3（"Minnlabic"と相性のよいスパイスを含むブレンドが0で上記5種を除く）とすると，それぞれのGroupを含む皿のカレーのブレンドのポイントは以下のようになる。
Group1 Group2 Group3
Group1 936
Group2 12 775
Group3 0 11 448
  • この際，2皿のポイントの合計を，3つのグループのうちどの1つを別の皿にするかで検討すると，Group3を別の皿にする場合がブレンドリスト中の11の組み合わせだけが達成されないことになり，最大である。
  • よってGroup1とGroup2に含まれるスパイスの名前（並べ替えた表の[B2:CQ2]）を別のエクセルに貼り付けて昇順に並べ替え，テキストファイルで保存。
  • 保存したテキストファイルのタブをスペースに置換。

数値標高モデルデータを実際に使ってみる(2)

GIS Omap

手順

  1. 基盤地図情報2500のダウンロード
  2. シェープファイルへの変換（公共測量ビューア・コンバータ）
  3. 基盤地図情報数値標高モデルのダウンロード
  4. GeoTiffへの変換（基盤地図情報DEMコンバータVer1.4）
  5. QGISでの表示（緯度経度）
  6. 等高線シェープファイルの作成
  7. 計曲線の色分け
  8. 平面直角座標への変換
  9. GPSログの読み込みとシェープファイルへの変換
  10. OCADへのシェープファイルの読み込み（平面直角座標の設定）
  11. OCADでの編集

（続く）

数値標高モデルデータを実際に使ってみる

GIS Omap

オリエンテーリング地図(Omap)は，国際オリエンテーリング連盟の定める規定（http://orienteering.org/resources/mapping/）に従って作成された，ナビゲーションのための地図です。Omapは都市計画図などを基礎図として，トレイル，微地形，林内の通行可能度など，地図を持って山中を走るときに必要な情報を現地踏査により書き加えて作ります。

比較

実際に基盤地図情報の数値標高モデル（http://fgd.gsi.go.jp/download/）を使った地図がOmapの「調査原図」としてどのように使えるか，一つの仮想例としてすでにOmapがあり，各種の数値標高モデル（5mメッシュ標高［航空レーザ測量，写真測量］，10mメッシュ標高［1:25000地形図ベース］）がそろっている「京都地区」で例示してみます。

Omap「奥大文字・山紫水明東山」 http://www.o-news.net/wr/2010/100530.jpg
この地図を使ったオリエンテーリング大会の記事 http://www.o-news.net/2010/06/z8.php


Omap「奥大文字山紫水明東山」の一部分。磁北方向が上になっているため，下の地図と比較すると約7°左へ回転している。



5mメッシュ標高からGDAL（d:id:kooi:20110713）で5m間隔の等高線を作成し，基盤地図情報縮尺レベル2500の道路縁，建物，水涯線と合わせて表示したもの。等高線が赤茶色の部分が「航空レーザ測量」DEM5A，薄茶色が「写真測量」DEM5Bのデータ部分。基盤地図情報では，同じ地区のDEM5AとDEM5Bは重複しては公開されていない。

基盤地図情報で全国のデータが作成されている10mメッシュ標高から，上記と同様に5m間隔の等高線を作成して表示したもの。

いかがですか。等高線による地形表現をOmapと比較してみると，下図のA，B地点では微小な沢地形が表現されていて航空レーザ測量の精度の高さをうかがい知ることができます。一方，C，D地点ではむしろ10mメッシュ標高の等高線の方が現実の地形に近いと考えられ，写真測量による5mメッシュDEMは（場所によっては）Omapの原図としての精度が10mメッシュ標高より劣ると考えられる場合もあるようです。

数値標高モデルを使った詳細地図

GIS

この項はオリエンテーリング地図の作製も念頭に置いて書いています。
基盤地図情報2500には等高線情報は含まれていません。そこで，数値標高モデル（5mメッシュ）から等高線を作成し，基盤地図情報2500とあわせて野外調査の現地野帳として使用できる詳細地図をつくります。ただし基盤地図情報2500はまだ公開されていないエリアも多く，数値標高モデルも10mメッシュは全国で作成されていますが，5mメッシュはまだ一部です。

必要なソフトウェア

GISはオープンソースのQGIS（http://www.qgis.org/wiki/Download）を使用しています（過去エントリ参照）。今回はプラグインのGdalToolsを使用するのでプラグインマネージャで使用可にしておいてください。
数値標高モデルをGeoTiff形式に変換するには，三匹のウリボウさんの「基盤地図情報DEMコンバータVer1.4」（http://space.geocities.jp/bischofia_vb/）を使用させていただきます。
基盤地図情報2500のシェープファイルへの変換は，公共測量ビューア・コンバータ（http://psgsv.gsi.go.jp/koukyou/public/sien/pindex.html）を使います。

データ（基盤地図情報）

基盤地図情報ダウンロードサービス（http://fgd.gsi.go.jp/download/）から，基盤地図情報2500と数値標高モデルの必要な地域のファイルをダウンロードします。

作業手順

1) 基盤地図情報2500は市区町村単位でダウンロードできるので必要な地域をダウンロードし，「公共測量ビューア・コンバータ」でシェープファイルに変換します。
2) 数値標高モデルは（世界測地系）2次メッシュ単位でダウンロードできるので，必要な地域のファイルをダウンロードし，「基盤地図情報DEMコンバータ」で開いて，GeoTiff形式で保存します。
3) QGISで，新しいプロジェクト（座標系はJGD2000）にDEMをコンバートしたGeoTiff形式のファイルを追加し，メニューの「ラスタ」（GdalTools）の「等高線」で，等高線のシェープファイルを作成します。

4) できたシェープファイルをQGISの新しいプロジェクト（座標系は平面直角座標で，オンザフライ座標変換をonにする）に追加します。
5) 必要なシェープファイルを「ベクタ→データマネージメントツール→新しい投影法へエクスポートする」で平面直角座標に変換します。
6) オリエンテーリング地図（OCAD使用 http://www.ocad.com/en/index.htm）へは平面直角座標系に変換したシェープファイルをインポートし，ファイル（レイヤ）ごとに適切な記号をあてます。（測量法で定められた利用申請はしてください）

参考

GDALを使って、基盤地図情報の標高データから等高線のshpファイルを作成する方法 - 自然環境保全のための周辺技術 http://d.hatena.ne.jp/tmizu23/20100215/1266229611
基盤地図情報（標高）をGDALやGRASSで使いたい - 地図はたいへん http://d.hatena.ne.jp/vec2ras/20100107/1262833221

基盤地図情報25000の等高線で計曲線・主曲線・補助曲線を描き分ける

GIS

自然環境調査などで使う印刷用3次メッシュマップのQGISを使った作り方（http://d.hatena.ne.jp/kooi/20110619）の補足です。
QGISの地図プロジェクトに読み込んだベクタレイヤの記号は今のところ「共通シンボル」で設定されています。これを「固有値」を使って属性によって記号を変化させるという話です。
基盤地図情報25000の等高線データ(Cntr.shp)の属性ファイル(Cntr.dbf)には「標高」という列があります。この標高フィールドを使って等高線を色分けしてみます。プロパティのシンボルで凡例タイプに「固有値」を選び，分類フィールドを「標高」にします。「分類」のボタンを押すと自動で全ての属性値に個別の凡例をつけることができるようになります。（「標高」のフィールドは数値ではなくテキストでデータが入っています。）

しかしこれを1つ1つ設定するのはかなり面倒です。10，20，30，40，60…と主曲線にあたる標高は「クラスを削除」してデフォルトで一括指定するとしても，50mおきの計曲線だけでも20本以上あります。緩傾斜のところで5mごとに入ってくる補助曲線も何本もあり，それぞれ定義していたらやはり大変です。そこで，ちょうど属性ファイルの「表示区分」フィールドは空欄になっているようなので，ここに等高線の種類を入れて表示区分で凡例を分類することにします。

  • QGISの属性テーブルを編集モードにしたときに使える「フィールド計算機」では，まだ複雑な計算は難しいようです。http://www.osgeo.jp/user_guide/user_guidech3.html#x10-1750003.7
  • Excel2003まではdbf形式を扱うことができたので，Excelで開いて表示区分のG2セルに以下の式を入れて列へコピーすることが考えられます。しかし，Excel2003までは65535行しか扱えないので，福岡県全部（約29万行）だとオーバーフローします。Excel2007ではdbf形式で保存できないのでこの方法はできません。基盤地図情報閲覧・コンバートソフトでシェープファイルにするときに少しずつ（ダウンロードしたzipファイル1つずつ）別のシェープファイルに変換すればExcel2003でオーバーフローせずにできるでしょう。

=if(right(I2,1)="5"),"2",if(or(right(I2,2)="00",right(I2,2)="50"),"1","0"))

今回は，約29万行のdbfファイルを一発で扱うために，またしても統計解析ソフトウェアRを使います。コードは以下のようになります。dbfはいちおう別名で保存していますので，できたdbfファイルを確認してから，名前をCntr.dbf→Cntr_ori.dbf，Cntr_R.dbf→Cntr.dbfというように変更してください。

library(foreign)
CntrType <- function(x) {
y <- as.numeric(as.character(x))
ifelse(floor(y/50)==y/50,"1",ifelse(floor(y/10)==y/10,"0","2"))
}
dd <- read.dbf("C:/GIS/Cntr.dbf")
dd$表示区分 <- CntrType(dd$標高)
write.dbf(dd, "C:/GIS/Cntr_R.dbf")

Cntr.dbfを置き換えたら，表示区分のフィールドに主曲線は"0"，計曲線（50mごと）は"1"，補助曲線(5m）は"2"が入っているので，それぞれ適切な記号を定義します。

  • R（参考RjpWiki http://www.okada.jp.org/RWiki/）のベクトル演算は強力で，上記のコードではデータの読み込みや書き出しに比べたらあっという間に計算は終わります。関数中でif文ではなくifelseを使ってベクトル演算できるようにしておくのがミソです。
  • 個人的には一番手間取ったのはdbfを読み込んだデータフレームで標高の数値がただのテキストではなくてfactorだったので，as.characterをかまさないと期待した結果を返してくれなかったことでした。

このようにして作った地図（pdfに出力したもの）の一部です。

インテルメッツォ：地図データを使うのに申請は必要か

GIS

国土地理院の地図を使うのには「測量法」できめられた申請が必要ですが，私的な利用等へは申請不要となっています。書籍やインターネットには，限られた量の掲載であれば「出所の明示」で申請不要とされています。
基盤地図情報に関するFAQ。5番目の項目が地図の複製に関することです。
http://www.gsi.go.jp/kiban/faq.html#52
そこでリンクされている「申請フロー」（http://www.gsi.go.jp/LAW/2930-flow.html）を見ると，以下のように書かれています（要約）。

  • 私的な使用，教育機関での使用は申請不要
  • 学術論文に掲載する場合，刊行物に少量の地図を掲載（挿入）する場合は「出所の明示」で申請不要
  • 使用目的により「複製」と「使用」を区別
  • 「複製」の場合，サークル内での使用や特定の者へ提出する報告書では申請不要

インターネットへ公開する場合は，300×400ピクセル以下の画像であれば「少量の地図を挿入」と見なされるようです。それ以上1画面範囲内（スクロール不要のサイズ）であれば5枚以内とされています。
http://www.gsi.go.jp/LAW/2930-qa.html
はてなfotofileは画像の長辺が450ピクセルになります。私がこれまであげてきたGISの使い方記事の画像は地図面だけでいえば300×400ピクセル以下のものがほとんどですので，「出所の明示」で申請不要と思っています。