算法学习笔记(5):Dijkstra

在刷Hackerrank的时候碰到了这道题,其本身倒是十分普通的一道题,照本宣科的把Dijkstra实现一遍就可以做出来了。在讨论区和别人讨论发现的时间复杂度也是会被Accept的。

如果使用数组储存所有unvisited点的距离,在找当前最近点的时候,会需要遍历整个数组,而这是的时间复杂度。如果使用一个最小优先级队列,那么只需要将队首的元素取出,然后min-heapify一下,保持最小堆的性质,而这是的时间复杂度。能够显著减少时间消耗。

但是如果使用最小堆的话,降低堆中元素的优先级的操作就很麻烦了。首先如果想要找到堆中的某个元素就需要的时间复杂度,而修改了这个元素的优先级之后,想要维持最小堆的性质又要的时间复杂度。所以整体而言就变得和没有使用最小堆一样了。

为了解决这个问题有一个弥补的方法。就是在一个散列表中存储优先级队列中对应点的元素的handle。

如果想要降低某个点的优先级,那么可以通过查询散列表,获得最小堆中对应元素的handle,然后直接标志该点已经被删除(但是不实际从堆中删除),然后在储存堆的数组的最后新添加一个储存新优先级和该点的元素。标志删除的操作因为没有破坏最小堆的性质,只是 的时间复杂度,而最小堆中插入一个元素是 的时间复杂度。最后就实现了 的降低优先级的操作。

如果使用python实现的话,python的官方文档中简单介绍了使用heapq实现最小优先级队列的方法。

不过残念的是,add_task()方法的第一行就出现了个遍历dict,而这是的,所以推荐额外再使用一个set来维护entry_finder.keys(),这样查找和删除操作就都是 了。

另外一个有意思的事情是,heappush()原地的,所以pq会在原地修改。而pq中保存的是[distance, vertex]的list,而list是mutable的。所以在remove_vertex()把entry[-1]也就是vertex标志城REMOVED了之后,pq中的对应entry会自动修改。

最后挂上自己的代码。

# Enter your code here. Read input from STDIN. Print output to STDOUT
from heapq import *
 
def add_vertex(vertex, distance=0):
    if vertex in entry_set:
        remove_vertex(vertex)
    entry = [distance, vertex]
    entry_finder[vertex] = entry
    entry_set.add(vertex)
    heappush(pq, entry)
 
def remove_vertex(vertex):
    entry = entry_finder.pop(vertex)
    entry[-1] = REMOVED
    entry_set.remove(vertex)
 
def pop_vertex():
    while pq:
        priority, vertex = heappop(pq)
        if vertex is not REMOVED:
            del entry_finder[vertex]
            return priority, vertex
 
T = int(raw_input())
for T_ in range(T):
    pq = []
    entry_finder = {}
    REMOVED = "REMOVED"
 
    N, M = [int(i) for i in raw_input().split(" ")]
    connection = [[] for conn_ in range(N)]
    for M_ in range(M):
        x, y, r = [int(i)-1 for i in raw_input().split(" ")]
        connection[x].append([y, r+1])
        connection[y].append([x, r+1])
    S = int(raw_input())-1
     
    entry_set = set()
    unvisited = set()
    distance = {}
 
    for node in range(M):
        unvisited.add(node)
        if node!=S:
            distance[node] = float("Inf")
            add_vertex(vertex = node, distance = float("Inf"))
        else:
            distance[node] = 0
            add_vertex(vertex = node, distance = 0)
 
    while len(unvisited)!=0:
        current_distance, current_node = pop_vertex()
        if current_distance==float("Inf"):
            break
        for node, edge in connection[current_node]:
            if node in unvisited:
                new_distance = current_distance + edge
                if new_distance < distance[node]:
                    distance[node] = new_distance
                    add_vertex(node, new_distance)
                distance[node] = new_distance if new_distance < distance[node] else distance[node]
        unvisited.remove(current_node)
 
    for node in range(N):
        if node!=S:
            if distance[node]==float("Inf"):
                print -1,
            else:
                print distance[node],
    print

GMM的EM算法实现

GMM的全名是Gaussian Mixture Model,即高斯混合模型,这个名字挺直观。

如果给了一组数据点,假设这些数据点符合iid,那么很容就能根据数据拟合出一个最佳的高斯分布。

但是如果给了一堆数据点,而他们却不符合iid,只能做到互相独立,并不能做到同分布。相反的,数据点们有可能来自不同的高斯分布,这也就是『高斯混合』模型名字的意思。这样如何能拟合出最佳的高斯分布们呢?

EM算法提供了一种方法,这种方法的思路有些像牛拉法。先是确定一个初值,然后看看当前值距离目标还有多远,然后逐渐的缩小差距(不断修正membership weight)。

具体原理见这里

这篇小短文大概记录下这次作业的内容。作业题目是生成了2个二维高斯分布混合的数据点,共2000个,然后用这两千个点拟合出2个二维高斯分布的参数。

clear;
%产生2个二维正态数据
MU1    = [1 2];
SIGMA1 = [1 0; 0 0.5];
MU2    = [-1 -1];
SIGMA2 = [1 0; 0 1];
X      = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];
scatter(X(:,1),X(:,2),10,'.');

%initial values
h_mu1 = X(randi(2000),:);
h_sigma1 = cov(X);
h_mu2 = X(randi(2000),:);
h_sigma2 = cov(X);
h_alpha1 = 0.5;
h_alpha2 = 0.5;

%simple iteration
for iter=1:100
    %E step
    w1 = mvnpdf(X,h_mu1,h_sigma1)*h_alpha1./...
    (mvnpdf(X,h_mu1,h_sigma1)*h_alpha1 + mvnpdf(X,h_mu2,h_sigma2)*h_alpha2);
    w2 = mvnpdf(X,h_mu2,h_sigma2)*h_alpha2./...
    (mvnpdf(X,h_mu1,h_sigma1)*h_alpha1 + mvnpdf(X,h_mu2,h_sigma2)*h_alpha2);

    %M step
    h_alpha1 = sum(w1)/2000;
    h_alpha2 = sum(w2)/2000;
    
    h_mu1 = sum([w1.*X(:,1) w1.*X(:,2)])/sum(w1);
    h_mu2 = sum([w2.*X(:,2) w2.*X(:,2)])/sum(w2);

    M1 = X-ones(2000,1)*h_mu1;
    h_sigma1 = ([sum(w1.*M1(:,1).^2) 0;0 sum(w1.*M1(:,2).^2)]+...
        sum(w1.*(M1(:,1).*M1(:,2)))*[0 1;1 0])/sum(w1);
    M2 = X-ones(2000,1)*h_mu2;
    h_sigma2 = ([sum(w2.*M2(:,1).^2) 0;0 sum(w2.*M2(:,2).^2)]+...
        sum(w2.*(M2(:,1).*M2(:,2)))*[0 1;1 0])/sum(w2);
end

figure;
[x,y] = meshgrid(linspace(-5,5,75));
p1 = mvnpdf([x(:),y(:)],h_mu1,h_sigma1);
p2 = mvnpdf([x(:),y(:)],h_mu2,h_sigma2);
surf(x,y,reshape(p1,75,75));hold on
surf(x,y,reshape(p2,75,75));

options = statset('Display','final');
obj = gmdistribution.fit(X,2,'Options',options);
figure,h = ezmesh(@(x,y)pdf(obj,[x,y]),[-8 6], [-8 6]);

绝大多数时候都没有什么问题,不过由于初始点是从X里随便选的一个数据点,极偶然会跑偏。

以下分别是数据点的分布、自行拟合出的结果、利用matlab自带拟合函数拟合出的结果。

originalgmm

根据给定的均值、协方差矩阵画出的GMM的图

xscatter

数据点的散点图,保存之后可能有点糊

fitbymyself

自己写的脚本拟合出的GMM的图

fitbymatlabtoolbox

matlab提供的拟合函数拟合出的图

迭代100次之后,高斯分布1的均值分别偏差0.2776%、4.4462%,协方差分别偏差了2.4458%、5.0146%;高斯分布2的均值偏差了2.6981%、1.2911%,协方差偏差了3.8498%、1.5517%。精度凑凑活活把...

对比之下,matlab自带的拟合函数拟合出来的高斯分布1的均值偏差2.2415%、2.2415%,协方差偏差为4.02%、1.32%;高斯分布2的均值偏差了0.4165%、0.4165%,协方差偏差为2.39%、4.62%。精度也没高到哪里去...不过均值偏差一样让我很是在意...

然而看图1、图3和图4的话,明显还是自己写的代码拟合效果更好一些啊,有些不科学。

算法学习笔记(4):计数排序

在宿舍呆着无聊接着在hackerrank上刷水题。这道题由于排序的依据是一定范围内的整数,所以用计数排序效率应当不错。

说到计数排序,实际上思路也很常见。最常见的一个例子就是比赛或者是考试的时候要进行排名。

比如某班有3名同学考了100分,2名同学考了99分,5名同学考了98分。那么很明显,比99分高的同学一共有3名,所以考99分的同学只能是第4名了。

计数排序大概就是做了个数数每个元素前面有多少元素这个事情。先统计最高分多少分,然后依次统计各个分数的学生个数,然后以此叠加就可以统计出前面有多少个学生了。

而且由于不同分数段的学生是按照原顺序放置到新数组中的,所以在原数组中的顺序并不会被打乱,即这个算法是稳定的。这也是为什么这道题可以用计数排序解,因为并不会打乱文本的顺序。

顺便一提,这道题中的最大值max()的复杂度应该是,所以总运行时间应该是。因为计数排序并不涉及到不同元素之间的比较,即不是比较排序,所以的下界并不适用于计数排序。

最后附上这道题的一个无脑解。

# Enter your code here. Read input from STDIN. Print output to STDOUT
n = int(raw_input())
ar = []
idx = []
out = [None]*n
for _ in range(n):
    ar.append(raw_input().split())
    ar[-1][0] = int(ar[-1][0])
    idx.append(ar[-1][0])
    if _<n/2:
        ar[-1][1] = "-"
countlist = [0]*(max(idx)+1)
for item in ar:
    countlist[item[0]]+=1
for i in range(len(countlist)-1):
    countlist[i+1]+=countlist[i]
for item in reversed(ar):
    out[countlist[item[0]]-1]=item
    countlist[item[0]]-=1
for item in out:
    print item[1],

算法学习笔记(3):插入排序、归并排序与逆序数

前些天在挨个做hackerrank上的题目的时候,碰到了这么一道题

题目的完成率并不高,只有百分之三十六点多,我估计了下难度,自己应该是做不出来的,最后果然没做出来。不过在Discussion下面还是看到了些有用的帮助,有人说用mergesort能做出来。于是又仔细的想了想,还是花了些力气做出来了...

题目的内容是计算插入排序一个序列的移动元素的总次数。实际上是计算的序列的逆序数

逆序数的定义上面那个wiki的链接里说的挺清楚的了,数的就是反序排列的数字对数。

比如[1,3,2]的逆序数是1,在插入排序的时候,2会向左移动一次,整个序列变成了[1,2,3]。

[2,3,1]的逆序数是2,在插入排序的时候,1会向左移动两次,整个序列变成[1,2,3]。

hackerrank在排序这一章前面也有一道题是计算插入排序移动次数的,当时按照题目的提示,使用的直接在插入排序的时候,移动一次,加一次,最后依次数出来移动的次数。当然这个的时间复杂度就是,在序列比较长的情况下很容易超时。这种比较『暴力』的方法,显然不适合本题。

mergesort是一种比较简单的方法。可以很容易的想象出,在merge过程中,如果『小序列』非空,每当『大序列』的最小元素被选出来,就意味着逆序数需要自增目前『小序列』中元素个数那么多,因为这个时候『小序列』中的所有元素都比位于他们之后的『大序列』中挑出的最小元素大。

而且mergesort的时间复杂度是,已经是比较排序的最优了,所以也理应在本题中不会超时。

知道这个关键点,代码实现就很容易了,就是增添个自增的语句而已,代码如下所示。

def merge(a,b):
    global answer
    a_flag = 0
    b_flag = 0
    result = []
    while a_flag<len(a) and b_flag<len(b):
        if a[a_flag]>b[b_flag]:
            result.append(b[b_flag])
            answer+=(len(a)-a_flag)
            b_flag+=1
        else:
            result.append(a[a_flag])
            a_flag+=1
    if a_flag==len(a):
        result.extend(b[b_flag:])
    else:
        result.extend(a[a_flag:])
    return result

def mergesort(l):
    if len(l)>1:
        mid = len(l)/2
        a = l[:mid]
        b = l[mid:]
        a = mergesort(a)
        b = mergesort(b)
        return merge(a,b)
    else:
        return l

n = input()
for iterate in range( n ):
    x = input()
    a = [ int( i ) for i in raw_input().strip().split() ]
    answer = 0
    mergesort(a)
    print answer

算法学习笔记(2):快排

前言:最近两天又改变了下学习方法。直接看《算法导论》,然后看完立刻实现感觉太直接,自己不需要动脑子,引导性也不是很强。于是去hackerrank上一点点学。昨天把warmup和implement完成了,熟悉了下大概环境。今天把插排和快排看完了。

网上之前一直流传的一个图书馆大妈给书排序的段子形象的解释了快排的思路。

图书馆看见两个志愿者需要把还回来的一堆书按顺序入架, 管理员大妈给他们教的时候说:『你先在这堆书里拉出一本来,把比它号小的扔到一边,比它大的扔到另一边,然后剩下两堆继续这样整,这样排的快!』

从分治法的思维来看,快排和堆排的路数差不多。堆排中是假设左子树和右子树是已经排好了的树,而快排中则是假设左侧的序列和右侧的序列分别是排好了的序列。其实就是每次排序之前都能保证小的在左面,大的在右面。思路很简朴。

对比之下,由于堆排中使用了二叉树的结构,所以在排序这个事情上徒增了复杂度,所以导致了其比快排效率低。

过来人也都对这种『凡事不能搞得太复杂』的思想也有过阐述。

Everything should be made as simple as possible, but no simpler.

——爱因斯坦

万物之始,大道至简,衍化至繁。

——老子《道德经》

说回图书馆大妈的段子。毕竟图书馆大妈不是程序员,只能想到最简朴的快排的方法。事实上,在把书分成三摞(即『比Pivot小的一摞』、『Pivot一摞』、『比Pivot大的一摞』)的不停迭代的过程中,实际上是个把最初的一整摞书摊开的过程。形象的说就是原来一摞书撑死只占A4纸那么大的面积,然而摊开之后会占许多倍原来的面积。即段子中的排序实际上是会多占空间的。

对应的也有一种『原地快排』的算法,并不需要把书按着大小摊开,只需要把书从那一摞里『抽出来』再『塞回去』(对应着其实就是数组元素之间的互相调换)即可以完成排序。《算法导论》上杭也介绍的是这种快排方法。

这种原地快排的代码也不是很复杂,书上写的也很清晰。我觉得一个挺巧妙的地方是《算法导论》中使用两个下标来圈定了『大于Pivot的右侧子数组』。这样处理的时候很直观。

别的关于快排也并没有什么了,毕竟学的时候感觉一切都是那么顺理成章。又想起来小波课上老师说的,『一看都不对的结论肯定有问题』。

算法学习笔记(1)——堆排序

前记:最近回家除了干干所里剩下的活,手头剩了不少时间。于是趁着这个寒假比较长,准备把算法导论看完了。不过隔了很久突然看起,有些生疏,于是准备再实现一遍。选择语言秉承了「人生苦短,我用python」的原则,果断的选择了python。

虽然把代码提交到了github上,不过感觉不能看着书把伪代码翻译一遍就算完了,写完了之后好好感悟一下,记录一下自己的感想也是很重要的,于是准备在博客上大概记录下自己的学习笔记和感想。

所以重点不在于名词定义和原理解释,毕竟wiki上也好,书上也好,都是一翻到处都有。我想记下来的更多的还是自己个人在学习过程中觉得这些算法构思巧妙的地方。

第一天的复习内容是堆排序,我认为巧妙的地方有以下几点。

1.首先这里的堆使用的数据结构是二叉树,在顺序存储二叉树的时候,父节点、左子节点和右子节点的下标之间具有天然的简单倍数、邻近关系。这直接决定了在访问数据时,不需要很大的代价,客观上显著降低了复杂度。

访问各个节点的代码如下所示,为了简便并没有考虑下标越界的异常处理,这也是为了突出原理。用书上的话说就是:

数据抽象、模块化和错误处理等问题往往都被忽略掉了,以便更简练的表达算法的核心内容。

def get_parent(self, index):
	return self.heaplist[index/2+1]

def get_leftchild(self, index):
	return self.heaplist[2*index+1]

def get_rightchild(self, index):
	return self.heaplist[2*index+2]

2.其次,堆排序将一个复杂的工作分成了许多层次,即所谓的分治法(divide-and-conquer)。

逆向思考的话,就是将一个排序问题分解成了「提取最大元素」、「构建最大堆」、「使用两个最大堆一个节点组成一个最大堆」三个层次。

而如果正向思考的话,由于最大堆的最大的节点就是根节点,所以实际问题就变成了「如何将最小的元素都放在叶节点上」、「如何将根节点和某个叶节点交换之后,再次构建最大堆」两步。

虽然实际上使用的方法都是相同的,即从叶节点开始(由于是顺序存储结构,只需要从数组的中间开始倒序进行)构建最大堆,最后自然整个树会变成最大堆;而后将根节点拿出来(《算法导论》中是将根节点放到了数组末尾,由于放到数组末尾之后并不是二叉树的组成部分,其下标与其他节点的下标没有关系,实际就是放到了另一个数组里)。最后自然形成了一定的大小顺序。

从这种角度来看,其实这种排序方法与插排的思想差不多。就是依次为各个元素找到合适的位置,而由于在元素比较多的时候,第一次构建最大堆的时候,所有的元素形成了一定的有序性。所以一般堆排序比插排的效率要高,不过也不是很高,因为受到了算法设计思维(单个、依次操作)的限制。

实现的代码如下所示。

def max_heapify(self, i=0):
	l = 2*i+1
	r = 2*i+2
	if l <= self.length-1 and self.nodelist[l] > self.nodelist[i]:
		largest = l
	else:
		largest = i
	if r <= self.length-1 and self.nodelist[r] > self.nodelist[largest]:
		largest = r
	if largest != i:
		self.nodelist[i], self.nodelist[largest] = self.nodelist[largest], self.nodelist[i]
		self.max_heapify(largest)

def build_max_heap(self):
	for i in reversed(range(0, self.length-2-1)):
		self.max_heapify(i)

def heapsort(self):
	self.build_max_heap()
	for index in reversed(range(1, self.length)):
		self.sortlist.append(self.nodelist[0])
		self.nodelist[index], self.nodelist[0] = self.nodelist[0], self.nodelist[index]
		del self.nodelist[index]
		self.length = len(self.nodelist)
		self.max_heapify()
	self.nodelist = self.heaplist
	self.length = len(self.nodelist)

最终所有的代码在heapsort.py