Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-04-09 07:49:45

0001 #pragma once
0002 
0003 #if defined(__CUDACC__) || defined(__CUDABE__)
0004    #define QAT4_METHOD __device__ __host__ __forceinline__
0005    #define QAT4_FUNCTION  __device__ __host__ __forceinline__
0006 #else
0007    #define QAT4_METHOD
0008    #define QAT4_FUNCTION inline
0009 #endif
0010 
0011 #if defined(__CUDACC__) || defined(__CUDABE__)
0012 #else
0013    #include <iostream>
0014    #include <iomanip>
0015    #include <vector>
0016    #include <sstream>
0017    #include <cstring>
0018    #include <algorithm>
0019 #endif
0020 
0021 #include "squad.h"
0022 
0023 struct qat4
0024 {
0025     quad q0, q1, q2, q3 ;
0026 
0027     QAT4_METHOD void zero()
0028     {
0029         q0.f.x = 0.f ; q0.f.y = 0.f ; q0.f.z = 0.f ; q0.f.w = 0.f ;
0030         q1.f.x = 0.f ; q1.f.y = 0.f ; q1.f.z = 0.f ; q1.f.w = 0.f ;
0031         q2.f.x = 0.f ; q2.f.y = 0.f ; q2.f.z = 0.f ; q2.f.w = 0.f ;
0032         q3.f.x = 0.f ; q3.f.y = 0.f ; q3.f.z = 0.f ; q3.f.w = 0.f ;
0033     }
0034 
0035     QAT4_METHOD void init()
0036     {
0037         q0.f.x = 1.f ; q0.f.y = 0.f ; q0.f.z = 0.f ; q0.f.w = 0.f ;
0038         q1.f.x = 0.f ; q1.f.y = 1.f ; q1.f.z = 0.f ; q1.f.w = 0.f ;
0039         q2.f.x = 0.f ; q2.f.y = 0.f ; q2.f.z = 1.f ; q2.f.w = 0.f ;
0040         q3.f.x = 0.f ; q3.f.y = 0.f ; q3.f.z = 0.f ; q3.f.w = 1.f ;
0041     }
0042 
0043 
0044     // notice that with *right_multiply* the .w column (with identity info) is not used
0045     QAT4_METHOD float3 right_multiply( const float3& v, const float w ) const
0046     {
0047         float3 ret;
0048         ret.x = q0.f.x * v.x + q1.f.x * v.y + q2.f.x * v.z + q3.f.x * w ;
0049         ret.y = q0.f.y * v.x + q1.f.y * v.y + q2.f.y * v.z + q3.f.y * w ;
0050         ret.z = q0.f.z * v.x + q1.f.z * v.y + q2.f.z * v.z + q3.f.z * w ;
0051         return ret;
0052     }
0053     QAT4_METHOD void right_multiply_inplace( float4& v, const float w ) const   // v.w is ignored
0054     {
0055         float x = q0.f.x * v.x + q1.f.x * v.y + q2.f.x * v.z + q3.f.x * w ;
0056         float y = q0.f.y * v.x + q1.f.y * v.y + q2.f.y * v.z + q3.f.y * w ;
0057         float z = q0.f.z * v.x + q1.f.z * v.y + q2.f.z * v.z + q3.f.z * w ;
0058         v.x = x ;
0059         v.y = y ;
0060         v.z = z ;
0061     }
0062     QAT4_METHOD void right_multiply_inplace( float3& v, const float w ) const
0063     {
0064         float x = q0.f.x * v.x + q1.f.x * v.y + q2.f.x * v.z + q3.f.x * w ;
0065         float y = q0.f.y * v.x + q1.f.y * v.y + q2.f.y * v.z + q3.f.y * w ;
0066         float z = q0.f.z * v.x + q1.f.z * v.y + q2.f.z * v.z + q3.f.z * w ;
0067         v.x = x ;
0068         v.y = y ;
0069         v.z = z ;
0070     }
0071 
0072 
0073     /**
0074     qat4::left_multiply
0075     ---------------------
0076 
0077     Note that with *left_multiply* the .w column is used.
0078     SO MUST CLEAR IDENTITY INFO.
0079 
0080     POTENTIAL FOR SLEEPER BUG HERE AS NORMALLY THE FLOAT
0081     VIEW OF INTEGER IDENTITY COLUMNS GIVES VERY SMALL VALUES
0082     SO IT WOULD NOT BE VISIBLE
0083 
0084     This gets used less that *right_multiply* but it is
0085     used for normals in CSG/intersect_leaf.h
0086 
0087     **/
0088     QAT4_METHOD float3 left_multiply( const float3& v, const float w ) const
0089     {
0090         float3 ret;
0091         ret.x = q0.f.x * v.x + q0.f.y * v.y + q0.f.z * v.z + q0.f.w * w ;
0092         ret.y = q1.f.x * v.x + q1.f.y * v.y + q1.f.z * v.z + q1.f.w * w ;
0093         ret.z = q2.f.x * v.x + q2.f.y * v.y + q2.f.z * v.z + q2.f.w * w ;
0094         return ret;
0095     }
0096     QAT4_METHOD void left_multiply_inplace( float4& v, const float w ) const
0097     {
0098         float x = q0.f.x * v.x + q0.f.y * v.y + q0.f.z * v.z + q0.f.w * w ;
0099         float y = q1.f.x * v.x + q1.f.y * v.y + q1.f.z * v.z + q1.f.w * w ;
0100         float z = q2.f.x * v.x + q2.f.y * v.y + q2.f.z * v.z + q2.f.w * w ;
0101         v.x = x ;
0102         v.y = y ;
0103         v.z = z ;
0104     }
0105 
0106     /**
0107     qat4::copy_columns_3x3
0108     -----------------------
0109 
0110     Canonical usage from CSGOptiX/IAS_Builder::Build
0111     Note that the 4th column is not read, so there
0112     is no need to clear any 4th ".w" column identity info.
0113 
0114     ::
0115 
0116                x  y  z  w
0117           q0   0  4  8  -
0118           q1   1  5  9  -
0119           q2   2  6 10  -
0120           q3   3  7 11  -
0121 
0122     **/
0123     QAT4_METHOD void copy_columns_3x4( float* dst ) const
0124     {
0125          dst[0]  = q0.f.x ;
0126          dst[1]  = q1.f.x ;
0127          dst[2]  = q2.f.x ;
0128          dst[3]  = q3.f.x ;
0129 
0130          dst[4]  = q0.f.y ;
0131          dst[5]  = q1.f.y ;
0132          dst[6]  = q2.f.y ;
0133          dst[7]  = q3.f.y ;
0134 
0135          dst[8]  = q0.f.z ;
0136          dst[9]  = q1.f.z ;
0137          dst[10] = q2.f.z ;
0138          dst[11] = q3.f.z ;
0139     }
0140 
0141     QAT4_METHOD qat4(const quad6& gs)   // ctor from genstep, needed for TORCH genstep generation
0142     {
0143         q0.f.x = gs.q2.f.x ;  q0.f.y = gs.q2.f.y ;   q0.f.z = gs.q2.f.z ;  q0.f.w = gs.q2.f.w ;
0144         q1.f.x = gs.q3.f.x ;  q1.f.y = gs.q3.f.y ;   q1.f.z = gs.q3.f.z ;  q1.f.w = gs.q3.f.w ;
0145         q2.f.x = gs.q4.f.x ;  q2.f.y = gs.q4.f.y ;   q2.f.z = gs.q4.f.z ;  q2.f.w = gs.q4.f.w ;
0146         q3.f.x = gs.q5.f.x ;  q3.f.y = gs.q5.f.y ;   q3.f.z = gs.q5.f.z ;  q3.f.w = gs.q5.f.w ;
0147     }
0148 
0149     QAT4_METHOD void write( quad6& gs ) const
0150     {
0151          gs.q2.f.x = q0.f.x ; gs.q2.f.y = q0.f.y ; gs.q2.f.z = q0.f.z ; gs.q2.f.w = q0.f.w ;
0152          gs.q3.f.x = q1.f.x ; gs.q3.f.y = q1.f.y ; gs.q3.f.z = q1.f.z ; gs.q3.f.w = q1.f.w ;
0153          gs.q4.f.x = q2.f.x ; gs.q4.f.y = q2.f.y ; gs.q4.f.z = q2.f.z ; gs.q4.f.w = q2.f.w ;
0154          gs.q5.f.x = q3.f.x ; gs.q5.f.y = q3.f.y ; gs.q5.f.z = q3.f.z ; gs.q5.f.w = q3.f.w ;
0155     }
0156 
0157 
0158 #if defined(__CUDACC__) || defined(__CUDABE__)
0159 #else
0160     QAT4_METHOD qat4()
0161     {
0162         q0.f.x = 1.f ;  q0.f.y = 0.f ;   q0.f.z = 0.f ;  q0.f.w = 0.f ;
0163         q1.f.x = 0.f ;  q1.f.y = 1.f ;   q1.f.z = 0.f ;  q1.f.w = 0.f ;
0164         q2.f.x = 0.f ;  q2.f.y = 0.f ;   q2.f.z = 1.f ;  q2.f.w = 0.f ;
0165         q3.f.x = 0.f ;  q3.f.y = 0.f ;   q3.f.z = 0.f ;  q3.f.w = 1.f ;
0166     }
0167 
0168     QAT4_METHOD bool is_identity() const
0169     {
0170         return
0171             q0.f.x == 1.f && q0.f.y == 0.f && q0.f.z == 0.f && q0.f.w == 0.f  &&
0172             q1.f.x == 0.f && q1.f.y == 1.f && q1.f.z == 0.f && q1.f.w == 0.f  &&
0173             q2.f.x == 0.f && q2.f.y == 0.f && q2.f.z == 1.f && q2.f.w == 0.f  &&
0174             q3.f.x == 0.f && q3.f.y == 0.f && q3.f.z == 0.f && q3.f.w == 1.f ;
0175     }
0176 
0177     QAT4_METHOD bool is_identity(float eps) const
0178     {
0179         return
0180             std::abs(q0.f.x-1.f)<eps && std::abs(q0.f.y)    <eps && std::abs(q0.f.z) < eps     && std::abs(q0.f.w) < eps  &&
0181             std::abs(q1.f.x)<eps     && std::abs(q1.f.y-1.f)<eps && std::abs(q1.f.z) < eps     && std::abs(q1.f.w) < eps  &&
0182             std::abs(q2.f.x)<eps     && std::abs(q2.f.y)<eps     && std::abs(q2.f.z-1.f) < eps && std::abs(q2.f.w) < eps  &&
0183             std::abs(q3.f.x)<eps     && std::abs(q3.f.y)<eps     && std::abs(q3.f.z) < eps     && std::abs(q3.f.w-1.f) < eps ;
0184     }
0185 
0186 
0187     QAT4_METHOD bool is_zero() const
0188     {
0189         return
0190             q0.f.x == 0.f && q0.f.y == 0.f && q0.f.z == 0.f && q0.f.w == 0.f  &&
0191             q1.f.x == 0.f && q1.f.y == 0.f && q1.f.z == 0.f && q1.f.w == 0.f  &&
0192             q2.f.x == 0.f && q2.f.y == 0.f && q2.f.z == 0.f && q2.f.w == 0.f  &&
0193             q3.f.x == 0.f && q3.f.y == 0.f && q3.f.z == 0.f && q3.f.w == 0.f ;
0194     }
0195 
0196     QAT4_METHOD qat4(float tx, float ty, float tz)
0197     {
0198         q0.f.x = 1.f ;  q0.f.y = 0.f ;   q0.f.z = 0.f ;  q0.f.w = 0.f ;
0199         q1.f.x = 0.f ;  q1.f.y = 1.f ;   q1.f.z = 0.f ;  q1.f.w = 0.f ;
0200         q2.f.x = 0.f ;  q2.f.y = 0.f ;   q2.f.z = 1.f ;  q2.f.w = 0.f ;
0201         q3.f.x = tx  ;  q3.f.y = ty  ;   q3.f.z = tz  ;  q3.f.w = 1.f ;
0202     }
0203 
0204     QAT4_METHOD qat4(const float* v)
0205     {
0206         if(v)
0207         {
0208             read(v);
0209         }
0210         else
0211         {
0212             init();
0213         }
0214     }
0215     QAT4_METHOD qat4(const double* v) // narrowing
0216     {
0217         if(v)
0218         {
0219             read_narrow(v);
0220         }
0221         else
0222         {
0223             init();
0224         }
0225     }
0226 
0227 
0228     QAT4_METHOD void read(const float* v)
0229     {
0230         q0.f.x = *(v+0)  ;  q0.f.y = *(v+1)  ;   q0.f.z = *(v+2)  ;  q0.f.w = *(v+3) ;
0231         q1.f.x = *(v+4)  ;  q1.f.y = *(v+5)  ;   q1.f.z = *(v+6)  ;  q1.f.w = *(v+7) ;
0232         q2.f.x = *(v+8)  ;  q2.f.y = *(v+9)  ;   q2.f.z = *(v+10) ;  q2.f.w = *(v+11) ;
0233         q3.f.x = *(v+12) ;  q3.f.y = *(v+13) ;   q3.f.z = *(v+14) ;  q3.f.w = *(v+15) ;
0234     }
0235 
0236 
0237     QAT4_METHOD void read_narrow(const double* v)
0238     {
0239         q0.f.x = float(*(v+0))  ;  q0.f.y = float(*(v+1))  ;   q0.f.z = float(*(v+2))  ;  q0.f.w = float(*(v+3)) ;
0240         q1.f.x = float(*(v+4))  ;  q1.f.y = float(*(v+5))  ;   q1.f.z = float(*(v+6))  ;  q1.f.w = float(*(v+7)) ;
0241         q2.f.x = float(*(v+8))  ;  q2.f.y = float(*(v+9))  ;   q2.f.z = float(*(v+10)) ;  q2.f.w = float(*(v+11)) ;
0242         q3.f.x = float(*(v+12)) ;  q3.f.y = float(*(v+13)) ;   q3.f.z = float(*(v+14)) ;  q3.f.w = float(*(v+15)) ;
0243     }
0244 
0245 
0246     QAT4_METHOD float* data()
0247     {
0248         return &q0.f.x ;
0249     }
0250 
0251     QAT4_METHOD const float* cdata() const
0252     {
0253         return &q0.f.x ;
0254     }
0255 
0256     QAT4_METHOD qat4* copy() const
0257     {
0258         return new qat4(cdata());
0259     }
0260 
0261     static QAT4_METHOD qat4* identity()
0262     {
0263         qat4* q = new qat4 ;
0264         q->init() ;
0265         return q ;
0266     }
0267 
0268     static QAT4_METHOD void copy(qat4& b, const qat4& a )
0269     {
0270         b.q0.f.x = a.q0.f.x ; b.q0.f.y = a.q0.f.y ; b.q0.f.z = a.q0.f.z ; b.q0.f.w = a.q0.f.w ;
0271         b.q1.f.x = a.q1.f.x ; b.q1.f.y = a.q1.f.y ; b.q1.f.z = a.q1.f.z ; b.q1.f.w = a.q1.f.w ;
0272         b.q2.f.x = a.q2.f.x ; b.q2.f.y = a.q2.f.y ; b.q2.f.z = a.q2.f.z ; b.q2.f.w = a.q2.f.w ;
0273         b.q3.f.x = a.q3.f.x ; b.q3.f.y = a.q3.f.y ; b.q3.f.z = a.q3.f.z ; b.q3.f.w = a.q3.f.w ;
0274     }
0275 
0276     /**
0277 
0278 
0279             a.q0.f.x   a.q0.f.y   a.q0.f.z   a.q0.f.w         b.q0.f.x   b.q0.f.y   b.q0.f.z   b.q0.f.w
0280 
0281             a.q1.f.x   a.q1.f.y   a.q1.f.z   a.q1.f.w         b.q1.f.x   b.q1.f.y   b.q1.f.z   b.q1.f.w
0282 
0283             a.q2.f.x   a.q2.f.y   a.q2.f.z   a.q2.f.w         b.q2.f.x   b.q2.f.y   b.q2.f.z   b.q2.f.w
0284 
0285             a.q3.f.x   a.q3.f.y   a.q3.f.z   a.q3.f.w         b.q3.f.x   b.q3.f.y   b.q3.f.z   b.q3.f.w
0286 
0287     **/
0288 
0289     QAT4_METHOD qat4(const qat4& a_, const qat4& b_ , bool flip)  // flip=true matches glm A*B
0290     {
0291         const qat4& a = flip ? b_ : a_ ;
0292         const qat4& b = flip ? a_ : b_ ;
0293 
0294         q0.f.x = a.q0.f.x * b.q0.f.x  +  a.q0.f.y * b.q1.f.x  +  a.q0.f.z * b.q2.f.x  + a.q0.f.w * b.q3.f.x ;
0295         q0.f.y = a.q0.f.x * b.q0.f.y  +  a.q0.f.y * b.q1.f.y  +  a.q0.f.z * b.q2.f.y  + a.q0.f.w * b.q3.f.y ;
0296         q0.f.z = a.q0.f.x * b.q0.f.z  +  a.q0.f.y * b.q1.f.z  +  a.q0.f.z * b.q2.f.z  + a.q0.f.w * b.q3.f.z ;
0297         q0.f.w = a.q0.f.x * b.q0.f.w  +  a.q0.f.y * b.q1.f.w  +  a.q0.f.z * b.q2.f.w  + a.q0.f.w * b.q3.f.w ;
0298 
0299         q1.f.x = a.q1.f.x * b.q0.f.x  +  a.q1.f.y * b.q1.f.x  +  a.q1.f.z * b.q2.f.x  + a.q1.f.w * b.q3.f.x ;
0300         q1.f.y = a.q1.f.x * b.q0.f.y  +  a.q1.f.y * b.q1.f.y  +  a.q1.f.z * b.q2.f.y  + a.q1.f.w * b.q3.f.y ;
0301         q1.f.z = a.q1.f.x * b.q0.f.z  +  a.q1.f.y * b.q1.f.z  +  a.q1.f.z * b.q2.f.z  + a.q1.f.w * b.q3.f.z ;
0302         q1.f.w = a.q1.f.x * b.q0.f.w  +  a.q1.f.y * b.q1.f.w  +  a.q1.f.z * b.q2.f.w  + a.q1.f.w * b.q3.f.w ;
0303 
0304         q2.f.x = a.q2.f.x * b.q0.f.x  +  a.q2.f.y * b.q1.f.x  +  a.q2.f.z * b.q2.f.x  + a.q2.f.w * b.q3.f.x ;
0305         q2.f.y = a.q2.f.x * b.q0.f.y  +  a.q2.f.y * b.q1.f.y  +  a.q2.f.z * b.q2.f.y  + a.q2.f.w * b.q3.f.y ;
0306         q2.f.z = a.q2.f.x * b.q0.f.z  +  a.q2.f.y * b.q1.f.z  +  a.q2.f.z * b.q2.f.z  + a.q2.f.w * b.q3.f.z ;
0307         q2.f.w = a.q2.f.x * b.q0.f.w  +  a.q2.f.y * b.q1.f.w  +  a.q2.f.z * b.q2.f.w  + a.q2.f.w * b.q3.f.w ;
0308 
0309         q3.f.x = a.q3.f.x * b.q0.f.x  +  a.q3.f.y * b.q1.f.x  +  a.q3.f.z * b.q2.f.x  + a.q3.f.w * b.q3.f.x ;
0310         q3.f.y = a.q3.f.x * b.q0.f.y  +  a.q3.f.y * b.q1.f.y  +  a.q3.f.z * b.q2.f.y  + a.q3.f.w * b.q3.f.y ;
0311         q3.f.z = a.q3.f.x * b.q0.f.z  +  a.q3.f.y * b.q1.f.z  +  a.q3.f.z * b.q2.f.z  + a.q3.f.w * b.q3.f.z ;
0312         q3.f.w = a.q3.f.x * b.q0.f.w  +  a.q3.f.y * b.q1.f.w  +  a.q3.f.z * b.q2.f.w  + a.q3.f.w * b.q3.f.w ;
0313     }
0314 
0315 
0316     // .w column looks to be in use : potential bug ? : TODO: compare with longhand approach
0317     QAT4_METHOD void transform_aabb_inplace( float* aabb ) const
0318     {
0319         float4 xa = q0.f * *(aabb+0) ;
0320         float4 xb = q0.f * *(aabb+3) ;
0321         float4 xmi = fminf(xa, xb);
0322         float4 xma = fmaxf(xa, xb);
0323 
0324         float4 ya = q1.f * *(aabb+1) ;
0325         float4 yb = q1.f * *(aabb+4) ;
0326         float4 ymi = fminf(ya, yb);
0327         float4 yma = fmaxf(ya, yb);
0328 
0329         float4 za = q2.f * *(aabb+2) ;
0330         float4 zb = q2.f * *(aabb+5) ;
0331         float4 zmi = fminf(za, zb);
0332         float4 zma = fmaxf(za, zb);
0333 
0334         float4 tmi = xmi + ymi + zmi + q3.f ;
0335         float4 tma = xma + yma + zma + q3.f ;
0336 
0337         *(aabb + 0) = tmi.x ;
0338         *(aabb + 1) = tmi.y ;
0339         *(aabb + 2) = tmi.z ;
0340         *(aabb + 3) = tma.x ;
0341         *(aabb + 4) = tma.y ;
0342         *(aabb + 5) = tma.z ;
0343     }
0344 
0345     QAT4_METHOD void getIdentity(int& ins_idx, int& gas_idx, int& sensor_identifier, int& sensor_index ) const
0346     {
0347         ins_idx           = q0.i.w ;   // formerly used unsigned and "- 1"
0348         gas_idx           = q1.i.w ;
0349         sensor_identifier = q2.i.w ;
0350         sensor_index      = q3.i.w  ;
0351     }
0352 
0353     /**
0354     sqat4::get_IAS_OptixInstance_instanceId
0355     ----------------------------------------
0356 
0357     Canonical use by SBT::collectInstances
0358 
0359     * July 2023 : QPMT needs lpmtid GPU side so changed to sensor_identifier from ins_idx
0360     * FORMER BUG: -1 identifier when unsigned became ~0u exceeding permissable bits,
0361       causing change to +1 and using zero to mean not-a-sensor
0362 
0363     YEP: that ~0u caused notes/issues/unspecified_launch_failure_with_simpleLArTPC.rst
0364 
0365     **/
0366 
0367     QAT4_METHOD int get_IAS_OptixInstance_instanceId() const
0368     {
0369         const int& sensor_identifier_1 = q2.i.w ;
0370         assert( sensor_identifier_1 >= 0 );  // 0 means not a sensor GPU side
0371         return sensor_identifier_1 ;    // NB this is +1,  subtract 1 to get original sensorId of sensors
0372     }
0373 
0374     /**
0375     sqat4::setIdentity
0376     -------------------
0377 
0378     Canonical usage from CSGFoundry::addInstance  where sensor_identifier gets +1
0379     with 0 meaning not a sensor.
0380     **/
0381 
0382     QAT4_METHOD void setIdentity(int ins_idx, int gas_idx, int sensor_identifier_1, int sensor_index )
0383     {
0384         assert( sensor_identifier_1 >= 0 );
0385 
0386         q0.i.w = ins_idx ;               // formerly unsigned and "+ 1"
0387         q1.i.w = gas_idx ;
0388         q2.i.w = sensor_identifier_1 ;   // now +1 with 0 meaning not-a-sensor
0389         q3.i.w = sensor_index ;
0390     }
0391 
0392     /**
0393     sqat4::incrementSensorIdentifier
0394     ---------------------------------
0395 
0396     Canonical usage from CSGFoundry::addInstanceVector to adjust from
0397 
0398     * CPU -1:not-a-sensor
0399     * GPU  0:not-a-sensor
0400 
0401     **/
0402     QAT4_METHOD void incrementSensorIdentifier()
0403     {
0404         assert( q2.i.w >= -1 );
0405         q2.i.w += 1 ;
0406         assert( q2.i.w >=  0 );
0407     }
0408 
0409 
0410     QAT4_METHOD void clearIdentity() // prepare for matrix multiply by clearing any auxiliary info in the "spare" 4th column
0411     {
0412         q0.f.w = 0.f ;
0413         q1.f.w = 0.f ;
0414         q2.f.w = 0.f ;
0415         q3.f.w = 1.f ;
0416     }
0417 
0418     static QAT4_METHOD bool IsDiff( const qat4& a, const qat4& b )
0419     {
0420         return false ;
0421     }
0422 
0423     // collects unique ins/gas/ias indices found in qv vector of instances
0424     static QAT4_METHOD void find_unique(const std::vector<qat4>& qv,
0425               std::vector<int>& ins,
0426               std::vector<int>& gas,
0427               std::vector<int>& s_ident,
0428               std::vector<int>& s_index )
0429     {
0430          for(unsigned i=0 ; i < qv.size() ; i++)
0431          {
0432              const qat4& q = qv[i] ;
0433              int ins_idx,  gas_idx, sensor_identifier, sensor_index ;
0434              q.getIdentity(ins_idx,  gas_idx, sensor_identifier, sensor_index );
0435 
0436              if(std::find(ins.begin(), ins.end(), ins_idx) == ins.end() ) ins.push_back(ins_idx);
0437              if(std::find(gas.begin(), gas.end(), gas_idx) == gas.end() ) gas.push_back(gas_idx);
0438              if(std::find(s_ident.begin(), s_ident.end(), sensor_identifier) == s_ident.end() ) s_ident.push_back(sensor_identifier);
0439              if(std::find(s_index.begin(), s_index.end(), sensor_index)      == s_index.end() ) s_index.push_back(sensor_index);
0440 
0441          }
0442     }
0443 
0444     static QAT4_METHOD void find_unique_gas(const std::vector<qat4>& qv, std::vector<int>& gas )
0445     {
0446          for(unsigned i=0 ; i < qv.size() ; i++)
0447          {
0448              const qat4& q = qv[i] ;
0449              int ins_idx,  gas_idx, sensor_identifier, sensor_index ;
0450              q.getIdentity(ins_idx,  gas_idx, sensor_identifier, sensor_index );
0451              if(std::find(gas.begin(), gas.end(), gas_idx) == gas.end() ) gas.push_back(gas_idx);
0452          }
0453     }
0454 
0455 
0456     // count the number of instances with the provided ias_idx, that are among the emm if that is non-zero
0457     static QAT4_METHOD unsigned count_ias( const std::vector<qat4>& qv , int /*ias_idx*/, unsigned long long emm)
0458     {
0459         unsigned count = 0 ;
0460         for(unsigned i=0 ; i < qv.size() ; i++)
0461         {
0462             const qat4& q = qv[i] ;
0463             int ins_idx,  gas_idx, sensor_identifier, sensor_index  ;
0464             q.getIdentity(ins_idx,  gas_idx, sensor_identifier, sensor_index );
0465             unsigned long long gas_idx_ull = gas_idx ;
0466             bool gas_enabled = emm == 0ull ? true : emm & ( 0x1ull << gas_idx_ull ) ;
0467             if( gas_enabled ) count += 1 ;
0468         }
0469         return count ;
0470     }
0471     // select instances with the provided ias_idx, that are among the emm if that is non-zero,  ordered as they are found
0472     static QAT4_METHOD void select_instances_ias(const std::vector<qat4>& qv, std::vector<qat4>& select_qv, int /*ias_idx*/, unsigned long long emm  )
0473     {
0474         for(unsigned i=0 ; i < qv.size() ; i++)
0475         {
0476             const qat4& q = qv[i] ;
0477             int ins_idx,  gas_idx, sensor_identifier, sensor_index  ;
0478             q.getIdentity(ins_idx,  gas_idx, sensor_identifier, sensor_index );
0479             unsigned long long gas_idx_ull = gas_idx ;
0480             bool gas_enabled = emm == 0ull ? true : emm & ( 0x1ull << gas_idx_ull ) ;
0481             if( gas_enabled ) select_qv.push_back(q) ;
0482         }
0483     }
0484 
0485 
0486     // count the number of instances with the provided gas_idx
0487     static QAT4_METHOD unsigned count_gas( const std::vector<qat4>& qv , int gas_idx_ )
0488     {
0489         unsigned count = 0 ;
0490         for(unsigned i=0 ; i < qv.size() ; i++)
0491         {
0492             const qat4& q = qv[i] ;
0493             int ins_idx,  gas_idx, sensor_identifier, sensor_index  ;
0494             q.getIdentity(ins_idx,  gas_idx, sensor_identifier, sensor_index );
0495             if( gas_idx_ == gas_idx ) count += 1 ;
0496         }
0497         return count ;
0498     }
0499     // select instances with the provided gas_idx, ordered as they are found
0500     static QAT4_METHOD void select_instances_gas(const std::vector<qat4>& qv, std::vector<qat4>& select_qv, int gas_idx_ )
0501     {
0502         for(unsigned i=0 ; i < qv.size() ; i++)
0503         {
0504             const qat4& q = qv[i] ;
0505             int ins_idx,  gas_idx, sensor_identifier, sensor_index  ;
0506             q.getIdentity(ins_idx,  gas_idx, sensor_identifier, sensor_index );
0507             if( gas_idx_ == gas_idx ) select_qv.push_back(q) ;
0508         }
0509     }
0510 
0511     // select instance pointers with the provided gas_idx, ordered as they are found
0512     static QAT4_METHOD void select_instance_pointers_gas(const std::vector<qat4>& qv, std::vector<const qat4*>& select_qi, int gas_idx_ )
0513     {
0514         for(unsigned i=0 ; i < qv.size() ; i++)
0515         {
0516             const qat4* qi = qv.data() + i ;
0517             int ins_idx,  gas_idx, sensor_identifier, sensor_index  ;
0518             qi->getIdentity(ins_idx,  gas_idx, sensor_identifier, sensor_index );
0519             if( gas_idx_ == gas_idx ) select_qi.push_back(qi) ;
0520         }
0521     }
0522 
0523     // return absolute instance index of the ordinal-th instance with the provided gas_idx or -1 if not found
0524     static QAT4_METHOD int find_instance_gas(const std::vector<qat4>& qv, int gas_idx_, unsigned ordinal  )
0525     {
0526         unsigned count = 0 ;
0527         int index = -1 ;
0528         for(unsigned i=0 ; i < qv.size() ; i++)
0529         {
0530             const qat4& q = qv[i] ;
0531             int ins_idx,  gas_idx, sensor_identifier, sensor_index  ;
0532             q.getIdentity(ins_idx,  gas_idx, sensor_identifier, sensor_index );
0533             if( gas_idx_ == gas_idx )
0534             {
0535                 if( count == ordinal ) index = i  ;
0536                 count += 1 ;
0537             }
0538         }
0539         return index ;
0540     }
0541 
0542     QAT4_METHOD void right_multiply_inplace( std::vector<float3>& points, const float w ) const
0543     {
0544         for(int i=0 ; i < int(points.size()) ; i++)
0545         {
0546             float3& p = points[i] ;
0547             right_multiply_inplace( p, w );
0548         }
0549     }
0550 
0551 
0552     QAT4_METHOD float element(unsigned j, unsigned k) const
0553     {
0554         const float4& f = j == 0 ? q0.f : ( j == 1 ? q1.f : ( j == 2 ? q2.f : q3.f ) ) ;
0555         return   k == 0  ? f.x : ( k == 1 ? f.y : ( k == 2 ? f.z : f.w )) ;
0556     }
0557 
0558     QAT4_METHOD void add_translate(float tx, float ty, float tz)
0559     {
0560         q3.f.x += tx ;
0561         q3.f.y += ty ;
0562         q3.f.z += tz ;
0563     }
0564 
0565 
0566     static QAT4_METHOD int compare(const qat4& a, const qat4& b, float epsilon=1e-6 )
0567     {
0568         int rc = 0 ;
0569         if( std::abs(a.q0.f.x - b.q0.f.x) > epsilon ) rc += 1 ;
0570         if( std::abs(a.q0.f.y - b.q0.f.y) > epsilon ) rc += 1 ;
0571         if( std::abs(a.q0.f.z - b.q0.f.z) > epsilon ) rc += 1 ;
0572         if( std::abs(a.q0.f.w - b.q0.f.w) > epsilon ) rc += 1 ;
0573 
0574         if( std::abs(a.q1.f.x - b.q1.f.x) > epsilon ) rc += 1 ;
0575         if( std::abs(a.q1.f.y - b.q1.f.y) > epsilon ) rc += 1 ;
0576         if( std::abs(a.q1.f.z - b.q1.f.z) > epsilon ) rc += 1 ;
0577         if( std::abs(a.q1.f.w - b.q1.f.w) > epsilon ) rc += 1 ;
0578 
0579         if( std::abs(a.q2.f.x - b.q2.f.x) > epsilon ) rc += 1 ;
0580         if( std::abs(a.q2.f.y - b.q2.f.y) > epsilon ) rc += 1 ;
0581         if( std::abs(a.q2.f.z - b.q2.f.z) > epsilon ) rc += 1 ;
0582         if( std::abs(a.q2.f.w - b.q2.f.w) > epsilon ) rc += 1 ;
0583 
0584         if( std::abs(a.q3.f.x - b.q3.f.x) > epsilon ) rc += 1 ;
0585         if( std::abs(a.q3.f.y - b.q3.f.y) > epsilon ) rc += 1 ;
0586         if( std::abs(a.q3.f.z - b.q3.f.z) > epsilon ) rc += 1 ;
0587         if( std::abs(a.q3.f.w - b.q3.f.w) > epsilon ) rc += 1 ;
0588 
0589         return rc ;
0590     }
0591 
0592 
0593     static QAT4_METHOD qat4* from_string(const char* s0, const char* replace="[]()," )
0594     {
0595         char* s1 = strdup(s0);
0596         for(unsigned i=0 ; i < strlen(s1) ; i++) if(strchr(replace, s1[i]) != nullptr) s1[i] = ' ' ;
0597 
0598         std::stringstream ss(s1);
0599         std::string s ;
0600         std::vector<float> vv ;
0601         float v ;
0602 
0603         while(std::getline(ss, s, ' '))
0604         {
0605             if(strlen(s.c_str()) == 0 ) continue;
0606             std::stringstream tt(s);
0607             tt >> v ;
0608             vv.push_back(v);
0609         }
0610 
0611         unsigned num_vv = vv.size() ;
0612         qat4* q = nullptr ;
0613 
0614         if( num_vv == 16 )
0615         {
0616             q = new qat4(vv.data()) ;
0617         }
0618         else if( num_vv == 3 || num_vv == 4 )
0619         {
0620             float tx = vv[0] ;
0621             float ty = vv[1] ;
0622             float tz = vv[2] ;
0623 
0624             q = new qat4(tx, ty, tz) ;
0625         }
0626         return q ;
0627     }
0628 
0629 
0630     QAT4_METHOD std::string desc(char mat='t', unsigned wid=8, unsigned prec=3, bool id=true) const
0631     {
0632         std::stringstream ss ;
0633         ss << desc_(mat, wid, prec)
0634            << " "
0635            << ( id ? descId() : "" )
0636            ;
0637         std::string s = ss.str() ;
0638         return s ;
0639     }
0640     QAT4_METHOD std::string desc_(char mat='t', unsigned wid=8, unsigned prec=3) const
0641     {
0642         std::stringstream ss ;
0643         ss << mat << ":" ;
0644         for(int j=0 ; j < 4 ; j++ )
0645         {
0646             ss << "[" ;
0647             for(int k=0 ; k < 4 ; k++ ) ss << std::setw(wid) << std::fixed << std::setprecision(prec) << element(j,k) << " " ;
0648             ss << "]" ;
0649         }
0650         std::string s = ss.str() ;
0651         return s ;
0652     }
0653 
0654     QAT4_METHOD std::string descId() const
0655     {
0656         int ins_idx,  gas_idx, sensor_identifier, sensor_index ;
0657         getIdentity(ins_idx,  gas_idx, sensor_identifier, sensor_index );
0658 
0659         std::stringstream ss ;
0660         ss
0661             << "( i/g/si/sx "
0662             << std::setw(7) << ins_idx
0663             << std::setw(3) << gas_idx
0664             << std::setw(7) << sensor_identifier
0665             << std::setw(5) << ( sensor_index == 1065353216 ? -1 : sensor_index )  // 1065353216 is 1.f viewed as int
0666             << " )"
0667             ;
0668 
0669         std::string s = ss.str() ;
0670         return s ;
0671     }
0672 
0673     static QAT4_METHOD void dump(const std::vector<qat4>& qv);
0674 
0675 #endif
0676 
0677 };
0678 
0679 
0680 
0681 #if defined(__CUDACC__) || defined(__CUDABE__)
0682 #else
0683 
0684 inline std::ostream& operator<<(std::ostream& os, const qat4& v)
0685 {
0686     os
0687        << v.q0.f
0688        << v.q1.f
0689        << v.q2.f
0690        << v.q3.f
0691        ;
0692     return os;
0693 }
0694 #endif
0695 
0696 
0697 #if defined(__CUDACC__) || defined(__CUDABE__)
0698 #else
0699 
0700 QAT4_FUNCTION void qat4::dump(const std::vector<qat4>& qv)
0701 {
0702     for(unsigned i=0 ; i < qv.size() ; i++)
0703     {
0704         const qat4& q = qv[i] ;
0705         std::cout
0706             << " i " << std::setw(4) << i
0707             << " " << q.descId()
0708             << " " << q.desc()
0709             << std::endl
0710             ;
0711 
0712     }
0713 }
0714 
0715 #endif
0716