<commit_before># -*- coding: utf-8 -*-
# /***************************************************************************
# Irmt
#                                 A QGIS plugin
# OpenQuake Integrated Risk Modelling Toolkit
#                              -------------------
#        begin                : 2013-10-24
#        copyright            : (C) 2014 by GEM Foundation
#        email                : <EMAIL>
# ***************************************************************************/
#
# OpenQuake is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# OpenQuake is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake.  If not, see <http://www.gnu.org/licenses/>.

import csv
import os
import time
import matplotlib.pyplot as plt
from collections import defaultdict

from PyQt4.QtCore import pyqtSlot, QSettings, QDir
from PyQt4.QtGui import (QDialog,
                         QFileDialog,
                         QDialogButtonBox)
from qgis.core import QgsMapLayer
from qgis.gui import QgsMessageBar
from svir.calculations.aggregate_loss_by_zone import add_zone_id_to_points

from svir.utilities.utils import (get_ui_class,
                                  read_config_file,
                                  reload_layers_in_cbx,
                                  reload_attrib_cbx,
                                  WaitCursorManager,
                                  tr,
                                  create_progress_message_bar,
                                  clear_progress_message_bar,
                                  )
from svir.recovery_modeling.building import Building

FORM_CLASS = get_ui_class('ui_recovery_modeling.ui')

DAYS_BEFORE_EVENT = 200
SVI_WEIGHT_COEFF = 1  # FIXME: Let the user set this parameter


class RecoveryModelingDialog(QDialog, FORM_CLASS):
    """
    Modal dialog giving to perform recovery modeling analysis.

    Object-oriented Programming(OOP) is used to model post-earthquake community
    recovery of residential community.

    PAOLO: FIXME In the last version, we only have a Building class
    In this Programming, four classes are
    adopted, namely Building Class, Ground Motion FIelds Class, Constructuon &
    Enigneering Service Class and Community Class.

    In each class, a set of data and functions work together to complete a task
    which we call attributes and functions.  Attribute is used to specify if
    and how object properties are accessed within the programming environment
    Function is used to create/altere to achieve task or functionality and
    return the attributes of an object

    Procedure of this OOP Programming:
    Step 1: Define attributes of objects in each class
    Step 2: Set up functions in each class
    Step 3: Incorporate Napa Data to community recovery model

    Methodology:
    Time-based method is utilized which characterize a probability density
    function of the time it takes to a higher or lower functioning state given
    a set of explanatory variables such as the extent of damage to the
    building.
    """
    def __init__(self, iface):
        QDialog.__init__(self)
        self.iface = iface
        # Set up the user interface from Designer.
        self.setupUi(self)
        self.ok_button = self.buttonBox.button(QDialogButtonBox.Ok)
        self.populate_approach_cbx()
        reload_layers_in_cbx(
            self.dmg_by_asset_layer_cbx, [QgsMapLayer.VectorLayer])
        reload_layers_in_cbx(self.svi_layer_cbx, [QgsMapLayer.VectorLayer])
        recovery_modeling_dir = os.path.join(
            os.path.dirname(os.path.abspath(__file__)),
            '..', 'recovery_modeling')
        self.input_data_dir = os.path.join(recovery_modeling_dir, 'input_data')
        self.restoreState()
        self.set_ok_button()

    def restoreState(self):
        """
        Reinstate the options based on the user's stored session info.
        """
        mySettings = QSettings()
        output_data_dir = mySettings.value('irmt/recovery_output_data_dir', '')
        # hack for strange mac behaviour
        if not output_data_dir:
            output_data_dir = ''
        self.output_data_dir_le.setText(output_data_dir)

    def saveState(self):
        """
        Store the options into the user's stored session info.
        """
        mySettings = QSettings()
        mySettings.setValue('irmt/recovery_output_data_dir',
                            self.output_data_dir_le.text())

    def set_ok_button(self):
        self.ok_button.setEnabled(
            os.path.isdir(self.output_data_dir_le.text())
            and self.approach_cbx.currentIndex != -1)
        # and self.dmg_by_asset_layer_cbx.currentIndex != -1
        # and self.svi_layer_cbx.currentIndex != -1
        # and self.svi_field_name_cbx.currentIndex != -1
        # and self.zone_field_name_cbx.currentIndex != -1)

    def populate_approach_cbx(self):
        self.approach_cbx.addItems(['Aggregate', 'Disaggregate'])

    @pyqtSlot(str)
    def on_output_data_dir_le_textChanged(self, text):
        self.output_data_dir = self.output_data_dir_le.text()
        self.set_ok_button()

    @pyqtSlot(str)
    def on_approach_cbx_currentIndexChanged(self, selected_text):
        # TODO: we might need to ask the user to provide the necessary files
        self.set_ok_button()

    @pyqtSlot(int)
    def on_dmg_by_asset_layer_cbx_currentIndexChanged(self, selected_index):
        self.dmg_by_asset_layer = self.dmg_by_asset_layer_cbx.itemData(
            selected_index)

    @pyqtSlot(int)
    def on_svi_layer_cbx_currentIndexChanged(self, selected_index):
        self.svi_layer = self.svi_layer_cbx.itemData(selected_index)
        # FIXME self.svi_field_name is temporarily ignored
        # reload_attrib_cbx(self.svi_field_name_cbx, self.svi_layer)
        reload_attrib_cbx(self.zone_field_name_cbx, self.svi_layer)

    @pyqtSlot()
    def on_output_data_dir_btn_clicked(self):
        path = QFileDialog.getExistingDirectory(
            self, self.tr('Choose output directory'), QDir.homePath())
        if path:
            self.output_data_dir_le.setText(path)

    def read_all_configuration_files(self):
        inspectionTimes = read_config_file('InspectionTimes.txt')
        assessmentTimes = read_config_file('AssessmentTimes.txt')
        mobilizationTimes = read_config_file('MobilizationTimes.txt')
        repairTimes = read_config_file('RepairTimes.txt')
        recoveryTimes = read_config_file('RecoveryTimes.txt')
        leadTimeDispersion = read_config_file(
            'LeadTimeDispersion.txt', float)
        repairTimeDispersion = read_config_file(
            'RepairTimeDispersion.txt', float)
        numberOfDamageSimulations = read_config_file(
            'NumberOfDamageSimulations.txt', int)[0]
        return (inspectionTimes, assessmentTimes, mobilizationTimes,
                repairTimes, recoveryTimes, leadTimeDispersion,
                repairTimeDispersion, numberOfDamageSimulations)

    def generate_community_level_recovery_curve(self, integrate_svi=True):
        # Developed By: <NAME>
        # Edited by: <NAME>
        # Reimplemented for this plugin by: <NAME> and <NAME>
        # Objective: GenerateCommunityLevelRecoveryCurve
        # Initial date: August 26, 2016

        if integrate_svi:
            self.svi_layer = self.svi_layer_cbx.itemData(
                    self.svi_layer_cbx.currentIndex())
            # FIXME self.svi_field_name is temporarily ignored
            # self.svi_field_name = self.svi_field_name_cbx.currentText()
            self.zone_field_name = self.zone_field_name_cbx.currentText()

        start = time.clock()  # FIXME

        # Incorporate Napa Data to community recovery model

        # build dictionary zone_id -> dmg_by_asset_probs
        zonal_dmg_by_asset_probs = defaultdict(list)
        zonal_asset_refs = defaultdict(list)
        if integrate_svi:
            # FIXME self.svi_field_name is temporarily ignored
            # svi_by_zone = dict()
            for zone_feat in self.svi_layer.getFeatures():
                zone_id = str(zone_feat[self.zone_field_name])
                # FIXME self.svi_field_name is temporarily ignored
                # svi_value = zone_feat[self.svi_field_name]
                # svi_by_zone[zone_id] = svi_value
            msg = 'Reading damage state probabilities...'
            msg_bar_item, progress = create_progress_message_bar(
                self.iface.messageBar(), msg)
            tot_features = self.dmg_by_asset_layer.featureCount()
            for feat_idx, dmg_by_asset_feat in enumerate(
                    self.dmg_by_asset_layer.getFeatures(), start=1):
                zone_id = dmg_by_asset_feat[self.zone_field_name]
                # FIXME: hack to handle case in which the zone id is an integer
                # but it is stored as Real
                try:
                    zone_id = str(int(zone_id))
                except:
                    zone_id = str(zone_id)
                # FIXME: same hack as above
                asset_ref = dmg_by_asset_feat['asset_ref']
                try:
                    asset_ref = str(int(asset_ref))
                except:
                    asset_ref = str(asset_ref)
                # select fields that contain probabilities
                # i.e., ignore asset id and taxonomy (first 2 items)
                # and get only columns containing means, discarding
                # those containing stddevs, therefore getting one item
                # out of two for the remaining columns
                # Also discard the last field, containing zone ids
                dmg_by_asset_probs = dmg_by_asset_feat.attributes()[
                    2:-1:2]
                zonal_dmg_by_asset_probs[zone_id].append(dmg_by_asset_probs)
                zonal_asset_refs[zone_id].append(asset_ref)
                progress_perc = feat_idx / float(tot_features) * 100
                progress.setValue(progress_perc)
            clear_progress_message_bar(self.iface.messageBar(), msg_bar_item)
        else:  # ignore svi
            msg = 'Reading damage state probabilities...'
            msg_bar_item, progress = create_progress_message_bar(
                self.iface.messageBar(), msg)
            tot_features = self.dmg_by_asset_layer.featureCount()
            for idx, dmg_by_asset_feat in enumerate(
                    self.dmg_by_asset_layer.getFeatures(), start=1):
                # we don't have any field containing zone ids, to be discarded
                dmg_by_asset_probs = dmg_by_asset_feat.attributes()[2::2]
                asset_ref = dmg_by_asset_feat['asset_ref']
                zonal_dmg_by_asset_probs['ALL'].append(dmg_by_asset_probs)
                zonal_asset_refs['ALL'].append(asset_ref)
                progress_perc = idx / float(tot_features) * 100
                progress.setValue(progress_perc)
            clear_progress_message_bar(self.iface.messageBar(), msg_bar_item)

        tot_zones = len(zonal_dmg_by_asset_probs)
        msg = 'Calculating zone-level recovery curves...'
        msg_bar_item, progress = create_progress_message_bar(
            self.iface.messageBar(), msg)
        # for each zone, calculate a zone-level recovery function
        for idx, zone_id in enumerate(zonal_dmg_by_asset_probs.keys(),
                                      start=1):
            # TODO: use svi_by_zone[zone_id] to adjust recovery times (how?)

            dmg_by_asset_probs = zonal_dmg_by_asset_probs[zone_id]
            asset_refs = zonal_asset_refs[zone_id]

            (LossBasedDamageStateProbabilities,
             RecoveryBasedDamageStateProbabilities,
             fractionCollapsedAndIrreparableBuildings) = \
                self.loss_based_to_recovery_based_probs(dmg_by_asset_probs)

            # FIXME self.svi_field_name is temporarily ignored
            # svi_value = svi_by_zone[zone_id] if integrate_svi else None

            # FIXME If we don't read files again for each zone, time increases
            # across zones. This is not optimal, but configuration files are
            # very small, so reading them is almost instantaneous.
            (inspectionTimes, assessmentTimes, mobilizationTimes, repairTimes,
             recoveryTimes, leadTimeDispersion, repairTimeDispersion,
             numberOfDamageSimulations) = self.read_all_configuration_files()

            # FIXME - when aggregating by zone we are constantly increasing
            # the times with this approach
            (timeList, inspectionTimes,
             assessmentTimes, mobilizationTimes) = self.calculate_times(
                fractionCollapsedAndIrreparableBuildings, inspectionTimes,
                assessmentTimes, mobilizationTimes, repairTimes)
            # assessmentTimes, mobilizationTimes, repairTimes, svi_value)

            # Initialize community recovery function
            communityRecoveryFunction = [0 for x in range(len(timeList))]
            New_communityRecoveryFunction = [
                0 for x in range(len(timeList)+DAYS_BEFORE_EVENT)]

            # Looping over all damage simulations
            for sim in range(numberOfDamageSimulations):
                simulationRecoveryFunction = \
                    self.generate_simulation_recovery_curve(
                        timeList, LossBasedDamageStateProbabilities,
                        RecoveryBasedDamageStateProbabilities, inspectionTimes,
                        recoveryTimes, repairTimes, leadTimeDispersion,
                        repairTimeDispersion, assessmentTimes,
                        mobilizationTimes, zone_id, asset_refs)
                # Sum up all building level recovery function
                # TODO: use enumerate instead
                for timePoint in range(len(timeList)):
                    communityRecoveryFunction[timePoint] \
                        += simulationRecoveryFunction[timePoint]

            # PAOLO: instead of calculating the community level recovery
            # function on all points, we should aggregate points by the same
            # zones defined for the socioeconomic dataset, and then we should
            # produce a community recovery function for each zone.
            # This has to be done on the damage by asset layer
            # (For the aggregation we can use SAGA:
            #  "Add Polygon Attributes to Points", i.e.
            #  processing.runalg('saga:addpolygonattributestopoints', input,
            #                    polygons, field, output))

            # Calculate community level recovery function
            # TODO: use enumerate instead
            for timePoint in range(len(timeList)):
                communityRecoveryFunction[timePoint] /= \
                    numberOfDamageSimulations

            # PAOLO: should we plot this?
            # Plot community level recovery curve
            # plt.plot(timeList, communityRecoveryFunction)
            # plt.show()

            # Plot community level recovery curve which can presents the number
            # of occupants before earthquake
            New_timeList = [x for x in range(len(timeList)+DAYS_BEFORE_EVENT)]
            # TODO: use enumerate instead
            for i in range(len(timeList)+DAYS_BEFORE_EVENT):
                if i < DAYS_BEFORE_EVENT:
                    New_communityRecoveryFunction[i] = 1
                else:
                    New_communityRecoveryFunction[i] = (
                        communityRecoveryFunction[i - DAYS_BEFORE_EVENT]
                        / len(LossBasedDamageStateProbabilities))

            fig = plt.figure()
            # highlight values at observation days (after 6, 12 and 18 months)
            obs_days = [DAYS_BEFORE_EVENT + day for day in (180, 360, 540)]
            plt.plot(New_timeList, New_communityRecoveryFunction)
            for obs_day in obs_days:
                plt.axvline(x=obs_day, linestyle='dotted')
                i = obs_day
                j = New_communityRecoveryFunction[obs_day]
                plt.annotate(
                    '%.3f' % j,
                    xy=(i, j),
                    xytext=(0, 30),
                    textcoords='offset points')
            plt.xlabel('Time (days)')
            plt.ylabel('Normalized recovery level')
            plt.title('Community level recovery curve for zone %s' % zone_id)
            plt.ylim((0.0, 1.2))
            # plt.show()
            filestem = os.path.join(
                self.output_data_dir, "recovery_function_zone_%s" % zone_id)
            fig.savefig(filestem + '.png')

            # Save community recovery function
            f3 = open(filestem + '.txt', "w")
            f3.write(str(New_communityRecoveryFunction))
            f3.close()

            end = time.clock()
            print (end - start)
            progress_perc = idx / float(tot_zones) * 100
            progress.setValue(progress_perc)

        clear_progress_message_bar(self.iface.messageBar(), msg_bar_item)

    def generate_simulation_recovery_curve(
            self, timeList, LossBasedDamageStateProbabilities,
            RecoveryBasedDamageStateProbabilities, inspectionTimes,
            recoveryTimes, repairTimes, leadTimeDispersion,
            repairTimeDispersion, assessmentTimes, mobilizationTimes,<commit_msg>Save a csv with recovery predictions at observed days<commit_after># -*- coding: utf-8 -*-
# /***************************************************************************
# Irmt
#                                 A QGIS plugin
# OpenQuake Integrated Risk Modelling Toolkit
#                              -------------------
#        begin                : 2013-10-24
#        copyright            : (C) 2014 by GEM Foundation
#        email                : <EMAIL>
# ***************************************************************************/
#
# OpenQuake is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# OpenQuake is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake.  If not, see <http://www.gnu.org/licenses/>.

import csv
import os
import time
import matplotlib.pyplot as plt
from collections import defaultdict

from PyQt4.QtCore import pyqtSlot, QSettings, QDir
from PyQt4.QtGui import (QDialog,
                         QFileDialog,
                         QDialogButtonBox)
from qgis.core import QgsMapLayer
from qgis.gui import QgsMessageBar
from svir.calculations.aggregate_loss_by_zone import add_zone_id_to_points

from svir.utilities.utils import (get_ui_class,
                                  read_config_file,
                                  reload_layers_in_cbx,
                                  reload_attrib_cbx,
                                  WaitCursorManager,
                                  tr,
                                  create_progress_message_bar,
                                  clear_progress_message_bar,
                                  )
from svir.recovery_modeling.building import Building

FORM_CLASS = get_ui_class('ui_recovery_modeling.ui')

DAYS_BEFORE_EVENT = 200
SVI_WEIGHT_COEFF = 1  # FIXME: Let the user set this parameter


class RecoveryModelingDialog(QDialog, FORM_CLASS):
    """
    Modal dialog giving to perform recovery modeling analysis.

    Object-oriented Programming(OOP) is used to model post-earthquake community
    recovery of residential community.

    PAOLO: FIXME In the last version, we only have a Building class
    In this Programming, four classes are
    adopted, namely Building Class, Ground Motion FIelds Class, Constructuon &
    Enigneering Service Class and Community Class.

    In each class, a set of data and functions work together to complete a task
    which we call attributes and functions.  Attribute is used to specify if
    and how object properties are accessed within the programming environment
    Function is used to create/altere to achieve task or functionality and
    return the attributes of an object

    Procedure of this OOP Programming:
    Step 1: Define attributes of objects in each class
    Step 2: Set up functions in each class
    Step 3: Incorporate Napa Data to community recovery model

    Methodology:
    Time-based method is utilized which characterize a probability density
    function of the time it takes to a higher or lower functioning state given
    a set of explanatory variables such as the extent of damage to the
    building.
    """
    def __init__(self, iface):
        QDialog.__init__(self)
        self.iface = iface
        # Set up the user interface from Designer.
        self.setupUi(self)
        self.ok_button = self.buttonBox.button(QDialogButtonBox.Ok)
        self.populate_approach_cbx()
        reload_layers_in_cbx(
            self.dmg_by_asset_layer_cbx, [QgsMapLayer.VectorLayer])
        reload_layers_in_cbx(self.svi_layer_cbx, [QgsMapLayer.VectorLayer])
        recovery_modeling_dir = os.path.join(
            os.path.dirname(os.path.abspath(__file__)),
            '..', 'recovery_modeling')
        self.input_data_dir = os.path.join(recovery_modeling_dir, 'input_data')
        self.restoreState()
        self.set_ok_button()

    def restoreState(self):
        """
        Reinstate the options based on the user's stored session info.
        """
        mySettings = QSettings()
        output_data_dir = mySettings.value('irmt/recovery_output_data_dir', '')
        # hack for strange mac behaviour
        if not output_data_dir:
            output_data_dir = ''
        self.output_data_dir_le.setText(output_data_dir)

    def saveState(self):
        """
        Store the options into the user's stored session info.
        """
        mySettings = QSettings()
        mySettings.setValue('irmt/recovery_output_data_dir',
                            self.output_data_dir_le.text())

    def set_ok_button(self):
        self.ok_button.setEnabled(
            os.path.isdir(self.output_data_dir_le.text())
            and self.approach_cbx.currentIndex != -1)
        # and self.dmg_by_asset_layer_cbx.currentIndex != -1
        # and self.svi_layer_cbx.currentIndex != -1
        # and self.svi_field_name_cbx.currentIndex != -1
        # and self.zone_field_name_cbx.currentIndex != -1)

    def populate_approach_cbx(self):
        self.approach_cbx.addItems(['Aggregate', 'Disaggregate'])

    @pyqtSlot(str)
    def on_output_data_dir_le_textChanged(self, text):
        self.output_data_dir = self.output_data_dir_le.text()
        self.set_ok_button()

    @pyqtSlot(str)
    def on_approach_cbx_currentIndexChanged(self, selected_text):
        # TODO: we might need to ask the user to provide the necessary files
        self.set_ok_button()

    @pyqtSlot(int)
    def on_dmg_by_asset_layer_cbx_currentIndexChanged(self, selected_index):
        self.dmg_by_asset_layer = self.dmg_by_asset_layer_cbx.itemData(
            selected_index)

    @pyqtSlot(int)
    def on_svi_layer_cbx_currentIndexChanged(self, selected_index):
        self.svi_layer = self.svi_layer_cbx.itemData(selected_index)
        # FIXME self.svi_field_name is temporarily ignored
        # reload_attrib_cbx(self.svi_field_name_cbx, self.svi_layer)
        reload_attrib_cbx(self.zone_field_name_cbx, self.svi_layer)

    @pyqtSlot()
    def on_output_data_dir_btn_clicked(self):
        path = QFileDialog.getExistingDirectory(
            self, self.tr('Choose output directory'), QDir.homePath())
        if path:
            self.output_data_dir_le.setText(path)

    def read_all_configuration_files(self):
        inspectionTimes = read_config_file('InspectionTimes.txt')
        assessmentTimes = read_config_file('AssessmentTimes.txt')
        mobilizationTimes = read_config_file('MobilizationTimes.txt')
        repairTimes = read_config_file('RepairTimes.txt')
        recoveryTimes = read_config_file('RecoveryTimes.txt')
        leadTimeDispersion = read_config_file(
            'LeadTimeDispersion.txt', float)
        repairTimeDispersion = read_config_file(
            'RepairTimeDispersion.txt', float)
        numberOfDamageSimulations = read_config_file(
            'NumberOfDamageSimulations.txt', int)[0]
        return (inspectionTimes, assessmentTimes, mobilizationTimes,
                repairTimes, recoveryTimes, leadTimeDispersion,
                repairTimeDispersion, numberOfDamageSimulations)

    def generate_community_level_recovery_curve(self, integrate_svi=True):
        # Developed By: <NAME>
        # Edited by: <NAME>
        # Reimplemented for this plugin by: <NAME> and <NAME>
        # Objective: GenerateCommunityLevelRecoveryCurve
        # Initial date: August 26, 2016

        if integrate_svi:
            self.svi_layer = self.svi_layer_cbx.itemData(
                    self.svi_layer_cbx.currentIndex())
            # FIXME self.svi_field_name is temporarily ignored
            # self.svi_field_name = self.svi_field_name_cbx.currentText()
            self.zone_field_name = self.zone_field_name_cbx.currentText()

        start = time.clock()  # FIXME

        # Incorporate Napa Data to community recovery model

        # build dictionary zone_id -> dmg_by_asset_probs
        zonal_dmg_by_asset_probs = defaultdict(list)
        zonal_asset_refs = defaultdict(list)
        if integrate_svi:
            # FIXME self.svi_field_name is temporarily ignored
            # svi_by_zone = dict()
            for zone_feat in self.svi_layer.getFeatures():
                zone_id = str(zone_feat[self.zone_field_name])
                # FIXME self.svi_field_name is temporarily ignored
                # svi_value = zone_feat[self.svi_field_name]
                # svi_by_zone[zone_id] = svi_value
            msg = 'Reading damage state probabilities...'
            msg_bar_item, progress = create_progress_message_bar(
                self.iface.messageBar(), msg)
            tot_features = self.dmg_by_asset_layer.featureCount()
            for feat_idx, dmg_by_asset_feat in enumerate(
                    self.dmg_by_asset_layer.getFeatures(), start=1):
                zone_id = dmg_by_asset_feat[self.zone_field_name]
                # FIXME: hack to handle case in which the zone id is an integer
                # but it is stored as Real
                try:
                    zone_id = str(int(zone_id))
                except:
                    zone_id = str(zone_id)
                # FIXME: same hack as above
                asset_ref = dmg_by_asset_feat['asset_ref']
                try:
                    asset_ref = str(int(asset_ref))
                except:
                    asset_ref = str(asset_ref)
                # select fields that contain probabilities
                # i.e., ignore asset id and taxonomy (first 2 items)
                # and get only columns containing means, discarding
                # those containing stddevs, therefore getting one item
                # out of two for the remaining columns
                # Also discard the last field, containing zone ids
                dmg_by_asset_probs = dmg_by_asset_feat.attributes()[
                    2:-1:2]
                zonal_dmg_by_asset_probs[zone_id].append(dmg_by_asset_probs)
                zonal_asset_refs[zone_id].append(asset_ref)
                progress_perc = feat_idx / float(tot_features) * 100
                progress.setValue(progress_perc)
            clear_progress_message_bar(self.iface.messageBar(), msg_bar_item)
        else:  # ignore svi
            msg = 'Reading damage state probabilities...'
            msg_bar_item, progress = create_progress_message_bar(
                self.iface.messageBar(), msg)
            tot_features = self.dmg_by_asset_layer.featureCount()
            for idx, dmg_by_asset_feat in enumerate(
                    self.dmg_by_asset_layer.getFeatures(), start=1):
                # we don't have any field containing zone ids, to be discarded
                dmg_by_asset_probs = dmg_by_asset_feat.attributes()[2::2]
                asset_ref = dmg_by_asset_feat['asset_ref']
                zonal_dmg_by_asset_probs['ALL'].append(dmg_by_asset_probs)
                zonal_asset_refs['ALL'].append(asset_ref)
                progress_perc = idx / float(tot_features) * 100
                progress.setValue(progress_perc)
            clear_progress_message_bar(self.iface.messageBar(), msg_bar_item)

        tot_zones = len(zonal_dmg_by_asset_probs)
        msg = 'Calculating zone-level recovery curves...'
        msg_bar_item, progress = create_progress_message_bar(
            self.iface.messageBar(), msg)
        at_observed_days_filename = os.path.join(
            self.output_data_dir, 'at_observed_days.csv')
        at_observed_days = open(at_observed_days_filename, 'w')
        writer = csv.writer(at_observed_days)
        header = ['zone_id', '6months', '12months', '18months']
        writer.writerow(header)
        # for each zone, calculate a zone-level recovery function
        for idx, zone_id in enumerate(zonal_dmg_by_asset_probs.keys(),
                                      start=1):
            # TODO: use svi_by_zone[zone_id] to adjust recovery times (how?)

            dmg_by_asset_probs = zonal_dmg_by_asset_probs[zone_id]
            asset_refs = zonal_asset_refs[zone_id]

            (LossBasedDamageStateProbabilities,
             RecoveryBasedDamageStateProbabilities,
             fractionCollapsedAndIrreparableBuildings) = \
                self.loss_based_to_recovery_based_probs(dmg_by_asset_probs)

            # FIXME self.svi_field_name is temporarily ignored
            # svi_value = svi_by_zone[zone_id] if integrate_svi else None

            # FIXME If we don't read files again for each zone, time increases
            # across zones. This is not optimal, but configuration files are
            # very small, so reading them is almost instantaneous.
            (inspectionTimes, assessmentTimes, mobilizationTimes, repairTimes,
             recoveryTimes, leadTimeDispersion, repairTimeDispersion,
             numberOfDamageSimulations) = self.read_all_configuration_files()

            # FIXME - when aggregating by zone we are constantly increasing
            # the times with this approach
            (timeList, inspectionTimes,
             assessmentTimes, mobilizationTimes) = self.calculate_times(
                fractionCollapsedAndIrreparableBuildings, inspectionTimes,
                assessmentTimes, mobilizationTimes, repairTimes)
            # assessmentTimes, mobilizationTimes, repairTimes, svi_value)

            # Initialize community recovery function
            communityRecoveryFunction = [0 for x in range(len(timeList))]
            New_communityRecoveryFunction = [
                0 for x in range(len(timeList)+DAYS_BEFORE_EVENT)]

            # Looping over all damage simulations
            for sim in range(numberOfDamageSimulations):
                simulationRecoveryFunction = \
                    self.generate_simulation_recovery_curve(
                        timeList, LossBasedDamageStateProbabilities,
                        RecoveryBasedDamageStateProbabilities, inspectionTimes,
                        recoveryTimes, repairTimes, leadTimeDispersion,
                        repairTimeDispersion, assessmentTimes,
                        mobilizationTimes, zone_id, asset_refs)
                # Sum up all building level recovery function
                # TODO: use enumerate instead
                for timePoint in range(len(timeList)):
                    communityRecoveryFunction[timePoint] \
                        += simulationRecoveryFunction[timePoint]

            # PAOLO: instead of calculating the community level recovery
            # function on all points, we should aggregate points by the same
            # zones defined for the socioeconomic dataset, and then we should
            # produce a community recovery function for each zone.
            # This has to be done on the damage by asset layer
            # (For the aggregation we can use SAGA:
            #  "Add Polygon Attributes to Points", i.e.
            #  processing.runalg('saga:addpolygonattributestopoints', input,
            #                    polygons, field, output))

            # Calculate community level recovery function
            # TODO: use enumerate instead
            for timePoint in range(len(timeList)):
                communityRecoveryFunction[timePoint] /= \
                    numberOfDamageSimulations

            # PAOLO: should we plot this?
            # Plot community level recovery curve
            # plt.plot(timeList, communityRecoveryFunction)
            # plt.show()

            # Plot community level recovery curve which can presents the number
            # of occupants before earthquake
            New_timeList = [x for x in range(len(timeList)+DAYS_BEFORE_EVENT)]
            # TODO: use enumerate instead
            for i in range(len(timeList)+DAYS_BEFORE_EVENT):
                if i < DAYS_BEFORE_EVENT:
                    New_communityRecoveryFunction[i] = 1
                else:
                    New_communityRecoveryFunction[i] = (
                        communityRecoveryFunction[i - DAYS_BEFORE_EVENT]
                        / len(LossBasedDamageStateProbabilities))

            fig = plt.figure()
            # highlight values at observation days (after 6, 12 and 18 months)
            obs_days = [DAYS_BEFORE_EVENT + day for day in (180, 360, 540)]
            plt.plot(New_timeList, New_communityRecoveryFunction)
            row = [zone_id]
            for obs_day in obs_days:
                value_at_obs_day = New_communityRecoveryFunction[obs_day]
                row.append(value_at_obs_day)
                plt.axvline(x=obs_day, linestyle='dotted')
                i = obs_day
                j = value_at_obs_day
                plt.annotate(
                    '%.3f' % j,
                    xy=(i, j),
                    xytext=(0, 30),
                    textcoords='offset points')
            writer.writerow(row)
            plt.xlabel('Time (days)')
            plt.ylabel('Normalized recovery level')
            plt.title('Community level recovery curve for zone %s' % zone_id)
            plt.ylim((0.0, 1.2))
            # plt.show()
            filestem = os.path.join(
                self.output_data_dir, "recovery_function_zone_%s" % zone_id)
            fig.savefig(filestem + '.png')

            # Save community recovery function
            f3 = open(filestem + '.txt', "w")
            f3.write(str(New_communityRecoveryFunction))
            f3.close()

            end = time.clock()
            print (end - start)
            progress_perc = idx / float(tot_zones) * 100
            progress.setValue(progress_perc)

        clear_progress_message_bar(self.iface.messageBar(), msg_bar_item)
        at_observed_days.close()

    def generate_simulation_recovery_curve(
            self, timeList, LossBasedDamageStateProbabilities,
            RecoveryBasedDamageStateProbabilities, inspectionTimes,
            recoveryTimes, repairTimes, leadTimeDispersion,
            repairTimeDispersion, assessmentTimes, mobilizationTimes,
{'fraction_non_alphanumeric': 0.07570593288441821, 'fraction_numerical': 0.0046277052678249785, 'mean_word_length': 2.2192685102586975, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2637038', 'n_tokens_mistral': 10121, 'n_tokens_neox': 9647, 'n_words': 2433}