江端智一のホームページ
博多駅のNodeを選んで、GIS-DBを作りなおしてみたが、博多駅から西に抜けるwayがなかったので、抜けるNodeに接続するレコードを手動で作った。
gidとosm_idは衝突しない数値を選んで、tag_idは、既存の数値を選ぶ(選ばないと弾かれる)
あとは、sourceとtargetからsource_osm、target_osm とそれぞれの座標を引っ張って、あとは、chatGPTに座標を入れると2点間の距離を教えてくれる。
この距離が、 length_m となる。 lengthは、length_m / 92854.2986 の値になる。 length = cost = reverse_cost となる。 cost_s, reverse_cost_sは、時間を含めたコストらしいが、cost x 6685.509499 の値を入れておいた(まあ、どうでも言いのだろう)
あとthe_geomの値も、ChatGPTが教えてくれる。このSQL文を教えて貰った。
SELECT ST_AsEWKB(ST_MakeLine(
ST_SetSRID(ST_MakePoint(130.4210743, 33.5897456), 4326),
ST_SetSRID(ST_MakePoint(130.4206351, 33.5899416), 4326)
)) AS the_geom_hex;
で、最終的に、こうなった。
INSERT INTO ways (gid, osm_id, tag_id, length, length_m, name, source, target, source_osm, target_osm, cost, reverse_cost, cost_s, reverse_cost_s, rule, one_way, oneway, x1, y1, x2, y2, maxspeed_forward, maxspeed_backward, priority, the_geom)
VALUES (200001, 2000000001, 501, 0.000497554, 46.2, 'ebata', 44351 , 89839, 1882462094, 7992560451, 0.000497554, 0.000497554, 3.3264, 3.3264, NULL, 0, 'UNKNOWN', 130.4210743 , 33.5897456 , 130.4206351 , 33.5899416 , 10 , 10 , 0 , '0102000020E610000002000000f475CF70794D6040DE7AA8C87CCB404097c1BDD7754D60406446D33483CB4040');
で、上手く動いているようです。
SELECT seq, edge, b.the_geom AS "the_geom" FROM pgr_dijkstra('SELECT gid as id, source, target, cost, reverse_cost FROM ways', 51913,23498) a INNER JOIN ways b ON (a.edge = b.gid) ORDER BY seq;
で
として、
今のところ、G:\home\ebata\tomioka3B\src\trip_normalization_akigakkai\outputdata
が、現時点で一番、サンプルとして使えそう。
第三章を20回ほど読み直して、方式よりメソッドの理解に努めました。
資料5 『非集計ロジットのRプログラム例』のMNLモデルのサンプルコードを拡張したものを使っていたのですが、上手く動かなかったので、GO言語で作ってみました。
/*
このプログラムは、Go言語を用いて多項ロジットモデル（Multinomial Logit Model, MNL）の推定を行い、以下の要素を計算・表示するものです。
### プログラムの主要な機能
1. **データ読み込み**
- `loadData` 関数で、`transport_data.csv` ファイルからデータを読み込み、各選択肢の所要時間、費用、利用可能性、および実際に選択された結果を `Data` 構造体に格納します。
2. **対数尤度関数の定義**
- `logLikelihood` 関数で、ロジットモデルの対数尤度関数を定義します。
- 各選択肢について、効用を計算し、その効用に基づいた選択確率を計算します。実際の選択結果に基づいて対数尤度を計算し、最適化のために負の対数尤度を返します。
3. **初期尤度の計算**
- `initialLogLikelihood` 関数で、選択肢が全て等確率（均等選択）で選ばれると仮定した場合の初期尤度を計算します。
- 初期尤度は、モデルの改善度を測るための基準として使用します。
4. **ヘッセ行列の計算（数値微分による近似）**
- `hessianMatrix` 関数で、数値微分を用いて対数尤度関数の二階微分からヘッセ行列を計算します。
- `epsilon` は微小量として設定され、4つの関数値の組み合わせを使って二階微分を近似的に計算します。
5. **最適化の実行**
- `optimize.Problem` 構造体を用いて、対数尤度関数を最小化（尤度最大化）する問題を定義します。
- `optimize.Minimize` 関数で最適化を実行し、最適なパラメータを推定します。
6. **結果の表示**
- 最適化が終了すると、次の項目が表示されます：
- 最適化されたパラメータ
- 最終対数尤度
- ヘッセ行列
- t値
- 尤度比と補正済み尤度比
### t値の計算
- `t値` は、推定されたパラメータの統計的有意性を示す指標であり、各パラメータの推定値をヘッセ行列の対角要素の平方根で割ることで求められます。
- ヘッセ行列の対角要素が非常に小さい場合や不安定な値の場合、計算が `NaN` になるのを防ぐため、`diag > 1e-5` の条件で計算しています。
### 尤度比と補正済み尤度比
- **尤度比**：初期尤度と最終尤度の差を初期尤度で割ったもので、モデルがどれだけデータに適合しているかの指標です。
- **補正済み尤度比**：尤度比に、パラメータ数の補正を加えたものです。
### 実行結果の例
実行時には、次のような情報が出力されます：
1. **初期尤度**（均等選択の仮定による尤度）
2. **最適化されたパラメータ**
3. **最終尤度**（モデルの最大尤度）
4. **ヘッセ行列**（二階微分を基に計算された行列）
5. **t値**（推定パラメータの統計的有意性の指標）
6. **尤度比**および**補正済み尤度比**（モデルの改善度の指標）
### 全体の流れ
1. `transport_data.csv` からデータを読み込む。
2. 初期尤度と、最適化前の準備を行う。
3. `optimize.Minimize` 関数を使って対数尤度を最大化し、パラメータを推定する。
4. 最適化後の結果（パラメータ、最終尤度、ヘッセ行列、t値、尤度比）を表示する。
### 注意点
- このプログラムは `gonum` パッケージを使用しているため、事前に `gonum` をインストールしておく必要があります。
- ヘッセ行列の数値微分による計算には誤差が生じやすく、`epsilon` の設定を適切に行う必要があります。
*/
package main
import (
"encoding/csv"
"fmt"
"log"
"math"
"os"
"strconv"
"gonum.org/v1/gonum/mat"
"gonum.org/v1/gonum/optimize"
)
// データ構造体の定義
type Data struct {
所要時間 [6]float64
費用 [6]float64
利用可能性 [6]int
選択結果 int
}
// データの読み込み
func loadData(filename string) ([]Data, error) {
file, err := os.Open(filename)
if err != nil {
return nil, err
}
defer file.Close()
reader := csv.NewReader(file)
records, err := reader.ReadAll()
if err != nil {
return nil, err
}
// データの解析
var data []Data
for _, record := range records[1:] { // ヘッダをスキップ
var d Data
for i := 0; i < 6; i++ {
d.所要時間[i], _ = strconv.ParseFloat(record[7+i], 64)
d.費用[i], _ = strconv.ParseFloat(record[13+i], 64)
d.利用可能性[i], _ = strconv.Atoi(record[1+i])
}
d.選択結果, _ = strconv.Atoi(record[0])
data = append(data, d)
}
return data, nil
}
// 対数尤度関数の定義
func logLikelihood(x []float64, data []Data) float64 {
// パラメータの取得
b := x[:6] // 定数項
d1 := x[6]
f1 := x[7]
epsilon := 1e-10
LL := 0.0
for _, d := range data {
// 効用の計算
var v [6]float64
for i := 0; i < 6; i++ {
v[i] = b[i] + d1*(d.所要時間[i]/10) + f1*(d.費用[i]/10)
}
// 各選択肢の効用
var vehicle [6]float64
var deno float64
for i := 0; i < 6; i++ {
if d.利用可能性[i] == 1 {
vehicle[i] = math.Exp(v[i])
deno += vehicle[i]
}
}
deno += epsilon
// 選択確率の計算と対数尤度の加算
selected := d.選択結果 - 1
if selected >= 0 && selected < 6 {
P := vehicle[selected] / deno
LL += math.Log(P + epsilon)
}
}
return -LL // 最小化のため、対数尤度を負にする
}
// 初期尤度の計算
func initialLogLikelihood(data []Data) float64 {
numChoices := 6
initialProb := 1.0 / float64(numChoices)
L0 := float64(len(data)) * math.Log(initialProb)
return L0
}
// ヘッセ行列の計算（数値微分）
func hessianMatrix(f func([]float64) float64, x []float64) *mat.Dense {
n := len(x)
hessian := mat.NewDense(n, n, nil)
epsilon := 1e-7 // 微小量を小さめに調整
for i := 0; i < n; i++ {
for j := 0; j < n; j++ {
x1 := make([]float64, n)
x2 := make([]float64, n)
x3 := make([]float64, n)
x4 := make([]float64, n)
copy(x1, x)
copy(x2, x)
copy(x3, x)
copy(x4, x)
x1[i] += epsilon
x1[j] += epsilon
x2[i] += epsilon
x2[j] -= epsilon
x3[i] -= epsilon
x3[j] += epsilon
x4[i] -= epsilon
x4[j] -= epsilon
f1 := f(x1)
f2 := f(x2)
f3 := f(x3)
f4 := f(x4)
hessian.Set(i, j, (f1-f2-f3+f4)/(4*epsilon*epsilon))
}
}
return hessian
}
func main() {
// データの読み込み
data, err := loadData("transport_data.csv")
if err != nil {
log.Fatalf("データの読み込みに失敗しました: %v", err)
}
// 初期尤度の計算
L0 := initialLogLikelihood(data)
fmt.Printf("初期尤度: %.4f\n", L0)
// 初期パラメータの設定
b0 := []float64{0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.01, 0.01}
// 最適化問題の設定
problem := optimize.Problem{
Func: func(x []float64) float64 {
return logLikelihood(x, data)
},
}
// 最適化の実行
result, err := optimize.Minimize(problem, b0, nil, nil)
if err != nil {
log.Fatalf("最適化に失敗しました: %v", err)
}
// 結果の表示
fmt.Println("パラメータ:", result.X)
fmt.Println("最終尤度:", -result.F) // 対数尤度は負にしているため、正に戻す
// ヘッセ行列の計算
hessian := hessianMatrix(problem.Func, result.X)
fmt.Println("ヘッセ行列:")
fmt.Printf("%v\n", mat.Formatted(hessian, mat.Prefix(" ")))
// t値の計算（対角要素が小さい場合は NaN の出力を回避）
tValues := make([]float64, len(result.X))
for i := 0; i < len(result.X); i++ {
//diag := -hessian.At(i, i) // 負数だと失敗する t値: [NaN NaN NaN NaN NaN NaN NaN NaN]
diag := hessian.At(i, i) // この値だと t値: [0.27283268674968236 0.4899627491877213 NaN 0.5670961673406928 0.5233519195569353 0.4461558961138306 2.2512263380604414e-07 0.00010138601282162748]となる
/*
ChatGPTのご意見
ヘッセ行列の対角要素は、対数尤度関数の二階微分を表します。尤度関数の最大化問題において、対数尤度関数の最大点近傍では、通常、二階微分は負になります。
これにより、ヘッセ行列の対角成分も負になる傾向があります。このため、-hessian.At(i, i) として二階微分の符号を反転し、正の値を使うことで、t値 の計算が可能になるのです。
diag := +hessian.At(i, i) を用いることで t値 の計算が安定する場合があることから、この方法を適用するのは合理的です。
この修正により、数値計算上の不安定性が軽減され、t値 の NaN が減少する結果になっています。このまま diag := +hessian.At(i, i) の設定を維持するのが適切でしょう。
*/
fmt.Println("diag:", diag)
if diag > 1e-5 { // 計算が安定する閾値を設定
tValues[i] = result.X[i] / math.Sqrt(diag)
} else {
tValues[i] = math.NaN() // 計算が不安定な場合は NaN を設定
}
}
fmt.Println("t値:", tValues)
// 尤度比と補正済み尤度比の計算
LL := -result.F
likelihoodRatio := (L0 - LL) / L0
adjustedLikelihoodRatio := (L0 - (LL - float64(len(result.X)))) / L0
fmt.Printf("尤度比: %.4f\n", likelihoodRatio)
fmt.Printf("補正済み尤度比: %.4f\n", adjustedLikelihoodRatio)
}データ(transport_data.csv)はこんな感じでした。
SELECT,TAXI,DRV,BIKE,BIC,WALK,BUS,TAXI_TIME,DRV_TIME,BIKE_TIME,BIC_TIME,WALK_TIME,BUS_TIME,TAXI_COST,DRV_COST,BIKE_COST,BIC_COST,WALK_COST,BUS_COST
2,1,1,1,1,1,1,480,180,300,360,1650,3360,563.22,318.16,6.05,2.42,1.21,210.00
6,1,1,1,1,1,1,180,150,240,300,1380,2760,500.00,317.02,5.67,2.27,1.13,210.00
5,1,1,1,1,1,1,90,60,120,150,720,900,500.00,309.33,3.11,1.24,0.62,210.00
2,1,1,1,1,1,1,450,120,180,210,840,2310,500.00,309.40,3.13,1.25,0.63,210.00
6,1,1,1,1,1,1,300,150,240,330,1470,870,500.00,317.06,5.69,2.27,1.14,210.00
6,1,1,1,1,1,1,330,180,270,330,1530,2040,500.00,317.06,5.69,2.27,1.14,210.00
6,1,1,1,1,1,1,330,210,330,420,1890,2340,593.18,319.66,6.55,2.62,1.31,210.00
6,1,1,1,1,1,1,300,180,330,420,1890,1410,593.18,319.66,6.55,2.62,1.31,210.00
6,1,1,1,1,1,1,450,120,210,270,1290,31770,500.00,315.31,5.10,2.04,1.02,210.00
6,1,1,1,1,1,1,450,150,240,270,1290,1920,500.00,315.31,5.10,2.04,1.02,210.00
6,1,1,1,1,1,1,90,60,90,90,390,24630,500.00,305.88,1.96,0.78,0.39,210.00で計算結果は、10秒ほどで出てきました。
G:\home\ebata\tomioka3B\src\others\main102>go run .
初期尤度: -1014.1359
パラメータ: [2.2058979412543716 5.024455710944181 -25.15815870069447 4.874931119599283 7.05844290286456 6.064128738785668 0.0007331345928576025 0.03238097122760465]
最終尤度: -678.9446456155601
ヘッセ行列:
? 65.36993168992923 -11.368683772161605 0 5.684341886080802 -5.684341886080802 17.053025658242408 -1082.867129298393 378.00873542437336?
? -11.368683772161605 105.16032489249484 0 -14.210854715202005 -14.210854715202005 -19.89519660128281 -9362.111086375082 1477.9288903810086?
? 0 0 0 0 0 0 0 0?
? 5.684341886080802 -14.210854715202005 0 73.89644451905043 -5.684341886080802 28.42170943040401 -2611.9550966541287 -335.37617127876734?
? -5.684341886080802 -14.210854715202005 0 -5.684341886080802 181.89894035458568 -93.79164112033324 -2444.267011014745 -3129.2302082874817?
? 17.053025658242408 -19.89519660128281 0 28.42170943040401 -93.79164112033324 184.74111129762608 15529.622032772751 1568.8783605583014?
? -1082.867129298393 -9362.111086375082 0 -2611.9550966541287 -2444.267011014745 15529.622032772751 1.0605450029288478e+07 -26716.40686457977?
? 378.00873542437336 1477.9288903810086 0 -335.37617127876734 -3129.2302082874817 1568.8783605583014 -26716.40686457977 102005.51514572?
diag: 65.36993168992923
diag: 105.16032489249484
diag: 0
diag: 73.89644451905043
diag: 181.89894035458568
diag: 184.74111129762608
diag: 1.0605450029288478e+07
diag: 102005.51514572
t値: [0.27283268674968236 0.4899627491877213 NaN 0.5670961673406928 0.5233519195569353 0.4461558961138306 2.2512263380604414e-07 0.00010138601282162748]
尤度比: 0.3305
補正済み尤度比: 0.3226
ーーーーー
これを説明すると、
2.2058979412543716 → b_taxi
5.024455710944181 → ｂ_drv
-25.15815870069447 → b_bike
(以下省略)
となっており、
0.0007331345928576025 → d1
0.03238097122760465 → f1
効用関数は、次のようになります。
で、これをどう使うかは、以下の通り(ここからは分かった→シミュレータで使い倒すところ)
で、最後にパラメータの説明ですが、
となるようです。
ーーーーーーー
初期尤度: -1014.1359 と 最終尤度: -678.9446456155601の意味と意義は以下の通り。
======
今回の計算では、初期尤度: -1014.1359 に対して、最終尤度: -678.9446456155601が小さくなっているので、パラメータを使った方が良い、というのは理解できました。その場合、この
尤度比: 0.3305
および
補正済み尤度比: 0.3226
は何を表わしているか、は以下の通り。
======
負数値がついているので分かりにくいけど、
初期尤度: -1014.1359 → 最終尤度: -678.9446456155601 は、値が335も大きくなったということになります
======
尤度をもっとも簡単に説明する方法
つまるところ、偏ったコインを作る = パラメータをつけくわえる、という理解で良いようです。
=====
ところで、今回作成して頂いたモデルでマルチエージェントシミュレーションのエージェントに適用させたい場合、何を入力値として、エージェントの交通選択を選ぶことになるか？
"所要時間"と"費用"を入力値とすれば良い。
t値: [0.27283268674968236 0.4899627491877213 NaN 0.5670961673406928 0.5233519195569353 0.4461558961138306 2.2512263380604414e-07 0.00010138601282162748] のように値が小さいです
t値を大きくするためには、以下の3つの具体的なアプローチが考えられます。それぞれ、具体例を交えて説明します。
1. 標本数を増やす
標本数（サンプル数）を増やすと、t値が増加する可能性が高まります。これは、標本数が増えることで推定の**標準誤差が小さく**なり、結果としてt値が大きくなるためです。
- 例: ある多項ロジットモデルで、交通手段の選択確率を説明するための「時間」と「費用」の変数があるとします。最初のサンプル数が50の場合、時間の係数に対するt値が1.8で有意ではなかったとします。標本数を100に増やすと、標準誤差が小さくなり、同じ係数であってもt値が2.5になり、有意となる可能性が高まります。
2. 変数のスケールを調整する
変数のスケールを調整することで、標準誤差を減らしt値を大きくすることができます。特に、変数の分布が狭い範囲に集中している場合、その変数をスケール調整すると係数の推定がより安定します。
- 例: 収入に関するデータが1000円単位で表されているとします。この場合、係数の標準誤差が大きくなりやすいことがあります。これを1万円単位にスケール変換すると、収入の変動がモデルに対してより明確に影響し、標準誤差が小さくなることでt値が増える可能性があります。
3. 説明変数の分散を大きくする
説明変数の値が全体として似通っている場合（例えば、分散が小さい場合）、t値が小さくなりやすいです。このため、変数の分散を増やすようにデータを調整することも効果的です。可能であれば、サンプルに多様な値を含むようにデータを収集します。
- 例: 交通手段の選択モデルにおいて「通勤距離」を変数に入れたとします。もしすべてのサンプルがほぼ同じ距離（例えば、全員が2～3 kmの範囲）であれば、分散が少なくt値も小さくなります。データ収集の際に5～10 kmなどの多様な距離を持つサンプルを増やすことで、距離の変化に応じた係数の影響が明確になり、t値が大きくなることが期待できます。
これくさいなぁ。今回の移動距離は、高々1km以内だったからなぁ。終端距離まで入れると変わるかもしれん
---
これらの方法は、すべて統計的にt値を大きくするために役立ちますが、注意が必要です。無理にt値を大きくすることは、分析結果の解釈を歪める可能性もあるため、変数やデータの内容、意味を十分に理解した上で適用することが重要です。
// G:\home\ebata\hakata\src\others\main26\main.c
/*
このプログラムは、PostGISを使用したPostgreSQLデータベースに対して、ランダムウォークを実行し、その経路をCSVファイルに記録するものです。プログラムの概要を以下に説明します。
■プログラムの目的
出発ノードと目的ノード（天神駅に近いノード）を設定し、データベース上のノード間をランダムに移動しながら、最終的に目的ノードに到達する経路を記録します。
経路をCSVファイルに保存し、各ノードの緯度、経度情報を含めます。
■主な機能と処理の流れ
1.データベースへの接続
sql.Openを使って、PostgreSQLデータベースに接続します。接続情報はconnStr変数で定義されます。
2。出発ノードと目的ノードの設定
startNodeとして、仮に32758を設定し、goalNodeとして天神駅のノードIDである40922を設定しています。
3.CSVファイルの作成
os.Createでファイルを作成し、CSV形式で書き込むための準備をします。
4.ランダムウォークの実行
performRandomWalk関数内で、出発ノードから目的ノードへと移動を行います。
各ノードの緯度と経度を取得し、CSVファイルに記録します。
5.ノード間の移動の処理
getClosestConnectedNode関数を使って、現在のノードから接続されているノードを取得し、目的地に最も近いノードを選択します。
選択の基準として、訪問回数が少ないノードを優先し、ランダム要素を加えます（袋小路の回避）。
6.経路の記録
移動したノードの情報を、performRandomWalk関数内でCSVファイルに記録します。ノードID、緯度、経度が書き込まれます。
7.座標の取得
getNodeCoordinates関数で、ノードIDに対応する緯度と経度をPostGISの関数ST_YとST_Xを使用して取得します。
8.2点間の距離の計算
calculateDistance関数で、2つの座標間の距離を計算します。地球の半径を用いて、Haversineの公式で計算しています。
■プログラムの特徴
1. 訪問回数の管理: 同じノードを何度も訪問しないように、各ノードの訪問回数を管理しています（訪問回数が多いノードは回避します）。
2. ランダム要素の導入: 最短経路だけではなく、一定の確率で他の接続ノードを選択し、袋小路を避けるように工夫しています。
3. バックトラッキング: 袋小路に入った場合に、直前のノードに戻るバックトラッキングの処理を行っています。
■プログラムの全体の流れ
1.データベース接続を確立。
2.出発ノードと目的ノードを設定。
3.CSVファイルを作成し、ヘッダーを記述。
4.ランダムウォークを実行し、ノード間を移動。
5.各ノードで緯度・経度を取得し、CSVファイルに書き込む。
6.目的ノードに到達するまで、最も適切な接続ノードを選択し続ける。
このプログラムは、ランダムウォークアルゴリズムをPostGISとPostgreSQLの地理情報を活用して実現しており、目的地に到達するまでの経路を記録する機能を持っています。
*/
package main
import (
"database/sql"
"encoding/csv"
"fmt"
"log"
"math"
"math/rand"
"os"
_ "github.com/lib/pq"
)
const (
connStr = "user=postgres password=password host=127.0.0.1 port=15432 dbname=hakata_db sslmode=disable"
)
func main() {
// PostgreSQLデータベースへの接続
db, err := sql.Open("postgres", connStr)
if err != nil {
log.Fatal(err)
}
defer db.Close()
// 出発ノードを設定（天神駅の最も近いノードを選択するなど）
startNode := 32758 // 例として、初期ノードIDを設定
goalNode := 40922 // 目的ノードIDを天神駅に設定
// CSVファイルを作成
csvFile, err := os.Create("random_walk_result.csv")
if err != nil {
log.Fatal(err)
}
defer csvFile.Close()
writer := csv.NewWriter(csvFile)
defer writer.Flush()
// CSVのヘッダーを作成
err = writer.Write([]string{"Node", "Latitude", "Longitude"})
if err != nil {
log.Fatal(err)
}
// ランダムウォークの実行
err = performRandomWalk(db, startNode, goalNode, writer)
if err != nil {
log.Fatal(err)
}
fmt.Println("Random walk completed and saved to CSV.")
}
// performRandomWalk ランダムウォークを実行し、結果をCSVに保存
func performRandomWalk(db *sql.DB, startNode, goalNode int, writer *csv.Writer) error {
currentNode := startNode
visited := make(map[int]int) // ノードの訪問回数を追跡
pathStack := []int{} // バックトラッキング用のスタック
for currentNode != goalNode {
// 現在のノードの緯度と経度を取得し、CSVに書き込み
lat, lon, err := getNodeCoordinates(db, currentNode)
if err != nil {
return err
}
fmt.Printf("Node %d: Latitude: %f, Longitude: %f\n", currentNode, lat, lon)
err = writer.Write([]string{
fmt.Sprintf("%d", currentNode),
fmt.Sprintf("%f", lat),
fmt.Sprintf("%f", lon),
})
if err != nil {
return err
}
writer.Flush()
// 現在のノードを訪問済みとして記録
visited[currentNode]++
// 次のノードを取得（目的地に近づくノードを優先）
nextNode, err := getClosestConnectedNode(db, currentNode, goalNode, visited)
if err != nil || visited[nextNode] > 2 { // 訪問回数が多い場合バックトラック
if len(pathStack) > 0 {
// スタックから一つ戻る
fmt.Printf("Backtracking from Node %d\n", currentNode)
currentNode = pathStack[len(pathStack)-1]
pathStack = pathStack[:len(pathStack)-1]
continue
}
return fmt.Errorf("no viable path found")
}
// 次のノードをスタックに追加
pathStack = append(pathStack, currentNode)
currentNode = nextNode
}
// ゴールノードの座標をCSVに保存
lat, lon, err := getNodeCoordinates(db, goalNode)
if err != nil {
return err
}
fmt.Printf("Arrived at Goal Node %d: Latitude: %f, Longitude: %f\n", goalNode, lat, lon)
err = writer.Write([]string{
fmt.Sprintf("%d", goalNode),
fmt.Sprintf("%f", lat),
fmt.Sprintf("%f", lon),
})
if err != nil {
return err
}
writer.Flush()
return nil
}
// getClosestConnectedNode 現在のノードから接続されているノードのうち、目的地に最も近いノードを取得
func getClosestConnectedNode(db *sql.DB, currentNode, goalNode int, visited map[int]int) (int, error) {
// 現在のノードから接続されているノードを取得
query := `
SELECT target
FROM ways
WHERE source = $1
UNION
SELECT source
FROM ways
WHERE target = $1;
`
rows, err := db.Query(query, currentNode)
if err != nil {
return 0, err
}
defer rows.Close()
// 接続ノードのリストを作成
type NodeDistance struct {
ID int
Distance float64
}
var connectedNodes []NodeDistance
for rows.Next() {
var node int
if err := rows.Scan(&node); err != nil {
return 0, err
}
// ノードの座標を取得
lat, lon, err := getNodeCoordinates(db, node)
if err != nil {
return 0, err
}
// 目的地との距離を計算
goalLat, goalLon, err := getNodeCoordinates(db, goalNode)
if err != nil {
return 0, err
}
distance := calculateDistance(lat, lon, goalLat, goalLon)
// 訪問回数が少ないノードを優先
if visited[node] < 3 {
connectedNodes = append(connectedNodes, NodeDistance{
ID: node,
Distance: distance,
})
}
}
if len(connectedNodes) == 0 {
return 0, fmt.Errorf("no connected nodes found")
}
// 目的地に最も近いノードを優先して選択
var closestNode NodeDistance
minDistance := math.MaxFloat64
for _, node := range connectedNodes {
// 最も目的地に近いノードを選択
if node.Distance < minDistance {
closestNode = node
minDistance = node.Distance
}
}
// ランダム選択の要素を追加（袋小路を避けるため）
if len(connectedNodes) > 1 && rand.Float64() < 0.2 {
return connectedNodes[rand.Intn(len(connectedNodes))].ID, nil
}
return closestNode.ID, nil
}
// getNodeCoordinates ノードIDに対応する緯度と経度を取得
func getNodeCoordinates(db *sql.DB, nodeID int) (float64, float64, error) {
var lat, lon float64
query := `
SELECT ST_Y(the_geom) AS lat, ST_X(the_geom) AS lon
FROM ways_vertices_pgr
WHERE id = $1;
`
err := db.QueryRow(query, nodeID).Scan(&lat, &lon)
if err != nil {
return 0, 0, err
}
return lat, lon, nil
}
// calculateDistance 2つの座標間の距離を計算
func calculateDistance(lat1, lon1, lat2, lon2 float64) float64 {
rad := func(deg float64) float64 {
return deg * math.Pi / 180
}
lat1Rad, lon1Rad := rad(lat1), rad(lon1)
lat2Rad, lon2Rad := rad(lat2), rad(lon2)
dlat := lat2Rad - lat1Rad
dlon := lon2Rad - lon1Rad
a := math.Sin(dlat/2)*math.Sin(dlat/2) + math.Cos(lat1Rad)*math.Cos(lat2Rad)*math.Sin(dlon/2)*math.Sin(dlon/2)
c := 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a))
const EarthRadius = 6371e3 // 地球の半径（メートル）
return EarthRadius * c
}
上記のDBはpostGISで構成されており、ways, ways_vertices_pgr があり、SQLでダイクストラ計算が行えます。
ここで以下のようなステップを行うランダムウォークをpostGISを使ったpostgreSQLを使ってGo言語で作成してみました。
(Step.1) ランダムォークの終了地点は、天神駅とする
(Step.2) 出発地点を乱数デ適当に選んで下さい。出発地点は、天神駅から1km程度の地点を適当に抽出する
(Step.3)天神駅の座標から最も近いNode(以下、これを天神ノードといいます)と、出発地点から最も近いNode(以下、これを出発ノードと言います)を抽出する
(Step.4) 出発ノードと天神ノードを直線で繋いで、出発ノードから-90度から+90度の範囲のノードを一つ選ぶ
(Step.5)上記で見つかったノードを次の出発ノードとする。正し、すでに出発ノードとしたノードを、新しい出発ノードにはできないようにする
(Step.6)上記(Step.4)に戻って、天神ノードに到着するまで、出発ノードを更新し続ける。なお更新される出発ノードは全てノド番号と緯度と座標を表示させる
というようなアルゴリズムで、ランダムウォークを作ろうとしたのですが、ダメでした。
一言でいうと「どんつき」の路地に入る→袋小路に入ってしまって、そこで動けなくなるからです。
これを回避するアルゴリズムを考えるの面倒になって、pgr_ksp()を使うことにしました。
pgr_ksp()は、最短経路問題の1、2、3、4番目を出してくれるので、これで1000番目とかを仕えば、そこそこランダムウォークになるかな、と思いまして。
// G:\home\ebata\hakata\src\others\main25\main.go
package main
import (
"database/sql"
"encoding/csv"
"fmt"
"log"
"os"
_ "github.com/lib/pq"
)
const (
connStr = "user=postgres password=password host=127.0.0.1 port=15432 dbname=hakata_db sslmode=disable"
)
type PathResult struct {
PathID int
Node int
Lat float64
Lon float64
}
func main() {
// PostgreSQLデータベースへの接続
db, err := sql.Open("postgres", connStr)
if err != nil {
log.Fatal(err)
}
defer db.Close()
// すべてのパス結果を一度に取得
results, err := getPathResults(db)
if err != nil {
log.Fatal(err)
}
// 各パスIDの結果をCSVに保存
pathIDs := []int{1, 5, 10, 100, 1000, 10000}
for _, pathID := range pathIDs {
filename := fmt.Sprintf("path_id_%d.csv", pathID)
err := saveResultsToCSV(results, pathID, filename)
if err != nil {
log.Fatal(err)
}
}
fmt.Println("CSV files created successfully.")
}
// getPathResults すべてのパス結果を取得
func getPathResults(db *sql.DB) ([]PathResult, error) {
// クエリで指定のpath_idの結果をすべて取得
query := `
WITH ksp_result AS (
SELECT * FROM pgr_ksp(
'SELECT gid AS id, source, target, cost FROM ways',
35382, -- 出発ノードID
40922, -- 目的ノードID
10000, -- 上位10000本の経路を取得
false -- エッジの向きを考慮しない
)
)
SELECT ksp.path_id, ksp.node, ST_Y(wv.the_geom) AS lat, ST_X(wv.the_geom) AS lon
FROM ksp_result AS ksp
JOIN ways_vertices_pgr AS wv ON ksp.node = wv.id
WHERE ksp.path_seq > 0
ORDER BY ksp.path_id, ksp.path_seq;
`
// クエリの実行
rows, err := db.Query(query)
if err != nil {
return nil, err
}
defer rows.Close()
// 結果をメモリに保存
var results []PathResult
for rows.Next() {
var result PathResult
err := rows.Scan(&result.PathID, &result.Node, &result.Lat, &result.Lon)
if err != nil {
return nil, err
}
results = append(results, result)
}
return results, nil
}
// saveResultsToCSV 指定されたpath_idの結果をCSVに保存
func saveResultsToCSV(results []PathResult, pathID int, filename string) error {
// CSVファイルの作成
file, err := os.Create(filename)
if err != nil {
return err
}
defer file.Close()
writer := csv.NewWriter(file)
defer writer.Flush()
// CSVのヘッダーを作成
err = writer.Write([]string{"Node", "Latitude", "Longitude"})
if err != nil {
return err
}
// 結果をCSVに書き込む
for _, result := range results {
if result.PathID == pathID {
err := writer.Write([]string{
fmt.Sprintf("%d", result.Node),
fmt.Sprintf("%f", result.Lat),
fmt.Sprintf("%f", result.Lon),
})
if err != nil {
return err
}
}
}
return nil
}
1番目、10番目、1000番目を表示してみました。
まあ、ランダムウォークという感じではないですが、私たちが移動するときも、まあ、こんな感じだと思うので、とりあえずこれで進めようかなと思います。
<!-- Way: Dummy Railway Line for JR篠栗線 -->
<way id="123456800" version="1" timestamp="2024-08-22T06:13:26Z" changeset="99999999" uid="999999" user="kobore.net">
<nd ref="1882462094" /> <!-- 博多駅 -->
<nd ref="1239476659" /> <!-- 吉塚駅 -->
<nd ref="1843493667" /> <!-- 柚須駅 -->
<nd ref="1843644210" /> <!-- 伊賀駅 -->
<nd ref="1842266870" /> <!-- 原町駅 -->
<nd ref="5445230442" /> <!-- 長者原駅 -->
<nd ref="6072021673" /> <!-- 門松駅 -->
<nd ref="796761953" /> <!-- 篠栗駅 -->
<nd ref="796761953" /> <!-- 筑前山手駅 -->
<nd ref="796761953" /> <!-- 城戸南蔵院前駅 -->
<tag k="railway" v="dummyrail" />
<tag k="name" v="Dummy Railway Line JR篠栗線" />
</way>
このデータを、
<relation id="17960683" version="1" timestamp="2024-08-22T06:13:26Z" changeset="155587618" uid="109705" user="gscholz">
<member type="way" ref="1309534187" role="outer"/>
<member type="way" ref="1309192923" role="inner"/>
<member type="way" ref="1309192922" role="inner"/>
<member type="way" ref="1309192926" role="inner"/>
<member type="way" ref="1309192929" role="inner"/>
<member type="way" ref="1309192924" role="inner"/>
<member type="way" ref="1309192925" role="inner"/>
<member type="way" ref="1309192928" role="inner"/>
<member type="way" ref="733239311" role="inner"/>
<member type="way" ref="1309192927" role="inner"/>
<tag k="highway" v="pedestrian"/>
<tag k="surface" v="gravel"/>
<tag k="type" v="multipolygon"/>
</relation>
この辺に置いた
</osm>
こうしたら、
G:\home\ebata\hakata\hakata_db>osm2pgrouting -f hakata_ex.osm -c mapconfig_for_ex.xml -d hakata_ex_db -U postgres -h localhost -p 15432 -W password
でDBの作成に失敗しました。
wayのタグはwayのタグの集まっているところに置いておかないと、DBの作成で失敗するようです。
(Step 1)地図をスキャンする
(Step 2)SAIで地図を切り取って、さらに追加したレイヤに点を打つ
(Step 3)点を打ったレイヤだけをjpgでセーブする
page61.jpgの内容
(Step 4)このjpegを左上の座標と右上の座標を指定して、ChatGPTに添付して抽出をお願いする。
添付したファイルは地図のある建物の位置を示しているものです。
左上の経度、緯度は、130.41636503641814, 33.59644753614503
右下の経度、緯度は、130.42569034464057, 33.58606144150355
となっています。
添付したファイルに記載された点の座標の経度、緯度を算出して下さい。
(Step 5)で、色々失敗するので、何度か修正をお願いする。たとえば以下のような感じ。
左上の経度、緯度は、130.41636503641814, 33.59644753614503
右下の経度、緯度は、130.42569034464057, 33.58606144150355
ファイルには、赤色、緑色の点が記載されています。この点の座標を、上記の左上の緯度経度と右下の移動経度から、線形計算をして算出して下さい。
私はこの添付ファイルに20～30の点しか付けていませんが、頂いた点は非常に多いです。見直して下さい。
同じ点が2～3個で表示されているようです。点と点の距離が20メートル以内であれば、点を1つにして下さい。
(Step 6)ChatGPTが作ってくれたcsvファイルをQGIS等で表示して確認する。
まあ、とりあえず、位置情報だけは、ChatGPTの協力を得て、取り出すことができました。
しかし、店舗名の紐づけは、どうしようもありませんね。これを何とかできる人は、私を助けて下さい。
======
と思っていたら、以外なほど簡単にできました。
QGISで上記のCSVファイルを示すことができましたが、csvファイルとの対応が分かりません。
例えば、
#,Latitude,Longitude
1,33.59181234, 130.4199366
2.33.59123927, 130.4190329
3.33.59106178 130.4217193
のようにして、QGISで表示される点に#の番号を表示させるにはどうしたら良いですか？
とChatGPTに聞いたら教えてくれました。
まず、普通にcsvを取り込む
点が見えにくいので、大きな赤丸にする。
csvのレイヤを選択して → プロパティ → シンポロジで "dot red" を選ぶ(趣味の問題だが)
で、再び、
csvのレイヤを選択して → プロパティ → ラベル → (上のほうに"なし"と記載されているメニューから) 「単一定義」を選択。あとはよく分からんが、適用を押す。
おお！ ちゃんと番号が付与されている。
これで、対応付けできるし、なんなら地物名称も入れれるかもしれない。
Teamsの会議でビデオの音声を共有するには、以下の手順を実行してください：
この設定により、Teams会議の参加者全員があなたのコンピューターの音声を聞くことができるようになります。
「Redisをインストールすることなく、Windows10で使う方法」がある、と、ChatGPTに言われて「ホンマかいな？」と思いながらやってみたのですが、さっくり動きました。
$ docker run --name redis -p 6379:6379 -d redis
これだけで、Windows上で動かしているGoプログラムの中のredisが動いています(理屈は分かりません)。
Redisの動作確認
Dockerコンテナ内のRedisが動作しているか確認するためには、以下のコマンドを実行します。
$ docker exec -it redis redis-cli ping
-----
また、「docker run --name redis -p 6379:6379 -d redis を次回再起動するには、どうしたら良いですか」と聞いてみたところ、
docker start redis
docker psコマンドを使って、Redisコンテナが正しく動作しているか確認できます：
docker ps
Redisコンテナが実行中であれば、一覧に表示されます。
だ、そうです。ちなみに、
もし、Dockerが起動するたびにRedisコンテナを自動的に再起動させたい場合は、--restartオプションを使って自動再起動ポリシーを設定できます。次回docker runを実行するとき、--restartオプションを追加しておくと、自動再起動が設定されます。
$ docker run --name redis -p 6379:6379 -d --restart unless-stopped redis
私は、redisが必要のない時に立ち上げたくはないので、上記の設定はしていません。
以上