main.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c511dc3c-074e-457d-b2ea-b2be5b3a1d2b",
   "metadata": {},
   "source": [
    "# Main - Charcrete LCA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00fd2deb-66f0-4de4-a91b-03617d9294bd",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 1. Set the project and import Charcrete database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1a0b7cb1-236c-48f1-8cae-c5eba833b87c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# === Switch to kernel bw25! ===\n",
    "\n",
    "# Brightway libraries\n",
    "import bw2data as bd\n",
    "import bw2io as bi\n",
    "import bw2calc as bc\n",
    "import bw2analyzer as ba\n",
    "import bw_processing as bwp\n",
    "\n",
    "# General libraries\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from fs.zipfs import ZipFS\n",
    "import json                        # Library for working with json files\n",
    "from pathlib import Path           # Library for working with paths in different OS     \n",
    "import matplotlib.pyplot as plt    # Library for creating plots\n",
    "import os\n",
    "\n",
    "# Our own functions\n",
    "%run utils.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "01f01d0c-536a-4fa0-9559-fdc815fc2dfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "bd.projects.set_current(\"ei38-teaching-25\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c78cbc52-2e04-493c-88f1-7db9cbc75b33",
   "metadata": {},
   "outputs": [],
   "source": [
    "wdir = os.getcwd()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8080533-ffa5-4061-a834-13bf8bb68ed8",
   "metadata": {},
   "source": [
    "### Importing Charcrete database Excel-file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "66d2ebb1-737c-4eda-aa46-61651dad718d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Extracted 2 worksheets in 0.01 seconds\n",
      "Applying strategy: csv_restore_tuples\n",
      "Applying strategy: csv_restore_booleans\n",
      "Applying strategy: csv_numerize\n",
      "Applying strategy: csv_drop_unknown\n",
      "Applying strategy: csv_add_missing_exchanges_section\n",
      "Applying strategy: normalize_units\n",
      "Applying strategy: normalize_biosphere_categories\n",
      "Applying strategy: normalize_biosphere_names\n",
      "Applying strategy: strip_biosphere_exc_locations\n",
      "Applying strategy: set_code_by_activity_hash\n",
      "Applying strategy: link_iterable_by_fields\n",
      "Applying strategy: assign_only_product_as_production\n",
      "Applying strategy: link_technosphere_by_activity_hash\n",
      "Applying strategy: drop_falsey_uncertainty_fields_but_keep_zeros\n",
      "Applying strategy: convert_uncertainty_types_to_integers\n",
      "Applying strategy: convert_activity_parameters_to_list\n",
      "Applied 16 strategies in 7.10 seconds\n",
      "Applying strategy: link_iterable_by_fields\n",
      "Applying strategy: link_iterable_by_fields\n",
      "7 datasets\n",
      "43 exchanges\n",
      "0 unlinked exchanges\n",
      "  \n",
      "Warning: No valid output stream.\n",
      "Title: Writing activities to SQLite3 database:\n",
      "  Started: 10/28/2022 08:59:53\n",
      "  Finished: 10/28/2022 08:59:53\n",
      "  Total time elapsed: 00:00:00\n",
      "  CPU %: 0.00\n",
      "  Memory %: 0.25\n",
      "Created database: Charcrete_38\n"
     ]
    }
   ],
   "source": [
    "if \"Charcrete_38\" in bd.databases:\n",
    "    del bd.databases[\"Charcrete_38\"]\n",
    "\n",
    "# 1. Specify filepath to your foreground inventories.\n",
    "charcrete_path = wdir+\"/data/Charcrete_database_38.xlsx\"\n",
    "#charcrete_path = \"/home/jupyter-mija.frossard/Charcrete-LCA/Water_database_38.xlsx\"\n",
    "# 2. Create an instance of a class that contains basic methods for importing a database from an excel file.\n",
    "charcrete = bi.ExcelImporter(charcrete_path)  \n",
    "# 3. `apply_strategies` is one of such basic methods, it makes sure units, locations, etc are in correct format.\n",
    "charcrete.apply_strategies()\n",
    "# 4. Next step is to link your foreground exchanges to existing databases by matching relevant exchanges fields.\n",
    "charcrete.match_database(\"biosphere3\", fields=(\"name\", \"unit\", \"categories\"))\n",
    "charcrete.match_database(\"ei 3.8 cutoff\", fields=(\"name\", \"location\", \"unit\"))\n",
    "#charcrete.metadata.pop(None)  # Remove metadata None entry. TODO\n",
    "# 5. If everything is linked, write database so that it is saved in your project.\n",
    "if charcrete.all_linked:\n",
    "    charcrete.write_database()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d3b2f6b6-077d-4e0b-9cca-d3d6f1fc3b99",
   "metadata": {},
   "outputs": [],
   "source": [
    "for unlinked in charcrete.unlinked:\n",
    "    print(unlinked)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5df72dbb-0169-48b3-beab-cc234ddefcb6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Databases dictionary with 7 object(s):\n",
       "\tCharcrete_38\n",
       "\tEXIOBASE 3.8.1 2017 monetary\n",
       "\tEXIOBASE 3.8.1 2017 monetary biosphere\n",
       "\tKorean Motors Chaebol\n",
       "\tWater_38\n",
       "\tbiosphere3\n",
       "\tei 3.8 cutoff"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bd.databases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c54fe20e-eaef-4fef-9303-4570f9e3dc7b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(bd.Database(\"Charcrete_38\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bab9c9c6-1934-49c9-ad6e-870b876e64aa",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. Basic calculations and comparisons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "cbb1175a-0e23-4901-8537-6a84da3881bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "bd.projects.set_current(\"ei38-teaching-25\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1c9c7eb0-6c08-4848-914f-e1c8e0ed1cde",
   "metadata": {},
   "outputs": [],
   "source": [
    "char_db = bd.Database(\"Charcrete_38\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "98760784-3c4e-4bf5-9d2f-41612c33de75",
   "metadata": {},
   "outputs": [],
   "source": [
    "ei_db = bd.Database(\"ei 3.8 cutoff\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "edaa2772-804a-4dc4-9f2b-e58abaed2966",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Charcrete production 5%' (cubic meter, SE, None)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "charcrete = [act for act in char_db if \"Charcrete production 5%\" in act[\"name\"]][0]\n",
    "charcrete_7 = [act for act in char_db if \"Charcrete production 7.5%\" in act[\"name\"]][0]\n",
    "charcrete_10 = [act for act in char_db if \"Charcrete production 10%\" in act[\"name\"]][0]\n",
    "charcrete"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b5307f7f-da4a-4083-a3be-0874a0aaf3d6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['concrete production 20MPa' (cubic meter, RoW, None),\n",
       " 'concrete production 20MPa' (cubic meter, ZA, None),\n",
       " 'concrete production 20MPa' (cubic meter, North America without Quebec, None),\n",
       " 'concrete production 20MPa' (cubic meter, CA-QC, None)]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "concrete_ei = [act for act in ei_db if \"concrete production 20MPa\" in act[\"name\"]]\n",
    "concrete_ei"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57f6a7a8-d933-45b0-a0bd-2193831b55c4",
   "metadata": {},
   "source": [
    "#### Choice of an LCIA method: IPCC 2013 GWP 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "37b7a808-05b6-40fb-8e2e-3453738651c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "ipcc = ('IPCC 2013', 'climate change', 'GWP 100a')\n",
    "method = bd.Method(ipcc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "901dfde0-ad13-447a-903a-6394e6e5cff3",
   "metadata": {},
   "source": [
    "#### Analysis of standard concrete production"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5430e209-74a0-4b50-857d-f5515fa808bb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "'concrete production 20MPa' (cubic meter, RoW, None) 193.8539105051557\n",
      "'concrete production 20MPa' (cubic meter, ZA, None) 229.71912735575796\n",
      "'concrete production 20MPa' (cubic meter, North America without Quebec, None) 217.5524366506636\n",
      "'concrete production 20MPa' (cubic meter, CA-QC, None) 234.18799189342727\n"
     ]
    }
   ],
   "source": [
    "concretes_gwp = [calculateLCA(act, 1, ipcc) for act in concrete_ei]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "64aff395-9d6f-4efc-ac3b-24f601bb7066",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>activity</th>\n",
       "      <th>location</th>\n",
       "      <th>GWP</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>concrete production 20MPa</td>\n",
       "      <td>RoW</td>\n",
       "      <td>193.853911</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>concrete production 20MPa</td>\n",
       "      <td>ZA</td>\n",
       "      <td>229.719127</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>concrete production 20MPa</td>\n",
       "      <td>North America without Quebec</td>\n",
       "      <td>217.552437</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>concrete production 20MPa</td>\n",
       "      <td>CA-QC</td>\n",
       "      <td>234.187992</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                    activity                      location         GWP\n",
       "0  concrete production 20MPa                           RoW  193.853911\n",
       "1  concrete production 20MPa                            ZA  229.719127\n",
       "2  concrete production 20MPa  North America without Quebec  217.552437\n",
       "3  concrete production 20MPa                         CA-QC  234.187992"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "first_lca = pd.DataFrame(list(zip([act[\"name\"] for act in concrete_ei],[act[\"location\"] for act in concrete_ei], concretes_gwp)), columns=[\"activity\", \"location\", \"GWP\"])\n",
    "first_lca"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "1225583b-6443-4db8-ad57-3e1bf8336031",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_concrete_idx = first_lca[['GWP']].idxmin()[0]\n",
    "worst_concrete_idx = first_lca[['GWP']].idxmax()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "b53a46aa-f80c-48f4-b245-6b1ccb0b8154",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fraction of score | Absolute score | Amount | Activity\n",
      "0001 | 193.9 |     1 | 'concrete production 20MPa' (cubic meter, RoW, None)\n",
      "  0.0214 | 4.153 | 858.8 | 'market for gravel, round' (kilogram, CH, None)\n",
      "    0.0118 | 2.283 |  17.4 | 'market for transport, freight, lorry, unspecified' (ton kilometer, RE\n",
      "  0.839 | 162.6 | 256.6 | 'market for cement, portland fly ash cement 21-35%' (kilogram, RoW, No\n",
      "    00.8 |   155 | 256.6 | 'cement production, fly ash 21-35%' (kilogram, RoW, None)\n",
      "      0.0296 | 5.739 | 6.118 | 'market group for electricity, medium voltage' (kilowatt hour, RAS, No\n",
      "      0.146 | 28.27 | 30.21 | 'market for clinker' (kilogram, Europe without Switzerland, None)\n",
      "      0.0303 | 5.868 | 6.363 | 'market for clinker' (kilogram, BR, None)\n",
      "      0.129 | 24.98 | 28.23 | 'market for clinker' (kilogram, IN, None)\n",
      "      0.43 | 83.42 | 90.59 | 'market for clinker' (kilogram, RoW, None)\n",
      "    0.0318 | 6.173 | 44.72 | 'market for transport, freight, lorry, unspecified' (ton kilometer, Ro\n",
      "      0.0133 | 2.573 | 18.16 | 'transport, freight, lorry, all sizes, EURO3 to generic market for tra\n",
      "      0.0105 | 2.045 | 15.35 | 'transport, freight, lorry, all sizes, EURO4 to generic market for tra\n",
      "  0.0167 | 3.237 | 1.752 | 'market for ethylene oxide' (kilogram, RER, None)\n",
      "    0.0167 | 3.237 | 1.752 | 'ethylene oxide production' (kilogram, RER, None)\n",
      "      0.0108 | 2.101 | 1.445 | 'market for ethylene' (kilogram, RER, None)\n",
      "  0.0106 |  2.06 | 2.196 | 'market group for electricity, medium voltage' (kilowatt hour, RAS, No\n",
      "  0.0792 | 15.36 |  1142 | 'market for sand' (kilogram, BR, None)\n",
      "    0.0116 | 2.241 | 8.713 | 'market for transport, freight, lorry 3.5-7.5 metric ton, EURO3' (ton \n",
      "      0.0116 | 2.241 | 8.713 | 'transport, freight, lorry 3.5-7.5 metric ton, EURO3' (ton kilometer, \n",
      "    0.0271 | 5.252 | 40.44 | 'market for transport, freight, lorry 16-32 metric ton, EURO3' (ton ki\n",
      "      0.0271 | 5.252 | 40.44 | 'transport, freight, lorry 16-32 metric ton, EURO3' (ton kilometer, BR\n",
      "    0.0195 | 3.789 | 799.6 | 'sand quarry operation, extraction from river bed' (kilogram, BR, None\n",
      "      0.019 | 3.683 | 39.91 | 'market for diesel, burned in building machine' (megajoule, GLO, None)\n",
      "    0.0124 |  2.41 | 342.7 | 'sand quarry operation, open pit mine' (kilogram, BR, None)\n",
      "      0.0122 | 2.365 | 25.63 | 'market for diesel, burned in building machine' (megajoule, GLO, None)\n"
     ]
    }
   ],
   "source": [
    "ba.print_recursive_calculation(concrete_ei[best_concrete_idx], ipcc, cutoff=0.010)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "d9f09a18-1e10-4c69-9843-0089f39f71ea",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fraction of score | Absolute score | Amount | Activity\n",
      "0001 | 234.2 |     1 | 'concrete production 20MPa' (cubic meter, CA-QC, None)\n",
      "  0.0483 |  11.3 |   950 | 'market for sand' (kilogram, RoW, None)\n",
      "    0.0181 | 4.235 | 30.68 | 'market for transport, freight, lorry, unspecified' (ton kilometer, Ro\n",
      "    0.0121 | 2.824 | 672.1 | 'gravel and sand quarry operation' (kilogram, RoW, None)\n",
      "  0.0477 | 11.17 | 960.1 | 'market for gravel, round' (kilogram, RoW, None)\n",
      "    0.0172 | 4.034 | 960.1 | 'gravel and sand quarry operation' (kilogram, RoW, None)\n",
      "    0.0159 | 3.714 | 26.91 | 'market for transport, freight, lorry, unspecified' (ton kilometer, Ro\n",
      "  0.878 | 205.6 | 239.5 | 'market for cement, Portland' (kilogram, CA-QC, None)\n",
      "    0.0151 | 3.534 |  25.6 | 'market for transport, freight, lorry, unspecified' (ton kilometer, Ro\n",
      "    0.855 | 200.2 | 239.5 | 'cement production, Portland' (kilogram, CA-QC, None)\n",
      "      0.85 |   199 | 220.4 | 'clinker production' (kilogram, CA-QC, None)\n"
     ]
    }
   ],
   "source": [
    "ba.print_recursive_calculation(concrete_ei[worst_concrete_idx], ipcc, cutoff=0.010)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40ea5e88-a5cd-4979-992e-2354e0a1e596",
   "metadata": {},
   "source": [
    "Concrete with the least impact has \"market for cement, portland fly ash cement 21-35%, RoW\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "384c4b51-64cd-4d72-aa78-cf9227484df5",
   "metadata": {},
   "outputs": [],
   "source": [
    "concrete_ei = [*concrete_ei, charcrete, charcrete_7, charcrete_10]\n",
    "first_lca = pd.DataFrame(list(zip([act[\"name\"] for act in concrete_ei],[act[\"location\"] for act in concrete_ei], concretes_gwp)), columns=[\"activity\", \"location\", \"GWP\"])\n",
    "first_lca\n",
    "first_lca.to_csv(wdir+\"/outputs/first_lca_results_concrete_1_m3\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc6a02f4-058b-4f43-8f4e-36b8878885a0",
   "metadata": {},
   "source": [
    "#### Analysis over the cement production"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "cdbac957-0fb4-4bfa-b08b-343c89005398",
   "metadata": {},
   "outputs": [],
   "source": [
    "cement = [act for act in ei_db if 'market for cement,' in act[\"name\"] and \"kilogram\" in act[\"unit\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "12338791-69ed-41bf-9e44-2b1a9c07c2d7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>activity</th>\n",
       "      <th>location</th>\n",
       "      <th>GWP</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>market for cement, pozzolana and fly ash 5-15%</td>\n",
       "      <td>RoW</td>\n",
       "      <td>0.848432</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>market for cement, alternative constituents 21...</td>\n",
       "      <td>CH</td>\n",
       "      <td>0.616254</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>market for cement, pozzolana and fly ash 36-55%</td>\n",
       "      <td>RoW</td>\n",
       "      <td>0.600784</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>market for cement, pozzolana and fly ash 11-35%</td>\n",
       "      <td>Europe without Switzerland</td>\n",
       "      <td>0.714784</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>market for cement, blast furnace slag 81-95%</td>\n",
       "      <td>Europe without Switzerland</td>\n",
       "      <td>0.261332</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>93</th>\n",
       "      <td>market for cement, unspecified</td>\n",
       "      <td>PE</td>\n",
       "      <td>0.771636</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>94</th>\n",
       "      <td>market for cement, pozzolana and fly ash 11-35%</td>\n",
       "      <td>CH</td>\n",
       "      <td>0.603359</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>95</th>\n",
       "      <td>market for cement, blast furnace slag 66-80%</td>\n",
       "      <td>RoW</td>\n",
       "      <td>0.443241</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>96</th>\n",
       "      <td>market for cement, unspecified</td>\n",
       "      <td>Europe without Switzerland</td>\n",
       "      <td>0.788553</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>97</th>\n",
       "      <td>market for cement, blast furnace slag 25-70%</td>\n",
       "      <td>US</td>\n",
       "      <td>0.599254</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>98 rows × 3 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                             activity  \\\n",
       "0      market for cement, pozzolana and fly ash 5-15%   \n",
       "1   market for cement, alternative constituents 21...   \n",
       "2     market for cement, pozzolana and fly ash 36-55%   \n",
       "3     market for cement, pozzolana and fly ash 11-35%   \n",
       "4        market for cement, blast furnace slag 81-95%   \n",
       "..                                                ...   \n",
       "93                     market for cement, unspecified   \n",
       "94    market for cement, pozzolana and fly ash 11-35%   \n",
       "95       market for cement, blast furnace slag 66-80%   \n",
       "96                     market for cement, unspecified   \n",
       "97       market for cement, blast furnace slag 25-70%   \n",
       "\n",
       "                      location       GWP  \n",
       "0                          RoW  0.848432  \n",
       "1                           CH  0.616254  \n",
       "2                          RoW  0.600784  \n",
       "3   Europe without Switzerland  0.714784  \n",
       "4   Europe without Switzerland  0.261332  \n",
       "..                         ...       ...  \n",
       "93                          PE  0.771636  \n",
       "94                          CH  0.603359  \n",
       "95                         RoW  0.443241  \n",
       "96  Europe without Switzerland  0.788553  \n",
       "97                          US  0.599254  \n",
       "\n",
       "[98 rows x 3 columns]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#cement_gwp = [calculateLCA(act, 1, ipcc) for act in cement]\n",
    "fu, data_objs, _ = bd.prepare_lca_inputs({cement[0]: 1}, method=ipcc)\n",
    "lca = bc.LCA(fu, data_objs=data_objs)\n",
    "lca.lci()\n",
    "lca.lcia()\n",
    "results_ciment = []\n",
    "for act in cement:\n",
    "    lca.redo_lcia({act.id: 1})\n",
    "    results_ciment.append((act[\"name\"], act[\"location\"], lca.score))\n",
    "\n",
    "cements_lca = pd.DataFrame(results_ciment, columns=[\"activity\", \"location\", \"GWP\"])\n",
    "cements_lca.to_csv(wdir+\"/outputs/gwp_results_cement_1_kg.csv\")\n",
    "cements_lca"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "6272bf6f-471d-4a63-987b-41140de4261f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "activity    32\n",
       "location    10\n",
       "GWP         97\n",
       "dtype: int64"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cements_lca.nunique()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c85170ed-c5da-40fa-9d90-dd1ef11edad7",
   "metadata": {},
   "source": [
    "#### Adding the charcrete to the comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "a8ab9468-19dc-4cfe-a01c-d8c71e1bc37c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['concrete production 20MPa' (cubic meter, RoW, None), 'concrete production 20MPa' (cubic meter, ZA, None), 'concrete production 20MPa' (cubic meter, North America without Quebec, None), 'concrete production 20MPa' (cubic meter, CA-QC, None), 'Charcrete production 5%' (cubic meter, SE, None), 'Charcrete production 7.5%' (cubic meter, SE, None), 'Charcrete production 10%' (cubic meter, SE, None), 'Charcrete production 5%' (cubic meter, SE, None), 'Charcrete production 7.5%' (cubic meter, SE, None), 'Charcrete production 10%' (cubic meter, SE, None)]\n",
      "'concrete production 20MPa' (cubic meter, RoW, None) 193.8539105051557\n",
      "'concrete production 20MPa' (cubic meter, ZA, None) 229.71912735575796\n",
      "'concrete production 20MPa' (cubic meter, North America without Quebec, None) 217.5524366506636\n",
      "'concrete production 20MPa' (cubic meter, CA-QC, None) 234.18799189342727\n",
      "'Charcrete production 5%' (cubic meter, SE, None) 177.64562302249578\n",
      "'Charcrete production 7.5%' (cubic meter, SE, None) 173.55581976132007\n",
      "'Charcrete production 10%' (cubic meter, SE, None) 169.3747240667536\n",
      "'Charcrete production 5%' (cubic meter, SE, None) 177.64562302249578\n",
      "'Charcrete production 7.5%' (cubic meter, SE, None) 173.55581976132007\n",
      "'Charcrete production 10%' (cubic meter, SE, None) 169.3747240667536\n"
     ]
    }
   ],
   "source": [
    "concrete_ei = [*concrete_ei, charcrete, charcrete_7, charcrete_10]\n",
    "print(concrete_ei)\n",
    "concretes_gwp = [calculateLCA(act, 1, ipcc) for act in concrete_ei]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "1ae61a45-2914-4df4-ab89-54b38bfc2802",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>activity</th>\n",
       "      <th>location</th>\n",
       "      <th>GWP</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>concrete production 20MPa</td>\n",
       "      <td>RoW</td>\n",
       "      <td>193.853911</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>concrete production 20MPa</td>\n",
       "      <td>ZA</td>\n",
       "      <td>229.719127</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>concrete production 20MPa</td>\n",
       "      <td>North America without Quebec</td>\n",
       "      <td>217.552437</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>concrete production 20MPa</td>\n",
       "      <td>CA-QC</td>\n",
       "      <td>234.187992</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Charcrete production 5%</td>\n",
       "      <td>SE</td>\n",
       "      <td>177.645623</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Charcrete production 7.5%</td>\n",
       "      <td>SE</td>\n",
       "      <td>173.555820</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Charcrete production 10%</td>\n",
       "      <td>SE</td>\n",
       "      <td>169.374724</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Charcrete production 5%</td>\n",
       "      <td>SE</td>\n",
       "      <td>177.645623</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Charcrete production 7.5%</td>\n",
       "      <td>SE</td>\n",
       "      <td>173.555820</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Charcrete production 10%</td>\n",
       "      <td>SE</td>\n",
       "      <td>169.374724</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                    activity                      location         GWP\n",
       "0  concrete production 20MPa                           RoW  193.853911\n",
       "1  concrete production 20MPa                            ZA  229.719127\n",
       "2  concrete production 20MPa  North America without Quebec  217.552437\n",
       "3  concrete production 20MPa                         CA-QC  234.187992\n",
       "4    Charcrete production 5%                            SE  177.645623\n",
       "5  Charcrete production 7.5%                            SE  173.555820\n",
       "6   Charcrete production 10%                            SE  169.374724\n",
       "7    Charcrete production 5%                            SE  177.645623\n",
       "8  Charcrete production 7.5%                            SE  173.555820\n",
       "9   Charcrete production 10%                            SE  169.374724"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "first_lca = pd.DataFrame(list(zip([act[\"name\"] for act in concrete_ei],[act[\"location\"] for act in concrete_ei], concretes_gwp)), columns=[\"activity\", \"location\", \"GWP\"])\n",
    "first_lca.to_csv(wdir+\"/outputs/first_lca_results_concrete_1_m3.csv\")\n",
    "first_lca"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "f3ebc5e7-c433-41bd-adb1-d4a0a1e91575",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fraction of score | Absolute score | Amount | Activity\n",
      "0001 | 177.6 |     1 | 'Charcrete production 5%' (cubic meter, SE, None)\n",
      "  0.0123 | 2.184 |  0.01 | 'Biochar production' (ton, SE, None)\n",
      "  0.0346 | 6.151 |  1000 | 'market for gravel, crushed' (kilogram, CH, None)\n",
      "    0.015 | 2.659 | 20.26 | 'market for transport, freight, lorry, unspecified' (ton kilometer, RE\n",
      "    0.0194 | 3.437 |  1000 | 'gravel production, crushed' (kilogram, CH, None)\n",
      "  0.0255 | 4.536 |   938 | 'market for sand' (kilogram, CH, None)\n",
      "    0.014 | 2.494 |    19 | 'market for transport, freight, lorry, unspecified' (ton kilometer, RE\n",
      "    0.0112 | 1.991 |   938 | 'gravel and sand quarry operation' (kilogram, CH, None)\n",
      "  0.024 | 4.261 |   3.2 | 'market for plasticiser, for concrete, based on sulfonated melamine fo\n",
      "    0.0231 |   4.1 |   3.2 | 'plasticiser production, for concrete, based on sulfonated melamine fo\n",
      "      0.0121 | 2.147 | 1.075 | 'market for chemical, organic' (kilogram, GLO, None)\n",
      "  0.92 | 163.5 | 186.8 | 'market for cement, Portland' (kilogram, Europe without Switzerland, N\n",
      "    0.0148 |  2.62 | 19.97 | 'market for transport, freight, lorry, unspecified' (ton kilometer, RE\n",
      "    0.898 | 159.5 | 186.8 | 'cement production, Portland' (kilogram, Europe without Switzerland, N\n",
      "      0.0154 | 2.744 | 7.024 | 'market group for electricity, medium voltage' (kilowatt hour, Europe \n",
      "      0.878 |   156 | 168.6 | 'clinker production' (kilogram, Europe without Switzerland, None)\n"
     ]
    }
   ],
   "source": [
    "ba.print_recursive_calculation(charcrete, ipcc, cutoff=0.010)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "952e4578-f5c9-4a8b-b3f0-483bb80f4c56",
   "metadata": {},
   "source": [
    "**Charcrete** is not that good I guess... Need to investigate the formulation of the concrete and other parameters with more data/uncertainty analysis on : \n",
    "- cement production process and localisation\n",
    "- formulation of cement/biochar/other..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf41cd9b-9aed-4d27-a05c-7a43ad565a2e",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Replacing stuff in the activity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "1611b390-3c56-46df-a615-23ef4cf5775c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup\n",
    "concrete_market = bd.get_activity(database=\"ei 3.8 cutoff\", name=\"market group for concrete, normal\", location = \"GLO\")\n",
    "concrete_RoW = bd.get_activity(database=\"ei 3.8 cutoff\", name='market for concrete, normal', location = 'RoW')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "f1e6ad90-0f37-4677-b2d8-fd2ef75ba599",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create missing column\n",
    "concrete_market.new_edge(input=charcrete, amount=0, type=\"technosphere\").save() # declare dependency on new database by drawing edges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "85c81825-9baf-478e-b5df-2d5f6f552752",
   "metadata": {},
   "outputs": [],
   "source": [
    "RoW_results = replace_stuff(concrete_market, concrete_RoW, charcrete, ipcc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "a984a530-7a58-4b94-bef4-ac8e1227c847",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGxCAYAAACEFXd4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+0klEQVR4nO3deVwW5f7/8fetLCICgguLEmJuJW6puZW4pVluZWnZKW3xeHI7pKV5SsVSKcvtaNnmWpn+yq3MVFxQO2guZa6V5X4UUVMQZVG4fn/45T7eAsqNN0HT6/l4zOPhfc01c3/mZhjezlwzt80YYwQAAGBRJYq6AAAAgMJE2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2IHLzJkzRzabTTabTXFxcTnmG2NUrVo12Ww2tWrVqlBriYuLy7MOV8re5sOHDzu9bHx8vKKjo3X+/HmX19WnTx9VqVLF5estDNHR0bLZbDpz5swtr+vSpUuKjo52+ue+adMmeXp66siRI/a2Vq1a2fdnm82mUqVK6c4779TYsWOVkZHhdG2dOnWSj4+Prly54tD+ww8/yGazKTg4ONe6bDab/v3vfzv9foXlq6++UufOnRUYGCgPDw8FBASobdu2+vTTT3X58mVJ0uHDh2Wz2fT2228XcbWuNX78eC1durTAy587d05ly5a9pXWgYAg7cDkfHx/NnDkzR/uGDRv022+/ycfHp9BruOuuu7R582bdddddhf5eBRUfH68xY8YUStj5q7p06ZLGjBnjVNgxxigqKkp9+/ZVWFiYw7yqVatq8+bN2rx5sz7//HNVr15dI0eO1MCBA52urXXr1kpJSdH27dsd2uPi4uTt7a2EhAT99NNPOeZlL1vUjDF6+umn1aVLF2VlZWnSpElas2aN5s6dq3r16ql///569913i7rMQnWrYcff318vvPCCXnrppQIFZhQcYQcu17NnTy1atEjJyckO7TNnzlSzZs102223FXoNvr6+atq0qXx9fQv9vVAwqampKg5fzbdy5Up9//33GjRoUI55Xl5eatq0qZo2barOnTtr0aJFql69uubOnau0tDSn3ic7sFwfxOLi4tS1a1cFBwdr/fr1OeaVL19eERERzm3U/zlw4IDLPuO33npLc+bM0ZgxY/T111/riSeeUMuWLdW5c2dNmjRJP//8sxo2bOiS98qvW92HLl265MJq8ucf//iHDh8+rC+++OIPf++/MsIOXO7xxx+XJH322Wf2tqSkJC1atEjPPPNMrsuMGTNGTZo0UUBAgHx9fXXXXXdp5syZOQ5kVapUUadOnbRy5Urddddd8vLyUq1atTRr1iyHfrldxurTp4/KlCmjvXv3qm3btvL29laFChU0cOBAh4Ne9in4OXPm5KjTZrMpOjr6htsfGxurrl27qnLlyipVqpSqVaumfv36OVymiY6O1ksvvSRJCg8Pz/Xy38KFC9WsWTN5e3urTJky6tChg3744Ycc7zdnzhzVrFlTnp6euuOOOzRv3rwb1net7M9z+fLlatCggby8vHTHHXdo+fLl9nXfcccd8vb21t13353jrMT27dv12GOPqUqVKvLy8lKVKlX0+OOPO1wOyl6PzWbT6tWr9cwzz6hChQoqXbq00tPTc63rp59+UtWqVdWkSRMlJiZKkhISEtSvXz9VrlxZHh4eCg8P15gxY+yXhQ4fPqwKFSpIuro/ZX+mffr0ueFnMGPGDDVu3Fg1a9a86efl5uam+vXrKyMjw+GMXFpamkaMGKHw8HB5eHioUqVKGjBggEOf+vXry9/f3+FnnJWVpU2bNqlVq1aKjIx0CDsZGRnavHmz/XJaQWSfrRo2bFiu+05+Xb58WW+++aZq1aqlkSNH5tonKChI99xzT472SZMmKTw8XGXKlFGzZs20ZcsWh/mu2ofmz5+vZs2aqUyZMipTpozq16/vcIa5VatWioiI0MaNG9W8eXOVLl3afjxKTk7Wiy++6PDzi4qK0sWLF+3L22w2Xbx4UXPnzrXvW9dejr/Z/pktMDBQ9913n9577718fPJwGQO4yOzZs40ks23bNvPkk0+au+++2z5vxowZxtvb2yQnJ5vatWubyMhIh2X79OljZs6caWJjY01sbKx5/fXXjZeXlxkzZoxDv7CwMFO5cmVz5513mnnz5plVq1aZRx991EgyGzZssPdbv369kWTWr19vb+vdu7fx8PAwt912mxk3bpxZvXq1iY6ONm5ubqZTp072focOHTKSzOzZs3NsoyQzevToHNt86NAhh22NiYkxX375pdmwYYOZO3euqVevnqlZs6bJyMgwxhhz7NgxM2jQICPJLF682GzevNls3rzZJCUlGWOMGTdunLHZbOaZZ54xy5cvN4sXLzbNmjUz3t7eZu/evTnev2vXruarr74yn3zyialWrZoJDQ01YWFhN/uR2T/PiIgI89lnn5kVK1aYJk2aGHd3dzNq1CjTokULs3jxYrNkyRJTo0YNExgYaC5dumRf/vPPPzejRo0yS5YsMRs2bDALFiwwkZGRpkKFCub06dM56qxUqZL5+9//br755hvzxRdfmCtXrpjRo0cbSfb+cXFxxt/f33Tt2tVcvHjRGGPMyZMn7dv0/vvvmzVr1pjXX3/deHp6mj59+hhjjElLSzMrV640ksyzzz5r/0x//fXXPLc/PT3deHl5mWHDhuWYFxkZaWrXrp2jvVGjRqZs2bLmypUrxhhjsrKyTIcOHYybm5sZOXKkWb16tXn77beNt7e3adCggUlLS7Mv27VrV+Pt7W0uX75sjDFmx44dRpL5+eefzYwZM0zFihXtfTds2GAkmXfeeeemP8e87N2714wYMcJUrVrVSDK1atUyY8aMMb/88otT64mPjzeSzPDhw/PVP/t3qEqVKub+++83S5cuNUuXLjV16tQx/v7+5vz58/a+rtiHRo4caSSZhx9+2Hz++edm9erVZtKkSWbkyJH25SMjI01AQIAJDQ0106ZNM+vXrzcbNmwwFy9eNPXr1zfly5c3kyZNMmvWrDFTp041fn5+pk2bNiYrK8sYY8zmzZuNl5eXeeCBB+z7VvbvYn72z2u9+eabpkSJEubcuXNO/RxQcIQduMy1YSc7bOzZs8cYY0zjxo3tv/S5hZ1rZWZmmsuXL5vXXnvNlCtXzn6wMebqH+dSpUqZI0eO2NtSU1NNQECA6devn70tr7AjyUydOtXh/caNG2ckmW+//dYYc+th51pZWVnm8uXL5siRI0aSWbZsmX3eW2+9leuyR48eNW5ubmbQoEEO7RcuXDBBQUGmR48e9s8pJCTE3HXXXQ6f0eHDh427u3u+w46Xl5c5fvy4vW3nzp1GkgkODraHDWOMWbp0qZFkvvzyyzzXd+XKFZOSkmK8vb0dPufsz+mpp57Kscy1Yefjjz82Hh4eZvDgwSYzM9Pep1+/fqZMmTIOP3djjHn77beNJPsfndOnT+f4Gd3Id999ZySZBQsW5JiXHXYuX75sLl++bE6ePGlGjRplJJn33nvP3i87YE2YMMFh+YULFxpJ5oMPPrC3TZkyxUgy8fHxxhhjJk6caIKDg40xxuzbt8/hd2bMmDFGktm3b1++tuVmtm7daoYMGWIqV65sJJmGDRuaiRMnOvzs87JgwYIc230j2b9DderUsYfC7Bokmc8++yzPZZ3dhw4ePGhKlixpnnjiiRvWFBkZaSSZtWvXOrTHxMSYEiVKmG3btjm0f/HFF0aSWbFihb3N29vb9O7dO8e687t/ZouNjTWSzDfffHPDmuE6XMZCoYiMjNTtt9+uWbNmaffu3dq2bVuel7Akad26dWrXrp38/PxUsmRJubu7a9SoUTp79qz9Mka2+vXrO4z7KVWqlGrUqJHjtHdennjiCYfXvXr1kqQc4yUKKjExUf/4xz8UGhoqNzc3ubu72we+7t+//6bLr1q1SleuXNFTTz2lK1eu2KdSpUopMjLSfhnk559/1okTJ9SrVy+HyxxhYWFq3rx5vuutX7++KlWqZH99xx13SLp62r906dI52q/9nFNSUjR8+HBVq1ZNbm5ucnNzU5kyZXTx4sVct7V79+551jFu3Dj16dNHb7zxhqZOnaoSJf53eFq+fLlat26tkJAQh8+kY8eOkq4Ofi+IEydOSJIqVqyY6/y9e/fK3d1d7u7uCg4O1muvvaYRI0aoX79+9j7r1q2TpByXyx599FF5e3tr7dq19rbrx+3ExcUpMjJS0tXPt2LFivb9MC4uToGBgfbPPS+ZmZkOn0lWVlau/Ro3bqyJEyfq6NGj2rhxo5o2bao333xTt912m4YOHXrD9yioBx98UCVLlrS/rlu3riTX7kOxsbHKzMzUgAEDblqPv7+/2rRp49C2fPlyRUREqH79+g6fY4cOHfJ9R6ez+2f2/vbf//73puuGa7gVdQGwJpvNpqefflr//ve/lZaWpho1aujee+/Nte/WrVvVvn17tWrVSh9++KH9mvfSpUs1btw4paamOvQvV65cjnV4enrm6JcbNze3HMsHBQVJks6ePZvfzctTVlaW2rdvrxMnTmjkyJGqU6eOvL29lZWVpaZNm+arxlOnTkm6+scpN9khILve7PqvFRQUlO/b4QMCAhxee3h43LD92oG5vXr10tq1azVy5Eg1btxYvr6+stlseuCBB3Ld1txur872ySefqFKlSnrsscdyzDt16pS++uorubu757psQW9bz66xVKlSuc6//fbbtWDBAhljdOTIEY0dO1YxMTGqW7euvc6zZ8/Kzc3NPl4om81mU1BQkMN+VadOHZUvX17r16/X8OHDtWnTJsXExNjnt2zZUnFxcerbt682b96sbt263XQb2rZt6/DHtHfv3rmON8t2+fJlJSUl6fz580pNTZWHh8dN75DM/s/FoUOHblrPta7/XfP09JQkh33jVveh06dPS5IqV65803py2/9OnTqlX3/99Zb2LWf3z+z9LT/HA7gGYQeFpk+fPho1apTee+89jRs3Ls9+CxYskLu7u5YvX+7wR6cwnkVx5coVnT171uEgnJCQIOl/B+bsGq4fPJufMLRnzx79+OOPmjNnjnr37m1v//XXX/NdY/ny5SVJX3zxRY5boa+VXW92/dfKrc3VkpKStHz5co0ePVovv/yyvT09PV2///57rsvcaKDtypUr1bNnT917771au3atw7aXL19edevWzXM/CgkJKdA2ZH/WedVbqlQpNWrUSNLV8Nm6dWvVrl1bUVFR6tSpk8qUKaNy5crpypUrOn36tEPgMcYoISHBIbTabDZFRkZq5cqV2rp1q86fP28/syNdPSMaHR2tzZs3Ky0tLV+3nL///vu6cOFCjm261pUrV7R27VotXLhQS5Ys0YULF9SmTRtNnTpV3bt3v+ldi40aNVJAQICWLVummJiYAg+Yvp4r9qHsz/z48eMKDQ294fvlVnf58uXl5eWV4yaHa+ffjLP7Z/a25WfdcA3CDgpNpUqV9NJLL+mnn35y+MN/PZvNJjc3N4fT3ampqfr4448Lpa5PP/1UgwcPtr+eP3++JNnvrAgMDFSpUqW0a9cuh+WWLVt203VnH0yz/web7f3338/RN7f/5UpShw4d5Obmpt9+++2Gl31q1qyp4OBgffbZZxoyZIj9vY8cOaL4+PgCB4D8stlsMsbk2NaPPvpImZmZTq8vLCxMmzZtUrt27eyBp3r16pKuPpBvxYoVuv322+Xv75/nOvL6TPOSfYnot99+y1f/cuXK6Y033tDTTz+tadOmacSIEWrbtq0mTJigTz75RC+88IK976JFi3Tx4kW1bdvWYR2tW7fWokWL9NZbb6lixYoOl6kiIyN19uxZTZs2zd73Zm50F9nmzZs1d+5cLVq0SGfOnFHTpk01ZswY9ezZU4GBgfnaZklyd3fX8OHDNXz4cL3++usaNWpUjj6JiYk6cOCAWrRoke/1umIfat++vUqWLKkZM2aoWbNm+X7vbJ06ddL48eNVrlw5hYeH37BvXmeQ87t/Zjt48KAk6c4773S6XhQMYQeF6o033rhpnwcffFCTJk1Sr1699Pe//11nz57V22+/neMA6AoeHh6aOHGiUlJS1LhxY8XHx2vs2LHq2LGj/bZZm82mv/3tb5o1a5Zuv/121atXT1u3brWHohupVauWbr/9dr388ssyxiggIEBfffWVYmNjc/StU6eOJGnq1Knq3bu33N3dVbNmTVWpUkWvvfaaXnnlFR08eFD333+//P39derUKW3dulXe3t4aM2aMSpQooddff13PPfecHnroIfXt21fnz59XdHR0rpe2XM3X11ctW7bUW2+9pfLly6tKlSrasGGDZs6cqbJlyxZoncHBwdqwYYM6dOigli1bKjY2VhEREXrttdcUGxur5s2ba/DgwapZs6bS0tJ0+PBhrVixQu+9954qV64sHx8fhYWFadmyZWrbtq0CAgLsteWmcuXKqlq1qrZs2eIQgG/kqaee0qRJk/T2229rwIABuu+++9ShQwcNHz5cycnJatGihXbt2qXRo0erQYMGevLJJx2Wzw4wS5Ys0SOPPOIwLyIiQuXKldOSJUtUqVIle9grqBEjRujMmTOKiopSr169bvrH/EZeeukl7d+/X6NHj9bWrVvVq1cvhYaGKikpSRs3btQHH3ygMWPGOBV2XLEPValSRf/617/0+uuvKzU1VY8//rj8/Py0b98+nTlzRmPGjLnh8lFRUVq0aJFatmypF154QXXr1lVWVpaOHj2q1atXa+jQoWrSpImkq7+zcXFx+uqrrxQcHCwfHx/VrFkz3/tnti1btqhcuXL2YwD+AEU6PBqWcu3dWDeS291Ys2bNMjVr1jSenp6matWqJiYmxsycOTPH3UphYWHmwQcfzLHOyMhIh3XmdTeWt7e32bVrl2nVqpXx8vIyAQEB5vnnnzcpKSkO60tKSjLPPfecCQwMNN7e3qZz587m8OHD+boba9++fea+++4zPj4+xt/f3zz66KPm6NGjud4lNGLECBMSEmJKlCiRo96lS5ea1q1bG19fX+Pp6WnCwsLMI488YtasWeOwjo8++shUr17deHh4mBo1aphZs2aZ3r175/turNw+T0lmwIABDm3Zd9i89dZb9rbjx4+b7t27G39/f+Pj42Puv/9+s2fPHhMWFuZw18qN9o3rbz03xpjz58+bFi1amICAAPsyp0+fNoMHDzbh4eHG3d3dBAQEmIYNG5pXXnnF4ee3Zs0a06BBA+Pp6Wkk5Xr3zLVGjhxp/P39HW4RNybvW8+NMebrr782kuyPRkhNTTXDhw83YWFhxt3d3QQHB5vnn38+z1uLg4KCjCQzffr0HPO6detmJN307qL8OHHixC2v43rLli0zDz74oKlQoYJxc3Mz/v7+pnXr1ua9994z6enpxpjc95Vs1/8euGIfMsaYefPmmcaNG5tSpUqZMmXKmAYNGjjcUXmjn2dKSop59dVXTc2aNY2Hh4fx8/MzderUMS+88IJJSEiw99u5c6dp0aKFKV26tJHkcMzJ7/6ZlZVlwsLCctxticJlM6YYPMIU+AP06dNHX3zxhVJSUoq6FBQjJ06cUHh4uObNm6eePXsWdTmwuLVr16p9+/bau3evatWqVdTl/GVw6zmAv7SQkBBFRUVp3Lhxed62DbjK2LFj9cwzzxB0/mCM2QHwl/fqq6+qdOnS+u9//3vTO3qAgjp37pwiIyPVv3//oi7lL4fLWAAAwNK4jAUAACyNsAMAACyNsAMAACyNAcq6+n1GJ06ckI+Pj8segw4AAAqXMUYXLlxQSEiIw5cHX4+wo6vP2eAODAAA/pyOHTt2wy+DJexI9m/8PXbs2E2/EA8AABQPycnJCg0Ntf8dzwthR//78kZfX1/CDgAAfzI3G4LCAGUAAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBpbkVdgNX9+sL0oi4BLlZt8sCiLgEA4ATO7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsr0rATExOjxo0by8fHRxUrVlS3bt30888/O/Tp06ePbDabw9S0aVOHPunp6Ro0aJDKly8vb29vdenSRcePH/8jNwUAABRTRRp2NmzYoAEDBmjLli2KjY3VlStX1L59e128eNGh3/3336+TJ0/apxUrVjjMj4qK0pIlS7RgwQJ9++23SklJUadOnZSZmflHbg4AACiG3IryzVeuXOnwevbs2apYsaJ27Nihli1b2ts9PT0VFBSU6zqSkpI0c+ZMffzxx2rXrp0k6ZNPPlFoaKjWrFmjDh06FN4GAACAYq9YjdlJSkqSJAUEBDi0x8XFqWLFiqpRo4b69u2rxMRE+7wdO3bo8uXLat++vb0tJCREERERio+Pz/V90tPTlZyc7DABAABrKjZhxxijIUOG6J577lFERIS9vWPHjvr000+1bt06TZw4Udu2bVObNm2Unp4uSUpISJCHh4f8/f0d1hcYGKiEhIRc3ysmJkZ+fn72KTQ0tPA2DAAAFKkivYx1rYEDB2rXrl369ttvHdp79uxp/3dERIQaNWqksLAwff3113r44YfzXJ8xRjabLdd5I0aM0JAhQ+yvk5OTCTwAAFhUsTizM2jQIH355Zdav369KleufMO+wcHBCgsL04EDByRJQUFBysjI0Llz5xz6JSYmKjAwMNd1eHp6ytfX12ECAADWVKRhxxijgQMHavHixVq3bp3Cw8NvuszZs2d17NgxBQcHS5IaNmwod3d3xcbG2vucPHlSe/bsUfPmzQutdgAA8OdQpJexBgwYoPnz52vZsmXy8fGxj7Hx8/OTl5eXUlJSFB0dre7duys4OFiHDx/Wv/71L5UvX14PPfSQve+zzz6roUOHqly5cgoICNCLL76oOnXq2O/OAgAAf11FGnZmzJghSWrVqpVD++zZs9WnTx+VLFlSu3fv1rx583T+/HkFBwerdevWWrhwoXx8fOz9J0+eLDc3N/Xo0UOpqalq27at5syZo5IlS/6RmwMAAIohmzHGFHURRS05OVl+fn5KSkpy+fidX1+Y7tL1oehVmzywqEsAACj/f7+LxQBlAACAwkLYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAllakYScmJkaNGzeWj4+PKlasqG7duunnn3926GOMUXR0tEJCQuTl5aVWrVpp7969Dn3S09M1aNAglS9fXt7e3urSpYuOHz/+R24KAAAopoo07GzYsEEDBgzQli1bFBsbqytXrqh9+/a6ePGivc+ECRM0adIkTZ8+Xdu2bVNQUJDuu+8+Xbhwwd4nKipKS5Ys0YIFC/Ttt98qJSVFnTp1UmZmZlFsFgAAKEZsxhhT1EVkO336tCpWrKgNGzaoZcuWMsYoJCREUVFRGj58uKSrZ3ECAwP15ptvql+/fkpKSlKFChX08ccfq2fPnpKkEydOKDQ0VCtWrFCHDh1u+r7Jycny8/NTUlKSfH19XbpNv74w3aXrQ9GrNnlgUZcAAFD+/34XqzE7SUlJkqSAgABJ0qFDh5SQkKD27dvb+3h6eioyMlLx8fGSpB07dujy5csOfUJCQhQREWHvc7309HQlJyc7TAAAwJqKTdgxxmjIkCG65557FBERIUlKSEiQJAUGBjr0DQwMtM9LSEiQh4eH/P398+xzvZiYGPn5+dmn0NBQV28OAAAoJopN2Bk4cKB27dqlzz77LMc8m83m8NoYk6PtejfqM2LECCUlJdmnY8eOFbxwAABQrBWLsDNo0CB9+eWXWr9+vSpXrmxvDwoKkqQcZ2gSExPtZ3uCgoKUkZGhc+fO5dnnep6envL19XWYAACANRVp2DHGaODAgVq8eLHWrVun8PBwh/nh4eEKCgpSbGysvS0jI0MbNmxQ8+bNJUkNGzaUu7u7Q5+TJ09qz5499j4AAOCvy60o33zAgAGaP3++li1bJh8fH/sZHD8/P3l5eclmsykqKkrjx49X9erVVb16dY0fP16lS5dWr1697H2fffZZDR06VOXKlVNAQIBefPFF1alTR+3atSvKzQMAAMVAkYadGTNmSJJatWrl0D579mz16dNHkjRs2DClpqaqf//+OnfunJo0aaLVq1fLx8fH3n/y5Mlyc3NTjx49lJqaqrZt22rOnDkqWbLkH7UpAACgmCpWz9kpKjxnB87gOTsAUDz8KZ+zAwAA4GqEHQAAYGmEHQAAYGmEHQAAYGmEHQAAYGlOh502bdro/PnzOdqTk5PVpk0bV9QEAADgMk6Hnbi4OGVkZORoT0tL06ZNm1xSFAAAgKvk+6GCu3btsv973759Dt9XlZmZqZUrV6pSpUqurQ4A8JfD88msp6ifT5bvsFO/fn3ZbDbZbLZcL1d5eXlp2rRpLi0OAADgVuU77Bw6dEjGGFWtWlVbt25VhQoV7PM8PDxUsWJFvp4BAAAUO/kOO2FhYZKkrKysQisGAADA1ZweoBwTE6NZs2blaJ81a5befPNNlxQFAADgKk6Hnffff1+1atXK0V67dm299957LikKAADAVZwOOwkJCQoODs7RXqFCBZ08edIlRQEAALiK02EnNDRU//nPf3K0/+c//1FISIhLigIAAHCVfA9Qzvbcc88pKipKly9ftt+CvnbtWg0bNkxDhw51eYEAAAC3wumwM2zYMP3+++/q37+//UnKpUqV0vDhwzVixAiXFwgAAHArnA47NptNb775pkaOHKn9+/fLy8tL1atXl6enZ2HUBwAAcEucDjvZypQpo8aNG7uyFgAAAJcrUNjZtm2bPv/8cx09ejTHl4IuXrzYJYUBAAC4gtN3Yy1YsEAtWrTQvn37tGTJEl2+fFn79u3TunXr5OfnVxg1AgAAFJjTYWf8+PGaPHmyli9fLg8PD02dOlX79+9Xjx49dNtttxVGjQAAAAXmdNj57bff9OCDD0qSPD09dfHiRdlsNr3wwgv64IMPXF4gAADArXA67AQEBOjChQuSpEqVKmnPnj2SpPPnz+vSpUuurQ4AAOAWOT1A+d5771VsbKzq1KmjHj166J///KfWrVun2NhYtW3btjBqBAAAKDCnw8706dOVlpYmSRoxYoTc3d317bff6uGHH9bIkSNdXiAAAMCtcDrsBAQE2P9dokQJDRs2TMOGDXNpUQAAAK5SoOfsZGZmasmSJdq/f79sNpvuuOMOde3aVW5uBX5GIQAAQKFwOp3s2bNHXbt2VUJCgmrWrClJ+uWXX1ShQgV9+eWXqlOnjsuLBAAAKCin78Z67rnnVLt2bR0/flzff/+9vv/+ex07dkx169bV3//+98KoEQAAoMCcPrPz448/avv27fL397e3+fv7a9y4cXxXFgAAKHacPrNTs2ZNnTp1Kkd7YmKiqlWr5pKiAAAAXKVAXxcxePBgffHFFzp+/LiOHz+uL774QlFRUXrzzTeVnJxsnwAAAIqa05exOnXqJEnq0aOHbDabJMkYI0nq3Lmz/bXNZlNmZqar6gQAACgQp8PO+vXrC6MOAACAQuF02ImMjCyMOgAAAApFgZ4CmJaWpl27dikxMVFZWVkO87p06eKSwgAAAFzB6bCzcuVKPfXUUzpz5kyOeYzTAQAAxY3Td2MNHDhQjz76qE6ePKmsrCyHiaADAACKG6fDTmJiooYMGaLAwMDCqAcAAMClnA47jzzyiOLi4gqhFAAAANdzeszO9OnT9eijj2rTpk2qU6eO3N3dHeYPHjzYZcUBAADcKqfDzvz587Vq1Sp5eXkpLi7O/mBB6eoAZcIOAAAoTpwOO6+++qpee+01vfzyyypRwumrYAAAAH8op9NKRkaGevbsSdABAAB/Ck4nlt69e2vhwoWFUQsAAIDLOX0ZKzMzUxMmTNCqVatUt27dHAOUJ02a5LLiAAAAbpXTYWf37t1q0KCBJGnPnj0O864drAwAAFAc8K3nAADA0hhlDAAALC3fZ3YefvjhfPVbvHhxgYsBAABwtXyf2fHz88vX5IyNGzeqc+fOCgkJkc1m09KlSx3m9+nTRzabzWFq2rSpQ5/09HQNGjRI5cuXl7e3t7p06aLjx487VQcAALCufJ/ZmT17tsvf/OLFi6pXr56efvppde/ePdc+999/v8N7e3h4OMyPiorSV199pQULFqhcuXIaOnSoOnXqpB07dqhkyZIurxkAAPy5OD1A2ZU6duyojh073rCPp6engoKCcp2XlJSkmTNn6uOPP1a7du0kSZ988olCQ0O1Zs0adejQweU1AwCAP5diP0A5Li5OFStWVI0aNdS3b18lJiba5+3YsUOXL19W+/bt7W0hISGKiIhQfHx8nutMT09XcnKywwQAAKypWIedjh076tNPP9W6des0ceJEbdu2TW3atFF6erokKSEhQR4eHvL393dYLjAwUAkJCXmuNyYmxmGcUWhoaKFuBwAAKDpFehnrZnr27Gn/d0REhBo1aqSwsDB9/fXXN7w7zBhzwwccjhgxQkOGDLG/Tk5OJvAAAGBRxfrMzvWCg4MVFhamAwcOSJKCgoKUkZGhc+fOOfRLTExUYGBgnuvx9PSUr6+vwwQAAKypQGHn448/VosWLRQSEqIjR45IkqZMmaJly5a5tLjrnT17VseOHVNwcLAkqWHDhnJ3d1dsbKy9z8mTJ7Vnzx41b968UGsBAAB/Dk6HnRkzZmjIkCF64IEHdP78eWVmZkqSypYtqylTpji1rpSUFO3cuVM7d+6UJB06dEg7d+7U0aNHlZKSohdffFGbN2/W4cOHFRcXp86dO6t8+fJ66KGHJF199s+zzz6roUOHau3atfrhhx/0t7/9TXXq1LHfnQUAAP7anA4706ZN04cffqhXXnnF4Tk2jRo10u7du51a1/bt29WgQQP7F4sOGTJEDRo00KhRo1SyZEnt3r1bXbt2VY0aNdS7d2/VqFFDmzdvlo+Pj30dkydPVrdu3dSjRw+1aNFCpUuX1ldffcUzdgAAgKQCDFA+dOiQPZxcy9PTUxcvXnRqXa1atZIxJs/5q1atuuk6SpUqpWnTpmnatGlOvTcAAPhrcPrMTnh4uP2y07W++eYb3Xnnna6oCQAAwGWcPrPz0ksvacCAAUpLS5MxRlu3btVnn32mmJgYffTRR4VRIwAAQIE5HXaefvppXblyRcOGDdOlS5fUq1cvVapUSVOnTtVjjz1WGDUCAAAUWIEeKti3b1/17dtXZ86cUVZWlipWrOjqugAAAFzC6TE7bdq00fnz5yVJ5cuXtwed5ORktWnTxqXFAQAA3Cqnw05cXJwyMjJytKelpWnTpk0uKQoAAMBV8n0Za9euXfZ/79u3z+GLNjMzM7Vy5UpVqlTJtdUBAADconyHnfr168tms8lms+V6ucrLy4tn3QAAgGIn32Hn0KFDMsaoatWq2rp1qypUqGCf5+HhoYoVK/LUYgAAUOzkO+yEhYVJkrKysgqtGAAAAFcr0K3n0tVxO0ePHs0xWLlLly63XBQAAICrOB12Dh48qIceeki7d++WzWazf7eVzWaTJPu3oAMAABQHTt96/s9//lPh4eE6deqUSpcurb1792rjxo1q1KiR4uLiCqFEAACAgnP6zM7mzZu1bt06VahQQSVKlFCJEiV0zz33KCYmRoMHD9YPP/xQGHUCAAAUiNNndjIzM1WmTBlJV5+gfOLECUlXBzD//PPPrq0OAADgFjl9ZiciIkK7du1S1apV1aRJE02YMEEeHh764IMPVLVq1cKoEQAAoMCcDjuvvvqqLl68KEkaO3asOnXqpHvvvVflypXTwoULXV4gAADArXA67HTo0MH+76pVq2rfvn36/fff5e/vb78jCwAAoLgo8HN2rhUQEOCK1QAAALic02EnLS1N06ZN0/r165WYmJjjicrff/+9y4oDAAC4VU6HnWeeeUaxsbF65JFHdPfdd3PpCgAAFGtOh52vv/5aK1asUIsWLQqjHgAAAJdy+jk7lSpVko+PT2HUAgAA4HJOh52JEydq+PDhOnLkSGHUAwAA4FJOX8Zq1KiR0tLSVLVqVZUuXVru7u4O83///XeXFQcAAHCrnA47jz/+uP773/9q/PjxCgwMZIAyAAAo1pwOO/Hx8dq8ebPq1atXGPUAAAC4lNNjdmrVqqXU1NTCqAUAAMDlnA47b7zxhoYOHaq4uDidPXtWycnJDhMAAEBx4vRlrPvvv1+S1LZtW4d2Y4xsNpsyMzNdUxkAAIALOB121q9fXxh1AAAAFAqnw05kZGRh1AEAAFAo8hV2du3apYiICJUoUUK7du26Yd+6deu6pDAAAABXyFfYqV+/vhISElSxYkXVr19fNptNxpgc/RizAwAAipt8hZ1Dhw6pQoUK9n8DAAD8WeQr7ISFheX6bwAAgOIuX2Hnyy+/zPcKu3TpUuBiAAAAXC1fYadbt275WhljdgAAQHGTr7CTlZVV2HUAAAAUCqe/LgIAAODPJN9hZ926dbrzzjtz/f6rpKQk1a5dWxs3bnRpcQAAALcq32FnypQp6tu3r3x9fXPM8/PzU79+/TR58mSXFgcAAHCr8h12fvzxR/uXgOamffv22rFjh0uKAgAAcJV8h51Tp07J3d09z/lubm46ffq0S4oCAABwlXyHnUqVKmn37t15zt+1a5eCg4NdUhQAAICr5DvsPPDAAxo1apTS0tJyzEtNTdXo0aPVqVMnlxYHAABwq/L1nB1JevXVV7V48WLVqFFDAwcOVM2aNWWz2bR//3698847yszM1CuvvFKYtQIAADgt32EnMDBQ8fHxev755zVixAj7t57bbDZ16NBB7777rgIDAwutUAAAgILId9iRrn4J6IoVK3Tu3Dn9+uuvMsaoevXq8vf3L6z6AAAAbolTYSebv7+/Gjdu7OpaAAAAXI6viwAAAJZWpGFn48aN6ty5s0JCQmSz2bR06VKH+cYYRUdHKyQkRF5eXmrVqpX27t3r0Cc9PV2DBg1S+fLl5e3trS5duuj48eN/4FYAAIDirEjDzsWLF1WvXj1Nnz491/kTJkzQpEmTNH36dG3btk1BQUG67777dOHCBXufqKgoLVmyRAsWLNC3336rlJQUderUSZmZmX/UZgAAgGKsQGN2XKVjx47q2LFjrvOMMZoyZYpeeeUVPfzww5KkuXPnKjAwUPPnz1e/fv2UlJSkmTNn6uOPP1a7du0kSZ988olCQ0O1Zs0adejQ4Q/bFgAAUDwV2zE7hw4dUkJCgtq3b29v8/T0VGRkpOLj4yVJO3bs0OXLlx36hISEKCIiwt4nN+np6UpOTnaYAACANRXbsJOQkCBJOZ7dExgYaJ+XkJAgDw+PHLe+X9snNzExMfLz87NPoaGhLq4eAAAUF8U27GSz2WwOr40xOdqud7M+I0aMUFJSkn06duyYS2oFAADFT7ENO0FBQZKU4wxNYmKi/WxPUFCQMjIydO7cuTz75MbT01O+vr4OEwAAsKZiG3bCw8MVFBSk2NhYe1tGRoY2bNig5s2bS5IaNmwod3d3hz4nT57Unj177H0AAMBfW5HejZWSkqJff/3V/vrQoUPauXOnAgICdNtttykqKkrjx49X9erVVb16dY0fP16lS5dWr169JEl+fn569tlnNXToUJUrV04BAQF68cUXVadOHfvdWQAA4K+tSMPO9u3b1bp1a/vrIUOGSJJ69+6tOXPmaNiwYUpNTVX//v117tw5NWnSRKtXr5aPj499mcmTJ8vNzU09evRQamqq2rZtqzlz5qhkyZJ/+PYAAIDix2ayv778Lyw5OVl+fn5KSkpy+fidX1/I/YGJ+POqNnlgUZcAWBrHTesprONmfv9+F9sxOwAAAK5A2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZWrMNOdHS0bDabwxQUFGSfb4xRdHS0QkJC5OXlpVatWmnv3r1FWDEAAChuinXYkaTatWvr5MmT9mn37t32eRMmTNCkSZM0ffp0bdu2TUFBQbrvvvt04cKFIqwYAAAUJ8U+7Li5uSkoKMg+VahQQdLVszpTpkzRK6+8oocfflgRERGaO3euLl26pPnz5xdx1QAAoLgo9mHnwIEDCgkJUXh4uB577DEdPHhQknTo0CElJCSoffv29r6enp6KjIxUfHz8DdeZnp6u5ORkhwkAAFhTsQ47TZo00bx587Rq1Sp9+OGHSkhIUPPmzXX27FklJCRIkgIDAx2WCQwMtM/LS0xMjPz8/OxTaGhooW0DAAAoWsU67HTs2FHdu3dXnTp11K5dO3399deSpLlz59r72Gw2h2WMMTnarjdixAglJSXZp2PHjrm+eAAAUCwU67BzPW9vb9WpU0cHDhyw35V1/VmcxMTEHGd7rufp6SlfX1+HCQAAWNOfKuykp6dr//79Cg4OVnh4uIKCghQbG2ufn5GRoQ0bNqh58+ZFWCUAAChO3Iq6gBt58cUX1blzZ912221KTEzU2LFjlZycrN69e8tmsykqKkrjx49X9erVVb16dY0fP16lS5dWr169irp0AABQTBTrsHP8+HE9/vjjOnPmjCpUqKCmTZtqy5YtCgsLkyQNGzZMqamp6t+/v86dO6cmTZpo9erV8vHxKeLKAQBAcVGsw86CBQtuON9msyk6OlrR0dF/TEEAAOBP5081ZgcAAMBZhB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBphB0AAGBplgk77777rsLDw1WqVCk1bNhQmzZtKuqSAABAMWCJsLNw4UJFRUXplVde0Q8//KB7771XHTt21NGjR4u6NAAAUMQsEXYmTZqkZ599Vs8995zuuOMOTZkyRaGhoZoxY0ZRlwYAAIqYW1EXcKsyMjK0Y8cOvfzyyw7t7du3V3x8fK7LpKenKz093f46KSlJkpScnOzy+i6kp7p8nShahbGfAPgfjpvWU1jHzez1GmNu2O9PH3bOnDmjzMxMBQYGOrQHBgYqISEh12ViYmI0ZsyYHO2hoaGFUiMsZsawoq4AAP5cCvm4eeHCBfn5+eU5/08fdrLZbDaH18aYHG3ZRowYoSFDhthfZ2Vl6ffff1e5cuXyXAY3lpycrNDQUB07dky+vr5FXQ4AFHscN2+dMUYXLlxQSEjIDfv96cNO+fLlVbJkyRxncRITE3Oc7cnm6ekpT09Ph7ayZcsWVol/Kb6+vvzSAoATOG7emhud0cn2px+g7OHhoYYNGyo2NtahPTY2Vs2bNy+iqgAAQHHxpz+zI0lDhgzRk08+qUaNGqlZs2b64IMPdPToUf3jH/8o6tIAAEARs0TY6dmzp86ePavXXntNJ0+eVEREhFasWKGwsLCiLu0vw9PTU6NHj85xeRAAkDuOm38cm7nZ/VoAAAB/Yn/6MTsAAAA3QtgBAACWRtgBAACWRtgBAACWRtgBAOAGbDabli5dWtRl4BYQdiwqISFBgwYNUtWqVeXp6anQ0FB17txZa9euLerSbujw4cOy2WzauXNnUZcC4C/iz3q8zEt0dLTq169f1GUUK5Z4zg4cHT58WC1atFDZsmU1YcIE1a1bV5cvX9aqVas0YMAA/fTTT0VS1+XLl+Xu7l4k7w0AuSmq42VGRoY8PDycWoZj6C0wsJyOHTuaSpUqmZSUlBzzzp07Z4wx5siRI6ZLly7G29vb+Pj4mEcffdQkJCTY+40ePdrUq1fPzJs3z4SFhRlfX1/Ts2dPk5ycbO+TmZlp3njjDXP77bcbDw8PExoaasaOHWuMMebQoUNGklm4cKGJjIw0np6eZtasWcYYY2bNmmVq1aplPD09Tc2aNc0777xjX6ckhykyMtI+70bLAUBB5Od4Kcl8+OGHplu3bsbLy8tUq1bNLFu2zN7vypUr5plnnjFVqlQxpUqVMjVq1DBTpkxxWFfv3r1N165dzfjx401wcLAJCwszxhhz7Ngx07NnT+Pv729Kly5tGjZsaLZs2WKM+d9xeObMmSY8PNzYbDaTlZVlzp8/b/r27WsqVKhgfHx8TOvWrc3OnTuNMcbMnj07x3F09uzZxhhzw+WsjrBjMWfPnjU2m82MHz8+zz5ZWVmmQYMG5p577jHbt283W7ZsMXfddZdDsBg9erQpU6aMefjhh83u3bvNxo0bTVBQkPnXv/5l7zNs2DDj7+9v5syZY3799VezadMm8+GHHxpj/hd2qlSpYhYtWmQOHjxo/vvf/5oPPvjABAcH29sWLVpkAgICzJw5c4wxxmzdutVIMmvWrDEnT540Z8+eNcaYmy4HAM7Kz/HSmKthp3Llymb+/PnmwIEDZvDgwaZMmTL241NGRoYZNWqU2bp1qzl48KD55JNPTOnSpc3ChQvt6+jdu7cpU6aMefLJJ82ePXvM7t27zYULF0zVqlXNvffeazZt2mQOHDhgFi5caOLj440xV4/D3t7epkOHDub77783P/74o8nKyjItWrQwnTt3Ntu2bTO//PKLGTp0qClXrpw5e/asuXTpkhk6dKipXbu2OXnypDl58qS5dOnSTZezOsKOxXz33XdGklm8eHGefVavXm1Klixpjh49am/bu3evkWS2bt1qjLn6S1a6dGmHMzkvvfSSadKkiTHGmOTkZOPp6WkPN9fLDjvX/+8mNDTUzJ8/36Ht9ddfN82aNXNY7ocffnBqOQBwVn6Ol8ZcDTuvvvqq/XVKSoqx2Wzmm2++yXOZ/v37m+7du9tf9+7d2wQGBpr09HR72/vvv298fHzyDBujR4827u7uJjEx0d62du1a4+vra9LS0hz63n777eb999+3L1evXj2H+flZzsoYs2Mx5v++/cNms+XZZ//+/QoNDVVoaKi97c4771TZsmW1f/9+NW7cWJJUpUoV+fj42PsEBwcrMTHRvo709HS1bdv2hvU0atTI/u/Tp0/r2LFjevbZZ9W3b197+5UrV+Tn55fnOgq6HADcSH6Ol9nq1q1r/7e3t7d8fHzsx0NJeu+99/TRRx/pyJEjSk1NVUZGRo5BwnXq1HEYp7Nz5041aNBAAQEBeb5vWFiYKlSoYH+9Y8cOpaSkqFy5cg79UlNT9dtvv+W5noIuZxWEHYupXr26bDab9u/fr27duuXaxxiT6y/39e3XD4Sz2WzKysqSJHl5eeWrHm9vb/u/s5f98MMP1aRJE4d+JUuWzHMdBV0OAG4kP8fLbDc6Hv6///f/9MILL2jixIlq1qyZfHx89NZbb+m7775zWOba46GUv+Po9ctkZWUpODhYcXFxOfqWLVs2z/UUdDmrIOxYTEBAgDp06KB33nlHgwcPzvGLcv78ed155506evSojh07Zj+7s2/fPiUlJemOO+7I1/tUr15dXl5eWrt2rZ577rl8LRMYGKhKlSrp4MGDeuKJJ3Ltk/2/nszMTKeWAwBn5ed4mZ8gsGnTJjVv3lz9+/e3t+XnbEndunX10Ucf6ffff7/h2Z1r3XXXXUpISJCbm5uqVKmSax8PDw+HY2h+l7MynrNjQe+++64yMzN19913a9GiRTpw4ID279+vf//732rWrJnatWununXr6oknntD333+vrVu36qmnnlJkZKTDZacbKVWqlIYPH65hw4Zp3rx5+u2337RlyxbNnDnzhstFR0crJiZGU6dO1S+//KLdu3dr9uzZmjRpkiSpYsWK8vLy0sqVK3Xq1CklJSXlazkAKIibHS/zo1q1atq+fbtWrVqlX375RSNHjtS2bdtuutzjjz+uoKAgdevWTf/5z3908OBBLVq0SJs3b85zmXbt2qlZs2bq1q2bVq1apcOHDys+Pl6vvvqqtm/fLunqEIRDhw5p586dOnPmjNLT0/O1nKUV7ZAhFJYTJ06YAQMGmLCwMOPh4WEqVapkunTpYtavX2+Myf+t59eaPHmy/XZJY67eej527FgTFhZm3N3dzW233Wa/qyGvgcbGGPPpp5+a+vXrGw8PD+Pv729atmzpMEDwww8/NKGhoaZEiRIOd4jdbDkAKIibHS8lmSVLljgs4+fnZ7+lOy0tzfTp08f4+fmZsmXLmueff968/PLLDsfQ7FvPr3f48GHTvXt34+vra0qXLm0aNWpkvvvuO2NM7sdhY67eIDJo0CATEhJi3N3dTWhoqHniiSfsN52kpaWZ7t27m7Jlyzrcen6z5azMZsz/jdACAACwIC5jAQAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAASyPsAAAAS/v/yZLEXajgbvQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "labels = ['Concrete', 'Charcrete']\n",
    "width = 0.35       # the width of the bars: can also be len(x) sequence\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.bar(labels, RoW_results, width, color='palevioletred')\n",
    "\n",
    "ax.set_ylabel('Climate Impact')\n",
    "ax.set_title('Manipulated market (RoW -> Charcrete)')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "93773418-2043-4dff-a853-a7a7def8ba74",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.savefig(wdir+\"/figures/market_concrete_charcrete\") "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb3d7845-9af7-4ba4-a8e9-00fc1e519580",
   "metadata": {},
   "source": [
    "#### Replacing cement in charcrete"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "869b85ca-252e-4dc5-8f32-b8a6d309ab1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "actual_cement = bd.get_activity(database=\"ei 3.8 cutoff\", name=\"market for cement, Portland\", location = \"Europe without Switzerland\")\n",
    "best_cement = bd.get_activity(database=\"ei 3.8 cutoff\", name=\"market for cement, blast furnace slag 81-95%\", location = \"CH\")\n",
    "worst_cement = bd.get_activity(database=\"ei 3.8 cutoff\", name=\"market for cement, Portland\", location = \"ZA\")\n",
    "\n",
    "charcrete_best = replace_stuff(charcrete, actual_cement, best_cement, ipcc)\n",
    "charcrete_worst = replace_stuff(charcrete, actual_cement, worst_cement, ipcc)\n",
    "\n",
    "results = [*charcrete_best, charcrete_worst[1]]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90e87bf7-6c8e-4a19-95ba-7f2181448a2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "charcrete_best\n",
    "charcrete_worst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56a16d70-be43-487e-9521-c9aa1e60d76b",
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = ['Base case (market for cement, Portland, EU w/ CH)', 'Best (market for cement, blast furnace slag 81-95%, CH)', 'Worst (market for cement, Portland, ZA)' ]\n",
    "width = 0.35       # the width of the bars: can also be len(x) sequence\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.bar(labels, results, width, color='palevioletred')\n",
    "\n",
    "ax.set_ylabel('Climate Impact')\n",
    "ax.set_title('Change of cement')\n",
    "\n",
    "plt.xticks(rotation=45)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfa37955-dd04-4793-8b09-afaf0c45e6ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.savefig(wdir+\"/figures/cement_in_charcrete\") "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d6df822-0540-46ff-b3e8-2076a5dccdeb",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Scenarios of charcrete"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b6a49a2-64fe-4e88-a580-bff48e14363d",
   "metadata": {},
   "outputs": [],
   "source": [
    "scenario_df = pd.read_excel(wdir+\"/data/Charcrete_database_38.xlsx\", sheet_name=\"scenario\")\n",
    "scenario_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65822e04-2aa2-4571-9963-f543ba9fc5e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_id(x):\n",
    "    if x[\"type\"] == \"technosphere\" or x[\"type\"] == \"production\" :\n",
    "        try:\n",
    "            activity = bd.get_activity(database=x[\"database\"], name=x[\"Activity name\"], location=x[\"location\"])\n",
    "            id_val = activity.id\n",
    "        except:\n",
    "            print(\"erreur : \"+x[\"Activity name\"])\n",
    "    if x[\"type\"] == \"biosphere\":\n",
    "            act = bd.get_activity(database=\"Charcrete_38\", name=\"Charcrete production 5%\", location=\"SE\")\n",
    "            activity = [flow for flow in act.biosphere() if x[\"Activity name\"] in flow[\"name\"]][0]\n",
    "            id_val = activity.input.id\n",
    "    return id_val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9284542b-a2c7-47f1-bf9f-7f00fbe6c027",
   "metadata": {},
   "outputs": [],
   "source": [
    "scenario_df[\"id\"] = scenario_df.apply(find_id, axis=1)\n",
    "scenario_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bce16b6-26d4-4f78-acc9-8ac591b3558b",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = scenario_df[(scenario_df[\"type\"]==\"technosphere\") | (scenario_df[\"type\"]==\"production\")][[\"Charcrete production 5%\", \"Charcrete production 7.5%\", \"Charcrete production 10%\"]].values.tolist()\n",
    "t_data_array = np.array(data)\n",
    "t_data_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29a08fcc-530c-4b51-b245-68900d5f07bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "[flow for flow in charcrete.biosphere()][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6493e4a6-ba69-4066-96fd-8460ba020191",
   "metadata": {},
   "outputs": [],
   "source": [
    "ba.print_recursive_supply_chain(charcrete)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98b226c9-9b57-44a7-94b6-077ba591c71e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ind_list = scenario_df[(scenario_df[\"type\"]==\"technosphere\") | (scenario_df[\"type\"]==\"production\")][\"id\"].values.tolist()\n",
    "indices = [(ind, 33105) for ind in ind_list]\n",
    "t_indices = np.array(indices, dtype=bwp.INDICES_DTYPE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91d2dea2-781d-4c70-8626-efe62c3b2f3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_flip = np.array([False]+[True for i in range(len(scenario_df[(scenario_df[\"type\"]==\"technosphere\")]))])\n",
    "t_flip"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d9c76c6-70e3-45f1-b923-d37a4371cd87",
   "metadata": {},
   "outputs": [],
   "source": [
    "ind_list = scenario_df[(scenario_df[\"type\"]==\"biosphere\")][\"id\"].values.tolist()\n",
    "indices = [(ind, 33105) for ind in ind_list]\n",
    "b_indices = np.array(indices, dtype=bwp.INDICES_DTYPE)\n",
    "b_indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30ec88f3-e81d-440e-9806-7a4bd791290a",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = scenario_df[(scenario_df[\"type\"]==\"biosphere\")][[\"Charcrete production 5%\", \"Charcrete production 7.5%\", \"Charcrete production 10%\"]].values.tolist()\n",
    "b_data = np.array(data)\n",
    "b_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23856ba5-9160-44dc-b9e6-f0a605e373e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_arrays = bwp.create_datapackage(sequential=True, seed=25323) # sequential sampling of columns = [0,1,2,0,1,...]\n",
    "dp_arrays.add_persistent_array(\n",
    "    matrix='technosphere_matrix',\n",
    "    indices_array=t_indices,\n",
    "    data_array=t_data_array,\n",
    "    flip_array=t_flip,\n",
    ")\n",
    "dp_arrays.add_persistent_array(\n",
    "    matrix='biosphere_matrix',\n",
    "    indices_array=b_indices,\n",
    "    data_array=b_data,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e33fe9f-bb99-4b57-86b1-135e47928512",
   "metadata": {},
   "outputs": [],
   "source": [
    "charcrete"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0643a3aa-44d8-499c-ad30-9f8e8ff3e11d",
   "metadata": {},
   "outputs": [],
   "source": [
    "fu, data_objs, _ = bd.prepare_lca_inputs({charcrete: 1}, method=ipcc)\n",
    "lca_b = bc.LCA(\n",
    "    fu,\n",
    "    data_objs=[*data_objs, dp_arrays],\n",
    "    use_distributions=False,\n",
    "    use_arrays=True, \n",
    "#     seed_override=42,  # Seed should not be used it randoms the columns order\n",
    ")\n",
    "lca_b.lci()\n",
    "lca_b.lcia()\n",
    "lca_b.score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5c15875-5175-4b14-b58b-ca3d56b8d474",
   "metadata": {},
   "outputs": [],
   "source": [
    "lca_b.keep_first_iteration()\n",
    "iterations = 3\n",
    "scores_b = [lca_b.score for _ in zip(range(iterations), lca_b)]\n",
    "scores_b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "215d48e5-b6d7-4439-b2d7-8ceb40d6fea8",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Correlated uncertainties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bdadbb3-a32f-4766-8d4d-57470b5dd3b0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "1de92aa4-d867-4149-b98f-86d6223738f6",
   "metadata": {},
   "source": [
    "## Sensitivity analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80afe3da-9600-45f0-8d11-768e1d595744",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "me = bd.Method(ipcc).datapackage()\n",
    "iterations_mc = 10\n",
    "seed_mc = 5555  # specifying random seeds is needed for reproducibility of results\n",
    "\n",
    "# Run MC simulations for all cases\n",
    "case = \"foreground\"\n",
    "fu, data_objs, _ = bd.prepare_lca_inputs({charcrete:1}, method=ipcc)\n",
    "\n",
    "kr_nounct = data_objs[0].exclude({\"kind\": \"distributions\"})\n",
    "ei_nounct = data_objs[2].exclude({\"kind\": \"distributions\"})\n",
    "char_dp_nounct = data_objs[3].exclude({\"kind\": \"distributions\"})\n",
    "# Case 1: uncertainties are present in both background and foreground\n",
    "dps_fb = [data_objs[0], data_objs[1], data_objs[2], data_objs[3], me]\n",
    "\n",
    "# Case 2: uncertainties are present ONLY in FOREground\n",
    "dps_fg = [kr_nounct, data_objs[1], ei_nounct, data_objs[3], me]\n",
    "\n",
    "# Case 3: uncertainties are present ONLY in BACKground\n",
    "dps_bg = [data_objs[0], data_objs[1], data_objs[2], char_dp_nounct, me]\n",
    "\n",
    "cases = {\n",
    "    \"foreground_background\": {\"datapackages\": dps_fb}, \n",
    "    \"foreground\": {\"datapackages\": dps_fg}, \n",
    "    \"background\": {\"datapackages\": dps_bg}\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31606180-f82f-474b-be7a-4b6f0b8a555a",
   "metadata": {},
   "outputs": [],
   "source": [
    "lca_mc = bc.LCA(\n",
    "    fu, \n",
    "    data_objs=cases[\"background\"][\"datapackages\"],\n",
    "    use_distributions=True, \n",
    "    seed_override=seed_mc,\n",
    ")\n",
    "lca_mc.lci()\n",
    "lca_mc.lcia()\n",
    "\n",
    "# Run Monte Carlo simulations\n",
    "lca_scores_mc = []\n",
    "for i in range(iterations_mc):\n",
    "        next(lca_mc)\n",
    "        lca_scores_mc.append(lca_mc.score)\n",
    "lca_scores_mc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14001119-b547-45bc-a5f8-a4695ad33595",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "iterations_mc = 100\n",
    "seed_mc = 5555  # specifying random seeds is needed for reproducibility of results\n",
    "\n",
    "# Run MC simulations for all cases\n",
    "for case, dict_ in cases.items():\n",
    "    print(case)\n",
    "    dps = dict_[\"datapackages\"]\n",
    "    lca_mc = bc.LCA(\n",
    "        {charcrete.id: 1}, \n",
    "        data_objs=dps,\n",
    "        use_distributions=True, \n",
    "        seed_override=seed_mc,\n",
    "    )\n",
    "    lca_mc.lci()\n",
    "    lca_mc.lcia()\n",
    "\n",
    "    # Since MC takes time, make sure to save the results!\n",
    "    # Also, read lca_scores from file if you already computed them once\n",
    "    output_dir = Path(wdir+\"/outputs\")\n",
    "    fp_mc = output_dir / f'lca_scores_{seed_mc}_{iterations_mc}_{case}.json'\n",
    "    replace = True\n",
    "    if replace == False:\n",
    "        # Read LCIA scores\n",
    "        with open(fp_mc, 'r') as f:\n",
    "            lca_scores_mc = json.load(f)\n",
    "    else:\n",
    "        # Run Monte Carlo simulations\n",
    "        lca_scores_mc = []\n",
    "        for i in range(iterations_mc):\n",
    "            next(lca_mc)\n",
    "            lca_scores_mc.append(lca_mc.score)\n",
    "        print(\"Save LCIA scores\")\n",
    "        with open(fp_mc, 'w') as f:\n",
    "            json.dump(lca_scores_mc, f)\n",
    "\n",
    "    dict_['lca_scores_mc'] = lca_scores_mc\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28ead0f5-2d80-42b5-a442-507c3eaea753",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_bins = 100\n",
    "lca_scores_mc_all = np.hstack([dict_['lca_scores_mc'] for dict_ in cases.values()])\n",
    "lca_scores_mc_all = lca_scores_mc_all[\n",
    "    np.logical_and(\n",
    "        lca_scores_mc_all > np.percentile(lca_scores_mc_all, 0),\n",
    "        lca_scores_mc_all < np.percentile(lca_scores_mc_all, 99),\n",
    "    )\n",
    "]\n",
    "bins = np.linspace(min(lca_scores_mc_all), max(lca_scores_mc_all), num_bins, endpoint=True)\n",
    "midbins = (bins[:-1] + bins[1:]) / 2\n",
    "width = (bins[1]-bins[0])*0.8\n",
    "\n",
    "plt.figure(figsize=(10, 6), dpi=80)\n",
    "\n",
    "# Uncertainty distributions for all cases\n",
    "for case, dict_ in cases.items():\n",
    "    lca_scores_mc = dict_['lca_scores_mc']\n",
    "    freq, _ = np.histogram(lca_scores_mc, bins=bins, density=False)\n",
    "    plt.bar(midbins, freq, width=width, label=case, alpha=0.6)\n",
    "\n",
    "plt.xlabel(\"LCIA scores, kg CO2-eq\")\n",
    "plt.ylabel(\"Count\")\n",
    "plt.legend(loc=\"upper right\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d1fe527-8497-4151-aee3-c87a3e181234",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c61aa94b-97c5-4cae-aeb2-f6a7be825683",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96985b8e-d26b-4956-a5cf-dc648b9415af",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fe251c1-716b-4a05-a0e1-fff437d37d85",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:bw25]",
   "language": "python",
   "name": "conda-env-bw25-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}