git-annex in git@b0a251890b7e:/data/tmp/appdata/tmp/local-repo/895

functions/MakePSR04.m

@@ -1 +1,49 @@
-/annex/objects/MD5-s1418--297e275b9ef65bd6556f94c1e4bf4b3f
+function  [PSRout Nout]=MakePSR04(Nrast,Erast,Win,PRMS)
+%Win is the time window around the event
+%PRMS=2 give a PSR cell {neurons,trials} and PRMS=1 collapses all the trials >
+%cell{neurons,1}
+%N is the number of trials
+
+
+
+Otype=PRMS{1};
+FLAG=0;Bmax=0;
+if length(PRMS)>1
+    FLAG=1;
+end
+
+if Otype==1
+    PSRout=cell(size(Nrast,1),1);
+    for i=1:length(Nrast)
+        
+        for k=1:size(Erast,1)
+            tmp=Nrast{i}-Erast(k);
+            if FLAG==0
+                PSRout{i}=cat(1,PSRout{i},tmp(tmp>Win(1)-Bmax & tmp<Win(2)+Bmax));
+            else
+                PSRout{i}=cat(1,PSRout{i},tmp(find(tmp>Win(1)-Bmax & tmp<Win(2)+Bmax,1,PRMS{2})));
+            end
+        end
+        
+    end
+elseif Otype==2
+    
+    PSRout=cell(size(Nrast,1),size(Erast,1)); %PSRout is a cell(neurons,trials) 
+    for i=1:length(Nrast) %loops through the number of neurons
+        
+        for k=1:size(Erast,1) %loops through the trials
+            tmp=Nrast{i}-Erast(k);
+            if ~isempty(find(tmp>Win(1) & tmp<Win(2), 1))
+                if FLAG==0
+                    PSRout{i,k}=tmp(tmp>Win(1)-Bmax & tmp<Win(2)+Bmax);
+                else
+                    PSRout{i,k}=tmp(find(tmp>Win(1)-Bmax & tmp<Win(2)+Bmax,1,PRMS{2}));
+                end
+            else
+                PSRout{i,k}=NaN;
+            end
+        end
+        
+    end
+end
+Nout=size(Erast,1);

functions/MakePSR05.m

@@ -1 +1,51 @@
-/annex/objects/MD5-s1548--a6dc3ec78eb470f813235a7751d67e95
+function  [PSRout Nout]=MakePSR04(Nrast,Erast,Win,PRMS)
+%Win is the time window around the event
+%PRMS=2 give a PSR cell {neurons,trials} and PRMS=1 collapses all the trials >
+%cell{neurons,1}
+%N is the number of trials
+
+
+Otype=PRMS{1};
+FLAG=0;Bmax=0;
+if length(PRMS)>1
+    FLAG=1;
+end
+
+if Otype==1
+    PSRout=cell(size(Nrast,1),1);
+    for i=1:length(Nrast)
+        
+        for k=1:size(Erast,1)
+            tmp=Nrast{i}-Erast(k);
+            if ~isempty(find(tmp>Win(1) & tmp<Win(2)))
+                if FLAG==0
+                    PSRout{i}=cat(1,PSRout{i},tmp(tmp>Win(1)-Bmax & tmp<Win(2)+Bmax));
+                else
+                    PSRout{i}=cat(1,PSRout{i},tmp(find(tmp>Win(1)-Bmax & tmp<Win(2)+Bmax,1,PRMS{2})));
+                end
+            else PSRout{i}=NaN;
+            end         
+        end
+        
+    end
+elseif Otype==2
+    
+    PSRout=cell(size(Nrast,1),size(Erast,1)); %PSRout is a cell(neurons,trials) 
+    for i=1:length(Nrast) %loops through the number of neurons
+        
+        for k=1:size(Erast,1) %loops through the trials
+            tmp=Nrast{i}-Erast(k);
+            if ~isempty(find(tmp>Win(1) & tmp<Win(2)))
+                if FLAG==0
+                    PSRout{i,k}=tmp(tmp>Win(1)-Bmax & tmp<Win(2)+Bmax);
+                else
+                    PSRout{i,k}=tmp(find(tmp>Win(1)-Bmax & tmp<Win(2)+Bmax,1,PRMS{2}));
+                end
+            else
+                PSRout{i,k}=NaN;
+            end
+        end
+        
+    end
+end
+Nout=size(Erast,1);

functions/MakePSR0int.m

@@ -1 +1,55 @@
-/annex/objects/MD5-s1630--8e2c1efa8916bcdf2b9c74112967bb8f
+function  [PSRout, Nout]=MakePSR0int(Nrast,EintS,EintE,PRMS,Int)
+%Eint is the interval start and end
+%PRMS=2 give a PSR cell {neurons,trials} and PRMS=1 collapses all the trials >
+%cell{neurons,1}
+%N is the number of interval
+%Int: 1=iti. 2=trial
+
+inttype=Int{1};
+if inttype==1
+    Timeshift=10;
+elseif inttype==2
+    Timeshift=0;
+end
+
+Otype=PRMS{1};
+FLAG=0;
+if length(PRMS)>1
+    FLAG=1;
+end
+
+if Otype==1
+    PSRout=cell(size(Nrast,1),1);
+    for i=1:length(Nrast)
+        
+        for k=1:size(EintS,1)
+            if FLAG==0
+                tmp=Nrast{i};
+                PSRout{i}=cat(1,PSRout{i},tmp(tmp>EintS(k,1)+Timeshift & tmp<EintE(k,1)-Timeshift));
+            else
+                PSRout{i}=cat(1,PSRout{i},tmp(find(tmp>EintS(k,1)+Timeshift & tmp<EintE(k,1)-Timeshift),1,PRMS{2}));
+            end
+        end
+        
+    end
+elseif Otype==2
+    
+    PSRout=cell(size(Nrast,1),size(EintS,1)); %PSRout is a cell(neurons,intervals) 
+    for i=1:length(Nrast) %loops through the number of neurons
+        
+        for k=1:size(EintS,1) %loops through interval, exclude the 1st interval
+            tmp=Nrast{i};
+            if ~isempty(find(tmp>EintS(k,1) & tmp<EintE(k,1), 1))
+                if FLAG==0
+                    PSRout{i,k}=tmp(tmp>EintS(k,1)+Timeshift & tmp<EintE(k,1)-Timeshift);
+                else
+                    PSRout{i,k}=tmp(find(tmp(tmp>EintS(k,1)+Timeshift & tmp<EintE(k,1)-Timeshift),1,PRMS{2}));
+                end
+            else
+                PSRout{i,k}=NaN;
+            end
+        end
+        
+    end
+end
+Nout=size(EintS,1)-1;

functions/MakePTH07.m

@@ -1 +1,118 @@
-/annex/objects/MD5-s4609--098276466f9a21869f8cf6502f58b7fd
+function [PTHout,BWout,MSRout]=MakePTH06(PSRin,Numin,PRMS)
+%!!!!!!!!CAUTION !!!!!
+%this version was modified by Fred and Ali to use DP and NOT Optimal BW
+%
+%INPUTS:
+%PSRin can be collapsed  or not. if not, the output is uncollapsed as well
+%Numin= number of trials
+%PRMS{1}==1 for making baseline-only PSTH and ==2 for normal, main PSTH %(i.e. -15 to 15s)
+%PRMS{2}==0 for fixed binwidth and ==1 for DP binwidth
+%PRMS{3}=binwidth used
+%OUTPUTS:
+%PTHout
+%BWout is a vector stores the DP BW
+%MSRout is MSR=f(BW)
+
+
+
+
+global Tm Tbase
+
+if PRMS{1}==1
+    T0=Tbase;
+elseif PRMS{1}==2
+    T0=Tm;
+end
+BWx=[0.01:0.005:1]; % Range of binsizes considered
+
+
+PTH=[];
+if size(PSRin,2)==1 %if TRUE, all trials are collapsed
+    if PRMS{2}==0 %fixed/assigned binwidth for all
+        if numel(PRMS{3})==1,
+            BW=repmat(PRMS{3},length(PSRin),1);
+        else
+            BW=PRMS{3};
+        end
+        for k=1:length(PSRin)
+            Tn=[-BW(k)/2:-BW(k):T0(1)-BW(k)]; Tp=[BW(k)/2:BW(k):T0(end)+BW(k)];T=[Tn(end:-1:1),Tp];%making PSTH centered on zero
+            %T=[T0(1)-BW(k):BW(k):T0(end)+BW(k)];
+            TMP=hist(PSRin{k},T)/(BW(k)*Numin(k)); %  make PSTH and normalize by number of refs used (Hz)
+            %deals with low firing in which there is too few DP bins to do
+            %interppolation
+            if length(T)>4
+                PTH(k,:)=interp1(T(2:end-1),TMP(2:end-1),T0,'nearest','extrap');BWout(k,1)=BW(k);
+            else
+                PTH(k,:)=repmat(mean(TMP),1,length(T0));BWout(k,1)=T0(end)-T0(1);
+            end
+            
+        end
+        PTHout=PTH;MSRout=[];
+    else %DP binwidth for all
+        for k=1:length(PSRin)
+            [~,~,MSR]=optimalBW05(PSRin(k),T0); %  make PSTH and normalize by number of refs used (Hz)
+            [~,I]=findDP(BWx,MSR,{2,.9});%{DP=1,DEG/TSH=2,TSHOLD}
+            Tn=[-BWx(I)/2:-BWx(I):T0(1)-BWx(I)]; Tp=[BWx(I)/2:BWx(I):T0(end)+BWx(I)];T=[Tn(end:-1:1),Tp];%making PSTH centered on zero
+            
+            TMP=hist(PSRin{k},T)/(BWx(I)*Numin(k)); %  make PSTH and normalize by number of refs used (Hz)
+            %deals with low firing in which there is too few DP bins to do
+            %interppolation
+            if length(T)>4
+                PTH(k,:)=interp1(T(2:end-1),TMP(2:end-1),T0,'nearest','extrap'); BWout(k,1)=BWx(I);
+            else
+                PTH(k,:)=repmat(mean(TMP),1,length(T0));BWout(k,1)=T0(end)-T0(1);
+            end
+            
+           MSRout(k,:)=MSR;
+        end
+        PTHout=PTH;
+    end
+else % trial-by-trial PTH
+    if PRMS{2}==0 %fixed/assigned binwidth for all
+        if numel(PRMS{3})==1,
+            BW=repmat(PRMS{3},length(PSRin),1);
+        else
+            BW=PRMS{3};
+        end
+        for k=1:size(PSRin,1) %neuron
+            Tn=[-BW(k)/2:-BW(k):T0(1)-BW(k)]; Tp=[BW(k)/2:BW(k):T0(end)+BW(k)];T=[Tn(end:-1:1),Tp];%making PSTH centered on zero
+            %  T=[T0(1)-BW(k):BW(k):T0(end)+BW(k)];
+           
+            for n=1:size(PSRin,2)%trial
+                TMP=hist(PSRin{k,n},T)/BW(k); %  make PSTH and normalize by number of refs used (Hz)
+                %deals with low firing in which there is too few DP bins to do
+                %interppolation
+                if length(T)>4
+                    PTH(k,:,n)=interp1(T(2:end-1),TMP(2:end-1),T0,'nearest','extrap');BWout(k,1)=BW(k);
+                else
+                    PTH(k,:,n)=repmat(mean(TMP),1,length(T0));BWout(k,1)=T0(end)-T0(1);
+                end
+                
+            end
+        end
+        PTHout=PTH;MSRout=[];
+        
+    else %DP binwidth for all
+        for k=1:size(PSRin,1)%neuron
+            [~,~,MSR]=optimalBW05(PSRin(k,:),T0);
+            [~,I]=findDP(BWx,MSR,{2,.9});%{DP=1,DEG/TSH=2,TSHOLD}
+            Tn=[-BWx(I)/2:-BWx(I):T0(1)-BWx(I)]; Tp=[BWx(I)/2:BWx(I):T0(end)+BWx(I)];T=[Tn(end:-1:1),Tp];%making PSTH centered on zero
+            for n=1:size(PSRin,2)%trial
+                TMP=hist(PSRin{k,n},T)/BWx(I); %  make PSTH and normalize by number of refs used (Hz)
+                %deals with low firing in which there is too few DP bins to do
+                %interppolation
+                if length(T)>4
+                    PTH(k,:,n)=interp1(T(2:end-1),TMP(2:end-1),T0,'nearest','extrap');BWout(k,1)=BWx(I);
+                else
+                    PTH(k,:,n)=repmat(mean(TMP),1,length(T0));BWout(k,1)=T0(end)-T0(1);
+                end
+                
+            end
+           MSRout(k,:)=MSR;
+        end
+        PTHout=PTH;
+    end
+    
+end
+
+

functions/classRF_predict.m

@@ -1 +1,62 @@
-/annex/objects/MD5-s2166--7e026fb9b31f99feae58d36b9cf6c2e0
+%**************************************************************
+%* mex interface to Andy Liaw et al.'s C code (used in R package randomForest)
+%* Added by Abhishek Jaiantilal ( abhishek.jaiantilal@colorado.edu )
+%* License: GPLv2
+%* Version: 0.02
+%
+% Calls Classification Random Forest
+% A wrapper matlab file that calls the mex file
+% This does prediction given the data and the model file
+% Options depicted in predict function in http://cran.r-project.org/web/packages/randomForest/randomForest.pdf
+%**************************************************************
+%function [Y_hat votes] = classRF_predict(X,model, extra_options)
+% requires 2 arguments
+% X: data matrix
+% model: generated via classRF_train function
+% extra_options.predict_all = predict_all if set will send all the prediction. 
+%
+%
+% Returns
+% Y_hat - prediction for the data
+% votes - unnormalized weights for the model
+% prediction_per_tree - per tree prediction. the returned object .
+%           If predict.all=TRUE, then the individual component of the returned object is a character
+%           matrix where each column contains the predicted class by a tree in the forest.
+%
+%
+% Not yet implemented
+% proximity
+
+function [Y_new, votes, prediction_per_tree] = classRF_predict(X,model, extra_options)
+    
+    if nargin<2
+		error('need atleast 2 parameters,X matrix and model');
+    end
+    
+    if exist('extra_options','var')
+        if isfield(extra_options,'predict_all') 
+            predict_all = extra_options.predict_all;
+        end
+    end
+    
+    if ~exist('predict_all','var'); predict_all=0;end
+            
+        
+    
+	[Y_hat,prediction_per_tree,votes] = mexClassRF_predict(X',model.nrnodes,model.ntree,model.xbestsplit,model.classwt,model.cutoff,model.treemap,model.nodestatus,model.nodeclass,model.bestvar,model.ndbigtree,model.nclass, predict_all);
+	%keyboard
+    votes = votes';
+    
+    clear mexClassRF_predict
+    
+    Y_new = double(Y_hat);
+    new_labels = model.new_labels;
+    orig_labels = model.orig_labels;
+    
+    for i=1:length(orig_labels)
+        Y_new(find(Y_hat==new_labels(i)))=Inf;
+        Y_new(isinf(Y_new))=orig_labels(i);
+    end
+    
+    1;
+    

functions/classRF_train.m

@@ -1 +1,383 @@
-/annex/objects/MD5-s14829--82a321d0a7c77f33b104acec4394c6ee
+%**************************************************************
+%* mex interface to Andy Liaw et al.'s C code (used in R package randomForest)
+%* Added by Abhishek Jaiantilal ( abhishek.jaiantilal@colorado.edu )
+%* License: GPLv2
+%* Version: 0.02
+%
+% Calls Classification Random Forest
+% A wrapper matlab file that calls the mex file
+% This does training given the data and labels 
+% Documentation copied from R-packages pdf 
+% http://cran.r-project.org/web/packages/randomForest/randomForest.pdf 
+% Tutorial on getting this working in tutorial_ClassRF.m
+%**************************************************************
+% function model = classRF_train(X,Y,ntree,mtry, extra_options)
+% 
+%___Options
+% requires 2 arguments and the rest 3 are optional
+% X: data matrix
+% Y: target values 
+% ntree (optional): number of trees (default is 500). also if set to 0
+%           will default to 500
+% mtry (default is floor(sqrt(size(X,2))) D=number of features in X). also if set to 0
+%           will default to 500
+%
+%
+% Note: TRUE = 1 and FALSE = 0 below
+% extra_options represent a structure containing various misc. options to
+%      control the RF
+%  extra_options.replace = 0 or 1 (default is 1) sampling with or without
+%                           replacement
+%  extra_options.classwt = priors of classes. Here the function first gets
+%                       the labels in ascending order and assumes the
+%                       priors are given in the same order. So if the class
+%                       labels are [-1 1 2] and classwt is [0.1 2 3] then
+%                       there is a 1-1 correspondence. (ascending order of
+%                       class labels). Once this is set the freq of labels in
+%                       train data also affects.
+%  extra_options.cutoff (Classification only) = A vector of length equal to number of classes. The ?winning?
+%                       class for an observation is the one with the maximum ratio of proportion
+%                       of votes to cutoff. Default is 1/k where k is the number of classes (i.e., majority
+%                       vote wins). 
+%  extra_options.strata = (not yet stable in code) variable that is used for stratified
+%                       sampling. I don't yet know how this works. Disabled
+%                       by default
+%  extra_options.sampsize =  Size(s) of sample to draw. For classification, 
+%                   if sampsize is a vector of the length the number of strata, then sampling is stratified by strata, 
+%                   and the elements of sampsize indicate the numbers to be
+%                   drawn from the strata. 
+%  extra_options.nodesize = Minimum size of terminal nodes. Setting this number larger causes smaller trees
+%                   to be grown (and thus take less time). Note that the default values are different
+%                   for classification (1) and regression (5).
+%  extra_options.importance =  Should importance of predictors be assessed?
+%  extra_options.localImp = Should casewise importance measure be computed? (Setting this to TRUE will
+%                   override importance.)
+%  extra_options.proximity = Should proximity measure among the rows be calculated?
+%  extra_options.oob_prox = Should proximity be calculated only on 'out-of-bag' data?
+%  extra_options.do_trace = If set to TRUE, give a more verbose output as randomForest is run. If set to
+%                   some integer, then running output is printed for every
+%                   do_trace trees.
+%  extra_options.keep_inbag Should an n by ntree matrix be returned that keeps track of which samples are
+%                   'in-bag' in which trees (but not how many times, if sampling with replacement)
+%
+% Options eliminated
+% corr_bias which happens only for regression ommitted
+% norm_votes - always set to return total votes for each class.
+%
+%___Returns model which has
+% importance =  a matrix with nclass + 2 (for classification) or two (for regression) columns.
+%       For classification, the first nclass columns are the class-specific measures
+%       computed as mean decrease in accuracy. The nclass + 1st column is the
+%       mean decrease in accuracy over all classes. The last column is the mean decrease
+%       in Gini index. For Regression, the first column is the mean decrease in
+%       accuracy and the second the mean decrease in MSE. If importance=FALSE,
+%       the last measure is still returned as a vector.
+% importanceSD = The ?standard errors? of the permutation-based importance measure. For classification,
+%       a p by nclass + 1 matrix corresponding to the first nclass + 1
+%       columns of the importance matrix. For regression, a length p vector.
+% localImp = a p by n matrix containing the casewise importance measures, the [i,j] element
+%       of which is the importance of i-th variable on the j-th case. NULL if
+%       localImp=FALSE.
+% ntree = number of trees grown.
+% mtry  = number of predictors sampled for spliting at each node.
+% votes (classification only) a matrix with one row for each input data point and one
+%       column for each class, giving the fraction or number of ?votes? from the random
+%       forest.
+% oob_times number of times cases are 'out-of-bag' (and thus used in computing OOB error
+%       estimate)
+% proximity if proximity=TRUE when randomForest is called, a matrix of proximity
+%       measures among the input (based on the frequency that pairs of data points are
+%       in the same terminal nodes).
+% errtr = first column is OOB Err rate, second is for class 1 and so on
+
+function model=classRF_train(X,Y,ntree,mtry, extra_options)
+    DEFAULTS_ON =0;
+    %DEBUG_ON=0;
+
+    TRUE=1;
+    FALSE=0;
+    
+    orig_labels = sort(unique(Y));
+    Y_new = Y;
+    new_labels = 1:length(orig_labels);
+    
+    for i=1:length(orig_labels)
+        Y_new(find(Y==orig_labels(i)))=Inf;
+        Y_new(isinf(Y_new))=new_labels(i);
+    end
+    
+    Y = Y_new;
+    
+    if exist('extra_options','var')
+        if isfield(extra_options,'DEBUG_ON');  DEBUG_ON = extra_options.DEBUG_ON;    end
+        if isfield(extra_options,'replace');  replace = extra_options.replace;       end
+        if isfield(extra_options,'classwt');  classwt = extra_options.classwt;       end
+        if isfield(extra_options,'cutoff');  cutoff = extra_options.cutoff;       end
+        if isfield(extra_options,'strata');  strata = extra_options.strata;       end
+        if isfield(extra_options,'sampsize');  sampsize = extra_options.sampsize;       end
+        if isfield(extra_options,'nodesize');  nodesize = extra_options.nodesize;       end
+        if isfield(extra_options,'importance');  importance = extra_options.importance;       end
+        if isfield(extra_options,'localImp');  localImp = extra_options.localImp;       end
+        if isfield(extra_options,'nPerm');  nPerm = extra_options.nPerm;       end
+        if isfield(extra_options,'proximity');  proximity = extra_options.proximity;       end
+        if isfield(extra_options,'oob_prox');  oob_prox = extra_options.oob_prox;       end
+        %if isfield(extra_options,'norm_votes');  norm_votes = extra_options.norm_votes;       end
+        if isfield(extra_options,'do_trace');  do_trace = extra_options.do_trace;       end
+        %if isfield(extra_options,'corr_bias');  corr_bias = extra_options.corr_bias;       end
+        if isfield(extra_options,'keep_inbag');  keep_inbag = extra_options.keep_inbag;       end
+    end
+    keep_forest=1; %always save the trees :)
+    
+    %set defaults if not already set
+    if ~exist('DEBUG_ON','var')     DEBUG_ON=FALSE; end
+    if ~exist('replace','var');     replace = TRUE; end
+    %if ~exist('classwt','var');     classwt = []; end %will handle these three later
+    %if ~exist('cutoff','var');      cutoff = 1; end    
+    %if ~exist('strata','var');      strata = 1; end
+    if ~exist('sampsize','var');    
+        if (replace) 
+            sampsize = size(X,1); 
+        else
+            sampsize = ceil(0.632*size(X,1));
+        end; 
+    end
+    if ~exist('nodesize','var');    nodesize = 1; end %classification=1, regression=5
+    if ~exist('importance','var');  importance = FALSE; end
+    if ~exist('localImp','var');    localImp = FALSE; end
+    if ~exist('nPerm','var');       nPerm = 1; end
+    %if ~exist('proximity','var');   proximity = 1; end  %will handle these two later
+    %if ~exist('oob_prox','var');    oob_prox = 1; end
+    %if ~exist('norm_votes','var');    norm_votes = TRUE; end
+    if ~exist('do_trace','var');    do_trace = FALSE; end
+    %if ~exist('corr_bias','var');   corr_bias = FALSE; end
+    if ~exist('keep_inbag','var');  keep_inbag = FALSE; end
+    
+
+    if ~exist('ntree','var') | ntree<=0
+		ntree=500;
+        DEFAULTS_ON=1;
+    end
+    if ~exist('mtry','var') | mtry<=0 | mtry>size(X,2)
+        mtry =floor(sqrt(size(X,2)));
+    end
+    
+    addclass =isempty(Y);
+    
+    if (~addclass && length(unique(Y))<2)
+        error('need atleast two classes for classification');
+    end
+    [N D] = size(X);
+    
+    if N==0; error(' data (X) has 0 rows');end
+    
+    if (mtry <1 || mtry > D)
+        DEFAULTS_ON=1;
+    end
+    
+    mtry = max(1,min(D,round(mtry)));
+    
+    if DEFAULTS_ON
+        fprintf('\tSetting to defaults %d trees and mtry=%d\n',ntree,mtry);
+    end
+    
+    if ~isempty(Y)
+        if length(Y)~=N,    
+            error('Y size is not the same as X size');  
+        end
+        addclass = FALSE;
+    else
+        if ~addclass, 
+            addclass=TRUE;
+        end
+        error('have to fill stuff here')
+    end
+    
+    if ~isempty(find(isnan(X)));  error('NaNs in X');   end
+    if ~isempty(find(isnan(Y)));  error('NaNs in Y');   end
+    
+    %now handle categories. Problem is that categories in R are more
+    %enhanced. In this i ask the user to specify the column/features to
+    %consider as categories, 1 if all the values are real values else
+    %specify the number of categories here
+    if exist ('extra_options','var') && isfield(extra_options,'categories')
+        ncat = extra_options.categories;      
+    else
+        ncat = ones(1,D);
+    end
+    
+    maxcat = max(ncat);
+    if maxcat>32
+        error('Can not handle categorical predictors with more than 32 categories');
+    end
+
+    %classRF - line 88 in randomForest.default.R
+    nclass = length(unique(Y));
+    if ~exist('cutoff','var') 
+        cutoff = ones(1,nclass)* (1/nclass);
+    else
+        if sum(cutoff)>1 || sum(cutoff)<0 || length(find(cutoff<=0))>0 || length(cutoff)~=nclass
+            error('Incorrect cutoff specified');
+        end
+    end
+    if ~exist('classwt','var')
+        classwt = ones(1,nclass);
+        ipi=0;
+    else
+        if length(classwt)~=nclass
+            error('Length of classwt not equal to the number of classes')
+        end
+        if ~isempty(find(classwt<=0))
+            error('classwt must be positive');
+        end
+        ipi=1;
+    end
+
+    if ~exist('proximity','var')
+        proximity = addclass;
+        oob_prox = proximity;
+    end
+    
+    if ~exist('oob_prox','var')
+        oob_prox = proximity;
+    end
+    
+    %i handle the below in the mex file
+%     if proximity
+%         prox = zeros(N,N);
+%         proxts = 1;
+%     else
+%         prox = 1;
+%         proxts = 1;
+%     end
+    
+    %i handle the below in the mex file
+    if localImp
+        importance = TRUE;
+%        impmat = zeors(D,N);
+    else
+%        impmat = 1;
+    end
+    
+    if importance
+        if (nPerm<1)
+            nPerm = int32(1);
+        else
+            nPerm = int32(nPerm);
+        end
+        
+        %classRF
+%        impout = zeros(D,nclass+2);
+%        impSD  = zeros(D,nclass+1);
+    else
+%        impout = zeros(D,1);
+%        impSD =  1;
+    end
+    
+    %i handle the below in the mex file
+    %somewhere near line 157 in randomForest.default.R
+    if addclass
+%        nsample = 2*n;
+    else
+%        nsample = n;
+    end
+    
+    Stratify = (length(sampsize)>1);
+    if (~Stratify && sampsize>N) 
+        error('Sampsize too large')
+    end
+    
+    if Stratify
+        if ~exist('strata','var')
+            strata = Y;
+        end
+        nsum = sum(sampsize);
+        if ( ~isempty(find(sampsize<=0)) || nsum==0)
+            error('Bad sampsize specification');
+        end
+    else
+        nsum = sampsize;
+    end
+    %i handle the below in the mex file
+    %nrnodes = 2*floor(nsum/nodesize)+1;
+    %xtest = 1;
+    %ytest = 1;
+    %ntest = 1;
+    %labelts = FALSE;
+    %nt = ntree;
+    
+    
+    
+    
+	%[ldau,rdau,nodestatus,nrnodes,upper,avnode,mbest,ndtree]=
+    %keyboard
+    
+    
+    
+    if Stratify
+        strata = int32(strata);
+    else
+        strata = int32(1);
+    end
+    
+    Options = int32([addclass, importance, localImp, proximity, oob_prox, do_trace, keep_forest, replace, Stratify, keep_inbag]);
+
+    
+    if DEBUG_ON
+        %print the parameters that i am sending in
+        fprintf('size(x) %d\n',size(X));
+        fprintf('size(y) %d\n',size(Y));
+        fprintf('nclass %d\n',nclass);
+        fprintf('size(ncat) %d\n',size(ncat));
+        fprintf('maxcat %d\n',maxcat);
+        fprintf('size(sampsize) %d\n',size(sampsize));
+        fprintf('sampsize[0] %d\n',sampsize(1));
+        fprintf('Stratify %d\n',Stratify);
+        fprintf('Proximity %d\n',proximity);
+        fprintf('oob_prox %d\n',oob_prox);
+        fprintf('strata %d\n',strata);
+        fprintf('ntree %d\n',ntree);
+        fprintf('mtry %d\n',mtry);
+        fprintf('ipi %d\n',ipi);
+        fprintf('classwt %f\n',classwt);
+        fprintf('cutoff %f\n',cutoff);
+        fprintf('nodesize %f\n',nodesize);
+    end    
+    
+    
+    [nrnodes,ntree,xbestsplit,classwt,cutoff,treemap,nodestatus,nodeclass,bestvar,ndbigtree,mtry ...
+        outcl, counttr, prox, impmat, impout, impSD, errtr, inbag] ...
+        = mexClassRF_train(X',int32(Y_new),length(unique(Y)),ntree,mtry,int32(ncat), ... 
+                           int32(maxcat), int32(sampsize), strata, Options, int32(ipi), ...
+                           classwt, cutoff, int32(nodesize),int32(nsum));
+ 	model.nrnodes=nrnodes;
+ 	model.ntree=ntree;
+ 	model.xbestsplit=xbestsplit;
+ 	model.classwt=classwt;
+ 	model.cutoff=cutoff;
+ 	model.treemap=treemap;
+ 	model.nodestatus=nodestatus;
+ 	model.nodeclass=nodeclass;
+ 	model.bestvar = bestvar;
+ 	model.ndbigtree = ndbigtree;
+    model.mtry = mtry;
+    model.orig_labels=orig_labels;
+    model.new_labels=new_labels;
+    model.nclass = length(unique(Y));
+    model.outcl = outcl;
+    model.counttr = counttr;
+    if proximity
+        model.proximity = prox;
+    else
+        model.proximity = [];
+    end
+    model.localImp = impmat;
+    model.importance = impout;
+    model.importanceSD = impSD;
+    model.errtr = errtr';
+    model.inbag = inbag;
+    model.votes = counttr';
+    model.oob_times = sum(counttr)';
+ 	clear mexClassRF_train
+    %keyboard
+    1;
+

functions/myfft.m

@@ -1 +1,95 @@
-/annex/objects/MD5-s2887--a7d792929f2c1d67cad4972bb682b41c
+function [X,f,P,Nt]=myfft(x,dtORfs,plotSpec)
+% [X,f,P,Nt]=myfft(x,dtORfs,plotSpec)
+% 
+% Inputs
+% x       : time domain sigmal
+% dtORfs  : sampling period or frequency. Currently assumes that dtORfs <=1
+%           is a sampling period, whereas dtORfs >1 is a sampling frequency. 
+% plotSpec: =1 => show plot; not specified, or =0 => don't plot
+%
+% Outputs
+% X : fft(x)
+% f : frequency range (Hz)
+% P : power spectrum of x
+% Nt: time stamp for x (based on your input dtORfs).
+%
+% Example: 
+% dt=0.01; t=0:dt:5;
+% x=7*sin(2*pi*10*t) + 19*sin(2*pi*23.5*t);
+% [X,f,P,Nt]=myfft(x,dtORfs,1);
+% 
+% Note in this example that the multiplier of t inside the sin function (i.e,. 2*pi*10)
+% is frequency in radians/sec (omega). Recall that f = 2*pi*omega. f is in Hz.
+% So in the example above, the two frequencies are 10Hz or 2*pi*10 rad/s;
+% and 23.5Hz or 2*pi*23.5 rad/s
+%
+% Shreesh P. Mysore
+% shreesh@stanford.edu
+% 2.10.07
+%
+
+
+xhat=[];
+if nargin<3, plotSpec=0; end
+
+%-- dt or fs?
+if dtORfs > 1, Fs = dtORfs; dt = 1/Fs; else dt = dtORfs; Fs =1/dt; end 
+% min sampling frequency is 1Hz just as a way to identify what is being inputted.
+
+%-- fft
+N=length(x); Nt = dt*(1:N);
+Nfft = 2^nextpow2(N); %N=2^round(log2(N));
+%X=fft(x,Nfft)/N;Pxx=2*abs(X);P=Pxx(1:Nfft/2);
+
+X=fft(x,N)/N;Pxx=2*abs(X);
+%<==> 
+%X = fft(x);Pxx = 2*sqrt(X.*conj(X))/N; %Pxx = X.*conj(X)/N;
+
+%-- power spectrum %magnitude of complex Fourier coefficients)
+P=Pxx(1:floor(N/2)+1,:);
+%P=0.5*Pxx; % all points in Pxx
+
+% %-- phase of coefficients
+phaseAngle = unwrap (angle (X(1:floor(N/2)+1))); %in radians
+phaseAngle = phaseAngle*180/pi; %in degrees
+
+
+%-- frequency range
+f = Fs/2*linspace(0,1,Nfft/2);
+f = Fs/2*linspace(0,1,floor(N/2)+1);
+f = ((0: 1/(floor(N/2)+1): 1-1/(floor(N/2)+1))*Fs/2).'; 
+%[length(f) N]
+% pause
+%f = ((0: 1/N: 1-1/N)*Fs).'; %all points in Pxx 
+
+
+%-- plotting
+if plotSpec==1,
+    figure
+    %  ind = find(f>0 & f~=10); P(ind) = zeros(size(ind));
+    %  X(ind) = zeros(size(ind));
+    %  X = [X(1:Nfft/2) fliplr(X(1:Nfft/2))];
+    %xhat = ifft(N*X);
+    
+    subplot(3,1,1), plot(Nt,x);xh=xlabel('time (seconds)');
+    set(gca,'fontsize',12)
+    yh=ylabel('x'); set(xh,'fontsize',15);set(yh,'fontsize',15)
+    axis tight;
+    
+    subplot(3,1,2), plot(f,P,'.-'), 
+    set(gca,'fontsize',12)
+    xh=xlabel('frequency (Hz)'); yh=ylabel('| fft(x) |'); %th = title('Frequency spectrum of X');
+    set(xh,'fontsize',15);set(yh,'fontsize',15); %set(th,'fontsize',20);
+    %subplot(3,1,3), plot(Nt,xhat(1:N));title('x'), xlabel('time (seconds)')
+    set(gcf,'color','w')
+    
+%     subplot(3,1,3), 
+%     plot(f,phaseAngle,'-'), 
+%     set(gca,'fontsize',12); set(gca,'ylim',[-180 180])
+%     xh=xlabel('frequency (Hz)'); yh=ylabel('phase (deg)');
+%     set(xh,'fontsize',15);set(yh,'fontsize',15); 
+%     %subplot(3,1,3), plot(Nt,xhat(1:N));title('x'), xlabel('time (seconds)')
+%     set(gcf,'color','w')
+   
+end
+

functions/myfind.m

@@ -1 +1,21 @@
-/annex/objects/MD5-s388--41b148963ab68d2bec4fa6b8a3b2234b
+function Eind=myfind(EVENTS,Ename)
+Eind=[];
+if isfield(EVENTS{1},'name')
+    for k=1:length(EVENTS)
+        
+        if strcmpi(EVENTS{k}.name,Ename)
+            Eind=[Eind,k];
+        end
+        
+    end
+elseif iscell(EVENTS)
+    for k=1:length(EVENTS)
+        
+        if strcmpi(EVENTS{k},Ename)
+            Eind=[Eind,k];
+        end
+        
+    end
+    
+end
+end

functions/normalize.m

@@ -1 +1,23 @@
-/annex/objects/MD5-s748--4597c5d01587b0300d9c4137e49113a6
+function Mnormal=normalize(M,Mbase,MD)
+
+%MD is the type of normalization. 
+%   0 = Z-score  >>(X-mean)/SD
+%   1 = normalized to response max
+%   2 = Baseline substracted
+%   3 = normalized by SD only >> X/SD
+
+%Added April 2012 to handle vectors
+if size(M,2)==1, M=M'; Mbase=Mbase'; end
+Mnormal=zeros(size(M));
+for k=1:size(M,1)
+    if MD==0
+        Mnormal(k,:)=(M(k,:)-nanmean(Mbase(k,:)))/nanstd(Mbase(k,:)); % Zscore
+    elseif MD==1
+        Tmp=M(k,:)-nanmean(Mbase(k,:));
+        Mnormal(k,:)=Tmp/max(abs(Tmp));% mean corrected normalized to RESPONSE MAX
+    elseif MD==2
+       Mnormal(k,:)=(M(k,:)-nanmean(Mbase(k,:)));
+    elseif MD==3
+            Mnormal(k,:)=M(k,:)/nanstd(Mbase(k,:)); % Zscore no demean
+    end
+end