From 4ffcfad909820da892fb04668ece2e9df1661982 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Wed, 1 Dec 2021 14:37:12 +0100 Subject: [PATCH 01/63] black changes, vlba as default --- vipy/fits/writer.py | 13 +++--- vipy/simulation/scan.py | 88 ++++++++++++++++++++++------------------ vipy/simulation/utils.py | 2 +- 3 files changed, 57 insertions(+), 46 deletions(-) diff --git a/vipy/fits/writer.py b/vipy/fits/writer.py index 859e7c2..4461edc 100644 --- a/vipy/fits/writer.py +++ b/vipy/fits/writer.py @@ -7,7 +7,7 @@ import astropy.constants as const -def create_vis_hdu(data, conf, layout="EHT", source_name="sim-source-0"): +def create_vis_hdu(data, conf, layout="vlba", source_name="sim-source-0"): u = data.u v = data.v @@ -32,8 +32,11 @@ def create_vis_hdu(data, conf, layout="EHT", source_name="sim-source-0"): vis = np.swapaxes( np.swapaxes( np.stack([values.real, values.imag, np.ones(values.shape)], axis=1), 1, 2 - ),0,1,).reshape(-1, 1, 1, 1, 1, 4, 3) - DATA = vis + ), + 0, + 1, + ).reshape(-1, 1, 1, 1, 1, 4, 3) + DATA = vis # in dim 4 = IFs , dim = 1, dim 4 = number of jones, 3 = real, imag, weight # wcs @@ -232,7 +235,7 @@ def create_frequency_hdu(conf): return hdu_freq -def create_antenna_hdu(layout_txt, conf, layout="EHT"): +def create_antenna_hdu(layout_txt, conf, layout="vlba"): array = pd.read_csv(layout_txt, sep=" ") ANNAME = np.chararray(len(array), itemsize=8, unicode=True) @@ -370,7 +373,7 @@ def create_antenna_hdu(layout_txt, conf, layout="EHT"): return hdu_ant -def create_hdu_list(data, conf, path="../vipy/layouts/eht.txt"): +def create_hdu_list(data, conf, path="../vipy/vipy/layouts/vlba.txt"): vis_hdu = create_vis_hdu(data, conf) time_hdu = create_time_hdu(data) freq_hdu = create_frequency_hdu(conf) diff --git a/vipy/simulation/scan.py b/vipy/simulation/scan.py index 7bae824..2e79749 100644 --- a/vipy/simulation/scan.py +++ b/vipy/simulation/scan.py @@ -202,12 +202,16 @@ def rd_grid(fov, samples, src_crd): 3d array Returns a 3d array with every pixel containing a RA and Dec value """ - res = fov/samples + res = fov / samples - rd_grid = np.zeros((samples,samples,2)) + rd_grid = np.zeros((samples, samples, 2)) for i in range(samples): - rd_grid[i,:,0] = np.array([(i-samples/2)*res + src_crd.ra.rad for i in range(samples)]) - rd_grid[:,i,1] = np.array([-(i-samples/2)*res + src_crd.dec.rad for i in range(samples)]) + rd_grid[i, :, 0] = np.array( + [(i - samples / 2) * res + src_crd.ra.rad for i in range(samples)] + ) + rd_grid[:, i, 1] = np.array( + [-(i - samples / 2) * res + src_crd.dec.rad for i in range(samples)] + ) return rd_grid @@ -228,9 +232,13 @@ def lm_grid(rd_grid, src_crd): Returns a 3d array with every pixel containing a l and m value """ lm_grid = np.zeros(rd_grid.shape) - lm_grid[:,:,0] = np.cos(rd_grid[:,:,1]) * np.sin(rd_grid[:,:,0] - src_crd.ra.rad) - lm_grid[:,:,1] = np.sin(rd_grid[:,:,1]) * np.cos(src_crd.dec.rad) - np.cos(src_crd.dec.rad) * np.sin(src_crd.dec.rad) * np.cos(rd_grid[:,:,0] - src_crd.ra.rad) - + lm_grid[:, :, 0] = np.cos(rd_grid[:, :, 1]) * np.sin( + rd_grid[:, :, 0] - src_crd.ra.rad + ) + lm_grid[:, :, 1] = np.sin(rd_grid[:, :, 1]) * np.cos(src_crd.dec.rad) - np.cos( + src_crd.dec.rad + ) * np.sin(src_crd.dec.rad) * np.cos(rd_grid[:, :, 0] - src_crd.ra.rad) + return lm_grid @@ -257,12 +265,11 @@ def uncorrupted(lm, baselines, wave, time, src_crd, array_layout, I): Returns ------- 4d array - Returns visibility for every lm and baseline + Returns visibility for every lm and baseline """ stat_num = array_layout.st_num.shape[0] base_num = int(stat_num * (stat_num - 1) / 2) - vectorized_num = np.vectorize(lambda st: st.st_num, otypes=[int]) st1, st2 = get_valid_baselines(baselines, base_num) @@ -275,12 +282,12 @@ def uncorrupted(lm, baselines, wave, time, src_crd, array_layout, I): B = np.zeros((lm.shape[0], lm.shape[1], 2, 2), dtype=complex) - B[:,:,0,0] = I[:,:,0]+I[:,:,1] - B[:,:,0,1] = I[:,:,2]+1j*I[:,:,3] - B[:,:,1,0] = I[:,:,2]-1j*I[:,:,3] - B[:,:,1,1] = I[:,:,0]-I[:,:,1] + B[:, :, 0, 0] = I[:, :, 0] + I[:, :, 1] + B[:, :, 0, 1] = I[:, :, 2] + 1j * I[:, :, 3] + B[:, :, 1, 0] = I[:, :, 2] - 1j * I[:, :, 3] + B[:, :, 1, 1] = I[:, :, 0] - I[:, :, 1] - X = torch.einsum('lmij,lmb->lmbij', torch.tensor(B), K) + X = torch.einsum("lmij,lmb->lmbij", torch.tensor(B), K) return X @@ -310,12 +317,11 @@ def corrupted(lm, baselines, wave, time, src_crd, array_layout, I, rd): Returns ------- 4d array - Returns visibility for every lm and baseline + Returns visibility for every lm and baseline """ stat_num = array_layout.st_num.shape[0] base_num = int(stat_num * (stat_num - 1) / 2) - vectorized_num = np.vectorize(lambda st: st.st_num, otypes=[int]) st1, st2 = get_valid_baselines(baselines, base_num) st1_num = vectorized_num(st1) @@ -327,13 +333,13 @@ def corrupted(lm, baselines, wave, time, src_crd, array_layout, I, rd): B = np.zeros((lm.shape[0], lm.shape[1], 2, 2), dtype=complex) - B[:,:,0,0] = I[:,:,0]+I[:,:,1] - B[:,:,0,1] = I[:,:,2]+1j*I[:,:,3] - B[:,:,1,0] = I[:,:,2]-1j*I[:,:,3] - B[:,:,1,1] = I[:,:,0]-I[:,:,1] + B[:, :, 0, 0] = I[:, :, 0] + I[:, :, 1] + B[:, :, 0, 1] = I[:, :, 2] + 1j * I[:, :, 3] + B[:, :, 1, 0] = I[:, :, 2] - 1j * I[:, :, 3] + B[:, :, 1, 1] = I[:, :, 0] - I[:, :, 1] # coherency - X = torch.einsum('lmij,lmb->lmbij', torch.tensor(B), K) + X = torch.einsum("lmij,lmb->lmbij", torch.tensor(B), K) del K # telescope response @@ -344,10 +350,10 @@ def corrupted(lm, baselines, wave, time, src_crd, array_layout, I, rd): E2 = torch.tensor(E_st[:, :, st2_num], dtype=torch.cdouble) # EX = torch.einsum('lmbij,lmbjk->lmbik',E1,X) - EX = torch.einsum('lmb,lmbij->lmbij',E1,X) + EX = torch.einsum("lmb,lmbij->lmbij", E1, X) del E1, X # EXE = torch.einsum('lmbij,lmbjk->lmbik',EX,torch.transpose(torch.conj(E2),3,4)) - EXE = torch.einsum('lmbij,lmb->lmbij',EX,torch.conj(E2)) + EXE = torch.einsum("lmbij,lmb->lmbij", EX, torch.conj(E2)) del EX, E2 # P matrix # parallactic angle @@ -362,12 +368,14 @@ def corrupted(lm, baselines, wave, time, src_crd, array_layout, I, rd): tsob = time_step_of_baseline(baselines, base_num) b1 = np.array([beta[st1_num[i], tsob[i]] for i in range(st1_num.shape[0])]) b2 = np.array([beta[st2_num[i], tsob[i]] for i in range(st2_num.shape[0])]) - P1 = torch.tensor(getP(b1),dtype=torch.cdouble) - P2 = torch.tensor(getP(b2),dtype=torch.cdouble) + P1 = torch.tensor(getP(b1), dtype=torch.cdouble) + P2 = torch.tensor(getP(b2), dtype=torch.cdouble) - PEXE = torch.einsum('bij,lmbjk->lmbik',P1,EXE) + PEXE = torch.einsum("bij,lmbjk->lmbik", P1, EXE) del EXE - PEXEP = torch.einsum('lmbij,bjk->lmbik',PEXE,torch.transpose(torch.conj(P2),1,2)) + PEXEP = torch.einsum( + "lmbij,bjk->lmbik", PEXE, torch.transpose(torch.conj(P2), 1, 2) + ) del PEXE return PEXEP @@ -389,19 +397,19 @@ def integrate(X1, X2): """ X_f = torch.stack((X1, X2)) - int_m = torch.sum(X_f,dim=2) + int_m = torch.sum(X_f, dim=2) del X_f - int_l = torch.sum(int_m,dim=1) + int_l = torch.sum(int_m, dim=1) del int_m - int_f = 0.5*torch.sum(int_l, dim=0) + int_f = 0.5 * torch.sum(int_l, dim=0) del int_l X_t = torch.stack(torch.split(int_f, int(int_f.shape[0] / 2), dim=0)) del int_f - int_t = 0.5*torch.sum(X_t, dim=0) + int_t = 0.5 * torch.sum(X_t, dim=0) del X_t - return int_t + return int_t def getE(rd, array_layout, wave, src_crd): @@ -429,11 +437,11 @@ def getE(rd, array_layout, wave, src_crd): # get diameters of all stations and do vectorizing stuff diameters = array_layout.diam - theta = angularDistance(rd,src_crd) + theta = angularDistance(rd, src_crd) - x = 2*np.pi/wave * np.einsum('s,rd->rds',diameters,theta) + x = 2 * np.pi / wave * np.einsum("s,rd->rds", diameters, theta) - E[:,:,:] = jinc(x) + E[:, :, :] = jinc(x) # E[:,:,:,0,0] = jinc(x) # E[..., 1, 1] = E[..., 0, 0] @@ -455,10 +463,10 @@ def angularDistance(rd, src_crd): 2d array Returns angular Distance for every pixel in rd grid with respect to source position """ - r = rd[:,:,0] - src_crd.ra.rad - d = rd[:,:,1] - src_crd.dec.rad + r = rd[:, :, 0] - src_crd.ra.rad + d = rd[:, :, 1] - src_crd.dec.rad - theta = np.arcsin(np.sqrt(r**2+d**2)) + theta = np.arcsin(np.sqrt(r ** 2 + d ** 2)) return theta @@ -526,11 +534,11 @@ def getK(baselines, lm, wave, base_num): l = torch.tensor(lm[:, :, 0]) m = torch.tensor(lm[:, :, 1]) - n = torch.sqrt(1-l**2-m**2) + n = torch.sqrt(1 - l ** 2 - m ** 2) ul = torch.einsum("b,ij->ijb", torch.tensor(u_cmplt), l) vm = torch.einsum("b,ij->ijb", torch.tensor(v_cmplt), m) - wn = torch.einsum('b,ij->ijb', torch.tensor(w_cmplt), (n-1)) + wn = torch.einsum("b,ij->ijb", torch.tensor(w_cmplt), (n - 1)) K = torch.exp(-2 * np.pi * 1j * (ul + vm + wn)) diff --git a/vipy/simulation/utils.py b/vipy/simulation/utils.py index c12338c..d27fac9 100644 --- a/vipy/simulation/utils.py +++ b/vipy/simulation/utils.py @@ -24,7 +24,7 @@ def read_config(conf): sim_conf["src_coord"] = SkyCoord( ra=config["sampling_options"]["fov_center_ra"], dec=config["sampling_options"]["fov_center_dec"], - unit=(u.hourangle, u.deg), + unit=(u.deg, u.deg), ) sim_conf["fov_size"] = config["sampling_options"]["fov_size"] sim_conf["corr_int_time"] = config["sampling_options"]["corr_int_time"] From fc84e526d7167d578fa82b84ba0f13f07903d5cb Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Wed, 1 Dec 2021 14:38:35 +0100 Subject: [PATCH 02/63] bump version number --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 049f5f9..33f274a 100644 --- a/setup.py +++ b/setup.py @@ -2,7 +2,7 @@ setup( name="vipy", - version="0.0.3", + version="0.0.4", description="Simulate radio interferometer observations and visibility generation.", url="https://github.com/radionets-project/vipy", author="Kevin Schmidt, Felix Geyer, Stefan Fröse", From 8c5d32b81347346d01a1c5a36cbba91f729a3b98 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Wed, 1 Dec 2021 14:44:17 +0100 Subject: [PATCH 03/63] rename to pyvisgen --- .github/workflows/ci.yml | 2 +- README.md | 2 +- examples/baselines.ipynb | 10 +++++----- examples/fits_tests.ipynb | 4 ++-- examples/visibility.ipynb | 14 +++++++------- {vipy => pyvisgen}/__init__.py | 0 {vipy => pyvisgen}/fits/writer.py | 2 +- {vipy => pyvisgen}/layouts/eht.txt | 0 {vipy => pyvisgen}/layouts/layouts.py | 0 {vipy => pyvisgen}/layouts/vlba.txt | 0 {vipy => pyvisgen}/layouts/vlba_light.txt | 0 {vipy => pyvisgen}/simulation/scan.py | 4 ++-- {vipy => pyvisgen}/simulation/utils.py | 0 {vipy => pyvisgen}/simulation/uv_plots.py | 0 scripts/cluster.py | 10 +++++----- setup.py | 4 ++-- tests/test_layouts.py | 2 +- tests/test_utils.py | 12 ++++++------ 18 files changed, 33 insertions(+), 33 deletions(-) rename {vipy => pyvisgen}/__init__.py (100%) rename {vipy => pyvisgen}/fits/writer.py (99%) rename {vipy => pyvisgen}/layouts/eht.txt (100%) rename {vipy => pyvisgen}/layouts/layouts.py (100%) rename {vipy => pyvisgen}/layouts/vlba.txt (100%) rename {vipy => pyvisgen}/layouts/vlba_light.txt (100%) rename {vipy => pyvisgen}/simulation/scan.py (99%) rename {vipy => pyvisgen}/simulation/utils.py (100%) rename {vipy => pyvisgen}/simulation/uv_plots.py (100%) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index cddbf95..cd4e4d3 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -25,6 +25,6 @@ jobs: pip install . - name: tests - run: python -m pytest --cov vipy --cov-report=xml + run: python -m pytest --cov pyvisgen --cov-report=xml - uses: codecov/codecov-action@v1 diff --git a/README.md b/README.md index cd3506a..87df06d 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# vipy [![Actions Status](https://github.com/radionets-project/vipy/workflows/CI/badge.svg)](https://github.com/radionets-project/vipy/actions) [![codecov](https://codecov.io/gh/radionets-project/vipy/branch/main/graph/badge.svg?token=84ATCQQAMN)](https://codecov.io/gh/radionets-project/vipy) +# pyvisgen [![Actions Status](https://github.com/radionets-project/pyvisgen/workflows/CI/badge.svg)](https://github.com/radionets-project/pyvisgen/actions) [![codecov](https://codecov.io/gh/radionets-project/pyvisgen/branch/main/graph/badge.svg?token=84ATCQQAMN)](https://codecov.io/gh/radionets-project/pyvisgen) Possible names: * pyvisgen diff --git a/examples/baselines.ipynb b/examples/baselines.ipynb index 2ebcf4b..e9ede4e 100644 --- a/examples/baselines.ipynb +++ b/examples/baselines.ipynb @@ -6,7 +6,7 @@ "metadata": {}, "outputs": [], "source": [ - "from vipy.simulation.utils import read_config" + "from pyvisgen.simulation.utils import read_config" ] }, { @@ -61,7 +61,7 @@ "metadata": {}, "outputs": [], "source": [ - "from vipy.layouts.layouts import get_array_layout" + "from pyvisgen.layouts.layouts import get_array_layout" ] }, { @@ -101,8 +101,8 @@ "metadata": {}, "outputs": [], "source": [ - "from vipy.simulation.scan import get_baselines\n", - "from vipy.simulation.utils import calc_time_steps\n", + "from pyvisgen.simulation.scan import get_baselines\n", + "from pyvisgen.simulation.utils import calc_time_steps\n", "import matplotlib.pyplot as plt" ] }, @@ -215,7 +215,7 @@ "metadata": {}, "outputs": [], "source": [ - "from vipy.simulation.scan import create_bgrid" + "from pyvisgen.simulation.scan import create_bgrid" ] }, { diff --git a/examples/fits_tests.ipynb b/examples/fits_tests.ipynb index f9c3eb2..97545db 100644 --- a/examples/fits_tests.ipynb +++ b/examples/fits_tests.ipynb @@ -188,7 +188,7 @@ "outputs": [], "source": [ "from astropy.time import Time\n", - "from vipy.simulation.utils import read_config" + "from pyvisgen.simulation.utils import read_config" ] }, { @@ -547,7 +547,7 @@ "metadata": {}, "outputs": [], "source": [ - "def create_hdu_list(conf, path=\"../vipy/layouts/eht.txt\"):\n", + "def create_hdu_list(conf, path=\"../pyvisgen/layouts/eht.txt\"):\n", " vis_hdu = create_vis_hdu(None, conf)\n", " time_hdu = create_time_hdu(None)\n", " freq_hdu = create_frequency_hdu(conf)\n", diff --git a/examples/visibility.ipynb b/examples/visibility.ipynb index a16133b..9d4e5b1 100644 --- a/examples/visibility.ipynb +++ b/examples/visibility.ipynb @@ -23,8 +23,8 @@ "metadata": {}, "outputs": [], "source": [ - "import vipy.simulation.utils as ut\n", - "import vipy.layouts.layouts as layouts\n", + "import pyvisgen.simulation.utils as ut\n", + "import pyvisgen.layouts.layouts as layouts\n", "import astropy.constants as const\n", "from astropy import units as un\n", "import time as t\n", @@ -41,7 +41,7 @@ "metadata": {}, "outputs": [], "source": [ - "# rc = ut.read_config('/net/nfshome/home/sfroese/vipy/config/default.toml')\n", + "# rc = ut.read_config('/net/nfshome/home/sfroese/pyvisgen/config/default.toml')\n", "rc = ut.read_config(\"../config/default.toml\")\n", "array_layout = layouts.get_array_layout('eht')\n", "src_crd = rc['src_coord']\n", @@ -65,7 +65,7 @@ "metadata": {}, "outputs": [], "source": [ - "import vipy.simulation.scan as scan\n", + "import pyvisgen.simulation.scan as scan\n", "import matplotlib.pyplot as plt" ] }, @@ -457,7 +457,7 @@ } ], "source": [ - "from vipy.simulation.scan import get_valid_baselines\n", + "from pyvisgen.simulation.scan import get_valid_baselines\n", "import astropy.units as un\n", "from astropy.time import Time\n", "visibilities = Visibilities([], [], [], [], [], [], [], [], [], [], [], [])\n", @@ -545,8 +545,8 @@ "metadata": {}, "outputs": [], "source": [ - "import vipy.simulation.utils as ut\n", - "import vipy.fits.writer as writer" + "import pyvisgen.simulation.utils as ut\n", + "import pyvisgen.fits.writer as writer" ] }, { diff --git a/vipy/__init__.py b/pyvisgen/__init__.py similarity index 100% rename from vipy/__init__.py rename to pyvisgen/__init__.py diff --git a/vipy/fits/writer.py b/pyvisgen/fits/writer.py similarity index 99% rename from vipy/fits/writer.py rename to pyvisgen/fits/writer.py index 4461edc..6b0c684 100644 --- a/vipy/fits/writer.py +++ b/pyvisgen/fits/writer.py @@ -373,7 +373,7 @@ def create_antenna_hdu(layout_txt, conf, layout="vlba"): return hdu_ant -def create_hdu_list(data, conf, path="../vipy/vipy/layouts/vlba.txt"): +def create_hdu_list(data, conf, path="../pyvisgen/pyvisgen/layouts/vlba.txt"): vis_hdu = create_vis_hdu(data, conf) time_hdu = create_time_hdu(data) freq_hdu = create_frequency_hdu(conf) diff --git a/vipy/layouts/eht.txt b/pyvisgen/layouts/eht.txt similarity index 100% rename from vipy/layouts/eht.txt rename to pyvisgen/layouts/eht.txt diff --git a/vipy/layouts/layouts.py b/pyvisgen/layouts/layouts.py similarity index 100% rename from vipy/layouts/layouts.py rename to pyvisgen/layouts/layouts.py diff --git a/vipy/layouts/vlba.txt b/pyvisgen/layouts/vlba.txt similarity index 100% rename from vipy/layouts/vlba.txt rename to pyvisgen/layouts/vlba.txt diff --git a/vipy/layouts/vlba_light.txt b/pyvisgen/layouts/vlba_light.txt similarity index 100% rename from vipy/layouts/vlba_light.txt rename to pyvisgen/layouts/vlba_light.txt diff --git a/vipy/simulation/scan.py b/pyvisgen/simulation/scan.py similarity index 99% rename from vipy/simulation/scan.py rename to pyvisgen/simulation/scan.py index 2e79749..4815fa1 100644 --- a/vipy/simulation/scan.py +++ b/pyvisgen/simulation/scan.py @@ -6,8 +6,8 @@ import scipy.constants as const import scipy.signal as sig from astroplan import Observer -from vipy.simulation.utils import single_occurance, get_pairs -from vipy.layouts import layouts +from pyvisgen.simulation.utils import single_occurance, get_pairs +from pyvisgen.layouts import layouts import torch import itertools diff --git a/vipy/simulation/utils.py b/pyvisgen/simulation/utils.py similarity index 100% rename from vipy/simulation/utils.py rename to pyvisgen/simulation/utils.py diff --git a/vipy/simulation/uv_plots.py b/pyvisgen/simulation/uv_plots.py similarity index 100% rename from vipy/simulation/uv_plots.py rename to pyvisgen/simulation/uv_plots.py diff --git a/scripts/cluster.py b/scripts/cluster.py index 225efb0..c57fbaa 100644 --- a/scripts/cluster.py +++ b/scripts/cluster.py @@ -1,15 +1,15 @@ -import vipy.simulation.utils as ut -import vipy.layouts.layouts as layouts +import pyvisgen.simulation.utils as ut +import pyvisgen.layouts.layouts as layouts import astropy.constants as const from astropy import units as un import time as t import numpy as np -import vipy.simulation.scan as scan +import pyvisgen.simulation.scan as scan import torch from astropy.io import fits -from vipy.simulation.scan import get_valid_baselines +from pyvisgen.simulation.scan import get_valid_baselines from astropy.time import Time -import vipy.fits.writer as writer +import pyvisgen.fits.writer as writer from dataclasses import dataclass import sys diff --git a/setup.py b/setup.py index 33f274a..a094a11 100644 --- a/setup.py +++ b/setup.py @@ -1,10 +1,10 @@ from setuptools import setup, find_packages setup( - name="vipy", + name="pyvisgen", version="0.0.4", description="Simulate radio interferometer observations and visibility generation.", - url="https://github.com/radionets-project/vipy", + url="https://github.com/radionets-project/pyvisgen", author="Kevin Schmidt, Felix Geyer, Stefan Fröse", author_email="kevin3.schmidt@tu-dortmund.de", license="MIT", diff --git a/tests/test_layouts.py b/tests/test_layouts.py index fe05876..41b69e0 100644 --- a/tests/test_layouts.py +++ b/tests/test_layouts.py @@ -2,7 +2,7 @@ def test_get_array_layout(): - from vipy.layouts.layouts import get_array_layout + from pyvisgen.layouts.layouts import get_array_layout layout = get_array_layout("eht") diff --git a/tests/test_utils.py b/tests/test_utils.py index 9f109b7..5305a81 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -2,7 +2,7 @@ def test_read_config(): - from vipy.simulation.utils import read_config + from pyvisgen.simulation.utils import read_config conf = read_config("config/default.toml") @@ -20,7 +20,7 @@ def test_read_config(): def test_single_occurance(): - from vipy.simulation.utils import single_occurance + from pyvisgen.simulation.utils import single_occurance arr = np.array([1, 2, 2, 3, 4, 4]) @@ -30,8 +30,8 @@ def test_single_occurance(): def test_get_pairs(): - from vipy.layouts.layouts import get_array_layout - from vipy.simulation.utils import get_pairs + from pyvisgen.layouts.layouts import get_array_layout + from pyvisgen.simulation.utils import get_pairs layout = get_array_layout("eht") delta_x, delta_y, delta_z = get_pairs(layout) @@ -41,8 +41,8 @@ def test_get_pairs(): def test_calc_time_steps(): - from vipy.simulation.utils import read_config - from vipy.simulation.utils import calc_time_steps + from pyvisgen.simulation.utils import read_config + from pyvisgen.simulation.utils import calc_time_steps conf = read_config("config/default.toml") time = calc_time_steps(conf) From dc6049bb29bb1a0ff049cc1a295b8dab23536e21 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Wed, 1 Dec 2021 15:19:47 +0100 Subject: [PATCH 04/63] black --- pyvisgen/simulation/scan.py | 68 +++++++++++-------------------------- 1 file changed, 20 insertions(+), 48 deletions(-) diff --git a/pyvisgen/simulation/scan.py b/pyvisgen/simulation/scan.py index 3647d45..91060d5 100644 --- a/pyvisgen/simulation/scan.py +++ b/pyvisgen/simulation/scan.py @@ -11,9 +11,10 @@ import torch import itertools import time as t -import numexpr as ne # fast exponential +import numexpr as ne # fast exponential from einsumt import einsumt as einsum + @dataclass class Baselines: name: [str] @@ -321,7 +322,7 @@ def corrupted(lm, baselines, wave, time, src_crd, array_layout, I, rd): 4d array Returns visibility for every lm and baseline """ - + stat_num = array_layout.st_num.shape[0] base_num = int(stat_num * (stat_num - 1) / 2) @@ -332,11 +333,8 @@ def corrupted(lm, baselines, wave, time, src_crd, array_layout, I, rd): if st1_num.shape[0] == 0: return torch.zeros(1) - K = getK(baselines, lm, wave, base_num) - - B = np.zeros((lm.shape[0], lm.shape[1], 2, 2), dtype=complex) B[:, :, 0, 0] = I[:, :, 0] + I[:, :, 1] @@ -345,15 +343,10 @@ def corrupted(lm, baselines, wave, time, src_crd, array_layout, I, rd): B[:, :, 1, 1] = I[:, :, 0] - I[:, :, 1] # coherency -<<<<<<< HEAD:pyvisgen/simulation/scan.py X = torch.einsum("lmij,lmb->lmbij", torch.tensor(B), K) -======= - X = torch.einsum('lmij,lmb->lmbij', torch.tensor(B), K) # X = np.einsum('lmij,lmb->lmbij', B, K, optimize=True) # X = torch.tensor(B)[:,:,None,:,:] * K[:,:,:,None,None] - ->>>>>>> origin/main:vipy/simulation/scan.py del K # telescope response @@ -363,20 +356,11 @@ def corrupted(lm, baselines, wave, time, src_crd, array_layout, I, rd): E1 = torch.tensor(E_st[:, :, st1_num], dtype=torch.cdouble) E2 = torch.tensor(E_st[:, :, st2_num], dtype=torch.cdouble) -<<<<<<< HEAD:pyvisgen/simulation/scan.py - # EX = torch.einsum('lmbij,lmbjk->lmbik',E1,X) EX = torch.einsum("lmb,lmbij->lmbij", E1, X) - del E1, X - # EXE = torch.einsum('lmbij,lmbjk->lmbik',EX,torch.transpose(torch.conj(E2),3,4)) - EXE = torch.einsum("lmbij,lmb->lmbij", EX, torch.conj(E2)) -======= - - EX = torch.einsum('lmb,lmbij->lmbij',E1,X) del E1, X # EXE = torch.einsum('lmbij,lmbjk->lmbik',EX,torch.transpose(torch.conj(E2),3,4)) - EXE = torch.einsum('lmbij,lmb->lmbij',EX,E2) ->>>>>>> origin/main:vipy/simulation/scan.py + EXE = torch.einsum("lmbij,lmb->lmbij", EX, E2) del EX, E2 # P matrix @@ -396,13 +380,7 @@ def corrupted(lm, baselines, wave, time, src_crd, array_layout, I, rd): P1 = torch.tensor(getP(b1), dtype=torch.cdouble) P2 = torch.tensor(getP(b2), dtype=torch.cdouble) -<<<<<<< HEAD:pyvisgen/simulation/scan.py PEXE = torch.einsum("bij,lmbjk->lmbik", P1, EXE) -======= - - - PEXE = torch.einsum('bij,lmbjk->lmbik',P1,EXE) ->>>>>>> origin/main:vipy/simulation/scan.py del EXE PEXEP = torch.einsum( "lmbij,bjk->lmbik", PEXE, torch.transpose(torch.conj(P2), 1, 2) @@ -411,6 +389,7 @@ def corrupted(lm, baselines, wave, time, src_crd, array_layout, I, rd): return PEXEP + def direction_independent(lm, baselines, wave, time, src_crd, array_layout, I, rd): """Calculates direction independet visibility @@ -436,13 +415,11 @@ def direction_independent(lm, baselines, wave, time, src_crd, array_layout, I, r Returns ------- 4d array - Returns visibility for every lm and baseline + Returns visibility for every lm and baseline """ - stat_num = array_layout.st_num.shape[0] base_num = int(stat_num * (stat_num - 1) / 2) - vectorized_num = np.vectorize(lambda st: st.st_num, otypes=[int]) st1, st2 = get_valid_baselines(baselines, base_num) st1_num = vectorized_num(st1) @@ -450,23 +427,18 @@ def direction_independent(lm, baselines, wave, time, src_crd, array_layout, I, r if st1_num.shape[0] == 0: return torch.zeros(1) - K = getK(baselines, lm, wave, base_num) - - B = np.zeros((lm.shape[0], lm.shape[1], 2, 2), dtype=complex) - B[:,:,0,0] = I[:,:,0]+I[:,:,1] - B[:,:,0,1] = I[:,:,2]+1j*I[:,:,3] - B[:,:,1,0] = I[:,:,2]-1j*I[:,:,3] - B[:,:,1,1] = I[:,:,0]-I[:,:,1] + B[:, :, 0, 0] = I[:, :, 0] + I[:, :, 1] + B[:, :, 0, 1] = I[:, :, 2] + 1j * I[:, :, 3] + B[:, :, 1, 0] = I[:, :, 2] - 1j * I[:, :, 3] + B[:, :, 1, 1] = I[:, :, 0] - I[:, :, 1] # coherency - X = torch.einsum('lmij,lmb->lmbij', torch.tensor(B), K) - + X = torch.einsum("lmij,lmb->lmbij", torch.tensor(B), K) - del K # telescope response @@ -475,16 +447,16 @@ def direction_independent(lm, baselines, wave, time, src_crd, array_layout, I, r E1 = torch.tensor(E_st[:, :, st1_num], dtype=torch.cdouble) E2 = torch.tensor(E_st[:, :, st2_num], dtype=torch.cdouble) - - EX = torch.einsum('lmb,lmbij->lmbij',E1,X) + EX = torch.einsum("lmb,lmbij->lmbij", E1, X) del E1, X - EXE = torch.einsum('lmbij,lmb->lmbij',EX,E2) + EXE = torch.einsum("lmbij,lmb->lmbij", EX, E2) del EX, E2 return EXE + def integrate(X1, X2): """Summation over l and m and avering over time and freq @@ -619,7 +591,7 @@ def getK(baselines, lm, wave, base_num): Shape is given by lm axes and baseline axis """ # new valid baseline calculus. for details see function get_valid_baselines() - + valid = baselines.valid.reshape(-1, base_num) mask = np.array(valid[:-1]).astype(bool) & np.array(valid[1:]).astype(bool) @@ -640,15 +612,15 @@ def getK(baselines, lm, wave, base_num): l = torch.tensor(lm[:, :, 0]) m = torch.tensor(lm[:, :, 1]) - n = torch.sqrt(1-l**2-m**2) - + n = torch.sqrt(1 - l ** 2 - m ** 2) + ul = torch.einsum("b,ij->ijb", torch.tensor(u_cmplt), l) vm = torch.einsum("b,ij->ijb", torch.tensor(v_cmplt), m) - wn = torch.einsum("b,ij->ijb", torch.tensor(w_cmplt), (n-1)) - + wn = torch.einsum("b,ij->ijb", torch.tensor(w_cmplt), (n - 1)) + pi = np.pi test = ul + vm + wn - K = ne.evaluate("exp(-2 * pi * 1j * (ul + vm + wn))") #-0.4 secs for vlba + K = ne.evaluate("exp(-2 * pi * 1j * (ul + vm + wn))") # -0.4 secs for vlba return torch.tensor(K) From e469aee052976248f2afcf77aae704be33ecd460 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Thu, 2 Dec 2021 17:10:19 +0100 Subject: [PATCH 05/63] black --- scripts/cluster.py | 86 ++++++++++++++++++++++++++++++---------------- 1 file changed, 57 insertions(+), 29 deletions(-) diff --git a/scripts/cluster.py b/scripts/cluster.py index c57fbaa..044e485 100644 --- a/scripts/cluster.py +++ b/scripts/cluster.py @@ -15,8 +15,10 @@ import sys from PIL import Image + # jobid = sys.getenv('SLURM_ARRAY_TASK_ID') + @dataclass class Visibilities: I: [complex] @@ -51,7 +53,7 @@ def __getitem__(self, i): def get_values(self): return np.array([self.I, self.Q, self.U, self.V]) - + def add(self, visibilities): self.I = np.concatenate([self.I, visibilities.I]) self.Q = np.concatenate([self.Q, visibilities.Q]) @@ -65,7 +67,7 @@ def add(self, visibilities): self.w = np.concatenate([self.w, visibilities.w]) self.date = np.concatenate([self.date, visibilities.date]) self._date = np.concatenate([self._date, visibilities._date]) - + @dataclass class Vis: @@ -85,32 +87,55 @@ class Vis: # set number of threads + def loop(): torch.set_num_threads(48) # read config rc = ut.read_config("../config/vlba.toml") - array_layout = layouts.get_array_layout('vlba') - src_crd = rc['src_coord'] - - wave1 = const.c/((float(rc['channel'].split(':')[0])-float(rc['channel'].split(':')[1])/2)*10**6/un.second)/un.meter - wave2 = const.c/((float(rc['channel'].split(':')[0])+float(rc['channel'].split(':')[1])/2)*10**6/un.second)/un.meter + array_layout = layouts.get_array_layout("vlba") + src_crd = rc["src_coord"] + + wave1 = ( + const.c + / ( + ( + float(rc["channel"].split(":")[0]) + - float(rc["channel"].split(":")[1]) / 2 + ) + * 10 ** 6 + / un.second + ) + / un.meter + ) + wave2 = ( + const.c + / ( + ( + float(rc["channel"].split(":")[0]) + + float(rc["channel"].split(":")[1]) / 2 + ) + * 10 ** 6 + / un.second + ) + / un.meter + ) # calculate rd, lm - rd = scan.rd_grid(rc['fov_size']*np.pi/(3600*180),256, src_crd) + rd = scan.rd_grid(rc["fov_size"] * np.pi / (3600 * 180), 256, src_crd) lm = scan.lm_grid(rd, src_crd) # calculate time steps time = ut.calc_time_steps(rc) # open image - img = np.asarray(Image.open('/net/nfshome/home/sfroese/150.jpg')) + img = np.asarray(Image.open("/net/nfshome/home/sfroese/150.jpg")) norm = np.sum(img) - img = img/norm + img = img / norm img = torch.tensor(img) - I = torch.zeros((img.shape[0],img.shape[1],4), dtype=torch.cdouble) - I[...,0] = img - I[...,1] = img + I = torch.zeros((img.shape[0], img.shape[1], 4), dtype=torch.cdouble) + I[..., 0] = img + I[..., 1] = img stat_num = array_layout.st_num.shape[0] base_num = int(stat_num * (stat_num - 1) / 2) @@ -118,12 +143,12 @@ def loop(): # calculate vis visibilities = Visibilities([], [], [], [], [], [], [], [], [], [], [], []) vis_num = np.zeros(1) -#i in total number of scans + # i in total number of scans for i in range(72): # i=30 - t = time[i*31:(i+1)*31] + t = time[i * 31 : (i + 1) * 31] baselines = scan.get_baselines(src_crd, t, array_layout) - + valid = baselines.valid.reshape(-1, base_num) mask = np.array(valid[:-1]).astype(bool) & np.array(valid[1:]).astype(bool) u = baselines.u.reshape(-1, base_num) @@ -133,41 +158,44 @@ def loop(): u_valid = u[:-1][mask] v_valid = v[:-1][mask] w_valid = w[:-1][mask] - date = np.repeat((t[:-1]+rc['corr_int_time']*un.second/2).jd.reshape(-1, 1), base_num, axis=1)[mask] - _date = np.zeros(len(u_valid)) - - + date = np.repeat( + (t[:-1] + rc["corr_int_time"] * un.second / 2).jd.reshape(-1, 1), + base_num, + axis=1, + )[mask] + _date = np.zeros(len(u_valid)) + X1 = scan.corrupted(lm, baselines, wave1, time, src_crd, array_layout, I, rd) if X1.shape[0] == 1: continue X2 = scan.corrupted(lm, baselines, wave1, time, src_crd, array_layout, I, rd) - - vis_num = np.arange(X1.shape[2]//2) + 1 + vis_num.max() - + + vis_num = np.arange(X1.shape[2] // 2) + 1 + vis_num.max() int_values = scan.integrate(X1, X2) del X1, X2 - int_values = int_values.reshape(-1,4) - + int_values = int_values.reshape(-1, 4) + vis = Visibilities( int_values[:, 0], int_values[:, 1], int_values[:, 2], int_values[:, 3], vis_num, - np.repeat(i+1, len(vis_num)), + np.repeat(i + 1, len(vis_num)), np.array([baselines[i].baselineNum() for i in base_valid]), u_valid, v_valid, w_valid, date, - _date, + _date, ) - + visibilities.add(vis) # break return visibilities + # save vis hdu_list = writer.create_hdu_list(loop(), rc) -hdu_list.writeto('test07_vlba.fits', overwrite=True) \ No newline at end of file +hdu_list.writeto("test07_vlba.fits", overwrite=True) From fc5223e9bb90976c38ea1c0d02d3472dc8541286 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Thu, 2 Dec 2021 17:10:46 +0100 Subject: [PATCH 06/63] sampling mask simulation with multi channel --- pyvisgen/simulation/sampling_masks.py | 76 +++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 pyvisgen/simulation/sampling_masks.py diff --git a/pyvisgen/simulation/sampling_masks.py b/pyvisgen/simulation/sampling_masks.py new file mode 100644 index 0000000..1b3c830 --- /dev/null +++ b/pyvisgen/simulation/sampling_masks.py @@ -0,0 +1,76 @@ +import re +import h5py +import numpy as np +from pathlib import Path +import matplotlib.pyplot as plt +from pyvisgen.simulation.scan import get_baselines +from pyvisgen.layouts.layouts import get_array_layout +from pyvisgen.simulation.utils import read_config, calc_time_steps + + +def create_sampling_mask( + layout="vlba", + size=256, + multi_channel=True, + bandwidths=4, + base_freq=15.21e10, + frequsel=[0, 8e7, 1.44e8, 2.08e8], +): + conf = read_config("../../config/vlba.toml") + time = calc_time_steps(conf) + layout = get_array_layout(layout) + baselines = get_baselines(conf["src_coord"], time, layout) + u = np.concatenate( + (baselines.u[baselines.valid == True], -baselines.u[baselines.valid == True]) + ) + v = np.concatenate( + (baselines.v[baselines.valid == True], -baselines.v[baselines.valid == True]) + ) + + if multi_channel: + u = np.repeat(u[None], bandwidths, axis=0) + v = np.repeat(v[None], bandwidths, axis=0) + scales = np.array(frequsel) + base_freq + u /= scales[:, None] + v /= scales[:, None] + else: + u /= base_freq + v /= base_freq + u = u[None] + v = v[None] + + uv_hist, _, _ = np.histogram2d( + u.ravel(), + v.ravel(), + bins=size, + ) + mask = uv_hist > 0 + return mask + + +def open_sky_distributions(path): + with h5py.File(path, "r") as hf: + images = np.array([hf["sky" + str(i)] for i in range(int(len(hf) / 3))]) + return images + + +def get_data_paths(path): + data_dir = Path(path) + bundles = np.array([x for x in data_dir.iterdir()]) + data_paths = np.sort([path for path in bundles if re.findall("source_", path.name)]) + return data_paths + + +def sampling(path="../../../radiosim/build/example_data"): + data_paths = get_data_paths(path) + for data in data_paths: + images = open_sky_distributions(data) + mask = np.array([create_sampling_mask() for i in range(len(images))]) + sampled = images.copy() + sampled[~mask.astype(bool)] = 0 + plt.imshow(mask[0].astype(bool)) + plt.show() + + +if __name__ == "__main__": + sampling() From cb37098e9aa44337dd8b5c6a2e558f49cfcf98b4 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Thu, 2 Dec 2021 18:12:56 +0100 Subject: [PATCH 07/63] add sky dist data set --- pyvisgen/utils/sky_distributions.py | 50 +++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 pyvisgen/utils/sky_distributions.py diff --git a/pyvisgen/utils/sky_distributions.py b/pyvisgen/utils/sky_distributions.py new file mode 100644 index 0000000..2951d03 --- /dev/null +++ b/pyvisgen/utils/sky_distributions.py @@ -0,0 +1,50 @@ +import h5py +import numpy as np + + +class h5_sky_distributions: + def __init__(self, bundle_paths): + """ + Save the bundle paths and the number of bundles in one file. + """ + self.bundles = bundle_paths + self.num_img = len(self.open_bundle(self.bundles[0])) // 3 + + def __call__(self): + return print("This is the h5_sky_distributions class.") + + def __len__(self): + """ + Returns the total number of images in this dataset. + """ + return len(self.bundles) * self.num_img + + def __getitem__(self, i): + sky = self.open_image("sky", i) + return sky + + def open_bundle(self, bundle_path): + bundle = h5py.File(bundle_path, "r") + return bundle + + def open_image(self, var, i): + if isinstance(i, int): + i = np.array([i]) + indices = np.sort(i) + bundle = indices // self.num_img + image = indices - bundle * self.num_img + bundle_unique = np.unique(bundle) + bundle_files = [ + h5py.File(self.bundles[bundle], "r") for bundle in bundle_unique + ] + bundle_files_str = list(map(str, bundle_files)) + data = np.array( + [ + bund["sky" + str(img)] + for bund, bund_str in zip(bundle_files, bundle_files_str) + for img in image[ + bundle == bundle_unique[bundle_files_str.index(bund_str)] + ] + ] + ) + return data From 15fad6b1444cc4d9eae7f0074001ab1995a0015f Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Thu, 2 Dec 2021 18:25:10 +0100 Subject: [PATCH 08/63] rename --- pyvisgen/utils/{sky_distributions.py => data.py} | 9 +++++++++ 1 file changed, 9 insertions(+) rename pyvisgen/utils/{sky_distributions.py => data.py} (85%) diff --git a/pyvisgen/utils/sky_distributions.py b/pyvisgen/utils/data.py similarity index 85% rename from pyvisgen/utils/sky_distributions.py rename to pyvisgen/utils/data.py index 2951d03..431f664 100644 --- a/pyvisgen/utils/sky_distributions.py +++ b/pyvisgen/utils/data.py @@ -1,5 +1,14 @@ +import re import h5py import numpy as np +from pathlib import Path + + +def get_data_paths(path): + data_dir = Path(path) + bundles = np.array([x for x in data_dir.iterdir()]) + data_paths = np.sort([path for path in bundles if re.findall("source_", path.name)]) + return data_paths class h5_sky_distributions: From 70caa05f7d5a742d61aed29f641ae447b1f5fab8 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Thu, 2 Dec 2021 18:25:34 +0100 Subject: [PATCH 09/63] keep only basic sampling --- pyvisgen/simulation/sampling_masks.py | 55 +++++++++------------------ 1 file changed, 17 insertions(+), 38 deletions(-) diff --git a/pyvisgen/simulation/sampling_masks.py b/pyvisgen/simulation/sampling_masks.py index 1b3c830..05408af 100644 --- a/pyvisgen/simulation/sampling_masks.py +++ b/pyvisgen/simulation/sampling_masks.py @@ -1,7 +1,4 @@ -import re -import h5py import numpy as np -from pathlib import Path import matplotlib.pyplot as plt from pyvisgen.simulation.scan import get_baselines from pyvisgen.layouts.layouts import get_array_layout @@ -11,8 +8,6 @@ def create_sampling_mask( layout="vlba", size=256, - multi_channel=True, - bandwidths=4, base_freq=15.21e10, frequsel=[0, 8e7, 1.44e8, 2.08e8], ): @@ -27,17 +22,11 @@ def create_sampling_mask( (baselines.v[baselines.valid == True], -baselines.v[baselines.valid == True]) ) - if multi_channel: - u = np.repeat(u[None], bandwidths, axis=0) - v = np.repeat(v[None], bandwidths, axis=0) - scales = np.array(frequsel) + base_freq - u /= scales[:, None] - v /= scales[:, None] - else: - u /= base_freq - v /= base_freq - u = u[None] - v = v[None] + u = np.repeat(u[None], len(frequsel), axis=0) + v = np.repeat(v[None], len(frequsel), axis=0) + scales = np.array(frequsel) + base_freq + u /= scales[:, None] + v /= scales[:, None] uv_hist, _, _ = np.histogram2d( u.ravel(), @@ -48,28 +37,18 @@ def create_sampling_mask( return mask -def open_sky_distributions(path): - with h5py.File(path, "r") as hf: - images = np.array([hf["sky" + str(i)] for i in range(int(len(hf) / 3))]) - return images - - -def get_data_paths(path): - data_dir = Path(path) - bundles = np.array([x for x in data_dir.iterdir()]) - data_paths = np.sort([path for path in bundles if re.findall("source_", path.name)]) - return data_paths - - -def sampling(path="../../../radiosim/build/example_data"): - data_paths = get_data_paths(path) - for data in data_paths: - images = open_sky_distributions(data) - mask = np.array([create_sampling_mask() for i in range(len(images))]) - sampled = images.copy() - sampled[~mask.astype(bool)] = 0 - plt.imshow(mask[0].astype(bool)) - plt.show() +def sampling(config, sky_dist): + mask = create_sampling_mask( + layout=config["layout"], + size=sky_dist.shape[-1], + base_freq=config["base_freq"], + frequsel=config["frequsel"], + ) + sampled = sky_dist.copy() + sampled[~mask.astype(bool)] = 0 + plt.imshow(mask[0].astype(bool)) + plt.show() + return sampled if __name__ == "__main__": From 353d037f20d0914168301487e3bd7a0ebdf6176b Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Wed, 8 Dec 2021 12:10:00 +0100 Subject: [PATCH 10/63] black --- pyvisgen/simulation/scan.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyvisgen/simulation/scan.py b/pyvisgen/simulation/scan.py index 91060d5..398af8f 100644 --- a/pyvisgen/simulation/scan.py +++ b/pyvisgen/simulation/scan.py @@ -538,7 +538,8 @@ def angularDistance(rd, src_crd): Returns ------- 2d array - Returns angular Distance for every pixel in rd grid with respect to source position + Returns angular Distance for every pixel in rd grid with respect + to source position """ r = rd[:, :, 0] - src_crd.ra.rad d = rd[:, :, 1] - src_crd.dec.rad From 7e179ccb4262dd548a5e923b442d48daec6c0191 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Wed, 8 Dec 2021 12:11:43 +0100 Subject: [PATCH 11/63] new time format in conf, units as un --- pyvisgen/simulation/utils.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pyvisgen/simulation/utils.py b/pyvisgen/simulation/utils.py index d27fac9..7ab4193 100644 --- a/pyvisgen/simulation/utils.py +++ b/pyvisgen/simulation/utils.py @@ -1,6 +1,6 @@ import toml from astropy.coordinates import SkyCoord -import astropy.units as u +import astropy.units as un import numpy as np from astropy.time import Time @@ -24,7 +24,7 @@ def read_config(conf): sim_conf["src_coord"] = SkyCoord( ra=config["sampling_options"]["fov_center_ra"], dec=config["sampling_options"]["fov_center_dec"], - unit=(u.deg, u.deg), + unit=(un.deg, un.deg), ) sim_conf["fov_size"] = config["sampling_options"]["fov_size"] sim_conf["corr_int_time"] = config["sampling_options"]["corr_int_time"] @@ -67,7 +67,7 @@ def get_pairs(array_layout): def calc_time_steps(conf): - start_time = Time(conf["scan_start"], format="yday") + start_time = Time(conf["scan_start"].isoformat(), format='isot') interval = conf["interval_length"] integration_time = conf["corr_int_time"] num_scans = conf["scans"] @@ -75,7 +75,7 @@ def calc_time_steps(conf): int_time = conf["corr_int_time"] time_lst = [ - start_time + interval * i * u.second + j * integration_time * u.second + start_time + interval * i * un.second + j * integration_time * un.second for i in range(num_scans) for j in range( int(scan_duration / int_time) + 1 From ac1c94dc7f943f34fe01a5c128d1ff585b55ec11 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Wed, 8 Dec 2021 12:12:38 +0100 Subject: [PATCH 12/63] adjust to new config format --- pyvisgen/fits/writer.py | 35 +++++++++++++++-------------------- 1 file changed, 15 insertions(+), 20 deletions(-) diff --git a/pyvisgen/fits/writer.py b/pyvisgen/fits/writer.py index 6b0c684..0decf78 100644 --- a/pyvisgen/fits/writer.py +++ b/pyvisgen/fits/writer.py @@ -18,10 +18,9 @@ def create_vis_hdu(data, conf, layout="vlba", source_name="sim-source-0"): data.date.min() ) # placeholder, julian date of vis, central time in the integration period - # I think this is not really needed, but dunno, documentation is again insane _DATE = ( data._date - ) # relative julian date for the observation day??, central time in the integration period + ) # central time in the integration period BASELINE = data.base_num @@ -31,7 +30,7 @@ def create_vis_hdu(data, conf, layout="vlba", source_name="sim-source-0"): values = data.get_values() vis = np.swapaxes( np.swapaxes( - np.stack([values.real, values.imag, np.ones(values.shape)], axis=1), 1, 2 + np.stack([values.real, values.imag, np.ones(values.shape)-0.8], axis=1), 1, 2 ), 0, 1, @@ -40,11 +39,11 @@ def create_vis_hdu(data, conf, layout="vlba", source_name="sim-source-0"): # in dim 4 = IFs , dim = 1, dim 4 = number of jones, 3 = real, imag, weight # wcs - ra = conf["src_coord"].ra.value - dec = conf["src_coord"].dec.value - freq, freq_d = freq, freq_d = ( - (np.array(conf["channel"].split(":")).astype("int") * un.MHz).to(un.Hz).value - ) + ra = conf["fov_center_ra"] + dec = conf["fov_center_dec"] + freq = (conf["base_freq"] * un.Hz).value + freq_d = (conf["bandwidths"][0] * un.Hz).value + ws = wcs.WCS(naxis=7) ws.wcs.crpix = [1, 1, 1, 1, 1, 1, 1] ws.wcs.cdelt = np.array([1, 1, -1, freq_d, 1, 1, 1]) @@ -91,9 +90,8 @@ def create_vis_hdu(data, conf, layout="vlba", source_name="sim-source-0"): hdu_vis.header.comments["PTYPE6"] = "Relative Julian date ?" hdu_vis.header.comments["PTYPE7"] = "Integration time" - date_obs = Time(conf["scan_start"], format="yday").to_value( - format="iso", subfmt="date" - ) + date_obs = conf["scan_start"].date().strftime("%Y-%m-%d") + date_map = Time.now().to_value(format="iso", subfmt="date") # add additional keys @@ -179,10 +177,9 @@ def create_time_hdu(data): def create_frequency_hdu(conf): - freq, freq_d = freq, freq_d = ( - (np.array(conf["channel"].split(":")).astype("int") * un.MHz).to(un.Hz).value - ) - num_ifs = 1 # at the moment only 1 possible + freq_d = (conf["bandwidths"][0] * un.Hz).value + + # num_ifs = 1 # at the moment only 1 possible FRQSEL = np.array([1], dtype=">i4") col1 = fits.Column(name="FRQSEL", format="1J", unit=" ", array=FRQSEL) @@ -303,10 +300,8 @@ def create_antenna_hdu(layout_txt, conf, layout="vlba"): ) hdu_ant = fits.BinTableHDU.from_columns(coldefs_ant) - freq, freq_d = freq, freq_d = ( - (np.array(conf["channel"].split(":")).astype("int") * un.MHz).to(un.Hz).value - ) - ref_date = Time(conf["scan_start"], format="yday") + freq = (conf["base_freq"] * un.Hz).value + ref_date = Time(conf["scan_start"].isoformat(), format="isot") # add additional keywords hdu_ant.header["EXTNAME"] = ("AIPS AN", "AIPS table file") @@ -373,7 +368,7 @@ def create_antenna_hdu(layout_txt, conf, layout="vlba"): return hdu_ant -def create_hdu_list(data, conf, path="../pyvisgen/pyvisgen/layouts/vlba.txt"): +def create_hdu_list(data, conf, path="../layouts/vlba.txt"): vis_hdu = create_vis_hdu(data, conf) time_hdu = create_time_hdu(data) freq_hdu = create_frequency_hdu(conf) From c51ae1ad7eef05b493e68fcf0863874a1c77e4a8 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Wed, 8 Dec 2021 12:13:51 +0100 Subject: [PATCH 13/63] add more build directories --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 57e7c58..474fe7a 100644 --- a/.gitignore +++ b/.gitignore @@ -11,6 +11,7 @@ # make *_done */build +*/*/build *.ipynb_checkpoints *.pyc From 3c96e46b8aea4698ef608c83b26c9fc385794834 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Wed, 8 Dec 2021 12:14:14 +0100 Subject: [PATCH 14/63] add config to create data sets --- config/data_set.toml | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) create mode 100644 config/data_set.toml diff --git a/config/data_set.toml b/config/data_set.toml new file mode 100644 index 0000000..4be20ab --- /dev/null +++ b/config/data_set.toml @@ -0,0 +1,24 @@ +[sampling_options] +mode = "basic" +layout = "vlba" +img_size = 256 +fov_center_ra = [00:00:00.0, 23:59:59.59] +fov_center_dec = [00:00:00.0, 11:59:59.59] +fov_size = 25.6 # 256 pixels, max res 0.1 +corr_int_time = 10.0 +scan_start = ["01-01-2020 00:00:01", "31-12-2021 23:59:59"] +scan_duration = [50, 300] +scans = [30, 72] +interval_length = 1200 +base_freq = 15.21e9 +frequsel = [0e8, 0.8e8, 1.44e8, 2.08e8] +bandwidths = [6.4e7, 6.4e7, 6.4e7, 6.4e7] + +[bundle_options] +num_bundles = 5 +size_bundles = 10 +in_path = "./build" +out_path = "./build" + +[cluster_options] +num_jobs = 4 \ No newline at end of file From 792ca1d3d2a68dcc5c9095a4dd410bd8694337c4 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Wed, 8 Dec 2021 12:14:45 +0100 Subject: [PATCH 15/63] functions for data set creation --- pyvisgen/simulation/data_set.py | 101 ++++++++++++++++++++++++++++++++ 1 file changed, 101 insertions(+) create mode 100644 pyvisgen/simulation/data_set.py diff --git a/pyvisgen/simulation/data_set.py b/pyvisgen/simulation/data_set.py new file mode 100644 index 0000000..e3176ae --- /dev/null +++ b/pyvisgen/simulation/data_set.py @@ -0,0 +1,101 @@ +import numpy as np +import pandas as pd +from pathlib import Path +from datetime import datetime +from pyvisgen.utils.config import read_data_set_conf +from pyvisgen.simulation.visibility import vis_loop +import pyvisgen.fits.writer as writer + + +def simulate_data_set(config, source_idx=0): + conf = read_data_set_conf(config) + out_path = Path(conf["out_path"]) + + out = out_path / Path("vis_" + str(source_idx) + ".fits") + print(out) + + samp_ops = create_sampling_rc(conf) + + hdu_list = writer.create_hdu_list(vis_loop(samp_ops, img_idx=0), samp_ops) + hdu_list.writeto(out, overwrite=True) + + +def create_sampling_rc(conf): + sampling_opts = draw_sampling_opts(conf) + samp_ops = { + "mode": sampling_opts[0], + "layout": sampling_opts[1], + "img_size": sampling_opts[2], + "fov_center_ra": sampling_opts[3], + "fov_center_dec": sampling_opts[4], + "fov_size": sampling_opts[5], + "corr_int_time": sampling_opts[6], + "scan_start": sampling_opts[7], + "scan_duration": sampling_opts[8], + "scans": sampling_opts[9], + "interval_length": sampling_opts[10], + "base_freq": sampling_opts[11], + "frequsel": sampling_opts[12], + "bandwidths": sampling_opts[13], + } + return samp_ops + + +def draw_sampling_opts(conf): + date_str_ra = pd.date_range( + conf["fov_center_ra"][0][0].strftime("%H:%M:%S"), + conf["fov_center_ra"][0][1].strftime("%H:%M:%S"), + freq="1min", + ).strftime("%H:%M:%S") + times_ra = [ + datetime.time(datetime.strptime(date, "%H:%M:%S")) for date in date_str_ra + ] + fov_center_ra = np.random.choice(times_ra) + date_str_dec = pd.date_range( + conf["fov_center_dec"][0][0].strftime("%H:%M:%S"), + conf["fov_center_dec"][0][1].strftime("%H:%M:%S"), + freq="1min", + ).strftime("%H:%M:%S") + times_dec = [ + datetime.time(datetime.strptime(date, "%H:%M:%S")) for date in date_str_dec + ] + fov_center_dec = np.random.choice(times_dec) + start_time_l = datetime.strptime(conf["scan_start"][0], "%d-%m-%Y %H:%M:%S") + start_time_h = datetime.strptime(conf["scan_start"][1], "%d-%m-%Y %H:%M:%S") + start_times = pd.date_range( + start_time_l, + start_time_h, + freq="1h", + ).strftime("%d-%m-%Y %H:%M:%S") + scan_start = np.random.choice( + [datetime.strptime(time, "%d-%m-%Y %H:%M:%S") for time in start_times] + ) + scan_duration = ( + np.random.randint(conf["scan_duration"][0] / 10, conf["scan_duration"][1] / 10) + * conf["corr_int_time"] + ) + scans = np.random.randint(conf["scans"][0], conf["scans"][1]) + opts = np.array( + [ + conf["mode"][0], + conf["layout"][0], + conf["img_size"][0], + fov_center_ra, + fov_center_dec, + conf["fov_size"], + conf["corr_int_time"], + scan_start, + scan_duration, + scans, + conf["interval_length"], + conf["base_freq"], + conf["frequsel"], + conf["bandwidths"], + ], + dtype="object", + ) + return opts + + +if __name__ == "__main__": + simulate_data_set("../../config/data_set.toml") From a8c2200a50e1189735803a1a1822a328d5ce8de5 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Wed, 8 Dec 2021 12:15:03 +0100 Subject: [PATCH 16/63] read data set config --- pyvisgen/utils/config.py | 74 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) create mode 100644 pyvisgen/utils/config.py diff --git a/pyvisgen/utils/config.py b/pyvisgen/utils/config.py new file mode 100644 index 0000000..d6cb2a0 --- /dev/null +++ b/pyvisgen/utils/config.py @@ -0,0 +1,74 @@ +import toml + + +def read_data_set_conf(conf_toml): + """Read toml data set configuration file and convert it into a dictionary. + + Parameters + ---------- + conf : toml file + path to config file + + Returns + ------- + sim_conf : dictionary + simulation configuration + """ + config = toml.load(conf_toml) + conf = {} + + conf["mode"] = config["sampling_options"]["mode"], + conf["layout"] = config["sampling_options"]["layout"], + conf["img_size"] = config["sampling_options"]["img_size"], + conf["fov_center_ra"] = config["sampling_options"]["fov_center_ra"], + conf["fov_center_dec"] = config["sampling_options"]["fov_center_dec"], + conf["fov_size"] = config["sampling_options"]["fov_size"] + conf["corr_int_time"] = config["sampling_options"]["corr_int_time"] + conf["scan_start"] = config["sampling_options"]["scan_start"] + conf["scan_duration"] = config["sampling_options"]["scan_duration"] + conf["scans"] = config["sampling_options"]["scans"] + conf["interval_length"] = config["sampling_options"]["interval_length"] + conf["base_freq"] = config["sampling_options"]["base_freq"] + conf["frequsel"] = config["sampling_options"]["frequsel"] + conf["bandwidths"] = config["sampling_options"]["bandwidths"] + + conf["num_bundles"] = config["bundle_options"]["num_bundles"] + conf["size_bundles"] = config["bundle_options"]["size_bundles"] + conf["num_images"] = conf["num_bundles"] * conf["size_bundles"] + conf["in_path"] = config["bundle_options"]["in_path"] + conf["out_path"] = config["bundle_options"]["out_path"] + + conf["num_jobs"] = config["cluster_options"]["num_jobs"] + + return conf + + +def read_config(conf): + """Read toml simulation configuration file and convert it into a dictionary. + + Parameters + ---------- + conf : toml file + path to config file + + Returns + ------- + sim_conf : dictionary + simulation configuration + """ + config = toml.load(conf) + sim_conf = {} + + sim_conf["src_coord"] = SkyCoord( + ra=config["sampling_options"]["fov_center_ra"], + dec=config["sampling_options"]["fov_center_dec"], + unit=(u.deg, u.deg), + ) + sim_conf["fov_size"] = config["sampling_options"]["fov_size"] + sim_conf["corr_int_time"] = config["sampling_options"]["corr_int_time"] + sim_conf["scan_start"] = config["sampling_options"]["scan_start"] + sim_conf["scan_duration"] = config["sampling_options"]["scan_duration"] + sim_conf["scans"] = config["sampling_options"]["scans"] + sim_conf["channel"] = config["sampling_options"]["channel"] + sim_conf["interval_length"] = config["sampling_options"]["interval_length"] + return sim_conf From 2317b601b900952fa2b10663c78451c6c08faf94 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Wed, 8 Dec 2021 12:15:30 +0100 Subject: [PATCH 17/63] loop to calc visibilities --- pyvisgen/simulation/visibility.py | 176 ++++++++++++++++++++++++++++++ 1 file changed, 176 insertions(+) create mode 100644 pyvisgen/simulation/visibility.py diff --git a/pyvisgen/simulation/visibility.py b/pyvisgen/simulation/visibility.py new file mode 100644 index 0000000..5f31bb5 --- /dev/null +++ b/pyvisgen/simulation/visibility.py @@ -0,0 +1,176 @@ +import numpy as np +from dataclasses import dataclass +import pyvisgen.simulation.utils as ut +import pyvisgen.layouts.layouts as layouts +import astropy.constants as const +from astropy import units as un +import pyvisgen.simulation.scan as scan +import torch +from astropy.coordinates import SkyCoord +from pyvisgen.utils.data import get_data_paths, h5_sky_distributions + + +@dataclass +class Visibilities: + SI: [complex] + SQ: [complex] + SU: [complex] + SV: [complex] + num: [int] + scan: [int] + base_num: [int] + u: [float] + v: [float] + w: [float] + date: [float] + _date: [float] + + def __getitem__(self, i): + baseline = Vis( + self.SI[i], + self.SQ[i], + self.SU[i], + self.SV[i], + self.num[i], + self.scan[i], + self.base_num[i], + self.u[i], + self.v[i], + self.w[i], + self.date[i], + self._date[i], + ) + return baseline + + def get_values(self): + return np.array([self.SI, self.SQ, self.SU, self.SV]) + + def add(self, visibilities): + self.SI = np.concatenate([self.SI, visibilities.SI]) + self.SQ = np.concatenate([self.SQ, visibilities.SQ]) + self.SU = np.concatenate([self.SU, visibilities.SU]) + self.SV = np.concatenate([self.SV, visibilities.SV]) + self.num = np.concatenate([self.num, visibilities.num]) + self.scan = np.concatenate([self.scan, visibilities.scan]) + self.base_num = np.concatenate([self.base_num, visibilities.base_num]) + self.u = np.concatenate([self.u, visibilities.u]) + self.v = np.concatenate([self.v, visibilities.v]) + self.w = np.concatenate([self.w, visibilities.w]) + self.date = np.concatenate([self.date, visibilities.date]) + self._date = np.concatenate([self._date, visibilities._date]) + + +@dataclass +class Vis: + SI: complex + SQ: complex + SU: complex + SV: complex + num: int + scan: int + base_num: int + u: float + v: float + w: float + date: float + _date: float + + +def vis_loop(rc, img_idx=0, num_threads=48): + # torch.set_num_threads(num_threads) + + # read config + array_layout = layouts.get_array_layout(rc["layout"]) + src_crd = SkyCoord( + ra=170, # rc["fov_center_ra"].strftime("%H:%M:%S"), + dec=22, # rc["fov_center_dec"].strftime("%H:%M:%S"), + unit=(un.deg, un.deg), + ) + rc["fov_center_ra"] = 170 + rc["fov_center_dec"] = 22 + + wave = np.array( + [ + const.c / ((rc["base_freq"] + float(freq)) / un.second) / un.meter + for freq in rc["frequsel"] + ] + ) + + # calculate rd, lm + rd = scan.rd_grid(rc["fov_size"] * np.pi / (3600 * 180), rc["img_size"], src_crd) + lm = scan.lm_grid(rd, src_crd) + + # calculate time steps + time = ut.calc_time_steps(rc) + + # open image + path = "../../../radiosim/build/example_data" + data_paths = get_data_paths(path) + h5 = h5_sky_distributions(data_paths) + img = torch.tensor(h5[img_idx]) + SI = torch.zeros((img.shape[1], img.shape[2], 4), dtype=torch.cdouble) + print(SI.shape) + SI[..., 0] = img + SI[..., 1] = img + + # number statiosn, number baselines + stat_num = array_layout.st_num.shape[0] + base_num = int(stat_num * (stat_num - 1) / 2) + + # calculate vis + visibilities = Visibilities([], [], [], [], [], [], [], [], [], [], [], []) + vis_num = np.zeros(1) + # i in total number of scans + for i in range(rc["scans"]): + end_idx = int((rc["scan_duration"] / rc["corr_int_time"]) + 1) + t = time[i * end_idx : (i + 1) * end_idx] + + baselines = scan.get_baselines(src_crd, t, array_layout) + + valid = baselines.valid.reshape(-1, base_num) + mask = np.array(valid[:-1]).astype(bool) & np.array(valid[1:]).astype(bool) + u = baselines.u.reshape(-1, base_num) + v = baselines.v.reshape(-1, base_num) + w = baselines.w.reshape(-1, base_num) + base_valid = np.arange(len(baselines.u)).reshape(-1, base_num)[:-1][mask] + u_valid = u[:-1][mask] + v_valid = v[:-1][mask] + w_valid = w[:-1][mask] + date = np.repeat( + (t[:-1] + rc["corr_int_time"] * un.second / 2).jd.reshape(-1, 1), + base_num, + axis=1, + )[mask] + print(len(base_valid)) + + _date = np.zeros(len(u_valid)) + + print(wave) + X1 = scan.uncorrupted(lm, baselines, wave[0], t, src_crd, array_layout, SI) + if X1.shape[0] == 1: + continue + X2 = scan.uncorrupted(lm, baselines, wave[0], t, src_crd, array_layout, SI) + + vis_num = np.arange(X1.shape[2] // 2) + 1 + vis_num.max() + + int_values = scan.integrate(X1, X2) + del X1, X2 + int_values = int_values.reshape(-1, 4) + + vis = Visibilities( + int_values[:, 0], + int_values[:, 1], + int_values[:, 2], + int_values[:, 3], + vis_num, + np.repeat(i + 1, len(vis_num)), + np.array([baselines[i].baselineNum() for i in base_valid]), + u_valid, + v_valid, + w_valid, + date, + _date, + ) + + visibilities.add(vis) + return visibilities From a65531dfcd168fca80f922ac58fed13c992bcad9 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Fri, 4 Mar 2022 03:17:29 +0100 Subject: [PATCH 18/63] add config to create data set --- config/data_set.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/config/data_set.toml b/config/data_set.toml index 4be20ab..c2b6980 100644 --- a/config/data_set.toml +++ b/config/data_set.toml @@ -1,10 +1,10 @@ [sampling_options] mode = "basic" layout = "vlba" -img_size = 256 +img_size = 64 fov_center_ra = [00:00:00.0, 23:59:59.59] fov_center_dec = [00:00:00.0, 11:59:59.59] -fov_size = 25.6 # 256 pixels, max res 0.1 +fov_size = 0.0064 # max res 0.1 corr_int_time = 10.0 scan_start = ["01-01-2020 00:00:01", "31-12-2021 23:59:59"] scan_duration = [50, 300] From 23b59d0c946c99dfc9106677fa5bcf59098fa885 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Fri, 4 Mar 2022 03:17:51 +0100 Subject: [PATCH 19/63] add script for data set creation --- pyvisgen/simulation/data_set.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/pyvisgen/simulation/data_set.py b/pyvisgen/simulation/data_set.py index e3176ae..fa84e43 100644 --- a/pyvisgen/simulation/data_set.py +++ b/pyvisgen/simulation/data_set.py @@ -1,22 +1,29 @@ import numpy as np import pandas as pd +import torch from pathlib import Path from datetime import datetime from pyvisgen.utils.config import read_data_set_conf from pyvisgen.simulation.visibility import vis_loop import pyvisgen.fits.writer as writer +from radiosim.data import radiosim_data def simulate_data_set(config, source_idx=0): + np.random.seed(42) conf = read_data_set_conf(config) out_path = Path(conf["out_path"]) out = out_path / Path("vis_" + str(source_idx) + ".fits") - print(out) samp_ops = create_sampling_rc(conf) - hdu_list = writer.create_hdu_list(vis_loop(samp_ops, img_idx=0), samp_ops) + # open image + path = "../../../test_radiosim/build/test_data" + data = radiosim_data(path) # conf["in_path"] + SI = torch.tensor(data[0][0][0], dtype=torch.cdouble) + + hdu_list = writer.create_hdu_list(vis_loop(samp_ops, SI), samp_ops) hdu_list.writeto(out, overwrite=True) From 89a0ca014a538ff7fb607ae6fe6295498a808f73 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Fri, 4 Mar 2022 03:18:21 +0100 Subject: [PATCH 20/63] only calculate I --- pyvisgen/simulation/scan.py | 16 +++++++++------- pyvisgen/simulation/visibility.py | 22 +++++----------------- 2 files changed, 14 insertions(+), 24 deletions(-) diff --git a/pyvisgen/simulation/scan.py b/pyvisgen/simulation/scan.py index 398af8f..f24db42 100644 --- a/pyvisgen/simulation/scan.py +++ b/pyvisgen/simulation/scan.py @@ -245,7 +245,7 @@ def lm_grid(rd_grid, src_crd): return lm_grid -def uncorrupted(lm, baselines, wave, time, src_crd, array_layout, I): +def uncorrupted(lm, baselines, wave, time, src_crd, array_layout, SI): """Calculates uncorrupted visibility Parameters @@ -283,14 +283,16 @@ def uncorrupted(lm, baselines, wave, time, src_crd, array_layout, I): K = getK(baselines, lm, wave, base_num) - B = np.zeros((lm.shape[0], lm.shape[1], 2, 2), dtype=complex) + B = np.zeros((lm.shape[0], lm.shape[1], 1), dtype=complex) - B[:, :, 0, 0] = I[:, :, 0] + I[:, :, 1] - B[:, :, 0, 1] = I[:, :, 2] + 1j * I[:, :, 3] - B[:, :, 1, 0] = I[:, :, 2] - 1j * I[:, :, 3] - B[:, :, 1, 1] = I[:, :, 0] - I[:, :, 1] + B[:, :, 0] = SI + SI + # # only calculate without polarization for the moment + # B[:, :, 0, 0] = SI[:, :, 0] + SI[:, :, 1] + # B[:, :, 0, 1] = SI[:, :, 2] + 1j * SI[:, :, 3] + # B[:, :, 1, 0] = SI[:, :, 2] - 1j * SI[:, :, 3] + # B[:, :, 1, 1] = SI[:, :, 0] - SI[:, :, 1] - X = torch.einsum("lmij,lmb->lmbij", torch.tensor(B), K) + X = torch.einsum("lmi,lmb->lmbi", torch.tensor(B), K) return X diff --git a/pyvisgen/simulation/visibility.py b/pyvisgen/simulation/visibility.py index 5f31bb5..1e500d0 100644 --- a/pyvisgen/simulation/visibility.py +++ b/pyvisgen/simulation/visibility.py @@ -76,7 +76,7 @@ class Vis: _date: float -def vis_loop(rc, img_idx=0, num_threads=48): +def vis_loop(rc, SI, num_threads=48): # torch.set_num_threads(num_threads) # read config @@ -103,16 +103,6 @@ def vis_loop(rc, img_idx=0, num_threads=48): # calculate time steps time = ut.calc_time_steps(rc) - # open image - path = "../../../radiosim/build/example_data" - data_paths = get_data_paths(path) - h5 = h5_sky_distributions(data_paths) - img = torch.tensor(h5[img_idx]) - SI = torch.zeros((img.shape[1], img.shape[2], 4), dtype=torch.cdouble) - print(SI.shape) - SI[..., 0] = img - SI[..., 1] = img - # number statiosn, number baselines stat_num = array_layout.st_num.shape[0] base_num = int(stat_num * (stat_num - 1) / 2) @@ -141,11 +131,9 @@ def vis_loop(rc, img_idx=0, num_threads=48): base_num, axis=1, )[mask] - print(len(base_valid)) _date = np.zeros(len(u_valid)) - print(wave) X1 = scan.uncorrupted(lm, baselines, wave[0], t, src_crd, array_layout, SI) if X1.shape[0] == 1: continue @@ -155,13 +143,13 @@ def vis_loop(rc, img_idx=0, num_threads=48): int_values = scan.integrate(X1, X2) del X1, X2 - int_values = int_values.reshape(-1, 4) + int_values = int_values vis = Visibilities( int_values[:, 0], - int_values[:, 1], - int_values[:, 2], - int_values[:, 3], + torch.zeros(int_values.shape[0], dtype=torch.complex128), + torch.zeros(int_values.shape[0], dtype=torch.complex128), + torch.zeros(int_values.shape[0], dtype=torch.complex128), vis_num, np.repeat(i + 1, len(vis_num)), np.array([baselines[i].baselineNum() for i in base_valid]), From 7817a2e083edc858e5a4905ca2a8c0a7861a3f3a Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Fri, 4 Mar 2022 03:18:49 +0100 Subject: [PATCH 21/63] black --- pyvisgen/fits/writer.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/pyvisgen/fits/writer.py b/pyvisgen/fits/writer.py index 0decf78..11b4b5f 100644 --- a/pyvisgen/fits/writer.py +++ b/pyvisgen/fits/writer.py @@ -18,9 +18,7 @@ def create_vis_hdu(data, conf, layout="vlba", source_name="sim-source-0"): data.date.min() ) # placeholder, julian date of vis, central time in the integration period - _DATE = ( - data._date - ) # central time in the integration period + _DATE = data._date # central time in the integration period BASELINE = data.base_num @@ -28,9 +26,12 @@ def create_vis_hdu(data, conf, layout="vlba", source_name="sim-source-0"): # visibility data values = data.get_values() + vis = np.swapaxes( np.swapaxes( - np.stack([values.real, values.imag, np.ones(values.shape)-0.8], axis=1), 1, 2 + np.stack([values.real, values.imag, np.ones(values.shape) - 0.8], axis=1), + 1, + 2, ), 0, 1, From 6902b668e02ad8fa0cb5c679c2b8f05e2ae37b3f Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Fri, 4 Mar 2022 04:08:38 +0100 Subject: [PATCH 22/63] use data routines from radiosim --- pyvisgen/utils/data.py | 59 ------------------------------------------ 1 file changed, 59 deletions(-) delete mode 100644 pyvisgen/utils/data.py diff --git a/pyvisgen/utils/data.py b/pyvisgen/utils/data.py deleted file mode 100644 index 431f664..0000000 --- a/pyvisgen/utils/data.py +++ /dev/null @@ -1,59 +0,0 @@ -import re -import h5py -import numpy as np -from pathlib import Path - - -def get_data_paths(path): - data_dir = Path(path) - bundles = np.array([x for x in data_dir.iterdir()]) - data_paths = np.sort([path for path in bundles if re.findall("source_", path.name)]) - return data_paths - - -class h5_sky_distributions: - def __init__(self, bundle_paths): - """ - Save the bundle paths and the number of bundles in one file. - """ - self.bundles = bundle_paths - self.num_img = len(self.open_bundle(self.bundles[0])) // 3 - - def __call__(self): - return print("This is the h5_sky_distributions class.") - - def __len__(self): - """ - Returns the total number of images in this dataset. - """ - return len(self.bundles) * self.num_img - - def __getitem__(self, i): - sky = self.open_image("sky", i) - return sky - - def open_bundle(self, bundle_path): - bundle = h5py.File(bundle_path, "r") - return bundle - - def open_image(self, var, i): - if isinstance(i, int): - i = np.array([i]) - indices = np.sort(i) - bundle = indices // self.num_img - image = indices - bundle * self.num_img - bundle_unique = np.unique(bundle) - bundle_files = [ - h5py.File(self.bundles[bundle], "r") for bundle in bundle_unique - ] - bundle_files_str = list(map(str, bundle_files)) - data = np.array( - [ - bund["sky" + str(img)] - for bund, bund_str in zip(bundle_files, bundle_files_str) - for img in image[ - bundle == bundle_unique[bundle_files_str.index(bund_str)] - ] - ] - ) - return data From d7a27640ad238661edb0a0aba8ee3a9c70698701 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Fri, 4 Mar 2022 04:11:33 +0100 Subject: [PATCH 23/63] use in_path from config file --- config/data_set.toml | 2 +- pyvisgen/simulation/data_set.py | 3 +-- pyvisgen/simulation/visibility.py | 1 - 3 files changed, 2 insertions(+), 4 deletions(-) diff --git a/config/data_set.toml b/config/data_set.toml index c2b6980..d1f77c5 100644 --- a/config/data_set.toml +++ b/config/data_set.toml @@ -17,7 +17,7 @@ bandwidths = [6.4e7, 6.4e7, 6.4e7, 6.4e7] [bundle_options] num_bundles = 5 size_bundles = 10 -in_path = "./build" +in_path = "../../../test_radiosim/build/test_data" out_path = "./build" [cluster_options] diff --git a/pyvisgen/simulation/data_set.py b/pyvisgen/simulation/data_set.py index fa84e43..c8ecb76 100644 --- a/pyvisgen/simulation/data_set.py +++ b/pyvisgen/simulation/data_set.py @@ -19,8 +19,7 @@ def simulate_data_set(config, source_idx=0): samp_ops = create_sampling_rc(conf) # open image - path = "../../../test_radiosim/build/test_data" - data = radiosim_data(path) # conf["in_path"] + data = radiosim_data(conf["in_path"]) SI = torch.tensor(data[0][0][0], dtype=torch.cdouble) hdu_list = writer.create_hdu_list(vis_loop(samp_ops, SI), samp_ops) diff --git a/pyvisgen/simulation/visibility.py b/pyvisgen/simulation/visibility.py index 1e500d0..3cd1790 100644 --- a/pyvisgen/simulation/visibility.py +++ b/pyvisgen/simulation/visibility.py @@ -7,7 +7,6 @@ import pyvisgen.simulation.scan as scan import torch from astropy.coordinates import SkyCoord -from pyvisgen.utils.data import get_data_paths, h5_sky_distributions @dataclass From 120a4a4d351c328c321b962f26511825f7c0f87e Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Fri, 4 Mar 2022 04:18:58 +0100 Subject: [PATCH 24/63] add missing import --- pyvisgen/utils/config.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/pyvisgen/utils/config.py b/pyvisgen/utils/config.py index d6cb2a0..ad5ac3d 100644 --- a/pyvisgen/utils/config.py +++ b/pyvisgen/utils/config.py @@ -1,4 +1,6 @@ import toml +import astropy.units as un +from astropy.coordinates import SkyCoord def read_data_set_conf(conf_toml): @@ -17,11 +19,11 @@ def read_data_set_conf(conf_toml): config = toml.load(conf_toml) conf = {} - conf["mode"] = config["sampling_options"]["mode"], - conf["layout"] = config["sampling_options"]["layout"], - conf["img_size"] = config["sampling_options"]["img_size"], - conf["fov_center_ra"] = config["sampling_options"]["fov_center_ra"], - conf["fov_center_dec"] = config["sampling_options"]["fov_center_dec"], + conf["mode"] = (config["sampling_options"]["mode"],) + conf["layout"] = (config["sampling_options"]["layout"],) + conf["img_size"] = (config["sampling_options"]["img_size"],) + conf["fov_center_ra"] = (config["sampling_options"]["fov_center_ra"],) + conf["fov_center_dec"] = (config["sampling_options"]["fov_center_dec"],) conf["fov_size"] = config["sampling_options"]["fov_size"] conf["corr_int_time"] = config["sampling_options"]["corr_int_time"] conf["scan_start"] = config["sampling_options"]["scan_start"] @@ -62,7 +64,7 @@ def read_config(conf): sim_conf["src_coord"] = SkyCoord( ra=config["sampling_options"]["fov_center_ra"], dec=config["sampling_options"]["fov_center_dec"], - unit=(u.deg, u.deg), + unit=(un.deg, un.deg), ) sim_conf["fov_size"] = config["sampling_options"]["fov_size"] sim_conf["corr_int_time"] = config["sampling_options"]["corr_int_time"] From fdbb023c75fcbca0654cf087c2eb820ff6c2024d Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Fri, 4 Mar 2022 04:19:20 +0100 Subject: [PATCH 25/63] loop over whole data set --- pyvisgen/simulation/data_set.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pyvisgen/simulation/data_set.py b/pyvisgen/simulation/data_set.py index c8ecb76..7ca797b 100644 --- a/pyvisgen/simulation/data_set.py +++ b/pyvisgen/simulation/data_set.py @@ -20,10 +20,11 @@ def simulate_data_set(config, source_idx=0): # open image data = radiosim_data(conf["in_path"]) - SI = torch.tensor(data[0][0][0], dtype=torch.cdouble) + for i in range(len(data)): + SI = torch.tensor(data[i][0][0], dtype=torch.cdouble) - hdu_list = writer.create_hdu_list(vis_loop(samp_ops, SI), samp_ops) - hdu_list.writeto(out, overwrite=True) + hdu_list = writer.create_hdu_list(vis_loop(samp_ops, SI), samp_ops) + hdu_list.writeto(out, overwrite=True) def create_sampling_rc(conf): From 6982bf899652316c34169bc20b460ec862a5e813 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Fri, 4 Mar 2022 06:10:15 +0100 Subject: [PATCH 26/63] drop unused import --- pyvisgen/simulation/scan.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyvisgen/simulation/scan.py b/pyvisgen/simulation/scan.py index f24db42..64c24bc 100644 --- a/pyvisgen/simulation/scan.py +++ b/pyvisgen/simulation/scan.py @@ -12,7 +12,6 @@ import itertools import time as t import numexpr as ne # fast exponential -from einsumt import einsumt as einsum @dataclass From 6ffb5a67cb128aaa0045325b53ce976f122edb93 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Fri, 4 Mar 2022 06:10:49 +0100 Subject: [PATCH 27/63] add console script --- pyvisgen/simulation/scripts/create_dataset.py | 23 +++++++++++++++++++ setup.py | 5 ++++ 2 files changed, 28 insertions(+) create mode 100644 pyvisgen/simulation/scripts/create_dataset.py diff --git a/pyvisgen/simulation/scripts/create_dataset.py b/pyvisgen/simulation/scripts/create_dataset.py new file mode 100644 index 0000000..1030385 --- /dev/null +++ b/pyvisgen/simulation/scripts/create_dataset.py @@ -0,0 +1,23 @@ +import click +from pyvisgen.simulation.data_set import simulate_data_set + + +@click.command() +@click.argument("configuration_path", type=click.Path(exists=True, dir_okay=False)) +# @click.option( +# "--mode", +# type=click.Choice( +# [ +# "simulate", +# "overview", +# ], +# case_sensitive=False, +# ), +# default="simulate", +# ) +def main(configuration_path): + simulate_data_set(configuration_path) + + +if __name__ == "__main__": + main() diff --git a/setup.py b/setup.py index a094a11..19791a9 100644 --- a/setup.py +++ b/setup.py @@ -27,6 +27,11 @@ setup_requires=["pytest-runner"], tests_require=["pytest"], zip_safe=False, + entry_points={ + "console_scripts": [ + "pyvisgen_create_dataset = pyvisgen.simulation.scripts.create_dataset:main", + ], + }, classifiers=[ "Development Status :: 3 - Alpha", "Intended Audience :: Science/Research", From c66b405b579f8286f6461dd2da422cd58e4e0806 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Fri, 4 Mar 2022 06:31:50 +0100 Subject: [PATCH 28/63] fix layout loading in writer --- pyvisgen/fits/writer.py | 12 ++++++------ pyvisgen/layouts/layouts.py | 7 +++++-- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/pyvisgen/fits/writer.py b/pyvisgen/fits/writer.py index 11b4b5f..0f97f4b 100644 --- a/pyvisgen/fits/writer.py +++ b/pyvisgen/fits/writer.py @@ -3,8 +3,8 @@ from astropy import wcs import astropy.units as un from astropy.time import Time -import pandas as pd import astropy.constants as const +import pyvisgen.layouts.layouts as layouts def create_vis_hdu(data, conf, layout="vlba", source_name="sim-source-0"): @@ -233,8 +233,8 @@ def create_frequency_hdu(conf): return hdu_freq -def create_antenna_hdu(layout_txt, conf, layout="vlba"): - array = pd.read_csv(layout_txt, sep=" ") +def create_antenna_hdu(conf): + array = layouts.get_array_layout(conf["layout"], writer=True) ANNAME = np.chararray(len(array), itemsize=8, unicode=True) ANNAME[:] = array["station_name"].values @@ -334,7 +334,7 @@ def create_antenna_hdu(layout_txt, conf, layout="vlba"): hdu_ant.header["UT1UTC"] = (0, "UT1 - UTC (sec)") # missing hdu_ant.header["DATUTC"] = (0, "time system - UTC (sec)") # missing hdu_ant.header["TIMSYS"] = ("UTC", "Time system") - hdu_ant.header["ARRNAM"] = (layout, "Array name") + hdu_ant.header["ARRNAM"] = (conf["layout"], "Array name") hdu_ant.header["XYZHAND"] = ("RIGHT", "Handedness of station coordinates") hdu_ant.header["FRAME"] = ("????", "Coordinate frame") hdu_ant.header["NUMORB"] = (0, "Number orbital parameters in table (n orb)") @@ -369,10 +369,10 @@ def create_antenna_hdu(layout_txt, conf, layout="vlba"): return hdu_ant -def create_hdu_list(data, conf, path="../layouts/vlba.txt"): +def create_hdu_list(data, conf): vis_hdu = create_vis_hdu(data, conf) time_hdu = create_time_hdu(data) freq_hdu = create_frequency_hdu(conf) - ant_hdu = create_antenna_hdu(path, conf) + ant_hdu = create_antenna_hdu(conf) hdu_list = fits.HDUList([vis_hdu, time_hdu, freq_hdu, ant_hdu]) return hdu_list diff --git a/pyvisgen/layouts/layouts.py b/pyvisgen/layouts/layouts.py index f6f49e7..915a723 100644 --- a/pyvisgen/layouts/layouts.py +++ b/pyvisgen/layouts/layouts.py @@ -56,7 +56,7 @@ class Station: altitude: float -def get_array_layout(array_name): +def get_array_layout(array_name, writer=False): """Reads telescope layout txt file and converts it into a dataclass. Available arrays: - EHT @@ -85,4 +85,7 @@ def get_array_layout(array_name): array["SEFD"].values, array["altitude"].values, ) - return stations + if writer: + return array + else: + return stations From cb767833c7b0c3ae0dab757490c4705d5401fa54 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Fri, 4 Mar 2022 06:32:07 +0100 Subject: [PATCH 29/63] fix out_dir creation --- pyvisgen/simulation/data_set.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pyvisgen/simulation/data_set.py b/pyvisgen/simulation/data_set.py index 7ca797b..2b87e22 100644 --- a/pyvisgen/simulation/data_set.py +++ b/pyvisgen/simulation/data_set.py @@ -13,6 +13,7 @@ def simulate_data_set(config, source_idx=0): np.random.seed(42) conf = read_data_set_conf(config) out_path = Path(conf["out_path"]) + out_path.mkdir(exist_ok=True) out = out_path / Path("vis_" + str(source_idx) + ".fits") From a933c97f9edff4299b0458d8c16dd4adc5a48586 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Fri, 4 Mar 2022 06:38:49 +0100 Subject: [PATCH 30/63] fix out_paths --- pyvisgen/simulation/data_set.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pyvisgen/simulation/data_set.py b/pyvisgen/simulation/data_set.py index 2b87e22..0b9a531 100644 --- a/pyvisgen/simulation/data_set.py +++ b/pyvisgen/simulation/data_set.py @@ -15,13 +15,12 @@ def simulate_data_set(config, source_idx=0): out_path = Path(conf["out_path"]) out_path.mkdir(exist_ok=True) - out = out_path / Path("vis_" + str(source_idx) + ".fits") - samp_ops = create_sampling_rc(conf) # open image data = radiosim_data(conf["in_path"]) for i in range(len(data)): + out = out_path / Path("vis_" + str(i) + ".fits") SI = torch.tensor(data[i][0][0], dtype=torch.cdouble) hdu_list = writer.create_hdu_list(vis_loop(samp_ops, SI), samp_ops) From 75cbb8e27b62689efca61232dea74dfcca241c97 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Sat, 5 Mar 2022 06:07:49 +0100 Subject: [PATCH 31/63] add basic gridding routine --- pyvisgen/gridding/gridder.py | 123 +++++++++++++++++++++++++++++++++++ 1 file changed, 123 insertions(+) create mode 100644 pyvisgen/gridding/gridder.py diff --git a/pyvisgen/gridding/gridder.py b/pyvisgen/gridding/gridder.py new file mode 100644 index 0000000..8bc291d --- /dev/null +++ b/pyvisgen/gridding/gridder.py @@ -0,0 +1,123 @@ +import numpy as np +from tqdm import tqdm +import matplotlib.pyplot as plt +from pyvisgen.fits.data import fits_data + + +def create_gridded_dataset(fits_data_path): + fits_files = fits_data(fits_data_path) + size = len(fits_files) + + for i in tqdm(range(1)): + freq_data, base_freq = fits_files.get_freq_data(i) + uv_data = fits_files.get_uv_data(i) + + gridded_data = grid_data(uv_data, base_freq) + + +def grid_data(uv_data, base_freq): + freq = base_freq + + cmplx = uv_data["DATA"] + real = np.squeeze(cmplx[..., 0, 0]) + imag = np.squeeze(cmplx[..., 0, 1]) + # weight = np.squeeze(cmplx[..., 0, 2]) + ap = np.sqrt(real ** 2 + imag ** 2) + ph = np.angle(real + 1j * imag) + samps = np.array( + [ + np.append(uv_data["UU--"] * freq, -uv_data["UU--"] * freq), + np.append(uv_data["VV--"] * freq, -uv_data["VV--"] * freq), + np.append(ap, ap), + np.append(ph, -ph), + ] + ) + + # Generate Mask + u_0 = samps[0] + v_0 = samps[1] + N = 256 # image size, needs to be from config file + mask = np.zeros((N, N, u_0.shape[0])) + fov = ( + 0.0256 * np.pi / (3600 * 180) + ) # hard code #default 0.00018382, FoV from VLBA 163.7 + print(fov) + delta_u = 1 / (fov) + print(delta_u) + for i in range(N): + for j in range(N): + u_cell = (j - N / 2) * delta_u + v_cell = (i - N / 2) * delta_u + mask[i, j] = ((u_cell <= u_0) & (u_0 <= u_cell + delta_u)) & ( + (v_cell <= v_0) & (v_0 <= v_cell + delta_u) + ) + + mask = np.flip(mask, [0]) + points = np.sum(mask, 2) + points[points == 0] = 1 + gridded_vis = np.zeros((2, N, N)) + gridded_vis[0] = np.matmul(mask, samps[2].T) / points + gridded_vis[0] = (np.log10(gridded_vis[0] + 1e-10) / 10) + 1 + gridded_vis[1] = np.matmul(mask, samps[3].T) / points + + print(gridded_vis.shape) + + plt.imshow(gridded_vis[1]) + plt.colorbar() + plt.show() + + return 1 + + +if __name__ == "__main__": + create_gridded_dataset("/net/big-tank/POOL/projects/radio/test_rime/build/uvfits") + + # img = np.zeros((size, 256, 256)) + # samps = np.zeros((size, 4, 21000)) # hard code + # for i in tqdm(range(size)): + # sampled = bundle_paths_input[i] + # target = bundle_paths_target[i] + + # img[i] = np.asarray(Image.open(str(target))) + # # img[i] = img[i]/np.sum(img[i]) + + # print(f"\n Gridding VLBI data set.\n") + + # # Generate Mask + # u_0 = samps[0][0] + # v_0 = samps[0][1] + # N = 63 # hard code + # mask = np.zeros((N, N, u_0.shape[0])) + # fov = 0.00018382 * np.pi / (3600 * 180) # hard code #default 0.00018382 + # # delta_u = 1/(fov*N/256) # hard code + # delta_u = 1 / (fov) + # for i in range(N): + # for j in range(N): + # u_cell = (j - N / 2) * delta_u + # v_cell = (i - N / 2) * delta_u + # mask[i, j] = ((u_cell <= u_0) & (u_0 <= u_cell + delta_u)) & ( + # (v_cell <= v_0) & (v_0 <= v_cell + delta_u) + # ) + + # mask = np.flip(mask, [0]) + # points = np.sum(mask, 2) + # points[points == 0] = 1 + # samp_img = np.zeros((size, 2, N, N)) + # img_resized = np.zeros((size, N, N)) + # for i in tqdm(range(samps.shape[0])): + # samp_img[i][0] = np.matmul(mask, samps[i][2].T) / points + # samp_img[i][0] = (np.log10(samp_img[i][0] + 1e-10) / 10) + 1 + # samp_img[i][1] = np.matmul(mask, samps[i][3].T) / points + # img_resized[i] = cv2.resize(img[i], (N, N)) + # img_resized[i] = img_resized[i] / np.sum(img_resized[i]) + + # # truth_fft = np.array([np.fft.fft2(np.fft.fftshift(img)) for im in img_resized]) + # truth_fft = np.fft.fftshift( + # np.fft.fft2(np.fft.fftshift(img_resized, axes=(1, 2)), axes=(1, 2)), axes=(1, 2) + # ) + # fft_scaled_truth = prepare_fft_images(truth_fft, True, False) + + # out = data_path + "/samp_train0.h5" + # save_fft_pair(out, samp_img[:2300], fft_scaled_truth[:2300]) + # out = data_path + "/samp_valid0.h5" + # save_fft_pair(out, samp_img[2300:], fft_scaled_truth[2300:]) \ No newline at end of file From 8325c272efdad1655a49f1423d258c1af06e43eb Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Sat, 5 Mar 2022 06:08:05 +0100 Subject: [PATCH 32/63] add fits data set class --- pyvisgen/fits/data.py | 53 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) create mode 100644 pyvisgen/fits/data.py diff --git a/pyvisgen/fits/data.py b/pyvisgen/fits/data.py new file mode 100644 index 0000000..9c78b83 --- /dev/null +++ b/pyvisgen/fits/data.py @@ -0,0 +1,53 @@ +import numpy as np +from astropy.io import fits +from pathlib import Path + + +class fits_data: + def __init__(self, data_path): + """Handle h5 files simulated with radiosim package. + + Parameters + ---------- + data_path: str + path to fits data directory + + """ + self.path = data_path + self.files = self.get_files(Path(data_path)) + self.num_files = len(self.files) + + def __call__(self): + return print("This is the pyvisgen fits files data set class.") + + def __len__(self): + """ + Returns the total number of fits files in this dataset. + """ + return self.num_files + + def __getitem__(self, i): + return self.open_file(i) + + def get_files(self, path): + return np.array([x for x in path.iterdir()]) + + def get_uv_data(self, i): + with fits.open(self.files[i]) as hdul: + uv_data = hdul[0].data + return uv_data + + def get_freq_data(self, i): + with fits.open(self.files[i]) as hdul: + base_freq = hdul[0].header["CRVAL4"] + freq_data = hdul[2].data + return freq_data, base_freq + + def open_file(self, i): + if isinstance(i, np.ndarray): + print("Only one file can be open at the same time!") + return + + with fits.open(self.files[i]) as hdul: + fits_file = hdul + return fits_file.info() From 3c9e8c6771629c289008e30692f3b01f55bbc390 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Sat, 5 Mar 2022 06:09:01 +0100 Subject: [PATCH 33/63] add tqdm --- pyvisgen/simulation/data_set.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyvisgen/simulation/data_set.py b/pyvisgen/simulation/data_set.py index 0b9a531..dceb9a7 100644 --- a/pyvisgen/simulation/data_set.py +++ b/pyvisgen/simulation/data_set.py @@ -1,6 +1,7 @@ import numpy as np import pandas as pd import torch +from tqdm import tqdm from pathlib import Path from datetime import datetime from pyvisgen.utils.config import read_data_set_conf @@ -19,7 +20,7 @@ def simulate_data_set(config, source_idx=0): # open image data = radiosim_data(conf["in_path"]) - for i in range(len(data)): + for i in tqdm(range(len(data))): out = out_path / Path("vis_" + str(i) + ".fits") SI = torch.tensor(data[i][0][0], dtype=torch.cdouble) From 3f1895070525a1ddbd3149325b241efed22c9fdd Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Sun, 6 Mar 2022 06:58:06 +0100 Subject: [PATCH 34/63] set up bundle options --- pyvisgen/gridding/gridder.py | 94 +++++++++++++++++++++++++----------- 1 file changed, 66 insertions(+), 28 deletions(-) diff --git a/pyvisgen/gridding/gridder.py b/pyvisgen/gridding/gridder.py index 8bc291d..b61d0d2 100644 --- a/pyvisgen/gridding/gridder.py +++ b/pyvisgen/gridding/gridder.py @@ -1,49 +1,85 @@ import numpy as np from tqdm import tqdm import matplotlib.pyplot as plt +from pathlib import Path from pyvisgen.fits.data import fits_data +from pyvisgen.utils.config import read_data_set_conf -def create_gridded_dataset(fits_data_path): - fits_files = fits_data(fits_data_path) - size = len(fits_files) - - for i in tqdm(range(1)): - freq_data, base_freq = fits_files.get_freq_data(i) - uv_data = fits_files.get_uv_data(i) - - gridded_data = grid_data(uv_data, base_freq) - +def create_gridded_data_set(config): + conf = read_data_set_conf(config) + out_path_fits = Path(conf["out_path"]) + out_path = out_path_fits.parent / "gridded_data/" + ###### + out_path = Path( + "/net/big-tank/POOL/projects/radio/test_rime/build_test/gridded_data" + ) + out_path_fits = Path( + "/net/big-tank/POOL/projects/radio/test_rime/build_test/uvfits" + ) + ###### + print(out_path) + out_path.mkdir(parents=True, exist_ok=True) -def grid_data(uv_data, base_freq): - freq = base_freq + fits_files = fits_data(out_path_fits) + size = len(fits_files) + num_bundles = int(size // conf["bundle_size"]) + if conf["num_test_images"] > 0: + budles_test = int(conf["num_test_images"] // conf["bundle_size"]) + size -= conf["num_test_images"] + else: + bundle_test = 0 + + size_valid = conf["train_valid_split"] * size + size_train = size - size_valid + bundle_train = int(size_train // conf["bundle_size"]) + bundle_valid = int(size_valid // conf["bundle_size"]) + + for i in tqdm(range(bundle_train)): + uv_data_train = np.array( + [ + fits_files.get_uv_data(n) + for n in np.arange( + i * conf["bundle_size"], + (i * conf["bundle_size"]) + conf["bundle_size"], + ) + ], + dtype="object", + ) + gridded_data_train = np.array( + [grid_data(data, conf) for data in tqdm(uv_data_train)] + ) + print(gridded_data_train.shape) + print(gridded_data_train[0].shape) + + +def grid_data(uv_data, conf): + freq = conf["base_freq"] cmplx = uv_data["DATA"] real = np.squeeze(cmplx[..., 0, 0]) imag = np.squeeze(cmplx[..., 0, 1]) # weight = np.squeeze(cmplx[..., 0, 2]) - ap = np.sqrt(real ** 2 + imag ** 2) - ph = np.angle(real + 1j * imag) + # ap = np.sqrt(real ** 2 + imag ** 2) + # ph = np.angle(real + 1j * imag) samps = np.array( [ np.append(uv_data["UU--"] * freq, -uv_data["UU--"] * freq), np.append(uv_data["VV--"] * freq, -uv_data["VV--"] * freq), - np.append(ap, ap), - np.append(ph, -ph), + np.append(real, real), + np.append(imag, -imag), ] ) # Generate Mask u_0 = samps[0] v_0 = samps[1] - N = 256 # image size, needs to be from config file + N = conf["grid_size"] # image size, needs to be from config file mask = np.zeros((N, N, u_0.shape[0])) fov = ( - 0.0256 * np.pi / (3600 * 180) - ) # hard code #default 0.00018382, FoV from VLBA 163.7 - print(fov) + conf["fov_size"] * np.pi / (3600 * 180) + ) # hard code #default 0.00018382, FoV from VLBA 163.7 <- wrong! depends on setting of simulations delta_u = 1 / (fov) - print(delta_u) for i in range(N): for j in range(N): u_cell = (j - N / 2) * delta_u @@ -57,20 +93,22 @@ def grid_data(uv_data, base_freq): points[points == 0] = 1 gridded_vis = np.zeros((2, N, N)) gridded_vis[0] = np.matmul(mask, samps[2].T) / points - gridded_vis[0] = (np.log10(gridded_vis[0] + 1e-10) / 10) + 1 + # gridded_vis[0] = (np.log10(gridded_vis[0] + 1e-10) / 10) + 1 gridded_vis[1] = np.matmul(mask, samps[3].T) / points - print(gridded_vis.shape) + # print(gridded_vis.shape) - plt.imshow(gridded_vis[1]) - plt.colorbar() - plt.show() + # plt.imshow(gridded_vis[1]) + # plt.colorbar() + # plt.show() - return 1 + return gridded_vis if __name__ == "__main__": - create_gridded_dataset("/net/big-tank/POOL/projects/radio/test_rime/build/uvfits") + create_gridded_data_set( + "/net/big-tank/POOL/projects/radio/test_rime/create_dataset.toml" + ) # img = np.zeros((size, 256, 256)) # samps = np.zeros((size, 4, 21000)) # hard code From fd457de82646b0e21d471f3abe868bbc019dc763 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Sun, 6 Mar 2022 06:59:01 +0100 Subject: [PATCH 35/63] new sampling rc for each source --- pyvisgen/simulation/data_set.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/pyvisgen/simulation/data_set.py b/pyvisgen/simulation/data_set.py index dceb9a7..b829c0a 100644 --- a/pyvisgen/simulation/data_set.py +++ b/pyvisgen/simulation/data_set.py @@ -10,13 +10,11 @@ from radiosim.data import radiosim_data -def simulate_data_set(config, source_idx=0): +def simulate_data_set(config): np.random.seed(42) conf = read_data_set_conf(config) out_path = Path(conf["out_path"]) - out_path.mkdir(exist_ok=True) - - samp_ops = create_sampling_rc(conf) + out_path.mkdir(parents=True, exist_ok=True) # open image data = radiosim_data(conf["in_path"]) @@ -24,6 +22,7 @@ def simulate_data_set(config, source_idx=0): out = out_path / Path("vis_" + str(i) + ".fits") SI = torch.tensor(data[i][0][0], dtype=torch.cdouble) + samp_ops = create_sampling_rc(conf) hdu_list = writer.create_hdu_list(vis_loop(samp_ops, SI), samp_ops) hdu_list.writeto(out, overwrite=True) From 56b80b87c2686bfb4feb79886d7b7ea0f2f8e169 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Sun, 6 Mar 2022 06:59:31 +0100 Subject: [PATCH 36/63] add gridding options --- pyvisgen/utils/config.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pyvisgen/utils/config.py b/pyvisgen/utils/config.py index ad5ac3d..3358a68 100644 --- a/pyvisgen/utils/config.py +++ b/pyvisgen/utils/config.py @@ -34,9 +34,10 @@ def read_data_set_conf(conf_toml): conf["frequsel"] = config["sampling_options"]["frequsel"] conf["bandwidths"] = config["sampling_options"]["bandwidths"] - conf["num_bundles"] = config["bundle_options"]["num_bundles"] - conf["size_bundles"] = config["bundle_options"]["size_bundles"] - conf["num_images"] = conf["num_bundles"] * conf["size_bundles"] + conf["num_test_images"] = config["bundle_options"]["num_test_images"] + conf["bundle_size"] = config["bundle_options"]["bundle_size"] + conf["train_valid_split"] = config["bundle_options"]["train_valid_split"] + conf["grid_size"] = config["bundle_options"]["grid_size"] conf["in_path"] = config["bundle_options"]["in_path"] conf["out_path"] = config["bundle_options"]["out_path"] From 705809fbd89978239a17a542a4775f80f39fafb2 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Sun, 6 Mar 2022 06:59:54 +0100 Subject: [PATCH 37/63] add gridding mode --- pyvisgen/simulation/scripts/create_dataset.py | 31 +++++++++++-------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/pyvisgen/simulation/scripts/create_dataset.py b/pyvisgen/simulation/scripts/create_dataset.py index 1030385..e664a99 100644 --- a/pyvisgen/simulation/scripts/create_dataset.py +++ b/pyvisgen/simulation/scripts/create_dataset.py @@ -1,22 +1,27 @@ import click from pyvisgen.simulation.data_set import simulate_data_set +from pyvisgen.gridding.gridder import create_gridded_data_set @click.command() @click.argument("configuration_path", type=click.Path(exists=True, dir_okay=False)) -# @click.option( -# "--mode", -# type=click.Choice( -# [ -# "simulate", -# "overview", -# ], -# case_sensitive=False, -# ), -# default="simulate", -# ) -def main(configuration_path): - simulate_data_set(configuration_path) +@click.option( + "--mode", + type=click.Choice( + [ + "simulate", + "gridding", + ], + case_sensitive=False, + ), + default="simulate", +) +def main(configuration_path, mode): + if mode == "simulate": + print("hi") + simulate_data_set(configuration_path) + if mode == "gridding": + create_gridded_data_set(configuration_path) if __name__ == "__main__": From 94cc7c2e952ed4d8d24959dc64e2c01b3b85f75f Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Mon, 7 Mar 2022 06:27:58 +0100 Subject: [PATCH 38/63] sort filenames --- pyvisgen/fits/data.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyvisgen/fits/data.py b/pyvisgen/fits/data.py index 9c78b83..dd8c8de 100644 --- a/pyvisgen/fits/data.py +++ b/pyvisgen/fits/data.py @@ -30,7 +30,7 @@ def __getitem__(self, i): return self.open_file(i) def get_files(self, path): - return np.array([x for x in path.iterdir()]) + return np.sort(np.array([x for x in path.iterdir()])) def get_uv_data(self, i): with fits.open(self.files[i]) as hdul: From 89da75e6c72a5426821396e49896429d23a0a4d9 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Mon, 7 Mar 2022 06:28:45 +0100 Subject: [PATCH 39/63] gridder v1 --- pyvisgen/gridding/gridder.py | 196 +++++++++++++++++++++-------------- 1 file changed, 120 insertions(+), 76 deletions(-) diff --git a/pyvisgen/gridding/gridder.py b/pyvisgen/gridding/gridder.py index b61d0d2..d25b2d4 100644 --- a/pyvisgen/gridding/gridder.py +++ b/pyvisgen/gridding/gridder.py @@ -1,9 +1,11 @@ +import cv2 import numpy as np from tqdm import tqdm -import matplotlib.pyplot as plt from pathlib import Path +from radiosim.data import radiosim_data from pyvisgen.fits.data import fits_data from pyvisgen.utils.config import read_data_set_conf +from radionets.dl_framework.data import save_fft_pair def create_gridded_data_set(config): @@ -11,24 +13,27 @@ def create_gridded_data_set(config): out_path_fits = Path(conf["out_path"]) out_path = out_path_fits.parent / "gridded_data/" ###### - out_path = Path( - "/net/big-tank/POOL/projects/radio/test_rime/build_test/gridded_data" - ) - out_path_fits = Path( - "/net/big-tank/POOL/projects/radio/test_rime/build_test/uvfits" - ) + # out_path = Path( + # "/net/big-tank/POOL/projects/radio/test_rime/build_test/gridded_data" + # ) + # out_path_fits = Path( + # "/net/big-tank/POOL/projects/radio/test_rime/build_test/uvfits" + # ) ###### - print(out_path) out_path.mkdir(parents=True, exist_ok=True) + # conf[ + # "in_path" + # ] = "/net/big-tank/POOL/projects/radio/test_rime/build_test/sky_simulations" + sky_dist = radiosim_data(conf["in_path"]) fits_files = fits_data(out_path_fits) size = len(fits_files) - num_bundles = int(size // conf["bundle_size"]) - if conf["num_test_images"] > 0: - budles_test = int(conf["num_test_images"] // conf["bundle_size"]) - size -= conf["num_test_images"] - else: - bundle_test = 0 + # num_bundles = int(size // conf["bundle_size"]) + # if conf["num_test_images"] > 0: + # bundle_test = int(conf["num_test_images"] // conf["bundle_size"]) + # size -= conf["num_test_images"] + # else: + # bundle_test = 0 size_valid = conf["train_valid_split"] * size size_train = size - size_valid @@ -49,8 +54,91 @@ def create_gridded_data_set(config): gridded_data_train = np.array( [grid_data(data, conf) for data in tqdm(uv_data_train)] ) - print(gridded_data_train.shape) - print(gridded_data_train[0].shape) + sky_dist_train = np.array( + [ + cv2.resize( + sky_dist[int(n)][0][0], (conf["grid_size"], conf["grid_size"]) + ) + for n in np.arange( + i * conf["bundle_size"], + (i * conf["bundle_size"]) + conf["bundle_size"], + ) + ] + ) + truth_fft_train = np.fft.fftshift( + np.fft.fft2(np.fft.fftshift(sky_dist_train, axes=(1, 2)), axes=(1, 2)), + axes=(1, 2), + ) + + # sim_real_imag_train = np.array( + # (gridded_data_train[:, 0] + 1j * gridded_data_train[:, 1]) + # ) + # dirty_image_train = np.abs( + # np.fft.fftshift( + # np.fft.fft2( + # np.fft.fftshift(sim_real_imag_train, axes=(1, 2)), axes=(1, 2) + # ), + # axes=(1, 2), + # ) + # ) + + if conf["amp_phase"]: + gridded_data_train = convert_amp_phase(gridded_data_train) + truth_amp_phase_train = convert_amp_phase(truth_fft_train, sky_sim=True) + + out = out_path / Path("bundle_train" + str(i) + ".h5") + save_fft_pair(out, gridded_data_train, truth_amp_phase_train) + train_index_last = i + + for i in tqdm(range(bundle_valid)): + i += train_index_last + uv_data_valid = np.array( + [ + fits_files.get_uv_data(n) + for n in np.arange( + i * conf["bundle_size"], + (i * conf["bundle_size"]) + conf["bundle_size"], + ) + ], + dtype="object", + ) + gridded_data_valid = np.array( + [grid_data(data, conf) for data in tqdm(uv_data_valid)] + ) + sky_dist_valid = np.array( + [ + cv2.resize( + sky_dist[int(n)][0][0], (conf["grid_size"], conf["grid_size"]) + ) + for n in np.arange( + i * conf["bundle_size"], + (i * conf["bundle_size"]) + conf["bundle_size"], + ) + ] + ) + truth_fft_valid = np.fft.fftshift( + np.fft.fft2(np.fft.fftshift(sky_dist_valid, axes=(1, 2)), axes=(1, 2)), + axes=(1, 2), + ) + + # sim_real_imag_valid = np.array( + # (gridded_data_valid[:, 0] + 1j * gridded_data_valid[:, 1]) + # ) + # dirty_image_valid = np.abs( + # np.fft.fftshift( + # np.fft.fft2( + # np.fft.fftshift(sim_real_imag_valid, axes=(1, 2)), axes=(1, 2) + # ), + # axes=(1, 2), + # ) + # ) + + if conf["amp_phase"]: + gridded_data_valid = convert_amp_phase(gridded_data_valid) + truth_amp_phase_valid = convert_amp_phase(truth_fft_valid, sky_sim=True) + + out = out_path / Path("bundle_valid" + str(i) + ".h5") + save_fft_pair(out, gridded_data_valid, truth_amp_phase_valid) def grid_data(uv_data, conf): @@ -60,8 +148,7 @@ def grid_data(uv_data, conf): real = np.squeeze(cmplx[..., 0, 0]) imag = np.squeeze(cmplx[..., 0, 1]) # weight = np.squeeze(cmplx[..., 0, 2]) - # ap = np.sqrt(real ** 2 + imag ** 2) - # ph = np.angle(real + 1j * imag) + samps = np.array( [ np.append(uv_data["UU--"] * freq, -uv_data["UU--"] * freq), @@ -78,7 +165,9 @@ def grid_data(uv_data, conf): mask = np.zeros((N, N, u_0.shape[0])) fov = ( conf["fov_size"] * np.pi / (3600 * 180) - ) # hard code #default 0.00018382, FoV from VLBA 163.7 <- wrong! depends on setting of simulations + ) # hard code #default 0.00018382, FoV from VLBA 163.7 <- wrong! + # depends on setting of simulations + delta_u = 1 / (fov) for i in range(N): for j in range(N): @@ -93,69 +182,24 @@ def grid_data(uv_data, conf): points[points == 0] = 1 gridded_vis = np.zeros((2, N, N)) gridded_vis[0] = np.matmul(mask, samps[2].T) / points - # gridded_vis[0] = (np.log10(gridded_vis[0] + 1e-10) / 10) + 1 gridded_vis[1] = np.matmul(mask, samps[3].T) / points + return gridded_vis - # print(gridded_vis.shape) - - # plt.imshow(gridded_vis[1]) - # plt.colorbar() - # plt.show() - return gridded_vis +def convert_amp_phase(data, sky_sim=False): + if sky_sim: + amp = np.abs(data) + phase = np.angle(data) + amp = (np.log10(amp + 1e-10) / 10) + 1 + data = np.stack((amp, phase), axis=1) + else: + data[:, 0] = np.sqrt(data[:, 0] ** 2 + data[:, 1] ** 2) + data[:, 1] = np.angle(data[:, 0] + 1j * data[:, 1]) + data[:, 0] = (np.log10(data[:, 0] + 1e-10) / 10) + 1 + return data if __name__ == "__main__": create_gridded_data_set( "/net/big-tank/POOL/projects/radio/test_rime/create_dataset.toml" ) - - # img = np.zeros((size, 256, 256)) - # samps = np.zeros((size, 4, 21000)) # hard code - # for i in tqdm(range(size)): - # sampled = bundle_paths_input[i] - # target = bundle_paths_target[i] - - # img[i] = np.asarray(Image.open(str(target))) - # # img[i] = img[i]/np.sum(img[i]) - - # print(f"\n Gridding VLBI data set.\n") - - # # Generate Mask - # u_0 = samps[0][0] - # v_0 = samps[0][1] - # N = 63 # hard code - # mask = np.zeros((N, N, u_0.shape[0])) - # fov = 0.00018382 * np.pi / (3600 * 180) # hard code #default 0.00018382 - # # delta_u = 1/(fov*N/256) # hard code - # delta_u = 1 / (fov) - # for i in range(N): - # for j in range(N): - # u_cell = (j - N / 2) * delta_u - # v_cell = (i - N / 2) * delta_u - # mask[i, j] = ((u_cell <= u_0) & (u_0 <= u_cell + delta_u)) & ( - # (v_cell <= v_0) & (v_0 <= v_cell + delta_u) - # ) - - # mask = np.flip(mask, [0]) - # points = np.sum(mask, 2) - # points[points == 0] = 1 - # samp_img = np.zeros((size, 2, N, N)) - # img_resized = np.zeros((size, N, N)) - # for i in tqdm(range(samps.shape[0])): - # samp_img[i][0] = np.matmul(mask, samps[i][2].T) / points - # samp_img[i][0] = (np.log10(samp_img[i][0] + 1e-10) / 10) + 1 - # samp_img[i][1] = np.matmul(mask, samps[i][3].T) / points - # img_resized[i] = cv2.resize(img[i], (N, N)) - # img_resized[i] = img_resized[i] / np.sum(img_resized[i]) - - # # truth_fft = np.array([np.fft.fft2(np.fft.fftshift(img)) for im in img_resized]) - # truth_fft = np.fft.fftshift( - # np.fft.fft2(np.fft.fftshift(img_resized, axes=(1, 2)), axes=(1, 2)), axes=(1, 2) - # ) - # fft_scaled_truth = prepare_fft_images(truth_fft, True, False) - - # out = data_path + "/samp_train0.h5" - # save_fft_pair(out, samp_img[:2300], fft_scaled_truth[:2300]) - # out = data_path + "/samp_valid0.h5" - # save_fft_pair(out, samp_img[2300:], fft_scaled_truth[2300:]) \ No newline at end of file From b6ec77f477a15de2062c7cbcf8f990eecaa63e7f Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Mon, 7 Mar 2022 06:29:18 +0100 Subject: [PATCH 40/63] adjust declination angle --- pyvisgen/simulation/data_set.py | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/pyvisgen/simulation/data_set.py b/pyvisgen/simulation/data_set.py index b829c0a..db593d8 100644 --- a/pyvisgen/simulation/data_set.py +++ b/pyvisgen/simulation/data_set.py @@ -58,15 +58,12 @@ def draw_sampling_opts(conf): datetime.time(datetime.strptime(date, "%H:%M:%S")) for date in date_str_ra ] fov_center_ra = np.random.choice(times_ra) - date_str_dec = pd.date_range( - conf["fov_center_dec"][0][0].strftime("%H:%M:%S"), - conf["fov_center_dec"][0][1].strftime("%H:%M:%S"), - freq="1min", - ).strftime("%H:%M:%S") - times_dec = [ - datetime.time(datetime.strptime(date, "%H:%M:%S")) for date in date_str_dec - ] - fov_center_dec = np.random.choice(times_dec) + angles_dec = np.arange( + conf["fov_center_dec"][0][0], + conf["fov_center_dec"][0][1], + step=0.1, + ) + fov_center_dec = np.random.choice(angles_dec) start_time_l = datetime.strptime(conf["scan_start"][0], "%d-%m-%Y %H:%M:%S") start_time_h = datetime.strptime(conf["scan_start"][1], "%d-%m-%Y %H:%M:%S") start_times = pd.date_range( @@ -105,4 +102,4 @@ def draw_sampling_opts(conf): if __name__ == "__main__": - simulate_data_set("../../config/data_set.toml") + simulate_data_set("/net/big-tank/POOL/projects/radio/test_rime/create_dataset.toml") From 736ad35bb0202eefa37418a6dde0e4b36164b88b Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Mon, 7 Mar 2022 06:29:51 +0100 Subject: [PATCH 41/63] direction independent ony for I --- pyvisgen/simulation/scan.py | 59 ++++++++++++++++++++----------------- 1 file changed, 32 insertions(+), 27 deletions(-) diff --git a/pyvisgen/simulation/scan.py b/pyvisgen/simulation/scan.py index 64c24bc..2f96cbc 100644 --- a/pyvisgen/simulation/scan.py +++ b/pyvisgen/simulation/scan.py @@ -296,7 +296,7 @@ def uncorrupted(lm, baselines, wave, time, src_crd, array_layout, SI): return X -def corrupted(lm, baselines, wave, time, src_crd, array_layout, I, rd): +def corrupted(lm, baselines, wave, time, src_crd, array_layout, SI, rd): """Calculates corrupted visibility Parameters @@ -336,15 +336,16 @@ def corrupted(lm, baselines, wave, time, src_crd, array_layout, I, rd): K = getK(baselines, lm, wave, base_num) - B = np.zeros((lm.shape[0], lm.shape[1], 2, 2), dtype=complex) + B = np.zeros((lm.shape[0], lm.shape[1], 1), dtype=complex) - B[:, :, 0, 0] = I[:, :, 0] + I[:, :, 1] - B[:, :, 0, 1] = I[:, :, 2] + 1j * I[:, :, 3] - B[:, :, 1, 0] = I[:, :, 2] - 1j * I[:, :, 3] - B[:, :, 1, 1] = I[:, :, 0] - I[:, :, 1] + B[:, :, 0] = SI + SI + # B[:, :, 0, 0] = I[:, :, 0] + I[:, :, 1] + # B[:, :, 0, 1] = I[:, :, 2] + 1j * I[:, :, 3] + # B[:, :, 1, 0] = I[:, :, 2] - 1j * I[:, :, 3] + # B[:, :, 1, 1] = I[:, :, 0] - I[:, :, 1] # coherency - X = torch.einsum("lmij,lmb->lmbij", torch.tensor(B), K) + X = torch.einsum("lmi,lmb->lmbi", torch.tensor(B), K) # X = np.einsum('lmij,lmb->lmbij', B, K, optimize=True) # X = torch.tensor(B)[:,:,None,:,:] * K[:,:,:,None,None] @@ -357,13 +358,15 @@ def corrupted(lm, baselines, wave, time, src_crd, array_layout, I, rd): E1 = torch.tensor(E_st[:, :, st1_num], dtype=torch.cdouble) E2 = torch.tensor(E_st[:, :, st2_num], dtype=torch.cdouble) - EX = torch.einsum("lmb,lmbij->lmbij", E1, X) + EX = torch.einsum("lmb,lmbi->lmbi", E1, X) del E1, X # EXE = torch.einsum('lmbij,lmbjk->lmbik',EX,torch.transpose(torch.conj(E2),3,4)) - EXE = torch.einsum("lmbij,lmb->lmbij", EX, E2) + EXE = torch.einsum("lmbi,lmb->lmbi", EX, E2) del EX, E2 + # return EXE + # P matrix # parallactic angle @@ -381,17 +384,17 @@ def corrupted(lm, baselines, wave, time, src_crd, array_layout, I, rd): P1 = torch.tensor(getP(b1), dtype=torch.cdouble) P2 = torch.tensor(getP(b2), dtype=torch.cdouble) - PEXE = torch.einsum("bij,lmbjk->lmbik", P1, EXE) + print("P", P1.shape) + + PEXE = torch.einsum("bi,lmbj->lmbi", P1, EXE) del EXE - PEXEP = torch.einsum( - "lmbij,bjk->lmbik", PEXE, torch.transpose(torch.conj(P2), 1, 2) - ) + PEXEP = torch.einsum("lmbi,bk->lmbik", PEXE, torch.transpose(torch.conj(P2), 1, 2)) del PEXE return PEXEP -def direction_independent(lm, baselines, wave, time, src_crd, array_layout, I, rd): +def direction_independent(lm, baselines, wave, time, src_crd, array_layout, SI, rd): """Calculates direction independet visibility Parameters @@ -430,15 +433,16 @@ def direction_independent(lm, baselines, wave, time, src_crd, array_layout, I, r K = getK(baselines, lm, wave, base_num) - B = np.zeros((lm.shape[0], lm.shape[1], 2, 2), dtype=complex) + B = np.zeros((lm.shape[0], lm.shape[1], 1), dtype=complex) - B[:, :, 0, 0] = I[:, :, 0] + I[:, :, 1] - B[:, :, 0, 1] = I[:, :, 2] + 1j * I[:, :, 3] - B[:, :, 1, 0] = I[:, :, 2] - 1j * I[:, :, 3] - B[:, :, 1, 1] = I[:, :, 0] - I[:, :, 1] + B[:, :, 0] = SI + SI + # B[:, :, 0, 0] = I[:, :, 0] + I[:, :, 1] + # B[:, :, 0, 1] = I[:, :, 2] + 1j * I[:, :, 3] + # B[:, :, 1, 0] = I[:, :, 2] - 1j * I[:, :, 3] + # B[:, :, 1, 1] = I[:, :, 0] - I[:, :, 1] # coherency - X = torch.einsum("lmij,lmb->lmbij", torch.tensor(B), K) + X = torch.einsum("lmi,lmb->lmbi", torch.tensor(B), K) del K @@ -448,11 +452,11 @@ def direction_independent(lm, baselines, wave, time, src_crd, array_layout, I, r E1 = torch.tensor(E_st[:, :, st1_num], dtype=torch.cdouble) E2 = torch.tensor(E_st[:, :, st2_num], dtype=torch.cdouble) - EX = torch.einsum("lmb,lmbij->lmbij", E1, X) + EX = torch.einsum("lmb,lmbi->lmbi", E1, X) del E1, X - EXE = torch.einsum("lmbij,lmb->lmbij", EX, E2) + EXE = torch.einsum("lmbi,lmb->lmbi", EX, E2) del EX, E2 return EXE @@ -511,6 +515,7 @@ def getE(rd, array_layout, wave, src_crd): # calculate matrix E for every point in grid # E = np.zeros((rd.shape[0], rd.shape[1], array_layout.st_num.shape[0], 2, 2)) E = np.zeros((rd.shape[0], rd.shape[1], array_layout.st_num.shape[0])) + # print("E", E.shape) # get diameters of all stations and do vectorizing stuff diameters = array_layout.diam @@ -565,12 +570,12 @@ def getP(beta): Shape is given by beta axis and (2,2) Jones matrix axes """ # calculate matrix P with parallactic angle beta - P = np.zeros((beta.shape[0], 2, 2)) + P = np.zeros((beta.shape[0], 1)) - P[:, 0, 0] = np.cos(beta) - P[:, 0, 1] = -np.sin(beta) - P[:, 1, 0] = np.sin(beta) - P[:, 1, 1] = np.cos(beta) + P[:, 0] = np.cos(beta) + # P[:, 0, 1] = -np.sin(beta) + # P[:, 1, 0] = np.sin(beta) + # P[:, 1, 1] = np.cos(beta) return P From e4ecf2018228147e433093a477c245222d43ab74 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Mon, 7 Mar 2022 06:30:12 +0100 Subject: [PATCH 42/63] delete print --- pyvisgen/simulation/scripts/create_dataset.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyvisgen/simulation/scripts/create_dataset.py b/pyvisgen/simulation/scripts/create_dataset.py index e664a99..85f4d20 100644 --- a/pyvisgen/simulation/scripts/create_dataset.py +++ b/pyvisgen/simulation/scripts/create_dataset.py @@ -18,7 +18,6 @@ ) def main(configuration_path, mode): if mode == "simulate": - print("hi") simulate_data_set(configuration_path) if mode == "gridding": create_gridded_data_set(configuration_path) From 32cda82d1a6addc37979d844ac1baa9ea9db3660 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Mon, 7 Mar 2022 06:31:01 +0100 Subject: [PATCH 43/63] use new ra dec, direction independent --- pyvisgen/simulation/visibility.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pyvisgen/simulation/visibility.py b/pyvisgen/simulation/visibility.py index 3cd1790..fe1cdf6 100644 --- a/pyvisgen/simulation/visibility.py +++ b/pyvisgen/simulation/visibility.py @@ -81,12 +81,12 @@ def vis_loop(rc, SI, num_threads=48): # read config array_layout = layouts.get_array_layout(rc["layout"]) src_crd = SkyCoord( - ra=170, # rc["fov_center_ra"].strftime("%H:%M:%S"), - dec=22, # rc["fov_center_dec"].strftime("%H:%M:%S"), - unit=(un.deg, un.deg), + ra=rc["fov_center_ra"].strftime("%H:%M:%S"), + dec=rc["fov_center_dec"], + unit=(un.hourangle, un.deg), ) - rc["fov_center_ra"] = 170 - rc["fov_center_dec"] = 22 + rc["fov_center_ra"] = src_crd.ra.value + rc["fov_center_dec"] = src_crd.dec.value wave = np.array( [ @@ -133,10 +133,10 @@ def vis_loop(rc, SI, num_threads=48): _date = np.zeros(len(u_valid)) - X1 = scan.uncorrupted(lm, baselines, wave[0], t, src_crd, array_layout, SI) + X1 = scan.direction_independent(lm, baselines, wave[0], t, src_crd, array_layout, SI, rd) if X1.shape[0] == 1: continue - X2 = scan.uncorrupted(lm, baselines, wave[0], t, src_crd, array_layout, SI) + X2 = scan.direction_independent(lm, baselines, wave[0], t, src_crd, array_layout, SI, rd) vis_num = np.arange(X1.shape[2] // 2) + 1 + vis_num.max() From 0691a42752199f526360aca219902d25c3ed90d2 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Mon, 7 Mar 2022 06:31:24 +0100 Subject: [PATCH 44/63] add amp phase options --- pyvisgen/utils/config.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyvisgen/utils/config.py b/pyvisgen/utils/config.py index 3358a68..3eea4ad 100644 --- a/pyvisgen/utils/config.py +++ b/pyvisgen/utils/config.py @@ -38,6 +38,8 @@ def read_data_set_conf(conf_toml): conf["bundle_size"] = config["bundle_options"]["bundle_size"] conf["train_valid_split"] = config["bundle_options"]["train_valid_split"] conf["grid_size"] = config["bundle_options"]["grid_size"] + conf["amp_phase"] = config["bundle_options"]["amp_phase"] + conf["target"] = config["bundle_options"]["target"] conf["in_path"] = config["bundle_options"]["in_path"] conf["out_path"] = config["bundle_options"]["out_path"] From 51ca6414345cd3210e884e1e1a22a13710ff6189 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Tue, 8 Mar 2022 01:56:01 +0100 Subject: [PATCH 45/63] rename out to samp --- pyvisgen/gridding/gridder.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyvisgen/gridding/gridder.py b/pyvisgen/gridding/gridder.py index d25b2d4..1c50973 100644 --- a/pyvisgen/gridding/gridder.py +++ b/pyvisgen/gridding/gridder.py @@ -86,7 +86,7 @@ def create_gridded_data_set(config): gridded_data_train = convert_amp_phase(gridded_data_train) truth_amp_phase_train = convert_amp_phase(truth_fft_train, sky_sim=True) - out = out_path / Path("bundle_train" + str(i) + ".h5") + out = out_path / Path("samp_train" + str(i) + ".h5") save_fft_pair(out, gridded_data_train, truth_amp_phase_train) train_index_last = i @@ -137,7 +137,7 @@ def create_gridded_data_set(config): gridded_data_valid = convert_amp_phase(gridded_data_valid) truth_amp_phase_valid = convert_amp_phase(truth_fft_valid, sky_sim=True) - out = out_path / Path("bundle_valid" + str(i) + ".h5") + out = out_path / Path("samp_valid" + str(i - train_index_last) + ".h5") save_fft_pair(out, gridded_data_valid, truth_amp_phase_valid) From 7b09dd267ba8594c0b2ddc297278a8d63ee0c816 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Tue, 8 Mar 2022 04:09:51 +0100 Subject: [PATCH 46/63] suppress warnings --- pyvisgen/fits/writer.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyvisgen/fits/writer.py b/pyvisgen/fits/writer.py index 0f97f4b..09aaf7d 100644 --- a/pyvisgen/fits/writer.py +++ b/pyvisgen/fits/writer.py @@ -5,6 +5,7 @@ from astropy.time import Time import astropy.constants as const import pyvisgen.layouts.layouts as layouts +import warnings def create_vis_hdu(data, conf, layout="vlba", source_name="sim-source-0"): @@ -370,6 +371,7 @@ def create_antenna_hdu(conf): def create_hdu_list(data, conf): + warnings.filterwarnings("ignore", module="astropy.io.fits") vis_hdu = create_vis_hdu(data, conf) time_hdu = create_time_hdu(data) freq_hdu = create_frequency_hdu(conf) From 35f4d557e326e3c7b94d0d6eac070894551c6ee1 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Tue, 8 Mar 2022 04:11:37 +0100 Subject: [PATCH 47/63] avoid vis with less than 10k entries --- pyvisgen/simulation/data_set.py | 7 ++++++- pyvisgen/simulation/visibility.py | 10 ++++++++-- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/pyvisgen/simulation/data_set.py b/pyvisgen/simulation/data_set.py index db593d8..fb3e7d9 100644 --- a/pyvisgen/simulation/data_set.py +++ b/pyvisgen/simulation/data_set.py @@ -23,7 +23,12 @@ def simulate_data_set(config): SI = torch.tensor(data[i][0][0], dtype=torch.cdouble) samp_ops = create_sampling_rc(conf) - hdu_list = writer.create_hdu_list(vis_loop(samp_ops, SI), samp_ops) + vis_data = vis_loop(samp_ops, SI) + while vis_data == 0: + print("Draw new config!") + samp_ops = create_sampling_rc(conf) + vis_data = vis_loop(samp_ops, SI) + hdu_list = writer.create_hdu_list(vis_data, samp_ops) hdu_list.writeto(out, overwrite=True) diff --git a/pyvisgen/simulation/visibility.py b/pyvisgen/simulation/visibility.py index fe1cdf6..d90ddb3 100644 --- a/pyvisgen/simulation/visibility.py +++ b/pyvisgen/simulation/visibility.py @@ -133,10 +133,14 @@ def vis_loop(rc, SI, num_threads=48): _date = np.zeros(len(u_valid)) - X1 = scan.direction_independent(lm, baselines, wave[0], t, src_crd, array_layout, SI, rd) + X1 = scan.direction_independent( + lm, baselines, wave[0], t, src_crd, array_layout, SI, rd + ) if X1.shape[0] == 1: continue - X2 = scan.direction_independent(lm, baselines, wave[0], t, src_crd, array_layout, SI, rd) + X2 = scan.direction_independent( + lm, baselines, wave[0], t, src_crd, array_layout, SI, rd + ) vis_num = np.arange(X1.shape[2] // 2) + 1 + vis_num.max() @@ -160,4 +164,6 @@ def vis_loop(rc, SI, num_threads=48): ) visibilities.add(vis) + if visibilities.get_values().shape[1] < 10000: + return 0 return visibilities From 385f605c99c4a71e0c21d58dcc9935ffc84f4957 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Tue, 8 Mar 2022 06:52:56 +0100 Subject: [PATCH 48/63] add slurm options --- pyvisgen/simulation/data_set.py | 28 +++++++++++++------ pyvisgen/simulation/scripts/create_dataset.py | 7 ++++- 2 files changed, 26 insertions(+), 9 deletions(-) diff --git a/pyvisgen/simulation/data_set.py b/pyvisgen/simulation/data_set.py index fb3e7d9..a48f15e 100644 --- a/pyvisgen/simulation/data_set.py +++ b/pyvisgen/simulation/data_set.py @@ -10,27 +10,39 @@ from radiosim.data import radiosim_data -def simulate_data_set(config): - np.random.seed(42) +def simulate_data_set(config, slurm=False, job_id=None, n=None): conf = read_data_set_conf(config) out_path = Path(conf["out_path"]) out_path.mkdir(parents=True, exist_ok=True) - # open image - data = radiosim_data(conf["in_path"]) - for i in tqdm(range(len(data))): - out = out_path / Path("vis_" + str(i) + ".fits") - SI = torch.tensor(data[i][0][0], dtype=torch.cdouble) + if slurm: + job_id = int(job_id + n * 1000) + data = radiosim_data(conf["in_path"]) + out = out_path / Path("vis_" + str(job_id) + ".fits") + SI = torch.tensor(data[job_id][0][0], dtype=torch.cdouble) samp_ops = create_sampling_rc(conf) vis_data = vis_loop(samp_ops, SI) while vis_data == 0: - print("Draw new config!") samp_ops = create_sampling_rc(conf) vis_data = vis_loop(samp_ops, SI) hdu_list = writer.create_hdu_list(vis_data, samp_ops) hdu_list.writeto(out, overwrite=True) + else: + data = radiosim_data(conf["in_path"]) + for i in tqdm(range(len(data))): + out = out_path / Path("vis_" + str(i) + ".fits") + SI = torch.tensor(data[i][0][0], dtype=torch.cdouble) + + samp_ops = create_sampling_rc(conf) + vis_data = vis_loop(samp_ops, SI) + while vis_data == 0: + samp_ops = create_sampling_rc(conf) + vis_data = vis_loop(samp_ops, SI) + hdu_list = writer.create_hdu_list(vis_data, samp_ops) + hdu_list.writeto(out, overwrite=True) + def create_sampling_rc(conf): sampling_opts = draw_sampling_opts(conf) diff --git a/pyvisgen/simulation/scripts/create_dataset.py b/pyvisgen/simulation/scripts/create_dataset.py index 85f4d20..622b7fa 100644 --- a/pyvisgen/simulation/scripts/create_dataset.py +++ b/pyvisgen/simulation/scripts/create_dataset.py @@ -5,20 +5,25 @@ @click.command() @click.argument("configuration_path", type=click.Path(exists=True, dir_okay=False)) +@click.option('--job_id', required=False, type=int, default=None) +@click.option('--n', required=False, type=int, default=None) @click.option( "--mode", type=click.Choice( [ "simulate", "gridding", + "slurm", ], case_sensitive=False, ), default="simulate", ) -def main(configuration_path, mode): +def main(configuration_path, mode, job_id=None, n=None): if mode == "simulate": simulate_data_set(configuration_path) + if mode == "slurm": + simulate_data_set(configuration_path, slurm=True, job_id=job_id, n=n) if mode == "gridding": create_gridded_data_set(configuration_path) From 58640c620a14d1fb2dd5d5bea59f18a04bee41a9 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Tue, 5 Apr 2022 17:11:17 +0200 Subject: [PATCH 49/63] fix index when creating test images --- pyvisgen/gridding/gridder.py | 126 +++++++++++++++++++++++++++-------- 1 file changed, 100 insertions(+), 26 deletions(-) diff --git a/pyvisgen/gridding/gridder.py b/pyvisgen/gridding/gridder.py index 1c50973..590df92 100644 --- a/pyvisgen/gridding/gridder.py +++ b/pyvisgen/gridding/gridder.py @@ -11,7 +11,7 @@ def create_gridded_data_set(config): conf = read_data_set_conf(config) out_path_fits = Path(conf["out_path"]) - out_path = out_path_fits.parent / "gridded_data/" + out_path = out_path_fits.parent / "gridded_data_large/" ###### # out_path = Path( # "/net/big-tank/POOL/projects/radio/test_rime/build_test/gridded_data" @@ -28,12 +28,70 @@ def create_gridded_data_set(config): sky_dist = radiosim_data(conf["in_path"]) fits_files = fits_data(out_path_fits) size = len(fits_files) + print(size) # num_bundles = int(size // conf["bundle_size"]) - # if conf["num_test_images"] > 0: - # bundle_test = int(conf["num_test_images"] // conf["bundle_size"]) - # size -= conf["num_test_images"] - # else: - # bundle_test = 0 + + ################### + # test + if conf["num_test_images"] > 0: + bundle_test = int(conf["num_test_images"] // conf["bundle_size"]) + size -= conf["num_test_images"] + + for i in tqdm(range(bundle_test)): + uv_data_test = np.array( + [ + fits_files.get_uv_data(n) + for n in np.arange( + i * conf["bundle_size"], + (i * conf["bundle_size"]) + conf["bundle_size"], + ) + ], + dtype="object", + ) + gridded_data_test = np.array( + [grid_data(data, conf) for data in tqdm(uv_data_test)] + ) + sky_dist_test = np.array( + [ + cv2.resize( + sky_dist[int(n)][0][0], (conf["grid_size"], conf["grid_size"]) + ) + for n in np.arange( + i * conf["bundle_size"], + (i * conf["bundle_size"]) + conf["bundle_size"], + ) + ] + ) + + # norm = np.sum(np.sum(sky_dist_test, keepdims=True, axis=1), axis=2) + # sky_dist_test = np.expand_dims(sky_dist_test, -1) / norm[:, None, None] + truth_fft_test = np.fft.fftshift( + np.fft.fft2( + np.fft.fftshift(sky_dist_test / 1e3, axes=(1, 2)), axes=(1, 2) + ), + axes=(1, 2), + ) + + # sim_real_imag_test = np.array( + # (gridded_data_test[:, 0] + 1j * gridded_data_test[:, 1]) + # ) + # dirty_image_test = np.abs( + # np.fft.fftshift( + # np.fft.fft2( + # np.fft.fftshift(sim_real_imag_test, axes=(1, 2)), axes=(1, 2) + # ), + # axes=(1, 2), + # ) + # ) + + if conf["amp_phase"]: + gridded_data_test = convert_amp_phase(gridded_data_test) + truth_amp_phase_test = convert_amp_phase(truth_fft_test, sky_sim=True) + + out = out_path / Path("samp_test" + str(i) + ".h5") + save_fft_pair(out, gridded_data_test, truth_amp_phase_test) + # + ######### size_valid = conf["train_valid_split"] * size size_train = size - size_valid @@ -41,6 +99,7 @@ def create_gridded_data_set(config): bundle_valid = int(size_valid // conf["bundle_size"]) for i in tqdm(range(bundle_train)): + i += bundle_test uv_data_train = np.array( [ fits_files.get_uv_data(n) @@ -65,8 +124,13 @@ def create_gridded_data_set(config): ) ] ) + + # norm = np.sum(np.sum(sky_dist_train, keepdims=True, axis=1), axis=2) + # sky_dist_train = np.expand_dims(sky_dist_train, -1) / norm[:, None, None] truth_fft_train = np.fft.fftshift( - np.fft.fft2(np.fft.fftshift(sky_dist_train, axes=(1, 2)), axes=(1, 2)), + np.fft.fft2( + np.fft.fftshift(sky_dist_train / 1e3, axes=(1, 2)), axes=(1, 2) + ), axes=(1, 2), ) @@ -116,8 +180,13 @@ def create_gridded_data_set(config): ) ] ) + + # norm = np.sum(np.sum(sky_dist_valid, keepdims=True, axis=1), axis=2) + # sky_dist_valid = np.expand_dims(sky_dist_valid, -1) / norm[:, None, None] truth_fft_valid = np.fft.fftshift( - np.fft.fft2(np.fft.fftshift(sky_dist_valid, axes=(1, 2)), axes=(1, 2)), + np.fft.fft2( + np.fft.fftshift(sky_dist_valid / 1e3, axes=(1, 2)), axes=(1, 2) + ), axes=(1, 2), ) @@ -159,30 +228,35 @@ def grid_data(uv_data, conf): ) # Generate Mask - u_0 = samps[0] - v_0 = samps[1] - N = conf["grid_size"] # image size, needs to be from config file - mask = np.zeros((N, N, u_0.shape[0])) + N = conf["grid_size"] # image size fov = ( conf["fov_size"] * np.pi / (3600 * 180) ) # hard code #default 0.00018382, FoV from VLBA 163.7 <- wrong! # depends on setting of simulations + delta = 1 / fov - delta_u = 1 / (fov) - for i in range(N): - for j in range(N): - u_cell = (j - N / 2) * delta_u - v_cell = (i - N / 2) * delta_u - mask[i, j] = ((u_cell <= u_0) & (u_0 <= u_cell + delta_u)) & ( - (v_cell <= v_0) & (v_0 <= v_cell + delta_u) - ) + bins = np.arange(start=-(N / 2) * delta, stop=(N / 2 + 1) * delta, step=delta) + + mask, *_ = np.histogram2d(samps[0], samps[1], bins=[bins, bins], normed=False) + mask[mask == 0] = 1 + mask = np.rot90(mask) + + mask_real, x_edges, y_edges = np.histogram2d( + samps[0], samps[1], bins=[bins, bins], weights=samps[2], normed=False + ) + mask_imag, x_edges, y_edges = np.histogram2d( + samps[0], samps[1], bins=[bins, bins], weights=samps[3], normed=False + ) + + mask_real = np.rot90(mask_real) + mask_imag = np.rot90(mask_imag) + + mask_real /= mask + mask_imag /= mask - mask = np.flip(mask, [0]) - points = np.sum(mask, 2) - points[points == 0] = 1 gridded_vis = np.zeros((2, N, N)) - gridded_vis[0] = np.matmul(mask, samps[2].T) / points - gridded_vis[1] = np.matmul(mask, samps[3].T) / points + gridded_vis[0] = mask_real + gridded_vis[1] = mask_imag return gridded_vis @@ -193,7 +267,7 @@ def convert_amp_phase(data, sky_sim=False): amp = (np.log10(amp + 1e-10) / 10) + 1 data = np.stack((amp, phase), axis=1) else: - data[:, 0] = np.sqrt(data[:, 0] ** 2 + data[:, 1] ** 2) + data[:, 0] = np.sqrt(data[:, 0] ** 2 + data[:, 1] ** 2) / 1e3 data[:, 1] = np.angle(data[:, 0] + 1j * data[:, 1]) data[:, 0] = (np.log10(data[:, 0] + 1e-10) / 10) + 1 return data From 40ef4ee37951d3e4717c77650e573024bfd3edbc Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Tue, 5 Apr 2022 17:12:36 +0200 Subject: [PATCH 50/63] add multichannel, use infos from config --- pyvisgen/fits/writer.py | 32 ++++++++++++++++++++++---------- 1 file changed, 22 insertions(+), 10 deletions(-) diff --git a/pyvisgen/fits/writer.py b/pyvisgen/fits/writer.py index 09aaf7d..8329dc9 100644 --- a/pyvisgen/fits/writer.py +++ b/pyvisgen/fits/writer.py @@ -28,17 +28,27 @@ def create_vis_hdu(data, conf, layout="vlba", source_name="sim-source-0"): # visibility data values = data.get_values() + print(values.shape) + # print(values.real.shape) + # stack = np.stack([values.real, values.imag, np.ones(values.shape)], axis=2) + # print("stack", stack.shape) + # print(stack[:, :, :, 0:2]) + + num_ifs = values.shape[0] + vis = np.swapaxes( np.swapaxes( - np.stack([values.real, values.imag, np.ones(values.shape) - 0.8], axis=1), - 1, - 2, + np.stack([values.real, values.imag, np.ones(values.shape)], axis=2), + 0, + 3, ), - 0, - 1, - ).reshape(-1, 1, 1, 1, 1, 4, 3) + 2, + 3, + ).reshape(-1, 1, 1, num_ifs, 1, 4, 3) DATA = vis # in dim 4 = IFs , dim = 1, dim 4 = number of jones, 3 = real, imag, weight + print(DATA) + print(DATA.shape) # wcs ra = conf["fov_center_ra"] @@ -179,24 +189,26 @@ def create_time_hdu(data): def create_frequency_hdu(conf): - freq_d = (conf["bandwidths"][0] * un.Hz).value + # freq_d = (conf["bandwidths"][0] * un.Hz).value # num_ifs = 1 # at the moment only 1 possible FRQSEL = np.array([1], dtype=">i4") col1 = fits.Column(name="FRQSEL", format="1J", unit=" ", array=FRQSEL) - IF_FREQ = np.array([[0.00e00]], dtype=">f8") # start with 0, add ch_with per IF + IF_FREQ = np.array( + [conf["frequsel"]], dtype=">f8" + ) # start with 0, add ch_with per IF col2 = fits.Column( name="IF FREQ", format=str(IF_FREQ.shape[-1]) + "D", unit="Hz", array=IF_FREQ ) - CH_WIDTH = np.repeat(np.array([[freq_d]], dtype=">f4"), 1, axis=1) + CH_WIDTH = np.repeat(np.array([conf["bandwidths"]], dtype=">f4"), 1, axis=1) col3 = fits.Column( name="CH WIDTH", format=str(CH_WIDTH.shape[-1]) + "E", unit="Hz", array=CH_WIDTH ) - TOTAL_BANDWIDTH = np.repeat(np.array([[freq_d]], dtype=">f4"), 1, axis=1) + TOTAL_BANDWIDTH = np.repeat(np.array([conf["bandwidths"]], dtype=">f4"), 1, axis=1) col4 = fits.Column( name="TOTAL BANDWIDTH", format=str(TOTAL_BANDWIDTH.shape[-1]) + "E", From e43a36bb51bbabb31773a1a48594225ab470ac9b Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Tue, 5 Apr 2022 17:13:06 +0200 Subject: [PATCH 51/63] change random seed --- pyvisgen/simulation/data_set.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyvisgen/simulation/data_set.py b/pyvisgen/simulation/data_set.py index fb3e7d9..5d3fcbd 100644 --- a/pyvisgen/simulation/data_set.py +++ b/pyvisgen/simulation/data_set.py @@ -11,7 +11,7 @@ def simulate_data_set(config): - np.random.seed(42) + np.random.seed(1) conf = read_data_set_conf(config) out_path = Path(conf["out_path"]) out_path.mkdir(parents=True, exist_ok=True) From 8661081b7bef06c50ed23853318405a01b0d3b09 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Tue, 5 Apr 2022 17:13:45 +0200 Subject: [PATCH 52/63] delete unused imports --- pyvisgen/simulation/scan.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/pyvisgen/simulation/scan.py b/pyvisgen/simulation/scan.py index 2f96cbc..6ac7d18 100644 --- a/pyvisgen/simulation/scan.py +++ b/pyvisgen/simulation/scan.py @@ -1,17 +1,13 @@ from dataclasses import dataclass from astropy import units as un -from astropy.coordinates import SkyCoord, EarthLocation, AltAz, Angle +from astropy.coordinates import EarthLocation, AltAz, Angle import numpy as np from scipy.special import j1 -import scipy.constants as const -import scipy.signal as sig from astroplan import Observer from pyvisgen.simulation.utils import single_occurance, get_pairs -from pyvisgen.layouts import layouts import torch import itertools -import time as t -import numexpr as ne # fast exponential +import numexpr as ne @dataclass From a8d5ecfb4b80f2d4cd104f218eb01f767f33dd07 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Tue, 5 Apr 2022 17:14:31 +0200 Subject: [PATCH 53/63] add multichannel --- pyvisgen/simulation/visibility.py | 77 +++++++++++++++++++++++-------- 1 file changed, 57 insertions(+), 20 deletions(-) diff --git a/pyvisgen/simulation/visibility.py b/pyvisgen/simulation/visibility.py index d90ddb3..987a92b 100644 --- a/pyvisgen/simulation/visibility.py +++ b/pyvisgen/simulation/visibility.py @@ -88,7 +88,7 @@ def vis_loop(rc, SI, num_threads=48): rc["fov_center_ra"] = src_crd.ra.value rc["fov_center_dec"] = src_crd.dec.value - wave = np.array( + waves = np.array( [ const.c / ((rc["base_freq"] + float(freq)) / un.second) / un.meter for freq in rc["frequsel"] @@ -107,7 +107,20 @@ def vis_loop(rc, SI, num_threads=48): base_num = int(stat_num * (stat_num - 1) / 2) # calculate vis - visibilities = Visibilities([], [], [], [], [], [], [], [], [], [], [], []) + visibilities = Visibilities( + np.empty(shape=[0] + [len(waves)]), + np.empty(shape=[0] + [len(waves)]), + np.empty(shape=[0] + [len(waves)]), + np.empty(shape=[0] + [len(waves)]), + [], + [], + [], + [], + [], + [], + [], + [], + ) vis_num = np.zeros(1) # i in total number of scans for i in range(rc["scans"]): @@ -133,26 +146,34 @@ def vis_loop(rc, SI, num_threads=48): _date = np.zeros(len(u_valid)) - X1 = scan.direction_independent( - lm, baselines, wave[0], t, src_crd, array_layout, SI, rd + int_values = np.array( + [ + calc_vis( + lm, + baselines, + wave, + t, + src_crd, + array_layout, + SI, + rd, + vis_num, + ) + for wave in waves + ] ) - if X1.shape[0] == 1: + if int_values.shape[1] == 1: continue - X2 = scan.direction_independent( - lm, baselines, wave[0], t, src_crd, array_layout, SI, rd - ) - - vis_num = np.arange(X1.shape[2] // 2) + 1 + vis_num.max() - - int_values = scan.integrate(X1, X2) - del X1, X2 - int_values = int_values + int_values = np.swapaxes(int_values, 0, 1) + vis_num = np.arange(int_values.shape[0]) + 1 + vis_num.max() + print(int_values.shape) + print(int_values[:, :, 0].shape) vis = Visibilities( - int_values[:, 0], - torch.zeros(int_values.shape[0], dtype=torch.complex128), - torch.zeros(int_values.shape[0], dtype=torch.complex128), - torch.zeros(int_values.shape[0], dtype=torch.complex128), + torch.tensor(int_values[:, :, 0]), + torch.zeros(int_values[:, :, 0].shape, dtype=torch.complex128), + torch.zeros(int_values[:, :, 0].shape, dtype=torch.complex128), + torch.zeros(int_values[:, :, 0].shape, dtype=torch.complex128), vis_num, np.repeat(i + 1, len(vis_num)), np.array([baselines[i].baselineNum() for i in base_valid]), @@ -164,6 +185,22 @@ def vis_loop(rc, SI, num_threads=48): ) visibilities.add(vis) - if visibilities.get_values().shape[1] < 10000: - return 0 + # if visibilities.get_values().shape[1] < 10000: + # return 0 return visibilities + + +def calc_vis(lm, baselines, wave, t, src_crd, array_layout, SI, rd, vis_num): + X1 = scan.direction_independent( + lm, baselines, wave, t, src_crd, array_layout, SI, rd + ) + if X1.shape[0] == 1: + return -1 + X2 = scan.direction_independent( + lm, baselines, wave, t, src_crd, array_layout, SI, rd + ) + + int_values = scan.integrate(X1, X2).numpy() + del X1, X2 + int_values = int_values + return int_values From 99e35f0891c9dfb52a9274992a46698a0551209f Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Tue, 5 Apr 2022 17:54:06 +0200 Subject: [PATCH 54/63] fix tests --- config/default.toml | 2 +- tests/test_utils.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/config/default.toml b/config/default.toml index aceb24c..f13052b 100644 --- a/config/default.toml +++ b/config/default.toml @@ -3,7 +3,7 @@ fov_center_ra = "12:30:49.423382" fov_center_dec = "12:23:28.04366" fov_size = 0.00018382 corr_int_time = 10.0 -scan_start = "2016:95:00:00:00" +scan_start = "22-05-2021 04:00:01" scan_duration = 300 scans = 72 interval_length = 1200 diff --git a/tests/test_utils.py b/tests/test_utils.py index 5305a81..4e9cd4a 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -47,4 +47,4 @@ def test_calc_time_steps(): conf = read_config("config/default.toml") time = calc_time_steps(conf) - assert time.shape == (2232,) \ No newline at end of file + assert time.shape == (2232,) From 7de2537a9d1f9e43b3d6b2e05feba842b9b01a31 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Tue, 5 Apr 2022 17:54:23 +0200 Subject: [PATCH 55/63] fix tests --- pyvisgen/simulation/utils.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/pyvisgen/simulation/utils.py b/pyvisgen/simulation/utils.py index 7ab4193..cf591e9 100644 --- a/pyvisgen/simulation/utils.py +++ b/pyvisgen/simulation/utils.py @@ -3,6 +3,7 @@ import astropy.units as un import numpy as np from astropy.time import Time +from datetime import datetime def read_config(conf): @@ -28,7 +29,9 @@ def read_config(conf): ) sim_conf["fov_size"] = config["sampling_options"]["fov_size"] sim_conf["corr_int_time"] = config["sampling_options"]["corr_int_time"] - sim_conf["scan_start"] = config["sampling_options"]["scan_start"] + sim_conf["scan_start"] = datetime.strptime( + config["sampling_options"]["scan_start"], "%d-%m-%Y %H:%M:%S" + ) sim_conf["scan_duration"] = config["sampling_options"]["scan_duration"] sim_conf["scans"] = config["sampling_options"]["scans"] sim_conf["channel"] = config["sampling_options"]["channel"] @@ -67,7 +70,7 @@ def get_pairs(array_layout): def calc_time_steps(conf): - start_time = Time(conf["scan_start"].isoformat(), format='isot') + start_time = Time(conf["scan_start"].isoformat(), format="isot") interval = conf["interval_length"] integration_time = conf["corr_int_time"] num_scans = conf["scans"] @@ -77,10 +80,10 @@ def calc_time_steps(conf): time_lst = [ start_time + interval * i * un.second + j * integration_time * un.second for i in range(num_scans) - for j in range( - int(scan_duration / int_time) + 1 - ) # +1 because t_1 is the stop time of t_0. in order to save computing power we take one time more to complete interval + for j in range(int(scan_duration / int_time) + 1) ] + # +1 because t_1 is the stop time of t_0 + # in order to save computing power we take one time more to complete interval time = Time(time_lst) return time From 6a1d57a7424330925c465559b145426d81ee6484 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Tue, 5 Apr 2022 17:54:41 +0200 Subject: [PATCH 56/63] flake8 changes --- pyvisgen/simulation/scan.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/pyvisgen/simulation/scan.py b/pyvisgen/simulation/scan.py index 6ac7d18..643b746 100644 --- a/pyvisgen/simulation/scan.py +++ b/pyvisgen/simulation/scan.py @@ -650,7 +650,8 @@ def jinc(x): def get_valid_baselines(baselines, base_num): - """Calculates all valid baselines. This depens on the baselines that are visible at start and stop times + """Calculates all valid baselines. This depens on the baselines that are visible at + start and stop times. Parameters ---------- @@ -668,7 +669,8 @@ def get_valid_baselines(baselines, base_num): valid = baselines.valid.reshape(-1, base_num) # generate a mask to only take baselines that are visible at start and stop time - # example: telescope is visible at time t_0 but not visible at time t_1, therefore throw away baseline + # example: telescope is visible at time t_0 but not visible at time t_1, therefore + # throw away baseline # this is checked for every pair of time: t_0-t_1, t_1-t_2,... # t_0<-mask[0]->t_1, t_1<-mask[1]->t_2,... mask = np.array(valid[:-1]).astype(bool) & np.array(valid[1:]).astype(bool) @@ -704,13 +706,15 @@ def time_step_of_baseline(baselines, base_num): Returns ------- 1d array - Return array with every time step repeated N times, where N is the number of valid baselines per time step + Return array with every time step repeated N times, where N is the number of + valid baselines per time step """ # reshape valid mask to (time, total baselines per time) valid = baselines.valid.reshape(-1, base_num) # generate a mask to only take baselines that are visible at start and stop time - # example: telescope is visible at time t_0 but not visible at time t_1, therefore throw away baseline + # example: telescope is visible at time t_0 but not visible at time t_1, therefore + # throw away baseline # this is checked for every pair of time: t_0-t_1, t_1-t_2,... # t_0<-mask[0]->t_1, t_1<-mask[1]->t_2,... mask = np.array(valid[:-1]).astype(bool) & np.array(valid[1:]).astype(bool) From 5b2b2c86c2b03c1ae36f3e88481a423c0e56b229 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Mon, 8 Aug 2022 16:46:45 +0200 Subject: [PATCH 57/63] close h5 files --- pyvisgen/fits/data.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pyvisgen/fits/data.py b/pyvisgen/fits/data.py index dd8c8de..83acb8c 100644 --- a/pyvisgen/fits/data.py +++ b/pyvisgen/fits/data.py @@ -35,12 +35,14 @@ def get_files(self, path): def get_uv_data(self, i): with fits.open(self.files[i]) as hdul: uv_data = hdul[0].data + hdul.close() return uv_data def get_freq_data(self, i): with fits.open(self.files[i]) as hdul: base_freq = hdul[0].header["CRVAL4"] freq_data = hdul[2].data + hdul.close() return freq_data, base_freq def open_file(self, i): @@ -50,4 +52,5 @@ def open_file(self, i): with fits.open(self.files[i]) as hdul: fits_file = hdul + hdul.close() return fits_file.info() From 978cee1beba98ee4994b4ac1bef1d90bf7e459e6 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Mon, 8 Aug 2022 16:48:36 +0200 Subject: [PATCH 58/63] implement multichannel --- pyvisgen/fits/writer.py | 22 +++++----------------- 1 file changed, 5 insertions(+), 17 deletions(-) diff --git a/pyvisgen/fits/writer.py b/pyvisgen/fits/writer.py index 8329dc9..9f04990 100644 --- a/pyvisgen/fits/writer.py +++ b/pyvisgen/fits/writer.py @@ -17,7 +17,7 @@ def create_vis_hdu(data, conf, layout="vlba", source_name="sim-source-0"): DATE = data.date - int( data.date.min() - ) # placeholder, julian date of vis, central time in the integration period + ) _DATE = data._date # central time in the integration period @@ -26,20 +26,14 @@ def create_vis_hdu(data, conf, layout="vlba", source_name="sim-source-0"): INTTIM = np.repeat(np.array(conf["corr_int_time"], dtype=">f4"), len(u)) # visibility data - values = data.get_values() + values = np.swapaxes(data.get_values(), 0, 1) - print(values.shape) - # print(values.real.shape) - # stack = np.stack([values.real, values.imag, np.ones(values.shape)], axis=2) - # print("stack", stack.shape) - # print(stack[:, :, :, 0:2]) - - num_ifs = values.shape[0] + num_ifs = values.shape[1] vis = np.swapaxes( np.swapaxes( np.stack([values.real, values.imag, np.ones(values.shape)], axis=2), - 0, + 1, 3, ), 2, @@ -47,8 +41,6 @@ def create_vis_hdu(data, conf, layout="vlba", source_name="sim-source-0"): ).reshape(-1, 1, 1, num_ifs, 1, 4, 3) DATA = vis # in dim 4 = IFs , dim = 1, dim 4 = number of jones, 3 = real, imag, weight - print(DATA) - print(DATA.shape) # wcs ra = conf["fov_center_ra"] @@ -63,7 +55,7 @@ def create_vis_hdu(data, conf, layout="vlba", source_name="sim-source-0"): ws.wcs.ctype = ["", "COMPLEX", "STOKES", "FREQ", "IF", "RA", "DEC"] h = ws.to_header() - scale = 1 # / freq + scale = 1 u_scale = u / const.c v_scale = v / const.c w_scale = w / const.c @@ -189,10 +181,6 @@ def create_time_hdu(data): def create_frequency_hdu(conf): - # freq_d = (conf["bandwidths"][0] * un.Hz).value - - # num_ifs = 1 # at the moment only 1 possible - FRQSEL = np.array([1], dtype=">i4") col1 = fits.Column(name="FRQSEL", format="1J", unit=" ", array=FRQSEL) From 9f950e2c971dc51cbd36cb04b05880b23136e128 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Mon, 8 Aug 2022 16:53:30 +0200 Subject: [PATCH 59/63] define basic and ducc0 based gridder --- pyvisgen/gridding/alt_gridder.py | 124 ++++++++++++ pyvisgen/gridding/gridder.py | 315 ++++++++++++++++--------------- 2 files changed, 288 insertions(+), 151 deletions(-) create mode 100644 pyvisgen/gridding/alt_gridder.py diff --git a/pyvisgen/gridding/alt_gridder.py b/pyvisgen/gridding/alt_gridder.py new file mode 100644 index 0000000..7f95954 --- /dev/null +++ b/pyvisgen/gridding/alt_gridder.py @@ -0,0 +1,124 @@ +from scipy.special.orthogonal import p_roots +import scipy +import numpy as np +import astropy.constants as const + +speedoflight = const.c.value + + +def init( + uvw, freq, nxdirty, nydirty, pixsizex, pixsizey, epsilon, do_wgridding, ms=None +): + ofac = 2 + x, y = np.meshgrid( + *[-ss / 2 + np.arange(ss) for ss in (nxdirty, nydirty)], indexing="ij" + ) + x *= pixsizex + y *= pixsizey + eps = x**2 + y**2 + nm1 = -eps / (np.sqrt(1.0 - eps) + 1.0) + ng = ofac * nxdirty, ofac * nydirty + supp = int(np.ceil(np.log10(1 / epsilon * (3 if do_wgridding else 2)))) + 1 + kernel = es_kernel(supp, 2.3 * supp) + uvw = np.transpose((uvw[..., None] * freq / speedoflight), (0, 2, 1)).reshape(-1, 3) + conjind = uvw[:, 2] < 0 + uvw[conjind] *= -1 + u, v, w = uvw.T + if do_wgridding: + wmin, wmax = np.min(w), np.max(w) + dw = 1 / ofac / np.max(np.abs(nm1)) / 2 + nwplanes = int(np.ceil((wmax - wmin) / dw + supp)) if do_wgridding else 1 + w0 = (wmin + wmax) / 2 - dw * (nwplanes - 1) / 2 + else: + nwplanes, w0, dw = 1, None, None + gridcoord = [np.linspace(-0.5, 0.5, nn, endpoint=False) for nn in ng] + slc0, slc1 = slice(nxdirty // 2, nxdirty * 3 // 2), slice( + nydirty // 2, nydirty * 3 // 2 + ) + u *= pixsizex + v *= pixsizey + if ms is not None: + ms[conjind] = ms[conjind].conjugate() + return u, v, w, w0, dw, nwplanes, nm1, kernel, gridcoord, slc0, slc1, ng, ms + return u, v, w, w0, dw, nwplanes, nm1, kernel, gridcoord, slc0, slc1, ng, conjind + + +class Kernel: + def __init__(self, supp, func): + self._func = func + self._supp = supp + + def ft(self, x): + x = x * self._supp * np.pi + nroots = 2 * self._supp + if self._supp % 2 == 0: + nroots += 1 + q, weights = p_roots(nroots) + ind = q > 0 + weights = weights[ind] + q = q[ind] + kq = np.outer(x, q) if len(x.shape) == 1 else np.einsum("ij,k->ijk", x, q) + arr = np.sum(weights * self._raw(q) * np.cos(kq), axis=-1) + return self._supp * arr + + def __call__(self, x): + return self._raw(x / self._supp * 2) + + def _raw(self, x): + ind = np.logical_and(x <= 1, x >= -1) + res = np.zeros_like(x) + res[ind] = self._func(x[ind]) + return res + + +def es_kernel(supp, beta): + return Kernel(supp, lambda x: np.exp(beta * (pow((1 - x) * (1 + x), 0.5) - 1))) + + +def ms2dirty_python_fast( + uvw, freq, vis, nxdirty, nydirty, pixsizex, pixsizey, epsilon, do_wgridding +): + u, v, w, w0, dw, nwplanes, nm1, kernel, gridcoord, slc0, slc1, ng, vis = init( + uvw, freq, nxdirty, nydirty, pixsizex, pixsizey, epsilon, do_wgridding, vis + ) + im = np.zeros((nxdirty, nydirty)) + supp = kernel._supp + for ii in range(nwplanes): + grid = np.zeros(ng, dtype=vis.dtype) + for uu, vv, vis in zip( + u, v, vis * kernel(ii - (w - w0) / dw) if do_wgridding else vis + ): + if vis == 0: + continue + ratposx = (uu * ng[0]) % ng[0] + ratposy = (vv * ng[1]) % ng[1] + xle = int(np.round(ratposx)) - supp // 2 + yle = int(np.round(ratposy)) - supp // 2 + # pos = np.arange(0, supp) + # xkernel = kernel(pos - ratposx + xle) + # ykernel = kernel(pos - ratposy + yle) + for xx in range(supp): + foo = vis #* xkernel[xx] + myxpos = (xle + xx) % ng[0] + for yy in range(supp): + myypos = (yle + yy) % ng[1] + grid[myxpos, myypos] += foo #* ykernel[yy] + # loopim = np.fft.fftshift(np.fft.ifft2(grid)*np.prod(ng)) + # loopim = loopim[slc0, slc1] + # if do_wgridding: + # loopim *= np.exp(-2j*np.pi*nm1*(w0+ii*dw)) + # im += loopim.real + # im /= kernel.ft(gridcoord[0][slc0])[:, None] + # im /= kernel.ft(gridcoord[1][slc1]) + # if do_wgridding: + # im /= (nm1+1)*kernel.ft(nm1*dw) + return grid + + +def get_npixdirty(uvw, freq, fov_deg, mask): + speedOfLight = 299792458. + bl = np.sqrt(uvw[:,0]**2+uvw[:,1]**2+uvw[:,2]**2) + bluvw = bl.reshape((-1,1))*freq.reshape((1,-1))/speedOfLight + maxbluvw = np.max(bluvw*mask) + minsize = int((2*fov_deg*np.pi/180*maxbluvw)) + 1 + return minsize+(minsize%2) # make even diff --git a/pyvisgen/gridding/gridder.py b/pyvisgen/gridding/gridder.py index 590df92..b6033f5 100644 --- a/pyvisgen/gridding/gridder.py +++ b/pyvisgen/gridding/gridder.py @@ -6,30 +6,19 @@ from pyvisgen.fits.data import fits_data from pyvisgen.utils.config import read_data_set_conf from radionets.dl_framework.data import save_fft_pair +import astropy.constants as const +from pyvisgen.gridding.alt_gridder import ms2dirty_python_fast, get_npixdirty def create_gridded_data_set(config): conf = read_data_set_conf(config) out_path_fits = Path(conf["out_path"]) - out_path = out_path_fits.parent / "gridded_data_large/" - ###### - # out_path = Path( - # "/net/big-tank/POOL/projects/radio/test_rime/build_test/gridded_data" - # ) - # out_path_fits = Path( - # "/net/big-tank/POOL/projects/radio/test_rime/build_test/uvfits" - # ) - ###### + out_path = out_path_fits.parent / "gridded/" out_path.mkdir(parents=True, exist_ok=True) - # conf[ - # "in_path" - # ] = "/net/big-tank/POOL/projects/radio/test_rime/build_test/sky_simulations" sky_dist = radiosim_data(conf["in_path"]) fits_files = fits_data(out_path_fits) size = len(fits_files) - print(size) - # num_bundles = int(size // conf["bundle_size"]) ################### # test @@ -38,101 +27,40 @@ def create_gridded_data_set(config): size -= conf["num_test_images"] for i in tqdm(range(bundle_test)): - uv_data_test = np.array( - [ - fits_files.get_uv_data(n) - for n in np.arange( - i * conf["bundle_size"], - (i * conf["bundle_size"]) + conf["bundle_size"], - ) - ], - dtype="object", - ) - gridded_data_test = np.array( - [grid_data(data, conf) for data in tqdm(uv_data_test)] - ) - sky_dist_test = np.array( - [ - cv2.resize( - sky_dist[int(n)][0][0], (conf["grid_size"], conf["grid_size"]) - ) - for n in np.arange( - i * conf["bundle_size"], - (i * conf["bundle_size"]) + conf["bundle_size"], - ) - ] - ) + ( + uv_data_test, + freq_data_test, + gridded_data_test, + sky_dist_test, + ) = open_data(fits_files, sky_dist, conf, i) - # norm = np.sum(np.sum(sky_dist_test, keepdims=True, axis=1), axis=2) - # sky_dist_test = np.expand_dims(sky_dist_test, -1) / norm[:, None, None] - truth_fft_test = np.fft.fftshift( - np.fft.fft2( - np.fft.fftshift(sky_dist_test / 1e3, axes=(1, 2)), axes=(1, 2) - ), - axes=(1, 2), - ) - - # sim_real_imag_test = np.array( - # (gridded_data_test[:, 0] + 1j * gridded_data_test[:, 1]) - # ) - # dirty_image_test = np.abs( - # np.fft.fftshift( - # np.fft.fft2( - # np.fft.fftshift(sim_real_imag_test, axes=(1, 2)), axes=(1, 2) - # ), - # axes=(1, 2), - # ) - # ) + truth_fft_test = calc_truth_fft(sky_dist_test) if conf["amp_phase"]: - gridded_data_test = convert_amp_phase(gridded_data_test) + gridded_data_test = convert_amp_phase(gridded_data_test, sky_sim=False) truth_amp_phase_test = convert_amp_phase(truth_fft_test, sky_sim=True) + assert gridded_data_test.shape[1] == 2 + out = out_path / Path("samp_test" + str(i) + ".h5") save_fft_pair(out, gridded_data_test, truth_amp_phase_test) - # - ######### + # + ################### size_valid = conf["train_valid_split"] * size size_train = size - size_valid bundle_train = int(size_train // conf["bundle_size"]) bundle_valid = int(size_valid // conf["bundle_size"]) + ################### + # train for i in tqdm(range(bundle_train)): i += bundle_test - uv_data_train = np.array( - [ - fits_files.get_uv_data(n) - for n in np.arange( - i * conf["bundle_size"], - (i * conf["bundle_size"]) + conf["bundle_size"], - ) - ], - dtype="object", - ) - gridded_data_train = np.array( - [grid_data(data, conf) for data in tqdm(uv_data_train)] - ) - sky_dist_train = np.array( - [ - cv2.resize( - sky_dist[int(n)][0][0], (conf["grid_size"], conf["grid_size"]) - ) - for n in np.arange( - i * conf["bundle_size"], - (i * conf["bundle_size"]) + conf["bundle_size"], - ) - ] + uv_data_train, freq_data_train, gridded_data_train, sky_dist_train = open_data( + fits_files, sky_dist, conf, i ) - # norm = np.sum(np.sum(sky_dist_train, keepdims=True, axis=1), axis=2) - # sky_dist_train = np.expand_dims(sky_dist_train, -1) / norm[:, None, None] - truth_fft_train = np.fft.fftshift( - np.fft.fft2( - np.fft.fftshift(sky_dist_train / 1e3, axes=(1, 2)), axes=(1, 2) - ), - axes=(1, 2), - ) + truth_fft_train = calc_truth_fft(sky_dist_train) # sim_real_imag_train = np.array( # (gridded_data_train[:, 0] + 1j * gridded_data_train[:, 1]) @@ -147,81 +75,156 @@ def create_gridded_data_set(config): # ) if conf["amp_phase"]: - gridded_data_train = convert_amp_phase(gridded_data_train) + gridded_data_train = convert_amp_phase(gridded_data_train, sky_sim=False) truth_amp_phase_train = convert_amp_phase(truth_fft_train, sky_sim=True) out = out_path / Path("samp_train" + str(i) + ".h5") save_fft_pair(out, gridded_data_train, truth_amp_phase_train) train_index_last = i + # + ################### + ################### + # valid for i in tqdm(range(bundle_valid)): i += train_index_last - uv_data_valid = np.array( - [ - fits_files.get_uv_data(n) - for n in np.arange( - i * conf["bundle_size"], - (i * conf["bundle_size"]) + conf["bundle_size"], - ) - ], - dtype="object", - ) - gridded_data_valid = np.array( - [grid_data(data, conf) for data in tqdm(uv_data_valid)] - ) - sky_dist_valid = np.array( - [ - cv2.resize( - sky_dist[int(n)][0][0], (conf["grid_size"], conf["grid_size"]) - ) - for n in np.arange( - i * conf["bundle_size"], - (i * conf["bundle_size"]) + conf["bundle_size"], - ) - ] + uv_data_valid, freq_data_valid, gridded_data_valid, sky_dist_valid = open_data( + fits_files, sky_dist, conf, i ) - # norm = np.sum(np.sum(sky_dist_valid, keepdims=True, axis=1), axis=2) - # sky_dist_valid = np.expand_dims(sky_dist_valid, -1) / norm[:, None, None] - truth_fft_valid = np.fft.fftshift( - np.fft.fft2( - np.fft.fftshift(sky_dist_valid / 1e3, axes=(1, 2)), axes=(1, 2) - ), - axes=(1, 2), - ) - - # sim_real_imag_valid = np.array( - # (gridded_data_valid[:, 0] + 1j * gridded_data_valid[:, 1]) - # ) - # dirty_image_valid = np.abs( - # np.fft.fftshift( - # np.fft.fft2( - # np.fft.fftshift(sim_real_imag_valid, axes=(1, 2)), axes=(1, 2) - # ), - # axes=(1, 2), - # ) - # ) + truth_fft_valid = calc_truth_fft(sky_dist_valid) if conf["amp_phase"]: - gridded_data_valid = convert_amp_phase(gridded_data_valid) + gridded_data_valid = convert_amp_phase(gridded_data_valid, sky_sim=False) truth_amp_phase_valid = convert_amp_phase(truth_fft_valid, sky_sim=True) out = out_path / Path("samp_valid" + str(i - train_index_last) + ".h5") save_fft_pair(out, gridded_data_valid, truth_amp_phase_valid) + # + ################### + + +def open_data(fits_files, sky_dist, conf, i): + uv_data = np.array( + [ + fits_files.get_uv_data(n).copy() + for n in np.arange( + i * conf["bundle_size"], + (i * conf["bundle_size"]) + conf["bundle_size"], + ) + ], + dtype="object", + ) + freq_data = np.array( + [ + fits_files.get_freq_data(n) + for n in np.arange( + i * conf["bundle_size"], + (i * conf["bundle_size"]) + conf["bundle_size"], + ) + ], + dtype="object", + ) + gridded_data = np.array( + [ + ducc0_gridding(data, freq).copy() + for data, freq in tqdm(zip(uv_data, freq_data)) + ] + ) + gridded_truth = np.array( + [ + sky_dist[int(n)][0][0] + for n in np.arange( + i * conf["bundle_size"], + (i * conf["bundle_size"]) + conf["bundle_size"], + ) + ] + ) + return uv_data, freq_data, gridded_data, gridded_truth + +def calc_truth_fft(sky_dist): + # norm = np.sum(np.sum(sky_dist_test, keepdims=True, axis=1), axis=2) + # sky_dist_test = np.expand_dims(sky_dist_test, -1) / norm[:, None, None] + truth_fft = np.fft.fftshift( + np.fft.fft2(np.fft.fftshift(sky_dist, axes=(1, 2)), axes=(1, 2)), + axes=(1, 2), + ) + return truth_fft + + +def ducc0_gridding(uv_data, freq_data): + vis_ = uv_data["DATA"] + vis = np.array([vis_[:, 0, 0, 0, 0, 0, 0] + 1j * vis_[:, 0, 0, 0, 0, 0, 1]]).T + vis_compl = np.array([vis_[:, 0, 0, 0, 0, 0, 0] + 1j * vis_[:, 0, 0, 0, 0, 0, 1]]).T + uu = np.array(uv_data["UU--"], dtype=np.float64) + uu_compl = np.array(-uv_data["UU--"], dtype=np.float64) + vv = np.array(uv_data["VV--"], dtype=np.float64) + vv_compl = np.array(-uv_data["VV--"], dtype=np.float64) + ww = np.array(uv_data["WW--"], dtype=np.float64) + ww_compl = np.array(uv_data["WW--"], dtype=np.float64) + uvw = np.stack([uu, vv, ww]).T + uvw_compl = np.stack([uu_compl, vv_compl, ww_compl]).T + uvw *= const.c.value + uvw_compl *= const.c.value + # complex conjugated + uvw = np.append(uvw, uvw_compl, axis=0) + vis = np.append(vis, vis_compl) + + freq = freq_data[1] + freq = (freq_data[0]["IF FREQ"] + freq).reshape(-1, 1)[0] + wgt = np.ones((vis.shape[0], 1)) + mask = None + + wgt[vis == 0] = 0 + if mask is None: + mask = np.ones(wgt.shape, dtype=np.uint8) + mask[wgt == 0] = False + + DEG2RAD = np.pi / 180 + nthreads = 4 + epsilon = 1e-4 + do_wgridding = False + verbosity = 1 + + do_sycl = False # True + do_cng = False # True + + ntries = 1 + + fov_deg = 0.02 # 1e-5 # 3.3477833333331884e-5 + + npixdirty = 64 # get_npixdirty(uvw, freq, fov_deg, mask) + pixsize = fov_deg / npixdirty * DEG2RAD + + mintime = 1e300 + + grid = ms2dirty_python_fast( + uvw, freq, vis, npixdirty, npixdirty, pixsize, pixsize, epsilon, False + ) + grid = np.rot90(np.fft.fftshift(grid)) + # assert grid.shape[0] == 256 + return grid -def grid_data(uv_data, conf): - freq = conf["base_freq"] +def grid_data(uv_data, freq_data, conf): cmplx = uv_data["DATA"] - real = np.squeeze(cmplx[..., 0, 0]) - imag = np.squeeze(cmplx[..., 0, 1]) + real = np.squeeze(cmplx[..., 0, 0]).ravel() + imag = np.squeeze(cmplx[..., 0, 1]).ravel() # weight = np.squeeze(cmplx[..., 0, 2]) + freq = freq_data[1] + IF_bands = (freq_data[0]["IF FREQ"] + freq).reshape(-1, 1) + + u = np.repeat([uv_data["UU--"]], np.squeeze(cmplx[..., 0, 0]).shape[1], axis=0) + v = np.repeat([uv_data["VV--"]], np.squeeze(cmplx[..., 0, 0]).shape[1], axis=0) + u = (u * IF_bands).T.ravel() + v = (v * IF_bands).T.ravel() + samps = np.array( [ - np.append(uv_data["UU--"] * freq, -uv_data["UU--"] * freq), - np.append(uv_data["VV--"] * freq, -uv_data["VV--"] * freq), + np.append(u, -u), + np.append(v, -v), np.append(real, real), np.append(imag, -imag), ] @@ -239,7 +242,9 @@ def grid_data(uv_data, conf): mask, *_ = np.histogram2d(samps[0], samps[1], bins=[bins, bins], normed=False) mask[mask == 0] = 1 - mask = np.rot90(mask) + # flip or rot90? Stefan used flip -> write test!!! + # mask = np.flip(mask, [1]) + mask = np.rot90(mask, 1) mask_real, x_edges, y_edges = np.histogram2d( samps[0], samps[1], bins=[bins, bins], weights=samps[2], normed=False @@ -248,8 +253,8 @@ def grid_data(uv_data, conf): samps[0], samps[1], bins=[bins, bins], weights=samps[3], normed=False ) - mask_real = np.rot90(mask_real) - mask_imag = np.rot90(mask_imag) + mask_real = np.rot90(mask_real, 1) + mask_imag = np.rot90(mask_imag, 1) mask_real /= mask mask_imag /= mask @@ -260,16 +265,24 @@ def grid_data(uv_data, conf): return gridded_vis -def convert_amp_phase(data, sky_sim=False): +def convert_amp_phase(data, sky_sim=False, rescale=False): if sky_sim: amp = np.abs(data) phase = np.angle(data) - amp = (np.log10(amp + 1e-10) / 10) + 1 + # get rescale clear + # amp = (np.log10(amp + 1e-10) / 10) + 1 + if rescale: + amp = np.array([cv2.resize(a, (128, 128)) for a in amp]) + phase = np.array([cv2.resize(p, (128, 128)) for p in phase]) data = np.stack((amp, phase), axis=1) else: - data[:, 0] = np.sqrt(data[:, 0] ** 2 + data[:, 1] ** 2) / 1e3 - data[:, 1] = np.angle(data[:, 0] + 1j * data[:, 1]) - data[:, 0] = (np.log10(data[:, 0] + 1e-10) / 10) + 1 + amp = np.abs(data) + phase = np.angle(data) + # amp = (np.log10(amp + 1e-10) / 10) + 1 + if rescale: + amp = np.array([cv2.resize(a, (128, 128)) for a in amp]) + phase = np.array([cv2.resize(p, (128, 128)) for p in phase]) + data = np.stack((amp, phase), axis=1) return data From 7c49a8210c0a0dff5dbf2eea1628e3206c953c53 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Mon, 8 Aug 2022 16:55:34 +0200 Subject: [PATCH 60/63] add vla layout --- pyvisgen/layouts/vla.txt | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) create mode 100644 pyvisgen/layouts/vla.txt diff --git a/pyvisgen/layouts/vla.txt b/pyvisgen/layouts/vla.txt new file mode 100644 index 0000000..fae9dd9 --- /dev/null +++ b/pyvisgen/layouts/vla.txt @@ -0,0 +1,29 @@ +station_name X Y Z dish_dia el_low el_high SEFD altitude +W32 1640.02760601 -4329.92665085 -2416.70580433 25 15 85 110 5030 +N24 -1660.49128127 -259.39960505 2454.41532549 25 15 85 110 5030 +OUT 0. 0. 0. 25 15 85 110 5030 +W20 733.34798474 -1932.98013829 -1078.10923116 25 15 85 110 5030 +E24 765.21097819 2889.44852033 -1108.8777769 25 15 85 110 5030 +N32 -2629.09218269 -410.65158466 3885.60273242 25 15 85 110 5030 +E36 1534.56120139 5793.91275687 -2223.48335235 25 15 85 110 5030 +N16 -801.39819671 -124.9731867 1182.11742906 25 15 85 110 5030 +W8 152.75953072 -401.27175215 -223.40441888 25 15 85 110 5030 +W28 1316.44807344 -3443.3076659 -1913.52714923 25 15 85 110 5030 +W12 306.16781395 -804.58239024 -448.07786236 25 15 85 110 5030 +W36 2000.06734992 -5299.79581739 -2962.89129891 25 15 85 110 5030 +N28 -2091.44294913 -326.5939248 3089.41066984 25 15 85 110 5030 +E32 1253.24728238 4733.61503622 -1816.91333411 25 15 85 110 5030 +N36 -3217.58422497 -502.66947861 4756.12344895 25 15 85 110 5030 +E28 998.66920939 3764.30678564 -1443.47495945 25 15 85 110 5030 +N4 -74.82206165 -11.7550797 111.62669079 25 15 85 110 5030 +E4 35.60493528 133.64967044 -51.10020632 25 15 85 110 5030 +W4 46.92150279 -122.02228857 -67.60757152 25 15 85 110 5030 +W16 499.83560566 -1317.9970782 -735.20513104 25 15 85 110 5030 +N8 -243.60167375 -38.04616288 360.03532459 25 15 85 110 5030 +E16 376.98962927 1440.98974937 -556.13018294 25 15 85 110 5030 +E8 114.43465698 438.68781542 -169.49232339 25 15 85 110 5030 +E20 560.1047806 2113.23643142 -810.69974416 25 15 85 110 5030 +E12 229.46680179 879.59799397 -339.87292654 25 15 85 110 5030 +N20 -1174.33720721 -183.30338162 1734.23854793 25 15 85 110 5030 +N12 -489.30549429 -76.3090785 721.52284041 25 15 85 110 5030 +W24 1005.42662888 -2642.99130871 -1472.20242709 25 15 85 110 5030 \ No newline at end of file From 377558dcb4a2e2d05e7bc5f28fc60b7d9ccf1939 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Mon, 8 Aug 2022 16:56:19 +0200 Subject: [PATCH 61/63] delete division by 10 --- pyvisgen/simulation/data_set.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyvisgen/simulation/data_set.py b/pyvisgen/simulation/data_set.py index eeebe5e..bf5adde 100644 --- a/pyvisgen/simulation/data_set.py +++ b/pyvisgen/simulation/data_set.py @@ -93,7 +93,7 @@ def draw_sampling_opts(conf): [datetime.strptime(time, "%d-%m-%Y %H:%M:%S") for time in start_times] ) scan_duration = ( - np.random.randint(conf["scan_duration"][0] / 10, conf["scan_duration"][1] / 10) + np.random.randint(conf["scan_duration"][0], conf["scan_duration"][1]) * conf["corr_int_time"] ) scans = np.random.randint(conf["scans"][0], conf["scans"][1]) From 955e07089b2f6a3ff6916e0bc6076b15b1d6aaa3 Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Mon, 8 Aug 2022 17:06:41 +0200 Subject: [PATCH 62/63] use uncorrupted --- pyvisgen/simulation/visibility.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/pyvisgen/simulation/visibility.py b/pyvisgen/simulation/visibility.py index 987a92b..e624d82 100644 --- a/pyvisgen/simulation/visibility.py +++ b/pyvisgen/simulation/visibility.py @@ -122,7 +122,6 @@ def vis_loop(rc, SI, num_threads=48): [], ) vis_num = np.zeros(1) - # i in total number of scans for i in range(rc["scans"]): end_idx = int((rc["scan_duration"] / rc["corr_int_time"]) + 1) t = time[i * end_idx : (i + 1) * end_idx] @@ -162,12 +161,10 @@ def vis_loop(rc, SI, num_threads=48): for wave in waves ] ) - if int_values.shape[1] == 1: + if int_values.dtype != np.complex128: continue int_values = np.swapaxes(int_values, 0, 1) vis_num = np.arange(int_values.shape[0]) + 1 + vis_num.max() - print(int_values.shape) - print(int_values[:, :, 0].shape) vis = Visibilities( torch.tensor(int_values[:, :, 0]), @@ -185,19 +182,21 @@ def vis_loop(rc, SI, num_threads=48): ) visibilities.add(vis) - # if visibilities.get_values().shape[1] < 10000: - # return 0 + # workaround to guarantee min number of visibilities + # when num vis is below N sampling is redone + # if visibilities.get_values().shape[1] < 3500: + # return 0 return visibilities def calc_vis(lm, baselines, wave, t, src_crd, array_layout, SI, rd, vis_num): - X1 = scan.direction_independent( - lm, baselines, wave, t, src_crd, array_layout, SI, rd + X1 = scan.uncorrupted( + lm, baselines, wave, t, src_crd, array_layout, SI#, rd ) if X1.shape[0] == 1: return -1 - X2 = scan.direction_independent( - lm, baselines, wave, t, src_crd, array_layout, SI, rd + X2 = scan.uncorrupted( + lm, baselines, wave, t, src_crd, array_layout, SI#, rd ) int_values = scan.integrate(X1, X2).numpy() From 4f57e0d48c93be5f3e0349a3a8c104e59d0c434a Mon Sep 17 00:00:00 2001 From: Kevin Schmidt Date: Mon, 8 Aug 2022 17:07:39 +0200 Subject: [PATCH 63/63] start description of sim chain --- examples/simulation_chain.ipynb | 1217 +++++++++++++++++++++++++++++++ 1 file changed, 1217 insertions(+) create mode 100644 examples/simulation_chain.ipynb diff --git a/examples/simulation_chain.ipynb b/examples/simulation_chain.ipynb new file mode 100644 index 0000000..0dcc085 --- /dev/null +++ b/examples/simulation_chain.ipynb @@ -0,0 +1,1217 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "organizational-boutique", + "metadata": {}, + "outputs": [], + "source": [ + "from radiosim.data import radiosim_data\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "standard-crack", + "metadata": {}, + "outputs": [], + "source": [ + "data = radiosim_data(\"../../test_radiosim/build/test_data/\")" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "processed-environment", + "metadata": {}, + "outputs": [], + "source": [ + "skysim = data[np.arange(10)][:][0]" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "corporate-stephen", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAASoAAAD7CAYAAADdL9kRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAiD0lEQVR4nO2df4xc13Xfv9+Z2V8kl1z+Fi2xppTQsh1BllRFkaPWpSUrURTB8h91YBcKWNcFgcJx5caGLdVojKAIoDaFawMpWhCxIrZWkyiyFBmGYZlgrboBElrUL0syKeuHaYoipSVXJJf7e2be6R/zdu+9b3dm3sy+nXmz/H6Ih71v7nv3nR3Onrnn3HPPoZlBCCHyTKHbAgghRDOkqIQQuUeKSgiRe6SohBC5R4pKCJF7pKiEELlnWYqK5B0kXyH5Gsn7shJKCCF82G4cFckigJ8DuB3ASQBPA/i0mf0sO/GEEAIoLePemwC8ZmZvAADJvwJwN4C6ioqkokuFWGHMjMu5/7fvuNbGzk6kuvaZZ37xpJndsZznpWE5iupyAG965ycB/Ebz24rLeKQQojHVZY8wdnYCh4/8x1TXlnjPlmU/MAXL8VEtpbUXzZhI7iN5hOSRZTxLCNEhDIYoqqY6mkHyQZKjJF9KvP752L/9Msn/3Gyc5cyoTgLY6Z1fAeBU8iIz2w9gfyycTD8h8o4Zomg2q9EeAvBnAP7n/AskP4qam+haM5slua3ZIMtRVE8D2E3ySgBvAfgUgH+xjPGEEDnAYIisks1YZj8muSvx8r8B8ICZzcbXjDYbp23Tz8wqAP4AwJMAjgJ4xMxebnc8IUReMJhVUh1t8j4A/5TkYZL/l+SvN7thOTMqmNn3AXx/OWMIIfKGtaKEtiT8z/tjd08jSgA2ArgZwK8DeITkVdYgVmpZikoIsQoxg0WpFdVZM7uxxSecBPBYrJh+QjICsAXAmXo3aAuNEGIxVkl3tMffArgVAEi+D0A/gLONbtCMSgiRIIJVpzMZieRfAtiDmol4EsDXADwI4ME4ZGEOwN5GZh8gRSWESGDWko+q2VifrtN1TyvjSFEJIRIYkN5H1RGkqIQQISZFJYToBTIy/bJCikoIEUCLwMpMt8UIkKISQiSQ6SeEyD0GyvQTQuQaA5AihUsnkaISQiQwUKafECLfmGZUQoicYwZWMkuclwlSVEKIENOMSgjRA1CKSgiRbzSjEkLkHJppRiWEyDlmYGWu21IESFEJIRajGZUQIt8YGEXdFiJAikoIEZLDLTQq7iCESBCv+qU5mlCvpHvc9yWSRnJLs3GkqIQQi6BFqY4UPATgjkXjkzsB3A7gRJpBpKiEECFmQKWc7mg6lP0YwLtLdP1XAF9GzdBsSlNFtdTUjeQmkgdJvhr/3JjmYUKIHsAMiKJ0RxuQ/DiAt8zshbT3pJlRPYTFU7f7ABwys90ADsXnQohVAqNqqgNxSXfv2NdwXHINgK8C+KNW5Gm66mdmPya5K/Hy3agVFQSAAwCeAvCVVh4shMgr1spsqdWS7r8C4EoAL5AEgCsAPEvyJjN7u95N7YYnbDez0wBgZqdJbmtzHCFE3jC0bdY1HdrsRQAL+oLkcQA3mlnDku4r7kwnuW9+WrjSzxJCZEF2Pqq4pPvfA7ia5EmSn21HonZnVO+Q3BHPpnYAGK13oZntB7AfAEim8vALIboHzcAUK3ppaFDSfb5/V5px2p1RfRfA3ri9F8ATbY4jhMgjK7jq1w5NZ1Tx1G0Pat79kwC+BuABAI/E07gTAD65kkIKITrICvqo2iXNql+9qdttGcsihMgFLa36dQRtShZChBiAKF/uZCkqIUQCAyqq6yeEyDOaUQkheoJ0mRE6hhSVECKBaUYlhMg5Mv2EED2BFJUQIs+YAVaRohJC5BkDkC9fuhSVEGIJpKiEELknX5afFJUQIoEBFrHbUgRIUQkhFiPTT4g8kNWMIWc2UhYYYJV8VdLLlzRCiBxAIEp5NBtp6XJ7f0ryGMmfknyc5EizcaSohBCLMaY7mvMQFpfbOwjgGjO7FsDPAdzfbBApKiFESOxMT3M0HWqJSslm9kMzm88j8w+olcxqiHxUYpXR6I+n/vcyU/qsrKFPqpEHusd8WVHH5jD/CsBfN7tIikoIEWKEVVMrqi2JUnj748pTTSH5VQAVAA83u1aKSgixmPQzqlYrJQMASO4FcBeA28ys6XRTikr0IOnNu8CkY8l7vdhg+MQfqZdELvlkQ9U7qXivJ//2escstBUO+CR5B4CvAPhnZjaV5h4504UQCTINT1iqUvKfARgGcJDk8yT/R7NxNKMSQizC0oUepBhnyXJ732p1HCkqIUSIoZOrfqmQohI5pt63esIPRedvIvoSfaUl2wWG1xUK6f4UoigsIxVZeaFtno8qstnwxtT+qzz4q4go/apfR2gqDcmdJH9E8ijJl0neG7++ieRBkq/GPzeuvLhCiBVnfkaV5ugQaZ5UAfBFM/sAgJsBfI7kBwHcB+CQme0GcCg+F0KsArKKTM+KpvNdMzsN4HTcvkjyKIDLAdwNYE982QEAT6G25ChEmyQ/+AWvZ+kwAwAocMC1C/1BX6kwtNDuC9prguuKninIxPe3eaZZJWHSlaMprz3trquGMkbRjH9XOH5g7XXfDDRk50zPipZ8VCR3AbgewGEA22MlBjM7TXJb9uIJITqOsXed6STXAfgOgC+Y2TiZTuOS3AdgX3viCSG6QU9m+CTZh5qSetjMHotffofkjng2tQPA6FL3xvt+9sfj5GFJQwjRCCOs2iByvws0VVSsTZ2+BeComX3d6/ougL0AHoh/PrEiEopVRv1v6kXbWrytLPT8UMXCYHCZ74caKK4P+oYKG1wbrm+trQuu64ucj6qQ8FFFnt9oujAd9E2ULiy0p6Jz7jqeC66b86IaokW7adwLefkm78UZ1S0Afh/AiySfj1/796gpqEfikPgTAD65IhIKITpKTzrTzezvUP9r8LZsxRFCdB1VoRGXJv6HPn1UuR9qUCquXWgnzbs1BRdrPIxNQd+GaGShvZ7OZFxXDD/6gwUnY3KdaC5yBtlEJTQZL0TDC+0xzwRlKRHi4GVgKCcyKVQjPwODb/wlbcROGYaEWY+u+gkhLh1aSJzXEaSohBAhMv3EpYNv0vlm1UBwlW/eFT3TCQD6PXNvqOjMu/XYGly3MXLm3uZExPnIgJNjpN/JMdwXmlH9hfpm1UzV3Xd+LjRPh2a9lUpvZa9cTESwF91qYTWaC/oiejcGm5e7g8n0E0L0AppRCSHyjfVgeIIQ4tJDikqsUsKocnrZCHzfU9IPNVB0y/uDxQ1B3xovxdmI54faiuHgus0D7mO8dZCJPrf0v77PtdcUq8F1pYILBagm/kgnyu53GUgk2Kt6vpypqgt/GLdQxknv955NjEFvA7D5hSUWFWepohMYiCijLTQkH0St2syomV0Tv7YJtVp+uwAcB/B7Znau3hiAijsIIZJkWCkZS5d0bzmXnRSVEGIRZkx1NB9ncUl31HLZHYjbBwB8otk4Mv1ECzRIbJfIQV7yTLqhkjPb1hXDtGUbsMX1RWHUtx9JPlJy4yfNu8uGnEm0dSBMSreh34UCrO1zfX2F0IyKvD+62Ur9P4upRCBkf8GZSH1eSHvJwjEaJebLIyvso2o5l50UlRAixFpKM9x2SfdWkKISQgQYgChK7Uxvp6R7qlx2PvmfgwohOk5kTHW0yXwuOyBlLjvNqMQSLJ3tgEkflVdkoVhYG3St7du+0N5e+JWF9hUWuiO29jvfzXBfOP6Q96U+3OfCBzYn/FBbB9x2lU2DU0HfOq+v6IcgJHKCT5fdVp6kj8r315QTJlHF/LY7qTCs/1f1CjrYoqwIOaM1068hcUn3PaiZiCcBfA1t5LKTohJCBGSZOK9OSXegxVx2UlRCiEUoMl3kBC+jwaJc5X7pc5ftgAzNJb9vqC9MWOebex8ovGeh/f4N4R/AzjXONBtKRIsXvFogfl2QtaVECMKgq5m3bmAm6OsrOhOs4kVbzyXMu+mKM0HPz4UZHsZmXd+5ufA9uDDn5Lro1fybLoQmaKXq+pJl4QNT0HyzsHsmohSVECLXmBHVXqtCI4S49NCMSmRAow9RuogTP1d5csWuv+SVmPIS1pUSSe98hrE5OH+vuVW/Xxtxr9+0ZSy8bosLoSklTL+JKZcEb3LWRalXE0ndSoX6m3VnPZNuyjPpzk+Hm6PPzLjzt2fCsvCjMwWvHZpjo2Vn0o0V3E6RCTsbXDcXTSy0o0RZeDN/RTAfBbOkqIQQ+cawnBipFUGKSggRUEtFLEUlhMg5UlSXPGk/APV9Tb5/aVGdPD/UgPXHKHl+qeH+HUHfTrx/of0euhp6/YVwvMlqfd/QzjXON7R72Pln3n/Fm8F123adXGgn/zjOnXJR7JVR176Q8C/NVJz/6uJsWO59zlu9Gvd8VGdmk36oktcO5Tgz68IJxqIw7OCc55cax5mF9nTlfHBduTq50I4SxR38gg5hSEL3/FXJyP1u01QakoMkf0LyBZIvk/zj+PVNJA+SfDX+ubHZWEKI/GOWXT6qrEijNmcB3GpmHwJwHYA7SN6MNrL0CSF6gXQbkjvpcG9q+pmZAZifu/fFh6GWpW9P/PoBAE8B+ErmEvYE6cMFFm3sDTr9a+snpSvQL3UemkHFOiEEhYQZ6Icd7I5+Lei7eaMb89qNzmzrL4RL87+ccObjW9PhR2m45MyWDd7G4LWeGQgA/RvceXUmlN3fGDvhmXSnp8IEe5NelHklYbL4NfkulJ0ZODYb/j+cnXW/29nKdND3buG8G6MQhh1MRS7V92x1fKHtm3oAEEUuYt4QRtbnJSTBJ28+qlSGKMkiyedRyxtz0MwOI5GlD0DTLH1CiN4gb6ZfKme6mVUBXEdyBMDjJK9J+wCS+wDsa088IUSnMcufM72lVT8zO0/yKdSqSqTK0henJd0PAPR3lgohckpn/U9paKqoSG4FUI6V1BCAjwH4T3BZ+h5Ayix9vUc639PihHJ+X1+DvoT/yuvzMxOUiuGSu+9fGi6EFveQuW0nBW+rScTQvzQSubCDf7xhTdB395W/XGhf/5HDbry+cNf/64evXWj/v1/8atD37pz7vafL7mN24VxYu6/vlPNfzU2Fv+epM+53e+PCyEL7+GToyxovu/d/JhExMVN1342TFfcenKuGIQLneMGNVwxLzE1GbtvPbDX0sfm+qGrkfFtmiQwJ3raZxT6pfIQk+BJY6jCazpBmRrUDwAHWgncKAB4xs++R/Hu0mKVPCNEb5M2ZnmbV76cArl/i9TG0mKVPCNEbZGn6kfx3AP41apO1FwF8xsxmGt8Vosj0dkMLvORyjcIHioXQTCkWvL5F97nzPq8E+DpuCa7bXnWZCXYwzHwwMrC0E9Q3gQBgbcn9LjdsCs2Z637z6YW2fe72hXalFP4uv1r8zkL79dEwuv3taWeeHp90Mva9tTO4bvM5l3BvqhxGi78+7szEV8bds9+cDO27C1W33D+bWPqfpTufposqnyxeCK4LwgzK9c27ZOYDP8o8CDuw0NS2huXYu2/uhWS3okfycgD/FsAHzWya5CMAPoVaBeXUSFEJIQJWYNWvBGCIZBnAGgCnWh0gX2uQQohcEIGpjmaY2VsA/gtqfuzTAC6Y2Q9blecSmlHVe1PTmXdAaOIVC26FqlQIo8P7iy5yeqAQRlEP0Ouz8L4Bc2MORa5vSyKx3a5hJ9fu4dDU2T7kRXp70/ex2dBs88s+XT48HvQN/iMXfV14z2+jHtHI/1poW8J6OTPr3tfJST/xXLgldK2XpM+PIgeA016A+IkZd3KmECbf8824OQs3DVe9pHTlquurJEy4StW5TPzVOyBcwUtGlfsmXrial8x3njfzrj4tVqFpWCk53gN8N4ArAZwH8Dck7zGzb7ci0yWkqIQQ6WgpjqpZpeSPAfiFmZ0BAJKPAfhNAFJUQojlkWF4wgkAN5NcA2AatUiBI41vWYwUlRAiwJBdoS4zO0zyUQDPAqgAeA7xTpVWWMWKKvmNUCeSPOGH8iPCC4Vwudz3Rfl+qDWF0O8yDLfkviEaCfrWwY2/phg+e7DoZFzjhQ9cNhT6N65e75bLr9n6dtC3fZtL3lbxIsLffOey4Lq3LjrfUPLbc+aXLiK8/5lvuI4oXGIff/HKhfbxieGg7/ik8+ucjZy8A9NhSEaf9/9STvx5+FkLzhfcDq3JapjBoOz5lKJERHjkF07w+pLX+QUWkn+mta2uS/eF9I4fqiEZr/qZ2ddQK+XeNqtYUQkh2qUXt9AIIS4hrBc3JfcW/pub3PBb9K5y5kehEG6E9RPRDRTXB32DBXe+AVsX2huj0PTb7JmIm4fCt3jEsybXlkJTYbDozIohr719MFxKv9qrhXfV+14P+obf60zByoSTo1wOTa6xKRfycCphtv30H25w47993I2XGOMnx9xiz3PvhmbyGzjhnsWTqIcfnV+1cOl/ztsAPFdn8y+QjA5vFAFer3R6s9CC8MpLgShnv+YqU1RCiCyQ6SeEyDW1LTRSVEKInJNme0wn6XFF1SAEIah9B9ALO/C3v/Qn/FBDpZGF9vpEGvgRzxe12dsKs3UwfBu3DXqhBYPhMvjGAedPGS6FPpk+r3hC0WtvGAgzYmza4LaMDKwPiwgUhpw/i9POb5QMQTjv1bh7fSLcyvPqRZdIb+TU5QvtmcSS9c/H3Xv87EwYJnE6Oubu82rcmTXy/4R9fghBED6QDC1I7TeS7ykNLW6h6Qg9rqiEENmjVT8hRA+Qt/nlqlJUfsR5Mld5PXNvXWlrcN0muARwW6NNQd+WPmdKbRt0ZtBlQ+GS+GWDzrzbNhTu5h8ZdEvrQ/1h3m7/W8yPDC4m6umVvYjzybMjQV9lxsk49a6LPn/zbGjGnph078exC+HH8t2KkyvyPrLJpHRnii5MYix6I+ibmnN9kXnhBA1NvxCFDHQHQ7YZPrNgVSkqIUQ2VKWohBB5xkwzqgyoH30elKJKbDb2NxT7K3sj2B5ct90z9y4bCJPNvcdbHNs+5Mygy4bCVbktg87c27Q2XJVb45l+SZNuds6ZbZNeCXO/DQBT3orduxfDVcuCVzpxbNKtTL58fiS47ti4ex9frb4T9J0pvLnQnovc71JORISXy+53q1QvBn1heahG0eKNkEnXLfL2zvegohJCrDSaUQkhck2W+aiyQopKCJEgu3JZWbHKFFUjH5Xz6wzRLduvr4Y+nk1BCEL4n7Wjjl9qy1Doh1o/5Hw5g4kQBHo+pHIlDKHwfVHnvOwGYzNh5PjFivvd5qqhn86PHh+bdZHjJybD78jXKq5Awmm8FvRNzDqfVcMy5fCT0iX9UO36pUS3MeRv1S91Gj+SRZLPkfxefL6J5EGSr8Y/NzYbQwjRG0SW7ugUreQbvRfAUe/8PgCHzGw3gEPxuRBiFWApjzSQHCH5KMljJI+S/HCr8qQy/UheAeB3AfwJgD+MX74bwJ64fQDAUwC+0qoAy4GJTcmEM3UKiXLpRW9Tcp95ZiDChG9rS053D/dFiT5nzgx5G4pLxWSObSfXzFw4frnqzLjJRK29c56J97bXHp1JJL2bdeOPz4Ufl4mqk+W8V8fuTCHMM37WS2w3MReGJ5Sr591JUNOul8qSi3ZZgTiqbwL4gZn9c5L9qFVLbom0PqpvAPgyAD8V5HYzOw0AZnaa5LalbhRC9B5ZrfqRXA/gIwD+JQCY2RyAuUb3LEVT04/kXQBGzeyZVgeP799H8kiimqoQIqfMO9PTHIgrJXvHvsRwVwE4A+AvYh/3n5Ncm3xmM9LMqG4B8HGSdwIYBLCe5LcBvENyRzyb2gFgdKmb4/LO+wGA/pKXECK3WPq/1GaVkksAbgDw+bjG3zdR82f/h1bkaaqozOx+APcDAMk9AL5kZveQ/FMAewE8EP98opUH5xV/JSO5ROufz1WdP2wq4Yea9oogVBLJ5ibL7tp3E/eNzbr73p52943OhBPxsxUXGnGB4daVi4XxhfaE55ea9pLXAcCsd16NwvCK+jXu9D1zqZBhhs+TAE6a2eH4/FG0sfC2nCqDDwC4neSrAG6Pz4UQPU4tzUs24Qlm9jaAN0leHb90G4CftSpTSwGfZvYUaqt7MLOx+KFCiFVGC6ZfGj4P4OF4xe8NAJ9pdYCejkxfnCvbLZ9HiRpxVW83f5muPZ1YgJisOHNsvBxOOAe9Euxl8+r/zYVhBv7S7lQicny87EzGsdmwb2zW/T5n55yMZ3khuO5c0ZVtn4rOBX0zZXdt2auF59e+q8noR5wnww5k7l3aMNPiDmb2PIBGfqym9LSiEkJkjxlQzdn3kxSVEGIRqpS8ojiTJbmBthJ55p458+hCIQzpGPA2ChemQ5NuxlvpW1ty7WJillz2LKfphFV1wYskP18JzbFz5ky1d4tu0/DFKIz8mCm7lb256njQ55t4kWfuLi4xVQ3OhPDJ2ydilSkqIcRyUXEHIURPkPGq37KRohJCBNS20HRbipAeVFQNar1501VD0kflluP9KO3zpbD0e+QVXJgth3X9zpW9rAuF8L7gWd7X0VTCNzRBJ8c4zyf6nF9quuLCDuaqYeS4n8yumii4EGY7aFQXL2efRJErlIpYCJFvOpwULw1SVEKIgFaS4nWKVaWoQlMnEZkeuY28flBAMhlcpeiW9KcK4dL/IFydvL5Ewr3gWXTm12whNM380Ii56kTQNxf5dfLcfZGFYQxBzTxFlYsVQDMqIUTu0aqfECLXGICKFJUQIu/kTE/1uqJKvp3+Fprktc6vU63W32pT9bagzBTCrAUlr0AEWT88wfcbVfxtLInxq1HY5/uiguR1i7a/KOxArBzz+ajyRI8rKiFE5ph8VEKIHkABnytKfZMo/IZwfdUoXN73sw9UEmXh59j627WoDHqjOnmKKhc5QKafEKInyHKvH2sO3SMA3jKzu9oZQ4pKCLGIjH1U9wI4CmB9uwOsYkXVYEUwuCx5nW9+hdHt4Y2NCvg0sPDNl6O+jA0eLMSKYsi0UvIVAH4XwJ8A+MN2x1nFikoI0S5RdlOqbwD4MoDh5QyynLp+QohViqU80KCkO8m7AIya2TPLlUczKiFEgJmhmn5G1aik+y0APk7yTgCDANaT/LaZ3dOqTJeQoqr3xkf1r1qUmcCfgCb76o/ZukxCdI+swhPM7H4A9wMAyT0AvtSOkgJSKiqSxwFcRO2vs2JmN5LcBOCvAewCcBzA75nZuXpjCCF6h7wFfLbio/qomV3nTfPuA3DIzHYDOBSfCyFWAWaW6mhhvKfajaECludMvxvAgbh9AMAnljFWF0nrNjTUJpRpjhZckULkjPnwhDRHp0irqAzAD0k+43n1t5vZaQCIf25bCQGFEJ0n6xnVcknrTL/FzE6R3AbgIMljaR8QK7Z9TS8UQuSCWuK8fM34U82ozOxU/HMUwOMAbgLwDskdABD/HK1z734zu7HBEqYQImdYyn+doqmiIrmW5PB8G8BvAXgJwHcB7I0v2wvgiZUSUgjRWfLmo0pj+m0H8DjJ+ev/t5n9gOTTAB4h+VkAJwB8cuXEFEJ0CoMhytliT1NFZWZvAPjQEq+PAbhtJYQSQnQRy3SvXyZcQpHpQoi0dNL/lAYpKiFEgAGo5Cw2XYpKCJGgsyt6aZCiEkIE1CLTpaiEEHmGQESZfkKInKMZlRAi1xgM1Yb51jqPFJUQYhEy/YQQuaYWmS5FJYTIOXlTVKpCI4RIML/br/nRDJI7Sf6I5FGSL5O8tx2JNKMSQgQYMvVRVQB80cyejbOwPEPyoJn9rJVBpKiEEAkM1WSV8HZHqmX/nc8EfJHkUQCXA5CiEkK0z0o500nuAnA9gMOt3itFJYRYRAuKagvJI975fjPbn7yI5DoA3wHwBTMbb1UeKSohRIJayGdKGlVKBgCQ7ENNST1sZo+1I5EUlRAioLYpORvTj7XUwN8CcNTMvt7uOApPEEIsIqvwBAC3APh9ALeSfD4+7mxVHs2ohBABlu2q398B4HLHkaISQiQwRKZNyUKInJPSrOsYUlRCiAQtrfp1BCkqIUSAAYhMMyohRJ4xQ2TZONOzQopKCBGQx3xUqeKoSI6QfJTksThdw4dJbiJ5kOSr8c+NKy2sEKIzmEWpjk6RNuDzmwB+YGbvR628+1EA9wE4ZGa7ARyKz4UQPU/NmZ7m6BRNFRXJ9QA+gloYPMxszszOA7gbwIH4sgMAPrEyIgohOk0vzqiuAnAGwF+QfI7kn5NcC2B7nGtmPufMthWUUwjRMbLL8JkVaRRVCcANAP67mV0PYBItmHkk95E8kkgFIYTIKQZDFJVTHZ0ijaI6CeCkmc0nu3oUNcX1DskdABD/HF3qZjPbb2Y3NksFIYTIDz03ozKztwG8SfLq+KXbUEsj+l0Ae+PX9gJ4YkUkFEJ0FsufjyptHNXnATxMsh/AGwA+g5qSe4TkZwGcAPDJlRFRCNFZrDf3+pnZ8wCWMt1uy1QaIUTXMQCm7AlCiHxTy/GZJ6SohBAJDJFVui1EgBSVEGIJ8jWjUs50IcRiLEp3pIDkHSRfIfkayba22klRCSESZBeZTrII4L8B+B0AHwTwaZIfbFUiKSohxBJEKY+m3ATgNTN7w8zmAPwVavuEW0I+KiFEAssymPNyAG965ycB/Earg3RaUZ0Fqr8EsKXW7jqSI0RyhORBjlZleG8Gz3wSqGxJee1gk5LuS5XKslYF6qiiMrOtAEDySB72/kkOyZF3Obohg5ndkeFwJwHs9M6vAHCq1UHkoxJCrCRPA9hN8sp4C96nUNsn3BLyUQkhVgwzq5D8AwBPAigCeNDMXm51nG4pqv3NL+kIkiNEcoTkQY48yLAszOz7AL6/nDFo1rJfSwghOop8VEKI3NNRRZVFKH2bz32Q5CjJl7zXOl7ui+ROkj+KS469TPLebshCcpDkT0i+EMvxx92Qw5OnGOfj/1635CB5nOSLJJ+fX27vkhwqTbcEHVNUWYXSt8lDAJJLrt0o91UB8EUz+wCAmwF8Ln4POi3LLIBbzexDAK4DcAfJm7sgxzz3olaCbZ5uyfFRM7vOCwfohhwqTbcUZtaRA8CHATzpnd8P4P4OPn8XgJe881cA7IjbOwC80ilZPBmeAHB7N2UBsAbAs6hFC3dcDtTiag4BuBXA97r1fwPgOIAtidc6KgeA9QB+gdh33C058nh00vRbKpT+8g4+P0lXy32R3AXgegCHuyFLbG49j1pRjoNWK97RjffkGwC+jHDjWDfkMAA/JPkMyX1dkkOl6erQSUWVSSj9aoDkOgDfAfAFMxvvhgxmVjWz61Cb0dxE8ppOy0DyLgCjZvZMp5+9BLeY2Q2ouSY+R/IjXZBhWaXpVjOdVFSZhNJnSKpyX1lDsg81JfWwmT3WTVkAwGpVr59CzYfXaTluAfBxksdR21V/K8lvd0EOmNmp+OcogMdR2/XfaTmWVZpuNdNJRZVJKH2GdLzcF0kC+BaAo2b29W7JQnIryZG4PQTgYwCOdVoOM7vfzK4ws12ofR7+j5nd02k5SK4lOTzfBvBbAF7qtBym0nT16aRDDMCdAH4O4HUAX+3gc/8SwGkAZdS+tT4LYDNqTtxX45+bOiDHP0HN3P0pgOfj485OywLgWgDPxXK8BOCP4tc7/p54Mu2Bc6Z3+v24CsAL8fHy/GezS5+R6wAcif9v/hbAxm7+v+TlUGS6ECL3KDJdCJF7pKiEELlHikoIkXukqIQQuUeKSgiRe6SohBC5R4pKCJF7pKiEELnn/wOS4GosRdyuGwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# just for inspection\n", + "plt.imshow(skysim[0].T, cmap=\"inferno\")\n", + "plt.colorbar()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "abroad-afghanistan", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "labeled-bouquet", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "reverse-award", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 54, + "id": "6dfb3f84", + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "from datetime import datetime" + ] + }, + { + "cell_type": "code", + "execution_count": 78, + "id": "50ffd1f6", + "metadata": {}, + "outputs": [], + "source": [ + "date_str = pd.date_range(\"0:00\", \"23:59\", freq=\"1min\").strftime('%H:%M:%S')" + ] + }, + { + "cell_type": "code", + "execution_count": 76, + "id": "fb480fad", + "metadata": {}, + "outputs": [], + "source": [ + "times = [datetime.time(datetime.strptime(date, \"%H:%M:%S\")) for date in date_str]" + ] + }, + { + "cell_type": "code", + "execution_count": 77, + "id": "6a034205", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[datetime.time(0, 0),\n", + " datetime.time(0, 1),\n", + " datetime.time(0, 2),\n", + " datetime.time(0, 3),\n", + " datetime.time(0, 4),\n", + " datetime.time(0, 5),\n", + " datetime.time(0, 6),\n", + " datetime.time(0, 7),\n", + " datetime.time(0, 8),\n", + " datetime.time(0, 9),\n", + " datetime.time(0, 10),\n", + " datetime.time(0, 11),\n", + " datetime.time(0, 12),\n", + " datetime.time(0, 13),\n", + " datetime.time(0, 14),\n", + " datetime.time(0, 15),\n", + " datetime.time(0, 16),\n", + " datetime.time(0, 17),\n", + " datetime.time(0, 18),\n", + " datetime.time(0, 19),\n", + " datetime.time(0, 20),\n", + " datetime.time(0, 21),\n", + " datetime.time(0, 22),\n", + " datetime.time(0, 23),\n", + " datetime.time(0, 24),\n", + " datetime.time(0, 25),\n", + " datetime.time(0, 26),\n", + " datetime.time(0, 27),\n", + " datetime.time(0, 28),\n", + " datetime.time(0, 29),\n", + " datetime.time(0, 30),\n", + " datetime.time(0, 31),\n", + " datetime.time(0, 32),\n", + " datetime.time(0, 33),\n", + " datetime.time(0, 34),\n", + " datetime.time(0, 35),\n", + " datetime.time(0, 36),\n", + " datetime.time(0, 37),\n", + " datetime.time(0, 38),\n", + " datetime.time(0, 39),\n", + " datetime.time(0, 40),\n", + " datetime.time(0, 41),\n", + " datetime.time(0, 42),\n", + " datetime.time(0, 43),\n", + " datetime.time(0, 44),\n", + " datetime.time(0, 45),\n", + " datetime.time(0, 46),\n", + " datetime.time(0, 47),\n", + " datetime.time(0, 48),\n", + " datetime.time(0, 49),\n", + " datetime.time(0, 50),\n", + " datetime.time(0, 51),\n", + " datetime.time(0, 52),\n", + " datetime.time(0, 53),\n", + " datetime.time(0, 54),\n", + " datetime.time(0, 55),\n", + " datetime.time(0, 56),\n", + " datetime.time(0, 57),\n", + " datetime.time(0, 58),\n", + " datetime.time(0, 59),\n", + " datetime.time(1, 0),\n", + " datetime.time(1, 1),\n", + " datetime.time(1, 2),\n", + " datetime.time(1, 3),\n", + " datetime.time(1, 4),\n", + " datetime.time(1, 5),\n", + " datetime.time(1, 6),\n", + " datetime.time(1, 7),\n", + " datetime.time(1, 8),\n", + " datetime.time(1, 9),\n", + " datetime.time(1, 10),\n", + " datetime.time(1, 11),\n", + " datetime.time(1, 12),\n", + " datetime.time(1, 13),\n", + " datetime.time(1, 14),\n", + " datetime.time(1, 15),\n", + " datetime.time(1, 16),\n", + " datetime.time(1, 17),\n", + " datetime.time(1, 18),\n", + " datetime.time(1, 19),\n", + " datetime.time(1, 20),\n", + " datetime.time(1, 21),\n", + " datetime.time(1, 22),\n", + " datetime.time(1, 23),\n", + " datetime.time(1, 24),\n", + " datetime.time(1, 25),\n", + " datetime.time(1, 26),\n", + " datetime.time(1, 27),\n", + " datetime.time(1, 28),\n", + " datetime.time(1, 29),\n", + " datetime.time(1, 30),\n", + " datetime.time(1, 31),\n", + " datetime.time(1, 32),\n", + " datetime.time(1, 33),\n", + " datetime.time(1, 34),\n", + " datetime.time(1, 35),\n", + " datetime.time(1, 36),\n", + " datetime.time(1, 37),\n", + " datetime.time(1, 38),\n", + " datetime.time(1, 39),\n", + " datetime.time(1, 40),\n", + " datetime.time(1, 41),\n", + " datetime.time(1, 42),\n", + " datetime.time(1, 43),\n", + " datetime.time(1, 44),\n", + " datetime.time(1, 45),\n", + " datetime.time(1, 46),\n", + " datetime.time(1, 47),\n", + " datetime.time(1, 48),\n", + " datetime.time(1, 49),\n", + " datetime.time(1, 50),\n", + " datetime.time(1, 51),\n", + " datetime.time(1, 52),\n", + " datetime.time(1, 53),\n", + " datetime.time(1, 54),\n", + " datetime.time(1, 55),\n", + " datetime.time(1, 56),\n", + " datetime.time(1, 57),\n", + " datetime.time(1, 58),\n", + " datetime.time(1, 59),\n", + " datetime.time(2, 0),\n", + " datetime.time(2, 1),\n", + " datetime.time(2, 2),\n", + " datetime.time(2, 3),\n", + " datetime.time(2, 4),\n", + " datetime.time(2, 5),\n", + " datetime.time(2, 6),\n", + " datetime.time(2, 7),\n", + " datetime.time(2, 8),\n", + " datetime.time(2, 9),\n", + " datetime.time(2, 10),\n", + " datetime.time(2, 11),\n", + " datetime.time(2, 12),\n", + " datetime.time(2, 13),\n", + " datetime.time(2, 14),\n", + " datetime.time(2, 15),\n", + " datetime.time(2, 16),\n", + " datetime.time(2, 17),\n", + " datetime.time(2, 18),\n", + " datetime.time(2, 19),\n", + " datetime.time(2, 20),\n", + " datetime.time(2, 21),\n", + " datetime.time(2, 22),\n", + " datetime.time(2, 23),\n", + " datetime.time(2, 24),\n", + " datetime.time(2, 25),\n", + " datetime.time(2, 26),\n", + " datetime.time(2, 27),\n", + " datetime.time(2, 28),\n", + " datetime.time(2, 29),\n", + " datetime.time(2, 30),\n", + " datetime.time(2, 31),\n", + " datetime.time(2, 32),\n", + " datetime.time(2, 33),\n", + " datetime.time(2, 34),\n", + " datetime.time(2, 35),\n", + " datetime.time(2, 36),\n", + " datetime.time(2, 37),\n", + " datetime.time(2, 38),\n", + " datetime.time(2, 39),\n", + " datetime.time(2, 40),\n", + " datetime.time(2, 41),\n", + " datetime.time(2, 42),\n", + " datetime.time(2, 43),\n", + " datetime.time(2, 44),\n", + " datetime.time(2, 45),\n", + " datetime.time(2, 46),\n", + " datetime.time(2, 47),\n", + " datetime.time(2, 48),\n", + " datetime.time(2, 49),\n", + " datetime.time(2, 50),\n", + " datetime.time(2, 51),\n", + " datetime.time(2, 52),\n", + " datetime.time(2, 53),\n", + " datetime.time(2, 54),\n", + " datetime.time(2, 55),\n", + " datetime.time(2, 56),\n", + " datetime.time(2, 57),\n", + " datetime.time(2, 58),\n", + " datetime.time(2, 59),\n", + " datetime.time(3, 0),\n", + " datetime.time(3, 1),\n", + " datetime.time(3, 2),\n", + " datetime.time(3, 3),\n", + " datetime.time(3, 4),\n", + " datetime.time(3, 5),\n", + " datetime.time(3, 6),\n", + " datetime.time(3, 7),\n", + " datetime.time(3, 8),\n", + " datetime.time(3, 9),\n", + " datetime.time(3, 10),\n", + " datetime.time(3, 11),\n", + " datetime.time(3, 12),\n", + " datetime.time(3, 13),\n", + " datetime.time(3, 14),\n", + " datetime.time(3, 15),\n", + " datetime.time(3, 16),\n", + " datetime.time(3, 17),\n", + " datetime.time(3, 18),\n", + " datetime.time(3, 19),\n", + " datetime.time(3, 20),\n", + " datetime.time(3, 21),\n", + " datetime.time(3, 22),\n", + " datetime.time(3, 23),\n", + " datetime.time(3, 24),\n", + " datetime.time(3, 25),\n", + " datetime.time(3, 26),\n", + " datetime.time(3, 27),\n", + " datetime.time(3, 28),\n", + " datetime.time(3, 29),\n", + " datetime.time(3, 30),\n", + " datetime.time(3, 31),\n", + " datetime.time(3, 32),\n", + " datetime.time(3, 33),\n", + " datetime.time(3, 34),\n", + " datetime.time(3, 35),\n", + " datetime.time(3, 36),\n", + " datetime.time(3, 37),\n", + " datetime.time(3, 38),\n", + " datetime.time(3, 39),\n", + " datetime.time(3, 40),\n", + " datetime.time(3, 41),\n", + " datetime.time(3, 42),\n", + " datetime.time(3, 43),\n", + " datetime.time(3, 44),\n", + " datetime.time(3, 45),\n", + " datetime.time(3, 46),\n", + " datetime.time(3, 47),\n", + " datetime.time(3, 48),\n", + " datetime.time(3, 49),\n", + " datetime.time(3, 50),\n", + " datetime.time(3, 51),\n", + " datetime.time(3, 52),\n", + " datetime.time(3, 53),\n", + " datetime.time(3, 54),\n", + " datetime.time(3, 55),\n", + " datetime.time(3, 56),\n", + " datetime.time(3, 57),\n", + " datetime.time(3, 58),\n", + " datetime.time(3, 59),\n", + " datetime.time(4, 0),\n", + " datetime.time(4, 1),\n", + " datetime.time(4, 2),\n", + " datetime.time(4, 3),\n", + " datetime.time(4, 4),\n", + " datetime.time(4, 5),\n", + " datetime.time(4, 6),\n", + " datetime.time(4, 7),\n", + " datetime.time(4, 8),\n", + " datetime.time(4, 9),\n", + " datetime.time(4, 10),\n", + " datetime.time(4, 11),\n", + " datetime.time(4, 12),\n", + " datetime.time(4, 13),\n", + " datetime.time(4, 14),\n", + " datetime.time(4, 15),\n", + " datetime.time(4, 16),\n", + " datetime.time(4, 17),\n", + " datetime.time(4, 18),\n", + " datetime.time(4, 19),\n", + " datetime.time(4, 20),\n", + " datetime.time(4, 21),\n", + " datetime.time(4, 22),\n", + " datetime.time(4, 23),\n", + " datetime.time(4, 24),\n", + " datetime.time(4, 25),\n", + " datetime.time(4, 26),\n", + " datetime.time(4, 27),\n", + " datetime.time(4, 28),\n", + " datetime.time(4, 29),\n", + " datetime.time(4, 30),\n", + " datetime.time(4, 31),\n", + " datetime.time(4, 32),\n", + " datetime.time(4, 33),\n", + " datetime.time(4, 34),\n", + " datetime.time(4, 35),\n", + " datetime.time(4, 36),\n", + " datetime.time(4, 37),\n", + " datetime.time(4, 38),\n", + " datetime.time(4, 39),\n", + " datetime.time(4, 40),\n", + " datetime.time(4, 41),\n", + " datetime.time(4, 42),\n", + " datetime.time(4, 43),\n", + " datetime.time(4, 44),\n", + " datetime.time(4, 45),\n", + " datetime.time(4, 46),\n", + " datetime.time(4, 47),\n", + " datetime.time(4, 48),\n", + " datetime.time(4, 49),\n", + " datetime.time(4, 50),\n", + " datetime.time(4, 51),\n", + " datetime.time(4, 52),\n", + " datetime.time(4, 53),\n", + " datetime.time(4, 54),\n", + " datetime.time(4, 55),\n", + " datetime.time(4, 56),\n", + " datetime.time(4, 57),\n", + " datetime.time(4, 58),\n", + " datetime.time(4, 59),\n", + " datetime.time(5, 0),\n", + " datetime.time(5, 1),\n", + " datetime.time(5, 2),\n", + " datetime.time(5, 3),\n", + " datetime.time(5, 4),\n", + " datetime.time(5, 5),\n", + " datetime.time(5, 6),\n", + " datetime.time(5, 7),\n", + " datetime.time(5, 8),\n", + " datetime.time(5, 9),\n", + " datetime.time(5, 10),\n", + " datetime.time(5, 11),\n", + " datetime.time(5, 12),\n", + " datetime.time(5, 13),\n", + " datetime.time(5, 14),\n", + " datetime.time(5, 15),\n", + " datetime.time(5, 16),\n", + " datetime.time(5, 17),\n", + " datetime.time(5, 18),\n", + " datetime.time(5, 19),\n", + " datetime.time(5, 20),\n", + " datetime.time(5, 21),\n", + " datetime.time(5, 22),\n", + " datetime.time(5, 23),\n", + " datetime.time(5, 24),\n", + " datetime.time(5, 25),\n", + " datetime.time(5, 26),\n", + " datetime.time(5, 27),\n", + " datetime.time(5, 28),\n", + " datetime.time(5, 29),\n", + " datetime.time(5, 30),\n", + " datetime.time(5, 31),\n", + " datetime.time(5, 32),\n", + " datetime.time(5, 33),\n", + " datetime.time(5, 34),\n", + " datetime.time(5, 35),\n", + " datetime.time(5, 36),\n", + " datetime.time(5, 37),\n", + " datetime.time(5, 38),\n", + " datetime.time(5, 39),\n", + " datetime.time(5, 40),\n", + " datetime.time(5, 41),\n", + " datetime.time(5, 42),\n", + " datetime.time(5, 43),\n", + " datetime.time(5, 44),\n", + " datetime.time(5, 45),\n", + " datetime.time(5, 46),\n", + " datetime.time(5, 47),\n", + " datetime.time(5, 48),\n", + " datetime.time(5, 49),\n", + " datetime.time(5, 50),\n", + " datetime.time(5, 51),\n", + " datetime.time(5, 52),\n", + " datetime.time(5, 53),\n", + " datetime.time(5, 54),\n", + " datetime.time(5, 55),\n", + " datetime.time(5, 56),\n", + " datetime.time(5, 57),\n", + " datetime.time(5, 58),\n", + " datetime.time(5, 59),\n", + " datetime.time(6, 0),\n", + " datetime.time(6, 1),\n", + " datetime.time(6, 2),\n", + " datetime.time(6, 3),\n", + " datetime.time(6, 4),\n", + " datetime.time(6, 5),\n", + " datetime.time(6, 6),\n", + " datetime.time(6, 7),\n", + " datetime.time(6, 8),\n", + " datetime.time(6, 9),\n", + " datetime.time(6, 10),\n", + " datetime.time(6, 11),\n", + " datetime.time(6, 12),\n", + " datetime.time(6, 13),\n", + " datetime.time(6, 14),\n", + " datetime.time(6, 15),\n", + " datetime.time(6, 16),\n", + " datetime.time(6, 17),\n", + " datetime.time(6, 18),\n", + " datetime.time(6, 19),\n", + " datetime.time(6, 20),\n", + " datetime.time(6, 21),\n", + " datetime.time(6, 22),\n", + " datetime.time(6, 23),\n", + " datetime.time(6, 24),\n", + " datetime.time(6, 25),\n", + " datetime.time(6, 26),\n", + " datetime.time(6, 27),\n", + " datetime.time(6, 28),\n", + " datetime.time(6, 29),\n", + " datetime.time(6, 30),\n", + " datetime.time(6, 31),\n", + " datetime.time(6, 32),\n", + " datetime.time(6, 33),\n", + " datetime.time(6, 34),\n", + " datetime.time(6, 35),\n", + " datetime.time(6, 36),\n", + " datetime.time(6, 37),\n", + " datetime.time(6, 38),\n", + " datetime.time(6, 39),\n", + " datetime.time(6, 40),\n", + " datetime.time(6, 41),\n", + " datetime.time(6, 42),\n", + " datetime.time(6, 43),\n", + " datetime.time(6, 44),\n", + " datetime.time(6, 45),\n", + " datetime.time(6, 46),\n", + " datetime.time(6, 47),\n", + " datetime.time(6, 48),\n", + " datetime.time(6, 49),\n", + " datetime.time(6, 50),\n", + " datetime.time(6, 51),\n", + " datetime.time(6, 52),\n", + " datetime.time(6, 53),\n", + " datetime.time(6, 54),\n", + " datetime.time(6, 55),\n", + " datetime.time(6, 56),\n", + " datetime.time(6, 57),\n", + " datetime.time(6, 58),\n", + " datetime.time(6, 59),\n", + " datetime.time(7, 0),\n", + " datetime.time(7, 1),\n", + " datetime.time(7, 2),\n", + " datetime.time(7, 3),\n", + " datetime.time(7, 4),\n", + " datetime.time(7, 5),\n", + " datetime.time(7, 6),\n", + " datetime.time(7, 7),\n", + " datetime.time(7, 8),\n", + " datetime.time(7, 9),\n", + " datetime.time(7, 10),\n", + " datetime.time(7, 11),\n", + " datetime.time(7, 12),\n", + " datetime.time(7, 13),\n", + " datetime.time(7, 14),\n", + " datetime.time(7, 15),\n", + " datetime.time(7, 16),\n", + " datetime.time(7, 17),\n", + " datetime.time(7, 18),\n", + " datetime.time(7, 19),\n", + " datetime.time(7, 20),\n", + " datetime.time(7, 21),\n", + " datetime.time(7, 22),\n", + " datetime.time(7, 23),\n", + " datetime.time(7, 24),\n", + " datetime.time(7, 25),\n", + " datetime.time(7, 26),\n", + " datetime.time(7, 27),\n", + " datetime.time(7, 28),\n", + " datetime.time(7, 29),\n", + " datetime.time(7, 30),\n", + " datetime.time(7, 31),\n", + " datetime.time(7, 32),\n", + " datetime.time(7, 33),\n", + " datetime.time(7, 34),\n", + " datetime.time(7, 35),\n", + " datetime.time(7, 36),\n", + " datetime.time(7, 37),\n", + " datetime.time(7, 38),\n", + " datetime.time(7, 39),\n", + " datetime.time(7, 40),\n", + " datetime.time(7, 41),\n", + " datetime.time(7, 42),\n", + " datetime.time(7, 43),\n", + " datetime.time(7, 44),\n", + " datetime.time(7, 45),\n", + " datetime.time(7, 46),\n", + " datetime.time(7, 47),\n", + " datetime.time(7, 48),\n", + " datetime.time(7, 49),\n", + " datetime.time(7, 50),\n", + " datetime.time(7, 51),\n", + " datetime.time(7, 52),\n", + " datetime.time(7, 53),\n", + " datetime.time(7, 54),\n", + " datetime.time(7, 55),\n", + " datetime.time(7, 56),\n", + " datetime.time(7, 57),\n", + " datetime.time(7, 58),\n", + " datetime.time(7, 59),\n", + " datetime.time(8, 0),\n", + " datetime.time(8, 1),\n", + " datetime.time(8, 2),\n", + " datetime.time(8, 3),\n", + " datetime.time(8, 4),\n", + " datetime.time(8, 5),\n", + " datetime.time(8, 6),\n", + " datetime.time(8, 7),\n", + " datetime.time(8, 8),\n", + " datetime.time(8, 9),\n", + " datetime.time(8, 10),\n", + " datetime.time(8, 11),\n", + " datetime.time(8, 12),\n", + " datetime.time(8, 13),\n", + " datetime.time(8, 14),\n", + " datetime.time(8, 15),\n", + " datetime.time(8, 16),\n", + " datetime.time(8, 17),\n", + " datetime.time(8, 18),\n", + " datetime.time(8, 19),\n", + " datetime.time(8, 20),\n", + " datetime.time(8, 21),\n", + " datetime.time(8, 22),\n", + " datetime.time(8, 23),\n", + " datetime.time(8, 24),\n", + " datetime.time(8, 25),\n", + " datetime.time(8, 26),\n", + " datetime.time(8, 27),\n", + " datetime.time(8, 28),\n", + " datetime.time(8, 29),\n", + " datetime.time(8, 30),\n", + " datetime.time(8, 31),\n", + " datetime.time(8, 32),\n", + " datetime.time(8, 33),\n", + " datetime.time(8, 34),\n", + " datetime.time(8, 35),\n", + " datetime.time(8, 36),\n", + " datetime.time(8, 37),\n", + " datetime.time(8, 38),\n", + " datetime.time(8, 39),\n", + " datetime.time(8, 40),\n", + " datetime.time(8, 41),\n", + " datetime.time(8, 42),\n", + " datetime.time(8, 43),\n", + " datetime.time(8, 44),\n", + " datetime.time(8, 45),\n", + " datetime.time(8, 46),\n", + " datetime.time(8, 47),\n", + " datetime.time(8, 48),\n", + " datetime.time(8, 49),\n", + " datetime.time(8, 50),\n", + " datetime.time(8, 51),\n", + " datetime.time(8, 52),\n", + " datetime.time(8, 53),\n", + " datetime.time(8, 54),\n", + " datetime.time(8, 55),\n", + " datetime.time(8, 56),\n", + " datetime.time(8, 57),\n", + " datetime.time(8, 58),\n", + " datetime.time(8, 59),\n", + " datetime.time(9, 0),\n", + " datetime.time(9, 1),\n", + " datetime.time(9, 2),\n", + " datetime.time(9, 3),\n", + " datetime.time(9, 4),\n", + " datetime.time(9, 5),\n", + " datetime.time(9, 6),\n", + " datetime.time(9, 7),\n", + " datetime.time(9, 8),\n", + " datetime.time(9, 9),\n", + " datetime.time(9, 10),\n", + " datetime.time(9, 11),\n", + " datetime.time(9, 12),\n", + " datetime.time(9, 13),\n", + " datetime.time(9, 14),\n", + " datetime.time(9, 15),\n", + " datetime.time(9, 16),\n", + " datetime.time(9, 17),\n", + " datetime.time(9, 18),\n", + " datetime.time(9, 19),\n", + " datetime.time(9, 20),\n", + " datetime.time(9, 21),\n", + " datetime.time(9, 22),\n", + " datetime.time(9, 23),\n", + " datetime.time(9, 24),\n", + " datetime.time(9, 25),\n", + " datetime.time(9, 26),\n", + " datetime.time(9, 27),\n", + " datetime.time(9, 28),\n", + " datetime.time(9, 29),\n", + " datetime.time(9, 30),\n", + " datetime.time(9, 31),\n", + " datetime.time(9, 32),\n", + " datetime.time(9, 33),\n", + " datetime.time(9, 34),\n", + " datetime.time(9, 35),\n", + " datetime.time(9, 36),\n", + " datetime.time(9, 37),\n", + " datetime.time(9, 38),\n", + " datetime.time(9, 39),\n", + " datetime.time(9, 40),\n", + " datetime.time(9, 41),\n", + " datetime.time(9, 42),\n", + " datetime.time(9, 43),\n", + " datetime.time(9, 44),\n", + " datetime.time(9, 45),\n", + " datetime.time(9, 46),\n", + " datetime.time(9, 47),\n", + " datetime.time(9, 48),\n", + " datetime.time(9, 49),\n", + " datetime.time(9, 50),\n", + " datetime.time(9, 51),\n", + " datetime.time(9, 52),\n", + " datetime.time(9, 53),\n", + " datetime.time(9, 54),\n", + " datetime.time(9, 55),\n", + " datetime.time(9, 56),\n", + " datetime.time(9, 57),\n", + " datetime.time(9, 58),\n", + " datetime.time(9, 59),\n", + " datetime.time(10, 0),\n", + " datetime.time(10, 1),\n", + " datetime.time(10, 2),\n", + " datetime.time(10, 3),\n", + " datetime.time(10, 4),\n", + " datetime.time(10, 5),\n", + " datetime.time(10, 6),\n", + " datetime.time(10, 7),\n", + " datetime.time(10, 8),\n", + " datetime.time(10, 9),\n", + " datetime.time(10, 10),\n", + " datetime.time(10, 11),\n", + " datetime.time(10, 12),\n", + " datetime.time(10, 13),\n", + " datetime.time(10, 14),\n", + " datetime.time(10, 15),\n", + " datetime.time(10, 16),\n", + " datetime.time(10, 17),\n", + " datetime.time(10, 18),\n", + " datetime.time(10, 19),\n", + " datetime.time(10, 20),\n", + " datetime.time(10, 21),\n", + " datetime.time(10, 22),\n", + " datetime.time(10, 23),\n", + " datetime.time(10, 24),\n", + " datetime.time(10, 25),\n", + " datetime.time(10, 26),\n", + " datetime.time(10, 27),\n", + " datetime.time(10, 28),\n", + " datetime.time(10, 29),\n", + " datetime.time(10, 30),\n", + " datetime.time(10, 31),\n", + " datetime.time(10, 32),\n", + " datetime.time(10, 33),\n", + " datetime.time(10, 34),\n", + " datetime.time(10, 35),\n", + " datetime.time(10, 36),\n", + " datetime.time(10, 37),\n", + " datetime.time(10, 38),\n", + " datetime.time(10, 39),\n", + " datetime.time(10, 40),\n", + " datetime.time(10, 41),\n", + " datetime.time(10, 42),\n", + " datetime.time(10, 43),\n", + " datetime.time(10, 44),\n", + " datetime.time(10, 45),\n", + " datetime.time(10, 46),\n", + " datetime.time(10, 47),\n", + " datetime.time(10, 48),\n", + " datetime.time(10, 49),\n", + " datetime.time(10, 50),\n", + " datetime.time(10, 51),\n", + " datetime.time(10, 52),\n", + " datetime.time(10, 53),\n", + " datetime.time(10, 54),\n", + " datetime.time(10, 55),\n", + " datetime.time(10, 56),\n", + " datetime.time(10, 57),\n", + " datetime.time(10, 58),\n", + " datetime.time(10, 59),\n", + " datetime.time(11, 0),\n", + " datetime.time(11, 1),\n", + " datetime.time(11, 2),\n", + " datetime.time(11, 3),\n", + " datetime.time(11, 4),\n", + " datetime.time(11, 5),\n", + " datetime.time(11, 6),\n", + " datetime.time(11, 7),\n", + " datetime.time(11, 8),\n", + " datetime.time(11, 9),\n", + " datetime.time(11, 10),\n", + " datetime.time(11, 11),\n", + " datetime.time(11, 12),\n", + " datetime.time(11, 13),\n", + " datetime.time(11, 14),\n", + " datetime.time(11, 15),\n", + " datetime.time(11, 16),\n", + " datetime.time(11, 17),\n", + " datetime.time(11, 18),\n", + " datetime.time(11, 19),\n", + " datetime.time(11, 20),\n", + " datetime.time(11, 21),\n", + " datetime.time(11, 22),\n", + " datetime.time(11, 23),\n", + " datetime.time(11, 24),\n", + " datetime.time(11, 25),\n", + " datetime.time(11, 26),\n", + " datetime.time(11, 27),\n", + " datetime.time(11, 28),\n", + " datetime.time(11, 29),\n", + " datetime.time(11, 30),\n", + " datetime.time(11, 31),\n", + " datetime.time(11, 32),\n", + " datetime.time(11, 33),\n", + " datetime.time(11, 34),\n", + " datetime.time(11, 35),\n", + " datetime.time(11, 36),\n", + " datetime.time(11, 37),\n", + " datetime.time(11, 38),\n", + " datetime.time(11, 39),\n", + " datetime.time(11, 40),\n", + " datetime.time(11, 41),\n", + " datetime.time(11, 42),\n", + " datetime.time(11, 43),\n", + " datetime.time(11, 44),\n", + " datetime.time(11, 45),\n", + " datetime.time(11, 46),\n", + " datetime.time(11, 47),\n", + " datetime.time(11, 48),\n", + " datetime.time(11, 49),\n", + " datetime.time(11, 50),\n", + " datetime.time(11, 51),\n", + " datetime.time(11, 52),\n", + " datetime.time(11, 53),\n", + " datetime.time(11, 54),\n", + " datetime.time(11, 55),\n", + " datetime.time(11, 56),\n", + " datetime.time(11, 57),\n", + " datetime.time(11, 58),\n", + " datetime.time(11, 59),\n", + " datetime.time(12, 0),\n", + " datetime.time(12, 1),\n", + " datetime.time(12, 2),\n", + " datetime.time(12, 3),\n", + " datetime.time(12, 4),\n", + " datetime.time(12, 5),\n", + " datetime.time(12, 6),\n", + " datetime.time(12, 7),\n", + " datetime.time(12, 8),\n", + " datetime.time(12, 9),\n", + " datetime.time(12, 10),\n", + " datetime.time(12, 11),\n", + " datetime.time(12, 12),\n", + " datetime.time(12, 13),\n", + " datetime.time(12, 14),\n", + " datetime.time(12, 15),\n", + " datetime.time(12, 16),\n", + " datetime.time(12, 17),\n", + " datetime.time(12, 18),\n", + " datetime.time(12, 19),\n", + " datetime.time(12, 20),\n", + " datetime.time(12, 21),\n", + " datetime.time(12, 22),\n", + " datetime.time(12, 23),\n", + " datetime.time(12, 24),\n", + " datetime.time(12, 25),\n", + " datetime.time(12, 26),\n", + " datetime.time(12, 27),\n", + " datetime.time(12, 28),\n", + " datetime.time(12, 29),\n", + " datetime.time(12, 30),\n", + " datetime.time(12, 31),\n", + " datetime.time(12, 32),\n", + " datetime.time(12, 33),\n", + " datetime.time(12, 34),\n", + " datetime.time(12, 35),\n", + " datetime.time(12, 36),\n", + " datetime.time(12, 37),\n", + " datetime.time(12, 38),\n", + " datetime.time(12, 39),\n", + " datetime.time(12, 40),\n", + " datetime.time(12, 41),\n", + " datetime.time(12, 42),\n", + " datetime.time(12, 43),\n", + " datetime.time(12, 44),\n", + " datetime.time(12, 45),\n", + " datetime.time(12, 46),\n", + " datetime.time(12, 47),\n", + " datetime.time(12, 48),\n", + " datetime.time(12, 49),\n", + " datetime.time(12, 50),\n", + " datetime.time(12, 51),\n", + " datetime.time(12, 52),\n", + " datetime.time(12, 53),\n", + " datetime.time(12, 54),\n", + " datetime.time(12, 55),\n", + " datetime.time(12, 56),\n", + " datetime.time(12, 57),\n", + " datetime.time(12, 58),\n", + " datetime.time(12, 59),\n", + " datetime.time(13, 0),\n", + " datetime.time(13, 1),\n", + " datetime.time(13, 2),\n", + " datetime.time(13, 3),\n", + " datetime.time(13, 4),\n", + " datetime.time(13, 5),\n", + " datetime.time(13, 6),\n", + " datetime.time(13, 7),\n", + " datetime.time(13, 8),\n", + " datetime.time(13, 9),\n", + " datetime.time(13, 10),\n", + " datetime.time(13, 11),\n", + " datetime.time(13, 12),\n", + " datetime.time(13, 13),\n", + " datetime.time(13, 14),\n", + " datetime.time(13, 15),\n", + " datetime.time(13, 16),\n", + " datetime.time(13, 17),\n", + " datetime.time(13, 18),\n", + " datetime.time(13, 19),\n", + " datetime.time(13, 20),\n", + " datetime.time(13, 21),\n", + " datetime.time(13, 22),\n", + " datetime.time(13, 23),\n", + " datetime.time(13, 24),\n", + " datetime.time(13, 25),\n", + " datetime.time(13, 26),\n", + " datetime.time(13, 27),\n", + " datetime.time(13, 28),\n", + " datetime.time(13, 29),\n", + " datetime.time(13, 30),\n", + " datetime.time(13, 31),\n", + " datetime.time(13, 32),\n", + " datetime.time(13, 33),\n", + " datetime.time(13, 34),\n", + " datetime.time(13, 35),\n", + " datetime.time(13, 36),\n", + " datetime.time(13, 37),\n", + " datetime.time(13, 38),\n", + " datetime.time(13, 39),\n", + " datetime.time(13, 40),\n", + " datetime.time(13, 41),\n", + " datetime.time(13, 42),\n", + " datetime.time(13, 43),\n", + " datetime.time(13, 44),\n", + " datetime.time(13, 45),\n", + " datetime.time(13, 46),\n", + " datetime.time(13, 47),\n", + " datetime.time(13, 48),\n", + " datetime.time(13, 49),\n", + " datetime.time(13, 50),\n", + " datetime.time(13, 51),\n", + " datetime.time(13, 52),\n", + " datetime.time(13, 53),\n", + " datetime.time(13, 54),\n", + " datetime.time(13, 55),\n", + " datetime.time(13, 56),\n", + " datetime.time(13, 57),\n", + " datetime.time(13, 58),\n", + " datetime.time(13, 59),\n", + " datetime.time(14, 0),\n", + " datetime.time(14, 1),\n", + " datetime.time(14, 2),\n", + " datetime.time(14, 3),\n", + " datetime.time(14, 4),\n", + " datetime.time(14, 5),\n", + " datetime.time(14, 6),\n", + " datetime.time(14, 7),\n", + " datetime.time(14, 8),\n", + " datetime.time(14, 9),\n", + " datetime.time(14, 10),\n", + " datetime.time(14, 11),\n", + " datetime.time(14, 12),\n", + " datetime.time(14, 13),\n", + " datetime.time(14, 14),\n", + " datetime.time(14, 15),\n", + " datetime.time(14, 16),\n", + " datetime.time(14, 17),\n", + " datetime.time(14, 18),\n", + " datetime.time(14, 19),\n", + " datetime.time(14, 20),\n", + " datetime.time(14, 21),\n", + " datetime.time(14, 22),\n", + " datetime.time(14, 23),\n", + " datetime.time(14, 24),\n", + " datetime.time(14, 25),\n", + " datetime.time(14, 26),\n", + " datetime.time(14, 27),\n", + " datetime.time(14, 28),\n", + " datetime.time(14, 29),\n", + " datetime.time(14, 30),\n", + " datetime.time(14, 31),\n", + " datetime.time(14, 32),\n", + " datetime.time(14, 33),\n", + " datetime.time(14, 34),\n", + " datetime.time(14, 35),\n", + " datetime.time(14, 36),\n", + " datetime.time(14, 37),\n", + " datetime.time(14, 38),\n", + " datetime.time(14, 39),\n", + " datetime.time(14, 40),\n", + " datetime.time(14, 41),\n", + " datetime.time(14, 42),\n", + " datetime.time(14, 43),\n", + " datetime.time(14, 44),\n", + " datetime.time(14, 45),\n", + " datetime.time(14, 46),\n", + " datetime.time(14, 47),\n", + " datetime.time(14, 48),\n", + " datetime.time(14, 49),\n", + " datetime.time(14, 50),\n", + " datetime.time(14, 51),\n", + " datetime.time(14, 52),\n", + " datetime.time(14, 53),\n", + " datetime.time(14, 54),\n", + " datetime.time(14, 55),\n", + " datetime.time(14, 56),\n", + " datetime.time(14, 57),\n", + " datetime.time(14, 58),\n", + " datetime.time(14, 59),\n", + " datetime.time(15, 0),\n", + " datetime.time(15, 1),\n", + " datetime.time(15, 2),\n", + " datetime.time(15, 3),\n", + " datetime.time(15, 4),\n", + " datetime.time(15, 5),\n", + " datetime.time(15, 6),\n", + " datetime.time(15, 7),\n", + " datetime.time(15, 8),\n", + " datetime.time(15, 9),\n", + " datetime.time(15, 10),\n", + " datetime.time(15, 11),\n", + " datetime.time(15, 12),\n", + " datetime.time(15, 13),\n", + " datetime.time(15, 14),\n", + " datetime.time(15, 15),\n", + " datetime.time(15, 16),\n", + " datetime.time(15, 17),\n", + " datetime.time(15, 18),\n", + " datetime.time(15, 19),\n", + " datetime.time(15, 20),\n", + " datetime.time(15, 21),\n", + " datetime.time(15, 22),\n", + " datetime.time(15, 23),\n", + " datetime.time(15, 24),\n", + " datetime.time(15, 25),\n", + " datetime.time(15, 26),\n", + " datetime.time(15, 27),\n", + " datetime.time(15, 28),\n", + " datetime.time(15, 29),\n", + " datetime.time(15, 30),\n", + " datetime.time(15, 31),\n", + " datetime.time(15, 32),\n", + " datetime.time(15, 33),\n", + " datetime.time(15, 34),\n", + " datetime.time(15, 35),\n", + " datetime.time(15, 36),\n", + " datetime.time(15, 37),\n", + " datetime.time(15, 38),\n", + " datetime.time(15, 39),\n", + " datetime.time(15, 40),\n", + " datetime.time(15, 41),\n", + " datetime.time(15, 42),\n", + " datetime.time(15, 43),\n", + " datetime.time(15, 44),\n", + " datetime.time(15, 45),\n", + " datetime.time(15, 46),\n", + " datetime.time(15, 47),\n", + " datetime.time(15, 48),\n", + " datetime.time(15, 49),\n", + " datetime.time(15, 50),\n", + " datetime.time(15, 51),\n", + " datetime.time(15, 52),\n", + " datetime.time(15, 53),\n", + " datetime.time(15, 54),\n", + " datetime.time(15, 55),\n", + " datetime.time(15, 56),\n", + " datetime.time(15, 57),\n", + " datetime.time(15, 58),\n", + " datetime.time(15, 59),\n", + " datetime.time(16, 0),\n", + " datetime.time(16, 1),\n", + " datetime.time(16, 2),\n", + " datetime.time(16, 3),\n", + " datetime.time(16, 4),\n", + " datetime.time(16, 5),\n", + " datetime.time(16, 6),\n", + " datetime.time(16, 7),\n", + " datetime.time(16, 8),\n", + " datetime.time(16, 9),\n", + " datetime.time(16, 10),\n", + " datetime.time(16, 11),\n", + " datetime.time(16, 12),\n", + " datetime.time(16, 13),\n", + " datetime.time(16, 14),\n", + " datetime.time(16, 15),\n", + " datetime.time(16, 16),\n", + " datetime.time(16, 17),\n", + " datetime.time(16, 18),\n", + " datetime.time(16, 19),\n", + " datetime.time(16, 20),\n", + " datetime.time(16, 21),\n", + " datetime.time(16, 22),\n", + " datetime.time(16, 23),\n", + " datetime.time(16, 24),\n", + " datetime.time(16, 25),\n", + " datetime.time(16, 26),\n", + " datetime.time(16, 27),\n", + " datetime.time(16, 28),\n", + " datetime.time(16, 29),\n", + " datetime.time(16, 30),\n", + " datetime.time(16, 31),\n", + " datetime.time(16, 32),\n", + " datetime.time(16, 33),\n", + " datetime.time(16, 34),\n", + " datetime.time(16, 35),\n", + " datetime.time(16, 36),\n", + " datetime.time(16, 37),\n", + " datetime.time(16, 38),\n", + " datetime.time(16, 39),\n", + " ...]" + ] + }, + "execution_count": 77, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "times" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12ed775b", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "603a33a2", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n", + " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n", + " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n", + " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n", + " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n", + " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n", + " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n", + " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n", + " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n", + " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n", + " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n", + " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n", + " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j])" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.zeros((100), dtype=\"complex\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "combined-ethiopia", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}