Commit a67a8401933ce41a867fed5558beebf03d6112b5

Authored by Eduardo Eloy
1 parent 33831508
Exists in master

python script to plot 3d pareto frontiers

Showing 1 changed file with 110 additions and 0 deletions   Show diff stats
MODFIRE-Prototype/plot3d.py 0 → 100644
... ... @@ -0,0 +1,110 @@
  1 +import numpy as np
  2 +import matplotlib.pyplot as plt
  3 +
  4 +
  5 +#val = input("1:WoodYield/SoilLoss \n 2:WoodYield/RiskPercentile \n 3:SoilLoss/RiskPercentile \n 4:All three
  6 +
  7 +f=open("paretoWSR.csv","r")
  8 +#lines=sorted(list(map(str, f.readlines())))
  9 +
  10 +f2 =open("nonParetoWSR.csv","r")
  11 +#lines2 = sorted(list(map(str, f2.readlines())))
  12 +
  13 +lines=f.readlines()
  14 +lines2 = f2.readlines()
  15 +
  16 +
  17 +paretoPoints = []
  18 +dominatedPoints = []
  19 +x = []
  20 +y = []
  21 +z = []
  22 +lst = []
  23 +scatterx = []
  24 +scattery = []
  25 +
  26 +for i in lines:
  27 + result = i.split(',')
  28 + x.append(int(result[0]))
  29 + y.append(int(result[1]))
  30 + paretoPoints.append([int(result[0]),int(result[1]),int(result[2])])
  31 +f.close()
  32 +
  33 +for i in lines2:
  34 + result = i.split(',')
  35 + x.append(int(result[0]))
  36 + y.append(int(result[1]))
  37 + z.append(int(result[2]))
  38 + if([int(result[0]),int(result[1]),int(result[2])] not in paretoPoints):
  39 + dominatedPoints.append([int(result[0]),int(result[1]),int(result[2])])
  40 +f2.close()
  41 +
  42 +#print(dominatedPoints)
  43 +
  44 +def simple_cull(inputPoints, dominates):
  45 + paretoPoints = set()
  46 + candidateRowNr = 0
  47 + dominatedPoints = set()
  48 + while True:
  49 + candidateRow = inputPoints[candidateRowNr]
  50 + inputPoints.remove(candidateRow)
  51 + rowNr = 0
  52 + nonDominated = True
  53 + while len(inputPoints) != 0 and rowNr < len(inputPoints):
  54 + row = inputPoints[rowNr]
  55 + if dominates(candidateRow, row):
  56 + # If it is worse on all features remove the row from the array
  57 + inputPoints.remove(row)
  58 + dominatedPoints.add(tuple(row))
  59 + elif dominates(row, candidateRow):
  60 + nonDominated = False
  61 + dominatedPoints.add(tuple(candidateRow))
  62 + rowNr += 1
  63 + else:
  64 + rowNr += 1
  65 +
  66 + if nonDominated:
  67 + # add the non-dominated point to the Pareto frontier
  68 + paretoPoints.add(tuple(candidateRow))
  69 +
  70 + if len(inputPoints) == 0:
  71 + break
  72 + return paretoPoints, dominatedPoints
  73 +
  74 +def dominates(row, candidateRow):
  75 + return sum([row[x] >= candidateRow[x] for x in range(len(row))]) == len(row)
  76 +
  77 +import random
  78 +inputPoints = [[random.randint(70,100) for i in range(3)] for j in range(500)]
  79 +#paretoPoints, dominatedPoints = simple_cull(inputPoints, dominates)
  80 +
  81 +#print()
  82 +
  83 +#print(paretoPoints)
  84 +
  85 +import matplotlib.pyplot as plt
  86 +from mpl_toolkits.mplot3d import Axes3D
  87 +
  88 +fig = plt.figure()
  89 +ax = fig.add_subplot(111, projection='3d')
  90 +dp = np.array(dominatedPoints)
  91 +pp = np.array(paretoPoints)
  92 +#print(dp)
  93 +ax.scatter(dp[:,0],dp[:,1],dp[:,2])
  94 +ax.scatter(pp[:,0],pp[:,1],pp[:,2],color='red')
  95 +
  96 +import matplotlib.tri as mtri
  97 +triang = mtri.Triangulation(pp[:,0],pp[:,1])
  98 +ax.plot_trisurf(triang,pp[:,2],color='red')
  99 +
  100 +ax.set_xlabel("Wood Yield",linespacing=5)
  101 +ax.set_ylabel("Soil Loss",linespacing=5)
  102 +ax.set_zlabel("Fire Risk Protection",linespacing=5)
  103 +ax.xaxis.labelpad=7
  104 +ax.yaxis.labelpad=7
  105 +ax.zaxis.labelpad=7
  106 +#plt.xlabel("Wood Yield")
  107 +#plt.ylabel("Soil Loss")
  108 +#plt.clabel("Fire Risk Protection")
  109 +
  110 +plt.show()
0 111 \ No newline at end of file
... ...