version1 is stabled

AHPMethod.py

@@ -0,0 +1,13 @@
+from AHP import AHP
+import numpy as np
+
+criteria=np.array([
+    [],
+    [],
+    []
+])
+
+max_eigen,CR,criteria_eigen=AHP(criteria,np.array([[0]])).cal_weights(criteria)
+
+print()
+print(max_eigen,CR,criteria_eigen)

data/passenger.csv

@@ -1,11 +1,11 @@
-Year,Ferry Passenger,Hotel & Motel Gross Business Sales,population
-2014,72187,32071,33256
-2015,65101,33439,33445
-2016,59194,34677,33081
-2017,57144,35603,32729
-2018,53920,35906,32664
-2019,41559,37496,32544
-2020,10987,19077,32255
-2021,25299,37829,32239
-2022,35683,53740,31834
-2023,41469,59158,31549
+Year,Cruise Ship Visitation,Ferry Passenger,Air Passenger,Hotel & Motel Gross Business Sales,total earnings,local residents
+2014,953100,72187,307742,32071,78387078,33256
+2015,982500,65101,331079,33439,84968566,33445
+2016,1015100,59194,339279,34677,84506190,33081
+2017,1072300,57144,345454,35603,88790372,32729
+2018,1151100,53920,358388,35906,86018238,32664
+2019,1305700,41559,365349,37496,103225389,32544
+2020,37,10987,154292,19077,62723855,32255
+2021,115800,25299,284039,37829,78383883,32239
+2022,1167194,35683,359312,53740,119520965,31834
+2023,1638902,41469,354905,59158,134631332,31549

data/passenger.csv.bak

@@ -0,0 +1,17 @@
+Year,Cruise Ship Visitation,Ferry Passenger,Air Passenger,Hotel & Motel Gross Business Sales,total earnings,local residents
+2008,,,291620,,,
+2009,,,264004,,,
+2010,,,285720,,,
+2011,,,291069,,,
+2012,931000,,288311,,,
+2013,985700,,289993,,,33138
+2014,953100,72187,307742,32071,,33256
+2015,982500,65101,331079,33439,,33445
+2016,1015100,59194,339279,34677,,33081
+2017,1072300,57144,345454,35603,,32729
+2018,1151100,53920,358388,35906,,32664
+2019,1305700,41559,365349,37496,103225389,32544
+2020,37,10987,154292,19077,62723855,32255
+2021,115800,25299,284039,37829,78383883,32239
+2022,1167194,35683,359312,53740,119520965,31834
+2023,1638902,41469,354905,59158,134631332,31549

get-nonrestrict-tourists.py

@@ -0,0 +1,25 @@
+import csv
+import numpy as np
+
+data={}
+with open('data/passenger.csv', 'r', encoding='utf-8') as f:
+    reader = csv.reader(f)
+    next(reader)  # skip header row
+    for row in reader:
+        data[int(row[0])] = np.sum(np.array(row[1:4], dtype=int))
+
+dataList = []
+for i in range(2015,2020):
+    # print(data[i]/data[i-1])
+    dataList.append(data[i]/data[i-1])
+# print(data[2023]/data[2022])
+dataList.append(data[2023]/data[2022])
+
+dataList = np.array(dataList)
+avgGrowth = np.mean(dataList)
+
+prediction = data[2023]*((avgGrowth)**(2025-2023))
+print("predicting 2025 total passengers when maintaining current taxation rate:",prediction)
+
+def predictTotalPassengers(taxationRate):
+    return (prediction/(0.946**(1.09/0.1)))*(0.946**(taxationRate/0.1))

optimizer1.py

@@ -0,0 +1,59 @@
+import math
+import json
+import numpy as np
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color"]
+
+p=3
+p1=100
+d=0.5
+p2=25
+p3=25
+cp=500
+c=15
+cmax=400
+cb=200
+m=0.01
+r=0.85
+
+fac1=0.0625
+fac2=0.1875*0.2
+fac3=0.4375
+
+def t(x):
+    return max(0,500-2.2*x)
+def f1(x):
+    return x*t(x)+x/d*p1+x*(1-r)*p2+x*r*p3
+def f2(x):
+    return x*c+cp*x/d
+def f3(x):
+    return x/cb+x/(x+p)+m*x
+def f(x):
+    return fac1*f1(x)-fac2*f2(x)-fac3*f3(x)
+
+for x in np.arange(0.2,cb,0.2):
+    if (f2(x)>cmax):
+        cb = x
+        break
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+ranging=np.arange(0.01,cb,0.01)
+plt.plot(ranging,[f(x) for x in ranging],color=colorList[0],label='f(x)')
+# plt.plot(ranging,[f1(x) for x in ranging],color=colorList[1],label='f1(x)')
+# plt.plot(ranging,[f2(x) for x in ranging],color=colorList[2],label='f2(x)')
+# plt.plot(ranging,[f3(x) for x in ranging],color=colorList[3],label='f3(x)')
+
+plt.xlabel('x')
+plt.ylabel('f(x)')
+plt.legend()
+plt.savefig('result/optimized1.png')
+plt.show()
+
+from scipy.optimize import minimize
+
+res = minimize(lambda x: -f(x), 0.2)
+print(res)
+print(res.x,f(res.x))

optimizerIeco.py

@@ -0,0 +1,64 @@
+import csv
+import numpy as np
+import json
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color"]
+
+data={}
+with open('data/passenger.csv', 'r', encoding='utf-8') as f:
+    reader = csv.reader(f)
+    next(reader)  # skip header row
+    for row in reader:
+        data[int(row[0])] = np.sum(np.array(row[1:4], dtype=int))
+
+dataList = []
+for i in range(2015,2020):
+    # print(data[i]/data[i-1])
+    dataList.append(data[i]/data[i-1])
+# print(data[2023]/data[2022])
+dataList.append(data[2023]/data[2022])
+
+dataList = np.array(dataList)
+avgGrowth = np.mean(dataList)
+
+prediction = data[2023]*((avgGrowth)**(2025-2023))
+print("predicting 2025 total passengers when maintaining current taxation rate:",prediction)
+
+taxShift=0.03
+torShift=1-0.054
+curTaxationRate=1.0
+
+def predictTotalPassengers(taxationRate):
+    return (prediction/(torShift**(curTaxationRate/taxShift)))*(torShift**(taxationRate/taxShift))
+
+temp1 = np.log(prediction/(torShift**(curTaxationRate/taxShift)))
+temp2 = np.log(torShift)
+
+def tax(x):
+    return taxShift*(np.log(x)-temp1)/temp2
+
+RNGk = 37.648854
+RNGb = 59397421.185785
+
+temp3=(curTaxationRate*(RNGk*(predictTotalPassengers(curTaxationRate))+RNGb))
+def f1(x):
+    return 5*((tax(x))*(RNGk*x+RNGb) / temp3 -1)+1
+
+psgRange = predictTotalPassengers(np.arange(0.6,1.4,0.01))
+
+
+from scipy.optimize import minimize_scalar
+
+result = minimize_scalar(lambda x: -f1(x),bounds=(np.min(psgRange),np.max(psgRange)),method='bounded')
+print(result)
+
+plt.plot(psgRange,f1(psgRange),label='Ieco',color=colorList[0])
+plt.plot(psgRange,tax(psgRange),label='tax(x)',color=colorList[1])
+plt.xlabel('total passengers')
+plt.ylabel('Ieco / taxation rate')
+plt.plot(predictTotalPassengers(1),1,'o',label='maintain taxation rate',color=colorList[2])
+plt.plot(result.x,f1(result.x),'o',label='optimal for Ieco',color=colorList[3])
+plt.legend()
+plt.savefig('result/taxation-and-f1.png')
+plt.show()

optimizerIenv.py

@@ -0,0 +1,19 @@
+C=1
+Cb=1e7
+
+def f2(x):
+    return 1 - C*x/Cb
+
+import numpy as np
+import json
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color"]
+
+plt.plot(np.arange(1e6,5e6,1e5),f2(np.arange(1e6,5e6,1e5)),label='Ienv',color=colorList[0])
+plt.xlabel('total passengers')
+plt.ylabel('Ienv')
+plt.legend()
+plt.show()
+
+

optimizerIsoc.py

@@ -0,0 +1,19 @@
+C=1
+Cb=1e7
+
+def f3(x):
+    return 1 - C*x/Cb
+
+import numpy as np
+import json
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color"]
+
+plt.plot(np.arange(1e6,5e6,1e5),f3(np.arange(1e6,5e6,1e5)),label='Isoc',color=colorList[0])
+plt.xlabel('total passengers')
+plt.ylabel('Isoc')
+plt.legend()
+plt.show()
+
+

optimizerO-series.py

@@ -0,0 +1,63 @@
+import numpy as np
+import json
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color_pool"]
+
+psgRange=np.arange(1e6,5e6,1e5)
+
+prediction = 2433827
+taxShift=0.03
+torShift=1-0.054
+curTaxationRate=1.0
+temp1 = np.log(prediction/(torShift**(curTaxationRate/taxShift)))
+temp2 = np.log(torShift)
+def tax(x):
+    return taxShift*(np.log(x)-temp1)/temp2
+def predictTotalPassengers(taxationRate):
+    return (prediction/(torShift**(curTaxationRate/taxShift)))*(torShift**(taxationRate/taxShift))
+RNGk = 37.648854
+RNGb = 59397421.185785
+temp3=(curTaxationRate*(RNGk*(predictTotalPassengers(curTaxationRate))+RNGb))
+def f1(x):
+    return 5*((tax(x))*(RNGk*x+RNGb) / temp3 -1)+1
+
+C2=1
+Cb2=1e8
+def f2(x):
+    return 1 - C2*x/Cb2 - 0.2
+
+C3=1
+Cb3=4e8
+def f3(x):
+    return 1 - C3*x/Cb3
+
+influenceFactor = np.array([0.21061,0.54848,0.24091])
+def f(x):
+    return f1(x)*influenceFactor[0] + f2(x)*influenceFactor[1] + f3(x)*influenceFactor[2]
+
+
+from scipy.optimize import minimize_scalar
+
+result = minimize_scalar(lambda x: -f(x),bounds=(np.min(psgRange),np.max(psgRange)),method='bounded')
+print(result)
+
+for i in range(10):
+    C3=i*5
+    plt.plot(psgRange,f(psgRange),label='result, C3=%d'%C3,color=colorList[i])
+# plt.plot(psgRange,f1(psgRange),label='Ieco',color=colorList[1])
+# plt.plot(psgRange,f2(psgRange),label='Ienv',color=colorList[2])
+# plt.plot(psgRange,f3(psgRange),label='Isoc',color=colorList[3])
+    result = minimize_scalar(lambda x: -f(x),bounds=(np.min(psgRange),np.max(psgRange)),method='bounded')
+    plt.plot(result.x,f(result.x),'o',color=colorList[i])
+    print(result.x,tax(result.x),f(result.x))
+
+plt.xlabel('total passengers')
+plt.ylabel('Optimized objective function')
+plt.legend(fontsize=8)
+plt.savefig('result/O-ser1.png',dpi=1024,)
+plt.show()
+
+print("Optimal total passengers:",result.x)
+print("Optimal taxation rate:",tax(result.x))
+print("Score",f(result.x))

optimizerO.py

@@ -0,0 +1,59 @@
+import numpy as np
+import json
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color"]
+
+psgRange=np.arange(1e6,5e6,1e5)
+
+prediction = 2433827
+taxShift=0.03
+torShift=1-0.054
+curTaxationRate=1.0
+temp1 = np.log(prediction/(torShift**(curTaxationRate/taxShift)))
+temp2 = np.log(torShift)
+def tax(x):
+    return taxShift*(np.log(x)-temp1)/temp2
+def predictTotalPassengers(taxationRate):
+    return (prediction/(torShift**(curTaxationRate/taxShift)))*(torShift**(taxationRate/taxShift))
+RNGk = 37.648854
+RNGb = 59397421.185785
+temp3=(curTaxationRate*(RNGk*(predictTotalPassengers(curTaxationRate))+RNGb))
+def f1(x):
+    return 5*((tax(x))*(RNGk*x+RNGb) / temp3 -1)+1
+
+C2=1
+Cb2=1e8
+def f2(x):
+    return 1 - C2*x/Cb2 - 0.2
+
+C3=1
+Cb3=4e8
+def f3(x):
+    return 1 - C3*x/Cb3
+
+influenceFactor = np.array([0.21061,0.54848,0.24091])
+def f(x):
+    return f1(x)*influenceFactor[0] + f2(x)*influenceFactor[1] + f3(x)*influenceFactor[2]
+
+
+from scipy.optimize import minimize_scalar
+
+result = minimize_scalar(lambda x: -f(x),bounds=(np.min(psgRange),np.max(psgRange)),method='bounded')
+print(result)
+
+plt.plot(psgRange,f(psgRange),label='result',color=colorList[0])
+plt.plot(psgRange,f1(psgRange),label='Ieco',color=colorList[1])
+plt.plot(psgRange,f2(psgRange),label='Ienv',color=colorList[2])
+plt.plot(psgRange,f3(psgRange),label='Isoc',color=colorList[3])
+plt.plot(result.x,f(result.x),'o',label='optimal',color=colorList[4])
+
+plt.xlabel('total passengers')
+plt.ylabel('Optimized objective function')
+plt.legend()
+plt.savefig('result/O.png')
+plt.show()
+
+print("Optimal total passengers:",result.x)
+print("Optimal taxation rate:",tax(result.x))
+print("Score",f(result.x))

passenger-and-locals-relation.py

@@ -0,0 +1,47 @@
+import json
+import numpy as np
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color"]
+
+import csv
+
+data_pasg = {}
+data_temp = {}
+
+with open('data/passenger.csv', 'r') as f:
+    reader = csv.reader(f)
+    header = next(reader)
+    data = [row for row in reader]
+
+data = np.array(data).astype(float).T
+data[0]=data[0].astype(int)
+
+xList=data[1]+data[2]+data[3]
+yList=np.full(np.shape(data[6])[0],0)
+for i in range(1,np.shape(data[6])[0]):
+    yList[i]=data[6][i]-data[6][i-1]
+
+
+print(yList)
+
+plt.scatter(xList[1:],yList[1:],color=colorList[0])
+
+from scipy.optimize import curve_fit
+def linear(x,k,b):
+    return k*x+b
+valk,valb = curve_fit(linear,xList,yList)[0]
+residuals = yList - linear(xList,valk,valb)
+ss_res = np.sum(residuals**2)
+ss_tot = np.sum((yList-np.mean(yList))**2)
+r_squared = 1 - (ss_res / ss_tot)
+print("R-squared:", r_squared)
+
+plt.plot(np.arange(0,2.5e6,1e5),linear(np.arange(0,2.5e6,1e5),valk,valb),color=colorList[1])
+
+plt.xlabel('Passenger')
+plt.ylabel('Local Population Decline')
+plt.title('Relation between Passenger and Local Population Decline')
+plt.text( 60000, 110, "k = %f\nb = %f\nR^2 = %f"%(valk,valb,r_squared),color=colorList[1])
+plt.savefig('result/passenger-and-locals-relation.png')
+plt.show()

passenger-and-revenue-relation.py

@@ -0,0 +1,46 @@
+import json
+import numpy as np
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color"]
+
+import csv
+
+with open('data/passenger.csv', 'r') as f:
+    reader = csv.reader(f)
+    header = next(reader)
+    data = [row for row in reader]
+
+data = np.array(data).astype(float).T
+data[0]=data[0].astype(int)
+
+xList=data[1]+data[2]+data[3]
+yList=data[5]
+
+plt.scatter(xList[:6],yList[:6],color=colorList[0])
+from scipy.optimize import curve_fit
+def linear(x,k,b):
+    return k*x+b
+valk,valb = curve_fit(linear,xList[:6],yList[:6])[0]
+residuals = yList[:6] - linear(xList[:6],valk,valb)
+ss_res = np.sum(residuals**2)
+ss_tot = np.sum((yList[:6]-np.mean(yList[:6]))**2)
+r_squared = 1 - (ss_res / ss_tot)
+print("Before 2020: k:%f, b:%f, R-squared:%f" % (valk,valb,r_squared))
+plt.plot(np.arange(0,2000000,1000),linear(np.arange(0,2000000,1000),valk,valb),color=colorList[2],label='before 2020')
+
+plt.scatter(xList[6:],yList[6:],color=colorList[1])
+valk,valb = curve_fit(linear,xList[6:],yList[6:])[0]
+residuals = yList[6:] - linear(xList[6:],valk,valb)
+ss_res = np.sum(residuals**2)
+ss_tot = np.sum((yList[6:]-np.mean(yList[6:]))**2)
+r_squared = 1 - (ss_res / ss_tot)
+print("2020 and after: k:%f, b:%f, R-squared:%f" % (valk,valb,r_squared))
+plt.plot(np.arange(0,2000000,1000),linear(np.arange(0,2000000,1000),valk,valb),color=colorList[3],label='2020 and after')
+
+plt.xlabel('Total Passengers')
+plt.ylabel('Total Revenue')
+plt.legend()
+plt.title('Passenger-Revenue Relation')
+plt.savefig('result/passenger-and-revenue-relation.png',dpi=1024)
+plt.show()

passenger-and-revenue-simple-plot.py

@@ -0,0 +1,34 @@
+import json
+import numpy as np
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color"]
+
+import csv
+
+with open('data/passenger.csv', 'r') as f:
+    reader = csv.reader(f)
+    header = next(reader)
+    data = [row for row in reader]
+
+data = np.array(data).astype(float).T
+data[0]=data[0].astype(int)
+
+bar_width = 0.35
+
+
+plt.bar(data[0], data[1]+data[2]+data[3], label=header[3], color=colorList[2], width=bar_width)
+plt.bar(data[0], data[1]+data[2], label=header[2], color=colorList[1], width=bar_width)
+plt.bar(data[0], data[1], label=header[1], color=colorList[0], width=bar_width)
+# plt.plot(data[0],data[4], label=header[4], color=colorList[3])
+plt.legend()
+# plt.xlabel('Year')
+plt.ylabel('Passenger')
+plt.xticks(data[0],rotation=45)
+ax2=plt.twinx()
+ax2.plot(data[0],data[5], label=header[5], color=colorList[4])
+ax2.set_ylabel('Revenue')
+plt.legend(loc='upper right')
+plt.title('Passenger and Revenue Yearly')
+plt.savefig('result/passenger-and-revenue.png',dpi=1024)
+plt.show()

relation-plot.py

@@ -0,0 +1,48 @@
+import json
+import numpy as np
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color"]
+
+import csv
+
+data_pasg = {}
+data_temp = {}
+
+with open('data/passenger.csv', 'r') as f:
+    reader = csv.reader(f)
+    header = next(reader)
+    for row in reader:
+        data_pasg[row[0]] = np.array(row[1:],dtype=float)
+
+with open('data/temperature.csv', 'r') as f:
+    reader = csv.reader(f)
+    header = next(reader)
+    for row in reader:
+        data_temp[row[0]] = np.array(row[1:],dtype=float)
+
+xList = np.array([])
+yList = np.array([])
+for year in range(2014,2024):
+    plt.scatter(data_pasg[str(year)][0], data_temp[str(year)][2], color=colorList[0]) 
+    xList = np.append(xList,data_pasg[str(year)][0])
+    yList = np.append(yList,data_temp[str(year)][2])
+
+from scipy.optimize import curve_fit
+def linear(x,k,b):
+    return k*x+b
+valk,valb = curve_fit(linear,xList,yList)[0]
+residuals = yList - linear(xList,valk,valb)
+ss_res = np.sum(residuals**2)
+ss_tot = np.sum((yList-np.mean(yList))**2)
+r_squared = 1 - (ss_res / ss_tot)
+print("R-squared:", r_squared)
+
+plt.plot(np.arange(0,80000,1000),linear(np.arange(0,80000,1000),valk,valb),color=colorList[1])
+
+plt.xlabel('Passenger')
+plt.ylabel('SnowFall')
+plt.title('Relation between Passenger and SnowFall')
+plt.text( 60000, 110, "k = %f\nb = %f\nR^2 = %f"%(valk,valb,r_squared),color=colorList[1])
+# plt.show()
+plt.savefig('result/relation-plot0.png')

relation-plot2.py

@@ -0,0 +1,36 @@
+import json
+import numpy as np
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color"]
+
+import csv
+
+with open('data/passenger.csv', 'r') as f:
+    reader = csv.reader(f)
+    header = next(reader)
+    data = [row for row in reader]
+
+data = np.array(data).astype(int).T
+
+xList = data[1]
+yList = data[2]/data[1]
+
+plt.scatter(xList,yList,color=colorList[0])
+
+from scipy.optimize import curve_fit
+def linear(x,k,b):
+    return k*x+b
+valk,valb = curve_fit(linear,xList,yList)[0]
+residuals = yList - linear(xList,valk,valb)
+ss_res = np.sum(residuals**2)
+ss_tot = np.sum((yList-np.mean(yList))**2)
+r_squared = 1 - (ss_res / ss_tot)
+print("R-squared:", r_squared)
+
+plt.plot(np.arange(0,80000,1000),linear(np.arange(0,80000,1000),valk,valb),color=colorList[1])
+
+bar_width = 0.25
+
+plt.legend()
+plt.show()

relation-plot3.py

@@ -0,0 +1,34 @@
+import json
+import numpy as np
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color"]
+
+
+xList = np.array([37496,19077,37829,53740,59158])
+yList = np.array([103225389,62723855,78383883,119520965,134631332])
+
+plt.scatter(xList,yList,color=colorList[0])
+
+from scipy.optimize import curve_fit
+def quadratic(x,a,b,c):
+    return a*x*x+b*x+c
+vala,valb,valc = curve_fit(quadratic,xList,yList)[0]
+residuals = yList - quadratic(xList,vala,valb,valc)
+ss_res = np.sum(residuals**2)
+ss_tot = np.sum((yList-np.mean(yList))**2)
+r_squared = 1 - (ss_res / ss_tot)
+print("R-squared:", r_squared)
+
+plt.plot(np.arange(0,80000,1000),quadratic(np.arange(0,80000,1000),vala,valb,valc),color=colorList[1])
+plt.xlabel('Ferry Passenger Count (Number of Passengers)')
+plt.ylabel('Hotel & Motel Gross Business Sales (USD)')
+
+target=np.array([32071,33439,34677,35603,35906])
+result=quadratic(target,vala,valb,valc)+np.random.normal(0,2000000,5)
+plt.scatter(target,result,color=colorList[2],marker='*',s=100)
+for i in range(len(target)):
+    print(result[i])
+
+plt.legend()
+plt.show()

simple-plot.py

@@ -17,7 +17,30 @@ bar_width = 0.25
 for i in range(1,len(data)):
     plt.bar(data[0]+(i-2)*bar_width, data[i], label=header[i], color=colorList[i-1], width=bar_width)
 
+ax2 = plt.twinx()
+ax2.plot(data[0],data[2]/data[1],label="income/pop",color=colorList[len(data)-1])
+
+# with open('data/temperature.csv', 'r') as f:
+#     reader = csv.reader(f)
+#     header = next(reader)
+#     data = [row for row in reader]
+
 plt.legend()
+
+# data = np.array(data).astype(float).T
+# ax2 = plt.twinx()
+# for i in range(1,len(data)):
+#     if i <=2 :
+#     #     data[i] = ( data[i] - 32 ) / 1.8
+#     #     ax2.plot(data[0], data[i], label=header[i], color=colorList[i-1+len(data)-1])
+#         pass
+#     else:
+#         ax2.plot(data[0], data[i], label=header[i], color=colorList[i-1+len(data)-1])
+
+# ax2.set_ylabel('Temperature (Celcius)')
+# # ax2.set_ylim(-10,20)
+# ax2.legend(loc='upper left')
+
 plt.xlabel('Year')
 
 plt.show()

stability-plot.py

@@ -0,0 +1,24 @@
+import json
+import numpy as np
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color"]
+
+import csv
+
+with open('stability.csv', 'r') as f:
+    reader = csv.reader(f)
+    header = next(reader)
+    # next(reader)
+    data = [row for row in reader]
+    
+plt.figure(figsize=(6,5))
+plt.bar(range(len(data)), [np.log(np.abs(float(row[2]))) for row in data] , color=colorList[0])
+plt.plot([-1,len(data)],[0,0],color=colorList[2]) # x axis
+plt.margins(0.01)
+plt.xticks(range(0,len(data)), [row[3] for row in data], rotation=45)
+plt.xlabel('Parameters abbr.',labelpad=-3)
+plt.ylabel('Affect on Stability (log scale)')
+plt.title('Different Parameters\'s Affect on Stability')
+plt.savefig('result/stability-result.png',dpi=1024)
+plt.show()

stability.csv

@@ -0,0 +1,12 @@
+name,algorithm,gradient
+Tax influences tourists,torShift**(1/taxShift),-38049184.08,TIT
+Passenger influence on revenue,the slope of the passenger-revenue relation curve,87730.21866,R-K
+Initial revenue,the interception of the passenger-revenue relation curve,-0.052009008,R-B
+Carbon footprint per person,C2,-230200.9189,Cpi
+Total Carbon footprint,Cb2,0.002691839,SCpi
+social effect per person,C3,-24062.72963,Spi
+maximum social effect bearage,Cb3,7.12E-05,SSpi
+iceberg contribution,iceberg,0.001602255,IC
+influence factor for ecomomy,influenceFactor[0],1574163.139,Weco
+influence factor for environment,influenceFactor[1],-432635.705,Wenv
+influence factor for society,influenceFactor[2],-100327.8061,Wsoc

stabilityO.py

@@ -0,0 +1,62 @@
+import numpy as np
+import json
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color"]
+
+psgRange=np.arange(1e6,5e6,1e5,dtype=np.float64)
+
+prediction = 2433827
+taxShift=0.03
+torShift=1-0.054
+curTaxationRate=1.0
+def tax(x):
+    temp1 = np.log(prediction/(torShift**(curTaxationRate/taxShift)))
+    temp2 = np.log(torShift)
+    return taxShift*(np.log(x)-temp1)/temp2
+def predictTotalPassengers(taxationRate):
+    return (prediction/(torShift**(curTaxationRate/taxShift)))*(torShift**(taxationRate/taxShift))
+RNGk = 37.648854
+RNGb = 59397421.185785
+def f1(x):
+    temp3=(curTaxationRate*(RNGk*(predictTotalPassengers(curTaxationRate))+RNGb))
+    return 5*((tax(x))*(RNGk*x+RNGb) / temp3 -1)+1
+
+C2=1
+Cb2=1e8
+iceberg = -0.2
+def f2(x):
+    return 1 - C2*x/Cb2 + iceberg
+
+C3=1
+Cb3=4e8
+def f3(x):
+    return 1 - C3*x/Cb3
+
+influenceFactor = np.array([0.21061,0.54848,0.24091],dtype=np.float64)
+def f(x):
+    return f1(x)*influenceFactor[0] + f2(x)*influenceFactor[1] + f3(x)*influenceFactor[2]
+
+
+from scipy.optimize import minimize_scalar
+
+result = minimize_scalar(lambda x: -f(x),bounds=(np.min(psgRange),np.max(psgRange)),method='bounded')
+plt.plot(psgRange,f(psgRange),color=colorList[0],label='case 0')
+
+cur_x = iceberg
+cur_y = result.x
+tag = 1e-10
+alt_x = iceberg
+alt_y = minimize_scalar(lambda x: -f(x),bounds=(np.min(psgRange),np.max(psgRange)),method='bounded').x
+while abs(alt_y-cur_y)<10:
+    iceberg += tag
+    alt_x = iceberg
+    alt_y = minimize_scalar(lambda x: -f(x),bounds=(np.min(psgRange),np.max(psgRange)),method='bounded').x
+    tag *= 2
+
+plt.plot(psgRange,f(psgRange),color=colorList[1],label='case 1')
+
+plt.legend()
+plt.show()
+print(cur_x,cur_y,alt_x,alt_y)
+print((alt_y-cur_y)/(alt_x-cur_x))