3.250,3.430-Colors,SensorFailure.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\n",
    "using Statistics\n",
    "using LinearAlgebra\n",
    "using JSON\n",
    "\n",
    "include(\"readclassjson.jl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ; wget http://ee263.stanford.edu/hw/color_perception_data.json -P data/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ; wget http://ee263.stanford.edu/hw/one_bad_sensor.json -P data/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.250 \n",
    "\n",
    "*Human color perception is based on the responses of three different types of color light receptors, called cones. The three types of cones have different spectral-response characteristics, and are called L, M, and, S because they respond mainly to long, medium, and short wavelengths, respectively. In this problem we will divide the visible spectrum into 20 bands, and model the cones’ responses as follows:*\n",
    "\n",
    "$L_{\\text {cone }}=\\sum_{i=1}^{20} l_{i} p_{i}, \\quad M_{\\text {cone }}=\\sum_{i=1}^{20} m_{i} p_{i}, \\quad S_{\\text {cone }}=\\sum_{i=1}^{20} s_{i} p_{i}$\n",
    "\n",
    "*where $p_{i}$ is the incident power in the $i$ th wavelength band, and $l_{i}, m_{i}$ and $s_{i}$ are nonnegative constants that describe the spectral responses of the different cones. The perceived color is a complex function of the three cone responses, i.e., the vector $\\left(L_{\\text {cone }}, M_{\\text {cone }}, S_{\\text {cone }}\\right),$ with different cone response vectors perceived as different colors. (Actual color perception is a bit more complicated than this, but the basic idea is right.)*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### a)  Metamers\n",
    "*When are two light spectra, $p$ and $\\tilde{p},$ visually indistinguishable? (Visually identical lights with different spectral power compositions are called metamers.)*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the vector $\\left(L_{\\text {cone }}, M_{\\text {cone }}, S_{\\text {cone }}\\right),$ determines the perception of the color, two light spectra are indistinguishable when they produce the same values in this vector. If we construct a matrix $A = (l, m, s)^T$, then we are saying that A is not 1-1 and has a null space larger than {0} since it can map multiple input vectors $p$ to the same output in the equation y = Ap. Essentially, p and $\\tilde{p}$ will be visually indistinguishable when $(\\tilde{p}-p)\\in null(A)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### b) Visual color matching\n",
    "*In a color matching problem, an observer is shown a test light, and is asked to change the intensities of three primary lights until the sum of the primary lights looks like the test light. In other words, the observer is asked the find a spectrum of the form $p_{\\text {match }}=a_{1} u+a_{2} v+a_{3} w$ where $u, v, w$ are the spectra of the primary lights, and $a_{i}$ are the intensities to be found, that is visually indistinguishable from a given test light spectrum $p_{\\text {test. }}$ Can this always be done? Discuss briefly.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a control problem as we are adjusting inputs to match a specific desired output. If we write the matching vector as a matrix equation $p_{match}  = \\underset{20 \\times 3}{B}\\quad\\underset{3 \\times 1}{a} =  \\begin{bmatrix} u & v & w\\end{bmatrix}\\begin{bmatrix} a_1 \\\\ a_2 \\\\ a_3 \\end{bmatrix}$ then we can see that we can only achieve $p_{match}$ that are in the range of B. Similarly, we can say $p_{match}$ will be perceived as $Ap_{match}$ for A defined as above. Then we can achieve visual perceptions that are in the range of $AB$. We saw that $range(AB) \\subset range(A)$, so we can always match the test light if $R(A) = R(AB)$ since the test light will be perceived as $Ap_{test}$ and primary lights $ABa$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### c) Visual matching with phosphors\n",
    "*A computer monitor has three phosphors, R, G, and B.  It is desired to adjust the phosphor intensities to create a color that looks like a reference test light.  Find weights that achieve the match or explain why no such weights exist. The data for this problem is in color_perception_data.json, which contains the vectors wavelength,B_phosphor,G_phosphor,R_phosphor,L_coefficients,M_coefficients,S_coefficients, and test_light.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dict{Any,Any} with 10 entries:\n",
       "  \"test_light\"     => [58.2792, 42.3496, 51.5512, 33.3951, 43.2907, 22.595, 57.…\n",
       "  \"wavelength\"     => [380.0, 400.0, 420.0, 440.0, 460.0, 480.0, 500.0, 520.0, …\n",
       "  \"tungsten\"       => [20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 1…\n",
       "  \"L_coefficients\" => [0.12589, 0.079433, 0.1, 0.15849, 0.19953, 0.25119, 0.316…\n",
       "  \"M_coefficients\" => [0.15849, 0.14125, 0.14125, 0.25119, 0.31623, 0.39811, 0.…\n",
       "  \"G_phosphor\"     => [21.0, 23.0, 26.0, 30.0, 35.0, 45.0, 75.0, 90.0, 100.0, 9…\n",
       "  \"R_phosphor\"     => [15.0, 16.0, 17.0, 18.0, 19.0, 21.0, 23.0, 26.0, 30.0, 35…\n",
       "  \"B_phosphor\"     => [30.0, 35.0, 45.0, 75.0, 90.0, 100.0, 90.0, 75.0, 45.0, 3…\n",
       "  \"S_coefficients\" => [0.39811, 0.50119, 0.63096, 0.79433, 0.50119, 0.25119, 0.…\n",
       "  \"sunlight\"       => [40.0, 70.0, 100.0, 160.0, 240.0, 220.0, 180.0, 160.0, 18…"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "colors = readclassjson(\"data/color_perception_data.json\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define relevant matrices and vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = colors[\"test_light\"]\n",
    "w = colors[\"wavelength\"]\n",
    "rgb = [colors[\"R_phosphor\"] colors[\"G_phosphor\"] colors[\"B_phosphor\"]]\n",
    "lms = [colors[\"L_coefficients\"] colors[\"M_coefficients\"] colors[\"S_coefficients\"]];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The vector y is the perceived color of the test light"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 191.16611560329102\n",
       " 242.7872917703638\n",
       " 135.92624672925962"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = transpose(lms) * test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matrix A is the perceived colors from the RGB spectra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " 145.555   254.398   197.777\n",
       " 148.341   333.869   278.369\n",
       "  55.7972   95.5058  194.713"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = transpose(lms) * rgb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By taking the inverse of A, we can reconstruct the intensity pattern needed to produce y exactly. The vector x is what we set out to find - weights that achieve the match"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 0.4225929938835502\n",
       " 0.09874256378706159\n",
       " 0.5285525473174244"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "invA = inv(A)\n",
    "x = invA * y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having constructed x from the inverse, $Ax$ then should give us y exactly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 191.16611560329105\n",
       " 242.78729177036385\n",
       " 135.92624672925962"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A*x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testing for equivalence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isapprox(A*x,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### d)\n",
    "Effects of illumination. An object's surface can be characterized by its reflectance ( $i . e .,$ the fraction of light it reflects) for each band of wavelengths. If the object is illuminated with a light spectrum characterized by $I_{i},$ and the reflectance of the object is $r_{i}($ which is between 0 and 1 ), then the reflected light spectrum is given by $I_{i} r_{i},$ where $i=1, \\ldots, 20$ denotes the wavelength band. Now consider two objects illuminated (at different times) by two different light sources, say an incandescent bulb and sunlight. Sally argues that if the two objects look identical when illuminated by a tungsten bulb, then they will look identical when illuminated by sunlight. Beth disagrees: she says that two objects can appear identical when illuminated by a tungsten bulb, but look different when lit by sunlight. Who is right? If Sally is right, explain why. If Beth is right give an example of two objects that appear identical under one light source and different under another. You can use the vectors sunlight and tungsten defined in the data file as the light\n",
    "sources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "tungsten = colors[\"tungsten\"]\n",
    "sunlight = colors[\"sunlight\"];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "B = transpose(lms);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start with some reflectance vector $r1$. By taking the elementwise product with the tungsten illumination, we produce $I_tr_1$, the reflected light vector for our object 1 with reflectance of $r1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "r1 = [repeat([0.5], 10) ; repeat([0.1], 10)]\n",
    "IR1 = (tungsten .* r1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "nullVector = nullspace(B)[:,1];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By finding vectors in the null space of B, we can add these vectors to our input to achieve the same output. Therefore we are guaranteed to construct another object 2 with different reflectance $r2$ that after an elementwise product with tungsten, B maps to the same output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "r2 = (IR1 - 50 .* nullVector) ./ tungsten\n",
    "IR2 = (tungsten .* r2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reflectances are different"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "false"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isapprox(r1, r2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reflected light spectra under tungsten are also different"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "false"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isapprox(IR2, IR1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate how reflected light spectra are perceived under tungsten"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 242,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 124.82236869999998\n",
       " 173.2441382523\n",
       "  71.23471167566724"
      ]
     },
     "execution_count": 242,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tungstenLook1 = B * IR1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 243,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 124.82236869999997\n",
       " 173.2441382523\n",
       "  71.23471167566724"
      ]
     },
     "execution_count": 243,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tungstenLook2 = B * IR2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They are perceived as the same"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 244,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isapprox(tungstenLook1, tungstenLook2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now turn to sunlight. We take the same reflexivity properties we created before, but now calculate their reflected spectra under the new light source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 245,
   "metadata": {},
   "outputs": [],
   "source": [
    "sun_IR1 = (sunlight .* r1);\n",
    "sun_IR2 = (sunlight .* r2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we use B to map the reflected spectra to their perceptions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 247,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 247.85898313\n",
       " 359.8846318979\n",
       " 215.66132988271158"
      ]
     },
     "execution_count": 247,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sunLook1 = B * sun_IR1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 248,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 246.60513901099375\n",
       " 357.20225222205795\n",
       " 196.89366873222457"
      ]
     },
     "execution_count": 248,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sunLook2 = B * sun_IR2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 249,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "false"
      ]
     },
     "execution_count": 249,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isapprox(sunLook1, sunLook2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The objects that looked the same under tungsten do not appear the same under sunlight. Using objects $r1, r2$ we can see Beth is correct, as they appear identical under one light source and different under another."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.430. Single sensor failure detection and identification\n",
    "Single sensor failure detection and identification. We have $y=A x,$ where $A \\in \\mathbb{R}^{m \\times n}$ is known, and $x \\in \\mathbb{R}^{n}$ is to be found. Unfortunately, up to one sensor may have failed (but you don't know which one has failed, or even whether any has failed). You are given $\\tilde{y}$ and not $y,$ where $\\tilde{y}$ is the same as $y$ in all entries except, possibly, one (say, the $k$ th entry). If all sensors are operating correctly, we have $y=\\tilde{y} .$ If the $k$ th sensor fails, we have $\\tilde{y}_{i}=y_{i}$ for all $i \\neq k$\n",
    "The file one_bad_sensor.m, available on the course web site, defines $A$ and $\\tilde{y}$ (as A and ytilde). Determine which sensor has failed (or if no sensors have failed). You must explain your method, and submit your code.\n",
    "\n",
    "For this exercise, you can use the matlab code rank ( $[\\mathrm{F}$ g] ) = rank (F) to check if $g \\in$ range $(F) .$ (We will see later a much better way to check if $g \\in \\operatorname{range}(F) .)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sensor = readclassjson(\"data/one_bad_sensor.json\")\n",
    "A = sensor[\"A\"]\n",
    "ytilde = sensor[\"ytilde\"];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(15, 10)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "size(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we can check whether there are any bad sensors. We see below that the rank of A differs from the rank of $\\begin{bmatrix} A & ytilde\\end{bmatrix}$. This tells us that ytilde is not a linear combination of the columns of A, since otherwise it would not change the rank when appended to A. This tells us the current measurements in ytilde are not in the range of A and are an invalid output of the system y = Ax. Therefore we have a bad sensor and need to find it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rank(A) = 10\n",
      "rank(hcat(A, ytilde)) = 11\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "11"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@show rank(A)\n",
    "@show rank(hcat(A,ytilde))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our method to find the bad sensor is simple. First, we know there is only a single bad sensor. Second, we know that each row of the vector y is mapped from the corresponding row in the matrix A. No other rows are used to produce that entry. Therefore we can find the bad sensor by removing the same indexed row from y and A and checking to see if the reduced y, y_prime, is in the range of the reduced A, A_prime. If $y_{prime} \\notin R(A_{prime})$ then the row we removed was not the bad sensor. If however $y_{prime} \\in R(A_{prime})$, then the row we removed corresponds to the bad sensor, since y_prime is now a feasible output as a linear combination of the columns of A_prime. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11\n"
     ]
    }
   ],
   "source": [
    "for i = 1:size(A)[1]\n",
    "    \n",
    "    idx = 1:15 .!= i\n",
    "    A_prime = A[idx, :]\n",
    "    y_prime = ytilde[idx]\n",
    "\n",
    "    if rank(hcat(A_prime,y_prime)) == rank(A_prime)\n",
    "        println(i)\n",
    "    end\n",
    "end\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the 11th row of y and A were removed, $y_{prime} \\in R(A_{prime})$ meaning the 11th sensor is the producing the bad readings."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.1",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}