% Interactive Dynamic Neuron Simulator
% Sebastian Schneegans
% Ruhr-Universitaet Bochum, 2010



function singleNeuronInteractiveSim

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model parameters (default values) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tau_u = 20; tau_v = 20; % time constants
h_u = -5; h_v = -5; % resting levels of fields
beta_u = 5; beta_v = 5; % steepness of sigmoid output functions

c_uu = 0;
c_vv = 0;
c_vu = 0;
c_uv = 0;

q_u = 0.0; q_v = 0.0;

i_u = 0;
i_v = 0;


%%%%%%%%%%%%%%%%%%%%
% program settings %
%%%%%%%%%%%%%%%%%%%%

% define parameters that will be saved to or loaded from .mat file
saveParamNames = {'tau_u','tau_v','beta_u','beta_v','h_u','h_v',...
  'c_uu','c_vv','c_vu','c_uv','q_u', 'q_v','i_u','i_v'};
saveParamStr = cellToString(saveParamNames);

paramFile = 0;
paramPath = 0;
advParamsChanged = false;
initialize = true;
loop = false;
sliderChanged = [];
quit = false;

historyRange = 500;
changePlotRange = [-10, 10];
changePlotStep = 0.01;

%%%%%%%%%%%%%%
% create GUI %
%%%%%%%%%%%%%%

% note: some changes to the sliders must be done during the initialization
% (setting slider range for fieldSize-dependent values ip1 and ip2)

fig = figure('Position', [50, 100, 1100, 650], 'Name', 'Neural Field', ...
  'Color', 'w', 'NumberTitle', 'off', 'MenuBar', 'none');

% create sliders for model parameters
controlFieldHeight = 0.05;
controlFieldWidth = 0.175;
controlFieldOffsetX = 0.825;
controlFieldOffsetY = 0.05;
sliderWidth = 0.1;
gapWidth = 0.01;
textWidth = controlFieldWidth - sliderWidth - gapWidth;

controlParamNames = {'h_u','c_uu','c_uv','q_u','i_u','h_v','c_vv','c_vu','q_v','i_v'};
controlToolTips = {'resting level of neuron u', 'self-excitation of neuron u', 'connection weight from v to u', ...
  'noise level of field u', 'strength of input to neuron u', 'resting level of neuron v', ...
  'self-excitation of neuron v', 'connection weight from u to v', ...
  'noise level of field v', 'strength of input to neuron v'};
controlPosX = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0] * controlFieldWidth + controlFieldOffsetX;
controlPosY = [17, 16, 15, 14, 13, 7, 6, 5, 4, 3] * controlFieldHeight + controlFieldOffsetY;
controlMin = [-10, 0, -10, 0, 0, -10, 0, -10, 0, 0];
controlMax = [0, 10, 10, 1, 10, 0, 10, 10, 1, 10];
textFormat = {'%0.1f', '%0.1f', '%0.1f', '%0.2f', '%0.1f', ...
  '%0.1f', '%0.1f', '%0.1f', '%0.2f', '%0.1f'};

nControlParams = length(controlParamNames);
sliders = zeros(nControlParams, 1);
textFields = zeros(nControlParams, 1);
for i = 1 : nControlParams
  eval(['tmp = ' controlParamNames{i} ';']);
  sliders(i) = uicontrol(fig, 'Style', 'Slider', 'Units', 'Norm', 'Position', ...
    [controlPosX(i)+textWidth+gapWidth, controlPosY(i), sliderWidth, controlFieldHeight], ...
    'Value', tmp, 'Min', controlMin(i), 'Max', controlMax(i), 'Tooltip', controlToolTips{i}, ...
    'Callback', @sliderCallback);
  textFields(i) = uicontrol(fig,'Style','Text','Units','Norm','HorizontalAlignment', 'left', ...
    'String',[controlParamNames{i} '=' num2str(tmp, textFormat{i})], 'BackgroundColor', 'w',...
    'Position',[controlPosX(i)+gapWidth, controlPosY(i) textWidth controlFieldHeight], ...
    'Tooltip', controlToolTips{i});
end

% create buttons for general functions
buttonWidth = controlFieldWidth/2;
buttonHeight = controlFieldHeight;
buttonOffsetX = controlFieldOffsetX;
buttonOffsetY = 0; %buttonHeight/2;

buttonLabels = {'Advanced', 'Reset', 'Load', 'Save', 'Quit'};
buttonToolTips = {'Advanced model parameters', 'Reset neuron activities', 'Load parameters', 'Save parameters', 'Quit'};
buttonPosX = [0 1 0 1 1] * buttonWidth + buttonOffsetX;
buttonPosY = [2 2 1 1 0] * buttonHeight + buttonOffsetY;

nButtons = length(buttonLabels);
buttons = zeros(nButtons, 1);
for i = 1 : nButtons
  buttons(i) = uicontrol(fig,'Style','Pushbutton','String',buttonLabels{i},...
    'Units','Norm','Position', [buttonPosX(i), buttonPosY(i), buttonWidth, buttonHeight], ...
    'Callback',@buttonCallback,'ToolTip',buttonToolTips{i});
end


% prepare GUI for advanced model parameters
advButtonHeight = 0.25;
advControlFieldHeight = (1 - advButtonHeight) / 5;
advControlFieldWidth = 0.5;
advEditWidth = advControlFieldWidth / 2;
advTextPadding = advControlFieldWidth / 8;
advTextWidth = advControlFieldWidth  - advEditWidth - advTextPadding;

advControlParams = {'tau_u', 'tau_v', 'beta_u', 'beta_v'};
advControlToolTips = {'time constant of neuron u', 'time constant of neuron v', ...
  'steepness of output function for neuron u', 'steepness of output function for neuron v'};
advControlPosX = [0 1 0 1] * advControlFieldWidth;
advControlPosY = [2 2 1 1] * advControlFieldHeight;
advParamsMin = [1 1 0 0]; % minimum values for parameters

nAdvControlParams = length(advControlParams);
advEditFields = zeros(nAdvControlParams);
advTextFields = zeros(nAdvControlParams);
advFig = 0;
advCancelButton = 0;
advOkButton = 0;


% axes parameters
xlim_history = [-historyRange, historyRange/5];
ylim_history = [-10, 10];
xlim_change = ylim_history;
ylim_change = [-0.5, 0.5];


% create axes for field plots
phaseAxes = axes('Position', [0.0375, 0.3, 0.225, 0.4], 'XLim', ylim_history, 'YLim', ylim_history, ...
  'Box', 'on', 'XGrid', 'on', 'YGrid', 'on', 'NextPlot', 'add');
historyAxes_u = axes('Position', [0.3125 0.55 0.2125 0.4], 'XLim', xlim_history, 'YLim', ylim_history, ...
  'Box', 'on', 'XGrid', 'on', 'YGrid', 'on', 'NextPlot', 'add');
changeAxes_u = axes('Position', [0.575 0.55 0.225 0.4], 'XLim', xlim_change, 'YLim', ylim_change, ...
  'Box', 'on', 'XGrid', 'on', 'YGrid', 'on', 'NextPlot', 'add');
historyAxes_v = axes('Position', [0.3125 0.05 0.2125 0.4], 'XLim', xlim_history, 'YLim', ylim_history, ...
  'Box', 'on', 'XGrid', 'on', 'YGrid', 'on', 'NextPlot', 'add');
changeAxes_v = axes('Position', [0.575 0.05 0.225 0.4], 'XLim', xlim_change, 'YLim', ylim_change, ...
  'Box', 'on', 'XGrid', 'on', 'YGrid', 'on', 'NextPlot', 'add');

% plot graphs
histPlot_uv = plot(phaseAxes, zeros(1, historyRange+1)+h_u, zeros(1, historyRange+1)+h_v, '-b');
statePlot_uv = plot(phaseAxes, h_u, h_v, 'ob', 'MarkerFaceColor', 'b', 'MarkerSize', 5, 'LineWidth', 2);
% attPlot_uv = plot(phaseAxes, h_u, h_v, 'dc', 'MarkerSize', 7);
ylabel(phaseAxes, 'activity v');
xlabel(phaseAxes, 'activity u');

inPlot_u = plot(historyAxes_u, (-historyRange:0), zeros(1, historyRange+1)+h_u, '-m');
histPlot_u = plot(historyAxes_u, (-historyRange:0), zeros(1, historyRange+1)+h_u, '-b', 'LineWidth', 2);
actPlot_u = plot(historyAxes_u, 0, h_u, 'ob', 'MarkerFaceColor', 'b', 'MarkerSize', 7);
ylabel(historyAxes_u, 'activity u');
xlabel(historyAxes_u, 'time');

changePlot_u = plot(changeAxes_u, changePlotRange(1):changePlotStep:changePlotRange(2), ...
  -(changePlotRange(1):changePlotStep:changePlotRange(2))+h_u, '-b', 'LineWidth', 2);
attPlot_u = plot(changeAxes_u, h_u, 0, 'dr', 'MarkerSize', 9, 'MarkerFaceColor', 'r');
statePlot_u = plot(changeAxes_u, h_u, 0, 'ob', 'MarkerFaceColor', 'b', 'MarkerSize', 7);
ylabel(changeAxes_u, 'rate of change u');
xlabel(changeAxes_u, 'activity u');

inPlot_v = plot(historyAxes_v, (-historyRange:0), zeros(1, historyRange+1)+h_v, '-m');
histPlot_v = plot(historyAxes_v, (-historyRange:0), zeros(1, historyRange+1)+h_v, '-g', 'LineWidth', 2);
actPlot_v = plot(historyAxes_v, 0, h_v, 'og', 'MarkerFaceColor', 'g', 'MarkerSize', 5);
ylabel(historyAxes_v, 'activity v');
xlabel(historyAxes_v, 'time');

changePlot_v = plot(changeAxes_v, changePlotRange(1):changePlotStep:changePlotRange(2), ...
  -(changePlotRange(1):changePlotStep:changePlotRange(2))+h_v, '-g', 'LineWidth', 2);
attPlot_v = plot(changeAxes_v, h_v, 0, 'dr', 'MarkerFaceColor', 'r', 'MarkerSize', 9);
statePlot_v = plot(changeAxes_v, h_v, 0, 'og', 'MarkerFaceColor', 'g', 'MarkerSize', 7);
ylabel(changeAxes_v, 'rate of change v');
xlabel(changeAxes_v, 'activity u');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% main loop with inialization and simulation %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

while ~quit
  
  % load parameters from file
  if paramFile
    fileParamNames = who('-file', [paramPath paramFile]);
    missingParamNames = setdiff(saveParamNames, fileParamNames);
    if ~isempty(missingParamNames)
      msgbox(['The follwing model parameters were not found in the file: ' ...
        cellToString(missingParamNames) ' Previously assigned values will be used for these parameters.'],...
        'Warning', 'warn');
    end
    loadParamNames = intersect(saveParamNames, fileParamNames);
    if ~isempty(loadParamNames)
      eval(['load ''' paramPath paramFile ''' ' cellToString(loadParamNames)]);
    end
    paramFile = 0;
    initialize = 1;
  end
  
  
  % change advanced parameters
  if advParamsChanged
    paramsOk = true;
    for i = 1 : nAdvControlParams
      tmp = str2double(get(advEditFields(i), 'String'));
      if isnan(tmp) || tmp < advParamsMin(i)
        msgbox(['Value for parameter ' advControlParams{i} ' is outside the valid range.'],...
          'Error', 'error');
        paramsOk = false;
      end
    end
    if paramsOk
      for i = 1 : nAdvControlParams
        tmp = str2double(get(advEditFields(i), 'String'));
        eval([advControlParams{i} ' = tmp;']);
      end
      close(advFig);
      initialize = true;
    end
    advParamsChanged = false;
  end
 
  
  % model initialization
  if initialize
    
    % set values of sliders and text fields
    for i = 1 : nControlParams
      if strcmp(controlParamNames{i}, 'ip1') || strcmp(controlParamNames{i}, 'ip2')
        set(sliders(i), 'Max', fieldSize-1);
      end
      
      eval(['tmp = ' controlParamNames{i} ';']);
      if get(sliders(i), 'Min') > tmp
        set(sliders(i), 'Min', floor(tmp));
      elseif get(sliders(i), 'Max') < tmp
        set(sliders(i), 'Max', ceil(tmp));
      end
      set(sliders(i), 'Value', tmp);
      set(textFields(i), 'String', [controlParamNames{i} '=' num2str(tmp, textFormat{i})]);
    end
    
    % initialization of dynamic neurons
    neuron_u = h_u;
    neuron_v = h_v;
    
    history_u = zeros(1, historyRange+1) + h_u;
    inputHistory_u = zeros(1, historyRange+1) + h_u;
    history_v = zeros(1, historyRange+1) + h_v;
    inputHistory_v = zeros(1, historyRange+1) + h_v;
    
    changePlotStates = changePlotRange(1):changePlotStep:changePlotRange(2);
    changeOutputs_u = sigmoid(changePlotStates, beta_u, 0);
    changeOutputs_v = sigmoid(changePlotStates, beta_v, 0);
            
    initialize = false;
    loop = true;
  end
  
  pause(0.1);
  
  % loop over time steps
  while loop
    output_u = sigmoid(neuron_u, beta_u, 0);
    output_v = sigmoid(neuron_v, beta_v, 0);

    neuron_u = neuron_u + 1/tau_u * (-neuron_u + h_u + i_u + c_uu * output_u + c_uv * output_v) + q_u * randn;
    neuron_v = neuron_v + 1/tau_v * (-neuron_v + h_v + i_v + c_vv * output_v + c_vu * output_u) + q_v * randn;
    
    history_u = [history_u(2:end), neuron_u];
    inputHistory_u = [inputHistory_u(2:end), i_u+h_u];
    history_v = [history_v(2:end), neuron_v];
    inputHistory_v = [inputHistory_v(2:end), i_v+h_v];
    
    change_u = 1/ tau_u * ( -neuron_u + h_u + i_u + c_uu * sigmoid(neuron_u, beta_u, 0) + c_uv * output_v);
    change_v =  1/tau_v * (  -neuron_v + h_v + i_v + c_vv * sigmoid(neuron_v, beta_v, 0) + c_vu * output_u);
    
    changeRates_u = 1/ tau_u * ( -changePlotStates + h_u + i_u + c_uu*changeOutputs_u + c_uv*output_v);
    changeRates_v = 1/tau_v * ( -changePlotStates + h_v + i_v + c_vv*changeOutputs_v + c_vu*output_u);
    
    attractors_u = changePlotStates(diff(changeRates_u<0)==1)+0.5*changePlotStep;
    attractors_v = changePlotStates(diff(changeRates_v<0)==1)+0.5*changePlotStep;
%     attractors_u_exp = repmat(attractors_u, [1, length(attractors_v)]);
%     attractors_v_exp = reshape(repmat(attractors_v, [length(attractors_u), 1]), 1, []);
    
    set(histPlot_uv, 'XData', history_u, 'YData', history_v);
    set(statePlot_uv, 'XData', neuron_u, 'YData', neuron_v);
%     set(attPlot_uv, 'XData', attractors_u_exp, 'YData', attractors_v_exp);
    
    set(inPlot_u, 'YData', inputHistory_u);
    set(histPlot_u, 'YData', history_u);
    set(actPlot_u, 'YData', neuron_u);
    set(changePlot_u, 'YData', changeRates_u);
    set(attPlot_u, 'XData', attractors_u, 'YData', zeros(size(attractors_u)));
    set(statePlot_u, 'XData', neuron_u, 'YData', change_u);
    
    set(inPlot_v, 'YData', inputHistory_v);
    set(histPlot_v, 'YData', history_v);
    set(actPlot_v, 'YData', neuron_v);
    set(changePlot_v, 'YData', changeRates_v);
    set(attPlot_v, 'XData', attractors_v, 'YData', zeros(size(attractors_v)));
    set(statePlot_v, 'XData', neuron_v, 'YData', change_v);


    drawnow;
    pause(0.01);
  end
end
close(fig);


%%%%%%%%%%%%%%%%%%%%%%
% callback functions %
%%%%%%%%%%%%%%%%%%%%%%

function buttonCallback(hObject, eventdata) %#ok<INUSD>
  buttonPressed = find(hObject == buttons);

  switch buttonPressed
    case {1} % Advanced
      loop = false;
      for j = 1 : nControlParams
        set(sliders(j), 'Enable', 'off');
      end
      for j = 1 : nButtons
        set(buttons(j), 'Enable', 'off');
      end
     
      advFig = figure('Position',[150,250,600,120],'Name','Advanced model parameters',...
        'Color','w','NumberTitle','off','MenuBar','none','DeleteFcn',@advCloseCallback);
      for j = 1 : nAdvControlParams
        eval(['tmp = ' advControlParams{j} ';']);
        advEditFields(j) = uicontrol(advFig, 'Style', 'Edit', 'Units', 'Norm', 'Position', ...
          [advControlPosX(j)+advTextWidth+advTextPadding, advControlPosY(j)+advButtonHeight, advEditWidth, advControlFieldHeight], ...
          'String',  num2str(tmp), 'Min', 0, 'Max', 0, 'Tooltip', advControlToolTips{j});
        advTextFields(j) = uicontrol(advFig,'Style','Text','Units','Norm','HorizontalAlignment', 'left', ...
          'String',[advControlParams{j} ':'], 'BackgroundColor', 'w', 'Tooltip', advControlToolTips{j},...
          'Position',[advControlPosX(j)+advTextPadding, advControlPosY(j)+advButtonHeight, advTextWidth, advControlFieldHeight]); %#ok<SETNU>
      end
      advCancelButton = uicontrol(advFig, 'Style', 'Pushbutton', 'String', 'Cancel', 'Units', 'Norm', ...
        'Position', [0.25, 0, 0.25, advButtonHeight], 'Callback', @advCancelCallback); %#ok<SETNU>
      advOkButton = uicontrol(advFig, 'Style', 'Pushbutton', 'String', 'Ok', 'Units', 'Norm', ...
        'Position', [0.5, 0, 0.25, advButtonHeight], 'Callback', @advOkCallback);   %#ok<SETNU>
    case {2} % Reset
      neuron_u = h_u;
      neuron_v = h_v;
    case {3} % Load
      [paramFile, paramPath] = uigetfile('*.mat', 'Load parameter file');
      loop = (paramFile == 0);
    case {4} % Save
      [paramFile, paramPath] = uiputfile('*.mat', 'Save parameter file');
      if paramFile
        eval(['save ''' paramPath paramFile ''' ' saveParamStr ' -v6']);
        paramFile = 0;
      end
    case {5} % Quit
      quit = true;
      loop = false;
  end
end


% update paramter values after slider changed
function sliderCallback(hObject, eventdata)  %#ok<INUSD>
  sliderChanged = find(hObject == sliders);

  paramName = controlParamNames{sliderChanged};
  tmp = get(sliders(sliderChanged), 'Value');
  set(textFields(sliderChanged), 'String', [paramName '=' num2str(tmp, textFormat{sliderChanged})]);
  eval([paramName '= tmp;']);
end


% close advanced parameters windows after cancel
function advCancelCallback(hObject, eventdata) %#ok<INUSD>
  close(advFig);
end


% continue simulation after advanced parameter setting has been cancelled
function advCloseCallback(hObject, eventdata) %#ok<INUSD>
  for j = 1 : nControlParams
    set(sliders(j), 'Enable', 'on');
  end
  for j = 1 : nButtons
    set(buttons(j), 'Enable', 'on');
  end

  loop = true;
end



% update parameters after advanced parameter setting has been confirmed
function advOkCallback(hObject, eventdata) %#ok<INUSD>
  advParamsChanged = true;
end

end


%%%%%%%%%%%%%%%%%%%%%%%
% auxiliary functions %
%%%%%%%%%%%%%%%%%%%%%%%

% create single string from cell array of strings
% elements of cell array are seperated by blanks
function str = cellToString(c)
  str = '';
  for i = 1 : length(c)
    str = [str c{i} ' ']; %#ok<AGROW>
  end
  deblank(str);
end


