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