{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# BEE 4750 Homework 1: Introduction to Using Julia\n",
        "\n",
        "**Name**:\n",
        "\n",
        "**ID**:\n",
        "\n",
        "> **Due Date**\n",
        ">\n",
        "> Thursday, 9/03/26, 9:00pm\n",
        "\n",
        "## Overview\n",
        "\n",
        "### Load Environment\n",
        "\n",
        "The following code loads the environment and makes sure all needed\n",
        "packages are installed. This should be at the start of most Julia\n",
        "scripts."
      ],
      "id": "1ffa5682-32fd-4954-945a-0a7758f69158"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "import Pkg\n",
        "Pkg.activate(@__DIR__)\n",
        "Pkg.instantiate()"
      ],
      "id": "2"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Standard Julia practice is to load all needed packages at the top of a\n",
        "file. If you need to load any additional packages in any assignments\n",
        "beyond those which are loaded by default, feel free to add a `using`\n",
        "statement, though [you may need to install the\n",
        "package](https://viveks.me/environmental-systems-analysis/tutorials/julia-basics.html#package-management)."
      ],
      "id": "9c09e58b-3a42-4116-b4c1-a912bce8398f"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "using Random # allows for random seed generation\n",
        "using Plots # basic plotting: can also install and use other packages\n",
        "using Graphs # for making graphs and networks\n",
        "using GraphRecipes # basic graph generation\n",
        "using LaTeXStrings # allows for LaTeX formatting in plots\n",
        "using Distributions # sampling and fitting of probability distributions\n",
        "using CSV # file I/O for CSV files\n",
        "using DataFrames # data structure for tabular data"
      ],
      "id": "4"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "# this sets a random seed, which ensures reproducibility of random number generation. You should always set a seed when working with random numbers.\n",
        "Random.seed!(1)"
      ],
      "id": "6"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Problems (Total: 50 Points)\n",
        "\n",
        "### Problem 1 (10)\n",
        "\n",
        "***Systems Modeling***:\n",
        "\n",
        "Three factories are discharging a chemical called Chlororadiated\n",
        "Ureadicarboxyl (or CRUD) into the Riley River. The river inflow has a\n",
        "volume of 500,000 m$^3$/d and a CRUD concentration of 0.2 mg/L. The\n",
        "river velocity is 25 km/d.\n",
        "\n",
        "- Factory 1 discharges $100,000 \\text{m}^3$/ of effluent with a CRUD\n",
        "  concentration of 10 mg/L;\n",
        "- Factory 2 discharges $60,000 \\text{m}^3$/ of effluent with a CRUD\n",
        "  concentration of 20 mg/L;\n",
        "- Factory 3 discharges $200,000 \\text{m}^3$/ of effluent with a CRUD\n",
        "  concentration of 8 mg/L.\n",
        "\n",
        "CRUD decays in the river with a first order reaction rate\n",
        "$k=0.45 \\text{d}^{-1}$.\n",
        "\n",
        "Environmental authorities have taken samples downriver of the three\n",
        "factories and found that the CRUD concentration exceeds the regulatory\n",
        "limit of 1 mg/L. CRUD can be removed at the source discharge at a cost\n",
        "of $50E^2$ per $1000 \\text{m}^3$ treated, where $E$ is the efficiency of\n",
        "removal (the fraction of CRUD removed).\n",
        "\n",
        "#### Problem 1.1 (2)\n",
        "\n",
        "Draw a diagram of the river and factory system. Include definitions for\n",
        "any notation you will use.\n",
        "\n",
        "#### Problem 1.2 (5)\n",
        "\n",
        "Write down a set of conditions on the removal efficiencies $E_1$, $E_2$,\n",
        "and $E_3$ (where the subscript corresponds to the factory where the\n",
        "effluent is being treated) to guarantee compliance.\n",
        "\n",
        "#### Problem 2.3 (3)\n",
        "\n",
        "Without finding an “optimal” solution, reason about which factories\n",
        "ought to treat their effluent at a higher level than the others. This\n",
        "can draw on considerations such as cost, the sequencing of discharges,\n",
        "and the amount of CRUD discharged by each factory, or others.\n",
        "\n",
        "### Problem 2 (6)\n",
        "\n",
        "***Systems Perturbations***:\n",
        "\n",
        "Consider the perturbed phosphorous cycle example from Lecture 2. Now\n",
        "suppose that the initial perturbation at $t=0$ is a 10% reduction in the\n",
        "rate constant for the conversion of phosphorous from organic to\n",
        "inorganic form (through an inhibition of phosphatase activity due to\n",
        "e.g. acidic rain fall). When the new steady state is reached, how much\n",
        "phosphorous will be in each compartment?\n",
        "\n",
        "### Problem 3 (9)\n",
        "\n",
        "***Debugging Code***:\n",
        "\n",
        "The following subproblems all involve code snippets that require\n",
        "debugging. You are encouraged to use online resources (*e.g.* Julia\n",
        "documentation and forums, Stack Overflow, Reddit, etc) to help find\n",
        "diagnose error messages and find solutions, but make sure that you\n",
        "clearly document which resources you used and how you used them.\n",
        "\n",
        "#### Problem 3.1 (3)\n",
        "\n",
        "You’ve been tasked with writing code to identify the minimum value in an\n",
        "array. You cannot use a predefined function. Your colleague suggested\n",
        "the function below, but it does not return the minimum value."
      ],
      "id": "2bbf6031-76ab-44e8-a7be-5ff14ccdd46f"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "minimum(array_values) = 0"
          ]
        }
      ],
      "source": [
        "function minimum(array)\n",
        "    # initialize the minimum value counter\n",
        "    min_value = 0\n",
        "    # update minimum values\n",
        "    for i in 1:length(array)\n",
        "        if array[i] < min_value\n",
        "            min_value = array[i]\n",
        "        end\n",
        "    end\n",
        "    # return found minimum\n",
        "    return min_value\n",
        "end\n",
        "\n",
        "array_values = [89, 90, 95, 100, 100, 78, 99, 98, 100, 95]\n",
        "@show minimum(array_values);"
      ],
      "id": "8"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "What is the logical flaw in the code? Fix it (describe what you changed\n",
        "rather than reproducing the entire code snippet) and show that your\n",
        "fixed code solves the problem.\n",
        "\n",
        "#### Problem 3.2 (3)\n",
        "\n",
        "Your team wants to know the expected payout of an old Italian dice game\n",
        "called *passadieci* (which was analyzed by Galileo as one of the first\n",
        "examples of a rigorous study of probability). The goal of passadieci is\n",
        "to get at least an 11 from rolling three fair, six-sided dice. Your\n",
        "strategy is to compute the average wins from 1,000 trials, but the code\n",
        "you’ve written below produces a `MethodError`."
      ],
      "id": "e5b3bf70-5b76-45dc-b141-a8dd30d9b6d7"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "# function to simulate a passadieci roll: 3 6-sided dice\n",
        "function passadieci()\n",
        "    # this rand() call samples 3 values from the vector [1, 6]\n",
        "    roll = rand(1:6, 3) \n",
        "    return roll\n",
        "end\n",
        "# set number of trials and initialize outcome vector\n",
        "n_trials = 1_000\n",
        "outcomes = zero(n_trials)\n",
        "# simulate number of passadieci rolls and count wins\n",
        "for i = 1:n_trials\n",
        "    outcomes[i] = (sum(passadieci()) > 11)\n",
        "end\n",
        "win_prob = sum(outcomes) / n_trials # compute average number of wins\n",
        "@show win_prob;"
      ],
      "id": "10"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "What does this `MethodError` mean? Why does this code produce one? Fix\n",
        "the error (describe what you changed rather than reproducing the entire\n",
        "code snippet) and solve the problem.\n",
        "\n",
        "#### Problem 3.3 (3)\n",
        "\n",
        "You’re interested in writing some code to remove the mean of a vector\n",
        "from all of its components. You’ve written the following code and tried\n",
        "to test it on a random vector, but your code returns a `MethodError`."
      ],
      "id": "099b4ef8-c6a5-4b07-a1b0-316ec56ce063"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "# function to remove mean from a vector\n",
        "function remove_mean(vect)\n",
        "    # fucntion to compute the mean\n",
        "    function compute_mean(vect)\n",
        "        element_sum = 0 # initialize sum\n",
        "        # compute mean and return\n",
        "        for v in vect\n",
        "            element_sum += v\n",
        "        end\n",
        "        return element_sum / length(vect)\n",
        "    end\n",
        "\n",
        "    m = compute_mean(vect) # compute mean\n",
        "    # return demeaned vector\n",
        "    return vect - m\n",
        "end\n",
        "\n",
        "random_vect = rand(1_000)\n",
        "@show remove_mean(random_vect)"
      ],
      "id": "12"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "What does this `MethodError` mean? Why does this code produce one? Fix\n",
        "the error (describe what you changed rather than reproducing the entire\n",
        "code snippet) and solve the problem.\n",
        "\n",
        "### Problem 4 (5)\n",
        "\n",
        "***Loading and Manipulating Data***\n",
        "\n",
        "The file `fha.csv` contains data[1] on the size of American metro areas\n",
        "and the number of miles people drive in each metro area each day. This\n",
        "problem asks you to load and work with this data to get practice with\n",
        "data manipulation in Julia.\n",
        "\n",
        "#### Problem 4.1 (2)\n",
        "\n",
        "Load `fha.csv` into a Julia `DataFrame`. How many rows and columns does\n",
        "the data have?\n",
        "\n",
        "#### Problem 4.2 (3)\n",
        "\n",
        "Calculate the number of miles driven per person per day for every metro\n",
        "area. Make a scatterplot of this on the $y$-axis and the population of\n",
        "the city on the $x$-axis. Make sure to label your axes and include a\n",
        "caption clearly describing the plot. Highlight where the New York City\n",
        "and Houston metro areas are on the plot and make a legend clearly\n",
        "labelling any plot features.\n",
        "\n",
        "### Problem 5 (20)\n",
        "\n",
        "***Computational Modeling***:\n",
        "\n",
        "Cheap Plastic Products, Inc. is operating a plant that produces\n",
        "$100 \\text{m}^3\\text{/day}$ of wastewater that is discharged into\n",
        "Pristine Brook. The wastewater contains $1 \\text{kg/m}^3$ of YUK, a\n",
        "toxic substance. The US Environmental Protection Agency has imposed an\n",
        "effluent standard on the plant prohibiting discharge of more than\n",
        "$20 \\text{kg/day}$ of YUK into Pristine Brook.\n",
        "\n",
        "Cheap Plastic Products has analyzed two methods for reducing its\n",
        "discharges of YUK. Method 1 is land disposal, which costs $X_1^2/20$\n",
        "dollars per day, where $X_1$ is the amount of wastewater disposed of on\n",
        "the land ($\\text{m}^3\\text{/day}$). With this method, 20% of the YUK\n",
        "applied to the land will eventually drain into the stream (*i.e.*, 80%\n",
        "of the YUK is removed by the soil).\n",
        "\n",
        "Method 2 is a chemical treatment procedure which costs \\$1.50 per\n",
        "$\\text{m}^3$ of wastewater treated. The chemical treatment has an\n",
        "efficiency of $e= 1 - 0.005X_2$, where $X_2$ is the quantity of\n",
        "wastewater ($\\text{m}^3\\text{/day}$) treated. For example, if\n",
        "$X_2 = 50 \\text{m}^3\\text{/day}$, then $e = 1 - 0.005(50) = 0.75$, so\n",
        "that 75% of the YUK is removed.\n",
        "\n",
        "Cheap Plastic Products is wondering how to allocate their wastewater\n",
        "between these three disposal and treatment methods (land disposal,\n",
        "chemical treatment, and direct disposal) to meet the effluent standard\n",
        "while keeping costs manageable. The flow of wastewater through this\n",
        "treatment system is shown in\n",
        "<a href=\"#fig-wastewater\" class=\"quarto-xref\">Figure 1</a>.\n",
        "\n",
        "[1] From the [2023 Federal Highway Administration’s Highway Statistics\n",
        "report](https://www.fhwa.dot.gov/policyinformation/statistics/2023/pdf/hm72.pdf)"
      ],
      "id": "2b050b77-6665-4976-9ac1-5e3dab292458"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "A = [0 1 1 1;\n",
        "    0 0 0 1;\n",
        "    0 0 0 1;\n",
        "    0 0 0 0]\n",
        "\n",
        "names = [\"Plant\", \"Land Treatment\", \"Chem Treatment\", \"Pristine Brook\"]\n",
        "# modify this dictionary to add labels\n",
        "edge_labels = Dict((1, 2) => \"\", (1,3) => \"\", (1, 4) => \"\",(2, 4) => \"\",(3, 4) => \"\")\n",
        "shapes=[:hexagon, :rect, :rect, :hexagon]\n",
        "xpos = [0, -1.5, -0.25, 1]\n",
        "ypos = [1, 0, 0, -1]\n",
        "\n",
        "p = graphplot(A, names=names,edgelabel=edge_labels, markersize=0.15, markershapes=shapes, markercolor=:white, x=xpos, y=ypos)\n",
        "display(p)"
      ],
      "id": "cell-fig-wastewater"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "#### Problem 5.1 (2)\n",
        "\n",
        "Make a figure showing notation for how YUK moves through the system for\n",
        "a given treatment allocation. Make sure to define your variables. You\n",
        "can make this by modifying the graph above or with an independent figure\n",
        "included in your solution.\n",
        "\n",
        "You can modify the edge labels (by editing the `edge_labels` dictionary\n",
        "in the code) in\n",
        "<a href=\"#fig-wastewater\" class=\"quarto-xref\">Figure 1</a> to show how\n",
        "the wastewater allocations result in the final YUK discharge into\n",
        "Pristine Brook. For the `edge_label` dictionary, the tuple $(i, j)$\n",
        "corresponds to the arrow going from node $i$ to node $j$. The syntax for\n",
        "any entry is `(i, j) => \"label text\"`, and the label text can include\n",
        "mathematical notation if the string is prefaced with an `L`, as in\n",
        "`L\"x_1\"` will produce $x_1$.\n",
        "\n",
        "#### Problem 5.2 (4)\n",
        "\n",
        "Formulate a mathematical model for the treatment cost and the amount of\n",
        "YUK that will be discharged into Pristine Brook based on the wastewater\n",
        "allocations. This is best done with some equations and supporting text\n",
        "explaining the derivation. Make sure you include, as additional\n",
        "equations in the model, any needed constraints on relevant values.\n",
        "\n",
        "> **None**\n",
        ">\n",
        "> You can find some basics on writing mathematical equations using the\n",
        "> LaTeX typesetting syntax\n",
        "> [here](https://viveks.me/environmental-systems-analysis/tutorials/latex-notebook.qmd),\n",
        "> and a cheatsheet with LaTeX commands can be found on the course\n",
        "> website’s [Resources\n",
        "> page](https://viveks.me/environmental-systems-analysis/resources/markdown.qmd).\n",
        "\n",
        "#### Problem 5.3 (3)\n",
        "\n",
        "Implement your systems model as a Julia function which computes the\n",
        "resulting YUK discharge and cost for a particular treatment plan. You\n",
        "can return multiple values from a function with a\n",
        "[tuple](https://docs.julialang.org/en/v1/manual/functions/#Tuples-1), as\n",
        "in:"
      ],
      "id": "89b23b68-626a-44cb-8238-bac81392f5b1"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "a = 7\n",
            "b = 10"
          ]
        }
      ],
      "source": [
        "function multiple_return_values(x, y)\n",
        "    return (x+y, x*y)\n",
        "end\n",
        "\n",
        "a, b = multiple_return_values(2, 5)\n",
        "@show a;\n",
        "@show b;"
      ],
      "id": "16"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To evalute the function over vectors of inputs, you can *broadcast* the\n",
        "function by adding a decimal `.` before the function arguments and\n",
        "accessing the resulting values by writing a *comprehension* to loop over\n",
        "the individual outputs in the vector:"
      ],
      "id": "aa70b269-e2df-4fc2-b595-58444c651730"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "a = [7, 9, 11, 13, 15]\n",
            "b = [6, 14, 24, 36, 50]"
          ]
        }
      ],
      "source": [
        "x = [1, 2, 3, 4, 5]\n",
        "y = [6, 7, 8, 9, 10]\n",
        "\n",
        "output = multiple_return_values.(x, y)\n",
        "a = [out[1] for out in output]\n",
        "b = [out[2] for out in output]\n",
        "@show a;\n",
        "@show b;"
      ],
      "id": "18"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Make sure you comment your code appropriately to make it clear what is\n",
        "going on and why.\n",
        "\n",
        "Use your function to experiment with 1,000 different combinations of\n",
        "wastewater discharge and treatment. You can do this with either a grid\n",
        "search or by sampling from a [Dirichlet\n",
        "distribution](https://en.wikipedia.org/wiki/Dirichlet_distribution) (a\n",
        "$\\text{Dirichlet}(1, n)$ distribution will generate uniformly-weighted\n",
        "$n$-dimensional vectors whose components add up to 1; see the\n",
        "[`Distributions.jl`\n",
        "documentation](https://juliastats.org/Distributions.jl/stable/starting/)\n",
        "for how to sample from probability distributions in Julia).\n",
        "\n",
        "Plot the results of these experiments by making a scatterplot of YUK\n",
        "discharge into the lake vs. treatment cost. Make sure to label your axes\n",
        "and include a legend, as well as a caption clearly describing the plot.\n",
        "\n",
        "#### Problem 5.4 (3)\n",
        "\n",
        "Do any satisfy the YUK effluent standard (plot this as well as a dashed\n",
        "red line). What was the cost of solutions satisfying the standard?\n",
        "\n",
        "#### Problem 5.5 (4)\n",
        "\n",
        "What can you say about the tradeoff between treatment cost and YUK\n",
        "discharge? You don’t have to find an “optimal” solution to this problem,\n",
        "but what do you think would be needed to find a better solution?\n",
        "\n",
        "#### Problem 5.6 (4)\n",
        "\n",
        "Find the strategies which minimize cost and YUK discharge (these will be\n",
        "different strategies) analytically and find the values of the objective\n",
        "metrics. Plot these values in the plot that you created for Problem 4.3.\n",
        "How do their values compare to the spread of values that you found in\n",
        "that problem?\n",
        "\n",
        "## References\n",
        "\n",
        "List any external references consulted, including classmates."
      ],
      "id": "9fa773ce-749d-4b38-ba7f-8d84ae544e15"
    }
  ],
  "nbformat": 4,
  "nbformat_minor": 5,
  "metadata": {
    "kernel_info": {
      "name": "julia"
    },
    "kernelspec": {
      "name": "julia",
      "display_name": "Julia",
      "language": "julia"
    },
    "language_info": {
      "name": "julia",
      "codemirror_mode": "julia",
      "version": "1.11.5"
    }
  }
}