#!/usr/bin/env python3
# -*- coding: utf-8 -*-

#%%
import numpy as np
import matplotlib.pyplot as plt

def f(x):
    ''' Paraboloid '''
    return x[0]**2 + x[1]**2

def H(x): 
    ''' Indikatorfunktion für die Menge {x : f(x) < 1} '''
    return f(x) < 1

# 1 Hoehenlinie von f
x0 = np.linspace(-1, 1)
x1 = np.linspace(-1, 1)
X0, X1 = np.meshgrid(x0, x1) # Gitterpunkte im Quadrate [-1,1] x [-1,1] 
Z = f([X0, X1])
plt.figure(1)
plt.contour(x0, x1, Z, [0.1, 0.5, 0.9, 1])  #  Die Liste [1] sorgt dafür, dass Höhenlinien zu den Werten in der Liste gezeichnet werden
plt.axis('equal')


#%%
def MCquad(Hfh, Ns, xmin, xmax):
    ''' 
    Monte Carlo Quadratur für Indikatorfunktion H 
    Hfh        function handle zur Auswertung von H
    N          integer Anzahl der zufällig gewählten Werte
    xmin, xmax Listen die die linke untere und rechte obere Ecke des Quaders beschreiben 
        in den H nicht Null ist.
    '''
    
    rx = [ xmin[i] + (xmax[i]-xmin[i])*np.random.rand(Ns) for i in range(2) ]
    rH = Hfh(rx)
    ind = rH == 1

    plt.figure(Ns) # das geht so nur in 2D
    plt.plot(rx[0][ind], rx[1][ind], 'r.')
    plt.plot(rx[0][~ind], rx[1][~ind], 'b.')
    plt.axis('equal')
    plt.show()

    return sum(ind)/Ns*np.prod([ xmax[i]-xmin[i] for i in range(2) ])


for N in [100, 1000, 10000, 1000000]:
    print(MCquad(H, N, [-1, -1], [1, 1]) - np.pi)