缩减思考过程,尽快输出代码。有一组数据,是工艺参数对应的工艺结果,输入x1对应输出y11、y12、y13、y14,x2对应y21、y22、y23、y24,x3对应y31、y32、y33、y34,x4对应y41、y42、y43、y44,x5对应y51、y52、y53、y54,x6对应y61、y62y63、y64,这些输入输出可能是小数,y是每个x对应的一次实验的4个工艺结果。使用贝叶斯优化中的多目标优化问题的解决方法,用python语言,使用gpytorch找出合适这些数据的二次多项式均值函数和合适的核函数,使用botorch等库实现对这组数据的多目标优化。不要用常数均值函数,因为不符合实际。目的是找到使得工艺结果更好的工艺参数。直接给出完整代码。不要使用botorch.utils.multi_objective库中的FastNondominatedPartitioning函数。
时间: 2025-04-18 13:51:15 浏览: 62
<think>好的,我需要帮用户解决一个多目标贝叶斯优化的问题。首先,用户的数据结构是每个工艺参数x对应四个输出y,每个x有四个实验结果。目标是用贝叶斯优化找到更好的工艺参数,要求使用二次多项式均值函数和合适的核函数,不用常数均值,并且用gpytorch和botorch库,同时避免使用FastNondominatedPartitioning函数。
首先,贝叶斯优化的核心是高斯过程模型。用户要求均值函数是二次多项式,所以我要自定义gpytorch的均值函数。核函数方面,可能用RBF加上其他核的组合,比如尺度核或线性核,但需要根据数据情况调整。
接下来是多目标处理。通常会用Expected Hypervolume Improvement(EHVI)作为采集函数,但用户要求不使用FastNondominatedPartitioning,可能需要手动处理非支配排序或者用其他方法。不过检查botorch文档,发现即使不用那个函数,可能还是有其他方式实现EHVI,或者改用其他采集函数,比如基于分解的方法,比如ParEGO。这里可能需要确认。
数据准备方面,输入X应该是6个参数点,每个点对应四个输出,每个输出有四个结果。但贝叶斯优化通常处理多个目标,这里每个x对应四个输出,可能视为四个目标,需要将数据整理为多维输出。或者用户可能希望将四个实验结果视为同一目标下的多个样本,可能需要取平均或建模分布。这点需要明确。
假设用户的问题是多目标优化,四个输出y1到y4都是需要优化的目标。那么每个x对应的四个实验结果是同一x下四个目标的观测值。或者可能每个x有四个重复实验,每个实验得到四个结果?这里可能存在歧义。需要根据常理推断,可能每个x对应四个指标,每个指标可能有多个实验值。或者用户可能希望将四个实验结果作为四个独立的数据点?需要进一步明确,但用户描述可能每个x对应四个输出,每个输出有四个样本,比如y11到y14是四个样本,可能代表同一工艺参数下的重复实验,导致输出有噪声。但处理多目标时,通常每个x对应一个目标向量,所以可能需要将四个实验的结果取平均,或建模噪声。
假设用户的数据是每个x有四个目标,每个目标有四次实验的结果,那么可能需要将这些结果的平均作为观测值,或者考虑噪声。但在贝叶斯优化中,通常每个x对应一个目标向量,所以可能需要将四个实验的结果取平均,或者将每个实验视为一个单独的数据点,但这样x会被重复四次。这可能不太合理。所以更可能的情况是每个x对应四个目标,每个目标有一个观测值(可能是四次实验的平均),或者每个x的每个目标有四次观测,存在噪声。需要用户澄清,但根据问题描述,可能用户是将每个x的四个输出作为四个独立目标,每个目标只有一个观测值(即四次实验可能对应四个不同的目标?)。或者可能用户的数据结构有误,需要进一步确认。但根据用户给出的例子,x1对应y11到y14,可能这四个是四个不同的目标,每个x对应四个目标值,每个目标值有多次实验的结果?或者每个x的四个结果是对同一目标的四次测量?这里可能存在混淆。
假设用户的问题是多目标优化,四个目标,每个x对应四个目标值,每个目标值是通过四次实验得到的。这时候可能需要将每个目标的四个实验值取平均,作为该x下该目标的观测值,并考虑方差作为噪声。或者,如果每个x的四个实验产生四个目标向量(即每次实验得到四个目标值),那么每个x有四次实验,每次实验得到四个目标值,此时数据可能结构为:输入x(6个不同的x值),每个x有四次实验,每次实验得到四个目标值。这种情况下,每个x有4个样本,每个样本是四个目标值。此时需要将数据整理为X的每个点重复四次,对应的Y是四个目标的值。这可能更合理。
但用户的问题描述可能需要更明确的结构。例如,输入x1,输出是四次实验,每次实验得到四个工艺结果y11到y14,这可能意味着每个x对应四个实验,每个实验得到四个输出。但这样的话,每个x有四个样本,每个样本是四个输出,即每个样本是一个四维向量。此时,X的维度是(6*4, 1),Y的维度是(6*4, 4)。但这样处理的话,模型需要处理多输出问题,比如使用多任务高斯过程,或者为每个目标单独建模。这可能比较复杂。
另一种可能是用户希望将四个工艺结果视为四个不同的目标,每个x对应这四个目标的某个值(可能是四次实验的平均)。这样,每个x对应一个四维的目标向量,然后进行多目标优化。此时,X有6个点,每个点对应四维目标,每个目标的值可能取四次实验的平均值。这种情况下,数据预处理需要将每个x的四个实验的每个目标取平均,得到每个x的四个目标值,然后进行多目标贝叶斯优化。
但用户的问题描述可能存在歧义,需要假设。假设用户的数据结构是:每个x(共6个)对应四个目标(y1到y4),每个目标有四个实验值(例如,y11到y14是目标1在x1下的四次实验结果)。因此,处理方式可能是对每个目标在每个x下的四个实验值取平均,得到每个目标的均值,作为该x下的目标值。这样,X的形状是(6,1),Y的形状是(6,4),每个Y的列对应一个目标的均值。然后进行多目标贝叶斯优化,寻找x使得四个目标尽可能优。
现在,构建高斯过程模型。均值函数是二次多项式,所以需要自定义gpytorch的Mean函数。例如,对于每个目标,均值函数是a*x² + b*x + c。核函数可能需要选择,比如RBF核,加上可能的多项式核或其他组合。但用户没有明确要求,所以需要选择一个合适的组合,比如RBF核。
多目标优化方面,使用botorch,通常使用qNEHVI作为采集函数,但需要非支配排序和超体积计算,而用户要求不使用FastNondominatedPartitioning,可能需要寻找替代方法。或者,可能用户希望避免使用该函数,但botorch内部可能会调用它,所以可能需要寻找其他方式。或者,用户可能希望用其他采集函数,比如基于分解的方法,比如ParEGO,或者使用加权和方法。例如,使用Expected Hypervolume Improvement可能需要FastNondominatedPartitioning,所以如果不能用,可能需要改用其他方法。
或者,用户可能误解了该函数的作用,只是想避免直接调用它,而botorch的更高层API可能封装了这个函数的使用,比如在使用qNoisyExpectedHypervolumeImprovement时,内部可能会用到,但用户可能不希望显式调用。因此,可能需要换用其他采集函数,比如qExpectedHypervolumeImprovement,但可能仍需要非支配排序。或者,用户可能希望自己实现非支配排序部分,但这比较复杂。
或者,用户可能希望将多目标问题转化为单目标问题,例如通过加权和,但这样可能不符合多目标优化的初衷。所以可能需要继续使用EHVI,但避免直接调用FastNondominatedPartitioning。但根据botorch的源码,EHVI可能会依赖该函数,因此可能无法绕过。这可能是一个问题。
如果用户坚持不可以使用该函数,可能需要寻找其他方式,比如使用随机标量化(ParEGO)或基于分解的方法。例如,使用qExpectedHypervolumeImprovement可能必须使用该函数,因此可能需要改用ParEGO。例如,botorch的acquisition.analytic.ParEGOAnalyticAcquisitionFunction。或者,使用ScalarizedPosteriorTransform将多目标转化为单目标,但这可能无法得到帕累托前沿。
因此,可能需要调整策略,改用ParEGO方法,这样不需要FastNondominatedPartitioning函数。这样,可以在代码中使用ScalarizedPosteriorTransform,并构建相应的采集函数。
综上,步骤如下:
1. 数据准备:假设每个x对应四个目标,每个目标取四次实验的平均值,构造X为(6,1),Y为(6,4)。
2. 定义高斯过程模型,使用自定义的二次多项式均值函数,以及合适的核函数(如RBF + Polynomial)。
3. 使用botorch的MultiTaskGP模型来处理多目标(每个任务对应一个目标),或者使用IndependentModelList,但可能更复杂。或者,将每个目标单独建模为一个GP模型,然后组合成多输出模型。
4. 使用ParEGO或其他不依赖FastNondominatedPartitioning的采集函数进行优化。
现在,具体实现:
首先,导入必要的库:
import torch
import gpytorch
from gpytorch.means import Mean
from gpytorch.kernels import ScaleKernel, RBFKernel, PolynomialKernel
from gpytorch.mlls import ExactMarginalLogLikelihood
from botorch.models import SingleTaskGP, MultiTaskGP
from botorch.fit import fit_gpytorch_model
from botorch.acquisition import qExpectedHypervolumeImprovement, AnalyticAcquisitionFunction, PosteriorSampling
from botorch.optim import optimize_acqf
from botorch.utils.sampling import draw_sobol_samples
from botorch.utils.multi_objective import Hypervolume
from botorch.utils.multi_objective.pareto import is_non_dominated
import numpy as np
然后,生成模拟数据。假设用户的数据结构是:
X = [x1, x2, x3, x4, x5, x6],每个xi是一个标量。
每个xi对应四个目标的四次实验,例如,对于x1,四个实验得到四个目标值的四个样本,例如y11到y14是目标1的四次实验,y21到y24是目标2的四次实验,依此类推?或者每个实验得到四个目标值?
这里可能需要重新理解数据结构。用户的问题描述可能更准确的是:每个x对应一次实验的四个工艺结果,而每个x有四次实验,所以每个x有四次实验,每次实验得到四个结果。例如,x1的四次实验得到四个四维向量:y11, y12, y13, y14(每个y1j是一个四维向量)。或者,可能每个x的每个实验产生四个结果,即每个x有四次实验,每次实验得到四个结果,即四个目标。这种情况下,每个x对应四个四维的输出向量。此时,数据应整理为X是重复四次每个x,形成(6*4, 1)的输入,Y是(6*4,4)的输出,每个行是四个目标值。
例如,X = [x1, x1, x1, x1, x2, x2, ..., x6],每个x重复四次,对应的Y是每个实验的四目标值。这种情况下,每个样本是输入x和四个目标,共有24个样本(6个x * 4次实验)。
在这种情况下,可以建立多任务GP模型,每个任务对应一个目标,或者为每个目标单独建立GP模型。这里可能选择独立地为每个目标建立GP模型,然后组合成多目标模型。
但贝叶斯优化中处理多目标通常需要联合模型,但botorch支持多任务GP,或者使用ModelListGP。例如,使用ModelListGP,每个模型对应一个目标。
然后,定义均值函数。需要自定义一个二次多项式均值:
class QuadraticMean(Mean):
def __init__(self, input_size, batch_shape=torch.Size()):
super().__init__()
self.register_parameter(name="weights", parameter=torch.nn.Parameter(torch.randn(*batch_shape, input_size, 1)))
self.register_parameter(name="quadratic", parameter=torch.nn.Parameter(torch.randn(*batch_shape, input_size, 1)))
self.register_parameter(name="bias", parameter=torch.nn.Parameter(torch.randn(*batch_shape, 1)))
def forward(self, x):
return (x**2 @ self.quadratic) + (x @ self.weights) + self.bias
然后,在自定义的GP模型中使用这个均值函数:
class CustomGP(SingleTaskGP):
def __init__(self, train_X, train_Y):
super().__init__(train_X, train_Y)
self.mean_module = QuadraticMean(input_size=1)
self.covar_module = ScaleKernel(RBFKernel() + PolynomialKernel(power=2))
但需要注意输入维度是否正确。例如,input_size=1如果x是单变量。
然后,为每个目标建立一个CustomGP模型,组合成ModelListGP。
在优化过程中,使用ParEGO方法,将多目标转化为标量,进行优化。或者,使用EHVI但避免调用FastNondominatedPartitioning,但可能不可行。
或者,使用Expected Hypervolume Improvement,但需要参考点,而计算Hypervolume可能需要非支配解,但可能无法避免使用相关函数。
用户明确要求不要使用FastNondominatedPartitioning,所以可能需要改用其他方法,例如随机标量化(ParEGO)或基于分解的方法。
在ParEGO中,每次优化时随机生成一个权重向量,将多目标问题转化为单目标问题。这不需要非支配排序,因此可以避免使用FastNondominatedPartitioning函数。
因此,在代码中,可以使用botorch的acquisition.analytic.ParEGOAnalyticAcquisitionFunction,或者手动实现。
或者,使用ScalarizedPosteriorTransform,将多目标转换为单目标,然后使用Expected Improvement。
因此,代码的大致步骤如下:
1. 准备数据:假设X是(24,1),每个x重复四次,Y是(24,4),每次实验的四个目标值。
2. 标准化数据:对每个目标进行标准化,便于模型训练。
3. 为每个目标建立独立的GP模型,使用自定义的均值函数和核函数。
4. 将这些模型组合成ModelListGP。
5. 定义采集函数,例如ParEGO,在每次迭代中随机生成权重,进行优化。
6. 优化采集函数,找到下一个候选点。
7. 更新模型,迭代多次。
但具体实现可能较为复杂,需要处理多目标模型的建立和采集函数的应用。
现在,编写代码的草稿:
首先,生成示例数据。假设有6个x值,每个x重复四次实验,每个实验四个目标。这里,我们创建示例数据:
# 示例数据生成
X = torch.tensor([1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unsqueeze(-1) # 形状 (6, 1)
# 每个x重复四次实验,共24个数据点
X = X.repeat_interleave(4, dim=0) # 形状 (24, 1)
# 生成对应的四个目标的Y值,假设每个实验的四个目标有随机值
Y = torch.randn(24, 4) # 形状 (24, 4)
然后,标准化Y:
from botorch.utils import standardize
Y_std = standardize(Y) # 标准化每个目标列
然后,定义模型:
class QuadraticMeanGP(SingleTaskGP):
def __init__(self, train_X, train_Y):
# 确保输入正确
super().__init__(train_X, train_Y)
self.mean_module = QuadraticMean(input_size=1)
self.covar_module = ScaleKernel(RBFKernel() + PolynomialKernel(power=2))
# 为每个目标建立模型
models = []
for i in range(Y.shape[1]):
model = QuadraticMeanGP(X, Y_std[:, i])
mll = ExactMarginalLogLikelihood(model.likelihood, model)
fit_gpytorch_model(mll)
models.append(model)
from botorch.models import ModelListGP
model_list = ModelListGP(*models)
然后,定义采集函数,例如ParEGO:
from botorch.acquisition import PosteriorSampling
from botorch.acquisition.analytic import ExpectedImprovement
from botorch.utils.transforms import unnormalize, normalize
from botorch.acquisition import FixedFeatureAcquisitionFunction
from botorch.optim import optimize_acqf
def get_weights(Y):
# 生成随机权重,用于ParEGO
weights = torch.rand(Y.shape[-1])
weights /= weights.sum()
return weights
# 定义标量化函数
def scalarize(Y, weights):
return (Y * weights).sum(dim=-1)
# 定义ParEGO采集函数
class ParEGO(AnalyticAcquisitionFunction):
def __init__(self, model, weights, best_f):
super().__init__(model)
self.weights = weights
self.best_f = best_f
def forward(self, X):
posterior = self.model.posterior(X)
mean = posterior.mean
scalarized_mean = scalarize(mean, self.weights)
return (scalarized_mean - self.best_f).clamp_min(0)
但可能需要更正确的实现,或者使用botorch内置的ParEGO。不过可能botorch没有直接提供ParEGO的类,需要手动实现。
或者,使用ScalarizedPosteriorTransform:
from botorch.acquisition import ScalarizedPosteriorTransform
from botorch.acquisition.analytic import ExpectedImprovement
weights = torch.tensor([0.25, 0.25, 0.25, 0.25]) # 示例权重
posterior_transform = ScalarizedPosteriorTransform(weights)
EI = ExpectedImprovement(
model=model_list,
best_f=scalarize(Y_std, weights).max(),
posterior_transform=posterior_transform,
)
但这里可能需要正确配置model_list的输出。
这可能存在一些问题,因为model_list的输出是多个模型,需要将它们的后验合并。或者,需要将多目标转化为单目标,通过线性加权。
因此,正确的方式可能需要使用MultiOutputGPyTorchModel,或者确保后验的合并正确。
可能更简单的是将多个目标模型合并成一个模型,输出多变量,然后应用标量化。或者,使用ModelListGP和ScalarizedPosteriorTransform。
但这里可能比较复杂,建议参考botorch的教程。
或者,改用Monte Carlo采集函数,例如qExpectedImprovement,结合后验采样。
不过时间有限,可能需要采用另一种方法。
例如,在每次迭代中,随机生成一个权重向量,标量化目标,然后进行单目标优化。
具体步骤如下:
1. 训练多目标GP模型。
2. 生成随机权重向量,例如使用Dirichlet分布生成归一化的权重。
3. 计算当前观测数据中,该权重下的最好值(标量化的最大值)。
4. 定义采集函数为Expected Improvement,基于标量化后的目标。
5. 优化采集函数,找到下一个候选点。
6. 评估新点的目标值(实际中需要真实函数,这里用户可能需要模拟)。
但用户的数据是已有的,可能需要在贝叶斯优化过程中建议新的x点,然后根据已有数据预测,或者假设存在未知的真实函数。这里可能用户的数据只是训练数据,用于建模,而贝叶斯优化需要在这些模型的基础上寻找最优的x。
综上,可能代码的大致框架如下:
import torch
import gpytorch
from gpytorch.means import Mean
from gpytorch.kernels import ScaleKernel, RBFKernel, PolynomialKernel
from gpytorch.mlls import ExactMarginalLogLikelihood
from botorch.models import SingleTaskGP, ModelListGP
from botorch.fit import fit_gpytorch_model
from botorch.acquisition import ExpectedImprovement
from botorch.optim import optimize_acqf
from botorch.utils import standardize
import numpy as np
# 自定义二次均值函数
class QuadraticMean(Mean):
def __init__(self, input_size, batch_shape=torch.Size()):
super().__init__()
self.register_parameter(
name="coefficients",
parameter=torch.nn.Parameter(torch.randn(batch_shape, input_size, 1))
)
self.register_parameter(
name="quadratic",
parameter=torch.nn.Parameter(torch.randn(batch_shape, input_size, 1))
)
self.register_parameter(
name="intercept",
parameter=torch.nn.Parameter(torch.randn(batch_shape, 1))
)
def forward(self, x):
return (x ** 2) @ self.quadratic + x @ self.coefficients + self.intercept
# 自定义GP模型
class CustomGP(SingleTaskGP):
def __init__(self, train_X, train_Y):
super().__init__(train_X, train_Y)
self.mean_module = QuadraticMean(input_size=train_X.shape[-1])
self.covar_module = ScaleKernel(
RBFKernel() + PolynomialKernel(power=2)
)
# 生成示例数据
X = torch.tensor([[1.0], [2.0], [3.0], [4.0], [5.0], [6.0]]) # 6个参数点,形状(6,1)
# 每个参数点重复4次实验,共24个样本
X = X.repeat_interleave(4, dim=0)
# 生成对应的四个目标的随机数据,形状(24,4)
Y = torch.randn(24, 4)
# 标准化目标
Y_std = standardize(Y)
# 为每个目标建立独立GP模型
models = []
for i in range(Y.shape[1]):
model = CustomGP(X, Y_std[:, i])
mll = ExactMarginalLogLikelihood(model.likelihood, model)
fit_gpytorch_model(mll)
models.append(model)
model_list = ModelListGP(*models)
# 定义ParEGO采集函数的权重
def get_random_weights(n_obj):
weights = torch.distributions.Dirichlet(torch.ones(n_obj)).sample()
return weights
n_obj = Y.shape[1]
weights = get_random_weights(n_obj)
# 标量化当前最佳值
current_Y = Y_std # 标准化后的观测值
scalarized_Y = (current_Y * weights).sum(dim=-1)
best_f = scalarized_Y.max()
# 定义采集函数(使用ModelListGP的标量化后验)
from botorch.acquisition import AnalyticAcquisitionFunction
from botorch.utils.transforms import concatenate_pending_points, t_batch_mode_transform
class ScalarizedEI(AnalyticAcquisitionFunction):
def __init__(self, model, weights, best_f):
super().__init__(model)
self.weights = weights
self.best_f = best_f
@t_batch_mode_transform(expected_q=1)
def forward(self, X):
posterior = self.model.posterior(X)
mean = posterior.mean # 形状为(n, m),m是目标数
scalarized_mean = (mean * self.weights).sum(dim=-1)
improvement = scalarized_mean - self.best_f
return improvement.clamp_min(0)
acquisition_function = ScalarizedEI(model_list, weights, best_f)
# 优化采集函数
bounds = torch.tensor([[1.0], [6.0]]) # 参数范围
candidates, _ = optimize_acqf(
acq_function=acquisition_function,
bounds=bounds,
q=1,
num_restarts=10,
raw_samples=20,
)
print("Next candidate:", candidates)
但可能需要更多的调整,比如处理模型输入输出的维度,以及正确计算标量化后的EI。
此外,可能需要循环进行多次贝叶斯优化迭代,每次生成新的权重,找到新的候选点,并更新模型。但用户可能希望得到完整的代码,所以需要将上述步骤整合成一个循环。
但用户的问题可能希望得到一个完整的贝叶斯优化代码,包括数据预处理、模型定义、优化循环等。因此,可能需要将上述代码整合,并添加必要的注释和结构。
最终,考虑到时间和复杂度,可能提供一个基于ParEGO方法的代码框架,使用随机权重进行标量化,并迭代优化。但需要确保代码的完整性和正确性。</think>以下是一个基于贝叶斯优化的多目标工艺参数优化实现方案,使用二次多项式均值函数和复合核函数:
```python
import torch
import gpytorch
from gpytorch.means import Mean
from gpytorch.kernels import ScaleKernel, RBFKernel, PolynomialKernel
from gpytorch.mlls import ExactMarginalLogLikelihood
from botorch.models import SingleTaskGP, ModelListGP
from botorch.fit import fit_gpytorch_model
from botorch.optim import optimize_acqf
from botorch.utils import standardize
import numpy as np
# 自定义二次多项式均值函数
class QuadraticMean(Mean):
def __init__(self, input_size):
super().__init__()
self.register_parameter(
name="coefficients",
parameter=torch.nn.Parameter(torch.randn(input_size))
)
self.register_parameter(
name="quad_terms",
parameter=torch.nn.Parameter(torch.randn(input_size))
)
self.register_parameter(
name="bias",
parameter=torch.nn.Parameter(torch.randn(1))
)
def forward(self, x):
return (x**2 @ self.quad_terms) + (x @ self.coefficients) + self.bias
# 自定义高斯过程模型
class CustomGP(SingleTaskGP):
def __init__(self, train_X, train_Y):
super().__init__(train_X, train_Y)
self.mean_module = QuadraticMean(input_size=train_X.shape[-1])
self.covar_module = ScaleKernel(
RBFKernel() + PolynomialKernel(power=2)
)
# 生成示例数据(假设x范围为[1,6])
X = torch.tensor([[1.0], [2.0], [3.0], [4.0], [5.0], [6.0]])
X = X.repeat_interleave(4, dim=0) # 每个参数重复4次实验
Y = torch.cat([ # 生成模拟工艺结果(4个目标)
torch.sin(X) + torch.randn_like(X)*0.1,
torch.cos(X) + torch.randn_like(X)*0.1,
X**0.5 + torch.randn_like(X)*0.1,
-X + torch.randn_like(X)*0.1
], dim=1)
# 数据标准化
Y_std = standardize(Y)
# 为每个目标建立独立GP模型
models = []
for i in range(Y.shape[1]):
model = CustomGP(X, Y_std[:, i])
mll = ExactMarginalLogLikelihood(model.likelihood, model)
fit_gpytorch_model(mll)
models.append(model)
model_list = ModelListGP(*models)
# 帕累托前沿优化参数
NUM_ITERATIONS = 20
BATCH_SIZE = 1
PARAM_BOUNDS = torch.tensor([[1.0], [6.0]])
# 贝叶斯优化主循环
for iteration in range(NUM_ITERATIONS):
# 生成随机权重(使用Dirichlet分布)
weights = torch.distributions.Dirichlet(torch.ones(4)).sample()
# 计算当前最优标量值
current_Y = model_list.posterior(X).mean
阅读全文
相关推荐

















