New MCMC script

scenariosAleatoires/MCMC.py

+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+import pandas as pd
+import geopandas as gpd
+import os, sys
+import yaml
+import numpy as np
+from proximite import Proximite
+from resilience_list import Resilience
+from productivite import Productivity
+from indice_biodiversite_2 import Biodiversity
+from social import Social
+from tqdm import tqdm
+import hashlib
+
+class Scenario:
+    def __init__(self, patches):
+        self.patches = patches
+
+    def load_shp(shpfilename):
+        '''
+        Return an instance of class Patches by loading a shapefile of initial patches.
+        '''
+        scenario = Scenario(gpd.GeoDataFrame.from_file(shpfilename, encoding='utf-8'))
+        scenario.patches['init_cult'] = scenario.patches['cultgeopat']
+        return scenario
+
+    def reallocate(self, rng, targetPAT, ratioNbPatches):
+        nbPatches = int(len(self.patches)*ratioNbPatches)
+        surfDelta = targetPAT - self.patches.groupby('cultgeopat')['SURF_PARC'].sum()
+        cult_to_decrease = surfDelta[surfDelta<0].sort_values(ascending=True).keys().tolist()
+        cult_to_increase = surfDelta[surfDelta>0].sort_values(ascending=False).keys().tolist()
+        # Sampling the patches to reallocate
+        samples = self.patches[self.patches['cultgeopat'].isin(cult_to_decrease)].sample(n=nbPatches, random_state=rng)#.reset_index(drop=True)
+        # Building the new culture reallocated
+        factors = surfDelta[cult_to_increase]
+        factors = (factors*len(samples)/factors.sum()).map(round) # normalize on nb samples
+        newCult = pd.Series(cult_to_increase).repeat(factors)
+        if len(newCult) < len(samples): # may be due to factors rounding
+            newCult = newCult.append(newCult.sample(n=len(samples)-len(newCult), random_state=rng), ignore_index=True)
+        newCult = newCult.sample(frac=1, random_state=rng)[:len(samples)].reset_index(drop=True) # shuffle and cut extra elements
+        # Doing the reallocation
+        self.patches.loc[samples.index.values,'cultgeopat'] = newCult.values
+
+class Indicators:
+    def __init__(self, config, initial_patches, patches_md5sum, targetPAT):
+        self._proximity = Proximite(initial_patches, targetPAT['Fruits et légumes'])
+        self._resilience = Resilience(config['resilience'], initial_patches)
+        self._productivity = Productivity()
+        biodiv_config = config['biodiversity']
+        matrixfilename = biodiv_config['matrixfilename']
+        self._biodiversity = Biodiversity(initial_patches, patches_md5sum, biodiv_config['dmax'], biodiv_config['epsilon'])
+        if not os.path.isfile(matrixfilename):
+        	self._biodiversity.compute_dispersal_probabilities(initial_patches)
+        	self._biodiversity.save_probabilities(matrixfilename)
+        self._biodiversity.load_probabilities(matrixfilename)
+        social_config = config['social']
+        self._social = Social(initial_patches, patches_md5sum, social_config['cost_matrix_filename'],
+            md5sum(social_config['cost_matrix_filename']), social_config['bdv_threshold'])
+        patches_costs_filename = social_config['patches_costs_filename']
+        if not os.path.isfile(patches_costs_filename):
+            self._social.compute_patches_transition_cost(initial_patches)
+            self._social.save_patches_transition_cost(patches_costs_filename)
+        self._social.load_patches_transition_cost(patches_costs_filename)
+        self.indicators_names = ['proximity', 'resilience', 'productivity', 'biodiversity', 'social']
+
+    def compute_indicators(self, patches):
+        return [self.proximity(patches), self.resilience(patches),
+                self.productivity(patches), self.biodiversity(patches), self.social(patches)]
+
+    def compute_indicators_pool(self, scenarios):
+        rows=[]
+        for patches in scenarios:
+            rows.append(self.compute_indicators(patches))
+        return pd.DataFrame(rows, columns=self.indicators_names)
+
+    def proximity(self, patches):
+        return self._proximity.compute_indicator(patches, False)
+
+    def resilience(self, patches):
+	       return self._resilience.compute_indicator(patches, False)
+
+    def productivity(self, patches):
+	       return self._productivity.compute_indicator(patches)
+
+    def biodiversity(self, patches):
+	       return self._biodiversity.compute_indicator(patches)
+
+    def social(self, patches):
+        return self._social.compute_indicator(patches)
+
+def md5sum(filename):
+    hash_md5 = hashlib.md5()
+    with open(filename, "rb") as f:
+        for chunk in iter(lambda: f.read(4096), b""):
+            hash_md5.update(chunk)
+    return hash_md5.hexdigest()
+
+class MCMC:
+    def __init__(self, mcmc_config_filename):
+        if not os.path.isfile(mcmc_config_filename):
+        	print('Error: file not found "{}"'.format(mcmc_config_filename))
+        	print('Please copy the template file "MCMC_config.sample.yml" and adjust to your settings and run again this program')
+        	sys.exit(1)
+        self.mcmc_config = yaml.load(open(mcmc_config_filename,'r'))
+        if 'rng_seed' in self.mcmc_config:
+            self.rng = np.random.RandomState(self.mcmc_config['rng_seed'])
+        else:
+            self.rng = np.random.RandomState(42)
+            print('MCMC initialized with default seed') # self.rng.get_state()
+        self.patches_md5sum = md5sum(self.mcmc_config['patches'])
+        self.patches = gpd.GeoDataFrame.from_file(self.mcmc_config['patches'], encoding='utf-8')
+        self.patches = self.patches[self.patches['cultgeopat']!='Non Considérée']
+        self.patches['init_cult'] = self.patches['cultgeopat']
+        self.target = pd.read_csv(self.mcmc_config['target'], sep=';',index_col=0)
+        targetRatio = (self.target['2050']-self.target['2016'])/self.target['2016']
+        self.targetPAT = self.patches.groupby('cultgeopat')['SURF_PARC'].sum()*(1+targetRatio)
+        self.indicators = Indicators(self.mcmc_config['indicators_config'], self.patches, self.patches_md5sum, self.targetPAT)
+
+    def run(self):
+        # Initial sampling and evaluation
+        scores = []
+        for i in tqdm(range(self.mcmc_config['initial_nb_particles'])):
+            scenario = Scenario(self.patches.copy()) # 0.2 ms
+            scenario.reallocate(self.rng, self.targetPAT, self.mcmc_config['ratio_patches_to_modify']) # 3.8 ms
+            scores.append(self.indicators.compute_indicators(scenario.patches))
+        scores = pd.DataFrame(scores, columns=self.indicators.indicators_names)
+        scores.to_csv('../output/mcmc2.csv')
+        # TODO
+        # Storing variation of indicators
+        init_var = scores.std()
+        # Selecting particles
+        # sequential optimization loop
+
+if __name__ == '__main__':
+    # scenario = Scenario.load_shp('../output/PAT_patches/PAT_patches.shp')
+    # target = pd.read_csv('../resources/targetPAT.csv', sep=';',index_col=0)
+    # targetRatio = (target['2050']-target['2016'])/target['2016']
+    # targetPAT = scenario.patches.groupby('cultgeopat')['SURF_PARC'].sum()*(1+targetRatio)
+    # rng = np.random.RandomState()
+    # scenario.reallocate(rng, targetPAT, 50)
+    mcmc = MCMC('MCMC_config.yml')
+    mcmc.run()
+    # print(mcmc.indicators.biodiversity(mcmc.patches))
+    # print(mcmc.indicators.proximity(mcmc.patches))

scenariosAleatoires/MCMC_config.yml

+patches: ../output/PAT_patches/PAT_patches.shp
+target: ../resources/targetPAT.csv
+initial_nb_particles: 10
+ratio_patches_to_modify: 0.05
+indicators_config:
+ proximite:
+   surf_totale_cible: 20165939.605135553
+ resilience: Parcelle_PAT/Grille_resilience.shp
+ biodiversity:
+   dmax: 1000
+   epsilon: 0.001
+   matrixfilename: ../output/matrix_biodiversity.npz
+ social:
+   cost_matrix_filename: ../resources/matrice_transition.csv
+   patches_costs_filename: ../output/patches_transition_costs.npz
+   bdv_threshold: 0.1