Files
2026-07-13 12:47:05 +08:00

571 lines
23 KiB
C++

/* ******************************************************************************
*
*
* This program and the accompanying materials are made available under the
* terms of the Apache License, Version 2.0 which is available at
* https://www.apache.org/licenses/LICENSE-2.0.
*
* See the NOTICE file distributed with this work for additional
* information regarding copyright ownership.
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
* License for the specific language governing permissions and limitations
* under the License.
*
* SPDX-License-Identifier: Apache-2.0
******************************************************************************/
#ifndef LIBND4J_SPECIAL_RANDOM_OPS_H
#define LIBND4J_SPECIAL_RANDOM_OPS_H
#include <execution/Threads.h>
#include <graph/RandomGenerator.h>
#include <helpers/shape.h>
#include <ops/random_ops.h>
#include <ops/specials_cuda.h>
namespace randomOps {
//////////////////////////////////////////////////////////////////////
template <typename T>
class Choice {
public:
method_idx method_X method_XY
static const bool requiresSpecial = true;
#ifdef __CUDACC__
static SD_INLINE SD_DEVICE void specialOpCuda(sd::Pointer state, T const *x, sd::LongType const *xShapeBuffer,
T const *y, sd::LongType const *yShapeBuffer, T *z,
sd::LongType const *zShapeBuffer, T *extraArguments) {
// ... (CUDA implementation remains unchanged)
}
#endif
static inline void specialOp(sd::Pointer state, const T *x, const sd::LongType *xShapeBuffer, const T *y,
const sd::LongType *yShapeBuffer, T *z, const sd::LongType *zShapeBuffer,
T *extraArguments) {
sd::LongType zLength = shape::length(zShapeBuffer);
sd::LongType yLength = shape::length(yShapeBuffer);
int elementsPerThread = zLength / TAD_THRESHOLD;
int _threads = sd::math::sd_max<int>(1, elementsPerThread);
_threads = sd::math::sd_min<int>(_threads, sd::Environment::getInstance().maxThreads());
sd::LongType zRank = shape::rank(zShapeBuffer);
sd::LongType *zShape = shape::shapeOf(zShapeBuffer);
sd::LongType *zStride = shape::stride(zShapeBuffer);
sd::LongType yRank = shape::rank(yShapeBuffer);
sd::LongType *yShape = shape::shapeOf(yShapeBuffer);
sd::LongType *yStride = shape::stride(yShapeBuffer);
sd::LongType *xShape = shape::shapeOf(xShapeBuffer);
sd::LongType xRank = shape::rank(xShapeBuffer);
sd::LongType *xStride = shape::stride(xShapeBuffer);
sd::graph::RandomGenerator *rng = reinterpret_cast<sd::graph::RandomGenerator *>(state);
auto func = PRAGMA_THREADS_FOR {
for (auto e = start; e < stop; e++) {
sd::LongType coords[SD_MAX_RANK];
INDEX2COORDS(e, zRank, zShape, coords);
sd::LongType zOffset;
COORDS2INDEX(zRank, zStride, coords, zOffset);
T prob = rng->relativeT<T>(e);
T cumProb = (T)0.0f;
for (sd::LongType f = 0; f < yLength; f++) {
sd::LongType yCoords[SD_MAX_RANK];
INDEX2COORDS(f, yRank, yShape, yCoords);
sd::LongType yOffset;
COORDS2INDEX(yRank, yStride, yCoords, yOffset);
T relProb = y[yOffset];
cumProb += relProb;
if (prob <= cumProb || f == yLength - 1) {
sd::LongType xCoords[SD_MAX_RANK];
INDEX2COORDS(f, xRank, xShape, xCoords);
sd::LongType xOffset;
COORDS2INDEX(xRank,xStride , xCoords, xOffset);
z[zOffset] = x[xOffset];
break;
}
}
}
};
samediff::Threads::parallel_for(func, 0, zLength, 1, _threads);
}
};
//////////////////////////////////////////////////////////////////////
template <typename T>
class GaussianDistribution {
public:
method_XY method_X method_idx
static const bool requiresSpecial = true;
#ifdef __CUDACC__
static SD_INLINE SD_DEVICE void specialOpCuda(sd::Pointer state, T const *x, sd::LongType const *xShapeBuffer,
T const *y, sd::LongType const *yShapeBuffer, T *z,
sd::LongType const *zShapeBuffer, T *extraArguments) {
// ... (CUDA implementation remains unchanged)
}
#endif
static inline void specialOp(sd::Pointer state, const T *x, const sd::LongType *xShapeBuffer, const T *y,
const sd::LongType *yShapeBuffer, T *z, const sd::LongType *zShapeBuffer,
T *extraArguments) {
const T two_pi = static_cast<T>(2.0f) * static_cast<T>(3.14159265358979323846);
sd::LongType zLength = shape::length(zShapeBuffer);
auto middle = zLength % 2 + zLength / 2;
int elementsPerThread = middle / TAD_THRESHOLD;
int _threads = sd::math::sd_max<int>(1, elementsPerThread);
_threads = sd::math::sd_min<int>(_threads, sd::Environment::getInstance().maxThreads());
sd::graph::RandomGenerator *rng = reinterpret_cast<sd::graph::RandomGenerator *>(state);
const T mean = extraArguments[0];
const T stddev = extraArguments[1];
const T epsilon = static_cast<T>(1e-5);
sd::LongType zRank = shape::rank(zShapeBuffer);
sd::LongType *zShape = shape::shapeOf(zShapeBuffer);
sd::LongType *zStride = shape::stride(zShapeBuffer);
sd::LongType yRank = shape::rank(yShapeBuffer);
sd::LongType *yShape = shape::shapeOf(yShapeBuffer);
sd::LongType *yStride = shape::stride(yShapeBuffer);
sd::LongType *xShape = shape::shapeOf(xShapeBuffer);
sd::LongType xRank = shape::rank(xShapeBuffer);
sd::LongType *xStride = shape::stride(xShapeBuffer);
auto func = PRAGMA_THREADS_FOR {
for (auto e = start; e < stop; e++) {
sd::LongType coords[SD_MAX_RANK];
INDEX2COORDS(e, zRank, zShape, coords);
sd::LongType zOffset;
COORDS2INDEX(zRank, zStride, coords, zOffset);
auto epm = e + middle;
// we need to get random values
T r0 = rng->relativeT<T>(e, epsilon, static_cast<T>(1.0f));
T r1 = rng->relativeT<T>(epm, epsilon, static_cast<T>(1.0f));
sd::LongType yOffset;
COORDS2INDEX(yRank, yStride, coords, yOffset);
T realMean0 = y == z ? mean : y[yOffset];
z[zOffset] = (sd::math::sd_sqrt<T, T>(static_cast<T>(-2.0f) * sd::math::sd_log<T, T>(r0)) *
sd::math::sd_cos<T, T>(two_pi * r1)) * stddev + realMean0;
if (epm < zLength) {
INDEX2COORDS(epm, zRank, zShape, coords);
COORDS2INDEX(zRank, zStride, coords, zOffset);
COORDS2INDEX(yRank, yStride, coords, yOffset);
T realMean1 = y == z ? mean : y[yOffset];
z[zOffset] = (sd::math::sd_sqrt<T, T>(static_cast<T>(-2.0f) * sd::math::sd_log<T, T>(r0)) *
sd::math::sd_sin<T, T>(two_pi * r1)) * stddev + realMean1;
}
}
};
samediff::Threads::parallel_for(func, 0, middle, 1, _threads);
}
};
//////////////////////////////////////////////////////////////////////
template <typename T>
class BinomialDistribution {
public:
method_XY method_X method_idx
static const bool requiresSpecial = true;
#ifdef __CUDACC__
static SD_INLINE SD_DEVICE void specialOpCuda(sd::Pointer state, T const *x, sd::LongType const *xShapeBuffer,
T const *y, sd::LongType const *yShapeBuffer, T *z,
sd::LongType const *zShapeBuffer, T *extraArguments) {
// ... (CUDA implementation remains unchanged)
}
#endif
static inline void specialOp(sd::Pointer state, const T *x, const sd::LongType *xShapeBuffer, const T *y,
const sd::LongType *yShapeBuffer, T *z, const sd::LongType *zShapeBuffer,
T *extraArguments) {
int trials = (int)extraArguments[0];
sd::LongType zLength = shape::length(zShapeBuffer);
int elementsPerThread = zLength / TAD_THRESHOLD;
int _threads = sd::math::sd_max<int>(1, elementsPerThread);
_threads = sd::math::sd_min<int>(_threads, sd::Environment::getInstance().maxThreads());
T prob = extraArguments[1];
sd::LongType zRank = shape::rank(zShapeBuffer);
sd::LongType *zShape = shape::shapeOf(zShapeBuffer);
sd::LongType *zStride = shape::stride(zShapeBuffer);
sd::LongType yRank = shape::rank(yShapeBuffer);
sd::LongType *yShape = shape::shapeOf(yShapeBuffer);
sd::LongType *yStride = shape::stride(yShapeBuffer);
sd::LongType *xShape = shape::shapeOf(xShapeBuffer);
sd::LongType xRank = shape::rank(xShapeBuffer);
sd::LongType *xStride = shape::stride(xShapeBuffer);
sd::graph::RandomGenerator *rng = reinterpret_cast<sd::graph::RandomGenerator *>(state);
auto func = PRAGMA_THREADS_FOR {
for (auto e = start; e < stop; e++) {
sd::LongType coords[SD_MAX_RANK];
INDEX2COORDS(e,zRank, zShape, coords);
sd::LongType zOffset;
COORDS2INDEX(zRank, zStride, coords, zOffset);
int success = 0;
for (int t = 1; t <= trials; t++) {
T randVal = rng->relativeT<T>((e + 1) * t);
if (y != z) {
// we're using external probs
sd::LongType yOffset;
COORDS2INDEX(yRank,yStride, coords, yOffset);
prob = y[yOffset];
}
if (randVal < prob) success++;
}
// if trials is set to 0, effectively we just have successful memset
z[zOffset] = static_cast<T>(success);
}
};
samediff::Threads::parallel_for(func, 0, zLength, 1, _threads);
}
};
//////////////////////////////////////////////////////////////////////
template <typename T>
class BinomialDistributionEx {
public:
method_XY method_X method_idx
static const bool requiresSpecial = true;
#ifdef __CUDACC__
static SD_INLINE SD_DEVICE void specialOpCuda(sd::Pointer state, T const *x, sd::LongType const *xShapeBuffer,
T const *y, sd::LongType const *yShapeBuffer, T *z,
sd::LongType const *zShapeBuffer, T *extraArguments) {
// ... (CUDA implementation remains unchanged)
}
#endif
static inline void specialOp(sd::Pointer state, const T *x, const sd::LongType *xShapeBuffer, const T *y,
const sd::LongType *yShapeBuffer, T *z, const sd::LongType *zShapeBuffer,
T *extraArguments) {
int trials = (int)extraArguments[0];
sd::LongType zLength = shape::length(zShapeBuffer);
int elementsPerThread = zLength / TAD_THRESHOLD;
int _threads = sd::math::sd_max<int>(1, elementsPerThread);
_threads = sd::math::sd_min<int>(_threads, sd::Environment::getInstance().maxThreads());
sd::LongType zRank = shape::rank(zShapeBuffer);
sd::LongType *zShape = shape::shapeOf(zShapeBuffer);
sd::LongType *zStride = shape::stride(zShapeBuffer);
sd::LongType yRank = shape::rank(yShapeBuffer);
sd::LongType *yShape = shape::shapeOf(yShapeBuffer);
sd::LongType *yStride = shape::stride(yShapeBuffer);
sd::LongType *xShape = shape::shapeOf(xShapeBuffer);
sd::LongType xRank = shape::rank(xShapeBuffer);
sd::LongType *xStride = shape::stride(xShapeBuffer);
T prob = extraArguments[1];
auto rng = reinterpret_cast<sd::graph::RandomGenerator *>(state);
auto func = PRAGMA_THREADS_FOR {
for (auto e = start; e < stop; e++) {
sd::LongType coords[SD_MAX_RANK];
INDEX2COORDS(e,zRank, zShape, coords);
sd::LongType zOffset;
COORDS2INDEX(zRank, zStride, coords, zOffset);
int success = 0;
for (int t = 1; t <= trials; t++) {
T randVal = rng->relativeT<T>((e + 1) * t);
if (y != z) {
// we're using external probs
sd::LongType yOffset;
COORDS2INDEX(shape::rank(yShapeBuffer), shape::stride(yShapeBuffer), coords, yOffset);
prob = y[yOffset];
}
if (randVal < prob) success++;
}
// if trials is set to 0, effectively we just have successful memset
z[zOffset] = static_cast<T>(success);
}
};
samediff::Threads::parallel_for(func, 0, zLength, 1, _threads);
}
};
//////////////////////////////////////////////////////////////////////
template <typename T>
class TruncatedNormalDistribution {
private:
static SD_INLINE SD_HOST_DEVICE T step(sd::graph::RandomGenerator *rng, T mean, T stddev, sd::LongType e,
sd::LongType middle, T &z) {
auto epm = e + middle;
const T two_pi = static_cast<T>(2.0f) * static_cast<T>(3.14159265358979323846);
const T epsilon = static_cast<T>(1.e-5f);
// we need to get random values
T r0 = rng->relativeT<T>(e, epsilon, static_cast<T>(1.0f));
T r1 = rng->relativeT<T>(epm, epsilon, static_cast<T>(1.0f));
T realMean0 = mean;
auto z0 = (sd::math::sd_sqrt<T, T>(static_cast<T>(-2.0f) * sd::math::sd_log<T, T>(r0)) *
sd::math::sd_cos<T, T>(two_pi * r1)) *
stddev +
realMean0;
z = z0;
if (epm < middle) {
T realMean1 = mean;
auto z1 = (sd::math::sd_sqrt<T, T>(static_cast<T>(-2.0f) * sd::math::sd_log<T, T>(r0)) *
sd::math::sd_sin<T, T>(two_pi * r1)) *
stddev +
realMean1;
z = z1;
}
return z;
}
public:
method_XY method_X method_idx
static const bool requiresSpecial = true;
#ifdef __CUDACC__
static SD_INLINE SD_DEVICE void specialOpCuda(sd::Pointer state, T const *x, sd::LongType const *xShapeBuffer,
T const *y, sd::LongType const *yShapeBuffer, T *z,
sd::LongType const *zShapeBuffer, T *extraArguments) {
// ... (CUDA implementation remains unchanged)
}
#endif
static inline void specialOp(sd::Pointer state, const T *x, const sd::LongType *xShapeBuffer, const T *y,
const sd::LongType *yShapeBuffer, T *z, const sd::LongType *zShapeBuffer,
T *extraArguments) {
GaussianDistribution<T>::specialOp(state, x, xShapeBuffer, y, yShapeBuffer, z, zShapeBuffer, extraArguments);
sd::LongType zLength = shape::length(zShapeBuffer);
auto rng = reinterpret_cast<sd::graph::RandomGenerator *>(state);
T mean = extraArguments[0];
T stddev = extraArguments[1];
T ds = sd::math::sd_abs<T,T>(stddev) * (T)2.0f;
sd::LongType middle = zLength / 2 + (zLength % 2);
int elementsPerThread = middle / TAD_THRESHOLD;
int _threads = sd::math::sd_max<int>(1, elementsPerThread);
_threads = sd::math::sd_min<int>(_threads, sd::Environment::getInstance().maxThreads());
sd::LongType zRank = shape::rank(zShapeBuffer);
sd::LongType *zShape = shape::shapeOf(zShapeBuffer);
sd::LongType *zStride = shape::stride(zShapeBuffer);
sd::LongType yRank = shape::rank(yShapeBuffer);
sd::LongType *yShape = shape::shapeOf(yShapeBuffer);
sd::LongType *yStride = shape::stride(yShapeBuffer);
sd::LongType *xShape = shape::shapeOf(xShapeBuffer);
sd::LongType xRank = shape::rank(xShapeBuffer);
sd::LongType *xStride = shape::stride(xShapeBuffer);
auto func = PRAGMA_THREADS_FOR {
for (auto e = start; e < stop; e++) {
sd::LongType coords[SD_MAX_RANK];
INDEX2COORDS(e, zRank,zShape, coords);
sd::LongType zOffset;
COORDS2INDEX(zRank, zStride, coords, zOffset);
if (z[zOffset] > mean + ds || z[zOffset] < mean - ds) {
z[zOffset] = step(rng, mean, stddev, e, middle, z[zOffset]);
if (z[zOffset] > mean + ds || z[zOffset] < mean - ds) z[zOffset] = mean + sd::DataTypeUtils::min_positive<T>();
}
}
};
samediff::Threads::parallel_for(func, 0, zLength, 1, _threads);
}
};
//////////////////////////////////////////////////////////////////////
template <typename T>
class LogNormalDistribution {
public:
method_XY method_X method_idx
static const bool requiresSpecial = true;
#ifdef __CUDACC__
static SD_INLINE SD_DEVICE void specialOpCuda(sd::Pointer state, T const *x, sd::LongType const *xShapeBuffer,
T const *y, sd::LongType const *yShapeBuffer, T *z,
sd::LongType const *zShapeBuffer, T *extraArguments) {
__shared__ T epsilon;
__shared__ T two_pi;
__shared__ sd::LongType zLength;
__shared__ T mean;
__shared__ T stddev;
__shared__ int step;
__shared__ T *tZ;
__shared__ sd::graph::RandomGenerator *rng;
__shared__ unsigned char *cB;
__shared__ unsigned char *dB;
__shared__ sd::graph::RandomGenerator *devRng;
__shared__ sd::LongType yRank;
__shared__ sd::LongType *yShape;
__shared__ sd::LongType *yStride;
__shared__ sd::LongType *xShape;
__shared__ sd::LongType xRank;
__shared__ sd::LongType *xStride;
__shared__ sd::LongType *zShape;
__shared__ sd::LongType *zStride;
__shared__ sd::LongType zRank;
if (threadIdx.x == 0) {
extern __shared__ unsigned char shmem[];
rng = reinterpret_cast<sd::graph::RandomGenerator *>(state);
cB = shmem;
devRng = reinterpret_cast<sd::graph::RandomGenerator *>(state);
dB = reinterpret_cast<unsigned char *>(state);
tZ = reinterpret_cast<T *>(shmem + sizeof(sd::graph::RandomGenerator));
zLength = shape::length(zShapeBuffer);
epsilon = static_cast<T>(1e-5);
two_pi = static_cast<T>(2.0f) * static_cast<T>(3.14159265358979323846);
mean = extraArguments[0];
stddev = extraArguments[1];
step = (blockDim.x * gridDim.x);
xRank = shape::rank(xShapeBuffer);
xShape = shape::shapeOf(xShapeBuffer);
xStride = shape::stride(xShapeBuffer);
yRank = shape::rank(yShapeBuffer);
yShape = shape::shapeOf(yShapeBuffer);
yStride = shape::stride(yShapeBuffer);
zRank = shape::rank(zShapeBuffer);
zShape = shape::shapeOf(zShapeBuffer);
zStride = shape::stride(zShapeBuffer);
}
__syncthreads();
// using this loop instead of memcpy
for (int e = threadIdx.x; e < sizeof(sd::graph::RandomGenerator); e += blockDim.x) cB[e] = dB[e];
__syncthreads();
int tid = blockIdx.x * blockDim.x + threadIdx.x;
int middle = zLength % 2 == 0 ? zLength / 2 : zLength / 2 + 1;
for (sd::LongType e = tid; e < middle; e += step) {
auto epm = e + middle;
// we need to get random values
T r0 = rng->relativeT<T>(e, epsilon, static_cast<T>(1.0f));
T r1 = rng->relativeT<T>(epm, epsilon, static_cast<T>(1.0f));
sd::LongType coords[SD_MAX_RANK];
sd::LongType yCoords[SD_MAX_RANK];
INDEX2COORDS(e, yRank, yShape, yCoords);
sd::LongType yOffset;
COORDS2INDEX(yRank, yStride, yCoords, yOffset);
sd::LongType zOffset;
INDEX2COORDS(e, zRank, zShape, coords);
COORDS2INDEX(zRank, zStride, coords, zOffset);
T realMean = y == z ? mean : y[yOffset];
z[zOffset] =
sd::math::sd_exp<T, T>((sd::math::sd_sqrt<T, T>(static_cast<T>(-2.0f) * sd::math::sd_log<T, T>(r0)) *
sd::math::sd_cos<T, T>(two_pi * r1)) *
stddev +
realMean);
if (epm < zLength) {
realMean = y == z ? mean : y[epm + yOffset];
z[epm + zOffset] =
sd::math::sd_exp<T, T>((sd::math::sd_sqrt<T, T>(static_cast<T>(-2.0f) * sd::math::sd_log<T, T>(r0)) *
sd::math::sd_sin<T, T>(two_pi * r1)) *
stddev +
realMean);
}
}
}
#endif
static inline void specialOp(sd::Pointer state, const T *x, const sd::LongType *xShapeBuffer, const T *y,
const sd::LongType *yShapeBuffer, T *z, const sd::LongType *zShapeBuffer,
T *extraArguments) {
const T two_pi = static_cast<T>(2.0f) * static_cast<T>(3.14159265358979323846);
sd::LongType zLength = shape::length(zShapeBuffer);
auto middle = zLength % 2 == 0 ? zLength / 2 : zLength / 2 + 1;
int elementsPerThread = middle / TAD_THRESHOLD;
int _threads = sd::math::sd_max<int>(1, elementsPerThread);
_threads = sd::math::sd_min<int>(_threads, sd::Environment::getInstance().maxThreads());
auto rng = reinterpret_cast<sd::graph::RandomGenerator *>(state);
const T mean = extraArguments[0];
const T stddev = extraArguments[1];
const T epsilon = static_cast<T>(1e-5);
sd::LongType zRank = shape::rank(zShapeBuffer);
sd::LongType *zShape = shape::shapeOf(zShapeBuffer);
sd::LongType *zStride = shape::stride(zShapeBuffer);
sd::LongType yRank = shape::rank(yShapeBuffer);
sd::LongType *yShape = shape::shapeOf(yShapeBuffer);
sd::LongType *yStride = shape::stride(yShapeBuffer);
sd::LongType *xShape = shape::shapeOf(xShapeBuffer);
sd::LongType xRank = shape::rank(xShapeBuffer);
sd::LongType *xStride = shape::stride(xShapeBuffer);
auto func = PRAGMA_THREADS_FOR {
for (auto e = start; e < stop; e++) {
sd::LongType coords[SD_MAX_RANK];
INDEX2COORDS(e, zRank, zShape, coords);
sd::LongType zOffset;
COORDS2INDEX(zRank, zStride, coords, zOffset);
auto epm = e + middle;
// we need to get random values
T r0 = rng->relativeT<T>(e, epsilon, static_cast<T>(1.0f));
T r1 = rng->relativeT<T>(epm, epsilon, static_cast<T>(1.0f));
sd::LongType yOffset;
COORDS2INDEX(yRank, yStride, coords, yOffset);
T realMean = y == z ? mean : y[yOffset];
z[zOffset] =
sd::math::sd_exp<T, T>((sd::math::sd_sqrt<T, T>(static_cast<T>(-2.0f) * sd::math::sd_log<T, T>(r0)) *
sd::math::sd_cos<T, T>(two_pi * r1)) *
stddev +
realMean);
if (epm < zLength) {
INDEX2COORDS(epm,zRank, zShape, coords);
COORDS2INDEX(zRank, zStride, coords, zOffset);
COORDS2INDEX(yRank, yStride, coords, yOffset);
realMean = y == z ? mean : y[yOffset];
z[zOffset] =
sd::math::sd_exp<T, T>((sd::math::sd_sqrt<T, T>(static_cast<T>(-2.0f) * sd::math::sd_log<T, T>(r0)) *
sd::math::sd_sin<T, T>(two_pi * r1)) *
stddev +
realMean);
}
}
};
samediff::Threads::parallel_for(func, 0, middle, 1, _threads);
}
};
} // namespace randomOps
#endif // LIBND4J_SPECIAL_RANDOM_OPS_H