Commit e06ca06f authored by Martin Reinecke's avatar Martin Reinecke

not yet working

parent 2f2bbc3e
Pipeline #78085 failed with stages
in 14 minutes and 13 seconds
......@@ -4,6 +4,7 @@
#include "ducc0/math/pointing.cc"
#include "ducc0/math/geom_utils.cc"
#include "ducc0/math/space_filling.cc"
#include "ducc0/math/least_misfit.cc"
#include "ducc0/sharp/sharp.cc"
#include "ducc0/sharp/sharp_almhelpers.cc"
#include "ducc0/sharp/sharp_core.cc"
......
......@@ -36,8 +36,8 @@
#include "ducc0/infra/useful_macros.h"
#include "ducc0/infra/mav.h"
#include "ducc0/infra/simd.h"
#include "ducc0/math/es_kernel.h"
#include "ducc0/math/gridding_kernel.h"
#include "ducc0/math/least_misfit.h"
namespace ducc0 {
......@@ -202,15 +202,18 @@ template<typename T> class GridderConfig
{
protected:
size_t nx_dirty, ny_dirty, nu, nv;
double ofactor, eps, psx, psy;
double ofactor;
public:
shared_ptr<GriddingKernel<T>> krn;
protected:
double psx, psy;
size_t supp, nsafe;
double beta;
size_t nthreads;
double ushift, vshift;
int maxiu0, maxiv0;
public:
shared_ptr<GriddingKernel<T>> krn;
protected:
complex<T> wscreen(T x, T y, T w, bool adjoint) const
{
......@@ -228,14 +231,12 @@ template<typename T> class GridderConfig
double epsilon, double pixsize_x, double pixsize_y, size_t nthreads_)
: nx_dirty(nxdirty), ny_dirty(nydirty), nu(nu_), nv(nv_),
ofactor(min(double(nu)/nxdirty, double(nv)/nydirty)),
eps(epsilon),
krn(make_shared<HornerKernel<T>>(selectLeastMisfitKernel(ofactor, epsilon))),
psx(pixsize_x), psy(pixsize_y),
supp(ES_Kernel::get_supp(epsilon, ofactor)), nsafe((supp+1)/2),
beta(ES_Kernel::get_beta(supp, ofactor)*supp),
supp(krn->support()), nsafe((supp+1)/2),
nthreads(nthreads_),
ushift(supp*(-0.5)+1+nu), vshift(supp*(-0.5)+1+nv),
maxiu0((nu+nsafe)-supp), maxiv0((nv+nsafe)-supp),
krn(make_shared<HornerKernel<T>>(supp, supp+3, [this](double v){return double(esk(v,double(beta)));}))
maxiu0((nu+nsafe)-supp), maxiv0((nv+nsafe)-supp)
{
MR_assert(nu>=2*nsafe, "nu too small");
MR_assert(nv>=2*nsafe, "nv too small");
......@@ -246,8 +247,6 @@ template<typename T> class GridderConfig
MR_assert(epsilon>0, "epsilon must be positive");
MR_assert(pixsize_x>0, "pixsize_x must be positive");
MR_assert(pixsize_y>0, "pixsize_y must be positive");
MR_assert(ofactor>=1.175,
"oversampling factor too small (>=1.2 recommended)");
}
GridderConfig(size_t nxdirty, size_t nydirty,
double epsilon, double pixsize_x, double pixsize_y, size_t nthreads_)
......@@ -256,7 +255,6 @@ template<typename T> class GridderConfig
pixsize_y, nthreads_) {}
size_t Nxdirty() const { return nx_dirty; }
size_t Nydirty() const { return ny_dirty; }
double Epsilon() const { return eps; }
double Pixsize_x() const { return psx; }
double Pixsize_y() const { return psy; }
size_t Nu() const { return nu; }
......@@ -264,7 +262,6 @@ template<typename T> class GridderConfig
size_t Supp() const { return supp; }
size_t Nsafe() const { return nsafe; }
size_t Nthreads() const { return nthreads; }
double Ofactor() const{ return ofactor; }
void grid2dirty_post(mav<T,2> &tmav,
mav<T,2> &dirty) const
......
......@@ -52,60 +52,60 @@ def explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, xpixsize, ypixsize,
res += (ms[row, chan]*wgt[row, chan]
* np.exp(2j*np.pi*phase)).real
return res/n
#
#
# @pmp("nxdirty", (30, 128))
# @pmp("nydirty", (128, 250))
# @pmp("ofactor", (1.2, 1.5, 1.7, 2.0))
# @pmp("nrow", (2, 27))
# @pmp("nchan", (1, 5))
# @pmp("epsilon", (1e-1, 1e-3, 1e-5))
# @pmp("singleprec", (True, False))
# @pmp("wstacking", (True, False))
# @pmp("use_wgt", (True, False))
# @pmp("nthreads", (1, 2))
# def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
# singleprec, wstacking, use_wgt, nthreads):
# if singleprec and epsilon < 5e-5:
# return
# rng = np.random.default_rng(42)
# pixsizex = np.pi/180/60/nxdirty*0.2398
# pixsizey = np.pi/180/60/nxdirty
# speedoflight, f0 = 299792458., 1e9
# freq = f0 + np.arange(nchan)*(f0/nchan)
# uvw = (rng.random((nrow, 3))-0.5)/(pixsizey*f0/speedoflight)
# ms = rng.random((nrow, nchan))-0.5 + 1j*(rng.random((nrow, nchan))-0.5)
# wgt = rng.random((nrow, nchan)) if use_wgt else None
# dirty = rng.random((nxdirty, nydirty))-0.5
# nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1
# if nu & 1:
# nu += 1
# if nv & 1:
# nv += 1
# if singleprec:
# ms = ms.astype("c8")
# dirty = dirty.astype("f4")
# if wgt is not None:
# wgt = wgt.astype("f4")
# dirty2 = ng.ms2dirty(uvw, freq, ms, wgt, nxdirty, nydirty, pixsizex,
# pixsizey, nu, nv, epsilon, wstacking, nthreads, 0).astype("f8")
# ms2 = ng.dirty2ms(uvw, freq, dirty, wgt, pixsizex, pixsizey, nu, nv, epsilon,
# wstacking, nthreads, 0).astype("c16")
# tol = 5e-4 if singleprec else 5e-11
# assert_allclose(np.vdot(ms, ms2).real, np.vdot(dirty2, dirty), rtol=tol)
@pmp("nxdirty", (30, 128))
@pmp("nydirty", (128, 250))
@pmp("ofactor", (1.2, 1.5, 1.7, 2.0))
@pmp("nrow", (2, 27))
@pmp("nchan", (1, 5))
@pmp("epsilon", (1e-1, 1e-3, 1e-5))
@pmp("singleprec", (True, False))
@pmp("wstacking", (True, False))
@pmp("use_wgt", (True, False))
@pmp("nthreads", (1, 2))
def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
singleprec, wstacking, use_wgt, nthreads):
if singleprec and epsilon < 5e-5:
return
rng = np.random.default_rng(42)
pixsizex = np.pi/180/60/nxdirty*0.2398
pixsizey = np.pi/180/60/nxdirty
speedoflight, f0 = 299792458., 1e9
freq = f0 + np.arange(nchan)*(f0/nchan)
uvw = (rng.random((nrow, 3))-0.5)/(pixsizey*f0/speedoflight)
ms = rng.random((nrow, nchan))-0.5 + 1j*(rng.random((nrow, nchan))-0.5)
wgt = rng.random((nrow, nchan)) if use_wgt else None
dirty = rng.random((nxdirty, nydirty))-0.5
nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1
if nu & 1:
nu += 1
if nv & 1:
nv += 1
if singleprec:
ms = ms.astype("c8")
dirty = dirty.astype("f4")
if wgt is not None:
wgt = wgt.astype("f4")
dirty2 = ng.ms2dirty(uvw, freq, ms, wgt, nxdirty, nydirty, pixsizex,
pixsizey, nu, nv, epsilon, wstacking, nthreads, 0).astype("f8")
ms2 = ng.dirty2ms(uvw, freq, dirty, wgt, pixsizex, pixsizey, nu, nv, epsilon,
wstacking, nthreads, 0).astype("c16")
tol = 5e-4 if singleprec else 5e-11
assert_allclose(np.vdot(ms, ms2).real, np.vdot(dirty2, dirty), rtol=tol)
@pmp('nxdirty', [16, 64])
@pmp('nxdirty', [16])
@pmp('nydirty', [64])
@pmp('ofactor', [1.2, 1.4, 1.7, 2])
@pmp("nrow", (2, 27))
@pmp("nchan", (1, 5))
@pmp("epsilon", (1e-2, 1e-3, 1e-4, 1e-7))
@pmp('ofactor', [1.2, 2])
@pmp("nrow", (2,))
@pmp("nchan", (1,))
@pmp("epsilon", ( 1e-2, 1e-4, 1e-7, 1e-9, 1e-12))
@pmp("singleprec", (False,))
@pmp("wstacking", (True,))
@pmp("use_wgt", (True,))
@pmp("nthreads", (1, 2))
@pmp("fov", (1., 20.))
@pmp("wstacking", (False,))
@pmp("use_wgt", (False,))
@pmp("nthreads", (1, ))
@pmp("fov", (1., ))
def test_ms2dirty_against_wdft2(nxdirty, nydirty, ofactor, nrow, nchan, epsilon, singleprec, wstacking, use_wgt, fov, nthreads):
if singleprec and epsilon < 5e-5:
return
......
......@@ -21,7 +21,7 @@ include_dirs = ['.', './src/',
pybind11.get_include(True),
pybind11.get_include(False)]
extra_compile_args = ['--std=c++17', '-march=native', '-ffast-math', '-O3']
extra_compile_args = ['--std=c++17', '-march=native', '-ffast-math', '-O3', '-g']
python_module_link_args = []
......
......@@ -25,6 +25,7 @@
#include <vector>
#include "ducc0/infra/simd.h"
#include "ducc0/math/gl_integrator.h"
#include "ducc0/math/least_misfit.h"
namespace ducc0 {
......@@ -84,7 +85,9 @@ template<typename T> class GriddingKernel
{
public:
virtual ~GriddingKernel() {}
// virtual size_t W() const = 0;
virtual size_t support() const = 0;
/*! Returns the function approximation at W different locations with the
abscissas x, x+2./W, x+4./W, ..., x+(2.*W-2)/W.
x must lie in [-1; -1+2./W].
......@@ -143,7 +146,6 @@ class KernelCorrection
while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i)
res[i] = corfunc(i*dx);
});
cout <<"bleep" << n << " " << res.size() << endl;
return res;
}
};
......@@ -216,6 +218,33 @@ template<typename T> class HornerKernel: public GriddingKernel<T>
}
return coeff;
}
static vector<Tsimd> coeffFromPolyKernel(const PolyKernel &krn)
{
auto W = krn.W;
auto D = krn.coeff.size()/((W+1)/2)-1;
auto nvec = ((W+vlen-1)/vlen);
vector<Tsimd> coeff(nvec*(D+1), 0);
vector<double> coeff_raw(W*(D+1), 0);
size_t Whalf = W/2;
for (size_t j=0; j<=D; ++j)
{
double flip = (j&1) ? -1:1;
for (size_t i=Whalf; i<W; ++i)
{
double val = krn.coeff[(i-Whalf)*(D+1)+j];
coeff_raw[j*W+i] = val;
coeff_raw[j*W+W-1-i] = flip*val;
}
}
for (size_t j=0; j<=D; ++j)
{
for (size_t i=0; i<W; ++i)
coeff[(D-j)*nvec + i/vlen][i%vlen] = T(coeff_raw[j*W+i]);
for (size_t i=W; i<vlen*nvec; ++i)
coeff[(D-j)*nvec + i/vlen][i%vlen] = T(0);
}
return coeff;
}
public:
using GriddingKernel<T>::eval;
......@@ -228,6 +257,14 @@ template<typename T> class HornerKernel: public GriddingKernel<T>
corr(W_, [this](T v){return eval_single(v);})
{}
HornerKernel(const PolyKernel &krn)
: W(krn.W), D(krn.coeff.size()/((W+1)/2)-1), nvec((W+vlen-1)/vlen),
coeff(coeffFromPolyKernel(krn)), evalfunc(evfhelper1<1>()),
corr(W, [this](T v){return eval_single(v);})
{}
virtual size_t support() const { return W; }
virtual void eval(T x, native_simd<T> *res) const
{ (this->*evalfunc)(x, res); }
/*! Returns the function approximation at location x.
......
This source diff could not be displayed because it is too large. You can view the blob instead.
#ifndef DUCC0_LEAST_MISFIT_H
#define DUCC0_LEAST_MISFIT_H
/*
* This file is part of nifty_gridder.
*
* nifty_gridder is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* nifty_gridder is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with nifty_gridder; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/* Copyright (C) 2020 Max-Planck-Society
Author: Martin Reinecke */
#include <vector>
//#include "least_misfit.h"
namespace ducc0 {
namespace detail_least_misfit {
using namespace std;
struct PolyKernel
{
size_t W;
double ofactor;
size_t D;
double epsilon;
vector<double> coeff;
};
const PolyKernel &selectLeastMisfitKernel(double ofactor, double epsilon);
} // namespace detail_least_misfit
using detail_least_misfit::PolyKernel;
using detail_least_misfit::selectLeastMisfitKernel;
} // namespace ducc0
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment