text: # -*- coding: utf-8 -*-
  #
  # Copyright (c) 2018 Leland Stanford Junior University
  # Copyright (c) 2018 The Regents of the University of California
  #
  # This file is part of pelicun.
  # 
  # Redistribution and use in source and binary forms, with or without 
  # modification, are permitted provided that the following conditions are met:
  #
  # 1. Redistributions of source code must retain the above copyright notice, 
  # this list of conditions and the following disclaimer.
  #
  # 2. Redistributions in binary form must reproduce the above copyright notice, 
  # this list of conditions and the following disclaimer in the documentation 
  # and/or other materials provided with the distribution.
  #
  # 3. Neither the name of the copyright holder nor the names of its contributors 
  # may be used to endorse or promote products derived from this software without 
  # specific prior written permission.
  #
  # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
  # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
  # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
  # ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE 
  # LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
  # CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
  # SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
  # INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
  # CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
  # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
  # POSSIBILITY OF SUCH DAMAGE.
  # 
  # You should have received a copy of the BSD 3-Clause License along with 
  # pelicun. If not, see <http://www.opensource.org/licenses/>.
  #
  # Contributors:
  # Adam Zsarnóczay
  
  """
  This module has classes and methods that control the loss assessment.
  
  .. rubric:: Contents
  
  .. autosummary::
  
      Assessment
      FEMA_P58_Assessment
  
  """
  
  from .base import *
  from .uq import *
  from .model import *
  from .file_io import *
  
  class Assessment(object):
      """
      A high-level class that collects features common to all supported loss
      assessment methods. This class will only rarely be called directly when
      using pelicun.
      """
      
      def __init__(self):
          
          # initialize the basic data containers      
          # inputs
          self._AIM_in = None
          self._EDP_in = None
          self._POP_in = None
          self._FG_in = None
          
          # random variables and loss model
          self._RV_dict = None # dictionary to store random variables
          self._EDP_dict = None
          self._FG_dict = None
          
          # results
          self._TIME = None
          self._POP = None
          self._COL = None
          self._ID_dict = None
          self._DMG = None
          self._DV_dict = None
          self._SUMMARY = None
          
      def read_inputs(self, path_DL_input, path_EDP_input, verbose=False):
          """
          Read and process the input files to describe the loss assessment task.
          
          Parameters
          ----------
          path_DL_input: string
              Location of the Damage and Loss input file. The file is expected to
              be a JSON with data stored in a standard format described in detail
              in the Input section of the documentation.
          path_EDP_input: string
              Location of the EDP input file. The file is expected to follow the 
              output formatting of Dakota. The Input section of the documentation
              provides more information about the expected formatting.
          verbose: boolean, default: False
              If True, the method echoes the information read from the files.
              This can be useful to ensure that the information in the file is 
              properly read by the method.
          """
          
          # read SimCenter inputs -----------------------------------------------
          # BIM file
          self._AIM_in = read_SimCenter_DL_input(path_DL_input, verbose=verbose)
  
          # EDP file
          self._EDP_in = read_SimCenter_EDP_input(path_EDP_input, 
                                                  verbose=verbose)
      
      def define_random_variables(self):
          """
          Define the random variables used for loss assessment.
  
          """
          pass
      
      def define_loss_model(self):
          """
          Create the stochastic loss model based on the inputs provided earlier.
  
          """
          pass
      
      def calculate_damage(self):
          """
          Characterize the damage experienced in each random event realization.
  
          """
          self._ID_dict = {}
      
      def calculate_losses(self):
          """
          Characterize the consequences of damage in each random event realization.
  
          """
          self._DV_dict = {}
      
      def write_outputs(self):
          """
          Export the results.
  
          """
          pass
      
  class FEMA_P58_Assessment(Assessment):
      """
      An Assessment class that implements the loss assessment method in FEMA P58.
      """
      def __init__(self, inj_lvls = 2):
          super(FEMA_P58_Assessment, self).__init__()
          
          # constants for the FEMA-P58 methodology
          self._inj_lvls = inj_lvls
          
      @property
      def beta_additional(self):
          """
          Calculate the total additional uncertainty for post processing.
          
          The total additional uncertainty is the squared root of sum of squared
          uncertainties corresponding to ground motion and modeling.
          
          Returns
          -------
          beta_tot: float
              The total uncertainty (logarithmic EDP standard deviation) to add
              to the EDP distribution.
          """
          
          AU = self._AIM_in['general']['added_uncertainty']
          
          beta_tot = 0.
          if AU['beta_m'] is not None:
              beta_tot += AU['beta_m']**2.
          if AU['beta_gm'] is not None:
              beta_tot += AU['beta_gm']**2.
              
          # if no uncertainty is assigned, we consider a minimum of 10^-4
          if beta_tot == 0:
              beta_tot = 1e-4
          else:
              beta_tot = np.sqrt(beta_tot)
          
          return beta_tot
  
      def read_inputs(self, path_DL_input, path_EDP_input, 
                      path_CMP_data=None, path_POP_data=None, verbose=False):
          """
          Read and process the input files to describe the loss assessment task.
          
          Parameters
          ----------
          path_DL_input: string
              Location of the Damage and Loss input file. The file is expected to
              be a JSON with data stored in a standard format described in detail
              in the Input section of the documentation.
          path_EDP_input: string
              Location of the EDP input file. The file is expected to follow the 
              output formatting of Dakota. The Input section of the documentation
              provides more information about the expected formatting.
          path_CMP_data: string, default: None
              Location of the folder with component damage and loss data files. 
              The default None value triggers the use of the FEMA P58 first
              edition data from pelicun/resources/component DL/FEMA P58 first 
              edition/.
          path_POP_data: string, default: None
              Location of the JSON file that describes the temporal distribution
              of the population per the FEMA P58 method. The default None value
              triggers the use of the FEMA P58 first edition distribution from
              pelicun/resources/population.json. 
          verbose: boolean, default: False
              If True, the method echoes the information read from the files.
              This can be useful to ensure that the information in the file is 
              properly read by the method.
  
          """
          
          super(FEMA_P58_Assessment, self).read_inputs(path_DL_input, path_EDP_input, verbose)
          
          # check if the component data path is provided by the user
          if path_CMP_data is None:
              warnings.warn(UserWarning(
                  "The component database is not specified; using the default "
                  "FEMA P58 first edition data."
              ))
              #path_CMP_data = '../../resources/component DL/FEMA P58 first edition/'
              path_CMP_data = pelicun_path
              path_CMP_data += '/resources/component DL/FEMA P58 first edition/'
  
          if path_POP_data is None:
              warnings.warn(UserWarning(
                  "The population distribution is not specified; using the default "
                  "FEMA P58 first edition data."
              ))
              #path_POP_data = '../../resources/population.json'
              path_POP_data = pelicun_path
              path_POP_data += '/resources/population.json'
  
          # assume that the asset is a building
          # TODO: If we want to apply FEMA-P58 to non-building assets, several parts of this methodology need to be extended.
          BIM = self._AIM_in
  
          # read component and population data ----------------------------------
          # components
          self._FG_in = read_component_DL_data(path_CMP_data, BIM['components'],
                                               verbose=verbose)
          
          # population
          POP = read_population_distribution(
              path_POP_data, 
              BIM['general']['occupancy_type'], 
              verbose=verbose)
  
          POP['peak'] = BIM['general']['population']          
          self._POP_in = POP              
  
      def define_random_variables(self):
          """
          Define the random variables used for loss assessment.
          
          Following the FEMA P58 methodology, the groups of parameters below are 
          considered random. Simple correlation structures within each group can
          be specified through the DL input file. The random decision variables
          are only created and used later if those particular decision variables
          are requested in the input file.
          
          1. Demand (EDP) distribution
          
          Describe the uncertainty in the demands. Unlike other random variables,
          the EDPs are characterized by the EDP input data provided earlier. All
          EDPs are handled in one multivariate lognormal distribution. If more 
          than one sample is provided, the distribution is fit to the EDP data.
          Otherwise, the provided data point is assumed to be the median value
          and the additional uncertainty prescribed describes the dispersion. See
          _create_RV_demands() for more details.
          
          2. Component quantities
          
          Describe the uncertainty in the quantity of components in each 
          Performance Group. All Fragility Groups are handled in the same 
          multivariate distribution. Consequently, correlation between various 
          groups of component quantities can be specified. See 
          _create_RV_quantities() for details. 
          
          3. Fragility EDP limits
          
          Describe the uncertainty in the EDP limit that corresponds to 
          exceedance of each Damage State. EDP limits are grouped by Fragility
          Groups. Consequently, correlation between fragility limits are 
          currently are limited within Fragility Groups. See 
          _create_RV_fragilities() for details.
          
          4. Reconstruction cost and time
          
          Describe the uncertainty in the cost and duration of reconstruction of 
          each component conditioned on the damage state of the component. All
          Fragility Groups are handled in the same multivariate distribution.
          Consequently, correlation between various groups of component 
          reconstruction time and cost estimates can be specified. See 
          _create_RV_repairs() for details.
          
          5. Damaged component proportions that trigger a red tag
          
          Describe the uncertainty in the amount of damaged components needed to
          trigger a red tag for the building. All Fragility Groups are handled in
          the same multivariate distribution. Consequently, correlation between
          various groups of component proportion limits can be specified. See
          _create_RV_red_tags() for details.
          
          6. Injuries
          
          Describe the uncertainty in the proportion of people in the affected
          area getting injuries exceeding a certain level of severity. FEMA P58
          uses two severity levels: injury and fatality. Both levels for all
          Fragility Groups are handled in the same multivariate distribution. 
          Consequently, correlation between various groups of component injury
          expectations can be specified. See _create_RV_injuries() for details. 
          """
          super(FEMA_P58_Assessment, self).define_random_variables()
  
          # create the random variables -----------------------------------------
          DEP = self._AIM_in['dependencies']
          
          self._RV_dict = {}
  
          # quantities 100
          self._RV_dict.update({'QNT': 
                              self._create_RV_quantities(DEP['quantities'])})
          
          # fragilities 300
          s_fg_keys = sorted(self._FG_in.keys())
          for c_id, c_name in enumerate(s_fg_keys):
              comp = self._FG_in[c_name]
              
              self._RV_dict.update({
                  'FR-' + c_name: 
                      self._create_RV_fragilities(c_id, comp, 
                                                  DEP['fragilities'])})
              
          # consequences 400
          DVs = self._AIM_in['decision_variables']
          
          if DVs['red_tag']:
              self._RV_dict.update({'DV_RED': 
                                  self._create_RV_red_tags(DEP['red_tags'])})
          if DVs['rec_time'] or DVs['rec_cost']:
              self._RV_dict.update({'DV_REP': 
                                  self._create_RV_repairs(
                                      DEP['rec_costs'], 
                                      DEP['rec_times'],
                                      DEP['cost_and_time'])})
          if DVs['injuries']:
              self._RV_dict.update({'DV_INJ': 
                                  self._create_RV_injuries(
                                      DEP['injuries'],
                                      DEP['injury_lvls'])})
  
          # demands 200
          self._RV_dict.update({'EDP': self._create_RV_demands(
              self.beta_additional)})
  
          # sample the random variables -----------------------------------------
          s_rv_keys = sorted(self._RV_dict.keys())
          for r_i in s_rv_keys:
              rv = self._RV_dict[r_i]
              if rv is not None:
                  rv.sample_distribution(self._AIM_in['general']['realizations'])
  
      def define_loss_model(self):
          """
          Create the stochastic loss model based on the inputs provided earlier.
          
          Following the FEMA P58 methodology, the components specified in the 
          Damage and Loss input file are used to create Fragility Groups. Each
          Fragility Group corresponds to a component that might be present in 
          the building at several locations. See _create_fragility_groups() for
          more details about the creation of Fragility Groups.
  
          """
          super(FEMA_P58_Assessment, self).define_loss_model()
          
          # fragility groups
          self._FG_dict = self._create_fragility_groups()
  
          # demands
          self._EDP_dict = dict(
              [(tag, RandomVariableSubset(self._RV_dict['EDP'],tags=tag))
               for tag in self._RV_dict['EDP']._dimension_tags])
          
      def calculate_damage(self):
          """
          Characterize the damage experienced in each random event realization.
          
          First, the time of the event (month, weekday/weekend, hour) is randomly
          generated for each realization. Given the event time, the number of 
          people present at each floor of the building is calculated.
          
          Second, the realizations that led to collapse are filtered. See 
          _calc_collapses() for more details on collapse estimation.
          
          Finally, the realizations that did not lead to building collapse are 
          further investigated and the quantities of components in each damage
          state are estimated. See _calc_damage() for more details on damage 
          estimation.
  
          """
          super(FEMA_P58_Assessment, self).calculate_damage()
  
          # event time - month, weekday, and hour realizations
          self._TIME = self._sample_event_time()
          
          # get the population conditioned on event time
          self._POP = self._get_population()
          
          # collapses
          self._COL, collapsed_IDs = self._calc_collapses()
          self._ID_dict.update({'collapse':collapsed_IDs})
  
          # select the non-collapse cases for further analyses
          non_collapsed_IDs = self._POP[
              ~self._POP.index.isin(collapsed_IDs)].index.values.astype(int)
          self._ID_dict.update({'non-collapse': non_collapsed_IDs})
          
          # damage in non-collapses
          self._DMG = self._calc_damage()
  
      def calculate_losses(self):
          """
          Characterize the consequences of damage in each random event realization.
          
          For the sake of efficiency, only the decision variables requested in 
          the input file are estimated. The following consequences are handled by 
          this method:
          
          Reconstruction time and cost
          Estimate the irrepairable cases based on residual drift magnitude and
          the provided irrepairable drift limits. Realizations that led to 
          irrepairable damage or collapse are assigned the replacement cost and
          time of the building when reconstruction cost and time is estimated. 
          Repairable cases get a cost and time estimate for each Damage State in 
          each Performance Group. For more information about estimating 
          irrepairability see _calc_irrepairable() and reconstruction cost and
          time see _calc_repair_cost_and_time() methods.
          
          Injuries
          Collapse-induced injuries are based on the collapse modes and 
          corresponding injury characterization. Injuries conditioned on no 
          collapse are based on the affected area and the probability of 
          injuries of various severity specified in the component data file. For
          more information about estimating injuries conditioned on collapse and
          no collapse, see _calc_collapse_injuries() and 
          _calc_non_collapse_injuries, respecitvely.
          
          Red Tag
          The probability of getting an unsafe placard or red tag is a function
          of the amount of damage experienced in various Damage States for each
          Performance Group. The damage limits that trigger an unsafe placard are
          specified in the component data file. For more information on 
          assigning red tags to realizations see the _calc_red_tag() method.
          
          """
          super(FEMA_P58_Assessment, self).calculate_losses()
          DVs = self._AIM_in['decision_variables']
          
          # red tag probability
          if DVs['red_tag']:
              DV_RED = self._calc_red_tag()
              
              self._DV_dict.update({'red_tag': DV_RED})
  
          # reconstruction cost and time
          if DVs['rec_cost'] or DVs['rec_time']:
              # irrepairable cases
              irrepairable_IDs = self._calc_irrepairable()
      
              # collect the IDs of repairable realizations
              P_NC = self._POP.loc[self._ID_dict['non-collapse']]
              repairable_IDs = P_NC[
                  ~P_NC.index.isin(irrepairable_IDs)].index.values.astype(int)
      
              self._ID_dict.update({'repairable': repairable_IDs})
              self._ID_dict.update({'irrepairable': irrepairable_IDs})
      
              # reconstruction cost and time for repairable cases
              DV_COST, DV_TIME = self._calc_repair_cost_and_time()
  
              self._DV_dict.update({
                  'rec_cost': DV_COST,
                  'rec_time': DV_TIME,
              })
              
          # injuries due to collapse
          if DVs['injuries']:
              COL_INJ = self._calc_collapse_injuries()
      
              # injuries in non-collapsed cases
              DV_INJ_dict = self._calc_non_collapse_injuries()
              
              # store results
              if COL_INJ is not None:
                  self._COL = pd.concat([self._COL, COL_INJ], axis=1)
  
              self._DV_dict.update({'injuries': DV_INJ_dict})
          
      def aggregate_results(self):
          """
          
          Returns
          -------
  
          """
  
          DVs = self._AIM_in['decision_variables']
  
          MI_raw = [
              ('event time', 'month'),
              ('event time', 'weekday?'),
              ('event time', 'hour'),
              ('inhabitants', ''),
              ('collapses', 'collapsed?'),
              ('collapses', 'mode'),
              ('red tagged?', ''),
              ('reconstruction', 'irrepairable?'),
              ('reconstruction', 'cost impractical?'),
              ('reconstruction', 'cost'),
              ('reconstruction', 'time impractical?'),
              ('reconstruction', 'time-sequential'),
              ('reconstruction', 'time-parallel'),
              ('injuries', 'casualties'),
              ('injuries', 'fatalities'),
          ]
          
          ncID = self._ID_dict['non-collapse']
          colID = self._ID_dict['collapse']
          if DVs['rec_cost'] or DVs['rec_time']:
              repID = self._ID_dict['repairable']
              irID = self._ID_dict['irrepairable']
  
          MI = pd.MultiIndex.from_tuples(MI_raw)
  
          SUMMARY = pd.DataFrame(np.empty((
              self._AIM_in['general']['realizations'], 
              len(MI))), columns=MI)
          SUMMARY[:] = np.NaN
  
          # event time
          for prop in ['month', 'weekday?', 'hour']:
              offset = 0
              if prop == 'month':
                  offset = 1
              SUMMARY.loc[:, ('event time', prop)] = \
                  self._TIME.loc[:, prop] + offset
  
          # inhabitants
          SUMMARY.loc[:, ('inhabitants', '')] = self._POP.sum(axis=1)
  
          # collapses
          SUMMARY.loc[:, ('collapses', 'collapsed?')] = self._COL.iloc[:, 0]        
  
          # red tag
          if DVs['red_tag']:
              SUMMARY.loc[ncID, ('red tagged?', '')] = \
                  self._DV_dict['red_tag'].max(axis=1)
  
          # reconstruction cost
          if DVs['rec_cost']:
              SUMMARY.loc[ncID, ('reconstruction', 'cost')] = \
                  self._DV_dict['rec_cost'].sum(axis=1)
      
              repl_cost = self._AIM_in['general']['replacement_cost']
              SUMMARY.loc[colID, ('reconstruction', 'cost')] = repl_cost
  
          if DVs['rec_cost'] or DVs['rec_time']:
              SUMMARY.loc[ncID, ('reconstruction', 'irrepairable?')] = 0
              SUMMARY.loc[irID, ('reconstruction', 'irrepairable?')] = 1
  
          if DVs['rec_cost']:
              SUMMARY.loc[irID, ('reconstruction', 'cost')] = repl_cost
  
          
              repair_impractical_IDs = SUMMARY.loc[
                  SUMMARY.loc[:, ('reconstruction', 'cost')] > repl_cost].index
              SUMMARY.loc[repID, ('reconstruction', 'cost impractical?')] = 0
              SUMMARY.loc[repair_impractical_IDs,
                          ('reconstruction', 'cost impractical?')] = 1
              SUMMARY.loc[
                  repair_impractical_IDs, ('reconstruction', 'cost')] = repl_cost
  
          # reconstruction time
          if DVs['rec_time']:
              SUMMARY.loc[ncID, ('reconstruction', 'time-sequential')] = \
                  self._DV_dict['rec_time'].sum(axis=1)
              SUMMARY.loc[ncID, ('reconstruction', 'time-parallel')] = \
                  self._DV_dict['rec_time'].max(axis=1)
  
              rep_time = self._AIM_in['general']['replacement_time']
  
              for t_label in ['time-sequential', 'time-parallel']:
                  SUMMARY.loc[colID, ('reconstruction', t_label)] = rep_time
                  SUMMARY.loc[irID, ('reconstruction', t_label)] = rep_time
  
              repair_impractical_IDs = \
                  SUMMARY.loc[SUMMARY.loc[:, ('reconstruction',
                                              'time-parallel')] > rep_time].index
              SUMMARY.loc[repID, ('reconstruction', 'time impractical?')] = 0
              SUMMARY.loc[repair_impractical_IDs,('reconstruction', 
                                                  'time impractical?')] = 1
              SUMMARY.loc[repair_impractical_IDs, ('reconstruction', 
                                                   'time-parallel')] = rep_time
  
          # injuries
          if DVs['injuries']:
              if 'CM' in self._COL.columns:
                  SUMMARY.loc[colID, ('collapses', 'mode')] = self._COL.loc[:, 'CM']
              
                  SUMMARY.loc[colID, ('injuries', 'casualties')] = \
                      self._COL.loc[:, 'INJ-0']
                  SUMMARY.loc[colID, ('injuries', 'fatalities')] = \
                      self._COL.loc[:, 'INJ-1']
      
              SUMMARY.loc[ncID, ('injuries', 'casualties')] = \
                  self._DV_dict['injuries'][0].sum(axis=1)
              SUMMARY.loc[ncID, ('injuries', 'fatalities')] = \
                  self._DV_dict['injuries'][1].sum(axis=1)
          
          self._SUMMARY = SUMMARY.dropna(axis=1,how='all')
  
      def write_outputs(self):
          """
          
          Returns
          -------
  
          """
          super(FEMA_P58_Assessment, self).write_outputs()
  
      def _create_correlation_matrix(self, rho_target, c_target=-1, 
                                     include_CSG=False, 
                                     include_DSG=False, include_DS=False):
          """
          
          Parameters
          ----------
          rho_target
          c_target
          include_CSG
          include_DSG
          include_DS
  
          Returns
          -------
  
          """
  
          # set the correlation structure
          rho_FG, rho_PG, rho_LOC, rho_DIR, rho_CSG, rho_DS = np.zeros(6)
  
          if rho_target in ['FG', 'PG', 'DIR', 'LOC', 'CSG', 'ATC', 'DS']:
              rho_DS = 1.0
          if rho_target in ['FG', 'PG', 'DIR', 'LOC', 'CSG']:
              rho_CSG = 1.0
          if rho_target in ['FG', 'PG', 'DIR']:
              rho_DIR = 1.0
          if rho_target in ['FG', 'PG', 'LOC']:
              rho_LOC = 1.0
          if rho_target in ['FG', 'PG']:
              rho_PG = 1.0
          if rho_target == 'FG':
              rho_FG = 1.0
  
          L_D_list = []
          dims = []
          DS_list = []
          ATC_rho = []
          s_fg_keys = sorted(self._FG_in.keys())
          for c_id, c_name in enumerate(s_fg_keys):
              comp = self._FG_in[c_name]
              
              if ((c_target == -1) or (c_id == c_target)):
                  c_L_D_list = []
                  c_DS_list = []
                  ATC_rho.append(comp['correlation'])
  
                  if include_DSG:
                      DS_count = 0
                      s_dsg_keys = sorted(comp['DSG_set'].keys())
                      for dsg_i in s_dsg_keys:
                          DSG = comp['DSG_set'][dsg_i]
                          if include_DS:
                              DS_count += len(DSG['DS_set'])
                          else:
                              DS_count += 1
                  else:
                      DS_count = 1
  
                  for loc in comp['locations']:
                      if include_CSG:
                          u_dirs = comp['directions']
                      else:
                          u_dirs = np.unique(comp['directions'])
                      c_L_D_list.append([])
                      for dir_ in u_dirs:
                          c_DS_list.append(DS_count)
                          for ds_i in range(DS_count):
                              c_L_D_list[-1].append(dir_)
  
                  c_dims = sum([len(loc) for loc in c_L_D_list])
                  dims.append(c_dims)
                  L_D_list.append(c_L_D_list)
                  DS_list.append(c_DS_list)
  
          rho = np.ones((sum(dims), sum(dims))) * rho_FG
  
          f_pos_id = 0
          for c_id, (c_L_D_list, c_dims, c_DS_list) in enumerate(
              zip(L_D_list, dims, DS_list)):
              c_rho = np.ones((c_dims, c_dims)) * rho_PG
  
              # dependencies btw directions
              if rho_DIR != 0:
                  c_pos_id = 0
                  for loc_D_list in c_L_D_list:
                      l_dim = len(loc_D_list)
                      c_rho[c_pos_id:c_pos_id + l_dim,
                      c_pos_id:c_pos_id + l_dim] = rho_DIR
                      c_pos_id = c_pos_id + l_dim
  
              # dependencies btw locations
              if rho_LOC != 0:
                  flat_dirs = []
                  [[flat_dirs.append(dir_i) for dir_i in dirs] for dirs in
                   c_L_D_list]
                  flat_dirs = np.array(flat_dirs)
                  for u_dir in np.unique(flat_dirs):
                      dir_ids = np.where(flat_dirs == u_dir)[0]
                      for i in dir_ids:
                          for j in dir_ids:
                              c_rho[i, j] = rho_LOC
  
              if ((rho_CSG != 0) or (rho_target == 'ATC')):
                  c_pos_id = 0
                  if rho_target == 'ATC':
                      rho_to_use = float(ATC_rho[c_id])
                  else:
                      rho_to_use = rho_CSG
                  for loc_D_list in c_L_D_list:
                      flat_dirs = np.array(loc_D_list)
                      for u_dir in np.unique(flat_dirs):
                          dir_ids = np.where(flat_dirs == u_dir)[0]
                          for i in dir_ids:
                              for j in dir_ids:
                                  c_rho[c_pos_id + i, c_pos_id + j] = rho_to_use
                      c_pos_id = c_pos_id + len(loc_D_list)
  
              if rho_DS != 0:
                  c_pos_id = 0
                  for l_dim in c_DS_list:
                      c_rho[c_pos_id:c_pos_id + l_dim,
                            c_pos_id:c_pos_id + l_dim] = rho_DS
                      c_pos_id = c_pos_id + l_dim
  
              rho[f_pos_id:f_pos_id + c_dims,
                  f_pos_id:f_pos_id + c_dims] = c_rho
              f_pos_id = f_pos_id + c_dims
  
          np.fill_diagonal(rho, 1.0)
  
          return rho
  
      def _create_RV_quantities(self, rho_qnt):
          """
          
          Parameters
          ----------
          rho_qnt
  
          Returns
          -------
  
          """
  
          q_theta, q_sig, q_tag, q_dist = [np.array([]) for i in range(4)]
  
          # collect the parameters for each quantity dimension
          s_fg_keys = sorted(self._FG_in.keys())
          for c_id in s_fg_keys:
              comp = self._FG_in[c_id]
              u_dirs = np.unique(comp['directions'])
  
              dir_weights = comp['dir_weights']
              theta_list = []
              [[theta_list.append(qnt * dw) for dw in dir_weights]
               for qnt in comp['quantities']]
              q_theta = np.append(q_theta, theta_list)
  
              if comp['distribution_kind'] == 'normal':
                  q_sig = np.append(q_sig, (
                      comp['cov'] * np.asarray(theta_list)).tolist())
              else:
                  q_sig = np.append(q_sig, (
                      np.ones(len(theta_list)) * comp['cov']).tolist())
  
              q_tag = np.append(q_tag,
                                [[c_id + '-QNT-' + str(s_i) + '-' + str(d_i)
                                  for d_i in u_dirs]
                                 for s_i in comp['locations']])
              q_dist = np.append(q_dist,
                                 [[comp['distribution_kind'] for d_i in u_dirs]
                                  for s_i in comp['locations']])
  
          dims = len(q_theta)
          rho = self._create_correlation_matrix(rho_qnt)
          q_COV = np.outer(q_sig, q_sig) * rho
  
          # add lower limits to ensure only positive quantities
          # zero is probably too low, and it might make sense to introduce upper 
          # limits as well
          tr_lower = [0. for d in range(dims)]
          tr_upper = [None for d in range(dims)]
          # to avoid truncations affecting other dimensions when rho_QNT is large, 
          # assign a post-truncation correlation structure
          corr_ref = 'post'
  
          # create a random variable for component quantities in performance groups
          if q_tag.size > 0:
              quantity_RV = RandomVariable(ID=100,
                                           dimension_tags=q_tag,
                                           distribution_kind=q_dist,
                                           theta=q_theta,
                                           COV=q_COV,
                                           truncation_limits=[tr_lower, tr_upper],
                                           corr_ref=corr_ref)
          else:
              quantity_RV = None
  
          return quantity_RV
  
      def _create_RV_fragilities(self, c_id, comp, rho_fr):
          """
          
          Parameters
          ----------
          c_id
          comp
          rho_fr
  
          Returns
          -------
  
          """
  
          # prepare the basic multivariate distribution data for one component subgroup considering all damage states
          d_theta, d_sig, d_tag = [np.array([]) for i in range(3)]
  
          s_dsg_keys = sorted(comp['DSG_set'].keys())
          for d_id in s_dsg_keys:
              DSG = comp['DSG_set'][d_id]
              d_theta = np.append(d_theta, DSG['theta'])
              d_sig = np.append(d_sig, DSG['sig'])
              d_tag = np.append(d_tag, comp['ID'] + '-' + str(d_id))
          dims = len(d_theta)
  
          # get the total number of random variables for this fragility group
          rv_count = len(comp['locations']) * len(comp['directions']) * dims
  
          # create the (empty) input arrays for the RV
          c_theta = np.zeros(rv_count)
          c_tag = np.empty(rv_count, dtype=object)
          c_sig = np.zeros(rv_count)
  
          pos_id = 0
          for l_id in comp['locations']:
              # for each location-direction pair)
              for d_id, __ in enumerate(comp['directions']):
                  # for each component-subgroup
                  c_theta[pos_id:pos_id + dims] = d_theta
                  c_sig[pos_id:pos_id + dims] = d_sig
                  c_tag[pos_id:pos_id + dims] = [
                      t + '-LOC-{}-CSG-{}'.format(l_id, d_id) for t in d_tag]
                  pos_id += dims
  
          # create the covariance matrix
          c_rho = self._create_correlation_matrix(rho_fr, c_target=c_id, 
                                                  include_DSG=True,
                                                  include_CSG=True)
          c_COV = np.outer(c_sig, c_sig) * c_rho
  
          if c_tag.size > 0:
              fragility_RV = RandomVariable(ID=300 + c_id,
                                            dimension_tags=c_tag,
                                            distribution_kind='lognormal',
                                            theta=c_theta,
                                            COV=c_COV)
          else:
              fragility_RV = None
  
          return fragility_RV
  
      def _create_RV_red_tags(self, rho_target):
  
          f_theta, f_sig, f_tag = [np.array([]) for i in range(3)]
          s_fg_keys = sorted(self._FG_in.keys())
          for c_id, c_name in enumerate(s_fg_keys):
              comp = self._FG_in[c_name]
  
              d_theta, d_sig, d_tag = [np.array([]) for i in range(3)]
  
              s_dsg_keys = sorted(comp['DSG_set'].keys())
              for dsg_i in s_dsg_keys:
                  DSG = comp['DSG_set'][dsg_i]
                  s_ds_keys = sorted(DSG['DS_set'].keys())
                  for ds_i in s_ds_keys:
                      DS = DSG['DS_set'][ds_i]
                      theta = DS['red_tag']['theta']
                      d_theta = np.append(d_theta, theta)
                      d_sig = np.append(d_sig, DS['red_tag']['cov'])
                      d_tag = np.append(d_tag,
                                        comp['ID'] + '-' + str(dsg_i) + '-' + str(
                                            ds_i))
  
              for loc in comp['locations']:
                  for dir_ in np.unique(comp['directions']):
                      f_theta = np.append(f_theta, d_theta)
                      f_sig = np.append(f_sig, d_sig)
                      f_tag = np.append(f_tag,
                                        [t + '-LOC-{}-DIR-{}'.format(loc, dir_)
                                         for t in d_tag])
  
          rho = self._create_correlation_matrix(rho_target, c_target=-1, 
                                                include_DSG=True,
                                                include_DS=True)
  
          # remove the unnecessary fields
          to_remove = np.where(f_theta == 0)[0]
          rho = np.delete(rho, to_remove, axis=0)
          rho = np.delete(rho, to_remove, axis=1)
  
          f_theta, f_sig, f_tag = [np.delete(f_vals, to_remove)
                                   for f_vals in [f_theta, f_sig, f_tag]]
  
          f_COV = np.outer(f_sig, f_sig) * rho
  
          tr_upper = 1. + (1. - f_theta) / f_theta
  
          if f_tag.size > 0:
              red_tag_RV = RandomVariable(ID=400,
                                          dimension_tags=f_tag,
                                          distribution_kind='normal',
                                          theta=np.ones(len(f_theta)),
                                          COV=f_COV,
                                          corr_ref='post',
                                          truncation_limits=[np.zeros(len(f_theta)),
                                                             tr_upper])
          else:
              red_tag_RV = None
  
          return red_tag_RV
  
      def _create_RV_repairs(self, rho_cost, rho_time, rho_cNt):
  
          # prepare the cost and time parts of the data separately
          ct_sig, ct_tag, ct_dkind = [np.array([]) for i in range(3)]
          for rho_target, name in zip([rho_cost, rho_time], ['cost', 'time']):
              
              f_sig, f_tag, f_dkind = [np.array([]) for i in range(3)]
              
              s_fg_keys = sorted(self._FG_in.keys())
              for c_id, c_name in enumerate(s_fg_keys):
                  comp = self._FG_in[c_name]
  
                  d_sig, d_tag, d_dkind = [np.array([]) for i in range(3)]
  
                  s_dsg_keys = sorted(comp['DSG_set'].keys())
                  for dsg_i in s_dsg_keys:
                      DSG = comp['DSG_set'][dsg_i]
                      s_ds_keys = sorted(DSG['DS_set'].keys())
                      for ds_i in s_ds_keys:
                          DS = DSG['DS_set'][ds_i]
                          d_sig = np.append(d_sig,
                                            DS['repair_{}'.format(name)]['cov'])
                          d_dkind = np.append(d_dkind,
                                              DS['repair_{}'.format(name)][
                                                  'distribution_kind'])
                          d_tag = np.append(d_tag,
                                            comp['ID'] + '-' + str(
                                                dsg_i) + '-' + str(
                                                ds_i) + '-{}'.format(name))
  
                  for loc in comp['locations']:
                      for dir_ in np.unique(comp['directions']):
                          f_sig = np.append(f_sig, d_sig)
                          f_dkind = np.append(f_dkind, d_dkind)
                          f_tag = np.append(f_tag,
                                            [t + '-LOC-{}-DIR-{}'.format(loc,
                                                                         dir_)
                                             for t in d_tag])
              ct_sig = np.append(ct_sig, f_sig)
              ct_tag = np.append(ct_tag, f_tag)
              ct_dkind = np.append(ct_dkind, f_dkind)
  
          rho_c = self._create_correlation_matrix(rho_cost, c_target=-1,
                                            include_DSG=True,
                                            include_DS=True)
          rho_t = self._create_correlation_matrix(rho_time, c_target=-1,
                                            include_DSG=True,
                                            include_DS=True)
  
          dims = len(ct_tag)
          ct_rho = np.zeros((dims, dims))
  
          dims = dims // 2
          if rho_cNt == False:
              ct_rho[:dims, :dims] = rho_c
              ct_rho[dims:, dims:] = rho_t
          else:    
              # if correlation between cost and time is considered, then the 
              # envelope of rho_c and rho_t shall be applied
              rho_ct = np.maximum(rho_c, rho_t)
              ct_rho[:dims, :dims] = rho_ct
              ct_rho[dims:, dims:] = rho_ct
              
              # apply the same blocks to the off-diagonal positions
              ct_rho[:dims, dims:] = rho_ct
              ct_rho[dims:, :dims] = rho_ct
  
          ct_COV = np.outer(ct_sig, ct_sig) * ct_rho
  
          if ct_tag.size > 0:
              repair_RV = RandomVariable(ID=401,
                                         dimension_tags=ct_tag,
                                         distribution_kind=ct_dkind,
                                         theta=np.ones(len(ct_sig)),
                                         COV=ct_COV,
                                         corr_ref='post',
                                         truncation_limits=[np.zeros(len(ct_sig)),
                                                            None])
          else:
              repair_RV = None
  
          return repair_RV
  
      def _create_RV_injuries(self, rho_target, rho_lvls):
  
          inj_lvls = self._inj_lvls
  
          # prepare the cost and time parts of the data separately
          full_theta, full_sig, full_tag = [np.array([]) for i in range(3)]
          for i_lvl in range(inj_lvls):
  
              f_theta, f_sig, f_tag = [np.array([]) for i in range(3)]
              s_fg_keys = sorted(self._FG_in.keys())
              for c_id, c_name in enumerate(s_fg_keys):
                  comp = self._FG_in[c_name]
  
                  d_theta, d_sig, d_tag = [np.array([]) for i in range(3)]
  
                  s_dsg_keys = sorted(comp['DSG_set'].keys())
                  for dsg_i in s_dsg_keys:
                      DSG = comp['DSG_set'][dsg_i]
                      s_ds_keys = sorted(DSG['DS_set'].keys())
                      for ds_i in s_ds_keys:
                          DS = DSG['DS_set'][ds_i]
                          d_theta = np.append(d_theta, DS['injuries']['theta'][i_lvl])
                          d_sig = np.append(d_sig, DS['injuries']['cov'][i_lvl])
                          d_tag = np.append(d_tag,
                                            comp['ID'] + '-' + str(
                                                dsg_i) + '-' + str(
                                                ds_i) + '-{}'.format(i_lvl))
  
                  for loc in comp['locations']:
                      for dir_ in np.unique(comp['directions']):
                          f_theta = np.append(f_theta, d_theta)
                          f_sig = np.append(f_sig, d_sig)
                          f_tag = np.append(f_tag,
                                            [t + '-LOC-{}-DIR-{}'.format(loc,
                                                                         dir_)
                                             for t in d_tag])
  
              full_theta = np.append(full_theta, f_theta)
              full_sig = np.append(full_sig, f_sig)
              full_tag = np.append(full_tag, f_tag)
  
          dims = len(full_tag)
          full_rho = np.zeros((dims, dims))
          dims = dims // inj_lvls
  
          # if correlation between different levels of severity is considered, add that to the matrix
          if rho_lvls:
              rho_i = self._create_correlation_matrix(rho_target, c_target=-1,
                                                include_DSG=True,
                                                include_DS=True)            
              for i in range(inj_lvls):
                  for j in range(inj_lvls):
                      full_rho[i * dims:(i + 1) * dims,
                      j * dims:(j + 1) * dims] = rho_i
  
          # and now add the values around the main diagonal
          for i in range(inj_lvls):
              rho_i = self._create_correlation_matrix(rho_target, c_target=-1,
                                                include_DSG=True,
                                                include_DS=True)
              full_rho[i * dims:(i + 1) * dims, i * dims:(i + 1) * dims] = rho_i
  
              # finally, remove the unnecessary lines
          to_remove = np.where(full_theta == 0)[0]
          full_rho = np.delete(full_rho, to_remove, axis=0)
          full_rho = np.delete(full_rho, to_remove, axis=1)
  
          full_theta, full_sig, full_tag = [np.delete(f_vals, to_remove)
                                            for f_vals in
                                            [full_theta, full_sig, full_tag]]
  
          full_COV = np.outer(full_sig, full_sig) * full_rho
  
          tr_upper = 1. + (1. - full_theta) / full_theta
  
          if full_tag.size > 0:
              injury_RV = RandomVariable(ID=402,
                                         dimension_tags=full_tag,
                                         distribution_kind='normal',
                                         theta=np.ones(len(full_sig)),
                                         COV=full_COV,
                                         corr_ref='post',
                                         truncation_limits=[np.zeros(len(full_sig)),
                                                            tr_upper])
          else:
              injury_RV = None
  
          return injury_RV
  
      def _create_RV_demands(self, beta_added):
  
          # unlike other random variables, the demand RV is based on raw data
          # first, collect the raw values from the EDP dict
          demand_data = []
          d_tags = []
          detection_limits = []
          s_edp_keys = sorted(self._EDP_in.keys())
          for d_id in s_edp_keys:
              d_list = self._EDP_in[d_id]
              for i in range(len(d_list)):
                  demand_data.append(d_list[i]['raw_data'])
                  d_tags.append(str(d_id) + 
                                '-LOC-' + str(d_list[i]['location']) + 
                                '-DIR-' + str(d_list[i]['direction']))
                  det_lim = self._AIM_in['general']['detection_limits'][d_id]
                  if det_lim is None:
                      det_lim = np.inf
                  detection_limits.append([0., det_lim])
  
          detection_limits = np.transpose(np.asarray(detection_limits))
          demand_data = np.transpose(np.asarray(demand_data))
          
          # if more than one sample is provided
          if demand_data.shape[0] > 1:
              # get the number of censored samples
              EDP_filter = np.all([np.all(demand_data > detection_limits[0], axis=1),
                                   np.all(demand_data < detection_limits[1], axis=1)],
                                  axis=0)
              censored_count = len(EDP_filter) - sum(EDP_filter)
              demand_data = demand_data[EDP_filter]
              demand_data = np.transpose(demand_data)
      
              # create the random variable
              demand_RV = RandomVariable(ID=200, dimension_tags=d_tags,
                                         raw_data=demand_data,
                                         detection_limits=detection_limits,
                                         censored_count=censored_count
                                         )
      
              # fit a multivariate lognormal distribution to the censored raw data
              demand_RV.fit_distribution('lognormal')
          else:
              # Since we only have one data point, the best we can do is assume
              # it is the median of the multivariate distribution. The dispersion
              # is assumed to be negligible.
              dim = len(demand_data[0])
              if dim > 1:
                  sig = np.ones(dim)*1e-4
                  rho = np.zeros((dim,dim))
                  np.fill_diagonal(rho, 1.0)
                  COV = np.outer(sig,sig) * rho
              else:
                  COV = np.asarray((1e-4)**2.)
                  
              demand_RV = RandomVariable(ID=200, dimension_tags=d_tags,
                                         distribution_kind='lognormal',
                                         theta=demand_data[0],
                                         COV=COV)
  
          # if we want to add other sources of uncertainty, we will need to 
          # redefine the random variable
          if beta_added > 0.:
              # get the covariance matrix with added uncertainty
              COV_orig = demand_RV.COV
              if COV_orig.shape is not ():
                  sig_orig = np.sqrt(np.diagonal(COV_orig))
                  rho_orig = COV_orig / np.outer(sig_orig, sig_orig)
                  sig_mod = np.sqrt(sig_orig ** 2. + beta_added ** 2.)
                  COV_mod = np.outer(sig_mod, sig_mod) * rho_orig
              else:
                  COV_mod = np.sqrt(COV_orig**2. + beta_added**2.)
  
              # redefine the random variable
              demand_RV = RandomVariable(ID=200,
                                         dimension_tags=demand_RV.dimension_tags,
                                         distribution_kind='lognormal',
                                         theta=demand_RV.theta,
                                         COV=COV_mod)
  
          return demand_RV
  
      def _create_fragility_groups(self):
  
          RVd = self._RV_dict
          DVs = self._AIM_in['decision_variables']
  
          # create a list for the fragility groups
          FG_dict = dict()
  
          s_fg_keys = sorted(self._FG_in.keys())
          for c_id in s_fg_keys:
              comp = self._FG_in[c_id]
  
              FG_ID = len(FG_dict.keys())
  
              # create a list for the performance groups
              performance_groups = []
  
              # one group for each of the stories prescribed by the user
              PG_locations = comp['locations']
              PG_directions = np.unique(comp['directions'])
              for loc in PG_locations:
                  for dir_ in PG_directions:
                      PG_ID = 1000 + dir_*100 + FG_ID * 10 + loc
      
                      # get the quantity
                      QNT = RandomVariableSubset(
                          RVd['QNT'],
                          tags=[c_id + '-QNT-' + str(loc) + '-' + str(dir_), ])
      
                      # create the damage objects
                      # consequences are calculated on a performance group level
      
                      # create a list for the damage state groups and their tags
                      DSG_list = []
                      d_tags = []
                      s_dsg_keys = sorted(comp['DSG_set'].keys())
                      for dsg_i, DSG_ID in enumerate(s_dsg_keys):
                          DSG = comp['DSG_set'][DSG_ID]
                          d_tags.append(c_id + '-' + DSG_ID)
      
                          # create a list for the damage states
                          DS_set = []
                          
                          s_ds_keys = sorted(DSG['DS_set'].keys())
                          for ds_i, DS_ID in enumerate(s_ds_keys):
                              DS = DSG['DS_set'][DS_ID]
      
                              # create the consequence functions
                              if DVs['rec_cost']:
                                  data = DS['repair_cost']
                                  f_median = prep_bounded_linear_median_DV(
                                      **{k: data.get(k, None) for k in
                                         ('median_max', 'median_min',
                                          'quantity_lower', 'quantity_upper')})
                                  cf_tag = c_id + '-' + DSG_ID + '-' + DS_ID + \
                                           '-cost' + \
                                           '-LOC-{}-DIR-{}'.format(loc, dir_)
                                  CF_RV = RandomVariableSubset(RVd['DV_REP'],
                                                               tags=cf_tag)
                                  CF_cost = ConsequenceFunction(DV_median=f_median,
                                                                DV_distribution=CF_RV)
                              else:
                                  CF_cost = None
      
                              if DVs['rec_time']:
                                  data = DS['repair_time']
                                  f_median = prep_bounded_linear_median_DV(
                                      **{k: data.get(k, None) for k in
                                         ('median_max', 'median_min', 'quantity_lower',
                                          'quantity_upper')})
                                  cf_tag = c_id + '-' + DSG_ID + '-' + DS_ID + \
                                           '-time' + \
                                           '-LOC-{}-DIR-{}'.format(loc, dir_)
                                  CF_RV = RandomVariableSubset(RVd['DV_REP'],
                                                               tags=cf_tag)
                                  CF_time = ConsequenceFunction(DV_median=f_median,
                                                                DV_distribution=CF_RV)
                              else:
                                  CF_time = None
      
                              if DVs['red_tag']:
                                  data = DS['red_tag']
                                  if data['theta'] > 0:
                                      f_median = prep_constant_median_DV(data['theta'])
                                      cf_tag = c_id + '-' + DSG_ID + '-' + DS_ID + \
                                               '-LOC-{}-DIR-{}'.format(loc, dir_)
                                      CF_RV = RandomVariableSubset(RVd['DV_RED'],
                                                                   tags=cf_tag)
                                      CF_red_tag = ConsequenceFunction(DV_median=f_median,
                                                                       DV_distribution=CF_RV)
                                  else:
                                      CF_red_tag = None
                              else:
                                  CF_red_tag = None
      
                              if DVs['injuries']:
                                  data = DS['injuries']
                                  CF_inj_set = []
                                  for inj_i, theta in enumerate(data['theta']):
                                      if theta > 0.:
                                          f_median = prep_constant_median_DV(theta)
                                          cf_tag = c_id + '-' + DSG_ID + '-' + DS_ID + \
                                                   '-{}-LOC-{}-DIR-{}'.format(inj_i, loc, dir_)
                                          CF_RV = RandomVariableSubset(RVd['DV_INJ'],
                                                                       tags=cf_tag)
                                          CF_inj_set.append(ConsequenceFunction(
                                              DV_median=f_median,
                                              DV_distribution=CF_RV))
                                      else:
                                          CF_inj_set.append(None)
                              else:
                                  CF_inj_set = [None,]
      
                              # add the DS to the list
                              DS_set.append(DamageState(ID=ds_i + 1,
                                                        description=DS['description'],
                                                        weight=DS['weight'],
                                                        affected_area=DS[
                                                            'affected_area'],
                                                        repair_cost_CF=CF_cost,
                                                        reconstruction_time_CF=CF_time,
                                                        red_tag_CF=CF_red_tag,
                                                        injuries_CF_set=CF_inj_set
                                                        ))
      
                          # add the DSG to the list
                          DSG_list.append(DamageStateGroup(ID=dsg_i + 1,
                                                           DS_set=DS_set,
                                                           DS_set_kind=DSG[
                                                               'DS_set_kind'],
                                                           description=DSG[
                                                               'description']
                                                           ))
      
                      # create the fragility functions
                      FF_set = []
                      CSG_this = np.where(comp['directions']==dir_)[0]
                      PG_weights = np.asarray(comp['csg_weights'])[CSG_this]
                      # normalize the weights
                      PG_weights /= sum(PG_weights)
                      for csg_id in CSG_this:
                          # assign the appropriate random variable to the fragility 
                          # function
                          ff_tags = [t + '-LOC-{}-CSG-{}'.format(loc, csg_id)
                                     for t in d_tags]
                          EDP_limit = RandomVariableSubset(RVd['FR-' + c_id],
                                                           tags=ff_tags)
                          FF_set.append(FragilityFunction(EDP_limit))
      
                      # create the performance group
                      PG = PerformanceGroup(ID=PG_ID,
                                            location=loc,
                                            quantity=QNT,
                                            fragility_functions=FF_set,
                                            DSG_set=DSG_list,
                                            csg_weights=PG_weights,
                                            direction=dir_
                                            )
                      performance_groups.append(PG)
  
              # create the fragility group
              FG = FragilityGroup(ID=FG_ID,
                                  kind=comp['kind'],
                                  demand_type=comp['demand_type'],
                                  performance_groups=performance_groups,
                                  directional=comp['directional'],
                                  correlation=comp['correlation'],
                                  demand_location_offset=comp['offset'],
                                  incomplete=comp['incomplete'],
                                  name=str(FG_ID) + ' - ' + comp['ID'],
                                  description=comp['description']
                                  )
  
              FG_dict.update({comp['ID']:FG})
  
          return FG_dict
  
      def _sample_event_time(self):
          
          sample_count = self._AIM_in['general']['realizations']
          
          # month - uniform distribution over [0,11]
          month = np.random.randint(0, 12, size=sample_count)
  
          # weekday - binomial with p=5/7
          weekday = np.random.binomial(1, 5. / 7., size=sample_count)
  
          # hour - uniform distribution over [0,23]
          hour = np.random.randint(0, 24, size=sample_count)
  
          data = pd.DataFrame(data={'month'   : month,
                                    'weekday?': weekday,
                                    'hour'    : hour},
                              dtype=int)
  
          return data
  
      def _get_population(self):
          """
          Use the population characteristics to generate random population samples.
  
          Returns
          -------
          
          """
          POPin = self._POP_in
          TIME = self._TIME
          
          POP = pd.DataFrame(
              np.ones((len(TIME.index), len(POPin['peak']))) * POPin['peak'],
              columns=['LOC' + str(loc + 1)
                       for loc in range(len(POPin['peak']))])
  
          weekdays = TIME[TIME['weekday?'] == 1].index
          weekends = TIME[~TIME.index.isin(weekdays)].index
  
          for col in POP.columns.values:
              POP.loc[weekdays, col] = (
                  POP.loc[weekdays, col] *
                  np.array(POPin['weekday']['daily'])[
                      TIME.loc[weekdays, 'hour'].values.astype(int)] *
                  np.array(POPin['weekday']['monthly'])[
                      TIME.loc[weekdays, 'month'].values.astype(int)])
  
              POP.loc[weekends, col] = (
                  POP.loc[weekends, col] *
                  np.array(POPin['weekend']['daily'])[
                      TIME.loc[weekends, 'hour'].values.astype(int)] *
                  np.array(POPin['weekend']['monthly'])[
                      TIME.loc[weekends, 'month'].values.astype(int)])
  
          return POP
  
      def _calc_collapses(self):
  
          # filter the collapsed cases based on the demand samples
          collapsed_IDs = np.array([])
          s_edp_keys = sorted(self._EDP_dict.keys())
          for demand_ID in s_edp_keys:
              demand = self._EDP_dict[demand_ID]
              coll_df = pd.DataFrame()
              kind = demand_ID[:3]
              collapse_limit = self._AIM_in['general']['collapse_limits'][kind]
              if collapse_limit is not None:
                  EDP_samples = demand.samples
                  coll_df = EDP_samples[EDP_samples > collapse_limit]
              collapsed_IDs = np.concatenate(
                  (collapsed_IDs, coll_df.index.values))
  
          # get a list of IDs of the collapsed cases
          collapsed_IDs = np.unique(collapsed_IDs).astype(int)
  
          COL = pd.DataFrame(np.zeros(self._AIM_in['general']['realizations']), 
                             columns=['COL', ])
          COL.loc[collapsed_IDs, 'COL'] = 1
  
          return COL, collapsed_IDs
  
      def _calc_damage(self):
          
          ncID = self._ID_dict['non-collapse']
          NC_samples = len(ncID)
          DMG = pd.DataFrame()
  
          s_fg_keys = sorted(self._FG_dict.keys())
          for fg_id in s_fg_keys:
              FG = self._FG_dict[fg_id]
  
              PG_set = FG._performance_groups
  
              DS_list = []
              for DSG in PG_set[0]._DSG_set:
                  for DS in DSG._DS_set:
                      DS_list.append(str(DSG._ID) + '-' + str(DS._ID))
              d_count = len(DS_list)
  
              MI = pd.MultiIndex.from_product([[FG._ID, ],
                                               [pg._ID for pg in PG_set],
                                               DS_list],
                                              names=['FG', 'PG', 'DS'])
  
              FG_damages = pd.DataFrame(np.zeros((NC_samples, len(MI))),
                                        columns=MI,
                                        index=ncID)
  
              for pg_i, PG in enumerate(PG_set):
  
                  PG_ID = PG._ID
                  PG_qnt = PG._quantity.samples.loc[ncID]
  
                  # get the corresponding demands
                  demand_ID = (FG._demand_type + 
                               '-LOC-' + str(PG._location) +
                               '-DIR-' + str(PG._direction+1))
                  if demand_ID in self._EDP_dict.keys():
                      EDP_samples = self._EDP_dict[demand_ID].samples.loc[ncID]
                  else:
                      # If the required demand is not available, then we are most
                      # likely analyzing a 3D structure using results from a 2D 
                      # simulation. The best thing we can do in that particular
                      # case is to use the EDP from the 0 direction for all other
                      # directions.
                      demand_ID = (FG._demand_type + 
                                   '-LOC-' + str(PG._location) + '-DIR-1')
                      EDP_samples = self._EDP_dict[demand_ID].samples.loc[ncID]
                      
                  csg_w_list = PG._csg_weights
  
                  for csg_i, csg_w in enumerate(csg_w_list): 
                      DSG_df = PG._FF_set[csg_i].DSG_given_EDP(EDP_samples)
                      
                      for DSG in PG._DSG_set:
                          in_this_DSG = DSG_df[DSG_df.values == DSG._ID].index
                          if DSG._DS_set_kind == 'single':
                              DS = DSG._DS_set[0]
                              DS_tag = str(DSG._ID) + '-' + str(DS._ID)
                              FG_damages.loc[in_this_DSG,
                                             (FG._ID, PG_ID, DS_tag)] += csg_w
                          elif DSG._DS_set_kind == 'mutually exclusive':
                              DS_weights = [DS._weight for DS in DSG._DS_set]
                              DS_RV = RandomVariable(
                                  ID=-1, dimension_tags=['me_DS', ],
                                  distribution_kind='multinomial',
                                  p_set=DS_weights)
                              DS_df = DS_RV.sample_distribution(
                                  len(in_this_DSG)) + 1
                              for DS in DSG._DS_set:
                                  DS_tag = str(DSG._ID) + '-' + str(DS._ID)
                                  in_this_DS = DS_df[DS_df.values == DS._ID].index
                                  FG_damages.loc[in_this_DSG[in_this_DS],
                                                 (FG._ID, PG_ID, DS_tag)] += csg_w
                          else:
                              # TODO: simultaneous
                              print(DSG._DS_set_kind)
  
                  FG_damages.iloc[:,pg_i * d_count:(pg_i + 1) * d_count] = \
                      FG_damages.mul(PG_qnt.iloc[:, 0], axis=0)
  
              DMG = pd.concat((DMG, FG_damages), axis=1)
  
          DMG.index = ncID
          
          # sort the columns to enable index slicing later
          DMG = DMG.sort_index(axis=1, ascending=True)
  
          return DMG
  
      def _calc_red_tag(self):
          idx = pd.IndexSlice
  
          ncID = self._ID_dict['non-collapse']
          NC_samples = len(ncID)
          DV_RED = pd.DataFrame()
  
          s_fg_keys = sorted(self._FG_dict.keys())
          for fg_id in s_fg_keys:
              FG = self._FG_dict[fg_id]
  
              PG_set = FG._performance_groups
  
              DS_list = self._DMG.loc[:, idx[FG._ID, PG_set[0]._ID, :]].columns
              DS_list = DS_list.levels[2][DS_list.labels[2]].values
  
              MI = pd.MultiIndex.from_product([[FG._ID, ],
                                               [pg._ID for pg in PG_set],
                                               DS_list],
                                              names=['FG', 'PG', 'DS'])
  
              FG_RED = pd.DataFrame(np.zeros((NC_samples, len(MI))),
                                    columns=MI,
                                    index=ncID)
  
              for pg_i, PG in enumerate(PG_set):
  
                  PG_ID = PG._ID
                  PG_qnt = PG._quantity.samples.loc[ncID]
  
                  PG_DMG = self._DMG.loc[:, idx[FG._ID, PG_ID, :]].div(
                      PG_qnt.iloc[:, 0],
                      axis=0)
  
                  for d_i, d_tag in enumerate(DS_list):
                      dsg_i = int(d_tag[0]) - 1
                      ds_i = int(d_tag[-1]) - 1
  
                      DS = PG._DSG_set[dsg_i]._DS_set[ds_i]
  
                      if DS._red_tag_CF is not None:
                          RED_samples = DS.red_tag_dmg_limit(
                              sample_size=NC_samples)
                          RED_samples.index = ncID
  
                          is_red = PG_DMG.loc[:, (FG._ID, PG_ID, d_tag)].sub(
                              RED_samples, axis=0)
                          
                          FG_RED.loc[:, (FG._ID, PG_ID, d_tag)] = (
                              is_red > 0.).astype(int)
                      else:
                          FG_RED.drop(labels=[(FG._ID, PG_ID, d_tag), ], axis=1,
                                      inplace=True)
  
              if FG_RED.size > 0:
                  DV_RED = pd.concat((DV_RED, FG_RED), axis=1)
  
          # sort the columns to enable index slicing later
          DV_RED = DV_RED.sort_index(axis=1, ascending=True)
  
          return DV_RED
  
      def _calc_irrepairable(self):
  
          ncID = self._ID_dict['non-collapse']
          NC_samples = len(ncID)
  
          # determine which realizations lead to irrepairable damage
          # get the max residual drifts
          RED_max = None
          PID_max = None
          s_edp_keys = sorted(self._EDP_dict.keys())
          for demand_ID in s_edp_keys:
              demand = self._EDP_dict[demand_ID]
              kind = demand_ID[:3]
              if kind == 'RED':
                  r_max = demand.samples.loc[ncID].values
                  if RED_max is None:
                      RED_max = r_max
                  else:
                      RED_max = np.max((RED_max, r_max), axis=0)
              elif kind == 'PID':
                  d_max = demand.samples.loc[ncID].values
                  if PID_max is None:
                      PID_max = d_max
                  else:
                      PID_max = np.max((PID_max, d_max), axis=0)
  
          if (RED_max is None) and (PID_max is not None):
              # we need to estimate residual drifts based on peak drifts
              RED_max = np.zeros(NC_samples)
  
              # based on Appendix C in FEMA P-58
              delta_y = self._AIM_in['general']['yield_drift']
              small = PID_max < delta_y
              medium = PID_max < 4 * delta_y
              large = PID_max >= 4 * delta_y
  
              RED_max[large] = PID_max[large] - 3 * delta_y
              RED_max[medium] = 0.3 * (PID_max[medium] - delta_y)
              RED_max[small] = 0.
          else:
              # If no drift data is available, then we cannot provide an estimate
              # of irrepairability. We assume that all non-collapse realizations 
              # are repairable in this case.
              return np.array([])
              
          # get the probabilities of irrepairability
          irrep_frag = self._AIM_in['general']['irrepairable_res_drift']
          RV_irrep = RandomVariable(ID=-1, dimension_tags=['RED_irrep', ],
                                    distribution_kind='lognormal',
                                    theta=irrep_frag['Median'],
                                    COV=irrep_frag['Sig'] ** 2.
                                    )
          RED_irrep = RV_irrep.sample_distribution(NC_samples)['RED_irrep'].values
  
          # determine if the realizations are repairable
          irrepairable = RED_max > RED_irrep
          irrepairable_IDs = ncID[np.where(irrepairable)[0]]
  
          return irrepairable_IDs
  
      def _calc_repair_cost_and_time(self):
  
          idx = pd.IndexSlice
          DVs = self._AIM_in['decision_variables']
  
          DMG_by_FG_and_DS = self._DMG.groupby(level=[0, 2], axis=1).sum()
  
          repID = self._ID_dict['repairable']
          REP_samples = len(repID)
          DV_COST = pd.DataFrame(np.zeros((REP_samples, len(self._DMG.columns))),
                                 columns=self._DMG.columns, index=repID)
          DV_TIME = deepcopy(DV_COST)
  
          s_fg_keys = sorted(self._FG_dict.keys())
          for fg_id in s_fg_keys:
              FG = self._FG_dict[fg_id]
  
              PG_set = FG._performance_groups
  
              DS_list = self._DMG.loc[:, idx[FG._ID, PG_set[0]._ID, :]].columns
              DS_list = DS_list.levels[2][DS_list.labels[2]].values
  
              for pg_i, PG in enumerate(PG_set):
  
                  PG_ID = PG._ID
  
                  for d_i, d_tag in enumerate(DS_list):
                      dsg_i = int(d_tag[0]) - 1
                      ds_i = int(d_tag[-1]) - 1
  
                      DS = PG._DSG_set[dsg_i]._DS_set[ds_i]
  
                      TOT_qnt = DMG_by_FG_and_DS.loc[repID, (FG._ID, d_tag)]
                      PG_qnt = self._DMG.loc[repID, 
                                             (FG._ID, PG_ID, d_tag)]
  
                      # repair cost
                      if DVs['rec_cost']:
                          COST_samples = DS.unit_repair_cost(quantity=TOT_qnt)
                          DV_COST.loc[:,
                          (FG._ID, PG_ID, d_tag)] = COST_samples * PG_qnt
                      
                      if DVs['rec_time']:
                          # repair time
                          TIME_samples = DS.unit_reconstruction_time(quantity=TOT_qnt)
                          DV_TIME.loc[:,
                          (FG._ID, PG_ID, d_tag)] = TIME_samples * PG_qnt
  
          # sort the columns to enable index slicing later
          if DVs['rec_cost']:
              DV_COST = DV_COST.sort_index(axis=1, ascending=True)
          else:
              DV_COST = None
          if DVs['rec_time']:
              DV_TIME = DV_TIME.sort_index(axis=1, ascending=True)
          else:
              DV_TIME = None
  
          return DV_COST, DV_TIME
  
      def _calc_collapse_injuries(self):
  
          inj_lvls = self._inj_lvls
  
          # calculate casualties and injuries for the collapsed cases
          # generate collapse modes
          colID = self._ID_dict['collapse']
          C_samples = len(colID)
          
          if C_samples > 0:
              
              inj_lvls = self._inj_lvls
              coll_modes = self._AIM_in['collapse_modes']
              P_keys = [cmk for cmk in coll_modes.keys()]
              P_modes = [coll_modes[k]['w'] for k in P_keys]
      
              # create the DataFrame that collects the decision variables
              inj_cols = ['CM',]
              for i in range(inj_lvls):
                  inj_cols.append('INJ-{}'.format(i))
              COL_INJ = pd.DataFrame(np.zeros((C_samples, inj_lvls + 1)),
                                     columns=inj_cols, index=colID)
      
              CM_RV = RandomVariable(ID=-1, dimension_tags=['CM', ],
                                     distribution_kind='multinomial',
                                     p_set=P_modes)
              COL_INJ['CM'] = CM_RV.sample_distribution(C_samples).values
      
              # get the popoulation values corresponding to the collapsed cases
              P_sel = self._POP.loc[colID]
      
              # calculate the exposure of the popoulation
              for cm_i, cmk in enumerate(P_keys):
                  mode_IDs = COL_INJ[COL_INJ['CM'] == cm_i].index
                  CFAR = coll_modes[cmk]['affected_area']
                  INJ = coll_modes[cmk]['injuries']
                  for loc_i in range(len(CFAR)):
                      loc_label = 'LOC{}'.format(loc_i + 1)
                      if loc_label in P_sel.columns:
                          for inj_i in range(inj_lvls):
                              INJ_i = P_sel.loc[mode_IDs, loc_label] * CFAR[loc_i] * \
                                      INJ[inj_i]
                              COL_INJ.loc[mode_IDs, 'INJ-{}'.format(inj_i)] = (
                                  COL_INJ.loc[mode_IDs, 'INJ-{}'.format(inj_i)].add(INJ_i, axis=0).values)
  
              return COL_INJ
          
          else:
              return None
  
      def _calc_non_collapse_injuries(self):
  
          idx = pd.IndexSlice
  
          ncID = self._ID_dict['non-collapse']
          NC_samples = len(ncID)
          DV_INJ_dict = dict([(i, pd.DataFrame(np.zeros((NC_samples,
                                                         len(self._DMG.columns))),
                                               columns=self._DMG.columns,
                                               index=ncID))
                              for i in range(self._inj_lvls)])
          s_fg_keys = sorted(self._FG_dict.keys())
          for fg_id in s_fg_keys:
              FG = self._FG_dict[fg_id]
  
              PG_set = FG._performance_groups
  
              DS_list = self._DMG.loc[:, idx[FG._ID, PG_set[0]._ID, :]].columns
              DS_list = DS_list.levels[2][DS_list.labels[2]].values
  
              for pg_i, PG in enumerate(PG_set):
  
                  PG_ID = PG._ID
  
                  for d_i, d_tag in enumerate(DS_list):
                      dsg_i = int(d_tag[0]) - 1
                      ds_i = int(d_tag[-1]) - 1
  
                      DS = PG._DSG_set[dsg_i]._DS_set[ds_i]
  
                      if DS._affected_area > 0.:
                          P_affected = (self._POP.loc[ncID]
                                        * DS._affected_area / 
                                        self._AIM_in['general']['plan_area'])
  
                          QNT = self._DMG.loc[:, (FG._ID, PG_ID, d_tag)]
  
                          # estimate injuries
                          for i in range(self._inj_lvls):
                              INJ_samples = DS.unit_injuries(severity_level=i,
                                                             sample_size=NC_samples)
                              if INJ_samples is not None:
                                  INJ_samples.index = ncID
                                  P_aff_i = P_affected.loc[:,
                                            'LOC{}'.format(PG._location)]
                                  INJ_i = INJ_samples * P_aff_i * QNT
                                  DV_INJ_dict[i].loc[:,
                                  (FG._ID, PG_ID, d_tag)] = INJ_i
  
                                  # remove the useless columns from DV_INJ
          for i in range(self._inj_lvls):
              DV_INJ = DV_INJ_dict[i]
              DV_INJ_dict[i] = DV_INJ.loc[:, (DV_INJ != 0.0).any(axis=0)]
  
          # sort the columns to enable index slicing later
          for i in range(self._inj_lvls):
              DV_INJ_dict[i] = DV_INJ_dict[i].sort_index(axis=1, ascending=True)
  
          return DV_INJ_dict
  