MCMC_class.py 37.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#!/usr/bin/python
# -*- coding: utf-8 -*-

import fiona
import numpy as np
import re
from shapely.geometry import mapping, shape
from shapely.ops import cascaded_union
from fiona import collection
import geopandas as gpd
import random
from random import randint
import json
import gzip
import sqlite3
16
17
import sys, os
import yaml
18

19
20
21
from proximite import Proximite
from resilience_list import Resilience
from productivite import Productivity
22
#from indice_biodiversite import biodiversite
23
24
from indice_biodiversite_2 import Biodiversity
from social import Social
25
26
27
28
from EPCI_culture_repartition import EPCI_culture_repartition

from radar import scenarios_graph

29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
class Indicators:
	def __init__(self, mcmc_config_filename, initial_patches):
			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'))
			self._proximity = Proximite(initial_patches, self.mcmc_config['indicators_config']['proximite']['surf_totale_cible'])
			self._resilience = Resilience(self.mcmc_config['indicators_config']['resilience'])
			self._productivity = Productivity()
			biodiv_config = self.mcmc_config['indicators_config']['biodiversity']
			matrixfilename = biodiv_config['matrixfilename']
			self._biodiversity = Biodiversity(initial_patches, biodiv_config['dmax'], biodiv_config['epsilon'])
			self._biodiversity.load_probabilities(matrixfilename)
			social_config = self.mcmc_config['indicators_config']['social']
			self._social = Social(initial_patches, social_config['cost_matrix_filename'], social_config['bdv_threshold'])

	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)

61
62
63
64
class Scenario :
	'''
	This class is used to create a new randomized scenario from an initial scenario. A few functions are included to permit
	the manipulation of this object.
65

66
67
68
69
70
71
72
	The particularity of this class is the way he permit to chose a direction for the random. You have to chose the proportion
	of each kind of cultural patches and the algorithme will generate a new scenario who will try to approach this proportion.
	'''

	def return_culture_code(self, culture):
		'''
		Return the code associated to the kind of cultural type put as a parameter of the function
73

74
75
76
		:param culture: Take the value of a kind of cultural type (Cereales, Culture industrielles, Fruits et legumes
						Oleagineux, Prairies, Proteagineux)
		'''
77
		if culture == "Céréales":
78
79
			#print culture
			return 0
80

81
82
		if culture == "Cultures industrielles"  :
			return 1
83

84
85
		if culture == "Fourrages" :
			return 2
86

87
		if culture == "Fruits et légumes" :
88
89
			#print culture
			return 3
90

91
		if culture == "Oléagineux" :
92
93
			#print culture
			return 4
94

95
96
		if culture == "Prairies" :
			return 5
97

98
		if culture == "Protéagineux" :
99
			return 6
100

101
	def __init__(self, indicators, scenarioInitial, nbr_parcelles_a_modifier,  cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste, rng) :
102
103
			'''
			Create a new randomized scenario according to the parameters
104

105
106
			:param scenarioInitial: the new scenario will originate from this one, he will take his values then change them
			:param nbr_parcelles_a_modifier: the number of patches who will be changed from the initial scenario
107
			:param cereales, culturesI, fruits_legumes, oleagineux, prairies, proteagineux: the proportion of area of each
108
			kind of cultural patches we want to approach in the new scenario
109
			:param liste: if 'ScenarioInitial' is a filename, chose False and if it's a list of patches, chose True. It's helpful
110
111
112
113
114
115
116
117
118
119
120
			to manage both kind of initial scenario.
			'''
			self.shape = scenarioInitial
			self.nbr_parcelles_a_modifier = nbr_parcelles_a_modifier
			self.cereales = cereales
			self.culturesI = culturesI
			self.fourrages = fourrages
			self.fruits_legumes = fruits_legumes
			self.oleagineux = oleagineux
			self.prairies = prairies
			self.proteagineux = proteagineux
121
122

			self.indicators = indicators
123

124
125
			self.initial_patches = gpd.GeoDataFrame.from_file(scenarioInitial, encoding='utf-8')

126
127
128
			# If ScenarioInitial is a filename -> transform it in a list
			if liste == False :
				scenario = []
129
				with fiona.open(self.shape, 'r', encoding='utf-8') as input:
130
131
132
133
134
135
136
137
138
					for f in input :
						scenario.append(f)
					self.shape = scenario
			else:
				scenario = scenarioInitial
			# Sort the initial list because the work done later will change the order. This way, we will be able to order the final list too et not be confused when we use the FID.
			newList = sorted(scenario, key=lambda k: k['properties']['ID_PARCEL'])
			self.shape = newList
			input = self.shape
139

140
			resultat_final = []
141

142
			nombre_culture_a_modifie = self.nbr_parcelles_a_modifier
143

144
145
146
147
148
149
150
151
152
153
			#Area wanted for each kind of cultural patches for the simulation
			# 0 : Céréales, 1 : culture industrielle, 2 : fourrages, 3 : fruits et légumes, 4 : oléagineux, 5 : prairies, 6 : protéagineux
			pourcentage_surface = {}
			pourcentage_surface[0] = self.cereales #Céréales
			pourcentage_surface[1] = self.culturesI # Cultures industrielles
			pourcentage_surface[2] = self.fourrages # Fourrages
			pourcentage_surface[3] = self.fruits_legumes # Fruits et légumes
			pourcentage_surface[4] = self.oleagineux # Oléagineux
			pourcentage_surface[5] = self.prairies # Prairies
			pourcentage_surface[6] = self.proteagineux # Protéagineux
154

155
156
			#Initialize the area of each kind of cultural patches to 0. They will be incremented accordingly of the state of the scenario
			cultures = {}
157

158
159
160
			for i in range (8) :
				cultures[i] = { 'surf': 0.0, 'parcelles': []}

161

162
163
			#Total area of PAT patches
			total_culture = 0.0 #153256
164

165
			i = 0
166
			for f in input :
167
				total_culture += float(f["properties"]["SURF_PARC"])
168

169
				# Initialize the area of the initial scenario for each kind of culture and store the FID of patches accordingly of the culture of the patch
170
				culture = f["properties"]["cultgeopat"]
171
				if culture != None and len(culture)>0:
172
173
174
175
176
177
178
					cultures[self.return_culture_code(culture)]['surf'] += float(f["properties"]["SURF_PARC"])
					cultures[self.return_culture_code(culture)]['parcelles'].append(i)
				else :
					cultures[7]['parcelles'].append(i)
					cultures[7]['surf'] += float(f["properties"]["SURF_PARC"])
				#Increment the FID of patches
				i += 1
179

180
181
			liste_cultures_a_reduire = []
			liste_cultures_a_augmenter = []
182

183
184
185
186
187
			for i in range(7) :
				if cultures[i]['surf']*100/total_culture <= pourcentage_surface[i] :
					liste_cultures_a_augmenter.append(i)
				else :
					liste_cultures_a_reduire.append(i)
188

189
			culture_tab = ["Céréales", "Cultures industrielles", "Fourrages", "Fruits et légumes", "Oléagineux", "Prairies", "Protéagineux", None]
190

191
			cpt = 0
192
			while cpt < nombre_culture_a_modifie and len(liste_cultures_a_augmenter) + len(liste_cultures_a_reduire) > 1 :
193
194
195
						#Get a random patches who need to reduce in area
						if len(liste_cultures_a_reduire) != 0 :
							rand_cul = rng.choice(liste_cultures_a_reduire)
196
						else :
197
198
199
200
201
							# Case when the desired area is reached in all kind of cultural patches but there is still patches to modify
							# We add a random culture to increase in the list of culture we need to reduce, so the algorithme go on until he reach "nombre_culture_a_modifie"
							rand_cul = rng.choice(liste_cultures_a_augmenter)
							liste_cultures_a_augmenter.remove(rand_cul)
							liste_cultures_a_reduire.append(rand_cul)
202

203
204
205
206
						if len(liste_cultures_a_augmenter) == 0 :
							rand_augmenter = rng.choice(liste_cultures_a_reduire)
							liste_cultures_a_augmenter.append(rand_augmenter)
							liste_cultures_a_reduire.remove(rand_augmenter)
207
208


209
						random_parcelles = rng.choice(cultures[rand_cul]['parcelles'])
210

211
						code_culture_rand = input[random_parcelles]["properties"]["cultgeopat"]
212

213
214
215
216
217
218
219
220
						recW = {}
						recW["geometry"]=input[random_parcelles]["geometry"]
						recW["properties"] = {}
						recW["properties"]['ID_PARCEL']=input[random_parcelles]["properties"]["ID_PARCEL"]
						recW["properties"]['CODE_CULTU']=input[random_parcelles]["properties"]["CODE_CULTU"]
						recW["properties"]['SURF_PARC']=input[random_parcelles]["properties"]["SURF_PARC"]
						recW["properties"]['Bdv']=input[random_parcelles]["properties"]["Bdv"]
						recW["properties"]['INSEE_COM']=input[random_parcelles]["properties"]["INSEE_COM"]
221
						recW["properties"]['VALEUR_CAD']=input[random_parcelles]["properties"]["VALEUR_CAD"]
222

223
						succes = 0
224
						while succes == 0 :
225
							rand_cult = rng.choice(liste_cultures_a_augmenter)
226

227
							cultures[rand_cult]['surf'] += float(input[random_parcelles]["properties"]["SURF_PARC"] )
228
							recW["properties"]['cultgeopat'] = culture_tab[rand_cult]
229
230

							#Remove the patch 'random_parcelles' of 'liste_nombre'
231
							cultures[self.return_culture_code(input[random_parcelles]["properties"]["cultgeopat"])]['parcelles'].remove(random_parcelles)
232

233
							#Remove the area of the old value contained in 'culture[type_culture]
234
							cultures[self.return_culture_code(input[random_parcelles]["properties"]["cultgeopat"])]['surf'] -= float(input[random_parcelles]["properties"]["SURF_PARC"])
235
							resultat_final.append(recW)
236

237
238
239
240
241
							# If the proportion of the area of one kind of culture exceeds the proportion wanted
							# We remove it from the list of culture to increase and we add in the one list of the one to reduce
							if float(cultures[rand_cult]['surf']*100/total_culture) >= pourcentage_surface[rand_cult] :
								liste_cultures_a_augmenter.remove(rand_cult)
								liste_cultures_a_reduire.append(rand_cult)
242

243
							# Same than last block but in reverse, reduce_list to increase_list
244
245
246
							if float(cultures[self.return_culture_code(input[random_parcelles]["properties"]["cultgeopat"])]['surf']*100/total_culture) <= pourcentage_surface[self.return_culture_code(input[random_parcelles]["properties"]["cultgeopat"])] :
								liste_cultures_a_reduire.remove(self.return_culture_code(input[random_parcelles]["properties"]["cultgeopat"]))
								liste_cultures_a_augmenter.append(self.return_culture_code(input[random_parcelles]["properties"]["cultgeopat"]))
247

248
249
							succes = 1
							cpt += 1
250

251
252
253
							# Remove the comment to see the proportion change during the algorithm
							#for i in range(7) :
							#print cultures[i]['surf']*100/total_culture
254

255
256
257
258
						# If all patches of a kind of culture as already been modified, we remove this culture type from the list (increase_list or reduce_list)
						for i in range(7):
							if len(cultures[i]["parcelles"]) == 0 :
								if i in liste_cultures_a_reduire :
259
									liste_cultures_a_reduire.remove(i)
260
261
								elif i in liste_cultures_a_augmenter :
									liste_cultures_a_augmenter.remove(i)
262
263

			# Every patch remaining will be added in the new layer. Patches with an empty type will be added now
264
265
266
267
268
269
270
271
272
			for i in range(8) :
				for e in cultures[i]["parcelles"] :
					recW = {}
					recW["geometry"] = input[e]["geometry"]
					recW["properties"] = {}
					recW["properties"]['ID_PARCEL'] = input[e]["properties"]["ID_PARCEL"]
					recW["properties"]['CODE_CULTU'] = input[e]["properties"]["CODE_CULTU"]
					recW["properties"]['SURF_PARC'] = input[e]["properties"]["SURF_PARC"]
					recW["properties"]['Bdv'] = input[e]["properties"]["Bdv"]
273
					recW["properties"]['id_expl'] = input[e]["properties"]["id_expl"]
274
					recW["properties"]['INSEE_COM'] = input[e]["properties"]["INSEE_COM"]
275
276
					recW["properties"]['VALEUR_CAD'] = input[e]["properties"]["VALEUR_CAD"]
					recW["properties"]['cultgeopat'] = culture_tab[i]
277

278
					resultat_final.append(recW)
279

280
			#print
281
			#print ("Pourcentage culture apres MCMC : ")
282
283
			#for i in range(7) :
			#print (cultures[i]['surf']*100/total_culture)
284

285
286
287
			# Sort the final list to make it ordered as the same way of the initial list
			resultat_final = sorted(resultat_final, key=lambda k: k['properties']['ID_PARCEL'])
			self.resultat_final = resultat_final
288
			self.resultat_final_gpd = gpd.GeoDataFrame.from_features(resultat_final)
289
290
291



292
293
	def get_scenario(self):
		'''
294
		Return the object himself
295
296
		'''
		return self
297
298


299
300
301
302
303
	def get_resultat_final(self):
		'''
		Return all entities of the scenario in the list format (faster to work on than a shapefile)
		'''
		return self.resultat_final
304
305


306
307
308
309
310
	def write_in_shape(self, shape) :
		'''
		Write the scenario in a shapefile, permits to visualize the map
		:param shape: the name of the file which will contains the scenario
		'''
311
		schema = {'geometry': 'Polygon', 'properties': {'CODE_CULTU' : 'str', 'Bdv' : 'str', 'ID_PARCEL' : 'str', 'INSEE_COM' : 'str', 'SURF_PARC' : 'float', 'cultgeopat' : 'str'}}
312
		crs = {'lon_0': 3, 'ellps': 'GRS80', 'y_0': 6600000, 'no_defs': True, 'proj': 'lcc', 'x_0': 700000, 'units': 'm', 'lat_2': 44, 'lat_1': 49, 'lat_0': 46.5}
313
		driver = "ESRI Shapefile"
314

315
		with fiona.open(shape,'w',driver=driver, crs=crs, schema=schema, encoding='utf-8') as parcelles :
316
317
			for f in self.resultat_final :
				parcelles.write(f)
318

319
320
321
322
323
324
325
326
327
328
	def write_in_json(self):
		'''
		*** not used, can be removed/modified ***
		Write the scenario in a Json file
		'''
		data = self.resultat_final

		with gzip.GzipFile('./best_scenario.gz','w') as fid_gz:
			json_str = json.dumps(self.resultat_final)
			fid_gz.write(json_str)
329

330
331
332
333
334
335
336
	def write_indices_in_json(self):
		'''
		Write the values of indices of the initial scenario in a JSON file
		This function is used for display and statistics in the application
		'''
		dict_indices = {}
		dict_indices["proximite"] = self.indice_proximite()
337
338
339
340
		dict_indices["resilience"] = self.indice_resilience()
		dict_indices["productivite"] = self.indice_productivitee()
		dict_indices["biodiversite"] = self.indice_biodiversite()
		dict_indices["social"] = self.indice_social()
341

342
		with gzip.GzipFile('./initial_scenario_indices.gz','w') as fid_gz:
343

344
			json_str = json.dumps(dict_indices, ensure_ascii=False)
Dumoulin Nicolas's avatar
Dumoulin Nicolas committed
345
			fid_gz.write(json_str.encode('utf-8'))
346

347
348
	'''
	These functions returns values of indices of this scenario
349
350
	'''
	def indice_proximite(self) :
351
		return self.indicators.proximity(self.resultat_final_gpd)
352

353
	def indice_resilience(self):
354
		return self.indicators.resilience(self.resultat_final_gpd)
355

356
	def indice_productivitee(self):
357
		return self.indicators.productivity(self.resultat_final_gpd)
358

359
	def indice_biodiversite(self):
360
		return self.indicators.biodiversity(self.resultat_final_gpd)
361

362
	def indice_social(self):
363
		return self.indicators.social(self.resultat_final_gpd)
364

365
366
367
368
369
370
371
372
373
374
375
dict_scenario = {} # Contains all scenarios created by the function "générer_x_scenarios"

'''
These few function are used to get some data and make stats. They are not called in the main program but as they can be useful, we keep them.
'''
def generer_x_scenarios(x, scenario_initial, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste ): #Génère x scénarios et les ranges dans "dict_scenario"
	'''
	Generate X scenarios and stocks them in a dict. It will be used to compare all scenario in the dict and keep the betters one
	'''
	dict_scenario = {}
	for i in range(x):
376
		s = Scenario(scenario_initial, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste)
377
378
		s_r = s.get_scenario()
		dict_scenario[i] = s_r
379

380
381
382
383
384
385
386
387
388
389
390
391
392
	#comparer_2_scenarios(dict_scenario[0], dict_scenario[1])
	return dict_scenario


def x_all_indices(scenarInitial, x, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux):
	'''
	Create X numbers of scenario, each new scenario is created from an initial one.
	'''
	scenar = scenarInitial
	liste = False
	for i in range(x) :
			s = Scenario(scenar, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste)
			s.indice_proximite()
393
394
395
			s.indice_resilience()
			s.indice_productivitee()
			s.indice_biodiversite()
396
397
			scenar = s.get_resultat_final()
			liste = True
398

399
400
401
402
403
404
405
406
407
408
409

def x_proximite(scenarInitial, x, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux):
	'''
	Calculate X number of proximity indices
	Create an X numbers of scenario, each new scenario is created from an initial one. It will create 4 set of data and put them into a file readable as a csv.
	Each set of data change a different number of patches.
	'''
	scenar = scenarInitial
	liste = False
	fichier = open("./indice_proximite_300.txt", "w")
	for i in range(x) :
410
			print(i)
411
412
413
414
415
			list = [300, 1000, 5000, 10000]
			for nbr in list :
				nbr_parcelles_a_modifier = nbr
				s = Scenario(scenar, nbr, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste)
				ind = s.indice_proximite()
416

417
418
419
420
				if nbr != 10000:
					fichier.write(str(ind) + ";")
				else :
					fichier.write(str(ind) + "\n")
421
422


423
424
425
426
427
428
429
430
431
432
def x_resilience(scenarInitial, x, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux):
	'''
	Calculate X number of resiliency indices
	Create an X numbers of scenario, each new scenario is created from an initial one. It will create 4 set of data and put them into a file readable as a csv.
	Each set of data change a different number of patches.
	'''
	scenar = scenarInitial
	liste = False
	fichier = open("./indice_resilience.txt", "w")
	for i in range(x) :
433
			print(i)
434
435
436
437
438
			list = [300, 1000, 5000, 10000]
			for nbr in list :
				nbr_parcelles_a_modifier = nbr
				s = Scenario(scenar, nbr, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste)
				#print "Scenario " + str(i) + " : "
439
				ind = s.indice_resilience()
440
441
442
443
444
445
				#scenar = s.get_resultat_final()
				#liste = True
				if nbr != 10000:
					fichier.write(str(ind) + ";")
				else :
					fichier.write(str(ind) + "\n")
446
447


448
449
450
451
452
453
454
455
456
457
def x_productivite(scenarInitial, x, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux):
	'''
	Calculate X number of productivity indices
	Create an X numbers of scenario, each new scenario is created from an initial one. It will create 4 set of data and put them into a file readable as a csv.
	Each set of data change a different number of patches.
	'''
	scenar = scenarInitial
	liste = False
	fichier = open("./indice_productivite.txt", "w")
	for i in range(x) :
458
			print(i)
459
460
461
462
			list = [300, 1000, 5000, 10000]
			for nbr in list :
				s = Scenario(scenar, nbr, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste)
				#print "Scenario " + str(i) + " : "
463
				ind = s.indice_productivitee()
464
465
466

				#scenar = s.get_resultat_final()
				#liste = True
467

468
469
470
471
				if nbr != 10000:
					fichier.write(str(ind) + ";")
				else :
					fichier.write(str(ind) + "\n")
472
473


474
475
476
477
478
479
480
481
482
483
def x_biodiversite(scenarInitial, x, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux):
	'''
	Calculate X number of biodiversity indices
	Create an X numbers of scenario, each new scenario is created from an initial one. It will create 4 set of data and put them into a file readable as a csv.
	Each set of data change a different number of patches.
	'''
	scenar = scenarInitial
	liste = False
	fichier = open("./indice_biodiversite.txt", "w")
	for i in range(x) :
484
			print(i)
485
486
487
488
			list = [300, 1000, 5000, 10000]
			for nbr in list :
				s = Scenario(scenar, nbr, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste)
				#print "Scenario " + str(i) + " : "
489
				ind = s.indice_biodiversite()
490
491
492
493
494
495
				#scenar = s.get_resultat_final()
				#liste = True
				if nbr != 10000:
					fichier.write(str(ind) + ";")
				else :
					fichier.write(str(ind) + "\n")
496

497
498
499
500
501
502
503
504
505
506
def x_social(scenarInitial, x, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux):
	'''
	Calculate X number of social indices
	Create an X numbers of scenario, each new scenario is created from an initial one. It will create 4 set of data and put them into a file readable as a csv.
	Each set of data change a different number of patches.
	'''
	scenar = scenarInitial
	liste = False
	fichier = open("./indice_social.txt", "w")
	for i in range(x) :
507
			print(i)
508
509
510
511
			list = [300, 1000, 5000, 10000]
			for nbr in list :
				s = Scenario(scenar, nbr, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste)
				#print "Scenario " + str(i) + " : "
512
				ind = s1.indice_social()
513
514
515
516
517
518
519
				#scenar = s.get_resultat_final()
				#liste = True
				if nbr != 10000:
					fichier.write(str(ind) + ";")
				else :
					fichier.write(str(ind) + "\n")

520
521


522
523
524
525
526
527
528
529
def comparer_2_scenarios(scenar1, scenar2) :
	'''
	*** Not really useful anymore ***
	Compare 2 scenarios who are in the global variable 'dict_scenario'
	'''
	dico = {}
	dico[0] = {}
	dico[1] = {}
530

531
532
533
534
	s1_proximite = scenar1.indice_proximite()
	dico[0]["proximite"] = s1_proximite
	s2_proximite = scenar2.indice_proximite()
	dico[1]["proximite"] = s2_proximite
535

536
	s1_resilience = scenar1.indice_resilience()
537
	dico[0]["resilience"] = s1_resilience
538
	s2_resilience = scenar2.indice_resilience()
539
	dico[1]["resilience"] = s2_resilience
540

541
	s1_productivite = scenar1.indice_productivitee()
542
	dico[0]["productivite"] = s1_productivite
543
	s2_productivite = scenar2.indice_productivitee()
544
	dico[1]["productivite"] = s2_productivite
545

546
	s1_biodiversite = scenar1.indice_biodiversite()
547
	dico[0]["biodiversite"] = s1_biodiversite
548
	s2_biodiversite = scenar2.indice_biodiversite()
549
	dico[1]["biodiversite"] = s2_biodiversite
550

551
552
	dico[0]["social"] = 0
	dico[1]["social"] = 0
553
554

	print()
555
	if s1_proximite > s2_proximite :
556
557
558
559
560
		print("Le scénario 1 possède un meilleur indice de proximité que le scenario 2 de : " + str((s1_proximite - s2_proximite) / s1_proximite * 100) + "%")
	else :
		print("Le scénario 2 possède un meilleur indice de proximité que le scénario 1 de : " + str((s2_proximite - s1_proximite) / s2_proximite * 100) + "%")
	print()

561
	if s1_resilience > s2_resilience :
562
563
564
565
566
		print("Le scénario 1 possède un meilleur indice de résilience que le scenario 2 de : " + str((s1_resilience - s2_resilience) / s1_resilience * 100) + "%")
	else :
		print("Le scénario 2 possède un meilleur indice de résilience que le scénario 1 de : " + str((s2_resilience - s1_resilience) / s2_resilience * 100) + "%")
	print()

567
	if s1_productivite > s2_productivite :
568
569
570
571
		print("Le scénario 1 possède un meilleur indice de productivité que le scenario 2 de : " + str((s1_productivite - s2_productivite) / s1_productivite * 100) + "%")
	else :
		print("Le scénario 2 possède un meilleur indice de productivité que le scénario 1 de : " + str((s2_productivite - s1_productivite) / s2_productivite * 100) + "%")

572
	if s1_biodiversite > s2_biodiversite :
573
574
575
576
577
		print("Le scénario 1 possède un meilleur indice de biodiversite  que le scenario 2 de : " + str((s1_biodiversite - s2_biodiversite) / s1_biodiversite  * 100) + "%")
	else :
		print("Le scénario 2 possède un meilleur indice de biodiversite  que le scénario 1 de : " + str((s2_biodiversite - s1_biodiversite) / s2_biodiversite  * 100) + "%")


578
579
580
581
582
583
584
585
586
	scenarios_graph(2, dico)




#Use to stock in JSON file the data about the area of each EPCI for each best scenario
def EPCI_repartition(best_scenario):
	'''
	Function used to get the area of each kind of culture for each EPCI in the PAT territory
587

588
589
590
591
	:param best_scenario: dictionnary of the final best scenario produced by the MCMC function
	'''
	list_EPCI_scenario = []
	EPCI_layer = []
592
	with fiona.open("./EPCI_clean.shp", 'r', encoding='utf-8') as epci :
593
594
595
		for f in epci :
			EPCI_layer.append(f)

596

597
	for i in range(len(best_scenario)) :
598
		print(best_scenario[i]['proximite'])
599
		list_EPCI_scenario.append(EPCI_culture_repartition(EPCI_layer, best_scenario[i]['scenario']))
600

601
602
	with gzip.GzipFile('./EPCI_repartition.gz','w') as fid_gz:
		json_str = json.dumps(list_EPCI_scenario)
603
		fid_gz.write(json_str)
604
605
606
607
608
609
610
611




def get_best_scenario(dict_scenar, nbr_scenar_to_stock):
	'''
	Function used is the MCMC function
	Get the 'nbr_scenar_to_stock' best scenario for each kind of cultural patches from a dictionary 'dict_scenar'
612

613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
	:param dict_scenar: dictionary containing a definite by the MCMC function numbers of scenarios
	:param nbr_scenar_to_stock: number of scenario we want to keep for each indices
	'''
	dict_resultat = {}
	list_indices = ['proximite', 'resilience', 'productivite', 'biodiversite', 'social']
	for i in range(5):
		dict_resultat[i] = {}
		dict_resultat[i]['values'] = []
		dict_resultat[i]['scenario'] = []

	# Add nbr_scenar_to_stock number of scenarios in the dictionnary to start the sort of the best scenarios
	for y in range(nbr_scenar_to_stock) :
		for indices in range(5):
			dict_resultat[indices]['scenario'].append(dict_scenar[y])
			dict_resultat[indices]['values'].append(dict_scenar[y][list_indices[indices]])
628

629
630
631
632
	# Loop who select the best scenarios
	for i in range(len(dict_scenar)):
		for y in range(nbr_scenar_to_stock) :
			for indices in range(5):
633

634
635
636
637
				# Low = best : proximite         high = best : resilience, productivite, biodiversite, social
				# For the proximity indice, the lowest value is the best, a special loop is created for this indice
				if i == 0 :
					if dict_scenar[i][list_indices[indices]] < dict_resultat[indices]['values'][y] and not(dict_scenar[i][list_indices[indices]] in dict_resultat[indices]['values']) :
638

639
640
						del dict_resultat[indices]['values'][y]
						dict_resultat[indices]['values'].append(dict_scenar[i][list_indices[indices]])
641

642
643
						del dict_resultat[indices]['scenario'][y]
						dict_resultat[indices]['scenario'].append(dict_scenar[i])
644

645
				elif dict_scenar[i][list_indices[indices]] > dict_resultat[indices]['values'][y] and not(dict_scenar[i][list_indices[indices]] in dict_resultat[indices]['values']) :
646

647
648
						del dict_resultat[indices]['values'][y]
						dict_resultat[indices]['values'].append(dict_scenar[i][list_indices[indices]])
649

650
651
						del dict_resultat[indices]['scenario'][y]
						dict_resultat[indices]['scenario'].append(dict_scenar[i])
652

653
	return dict_resultat
654

655
656
657
658
659
660
661
662
663
dict_scenario = {}
dump_scenario = {}
dump_cpt = 0

def MCMC(scenarInitial, nbr_loop, nbr_scenario, nbr_scenario_to_stock, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, rng):
	'''
	Final function of the project. His purpose is to find X best scenario from a Y number of generated scenario multiplied by the number of loops we want.
	The first loop will create scenarios from the initial scenario put in the parameters. The others loops will used the already chosen best scenario to create new ones.
	This method will permit to analyse a high variety of different scenario and enhance the size of the spectrum of data created.
664

665
	The function save in an Sqlite database all the scenarios created to keep trace of them. Database name : scenario_dump      Table name : scenario
666

667
668
669
670
671
672
	:param scenarInitial: Initial scenario who will be used to create the scenarios in the first loop
	:param nbr_loop: The number of loop who will run in this algorithme. The first loop is counted as one, so the value need to be at least 1
	:param nbr_scenario: The number of scenario created in each loop
	:param nbr_scenario_to_stock: The number of best scenario saved in a dictionnary at the end of each loop for each indices.
	 The final number of best scenario saved is 'nbr_scenario_to_stock' * number of indices (5)
	:param nbr_parcelles_a_modifier: the number of patches who will be changed from the initial scenario
673
	:param cereales, culturesI, fruits_legumes, oleagineux, prairies, proteagineux: the proportion of area of each
674
	 kind of cultural patches we want to approach in the new scenario
675
676
	:param liste: if 'ScenarioInitial' is a filename, chose False and if it's a list of patches, chose True. It's helpful
	 to manage both kind of initial scenario.
677
678
679
	'''
	liste = True
	global dump_cpt
680

681
682
683
684
685
686
687
	# Put in a sqlite database all the scenario created during the algorithme
	conn = sqlite3.connect('./scenario_dump.db')
	cursor = conn.cursor()
	cursor.execute("""
	DROP TABLE IF EXISTS scenario
	""")
	conn.commit()
688

689
690
691
692
693
	cursor.execute("""
	CREATE TABLE scenario (
		ID integer primary key,
		scenario VARCHAR(10000000),
		indice_proximite double,
694
695
696
697
		indice_resilience double,
		indice_productivite double,
		indice_biodiversite double,
		indice_social double
698
699
700
	)
	""")
	conn.commit()
701
702


703
	scenario = []
704
	with fiona.open(scenarInitial, 'r', encoding='utf-8') as input:
705
706
		for f in input :
			scenario.append(f)
707

708
709
710
711
	# Get the values of indices of the initial scenario
	si = Scenario(scenario, 0, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, True, rng)
	si_indices = {}
	si_indices['proximite'] = si.indice_proximite()
712
713
714
715
	si_indices['resilience'] = si.indice_resilience()
	si_indices['productivite'] = si.indice_productivitee()
	si_indices['biodiversite'] = si.indice_biodiversite()
	si_indices['social'] = si.indice_social()
716

717
718
719
720
721
722
	s_new_indices_avg = {}
	s_new_indices_avg['proximite'] = 0
	s_new_indices_avg['resilience'] = 0
	s_new_indices_avg['productivite'] = 0
	s_new_indices_avg['biodiversite'] = 0
	s_new_indices_avg['social'] = 0
723

724
725
726
727
728
	for nbr_scenar in range(nbr_scenario):
		s = Scenario(scenario, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste, rng)
		dict_scenario[nbr_scenar] = {}
		dict_scenario[nbr_scenar]['scenario'] = s.get_scenario()
		dict_scenario[nbr_scenar]['proximite'] = s.indice_proximite()
729
730
731
732
		dict_scenario[nbr_scenar]['resilience'] = s.indice_resilience()
		dict_scenario[nbr_scenar]['productivite'] = s.indice_productivitee()
		dict_scenario[nbr_scenar]['biodiversite'] = s.indice_biodiversite()
		dict_scenario[nbr_scenar]['social'] = s.indice_social()
733
734


735
736
737
738
739
740
741
742
743
744
745
746
		s_new_indices_avg['proximite'] += dict_scenario[nbr_scenar]['proximite']
		s_new_indices_avg['resilience'] += dict_scenario[nbr_scenar]['resilience']
		s_new_indices_avg['productivite'] += dict_scenario[nbr_scenar]['productivite']
		s_new_indices_avg['biodiversite'] += dict_scenario[nbr_scenar]['biodiversite']
		s_new_indices_avg['social'] += dict_scenario[nbr_scenar]['social']
		#dump_scenario[dump_cpt] = s.get_resultat_final()
		cursor.execute("""
			INSERT INTO scenario(id, scenario, indice_proximite, indice_resilience, indice_productivite, indice_biodiversite, indice_social)
			VALUES(?, ?, ?, ?, ?, ?, ?)
		""",
		(dump_cpt ,str(s.get_resultat_final()), dict_scenario[nbr_scenar]['proximite'], dict_scenario[nbr_scenar]['resilience'], dict_scenario[nbr_scenar]['productivite'], dict_scenario[nbr_scenar]['biodiversite'], dict_scenario[nbr_scenar]['social'])
		)
747

748
		dump_cpt += 1
749
750
751

	conn.commit()

752
	list_indices = ['proximite', 'resilience', 'productivite', 'biodiversite', 'social']
753

754
755
756
757
758
759
760
761
762
	'''
	for i in list_indices :
		print i +" initiale : " + str(si_indices[i])
		print i + " average of new scenarios : " + str(s_new_indices_avg[i] / nbr_scenario)
	'''
	value_end_loop = {}
	for i in list_indices:
		value_end_loop[i] = abs(si_indices[i] - (s_new_indices_avg[i] / nbr_scenario)) / 100
	best_scenario = get_best_scenario(dict_scenario, nbr_scenario_to_stock)
763

764
765
766
767
768
769
770
	#Triggers of end condition
	trigger = {}
	trigger['proximite'] = 0
	trigger['resilience'] = 0
	trigger['productivite'] = 0
	trigger['biodiversite'] = 0
	trigger['social'] = 0
771
772
773

	cpt_loop = 0

774
775
776
	fichier = open("./difference_values.txt", "w")
	list_indices = ['proximite', 'resilience', 'productivite', 'biodiversite', 'social']

777
	while cpt_loop < nbr_loop-1 and (trigger['proximite'] + trigger['resilience'] + trigger['productivite'] + trigger['biodiversite'] + trigger['social']) < 5 :
778
779
780
781
782
783

		trigger['proximite'] = 0
		trigger['resilience'] = 0
		trigger['productivite'] = 0
		trigger['biodiversite'] = 0
		trigger['social'] = 0
784

785
786
787
788
		cpt = 0
		for i in best_scenario:
			for val in best_scenario[i]['scenario']:
				dict_scenario[cpt] = val
789

790
791
792
793
794
		for nbr_scenar in range(nbr_scenario - nbr_scenario_to_stock*5, nbr_scenario) :
			# Create new scenario using the best scenario chosen in the previous step
			s = Scenario(dict_scenario[nbr_scenar%(nbr_scenario_to_stock * 5)]['scenario'].get_resultat_final(), nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste, rng)
			dict_scenario[nbr_scenar]['scenario'] = s.get_scenario()
			dict_scenario[nbr_scenar]['proximite'] = s.indice_proximite()
795
796
797
798
			dict_scenario[nbr_scenar]['resilience'] = s.indice_resilience()
			dict_scenario[nbr_scenar]['productivite'] = s.indice_productivitee()
			dict_scenario[nbr_scenar]['biodiversite'] = s.indice_biodiversite()
			dict_scenario[nbr_scenar]['social'] = s.indice_social()
799

800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
		old_best_scenario = best_scenario
		best_scenario = get_best_scenario(dict_scenario, nbr_scenario_to_stock)
		cpt_loop += 1
		for ind in range(5):
			sum_old = 0
			sum_new = 0
			for y in range(nbr_scenario_to_stock) :
				sum_old += old_best_scenario[ind]['values'][y]
				sum_new += best_scenario[ind]['values'][y]
			difference = abs(((sum_old/nbr_scenario_to_stock) - (sum_new/nbr_scenario_to_stock)) / (sum_new/nbr_scenario_to_stock) * 100)
			fichier.write(str(difference) + ";")
			#print "Difference for old and new " + list_indices[ind] + " is : " + str(difference)
			#print si_indices[list_indices[ind]]
			#print value_end_loop[list_indices[ind]]
			#print
815

816
817
			if difference <= value_end_loop[list_indices[ind]] :
				trigger[list_indices[ind]] = 1
818
819


820
821
822
823
824
	scenario_json = {}
	for i in range(nbr_scenario_to_stock*5):
		scenario_json[i] = {}
	cpt_json = 0

825

826
827
	for i in range(len(best_scenario)):
			for y in range(nbr_scenario_to_stock):
828

829
830
831
832
833
834
				scenario_json[cpt_json]['scenario'] = best_scenario[i]['scenario'][y]['scenario'].get_resultat_final()
				scenario_json[cpt_json]['proximite'] = best_scenario[i]['scenario'][y]['proximite']
				scenario_json[cpt_json]['resilience'] = best_scenario[i]['scenario'][y]['resilience']
				scenario_json[cpt_json]['productivite'] = best_scenario[i]['scenario'][y]['productivite']
				scenario_json[cpt_json]['biodiversite'] = best_scenario[i]['scenario'][y]['biodiversite']
				scenario_json[cpt_json]['social'] = best_scenario[i]['scenario'][y]['social']
835

836
				cpt_json += 1
837
838
839


	EPCI_repartition(scenario_json)
840
841
	with gzip.GzipFile('./best_scenario.gz','w') as fid_gz:
		json_str = json.dumps(scenario_json)
842
843
		fid_gz.write(json_str)

844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
if __name__ == '__main__':
	# Modifie these values to define how many of each cultures you want in your scenario, the modifications will tends to approach these values
	cereales = 22
	culturesI = 5
	fourrages = 5
	fruits_legumes = 3
	oleagineux = 1
	prairies = 59
	proteagineux = 2

	# Change to select how many patches you will change from the initial scenario
	nbr_parcelles_a_modifier = 20000

	# Random seed
	seed = random.randrange(sys.maxsize)
	rng = random.Random(seed)
	print(("Seed was:", seed))

	fichier_seed = open("./lastest_seed.txt", "w")
	fichier_seed.write(str(seed))
	#s1 = Scenario("./Parcelle_PAT_MCMC.shp",nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, False)
	#s1.indice_proximite()
	#s1.indice_resilience('../Parcelle_PAT/Grille_resilience.shp')
	#s1.indice_productivitee("../Parcelle_PAT/Parcelle_PAT_valCad.shp")
	#s1.indice_biodiversite('./compute_dispersal_probabilities_1000.npz')
	#s1.indice_biodiversite('./matrix_biodiversity.npz')
	#s1.indice_social('./Bassin_de_vie_Pat.shp', './Parcelle_PAT_exploitant.shp')
	#s1.write_in_json()

	#s1.write_in_shape('./testSimulationSurf_MCMC.shp')


	#dict_scenario = generer_x_scenarios(2, "./Parcelle_PAT_MCMC.shp", nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, False)

	#x_proximite("./Parcelle_PAT_MCMC.shp", 100 , nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux)
	#x_resilience("./Parcelle_PAT_MCMC.shp", 100 , nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux)
	#x_productivite("./Parcelle_PAT_MCMC.shp", 100 , nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux)
	#x_biodiversite("./Parcelle_PAT_MCMC.shp", 1 , nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux)
	#x_social("./Parcelle_PAT_MCMC.shp", 1 , nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux)

	# scenarInitial, nbr_loop, nbr_scenario, nbr_scenario_to_stock, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux
	#scenario_initial = Scenario("./Parcelle_PAT_MCMC.shp",0, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, False, rng)
886
887
888
	patches = gpd.GeoDataFrame.from_file("../output/PAT_patches/PAT_patches.shp", encoding='utf-8')
	indicators = Indicators('MCMC_config.yml', patches)
	scenario_initial = Scenario(indicators,"../output/PAT_patches/PAT_patches.shp",0, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, False, rng)
889
890
891
	#scenario_initial = Scenario("../output/patches_RiomLimagneVolcans.shp",0, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, False, rng)
	scenario_initial.write_indices_in_json()
	#MCMC("./Parcelle_PAT_MCMC.shp", 2, 11, 2, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, rng)