Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Simon Perkins
ducc
Commits
dc7877d7
Commit
dc7877d7
authored
Jan 16, 2020
by
Martin Reinecke
Browse files
slight source code reorganization
parent
c9c3cd15
Changes
8
Hide whitespace changes
Inline
Side-by-side
libsharp2/sharp.cc
View file @
dc7877d7
...
@@ -197,137 +197,6 @@ struct ringhelper
...
@@ -197,137 +197,6 @@ struct ringhelper
}
}
};
};
sharp_standard_alm_info
::
sharp_standard_alm_info
(
size_t
lmax__
,
size_t
nm_
,
ptrdiff_t
stride_
,
const
size_t
*
mval__
,
const
ptrdiff_t
*
mstart
)
:
lmax_
(
lmax__
),
mval_
(
nm_
),
mvstart
(
nm_
),
stride
(
stride_
)
{
for
(
size_t
mi
=
0
;
mi
<
nm_
;
++
mi
)
{
mval_
[
mi
]
=
mval__
[
mi
];
mvstart
[
mi
]
=
mstart
[
mi
];
}
}
sharp_standard_alm_info
::
sharp_standard_alm_info
(
size_t
lmax__
,
size_t
mmax_
,
ptrdiff_t
stride_
,
const
ptrdiff_t
*
mstart
)
:
lmax_
(
lmax__
),
mval_
(
mmax_
+
1
),
mvstart
(
mmax_
+
1
),
stride
(
stride_
)
{
for
(
size_t
i
=
0
;
i
<=
mmax_
;
++
i
)
{
mval_
[
i
]
=
i
;
mvstart
[
i
]
=
mstart
[
i
];
}
}
ptrdiff_t
sharp_standard_alm_info
::
index
(
int
l
,
int
mi
)
{
return
mvstart
[
mi
]
+
stride
*
l
;
}
sharp_standard_geom_info
::
sharp_standard_geom_info
(
size_t
nrings
,
const
size_t
*
nph
,
const
ptrdiff_t
*
ofs
,
ptrdiff_t
stride_
,
const
double
*
phi0
,
const
double
*
theta
,
const
double
*
wgt
)
:
ring
(
nrings
),
stride
(
stride_
)
{
size_t
pos
=
0
;
nphmax_
=
0
;
for
(
size_t
m
=
0
;
m
<
nrings
;
++
m
)
{
ring
[
m
].
theta
=
theta
[
m
];
ring
[
m
].
cth
=
cos
(
theta
[
m
]);
ring
[
m
].
sth
=
sin
(
theta
[
m
]);
ring
[
m
].
weight
=
(
wgt
!=
nullptr
)
?
wgt
[
m
]
:
1.
;
ring
[
m
].
phi0
=
phi0
[
m
];
ring
[
m
].
ofs
=
ofs
[
m
];
ring
[
m
].
nph
=
nph
[
m
];
if
(
nphmax_
<
nph
[
m
])
nphmax_
=
nph
[
m
];
}
sort
(
ring
.
begin
(),
ring
.
end
(),[](
const
Tring
&
a
,
const
Tring
&
b
)
{
return
(
a
.
sth
<
b
.
sth
);
});
while
(
pos
<
nrings
)
{
pair_
.
push_back
(
Tpair
());
pair_
.
back
().
r1
=
pos
;
if
((
pos
<
nrings
-
1
)
&&
approx
(
ring
[
pos
].
cth
,
-
ring
[
pos
+
1
].
cth
,
1e-12
))
{
if
(
ring
[
pos
].
cth
>
0
)
// make sure northern ring is in r1
pair_
.
back
().
r2
=
pos
+
1
;
else
{
pair_
.
back
().
r1
=
pos
+
1
;
pair_
.
back
().
r2
=
pos
;
}
++
pos
;
}
else
pair_
.
back
().
r2
=
size_t
(
~
0
);
++
pos
;
}
sort
(
pair_
.
begin
(),
pair_
.
end
(),
[
this
]
(
const
Tpair
&
a
,
const
Tpair
&
b
)
{
if
(
ring
[
a
.
r1
].
nph
==
ring
[
b
.
r1
].
nph
)
return
(
ring
[
a
.
r1
].
phi0
<
ring
[
b
.
r1
].
phi0
)
?
true
:
((
ring
[
a
.
r1
].
phi0
>
ring
[
b
.
r1
].
phi0
)
?
false
:
(
ring
[
a
.
r1
].
cth
>
ring
[
b
.
r1
].
cth
));
return
ring
[
a
.
r1
].
nph
<
ring
[
b
.
r1
].
nph
;
});
}
/* This currently requires all m values from 0 to nm-1 to be present.
It might be worthwhile to relax this criterion such that holes in the m
distribution are permissible. */
size_t
sharp_standard_alm_info
::
mmax
()
const
{
//FIXME: if gaps are allowed, we have to search the maximum m in the array
auto
nm_
=
mval_
.
size
();
vector
<
bool
>
mcheck
(
nm_
,
false
);
for
(
auto
m_cur
:
mval_
)
{
MR_assert
(
m_cur
<
nm_
,
"not all m values are present"
);
MR_assert
(
mcheck
[
m_cur
]
==
false
,
"duplicate m value"
);
mcheck
[
m_cur
]
=
true
;
}
return
nm_
-
1
;
}
MRUTIL_NOINLINE
void
sharp_standard_geom_info
::
clear_map
(
double
*
map
)
const
{
for
(
const
auto
&
r
:
ring
)
{
if
(
stride
==
1
)
memset
(
&
map
[
r
.
ofs
],
0
,
r
.
nph
*
sizeof
(
double
));
else
for
(
size_t
i
=
0
;
i
<
r
.
nph
;
++
i
)
map
[
r
.
ofs
+
i
*
stride
]
=
0
;
}
}
MRUTIL_NOINLINE
void
sharp_standard_geom_info
::
clear_map
(
float
*
map
)
const
{
for
(
const
auto
&
r
:
ring
)
{
if
(
stride
==
1
)
memset
(
&
map
[
r
.
ofs
],
0
,
r
.
nph
*
sizeof
(
float
));
else
for
(
size_t
i
=
0
;
i
<
r
.
nph
;
++
i
)
map
[
r
.
ofs
+
i
*
stride
]
=
0
;
}
}
void
sharp_standard_alm_info
::
clear_alm
(
dcmplx
*
alm
)
const
{
for
(
size_t
mi
=
0
;
mi
<
mval_
.
size
();
++
mi
)
for
(
size_t
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
reinterpret_cast
<
dcmplx
*>
(
alm
)[
mvstart
[
mi
]
+
l
*
stride
]
=
0.
;
}
void
sharp_standard_alm_info
::
clear_alm
(
fcmplx
*
alm
)
const
{
for
(
size_t
mi
=
0
;
mi
<
mval_
.
size
();
++
mi
)
for
(
size_t
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
reinterpret_cast
<
fcmplx
*>
(
alm
)[
mvstart
[
mi
]
+
l
*
stride
]
=
0.
;
}
MRUTIL_NOINLINE
void
sharp_job
::
init_output
()
MRUTIL_NOINLINE
void
sharp_job
::
init_output
()
{
{
if
(
flags
&
SHARP_ADD
)
return
;
if
(
flags
&
SHARP_ADD
)
return
;
...
@@ -365,27 +234,6 @@ void sharp_job::alloc_almtmp (size_t lmax, vector<dcmplx> &data)
...
@@ -365,27 +234,6 @@ void sharp_job::alloc_almtmp (size_t lmax, vector<dcmplx> &data)
almtmp
=
data
.
data
();
almtmp
=
data
.
data
();
}
}
void
sharp_standard_alm_info
::
get_alm
(
size_t
mi
,
const
dcmplx
*
alm
,
dcmplx
*
almtmp
,
size_t
nalm
)
const
{
for
(
auto
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
almtmp
[
nalm
*
l
]
=
alm
[
mvstart
[
mi
]
+
l
*
stride
];
}
void
sharp_standard_alm_info
::
get_alm
(
size_t
mi
,
const
fcmplx
*
alm
,
dcmplx
*
almtmp
,
size_t
nalm
)
const
{
for
(
auto
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
almtmp
[
nalm
*
l
]
=
alm
[
mvstart
[
mi
]
+
l
*
stride
];
}
void
sharp_standard_alm_info
::
add_alm
(
size_t
mi
,
const
dcmplx
*
almtmp
,
dcmplx
*
alm
,
size_t
nalm
)
const
{
for
(
auto
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
alm
[
mvstart
[
mi
]
+
l
*
stride
]
+=
almtmp
[
nalm
*
l
];
}
void
sharp_standard_alm_info
::
add_alm
(
size_t
mi
,
const
dcmplx
*
almtmp
,
fcmplx
*
alm
,
size_t
nalm
)
const
{
for
(
auto
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
alm
[
mvstart
[
mi
]
+
l
*
stride
]
+=
fcmplx
(
almtmp
[
nalm
*
l
]);
}
MRUTIL_NOINLINE
void
sharp_job
::
alm2almtmp
(
size_t
lmax
,
size_t
mi
)
MRUTIL_NOINLINE
void
sharp_job
::
alm2almtmp
(
size_t
lmax
,
size_t
mi
)
{
{
if
(
type
!=
SHARP_MAP2ALM
)
if
(
type
!=
SHARP_MAP2ALM
)
...
@@ -439,39 +287,6 @@ MRUTIL_NOINLINE void sharp_job::almtmp2alm (size_t lmax, size_t mi)
...
@@ -439,39 +287,6 @@ MRUTIL_NOINLINE void sharp_job::almtmp2alm (size_t lmax, size_t mi)
ainfo
->
add_alm
(
mi
,
almtmp
+
i
,
reinterpret_cast
<
fcmplx
**>
(
alm
)[
i
],
nalm
);
ainfo
->
add_alm
(
mi
,
almtmp
+
i
,
reinterpret_cast
<
fcmplx
**>
(
alm
)[
i
],
nalm
);
}
}
//virtual
void
sharp_standard_geom_info
::
add_ring
(
bool
weighted
,
size_t
iring
,
const
double
*
ringtmp
,
double
*
map
)
const
{
double
*
MRUTIL_RESTRICT
p1
=&
map
[
ring
[
iring
].
ofs
];
double
wgt
=
weighted
?
ring
[
iring
].
weight
:
1.
;
for
(
size_t
m
=
0
;
m
<
ring
[
iring
].
nph
;
++
m
)
p1
[
m
*
stride
]
+=
ringtmp
[
m
]
*
wgt
;
}
//virtual
void
sharp_standard_geom_info
::
add_ring
(
bool
weighted
,
size_t
iring
,
const
double
*
ringtmp
,
float
*
map
)
const
{
float
*
MRUTIL_RESTRICT
p1
=&
map
[
ring
[
iring
].
ofs
];
double
wgt
=
weighted
?
ring
[
iring
].
weight
:
1.
;
for
(
size_t
m
=
0
;
m
<
ring
[
iring
].
nph
;
++
m
)
p1
[
m
*
stride
]
+=
float
(
ringtmp
[
m
]
*
wgt
);
}
//virtual
void
sharp_standard_geom_info
::
get_ring
(
bool
weighted
,
size_t
iring
,
const
double
*
map
,
double
*
ringtmp
)
const
{
const
double
*
MRUTIL_RESTRICT
p1
=&
map
[
ring
[
iring
].
ofs
];
double
wgt
=
weighted
?
ring
[
iring
].
weight
:
1.
;
for
(
size_t
m
=
0
;
m
<
ring
[
iring
].
nph
;
++
m
)
ringtmp
[
m
]
=
p1
[
m
*
stride
]
*
wgt
;
}
//virtual
void
sharp_standard_geom_info
::
get_ring
(
bool
weighted
,
size_t
iring
,
const
float
*
map
,
double
*
ringtmp
)
const
{
const
float
*
MRUTIL_RESTRICT
p1
=&
map
[
ring
[
iring
].
ofs
];
double
wgt
=
weighted
?
ring
[
iring
].
weight
:
1.
;
for
(
size_t
m
=
0
;
m
<
ring
[
iring
].
nph
;
++
m
)
ringtmp
[
m
]
=
p1
[
m
*
stride
]
*
wgt
;
}
MRUTIL_NOINLINE
void
sharp_job
::
ringtmp2ring
(
size_t
iring
,
MRUTIL_NOINLINE
void
sharp_job
::
ringtmp2ring
(
size_t
iring
,
const
vector
<
double
>
&
ringtmp
,
ptrdiff_t
rstride
)
const
vector
<
double
>
&
ringtmp
,
ptrdiff_t
rstride
)
{
{
...
...
libsharp2/sharp.h
View file @
dc7877d7
...
@@ -59,52 +59,6 @@ class sharp_geom_info
...
@@ -59,52 +59,6 @@ class sharp_geom_info
virtual
void
add_ring
(
bool
weighted
,
size_t
iring
,
const
double
*
ringtmp
,
float
*
map
)
const
=
0
;
virtual
void
add_ring
(
bool
weighted
,
size_t
iring
,
const
double
*
ringtmp
,
float
*
map
)
const
=
0
;
};
};
class
sharp_standard_geom_info
:
public
sharp_geom_info
{
private:
struct
Tring
{
double
theta
,
phi0
,
weight
,
cth
,
sth
;
ptrdiff_t
ofs
;
size_t
nph
;
};
std
::
vector
<
Tring
>
ring
;
std
::
vector
<
Tpair
>
pair_
;
ptrdiff_t
stride
;
size_t
nphmax_
;
public:
/*! Creates a geometry information from a set of ring descriptions.
All arrays passed to this function must have \a nrings elements.
\param nrings the number of rings in the map
\param nph the number of pixels in each ring
\param ofs the index of the first pixel in each ring in the map array
\param stride the stride between consecutive pixels
\param phi0 the azimuth (in radians) of the first pixel in each ring
\param theta the colatitude (in radians) of each ring
\param wgt the pixel weight to be used for the ring in map2alm
and adjoint map2alm transforms.
Pass nullptr to use 1.0 as weight for all rings.
*/
sharp_standard_geom_info
(
size_t
nrings
,
const
size_t
*
nph_
,
const
ptrdiff_t
*
ofs
,
ptrdiff_t
stride_
,
const
double
*
phi0_
,
const
double
*
theta_
,
const
double
*
wgt
);
virtual
size_t
nrings
()
const
{
return
ring
.
size
();
}
virtual
size_t
npairs
()
const
{
return
pair_
.
size
();
}
virtual
size_t
nph
(
size_t
iring
)
const
{
return
ring
[
iring
].
nph
;
}
virtual
size_t
nphmax
()
const
{
return
nphmax_
;
}
virtual
double
theta
(
size_t
iring
)
const
{
return
ring
[
iring
].
theta
;
}
virtual
double
cth
(
size_t
iring
)
const
{
return
ring
[
iring
].
cth
;
}
virtual
double
sth
(
size_t
iring
)
const
{
return
ring
[
iring
].
sth
;
}
virtual
double
phi0
(
size_t
iring
)
const
{
return
ring
[
iring
].
phi0
;
}
virtual
Tpair
pair
(
size_t
ipair
)
const
{
return
pair_
[
ipair
];
}
virtual
void
clear_map
(
double
*
map
)
const
;
virtual
void
clear_map
(
float
*
map
)
const
;
virtual
void
get_ring
(
bool
weighted
,
size_t
iring
,
const
double
*
map
,
double
*
ringtmp
)
const
;
virtual
void
get_ring
(
bool
weighted
,
size_t
iring
,
const
float
*
map
,
double
*
ringtmp
)
const
;
virtual
void
add_ring
(
bool
weighted
,
size_t
iring
,
const
double
*
ringtmp
,
double
*
map
)
const
;
virtual
void
add_ring
(
bool
weighted
,
size_t
iring
,
const
double
*
ringtmp
,
float
*
map
)
const
;
};
/*! \defgroup almgroup Helpers for dealing with a_lm */
/*! \defgroup almgroup Helpers for dealing with a_lm */
/*! \{ */
/*! \{ */
...
@@ -123,60 +77,6 @@ class sharp_alm_info
...
@@ -123,60 +77,6 @@ class sharp_alm_info
virtual
void
add_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
almtmp
,
std
::
complex
<
double
>
*
alm
,
size_t
nalm
)
const
=
0
;
virtual
void
add_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
almtmp
,
std
::
complex
<
double
>
*
alm
,
size_t
nalm
)
const
=
0
;
virtual
void
add_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
almtmp
,
std
::
complex
<
float
>
*
alm
,
size_t
nalm
)
const
=
0
;
virtual
void
add_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
almtmp
,
std
::
complex
<
float
>
*
alm
,
size_t
nalm
)
const
=
0
;
};
};
/*! \internal
Helper type for index calculation in a_lm arrays. */
class
sharp_standard_alm_info
:
public
sharp_alm_info
{
private:
/*! Maximum \a l index of the array */
size_t
lmax_
;
/*! Array with \a nm entries containing the individual m values */
std
::
vector
<
size_t
>
mval_
;
/*! Array with \a nm entries containing the (hypothetical) indices of
the coefficients with quantum numbers 0,\a mval[i] */
std
::
vector
<
ptrdiff_t
>
mvstart
;
/*! Stride between a_lm and a_(l+1),m */
ptrdiff_t
stride
;
public:
/*! Creates an a_lm data structure from the following parameters:
\param lmax maximum \a l quantum number (>=0)
\param mmax maximum \a m quantum number (0<= \a mmax <= \a lmax)
\param stride the stride between entries with identical \a m, and \a l
differing by 1.
\param mstart the index of the (hypothetical) coefficient with the
quantum numbers 0,\a m. Must have \a mmax+1 entries.
*/
sharp_standard_alm_info
(
size_t
lmax__
,
size_t
mmax_
,
ptrdiff_t
stride_
,
const
ptrdiff_t
*
mstart
);
/*! Creates an a_lm data structure which from the following parameters:
\param lmax maximum \a l quantum number (\a >=0)
\param nm number of different \a m (\a 0<=nm<=lmax+1)
\param stride the stride between entries with identical \a m, and \a l
differing by 1.
\param mval array with \a nm entries containing the individual m values
\param mvstart array with \a nm entries containing the (hypothetical)
indices of the coefficients with the quantum numbers 0,\a mval[i]
*/
sharp_standard_alm_info
(
size_t
lmax__
,
size_t
nm__
,
ptrdiff_t
stride_
,
const
size_t
*
mval__
,
const
ptrdiff_t
*
mvstart_
);
/*! Returns the index of the coefficient with quantum numbers \a l,
\a mval_[mi].
\note for a \a sharp_alm_info generated by sharp_make_alm_info() this is
the index for the coefficient with the quantum numbers \a l, \a mi. */
ptrdiff_t
index
(
int
l
,
int
mi
);
virtual
size_t
lmax
()
const
{
return
lmax_
;
}
virtual
size_t
mmax
()
const
;
virtual
size_t
nm
()
const
{
return
mval_
.
size
();
}
virtual
size_t
mval
(
size_t
i
)
const
{
return
mval_
[
i
];
}
virtual
void
clear_alm
(
std
::
complex
<
double
>
*
alm
)
const
;
virtual
void
clear_alm
(
std
::
complex
<
float
>
*
alm
)
const
;
virtual
void
get_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
alm
,
std
::
complex
<
double
>
*
almtmp
,
size_t
nalm
)
const
;
virtual
void
get_alm
(
size_t
mi
,
const
std
::
complex
<
float
>
*
alm
,
std
::
complex
<
double
>
*
almtmp
,
size_t
nalm
)
const
;
virtual
void
add_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
almtmp
,
std
::
complex
<
double
>
*
alm
,
size_t
nalm
)
const
;
virtual
void
add_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
almtmp
,
std
::
complex
<
float
>
*
alm
,
size_t
nalm
)
const
;
};
/*! \} */
/*! \} */
...
...
libsharp2/sharp_almhelpers.cc
View file @
dc7877d7
...
@@ -25,9 +25,87 @@
...
@@ -25,9 +25,87 @@
* \author Martin Reinecke
* \author Martin Reinecke
*/
*/
#include
"mr_util/error_handling.h"
#include
"libsharp2/sharp_almhelpers.h"
#include
"libsharp2/sharp_almhelpers.h"
using
namespace
std
;
using
namespace
std
;
using
dcmplx
=
complex
<
double
>
;
using
fcmplx
=
complex
<
float
>
;
sharp_standard_alm_info
::
sharp_standard_alm_info
(
size_t
lmax__
,
size_t
nm_
,
ptrdiff_t
stride_
,
const
size_t
*
mval__
,
const
ptrdiff_t
*
mstart
)
:
lmax_
(
lmax__
),
mval_
(
nm_
),
mvstart
(
nm_
),
stride
(
stride_
)
{
for
(
size_t
mi
=
0
;
mi
<
nm_
;
++
mi
)
{
mval_
[
mi
]
=
mval__
[
mi
];
mvstart
[
mi
]
=
mstart
[
mi
];
}
}
sharp_standard_alm_info
::
sharp_standard_alm_info
(
size_t
lmax__
,
size_t
mmax_
,
ptrdiff_t
stride_
,
const
ptrdiff_t
*
mstart
)
:
lmax_
(
lmax__
),
mval_
(
mmax_
+
1
),
mvstart
(
mmax_
+
1
),
stride
(
stride_
)
{
for
(
size_t
i
=
0
;
i
<=
mmax_
;
++
i
)
{
mval_
[
i
]
=
i
;
mvstart
[
i
]
=
mstart
[
i
];
}
}
void
sharp_standard_alm_info
::
clear_alm
(
dcmplx
*
alm
)
const
{
for
(
size_t
mi
=
0
;
mi
<
mval_
.
size
();
++
mi
)
for
(
size_t
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
reinterpret_cast
<
dcmplx
*>
(
alm
)[
mvstart
[
mi
]
+
l
*
stride
]
=
0.
;
}
void
sharp_standard_alm_info
::
clear_alm
(
fcmplx
*
alm
)
const
{
for
(
size_t
mi
=
0
;
mi
<
mval_
.
size
();
++
mi
)
for
(
size_t
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
reinterpret_cast
<
fcmplx
*>
(
alm
)[
mvstart
[
mi
]
+
l
*
stride
]
=
0.
;
}
void
sharp_standard_alm_info
::
get_alm
(
size_t
mi
,
const
dcmplx
*
alm
,
dcmplx
*
almtmp
,
size_t
nalm
)
const
{
for
(
auto
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
almtmp
[
nalm
*
l
]
=
alm
[
mvstart
[
mi
]
+
l
*
stride
];
}
void
sharp_standard_alm_info
::
get_alm
(
size_t
mi
,
const
fcmplx
*
alm
,
dcmplx
*
almtmp
,
size_t
nalm
)
const
{
for
(
auto
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
almtmp
[
nalm
*
l
]
=
alm
[
mvstart
[
mi
]
+
l
*
stride
];
}
void
sharp_standard_alm_info
::
add_alm
(
size_t
mi
,
const
dcmplx
*
almtmp
,
dcmplx
*
alm
,
size_t
nalm
)
const
{
for
(
auto
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
alm
[
mvstart
[
mi
]
+
l
*
stride
]
+=
almtmp
[
nalm
*
l
];
}
void
sharp_standard_alm_info
::
add_alm
(
size_t
mi
,
const
dcmplx
*
almtmp
,
fcmplx
*
alm
,
size_t
nalm
)
const
{
for
(
auto
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
alm
[
mvstart
[
mi
]
+
l
*
stride
]
+=
fcmplx
(
almtmp
[
nalm
*
l
]);
}
ptrdiff_t
sharp_standard_alm_info
::
index
(
int
l
,
int
mi
)
{
return
mvstart
[
mi
]
+
stride
*
l
;
}
/* This currently requires all m values from 0 to nm-1 to be present.
It might be worthwhile to relax this criterion such that holes in the m
distribution are permissible. */
size_t
sharp_standard_alm_info
::
mmax
()
const
{
//FIXME: if gaps are allowed, we have to search the maximum m in the array
auto
nm_
=
mval_
.
size
();
vector
<
bool
>
mcheck
(
nm_
,
false
);
for
(
auto
m_cur
:
mval_
)
{
MR_assert
(
m_cur
<
nm_
,
"not all m values are present"
);
MR_assert
(
mcheck
[
m_cur
]
==
false
,
"duplicate m value"
);
mcheck
[
m_cur
]
=
true
;
}
return
nm_
-
1
;
}
unique_ptr
<
sharp_standard_alm_info
>
sharp_make_triangular_alm_info
(
int
lmax
,
int
mmax
,
int
stride
)
unique_ptr
<
sharp_standard_alm_info
>
sharp_make_triangular_alm_info
(
int
lmax
,
int
mmax
,
int
stride
)
{
{
...
...
libsharp2/sharp_almhelpers.h
View file @
dc7877d7
...
@@ -32,6 +32,61 @@
...
@@ -32,6 +32,61 @@
#include
<memory>
#include
<memory>
#include
"libsharp2/sharp.h"
#include
"libsharp2/sharp.h"
/*! \internal
Helper type for index calculation in a_lm arrays. */
class
sharp_standard_alm_info
:
public
sharp_alm_info
{
private:
/*! Maximum \a l index of the array */
size_t
lmax_
;
/*! Array with \a nm entries containing the individual m values */
std
::
vector
<
size_t
>
mval_
;
/*! Array with \a nm entries containing the (hypothetical) indices of
the coefficients with quantum numbers 0,\a mval[i] */
std
::
vector
<
ptrdiff_t
>
mvstart
;
/*! Stride between a_lm and a_(l+1),m */
ptrdiff_t
stride
;
public:
/*! Creates an a_lm data structure from the following parameters:
\param lmax maximum \a l quantum number (>=0)
\param mmax maximum \a m quantum number (0<= \a mmax <= \a lmax)
\param stride the stride between entries with identical \a m, and \a l
differing by 1.
\param mstart the index of the (hypothetical) coefficient with the
quantum numbers 0,\a m. Must have \a mmax+1 entries.
*/
sharp_standard_alm_info
(
size_t
lmax__
,
size_t
mmax_
,
ptrdiff_t
stride_
,
const
ptrdiff_t
*
mstart
);
/*! Creates an a_lm data structure which from the following parameters:
\param lmax maximum \a l quantum number (\a >=0)
\param nm number of different \a m (\a 0<=nm<=lmax+1)
\param stride the stride between entries with identical \a m, and \a l
differing by 1.
\param mval array with \a nm entries containing the individual m values
\param mvstart array with \a nm entries containing the (hypothetical)
indices of the coefficients with the quantum numbers 0,\a mval[i]
*/
sharp_standard_alm_info
(
size_t
lmax__
,
size_t
nm__
,
ptrdiff_t
stride_
,
const
size_t
*
mval__
,
const
ptrdiff_t
*
mvstart_
);
/*! Returns the index of the coefficient with quantum numbers \a l,
\a mval_[mi].
\note for a \a sharp_alm_info generated by sharp_make_alm_info() this is
the index for the coefficient with the quantum numbers \a l, \a mi. */
ptrdiff_t
index
(
int
l
,
int
mi
);
virtual
size_t
lmax
()
const
{
return
lmax_
;
}
virtual
size_t
mmax
()
const
;
virtual
size_t
nm
()
const
{
return
mval_
.
size
();
}
virtual
size_t
mval
(
size_t
i
)
const
{
return
mval_
[
i
];
}
virtual
void
clear_alm
(
std
::
complex
<
double
>
*
alm
)
const
;
virtual
void
clear_alm
(
std
::
complex
<
float
>
*
alm
)
const
;
virtual
void
get_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
alm
,
std
::
complex
<
double
>
*
almtmp
,
size_t
nalm
)
const
;
virtual
void
get_alm
(
size_t
mi
,
const
std
::
complex
<
float
>
*
alm
,
std
::
complex
<
double
>
*
almtmp
,
size_t
nalm
)
const
;
virtual
void
add_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
almtmp
,
std
::
complex
<
double
>
*
alm
,
size_t
nalm
)
const
;
virtual
void
add_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
almtmp
,
std
::
complex
<
float
>
*
alm
,
size_t
nalm
)
const
;
};
/*! Initialises an a_lm data structure according to the scheme used by
/*! Initialises an a_lm data structure according to the scheme used by
Healpix_cxx.
Healpix_cxx.
\ingroup almgroup */
\ingroup almgroup */
...
...
libsharp2/sharp_geomhelpers.cc
View file @
dc7877d7
...
@@ -31,8 +31,118 @@
...
@@ -31,8 +31,118 @@
#include
"mr_util/gl_integrator.h"
#include
"mr_util/gl_integrator.h"
#include
"mr_util/fft.h"
#include
"mr_util/fft.h"
#include
"mr_util/error_handling.h"
#include
"mr_util/error_handling.h"
#include
"mr_util/math_utils.h"
using
namespace
std
;
using
namespace
std
;
using
namespace
mr
;
sharp_standard_geom_info
::
sharp_standard_geom_info
(
size_t
nrings
,
const
size_t
*
nph
,
const
ptrdiff_t
*
ofs
,
ptrdiff_t
stride_
,
const
double
*
phi0
,
const
double
*
theta
,
const
double
*
wgt
)
:
ring
(
nrings
),
stride
(
stride_
)
{
size_t
pos
=
0
;
nphmax_
=
0
;
for
(
size_t
m
=
0
;
m
<
nrings
;
++
m
)
{
ring
[
m
].
theta
=
theta
[
m
];
ring
[
m
].
cth
=
cos
(
theta
[
m
]);
ring
[
m
].
sth
=
sin
(
theta
[
m
]);
ring
[
m
].
weight
=
(
wgt
!=
nullptr
)
?
wgt
[
m
]
:
1.
;
ring
[
m
].
phi0
=
phi0
[
m
];
ring
[
m
].
ofs
=
ofs
[
m
];
ring
[
m
].
nph
=
nph
[
m
];
if
(
nphmax_
<
nph
[
m
])
nphmax_
=
nph
[
m
];
}
sort
(
ring
.
begin
(),
ring
.
end
(),[](
const
Tring
&
a
,
const
Tring
&
b
)
{
return
(
a
.
sth
<
b
.
sth
);
});
while
(
pos
<
nrings
)
{
pair_
.
push_back
(
Tpair
());
pair_
.
back
().
r1
=
pos
;
if
((
pos
<
nrings
-
1
)
&&
approx
(
ring
[
pos
].
cth
,
-
ring
[
pos
+
1
].
cth
,
1e-12
))
{
if
(
ring
[
pos
].
cth
>
0
)
// make sure northern ring is in r1
pair_
.
back
().
r2
=
pos
+
1
;
else
{
pair_
.
back
().
r1
=
pos
+
1
;
pair_
.
back
().
r2
=
pos
;
}
++
pos
;
}
else
pair_
.
back
().
r2
=
size_t
(
~
0
);
++
pos
;
}
sort
(
pair_
.
begin
(),
pair_
.
end
(),
[
this
]
(
const
Tpair
&
a
,
const
Tpair
&
b
)
{
if
(
ring
[
a
.
r1
].
nph
==
ring
[
b
.
r1
].
nph
)
return
(
ring
[
a
.
r1
].
phi0
<
ring
[
b
.
r1
].
phi0
)
?
true
:
((
ring
[
a
.
r1
].
phi0
>
ring
[
b
.
r1
].
phi0
)
?
false
:
(
ring
[
a
.
r1
].
cth
>
ring
[
b
.
r1
].
cth
));
return
ring
[
a
.
r1
].
nph
<
ring
[
b
.
r1
].
nph
;
});
}
void
sharp_standard_geom_info
::
clear_map
(
double
*
map
)
const
{
for
(
const
auto
&
r
:
ring
)
{
if
(
stride
==
1
)
memset
(
&
map
[
r
.
ofs
],
0
,
r
.
nph
*
sizeof
(
double
));
else
for
(
size_t
i
=
0
;
i
<
r
.
nph
;
++
i
)
map
[
r
.
ofs
+
i
*
stride
]
=
0
;
}
}
void
sharp_standard_geom_info
::
clear_map
(
float
*
map
)
const
{
for
(
const
auto
&
r
:
ring
)
{
if
(
stride
==
1
)
memset
(
&
map
[
r
.
ofs
],
0
,
r
.
nph
*
sizeof
(
float
));
else
for
(
size_t
i
=
0
;
i
<
r
.
nph
;
++
i
)
map
[
r
.
ofs
+
i
*
stride
]
=
0
;
}
}
//virtual
void
sharp_standard_geom_info
::
add_ring
(
bool
weighted
,
size_t
iring
,
const
double
*
ringtmp
,
double
*
map
)
const
{
double
*
MRUTIL_RESTRICT
p1
=&
map
[
ring
[
iring
].