江端さんの技術メモ

第三章を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も大きくなったということになります

======
尤度をもっとも簡単に説明する方法

つまるところ、偏ったコインを作る = パラメータをつけくわえる、という理解で良いようです。

=====

ところで、今回作成して頂いたモデルでマルチエージェントシミュレーションのエージェントに適用させたい場合、何を入力値として、エージェントの交通選択を選ぶことになるか？

"所要時間"と"費用"を入力値とすれば良い。

//  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の会議でビデオの音声を共有するには、以下の手順を実行してください：

1. 会議を開始するか参加：
   • Teamsの会議を開始する、もしくは既存の会議に参加してください。
2. 「画面を共有」ボタンを押す：
   • 会議ウィンドウの中央下部にある「共有」ボタン（四角と上向きの矢印アイコン）をクリックします。
3. 共有するウィンドウや画面を選択：
   • 共有したいビデオを再生しているウィンドウや画面、もしくはデスクトップ全体を選択します。
4. 「コンピューターの音声を含める」オプションを有効化：
   • 画面共有の際に、画面の上部に「コンピューターの音声を含める」というチェックボックスが表示されます。これにチェックを入れることで、ビデオの音声も一緒に共有できます。

1. ビデオを再生：
   • 画面共有を開始した後、会議参加者にビデオの映像と音声が共有されます。

この設定により、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 を次回再起動するには、どうしたら良いですか」と聞いてみたところ、

Redisコンテナを再起動するコマンド

1. コンテナを再起動する コンテナが停止している場合、次のコマンドで再起動します：
   docker start redis

2. コンテナが実行中かどうかを確認する 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が必要のない時に立ち上げたくはないので、上記の設定はしていません。

以上

(1)地図DBに接続されている状態であること。

(2)go.modが入っていること

(3)prune_isonodes.go の依存パッケージのインストール方法が変っていること

この地図のプロパティは、424 x 379

で、

set datafile separator "," # データの区切り文字をコンマに指定
set title "Cluster Data with Map" # グラフタイトル
set xlabel "X-axis" # X軸のラベル
set ylabel "Y-axis" # Y軸のラベル
set zlabel "Z-axis" # Z軸のラベル
set xrange [0:1] # X軸の範囲をデータに合わせて調整
set yrange [0:1] # Y軸の範囲をデータに合わせて調整
set zrange [0:1] # Z軸の範囲
set ticslevel 0 # Z軸を平面上に表示
set xyplane at 0 # XY平面をZ=0に設定
set view 40, 30, 1, 1 # 視点の角度とスケールを調整

# 地図画像をXY平面に表示 (画像のX軸、Y軸範囲を設定し、Z軸は0に固定)
splot "img_664a03bcd2376.png" binary filetype=png dx=1.0/424 dy=1.0/379 with rgbimage, \
'Cluster1.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" title 'Cluster 1', \
'Cluster2.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster3.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster4.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster5.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster6.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster7.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster8.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster9.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster10.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster11.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster12.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster13.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster14.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster15.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster16.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster17.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster18.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster19.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster20.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster21.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster22.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster23.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster24.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster25.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster26.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster27.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster28.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster29.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster30.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster31.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster32.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster33.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster34.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster35.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster36.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster37.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster38.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster39.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster40.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster41.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle

cluster1.csvの内容はこんな感じ
0.63, 0.65, 0.58
0.62, 0.65, 0.56
0.61, 0.65, 0.56
0.61, 0.65, 0.56
0.62, 0.65, 0.56
0.62, 0.66, 0.55
0.63, 0.65, 0.55
0.62, 0.65, 0.56
0.62, 0.66, 0.57
0.62, 0.65, 0.58
0.62, 0.65, 0.58
(後略)

chartjs-plugin-streaming.js は、今も便利に使わせて貰っています。

Golangから送り込んだデータを、chartjs-plugin-streaming(リアルタイムストリーミングデータ向け Chart.js プラグイン)で表示してみた件

で、今回、リアルタイムは必要なくて、単にグラフ表示が必要になりました。

このライブラリが必要です。https://www.papaparse.com/ からダウンロードできます。

sample.htmlという名前で以下のファイルをセーブします。

<!DOCTYPE html>
<html lang="en">
<head>
  <meta charset="UTF-8">
  <meta name="viewport" content="width=device-width, initial-scale=1.0">
  <title>CSV to Chart.js</title>
  <script src="./chart.js"></script>
  <script src="./papaparse.min.js"></script> <!-- 正しいPapaParseのCDN -->
</head>
<body>
  <canvas id="myChart" width="400" height="200"></canvas>
  <script>
    // サンプルのCSVファイルを読み込む処理
    fetch('sample.csv')
      .then(response => response.text())
      .then(csvData => {
        Papa.parse(csvData, {
          header: true,
          dynamicTyping: true,
          complete: function(results) {
            const data = results.data;
            const labels = [];
            const values = [];

            // CSVファイルからラベルとデータを抽出
            data.forEach(row => {
              labels.push(row['label_column']);  // 'label_column' はラベル用の列名
              values.push(row['value_column']);  // 'value_column' は値用の列名
            });

            // Chart.jsでグラフを描画
            const ctx = document.getElementById('myChart').getContext('2d');
            const chart = new Chart(ctx, {
              type: 'line',  // 折れ線グラフの例
              data: {
                labels: labels,
                datasets: [{
                  label: 'CSV Data',
                  data: values,
                  borderColor: 'rgba(75, 192, 192, 1)',
                  borderWidth: 1,
                  fill: false
                }]
              },
              options: {
                scales: {
                  x: {
                    beginAtZero: true
                  },
                  y: {
                    beginAtZero: true
                  }
                }
              }
            });
          }
        });
      });
  </script>
</body>
</html>

sample.csv は以下の通り。

label_column,value_column
January,10
February,15
March,20
April,25
May,30
June,35
July,40
August,45
September,50
October,55
November,60
December,65

あとは、

python3 -m http.server 8000

などとしてサーバを立ち上げておいて、ブラウザから、

http://localhost:8000/sample.html

とすれば、