SimCenter Forum

I am posting a discussion I had in email with Vaafoulay K., with their permission, to preserve it for others to see.

Vaafoulay: Can I obtain the +1 standard deviation costs and downtime directly from the DL summary files?

Adam: If you are interested in the mean+1 standard deviation, I suggest looking at another file, the DL_summary_stats.csv. That file provides statistics on the main outputs, including repair cost and repair time. If you take the value in the mean row and add the value in the std row, you'll get the mean + 1 standard deviation result.

Now, you might want to look at some statistics that are not included in that file. For example, it is quite common to assume that the distribution of the results is lognormal and we are interested in the median + 1 log standard deviation. Such results are not provided in the stats file, so you'd have to calculate them using the sample in the DL_summary.csv.  So, you'd open the DL summary file, take the column of interest (e.g., reconstruction/cost) and use an app of your choice (e.g., Excel, Matlab, Python scripts) to calculate the desired statistic based on the available realizations.

Hi Jiajun,

I am writing to let you know that the calibration notebook is on my list of todos and I'll get to it shortly.

Thank you for your patience.

Adam

Hi Andres,

Apologies for the late response.

I decided to make a few enhancements in Pelicun to streamline the process of adding a new distribution type. I've just released v3.1.b6 and also updated the Example notebook on DesignSafe with a lot of additional details and explanation on how Pelicun 3 works. I encourage you to take a look here: https://www.designsafe-ci.org/data/browser/public/designsafe.storage.published/PRJ-3411v5

Read on if you are still interested in adding the skew normal distribution.

First, you'll need to have the source code available locally. I also recommend linking this version of pelicun to your active Python on the local system to make testing easier. I am not sure how familiar you are with the steps to do this, so let me give you a few tips - please don't hesitate to ask me if you need more information:
- You should get pelicun from the main NHERI-SimCenter account on GitHub. The master branch always provides the latest stable release. Since you are interested in extending the code, you'll need the develop branch; you can find that here: https://github.com/NHERI-SimCenter/pelicun/tree/develop
- I suggest forking the NHERI-SimCenter/pelicun repo and then cloning the develop branch of your own fork to your local hard drive. This will give you the latest version of the code base.
- In your local system, you can store pelicun in any location and make sure Python finds it by adding it to the PYTHONPATH environment variable or by using the externals.pth file. I am happy to provide more information on either of these if you are not familiar with them.
- Once you set things up properly, you should see 3.1.b6 when you import pelicun and check pelicun.__ version__ in a python interpreter.
- Now you can make your edits, test the resulting code and, when you are happy with it, I would appreciate if you committed your contributions to the main repo so that everyone can use this new distribution. To do that, you'll need to first commit your changes to the develop branch on your fork of pelicun and then open a pull request from your develop branch to the develop branch of the main repo. When I see your pull request, I'll test your version and - if everything looks okay - I'll accept it.

That's it, that's how you can extend pelicun and share your work with the community.

Now, let's see how you would go about adding a new distribution. You'll only need to make changes in the uq module to do this. See the latest version of the script here for reference: https://github.com/NHERI-SimCenter/pelicun/blob/develop/pelicun/uq.py . Note that if the uq module gets updated by the time you are reading this, you can always click on history in the top right and go back to today's version so that the line numbers I give below will point to the right location.

Sampling an N dimensional multivariate distribution in Pelicun follows the logic below (see generate_sample method starting at line 1465):
- Sample an N dimensional uniform distribution using Monte Carlo or Latin Hypercube Sampling
- If there are variables across the N dimensions with non-zero correlation, apply the prescribed correlations assuming a Gaussian copula function. First, we try to use a fast Cholesky transformation; if this fails because the correlation matrix is not positive semidefinite, we use another method based on Singular Value Decomposition to apply correlations that preserve as much as possible from the prescribed correlation matrix.
- Perfect correlation can be handled very efficiently with a special 'anchor' feature, but that is almost surely outside of scope for your edits; let me know if you want to know more about it.
- Finally, we take each marginal and use inverse probability integral transformation to transform the sample from uniform distribution to the desired distribution.

Adding a new distribution only affects the last step of the above (as long as you don't want to have some special, non-Gaussian copula function also included - let me know if you do). Below is a list of locations where you'd need to make edits:
- Add the name of the new distribution and a description of its parameters in the documentation at the top of the RandomVariable class - line 791-827; add more info at least under 'distribution' and 'theta' parameters, but potentially also to others, if needed. Currently, Pelicun supports up to three parameters for distributions. That should be sufficient for the skew normal, but let me know if you need more than three.
- If the new distribution requires some special checks, you can add those in the init method (line 829). See, for example, that the multinomial distribution is checked to have probabilities sum up to at most 1.0.
- The cdf() method (line 994) returns the Cumulative Distribution Function ordinates for a given set of x values. You'll need to add your distribution here in an elif clause. Take a look at how the existing distributions are handled and try to mimic the same robustness:
   - You'll see that we get the parameters of the distribution in the theta vector and truncation limits in another vector. Truncation limits can be undefined which should lead to an unbounded distribution in one or both ends and not throw an error.
   - As for the parameters of the distribution, they are typically mandatory, but in some cases we can have rules set up to replace missing values - e.g., see how the limits of the uniform distribution are optional and replaced with infinite values if missing.
   - Make sure the input values are valid inputs to the CDF. For example, 0 and negative numbers are not part of the input domain for the lognormal CDF, so we need to make sure zeros are replaced by the smallest positive number (the computer can handle) in line 1050. This also shows that in general I try to make Pelicun work and handle issues gracefully rather than throwing error messages all the time when something out of the ordinary happens. So, someone feeding a zero to a lognormal CDF will get a probability of 0. I believe in most cases this behavior is preferred over the error message that would terminate execution.
   - I suggest starting with the implementation of the non-truncated version of the distribution. After testing and making sure that it works, you shall expand it by adding the truncation option. Let me know if you need help with this.
- The next method to edit is the inverse_transform() starting at line 1070. Here, you'll need to add an elif clause; I'd add it after line 1144 because the distribution you plan to add belongs to the normal family. The task is similar to the cdf method, but you are implementing an inverse cdf function here:
   - Note that the sample_size argument of the method is only used in special cases. If you implement a skew normal distribution, you should expect to have an array of values (which come from the [0,1] uniform distribution following the sampling logic I explained earlier) that your script will transform to the target distribution.
   - Pay attention to handling missing inputs and truncation limits - the advice I gave for the cdf() method applies here too.
- Finally, you'll need to edit the scale_distribution() method at the top of the uq module (line 70). This method is used by various models in Pelicun to scale input parameters defined in one unit to the SI unit that is used internally. Scaling distributions so far was straightforward as you'll see in the implementation - I hope the skew normal will not be an exception. Note that I define the normal distribution with the mean and coefficient of variation so that the second parameter is unitless and does not need to be scaled. It might be worth pulling similar tricks with the skew normal.

That's it. At this point you have the skew normal distribution implemented and you can sample it by creating a RandomVariable object, adding it to a RandomVariableRegistry and then calling the generate_sample() method of that registry. This is done all over the model module: https://github.com/NHERI-SimCenter/pelicun/blob/develop/pelicun/model.py For example, take a look at the generate_cmp_sample() method starting at line 1165. The _create_cmp_RVs() creates the RandomVariable objects and the registry (line 1139) and then at line 1180 we sample the registry.

Notice that RandomVariable objects are created by feeding the family and parameters of the distribution to the uq module automatically. Using the previous cmp example, take a look at line 1150 in the model module. As long as you feed in the correct name for the family and valid parameters, your new distribution should work immediately without making any changes to the model.py .

The only exception to this is the LossModel - currently, only normal and lognormal random variables are supported for probabilistic loss calculation. Let me know if you want to include the skewed normal there and I can help you set that up.

I hope you'll find the above helpful. Please let me know if you're working on this and don't hesitate to ask questions here if you run into any issues.

Hi Jiajun,

Are you familiar with Google Colab? I could prepare a notebook that shows you an example usage. If not, I can just post a series of commands here that you can use as a template.

As for the environmental impact, please don't feel any pressure about this. I will get the feature added before the end of the year anyway, so the only reason for you to get involved is if you want to have it available before then.

Let me know your preference on Colab and I'll share the template for demand calibration.

Adam

Hi Jiajun,

Do you want to use the Matlab code provided by FEMA P58 or you want to reproduce the method from Appendix G in Pelicun?
That method describes fitting a multivariate lognormal distribution to EDP data, increasing the log standard deviations, and then sampling it. Pelicun can certainly do this if you are interested.

 1 The values in the example are indeed different from the data provided in the tables that you pasted. I wanted to have EDPs for two directions and those tables only provide one value per floor. Furthermore, I could not find in the document what those values are. They might be the average, geometric mean, maximum, or some other combination of the results in the two directions.
The nonlinear analysis results from Figure 1-14 - 1-21 gave me the 10th percentile, median, and 90th percentile of EDPs in two directions on each floor at each intensity level. Those were sufficient information to fit a lognormal distribution and I considered them the raw results of the analysis. Note that these do not match the results in Table 1-35 - 1-42, so those must have been processed in some way.
If you have some suggestions of how to reconcile those differences, please let me know.

2. I have good news: adding environmental impact is high on our list. The framework is ready to support it, but the actual timeline depends on who gets the job done. The worst-case scenario is I have to do it, which means it will be available sometime around the end of this year. I might get some help from graduate students that could expedite this development by a few months.
If you are interested in adding this feature, that's also an option and I would be more than happy to guide you on the way. To give you an idea of what such a development would entail:
- You'd need to create a new consequence function table that describes the parameters of functions used for environmental impacts. This is similar to the repair cost/time table I already have prepared in a CSV file.
- Then, you'd need to add a new child class of the LossModel in the model.py module that handles the calculation of such impacts. Based on my superficial understanding of the environmental impact calculation, you might be able to use most of the BldgRepairModel and perhaps even have the EnvironmentalImpactModel be a child of that class to make the calculation simpler.
That's all we need to get this done. Let me know if you're interested.

Adam

Hi Pooya,

That's a great question.

The offset defined in FEMA P58 for D.30.31.023 is 0. I think this suggests the developers of P58 intended to assign it to the floor of a story - rather than to the ceiling - which makes perfect sense. Let me walk you through two approaches to handle the ground floor and roof cases:

1, Expand the component library
You could add another chiller component to the library (i.e., another row to the fragility_DB_FEMA_P58_2nd.csv file) that has an offset of 1 and otherwise is identical to the existing chiller. You could call this D.30.31.023_roof, for example. Then, if you want to put a chiller on the roof, you would assign this component to the top floor. If you want to put a chiller on any other floor, you would assign the original component to that floor.

2, Assign the component to a non-existing story
Pelicun will not have a problem with you assigning a component to the 5th floor of a 4-story building. You could keep using the original component from FEMA P58 and specify its location as 5. Then, internally, pelicun will see that the component uses PFA as demand and look for a PFA-4 to calculate the damage - which is exactly what you want.

Looking at the two options above, the second one seems more intuitive to me. I am inclined to revise the interpretation of the 'roof' keyword to be parsed as "number of stories + 1". That would make it easy for users to assign such components to the roof without having to worry about the above mechanism inside Pelicun. If you agree, I will get this implemented in the upcoming 3.1.b5

Let me know your thoughts.

Adam

Hi Aram,

I hope you could successfully import the fragility database.

As for the components you asked about earlier, I took a look at the NQET spreadsheet. As far as I understand, the spreadsheet is designed to provide you the quantity of those missing ("None found") component types, but FEMA P58 does not have a corresponding fragility function to use to describe their damage. I assume the spreadsheet and the fragilities were developed by different working groups in the project and not everything lined up perfectly in the end.

You can confirm this by looking at the "Normative Quantity Database" sheet where you'll see those broad component groups listed in Column J. The quantities assigned to various occupancy types are to the right and the components assigned to each group are to the left (column D). Now, the groups you highlighted all have NONE AVAILABLE in column D. Furthermore, if you take a note of the group ID listed in column C for them (that is, D302, D309, and D502, respectively), and go to the "Fragility Database" sheet, you'll see that those IDs are not available in column B.

If this was a consulting project, I would encourage you to look at the literature and try to find a way to model the fragility of those components. Since this is a homework assignment, I suggest checking with your instructor to make sure they're okay with you leaving these out.

Adam

Hi Pooya,

Thank you for these questions; I can assure you, they are very well received and welcome. Please don't hesitate to ask.

The short answer is that you do not need to worry about what the unit specified in the database as long as you use the same type of unit (i.e., length, area, volume, etc.).

The C10.11.001a is a perfect example of a case where I did not follow that rule: I should've used 'ft' for the partition walls. Pelicun cannot convert between area and length units because it does not know what the missing dimension of the partition wall should be. You can see all of the available units in the base.py module of Pelicun here: https://github.com/zsarnoczay/pelicun/blob/develop/pelicun/base.py#L576
I will fix the portion wall units and update the notebook.

In the databases, I kept the units used in the original document to show that a person who develops fragility or consequence functions does not need to worry about what units they used to define those functions because pelicun will automatically perform the conversions in the background.

As for the low-voltage switchgear, the units provided in the FEMA P58 documentation are somewhat misleading. I double-checked this with researchers who were involved in the development of the P58 methodology: all components that have TN, AP, CF, or KV units should be calculated as if they had EA. The rationale is that you would need to know how many of those equipment you have rather than providing a measure of their weight, capacity, or performance. That's why I replaced these measures with EA even in the supporting database.

I hope this helps and thanks for catching that error with the ft2.

Let me know if you have more questions.

Adam