def similarness(image1,image2):
    """
Return the correlation distance be1tween the histograms. This is 'normalized' so that
1 is a perfect match while -1 is a complete mismatch and 0 is no match.
"""
    # Open and resize images to 200x200
    i1 = Image.open(image1).resize((200,200))
    i2 = Image.open(image2).resize((200,200))

    # Get histogram and seperate into RGB channels
    i1hist = numpy.array(i1.histogram()).astype('float32')
    i1r, i1b, i1g = i1hist[0:256], i1hist[256:256*2], i1hist[256*2:]
    # Re bin the histogram from 256 bins to 48 for each channel
    i1rh = numpy.array([sum(i1r[i*16:16*(i+1)]) for i in range(16)]).astype('float32')
    i1bh = numpy.array([sum(i1b[i*16:16*(i+1)]) for i in range(16)]).astype('float32')
    i1gh = numpy.array([sum(i1g[i*16:16*(i+1)]) for i in range(16)]).astype('float32')
    # Combine all the channels back into one array
    i1histbin = numpy.ravel([i1rh, i1bh, i1gh]).astype('float32')

    # Same steps for the second image
    i2hist = numpy.array(i2.histogram()).astype('float32')
    i2r, i2b, i2g = i2hist[0:256], i2hist[256:256*2], i2hist[256*2:]
    i2rh = numpy.array([sum(i2r[i*16:16*(i+1)]) for i in range(16)]).astype('float32')
    i2bh = numpy.array([sum(i2b[i*16:16*(i+1)]) for i in range(16)]).astype('float32')
    i2gh = numpy.array([sum(i2g[i*16:16*(i+1)]) for i in range(16)]).astype('float32')
    i2histbin = numpy.ravel([i2rh, i2bh, i2gh]).astype('float32')

    return cv2.compareHist(i1histbin, i2histbin, 0)

def kdtree_fast(latvar,lonvar,lat0,lon0):
    '''
    :param latvar:
    :param lonvar:
    :param lat0:
    :param lon0:
    :return:
    '''
    rad_factor = pi/180.0 # for trignometry, need angles in radians
    # Read latitude and longitude from file into numpy arrays
    latvals = latvar[:] * rad_factor
    lonvals = lonvar[:] * rad_factor
    ny,nx = latvals.shape
    clat,clon = cos(latvals),cos(lonvals)
    slat,slon = sin(latvals),sin(lonvals)
    # Build kd-tree from big arrays of 3D coordinates
    triples = list(zip(ravel(clat*clon), ravel(clat*slon), ravel(slat)))
    kdt = cKDTree(triples)
    lat0_rad = lat0 * rad_factor
    lon0_rad = lon0 * rad_factor
    clat0,clon0 = cos(lat0_rad),cos(lon0_rad)
    slat0,slon0 = sin(lat0_rad),sin(lon0_rad)
    dist_sq_min, minindex_1d = kdt.query([clat0*clon0, clat0*slon0, slat0])
    iy_min, ix_min = unravel_index(minindex_1d, latvals.shape)
    return iy_min,ix_min

 def counts(self, a):
     """Returns array containing counts of each item in a.
     
     For example, on the enumeration 'UCAG', the sequence 'CCUG' would
     return the array [1,2,0,1] reflecting one count for the first item
     in the enumeration ('U'), two counts for the second item ('C'), no
     counts for the third item ('A'), and one count for the last item ('G').
     
     The result will always be a vector of Int with length equal to
     the length of the enumeration. We return Int and non an unsigned
     type because it's common to subtract counts, which produces surprising
     results on unit types (i.e. wrapraround to maxint) unless the type
     is explicitly coerced by the user.
     
     Sliently ignores any unrecognized indices, e.g. if your enumeration
     contains 'TCAG' and you get an 'X', the 'X' will be ignored because
     it has no index in the enumeration.
     """
     try:
         data = ravel(a)
     except ValueError:  #ravel failed; try coercing to array
         try:
             data = ravel(array(a))
         except ValueError: #try mapping to string
             data = ravel(array(map(str, a)))
     return sum(asarray(self._allowed_range == data, Int), axis=-1)

def create_edisp(event_class, event_type, erec, egy, cth):
    """Create an array of energy response values versus energy and
    inclination angle.

    Parameters
    ----------
    egy : `~numpy.ndarray`
        Energy in MeV.

    cth : `~numpy.ndarray`
        Cosine of the incidence angle.

    """
    irf = create_irf(event_class, event_type)
    theta = np.degrees(np.arccos(cth))
    v = np.zeros((len(erec), len(egy), len(cth)))
    m = (erec[:,None] / egy[None,:] < 3.0) & (erec[:,None] / egy[None,:] > 0.33333)
    #    m |= ((erec[:,None] / egy[None,:] < 3.0) &
    #          (erec[:,None] / egy[None,:] > 0.5) & (egy[None,:] < 10**2.5))    
    m = np.broadcast_to(m[:,:,None], v.shape)

    try:    
        x = np.ones(v.shape)*erec[:,None,None]
        y = np.ones(v.shape)*egy[None,:,None]
        z = np.ones(v.shape)*theta[None,None,:]
        v[m] = irf.edisp().value(np.ravel(x[m]), np.ravel(y[m]), np.ravel(z[m]), 0.0)
    except:
        for i, x in enumerate(egy):
            for j, y in enumerate(theta):
                m = (erec / x < 3.0) & (erec / x > 0.333)
                v[m, i, j] = irf.edisp().value(erec[m], x, y, 0.0)
        
    return v

    def __eq__(self, other):
        if not isinstance(other, DenseMatrix) or self.numRows != other.numRows or self.numCols != other.numCols:
            return False

        self_values = np.ravel(self.toArray(), order="F")
        other_values = np.ravel(other.toArray(), order="F")
        return all(self_values == other_values)

    def __init__(self, con_id=None, onset=None, amplitude=None):
        """
        Parameters
        ----------
        con_id: array of shape (n_events), type = string, optional
               identifier of the events
        onset: array of shape (n_events), type = float, optional,
               onset time (in s.) of the events
        amplitude: array of shape (n_events), type = float, optional,
                   amplitude of the events (if applicable)
        """
        self.con_id = con_id
        self.onset = onset
        self.amplitude = amplitude
        self.n_event = 0
        if con_id is not None:
            self.n_events = len(con_id)
            try:
                # this is only for backward compatibility:
                # if con_id were integers, they become a string
                self.con_id = np.array(["c" + str(int(float(c))) for c in con_id])
            except:
                self.con_id = np.ravel(np.array(con_id)).astype("str")

        if onset is not None:
            if len(onset) != self.n_events:
                raise ValueError("inconsistent definition of ids and onsets")
            self.onset = np.ravel(np.array(onset)).astype(np.float)
        if amplitude is not None:
            if len(amplitude) != self.n_events:
                raise ValueError("inconsistent definition of amplitude")
            self.amplitude = np.ravel(np.array(amplitude))
        self.type = "event"
        self.n_conditions = len(np.unique(self.con_id))

 def _check_hessian(self):
     if self.ff.system.cell.nvec != 0:
         # external rotations should be implemented properly for periodic systems.
         # 1D -> one external rotation, 2D and 3D -> no external rotation
         raise NotImplementedError('The hessian test is only working for isolated systems')
     # compute hessian
     hessian = estimate_cart_hessian(self.ff)
     # construct basis of external/internal degrees (rows)
     x, y, z = self.ff.system.pos.T
     natom = self.ff.system.natom
     ext_basis = np.array([
         [1.0, 0.0, 0.0]*natom,
         [0.0, 1.0, 0.0]*natom,
         [0.0, 0.0, 1.0]*natom,
         # TODO: this assumes geometry is centered for good conditioning
         np.ravel(np.array([np.zeros(natom), z, -y]).T),
         np.ravel(np.array([-z, np.zeros(natom), x]).T),
         np.ravel(np.array([y, -x, np.zeros(natom)]).T),
     ]).T
     u, s, vt = np.linalg.svd(ext_basis, full_matrices=True)
     rank = (s > s.max()*1e-10).sum() # for linear and
     int_basis = u[:,rank:]
     # project hessian
     int_hessian = np.dot(int_basis.T, np.dot(hessian, int_basis))
     evals = np.linalg.eigvalsh(int_hessian)
     self.num_neg_evals = (evals < 0).sum()
     # call tamkin as double check
     import tamkin
     system = self.ff.system
     mol = tamkin.Molecule(system.numbers, system.pos, system.masses, self.energy, self.gpos, hessian)
     nma = tamkin.NMA(mol, tamkin.ConstrainExt())
     invcm = lightspeed/centimeter
     #print nma.freqs/invcm
     self.num_neg_evals = (nma.freqs < 0).sum()

def make_kernel_grid(freq, kernel_size, n_pix, placement_grid):
    """
    make_kernel_grid(freq,kernel_size,n_pix,placement_grid)
    
    freq ~ cyc/n_pix
    
    kernel_size ~ pix. the fwhm of the gaussian envelope, effectively the kernel radius.
    
    n_pix ~ pixels per side of square image
    
    placement_grid = (X,Y) grid of kernel centers, as from meshgrid
   
    return:
      kernel_set  = 3D numpy array of complex ripple filters ~ [number_of_filters] x [n_pix] x [n_pix]
    
    """
    iter_x = np.ravel(placement_grid[0])
    iter_y = np.ravel(placement_grid[1])

    kernel_set = np.zeros((len(iter_x), n_pix, n_pix)).astype(complex)
    count = 0
    print "constructing %d filters" % (len(iter_x))
    for x, y in zip(iter_x, iter_y):
        kernel_set[count, :, :] = complex_ripple_filter(freq, (x, y), kernel_size, n_pix)
        count += 1
    return kernel_set

    def on_epoch_end(self, epoch, logs={}):
        model.save_weights(weightSavePath + "bestWeights_regressMOS_smallNetwork_latestModel.h5",overwrite=True)
        logging.info(" -- Epoch "+str(epoch)+" done, loss : "+ str(logs.get('loss')))

        predictedScoresVal = np.ravel(model.predict(valData,batch_size=batchSize))
        predictedScoresTest = np.ravel(model.predict(testData,batch_size=batchSize))
        sroccVal = scipy.stats.spearmanr(predictedScoresVal, valLabels)
        plccVal =  scipy.stats.pearsonr(predictedScoresVal, valLabels)
        sroccTest = scipy.stats.spearmanr(predictedScoresTest, testLabels)
        plccTest =  scipy.stats.pearsonr(predictedScoresTest, testLabels)
        t_str_val = '\nSpearman corr for validation set is ' + str(sroccVal[0]) + '\nPearson corr for validation set is '+ str(plccVal[0]) + '\nMean absolute error for validation set is ' + str(np.mean(np.abs(predictedScoresVal-valLabels)))
        t_str_test = '\nSpearman corr for test set is ' + str(sroccTest[0]) + '\nPearson corr for test set is '+ str(plccTest[0]) + '\nMean absolute error for test set is ' + str(np.mean(np.abs(predictedScoresTest-testLabels)))
        print t_str_val
        print t_str_test

        mean_corr = sroccVal[0] + plccVal[0]
        if mean_corr > self.best_mean_corr:
            self.best_mean_corr = mean_corr
            model.save_weights(weightSavePath + "bestWeights_regressMOS_smallNetwork_bestCorr.h5",overwrite=True)
            printing("Best correlation loss model saved at Epoch " + str(epoch) + "\n")

        self.metric.append(logs.get("val_loss"))
        if epoch % 5 == 0:
            model.optimizer.lr.set_value(round(Decimal(0.8*model.optimizer.lr.get_value()),8))
            learningRate = model.optimizer.lr.get_value()
            printing("")
            printing("The current learning rate is: " + str(learningRate))

    def _binopt(self, other, op, in_shape=None, out_shape=None):
        """apply the binary operation fn to two sparse matrices"""

        # ideally we'd take the GCDs of the blocksize dimensions
        # and explode self and other to match
        other = self.__class__(other, blocksize=self.blocksize)

        # e.g. bsr_plus_bsr, etc.
        fn = getattr(sparsetools, self.format + op + self.format)

        R,C = self.blocksize

        max_bnnz = len(self.data) + len(other.data)
        indptr = np.empty_like(self.indptr)
        indices = np.empty(max_bnnz, dtype=np.intc)
        data = np.empty(R*C*max_bnnz, dtype=upcast(self.dtype,other.dtype))

        fn(self.shape[0]//R, self.shape[1]//C, R, C,
                self.indptr, self.indices, np.ravel(self.data),
                other.indptr, other.indices, np.ravel(other.data),
                indptr, indices, data)

        actual_bnnz = indptr[-1]
        indices = indices[:actual_bnnz]
        data = data[:R*C*actual_bnnz]

        if actual_bnnz < max_bnnz/2:
            indices = indices.copy()
            data = data.copy()

        data = data.reshape(-1,R,C)

        return self.__class__((data, indices, indptr), shape=self.shape)

def cost(params, Y, R, num_features, lambdas):
    Y = np.matrix(Y)  # (1682, 943)
    R = np.matrix(R)  # (1682, 943)
    num_movies = Y.shape[0]
    num_users = Y.shape[1]

    # reshape the parameter array into parameter matrices
    X = np.matrix(np.reshape(params[:num_movies * num_features], (num_movies, num_features)))  # (1682, 10)
    Theta = np.matrix(np.reshape(params[num_movies * num_features:], (num_users, num_features)))  # (943, 10)

    # initializations
    J = 0
    X_grad = np.zeros(X.shape)  # (1682, 10)
    Theta_grad = np.zeros(Theta.shape)  # (943, 10)


    # compute the cost
    error = np.multiply((X * Theta.T) - Y, R)  # (1682, 943)
    squared_error = np.power(error, 2)  # (1682, 943)
    J = (1. / 2) * np.sum(squared_error)

    # add the cost regularization
    J = J + ((lambdas / 2) * np.sum(np.power(Theta, 2)))
    J = J + ((lambdas / 2) * np.sum(np.power(X, 2)))

    # calculate the gradients with regularization
    X_grad = (error * Theta) + (lambdas * X)
    Theta_grad = (error.T * X) + (lambdas * Theta)

    # unravel the gradient matrices into a single array
    grad = np.concatenate((np.ravel(X_grad), np.ravel(Theta_grad)))

    return J, grad

    def objective_function(self, fps, fjac=None, **kwargs):
        """
        Function to minimize.

        Parameters
        ----------
        fps : list
            parameters returned by the fitter
        fjac : None or list
            parameters for which to compute the jacobian
        args : list
            [model, [weights], [input coordinates]]
        """
        status = 0
        model = kwargs['model']
        weights = kwargs['weights']
        model.parameters = fps
        meas = kwargs['err']
        if 'y' in kwargs:
            args = (kwargs['x'], kwargs['y'])
        else:
            args = (kwargs['x'],)
        r = [status]
        if weights is None:
            residuals = np.ravel(model(*args) - meas)
            r.append(residuals)
        else:
            residuals = np.ravel(weights * (model(*args) - meas))
            r.append(residuals)
        if fjac is not None:
            args = args + (meas,)
            fderiv = np.array(self._wrap_deriv(fps, model, weights, *args))
            r.append(fderiv)
        return r

def test_cvxopt():
    mycvxopt.solvers().qp(0,0,0,0,0,0)
    path = '/Users/Admin/Dropbox/ml/MachineLearning_CS6140'
    with open(os.path.join(path, 'cvxopt.pkl'), 'rb') as f:
        arr = pickle.load(f)
    print 'pickle loaded'
    P = arr[0]
    q = arr[1]
    G = arr[2]
    h = arr[3]
    A = arr[4]
    b = arr[5]
    print 'input assigned'
    #     pcost       dcost       gap    pres   dres
    #0: -6.3339e+03 -5.5410e+05  2e+06  2e+00  2e-14
    #1:  5.8332e+02 -3.1277e+05  5e+05  2e-01  2e-14
    #2:  1.3585e+03 -1.3003e+05  2e+05  7e-02  2e-14
    #return np.ravel(solution['x'])
    with open(os.path.join(path, 'cvxopt_solution.pkl'), 'rb') as f:
        solution = pickle.load(f)
    print 'solution pickle loaded'

    mysolution = cvxopt.solvers.qp(P, q, G, h, A, b)
    print 'convex optimizer solved'
    if np.allclose(np.ravel(mysolution['x']), np.ravel(solution['x'])):
        print 'EQUAL!!!'
    else:
        print 'WROng!!!'

            def producer():

                try:
                    # Load the data from HDF5 file
                    with h5py.File(self.hdf5_file, "r") as hf:
                        num_chan, height, width = self.X_shape[-3:]
                        # Select start_idx at random for the batch
                        idx_start = np.random.randint(0, self.X_shape[0] - self.batch_size)
                        idx_end = idx_start + self.batch_size
                        # Get X and y
                        X_batch_color = hf["%s_lab_data" % self.dset][idx_start: idx_end, :, :, :]

                        X_batch_black = X_batch_color[:, :1, :, :]
                        X_batch_ab = X_batch_color[:, 1:, :, :]
                        npts, c, h, w = X_batch_ab.shape
                        X_a = np.ravel(X_batch_ab[:, 0, :, :])
                        X_b = np.ravel(X_batch_ab[:, 1, :, :])
                        X_batch_ab = np.vstack((X_a, X_b)).T

                        Y_batch = self.get_soft_encoding(X_batch_ab, nn_finder, nb_q)
                        # Add prior weight to Y_batch
                        idx_max = np.argmax(Y_batch, axis=1)
                        weights = prior_factor[idx_max].reshape(Y_batch.shape[0], 1)
                        Y_batch = np.concatenate((Y_batch, weights), axis=1)
                        # # Reshape Y_batch
                        Y_batch = Y_batch.reshape((npts, h, w, nb_q + 1))

                        # Put the data in a queue
                        queue.put((X_batch_black, X_batch_color, Y_batch))
                except:
                    print("Nothing here")

 def _lf_acc(self, subset, lf_idx):
   gt = self.gt._gt_vec
   pred = np.ravel(self.lf_matrix.tocsc()[:,lf_idx].todense())
   has_label = np.where(pred != 0)
   has_gt = np.where(gt != 0)
   # Get labels/gt for candidates in dev set, with label, with gt
   gd_idxs = np.intersect1d(has_label, subset)
   gd_idxs = np.intersect1d(has_gt, gd_idxs)
   gt = np.ravel(gt[gd_idxs])
   pred_sub = np.ravel(pred[gd_idxs])
   n_neg = np.sum(pred_sub == -1)
   n_pos = np.sum(pred_sub == 1)
   if np.sum(pred == -1) == 0:
     neg_acc = -1
   elif n_neg == 0:
     neg_acc = 0
   else:
     neg_acc = float(np.sum((pred_sub == -1) * (gt == -1))) / n_neg
   if np.sum(pred == 1) == 0:
     pos_acc = -1
   elif n_pos == 0:
     pos_acc = 0
   else: 
     pos_acc = float(np.sum((pred_sub == 1) * (gt == 1))) / n_pos
   return (pos_acc, n_pos, neg_acc, n_neg)

def plot_3d_covariance(mean, cov):
    o,w,h = covariance_ellipse(cov,3)
    # rotate width and height to x,y axis
    wx = abs(w*np.cos(o) + h*np.sin(o))*1.2
    wy = abs(h*np.cos(o) - w*np.sin(o))*1.2
    # scale w 拿来重用
    if wx > wy:
        w = wx
    else:
        w = wy

    minx = mean[0] - w
    maxx = mean[0] + w
    miny = mean[1] - w
    maxy = mean[1] + w

    xs = np.arange(minx, maxx, (maxx-minx)/40.)
    ys = np.arange(miny, maxy, (maxy-miny)/40.)
    xv, yv = np.meshgrid (xs, ys)

    zs = np.array([100.* stats.multivariate_normal.pdf(np.array([x,y]),mean,cov) \
                   for x,y in zip(np.ravel(xv), np.ravel(yv))])
    zv = zs.reshape(xv.shape)

    ax = plt.figure().add_subplot(111, projection='3d')
    ax.plot_surface(xv, yv, zv, rstride=1, cstride=1, cmap=cm.autumn)

    ax.set_xlabel('X')
    ax.set_ylabel('Y')

    ax.contour(xv, yv, zv, zdir='x', offset=minx-1, cmap=cm.autumn)
    ax.contour(xv, yv, zv, zdir='y', offset=maxy, cmap=cm.BuGn)

    def add_ij(self):
        self.imat,self.jmat = np.meshgrid(np.arange(self.i2-self.i1),
                                          np.arange(self.j2-self.j1))
        self.kdijvec = np.vstack((np.ravel(self.imat),
                                  np.ravel(self.jmat))).T

        self._ijvec = np.arange(np.prod(self.imat.shape))

    def _returnXY(self):
        """
        Returns gridded points as a vector
        """
        X,Y = np.meshgrid(self.xgrd,self.ygrd)

        return np.column_stack((np.ravel(X),np.ravel(Y)))

    def reproject(self, nj_obj, field):
        """Reproject a field of another njord inst. to the current grid"""
        if not hasattr(self,'nj_ivec'):
            self.add_njijvec(nj_obj)
        field = getattr(nj_obj, field) if type(field) is str else field
        
        if hasattr(nj_obj, 'tvec') and (len(nj_obj.tvec) == field.shape[0]):
            newfield = np.zeros(nj_obj.tvec.shape + self.llat.shape)
            for tpos in range(len(nj_obj.tvec)):
                newfield[tpos,:,:] = self.reproject(nj_obj, field[tpos,...])
            return newfield

        di = self.i2 - self.i1
        dj = self.j2 - self.j1
        xy = np.vstack((self.nj_jvec, self.nj_ivec))
        if type(field) == str:
            weights = np.ravel(nj_obj.__dict__[field])[self.nj_mask]
        else:
            weights = np.ravel(field)[self.nj_mask]
        mask = ~np.isnan(weights) 
        flat_coord = np.ravel_multi_index(xy[:,mask],(dj, di))
        sums = np.bincount(flat_coord, weights[mask])
        cnts = np.bincount(flat_coord)
        fld = np.zeros((dj, di)) * np.nan
        fld.flat[:len(sums)] = sums.astype(np.float)/cnts
        try:
            self.add_landmask()
            fld[self.landmask] = np.nan
        except:
            print "Couldn't load landmask for %s" % self.projname
        return fld

    def predict(self, X, return_std=False):
        """
        Perform classification on test vectors X.

        Parameters
        ----------
        X : {array-like, object with finite length or shape}
            Training data, requires length = n_samples

        return_std : boolean, optional
            Whether to return the standard deviation of posterior prediction.
            All zeros in this case.

        Returns
        -------
        y : array, shape = [n_samples]  or [n_samples, n_outputs]
            Predicted target values for X.

        y_std : array, shape = [n_samples]  or [n_samples, n_outputs]
            Standard deviation of predictive distribution of query points.
        """
        check_is_fitted(self, "constant_")
        n_samples = _num_samples(X)

        y = np.full((n_samples, self.n_outputs_), self.constant_,
                    dtype=np.array(self.constant_).dtype)
        y_std = np.zeros((n_samples, self.n_outputs_))

        if self.n_outputs_ == 1 and not self.output_2d_:
            y = np.ravel(y)
            y_std = np.ravel(y_std)

        return (y, y_std) if return_std else y