# Write your code here

import numpy as np
import matplotlib.pyplot as plt

#Parameters of the distribution
mu = 0.1
sigma = 0.2
S0 = 50

#We want to simulate npaths=10000 paths of the r.v
#for tha we discretize the time from t=0 to T=10 in n=100 steps

n = 100
npaths=10000
T=10
dt =T/n

np.random.seed(1)

s = np.exp(
    (mu - sigma ** 2 / 2) * dt
    + sigma * np.random.normal(0, np.sqrt(dt), size=(npaths,n)).T
)
    
s = np.vstack([np.ones(npaths), s])
s = S0 * s.cumprod(axis=0)

#This is for scaling the x axis correctly to 0 to T=10 in 100 steps
time=np.linspace(0,T,n+1)
t=np.full(shape=(npaths,n+1),fill_value=time).T

#Plot the paths
plt.plot(t,s)
plt.xlabel("$t$")
plt.ylabel("$x$")
plt.title(
    "10000 paths of GBM for $\mu=0.1$ and $\sigma=0.2$"
)
plt.show()

#We calculate the expectation by averaging over the samples
expectation_numerical=sum(s[100])/npaths
expectation_numerical

# E_ST
#We compare with the theoretical result given by the analytical solution
#E=S_0*exp(mu*t)

expectation_analytical=S0*np.exp(mu*10)
expectation_analytical

#As we can see we get quite similar results 




%matplotlib inline


from qiskit import QuantumCircuit,  IBMQ
from qiskit import Aer
from qiskit_ibm_runtime import QiskitRuntimeService, Session, Estimator, Options
from qiskit.algorithms import IterativeAmplitudeEstimation, EstimationProblem, MaximumLikelihoodAmplitudeEstimation
from qiskit.circuit.library import LinearAmplitudeFunction
from qiskit_aer.primitives import Sampler
from qiskit_finance.circuit.library import LogNormalDistribution


# Loading your IBM Q account(s)
IBMQ.save_account('a68087d9c262053b5e20eea6f7d5dd9a02252efbc722cc8774bde7a6b69f71b9b077e3cd8541be84c5841c130f3d8f1c4d3940bf42ccbd3f1caae0d843c42976')

# Write your code here

# number of qubits to represent the uncertainty
num_samples_qubits = 4

# parameters for considered random distribution
S0 = 50  
T = 10

#parameters for log-normal distribution
mu = 0.1
sigma = 0.2

mean = np.exp(mu + sigma**2 / 2)
variance = (np.exp(sigma*2) - 1) * np.exp(2 * mu + sigma*2)
stddev = np.sqrt(variance)


#lowest and highest value considered 
low = np.maximum(0, mean - 3 * stddev)
high = mean + 3* stddev

#We construct a circuit to load a log-normal random distribution into a quantum state.
preparation_model = LogNormalDistribution(
    num_samples_qubits, mu=mu, sigma=sigma**2, bounds=(low, high)
)

preparation_model.num_qubits

# plot probability distribution
x = preparation_model.values
y = preparation_model.probabilities
plt.bar(x, y, width=0.2)
plt.xticks(x, size=15, rotation=90)
plt.yticks(size=15)
plt.grid()
plt.xlabel("$S_T$", size=15)
plt.ylabel("Probability", size=15)
plt.show()

# Write your code here

sim=Aer.get_backend('qasm_simulator')



num_samples_qubits = 4 #this number should be a power of 2
num_qubits = num_samples_qubits + 1 #add the target qubit that we will measure

circuit = QuantumCircuit(num_qubits)
circuit.append(preparation_model, range(num_samples_qubits))

def controlled_rotations(num_control_qubits, angle):
    qc = QuantumCircuit(num_control_qubits+1)
    #do the  controlled rotations
    for i in range(num_control_qubits-1, -1, -1):
        qc.cry(2**(num_control_qubits-1-i)*angle, control_qubit=i,target_qubit=num_control_qubits)
    return qc.to_gate()

#The bounds of the support of f(x)
x_l = 0.01
x_u = 0.3*sigma

num_fourier_terms = 4
P = (x_u-x_l) #period of the periodic function
delta = (x_u-x_l)/(2**num_samples_qubits-1)  #step in the discretized space

#Initialize sampler to run circuits
#service = QiskitRuntimeService(channel='ibm_quantum')
#backend = service.backend('ibmq_qasm_simulator') #choose the simulator or device
sampler = Sampler()


#Arrays where fourier terms are stored
fourier_terms_cos = []
fourier_terms_sin = []

for n in range(num_fourier_terms):
    w = 2*np.pi/P  #phase of the fourier terms
    theta = n*w*delta  

    alpha_cos=2*(np.pi)*x_l
    alpha_sin=2*(np.pi)*x_l+np.pi/2

    #circuit must be reinitialized everytime
    circuit_fourier=circuit.copy()

    #distinguish two cases:

    #If we want to estimate the cos term:
    circuit_sin=circuit_fourier.copy()
    circuit_sin.ry(alpha_cos,num_samples_qubits)
    circuit_sin.append(controlled_rotations(num_samples_qubits,theta),range(num_qubits))
    #run QAE for the cos term
    problem = EstimationProblem(circuit_sin, objective_qubits = num_samples_qubits)
    mlae = MaximumLikelihoodAmplitudeEstimation(evaluation_schedule = 3, sampler = sampler)
    mlae_result = mlae.estimate(problem)
    fourier_terms_cos.append(1-2*mlae_result.estimation)
    print("Estimate of a_",n+1,"is",1-2*mlae_result.estimation)

    
    #If we want to estimate the sin term:
    circuit_cos=circuit_fourier.copy()
    circuit_cos.ry(alpha_sin,num_samples_qubits)
    circuit_cos.append(controlled_rotations(num_samples_qubits,theta),range(num_qubits))
    #run QAE for the sin term
    problem = EstimationProblem(circuit_cos, num_samples_qubits)
    mlae = MaximumLikelihoodAmplitudeEstimation(evaluation_schedule = 3, sampler = sampler)
    mlae_result = mlae.estimate(problem)
    fourier_terms_sin.append(1-2*mlae_result.estimation)
    print("Estimate of b_",n+1,"is",1-2*mlae_result.estimation,"\n")

circuit_sin.draw(output='mpl', style={'backgroundcolor': '#EEEEEE'})

#New version to compute coefficients
import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt

#x = np.linspace(45, 144, 100)
#y = np.piecewise(x, [((x >= 45) & (x <= 100)), ((x >= 100) & (x < 145))], [0, lambda x: x-100])
#plt.plot(x, y)
def f(x):
    return np.piecewise(x, [((x >= 45) & (x <= 100)), ((x >= 100) & (x < 145))], [0, lambda x: x-100])
    

#t is the independent variable
x_l = 50
x_u = 107
P = (x_u-x_l) #period value

BT=0. #initian value of t (begin time)
ET=10. #final value of t (end time)
FS=100 #number of discrete values of t between BT and ET

#all discrete values of t in the interval from BT and ET
t_range = np.linspace(BT, ET, FS)
y_true = f(t_range) #the true f(t)

#function that computes the real fourier couples of coefficients (a0, 0), (a1, b1)...(aN, bN)
def compute_real_fourier_coeffs(func, N):
    result = []
    for n in range(N+1):
        an = (2./P) * spi.quad(lambda t: func(t) * np.cos(2 * np.pi * n * t / P), 0, P)[0]
        bn = (2./P) * spi.quad(lambda t: func(t) * np.sin(2 * np.pi * n * t / P), 0, P)[0]
        result.append((an, bn))
    return np.array(result)

def fit_func_by_fourier_series_with_real_coeffs(t, AB):
    result = 0.
    A = AB[:,0]
    B = AB[:,1]
    for n in range(0, len(AB)):
        if n > 0:
            result +=  A[n] * fourier_terms_cos[n] + B[n]* fourier_terms_sin[n]
        else:
            result +=  A[0]
    return result



def f(t):
    return max(np.exp(mu*t) - 100, 0)

AB = compute_real_fourier_coeffs(f,num_fourier_terms-1 )
est = fit_func_by_fourier_series_with_real_coeffs(0, AB)
print(est)





#The results for different values of n (precission of Fourier)

#44.698394965554456 for n=4
#41.10875168637083 for n=8 
#48.88444858978799 for n=12
#37.715743030433934 for n=16 
#45.077490680927184 for n=20
#48.90604298278019 for n=32



#We choose the Maximum Likelihood Amplitude estimation because it is one of the fastest ones and it is
#also subject to paralelization for a more efficient use



# Write your code here

# Write your code here

# Build your presentation and rehearse your pitch!