%{

Description:
This script opens example data and estimates the scale size u* and T* from
given inputs. After that, this program creates four simple plots to 
visualise the results given by the neural network.

Required files: 'DE-KaN.dat', 'net_6_3_2.mat' (located in same directory as
    this run.m script.)

Feature explanation:
'dTheta_dz'     -> difference quotient of potential temperature theta(z)
                    dTheta_dz = (theta(z2)-theta(z1))/z2-z1 , z2>z1
'dU_dz'         -> difference quotient of horizontal wind velocity u(z)
                    dU_dz = (u(z2)-u(z1))/z2-z1 , z2>z1
'Tmean'         -> mean potential temperature of two levels
                    Tmean = 0.5*(theta(z2)+theta(z1))
'cveg'          -> vegetation classifier
                    cveg = 1 if vegetation is low (e.g. grasslands)
                    cveg = 2 if vegetation is tall (e.g. forest)
'dThetadz_dUdz' -> ratio of difference quotients
                    dThetadz_dUdz = dTheta_dz / dU_dz
'umean'         -> mean wind velocity of two levels
                    umean = 0.5*(u(z2)+u(z1))

output layer:
netForecast(:,1)-> network's estimate for u*
netForecast(:,2)-> network's estimate for T*

%}

% set import file name
fileName = 'DE-KaN.dat';
% set input features
featureNames = {'dTheta_dz', 'dU_dz','Tmean','cveg','dThetadz_dUdz','umean'};
% set targets
targetNames = {'USTAR', 'TSTAR'};
% load inputs, targets and date array
[inputs, targets, date] = importFunction(fileName,featureNames,targetNames);

% load net
[net] = loadNet('6_3_2','./');

% do forecast
netForecast = net(inputs')';

% simple plot
plotScatter(netForecast, targets);
plotTimeseries(netForecast, targets, date);


function [] = plotTimeseries(netForecast, targets, date)
% plot u*
figure()
plot(date, targets(:,1),'DisplayName','observations');
hold on;
plot(date, netForecast(:,1),'DisplayName','neural network');
xlabel('date');
ylabel('u* [ms-1]');
legend();
% plot T*
figure()
plot(date, targets(:,2),'DisplayName','observations');
hold on;
plot(date, netForecast(:,2),'DisplayName','neural network');
xlabel('date');
ylabel('T* [K]');
legend();
end

function [] = plotScatter(netForecast, targets)
% plot u*
figure();
scatter(targets(:,1),netForecast(:,1));
hold on;
plot([0,1],[0,1])
xlabel('u* observations [ms-1]');
ylabel('u* neural network [ms-1]');
% plot T*
figure();
scatter(targets(:,2),netForecast(:,2));
hold on;
plot([-0.8,0.3],[-0.8,0.3])
xlabel('T* observations [K]');
ylabel('T* neural network [K]');
end

function [returnNet] = loadNet(netName, path)
% load net
tmp = load(sprintf([path,'net_%s.mat'], netName));
returnNet = tmp.net;
% show net structure
view(returnNet);
end

function [inputValues, targetValues, dateArray] = importFunction(fileName, inputNames, targetNames)

% open data file
dataSet = readtable(sprintf(['./',fileName]));

% check if all requested input variables are available
[bool, loc] = ismember(inputNames, dataSet.Properties.VariableNames);
if sum(bool)<length(bool)
    msg = mat2str(inputNames(~bool));
    error('Missing input variables in given file: %s', msg);
end
% set inputs
inputValues = dataSet{:,loc};

% check if all requested target variables are available
[bool, loc] = ismember(targetNames, dataSet.Properties.VariableNames);
if sum(bool)<length(bool)
    msg = mat2str(targetNames(~bool));
    error('Missing target variables in given file: %s', msg);
end
% set targets
targetValues = dataSet{:,loc};

% extract date array
dateArray = dataSet{:,1};

end