下面是题目,需要作业答案请联系本人微信scratch8
"""
## About This Homework 关于家庭作业
The homework consists of 12 questions, for a total of 86 points. 家作由12个问题组成，共86个要点。
In this module, we will write code to represent regions of space. 在本模块中，我们要写代码，表示空间中的区域。
We will represent a _region_ as the union of rectangles whose sides are parallel to the Cartesian axes.我们将使用一些和笛卡尔坐标系平等的矩形来表示区域。
Of course, this does not suffice to represent all possible regions of space,当然这并不能表示所有的空间区域。
but by using lots of small rectangles, we can at least approximate most (continuous, etc) regions. 但是靠细分矩形，我们能近似的表示大多数区域。
Since a region is the union of rectangles, let us turn our attention to rectangles. 由于区域由一些矩形组成，我们专注于矩形。
## Rectangles 矩形
A 1D rectangle is simply an interval. 一维矩形只是个间隔。
We can think of a 2D rectangle, rather than as a collection of vertices,我们认为2维矩形是间隔的集合。
as the intersection of two intervals: one for the x-axis, and one for the y-axis.由于是两个间隔，一个为x轴， 一个为y轴。
Similarly, a 3D rectangle can be thought of as the intersection of three intervals, on the x, y, and z axes.三维矩形就是在x,y,z三个轴上的间隔。
To compute the intersection of two rectangles, 计算矩形重叠，
we just need to compute the intersection of their intervals，我们只要计算间隔的重叠即可。
on the respective axes. Thus, intervals provide a simple representation that generalizes well to multiple dimensions. 在各自的轴上，表示一个间隔非常简单，可以很好地推广到多个维度。
We start building our representation from intervals. We then use lists of intervals to represent rectangles,我们将从如何表示一个间隔开始，用一系列间隔，来表示矩形。
and unions of rectangles to represent regions of space. 用多个矩形来表示一个区域。
## Intervals 间隔
An interval is defined by its two endpoints. 一个间隔用两个端点来定义。
We keep them sorted, which will make it (a lot) easier to operate on them. 这两个端点将进行排序，这样以后更方便操作。
"""
class Interval(object):
def __init__(self, x0, x1):
# 对端点进行排序，确保 x0 <= x1.
x0, x1 = (x0, x1) if x0 < x1 else (x1, x0)
assert x0 < x1 # No point intervals. self.x0 = x0 self.x1 = x1 @property def length(self): return self.x1 - self.x0 def endpoints(self): return (self.x0, self.x1) def __getitem__(self, i): """定义访问端点的方法.""" if i == 0: return self.x0 elif i == 1: return self.x1 raise KeyError() def __repr__(self): return "[{},{}]".format(self.x0, self.x1) """The main operations we need on intervals, to do anything interesting, are: “间隔”类需要的主要方法，在下面一一说明: * **Intersection.** 交叉 Given two intervals `i` and `j`, 给定两个间隔，i和j。 we want to define `__and__` for an interval so that `i & j` will be either `None`,if `i` and `j` have no intersection,如果没有交叉，i&j返回None。 or the interval corresponding to the intersection of `i` and `j`. 如果有交叉，就返回新的交叉区域的间隔。 * **Union.** 联合 The union of two intervals is not necessarily an interval; it could also be _two_ intervals, 两个间隔的联合不一定是一个间隔，或许是两个间隔。 if the original intervals are disjoint and there is a gap in between. 如果两个间隔是没有关联的，中间有间隙。 Thus, we define `__or__` so that `i | j` returns a _list_ consisting of 1 or 2 intervals. 那么，我们定义__or__方法，让i|j返回一个或者两个间隔的列表。 * **Difference.** 差集 The difference `i - j` is the portion of `i` that is not in `j`. 让i减去j，结果是保留i的，但去除i和j所相同的间隔。 The result is a _list_ of intervals, containing 0 intervals (if `j` includes `i`), 结果是个列表，如果j完全在i中，则是空列表。 one interval (if `j` does not overlap `i`, or if it overlaps only from one side of `i`),如果j不和i重叠，或者j在另一端和i重叠，则列表中就只有一个间隔。 or two intervals (if `j` falls in the middle of `i`). 如果j在i中间，则列表中有两个间隔。 * **Equality.** 相等 Two intervals are equal if, well, they are equal. 如果两个间隔(端点)都相等，则它们是相等的。 * **Membership.** 成员 We test for a point belonging to an interval. 我们测试一个点是否在间隔之间 In defining these operations, we disregard isolated points, and we blur the distinction between open and closed intervals. After all, we only care about representing regions of space, and so isolated points and things that have no extension, or no volume, are not a concern of ours. So for instance, if we subtract the interval $[3, 5]$ from $[0, 4]$, the result will be simply the interval $[4, 5]$: we do not track whether the interval is open or closed at 5. Likewise, point-wise intervals such as $[5, 5]$ are simply not considered. 在定义这些操作时，我们忽略了孤立点，模糊了开区间和闭区间的区别。毕竟，我们只关心表示空间的区域，因此没有扩展或没有体积的孤立点和事物不是我们关心的问题。例如，如果我们从$[0，4]$中减去区间$[3，5]$，结果就是区间$[4，5]$：我们不跟踪区间是在5处打开还是关闭。同样地，像$[5，5]$这样的逐点区间也不被考虑。 ## Question 1: Interval Equality 问题1，内部相等 Let us start by implementing equality. We leave this to you.让我们从实现相等开始，我们把这个问题留给你做。 """ ### Defining Equality def interval_equality(self, other): """如果内部数据相等，则返回真，否则返回假""" ### YOUR CODE HERE (下面是我写的代码，下同) Interval.__eq__ = interval_equality i = Interval(3, 5) j = Interval(4, 5) i == j # Tests for equality. 5 points. i = Interval(3, 5) j = Interval(4, 5) assert i != j assert Interval(5, 7) == Interval(5, 7) assert Interval(2.3, 3.4) == Interval(2.3, 3.4) """### Union We define union for you, to give you an example. 我们给你定义好了，这是给你一个例子 """ def interval_or(self, other): """Union of self and other. Returns a list of 1 or 2 non-overlapping intervals.""" Interval.__or__ = interval_or # The union of these two intervals is a list of two intervals. Interval(3, 5) | Interval(7, 10) # The union of these two intervals is a single interval. Interval(3, 5) | Interval(4, 10) """## Question 2: Interval Intersection The intersection of two intervals `i` and `j` consists either of a single interval, or `None`, if the two intervals have no intersection. We leave it to you to implement it. """ ### Interval intersection def interval_and(self, other): """Intersection; returns an interval, or None.""" ### YOUR CODE HERE Interval.__and__ = interval_and # These two intervals should have empty intersection assert Interval(3, 4) & Interval(5, 6) is None # These two intervals should have non-empty intersection. assert Interval(3, 10) & Interval(6, 20) == Interval(6, 10) # 5 points: tests for intersection. assert Interval(3, 4) & Interval(5, 6) is None assert Interval(3, 10) & Interval(6, 20) == Interval(6, 10) assert Interval(-3, 10) & Interval(4, 5) == Interval(4, 5) """## Question 3: Interval Membership Given an interval `i`, and a floating point number `x`, we can write a method `__contains__` of an interval, which checks if `x` belongs to the interval. In this way, writing `x in i` will return `True` if `x` belongs to `i`, and `False` otherwise. For the purpose of this method, you can consider an interval closed, so that 3 in Interval(3, 5) returns `True`. """ ### Membership of a point in an interval def interval_contains(self, x): ### YOUR CODE HERE Interval.__contains__ = interval_contains assert 3 in Interval(3, 5) assert not (1 in Interval(3, 5)) # 5 points: tests for interval membership. assert 3 in Interval(3, 5) assert 2 not in Interval(3, 5) assert 8 not in Interval(3, 5) """## Question 4: Interval Difference 内部差 For intervals `i`, `j`, the difference of `i - j` consists of 0, 1, or 2 non-overlapping intervals. Again, we leave the implementation to you. """ ### Interval difference def interval_sub(self, other): """Subtracts from this interval the interval other, returning a possibly empty list of intervals.""" ### YOUR CODE HERE Interval.__sub__ = interval_sub assert Interval(4, 6) - Interval(5, 8) == [Interval(4, 5)] assert Interval(0, 10) - Interval(4, 5) == [Interval(0, 4), Interval(5, 10)] # 5 points: tests for interval difference. assert Interval(4, 6) - Interval(5, 8) == [Interval(4, 5)] assert Interval(0, 10) - Interval(4, 5) == [Interval(0, 4), Interval(5, 10)] assert Interval(0, 2) - Interval(-3, 6) == [] assert Interval(0, 10) - Interval(0, 5) == [Interval(5, 10)] assert Interval(-4, -2) - Interval(-3, -2) == [Interval(-4, -3)] assert Interval(4, 5) - Interval(4, 5) == [] """Another way of testing this code is the following. Let's generate many random intervals $I_1$ and $I_2$. Denoting with $-$ the difference of intervals and with $\cap$ their intersection, and denoting the length of an interval $I$ by $l(I)$, the following invariant must be true: $$ l(I_1 - I_2) + l(I_2 - I_1) + 2l(I_1 \cap I_2) = l(I_1) + l(I_2) $$ To verify this, let us start by defining this total length function precisely. """ import numpy as np def total_length(x): if x is None: return 0. elif type(x) == list: return np.sum([i.length for i in x]) else: return x.length print(total_length(None)) i1 = Interval(0, 1) i2 = Interval(3, 5) print("i1:", total_length(i1)) print("i2:", total_length(i2)) print("i1+i2:", total_length([i1, i2])) # 5 points: more tests for interval difference. import random def test_random(): i1 = Interval(random.random(), random.random()) i2 = Interval(random.random(), random.random()) d1 = i1 - i2 d2 = i2 - i1 inters = i1 & i2 assert (total_length(d1) + total_length(d2) + 2. * total_length(inters) == i1.length + i2.length) for _ in range(100): test_random() """### Rectangles 现在，让我们来开发个小项目来表示矩形，用区间的交集表示。 Let us now develop a representation of a rectangle, in terms of intersection of intervals. 我们定义的程序能在任何维度上都能正常工作，而不是(打个比方说),只在3维空间有效。 We will phrase the definition in such a way that it works in any number of dimensions, storing the intervals as a list, as opposed to (say) storing the three intervals separately for 3D. """ import string class Rectangle(object): def __init__(self, *intervals, name=None): """A rectangle is initialized with a list, whose elements are either Interval, or a pair of numbers. It would be perhaps cleaner to accept only list of intervals, but specifying rectangles via a list of pairs, with each pair defining an interval, makes for a concise shorthand that will be useful in tests. Every rectangle has a name, used to depict it. If no name is provided, we invent a random one.""" self.intervals = [] for i in intervals: self.intervals.append(i if type(i) == Interval else Interval(*i)) # I want each rectangle to have a name. if name is None: self.name = ''.join( random.choices(string.ascii_letters + string.digits, k=8)) else: self.name = name def __repr__(self): """Function used to print a rectangle.""" s = "Rectangle " + self.name + ": " s += repr([(i.x0, i.x1) for i in self.intervals]) return s def clone(self, name=None): """Returns a clone of itself, with a given name.""" name = name or self.name + "'" return Rectangle(*self.intervals, name=name) def __len__(self): """Returns the number of dimensions of the rectangle (not the length of the edges). This is used with __getitem__ below, to get the interval along a dimension.""" return len(self.intervals) def __getitem__(self, n): """Returns the interval along the n-th dimension""" return self.intervals[n] def __setitem__(self, n, i): """Sets the interval along the n-th dimension to be i""" self.intervals[n] = i @property def ndims(self): """Returns the number of dimensions of the interval.""" return len(self.intervals) @property def volume(self): return np.prod([i.length for i in self.intervals]) print(Rectangle(Interval(3., 4.), Interval(1., 4.))) r = Rectangle(Interval(1., 2.), (5., 6.), name="my_rectangle") print(r) print(r.clone()) """### Drawing rectangles Before we go much further, it is useful to be able to draw rectangles. Rectangles can have any number of dimensions, and we will write here code to draw them on 2D, projecting away all other dimensions. """ import matplotlib import matplotlib.pyplot as plt import matplotlib.path as mpath import matplotlib.patches as mpatches from matplotlib.collections import PatchCollection matplotlib.rcParams['figure.figsize'] = (6.0, 4.0) def draw_rectangles(*rectangles, prefix=""): """Here, rectangles is a rectangle iterator; it could be a list, for instance.""" fig, ax = plt.subplots() patches = [] # We keep track of the limits. lo_x, hi_x = [], [] lo_y, hi_y = [], [] for r in rectangles: x0, x1 = r[0].endpoints() y0, y1 = r[1].endpoints() lo_x.append(x0) hi_x.append(x1) lo_y.append(y0) hi_y.append(y1) # Prepares the "patch" for the rectangle, see # https://matplotlib.org/api/_as_gen/matplotlib.patches.Rectangle.html p = mpatches.Rectangle((x0, y0), x1 - x0, y1 - y0) y = (y0 + y1) / 2. - 0.0 x = (x0 + x1) / 2. - 0.0 plt.text(x, y, prefix + r.name, ha="center", family='sans-serif', size=12) patches.append(p) # Draws the patches. colors = np.linspace(0, 1, len(patches) + 1) collection = PatchCollection(patches, cmap=plt.cm.hsv, alpha=0.3) collection.set_array(np.array(colors)) ax.add_collection(collection) # Computes nice ax limits. Note that I need to take care of the case # in which the rectangle lists are empty. lox, hix = (min(lo_x), max(hi_x)) if len(lo_x) > 0 else (0., 1.)
loy, hiy = (min(lo_y), max(hi_y)) if len(lo_y) > 0 else (0., 1.)
sx, sy = hix - lox, hiy - loy
lox -= 0.2 * sx
hix += 0.2 * sx
loy -= 0.2 * sy
hiy += 0.2 * sy
ax.set_xlim(lox, hix)
ax.set_ylim(loy, hiy)
plt.gca().set_aspect('equal', adjustable='box')
plt.grid()
plt.show()
r1 = Rectangle((3., 5.), (1., 4.), name="A")
r2 = Rectangle((1., 4.), (2., 6.), name="B")
r3 = Rectangle((2., 3.5), (1.5, 5.), name="C")
draw_rectangles(r1, r2, r3)
"""There are three main operations on rectangles: intersection, union, and difference.
Among them, only intersection is guaranteed to return another rectangle.
In general, the union of two rectangles is ... two rectangles,
and the difference between two rectangles is ... a whole lot of rectangles,
as we will see.
We let you implement rectangle equality, intersection, and membership of a point in a rectangle.
**Equality:** Two rectangles $R$ and $T$ are equal if they have the same number of dimensions,
and if for every dimension $k$, the interval of $R$ along $k$ is equal to the interval of $T$ along $k$. For example,
Rectangle((2, 3), (4, 5)) == Rectangle((2, 3), (4, 5))
Rectangle((2, 3), (4, 5)) != Rectangle((4, 5), (2, 3))
Rectangle((2, 3), (4, 5)) != Rectangle((2, 3), (4, 5), (6, 7))
**Intersection:** The intersection is defined only if the rectangles have the same number of dimensions.
The intersection is computed by taking the intersection of the intervals of the two rectangles for corresponding dimensions.
**Membership:** For an $n$-dimensional point $(x_0, x_1, \ldots, x_n)$ and an $n$-dimensional rectangle $R$, we have $(x_0, x_1, \ldots, x_n) \in R$ if the point is in the region $R$. For instance:
(2.5, 4.5) in Rectangle((2, 3), (4, 5))
(2, 3) in Rectangle((2, 3), (4, 5))
(1, 3) not in Rectangle((2, 3), (4, 5))
If the point and the rectangle have different dimensions, you can raise a `TypeError`.
## Question 5: Rectangle Equality
"""
def rectangle_eq(self, other):
### YOUR CODE HERE
Rectangle.__eq__ = rectangle_eq
# 5 points: tests for rectangle equality.
assert Rectangle((2, 3), (4, 5)) == Rectangle((2, 3), (4, 5))
assert Rectangle((2, 3), (4, 5)) != Rectangle((4, 5), (2, 3))
assert Rectangle((2, 3), (4, 5)) != Rectangle((2, 3), (4, 5), (6, 7))
"""## Question 6: Rectangle Intersection """
### Rectangle intersection
def rectangle_and(self, other):
if self.ndims != other.ndims:
raise TypeError("The rectangles have different dimensions: {} and {}".format(
self.ndims, other.ndims
))
# Challenge: can you write this as a one-liner shorter than this comment is?
# There are no bonus points, note. Just for the fun.
### YOUR CODE HERE
Rectangle.__and__ = rectangle_and
# Let's see how your rectangle intersection works.
r1 = Rectangle((2, 3), (0, 4))
r2 = Rectangle((0, 4), (1, 3))
draw_rectangles(r1, r2)
draw_rectangles(r1 & r2)
# 10 points: tests for rectangle intersection.
r1 = Rectangle((2, 3), (0, 4))
r2 = Rectangle((0, 4), (1, 3))
assert r1 & r2 == Rectangle((2, 3), (1, 3))
r1 = Rectangle((2, 3), (0, 4))
r2 = Rectangle((0, 4), (1, 5))
assert r1 & r2 == Rectangle((2, 3), (1, 4))
r1 = Rectangle((-1, 5), (0, 6))
r2 = Rectangle((0, 4), (-1, 3))
assert r1 & r2 == Rectangle((0, 4), (0, 3))
r1 = Rectangle((2, 6), (0, 4))
r2 = Rectangle((0, 6), (0, 3))
assert r1 & r2 == Rectangle((2, 6), (0, 3))
"""## Question 7: Point Membership in a Rectangle"""
### Membership of a point in a rectangle.
def rectangle_contains(self, p):
# The point is a tuple with one element per dimension of the rectangle.
if len(p) != self.ndims:
raise TypeError()
### YOUR CODE HERE (下面的代码只适合二维，需要改成和维数无关，请参考上面的自行修改)
Rectangle.__contains__ = rectangle_contains
# 5 points: tests for membership.
assert (2, 3) in Rectangle((0, 4), (1, 5))
assert (0, 4) in Rectangle((0, 4), (4, 5))
assert (4, 5) in Rectangle((0, 4), (4, 5))
assert (0, 0, 0) not in Rectangle((3, 4), (0, 3), (0, 8))
"""## Regions
The problem with rectangles is that they are not closed under union: the union of two rectangles is not necessarily a rectangle.
We want a representation for objects in space that is closed under union, intersection, and difference.
To this end, we introduce _regions_, which are unions of rectangles.
"""
class Region(object):
def __init__(self, *rectangles, name=None):
"""A region is initialized via a set of rectangles."""
self.rectangles = list(rectangles)
if name is None:
self.name = ''.join(
random.choices(string.ascii_letters + string.digits, k=8))
else:
self.name = name
def draw(self):
draw_rectangles(*self.rectangles, prefix=self.name + ":")
def __or__(self, other):
"""Union of regions."""
return Region(*(self.rectangles + other.rectangles), name=self.name + "_union_" + other.name)
# Let us try.
r = Rectangle((0., 4.), (0., 4.), name="R")
t = Rectangle((1.5, 3.5), (1., 5.), name="T")
reg1 = Region(r, name="Reg1")
reg2 = Region(t, name="Reg2")
(reg1 | reg2).draw()
"""## Question 8: Membership of a Point in a Region
A point belongs into a region if it belongs into some rectangle of the region. We let you implement this.
"""
### Membership of a point in a region
def region_contains(self, p):
### YOUR CODE HERE
Region.__contains__ = region_contains
assert (2, 1) in Region(Rectangle((0, 2), (0, 3)),
Rectangle((4, 6), (5, 8)))
assert (2, 1) not in Region(Rectangle((0, 1), (0, 3)),
Rectangle((4, 6), (5, 8)))
"""## Monte-Carlo Methods
There are some obvious things we might want to do with a region, namely, compute its volume,
compute whether two regions are equal, and compute the center of mass of a region.
There are two approaches to this.
One is to develop a precise approach. The problem in computing the volume of a region is that the rectangles in it might overlap.
To solve this, one can use our method for computing disjoing differences to put regions in _normal form_,consisting of non-overlapping rectangles.
The idea is to keep a region as a list of non-overlapping rectangles.
When we add a rectangle $S$ from a region consisting of non-overlapping rectangles $R_1, \ldots, R_n$,
we first subtract from $S$ each of $R_1, \ldots, R_n$ in turn, getting as result a bunch of subrectangles of $S$;
we can then add these subrectangles to the region.
But this sounds like work!
An alternative is to develop a _randomized_ approach.
### A Monte-Carlo algorithm for region area
We can develop a randomized approach to measuring the area of a region as follows. First,
we compute a _bounding box_ around it, which is simply the smallest rectangle guaranteed to contain the region.
We simply take, for each coordinate, the min and max values of that coordinate of any rectangle in the region.
Once we have a bounding box $B$ for a region $\cal R$, we simply pick at random a lot of points $x \in B$,
using our Python `random` function. We can use our test $x \in \cal R$, written in code as `x in my_region`,
to check whether a point $x$ belongs to region $\cal R$. Let $N$ be the number of points we generate,
and $M$ be the number of points that end up in $\cal R$. The volume $V_{\cal R}$ of the region $\cal R$ can be simply written as:
$$
V_{\cal R} = \frac{M}{N} \cdot V_B \; ,
$$
where $V_B$ is the volume of the bounding box. We will lead you to implement this code in steps.
This method is an example of a [Monte Carlo method](https://en.wikipedia.org/wiki/Monte_Carlo_method),
a method which gives an answer to a question via repeated randomized experiments, rather than via mathematical computation,
which may be complex or unfeasible.
## Question 9: Compute Bounding Boxes
First, write a method `bounding_box` of a region, which returns the bounding box as a rectangle.
The bounding box is the smallest rectangle that contains the region.
"""
### Compute the bounding box of a region
def region_bounding_box(self):
"""返回区域最小包围矩形(绑定盒)
Returns the bounding box of the region, as a rectangle.
This returns None if the region does not contain any rectangle."""
if len(self.rectangles) == 0:
return None
### YOUR CODE HERE
Region.bounding_box = property(region_bounding_box)
# 10 points: tests for bounding boxes
reg = Region(Rectangle((0, 2), (1, 3)), Rectangle((4, 6), (5, 8)))
assert reg.bounding_box == Rectangle((0, 6), (1, 8))
reg = Region(
Rectangle((0, 5), (4, 5), (1, 9)),
Rectangle((4, 20), (-2, 3), (4, 21)),
Rectangle((7, 99), (3, 7), (2, 3))
)
assert reg.bounding_box == Rectangle((0, 99), (-2, 7), (1, 21))
"""### Select random points from a rectangle
Next, we write a Rectangle method `random_point`, which returns a point of a rectangle chosen uniformly at random each time it is called.
To this end, it is easier first to write the corresponding method for an interval. [In Python](https://docs.python.org/3/library/random.html#real-valued-distributions),
random.random()
returns a random value uniformly distributed between 0 and 1, and
random.uniform(a, b)
returns a random value uniformly distributed between a and b. We can use this to define the interval method:
"""
import random
def interval_random_point(self):
return random.uniform(self.x0, self.x1)
Interval.random_point = interval_random_point
# Or if we wanted to be concise, we could just have written:
Interval.random_point = lambda self : random.uniform(self.x0, self.x1)
"""## Question 10: Random Point in a Rectangle
To select a random point from a rectangle, we just need to return a tuple formed by choosing a random point from each of the rectangle's intervals.
We leave this to you. Remember that the intervals of a rectangle `self` are in `self.intervals`.
"""
### Random point of a rectangle
def rectangle_random_point(self):
### YOUR CODE HERE
Rectangle.random_point = rectangle_random_point
# 3 points: random point of a rectangle.
r = Rectangle((0, 2), (1, 3))
for i in range(5):
p = r.random_point()
assert isinstance(p, tuple)
assert len(p) == 2
assert p in r
print(p)
# 3 points: random point of a rectangle.
import numpy as np
r = Rectangle((1, 2), (1, 6))
xs, ys = [], []
for _ in range(10000):
p = r.random_point()
assert p in r
xs.append(p[0])
ys.append(p[1])
assert np.std(xs) * 4.9 < np.std(ys) < np.std(xs) * 5.1
"""## Question 11: Volume via Monte Carlo
We are now ready to compute the volume of a region $\cal R$ using Monte Carlo methods. The form of your code is:
* Compute the bounding box $B$
* Pick $N$ points at random from $B$, and count the number $M$ of them that fall in $\cal R$.
* Return $B.volume \cdot (M/N)$.
"""
### Monte carlo Volume
def region_montecarlo_volume(self, n=1000):
"""Computes the volume of a region, using Monte Carlo approximation
with n samples."""
# The solution, written without any particular trick, takes 7 lines.
# If you write a much longer solution, you are on the wrong track.
### YOUR CODE HERE
Region.montecarlo_volume = region_montecarlo_volume
"""The approximation becomes the more precise, the more samples we have. """
reg = Region(Rectangle((0, 4), (0, 2)), Rectangle((1, 2), (0, 4)))
reg.draw()
print(" 10 samples:", reg.montecarlo_volume(n=10))
print(" 100 samples:", reg.montecarlo_volume(n=100))
print(" 1000 samples:", reg.montecarlo_volume(n=1000))
print("10000 samples:", reg.montecarlo_volume(n=10000))
# 10 points: Volume via Monte Carlo
reg = Region(Rectangle((0, 1), (0, 4)), Rectangle((0, 4), (0, 1)))
reg.draw()
print(" 10 samples:", reg.montecarlo_volume(n=10))
print(" 100 samples:", reg.montecarlo_volume(n=100))
print(" 1000 samples:", reg.montecarlo_volume(n=1000))
print("10000 samples:", reg.montecarlo_volume(n=10000))
v = reg.montecarlo_volume(n=10000)
assert 6.2 < v < 7.8
"""We could quantify the standard deviation of the result, but it is beyond the scope of this class.
**Exercise:** Develop a Monte-Carlo method for computing the center of mass of a region.
The idea consists in sampling uniformly at random from the bounding box, retaining only the points that are in the region.
The center of mass of the sampled points in the region provides an approximation for the center of mass of the region.
## Question 12: A Monte-Carlo method for region equality
We can apply Monte Carlo methods also to the question of deciding region equality.
One way to decide whether two regions $\cal R_1$ and $\cal R_2$ are equal consists in subtracting $\cal R_2$ from $\cal R_1$ and
checking that the result is empty, and then subtracting $\cal R_1$ from $\cal R_2$, and checking that it is also empty.
But again, this sounds like work, and why work if we can just guess?
The idea is to compute the bounding box $B$ of $\cal R_1 \cup \cal R_2$, and to sample points from $B$.
If we find a point $p$ that belongs to one region but not the other, we declare the regions distinct.
If we do not find such "distinguishing" point after $N$ trials, we declare the regions identical,
and the point serves as a witness to their difference.
We leave the implementation to you.
"""
### Monte Carlo difference and equality between regions
def region_montecarlo_difference(self, other, n=1000):
"""Checks whether a region self is different from a region other, using
a Monte Carlo method with n samples. It returns either a point p that
witnesses the difference of the regions, or None, if no such point is found."""
# This can be done without hurry in 6 lines of code.
### YOUR CODE HERE
Region.montecarlo_difference = region_montecarlo_difference
def region_montecarlo_equality(self, other, n=1000):
return self.montecarlo_difference(other, n=n) is None
Region.montecarlo_equality = region_montecarlo_equality
reg1 = Region(Rectangle((0, 4), (0, 2)), Rectangle((1, 2), (0, 4)), name="reg1")
reg2 = Region(Rectangle((0, 4), (0, 2)), Rectangle((1, 2), (1, 4)), name="reg2")
reg3 = Region(Rectangle((0, 4), (0, 2)), Rectangle((1, 2), (0, 3)), name="reg3")
reg1.draw()
reg2.draw()
reg3.draw()
print("reg1 vs reg2", reg1.montecarlo_equality(reg2))
print("reg1 vs reg2", reg1.montecarlo_equality(reg3))
print("reg1 vs reg3", reg1.montecarlo_difference(reg3))
# 10 points: Equality of regions via Monte Carlo.
reg1 = Region(Rectangle((0, 4), (0, 2)), Rectangle((1, 2), (0, 4)), name="reg1")
reg2 = Region(Rectangle((0, 4), (0, 2)), Rectangle((1, 2), (1, 4)), name="reg2")
reg3 = Region(Rectangle((0, 4), (0, 2)), Rectangle((1, 2), (0, 3)), name="reg3")
assert reg1.montecarlo_equality(reg2)
assert not reg1.montecarlo_equality(reg3)
