def classifykNN():
    print 'Classify kNN'
    target_names = ['unacc', 'acc','good','v-good']
    df = pd.read_csv("data/cars-cleaned.txt", delimiter=",");    
    print df
    print df.dtypes
    df_y = df['accept']
    df_x = df.ix[:,:-1]

    #print df_y
    #print df_x
    train_y, test_y, train_x, test_x = train_test_split(df_y, df_x, test_size = 0.3, random_state=33)
    
    clf = KNeighborsClassifier(n_neighbors=3)
    tstart=time.time()
    model = clf.fit(train_x, train_y)
    print "training time:", round(time.time()-tstart, 3), "seconds"
    y_predictions = model.predict(test_x)
    print "Accuracy : " , model.score(test_x, test_y)
    #print y_predictions
    c_matrix = confusion_matrix(test_y,y_predictions)
    print "confusion matrix:"
    print c_matrix

    print "Nearest Neighbors probabilities"
    print model.predict_proba(test_x)
    
    plt.matshow(c_matrix)
    plt.colorbar();
    tick_marks = np.arange(len(target_names))
    plt.xticks(tick_marks, target_names, rotation=45)
    plt.yticks(tick_marks, target_names)    
    plt.ylabel('true label')
    plt.xlabel('predicted label')
    plt.show()

def classify():
    print 'Classify SVM'
    target_names = ['unacc', 'acc','good','v-good']
    df = pd.read_csv("data/cars-cleaned.txt", delimiter=",");    
    print df
    print df.dtypes
    df_y = df['accept']
    df_x = df.ix[:,:-1]
   
    train_y, test_y, train_x, test_x = train_test_split(df_y, df_x, test_size = 0.3, random_state=33)
    
    clf = svm.SVC(kernel="linear", C=0.01)
    tstart=time.time()
    model = clf.fit(train_x, train_y)
    print "training time:", round(time.time()-tstart, 3), "seconds"
    y_predictions = model.predict(test_x)
    print "Accuracy : " , model.score(test_x, test_y)
    c_matrix = confusion_matrix(test_y,y_predictions)
    print "confusion matrix:"
    print c_matrix

    plt.matshow(c_matrix)
    plt.colorbar();
    tick_marks = np.arange(len(target_names))
    plt.xticks(tick_marks, target_names, rotation=45)
    plt.yticks(tick_marks, target_names)    
    plt.ylabel('true label')
    plt.xlabel('predicted label')
    plt.show()

def export_pdf(data, output, gradient=True):
    """Exports the data as a heatmap to the file specified by output (.pdf). The gradient
	marker specifies if the color in the heatmap should be a gradient that shows the
	fold change score, or as binary red/white colors with a cutoff for maximum fold-change
	score to be considered a hit.
	"""
    bacteria = sorted(data.values()[0].keys())
    matrix = np.zeros((len(data), len(bacteria)))
    rows = sorted(
        data.keys(),
        key=lambda w: np.sum([d ** 3 for d in data[w].values()]) / (len([d for d in data[w].values() if d < 0.3]) + 1),
    )
    for wind, well in enumerate(rows):
        for bind, bacterium in enumerate(bacteria):
            if gradient:
                matrix[(wind, bind)] = data[well][bacterium]
            else:
                matrix[(wind, bind)] = 0 if data[well][bacterium] < 0.3 else 0.5
    plt.matshow(matrix, cmap=plt.get_cmap("RdYlGn"), vmin=0, vmax=1)
    plt.xticks(range(len(bacteria)), bacteria, rotation=90)
    plt.yticks(range(len(rows)), rows)
    ax = plt.gca()
    for posi in ax.spines:
        ax.spines[posi].set_color("none")
    ax.tick_params(labelcolor="k", top="off", bottom="off", left="off", right="off")
    fig = plt.gcf()
    fig.set_size_inches(10, 100)
    plt.savefig(output + ".pdf", bbox_inches="tight", dpi=200)
    plt.close()

def plot_confusion_matrix(cm):
	pl.matshow(cm)
	pl.title('Confusion matrix')
	pl.colorbar()
	pl.ylabel('True label')
	pl.xlabel('Predicted label')
	pl.show()

def generate_single_funnel_test_data( excitation_angles, emission_angles, \
                                          md_ex=0, md_fu=1, \
                                          phase_ex=0, phase_fu=0, \
                                          gr=1.0, et=1.0 ):

    ex, em = np.meshgrid( excitation_angles, emission_angles )

    alpha = 0.5 * np.arccos( .5*(((gr+2)*md_ex)-gr) )

    ph_ii_minus = phase_ex - alpha
    ph_ii_plus  = phase_ex + alpha
    
    print ph_ii_minus
    print ph_ii_plus

    Fnoet  =    np.cos( ex-ph_ii_minus )**2 * np.cos( em-ph_ii_minus )**2
    Fnoet += gr*np.cos( ex-phase_ex )**2    * np.cos( em-phase_ex )**2
    Fnoet +=    np.cos( ex-ph_ii_plus )**2  * np.cos( em-ph_ii_plus )**2
        
    Fnoet /= (2+gr)
    
    Fet   = .25 * (1+md_ex*np.cos(2*(ex-phase_ex))) \
        * (1+md_fu*np.cos(2*(em-phase_fu-phase_ex)))
    
    
    Fem = et*Fet + (1-et)*Fnoet


    import matplotlib.pyplot as plt
    plt.interactive(True)
    plt.matshow( Fem, origin='bottom' )
    plt.colorbar()

def test_initialize_at_truth():
    global alpha, beta, num_topics, num_vocab, document_lengths, \
            doc_topic, topic_word, docs, model
    alpha = 5.
    beta = 20.
    num_topics = 20
    num_vocab = 1000
    document_lengths = [100]*1000

    doc_topic, topic_word, docs = generate_synthetic(alpha,beta,
            num_topics,num_vocab,document_lengths)

    model = lda.CollapsedSampler(alpha,beta,num_topics,num_vocab)
    model.add_documents_spmat(docs)

    # initialize at truth
    model.document_topic_counts = (model.document_topic_counts.sum(1)[:,None] * doc_topic).round()
    model.topic_word_counts = (model.topic_word_counts.sum(1)[:,None] * topic_word).round()

    model.resample(1000)

    plt.matshow(topic_word[:20,:20])
    plt.title('true topic_word on first 20 words')
    plt.matshow(model.topic_word_counts[:20,:20])
    plt.title('topic_word counts on first 20 words')

def TestSVM(features, labels, silence=True):
    X_train = features[0:1600,:]
    Y_train = labels[0:1600]

    X_test = features[1600:,:]
    Y_test = labels[1600:]

    clf = SVM.SVC()
    clf.fit(X_train, Y_train)

    predictions = clf.predict(X_test)  

    error = np.mean(abs(predictions-Y_test))

    cm = confusion_matrix(Y_test, predictions)

    cm_sum = np.sum(cm, axis=1)

    cm_mean = cm.T / cm_sum
    cm_mean = cm_mean.T
    
    if silence==False:
        plt.matshow(cm_mean)
        plt.title('Confusion matrix')
        plt.colorbar()
        plt.ylabel('True label')
        plt.xlabel('Predicted label')
        plt.show()
        
    return error, cm_mean

def compare_autoencoder_outputs(imgs, model, indices=[0], img_dim=(28, 28)):
    pred = model.predict(imgs)
    for i in indices:
        tup = (imgs[i].reshape(img_dim), pred[i].reshape(img_dim))
        plt.matshow(tup[0])
        plt.matshow(tup[1])
    plt.show()

def makeConfusionMatrix(n=100):
	# import some data to play with
	trainingdata = sio.loadmat('train.mat')
	X = np.swapaxes(trainingdata['train_images'].reshape(784,60000), 0, 1)
	y = np.array(trainingdata['train_labels']).transpose()[0]
	# Split the data into a training set and a test set
	X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=n/60000.0)

	# Run classifier
	classifier = svm.SVC(kernel='linear')
	y_pred = classifier.fit(X_train, y_train).predict(X_test)
	# Compute confusion matrix
	cm = confusion_matrix(y_test, y_pred)

	print(cm)

	# Show confusion matrix in a separate window
	plt.matshow(cm)
	plt.title('Confusion matrix')
	plt.colorbar()
	plt.ylabel('True label')
	plt.xlabel('Predicted label')
	plt.savefig('matrix'+str(n))
	plt.show()
	return

def verify_gradient(f, x, eps=1e-4, tol=1e-6, **kwargs):
    """
    Compares the numerical and analytical gradients.
    """
    # print
    fval, fgrad = f(x=x, **kwargs)
    # print fval, fgrad.shape
    # bbbbbbbb
    ngrad = numerical_gradient(f=f, x=x, eps=eps, tol=tol, **kwargs)

    fgradnorm = numpy.sqrt(numpy.sum(fgrad**2))
    ngradnorm = numpy.sqrt(numpy.sum(ngrad**2))
    diffnorm = numpy.sqrt(numpy.sum((fgrad-ngrad)**2))
    # print fval.shape
    plt.matshow(fgrad)
    plt.show()
    plt.matshow(ngrad)
    plt.show()
    if fgradnorm > 0 or ngradnorm > 0:
        norm = numpy.maximum(fgradnorm, ngradnorm)
        if not (diffnorm < tol or diffnorm/norm < tol):
            raise Exception("Numerical and analytical gradients "
                            "are different: %s != %s!" % (ngrad, fgrad))
    else:
        if not (diffnorm < tol):
            raise Exception("Numerical and analytical gradients "
                            "are different: %s != %s!" % (ngrad, fgrad))
    return True

def coloc(dataR, dataG):
    #returns heatmap for colocalization based on the angle in the red- green value plot
    #regions with low intensity are filterd out

    if dataR.shape != dataG.shape:
        print('images must have same shape')
        return 0

    tol=0.02
    dataB = np.zeros(dataR.shape)
    dataB[...,0]=np.tan(dataR[...,2]/dataG[...,1])
    dataB[...,0]=np.where((dataB[...,0]-np.pi/2.)**2>tol,0,dataB[...,0])
    maskG=np.where(dataG[...,1]<np.mean(dataG[...,1]),0,dataG[...,1])
    maskR=np.where(dataR[...,2]<np.mean(dataR[...,2]),0,dataR[...,2])
    from matplotlib import pyplot
    #pyplot.matshow(maskR)
    #pyplot.show()
    #pyplot.matshow(maskG)
    #pyplot.show()

    dataB[...,0]=dataB[...,0]*maskG*maskR
    dataB=dataB*255/np.max(dataB)
    dataB = np.array(dataB, dtype=np.uint8)
    print(np.mean(dataB[...,0]))


    plot.matshow(dataB[...,0])
    plot.show()

    return dataB

def test(training_file, testing_file):
    X,y = train.process_training_examples(training_file)
    X_test, y_test = train.process_training_examples(testing_file)

    lin_svm = train.train_linear_svm(X, y)
    linear_svm_accuracy = test_with(lin_svm, X_test, y_test)
    print("LinearSVM has classification accuracy of {}%".format(100 * linear_svm_accuracy))

    rbf_svm = train.train_rbf_svm(X, y)
    rbf_svm_accuracy = test_with(rbf_svm, X_test, y_test)
    print("RBF-SVM has classification accuracy of {}%".format(100 * rbf_svm_accuracy))

    nbc = train.train_naive_bayes(X, y)
    nb_accuracy = test_with(nbc, X_test, y_test)
    print("Multinomial Naive Bayes has classification accuracy of {}%".format(100 * nb_accuracy))

    lda = train.train_lda(X, y)
    lda_accuracy = test_with(lda, X_test, y_test)
    print("LDA has classification accuracy of {}%".format(100 * lda_accuracy))

    #Print SVM confusion matrix
    y_pred = lin_svm.predict(X_test)
    cm = confusion_matrix(y_test, y_pred)
    plt.matshow(cm)
    plt.title('Confusion matrix for SVM Classification of File Fragment Types')
    #plt.colorbar()
    plt.ylabel('True File Type')
    plt.xlabel('Predicted File Type')
    plt.show()

def show_confusion_matrix(X, y):
    """docstring for show_confusion_matrix"""

    print "show matrix..."
    # Split the data into a training set and a test set
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0, test_size=.1)

    print "running classifier...."
    # Run classifier
    classifier = svm.SVC()
    y_pred = classifier.fit(X_train, y_train).predict(X_test)

    print "compute confusion matrix..."
    # Compute confusion matrix
    cm = confusion_matrix(y_test, y_pred)

    cm_sum = np.sum(cm, axis=1).T

    cm_ratio = cm / cm_sum.astype(float)[:, np.newaxis]

    print(cm_ratio)
    print cm
    print cm_sum

    print "plot matrix..."
    # Show confusion matrix in a separate window
    plt.matshow(cm_ratio)
    plt.title('Confusion matrix')
    plt.colorbar()
    plt.ylabel('True label')
    plt.xlabel('Predicted label')
    plt.show()

 def calc_field_OneDistance(self, distance, n_z_points, temp_field=None,
                            plot_matrix=False):
     if temp_field is not None:
         old_field = [self.do_static, self.do_induction, self.do_radiation]
         self.set_fields(temp_field)
         
     Z_array, dz = np.linspace(0, self.channel_height, n_z_points,
                               retstep=True)
     
     min_t = self.min_starttime + distance/C
     max_t = self.max_endtime + np.sqrt(distance*distance +
                              self.channel_height*self.channel_height)/C
     n_t_points = int((max_t-min_t)/self.dt)+1
     
     T_array = min_t + np.arange(n_t_points)*self.dt
     
     field_matrix=np.zeros((n_t_points,n_z_points))
     
     for z_i in range(n_z_points):
         self.integrand(Z_array[z_i], distance, min_t, field_matrix[:,z_i])
         
     Es=-simps(field_matrix, dx=dz)/two_pi_e0
     
     if temp_field!=None:
         self.do_static, self.do_induction, self.do_radiation=old_field
         
     if plot_matrix:
         plt.matshow(-field_matrix/two_pi_e0)
         plt.colorbar()
         plt.show()
     
     return T_array, Es

def marg_mult(model, db, samples, burn=0, filename=None, n5=False):
    """
    generates histogram for marginal distribution of posterior multiplicities.

    :param model: TorsionFitModel
    :param db: pymc.database for model
    :param samples: length of trace
    :param burn: int. number of steps to skip
    :param filename: filename for plot to save
    """
    if n5:
        multiplicities = tuple(range(1, 7))
    else:
        multiplicities = (1, 2, 3, 4, 6)
    mult_bitstring = []
    for i in model.pymc_parameters.keys():
        if i.split('_')[-1] == 'bitstring':
            mult_bitstring.append(i)
    if n5:
        histogram = np.zeros((len(mult_bitstring), samples, 5))
    else:
        histogram = np.zeros((len(mult_bitstring), samples, 5))

    for m, torsion in enumerate(mult_bitstring):
        for i, j in enumerate(db.trace('%s' % torsion)[burn:]):
            for k, l in enumerate(multiplicities):
                if 2**(l-1) & int(j):
                    histogram[m][i][k] = 1

    plt.matshow(histogram.sum(1), cmap='Blues',  extent=[0, 5, 0, 20]), plt.colorbar()
    plt.yticks([])
    plt.xlabel('multiplicity term')
    plt.ylabel('torsion')
    if filename:
        plt.savefig(filename)

def confusion(y,y_auto,l,set,method,plot=False,output=False):
  """
  Computes the confusion matrix
  """
  from sklearn.metrics import confusion_matrix
  cmat = confusion_matrix(y.values[:,0],y_auto)
  cmat = np.array(cmat,dtype=float)
  for i in range(cmat.shape[0]):
    cmat[i,:] = cmat[i,:]*100./len(np.where(y.values[:,0]==i)[0])

  if plot:
    plt.matshow(cmat,cmap=plt.cm.gray_r)
    for i in range(cmat.shape[0]):
      for j in range(cmat.shape[0]):
        if cmat[j,i] >= np.max(cmat)/2. or cmat[j,i] > 50:
          col = 'w'
        else:
          col = 'k'
        if cmat.shape[0] <= 4:
          plt.text(i,j,"%.2f"%cmat[j,i],color=col)
        else:
          plt.text(i,j,"%d"%np.around(cmat[j,i]),color=col)
    plt.title('%s set - %s'%(set,method.upper()))
    plt.xlabel('Prediction')
    plt.ylabel('Observation')
    if len(l) <= 4:
      plt.xticks(range(len(l)),l)
      plt.yticks(range(len(l)),l)
  if output:
    return cmat

def main():
    np.set_printoptions(threshold=np.nan)
    testPixelRow = 24
    testPixelCol = 17
    #obs_20120919-131142.h5,obs_20120919-131346.h5
    #create a cal file from a twilight flat
    cal = FlatCal('../../params/flatCal.dict')
    #open another twilight flat as an observation and apply a wavelength cal and the new flat cal
#    run='LICK2012'
#    obsFileName = FileName(run=run,date='20120918',tstamp='20120919-131142').flat()
#    flatCalFileName = FileName(run=run,date='20120918',tstamp='20120919-131448').flatSoln()
#    wvlCalFileName = FileName(run=run,date='20120916',tstamp='20120917-072537').calSoln()

    run = 'PAL2012'
    obsFileName = FileName(run=run,date='20121211',tstamp='20121212-140003').obs()
    flatCalFileName = FileName(run=run,date='20121210',tstamp='').flatSoln()
    wvlCalFileName = FileName(run=run,date='20121210',tstamp='20121211-133056').calSoln()
    flatCalPath = os.path.dirname(flatCalFileName)

    ob = ObsFile(obsFileName)#('obs_20120919-131142.h5')
    ob.loadWvlCalFile(wvlCalFileName)#('calsol_20120917-072537.h5')
    ob.loadFlatCalFile(flatCalFileName)#('flatsol_20120919-131142.h5')

    #plot some uncalibrated and calibrated spectra for one pixel
    fig = plt.figure()
    ax = fig.add_subplot(211)
    ax2 = fig.add_subplot(212)
    print ob.getPixelCount(testPixelRow,testPixelCol)

    #flatSpectrum,wvlBinEdges = ob.getPixelSpectrum(testPixelRow,testPixelCol,weighted=False)
    spectrum,wvlBinEdges = ob.getPixelSpectrum(testPixelRow,testPixelCol,wvlStart=cal.wvlStart,wvlStop=cal.wvlStop,wvlBinWidth=cal.wvlBinWidth,weighted=False,firstSec=0,integrationTime=-1)

    weightedSpectrum,wvlBinEdges = ob.getPixelSpectrum(testPixelRow,testPixelCol,weighted=True)
    #flatSpectrum,wvlBinEdges = cal.flatFile.getPixelSpectrum(testPixelRow,testPixelCol,wvlStart=cal.wvlStart,wvlStop=cal.wvlStop,wvlBinWidth=cal.wvlBinWidth,weighted=False,firstSec=0,integrationTime=-1)
    flatSpectrum = cal.spectra[testPixelRow,testPixelCol]
    x = wvlBinEdges[0:-1]
    ax.plot(x,cal.wvlMedians,label='median spectrum',alpha=.5)
    ax2.plot(x,cal.flatFactors[testPixelRow,testPixelCol,:],label='pixel weights',alpha=.5)
    ax2.set_title('flat weights for pixel %d,%d'%(testPixelRow,testPixelCol))
    ax.plot(x,spectrum+20,label='unweighted spectrum for pixel %d,%d'%(testPixelRow,testPixelCol),alpha=.5)
    ax.plot(x,weightedSpectrum+10,label='weighted %d,%d'%(testPixelRow,testPixelCol),alpha=.5)
    ax.plot(x,flatSpectrum+30,label='flatFile %d,%d'%(testPixelRow,testPixelCol),alpha=.5)

    ax.legend(loc='upper center', bbox_to_anchor=(0.5, 1.3),fancybox=True,ncol=3)
    plt.show()

    #display a time-flattened image of the twilight flat as it is and after using itself as it's flat cal
    #cal.flatFile.loadFlatCalFile(flatCalFileName)#('flatsol_20120919-131142.h5')
    #cal.flatFile.displaySec(weighted=True,integrationTime=-1)
    #ob.displaySec(integrationTime=-1)
    #ob.displaySec(weighted=True,integrationTime=-1)


    for idx in range(0,100,20):
        factors10 = cal.flatFactors[:,:,idx]
        plt.matshow(factors10,vmax=np.mean(factors10)+1.5*np.std(factors10))
        plt.title('Flat weights at %d'%cal.wvlBinEdges[idx])
        plt.colorbar()
        plt.savefig('plots/factors%d.png'%idx)
        plt.show()

def on_slic_superpixels():
    data = load_data('train', independent=True)
    probs = get_kraehenbuehl_pot_sp(data)
    results = eval_on_pixels(data, [np.argmax(prob, axis=-1) for prob in
                                    probs])
    plt.matshow(results['confusion'])
    plt.show()

def benchmark(clf_class, params, name):
    print("parameters:", params)
    t0 = time()
    clf = clf_class(**params).fit(X_train, y_train)
    print("done in %fs" % (time() - t0))

    if hasattr(clf, 'coef_'):
        print("Percentage of non zeros coef: %f"
              % (np.mean(clf.coef_ != 0) * 100))
    print("Predicting the outcomes of the testing set")
    t0 = time()
    pred = clf.predict(X_test)
    print("done in %fs" % (time() - t0))

    print("Classification report on test set for classifier:")
    print(clf)
    print()
    print(classification_report(y_test, pred,
                                target_names=news_test.target_names))

    cm = confusion_matrix(y_test, pred)
    print("Confusion matrix:")
    print(cm)

    # Show confusion matrix
    plt.matshow(cm)
    plt.title('Confusion matrix of the %s classifier' % name)
    plt.colorbar()

def makeConfMat(estClasses, gtClasses, outFilename, numClasses = None, plotLabels = False):
   #If not defined, find number of unique numbers in gtClasses
   if numClasses == None:
      numClasses = len(np.unique(gtClasses))

   #X axis is est, y axis is gt
   #First index is y, second is x
   confMat = np.zeros((numClasses, numClasses))
   numInstances = len(estClasses)
   for (gtIdx, estIdx) in zip(gtClasses.astype(int), estClasses.astype(int)):
      confMat[gtIdx, estIdx]  += 1

   plt.matshow(confMat)
   plt.colorbar()
   plt.xlabel("Est class")
   plt.ylabel("True class")
   plt.title("Confusion matrix")
   ax = plt.gca()
   ax.xaxis.set_ticks_position('bottom')

   #Plot labels for each field
   if plotLabels:
      for i in range(numClasses):
         for j in range(numClasses):
            labelStr = generateStatString(confMat, i, j)
            #text receives x, y coord of plot
            ax.text(j, i, labelStr, fontweight='bold',
                  horizontalalignment='center', verticalalignment='center',
                  bbox={'facecolor':'white'})

   #plt.show()
   plt.savefig(outFilename)