macromultimodelling:fix_real_data_and_add_ecoWeight_param(diff_from_ecoMult_in_fact)

Author: JusteRaimbault

library/src/main/scala/org/openmole/spatialdata/application/multimodelmacro/RealData.scala

@@ -1,7 +1,8 @@
 package org.openmole.spatialdata.application.multimodelmacro
 
 import org.openmole.spatialdata.model.urbandynamics.{Coevolution, EconomicExchanges, Innovation, MacroResultFit, MultiMacroModel}
-import org.openmole.spatialdata.model.urbandynamics.MultiMacroModel.MultiMacroResult
+import org.openmole.spatialdata.model.urbandynamics.MultiMacroModel.{MultiMacroResult, MultiMacroState}
+import org.openmole.spatialdata.utils
 import org.openmole.spatialdata.utils.io.CSV
 import org.openmole.spatialdata.utils.math.{DenseMatrix, EmptyMatrix, Matrix, Statistics}
 import org.openmole.spatialdata.utils.math.Matrix.MatrixImplementation
@@ -21,14 +22,15 @@ object RealData {
                      result: MultiMacroResult,
                      fit: MacroResultFit
                    ){
-    def logmse: Double = 0.0
-    def mselog: Double = 0.0
+    def logmse: Double = fit.logmse
+    def mselog: Double = fit.mselog
   }
 
   object Result {
     def apply(result: MultiMacroResult): Result = {
-      val targetPop = result.states.head.populations
-      val simulatedStates = result.states.toVector
+      val targetPop = result.realPopulations
+      val simulatedStates: Vector[MultiMacroState] = result.states.toVector
+      utils.log(s"Computing real fit: total target pops ${targetPop.colSum.toSeq}; simulated pops: ${simulatedStates.map(_.populations.sum)}")
       val fit = MacroResultFit(targetPop, simulatedStates)
       Result(result, fit)
     }
@@ -48,6 +50,8 @@ object RealData {
     val rawdates: Seq[String] = CSV.readCSV(datesFile, withHeader = false).values.toSeq.head
     val dates: Array[Double] = rawdates.map(_.toDouble).toArray
 
+    utils.log(s"Real data setup : total populations ${populationMatrix.colSum.toSeq} ; average distance ${distancesMatrix.flatValues.sum/(distancesMatrix.nrows*distancesMatrix.nrows)}")
+
     (populationMatrix, distancesMatrix, dates)
   }
 

library/src/main/scala/org/openmole/spatialdata/application/multimodelmacro/RunSDG.scala

@@ -12,51 +12,55 @@ object RunSDG extends App {
 
   implicit val rng: Random = new Random
 
-  // synthetic setup
-  /*val initialState = SDG.setupSynthetic(syntheticCities = 30,
-    syntheticHierarchy = 1.0,
-    syntheticMaxPop =  100000.0,
-    finalTime = 20
-  )
-  */
+  val synthetic = false
 
-  // real data setup
-  val initialState = RealData.setupRealData(
-    "","","" //TODO
-  )
+  val initialState =
+    if(synthetic) {
+      // synthetic setup
+      SDG.setupSynthetic(syntheticCities = 30,
+        syntheticHierarchy = 1.0,
+        syntheticMaxPop = 100000.0,
+        finalTime = 20, // limit to avoid huge values: 50 -> up to larger than 1e9 !
+        seed = 42L
+      )
+    }else {
+      // real data setup
+      RealData.setupRealData("data/urbansystems/EU_pops.csv", "data/urbansystems/EU_dist.csv", "data/urbansystems/EU_dates.csv")
+    }
 
-  val result = SDG.runMultiModelMacro(
+  val result: SDG.Result = SDG.runMultiModelMacro(
     SDG.setupMultiModelMacro(
       initialState._1, initialState._2, initialState._3,
-      finalTime = 20, // limit to avoid huge values: 50 -> up to larger than 1e9 !
+      finalTime = initialState._3.length,
       seed = rng.nextInt(),
-      innovationWeight = 0.005,
-      innovationGravityDecay = 1.0,//0.5,
-      innovationInnovationDecay = 0.3,
+      innovationWeight = 1e-2, // 0.005
+      innovationGravityDecay = 1000.0,//0.5,
+      innovationInnovationDecay = 1000.0,
       innovationMutationRate = 0.2,
       innovationNewInnovationHierarchy = 0.5,
       innovationEarlyAdoptersRate = 0.2,
       innovationUtilityStd = 1.0,
       innovationUtilityDistrib = "log-normal",
-      ecoWeight = 0.05,
+      ecoWeight = 1e-2,
+      ecoMultiplier = 2.0,
       ecoSizeEffectOnDemand = 0.1,
       ecoSizeEffectOnSupply = 0.2,
-      ecoGravityDecay = 0.1,//0.5,
-      ecoWealthToPopulationExponent = 1.5,
-      ecoPopulationToWealthExponent = 1.5,
-      coevolWeight = 0.005,
+      ecoGravityDecay = 1000.0,//0.5,
+      ecoWealthToPopulationExponent = 1.0,
+      ecoPopulationToWealthExponent = 1.0,
+      coevolWeight = 1e-2,
       coevolGamma = 1.0,
-      coevolGravityDecay = 0.5,
-      coevolNetworkGmax = 0.05,
+      coevolGravityDecay = 1000.0,
+      coevolNetworkGmax = 0.0001,
       coevolNetworkExponent = 1.0,
       coevolNetworkThresholdQuantile = 0.5
     )
   )
 
   utils.log(s"Indicators: Emissions = ${SDG.cumulatedFlows(result)}; Innovation = ${- SDG.averageUtility(result)}; Infrastructure = ${SDG.averageDistance(result)}; Eco inequality = ${SDG.giniEconomicWealth(result)}; Wealth = ${SDG.averageWealth(result)}")
 
-  val realresult: RealData.Result = RealData.Result(result.asInstanceOf[MultiMacroResult])
+  val realresult: RealData.Result = RealData.Result.apply(result.result)
 
-  utils.log(s"Fit on populations : logmse =  ${realresult.logmse} ; mselog = ${realresult.mselog}")
+  utils.log(s"Fit on populations (real setup): logmse =  ${realresult.logmse} ; mselog = ${realresult.mselog}")
 
 }

library/src/main/scala/org/openmole/spatialdata/application/multimodelmacro/SDG.scala

@@ -4,7 +4,7 @@ import org.openmole.spatialdata.model.urbandynamics.Coevolution.CoevolutionState
 import org.openmole.spatialdata.model.urbandynamics.EconomicExchanges.EconomicExchangesState
 import org.openmole.spatialdata.model.urbandynamics.Innovation.{InnovationState, InnovationUtilityLogNormalDistribution, InnovationUtilityNormalDistribution, mutationInnovation}
 import org.openmole.spatialdata.model.urbandynamics.MultiMacroModel.MultiMacroResult
-import org.openmole.spatialdata.model.urbandynamics.{Coevolution, EconomicExchanges, Innovation, MultiMacroModel}
+import org.openmole.spatialdata.model.urbandynamics.{Coevolution, EconomicExchanges, Innovation, MacroResult, MultiMacroModel}
 import org.openmole.spatialdata.utils.math.Matrix.MatrixImplementation
 import org.openmole.spatialdata.utils.math.{DenseMatrix, EmptyMatrix, Matrix, Statistics}
 import org.openmole.spatialdata.vector.measures.Spatstat
@@ -77,6 +77,7 @@ object SDG {
                            innovationUtilityStd: Double,
                            innovationUtilityDistrib: String,
                            ecoWeight: Double,
+                           ecoMultiplier: Double,
                            ecoSizeEffectOnDemand: Double,
                            ecoSizeEffectOnSupply: Double,
                            ecoGravityDecay: Double,
@@ -102,7 +103,7 @@ object SDG {
     // eco exchanges - ! add gravitypotentials in eco state (total emission: cumulated interactions across models)
     // ! issue: not same unit for innov / eco (supply - demand gravity potentials) -> should compute something equivalent?
     // ! marius weight? - rescale pop updates?
-    val ecoModel = EconomicExchanges(population, distances, dates, ecoWeight, ecoSizeEffectOnDemand, ecoSizeEffectOnSupply, ecoGravityDecay, ecoWealthToPopulationExponent, ecoPopulationToWealthExponent)
+    val ecoModel = EconomicExchanges(population, distances, dates, ecoWeight, ecoMultiplier, ecoSizeEffectOnDemand, ecoSizeEffectOnSupply, ecoGravityDecay, ecoWealthToPopulationExponent, ecoPopulationToWealthExponent)
     val ecoInitialState = EconomicExchanges.initialState(ecoModel)
 
     // coevolution
@@ -111,10 +112,10 @@ object SDG {
       0.0, 1.0, 1.0, coevolNetworkGmax, coevolNetworkExponent, coevolNetworkThresholdQuantile)
     val coevolInitialState = Coevolution.initialState(coevolModel)
 
-    MultiMacroModel(Seq(innovModel, ecoModel, coevolModel), Seq(innovInitialState, ecoInitialState, coevolInitialState))
+    MultiMacroModel(population.clone,Seq(innovModel, ecoModel, coevolModel), Seq(innovInitialState, ecoInitialState, coevolInitialState))
   }
 
-  def runMultiModelMacro(model: MultiMacroModel)(implicit rng: Random): Result = Result(model.run.asInstanceOf[MultiMacroResult])
+  def runMultiModelMacro(model: MultiMacroModel)(implicit rng: Random): Result = Result(model.run)
 
 
   /**
@@ -151,7 +152,7 @@ object SDG {
     */
   def averageDistance(res: Result): Double = res.result.states.map{s =>
     val ptot = s.populations.sum
-    val diagpop = DenseMatrix.diagonal(s.populations.flatValues.map( _ / ptot))
+    val diagpop = DenseMatrix.diagonal(s.populations.map( _ / ptot).toArray)
     (diagpop %*% s.distanceMatrix %*% diagpop).sum
   }.sum / res.result.states.length
 

library/src/main/scala/org/openmole/spatialdata/model/urbandynamics/Coevolution.scala

@@ -49,7 +49,7 @@ case class Coevolution(
 
   override def run: MacroResultFit = Coevolution.run(model = this)
 
-  override def nextStep(state: MacroState, populations: Matrix, distanceMatrix: Matrix): MacroState = {
+  override def nextStep(state: MacroState, populations: Seq[Double], distanceMatrix: Matrix): MacroState = {
     val updatedState = Coevolution.updateState(state.asInstanceOf[CoevolutionState], populations, distanceMatrix)
     Coevolution.nextState(this, updatedState)
   }
@@ -71,7 +71,7 @@ object Coevolution {
 
   case class CoevolutionState(
                              time: Int,
-                             populations: Matrix,
+                             populations: Seq[Double],
                              distanceMatrix: Matrix,
                              interactionPotentials: Matrix,
                              flows: Matrix
@@ -155,7 +155,7 @@ object Coevolution {
     )
   }
 
-  def updateState(state: CoevolutionState, populations: Matrix, distanceMatrix: Matrix): CoevolutionState = {
+  def updateState(state: CoevolutionState, populations: Seq[Double], distanceMatrix: Matrix): CoevolutionState = {
     state.copy(populations = populations, distanceMatrix = distanceMatrix)
   }
 
@@ -168,7 +168,7 @@ object Coevolution {
     val interactionPotentials = computeFlows(pops, gravityDistanceWeights, model.gravityGamma, {case (p,g) =>math.pow(p / totalpop,g) })
     val flows =  computeFlows(pops, gravityDistanceWeights, 1.0, {case (p,_) =>p / totalpop })
 
-    CoevolutionState(0, pops, dmat, interactionPotentials, flows)
+    CoevolutionState(0, pops.flatValues.toSeq, dmat, interactionPotentials, flows)
   }
 
   /**
@@ -188,12 +188,12 @@ object Coevolution {
     val gravityDistanceWeights = updatedPopulationState.distanceMatrix.map{ d => Math.exp(-d / model.gravityDecay) }
 
     val totalpop = updatedPopulationState.populations.sum
-    val interactionPotentials = computeFlows(updatedPopulationState.populations, gravityDistanceWeights, model.gravityGamma, {case (p,g) =>math.pow(p / totalpop,g) })
+    val interactionPotentials = computeFlows(RealMatrix(updatedPopulationState.populations.toArray,row=false), gravityDistanceWeights, model.gravityGamma, {case (p,g) =>math.pow(p / totalpop,g) })
 
     val updatedNetwork = updateNetwork(updatedPopulationState.distanceMatrix, interactionPotentials, networkGmax, networkExponent, networkThresholdQuantile)
 
     val updatedGravityDistanceWeight = updatedNetwork.map{ d => Math.exp(-d / model.gravityDecay) }
-    val flows = computeFlows(updatedPopulationState.populations, updatedGravityDistanceWeight, 1.0, {case (p,_) => p / totalpop})
+    val flows = computeFlows(RealMatrix(updatedPopulationState.populations.toArray,row=false), updatedGravityDistanceWeight, 1.0, {case (p,_) => p / totalpop})
 
     utils.log(s"Delta distance = ${updatedPopulationState.distanceMatrix.flatValues.zip(updatedNetwork.flatValues).map{case (d1,d2) =>math.abs(d1-d2)}.sum}")
 
@@ -244,7 +244,7 @@ object Coevolution {
     def step(state: (Vector[MacroStateGen],Vector[Double])): (Vector[MacroStateGen],Vector[Double]) =
       (Vector(MacroModel.macroStep(state._1.head,growthRates(model),state._2.head))++state._1,state._2.tail)
 
-    val initState = MacroStateGen(0,populationMatrix.getSubmat(0, 0, nrows = n, ncols= 1), distancesMatrices(0))
+    val initState = MacroStateGen(0,populationMatrix.getSubmat(0, 0, nrows = n, ncols= 1).flatValues.toSeq, distancesMatrices(0))
     val deltats = dates.tail.zip(dates.dropRight(1)).map{case (next,prev)=> next-prev}.toVector
 
     val finalState = Iterator.iterate((Vector(initState),deltats))(step).takeWhile(_._2.nonEmpty).toVector.last._1