if_gui.m

ifnn

function varargout = IF_gui(varargin)

global NPOINTS_TRANSFER_FUNCTION; 
global TRANSFER_FUNCTION_IMAX; 
global PHASE_BIN_WIDTH; 
global ISI_BIN_WIDTH; 
global POWER_SPECTRUM_NBINS; 
global STIM_PARAMS;
global MODEL_PARAMS;
global SCREENSIZE; global FIGSIZE;

if nargin == 0  % LAUNCH GUIf
    
    NPOINTS_TRANSFER_FUNCTION = 2000; %The number of points used in plots of the transfer function
    TRANSFER_FUNCTION_IMAX = 2; %The maximum value of I in the plot
    PHASE_BIN_WIDTH = 0.2; %The width of a bin in the phase histogram
    ISI_BIN_WIDTH = 1; %The width of the bins used in the ISI histogram
    POWER_SPECTRUM_NBINS = 666; %determines whether the whole power spectrum is drawn with the signal-to-noise ratio
    STIM_PARAMS = get_basestimparams;
    MODEL_PARAMS = get_basemodelparams;
    SCREENSIZE = get(0,'screensize');
    
    fig = openfig(mfilename,'reuse');
    curunits = get(fig,'units'); set(fig,'units','pixels');
    cursize = get(fig,'position');
    set(fig,'units',curunits);
    FIGSIZE = [10, 40, 0.95*(SCREENSIZE(3)-10), (cursize(4) < 0.8*SCREENSIZE(4))*cursize(4) + (cursize(4) >= 0.8*SCREENSIZE(4))*0.8*SCREENSIZE(4)];

	% Generate a structure of handles to pass to callbacks, and store it. 
	handles = guihandles(fig);
	guidata(fig, handles);
    curunits = get(fig,'units'); set(fig,'units','pixels','position',FIGSIZE); set(fig,'units',curunits);
    FIGSIZE = get(fig,'position');
    setaxespos(handles,setpositions(handles));
    if(get(handles.popup_stimulus,'value') ~= 2) 
        set(handles.draw_phase,'enable','off');
        set(handles.draw_SNR,'enable','off');
    end;

    if nargout > 0
		varargout{1} = fig;
	end
    
    
elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

	try
		if (nargout)
			[varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
		else
			feval(varargin{:}); % FEVAL switchyard
		end
	catch
		disp(lasterr);
	end

end

function setaxespos(handles,pos)
global FIGSIZE;    
    xgap = 0.07*FIGSIZE(3); ygap = 0.1*FIGSIZE(4);
    AXESHEIGHT = 0.4*FIGSIZE(4); AXESWIDTH = pos-1.5*xgap; AXESSEP = 0.05*FIGSIZE(4);
    set(handles.Iaxes,'visible','on','position',[xgap, ygap, AXESWIDTH, AXESHEIGHT]);
    set(handles.vaxes,'visible','on','position',[xgap, ygap + AXESHEIGHT + AXESSEP, AXESWIDTH,  AXESHEIGHT]);
    set(get(handles.Iaxes,'xlabel'),'string','time (ms)');
    set(get(handles.vaxes,'ylabel'),'string','membrane voltage (mV)');
    set(get(handles.Iaxes,'ylabel'),'string','current (uA)');
    set(handles.slider2,'position',[xgap,min(0.8,ygap/4),AXESWIDTH,min(1.385,ygap/6)]);
    set(handles.slider3,'position',[min(1.8,xgap/4),ygap,min(4,xgap/6),2*AXESHEIGHT+AXESSEP]);
%| ABOUT CALLBACKS:
%| GUIDE automatically appends subfunction prototypes to this file, and 
%| sets objects' callback properties to call them through the FEVAL 
%| switchyard above. This comment describes that mechanism.
%|
%| Each callback subfunction declaration has the following form:
%| <SUBFUNCTION_NAME>(H, EVENTDATA, HANDLES, VARARGIN)
%|
%| The subfunction name is composed using the object's Tag and the 
%| callback type separated by '_', e.g. 'slider2_Callback',
%| 'figure1_CloseRequestFcn', 'axis1_ButtondownFcn'.
%|
%| H is the callback object's handle (obtained using GCBO).
%|
%| EVENTDATA is empty, but reserved for future use.
%|
%| HANDLES is a structure containing handles of components in GUI using
%| tags as fieldnames, e.g. handles.figure1, handles.slider2. This
%| structure is created at GUI startup using GUIHANDLES and stored in
%| the figure's application data using GUIDATA. A copy of the structure
%| is passed to each callback.  You can store additional information in
%| this structure at GUI startup, and you can change the structure
%| during callbacks.  Call guidata(h, handles) after changing your
%| copy to replace the stored original so that subsequent callbacks see
%| the updates. Type "help guihandles" and "help guidata" for more
%| information.
%|
%| VARARGIN contains any extra arguments you have passed to the
%| callback. Specify the extra arguments by editing the callback
%| property in the inspector. By default, GUIDE sets the property to:
%| <MFILENAME>('<SUBFUNCTION_NAME>', gcbo, [], guidata(gcbo))
%| Add any extra arguments after the last argument, before the final
%| closing parenthesis.

% --------------------------------------------------------------------
function varargout = apply_Callback(h, eventdata, handles, varargin)

[sim membrane noise stimulus showspikes drawphase drawISI drawS drawSNR drawpower] = getparams(handles);
N = ceil(sim.duration/sim.dt);
t = 0;
v = membrane.vr;
I = 0;
hideplots = get(handles.hideplots,'value');
if(~hideplots)
    set(handles.vaxes,'Nextplot','add');
    set(handles.Iaxes,'Nextplot','add');
    lv = plot(t,v,'parent',handles.vaxes,'erasemode','none','color','b');
    lI = plot(t,I,'parent',handles.Iaxes,'erasemode','none','color','k');
end;
if(noise.tnoise) vth = membrane.vth; end;
if(noise.rnoise) vr = membrane.vr; end;

if(drawphase | drawISI | drawSNR | drawpower | drawS) 
    recordspikes = 1;
    tf = [];
else recordspikes = 0;
end;


switch(stimulus.type)
case 1 %Step pulse
    Iext_func = @Iext_step;
    Iext_args = {stimulus};
case 2 %Sine wave
    Iext_func = @Iext_sine;
    Iext_args = {stimulus};
case 3 %Spike train, uniform
    Iext_func = @Iext_uspike;
    Iext_args = {stimulus,sim.dt,membrane.C};
end;

%Sets axes and ticks
TRANS = (get(handles.slider3,'value')-0.5);
ZOOM = 1/(1+1.8*(get(handles.slider2,'value')-0.5));
dE = (membrane.vth-membrane.vr);
switch(stimulus.type)
    case {1,2}
        VMAX = min((membrane.R*stimulus.param1/(1 + 0.33*membrane.R*membrane.C/stimulus.param2)),membrane.vth);
        set(handles.vaxes,'xlim',[0 sim.duration],'ylim',2*dE*TRANS + ...
            ZOOM*[min(VMAX -dE,membrane.vr-noise.rnoise*noise.ramp-dE), max(membrane.vth+noise.tnoise*noise.tamp+dE/2,VMAX + dE/2)]);
        
        set(handles.Iaxes,'xlim',[0 sim.duration],'ylim',(stimulus.type==2)*stimulus.param3 + ...
        stimulus.param1*TRANS + ...
        ZOOM*[-1-((stimulus.param1 < 0)*1.5*abs(stimulus.param1)+ ... 
        noise.dnoise*sqrt(noise.damp)),...
    1+1.5*(stimulus.param1 > 0)*abs(stimulus.param1)+noise.dnoise*sqrt(noise.damp)]);
    case 3
        M = (stimulus.param1*stimulus.param2 - stimulus.param3*stimulus.param4);
        std = sqrt(stimulus.param1^2*stimulus.param2^2 + stimulus.param3^2*stimulus.param4^2)/sim.dt/sqrt(1e3);

        set(handles.vaxes,'xlim',[0 sim.duration],'ylim',2*dE*TRANS + min(membrane.R*M,membrane.vth) + ...
            ZOOM*[membrane.vr-noise.rnoise*noise.ramp-std/sqrt(membrane.R*membrane.C) - dE, ...
            max(membrane.vth+noise.tnoise*noise.tamp+dE/2 + std/sqrt(membrane.R*membrane.C),membrane.vth + dE/2)]);

        
        set(handles.Iaxes,'xlim',[0 sim.duration],'ylim',M*TRANS + ...
        ZOOM*[M-20*std,M+20*std]);       
end;
tick = get(handles.vaxes,'xtick');
ticklabel = {};
for i = 1:length(tick)
    ticklabel{i} = num2str(tick(i));
end;
set(handles.vaxes,'xticklabel',ticklabel);

tick = get(handles.Iaxes,'xtick');
ticklabel = {};
for i = 1:length(tick)
    ticklabel{i} = num2str(tick(i));
end;
set(handles.Iaxes,'xticklabel',ticklabel);
tau = membrane.R*membrane.C;
tf = [];
switch(membrane.type)
case 1 %Integrate-and-fire model

%%%%%%%%%%%%%%%%%%%%%%%%%
%    hh = figure; aa = axes('parent',hh); hold on;
%    lv2 = plot(t,v,'parent',aa,'color','b');
%    lI2 = plot(t,I,'parent',aa,'color','k');
%%%%%%%%%%%%%%%%%%%%%%%%%

    a = 1 - sim.dt/tau;
    tflast = -2*membrane.trefr;
    if(noise.enoise) spikefunc = @hazard; varargspikefunc = {noise.eamp,sim.dt}; else spikefunc = @hardthreshold; varargspikefunc = {}; end; %Sets the method for spiek detection
    switch(noise.dnoise)
        case 0 %No white noise
            for i = 1:N
            tnext = i*sim.dt;
            if(tnext - tflast < membrane.trefr) vnext = membrane.vr; 
            else 
                %Checks for spikes
                if(feval(spikefunc,v,membrane.vth,varargspikefunc{:}))
                tflast = t;
                if(~hideplots) 
                    if(showspikes(1)) set(lv,'Xdata',[t t],'Ydata',[v+50 membrane.vr]); end;
                    if(showspikes(2)) set(lI,'Xdata',[t t],'Ydata',[-10 10]); end;                
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
                    %if(showspikes(1)) set(lv2,'Xdata',cat(2,get(lv2,'xdata'),[t t]),'Ydata',cat(2,get(lv2,'ydata'),[v+50 v-(membrane.vth -membrane.vr)])); end;
                    %if(showspikes(2)) set(lI2,'Xdata',cat(2,get(lI2,'xdata'),[t t]),'Ydata',cat(2,get(lI2,'ydata'),[-10 10])); end;                
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                end;
                v = membrane.vr; %%v-(membrane.vth -membrane.vr); %%The second possibility interpolates a bit but is not appropriate for escape noise
                if(noise.tnoise) membrane.vth = vth + noise.tamp*randn(1); end;
                if(noise.rnoise) membrane.vr = vr + noise.ramp*randn(1); end;
                if(recordspikes) tf = cat(2,tf,t); end;
                end;
                %RK2 algorithm
                F1 = FIF(v,tau,membrane.C,Iext_func,t,Iext_args{:});
                vnext = v + (F1 + FIF(v + sim.dt*F1,tau,membrane.C,Iext_func,t+sim.dt/2,Iext_args{:}))*sim.dt/2;
            end;
            Inext = feval(Iext_func,tnext,Iext_args{:});
            if(~hideplots)
                set(lv,'Xdata',[t tnext],'Ydata',[v vnext]);
                set(lI,'Xdata',[t tnext],'Ydata',[I Inext]); drawnow;
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %set(lv2,'Xdata',cat(2,get(lv2,'xdata'),tnext),'Ydata',cat(2,get(lv2,'ydata'),vnext));
                %set(lI2,'Xdata',cat(2,get(lI2,'xdata'),tnext),'Ydata',cat(2,get(lI2,'ydata'),Inext)); drawnow;
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            end;
            v = vnext; t = tnext; I = Inext;
            end;
        case 1 %White noise
        b = noise.damp*sqrt(sim.dt/tau);
        for i = 1:N
            tnext = i*sim.dt;
            if(tnext - tflast < membrane.trefr) vnext = membrane.vr; 
            else 
                %Checks for spikes
                if(feval(spikefunc,v,varargspikefunc{:}))
                tflast = t;
                if(~hideplots)
                    if(showspikes(1)) set(lv,'Xdata',[t t],'Ydata',[v+50 membrane.vr]); end;
                    if(showspikes(2)) set(lI,'Xdata',[t t],'Ydata',[-10 10]); end;                
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
                    %if(showspikes(1)) set(lv2,'Xdata',cat(2,get(lv2,'xdata'),[t t]),'Ydata',cat(2,get(lv2,'ydata'),[v+50 v-(membrane.vth -membrane.vr)])); end;
                    %if(showspikes(2)) set(lI2,'Xdata',cat(2,get(lI2,'xdata'),[t t]),'Ydata',cat(2,get(lI2,'ydata'),[-10 10])); end;                
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                end;
                v = membrane.vr; %%v-(membrane.vth -membrane.vr); %%The second possibility interpolates a bit but is not appropriate for escape noise
                if(noise.tnoise) membrane.vth = vth + noise.tamp*randn(1); end;
                if(noise.rnoise) membrane.vr = vr + noise.ramp*randn(1); end;
                if(recordspikes) tf = cat(2,tf,t); end;
                end;
                psi = b*randn;
                %SRK2 algorithm
                F1 = FIF(v,tau,membrane.C,Iext_func,t,Iext_args{:});
                vnext = v + (F1 + FIF(v + sim.dt*F1 + psi,tau,membrane.C,Iext_func,t+sim.dt/2,Iext_args{:}))*sim.dt/2 + psi;
            end;
            Inext = feval(Iext_func,tnext,Iext_args{:}) + psi;
            if(~hideplots)
                set(lv,'Xdata',[t tnext],'Ydata',[v vnext]);
                set(lI,'Xdata',[t tnext],'Ydata',[I Inext]); drawnow;
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %set(lv2,'Xdata',cat(2,get(lv2,'xdata'),tnext),'Ydata',cat(2,get(lv2,'ydata'),vnext));
                %set(lI2,'Xdata',cat(2,get(lI2,'xdata'),tnext),'Ydata',cat(2,get(lI2,'ydata'),Inext)); drawnow;
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            end;
            v = vnext; t = tnext; I = Inext;
        end;
        end;

case 2 %Poisson neuron
lambda = membrane.C*sim.dt/1000;
trefr = membrane.R;
if(noise.dnoise) b = noise.damp*sqrt(sim.dt/tau); end;
    tflast = -2*trefr;
    for i = 1:N
    tnext = i*sim.dt;
    if((tnext - tflast) > trefr)
        if(trefr>0) nspikes = min(poissrnd(lambda),1); else nspikes = poissrnd(lambda); end;
    else
        nspikes = 0;
    end;
    if(nspikes>0)
        tflast = tnext;
        for j = 1:length(nspikes)
            if(recordspikes) tf = cat(2,tf,tnext); end;
        end
    end;
    if(~hideplots)        
        if(nspikes>0)
            if(showspikes(1)) set(lv,'Xdata',[t t],'Ydata',[v+50 membrane.vr]); end;
            if(showspikes(2)) set(lI,'Xdata',[t t],'Ydata',[-10 10]); end;                
        end;
            set(lv,'Xdata',[t tnext],'Ydata',[v v]);
            if(noise.dnoise) Inext = feval(Iext_func,tnext,Iext_args{:}) + b*randn; else Inext = feval(Iext_func,tnext,Iext_args{:}); end;
            set(lI,'Xdata',[t tnext],'Ydata',[I  Inext]); drawnow;
            I = Inext;
    end;
            t = tnext;        
    end;

end;

%Draws the histograms and plots
if(length(tf)>1) 
    if(strcmp(get(handles.cumulative,'checked'),'on')) cumulative = 1; else cumulative = 0; end;
    if(drawphase) draw_phase_histogram(tf,stimulus.param2,cumulative); end;
    if(drawISI) draw_ISI_histogram(tf,cumulative); end;
    if(drawS) draw_survivor_function(tf,sim.dt,cumulative); end;
    if(drawSNR)
        %Here there should be only one type of noise
        sigma = 0 ;
        if(noise.rnoise) sigma = noise.ramp; end%/(membrane.vth - stimulus.param3*membrane.R); end;
        if(noise.tnoise) sigma = noise.tamp; end%/(membrane.vth - stimulus.param3*membrane.R); end;
        if(noise.dnoise) sigma = noise.damp; end%/(membrane.vth - stimulus.param3*membrane.R); end;