-
Adriaens, Ines authoredAdriaens, Ines authored
S8b_cpnts.m 18.82 KiB
%% S8_cpnts
% visualisation of the lying bouts + changepoints for development purposes
clear variables
close all
clc
% results directory
init_.resdir = ['C:\Users\adria036\OneDrive - Wageningen '...
'University & Research\iAdriaens_doc\Projects\' ...
'eEllen\B4F_indTracking_cattle\B4F_results\'];
% load data
load([init_.resdir 'D3_sdata_IQ_UWB.mat'])
load([init_.resdir 'D2_data_meta.mat'])
%% set constants
% delete cow 280
sdata = rmfield(sdata,'cow_280');
fields_ = fieldnames(sdata);
for f = 1:length(fields_)
% set cow
cow = fields_{f};
% cow = 'cow_99';
% set technique
cpts.Zstat = 'mean'; % technique based on which cpts are found for z
cpts.CDstat = 'std'; % technique based on which cpts are found for CD
% gapsize for new segment for changepoint analysis
cpts.gapsize = 300; % 5 minutes
% segment length
cpts.segmentlength = 600; % 10 minutes
% set windowsize in days
cpts.window = 600/(60*60*24); % windowsize in which cpts are detected
% set settings
cpts.MinDistance = 1200; % 20 min
cpts.NumChanges = 2; % number of changes per hour of data
% figure settings
cpts.Units = 'centimeters'; % set units
cpts.figPos = [-35 1 30 22]; % set position
%% select data and plot lying bouts
% select data of one day per cow
clear data
close all
% data contains "z" and "CD" and runs from first lying to last lying bout
data = sdata.(cow)(:,[2 3 7 13 17 18 15]);
% reference date for numeric time
refdate = datenum(2019,7,3);
% summary of the lying bouts
data.islying_prev(1) = 0;
data.islying_prev(2:end) = data.islying(1:end-1);
data.lychange = ((data.islying == 1 & data.islying_prev == 0) | ...
(data.islying == 0 & data.islying_prev == 1));
data_ld = data(data.lychange == 1,:);
data_ld.dateup(:,1) = NaT;
data_ld.dateup(1:2:end-1) = data_ld.date(2:2:end);
data_ld(isnat(data_ld.dateup),:) = [];
data_ld = movevars(data_ld,'dateup','After','date');
data_ld.numdatedown = datenum(data_ld.date) - refdate;data_ld.numdateup = datenum(data_ld.dateup) - refdate;
% visualise the data and the lying bouts (data_meta.selection)
h = figure('Units',cpts.Units,...
'OuterPosition',cpts.figPos); % figure
subplot(2,1,1); hold on; box on;title('Z'); % prepare z plot
ylim([0 2.5])
subplot(2,1,2); hold on; box on;title('CD'); % prepare CD plot
ylim([0 12])
% plot lying bouts
for i = 1:2
subplot(2,1,i);
for j = 1:height(data_ld)
x = [datenum(data_ld.date)-refdate ...
datenum(data_ld.dateup)-refdate]; % numeric time
y = [0.1 0.1]; % y variable = thick line
plot(x,y,'LineWidth',4.5,'Color',[204/255 0 102/255])
x = [datenum(data_ld.date)-refdate ... % numeric time
datenum(data_ld.date)-refdate];
y = [0 15]; % vertical line begin
plot(x,y,'LineWidth',0.8,'Color',[204/255 0 102/255])
x = [datenum(data_ld.dateup)-refdate ...
datenum(data_ld.dateup)-refdate]; % vertical line end
plot(x,y,'LineWidth',0.8,'Color',[204/255 0 102/255])
end
end
% plot z data and implement changepoint analysis
% stap 1: determine data with gaps not too long and plot gap
% step 2: split data according to gaps
% step 3: determine max amount of CP based on data
% step 4: CP analysis
% step 5: visualise result
data2 = data(isnan(data.avg_z_sm)~= 1 & isnan(data.centerdist) ~= 1,:);
data2.numdate = datenum(data2.date)-refdate;
% plot z data
subplot(2,1,1); plot(data2.numdate,data2.avg_z_sm, ...
'LineWidth',0.9,'Color',[60/255 179/255 113/255]);
xlim([min(data2.numdate) max(data2.numdate)])
% plot cd data
subplot(2,1,2); plot(data2.numdate,data2.centerdist, ...
'LineWidth',0.9,'Color',[60/255 179/255 113/255]);
xlim([min(data2.numdate) max(data2.numdate)])
% calculate gapsize
data2.gapsize(end) = 0;
data2.gapsize(1:end-1) = diff(data2.numdate)*60*60*24;
% select data segments and characterize the segments
ind = find(data2.gapsize > cpts.gapsize);
segments = data2(ind,:); % segments
segments.segstart(1) = data2.numdate(1); % segment start
segments.segstart(2:end) = data2.numdate(ind(1:end-1)+1); % segment start
segments.dur = (segments.numdate - segments.segstart)*60*60*24; % dur in s
% find the segment where transitions happens
for i = 1:height(data_ld)
% find whether the lie down event is in a segment
ind = find(segments.segstart < data_ld.numdatedown(i) & ...
segments.numdate > data_ld.numdatedown(i));
if ~isempty(ind)
data_ld.segdown(i) = ind;
else
data_ld.segdown(i) = NaN;
end
% find whether the get up event is in a segment
ind = find(segments.segstart < data_ld.numdateup(i) & ...
segments.numdate > data_ld.numdateup(i));
if ~isempty(ind)
data_ld.segup(i) = ind;
else
data_ld.segup(i) = NaN;
end
end
% implement changepoint analysis on the segments that are long enough, and
% that have at least one getup/liedown event
% segment number
segments.number = (1:height(segments))';
% segments with an event get up or lie down
seglist = unique([data_ld.segdown;data_ld.segup]);
seglist(isnan(seglist)) = [];
% add to segments whether there is a get up or lie down event
segments.liedown = ismember(segments.number, seglist);
%%
%---------------------------------------------------------------------
% changepoint analysis - ATTEMPT 2 - take full day independent of segment
% length and a fixed number of max changepoints
% --> evaluate validity of detection only
% if there is a getup or liedown event
days = unique(floor(data2.numdate));
for i = 1:length(days)
% select daydata
data_day = data2(floor(data2.numdate) == days(i),:);
if height(data_day) > 12*60*60
% prepare results
day = sprintf('day_%d',days(i));
% prepare figure
k = figure('Units',cpts.Units,...
'OuterPosition',cpts.figPos); % figure
subplot(2,1,1); hold on; box on; % prepare z plot
title([cow ', day = ' num2str(days(i)) ', Z'], ...
'Interpreter','none');
ylim([0 2.5])
subplot(2,1,2); hold on; box on;title('CD'); % prepare CD plot
ylim([0 12])
% plot z and cd of the segment
subplot(2,1,1);
plot(data_day.numdate, data_day.avg_z_sm, ...
'LineWidth',0.9,'Color',[60/255 179/255 113/255])
subplot(2,1,2);
plot(data_day.numdate, data_day.centerdist, ...
'LineWidth',0.9,'Color',[60/255 179/255 113/255])
% select events
downevents = data_ld(floor(data_ld.numdatedown) == days(i),:); % all down
upevents = data_ld(floor(data_ld.numdatedown) == days(i),:); % all up
% add detection window
downevents.winstart = downevents.numdatedown - cpts.window/2;
downevents.winend = downevents.numdatedown + cpts.window/2;
upevents.winstart = upevents.numdateup - cpts.window/2;
upevents.winend = upevents.numdateup + cpts.window/2;
% plot events as ground thruth + window for j = 1:height(downevents)
% z subplot
subplot(2,1,1);
plot([downevents.numdatedown(j) downevents.numdatedown(j)], ...
[0 2.5],'LineWidth',0.8,'Color',[204/255 0 102/255])
plot([downevents.winstart(j) downevents.winend(j)],...
[1 1],'LineWidth',4.5,'Color',[204/255 0 102/255])
% cd subplot
subplot(2,1,2);
plot([downevents.numdatedown(j) downevents.numdatedown(j)], ...
[0 15],'LineWidth',0.8,'Color',[204/255 0 102/255])
plot([downevents.winstart(j) downevents.winend(j)],...
[6 6],'LineWidth',4.5,'Color',[204/255 0 102/255])
end
for j = 1:height(upevents)
% z subplot
subplot(2,1,1); xlim([days(i) days(i)+1])
plot([upevents.numdateup(j) upevents.numdateup(j)], ...
[0 2.5],'LineWidth',0.8,'Color',[204/255 0 102/255])
plot([upevents.winstart(j) upevents.winend(j)],...
[1 1],'LineWidth',4.5,'Color',[204/255 0 102/255])
% cd subplot
subplot(2,1,2);
xlim([days(i) days(i)+1])
plot([upevents.numdateup(j) upevents.numdateup(j)], ...
[0 15],'LineWidth',0.8,'Color',[204/255 0 102/255])
plot([upevents.winstart(j) upevents.winend(j)],...
[6 6],'LineWidth',4.5,'Color',[204/255 0 102/255])
end
% changepoints avg_z
[cptsr.(cow).(day).ipoints, cptsr.(cow).(day).resid] = ...
findchangepts(data_day.avg_z_sm,...
'Statistic',cpts.Zstat,...
'MinDistance',cpts.MinDistance,...
'MaxNumChanges',cpts.NumChanges*24);
cptsr.(cow).(day).ipoints = [[1;...
cptsr.(cow).(day).ipoints(1:end-1);...
cptsr.(cow).(day).ipoints(end)],...
[cptsr.(cow).(day).ipoints;height(data_day)]];
% calculate mean value + std
for j = 1:size(cptsr.(cow).(day).ipoints,1)
cptsr.(cow).(day).ipoints(j,3) = ...
mean(data_day.avg_z_sm(...
cptsr.(cow).(day).ipoints(j,1):...
cptsr.(cow).(day).ipoints(j,2)));
cptsr.(cow).(day).ipoints(j,4) = ...
std(data_day.avg_z_sm(...
cptsr.(cow).(day).ipoints(j,1):...
cptsr.(cow).(day).ipoints(j,2)));
% plot changepoints avg_z_sm
subplot(2,1,1)
plot([data_day.numdate(cptsr.(cow).(day).ipoints(j,2)) ...
data_day.numdate(cptsr.(cow).(day).ipoints(j,2))],...
[0 3],'b','LineWidth',1)
plot([data_day.numdate(cptsr.(cow).(day).ipoints(j,1)) ...
data_day.numdate(cptsr.(cow).(day).ipoints(j,2))],...
[cptsr.(cow).(day).ipoints(j,3) ...
cptsr.(cow).(day).ipoints(j,3)],'b','LineWidth',2)
end
% changepoints cd
[cptsr.(cow).(day).ipointsCD, cptsr.(cow).(day).residCD] = ...
findchangepts(data_day.centerdist,...
'Statistic',cpts.CDstat,...
'MinDistance',cpts.MinDistance,...
'MaxNumChanges',cpts.NumChanges*24);
cptsr.(cow).(day).ipointsCD = [[1;... cptsr.(cow).(day).ipointsCD(1:end-1);...
cptsr.(cow).(day).ipointsCD(end)],...
[cptsr.(cow).(day).ipointsCD;height(data_day)]];
% calculate mean value + std
for j = 1:size(cptsr.(cow).(day).ipointsCD,1)
cptsr.(cow).(day).ipointsCD(j,3) = mean(data_day.centerdist(...
cptsr.(cow).(day).ipointsCD(j,1):...
cptsr.(cow).(day).ipointsCD(j,2)));
cptsr.(cow).(day).ipointsCD(j,4) = std(data_day.centerdist(...
cptsr.(cow).(day).ipointsCD(j,1):...
cptsr.(cow).(day).ipointsCD(j,2)));
% plot changepoints centerdist
subplot(2,1,2)
plot([data_day.numdate(cptsr.(cow).(day).ipointsCD(j,2)) ...
data_day.numdate(cptsr.(cow).(day).ipointsCD(j,2))],...
[0 15],'b','LineWidth',1)
plot([data_day.numdate(cptsr.(cow).(day).ipointsCD(j,1)) ...
data_day.numdate(cptsr.(cow).(day).ipointsCD(j,2))],...
[cptsr.(cow).(day).ipointsCD(j,3)
cptsr.(cow).(day).ipointsCD(j,3)],'b',...
'LineWidth',2)
end
saveas(k,[init_.resdir 'figb_cpts_' cow '_' day '.fig'])
saveas(k,[init_.resdir 'figb_cpts_' cow '_' day '.tif'])
end
end
end
% save results
save([init_.resdir 'D5b_cpts_res.mat'],'cptsr','cpts')
% clear variables
clear ans cow data data_day daat_ld day days downevents f j i h ind
clear x y upevents seglist data2 k
%% Postprocessing - filter and combine results
% >>> SENSITIVITY
% >>> SPECIFICITY
% - change = in segment
% - difference in mean value (~std of the signal)
% - CD: difference at changepoints in level > 0.25m
% - CD: difference at changepoints in std >
% load results
% load([init_.resdir 'D4_cpts_res.mat'])
% cows
cows = fieldnames(cptsr);
for i = 1:length(fields_) % all cows with results
T=1; % count cpts for cd
T2=1; % count cpts for z
% all days for this cow
days = fieldnames(cptsr.(cows{i}));
close all
% add to changepoints the information for postprocessing
for j = 1:length(days)
% ref day numeric
day = days{j}; day = str2double(day(5:end)); % num day (~refday)
% select and store/split data
data.(cows{i}).(days{j}) = sdata.(cows{i})(:,[2 3 7 13 17 18 15]);
data.(cows{i}).(days{j})(isnan(data.(cows{i}).(days{j}).avg_z_sm) | ...
isnan(data.(cows{i}).(days{j}).centerdist),:) = [];
data.(cows{i}).(days{j}).numdate = ...
datenum(data.(cows{i}).(days{j}).date)-refdate; % numeric date
data.(cows{i}).(days{j}) = data.(cows{i}).(days{j})(...
floor(data.(cows{i}).(days{j}).numdate) == day,:);
% load figure and plot
% k = openfig([init_.resdir 'fig_cpts_' cows{i} '_' days{j} '.fig']);
% results z
rzdata = array2table(cptsr.(cows{i}).(days{j}).ipoints,...
'VariableNames',{'istart','iend','avg','std'});
% difference in mean for CD
rzdata.avgdif(1) = NaN;
rzdata.avgdif(2:end) = abs(diff(rzdata.avg));
% idx = find(rzdata.avgdif <= 0.25); % less than 0.25 meter diff
% if ~isempty(idx)
% rzdata.iend(idx-1) = rzdata.iend(idx);
% rzdata(idx,:) = [];
% end
% calculate mean value + std
for l = 1:height(rzdata)
rzdata.avgnew(l) = mean(...
data.(cows{i}).(days{j}).avg_z_sm(...
rzdata.istart(l):...
rzdata.iend(l)));
rzdata.stdnew(l) = std(...
data.(cows{i}).(days{j}).avg_z_sm(...
rzdata.istart(l):...
rzdata.iend(l)));
rzdata.numdatestart(l) = data.(cows{i}).(days{j}).numdate(...
rzdata.istart(l)); % start
rzdata.numdateend(l) = data.(cows{i}).(days{j}).numdate(...
rzdata.iend(l)); % end
% plot changepoints centerdist
% subplot(2,1,1)
% plot([data.(cows{i}).(days{j}).numdate(rzdata.iend(l)) ...
% data.(cows{i}).(days{j}).numdate(rzdata.iend(l))],...
% [0 2.5],'Color',[1 185/255 0],'LineWidth',1)
% plot([data.(cows{i}).(days{j}).numdate(rzdata.istart(l)) ...
% data.(cows{i}).(days{j}).numdate(rzdata.istart(l))],...
% [0 2.5],'Color',[1 185/255 0],'LineWidth',1)
% plot([data.(cows{i}).(days{j}).numdate(rzdata.istart(l)) ...
% data.(cows{i}).(days{j}).numdate(rzdata.iend(l))],...
% [rzdata.avgnew(l) rzdata.avgnew(l)],...
% 'Color',[1 185/255 0],'LineWidth',2.5)
%
end
% add that it's a changepoint to data
data.(cows{i}).(days{j}).ischangeZ(1) = T2;
data.(cows{i}).(days{j}).ischangeZ(rzdata.iend) = ...
(T2+1:T2+height(rzdata))';
% save new changepoints
fie = ['z_' days{j}];
data.(cows{i}).(fie) = rzdata;
T2 = T2+height(rzdata);
% results cd
rcpdata = array2table(cptsr.(cows{i}).(days{j}).ipointsCD,...
'VariableNames',{'istart','iend','avg','std'});
% difference in mean for CD
rcpdata.avgdif(1) = NaN;
rcpdata.avgdif(2:end) = abs(diff(rcpdata.avg));
% idx = find(rcpdata.avgdif <= 0.5);% less than 0.5 meter diff
% if ~isempty(idx)
% rcpdata.iend(idx-1) = rcpdata.iend(idx);
% rcpdata(idx,:) = [];
% end
% calculate mean value + std
for l = 1:height(rcpdata)
rcpdata.avgnew(l) = mean(...
data.(cows{i}).(days{j}).centerdist(...
rcpdata.istart(l):... rcpdata.iend(l)));
rcpdata.stdnew(l) = std(...
data.(cows{i}).(days{j}).centerdist(...
rcpdata.istart(l):...
rcpdata.iend(l)));
rcpdata.numdatestart(l) = data.(cows{i}).(days{j}).numdate(...
rcpdata.istart(l)); % start
rcpdata.numdateend(l) = data.(cows{i}).(days{j}).numdate(...
rcpdata.iend(l)); % end
% plot changepoints centerdist
% subplot(2,1,2)
% plot([data.(cows{i}).(days{j}).numdate(rcpdata.iend(l)) ...
% data.(cows{i}).(days{j}).numdate(rcpdata.iend(l))],...
% [0 15],'Color',[1 185/255 0],'LineWidth',1)
% plot([data.(cows{i}).(days{j}).numdate(rcpdata.istart(l)) ...
% data.(cows{i}).(days{j}).numdate(rcpdata.istart(l))],...
% [0 15],'Color',[1 185/255 0],'LineWidth',1)
% plot([data.(cows{i}).(days{j}).numdate(rcpdata.istart(l)) ...
% data.(cows{i}).(days{j}).numdate(rcpdata.iend(l))],...
% [rcpdata.avgnew(l) rcpdata.avgnew(l)],...
% 'Color',[1 185/255 0],'LineWidth',2.5)
end
% add that it's a changepoint to data
data.(cows{i}).(days{j}).ischangeCD(1) = T;
data.(cows{i}).(days{j}).ischangeCD(rcpdata.iend) = ...
(T+1:T+height(rcpdata))';
T = T+height(rcpdata);
% save new changepoints
fie = ['cd_' days{j}];
data.(cows{i}).(fie) = rcpdata;
% save figure
% saveas(k, [init_.resdir 'nfig_cpts_' cows{i} '_' days{j} '.fig'])
% saveas(k, [init_.resdir 'nfig_cpts_' cows{i} '_' days{j} '.tif'])
end
end
save([init_.resdir 'D5b_dataday.mat'],'data', '-v7.3')