added find_dark2 (improved find_dark method)

doc/find_dark.ipynb

@@ -276,8 +276,8 @@
     "def estimate_dark_lens_fluo(I, A, P, R, f_1, F_1):\n",
     "    lenses = compute_lenses(A, P, f_1)\n",
     "    fluos = estimate_fluos(f_1, F_1, R)\n",
-    "    lens = lenses[0]\n",
+    "    lens = np.mean(lenses, axis=0)\n",
-    "    fluo = fluos[:,0]\n",
+    "    fluo = np.mean(fluos, axis=1)\n",
     "    dark = I - lens.reshape(1, len(lens)) * fluo.reshape(len(fluo), 1)\n",
     "    return dark, lens, fluo\n",
     "\n",
@@ -326,8 +326,9 @@
     "        dark_matrix, lens, fluo = estimate_dark_lens_fluo(I, A, P, R, f_1, F_1)\n",
     "        error = np.var(dark_matrix)\n",
     "        return error\n",
-    "    \n",
+    "    start_f_1 = -0.7\n",
-    "    minimize_result = optimize.minimize(fun=function_to_minimize, x0=(-0.7, 0.5), method = 'Nelder-Mead')\n",
+    "    start_F_1 = 0.5\n",
+    "    minimize_result = optimize.minimize(fun=function_to_minimize, x0=(start_f_1, start_F_1), method = 'Nelder-Mead')\n",
     "    if not minimize_result.success:\n",
     "        print(type(minimize_result))\n",
     "        print(minimize_result)\n",
@@ -370,7 +371,111 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Even if the values of $f_1$ and $F_1$ are not correct, the exact value of for the dark $D$ has been found, which proves that this method is able to estimate the dark. We can also observe that $F$ and $L$ signals are recovered up to a scale factor."
+    "Even if the values of $f_1$ and $F_1$ are not correct, the exact value of for the dark $D$ has been found, which proves that this method is able to estimate the dark. We can also observe that $F$ and $L$ signals are recovered up to a scale factor:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lens_ratio = estimated_lens / lens\n",
+    "fluo_ratio = estimated_fluo / fluo\n",
+    "print('lens_ratio:',lens_ratio)\n",
+    "print('fluo_ratio:',fluo_ratio)\n",
+    "print('lens_ratio*fluo_ratio:', np.mean(lens_ratio)*np.mean(fluo_ratio))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## adding the constraint `F_1=fluo[0]`\n",
+    "\n",
+    "When we think about it, it makes sense $F$ and $L$ signals are recovered up to a scale factor only; indeed scaling $L$ up a factor $r$ has the same effect on $I$ as scaling $F$ up a factor $r$. In other words, the couple `(lens, fluo)` produce the same image as the couple `(estimated_lens, estimated_fluo)`. As a result, the sclaing factor $r$ is arbitrary. So, we need to reduce the degree of freedom. As a convention we then decide to fix the value of $F_1=1$. The corresponding code becomes:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def find_dark2(I, F_1=1.0):\n",
+    "    (num_times, num_sites) = I.shape\n",
+    "    A = I[1:num_times] - I[0:num_times-1]\n",
+    "    R = A[1:num_times-1] / A[0:num_times-2]\n",
+    "    P = np.multiply.accumulate(R)\n",
+    "    \n",
+    "    def function_to_minimize(x):\n",
+    "        f_1 = x[0]\n",
+    "        dark_matrix, lens, fluo = estimate_dark_lens_fluo(I, A, P, R, f_1, F_1)\n",
+    "        error = np.var(dark_matrix)\n",
+    "        return error\n",
+    "    start_f_1 = -0.7\n",
+    "    minimize_result = optimize.minimize(fun=function_to_minimize, x0=(start_f_1), method = 'Nelder-Mead')\n",
+    "    if not minimize_result.success:\n",
+    "        print(type(minimize_result))\n",
+    "        print(minimize_result)\n",
+    "        raise Exception(minimize_result.message)\n",
+    "    f_1 = minimize_result.x[0]\n",
+    "    print('estimated value of f_1 : %f' % f_1)\n",
+    "    dark, lens, fluo = estimate_dark_lens_fluo(I, A, P, R, f_1, F_1)\n",
+    "    return dark, lens, fluo\n",
+    "\n",
+    "print('exact value of f_1 : %f' % (fluo[1]-fluo[0]))\n",
+    "\n",
+    "estimated_dark, estimated_lens, estimated_fluo = find_dark2(i, fluo[0])\n",
+    "print('estimated dark value : ', estimated_dark)\n",
+    "print('exact lens : ', lens)\n",
+    "print('estimated lens : ', estimated_lens)\n",
+    "pyplot.plot(lens)\n",
+    "pyplot.plot(estimated_lens)\n",
+    "pyplot.xlabel('$s$ (pixel position)')\n",
+    "pyplot.ylabel('$L(s)$')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print('exact fluo :', fluo)\n",
+    "print('estimated fluo', estimated_fluo)\n",
+    "pyplot.plot(fluo)\n",
+    "pyplot.plot(estimated_fluo)\n",
+    "pyplot.xlabel('$t$ (time)')\n",
+    "pyplot.ylabel('$F(s)$')\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now that `F_1` is fixed to the exact value `fluo[0]`, both $L$ and $F$ estimations are the same as the ones that were used to construct $I$.\n",
+    "\n",
+    "Let's see how `find_dark2` is tolerant to noise:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "estimated_dark, estimated_lens, estimated_fluo = find_dark2(i + np.random.rand(i.shape[0], i.shape[1]) * 0.1, F_1=fluo[0])\n",
+    "print('estimated dark:', np.mean(estimated_dark))\n",
+    "print('exact dark:', real_dark)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## conclusion\n",
+    "`find_dark2` function is able to estimate $D$, $F$ and $L$ from $I$, which is the goal of this study."
    ]
   },
   {