{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# edapy_08_00 clear all variables and import various 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( N==0 ):\n",
    "            continue; # skip vector of zero length\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",
    "\n",
    "def eda_kmeans(S, mu):\n",
    "    # k-mean cluster analsis of NxM sample natrix S\n",
    "    # starting from initial guess KxM mu of cluster means\n",
    "    # returns:\n",
    "    # newmu: new cluster means\n",
    "    # k: Nx1 vector of cluster assignments\n",
    "    # dist2: Nx1 vector of squared distance to cluster center\n",
    "    # count: Kx1 vector of samples in cluster\n",
    "    # Kdist2: KxK matrix of inter-cluster squared distances\n",
    "    N, M = np.shape(S);   # N samples with M elements\n",
    "    K, M2 = np.shape(mu); # K cluster meabs\n",
    "    k = np.zeros((N,1),dtype=np.int); # assoc sample w cluster\n",
    "    dist2 = np.zeros((N,1)); # sq distance from cluster mean\n",
    "    Kdist2 = np.zeros((K,K)); # cluster to nearest cluster sq distance\n",
    "    oldmu = np.copy(mu);  # old cluster means\n",
    "    \n",
    "    MAXP = 2*K*M;  # max iterations\n",
    "    for p in range(MAXP):  # iterate\n",
    "    \n",
    "        # cluster assignment step\n",
    "        changes = 0; # number of changes\n",
    "        for i in range(N):     # loop over samples\n",
    "            for j in range(K): # loop over cluster means\n",
    "                # deviation of sample from mean\n",
    "                delta = (S[i:i+1,0:M] - oldmu[j:j+1,0:M]).T;\n",
    "                # squared distance\n",
    "                r2 = np.matmul(delta.T,delta);\n",
    "                if( j==0 ): # always accept first distance\n",
    "                    r2min = r2;\n",
    "                    newk = j;\n",
    "                elif r2<r2min: # accept if smaller\n",
    "                    r2min = r2;\n",
    "                    newk = j;\n",
    "            if k[i,0] != newk: # count reassignements\n",
    "                changes=changes+1;\n",
    "            k[i,0] = newk;     # reassigned cluster\n",
    "            dist2[i,0] = r2min; # sq distance to that cluster\n",
    "\n",
    "        # end loop over iterations on no further changes\n",
    "        if( (p!=0) and (changes==0) ):\n",
    "            break;\n",
    "\n",
    "        # means step\n",
    "        newmu = np.zeros((K,M));  # new cluster means\n",
    "        count = np.zeros((K,1),dtype=np.int); # samples in cluster\n",
    "        for i in range(N): # loop over samples\n",
    "            for j in range(M): # loop over elements\n",
    "                # cluster sum\n",
    "                newmu[k[i,0],j] = newmu[k[i,0],j] + S[i,j];\n",
    "            # count samples in cluster\n",
    "            count[k[i,0],0] = count[k[i,0],0] + 1;\n",
    "        for j in range(K): # loop over clusters\n",
    "            if( count[j,0] == 0 ): # cluser with no samples\n",
    "                # no change in cluster mean\n",
    "                newmu[j,0:M] = oldmu[j,0:M];\n",
    "            else:\n",
    "                # new cluster mean\n",
    "                newmu[j,0:M] = newmu[j,0:M]/count[j,0]\n",
    "        oldmu=newmu; # update cluster mean\n",
    "        \n",
    "    # inter-cluster distances\n",
    "    for i in range(K-1):     # loop over cluster means\n",
    "        for j in range(1,K): # loop over cluster means\n",
    "            delta = (newmu[i:i+1,0:M] - newmu[j:j+1,0:M]).T;\n",
    "            r2 = np.matmul(delta.T,delta);\n",
    "            Kdist2[i,j] = r2; # sq distance between clusters\n",
    "            Kdist2[j,i] = r2; # sq distance between clusters\n",
    "    return newmu, k, dist2, count, Kdist2;\n",
    "\n",
    "def eda_kmeans_mod(S, mu):\n",
    "    # k-mean cluster analsis of NxM sample natrix S\n",
    "    # modified per Bradley & Fayyad (1998) never to\n",
    "    # return cluster with zero samples.\n",
    "    # Starts with initial guess KxM mu of cluster means\n",
    "    # returns:\n",
    "    # mymu2: new cluster means\n",
    "    # myk2: Nx1 vector of cluster assignments\n",
    "    # mydist2: Nx1 vector of distance to cluster center\n",
    "    # mycount: Kx1 vector of samples in cluster\n",
    "    initmu = np.copy(mu);\n",
    "    mymu, myk, mydist, mycount, myKdist = eda_kmeans(S, initmu);\n",
    "    while np.min(mycount)==0:  # while cluster with no samples\n",
    "        ic=np.argmin(mycount); # index of cluster\n",
    "        ir=np.argmax(mydist);  # index of most distant sample\n",
    "        ik=myk[ir,0];          # cluster of this sample\n",
    "        mymu[ic:ic+1,0:M] = S[ir:ir+1,0:M]; # new mean\n",
    "        mycount[ic,0]=1;       # reset cluster count\n",
    "        myk[ir,0]=ic;          # reset sample membership\n",
    "        mydist[ir,0]=0.0;      # reset sample distance\n",
    "        # note: in rare cases, next line could create\n",
    "        # infinite loop, so don't do it.\n",
    "        # mycount[ik,0]=mycount[ik,0]-1;\n",
    "    # redo k-means\n",
    "    mymu2, myk2, mydist2, mycount2, myKdist2 = eda_kmeans(S, mymu);  \n",
    "    return mymu2, myk2, mydist2, mycount2, myKdist2;\n",
    "\n",
    "def eda_forgy_mean(S,K):\n",
    "    # create initial k-mean cluster means, from samples\n",
    "    # randomly drawn NxM sample matrix S\n",
    "    # returns KxM matrix mu of means\n",
    "    N,M=np.shape(S);\n",
    "    k = np.random.randint(low=0,high=N,size=K,dtype=np.int);\n",
    "    mu = S[k,0:M];\n",
    "    return mu;\n",
    "\n",
    "def eda_random_partition_mean(S,K):\n",
    "    # create initial k-mean cluster means, by K random\n",
    "    # partitions of NxM sample matrix S\n",
    "    # returns KxM matrix mu of means\n",
    "    N,M=np.shape(S);\n",
    "    # random cluster-assignments\n",
    "    k = np.random.randint(low=0,high=K,size=N,dtype=np.int);\n",
    "    mu = np.zeros((K,M));\n",
    "    for i in range(K): # loop over clusters\n",
    "        r = np.where(k==i); # samples in cluster i\n",
    "        kk = r[0];\n",
    "        mu[i:i+1,0:M] = np.mean(S[kk,0:M],axis=0); # their means\n",
    "    return mu;\n",
    "\n",
    "def eda_bradley_fayyad_mean(S,K,Nsets,fraction):\n",
    "    # create initial k-mean cluster means, by\n",
    "    # the Bradley & Fayyad (1998) method, which\n",
    "    # is based on clustering clusters. The NxM\n",
    "    # sample matrix S is randomly subsampled Nsets\n",
    "    # times, with each set containing (fraction*N+2*K)\n",
    "    # samples.\n",
    "    # returns KxM matrix mu of means\n",
    "    \n",
    "    N,M=np.shape(S);\n",
    "    initmu = eda_random_partition_mean(S,K);\n",
    "\n",
    "    KiMAX = K*Nsets;  # maximum number of means mui\n",
    "    Ni = floor(fraction*N)+2*K; # size of small sets\n",
    "\n",
    "    # Step 1, create KixM table CM of means of subsampled data\n",
    "    CM = np.zeros((KiMAX,M));\n",
    "    Ki=0; # counter for rows of CM\n",
    "    for j in range(Nsets):\n",
    "        k = np.random.randint(low=0,high=N,size=Ni,dtype=np.int);\n",
    "        Si = S[k,0:M];\n",
    "        mymu, myk, mydist, mycount, myKdist = eda_kmeans_mod(Si, initmu);\n",
    "        CM[Ki:Ki+K,0:M] = mymu[0:K,0:M];\n",
    "        Ki=Ki+K;\n",
    "\n",
    "    # Step 2, create KixM table FMS of means of cluster means\n",
    "    FMS = np.zeros((KiMAX,M));\n",
    "    Ki=0; # counter for rows of CM\n",
    "    for j in range(Nsets):\n",
    "        mymu, myk, mydist, mycount, myKdist = eda_kmeans(CM, CM[3*j:3*j+K,0:M]);\n",
    "        FMS[Ki:Ki+K,0:M] = mymu[0:K,0:M];\n",
    "        Ki=Ki+K;\n",
    "\n",
    "    # Step 3A: associate each KxM member FMi of Fm\n",
    "    #          with all individual 1xM elements of CM\n",
    "    #          and compute the overall error\n",
    "    error = np.zeros((Nsets,1));\n",
    "    Ki=0;\n",
    "    for i in range(Nsets):\n",
    "        FMi = FMS[Ki:Ki+K,0:M] # one group of means\n",
    "        Ki=Ki+K;\n",
    "        for j in range(KiMAX):\n",
    "            CMj = CM[j:j+1,0:M]; # one member of CM\n",
    "            for m in range(K):\n",
    "                delta = (CMj[0:1,0:M] - FMi[m:m+1,0:M]).T;\n",
    "                r2 = np.matmul(delta.T,delta);\n",
    "                if( m==0 ): # always acceot first distance\n",
    "                    r2min = r2;\n",
    "                elif r2<r2min: # accept if smaller\n",
    "                    r2min = r2;\n",
    "        error[i,0] = error[i,0] + r2min;\n",
    "\n",
    "    # Step 3B: select the FMi with lowest error\n",
    "    i = np.argmin(error);\n",
    "    mu = np.copy( FMS[K*i:K*i+K,0:M] );\n",
    "    return(mu);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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": [
    "# edaps_08_01, 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": 3,
   "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": [
    "# edapy_08_02, 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": 4,
   "metadata": {
    "scrolled": true
   },
   "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": [
    "# edapy_08_03, alternative 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": 5,
   "metadata": {
    "scrolled": true
   },
   "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": [
    "# edapy_08_04, factor analysis of Atlantic Rocks dataset\n",
    "# using singular value decomposition\n",
    "\n",
    "# load data\n",
    "S = np.genfromtxt('../Data/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": 6,
   "metadata": {
    "scrolled": true
   },
   "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": [
    "# edapy_08_05: varimax procedure applied 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",
    "A = 2*M*np.matmul(u.T,v);\n",
    "B = np.sum(u)*np.sum(v);\n",
    "top = A - B;\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",
    "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": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max relative error in S-CF   0.05663026082140537\n",
      "max relative error in S-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('../Data/rocks.txt', delimiter='\\t')\n",
    "N, M = np.shape(S);\n",
    "Smax = np.max(np.abs(S));\n",
    "\n",
    "# compute factors and factor loadings using singular value decompostion\n",
    "[U, sigma, VT] = la.svd(S,full_matrices=False);\n",
    "\n",
    "# keep only first five singular values\n",
    "P=5;\n",
    "F = np.copy(VT[0:P,:]);  # P factors\n",
    "C = np.matmul( U[:,0:P], np.diag(sigma[0:P])); # P loadings\n",
    "SP = np.matmul(C,F); # reconstruction of S with P factors\n",
    "\n",
    "\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));   # fator A of a pair\n",
    "fB = np.zeros((M,1));   # fator B of a pair\n",
    "fAp = np.zeros((M,1));  # fator A of a pair after rotation\n",
    "fBp = np.zeros((M,1));  # fator B of a pair after rotation\n",
    "mA = np.zeros((P,1));   # col of rotation matrix associated with A\n",
    "mB = np.zeros((P,1));   # col of rotation matrix associated with B\n",
    "mAp = np.zeros((P,1));  # col A of rotation matrix after rotation\n",
    "mBp = np.zeros((P,1));  # col B of rotation matrix after rotation\n",
    "\n",
    "# spike these factors using the varimax procedure\n",
    "k = [1, 2, 3, 4];\n",
    "Nk = len(k);\n",
    "\n",
    "for iter in range(3):  # three iterations almost always sufficient\n",
    "    for ii in range(Nk): # first factor of pair\n",
    "        for jj in range(ii+1,Nk): # second factor of pair\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:i+1,0:M] = fAp[0:M,0];\n",
    "            FP[j:j+1,0:M] = fBp[0:M,0];\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:i+1,0:P] = mAp[0:P,0];\n",
    "            MR[j:j+1,0:P] = mBp[0:P,0];\n",
    "\n",
    "# new factor loadings\n",
    "CP = np.matmul(C,MR.T);\n",
    "\n",
    "# reconstruction using new factors\n",
    "SPP = np.matmul(CP,FP);\n",
    "\n",
    "# check error of reconstruction\n",
    "Smax = np.max(np.abs(S));\n",
    "Emax1 = np.max(np.abs(S - SP));\n",
    "Emax2 = np.max(np.abs(S - SPP));\n",
    "print('max relative error in S-CF  ', Emax1/Smax);\n",
    "print('max relative error in S-CPFP', Emax2/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": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "witout weighting, max relative error in CF   0.18914537524554229\n",
      "with   weighting, max relative error in CF   0.18145117130234192\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 23492 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 17500 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": 12,
   "metadata": {},
   "outputs": [
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: 'coastdata.txt'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-12-9c9a667c3bd2>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m    156\u001b[0m \u001b[1;31m# each list deliited with a \">\".\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    157\u001b[0m \u001b[0mcoastfile\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'coastdata.txt'\u001b[0m\u001b[1;33m;\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 158\u001b[1;33m \u001b[0mfd\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mopen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcoastfile\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m\"r\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m;\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    159\u001b[0m \u001b[0mNN\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m;\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    160\u001b[0m \u001b[0mNNMAX\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m10000\u001b[0m\u001b[1;33m;\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'coastdata.txt'"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 81700 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('../Data/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 = '../Data/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 = '../Data/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": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Images: 399\n",
      "start year 1970 month 1\n",
      "end   year 2003 month 3\n",
      "6 blocks of years\n",
      "lastmonth 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, 193392 characters]\n",
      "text/plain": [
       "<Figure size 1152x720 with 39 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('../Data/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;                         # months in year\n",
    "YEARS=ceil(IMAGES/12);             # years of images\n",
    "YBLOCKSIZE = 6;                    # 6 years per plot block of years\n",
    "YBLOCKS = ceil(YEARS/YBLOCKSIZE);  # number of blocks\n",
    "\n",
    "print(\"%d blocks of years\" % (YBLOCKS) );\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",
    "print(\"lastmonth %d\" % (lastmonth));\n",
    "done=0;\n",
    "for yb in range(YBLOCKS):\n",
    "    fig = plt.figure(yb+2,figsize=(16,10));\n",
    "    for y in range(YBLOCKSIZE):\n",
    "        for m in range(MONTHS):\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-1)))  ):\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": 14,
   "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('../Data/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": 15,
   "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('../Data/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": 3,
   "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"
    },
    {
     "data": {
      "image/png": "[encoded data, 285432 characters]\n",
      "text/plain": [
       "<Figure size 1152x720 with 72 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 284424 characters]\n",
      "text/plain": [
       "<Figure size 1152x720 with 72 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 295252 characters]\n",
      "text/plain": [
       "<Figure size 1152x720 with 72 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 19516 characters]\n",
      "text/plain": [
       "<Figure size 1152x720 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda08_12, approximate SST's using lowest factors\n",
    "# CAC Pacific Sea Surface Temperature (SST) dataset\n",
    "\n",
    "# load the data\n",
    "D = np.genfromtxt('../Data/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": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 13276 characters]\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda08_13, create a dataset for cluster analysis\n",
    "\n",
    "# Note: In order to prevent ooverwriting of the originally\n",
    "# distributed dataset, filenames have been changed from\n",
    "# \"../Data/threeclusters.txt\" to \"../Data/mythreeclusters.txt\"\n",
    "# \"../Data/threeclustermeans.txt\" to \"../Data/threeclustermeans.txt\"\n",
    "# change them back if you want ...\n",
    "\n",
    "K=3; # three clusters\n",
    "# 100 samples per cluster\n",
    "N1=100;\n",
    "N2=100;\n",
    "N3=100;\n",
    "N=N1+N2+N3;\n",
    "\n",
    "M=2; # two elements\n",
    "\n",
    "# means set by hand\n",
    "mu = np.zeros((K,M));\n",
    "mu[0,0] = 1.0;\n",
    "mu[0,1] = 2.0;\n",
    "mu[1,0] = 3.0;\n",
    "mu[1,1] = 4.0;\n",
    "mu[2,0] = 4.0;\n",
    "mu[2,1] = 1.0;\n",
    "\n",
    "# sample matrix randomly assigned\n",
    "S = np.zeros((N,M));\n",
    "sc = 0.5;\n",
    "S[0:N1,0:1] = np.random.normal( loc=mu[0,0], scale=sc, size=(N1,1) );\n",
    "S[0:N1,1:2] = np.random.normal( loc=mu[0,1], scale=sc, size=(N1,1) );\n",
    "S[N1:N1+N2,0:1] = np.random.normal( loc=mu[1,0], scale=sc, size=(N2,1) );\n",
    "S[N1:N1+N2,1:2] = np.random.normal( loc=mu[1,1], scale=sc, size=(N2,1) );\n",
    "S[N1+N2:N,0:1] = np.random.normal( loc=mu[2,0], scale=sc, size=(N3,1) );\n",
    "S[N1+N2:N,1:2] = np.random.normal( loc=mu[2,1], scale=sc, size=(N3,1) );\n",
    "\n",
    "# cluster membership\n",
    "k1 = np.zeros((N1,1),dtype=np.int);\n",
    "k2 = np.ones((N2,1),dtype=np.int);\n",
    "k3 = 2*np.ones((N3,1),dtype=np.int);\n",
    "myk = np.concatenate((k1,k2,k3),axis=0);\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(5,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, 5, 0, 5] );\n",
    "for i in range(N):\n",
    "    if myk[i,0]==0:\n",
    "        plt.plot( S[i,0], S[i,1], 'k.', lw=1 );\n",
    "    elif myk[i,0]==1:\n",
    "        plt.plot( S[i,0], S[i,1], 'k+', lw=1 );\n",
    "    else:\n",
    "        plt.plot( S[i,0], S[i,1], 'ko', lw=1 );\n",
    "plt.plot( mu[0:K,0:1], mu[0:K,1:2], 'b^', lw=2 );\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.show();\n",
    "\n",
    "Dm = np.concatenate( (S,myk),axis=1);\n",
    "np.savetxt(\"../Data/mythreeclusters.txt\",Dm, delimiter=\"\\t\");\n",
    "np.savetxt(\"../Data/mythreeclustermeans.txt\",mu, delimiter=\"\\t\");\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cluster means\n",
      "[[1.07859747 2.02539117]\n",
      " [3.04311042 4.05478594]\n",
      " [4.05072009 1.03414835]]\n",
      "mean sample distance\n",
      "0.6470009786672823\n",
      "cluster counts\n",
      "[[100]\n",
      " [100]\n",
      " [100]]\n",
      "mean inter-cluster distance\n",
      "2.0315147021424114\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13940 characters]\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_08_14, cluster analysis w/ prior mean initialization\n",
    "\n",
    "# load the data\n",
    "K = 3;\n",
    "D = np.genfromtxt(\"../Data/threeclusters.txt\", delimiter='\\t');\n",
    "N, M = np.shape(D);\n",
    "M = M-1;\n",
    "S = np.copy( D[0:N,0:2] )\n",
    "myk = np.ndarray.astype(D[0:N,2:3], np.int);\n",
    "mu = np.genfromtxt(\"../Data/threeclustermeans.txt\", delimiter='\\t');\n",
    "\n",
    "# initial means are the prior means\n",
    "initmu = np.zeros((K,M));\n",
    "initmu[0,0] = 1.2;\n",
    "initmu[0,1] = 2.2;\n",
    "initmu[1,0] = 3.4;\n",
    "initmu[1,1] = 4.4;\n",
    "initmu[2,0] = 4.1;\n",
    "initmu[2,1] = 1.1;\n",
    "\n",
    "# k-mean analysis\n",
    "mymu, myk, mydist2, mycount, myKdist2 = eda_kmeans(S, initmu);\n",
    "\n",
    "print('cluster means');\n",
    "print(mymu);\n",
    "print('mean sample distance');\n",
    "print(np.mean(np.sqrt(mydist2)));\n",
    "print('cluster counts');\n",
    "print(mycount);\n",
    "print('mean inter-cluster distance');\n",
    "print(np.mean(np.sqrt(myKdist2)));\n",
    "\n",
    "fig2 = plt.figure(2,figsize=(5,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, 5, 0, 5] );\n",
    "plt.plot( initmu[0:K,0:1], initmu[0:K,1:2], 'r^', lw=2 )\n",
    "for i in range(N):\n",
    "    if myk[i,0]==0:\n",
    "        plt.plot( S[i,0], S[i,1], 'k.', lw=1 );\n",
    "    elif myk[i,0]==1:\n",
    "        plt.plot( S[i,0], S[i,1], 'k+', lw=1 );\n",
    "    else:\n",
    "        plt.plot( S[i,0], S[i,1], 'ko', lw=1 );\n",
    "plt.plot( mymu[0:K,0:1], mymu[0:K,1:2], 'b^', lw=2 )\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cluster means\n",
      "[[3.04311042 4.05478594]\n",
      " [4.05072009 1.03414835]\n",
      " [1.07859747 2.02539117]]\n",
      "mean sample distance\n",
      "0.6470009786672823\n",
      "cluster counts\n",
      "[[100]\n",
      " [100]\n",
      " [100]]\n",
      "mean inter-cluster distance\n",
      "2.0315147021424114\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13708 characters]\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_08_15, cluster analysis w/ random-partition initialization\n",
    "\n",
    "# load the data\n",
    "K = 3;\n",
    "D = np.genfromtxt(\"../Data/threeclusters.txt\", delimiter='\\t');\n",
    "N, M = np.shape(D);\n",
    "M = M-1;\n",
    "S = np.copy( D[0:N,0:2] )\n",
    "myk = np.ndarray.astype(D[0:N,2:3], np.int);\n",
    "mu = np.genfromtxt(\"../Data/threeclustermeans.txt\", delimiter='\\t');\n",
    "\n",
    "# initial means are from random partition of the data set\n",
    "initmu = eda_random_partition_mean(S,K);\n",
    "\n",
    "# k-mean analysis\n",
    "mymu, myk, mydist, mycount, myKdist = eda_kmeans(S, initmu);\n",
    "\n",
    "print('cluster means');\n",
    "print(mymu);\n",
    "print('mean sample distance');\n",
    "print(np.mean(np.sqrt(mydist)));\n",
    "print('cluster counts');\n",
    "print(mycount);\n",
    "print('mean inter-cluster distance');\n",
    "print(np.mean(np.sqrt(myKdist)));\n",
    "\n",
    "fig3 = plt.figure(3,figsize=(5,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, 5, 0, 5] );\n",
    "plt.plot( initmu[0:K,0:1], initmu[0:K,1:2], 'r^', lw=4 )\n",
    "for i in range(N):\n",
    "    if myk[i,0]==0:\n",
    "        plt.plot( S[i,0], S[i,1], 'k.', lw=1 );\n",
    "    elif myk[i,0]==1:\n",
    "        plt.plot( S[i,0], S[i,1], 'k+', lw=1 );\n",
    "    else:\n",
    "        plt.plot( S[i,0], S[i,1], 'ko', lw=1 );\n",
    "plt.plot( mymu[0:K,0:1], mymu[0:K,1:2], 'b^', lw=4 );\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cluster means\n",
      "[[3.04311042 4.05478594]\n",
      " [1.07859747 2.02539117]\n",
      " [4.05072009 1.03414835]]\n",
      "mean sample distance\n",
      "0.6470009786672823\n",
      "cluster counts\n",
      "[[100]\n",
      " [100]\n",
      " [100]]\n",
      "mean inter-cluster distance\n",
      "2.0315147021424114\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13528 characters]\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_08_16, cluster analysis w/ forgy initialization\n",
    "\n",
    "# load the data\n",
    "K = 3;\n",
    "D = np.genfromtxt(\"../Data/threeclusters.txt\", delimiter='\\t');\n",
    "N, M = np.shape(D);\n",
    "M = M-1;\n",
    "S = np.copy( D[0:N,0:2] )\n",
    "myk = np.ndarray.astype(D[0:N,2:3], np.int);\n",
    "mu = np.genfromtxt(\"../Data/threeclustermeans.txt\", delimiter='\\t');\n",
    "    \n",
    "# initial means are randomly-chosen samples\n",
    "initmu = eda_forgy_mean(S,K);\n",
    "\n",
    "# k-mean analysis\n",
    "mymu, myk, mydist2, mycount, myKdist2 = eda_kmeans(S, initmu);\n",
    "\n",
    "print('cluster means');\n",
    "print(mymu);\n",
    "print('mean sample distance');\n",
    "print(np.mean(np.sqrt(mydist2)));\n",
    "print('cluster counts');\n",
    "print(mycount);\n",
    "print('mean inter-cluster distance');\n",
    "print(np.mean(np.sqrt(myKdist2)));\n",
    "\n",
    "fig4 = plt.figure(4,figsize=(5,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, 5, 0, 5] );\n",
    "plt.plot( initmu[0:K,0:1], initmu[0:K,1:2], 'r^', lw=4 )\n",
    "for i in range(N):\n",
    "    if myk[i,0]==0:\n",
    "        plt.plot( S[i,0], S[i,1], 'k.', lw=1 );\n",
    "    elif myk[i,0]==1:\n",
    "        plt.plot( S[i,0], S[i,1], 'k+', lw=1 );\n",
    "    else:\n",
    "        plt.plot( S[i,0], S[i,1], 'ko', lw=1 );\n",
    "plt.plot( mymu[0:K,0:1], mymu[0:K,1:2], 'b^', lw=4 )\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cluster means\n",
      "[[3.04311042 4.05478594]\n",
      " [4.05072009 1.03414835]\n",
      " [1.07859747 2.02539117]]\n",
      "mean sample distance\n",
      "0.6470009786672823\n",
      "cluster counts\n",
      "[[100]\n",
      " [100]\n",
      " [100]]\n",
      "mean inter-cluster distance\n",
      "2.0315147021424114\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13700 characters]\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_08_17, cluster analysis w/ Bradley-Fayyad initialization\n",
    "\n",
    "# load the data\n",
    "K = 3;\n",
    "D = np.genfromtxt(\"../Data/threeclusters.txt\", delimiter='\\t');\n",
    "N, M = np.shape(D);\n",
    "M = M-1;\n",
    "S = np.copy( D[0:N,0:2] )\n",
    "myk = np.ndarray.astype(D[0:N,2:3], np.int);\n",
    "mu = np.genfromtxt(\"../Data/threeclustermeans.txt\", delimiter='\\t');\n",
    "\n",
    "# bradley-fayyad parameters\n",
    "Nsets = 100;     # number of small sets\n",
    "fraction = 0.05; # size of small sets = fraction of N\n",
    "\n",
    "# initial means are bradley-fayyad means\n",
    "initmu = eda_bradley_fayyad_mean(S,K,Nsets,fraction);\n",
    "\n",
    "# k-mean analysis\n",
    "mymu, myk, mydist2, mycount, myKdist2 = eda_kmeans(S, initmu);\n",
    "\n",
    "print('cluster means');\n",
    "print(mymu);\n",
    "print('mean sample distance');\n",
    "print(np.mean(np.sqrt(mydist2)));\n",
    "print('cluster counts');\n",
    "print(mycount);\n",
    "print('mean inter-cluster distance');\n",
    "print(np.mean(np.sqrt(myKdist2)));\n",
    "\n",
    "fig5 = plt.figure(5,figsize=(5,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, 5, 0, 5] );\n",
    "plt.plot( initmu[0:K,0:1], initmu[0:K,1:2], 'r^', lw=4 )\n",
    "for i in range(N):\n",
    "    if myk[i,0]==0:\n",
    "        plt.plot( S[i,0], S[i,1], 'k.', lw=1 );\n",
    "    elif myk[i,0]==1:\n",
    "        plt.plot( S[i,0], S[i,1], 'k+', lw=1 );\n",
    "    else:\n",
    "        plt.plot( S[i,0], S[i,1], 'ko', lw=1 );\n",
    "plt.plot( mymu[0:K,0:1], mymu[0:K,1:2], 'b^', lw=4 )\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cluster means\n",
      "[[4.05072009 1.03414835]\n",
      " [2.59186649 4.05694587]\n",
      " [1.15537449 1.29173102]\n",
      " [1.36411249 2.41193976]\n",
      " [3.51277248 4.05253784]\n",
      " [0.60469989 2.03226719]]\n",
      "mean sample distance\n",
      "0.53212144710493\n",
      "cluster counts\n",
      "[[100]\n",
      " [ 51]\n",
      " [ 24]\n",
      " [ 45]\n",
      " [ 49]\n",
      " [ 31]]\n",
      "mean inter-cluster distance\n",
      "2.0915969746599563\n",
      "inter-cluster distances\n",
      "[[0.         3.35642052 2.90678087 3.01929952 3.06595217 3.58765891]\n",
      " [3.35642052 0.         3.11607484 2.0526629  0.92091654 2.83692701]\n",
      " [2.90678087 3.11607484 0.         1.13949075 3.63034152 0.92284145]\n",
      " [3.01929952 2.0526629  1.13949075 0.         2.70338709 0.84903402]\n",
      " [3.06595217 0.92091654 3.63034152 2.70338709 0.         3.54095745]\n",
      " [3.58765891 2.83692701 0.92284145 0.84903402 3.54095745 0.        ]]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 15684 characters]\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_08_18, more clusters than actually occur\n",
    "\n",
    "# load the data\n",
    "K = 3;\n",
    "D = np.genfromtxt(\"../Data/threeclusters.txt\", delimiter='\\t');\n",
    "N, M = np.shape(D);\n",
    "M = M-1;\n",
    "S = np.copy( D[0:N,0:2] )\n",
    "myk = np.ndarray.astype(D[0:N,2:3], np.int);\n",
    "mu = np.genfromtxt(\"../Data/threeclustermeans.txt\", delimiter='\\t');\n",
    "\n",
    "\n",
    "# six clusers\n",
    "K6=6;\n",
    "\n",
    "# initial means are forgy means\n",
    "initmu = eda_forgy_mean(S,K6);\n",
    "\n",
    "# k-mean analysis\n",
    "mymu, myk, mydist2, mycount, myKdist2 = eda_kmeans(S, initmu);\n",
    "\n",
    "print('cluster means');\n",
    "print(mymu);\n",
    "print('mean sample distance');\n",
    "print(np.mean(np.sqrt(mydist2)));\n",
    "print('cluster counts');\n",
    "print(mycount);\n",
    "print('mean inter-cluster distance');\n",
    "print(np.mean(np.sqrt(myKdist2)));\n",
    "print('inter-cluster distances');\n",
    "print(np.sqrt(myKdist2));\n",
    "\n",
    "fig4 = plt.figure(4,figsize=(5,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, 5, 0, 5] );\n",
    "for i in range(N):\n",
    "    if myk[i,0]==0:\n",
    "        plt.plot( S[i,0], S[i,1], 'k.', lw=1 );\n",
    "    elif myk[i,0]==1:\n",
    "        plt.plot( S[i,0], S[i,1], 'k+', lw=1 );\n",
    "    elif myk[i,0]==2:\n",
    "        plt.plot( S[i,0], S[i,1], 'k^', lw=1 );\n",
    "    elif myk[i,0]==3:\n",
    "        plt.plot( S[i,0], S[i,1], 'kx', lw=1 );\n",
    "    elif myk[i,0]==4:\n",
    "        plt.plot( S[i,0], S[i,1], 'kv', lw=1 );\n",
    "    else:\n",
    "        plt.plot( S[i,0], S[i,1], 'ks', lw=1 );\n",
    "plt.plot( mymu[0:K6,0:1], mymu[0:K6,1:2], 'r^', lw=2 )\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K mean_dist\n",
      "1 2.450872469766111\n",
      "2 2.159369027284355\n",
      "3 1.87730362348076\n",
      "4 1.6844961480997582\n",
      "5 1.5978670360338794\n",
      "6 1.5629487428459847\n",
      "7 1.4514984830833648\n",
      "8 1.4064497014366697\n",
      "9 1.3781130639717087\n",
      "10 1.3463808953032275\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 12700 characters]\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_08_19, clusters in the Atlantic Rocks dataset\n",
    "\n",
    "# read data\n",
    "S = np.genfromtxt('../Data/rocks.txt', delimiter='\\t');\n",
    "N, M = np.shape(S);\n",
    "\n",
    "# examine mean-sample distances for sequence of number\n",
    "# of clusters\n",
    "KMAX = 11;\n",
    "Klist = np.zeros((KMAX,1));\n",
    "Dlist = np.zeros((KMAX,1));\n",
    "print('K','mean_dist');\n",
    "for K in range(1,KMAX):\n",
    "    initmu = eda_forgy_mean(S,K);\n",
    "    mymu, myk, mydist2, mycount, myKdist2 = eda_kmeans(S, initmu);\n",
    "    Klist[K,0]=K;\n",
    "    Dlist[K,0]=np.mean(np.sqrt(mydist2));\n",
    "    print(K,Dlist[K,0]);\n",
    "\n",
    "# plot of distances\n",
    "fig1 = plt.figure(1,figsize=(10,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, KMAX, 0, 3] );\n",
    "plt.plot( Klist[1:KMAX,0:1], Dlist[1:KMAX,0:1], 'k-', lw=2 )\n",
    "plt.plot( Klist[1:KMAX,0:1], Dlist[1:KMAX,0:1], 'ko', lw=2 )\n",
    "plt.xlabel('K');\n",
    "plt.ylabel('mean Rij');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cluster\tSiO2\tTiO2\tAl203\tFeO-t\tMgO\tCaO\tNa2O\tK2O\tcount\n",
      "0\t50.50\t1.17\t15.79\t8.90\t8.12\t12.17\t2.37\t0.20\t1540\n",
      "1\t50.31\t1.60\t15.47\t9.89\t7.60\t11.16\t2.85\t0.18\t1590\n",
      "2\t48.11\t1.04\t17.58\t8.10\t6.76\t13.38\t2.23\t0.19\t360\n",
      "3\t51.03\t1.27\t14.77\t10.21\t7.54\t12.06\t2.29\t0.15\t1260\n",
      "4\t52.00\t1.79\t15.23\t9.83\t6.36\t10.34\t3.03\t0.45\t419\n",
      "5\t48.29\t0.84\t16.11\t9.05\t10.04\t11.93\t2.17\t0.11\t314\n",
      "6\t50.69\t1.62\t14.08\t12.27\t6.77\t11.40\t2.32\t0.16\t819\n",
      "7\t44.90\t0.06\t2.51\t8.20\t40.22\t1.52\t0.13\t0.03\t54\n"
     ]
    }
   ],
   "source": [
    "# edapy_08_20, print 8 clusters in the Atlantic Rocks dataset\n",
    "\n",
    "# read the data\n",
    "S = np.genfromtxt('../Data/rocks.txt', delimiter='\\t')\n",
    "N, M = np.shape(S);\n",
    "\n",
    "K=8; # 8 clusters\n",
    "initmu = eda_forgy_mean(S,K); # forgy initial means\n",
    "\n",
    "# k-mean analysis\n",
    "mymu, myk, mydist, mycount, myKdist = eda_kmeans(S, initmu);\n",
    "\n",
    "# print table of cluster means\n",
    "names = [\"SiO2\", \"TiO2\", \"Al203\", \"FeO-t\", \"MgO\", \"CaO\", \"Na2O\", \"K2O\" ];\n",
    "print('cluster\\t', end='');\n",
    "for j in range(M):\n",
    "    print(\"%s\\t\" % (names[j]), end='');\n",
    "print('count');\n",
    "for i in range(K):\n",
    "    print(\"%d\\t\" % (i), end='');\n",
    "    for j in range(M):\n",
    "        print(\"%.2f\\t\" % (mymu[i,j]), end='');\n",
    "    print(\"%d\" % (mycount[i,0]));\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 47036 characters]\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_08_21, 3D plot of clusters in Atlantic Rocks dataset\n",
    "\n",
    "# read the data\n",
    "S = np.genfromtxt('../Data/rocks.txt', delimiter='\\t')\n",
    "N, M = np.shape(S);\n",
    "names = [\"SiO2\", \"TiO2\", \"Al203\", \"FeO-t\", \"MgO\", \"CaO\", \"Na2O\", \"K2O\" ];\n",
    "\n",
    "K=8; # 8 clusters\n",
    "initmu = eda_forgy_mean(S,K); # forgy initial means\n",
    "\n",
    "# k-mean analysis\n",
    "mymu, myk, mydist, mycount, myKdist = eda_kmeans(S, initmu);\n",
    "\n",
    "# elements to plot (and their scale)\n",
    "\n",
    "m1 = 2; # Al203\n",
    "xmin = 0.0\n",
    "xmax = 50.0;\n",
    "\n",
    "m2 = 4; # MgO\n",
    "ymin = 0.0\n",
    "ymax = 50.0;\n",
    "\n",
    "m3 = 5; # Ca0\n",
    "zmin = 0.0\n",
    "zmax = 20.0;\n",
    "\n",
    "# 3D plot\n",
    "\n",
    "fig2 = plt.figure(2,figsize=(10,5));\n",
    "ax1 = fig2.add_subplot(121, projection='3d')\n",
    "ax1.axes.set_xlim3d(left=xmin, right=xmax);\n",
    "ax1.axes.set_ylim3d(bottom=ymin, top=ymax);\n",
    "ax1.axes.set_zlim3d(bottom=zmin, top=zmax);\n",
    "ax1.view_init(1, 30);\n",
    "\n",
    "ax1.scatter( xs=S[:,m1], ys=S[:,m2], zs=S[:,m3], color=(0,0,0) );\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymin,ymin], zs=[zmin,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymin,ymax], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmin,xmax], ys=[ymin,ymin], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymax,ymax], zs=[zmax,zmin], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymax,ymin], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmin], ys=[ymax,ymax], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmax,xmin], ys=[ymin,ymin], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymin,ymin], zs=[zmax,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymax,ymin], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymax,ymax], zs=[zmax,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymax,ymin], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmin], ys=[ymax,ymax], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.set_xlabel(names[m1]);\n",
    "ax1.set_ylabel(names[m2]);\n",
    "ax1.set_zlabel(names[m3]);\n",
    "\n",
    "ax1 = fig2.add_subplot(122, projection='3d')\n",
    "ax1.axes.set_xlim3d(left=xmin, right=xmax);\n",
    "ax1.axes.set_ylim3d(bottom=ymin, top=ymax);\n",
    "ax1.axes.set_zlim3d(bottom=zmin, top=zmax);\n",
    "ax1.view_init(1, 30);\n",
    "\n",
    "ax1.scatter( xs=mymu[:,m1], ys=mymu[:,m2], zs=mymu[:,m3], color=(0,0,0) );\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymin,ymin], zs=[zmin,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymin,ymax], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmin,xmax], ys=[ymin,ymin], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymax,ymax], zs=[zmax,zmin], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymax,ymin], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmin], ys=[ymax,ymax], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmax,xmin], ys=[ymin,ymin], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymin,ymin], zs=[zmax,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymax,ymin], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymax,ymax], zs=[zmax,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymax,ymin], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmin], ys=[ymax,ymax], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.set_xlabel(names[m1]);\n",
    "ax1.set_ylabel(names[m2]);\n",
    "ax1.set_zlabel(names[m3]);\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "true patterns\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 40712 characters]\n",
      "text/plain": [
       "<Figure size 720x360 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "amplitudes of true patterns\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 32884 characters]\n",
      "text/plain": [
       "<Figure size 720x360 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 120496 characters]\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 10792 characters]\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 85972 characters]\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EOFs\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 49676 characters]\n",
      "text/plain": [
       "<Figure size 720x360 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "amplitudes of EOFs\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 60560 characters]\n",
      "text/plain": [
       "<Figure size 720x360 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note that the factors do not match the original patterns,\n",
      "although they do bear some resemblence to them.\n",
      "Thats because any linear combination of factors will do.\n",
      "Note also that fv4 and cu4 are very noisy. Thats because the\n",
      "p=4 singular value is close to the noise level.\n"
     ]
    }
   ],
   "source": [
    "# edama_08_22, Script to support Cribsheet 8.2 on EOF's\n",
    "\n",
    "# standard set-up of v-axis\n",
    "M=100;\n",
    "Dv = 1.0;\n",
    "vmin = 0.0;\n",
    "vmax = vmin + Dv*(M-1);\n",
    "v = eda_cvec(np.linspace(vmin,vmax,M));\n",
    "\n",
    "# standard set-up of u-axis\n",
    "N=1000;\n",
    "Du = 1.0;\n",
    "umin = 0.0;\n",
    "umax = umin + Du*(N-1);\n",
    "u = eda_cvec(np.linspace(umin,umax,N));\n",
    "\n",
    "# gre-scale colormap\n",
    "bw = np.zeros((256,4));\n",
    "vv = (255 - np.linspace( 0, 255, 256 ))/255;\n",
    "bw[:,0] = vv;\n",
    "bw[:,1] = vv;\n",
    "bw[:,2] = vv;\n",
    "bw[:,3] = np.ones(256);\n",
    "bwcmap = ListedColormap(bw);\n",
    "\n",
    "# Part 1: make synthetic data\n",
    "\n",
    "# four true patterns\n",
    "fv1 = np.mod(v,10)/10;\n",
    "v0 = (vmin+vmax)/2;\n",
    "s2 = (vmax-vmin)/5;\n",
    "fv2 = np.exp(-0.5*np.power((v-v0)/s2,2));\n",
    "v1 = (vmin+vmax)/4;\n",
    "v2 = 3*(vmin+vmax)/4;\n",
    "s2 = (vmax-vmin)/10;\n",
    "fv3 = np.exp(-0.5*np.power((v-v1)/s2,2))-np.exp(-0.5*np.power((v-v2),2)/s2);\n",
    "L = vmax-vmin;\n",
    "n=10;\n",
    "fv4 = np.cos(n*pi*v/L);\n",
    "\n",
    "# plot u dependence of true patterns\n",
    "print('true patterns');\n",
    "fig1 = plt.figure(1,figsize=(10,5));\n",
    "fvlist = [fv1, fv2, fv3, fv4];\n",
    "for i in range(4):\n",
    "    fvi = fvlist[i];\n",
    "    ax1 = plt.subplot(4,1,i+1);\n",
    "    fmax = np.max(np.abs(fvi));\n",
    "    plt.axis( [vmin, vmax, -fmax, fmax] );\n",
    "    plt.plot( v, fvi, 'k-', lw=2 );\n",
    "    plt.xlabel('v');\n",
    "    plt.ylabel(\"fv%d\" % (i));\n",
    "plt.show()\n",
    "\n",
    "# true amplitudes of patterns\n",
    "L = umax-umin;\n",
    "cu1 = (u-umin)/L;\n",
    "cu2 = np.power((u-umin),0.25)/(L**0.25);\n",
    "n=10;\n",
    "cu3 = np.cos(n*pi*u/L);\n",
    "cu4 = (100-np.mod(u,100))/100;\n",
    "\n",
    "# plot true v dependence of amplitudes\n",
    "print('amplitudes of true patterns');\n",
    "fig2 = plt.figure(2,figsize=(10,5));\n",
    "culist = [cu1, cu2, cu3, cu4];\n",
    "for i in range(4):\n",
    "    cui = culist[i];\n",
    "    ax1 = plt.subplot(4,1,i+1);\n",
    "    cmax = np.max(np.abs(cui));\n",
    "    plt.axis( [umin, umax, -cmax, cmax] );\n",
    "    plt.plot( u, cui, 'k-', lw=2 );\n",
    "    plt.xlabel('u');\n",
    "    plt.ylabel(\"cu%d\" % (i));\n",
    "plt.show()\n",
    "\n",
    "# build hypothetical sample matric\n",
    "S = np.matmul(\n",
    "        np.concatenate( (cu1, cu2, cu3, cu4), axis=1),\n",
    "        np.concatenate( (fv1.T, fv2.T, fv3.T, fv4.T), axis=0) );\n",
    "Smax = np.max(np.abs(S));\n",
    "S = S + np.random.normal(loc=0,scale=(0.1*Smax),size=(N,M));\n",
    "\n",
    "# plot S\n",
    "fig3 = plt.figure(3,figsize=(10,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "Smax = np.max(np.abs(S));\n",
    "plt.axis( [vmin,vmax,umin,umax] );\n",
    "plt.imshow(S,cmap=bwcmap,vmin=-Smax,vmax=Smax,extent=(vmin,vmax,umin,umax),aspect='auto');\n",
    "plt.xlabel('v');\n",
    "plt.ylabel('u');\n",
    "plt.title('S');\n",
    "plt.show();\n",
    "\n",
    "# Part 2: singular value decompsition\n",
    "# F is a matrix of the EOF's and\n",
    "# C is a matrix of their amplitudes\n",
    "\n",
    "# factors F and loadings F\n",
    "[U, sigma, VT] = la.svd(S,full_matrices=False);\n",
    "F = VT; # EOS's (factors)\n",
    "C = np.matmul( U, np.diag(sigma) ); # amplitudes (loadings)\n",
    "\n",
    "fig4 = plt.figure(4,figsize=(10,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, M, 0, np.max(sigma)] );\n",
    "plt.plot( np.linspace(1,M,M), sigma, 'k-', lw=1 );\n",
    "plt.plot( np.linspace(1,M,M), sigma, 'ko', lw=2 );\n",
    "plt.xlabel('i');\n",
    "plt.ylabel('sigma i');\n",
    "plt.title('singular values');\n",
    "\n",
    "# Part 3: selection of number P of significant EOS's\n",
    "\n",
    "# On the basis of the last plot, P=4\n",
    "P = 4;\n",
    "FP = np.copy( F[0:P,:] );\n",
    "CP = np.copy( C[:,0:P] );\n",
    "SP = np.matmul( CP, FP );\n",
    "\n",
    "# plot SP\n",
    "fig5 = plt.figure(5,figsize=(10,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "SPmax = np.max(np.abs(SP));\n",
    "plt.axis( [vmin,vmax,umin,umax] );\n",
    "plt.imshow(SP,cmap=bwcmap,vmin=-Smax,vmax=Smax,extent=(vmin,vmax,umin,umax),aspect='auto');\n",
    "plt.xlabel('v');\n",
    "plt.ylabel('u');\n",
    "plt.title('SP');\n",
    "plt.show();\n",
    "\n",
    "# Part 3: plot EOF's as a function of v\n",
    "# and their amplitudes as a function of u\n",
    "\n",
    "# plot EOF's\n",
    "print('EOFs');\n",
    "fig6 = plt.figure(6,figsize=(10,5));\n",
    "for i in range(P):\n",
    "    ax1 = plt.subplot(4,1,i+1);\n",
    "    fvi = np.copy(FP[i:i+1,0:M]).T;\n",
    "    plt.axis( [vmin, vmax, -max(abs(fvi)), max(abs(fvi))])\n",
    "    plt.plot( v, fvi, 'k-', lw=2);\n",
    "    plt.xlabel('v');\n",
    "    plt.ylabel(\"fv%d\" % (i));\n",
    "plt.show();\n",
    "\n",
    "# plot amplitudes\n",
    "print('amplitudes of EOFs');\n",
    "fig7 = plt.figure(7,figsize=(10,5));\n",
    "for i in range(P):\n",
    "    ax1 = plt.subplot(4,1,i+1);\n",
    "    cui=np.copy(CP[0:N,i:i+1]);\n",
    "    plt.axis( [umin, umax, -max(abs(cui)), max(abs(cui))]);\n",
    "    plt.plot( u, cui, 'k-', lw=2);\n",
    "    plt.xlabel('u');\n",
    "    plt.ylabel(\"cu%d\" % (i));\n",
    "plt.show();\n",
    "\n",
    "print('Note that the factors do not match the original patterns,');\n",
    "print('although they do bear some resemblence to them.');\n",
    "print('That''s because any linear combination of factors will do.');\n",
    "print('Note also that fv4 and cu4 are very noisy. That''s because the');\n",
    "print('p=4 singular value is close to the noise level.');"
   ]
  },
  {
   "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
}
