MRI.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MRI\n",
    "\n",
    "The following problem was posted on the [Julia Discourse](https://discourse.julialang.org/t/1e6-constrained-optimization-problems/):\n",
    "\n",
    "> So I have three magnetic resonance (MR) images: $X, Y, Z$. They come from three different [pulse sequences](https://en.wikipedia.org/wiki/MRI_sequence) which interact with [NMR relaxation parameters ](https://en.wikipedia.org/wiki/Relaxation_(NMR)) in known ways. That is, we have equations which given the underlying NMR parameters (specifically, T1, T2 and PD) and the scanner parameters $p^X, p^Y, p^Z$ corresponding to each image $X,Y,Z$, we can model the output intensity value observed in the acquired image at some voxel [$(i, j, k)$, e.g., $X_{i,j,k}$, since we are in 3D].\n",
    "\n",
    "> However, we don’t know the $T_1$, $T_2$ and $PD$ parameters but we do have $X,Y,Z$ and we can separately solve for each $p^{(\\cdot)}$. If you haven’t followed any of the previous mediocre explanation, then the problem simply reduces to the following system of equations, for which I am trying to solve $T_1$, $T_2$, and $P_D$ where everything else is given.\n",
    "\n",
    "> $$\n",
    " \\begin{align}\n",
    " \\log(X_{i,j,k}) &= \\log(PD) + p_0^X + p_1^X T_1 + p_2^X T_1^2 \\\\\n",
    "\\log(Y_{i,j,k}) &= \\log(PD) + p_0^Y + p_1^Y T_1 + p_2^Y / T_2 \\\\\n",
    "\\log(Z_{i,j,k}) &= \\log(PD) + p_0^Z + p_1^Z T_1 + p_2^Z / T_2 \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    ">In previous iterations of this project, the person who built the code used non-linear least squares (in MATLAB, FWIW) to solve for the values of $T_1$, $T_2$ and $PD$. However, the person who built the code left our lab, and with that has gone the knowledge of why that was chosen.\n",
    "\n",
    ">Previously I was using python’s [scipy implementation of least squares ](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html) to solve the problem at every voxel; however, even with parallelization this *takes hours*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since these are polynomial equations we can use homotopy continuation to solve the problem. Let's setup the system first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Polynomial{true,Int64},1}:\n",
       " p₁₋₃T₁² + p₁₋₂T₁ - x + PD + p₁₋₁     \n",
       " p₂₋₂T₁T₂ - yT₂ + PDT₂ + p₂₋₁T₂ + p₂₋₃\n",
       " p₃₋₂T₁T₂ - zT₂ + PDT₂ + p₃₋₁T₂ + p₃₋₃"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using HomotopyContinuation, DynamicPolynomials\n",
    "@polyvar x y z PD p[1:3,1:3] T[1:2]\n",
    "\n",
    "f1 = PD + p[1,1] + p[1,2]*T[1] + p[1,3] * T[1]^2 - x\n",
    "# note that we eliminate the denominator in equation 2 and 3\n",
    "f2 = (PD + p[2,1] + p[2,2]*T[1])*T[2] + p[2,3] - y * T[2]\n",
    "f3 = (PD + p[3,1] + p[3,2]*T[1])*T[2] + p[3,3] - z * T[2]\n",
    "\n",
    "f = [f1, f2, f3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we need to solve the same instance many many times it makes sense to first compute the generic number of solutions of this polynomial system and then to use a [parameter homotopy](https://www.juliahomotopycontinuation.org/guides/parameter-homotopies/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need some generic points to generate a generic instance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Polynomial{true,Complex{Float64}},1}:\n",
       " (0.06560018192766423 + 1.226509609243317im)T₁² + PD + (-0.40076878384840564 + 0.3936976384006835im)T₁ + (-0.4034430867072931 + 2.130445251499418im)   \n",
       " PDT₂ + (0.17548682027400092 + 1.0868781469710116im)T₁T₂ + (1.0068779513453254 - 0.9349141017423815im)T₂ + (0.6579055458504889 + 0.11880590069556508im)\n",
       " PDT₂ + (-0.6911103187839271 - 1.1844653755707422im)T₁T₂ + (0.216353296242008 - 0.3087205812565972im)T₂ + (-1.3312172733617187 + 1.9140562766862028im) "
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params = [vec(p); x; y; z]\n",
    "generic_params = randn(ComplexF64, length(params)) # Note we should sample them randomly\n",
    "f_generic = subs.(f, Ref(params => generic_params))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AffineResult with 8 tracked paths\n",
       "==================================\n",
       "• 2 non-singular finite solutions (0 real)\n",
       "• 0 singular finite solutions (0 real)\n",
       "• 6 solutions at infinity\n",
       "• 0 failed paths\n",
       "• random seed: 808747\n"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "generic_result = solve(f_generic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the system has only 2 solutions (instead the maximal possible 8 solutions)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Array{Complex{Float64},1},1}:\n",
       " [0.100221+0.505434im, -0.256837+1.29393im, 0.809632-0.787296im]\n",
       " [-2.12345+1.73512im, 0.438372-2.01187im, -0.398561+0.217356im] "
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "generic_solutions = solutions(generic_result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's prepare to solve our actual problems:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PathTracker()"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Construct a path tracker to track the two solutions to our specific instances\n",
    "# p₁ and \n",
    "tracker = pathtracker(f, generic_solutions;\n",
    "            parameters=params,\n",
    "            # we just use som dummy parameters for now\n",
    "            targetparameters=generic_params,\n",
    "            startparameters=generic_params,\n",
    "            affine=true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Out of lack of real world data let's sample some random parameters again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "specific_params = randn(length(params));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the parameters order is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12-element Array{PolyVar{true},1}:\n",
       " p₁₋₁\n",
       " p₂₋₁\n",
       " p₃₋₁\n",
       " p₁₋₂\n",
       " p₂₋₂\n",
       " p₃₋₂\n",
       " p₁₋₃\n",
       " p₂₋₃\n",
       " p₃₋₃\n",
       " x   \n",
       " y   \n",
       " z   "
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PathTrackerResult{Array{Complex{Float64},1}}:\n",
       " • returncode → success\n",
       " • x → Complex{Float64}[0.454103+9.56595e-24im, 1.3872-3.60672e-23im, 0.132065+1.31272e-24im]\n",
       " • t → 0.0 + 0.0im\n",
       " • accuracy → 2.7287366859420435e-15\n",
       " • accepted_steps → 17\n",
       " • rejected_steps → 2\n"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Setup the new parameters\n",
    "set_parameters!(tracker.homotopy, generic_params, specific_params);\n",
    "\n",
    "# track to the two solutions\n",
    "r1 = track(tracker, generic_solutions[1])\n",
    "r2 = track(tracker, generic_solutions[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Lets check that the tracking was successfull for both\n",
    "r1.returncode == PathTrackerStatus.success && r2.returncode == PathTrackerStatus.success"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we get our two specific solutions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Complex{Float64},1}:\n",
       "  -0.4432752095469262 - 1.2672493700951796e-19im\n",
       " -0.18955251455468036 + 9.362841922611108e-19im \n",
       "  0.20328610096481764 - 7.293918355685165e-20im "
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r1.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Complex{Float64},1}:\n",
       " 0.45410326818290414 + 9.565951274599459e-24im \n",
       "  1.3871995772087389 - 3.60671834508408e-23im  \n",
       "  0.1320653938006548 + 1.3127235321293564e-24im"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r2.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also see that these are complex vectors, but we are interested in the real solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First solution: [-0.443275, -0.189553, 0.203286]\n",
      "Second solution: [0.454103, 1.3872, 0.132065]\n"
     ]
    }
   ],
   "source": [
    "if maximum(abs ∘ imag, r1.x) < 1e-8\n",
    "    println(\"First solution: \", real.(r1.x))\n",
    "end\n",
    "if maximum(abs ∘ imag, r2.x) < 1e-8\n",
    "    println(\"Second solution: \", real.(r2.x))\n",
    "end"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.1.0",
   "language": "julia",
   "name": "julia-1.1"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}