Merge pull request llr-cta#154 from sfegan/main

Work on simulations

.github/workflows/native-build-and-test.yml

@@ -59,7 +59,7 @@ jobs:
       CXX: ${{ matrix.compiler[1] }}
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v4
     - name: Install prerequisites
       shell: bash
       run: |

.github/workflows/singularity-build-and-deploy.yml

@@ -38,7 +38,7 @@ jobs:
     runs-on: ubuntu-20.04
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v4
 
     - name: Install prerequisites
       shell: bash
@@ -87,7 +87,7 @@ jobs:
         ls -l /tmp/calin_installed.tgz
 
     - name: Set up Go 1.13
-      uses: actions/setup-go@v1
+      uses: actions/setup-go@v4
       with:
         go-version: 1.13
       id: go

.gitignore

@@ -3,3 +3,4 @@
 .#*
 #*#
 .ipynb_checkpoints/
+.vscode/

include/math/fftw_util.hpp

@@ -91,6 +91,103 @@ void hcvec_scale_and_multiply(T* ovec, const T* ivec1,
   if(ro==co)(*ro) = (*ri1) * (*ri2) * scale;
 }
 
+template<typename VCLReal>
+void hcvec_scale_and_multiply_vcl(typename VCLReal::real_t* ovec,
+  const typename VCLReal::real_t* ivec1, const typename VCLReal::real_t* ivec2,
+  unsigned nsample, typename VCLReal::real_t scale = 1.0)
+{
+  typename VCLReal::real_t* ro = ovec;
+  typename VCLReal::real_t* co = ovec + nsample;
+  const typename VCLReal::real_t* ri1 = ivec1;
+  const typename VCLReal::real_t* ci1 = ivec1 + nsample;
+  const typename VCLReal::real_t* ri2 = ivec2;
+  const typename VCLReal::real_t* ci2 = ivec2 + nsample;
+
+  (*ro++) = (*ri1++) * (*ri2++) * scale;
+
+  while(co - ro >= 2*VCLReal::num_real)
+  {
+    co -= VCLReal::num_real;
+    ci1 -= VCLReal::num_real;
+    ci2 -= VCLReal::num_real;
+    typename VCLReal::real_vt vri1; vri1.load(ri1);
+    typename VCLReal::real_vt vci1; vci1.load(ci1);
+    typename VCLReal::real_vt vri2; vri2.load(ri2);
+    typename VCLReal::real_vt vci2; vci2.load(ci2);
+
+    typename VCLReal::real_vt vro = (vri1*vri2 - calin::util::vcl::reverse(vci1*vci2)) * scale;
+    typename VCLReal::real_vt vco = (calin::util::vcl::reverse(vri1)*vci2 + vci1*calin::util::vcl::reverse(vri2)) * scale;
+
+    vro.store(ro);
+    vco.store(co);
+
+    ro += VCLReal::num_real;
+    ri1 += VCLReal::num_real;
+    ri2 += VCLReal::num_real;
+  }
+
+  --co;
+  --ci1;
+  --ci2;
+  while(ro < co)
+  {
+    double vri1 = *ri1++;
+    double vci1 = *ci1--;
+    double vri2 = *ri2++;
+    double vci2 = *ci2--;
+    (*ro++) = (vri1*vri2 - vci1*vci2) * scale;
+    (*co--) = (vri1*vci2 + vci1*vri2) * scale;
+  }
+  if(ro==co)(*ro) = (*ri1) * (*ri2) * scale;
+}
+// NOTE : This function overrides the template for systems with AVX !!!
+#if INSTRSET >= 7
+void hcvec_scale_and_multiply(double* ovec, const double* ivec1,
+  const double* ivec2, unsigned nsample, double scale = 1.0);
+#endif
+
+#if 0
+template<typename T>
+void hcvec_scale_and_multiply_with_stride(T* ovec, int ostride,
+  const T* ivec1, int istride1, const T* ivec2, int istride2, unsigned nsample, T scale = 1.0)
+{
+  T *ro = ovec;
+  T *co = ovec + (nsample-1)*ostride;
+  const T *ri1 = ivec1;
+  const T *ci1 = ivec1 + (nsample-1)*istride1;
+  const T *ri2 = ivec2;
+  const T *ci2 = ivec2 + (nsample-1)*istride2;
+
+  *ro = *ri1 * *ri2 * scale;
+  ro += ostride, ri1 += istride1, ri2 += istride2;
+  if(ro==ri1 or ro==ri2)
+  {
+    while(ro < co)
+    {
+      T vri1 = *ri1;
+      T vci1 = *ci1;
+      T vri2 = *ri2;
+      T vci2 = *ci2;
+      *ro = (vri1*vri2 - vci1*vci2)*scale;
+      *co = (vri1*vci2 + vci1*vri2)*scale;
+      ro += ostride, ri1 += istride1, ri2 += istride2;
+      co -= ostride, ci1 -= istride1, ci2 -= istride2;
+    }
+   }
+  else
+  {
+    while(ro < co)
+    {
+      *ro = ((*ri1)*(*ri2) - (*ci1)*(*ci2)) * scale;
+      *co = ((*ri1)*(*ci2) + (*ci1)*(*ri2)) * scale;
+      ro += ostride, ri1 += istride1, ri2 += istride2;
+      co -= ostride, ci1 -= istride1, ci2 -= istride2;
+    }
+  }
+  if(ro==co)(*ro) = (*ri1) * (*ri2) * scale;
+}
+#endif
+
 template<typename T>
 void hcvec_scale_and_multiply_conj(T* ovec, const T* ivec1,
   const T* ivec2_conj, unsigned nsample, T scale = 1.0)

include/math/geometry.hpp

@@ -411,7 +411,7 @@ scatter_direction_in_place(Eigen::Vector3d& v, double dispersion_per_axis,
   calin::math::rng::RNG& rng)
 {
   Eigen::Vector3d x;
-  rng.normal_two_bm(x.x(), x.y());
+  rng.normal_pair_bm(x.x(), x.y());
   x.x() *= dispersion_per_axis;
   x.y() *= dispersion_per_axis;
   x.z() = sqrt(1.0-x.x()*x.x()-x.y()*x.y());
@@ -559,6 +559,12 @@ inline Eigen::Vector3d norm_of_polynomial_surface(double x, double z, const Eige
   return norm_of_polynomial_surface(x, z, p.data(), p.size());
 }
 
+inline Eigen::Vector3d norm_and_y_of_polynomial_surface(double& y_out, double x, double z,
+  const Eigen::VectorXd& p)
+{
+  return norm_and_y_of_polynomial_surface(y_out, x, z, p.data(), p.size());
+}
+
 inline int find_square_grid_site(double x, double y, double pitch_inv, unsigned nside,
   double xc = 0, double yc = 0, double dead_space_fraction = 0)
 {

include/math/geometry_vcl.hpp

@@ -321,7 +321,7 @@ template<typename VCLReal> struct alignas(VCLReal::vec_bytes) VCL: public VCLRea
     calin::math::rng::VCLRealRNG<VCLReal>& rng)
   {
     vec3_vt x;
-    rng.normal_two_bm(x.x(), x.y());
+    rng.normal_pair_bm(x.x(), x.y());
     x.x() *= dispersion_per_axis;
     x.y() *= dispersion_per_axis;
     // x.z() = sqrt(1.0-x.x()*x.x()-x.y()*x.y());

include/math/nspace.hpp

@@ -29,6 +29,7 @@
 
 #include <Eigen/Dense>
 #include <calin_global_definitions.hpp>
+#include <math/nspace.pb.h>
 
 namespace calin { namespace math { namespace nspace {
 
@@ -46,6 +47,7 @@ class TreeSparseNSpace
   TreeSparseNSpace(const Eigen::VectorXd& xlo, const Eigen::VectorXd& xhi,
     const Eigen::VectorXi& n);
   TreeSparseNSpace(const std::vector<Axis>& axes);
+  TreeSparseNSpace(const calin::ix::math::nspace::NSpaceData& proto);
 
   void clear() { bins_.clear(); }
   void injest(const TreeSparseNSpace& o);
@@ -187,6 +189,11 @@ class TreeSparseNSpace
   Eigen::MatrixXd covar_mean_and_total_weight(Eigen::VectorXd& w1, double& w0) const;
   Eigen::MatrixXd covar() const;
 
+  void save_to_proto(calin::ix::math::nspace::NSpaceData* proto) const;
+  calin::ix::math::nspace::NSpaceData* as_proto() const;
+  void accumulate_from_proto(const calin::ix::math::nspace::NSpaceData& proto);
+  static TreeSparseNSpace* create_from_proto(const calin::ix::math::nspace::NSpaceData& proto);
+
 private:
   Eigen::VectorXd xlo_;
   Eigen::VectorXd xhi_;
@@ -204,6 +211,7 @@ class BlockSparseNSpace
   BlockSparseNSpace(const Eigen::VectorXd& xlo, const Eigen::VectorXd& xhi,
     const Eigen::VectorXi& n, unsigned log2_block_size = 0);
   BlockSparseNSpace(const std::vector<Axis>& axes, unsigned log2_block_size = 0);
+  BlockSparseNSpace(const calin::ix::math::nspace::NSpaceData& proto, unsigned log2_block_size = 0);
 
   ~BlockSparseNSpace();
 
@@ -233,6 +241,9 @@ class BlockSparseNSpace
   bool index_of_bin(const Eigen::VectorXi& ix, int64_t& array_index, int64_t& block_index) const;
   bool bin_coords(Eigen::VectorXi& ix_out, int64_t array_index, int64_t block_index) const;
 
+  int64_t map_index(const Eigen::VectorXi& ix) const;
+  void map_bin_coords(Eigen::VectorXi& ix_out, int64_t index) const;
+
   uint64_t size() const { return N_; }
   uint64_t block_size() const { return block_size_; }
   uint64_t block_array_size() const { return array_.size(); }
@@ -276,6 +287,12 @@ class BlockSparseNSpace
   Eigen::MatrixXd covar_mean_and_total_weight(Eigen::VectorXd& w1, double& w0) const;
   Eigen::MatrixXd covar() const;
 
+  void save_to_proto(calin::ix::math::nspace::NSpaceData* proto) const;
+  calin::ix::math::nspace::NSpaceData* as_proto() const;
+  void accumulate_from_proto(const calin::ix::math::nspace::NSpaceData& proto);
+  static BlockSparseNSpace* create_from_proto(const calin::ix::math::nspace::NSpaceData& proto, 
+    unsigned log2_block_size = 0);
+
 #if 0
   void subspace_covar_mean_and_total_weight(const Eigen::VectorXi& subspace_axes,
     BlockSparseNSpace** w2_space, BlockSparseNSpace** w1_space, BlockSparseNSpace** w0_space);

include/math/ray.hpp

@@ -206,6 +206,10 @@ class Ray
       time_reversal_ok, n, ray_is_close_to_surface, tol, niter);
   }
 
+  void scatter_direction(double dispersion_per_axis, calin::math::rng::RNG& rng) {
+    calin::math::geometry::scatter_direction_in_place(dir_, dispersion_per_axis, rng);
+  }
+
 #ifndef SWIG
   inline Eigen::Vector3d norm_of_polynomial_surface(const double* p, unsigned np) const
   {
@@ -215,6 +219,16 @@ class Ray
   double reflect_from_polynomial_surface(const double* p, unsigned np);
   double reflect_from_rough_polynomial_surface(const double* p, unsigned np,
     double roughness, calin::math::rng::RNG& rng);
+
+  void refract_at_polynomial_surface_in(const double* p, unsigned np, double n)
+  {
+    refract_at_surface_in(norm_of_polynomial_surface(p,np), n);
+  }
+
+  void refract_at_polynomial_surface_out(const double* p, unsigned np, double n)
+  {
+    refract_at_surface_out(norm_of_polynomial_surface(p,np), n);
+  }
 #endif
 
   Eigen::Vector3d norm_of_polynomial_surface(const Eigen::VectorXd& p) const
@@ -233,6 +247,16 @@ class Ray
     return reflect_from_rough_polynomial_surface(p.data(), p.size(), roughness, rng);
   }
 
+  void refract_at_polynomial_surface_in(const Eigen::VectorXd& p, double n)
+  {
+    refract_at_polynomial_surface_in(p.data(), p.size(), n);
+  }
+
+  void refract_at_polynomial_surface_out(const Eigen::VectorXd& p, double n)
+  {
+    refract_at_polynomial_surface_out(p.data(), p.size(), n);
+  }
+
 private:
   Eigen::Vector3d pos_ = Eigen::Vector3d::Zero();
   Eigen::Vector3d dir_ = Eigen::Vector3d::Zero();

include/math/ray_vcl.hpp

@@ -124,6 +124,10 @@ template<typename VCLReal> class alignas(VCLReal::vec_bytes) VCLRay: public VCLR
         pos_, ux_inv_, uy_inv_, uz_inv_);
   }
 
+  void scatter_direction(real_vt dispersion_per_axis, calin::math::rng::VCLRealRNG<VCLReal>& rng) {
+    calin::math::geometry::scatter_direction_in_place(dir_, dispersion_per_axis, rng);
+  }
+
   void reflect_from_surface_with_mask(const bool_vt& mask, const vec3_vt& surface_norm) {
     clear_dir_inv();
     dir_ -= surface_norm * select(mask, 2.0*(dir_.dot(surface_norm)), 0);

include/math/rng.hpp

@@ -132,11 +132,22 @@ class RNG
   void uniform_by_type(float& x) { x = uniform_float(); }
 
   double uniform() { return uniform_double(); }
-  double exponential() { return -std::log(uniform()); }
-  double exponential(double mean) { return -mean*std::log(uniform()); }
-  double normal();
+
+  double exponential() { return exponential_ziggurat(); }
+  double exponential(double mean) { return mean*exponential(); }
+  double exponential_transformation() { return -std::log(uniform()); }
+  double exponential_ziggurat();
+
+  double normal() { return normal_ziggurat(); }
   double normal(double mean, double sigma) { return mean+normal()*sigma; }
-  void normal_two_bm(double& x, double& y);
+  double normal_bm();
+  void normal_pair_bm(double& x, double& y);
+  double normal_ziggurat();
+
+  double x_exp_minus_x_squared() { return x_exp_minus_x_squared_ziggurat(); }
+  double x_exp_minus_x_squared_transformation() { return std::sqrt(-std::log(uniform())); }
+  double x_exp_minus_x_squared_ziggurat();
+
   double gamma_by_alpha_and_beta(double alpha, double beta);
   double gamma_by_mean_and_sigma(const double mean, const double sigma) {
     double b = mean/(sigma*sigma);
@@ -271,24 +282,24 @@ class Ranlux48RNGCore: public RNGCore
     res <<= 16;
     switch(dev_blocks_)
     {
-      case 0:
-        dev_ = gen_();
-        gen_calls_++;
-        res |= dev_&0x000000000000FFFFULL;
-        dev_ >>= 16;
-        dev_blocks_ = 2;
-        break;
-      case 1:
-        res |= dev_&0x000000000000FFFFULL;
-        dev_blocks_ = 0;
-        break;
-      case 2:
-        res |= dev_&0x000000000000FFFFULL;
-        dev_ >>= 16;
-        dev_blocks_ = 1;
-        break;
-      default:
-        assert(0);
+    case 0:
+      dev_ = gen_();
+      gen_calls_++;
+      res |= dev_&0x000000000000FFFFULL;
+      dev_ >>= 16;
+      dev_blocks_ = 2;
+      break;
+    case 1:
+      res |= dev_&0x000000000000FFFFULL;
+      dev_blocks_ = 0;
+      break;
+    case 2:
+      res |= dev_&0x000000000000FFFFULL;
+      dev_ >>= 16;
+      dev_blocks_ = 1;
+      break;
+    default:
+      assert(0);
     }
     return res;
   }
@@ -477,4 +488,33 @@ class NR3_EmulateSIMD_RNGCore: public RNGCore
   unsigned ndev_ = NSTREAM;
 };
 
+#ifndef SWIG
+namespace gaussian_ziggurat {
+  extern double v;
+  extern double r;
+  extern double r_inv;
+  extern double xi[257];
+  extern double yi[257];
+} // namespace calin::math::rng::gaussian_ziggurat
+
+namespace exponential_ziggurat {
+  extern double v;
+  extern double r;
+  extern double exp_minus_r;
+  extern double xi[257];
+  extern double yi[257];
+} // namespace calin::math::rng::exponential_ziggurat
+
+namespace x_exp_minus_x_squared_ziggurat {
+  extern double v;
+  extern double xpeak;
+  extern double r;
+  extern double exp_minus_r_squared;
+  extern double xli[257];
+  extern double xri[257];
+  extern double yi[257];
+} // namespace calin::math::rng::exponential_ziggurat
+
+#endif
+
 } } } // namespace calin::math::rng