spatiotemporal_pad0.m

%cell_ids = [348 400 498 519];
%cell_ids = [45 54 66 69];
cell_ids = [498];
numZoness = [6]; % Number of neighbour cells to incorporate
maxTimestamp = 1440; % Number of time periods in a day(depend on sampling interval)
samplingSteps = [10];

% Construct distance matrix between cells
% columns are cellid, centroid_x, centroid_y, etc
celldata = csvread('./Data/cell_info_1.0km_priority.csv',2,0);
numCells = rows(celldata);
distance = zeros(numCells);

for i = 1:numCells
    if any(i == cell_ids)
      for j = 1:numCells          
          if i == j
              ang(i,j) = 10000;
              distance(i,j) = 10000;
          else
              distance(i,j) = sqrt((celldata(i,2)-celldata(j,2))^2 + (celldata(i,3) - celldata(j,3))^2);
              dx = celldata(j,2) - celldata(i,2);
              dy = celldata(j,3) - celldata(i,3);
              h = distance(i,j);
              if dy > 0
                ang(i,j) = acos(dx/h);
              else
                ang(i,j) = 2*pi() - acos(dx/h);
              end      
          end   
      end
    end
end

mseall = [];
for q = 1:columns(cell_ids)
  mse = zeros(columns(numZoness),columns(samplingSteps));
  % Retrieve indexes of the nearest numZones cells
  minIndex = zeros(numCells,max(numZoness));
  for i = 1:numCells
      if any(cell_ids==i)
        [~,d] = sort(distance(i,:));
        minIndex(i,:) = d(1:max(numZoness));
      else
        minIndex(i,:) = zeros(1,max(numZoness));
      end
  end
  neighbourcell = minIndex(cell_id,:)
  for r = 1:columns(numZoness)
    for s = 1:columns(samplingSteps)

        cell_id = cell_ids(q);
        numZones = numZoness(r);
        samplingStep = samplingSteps(s);
        numTimePeriod = maxTimestamp / samplingStep;

        numDays = 4; % Number of days of data
        trainDays = 3; % Number of days used as training data
        weekend = 0;
        
        traj = [];
        time = [];
        trajNeigh = [];
        time2 = [];
        dow = [];
  
        % Read in Data
        for date = 1:numDays

            if (date < 10) 
                d = strcat('0',num2str(date));
            else
                d = num2str(date);
            end
            
            % columns are cell_id,time_id,start_time,end_time,num_traj,inflow,outflow,avg_traveldist_km,avg_traveltime_min,avg_speed_kph
            data = csvread(strcat('./Data/node_measure_1km_1min_1min_201507',d,'.csv'),2,1);
            
            % Take out weekends (optional depending on needs)
            if weekday(strcat('2015-07-',d)) == 1 || weekday(strcat('2015-07-',d)) == 7
                weekend = weekend + 1;
                if (date <= trainDays+weekend)
                    trainDays = trainDays - 1;
                end
                numDays = numDays - 1;
                continue
            end
            
            % Call GetTrajs to get the flows of the neighbouring cells
            temp = [];       
            for i = 1:numZones
                dataneighbour = [];
                while (rows(dataneighbour) == 0)
                    dataneighbour = data(data(:,1) == neighbourcell(i),:);
                end
                [a, ~] = GetTrajs(dataneighbour,d,maxTimestamp,samplingStep);
                temp = [temp a];
            end 
            trajNeigh = [trajNeigh; temp];
            % Get information of the interested cell
            data = data(data(:,1) == cell_id,:);
            [a, b] = GetTrajs(data,d,maxTimestamp,samplingStep);
            traj = [traj; a];
            time = [time; b];
            
        end

        if (numDays == trainDays)
            disp("Cannot predict weekends")
            break
        end
        
        coordx = repmat(celldata(cell_id,2),numTimePeriod*numDays,1);
        coordy = repmat(celldata(cell_id,3),numTimePeriod*numDays,1);
        angles = repmat(pi()/2, numTimePeriod*numDays,1);
        for i = 1:numZones
            coordx = [coordx;repmat(celldata(neighbourcell(i),2),numTimePeriod*numDays,1)];
            coordy = [coordy;repmat(celldata(neighbourcell(i),3),numTimePeriod*numDays,1)];
            angles = [angles;repmat(ang(cell_id, neighbourcell(i)),numTimePeriod*numDays,1)];
        end
        
        x_temporal = repmat(time,numZones+1,1); % Prediction only works on continuous time stamps (1,2,3,...) not (0,5,10,...)
        y = traj;
        for i = 1:numZones
            y = [y; trajNeigh(:,i)];
        end
        x = [coordx coordy angles x_temporal];
        x = Normalize(x);
        xtrain = x(1:numTimePeriod*trainDays,:);
        ytrain = y(1:numTimePeriod*trainDays);
        for i = 1:numZones
            xtrain = [xtrain; x(i*numTimePeriod*numDays+1:i*numTimePeriod*numDays+numTimePeriod*trainDays, :)];
            ytrain = [ytrain; y(i*numTimePeriod*numDays+1:i*numTimePeriod*numDays+numTimePeriod*trainDays)];
        end
        
        meanfunc = {};
        covfunc = {@covProd,{{@covMask,{[1 1 1 0],@covSEiso}},{@covMask,{[0 0 0 1],{'covProd',{'covPeriodic',{'covProd',{'covRQiso','covPeriodic'}}}}}}}};
        %covfunc = {@covMask,{[0 0 1],{'covProd',{'covPeriodic',{'covProd',{'covRQiso','covPeriodic'}}}}}};
        %covfunc = {@covProd,{{@covMask,{[1 1 0],@covSEiso}}}};
        likfunc = @likGauss;
        hyp = struct('mean',[],'cov', [0 0 0 0 0 0 0 0 0 0 0], 'lik', -1);
        hyp2 = minimize(hyp, @gp, -100, @infExact, meanfunc, covfunc, likfunc, xtrain, ytrain);

        % Predict and Plot
        plotdays = linspace(trainDays, numDays, numDays-trainDays); % the dates to predict
        for i = 1:rows(plotdays)
            xtest = x(numTimePeriod*(plotdays(i)-1)+1:numTimePeriod*plotdays(i), :);
            ytest = y(numTimePeriod*(plotdays(i)-1)+1:numTimePeriod*plotdays(i));
            for j = 1:numZones
                xtest = [xtest;x(numTimePeriod*(plotdays(i)-1)+j*numTimePeriod*numDays+1:numTimePeriod*plotdays(i)+j*numTimePeriod*numDays, :)];
                ytest = [ytest;y(numTimePeriod*(plotdays(i)-1)+j*numTimePeriod*numDays+1:numTimePeriod*plotdays(i)+j*numTimePeriod*numDays)];
            end

            [mu, s2] = gp(hyp2, @infExact, meanfunc, covfunc, likfunc, xtrain, ytrain, xtest);
            mse(r,s) = sum((mu(1:numTimePeriod)-ytest(1:numTimePeriod)).^2)/numTimePeriod;
            fig = figure;
            plot(ytest(1:numTimePeriod),'b', 'LineWidth',2);
            hold on;
            plot(mu(1:numTimePeriod),'r');
            legend('Actual Flow', 'Prediction');
            hold off;
            saveas(fig, strcat('07-',num2str(plotdays(i)),'-',num2str(cell_id),'-1km-',num2str(numZones),'-',num2str(samplingStep),'-angle.png'));
            close;
        end
    end
  end
  mseall = [mseall mse];
end