%% Data processing of 50 subjects included for feasibility of fever assessment by wearable temperature sensor

%% Color map definition
Sm = -4:0.01:4;
smarrR = [(0.5:0.005:1) (1:-0.002:0.004) (0:0.005:1)];%{0.5+sin(Sm./1-3)./2;%(801:-1:1)./801;%(1.*Sm.^2+0.5.*Sm+0)./20+0.005; %0.5+sin(Sm./1.5-1.5)./2;
smarrG = [zeros(1,550) (0:0.004:1)]; %(3.^Sm)./81;%1./(1+2.^-Sm);
smarrB = [zeros(1,100) (0:0.0016:1) ones(1,75)];%(1.4.^Sm)./4;%(1:1:801)./801;
smarr = [smarrB' smarrG' smarrR'];

SickC = [0 0 0.5; 0.9 0 0; 0 0.4 1; 0.3 0.7 0]; %consistent colors for all four modalities

%% loading subject data
load('PubData.mat')
%loaded as 45 days prior to symptom report + 20 days after (including day 1); Symptom report onset is the 46th day.
% Symptom duration is in cell column 2
% Days prior successfully loaded is in cell column 3. this is for averaging

%% Cleaning and organizing matrices for aligned comparison by variable.
Ln = length(Cases50);
% Variable naming convention: T for temperature; H for HR, V for HRV, and B for RR (Breathing rate)
AggT = zeros(93600,Ln); %Agg for aggregated data
for i = 1:Ln
    a = Cases50{i,1}(:,1); % load the data for case i
    if Cases50{i,3} < 45 % checks if buffer is needed
        b = zeros((45-Cases50{i,3}).*1440,1); %if so, adds buffer to fill days specificied in Cases50{i,3};
        AggT(:,i) = [b' a']';
    else
        AggT(:,i) = a;
    end
end

AggH = zeros(93600,Ln);
for i = 1:Ln
    a = Cases50{i,1}(:,2); % load the data for case i
    if Cases50{i,3} < 45 % checks if buffer is needed
        b = zeros((45-Cases50{i,3}).*1440,1);
        AggH(:,i) = [b' a']';
    else
        AggH(:,i) = a;
    end
end 

AggV = zeros(93600,Ln);
for i = 1:Ln
    a = Cases50{i,1}(:,3); % load the data for case i
    if Cases50{i,3} < 45 % checks if buffer is needed
        b = zeros((45-Cases50{i,3}).*1440,1);
        AggV(:,i) = [b' a']';
    else
        AggV(:,i) = a;
    end
end 

AggB = zeros(93600,Ln);
for i = 1:Ln
    a = Cases50{i,1}(:,4); % load the data for case i
    if Cases50{i,3} < 45 % checks if buffer is needed
        b = zeros((45-Cases50{i,3}).*1440,1);
        AggB(:,i) = [b' a']';
    else
        AggB(:,i) = a;
    end
end 
clear a b

AggBIDfix = [6 85221; 9 82045; 18 75170; 21 89366; 22 83756; 27 82143; 29 88068; 36 87222; 44 75000; 47 80826; 49 91316];
for i = 1:length(AggBIDfix)
    AggB(AggBIDfix(i,2):end,AggBIDfix(i,1)) = NaN;
end
AggB(1:59117,28) = NaN;
AggT(AggT==0) = NaN;
AggH(AggH==0) = NaN;
AggV(AggV==0) = NaN;
AggB(AggB==0) = NaN;

%% parameters for comparisons
for i = 1:Ln
    SymDur(i,1) = Cases50{i,2}; %extracting subject symptom window durations
end
SymOn = 1440*45.3; SymOff = SymOn+1440*mean(SymDur); %defining onset and offset for individual symptom windows
SymD20 = SymDur; SymD20(SymDur>20)=20; %limits displayed symptom window to the end of the data

for i = 1:Ln
    AggTz(:,i) = (AggT(:,i)-nanmean(AggT(1:1440*40,i)))./nanstd(AggT(1:1440*40,i));
    AggHz(:,i) = (AggH(:,i)-nanmean(AggH(1:1440*40,i)))./nanstd(AggH(1:1440*40,i));
    AggVz(:,i) = (AggV(:,i)-nanmean(AggV(1:1440*40,i)))./nanstd(AggV(1:1440*40,i));
    AggBz(:,i) = (AggB(:,i)-nanmean(AggB(1:1440*40,i)))./nanstd(AggB(1:1440*40,i));
end
AggTz0 = AggTz; AggHz0 = AggHz; AggVz0 = AggVz; AggBz0 = AggBz;
AggTz0(isnan(AggTz0)) = 0; AggHz0(isnan(AggHz0)) = 0; AggVz0(isnan(AggVz0)) = 0; AggBz0(isnan(AggBz0)) = 0;

%% Individual ranges - Figure 1A
for i = 1:Ln
CaseBMs(1,i) = nanmean(AggT(1:1440*40,i)); %BMs = Baseline Means;
CaseBSs(1,i) = nanstd(AggT(1:1440*40,i)); %BSs = Baseline Standard deviations
CaseSMs(1,i) = nanmean(AggT(44*1440:1440*(44+SymD20(i)),i)); %SMs = Symptom window Means
CaseSSs(1,i) = nanstd(AggT(44*1440:1440*(44+SymD20(i)),i)); %SSs = Symptom window Standard deviations
end
[Ts Tsi] = sort(CaseBMs(1,:)); %sorting individuals for Figure 1A.
[TSs TSsi] = sort(CaseSMs(1,:));

p = ranksum(CaseBMs(1,:),CaseSMs(1,:))% testing change in mean T from baseline to symptom window, within individual

pinetrees = figure;
hold on; errorbar([2:1:51],Ts,CaseBSs(1,Tsi),'LineStyle','none','Color',SickC(1,:)); scatter([2:1:51],Ts,30,SickC(1,:),'filled'); set(gca,'TickDir','out'); 
hold on; errorbar([1.9:1:50.9],Ts,CaseSSs(1,TSsi),'LineStyle','none','Color',SickC(1,:)./3); scatter([1.9:1:50.9],TSs,30,SickC(1,:)./3); set(gca,'TickDir','out'); 

%% Generating daily min and max for Figure 1B-D
% ID:  Subject ID (internal to this analysis only - not connected to other data) for figure labels
ID   =  [01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50];
% Fever: Has fever seld report = 01; does not = 02 (00 caused errors in some places, which is lame)
Fever = [02 01 01 01 01 01 01 02 02 01 01 01 02 01 02 01 02 01 01 02 01 01 01 01 01 01 01 01 01 02 01 01 02 02 01 01 01 01 01 02 01 01 01 01 01 02 01 01 01 01]; % did subject report fever as a symptom
% FevBot: fever by robot: Has fever-like day detected by digital biomarkers = 01; does not = 02
FevBot =[02 01 01 01 01 02 01 02 02 01 01 01 02 01 01 01 02 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 02 01 01 01 01 01 01 02 01 01 01 01 01 01]; % rule: max > 1.2 &| min > -0.2 in normalized zscores, excluding artifacts
% PrSymF: presymptom fever detected = 01; not = 02;
PrSymF =[01 01 01 01 01 02 02 01 01 01 01 02 02 01 01 01 02 01 01 01 02 02 01 02 01 01 01 01 01 01 01 01 01 01 01 02 02 01 02 01 01 01 01 02 01 01 01 01 01 01]; % was there sickness (by rule) before symptom report onset

for i = 1:Ln
    for j = 1:65
        A = sort(movmean(AggTz((j-1)*1440+1:1440*j,i),360));
        AggMinTz(j,i) = median(A(1:360)); %generates min and max from Zscore, but Zscores are messed up still, so then uses these min and maxes to renormalize
        AggMaxTz(j,i) = median(A(1320:1440));
        D = sort(movmean(AggHz((j-1)*1440+1:1440*j,i),360));
        AggMinHz(j,i) = median(D(1:360)); %generates min and max from Zscore, but Zscores are messed up still, so then uses these min and maxes to renormalize
        AggMaxHz(j,i) = median(D(1320:1440));
        G = sort(movmean(AggVz((j-1)*1440+1:1440*j,i),360));
        AggMinVz(j,i) = median(G(1:360)); %generates min and max from Zscore, but Zscores are messed up still, so then uses these min and maxes to renormalize
        AggMaxVz(j,i) = median(G(1320:1440));
        J = sort(movmean(AggBz((j-1)*1440+1:1440*j,i),360));
        AggMinBz(j,i) = median(J(1:360)); %generates min and max from Zscore, but Zscores are messed up still, so then uses these min and maxes to renormalize
        AggMaxBz(j,i) = median(J(1320:1440));
    clear A D G J
    end
    A = nanmedian(AggMinTz(:,i)); B = nanmedian(AggMaxTz(:,i));
    AggMTzN(:,i) = (AggTz(:,i)-A)./(B-A).*2-1; % normalized from -1 to 1 by Zscore min and maxes
    D = nanmedian(AggMinHz(:,i)); E = nanmedian(AggMaxHz(:,i));
    AggMHzN(:,i) = (AggHz(:,i)-D)./(E-D).*2-1; % normalized from -1 to 1 by Zscore min and maxes
    G = nanmedian(AggMinVz(:,i)); H = nanmedian(AggMaxVz(:,i));
    AggMVzN(:,i) = (AggVz(:,i)-G)./(H-G).*2-1; % normalized from -1 to 1 by Zscore min and maxes
    J = nanmedian(AggMinBz(:,i)); K = nanmedian(AggMaxBz(:,i));
    AggMBzN(:,i) = (AggBz(:,i)-J)./(K-J).*2-1; % normalized from -1 to 1 by Zscore min and maxes
    clear A B D E G H J K
    for k = 1:65
        A = sort(movmean(AggMTzN((k-1)*1440+1:1440*k,i),360));
        AggMinTn(k,i) = median(A(1:360)); %regenerating the min and maxes from the normalized Zscore data
        AggMaxTn(k,i) = median(A(1320:1440));
        D = sort(movmean(AggMHzN((k-1)*1440+1:1440*k,i),360));
        AggMinHn(k,i) = median(D(1:360)); %regenerating the min and maxes from the normalized Zscore data
        AggMaxHn(k,i) = median(D(1320:1440));
        G = sort(movmean(AggMVzN((k-1)*1440+1:1440*k,i),360));
        AggMinVn(k,i) = median(G(1:360)); %regenerating the min and maxes from the normalized Zscore data
        AggMaxVn(k,i) = median(G(1320:1440));
        J = sort(movmean(AggMBzN((k-1)*1440+1:1440*k,i),360));
        AggMinBn(k,i) = median(J(1:360)); %regenerating the min and maxes from the normalized Zscore data
        AggMaxBn(k,i) = median(J(1320:1440));
    end
end
% sorting out those with and without fever as either self-report or detected by digital biomarker
j = 0; k = 0; l=0; m=0;
for i = 1:Ln
    if Fever(1,i) == 1
        j = j+1;
        Fever1s(j) = i;
    else
        k = k+1;
    Fever0s(k) = i; %replaced all 0s with 2s in Fever because it couldn't find 0s. Why.
    end
    if FevBot(1,i) == 1
        l = l+1;
        FevBo1s(l) = i;
    else
        m = m+1;
    FevBo0s(m) = i; %replaced all 0s with 2s in Fever because it couldn't find 0s.
    end    
end
F1SymDur = SymDur(Fever1s); FB1SymDur = SymDur(FevBo1s);
F0SymDur = SymDur(Fever0s); FB0SymDur = SymDur(FevBo0s);
F1SymOn = 1440*45.3; F1SymOff = F1SymOn+1440*mean(F1SymDur); FB1SymOff = F1SymOn+1440*mean(FB1SymDur);
F0SymOn = 1440*45.3; F0SymOff = F0SymOn+1440*mean(F0SymDur); FB0SymOff = F1SymOn+1440*mean(FB0SymDur);

fevers = figure; %Figure 1B-D

subplot(2,2,[1 3]);
hold on
plot(nanmedian(AggMTzN(:,FevBo0s),2),'Color',SickC(1,:),'LineWidth',2); xlim([0 1440*65]);
plot(nanmedian(AggMTzN(:,FevBo0s),2),'w','LineWidth',1); xlim([0 1440*65]);
plot(nanmedian(AggMTzN(:,FevBo1s),2),'Color',SickC(1,:),'LineWidth',3); 
xtick([F0SymOn FB0SymOff FB1SymOff]); xticklabels({'Symptoms On', '1s Off', '0s Off'}); ylim([-1.5 1.2]); set(gca,'TickDir','out'); 

subplot(2,2,2);
hold on; ylabel('Temp')
%fill([F1SymOn FB1SymOff FB1SymOff F1SymOn],[-3 -3 3 3],[0.95 0.95 0.95],'LineStyle','none');
scatter([1:1:65].*1440,nanmedian(AggMaxTn(:,FevBo0s),2),40,SickC(1,:)); scatter([1:1:65].*1440,nanmedian(AggMinTn(:,Fever1s),2),36,SickC(1,:));
scatter([1:1:65].*1440,nanmedian(AggMaxTn(:,FevBo1s),2),36,SickC(1,:),'filled');
xtick([F0SymOn FB0SymOff FB1SymOff]); xticklabels({'Symptoms On', '1s Off', '0s Off'}); ylim([0.9 1.2]); set(gca,'TickDir','out'); 
plot([1:1:65].*1440,movmean(nanmedian(AggMaxTn(:,FevBo0s),2),3),'Color',SickC(1,:),'LineWidth',3);
plot([1:1:65].*1440,movmean(nanmedian(AggMaxTn(:,FevBo0s),2),3),'w','LineWidth',2)
plot([1:1:65].*1440,movmean(nanmedian(AggMaxTn(:,FevBo1s),2),3),'Color',SickC(1,:),'LineWidth',3);
xlim([0 1440*65]);

[p,h] = ranksum((nanmedian(AggMaxTn(1:40,FevBo0s),2)),(nanmedian(AggMaxTn(45:52,FevBo0s),2)),'tail','left'); % p = 0.4
[p,h] = ranksum((nanmedian(AggMaxTn(1:40,FevBo1s),2)),(nanmedian(AggMaxTn(45:52,FevBo1s),2)),'tail','left'); % p = 7.5E-6

subplot(2,2,4);
hold on; ylabel('Temp')
%fill([F1SymOn FB1SymOff FB1SymOff F1SymOn],[-3 -3 3 3],[0.95 0.95 0.95],'LineStyle','none');
scatter([1:1:65].*1440,nanmedian(AggMinTn(:,FevBo0s),2),40,SickC(1,:)); scatter([1:1:65].*1440,nanmedian(AggMinTn(:,Fever1s),2),36,SickC(1,:));
scatter([1:1:65].*1440,nanmedian(AggMinTn(:,FevBo1s),2),36,SickC(1,:),'filled');
xtick([F0SymOn FB0SymOff FB1SymOff]); xticklabels({'Symptoms On', '1s Off', '0s Off'}); ylim([-1.5 -0.5]); set(gca,'TickDir','out'); 
plot([1:1:65].*1440,movmean(nanmedian(AggMinTn(:,FevBo0s),2),3),'Color',SickC(1,:),'LineWidth',3);
plot([1:1:65].*1440,movmean(nanmedian(AggMinTn(:,FevBo0s),2),3),'w','LineWidth',2)
plot([1:1:65].*1440,movmean(nanmedian(AggMinTn(:,FevBo1s),2),3),'Color',SickC(1,:),'LineWidth',3);
xlim([0 1440*65]);

[p,h] = ranksum((nanmedian(AggMinTn(1:40,FevBo0s),2)),(nanmedian(AggMinTn(45:52,FevBo0s),2)),'tail','left'); % p = 0.4
[p,h] = ranksum((nanmedian(AggMinTn(1:40,FevBo1s),2)),(nanmedian(AggMinTn(45:52,FevBo1s),2)),'tail','left'); % p = 2E-5

%block for adding symptom windows as grey rectangles in different plots
Sf1 = (length(FevBo1s).^0.5); Sf0 = (length(FevBo0s).^0.5); %to be added by figure - change FevBo1s as appropriate based on group being plotted
hold on; ylabel('Temp')
fill([F1SymOn FB1SymOff FB1SymOff F1SymOn],[-3 -3 3 3],[0.9 0.9 0.9],'LineStyle','none');
fill([1:1:1440*65 1440*65:-1:1],[(nanmedian(AggMTzN(:,FevBo1s),2)+nanstd(AggMTzN(:,FevBo1s),[],2)./Sf1)' flipud(nanmedian(AggMTzN(:,FevBo1s),2)-nanstd(AggMTzN(:,FevBo1s),[],2)./Sf1)'],SickC(1,:),'LineStyle','none'); alpha(0.3) 
plot(nanmedian(AggMTzN(:,FevBo1s),2),'Color',SickC(1,:))
xtick([F1SymOn FB1SymOff]); xticklabels({'Symptoms On', 'Symptoms Off'}); ylim([-2 1.5]); set(gca,'TickDir','out'); 
plot([0 65*1440],[1.1 1.1],'k','LineStyle',':'); plot([0 65*1440],[-0.4 -0.4],'k','LineStyle',':')
scatter([1:1:65].*1440,nanmedian(AggMinTn(:,FevBo1s),2),40,'k'); scatter([1:1:65].*1440,nanmedian(AggMinTn(:,FevBo1s),2),36,[0.1 0.5 1],'filled');
scatter([1:1:65].*1440-480,nanmedian(AggMaxTn(:,FevBo1s),2),40,'k'); scatter([1:1:65].*1440-480,nanmedian(AggMaxTn(:,FevBo1s),2),36,[0 0.3 0.8],'filled');

% Figure 1 E-H are windows on ID 3 from Agg* of *= T or H or V or B

%% processing non-interpolated data (mostly at night, hence the name) for Figure 2
for i = 1:Ln
   NightSums(:,i) = sum(Cases50{i,4}(:,:),2);
     a = find(Cases50{i,4}(:,2) > 0);
    NightsT(a,i) = AggT(a,i); % Temp  
    NightsH(a,i) = AggV(a,i); % HR 
    NightsV(a,i) = AggH(a,i); % HRv 
    NightsB(a,i) = AggB(a,i); % BR 
    NightsTime(a,i) = a; % Time stamps
  clear a NightSums
end

NightsT(NightsT==0) = NaN;
NightsH(NightsH==0) = NaN;
NightsV(NightsV==0) = NaN;
NightsB(NightsB==0) = NaN;

for i = 1:Ln
NightsTz(:,i) = (NightsT(:,i)-nanmean(NightsT(1:40*1440,i)))./nanstd(NightsT(1:40*1440,i));
NightsHz(:,i) = (NightsH(:,i)-nanmean(NightsH(1:40*1440,i)))./nanstd(NightsH(1:40*1440,i));
NightsVz(:,i) = (NightsV(:,i)-nanmean(NightsV(1:40*1440,i)))./nanstd(NightsV(1:40*1440,i));
NightsBz(:,i) = (NightsB(:,i)-nanmean(NightsB(1:40*1440,i)))./nanstd(NightsB(1:40*1440,i));
end

for i = 1:Ln
NightsC{i,1} = NightsTz((all((~isnan(NightsTz(:,i))),2)),i);
NightsC{i,2} = NightsHz((all((~isnan(NightsHz(:,i))),2)),i);
NightsC{i,3} = NightsVz((all((~isnan(NightsVz(:,i))),2)),i);
NightsC{i,4} = NightsBz((all((~isnan(NightsBz(:,i))),2)),i);
end

% scatter made by sorting individuals by category, using FevBot: FevBo0 = negative, FevBo1 = positive
scatterfun = figure;
subplot(1,4,1:2); set(gca,'Color',[0.6 0.6 0.6]); axis([-6 12 -6 12]);
subplot(1,4,3:4); set(gca,'Color',[0.6 0.6 0.6]); axis([-6 12 -6 12]);
for i = 1:Ln
    if isempty(find(FevBo1s == i))
        j = find(FevBo0s == i);
            if FB0SymDur(j) > 19
                k = 19;
            else
                k = FB0SymDur(j);
            end
        subplot(1,4,3:4); hold on; scatter(NightsHz(F0SymOn-1440:F0SymOn+k*1440,i),NightsVz(F0SymOn-1440:F0SymOn+k*1440,i),(10.*NightsBz(F0SymOn-1440:F0SymOn+k*1440,i)).^2,...
            [0.2 0.2 0.2],'filled'); alpha(0.2)
    else
        j = find(FevBo1s == i);
            if FB1SymDur(j) > 19
                k = 19;
            else
                k = FB1SymDur(j);
            end
        subplot(1,4,1:2); hold on; scatter(NightsHz(F1SymOn-1440:F1SymOn+k*1440,i),NightsVz(F1SymOn-1440:F1SymOn+k*1440,i),(10.*NightsBz(F1SymOn-1440:F1SymOn+k*1440,i)).^2,...
            [0.2 0.2 0.2],'filled'); alpha(0.2);
    end
end
for i = 1:Ln
    if isempty(find(FevBo1s == i))
        j = find(FevBo0s == i);
            if FB0SymDur(j) > 19
                k = 19;
            else
                k = FB0SymDur(j);
            end
        subplot(1,4,1:2); hold on; scatter(NightsHz(1:40*1440,i),NightsVz(1:40*1440,i),(10.*NightsBz(1:40*1440,i)).^2,'w','filled'); alpha(0.2)
    else
        j = find(FevBo1s == i);
            if FB1SymDur(j) > 19
                k = 19;
            else
                k = FB1SymDur(j);
            end
        subplot(1,4,3:4); hold on; scatter(NightsHz(1:40*1440,i),NightsVz(1:40*1440,i),(10.*NightsBz(1:40*1440,i)).^2,'w','filled'); alpha(0.2)
    end
end
for i = 1:Ln
    if isempty(find(FevBo1s == i))
        j = find(FevBo0s == i);
            if FB0SymDur(j) > 19
                k = 19;
            else
                k = FB0SymDur(j);
            end
        ScatC = ((NightsTz(:,i)-min(NightsTz(:,i)))./(max(NightsTz(:,i))-min(NightsTz(:,i)))).*(max(AggMaxTn(44:45+k,i)./1.2)); ScatC(ScatC>1)=1; %normalizing to fever threshold 1.2
        subplot(1,4,1:2); hold on; scatter(NightsHz(F0SymOn-1440:F0SymOn+k*1440,i),NightsVz(F0SymOn-1440:F0SymOn+k*1440,i),(10.*NightsBz(F0SymOn-1440:F0SymOn+k*1440,i)).^2,...
            [ScatC(F0SymOn-1440:F0SymOn+k*1440).^10 abs(ScatC(F0SymOn-1440:F0SymOn+k*1440)-0.75).^1 1-(ScatC(F0SymOn-1440:F0SymOn+k*1440)).^10],'filled'); alpha(0.2)
    else
        j = find(FevBo1s == i);
            if FB1SymDur(j) > 19
                k = 19;
            else
                k = FB1SymDur(j);
            end
        ScatC = ((NightsTz(:,i)-min(NightsTz(:,i)))./(max(NightsTz(:,i))-min(NightsTz(:,i)))).*(max(AggMaxTn(44:45+k,i)./1.2)); ScatC(ScatC>1)=1;
        subplot(1,4,3:4); hold on; scatter(NightsHz(F1SymOn-1440:F1SymOn+k*1440,i),NightsVz(F1SymOn-1440:F1SymOn+k*1440,i),(10.*NightsBz(F1SymOn-1440:F1SymOn+k*1440,i)).^2,...
            [ScatC(F1SymOn-1440:F1SymOn+k*1440).^10 abs(ScatC(F0SymOn-1440:F0SymOn+k*1440)-0.75).^1 1-(ScatC(F1SymOn-1440:F1SymOn+k*1440)).^10],'filled'); alpha(0.2);
    end
end
set(gcf, 'InvertHardcopy', 'off');

%% Wavelets and peak determination for Figure 4 
% cirv v ult
% 24h: 336; 9h: 248; 5h: 194; 1h: 45.plot([336, 248, 194, 45],[24 9 5 1]);
x1 = [1 5 9 24]; y1 = [45 194 248 336];
%y = -5992*x^-0.01567+6037;
for i = 1:30; %Ys converts wave Y axis into hours (periods)
    Ys(i,1) = i; Ys(i,2) = round(-5992*i^-0.01567+6037);
end
%y scale scatter plot
%scatter(Ys(:,1),Ys(:,2))

%TEMPERATURE
%taking wavelet power of Temperature @ 1 / 1 min
NstepsPerHr=60; TimeStep=1/NstepsPerHr;
period1=24;  % period in hours of main oscillation
period2=24*56;% period in hours governing the varying frequency or amplitude
Nsteps=24*60*NstepsPerHr+1; t=(0:Nsteps-1)'/NstepsPerHr; T1=t(Nsteps)-t(1);
longestperiod=40;shortestperiod=1;  % range of periods to use in AWT
gamma=3;beta=5; % parameter values for Morse wavelet function
nvoice=64;
[fs,tau,qscaleArray] = CalcScaleForCWT(shortestperiod,longestperiod,T1,NstepsPerHr,nvoice);
taubnds=round(interp1(tau,1:numel(tau),[1 39]));

WtTn = AggMTzN;
WtTn(isnan(WtTn)) = 0;
for i = 1:Ln %downsampling to band max for circadian (C = 22-26h) and ultradian (U = 2-6 h)
     WavesTC(:,i) = max(WTTn(:,Ys(22,2):Ys(26,2),i),[],2);
end
%By hand, identifying peaks in reference to first day with Tmin above threshold of "fever-like day", and identifying artifacts and data holes
N = 22; %iterated by hand to ensure each one was not related to missing data or other artifacts for which we did not implement algorithmic cleaning
figy = figure;
subplot(3,1,1:2); hold on; plot(1:1:65, WTMaxTn(:,N),'b'); %plot(1:1:65, WTMinTn(:,N),'b');
plot(1:1:65, AggMinTn(:,N),'r'); %plot(1:1:65, AggMaxTn(:,N),'r');
xticks([0.3:7:42.3 45.3 49.3:7:70.3]); xlim([0 66]); yticks([-0.6 0 1.2]); grid on; xticklabels({[0:7:42],'SO',[49:7:70]});  title(N)
subplot(3,1,3); plot(AggTz(:,N)); %49 is crazy! It disintegrates into no daily rhythms!
xticks([0.3*1440:7*1440:42.3*1440 45.3*1440 49.3*1440:7*1440:70.3*1440]); yticks([-0.6 0 1.2]); grid on; xticklabels({[0:7:42],'SO',[49:7:70]}); xlim([0 66*1440]);
%IDs with primarily night data only; 10,16,20,30,35,36,39,40,43,47
wavepeak = [12,14,1; 22,26,1; 29,30,1; 33,37,1; 36,43,1; 21,25,2; 27,29,2; 33,36,2; 41,42,2; 45,50,2; 53,54,2; 9,13,3; 23,26,3; 34,38,3; 42,49,3; 53,56,3; 58,61,3; 42,46,4; 50,56,4;...
  8,10,5; 17,18,5; 24,26,5; 40,41,5; 42,43,5; 40,45,5; 1,7,6; 8,14,6; 13,17,6; 20,21,6; 36,42,6; 50,52,6; 54,58,6; 57,62,6; 5,10,7; 30,34,7; 39,40,7;...
   47,48,7; 56,61,7; 5,6,8; 12,18,8; 15,21,8; 20,23,8; 24,27,8; 28,35,8; 41,48,8; 51,53,8; 62,65,8; 21,25,9; 26,31,9; 44,48,9; 50,52,9; 59,61,9; 10,11,10;...
    8,13,11; 14,15,11; 21,26,11; 30,35,11; 35,42,11; 44,49,11; 48,52,11; 51,54,11; 55,56,11; 57,58,11; 42,44,12; 46,47,12; 53,57,12; 60,62,12; 63,64,12; 3,5,13; 20,24,13;...
     27,33,13; 35,38,13; 55,57,13; 4,5,14; 8,9,14; 11,12,14; 14,16,14; 18,22,14; 26,28,14; 37,41,14; 42,43,14; 46,49,14; 52,53,14; 55,56,14; 17,20,15; 46,47,15; 51,54,15; 5,7,17;...
      33,38,17; 44,46,17; 56,57,17; 60,62,17; 16,23,18; 25,30,18; 33,34,18; 36,38,18; 47,50,18; 57,58,18; 3,8,19; 15,21,19; 28,35,19; 33,40,19; 42,43,19; 48,51,19; 36,43,21;...
       45,47,21; 21,22,22; 40,46,22; 43,50,22; 49,56,22; 23,27,23; 30,32,23; 37,38,23; 39,41,23; 46,47,23; 7,8,24; 41,43,24; 44,46,24; 32,36,25; 49,51,25; 56,58,25; 60,64,25; 15,16,26;...
      27,30,26; 35,37,26; 42,44,26; 51,55,26; 63,65,26; 3,6,27; 19,21,27; 24,26,27; 28,29,27; 35,36,27; 37,40,27; 45,46,27; 54,58,27; 59,60,27; 45,47,28; 53,54,28; 55,56,28;...
     8,13,29; 18,23,29; 35,39,29; 43,44,29; 50,51,29; 58,60,29; 14,17,31; 19,20,31; 35,37,31; 42,45,31; 52,54,31; 62,65,31; 16,20,32; 25,28,32; 29,30,32; 33,34,32; 43,46,32;...
    49,51,32; 53,54,32; 59,63,32; 64,65,32; 12,18,33; 17,22,33; 33,38,33; 40,42,33; 45,48,33; 35,37,34; 42,44,34; 53,54,34; 2,5,36; 8,10,36; 16,17,36; 23,24,36; 51,52,37;...
   36,38,38; 59,61,38; 16,19,41; 20,26,41; 24,29,41; 32,33,41; 35,37,41; 39,41,41; 46,48,41; 50,55,41; 35,40,42; 38,45,42; 45,47,42; 46,49,42; 48,51,42; 15,20,44; 21,25,44;...
  30,32,44; 36,40,44; 57,60,44; 12,16,45; 26,27,45; 43,47,45; 62,65,45; 40,42,47; 43,44,47; 45,47,47; 49,51,47; 53,55,47; 58,59,47; 61,62,47; 3,5,48; 7,8,48; 11,14,48;...
 15,16,48; 19,22,48; 28,30,48; 35,37,48; 38,43,48; 47,51,48; 12,14,49; 18,23,49; 29,30,49; 35,37,49; 39,40,49; 43,44,49; 59,62,49; 2,6,50; 37,39,50; 47,52,50; 53,55,50; 64,65,50];
% day of WT peak, day of day fever AggMin, ID - picked to preceed day fevers; unique peaks w/in 7 days from 1st day of fever above threshold (-0.6)); exclude too early to see peaks1
orphans = [21,5; 48,5; 64,7; 36,18; 58,19; 24,22; 27,22; 35,23; 42,27; 44,27; 55,29; 40,31; 40,33; 46,34; 29,44; 62,44; 42,49; 57,50]; %day time fevers from AggMinTn without unique peaks within 1 week prior in WTMaxTn; day, ID 
% 42 & 44,27 is in a long fever where power drops for days; same 46,34; same as 62,44; 
wavepeakDiff = wavepeak(:,2)-wavepeak(:,1);
for i = 1:length(wavepeak)
    if wavepeak(i,1) < 20
        a = 20-wavepeak(i,1);  b = nan(a,1);  c = [b' WTMaxTn(:,wavepeak(i,3))']';
       d = 20-wavepeak(i,2); e = nan(d,1); f = [e' AggMinTn(:,wavepeak(i,3))']';
        WPTMa(:,i) = c(wavepeak(i,1)+a-9:wavepeak(i,1)+a+10); %Wave Peak Temp Max Aligned
        WPTCma(:,i) = f(wavepeak(i,2)+d-12:wavepeak(i,2)+d+7); %Wave Peak Temp Circ Min Aligned
    elseif 65-wavepeak(i,1) < 17
        a = 65-wavepeak(i,1);  b = nan(20,1);  c = [WTMaxTn(:,wavepeak(i,3))' b']';
       d = 65-wavepeak(i,2);  e = nan(20,1); f = [AggMinTn(:,wavepeak(i,3))' e']';
        WPTMa(:,i) = c(wavepeak(i,1)-9:wavepeak(i,1)+10); %Wave Peak Temp Max Aligned
        WPTCma(:,i) = f(wavepeak(i,2)-12:wavepeak(i,2)+7); %Wave Peak Temp Circ Min Aligned
    else
        c = WTMaxTn(:,wavepeak(i,3))';
       d = AggMinTn(:,wavepeak(i,3))';
        WPTMa(:,i) = c(wavepeak(i,1)-9:wavepeak(i,1)+10); %Wave Peak Temp Max Aligned
        WPTCma(:,i) = d(wavepeak(i,2)-12:wavepeak(i,2)+7); %Wave Peak Temp Circ Min Aligned
    end
    clear a b c d e f
end

%Zscores of peak wavelet values, modified so 0 excludes a 10 day window surrounding the peak to allow comparison of peak change pre- and post-peak power for all peaks
for i = 1:length(WPTMa)
WPTMaZ(:,i) = (WPTMa(:,i)-nanmean(WPTMa([1:6 15:20] ,i)))./(nanstd(WPTMa([1:6 15:20],i))); 
end
for i = 1:length(wavepeak)
    WPTamp(i) = WTMaxTn(wavepeak(i,1),wavepeak(i,3)); %peak height from WT
    WPTCamp(i) = AggMinTn(wavepeak(i,2),wavepeak(i,3)); %peak height from temp min
end

%mean(wavepeak(:,2)-wavepeak(:,1)) = 2.84;
a = nanmean(WPTMaZ,2)+nanstd(WPTMaZ,[],2); b = nanmean(WPTMaZ,2)-nanstd(WPTMaZ,[],2);
c = nanmean(WPTCma,2)+nanstd(WPTCma,[],2); d = nanmean(WPTCma,2)-nanstd(WPTCma,[],2);

WTtrace = figure; %aligned peaks in Tmin and daily wavelet band power in Figure 3B,D
subplot(1,6,1:3); hold on; fill([1:1:20 20:-1:1],[a' fliplr(b')],'b','LineStyle','none'); alpha(0.1);
xlim([0 21]); xtick([10 13]); xticklabels({'WT peak','Daytime fever onset'})
subplot(1,6,4:6); hold on; fill([1:1:20 20:-1:1],[c' fliplr(d')],'r','LineStyle','none'); alpha(0.1); 
xlim([0 21]); xtick([10 13]); xticklabels({'WT peak','Daytime fever onset'})
for i = 1:length(wavepeak)
    if wavepeakDiff(i) >= 6
        subplot(1,6,1:3); hold on; plot(WPTMaZ(:,i),'Color',[0 0.2 1],'LineWidth',1);
        subplot(1,6,4:6); hold on; plot(WPTCma(:,i),'Color',[0.9 0 0],'LineWidth',1);
    elseif wavepeakDiff(i) >= 3 && wavepeakDiff(i) < 6
        subplot(1,6,1:3); hold on; plot(WPTMaZ(:,i),'Color',[0 0.2 1],'LineWidth',0.5);
        subplot(1,6,4:6); hold on; plot(WPTCma(:,i),'Color',[0.9 0 0],'LineWidth',0.5);
    elseif wavepeakDiff(i) >= 1 && wavepeakDiff(i) < 3
        subplot(1,6,1:3); hold on; plot(WPTMaZ(:,i),'Color',[0 0.2 1],'LineWidth',0.1);
        subplot(1,6,4:6); hold on; plot(WPTCma(:,i),'Color',[0.9 0 0],'LineWidth',0.1);    
    end
end
subplot(1,6,1:3); plot(nanmean(WPTMaZ,2),'w','LineWidth',8); plot(nanmean(WPTMaZ,2),'Color',[0 0.1 0.75],'LineWidth',4);
subplot(1,6,4:6); plot(nanmean(WPTCma,2),'w','LineWidth',8); plot(nanmean(WPTCma,2),'Color',[0.75 0 0],'LineWidth',4);

xlim([0 8]); ylim([-3 15]); set(gca,'TickDir','out');
%Scatter for Figure 3C and correlation
for i = 1:length(WPTamp)
    hold on
    if isempty(find(FevBo1s == wavepeak(i,3)))
    scatter(wavepeakDiff(i),WPTamp(i),200,'w','filled'); alpha(0.4)
    scatter(wavepeakDiff(i),WPTamp(i),200,SickC(1,:)); 
    else
    scatter(wavepeakDiff(i),WPTamp(i),200,SickC(1,:),'filled'); 
    end
end
xlim([0 8]); ylim([-3 15]); set(gca,'TickDir','out');
% 18 orphans, 226 peak diffs, 244 total; 7% miss rate
[r,p] = corrcoef(wavepeakDiff,WPTamp); %r = 0.35; p = 1e-7

%% Supplemental Figure 1
%generating differences for all variables
%Label conventions: Difference of Medians for Variable (T,H,V,B) by Fever (self report) or FevBot (digital biomarker); 0 = not reported/detected; 1 = reported/detected
DMedTFr0 = nanmean(AggMTzN(F0SymOn:1440*55,Fever0s))-nanmedian(AggMTzN(1:1440*40,Fever0s));
DMedHFr0 = nanmean(AggMHzN(F0SymOn:1440*55,Fever0s))-nanmedian(AggMHzN(1:1440*40,Fever0s));
DMedVFr0 = nanmean(AggMVzN(F0SymOn:1440*55,Fever0s))-nanmedian(AggMVzN(1:1440*40,Fever0s));
DMedBFr0 = nanmean(AggMBzN(F0SymOn:1440*55,Fever0s))-nanmedian(AggMBzN(1:1440*40,Fever0s));

DMedTFB0 = nanmean(AggMTzN(F0SymOn:1440*55,FevBo0s))-nanmedian(AggMTzN(1:1440*40,FevBo0s));
DMedHFB0 = nanmean(AggMHzN(F0SymOn:1440*55,FevBo0s))-nanmedian(AggMHzN(1:1440*40,FevBo0s));
DMedVFB0 = nanmean(AggMVzN(F0SymOn:1440*55,FevBo0s))-nanmedian(AggMVzN(1:1440*40,FevBo0s));
DMedBFB0 = nanmean(AggMBzN(F0SymOn:1440*55,FevBo0s))-nanmedian(AggMBzN(1:1440*40,FevBo0s));

DMedTFr1 = nanmean(AggMTzN(F1SymOn:1440*55,Fever1s))-nanmedian(AggMTzN(1:1440*40,Fever1s));
DMedHFr1 = nanmean(AggMHzN(F1SymOn:1440*55,Fever1s))-nanmedian(AggMHzN(1:1440*40,Fever1s));
DMedVFr1 = nanmean(AggMVzN(F1SymOn:1440*55,Fever1s))-nanmedian(AggMVzN(1:1440*40,Fever1s));
DMedBFr1 = nanmean(AggMBzN(F1SymOn:1440*55,Fever1s))-nanmedian(AggMBzN(1:1440*40,Fever1s));

DMedTFB1 = nanmean(AggMTzN(F1SymOn:1440*55,FevBo1s))-nanmedian(AggMTzN(1:1440*40,FevBo1s));
DMedHFB1 = nanmean(AggMHzN(F1SymOn:1440*55,FevBo1s))-nanmedian(AggMHzN(1:1440*40,FevBo1s));
DMedVFB1 = nanmean(AggMVzN(F1SymOn:1440*55,FevBo1s))-nanmedian(AggMVzN(1:1440*40,FevBo1s));
DMedBFB1 = nanmean(AggMBzN(F1SymOn:1440*55,FevBo1s))-nanmedian(AggMBzN(1:1440*40,FevBo1s));

DErTFr0 = nanstd(DMedTFr0); DErTFB0 = nanstd(DMedTFB0); DErTFr1 = nanstd(DMedTFr1); DErTFB1 = nanstd(DMedTFB1);
DErHFr0 = nanstd(DMedHFr0); DErHFB0 = nanstd(DMedHFB0); DErHFr1 = nanstd(DMedHFr1); DErHFB1 = nanstd(DMedHFB1);
DErVFr0 = nanstd(DMedVFr0); DErVFB0 = nanstd(DMedVFB0); DErVFr1 = nanstd(DMedVFr1); DErVFB1 = nanstd(DMedVFB1);
DErBFr0 = nanstd(DMedTFr0); DErBFB0 = nanstd(DMedBFB0); DErBFr1 = nanstd(DMedBFr1); DErBFB1 = nanstd(DMedBFB1);
StDc = [length(Fever0s).^0.5, length(FevBo0s).^0.5, length(Fever1s).^0.5, length(FevBo1s).^0.5]; 

%Supp Figure 1
subplot(3,1,1)
errorbar([nanmean(DMedHFr0) nanmean(DMedHFB0) nanmean(DMedHFr1) nanmean(DMedHFB1)],[DErHFr0/StDc(1) DErHFB0/StDc(2) DErHFr1/StDc(3) DErHFB1/StDc(4)],'LineStyle','none');
hold on; bar([nanmean(DMedHFr0) nanmean(DMedHFB0) nanmean(DMedHFr1) nanmean(DMedHFB1)],'FaceColor',SickC(2,:));
xticklabels({'Fev Rep 0', 'Fev DBM 0', 'Fev Rep 1', 'Fev DBM 1'});  title('H')

subplot(3,1,2)
errorbar([nanmean(DMedVFr0) nanmean(DMedVFB0) nanmean(DMedVFr1) nanmean(DMedVFB1)],[DErVFr0/StDc(1) DErVFB0/StDc(2) DErVFr1/StDc(3) DErVFB1/StDc(4)],'LineStyle','none');
hold on; bar([nanmean(DMedVFr0) nanmean(DMedVFB0) nanmean(DMedVFr1) nanmean(DMedVFB1)],'FaceColor',SickC(3,:));
xticklabels({'Fev Rep 0', 'Fev DBM 0', 'Fev Rep 1', 'Fev DBM 1'});  title('V')

subplot(3,1,3)
errorbar([nanmean(DMedBFr0) nanmean(DMedBFB0) nanmean(DMedBFr1) nanmean(DMedBFB1)],[DErBFr0/StDc(1) DErBFB0/StDc(2) DErBFr1/StDc(3) DErBFB1/StDc(4)],'LineStyle','none');
hold on; bar([nanmean(DMedBFr0) nanmean(DMedBFB0) nanmean(DMedBFr1) nanmean(DMedBFB1)],'FaceColor',SickC(4,:));
xticklabels({'Fev Rep 0', 'Fev DBM 0', 'Fev Rep 1', 'Fev DBM 1'});  title('B')

%for post-hoc comparisons, only 2 comparisons are considered - 1v3 and 2v4, so correction is redone with bonferroni rather than the automatic, which corrects for all 6 comparisons, though 4 are not relevant here.
[P,ANOVATAB,STATS] = kruskalwallis([DMedHFr0 DMedHFB0 DMedHFr1 DMedHFB1],[ones(1,length(Fever0s))+1 ones(1,length(FevBo0s))+2 ones(1,length(Fever1s))+3 ones(1,length(FevBo1s))+4])
% Source     SS      df     MS      Chi-sq   Prob>Chi-sq - posthoc: Tukey-Kramer: difference from Fr0-Fr1 NS, FB0-FB1, p=0.052
% ------------------------------------------------------
% Groups    8363.3    3   2787.78    9.94      0.0191   
% Error    74936.7   96    780.59                       
% Total    83300     99                                                               
%    1v3 p = 0.1288
%    2v4 p = 0.0213

[P,ANOVATAB,STATS] = kruskalwallis([DMedVFr0 DMedVFB0 DMedVFr1 DMedVFB1],[ones(1,length(Fever0s))+1 ones(1,length(FevBo0s))+2 ones(1,length(Fever1s))+3 ones(1,length(FevBo1s))+4])
% Source     SS      df     MS      Chi-sq   Prob>Chi-sq - posthoc: Tukey-Kramer: difference from Fr0-Fr1 NS, FB0-FB1, p=0.073
% ------------------------------------------------------
% Groups    6508.1    3   2169.37    7.73      0.0518   
% Error    76791.9   96    799.92                       
% Total    83300     99                                                                
%    1v3 p = 0.3334
%    2v4 p = 0.0310

[P,ANOVATAB,STATS] = kruskalwallis([DMedBFr0 DMedBFB0 DMedBFr1 DMedBFB1],[ones(1,length(Fever0s))+1 ones(1,length(FevBo0s))+2 ones(1,length(Fever1s))+3 ones(1,length(FevBo1s))+4])
% Source      SS      df      MS      Chi-sq   Prob>Chi-sq - posthoc: Tukey-Kramer: difference from Fr0-Fr1 sig, FB0-FB1, sig
% --------------------------------------------------------
% Groups   14766.01    3   4922.004   18.27      0.0004   
% Error    63633.99   94    676.957                       
% Total    78400      97                                  
%    1v3 p = 0.0017
%    2v4 p = 0.0154; bonferonni, 2 comparisons

%% Supplemental Figure 2.
%insets of ID 23,33; ploting AggMax[variable]n for all 4 variables.
%Correlation values:
StatAgg = {AggMaxTn(:,FevBo1s), AggMaxHn(:,FevBo1s), AggMaxVn(:,FevBo1s), AggMaxBn(:,FevBo1s), AggMinTn(:,FevBo1s), AggMinHn(:,FevBo1s), AggMinVn(:,FevBo1s), AggMinBn(:,FevBo1s)}; 
for k = 1:length(FevBo1s)
    for i = 1:8
        for j = 1:8
            [r,p] = corrcoef(StatAgg{1,i}(:,k),StatAgg{1,j}(:,k),'rows','complete');
            AggRs(i,j,k) = r(1,2);
        end
    end
end
AggRa = nanmean(AggRs,3); %average R per relationship


