{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# eda08_00 clear all variables and import vatious modules\n",
    "%reset -f\n",
    "import os\n",
    "from datetime import date\n",
    "from math import exp, pi, sin, cos, sqrt, floor, ceil, atan2\n",
    "import numpy as np\n",
    "import scipy.sparse.linalg as las\n",
    "from scipy import sparse\n",
    "import scipy.linalg as la\n",
    "import scipy.signal as sg\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import patches as pch\n",
    "import matplotlib.cm as cm\n",
    "from matplotlib.colors import ListedColormap\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "# eda_draw function makes a \"pictorial matrix equation\"\n",
    "# arguments are vectors, matrices and strings\n",
    "# which are plotted in the order that the appear\n",
    "# except that strings starting with 'title ' are plotted\n",
    "# under the subseqeunt matrix or vector\n",
    "# always returns a status of 1\n",
    "\n",
    "def eda_draw(*argv):\n",
    "    bw = np.zeros((256,4));\n",
    "    v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;\n",
    "    bw[:,0] = v;\n",
    "    bw[:,1] = v;\n",
    "    bw[:,2] = v;\n",
    "    bw[:,3] = np.ones(256);\n",
    "    bwcmap = ListedColormap(bw);\n",
    "    # size of plot\n",
    "    W = 16;\n",
    "    H = 4;\n",
    "    fig1 = plt.figure(1);\n",
    "    # figsize width and height in inches\n",
    "    fig1.set_size_inches(W,H);\n",
    "    ax1 = plt.subplot(1,1,1);\n",
    "    plt.axis([0, W, -H/2, H/2]);\n",
    "    plt.axis('off');\n",
    "    LM = W/6;    # matrix width and heoght\n",
    "    LV = W/40;   # vector width\n",
    "    FS = 0.12;    # character width\n",
    "    TO = 0.4;    # title vertical offset\n",
    "    SP = 0.2;    # space between objects\n",
    "    LS = 0.2;    # leading space\n",
    "    p = LS; # starting x-position\n",
    "    istitle=0; # flags presence of a title\n",
    "    for a in argv:\n",
    "        if isinstance(a,np.ndarray):\n",
    "            sh = np.shape(a);\n",
    "            if len(sh) == 1:  # conversion to nx1 array\n",
    "                n = sh[0];\n",
    "                m = 1;\n",
    "                ap = a;\n",
    "                a = np.zeros((n,1));\n",
    "                a[:,0] = ap;\n",
    "            else:\n",
    "                n = sh[0];\n",
    "                m = sh[1];\n",
    "            if m==1:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LV;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LV;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "            else:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LM;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LM;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "        elif isinstance(a,str):\n",
    "            ns = len(a);\n",
    "            istitle=0;\n",
    "            if( ns>=6 ):\n",
    "                if 'title ' in a[0:6]:\n",
    "                    istitle=1;\n",
    "                    titlestr=a[6:];\n",
    "            if( istitle != 1):\n",
    "                plt.text(p,0,a);\n",
    "                p = p + ns*FS + SP;\n",
    "    plt.show();\n",
    "    return 1;\n",
    "\n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec1(v) from entity v that is array-like,\n",
    "# including a number, a list of numbers, a tuple\n",
    "# of numbers, an N-by-0 np array and a N-by1 np array.\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it.\n",
    "# This version takes just one input argment.\n",
    "# whereas c=eda_cvec(v1,v2,...) concatenates\n",
    "# many argiments.\n",
    "def eda_cvec1(v):\n",
    "    if isinstance(v, int):\n",
    "        w = np.zeros((1,1),dtype=int);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, float):\n",
    "        w = np.zeros((1,1),dtype=float);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, complex):\n",
    "        w = np.zeros((1,1),dtype=complex);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, np.ndarray):\n",
    "        s = np.shape(v);\n",
    "        if len(s) == 1:\n",
    "            return np.copy(np.reshape(v,(s[0],1)));\n",
    "        else:\n",
    "            [r,c]=s;\n",
    "            if( c==1 ):\n",
    "                return(np.copy(v));\n",
    "            elif(r==1):\n",
    "                return(np.copy(v.T));\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: %d by %d ndarray not allowed\" % (r, c));\n",
    "    elif isinstance(v, list):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: list contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, tuple):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: tuple contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    else:\n",
    "        raise TypeError(\"eda_cvec: %s not supported\" % type(v));\n",
    "        \n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec(v1, v2, ...) from a list of several\n",
    "# array-like entities v1, v2, including a number\n",
    "# a list of numbers, a tuple of numbers, an N-by-0 np array\n",
    "# and a N-by-1 np array. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context.\n",
    "# This version concatenates many argments,\n",
    "# whereas c=eda_cvec1(v1) takes just one argiment.\n",
    "# I recommend always using eda_cvec(v1, v2, ...)\n",
    "def eda_cvec(*argv):\n",
    "    t = int;\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        Nt = Nt + N;\n",
    "        if (t==int) and isinstance(v[0,0],float):\n",
    "            t=float;\n",
    "        elif isinstance(v[0,0],complex):\n",
    "            t=complex;\n",
    "    w = np.zeros((Nt,1),dtype=t);\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        w[Nt:Nt+N,0] = v[0:N,0];\n",
    "        Nt = Nt + N;\n",
    "    return w;\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 36764 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda08_01\n",
    "# factors shown on terniary diagram\n",
    "\n",
    "# vectrs of constants needed to make ternary diagram\n",
    "a = eda_cvec( [0.0,0.0] );\n",
    "b = eda_cvec( [1.0,0.0] );\n",
    "c = eda_cvec( [0.5,0.5*sqrt(3)] );\n",
    "\n",
    "# factor 1\n",
    "f1 = eda_cvec( [0.15, 0.25, 0.60] );\n",
    "f1 = f1/np.sum(f1); # normalize factor  1\n",
    "\n",
    "# factor 2\n",
    "f2 = eda_cvec( [0.25, 0.45, 0.30] );\n",
    "f2 = f2/np.sum(f2); # normalize factor 2\n",
    "\n",
    "# terniary diagram\n",
    "fig1 = plt.figure(1,figsize=(6,6));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, 1, 0, 1] );\n",
    "plt.axis('off');\n",
    "# bounding triangle\n",
    "plt.plot( [a[0,0], b[0,0]], [a[1,0], b[1,0]], 'k-');\n",
    "plt.plot( [b[0,0], c[0,0]], [b[1,0], c[1,0]], 'k-');\n",
    "plt.plot( [c[0,0], a[0,0]], [c[1,0], a[1,0]], 'k-');\n",
    "plt.text(a[0,0]-0.05, a[1,0],'A',horizontalalignment='center');\n",
    "plt.text(b[0,0]+0.05, b[1,0],'B',horizontalalignment='center');\n",
    "plt.text(c[0,0], c[1,0]+0.05,'C',horizontalalignment='center');\n",
    "# first set of grid lines\n",
    "L=10;\n",
    "x = np.zeros((2,2))\n",
    "for j in range(L):\n",
    "    for k in range(2):\n",
    "        v1 = float(j+1)/float(L);\n",
    "        v2 = (1.0-v1)*float(k);\n",
    "        v3 = 1.0 - (v1+v2);\n",
    "        x[:,k]=(v1*a+v2*b+v3*c).ravel();\n",
    "    plt.plot( [x[0,0], x[0,1]], [x[1,0], x[1,1]], 'k:' );\n",
    "# second set of grid lines\n",
    "L=10;\n",
    "for j in range(L):\n",
    "    for k in range(2):\n",
    "        v3 = float(j+1)/float(L);\n",
    "        v1 =(1-v3)*float(k);\n",
    "        v2 = 1.0 - (v1+v3);\n",
    "        x[:,k]=(v1*a+v2*b+v3*c).ravel();\n",
    "    plt.plot( [x[0,0], x[0,1]], [x[1,0], x[1,1]], 'k:' );\n",
    "#third set of grid lines\n",
    "L=10;\n",
    "for j in range(L):\n",
    "    for k in range(2):\n",
    "        v2 = float(j+1)/float(L);\n",
    "        v1 = (1-v2)*float(k);\n",
    "        v3 = 1.0 - (v1+v2);\n",
    "        x[:,k]=(v1*a+v2*b+v3*c).ravel();\n",
    "    plt.plot( [x[0,0], x[0,1]], [x[1,0], x[1,1]], 'k:' );\n",
    "    \n",
    "# plot factors\n",
    "y = f1[0]*a+f1[1]*b+f1[2]*c;\n",
    "plt.plot( y[0], y[1], 'k*' );\n",
    "plt.text(y[0], y[1]+0.025,'f1',horizontalalignment='center');\n",
    "y = f2[0]*a+f2[1]*b+f2[2]*c;\n",
    "plt.plot( y[0], y[1], 'k*' );\n",
    "plt.text(y[0], y[1]-0.05,'f2',horizontalalignment='center');\n",
    "\n",
    "# create and plot samples lying between f1 and f2\n",
    "N=10;\n",
    "for j in range(N-1):\n",
    "    c1=float(j+1)/float(N);\n",
    "    c2=1-c1;\n",
    "    s = c1*f1 + c2*f2;\n",
    "    y = s[0]*a+s[1]*b+s[2]*c;\n",
    "    plt.plot( y[0], y[1], 'ko' );\n",
    "    \n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 12604 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda08_02\n",
    "# factors shown on terniary diagram\n",
    "# this is eda08_01 w/o the grid lines and w/o annotation of factors\n",
    "\n",
    "# constants needed to make ternary diagram\n",
    "a = eda_cvec( [0.0,0.0] );\n",
    "b = eda_cvec( [1.0,0.0] );\n",
    "c = eda_cvec( [0.5,0.5*sqrt(3)] );\n",
    "\n",
    "# factor 1\n",
    "f1 = eda_cvec( [0.15, 0.25, 0.60] );\n",
    "f1 = f1/np.sum(f1); # normalize factor 1\n",
    "\n",
    "# factor 2\n",
    "f2 = eda_cvec( [0.25, 0.45, 0.30] );\n",
    "f2 = f2/np.sum(f2); # normalize factor 2\n",
    "\n",
    "# terniary diagram\n",
    "fig1 = plt.figure(1,figsize=(6,6));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, 1, 0, 1] );\n",
    "plt.axis('off');\n",
    "# bounding triangle\n",
    "plt.plot( [a[0,0], b[0,0]], [a[1,0], b[1,0]], 'k-');\n",
    "plt.plot( [b[0,0], c[0,0]], [b[1,0], c[1,0]], 'k-');\n",
    "plt.plot( [c[0,0], a[0,0]], [c[1,0], a[1,0]], 'k-');\n",
    "plt.text(a[0,0]-0.05, a[1,0],'A',horizontalalignment='center');\n",
    "plt.text(b[0,0]+0.05, b[1,0],'B',horizontalalignment='center');\n",
    "plt.text(c[0,0], c[1,0]+0.05,'C',horizontalalignment='center');\n",
    "\n",
    "# plot factors\n",
    "y = f1[0]*a+f1[1]*b+f1[2]*c;\n",
    "plt.plot( y[0], y[1], 'k*' );\n",
    "y = f2[0]*a+f2[1]*b+f2[2]*c;\n",
    "plt.plot( y[0], y[1], 'k*' );\n",
    "\n",
    "# create and plot samples lying between f1 and f2\n",
    "N=10;\n",
    "for j in range(N-1):\n",
    "    c1=float(j+1)/float(N);\n",
    "    c2=1-c1;\n",
    "    s = c1*f1 + c2*f2;\n",
    "    y = s[0]*a+s[1]*b+s[2]*c;\n",
    "    plt.plot( y[0], y[1], 'ko' );\n",
    "    \n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 12724 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda08_03\n",
    "# factors shown on terniary diagram\n",
    "# this is eda08_02 with a different choice of factors\n",
    "\n",
    "# constants needed to make ternary diagram\n",
    "a = eda_cvec( [0.0,0.0] );\n",
    "b = eda_cvec( [1.0,0.0] );\n",
    "c = eda_cvec( [0.5,0.5*sqrt(3)] );\n",
    "\n",
    "# factor 1\n",
    "f1 = eda_cvec( [0.15, 0.25, 0.60] );\n",
    "f1 = f1/np.sum(f1); # normalize factor  1\n",
    "\n",
    "# factor 2\n",
    "f2 = eda_cvec( [0.25, 0.45, 0.30] );\n",
    "f2 = f2/np.sum(f2); # normalie factor 2\n",
    "\n",
    "f1a = 0.5*(f1+f2); # factor 1a\n",
    "f2a = f1+0.5*(f1-f2); # factor 2a \n",
    "\n",
    "# terniary diagram\n",
    "fig1 = plt.figure(1,figsize=(6,6));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, 1, 0, 1] );\n",
    "plt.axis('off');\n",
    "# bounding triangle\n",
    "plt.plot( [a[0,0], b[0,0]], [a[1,0], b[1,0]], 'k-');\n",
    "plt.plot( [b[0,0], c[0,0]], [b[1,0], c[1,0]], 'k-');\n",
    "plt.plot( [c[0,0], a[0,0]], [c[1,0], a[1,0]], 'k-');\n",
    "plt.text(a[0,0]-0.05, a[1,0],'A',horizontalalignment='center');\n",
    "plt.text(b[0,0]+0.05, b[1,0],'B',horizontalalignment='center');\n",
    "plt.text(c[0,0], c[1,0]+0.05,'C',horizontalalignment='center');\n",
    "\n",
    "# create and plot samples lying between f1 and f2\n",
    "N=10;\n",
    "for j in range(N-1):\n",
    "    c1=float(j+1)/float(N);\n",
    "    c2=1-c1;\n",
    "    s = c1*f1 + c2*f2;\n",
    "    y = s[0]*a+s[1]*b+s[2]*c;\n",
    "    plt.plot( y[0], y[1], 'ko' );\n",
    "    \n",
    "# plot factors\n",
    "y = f1a[0]*a+f1a[1]*b+f1a[2]*c;\n",
    "plt.plot( y[0], y[1], 'k*', color=(0.7,0.7,0.7) );\n",
    "y = f2a[0]*a+f2a[1]*b+f2a[2]*c;\n",
    "plt.plot( y[0], y[1], 'k*', color=(0.7,0.7,0.7) );\n",
    "    \n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max relative error in CF   2.9341870488937076e-15\n",
      "factor 0\n",
      "    SiO2 -0.908829\n",
      "    TiO2 -0.024638\n",
      "    Al203 -0.275168\n",
      "    FeO-total -0.177851\n",
      "    MgO -0.141341\n",
      "    CaO -0.209989\n",
      "    Na2O -0.044611\n",
      "    K2O -0.003430\n",
      "factor 1\n",
      "    SiO2 0.007684\n",
      "    TiO2 -0.037474\n",
      "    Al203 -0.301583\n",
      "    FeO-total -0.018421\n",
      "    MgO 0.923193\n",
      "    CaO -0.226917\n",
      "    Na2O -0.058457\n",
      "    K2O -0.007204\n",
      "factor 2\n",
      "    SiO2 -0.161689\n",
      "    TiO2 -0.126343\n",
      "    Al203 0.567828\n",
      "    FeO-total -0.659205\n",
      "    MgO 0.255748\n",
      "    CaO 0.365682\n",
      "    Na2O -0.041738\n",
      "    K2O -0.006464\n",
      "factor 3\n",
      "    SiO2 0.209819\n",
      "    TiO2 0.151367\n",
      "    Al203 0.176021\n",
      "    FeO-total -0.427461\n",
      "    MgO -0.118643\n",
      "    CaO -0.780043\n",
      "    Na2O 0.302367\n",
      "    K2O 0.073403\n",
      "factor 4\n",
      "    SiO2 0.309495\n",
      "    TiO2 -0.100476\n",
      "    Al203 -0.670083\n",
      "    FeO-total -0.585155\n",
      "    MgO -0.195193\n",
      "    CaO 0.207980\n",
      "    Na2O -0.145318\n",
      "    K2O 0.015035\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 14928 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 54168 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda08_04\n",
    "# factor analysis on Atlantic Rocks dataset\n",
    "# using singular value decomposition\n",
    "\n",
    "# load data\n",
    "S = np.genfromtxt('rocks.txt', delimiter='\\t')\n",
    "N, M = np.shape(S);\n",
    "Smax = np.amax(S);\n",
    "\n",
    "# break out data into vectors\n",
    "sio2 = eda_cvec( S[:,0] );   # SiO2\n",
    "tio2 = eda_cvec( S[:,1] );   # TiO2\n",
    "als03 = eda_cvec( S[:,2] );  # Al203\n",
    "feot = eda_cvec( S[:,3] );   # FeO-total\n",
    "mgo = eda_cvec( S[:,4] );    # MgO\n",
    "cao = eda_cvec( S[:,5] );    # CaO\n",
    "na20 = eda_cvec( S[:,6] );   # Na2O\n",
    "k20 = eda_cvec( S[:,7] );    # K2O\n",
    "\n",
    "# compute factors F and factor loadings C using singular value decompostion\n",
    "[U, sigma, VT] = la.svd(S,full_matrices=False);\n",
    "sh = np.shape(sigma);\n",
    "Ns = sh[0];\n",
    "F = np.copy(VT);\n",
    "C = np.matmul( U, np.diag(sigma) );\n",
    "\n",
    "# check error of reconstruction\n",
    "Emax = np.amax(S - np.matmul(C,F));\n",
    "print('max relative error in CF  ', Emax/Smax);\n",
    "\n",
    "# plot singular values\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, Ns+1, 0, np.max(sigma)] );\n",
    "plt.plot( np.linspace(1,Ns,Ns), sigma, 'k-' );\n",
    "plt.plot( np.linspace(1,Ns,Ns), sigma, 'ko' );\n",
    "plt.title('singular values, s(i)');\n",
    "plt.xlabel('index, i');\n",
    "plt.ylabel('s(i)');\n",
    "\n",
    "# print first five factors\n",
    "for j in range(5):\n",
    "    fj = np.zeros((M,1));\n",
    "    fj[:,0] = F[j,:];\n",
    "    printstr = \"factor %d\" % (j);\n",
    "    print(printstr);\n",
    "    printstr = \"    SiO2 %f\"      % (fj[0]);\n",
    "    print(printstr);\n",
    "    printstr = \"    TiO2 %f\"      % (fj[1]);\n",
    "    print(printstr);\n",
    "    printstr = \"    Al203 %f\"     % (fj[2]);\n",
    "    print(printstr);\n",
    "    printstr = \"    FeO-total %f\" % (fj[3]);\n",
    "    print(printstr);\n",
    "    printstr = \"    MgO %f\"       % (fj[4]);\n",
    "    print(printstr);\n",
    "    printstr = \"    CaO %f\"       % (fj[5]);\n",
    "    print(printstr);\n",
    "    printstr = \"    Na2O %f\"      % (fj[6]);\n",
    "    print(printstr);\n",
    "    printstr = \"    K2O %f\"       % (fj[7]);\n",
    "    print(printstr);\n",
    "print(' ');\n",
    "\n",
    "cmin=(-50);\n",
    "cmax=50;\n",
    "fig2 = plt.figure(2,figsize=(6,6));\n",
    "ax1 = fig2.add_subplot(111, projection='3d')\n",
    "ax1.axes.set_xlim3d(left=cmin, right=cmax);\n",
    "ax1.axes.set_ylim3d(bottom=cmin, top=cmax);\n",
    "ax1.axes.set_zlim3d(bottom=cmin, top=cmax);\n",
    "ax1.view_init(30, 30)\n",
    "\n",
    "ax1.scatter( xs=C[:,1], ys=C[:,2], zs=C[:,3], color=(0,0,0) );\n",
    "ax1.plot( xs=[cmin,cmin], ys=[cmin,cmin], zs=[cmin,cmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[cmin,cmin], ys=[cmin,cmax], zs=[cmin,cmin], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[cmin,cmax], ys=[cmin,cmin], zs=[cmin,cmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[cmax,cmax], ys=[cmax,cmax], zs=[cmax,cmin], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[cmax,cmax], ys=[cmax,cmin], zs=[cmax,cmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[cmax,cmin], ys=[cmax,cmax], zs=[cmax,cmax], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[cmax,cmin], ys=[cmin,cmin], zs=[cmax,cmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[cmax,cmax], ys=[cmin,cmin], zs=[cmax,cmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[cmin,cmin], ys=[cmax,cmin], zs=[cmax,cmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[cmin,cmin], ys=[cmax,cmax], zs=[cmax,cmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[cmax,cmax], ys=[cmax,cmin], zs=[cmin,cmin], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[cmax,cmin], ys=[cmax,cmax], zs=[cmin,cmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.set_xlabel('C(2)');\n",
    "ax1.set_ylabel('C(3)');\n",
    "ax1.set_zlabel('C(4)');\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fA  0.500 0.500. 0.500 0.500\n",
      "fB  0.500 -0.500. 0.500 -0.500\n",
      " \n",
      "fAp 0.707 0.000. 0.707 0.000\n",
      "FBp 0.000 -0.707. 0.000 -0.707\n"
     ]
    }
   ],
   "source": [
    "# eda08_05\n",
    "# simple example of varimax procedure applued to 2 factors\n",
    "\n",
    "# define two non-spiky factors, fA and fB\n",
    "M=4;\n",
    "fA = eda_cvec( [ 1.0,  1.0,  1.0,  1.0 ] );\n",
    "fB = eda_cvec( [ 1.0, -1.0,  1.0, -1.0 ]);\n",
    "fA = fA / sqrt( np.matmul(fA.T, fA )); # normalize to unit length\n",
    "fB = fB / sqrt( np.matmul(fB.T, fB )); # normalize to unit length\n",
    "\n",
    "# standard varimax procedure to determine rotation angle q\n",
    "u = np.power(fA,2) - np.power(fB,2);\n",
    "v = 2.0 * np.multiply(fA, fB);\n",
    "\n",
    "A = 2*M*np.matmul(u.T,v);\n",
    "B = np.sum(u)*np.sum(v);\n",
    "top = A - B;\n",
    "\n",
    "C = M*(np.matmul(u.T,u)-np.matmul(v.T,v));\n",
    "D = (np.sum(u)**2)-(np.sum(v)**2);\n",
    "bot =  C - D;\n",
    "\n",
    "q = 0.25 * atan2(top,bot);\n",
    "\n",
    "# rotate factors\n",
    "cq = cos(q);\n",
    "sq = sin(q);\n",
    "fAp =  cq*fA + sq*fB;\n",
    "fBp = -sq*fA + cq*fB;\n",
    "\n",
    "# display results\n",
    "print(\"fA  %.3f %.3f. %.3f %.3f\" % (fA[0], fA[1], fA[2], fA[3]) );\n",
    "print(\"fB  %.3f %.3f. %.3f %.3f\" % (fB[0], fB[1], fB[2], fB[3]) );\n",
    "print(' ');\n",
    "print(\"fAp %.3f %.3f. %.3f %.3f\" % (fAp[0], fAp[1], fAp[2], fAp[3]) );\n",
    "print(\"FBp %.3f %.3f. %.3f %.3f\" % (fBp[0], fBp[1], fBp[2], fBp[3]) );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max relative error in CF   0.05663026082140537\n",
      "max relative error in CPFP 0.056630260821405376\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 5104 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda08_06\n",
    "# factor analysis on Atlantic Rocks dataset\n",
    "# using singular value decomposition\n",
    "# and the varimax procedure\n",
    "\n",
    "# load data\n",
    "S = np.genfromtxt('rocks.txt', delimiter='\\t')\n",
    "N, M = np.shape(S);\n",
    "Smax = np.amax(S);\n",
    "\n",
    "# compute factors and factor loadings using singular value decompostion\n",
    "[U, s, VT] = la.svd(S,full_matrices=False);\n",
    "\n",
    "# check error of reconstruction\n",
    "# Emax = np.amax(S - np.matmul(np.matmul(U,np.diag(s)),VT));\n",
    "# print('max relative error in USVT', Emax/Smax);\n",
    "\n",
    "# keep only first five singular values\n",
    "P=5;\n",
    "F = np.copy(VT[0:P,:]);\n",
    "C = np.matmul( U[:,0:P], np.diag(s[0:P]));\n",
    "\n",
    "# check error of reconstruction\n",
    "Emax = np.amax(S - np.matmul(C,F));\n",
    "print('max relative error in CF  ', Emax/Smax);\n",
    "\n",
    "# spike these factors using the varimax procedure\n",
    "k = [1, 2, 3, 4];\n",
    "Nk = len(k);\n",
    "\n",
    "# initial rotated factor matrix, FP, abd rotation matrix, MR\n",
    "MR=np.identity(P);\n",
    "FP=np.copy(F);\n",
    "\n",
    "fA = np.zeros((M,1));\n",
    "fB = np.zeros((M,1));\n",
    "fAp = np.zeros((M,1));\n",
    "fBp = np.zeros((M,1));\n",
    "mA = np.zeros((P,1));\n",
    "mB = np.zeros((P,1));\n",
    "mAp = np.zeros((P,1));\n",
    "mBp = np.zeros((P,1));\n",
    "for iter in range(3):\n",
    "    for ii in range(Nk):\n",
    "        for jj in range(ii+1,Nk):\n",
    "            # spike factors i and j\n",
    "            i=k[ii];\n",
    "            j=k[jj];\n",
    "            # copy factors from matrix to vectors\n",
    "            fA[:,0] = np.copy(FP[i,:]);\n",
    "            fB[:,0] = np.copy(FP[j,:]);\n",
    "            # standard varimax procedure to determine rotation angle q\n",
    "            u = np.power(fA,2) - np.power(fB,2);\n",
    "            v = 2.0 * np.multiply(fA, fB);\n",
    "            AA = 2*M*np.matmul(u.T,v);\n",
    "            BB = np.sum(u)*np.sum(v);\n",
    "            top = AA - BB;\n",
    "            CC = M*(np.matmul(u.T,u)-np.matmul(v.T,v));\n",
    "            DD = (np.sum(u)**2)-(np.sum(v)**2);\n",
    "            bot =  CC - DD;\n",
    "            q = 0.25 * atan2(top,bot);\n",
    "            # rotate factors\n",
    "            cq = cos(q);\n",
    "            sq = sin(q);\n",
    "            fAp =  cq*fA + sq*fB;\n",
    "            fBp = -sq*fA + cq*fB;\n",
    "            # put back into factor matrix, FP\n",
    "            FP[i,:] = fAp.ravel();\n",
    "            FP[j,:] = fBp.ravel();\n",
    "            # accumulate rotation\n",
    "            mA[:,0] = np.copy(MR[i,:]);\n",
    "            mB[:,0] = np.copy(MR[j,:]);\n",
    "            mAp =  cq*mA + sq*mB;\n",
    "            mBp = -sq*mA + cq*mB;\n",
    "            MR[i,:] = mAp.ravel();\n",
    "            MR[j,:] = mBp.ravel();\n",
    "\n",
    "# new factor loadings\n",
    "CP = np.matmul(C,MR.T);\n",
    "# check error of reconstruction\n",
    "Emax = np.amax(S - np.matmul(CP,FP));\n",
    "print('max relative error in CPFP', Emax/Smax);\n",
    "\n",
    "Fa = np.abs(F);\n",
    "FPa = np.abs(FP);\n",
    "eda_draw( 'title f1', Fa[1,:], 'title f2', Fa[2,:], 'title f3', Fa[3,:], 'title f4', Fa[4,:], 'becomes', 'title fp1', FPa[1,:], 'title fp2', FPa[2,:], 'title fp3', FPa[3,:], 'title fp4', FPa[4,:]);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "witout weighting, max relative error in CF   0.18873689079879263\n",
      "with   weighting, max relative error in CF   0.1813415848375578\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 23664 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 17848 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda08_07\n",
    "# comparison of weighted and unweighted factor analysis\n",
    "\n",
    "# Part 1: Make synthetic chemical data with a wide range of concentrations\n",
    "\n",
    "# 7 elements, three factors\n",
    "M = 7;\n",
    "P = 3;\n",
    "\n",
    "# typical analytic errors\n",
    "#              [Fe        Cu        Zn        Pb        As        Ag        Au]\n",
    "se = eda_cvec( [0.010000, 0.010000, 0.010000, 0.001000, 0.001000, 0.0000010, 0.0000001 ] );\n",
    "\n",
    "# factors as column vectors\n",
    "v1 = np.zeros((M,1));\n",
    "v2 = np.zeros((M,1));\n",
    "v3 = np.zeros((M,1));\n",
    "#              [Fe        Cu        Zn        Pb        As        Ag         Au]\n",
    "v1 = eda_cvec( [0.200000, 0.004000, 0.010000, 0.010000, 0.001000, 0.0000100, 0.0000010 ] );\n",
    "v2 = eda_cvec( [0.050000, 0.180000, 0.180000, 0.030000, 0.000200, 0.0000400, 0.0000005 ] );\n",
    "v3 = eda_cvec( [0.100000, 0.080000, 0.110000, 0.010000, 0.000700, 0.0000100, 0.0000001 ] );\n",
    "\n",
    "# factor matrix\n",
    "F = np.concatenate( (v1.T, v2.T, v3.T), axis=0 );\n",
    "\n",
    "# randomly chosen loadings\n",
    "N = 10000;\n",
    "c1 = np.random.uniform(0.0,1.0,(N,1));\n",
    "c2 = np.random.uniform(0.0,1.0,(N,1));\n",
    "c3 = 1.0-(c1+c2);\n",
    "C = np.concatenate( (c1, c2, c3), axis=1 )\n",
    "\n",
    "# true sample matrix\n",
    "Strue = np.matmul(C,F);\n",
    "Sobs = np.zeros((N,M));\n",
    "\n",
    "# observed sample matrix is true plus Normally-distributed random noise\n",
    "for i in range(M):\n",
    "    Sobs[:,i] = Strue[:,i] + np.random.normal(0,se[i,0],N);\n",
    "Smax = np.amax(Sobs);\n",
    "              \n",
    "# histogram limts\n",
    "Lh = 100;\n",
    "Fe_min = 0;\n",
    "Fe_max = 5;\n",
    "Au_min = 0;\n",
    "Au_max = 5;\n",
    "\n",
    "# Part two factor analysis without weighting\n",
    "\n",
    "# factors F and loadings F\n",
    "[U1, s1, V1T] = la.svd(Sobs,full_matrices=False);\n",
    "Fpre1 = V1T; # factors\n",
    "Cpre1 = np.matmul( U1, np.diag(s1) ); # loadings\n",
    "\n",
    "# P significant factors\n",
    "Fpre1P = np.copy( Fpre1[0:P,:] );\n",
    "Cpre1P = np.copy( Cpre1[:,0:P] );\n",
    "Spre1P = np.matmul( Cpre1P, Fpre1P );\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "\n",
    "# error\n",
    "e1 = np.abs(Strue-Spre1P);  # absolute values of individual errors\n",
    "Emax = np.max(e1);          # maximum error\n",
    "print('witout weighting, max relative error in CF  ', Emax/Smax);\n",
    "\n",
    "c, e = np.histogram( e1[:,0]/se[0,0], Lh, (Fe_min,Fe_max));\n",
    "Nc = len(c);\n",
    "Ne = len(e);\n",
    "count = eda_cvec( c );\n",
    "edges = eda_cvec( e[0:Ne-1] );\n",
    "\n",
    "ax1 = plt.subplot(1,2,1);\n",
    "plt.axis( [Fe_min, Fe_max, 0, 1] );\n",
    "plt.plot(edges,count/np.max(count),'k-');\n",
    "plt.title('error in Fe');\n",
    "plt.xlabel('error/se');\n",
    "plt.ylabel('count');\n",
    "\n",
    "c, e = np.histogram( e1[:,M-1]/se[M-1,0], Lh,(Au_min,Au_max));\n",
    "Nc = len(c);\n",
    "Ne = len(e);\n",
    "count = eda_cvec( c );\n",
    "edges = eda_cvec( e[0:Ne-1] );\n",
    "\n",
    "ax1 = plt.subplot(1,2,2);\n",
    "plt.axis( [Au_min, Au_max, 0, 1] );\n",
    "plt.plot(edges,count/np.max(count),'k-');\n",
    "plt.title('Error in Au');\n",
    "plt.xlabel('error/se');\n",
    "plt.ylabel('count');\n",
    "\n",
    "# Part three factor analysis with weighting by analytic precision\n",
    "\n",
    "# factors F and loadings F\n",
    "w = np.reciprocal(se);\n",
    "[U2, s2, V2T] = la.svd( np.matmul(Sobs,np.diag(w.ravel())),full_matrices=False);\n",
    "\n",
    "Fpre2 = np.matmul(V2T,np.diag(se.ravel())); # factors\n",
    "Cpre2 = np.matmul( U2, np.diag(s2) ); # loadings\n",
    "\n",
    "# P significant factors\n",
    "Fpre2P = np.copy( Fpre2[0:P,:] );\n",
    "Cpre2P = np.copy( Cpre2[:,0:P] );\n",
    "Spre2P = np.matmul( Cpre2P, Fpre2P );\n",
    "\n",
    "fig2 = plt.figure(2,figsize=(8,4));\n",
    "\n",
    "# check error of reconstruction\n",
    "# Emax is absolute vale of difference between true and reconstructed sample matrices\n",
    "e2 = np.abs(Strue-Spre2P);  \n",
    "Emax = np.max(e2);\n",
    "print('with   weighting, max relative error in CF  ', Emax/Smax);\n",
    "\n",
    "c, e = np.histogram( e2[:,0]/se[0,0], Lh, (Fe_min,Fe_max));\n",
    "Nc = len(c);\n",
    "Ne = len(e);\n",
    "count = eda_cvec( c );\n",
    "edges = eda_cvec( e[0:Ne-1] );\n",
    "\n",
    "ax1 = plt.subplot(1,2,1);\n",
    "plt.axis( [Fe_min, Fe_max, 0, 1] );\n",
    "plt.plot(edges,count/np.max(count),'k-');\n",
    "plt.title('error in Fe');\n",
    "plt.xlabel('error/se');\n",
    "plt.ylabel('count');\n",
    "\n",
    "c, e = np.histogram( e2[:,M-1]/se[M-1,0], Lh,(Au_min,Au_max));\n",
    "Nc = len(c);\n",
    "Ne = len(e);\n",
    "count = eda_cvec( c );\n",
    "edges = eda_cvec( e[0:Ne-1] );\n",
    "\n",
    "ax1 = plt.subplot(1,2,2);\n",
    "plt.axis( [Au_min, Au_max, 0, 1] );\n",
    "plt.plot(edges,count/np.max(count),'k-');\n",
    "plt.title('Error in Au');\n",
    "plt.xlabel('error/se');\n",
    "plt.ylabel('count');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 93052 characters]\n",
      "text/plain": [
       "<Figure size 1152x720 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda08_08\n",
    "# factor analysis on Mid Atlantic Ridge Rocks dataset\n",
    "# Q-mode factors with geographical plots\n",
    "\n",
    "# load data\n",
    "D = np.genfromtxt('rocks_with_lat_lon.txt', delimiter='\\t')\n",
    "N, K = np.shape(D);\n",
    "M=K-2;\n",
    "\n",
    "# break out (lat,lon) coordinates and sample matrix S into separate vctors\n",
    "lat = eda_cvec( D[:,8] );    # lat\n",
    "lon = eda_cvec( D[:,9] );    # lon\n",
    "S = np.copy(D[:,0:M]);       # sample matric S\n",
    "\n",
    "# svd\n",
    "[U, s, VT] = la.svd(S.T,full_matrices=False);\n",
    "\n",
    "# keep only first few singular values\n",
    "P=4;\n",
    "FP = np.copy(VT[0:P,:]); # factor matrix\n",
    "CP = np.matmul( U[:,0:P], np.diag(s[0:P])); # loading matrix\n",
    "\n",
    "# map bounds\n",
    "left = -60.0;\n",
    "right = 10.0;\n",
    "bottom = -70.0;\n",
    "top = 70.0;\n",
    "\n",
    "# plot\n",
    "fig1 = plt.figure(1,figsize=(16,10));\n",
    "\n",
    "# PART 1: factor 1\n",
    "ax1 = plt.subplot(1,3,1);\n",
    "plt.axis( [left, right, bottom, top] );\n",
    "\n",
    "# plot one rectangle per sample i, with color scaled according to\n",
    "# the strength of the i-th element in factor nf=1\n",
    "nf = 1;\n",
    "FPmean = np.mean(FP[nf,:]);\n",
    "FPstd = np.std(FP[nf,:]);\n",
    "fmin = -2*FPstd\n",
    "fmax = 2*FPstd\n",
    "for i in range(N):\n",
    "    h = 3;\n",
    "    c = (FP[nf,i]-fmin)/(fmax-fmin);\n",
    "    if( c<0.0 ):\n",
    "        c=0.0;\n",
    "    elif ( c>1.0 ):\n",
    "        c=1.0;\n",
    "    ax1.add_patch( pch.Rectangle( (lon[i,0], lat[i,0]), h, h, fc=(c,c,c), ec=(0.0,0.0,0.0) ) );\n",
    "\n",
    "plt.xlabel('longitude (deg)');\n",
    "plt.ylabel('latitude (deg)');\n",
    "plt.title('Factor 1');\n",
    "\n",
    "# plot coastlines on maps, using coarse coastlines from the file coastdata.txt\n",
    "# this file consists of lists of (lat,lons) along segments of coastlines, with\n",
    "# each list deliited with a \">\".\n",
    "coastfile = 'coastdata.txt';\n",
    "fd = open(coastfile,\"r\");\n",
    "NN = 0;\n",
    "NNMAX = 10000;\n",
    "mylat = np.zeros((NNMAX,1));\n",
    "mylon = np.zeros((NNMAX,1));\n",
    "for i in range(100000):\n",
    "    line = fd.readline();\n",
    "    if not line:\n",
    "        break;\n",
    "    if( line[0] == '>' ):\n",
    "        if( NN == 0):\n",
    "            continue;\n",
    "        plt.plot( mylon[0:NN,0], mylat[0:NN,0], 'k.', color=(0.7,0.7,0.7) );\n",
    "        NN=0;\n",
    "    else:\n",
    "        latlonstr=line.split();\n",
    "        mylat[NN,0] = float(latlonstr[1]);\n",
    "        mylon[NN,0] = float(latlonstr[0]);\n",
    "        NN=NN+1;\n",
    "fd.close();\n",
    "\n",
    "# PART 1: factor 2\n",
    "ax2 = plt.subplot(1,3,2);\n",
    "plt.axis( [left, right, bottom, top] );\n",
    "\n",
    "# plot one rectangle per sample i, with color scaled according to\n",
    "# the strength of the i-th element in factor nf=2\n",
    "nf = 2;\n",
    "FPmean = np.mean(FP[nf,:]);\n",
    "FPstd = np.std(FP[nf,:]);\n",
    "fmin = -2*FPstd\n",
    "fmax = 2*FPstd\n",
    "for i in range(N):\n",
    "    h = 3;\n",
    "    c = (FP[nf,i]-fmin)/(fmax-fmin);\n",
    "    if( c<0.0 ):\n",
    "        c=0.0;\n",
    "    elif ( c>1.0 ):\n",
    "        c=1.0;\n",
    "    ax2.add_patch( pch.Rectangle( (lon[i,0], lat[i,0]), h, h, fc=(c,c,c), ec=(0.0,0.0,0.0) ) );\n",
    "\n",
    "plt.xlabel('lonomgitude (deg)');\n",
    "plt.ylabel('latitude (deg)');\n",
    "plt.title('Factor 2');\n",
    "\n",
    "# plot coastlines on maps, using coarse coastlines from the file coastdata.txt\n",
    "# this file consists of lists of (lat,lons) along segments of coastlines, with\n",
    "# each list deliited with a \">\".\n",
    "coastfile = 'coastdata.txt';\n",
    "fd = open(coastfile,\"r\");\n",
    "NN = 0;\n",
    "NNMAX = 10000;\n",
    "mylat = np.zeros((NNMAX,1));\n",
    "mylon = np.zeros((NNMAX,1));\n",
    "for i in range(100000):\n",
    "    line = fd.readline();\n",
    "    if not line:\n",
    "        break;\n",
    "    if( line[0] == '>' ):\n",
    "        if( NN == 0):\n",
    "            continue;\n",
    "        plt.plot( mylon[0:NN,0], mylat[0:NN,0], 'k.', color=(0.7,0.7,0.7) );\n",
    "        NN=0;\n",
    "    else:\n",
    "        latlonstr=line.split();\n",
    "        mylat[NN,0] = float(latlonstr[1]);\n",
    "        mylon[NN,0] = float(latlonstr[0]);\n",
    "        NN=NN+1;\n",
    "fd.close();\n",
    "\n",
    "# PART 3: factor 3\n",
    "ax3 = plt.subplot(1,3,3);\n",
    "plt.axis( [left, right, bottom, top] );\n",
    "\n",
    "# plot one rectangle per sample i, with color scaled according to\n",
    "# the strength of the i-th element in factor nf=3\n",
    "nf = 3;\n",
    "FPmean = np.mean(FP[nf,:]);\n",
    "FPstd = np.std(FP[nf,:]);\n",
    "fmin = -2*FPstd\n",
    "fmax = 2*FPstd\n",
    "for i in range(N):\n",
    "    h = 3;\n",
    "    c = (FP[nf,i]-fmin)/(fmax-fmin);\n",
    "    if( c<0.0 ):\n",
    "        c=0.0;\n",
    "    elif ( c>1.0 ):\n",
    "        c=1.0;\n",
    "    ax3.add_patch( pch.Rectangle( (lon[i,0], lat[i,0]), h, h, fc=(c,c,c), ec=(0.0,0.0,0.0) ) );\n",
    "\n",
    "plt.xlabel('lonomgitude (deg)');\n",
    "plt.ylabel('latitude (deg)');\n",
    "plt.title('Factor 3');\n",
    "\n",
    "# plot coastlines on maps, using coarse coastlines from the file coastdata.txt\n",
    "# this file consists of lists of (lat,lons) along segments of coastlines, with\n",
    "# each list deliited with a \">\".\n",
    "coastfile = 'coastdata.txt';\n",
    "fd = open(coastfile,\"r\");\n",
    "NN = 0;\n",
    "NNMAX = 10000;\n",
    "mylat = np.zeros((NNMAX,1));\n",
    "mylon = np.zeros((NNMAX,1));\n",
    "for i in range(100000):\n",
    "    line = fd.readline();\n",
    "    if not line:\n",
    "        break;\n",
    "    if( line[0] == '>' ):\n",
    "        if( NN == 0):\n",
    "            continue;\n",
    "        plt.plot( mylon[0:NN,0], mylat[0:NN,0], 'k.', color=(0.7,0.7,0.7) );\n",
    "        NN=0;\n",
    "    else:\n",
    "        latlonstr=line.split();\n",
    "        mylat[NN,0] = float(latlonstr[1]);\n",
    "        mylon[NN,0] = float(latlonstr[0]);\n",
    "        NN=NN+1;\n",
    "fd.close();\n",
    "\n",
    "\n",
    "plt.show();\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Images: 399\n",
      "start year 1970 month 1\n",
      "end   year 2003 month 3\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 327192 characters]\n",
      "text/plain": [
       "<Figure size 1152x720 with 72 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 350092 characters]\n",
      "text/plain": [
       "<Figure size 1152x720 with 72 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 338992 characters]\n",
      "text/plain": [
       "<Figure size 1152x720 with 72 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 336632 characters]\n",
      "text/plain": [
       "<Figure size 1152x720 with 72 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 333040 characters]\n",
      "text/plain": [
       "<Figure size 1152x720 with 72 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 22056 characters]\n",
      "text/plain": [
       "<Figure size 1152x720 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda08_09\n",
    "# CAC Pacific Sea Surface Temperature (SST) dataset\n",
    "# read in and plot the data\n",
    "\n",
    "# load the data\n",
    "D = np.genfromtxt('cac_sst.txt', delimiter='\\t')\n",
    "\n",
    "# set up sizes of arrays\n",
    "I, J = np.shape(D);\n",
    "NBLOCK=85;\n",
    "NLON = NBLOCK-1;\n",
    "NLAT = J-1;\n",
    "IMAGES = floor(I/NBLOCK);\n",
    "Dt = 1/12;\n",
    "print('Images:',IMAGES)\n",
    "\n",
    "# arrays for SST, date, lats, lons\n",
    "SST = np.zeros((NLON,NLAT,IMAGES));\n",
    "thedate = np.zeros((IMAGES,1)); # in month.year format\n",
    "theyear = np.zeros((IMAGES,1)); # year extracted from thedate\n",
    "themonth = np.zeros((IMAGES,1)); # month extracted from thedate\n",
    "lats=np.zeros((NLAT,1));\n",
    "lons=np.zeros((NLON,1));\n",
    "\n",
    "# cut up data into SST arrays\n",
    "lats = np.copy(D[0,1:J]).T;\n",
    "lons = np.copy(D[1:NBLOCK,0]);\n",
    "for j in range(IMAGES):\n",
    "    k1 = j*NBLOCK;\n",
    "    k2 = k1+NLON;\n",
    "    thedate[j,0]=D[k1,0];\n",
    "    themonth[j,0]=floor(thedate[j,0]);\n",
    "    theyear[j,0]=floor(10000*(thedate[j,0]-themonth[j,0])+0.1);\n",
    "    SST[:,:,j] = np.copy(D[k1+1:k2+1,1:J]);\n",
    "    \n",
    "print(\"start year %d month %d\" % (theyear[0,0], themonth[0,0]));\n",
    "print(\"end   year %d month %d\" % (theyear[IMAGES-1,0], themonth[IMAGES-1,0]));\n",
    "    \n",
    "# some definitions related to plotting\n",
    "MONTHS=12;\n",
    "YEARS=ceil(IMAGES/12);\n",
    "YBLOCKSIZE = 6;\n",
    "YBLOCKS = ceil(YEARS/YBLOCKSIZE);\n",
    "\n",
    "# grey-scale colormap\n",
    "bw = np.zeros((256,4));\n",
    "v = (255 - np.linspace( 0, 255, 256 ))/255;\n",
    "bw[:,0] = v;\n",
    "bw[:,1] = v;\n",
    "bw[:,2] = v;\n",
    "bw[:,3] = np.ones(256);\n",
    "bwcmap = ListedColormap(bw);\n",
    "\n",
    "# plot data, 12 months by six years of maps per figure\n",
    "lastmonth = (IMAGES % YBLOCKSIZE*12) % YBLOCKSIZE ;\n",
    "done=0;\n",
    "for yb in range(YBLOCKS):\n",
    "    fig = plt.figure(yb+2,figsize=(16,10));\n",
    "    for m in range(MONTHS):\n",
    "        for y in range(YBLOCKSIZE):\n",
    "            j = m+MONTHS*y+MONTHS*YBLOCKSIZE*yb;\n",
    "            k = y + m*YBLOCKSIZE;\n",
    "            if( j>=IMAGES):\n",
    "                done=1;\n",
    "                break;\n",
    "            ax = plt.subplot(MONTHS, YBLOCKSIZE, k+1);\n",
    "            left=0;\n",
    "            right=3.0;\n",
    "            bottom=0.0;\n",
    "            top=1.0;\n",
    "            plt.axis( [left, right, bottom, top] );\n",
    "            ax.xaxis.set_ticks([]);\n",
    "            ax.yaxis.set_ticks([]);\n",
    "            MAP=np.flipud(SST[:,:,j].T); # reorient so it plots as map\n",
    "            MAPmin = np.min(MAP);\n",
    "            MAPmax = np.max(MAP);\n",
    "            MAPrange = MAPmax-MAPmin;\n",
    "            if( MAPrange<=0.0 ):\n",
    "                MAPrange=1;\n",
    "            plt.imshow( (MAP-MAPmin)/MAPrange, cmap=bwcmap, vmin=0, vmax=1, extent=(left,right,bottom,top) );\n",
    "            if( y==0 ):\n",
    "                plt.ylabel(\"%d\" % (themonth[j,0]) );\n",
    "            if( (m==11) or ((yb==(YBLOCKS-1)) and (m==lastmonth))  ):\n",
    "                plt.xlabel(\"%d\" % (theyear[j,0]) );\n",
    "        if( done==1 ):\n",
    "            break;\n",
    "    if( done==1 ):\n",
    "            break;\n",
    "    plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 14788 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#  eda08_10\n",
    "# CAC Pacific Sea Surface Temperature (SST) dataset\n",
    "# plot singular values\n",
    "\n",
    "# load the data\n",
    "D = np.genfromtxt('cac_sst.txt', delimiter='\\t')\n",
    "\n",
    "# set up sizes of arrays\n",
    "I, J = np.shape(D);\n",
    "NBLOCK=85;\n",
    "NLON = NBLOCK-1;\n",
    "NLAT = J-1;\n",
    "IMAGES = floor(I/NBLOCK);\n",
    "Dt = 1/12;\n",
    "\n",
    "# arrays for SST, date, lats, lons\n",
    "SST = np.zeros((NLON,NLAT,IMAGES));\n",
    "thedate = np.zeros((IMAGES,1)); # in month.year format\n",
    "theyear = np.zeros((IMAGES,1)); # year extracted from thedate\n",
    "themonth = np.zeros((IMAGES,1)); # month extracted from thedate\n",
    "lats=np.zeros((NLAT,1));\n",
    "lons=np.zeros((NLON,1));\n",
    "\n",
    "# cut up data into SST arrays\n",
    "lats = np.copy(D[0,1:J]).T;\n",
    "lons = np.copy(D[1:NBLOCK,0]);\n",
    "for j in range(IMAGES):\n",
    "    k1 = j*NBLOCK;\n",
    "    k2 = k1+NLON;\n",
    "    thedate[j,0]=D[k1,0];\n",
    "    themonth[j,0]=floor(thedate[j,0]);\n",
    "    theyear[j,0]=floor(10000*(thedate[j,0]-themonth[j,0])+0.1);\n",
    "    SST[:,:,j] = np.copy(D[k1+1:k2+1,1:J]);\n",
    "    \n",
    "# some definitions related to plotting\n",
    "MONTHS=12;\n",
    "YEARS=ceil(IMAGES/12);\n",
    "YBLOCKSIZE = 6;\n",
    "YBLOCKS = ceil(YEARS/YBLOCKSIZE);\n",
    "\n",
    "# fold out the images into the sample array\n",
    "SSTV = np.zeros( (IMAGES, NLAT*NLON) );\n",
    "for j in range(IMAGES):\n",
    "    for p in range(NLAT):\n",
    "        for q in range(NLON):\n",
    "            k = p + q*NLAT;\n",
    "            SSTV[j,k] = SST[q, p, j];\n",
    "\n",
    "# singular value decomposition\n",
    "[U, sigma, VT] = la.svd(SSTV,full_matrices=False); # singular values sigma\n",
    "sh = np.shape(sigma);\n",
    "Ns = sh[0];\n",
    "F = np.copy(VT); # factors\n",
    "C = np.matmul( U, np.diag(sigma) ); # factor loadings\n",
    "\n",
    "# plot singular values\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, Ns+1, 0, 1.1*np.max(sigma)] );\n",
    "plt.plot( np.linspace(1,Ns,Ns), sigma, 'k-');\n",
    "plt.title('Singular Values of the CAC dataset');\n",
    "plt.xlabel('index i');\n",
    "plt.ylabel('singular value i');\n",
    "\n",
    "plt.show();\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 116496 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 12 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 132248 characters]\n",
      "text/plain": [
       "<Figure size 1152x864 with 12 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda08_11\n",
    "# CAC Pacific Sea Surface Temperature (SST) dataset\n",
    "# plot factors as maps and loadings as timeseries\n",
    "\n",
    "# load the data\n",
    "D = np.genfromtxt('cac_sst.txt', delimiter='\\t')\n",
    "\n",
    "# set up sizes of arrays\n",
    "I, J = np.shape(D);\n",
    "NBLOCK=85;\n",
    "NLON = NBLOCK-1;\n",
    "NLAT = J-1;\n",
    "IMAGES = floor(I/NBLOCK);\n",
    "Dt = 1/12;\n",
    "\n",
    "# arrays for SST, date, lats, lons\n",
    "SST = np.zeros((NLON,NLAT,IMAGES));\n",
    "thedate = np.zeros((IMAGES,1)); # in month.year format\n",
    "theyear = np.zeros((IMAGES,1)); # year extracted from thedate\n",
    "themonth = np.zeros((IMAGES,1)); # month extracted from thedate\n",
    "lats=np.zeros((NLAT,1));\n",
    "lons=np.zeros((NLON,1));\n",
    "\n",
    "# cut up data into SST arrays\n",
    "lats = np.copy(D[0,1:J]).T;\n",
    "lons = np.copy(D[1:NBLOCK,0]);\n",
    "for j in range(IMAGES):\n",
    "    k1 = j*NBLOCK;\n",
    "    k2 = k1+NLON;\n",
    "    thedate[j,0]=D[k1,0];\n",
    "    themonth[j,0]=floor(thedate[j,0]);\n",
    "    theyear[j,0]=floor(10000*(thedate[j,0]-themonth[j,0])+0.1);\n",
    "    SST[:,:,j] = np.copy(D[k1+1:k2+1,1:J]);\n",
    "    \n",
    "    \n",
    "# some definitions related to plotting\n",
    "MONTHS=12;\n",
    "YEARS=ceil(IMAGES/12);\n",
    "YBLOCKSIZE = 6;\n",
    "YBLOCKS = ceil(YEARS/YBLOCKSIZE);\n",
    "\n",
    "# fold out the images into the sample matrix\n",
    "SSTV = np.zeros( (IMAGES, NLAT*NLON) );\n",
    "for j in range(IMAGES):\n",
    "    for p in range(NLAT):\n",
    "        for q in range(NLON):\n",
    "            k = p + q*NLAT;\n",
    "            SSTV[j,k] = SST[q, p, j];\n",
    "\n",
    "# singular value decomposition\n",
    "[U, sigma, VT] = la.svd(SSTV,full_matrices=False);\n",
    "sh = np.shape(sigma); # sigma is the singular values\n",
    "Ns = sh[0];\n",
    "\n",
    "# keep only P singular values\n",
    "P=Ns;\n",
    "F = np.copy(VT[0:P,:]); # factors\n",
    "C = np.matmul( U[:,0:P], np.diag(sigma[0:P])); # factor loadings\n",
    "\n",
    "# fold EOF's back into images (maps)\n",
    "Nf = 12;  # just do a few EOF's\n",
    "SSTF = np.zeros((NLON,NLAT,Nf));\n",
    "for j in range(Nf):\n",
    "    for p in range(NLAT):\n",
    "        for q in range(NLON):\n",
    "            k = p + q*NLAT;\n",
    "            SSTF[ q, p, j ] = F[j,k];\n",
    "            \n",
    "# grey-scale colormap\n",
    "bw = np.zeros((256,4));\n",
    "v = (255 - np.linspace( 0, 255, 256 ))/255;\n",
    "bw[:,0] = v;\n",
    "bw[:,1] = v;\n",
    "bw[:,2] = v;\n",
    "bw[:,3] = np.ones(256);\n",
    "bwcmap = ListedColormap(bw);\n",
    "\n",
    "# plot these EOF's\n",
    "fig1 = plt.figure(1,figsize=(16,6));\n",
    "for f in range(Nf):\n",
    "    ax = plt.subplot(3, 4, f+1);\n",
    "    left=0;\n",
    "    right=3.0;\n",
    "    bottom=0.0;\n",
    "    top=1.0;\n",
    "    plt.axis( [left, right, bottom, top] );\n",
    "    ax.xaxis.set_ticks([]);\n",
    "    ax.yaxis.set_ticks([]);\n",
    "    MAP=np.flipud(SSTF[:,:,f].T); # reorient so it plots as map\n",
    "    MAPmin = np.min(MAP);\n",
    "    MAPmax = np.max(MAP);\n",
    "    MAPrange = MAPmax-MAPmin;\n",
    "    if( MAPrange<=0.0 ):\n",
    "        MAPrange=1;\n",
    "    plt.imshow( (MAP-MAPmin)/MAPrange, cmap=bwcmap, vmin=0, vmax=1, extent=(left,right,bottom,top) );\n",
    "    plt.title(\"EOF %d\" % (f+1) );\n",
    "\n",
    "plt.show();\n",
    "\n",
    "# plot corresponding amplitude timeseries\n",
    "fig2 = plt.figure(2,figsize=(16,12));\n",
    "for f in range(Nf):\n",
    "    ax = plt.subplot(3, 4, f+1);\n",
    "    plt.axis( [theyear[0,0], theyear[IMAGES-1,0]+1, np.min(C[:,0]), np.max(C[:,0]) ] );\n",
    "    plt.plot( theyear[0,0]+np.linspace(1,IMAGES,IMAGES)/12, C[:,f], 'k-' );\n",
    "    plt.title(\"C(%d,t)\" % (f) );\n",
    "\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Images: 399\n",
      "start year 1970 month 1\n",
      "end   year 2003 month 3\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 290608 characters]\n",
      "text/plain": [
       "<Figure size 1152x720 with 72 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 293700 characters]\n",
      "text/plain": [
       "<Figure size 1152x720 with 72 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda08_12\n",
    "# CAC Pacific Sea Surface Temperature (SST) dataset\n",
    "# plot approximation to SST's using just lowest factors\n",
    "\n",
    "# load the data\n",
    "D = np.genfromtxt('cac_sst.txt', delimiter='\\t')\n",
    "\n",
    "# set up sizes of arrays\n",
    "I, J = np.shape(D);\n",
    "NBLOCK=85;\n",
    "NLON = NBLOCK-1;\n",
    "NLAT = J-1;\n",
    "IMAGES = floor(I/NBLOCK);\n",
    "Dt = 1/12;\n",
    "print('Images:',IMAGES)\n",
    "\n",
    "# arrays for SST, date, lats, lons\n",
    "SST = np.zeros((NLON,NLAT,IMAGES));\n",
    "thedate = np.zeros((IMAGES,1)); # in month.year format\n",
    "theyear = np.zeros((IMAGES,1)); # year extracted from thedate\n",
    "themonth = np.zeros((IMAGES,1)); # month extracted from thedate\n",
    "lats=np.zeros((NLAT,1));\n",
    "lons=np.zeros((NLON,1));\n",
    "\n",
    "# cut up data into SST arrays\n",
    "lats = np.copy(D[0,1:J]).T;\n",
    "lons = np.copy(D[1:NBLOCK,0]);\n",
    "for j in range(IMAGES):\n",
    "    k1 = j*NBLOCK;\n",
    "    k2 = k1+NLON;\n",
    "    thedate[j,0]=D[k1,0];\n",
    "    themonth[j,0]=floor(thedate[j,0]);\n",
    "    theyear[j,0]=floor(10000*(thedate[j,0]-themonth[j,0])+0.1);\n",
    "    SST[:,:,j] = np.copy(D[k1+1:k2+1,1:J]);\n",
    "    \n",
    "print(\"start year %d month %d\" % (theyear[0,0], themonth[0,0]));\n",
    "print(\"end   year %d month %d\" % (theyear[IMAGES-1,0], themonth[IMAGES-1,0]));\n",
    "    \n",
    "# some definitions related to plotting\n",
    "MONTHS=12;\n",
    "YEARS=ceil(IMAGES/12);\n",
    "YBLOCKSIZE = 6;\n",
    "YBLOCKS = ceil(YEARS/YBLOCKSIZE);\n",
    "\n",
    "# fold out images into the sample matrix\n",
    "SSTV = np.zeros( (IMAGES, NLAT*NLON) );\n",
    "for j in range(IMAGES):\n",
    "    for p in range(NLAT):\n",
    "        for q in range(NLON):\n",
    "            k = p + q*NLAT;\n",
    "            SSTV[j,k] = SST[q, p, j];\n",
    "\n",
    "# singular value decomposition\n",
    "[U, sigma, VT] = la.svd(SSTV,full_matrices=False);\n",
    "sh = np.shape(sigma); # sigma is the singular values\n",
    "Ns = sh[0];\n",
    "\n",
    "# keep only first few singular values\n",
    "P=5;\n",
    "FP = np.copy(VT[0:P,:]); # factors\n",
    "CP = np.matmul( U[:,0:P], np.diag(sigma[0:P])); # factor loadings\n",
    "\n",
    "# reconstruct with a small number of EOFs\n",
    "SSTVp = np.matmul(CP, FP);\n",
    "\n",
    "# fold reconstructed sample matrix back into images (maps)\n",
    "Nf = IMAGES;\n",
    "SSTp = np.zeros((NLON,NLAT,Nf));\n",
    "for j in range(Nf):\n",
    "    for p in range(NLAT):\n",
    "        for q in range(NLON):\n",
    "            k = p + q*NLAT;\n",
    "            SSTp[ q, p, j ] = SSTVp[j,k];\n",
    "\n",
    "# gre-scale colormap\n",
    "bw = np.zeros((256,4));\n",
    "v = (255 - np.linspace( 0, 255, 256 ))/255;\n",
    "bw[:,0] = v;\n",
    "bw[:,1] = v;\n",
    "bw[:,2] = v;\n",
    "bw[:,3] = np.ones(256);\n",
    "bwcmap = ListedColormap(bw);\n",
    "\n",
    "# plor data, 12 months by six years of maps per figure\n",
    "lastmonth = (IMAGES % YBLOCKSIZE*12) % YBLOCKSIZE ;\n",
    "done=0;\n",
    "for yb in range(YBLOCKS):\n",
    "    fig = plt.figure(yb+2,figsize=(16,10));\n",
    "    for m in range(MONTHS):\n",
    "        for y in range(YBLOCKSIZE):\n",
    "            j = m+MONTHS*y+MONTHS*YBLOCKSIZE*yb;\n",
    "            k = y + m*YBLOCKSIZE;\n",
    "            if( j>=IMAGES):\n",
    "                done=1;\n",
    "                break;\n",
    "            ax = plt.subplot(MONTHS, YBLOCKSIZE, k+1);\n",
    "            left=0;\n",
    "            right=3.0;\n",
    "            bottom=0.0;\n",
    "            top=1.0;\n",
    "            plt.axis( [left, right, bottom, top] );\n",
    "            ax.xaxis.set_ticks([]);\n",
    "            ax.yaxis.set_ticks([]);\n",
    "            MAP=np.flipud(SSTp[:,:,j].T); # reorient so it plots as map\n",
    "            MAPmin = np.min(MAP);\n",
    "            MAPmax = np.max(MAP);\n",
    "            MAPrange = MAPmax-MAPmin;\n",
    "            if( MAPrange<=0.0 ):\n",
    "                MAPrange=1;\n",
    "            plt.imshow( (MAP-MAPmin)/MAPrange, cmap=bwcmap, vmin=0, vmax=1, extent=(left,right,bottom,top) );\n",
    "            if( y==0 ):\n",
    "                plt.ylabel(\"%d\" % (themonth[j,0]) );\n",
    "            if( (m==11) or ((yb==(YBLOCKS-1)) and (m==lastmonth))  ):\n",
    "                plt.xlabel(\"%d\" % (theyear[j,0]) );\n",
    "        if( done==1 ):\n",
    "            break;\n",
    "    if( done==1 ):\n",
    "            break;\n",
    "    plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
