diff --git a/benchmarks/benchmarks.py b/benchmarks/benchmarks.py index 562750c..3d9bd9d 100644 --- a/benchmarks/benchmarks.py +++ b/benchmarks/benchmarks.py @@ -5,337 +5,345 @@ import matplotlib.colors as pltcolors import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from sklearn.cross_validation import KFold -import Measures +from pyFTS.benchmarks import Measures from pyFTS.partitioners import Grid -from pyFTS.common import Membership,FuzzySet,FLR,Transformations +from pyFTS.common import Membership, FuzzySet, FLR, Transformations + + +def getIntervalStatistics(original, models): + ret = "Model & RMSE & MAPE & Sharpness & Resolution & Coverage \\ \n" + for fts in models: + forecasts = fts.forecast(original) + ret = ret + fts.shortname + " & " + ret = ret + str(round(Measures.rmse_interval(original[fts.order - 1:], forecasts), 2)) + " & " + ret = ret + str(round(Measures.mape_interval(original[fts.order - 1:], forecasts), 2)) + " & " + ret = ret + str(round(Measures.sharpness(forecasts), 2)) + " & " + ret = ret + str(round(Measures.resolution(forecasts), 2)) + " & " + ret = ret + str(round(Measures.coverage(original[fts.order - 1:], forecasts), 2)) + " \\ \n" + return ret + -def getIntervalStatistics(original,models): - ret = "Model & RMSE & MAPE & Sharpness & Resolution & Coverage \\ \n" - for fts in models: - forecasts = fts.forecast(original) - ret = ret + fts.shortname + " & " - ret = ret + str( round(Measures.rmse_interval(original[fts.order-1 :],forecasts),2)) + " & " - ret = ret + str( round(Measures.mape_interval(original[fts.order-1 :],forecasts),2)) + " & " - ret = ret + str( round(Measures.sharpness(forecasts),2)) + " & " - ret = ret + str( round(Measures.resolution(forecasts),2)) + " & " - ret = ret + str( round(Measures.coverage(original[fts.order-1 :],forecasts),2)) + " \\ \n" - return ret - def plotDistribution(dist): - for k in dist.index: - alpha = np.array([dist[x][k] for x in dist])*100 - x = [k for x in np.arange(0,len(alpha))] - y = dist.columns - plt.scatter(x,y,c=alpha,marker='s',linewidths=0,cmap='Oranges',norm=pltcolors.Normalize(vmin=0,vmax=1),vmin=0,vmax=1,edgecolors=None) - -def plotComparedSeries(original,models, colors): - fig = plt.figure(figsize=[25,10]) - ax = fig.add_subplot(111) - - mi = [] - ma = [] - - ax.plot(original,color='black',label="Original") - count = 0 - for fts in models: - forecasted = fts.forecast(original) - - if fts.isInterval: - lower = [kk[0] for kk in forecasted] - upper = [kk[1] for kk in forecasted] - mi.append(min(lower)) - ma.append(max(upper)) - for k in np.arange(0,fts.order): - lower.insert(0,None) - upper.insert(0,None) - ax.plot(lower,color=colors[count],label=fts.shortname) - ax.plot(upper,color=colors[count]) - - else: - mi.append(min(forecasted)) - ma.append(max(forecasted)) - forecasted.insert(0,None) - ax.plot(forecasted,color=colors[count],label=fts.shortname) - - handles0, labels0 = ax.get_legend_handles_labels() - ax.legend(handles0,labels0) - count = count + 1 - #ax.set_title(fts.name) - ax.set_ylim([min(mi),max(ma)]) - ax.set_ylabel('F(T)') - ax.set_xlabel('T') - ax.set_xlim([0,len(original)]) + for k in dist.index: + alpha = np.array([dist[x][k] for x in dist]) * 100 + x = [k for x in np.arange(0, len(alpha))] + y = dist.columns + plt.scatter(x, y, c=alpha, marker='s', linewidths=0, cmap='Oranges', norm=pltcolors.Normalize(vmin=0, vmax=1), + vmin=0, vmax=1, edgecolors=None) -def plotComparedIntervalsAhead(original,models, colors, distributions, time_from, time_to): - fig = plt.figure(figsize=[25,10]) - ax = fig.add_subplot(111) - - mi = [] - ma = [] - - count = 0 - for fts in models: - if fts.isDensity and distributions[count]: - density = fts.forecastDistributionAhead(original[:time_from],time_to,25) - for k in density.index: - alpha = np.array([density[x][k] for x in density])*100 - x = [time_from + fts.order + k for x in np.arange(0,len(alpha))] - y = density.columns - ax.scatter(x,y,c=alpha,marker='s',linewidths=0,cmap='Oranges', - norm=pltcolors.Normalize(vmin=0,vmax=1),vmin=0,vmax=1,edgecolors=None) - - if fts.isInterval: - forecasts = fts.forecastAhead(original[:time_from],time_to) - lower = [kk[0] for kk in forecasts] - upper = [kk[1] for kk in forecasts] - mi.append(min(lower)) - ma.append(max(upper)) - for k in np.arange(0,time_from): - lower.insert(0,None) - upper.insert(0,None) - ax.plot(lower,color=colors[count],label=fts.shortname) - ax.plot(upper,color=colors[count]) - - else: - forecasts = fts.forecast(original) - mi.append(min(forecasts)) - ma.append(max(forecasts)) - for k in np.arange(0,time_from): - forecasts.insert(0,None) - ax.plot(forecasts,color=colors[count],label=fts.shortname) - - handles0, labels0 = ax.get_legend_handles_labels() - ax.legend(handles0,labels0) - count = count + 1 - ax.plot(original,color='black',label="Original") - #ax.set_title(fts.name) - ax.set_ylim([min(mi),max(ma)]) - ax.set_ylabel('F(T)') - ax.set_xlabel('T') - ax.set_xlim([0,len(original)]) +def plotComparedSeries(original, models, colors): + fig = plt.figure(figsize=[15, 5]) + ax = fig.add_subplot(111) + + mi = [] + ma = [] + + ax.plot(original, color='black', label="Original") + count = 0 + for fts in models: + if fts.hasPointForecasting: + forecasted = fts.forecast(original) + mi.append(min(forecasted)) + ma.append(max(forecasted)) + for k in np.arange(0, fts.order): + forecasted.insert(0, None) + ax.plot(forecasted, color=colors[count], label=fts.shortname, ls="-") + + if fts.hasIntervalForecasting: + forecasted = fts.forecastInterval(original) + lower = [kk[0] for kk in forecasted] + upper = [kk[1] for kk in forecasted] + mi.append(min(lower)) + ma.append(max(upper)) + for k in np.arange(0, fts.order): + lower.insert(0, None) + upper.insert(0, None) + ax.plot(lower, color=colors[count], label=fts.shortname,ls="--") + ax.plot(upper, color=colors[count],ls="--") + + handles0, labels0 = ax.get_legend_handles_labels() + ax.legend(handles0, labels0, loc=2) + count = count + 1 + # ax.set_title(fts.name) + ax.set_ylim([min(mi), max(ma)]) + ax.set_ylabel('F(T)') + ax.set_xlabel('T') + ax.set_xlim([0, len(original)]) -def plotCompared(original,forecasts,labels,title): - fig = plt.figure(figsize=[13,6]) - ax = fig.add_subplot(111) - ax.plot(original,color='k',label="Original") - for c in range(0,len(forecasts)): - ax.plot(forecasts[c],label=labels[c]) - handles0, labels0 = ax.get_legend_handles_labels() - ax.legend(handles0,labels0) - ax.set_title(title) - ax.set_ylabel('F(T)') - ax.set_xlabel('T') - ax.set_xlim([0,len(original)]) - ax.set_ylim([min(original),max(original)]) - -def SelecaoKFold_MenorRMSE(original,parameters,modelo): - nfolds = 5 - ret = [] - errors = np.array([[0 for k in parameters] for z in np.arange(0,nfolds)]) - forecasted_best = [] - print("Série Original") - fig = plt.figure(figsize=[18,10]) - fig.suptitle("Comparação de modelos ") - ax0 = fig.add_axes([0, 0.5, 0.65, 0.45]) #left, bottom, width, height - ax0.set_xlim([0,len(original)]) - ax0.set_ylim([min(original),max(original)]) - ax0.set_title('Série Temporal') - ax0.set_ylabel('F(T)') - ax0.set_xlabel('T') - ax0.plot(original,label="Original") - min_rmse_fold = 100000.0 - best = None - fc = 0 #Fold count - kf = KFold(len(original), n_folds=nfolds) - for train_ix, test_ix in kf: - train = original[train_ix] - test = original[test_ix] - min_rmse = 100000.0 - best_fold = None - forecasted_best_fold = [] - errors_fold = [] - pc = 0 #Parameter count - for p in parameters: - sets = Grid.GridPartitionerTrimf(train,p) - fts = modelo(str(p)+ " particoes") - fts.train(train,sets) - forecasted = [fts.forecast(xx) for xx in test] - error = Measures.rmse(np.array(forecasted),np.array(test)) - errors_fold.append(error) - print(fc, p, error) - errors[fc,pc] = error - if error < min_rmse: - min_rmse = error - best_fold = fts - forecasted_best_fold = forecasted - pc = pc + 1 - forecasted_best_fold = [best_fold.forecast(xx) for xx in original] - ax0.plot(forecasted_best_fold,label=best_fold.name) - if np.mean(errors_fold) < min_rmse_fold: - min_rmse_fold = np.mean(errors) - best = best_fold - forecasted_best = forecasted_best_fold - fc = fc + 1 - handles0, labels0 = ax0.get_legend_handles_labels() - ax0.legend(handles0, labels0) - ax1 = Axes3D(fig, rect=[0.7, 0.5, 0.3, 0.45], elev=30, azim=144) - #ax1 = fig.add_axes([0.6, 0.0, 0.45, 0.45], projection='3d') - ax1.set_title('Comparação dos Erros Quadráticos Médios') - ax1.set_zlabel('RMSE') - ax1.set_xlabel('K-fold') - ax1.set_ylabel('Partições') - X,Y = np.meshgrid(np.arange(0,nfolds),parameters) - surf = ax1.plot_surface(X, Y, errors.T, rstride=1, cstride=1, antialiased=True) - ret.append(best) - ret.append(forecasted_best) +def plotComparedIntervalsAhead(original, models, colors, distributions, time_from, time_to): + fig = plt.figure(figsize=[25, 10]) + ax = fig.add_subplot(111) - # Modelo diferencial - print("\nSérie Diferencial") - errors = np.array([[0 for k in parameters] for z in np.arange(0,nfolds)]) - forecastedd_best = [] - ax2 = fig.add_axes([0, 0, 0.65, 0.45]) #left, bottom, width, height - ax2.set_xlim([0,len(original)]) - ax2.set_ylim([min(original),max(original)]) - ax2.set_title('Série Temporal') - ax2.set_ylabel('F(T)') - ax2.set_xlabel('T') - ax2.plot(original,label="Original") - min_rmse = 100000.0 - min_rmse_fold = 100000.0 - bestd = None - fc = 0 - diff = Transformations.differential(original) - kf = KFold(len(original), n_folds=nfolds) - for train_ix, test_ix in kf: - train = diff[train_ix] - test = diff[test_ix] - min_rmse = 100000.0 - best_fold = None - forecasted_best_fold = [] - errors_fold = [] - pc = 0 - for p in parameters: - sets = Grid.GridPartitionerTrimf(train,p) - fts = modelo(str(p)+ " particoes") - fts.train(train,sets) - forecasted = [fts.forecastDiff(test,xx) for xx in np.arange(len(test))] - error = Measures.rmse(np.array(forecasted),np.array(test)) - print(fc, p,error) - errors[fc,pc] = error - errors_fold.append(error) - if error < min_rmse: - min_rmse = error - best_fold = fts - pc = pc + 1 - forecasted_best_fold = [best_fold.forecastDiff(original, xx) for xx in np.arange(len(original))] - ax2.plot(forecasted_best_fold,label=best_fold.name) - if np.mean(errors_fold) < min_rmse_fold: - min_rmse_fold = np.mean(errors) - best = best_fold - forecasted_best = forecasted_best_fold - fc = fc + 1 - handles0, labels0 = ax2.get_legend_handles_labels() - ax2.legend(handles0, labels0) - ax3 = Axes3D(fig, rect=[0.7, 0, 0.3, 0.45], elev=30, azim=144) - #ax1 = fig.add_axes([0.6, 0.0, 0.45, 0.45], projection='3d') - ax3.set_title('Comparação dos Erros Quadráticos Médios') - ax3.set_zlabel('RMSE') - ax3.set_xlabel('K-fold') - ax3.set_ylabel('Partições') - X,Y = np.meshgrid(np.arange(0,nfolds),parameters) - surf = ax3.plot_surface(X, Y, errors.T, rstride=1, cstride=1, antialiased=True) - ret.append(best) - ret.append(forecasted_best) - return ret - -def SelecaoSimples_MenorRMSE(original,parameters,modelo): - ret = [] - errors = [] - forecasted_best = [] - print("Série Original") - fig = plt.figure(figsize=[20,12]) - fig.suptitle("Comparação de modelos ") - ax0 = fig.add_axes([0, 0.5, 0.65, 0.45]) #left, bottom, width, height - ax0.set_xlim([0,len(original)]) - ax0.set_ylim([min(original),max(original)]) - ax0.set_title('Série Temporal') - ax0.set_ylabel('F(T)') - ax0.set_xlabel('T') - ax0.plot(original,label="Original") - min_rmse = 100000.0 - best = None - for p in parameters: - sets = Grid.GridPartitionerTrimf(original,p) - fts = modelo(str(p)+ " particoes") - fts.train(original,sets) - #print(original) - forecasted = fts.forecast(original) - forecasted.insert(0,original[0]) - #print(forecasted) - ax0.plot(forecasted,label=fts.name) - error = Measures.rmse(np.array(forecasted),np.array(original)) - print(p,error) - errors.append(error) - if error < min_rmse: - min_rmse = error - best = fts - forecasted_best = forecasted - handles0, labels0 = ax0.get_legend_handles_labels() - ax0.legend(handles0, labels0) - ax1 = fig.add_axes([0.7, 0.5, 0.3, 0.45]) #left, bottom, width, height - ax1.set_title('Comparação dos Erros Quadráticos Médios') - ax1.set_ylabel('RMSE') - ax1.set_xlabel('Quantidade de Partições') - ax1.set_xlim([min(parameters),max(parameters)]) - ax1.plot(parameters,errors) - ret.append(best) - ret.append(forecasted_best) - # Modelo diferencial - print("\nSérie Diferencial") - difffts = Transformations.differential(original) - errors = [] - forecastedd_best = [] - ax2 = fig.add_axes([0, 0, 0.65, 0.45]) #left, bottom, width, height - ax2.set_xlim([0,len(difffts)]) - ax2.set_ylim([min(difffts),max(difffts)]) - ax2.set_title('Série Temporal') - ax2.set_ylabel('F(T)') - ax2.set_xlabel('T') - ax2.plot(difffts,label="Original") - min_rmse = 100000.0 - bestd = None - for p in parameters: - sets = Grid.GridPartitionerTrimf(difffts,p) - fts = modelo(str(p)+ " particoes") - fts.train(difffts,sets) - forecasted = fts.forecast(difffts) - forecasted.insert(0,difffts[0]) - ax2.plot(forecasted,label=fts.name) - error = Measures.rmse(np.array(forecasted),np.array(difffts)) - print(p,error) - errors.append(error) - if error < min_rmse: - min_rmse = error - bestd = fts - forecastedd_best = forecasted - handles0, labels0 = ax2.get_legend_handles_labels() - ax2.legend(handles0, labels0) - ax3 = fig.add_axes([0.7, 0, 0.3, 0.45]) #left, bottom, width, height - ax3.set_title('Comparação dos Erros Quadráticos Médios') - ax3.set_ylabel('RMSE') - ax3.set_xlabel('Quantidade de Partições') - ax3.set_xlim([min(parameters),max(parameters)]) - ax3.plot(parameters,errors) - ret.append(bestd) - ret.append(forecastedd_best) - return ret - -def compareModelsPlot(original,models_fo,models_ho): - fig = plt.figure(figsize=[13,6]) + mi = [] + ma = [] + + count = 0 + for fts in models: + if fts.hasDistributionForecasting and distributions[count]: + density = fts.forecastDistributionAhead(original[:time_from], time_to, 25) + for k in density.index: + alpha = np.array([density[x][k] for x in density]) * 100 + x = [time_from + fts.order + k for x in np.arange(0, len(alpha))] + y = density.columns + ax.scatter(x, y, c=alpha, marker='s', linewidths=0, cmap='Oranges', + norm=pltcolors.Normalize(vmin=0, vmax=1), vmin=0, vmax=1, edgecolors=None) + + if fts.hasIntervalForecasting: + forecasts = fts.forecastAhead(original[:time_from], time_to) + lower = [kk[0] for kk in forecasts] + upper = [kk[1] for kk in forecasts] + mi.append(min(lower)) + ma.append(max(upper)) + for k in np.arange(0, time_from): + lower.insert(0, None) + upper.insert(0, None) + ax.plot(lower, color=colors[count], label=fts.shortname) + ax.plot(upper, color=colors[count]) + + else: + forecasts = fts.forecast(original) + mi.append(min(forecasts)) + ma.append(max(forecasts)) + for k in np.arange(0, time_from): + forecasts.insert(0, None) + ax.plot(forecasts, color=colors[count], label=fts.shortname) + + handles0, labels0 = ax.get_legend_handles_labels() + ax.legend(handles0, labels0) + count = count + 1 + ax.plot(original, color='black', label="Original") + # ax.set_title(fts.name) + ax.set_ylim([min(mi), max(ma)]) + ax.set_ylabel('F(T)') + ax.set_xlabel('T') + ax.set_xlim([0, len(original)]) + + +def plotCompared(original, forecasts, labels, title): + fig = plt.figure(figsize=[13, 6]) + ax = fig.add_subplot(111) + ax.plot(original, color='k', label="Original") + for c in range(0, len(forecasts)): + ax.plot(forecasts[c], label=labels[c]) + handles0, labels0 = ax.get_legend_handles_labels() + ax.legend(handles0, labels0) + ax.set_title(title) + ax.set_ylabel('F(T)') + ax.set_xlabel('T') + ax.set_xlim([0, len(original)]) + ax.set_ylim([min(original), max(original)]) + + +def SelecaoKFold_MenorRMSE(original, parameters, modelo): + nfolds = 5 + ret = [] + errors = np.array([[0 for k in parameters] for z in np.arange(0, nfolds)]) + forecasted_best = [] + print("Série Original") + fig = plt.figure(figsize=[18, 10]) fig.suptitle("Comparação de modelos ") - ax0 = fig.add_axes([0, 0, 1, 1]) #left, bottom, width, height + ax0 = fig.add_axes([0, 0.5, 0.65, 0.45]) # left, bottom, width, height + ax0.set_xlim([0, len(original)]) + ax0.set_ylim([min(original), max(original)]) + ax0.set_title('Série Temporal') + ax0.set_ylabel('F(T)') + ax0.set_xlabel('T') + ax0.plot(original, label="Original") + min_rmse_fold = 100000.0 + best = None + fc = 0 # Fold count + kf = KFold(len(original), n_folds=nfolds) + for train_ix, test_ix in kf: + train = original[train_ix] + test = original[test_ix] + min_rmse = 100000.0 + best_fold = None + forecasted_best_fold = [] + errors_fold = [] + pc = 0 # Parameter count + for p in parameters: + sets = Grid.GridPartitionerTrimf(train, p) + fts = modelo(str(p) + " particoes") + fts.train(train, sets) + forecasted = [fts.forecast(xx) for xx in test] + error = Measures.rmse(np.array(forecasted), np.array(test)) + errors_fold.append(error) + print(fc, p, error) + errors[fc, pc] = error + if error < min_rmse: + min_rmse = error + best_fold = fts + forecasted_best_fold = forecasted + pc = pc + 1 + forecasted_best_fold = [best_fold.forecast(xx) for xx in original] + ax0.plot(forecasted_best_fold, label=best_fold.name) + if np.mean(errors_fold) < min_rmse_fold: + min_rmse_fold = np.mean(errors) + best = best_fold + forecasted_best = forecasted_best_fold + fc = fc + 1 + handles0, labels0 = ax0.get_legend_handles_labels() + ax0.legend(handles0, labels0) + ax1 = Axes3D(fig, rect=[0.7, 0.5, 0.3, 0.45], elev=30, azim=144) + # ax1 = fig.add_axes([0.6, 0.0, 0.45, 0.45], projection='3d') + ax1.set_title('Comparação dos Erros Quadráticos Médios') + ax1.set_zlabel('RMSE') + ax1.set_xlabel('K-fold') + ax1.set_ylabel('Partições') + X, Y = np.meshgrid(np.arange(0, nfolds), parameters) + surf = ax1.plot_surface(X, Y, errors.T, rstride=1, cstride=1, antialiased=True) + ret.append(best) + ret.append(forecasted_best) + + # Modelo diferencial + print("\nSérie Diferencial") + errors = np.array([[0 for k in parameters] for z in np.arange(0, nfolds)]) + forecastedd_best = [] + ax2 = fig.add_axes([0, 0, 0.65, 0.45]) # left, bottom, width, height + ax2.set_xlim([0, len(original)]) + ax2.set_ylim([min(original), max(original)]) + ax2.set_title('Série Temporal') + ax2.set_ylabel('F(T)') + ax2.set_xlabel('T') + ax2.plot(original, label="Original") + min_rmse = 100000.0 + min_rmse_fold = 100000.0 + bestd = None + fc = 0 + diff = Transformations.differential(original) + kf = KFold(len(original), n_folds=nfolds) + for train_ix, test_ix in kf: + train = diff[train_ix] + test = diff[test_ix] + min_rmse = 100000.0 + best_fold = None + forecasted_best_fold = [] + errors_fold = [] + pc = 0 + for p in parameters: + sets = Grid.GridPartitionerTrimf(train, p) + fts = modelo(str(p) + " particoes") + fts.train(train, sets) + forecasted = [fts.forecastDiff(test, xx) for xx in np.arange(len(test))] + error = Measures.rmse(np.array(forecasted), np.array(test)) + print(fc, p, error) + errors[fc, pc] = error + errors_fold.append(error) + if error < min_rmse: + min_rmse = error + best_fold = fts + pc = pc + 1 + forecasted_best_fold = [best_fold.forecastDiff(original, xx) for xx in np.arange(len(original))] + ax2.plot(forecasted_best_fold, label=best_fold.name) + if np.mean(errors_fold) < min_rmse_fold: + min_rmse_fold = np.mean(errors) + best = best_fold + forecasted_best = forecasted_best_fold + fc = fc + 1 + handles0, labels0 = ax2.get_legend_handles_labels() + ax2.legend(handles0, labels0) + ax3 = Axes3D(fig, rect=[0.7, 0, 0.3, 0.45], elev=30, azim=144) + # ax1 = fig.add_axes([0.6, 0.0, 0.45, 0.45], projection='3d') + ax3.set_title('Comparação dos Erros Quadráticos Médios') + ax3.set_zlabel('RMSE') + ax3.set_xlabel('K-fold') + ax3.set_ylabel('Partições') + X, Y = np.meshgrid(np.arange(0, nfolds), parameters) + surf = ax3.plot_surface(X, Y, errors.T, rstride=1, cstride=1, antialiased=True) + ret.append(best) + ret.append(forecasted_best) + return ret + + +def SelecaoSimples_MenorRMSE(original, parameters, modelo): + ret = [] + errors = [] + forecasted_best = [] + print("Série Original") + fig = plt.figure(figsize=[20, 12]) + fig.suptitle("Comparação de modelos ") + ax0 = fig.add_axes([0, 0.5, 0.65, 0.45]) # left, bottom, width, height + ax0.set_xlim([0, len(original)]) + ax0.set_ylim([min(original), max(original)]) + ax0.set_title('Série Temporal') + ax0.set_ylabel('F(T)') + ax0.set_xlabel('T') + ax0.plot(original, label="Original") + min_rmse = 100000.0 + best = None + for p in parameters: + sets = Grid.GridPartitionerTrimf(original, p) + fts = modelo(str(p) + " particoes") + fts.train(original, sets) + # print(original) + forecasted = fts.forecast(original) + forecasted.insert(0, original[0]) + # print(forecasted) + ax0.plot(forecasted, label=fts.name) + error = Measures.rmse(np.array(forecasted), np.array(original)) + print(p, error) + errors.append(error) + if error < min_rmse: + min_rmse = error + best = fts + forecasted_best = forecasted + handles0, labels0 = ax0.get_legend_handles_labels() + ax0.legend(handles0, labels0) + ax1 = fig.add_axes([0.7, 0.5, 0.3, 0.45]) # left, bottom, width, height + ax1.set_title('Comparação dos Erros Quadráticos Médios') + ax1.set_ylabel('RMSE') + ax1.set_xlabel('Quantidade de Partições') + ax1.set_xlim([min(parameters), max(parameters)]) + ax1.plot(parameters, errors) + ret.append(best) + ret.append(forecasted_best) + # Modelo diferencial + print("\nSérie Diferencial") + difffts = Transformations.differential(original) + errors = [] + forecastedd_best = [] + ax2 = fig.add_axes([0, 0, 0.65, 0.45]) # left, bottom, width, height + ax2.set_xlim([0, len(difffts)]) + ax2.set_ylim([min(difffts), max(difffts)]) + ax2.set_title('Série Temporal') + ax2.set_ylabel('F(T)') + ax2.set_xlabel('T') + ax2.plot(difffts, label="Original") + min_rmse = 100000.0 + bestd = None + for p in parameters: + sets = Grid.GridPartitionerTrimf(difffts, p) + fts = modelo(str(p) + " particoes") + fts.train(difffts, sets) + forecasted = fts.forecast(difffts) + forecasted.insert(0, difffts[0]) + ax2.plot(forecasted, label=fts.name) + error = Measures.rmse(np.array(forecasted), np.array(difffts)) + print(p, error) + errors.append(error) + if error < min_rmse: + min_rmse = error + bestd = fts + forecastedd_best = forecasted + handles0, labels0 = ax2.get_legend_handles_labels() + ax2.legend(handles0, labels0) + ax3 = fig.add_axes([0.7, 0, 0.3, 0.45]) # left, bottom, width, height + ax3.set_title('Comparação dos Erros Quadráticos Médios') + ax3.set_ylabel('RMSE') + ax3.set_xlabel('Quantidade de Partições') + ax3.set_xlim([min(parameters), max(parameters)]) + ax3.plot(parameters, errors) + ret.append(bestd) + ret.append(forecastedd_best) + return ret + + +def compareModelsPlot(original, models_fo, models_ho): + fig = plt.figure(figsize=[13, 6]) + fig.suptitle("Comparação de modelos ") + ax0 = fig.add_axes([0, 0, 1, 1]) # left, bottom, width, height rows = [] for model in models_fo: fts = model["model"] @@ -345,29 +353,30 @@ def compareModelsPlot(original,models_fo,models_ho): ax0.plot(model["forecasted"], label=model["name"]) handles0, labels0 = ax0.get_legend_handles_labels() ax0.legend(handles0, labels0) - -def compareModelsTable(original,models_fo,models_ho): - fig = plt.figure(figsize=[12,4]) + + +def compareModelsTable(original, models_fo, models_ho): + fig = plt.figure(figsize=[12, 4]) fig.suptitle("Comparação de modelos ") - columns = ['Modelo','Ordem','Partições','RMSE','MAPE (%)'] + columns = ['Modelo', 'Ordem', 'Partições', 'RMSE', 'MAPE (%)'] rows = [] for model in models_fo: fts = model["model"] - error_r = Measures.rmse(model["forecasted"],original) - error_m = round(Measures.mape(model["forecasted"],original)*100,2) - rows.append([model["name"],fts.order,len(fts.sets),error_r,error_m]) + error_r = Measures.rmse(model["forecasted"], original) + error_m = round(Measures.mape(model["forecasted"], original) * 100, 2) + rows.append([model["name"], fts.order, len(fts.sets), error_r, error_m]) for model in models_ho: fts = model["model"] - error_r = Measures.rmse(model["forecasted"][fts.order:],original[fts.order:]) - error_m = round(Measures.mape(model["forecasted"][fts.order:],original[fts.order:])*100,2) - rows.append([model["name"],fts.order,len(fts.sets),error_r,error_m]) - ax1 = fig.add_axes([0, 0, 1, 1]) #left, bottom, width, height + error_r = Measures.rmse(model["forecasted"][fts.order:], original[fts.order:]) + error_m = round(Measures.mape(model["forecasted"][fts.order:], original[fts.order:]) * 100, 2) + rows.append([model["name"], fts.order, len(fts.sets), error_r, error_m]) + ax1 = fig.add_axes([0, 0, 1, 1]) # left, bottom, width, height ax1.set_xticks([]) ax1.set_yticks([]) ax1.table(cellText=rows, - colLabels=columns, - cellLoc='center', - bbox=[0,0,1,1]) + colLabels=columns, + cellLoc='center', + bbox=[0, 0, 1, 1]) sup = "\\begin{tabular}{" header = "" body = "" @@ -379,113 +388,115 @@ def compareModelsTable(original,models_fo,models_ho): header = header + " & " header = header + "\\textbf{" + c + "} " sup = sup + "|} \\hline\n" - header = header + "\\\\ \\hline \n" - + header = header + "\\\\ \\hline \n" + for r in rows: lin = "" for c in r: if len(lin) > 0: lin = lin + " & " lin = lin + str(c) - - body = body + lin + "\\\\ \\hline \n" - + + body = body + lin + "\\\\ \\hline \n" + return sup + header + body + "\\end{tabular}" + from pyFTS import hwang -def HOSelecaoSimples_MenorRMSE(original,parameters,orders): - ret = [] - errors = np.array([[0 for k in range(len(parameters))] for kk in range(len(orders))]) - forecasted_best = [] - print("Série Original") - fig = plt.figure(figsize=[20,12]) - fig.suptitle("Comparação de modelos ") - ax0 = fig.add_axes([0, 0.5, 0.6, 0.45]) #left, bottom, width, height - ax0.set_xlim([0,len(original)]) - ax0.set_ylim([min(original),max(original)]) - ax0.set_title('Série Temporal') - ax0.set_ylabel('F(T)') - ax0.set_xlabel('T') - ax0.plot(original,label="Original") - min_rmse = 100000.0 - best = None - pc = 0 - for p in parameters: - oc = 0 - for o in orders: - sets = Grid.GridPartitionerTrimf(original,p) - fts = hwang.HighOrderFTS(o,"k = " + str(p)+ " w = " + str(o)) - fts.train(original,sets) - forecasted = [fts.forecast(original, xx) for xx in range(o,len(original))] - error = Measures.rmse(np.array(forecasted),np.array(original[o:])) - for kk in range(o): - forecasted.insert(0,None) - ax0.plot(forecasted,label=fts.name) - print(o,p,error) - errors[oc,pc] = error - if error < min_rmse: - min_rmse = error - best = fts - forecasted_best = forecasted - oc = oc + 1 - pc = pc + 1 - handles0, labels0 = ax0.get_legend_handles_labels() - ax0.legend(handles0, labels0) - ax1 = Axes3D(fig, rect=[0.6, 0.5, 0.45, 0.45], elev=30, azim=144) - #ax1 = fig.add_axes([0.6, 0.5, 0.45, 0.45], projection='3d') - ax1.set_title('Comparação dos Erros Quadráticos Médios por tamanho da janela') - ax1.set_ylabel('RMSE') - ax1.set_xlabel('Quantidade de Partições') - ax1.set_zlabel('W') - X,Y = np.meshgrid(parameters,orders) - surf = ax1.plot_surface(X, Y, errors, rstride=1, cstride=1, antialiased=True) - ret.append(best) - ret.append(forecasted_best) + +def HOSelecaoSimples_MenorRMSE(original, parameters, orders): + ret = [] + errors = np.array([[0 for k in range(len(parameters))] for kk in range(len(orders))]) + forecasted_best = [] + print("Série Original") + fig = plt.figure(figsize=[20, 12]) + fig.suptitle("Comparação de modelos ") + ax0 = fig.add_axes([0, 0.5, 0.6, 0.45]) # left, bottom, width, height + ax0.set_xlim([0, len(original)]) + ax0.set_ylim([min(original), max(original)]) + ax0.set_title('Série Temporal') + ax0.set_ylabel('F(T)') + ax0.set_xlabel('T') + ax0.plot(original, label="Original") + min_rmse = 100000.0 + best = None + pc = 0 + for p in parameters: + oc = 0 + for o in orders: + sets = Grid.GridPartitionerTrimf(original, p) + fts = hwang.HighOrderFTS(o, "k = " + str(p) + " w = " + str(o)) + fts.train(original, sets) + forecasted = [fts.forecast(original, xx) for xx in range(o, len(original))] + error = Measures.rmse(np.array(forecasted), np.array(original[o:])) + for kk in range(o): + forecasted.insert(0, None) + ax0.plot(forecasted, label=fts.name) + print(o, p, error) + errors[oc, pc] = error + if error < min_rmse: + min_rmse = error + best = fts + forecasted_best = forecasted + oc = oc + 1 + pc = pc + 1 + handles0, labels0 = ax0.get_legend_handles_labels() + ax0.legend(handles0, labels0) + ax1 = Axes3D(fig, rect=[0.6, 0.5, 0.45, 0.45], elev=30, azim=144) + # ax1 = fig.add_axes([0.6, 0.5, 0.45, 0.45], projection='3d') + ax1.set_title('Comparação dos Erros Quadráticos Médios por tamanho da janela') + ax1.set_ylabel('RMSE') + ax1.set_xlabel('Quantidade de Partições') + ax1.set_zlabel('W') + X, Y = np.meshgrid(parameters, orders) + surf = ax1.plot_surface(X, Y, errors, rstride=1, cstride=1, antialiased=True) + ret.append(best) + ret.append(forecasted_best) # Modelo diferencial - print("\nSérie Diferencial") - errors = np.array([[0 for k in range(len(parameters))] for kk in range(len(orders))]) - forecastedd_best = [] - ax2 = fig.add_axes([0, 0, 0.6, 0.45]) #left, bottom, width, height - ax2.set_xlim([0,len(original)]) - ax2.set_ylim([min(original),max(original)]) - ax2.set_title('Série Temporal') - ax2.set_ylabel('F(T)') - ax2.set_xlabel('T') - ax2.plot(original,label="Original") - min_rmse = 100000.0 - bestd = None - pc = 0 - for p in parameters: - oc = 0 - for o in orders: - sets = Grid.GridPartitionerTrimf(Transformations.differential(original),p) - fts = hwang.HighOrderFTS(o,"k = " + str(p)+ " w = " + str(o)) - fts.train(original,sets) - forecasted = [fts.forecastDiff(original, xx) for xx in range(o,len(original))] - error = Measures.rmse(np.array(forecasted),np.array(original[o:])) - for kk in range(o): - forecasted.insert(0,None) - ax2.plot(forecasted,label=fts.name) - print(o,p,error) - errors[oc,pc] = error - if error < min_rmse: - min_rmse = error - bestd = fts - forecastedd_best = forecasted - oc = oc + 1 - pc = pc + 1 - handles0, labels0 = ax2.get_legend_handles_labels() - ax2.legend(handles0, labels0) - ax3 = Axes3D(fig, rect=[0.6, 0.0, 0.45, 0.45], elev=30, azim=144) - #ax3 = fig.add_axes([0.6, 0.0, 0.45, 0.45], projection='3d') - ax3.set_title('Comparação dos Erros Quadráticos Médios') - ax3.set_ylabel('RMSE') - ax3.set_xlabel('Quantidade de Partições') - ax3.set_zlabel('W') - X,Y = np.meshgrid(parameters,orders) - surf = ax3.plot_surface(X, Y, errors, rstride=1, cstride=1, antialiased=True) - ret.append(bestd) - ret.append(forecastedd_best) - return ret + print("\nSérie Diferencial") + errors = np.array([[0 for k in range(len(parameters))] for kk in range(len(orders))]) + forecastedd_best = [] + ax2 = fig.add_axes([0, 0, 0.6, 0.45]) # left, bottom, width, height + ax2.set_xlim([0, len(original)]) + ax2.set_ylim([min(original), max(original)]) + ax2.set_title('Série Temporal') + ax2.set_ylabel('F(T)') + ax2.set_xlabel('T') + ax2.plot(original, label="Original") + min_rmse = 100000.0 + bestd = None + pc = 0 + for p in parameters: + oc = 0 + for o in orders: + sets = Grid.GridPartitionerTrimf(Transformations.differential(original), p) + fts = hwang.HighOrderFTS(o, "k = " + str(p) + " w = " + str(o)) + fts.train(original, sets) + forecasted = [fts.forecastDiff(original, xx) for xx in range(o, len(original))] + error = Measures.rmse(np.array(forecasted), np.array(original[o:])) + for kk in range(o): + forecasted.insert(0, None) + ax2.plot(forecasted, label=fts.name) + print(o, p, error) + errors[oc, pc] = error + if error < min_rmse: + min_rmse = error + bestd = fts + forecastedd_best = forecasted + oc = oc + 1 + pc = pc + 1 + handles0, labels0 = ax2.get_legend_handles_labels() + ax2.legend(handles0, labels0) + ax3 = Axes3D(fig, rect=[0.6, 0.0, 0.45, 0.45], elev=30, azim=144) + # ax3 = fig.add_axes([0.6, 0.0, 0.45, 0.45], projection='3d') + ax3.set_title('Comparação dos Erros Quadráticos Médios') + ax3.set_ylabel('RMSE') + ax3.set_xlabel('Quantidade de Partições') + ax3.set_zlabel('W') + X, Y = np.meshgrid(parameters, orders) + surf = ax3.plot_surface(X, Y, errors, rstride=1, cstride=1, antialiased=True) + ret.append(bestd) + ret.append(forecastedd_best) + return ret diff --git a/common/FLR.py b/common/FLR.py index caacb86..3bd9338 100644 --- a/common/FLR.py +++ b/common/FLR.py @@ -1,10 +1,14 @@ +import numpy as np + + class FLR: def __init__(self, LHS, RHS): self.LHS = LHS self.RHS = RHS def __str__(self): - return str(self.LHS) + " -> " + str(self.RHS) + return self.LHS.name + " -> " + self.RHS.name + def generateNonRecurrentFLRs(fuzzyData): flrs = {} diff --git a/common/FuzzySet.py b/common/FuzzySet.py index 908eebd..e4d6843 100644 --- a/common/FuzzySet.py +++ b/common/FuzzySet.py @@ -20,7 +20,7 @@ class FuzzySet: return self.mf(x, self.parameters) def __str__(self): - return self.name + ": " + str(self.mf) + "(" + str(self.parameters) + ")" + return self.name + ": " + str(self.mf.__name__) + "(" + str(self.parameters) + ")" def fuzzyInstance(inst, fuzzySets): diff --git a/fts.py b/fts.py index 877add9..1ede686 100644 --- a/fts.py +++ b/fts.py @@ -10,9 +10,11 @@ class FTS: self.shortname = name self.name = name self.detail = name - self.isSeasonal = False - self.isInterval = False - self.isDensity = False + self.hasSeasonality = False + self.hasPointForecasting = True + self.hasIntervalForecasting = False + self.hasDistributionForecasting = False + self.dump = False def fuzzy(self, data): best = {"fuzzyset": "", "membership": 0.0} @@ -28,6 +30,21 @@ class FTS: def forecast(self, data): pass + def forecastInterval(self, data): + pass + + def forecastDistribution(self, data): + pass + + def forecastAhead(self, data, steps): + pass + + def forecastAheadInterval(self, data, steps): + pass + + def forecastAheadDistribution(self, data, steps): + pass + def train(self, data, sets): pass diff --git a/hofts.py b/hofts.py index b3db5d5..ec0bd6c 100644 --- a/hofts.py +++ b/hofts.py @@ -1,6 +1,6 @@ import numpy as np from pyFTS.common import FuzzySet,FLR -import fts +from pyFTS import fts class HighOrderFLRG: @@ -18,7 +18,7 @@ class HighOrderFLRG: if len(self.strlhs) == 0: for c in self.LHS: if len(self.strlhs) > 0: - self.strlhs = self.strlhs + ", " + self.strlhs += ", " self.strlhs = self.strlhs + c.name return self.strlhs @@ -63,7 +63,7 @@ class HighOrderFTS(fts.FTS): self.sets = sets for s in self.sets: self.setsDict[s.name] = s tmpdata = FuzzySet.fuzzySeries(data, sets) - flrs = FuzzySet.generateRecurrentFLRs(tmpdata) + flrs = FLR.generateRecurrentFLRs(tmpdata) self.flrgs = self.generateFLRG(flrs) def getMidpoints(self, flrg): diff --git a/hwang.py b/hwang.py index 81e660b..736b16a 100644 --- a/hwang.py +++ b/hwang.py @@ -1,6 +1,6 @@ import numpy as np from pyFTS.common import FuzzySet,FLR,Transformations -import fts +from pyFTS import fts class HighOrderFTS(fts.FTS): diff --git a/ifts.py b/ifts.py index 0bd7dfe..b6f98a4 100644 --- a/ifts.py +++ b/ifts.py @@ -1,6 +1,6 @@ import numpy as np from pyFTS.common import FuzzySet,FLR -import hofts, fts, tree +from pyFTS import hofts, fts, tree class IntervalFTS(hofts.HighOrderFTS): @@ -10,7 +10,8 @@ class IntervalFTS(hofts.HighOrderFTS): self.name = "Interval FTS" self.detail = "Silva, P.; Guimarães, F.; Sadaei, H. (2016)" self.flrgs = {} - self.isInterval = True + self.hasPointForecasting = False + self.hasIntervalForecasting = True def getUpper(self, flrg): if flrg.strLHS() in self.flrgs: diff --git a/partitioners/Grid.py b/partitioners/Grid.py index a34624d..91bd9af 100644 --- a/partitioners/Grid.py +++ b/partitioners/Grid.py @@ -10,17 +10,21 @@ from pyFTS.common import FuzzySet, Membership def GridPartitionerTrimf(data, npart, names=None, prefix="A"): sets = [] dmax = max(data) - dmax += dmax * 0.10 + dmax += dmax * 0.1 + print(dmax) dmin = min(data) - dmin -= dmin * 0.10 + dmin -= dmin * 0.1 + print(dmin) dlen = dmax - dmin partlen = math.ceil(dlen / npart) - partition = math.ceil(dmin) - for c in range(npart): + #p2 = partlen / 2 + #partition = dmin #+ partlen + count = 0 + for c in np.arange(dmin, dmax, partlen): sets.append( - FuzzySet.FuzzySet(prefix + str(c), Membership.trimf, [round(partition - partlen, 3), partition, partition + partlen], - partition)) - partition += partlen + FuzzySet.FuzzySet(prefix + str(count), Membership.trimf, [c - partlen, c, c + partlen],c)) + count += 1 + #partition += partlen return sets diff --git a/pfts.py b/pfts.py index 2ddce65..86deda2 100644 --- a/pfts.py +++ b/pfts.py @@ -2,21 +2,21 @@ import numpy as np import pandas as pd import math from pyFTS.common import FuzzySet, FLR -import hofts, ifts, tree +from pyFTS import hofts, ifts, tree class ProbabilisticFLRG(hofts.HighOrderFLRG): def __init__(self, order): super(ProbabilisticFLRG, self).__init__(order) self.RHS = {} - self.frequencyCount = 0 + self.frequencyCount = 0.0 def appendRHS(self, c): - self.frequencyCount = self.frequencyCount + 1 + self.frequencyCount += 1 if c.name in self.RHS: - self.RHS[c.name] = self.RHS[c.name] + 1 + self.RHS[c.name] += 1 else: - self.RHS[c.name] = 1 + self.RHS[c.name] = 1.0 def getProbability(self, c): return self.RHS[c] / self.frequencyCount @@ -38,23 +38,27 @@ class ProbabilisticFTS(ifts.IntervalFTS): self.detail = "Silva, P.; Guimarães, F.; Sadaei, H." self.flrgs = {} self.globalFrequency = 0 - self.isInterval = True - self.isDensity = True + self.hasPointForecasting = True + self.hasIntervalForecasting = True + self.hasDistributionForecasting = True def generateFLRG(self, flrs): flrgs = {} l = len(flrs) - for k in np.arange(self.order + 1, l): + for k in np.arange(self.order, l+1): + if self.dump: print("FLR: " + str(k)) flrg = ProbabilisticFLRG(self.order) for kk in np.arange(k - self.order, k): flrg.appendLHS(flrs[kk].LHS) + if self.dump: print("LHS: " + str(flrs[kk])) if flrg.strLHS() in flrgs: - flrgs[flrg.strLHS()].appendRHS(flrs[k].RHS) + flrgs[flrg.strLHS()].appendRHS(flrs[k-1].RHS) else: flrgs[flrg.strLHS()] = flrg; - flrgs[flrg.strLHS()].appendRHS(flrs[k].RHS) + flrgs[flrg.strLHS()].appendRHS(flrs[k-1].RHS) + if self.dump: print("RHS: " + str(flrs[k-1])) self.globalFrequency = self.globalFrequency + 1 return (flrgs) @@ -68,9 +72,9 @@ class ProbabilisticFTS(ifts.IntervalFTS): def getMidpoints(self, flrg): if flrg.strLHS() in self.flrgs: tmp = self.flrgs[flrg.strLHS()] - ret = sum(np.array([tmp.getProbability(s) * self.setsDict[s].midpoint for s in tmp.RHS])) + ret = sum(np.array([tmp.getProbability(s) * self.setsDict[s].centroid for s in tmp.RHS])) else: - ret = sum(np.array([0.33 * s.midpoint for s in flrg.LHS])) + ret = sum(np.array([0.33 * s.centroid for s in flrg.LHS])) return ret def getUpper(self, flrg): diff --git a/sfts.py b/sfts.py index 2d38ddd..6565143 100644 --- a/sfts.py +++ b/sfts.py @@ -27,7 +27,7 @@ class SeasonalFTS(fts.FTS): self.name = "Seasonal FTS" self.detail = "Chen" self.seasonality = 1 - self.isSeasonal = True + self.hasSeasonality = True def generateFLRG(self, flrs): flrgs = []