import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import solve

# Marine Food Web Model

# Define the ecological efficiencies
E = [0.92, 0.24, 0.20, 0.11, 0.13, 0.1]  
nodes = len(E);  # number of nodes in the food web

#Build the food-web model as a 6x6 matrix
Diet = np.array(
        [[0,0,0,0,0,0],
        [0,0,0,0,0,0],
        [0,0,0,0,0,0],
        [0,0,0,0,0,0],
        [0,0,0,0,0,0],
        [0,0,0,0,0,0]])

Diet[1,0] = 60
Diet[2,0] = 40
Diet[4,1] = 100
Diet[3,2] = 100
Diet[5,3] = 100
Diet[5,4] = 100

# columns should sum to 100 except the last one
Diet.sum(axis=0)

# Build the food-web model
print('Bottom-up Food-web model')
M = -0.01*Diet

for n in range(0,nodes): #update the diagonal
    M[n,n] += 1/E[n]
    
print(M)

# Define the vector of constants
C = np.array([0,0,0,0,0,0])
C[0]=1160 #initialize with primary productivity

print(C)

print('Solve for the fluxes')
#given the linear system x*C = M, solve for x
F = solve(M, C)

print(F)

objects = ('Phyto', 'Zoop', 'Benth', 'DemFish', 'PelFish', 'Pisc')
y_pos = np.arange(nodes)
 
plt.bar(y_pos, F, align='center', alpha=0.5)
plt.xticks(y_pos, objects)
plt.ylabel('Production')
plt.xlabel('Trophic compartment')
plt.title('Simple Food Web')
 
plt.show()