[Debian-astro-commits] [gyoto] 48/221: Support di and dj arguments, route yorick Scenery ray-tracing through Scenery::rayTrace when possible. The case that remains is when i_idx or j_idx is a list.
Thibaut Jean-Claude Paumard
thibaut at moszumanska.debian.org
Fri May 22 20:52:32 UTC 2015
This is an automated email from the git hooks/post-receive script.
thibaut pushed a commit to branch master
in repository gyoto.
commit 8c1a676e776387e3617c7468f17f7c361d244dfe
Author: Thibaut Paumard <paumard at users.sourceforge.net>
Date: Tue Oct 14 16:19:52 2014 +0200
Support di and dj arguments, route yorick Scenery ray-tracing through Scenery::rayTrace when possible. The case that remains is when i_idx or j_idx is a list.
---
bin/gyoto.1 | 7 +++-
bin/gyoto.C | 13 ++++++
include/GyotoAstrobj.h | 7 ++++
lib/Astrobj.C | 4 ++
lib/Scenery.C | 107 +++++++++++++++++++++++++++++++++++--------------
yorick/gyoto_Idx.C | 35 +++++++++++-----
yorick/gyoto_Scenery.C | 10 +++++
yorick/ygyoto_idx.h | 22 +++++-----
8 files changed, 156 insertions(+), 49 deletions(-)
diff --git a/bin/gyoto.1 b/bin/gyoto.1
index cece1d0..b4af1c6 100644
--- a/bin/gyoto.1
+++ b/bin/gyoto.1
@@ -6,7 +6,8 @@
Gyoto \- the General relativitY Orbit Tracer of Observatoire de Paris
.SH SYNOPSIS
gyoto [\fB\-\-silent\fR|\fB\-\-quiet\fR|\fB\-\-verbose\fR[=\fIN\fR]|\fB\-\-debug\fR]
- [\fB\-\-imin\fR=\fIi0\fR] [\fB\-\-imax\fR=\fIi1\fR] [\fB\-\-jmin\fR=\fIj0\fR] [\fB\-\-jmax\fR=\fIj1\fR]
+ [\fB\-\-imin\fR=\fIi0\fR] [\fB\-\-imax\fR=\fIi1\fR] [\fB\-\-di\fR=\fIdi\fR]
+ [\fB\-\-jmin\fR=\fIj0\fR] [\fB\-\-jmax\fR=\fIj1\fR] [\fB\-\-dj\fR=\fIdj\fR]
[\fB\-\-time\fR=\fItobs\fR] [\fB\-\-tmin\fR=\fItmin\fR]
[\fB\-\-fov\fR=\fIangle\fR] [\fB\-\-resolution\fR=\fInpix\fR] [\fB\-\-distance\fR=\fIdist\fR]
[\fB\-\-paln\fR=\fIOmega\fR] [\fB\-\-inclination\fR=\fIi\fR] [\fB\-\-argument\fR=\fItheta\fR]
@@ -63,10 +64,14 @@ pixel of the complete image has coordinates i=1 and j=1.
Default value: 1.
.IP \fB\-\-imax\fR=\fIi1
Default value: \fInpix\fR (see option \fB\-\-resolution\fR below).
+.IP \fB\-\-di\fR=\fIdi
+Default value:1.
.IP \fB\-\-jmin\fR=\fIj0
Default value: 1.
.IP \fB\-\-jmax\fR=\fIj1
Default value: \fInpix\fR (see option \fB\-\-resolution\fR below).
+.IP \fB\-\-dj\fR=\fIdj
+Default value:1.
.SS Setting the camera position
The following parameters are normally provided in the Screen section
diff --git a/bin/gyoto.C b/bin/gyoto.C
index 8fe4142..b3433f9 100644
--- a/bin/gyoto.C
+++ b/bin/gyoto.C
@@ -89,6 +89,7 @@ int main(int argc, char** argv) {
string param;
size_t imin=1, imax=1000000000, jmin=1, jmax=1000000000;
+ ptrdiff_t di=1, dj=1;
// double tobs, tmin, fov, dist, paln, incl, arg;
double tobs=0., tmin=0., fov=0., dist=0., paln=0., incl=0., arg=0.;
size_t res=0, nthreads=0, nprocs=0;
@@ -148,6 +149,14 @@ int main(int argc, char** argv) {
return 1;
}
}
+ else if (param.substr(0,5)=="--di=") {
+ di=atoi(param.substr(5).c_str());
+ if (di==0) di=1;
+ }
+ else if (param.substr(0,5)=="--dj=") {
+ dj=atoi(param.substr(5).c_str());
+ if (dj==0) dj=1;
+ }
else if (param.substr(0,15)=="--impact-coords") {
if (param.size() > 16 && param.substr(15,1)=="=")
ipctfile=param.substr(16);
@@ -324,6 +333,10 @@ int main(int argc, char** argv) {
// Allocate space for the output data
data = new Astrobj::Properties();
+ data->alloc=
+ Astrobj::Properties::entire_field |
+ Astrobj::Properties::all_rows;
+ data->di=di; data->dj=dj;
size_t curquant=0;
size_t offset=res*res;
diff --git a/include/GyotoAstrobj.h b/include/GyotoAstrobj.h
index ef258ce..c99aa2c 100644
--- a/include/GyotoAstrobj.h
+++ b/include/GyotoAstrobj.h
@@ -737,6 +737,13 @@ class Gyoto::Astrobj::Properties : protected Gyoto::SmartPointee {
*/
Gyoto::SmartPointer<Gyoto::Units::Converter> binspectrum_converter_ ;
# endif
+
+ typedef int alloc_t;
+ enum alloc_values {minimal=0, entire_field=1, all_rows=2};
+ alloc_t alloc;
+ ptrdiff_t di;
+ ptrdiff_t dj;
+
public:
Properties(); ///< Default constructor (everything is set to NULL);
Properties (double*, double*); ///<< Set intensity and time pointers.
diff --git a/lib/Astrobj.C b/lib/Astrobj.C
index ae382d6..0a3da93 100644
--- a/lib/Astrobj.C
+++ b/lib/Astrobj.C
@@ -487,6 +487,8 @@ Astrobj::Properties::Properties() :
, intensity_converter_(NULL), spectrum_converter_(NULL),
binspectrum_converter_(NULL)
# endif
+ , alloc(Properties::entire_field | Properties::all_rows)
+ , di(1), dj(1)
{}
Astrobj::Properties::Properties( double * I, double * t) :
@@ -499,6 +501,8 @@ Astrobj::Properties::Properties( double * I, double * t) :
, intensity_converter_(NULL), spectrum_converter_(NULL),
binspectrum_converter_(NULL)
# endif
+ , alloc(Properties::entire_field | Properties::all_rows)
+ , di(1), dj(1)
{}
void Astrobj::Properties::init(size_t nbnuobs) {
diff --git a/lib/Scenery.C b/lib/Scenery.C
index 1fcb55e..4c2ec17 100644
--- a/lib/Scenery.C
+++ b/lib/Scenery.C
@@ -132,7 +132,7 @@ typedef struct SceneryThreadWorkerArg {
pthread_mutex_t * mutex;
pthread_t * parent;
#endif
- size_t i, j, imin, imax, jmin, jmax;
+ size_t i, j, imin, imax, jmin, jmax, eol_offset;
Scenery *sc;
Photon * ph;
Astrobj::Properties *data;
@@ -160,8 +160,11 @@ static void * SceneryThreadWorker (void *arg) {
// local variables to store our parameters
size_t i, j;
- size_t eol_offset =
- larg->sc->screen()->resolution() - larg->imax + larg->imin -1;
+
+ // end-of-line offset, depends on whether memory was allocated for
+ // the full field or only for the imin-imax, jmin:jmax window
+ size_t eol_offset = larg -> eol_offset;
+
Astrobj::Properties data;
double * impactcoords = NULL;
@@ -189,13 +192,13 @@ static void * SceneryThreadWorker (void *arg) {
data = *larg->data;
// update i, j and pointer
- ++larg->i;
- ++(*larg->data);
+ larg->i += data.di;
+ *larg->data += data.alloc&Astrobj::Properties::all_rows?data.di:1;
if (larg->impactcoords) {
- impactcoords = larg->impactcoords; larg->impactcoords+=16;
+ impactcoords = larg->impactcoords; larg->impactcoords+=16*data.di;
}
if (larg->i > larg->imax) {
- ++larg->j; larg->i=larg->imin;
+ larg->j+=data.dj; larg->i=larg->imin;
(*larg->data) += eol_offset;
if (larg->impactcoords) larg->impactcoords+=16*eol_offset;
}
@@ -269,6 +272,33 @@ void Scenery::rayTrace(size_t imin, size_t imax,
if (data) setPropertyConverters(data);
+ // curcell will be the current cell, i.e. where to store the
+ // result of the next computation. It's initial value depends on
+ // how memory was allocated in data:
+ //
+ // data->alloc & Properties::entire_field means memory was
+ // allocated for the entire field. Else, memory was only
+ // allocated for the imin:imax, jmin:jmax portion.
+ //
+ // eol_offset is the additional offset at the end of a line. It is 0
+ // if memory was allocated only for the window.
+ size_t curcell=0, eol_offset=0;
+ if (data && data->alloc & Astrobj::Properties::entire_field) {
+ curcell = (jmin-1)*npix+imin-1;
+ eol_offset = data->dj*npix - data->di*((imax-imin)/data->di+1);
+ }
+
+ if (data) {
+ int toto=debug();
+ debug(1);
+ GYOTO_DEBUG_EXPR(data->di);
+ GYOTO_DEBUG_EXPR(data->dj);
+ GYOTO_DEBUG_EXPR(data->alloc);
+ GYOTO_DEBUG_EXPR(curcell);
+ GYOTO_DEBUG_EXPR(eol_offset);
+ GYOTO_DEBUG_EXPR((imax-imin)%data->di);
+ debug(toto);
+ }
#ifdef HAVE_MPI
if (mpi_team_) {
// We are in an MPI content, either the manager or a worker.
@@ -359,7 +389,7 @@ void Scenery::rayTrace(size_t imin, size_t imax,
// First tell the workers to join our task force
// The corresponding recv is in gyoto-scenery-worker.c
- size_t *ijr = new size_t[working*2+2];
+ vector<size_t> cell(working+1);
while (working) {
// receive one result, need to track back where it belongs and
@@ -374,33 +404,32 @@ void Scenery::rayTrace(size_t imin, size_t imax,
w = s.source();
if (s.tag()==Scenery::raytrace_done && data) {
- size_t cell=(ijr[2*w+1]-1)*npix+ijr[2*w]-1;
// Copy each relevant quantity, performing conversion if needed
if (data->intensity)
- data->intensity[cell]=data->intensity_converter_?
+ data->intensity[cell[w]]=data->intensity_converter_?
(*data->intensity_converter_)(*locdata->intensity):
*locdata->intensity;
- if (data->time) data->time[cell]=*locdata->time;
- if (data->distance) data->distance[cell]=*locdata->distance;
- if (data->first_dmin) data->first_dmin[cell]=*locdata->first_dmin;
- if (data->redshift) data->redshift[cell]=*locdata->redshift;
+ if (data->time) data->time[cell[w]]=*locdata->time;
+ if (data->distance) data->distance[cell[w]]=*locdata->distance;
+ if (data->first_dmin) data->first_dmin[cell[w]]=*locdata->first_dmin;
+ if (data->redshift) data->redshift[cell[w]]=*locdata->redshift;
if (data->impactcoords)
for (size_t k=0; k<16; ++k)
- data->impactcoords[cell*16+k]=locdata->impactcoords[k];
- if (data->user1) data->user1[cell]=*locdata->user1;
- if (data->user2) data->user2[cell]=*locdata->user2;
- if (data->user3) data->user3[cell]=*locdata->user3;
- if (data->user4) data->user4[cell]=*locdata->user4;
- if (data->user5) data->user5[cell]=*locdata->user5;
+ data->impactcoords[cell[w]*16+k]=locdata->impactcoords[k];
+ if (data->user1) data->user1[cell[w]]=*locdata->user1;
+ if (data->user2) data->user2[cell[w]]=*locdata->user2;
+ if (data->user3) data->user3[cell[w]]=*locdata->user3;
+ if (data->user4) data->user4[cell[w]]=*locdata->user4;
+ if (data->user5) data->user5[cell[w]]=*locdata->user5;
if (data->spectrum)
for (size_t c=0; c<nbnuobs; ++c)
- data->spectrum[cell+c*data->offset]=
+ data->spectrum[cell[w]+c*data->offset]=
data->spectrum_converter_?
(*data->spectrum_converter_)(locdata->spectrum[c]):
locdata->spectrum[c];
if (data->binspectrum)
for (size_t c=0; c<nbnuobs; ++c)
- data->binspectrum[cell+c*data->offset]=
+ data->binspectrum[cell[w]+c*data->offset]=
data->binspectrum_converter_?
(*data->binspectrum_converter_)(locdata->binspectrum[c]):
locdata->binspectrum[c];
@@ -408,25 +437,42 @@ void Scenery::rayTrace(size_t imin, size_t imax,
// give new task or decrease working counter
if (ij[0]<=imax) {
- ijr[2*w]=ij[0];
- ijr[2*w+1]=ij[1];
- size_t cell=(ijr[2*w+1]-1)*npix+ijr[2*w]-1;
+ cell[w]=curcell; // store curcell
mpi_team_ -> send(w, raytrace, ij, 2);
+ if (data) {
+ int toto=debug();
+ debug(1);
+ GYOTO_DEBUG_EXPR(data->di);
+ GYOTO_DEBUG_EXPR(data->dj);
+ GYOTO_DEBUG_EXPR(data->alloc);
+ GYOTO_DEBUG_EXPR(curcell);
+ GYOTO_DEBUG_EXPR(eol_offset);
+ GYOTO_DEBUG_EXPR((imax-imin)%data->di);
+ debug(toto);
+ }
+
if (impactcoords) {
- mpi_team_ -> send(w, Scenery::impactcoords, impactcoords+cell*16, 16);
+ mpi_team_ -> send(w, Scenery::impactcoords, impactcoords+cell[w]*16, 16);
}
- if (++ij[0]>imax) {
- cout << "\rj = " << ij[1] << " / " << jmax << " " << flush;
- if (++ij[1]<=jmax) ij[0]=imin;
+ curcell += (data&&data->alloc&Astrobj::Properties::all_rows)?
+ data->di:1;
+
+ if ( (ij[0]+=data?data->di:1) >imax) {
+ if (verbose())
+ cout << "\rj = " << ij[1] << " / " << jmax << " " << flush;
+ if ( (ij[1]+=data?data->dj:1) <=jmax) {
+ ij[0]=imin;
+ curcell += eol_offset;
+ }
}
} else {
mpi_team_ -> send(w, raytrace_done, ij, 2);
--working;
}
}
- delete [] ijr;
+ if (verbose()) cout << endl;
} else {
// We are a worker, do we need to query for impactcoords?
double ipct[16];
@@ -473,6 +519,7 @@ void Scenery::rayTrace(size_t imin, size_t imax,
larg.imax=imax;
larg.jmin=jmin;
larg.jmax=jmax;
+ larg.eol_offset=eol_offset;
struct timeval tim;
double start, end;
diff --git a/yorick/gyoto_Idx.C b/yorick/gyoto_Idx.C
index ca6ca46..cbbea97 100644
--- a/yorick/gyoto_Idx.C
+++ b/yorick/gyoto_Idx.C
@@ -19,20 +19,23 @@
#include "ygyoto_idx.h"
#include "yapi.h"
+#include "GyotoError.h"
using namespace YGyoto;
+using namespace Gyoto;
-int YGyoto::Idx::isNuller() {return _is_nuller;}
-long YGyoto::Idx::getNElements() {return _nel;}
-long YGyoto::Idx::current() {
+int YGyoto::Idx::isNuller() const {return _is_nuller;}
+long YGyoto::Idx::getNElements() const {return _nel;}
+long YGyoto::Idx::current() const {
if (_is_list) return _idx[_cur];
return _cur;
}
-double YGyoto::Idx::getDVal() {return _is_double?_dval:_range[0];}
-int YGyoto::Idx::isDouble() {return _is_double;}
-int YGyoto::Idx::isFirst() {return _is_first;}
+double YGyoto::Idx::getDVal() const {return _is_double?_dval:_range[0];}
+int YGyoto::Idx::isDouble() const {return _is_double;}
+int YGyoto::Idx::isRangeOrScalar() const {return _is_range || _is_scalar;}
+int YGyoto::Idx::isFirst() const {return _is_first;}
-int YGyoto::Idx::isLast() {
+int YGyoto::Idx::isLast() const {
if (_is_range) return _cur+_range[2] > _range[1];
if (_is_scalar) return 1;
if (_is_list) return _cur >= _nel;
@@ -47,7 +50,7 @@ long YGyoto::Idx::first() {
return 0;
}
-int YGyoto::Idx::valid() {
+int YGyoto::Idx::valid() const {
if (_is_range) return _cur<=_range[1];
if (_is_scalar) return _cur==_range[0];
if (_is_list) return _cur<_nel;
@@ -64,6 +67,20 @@ long YGyoto::Idx::next() {
return 0;
}
+long YGyoto::Idx::range_min() const {
+ if (!(_is_range || _is_scalar)) throwError("BUG: not a range");
+ return _range[0];
+}
+
+long YGyoto::Idx::range_max() const {
+ if (!(_is_range || _is_scalar)) throwError("BUG: not a range");
+ return _range[1];
+}
+
+long YGyoto::Idx::range_dlt() const {
+ if (!(_is_range || _is_scalar)) throwError("BUG: not a range");
+ return _range[2];
+}
YGyoto::Idx::Idx(int iarg, int res) :
_is_nuller(0), _is_range(0), _is_list(0), _is_scalar(0), _is_double(0)
@@ -128,7 +145,7 @@ YGyoto::Idx::Idx(int iarg, int res) :
y_error("unsupported range syntax");
}
-int YGyoto::Idx::getNDims() {
+int YGyoto::Idx::getNDims() const {
if (_is_range) return 1;
if (_is_list) return 1;
if (_is_scalar) return 0;
diff --git a/yorick/gyoto_Scenery.C b/yorick/gyoto_Scenery.C
index 4d48b52..59c48a9 100644
--- a/yorick/gyoto_Scenery.C
+++ b/yorick/gyoto_Scenery.C
@@ -365,6 +365,7 @@ extern "C" {
double * data=ypush_d(dims);
Astrobj::Properties prop;
+ prop.alloc=Astrobj::Properties::minimal;
SmartPointer<Screen> screen = (*OBJ) -> screen();
# ifdef HAVE_UDUNITS
if (data) (*OBJ)->setPropertyConverters(&prop);
@@ -441,7 +442,16 @@ extern "C" {
ph->setInitCoord(coord, -1);
ph->delta((*OBJ)->delta());
ph->hit(&prop);
+ } else if (i_idx.isRangeOrScalar() && j_idx.isRangeOrScalar()){
+ // i and j can be trated as ranges
+ // We will use Scenery::rayTrace()
+ prop.di=i_idx.range_dlt();
+ prop.dj=j_idx.range_dlt();
+ (*OBJ)->rayTrace(i_idx.range_min(), i_idx.range_max(),
+ j_idx.range_min(), j_idx.range_max(),
+ &prop, impactcoords);
} else {
+ // i and j specify integers; at least one is a list.
screen -> computeBaseVectors();
double coord[8];
screen -> getRayCoord(size_t(i_idx.first()),
diff --git a/yorick/ygyoto_idx.h b/yorick/ygyoto_idx.h
index 7c1cb71..8e7d774 100644
--- a/yorick/ygyoto_idx.h
+++ b/yorick/ygyoto_idx.h
@@ -40,17 +40,21 @@ private:
public:
Idx(int iarg, int res);
- int isNuller();
- int getNDims();
- long getNElements();
+ int isNuller() const;
+ int getNDims() const;
+ long getNElements() const;
long first();
- int valid();
+ int valid() const;
long next();
- long current();
- double getDVal();
- int isDouble();
- int isFirst();
- int isLast();
+ long current() const;
+ double getDVal() const;
+ int isDouble() const;
+ int isRangeOrScalar() const;
+ int isFirst() const;
+ int isLast() const;
+ long range_min() const;
+ long range_max() const;
+ long range_dlt() const;
};
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-astro/packages/gyoto.git
More information about the Debian-astro-commits
mailing list