from scipy.integrate import quad
import numpy as np
from prettytable import PrettyTable


def integral_solution(interval_bounds, family):
    answers = []

    for function in family:
        answers.append(quad(lambda x: eval(function), interval_bounds[0], interval_bounds[1])[0])

    return answers


def midpoint_rule(interval_bounds, partitions, family):
    intervals = list()  # array of all dx's
    dx = ((interval_bounds[1] - interval_bounds[0]) / partitions)
    step = interval_bounds[0]

    while step <= interval_bounds[1]:
        intervals.append(step)
        step += dx

    answers = [0.0] * len(family)
    i = 0
    while i < len(intervals) - 1:
        x = ((intervals[i] + intervals[i+1]) / 2.0)  # evaluate each with x
        for j in range(0, len(family)):
            answers[j] += 2*(eval(family[j]))
        i += 1
    return answers


def trapezoidal_rule(interval_bounds, partitions, family):
    x = list()  # array of all dx's
    dx = ((interval_bounds[1] - interval_bounds[0]) / partitions)
    step = interval_bounds[0]

    while step <= interval_bounds[1]:
        x.append(step)
        step += dx

    x = np.array(x)

    answers = []
    for function in family:
        if function.isdigit():
            answers.append(np.trapz(int(function) + (0*x), x))  # check if constant
        else:
            answers.append(np.trapz(eval(function), x))
    return answers


def find_weight(midpoint_rule_value, trapezoidal_rule_value, integral_value, family):
    answers = [0] * len(family)
    i = 0
    while i <= len(answers) - 1:
        if trapezoidal_rule_value[i] - midpoint_rule_value[i] == 0:
            answers[i] = 0
        else:
            answers[i] = (trapezoidal_rule_value[i] - integral_value[i]) / \
                         (trapezoidal_rule_value[i] - midpoint_rule_value[i])
        i += 1

    return answers


def table(family, midpoint_value, trapezoidal_value, integral_value, weights):
    table = PrettyTable()
    table.add_column("Functions", family)
    table.add_column("Midpoint Values", midpoint_value)
    table.add_column("Trapezoidal Values", trapezoidal_value)
    table.add_column("Integral Values", integral_value)
    table.add_column("Scalar Constants", weights)

    return table

def main():
    family = []
    print("Warning: Use ** for exponents instead of ^")
    while True:
        function = input("Please enter your next function or press q to finish: ")
        if function == "q" or function == "Q":
            break
        else:
            family.append(function)
    partitions = int(input("How many partitions? "))
    start = int(input("Interval start? "))
    end = int(input("Interval end? "))
    interval_bounds = [start, end]

    integral_value = integral_solution(interval_bounds, family)
    midpoint_value = midpoint_rule(interval_bounds, partitions, family)
    trapezoidal_value = trapezoidal_rule(interval_bounds, partitions*2, family)
    weights = find_weight(midpoint_value, trapezoidal_value, integral_value, family)
    print(table(family, midpoint_value, trapezoidal_value, integral_value, weights))

if __name__ == main():
    main()
