Merge pull request #10124 from gem/loss_types

Restored multiple loss types in damage calculations

debian/changelog

@@ -14,7 +14,6 @@
   * Raised an error in case of non-invertible hazard curves, affecting
     disaggregation and site-specific hazard spectrum calculations
   * Changed the ruptures.csv exporter to also export the source IDs
-  * Restricted damage calculations to a single loss type
   * Added support for consequence=losses for liquefaction and landslides
   * Added a check for missing secondary perils
   * Added loss types liquefaction and landslide

openquake/baselib/general.py

@@ -431,7 +431,8 @@ def gettemp(content=None, dir=None, prefix="tmp", suffix="tmp", remove=True):
     :param bool remove:
         True by default, meaning the file will be automatically removed
         at the exit of the program
-    :returns: a string with the path to the temporary file
+    :returns:
+        a string with the path to the temporary file
     """
     if dir is not None:
         if not os.path.exists(dir):

openquake/calculators/classical_damage.py

@@ -39,20 +39,21 @@ def classical_damage(riskinputs, param, monitor):
         dictionaries asset_ordinal -> damage(R, L, D)
     """
     crmodel = monitor.read('crmodel')
-    [loss_type] = crmodel.oqparam.loss_types
+    L = crmodel.oqparam.L
     mon = monitor('getting hazard', measuremem=False)
     for ri in riskinputs:
         R = ri.hazard_getter.R
         D = len(crmodel.damage_states)
-        result = AccumDict(accum=numpy.zeros((R, D), F32))
+        result = AccumDict(accum=numpy.zeros((R, L, D), F32))
         with mon:
             haz = ri.hazard_getter.get_hazard()
         for taxo, assets in ri.asset_df.groupby('taxonomy'):
             for rlz in range(R):
                 hcurve = haz[:, rlz]
                 [out] = crmodel.get_outputs(assets, hcurve)
-                for a, frac in zip(assets.ordinal, out[loss_type]):
-                    result[a][rlz] += frac
+                for li, lt in enumerate(crmodel.oqparam.loss_types):
+                    for a, frac in zip(assets.ordinal, out[lt]):
+                        result[a][rlz, li] += frac
         yield result
 
 
@@ -72,12 +73,13 @@ def post_execute(self, result):
             a dictionary asset_ordinal -> array(R, D)
         """
         D = len(self.crmodel.damage_states)
-        damages = numpy.zeros((self.A, self.R, D), numpy.float32)
+        damages = numpy.zeros((self.A, self.R, self.L, D), numpy.float32)
         for a in result:
             damages[a] = result[a]
         self.datastore['damages-rlzs'] = damages
         stats.set_rlzs_stats(self.datastore, 'damages-rlzs',
                              assets=self.assetcol['id'],
+                             loss_type=self.oqparam.loss_types,
                              dmg_state=self.crmodel.damage_states)
         dmg = views.view('portfolio_damage', self.datastore)
         logging.info('\n' + views.text_table(dmg, ext='org'))

openquake/calculators/event_based_damage.py

@@ -50,13 +50,13 @@ class Dparam:
     rng: scientific.MultiEventRNG
 
 
-def zero_dmgcsq(A, R, crmodel):
+def zero_dmgcsq(A, R, L, crmodel):
     """
-    :returns: an array of zeros of shape (A, R, Dc)
+    :returns: an array of zeros of shape (A, R, L, Dc)
     """
     dmg_csq = crmodel.get_dmg_csq()
     Dc = len(dmg_csq) + 1  # damages + consequences
-    return numpy.zeros((A, R, Dc), F32)
+    return numpy.zeros((A, R, L, Dc), F32)
 
 
 def damage_from_gmfs(gmfslices, oqparam, dstore, monitor):
@@ -78,13 +78,13 @@ def damage_from_gmfs(gmfslices, oqparam, dstore, monitor):
     return event_based_damage(df, oqparam, dstore, monitor)
 
 
-def _gen_dd3(asset_df, gmf_df, crmodel, dparam, mon):
-    # yields (aids, dd3) triples 
+def _gen_dd4(asset_df, gmf_df, crmodel, dparam, mon):
+    # yields (aids, dd4) triples 
     # this part is quite slow for discrete damage distributions
     oq = crmodel.oqparam
     E = len(dparam.eids)
     P = len(crmodel.perils)
-    [lt] = oq.loss_types  # assume single loss type
+    L = len(oq.loss_types)
     for taxo, adf in asset_df.groupby('taxonomy'):
         with mon:
             outs = crmodel.get_outputs(adf, gmf_df)  # dicts loss_type -> array
@@ -95,33 +95,36 @@ def _gen_dd3(asset_df, gmf_df, crmodel, dparam, mon):
                 number = assets['value-number']
             else:
                 number = assets['value-number'] = U32(assets['value-number'])
-            dd4 = numpy.zeros((P, A, E, dparam.Dc), F32)
+            dd5 = numpy.zeros((P, A, E, L, dparam.Dc), F32)
             D = dparam.D
             for p, out in enumerate(outs):
-                fractions = out[lt]
-                if oq.float_dmg_dist:
-                    for a in range(A):
-                        dd4[p, a, :, :D] = fractions[a] * number[a]
-                else:
-                    # this is a performance distaster; for instance
-                    # the Messina test in oq-risk-tests becomes 12x
-                    # slower even if it has only 25_736 assets
-                    dd4[p, :, :, :D] = dparam.rng.discrete_dmg_dist(
-                        dparam.eids, fractions, number)
+                for li, lt in enumerate(oq.loss_types):
+                    fractions = out[lt]  # shape (A, E, Dc)
+                    if oq.float_dmg_dist:
+                        for a in range(A):
+                            dd5[p, a, :, li, :D] = fractions[a] * number[a]
+                    else:
+                        # this is a performance distaster; for instance
+                        # the Messina test in oq-risk-tests becomes 12x
+                        # slower even if it has only 25_736 assets
+                        dd5[p, :, :, li, :D] = dparam.rng.discrete_dmg_dist(
+                            dparam.eids, fractions, number)
 
             # compose damage distributions
             if P > 1:
-                dd3 = numpy.empty(dd4.shape[1:])
-                for a in range(A):
-                    for e in range(E):
-                        dd3[a, e, :D] = scientific.compose_dds(dd4[:, a, e, :D])
+                dd4 = numpy.empty(dd5.shape[1:])
+                for li in range(L):
+                    for a in range(A):
+                        for e in range(E):
+                            dd4[a, e, li, :D] = scientific.compose_dds(
+                                dd5[:, a, e, li, :D])
             else:
-                dd3 = dd4[0]
+                dd4 = dd5[0]
             df = crmodel.tmap_df[crmodel.tmap_df.taxi == assets[0]['taxonomy']]
-            csq = crmodel.compute_csq(assets, dd4[:, :, :, :D], df, oq)
-            for name, values in csq.items():
-                dd3[:, :, dparam.csqidx[name]] = values
-        yield aids, dd3  # dd3 has shape (A, E, Dc)
+            csq = crmodel.compute_csq(assets, dd5[:, :, :, :, :D], df, oq)
+            for (name, li), values in csq.items():
+                dd4[:, :, li, dparam.csqidx[name]] = values
+        yield aids, dd4  # dd4 has shape (A, E, L, Dc)
 
 
 def event_based_damage(df, oq, dstore, monitor):
@@ -143,42 +146,38 @@ def event_based_damage(df, oq, dstore, monitor):
         aggids = monitor.read('aggids')
     dmg_csq = crmodel.get_dmg_csq()
     csqidx = {dc: i + 1 for i, dc in enumerate(dmg_csq)}
-    dmgcsq = zero_dmgcsq(len(assetcol), oq.R, crmodel)
-    _A, R, Dc = dmgcsq.shape
+    dmgcsq = zero_dmgcsq(len(assetcol), oq.R, oq.L, crmodel)
+    _A, R, L, Dc = dmgcsq.shape
     D = Dc - len(crmodel.get_consequences())
     rlzs = dstore['events']['rlz_id']
-    dddict = general.AccumDict(accum=numpy.zeros(Dc, F32))  # eid, kid
+    dddict = general.AccumDict(accum=numpy.zeros((L, Dc), F32))  # eid, kid
     for sid, asset_df in assetcol.to_dframe().groupby('site_id'):
         # working one site at the time
         gmf_df = df[df.sid == sid]
         if len(gmf_df) == 0:
             continue
-        oq = crmodel.oqparam
         eids = gmf_df.eid.to_numpy()
         if not oq.float_dmg_dist:
             rng = scientific.MultiEventRNG(
                 oq.master_seed, numpy.unique(eids))
         else:
             rng = None
         dparam = Dparam(eids, aggids, csqidx, D, Dc, rng)
-        for aids, dd3 in _gen_dd3(asset_df, gmf_df, crmodel, dparam, mon):
+        for aids, dd4 in _gen_dd4(asset_df, gmf_df, crmodel, dparam, mon):
+            # dd4 has shape (A, E, L, Dc)
             if R == 1:  # possibly because of collect_rlzs
-                dmgcsq[aids, 0] += dd3.sum(axis=1)
+                dmgcsq[aids, 0] += dd4.sum(axis=1)
             else:
                 for e, rlz in enumerate(rlzs[eids]):
-                    dmgcsq[aids, rlz] += dd3[:, e]
-            tot = dd3.sum(axis=0)  # sum on the assets
+                    dmgcsq[aids, rlz] += dd4[:, e]
+            tot = dd4.sum(axis=0)  # (E, L, Dc)
             for e, eid in enumerate(eids):
                 dddict[eid, oq.K] += tot[e]
                 if oq.K:
                     for kids in dparam.aggids:
                         for a, aid in enumerate(aids):
-                            dddict[eid, kids[aid]] += dd3[a, e]
-    try:
-        [lt] = oq.loss_types
-    except ValueError:
-        lt = oq.total_losses
-    return _dframe(dddict, csqidx, [lt]), dmgcsq
+                            dddict[eid, kids[aid]] += dd4[a, e]
+    return _dframe(dddict, csqidx, oq.loss_types), dmgcsq
 
 
 def _dframe(dddic, csqidx, loss_types):
@@ -190,7 +189,7 @@ def _dframe(dddic, csqidx, loss_types):
             dic['event_id'].append(eid)
             dic['loss_id'].append(scientific.LOSSID[lt])
             for cname, ci in csqidx.items():
-                dic[cname].append(dd[ci])
+                dic[cname].append(dd[li, ci])
     fix_dtypes(dic)
     return pandas.DataFrame(dic)
 
@@ -228,7 +227,7 @@ def execute(self):
             oq.parentdir = os.path.dirname(self.datastore.ppath)
         if oq.investigation_time:  # event based
             self.builder = get_loss_builder(self.datastore, oq)  # check
-        self.dmgcsq = zero_dmgcsq(len(self.assetcol), self.R, self.crmodel)
+        self.dmgcsq = zero_dmgcsq(len(self.assetcol), self.R, oq.L, self.crmodel)
         if oq.K:
             aggids, _ = self.assetcol.build_aggids(
                 oq.aggregate_by, oq.max_aggregations)
@@ -266,7 +265,7 @@ def post_execute(self, dummy):
         """
         oq = self.oqparam
         # no damage check, perhaps the sites where disjoint from gmf_data
-        if self.dmgcsq[:, :, 1:].sum() == 0:
+        if self.dmgcsq[:, :, :, 1:].sum() == 0:
             haz_sids = self.datastore['gmf_data/sid'][:]
             count = numpy.isin(haz_sids, self.sitecol.sids).sum()
             if count == 0:
@@ -284,28 +283,29 @@ def post_execute(self, dummy):
         with prc.datastore:
             prc.run(exports='')
 
-        _A, _R, _Dc = self.dmgcsq.shape
+        _A, _R, L, _Dc = self.dmgcsq.shape
         D = len(self.crmodel.damage_states)
         # fix no_damage distribution for events with zero damage
         number = self.assetcol['value-number']
         for r in range(self.R):
             self.dmgcsq[:, r] /= prc.num_events[r]
-            self.dmgcsq[:, r, 0] = number - self.dmgcsq[:, r, 1:D].sum(axis=1)
+            ndamaged = self.dmgcsq[:, r, :, 1:D].sum(axis=2)  # shape (A, L)
+            for li in range(L):
+                # set no_damage
+                self.dmgcsq[:, r, li, 0] = number - ndamaged[:, li]
 
         assert (self.dmgcsq >= 0).all()  # sanity check
         self.datastore['damages-rlzs'] = self.dmgcsq
         set_rlzs_stats(self.datastore,
                        'damages-rlzs',
                        asset_id=self.assetcol['id'],
-                       rlz=numpy.arange(self.R),
+                       loss_type=oq.loss_types,
                        dmg_state=['no_damage'] + self.crmodel.get_dmg_csq())
 
-        if (hasattr(oq, 'infrastructure_connectivity_analysis')
-                and oq.infrastructure_connectivity_analysis):
-
+        if oq.infrastructure_connectivity_analysis:
             logging.info('Running connectivity analysis')
-            conn_results = connectivity.analysis(self.datastore)
-            self._store_connectivity_analysis_results(conn_results)
+            results = connectivity.analysis(self.datastore)
+            self._store_connectivity_analysis_results(results)
 
     def _store_connectivity_analysis_results(self, conn_results):
         avg_dict = {}
@@ -318,58 +318,41 @@ def _store_connectivity_analysis_results(self, conn_results):
         if 'avg_connectivity_loss_ccl' in conn_results:
             avg_dict['ccl'] = [conn_results['avg_connectivity_loss_ccl']]
         if avg_dict:
-            avg_df = pandas.DataFrame(data=avg_dict)
             self.datastore.create_df(
-                'infra-avg_loss', avg_df,
+                'infra-avg_loss', pandas.DataFrame(data=avg_dict),
                 display_name=DISPLAY_NAME['infra-avg_loss'])
-            logging.info(
-                'Stored avarage connectivity loss (infra-avg_loss)')
         if 'event_connectivity_loss_eff' in conn_results:
             self.datastore.create_df(
                 'infra-event_efl',
                 conn_results['event_connectivity_loss_eff'],
                 display_name=DISPLAY_NAME['infra-event_efl'])
-            logging.info(
-                'Stored efficiency loss by event (infra-event_efl)')
         if 'event_connectivity_loss_pcl' in conn_results:
             self.datastore.create_df(
                 'infra-event_pcl',
                 conn_results['event_connectivity_loss_pcl'],
                 display_name=DISPLAY_NAME['infra-event_pcl'])
-            logging.info(
-                'Stored partial connectivity loss by event (infra-event_pcl)')
         if 'event_connectivity_loss_wcl' in conn_results:
             self.datastore.create_df(
                 'infra-event_wcl',
                 conn_results['event_connectivity_loss_wcl'],
                 display_name=DISPLAY_NAME['infra-event_wcl'])
-            logging.info(
-                'Stored weighted connectivity loss by event (infra-event_wcl)')
         if 'event_connectivity_loss_ccl' in conn_results:
             self.datastore.create_df(
                 'infra-event_ccl',
                 conn_results['event_connectivity_loss_ccl'],
                 display_name=DISPLAY_NAME['infra-event_ccl'])
-            logging.info(
-                'Stored complete connectivity loss by event (infra-event_ccl)')
         if 'taz_cl' in conn_results:
             self.datastore.create_df(
                 'infra-taz_cl',
                 conn_results['taz_cl'],
                 display_name=DISPLAY_NAME['infra-taz_cl'])
-            logging.info(
-                'Stored connectivity loss of TAZ nodes (taz_cl)')
         if 'dem_cl' in conn_results:
             self.datastore.create_df(
                 'infra-dem_cl',
                 conn_results['dem_cl'],
                 display_name=DISPLAY_NAME['infra-dem_cl'])
-            logging.info(
-                'Stored connectivity loss of demand nodes (dem_cl)')
         if 'node_el' in conn_results:
             self.datastore.create_df(
                 'infra-node_el',
                 conn_results['node_el'],
                 display_name=DISPLAY_NAME['infra-node_el'])
-            logging.info(
-                'Stored efficiency loss of nodes (node_el)')

openquake/calculators/export/risk.py

@@ -380,10 +380,13 @@ def export_loss_maps_npz(ekey, dstore):
 
 def modal_damage_array(data, damage_dt):
     # determine the damage state with the highest probability
-    A, _D = data.shape
-    dmgstate = damage_dt.names
-    arr = numpy.zeros(A, [('modal-ds', hdf5.vstr)])
-    arr['modal-ds'] = [dmgstate[data[a].argmax()] for a in range(A)]
+    A, _L, _D = data.shape                          
+    dmgstate = damage_dt['structural'].names
+    arr = numpy.zeros(A, [('modal-ds-' + lt, hdf5.vstr)
+                          for lt in damage_dt.names])
+    for li, loss_type in enumerate(damage_dt.names):
+        arr['modal-ds-' + loss_type] = [dmgstate[data[a, li].argmax()]                  
+                                        for a in range(A)]
     return arr
 
 
@@ -392,9 +395,9 @@ def modal_damage_array(data, damage_dt):
 def export_damages_csv(ekey, dstore):
     oq = dstore['oqparam']
     ebd = oq.calculation_mode == 'event_based_damage'
-    dmg_dt = build_damage_dt(dstore)
     rlzs = dstore['full_lt'].get_realizations()
-    orig = dstore[ekey[0]][:]  # shape (A, R, D)
+    orig = dstore[ekey[0]][:]  # shape (L, A, R, D)
+    dmg_dt = build_damage_dt(dstore)
     writer = writers.CsvWriter(fmt='%.6E')
     assets = get_assets(dstore)
     md = dstore.metadata
@@ -415,13 +418,14 @@ def export_damages_csv(ekey, dstore):
         if ebd:  # export only the consequences from damages-rlzs, i == 0
             rate = len(dstore['events']) * oq.time_ratio / len(rlzs)
             data = orig[:, i] * rate
-            A, Dc = data.shape
+            A, _L, Dc = data.shape
             if Dc == D:  # no consequences, export nothing
                 return []
             csq_dt = build_csq_dt(dstore)
-            damages = numpy.zeros(A, csq_dt)
-            for a in range(A):
-                damages[a] = tuple(data[a, D:Dc])
+            damages = numpy.zeros(A, [(lt, csq_dt) for lt in oq.loss_types])
+            for a in range(A):                                                        
+                for li, lt in enumerate(oq.loss_types):
+                    damages[lt][a] = tuple(data[a, li, D:Dc])
             fname = dstore.build_fname('avg_risk', ros, ekey[1])
         else:  # scenario_damage, classical_damage
             if oq.modal_damage_state: