Monte Carlo: Bumpy Sphere,A Tumor Model and Volume

Computer Science Level 2

Given the spherical equation of a (Bumpy Sphere)

ρ ( t , v ) = 1 + 0.2 ∗ s i n ( 6 ∗ t ) ∗ s i n ( 5 ∗ v ) , 0 < = t < = 2 ∗ π , 0 < = v < = π \rho(t,v) = 1 + 0.2*sin( 6*t )*sin( 5*v ) ,0 <=t<=2*\pi ,0< =v<=\pi

find its volume using the Monte Carlo method in conjunction with the GeoGebra app or any other graphical software.


The answer is 4.3.

Frank Petiprin by Frank Petiprin , 81, USA

This section requires Javascript.
You are seeing this because something didn't load right. We suggest you, (a) try refreshing the page, (b) enabling javascript if it is disabled on your browser and, finally, (c) loading the non-javascript version of this page . We're sorry about the hassle.

2 solutions

Piotr Idzik
May 22, 2020

Below you will find my C++ code using multithreading. 

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
 
#include <iostream>
#include <cmath>
#include <random>
#include <vector>
#include <chrono>
#include <thread>
#include <ctime>
#include <algorithm>

template<std::size_t Dimension, typename ValueType, typename RandomEngineType, typename IsInsideType>
class McWorker
{
    private:
        const IsInsideType isInside;
        const std::uintmax_t localIter;
        std::array<std::uniform_real_distribution<ValueType>, Dimension> pointGenerator;
        RandomEngineType randomEngine;
        ValueType amountInside;

    public:
        explicit McWorker(const std::array<std::pair<ValueType, ValueType>, Dimension>& _limitBox,
                          const IsInsideType& _isInside,
                          const std::uintmax_t _localIter,
                          const int _seed) :
        isInside(_isInside),
        localIter(_localIter),
        pointGenerator(),
        randomEngine(_seed),
        amountInside()
        {
            for(std::size_t curDim = 0; curDim < Dimension; ++curDim)
            {
                const auto curLimit = _limitBox[curDim];
                this->pointGenerator[curDim] = std::uniform_real_distribution<ValueType>(curLimit.first,
                                                                                         curLimit.second);
            }
        }
        void operator()()
        {
            for(std::uintmax_t curIter = 0; curIter < this->localIter; ++curIter)
            {
                const auto curPos = this->getRandomPoint();
                if(this->isInside(curPos))
                {
                    ++this->amountInside;
                }
            }
        }
        ValueType GetAmountInside() const
        {
            return this->amountInside;
        }
    private:
        std::array<ValueType, Dimension> getRandomPoint()
        {
            std::array<ValueType, Dimension> res;
            for(std::size_t curDim = 0; curDim < Dimension; ++curDim)
            {
                res[curDim] = this->pointGenerator[curDim](this->randomEngine);
            }
            return res;
        }
};

template<typename UnsignedIntType>
UnsignedIntType gcd(const UnsignedIntType _a, UnsignedIntType _b)
{
    UnsignedIntType a, b;
    if(_a > _b)
    {
        a = _a;
        b = _b;
    }
    else
    {
        a = _b;
        b = _a;
    }
    while(b > 0)
    {
        const UnsignedIntType div = a%b;
        a = b;
        b = div;
    }
    return a;
}

template <std::size_t Dimension, typename ValueType>
ValueType getTotalMeasure(const std::array<std::pair<ValueType, ValueType>, Dimension>& limitBox)
{
    auto prodFun = [](const ValueType tmpProd, const std::pair<ValueType, ValueType>& curDim)->ValueType
    {
        return tmpProd*(curDim.second-curDim.first);
    };
    return std::accumulate(limitBox.begin(), limitBox.end(), ValueType(1), prodFun);
}

template <std::size_t Dimension, typename ValueType, typename RandomEngineType, typename IsInsideType>
ValueType MeasureMcPar(const std::array<std::pair<ValueType, ValueType>, Dimension>& limitBox,
                       const IsInsideType& isInside,
                       const std::uintmax_t totalIter,
                       const unsigned numberOfThreads)
{
    if(totalIter%std::uintmax_t(numberOfThreads) != 0)
    {
        throw std::invalid_argument("totalIter must be a multiple of numberOfThreads.");
    }
    const std::uintmax_t localIter = totalIter/std::uintmax_t(numberOfThreads);
    std::vector<std::thread> threadVector;
    typedef McWorker<Dimension, ValueType, RandomEngineType, IsInsideType> WorkerType;
    threadVector.reserve(numberOfThreads);
    std::vector<WorkerType> workerVector;
    workerVector.reserve(numberOfThreads);
    const auto mainSeed = time(NULL);
    for(unsigned i = 0; i < numberOfThreads; ++i)
    {
        workerVector.push_back(WorkerType(limitBox, isInside, localIter, i+mainSeed));
        threadVector.push_back(std::thread(std::ref(workerVector[i])));
    }
    for(auto & th : threadVector)
    {
        th.join();
    }

    auto sumFun = [](const std::uintmax_t tmpSum, const WorkerType& inWorker)->std::uintmax_t
    {
        return tmpSum+inWorker.GetAmountInside();
    };
    const std::uintmax_t totalSum = std::accumulate(workerVector.begin(), workerVector.end(), 0ULL, sumFun);
    const auto div = gcd(totalSum, totalIter);
    const auto limitingMeasure = getTotalMeasure<Dimension, ValueType>(limitBox);
    return limitingMeasure*ValueType(totalSum/div)/ValueType(totalIter/div);
}

template <std::size_t Dimension, typename ValueType>
std::array<ValueType, Dimension> operator-(const std::array<ValueType, Dimension>& startPos,
                                           const std::array<ValueType, Dimension>& endPos)
{
    std::array<ValueType, Dimension> res;
    for(std::size_t i = 0; i < Dimension; ++i)
    {
        res[i] = endPos[i]-startPos[i];
    }
    return res;
}

template <std::size_t Dimension, typename ValueType>
ValueType normL2(const std::array<ValueType, Dimension>& inVec)
{
    auto sumSq = [](const ValueType& tmpSum, const ValueType& inVal)->ValueType
    {
        return tmpSum+inVal*inVal;
    };
    return std::sqrt(std::accumulate(inVec.begin(), inVec.end(), ValueType(0), sumSq));
}


template<typename ValueType>
class IsInsideBumpy
{
    public:
        bool operator()(const std::array<ValueType, 3>& inPos) const
        {
            bool res = true;
            const auto r = normL2(inPos);
            if(r > 0.8)
            {
                const auto x = inPos[0];
                const auto y = inPos[1];
                const auto z = inPos[2];
                const auto t = std::atan2(x, y);
                const auto v = std::acos(z/r);
                res = r <= 1.0L+0.2L*std::sin(6.0L*t)*std::sin(5.0L*v);
            }
            return res;
        }
};

int main()
{
    typedef long double ValueType;
    typedef IsInsideBumpy<ValueType> IsInsideType;
    const IsInsideType isInObj = IsInsideType();
    const ValueType limitSize = 1.2L;
    const std::array<std::pair<ValueType, ValueType>, 3> limitBox = {std::make_pair(-limitSize, limitSize),
                                                                     std::make_pair(-limitSize, limitSize),
                                                                     std::make_pair(-limitSize, limitSize)};
    const unsigned numberOfThreads = std::thread::hardware_concurrency();
    const std::uintmax_t totalIter = std::uintmax_t(numberOfThreads)*50000000ULL;

    const auto startTime = std::chrono::steady_clock::now();
    const auto res = MeasureMcPar<3, ValueType, std::mt19937_64, IsInsideType>(limitBox,
                                                                               isInObj,
                                                                               totalIter,
                                                                               numberOfThreads);
    const auto endTime = std::chrono::steady_clock::now();
    std::chrono::duration<double> runTime = endTime-startTime;
    std::cout<<"Runtime: "<<runTime.count()<<" [s]"<<std::endl;
    std::cout<<"Number of threads: "<<numberOfThreads<<std::endl;
    std::cout<<"Number of trials: "<<totalIter<<std::endl;
    std::cout<<"result: "<<res<<std::endl;
    return 0;
}

0 Helpful 0 Interesting 0 Brilliant 0 Confused
Frank Petiprin
Nov 1, 2019 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
 
python
#Bumpy Sphere  October 2019
'''
The program generates a random PT(x,y,z) and converts PT into a spherical
coordinate (Random_ro,t,v). Evaluates the function ro=1+.2*sin(6*t)*sin(5*v)
and if Randon_ro < ro then PT is contained in the Bumpy Sphere.
'''
import numpy as np
#main
seed=9089
#np.random.seed(seed) #Use as needed for Testing.
Trials,Att,TarS,TotalVol=10,100000000,1.25,0 #Att is size of column arrays (Rx,Ry,Rz)
#Used a GeoGebra graph to approximate the TargetBox size.
TargetBox=(2*TarS)**3
for j in range(1,Trials+1):
    Rx,Ry=np.random.uniform(-TarS,TarS,Att),np.random.uniform(-TarS,TarS,Att)
    Rz=np.random.uniform(-TarS,TarS,Att)
    Rro=np.sqrt(Rx*Rx+Ry*Ry+Rz*Rz)
    t=np.arctan(Ry/Rx)
    v=np.arccos(Rz/Rro)
    ro=1+.2*(np.sin(6*t)*np.sin(5*v)) # (ro,t,v) is on the Bumpy Sphere
    TorF=(Rro<=ro)
    TotHits=sum(TorF) #sums number of True
    print('total hits',TotHits,TotHits/Att,'Trial',j,"Attempts",Att)
    print('Volume',(TotHits/Att)*TargetBox)
    TotalVol+=(TotHits/Att)*TargetBox
print('Average Volume after ',Trials,' Trials is',TotalVol/Trials)

'''
total hits 27619258 0.27619258 Trial 1 Attempts 100000000
Volume 4.3155090625
total hits 27618813 0.27618813 Trial 2 Attempts 100000000
Volume 4.31543953125
total hits 27619790 0.2761979 Trial 3 Attempts 100000000
Volume 4.3155921875
total hits 27625644 0.27625644 Trial 4 Attempts 100000000
Volume 4.316506875
total hits 27615778 0.27615778 Trial 5 Attempts 100000000
Volume 4.3149653125
total hits 27625814 0.27625814 Trial 6 Attempts 100000000
Volume 4.3165334374999995
total hits 27616030 0.2761603 Trial 7 Attempts 100000000
Volume 4.3150046875
total hits 27624056 0.27624056 Trial 8 Attempts 100000000
Volume 4.31625875
total hits 27618524 0.27618524 Trial 9 Attempts 100000000
Volume 4.315394375
total hits 27620073 0.27620073 Trial 10 Attempts 100000000
Volume 4.31563640625
Average Volume after  10  Trials is 4.315684062500001

0 Helpful 0 Interesting 0 Brilliant 0 Confused

0 pending reports

×

Problem Loading...

Note Loading...

Set Loading...