A*算法lyh's blog

简要介绍两类启发式搜索算法：A*和IDA*算法。

众所周知，DFS和BFS是两类盲目搜索算法。所谓盲目，就是根据当前状态毫无目标性地进行扩展。想像一下，如果你正在操场上，你如何找到一条尽可能短的路走向旗杆？显然根据视觉上的提示，旗杆所在的方向就是行进的方向，很快就找到了一条路径。如果蒙上眼睛，你可能就会漫无目的的随处寻找了。这种视觉上的提示就称为“启发”，普通的搜索算法没有启发的引导，因此搜索过程冗杂，需要各种剪枝，效率比较低下。如果能够给搜索过程加上合适的启发方向，那么搜索过程将会大大简化，效率也会得到提升。这里所说的A*和IDA*算法就是两种启发式搜索算法。

A*算法

如何给状态设置搜索方向呢？A*算法引入了估价函数的概念。
对于状态n，我们可以求出从这个状态开始到目标状态的预估代价h(n)，这个h(n)函数称为估价函数。另外对于每一个状态，有到这个状态的实际代价g(n)。那么就可以求出从初始状态到目标状态的预估代价f(n)，它满足：

$f(n)=g(n)+h(n)$

估价函数的选取稍后再提，这里先介绍A*算法流程。
A*算法中存在一个开启列表和关闭列表，开启列表中储存当前转移的状态，关闭列表中储存当前已经经过转移的状态。一个状态只要进入关闭列表就不会再重新进入开启列表。A*算法步骤如下：

将初始状态加入开启列表，设置其g和h函数值。
从开启列表中选择f值最小的一个状态，从这个状态开始转移。如果转移后的状态在关闭列表中，则忽略，如果在开启列表中，则检查该状态的g值是否可以更新，若可以则更新g值。如果转移后的状态不在开启列表也不再关闭列表中，则将其加入开启列表，计算其g值和h值。
当从开启列表中选择出状态正好是目标状态时，结束算法。如果开启列表已空，则不存在到目标状态的路径。

如果需要记录路径，在每一个状态上加入前驱指标即可递归找到路径。由于需要在开启列表中找最小值，故开启列表常常用堆或优先队列实现。
容易看出，A*算法很像是将BFS中的队列换成优先队列，它的算法思路也是比较清晰的。如果存在路径，A*算法一定可以找到该路径，但最优性与估价函数选取有关，设实际最小代价为d(n)，则有以下规律：

若$h(n) < d(n)$，只要找到路径，该路径一定是最优解，$h(n)$越小效率越差。
若$h(n)=d(n)$，只要找到路径，该路径一定是最优解，此时效率最高。
若$h(n)>d(n)$，可以找到路径，但不保证最优。

        可见估价函数的正确选取对效率和最优性有着很大的影响。
        容易看出，对于转移代价恒定的情况，若h(n)恒为0，则A*算法退化为BFS，对于代价不定的情况，若h(n)恒为0，则退化为Dijkstra算法。
        既然A*算法可以比Dijkstra算法更好地求最短路，为何不用该算法代替Dijkstra算法呢？这是因为在给定的图中，预估当前结点到目标结点的转移代价是十分困难的，A*算法在这个问题上不实用。但是，对于方格图上转移（比如八数码问题），A*算法就大有用武之地。
        看这样一道题。
        很明显是一道搜索题，可以考虑用A*算法解决。首先注意到如果想让特殊方格移动，空白方格必须在其附近。那么我们可以先抽离这些状态先进行BFS预处理。预处理后，将会得到数组dist[x][y][i][j]表示特殊方格在坐标(x，y)处，空白方格从其i方向移动到其j方向的最小代价（方向定义为0上1下2右3左）。
        对于每一组询问，用BFS预处理将空白结点移动到特殊方格附近的转移代价（最多4个，分别对应4个方向），这几个对应的状态就是初始状态，将它们加入开启列表，计算h值和g值（g值就是BFS预处理的转移代价）。
        问题来了，如何计算h值？在这种问题上，我们通常用曼哈顿距离法来估计转移代价。对于平面上两点，它们的曼哈顿距离定义为：

$d=|x_2-x_2|+|y_1-y_2|$

        显然以曼哈顿距离作估价函数一定满足前文所说的$h(n)\leq d(n)$，可以保证得到最优解。
        下面就可以开始进行A*算法了。用小根堆维护开启列表，和Dijkstra算法手写堆版本相似，这里要构造一个哈希表来记录每一个状态在堆中的位置，若未入堆置为0，在关闭列表中置为-1。转移时最多有4个转移方向，若空白方格在特殊方格左侧，则它可以转移到上、下、右侧，也可以直接与特殊方格交换。
        示例代码，比较长：

#include<bits/stdc++.h>
//定义几个宏来简化代码
#define OK(x, y) ((x)>=1&&(x)<=n&&(y)>=1&&(y)<=m&&op[(x)][(y)])
#define H(x, y) (abs((x)-tx)+abs((y)-ty))
#define SUM(x) (g[(x).s][(x).sx][(x).sy]+h[(x).s][(x).sx][(x).sy])
#define RK(x) (rk[tree[x].s][tree[x].sx][tree[x].sy])
using namespace std;

struct Point {
    int x, y;

    Point(int x, int y) : x(x), y(y) {}
};

struct Node {
    int s, sx, sy;

    Node() {}

    Node(int a, int b, int c) {
        s = a, sx = b, sy = c;
    }

} tree[31 * 31 * 31 * 31];

queue<Point> que;
int n, m, q, op[35][35], tx, ty, TMP[35][35], dict[35][35][4][4], size = 0;
const int mx[4] = {-1, 1, 0, 0}, my[4] = {0, 0, 1, -1};
int rk[4][31][31], g[4][31][31], h[4][31][31];

inline void BFS(int ex, int ey, int sx, int sy, int d) {
    memset(TMP, -1, sizeof(TMP));
    TMP[ex][ey] = 0, TMP[sx][sy] = 1;
    que.push(Point(ex, ey));
    while (!que.empty()) {
        Point p = que.front();
        que.pop();
        for (int i = 0; i < 4; i++) {
            int xx = p.x + mx[i], yy = p.y + my[i];
            if (OK(xx, yy) && TMP[xx][yy] == -1)TMP[xx][yy] = TMP[p.x][p.y] + 1, que.push(Point(xx, yy));
        }
    }
    if (d == -1)return;
    for (int i = 0; i < 4; i++)dict[sx][sy][d][i] = TMP[sx + mx[i]][sy + my[i]];
}

void solve(int x) {
    int x1 = x << 1, x2 = x << 1 | 1, minn = x;
    if (x1 <= size && SUM(tree[x1]) < SUM(tree[minn]))minn = x1;
    if (x2 <= size && SUM(tree[x2]) < SUM(tree[minn]))minn = x2;
    if (minn != x)swap(tree[x], tree[minn]), swap(RK(x), RK(minn)), solve(minn);
}

void up(int x) {
    if (x == 1)return;
    if (SUM(tree[x >> 1]) > SUM(tree[x]))swap(tree[x >> 1], tree[x]), swap(RK(x >> 1), RK(x)), up(x >> 1);
}

inline void add(Node p) {
    tree[++size] = p, rk[p.s][p.sx][p.sy] = size, up(size);
}

inline Node top() {
    Node p = tree[1];
    rk[p.s][p.sx][p.sy] = -1;
    if (size > 1)tree[1] = tree[size--], RK(1) = 1, solve(1);
    else size = 0;
    return p;
}

inline void clear() {
    size = 0, memset(rk, 0, sizeof(rk));
    memset(g, 127, sizeof(g));//置为无穷大
}

int main() {
    ios::sync_with_stdio(false);
    cin >> n >> m >> q;
    for (int i = 1; i <= n; i++)for (int j = 1; j <= m; j++)cin >> op[i][j];
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= m; j++)
            for (int z = 0; z < 4; z++)if (OK(i, j) && OK(i + mx[z], j + my[z]))BFS(i + mx[z], j + my[z], i, j, z);
    for (int i = 1; i <= q; i++) {
        int ex, ey, sx, sy, tx, ty, flag = 0;
        cin >> ex >> ey >> sx >> sy >> tx >> ty;
        BFS(ex, ey, sx, sy, -1), clear();
        if (sx == tx && sy == ty) {
            cout << 0 << endl;
            continue;
        }
        for (int j = 0; j < 4; j++) {
            if (OK(sx + mx[j], sy + my[j]) && TMP[sx + mx[j]][sy + my[j]] != -1) {
                g[j][sx][sy] = TMP[sx + mx[j]][sy + my[j]], h[j][sx][sy] = H(sx, sy), add(Node(j, sx, sy));
            }
        }
        while (size > 0) {
            Node p = top();
            if (p.sx == tx && p.sy == ty) {
                cout << g[p.s][p.sx][p.sy] << endl, flag = 1;
                break;
            }
            for (int j = 0; j < 4; j++) {//不改变特殊方格位置的转移
                int xx = p.sx + mx[j], yy = p.sy + my[j], G;
                if (j == p.s || !OK(xx, yy) || dict[p.sx][p.sy][p.s][j] == -1)continue;
                G = dict[p.sx][p.sy][p.s][j] + g[p.s][p.sx][p.sy];
                if (G < g[j][p.sx][p.sy]) {
                    if (rk[j][p.sx][p.sy] == -1)continue;
                    if (!rk[j][p.sx][p.sy]) {
                        g[j][p.sx][p.sy] = G, h[j][p.sx][p.sy] = H(sx, sy), add(Node(j, p.sx, p.sy));
                    } else g[j][p.sx][p.sy] = G, up(rk[j][p.sx][p.sy]);
                }
            }
            int G = g[p.s][p.sx][p.sy] + 1;//与特殊方格交换的转移
            if (G < g[p.s ^ 1][p.sx + mx[p.s]][p.sy + my[p.s]]) {
                if (rk[p.s ^ 1][p.sx + mx[p.s]][p.sy + my[p.s]] == -1)continue;
                if (!rk[p.s ^ 1][p.sx + mx[p.s]][p.sy + my[p.s]]) {
                    g[p.s ^ 1][p.sx + mx[p.s]][p.sy + my[p.s]] = G;
                    h[p.s ^ 1][p.sx + mx[p.s]][p.sy + my[p.s]] = H(p.sx + mx[p.s], p.sy + my[p.s]);
                    add(Node(p.s ^ 1, p.sx + mx[p.s], p.sy + my[p.s]));
                } else {
                    g[p.s ^ 1][p.sx + mx[p.s]][p.sy + my[p.s]] = G;
                    up(rk[p.s ^ 1][p.sx + mx[p.s]][p.sy + my[p.s]]);
                }
            }
        }
        if (!flag)cout << -1 << endl;
    }
    return 0;
}

IDA*算法

        接下来就是IDA*算法，这其实是基于迭代加深的A*算法，个人感觉比A*算法简洁。
        IDA*算法类似于DFS，但是它的搜索深度有要求。所谓迭代加深就是每次搜索规定最大搜索深度，如果当前代价与预估代价之和高于最大搜索深度，则进行剪枝。IDA*算法的估价函数性质与A*算法类似。
        看这么一道题。
        答案深度最大为15，很有利于IDA*算法。首先从0到15枚举搜索深度，对于每一个深度进行搜索，搜索代码如下：

void IDA_star(int x, int y, int d) {
    if (d == maxDepth) {
        if (getH() == 0)ans = d;//找到答案，由于是从小到大枚举深度，故答案一定是在最大深度处取到的
        return;
    }
    if (ans != -1)return;
    for (int i = 0; i < 8; i++) {//枚举8个跳转方向
        int xx = x + mx[i], yy = y + my[i];
        if (OK(xx, yy)) {
            swap(op[x][y], op[xx][yy]);//交换
            if (d + getH() <= maxDepth) IDA_star(xx, yy, d + 1);//只要f(n)=g(n)+h(n)的值不大于最大深度，就进行转移
            swap(op[x][y], op[xx][yy]);
        }
    }
}

        只要可以保证h(n)始终不大于实际最小代价，就可以得到最优解。但是h(n)如何选取？
        在这个问题上，规定h(n)为当前状态与目标状态不同的方格数目，可以证明这种h(n)的选取方式满足上述条件。简单证明如下：
        观察到代码中的表达式d + getH()，其实这里的getH()计算的是已经进行了一次转移（就是上一行的swap）后的预估代价。如果实际最小代价为d(n)，能够证明getH()≤d(n)+1即可。
        getH()表示当前有多少方格与目标格不同，假设有p个。由于方格的交换必定是连续的（每一次都有空白方格的参与），故要想把p个方格的位置都移动到正确的位置，至少需要移动p-1次，那么有d(n)≥p-1，即d(n)+1≥p=h(n)，证毕。
        这样直接进行IDA*算法即可。

#include<bits/stdc++.h>

#define OK(x, y) (x>=1&&x<=5&&y>=1&&y<=5)
using namespace std;
const int ANS[6][6] = {
        {0, 0, 0, 0, 0, 0},
        {0, 1, 1, 1, 1, 1},
        {0, 0, 1, 1, 1, 1},
        {0, 0, 0, 2, 1, 1},
        {0, 0, 0, 0, 0, 1},
        {0, 0, 0, 0, 0, 0}
}, mx[] = {-1, -2, -2, -1, 1, 2, 2, 1}, my[] = {-2, -1, 1, 2, 2, 1, -1, -2};
int op[10][10], maxDepth, ans;

inline int getH() {
    int h = 0;
    for (int i = 1; i <= 5; i++) {
        for (int j = 1; j <= 5; j++)if (op[i][j] != ANS[i][j])h++;
    }
    return h;
}

void IDA_star(int x, int y, int d) {
    if (d == maxDepth) {
        if (getH() == 0)ans = d;
        return;
    }
    if (ans != -1)return;
    for (int i = 0; i < 8; i++) {
        int xx = x + mx[i], yy = y + my[i];
        if (OK(xx, yy)) {
            swap(op[x][y], op[xx][yy]);
            if (d + getH() <= maxDepth) IDA_star(xx, yy, d + 1);
            swap(op[x][y], op[xx][yy]);
        }
    }
}

int main() {
    ios::sync_with_stdio(false);
    int q, x, y;
    cin >> q;
    for (int z = 1; z <= q; z++) {
        char e;
        for (int i = 1; i <= 5; i++) {
            for (int j = 1; j <= 5; j++) {
                cin >> e;
                if (e == '0' || e == '1')op[i][j] = e - '0';
                else op[i][j] = 2, x = i, y = j;
            }
        }
        for (int i = 0; i < 16; i++) {
            maxDepth = i, ans = -1, IDA_star(x, y, 0);
            if (ans != -1) {
                cout << ans << endl;
                break;
            }
        }
        if (ans == -1)cout << -1 << endl;
    }
    return 0;
}

我们下面尝试用IDA*算法解决八数码问题。
设置最大深度1000（自行预估），进行IDA*算法即可。注意这里需要一处最优性剪枝：下一个状态的转移方向不能是上一个状态转移方向的反方向。举例来说，如果通过向下移动得到了一个状态，那么这个状态就不能再向上移动，因为这一定不是最优解。如果移动方向设置巧妙的话，可以用异或运算迅速求出反方向。

#include<bits/stdc++.h>

#define OK(x, y) (x>=1&&x<=3&&y>=1&&y<=3)
using namespace std;
int mp[5][5], depth, ans = -1;
const int mx[4] = {-1, 1, 0, 0}, my[4] = {0, 0, 1, -1};
const int ANS[4][4] = {
        {0, 0, 0, 0},
        {0, 1, 2, 3},
        {0, 8, 0, 4},
        {0, 7, 6, 5}
};

inline int getH() {
    int h = 0;
    for (int i = 1; i <= 3; i++)
        for (int j = 1; j <= 3; j++)if (mp[i][j] != ANS[i][j])h++;
    return h;
}

void IDA_star(int x, int y, int d, int pre) {//pre是上一个的移动方向
    if (d == depth) {
        if (getH() == 0)ans = d;
        return;
    }
    if (ans != -1)return;
    for (int i = 0; i < 4; i++) {
        int xx = x + mx[i], yy = y + my[i];
        if (OK(xx, yy) && pre != (i ^ 1)) {
            swap(mp[x][y], mp[xx][yy]);
            if (getH() + d <= depth)IDA_star(xx, yy, d + 1, i);
            swap(mp[x][y], mp[xx][yy]);
        }
    }
}

int main() {
    int x, y;
    for (int i = 0; i < 9; i++) {
        char e;
        cin >> e;
        mp[i / 3 + 1][i % 3 + 1] = e - '0';
        if (mp[i / 3 + 1][i % 3 + 1] == 0)x = i / 3 + 1, y = i % 3 + 1;
    }
    for (int i = 0; i <= 1000; i++) {
        ans = -1, depth = i, IDA_star(x, y, 0, -1);
        if (ans != -1) {
            cout << ans;
            return 0;
        }
    }
}

A*算法应用：图的第k短路

        A*算法的一个应用是求给定图的第k短路。
        求最短路是很容易的，用Dijkstra或者SPFA就能求出来，但第k短却不是那么容易。
        根据上面的讨论可以发现，如果h(n)≤d(n)，A*一定可以找到最有解（即最短路），当找到路时结束算法。如果不结束呢？如果对于出优先队列状态我们不将其加入结束列表，而是将其恢复成未进入过队列的状态以做到点可以重复入队，那这样找到最短路后再找的就是次短路，进而求出第k短路。当目标状态出队列达k次时，得到的就是第k短路。
        估价函数如何规定呢？为了找到最优解，必须满足h(n)≤d(n)，显然h(n)取这个点到目标点的最小距离即可。这个最小距离可以通过建反向图跑一遍Dijkstra或SPFA求出来。
        另外，边权非负的无向图一定有第k短路，有向图可能没有，这时根据队列已空但目标状态出队次数未到k这个指标来判断。
        模板题：POJ2446，注意坑点：s=t时k要自增1，因为0不算最短路。

//常年手写堆的我终于用优先队列了，因为这题手写堆过于复杂
#include<iostream>
#include<queue>
#include<cstring>

using namespace std;
struct Edge {
    int to, next, v;
} edge[100005], edge_opp[100005];

struct DNode {
    int t, dis;

    DNode(int x, int y) : t(x), dis(y) {}

    bool operator<(DNode x) const {
        return dis > x.dis;
    }
};

struct ANode {
    int t, g, h;

    ANode(int x, int y, int z) : t(x), g(y), h(z) {}

    bool operator<(ANode x) const {
        return g + h > x.g + x.h;
    }
};

int head[1005], head_opp[1005], cnt = 1, cnt_opp = 1, dict[1005], n, m, s, t, k, times[1005];


inline void add(int x, int y, int z) {
    edge[cnt].to = y, edge[cnt].next = head[x], edge[cnt].v = z, head[x] = cnt++;
}

inline void add_opp(int x, int y, int z) {
    edge_opp[cnt_opp].to = y, edge_opp[cnt_opp].next = head_opp[x], edge_opp[cnt_opp].v = z, head_opp[x] = cnt_opp++;
}

inline void Dijkstra() {
    priority_queue<DNode> Q;
    memset(dict, 127, sizeof(dict)), dict[t] = 0, Q.push(DNode(t, 0));
    while (!Q.empty()) {
        DNode t = Q.top();
        Q.pop();
        if (t.dis > dict[t.t])continue;
        for (int i = head_opp[t.t]; i; i = edge_opp[i].next) {
            if (dict[edge_opp[i].to] > dict[t.t] + edge_opp[i].v) {
                dict[edge_opp[i].to] = dict[t.t] + edge_opp[i].v;
                Q.push(DNode(edge_opp[i].to, dict[edge_opp[i].to]));
            }
        }
    }
}

inline int AStar(int a, int b) {
    priority_queue<ANode> Q;
    Q.push(ANode(a, 0, dict[a]));
    while (!Q.empty()) {
        ANode t = Q.top();
        Q.pop(), times[t.t]++;
        if (t.t == b && times[t.t] == k)return t.g;
        for (int i = head[t.t]; i; i = edge[i].next) {
            Q.push(ANode(edge[i].to, t.g + edge[i].v, dict[edge[i].to]));
        }
    }
    return -1;
}

int main() {
    ios::sync_with_stdio(false);
    cin >> n >> m;
    for (int i = 1; i <= m; i++) {
        int x, y, z;
        cin >> x >> y >> z;
        add(x, y, z), add_opp(y, x, z);
    }
    cin >> s >> t >> k;
    if (s == t)k++;
    Dijkstra(), cout << AStar(s, t);
    return 0;
}

更新：这里有一个优化。当求k短路时，若某一个点访问次数已经大于k，则可以直接continue，这是因为通过该点得到的路径一定不会排名在前k。这个优化可以减小时空消耗。