[Debian-astro-commits] [gyoto] 60/221: * New API: Screen::Coord1dSet & Coord2dSet * Change Scenery::rayTrace() API tu use a Coord2dSet instead of i/jmin/max * Make Astrobj::Properties::alloc a simple bool * Fix PATH in yorick/yorick1.in * Got rid of yorick-specific multi-thread code. * Support ray-tracing of a bunch of angle-specified directions, including in wYorick, and in parallel
Thibaut Jean-Claude Paumard
thibaut at moszumanska.debian.org
Fri May 22 20:52:33 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 fb035b8ec215a6e14bd87231645c70b9bd0729ed
Author: Thibaut Paumard <paumard at users.sourceforge.net>
Date: Thu Oct 16 16:15:01 2014 +0200
* New API: Screen::Coord1dSet & Coord2dSet
* Change Scenery::rayTrace() API tu use a Coord2dSet instead of i/jmin/max
* Make Astrobj::Properties::alloc a simple bool
* Fix PATH in yorick/yorick1.in
* Got rid of yorick-specific multi-thread code.
* Support ray-tracing of a bunch of angle-specified directions, including in wYorick, and in parallel
---
NEWS | 6 ++
bin/gyoto-mpi-worker.C | 7 +-
bin/gyoto.C | 20 ++--
include/GyotoAstrobj.h | 7 +-
include/GyotoScenery.h | 39 +++++--
include/GyotoScreen.h | 170 ++++++++++++++++++++++++++++-
lib/Astrobj.C | 6 +-
lib/Scenery.C | 196 ++++++++++++++++++----------------
lib/Screen.C | 151 ++++++++++++++++++++++++++
yorick/check-helpers.i | 2 +-
yorick/check-mpi.i | 23 ++++
yorick/check-scenery.i | 89 +++++++++++++++-
yorick/gyoto_Idx.C | 33 ++++--
yorick/gyoto_Scenery.C | 282 ++++++++++++++++---------------------------------
yorick/ygyoto_idx.h | 9 ++
yorick/yorick1.in | 1 +
16 files changed, 720 insertions(+), 321 deletions(-)
diff --git a/NEWS b/NEWS
index b047b24..c4977cd 100644
--- a/NEWS
+++ b/NEWS
@@ -1,3 +1,9 @@
+next future ABI 3
+ * New functionality: MPI-based parallelisation
+ * Change Scenery::rayTrace() API, replacing i/jmin and max with a
+ new concept Screen::Coord2dSet. This allows using the same
+ code-path for both the gyoto utility and all the yorick use cases.
+
0.2.1 2014/07/22 ABI 2.1
* ABI is backward compatible with 0.2.0.
* Support for the special values DBL_MAX, DBL_MIN, -DBL_MAX
diff --git a/bin/gyoto-mpi-worker.C b/bin/gyoto-mpi-worker.C
index 359fb0c..28cc6f4 100644
--- a/bin/gyoto-mpi-worker.C
+++ b/bin/gyoto-mpi-worker.C
@@ -23,6 +23,7 @@
#include "GyotoUtils.h"
#include "GyotoRegister.h"
#include "GyotoScenery.h"
+#include "GyotoScreen.h"
// FITS I/O
#include <fitsio.h>
@@ -67,7 +68,7 @@ int main(int argc, char** argv) {
*/
// For debug output
debug(0);
- // verbose(1);
+ verbose(GYOTO_SEVERE_VERBOSITY);
mpi::environment env(argc, argv);
@@ -88,6 +89,8 @@ int main(int argc, char** argv) {
sc = new Scenery();
sc -> mpi_team_ = &team;
+ Screen::Empty empty;
+
Scenery::mpi_tag task=Scenery::give_task;
Scenery::am_worker=true;
while (task != Scenery::terminate) {
@@ -105,7 +108,7 @@ int main(int argc, char** argv) {
case Scenery::raytrace:
curmsg = "In gyoto-mpi-worker.C: Error in ray-tracing: ";
curretval = 1;
- sc -> rayTrace(0, 0, 0, 0, NULL, NULL);
+ sc -> rayTrace(empty, NULL, NULL);
break;
case Scenery::terminate:
sc = NULL;
diff --git a/bin/gyoto.C b/bin/gyoto.C
index 0b79f90..65d5c7f 100644
--- a/bin/gyoto.C
+++ b/bin/gyoto.C
@@ -354,10 +354,7 @@ 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;
+ data->alloc=true;
size_t curquant=0;
size_t offset=res*res;
@@ -471,8 +468,19 @@ int main(int argc, char** argv) {
signal(SIGINT, sigint_handler);
curmsg = "In gyoto.C: Error during ray-tracing: ";
- scenery -> rayTrace(imin, imax, jmin, jmax, data,
- ipctdims[0]?impactcoords:NULL);
+
+ if (imax>res) imax=res;
+ if (jmax>res) jmax=res;
+
+
+ Screen::Range irange(imin, imax, di);
+ Screen::Range jrange(jmin, jmax, dj);
+ Screen::Grid grid(irange, jrange, "\rj = ");
+
+ if (verbose() >= GYOTO_QUIET_VERBOSITY)
+ cout << "j = " << imin << "/" << (imax-imin+1)/di << flush;
+
+ scenery -> rayTrace(grid, data, ipctdims[0]?impactcoords:NULL);
curmsg = "In gyoto.C: Error while saving: ";
if (verbose() >= GYOTO_QUIET_VERBOSITY)
diff --git a/include/GyotoAstrobj.h b/include/GyotoAstrobj.h
index c99aa2c..64408dd 100644
--- a/include/GyotoAstrobj.h
+++ b/include/GyotoAstrobj.h
@@ -738,11 +738,8 @@ 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;
+ /// True if buffers are allocated for entire field (npix*npix)
+ bool alloc;
public:
Properties(); ///< Default constructor (everything is set to NULL);
diff --git a/include/GyotoScenery.h b/include/GyotoScenery.h
index 9cd5014..338ad98 100644
--- a/include/GyotoScenery.h
+++ b/include/GyotoScenery.h
@@ -389,10 +389,9 @@ class Gyoto::Scenery : protected Gyoto::SmartPointee {
// Worker:
public:
- /// Perform ray-tracing for a square area on Screen
+ /// Perform ray-tracing
/**
- * For each Scenery::screen_ pixel in the square area limited by
- * imin, imax, jmin and jmax, launch a Photon back in time to
+ * For each directions specified, launch a Photon back in time to
* compute the various quantities.
*
* At this time, the computed quantities depend on on the pointers
@@ -407,7 +406,11 @@ class Gyoto::Scenery : protected Gyoto::SmartPointee {
* the various pointers in *data must be NULL or point to the first
* cell in an array of size at least Screen::npix_ squared.
*
- * If Scenery::nthreads_ is ≥2 and Gyoto has been compiled with
+ * If MPI support is built-in, MPI_Init() has been called, and
+ * nprocesses_ is ≥1, then rayTrace() will use several processes,
+ * launching them using mpiSpawn() if necessary.
+ *
+ * Else, if Scenery::nthreads_ is ≥2 and Gyoto has been compiled with
* pthreads support, rayTrace() will use Scenery::nthreads_ threads
* and launch photons in parallel. This works only if the
* Astrobj::Generic::clone() and Metric::Generic::clone() methods
@@ -415,9 +418,14 @@ class Gyoto::Scenery : protected Gyoto::SmartPointee {
* metric kind, and if they are both thread-safe. At the moment,
* unfortunately, Lorene metrics are known to not be thread-safe.
*
- * \param[in] imin, imax, jmin, jmax First and last rows and columns in
- * Scenery::screen_ to compute
-
+ * \param[in] ij Screen::Coord2dSet specification of rays to trace. e.g.:
+ *
+ * \code
+ * Screen::Range irange(imin, imax, di);
+ * Screen::Range jrange(jmin, jmax, dj);
+ * Screen::Grid ij(irange, jrange);
+ * \endcode
+ *
* \param[in, out] data Pointer to a preallocated
* Astrobj::Properties instance which sets which quantities must be
* computed and where to store the output.
@@ -431,9 +439,9 @@ class Gyoto::Scenery : protected Gyoto::SmartPointee {
* distributions are not. impactcoords can be computed using the
* ImpactCoords quantity.
*/
- void rayTrace(size_t imin, size_t imax, size_t jmin, size_t jmax,
- Astrobj::Properties* data, double * impactcoords = NULL);
-
+ void rayTrace(Screen::Coord2dSet & ij,
+ Astrobj::Properties *data,
+ double * impactcoords=NULL);
/// Ray-trace a single pixel in Scenery::screen_
/**
@@ -446,6 +454,17 @@ class Gyoto::Scenery : protected Gyoto::SmartPointee {
void operator() (size_t i, size_t j, Astrobj::Properties *data,
double * impactcoords = NULL, Photon * ph = NULL);
+ /// Ray-trace single direction
+ /**
+ * Almost identical to rayTrace(), but for a single direction.
+ *
+ * If ph is passed, it is assumed to have been properly initialized
+ * (with the right metric and astrobj etc.) already. Else, use
+ * &Scenery::ph_.
+ */
+ void operator() (double alpha, double delta, Astrobj::Properties *data,
+ Photon * ph = NULL);
+
#ifdef GYOTO_USE_XERCES
public:
/// Fill XML section
diff --git a/include/GyotoScreen.h b/include/GyotoScreen.h
index f40590f..fbb29a3 100644
--- a/include/GyotoScreen.h
+++ b/include/GyotoScreen.h
@@ -29,6 +29,13 @@
#include <iostream>
#include <fstream>
#include <string>
+#if defined HAVE_BOOST
+#include <boost/array.hpp>
+#define GYOTO_ARRAY boost::array
+#else
+#include <array>
+#define GYOTO_ARRAY std::array
+#endif
namespace Gyoto {
class Screen;
@@ -555,7 +562,168 @@ class Gyoto::Screen : protected Gyoto::SmartPointee {
static SmartPointer<Screen> Subcontractor(FactoryMessenger* fmp);
#endif
-
+ /// Enum to specify whether a coordinate set (Coord1dSet or Coord2dSet) holds pixel values or angles
+ enum CoordType_e {angle, pixel};
+
+ /// Set of 1-d coordinates: indices or angles
+ /**
+ * Acts like a container (array-like) of either size_t (pixel
+ * coordinate) or double (angle) values. This container can be
+ * iterated-through using the operator++(), derefenced using the
+ * operator*() (if containing pixel coordinates) or angle() (in
+ * containing angles).
+ */
+ class Coord1dSet {
+ public:
+ /// Whether this specifier represents angles or pixels
+ const CoordType_e kind;
+ public:
+ /// Set kind during initialization
+ Coord1dSet(CoordType_e k);
+ /// Reset specifier to point to the first value
+ virtual void begin() =0;
+ /// True if pointing to something, false if end has been reached.
+ virtual bool valid() =0;
+ /// Number of values in this container
+ virtual size_t size()=0;
+ /// Get size_t value crrently pointed to
+ virtual size_t operator*() const ;
+ /// Get double value currently pointed to
+ virtual double angle() const ;
+ /// Increment iterator (point to next value)
+ virtual Coord1dSet& operator++()=0;
+ };
+
+ /// Class to specify a set of points on the Screen
+ /**
+ * Container (array-like) holding several 2D points. Can be a 2D
+ * grid of pixel coordinates or a vector of floating-point (alpha,
+ * delta) pairs, for instance.
+ */
+ class Coord2dSet {
+ public:
+ /// Whether this set holds pixels or angle specifications
+ const CoordType_e kind;
+ /// Set kind at initialisation
+ Coord2dSet(CoordType_e k);
+ /// Increment pointer
+ virtual Coord2dSet& operator++() =0;
+ /// Get pixel coordinates
+ virtual GYOTO_ARRAY<size_t, 2> operator* () const;
+ /// Get angle coordinates
+ virtual GYOTO_ARRAY<double, 2> angles() const ;
+ /// Reset pointer
+ virtual void begin() =0;
+ /// Whether the end has not been passed
+ virtual bool valid() =0;
+ /// Number of positions contained
+ virtual size_t size()=0;
+ };
+
+ /// Class containing 2D-points organized in a grid
+ class Grid: public Coord2dSet {
+ protected:
+ protected:
+ /// If non-NULL, cout j each tims it is incremented.
+ const char * const prefix_;
+ Coord1dSet &iset_;
+ Coord1dSet &jset_;
+ public:
+ Grid(Coord1dSet &iset, Coord1dSet &jset, const char * const p=NULL);
+ virtual Coord2dSet& operator++();
+ virtual GYOTO_ARRAY<size_t, 2> operator* () const;
+ virtual void begin();
+ virtual bool valid();
+ virtual size_t size();
+ };
+
+ /// Class containing arbitrary 2D-points
+ /**
+ * ispec_ and jspec_ must be the same size.
+ */
+ class Bucket : public Coord2dSet {
+ protected:
+ Coord1dSet &alpha_;
+ Coord1dSet &delta_;
+ public:
+ Bucket(Coord1dSet &iset, Coord1dSet &jset);
+ virtual Coord2dSet& operator++();
+ virtual GYOTO_ARRAY<double, 2> angles() const;
+ virtual GYOTO_ARRAY<size_t, 2> operator*() const;
+ virtual void begin();
+ virtual bool valid();
+ virtual size_t size();
+ };
+
+ /// A dummy, empty 2D set.
+ class Empty: public Coord2dSet {
+ public:
+ Empty();
+ virtual Coord2dSet& operator++();
+ virtual void begin();
+ virtual bool valid();
+ virtual size_t size();
+ };
+
+ /// 1D coordinated specifier for a range
+ class Range : public Coord1dSet {
+ protected:
+ const size_t mi_, ma_, d_, sz_;
+ size_t cur_;
+ public:
+ /// Specify min, max and step of this range.
+ Range(size_t mi, size_t ma, size_t d);
+ void begin();
+ bool valid();
+ size_t size();
+ Coord1dSet& operator++();
+ size_t operator*() const ;
+ };
+
+ /// 1D specifier for an arbitrary pixel coordinate set.
+ class Indices : public Coord1dSet {
+ protected:
+ size_t const * const indices_;
+ size_t const sz_;
+ size_t i_;
+ public:
+ Indices (size_t const*const buf, size_t sz);
+ void begin();
+ bool valid();
+ size_t size();
+ Coord1dSet& operator++();
+ size_t operator*() const ;
+ };
+
+ /// 1D specifier for an arbitrary angle coordinate set.
+ class Angles : public Coord1dSet {
+ protected:
+ double const * const buf_;
+ size_t const sz_;
+ size_t i_;
+ public:
+ Angles (double const*const buf, size_t sz);
+ void begin();
+ bool valid();
+ size_t size();
+ Coord1dSet& operator++();
+ double angle() const ;
+ };
+
+ /// 1D specifier for an angle that is repeated.
+ class RepeatAngle : public Coord1dSet {
+ protected:
+ double const val_;
+ size_t const sz_;
+ size_t i_;
+ public:
+ RepeatAngle (double val, size_t sz);
+ void begin();
+ bool valid();
+ size_t size();
+ Coord1dSet& operator++();
+ double angle() const ;
+ };
};
#endif
diff --git a/lib/Astrobj.C b/lib/Astrobj.C
index 0a3da93..fe23119 100644
--- a/lib/Astrobj.C
+++ b/lib/Astrobj.C
@@ -487,8 +487,7 @@ Astrobj::Properties::Properties() :
, intensity_converter_(NULL), spectrum_converter_(NULL),
binspectrum_converter_(NULL)
# endif
- , alloc(Properties::entire_field | Properties::all_rows)
- , di(1), dj(1)
+ , alloc(false)
{}
Astrobj::Properties::Properties( double * I, double * t) :
@@ -501,8 +500,7 @@ 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)
+ , alloc(false)
{}
void Astrobj::Properties::init(size_t nbnuobs) {
diff --git a/lib/Scenery.C b/lib/Scenery.C
index 0f800da..3ff630f 100644
--- a/lib/Scenery.C
+++ b/lib/Scenery.C
@@ -35,6 +35,7 @@
#include <boost/mpi/collectives.hpp>
#include <string>
#include <boost/serialization/string.hpp>
+#include <boost/serialization/array.hpp>
namespace mpi = boost::mpi;
#endif
@@ -135,13 +136,23 @@ typedef struct SceneryThreadWorkerArg {
pthread_mutex_t * mutex;
pthread_t * parent;
#endif
- size_t i, j, imin, imax, jmin, jmax, eol_offset;
+ Screen::Coord2dSet & ij;
+ size_t cnt, npix;
Scenery *sc;
Photon * ph;
Astrobj::Properties *data;
double * impactcoords;
+ SceneryThreadWorkerArg(Screen::Coord2dSet & ijin);
+ bool is_pixel;
} SceneryThreadWorkerArg ;
+SceneryThreadWorkerArg::SceneryThreadWorkerArg(Screen::Coord2dSet & ijin)
+ :ij(ijin)
+{
+
+}
+
+
static void * SceneryThreadWorker (void *arg) {
/*
This is the real ray-tracing loop. It may be called by multiple
@@ -162,11 +173,8 @@ static void * SceneryThreadWorker (void *arg) {
#endif
// local variables to store our parameters
- size_t i, j;
-
- // 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;
+ GYOTO_ARRAY<size_t, 2> ijb;
+ GYOTO_ARRAY<double, 2> ad;
Astrobj::Properties data;
double * impactcoords = NULL;
@@ -182,9 +190,13 @@ static void * SceneryThreadWorker (void *arg) {
// lock mutex so we can safely read and update i, j et al.
if (larg->mutex) pthread_mutex_lock(larg->mutex);
#endif
- // copy i & j
- i = larg->i; j = larg->j;
- if (j > larg->jmax || (j==larg->jmax && i>larg->imax)) {
+
+ // copy i & j or alpha and delta
+
+ if (larg->is_pixel) ijb = *(larg->ij);
+ else ad = larg->ij.angles();
+
+ if (!larg->ij.valid()) {
// terminate, but first...
#ifdef HAVE_PTHREAD
// ...unlock mutex so our siblings can access i & j and terminate too
@@ -193,38 +205,26 @@ static void * SceneryThreadWorker (void *arg) {
break;
}
- data = *larg->data;
- // update i, j and pointer
- larg->i += data.di;
- *larg->data += data.alloc&Astrobj::Properties::all_rows?data.di:1;
- if (larg->impactcoords) {
- impactcoords = larg->impactcoords; larg->impactcoords+=16*data.di;
- }
- if (larg->i > larg->imax) {
- larg->j+=data.dj; larg->i=larg->imin;
- (*larg->data) += eol_offset;
- if (larg->impactcoords) larg->impactcoords+=16*eol_offset;
- }
+ size_t lcnt = larg->cnt++;
+
+ ++(larg->ij);
#ifdef HAVE_PTHREAD
// unlock mutex so our siblings can can access i, j et al. and procede
if (larg->mutex) pthread_mutex_unlock(larg->mutex);
#endif
- ////// 2- do the actual work.
- if (i==larg->imin && verbose() >= GYOTO_QUIET_VERBOSITY && !impactcoords) {
-# ifdef HAVE_PTHREAD
- if (larg->mutex) pthread_mutex_lock(larg->mutex);
-# endif
- cout << "\rj = " << j << " / " << larg->jmax << " " << flush;
-# ifdef HAVE_PTHREAD
- if (larg->mutex) pthread_mutex_unlock(larg->mutex);
-# endif
- }
-# if GYOTO_DEBUG_ENABLED
- GYOTO_DEBUG << "i = " << i << ", j = " << j << endl;
-# endif
- (*larg->sc)(i, j, &data, impactcoords, ph);
+ // Store current cell
+ data = *larg->data;
+ size_t cell=lcnt;
+ if (larg->is_pixel && data.alloc) cell=(ijb[1]-1)*larg->npix+ijb[0]-1;
+ data += cell;
+ impactcoords=larg->impactcoords?larg->impactcoords+16*cell:NULL;
+
+ if (larg->is_pixel)
+ (*larg->sc)(ijb[0], ijb[1], &data, impactcoords, ph);
+ else (*larg->sc)(ad[0], ad[1], &data, ph);
+
++count;
}
#ifdef HAVE_PTHREAD
@@ -238,10 +238,10 @@ static void * SceneryThreadWorker (void *arg) {
return NULL;
}
-void Scenery::rayTrace(size_t imin, size_t imax,
- size_t jmin, size_t jmax,
+void Scenery::rayTrace(Screen::Coord2dSet & ij,
Astrobj::Properties *data,
double * impactcoords) {
+ GYOTO_DEBUG_EXPR(am_worker);
#if defined HAVE_MPI
if (nprocesses_ && !mpi_team_) {
@@ -259,10 +259,10 @@ void Scenery::rayTrace(size_t imin, size_t imax,
- some housekeeping
*/
- const size_t npix = screen_->resolution();
- imax=(imax<=(npix)?imax:(npix));
- jmax=(jmax<=(npix)?jmax:(npix));
-
+ if (!screen_) {
+ if (am_worker) throwError("No screen, have you called mpiClone()?");
+ else throwError("Scenery::rayTrace() needs a Screen to work on");
+ }
screen_->computeBaseVectors();
// Necessary for KS integration, computes relation between
// observer's x,y,z coord and KS X,Y,Z coord. Will be used to
@@ -274,27 +274,17 @@ void Scenery::rayTrace(size_t imin, size_t imax,
ph_.spectrometer(spr);
ph_.freqObs(screen_->freqObs());
double coord[8];
- screen_ -> getRayCoord(imin,jmin, coord);
+ screen_ -> getRayCoord(size_t(1),size_t(1), coord);
ph_ . setInitCoord(coord, -1);
// delta is reset in operator()
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);
- }
+ GYOTO_ARRAY<size_t, 2> ijb;
+ GYOTO_ARRAY<double, 2> ad;
+
+ bool alloc=data?data->alloc:false;
+ size_t npix=screen_->resolution();
#ifdef HAVE_MPI
if (mpi_team_) {
@@ -306,14 +296,14 @@ void Scenery::rayTrace(size_t imin, size_t imax,
mpiTask(tag);
}
- size_t ij[2]={imin, jmin};
-
size_t nbnuobs=0;
Quantity_t quantities = (am_worker || !data)?GYOTO_QUANTITY_NONE:*data;
bool has_ipct=am_worker?false:bool(impactcoords);
+ bool is_pixel=(ij.kind==Screen::pixel);
mpi::broadcast(*mpi_team_, quantities, 0);
mpi::broadcast(*mpi_team_, has_ipct, 0);
+ mpi::broadcast(*mpi_team_, is_pixel, 0);
if (quantities & (GYOTO_QUANTITY_SPECTRUM | GYOTO_QUANTITY_BINSPECTRUM)) {
if (!spr) throwError("Spectral quantity requested but "
@@ -387,6 +377,7 @@ void Scenery::rayTrace(size_t imin, size_t imax,
// The corresponding recv is in gyoto-scenery-worker.c
vector<size_t> cell(working+1);
+ size_t cnt=0;
while (working) {
// receive one result, need to track back where it belongs and
@@ -403,30 +394,32 @@ void Scenery::rayTrace(size_t imin, size_t imax,
size_t cs=cell[w]; // remember where to store results
// give new task or decrease working counter
- if (ij[0]<=imax) {
- cell[w]=curcell; // store curcell
- mpi_team_ -> send(w, raytrace, ij, 2);
+ if (ij.valid()) {
+ cell[w]=cnt++;
+ if (is_pixel) {
+ ijb=*ij;
+ if (alloc) cell[w]=(ijb[1]-1)*npix+ijb[0]-1;
+ mpi_team_ -> send(w, raytrace, ijb);
+ } else {
+ ad = ij.angles();
+ mpi_team_ -> send(w, raytrace, ad);
+ }
if (impactcoords) {
mpi_team_ -> send(w, Scenery::impactcoords, impactcoords+cell[w]*16, 16);
}
- curcell += (data&&data->alloc&Astrobj::Properties::all_rows)?
- data->di:1;
+ ++ij;
- 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);
+ if (is_pixel)
+ mpi_team_ -> send(w, raytrace_done, ijb);
+ else
+ mpi_team_ -> send(w, raytrace_done, ad);
--working;
}
// Now that the worker is back to work, triage data it has just delivered
+
if (s.tag()==Scenery::raytrace_done && data) {
// Copy each relevant quantity, performing conversion if needed
if (data->intensity)
@@ -471,7 +464,10 @@ void Scenery::rayTrace(size_t imin, size_t imax,
mpi_team_->send(0, give_task, vect, nelt);
while (true) {
// Receive new coordinates to work on.
- s = mpi_team_->recv(0, mpi::any_tag, ij, 2);
+ if (is_pixel)
+ s = mpi_team_->recv(0, mpi::any_tag, ijb);
+ else
+ s = mpi_team_->recv(0, mpi::any_tag, ad);
if (s.tag()==raytrace_done) {
break;
}
@@ -479,7 +475,10 @@ void Scenery::rayTrace(size_t imin, size_t imax,
if (has_ipct)
s = mpi_team_->recv(0, Scenery::impactcoords, impactcoords, 16);
locdata->init(nbnuobs);
- (*this)(ij[0], ij[1], locdata, impactcoords, &ph_);
+ if (is_pixel)
+ (*this)(ijb[0], ijb[1], locdata, impactcoords, &ph_);
+ else
+ (*this)(ad[0], ad[1], locdata, &ph_);
// send result
mpi_team_->send(0, raytrace_done, vect, nelt);
}
@@ -490,24 +489,14 @@ void Scenery::rayTrace(size_t imin, size_t imax,
}
#endif
- if (data) {
- size_t first_index=(jmin-1)*npix + imin -1;
- (*data) += first_index;
- if (impactcoords) impactcoords += first_index * 16;
- }
-
- SceneryThreadWorkerArg larg;
+ SceneryThreadWorkerArg larg(ij);
larg.sc=this;
larg.ph=&ph_;
larg.data=data;
+ larg.cnt=0;
+ larg.npix=npix;
larg.impactcoords=impactcoords;
- larg.i=imin;
- larg.j=jmin;
- larg.imin=imin;
- larg.imax=imax;
- larg.jmin=jmin;
- larg.jmax=jmax;
- larg.eol_offset=eol_offset;
+ larg.is_pixel= (ij.kind==Screen::pixel);
struct timeval tim;
double start, end;
@@ -545,7 +534,7 @@ void Scenery::rayTrace(size_t imin, size_t imax,
gettimeofday(&tim, NULL);
end=double(tim.tv_sec)+(double(tim.tv_usec)/1000000.0);
- GYOTO_MSG << "\nRaytraced "<< (jmax-jmin+1) * (imax-imin+1)
+ GYOTO_MSG << "\nRaytraced "<< ij.size()
<< " photons in " << end-start
<< "s using " << nthreads_ << " thread(s)\n";
@@ -556,6 +545,7 @@ void Scenery::operator() (
Astrobj::Properties *data, double * impactcoords,
Photon *ph
) {
+
double coord[8];
SmartPointer<Spectrometer::Generic> spr = screen_->spectrometer();
size_t nbnuobs = spr() ? spr -> nSamples() : 0;
@@ -601,6 +591,32 @@ void Scenery::operator() (
}
}
+void Scenery::operator() (
+ double a, double d,
+ Astrobj::Properties *data,
+ Photon *ph
+ ) {
+
+ double coord[8];
+ SmartPointer<Spectrometer::Generic> spr = screen_->spectrometer();
+ size_t nbnuobs = spr() ? spr -> nSamples() : 0;
+
+ if (data) data -> init(nbnuobs);
+
+ if (!ph) {
+ ph = &ph_;
+ ph -> spectrometer(spr);
+ ph -> freqObs(screen_->freqObs());
+ }
+ // Always reset delta
+ ph -> delta(delta_);
+
+ screen_ -> getRayCoord(a, d, coord);
+ ph -> setInitCoord(coord, -1);
+ ph -> hit(data);
+
+}
+
SmartPointer<Photon> Scenery::clonePhoton() const {
return ph_.clone();
}
diff --git a/lib/Screen.C b/lib/Screen.C
index 77d68c3..7bf3ccb 100644
--- a/lib/Screen.C
+++ b/lib/Screen.C
@@ -1209,3 +1209,154 @@ SmartPointer<Screen> Screen::Subcontractor(FactoryMessenger* fmp) {
return scr;
}
#endif
+
+//////// COORD2DSET
+
+Screen::Coord2dSet::Coord2dSet(CoordType_e k) : kind(k) {}
+
+GYOTO_ARRAY<size_t, 2> Screen::Coord2dSet::operator* () const {
+ if (kind==pixel)
+ throwError("BUG: Coord2dSet of kind pixel should implement operator*");
+ else
+ throwError("Coord2dSet of kind angle should not be dereferenced");
+ // avoid warning
+ GYOTO_ARRAY<size_t, 2> a = {0, 0};
+ return a;
+}
+
+GYOTO_ARRAY<double, 2> Screen::Coord2dSet::angles () const {
+ if (kind==Screen::angle)
+ throwError("BUG: Coord2dSet of kind angle should implement angles()");
+ else
+ throwError("angles() should not be called on Coord2dSet of kind pixel");
+ // avoid warning
+ GYOTO_ARRAY<double, 2> a = {0., 0.};
+ return a;
+}
+
+//////// GRID
+
+Screen::Grid::Grid(Coord1dSet &iset, Coord1dSet &jset,
+ const char * const p)
+ : Coord2dSet(pixel),
+ prefix_(p),
+ iset_(iset), jset_(jset)
+{}
+
+GYOTO_ARRAY<size_t, 2> Screen::Grid::operator* () const {
+ GYOTO_ARRAY<size_t, 2> ij = {*iset_, *jset_};
+ return ij;
+}
+void Screen::Grid::begin() {iset_.begin(); jset_.begin();}
+bool Screen::Grid::valid() {return iset_.valid() && jset_.valid();}
+size_t Screen::Grid::size() {return (iset_.size())*(jset_.size());}
+
+Screen::Coord2dSet& Screen::Grid::operator++() {
+ if (!valid()) return *this;
+ ++iset_;
+ if (!iset_.valid()) {
+ iset_.begin();
+ ++jset_;
+ if (prefix_ && verbose() >= GYOTO_QUIET_VERBOSITY && jset_.valid())
+ cout << prefix_ << *jset_ << "/" << jset_.size() << flush;
+ }
+
+ return *this;
+}
+
+//////// COORD1DSET
+
+Screen::Coord1dSet::Coord1dSet(CoordType_e k) : kind(k) {}
+
+size_t Screen::Coord1dSet::operator* () const {
+ if (kind==pixel)
+ throwError("BUG: Coord1dSet of kind pixel should implement operator*");
+ else
+ throwError("Coord1dSet of kind angle should not be dereferenced");
+ // avoid warning
+ return 0;
+}
+
+double Screen::Coord1dSet::angle () const {
+ if (kind==Screen::angle)
+ throwError("BUG: Coord1dSet of kind angle should implement angle()");
+ else
+ throwError("angle() should not be called on Coord1dSet of kind pixel");
+ // avoid warning
+ return 0.;
+}
+
+///////
+
+Screen::Range::Range(size_t mi, size_t ma, size_t d)
+ : Coord1dSet(pixel), mi_(mi), ma_(ma), d_(d), cur_(mi), sz_((ma-mi+1)/d)
+{}
+
+void Screen::Range::begin() {cur_=mi_;}
+Screen::Coord1dSet& Screen::Range::operator++() {
+ cur_ += d_; return *this;
+}
+bool Screen::Range::valid() {return cur_ <= ma_;}
+size_t Screen::Range::size() {return sz_;}
+size_t Screen::Range::operator*() const {return cur_;}
+
+//////
+
+Screen::Indices::Indices (size_t const*const buf, size_t sz)
+ : Coord1dSet(pixel), indices_(buf), sz_(sz), i_(0)
+{}
+void Screen::Indices::begin() {i_=0;}
+bool Screen::Indices::valid() {return i_ < sz_;}
+size_t Screen::Indices::size(){return sz_;}
+Screen::Coord1dSet& Screen::Indices::operator++() {++i_; return *this;}
+size_t Screen::Indices::operator*() const {return indices_[i_];}
+
+
+/////
+
+Screen::Angles::Angles (double const*const buf, size_t sz)
+ : Coord1dSet(Screen::angle), buf_(buf), sz_(sz), i_(0)
+{}
+void Screen::Angles::begin() {i_=0;}
+bool Screen::Angles::valid() {return i_<sz_;}
+size_t Screen::Angles::size(){return sz_;}
+Screen::Coord1dSet& Screen::Angles::operator++(){++i_;}
+double Screen::Angles::angle() const {return buf_[i_];}
+
+Screen::RepeatAngle::RepeatAngle (double val, size_t sz)
+ : Coord1dSet(Screen::angle), val_(val), sz_(sz), i_(0)
+{}
+void Screen::RepeatAngle::begin() {i_=0;}
+bool Screen::RepeatAngle::valid() {return i_<sz_;}
+size_t Screen::RepeatAngle::size(){return sz_;}
+Screen::Coord1dSet& Screen::RepeatAngle::operator++(){++i_;}
+double Screen::RepeatAngle::angle() const {return val_;}
+
+Screen::Bucket::Bucket (Coord1dSet &alp, Coord1dSet &del)
+ : Coord2dSet(alp.kind), alpha_(alp), delta_(del)
+{
+ if (alp.kind != angle || alp.kind != angle)
+ throwError("both specifiers must be of same kind");
+ if (alp.size() != del.size())
+ throwError("alpha and delta should be of same size");
+}
+void Screen::Bucket::begin() {alpha_.begin(); delta_.begin();}
+bool Screen::Bucket::valid() {return alpha_.valid() && delta_.valid();}
+size_t Screen::Bucket::size(){return alpha_.size();}
+Screen::Coord2dSet& Screen::Bucket::operator++(){++alpha_; ++delta_;}
+GYOTO_ARRAY<double, 2> Screen::Bucket::angles() const {
+ GYOTO_ARRAY<double, 2> out {alpha_.angle(), delta_.angle()};
+ return out;
+}
+GYOTO_ARRAY<size_t, 2> Screen::Bucket::operator* () const {
+ GYOTO_ARRAY<size_t, 2> ij = {*alpha_, *delta_};
+ return ij;
+}
+
+///// Empty
+
+Screen::Empty::Empty () : Coord2dSet(pixel) {}
+Screen::Coord2dSet& Screen::Empty::operator++() {return *this;}
+void Screen::Empty::begin() {}
+bool Screen::Empty::valid() {return false;}
+size_t Screen::Empty::size() {return 0;}
diff --git a/yorick/check-helpers.i b/yorick/check-helpers.i
index 8335aba..8bbf0ec 100644
--- a/yorick/check-helpers.i
+++ b/yorick/check-helpers.i
@@ -158,7 +158,7 @@ func check_christoffels(gg, pos, tolerance=, eps=)
Gamma3=christoffel(gg, pos(,n), eps=eps);
if (anyof(Gamma1!=Gamma2))
error, "The two forms of the christoffel method don't yield the same result";
- if (max(abs(Gamma1(ind)-Gamma3(ind)))>tolerance)
+ if (max(abs(Gamma1-Gamma3))>tolerance)
error, "The Christoffels don't agree with their numerical estimate";
dot;
}
diff --git a/yorick/check-mpi.i b/yorick/check-mpi.i
index ee3226b..328228d 100644
--- a/yorick/check-mpi.i
+++ b/yorick/check-mpi.i
@@ -73,6 +73,29 @@ mdiff=max(abs(diff));
if (mdiff > 1e-6) error, "Results differ";
output, " OK (max rel. dif.: "+pr1(mdiff)+")";
+fov=sc.screen().fov();
+npix=sc.screen().resolution();
+delta= fov/double(npix);
+
+xx=((indgen(npix)-0.5)*delta-fov/2.)(, -:1:32);
+yy=transpose(xx);
+
+sc, mpispawn=8, mpiclone=;
+
+doing, "Integrating whole field, specifying angles, with MPI";
+data2 = sc(-xx, yy, );
+done;
+
+sc, mpispawn=0;
+
+doing, "Comparing results";
+diff=data-data2;
+ind=where(data);
+diff(ind)/=data(ind);
+mdiff=max(abs(diff));
+if (mdiff > 1e-6) error, "Results differ";
+output, " OK (max rel. dif.: "+pr1(mdiff)+")";
+
doing, "Deleting Scenery";
sc=[];
done;
diff --git a/yorick/check-scenery.i b/yorick/check-scenery.i
index ee26a7e..5432e96 100644
--- a/yorick/check-scenery.i
+++ b/yorick/check-scenery.i
@@ -18,6 +18,7 @@
*/
#include "check-helpers.i"
+restore, gyoto;
begin_section, "Scenery";
@@ -211,11 +212,97 @@ sc4;
ph = gyoto_Photon(initcoord=sc3, 6, 19);
ph.is_hit();
+doing, "Reading Scenery...";
+sc=Scenery("../doc/examples/example-complex-astrobj.xml");
+done;
+
+sc, nthreads=8, nprocesses=0, mpispawn=0;
+
+doing, "Integrating whole field...";
+data=sc();
+done;
+
+r1=8:25:4;
+r2=2:-2:3;
+v1=[1, 4, 16];
+v2=[15, 20, 22];
+
+doing, "Integrating subfield...";
+data2=sc(r1, r2, );
+done;
+
+doing, "Comparing...";
+if (anyof(data2 != data(r1, r2, ))) error, "result differ";
+done;
+
+doing, "Integrating subfield...";
+data2=sc(v1, v2, );
+done;
+
+doing, "Comparing...";
+if (anyof(data2 != data(v1, v2, ))) error, "result differ";
+done;
+
+doing, "Integrating subfield...";
+data2=sc(r1, v2, );
+done;
+
+doing, "Comparing...";
+if (anyof(data2 != data(r1, v2, ))) error, "result differ";
+done;
+
+doing, "Integrating subfield...";
+data2=sc(v1, r2, );
+done;
+
+doing, "Comparing...";
+if (anyof(data2 != data(v1, r2, ))) error, "result differ";
+done;
+
+fov=sc.screen().fov();
+npix=sc.screen().resolution();
+delta= fov/double(npix);
+
+xx=((indgen(npix)-0.5)*delta-fov/2.)(, -:1:32);
+yy=transpose(xx);
+
+data2=array(double, npix, npix, 1);
+doing, "integrating pix. by pix., specifying angles...\n";
+verbosity=verbose(0);
+for (j=1; j<=npix; ++j) {
+ write, format="\rj = %d/%d", j, npix;
+ for (i=1; i<=npix; ++i) {
+ data2(i, j, 1) = sc(-xx(i, j), yy(i, j), );
+ }
+ }
+write, format="%s\n", "";
+done;
+verbose, verbosity;
+
+doing, "Comparing results";
+diff=data-data2;
+ind=where(data);
+diff(ind)/=data(ind);
+mdiff=max(abs(diff));
+if (mdiff > 1e-6) error, "Results differ";
+output, " OK (max rel. dif.: "+pr1(mdiff)+")";
+
+doing, "integrating whole field, specifying angles...\n";
+data2 = sc(-xx, yy, );
+done;
+
+doing, "Comparing results";
+diff=data-data2;
+ind=where(data);
+diff(ind)/=data(ind);
+mdiff=max(abs(diff));
+if (mdiff > 1e-6) error, "Results differ";
+output, " OK (max rel. dif.: "+pr1(mdiff)+")";
if (batch()) {
// Free memory for easier checking with valgrind
- data=[];
+ xx=yy=data2=data=ind=[];
sc4=[];
sc3=[];
sc2=[];
diff --git a/yorick/gyoto_Idx.C b/yorick/gyoto_Idx.C
index cbbea97..8e7f34c 100644
--- a/yorick/gyoto_Idx.C
+++ b/yorick/gyoto_Idx.C
@@ -20,10 +20,12 @@
#include "ygyoto_idx.h"
#include "yapi.h"
#include "GyotoError.h"
-
+#include "GyotoDefs.h"
+#include "GyotoUtils.h"
+#include <iostream>
using namespace YGyoto;
using namespace Gyoto;
-
+using namespace std;
int YGyoto::Idx::isNuller() const {return _is_nuller;}
long YGyoto::Idx::getNElements() const {return _nel;}
long YGyoto::Idx::current() const {
@@ -31,7 +33,7 @@ long YGyoto::Idx::current() const {
return _cur;
}
double YGyoto::Idx::getDVal() const {return _is_double?_dval:_range[0];}
-int YGyoto::Idx::isDouble() const {return _is_double;}
+int YGyoto::Idx::isDouble() const {return _is_double || _is_dlist;}
int YGyoto::Idx::isRangeOrScalar() const {return _is_range || _is_scalar;}
int YGyoto::Idx::isFirst() const {return _is_first;}
@@ -83,7 +85,8 @@ long YGyoto::Idx::range_dlt() const {
}
YGyoto::Idx::Idx(int iarg, int res) :
- _is_nuller(0), _is_range(0), _is_list(0), _is_scalar(0), _is_double(0)
+ _is_nuller(0), _is_range(0), _is_list(0), _is_scalar(0), _is_double(0),
+ _buf(NULL), _is_dlist(0)
{
int flags = yget_range(iarg, _range), ndims=0;
if (flags) {
@@ -113,10 +116,11 @@ YGyoto::Idx::Idx(int iarg, int res) :
return;
}
if (yarg_rank(iarg) > 0) {
- _is_list=1;
- _nel=1;
- _idx = ygeta_l(iarg, &_nel, 0);
- return;
+ if (yarg_number(iarg)==1){
+ _is_list=1;
+ _idx = ygeta_l(iarg, &_nel, _dims);
+ return;
+ }
}
if (yarg_number(iarg)==1) {
_is_scalar=1;
@@ -131,7 +135,13 @@ YGyoto::Idx::Idx(int iarg, int res) :
if (yarg_number(iarg)==2) {
_is_scalar=1;
_is_double=1;
- _dval=ygets_d(iarg);
+ _buf=ygeta_d(iarg, &_nel, _dims);
+ _dval=_buf[0];
+ if (_dims[0]) _is_dlist=1;
+ else _is_scalar=1;
+ GYOTO_DEBUG_ARRAY(_dims, Y_DIMSIZE);
+ GYOTO_DEBUG_EXPR(_is_scalar);
+ GYOTO_DEBUG_EXPR(_is_dlist);
return;
}
if (iarg<0 || yarg_nil(iarg)) {
@@ -150,3 +160,8 @@ int YGyoto::Idx::getNDims() const {
if (_is_list) return 1;
if (_is_scalar) return 0;
}
+
+long const * YGyoto::Idx::getDims() const {return _dims;}
+
+long const * YGyoto::Idx::getBuffer() const {return _idx;}
+double const * YGyoto::Idx::getDoubleBuffer() const {return _buf;}
diff --git a/yorick/gyoto_Scenery.C b/yorick/gyoto_Scenery.C
index 72b2698..55a1489 100644
--- a/yorick/gyoto_Scenery.C
+++ b/yorick/gyoto_Scenery.C
@@ -41,111 +41,6 @@ using namespace std;
using namespace Gyoto;
using namespace YGyoto;
-typedef struct ySceneryThreadWorkerArg {
-#ifdef HAVE_PTHREAD
- pthread_mutex_t * mutex;
- pthread_t * parent;
-#endif
- Idx *i_idx, *j_idx;
- Scenery *sc;
- Photon * ph;
- Astrobj::Properties *data;
- double * impactcoords;
- size_t res;
-} ySceneryThreadWorkerArg ;
-
-static void * ySceneryThreadWorker (void *arg) {
- /*
- This is the real ray-tracing loop. It may be called by multiple
- threads in parallel, launched from ::rayTrace
- */
-
- ySceneryThreadWorkerArg *larg = static_cast<ySceneryThreadWorkerArg*>(arg);
-
- // Each thread needs its own Photon, clone cached Photon
- // it is assumed to be already initialized with spectrometer et al.
- Photon * ph = larg -> ph;
-#ifdef HAVE_PTHREAD
- if (larg->mutex) {
- GYOTO_DEBUG << "locking mutex\n";
- pthread_mutex_lock(larg->mutex);
- GYOTO_DEBUG << "mutex locked\n";
- ph = larg -> ph -> clone();
- GYOTO_DEBUG << "unlocking mutex\n";
- pthread_mutex_unlock(larg->mutex);
- GYOTO_DEBUG << "mutex unlocked\n";
- }
-#endif
-
- // local variables to store our parameters
- size_t i, j;
- Astrobj::Properties data;
- double * impactcoords = NULL;
-
- size_t count=0;
-
- while (1) {
- /////// 1- get input and output parameters and update them for next access
- //// i and j are input, data and impactcoords are where to store
- //// output. we must get them and increase them so that another
- //// thread can get the next values while we integrate.
-#ifdef HAVE_PTHREAD
- // lock mutex so we can safely read and update i, j et al.
- if (larg->mutex) pthread_mutex_lock(larg->mutex);
-#endif
- // copy i & j
- i = larg->i_idx->current(); j = larg->j_idx->current();
-
- // check whether there remains something to compute
- if (!larg->j_idx->valid()
- || (larg->j_idx->isLast() && !larg->i_idx->valid())) {
- // terminate, but first...
-#ifdef HAVE_PTHREAD
- // ...unlock mutex so our siblings can access i & j and terminate too
- if (larg->mutex) pthread_mutex_unlock(larg->mutex);
-#endif
- break;
- }
-
- // print some info
- if (larg->i_idx->isFirst() &&
- verbose() >= GYOTO_QUIET_VERBOSITY && !impactcoords) {
- cout << "\rRay-tracing scenery: j = " << j << flush ;
- }
-# if GYOTO_DEBUG_ENABLED
- GYOTO_DEBUG << "i = " << i << ", j = " << j << endl;
-# endif
-
- // update i & j
- larg->i_idx->next();
- if (!larg->i_idx->valid()) {
- larg->j_idx->next(); larg->i_idx->first();
- }
-
- // copy output pointers and update them
- data = *larg->data; ++(*larg->data);
- if (larg->impactcoords) impactcoords=larg->impactcoords+((j-1)*larg->res+i-1)*16;
-
-#ifdef HAVE_PTHREAD
- // unlock mutex so our siblings can can access i, j et al. and procede
- if (larg->mutex) pthread_mutex_unlock(larg->mutex);
-#endif
-
- ////// 2- do the actual work.
- (*larg->sc)(i, j, &data, impactcoords, ph);
- ++count;
- }
-#ifdef HAVE_PTHREAD
- if (larg->mutex) {
- delete ph;
- pthread_mutex_lock(larg->mutex);
- }
- GYOTO_MSG << "\nThread terminating after integrating " << count << " photons";
- if (larg->mutex) pthread_mutex_unlock(larg->mutex);
-# endif
- return NULL;
-}
-
YGYOTO_YUSEROBJ(Scenery, Scenery)
extern "C" {
@@ -277,9 +172,34 @@ extern "C" {
if (i_idx.isNuller()) return;
Idx j_idx (piargs[1], res);
if (j_idx.isNuller()) return;
- long ni=i_idx.getNElements();
- long nj=j_idx.getNElements();
- long nelem=ni*nj;
+
+ bool is_double=false;
+
+ long ndims, dims[Y_DIMSIZE]={0};
+ long nelem;
+
+ if (i_idx.isDouble() || j_idx.isDouble()) {
+ if (!i_idx.isDouble() || !j_idx.isDouble())
+ throwError("i and j must be of same type (double or long)");
+ is_double=true;
+ nelem=i_idx.getNElements();
+ long const *idims=i_idx.getDims();
+ for (int m=0; m<idims[0]+1; ++m) dims[m]=idims[m];
+ if (j_idx.getNElements() != nelem) {
+ if (nelem==1) {
+ nelem=j_idx.getNElements();
+ idims=j_idx.getDims();
+ for (int m=0; m<idims[0]+1; ++m) dims[m]=idims[m];
+ }
+ else throwError("alpha and delta must be conformable");
+ }
+ } else {
+ nelem=i_idx.getNElements()*j_idx.getNElements();
+ dims[0]=0;
+ if (precompute) dims[++dims[0]]=16;
+ if (i_idx.getNDims()) dims[++dims[0]]=i_idx.getNElements();
+ if (j_idx.getNDims()) dims[++dims[0]]=j_idx.getNElements();
+ }
long nk = 0; int rquant = 0; ystring_t * squant = NULL;
@@ -347,26 +267,19 @@ extern "C" {
if (has_bsp) nk+=nbnuobs-1;
if (!has_sp && !has_bsp) nbnuobs=0;
- long ndims=i_idx.getNDims()+j_idx.getNDims()
- + (precompute ? 1 : (((rquant>=1)||has_sp||has_bsp)));
- GYOTO_DEBUG << "i_idx.getNDims()=" << i_idx.getNDims()
- << ", j_idx.getNDims()" << j_idx.getNDims()
- << ", precompute ? 1 : (((rquant>=1)||has_sp||has_bsp))"
- << (precompute ? 1 : (((rquant>=1)||has_sp||has_bsp))) <<endl;
- long dims[4]={ndims};
- size_t offset=0;
- if (precompute) dims[++offset]=16;
- if (i_idx.getNDims()) dims[++offset]=ni;
- if (j_idx.getNDims()) dims[++offset]=nj;
- if ( !precompute && ((rquant>=1)||has_sp||has_bsp) ) dims[++offset]=nk;
+ if ( !precompute && ((rquant>=1)||has_sp||has_bsp) ) dims[++dims[0]]=nk;
+
GYOTO_DEBUG << "precompute=" << precompute << ", nk=" << nk
<< ", nbnuobs="<<nbnuobs << ", data=ypush_d({"<< dims[0]
<< ", " << dims[1] << ", " << dims[2] << ", " << dims[3]
<< "})\n";
+
+ GYOTO_DEBUG_ARRAY(dims, Y_DIMSIZE);
+
double * data=ypush_d(dims);
Astrobj::Properties prop;
- prop.alloc=Astrobj::Properties::minimal;
+ prop.alloc=false;
SmartPointer<Screen> screen = (*OBJ) -> screen();
# ifdef HAVE_UDUNITS
if (data) (*OBJ)->setPropertyConverters(&prop);
@@ -435,76 +348,52 @@ extern "C" {
} else y_errorq("unknown quantity: %s", squant[k]);
}
}
- if (i_idx.isDouble() ||j_idx.isDouble()) {
- prop.init(nbnuobs);
- SmartPointer<Photon> ph=(*OBJ)->clonePhoton();
- double coord[8];
- screen -> getRayCoord(i_idx.getDVal(), j_idx.getDVal(), coord);
- 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);
+
+ Screen::Coord2dSet *ijspec=NULL;
+ Screen::Coord1dSet *ispec=NULL, *jspec=NULL;
+
+ if (is_double) {
+ size_t sz = i_idx.getNElements() >= j_idx.getNElements()?
+ i_idx.getNElements():j_idx.getNElements();
+ if (i_idx.getNElements()==1)
+ ispec = new Screen::RepeatAngle(i_idx.getDVal(), sz);
+ else
+ ispec = new Screen::Angles(i_idx.getDoubleBuffer(), sz);
+ if (j_idx.getNElements()==1)
+ jspec = new Screen::RepeatAngle(j_idx.getDVal(), sz);
+ else
+ jspec = new Screen::Angles(j_idx.getDoubleBuffer(), sz);
+ ijspec = new Screen::Bucket(*ispec, *jspec);
} else {
- // i and j specify integers; at least one is a list.
- screen -> computeBaseVectors();
- double coord[8];
- screen -> getRayCoord(size_t(i_idx.first()),
- size_t(j_idx.first()),
- coord);
- SmartPointer<Photon> ph=(*OBJ)->clonePhoton();
- ph->setInitCoord(coord, -1);
- ph->spectrometer(screen->spectrometer());
- ph->freqObs(screen->freqObs());
- ph->delta((*OBJ)->delta());
-
- ySceneryThreadWorkerArg larg;
- larg.sc=(*OBJ);
- larg.ph=ph;
- larg.data=∝
- larg.impactcoords=impactcoords;
- larg.i_idx=&i_idx;
- larg.j_idx=&j_idx;
- larg.res=res;
-
-# ifdef HAVE_PTHREAD
- larg.mutex = NULL;
- pthread_mutex_t mumu = PTHREAD_MUTEX_INITIALIZER;
- pthread_t * threads = NULL;
- pthread_t pself = pthread_self();
- larg.parent = &pself;
- size_t nthreads = (*OBJ) -> nThreads();
- if (nthreads >= 2) {
- threads = new pthread_t[nthreads-1];
- larg.mutex = &mumu;
- for (size_t th=0; th < nthreads-1; ++th) {
- if (pthread_create(threads+th, NULL,
- ySceneryThreadWorker,
- static_cast<void*>(&larg)) < 0)
- y_error("Error creating thread");
- }
- }
-# endif
-
- // Call worker on the parent thread
- (*ySceneryThreadWorker)(static_cast<void*>(&larg));
-
-# ifdef HAVE_PTHREAD
- // Wait for the child threads
- if (nthreads>=2)
- for (size_t th=0; th < nthreads-1; ++th)
- pthread_join(threads[th], NULL);
-# endif
- if (impactcoords==NULL) GYOTO_MSG << endl;
+
+ if (i_idx.isRangeOrScalar())
+ ispec = new Screen::Range
+ (i_idx.range_min(), i_idx.range_max(), i_idx.range_dlt());
+ else
+ ispec = new Screen::Indices
+ (reinterpret_cast<size_t const*const>(i_idx.getBuffer()),
+ i_idx.getNElements());
+
+ if (j_idx.isRangeOrScalar())
+ jspec = new Screen::Range
+ (j_idx.range_min(), j_idx.range_max(), j_idx.range_dlt());
+ else
+ jspec = new Screen::Indices
+ (reinterpret_cast<size_t const*const>(j_idx.getBuffer()),
+ j_idx.getNElements());
+
+ ijspec = new Screen::Grid(*ispec, *jspec, "\rj = ");
+ if (verbose() >= GYOTO_QUIET_VERBOSITY)
+ cout << "\nj = 1/" << jspec->size() << flush;
}
- }
- }
+
+ (*OBJ)->rayTrace(*ijspec, &prop, impactcoords);
+
+ delete ispec;
+ delete jspec;
+ delete ijspec;
+ } // if (conditions for ray-tracing)
+ } // void gyoto_Scenery_eval(void *obj, int argc);
// Constructor
void Y_gyoto_Scenery(int argc) {
@@ -512,7 +401,7 @@ extern "C" {
gyoto_Scenery_eval(OBJ, argc);
}
- void Y_gyoto_Scenery_rayTrace(int argc) {
+void Y_gyoto_Scenery_rayTrace(int argc) {
size_t imin=1, imax=-1, jmin=1, jmax=-1;
if (argc<1) y_error("gyoto_Scenery_rayTrace takes at least 1 argument");
gyoto_Scenery * s_obj=(gyoto_Scenery*)yget_obj(argc-1, &gyoto_Scenery_obj);
@@ -526,6 +415,13 @@ extern "C" {
try {res=long(scenery->screen()->resolution());}
YGYOTO_STD_CATCH;
+ if (imax>res) imax=res;
+ if (jmax>res) jmax=res;
+
+ Screen::Range irange(imin, imax, 1);
+ Screen::Range jrange(jmin, jmax, 1);
+ Screen::Grid grid(irange, jrange, "\r j = ");
+
double * impactcoords = NULL;
int ipctout = 0;
if (argc>=6) {
@@ -598,8 +494,10 @@ extern "C" {
data.intensity=vect;
- try {scenery -> rayTrace(imin, imax, jmin, jmax, &data,
- ipctout?NULL:impactcoords);}
+ if (verbose() >= GYOTO_QUIET_VERBOSITY)
+ cout << endl << flush;
+
+ try {scenery -> rayTrace(grid, &data, ipctout?NULL:impactcoords);}
YGYOTO_STD_CATCH;
}
diff --git a/yorick/ygyoto_idx.h b/yorick/ygyoto_idx.h
index 8e7d774..bbf5d5d 100644
--- a/yorick/ygyoto_idx.h
+++ b/yorick/ygyoto_idx.h
@@ -20,6 +20,8 @@
#ifndef __YGYOTO_IDX_H
#define __YGYOTO_IDX_H
+#include "yapi.h"
+
namespace YGyoto { class Idx; }
class YGyoto::Idx {
@@ -29,10 +31,13 @@ private:
int _is_list;
int _is_scalar;
int _is_double;
+ int _is_dlist;
int _is_first;
long _range[3];
+ long _dims[Y_DIMSIZE];
double _dval;
+ double *_buf;
long *_idx;
long _nel;
long _cur;
@@ -49,12 +54,16 @@ public:
long current() const;
double getDVal() const;
int isDouble() const;
+ int isDoubleArray() const;
int isRangeOrScalar() const;
int isFirst() const;
int isLast() const;
long range_min() const;
long range_max() const;
long range_dlt() const;
+ long const * getBuffer() const;
+ double const * getDoubleBuffer() const;
+ long const * getDims() const;
};
diff --git a/yorick/yorick1.in b/yorick/yorick1.in
index 986148a..0ead7c4 100644
--- a/yorick/yorick1.in
+++ b/yorick/yorick1.in
@@ -26,6 +26,7 @@ abs_top_builddir=@abs_top_builddir@
YORICK=@YORICK@
export ${DYLIB_VAR}=${abs_top_builddir}/lib/.libs:${!DYLIB_VAR}
+export PATH=${abs_top_builddir}/bin:${PATH}
${YORICK} -i ${abs_top_builddir}/yorick/setpaths.i $@
--
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