%This program is useful to forecaste the ridership of Beijing 
%public transportation


% Get the input sampledata from the .txt file "sample.txt"
fid=fopen('sample.txt','rt');
originalData=fscanf(fid,'%f', [20,14]);
status=fclose(fid);

% transport sampledata into input training sampledata
inputSampledata=originalData';

% Get the output sampledata from the .txt file "goal.txt"
fid=fopen('goal.txt','rt');
outputData=fscanf(fid,'%f');
status=fclose(fid);
% transport sampledata into input training sampledata
outputSampledata=outputData';

% creating neural network and setting trainging parameters
gwwnet=newff(minmax(inputSampledata),[4,1],{'tansig','purelin'},'traingdm');
gwwnet.trainParam.show = 50;
gwwnet.trainParam.lr = 0.05;
gwwnet.trainParam.epochs = 8000;
gwwnet.trainParam.goal = 1e-3;

%data scaling (converting the network input and output data to the intervel [-1,1])
[input,mininput,maxinput,output,minoutput,maxoutput] = premnmx(inputSampledata,outputSampledata);

%training
[gwwnet,tr]=train(gwwnet,input,output);
y=sim(gwwnet,input);

%data offset (converting the network output data to it original unit)
nnoutput = postmnmx(y,minoutput,maxoutput);

%plot
time=1978:1:1997;
plot(time,outputSampledata,'-',time,nnoutput,'o');
legend('actual output','NN output');
xlabel('time');ylabel('Learning fitting curve');


%scenario1 forecasting process
column=10;
for i=1:column;    
      SceInput(1,i)=inputSampledata(1,20)*(1.0464^i);
      SceInput(2,i)=inputSampledata(2,20)*(1.0631^i);
      SceInput(3,i)=inputSampledata(3,20)*(1.0872^i);
      SceInput(4,i)=inputSampledata(4,20)*(1.2044^i);
      SceInput(5,i)=inputSampledata(5,20)*(1.2326^i);
      SceInput(6,i)=inputSampledata(6,20)*(1.0605^i);
      SceInput(7,i)=2*(1.01^i);
      SceInput(8,i)=42*(1.02^i);
      SceInput(9,i)=inputSampledata(9,20)*(1.1426^i);
      SceInput(10,i)=inputSampledata(10,20)*(1.017^i);
      SceInput(11,i)=inputSampledata(11,20)*(1.0205^i);
      SceInput(12,i)=inputSampledata(12,20)*(1.1336^i);
      SceInput(13,i)=inputSampledata(13,20)*(1.1599^i);
      SceInput(14,i)=inputSampledata(14,20)*(1.1783^i);
 end

for j=1:20;
   for i=1:14;
      recalldata(i,j)=inputSampledata(i,j);
   end
end

for j=21:30;
   for i=1:14;
      recalldata(i,j)=SceInput(i,j-20)  
  end
end

[alterinput,mininput,maxinput] = premnmx(recalldata);


%training
fvalue=sim(gwwnet,alterinput);

%data offset (converting the network output data to it original unit)
forecastvalue = postmnmx(fvalue,minoutput,maxoutput);

%plot

waitforbuttonpress;
clf;
time=1978:1:2007;
time1=1978:1:1997;
plot(time,forecastvalue,'o',time1,outputSampledata,'-');
legend('預(yù)測(cè)曲線','實(shí)際曲線');
title('客運(yùn)量曲線');
xlabel('時(shí)間');ylabel('公交客運(yùn)量');