MCMC.py

#!/usr/bin/python
# -*- coding: utf-8 -*-
import pandas as pd
import geopandas as gpd
import os, sys, time, shutil
import yaml
import numpy as np
from indicators import Indicators
from tqdm import tqdm
import patutils
from patutils import md5sum, load_pat_patches
import seaborn as sns
import matplotlib.pyplot as plt
from scenarios import ScenariosStack, Reallocator

class MCMC:
    def __init__(self, mcmc_config_filename, output=True):
        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.patches_md5sum = md5sum(self.mcmc_config['patches'])
        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()
        # Copying input data in output dir
        self.outputdir = self.mcmc_config['output_dir'] + '/' + time.strftime('%Y%m%d-%H%M%S')
        if output:
            if not os.path.exists(self.outputdir):
                os.makedirs(self.outputdir)
                print('All data will be written in {}'.format(self.outputdir))
            else:
                print('Output directory already exists! ({})'.format(self.outputdir))
                sys.exit(1)
            shutil.copy(mcmc_config_filename, self.outputdir)
            for f in [ self.mcmc_config['target'],
                    self.mcmc_config['indicators_config']['resilience'],
                    self.mcmc_config['indicators_config']['social']['cost_matrix_filename']
                    ]:
                shutil.copy(f, self.outputdir)
            with open(self.outputdir+'/seed.txt', 'w') as outfile:
                outfile.write('{}\n'.format(self.rng.get_state()))
            config_data = {
                'patches_md5sum':self.patches_md5sum,
                'biodiversity_matrix_md5sum':md5sum(self.mcmc_config['indicators_config']['biodiversity']['matrixfilename']),
                'social_patches_costs_md5sum':md5sum(self.mcmc_config['indicators_config']['social']['patches_costs_filename'])
            }
            with open(self.outputdir+'/config.yml', 'w') as outfile:
                yaml.dump(config_data, outfile, default_flow_style=False)
        # finishing init
        self.patches = load_pat_patches(self.mcmc_config['patches'])
        self.surfaces = self.patches['SURF_PARC']
        self.target = pd.read_csv(self.mcmc_config['target'], sep=';',index_col=0).rename(index=patutils.code_cultgeopat)
        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)
        self.reallocator = Reallocator(self.rng, self.patches, self.targetPAT, self.mcmc_config['ratio_patches_to_modify'])


    def is_pareto_efficient(costs, return_mask = True):
        """
        Find the pareto-efficient points
        :param costs: An (n_points, n_costs) array
        :param return_mask: True to return a mask
        :return: An array of indices of pareto-efficient points.
            If return_mask is True, this will be an (n_points, ) boolean array
            Otherwise it will be a (n_efficient_points, ) integer array of indices.
        Author:
        Source: https://stackoverflow.com/questions/32791911/fast-calculation-of-pareto-front-in-python
        """
        is_efficient = np.arange(costs.shape[0])
        n_points = costs.shape[0]
        next_point_index = 0  # Next index in the is_efficient array to search for
        while next_point_index<len(costs):
            nondominated_point_mask = np.any(costs<costs[next_point_index], axis=1)
            nondominated_point_mask[next_point_index] = True
            is_efficient = is_efficient[nondominated_point_mask]  # Remove dominated points
            costs = costs[nondominated_point_mask]
            next_point_index = np.sum(nondominated_point_mask[:next_point_index])+1
        if return_mask:
            is_efficient_mask = np.zeros(n_points, dtype = bool)
            is_efficient_mask[is_efficient] = True
            return is_efficient_mask
        else:
            return is_efficient

    def scenario_scores(patches, indicators):
        return indicators.compute_indicators(patches)

    def step(self, iter_nb, nb_particles, scenarios):
        '''
        Sample new scenarios, evaluate their scores and retain only pareto front.

        :param iter_nb: (int) number of this iteration.
        :param nb_particles: number of new scenarios to sample
        :param scenarios: (ScenariosStack) list of the scenarios used as base for sampling the new scenarios.
        '''
        # Loop of sampling and scoring
        for index, scenario in tqdm(scenarios.sample(nb_particles, self.rng).iterrows(), total=nb_particles):
            scenario = scenarios.reallocate(index)
            if scenario is None:
                scenarios.setFullReallocated(index, True)
            else:
                scenarios.append(iter_nb, scenario, index)
        # Computing pareto front
        start_time = time.time()
        pareto_mask = MCMC.is_pareto_efficient(scenarios.scores[['Resilience','Proximity','Productivity','Biodiversity','Social','TargetDelta']].values)
        scenarios.scores['pareto'] = pareto_mask
        elapsed_time = time.time() - start_time
        print(' Fully reallocated scenarios {}/{} - '.format(scenarios.nbFullReallocated(), len(scenarios.scores)), end="", flush=True)
        print('{:.2f}% on pareto front (computed in {}) - '.format(100*np.sum(pareto_mask)/len(pareto_mask), time.strftime("%M:%S", time.gmtime(elapsed_time))), end="", flush=True)
        # Writing output data
        write_times=[]
        start_time = time.time()
        scenarios.scores.to_csv(self.outputdir+'/scores_iter_{0:03d}.csv'.format(iter_nb), index=True)
        write_times.append(time.time() - start_time)
        start_time = time.time()
        scenarios.cultgeopat.to_hdf(self.outputdir+'/cult_iter_{0:03d}.h5'.format(iter_nb), 'cultgeopat', index=True, complevel=1)
        write_times.append(time.time() - start_time)
        if self.mcmc_config.get('write_data', False):
            start_time = time.time()
            try:
                pairplot = sns.pairplot(scenarios.scores, vars=['Resilience','Proximity','Productivity','Biodiversity','Social','TargetDelta'],
                        diag_kind="kde",hue="pareto"
                    )
                pairplot.savefig(self.outputdir+'/mcmc_iter_{0:03d}.png'.format(iter_nb))
                plt.close(pairplot.fig)
            except (np.linalg.linalg.LinAlgError, ValueError, AttributeError) as e:
                # when data doesn't vary enough on a dimension, it is impossible to generate the pareto density
                # and an error is triggered. Here, we ignore this error.
                pass
            write_times.append(time.time() - start_time)
            print('Scores written in {}/{}/{} '.format(*map(lambda x: time.strftime("%M:%S", time.gmtime(x)), write_times)), end="", flush=True)
        else:
            print('Scores written in {}/{} '.format(*map(lambda x: time.strftime("%M:%S", time.gmtime(x)), write_times)), end="", flush=True)
        print()
        # Retaining only optimal particules
        scenarios.retain(pareto_mask)
        return scenarios

    def write_data(self, iter_nb, scenarios):
        print(' Writing modified patches of scenarios')
        # Writing shapefiles
        start_time = time.time()
        shp_dir = self.outputdir + '/patches_iter_{0:03d}'.format(iter_nb)
        os.makedirs(shp_dir)
        counter = 0
        for index,scenario in tqdm(scenarios.cultgeopat.iterrows(), total=len(scenarios.cultgeopat)):
            #if scenarios.scores.loc[index,'full_reallocated']:
            counter += 1
            concat = pd.concat([self.patches, scenario.to_frame("newcult")], axis=1)
            concat = concat[concat['init_cult']!=concat['newcult']]
            if len(concat)>0:
                for c in ['newcult','cultgeopat','init_cult']:
                    patutils.decode_cultgeopat(concat, c)
                concat.to_file(shp_dir+'/{:03d}_{}'.format(iter_nb, index), encoding='utf-8')
        elapsed_time = time.time() - start_time
        #print('  {} shapefiles (only scenario fully reallocated) written in {}'.format(counter, time.strftime("%M:%S", time.gmtime(elapsed_time))))
        print('  {} shapefiles written in {}'.format(counter, time.strftime("%M:%S", time.gmtime(elapsed_time))))


    def run(self):
        # Initial sampling and evaluation
        nb_particles = self.mcmc_config['initial_nb_particles']
        # Start with initial scenario
        scenarios = ScenariosStack(self.indicators, self.reallocator, self.patches)
        # Iteration
        nb_iteration=0
        for i in range(self.mcmc_config['max_iterations']):
            nb_iteration=i+1
            print('Iteration #{}'.format(nb_iteration))
            # Write SHP only for the last iteration
            scenarios = self.step(nb_iteration, nb_particles, scenarios)
            if scenarios.allFullReallocated():
                print('All scenarios are fully reallocated. Simulation finished.')
                break
        self.write_data(nb_iteration, scenarios)
        return scenarios

if __name__ == '__main__':
    mcmc = MCMC('MCMC_config.yml')
    scenario = mcmc.run()