############################################################################
## CALCULUS - BIOINFORMATICS, 2018-19 T3
## Authors: Toni Guillamon, Alex Haro
############################################################################

##################################################################################################################
### Integrating tables of points, not necessary to know the underlying function f(x)                           ###
### https://docs.scipy.org/doc/scipy/reference/integrate.html#module-scipy.integrate                           ###
##################################################################################################################

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import*

### SOME PRELIMINARIES

### We work with double precision

###from decimal import *
###getcontext().prec = 16
###Decimal(1) / Decimal(7)

### Different ways to define functions
def fsin(x):
    y=np.sin(x**2)
    return y
  
### A) TRAPEZOIDS' AND SIMPSON'S RULES

## A.1) Trapezoids' rule from Scipy integration package
x = np.linspace(0, 2, num=17, endpoint=True)# Uniformly spaced, but including endpoint; it is n=16
y = fsin(x);
np.trapz (y,x)

## A.2) A self-made implementation of Trapezoids' rule.

#Warning: Remind to check that you have defined the function f3 that you want.
n=16 #number of intervals
a=0 #endpoints
b=2
h=(b-a)/n #integration step
s=fsin(a)+fsin(b) #we evaluate at the endpoints
for j in range(1,n): #loop to compute the inner components of the trapezoids' rule
    xj=a+h*j
    s=s+2*fsin(xj)

sT=s*h/2
sT

## A.3) Simpson's rule from Scipy integration package

x = np.linspace(0, 2, num=16+1, endpoint=True) # Uniformly spaced, but including endpoint; it is n=16
y = fsin(x)
simps(y,x) ### compare the agreement with Slide 8 from Lecture 9.

## A.4) A self-made implementation of Simpson's rule.

#Warning: Remind to check that you have defined the function f3 that you want.
n=16384 #number of intervals (must be EVEN)
a=0 #endpoints
b=2
h=(b-a)/n #integration step
s=fsin(a) #we evaluate at the endpoints
for j in range(1,int(n/2)+1): #loop to compute the inner components of the trapezoids' rule
    x1=a+h*(2*j-1)
    x2=a+h*2*j
    s=s+4*fsin(x1)+2*fsin(x2)

s=s-fsin(b) #we subtract an extra f(b) computed inside the loop
sS=s*h/3
sS


### B) ROMBERG'S METHOD (extrapolation)

#Default values: r=romberg(fsin, 0, 2, tol=1.48e-08, rtol=1.48e-08, show=False, divmax=10, vec_func=False)
r=romberg(fsin, 0, 2)
r

### For more details, try the function "RombergScheme" from file "Lab080-2019-RombergMethod.py" to retrieve the table "Extrapolation" given in "trapezoids.ods".