编程农场——入坑游玩第二天

合集 - 编程农场(2)

1.编程农场——入坑游玩第二天09-14

2.编程农场——入坑游玩第一天09-13

收起

书接上回

显然，我们可以知道，如果仅初始化图为线性规划的图，显然在模拟退火随机扰动时，会使其不同种类计数扰动，而这一扰动会导致后面完成退火后，得到的计数并不同初始图一样保证最优，我们选择第二种方案，我们可以将cnt1、cnt2、cnt3设为固定的，扰动的只有分布，但是我太菜了，没有想到这一扰动应该怎么实现，于是我又问了ai，他给出了一个十分巧妙的方法，我们不需要考虑如何在一个图里面维护三个种类的作物数量，而是先将初始化图设为满足最优数量的图，然后我们只需要随机交换里面的格子，就可以保证最优数量不会改变。

实际上在评估的过程中，我发现了一些反我的直觉的事情，我在通过线性规划预处理后的模拟退火分数甚至不如纯粹的模拟退火，不过经过一段时间的思考过后，我发现这个现象是显然存在的，因为线性规划并没有考虑到时序的问题，因此，线性规划的最优解在我目前算法设计中，可以退化成每个格子都有一台无人机，每当成熟立刻完成所有流程，本来我是想着，那么这个方案不行，算了吧，不过ai给了我一个新的方案，用于更好的提高SA质量，他重新设计了一个稳态函数代替原本线性规划的预处理，他假设了，每当我们遍历一次全图后，前面的都成熟了，所以加入了移动时间，但是依旧不考虑作物的成熟，这么做，让设定的最优数量更准确了。

第五版代码-结合一下

//Hard-version
#include 
#define ll long long
using namespace std;

// ===== 随机数 =====
// 用 mt19937 取代 rand()：MSVC 的 RAND_MAX 只有 32767，一轮 SA 要消耗约 37 万次
// 随机数，LCG 的周期和低比特质量都不够用。
// 默认用高精度时钟播种（每次运行结果不同）；需要复现时设环境变量 SA_SEED=<整数>。
static mt19937 rng;
inline int rnd_int(int lo, int hi){ return uniform_int_distribution<int>(lo, hi)(rng); }
inline double rnd01(){ return uniform_real_distribution<double>(0.0, 1.0)(rng); }

const double delta = 0.9112;
const int MAXN = 2002;        // 数组宽度（下标 1..n 有效，留出 [MAXN][MAXN] 的余地）
const int ROUNDS = 10;        // cal 内部蛇形遍历的轮数，越大越接近稳态、耗时也线性增加
int n,t;
int v;
int t_move,t_harvest,t_plant,t_till; //基础耗时（输入原值）
// 折算后的实际耗时：题目规定每个动作 Δt = ceil(基础耗时 / v)
int e_move,e_harvest,e_plant,e_till;
int g_grass,g_bush,g_corrot; //成熟耗时
int num_hay,num_wood,num_corrot; //作物数量
int perharvest_grass,perharvest_bush,perharvest_corrot; //单次收获量
int nowx,nowy; //当前无人机位置
int mp[MAXN][MAXN];           //当前解
int tmp_mp[MAXN][MAXN];       //候选解
int best_mp[MAXN][MAXN];      //历史最优解
ll time_mp[MAXN][MAXN];       //cal 用：各格上次被操作时间（仅 [1..n] 范围）
int vis_tag[MAXN][MAXN];      //扰动去重戳记
int cur_tag = 0;
int k0;
double ans = -1e9;            // 当前解 mp 的适应度，每轮 SA 开头会被校准为 cal(mp)
double best = -1e18;
unsigned int g_seed = 0;      // 本轮实际使用的随机种子，写进 meta.txt 便于复现
int cnt1,cnt2,cnt3;

//inline void solve(){ //线性规划找cnt
//	int t1 = e_harvest+g_grass; //草全流程耗时
//	int t2 = e_harvest+e_plant+g_bush; //灌木全流程耗时
//	int t3 = e_harvest+e_plant+e_till+g_corrot; //胡萝卜全流程耗时
//	
//	double a1 = 1.0*perharvest_grass/t1;
//	double a2 = 1.0*perharvest_bush/t2;
//	double a3 = 1.0*perharvest_corrot/t3;
//	double k = 1.0/t3;
//	
//	int num = pow(n,2);
//	double best_sum = -1e9;   // 目标值 = (u1*u2*u3)^(1/3)
//	for (int c3=0;c3<=num;c3++){
//		int R = num-c3;          // 分给草/灌木的格子
//		double K = k*c3;         // 胡萝卜对草、灌木的消耗速率
//		double u3 = c3*a3;
//		if (u3 <= 0 or a1 <= 0 or a2 <= 0) continue; // 几何平均要求三项都为正
//		
//		// u1>0 且 u2>0 对应的 c1 整数可行区间
//		int lo = (int)floor(K/a1) + 1;       // c1*a1 > K
//		int hi = R - (int)floor(K/a2) - 1;   // (R-c1)*a2 > K
//		if (lo < 0) lo = 0;
//		if (hi > R) hi = R;
//		if (lo > hi) continue;
//		
//		// 固定 c3 时 u3 是常数，只需最大化 u1*u2。
//		// u1*u2 = (c1*a1-K)*((R-c1)*a2-K) 是 c1 的开口向下抛物线，
//		// 顶点 x* = R/2 + K*(a2-a1)/(2*a1*a2)，整数最优必在 floor/ceil 处。
//		double xstar = R/2.0 + K*(a2-a1)/(2.0*a1*a2);
//		int cand[2] = {(int)floor(xstar), (int)ceil(xstar)};
//		for (int q=0;q<2;q++){
//			int c1 = cand[q];
//			if (c1 < lo) c1 = lo;
//			if (c1 > hi) c1 = hi;
//			int c2 = R-c1;
//			double u1 = c1*a1 - K;
//			double u2 = c2*a2 - K;
//			if (u1 <= 0 or u2 <= 0) continue;
//			double s = cbrt(u1*u2*u3);   // 几何平均
//			if (s > best_sum){
//				best_sum = s;
//				cnt1 = c1;
//				cnt2 = c2;
//				cnt3 = c3;
//			}
//		}
//	}
//	if (best_sum <= -1e8){
//		cerr << "警告：几何平均下无可行解（格子太少，放不下三种作物且都为正产出）" << endl;
//	}
//	cerr << "geometric mean = " << best_sum
//	     << "   (grass/bush/carrot = " << cnt1 << "/" << cnt2 << "/" << cnt3 << ")" << endl;
//}

inline void solve(){
    int n2 = n*n;
    // 每格每轮固定成本
    int t1 = e_harvest;                     // 草的收获开销
    int t2 = e_harvest + e_plant;           // 灌木
    int t3 = e_harvest + e_plant + e_till;  // 胡萝卜
    int move_cost = n2 * e_move;

    double best_sum = -1e9;
    cnt1 = cnt2 = cnt3 = 0;

    for (int c3 = 0; c3 <= n2; c3++){
        int R = n2 - c3;                    // 分给草/灌木
        if (R < 2) continue;

        // T = move_cost + c1*t1 + (R-c1)*t2 + c3*t3 = K0 - c1*(t2-t1)
        // 其中 K0 = move_cost + R*t2 + c3*t3
        long long K0 = (long long)move_cost + 1LL*R*t2 + 1LL*c3*t3;
        int dt = t2 - t1;                   // 一般 dt >= 0

        // 目标：maximize (c1-c3) * (R-c1-c3) * c3 / T^3
        // 先扫一遍 c1 的整数范围，因为 T 依赖 c1，抛物线顶点不是 O(1)
        // 范围 [c3+1, R-c3-1]（保证 u1,u2 > 0）
        int lo = c3 + 1;
        int hi = R - c3 - 1;
        if (lo > hi) continue;

        // 一维搜索：T 依赖 c1，但只线性依赖，抛物线顶点有显式解
        // d/dc1 [ ln(c1-c3) + ln(R-c1-c3) - 3 ln(K0 - c1*dt) ] = 0
        // 1/(c1-c3) - 1/(R-c1-c3) + 3*dt/(K0 - c1*dt) = 0
        // 用牛顿或二分解，这里为了简洁直接三分扫一遍区间（区间不大，n=30 时最多 900 点）
        double best_local = -1e9;
        int best_c1 = -1;
        for (int c1 = lo; c1 <= hi; c1++){
            int c2 = R - c1;
            long long T = K0 - 1LL*c1*dt;
            if (T <= 0) continue;
            double u1 = 1.0*(c1 - c3)/T;
            double u2 = 1.0*(c2 - c3)/T;
            double u3 = 1.0*c3/T;
            if (u1 <= 0 or u2 <= 0 or u3 <= 0) continue;
            double s = cbrt(u1*u2*u3);
            if (s > best_local){
                best_local = s;
                best_c1 = c1;
            }
        }
        if (best_c1 < 0) continue;
        if (best_local > best_sum){
            best_sum = best_local;
            cnt1 = best_c1;
            cnt2 = R - best_c1;
            cnt3 = c3;
        }
    }

    if (best_sum <= -1e8){
        cerr << "警告：稳态模型下无可行解" << endl;
    }
    cerr << "steady-state mean = " << best_sum
         << "   (grass/bush/carrot = " << cnt1 << "/" << cnt2 << "/" << cnt3 << ")" << endl;
}

inline void init(){  //初始化
	vector<int> vals;
	vals.reserve((size_t)n*n);
	for (int i=0; iint)vals.size()push_back(1);
	for (int i=0; iint)vals.size()push_back(2);
	for (int i=0; iint)vals.size()push_back(3);
	while ((int)vals.size() < (size_t)n*n) vals.push_back(1);   // 兜底
	shuffle(vals.begin(), vals.end(), rng);
	int idx = 0;
	for (int i=1;i<=n;i++)
		for (int j=1;j<=n;j++)
			mp[i][j] = vals[idx++];
}

// Δt = ceil(基础耗时 / v)，v 至少按 1 处理
inline int eff(int base){
	int vv = (v < 1) ? 1 : v;
	return (base + vv - 1) / vv;
}

// 计算布局 a 的适应度（纯函数：只读 a，不依赖任何全局状态）
// 无人机按蛇形顺序遍历 ROUNDS 轮，逐格判断成熟并执行操作，统计净产出率
inline double cal(int a[MAXN][MAXN]){
	// 只清理 [1..n] 范围的 time_mp，避免对 32MB 整块 memset
	for (int i=1;i<=n;i++){
		memset(time_mp[i]+1, 0, n*sizeof(ll));
	}
	int cnt1 = 0,cnt2 = 0,cnt3 = 0; //cnt1: hay cnt2: wood cnt3: corrot
	ll global_time = 0; //全局时间
	// 单次收获操作的耗时（已按 Δt=ceil(基础/v) 折算；成长时间由下方的"delta >= 成长阈值"判断隐式处理）
	int t1 = e_harvest;                          // 草：收割后原地重生，无需补种
	int t2 = e_harvest+e_plant;                  // 灌木：收割 + 补种
	int t3 = e_harvest+e_plant+e_till;           // 胡萝卜：收割 + 翻地 + 补种
	for (int k=1;k<=ROUNDS;k++){
		for (int i=1;i<=n;i++){
			// 蛇形扫描：奇数行正序、偶数行逆序
			int js, jt, stp;
			if (i%2 == 0){ js = n; jt = 0;  stp = -1; }
			else         { js = 1; jt = n+1; stp =  1; }
			for (int j=js; j!=jt; j+=stp){
				int type = a[i][j];
				ll delta_time = global_time - time_mp[i][j];
				if (type == 1){
					if (delta_time >= g_grass){
						cnt1++;
						global_time += t1;
						time_mp[i][j] = global_time;
					}
				}
				else if (type == 2){
					if (delta_time >= g_bush){
						cnt2++;
						global_time += t2;
						time_mp[i][j] = global_time;
					}
				}
				else if (type == 3){
					if (delta_time >= g_corrot){
						cnt3++;
						global_time += t3;
						time_mp[i][j] = global_time;
					}
				}
				global_time += e_move;
			}
		}
	}
	ll num1 = 1LL*perharvest_grass*cnt1; //草收获量
	ll num2 = 1LL*perharvest_bush*cnt2; //灌木收获量
	ll num3 = 1LL*perharvest_corrot*cnt3; //胡萝卜收获量
	// 每棵胡萝卜种植消耗 1 干草 + 1 木材，故从草/灌木的净产出里各扣除 cnt3
	double u1 = 1.0*num1/global_time - 1.0*cnt3/global_time; //单位时间草净产出
	double u2 = 1.0*num2/global_time - 1.0*cnt3/global_time; //单位时间灌木净产出
	double u3 = 1.0*num3/global_time; //单位时间胡萝卜产出

	if (u1 < 0 or u2 < 0){
		return -1e18; // 入不敷出：草/灌木产量撑不起胡萝卜消耗，方案不可行
	}

//	double sum = 3/(1/u1+1/u2+1/u3);
	double sum = cbrt(u1*u2*u3);
	return sum;
}

inline void SA(){
	k0 = n*n;
	// 把 ans 校准成当前解的真实适应度，避免历史残留造成 best_mp 与分数错位
	ans = cal(mp);

	// 自适应初始温度：采样若干单格扰动，让典型劣化移动在 T0 时的接受率约 0.5
	double T0;
	{
		const int SAMPLES = 30;
		double s = 0;
		for (int q=0; qint i1 = rnd_int(1, n), j1 = rnd_int(1, n);
			int i2 = rnd_int(1, n), j2 = rnd_int(1, n);
			int tries = 0;
			while ((i1 == i2 and j1 == j2)or(mp[i1][j1] == mp[i2][j2])){
				i1 = rnd_int(1,n);
				j1 = rnd_int(1,n);
				i2 = rnd_int(1,n);
				j2 = rnd_int(1,n);
				if (++tries > 1000){
					break;
				}
			}
			if (++tries > 1000) continue;
			swap(mp[i1][j1], mp[i2][j2]);
			double v = cal(mp);
			swap(mp[i1][j1],mp[i2][j2]);
			s += fabs(v - ans);
		}
		T0 = max(s/SAMPLES, 1e-6) / log(2.0);
		if (T0 < 1e-4) T0 = 1e-4;
	}
	double T = T0;
	double T_min = max(1e-9, T0 * 1e-6);
	while (T > T_min){
		int k = max(1,(int)(k0*(T/T0)));
		vectorint,4>> swaps;
		cur_tag++; // 本温度步的扰动去重戳记（等价于原来的 memset(vis,0,...)）
		for (int i=0;iint tries = 0;
			while (true){
				int i1 = rnd_int(1,n),j1 = rnd_int(1,n);
				int i2 = rnd_int(1,n),j2 = rnd_int(1,n);
				if (i1 == i2 and j1 == j2){
					if (++tries > 200) break;
					continue;
				}
				if (vis_tag[i1][j1] == cur_tag or vis_tag[i2][j2] == cur_tag){
					if (++tries > 200){
						break;
					}
					continue;
				}
				if (mp[i1][j1] == mp[i2][j2]){
					if (++tries > 200){
						break;
					}
					continue;
				}
				vis_tag[i1][j1] = vis_tag[i2][j2] = cur_tag;
				swaps.push_back({i1,j1,i2,j2});
				break;
			}
		}
		// 拷贝当前解到候选解
		for (int i=1;i<=n;i++){
			memcpy(tmp_mp[i]+1, mp[i]+1, n*sizeof(int));
		}
		for (auto [i1,j1,i2,j2] : swaps){
			swap(tmp_mp[i1][j1],tmp_mp[i2][j2]);
		}
		double now = cal(tmp_mp); // 评估候选解（关键修复点：原来是 cal(mp)，等于没评估扰动）
		double Delta = now-ans;
		if (Delta >= 0 or exp(Delta/T) > rnd01()){
			// 接受：覆盖 mp，并把 ans 同步到新 mp 的真实分数
			for (int i=1;i<=n;i++){
				memcpy(mp[i]+1, tmp_mp[i]+1, n*sizeof(int));
			}
			ans = now;
			if (ans > best){
				best = ans;
				for (int i=1;i<=n;i++){
					memcpy(best_mp[i]+1, mp[i]+1, n*sizeof(int));
				}
			}
		}
		T *= delta;
	}
}

signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	cin>>n>>t;
	cin>>v;
	cin>>t_move>>t_harvest>>t_plant>>t_till;
	cin>>g_grass>>g_bush>>g_corrot;
	cin>>num_hay>>num_wood>>num_corrot;
	cin>>perharvest_grass>>perharvest_bush>>perharvest_corrot;
	cin>>nowx>>nowy;
	// 按 Δt = ceil(基础耗时/v) 折算实际耗时（v 是速度因子，见题目）
	e_move    = eff(t_move);
	e_harvest = eff(t_harvest);
	e_plant   = eff(t_plant);
	e_till    = eff(t_till);
	// 播种：默认高精度时钟（每次运行结果不同）；想复现某次结果就设 SA_SEED=<整数>
	g_seed = (unsigned)chrono::high_resolution_clock::now().time_since_epoch().count();
	if (const char* env = getenv("SA_SEED")) g_seed = (unsigned)strtoul(env, nullptr, 10);
	rng.seed(g_seed);
	cerr << "seed = " << g_seed << endl;
	//3^(n*n) NP-hard?
	solve();
	init();
	int times = 100;
	while (times--){
		SA();
	}
	// 还原历史最优到 mp 并写出
	for (int i=1;i<=n;i++){
		memcpy(mp[i]+1, best_mp[i]+1, n*sizeof(int));
	}
	for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            cout << mp[i][j];
            if (j < n) cout << ',';
        }
        cout << '\n';
    }
	ofstream fout("layout.csv");
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            fout << mp[i][j];
            if (j < n) fout << ',';
        }
        fout << '\n';
    }
    fout.close();

    // ========== 输出参数到 meta.txt ==========
    ofstream meta("meta.txt");
    meta << "n=" << n << "\n";
    meta << "best_score=" << best << "\n";
    meta << "seed=" << g_seed << "\n";
    meta << "t_move=" << t_move << "\n";
    meta << "t_harvest=" << t_harvest << "\n";
    meta << "t_plant=" << t_plant << "\n";
    meta << "t_till=" << t_till << "\n";
    meta << "v=" << v << "\n";
    meta << "e_move=" << e_move << "\n";
    meta << "e_harvest=" << e_harvest << "\n";
    meta << "e_plant=" << e_plant << "\n";
    meta << "e_till=" << e_till << "\n";
    meta << "g_grass=" << g_grass << "\n";
    meta << "g_bush=" << g_bush << "\n";
    meta << "g_corrot=" << g_corrot << "\n";
    meta << "perharvest_grass=" << perharvest_grass << "\n";
    meta << "perharvest_bush=" << perharvest_bush << "\n";
    meta << "perharvest_corrot=" << perharvest_corrot << "\n";
    meta.close();

    // ========== 调用 Python 画图 ==========
#ifdef _WIN32
    int ret = system("python plot_layout.py");
#else
    int ret = system("python3 plot_layout.py");
#endif
    if (ret != 0) {
        cerr << "警告：Python 脚本执行失败，请检查 Python 环境和 plot_layout.py" << endl;
    }
	return 0;
}

到此为止我想应该没有什么可以优化的了。下一步就去解锁更多好东西吧。