Gamma and Poisson distributions (#27)
* Added implementation for random_gamma op. * Added implementation for random_poisson op and support classes. * Added helpers for random_poisson and random_gamma ops. * Implementation of random_poisson. The first working edition. * Implementation of random_poisson. Parallelized working edition. * Implementation of random_gamma. Parallelized working edition with alpha only. * Added cuda implementation for helper of poisson distribution. * Corrected shape calculation with random_gamma and tests. * Finished cpu implementation for gamma distribution. * Finished cuda implementation for random_gamma op. * Refactored cpu helpers for random_gamma and random_poisson ops. * Refactored cuda helpers for gamma and poisson distribution. * Refactored cuda helper for gamma distribution. * Refactored cpu helper for random_poisson op. * Refactored cpu helper for random_gamma op.master
shugeo
GitHub
parent
948ebef41c
commit
7b14a9f603
|
|
@ -323,7 +323,9 @@
|
|||
(11, TruncatedNormalDistribution) ,\
|
||||
(12, AlphaDropOut),\
|
||||
(13, ExponentialDistribution),\
|
||||
(14, ExponentialDistributionInv)
|
||||
(14, ExponentialDistributionInv), \
|
||||
(15, PoissonDistribution), \
|
||||
(16, GammaDistribution)
|
||||
|
||||
#define PAIRWISE_INT_OPS \
|
||||
(0, ShiftLeft), \
|
||||
|
|
|
|||
|
|
@ -0,0 +1,83 @@
|
|||
/*******************************************************************************
|
||||
* Copyright (c) 2015-2018 Skymind, Inc.
|
||||
*
|
||||
* 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.
|
||||
*
|
||||
* 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
|
||||
******************************************************************************/
|
||||
|
||||
//
|
||||
// @author George A. Shulinok <sgazeos@gmail.com>
|
||||
//
|
||||
|
||||
#include <op_boilerplate.h>
|
||||
#if NOT_EXCLUDED(OP_random_gamma)
|
||||
|
||||
#include <ops/declarable/headers/random.h>
|
||||
#include <ops/declarable/helpers/random.h>
|
||||
|
||||
namespace nd4j {
|
||||
namespace ops {
|
||||
CUSTOM_OP_IMPL(random_gamma, 2, 1, false, 0, 0) {
|
||||
// gamma distribution
|
||||
auto rng = block.randomGenerator();
|
||||
auto shape = INPUT_VARIABLE(0);
|
||||
auto alpha = INPUT_VARIABLE(1);
|
||||
NDArray* beta = nullptr;
|
||||
|
||||
if (block.width() > 2) {
|
||||
beta = INPUT_VARIABLE(2);
|
||||
REQUIRE_TRUE(ShapeUtils::areShapesBroadcastable(*alpha, *beta), 0, "random_gamma: alpha and beta shapes should be broadcastable.");
|
||||
}
|
||||
|
||||
auto output = OUTPUT_VARIABLE(0);
|
||||
auto seed = 0;
|
||||
|
||||
if (block.getIArguments()->size()) {
|
||||
seed = INT_ARG(0);
|
||||
}
|
||||
|
||||
rng.setSeed(seed);
|
||||
|
||||
helpers::fillRandomGamma(block.launchContext(), rng, alpha, beta, output);
|
||||
|
||||
return Status::OK();
|
||||
}
|
||||
|
||||
DECLARE_SHAPE_FN(random_gamma) {
|
||||
auto in = INPUT_VARIABLE(0);
|
||||
auto shape = in->template asVectorT<Nd4jLong>();
|
||||
auto alphaShape = inputShape->at(1);
|
||||
auto additionalShape = alphaShape;
|
||||
if (inputShape->size() > 2) {
|
||||
auto rest = inputShape->at(2); additionalShape = nullptr;
|
||||
REQUIRE_TRUE(ShapeUtils::areShapesBroadcastable(alphaShape, rest), 0, "random_gamma: alpha and beta shapes should be broadcastable.");
|
||||
ShapeUtils::evalBroadcastShapeInfo(alphaShape, rest, true, additionalShape, block.workspace());
|
||||
}
|
||||
auto lastDim = shape::sizeAt(alphaShape, 0);
|
||||
auto dtype = ArrayOptions::dataType(alphaShape);
|
||||
for (auto i = 0; i < shape::rank(additionalShape); i++)
|
||||
shape.push_back(shape::sizeAt(additionalShape, i));
|
||||
auto newShape = ConstantShapeHelper::getInstance()->createShapeInfo(dtype, 'c', shape);
|
||||
return SHAPELIST(newShape);
|
||||
}
|
||||
|
||||
DECLARE_TYPES(random_gamma) {
|
||||
getOpDescriptor()
|
||||
->setAllowedInputTypes(0, {ALL_INTS})
|
||||
->setAllowedInputTypes(1, {ALL_FLOATS})
|
||||
->setAllowedInputTypes(2, {ALL_FLOATS})
|
||||
->setAllowedOutputTypes({ALL_FLOATS});
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,67 @@
|
|||
/*******************************************************************************
|
||||
* Copyright (c) 2015-2018 Skymind, Inc.
|
||||
*
|
||||
* 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.
|
||||
*
|
||||
* 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
|
||||
******************************************************************************/
|
||||
|
||||
//
|
||||
// @author George A. Shulinok <sgazeos@gmail.com>
|
||||
//
|
||||
|
||||
#include <op_boilerplate.h>
|
||||
#if NOT_EXCLUDED(OP_random_poisson)
|
||||
|
||||
#include <ops/declarable/headers/random.h>
|
||||
#include <ops/declarable/helpers/random.h>
|
||||
|
||||
namespace nd4j {
|
||||
namespace ops {
|
||||
CUSTOM_OP_IMPL(random_poisson, 2, 1, false, 0, 0) {
|
||||
// gamma distribution
|
||||
auto rng = block.randomGenerator();
|
||||
auto shape = INPUT_VARIABLE(0);
|
||||
auto lambda = INPUT_VARIABLE(1);
|
||||
auto output = OUTPUT_VARIABLE(0);
|
||||
auto seed = 0;
|
||||
if (block.getIArguments()->size()) {
|
||||
seed = INT_ARG(0);
|
||||
}
|
||||
rng.setSeed(seed);
|
||||
helpers::fillRandomPoisson(block.launchContext(), rng, lambda, output);
|
||||
|
||||
return Status::OK();
|
||||
}
|
||||
|
||||
|
||||
DECLARE_SHAPE_FN(random_poisson) {
|
||||
auto in = INPUT_VARIABLE(0);
|
||||
auto shape = in->template asVectorT<Nd4jLong>();
|
||||
auto lambdaShape = inputShape->at(1);
|
||||
auto dtype = ArrayOptions::dataType(lambdaShape);
|
||||
for (auto d = 0; d < shape::rank(lambdaShape); ++d ) {
|
||||
shape.emplace_back(shape::sizeAt(lambdaShape, d));
|
||||
}
|
||||
auto newShape = ConstantShapeHelper::getInstance()->createShapeInfo(dtype, 'c', shape);
|
||||
return SHAPELIST(newShape);
|
||||
}
|
||||
|
||||
DECLARE_TYPES(random_poisson) {
|
||||
getOpDescriptor()
|
||||
->setAllowedInputTypes(0, {ALL_INTS})
|
||||
->setAllowedInputTypes(1, {ALL_FLOATS})
|
||||
->setAllowedOutputTypes({ALL_FLOATS});
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
|
|
@ -49,7 +49,23 @@ namespace nd4j {
|
|||
DECLARE_CUSTOM_OP(random_exponential, 1, 1, true, 1, 0);
|
||||
#endif
|
||||
|
||||
#if NOT_EXCLUDED(OP_random_crop)
|
||||
DECLARE_CUSTOM_OP(random_crop, 2, 1, false, 0, 0);
|
||||
#endif
|
||||
|
||||
/**
|
||||
* random_gamma op.
|
||||
*/
|
||||
#if NOT_EXCLUDED(OP_random_gamma)
|
||||
DECLARE_CUSTOM_OP(random_gamma, 2, 1, false, 0, 0);
|
||||
#endif
|
||||
|
||||
/**
|
||||
* random_poisson op.
|
||||
*/
|
||||
#if NOT_EXCLUDED(OP_random_poisson)
|
||||
DECLARE_CUSTOM_OP(random_poisson, 2, 1, false, 0, 0);
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,132 @@
|
|||
/*******************************************************************************
|
||||
* Copyright (c) 2015-2018 Skymind, Inc.
|
||||
*
|
||||
* 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.
|
||||
*
|
||||
* 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
|
||||
******************************************************************************/
|
||||
|
||||
//
|
||||
// @author sgazeos@gmail.com
|
||||
//
|
||||
|
||||
#include <ops/declarable/helpers/random.h>
|
||||
//#include <vector>
|
||||
#include <memory>
|
||||
//#include <graph/Context.h>
|
||||
#include <ShapeUtils.h>
|
||||
|
||||
namespace nd4j {
|
||||
namespace ops {
|
||||
namespace helpers {
|
||||
|
||||
template <typename T>
|
||||
void fillRandomGamma_(LaunchContext* context, graph::RandomGenerator& rng, NDArray* alpha, NDArray* beta, NDArray* output) {
|
||||
|
||||
Nd4jLong* broadcasted = nullptr;
|
||||
if (beta != nullptr)
|
||||
ShapeUtils::evalBroadcastShapeInfo(*alpha, *beta, true, broadcasted, context->getWorkspace());
|
||||
else
|
||||
broadcasted = alpha->shapeInfo();
|
||||
auto step = shape::length(broadcasted);
|
||||
auto shift = output->lengthOf() / step;
|
||||
|
||||
auto copyAlpha = alpha;
|
||||
auto copyBeta = beta;
|
||||
if (beta != nullptr) {
|
||||
NDArray alphaBroadcasted(broadcasted, alpha->dataType(), false, context);
|
||||
NDArray betaBroadcasted(broadcasted, beta->dataType(), false, context);
|
||||
|
||||
copyAlpha = (alphaBroadcasted.applyTrueBroadcast(BroadcastOpsTuple::Assign(), alpha));
|
||||
copyBeta = (betaBroadcasted.applyTrueBroadcast(BroadcastOpsTuple::Assign(), beta));
|
||||
|
||||
}
|
||||
// bool directAlpha = alpha->ews() == 1 && alpha->ordering() == 'c';
|
||||
bool directOutput = output->ews() == 1 && output->ordering() == 'c';
|
||||
T* outputBuf = output->dataBuffer()->primaryAsT<T>();
|
||||
|
||||
PRAGMA_OMP_PARALLEL_FOR
|
||||
for (auto k = 0; k < shift; k++) {
|
||||
auto pos = k * step;
|
||||
auto u = rng.relativeT<T>(k, 0., 1.);
|
||||
for (auto e = 0; e < step; e++)
|
||||
if (directOutput) {
|
||||
outputBuf[pos + e] = math::nd4j_igamma<T, T, T>(copyAlpha->t<T>(e),
|
||||
beta != nullptr ? copyBeta->t<T>(e) * u : u);
|
||||
}
|
||||
else {
|
||||
output->t<T>(pos + e) = math::nd4j_igamma<T, T, T>(copyAlpha->t<T>(e),
|
||||
beta != nullptr ? copyBeta->t<T>(e) * u : u);
|
||||
}
|
||||
}
|
||||
|
||||
if (beta != nullptr) {
|
||||
delete copyAlpha;
|
||||
delete copyBeta;
|
||||
//delete broadcasted;
|
||||
}
|
||||
}
|
||||
|
||||
void fillRandomGamma(LaunchContext* context, graph::RandomGenerator& rng, NDArray* alpha, NDArray* beta, NDArray* output) {
|
||||
BUILD_SINGLE_SELECTOR(output->dataType(), fillRandomGamma_, (context, rng, alpha, beta, output), FLOAT_NATIVE);
|
||||
}
|
||||
BUILD_SINGLE_TEMPLATE(template void fillRandomGamma_, (LaunchContext* context,
|
||||
graph::RandomGenerator& rng, NDArray* alpha, NDArray* beta, NDArray* output), FLOAT_NATIVE);
|
||||
|
||||
/*
|
||||
* algorithm Poisson generator based upon the inversion by sequential search:[48]:505
|
||||
init:
|
||||
Let x ← 0, p ← e−λ, s ← p.
|
||||
Generate uniform random number u in [0,1].
|
||||
while u > s do:
|
||||
x ← x + 1.
|
||||
p ← p * λ / x.
|
||||
s ← s + p.
|
||||
return x.
|
||||
* */
|
||||
template <typename T>
|
||||
void fillRandomPoisson_(LaunchContext* context, graph::RandomGenerator& rng, NDArray* lambda, NDArray* output) {
|
||||
auto shift = output->lengthOf() / lambda->lengthOf();
|
||||
auto step = lambda->lengthOf();
|
||||
T* lambdaBuf = lambda->dataBuffer()->primaryAsT<T>();
|
||||
T* outputBuf = output->dataBuffer()->primaryAsT<T>();
|
||||
bool directLa = lambda->ews() == 1 && lambda->ordering() == 'c';
|
||||
bool directOut = output->ews() == 1 && output->ordering() == 'c';
|
||||
PRAGMA_OMP_PARALLEL_FOR
|
||||
for (auto k = 0; k < shift; k++) {
|
||||
auto pos = k * step;
|
||||
auto u = rng.relativeT<T>(k, 0., 1.);
|
||||
for (auto e = 0; e < step; e++) {
|
||||
auto p = math::nd4j_exp<T, T>(-lambda->t<T>(e));
|
||||
auto s = p;
|
||||
auto x = T(0.f);
|
||||
while (u > s) {
|
||||
x += 1.f;
|
||||
p *= directLa?lambdaBuf[e]/x:lambda->t<T>(e) / x;
|
||||
s += p;
|
||||
}
|
||||
if (directOut)
|
||||
outputBuf[pos + e] = x;
|
||||
else
|
||||
output->t<T>(pos + e) = x;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void fillRandomPoisson(LaunchContext* context, graph::RandomGenerator& rng, NDArray* lambda, NDArray* output) {
|
||||
BUILD_SINGLE_SELECTOR(output->dataType(), fillRandomPoisson_, (context, rng, lambda, output), FLOAT_NATIVE);
|
||||
}
|
||||
BUILD_SINGLE_TEMPLATE(template void fillRandomPoisson_, (LaunchContext* context,
|
||||
graph::RandomGenerator& rng, NDArray* lambda, NDArray* output), FLOAT_TYPES);
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,186 @@
|
|||
/*******************************************************************************
|
||||
* Copyright (c) 2015-2018 Skymind, Inc.
|
||||
*
|
||||
* 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.
|
||||
*
|
||||
* 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
|
||||
******************************************************************************/
|
||||
|
||||
//
|
||||
// @author sgazeos@gmail.com
|
||||
//
|
||||
|
||||
#include <ops/declarable/helpers/random.h>
|
||||
//#include <NativeOps.h>
|
||||
#include <vector>
|
||||
#include <memory>
|
||||
#include <graph/Context.h>
|
||||
#include <helpers/RandomLauncher.h>
|
||||
#include <ShapeUtils.h>
|
||||
#include <NDArrayFactory.h>
|
||||
|
||||
namespace nd4j {
|
||||
namespace ops {
|
||||
namespace helpers {
|
||||
|
||||
/*
|
||||
* fillGammaKernel - fill up output with gamma distributed values
|
||||
*
|
||||
* uList - uniformly distributed values set
|
||||
* uLength - length of uList
|
||||
* alpha - alpha param
|
||||
* beta - beta param
|
||||
* output - distributed output.
|
||||
* */
|
||||
template <typename T>
|
||||
static __global__ void fillGammaKernel(T* uList, Nd4jLong uLength, T* alpha, Nd4jLong* alphaShape,
|
||||
T* beta, Nd4jLong* betaShape, T* output, Nd4jLong* outputShape) {
|
||||
// fill up
|
||||
__shared__ Nd4jLong aLength;
|
||||
if (threadIdx.x == 0) {
|
||||
aLength = shape::length(alphaShape);
|
||||
}
|
||||
__syncthreads();
|
||||
|
||||
for (auto k = blockIdx.x; k < (int)uLength; k += gridDim.x) {
|
||||
auto pos = k * aLength;
|
||||
auto u = uList[k]; // this is a vector
|
||||
for (auto e = threadIdx.x; e < (int)aLength; e += blockDim.x) {
|
||||
auto aIndex = shape::getIndexOffset(e, alphaShape);
|
||||
auto bIndex = betaShape?shape::getIndexOffset(e, betaShape):-1LL;
|
||||
auto betaV = T(beta != nullptr ? beta[bIndex] * u : u);
|
||||
auto zIndex = shape::getIndexOffset(e + pos, outputShape);
|
||||
|
||||
output[zIndex] = math::nd4j_igamma<T, T, T>(alpha[aIndex], betaV);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
static void fillRandomGamma_(LaunchContext* context, graph::RandomGenerator& rng, NDArray* alpha, NDArray* beta, NDArray* output) {
|
||||
// To fill up output need to broadcast alpha and beta to the same shape and in
|
||||
Nd4jLong* broadcasted = nullptr;
|
||||
if (beta != nullptr)
|
||||
ShapeUtils::evalBroadcastShapeInfo(*alpha, *beta, true, broadcasted, context->getWorkspace());
|
||||
else
|
||||
broadcasted = alpha->shapeInfo();
|
||||
auto step = shape::length(broadcasted);
|
||||
auto shift = output->lengthOf() / step;
|
||||
|
||||
auto copyAlpha = alpha;
|
||||
auto copyBeta = beta;
|
||||
if (beta != nullptr) {
|
||||
NDArray alphaBroadcasted(broadcasted, alpha->dataType(), true, context);
|
||||
NDArray betaBroadcasted(broadcasted, beta->dataType(), true, context);
|
||||
|
||||
copyAlpha = (alphaBroadcasted.applyTrueBroadcast(BroadcastOpsTuple::Assign(), alpha));
|
||||
copyBeta = (betaBroadcasted.applyTrueBroadcast(BroadcastOpsTuple::Assign(), beta));
|
||||
copyAlpha->tickWriteDevice(); copyBeta->tickWriteDevice();
|
||||
}
|
||||
|
||||
auto stream = context->getCudaStream();
|
||||
NDArray uniform = NDArrayFactory::create<T>('c', {shift}, context);
|
||||
uniform.syncToDevice();
|
||||
// fill up uniform with given length
|
||||
RandomLauncher::fillUniform(context, rng, &uniform, 0., 1.);
|
||||
|
||||
fillGammaKernel<T><<<128, 128, 256, *stream>>>(uniform.dataBuffer()->specialAsT<T>(), shift,
|
||||
copyAlpha->dataBuffer()->specialAsT<T>(), copyAlpha->specialShapeInfo(),
|
||||
beta?copyBeta->dataBuffer()->specialAsT<T>():(T*)nullptr,
|
||||
beta?copyBeta->specialShapeInfo():(Nd4jLong*)nullptr,
|
||||
output->dataBuffer()->specialAsT<T>(), output->specialShapeInfo());
|
||||
|
||||
if (beta != nullptr) {
|
||||
delete copyAlpha;
|
||||
delete copyBeta;
|
||||
//delete broadcasted;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void fillRandomGamma(LaunchContext* context, graph::RandomGenerator& rng, NDArray* alpha, NDArray* beta, NDArray* output) {
|
||||
if (beta)
|
||||
NDArray::prepareSpecialUse({output}, {alpha, beta});
|
||||
else
|
||||
NDArray::prepareSpecialUse({output}, {alpha});
|
||||
BUILD_SINGLE_SELECTOR(output->dataType(), fillRandomGamma_, (context, rng, alpha, beta, output), FLOAT_NATIVE);
|
||||
if (beta)
|
||||
NDArray::registerSpecialUse({output}, {alpha, beta});
|
||||
else
|
||||
NDArray::prepareSpecialUse({output}, {alpha});
|
||||
}
|
||||
BUILD_SINGLE_TEMPLATE(template void fillRandomGamma_, (LaunchContext* context, graph::RandomGenerator& rng, NDArray* alpha, NDArray* beta, NDArray* output), FLOAT_NATIVE);
|
||||
|
||||
|
||||
/*
|
||||
* algorithm Poisson generator based upon the inversion by sequential search
|
||||
*
|
||||
init:
|
||||
Let x ← 0, p ← e−λ, s ← p.
|
||||
using uniformly random sequence U (u in U) distributed at [0, 1].
|
||||
while u > s do:
|
||||
x ← x + 1.
|
||||
p ← p * λ / x.
|
||||
s ← s + p.
|
||||
return x.
|
||||
* */
|
||||
template <typename T>
|
||||
static __global__ void fillPoissonKernel(T* uList, Nd4jLong uLength, T* lambda, Nd4jLong* lambdaShape, T* output,
|
||||
Nd4jLong* outputShape) {
|
||||
|
||||
__shared__ Nd4jLong step;
|
||||
|
||||
if (threadIdx.x == 0) {
|
||||
step = shape::length(lambdaShape);
|
||||
}
|
||||
__syncthreads();
|
||||
|
||||
for (auto k = blockIdx.x; k < (int)uLength; k += gridDim.x) {
|
||||
auto pos = k * step;
|
||||
auto u = uList[k];
|
||||
for (auto e = threadIdx.x; e < step; e += blockDim.x) {
|
||||
auto p = math::nd4j_exp<T,T>(-lambda[e]);
|
||||
auto s = p;
|
||||
auto x = T(0.f);
|
||||
auto lIndex = shape::getIndexOffset(e, lambdaShape);
|
||||
auto zIndex = shape::getIndexOffset(e + pos, outputShape);
|
||||
while (u > s) {
|
||||
x += T(1.);
|
||||
p *= lambda[lIndex] / x;
|
||||
s += p;
|
||||
}
|
||||
output[zIndex] = x;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
static void fillRandomPoisson_(LaunchContext* context, graph::RandomGenerator& rng, NDArray* lambda, NDArray* output) {
|
||||
auto shift = output->lengthOf() / lambda->lengthOf();
|
||||
NDArray uniform('c', {shift}, output->dataType());
|
||||
auto stream = context->getCudaStream();
|
||||
// fill up uniform with given length
|
||||
RandomLauncher::fillUniform(context, rng, &uniform, 0., 1.);
|
||||
fillPoissonKernel<T><<<128, 256, 128, *stream>>>(uniform.dataBuffer()->specialAsT<T>(), uniform.lengthOf(),
|
||||
lambda->dataBuffer()->specialAsT<T>(), lambda->specialShapeInfo(),
|
||||
output->dataBuffer()->specialAsT<T>(), output->specialShapeInfo());
|
||||
}
|
||||
|
||||
void fillRandomPoisson(LaunchContext* context, graph::RandomGenerator& rng, NDArray* lambda, NDArray* output) {
|
||||
NDArray::prepareSpecialUse({output}, {lambda});
|
||||
BUILD_SINGLE_SELECTOR(output->dataType(), fillRandomPoisson_, (context, rng, lambda, output), FLOAT_NATIVE);
|
||||
NDArray::registerSpecialUse({output}, {lambda});
|
||||
}
|
||||
|
||||
BUILD_SINGLE_TEMPLATE(template void fillRandomPoisson_, (LaunchContext* context, graph::RandomGenerator& rng, NDArray* lambda, NDArray* output), FLOAT_NATIVE);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,40 @@
|
|||
/*******************************************************************************
|
||||
* Copyright (c) 2015-2018 Skymind, Inc.
|
||||
*
|
||||
* 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.
|
||||
*
|
||||
* 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
|
||||
******************************************************************************/
|
||||
|
||||
//
|
||||
// @author sgazeos@gmail.com
|
||||
//
|
||||
//
|
||||
// Declaration of distribution helpers
|
||||
//
|
||||
#ifndef __RANDOM_HELPERS__
|
||||
#define __RANDOM_HELPERS__
|
||||
#include <op_boilerplate.h>
|
||||
#include <NDArray.h>
|
||||
#include <helpers/helper_random.h>
|
||||
#include <graph/Context.h>
|
||||
|
||||
namespace nd4j {
|
||||
namespace ops {
|
||||
namespace helpers {
|
||||
|
||||
void fillRandomGamma(LaunchContext* context, graph::RandomGenerator& rng, NDArray* alpha, NDArray* beta, NDArray* output);
|
||||
void fillRandomPoisson(LaunchContext* context, graph::RandomGenerator& rng, NDArray* lambda, NDArray* output);
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
@ -129,6 +129,47 @@ namespace randomOps {
|
|||
}
|
||||
};
|
||||
|
||||
template <typename T>
|
||||
class PoissonDistribution {
|
||||
public:
|
||||
no_exec_special
|
||||
no_exec_special_cuda
|
||||
|
||||
method_XY
|
||||
|
||||
random_def T op(Nd4jLong idx, Nd4jLong length, nd4j::graph::RandomGenerator *helper, T *extraParams) {
|
||||
T lambda = extraParams[0];
|
||||
T x = helper->relativeT(idx, -nd4j::DataTypeUtils::template max<T>() / 10 , nd4j::DataTypeUtils::template max<T>() / 10);
|
||||
return x <= (T)0.f ? (T)0.f : nd4j::math::nd4j_igammac<T,T,T>(nd4j::math::nd4j_floor<T,T>(x), lambda);
|
||||
}
|
||||
|
||||
random_def T op(T valueX, Nd4jLong idx, Nd4jLong length, nd4j::graph::RandomGenerator *helper, T *extraParams) {
|
||||
T lambda = extraParams[0];
|
||||
return valueX <= (T)0.f ? (T)0.f : (T)nd4j::math::nd4j_igammac<T,T,T>(nd4j::math::nd4j_floor<T,T>(valueX), lambda);
|
||||
}
|
||||
};
|
||||
|
||||
template <typename T>
|
||||
class GammaDistribution {
|
||||
public:
|
||||
no_exec_special
|
||||
no_exec_special_cuda
|
||||
|
||||
method_XY
|
||||
|
||||
random_def T op(Nd4jLong idx, Nd4jLong length, nd4j::graph::RandomGenerator *helper, T *extraParams) {
|
||||
T alpha = extraParams[0];
|
||||
T beta = extraParams[1];
|
||||
T x = helper->relativeT(idx, -nd4j::DataTypeUtils::template max<T>() / 10 , nd4j::DataTypeUtils::template max<T>() / 10);
|
||||
return x <= (T)0.f ? (T)0.f : nd4j::math::nd4j_igamma<T,T,T>(alpha, x * beta);
|
||||
}
|
||||
|
||||
random_def T op(T valueX, Nd4jLong idx, Nd4jLong length, nd4j::graph::RandomGenerator *helper, T *extraParams) {
|
||||
T alpha = extraParams[0];
|
||||
T beta = extraParams[1];
|
||||
return valueX <= (T)0.f ? (T)0.f : nd4j::math::nd4j_igamma<T,T,T>(alpha, beta * valueX);
|
||||
}
|
||||
};
|
||||
|
||||
/**
|
||||
* Basic DropOut/DropConnect Op
|
||||
|
|
|
|||
|
|
@ -894,6 +894,10 @@ namespace nd4j {
|
|||
Z aim = nd4j_pow<X, X, Z>(x, a) / (nd4j_exp<X, Z>(x) * nd4j_gamma<Y, Z>(a));
|
||||
auto sum = Z(0.);
|
||||
auto denom = Z(1.);
|
||||
if (a <= X(0.000001))
|
||||
//throw std::runtime_error("Cannot calculate gamma for a zero val.");
|
||||
return Z(0);
|
||||
|
||||
for (int i = 0; Z(1./denom) > Z(1.0e-12); i++) {
|
||||
denom *= (a + i);
|
||||
sum += nd4j_pow<X, int, Z>(x, i) / denom;
|
||||
|
|
|
|||
|
|
@ -773,6 +773,88 @@ TEST_F(RNGTests, Test_ExponentialDistribution_2) {
|
|||
delete result;
|
||||
}
|
||||
|
||||
TEST_F(RNGTests, Test_PoissonDistribution_1) {
|
||||
auto x = NDArrayFactory::create<Nd4jLong>('c', {1}, {10});
|
||||
auto la = NDArrayFactory::create<float>('c', {2, 3});
|
||||
auto exp0 = NDArrayFactory::create<float>('c', {10, 2, 3});
|
||||
|
||||
la.linspace(1.0);
|
||||
|
||||
|
||||
nd4j::ops::random_poisson op;
|
||||
auto result = op.execute({&x, &la}, {}, {});
|
||||
ASSERT_EQ(Status::OK(), result->status());
|
||||
|
||||
auto z = result->at(0);
|
||||
// z->printIndexedBuffer("Poisson distribution");
|
||||
ASSERT_TRUE(exp0.isSameShape(z));
|
||||
ASSERT_FALSE(exp0.equalsTo(z));
|
||||
|
||||
delete result;
|
||||
}
|
||||
|
||||
TEST_F(RNGTests, Test_GammaDistribution_1) {
|
||||
auto x = NDArrayFactory::create<Nd4jLong>('c', {1}, {10});
|
||||
auto al = NDArrayFactory::create<float>('c', {2, 3});
|
||||
auto exp0 = NDArrayFactory::create<float>('c', {10, 2, 3});
|
||||
|
||||
al.linspace(1.0);
|
||||
|
||||
|
||||
nd4j::ops::random_gamma op;
|
||||
auto result = op.execute({&x, &al}, {}, {});
|
||||
ASSERT_EQ(Status::OK(), result->status());
|
||||
|
||||
auto z = result->at(0);
|
||||
// z->printIndexedBuffer("Gamma distribution");
|
||||
ASSERT_TRUE(exp0.isSameShape(z));
|
||||
ASSERT_FALSE(exp0.equalsTo(z));
|
||||
|
||||
delete result;
|
||||
}
|
||||
|
||||
TEST_F(RNGTests, Test_GammaDistribution_2) {
|
||||
auto x = NDArrayFactory::create<Nd4jLong>('c', {1}, {10});
|
||||
auto al = NDArrayFactory::create<float>('c', {2, 3});
|
||||
auto be = NDArrayFactory::create<float>('c', {2, 3});
|
||||
auto exp0 = NDArrayFactory::create<float>('c', {10, 2, 3});
|
||||
|
||||
al.linspace(1.0);
|
||||
be.assign(1.0);
|
||||
|
||||
nd4j::ops::random_gamma op;
|
||||
auto result = op.execute({&x, &al, &be}, {}, {});
|
||||
ASSERT_EQ(Status::OK(), result->status());
|
||||
|
||||
auto z = result->at(0);
|
||||
// z->printIndexedBuffer("Gamma distribution");
|
||||
ASSERT_TRUE(exp0.isSameShape(z));
|
||||
ASSERT_FALSE(exp0.equalsTo(z));
|
||||
|
||||
delete result;
|
||||
}
|
||||
|
||||
TEST_F(RNGTests, Test_GammaDistribution_3) {
|
||||
auto x = NDArrayFactory::create<Nd4jLong>('c', {1}, {10});
|
||||
auto al = NDArrayFactory::create<float>('c', {3, 1});
|
||||
auto be = NDArrayFactory::create<float>('c', {1, 2});
|
||||
auto exp0 = NDArrayFactory::create<float>('c', {10, 3, 2});
|
||||
|
||||
al.linspace(1.0);
|
||||
be.assign(2.0);
|
||||
|
||||
nd4j::ops::random_gamma op;
|
||||
auto result = op.execute({&x, &al, &be}, {}, {});
|
||||
ASSERT_EQ(Status::OK(), result->status());
|
||||
|
||||
auto z = result->at(0);
|
||||
// z->printIndexedBuffer("Gamma distribution");
|
||||
ASSERT_TRUE(exp0.isSameShape(z));
|
||||
ASSERT_FALSE(exp0.equalsTo(z));
|
||||
|
||||
delete result;
|
||||
}
|
||||
|
||||
namespace nd4j {
|
||||
namespace tests {
|
||||
static void fillList(Nd4jLong seed, int numberOfArrays, std::vector<Nd4jLong> &shape, std::vector<NDArray*> &list, nd4j::graph::RandomGenerator *rng) {
|
||||
|
|
|
|||
Loading…
Reference in New Issue