svm.py

   1 #
   2 # Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
   3 # Angel Farguell (angel.farguell@gmail.com)
   4 #
   5 # to install:
   6 #       conda install scikit-learn
   7 #       conda install scikit-image
   8
   9 from sklearn import svm
  10 from scipy import interpolate, spatial
  11 import matplotlib.pyplot as plt
  12 import matplotlib.font_manager
  13 import matplotlib.colors as colors
  14 from mpl_toolkits.mplot3d import axes3d
  15 from mpl_toolkits.mplot3d.art3d import Poly3DCollection
  16 import numpy as np
  17 from time import time
  18 from infrared_perimeters import process_infrared_perimeters
  19 import sys
  20 import saveload as sl
  21
  22 def preprocess_data_svm(lons, lats, U, L, T, scale, time_num_granules, C=None):
  23     """
  24     Preprocess satellite data from JPSSD and setup to use in Support Vector Machine
  25
  26     :param lons: longitud grid
  27     :param lats: latitde grid
  28     :param U: upper bound grid
  29     :param L: lower bound grid
  30     :param T: mask grid
  31     :param scale: time scales
  32     :param time_num_granules: times of the granules
  33     :return X: matrix of features for SVM
  34     :return y: vector of labels for SVM
  35
  36     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
  37     Angel Farguell (angel.farguell@gmail.com), 2019-04-01
  38     """
  39
  40     # Flatten coordinates
  41     lon = np.reshape(lons,np.prod(lons.shape)).astype(float)
  42     lat = np.reshape(lats,np.prod(lats.shape)).astype(float)
  43
  44     # Temporal scale to days
  45     tscale = 24*3600
  46     U = U/tscale
  47     L = L/tscale
  48
  49     # Ensuring U>=L always
  50     print 'U>L: ',(U>L).sum()
  51     print 'U<L: ',(U<L).sum()
  52     print 'U==L: ',(U==L).sum()
  53
  54     # Reshape to 1D
  55     uu = np.reshape(U,np.prod(U.shape))
  56     ll = np.reshape(L,np.prod(L.shape))
  57     tt = np.reshape(T,np.prod(T.shape))
  58
  59     # Maximum and minimums to NaN data
  60     uu[uu==uu.max()] = np.nan
  61     ll[ll==ll.min()] = np.nan
  62
  63     # Mask created during computation of lower and upper bounds
  64     mk = tt==scale[1]-scale[0]
  65     # Masking upper bounds outside the mask
  66     uu[mk] = np.nan
  67     # Creating maximum value considered of the upper bounds
  68     nuu = uu[~np.isnan(uu)]
  69     muu = nuu.max() # could be a different value like a mean value
  70     # Create a mask with lower bound less than the previous maximum upper bound value
  71     with np.errstate(invalid='ignore'):
  72         low = (ll <= muu)
  73     if low.sum() > 10000:
  74         # Create a mask with all False of low size
  75         mask = np.repeat(False,len(low[low == True]))
  76         # Take just a subset of the nodes
  77         clear_level = 50
  78         mask[0::clear_level] = True
  79         # Mask the subset
  80         low[low == True] = mask
  81     # Eliminate all the previous elements from the mask
  82     mk[low] = False
  83     # Masking lower bounds outside the mask
  84     ll[mk] = np.nan
  85
  86     # Values different than NaN in the upper and lower bounds
  87     um = np.array(~np.isnan(uu))
  88     lm = np.array(~np.isnan(ll))
  89     # Define all the x, y, and z components of upper and lower bounds
  90     ux = lon[um]
  91     uy = lat[um]
  92     uz = uu[um]
  93     lx = lon[lm]
  94     ly = lat[lm]
  95     lz = ll[lm]
  96
  97     # Create the data to call SVM3 function from svm3test.py
  98     X = np.c_[np.concatenate((lx,ux)),np.concatenate((ly,uy)),np.concatenate((lz,uz))]
  99     y = np.concatenate((-np.ones(len(lx)),np.ones(len(ux))))
 100     # Print the shape of the data
 101     print 'shape X: ', X.shape
 102     print 'shape y: ', y.shape
 103
 104     if C is None:
 105         c = 80*np.ones(y.shape)
 106     else:
 107         c = np.concatenate((C[0].ravel()[lm],C[1].ravel()[um]))
 108
 109     # Clean data if not in bounding box
 110     bbox = (lon.min(),lon.max(),lat.min(),lat.max(),time_num_granules)
 111     geo_mask = np.logical_and(np.logical_and(np.logical_and(X[:,0] >= bbox[0],X[:,0] <= bbox[1]), X[:,1] >= bbox[2]), X[:,1] <= bbox[3])
 112     btime = (0,(scale[1]-scale[0])/tscale)
 113     time_mask = np.logical_and(X[:,2] >= btime[0], X[:,2] <= btime[1])
 114     whole_mask = np.logical_and(geo_mask, time_mask)
 115     X = X[whole_mask,:]
 116     y = y[whole_mask]
 117     c = c[whole_mask]
 118
 119     return X,y,c
 120
 121 def make_fire_mesh(fxlon, fxlat, it, nt):
 122     """
 123     Create a mesh of points to evaluate the decision function
 124
 125     :param fxlon: data to base x-axis meshgrid on
 126     :param fxlat: data to base y-axis meshgrid on
 127     :param it: data to base z-axis meshgrid on
 128     :param nt: tuple of number of nodes at each direction, optional
 129     :param coarse: coarsening of the fire mesh
 130     :return xx, yy, zz: ndarray
 131
 132     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 133     Angel Farguell (angel.farguell@gmail.com), 2019-04-01
 134     """
 135
 136     xx = np.repeat(fxlon[:, :, np.newaxis], nt, axis=2)
 137     yy = np.repeat(fxlat[:, :, np.newaxis], nt, axis=2)
 138     tt = np.linspace(it[0],it[1],nt)
 139     zz = np.swapaxes(np.swapaxes(np.array([np.ones(fxlon.shape)*t for t in tt]),0,1),1,2)
 140
 141     return xx, yy, zz
 142
 143 def make_meshgrid(x, y, z, s=(50,50,50), exp=.1):
 144     """
 145     Create a mesh of points to evaluate the decision function
 146
 147     :param x: data to base x-axis meshgrid on
 148     :param y: data to base y-axis meshgrid on
 149     :param z: data to base z-axis meshgrid on
 150     :param s: tuple of number of nodes at each direction, optional
 151     :param exp: extra percentage of time steps in each direction (between 0 and 1), optional
 152     :return xx, yy, zz: ndarray
 153
 154     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 155     Angel Farguell (angel.farguell@gmail.com), 2019-02-20
 156     Modified version of:
 157     https://scikit-learn.org/stable/auto_examples/svm/plot_iris.html#sphx-glr-auto-examples-svm-plot-iris-py
 158     """
 159
 160     if not isinstance(s, tuple):
 161         print '>> FAILED <<'
 162         print 'The number of nodes at each direction is not a tuple: ', s
 163         sys.exit(1)
 164     # number of nodes in each direction
 165     sx, sy, sz = np.array(s).astype(int)
 166     # extra step sizes in each direction
 167     brx = int(sx * exp)
 168     bry = int(sy * exp)
 169     brz = int(sz * exp)
 170     # grid lengths in each directon
 171     lx = x.max() - x.min()
 172     ly = y.max() - y.min()
 173     lz = z.max() - z.min()
 174     # grid resolutions in each direction
 175     hx = lx / (sx - 2*brx - 1)
 176     hy = ly / (sy - 2*bry - 1)
 177     hz = lz / (sz - 2*brz - 1)
 178     # extrem values for each dimension
 179     x_min, x_max = x.min() - brx * hx, x.max() + brx * hx
 180     y_min, y_max = y.min() - bry * hy, y.max() + bry * hy
 181     z_min, z_max = z.min() - brz * hz, z.max() + brz * hz
 182     # generating the mesh grid
 183     xx, yy, zz = np.meshgrid(np.linspace(y_min, y_max, sy),
 184                              np.linspace(x_min, x_max, sx),
 185                              np.linspace(z_min, z_max, sz))
 186     return xx, yy, zz
 187
 188 def frontier(clf, xx, yy, zz, bal=.5, plot_decision = False, plot_poly=False, using_weights=False):
 189     """
 190     Compute the surface decision frontier for a classifier.
 191
 192     :param clf: a classifier
 193     :param xx: meshgrid ndarray
 194     :param yy: meshgrid ndarray
 195     :param zz: meshgrid ndarray
 196     :param bal: number between 0 and 1, balance between lower and upper bounds in decision function (in case not level 0)
 197     :param plot_decision: boolean of plotting decision volume
 198     :param plot_poly: boolean of plotting polynomial approximation
 199     :return F: 2D meshes with xx, yy coordinates and the hyperplane z which gives decision functon 0
 200
 201     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 202     Angel Farguell (angel.farguell@gmail.com), 2019-02-20
 203     Modified version of:
 204     https://www.semipol.de/2015/10/29/SVM-separating-hyperplane-3d-matplotlib.html
 205     """
 206
 207     # Creating the 3D grid
 208     XX = np.c_[xx.ravel(), yy.ravel(), zz.ravel()]
 209
 210     # Evaluating the decision function
 211     print '>> Evaluating the decision function...'
 212     sys.stdout.flush()
 213     t_1 = time()
 214     if using_weights:
 215         from libsvm_weights.python.svmutil import svm_predict
 216         _, _, p_vals = svm_predict([], XX, clf)
 217         ZZ = np.array([p[0] for p in p_vals])
 218     else:
 219         ZZ = clf.decision_function(XX)
 220     t_2 = time()
 221     print 'elapsed time: %ss.' % str(abs(t_2-t_1))
 222     hist = np.histogram(ZZ)
 223     print 'counts: ', hist[0]
 224     print 'values: ', hist[1]
 225     print 'decision function range: ', ZZ.min(), '~', ZZ.max()
 226
 227     # Reshaping decision function volume
 228     Z = ZZ.reshape(xx.shape)
 229     print 'decision function shape: ', Z.shape
 230
 231     if plot_decision:
 232         try:
 233             from skimage import measure
 234             from shiftcmap import shiftedColorMap
 235             verts, faces, normals, values = measure.marching_cubes_lewiner(Z, level=0, allow_degenerate=False)
 236             # Scale and transform to actual size of the interesting volume
 237             h = np.divide([xx.max()-xx.min(), yy.max() - yy.min(), zz.max() - zz.min()],np.array(xx.shape)-1)
 238             verts = verts * h
 239             verts = verts + [xx.min(), yy.min(), zz.min()]
 240             mesh = Poly3DCollection(verts[faces], facecolor='orange', alpha=.9)
 241             fig = plt.figure()
 242             ax = fig.gca(projection='3d')
 243             fig.suptitle("Decision volume")
 244             col = [(0, 0, 1), (.5, .5, .5), (1, 0, 0)]
 245             cm = colors.LinearSegmentedColormap.from_list('BuRd',col,N=100)
 246             midpoint = 1 - ZZ.max() / (ZZ.max() + abs(ZZ.min()))
 247             shiftedcmap = shiftedColorMap(cm, midpoint=midpoint, name='shifted')
 248             X = xx.ravel()
 249             Y = yy.ravel()
 250             T = zz.ravel()
 251             kk = 50
 252             p = ax.scatter(X[0::kk], Y[0::kk], T[0::kk], c=ZZ[0::kk], s=.1, alpha=.4, cmap=shiftedcmap)
 253             cbar = fig.colorbar(p)
 254             cbar.set_label('decision function value', rotation=270, labelpad=20)
 255             ax.add_collection3d(mesh)
 256             ax.set_zlim([xx.min(),xx.max()])
 257             ax.set_ylim([yy.min(),yy.max()])
 258             ax.set_zlim([zz.min(),zz.max()])
 259             ax.set_xlabel("Longitude normalized")
 260             ax.set_ylabel("Latitude normalized")
 261             ax.set_zlabel("Time normalized")
 262             plt.savefig('decision.png')
 263         except Exception as e:
 264             print 'Warning: something went wrong when plotting...'
 265             print e
 266
 267     if plot_poly:
 268         fig = plt.figure()
 269     # Computing fire arrival time from previous decision function
 270     print '>> Computing fire arrival time...'
 271     sys.stdout.flush()
 272     t_1 = time()
 273     # xx 2-dimensional array
 274     Fx = xx[:, :, 0]
 275     # yy 2-dimensional array
 276     Fy = yy[:, :, 0]
 277     # zz 1-dimensional array
 278     zr = zz[0, 0]
 279     # Initializing fire arrival time
 280     Fz = np.zeros(Fx.shape)
 281     # For each x and y
 282     for k1 in range(Fx.shape[0]):
 283         for k2 in range(Fx.shape[1]):
 284             # Approximate the vertical decision function by a piecewise polynomial (cubic spline interpolation)
 285             pz = interpolate.CubicSpline(zr, Z[k1,k2])
 286             # Compute the real roots of the the piecewise polynomial
 287             rr = pz.roots()
 288             # Just take the real roots between min(zz) and max(zz)
 289             realr = rr.real[np.logical_and(abs(rr.imag) < 1e-5, np.logical_and(rr.real > zr.min(), rr.real < zr.max()))]
 290             if len(realr) > 0:
 291                 # Take the minimum root
 292                 Fz[k1,k2] = realr.min()
 293                 # Plotting the approximated polynomial with the decision function
 294                 if plot_poly:
 295                     try:
 296                         plt.ion()
 297                         plt.plot(pz(zr),zr)
 298                         plt.plot(Z[k1,k2],zr,'+')
 299                         plt.plot(np.zeros(len(realr)),realr,'o',c='g')
 300                         plt.plot(0,Fz[k1,k2],'o',markersize=3,c='r')
 301                         plt.title('Polynomial approximation of decision_function(%f,%f,z)' % (Fx[k1,k2],Fy[k1,k2]))
 302                         plt.xlabel('decision function value')
 303                         plt.ylabel('Z')
 304                         plt.legend(['polynomial','decision values','roots','fire arrival time'])
 305                         plt.xlim([Z.min(),Z.max()])
 306                         plt.ylim([zz.min(),zz.max()])
 307                         plt.show()
 308                         plt.pause(.001)
 309                         plt.cla()
 310                     except Exception as e:
 311                         print 'Warning: something went wrong when plotting...'
 312                         print e
 313             else:
 314                 # If there is not a real root of the polynomial between zz.min() and zz.max(), just define as a Nan
 315                 Fz[k1,k2] = np.nan
 316     t_2 = time()
 317     print 'elapsed time: %ss.' % str(abs(t_2-t_1))
 318     F = [Fx,Fy,Fz]
 319
 320     return F
 321
 322 def SVM3(X, y, C=1., kgam=1., norm=True, fire_grid=None, weights=None):
 323     """
 324     3D SuperVector Machine analysis and plot
 325
 326     :param X: Training vectors, where n_samples is the number of samples and n_features is the number of features.
 327     :param y: Target values
 328     :param C: Weight to not having outliers (argument of svm.SVC class), optional
 329     :param kgam: Scalar multiplier for gamma (capture more details increasing it)
 330     :param norm: Normalize the data in the interval (0,1) in all the directions, optional
 331     :param fire_grid: The longitud and latitude grid where to have the fire arrival time
 332     :return F: tuple with (longitude grid, latitude grid, fire arrival time grid)
 333
 334     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 335     Angel Farguell (angel.farguell@gmail.com), 2019-02-20
 336     Modified version of:
 337     https://scikit-learn.org/stable/auto_examples/svm/plot_iris.html#sphx-glr-auto-examples-svm-plot-iris-py
 338     """
 339
 340     t_init = time()
 341
 342     # Plot options
 343     # plot original data
 344     plot_data = True
 345     # plot scaled data with artificial data
 346     plot_scaled = True
 347     # plot decision volume
 348     plot_decision = False
 349     # plot polynomial approximation
 350     plot_poly = False
 351     # plot full hyperplane vs detections with support vectors
 352     plot_supports = True
 353     # plot resulting fire arrival time vs detections
 354     plot_result = True
 355
 356     # Other options
 357     # number of vertical nodes per observation
 358     vN = 1
 359     # interpolate into the original fire mesh
 360     interp = False
 361     # if not Nans in the data are wanted (all Nans are going to be replaced by the maximum value)
 362     notnan = True
 363
 364     # Options better to not change
 365     # number of horizontal nodes per observation (if fire_grid==None)
 366     hN = 5
 367     # creation of under artificial lower bounds in the pre-processing
 368     artil = False
 369     # if artil = True: resolution of artificial upper bounds vertical to the fire detections
 370     hartil = .1
 371     # creation of over artificial upper bounds in the pre-processing
 372     artiu = True
 373     # if artiu = True: resolution of artificial lower bounds vertical to the ground detections
 374     hartiu = .1
 375     # creation of an artifitial mesh of down lower bounds
 376     downarti = True
 377     # if downarti = True: below min of z direction for lower bound artifitial creation
 378     dminz = .1
 379     # if downarti = True: confidence level of the artificial lower bounds
 380     confal = 100.
 381     # creation of an artifitial mesh of top upper bounds
 382     toparti = False
 383     # if toparti = True: proportion over max of z direction for upper bound artifitial creation
 384     dmaxz = .1
 385     # if toparti = True: confidence level of the artificial upper bounds
 386     confau = 100.
 387
 388     # using different weights for the data
 389     if isinstance(C,(list,tuple,np.ndarray)):
 390         using_weights = True
 391         from libsvm_weights.python.svm import svm_problem, svm_parameter
 392         from libsvm_weights.python.svmutil import svm_train
 393         from sklearn.utils import compute_class_weight
 394     else:
 395         using_weights = False
 396
 397     # Data inputs
 398     X = np.array(X).astype(float)
 399     y = np.array(y)
 400
 401     # Original data
 402     oX = np.array(X).astype(float)
 403     oy = np.array(y)
 404
 405     # Visualization of the data
 406     X0, X1, X2 = X[:, 0], X[:, 1], X[:, 2]
 407     if plot_data:
 408         try:
 409             fig = plt.figure()
 410             ax = fig.gca(projection='3d')
 411             fig.suptitle("Plotting the original data to fit")
 412             ax.scatter(X0, X1, X2, c=y, cmap=plt.cm.coolwarm, s=20, edgecolors='k', vmin=y.min(), vmax=y.max())
 413             ax.set_xlabel("Longitude")
 414             ax.set_ylabel("Latitude")
 415             ax.set_zlabel("Time (days)")
 416             plt.savefig('original_data.png')
 417         except Exception as e:
 418             print 'Warning: something went wrong when plotting...'
 419             print e
 420
 421     # Normalization of the data into [0,1]^3
 422     if norm:
 423         xmin = X0.min()
 424         xlen = X0.max() - X0.min()
 425         x0 = np.divide(X0 - xmin, xlen)
 426         ymin = X1.min()
 427         ylen = X1.max() - X1.min()
 428         x1 = np.divide(X1 - ymin, ylen)
 429         zmin = X2.min()
 430         zlen = X2.max() - X2.min()
 431         x2 = np.divide(X2 - zmin, zlen)
 432         X0, X1, X2 = x0, x1, x2
 433         X[:, 0] = X0
 434         X[:, 1] = X1
 435         X[:, 2] = X2
 436
 437     # Creation of fire and ground artificial detections
 438     if artil or artiu or toparti or downarti:
 439         # Extreme values at z direction
 440         minz = X[:, 2].min()
 441         maxz = X[:, 2].max()
 442         # Division of lower and upper bounds for data and confidence level
 443         fl = X[y==np.unique(y)[0]]
 444         fu = X[y==np.unique(y)[1]]
 445
 446         # Artifitial extensions of the lower bounds
 447         if artil:
 448             # Create artificial lower bounds
 449             flz = np.array([ np.unique(np.append(np.arange(f[2],minz,-hartil),f[2])) for f in fl ])
 450             # Definition of new ground detections after artificial detections added
 451             Xg = np.concatenate([ np.c_[(np.repeat(fl[k][0],len(flz[k])),np.repeat(fl[k][1],len(flz[k])),flz[k])] for k in range(len(flz)) ])
 452             if using_weights:
 453                 cl = C[y==np.unique(y)[0]]
 454                 Cg = np.concatenate([ np.repeat(cl[k],len(flz[k])) for k in range(len(flz)) ])
 455         else:
 456             Xg = fl
 457             if using_weights:
 458                 cl = C[y==np.unique(y)[0]]
 459                 Cg = cl
 460
 461         # Artifitial extensions of the upper bounds
 462         if artiu:
 463             # Create artificial upper bounds
 464             fuz = np.array([ np.unique(np.append(np.arange(f[2],maxz,hartiu),f[2])) for f in fu ])
 465             # Definition of new fire detections after artificial detections added
 466             Xf = np.concatenate([ np.c_[(np.repeat(fu[k][0],len(fuz[k])),np.repeat(fu[k][1],len(fuz[k])),fuz[k])] for k in range(len(fuz)) ])
 467             # Define new confidence levels
 468             if using_weights:
 469                 cu = C[y==np.unique(y)[1]]
 470                 Cf = np.concatenate([ np.repeat(cu[k],len(fuz[k])) for k in range(len(fuz)) ])
 471         else:
 472             Xf = fu
 473             if using_weights:
 474                 cu = C[y==np.unique(y)[1]]
 475                 Cf = cu
 476
 477         # Bottom artificial lower bounds
 478         if downarti:
 479             # Creation of the x,y new mesh of artificial lower bounds
 480             xn, yn = np.meshgrid(np.linspace(X[:, 0].min(), X[:, 0].max(), 20),
 481                 np.linspace(X[:, 1].min(), X[:, 1].max(), 20))
 482             # All the artificial new mesh are going to be below the data
 483             zng = np.repeat(minz-dminz,len(xn.ravel()))
 484             # Artifitial lower bounds
 485             Xga = np.c_[(xn.ravel(),yn.ravel(),zng.ravel())]
 486             # Definition of new ground detections after down artificial lower detections
 487             Xgn = np.concatenate((Xg,Xga))
 488             # Definition of new confidence level
 489             if using_weights:
 490                 Cga = np.ones(len(Xga))*confal
 491                 Cgn = np.concatenate((Cg,Cga))
 492         else:
 493             Xgn = Xg
 494             if using_weights:
 495                 Cgn = Cg
 496
 497         # Top artificial upper bounds
 498         if toparti:
 499             # Creation of the x,y new mesh of artificial upper bounds
 500             xn, yn = np.meshgrid(np.linspace(X[:, 0].min(), X[:, 0].max(), 20),
 501                                 np.linspace(X[:, 1].min(), X[:, 1].max(), 20))
 502             # All the artificial new mesh are going to be over the data
 503             znf = np.repeat(maxz+dmaxz,len(xn.ravel()))
 504             # Artifitial upper bounds
 505             Xfa = np.c_[(xn.ravel(),yn.ravel(),znf.ravel())]
 506             # Definition of new fire detections after top artificial upper detections
 507             Xfn = np.concatenate((Xf,Xfa))
 508             # Definition of new confidence level
 509             if using_weights:
 510                 Cfa = np.ones(len(Xfa))*confau
 511                 Cfn = np.concatenate((Cf,Cfa))
 512         else:
 513             Xfn = Xf
 514             if using_weights:
 515                 Cfn = Cf
 516
 517         # New definition of the training vectors
 518         X = np.concatenate((Xgn, Xfn))
 519         # New definition of the target values
 520         y = np.concatenate((np.repeat(np.unique(y)[0],len(Xgn)),np.repeat(np.unique(y)[1],len(Xfn))))
 521         # New definition of the confidence level
 522         if using_weights:
 523             C = np.concatenate((Cgn, Cfn))
 524         # New definition of each feature vector
 525         X0, X1, X2 = X[:, 0], X[:, 1], X[:, 2]
 526
 527     # Printing number of samples and features
 528     n0 = (y==np.unique(y)[0]).sum().astype(float)
 529     n1 = (y==np.unique(y)[1]).sum().astype(float)
 530     n_samples, n_features = X.shape
 531     print 'n_samples =', n_samples
 532     print 'n_samples_{-1} =', int(n0)
 533     print 'n_samples_{+1} =', int(n1)
 534     print 'n_features =', n_features
 535
 536     # Visualization of scaled data
 537     if plot_scaled:
 538         try:
 539             fig = plt.figure()
 540             ax = fig.gca(projection='3d')
 541             fig.suptitle("Plotting the data scaled to fit")
 542             ax.scatter(X0, X1, X2, c=y, cmap=plt.cm.coolwarm, s=20, edgecolors='k', vmin=y.min(), vmax=y.max())
 543             ax.set_xlabel("Longitude normalized")
 544             ax.set_ylabel("Latitude normalized")
 545             ax.set_zlabel("Time normalized")
 546             plt.savefig('scaled_data.png')
 547         except Exception as e:
 548             print 'Warning: something went wrong when plotting...'
 549             print e
 550
 551     # Reescaling gamma to include more detailed results
 552     gamma = kgam / (n_features * X.std())
 553     print 'gamma =', gamma
 554
 555     # Creating the SVM model
 556     print '>> Creating the SVM model...'
 557     sys.stdout.flush()
 558     if using_weights:
 559         # Compute class balanced weights
 560         cls, _ = np.unique(y, return_inverse=True)
 561         class_weight = compute_class_weight("balanced", cls, y)
 562         prob = svm_problem(C,y,X)
 563         arg = '-g %.15g -w%01d %.15g -w%01d %.15g -m 1000 -h 0' % (gamma, cls[0], class_weight[0],
 564                                             cls[1], class_weight[1])
 565         param = svm_parameter(arg)
 566     else:
 567         clf = svm.SVC(C=C, kernel="rbf", gamma=gamma, cache_size=1000, class_weight="balanced") # default kernel: exp(-gamma||x-x'||^2)
 568         print clf
 569
 570     # Fitting the data using Super Vector Machine technique
 571     print '>> Fitting the SVM model...'
 572     sys.stdout.flush()
 573     t_1 = time()
 574     if using_weights:
 575         clf = svm_train(prob,param)
 576     else:
 577         clf.fit(X, y)
 578     t_2 = time()
 579     print 'elapsed time: %ss.' % str(abs(t_2-t_1))
 580
 581     if not using_weights:
 582         # Check if the classification failed
 583         if clf.fit_status_:
 584             print '>> FAILED <<'
 585             print 'Failed fitting the data'
 586             sys.exit(1)
 587         print 'number of support vectors: ', clf.n_support_
 588         print 'score of trained data: ', clf.score(X,y)
 589
 590     # Creating the mesh grid to evaluate the classification
 591     print '>> Creating mesh grid to evaluate the classification...'
 592     nnodes = np.ceil(np.power(n_samples,1./n_features))
 593     if fire_grid is None:
 594         # Number of necessary nodes
 595         hnodes = hN*nnodes
 596         vnodes = vN*nnodes
 597         print 'number of horizontal nodes (%d meshgrid nodes for each observation): %d' % (hN,hnodes)
 598         print 'number of vertical nodes (%d meshgrid nodes for each observation): %d' % (vN,vnodes)
 599         # Computing resolution of the mesh to evaluate
 600         sdim = (hnodes,hnodes,vnodes)
 601         print 'grid_size = %dx%dx%d = %d' % (sdim[0],sdim[1],sdim[2],np.prod(sdim))
 602         t_1 = time()
 603         xx, yy, zz = make_meshgrid(X0, X1, X2, s=sdim)
 604         t_2 = time()
 605     else:
 606         fxlon = np.divide(fire_grid[0] - xmin, xlen)
 607         fxlat = np.divide(fire_grid[1] - ymin, ylen)
 608         it = (X2.min(),X2.max())
 609         vnodes = vN*nnodes
 610         sdim = (fxlon.shape[0],fxlon.shape[1],vnodes)
 611         print 'fire_grid_size = %dx%dx%d = %d' % (sdim + (np.prod(sdim),))
 612         t_1 = time()
 613         xx, yy, zz = make_fire_mesh(fxlon, fxlat, it, sdim[2])
 614         t_2 = time()
 615         print 'grid_created = %dx%dx%d = %d' % (zz.shape + (np.prod(zz.shape),))
 616     print 'elapsed time: %ss.' % str(abs(t_2-t_1))
 617
 618     # Computing the 2D fire arrival time, F
 619     print '>> Computing the 2D fire arrival time, F...'
 620     sys.stdout.flush()
 621     F = frontier(clf, xx, yy, zz, plot_decision=plot_decision, plot_poly=plot_poly, using_weights=using_weights)
 622
 623     print '>> Creating final results...'
 624     sys.stdout.flush()
 625     # Plotting the Separating Hyperplane of the SVM classification with the support vectors
 626     if plot_supports:
 627         try:
 628             if using_weights:
 629                 supp_ind = np.sort(clf.get_sv_indices())-1
 630                 supp_vec = X[supp_ind]
 631             else:
 632                 supp_ind = clf.support_
 633                 supp_vec = clf.support_vectors_
 634             fig = plt.figure()
 635             ax = fig.gca(projection='3d')
 636             fig.suptitle("Plotting the 3D Separating Hyperplane of an SVM")
 637             # plotting the separating hyperplane
 638             ax.plot_wireframe(F[0], F[1], F[2], color='orange')
 639             # computing the indeces where no support vectors
 640             rr = np.array(range(len(y)))
 641             ms = np.isin(rr,supp_ind)
 642             nsupp = rr[~ms]
 643             # plotting no-support vectors (smaller)
 644             ax.scatter(X0[nsupp], X1[nsupp], X2[nsupp], c=y[nsupp], cmap=plt.cm.coolwarm, s=.5, vmin=y.min(), vmax=y.max(), alpha=.1)
 645             # plotting support vectors (bigger)
 646             ax.scatter(supp_vec[:, 0], supp_vec[:, 1], supp_vec[:, 2], c=y[supp_ind], cmap=plt.cm.coolwarm, s=20, edgecolors='k', alpha=.2);
 647             ax.set_xlim(xx.min(),xx.max())
 648             ax.set_ylim(yy.min(),yy.max())
 649             ax.set_zlim(zz.min(),zz.max())
 650             ax.set_xlabel("Longitude normalized")
 651             ax.set_ylabel("Latitude normalized")
 652             ax.set_zlabel("Time normalized")
 653             plt.savefig('support.png')
 654         except Exception as e:
 655             print 'Warning: something went wrong when plotting...'
 656             print e
 657
 658     # Plot the fire arrival time resulting from the SVM classification normalized
 659     if plot_result:
 660         try:
 661             Fx, Fy, Fz = np.array(F[0]), np.array(F[1]), np.array(F[2])
 662             with np.errstate(invalid='ignore'):
 663                 Fz[Fz > X2.max()] = np.nan
 664             if notnan:
 665                 Fz[np.isnan(Fz)] = X2.max()
 666                 Fz = np.minimum(Fz, X2.max())
 667             fig = plt.figure()
 668             ax = fig.gca(projection='3d')
 669             fig.suptitle("Fire arrival time normalized")
 670             # plotting fire arrival time
 671             p = ax.plot_surface(Fx, Fy, Fz, cmap=plt.cm.coolwarm,
 672                            linewidth=0, antialiased=False)
 673             ax.set_xlim(xx.min(),xx.max())
 674             ax.set_ylim(yy.min(),yy.max())
 675             ax.set_zlim(zz.min(),zz.max())
 676             cbar = fig.colorbar(p)
 677             cbar.set_label('Fire arrival time normalized', labelpad=20, rotation=270)
 678             ax.set_xlabel("Longitude normalized")
 679             ax.set_ylabel("Latitude normalized")
 680             ax.set_zlabel("Time normalized")
 681             plt.savefig('tign_g.png')
 682         except Exception as e:
 683             print 'Warning: something went wrong when plotting...'
 684             print e
 685
 686     # Translate the result again into initial data scale
 687     if norm:
 688         f0 = F[0] * xlen + xmin
 689         f1 = F[1] * ylen + ymin
 690         f2 = F[2] * zlen + zmin
 691         FF = [f0,f1,f2]
 692
 693     # Set all the larger values at the end to be the same maximum value
 694     oX0, oX1, oX2 = oX[:, 0], oX[:, 1], oX[:, 2]
 695     FFx, FFy, FFz = FF[0], FF[1], FF[2]
 696
 697     with np.errstate(invalid='ignore'):
 698         FFz[FFz > oX2.max()] = np.nan
 699
 700         if notnan:
 701             FFz[np.isnan(FFz)] = oX2.max()
 702             FFz = np.minimum(FFz, oX2.max())
 703
 704     if (not fire_grid is None) and (interp):
 705         print '>> Interpolating the results in the fire mesh'
 706         Flon = fire_grid[0]
 707         Flat = fire_grid[1]
 708         points = np.c_[Fx.ravel(),Fy.ravel()]
 709         values = Fz.ravel()
 710         Ffire = interpolate.griddata(points,values,(Flon,Flat))
 711         FF = [Flon,Flat,Ffire]
 712     else:
 713         FF = [FFx,FFy,FFz]
 714
 715     # Plot the fire arrival time resulting from the SVM classification
 716     if plot_result:
 717         try:
 718             # Plotting the result
 719             fig = plt.figure()
 720             ax = fig.gca(projection='3d')
 721             fig.suptitle("Plotting the 3D graph function of a SVM")
 722             FFx, FFy, FFz = np.array(FF[0]), np.array(FF[1]), np.array(FF[2])
 723             # plotting original data
 724             ax.scatter(oX0, oX1, oX2, c=oy, cmap=plt.cm.coolwarm, s=2, vmin=y.min(), vmax=y.max())
 725             # plotting fire arrival time
 726             ax.plot_wireframe(FFx, FFy, FFz, color='orange', alpha=.5)
 727             ax.set_xlabel("Longitude")
 728             ax.set_ylabel("Latitude")
 729             ax.set_zlabel("Time (days)")
 730             plt.savefig('result.png')
 731         except Exception as e:
 732             print 'Warning: something went wrong when plotting...'
 733             print e
 734
 735     print '>> SUCCESS <<'
 736     t_final = time()
 737     print 'TOTAL elapsed time: %ss.' % str(abs(t_final-t_init))
 738     plt.close()
 739
 740     return FF
 741
 742
 743 if __name__ == "__main__":
 744     # Experiment to do
 745     exp = 1
 746
 747     # Defining ground and fire detections
 748     def exp1():
 749         Xg = [[0, 0, 0], [2, 2, 0], [2, 0, 0], [0, 2, 0]]
 750         Xf = [[0, 0, 1], [1, 1, 0], [2, 2, 1], [2, 0, 1], [0, 2, 1]]
 751         C = np.concatenate((10.*np.ones(len(Xg)),100.*np.ones(len(Xf))))
 752         kgam = 1.
 753         return Xg, Xf, C, kgam
 754     def exp2():
 755         Xg = [[0, 0, 0], [2, 2, 0], [2, 0, 0], [0, 2, 0],
 756             [4, 2, 0], [4, 0, 0], [2, 1, .5], [0, 1, .5],
 757             [4, 1, .5], [2, 0, .5], [2, 2, .5]]
 758         Xf = [[0, 0, 1], [1, 1, 0.25], [2, 2, 1], [2, 0, 1], [0, 2, 1], [3, 1, 0.25], [4, 2, 1], [4, 0, 1]]
 759         C = np.concatenate((np.array([50.,50.,50.,50.,50.,50.,
 760                         1000.,100.,100.,100.,100.]), 100.*np.ones(len(Xf))))
 761         kgam = 5.
 762         return Xg, Xf, C, kgam
 763
 764     # Creating the options
 765     options = {1 : exp1, 2 : exp2}
 766
 767     # Defining the option depending on the experiment
 768     Xg, Xf, C, kgam = options[exp]()
 769
 770     # Creating the data necessary to run SVM3 function
 771     X = np.concatenate((Xg, Xf))
 772     y = np.concatenate((-np.ones(len(Xg)), np.ones(len(Xf))))
 773
 774     # Running SVM classification
 775     SVM3(X,y,C=C,kgam=kgam)