Monte Carlo Integration OR throwing Random Darts(Ordered Pairs (x,y)) at a target.

Computer Science Level 4

Given two sets of functions

SetA:

r e d F ( x ) = x 5 + 3.5 ∗ x 4 − 2.5 x 3 − 12.5 x 2 + 1.5 x + 9 redF(x)=x^5 + 3.5*x^4 - 2.5x^3 - 12.5x^2 + 1.5x+9

b l u e F ( x ) = x 5 + x 4 − 3 x 3 − x 2 + 2 x blueF(x)=x^5 + x^4 - 3x^3 - x^2 + 2x

SetB:

g r e e n F ( x ) = 0.7 x 5 + x 4 − 3.5 x 3 − 5.2 x 2 + 6.5 x + 25 greenF(x)=0.7x^5+x^4-3.5x^3-5.2x^2+6.5x+25

p u r p l e F ( x ) = − 0.6 x 5 − x 4 + 2.5 x 3 + 3.8 x 2 − 4.1 x + 18.3 purpleF(x)= -0.6x^5-x^4+2.5x^3+3.8x^2-4.1x+18.3

Find the area between the curves of SetA on the interval [ -2 .1, 2 ] and the area between the curves of SetB on the interval [ -2.1 , 1.8 ] (Give the answer as the sum of the two areas)


The answer is 41.8.

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 21, 2020

My parallel C++ code: 

  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
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
 
#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(), 1.0L, 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 <typename ValueType>
bool isBetween(const ValueType& inVal, const ValueType& endA, const ValueType& endB)
{
    return (endA <= inVal && inVal <= endB) || (endB <= inVal && inVal <= endA);
}

template <typename FunTypeA, typename FunTypeB, typename ValueType>
class IsBetweenFunctions2D
{
    const FunTypeA funA;
    const FunTypeB funB;
    public:
        IsBetweenFunctions2D(const FunTypeA _funA, const FunTypeB _funB) :
            funA(_funA),
            funB(_funB)
        {}
        bool operator()(std::array<ValueType, 2> inPos) const
        {
            const auto x = inPos[0];
            const auto y = inPos[1];
            const auto valA = this->funA(x);
            const auto valB = this->funB(x);
            return isBetween<ValueType>(y, valA, valB);
        }
};

class SimpleFun
{
    private:
        long double (*const fun)(const long double);
    public:
        SimpleFun(long double (*const _fun)(const long double)) :
            fun(_fun)
        {}
        long double operator()(long double inVal) const
        {
            return this->fun(inVal);
        }
};

long double redF(const long double x)
{
    return ((((x+3.5)*x-2.5)*x-12.5)*x+1.5)*x+9;
}

long double blueF(const long double x)
{
    return (((((x+1)*x)-3)*x-1)*x+2)*x;
}

long double greenF(const long double x)
{
    return ((((0.7*x+1)*x-3.5)*x-5.2)*x+6.5)*x+25;
}

long double purpleF(const long double x)
{
    return ((((-0.6*x-1)*x+2.5)*x+3.8)*x-4.1)*x+18.3;
}

int main()
{
    typedef IsBetweenFunctions2D<SimpleFun, SimpleFun, long double> IsInsideType;
    const SimpleFun sRedF(&redF), sBlueF(&blueF), sGreenF(&greenF), sPurpleF(&purpleF);

    const IsInsideType isInSetA(sRedF, sBlueF);
    const std::array<std::pair<long double, long double>, 2> boxA = {std::make_pair(-2.1L, 2.0L), std::make_pair(-2.9L, 30.0L)};

    const IsInsideType isInSetB(sGreenF, sPurpleF);
    const std::array<std::pair<long double, long double>, 2> boxB = {std::make_pair(-2.1L, 1.8L), std::make_pair(11.6L, 26.7L)};


    const unsigned numberOfThreads = std::thread::hardware_concurrency();
    const std::uintmax_t totalIter = std::uintmax_t(numberOfThreads)*500000000ULL;

    const auto startTime = std::chrono::steady_clock::now();
    const auto resA = MeasureMcPar<2, long double, std::mt19937_64, IsInsideType>(boxA,
                                                                                  isInSetA,
                                                                                  totalIter,
                                                                                  numberOfThreads);
    const auto resB = MeasureMcPar<2, long double, std::mt19937_64, IsInsideType>(boxB,
                                                                                  isInSetB,
                                                                                  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<<"result: "<<resA+resB<<std::endl;
    return 0;
}

1 Helpful 0 Interesting 1 Brilliant 0 Confused
Frank Petiprin
May 6, 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
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
 
import random
#+++++++++++++++++ Using Pythonista on an Ipad Pro +++++++++++++++++++
def my_range5(start, end, step):
    while start <= end:
        yield start
        start += step

def redF(x):
      return x**5 + 3.5*x**4 - 2.5*x**3 - 12.5*x**2 + 1.5*x+9

def blueF(x):
      return x**5 + x**4 - 3*x**3 - x**2 + 2*x

def greenF(x):  
      return 0.7*x**5+1*x**4-3.5*x**3-5.2*x**2+6.5*x+25

def purpleF(x):
      return -0.6*x**5-1*x**4+2.5*x**3+3.8*x**2-4.1*x+18.3
#++++++++++++++++++++++++++
max,min,start,finish,step,eps=-9999,9999,-2.1,2,0.001,0.002
InterRBx=[0]*4 #ArrayStore x intercepts for redF(x) and blueF(x) 
index=0 #Array subscript for InterRB 

#Find the max and min value of the red function on interval [-2.1,2]
#The area of the targetArea is (max-min)*(2-(-2.1))
#The follwing code can be replaced by Python Routines to fine max,min and (x,y) intercepts

for x in my_range5(start,finish,step):
    red,blue=redF(x),blueF(x)
    purple,green=purpleF(x),greenF(x)
    if (red>max):max,maxx=red,x
    if(red<min):min,minx=red,x
#Find the x and y  intercepts of the two functions redF(x) and blueF(x).
    if(abs(red-blue)<=eps):
        InterRBx[index]=x
        index+=1
#Find the x and y  intercept of the two functions purpleF(x) and greenF(x).
    if(abs(purple-green)<=eps):InterPGx=x
targetArea=abs((max-min)*(2-(-2.1)))
#+++++++++++++++++++
seed=9007
#random.seed(seed)
Dartosses=100000000
TotArea=0
Trials=10
for i in range(1,Trials+1):
    Hits=0
    for j in range(1,Dartosses+1):
        x,y=random.uniform(-2.1,2),random.uniform(min,max)
        TFA0=((blueF(x)<=y<=redF(x)) and (-2.1<=x<=InterRBx[0]))
        TFA1=((redF(x)<=y<=blueF(x)) and (InterRBx[0]<x<=InterRBx[1]))
        TFA2=((blueF(x)<=y<=redF(x)) and (InterRBx[1]<x<=InterRBx[2]))
        TFA3=((redF(x)<=y<=blueF(x)) and (InterRBx[2]<InterRBx[3]))
        TFA4=((blueF(x)<=y<=redF(x)) and (InterRBx[3]<x<=2))
        TFA5=((greenF(x)<=y<=purpleF(x)) and (-2.1<x<=InterPGx))
        TFA6=((purpleF(x)<=y<=greenF(x)) and (InterPGx<x<=1.8))
        if(TFA0 or TFA1 or TFA2 or TFA3 or TFA4 or TFA5 or TFA6):Hits+=1
    Area=(Hits/Dartosses)*targetArea
    print('Area, Hits, Dart-Tosses',Area,Hits,Dartosses,i)
    TotArea+=Area
AvgArea=float(TotArea/Trials)
print('AVERAGE AREA, Tosses, Trials',AvgArea,Dartosses,Trials)
print('End Program')
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Area, Hits, Dart-Tosses 41.82497984989677 31088738 100000000 1
Area, Hits, Dart-Tosses 41.82021061308929 31085193 100000000 2
Area, Hits, Dart-Tosses 41.82454126845834 31088412 100000000 3
Area, Hits, Dart-Tosses 41.815718516699945 31081854 100000000 4
Area, Hits, Dart-Tosses 41.810033102102366 31077628 100000000 5
Area, Hits, Dart-Tosses 41.81585708690903 31081957 100000000 6
Area, Hits, Dart-Tosses 41.82009087766591 31085104 100000000 7
Area, Hits, Dart-Tosses 41.82759653976717 31090683 100000000 8
Area, Hits, Dart-Tosses 41.813583459206725 31080267 100000000 9
Area, Hits, Dart-Tosses 41.81961462665608 31084750 100000000 10
AVERAGE AREA , Tosses, Trials 41.81922259404516 100000000 10
End Program
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

1 Helpful 1 Interesting 1 Brilliant 0 Confused

0 pending reports

×

Problem Loading...

Note Loading...

Set Loading...