2 - Building and using matrices in bw2calc.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1c6c844b-622e-4e74-83b1-c8a1071aaaea",
   "metadata": {},
   "source": [
    "# 2 - Building and using matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7375df6-8326-4f36-a765-535e8f9b0576",
   "metadata": {},
   "outputs": [],
   "source": [
    "import bw_processing as bwp\n",
    "import matrix_utils as mu\n",
    "import bw2calc as bc\n",
    "import numpy as np\n",
    "import seaborn as sb\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a8c9830-d0ca-479a-9f6a-020b5dd8bc15",
   "metadata": {},
   "source": [
    "Before we dive into it, let's think about what we need to actually build a matrix. What specific data would you need? What don't you need?\n",
    "\n",
    "## Exercise\n",
    "\n",
    "Please think about the minimal set of information you would need to build a *sparse matrix* using [scipy.sparse.coo_matrix](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.coo_matrix.html) (sparse matrices store only non-zero values). Then, create this information as Numpy arrays and actually build a sparse matrix.\n",
    "\n",
    "Here is the matrix you should build:\n",
    "\n",
    "$$\\begin{bmatrix} 0 & 1 \\\\ 2 & 3 \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1635ef16-d33c-4663-9237-747b45a2ccde",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Hint"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b41ae23-bc12-4ddc-9292-49da20f33c32",
   "metadata": {
    "tags": []
   },
   "source": [
    "You will need three Numpy arrays: one for the data, one for the row indices, and one for the column indices."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41a9de30-1ac1-40af-9ba8-092564e31278",
   "metadata": {},
   "source": [
    "## Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47d29ce0-79a9-4afd-b1db-7c76faa40bf6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from scipy import sparse\n",
    "\n",
    "data = np.array([1, 2, 3])\n",
    "rows = np.array([0, 1, 1])\n",
    "cols = np.array([1, 0, 1])\n",
    "\n",
    "matrix = sparse.coo_matrix((data, (rows, cols)), (2, 2))\n",
    "matrix.toarray()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20404d5e-261f-4922-93e0-bd35ad0d0a15",
   "metadata": {},
   "source": [
    "## `bw_processing`\n",
    "\n",
    "We can run into difficulties when we want to store this data. The library `bw_processing` helps us create data packages, which can store this matrix-building data on variety of file systems. You can read the [`bw_processing` README](github.com/brightway-lca/bw_processing) for more information, and can see the [PyFilesystem2 Docs](https://docs.pyfilesystem.org/en/latest/) for more on the filesystems that can be used.\n",
    "\n",
    "Let's define this same matrix in `bw_processing`.\n",
    "\n",
    "Matrices by definition are two-dimensional, so we know that to build matrices we will always need to specify the row and column indices of the data. We combine these two arrays into a single Numpy [structured array](https://numpy.org/doc/stable/user/basics.rec.html), which uses the labels `row` and `col`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04ba1424-c124-4391-805e-68588143ef72",
   "metadata": {},
   "outputs": [],
   "source": [
    "import bw_processing as bwp\n",
    "import numpy as np\n",
    "\n",
    "indices_array = np.array([(0, 1), (1, 0), (1, 1)], dtype=bwp.INDICES_DTYPE)\n",
    "indices_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6396041c-17bf-4a0b-b72a-7f7d73a180f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "indices_array['row']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3dcba032-da1a-4818-8126-f4ab1e418c4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "bwp.INDICES_DTYPE"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1afb5959-335f-4ba0-a40b-91324cc408e0",
   "metadata": {},
   "source": [
    "The data array is the same as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b43f5478-c4aa-4b00-96c2-7d06f115f9b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_array = np.array([1, 2, 3])\n",
    "data_array"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d887832-21d2-4ca6-aba2-d4f9e27c2c07",
   "metadata": {},
   "source": [
    "This is all we need to create a data package:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9174a0c1-555f-4857-949f-6276d7d39bf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "dp = bwp.create_datapackage()\n",
    "\n",
    "dp.add_persistent_vector(\n",
    "    matrix=\"some name\",\n",
    "    data_array=data_array,\n",
    "    name=\"some name\",\n",
    "    indices_array=indices_array,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8081b317-4d2e-447e-a5dc-2e2485920030",
   "metadata": {},
   "source": [
    "But before this gets too abstract, let's do the same for our example system:\n",
    "\n",
    "<img src='images/simple-graph.png' width='400'>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c52f1897-f6c7-410c-9cfa-8fc0b0be51d6",
   "metadata": {},
   "source": [
    "We need integer indices for each concept. The numbers themselves are arbitrary, the point is that we will use them consistently in the datapackages, and maintain a mapping from our concepts to these numbers. This might be similar to other fundamental concepts such as a [basis](https://en.wikipedia.org/wiki/Basis_(linear_algebra))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a30f9eb-8914-4273-925f-3d78893cd7bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "natural_gas = 101\n",
    "carbon_fibre = 102\n",
    "bike = 103\n",
    "co2 = 201"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92f8c609-333e-4d9c-9829-1b5d6b9ed146",
   "metadata": {},
   "source": [
    "In our technosphere matrix, we will have three production exchanges (each of amount 1), and two consumption exchanges. Our biosphere matrix will only have one number, the emission of $CO_{2}$ from carbon fibre production.\n",
    "\n",
    "Our matrices should look like this:\n",
    "\n",
    "## Technosphere matrix\n",
    "\n",
    "Dimensions are products (rows) by activities (columns).\n",
    "\n",
    "$$\n",
    "\\left[\\begin{array}{ccc} \n",
    "1 & 0 & 0\\\\\n",
    "-237 & 1 & 0\\\\\n",
    "0 & -2.5 & 1\\\\\n",
    "\\end{array}\\right]\n",
    "$$\n",
    "\n",
    "## Biosphere matrix\n",
    "\n",
    "Dimensions are flows (rows) by activities (columns).\n",
    "\n",
    "$$\n",
    "\\left[\\begin{array}{ccc} \n",
    "0 & 26.6 & 0\\\\ \n",
    "\\end{array}\\right]\n",
    "$$ \n",
    "\n",
    "## Characterization matrix\n",
    "\n",
    "Dimensions are flows (rows) by flows (columns).\n",
    "\n",
    "$$\n",
    "\\left[\\begin{array}{c} \n",
    "1\\\\ \n",
    "\\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a2d45da-125b-40a8-a0de-a10f05af12d1",
   "metadata": {},
   "source": [
    "# A first LCA calculation\n",
    "\n",
    "To use datapackages, we need to know four thing for each point: the row and column ids, the data value, and the sign. We only need to store the non-zero points. Here is an example for the technosphere matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23951444-84a2-40b5-9687-e3cbd3053dc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_data = np.array([\n",
    "    1,   # production of natural gas\n",
    "    1,   # production of carbon fibre\n",
    "    1,   # production of bike\n",
    "    237, # input of natural gas\n",
    "    2.5, # input of carbon fibre\n",
    "])\n",
    "t_indices = np.array([\n",
    "    (natural_gas, natural_gas), # production of natural gas\n",
    "    (carbon_fibre, carbon_fibre), # production of carbon fibre\n",
    "    (bike, bike), # production of bike\n",
    "    (natural_gas, carbon_fibre), # input of natural gas\n",
    "    (carbon_fibre, bike), # input of carbon fibre\n",
    "    ], \n",
    "    dtype=bwp.INDICES_DTYPE\n",
    ")\n",
    "t_flip = np.array([False, False, False, True, True]) # Numerical sign of the inputs needs to be flipped negative"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b29dcb59-79eb-4ce4-8881-024b7b0892ab",
   "metadata": {},
   "source": [
    "And similarly for the other matrices (no need to flip signs, so we skip that part):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15eaabd8-0b91-4973-a7d6-cb6b5ed9f54e",
   "metadata": {},
   "outputs": [],
   "source": [
    "b_data = np.array([26.6])\n",
    "b_indices = np.array([\n",
    "    (co2, carbon_fibre), # emission of CO2\n",
    "    ], \n",
    "    dtype=bwp.INDICES_DTYPE\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e6efb7c-cac3-4536-b93e-da76fa2a49a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "c_data = np.array([1])\n",
    "c_indices = np.array([\n",
    "    (co2, co2), # CF of CO2\n",
    "    ], \n",
    "    dtype=bwp.INDICES_DTYPE\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53f00163-c997-4555-a563-8a767e75f2fa",
   "metadata": {},
   "source": [
    "We can now create our datapackage, and add the data for all three matrices to it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7fc1602-4fad-45c8-a84f-f47943bfe739",
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_static = bwp.create_datapackage()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46e9a1b2-6ec3-4627-924a-646a4e8e64a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_static.add_persistent_vector(\n",
    "    matrix='technosphere_matrix',\n",
    "    indices_array=t_indices,\n",
    "    data_array=t_data,\n",
    "    flip_array=t_flip,\n",
    ")\n",
    "dp_static.add_persistent_vector(\n",
    "    matrix='biosphere_matrix',\n",
    "    indices_array=b_indices,\n",
    "    data_array=b_data,\n",
    ")\n",
    "dp_static.add_persistent_vector(\n",
    "    matrix='characterization_matrix',\n",
    "    indices_array=c_indices,\n",
    "    data_array=c_data,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58cc32ba-d495-4a27-a7e6-a9c9ceaa7b85",
   "metadata": {},
   "source": [
    "This is already enough to calculate an LCA score:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4de9b4ab-ce03-4ab3-9e56-42420875356c",
   "metadata": {},
   "outputs": [],
   "source": [
    "lca = bc.LCA(\n",
    "    demand={bike: 1},\n",
    "    data_objs=[dp_static],\n",
    ")\n",
    "lca.lci()\n",
    "lca.lcia()\n",
    "lca.score"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c624a683-9374-47f3-b496-7af3b99be0ff",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Do the same for your steel bike. Make sure you get the same LCIA results as before."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f964f4d-28ca-4bbc-8ed2-18068c37cbe5",
   "metadata": {},
   "source": [
    "# A stochastic LCA\n",
    "\n",
    "This is pretty boring so far. The only marginally interesting thing is that we could put our datapackage on another computer, and do calculations on that computer - i.e. we are ready for cloud calculations. Indeed, the library `matrix_utils` allows us to specify datapackages on FTP sites and cloud storage such as Amazon S3.\n",
    "\n",
    "But we can make our example more interesting by adding uncertainty. To do this, we will use the standards in the [stats_arrays](https://stats-arrays.readthedocs.io/en/latest/) library.\n",
    "\n",
    "We will only add uncertainty to the inputs, and the CO2 emission. The other edges will have uncertainty type zero - no uncertainty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c82ee36-07d3-4a3e-aaa3-164682d68344",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_uncertainty = np.array([\n",
    "        (0, 1, np.NaN, np.NaN, np.NaN, np.NaN, False),\n",
    "        (0, 1, np.NaN, np.NaN, np.NaN, np.NaN, False),    \n",
    "        (0, 1, np.NaN, np.NaN, np.NaN, np.NaN, False),    \n",
    "        (5, 237, np.NaN, np.NaN, 200, 300, False), # triangular uncertainty from 200 to 300  \n",
    "        (5, 2.5, np.NaN, np.NaN, 2, 3, False), # triangular uncertainty from 2 to 3\n",
    "    ], \n",
    "    dtype=bwp.UNCERTAINTY_DTYPE\n",
    ")\n",
    "b_uncertainty = np.array([\n",
    "        (3, 26.6, 1.5, np.NaN, np.NaN, np.NaN, False), # normal uncertainty with std. dev. of 1.5\n",
    "    ], \n",
    "    dtype=bwp.UNCERTAINTY_DTYPE\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efc258ba-5ff1-4f87-bb7d-7b4d99e51e73",
   "metadata": {},
   "source": [
    "We then will create another data package, and use it to do Monte Carlo:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bde2e4dd-cf2e-4ad7-8f22-b23c7f7cc847",
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_stochastic = bwp.create_datapackage()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1613950-5464-4774-b09c-e06dabe6c57e",
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_stochastic.add_persistent_vector(\n",
    "    matrix='technosphere_matrix',\n",
    "    indices_array=t_indices,\n",
    "    data_array=t_data,\n",
    "    flip_array=t_flip,\n",
    "    distributions_array=t_uncertainty,\n",
    ")\n",
    "dp_stochastic.add_persistent_vector(\n",
    "    matrix='biosphere_matrix',\n",
    "    indices_array=b_indices,\n",
    "    data_array=b_data,\n",
    "    distributions_array=b_uncertainty,\n",
    ")\n",
    "dp_stochastic.add_persistent_vector(\n",
    "    matrix='characterization_matrix',\n",
    "    indices_array=c_indices,\n",
    "    data_array=c_data,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "984c71e8-dfcc-445e-9acd-9ba0bc8c9fc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "lca = bc.LCA(\n",
    "    demand={bike: 1},\n",
    "    data_objs=[dp_stochastic],\n",
    "    use_distributions=True,\n",
    ")\n",
    "lca.lci()\n",
    "lca.lcia()\n",
    "    \n",
    "mc_results = [lca.score for _ in zip(range(50), lca)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "666b3486-280f-4503-85e4-9e25ea5f3f4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "sb.displot(np.array(mc_results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bce3c9f5-cfd6-40f0-94af-5e36b8bcc3a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0e37bd1-bf63-4037-9916-06b9db6f9895",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "The numpy function for correlation coefficients is `corrcoef`. Can we do a basic sensitivity analysis on the Monte Carlo parameters? See if you can capture the matrix values during each iteration, and then calculate the correlation coefficients between this matrix values and the LCIA scores."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "688ed0a3-f5ff-476a-b1d6-cc1bc63abc6d",
   "metadata": {},
   "source": [
    "# Correlating exchanges\n",
    "\n",
    "If we sample each exchange independently, we lose some of the logic of our product system. In our example system, the $CO_{2}$ should be tied directly to the natural gas consumption, but it isn't:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b926b98-7cee-4c69-a846-35f6cbb97aa3",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = []\n",
    "\n",
    "for _ in range(100):\n",
    "    next(lca)\n",
    "    params.append({\n",
    "        'co2': lca.biosphere_matrix[lca.dicts.biosphere[co2], lca.dicts.activity[carbon_fibre]],\n",
    "        'ng': -1 * lca.technosphere_matrix[lca.dicts.product[natural_gas], lca.dicts.activity[carbon_fibre]],\n",
    "    })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c60872b-2309-4e15-b551-741266146c30",
   "metadata": {},
   "outputs": [],
   "source": [
    "sb.scatterplot(data=pd.DataFrame(params), x='ng', y='co2')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7347e5e4-776f-4a1b-8ac7-dab734c74801",
   "metadata": {},
   "source": [
    "We fix this by drawing correlated samples ahead of time - this is the idea behind the `presamples` package, whose functionality is now integrated directly into Brightway. We will make up a model to illustrate how this would work; we can assume that (again, totally made up) 80% of the natural gas goes to energy, the ratio stays the same, but with a little bit of noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6412f750-65bd-42f7-b290-0cef96ddd4be",
   "metadata": {},
   "outputs": [],
   "source": [
    "ng_samples = np.random.triangular(200, 237, 300, size=100)\n",
    "co2_samples = 26.6 / 237 * ng_samples * np.random.normal(loc=1, scale=0.025, size=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cf9ceb1-b7f3-4b8f-b13f-3d4c05f46553",
   "metadata": {},
   "source": [
    "We can now add another datapackage which will overwrite our previous values (but only where we tell it to). Note that we are now creating `arrays`, not `vectors`.\n",
    "\n",
    "We need to tell the datapackage that these two resources are correlated - otherwise there would be two RNGs used to samples them independently again. We can do this by either setting `sequential=True` (in which case they would start with column 0, then column 1, etc.), or by using the same RNG seed for both resources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18c956d5-a89c-423f-aad2-906012590a52",
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_correlated = bwp.create_datapackage(seed=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b13c9c8a-cecc-4909-80c9-d8d5b451861f",
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_correlated.add_persistent_array(\n",
    "    matrix='technosphere_matrix',\n",
    "    indices_array=np.array([(natural_gas, carbon_fibre)], dtype=bwp.INDICES_DTYPE),\n",
    "    data_array=ng_samples.reshape((1, -1)),\n",
    "    flip_array=np.array([True]) ,\n",
    ")\n",
    "dp_correlated.add_persistent_array(\n",
    "    matrix='biosphere_matrix',\n",
    "    indices_array=np.array([(co2, carbon_fibre)], dtype=bwp.INDICES_DTYPE),\n",
    "    data_array=co2_samples.reshape((1, -1)),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "192177b8-45e6-4f64-a201-c0fb0555f630",
   "metadata": {},
   "source": [
    "And now the use of our correlated samples reflects reality a lot better. Note that we justs add the new datapackage after the original one, and set `use_arrays=True`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ac9c5fb-a726-40c4-9e51-d1e422fe93a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "lca = bc.LCA(\n",
    "    demand={bike: 1},\n",
    "    data_objs=[dp_stochastic, dp_correlated],\n",
    "    use_distributions=True,\n",
    "    use_arrays=True,\n",
    ")\n",
    "lca.lci()\n",
    "lca.lcia()\n",
    "\n",
    "corr_params = []\n",
    "\n",
    "for _ in range(100):\n",
    "    next(lca)\n",
    "    corr_params.append({\n",
    "        'co2': lca.biosphere_matrix[lca.dicts.biosphere[co2], lca.dicts.activity[carbon_fibre]],\n",
    "        'ng': -1 * lca.technosphere_matrix[lca.dicts.product[natural_gas], lca.dicts.activity[carbon_fibre]],\n",
    "    })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d71fa519-7b1b-4122-be55-11cae2164437",
   "metadata": {},
   "outputs": [],
   "source": [
    "sb.scatterplot(data=pd.DataFrame(corr_params), x='ng', y='co2')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da36cc00-a224-4334-8da1-a55a24f05cf3",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Do the same for your steel bike - correlate iron ore resource consumption to mining."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8216209b-1a34-4e4c-86e6-a69862e463fe",
   "metadata": {},
   "source": [
    "# Using arrays for scenarios\n",
    "\n",
    "In addition to correlated samples, we can use arrays for scenarios. Let's imagine two different possibilities: a lightweight bike with 1.5 kilograms of carbon fibre, and a technology platform that allows for efficient bike sharing making each bike functionally equivalent to two bikes (I know, I find this silly as well). If we treat these as separate possibilities, we have four scenarios in total.\n",
    "\n",
    "For this to work, we will need to create **two** new arrays, one for each choice, and then tell the software to do combinatorial sampling:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c5af1e7-4063-46da-b3cc-3a25769e51c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_scenarios = bwp.create_datapackage(combinatorial=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4ad459f-29bd-4454-b51f-72dbc16a5356",
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_scenarios.add_persistent_array(\n",
    "    matrix='technosphere_matrix',\n",
    "    indices_array=np.array([(carbon_fibre, bike)], dtype=bwp.INDICES_DTYPE),\n",
    "    data_array=np.array([(2.5, 1.5)]),\n",
    "    flip_array=np.array([True]),\n",
    "    name='cf scenario'\n",
    ")\n",
    "dp_scenarios.add_persistent_array(\n",
    "    matrix='technosphere_matrix',\n",
    "    indices_array=np.array([(bike, bike)], dtype=bwp.INDICES_DTYPE),\n",
    "    data_array=np.array([(1, 2)]),\n",
    "    name='double bike'\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f75e957c-077a-481f-82e6-09dfd1007bbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "scenario_mapping = {\n",
    "    (0, 0): \"Original\",\n",
    "    (0, 1): \"Each bike counts double\",\n",
    "    (1, 0): \"Lightweight\",\n",
    "    (1, 1): \"Lightweight & each bike counts double\",\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13c7d614-db27-47bf-844f-090f5a5306db",
   "metadata": {},
   "outputs": [],
   "source": [
    "lca = bc.LCA(\n",
    "    demand={bike: 1},\n",
    "    data_objs=[dp_static, dp_scenarios],\n",
    "    use_arrays=True,\n",
    ")\n",
    "lca.lci()\n",
    "lca.lcia()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf559abe-39d1-410c-9410-9ea46c0c3ef2",
   "metadata": {},
   "outputs": [],
   "source": [
    "resource_group = next(grp for grp in lca.technosphere_mm.groups if grp.label == 'double bike').indexer.indexer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f85dcaef-f510-45c1-9a99-56436aa9c93b",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(lca.score, scenario_mapping[resource_group.index])\n",
    "\n",
    "for scenario_result in lca:\n",
    "    print(lca.score, scenario_mapping[resource_group.index])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f987fed9-e769-4940-845b-9780ea78b18c",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Can you build a half steel, half carbon fiber bike? Can you construct scenarios to vary the percentages (say, three scenarios)?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c511d688-6e26-4c86-a29a-ff3751eeebae",
   "metadata": {},
   "source": [
    "# Interfaces\n",
    "\n",
    "We can add some electricity consumption to the bike manufacturing process, and demonstrate how an `Interface` should work. For our example interface, we can use the Danish company Energinet's [API for the 5 minute interval CO2 intensity of Danish electricity](https://www.energidataservice.dk/tso-electricity/co2emis):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "005a5058-242c-4b96-874e-af37f840501f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "\n",
    "class DenmarkCO2Interface:\n",
    "    size = (1,)\n",
    "\n",
    "    def __next__(self):\n",
    "        # Supposed to work with ISO 8601 but T intervals raise error...\n",
    "        URL = \"https://api.energidataservice.dk/dataset/co2emis?start=now-P1D&end=now\"\n",
    "        print(\"Querying energidataservice.dk API\")\n",
    "        result = requests.get(URL).json()['records'][0]['CO2Emission']\n",
    "        return np.array([result]) / 1000 # g to kg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b20597e2-036c-4720-9585-ae7c60c18784",
   "metadata": {},
   "source": [
    "That was just the CO2 intensity, we also need to add a node for electricity consumption. We will give it the ID `104`, and use 100 kilowatt hours (way too high, I know)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68fcd3ca-5c52-457f-8827-f56c778fcb26",
   "metadata": {},
   "outputs": [],
   "source": [
    "electricity = 104"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d438e5fb-2a86-4dd6-a1d6-c02a4981d372",
   "metadata": {},
   "outputs": [],
   "source": [
    "dk_elec_dp = bwp.create_datapackage()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b38fc2a-288c-4081-88dd-89343aaac971",
   "metadata": {},
   "outputs": [],
   "source": [
    "dk_elec_dp.add_persistent_vector(\n",
    "    matrix='technosphere_matrix',\n",
    "    indices_array=np.array([\n",
    "            (electricity, electricity), # production of electricity\n",
    "            (electricity, bike), # use of electricity in bike production\n",
    "        ], \n",
    "        dtype=bwp.INDICES_DTYPE\n",
    "    ),\n",
    "    data_array=np.array([1, 100]),\n",
    "    flip_array=np.array([False, True]),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f610dba-8a85-41c5-a3f1-a81a7625bf08",
   "metadata": {},
   "source": [
    "We can now add our interface. We are not adding `persistent` data, but `dynamic` data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8bcd961-50de-4c6c-9dfa-5d16acc4ec44",
   "metadata": {},
   "outputs": [],
   "source": [
    "dk_elec_dp.add_dynamic_vector(\n",
    "    matrix='biosphere_matrix',\n",
    "    interface=DenmarkCO2Interface(),\n",
    "    indices_array=np.array([(co2, electricity)], dtype=bwp.INDICES_DTYPE),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccf82f27-477a-4513-a103-764b39ebb130",
   "metadata": {},
   "source": [
    "As before, we can just add the data packages together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee5e3e94-6833-48da-9c75-cdb501c07a92",
   "metadata": {},
   "outputs": [],
   "source": [
    "lca = bc.LCA(\n",
    "    demand={bike: 1},\n",
    "    data_objs=[dp_static, dk_elec_dp],\n",
    ")\n",
    "lca.lci()\n",
    "lca.lcia()\n",
    "lca.score"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8b5348e-d80d-451d-aacb-c4d7be6972ff",
   "metadata": {},
   "source": [
    "This is the value we got for CO2 intensity (in kg CO2-eq/kwh):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10bd8239-cd39-44ba-97eb-18be2ea27cda",
   "metadata": {},
   "outputs": [],
   "source": [
    "lca.biosphere_matrix[0, lca.dicts.activity[104]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0c9a1f2-eaa3-4c4c-ae00-7b9dcc998e18",
   "metadata": {},
   "source": [
    "Other countries also have similar APIs, for example France:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7001bbf3-44ee-45b2-bf9a-0234c6148947",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    The Interface below includes a parameter <code>verify=False</code> in the call to the <code>requests</code> library to pull the data.\n",
    "    \n",
    "    \n",
    "The TLS configuration from www.rte-france.com is incomplete, but the certificates are valid. You can verify the validity by pointing your browser to the <a href=\"https://www.rte-france.com/themes/swi/xml/power-co2-emission-fr.xml\"> API URL</a> <strong>. \n",
    "    We can safely use requests with: </strong> <code>\"verify=False\"</code>, but the <code>requests</code> library will still warn us about this.\n",
    "    \n",
    "More details on this at https://github.com/Depart-de-Sentier/try_brightway_dev/issues/3\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6b0ab2a-52f4-485d-badd-fb3bd4969fa4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from time import time\n",
    "import xml.etree.ElementTree as ET\n",
    "\n",
    "class FranceCO2Interface:\n",
    "    size = (1,)\n",
    "    \n",
    "    def __next__(self):\n",
    "        # Note: has timestamp but seems to always give latest data\n",
    "        URL_TEMPLATE = \"https://www.rte-france.com/themes/swi/xml/power-co2-emission-fr.xml?_={}\"\n",
    "        print(\"Querying RTE Eco2mix interface\")\n",
    "        URL = URL_TEMPLATE.format(int(time()))\n",
    "        # the TLS configuration from www.rte-france.com is incomplete, but the certificates are valid.\n",
    "        # To verify the validity, point your browser to the URL\n",
    "        # so, we can safely use requests with \"verify=False\"\n",
    "        # see issue on the repo of this notebook -> https://github.com/Depart-de-Sentier/try_brightway_dev/issues/3\n",
    "        # requests will still warn us about this.\n",
    "        resp = requests.get(URL, verify=False)        \n",
    "        root = ET.fromstring(resp.text)\n",
    "        return np.array([[float(child.text) for child in root.iter('valeur')][-1]]) / 1000 # g to kg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f2bcafa-e47b-44a5-80c9-1e662c7f4373",
   "metadata": {},
   "outputs": [],
   "source": [
    "fr_elec_dp = bwp.create_datapackage()\n",
    "\n",
    "fr_elec_dp.add_persistent_vector(\n",
    "    matrix='technosphere_matrix',\n",
    "    indices_array=np.array([\n",
    "            (electricity, electricity), # production of electricity\n",
    "            (electricity, bike), # use of electricity in bike production\n",
    "        ], \n",
    "        dtype=bwp.INDICES_DTYPE\n",
    "    ),\n",
    "    data_array=np.array([1, 100]),\n",
    "    flip_array=np.array([False, True]),\n",
    ")\n",
    "\n",
    "fr_elec_dp.add_dynamic_vector(\n",
    "    matrix='biosphere_matrix',\n",
    "    interface=FranceCO2Interface(),\n",
    "    indices_array=np.array([(co2, electricity)], dtype=bwp.INDICES_DTYPE),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32b67834-9770-4ffb-8b70-335a88584f02",
   "metadata": {},
   "source": [
    "As before, we can just add the data packages together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92075fec-0fdb-4e33-abad-ce53c71e54ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "lca = bc.LCA(\n",
    "    demand={bike: 1},\n",
    "    data_objs=[dp_static, fr_elec_dp],\n",
    ")\n",
    "lca.lci()\n",
    "lca.lcia()\n",
    "lca.score"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4375bc72-f325-44d1-b1fc-6ea8ccb025fa",
   "metadata": {},
   "source": [
    "This is the value we got for CO2 intensity (in kg CO2-eq/kwh):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b30507c-3481-4883-8bda-afdffc4795a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "lca.biosphere_matrix[0, lca.dicts.activity[electricity]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a2206a3-c55d-4089-ac2e-3c649a41a173",
   "metadata": {},
   "source": [
    "## Class Exercise\n",
    "\n",
    "How would we model a system where we could choose between Danish and French electricity as two scenarios?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}