From 28a19648e878d7c0784a0e2fc2daceaf65d71540 Mon Sep 17 00:00:00 2001
From: Leo Strobel <leo.strobel@uni-wuerzburg.de>
Date: Fri, 1 Mar 2024 10:42:55 +0100
Subject: [PATCH] [BUGFIX] Solved an issue with OD-Matrix calibration when
 transitions missing in the matrix where assumed to be zero in the
 uncalibrated model

---
 .../de/uniwuerzburg/omod/core/ODZone.kt       |  3 +-
 .../kotlin/de/uniwuerzburg/omod/core/Omod.kt  | 48 +++++++++++--------
 2 files changed, 31 insertions(+), 20 deletions(-)

diff --git a/src/main/kotlin/de/uniwuerzburg/omod/core/ODZone.kt b/src/main/kotlin/de/uniwuerzburg/omod/core/ODZone.kt
index 21af279f..bc3278f3 100644
--- a/src/main/kotlin/de/uniwuerzburg/omod/core/ODZone.kt
+++ b/src/main/kotlin/de/uniwuerzburg/omod/core/ODZone.kt
@@ -62,7 +62,8 @@ data class ODZone (
             for (entry in geoJson.features) {
                 val properties = entry.properties as GeoJsonODProperties
                 val originZone = nameMapping[properties.origin]!!
-                originZone.destinations = properties.destinations.map { Pair(nameMapping[it.key]!!, it.value) }
+                originZone.destinations = properties.destinations
+                    .filter { it.value > 0 }.map {Pair(nameMapping[it.key]!!, it.value)}
             }
 
             // Check if OD is valid
diff --git a/src/main/kotlin/de/uniwuerzburg/omod/core/Omod.kt b/src/main/kotlin/de/uniwuerzburg/omod/core/Omod.kt
index 7e8e7659..1049b74f 100644
--- a/src/main/kotlin/de/uniwuerzburg/omod/core/Omod.kt
+++ b/src/main/kotlin/de/uniwuerzburg/omod/core/Omod.kt
@@ -242,7 +242,11 @@ class Omod(
             // Calculate omod origin probability. For speed only on zone level.
             omodWeights[odZone] = odZone.aggLocs.sumOf { omodProbs[it]!! }
             // Calculate OD-Matrix origin probability.
-            odWeights[odZone] = odZone.destinations.filter { it.first.inFocusArea }.sumOf { it.second }
+            odWeights[odZone] = if (odZone.inFocusArea) {
+                odZone.destinations.sumOf { it.second }
+            } else {
+                odZone.destinations.filter { it.first.inFocusArea }.sumOf { it.second }
+            }
         }
 
         val weightSumOMOD = omodWeights.values.sum()
@@ -285,25 +289,31 @@ class Omod(
         require(activities == Pair(ActivityType.HOME, ActivityType.WORK)) {
             "Only OD-Matrices with Activities HOME->WORK are currently supported"
         }
-        
+
         val factors = mutableMapOf<Pair<ODZone, ODZone>, Double>()
         val priorProbs = calcOMODProbsAsMap(activities.first)
 
         for (originOdZone in odZones) {
-            val omodWeights = mutableMapOf<ODZone, Double>()
-            val odWeights = mutableMapOf<ODZone, Double>()
+            val omodWeights = odZones.associateWith { 0.0 } as MutableMap<ODZone, Double>
+            val odWeights = odZones.associateWith { 0.0 } as MutableMap<ODZone, Double>
+
+            // Calculate omod transition probability. For speed only on zone level.
+            for (origin in originOdZone.aggLocs) {
+                val pPriorLoc = priorProbs[origin]!!
+                if (pPriorLoc == 0.0) { continue }
+                val wDependentZone = getWeights(origin, zones, activities.second).sum()
+
+                for (destOdZone in odZones) {
+                    val wDependentLoc = getWeights(origin, destOdZone.aggLocs, activities.second).sum()
+                    omodWeights[destOdZone] = omodWeights[destOdZone]!! + (pPriorLoc * wDependentLoc / wDependentZone)
+                }
+            }
+
+            // Calculate OD-Matrix origin probability.
             for ((destOdZone, transitions) in originOdZone.destinations) {
-                // Calculate omod transition probability. For speed only on zone level.
-                var omodWeight = 0.0
-                for (origin in originOdZone.aggLocs) {
-                    val prior = priorProbs[origin]!!
-                    if (prior == 0.0) { continue }
-                    omodWeight +=  prior * getWeights(origin, destOdZone.aggLocs, activities.second).sum()
-                 }
-                omodWeights[destOdZone] = omodWeight
-
-                // Calculate OD-Matrix origin probability.
-                odWeights[destOdZone] = if (destOdZone.inFocusArea || originOdZone.inFocusArea) transitions else 0.0
+                if (destOdZone.inFocusArea || originOdZone.inFocusArea) {
+                    odWeights[destOdZone] = transitions
+                }
             }
 
             val weightSumOMOD = omodWeights.values.sum()
@@ -315,7 +325,7 @@ class Omod(
                     factors[Pair(originOdZone, destOdZone)] = 1.0
                 }
             } else {
-                for (destOdZone in originOdZone.destinations.map { it.first }) {
+                for (destOdZone in odZones) {
                     // Normalize
                     val omodProb = omodWeights[destOdZone]!! / weightSumOMOD
                     val odProb = odWeights[destOdZone]!! / weightSumOD
@@ -806,11 +816,11 @@ class Omod(
 
         // Work distribution
         val workProbs = DoubleArray(zones.size) { 0.0 }
-        for (zone in zones) {
+        for ((i, zone) in zones.withIndex()) {
             val workWeights = getWeights(zone, zones, ActivityType.WORK)
             val totalWorkWeight = workWeights.sum()
-            for (i in zones.indices) {
-                workProbs[i] += homeProbs[i] * workWeights[i] / totalWorkWeight
+            for (j in zones.indices) {
+                workProbs[j] += homeProbs[i] * workWeights[j] / totalWorkWeight
             }
         }
         return workProbs