close all
load output5.dat;
optfreq2;
data2 = output5;
%% DETERMINING DATES

 startDay = '13';                %Edit desired date and length of sample here
 startMonth = 'Dec';
 startYear = 2014;

 endDay = '1';
 endMonth = 'Jan';
 endYear = 2016;

 predDay = '1';
 predMonth = 'Jan';
 predYear = 2017;

startDay = '10';                %Edit desired date and length of sample here
startMonth = 'Feb';
startYear = 2016;

endDay = '27';
endMonth = 'Feb';
endYear = 2017;

predDay = '1';
predMonth = 'May';
predYear = 2017;

singlepredDay = '20';
singlepredMonth = 'Apr';
singlepredYear = 2017;

month = strcat(startDay, '-', startMonth, '-');
startDate = (datenum(strcat(month, num2str(startYear+0)))-725738)*24*6;

endmonth = strcat(endDay, '-', endMonth, '-');
endDate = (datenum(strcat(endmonth, num2str(endYear+0)))-725738)*24*6;

predNum = strcat(predDay, '-', predMonth, '-');
predDate = (datenum(strcat(predNum, num2str(predYear+0)))-725738)*24*6;

singlepredNum = strcat(singlepredDay, '-', singlepredMonth, '-');
singlepredDate = (datenum(strcat(singlepredNum, num2str(singlepredYear+0)))-725738)*24*6;

%% DEFINING IMPORTANT VARIABLES

fourierVec = data2(startDate:endDate);

if predDate < length(data2)
    predVec = data2(endDate:predDate);
end

predLength = length([endDate:predDate]);

lfV = length(fourierVec);

unit = 1/(6*24);
t = 0:unit:(lfV-1)*unit;

ampFilterMax = 200;            %Edit desired amplitude filtering here

frequencyDomain = 0: 1/(length(fourierVec)/(6*24)) : (length(fourierVec)-1)/(length(fourierVec)/(6*24));

%% CALCULATING FOURIER TRANSFORM AND FREQUENCIES

figure

fourierTransform = fft(fourierVec*24*6/lfV); %Creating and filtering the Fourier Transform

%% METHOD 1 (COMMENT ONE)

% fourierTransform(abs(fourierTransform) < ampFilterMax) = 0;
% freqVec = find(abs(fourierTransform) > 0);
% realfreqVec1 = [];
%
% for i =1:length(freqVec)
%     if(freqVec(i) < (1/2)*lfV)
%         realfreqVec1 = [realfreqVec1 2*pi*freqVec(i)/(lfV*unit)];
%     end
% end

%% METHOD 2 (COMMENT ONE)

% [peaks, freqVecAll] = findpeaks(abs(fourierTransform), 'MinPeakHeight', 100);
% freqVecAll = 2*pi*freqVecAll/(lfV*unit);
% realfreqVec2 = freqVecAll(1:1/2*length(freqVecAll))';

%%

realfreqVec = optfreq3;

display(length(realfreqVec));
absFourier = abs(fourierTransform);

%% MAKING APPROXIMATION

f=fourierVec-mean(fourierVec);
abSolve = [];
for i = 1:length(realfreqVec)
    abSolve=[abSolve; fminsearch(@(x) sum((f(1:length(f))' -x(1)*cos(realfreqVec(i)*t(1:length(f))) - x(2)*sin(realfreqVec(i)*t(1:length(f)))).^2),[-100,-100])];
    display(i)
end

a=[];
b=[];

for i=1:length(realfreqVec)
    a = [a abSolve(i,1)];
    b = [b abSolve(i,2)];
end

sumCosApprox = 0;

for i=1:length(a)
    sumCosApprox = sumCosApprox +a(i)*cos(realfreqVec(i)*t)+b(i)*sin(realfreqVec(i)*t);
end
sumCosApprox = sumCosApprox+mean(fourierVec);

t = 0:unit:(predLength-1)*unit;
sumCos = 0;

for i=1:length(a)
    sumCos = sumCos +a(i)*cos(realfreqVec(i)*t)+b(i)*sin(realfreqVec(i)*t);
end
sumCos = sumCos+mean(fourierVec);
%% PLOTTING THINGS

plot(frequencyDomain, absFourier);
axis([0 10 0 max(absFourier)]);
title('Filtered Fourier Transform');
xlabel('1/T');

figure

plot(data2(startDate:endDate), 'Color', 'r');
hold on
plot(sumCosApprox, 'Color', 'b')
title('Approximation using sines and cosines');
legend('Original Data', 'Approximation');
xlabel('Index of approximated data points');
ylabel('Water height in cm above NAP');
if max(data2)> max(sumCosApprox)
    if min(data2) < min(sumCosApprox)
        axis([0 lfV min(data2) max(data2)]);
    else
        axis([0 lfV min(sumCosApprox) max(data2)]);
    end
else
    if min(data2) > min(sumCosApprox)
        axis([0 lfV min(sumCosApprox) max(sumCosApprox)]);
    else
        axis([0 lfV min(data2) max(sumCosApprox)]);
    end
end
%% PREDICTING KNOWN DATA
errorSumCos = sumCos;

if predDate < length(data2)
    predVec = data2(endDate:predDate);
    errorApproximation = errorSumCos(1:length(errorSumCos)) - predVec(1:length(predVec))';

    figure

    plot(predVec);
    hold on
    plot(errorSumCos);

    axis([0 predLength min(predVec) max(predVec)])
    title('Prediction of known data');
    xlabel('Index of predicted data points');
    ylabel('Predicted water height in cm above NAP');
    legend('Original Data','Prediction');

    figure
    plot(errorApproximation);
    axis([0 predLength min(errorApproximation) max(errorApproximation)]);
    title('Error of predicting known data');
    xlabel('Index of predicted data');
    ylabel('Error in cm of predicted data')

    predictionTime = sumCos;
    predictionTime(find(sumCos<=160))=0;
    predictionTime(find(sumCos>160))=1;

    dataTime = data2(endDate:predDate);

    dataTime(find(dataTime<=150))=0;
    dataTime(find(dataTime>150))=1;

    figure

    predictionMistake = predictionTime(1:predLength)-dataTime(1:predLength)';
    plot(predictionMistake);
    axis([0 predLength -1.5 1.5]);
    title('Prediction accuracy of known data');
    xlabel('Index of predicted data points');
    ylabel('Error in prediction');

    predErrors = find(predictionMistake ~= 0);
    percentageWrong = length(predErrors)/predLength;
    display(percentageWrong);

    predErrorsUp = find(predictionMistake > 0);
    percentageUp = length(predErrorsUp)/predLength;
    display(percentageUp);

    predErrorsDown = find(predictionMistake < 0);
    percentageDown = length(predErrorsDown)/predLength;
    display(percentageDown);

    figure

    errorMean = mean(abs(errorApproximation));
    errorDeviation = std(errorApproximation);
    display(errorMean);
    display(errorDeviation);
    histogram(errorApproximation);
    title('Distribution of error of known data');
    xlabel('Error in cm');
    ylabel('Amount of data points');

    figure

    [F, y] = ecdf(errorApproximation);
    plot(y, F);
    title('Cumulative Distribution of error of known data');
    xlabel('Error in cm');
    ylabel('Cumulative probability');

    error = [];
    approx = [];
    for i=1:length(predVec)
        if predVec(i) >= 150
            approx = [approx errorSumCos(i)];
            error = [error errorApproximation(i)];
        end
    end
    figure
    histogram(error);
    title('Distribution of error looking at prediction of high water levels');
    xlabel('Error in cm');
    ylabel('Amount of data points');

    figure

    [G, z] = ecdf(error);
    plot(z, G);
    title('Cumulative Distribution of prediction at high water levels');
    xlabel('Error in cm');
    ylabel('Cumulative probability');

end

%% PREDICTING UNKNOWN DATA

if predDate > length(data2)

    figure

    plot(errorSumCos);
    axis([0 predLength min(errorSumCos) max(errorSumCos)]);
    title('Prediction');
    %axis([0 predictionLength min(errorSumCos) max(errorSumCos)]);
    xlabel('Index of predicted data points');
    ylabel('Predicted water height in cm above NAP');

    prediction_150 = errorSumCos;
    prediction_150(find(errorSumCos > 150)) = 1;
    prediction_150(find(errorSumCos <= 150)) = 0;

    figure

    plot(errorSumCos)
    axis([singlepredDate-endDate singlepredDate-endDate+24*6 min(errorSumCos) max(errorSumCos)]);
    title('Prediction of April 20th');
    xticks(singlepredDate-endDate:12:singlepredDate-endDate+24*6);
        xticklabels({'00:00', '02:00', '04:00', '06:00', '08:00', '10:00', '12:00', '14:00', '16:00', '18:00', '20:00', '22:00', '24:00'});
    xlabel('Time (hours)');
    ylabel('Predicted water height in cm above NAP');

    figure

    plot(prediction_150)
    axis([singlepredDate-endDate singlepredDate-endDate+24*6 -1.5 1.5]);
    xticks(singlepredDate-endDate:12:singlepredDate-endDate+24*6);
        xticklabels({'00:00', '02:00', '04:00', '06:00', '08:00', '10:00', '12:00', '14:00', '16:00', '18:00', '20:00', '22:00', '24:00'});
    title('Predicted intervals at which a ship can pass on April 20th');
    xlabel('Time (hours)');
    ylabel('Predicted water height in cm above NAP');

    if length(F)>0
        figure
        plot(y,F)

        figure
        errorplot = [];
        for i = 1:length(errorSumCos)
            errorplot(i) = max(F(y<=150-errorSumCos(i)));
        end
        plot(errorplot)
        hold on
        plot(prediction_150)
        axis([singlepredDate-endDate singlepredDate-endDate+24*6 0 1]);
        xticks(singlepredDate-endDate:12:singlepredDate-endDate+24*6);
        xticklabels({'00:00', '02:00', '04:00', '06:00', '08:00', '10:00', '12:00', '14:00', '16:00', '18:00', '20:00', '22:00', '24:00'});
        xlabel('Time (hours)');
        ylabel('Probability that a ship can pass through on April 20th');
    end

    if length(G)>0
        figure
        plot(z,G)

        figure
        errorplot = [];
        for i = 1:length(errorSumCos)
            errorplot(i) = max(G(z<=150-errorSumCos(i)));
        end
        plot(errorplot)
        hold on
        plot(prediction_150)
        axis([singlepredDate-endDate singlepredDate-endDate+24*6 0 1]);
        xticks(singlepredDate-endDate:12:singlepredDate-endDate+24*6);
        xticklabels({'00:00', '02:00', '04:00', '06:00', '08:00', '10:00', '12:00', '14:00', '16:00', '18:00', '20:00', '22:00', '24:00'});
        xlabel('Time (hours)');
        ylabel('Probability that a ship can pass through on April 30th - High points');
    end
end