working_points.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1e9aff8e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Changes to local modules might not update if Kernel is not restarted'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"Errors might remain if Kernel is not restarted\"\"\"\n",
    "\"\"\"Changes to local modules might not update if Kernel is not restarted\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "56e53226",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5 4\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "#### In the case of correspondence analysis, \n",
    "#### the coordinates representations of the points seems to be the right scalling\n",
    "Data = np.array([[21, 143, 15, 31],\n",
    "                 [10, 3, 156, 34],\n",
    "                 [250, 10, 14, 1],\n",
    "                 [120, 20, 1, 370],\n",
    "                 [5, 3, 12, 145]\n",
    "                ])\n",
    "M = Data.shape[0]\n",
    "N = Data.shape[1]\n",
    "print(M, N)\n",
    "joint_xy = Data/np.sum(Data)\n",
    "marj_x = np.sum(joint_xy, axis = 1)\n",
    "marj_y = np.sum(joint_xy, axis = 0)\n",
    "\n",
    "Dev_mat  =  joint_xy - marj_x[:, np.newaxis]*marj_y\n",
    "# Standardize the deviation matrix\n",
    "Dr = np.diag(1/np.sqrt(marj_x))\n",
    "Dc = np.diag(1/np.sqrt(marj_y))\n",
    "\n",
    "# standardized residuals\n",
    "D = Dr.dot(Dev_mat.dot(Dc))\n",
    "\n",
    "# Compute CA coords \n",
    "X = (Dr).dot(D)\n",
    "Y = (D.dot(Dc)).T # just bring them on the rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "84acc275",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function works pretty well from our previous projects\n",
    "B = Dr.dot(D.dot(Dc)) ## Pearson ratio(x_m, y_n) - 1 \n",
    "fXY = np.exp(- B) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4d4233d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy as sp\n",
    "# compute distance/proximity within sets of points (must be Euclidean distance matrices)\n",
    "DX = sp.spatial.distance.pdist(X)\n",
    "DX = sp.spatial.distance.squareform(DX)\n",
    "DY = sp.spatial.distance.pdist(Y)\n",
    "DY = sp.spatial.distance.squareform(DY)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9e018b1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is directly the norms of the points\n",
    "UX = np.linalg.norm(X, axis = 1)\n",
    "UY = np.linalg.norm(Y, axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ac9bf411",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5 1 2.0\n"
     ]
    }
   ],
   "source": [
    "# Provide parameters (c1, c2, c3) as dictionary, using conditions (i), (ii), and (iii) of Theorem 3.1 as a guideline\n",
    "c1, c2 = 1/2, 1 # or c2 = 2 \n",
    "a = 1. - 1./(M+N)\n",
    "b = 2.*c2/(M+N)\n",
    "c3 = ((2.*c1 + c2) - b)/a\n",
    "c = {\"c1\":c1, \"c2\":c2, \"c3\":c3} \n",
    "print(c1, c2, c3) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "58c72dab",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Is the cosine law matrix is symmetric: True\n",
      "[38.40708368 18.60563925 17.56353003 14.45702441 14.15110768 12.30342337\n",
      " 11.79638296 10.65062321 10.18650912  2.25931586  0.        ]\n"
     ]
    }
   ],
   "source": [
    "# Compute cosine law matrix\n",
    "import Methods.CosLM as CosLM\n",
    "COS_MAT, c1, c2, c3, zeta_f = CosLM.CosLM(DX, DY, UX, UY, fXY, c) # w1 = x1\n",
    "print(\"Is the cosine law matrix is symmetric:\", np.all(np.isclose(COS_MAT,COS_MAT.T)))\n",
    "sigma, U = sp.linalg.eigh(COS_MAT)\n",
    "sigma = np.real(sigma) # if COS_MAT is symmetric, thus imaginary number are supposed to be zero or numerical zeros\n",
    "sigma[np.isclose(sigma, np.zeros(len(sigma)))] = 0\n",
    "print(sigma[np.argsort(sigma)[::-1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a35e55ee",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Is COS_MAT PSD: True\n"
     ]
    }
   ],
   "source": [
    "print(\"Is COS_MAT PSD:\", np.all(sigma>=0))\n",
    "if (not np.all(sigma>=0)):\n",
    "    print(\"Check the conditions in Theorem 3.1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "15584a44",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Check theory on distances for the set X True\n",
      "Check theory on distances for the set Y True\n",
      "Check theory on distance to orgin for the set X: True\n",
      "Check theory on distance to orgin for the set Y: True\n",
      "Check theory on distance between points of sets: True\n"
     ]
    }
   ],
   "source": [
    "### check distances and proximities in Euclidean embedding ###\n",
    "W = U.dot(np.diag(np.sqrt(sigma)))\n",
    "EmbX = W[1:M+1, :]\n",
    "EmbY = W[M+2:, :]\n",
    "o = W[M+1, :]\n",
    "DMat = sp.spatial.distance.pdist(W[1:, :]) \n",
    "DMat = sp.spatial.distance.squareform(DMat)\n",
    "# Within sets \n",
    "DXe = DMat[:M, :M]\n",
    "DYe = DMat[M+1:, M+1:]\n",
    "\n",
    "print(\"Check theory on distances for the set X\", np.all(np.isclose(DXe, np.sqrt(DX**2 + c3*zeta_f) - np.sqrt(np.diag(c3*zeta_f*np.ones(DX.shape[0])))))) # diagonal elements of distance matrices are always 0s\n",
    "print(\"Check theory on distances for the set Y\", np.all(np.isclose(DYe, np.sqrt(DY**2 + c3*zeta_f) - np.sqrt(np.diag(c3*zeta_f*np.ones(DY.shape[0])))))) # diagonal elements of distance matrices are always 0s\n",
    "print(\"Check theory on distance to orgin for the set X:\", np.all(np.isclose(DMat[M, :M], np.sqrt(UX**2 + c3*zeta_f))))\n",
    "print(\"Check theory on distance to orgin for the set Y:\", np.all(np.isclose(DMat[M, M+1:], np.sqrt(UY**2 + c3*zeta_f))))\n",
    "\n",
    "# Between sets\n",
    "fXYe = DMat[:M, M+1:]\n",
    "print(\"Check theory on distance between points of sets:\", np.all(np.isclose(fXYe, np.sqrt(fXY**2 + c3*zeta_f))))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16ff1182",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}