update

Author: sunkafei

算法/二分搜索/二分查找.cpp

@@ -0,0 +1,47 @@
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <iostream>
+#include <algorithm>
+#include <functional>
+#include <array>
+using namespace std;
+const int dim = 2;
+using vec = array<double, dim>;
+/*
+* grad: 计算梯度的函数
+* pos: 搜索的起始点
+* lr: 学习率（移动的步长）
+* decay: 学习率的衰减指数
+* exit: 当学习率小于exit时算法结束
+* 注意该算法求出的是 函数的极**小**值点
+*/
+vec Adam(function<vec(vec)> grad, vec pos, double lr, double decay, double exit,
+	const double beta1 = 0.9, const double beta2 = 0.9) {
+	double pw1 = 1, pw2 = 1;
+	vec m = {}, v = {};
+	for (int t = 1; lr >= exit; ++t) {
+		vec g = grad(pos);
+		pw1 *= beta1;
+		pw2 *= beta2;
+		for (int i = 0; i < dim; ++i) {
+			m[i] = m[i] * beta1 + g[i] * (1 - beta1);
+			v[i] = v[i] * beta2 + g[i] * g[i] * (1 - beta2);
+			double update = m[i] / (1 - pw1);
+			double factor = lr / (sqrt(v[i] / (1 - pw2)) + 1e-8);
+			pos[i] -= update * factor;
+		}
+		lr *= decay;
+	}
+	return pos;
+}
+
+int main() {
+	auto grad = [](vec A) {
+		double x = A[0], y = A[1];
+		return vec{ 2 * y + -4 / (x * x), 2 * x + -4 / (y * y) };
+	};
+	auto ans = Adam(grad, vec{ 1, 1 }, 0.01, 0.999, 1e-8);
+	printf("%.10f %.10f\n", ans[0], ans[1]); //1.2599210498948732
+	return 0;
+}

算法/人工智能/Adam算法/Adam算法.cpp

@@ -0,0 +1,120 @@
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <iostream>
+#include <algorithm>
+#include <random>
+#include <tuple>
+#include <vector>
+#include <map>
+using namespace std;
+const int nrow = 10, ncol = 10, max_state = nrow * ncol, max_action = 4;
+const int dr[4] = { -1, 1, 0, 0 };
+const int dc[4] = { 0, 0, -1, 1 };
+struct CliffWalkingEnv {
+	char s[11][11];
+	int x, y;
+	vector<pair<int, int>> path;
+	CliffWalkingEnv() {
+		reset();
+	}
+	int reset() {
+		x = 0, y = 0;
+		path.clear();
+		path.emplace_back(x, y);
+		memset(s, 0, sizeof(s));
+		for (int i = 0; i < nrow; ++i) {
+			for (int j = 0; j < ncol; ++j) {
+				s[i][j] = '.';
+			}
+		}
+		for (int i = 0; i < 5; ++i) {
+			for (int j = 1; j < ncol - 1; ++j) {
+				s[i][j] = '#';
+			}
+		}
+		return x * ncol + y;
+	}
+	void print() {
+		char t[11][11];
+		memcpy(t, s, sizeof(s));
+		for (auto [x, y] : path) {
+			if (t[x][y] >= '1' && t[x][y] < '9')
+				t[x][y] += 1;
+			else
+				t[x][y] = '1';
+		}
+		for (int i = 0; i < nrow; ++i)
+			printf("%s\n", t[i]);
+		printf("\n");
+	}
+	tuple<int, int, int> step(int number) {
+		int nx = x + dr[number], ny = y + dc[number];
+		if (nx >= 0 && nx < nrow && ny >= 0 && ny < ncol)
+			x = nx, y = ny;
+		int reward = -1, done = false;
+		if (s[x][y] == '#') {
+			reward = -100;
+			done = true;
+		}
+		if (x == 0 && y == ncol - 1)
+			done = true;
+		path.emplace_back(x, y);
+		return make_tuple(x * ncol + y, reward, done);
+	}
+} env;
+
+const double alpha = 1e-1; /*学习率*/
+const double gamma = 0.9; /*折扣因子*/
+const double epsilon = 1e-2; /*epsilon-greedy*/
+const int N = 10; //Q-Planning的次数
+const int max_episode = 200; //训练多少轮
+uniform_real_distribution<double> p;
+uniform_int_distribution<int> d;
+default_random_engine e;
+double Q[max_state][max_action];
+map<pair<int, int>, int> id;
+vector<tuple<int, int, int, int>> model;
+int epsilon_greedy(int state) { /*基于epsilon-greedy选择动作*/
+	if (p(e) < epsilon)
+		return d(e) % max_action;
+	return max_element(Q[state], Q[state] + max_action) - Q[state];
+}
+void learn(int s0, int a0, int r, int s1) {
+	auto td_error = r + gamma * *max_element(Q[s1], Q[s1] + max_action) - Q[s0][a0];
+	Q[s0][a0] += alpha * td_error;
+}
+void update(int s0, int a0, int r, int s1) { 
+	learn(s0, a0, r, s1);
+	pair<int, int> pr(s0, a0);
+	if (!id.count(pr)) {
+		id[pr] = model.size();
+		model.emplace_back(s0, a0, r, s1);
+	}
+	else {
+		model[id[pr]] = make_tuple(s0, a0, r, s1);
+	}
+	for (int i = 0; i < N; ++i) {
+		int idx = d(e) % model.size();
+		auto [s0, a0, r, s1] = model[idx];
+		learn(s0, a0, r, s1);
+	}
+}
+void DynaQ() {
+	for (int i = 0; i < max_episode; ++i) {
+		int state = env.reset();
+		for (;;) {
+			int action = epsilon_greedy(state);
+			auto [next_state, reward, done] = env.step(action);
+			update(state, action, reward, next_state);
+			state = next_state;
+			if (done)
+				break;
+		}
+		env.print();
+	}
+}
+int main() {
+	DynaQ();
+	return 0;
+}

算法/人工智能/Q-learning算法/Dyna-Q.cpp

@@ -0,0 +1,126 @@
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <iostream>
+#include <algorithm>
+#include <vector>
+#include <array>
+#include <cmath>
+#include <ctime>
+using namespace std;
+const int maxn = 210000;
+const int maxdim = 1001;
+const double eps = 1e-8;
+using matrix = double[maxdim][maxdim];
+using vec = array<double, maxdim>;
+using pair_t = pair<double, vec>;
+struct PCA {
+	matrix A, V;
+	int column[maxdim], n;
+	void update(int r, int c, double v) {
+		A[r][c] = v;
+		if (column[r] == c || fabs(A[r][c]) > fabs(A[r][column[r]])) {
+			for (int i = 0; i < n; ++i) if (i != r)
+				if (fabs(A[r][i]) > fabs(A[r][column[r]]))
+					column[r] = i;
+		}
+	}
+	void Jacobi() {
+		for (int i = 0; i < n; ++i) {
+			for (int j = 0; j < n; ++j)
+				V[i][j] = 0;
+			V[i][i] = 1;
+			column[i] = (i == 0 ? 1 : 0);
+		}
+		for (int i = 0; i < n; ++i) 
+			for (int j = 0; j < n; ++j) 
+				if (j != i && fabs(A[i][j]) > fabs(A[i][column[i]]))
+					column[i] = j;
+		for (int T = 0; ; ++T) { //迭代次数限制
+			int x, y;
+			double val = 0;
+			for (int i = 0; i < n; ++i)
+				if (fabs(A[i][column[i]]) > val)
+					val = fabs(A[i][column[i]]), x = i, y = column[i];
+			if (val < eps) //精度限制
+				break;
+			double phi = atan2(-2 * A[x][y], A[y][y] - A[x][x]) / 2;
+			double sinp = sin(phi), cosp = cos(phi);
+			for (int i = 0; i < n; ++i) if (i != x && i != y) {
+				double a = A[x][i] * cosp + A[y][i] * sinp;
+				double b = A[x][i] * -sinp + A[y][i] * cosp;
+				update(x, i, a);
+				update(y, i, b);
+			}
+			for (int i = 0; i < n; ++i) if (i != x && i != y) {
+				double a = A[i][x] * cosp + A[i][y] * sinp;
+				double b = A[i][x] * -sinp + A[i][y] * cosp;
+				update(i, x, a);
+				update(i, y, b);
+			}
+			for (int i = 0; i < n; ++i) {
+				double a = V[i][x] * cosp + V[i][y] * sinp;
+				double b = V[i][x] * -sinp + V[i][y] * cosp;
+				V[i][x] = a, V[i][y] = b;
+			}
+			double a = A[x][x] * cosp * cosp + A[y][y] * sinp * sinp + 2 * A[x][y] * cosp * sinp;
+			double b = A[x][x] * sinp * sinp + A[y][y] * cosp * cosp - 2 * A[x][y] * cosp * sinp;
+			double tmp = (A[y][y] - A[x][x]) * sin(2 * phi) / 2 + A[x][y] * cos(2 * phi);
+			update(x, y, tmp);
+			update(y, x, tmp);
+			A[x][x] = a, A[y][y] = b;
+		}
+	}
+	//a 为输入向量组
+	//n 为向量的维数
+	//center 指针用来保存输入向量组的中心点
+	//返回特征值和特征向量的pair，按照特征值从大到小排序
+	//特征值是各个点在对应特征向量方向的坐标平方和，除以(a.size() - 1)为方差。
+	auto solve(vector<vec> a, int n, vec* center = nullptr) {
+		this->n = n;
+		vec s = {};
+		for (int i = 0; i < a.size(); ++i)
+			for (int j = 0; j < n; ++j)
+				s[j] += a[i][j];
+		for (int j = 0; j < n; ++j)
+			s[j] /= a.size();
+		for (int i = 0; i < a.size(); ++i)
+			for (int j = 0; j < n; ++j)
+				a[i][j] -= s[j];
+		if (center) *center = s;
+		for (int i = 0; i < n; ++i) {
+			for (int j = 0; j < n; ++j) {
+				A[i][j] = 0;
+				for (int k = 0; k < a.size(); ++k)
+					A[i][j] += a[k][i] * a[k][j];
+			}
+		}
+		Jacobi();
+		vector<pair_t> result;
+		for (int i = 0; i < n; ++i)
+			result.emplace_back(A[i][i], vec());
+		for (int i = 0; i < n; ++i)
+			for (int j = 0; j < n; ++j)
+				result[i].second[j] = V[j][i];
+		sort(result.begin(), result.end(), greater<pair_t>());
+		return result;
+	}
+}pca;
+int main() { //1329070654.526
+	freopen("in.txt", "r", stdin);
+	vector<vec> a;
+	int n, m;
+	scanf("%d %d", &n, &m);
+	for (int i = 0; i < n; ++i) {
+		vec v = {};
+		for (int j = 0; j < m; ++j)
+			scanf("%lf", &v[j]);
+		a.push_back(v);
+	}
+	auto now = clock();
+	auto result = pca.solve(a, m);
+	for (int i = 0; i < m; ++i)
+		printf("%.3f\n", result[i].first);
+	printf("time: %f\n", double(clock() - now) / CLOCKS_PER_SEC);
+	return 0;
+}