#!/usr/bin/env python

import math

from matplotlib import cm
import matplotlib.pyplot as plot
from mpl_toolkits.mplot3d import Axes3D
from numpy import arange

def combination(n, k):
    return math.factorial(n) / (math.factorial(k) * math.factorial(n - k))

la = 0.001
def reliability(p, m):
    r = []
    for t in range(4000):
        r_m = math.e ** (-la * t)
        r.append(0)
        for i in range(m, p + 1):
            r[-1] += combination(p, i) * (r_m ** i) * ((1 - r_m) ** (p - i))
    return r

def plot_reliability_threashold(p):
    r = []
    times = [arange(0, 4000, 1)] * (p - 2)
    threashold = []
    for m in range(2, p):
        r.append(reliability(p, m))
        threashold.append([m] * 4000)

    fig = plot.figure()
    ax = fig.gca(projection='3d')
    ax.plot_surface(times, threashold, r, cstride=100, cmap=cm.coolwarm)
    ax.set_yticks(range(2, p))
    ax.set_xlabel('Time')
    ax.set_ylabel('Voter Threashold')
    ax.set_zlabel('System Reliability')
    plot.title("Reliability of Self-Purging Redundancy ALU\n with {0} Modules over Time and Threashold".format(p))
    plot.savefig('reliability_{0}.pdf'.format(p), bbox_inches='tight')

def plot_reliability_modules():
    r = []
    times = [arange(0, 4000, 1)] * (9)
    num_modules = []
    for p in range(3, 12):
        r.append(reliability(p, 2))
        num_modules.append([p] * 4000)

    fig = plot.figure()
    ax = fig.gca(projection='3d')
    ax.plot_surface(times, num_modules, r, cstride=100, cmap=cm.coolwarm)
    ax.set_xlabel('Time')
    ax.set_ylabel('Number of Modules')
    ax.set_zlabel('System Reliability')
    plot.title("Reliability of Self-Purging Redundancy ALU\n with Voter Threashold of 2 over Time and Number of Modules".format(p))
    plot.savefig('reliability.pdf', bbox_inches='tight')

plot_reliability_threashold(5)
plot_reliability_threashold(11)

plot_reliability_modules()
