new version optimized with data cached

scenariosAleatoires/social.py

 #!/usr/bin/python3
 # -*- coding: utf-8 -*-
 import pandas as pd
+import numpy as np
+import os
 
 class Social:
 	'''
@@ -12,47 +14,97 @@ class Social:
 	These values were defined during reunions and with the review of our experts.
 	'''
 
-	def __init__(self, initial_patches, cost_matrix_filename, bdv_threshold):
+	def __init__(self, initial_patches, patches_md5sum, cost_matrix_filename, cost_matrix_md5sum, bdv_threshold):
 		'''
+		:param initial_patches Initial state of the patches (attributes cultugeopat and SURF_PARC required)
+		:param patches_md5sum md5 checksum of the input file of patches
 		:param cost_matrix_filename input file containing the matrix in CSV format (separator: \t)
 		:param bdv_threshold proportion of a cultural type needed for considering that this cultural type
 		is suffisantly present in the locality for ignoring the conversion cost.
 		'''
+		self.patches_md5sum = patches_md5sum
 		self.matrice_transition = pd.read_csv(cost_matrix_filename, sep='\t',index_col=0)
+		self.cost_matrix_md5sum = cost_matrix_md5sum
 		self.bdv_threshold = bdv_threshold
-		self.initial_patches = initial_patches[['ID_PARCEL','cultgeopat']]
+
+	def compute_patches_transition_cost(self, initial_patches):
+		costs=[]
+		for bdv,grouped in initial_patches.groupby('Bdv'):
+			grouped=grouped.copy()
+			surfaces_initial = grouped.groupby('cultgeopat')['SURF_PARC'].sum()
+			cult_threshold = surfaces_initial > surfaces_initial.sum() * self.bdv_threshold
+		for culture in self.matrice_transition.columns.get_values():
+			for row in initial_patches.itertuples():
+				cost = 0.0
+				# patches that have different farming that others in same Bdv
+				if culture not in cult_threshold or not cult_threshold[culture]:
+					# patches that have different farming that others of the same farmer
+					if culture not in initial_patches[initial_patches['id_expl']==row.id_expl]['cultgeopat']:
+						cost = self.matrice_transition.loc[row.cultgeopat,culture]
+				costs.append([row.ID_PARCEL, row.cultgeopat, culture, cost])
+		self.patches_transition_cost = pd.DataFrame(costs, columns=['ID_PARCEL','initialcult','cultgeopat','cost'])
+
+	def save_patches_transition_cost(self, filename):
+		'''
+		Save the matrix of costs computed for each patch and initial parameters in a file (npz format).
+		'''
+		m = self.patches_transition_cost.values
+		np.savez(filename, bdv_threshold=self.bdv_threshold, patches_md5sum=self.patches_md5sum, cost_matrix_md5sum=self.cost_matrix_md5sum,
+			data=m.data)
+
+	def load_patches_transition_cost(self, filename):
+		'''
+		Load the matrix of costs and initial parameters from a file (npz format).
+
+		The correspondance will be checked between
+		 class attributes and metadata stored in the file
+		'''
+		loader = np.load(filename)
+		error_message = ('Matrix metadata {0} does\'nt correspond to metadata given at indicator construction'
+			'\n\t{0} in metadata: {1} while asked: {2}'
+			'\n\tmatrix filename: ') + filename
+		for k1,k2 in [[self.bdv_threshold,'bdv_threshold'],[self.patches_md5sum,'patches_md5sum'],[self.cost_matrix_md5sum,'cost_matrix_md5sum']]:
+			if k2 not in loader or k1!=loader[k2]:
+				raise ValueError(error_message.format(k2, dict(loader.items()).get(k2, 'NA'), k1))
+		self.patches_transition_cost = pd.DataFrame(loader['data'], columns=['ID_PARCEL','initialcult','cultgeopat','cost'])
 
 	def compute_indicator(self, scenario_patches):
 		'''
-		:param initial_patches: Initial patches
 		:param scenario: The scenario to analyse
 		'''
 		cost = 0
-		filter1 = pd.merge(
-			scenario_patches[['ID_PARCEL','cultgeopat','Bdv','SURF_PARC','id_expl']],
-		 	self.initial_patches,
-			on='ID_PARCEL')
-		for bdv,grouped in filter1.groupby('Bdv'):
-			grouped=grouped.copy()
-			surfaces_initial = grouped.groupby('cultgeopat_y')['SURF_PARC'].sum()
-			cult_threshold = surfaces_initial > surfaces_initial.sum() * self.bdv_threshold
-			filter2 = grouped[grouped['cultgeopat_y'] != grouped['cultgeopat_x']]
-			for i,row in filter2.iterrows():
-				# patches that have different farming that others in same Bdv
-				if cult_threshold[row['cultgeopat_y']]:
-					# patches that have different farming that others of the same farmer
-					if row['cultgeopat_y'] not in filter1[filter1['id_expl']==row['id_expl']]['cultgeopat_x']:
-						cost += self.matrice_transition[row['cultgeopat_x']][row['cultgeopat_y']]
+		patches_altered = scenario_patches[scenario_patches['cultgeopat'] != scenario_patches['init_cult']]
+		cost = pd.merge(patches_altered, self.patches_transition_cost,on=['ID_PARCEL','cultgeopat'], how='left')['cost'].sum()
 		return 1/(1+cost)
 
 if __name__ == '__main__':
 	import geopandas as gpd
-	patches = gpd.GeoDataFrame.from_file("../output/PAT_patches/PAT_patches.shp", encoding='utf-8')
-	social = Social(patches, '../resources/matrice_transition.csv', 0.1)
+	import time
+	import os
+	import hashlib
+	from MCMC import md5sum
+	# loading and initializing patches data
+	patches_filename = "../output/PAT_patches/PAT_patches.shp"
+	matrix_filename = '../resources/matrice_transition.csv'
+	patches = gpd.GeoDataFrame.from_file(patches_filename, encoding='utf-8')
+	patches = patches[patches['cultgeopat']!='Non Considérée']
+	patches['init_cult'] = patches['cultgeopat']
+	# init indicator
+	social = Social(patches, md5sum(patches_filename), matrix_filename, md5sum(matrix_filename), 0.1)
+	costs_filename = '../output/patches_transition_costs.npz'
+	if not os.path.isfile(costs_filename):
+		social.compute_patches_transition_cost(patches)
+		social.save_patches_transition_cost(costs_filename)
+	social.load_patches_transition_cost(costs_filename)
+	# computing init score
 	scenario_patches = patches.copy()
 	print(social.compute_indicator(scenario_patches))
-	scenario_patches.ix[2,'cultgeopat']='Protéagineux'
-	print(social.compute_indicator(scenario_patches))
-	# should give 0.5 the cost for one transition to Protéagineux
-	scenario_patches['cultgeopat']='Fruits et légumes'
+	# small reallocation and computing score
+	scenario_patches.loc[1,'cultgeopat']='Fruits et légumes'
+	# should give 1/(1+1) the cost for one transition to Fruits et légumes
 	print(social.compute_indicator(scenario_patches))
+	# computing with more reallocation
+	for N in [2000,5000,10000]:
+		for i in range(N):
+			scenario_patches.ix[i,'cultgeopat']='Protéagineux'
+		print('{} {}'.format(N, social.compute_indicator(scenario_patches)))