Skip to content

Commit

Permalink
Merge pull request #258 from juliantaylor/perf-improvements
Browse files Browse the repository at this point in the history 
various coordinate performance improvements
  • Loading branch information
@dpetry
dpetry committed Oct 27, 2015
2 parents 8c0df8c + fdd1928 commit e7e765e
Show file tree
Hide file tree
Showing 12 changed files with 281 additions and 59 deletions.
58 changes: 45 additions & 13 deletions casa/Quanta/Euler.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,14 +38,44 @@ namespace casacore { //# NAMESPACE CASACORE - BEGIN

// Euler class

// simplistic Vector(3) cache to reduce allocation overhead for temporaries
#if defined(AIPS_CXX11) && !defined(__APPLE__)
thread_local size_t Euler::available = 0;
thread_local Euler::DataArrays Euler::arrays[50];
#endif

Euler::DataArrays Euler::get_arrays()
{
#if defined(AIPS_CXX11) && !defined(__APPLE__)
if (available > 0) {
return arrays[--available];
}
#endif
return std::make_pair(new Vector<Double>(3), new Vector<Int>(3));
}

void Euler::return_arrays(Euler::DataArrays array)
{
#if defined(AIPS_CXX11) && !defined(__APPLE__)
if (available < sizeof(arrays) / sizeof(arrays[0]) &&
array.first->size() == 3 && array.second->size() == 3 &&
array.first->nrefs() == 1 && array.second->nrefs() == 1) {
arrays[available++] = array;
return;
}
#endif
delete array.first;
delete array.second;
}

//# Constructors
Euler::Euler() : euler(3), axes(3) {
Euler::Euler() : data(get_arrays()), euler(*data.first), axes(*data.second) {
euler = Double(0.0);
indgen(axes,1,1);
}

Euler::Euler(const Euler &other) :
euler(3), axes(3) {
data(get_arrays()), euler(*data.first), axes(*data.second) {
euler = other.euler;
axes = other.axes;
}
Expand All @@ -59,7 +89,7 @@ Euler &Euler::operator=(const Euler &other) {
}

Euler::Euler(Double in0, Double in1, Double in2) :
euler(3), axes(3){
data(get_arrays()), euler(*data.first), axes(*data.second) {
euler(0) = in0;
euler(1) = in1;
euler(2) = in2;
Expand All @@ -69,7 +99,7 @@ euler(3), axes(3){

Euler::Euler(Double in0, uInt ax0, Double in1, uInt ax1, Double in2,
uInt ax2) :
euler(3), axes(3) {
data(get_arrays()), euler(*data.first), axes(*data.second) {
DebugAssert(ax0 <= 3 && ax1 <=3 && ax2 <=3, AipsError);
euler(0) = in0;
euler(1) = in1;
Expand All @@ -80,31 +110,31 @@ euler(3), axes(3) {
}

Euler::Euler(const Quantity &in0) :
euler(3), axes(3) {
data(get_arrays()), euler(*data.first), axes(*data.second) {
euler(0) = Euler::makeRad(in0);
euler(1) = 0;
euler(2) = 0;
indgen(axes,1,1);
}

Euler::Euler(const Quantity &in0, const Quantity &in1) :
euler(3), axes(3) {
data(get_arrays()), euler(*data.first), axes(*data.second) {
euler(0) = Euler::makeRad(in0);
euler(1) = Euler::makeRad(in1);
euler(2) = 0;
indgen(axes,1,1);
}

Euler::Euler(const Quantity &in0, const Quantity &in1, const Quantity &in2) :
euler(3), axes(3) {
data(get_arrays()), euler(*data.first), axes(*data.second) {
euler(0) = Euler::makeRad(in0);
euler(1) = Euler::makeRad(in1);
euler(2) = Euler::makeRad(in2);
indgen(axes,1,1);
}

Euler::Euler(const Quantity &in0, uInt ax0) :
euler(3), axes(3) {
data(get_arrays()), euler(*data.first), axes(*data.second) {
DebugAssert(ax0 <= 3, AipsError);
euler(0) = Euler::makeRad(in0);
euler(1) = 0;
Expand All @@ -114,7 +144,7 @@ euler(3), axes(3) {
axes(2) = 0;
}
Euler::Euler(const Quantity &in0, uInt ax0, const Quantity &in1, uInt ax1) :
euler(3), axes(3) {
data(get_arrays()), euler(*data.first), axes(*data.second) {
DebugAssert(ax0 <= 3 && ax1 <=3, AipsError);
euler(0) = Euler::makeRad(in0);
euler(1) = Euler::makeRad(in1);
Expand All @@ -125,7 +155,7 @@ euler(3), axes(3) {
}
Euler::Euler(const Quantity &in0, uInt ax0, const Quantity &in1, uInt ax1,
const Quantity &in2, uInt ax2) :
euler(3), axes(3) {
data(get_arrays()), euler(*data.first), axes(*data.second) {
DebugAssert(ax0 <= 3 && ax1 <=3 && ax2 <=3, AipsError);
euler(0) = Euler::makeRad(in0);
euler(1) = Euler::makeRad(in1);
Expand All @@ -136,7 +166,7 @@ euler(3), axes(3) {
}

Euler::Euler(const Quantum<Vector<Double> > &in) :
euler(3), axes(3) {
data(get_arrays()), euler(*data.first), axes(*data.second) {
Int i;
Vector<Double> tmp = Euler::makeRad(in);
Int j=tmp.size(); j=min(j,3);
Expand All @@ -150,7 +180,7 @@ euler(3), axes(3) {
}

Euler::Euler(const Quantum<Vector<Double> > &in, const Vector<uInt> &ax) :
euler(3), axes(3) {
data(get_arrays()), euler(*data.first), axes(*data.second) {
Vector<Double> tmp = Euler::makeRad(in);
Int j=tmp.size(); j=min(j,3); Int i=ax.size(); j=min(j,i);
for (i=0; i<j; i++) {
Expand All @@ -165,7 +195,9 @@ euler(3), axes(3) {
}

//# Destructor
Euler::~Euler() {}
Euler::~Euler() {
return_arrays(data);
}

//# Operators
Euler Euler::operator-() const {
Expand Down
18 changes: 14 additions & 4 deletions casa/Quanta/Euler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@
#include <casacore/casa/aips.h>
#include <casacore/casa/Arrays/Vector.h>
#include <casacore/casa/Quanta/Quantum.h>
#include <utility>

namespace casacore { //# NAMESPACE CASACORE - BEGIN

Expand Down Expand Up @@ -189,17 +190,26 @@ class Euler

private:
//# Data
// Actual vector with 3 Euler angles
Vector<Double> euler;
// Axes
Vector<Int> axes;
typedef std::pair<Vector<Double> *, Vector<Int> *> DataArrays;
// data container
DataArrays data;
// vector with 3 Euler angles (data.first)
Vector<Double> & euler;
// Axes (data.second)
Vector<Int> & axes;

//# Private Member Functions
// The makeRad functions check and convert the input Quantities to radians
// <group>
static Double makeRad(const Quantity &in);
static Vector<Double> makeRad(const Quantum<Vector<Double> > &in);
// </group>
DataArrays get_arrays();
void return_arrays(DataArrays array);
#if defined(AIPS_CXX11) && !defined(__APPLE__)
static thread_local DataArrays arrays[50];
static thread_local size_t available;
#endif
};


Expand Down
5 changes: 3 additions & 2 deletions casa/Quanta/MVEpoch.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@ namespace casacore { //# NAMESPACE CASACORE - BEGIN

//# Constants
const Double MVEpoch::secInDay(3600*24);
const Unit MVEpoch::unitDay("d");

//# Constructors
MVEpoch::MVEpoch() :
Expand Down Expand Up @@ -178,7 +179,7 @@ Double MVEpoch::get() const {
}

Quantity MVEpoch::getTime() const {
return (Quantity(get(), "d"));
return (Quantity(get(), unitDay));
}

Quantity MVEpoch::getTime(const Unit &unit) const {
Expand Down Expand Up @@ -249,7 +250,7 @@ Bool MVEpoch::putValue(const Vector<Quantum<Double> > &in) {

Double MVEpoch::makeDay(const Quantity &in) const {
in.assure(UnitVal::TIME);
return in.get("d").getValue();
return in.get(unitDay).getValue();
}

void MVEpoch::addTime(Double in) {
Expand Down
1 change: 1 addition & 0 deletions casa/Quanta/MVEpoch.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,6 +141,7 @@ class MVEpoch : public MeasValue {
//# General Member Functions
// Constants
static const Double secInDay;
static const Unit unitDay;

// Tell me your type
// <group>
Expand Down
70 changes: 50 additions & 20 deletions casa/Quanta/MVPosition.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,15 +46,44 @@ namespace casacore { //# NAMESPACE CASACORE - BEGIN
const Double MVPosition::loLimit = 743.568;
const Double MVPosition::hiLimit = 743.569;

// simplistic vector(3) cache to reduce allocation overhead for temporaries
#if defined(AIPS_CXX11) && !defined(__APPLE__)
static thread_local size_t available = 0;
static thread_local Vector<Double> * arrays[50];
#endif

Vector<Double> * get_array()
{
#if defined(AIPS_CXX11) && !defined(__APPLE__)
if (available > 0) {
return arrays[--available];
}
#endif
return new Vector<Double>(3);
}

void return_array(Vector<Double> * array)
{
#if defined(AIPS_CXX11) && !defined(__APPLE__)
if (available < sizeof(arrays) / sizeof(arrays[0]) &&
array->size() == 3 &&
array->nrefs() == 1) {
arrays[available++] = array;
return;
}
#endif
delete array;
}

//# Constructors
MVPosition::MVPosition() :
xyz(3) {
xyz(*get_array()) {
xyz = Double(0.0);
}

MVPosition::MVPosition(const MVPosition &other) :
MeasValue(),
xyz(3)
xyz(*get_array())
{
xyz = other.xyz;
}
Expand All @@ -67,27 +96,27 @@ MVPosition &MVPosition::operator=(const MVPosition &other) {
}

MVPosition::MVPosition(Double in) :
xyz(3) {
xyz(*get_array()) {
xyz = Double(0.0);
xyz(2) = in;
}

MVPosition::MVPosition(const Quantity &l) :
xyz(3) {
xyz(*get_array()) {
xyz = Double(0.0);
l.assure(UnitVal::LENGTH);
xyz(2) = l.getBaseValue();
}

MVPosition::MVPosition(Double in0, Double in1, Double in2) :
xyz(3) {
xyz(*get_array()) {
xyz(0) = in0;
xyz(1) = in1;
xyz(2) = in2;
}

MVPosition::MVPosition(const Quantity &l, Double angle0, Double angle1) :
xyz(3) {
xyz(*get_array()) {
Double loc = std::cos(angle1);
xyz(0) = std::cos(angle0)*loc;
xyz(1) = std::sin(angle0)*loc;
Expand All @@ -100,7 +129,7 @@ MVPosition::MVPosition(const Quantity &l, Double angle0, Double angle1) :

MVPosition::MVPosition(const Quantity &l, const Quantity &angle0,
const Quantity &angle1) :
xyz(3) {
xyz(*get_array()) {
Double loc = (cos(angle1)).getValue();
xyz(0) = ((cos(angle0)).getValue()) * loc;
xyz(1) = ((sin(angle0)).getValue()) * loc;
Expand All @@ -113,7 +142,7 @@ MVPosition::MVPosition(const Quantity &l, const Quantity &angle0,
}

MVPosition::MVPosition(const Quantum<Vector<Double> > &angle) :
xyz(3) {
xyz(*get_array()) {
uInt i; i = angle.getValue().nelements();
if (i > 3 ) {
throw (AipsError("Illegeal vector length in MVPosition constructor"));
Expand All @@ -139,7 +168,7 @@ MVPosition::MVPosition(const Quantum<Vector<Double> > &angle) :

MVPosition::MVPosition(const Quantity &l,
const Quantum<Vector<Double> > &angle) :
xyz(3) {
xyz(*get_array()) {
uInt i; i = angle.getValue().nelements();
if (i > 3 ) {
throw (AipsError("Illegal vector length in MVPosition constructor"));
Expand All @@ -166,7 +195,7 @@ MVPosition::MVPosition(const Quantity &l,
}

MVPosition::MVPosition(const Vector<Double> &other) :
xyz(3) {
xyz(*get_array()) {
uInt i; i = other.nelements();
if (i > 3 ) {
throw (AipsError("Illegal vector length in MVPosition constructor"));
Expand All @@ -190,14 +219,17 @@ MVPosition::MVPosition(const Vector<Double> &other) :
}

MVPosition::MVPosition(const Vector<Quantity> &other) :
xyz(3) {
xyz(*get_array()) {
if (!putValue(other)) {
throw (AipsError("Illegal quantum vector in MVPosition constructor"));
}
}

//# Destructor
MVPosition::~MVPosition() {}
MVPosition::~MVPosition()
{
return_array(&xyz);
}

//# Operators
Bool MVPosition::
Expand Down Expand Up @@ -273,14 +305,12 @@ MVPosition MVPosition::operator-(const MVPosition &right) const{
}

MVPosition &MVPosition::operator*=(const RotMatrix &right) {
MVPosition result;
for (Int i=0; i<3; i++) {
result(i) = 0;
for (Int j=0; j<3; j++) {
result(i) += xyz(j) * right(j,i);
}
}
*this = result;
Double x = xyz(0);
Double y = xyz(1);
Double z = xyz(2);
xyz(0) = x * right(0, 0) + y * right(1, 0) + z * right(2, 0);
xyz(1) = x * right(0, 1) + y * right(1, 1) + z * right(2, 1);
xyz(2) = x * right(0, 2) + y * right(1, 2) + z * right(2, 2);
return *this;
}

Expand Down
4 changes: 2 additions & 2 deletions casa/Quanta/MVPosition.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -155,7 +155,7 @@ class MVPosition : public MeasValue {
MVPosition &operator=(const MVPosition &other);

// Destructor
~MVPosition();
virtual ~MVPosition();

//# Operators
// Multiplication defined as in-product
Expand Down Expand Up @@ -279,7 +279,7 @@ class MVPosition : public MeasValue {
Double getLat(Double ln) const;
//# Data
// Position vector (in m)
Vector<Double> xyz;
Vector<Double> & xyz;
};

//# Global functions
Expand Down
Loading

0 comments on commit e7e765e

Please sign in to comment.