eic-opticks/sysrap/iexpand.h

0001 /*
0002  * Copyright (c) 2019 Opticks Team. All Rights Reserved.
0003  *
0004  * This file is part of Opticks
0005  * (see https://bitbucket.org/simoncblyth/opticks).
0006  *
0007  * Licensed under the Apache License, Version 2.0 (the "License");
0008  * you may not use this file except in compliance with the License.
0009  * You may obtain a copy of the License at
0010  *
0011  *   http://www.apache.org/licenses/LICENSE-2.0
0012  *
0013  * Unless required by applicable law or agreed to in writing, software
0014  * distributed under the License is distributed on an "AS IS" BASIS,
0015  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0016  * See the License for the specific language governing permissions and
0017  * limitations under the License.
0018  */
0019
0020 #pragma once
0021 /**
0022
0023 iexpand.h
0024 ===========
0025
0026 Adapted from  /usr/local/env/numerics/thrust/examples/expand.cu
0027
0028 Expand an input sequence of counts by replicating indices of each element the number
0029 of times specified by the count values.
0030
0031 The element counts are assumed to be non-negative integers.
0032
0033 Note that the length of the output is equal
0034 to the sum of the input counts.
0035
0036 For example::
0037
0038     iexpand([2,2,2]) -> [0,0,1,1,2,2]  2*0, 2*1, 2*2
0039     iexpand([3,0,1]) -> [0,0,0,2]      3*0, 0*1, 1*2
0040     iexpand([1,3,2]) -> [0,1,1,1,2,2]  1*0, 3*1, 2*2
0041
0042 NB the output device must be zeroed prior to calling iexpand.
0043 This is because the iexpand is implemented ending with an inclusive_scan
0044 to fill in the non-transition values which relies on initial zeroing.
0045
0046
0047 A more specific example:
0048
0049 Every optical photon generating genstep (Cerenkov or scintillation)
0050 specifies the number of photons it will generate.
0051 Applying iexpand to the genstep photon counts produces
0052 an array of genstep indices that is stored into the seed buffer
0053 and provides a reference back to the genstep that produced it.
0054 The seed values are used to translate from a photon index to a
0055 genstep index.
0056
0057 **/
0058
0059
0060 #include <thrust/device_vector.h>
0061 #include <thrust/reduce.h>
0062 #include <thrust/gather.h>
0063 #include <thrust/scan.h>
0064 #include <thrust/fill.h>
0065 #include <thrust/copy.h>
0066 #include <iostream>
0067
0068 #ifdef DEBUG
0069 #include <cassert>
0070 #include <iterator>
0071 template <typename Vector>
0072 void print(const std::string& s, const Vector& v)
0073 {
0074   typedef typename Vector::value_type T;
0075
0076   std::cout << s;
0077   thrust::copy(v.begin(), v.end(), std::ostream_iterator<T>(std::cout, " "));
0078   std::cout << std::endl;
0079 }
0080 #endif
0081
0082
0083 template <typename InputIterator,
0084           typename OutputIterator>
0085 void iexpand(InputIterator  counts_first,
0086              InputIterator  counts_last,
0087              OutputIterator output_first,
0088              OutputIterator output_last)
0089 {
0090   typedef typename thrust::iterator_difference<InputIterator>::type difference_type;
0091
0092   difference_type counts_size = thrust::distance(counts_first, counts_last);  // eg number of gensteps
0093   difference_type output_size = thrust::distance(output_first, output_last);  // eg number of photon "seeds" : back referencing genstep index
0094
0095 #ifdef DEBUG
0096   std::cout << "iexpand "
0097             << " counts_size " << counts_size
0098             << " output_size " << output_size
0099             << std::endl ;
0100 #endif
0101
0102
0103   thrust::device_vector<difference_type> output_offsets(counts_size, 0);
0104
0105   thrust::exclusive_scan(counts_first, counts_last, output_offsets.begin());
0106 #ifdef DEBUG
0107   print(
0108      " scan the counts to obtain output offsets for each input element \n"
0109      " exclusive_scan of input counts creating output_offsets of transitions \n"
0110      " exclusive_scan is a cumsum that excludes current value \n"
0111      " 1st result element always 0, last input value ignored  \n"
0112      " (output_offsets) \n"
0113    , output_offsets );
0114
0115   difference_type output_size2 = thrust::reduce(counts_first, counts_last);    // sum of input counts
0116   assert( output_size == output_size2 );
0117 #endif
0118
0119   // scatter indices into transition points of output
0120   thrust::scatter_if
0121     (thrust::counting_iterator<difference_type>(0),
0122      thrust::counting_iterator<difference_type>(counts_size),
0123      output_offsets.begin(),
0124      counts_first,
0125      output_first);
0126
0127 #ifdef DEBUG
0128   printf(
0129      " scatter the nonzero counts into their corresponding output positions \n"
0130      " scatter_if( first, last, map, stencil, output ) \n"
0131      "    conditionally copies elements from a source range (indices 0:N-1) into an output array according to a map \n"
0132      "    condition dictated by a stencil (input counts) which must be non-zero to be true \n"
0133      "    map provides indices of where to put the indice values in the output  \n"
0134    );
0135 #endif
0136
0137   thrust::inclusive_scan
0138     (output_first,
0139      output_last,
0140      output_first,
0141      thrust::maximum<difference_type>());
0142
0143 #ifdef DEBUG
0144   printf(
0145      " compute max-scan over the output indices, filling in the holes \n"
0146      " inclusive_scan is a cumsum that includes current value \n"
0147      " providing an binary operator (maximum) replaces the default of plus \n"
0148      "   because the empties are init to 0 this will fill in the empties to the right \n"
0149    );
0150 #endif
0151
0152
0153 }
0154
0155
0156