#!/bin/python from scipy.integrate import odeint import numpy as np from PyQt4.QtCore import * from PyQt4.QtGui import * from pyqtgraph.widgets.PlotWidget import PlotWidget from pyqtgraph.graphicsItems.InfiniteLine import InfiniteLine from pyqtgraph.graphicsItems.ScatterPlotItem import ScatterPlotItem from pyqtgraph import mkPen, mkBrush, setConfigOption from pyqtgraph.parametertree import ParameterTree, Parameter import pyqtgraph as pg from pylab import plot, show, xlabel, ylabel, quiver t_init, t_step, t_end = 0, 0.1, 10000 t = t_init t_values = np.arange(t_init, t_end, t_step) class MyWidget(QWidget): def __init__(self): super(MyWidget, self).__init__() self.parameter_tree = ParameterTree() self.parameter_tree.setMaximumWidth(300) self.k1, self.k2, self.k3, self.k4 = 0.4, 0.6, 0.1, 0.005 self.alpha, self.beta = 0.01, 0.1 self.Kf, self.Kp = 0, 0 self.d1, self.d2, self.d3 = 0.3, 0.2, 0.044 self.F, self.P, self.M = 100, 50, 2 self.parameter = Parameter(name='Parameters', children=[ { 'name':'Initial Conditions', 'type':'group', 'children': [{ 'name':'F', 'type':'float', 'value':self.F }, { 'name':'P', 'type':'float', 'value':self.P }, { 'name':'M', 'type':'float', 'value':self.M }] }, { 'name':'K', 'type':'group', 'children': [{ 'name':'k1', 'type':'float', 'value':self.k1 }, { 'name':'k2', 'type':'float', 'value':self.k2 }, { 'name':'k3', 'type':'float', 'value':self.k3 }, { 'name':'k4', 'type':'float', 'value':self.k4 }] }, { 'name':'D', 'type':'group', 'children': [{ 'name':'d1', 'type':'float', 'value':self.d1 }, { 'name':'d2', 'type':'float', 'value':self.d2 }, { 'name':'d3', 'type':'float', 'value':self.d3 }] }, { 'name':'Other paramters', 'type':'group', 'children': [{ 'name':'alpha', 'type':'float', 'value':self.alpha }, { 'name':'beta', 'type':'float', 'value':self.beta }, { 'name':'Kf', 'type':'float', 'value':self.Kf }, { 'name':'Kp', 'type':'float', 'value':self.Kp }, ] }] ) self.parameter_tree.setParameters(self.parameter) self.plot_widget = PlotWidget() self.plt = self.plot_widget.getPlotItem() l = self.plt.addLegend(size=(100, 100), offset=(500, 30)) self.plt.showGrid(True, True) res = self.solve() self.f_plt = self.plt.plot(t_values, res[0], name='Flowers', antialias=True) self.f_plt.setPen(mkPen('#EE02FF', width=3)) self.p_plt = self.plt.plot(t_values, res[1], name='Pollinators', antialias=True) self.p_plt.setPen(mkPen('#FFD800', width=3)) self.m_plt = self.plt.plot(t_values, res[2], name='Mantis', antialias=True) self.m_plt.setPen(mkPen('#FF0000', width=3)) print('Flowers {0}\nPollinators {1}\nMantis {2}'.format(res[0][-1], res[1][-1], res[2][-1])) self.fp_plt = PlotWidget() fp_plt = self.fp_plt.getPlotItem() fp_plt.showGrid(True, True) l = fp_plt.addLegend(size=(100, 100), offset=(500, 30)) self.fp_f_plt = fp_plt.plot([], [], name='Flowers', antialias=True) # , symbol='+') self.fp_f_plt.setPen(mkPen('#EE02FF', width=3)) self.fp_p_plt = fp_plt.plot([], [], name='Pollinators', antialias=True) # , symbol='s') self.fp_p_plt.setPen(mkPen('#FFD800', width=3)) self.fp_m_plt = fp_plt.plot([], [], name='Mantis', antialias=True) # , symbol='t') self.fp_m_plt.setPen(mkPen('#FF0000', width=3)) wid = QWidget() l2 = QVBoxLayout() self.min_value = QDoubleSpinBox() self.max_value = QDoubleSpinBox() self.parameter_name = QLineEdit('d3') run_button = QPushButton('Run') l3 = QHBoxLayout() l3.addWidget(QLabel('Parameter: ')) l3.addWidget(self.parameter_name) l3.addWidget(QLabel('Min value: ')) l3.addWidget(self.min_value) l3.addWidget(QLabel('Max value: ')) l3.addWidget(self.max_value) l3.addWidget(run_button) l2.addLayout(l3) l2.addWidget(self.fp_plt) wid.setLayout(l2) self.fpa_plt = PlotWidget() fpa_plt = self.fpa_plt.getPlotItem() fpa_plt.showGrid(True, True) yaxis = fpa_plt.getAxis('left') yaxis.setTicks([[(1, '0 0 0'), (2, 'CO'), (3, 'NM -'), (4, ('NM +'))]]) line1 = InfiniteLine(QPointF(0, 1), angle=0, pen=mkPen('#f00', width=2)) line2 = InfiniteLine(QPointF(0, 2), angle=0, pen=mkPen('#f00', width=2)) line3 = InfiniteLine(QPointF(0, 3), angle=0, pen=mkPen('#f00', width=2)) line4 = InfiniteLine(QPointF(0, 4), angle=0, pen=mkPen('#f00', width=2)) fpa_plt.addItem(line1) fpa_plt.addItem(line2) fpa_plt.addItem(line3) fpa_plt.addItem(line4) self.sc1 = ScatterPlotItem(antialias=True, brush=mkBrush('#00f'), size=20) self.sc2 = ScatterPlotItem(antialias=True, brush=mkBrush('#ff0'), size=20) self.index = 0 fpa_plt.addItem(self.sc1) fpa_plt.addItem(self.sc2) wid2 = QWidget() l4 = QVBoxLayout() l5 = QHBoxLayout() self.min_fxa = QDoubleSpinBox() self.max_fxa = QDoubleSpinBox() self.param_fxa = QLineEdit('d3') run_button_2 = QPushButton('Run') l5.addWidget(QLabel('Parameter: ')) l5.addWidget(self.param_fxa) l5.addWidget(QLabel('Min value: ')) l5.addWidget(self.min_fxa) l5.addWidget(QLabel('Max value: ')) l5.addWidget(self.max_fxa) l5.addWidget(run_button_2) l4.addLayout(l5) l4.addWidget(self.fpa_plt) wid2.setLayout(l4) lay = QHBoxLayout() self.tab = QTabWidget() l6 = QVBoxLayout() l6.addWidget(self.parameter_tree) phase_button = QPushButton('Plot phase diagram') l6.addWidget(phase_button) lay.addLayout(l6) lay.addWidget(self.tab) self.tab.addTab(self.plot_widget, 'Simulation') self.tab.addTab(wid, 'Fixed Points') self.tab.addTab(wid2, 'Fixed Points Analysis') self.setLayout(lay) self.parameter.sigTreeStateChanged.connect(self.upd) self.update_fp() run_button.clicked.connect(self.run_a_lot) run_button_2.clicked.connect(self.run_a_lot_of_points) phase_button.clicked.connect(self.plot_phase) def run_a_lot(self): interval = self.min_value.value(), self.max_value.value() values = np.linspace(interval[0], interval[1], 200) parameter_name = self.parameter_name.text() old_val = eval(str('self.' + parameter_name)) f, p, m = [], [], [] self.fp_plt.getPlotItem().setLabels(bottom=str(parameter_name), left='Fixed Points') for i in values: exec (str('self.' + parameter_name + '=' + str(i))) res = self.solve() f.append(res[0][-1]) p.append(res[1][-1]) m.append(res[2][-1]) l = len(f) print('Finished simulation for {0} = {1}'.format(parameter_name, eval(str('self.' + parameter_name)))) self.fp_f_plt.setData(values, f) self.fp_p_plt.setData(values, p) self.fp_m_plt.setData(values, m) self.fp_plt.update() exec (str('self.' + parameter_name + '=' + str(old_val))) def run_a_lot_of_points(self): interval = self.min_fxa.value(), self.max_fxa.value() values = np.linspace(interval[0], interval[1], 100) parameter_name = self.param_fxa.text() old_val = eval(str('self.' + parameter_name)) s1x, s1y, s2x, s2y = [], [], [], [] self.sc1.clear() self.sc2.clear() self.index = 0 for i in values: exec (str('self.' + parameter_name + '=' + str(i))) self.update_fp() res = self.solve() print('Finished simulation for {0} = {1}'.format(parameter_name, eval(str('self.' + parameter_name)))) f, p, m = res[0][-1], res[1][-1], res[2][-1] eps = 5 if np.abs(f - self.fp0[0]) < eps and np.abs(p - self.fp0[1]) < eps and np.abs(m - self.fp0[2]) < eps: val = 1 elif np.abs(f - self.fp1[0]) < eps and np.abs(p - self.fp1[1]) < eps and np.abs(m - self.fp1[2]) < eps: val = 2 elif np.abs(f - self.fp2[0]) < eps and np.abs(p - self.fp2[1]) < eps and np.abs(m - self.fp2[2]) < eps: val = 3 elif np.abs(f - self.fp3[0]) < eps and np.abs(p - self.fp3[1]) < eps and np.abs(m - self.fp3[2]) < eps: val = 4 else: val = 5 if ((self.k4 / self.d3) - (self.k1 / self.d1)) < self.beta \ and (self.k4 / self.d3) > self.beta \ and (self.alpha * self.k1 * self.d2) / (self.d1 * self.k2) \ < ((self.k4 / self.d3) - self.beta) * \ (self.beta + (self.k1 / self.d1) - (self.k4 / self.d3)): s2x.append(self.index) s2y.append(val) else: s1x.append(self.index) s1y.append(val) show = (4 * self.alpha * self.d1 * self.d2) / (self.k1 * self.k2) self.fpa_plt.setXRange(0, self.index, 0.5, True) self.index += 1 self.sc1.setData(s1x, s1y) self.sc2.setData(s2x, s2y) def upd(self, *args): exec ('self.' + args[1][0][0].name() + '=' + str(args[1][0][2])) self.update_fp() res = self.solve() print('Flowers {0}\nPollinators {1}\nMantis {2}'.format(res[0][-1], res[1][-1], res[2][-1])) if self.tab.currentIndex() == 0: self.f_plt.setData(t_values, res[0]) self.p_plt.setData(t_values, res[1]) self.m_plt.setData(t_values, res[2]) self.plot_widget.update() def update_fp(self): P0 = F0 = M0 = 0 P1 = self.d3 / (self.k4 - self.d3 * self.beta) F1 = (self.k1 * P1 - self.d1) / (self.alpha * P1 * self.d1) M1 = ((self.k2 * F1) / (1 + self.alpha * P1 * F1) - self.d2) * ((1 + self.beta * P1) / self.k3) M2 = 0 P2 = (self.k1 * self.k2 - np.sqrt(self.k1 * self.k2 * (-4 * self.alpha * self.d1 * self.d2 + self.k1 * self.k2))) \ / (2 * self.alpha * self.d2 * self.k1) F2 = (self.k1 * self.k2 - np.sqrt(self.k1 * self.k2 * (-4 * self.alpha * self.d1 * self.d2 + self.k1 * self.k2))) \ / (2 * self.alpha * self.d1 * self.k2) M3 = 0 P3 = (self.k1 * self.k2 + np.sqrt(self.k1 * self.k2 * (-4 * self.alpha * self.d1 * self.d2 + self.k1 * self.k2))) \ / (2 * self.alpha * self.d2 * self.k1) F3 = (self.k1 * self.k2 + np.sqrt(self.k1 * self.k2 * (-4 * self.alpha * self.d1 * self.d2 + self.k1 * self.k2))) \ / (2 * self.alpha * self.d1 * self.k2) self.fp0, self.fp1, self.fp2, self.fp3 = (F0, P0, M0), (F1, P1, M1), (F2, P2, M2), (F3, P3, M3) print('F0 = {0}\nP0 = {1}\nM0 = {2}\n' .format(F0, P0, M0)) print('F1 = {0}\nP1 = {1}\nM1 = {2}\n' .format(F1, P1, M1)) print('F2 = {0}\nP2 = {1}\nM2 = {2}\n' .format(F2, P2, M2)) print('F3 = {0}\nP3 = {1}\nM3 = {2}\n' .format(F3, P3, M3)) print() def solve(self): res = odeint(self.func(), (self.F, self.P, self.M), t_values).T return res def func(self): flower = lambda t, f, p, m:(self.k1 * p * f) / (1 + self.alpha * f * p) \ - self.d1 * f \ - self.Kf * f * f polli = lambda t, f, p, m:(self.k2 * p * f) / (1 + self.alpha * f * p) \ - (self.k3 * m * p) / (1 + self.beta * p) \ - self.d2 * p \ - self.Kp * p * p mantis = lambda t, f, p, m:((self.k4 * m * p)) / (1 + self.beta * p) \ - self.d3 * m return lambda y, t:(flower(t, *y), polli(t, *y), mantis(t, *y)) def plot_phase(self): # Equations without the mantis (simpler model for the phase space) flower = lambda t, f, p:(self.k1 * p * f) / (1 + self.alpha * f * p) \ - self.d1 * f polli = lambda t, f, p:(self.k2 * p * f) / (1 + self.alpha * f * p) \ - self.d2 * p f = lambda y, t:np.array([flower(t, *y), polli(t, *y)]) # res = odeint(f, (self.F, self.P), t_values).T # plot(res[0], res[1]) res = odeint(f, (20, 20), t_values).T plot(res[0], res[1]) res = odeint(f, (150, 20), t_values).T plot(res[0], res[1]) res = odeint(f, (150, 350), t_values).T plot(res[0], res[1]) res = odeint(f, (200, 20), t_values).T plot(res[0], res[1]) res = odeint(f, (150, 20), t_values).T plot(res[0], res[1]) res = odeint(f, (150, 0), t_values).T plot(res[0], res[1]) xlabel('Flowers') ylabel('Pollinators') res = odeint(f, (0, 50), t_values).T plot(res[0], res[1]) R, C = np.meshgrid(np.arange(-30, 200, 30), np.arange(-30, 330, 10)) dy = f(np.array([R, C]), 0) plot([self.fp0[0]], [self.fp0[1]], 'o') plot([self.fp2[0]], [self.fp2[1]], 'o') plot([self.fp3[0]], [self.fp3[1]], 'o') quiver(R, C, dy[0, :], dy[1, :], scale_units='xy', angles='xy') show() xlabel('Flowers') ylabel('Pollinators') res = odeint(f, (0.5, 0.3), t_values).T plot(res[0], res[1]) res = odeint(f, (0.3, 0.8), t_values).T plot(res[0], res[1]) res = odeint(f, (.1, 1), t_values).T plot(res[0], res[1]) res = odeint(f, (0.3, 0.7), t_values).T plot(res[0], res[1]) res = odeint(f, (0.2, 0.75), t_values).T plot(res[0], res[1]) res = odeint(f, (0.2, 0.7), t_values).T plot(res[0], res[1]) res = odeint(f, (0.334170, 0.75188), t_values).T plot(res[0], res[1]) plot([self.fp0[0]], [self.fp0[1]], 'o') plot([self.fp2[0]], [self.fp2[1]], 'o') R, C = np.meshgrid(np.arange(-1, 2, .2), np.arange(-1, 2, .2)) dy = f(np.array([R, C]), 0) quiver(R, C, dy[0, :], dy[1, :], scale_units='xy', angles='xy') show() app = QApplication([]) wid = MyWidget() wid.show() app.exec_()