fitMLEtoExperiment.m

%  Copyright 2010 Andrew Leifer et al <leifer@fas.harvard.edu>
%  This file is part of Mindcontrol-analysis.
% 
%  Mindcontrol-analysis is free software: you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation, either version 3 of the License, or
%  (at your option) any later version.
% 
%  Mindcontrol-analysis is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%  GNU General Public License for more details.
% 
%  You should have received a copy of the GNU General Public License
%  along with mindcontrol-analysis. If not, see <http://www.gnu.org/licenses/>.




function x=fitMLEtoExperiment(directory,thresh)
% x=fitMLEtoExperiment(directory)
%
% Go through a directory containing .mat files that were exported with
% aggregateReversals. There should be one .mat file for each "run" 
% (which itself may have once been composed of multiple recordings)
% 
% All of these runs should have been of the same experiment. 
%
% This script combines them all, graphs them and fits an exponential.
%
% 
%
% For the most up to date version of this software, see:
% http://github.com/samuellab/mindcontrol
% 
% NOTE: If you use any portion of this code in your research, kindly cite:
% Leifer, A.M., Fang-Yen, C., Gershow, M., Alkema, M., and Samuel A. D.T.,
%   "Optogenetic manipulation of neural activity with high spatial resolution
%   in freely moving Caenorhabditis elegans," Nature Methods, Submitted
%   (2010).
% 


disp('Welcome. Hit enter to close all windows and continue..');
pause;
close all;
files=ls([directory '\*.mat']);

h = waitbar(0,'Aggregating runs...');
steps = size(files,1);

if ~exist('thresh')
    thresh=-1;
end

figure(1);
xlabel('Time (s)');
ylabel('Change in Phase Velocity (worm length/s)');

figure(2);
xlabel('Time (s)');
ylabel('Binary Response');

T=[]; %time;
Q=[]; %response magnitue
R=[]; %binary response

n=0; %number of worms
for k=1:steps
    
    waitbar(k/ steps)

    
    disp(['Loading ' files(k,:) ' ...' ]);

    
    %Load the file
    load( [directory '\'  files(k,:)],'-mat')
    disp('loaded!');
    %Check to make sure we are dealing with a mat file that was generated
    %from the mindcontrol preview program
    if ( ~exist('t') || ~exist('short'))
        disp([directory '\'  files(k,:) '  does not seem to be a file generated by aggregateReversals.m as expected']  )
        continue
    end
    n=n+1; 
    
    
    %calculate binary response
    r=zeros(size(short));
    r(find(short<thresh))=1;
    
    %Concatenate this data to the existing data;
    T=[T t]; %time
    Q=[Q short]; %reversal magnitude
    R=[R r]; %binary response
    
    
    %plot the reversal magnitude
    figure(1);
    hold on;
    plot(t,short,'o','MarkerEdgeColor',[rand, rand, rand],'LineWidth',2);
    title(['Response Habituation (n=' num2str(n) ' worms)']);

    
    %plot binary response
    figure(2)
    hold on;
    plot(t,r,'o','MarkerEdgeColor',[rand, rand, rand],'LineWidth',2);
    title(['Response Habituation (n=' num2str(n) ' worms)']);
    
    
    clear('t','short','r');
end



%Beginning curve fitting
disp('beginning curve fitting');

a=.1;
b=.7;
c=.001;

x0=[a b c];  %initial conditions
lbound=[0,0,0];
ubound=[1,1,10];
 
 

f=@(x) sum(-logLikelihood(x(1),x(2),x(3),T,R),2);

%Add in constraint that a+b<1
coef=[1 1 0];
limit=1;

x = fmincon(f,x0,coef,limit,[],[],lbound,ubound)


figure(2)
hold on;
plot(sort(T),x(1)+x(2)*exp(-x(3)*sort(T)),'r');

%%%%%%%%%%%%
% Another way to visualize this is to plot the ration of responses to non
% responses in a sliding window.
figure(3)

%w=250; %bin width in seconds
nb=10; %number of time bins

%Sort the responses as a function of time
[Tsort Ix]=sort(T);
Rsort=R(Ix);

RsortSum=cumsum(Rsort);
m=0;
for k=1:nb
    %These indices correspond to the location of every 200seconds in Tsort
    TsortIxEven(k)=findClosest(Tsort,k*Tsort(end)/nb);
end
NumEventsPerBin=RsortSum(TsortIxEven)-RsortSum([1 TsortIxEven(1:end-1)]);
NumStimuliPerBin=TsortIxEven - [1 TsortIxEven(1:end-1)] ;
Ratio=NumEventsPerBin./(NumStimuliPerBin);
RatioTimeStamps=  ( Tsort(TsortIxEven)+Tsort([1 TsortIxEven(1:end-1)]) )./2;

%Calculate errorbars based on counting statistics
ebar= ( sqrt(NumEventsPerBin) .* (1-Ratio ) )./(NumStimuliPerBin);
errorbar( RatioTimeStamps, Ratio,ebar,'o'); hold on;

ylim([0 1]);
ylabel('Ratio of Response to Non-Response');
xlabel('Time (s)');
plot(sort(T),x(1)+x(2)*exp(-x(3)*sort(T)),'r');
title({['Ratio of responses to non responses'];[ ' bin size ' num2str(Tsort(end)/nb) ' seconds wide;(n=' num2str(n) ' worms); \tau=' num2str(1/x(3)/60) ' minutes']});




%Calculate errorbars based on counting statistics
figure(4)
ebar= ( sqrt(NumEventsPerBin) .* (1-Ratio ) )./(NumStimuliPerBin);
hold on;
h2=errorbar( RatioTimeStamps, Ratio-x(1),ebar,'o');
set(get(h2,'Parent'), 'YScale', 'log')
semilogy(sort(T),x(2)*exp(-x(3)*sort(T)),'r');



close(h)
disp('Goodbye.');