LeetCode-105-前序和中序构建二叉树

给定一个二叉树前序遍历数组A 124367,中序遍历得到数组B, 421637。我们的任务是从这两个结果中构造出原来二叉树。

根据前序遍历的顺序,可以知道数组是先存放根,然后是左子树,最后是右子树。

根据中序遍历的顺序,可以知道数组是先放左子树,然后是根,最后是右子树。

根据定义,前序遍历数组第一个位置是根节点,但左子树和右子树大小未知。而中序遍历一旦确定根的位置,就可以确定左右子树的大小。这两者刚好是互补关系。

因此具体过程是

  • 用前序数组的第一个位置确定中序数组的根位置
  • 将中序数组分成两个部分,得到左右子树的大小。
  • 用左右子树的大小确定前序数组的左右子树位置
  • 对于新生成的左右子树,用上面方法确定新的根位置和左右子树大小。

对于本问题,就是1将421637分成42和637,这意味着左子数大小为2,右子数大小为3.于是124367被分成1,24,367.

24和42是新的左子树问题,367和637是新的右子树问题。

理清思路之后,我们可以按照递归模版四步走起来

第零步:思考函数名和传递的参数和全局变量。根据我们的分析,我们最重要的就是子问题对应数组的区间,也就是我们需要记录前序数组的左右子树开始和结束(pb, pe), 中序数组的开始和结束(ib, ie), 函数名的话,就写dc吧

第一步:终止条件。这里终止条件是前序数组用完了(当然中序数组肯定也用完了)

1
if (pb > pe) return NULL;

第二步: 处理逻辑。每次用前序节点的第一个位置构建根节点,然后拆分中序节点,获取左右子树的大小,

1
2
3
4
TreeNode *root = new TreeNode(preorder[pb]);
int rootIndex = mp[root->val]; //获取中序的位置
int leftSize = rootIndex - ib; //左子树大小
int rightSize = ie - rootIndex; //右

第三步,下潜一层, 新的一层要带新的参数过去。

1
2
root->left = dc(preorder,pb+1,pb+leftSize, inorder,ib, rootIndex-1);
root->right = dc(preorder,pb+leftSize+1, pe, inorder, rootIndex+1, ie);

第四步:状态重置。此处不需要

最终返回root即可。

对应的完成代码见 https://github.com/xuzhougeng/algorithm006-class01/blob/master/Week_02/G20200343030533/LeetCode_105_533.cpp

LeetCode-042-接雨水

题目地址: https://leetcode-cn.com/problems/trapping-rain-water

给定 n 个非负整数表示每个宽度为 1 的柱子的高度图,计算按此排列的柱子,下雨之后能接多少雨水。

接雨水

上面是由数组 [0,1,0,2,1,0,1,3,2,1,2,1] 表示的高度图,在这种情况下,可以接 6 个单位的雨水(蓝色部分表示雨水)。 感谢 Marcos 贡献此图。

第一种方法是暴力求解,对于给定的一个柱子,我们分别找到它的左右边界,然后其中的最小值减除当前柱子的高度,就是它能接的水

i方法1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
int min(int a, int b){
return a > b ? b : a;
}

int max(int a, int b){
return a < b ? b : a;
}

int trap(int* height, int heightSize){

int left_max;
int right_max;
int ans=0;
//从1到n-1即可,因为第一个和最后一个不能有水
for (int i = 1; i < heightSize - 1; i++){
left_max = height[i];
right_max = height[i];
for (int j = 0; j < i; j++){
left_max = max(left_max, height[j]);
}
for (int j = heightSize - 1; j > i; j--){
right_max = max(right_max, height[j]);
}
//计算接水面积
ans += min(left_max, right_max) - height[i];

}
return ans;

}

这一种方法是两层循环,时间复杂度是O(n^2). 其中第二层的循环是实时计算当前柱子的左右边界。如果已知每个柱子的左右边界,我们就可以省掉嵌套循环,可以用 ans += min(left_max[i], right_max[i]) - height[i];直接计算面积。问题就在于如何提前计算好left_max[i]right_max[i]。我们只需要分别为左右遍历数组,在遍历过程中,比较当前值和前一个值,用其中较大值更新当前值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
class Solution {
public:
int trap(vector<int>& height)
{
if(height.size() == 0) return 0;
int ans = 0;
int size = height.size();
vector<int> left_max(size), right_max(size);
left_max[0] = height[0];
for (int i = 1; i < size; i++) {
left_max[i] = max(height[i], left_max[i - 1]);
}
right_max[size - 1] = height[size - 1];
for (int i = size - 2; i >= 0; i--) {
right_max[i] = max(height[i], right_max[i + 1]);
}
for (int i = 1; i < size - 1; i++) {
ans += min(left_max[i], right_max[i]) - height[i];
}
return ans;
}

};

这里用了2个额外的数组,避免了嵌套循环,用空间换时间。

方法3,前面是按列计算每一列的接水量,我们其实还可以按行计算每一行的接水量。

我们依次遍历柱子,如果当前柱子比之前的小,那就说明它有可能和之前的柱子组成一个凹型空间,也就是说它的左边界就是它的前一个柱子,但是右边界我们暂时不清楚。如果当前柱子比上一个元素大,就说明前一个元素的右边界确定了。下面是两个特殊的例子,下图左直到第五个柱子才确定了第四根柱子的左右边界,而下图右则到头都没有找到有边界

方法3
这种依次加入元素,然后从后往前取出元素的方式可以让我们联想到这种先进后出数据结构。具体算法逻辑如下

  • 读取最新元素,
  • 和栈顶元素对应的高度比较,
  • 如果低于栈顶元素对应高度,则直接入栈
  • 如果高于栈顶元素对应高度,则开始计算面积

面积的计算公式为,凹型区间两边柱子的较低者减去凹形区间最低元素高度,乘以两边柱子的距离。

代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
class Solution {
public:
int trap(vector<int>& height) {
int ans = 0;
stack<int> st;
int current = 0;
while( current < height.size()){
while( !st.empty() && height[current] > height[st.top()]){
//获取中间位置
int mid = st.top();
st.pop();
//获取左侧位置
if ( st.empty()) {
break;
}
int left = st.top();
//左右间距
int distance = current - left - 1;
//面积=高度差 x 距离
int bound_height = min(height[left], height[current] ) - height[mid];
ans += distance * bound_height;
}
st.push(current++);
}
return ans;
}
};

LeetCode-084-最大面积

题目地址: https://leetcode-cn.com/problems/largest-rectangle-in-histogram/

给定 n 个非负整数,用来表示柱状图中各个柱子的高度。每个柱子彼此相邻,且宽度为 1 。求在该柱状图中,能够勾勒出来的矩形的最大面积。

example

第一种方法,我们先确定左右两个边界,然后找边界中的最小值。比方说,我们左边界确定为(0,2),右边界确定为(4,2), 然后遍历中间元素,发现最小值是(1,1),那么面积就是(4-0) x 2 = 8

方法1

代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
int min(int a, int b){
return a > b ? b : a;
}
int max(int a, int b){
return a > b ? a: b;
}

int largestRectangleArea(int* heights, int heightsSize){

int max_area = 0;

for (int i = 0; i < heightsSize; i++){
int min_height = heights[i];
for (int j = i; j < heightsSize; j++){
//找到i-j中的最小高度
for ( int k = i; k <= j; k++){
min_height = min(min_height, heights[k]);
}
//计算面积
max_area = max(max_area, (j-i+1) * min_height);
}
}
return max_area;

}

我写代码时出错的两个地方

  • k <= j而不是k<j, 要包括最后一个位置。
  • 面积的计算公式为(j-i+1) * min_height, 我漏了+1,结果(0,1)我算成了0,而实际面积是(0-0+1)x1=1

上面方法是固定左右边界然后找中间的最小值,用到了三层循环,计算效率比较低。我们可以想办法省掉第三层循环,比如说当我们计算完(0,2)-(2,5)内的最小值后,对于(0,2)-(3,6)的最小值,其实只要比较当前高度和之前的最小值,如果比前面的小,就更新最小值,否则就用之前的最小值。代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
int largestRectangleArea(int* heights, int heightsSize){
int min_height = 0;
int max_area = 0;

for ( int i = 0 ; i < heightsSize; i++){
min_height = heights[i];
for ( int j = i ; j < heightsSize; j++){
min_height = min(min_height, heights[j]);
max_area = max(max_area, (j-i+1) * min_height);
}
}
return max_area;

}

或者我们可以换个思路,先把中间值固定住,然后找左右边界。这里的边界指的是左右出现的第一个比他小的棒子,比如说(1,1)就是最两边。

方法2

代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
int largestRectangleArea(int* heights, int heightsSize){
int max_area = 0;

for ( int i = 0; i < heightsSize; i++){
int left_border = 0;
int right_border = heightsSize;

for (int j = i; j < heightsSize; j++){
if (heights[j] < heights[i]){
right_border = j;
break;
}

}
for (int j = i; j >= 0; j--){
if (heights[j] < heights[i]){
left_border = j+1;
break;
}
}
max_area = max(max_area, (right_border-left_border) * heights[i]);

}
return max_area;

}

依旧需要两层循环,其中第二层循环的目的是就是找到左右边界。那有没有方法不需要用到循环就能找到左右的边界呢?

思路和第二种方法类似,通过记录已经出现的最小值来避免多余计算,只不过这里用栈处理,具体步骤为

  • 栈初始化,入栈-1
  • 对于每一个新元素,都和栈顶元素比较
  • 如果比栈顶元素大
    • 该元素就入栈
  • 如果比栈顶元素小
    • 先取出栈顶元素,计算栈顶元素对应的面积
    • 重复上面步骤,直到比栈顶元素大
    • 入栈
  • 遍历结束后,清空栈

我们以最特殊的两个数据来举例,对于[0,1,2,3,4,5,6,7], 每一个元素都比之前的小,那么每个元素入栈的时候,我们都只能确定它的左边界,也就是它的上一个元素,而无法确定它的右边界,比如说(1,1)的左边界就是(0,0), 而右边的位置必须等到所有元素都入栈之后才能确定,最后算出来的面积是(8-1)x1=7.

方法3-1

对于[7,6,5,4,3,2,1,0], 我们先入栈(0,7),然后入栈(1,6), 此时对于高度为7的棒子而言,它已经到头了,面积只肯能是7. 一波操作之后,栈内部元素为(1,6), 此时来了(2,5), 那么6也到头了,它的面积就是(1-(-1) x 6 = 12.,其中-1栈初始后第一个元素。

最终代码如下(大部分代码是实现栈的操作)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
int max(int a, int b){
return a > b ? a : b;
}

typedef struct {
int *arr;
int count;
} Stack;

Stack *create(int k){
Stack *st = malloc(sizeof(Stack));
st->arr = malloc(sizeof(int) * (k+2));
st->count = 0;
return st;
}


void push(Stack *st, int val){
st->arr[st->count] = val;
st->count++;
return ;
}

int peek(Stack *st){

return st->arr[st->count-1];

}
int pop(Stack *st){
st->count--;
return st->arr[st->count];
}

void destroy(Stack *st){
free(st->arr);
free(st);
return;
}


int largestRectangleArea(int* heights, int heightsSize){

Stack *st = create(heightsSize);
push(st, -1);

int max_area = 0;
for (int i = 0; i < heightsSize; i++){

//如果不为空,且当前元素小于栈顶元素
while ( && heights[i] < heights[peek(st)]){

max_area = max(max_area, heights[pop(st)] *( i - peek(st)-1 )) ;

}
push(st, i);
}
while ( peek(st) != -1){
max_area = max(max_area, heights[pop(st)] *( heightsSize - peek(st) - 1) );
}
destroy(st);
return max_area;
}

我写代码出错的两个地方

  • peek(st) != -1用于判断栈是否为空,因为至少有一个-1
  • 计算面积的代码是heights[pop(st)] *( i - peek(st) - 1), 这里要多减去一个1,因为此时的i多偏移了1个单位,见下图
方法3-2

最终借用额外到数据结构,时间复杂度从原来的O(n^3)优化到最后到O(n), 也就是利用空间换取了时间。

0-1背包问题

对于一组不同重量、不可分割的物品,我们需要选择一些装入背包,在满足背包最大重量限制前提下,背包中物品总重量的最大值是多少?假设此时是5个物品,2,2,4,6,3,然后背包最大承载两是9.

假如我们使用回溯算法解决该问题, 代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
int maxW = 0; //最大重量
int n = 5; //物品数目
int w = 9; // 背包最大重量
int weight[] = {2,2,4,6,3};// 物品重量,2,2,4,6,3

void rucksack(int i, int cw)
{
if ( cw == w || i == n ){
if ( cw > maxW ) maxW= w;
return ;
}
rucksack(i + 1, cw);//不装第i个物品
if ( cw + weight[i] <= w){ //如果装的下
rucksack(i + 1, cw + weight[i]);
}
}

如果将代码执行过程产生状态画成树,我们可以发现对于不加入物品2的选择f(1,0)和加入物品2的选择f(1,2), 在下一个选择时,他们有一个相同的状态,f(2,2)。

递归树

虽然得到它的过程不同,但是从这往后,大家都是一样的。如同你走迷宫,走到一个地方,你发现墙上有一张纸写着,”此路不通”,你就可以不用白费力气去探索了。因此,如果一个问题解决时存在重复子问题,我们可以通过记忆化的方式,避免重复运算,提高计算效率。

从动态规划的角度,我们可以将整个求解过程分为n个阶段,每个阶段都需要决策是否需要将物品放到背包中。每个物品的决策后,背包中物品的重量就有会有种情况,也就是达到了不同的状态,也就是递归树中的不同节点。在每一个层中,我们只记录不同的状态(比如说上图的第2层的两个f(2,2)就可以合并成一种情况,当然第四层就更多了)。这样一来,我们就保证了每一层的状态数就不会超过w个(w是背包的承载重量)。这种合并操作就可以认为是一种记忆化。

我们先用一个二维数组来记录每层可以达到的不同状态

初始状态

考场重量为2的物品后,会出现两种状态,0和2。

状态1

在上一个状态的基础再考察一个重量为2的物品,会有三种状态,一直不选择,先不选再选2,两次都选2.

状态2

继续考察重量为4的情况时,会出现五种情况,其中重量为4可能来源是0+0+4,2+2+0,这种重复状态就被合并成一种状态,因此减少了计算量。

状态3

最终状态如下

最终状态

将上面的思考过程翻译成代码就是,

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
int rucksackDp(int *weight, int num, int w)
{
//先定义一个二维数组
int **status = (int **)malloc( sizeof(int *) * num);
//初始化数组
for (int i = 0; i < num; i++){
status[i] = (int *)calloc( w+1, sizeof(int));
}
//初始化第一行的值
status[0][0] = 1; //不放
if ( weight[0] < w){
status[0][weight[0]] = 1; //放
}
//从第二个开始考虑
for (int i = 1; i < num; i++){
//不放的情况,直接将上面的结果复制给当前
for (int j = 0; j <= w; j++){
if (status[i-1][j] == 1) status[i][j] = status[i-1][j];
}
//上面循环可以写成
//memcpy(status[i], status[i-1], sizeof(int) * num);
//放: 在上一个状态基础上, 将增加后的重量对应位置设置为1
for ( int j = 0; j <= w -weight[i]; j++){
if ( status[i-1][j] == 1) status[i][j+weight[i]] = 1;
}
}
//输出结果
for ( int i = w; i >= 0; --i){
if (status[n-1][i] == 1) return i;
}
return 0;
}

上面我们使用的是二维数组用于保存所有状态,但实际上这里我们只需要一维数组维护前一个状态就可以推导出当前结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
int rucksackDp2(int *weight, int num, int w)
{
//先定义数组
int *status = (int *)calloc( num, sizeof(int ) );

//初始化第一行的值
status[0]= 1; //不放
if ( weight[0] < w){
status[weight[0]] = 1; //放
}
//从第二个开始考虑
for (int i = 1; i < num; i++){
//不放: 保持状态不变
//放: 在上一个状态基础上, 将增加后的重量对应位置设置为1
for ( int j = 0; j <= w -weight[i]; j++){
if ( status[j] == 1) status[j+weight[i]] = 1;
}
}
//输出结果
for ( int i = w; i >= 0; --i){
if (status[i] == 1) return i;
}
return 0;
}

写动态规划代码的关键在于状态定义和状态转移方程。在0-1背包问题中,我们定义的状态是status[i]就是当前决策结束后到达的重量,而转移方程就是if ( status[j] == 1) status[j+weight[i]] = 1;

用C语言实现单链表操作

用C写一个链表

链表(Linked List)是一种非连续的线性数据结构,相对于数组,它允许数据在内存中非连续存储,但是不支持随机读取。

链表

链表由一个个节点(Node)组成,每个节点除了记录数据以外,还需要记录下一个节点的位置(如果是双向链表,还需要记录上一个节点的位置)

1
2
3
4
5
6
7
struct _Node;
typedef struct _Node Node;

struct _Node{
int data; //记录整型数据
Node *next;
};

对于第一个节点,我们有一个指针指向它的地址,对于最后一个节点,它需要指向NULL,表示链表结束了。

1
2
3
4
5
typedef struct _List {
Node *head; //记录头地址
Node *tail; //记录尾巴地址
int num;
} List;

有了链表的数据结构后,我们需要定义三个基本函数,用于创建链表,往链表中加入数据和删除链表

创建链表比较简单,就是为链表分配内存,并将其赋值给一个指针,然后返回

1
2
3
4
5
6
7
8
// 创建链表
List *CreateList()
{
List *list;
list = (List*)malloc( sizeof(List) );
list->num = 0;
return list;
}

加入数据时,我们需要先声明两个节点指针,第一个用于记录当前节点的位置,第二个是记录新节点的位置。如果链表中没有节点,也就是head指向为NULL,那么直接插入新节点即可。如果链表中已经有了节点,那么获取最后第一个节点的位置, 然后在它的后面加入节点,同时将tail指向新的节点。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
bool AddNode(List *list, int data)
{
Node *node;
Node *new_node;

new_node = (Node *)malloc( sizeof(Node) );
if ( new_node == NULL) return false;
new_node->data = data;
new_node->next = NULL;

//获取链表head
node = list->head ;
//如果head指向NULL, 则直接插入到下一个
if ( node == NULL){
list->head = new_node;
list->tail = new_node;
list->num = 1;
return true;
}
// 否则在尾部插入节点
node = list->tail ;
node->next = new_node;
list->tail = new_node;
list->num+=1;

return true;

}

删除列表分为两步,先删除节点内容,然后删除列表这个结构。如果节点存放的数据是其他结构,那么还需要先删除节点存放的其他数据。

1
2
3
4
5
6
7
8
9
10
11
12
void DestroyList(List *list)
{
Node *current;
Node *next;
current = list->head;
while (current->next != NULL){
next = current->next;
free(current);
current = next;
}
free(list);
}

我们还可以定一个输出函数,将链表里存放的数据依次输出

1
2
3
4
5
6
7
8
9
//打印整个链表
void dump(List *list){
Node *node;
node = list->head;
while (node != NULL){
printf("%08d\n", node->data);
node = node->next;
}
}

有了上面的基本函数时候,我们就能够读取存放数字的文本,将其加入到链表中。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
int main(int argc, char const *argv[])
{
/* code */
if (argc == 1) exit(EXIT_FAILURE);

FILE *fp;
fp = fopen(argv[1], "r");
if (fp == NULL){
perror(argv[1]);
exit(EXIT_FAILURE);
}
int data;
List *list;
list = CreateList();
while (fscanf(fp, "%d", &data) != EOF){
AddNode(list, data);
}
dump(list)
return 0;

我们的链表还应该支持插入操作和删除操作。对于插入操作,我们要分为是插入到给定位置前,还是给定位置后。对于删除而言,也就是都是删除当前节点,而为了删除当前节点,我们需要前一个节点的位置。

无论是插入还是删除,我们都需要知道插入的位置和删除的位置,因此我们还需要一个搜索函数,用于搜索等于给定值的节点位置或者是上一个位置。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// 查找元素
// situ=true时, 返回当前位置, false, 则返回上一个位置
Node *Search(List *list, int data, bool situ)
{
Node *node;
node = list->head;
if ( situ ){
while ( node->next != NULL ){
if ( node->data == data)
return node;
node = node->next;
}
} else {
while ( node->next->next != NULL) {
if (node->next->data == data)
return node;
node = node->next;
}
}
return NULL;
}

我们先写一个删除操作, 用于删除等于给定的节点。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
//删除节点
bool DeleteNode( List* list, int data){

Node *node;
Node *tmp;
node = list->head;
// 判断这个节点是否是首节点
if ( node->data == data ){
free(list->head);
list->head = NULL;
list->tail = list->head;
list->num = 0;
return true;
}
// 查找给定节点的前一个节点
node = Search(list, data, false);
// 找不到节点
if ( node == NULL){
return false;
}
//删除
tmp = node->next->next;
free(node->next);
node->next = tmp;
return true;
}

然后将元素加入函数分为两种,一种是插入(当前位置前),一种是追加(当前位置后)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
//在给定元素前加节点
bool InsertNode( List* list, int query, int data){

Node *node;
Node *new_node;

// 为新节点分配内存
new_node = (Node *)malloc( sizeof(Node) );
if ( new_node == NULL) return false;
new_node->data = data;

node = list->head;
// 判断这个节点是否是首节点
if ( node->data == query ){
new_node->next = node->next ;
node->next = new_node;
return true;
}

// 查找给定节点的前一个节点
node = Search(list, query, false);
// 找不到节点
if ( node == NULL){
return false;
}
new_node->next = node->next ;
node->next = new_node;

return true;

}

//在给定元素后加
bool AppendNode( List* list, int query, int data){

Node *node;
Node *new_node;

// 为新节点分配内存
new_node = (Node *)malloc( sizeof(Node) );
if ( new_node == NULL) return false;
new_node->data = data;

// 查找给定节点的位置
node = Search(list, query, true);
// 找不到节点
if ( node == NULL){
return false;
}
new_node->next = node->next;
node->next = new_node;

return true;

}

进阶操作

上面都是链表的基础操作,创建、摧毁,增加,删除。下面几个则是考验对链表对深刻理解,

  • 单链表反转
  • 链表中环的检测
  • 两个有序链表的合并
  • 删除链表倒数第N个结点
  • 求链表的中间结点

单链表反转

如果要将单链表进行反转,每次移动的时候需要三个位置,前一个位置,当前位置和head。每次将head向后移动,记录了当前位置的下一个节点,然后将当前位置指向前一个位置。最后将前一个位置和当前位置向后移动。图解如下, 首先head指向链表第一个节点

head赋值

然后将cur设置到当前的head

赋值cur
接着将head往后移动一个位置, 保存了原本在cur后面的位置

head后移

然后将cur指向到res,也就是前面的位置

cur指向res

上面的操作后,就将res和cur的顺序反转了。接着就是将res和cur往后移动

移动cur和res

代码为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
List* reverseList(List* list){

Node *curr, *res;
res = NULL;
curr = list->head;
//尾巴是之前的开头
list->tail = list->head;
while ( curr ){
//移动head
list->head = list->head->next;
//将当前位置指向前一个位置
curr->next = res;
//依次向后移动res和curr
res = curr;
curr = list->head;
}
list->head = res;
return list;
}

中间节点

为了寻找中间节点,我们可以定义两个指针,快指针和慢指针。慢指针一次一步,快指针一次两步. 如果是偶数,那么快指针最后是NULL,如果是奇数,那么快指针的下一个是NULL。

1
2
3
4
5
6
7
8
9
10
11
12
Node *FindMidlle(List *list)
{
if (list->num == 0) return NULL;
Node *fast = list->head;
Node *slow = list->head;
while ( fast != NULL && fast->next != NULL){
slow = slow->next;
fast = fast->next->next;
}
return slow;

}

删除倒数第N个指针

同上,也是快慢两个指针,快指针先走N步,然后两个指针再一起走。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
bool RemoveLastN(List *list, int n)
{
//删除第一个
if ( list->num == n){
list->head = list->head->next;
return true;
}
Node *fast = list->head;
Node *slow = list->head;
Node *tmp;
while (n-- > 0){
fast = fast->next;
}
while (fast->next != NULL){
fast = fast->next;
slow = slow->next;
}
tmp = slow->next;
slow->next = slow->next->next;
free(tmp);
return true;
}

有序链表合并

假设两个有序链表分别为1->3->5->7->82->3->4->5->8, 那么合并之后应该是1->2->3->3->4->5->7->8.

我们需要创建一个新的链表用于存放两个链表排序的结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
//合并两个链表
List *MergeSortedList(List *list_a, List *list_b)
{
List *list_c;
list_c = CreateList();
Node *node_a, *node_b, *node_c;
node_a = list_a->head;
node_b = list_b->head;

//确定新列表的head
if ( node_a->data < node_b->data ){
list_c->head = node_a;
node_a = node_a->next;
} else {
list_c->head = node_b;
node_b = node_b->next;

}
node_c = list_c->head;
while( true ){
if (node_a->data < node_b->data){
node_c->next = node_a;
node_a = node_a->next;
if (node_a == NULL) break;
} else{
node_c->next = node_b;
node_b = node_b->next;
if (node_b == NULL) break;
}
node_c = node_c->next;
}
while ( node_a != NULL){
node_c->next = node_a;
node_a = node_a->next;
node_c = node_c->next;

}
while ( node_b != NULL){
node_c->next = node_b;
node_b = node_b->next;
node_c = node_c->next;
}
return list_c;
}

为了测试这个代码正确性,我写了一个测试函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
int MergeTest( const char *file1, const char *file2){
FILE *f1;
FILE *f2;

int data;
f1 = fopen(file1, "r");

List *list1;
list1 = CreateList();
//读取数据
while (fscanf(f1, "%d", &data) != EOF){
AddNode(list1, data);
}
dump(list1);
fclose(f1);

f2 = fopen(file2, "r");
if (f2 == NULL){
perror(file2);
exit(EXIT_FAILURE);
}
List *list2;
list2 = CreateList();

//读取数据
while (fscanf(f2, "%d", &data) != EOF){
AddNode(list2, data);
}
dump(list2);
fclose(f2);

List *res;
res = MergeSortedList(list1, list2);
dump(res);
return 0;

}

最终的代码在GitHub上https://github.com/xuzhougeng/learn-algo/blob/master/link_list.c

LeetCode和链表有关的几个题目

Klib之khash学习笔记

散列表(Hash table,也叫哈希表),是根据关键码值(Key value)而直接进行访问的数据结构。也就是说,它通过把关键码值映射到表中一个位置来访问记录,以加快查找的速度。这个映射函数叫做散列函数,存放记录的数组叫做散列表。
给定表M,存在函数f(key),对任意给定的关键字值key,代入函数后若能得到包含该关键字的记录在表中的地址,则称表M为哈希(Hash)表,函数f(key)为哈希(Hash) 函数。 –来源:百度百科

klib提供的khash.h的初始化方法分为两种数据结构,分别是SET和MAP。SET只有键,且键唯一,MAP有键和值,键唯一,而值不唯一。

SET和MAP分别有三种初始化方法,对应键的类型分别为INT,INT64STR,而哈希算法也分为数值和字符串两类

1
2
3
4
5
6
7
8
9
10
11
12
13
14
//SET
#define KHASH_SET_INIT_INT(name) \
KHASH_INIT(name, khint32_t, char, 0, kh_int_hash_func, kh_int_hash_equal)
#define KHASH_SET_INIT_INT64(name) \
KHASH_INIT(name, khint64_t, char, 0, kh_int64_hash_func, kh_int64_hash_equal)
#define KHASH_SET_INIT_STR(name) \
KHASH_INIT(name, kh_cstr_t, char, 0, kh_str_hash_func, kh_str_hash_equal)
//MAP
#define KHASH_MAP_INIT_INT(name, khval_t) \
KHASH_INIT(name, khint32_t, khval_t, 1, kh_int_hash_func, kh_int_hash_equal)
#define KHASH_MAP_INIT_INT64(name, khval_t) \
KHASH_INIT(name, khint64_t, khval_t, 1, kh_int64_hash_func, kh_int64_hash_equal)
#define KHASH_MAP_INIT_STR(name, khval_t) \
KHASH_INIT(name, kh_cstr_t, khval_t, 1, kh_str_hash_func, kh_str_hash_equal)

键值对中的kint32_tkhin64_t和系统有关,用于定义一个很大的取值范围。

1
2
3
4
5
6
7
8
9
10
11
#if UINT_MAX == 0xffffffffu
typedef unsigned int khint32_t;
#elif ULONG_MAX == 0xffffffffu
typedef unsigned long khint32_t;
#endif

#if ULONG_MAX == ULLONG_MAX
typedef unsigned long khint64_t;
#else
typedef unsigned long long khint64_t;
#endif

kh_cstr_t的定义是typedef const char *kh_cstr_t;, 是一个不会变的字符串。

这两种类型用于设置KHASH_INIT的参数khkey_tkhval_t, 用于初始化哈希表的结构定义

1
2
3
4
5
6
7
8
#define __KHASH_TYPE(name, khkey_t, khval_t) \
typedef struct kh_##name##_s { \
//桶的数目, 哈希表大小, 占用数, 上限
khint_t n_buckets, size, n_occupied, upper_bound; \
khint32_t *flags; \ //记录当前位置是否被使用
khkey_t *keys; \ //键的类型
khval_t *vals; \ //值的类型
} kh_##name##_t;

和哈希表操作有关的函数如下

  • kh_init(name): 初始化哈希表
  • kh_destroy(name, h); 删除哈希表
  • kh_clear(name, h): 保持哈希表大小不变,清空内容
  • kh_resize(name, h, s): 调整哈希表大小, 运行时它会被自动调用,用于扩容
  • kh_put(name, h, k, r): 将key放在哈希表中,并获取key的位置
  • kh_get(name, h, k): 获取key对应的位置
  • kh_del(name, h, k): 删除哈希表元素
  • kh_exist(h, x): 检查哈希表位置上是否有内容
  • kh_key(h, x): 获取哈希表中x对应的key
  • kh_value(h,x): 获取哈希表中键x的值
  • kh_begin(h): 获取哈希表的起始key
  • kh_end(h): 获取哈希表的最后key
  • kh_size(h): 获取哈希表的大小
  • kh_n_buckets(h): 哈希表中桶的数目
  • kh_foreach(h, kvar, vvar, code): 遍历哈希表,其中键赋值给kvar, 值赋值给vvar,运行code的代码
  • kh_foreach_value(h, vvar, code): 遍历哈希表,其中值赋予给vvar,运行code的代码

为了达到类似于Python的字典操作,例如d = {"abc": "aaa"}d["abc"],所需要写的代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#include <stdio.h>
#include <string.h>
#include "klib/khash.h"

KHASH_MAP_INIT_STR(dict, char*)

int main(int argc, char *argv[])
{

khiter_t k;

khash_t(dict) *h = kh_init(dict);

int ret;

k = kh_put(dict, h, "abc", &ret);
kh_key(h, k) = strdup("abc");
kh_value(h,k) = strdup("aaa");

k = kh_put(dict, h, "efg", &ret);
kh_key(h, k) = strdup("efg");
kh_value(h,k) = strdup("bbb");

k = kh_get(dict, h, "abc");
printf("%s\n", kh_value(h,k));

for ( k = kh_begin(h); k != kh_end(h) ; k++){
if (kh_exist(h,k)){
free((char*)kh_key(h, k));
free((char*)kh_value(h, k));
}
}

kh_destroy(dict, h);

}

对于字符串,哈希表结构中的keysvals并不存放实际的值,而是存放字符串的地址,因此如果没有专门内存用于存放键值对对字符串,那么用strdup在内存中获取字符串新的地址。 如果用一个字符串数组存放键值对字符串的地址,代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#include <stdio.h>
#include <string.h>
#include "klib/khash.h"

KHASH_MAP_INIT_STR(dict, char*)

int main(int argc, char *argv[])
{

khiter_t k;

khash_t(dict) *h = kh_init(dict);

int ret;

char *key[] = {"abc", "efg"};
char *value[] = {"aaa", "bbb"};

// h = {"abc":"efg"}
k = kh_put(dict, h, "abc", &ret);
kh_key(h, k) = key[0];
kh_value(h,k) = value[0];

// h["efg"] = "bbb"
k = kh_put(dict, h, "efg", &ret);
kh_key(h, k) = key[1];
kh_value(h,k) = value[1];

// h["abc"]
k = kh_get(dict, h, "abc");
printf("%s\n", kh_value(h,k));
kh_destroy(dict, h);

}

推荐阅读

数据结构之堆(heap)

极客时间的「数据结构与算法之美」的学习笔记,图片来源于「28 | 堆和堆排序:为什么说堆排序没有快速排序快?」

堆满足两个要求:

  1. 完全二叉树
  2. 父节点的元素大于(或小于)子节点的元素

堆的实现

为了实现一个堆,我们需要创造一个堆的数据结构,以及实现堆的插入和删除等操作函数。

堆的存储

由于堆是完全二叉树,因此可以用数组存放堆。第i个节点就放在数组的第i个位置上。它的左子节点是 2i, 它的右子节点是2i+1, 它的父节点是i/2.

堆的存放

这里定义了一个堆的结构,包含三个元素

1
2
3
4
5
typedef struct _heap {
int *heap; //指向存放堆对数组
int n; //堆的大小
int count; //堆目前的元素
} heap;

创建堆的函数一开始写的代码如下

1
2
3
4
5
6
7
8
9
heap*
createHeap(int n){

heap *h;
h->heap = (int*)malloc( sizeof(int) * (n+1));
h->n = n;
h->count = 0;
return h;
}

代码运行会出错,因为heap *h只是声明了一个变量,并没有分配一个内存空间用于构造一个heap结构,同时将h指向这个内存地址。 因此应该加一句,h = (heap*)malloc( sizeof(heap) );. 也就是

1
2
3
4
5
6
7
8
9
10
heap*
createHeap(int n){

heap *h;
h = (heap*)malloc( sizeof(heap) );
h->heap = (int*)malloc( sizeof(int) * (n+1));
h->n = n;
h->count = 0;
return h;
}

堆的操作

堆有两个最常用堆操作,插入元素和删除顶部元素。无论是何种操作,都需要保证操作之后的数据依旧满足堆的两个特性,也就是堆化(heapify)。

每次往堆里加入一个元素,其实是放在数据的最后一个位置,然后让加入元素继续满足堆的性质。

堆的插入

最开始写的插入代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
void
insertElement(heap *h, int e){
if (h->n >= h->count) return;

// 在数组最后加入新的元素
int *heap = h->heap;
int count = h->count;
heap[++count] = e;

printf("Insert %d at %d\n", e,count+1);
h->count = count;

int i = count;// 数组索引
//堆化,
while( (i/2) > 0 && heap[i] > heap[i/2] ){
swap(heap, i, i /2); //交换父子节点
}
}

上面代码在调试的时候发现,函数没有输出”Insert %d at %d\n”这一不会显示,也就是说前面的if语句就无法顺利运行。仔细检查发现是逻辑语句写反了, 应该是(h->count >= h->n)。更改此处错误之后,发现堆化依旧失败,原因是while语句中,每次循环中缺少一句i=i/2,导致循环之后结果不正确。

删除堆顶元素有两种方式。一种是直接删除第一个元素,然后开始堆化,但是写代码比较复杂,很可能产生一个非完全二叉树。第二个方式是删除第一个元素,并用最后一个元素替换。然后至上而下进行堆化

堆的删除操作

代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
void removeTop(heap *h){

if (h->count < 1) return ;

int count = h->count;
int *heap = h->heap;
//用最后一个元素替代第一个元素
heap[1] = heap[count];
//删除最后一个元素
h->count = --count;
// 自上而下堆化
int i = 1;
while ( true ){
int max_pos = i;
if ( i*2 <=count && heap[i] < heap[i*2]) max_pos = i*2; //和左子节点比较
if ( i*2+1<=count && heap[i] < heap[i*2+1]) max_pos = i*2+1; //和右子节点比较
if (max_pos == i) break; // 不再发生交换, 当前位置就是最大位置
swap(heap, i, max_pos); //将当前节点和子节点进行交换
i = max_pos ; //将i设置为子节点的索引
}

}

最后代码参考https://github.com/xuzhougeng/learn-algo/blob/master/heap.c

C语言实战课-klib库学习

对于Python和R而言,使用一个已有的工具,通常都是先安装,然后看帮助文档学会怎么调用函数。那么对于C语言来说,我应该如何使用一个已有的轮子呢?

前些日子,我翻译了李恒大神一篇关于seqtk代码介绍的博客,让我决定以seqtk为切入点,介绍下如何用别人造好的轮子。seqtk主要用到了两个头文件, khash.h和kseq.h, 同时我发现这两个文件在klib中,而klib提供了较为详细的介绍文档,因此我就直接去学klib了。

关于klib名字中的k, 我想到了之前英语老师说,在英语中k和c的发音类似,例如kindle读起来就像candle, 因此这里的k以及函数名,结构体名里的k,都可以认为是c,表示这是一个c语言库。

第一步:阅读文档

一个优秀的项目必然会有一个优秀的文档对这个项目进行介绍,否则就只能自娱自乐,不被广泛使用。因此,学习klib的第一步就是阅读它的文档。我为了让自己更好的理解,因此通过翻译的方式让自己能够认真阅读。

为了实现范型的容器,klib广泛的使用了C宏。为了使用这些数据结构,我们需要先扩展出宏中的实例方法。于是,源代码就显得有点难以阅读以及不利于后续调试。由于C语言没有模版特性,实现高效范型编程也只能用宏。只有在宏的帮助下,我们才能写出范型容器,它在初始化之后才能在效率上和类型特异的容器竞争。其他一些库,例如Glib, 使用void *类型实现容器。但是这种实现方法效率比较低,并且比klib使用更多内存。(参考这个评测结果

为了有效的使用klib,我们需要了解下它是如何实现范型编程。这里哈希表库(哈希表是一种键值对数据库,查询效率高)为例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#include "khash.h"
KHASH_MAP_INIT_INT(32, char) //利用宏构建函数和数据结构
int main() {
int ret, is_missing;
khiter_t k;
khash_t(32) *h = kh_init(32); //初始化
k = kh_put(32, h, 5, &ret); //插入key
kh_value(h, k) = 10; //根据key设置值
k = kh_get(32, h, 10); //根据值找key
is_missing = (k == kh_end(h));
k = kh_get(32, h, 5);
kh_del(32, h, k); //删除key
for (k = kh_begin(h); k != kh_end(h); ++k) //遍历哈希表
if (kh_exist(h, k)) kh_value(h, k) = 1; //判断桶是否有数据
kh_destroy(32, h); //删除哈希表
return 0;
}

在这里例子中,第二行实例化了一个键类型为unsigned,值类型为char的哈希表,其中m32是哈希表的名字。所有和该名字相关的函数和数据类型都是宏,之后会进行解释。宏kh_init()初始化哈希表,而kh_destroy()则负责释放内存。kh_put插入键,并返回一个哈希表的迭代器(或者说位置)。kh_get()kh_del()分别是获取值对应的键以及删除里面的元素。宏kh_exist()检测迭代器(或者说位置)是否在数据中。

看完代码之后,你会发现代码看起来不太像是有效的C程序,例如它缺少分号, 对一个函数赋值,m32没有预先定义。为了理解为什么代码是正确都,让我们更进一步都看下khash.h都源代码,代码骨架如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#define KHASH_INIT(name, SCOPE, key_t, val_t, is_map, _hashf, _hasheq) \
typedef struct { \
int n_buckets, size, n_occupied, upper_bound; \
unsigned *flags; \
key_t *keys; \
val_t *vals; \
} kh_##name##_t; \
SCOPE inline kh_##name##_t *init_##name() { \
return (kh_##name##_t*)calloc(1, sizeof(kh_##name##_t)); \
} \
SCOPE inline int get_##name(kh_##name##_t *h, key_t k) \
... \
SCOPE inline void destroy_##name(kh_##name##_t *h) { \
if (h) { \
free(h->keys); free(h->flags); free(h->vals); free(h); \
} \
}

#define _int_hf(key) (unsigned)(key)
#define _int_heq(a, b) (a == b)
#define khash_t(name) kh_##name##_t
#define kh_value(h, k) ((h)->vals[k])
#define kh_begin(h, k) 0
#define kh_end(h) ((h)->n_buckets)
#define kh_init(name) init_##name()
#define kh_get(name, h, k) get_##name(h, k)
#define kh_destroy(name, h) destroy_##name(h)
...
#define KHASH_MAP_INIT_INT(name, val_t) \
KHASH_INIT(name, static, unsigned, val_t, is_map, _int_hf, _int_heq)

KHASH_INIT是一个巨大都宏,他定义了所有都结构体方法。当这个宏被调用时,所有代码会通过C预处理器插入到程序中被调用都位置。如果一个宏被多次调用,那么代码会在程序中出现多个拷贝。为了避免哈希表中因为不同键-值类型导致都命名冲突,这个库使用token concatenation, 也就是预处理器能够使用宏里的参数替换符号中部分内容。最后C预处理器会产生下面的代码,然后传递给编译器。(由于kh_exist(h,k)比较复杂,这里就不展开介绍)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
typedef struct {
int n_buckets, size, n_occupied, upper_bound;
unsigned *flags;
unsigned *keys;
char *vals;
} kh_m32_t;
static inline kh_m32_t *init_m32() {
return (kh_m32_t*)calloc(1, sizeof(kh_m32_t));
}
static inline int get_m32(kh_m32_t *h, unsigned k)
...
static inline void destroy_m32(kh_m32_t *h) {
if (h) {
free(h->keys); free(h->flags); free(h->vals); free(h);
}
}

int main() {
int ret, is_missing;
khint_t k;
kh_m32_t *h = init_m32();
k = put_m32(h, 5, &ret);
if (!ret) del_m32(h, k);
h->vals[k] = 10;
k = get_m32(h, 10);
is_missing = (k == h->n_buckets);
k = get_m32(h, 5);
del_m32(h, k);
for (k = 0; k != h->n_buckets; ++k)
if (kh_exist(h, k)) h->vals[k] = 1;
destroy_m32(h);
return 0;
}

从这个例子中,我们可以看到宏和C预处理器是klib中非常重要的部分。为什么klib那么快,一部分原因就是编译器在编译过程中已经知道键值对类型,因此他就能将代码优化到类特异代码对水平。一个用void *写的范型库无法达到这种性能。

在实例化过程中会有大量代码插入,这也提示我们为什么C++编译速度慢以及使用STL/boost编译的二进制文件很大。Klib由于代码量很小,并且各个组件相对独立,因此表现更加优异。插入上百行代码不会让编译变得很慢。

第二步:运行案例

介绍文档给了哈希表作为案例,我们就可以复制代码进行测试

1
git clone https://github.com/attractivechaos/klib

新建一个hash_test.c, 粘贴下列内容

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#include "klib/khash.h"
KHASH_MAP_INIT_INT(m32, char) // instantiate structs and methods
int main() {
int ret, is_missing;
khint_t k;
khash_t(m32) *h = kh_init(m32); // allocate a hash table
k = kh_put(m32, h, 5, &ret); // insert a key to the hash table
if (!ret) kh_del(m32, h, k);
kh_value(h, k) = 10; // set the value
k = kh_get(m32, h, 10); // query the hash table
is_missing = (k == kh_end(h)); // test if the key is present
k = kh_get(m32, h, 5);
kh_del(m32, h, k); // remove a key-value pair
for (k = kh_begin(h); k != kh_end(h); ++k) // traverse
if (kh_exist(h, k)) // test if a bucket contains data
kh_value(h, k) = 1;
kh_destroy(m32, h); // deallocate the hash table
return 0;
}

编译和运行

1
2
gcc -o test hash_test.c
./test

这个代码在运行时不会有任何输出,下一步我们可以尝试改造代码,让代码能够输出信息。

第三步: 改造代码

学会写代码的最快路径就是多写代码。以前我学习的时候,总觉得看一遍就够了,没有必要真的写一遍。真正开始写代码的时候,发现自己啥都不会。因此,为了能够学会使用khash.h, 我们按照自己的想法对代码做一些改造,看看代码的表现是否符合预期。下面代码,我写了一个循环,用于设置键值对。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#include <stdio.h>
#include "klib/khash.h"
KHASH_MAP_INIT_INT(m32, char) // instantiate structs and methods
int main() {
int ret, is_missing;
khint_t k;
khash_t(m32) *h = kh_init(m32); // allocate a hash table

int i = 0;
for ( i = 0 ; i < 10 ; i++ ){
k = kh_put(m32, h, i, &ret); // insert a key to the hash table
if (!ret) kh_del(m32, h, k);
kh_value(h, k) = i * i; // set the value
}

for (k = kh_begin(h); k != kh_end(h); ++k) // traverse
if (kh_exist(h, k)) // test if a bucket contains data
printf("%d\n", kh_value(h, k) );
kh_destroy(m32, h); // deallocate the hash table
return 0;
}

假如你有个文件,也是键值对的关系,那么你可以读取这个文件,用文件里对内容对哈希表进行赋值。

第四步: 在项目中应用

学习klib的目的是为了解决我实际的问题,也就是写一个通用的fastq转换成fasta程序,这个程序能够自动解决行过长的问题。最终代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
#include <stdio.h>
#include <zlib.h>
#include "klib/kseq.h"

KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[])
{

gzFile fp;
kseq_t *seq;
int l;
if (argc == 1){
fprintf(stderr, "Usage: %s <in.fasta|in.fasta.gz>\n", argv[0]);
return 1;
}

fp = gzopen(argv[1], "r");
seq = kseq_init(fp); // 分配内存给seq
while( (l = kseq_read(seq)) >= 0){ //读取数据到seq中
printf(">%s\n", seq->name.s); //
printf("%s\n", seq->seq.s);
}

kseq_destroy(seq); //释放内存
gzclose(fp);

return 0;


}

这次的代码就特别的简洁,并且容易阅读。除了多了几行变量名定义和内存管理的代码外,其他时间就和写Python脚本一样,无非就是调用函数,kseq_read用于数据读取,kseq_t用于存放序列数据, 它的结构定义如下

1
2
3
4
5
6
7
8
9
typedef struct __kstring_t {
size_t l, m;
char *s;
} kstring_t;
typedef struct {
kstring_t name, comment, seq, qual;
int last_char;
kstream_t *f;
} kseq_t

这是我写的关于klib库的第一篇教程,往后还会不断用到这个库,用到的时候再写几篇进行介绍。

C语言实战课-让程序能够处理gz文件

之前的程序不能够处理压缩文件,而事实上为了节约空间,基本上fastq都会压缩成gz格式,因此这一课就是让程序能够支持压缩文件。

这里有两种思路,一种是利用管道,将之前的压缩文件通过zcat程序读取然后传递给我们的程序,另一种则是在程序中调用zlib库,让程序能够直接处理gz文件.

背景知识

针对思路1: 系统在每个C语言运行的时候都至少会提供三个流,标准输入(stdin),标准输出(stdout)和标准错误(stderr). 在之前的程序我们就用到了stdout,用于将结果输出到屏幕上, 即fprintf(stdout, "%s", line);。 同样,我们可以修改fgets(line, MAX_LINE_LENGTH, fi) 中的fistdin使得程序能够接受管道传递的数据。

针对思路2: 利用zlib读取gz文件并不复杂,只需要将原来的函数名前面加上或者改成gz。毕竟优秀的代码应该符合人的直觉,下面就是我们将要用到的几个函数

  • 声明水管: gzFile
  • 连接水管: gzopen(const char *path, const char *mode)
  • 读取一行: gzgets(gzFile file, char *buf, int len)
  • 输出到gz文件中: gzprintf(gzFile file, const char *format, ...)
  • 关闭水管: gzclose(gzFile file)

我们以一个非常简单的代码进行展示,它会从gz文件中按行读取然后输出

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#include <stdio.h>
#include <zlib.h>

int main(int argc, char *argv[])
{

gzFile fi;
fi = gzopen(argv[1], "r");
char string[1000];

while ( gzgets(fi, string, 1000) != NULL) {
printf("%s", string);
}
return 0;
}

它不是只能处理压缩文件,对于没有压缩的文件也是可以识别的,因此不需要写专门的逻辑判断语句,根据文件名来决定是否使用zlib相关的函数。

gzopen can be used to read a file which is not in gzip format; in this case gzread will directly read from the file without decompression. When reading, this will be detected automatically by looking for the magic two-byte gzip header

由于zlib不是标准库,我们在编译程序的时候需要额外加上-lz参数,用于链接动态库。

动手搞事情

有了以上的背景知识后,我们就可以来修改我们之前的代码了.

zlib在大部分服务器上都属于默认安装,因此只需要多加一行#include <zlib.h>. 但是还是有一小部分服务器居然没有安装zlib,因此我们需要先下载编译(root权限可以直接用apt/yum进行安装)。

1
2
3
wget https://www.zlib.net/zlib-1.2.11.tar.gz
tar xf zlib-1.2.11.tar.gz
cd zlib-1.2.11 && ./configure && make -j 8

修改后的fq2fa.c的代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "zlib.h"

#define MAX_LINE_LENGTH 1000

int main(int argc, char *argv[]){

gzFile fi;
gzFile fo;
if (argc < 2 ){
return -1;
}
fi = gzopen( argv[1], "rb");
if ( argc == 3 ){
fo = gzopen (argv[2], "wb");
}

char line[MAX_LINE_LENGTH];
int count = 0;
while ( gzgets(fi,line, MAX_LINE_LENGTH) != NULL) {
if ( count == 0) {
line[0] = '>';
}
if ( count < 2){
if ( argc == 2) {
fprintf(stdout, "%s", line);
}
if ( argc ==3 ){
gzprintf(fo, "%s", line);
}
}
if ( ++count > 3 ){
count = 0;
}

}
gzclose(fi);
if (argc == 3) gzclose(fo);
return 0;
}

编译方法为gcc -lz -Lzlib -o fq2fa fq2fa.c。关于代码的讲解和zlib的知识点可以看视频介绍, https://www.bilibili.com/video/av84768292/

代码目前存在的问题是行缓冲是1000字符,对于三代测序数据,按照宣传可以有几M,一般也是几百k,因此当前程序处理三代程序肯定不行。一个方法是将MAX_LINE_LENGTH设置为更大的值,比如说设置为10M,或者是将其设置为一个参数,可以手动指定。

不过下一次,将会学习一个新的轮子klib,用里面的函数来解决这个问题。

C语言实战课-将FASTQ转成FASTA

我这次学习C语言的目的非常单纯,就是尝试将C语言应用在日常的分析任务重,解决实际问题。既然如此,那么第一课就不是打印”hello world!” 了,毕竟说了那么多次世界你好,依旧写不好代码。

我们第一课直接就来处理一个实际的小需求,读取FASTQ,将其转成一个FASTA。处理这个问题和把大象放进冰箱里一样,都是分为三步,读取数据,处理数据,输出数据。其中第一步和第三步都是和文件打交道,而第二步考验的是对算法,数据结构和内存等有关知识对理解。

文件读写背景知识

我们应该如何读取数据呢?如果R语言,我会用readLines("input.fasta")直接读取所有的数据。如果是我以前用的python, 代码会是下面的样子

1
2
3
4
fi = open("input.fasta", "r")
for line in fi.readline():
print(line)
fi.close()

乍看起来,R好像比Python写的代码要少很多。但其实Python也可以一行完成,[line for line in open("input.fasta", "r"]。将数据输出到外部也差不多,基本上都能一行命令搞定。

总之,当你发现自己的代码可以少写了,其实是有人帮你简化了。

回归到C语言,所有的读写操作其实就是简单地将字节逐个移入程序中,将字节逐个从程序中移出而已,类似于水动的感觉。

我们要让程序读取数据和输出数据,就相当于要搞一个水管,把水引过来。下面代码就是用FILE *声明了fi和fo这两根水管

1
2
FILE *fi;
FILE *fo;

然后我们还需要把水管接到水流上。下面代码的fopen就是将我们之前对两根水管分别接到了输入文件和输出文件上,参数中的”r”表示读取(read), 而”w”表示写出(write)。

1
2
3
4
5
6
7
8
9
fi = fopen("input.txt", "r");
fo = fopen("output.txt", "w");

//确保文件能够正常打开
//否则退出
if ( fi == NULL ){
perror(argv[1]);
exit(EXIT_FAILURE);
}

水管接上了,那么如何打开水阀,让水不断地流入流出呢?一种是逐个字符操作,一种是逐行操作。前面的方法更加基础,因为逐行读取就是不断的读取一个个字符,当碰到换行符时,就把前面的字符合在一起,以字符串的形式传递,输出是以行为单位,每一行也是一个个字符的写出。

  • 逐字符函数: fgetc, fputc
  • 逐行函数: fgets, fputs

通常读取和写出都伴随着数据处理,因此这部分放到实际数据处理进一步介绍。

最后关闭水龙头,我们需要用到fclose()

1
2
fclose(fi);
fclose(fo);

动手搞事情

有了一些文件读写的基本概念之后,我们就可以开始动手写代码了。

根据fastq和fasta的格式定义,fasta是一行以>开头的序列表示,之后跟着N条序列。而fastq则是标准的4行,第一行是以@开头的序列标示符,第二行是序列,第三行以+开头,第四行是第二行序列对应的质量信息。处理过程描述如下

  • 读取第一行
  • 将第一个字符的@替换成>
  • 输出第一行
  • 读取第二行,输出第二行
  • 读取第三行,不输出
  • 读取第四行,不输出

下面实际的代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#include <stdio.h>
#include <stdlib.h>

#define MAX_LINE_LENGTH 1000

int main(int argc, char *argv[])
{
FILE *fi;
FILE *fo;

if ( argc == 2) {
fi = fopen(argv[1], "r");
} else if (argc == 3){
fo = fopen(argv[2], "w");
} else {
exit(EXIT_FAILURE);
}

char line[MAX_LINE_LENGTH];
int count=0;

while ( fgets(line, MAX_LINE_LENGTH, fi) != NULL ){
if (count == 0){
line[0] = '>';
}
if (count < 2){
if ( argc == 2){
fprintf(stdout, "%s", line);
} else{
fprintf(fo, "%s", line);
}
}
if (++count > 3){
count = 0;
}
}

fclose(fi);
if (argc == 3) fclose(fo);

return 0;

}

关于代码的讲解,我录制了专门的视频,https://www.bilibili.com/video/av84245773/

具体使用方法为

1
./fq2fa input.fq output.fa

目前代码代码比较简陋,至少存在下面这些问题

  • 不支持gz压缩
  • 不支持管道输入
  • 没有文件类型检验,不能提示正确的错误

根据这些需求,我们就有了学习的目标,也就是后面课程的内容。