def moments(self):
     """Calculate covariance and correlation matrices,
     trait, genotipic and ontogenetic means"""
     zs = np.array([ind["z"] for ind in self.pop])
     xs = np.array([ind["x"] for ind in self.pop])
     ys = np.array([ind["y"] for ind in self.pop])
     bs = np.array([ind["b"] for ind in self.pop])
     ymean = ys.mean(axis=0)
     zmean = zs.mean(axis=0)
     xmean = xs.mean(axis=0)
     ymean = ys.mean(axis=0)
     bmean = bs.mean(axis=0)
     phenotipic = np.cov(zs, rowvar=0, bias=1)
     genetic = np.cov(xs, rowvar=0, bias=1)
     heridability = genetic[np.diag_indices_from(genetic)] / phenotipic[np.diag_indices_from(phenotipic)]
     corr_phenotipic = np.corrcoef(zs, rowvar=0, bias=1)
     corr_genetic = np.corrcoef(xs, rowvar=0, bias=1)
     avgP = avg_ratio(corr_phenotipic, self.modules)
     avgG = avg_ratio(corr_genetic, self.modules)
     return {
         "y.mean": ymean,
         "b.mean": bmean,
         "z.mean": zmean,
         "x.mean": xmean,
         "P": phenotipic,
         "G": genetic,
         "h2": heridability,
         "avgP": avgP,
         "avgG": avgG,
         "corrP": corr_phenotipic,
         "corrG": corr_genetic,
     }

 def on_epoch_end(self, epoch, logs=None):
     
     if self.currentEpoch % self.freq == 0:
         self.results["epochs"].append(self.currentEpoch) # add the epoch's number
         
         evaluation = "prediction (r^2)"
         resultsText = ""
 
         if self.M is not None:
             yhatKeras = self.model.predict(self.M)
             yhatKeras += self.modelEpsilon # for numerical stability
             rSQ = np.corrcoef( self.y, yhatKeras, rowvar=0)[1,0]**2  # 0.1569    
             self.results["train_accuracy"].append(rSQ) 
             resultsText += "Training " +evaluation +":" +  str(rSQ) + " / "
 
         
         if self.M_validation is not None:
             yhatKeras = self.model.predict(self.M_validation)
             yhatKeras += self.modelEpsilon # for numerical stability
             rSQ = np.corrcoef( self.y_validation, yhatKeras, rowvar=0)[1,0]**2  # 0.1569  
             self.results["test_accuracy"].append(rSQ)
             resultsText += "Test " +evaluation +":" +  str(rSQ)
             
         print(resultsText, flush = True)
         
     self.currentEpoch += 1

def plotetc(x,y,stat,season):
	cc_all = np.corrcoef(x, y['All'])[0][1]
	cc_opt = np.corrcoef(x, y['Optimal'])[0][1]
	cc_b1 = np.corrcoef(x, y['b1'])[0][1]
	cc_b2 = np.corrcoef(x, y['b2'])[0][1]
	print "Correlation coefficients for scores with {0} NAO during {1}".format(stat, season)
	print "Optimal\tb1\tb2\tAll"
	print "{0:.3f}\t{1:.3f}\t{2:.3f}\t{3:.3f}\n".format(cc_opt, cc_b1, cc_b2, cc_all)
	# matplotlib.rcParams['axes.grid'] = True
# 	matplotlib.rcParams['legend.fancybox'] = True
# 	matplotlib.rcParams['figure.figsize'] = 18, 9
# 	matplotlib.rcParams['savefig.dpi'] = 300
# 	# Set figure name and number for pdf ploting
# 	pdfName = '{0}_{1}.pdf'.format(stat, season)
# 	pp1 = PdfPages(os.path.join('/Users/andrew/Google Drive/Work/MeltModel/Output/',pdfName))
# 	fig1 = plt.figure(1)
# 	ax1 = fig1.add_subplot(111)
# 	ax1.plot(x, y['Optimal'], 'ok', label='Optimum')
# 	ax1.plot(x, y['All'], 'or', label='All')
# 	ax1.plot(x, y['b1'], 'og', label='b1')
# 	ax1.plot(x, y['b2'], 'ob', label='b2')
# 	ax1.set_xlabel("NAO")
# 	ax1.set_xlim((-3,3))
# 	ax1.set_ylabel("Score")
# 	#
# 	#ax2 = ax1.twinx()
# 	#ax2.plot(x, y['AdjOptimal'], 'ok', label='Adjusted')
# 	#ax2.set_ylabel("Adjusted Score")
# 	plt.title(stat)
# 	plt.legend(loc='upper left')
# 	pp1.savefig(bbox_inches='tight')
# 	pp1.close()
# 	plt.close()
	return 0

def main():
    if len(sys.argv) < 2:
        print("Usage: ./bootstrap.py <project_dir>")
        sys.exit(-1)

    project_dir = sys.argv[1]
    project = Project(join(project_dir, "project.json"))

    # For each bootstraped model
    for (bootstrap_number, bootstrap) in enumerate(project.bootstraps):
        boot_dir = os.path.abspath(join(project_dir, "bootstrap{}-{}".format(bootstrap_number, type(bootstrap.base_model).__name__)))
        os.makedirs(boot_dir, exist_ok=True)

        # Interior
        f, ax = plot_covariance_matrix(np.corrcoef(bootstrap.internals, rowvar=0), ["fx", "fy", "ppx", "ppy", "ps"])
        savefigure(f, join(boot_dir, "covariance-interior"))

        # For each camera
        for (cam_number, cam) in enumerate(bootstrap.extract_cameras()):
            # Scatter and distribution plots
            f, ax = plot_scatter(cam[:,0], cam[:, 1])
            savefigure(f, join(boot_dir, "cam{}-xy".format(cam_number)))
            f, ax = plot_distribution(cam[:,2])
            ax.set_xlabel("Z")
            savefigure(f, join(boot_dir, "cam{}-z".format(cam_number)))

            # X, Y, Z and angles covariances matrices
            S = np.corrcoef(cam, rowvar=0)
            f, ax = plot_covariance_matrix(S[:3, :3], ["X", "Y", "Z"])
            savefigure(f, join(boot_dir, "covariance-cam{}-pos".format(cam_number)))
            f, ax = plot_covariance_matrix(S[3:, 3:], [r"$\Omega$", "$\phi$",
                r"$\kappa$"])
            savefigure(f, join(boot_dir, "covariance-cam{}-angles".format(cam_number)))

def test_c_within_and_c_between():

  # mocking the correlation values store
  test_c_values_store = {"test_network_1":{"test_roi_1":((0,0,0),(0,0,1)), "test_roi_2":((0,1,0),(0,1,1))},
                    "test_network_2":{"test_roi_3":((1,0,0),(1,0,1)), "test_roi_4":((1,1,0),(1,1,1))}}


  data = np.zeros((2,2,2,3))
  data[0,0,0] = [1,2,3]
  data[0,0,1] = [1,2,3]
  data[0,1,0] = [-1,-2,-3]
  data[0,1,1] = [-1,-2,-3]
  data[1,0,0] = [5,4,37]
  data[1,0,1] = [5,4,37]
  data[1,1,0] = [-3,-244,-1]
  data[1,1,1] = [-3,-244,-1]


  actual = connectivity_utils.c_within(data, test_c_values_store)

  # expected values are explicitly calculated according to the rules explained in the paper
  expected = {'test_network_1':(np.corrcoef([1,2,3],[-1,-2,-3])[1,0],), 'test_network_2': (np.corrcoef([5,4,37],[-3,-244,-1])[1,0],)}

  assert_almost_equal(expected['test_network_1'], expected['test_network_1'])
  assert_almost_equal(expected['test_network_2'], expected['test_network_2'])

  actual = connectivity_utils.c_between(data, test_c_values_store)

  # expected values are explicitly calculated according to the rules explained in the paper
  expected = [np.corrcoef([1,2,3],[5,4,37])[1,0], np.corrcoef([1,2,3],[-3,-244,-1])[1,0],np.corrcoef([-1,-2,-3],[5,4,37])[1,0], np.corrcoef([-1,-2,-3],[-3,-244,-1])[1,0]]

  assert_almost_equal(np.sort(expected), np.sort(actual['test_network_1-test_network_2']))

def corr_matrix(df, similar_type=ECoreCorrType.E_CORE_TYPE_PEARS, **kwargs):
    """
    与corr_xy的区别主要是，非两两corr计算，输入参数除类别外，只有一个矩阵的输入，且输入必须为pd.DataFrame对象 or np.array
    :param df: pd.DataFrame or np.array, 之所以叫df，是因为在内部会统一转换为pd.DataFrame
    :param similar_type: ECoreCorrType, 默认值ECoreCorrType.E_CORE_TYPE_PEARS
    :return: pd.DataFrame对象
    """
    if isinstance(df, np.ndarray):
        # 把np.ndarray转DataFrame，便统一处理
        df = pd.DataFrame(df)

    if not isinstance(df, pd.DataFrame):
        raise TypeError('df must pd.DataFrame object!!!')

    # FIXME 这里不应该支持ECoreCorrType.E_CORE_TYPE_PEARS.value，只严格按照ECoreCorrType对象相等
    if similar_type == ECoreCorrType.E_CORE_TYPE_PEARS or similar_type == ECoreCorrType.E_CORE_TYPE_PEARS.value:
        # 皮尔逊相关系数计算
        corr = np.corrcoef(df.T)
    elif similar_type == ECoreCorrType.E_CORE_TYPE_SPERM or similar_type == ECoreCorrType.E_CORE_TYPE_SPERM.value:
        # 斯皮尔曼相关系数计算, 使用自定义spearmanr，不计算p_value
        corr = spearmanr(df)
    elif similar_type == ECoreCorrType.E_CORE_TYPE_SIGN or similar_type == ECoreCorrType.E_CORE_TYPE_SIGN.value:
        # 序列＋－符号相关系数, 使用np.sign取符号后，再np.corrcoef计算
        corr = np.corrcoef(np.sign(df.T))
    elif similar_type == ECoreCorrType.E_CORE_TYPE_ROLLING or similar_type == ECoreCorrType.E_CORE_TYPE_ROLLING.value:
        # pop参数window，默认使用g_rolling_corr_window
        window = kwargs.pop('window', g_rolling_corr_window)
        corr = rolling_corr(df, window=window)
    else:
        # 还是给个默认的corr计算np.corrcoef(df.T)
        corr = np.corrcoef(df.T)
    # 将计算结果的corr转换为pd.DataFrame对象，行和列索引都使用df.columns
    corr = pd.DataFrame(corr, index=df.columns, columns=df.columns)
    return corr

def corr(x, y, reps=10**4, seed=None):
    r"""
    Simulate permutation p-value for Spearman correlation coefficient

    Parameters
    ----------
    x : array-like
    y : array-like
    reps : int
    seed : RandomState instance or {None, int, RandomState instance}
        If None, the pseudorandom number generator is the RandomState
        instance used by `np.random`;
        If int, seed is the seed used by the random number generator;
        If RandomState instance, seed is the pseudorandom number generator


    Returns
    -------
    tuple
        Returns test statistic, left-sided p-value,
        right-sided p-value, two-sided p-value, simulated distribution
    """
    prng = get_prng(seed)
    tst = np.corrcoef(x, y)[0, 1]
    sims = [np.corrcoef(prng.permutation(x), y)[0, 1] for i in range(reps)]
    left_pv = np.sum(sims <= tst) / reps
    right_pv = np.sum(sims >= tst) / reps
    two_sided_pv = np.min([1, 2 * np.min([left_pv, right_pv])])
    return tst, left_pv, right_pv, two_sided_pv, sims

def correlate(name, p_index, c_index, d_mean, d_sd):
    object = bpy.data.objects[name]
    uvs = getUVs(object, p_index)
    distances = getDistancesPerParticle(model.CONNECTION_RESULTS[c_index]['d'])

    uvs2 = []
    delays = []
    for index, ds in enumerate(distances):
        samples = []
        for i in range(1):
            delay_mm = max(delayModel_delayDistribLogNormal(d_mean, d_sd), 0.1)
            uvs2.append(list(uvs[index,:]))
            delays.append(max(ds * delay_mm, 0.1))
            #samples.append( max(ds * delay_mm, 0.1) )
        #delays.append(np.mean(samples))
    delays = np.array(delays)
    
    uvs2 = np.array(uvs2)
    print(len(uvs2))
    
    corr_dist_x = np.corrcoef(uvs[:,0], distances)
    corr_dist_y = np.corrcoef(uvs[:,1], distances)    
    corr_dela_x = np.corrcoef(uvs2[:,0], delays)
    corr_dela_y = np.corrcoef(uvs2[:,1], delays)    
    print('Correlation x with distance: %f' % corr_dist_x[0][1])
    print('Correlation x with delay: %f' % corr_dela_x[0][1])
    print('Correlation y with distance: %f' % corr_dist_y[0][1])
    print('Correlation y with delay: %f' % corr_dela_y[0][1])

def test():
    data = SimData(400, 4, 15)
    cor = np.nan_to_num(np.corrcoef(data.answers, rowvar=0)) # pearson metric
    cor = np.nan_to_num(np.corrcoef(cor))
    label1 = kmeans2(cor, 6, minit='points', iter=100)[1] # hack pocet komponent
    label2 = kmeans(cor, 6, True)

    xs, ys = mds(cor, euclid=True)
    plt.subplot(1, 2, 1)
    plt.title('kmeans2 ' + str(adjusted_rand_score(data.item_concept, label1)))
    plot_clustering(
        range(cor.shape[0]), xs, ys,
        labels=label1,
        shapes=data.item_concept,
    )

    plt.subplot(1, 2, 2)
    plt.title('Kmeans ' + str(adjusted_rand_score(data.item_concept, label2)))
    plot_clustering(
        range(cor.shape[0]), xs, ys,
        labels=label2,
        shapes=data.item_concept,
    )

    plt.show()

    def traverseplot(Xin,Yin,Field,name):
        string,nodeind,leaf,label=TraverseTree(regTree,Xin,Field)
        nband=Yin.shape[1]
        k=0
        for j in leaf:
            Ytemp=Yin[nodeind[j],:]
            Xtemp=Xin[nodeind[j],:]
            Yptemp=regTreeModel.predict(Xtemp)

            fitmodel=fitModelList[k]
            if predind.ndim==1:
                Ypnewtemp=fitmodel.predict(Xtemp[:,predind.astype(int)])
            else:
                Ypnewtemp=fitmodel.predict(Xtemp[:,predind[k,:].astype(int)])

            rmse,rmse_band=RMSECal(Yptemp,Ytemp)
            rmsenew,rmse_bandnew=RMSECal(Ypnewtemp,Ytemp)

            n=nband
            f, axarr = plt.subplots(int(np.ceil(n/2)), 2,figsize=(10,12))
            for i in range(n):
                pj=int(np.ceil(i/2))
                pi=int(i%2)
                axarr[pj, pi].plot(Yptemp[:,i],Ytemp[:,i],'.')
                axarr[pj, pi].plot(Ypnewtemp[:,i],Ytemp[:,i],'.r')
                axarr[pj, pi].set_title('cluster %s,\n cc=%.3f -> %.3f, r=%.3f -> %.3f'\
                                        %(i,np.corrcoef(Yptemp[:,i],Ytemp[:,i])[0,1],np.corrcoef(Ypnewtemp[:,i],Ytemp[:,i])[0,1],
                                          rmse_band[i],rmse_bandnew[i]))
                plotFun.plot121line(axarr[pj, pi])
            f.tight_layout()
            f.suptitle(string[j],fontsize=8)
            f.subplots_adjust(top=0.9)
            plt.savefig(savedir+name+"_node%i"%j)
            plt.close()
            k=k+1

def testSVM(linkSet, patterns = None):
    cont, ncont, vectors = [], [], []
    print "\nTesting\n"
    classifier = None
    if patterns == None:
		classifier = pickle.load(open("svm-classifier", "r"))
    else:
		classifier = pickle.load(open("svm-classifier2", "r"))
    
    for link in linkSet:	
		vec = getFeatures(link, patterns)
		vectors += [vec]
		result = classifier.predict(vec)
		if result == 1.0:
			if link.endswith(".htm"):
				cont += [link + 'l']
			else:
				cont += [link]
		else:
			ncont += [link]
	
    cont = [link for link in cont if checkBoilerplate(link)]
    #ones, zeros = clusterSimilar(cont)
    
    #cont = ones
    #ncont += zeros

    print "\nCorrelation Matrix\n"
    print numpy.corrcoef(numpy.transpose(numpy.array(vectors)))
    return sorted(cont, key=lambda x: len(x)), sorted(ncont, key=lambda x: len(x))

 def test_simulate_density(self):
     # generate a rings object both from an atomic and density model and
     # ensure the correlations match
     
     num_shots = 100
     num_phi   = 1024
     
     nq = 100 # number of q vectors
     q_values = [1.0, 2.0]
     
     # atomic model
     traj = mdtraj.load(ref_file('pentagon.pdb'))        
     r1 = xray.Rings.simulate(traj, 1, q_values, num_phi, num_shots)
                               
     # density model
     grid_dimensions = [151,] * 3
     grid_spacing = 1.0 # Angstroms
     grid = structure.atomic_to_density(traj, grid_dimensions, 
                                        grid_spacing)
                                        
     r2 = xray.Rings.simulate_density(grid, grid_spacing, num_shots, 
                                      q_values, num_phi)        
     
     # compute correlations & ensure match
     c1 = r1.correlate_intra(1.0, 1.0)
     c2 = r2.correlate_intra(1.0, 1.0)
     R = np.corrcoef(c1, c2)[0,1]
     assert R > 0.95
     
     c1 = r1.correlate_intra(2.0, 2.0)
     c2 = r2.correlate_intra(2.0, 2.0)
     R = np.corrcoef(c1, c2)[0,1]
     assert R > 0.95

def _region_features_for(histone, dna, region):
    pixels0 = histone[region].ravel()
    pixels1 = dna[region].ravel()
    bin0 = pixels0 > histone.mean()
    bin1 = pixels1 > dna.mean()
    overlap = [np.corrcoef(pixels0, pixels1)[0, 1], (bin0 & bin1).mean(), (bin0 | bin1).mean()]

    spi = mh.sobel(histone, just_filter=1)
    sp = spi[mh.erode(region)]
    sdi = mh.sobel(dna, just_filter=1)
    sd = sdi[mh.erode(region)]
    sobels = [
        np.dot(sp, sp) / len(sp),
        np.abs(sp).mean(),
        np.dot(sd, sd) / len(sd),
        np.abs(sd).mean(),
        np.corrcoef(sp, sd)[0, 1],
        np.corrcoef(sp, sd)[0, 1] ** 2,
        sp.std(),
        sd.std(),
    ]

    return np.concatenate(
        [
            [region.sum()],
            haralick(histone * region, ignore_zeros=True).mean(0),
            haralick(dna * region, ignore_zeros=True).mean(0),
            overlap,
            sobels,
            haralick(mh.stretch(sdi * region), ignore_zeros=True).mean(0),
            haralick(mh.stretch(spi * region), ignore_zeros=True).mean(0),
        ]
    )

    def correlation(self):
        keys_a = set(self.gdp.keys())
        keys_b = set(self.complaint_allstate.keys())
        intersection = keys_a & keys_b
        corr_dict = {}
        ax= []
        ay = []
        for v in intersection:
            y = self.gdp[v].values()
            x = self.complaint_allstate[v].values()
            ax.append(x)
            ay.append(y)
            '''
            if(len(x) != len(y)):
                continue
            else:
                corr_dict.update({v:np.corrcoef(x,y)[0,1]})'''
        if len(ax) != len(ay):
            if(len(ax)> len(ay)):
                ay = ay[:len(ax)]
            else:
                ax = ax[:len(ay)]
        print len(flatten(ax)),len(flatten(ay))
        print np.corrcoef(flatten(ax)[:735],flatten(ay))[0,1]
        corrdict = OrderedDict(sorted(corr_dict.items(), key=itemgetter(1)))
        #print corrdict

        '''

def correl():
	for eof in [ 1, 2 ]:
		cook=[]
		glue=[]
		for model in models:
			fmod = '{0}/run1/dtred/{0}.space{1}.txt'.format(model,eof)
			fobs = '../../sst-data/detrend/ersst.space{0}.txt'.format(eof)
			eof_mod = np.loadtxt(fmod)
			print fmod
			eof_obs = np.loadtxt(fobs)
			#print eof_mod[0:40]
			idm = np.where(eof_mod == 999.)
			ido = np.where(eof_obs == 999.)
			eof_mod = np.delete(eof_mod, idm)
			eof_obs = np.delete(eof_obs, idm)
			cook.append([ model, np.corrcoef(eof_mod, eof_obs)[0, 1]] )
			
			fmodpc = '{0}/run1/dtred/PC{1}.{0}.txt'.format(model,eof)
			fobspc = '../../sst-data/detrend/PC{0}.annual.txt'.format(eof)
			pc_mod = np.loadtxt(fmodpc)
			pc_obs = np.loadtxt(fobspc)
			#print pc_mod.shape, pc_obs.shape
			glue.append([model, np.corrcoef(pc_mod, pc_obs)[0, 1]] )

	# --- Writing spatial correlation from models and Observation - EOF
		npcook = np.array(cook)
		np.savetxt('eof{0}.ar4.correl.txt'.format(eof), npcook, fmt= '%s     %6s')

	# --- Writing time correlation from models and Observation - PC
		npglue = np.array(glue)
		np.savetxt('pc{0}.ar4.correl.txt'.format(eof), npglue, fmt= '%s     %6s')

def main():
    # Define matrix dimensions
    Nobs = 1000 # Number of observation
    Nvars = 50000 # Number of variables
    Ncomp = 100 # Number of components

    # Simulated true sources
    S_true = np.random.logistic(0,1,(Ncomp,Nvars))
    # Simulated true mixing
    A_true = np.random.normal(0,1,(Nobs,Ncomp))
    # X = AS
    X = np.dot(A_true,S_true)
    # add some noise
    X = X + np.random.normal(0,1,X.shape)
    # apply ICA on X and ask for 2 components

    model = ica1(Ncomp)
    
    start = time.time()
    A,S = model.fit(X)
    total = time.time() - start
    print('total time: {}'.format(total))
    # compare if our estimates are accurate
    # correlate A with Atrue and take 
    aCorr = np.abs(np.corrcoef(A.T,A_true.T)[:Ncomp,Ncomp:]).max(axis = 0).mean()
    sCorr = np.abs(np.corrcoef(S,S_true)[:Ncomp,Ncomp:]).max(axis = 0).mean()

    print("Accuracy of estimated sources: %.2f"%sCorr)
    print("Accuracy of estimated mixing: %.2f"%aCorr)

 def test_nancorr_pearson(self):
     targ0 = np.corrcoef(self.arr_float_2d, self.arr_float1_2d)[0, 1]
     targ1 = np.corrcoef(self.arr_float_2d.flat, self.arr_float1_2d.flat)[0, 1]
     self.check_nancorr_nancov_2d(nanops.nancorr, targ0, targ1, method="pearson")
     targ0 = np.corrcoef(self.arr_float_1d, self.arr_float1_1d)[0, 1]
     targ1 = np.corrcoef(self.arr_float_1d.flat, self.arr_float1_1d.flat)[0, 1]
     self.check_nancorr_nancov_1d(nanops.nancorr, targ0, targ1, method="pearson")

    def calc_correlation(self, Xs):
        
        for X in Xs:
            pass
            print np.corrcoef(X.T)
#             np.savetxt("correlations.csv",  np.corrcoef(X.T), delimiter=",")
            print 3

def PrintResults(all_ground_truth,all_b1_output,all_b2_output,all_b3_output,all_b4_output,all_combined_output):
   print 'Error on baseline 1: ', numpy.std(all_ground_truth - all_b1_output,axis=0), \
   numpy.mean(numpy.std(all_ground_truth - all_b1_output,axis=0))
   correlation_matrix = numpy.corrcoef(all_ground_truth.T,all_b1_output.T)
   print 'cur_rho: ', correlation_matrix[0,3], correlation_matrix[1,4], correlation_matrix[2,5], \
   (correlation_matrix[0,3]+correlation_matrix[1,4]+correlation_matrix[2,5])/3

   print 'Error on baseline 2: ', numpy.std(all_ground_truth - all_b2_output,axis=0), \
   numpy.mean(numpy.std(all_ground_truth - all_b2_output,axis=0))
   correlation_matrix = numpy.corrcoef(all_ground_truth.T,all_b2_output.T)
   print 'cur_rho: ', correlation_matrix[0,3], correlation_matrix[1,4], correlation_matrix[2,5], \
   (correlation_matrix[0,3]+correlation_matrix[1,4]+correlation_matrix[2,5])/3

   print 'Error on baseline 3: ', numpy.std(all_ground_truth - all_b3_output,axis=0), \
   numpy.mean(numpy.std(all_ground_truth - all_b3_output,axis=0))
   correlation_matrix = numpy.corrcoef(all_ground_truth.T,all_b3_output.T)
   print 'cur_rho: ', correlation_matrix[0,3], correlation_matrix[1,4], correlation_matrix[2,5], \
   (correlation_matrix[0,3]+correlation_matrix[1,4]+correlation_matrix[2,5])/3

   print 'Error on baseline 4: ', numpy.std(all_ground_truth - all_b4_output,axis=0), \
   numpy.mean(numpy.std(all_ground_truth - all_b4_output,axis=0))
   correlation_matrix = numpy.corrcoef(all_ground_truth.T,all_b4_output.T)
   print 'cur_rho: ', correlation_matrix[0,3], correlation_matrix[1,4], correlation_matrix[2,5], \
   (correlation_matrix[0,3]+correlation_matrix[1,4]+correlation_matrix[2,5])/3

   print 'Error on combined: ', numpy.std(all_ground_truth - all_combined_output,axis=0), \
   numpy.mean(numpy.std(all_ground_truth - all_combined_output,axis=0))
   correlation_matrix = numpy.corrcoef(all_ground_truth.T,all_combined_output.T)
   print 'cur_rho: ', correlation_matrix[0,3], correlation_matrix[1,4], correlation_matrix[2,5], \
   (correlation_matrix[0,3]+correlation_matrix[1,4]+correlation_matrix[2,5])/3

def corr(x, y, reps=10**4, prng=None):
    """
    Simulate permutation p-value for Spearman correlation coefficient

    Parameters
    ----------
    x : array-like
    y : array-like
    reps : int
    prng : RandomState instance or None, optional (default=None)
        If RandomState instance, prng is the pseudorandom number generator;
        If None, the pseudorandom number generator is the RandomState
        instance used by `np.random`.

    Returns
    -------
    tuple
        Returns test statistic, left-sided p-value,
        right-sided p-value, two-sided p-value, simulated distribution
    """
    if prng is None:
        prng = RandomState()
    tst = np.corrcoef(x, y)[0, 1]
    sims = [np.corrcoef(prng.permutation(x), y)[0, 1] for i in range(reps)]
    left_pv = np.sum(sims <= tst)/reps
    right_pv = np.sum(sims >= tst)/reps
    two_sided_pv = np.sum(np.abs(sims) >= np.abs(tst))/reps
    return tst, left_pv, right_pv, two_sided_pv, sims