behavioralSensitivity.m

% Behavioral Adjacency for all houses
% 8 April 2025
% Updated: 16 Sep 2025
% Dr. Kendall Parker
% The Behavioral Adjacency Model considers datasets adjacent if they
% differ only in the internal behavior of a single household over time,
% such as changes in consumption patterns. --> per-time-step changes in a
% household’s consumption.
% Sensitivity is therefore calibrated using a cycle based DCT assessment
% to capture the max difference at each time step of the most dominant cycle.
clc, clear all, close all

%% Import Data
load('dataOneYear93homes.mat')

%% Check if user wants to normalize the data
% Loop until the user enters a valid response
while true
    % Ask the user if they want to normalize the data
    normalizeData = input('Do you want to normalize the data? (y/n): ', 's');

    % Convert response to lowercase and remove spaces
    normalizeData = lower(strtrim(normalizeData));

    % Check user response
    if strcmp(normalizeData, 'y')
        % Normalize each variable using the function
        allConsump = f_normalize(allConsump);
        disp('Data has been normalized.');
        break; % Exit loop
    elseif strcmp(normalizeData, 'n')
        disp('Proceeding without normalization.');
        break; % Exit loop
    else
        % Display error message and ask again
        disp('Error: Please enter "y" or "n".');
    end
end

%% Preliminaries & Sanity check
[T, N] = size(allConsump); % T is number of samples, N is number of homes
samplingRateMin = 15; % [min]
samplingFreq = samplingRateMin/60; % [samples per hour]
fontSize = 12;
sampleUnitIdx = 10;
preservationAmt = 0.9; % For variance based thresholding
allTime = (0:1:T-1)' * 15 / 1440;

figure
plot(allTime, allConsump(:,1:10))
title('Normalized Energy Consumption - 10 Sample Homes')
xlabel('Time (days)','FontSize',fontSize)
ylabel('Consumption (kW)','FontSize',fontSize)
xticks(0:30:365)
xlim([0, 365])
legend()

%% Apply a Hamming Window
% For DCT, which assumes even symmetry, the 'symmetric' option ensures the 
% window tapers symmetrically about the midpoint — ideal for minimizing 
% edge effects.
w = hamming(T, 'symmetric'); 
allConsump = allConsump .* w;

%% For Loop to plot some of the homes and look at the cycles
numHomes = 10;
numDomCycles = 30;
fig = figure;
dctCoeffs = nan*ones(T,N);
ax = gobjects(numHomes,1);  % preallocate axes handles

for i = 1:numHomes
    % Compute DCT for household i
    dctCoeffs(:,i) = dct(allConsump(:,i));
    mags = abs(dctCoeffs(:,i));  % Magnitude of DCT coefficients

    % Find dominant coeff excluding DC component
    [sortedMags, sortIdx] = sort(mags(2:end), 'descend');
    sortIdx = sortIdx + 1; % Shift indices since we dropped the first one

    samplePeriodHours = samplingRateMin / 60; % 15 minutes → 0.25 hours
    
    % DCT freq index k corresponds to ~k/(2T) cycles per sample -> 2T/k samples per cycle
    cycleLengths = (2*T ./ sortIdx) * samplePeriodHours; % [hours]

    % if i == 1
    %     figure
    %     stem(cycleLengths(1:numDomCycles), sortedMags(1:numDomCycles), 'filled');
    %     xlabel('Cycle Length (hours)','FontSize',fontSize);
    %     ylabel('Magnitude','FontSize',fontSize);
    %     title(['(A) Household ' num2str(i) ' - Top Cycles']);
    %     set(gca,'XScale','log'); % optional: log scale helps spread short & long cycles
    %     set(gca, 'FontSize', fontSize);
    %     xlim([1, 5000]); % force start at 10^0 = 1
    %     hold on
    %     xline(4, '--k', '4h');
    %     xline(6, '--k', '6h');
    %     xline(8, '--r', '8h');
    %     xline(12, '--k', '12h');
    %     xline(24, '--k', '1d');
    %     xline(7*24, '--k', '7d');
    %     xline(14*24, '--k', '14d');
    %     xline(30*24,'--k','1m')
    %     xline(4*30*24,'--k','4m')
    % end

    % Create subplot and store handle
    figure
    ax(i) = subplot(numHomes, 1, i);

    % Plot results together
    % subplot(numHomes, 1, i);
    stem(cycleLengths(1:numDomCycles), sortedMags(1:numDomCycles), 'filled');
    set(gca,'XScale','log');
    xlim([1, 5000]);
    title(['Household ' num2str(i)],'FontSize',fontSize-3);
    set(gca,'FontSize',fontSize-4);

end
han = axes(fig,'visible','off');
han.XLabel.Visible = 'on';
han.YLabel.Visible = 'on';
xlabel(han,'Cycle Length (hours)','FontSize',fontSize);
ylabel(han,'DCT Magnitude','FontSize',fontSize);
title(han,'Top Cycles','FontSize',fontSize)

% --- Add vertical reference lines to every subplot ---
xlineVals = [1/6, 1/4, 1/3, 0.5, 1, 7, 14, 30, 4*30].*24;  % cycles to mark (hours)
xlineLabels = {'4h','6h','8h','12h','1d','7d','14d','30d','4m'};
for i = 1:numel(ax)
    for j = 1:numel(xlineVals)
        xl = xlineVals(j);
        if i == 1
            % add labeled line only on first subplot
            xline(ax(i), xl, '--r', xlineLabels{j}, ...
                'LabelVerticalAlignment', 'bottom', ...
                'LabelHorizontalAlignment', 'left');
        else
            % unlabeled lines on others
            xline(ax(i), xl, '--r');
        end
    end
end


% figure
% stem(cycleLengths(1:numDomCycles), sortedMags(1:numDomCycles), 'filled');
% xlabel('Cycle Length (hours)','FontSize',fontSize);
% ylabel('Magnitude','FontSize',fontSize);
% title(['(A) Household ' num2str(10) ' - Top Cycles']);
% set(gca,'XScale','log'); % optional: log scale helps spread short & long cycles
% set(gca, 'FontSize', fontSize);
% xlim([1, 5000]); % force start at 10^0 = 1
% hold on
% xline(12, '--r', '12 hours');
% xline(24, '--r', '1 day');
% xline(7*24, '--r', '1 week');
% xline(14*24, '--r', '2 weeks');
% xline(30*24,'--r','1 month')
% xline(3*30*24,'--r','3 months')

%% For loop to cycle through the houses and find sensitivity
maxNumHomes = 20;
cycleLengthsMin = xlineVals*60; % [min]
cycleLabels = xlineLabels;

numOfCycles = length(cycleLengthsMin);
allSensitivities = nan*ones(numOfCycles,maxNumHomes);

for ii = 1:maxNumHomes
    for jj = 1:numOfCycles % Select a cycle length
        
        % Convert cycle lengths to number of samples per cycle
        cycleLengthsSamples = cycleLengthsMin(jj) / samplingRateMin; %[samples]

        % Compute how many full cycles fit in T
        numCycles = floor(T / cycleLengthsSamples);

        % Truncate the vector so its length is divisible by cycleLength
        truncateVector = allConsump(1:cycleLengthsSamples*numCycles, ii);

        % Reshape: each column = one cycle
        houseConsump = reshape(truncateVector, cycleLengthsSamples, numCycles);

        % Initialize a vector to store the max difference per row
        timeMaxDiff = zeros(cycleLengthsSamples, 1);

        % Loop through each row
        for row = 1:cycleLengthsSamples
            maxDiff = 0; % Reset max difference for this row

            % Loop through unique column pairs
            for i = 1:size(houseConsump,2)
                for j = i+1:size(houseConsump,2)
                    % Compute absolute difference between columns i and j
                    colDiff = abs(houseConsump(row, i) - houseConsump(row, j));

                    % Update maxDiff for this row
                    maxDiff = max(maxDiff, colDiff);
                end
            end

            % Store the max difference for this row
            timeMaxDiff(row) = maxDiff;
        end
        % Take DCT of Perturbation
        housePertDct = dct(timeMaxDiff);
        allSensitivities(jj,ii) = norm(housePertDct(2:end),2);
    end
end

%% Sensitivity Heatmap for Appendix
heatMapPart1 = allSensitivities(1:round(numOfCycles/2),:);
figure;
imagesc(heatMapPart1); % rows = cycles, columns = households
colormap(parula); % other options 'hot', 'turbo' 
c = colorbar;
clim([min(heatMapPart1(:)), max(heatMapPart1(:))]); % for consistent scaling

% Set colorbar ticks and labels to reflect values in allSensitivities
c.Ticks = linspace(min(heatMapPart1(:)), max(heatMapPart1(:)), numOfCycles);
c.TickLabels = arrayfun(@(v) sprintf('%.2f', v),...
    linspace(min(heatMapPart1(:)), max(heatMapPart1(:)),...
    numOfCycles), 'UniformOutput', false);

% Axis formatting
xlabel('Household Index','FontSize',fontSize);
ylabel('Cycle Length','FontSize',fontSize);
title('Behavioral Sensitivity by Household and Cycle Length','FontSize',fontSize+1);

% Custom tick labels
yticks(1:length(cycleLabels));
yticklabels(cycleLabels(:,1:round(numOfCycles/2)));
xticks(1:maxNumHomes);
xticklabels(string(1:maxNumHomes));
set(gca, 'FontSize', fontSize);

% SECOND HEATMAP
heatMapPart2 = allSensitivities(numOfCycles/2+1:end,:);
figure;
imagesc(heatMapPart2); % rows = cycles, columns = households
colormap(parula); % You can try 'hot', 'viridis' (with add-on), or 'turbo' too
c = colorbar;
clim([min(heatMapPart2(:)), max(heatMapPart2(:))]); % for consistent scaling

% Set colorbar ticks and labels to reflect values in allSensitivities
c.Ticks = linspace(min(heatMapPart2(:)), max(heatMapPart2(:)), numOfCycles);
c.TickLabels = arrayfun(@(v) sprintf('%.2f', v),...
    linspace(min(heatMapPart2(:)), max(heatMapPart2(:)),...
    numOfCycles), 'UniformOutput', false);

% Axis formatting
xlabel('Household Index','FontSize',fontSize);
ylabel('Cycle Length','FontSize',fontSize);
title('Behavioral Sensitivity by Household and Cycle Length','FontSize',fontSize+1);

% Custom tick labels
yticks(1:length(cycleLabels));
yticklabels(cycleLabels(:,round(numOfCycles/2+1):end));
xticks(1:maxNumHomes);
xticklabels(string(1:maxNumHomes));
set(gca, 'FontSize', fontSize);

%% Select a house and plot cycle in time domain
sampleHouse = 1;
selectedCycleLength = 8*60; % [min]
cycleStr = '8h';

% Compute how many full cycles fit in T
numCycles = floor(T / selectedCycleLength);

% Truncate the vector so its length is divisible by cycleLength
truncateVector = allConsump(1:selectedCycleLength*numCycles, sampleHouse);

% Reshape: each column = one cycle
sampleHouseReshape = reshape(truncateVector, selectedCycleLength, numCycles);

timeAxis = linspace(0,selectedCycleLength/60,selectedCycleLength); %(0:samplingFreq:selectedCycleLength/60-1); % in hours

figure
plot(timeAxis, sampleHouseReshape)
xlabel('Time Within Cycle (hours)','FontSize',fontSize)
ylabel('Consumption (kW)','FontSize',fontSize)
title(['(B) House ' num2str(sampleHouse) ' Cycle Overlay: ' cycleStr],'FontSize',fontSize)
set(gca, 'FontSize', fontSize);
grid on

%% Visualize cycles for selected house
% Initialize a vector to store the max difference per row
timeMaxDiff = zeros(selectedCycleLength, 1);

% Loop through each row
for row = 1:selectedCycleLength
    maxDiff = 0; % Reset max difference for this row

    % Loop through unique column pairs
    for i = 1:size(sampleHouseReshape,2)
        for j = i+1:size(sampleHouseReshape,2)
            % Compute absolute difference between columns i and j
            colDiff = abs(sampleHouseReshape(row, i) - sampleHouseReshape(row, j));

            % Update maxDiff for this row
            maxDiff = max(maxDiff, colDiff);
        end
    end

    % Store the max difference for this row
    timeMaxDiff(row) = maxDiff;
end

timeVector = (0:selectedCycleLength-1)' * 15 / 1440;
% ((0:672-1)' * 15) / 1440;  % days

figure
plot(timeVector,timeMaxDiff,'.-','MarkerSize',10)
xlabel('Time (days)','FontSize',fontSize)
ylabel('Difference in Consumption at Time Stamp','FontSize',fontSize)
title(['(C) House ' num2str(sampleHouse) ' Perturbations for ' cycleStr ' Cycles'])
set(gca, 'FontSize', fontSize);
grid on; box on;

% Take DCT of Perturbation
% checkTimeNorm = norm(timeMaxDiff,2);
[singleHousePertDct, ~] = f_varianceThreshold(dct(timeMaxDiff),preservationAmt);

figure
plot(singleHousePertDct,'b-','LineWidth', 1.5)
title(['(D) DCT of Perturbations - House ' num2str(sampleHouse)],...
    'FontSize', fontSize+2, 'FontWeight', 'bold');
xlabel('Index', 'FontSize', fontSize);
ylabel('Magnitude', 'FontSize', fontSize);
set(gca, 'FontSize', fontSize);
grid on; box on;
% ax = gca; ax.FontSize = fontSize-2;
% ax.XMinorGrid = 'on'; ax.YMinorGrid = 'on';

householdDctS2 = norm(singleHousePertDct(2:end),2)