def column_stack(tup):
    """ Stack 1D arrays as columns into a 2D array

        Description:
            Take a sequence of 1D arrays and stack them as columns
            to make a single 2D array.  All arrays in the sequence
            must have the same first dimension.  2D arrays are
            stacked as-is, just like with hstack.  1D arrays are turned
            into 2D columns first.

        Arguments:
            tup -- sequence of 1D or 2D arrays.  All arrays must have the same
                   first dimension.
        Examples:
            >>> import numpy
            >>> a = array((1,2,3))
            >>> b = array((2,3,4))
            >>> numpy.column_stack((a,b))
            array([[1, 2],
                   [2, 3],
                   [3, 4]])

    """
    arrays = []
    for v in tup:
        arr = array(v,copy=False,subok=True)
        if arr.ndim < 2:
            arr = array(arr,copy=False,subok=True,ndmin=2).T
        arrays.append(arr)
    return _nx.concatenate(arrays,1)

 def testDistanceBetweenCurves(self):
   l1 = {'save_xd': array([[0.5,1,0], [1.5,1,0]]), 'lamb':array([0.0, 1.0])}
   l2 = {'save_xd': array([[0,0,0], [1,0,0], [2,0,0]])}
   x = arange(-1,1,0.005)
   line = array(zip(x,x,x)) #not actually needed for calcualtion, but dummy argument to residuals_cal for now
   residuals_calc = LPCResiduals(line, tube_radius = 0.2)
   dist = residuals_calc._distanceBetweenCurves(l1,l2)

def _replace_zero_by_x_arrays(sub_arys):
    for i in range(len(sub_arys)):
        if len(_nx.shape(sub_arys[i])) == 0:
            sub_arys[i] = _nx.array([])
        elif _nx.sometrue(_nx.equal(_nx.shape(sub_arys[i]),0)):
            sub_arys[i] = _nx.array([])
    return sub_arys

 def test_matrix_product(self):
     A = array( [[1,1],
                 [0,1]] )
     B = array( [[2,0],
                 [3,4]] )
     
     C = dot(A,B)
     numpy.testing.assert_array_equal(C,array([[5, 4],
                                               [3, 4]]))

 def test_elementwise_product(self):
     A = array( [[1,1],
                 [0,1]] )
     B = array( [[2,0],
                 [3,4]] )
     C = A*B                         # elementwise product
     
     numpy.testing.assert_array_equal(C, array([[2, 0],
                                                [0, 4]]))

def loadClassifierNormNum():
    # 加载数据
    datingDataMat, datingLabels = kNN.filedata2matrix('datingTestSet2.txt')
    normMat, ranges, minVals = kNN.autoNorm(datingDataMat)
    
    # 打印数据
    print normMat
    print ranges
    print minVals
    
    # 图形化显示数据
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(normMat[:, 0], normMat[:, 1], 15.0 * array(datingLabels), 15.0 * array(datingLabels))
    plt.show()

    def __init__(self,filepath):
        p = Path(filepath).resolve()
        self._filepath = p
        self._imgPtr = Image.open(str(p))
        self._imgArr = [array(self._imgPtr)] # cached copies of image levels

        self._levels = [LevelInfo(self._imgArr[0].shape[0],self._imgArr[0].shape[1],1,1,1)]

 def calculatePurity(self, curves, data_range, voxel_to_pdg_dictionary):
   '''NB - self._residuals_runner should have had calculateResiduals method called with calc_residuals = True before calling this method
   '''
   hit_tuples = voxel_to_pdg_dictionary.keys()
   if data_range is None:
     data_range = 1.0
   #rescales the truth data if necessary
   hits = array([[h[0], h[1], h[2]] for h in hit_tuples]) / data_range
   self._residuals_runner.setDataPoints(hits)
   self._residuals_runner.setLpcCurves(curves)
   self._residuals_runner.calculateResiduals(True, False, False)
   residuals = self._residuals_runner.getResiduals()
   tau_range = self._residuals_runner.getTauRange()
   purity = {}
   for tau in tau_range:
     pdg_code_energy_deposition = []
     for i in range(len(curves)):
       d = defaultdict(float)
       #NB This relies upon the the keys of voxel_to_pdg_dictionary (hit_tuple) being stored in the same order as the original data points in
       #lpc.Xi (as read in from file in lpcAnalyser), so the coverage_indices correctly index identify the hit.
       hit_labels = [voxel_to_pdg_dictionary[hit_tuples[i]] for i in residuals['curve_residuals'][i]['coverage_indices'][tau]]
       flattened_hit_labels = [pdg_code_weight for pdg_code_weight_list in hit_labels for pdg_code_weight in pdg_code_weight_list]
       for pdg_code_weight in flattened_hit_labels:
         d[pdg_code_weight[0]] += pdg_code_weight[1]
       pdg_code_energy_deposition.append(d)
     purity[tau] = pdg_code_energy_deposition
   return purity

    def test_linkage_to_d3_4_observations(self):
        Z = array([[ 1.        ,  3.        ,  0.45015331,  2.        ],   # arr[0], cluster4
            [ 0.        ,  2.        ,  1.29504919,  2.        ],   # arr[1], cluster5
            [ 4.        ,  5.        ,  1.55180264,  4.        ]])  # arr[2], cluster6

        expected = {
            "name": "cluster6",
            "children": [
                    {
                    "name": "cluster4",
                    "children": [
                            {"name": "cluster1", "size": 10},
                            {"name": "cluster3", "size": 10},
                    ]
                },
                {
                "name": "cluster5",
                "children": [
                        {"name": "cluster0", "size": 10},
                        {"name": "cluster2", "size": 10},
                ],
                },
            ]
        }

#        n = len(Z)+1
#        d3_dict = _do_linkage_to_d3(n, len(Z)+n-1, Z)
        d3_dict = linkage_to_d3(Z)
        self.assertDictEqual(expected, d3_dict)

def nanmin(a, axis=None):
    """Find the minimium over the given axis, ignoring NaNs.
    """
    y = array(a,subok=True)
    if not issubclass(y.dtype.type, _nx.integer):
        y[isnan(a)] = _nx.inf
    return y.min(axis)

def helixHeteroscedasticDiags():
  #Parameterise a helix (no noise)
  fig5 = plt.figure()
  t = arange(-1,1,0.0005)
  x = map(lambda x: x + gauss(0,0.001 + 0.001*sin(2*pi*x)**2), (1 - t*t)*sin(4*pi*t))
  y = map(lambda x: x + gauss(0,0.001 + 0.001*sin(2*pi*x)**2), (1 - t*t)*cos(4*pi*t))
  z = map(lambda x: x + gauss(0,0.001 + 0.001*sin(2*pi*x)**2), t)
  line = array(zip(x,y,z))
  lpc = LPCImpl(h = 0.1, t0 = 0.1, mult = 1, it = 500, scaled = False, cross = False)
  lpc_curve = lpc.lpc(X=line)
  ax = Axes3D(fig5)
  ax.set_title('helixHeteroscedastic')
  curve = lpc_curve[0]['save_xd']
  ax.scatter(x,y,z, c = 'red')
  ax.plot(curve[:,0],curve[:,1],curve[:,2])
  saveToPdf(fig5, '/tmp/helixHeteroscedastic.pdf')
  residuals_calc = LPCResiduals(line, tube_radius = 0.2, k = 20)
  residual_diags = residuals_calc.getPathResidualDiags(lpc_curve[0])
  fig6 = plt.figure()
  #plt.plot(lpc_curve[0]['lamb'][1:], residual_diags['line_seg_num_NN'], drawstyle = 'step', linestyle = '--')
  plt.plot(lpc_curve[0]['lamb'][1:], residual_diags['line_seg_mean_NN'])
  plt.plot(lpc_curve[0]['lamb'][1:], residual_diags['line_seg_std_NN'])
  saveToPdf(fig6, '/tmp/helixHeteroscedasticPathResiduals.pdf')
  coverage_graph = residuals_calc.getCoverageGraph(lpc_curve[0], arange(0.01, .052, 0.01))
  fig7 = plt.figure()
  plt.plot(coverage_graph[0],coverage_graph[1])
  saveToPdf(fig7, '/tmp/helixHeteroscedasticCoverage.pdf')
  residual_graph = residuals_calc.getGlobalResiduals(lpc_curve[0])
  fig8 = plt.figure()
  plt.plot(residual_graph[0], residual_graph[1])
  saveToPdf(fig8, '/tmp/helixHeteroscedasticResiduals.pdf')
  fig9 = plt.figure()
  plt.plot(range(len(lpc_curve[0]['lamb'])), lpc_curve[0]['lamb'])
  saveToPdf(fig9, '/tmp/helixHeteroscedasticPathLength.pdf')

def twoDisjointLinesWithMSClustering():
 
  t = arange(-1,1,0.002)
  x = map(lambda x: x + gauss(0,0.02)*(1-x*x), t)
  y = map(lambda x: x + gauss(0,0.02)*(1-x*x), t)
  z = map(lambda x: x + gauss(0,0.02)*(1-x*x), t)
  line1 = array(zip(x,y,z))
  line = vstack((line1, line1 + 3))
  lpc = LPCImpl(start_points_generator = lpcMeanShift(ms_h = 1), h = 0.05, mult = None, it = 200, cross = False, scaled = False, convergence_at = 0.001)
  lpc_curve = lpc.lpc(X=line)
  #Plot results
  fig = plt.figure()
  ax = Axes3D(fig)
  labels = lpc._startPointsGenerator._meanShift.labels_
  labels_unique = unique(labels)
  cluster_centers = lpc._startPointsGenerator._meanShift.cluster_centers_
  n_clusters = len(labels_unique)
  colors = cycle('bgrcmyk')
  for k, col in zip(range(n_clusters), colors):
    cluster_members = labels == k
    cluster_center = cluster_centers[k]
    ax.scatter(line[cluster_members, 0], line[cluster_members, 1], line[cluster_members, 2], c = col, alpha = 0.1)
    ax.scatter([cluster_center[0]], [cluster_center[1]], [cluster_center[2]], c = 'b', marker= '^')
    curve = lpc_curve[k]['save_xd']
    ax.plot(curve[:,0],curve[:,1],curve[:,2], c = col, linewidth = 3)
  plt.show()

def helixHeteroscedasticCrossingDemo():
  #Parameterise a helix (no noise)
  fig5 = plt.figure()
  t = arange(-1,1,0.001)
  x = map(lambda x: x + gauss(0,0.01 + 0.05*sin(8*pi*x)), (1 - t*t)*sin(4*pi*t))
  y = map(lambda x: x + gauss(0,0.01 + 0.05*sin(8*pi*x)), (1 - t*t)*cos(4*pi*t))
  z = map(lambda x: x + gauss(0,0.01 + 0.05*sin(8*pi*x)), t)
  line = array(zip(x,y,z))
  lpc = LPCImpl(h = 0.15, t0 = 0.1, mult = 2, it = 500, scaled = False)
  lpc_curve = lpc.lpc(line)
  ax = Axes3D(fig5)
  ax.set_title('helixHeteroscedasticWithCrossing')
  curve = lpc_curve[0]['save_xd']
  ax.scatter(x,y,z, c = 'red')
  ax.plot(curve[:,0],curve[:,1],curve[:,2])
  saveToPdf(fig5, '/tmp/helixHeteroscedasticWithCrossing.pdf')
  lpc.set_in_dict('cross', False, '_lpcParameters')
  fig6 = plt.figure()
  lpc_curve = lpc.lpc(X=line)
  ax = Axes3D(fig6)
  ax.set_title('helixHeteroscedasticWithoutCrossing')
  curve = lpc_curve[0]['save_xd']
  ax.scatter(x,y,z, c = 'red')
  ax.plot(curve[:,0],curve[:,1],curve[:,2])
  saveToPdf(fig6, '/tmp/helixHeteroscedasticWithoutCrossing.pdf')

def findCoords(gs, candidates=None):
    if candidates == None:
        candidates=[]
        # List all the possible z-level (heights)
        zRange = list(takewhile(lambda x : x < gs.boardSize[2], \
                 sort(unique(flatten(gs.heightMap())))))
        if zRange==[]:
            print "Board is full, cannot find legal coordinates !"
            return None
    else:
        zRange = sort(unique(map(third,candidates)))
    # Do we have a choice on the z-level ?
    if len(zRange)==1:
        z = zRange[0]
    else:
        print "\n",gs.boardToASCII(markedCubes=candidates)
        # Discard the z height max
        if zRange[-1]==gs.boardSize[2]:
            zRange = zRange[:-1]
        z = -1+input("Which z-level ? (%d-%d)\n> " \
                     % (zRange[0]+1,zRange[-1]+1))
    candidates = filter(lambda c: c[2]==z, candidates)
    if len(candidates)>1:
        # Display the z-level with xy coordinates as letter-number
        print '    '+''.join(chr(97+x) for x in xrange(gs.boardSize[0]))
        print '   +'+'-'*gs.boardSize[0]
        lines = gs.boardToASCII(zRange=[z],markedCubes=candidates)\
                .split('\n')
        for y in xrange(gs.boardSize[1]):
            print '%s |%s' % (str(y+1).zfill(2),lines[y])
        print "\n"
        xy = raw_input("Which xy coordinates ?\n> ")
        return array([ord(xy[0])-97,int(xy[1:])-1,z])
    else:
        return candidates[0]

def nanargmax(a, axis=None):
    """Find the maximum over the given axis ignoring NaNs.
    """
    y = array(a,subok=True)
    if not issubclass(y.dtype.type, _nx.integer):
        y[isnan(a)] = -_nx.inf
    return y.argmax(axis)

 def test_multi_itr_array(self):
     
     c = array( [ [[  0,  1,  2],
                   [ 10, 12, 13]],
                 
                 [[100,101,102],
                  [110,112,113]] ] )

 def testNoisyLine2(self):
   x = map(lambda x: x + gauss(0,0.005), arange(-1,1,0.005))
   y = map(lambda x: x + gauss(0,0.005), arange(-1,1,0.005))
   z = map(lambda x: x + gauss(0,0.005), arange(-1,1,0.005))
   line = array(zip(x,y,z))
   lpc = LPCImpl(h = 0.2, convergence_at = 0.001, mult = 2)
   lpc_curve = lpc.lpc(X = line) 

 def testNoisyLine1(self):
   x = map(lambda x: x + gauss(0,0.002), arange(-1,1,0.001))
   y = map(lambda x: x + gauss(0,0.002), arange(-1,1,0.001))
   z = map(lambda x: x + gauss(0,0.02), arange(-1,1,0.001))
   line = array(zip(x,y,z))
   lpc = LPCImpl(h = 0.2, mult = 2)
   lpc_curve = lpc.lpc(X = line)

 def getCoverageGraph(self, curves, tau_range):
   '''Return a 2*len(tau_range) array of the proportion of self._X points within tau (in tau_range, an array)
   of the curve segments, where 'curves' is a either a list of curve dictionaries as returned by LPCImpl.lpc, or an element thereof
   This should give graphs similar to the output of BOakley
   '''
   coverage = [1.0*len(self.calculateCoverageIndices(curves,tau))/len(self._X) for tau in tau_range]
   return array([tau_range,coverage])

 def calculatePurity(self, curves, data_range, voxel_to_pdg_dictionary):
   '''NB - self._residuals_runner should have had calculateResiduals method called with calc_residuals = True beofre calling this method
   '''
   hit_tuples = voxel_to_pdg_dictionary.keys()
   if data_range is None:
     data_range = 1.0
   #rescales the truth data if necessary
   hits = array([[h[0], h[1], h[2]] for h in hit_tuples]) / data_range
   self._residuals_runner.setDataPoints(hits)
   self._residuals_runner.setLpcCurves(curves)
   self._residuals_runner.calculateResiduals(True, False, False)
   residuals = self._residuals_runner.getResiduals()
   tau_range = self._residuals_runner.getTauRange()
   purity = {}
   for tau in tau_range:
     pdg_code_frequencies = []
     for i in range(len(curves)):
       d = defaultdict(int)
       hit_labels = [voxel_to_pdg_dictionary[hit_tuples[i]] for i in residuals['curve_residuals'][i]['coverage_indices'][tau]]
       flattened_hit_labels = [pdg_code for pdg_code_list in hit_labels for pdg_code in pdg_code_list]
       for pdg_code in flattened_hit_labels:
         d[pdg_code] += 1
       pdg_code_frequencies.append(d)
     purity[tau] = pdg_code_frequencies
   return purity