colab/flag_test.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "flag_test.ipynb",
      "provenance": [],
      "authorship_tag": "ABX9TyPsE0cUdNX5zFhv8yXu5XH/",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/dnguyend/ManNullRange/blob/master/colab/flag_test.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TTNItWhjbl1x",
        "colab_type": "text"
      },
      "source": [
        "$\\newcommand{\\KK}{\\mathbb{K}}$\n",
        "$\\newcommand{\\rN}{\\mathrm{N}}$\n",
        "$\\newcommand{\\rD}{\\mathrm{D}}$\n",
        "$\\newcommand{\\rK}{\\mathrm{K}}$\n",
        "$\\newcommand{\\cE}{\\mathcal{E}}$\n",
        "$\\newcommand{\\cB}{\\mathcal{B}}$\n",
        "$\\newcommand{\\cH}{\\mathcal{H}}$\n",
        "$\\newcommand{\\ft}{\\mathfrak{t}}$\n",
        "$\\newcommand{\\cEJ}{\\cE_{\\JJ}}$\n",
        "$\\newcommand{\\sfg}{\\mathtt{g}}$\n",
        "$\\newcommand{\\rhess}{\\mathtt{rhess}}$\n",
        "$\\newcommand{\\CTRL}{\\textsc{CTRL}}$\n",
        "$\\newcommand{\\Herm}[3]{\\mathrm{Sym}_{#1, #2, #3}}$\n",
        "$\\newcommand{\\AHerm}[3]{\\mathrm{Skew}_{#1, #2, #3}}$\n",
        "$\\newcommand{\\St}[3]{\\mathrm{St}_{#1, #2, #3}}$\n",
        "$\\newcommand{\\Sp}[3]{\\mathtt{S}^{+}_{#1, #2, #3}}$\n",
        "$\\newcommand{\\Sd}[2]{\\mathtt{S}^{+}_{#1, #2}}$\n",
        "$\\DeclareMathOperator{\\UU}{U}$\n",
        "$\\DeclareMathOperator{\\sym}{sym}$\n",
        "$\\DeclareMathOperator{\\asym}{skew}$\n",
        "$\\DeclareMathOperator{\\JJ}{J}$\n",
        "$\\DeclareMathOperator{\\Null}{Null}$\n",
        "$\\DeclareMathOperator{\\xtrace}{xtrace}$\n",
        "# Formulas used here:\n",
        "* $\\JJ_{rr}\\eta =\\alpha_{rr}Y_r^{\\ft}\\eta_r$\n",
        "* $\\JJ_{rs} = Y_r^{\\ft}\\eta_s + \\eta_r^{\\ft}Y_s$\n",
        "* $(\\sfg(Y)\\omega)_r = \\alpha_{r0}\\omega_r +\\sum_{j=1}^p (\\alpha_{rj}-\\alpha_{r0})Y_jY_j^{\\ft}\\omega_r\n",
        "$ a metric given as self-adjoint operator on ambient $\\cE$\n",
        "* $\\Pi_{\\sfg} = I - \\sfg^{-1}\\JJ^{\\ft}(\\JJ\\sfg^{-1}\\JJ^{\\ft})^{-1}\\JJ$ projection to $\\cH$, the horizontal space, identified with $T\\cB = T\\Sp{\\KK}{p}{n}$\n",
        "* Gradient is $\\Pi_{\\sfg}\\sfg^{-1}f_Y$ \n",
        "* $\\xtrace$: index raising operator for trace (Frobenius) inner product. Very simple for matrix expressions:\n",
        "   * $\\xtrace(AbC, b) = A^{\\ft}B^{\\ft}$ \n",
        "   * $\\xtrace(Ab^tC, b) = BA$ \n",
        "* $\\JJ^{\\ft}$ is evaluated by $\\xtrace$\n",
        "* $\\rK(\\xi, \\eta)$  Christoffel metric term  $\\frac{1}{2}((\\rD_{\\xi}\\sfg)\\eta + (\\rD_{\\eta})\\sfg\\xi-\\xtrace(\\langle(\\rD_\\phi\\sfg)\\xi, \\eta\\rangle_{\\cE}, \\phi))$. \\\\\n",
        "* $\\Gamma_{c}(\\xi, \\eta)$  Christoffel function $\\Pi_{\\sfg}\\sfg^{-1}\\rK(\\xi, \\eta)-(\\rD_{\\xi}\\Pi_{\\cH, \\sfg})\\eta$.\n",
        "* $\\nabla_{\\xi}\\eta = \\rD_{\\xi}\\eta + \\Gamma_{c}(\\xi, \\eta)$\n",
        "* $\\rhess^{02}(\\xi, \\eta) = f_{YY}-\\Gamma_c(\\xi, \\eta)$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4w49TJdo4joe",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 136
        },
        "outputId": "6ea45b1f-ac54-4173-e920-e8547e2f2505"
      },
      "source": [
        "!git clone https://github.com/dnguyend/ManNullRange.git"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Cloning into 'ManNullRange'...\n",
            "remote: Enumerating objects: 144, done.\u001b[K\n",
            "remote: Counting objects: 100% (144/144), done.\u001b[K\n",
            "remote: Compressing objects: 100% (89/89), done.\u001b[K\n",
            "remote: Total 144 (delta 95), reused 89 (delta 54), pack-reused 0\u001b[K\n",
            "Receiving objects: 100% (144/144), 165.37 KiB | 5.17 MiB/s, done.\n",
            "Resolving deltas: 100% (95/95), done.\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WLnBVbKM7R1b",
        "colab_type": "text"
      },
      "source": [
        "We adopt our framework for block matrices in this case, list_xtrace, list_DDR. list_spfy are block versions of xtrace, etc, for example x_trace is:\n",
        "Signature: list_xtrace(expr, v)\n",
        "Source:   \n",
        "def list_xtrace(expr, v):\n",
        "    return [xtrace(expr, vv) for vv in v]\n",
        "File:      /content/ManNullRange/symbolic/flag.py\n",
        "Type:      function\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "O7Qq4SAi5Gl5",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from collections import OrderedDict\n",
        "from IPython.display import display, Math\n",
        "from sympy import Integer\n",
        "from ManNullRange.symbolic.SymMat import (\n",
        "    latex_map, mat_latex, t, trace, mat_spfy, xtrace)\n",
        "from ManNullRange.symbolic.flag import (\n",
        "    make_stiefel_block, make_ambient_block, make_b, make_al_stiefel, make_al,\n",
        "    list_spfy, list_DDR, list_xtrace)\n",
        "\n",
        "\"\"\"\n",
        "from ManNullRange.symbolic import SymMat as sm\n",
        "from ManNullRange.symbolic.SymMat import (\n",
        "    matrices, t, scalars, mat_spfy, xtrace, trace, stiefels, DDR,\n",
        "    LatexMap, MatLatex)\n",
        "\"\"\"\n",
        "\n",
        "def pprint(expr):\n",
        "    display(Math(latex_map(mat_latex(expr), OrderedDict(\n",
        "        [('fYY', r'(f_YY})'), ('fY', '(f_Y)'), ('al', r'\\alpha')]))))\n",
        "    \n",
        "def get_al_idx(r, j):\n",
        "    return (r-1)*(p+1)+j\n",
        "\n",
        "def make_alpha_w_beta(alpha0, beta):\n",
        "    al = []\n",
        "    p = len(alpha0)    \n",
        "    for r in range(1, p+1):\n",
        "        aidx = get_al_idx(r, 0)\n",
        "        al.append(alpha0[r-1])\n",
        "        for s in range(1, p+1):\n",
        "          aidx = get_al_idx(r, s)\n",
        "          al.append((alpha0[r-1] + alpha0[s-1] + beta)/2)\n",
        "    return al\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "25p8D2iDzlET",
        "colab_type": "text"
      },
      "source": [
        "We define alpha, etc. The expression with symbolic alpha turns out to be quite long (we can turn it on by using al = make_al(p)), so we use the numerical alpha, beta. $Y$ is a manifold point, $eta$ is an ambient vector, $b$ is a vector in range of $J$:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1Buv3CEruHJQ",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 284
        },
        "outputId": "62f20014-d4b8-44ca-9f69-2665ea26c281"
      },
      "source": [
        "if True:\n",
        "    \"\"\"block Stiefel matrix with 3 blocks\n",
        "    \"\"\"\n",
        "    p = 3\n",
        "    Y = make_stiefel_block(p, 'Y')\n",
        "    eta = make_ambient_block(p, 'eta')\n",
        "    b = make_b(p)\n",
        "    # al = make_al(p)\n",
        "    # al = make_al_stiefel(p, [Integer(5), Integer(1)])    \n",
        "    al = make_alpha_w_beta(p*[Integer(1)], Integer(-1))\n",
        "    print(\"Alpha:\")\n",
        "    pprint(al)\n",
        "    # al = range(1, 1+p*(p+1)//2)\n",
        "\n",
        "    def get_al_idx(r, j):\n",
        "        return (r-1)*(p+1)+j\n",
        "    \n",
        "    def J(Y, eta):\n",
        "        a = []\n",
        "        for r in range(p):\n",
        "            aidx = get_al_idx(r+1, r+1)\n",
        "            a.append(al[aidx]*t(Y[r])*eta[r])\n",
        "            for s in range(r+1, p):\n",
        "                a.append(t(Y[r])*eta[s] + t(eta[r])*Y[s])\n",
        "        return a\n",
        "\n",
        "    def inner_ambient(v1, v2):\n",
        "        expr = Integer(0)\n",
        "        for j in range(p):\n",
        "            expr += trace(v1[j]*t(v2[j]))\n",
        "        return expr\n",
        "\n",
        "    def inner_range_J(a1, a2):\n",
        "        expr = Integer(0)\n",
        "        for idx in range(len(a1)):\n",
        "            expr += trace(a1[idx]*t(a2[idx]))\n",
        "        return expr\n",
        "    \n",
        "    def J_adj(Y, a):\n",
        "        dY = make_ambient_block(p, 'dY')\n",
        "        jaj = []\n",
        "        jsim = [mat_spfy(jj) for jj in J(Y, dY)]\n",
        "        expr = mat_spfy(inner_range_J(jsim, a))\n",
        "        for r in range(p):\n",
        "            jaj.append(xtrace(expr, dY[r]))\n",
        "        return jaj\n",
        "\n",
        "    def g(Y, eta):\n",
        "        expr = []\n",
        "        for r in range(1, p+1):\n",
        "            zidx = get_al_idx(r, 0)\n",
        "            expr.append(al[zidx]*eta[r-1])\n",
        "            for j in range(1, p+1):\n",
        "                rjidx = get_al_idx(r, j)\n",
        "                expr[r-1] += (al[rjidx]-al[zidx])*Y[j-1]*t(Y[j-1])*eta[r-1]\n",
        "        return expr\n",
        "\n",
        "    def g_inv(Y, eta):\n",
        "        expr = []\n",
        "        for r in range(1, p+1):\n",
        "            zidx = get_al_idx(r, 0)\n",
        "            expr.append(1/al[zidx]*eta[r-1])\n",
        "            for j in range(1, p+1):\n",
        "                rjidx = get_al_idx(r, j)\n",
        "                expr[r-1] += (1/al[rjidx]-1/al[zidx])*Y[j-1]*t(Y[j-1])*eta[r-1]\n",
        "        return expr\n",
        "    print(\"this is Y\")        \n",
        "    pprint(Y)\n",
        "    print(\"this is eta\")\n",
        "    pprint(eta)\n",
        "    print(\"this is b\")\n",
        "    pprint(b)\n",
        "\n",
        "    j1 = J(Y, eta)\n",
        "    j2 = J_adj(Y, b)\n",
        "    pprint(\"J\\\\text{ and }J^T:\")\n",
        "    pprint(j1)\n",
        "    pprint(j2)\n",
        "    J_giv_J_adj = list_spfy(J(Y, g_inv(Y, J_adj(Y, b))), Y, [])\n",
        "    pprint(\"Jg^{-1}J^T:\")\n",
        "    pprint(J_giv_J_adj)\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Alpha:\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "text/latex": "$$\\left [ 1, \\quad \\frac{1}{2}, \\quad \\frac{1}{2}, \\quad \\frac{1}{2}, \\quad 1, \\quad \\frac{1}{2}, \\quad \\frac{1}{2}, \\quad \\frac{1}{2}, \\quad 1, \\quad \\frac{1}{2}, \\quad \\frac{1}{2}, \\quad \\frac{1}{2}\\right ]$$",
            "text/plain": [
              "<IPython.core.display.Math object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "stream",
          "text": [
            "this is Y\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "text/latex": "$$\\left ( Y_{1}, \\quad Y_{2}, \\quad Y_{3}\\right )$$",
            "text/plain": [
              "<IPython.core.display.Math object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "stream",
          "text": [
            "this is eta\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "text/latex": "$$\\left ( \\eta_{1}, \\quad \\eta_{2}, \\quad \\eta_{3}\\right )$$",
            "text/plain": [
              "<IPython.core.display.Math object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "stream",
          "text": [
            "this is b\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "text/latex": "$$\\left [ b_{11}, \\quad b_{12}, \\quad b_{13}, \\quad b_{22}, \\quad b_{23}, \\quad b_{33}\\right ]$$",
            "text/plain": [
              "<IPython.core.display.Math object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "text/latex": "$$J\\text{ and }J^T:$$",
            "text/plain": [
              "<IPython.core.display.Math object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "text/latex": "$$\\left [ \\frac{Y_{1}^T \\eta_{1}}{2}, \\quad Y_{1}^T \\eta_{2} + \\eta_{1}^T Y_{2}, \\quad Y_{1}^T \\eta_{3} + \\eta_{1}^T Y_{3}, \\quad \\frac{Y_{2}^T \\eta_{2}}{2}, \\quad Y_{2}^T \\eta_{3} + \\eta_{2}^T Y_{3}, \\quad \\frac{Y_{3}^T \\eta_{3}}{2}\\right ]$$",
            "text/plain": [
              "<IPython.core.display.Math object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "text/latex": "$$\\left [ \\frac{Y_{1} b_{11}}{2} + Y_{2} b_{12}^T + Y_{3} b_{13}^T, \\quad Y_{1} b_{12} + \\frac{Y_{2} b_{22}}{2} + Y_{3} b_{23}^T, \\quad Y_{1} b_{13} + Y_{2} b_{23} + \\frac{Y_{3} b_{33}}{2}\\right ]$$",
            "text/plain": [
              "<IPython.core.display.Math object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "text/latex": "$$Jg^{-1}J^T:$$",
            "text/plain": [
              "<IPython.core.display.Math object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "text/latex": "$$\\left [ \\frac{b_{11}}{2}, \\quad 4 b_{12}, \\quad 4 b_{13}, \\quad \\frac{b_{22}}{2}, \\quad 4 b_{23}, \\quad \\frac{b_{33}}{2}\\right ]$$",
            "text/plain": [
              "<IPython.core.display.Math object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ei9NsbTC9FuQ",
        "colab_type": "text"
      },
      "source": [
        "From here, we have $Jg^{-1}J^T$ is diagonal, and we can take its inverse:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UqWTEBIq7ldJ",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 124
        },
        "outputId": "5e08a871-29e4-4c14-fcc3-88b1975ba0d9"
      },
      "source": [
        "def solve_JginvJadj(Y, b):\n",
        "    ret = []\n",
        "    cnt = 0\n",
        "    for r in range(1, p+1):\n",
        "        for s in range(r, p+1):\n",
        "            i_r_s = get_al_idx(r, s)\n",
        "\n",
        "            if r == s:\n",
        "                ret.append(1/al[i_r_s]*b[cnt])\n",
        "            else:\n",
        "                i_s_r = get_al_idx(s, r)\n",
        "                ret.append((al[i_r_s]*al[i_s_r]/(al[i_r_s]+al[i_s_r])) *\n",
        "                            b[cnt])\n",
        "            cnt += 1\n",
        "    return ret\n",
        "\n",
        "pprint(solve_JginvJadj(Y, J_giv_J_adj))\n",
        "\n",
        "def g_inv_J_adj(Y, a):\n",
        "    expr = g_inv(Y, list_spfy(J_adj(Y, a), Y, []))\n",
        "    return list_spfy(expr, Y, [])\n",
        "\n",
        "g_inv_J_adj_expr = g_inv_J_adj(Y, b)\n",
        "\n",
        "def g_inv_J_adj_func(Y, a):\n",
        "    rdict = dict((b[i], a[i]) for i in range(len(b)))\n",
        "    return [list_spfy(g_inv_J_adj_expr[i].xreplace(rdict), Y, [])\n",
        "            for i in range(len(a))]\n",
        "\n",
        "def proj(Y, omg):\n",
        "    w = g_inv_J_adj(Y, solve_JginvJadj(Y, list_spfy(J(Y, omg), Y, [])))\n",
        "    return [mat_spfy(omg[i] - w[i]) for i in range(len(omg))]\n",
        "proj_expr = list_spfy(proj(Y, eta), Y, [])\n",
        "\n",
        "def proj_func(Y, omg):\n",
        "    rdict = dict((eta[i], omg[i]) for i in range(len(omg)))\n",
        "    return [proj_expr[i].xreplace(rdict) for i in range(len(omg))]\n",
        "\n",
        "r_gradient_expr = list_spfy(\n",
        "    proj_func(Y, g_inv(Y, eta)), Y, [])\n",
        "                            \n",
        "def r_gradient(Y, omg):\n",
        "    rdict = dict((eta[i], omg[i]) for i in range(len(omg)))\n",
        "    return [r_gradient_expr[i].xreplace(rdict)\n",
        "            for i in range(len(omg))]\n",
        "print(\"This is the projection:\")\n",
        "pprint(proj_expr)\n",
        "print(\"This is the Riemannian gradient:\")\n",
        "pprint(r_gradient_expr)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/latex": "$$\\left [ b_{11}, \\quad b_{12}, \\quad b_{13}, \\quad b_{22}, \\quad b_{23}, \\quad b_{33}\\right ]$$",
            "text/plain": [
              "<IPython.core.display.Math object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "stream",
          "text": [
            "This is the projection:\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "text/latex": "$$\\left [ - Y_{1} Y_{1}^T \\eta_{1} - \\frac{Y_{2} Y_{2}^T}{2} \\eta_{1} - \\frac{Y_{2} \\eta_{2}^T}{2} Y_{1} - \\frac{Y_{3} Y_{3}^T}{2} \\eta_{1} - \\frac{Y_{3} \\eta_{3}^T}{2} Y_{1} + \\eta_{1}, \\quad - \\frac{Y_{1} Y_{1}^T}{2} \\eta_{2} - \\frac{Y_{1} \\eta_{1}^T}{2} Y_{2} - Y_{2} Y_{2}^T \\eta_{2} - \\frac{Y_{3} Y_{3}^T}{2} \\eta_{2} - \\frac{Y_{3} \\eta_{3}^T}{2} Y_{2} + \\eta_{2}, \\quad - \\frac{Y_{1} Y_{1}^T}{2} \\eta_{3} - \\frac{Y_{1} \\eta_{1}^T}{2} Y_{3} - \\frac{Y_{2} Y_{2}^T}{2} \\eta_{3} - \\frac{Y_{2} \\eta_{2}^T}{2} Y_{3} - Y_{3} Y_{3}^T \\eta_{3} + \\eta_{3}\\right ]$$",
            "text/plain": [
              "<IPython.core.display.Math object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "stream",
          "text": [
            "This is the Riemannian gradient:\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "text/latex": "$$\\left [ - Y_{1} Y_{1}^T \\eta_{1} - Y_{2} \\eta_{2}^T Y_{1} - Y_{3} \\eta_{3}^T Y_{1} + \\eta_{1}, \\quad - Y_{1} \\eta_{1}^T Y_{2} - Y_{2} Y_{2}^T \\eta_{2} - Y_{3} \\eta_{3}^T Y_{2} + \\eta_{2}, \\quad - Y_{1} \\eta_{1}^T Y_{3} - Y_{2} \\eta_{2}^T Y_{3} - Y_{3} Y_{3}^T \\eta_{3} + \\eta_{3}\\right ]$$",
            "text/plain": [
              "<IPython.core.display.Math object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9KXWXzy59RoQ",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 181
        },
        "outputId": "6c78eb90-6266-4c82-9109-fac077332b03"
      },
      "source": [
        "    xi = make_ambient_block(p, 'xi')\n",
        "    phi = make_ambient_block(p, 'phi')\n",
        "\n",
        "    # jo = list_spfy(g(Y, eta), Y, [])\n",
        "    # ee = list_DDR(jo, Y, phi)\n",
        "    \n",
        "    trilinear = mat_spfy(inner_ambient(list_DDR(g(Y, eta), Y, phi), xi))\n",
        "    xcross = list_xtrace(trilinear, phi)\n",
        "    \n",
        "    # right now code return none if there is no term. make it zero b/c when beta=0 we have constant metric.    \n",
        "    def replace_none(itm):\n",
        "      if itm is None:\n",
        "        return 0\n",
        "      else: return itm\n",
        "    xcross = [replace_none(itm) for itm in xcross]\n",
        "\n",
        "    ddr_xi = list_DDR(g(Y, eta), Y, xi)\n",
        "    ddr_eta = list_DDR(g(Y, xi), Y, eta)\n",
        "    \n",
        "    K1 = [(Integer(1)/Integer(2))*(ddr_eta[jj] + ddr_xi[jj] - xcross[jj])\n",
        "          for jj in range(len(xcross))]\n",
        "    K = list_spfy(K1, Y, [xi, eta])\n",
        "    print(\"This is the Christoffel metric term:\")\n",
        "    pprint(K)\n",
        "\n",
        "    d_proj_expr = list_DDR(proj_expr, Y, xi)\n",
        "    \n",
        "    def d_proj_func(Y, xi, omg):\n",
        "        rdict = dict((eta[i], omg[i]) for i in range(len(omg)))\n",
        "        return [aa.xreplace(rdict) for aa in d_proj_expr]\n",
        "\n",
        "    dp_xi_eta = d_proj_func(Y, xi, eta)\n",
        "    prK0 = r_gradient(Y, K)\n",
        "    prK1 = list_spfy(prK0, Y, [])\n",
        "    prK = list_spfy(prK1, Y, (xi, eta))\n",
        "    Gamma0 = [prK[jj] - dp_xi_eta[jj] for jj in range(len(prK))]\n",
        "    Gamma = mat_spfy(\n",
        "        list_spfy(Gamma0, Y, (xi, eta)))\n",
        "    # one more time\n",
        "    Gamma = mat_spfy(\n",
        "        list_spfy(Gamma, Y, (xi, eta)))\n",
        "    \n",
        "    print(\"This is the Christoffel function:\")\n",
        "    pprint(Gamma)\n",
        "    fY = make_ambient_block(p, 'fY')\n",
        "\n",
        "    # t(eta)*fYY*xi\n",
        "    rhess02 = mat_spfy(inner_ambient([-Gamma[jj] for jj in range(p)], fY))\n",
        "    rhess11_bf_gr = list_xtrace(rhess02, eta)\n",
        "    print(\"This is the Riemannian Hessian Vector Product:\")\n",
        "    pprint(rhess11_bf_gr)\n",
        "    \n",
        "\n",
        "    \n",
        "\n",
        "            \n",
        "    \n",
        "\n",
        "    "
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "This is the Christoffel metric term:\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "text/latex": "$$\\left [ - \\frac{Y_{1} \\eta_{1}^T}{4} \\xi_{1} - \\frac{Y_{1} \\xi_{1}^T}{4} \\eta_{1} - \\frac{Y_{2} \\eta_{2}^T}{4} \\xi_{1} - \\frac{Y_{2} \\xi_{2}^T}{4} \\eta_{1} - \\frac{Y_{3} \\eta_{3}^T}{4} \\xi_{1} - \\frac{Y_{3} \\xi_{3}^T}{4} \\eta_{1} - \\frac{\\eta_{2} Y_{2}^T}{2} \\xi_{1} - \\frac{\\eta_{3} Y_{3}^T}{2} \\xi_{1} - \\frac{\\xi_{2} Y_{2}^T}{2} \\eta_{1} - \\frac{\\xi_{3} Y_{3}^T}{2} \\eta_{1}, \\quad - \\frac{Y_{1} \\eta_{1}^T}{4} \\xi_{2} - \\frac{Y_{1} \\xi_{1}^T}{4} \\eta_{2} - \\frac{Y_{2} \\eta_{2}^T}{4} \\xi_{2} - \\frac{Y_{2} \\xi_{2}^T}{4} \\eta_{2} - \\frac{Y_{3} \\eta_{3}^T}{4} \\xi_{2} - \\frac{Y_{3} \\xi_{3}^T}{4} \\eta_{2} - \\frac{\\eta_{1} Y_{1}^T}{2} \\xi_{2} - \\frac{\\eta_{3} Y_{3}^T}{2} \\xi_{2} - \\frac{\\xi_{1} Y_{1}^T}{2} \\eta_{2} - \\frac{\\xi_{3} Y_{3}^T}{2} \\eta_{2}, \\quad - \\frac{Y_{1} \\eta_{1}^T}{4} \\xi_{3} - \\frac{Y_{1} \\xi_{1}^T}{4} \\eta_{3} - \\frac{Y_{2} \\eta_{2}^T}{4} \\xi_{3} - \\frac{Y_{2} \\xi_{2}^T}{4} \\eta_{3} - \\frac{Y_{3} \\eta_{3}^T}{4} \\xi_{3} - \\frac{Y_{3} \\xi_{3}^T}{4} \\eta_{3} - \\frac{\\eta_{1} Y_{1}^T}{2} \\xi_{3} - \\frac{\\eta_{2} Y_{2}^T}{2} \\xi_{3} - \\frac{\\xi_{1} Y_{1}^T}{2} \\eta_{3} - \\frac{\\xi_{2} Y_{2}^T}{2} \\eta_{3}\\right ]$$",
            "text/plain": [
              "<IPython.core.display.Math object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "stream",
          "text": [
            "This is the Christoffel function:\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "text/latex": "$$\\left [ \\frac{Y_{1} Y_{1}^T}{2} \\eta_{2} Y_{2}^T \\xi_{1} + \\frac{Y_{1} Y_{1}^T}{2} \\eta_{3} Y_{3}^T \\xi_{1} + \\frac{Y_{1} Y_{1}^T}{2} \\xi_{2} Y_{2}^T \\eta_{1} + \\frac{Y_{1} Y_{1}^T}{2} \\xi_{3} Y_{3}^T \\eta_{1} + Y_{1} \\xi_{1}^T \\eta_{1} + \\frac{Y_{2} Y_{2}^T}{2} \\eta_{3} Y_{3}^T \\xi_{1} + \\frac{Y_{2} Y_{2}^T}{2} \\xi_{3} Y_{3}^T \\eta_{1} + \\frac{Y_{2} \\eta_{2}^T}{2} \\xi_{1} + \\frac{Y_{2} \\xi_{2}^T}{2} \\eta_{1} + \\frac{Y_{3} Y_{3}^T}{2} \\eta_{2} Y_{2}^T \\xi_{1} + \\frac{Y_{3} Y_{3}^T}{2} \\xi_{2} Y_{2}^T \\eta_{1} + \\frac{Y_{3} \\eta_{3}^T}{2} \\xi_{1} + \\frac{Y_{3} \\xi_{3}^T}{2} \\eta_{1} - \\frac{\\eta_{2} Y_{2}^T}{2} \\xi_{1} - \\frac{\\eta_{3} Y_{3}^T}{2} \\xi_{1} - \\frac{\\xi_{2} Y_{2}^T}{2} \\eta_{1} - \\frac{\\xi_{3} Y_{3}^T}{2} \\eta_{1}, \\quad \\frac{Y_{1} Y_{1}^T}{2} \\eta_{3} Y_{3}^T \\xi_{2} + \\frac{Y_{1} Y_{1}^T}{2} \\xi_{3} Y_{3}^T \\eta_{2} + \\frac{Y_{1} \\eta_{1}^T}{2} \\xi_{2} + \\frac{Y_{1} \\xi_{1}^T}{2} \\eta_{2} + \\frac{Y_{2} Y_{2}^T}{2} \\eta_{1} Y_{1}^T \\xi_{2} + \\frac{Y_{2} Y_{2}^T}{2} \\eta_{3} Y_{3}^T \\xi_{2} + \\frac{Y_{2} Y_{2}^T}{2} \\xi_{1} Y_{1}^T \\eta_{2} + \\frac{Y_{2} Y_{2}^T}{2} \\xi_{3} Y_{3}^T \\eta_{2} + Y_{2} \\xi_{2}^T \\eta_{2} + \\frac{Y_{3} Y_{3}^T}{2} \\eta_{1} Y_{1}^T \\xi_{2} + \\frac{Y_{3} Y_{3}^T}{2} \\xi_{1} Y_{1}^T \\eta_{2} + \\frac{Y_{3} \\eta_{3}^T}{2} \\xi_{2} + \\frac{Y_{3} \\xi_{3}^T}{2} \\eta_{2} - \\frac{\\eta_{1} Y_{1}^T}{2} \\xi_{2} - \\frac{\\eta_{3} Y_{3}^T}{2} \\xi_{2} - \\frac{\\xi_{1} Y_{1}^T}{2} \\eta_{2} - \\frac{\\xi_{3} Y_{3}^T}{2} \\eta_{2}, \\quad \\frac{Y_{1} Y_{1}^T}{2} \\eta_{2} Y_{2}^T \\xi_{3} + \\frac{Y_{1} Y_{1}^T}{2} \\xi_{2} Y_{2}^T \\eta_{3} + \\frac{Y_{1} \\eta_{1}^T}{2} \\xi_{3} + \\frac{Y_{1} \\xi_{1}^T}{2} \\eta_{3} + \\frac{Y_{2} Y_{2}^T}{2} \\eta_{1} Y_{1}^T \\xi_{3} + \\frac{Y_{2} Y_{2}^T}{2} \\xi_{1} Y_{1}^T \\eta_{3} + \\frac{Y_{2} \\eta_{2}^T}{2} \\xi_{3} + \\frac{Y_{2} \\xi_{2}^T}{2} \\eta_{3} + \\frac{Y_{3} Y_{3}^T}{2} \\eta_{1} Y_{1}^T \\xi_{3} + \\frac{Y_{3} Y_{3}^T}{2} \\eta_{2} Y_{2}^T \\xi_{3} + \\frac{Y_{3} Y_{3}^T}{2} \\xi_{1} Y_{1}^T \\eta_{3} + \\frac{Y_{3} Y_{3}^T}{2} \\xi_{2} Y_{2}^T \\eta_{3} + Y_{3} \\xi_{3}^T \\eta_{3} - \\frac{\\eta_{1} Y_{1}^T}{2} \\xi_{3} - \\frac{\\eta_{2} Y_{2}^T}{2} \\xi_{3} - \\frac{\\xi_{1} Y_{1}^T}{2} \\eta_{3} - \\frac{\\xi_{2} Y_{2}^T}{2} \\eta_{3}\\right ]$$",
            "text/plain": [
              "<IPython.core.display.Math object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "stream",
          "text": [
            "This is the Riemannian Hessian Vector Product:\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "text/latex": "$$\\left [ - \\frac{Y_{2} Y_{2}^T}{2} (f_Y)_{2} \\xi_{2}^T Y_{1} - \\frac{Y_{2} Y_{2}^T}{2} (f_Y)_{3} \\xi_{3}^T Y_{1} - \\frac{Y_{2} \\xi_{2}^T}{2} Y_{1} Y_{1}^T (f_Y)_{1} - \\frac{Y_{2} \\xi_{2}^T}{2} Y_{3} Y_{3}^T (f_Y)_{1} + \\frac{Y_{2} \\xi_{2}^T}{2} (f_Y)_{1} - \\frac{Y_{3} Y_{3}^T}{2} (f_Y)_{2} \\xi_{2}^T Y_{1} - \\frac{Y_{3} Y_{3}^T}{2} (f_Y)_{3} \\xi_{3}^T Y_{1} - \\frac{Y_{3} \\xi_{3}^T}{2} Y_{1} Y_{1}^T (f_Y)_{1} - \\frac{Y_{3} \\xi_{3}^T}{2} Y_{2} Y_{2}^T (f_Y)_{1} + \\frac{Y_{3} \\xi_{3}^T}{2} (f_Y)_{1} + \\frac{(f_Y)_{2} \\xi_{2}^T}{2} Y_{1} + \\frac{(f_Y)_{3} \\xi_{3}^T}{2} Y_{1} - \\xi_{1} Y_{1}^T (f_Y)_{1} - \\frac{\\xi_{2} Y_{2}^T}{2} (f_Y)_{1} - \\frac{\\xi_{2} (f_Y)_{2}^T}{2} Y_{1} - \\frac{\\xi_{3} Y_{3}^T}{2} (f_Y)_{1} - \\frac{\\xi_{3} (f_Y)_{3}^T}{2} Y_{1}, \\quad - \\frac{Y_{1} Y_{1}^T}{2} (f_Y)_{1} \\xi_{1}^T Y_{2} - \\frac{Y_{1} Y_{1}^T}{2} (f_Y)_{3} \\xi_{3}^T Y_{2} - \\frac{Y_{1} \\xi_{1}^T}{2} Y_{2} Y_{2}^T (f_Y)_{2} - \\frac{Y_{1} \\xi_{1}^T}{2} Y_{3} Y_{3}^T (f_Y)_{2} + \\frac{Y_{1} \\xi_{1}^T}{2} (f_Y)_{2} - \\frac{Y_{3} Y_{3}^T}{2} (f_Y)_{1} \\xi_{1}^T Y_{2} - \\frac{Y_{3} Y_{3}^T}{2} (f_Y)_{3} \\xi_{3}^T Y_{2} - \\frac{Y_{3} \\xi_{3}^T}{2} Y_{1} Y_{1}^T (f_Y)_{2} - \\frac{Y_{3} \\xi_{3}^T}{2} Y_{2} Y_{2}^T (f_Y)_{2} + \\frac{Y_{3} \\xi_{3}^T}{2} (f_Y)_{2} + \\frac{(f_Y)_{1} \\xi_{1}^T}{2} Y_{2} + \\frac{(f_Y)_{3} \\xi_{3}^T}{2} Y_{2} - \\frac{\\xi_{1} Y_{1}^T}{2} (f_Y)_{2} - \\frac{\\xi_{1} (f_Y)_{1}^T}{2} Y_{2} - \\xi_{2} Y_{2}^T (f_Y)_{2} - \\frac{\\xi_{3} Y_{3}^T}{2} (f_Y)_{2} - \\frac{\\xi_{3} (f_Y)_{3}^T}{2} Y_{2}, \\quad - \\frac{Y_{1} Y_{1}^T}{2} (f_Y)_{1} \\xi_{1}^T Y_{3} - \\frac{Y_{1} Y_{1}^T}{2} (f_Y)_{2} \\xi_{2}^T Y_{3} - \\frac{Y_{1} \\xi_{1}^T}{2} Y_{2} Y_{2}^T (f_Y)_{3} - \\frac{Y_{1} \\xi_{1}^T}{2} Y_{3} Y_{3}^T (f_Y)_{3} + \\frac{Y_{1} \\xi_{1}^T}{2} (f_Y)_{3} - \\frac{Y_{2} Y_{2}^T}{2} (f_Y)_{1} \\xi_{1}^T Y_{3} - \\frac{Y_{2} Y_{2}^T}{2} (f_Y)_{2} \\xi_{2}^T Y_{3} - \\frac{Y_{2} \\xi_{2}^T}{2} Y_{1} Y_{1}^T (f_Y)_{3} - \\frac{Y_{2} \\xi_{2}^T}{2} Y_{3} Y_{3}^T (f_Y)_{3} + \\frac{Y_{2} \\xi_{2}^T}{2} (f_Y)_{3} + \\frac{(f_Y)_{1} \\xi_{1}^T}{2} Y_{3} + \\frac{(f_Y)_{2} \\xi_{2}^T}{2} Y_{3} - \\frac{\\xi_{1} Y_{1}^T}{2} (f_Y)_{3} - \\frac{\\xi_{1} (f_Y)_{1}^T}{2} Y_{3} - \\frac{\\xi_{2} Y_{2}^T}{2} (f_Y)_{3} - \\frac{\\xi_{2} (f_Y)_{2}^T}{2} Y_{3} - \\xi_{3} Y_{3}^T (f_Y)_{3}\\right ]$$",
            "text/plain": [
              "<IPython.core.display.Math object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "p7zUUgIRjjI2",
        "colab_type": "text"
      },
      "source": [
        "Optimization tests:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1ieU7gh0lba9",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 136
        },
        "outputId": "e10d14cd-aa81-4b62-9a00-08bfe34557a6"
      },
      "source": [
        "!git clone  https://github.com/pymanopt/pymanopt.git\n",
        "import sys\n",
        "sys.path.append(\"/content/pymanopt\")"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Cloning into 'pymanopt'...\n",
            "remote: Enumerating objects: 77, done.\u001b[K\n",
            "remote: Counting objects: 100% (77/77), done.\u001b[K\n",
            "remote: Compressing objects: 100% (62/62), done.\u001b[K\n",
            "remote: Total 4124 (delta 43), reused 40 (delta 15), pack-reused 4047\u001b[K\n",
            "Receiving objects: 100% (4124/4124), 907.13 KiB | 17.12 MiB/s, done.\n",
            "Resolving deltas: 100% (2863/2863), done.\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YGkEbDm4jmrQ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "from numpy.random import (randint, randn)\n",
        "from numpy import zeros, zeros_like, trace, allclose\n",
        "\n",
        "from ManNullRange.manifolds.RealFlag import RealFlag\n",
        "from ManNullRange.tests.test_tools import check_zero, make_sym_pos, random_orthogonal"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GQpct6YVlqGH",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "outputId": "f5484758-a1d6-41a3-f0c3-78677c0bf901"
      },
      "source": [
        "from pymanopt import Problem\n",
        "from pymanopt.solvers import TrustRegions\n",
        "from pymanopt.function import Callable\n",
        "\n",
        "n = 1000\n",
        "\n",
        "def make_simple_alpha(alpha0, beta):\n",
        "    \"\"\"Make the alpha_{rl} dependent on alpha_{t0} and beta only\n",
        "    \"\"\"\n",
        "    p = alpha0.shape[0]\n",
        "    alpha = zeros((p, p+1))\n",
        "    alpha[:, 0] = alpha0\n",
        "    for i in range(p):\n",
        "        for j in range(i, p):\n",
        "            alpha[i, j+1] = alpha0[i]/2 + alpha0[j]/2 + beta/2\n",
        "            alpha[j, i+1] = alpha0[i]/2 + alpha0[j]/2 + beta/2\n",
        "    return alpha\n",
        "\n",
        "# problem Tr(AXBX^T)\n",
        "for i in range(1):\n",
        "    dvec = np.array([0, 30, 2, 1])\n",
        "    dvec[0] = 1000 - dvec[1:].sum()\n",
        "    d = dvec[1:].sum()\n",
        "    D = randint(1, 10, n) * 0.02 + 1\n",
        "    OO = random_orthogonal(n)\n",
        "    A = OO @ np.diag(D) @ OO.T\n",
        "    \n",
        "    B = make_sym_pos(d)\n",
        "    p = dvec.shape[0]-1\n",
        "    alpha = randint(1, 10, (p, p+1)) * .1\n",
        "    alpha0 = randint(1, 10, (p))\n",
        "    # alpha0 = randint(1, 2, (p))\n",
        "    alpha = make_simple_alpha(alpha0, 0)\n",
        "    man = RealFlag(dvec, alpha=alpha)\n",
        "\n",
        "    @Callable\n",
        "    def cost(X):\n",
        "        return trace(A @ X @ B @ X.T)\n",
        "    \n",
        "    @Callable\n",
        "    def egrad(X):\n",
        "        return 2*A @ X @ B\n",
        "\n",
        "    @Callable\n",
        "    def ehess(X, H):\n",
        "        return 2*A @ H @ B\n",
        "\n",
        "    if False:\n",
        "        X = man.rand()\n",
        "        xi = man.randvec(X)\n",
        "        d1 = num_deriv(man, X, xi, cost)\n",
        "        d2 = trace(egrad(X) @ xi.T)\n",
        "        print(check_zero(d1-d2))\n",
        "    \n",
        "    prob = Problem(\n",
        "        man, cost, egrad=egrad)\n",
        "    XInit = man.rand()\n",
        "\n",
        "    prob = Problem(\n",
        "        man, cost, egrad=egrad, ehess=ehess)\n",
        "\n",
        "    solver = TrustRegions(maxtime=100000, maxiter=100)\n",
        "    opt = solver.solve(prob, x=XInit, Delta_bar=250)\n",
        "    print(cost(opt))\n",
        "\n",
        "    if False:\n",
        "        min_val = 1e190\n",
        "        # min_X = None\n",
        "        for i in range(100):\n",
        "            Xi = man.rand()\n",
        "            c = cost(Xi)\n",
        "            if c < min_val:\n",
        "                # min_X = Xi\n",
        "                min_val = c\n",
        "            if i % 1000 == 0:\n",
        "                print('i=%d min=%f' % (i, min_val))\n",
        "        print(min_val)\n",
        "    alpha_c = alpha.copy()\n",
        "    alpha_c[:] = 1\n",
        "    man1 = RealFlag(dvec, alpha=alpha_c)\n",
        "    prob = Problem(\n",
        "        man1, cost, egrad=egrad, ehess=ehess)\n",
        "\n",
        "    solver = TrustRegions(maxtime=100000, maxiter=100)\n",
        "    opt = solver.solve(prob, x=XInit, Delta_bar=250)\n",
        "    alpha_c5 = alpha_c.copy()\n",
        "    alpha_c5[:, 1:] = .5\n",
        "    man1 = RealFlag(dvec, alpha=alpha_c5)\n",
        "    prob = Problem(\n",
        "        man1, cost, egrad=egrad, ehess=ehess)\n",
        "\n",
        "    solver = TrustRegions(maxtime=100000, maxiter=100)\n",
        "    opt = solver.solve(prob, x=XInit, Delta_bar=250)\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Optimizing...\n",
            "                                            f: +1.227771e+03   |grad|: 1.574115e+01\n",
            "REJ TR-   k:     1     num_inner:     0     f: +1.227771e+03   |grad|: 1.574115e+01   negative curvature\n",
            "acc       k:     2     num_inner:     0     f: +1.180410e+03   |grad|: 1.315424e+01   negative curvature\n",
            "acc       k:     3     num_inner:     1     f: +1.156380e+03   |grad|: 9.529445e+00   exceeded trust region\n",
            "acc       k:     4     num_inner:     3     f: +1.146140e+03   |grad|: 4.879151e+00   exceeded trust region\n",
            "acc       k:     5     num_inner:     8     f: +1.142528e+03   |grad|: 2.690995e+00   exceeded trust region\n",
            "REJ TR-   k:     6     num_inner:     7     f: +1.142528e+03   |grad|: 2.690995e+00   exceeded trust region\n",
            "acc TR+   k:     7     num_inner:     2     f: +1.140744e+03   |grad|: 1.398314e+00   exceeded trust region\n",
            "acc       k:     8     num_inner:    14     f: +1.140236e+03   |grad|: 5.397423e-01   exceeded trust region\n",
            "acc       k:     9     num_inner:    26     f: +1.140117e+03   |grad|: 3.201076e-01   exceeded trust region\n",
            "acc       k:    10     num_inner:    49     f: +1.140094e+03   |grad|: 2.958882e-01   exceeded trust region\n",
            "REJ TR-   k:    11     num_inner:    53     f: +1.140094e+03   |grad|: 2.958882e-01   exceeded trust region\n",
            "acc TR+   k:    12     num_inner:    22     f: +1.140062e+03   |grad|: 5.145325e-02   exceeded trust region\n",
            "acc       k:    13     num_inner:   100     f: +1.140059e+03   |grad|: 2.037071e-02   exceeded trust region\n",
            "acc       k:    14     num_inner:   161     f: +1.140058e+03   |grad|: 2.335143e-02   exceeded trust region\n",
            "REJ TR-   k:    15     num_inner:   274     f: +1.140058e+03   |grad|: 2.335143e-02   exceeded trust region\n",
            "acc TR+   k:    16     num_inner:   108     f: +1.140058e+03   |grad|: 1.564020e-03   exceeded trust region\n",
            "acc       k:    17     num_inner:   650     f: +1.140058e+03   |grad|: 1.371189e-04   model increased\n",
            "REJ TR-   k:    18     num_inner:   441     f: +1.140058e+03   |grad|: 1.371189e-04   negative curvature\n",
            "REJ TR-   k:    19     num_inner:   441     f: +1.140058e+03   |grad|: 1.371189e-04   negative curvature\n",
            "REJ TR-   k:    20     num_inner:   441     f: +1.140058e+03   |grad|: 1.371189e-04   negative curvature\n",
            "acc       k:    21     num_inner:   441     f: +1.140058e+03   |grad|: 4.591028e-05   negative curvature\n",
            "REJ TR-   k:    22     num_inner:   479     f: +1.140058e+03   |grad|: 4.591028e-05   negative curvature\n",
            "acc TR+   k:    23     num_inner:    92     f: +1.140058e+03   |grad|: 1.547327e-05   exceeded trust region\n",
            "REJ TR-   k:    24     num_inner:   382     f: +1.140058e+03   |grad|: 1.547327e-05   negative curvature\n",
            "acc TR+   k:    25     num_inner:   380     f: +1.140058e+03   |grad|: 6.206452e-06   exceeded trust region\n",
            "REJ TR-   k:    26     num_inner:   329     f: +1.140058e+03   |grad|: 6.206452e-06   exceeded trust region\n",
            "acc TR+   k:    27     num_inner:   320     f: +1.140058e+03   |grad|: 2.101010e-06   exceeded trust region\n",
            "REJ TR-   k:    28     num_inner:   320     f: +1.140058e+03   |grad|: 2.101010e-06   exceeded trust region\n",
            "acc TR+   k:    29     num_inner:   314     f: +1.140058e+03   |grad|: 1.838027e-06   exceeded trust region\n",
            "acc       k:    30     num_inner:   249     f: +1.140058e+03   |grad|: 3.395067e-06   exceeded trust region\n",
            "acc TR+   k:    31     num_inner:   346     f: +1.140058e+03   |grad|: 2.634301e-06   exceeded trust region\n",
            "REJ TR-   k:    32     num_inner:   348     f: +1.140058e+03   |grad|: 2.634301e-06   exceeded trust region\n",
            "acc TR+   k:    33     num_inner:   341     f: +1.140058e+03   |grad|: 1.091577e-06   exceeded trust region\n",
            "acc       k:    34     num_inner:   270     f: +1.140058e+03   |grad|: 4.228640e-06   exceeded trust region\n",
            "acc TR+   k:    35     num_inner:   332     f: +1.140058e+03   |grad|: 1.833959e-06   exceeded trust region\n",
            "REJ TR-   k:    36     num_inner:   325     f: +1.140058e+03   |grad|: 1.833959e-06   exceeded trust region\n",
            "acc TR+   k:    37     num_inner:   294     f: +1.140058e+03   |grad|: 5.091123e-07   exceeded trust region\n",
            "Terminated - min grad norm reached after 37 iterations, 776.40 seconds.\n",
            "\n",
            "1140.0577849753927\n",
            "Optimizing...\n",
            "                                            f: +1.227771e+03   |grad|: 2.756721e+01\n",
            "REJ TR-   k:     1     num_inner:     0     f: +1.227771e+03   |grad|: 2.756721e+01   exceeded trust region\n",
            "acc TR-   k:     2     num_inner:     0     f: +1.185200e+03   |grad|: 2.918256e+01   exceeded trust region\n",
            "acc       k:     3     num_inner:     0     f: +1.161577e+03   |grad|: 1.496005e+01   exceeded trust region\n",
            "acc TR+   k:     4     num_inner:     1     f: +1.149300e+03   |grad|: 8.383889e+00   exceeded trust region\n",
            "acc       k:     5     num_inner:     4     f: +1.143293e+03   |grad|: 4.250105e+00   exceeded trust region\n",
            "acc       k:     6     num_inner:    10     f: +1.141014e+03   |grad|: 2.438730e+00   exceeded trust region\n",
            "acc       k:     7     num_inner:    20     f: +1.140441e+03   |grad|: 1.545523e+00   exceeded trust region\n",
            "acc       k:     8     num_inner:    32     f: +1.140288e+03   |grad|: 1.214079e+00   exceeded trust region\n",
            "REJ TR-   k:     9     num_inner:    30     f: +1.140288e+03   |grad|: 1.214079e+00   exceeded trust region\n",
            "acc TR+   k:    10     num_inner:    12     f: +1.140098e+03   |grad|: 2.861011e-01   exceeded trust region\n",
            "acc TR-   k:    11     num_inner:    37     f: +1.140086e+03   |grad|: 2.873049e-01   exceeded trust region\n",
            "acc TR+   k:    12     num_inner:    15     f: +1.140063e+03   |grad|: 1.243429e-01   exceeded trust region\n",
            "acc       k:    13     num_inner:    84     f: +1.140059e+03   |grad|: 5.341399e-02   exceeded trust region\n",
            "acc       k:    14     num_inner:   136     f: +1.140058e+03   |grad|: 9.181523e-02   exceeded trust region\n",
            "acc       k:    15     num_inner:    40     f: +1.140058e+03   |grad|: 8.283581e-03   reached target residual-theta (superlinear)\n",
            "acc       k:    16     num_inner:   358     f: +1.140058e+03   |grad|: 1.471655e-04   reached target residual-theta (superlinear)\n",
            "acc       k:    17     num_inner:   465     f: +1.140058e+03   |grad|: 3.019542e-08   reached target residual-theta (superlinear)\n",
            "Terminated - min grad norm reached after 17 iterations, 116.02 seconds.\n",
            "\n",
            "Optimizing...\n",
            "                                            f: +1.227771e+03   |grad|: 2.759455e+01\n",
            "REJ TR-   k:     1     num_inner:     0     f: +1.227771e+03   |grad|: 2.759455e+01   exceeded trust region\n",
            "acc TR-   k:     2     num_inner:     0     f: +1.185274e+03   |grad|: 2.920060e+01   exceeded trust region\n",
            "acc       k:     3     num_inner:     0     f: +1.161602e+03   |grad|: 1.496308e+01   exceeded trust region\n",
            "acc TR+   k:     4     num_inner:     1     f: +1.149327e+03   |grad|: 8.392623e+00   exceeded trust region\n",
            "acc       k:     5     num_inner:     4     f: +1.143301e+03   |grad|: 4.256184e+00   exceeded trust region\n",
            "acc       k:     6     num_inner:     9     f: +1.141040e+03   |grad|: 2.558507e+00   exceeded trust region\n",
            "acc       k:     7     num_inner:    16     f: +1.140667e+03   |grad|: 2.016349e+00   exceeded trust region\n",
            "REJ TR-   k:     8     num_inner:    18     f: +1.140667e+03   |grad|: 2.016349e+00   exceeded trust region\n",
            "acc TR+   k:     9     num_inner:     6     f: +1.140162e+03   |grad|: 6.866492e-01   exceeded trust region\n",
            "acc       k:    10     num_inner:    35     f: +1.140093e+03   |grad|: 4.067748e-01   exceeded trust region\n",
            "REJ TR-   k:    11     num_inner:    38     f: +1.140093e+03   |grad|: 4.067748e-01   exceeded trust region\n",
            "acc TR+   k:    12     num_inner:    19     f: +1.140065e+03   |grad|: 7.768211e-02   exceeded trust region\n",
            "acc       k:    13     num_inner:    53     f: +1.140061e+03   |grad|: 7.900577e-02   exceeded trust region\n",
            "REJ TR-   k:    14     num_inner:    52     f: +1.140061e+03   |grad|: 7.900577e-02   exceeded trust region\n",
            "acc TR+   k:    15     num_inner:    31     f: +1.140059e+03   |grad|: 2.341646e-02   exceeded trust region\n",
            "acc TR+   k:    16     num_inner:    51     f: +1.140058e+03   |grad|: 1.275293e-02   exceeded trust region\n",
            "acc       k:    17     num_inner:   126     f: +1.140058e+03   |grad|: 3.626688e-02   exceeded trust region\n",
            "REJ TR-   k:    18     num_inner:    58     f: +1.140058e+03   |grad|: 3.626688e-02   exceeded trust region\n",
            "acc       k:    19     num_inner:    49     f: +1.140058e+03   |grad|: 9.940384e-03   exceeded trust region\n",
            "acc TR+   k:    20     num_inner:   235     f: +1.140058e+03   |grad|: 2.013772e-03   exceeded trust region\n",
            "REJ TR-   k:    21     num_inner:   364     f: +1.140058e+03   |grad|: 2.013772e-03   exceeded trust region\n",
            "acc TR+   k:    22     num_inner:   323     f: +1.140058e+03   |grad|: 5.556020e-04   exceeded trust region\n",
            "REJ TR-   k:    23     num_inner:   471     f: +1.140058e+03   |grad|: 5.556020e-04   exceeded trust region\n",
            "acc TR+   k:    24     num_inner:   468     f: +1.140058e+03   |grad|: 2.372564e-04   exceeded trust region\n",
            "REJ TR-   k:    25     num_inner:   441     f: +1.140058e+03   |grad|: 2.372564e-04   negative curvature\n",
            "REJ TR-   k:    26     num_inner:   441     f: +1.140058e+03   |grad|: 2.372564e-04   negative curvature\n",
            "acc TR+   k:    27     num_inner:   439     f: +1.140058e+03   |grad|: 5.048154e-05   exceeded trust region\n",
            "REJ TR-   k:    28     num_inner:   381     f: +1.140058e+03   |grad|: 5.048154e-05   negative curvature\n",
            "acc       k:    29     num_inner:   379     f: +1.140058e+03   |grad|: 4.741951e-05   exceeded trust region\n",
            "REJ TR-   k:    30     num_inner:   376     f: +1.140058e+03   |grad|: 4.741951e-05   exceeded trust region\n",
            "acc TR+   k:    31     num_inner:   114     f: +1.140058e+03   |grad|: 4.674058e-06   exceeded trust region\n",
            "REJ TR-   k:    32     num_inner:   259     f: +1.140058e+03   |grad|: 4.674058e-06   exceeded trust region\n",
            "acc TR+   k:    33     num_inner:   253     f: +1.140058e+03   |grad|: 5.644118e-06   exceeded trust region\n",
            "REJ TR-   k:    34     num_inner:   264     f: +1.140058e+03   |grad|: 5.644118e-06   exceeded trust region\n",
            "acc TR+   k:    35     num_inner:   234     f: +1.140058e+03   |grad|: 1.085803e-06   exceeded trust region\n",
            "acc       k:    36     num_inner:   169     f: +1.140058e+03   |grad|: 5.603100e-06   exceeded trust region\n",
            "acc TR+   k:    37     num_inner:   262     f: +1.140058e+03   |grad|: 3.890321e-06   exceeded trust region\n",
            "REJ TR-   k:    38     num_inner:   227     f: +1.140058e+03   |grad|: 3.890321e-06   exceeded trust region\n",
            "acc TR+   k:    39     num_inner:   217     f: +1.140058e+03   |grad|: 1.857496e-06   exceeded trust region\n",
            "acc       k:    40     num_inner:   215     f: +1.140058e+03   |grad|: 4.284278e-06   exceeded trust region\n",
            "acc TR+   k:    41     num_inner:   276     f: +1.140058e+03   |grad|: 4.856227e-06   exceeded trust region\n",
            "REJ TR-   k:    42     num_inner:   252     f: +1.140058e+03   |grad|: 4.856227e-06   exceeded trust region\n",
            "acc TR+   k:    43     num_inner:   229     f: +1.140058e+03   |grad|: 1.779719e-06   exceeded trust region\n",
            "acc       k:    44     num_inner:   202     f: +1.140058e+03   |grad|: 4.171663e-06   exceeded trust region\n",
            "acc TR+   k:    45     num_inner:   255     f: +1.140058e+03   |grad|: 6.445384e-06   exceeded trust region\n",
            "REJ TR-   k:    46     num_inner:   242     f: +1.140058e+03   |grad|: 6.445384e-06   exceeded trust region\n",
            "acc TR+   k:    47     num_inner:   228     f: +1.140058e+03   |grad|: 1.666072e-06   exceeded trust region\n",
            "acc       k:    48     num_inner:   209     f: +1.140058e+03   |grad|: 5.417503e-06   exceeded trust region\n",
            "acc TR+   k:    49     num_inner:   240     f: +1.140058e+03   |grad|: 6.105868e-06   exceeded trust region\n",
            "REJ TR-   k:    50     num_inner:   254     f: +1.140058e+03   |grad|: 6.105868e-06   exceeded trust region\n",
            "acc TR+   k:    51     num_inner:   246     f: +1.140058e+03   |grad|: 2.572730e-06   exceeded trust region\n",
            "acc       k:    52     num_inner:   222     f: +1.140058e+03   |grad|: 5.444834e-06   exceeded trust region\n",
            "acc TR+   k:    53     num_inner:   235     f: +1.140058e+03   |grad|: 5.034121e-06   exceeded trust region\n",
            "REJ TR-   k:    54     num_inner:   255     f: +1.140058e+03   |grad|: 5.034121e-06   exceeded trust region\n",
            "acc TR+   k:    55     num_inner:   239     f: +1.140058e+03   |grad|: 2.755474e-06   exceeded trust region\n",
            "acc       k:    56     num_inner:   204     f: +1.140058e+03   |grad|: 4.727709e-06   exceeded trust region\n",
            "acc TR+   k:    57     num_inner:   267     f: +1.140058e+03   |grad|: 6.321116e-06   exceeded trust region\n",
            "REJ TR-   k:    58     num_inner:   223     f: +1.140058e+03   |grad|: 6.321116e-06   exceeded trust region\n",
            "acc TR+   k:    59     num_inner:   215     f: +1.140058e+03   |grad|: 2.613309e-06   exceeded trust region\n",
            "acc       k:    60     num_inner:   205     f: +1.140058e+03   |grad|: 5.331412e-06   exceeded trust region\n",
            "acc TR+   k:    61     num_inner:   229     f: +1.140058e+03   |grad|: 7.062096e-06   exceeded trust region\n",
            "REJ TR-   k:    62     num_inner:   201     f: +1.140058e+03   |grad|: 7.062096e-06   exceeded trust region\n",
            "acc TR+   k:    63     num_inner:   194     f: +1.140058e+03   |grad|: 2.952249e-06   exceeded trust region\n",
            "acc       k:    64     num_inner:   220     f: +1.140058e+03   |grad|: 3.737316e-06   exceeded trust region\n",
            "acc TR+   k:    65     num_inner:   256     f: +1.140058e+03   |grad|: 5.544505e-06   exceeded trust region\n",
            "REJ TR-   k:    66     num_inner:   228     f: +1.140058e+03   |grad|: 5.544505e-06   exceeded trust region\n",
            "acc TR+   k:    67     num_inner:   214     f: +1.140058e+03   |grad|: 2.988486e-06   exceeded trust region\n",
            "acc       k:    68     num_inner:   223     f: +1.140058e+03   |grad|: 8.390461e-06   exceeded trust region\n",
            "acc TR+   k:    69     num_inner:   257     f: +1.140058e+03   |grad|: 3.290934e-06   exceeded trust region\n",
            "REJ TR-   k:    70     num_inner:   264     f: +1.140058e+03   |grad|: 3.290934e-06   exceeded trust region\n",
            "acc TR+   k:    71     num_inner:   254     f: +1.140058e+03   |grad|: 1.228651e-06   exceeded trust region\n",
            "acc       k:    72     num_inner:   154     f: +1.140058e+03   |grad|: 9.587615e-06   exceeded trust region\n",
            "acc TR+   k:    73     num_inner:   270     f: +1.140058e+03   |grad|: 6.058291e-06   exceeded trust region\n",
            "REJ TR-   k:    74     num_inner:   215     f: +1.140058e+03   |grad|: 6.058291e-06   exceeded trust region\n",
            "acc TR+   k:    75     num_inner:   206     f: +1.140058e+03   |grad|: 2.304765e-06   exceeded trust region\n",
            "acc       k:    76     num_inner:   205     f: +1.140058e+03   |grad|: 8.254015e-06   exceeded trust region\n",
            "acc TR+   k:    77     num_inner:   225     f: +1.140058e+03   |grad|: 8.706789e-06   exceeded trust region\n",
            "REJ TR-   k:    78     num_inner:   256     f: +1.140058e+03   |grad|: 8.706789e-06   exceeded trust region\n",
            "acc TR+   k:    79     num_inner:   240     f: +1.140058e+03   |grad|: 1.666329e-06   exceeded trust region\n",
            "acc       k:    80     num_inner:   176     f: +1.140058e+03   |grad|: 8.262014e-06   exceeded trust region\n",
            "acc TR+   k:    81     num_inner:   259     f: +1.140058e+03   |grad|: 3.646822e-06   exceeded trust region\n",
            "REJ TR-   k:    82     num_inner:   253     f: +1.140058e+03   |grad|: 3.646822e-06   exceeded trust region\n",
            "acc TR+   k:    83     num_inner:   233     f: +1.140058e+03   |grad|: 1.794664e-06   exceeded trust region\n",
            "acc       k:    84     num_inner:   223     f: +1.140058e+03   |grad|: 4.688861e-06   exceeded trust region\n",
            "acc TR+   k:    85     num_inner:   260     f: +1.140058e+03   |grad|: 4.624598e-06   exceeded trust region\n",
            "REJ TR-   k:    86     num_inner:   252     f: +1.140058e+03   |grad|: 4.624598e-06   exceeded trust region\n",
            "acc TR+   k:    87     num_inner:   226     f: +1.140058e+03   |grad|: 1.556547e-06   exceeded trust region\n",
            "acc       k:    88     num_inner:   201     f: +1.140058e+03   |grad|: 6.138585e-06   exceeded trust region\n",
            "acc TR+   k:    89     num_inner:   269     f: +1.140058e+03   |grad|: 4.215399e-06   exceeded trust region\n",
            "REJ TR-   k:    90     num_inner:   243     f: +1.140058e+03   |grad|: 4.215399e-06   exceeded trust region\n",
            "acc TR+   k:    91     num_inner:   228     f: +1.140058e+03   |grad|: 1.671569e-06   exceeded trust region\n",
            "acc       k:    92     num_inner:   222     f: +1.140058e+03   |grad|: 4.677648e-06   exceeded trust region\n",
            "acc TR+   k:    93     num_inner:   270     f: +1.140058e+03   |grad|: 4.287151e-06   exceeded trust region\n",
            "REJ TR-   k:    94     num_inner:   263     f: +1.140058e+03   |grad|: 4.287151e-06   exceeded trust region\n",
            "acc TR+   k:    95     num_inner:   237     f: +1.140058e+03   |grad|: 1.255983e-06   exceeded trust region\n",
            "acc       k:    96     num_inner:   169     f: +1.140058e+03   |grad|: 5.206012e-06   exceeded trust region\n",
            "acc TR+   k:    97     num_inner:   242     f: +1.140058e+03   |grad|: 3.967505e-06   exceeded trust region\n",
            "REJ TR-   k:    98     num_inner:   233     f: +1.140058e+03   |grad|: 3.967505e-06   exceeded trust region\n",
            "acc TR+   k:    99     num_inner:   223     f: +1.140058e+03   |grad|: 2.641059e-06   exceeded trust region\n",
            "acc       k:   100     num_inner:   230     f: +1.140058e+03   |grad|: 7.610201e-06   exceeded trust region\n",
            "Terminated - max iterations reached after 1954.89 seconds.\n",
            "\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0y-X3MfAluBj",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 68
        },
        "outputId": "4e5d9c07-40ac-4a14-986f-ccd11a992435"
      },
      "source": [
        "#  TEST THE CHRISTOFFEL GAMMA\n",
        "dvec = np.array([0, 30, 2, 1])\n",
        "p = dvec.shape[0]-1\n",
        "alpha = randint(1, 10, (p, p+1)) * .1\n",
        "\n",
        "dvec[0] = 1000 - dvec[1:].sum()\n",
        "man = RealFlag(dvec, alpha=alpha)\n",
        "X = man.rand()\n",
        "eta = man.randvec(X)\n",
        "g2 = man.christoffel_gamma(X, eta, eta)\n",
        "\n",
        "egrad = man._rand_ambient()\n",
        "ehess = man._rand_ambient()\n",
        "ehess_val = man.base_inner_ambient(ehess, eta)\n",
        "print(man.rhess02_alt(X, eta, eta, egrad, ehess_val))\n",
        "print(man.rhess02(X, eta, eta, egrad, ehess))\n",
        "print(ehess_val - man.base_inner_ambient(g2, egrad))\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "-0.8043898571290166\n",
            "-0.8043898571290166\n",
            "-0.8043898571290166\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}