Skip to content

Commit a778938

Browse files
authored
Weak lensing: decouple from master, refactor descwl/anacal workflow (#1104)
* Refactor case_studies/weak_lensing * Fix refactor bug and move wl/cached_dataset * Update dc2 config * Allow count splits instead of percent splits * Remove deprecated files for generating encoder input * Script for running AnaCal * Update anacal script * Ignore anacal psf/mask/etc during NPE training * Update config * Move some descwl notebooks to deprecated folder * Deprecated scripts * Rename and reorganize descwl * Rename and reorganize DC2 * Rename descwl images notebook * Rename config * Update config paths * Move descwl ckpt * Update descwl NPE credible intervals script * Update descwl NPE credible intervals script * Speed up AnaCal and run on same test set as NPE * Get rid of as_completed in AnaCal multiprocessing * Update descwl credible intervals notebook * Fix and run all DC2 notebooks * Update descwl image notebook * Update descwl scatterplots * Update README * ruff + update gitignore
1 parent 39faef5 commit a778938

48 files changed

Lines changed: 5719 additions & 5875 deletions

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

.gitignore

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -32,10 +32,8 @@ data/dc2/merged_catalog
3232
case_studies/sdss_galaxies/models/simulated_blended_galaxies.pt
3333
case_studies/*/data/
3434
venv
35-
nohup.out
35+
*.out
3636
multirun
37-
DC2_*.out
3837
input_images/
39-
WeakLensingResults*
40-
ResultsDC2*
38+
Results*
4139
case_studies/weak_lensing/**/*.pt
Lines changed: 34 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,19 +1,44 @@
1-
### Neural posterior estimation for tomographic field-level weak lensing inference
2-
#### Tim White, Shreyas Chandrashekaran, Dingrui Tao, Camille Avestruz, and Jeffrey Regier
3-
#### with assistance from Steve Fan and Tahseen Younus
1+
## Neural posterior estimation for field-level weak lensing inference
2+
#### Tim White, Shreyas Chandrashekaran, Dingrui Tao, Camille Avestruz, and Jeffrey Regier, with assistance from Steve Fan and Tahseen Younus
43

5-
In this case study, we use neural posterior estimation to infer tomographic shear and convergence maps from LSST-like images.
4+
In this case study, we use neural posterior estimation (NPE) to infer weak lensing shear and convergence from LSST-like images. We use NPE to (1) infer tomographic mass maps for the [DC2 Simulated Sky Survey](https://data.lsstdesc.org/doc/dc2_sim_sky_survey) and (2) infer constant shear from images generated with the [`descwl-shear-sims` package](https://github.com/timwhite0/descwl-shear-sims).
65

7-
To train the encoder on DC2 images, run
6+
## DC2
7+
8+
### Generate catalog
9+
10+
```
11+
nohup python -u case_studies/weak_lensing/dc2/generate_catalog.py &> generate_catalog.out &
12+
```
13+
14+
### Generate mass maps and train MassMapEncoder
15+
16+
```
17+
nohup bliss -cp <path>/bliss/case_studies/weak_lensing/dc2 -cn config_train_npe.yaml mode=train &> train_dc2.out &
18+
```
19+
20+
### Notebooks
21+
22+
- **In `dc2/results`**: `credibleintervals.ipynb`, `posteriormeanmaps.ipynb`
23+
- **In `dc2/exploratory`**: `dc2imageandmaps.ipynb`, `ellipticity.ipynb`, `galaxyproperties.ipynb`, `twopoint.ipynb`
24+
25+
## descwl-shear-sims
26+
27+
### Train ScalarShearEncoder
828

929
```
10-
nohup bliss -cp <path>/bliss/case_studies/weak_lensing/dc2 -cn config_dc2.yaml mode=train &> train_on_dc2.out &
30+
nohup bliss -cp <path>/bliss/case_studies/weak_lensing/descwl -cn config_train_npe.yaml mode=train &> train_descwl.out &
1131
```
1232

13-
To train the encoder on descwl-shear-sims images, run
33+
### Run [AnaCal](https://github.com/mr-superonion/AnaCal)
34+
35+
Configure settings in `descwl/config_run_anacal.yaml`.
1436

1537
```
16-
nohup bliss -cp <path>/bliss/case_studies/weak_lensing/descwl -cn config_descwl.yaml mode=train &> train_on_descwl.out &
38+
nohup python -u case_studies/weak_lensing/descwl/run_anacal.py &> run_anacal.out &
1739
```
1840

19-
See `dc2/notebooks` and `descwl/notebooks` for some exploratory plots and our most recent results.
41+
### Notebooks
42+
43+
- **In `descwl/results`**: `compute_npe_credibleintervals.py`, `credibleintervals.ipynb`, `scatterplots.ipynb`
44+
- **In `descwl/exploratory`**: `images.ipynb`

case_studies/weak_lensing/bliss_simulator/angular_cl.npy renamed to case_studies/weak_lensing/bliss_simulator/deprecated_scripts/angular_cl.npy

File renamed without changes.

case_studies/weak_lensing/bliss_simulator/decoder.py renamed to case_studies/weak_lensing/bliss_simulator/deprecated_scripts/decoder.py

File renamed without changes.

case_studies/weak_lensing/bliss_simulator/generate_angular_cl.py renamed to case_studies/weak_lensing/bliss_simulator/deprecated_scripts/generate_angular_cl.py

File renamed without changes.

case_studies/weak_lensing/bliss_simulator/prior.py renamed to case_studies/weak_lensing/bliss_simulator/deprecated_scripts/prior.py

File renamed without changes.
Lines changed: 267 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,267 @@
1+
"""Cached dataset infrastructure for weak lensing case study.
2+
3+
This is a standalone copy of the data loading infrastructure from bliss.cached_dataset,
4+
modified to remove GlobalEnv dependency.
5+
"""
6+
7+
import functools
8+
import logging
9+
import math
10+
import operator
11+
import os
12+
import pathlib
13+
import random
14+
import re
15+
import warnings
16+
from typing import List
17+
18+
import pytorch_lightning as pl
19+
import torch
20+
from torch import distributed as dist
21+
from torch.utils.data import DataLoader, Dataset, DistributedSampler, Sampler
22+
from torchvision import transforms
23+
24+
# Suppress pytorch_lightning warnings
25+
warnings.filterwarnings(
26+
"ignore", ".*does not have many workers which may be a bottleneck.*", UserWarning
27+
)
28+
warnings.filterwarnings("ignore", ".*Total length of .* across ranks is zero.*", UserWarning)
29+
30+
31+
class ChunkingSampler(Sampler):
32+
"""Sampler that respects chunked data ordering."""
33+
34+
def __init__(self, dataset: Dataset) -> None:
35+
super().__init__()
36+
assert isinstance(dataset, ChunkingDataset), "dataset should be ChunkingDataset"
37+
self.dataset = dataset
38+
39+
def __len__(self):
40+
return len(self.dataset)
41+
42+
def __iter__(self):
43+
return iter(self.dataset.get_chunked_indices())
44+
45+
46+
class DistributedChunkingSampler(DistributedSampler):
47+
"""Distributed sampler that respects chunked data ordering."""
48+
49+
def __init__(
50+
self,
51+
dataset: Dataset,
52+
num_replicas: int | None = None,
53+
rank: int | None = None,
54+
shuffle: bool = False,
55+
seed: int = 0,
56+
drop_last: bool = False,
57+
) -> None:
58+
assert isinstance(dataset, ChunkingDataset), "dataset should be ChunkingDataset"
59+
assert not shuffle, "you should not use shuffle"
60+
super().__init__(dataset, num_replicas, rank, shuffle, seed, drop_last)
61+
62+
def __iter__(self):
63+
pre_indices = super().__iter__()
64+
chunked_indices = self.dataset.get_chunked_indices()
65+
return iter([chunked_indices[i] for i in pre_indices])
66+
67+
68+
class ChunkingDataset(Dataset):
69+
"""Dataset that loads data from chunked .pt files."""
70+
71+
def __init__(self, file_paths, shuffle=False, transform=None, seed=0) -> None:
72+
super().__init__()
73+
self.file_paths = file_paths
74+
self.shuffle = shuffle
75+
self.transform = transform
76+
self.seed = seed
77+
self.current_epoch = 0 # Updated by datamodule
78+
79+
self.accumulated_file_sizes = torch.zeros(len(self.file_paths), dtype=torch.int64)
80+
for i, file_path in enumerate(self.file_paths):
81+
file_size_match = re.search(r"size_(\d+)", file_path)
82+
if file_size_match:
83+
cached_data_len = int(file_size_match.group(1))
84+
else:
85+
if i == 0:
86+
logger = logging.getLogger("ChunkingDataset")
87+
warning_msg = (
88+
"WARNING: add postfix '_size_<chunk size>' to file name; "
89+
"otherwise it'll be very slow \n"
90+
)
91+
logger.warning(warning_msg)
92+
with open(file_path, "rb") as f:
93+
cached_data_len = len(torch.load(f, weights_only=False))
94+
95+
if i == 0:
96+
self.accumulated_file_sizes[i] = cached_data_len
97+
else:
98+
self.accumulated_file_sizes[i] = (
99+
self.accumulated_file_sizes[i - 1] + cached_data_len
100+
)
101+
102+
self.buffered_file_index = None
103+
self.buffered_data = None
104+
105+
def __len__(self):
106+
return self.accumulated_file_sizes[-1].item()
107+
108+
def __getitem__(self, index):
109+
converted_index = (self.accumulated_file_sizes <= index).sum().item()
110+
converted_sub_index = (index - self.accumulated_file_sizes[converted_index]).item()
111+
if self.buffered_file_index != converted_index:
112+
self.buffered_file_index = converted_index
113+
with open(self.file_paths[converted_index], "rb") as f:
114+
self.buffered_data = torch.load(f, weights_only=False)
115+
output_data = self.buffered_data[converted_sub_index]
116+
return self.transform(output_data)
117+
118+
def get_chunked_indices(self):
119+
"""Get indices respecting chunk boundaries, with optional shuffling."""
120+
accumulated_file_sizes_list = self.accumulated_file_sizes.tolist()
121+
122+
output_list = []
123+
if self.shuffle:
124+
# Use seed + epoch for reproducible shuffling
125+
epoch_seed = self.seed + self.current_epoch
126+
logger = logging.getLogger("ChunkingDataset")
127+
logger.info(
128+
"INFO: seed is %d; current epoch is %d; epoch_seed is set to %d",
129+
self.seed,
130+
self.current_epoch,
131+
epoch_seed,
132+
)
133+
right_shift_list = [0, *accumulated_file_sizes_list[:-1]]
134+
for start, end in zip(right_shift_list, accumulated_file_sizes_list, strict=True):
135+
rng = random.Random(epoch_seed)
136+
output_list.append(rng.sample(range(start, end), end - start))
137+
random.Random(epoch_seed).shuffle(output_list)
138+
return functools.reduce(operator.iadd, output_list, [])
139+
140+
return list(range(0, len(self)))
141+
142+
143+
class CachedSimulatedDataModule(pl.LightningDataModule):
144+
"""DataModule for loading cached simulation data from .pt files."""
145+
146+
def __init__(
147+
self,
148+
splits: str,
149+
batch_size: int,
150+
num_workers: int,
151+
cached_data_path: str,
152+
train_transforms: List,
153+
nontrain_transforms: List,
154+
subset_fraction: float = None,
155+
shuffle_file_order: bool = True,
156+
seed: int = 0,
157+
splits_type: str = "percent",
158+
):
159+
super().__init__()
160+
161+
self.splits = splits
162+
self.batch_size = batch_size
163+
self.num_workers = num_workers
164+
self.cached_data_path = pathlib.Path(cached_data_path)
165+
self.train_transforms = train_transforms
166+
self.nontrain_transforms = nontrain_transforms
167+
self.subset_fraction = subset_fraction
168+
self.shuffle_file_order = shuffle_file_order
169+
self.seed = seed
170+
self.splits_type = splits_type
171+
172+
self.file_paths = None
173+
self.slices = None
174+
self.train_dataset = None
175+
self.val_dataset = None
176+
self.test_dataset = None
177+
self.predict_dataset = None
178+
179+
def setup(self, stage: str) -> None:
180+
if self.file_paths is None or self.slices is None:
181+
self._load_file_paths_and_slices()
182+
183+
if stage == "fit":
184+
self.train_dataset = self._get_dataset(
185+
self.file_paths[self.slices[0]], self.train_transforms, shuffle=True
186+
)
187+
self.val_dataset = self._get_dataset(
188+
self.file_paths[self.slices[1]], self.nontrain_transforms
189+
)
190+
return None
191+
192+
if stage == "validate":
193+
if self.val_dataset is None:
194+
self.val_dataset = self._get_dataset(
195+
self.file_paths[self.slices[1]], self.nontrain_transforms
196+
)
197+
return None
198+
199+
if stage == "test":
200+
self.test_dataset = self._get_dataset(
201+
self.file_paths[self.slices[2]], self.nontrain_transforms
202+
)
203+
return None
204+
205+
if stage == "predict":
206+
self.predict_dataset = self._get_dataset(self.file_paths, self.nontrain_transforms)
207+
return None
208+
209+
raise RuntimeError(f"setup skips stage {stage}")
210+
211+
def _load_file_paths_and_slices(self):
212+
file_names = [
213+
f for f in sorted(os.listdir(str(self.cached_data_path))) if f.endswith(".pt")
214+
]
215+
if self.shuffle_file_order:
216+
random.shuffle(file_names)
217+
if self.subset_fraction:
218+
file_names = file_names[: math.ceil(len(file_names) * self.subset_fraction)]
219+
self.file_paths = [os.path.join(str(self.cached_data_path), f) for f in file_names]
220+
221+
self.slices = self.parse_slices(self.splits, len(self.file_paths), self.splits_type)
222+
223+
def _percent_to_idx(self, x, length):
224+
"""Converts string in percent to an integer index."""
225+
return int(float(x.strip()) / 100 * length) if x.strip() else None
226+
227+
def _count_to_idx(self, x):
228+
"""Converts string count to an integer index."""
229+
return int(x.strip()) if x.strip() else None
230+
231+
def parse_slices(self, splits: str, length: int, splits_type: str = "percent"):
232+
slices = [slice(0, 0) for _ in range(3)]
233+
for i, data_split in enumerate(splits.split("/")):
234+
if splits_type == "percent":
235+
slices[i] = slice(
236+
*(self._percent_to_idx(val, length) for val in data_split.split(":"))
237+
)
238+
else: # count
239+
slices[i] = slice(*(self._count_to_idx(val) for val in data_split.split(":")))
240+
return slices
241+
242+
def _get_dataset(self, sub_file_paths, defined_transforms, shuffle: bool = False):
243+
assert sub_file_paths, "No cached data found"
244+
transform = transforms.Compose(defined_transforms)
245+
return ChunkingDataset(sub_file_paths, shuffle=shuffle, transform=transform, seed=self.seed)
246+
247+
def _get_dataloader(self, my_dataset):
248+
distributed_is_used = dist.is_available() and dist.is_initialized()
249+
sampler_type = DistributedChunkingSampler if distributed_is_used else ChunkingSampler
250+
return DataLoader(
251+
my_dataset,
252+
batch_size=self.batch_size,
253+
num_workers=self.num_workers,
254+
sampler=sampler_type(my_dataset),
255+
)
256+
257+
def train_dataloader(self):
258+
return self._get_dataloader(self.train_dataset)
259+
260+
def val_dataloader(self):
261+
return self._get_dataloader(self.val_dataset)
262+
263+
def test_dataloader(self):
264+
return self._get_dataloader(self.test_dataset)
265+
266+
def predict_dataloader(self):
267+
return self._get_dataloader(self.predict_dataset)

0 commit comments

Comments
 (0)