bsoid_svm.m

function [labels,f_10fps_test] = bsoid_svm(data_test,fps,OF_mdl,smth_hstry,smth_futr)
%BSOID_SVM     Predicts mouse behavior based on the Support Vector Machine (SVM) classifier
%   
%   [LABELS,F_10FPS_TEST] = BSOID_SVM(DATA_TEST,FPS,MDL) outputs classified behaviors using the output from bsoid_mdl.
%
%   INPUTS:
%   DATA_TEST    Matrix of the positons of the 6-body parts outlining the rodent over time videotaped from the bottom looking up. 
%                Rows represents time.
%                Columns 1 & 2 tracks snout; columns 3 to 6 tracks the two front paws; Columns 7 to 10 tracks the two hind paws;
%                Columns 11 & 12 tracks the base of the tail. Tested on tracking data generated by DeepLabCut 2.0.
%   FPS    Rounded frame rate, can use VideoReader/ffmpeg(linux command) to automatically detect the input video fps. Try to match training dataset.
%   OF_MDL    Support Vector Machine Classifier Model. This is the output of bsoid_mdl.
%   SMTH_HSTRY    BOXCAR smoothing using number of frames from before. Default ~40ms before.
%   SMTH_FUTR    BOXCAR smoothing using number of frames from after. Default ~40ms after.
%
%   OUTPUTS:
%   LABELS    Predicted action based on the model, the group here matches the group number in the bsoid_gmm, 10 frame/second temporal resolution.
%   F_10FPS_TEST    The features collated for the test animal, 10 frame/second temporal resolution.
%
%   EXAMPLES:
%   clear data;
%   load test_mouse.mat OF_mdl.mat;
%   [labels,f_10fps_test] = bsoid_svm(data_test,60,OF_mdl);
%   
%   Created by Alexander Hsu, Date: 071919
%   Contact ahsu2@andrew.cmu.edu

    if nargin < 3
        error('Please input test dataset, frame rate, and the SVM model!')
    end
    if nargin < 4
        smth_hstry = round(0.05/(1/fps))-1;
        smth_futr = round(0.05/(1/fps))-1;
    end
    fprintf('Obtaining features from dataset... \n');
    for m = 1:length(data_test) % For each csv file you uploaded.
        %% Obtain features, 2 physical features and 5 time-varying signals
        clear fpd_norm cfp_pt_norm chp_pt_norm sn_pt_norm sn_pt_ang sn_disp pt_disp;
        fpd = data_test{m}(:,3:4)-data_test{m}(:,5:6); % Front Paw distance in x,y
        cfp = [mean([data_test{m}(:,3),data_test{m}(:,5)],2),mean([data_test{m}(:,4),data_test{m}(:,6)],2)]; % Center of front paws x,y
        cfp_pt = [cfp(:,1) - data_test{m}(:,11), cfp(:,2) - data_test{m}(:,12)]; % Center of front paw to proximal tail x,y
        chp = [mean([data_test{m}(:,7),data_test{m}(:,9)],2),mean([data_test{m}(:,8),data_test{m}(:,10)],2)]; % Center of hind paws x,y
        chp_pt = [chp(:,1) - data_test{m}(:,11), chp(:,2) - data_test{m}(:,12)]; % Center of hind paws to proximal tail x,y
        sn_pt = [data_test{m}(:,1) - data_test{m}(:,11), data_test{m}(:,2) - data_test{m}(:,12)]; % Snout to proximal tail x,y
        for i = 1:length(data_test{m}) % Euclidean distance of x,y, since position means nothing
            fpd_norm(i) = norm(data_test{m}(i,3:4)-data_test{m}(i,5:6)); % Front paw R to L euclidean distance
            cfp_pt_norm(i) = norm(cfp_pt(i,:)); % Center of front paws to proximal tail euclidean distance
            chp_pt_norm(i) = norm(chp_pt(i,:)); % Center of hind paws to proximal tail euclidean distance
            sn_pt_norm(i) = norm(sn_pt(i,:)); % Snout to proximal tail euclidean distance, i.e. body length
        end
        fpd_norm_smth{m} = movmean(fpd_norm,[smth_hstry,smth_futr]); % Reduce label noise
        sn_cfp_norm_smth{m} = movmean(sn_pt_norm-cfp_pt_norm,[smth_hstry,smth_futr]); % Reduce label noise
        sn_chp_norm_smth{m} = movmean(sn_pt_norm-chp_pt_norm,[smth_hstry,smth_futr]); % Reduce label noise
        sn_pt_norm_smth{m} = movmean(sn_pt_norm,[smth_hstry,smth_futr]); % Reduce label noise
        for k = 1:length(data_test{m})-1 % Velocity and angle over time
            b_3d = [sn_pt(k+1,:),0]; a_3d = [sn_pt(k,:),0]; c = cross(b_3d,a_3d);
            sn_pt_ang(k) = sign(c(3))*180/pi*atan2(norm(c),dot(sn_pt(k,:),sn_pt(k+1,:))); % Body angle, arctan between body
            sn_disp(k) = norm(data_test{m}(k+1,1:2)-data_test{m}(k,1:2)); % Snout displacement over time
            pt_disp(k) = norm(data_test{m}(k+1,11:12)-data_test{m}(k,11:12)); % Proximal tail displacement over time
        end
        sn_pt_ang_smth{m} = movmean(sn_pt_ang,[smth_hstry,smth_futr]); % Reduce label noise
        sn_disp_smth{m} = movmean(sn_disp,[smth_hstry,smth_futr]); % Reduce label noise
        pt_disp_smth{m} = movmean(pt_disp,[smth_hstry,smth_futr]); % Reduce label noise
        %% Collate 7 features. 
        feats_test{m} = [sn_cfp_norm_smth{m}(:,2:end); sn_chp_norm_smth{m}(:,2:end); fpd_norm_smth{m}(:,2:end); sn_pt_norm_smth{m}(:,2:end); ...
            sn_pt_ang_smth{m}(:,1:end); sn_disp_smth{m}(:,1:end); pt_disp_smth{m}(:,1:end)];
    end
    %% For each test dataset, we will predict the behavior based on the model parameters in MDL at a temporal resolution of 10fps.
    for n = 1:length(feats_test)
        feats1 = [];
        for k = fps/10:fps/10:length(feats_test{n})
            feats1(:,end+1) = [mean(feats_test{n}(1:4,k-fps/10+1:k),2);sum(feats_test{n}(5:7,k-fps/10+1:k),2)];
        end
        f_10fps_test{n} = feats1; % Store individual animal features so you can see what the machine grouping is focusing on.
        labels{n} = predict(OF_mdl,f_10fps_test{n}'); % Predict the labels based on your model.
    end
    
return