clc, clear all
getplotlyoffline('https://cdn.plot.ly/plotly-latest.min.js')

%% File Retrieval 
c1Cell = {'1.7','1.9','2.1','2.3','2.5','2.7','2.9','3.2'};
c1Cell = strcat(c1Cell,'e-3');
doCondHorz = 0;  doTermHorz = 0; 
doFactors = [1,3,7];
horizonDir = 'CSV_Data-Planned_01-Aug-2012_DOY213';      % Specific horizon directory to look at (3 weeks of data, determined 1 day of planned policy).
horizonPrefix = 'HorizonLogs_';

fprintf('[doCondHorz, doTermHorz] = [%d,%d]\n',doCondHorz,doTermHorz)
fprintf('For penalty factors: ')
fprintf('%s, ', c1Cell{doFactors})
fprintf('\n===============================================\n')
for penFactor = doFactors
    penaltyDir = [horizonPrefix,c1Cell{penFactor}];                            % Directory where all horizons are located for this penalty
    plannedDay = horizonDir(18:28);
    curdir = [pwd, '/', penaltyDir, '/', horizonDir]; 
    csvStruct = struct2cell(dir([curdir,'/*.csv']));
    csvNames = natsort(csvStruct(1,:));
    horizonData{penFactor} = [];
    fprintf('Sorting chronologically...\n')
    
    if doCondHorz == 0, startNum = 8; condOffset = 0; else startNum = 1; condOffset = -7; end 
    if doTermHorz == 0, endNum = 14; termOffset = 7; else endNum = length(csvNames); termOffset = 14; end
    for dayNum = startNum:endNum
           filename = [curdir,'/',csvNames{dayNum}];
           fprintf('Reading file: %s\n',filename);
           thisHorizon = csvread(filename,1,0);
           horizonData{penFactor} = [horizonData{penFactor};thisHorizon];
    end 
    fprintf('Done for %s. Number of files read: %d. For: horizonData{%d} \n\n',c1Cell{penFactor},length(csvNames),penFactor)
end 
allLabels = readHorizonHeader(filename, 1, 1); % All headers are the same anyway

% Column reading info
%=============================
readColumns = [[1:5],[12:23]];           % Columns of horizon data we're interested in 
readLabels = allLabels(readColumns);     % Their labels/header for legends
t1 = addtodate(datenum(plannedDay) , condOffset, 'day');   
t2 = addtodate(datenum(plannedDay) , termOffset, 'day');

disp('Truncating for labels.')
for x = doFactors
    horizonData{x} = horizonData{x}(:,readColumns);
end 
%********************

%% Plot 1: Tset and Person-Weighted Average Temperature
% Alongside this, the plot shows:
% -  $/kWh for the date range [t1,t2] (convert from pjmLMP.csv)
% -  Occupancy fraction (defined here)
% -  Total HVAC kW power usage
% -  Person-weighted PMV
interval = 5/24/60;
t = convertDate(t1:interval:t2); t = t(1:end-1);  % Hack off the end
%t = convertDate(linspace(t1,t2,size(horizonData{1},1)));
fig1 = figure('units','normalized','outerposition',[0 0 1 1]);
numSubplots = 5; 

% Get the occupancy information
% MF: (1 person / 18.58 m2)(volume per floor in each zone)(12 floors)(1/2.74m ceiling) 
% MF: ordering here and in the paper is N,S,E,W,Core, but 
% in EnergyPlus and NetPenalty.m it is S,E,N,W,Core.
totalPop = (1/18.58).*[860,860,554,554,6949].*(12/2.74);
fractionWeekday = [0, 8, 11, 12 ,13 ,14 ,18, 19, 24;   % Top row: Until: xx:00 (military time)
                   0, 0, 1,  .8, .4, .8, 1,  .5, 0];   % Bottom row: coefficient (%) of normal population present.fractionScheduleWeekend = [0,0,0.1,0.3,0.3,0,0,0.3,0.3,0.2,0.2,0.1,0];
fractionWeekend = [0, 8, 9,   11  ,13 ,16 ,18, 19, 24;   % Top row: Until: xx:00 (military time)
               0, 0, 0.2, 0.3, 0, .3, .2,  .1, 0];   % Bottom row: coefficient (%) of normal population present.
fractionSchedule = [];
totalDays = abs(endNum-startNum); % Counts "day 0"
for day = 0:totalDays
    today = addtodate(t1,day,'day');
    [~,dow] = weekday(today);
    if ismember(dow, {'Sat','Sun'})
        fractionSchedule = [fractionSchedule; genOccupyCoeff(fractionWeekend,12);];
    else
        fractionSchedule = [fractionSchedule; genOccupyCoeff(fractionWeekday,12);];
    end 
end 

% Time to do plotting.
%-------------------------
% Subplot 1: temperatures 'n stuff
ax(1) = subplot(numSubplots,1,1);
hold all
trace1 = {};
colormap = {'4c72b0','c44e52','55a868','8172b2','ccb974','64b5cd','cd8464','cc64cd'};
for i = doFactors
    TzoneAvg = sum(horizonData{i}(:,7:11).*kron(totalPop,ones(288*(totalDays+1),1)),2)./sum(totalPop);
    plot(t,horizonData{i}(:,1))         % Building-wide setpoints for each factor
    plot(t,TzoneAvg)                    % Weighted zone temperatures per factor
    TsetStruct = struct('x',t,'y',horizonData{i}(:,1),'xaxis','x','type', 'scatter','yaxis','y5','name',['Tset* ',c1Cell{i}]);
    TsetStruct.line.color = colormap{i};
    trace1 = [trace1,TsetStruct];
    zoneStruct = struct('x',t,'y',TzoneAvg,'xaxis','x','type', 'scatter','yaxis','y5','name',['E[Tzone] ',c1Cell{i}]);
    zoneStruct.line.dash = 'dot';
    zoneStruct.line.color = colormap{i};
    trace1 = [trace1,zoneStruct];
end 
plot(t,horizonData{i}(:,6))   % Single ambient temperature for all penalties
TambStruct = struct('x',t,'y',horizonData{i}(:,6),'xaxis','x','type', 'scatter','yaxis','y5','name','Tamb');
TambStruct.line.color = '212121'; 
TambStruct.opacity = 0.5;
trace1 = [trace1,TambStruct];


% Subplot 2: $/kWh
ax(2) = subplot(numSubplots,1,2);
hold all
doy_t1 = datevec(t1); doy_t1(1) = 2011; doy_t1 = datevec2doy(doy_t1); % Force no leap year
doy_t2 = datevec(t2); doy_t2(1) = 2011; doy_t2 = datevec2doy(doy_t2)-1;
lmpTab = importLambda('./pjmLMP.xlsx');
lmp = lmpTab(doy_t1:doy_t2,:)';    % Go to $/kWh
lmpVec = kron(lmp(:),ones(12,1));
plot(t,lmpVec)
trace2 = {};
trace2 = [trace2,struct('x',t,'y',lmpVec,'yaxis','y4','xaxis','x','type', 'scatter','name','$/MWh')];

% Subplot 3: fractional schedule
ax(3) = subplot(numSubplots,1,3);
hold all
stairs(t,fractionSchedule)
trace3 = {};
trace3 = [trace3,struct('x',t,'y',fractionSchedule,'yaxis','y3','xaxis','x','type', 'scatter','name','Light/Equip/Person Fraction')];

% Subplot 4: power usage (kW)
ax(4) = subplot(numSubplots,1,4);
hold all
trace4 = {};
for i = doFactors
    plot(t,horizonData{i}(:,12))    
    trace4 = [trace4,struct('x',t,'y',horizonData{i}(:,12),'yaxis','y2','xaxis','x','type', 'scatter','name',['Power Usage (W) ',c1Cell{i}])];
end 

% Subplot 5: person-weighted PMV
ax(5) = subplot(numSubplots,1,5);
hold all
trace5 = {};
for i = doFactors
    PMVweight = sum(horizonData{i}(:,13:17).*kron(totalPop,ones(288*(totalDays+1),1)),2)./sum(totalPop).*(fractionSchedule > 0);
    PMVweight(PMVweight == 0) = NaN;
    plot(t,PMVweight)                    % Weighted zone temperatures per factor
    trace5 = [trace5,struct('x',t,'y',PMVweight,'yaxis','y','xaxis','x','type', 'scatter','name',['Occupied PMV ',c1Cell{i}])];
end 

for handle = 1:5
subplot(numSubplots,1,handle)
thisax = gca;
set(thisax,'XTick',convertDate(t1:t2))
end

xaxis =  struct('domain', [0, 1],'nticks',totalDays,'ticks','inside','type',date,'side','bottom','title','Horizon Date'); % let the first x-axis span the first 45% of the plot width
yaxis =  struct('domain', [0.34, 0.41],'nticks',5, 'anchor','x','linecolor', '#636363','linewidth', 6,'gridcolor', '#bdbdbd','zeroline', true,'title','PMV'); % and let the first y-axis span the first 45% of the plot height
yaxis2 =  struct('domain', [0.43, 0.50],'nticks',6, 'anchor','x','linecolor', '#636363','linewidth', 6,'gridcolor', '#bdbdbd','zeroline', true,'title','Wattage'); % and let the first y-axis span the first 45% of the plot height
yaxis3 =  struct('domain', [0.52, 0.59],'nticks',8, 'anchor','x','linecolor', '#636363','linewidth', 6,'gridcolor', '#bdbdbd','zeroline', true,'title','Fraction'); % and let the first y-axis span the first 45% of the plot height
yaxis4 =  struct('domain', [0.61, 0.68],'nticks',5, 'anchor','x','linecolor', '#636363','linewidth', 6,'gridcolor', '#bdbdbd','zeroline', true,'title','$/MWh'); % and let the first y-axis span the first 45% of the plot height
yaxis5 =  struct('domain', [0.695, 1],'nticks',6, 'anchor','x','linecolor', '#636363','linewidth', 6,'gridcolor', '#bdbdbd','zeroline', true,'title','Temperature (C)'); % and let the first y-axis span the first 45% of the plot height

xaxis = struct('domain', [0, 1], ...
      'showgrid', true, ...
      'zeroline', true, ...
      'showline', true, ...
      'gridcolor', '#bdbdbd', ...
      'gridwidth', 2, ...
      'zerolinecolor', '#969696', ...
      'zerolinewidth', 4, ...
      'linecolor', '#636363', ...
      'linewidth', 6, ...
      'type', 'date', 'nticks',14, 'title','Simulated Time/Date'); 
layout = struct('xaxis',xaxis, 'yaxis',yaxis, 'yaxis2',yaxis2, 'yaxis3',yaxis3, 'yaxis4',yaxis4,'yaxis5',yaxis5);

layout.width = 2200; % required
layout.height = 1200; % required

p = plotlyfig; close all
p.data = {trace1{:},trace2{:},trace3{:},trace4{:},trace5{:}};
p.layout = layout;
p.PlotOptions.FileName = 'Plot_1';
html_file = plotlyoffline(p);
%r = plotly(p.data, struct('layout', layout));
disp(['View file: ', html_file])

%% Plot 2: Running Sums
clc,clear,close all
csvStruct = struct2cell(dir([pwd,'/Full_Logs/*.csv']));
csvNames = natsort(csvStruct(1,:));
fullLogs = {};
c1Cell = {'1.7','1.9','2.1','2.3','2.5','2.7','2.9','3.2','10','20','30','40'};
c1Cell = strcat(c1Cell,'e-3');
wtp = cellfun(@str2num,c1Cell);
for i = 1:length(csvNames)
        fullLogs = [fullLogs, openFullLog([pwd,'/Full_Logs/',csvNames{i}])];
        fullLogs{i} = [fullLogs{i},cumsum(fullLogs{i},1)];   % Attach cumulative sums
        sumPelec(:,i) = fullLogs{i}(:,5);
        sumabsDiscomfort(:,i) = fullLogs{i}(:,6)./wtp(i);
end

mycolors = distinguishable_colors(length(c1Cell));
figure
dumb75Pelec = kron(sumPelec(:,end-2),ones(1,length(c1Cell)));
dumb75Disc = kron(sumabsDiscomfort(:,end-2),ones(1,length(c1Cell)));
fracPelecSave = sumPelec-dumb75Pelec;
fracDiscSave = sumabsDiscomfort-dumb75Disc;

% fracPelecSave(:,length(c1Cell)-1) = fracPelecSave(:,length(c1Cell)-1) == -1;
% fracDiscSave(:,length(c1Cell)-1) = fracDiscSave(:,length(c1Cell)-1) == -1;

for i = 1:length(c1Cell)
    plot(1:31,fracPelecSave(:,i),'color',mycolors(i,:),'linewidth',3)
    hold all
end 
legend(c1Cell)
figure
for i = 1:length(c1Cell)
    plot(1:31,fracDiscSave(:,i),'color',mycolors(i,:),'linewidth',3)
    hold all
end 
legend(c1Cell)
plotlycolors = colorGradient([0,0,130],[51,204,255],8);
plotlycolors = rgb2hex(plotlycolors);
plotlycolors = [plotlycolors; ['#ff0000';'#ff0000';'#00e600';'#cc6699']];
plotlynames = {c1Cell{1:end-4},'Bounded w/ Term','Bounded w/o Term','Occupied 75F','Occupied Tideal=72.5F'};
for i = [1:length(c1Cell)]
   thisStruct = struct('x',1:31,'y',fracPelecSave(:,i),'type','scatter','name',plotlynames{i},'yaxis','y2');
   thisStruct.line.color = plotlycolors(i,:);
   if ismember(i,2:8), thisStruct.fill = 'tonexty';
   else thisStruct.fill = 'none'; end 
   if ismember(i,[1:8]), thisStruct.line.dash = 'dot'; 
   elseif i == 9, thisStruct.line.dash = 'dash'; end 
   trace1{i} = thisStruct;
end 

for i = [1:length(c1Cell)]
   thisStruct = struct('x',1:31,'y',fracDiscSave(:,i),'type','scatter','name',plotlynames{i},'yaxis','y');
   thisStruct.line.color = plotlycolors(i,:);
   if ismember(i,2:8), thisStruct.fill = 'tonexty';
   else thisStruct.fill = 'none'; end 
   if ismember(i,[1:8]), thisStruct.line.dash = 'dot';
   elseif i == 9, thisStruct.line.dash = 'dash'; end 
   trace2{i} = thisStruct;
end 

fig2 = figure('units','normalized','outerposition',[0 0 1 1]);
xaxis = struct('domain', [0, 0.9], ...
      'showgrid', true, ...
      'zeroline', true, ...
      'showline', true, ...
      'gridcolor', '#bdbdbd', ...
      'gridwidth', 2, ...
      'zerolinecolor', '#969696', ...
      'zerolinewidth', 4, ...
      'linecolor', '#636363', ...
      'linewidth', 6, ...
      'nticks',31, ...
      'title', 'Day of August', ...
      'range',[1,31]); 
yaxis =  struct('domain', [0, 0.47],'nticks',15, 'linecolor', '#636363','linewidth', 6, ...
            'gridcolor', '#bdbdbd','zeroline', true,'title','Cumulative Absolute Discomfort Difference','anchor','x','range',[-160e3,40e3]);          

yaxis2 =  struct('domain', [0.53, 1],'nticks',15, 'linecolor', '#636363','linewidth', 6, ...
                'gridcolor', '#bdbdbd','zeroline', true,'title','Cumulative Pelec Difference','range',[-150,250]); 

            
layout = struct('xaxis',xaxis, 'yaxis',yaxis,'yaxis2',yaxis2);
layout.width = 1500; % required
layout.height = layout.width/1.2; % required

p1 = plotlyfig; close all
p1.data = [trace1,trace2];
p1.layout = layout;
p1.PlotOptions.FileName = 'Plot_2';
html_file = plotlyoffline(p1);
r = plotly(p1.data, struct('layout', layout));
%% Plot 3: Pelec vs Absolute Discomfort
clc,clear,close all
csvStruct = struct2cell(dir([pwd,'/Full_Logs/*.csv']));
csvNames = natsort(csvStruct(1,:));
fullLogs = {};
c1Cell = {'1.7','1.9','2.1','2.3','2.5','2.7','2.9','3.2','10','20','30','40'};
c1Cell = strcat(c1Cell,'e-3');
wtp = cellfun(@str2num,c1Cell);
for i = 1:length(csvNames)
    fullLogs = [fullLogs, openFullLog([pwd,'/Full_Logs/',csvNames{i}])];
    fullLogs{i} = [fullLogs{i},cumsum(fullLogs{i},1)];   % Attach cumulative sums
    sumPelec(i) = fullLogs{i}(end,5);
    absDiscomfort(i) = fullLogs{i}(end,6)/wtp(i);
end
figure, hold all
plot(sumPelec(1:end-4),absDiscomfort(1:end-4),'b.-')
plot(sumPelec(end-3),absDiscomfort(end-3),'r*')
plot(sumPelec(end-2),absDiscomfort(end-2),'g*')
plot(sumPelec(end-1),absDiscomfort(end-1),'k*')
%plot(sumPelec(end),absDiscomfort(end),'m*')

% Linear regression stuff:
x = sumPelec(1:end-4)'; y = absDiscomfort(1:end-4)';
b1 = x\y;
yCalc1 = b1*x;
X = [ones(length(x),1) x];
b0 = X\y;
yCalc2 = X*b0;
plot(x,yCalc2,'--')
Rsq1 = 1 - sum((y - yCalc1).^2)/sum((y - mean(y)).^2);
Rsq2 = 1 - sum((y - yCalc2).^2)/sum((y - mean(y)).^2);

boundedInter = {[9500,sumPelec(end-3),sumPelec(end-3)],[absDiscomfort(end-3),absDiscomfort(end-3),0]};
unboundedInter = {[9500,sumPelec(end-2),sumPelec(end-2)],[absDiscomfort(end-2),absDiscomfort(end-2),0]};
baseInter = {[9500,sumPelec(end-1),sumPelec(end-1)],[absDiscomfort(end-1),absDiscomfort(end-1),0]};
idealInter = {[9500,sumPelec(end),sumPelec(end)],[absDiscomfort(end),absDiscomfort(end),0]};

c1text = strcat('w = ', c1Cell(1:end-4));
colormap = {'4c72b0','c44e52','8172b2','ccb974','64b5cd','cc64cd','cd8464','55a868'};
trace1 = struct('x',sumPelec(1:end-4),'y',absDiscomfort(1:end-4),'xaxis','x','mode', 'lines+markers+text','name','Dual-objective', ...
                'text', {c1text}, 'textposition','top right');
trace1.marker.size = 12;
trace2 = struct('x',sumPelec(end-3:end),'y',absDiscomfort(end-3:end),'xaxis','x','mode', 'markers+text', ...
                'name','Bounded w/ P(termination)', 'text',{{'Bounded w/ P(termination)','Bounded w/o P(termination)','Occupied 75F','Occupied Tideal=72.5F'}}, ...
                'textposition', {{'top right','top right','bottom left','top right'}});
trace2.marker.color = { colormap{7},colormap{end},colormap{2},colormap{3} };
trace2.marker.size = [12, 12, 12];

trace3 = struct('x',boundedInter{1},'y',boundedInter{2},'xaxis','x','type','line');
trace3.line.color = colormap{7};
trace3.line.dash = 'longdash';
trace3.marker.size = 0;

trace4 = struct('x',unboundedInter{1},'y',unboundedInter{2},'xaxis','x','type','line');
trace4.line.color = colormap{end};
trace4.line.dash = 'longdash';
trace4.marker.size = 0;

trace5 = struct('x',baseInter{1},'y',baseInter{2},'xaxis','x','type','line');
trace5.line.color = colormap{2};
trace5.line.dash = 'longdash';
trace5.marker.size = 0;

trace6 = struct('x',idealInter{1},'y',idealInter{2},'xaxis','x','type','line');
trace6.line.color = colormap{3};
trace6.line.dash = 'longdash';
trace6.marker.size = 0;

% 10e-3: using termination horizon, corbin
% 20e-3: no termination horizon, corbin
% 30e-3: dumb 75F
% 40e-3: dumb Tideal = 72.5F
fig3 = figure('units','normalized','outerposition',[0 0 1 1]);
xaxis = struct('domain', [0, 0.9], ...
      'showgrid', true, ...
      'zeroline', true, ...
      'showline', true, ...
      'gridcolor', '#bdbdbd', ...
      'gridwidth', 2, ...
      'zerolinecolor', '#969696', ...
      'zerolinewidth', 4, ...
      'linecolor', '#636363', ...
      'linewidth', 6, ...
      'nticks',14, ...
      'name', 'Summed Pelec', ...
      'range',[9.55e3,10.05e3]); 
yaxis =  struct('domain', [0, 0.9],'nticks',15, 'linecolor', '#636363','linewidth', 6, ...
                'gridcolor', '#bdbdbd','zeroline', true,'title','Summed absolute Discomfort [(delta F)^2*person*timestep]', ...
                'range',[200e3,420e3]); % and let the first y-axis span the first 45% of the plot height
layout = struct('xaxis',xaxis, 'yaxis',yaxis);
layout.width = 1200; % required
layout.height = 1000; % required

p = plotlyfig; close all
p.data = {trace1,trace2,trace3,trace4,trace5,trace6};
p.layout = layout;
p.PlotOptions.FileName = 'Plot_3';
html_file = plotlyoffline(p);
%r = plotly(p.data, struct('layout', layout));

disp(['View file: ', html_file])