|
| 1 | +# Learning to regularize with a variational autoencoder for hydrologic inverse analysis |
| 2 | + |
| 3 | +## 1.简介 |
| 4 | + |
| 5 | +本项目基于paddle框架复现。 |
| 6 | + |
| 7 | +论文主要点如下: |
| 8 | + |
| 9 | +* 作者提出了一种基于变分自动编码器 (VAE)的正则化方法; |
| 10 | +* 这种方法的优点1: 对来自VAE的潜在变量(此处指encoder的输出)执行正则化,可以简单地对其进行正则化; |
| 11 | +* 这种方法的优点2: VAE减少了优化问题中的变量数量,从而在伴随方法不可用时使基于梯度的优化在计算上更加高效。 |
| 12 | + |
| 13 | +本项目关键技术要点: |
| 14 | + |
| 15 | +* 实现paddle和julia混写代码梯度传递,避免大面积重写julia代码并能够调用优秀的julia代码; |
| 16 | +* 发现paddle minimize_lbfgs存在问题, 待提交issue确认。 |
| 17 | + |
| 18 | +实验结果要点: |
| 19 | + |
| 20 | +* 成功复现论文代码框架及全流程运行测试; |
| 21 | +* 本次复现精度因无法使用相同样本,无法与论文中数据进行一一比较。本项目给出了采用paddle编写的框架结果。 |
| 22 | + |
| 23 | +论文信息: |
| 24 | +O'Malley D, Golden J K, Vesselinov V V. Learning to regularize with a variational autoencoder for hydrologic inverse analysis[J]. arXiv preprint arXiv:1906.02401, 2019. |
| 25 | + |
| 26 | +参考GitHub地址: |
| 27 | +<https://github.com/madsjulia/RegAE.jl> |
| 28 | + |
| 29 | +项目aistudio地址: |
| 30 | +<https://aistudio.baidu.com/aistudio/projectdetail/5541961> |
| 31 | + |
| 32 | +模型结构 |
| 33 | + |
| 34 | + |
| 35 | +## 2.数据集 |
| 36 | + |
| 37 | +本项目数据集通过作者提供的julia代码生成,生成后保存为npz文件,已上传aistudio[数据集](https://aistudio.baidu.com/aistudio/datasetdetail/193595)并关联本项目。 |
| 38 | + |
| 39 | +以下为数据生成和数据保存代码的说明 |
| 40 | + |
| 41 | +(1)作者通过julia中的DPFEHM和GaussianRandomFields进行数据生成,代码可参考本项目/home/aistudio/RegAE.jl/examples/hydrology/ex_gaussian.jl,可根据其中参数进行修改; |
| 42 | + |
| 43 | +(2)数据保存代码。在/home/aistudio/RegAE.jl/examples/hydrology/ex.jl代码中增加以下代码,可将数据通过转换为numpy数据并保存为npz。 |
| 44 | + |
| 45 | +```julia |
| 46 | +using Distributed |
| 47 | +using PyCall # 增加此处引用 |
| 48 | + |
| 49 | +@everywhere variablename = "allloghycos" |
| 50 | +@everywhere datafilename = "$(results_dir)/trainingdata.jld2" |
| 51 | +if !isfile(datafilename) |
| 52 | + if nworkers() == 1 |
| 53 | + error("Please run in parallel: julia -p 32") |
| 54 | + end |
| 55 | + numsamples = 10^5 |
| 56 | + @time allloghycos = SharedArrays.SharedArray{Float32}(numsamples, ns[2], ns[1]; init=A -> samplehyco!(A; setseed=true)) |
| 57 | + # @time JLD2.@save datafilename allloghycos |
| 58 | + |
| 59 | + ########### 此处为增加部分 ########### |
| 60 | + p_trues = SharedArrays.SharedArray{Float32}(3, ns[2], ns[1]; init=samplehyco!) # 计算p_true |
| 61 | + |
| 62 | + np = pyimport("numpy") |
| 63 | + training_data = np.asarray(allloghycos) |
| 64 | + test_data = np.asarray(p_trues) |
| 65 | + |
| 66 | + np_coords = np.asarray(coords) |
| 67 | + np_neighbors = np.asarray(neighbors) |
| 68 | + np_areasoverlengths = np.asarray(areasoverlengths) |
| 69 | + np_dirichletnodes = np.asarray(dirichletnodes) |
| 70 | + np_dirichletheads = np.asarray(dirichletheads) |
| 71 | + |
| 72 | + np.savez("$(results_dir)/gaussian_train.npz", |
| 73 | + data=training_data, |
| 74 | + test_data=test_data, |
| 75 | + coords=np_coords, |
| 76 | + neighbors=np_neighbors, |
| 77 | + areasoverlengths=np_areasoverlengths, |
| 78 | + dirichletnodes=np_dirichletnodes, |
| 79 | + dirichletheads=np_dirichletheads) |
| 80 | +end |
| 81 | +``` |
| 82 | + |
| 83 | +## 数据标准化 |
| 84 | + |
| 85 | +* 数据标准化方式: $z = (x - \mu)/ \sigma$ |
| 86 | + |
| 87 | +```python |
| 88 | + class ScalerStd(object): |
| 89 | + """ |
| 90 | + Desc: Normalization utilities with std mean |
| 91 | + """ |
| 92 | + |
| 93 | + def __init__(self): |
| 94 | + self.mean = 0. |
| 95 | + self.std = 1. |
| 96 | + |
| 97 | + def fit(self, data): |
| 98 | + self.mean = np.mean(data) |
| 99 | + self.std = np.std(data) |
| 100 | + |
| 101 | + def transform(self, data): |
| 102 | + mean = paddle.to_tensor(self.mean).type_as(data).to( |
| 103 | + data.device) if paddle.is_tensor(data) else self.mean |
| 104 | + std = paddle.to_tensor(self.std).type_as(data).to( |
| 105 | + data.device) if paddle.is_tensor(data) else self.std |
| 106 | + return (data - mean) / std |
| 107 | + |
| 108 | + def inverse_transform(self, data): |
| 109 | + mean = paddle.to_tensor(self.mean) if paddle.is_tensor(data) else self.mean |
| 110 | + std = paddle.to_tensor(self.std) if paddle.is_tensor(data) else self.std |
| 111 | + return (data * std) + mean |
| 112 | +``` |
| 113 | + |
| 114 | +## 定义Dataset |
| 115 | + |
| 116 | +1. 通过读取预保存npz加载数据,当前数据类型为 [data_nums, 100, 100], 此处100为数据生成过程中指定 |
| 117 | +2. 数据reshape为 [data_nums, 10000] |
| 118 | +3. 数据划分为8:2用与train和test |
| 119 | +4. 通过对train数据得到标准化参数mean和std,并用此参数标准化train和test数据集 |
| 120 | +5. 通过dataloader得到的数据形式为 [batch_size, 10000] |
| 121 | + |
| 122 | +```python |
| 123 | + class CustomDataset(Dataset): |
| 124 | + def __init__(self, file_path, data_type="train"): |
| 125 | + """ |
| 126 | +
|
| 127 | + :param file_path: |
| 128 | + :param data_type: train or test |
| 129 | + """ |
| 130 | + super().__init__() |
| 131 | + all_data = np.load(file_path) |
| 132 | + data = all_data["data"] |
| 133 | + num, _, _ = data.shape |
| 134 | + data = data.reshape(num, -1) |
| 135 | + |
| 136 | + self.neighbors = all_data['neighbors'] |
| 137 | + self.areasoverlengths = all_data['areasoverlengths'] |
| 138 | + self.dirichletnodes = all_data['dirichletnodes'] |
| 139 | + self.dirichleths = all_data['dirichletheads'] |
| 140 | + self.Qs = np.zeros([all_data['coords'].shape[-1]]) |
| 141 | + self.val_data = all_data["test_data"] |
| 142 | + |
| 143 | + self.data_type = data_type |
| 144 | + |
| 145 | + self.train_len = int(num * 0.8) |
| 146 | + self.test_len = num - self.train_len |
| 147 | + |
| 148 | + self.train_data = data[:self.train_len] |
| 149 | + self.test_data = data[self.train_len:] |
| 150 | + |
| 151 | + self.scaler = ScalerStd() |
| 152 | + self.scaler.fit(self.train_data) |
| 153 | + |
| 154 | + self.train_data = self.scaler.transform(self.train_data) |
| 155 | + self.test_data = self.scaler.transform(self.test_data) |
| 156 | + |
| 157 | + def __getitem__(self, idx): |
| 158 | + if self.data_type == "train": |
| 159 | + return self.train_data[idx] |
| 160 | + else: |
| 161 | + return self.test_data[idx] |
| 162 | + |
| 163 | + def __len__(self): |
| 164 | + if self.data_type == "train": |
| 165 | + return self.train_len |
| 166 | + else: |
| 167 | + return self.test_len |
| 168 | +``` |
| 169 | + |
| 170 | +## 将数据转换为IterableNPZDataset的形式 |
| 171 | + |
| 172 | +```python |
| 173 | +np.savez("data.npz", p_train=train_data.train_data, p_test=train_data.test_data) |
| 174 | +``` |
| 175 | + |
| 176 | +## 3.环境依赖 |
| 177 | + |
| 178 | +本项目为julia和python混合项目。 |
| 179 | + |
| 180 | +### julia依赖 |
| 181 | + |
| 182 | +* DPFEHM |
| 183 | +* Zygote |
| 184 | + |
| 185 | +### python依赖 |
| 186 | + |
| 187 | +* paddle |
| 188 | +* julia (pip安装) |
| 189 | +* matplotlib |
| 190 | + |
| 191 | +本项目已经提供安装后压缩文档,可fork本项目后执行以下代码进行解压安装。 |
| 192 | + |
| 193 | +```python |
| 194 | +# 解压预下载文件和预编译文件 |
| 195 | +!tar zxf /home/aistudio/opt/curl-7.88.1.tar.gz -C /home/aistudio/opt # curl 预下载文件 |
| 196 | +!tar zxf /home/aistudio/opt/curl-7.88.1-build.tgz -C /home/aistudio/opt # curl 预编译文件 |
| 197 | +!tar zxf /home/aistudio/opt/julia-1.8.5-linux-x86_64.tar.gz -C /home/aistudio/opt # julia 预下载文件 |
| 198 | +!tar zxf /home/aistudio/opt/julia_package.tgz -C /home/aistudio/opt # julia依赖库文件 |
| 199 | +!tar zxf /home/aistudio/opt/external-libraries.tgz -C /home/aistudio/opt # pip依赖库文件 |
| 200 | +``` |
| 201 | + |
| 202 | +```python |
| 203 | +####### 以下指令需要时可参考执行,上述压缩包已经完成以下内容 ####### |
| 204 | + |
| 205 | +# curl 编译指令,当解压后无效使用 |
| 206 | +!mkdir -p /home/aistudio/opt/curl-7.88.1-build |
| 207 | +!/home/aistudio/opt/curl-7.88.1/configure --prefix=/home/aistudio/opt/curl-7.88.1-build --with-ssl --enable-tls-srp |
| 208 | +!make install -j4 |
| 209 | + |
| 210 | +# 指定curl预编译文件 |
| 211 | +!export LD_LIBRARY_PATH=/home/aistudio/opt/curl-7.88.1-build/lib:$LD_LIBRARY_PATH |
| 212 | +!export PATH=/home/aistudio/opt/curl-7.88.1-build/bin:$PATH |
| 213 | +!export CPATH=/home/aistudio/opt/curl-7.88.1-build/include:$CPATH |
| 214 | +!export LIBRARY_PATH=/home/aistudio/opt/curl-7.88.1-build/lib:$LIBRARY_PATH |
| 215 | + |
| 216 | +# 指定已经安装的julia包 |
| 217 | +!export JULIA_DEPOT_PATH=/home/aistudio/opt/julia_package |
| 218 | +# 指定julia使用清华源 |
| 219 | +!export JULIA_PKG_SERVER=https://mirrors.tuna.tsinghua.edu.cn/julia |
| 220 | +# julia 安装依赖库 |
| 221 | +# 需要先export JULIA_DEPOT_PATH 环境变量,否则安装位置为~/.julia, aistudio无法保存 |
| 222 | +!/home/aistudio/opt/julia-1.8.5/bin/julia -e "using Pkg; Pkg.add(\"DPFEHM\")" |
| 223 | +!/home/aistudio/opt/julia-1.8.5/bin/julia -e "using Pkg; Pkg.add(\"Zygote\")" |
| 224 | +!/home/aistudio/opt/julia-1.8.5/bin/julia -e "using Pkg; Pkg.add(\"PyCall\")" |
| 225 | +``` |
| 226 | + |
| 227 | +使用方法可以参考以下代码和julia导数传递测试.ipynb文件。 |
| 228 | + |
| 229 | +```python |
| 230 | +import paddle |
| 231 | +import os |
| 232 | +import sys |
| 233 | + |
| 234 | +# julia 依赖 |
| 235 | +os.environ['JULIA_DEPOT_PATH'] = '/home/aistudio/opt/julia_package' |
| 236 | +# pip 依赖 |
| 237 | +sys.path.append('/home/aistudio/opt/external-libraries') |
| 238 | + |
| 239 | +# julieries |
| 240 | +from julia.api import Julia |
| 241 | + |
| 242 | +jl = Julia(compiled_modules=False,runtime="/home/aistudio/opt/julia-1.8.5/bin/julia") |
| 243 | +# import julia |
| 244 | +from julia import Main |
| 245 | +``` |
| 246 | + |
| 247 | +## 4.快速开始 |
| 248 | + |
| 249 | +本项目运行分为两个步骤: |
| 250 | + |
| 251 | +* (1)训练步骤。通过运行train.ipynb文件,可以得到训练后的模型参数,具体代码请参考train.ipynb文件及其中注释说明; |
| 252 | +* (2)测试步骤。通过运行test.ipynb文件,应用训练后的模型参数,对latent domain进行优化。 |
| 253 | + |
| 254 | +## 5.代码结构与详细说明 |
| 255 | + |
| 256 | +```text |
| 257 | +├── data #预生成数据文件 |
| 258 | +│ └── data193595 |
| 259 | +├── main.ipynb #本说明文件 |
| 260 | +├── opt #环境配置文件,已压缩,解压即可使用 |
| 261 | +│ ├── curl-7.88.1 |
| 262 | +│ ├── curl-7.88.1-build |
| 263 | +│ ├── curl-7.88.1-build.tgz |
| 264 | +│ ├── curl-7.88.1.tar.gz |
| 265 | +│ ├── external-libraries |
| 266 | +│ ├── external-libraries.tgz |
| 267 | +│ ├── julia-1.8.5 |
| 268 | +│ ├── julia-1.8.5-linux-x86_64.tar.gz |
| 269 | +│ ├── julia_package |
| 270 | +│ └── julia_package.tgz |
| 271 | +├── params_vae_nz100 #模型参数文件 |
| 272 | +│ └── model.pdparams |
| 273 | +├── params_vae_nz200 |
| 274 | +│ └── model.pdparams |
| 275 | +├── params_vae_nz400 |
| 276 | +│ └── model.pdparams |
| 277 | +├── test.ipynb #测试文件 |
| 278 | +├── train.ipynb #训练文件 |
| 279 | +├── julia导数传递测试.ipynb #julia和python混合测试文件 |
| 280 | +``` |
| 281 | + |
| 282 | +### train文件和test文件关联性说明 |
| 283 | + |
| 284 | +我们依照论文作者的符号进行说明,$p$为数据输入,$\hat{p}$为数据输出,$loss=mse(p,\hat{p}) + loss_{kl}(\hat{p},N(0,1))$。 |
| 285 | + |
| 286 | +* (1)通过train能够得到训练后的Autoencoder(包含encoder和decoder); |
| 287 | +* (2)通过test调用训练后的encoder针对testdata生成latent_test,并得到latent_mean; |
| 288 | +* (3)针对新生成的样本$p_{new}$,通过LBFGS方法不断优化latent_mean,直到obj_fun最小,其中obj_fun = mse($p_{new}$,$\hat{p}_{new}$)+mse(sci_fun($p_{new}$),sci_fun($\hat{p}_{new}$)),sci_fun为任何其他科学计算模拟方法。 |
| 289 | + |
| 290 | +### paddle.incubate.optimizer.functional.minimize_lbfgs 问题 |
| 291 | + |
| 292 | +以下为paddle官方minimize_lbfgs API: |
| 293 | + |
| 294 | +```python |
| 295 | +paddle.incubate.optimizer.functional.minimize_lbfgs(objective_func, initial_position, history_size=100, max_iters=50, tolerance_grad=1e-08, tolerance_change=1e-08, initial_inverse_hessian_estimate=None, line_search_fn='strong_wolfe', max_line_search_iters=50, initial_step_length=1.0, dtype='float32', name=None) |
| 296 | +``` |
| 297 | + |
| 298 | +* (1)参数max_line_search_iters无效。虽然设置了此参数,但是内部没有传递对应参数; |
| 299 | +* (2)中wolfe条件1错误。line256处应为`phi_2 >= phi_1`,以下为paddle部分源码。 |
| 300 | + |
| 301 | +```python |
| 302 | + # 1. If phi(a2) > phi(0) + c_1 * a2 * phi'(0) or [phi(a2) >= phi(a1) and i > 1], |
| 303 | + # a_star= zoom(a1, a2) and stop; |
| 304 | + pred1 = ~done & ((phi_2 > phi_0 + c1 * a2 * derphi_0) | |
| 305 | + ((phi_2 >= phi_0) & (i > 1))) |
| 306 | +``` |
| 307 | + |
| 308 | +## 6.复现结果 |
| 309 | + |
| 310 | +### 不同latent维度对比 |
| 311 | + |
| 312 | + |
| 313 | + |
| 314 | +通过实验结果可以发现: |
| 315 | + |
| 316 | +* (1)不同样本之间存在差距,并不是所有样本都能优化得到良好的latent变量; |
| 317 | +* (2)随着模型latent维度的上升,模型效果逐渐提升。 |
| 318 | + |
| 319 | +### latent_random和latent_mean对比 |
| 320 | + |
| 321 | +本项目还增加了latent_random和latent_mean对生成结果的对比。此处对latent_random和latent_mean再次说明: |
| 322 | + |
| 323 | +* latent_random:通过paddle.randn生成的高斯噪声得到; |
| 324 | +* latent_mean:通过对所有testdata进行encoder结果平均得到。 |
| 325 | + |
| 326 | +以下为通过latent_random得到的实验结果 |
| 327 | + |
| 328 | + |
| 329 | +通过对比,可以发现latent_mean对优化结果重要影响。近似正确的latent变量能够得到更优的生成结果。 |
| 330 | + |
| 331 | +### LBFGS优化收敛情况 |
| 332 | + |
| 333 | +可以从如下图中看出,使用paddle minimize_lbfgs能够有效优化收敛。 |
| 334 | + |
| 335 | + |
| 336 | +## 7.延伸思考 |
| 337 | + |
| 338 | +如果深入思考本项目,会发现模型在test过程中是使用真实数据作为目标进行lbfgs优化,这种计算方式还有意义吗? |
| 339 | + |
| 340 | +回答是肯定的!有意义! |
| 341 | + |
| 342 | +以下为本人个人观点: |
| 343 | + |
| 344 | +* (1)通过实验对比latent_random和latent_mean的最终生成结果差距,可以发现一个良好的初值对模型的影响是巨大的。当前diffusion模型在sample生成过程对生成的高斯噪声不断做denoise操作,这其中生成的噪声数据如果经过预先优化,不仅能够加速diffusion的生成速度,而且能够提升数据的生成质量。 |
| 345 | +* (2)在域迁移等研究领域,可以使用这种latent逐渐生成中间过渡变量,达到不同域数据的迁移生成。 |
| 346 | + |
| 347 | +## 7.模型信息 |
| 348 | + |
| 349 | +| 信息 | 说明| |
| 350 | +| -------- | -------- | |
| 351 | +| 发布者 | 朱卫国 (DrownFish19) | |
| 352 | +| 发布时间 | 2023.03 | |
| 353 | +| 框架版本 | paddle 2.4.1 | |
| 354 | +| 支持硬件 | GPU、CPU | |
| 355 | +| aistudio | [notebook](https://aistudio.baidu.com/aistudio/projectdetail/5541961) | |
| 356 | + |
| 357 | +请点击[此处](https://ai.baidu.com/docs#/AIStudio_Project_Notebook/a38e5576)查看本环境基本用法. |
| 358 | +Please click [here](https://ai.baidu.com/docs#/AIStudio_Project_Notebook/a38e5576) for more detailed instructions. |
0 commit comments