From 954be97f644fc6838e484f9716169c036794c900 Mon Sep 17 00:00:00 2001
From: Yarmo Mackenbach <yarmo@protonmail.com>
Date: Wed, 1 Apr 2020 14:44:25 +0200
Subject: [PATCH] Added matlab files

---
 matlab/@pricklypear/anadata.m              |  59 +++++++++
 matlab/@pricklypear/basedir.m              |  18 +++
 matlab/@pricklypear/child.m                |  11 ++
 matlab/@pricklypear/datadir.m              |  15 +++
 matlab/@pricklypear/delete.m               |  91 +++++++++++++
 matlab/@pricklypear/directdir.m            |   9 ++
 matlab/@pricklypear/disp.m                 |  32 +++++
 matlab/@pricklypear/exist.m                |  73 +++++++++++
 matlab/@pricklypear/experiment.m           |   7 +
 matlab/@pricklypear/irec.m                 |   7 +
 matlab/@pricklypear/modeldir.m             |   9 ++
 matlab/@pricklypear/notes.m                |  54 ++++++++
 matlab/@pricklypear/parent.m               |  11 ++
 matlab/@pricklypear/plot.m                 |  83 ++++++++++++
 matlab/@pricklypear/pricklypear.m          | 146 +++++++++++++++++++++
 matlab/@pricklypear/run.m                  |  91 +++++++++++++
 matlab/@pricklypear/runSubthresholdModel.m |  83 ++++++++++++
 matlab/@pricklypear/save.m                 |  44 +++++++
 matlab/@pricklypear/spiketimes.m           |  25 ++++
 matlab/@pricklypear/strpad.m               |  60 +++++++++
 matlab/@pricklypear/struct.m               |  24 ++++
 matlab/@pricklypear/subsref.m              |  20 +++
 matlab/@pricklypear/x.m                    |   7 +
 matlab/flagParser.m                        |  31 +++++
 matlab/pp_dir.template.txt                 |   1 +
 25 files changed, 1011 insertions(+)
 create mode 100644 matlab/@pricklypear/anadata.m
 create mode 100644 matlab/@pricklypear/basedir.m
 create mode 100644 matlab/@pricklypear/child.m
 create mode 100644 matlab/@pricklypear/datadir.m
 create mode 100644 matlab/@pricklypear/delete.m
 create mode 100644 matlab/@pricklypear/directdir.m
 create mode 100644 matlab/@pricklypear/disp.m
 create mode 100644 matlab/@pricklypear/exist.m
 create mode 100644 matlab/@pricklypear/experiment.m
 create mode 100644 matlab/@pricklypear/irec.m
 create mode 100644 matlab/@pricklypear/modeldir.m
 create mode 100644 matlab/@pricklypear/notes.m
 create mode 100644 matlab/@pricklypear/parent.m
 create mode 100644 matlab/@pricklypear/plot.m
 create mode 100644 matlab/@pricklypear/pricklypear.m
 create mode 100644 matlab/@pricklypear/run.m
 create mode 100644 matlab/@pricklypear/runSubthresholdModel.m
 create mode 100644 matlab/@pricklypear/save.m
 create mode 100644 matlab/@pricklypear/spiketimes.m
 create mode 100644 matlab/@pricklypear/strpad.m
 create mode 100644 matlab/@pricklypear/struct.m
 create mode 100644 matlab/@pricklypear/subsref.m
 create mode 100644 matlab/@pricklypear/x.m
 create mode 100644 matlab/flagParser.m
 create mode 100644 matlab/pp_dir.template.txt

diff --git a/matlab/@pricklypear/anadata.m b/matlab/@pricklypear/anadata.m
new file mode 100644
index 0000000..08d503c
--- /dev/null
+++ b/matlab/@pricklypear/anadata.m
@@ -0,0 +1,59 @@
+function output = anadata(obj, varargin)
+%ANADATA Get the analog data from a PricklyPear run
+%   anadata(PP) returns the analog data (membrane potential) of the soma
+%   from a PricklyPear simulation run (therefore, PP needs an experiment
+%   and irec).
+%
+%   anadata(PP, channel) specifies the recording channel (c, d0d, ad...).
+%
+%   anadata(PP, channel, metric) specifies the metric (vm).
+%
+%   anadata(..., 'spt', spt) removes the spikes at spt (ms) by using
+%   linear interpolation over a [-1 1.2] ms window.
+%   
+%   anadata(..., 'spwin', [-1 1.2]) specifies the time window (ms) to be used
+%   for linear interpolation when suppressing spikes.
+
+    % Handle input
+    p = inputParser();
+    p.addOptional('channel', 'c', @ischar);
+    p.addOptional('metric', 'vm', @ischar);
+    p.addParamValue('spt', []);
+    p.addParamValue('spwin', [-1 1.2]);
+    p.parse(varargin{:});
+    p = p.Results;
+    
+    % The object must be a specific run
+    assert(~isempty(experiment(obj)), 'Please specify an experiment');
+    assert(~isempty(irec(obj)), 'Please specify a specific run');
+    
+    % Requested variables must be correct
+    assert(any(strcmp([obj.sim_channels {'all'}], p.channel)), sprintf('"%s" is not a valid channel', p.channel));
+    assert(any(strcmp([obj.sim_metrics {'all'}], p.metric)), sprintf('"%s" is not a valid metric', p.metric));
+    
+    % Handle output
+    output = {};
+    
+    % Generate the filename
+    if obj.irec == 0
+        fname = fullfile(obj.datadir(), sprintf('run_unsaved__vm_%s.txt', p.channel));
+    else
+        fname = fullfile(obj.datadir(), sprintf('run_%05i__vm_%s.txt', obj.irec, p.channel));
+    end
+    
+    % Read it if it exists
+    if exist(fname, 'file')
+        txt = fileread(fname);
+        data = textscan(txt, '%f', 'Delimiter', '\n');
+        output = data{1};
+        output(end) = [];
+    else
+        error('Data was not found');
+    end
+    
+    % Remove spikes
+    if ~isempty(p.spt)
+        output = cutFromWaveform(obj.dt, output, p.spt, p.spwin);
+    end
+    
+end
diff --git a/matlab/@pricklypear/basedir.m b/matlab/@pricklypear/basedir.m
new file mode 100644
index 0000000..562b9cd
--- /dev/null
+++ b/matlab/@pricklypear/basedir.m
@@ -0,0 +1,18 @@
+function output = basedir(obj)
+%DATADIR Get the pricklypear base directory
+%   directdir(PP)
+%
+%   See also pricklypear/datadir, pricklypear/modeldir
+    
+    fn = mfilename('fullpath');
+    fn = fileparts(fn);
+    fn = fileparts(fn);
+    fn = fullfile(fn, 'pp_dir.txt');
+    
+    fid = fopen(fn,'r');
+    txt = fscanf(fid, '%s');
+    fclose(fid);
+    
+    output = txt;
+
+end
diff --git a/matlab/@pricklypear/child.m b/matlab/@pricklypear/child.m
new file mode 100644
index 0000000..e13b937
--- /dev/null
+++ b/matlab/@pricklypear/child.m
@@ -0,0 +1,11 @@
+function output = child(obj, irec)
+%CHILD Spawn a child of a pricklypear object
+%   PP.irec -> []
+%   P = child(PP, 2)
+%   P.irec -> 2
+%
+%   See also pricklypear/parent
+
+    output = pricklypear(experiment(obj), irec);
+
+end
diff --git a/matlab/@pricklypear/datadir.m b/matlab/@pricklypear/datadir.m
new file mode 100644
index 0000000..b2718c3
--- /dev/null
+++ b/matlab/@pricklypear/datadir.m
@@ -0,0 +1,15 @@
+function output = datadir(obj)
+%DATADIR Get the directory where data is stored
+%   datadir(PP) returns the directory where data is stored if the object
+%   has a specified experiment, otherwise it returns the direct directory
+%   (where the model directly outputs its data).
+%
+%   See also pricklypear/directdir, pricklypear/modeldir
+
+    if isempty(obj.sim_experiment)
+        output = obj.directdir();
+    else
+        output = fullfile(obj.basedir, 'data', obj.sim_experiment);
+    end
+
+end
diff --git a/matlab/@pricklypear/delete.m b/matlab/@pricklypear/delete.m
new file mode 100644
index 0000000..62c2678
--- /dev/null
+++ b/matlab/@pricklypear/delete.m
@@ -0,0 +1,91 @@
+function output = delete(obj, varargin)
+%DELETE Delete the data of a run
+%   delete(PP) deletes the data of run PP if PP has a specified experiment
+%   and irec.
+
+    % Handle input
+    p = inputParser();
+    p.addParamValue('directMode', false);
+    p.parse(varargin{:});
+    p = p.Results;
+    
+    % Direct mode
+    fpath = obj.datadir();
+    irec = obj.irec;
+    if p.directMode
+        fpath = obj.directdir();
+        irec = obj.iter;
+    end
+
+    % The object must be a specific run
+    assert(~isempty(experiment(obj)), 'Please specify an experiment');
+    assert(~isempty(irec), 'Please specify a specific run');
+    
+    % Handle output
+    isDeleted = 0;
+    
+    % Variables
+    flExt = obj.sim_channels;
+    
+    % Delete every single channel known
+    for ii = 1:numel(flExt)
+        % Generate a filename
+        if irec == 0
+            fname = fullfile(fpath, sprintf('run_unsaved__vm_%s.txt', flExt{ii}));
+        else
+            fname = fullfile(fpath, sprintf('run_%05i__vm_%s.txt', irec, flExt{ii}));
+        end
+        
+        % Check if file exists
+        if exist(fname, 'file')
+            delete(fname);
+            isDeleted = 1;
+        end
+    end
+    % Delete other files
+    if irec == 0
+        fname = fullfile(fpath, sprintf('run_unsaved__evt_contra.txt'));
+        if exist(fname, 'file')
+            delete(fname);
+            isDeleted = 1;
+        end
+        fname = fullfile(fpath, sprintf('run_unsaved__evt_ipsi.txt'));
+        if exist(fname, 'file')
+            delete(fname);
+            isDeleted = 1;
+        end
+        fname = fullfile(fpath, sprintf('run_unsaved__properties.ppbin'));
+        if exist(fname, 'file')
+            delete(fname);
+            isDeleted = 1;
+        end
+    else
+        fname = fullfile(fpath, sprintf('run_%i__evt_contra.txt', irec));
+        if exist(fname, 'file')
+            delete(fname);
+            isDeleted = 1;
+        end
+        fname = fullfile(fpath, sprintf('run_%i__evt_ipsi.txt', irec));
+        if exist(fname, 'file')
+            delete(fname);
+            isDeleted = 1;
+        end
+        fname = fullfile(fpath, sprintf('run_%i__properties.ppbin', irec));
+        if exist(fname, 'file')
+            delete(fname);
+            isDeleted = 1;
+        end
+    end
+    
+    % Output
+    output = isDeleted;
+    if nargout == 0
+        switch isDeleted
+            case 0
+                fprintf('Run %i is not deleted\n', irec);
+            case 1
+                fprintf('Run %i is deleted\n', irec);
+        end
+    end
+
+end
diff --git a/matlab/@pricklypear/directdir.m b/matlab/@pricklypear/directdir.m
new file mode 100644
index 0000000..5f8fa66
--- /dev/null
+++ b/matlab/@pricklypear/directdir.m
@@ -0,0 +1,9 @@
+function output = directdir(obj)
+%DATADIR Get the directory where the model directly stored its output data
+%   directdir(PP)
+%
+%   See also pricklypear/datadir, pricklypear/modeldir
+
+    output = fullfile(obj.basedir, 'data', '_direct');
+
+end
diff --git a/matlab/@pricklypear/disp.m b/matlab/@pricklypear/disp.m
new file mode 100644
index 0000000..e7ae5ee
--- /dev/null
+++ b/matlab/@pricklypear/disp.m
@@ -0,0 +1,32 @@
+function disp(obj)
+%DISP Display information about a PricklyPear object
+%   disp(PP)
+
+    disp(' ');
+    disp('==================');
+    disp('PricklyPear object');
+    disp('==================');
+    disp(' ');
+    
+    fprintf('Experiment: %s\n', obj.sim_experiment);
+    fprintf('      Irec: %i\n', obj.sim_irec);
+    disp(' ');
+    
+    p = properties(obj);
+    pt = obj.strpad(p);
+    
+    disp('Properties');
+    disp('==========');
+    for ii = 1:numel(pt)
+        if any(strcmp(obj.template_exclude, p{ii}))
+            continue;
+        end
+        fprintf('    %s -> %g\n', pt{ii}, obj.(p{ii}));
+    end
+    disp(' ');
+    
+%     disp('Methods');
+%     disp('=======');
+%     fprintf('    run\n');
+
+end
diff --git a/matlab/@pricklypear/exist.m b/matlab/@pricklypear/exist.m
new file mode 100644
index 0000000..3b0ec62
--- /dev/null
+++ b/matlab/@pricklypear/exist.m
@@ -0,0 +1,73 @@
+function output = exist(obj, varargin)
+%EXIST Check if a run exists
+%   exist(PP) checks the existence of run PP if PP has a specified
+%   experiment and irec.
+%
+%   exist(PP, 'full', true) checks if all channels are actually present. If
+%   full is false, a single channel is sufficient to return true on whether
+%   the run exists. Eventtimes and other files are never checked. Default:
+%   false
+%
+%   exist(PP, 'directMode', true) checks if there is a run in the _direct
+%   folder, meaning it has just completed by the neuron program. Default:
+%   false
+
+    % Handle input
+    p = inputParser();
+    p.addParamValue('full', false);
+    p.addParamValue('directMode', false);
+    p.parse(varargin{:});
+    p = p.Results;
+    
+    % Direct mode
+    fpath = obj.datadir();
+    irec = obj.irec;
+    if p.directMode
+        fpath = obj.directdir();
+        irec = obj.iter;
+    end
+
+    % The object must be a specific run
+    assert(~isempty(experiment(obj)), 'Please specify an experiment');
+    assert(~isempty(irec), 'Please specify a specific run');
+    
+    % Handle output
+    doesExist = 0;
+    
+    % Variables
+    flExt = obj.sim_channels;
+    
+    % If a single channel exists, consider the run as existing
+    for ii = 1:numel(flExt)
+        % Generate a filename
+        if irec == 0
+            fname = fullfile(fpath, sprintf('run_unsaved__vm_%s.txt', flExt{ii}));
+        else
+            fname = fullfile(fpath, sprintf('run_%05i__vm_%s.txt', irec, flExt{ii}));
+        end
+        
+        % Check if file exists
+        if exist(fname, 'file')
+            doesExist = doesExist + 1;
+        end
+    end
+    
+    % Handle doesExist
+    if p.full
+        doesExist = doesExist >= numel(flExt);
+    else
+        doesExist = doesExist > 1;
+    end
+    
+    % Output
+    output = doesExist;
+    if nargout == 0
+        switch doesExist
+            case 0
+                fprintf('Run %i does not exist\n', irec);
+            case 1
+                fprintf('Run %i exists\n', irec);
+        end
+    end
+
+end
diff --git a/matlab/@pricklypear/experiment.m b/matlab/@pricklypear/experiment.m
new file mode 100644
index 0000000..828cc0c
--- /dev/null
+++ b/matlab/@pricklypear/experiment.m
@@ -0,0 +1,7 @@
+function output = experiment(obj)
+%EXPERIMENT Get experiment of pricklypear object
+%   experiment(PP) returns the experiment of PP
+    
+    output = obj.sim_experiment;
+    
+end
diff --git a/matlab/@pricklypear/irec.m b/matlab/@pricklypear/irec.m
new file mode 100644
index 0000000..27750f7
--- /dev/null
+++ b/matlab/@pricklypear/irec.m
@@ -0,0 +1,7 @@
+function output = irec(obj)
+%IREC Get irec of pricklypear object
+%   irec(PP) returns the irec of PP
+    
+    output = obj.sim_irec;
+    
+end
diff --git a/matlab/@pricklypear/modeldir.m b/matlab/@pricklypear/modeldir.m
new file mode 100644
index 0000000..943c1b9
--- /dev/null
+++ b/matlab/@pricklypear/modeldir.m
@@ -0,0 +1,9 @@
+function output = modeldir(obj)
+%MODELDIR Get the directory where the model is stored
+%   modeldir(PP)
+%
+%   See also pricklypear/datadir, pricklypear/directdir
+
+    output = fullfile(obj.basedir, 'neuron');
+
+end
diff --git a/matlab/@pricklypear/notes.m b/matlab/@pricklypear/notes.m
new file mode 100644
index 0000000..fae4ba9
--- /dev/null
+++ b/matlab/@pricklypear/notes.m
@@ -0,0 +1,54 @@
+function output = notes(obj, varargin)
+%NOTES See/edit notes for a pricklypear experiment or run
+%   notes(PP) for notes of an entire experiment (PP without irec)
+%   notes(PP) for notes of a single run (PP with irec)
+%   notes(..., 'edit') to edit the notes for PP
+
+    % Handle input
+    p = inputParser();
+    p.addOptional('mode', 'edit', @ischar);
+    p.parse(varargin{:});
+    p = p.Results;
+    
+    % The object must be an experiment or a saved run
+    assert(~isempty(experiment(obj)), 'Please specify an experiment');
+    
+    if ~isempty(obj.irec)
+        % The run must be saved
+        assert(obj.irec > 0, 'Only saved runs can have notes');
+        
+        % The run must exist
+        assert(exist(obj) == true, 'The run was not found');
+    else
+        % The experiment must exist
+        %assert(exists(obj) == true, 'The experiment was not found');
+    end
+    
+    % File name
+    if isempty(obj.irec) % Experiment
+        fname = fullfile(datadir(obj), 'experiment_notes.txt');
+    else % Run
+        fname = fullfile(datadir(obj), sprintf('run_%05i__notes.txt', obj.irec));
+    end
+    
+    % Mode
+    switch p.mode
+        case {'view' 'read'}
+            if ~exist(fname, 'file')
+                warning('pricklypear:noNotes', 'No notes found for this experiment/run');
+                return;
+            end
+            
+            txt = fileread(fname);
+            if nargout == 0
+                fprintf('\n%s\n\n', txt);
+            else
+                output = txt;
+            end
+        case {'edit' 'write'}
+            edit(fname);
+        otherwise
+            error('The provided action was not recognized (must be view/edit)');
+    end
+
+end
diff --git a/matlab/@pricklypear/parent.m b/matlab/@pricklypear/parent.m
new file mode 100644
index 0000000..642bd66
--- /dev/null
+++ b/matlab/@pricklypear/parent.m
@@ -0,0 +1,11 @@
+function output = parent(obj)
+%PARENT Spawn the parent of a pricklypear object
+%   P.irec -> 2
+%   PP = parent(P)
+%   PP.irec -> []
+%
+%   See also pricklypear/child
+
+    output = pricklypear(experiment(obj));
+
+end
diff --git a/matlab/@pricklypear/plot.m b/matlab/@pricklypear/plot.m
new file mode 100644
index 0000000..8eb3fd0
--- /dev/null
+++ b/matlab/@pricklypear/plot.m
@@ -0,0 +1,83 @@
+function plot(obj, varargin)
+%PLOT Plot the data of a PricklyPear object
+%   plot(PP)
+%
+%   See also pricklypear/anadata
+
+    % Handle input
+    p = inputParser();
+    p.addOptional('channel', '', @(x)ischar(x)||iscell(x));
+    p.addOptional('metric', '', @(x)ischar(x)||iscell(x));
+    p.parse(varargin{:});
+    p = p.Results;
+    
+    if ~isempty(p.channel) && ischar(p.channel)
+        p.channel = {p.channel};
+    end
+    if ~isempty(p.metric) && ischar(p.metric)
+        p.metric = {p.metric};
+    end
+    
+    % Plot
+    figureFullScreen(1000, 800);
+    hold on;
+
+    clrs = {'green_500' 'green_700' 'green_900'...
+            'blue_700'...
+            'green_900' 'green_700' 'green_500'...
+            'red_700' 'red_500'};
+    
+    % Select channels
+    channels = {'d0d' 'd0m' 'd0p' 'c' 'd1p' 'd1m' 'd1d' 'ap' 'ad'};
+    if ~isempty(p.channel)
+        [channels iChan] = intersect(channels, p.channel);
+        clrs = clrs(iChan);
+    end
+    
+    % Spiketimes
+    spkt = spiketimes(obj);
+    
+    % SubplotXY options
+    opts = struct();
+    opts.margin = [0.1 0.08 0.08 0.08];
+    opts.spaceBetweenPlots = 0;
+    opts.nX = 1;
+    opts.nY = numel(channels);
+    opts.iX = 1;
+    opts.iY = 1;
+    
+    % Other variables
+    h = zeros(numel(channels), 1);
+    yl = [-80 60];
+    
+    % For each channel
+    for ii = 1:numel(channels)
+        opts.iY = ii;
+        h(ii) = subplotXY(opts);
+        hold on;
+        
+        % Plot spiketimes
+        if numel(spkt) > 0
+            plot([spkt spkt zeros(size(spkt))*nan]',...
+                 [ones(size(spkt))*yl(1) ones(size(spkt))*yl(2) ones(size(spkt))*yl(2)]',...
+                 'Color', [0.75 0.75 0.75]);
+        end
+        
+        % Plot data
+        plot(obj.x, anadata(obj, channels{ii})+(ii-1), 'Color', mdc(clrs{ii}));
+        box off;
+        ylim(yl);
+        ylabel(channels{ii});
+        
+        if mod(ii,2) == 0
+            set(gca, 'Color', [1 1 1]*0.95);
+        end
+        
+        if ii < numel(channels)
+            set(gca, 'xtick', []);
+        end
+    end
+    
+    linkaxes(h, 'x');
+
+end
diff --git a/matlab/@pricklypear/pricklypear.m b/matlab/@pricklypear/pricklypear.m
new file mode 100644
index 0000000..2b83db9
--- /dev/null
+++ b/matlab/@pricklypear/pricklypear.m
@@ -0,0 +1,146 @@
+classdef pricklypear < handle
+    %PRICKYLPEAR Interface to the MSO neuron model from matlab
+    %   PP = pricklypear() creates a new MSO neuron model interface.
+    %   PP = pricklypear(expName) specifies the experiment name.
+    %   PP = pricklypear(expName, irec) specifies an irec.
+    %
+    %   See also pricklypear/run
+
+    properties
+        seed = 1                % RNG seed
+        iter = 99               % Used by model (do not touch!)
+
+        tstop = 2000            % Length of run
+        dt = 0.025              % Temporal resolution
+
+        v_init = -65            % Initial membrane potential
+        celsius = 38            % Temperature (celsius)
+
+        axon_soma_distance = 45 % Axon soma distance (on ipsi dendrite)
+        dend_n = 2              % Number of dendrites
+        dend_n_syn = 10         % Number of synapses per dendrite
+        dend_n_seg = 20         % Number of segments per dendrite
+        dend_length = 200       % Length of dendrites (um)
+        dend_syn_spread = 0.5   % Spread of synapses (index)
+        dend_syn_offset = 0.45  % Offset of synapses (index)
+
+        dend_0_exc_G = 11       % Excitatory conductance contra dendrite
+        dend_1_exc_G = 11       % Excitatory conductance ipsi dendrite
+
+        dend_0_inh_G = 0        % Inhibitory conductance contra dendrite
+        dend_1_inh_G = 0        % Inhibitory conductance ipsi dendrite
+
+        dend_0_exc_gain = 0.0008    % Excitatory gain contra dendrite
+        dend_1_exc_gain = 0.0008    % Excitatory gain ipsi dendrite
+
+        dend_0_inh_gain = 0.001     % Inhibitory gain contra dendrite
+        dend_1_inh_gain = 0.001     % Inhibitory gain ipsi dendrite
+
+        spk_thres_exc_contra    = 0 % Excitatory spike threshold contra
+        spk_thres_exc_ipsi      = 0	% Excitatory spike threshold ipsi
+        spk_thres_inh_contra    = 1	% Inhibitory spike threshold contra
+        spk_thres_inh_ipsi      = 1	% Inhibitory spike threshold ipsi
+
+        eventtimes_contra   = 10+(0:75:1500)	% Input spiketimes contra (ms)
+        eventtimes_ipsi     = 10+(0:100:1500)	% Input spiketimes ipsi (ms)
+        eventtimes_inh_contra   = []            % Inh input spiketimes contra (ms)
+        eventtimes_inh_ipsi     = []            % Inh input spiketimes ipsi (ms)
+        model_name          = 'default'         % Model name
+        user_data           = []                % Arbitrary user data
+        
+        % Private
+        sim_channels = {'c' 'd0p' 'd0m' 'd0d' 'd1p' 'd1m' 'd1d' 'ap' 'ad'}
+        sim_metrics = {'vm'}
+        sim_experiment = ''
+        sim_irec = []
+        template_exclude = {'eventtimes_contra' 'eventtimes_ipsi'...
+                            'eventtimes_inh_contra' 'eventtimes_inh_ipsi'...
+                            'template_exclude' 'model_name', 'user_data'...
+                            'sim_irec' 'sim_experiment' 'sim_channels' 'sim_metrics'}
+    end
+    
+    methods
+        % Class constructor
+        function obj = pricklypear(varargin)
+            % Handle input
+            p = inputParser();
+            p.addOptional('expname', '', @ischar);
+            p.addOptional('irec', []);
+            p.parse(varargin{:});
+            p = p.Results;
+            
+            % Check expname
+            assert(~strcmp(p.expname, '_direct'), '_direct is a protected experiment name, pick any other.');
+            
+            % Handle expname
+            if ~isempty(p.expname)
+                obj.sim_experiment = p.expname;
+            end
+            
+            % Handle irec
+            if ~isempty(p.expname) && ~isempty(p.irec)
+                obj.sim_irec = p.irec;
+                
+                % Check if this run exists
+                if exist(obj)
+                    % Get properties from old runs
+                    % Compatibility issues will likely arise here...
+                    props = properties(obj);
+                    if p.irec == 0
+                        storedProps = load(fullfile(obj.datadir(), sprintf('run_unsaved__properties.ppbin')), '-mat');
+                    else
+                        storedProps = load(fullfile(obj.datadir(), sprintf('run_%05i__properties.ppbin', p.irec)), '-mat');
+                    end
+                    
+                    for ii = 1:numel(props)
+                        prop = props{ii};
+                        prop_ext = prop;
+                        % Skip irec as it has already been dealt with
+                        if isequal(prop_ext, 'sim_irec')
+                            continue;
+                        end
+                        %%% Backwards compatibility
+                        % Skip a property if not stored
+                        if ~isfield(storedProps, prop)
+                            continue;
+                        end
+                        % Renamed spiketimes to eventtimes
+                        if ~isempty(strfind(prop_ext, 'spiketimes'))
+                            prop_ext = strrep(prop_ext, 'spiketimes', 'eventtimes');
+                        end
+                        % Overwrite the field
+                        obj.(prop) = storedProps.(prop_ext);
+                    end
+                else
+                    %error('Run %i does not exist', p.irec);
+                end
+                
+            end
+            
+        end
+        
+        disp(obj)
+        plot(obj, varargin)
+        obj = save(obj)
+        obj = run(obj)
+        obj = irec(obj)
+        obj = experiment(obj)
+        output = struct(obj)
+        output = exist(obj, varargin)
+        output = delete(obj, varargin)
+        output = anadata(obj, varargin)
+        output = spiketimes(obj, varargin)
+        output = notes(obj, varargin)
+        output = strpad(obj, varargin)
+        output = basedir(obj)
+        output = datadir(obj)
+        output = directdir(obj)
+        output = x(obj)
+        output = subsref(obj, irec)
+        output = parent(obj)
+        output = child(obj, irec)
+        output = runSubthresholdModel(obj, varargin)
+    end
+    
+end
+
diff --git a/matlab/@pricklypear/run.m b/matlab/@pricklypear/run.m
new file mode 100644
index 0000000..8f925a7
--- /dev/null
+++ b/matlab/@pricklypear/run.m
@@ -0,0 +1,91 @@
+function obj = run(obj)
+%RUN Run PricklyPear model simulation
+%   run(PP)
+
+    % Variables
+    cdir = pwd();
+    
+    % Read the base hoc file
+    fname = fullfile(obj.modeldir(), sprintf('msomodel_%s.hoc', obj.model_name));
+    fid = fopen(fname);
+    txt = textscan(fid, '%s', 'delimiter', '\n');
+    txt = txt{1};
+    fclose(fid);
+    
+    % Replace the main placeholder
+    phLine = find(strcmp('//%%%TEMPLATE%%%', txt));
+    
+    props = properties(obj);
+    phNew = cell(numel(props), 1);
+    for ii = 1:numel(props)
+        if any(strcmp(obj.template_exclude, props{ii}))
+            continue;
+        end
+        phNew{ii} = sprintf('%s = %g', props{ii}, obj.(props{ii}));
+    end
+    phNew = [{'// PricklyPear parameters'};...
+             phNew;...
+             {'// If you see this file, something went wrong during pricklypear.run(). Debug and remove this file'}];
+    
+    txt = [txt(1:phLine-1); phNew; txt(phLine+1:end)];
+    
+    % Replace the ending placeholder
+    phLine = strcmp('//%%%TEMPLATE_ENDING%%%', txt);
+    txt{phLine} = 'quit()';
+    
+    % Write the temporary hoc file
+    fname = fullfile(obj.modeldir(), 'msomodel__temp.hoc');
+    fid = fopen(fname, 'w');
+    fprintf(fid, '%s\n', txt{:});
+    fclose(fid);
+    
+    % Write the input event times
+    fname = fullfile(obj.modeldir(), 'eventtimes_contra.txt');
+    fid = fopen(fname, 'w');
+    fprintf(fid, '%i\r\n', obj.eventtimes_contra);
+    fclose(fid);
+    fname = fullfile(obj.modeldir(), 'eventtimes_ipsi.txt');
+    fid = fopen(fname, 'w');
+    fprintf(fid, '%i\n', obj.eventtimes_ipsi);
+    fclose(fid);
+    
+    % Remove existing data files
+    if exist(obj, 'directMode', true)
+        isDeleted = delete(obj, 'directMode', true);
+    end
+    parobj = parent(obj);
+    if exist(child(parobj, 0))
+        isDeleted = delete(child(parobj, 0));
+    end
+    
+    % Execute
+    cd(obj.modeldir());
+    system('neuron -nogui msomodel__temp.hoc');
+    cd(cdir);
+    
+    % Wait for simulation to finish
+    while ~exist(obj, 'directMode', true, 'full', true)
+        pause(4);
+    end
+    
+    % Move data
+    if ~isempty(obj.sim_experiment)
+        [unused unused] = mkdir(obj.datadir());
+        for ii = 1:numel(obj.sim_channels)
+            movefile(fullfile(obj.directdir(), sprintf('run_%05i__vm_%s.txt', obj.iter, obj.sim_channels{ii})),...
+                     fullfile(obj.datadir(), sprintf('run_unsaved__vm_%s.txt', obj.sim_channels{ii})));
+        end
+    end
+    copyfile(fullfile(obj.modeldir(), sprintf('eventtimes_contra.txt')),...
+             fullfile(obj.datadir(),  sprintf('run_unsaved__evt_contra.txt')));
+    copyfile(fullfile(obj.modeldir(), sprintf('eventtimes_ipsi.txt')),...
+             fullfile(obj.datadir(),  sprintf('run_unsaved__evt_ipsi.txt')));
+	props = struct(obj);
+    save(fullfile(obj.datadir(), 'run_unsaved__properties.ppbin'), '-struct', 'props');
+    
+    % Clean up
+    delete(fullfile(obj.modeldir(), 'msomodel__temp.hoc'));
+    delete(fullfile(obj.modeldir(), 'eventtimes_contra.txt'));
+    delete(fullfile(obj.modeldir(), 'eventtimes_ipsi.txt'));
+
+end
diff --git a/matlab/@pricklypear/runSubthresholdModel.m b/matlab/@pricklypear/runSubthresholdModel.m
new file mode 100644
index 0000000..4dd7b15
--- /dev/null
+++ b/matlab/@pricklypear/runSubthresholdModel.m
@@ -0,0 +1,83 @@
+function output = runSubthresholdModel(obj, S, varargin)
+%RUNSUBTHRESHOLDMODEL Run a simulation as defined by SubthresholdModel
+%   O = runSubthresholdModel(PP, S) uses the output of subthrmodel to run 
+%   new simulations mimicking that model as best as possible. The output O
+%   contains six PricklyPear instances corresponding to the 6 DZW
+%   conditions as well as all the spiketimes for each conditions.
+%
+%   runSubthresholdModel(..., 'axon_soma_distance', 45) sets the distance
+%   between axon and soma (on ipsi dendrite) in microns. Default: 45
+
+    % TODO Nicely implement this quickfix
+    add_delay = 200; % ms
+
+    % Handle input
+    p = inputParser();
+    p.addParamValue('axon_soma_distance', 45);
+    p.parse(varargin{:});
+    p = p.Results;
+
+    % Handle output
+    output = struct;
+    output.irec = [];
+    output.spt = struct;
+    
+    % The object must be an experiment (not a run)
+    assert(~isempty(experiment(obj)), 'Please specify an experiment');
+    assert(isempty(irec(obj)), 'Do not specify a run (irec)');
+    
+    % Log
+    fprintf('Starting runSubthresholdModel...\n');
+    
+    % For each of the conditions
+    for iCond = 1:numel(S.AllCond)
+        % Get the cond
+        cond = S.AllCond{iCond};
+        
+        % Log
+        fprintf('Condition %i: %s\n', iCond, cond);
+        fprintf('  - Setting parameters...\n');
+        
+        % Set PP parameters to S parameters
+        obj.axon_soma_distance = p.axon_soma_distance;
+        obj.eventtimes_contra = S.ev.(cond).C.evt + add_delay;
+        obj.eventtimes_ipsi   = S.ev.(cond).I.evt + add_delay;
+        obj.tstop = S.Param.Dur;
+        obj.user_data = struct;
+        obj.user_data.subthr = S;
+        obj.user_data.cond = cond;
+        
+        % Log
+        fprintf('  - Parameters set!\n');
+        fprintf('  - Simulation duration: %i seconds.\n', ceil(obj.tstop/1e3));
+        fprintf('  - Expected processing time: %i seconds.\n', ceil(obj.tstop/1e3)*5);
+        fprintf('  - Running simulation...\n');
+        
+        % Run simulation
+        run(obj);
+        
+        % Log
+        fprintf('  - Run finished!\n');
+        fprintf('  - Saving run...\n');
+        
+        % Save run
+        nextIrec = save(child(obj, 0));
+        output.irec(iCond) = nextIrec;
+        
+        % Log
+        fprintf('  - Run saved as run #%i!\n', nextIrec);
+        fprintf('  - Putting spiketimes in output...\n');
+        
+        % Put spiketimes in output struct
+        output.spt.(cond) = struct;
+        output.spt.(cond).spt = spiketimes(child(obj, nextIrec));
+        
+        % Log
+        fprintf('  - Done!\n');
+        
+    end
+    
+    % Log
+    fprintf('runSubthresholdModel is done!\n');
+
+end
diff --git a/matlab/@pricklypear/save.m b/matlab/@pricklypear/save.m
new file mode 100644
index 0000000..347098a
--- /dev/null
+++ b/matlab/@pricklypear/save.m
@@ -0,0 +1,44 @@
+function output = save(obj)
+%SAVE Save a PricklyPear run
+%   save(PP(0)) stores the unsaved run from experiment(PP) to the next
+%   available irec. Storing already saved runs is not allowed.
+
+    % The object must be the unsaved simulation (irec = 0)
+    assert(~isempty(experiment(obj)), 'Please specify an experiment');
+    assert(~isempty(obj.irec) && obj.irec == 0, 'Only the unsaved simulation (irec = 0) can be saved');
+    
+    % The unsaved experiment must exist
+    assert(exist(obj) == true, 'No unsaved simulation found');
+    
+    % What is the next available irec?
+    parobj = parent(obj);
+    nextIrec = [];
+    ii = 1;
+    while isempty(nextIrec)
+        if ~exist(child(parobj, ii), 'full', true)
+            nextIrec = ii;
+        else
+            ii = ii + 1;
+        end
+    end
+    
+    % Move data
+    [unused unused] = mkdir(obj.datadir());
+    for ii = 1:numel(obj.sim_channels)
+        movefile(fullfile(obj.datadir(), sprintf('run_unsaved__vm_%s.txt', obj.sim_channels{ii})),...
+                 fullfile(obj.datadir(), sprintf('run_%05i__vm_%s.txt', nextIrec, obj.sim_channels{ii})));
+    end
+    movefile(fullfile(obj.datadir(), sprintf('run_unsaved__evt_contra.txt')),...
+             fullfile(obj.datadir(), sprintf('run_%05i__evt_contra.txt', nextIrec)));
+    movefile(fullfile(obj.datadir(), sprintf('run_unsaved__evt_ipsi.txt')),...
+             fullfile(obj.datadir(), sprintf('run_%05i__evt_ipsi.txt', nextIrec)));
+    movefile(fullfile(obj.datadir(), sprintf('run_unsaved__properties.ppbin')),...
+             fullfile(obj.datadir(), sprintf('run_%05i__properties.ppbin', nextIrec)));
+    
+    % Notify the user or output the irec
+    output = nextIrec;
+    if nargout == 0
+        fprintf('Run saved as run %i\n', nextIrec);
+    end
+    
+end
diff --git a/matlab/@pricklypear/spiketimes.m b/matlab/@pricklypear/spiketimes.m
new file mode 100644
index 0000000..1bfec21
--- /dev/null
+++ b/matlab/@pricklypear/spiketimes.m
@@ -0,0 +1,25 @@
+function output = spiketimes(obj, varargin)
+%SPIKETIMES Spiketimes from pricklypear simulation
+%   spiketimes(PP)
+%   spiketimes(PP, thr)
+%   spiketimes(PP, thr, chan)
+
+    % Handle input
+    p = inputParser();
+    p.addOptional('threshold', 0);
+    p.addOptional('channel', 'ap', @ischar);
+    p.parse(varargin{:});
+    p = p.Results;
+    
+    % The object must be a specific run
+    assert(~isempty(experiment(obj)), 'Please specify an experiment');
+    assert(~isempty(irec(obj)), 'Please specify a specific run');
+    
+    % Get the raw recording
+    ad = anadata(obj, p.channel, 'vm');
+    
+    % Detect peaks
+    output = peakpicker(obj.dt, ad, p.threshold);
+    output = output{1};
+
+end
diff --git a/matlab/@pricklypear/strpad.m b/matlab/@pricklypear/strpad.m
new file mode 100644
index 0000000..1a85e0b
--- /dev/null
+++ b/matlab/@pricklypear/strpad.m
@@ -0,0 +1,60 @@
+function S = strpad(obj, S, varargin)
+%STRPAD String rightpadding
+%   S = strpad(CA) applies rightpadding of spaces to all elements of the
+%   cell array CA so that all have the same length and returns the cell
+%   array S.
+%
+%   S = strpad(CA, n) sets extra rules for the padding. If n = 0 (default),
+%   the longest string determines the length. If n > 0, n will be the
+%   length until which padding is added. Note that if a string is longer
+%   than n, it will not be truncated. If n < 0, -n is a margin, a number of
+%   spaces added to the longest string and that then decides the final
+%   length for all elements.
+%
+%   S = strpad(CA, n, '-truncate') truncates string if their length
+%   exceeds the value of n > 0 (see above).
+
+    % Handle flags
+    [flg varargin] = flagParser(varargin, 'truncate');
+
+    % Handle input
+    p = inputParser();
+    p.addRequired('S');
+    p.addOptional('n', 0);
+    p.parse(S, varargin{:});
+    p = p.Results;
+    
+    S = p.S;
+    n = p.n;
+    
+    % If S is an array, make it a cell array
+    if isnumeric(S)
+        S = cellfun(@num2str, num2cell(S), 'UniformOutput', false);
+    end
+    
+    % If S is a string, make it a cell array
+    if ischar(S)
+        S = {S};
+    end
+    
+    % Padding method depends on n
+    if n > 0
+        maxLength = p.n;
+    else
+        maxLength = max(max(cellfun(@length, S))) - n;
+    end
+    
+    % For each element, add padding
+    S = cellfun(@(s)([s char(zeros(1, maxLength-numel(s))+32)]), S, 'UniformOutput', false);
+    
+    % Truncation
+    if flg.truncate
+        S = cellfun(@(s)(s(1:maxLength)), S, 'UniformOutput', false);
+    end
+    
+    % Display mode TODO
+    if nargout > 0
+        return;
+    end
+
+end
diff --git a/matlab/@pricklypear/struct.m b/matlab/@pricklypear/struct.m
new file mode 100644
index 0000000..9cf5571
--- /dev/null
+++ b/matlab/@pricklypear/struct.m
@@ -0,0 +1,24 @@
+function output = struct(obj)
+%STRUCT Generate a struct from a PricklyPear object
+%   struct(PP)
+
+%     % Handle input
+%     p = inputParser();
+%     p.addOptional('mode', 'public', @ischar);
+%     p.parse(varargin{:});
+%     p = p.Results;
+
+    % Handle output
+    output = struct;
+    
+    % Handle the properties
+    P = properties(obj);
+    
+    % Transfer properties to output
+    for ii = 1:numel(P)
+        p = P{ii};
+        output.(p) = obj.(p);
+    end
+
+end
+
diff --git a/matlab/@pricklypear/subsref.m b/matlab/@pricklypear/subsref.m
new file mode 100644
index 0000000..f1c2904
--- /dev/null
+++ b/matlab/@pricklypear/subsref.m
@@ -0,0 +1,20 @@
+function output = subsref(obj, S)
+%SUBSREF Handle the subreffing of PricklyPear objects
+%   Not to be called directly
+
+    output = obj;
+
+    for ii = 1:numel(S)
+        s = S(ii);
+        
+        switch s.type
+            case '()'
+                if numel(s.subs) > 0
+                    output = obj.child(s.subs{1});
+                end
+            case '.'
+                output = builtin('subsref', output, s);
+        end
+    end
+
+end
diff --git a/matlab/@pricklypear/x.m b/matlab/@pricklypear/x.m
new file mode 100644
index 0000000..69fa188
--- /dev/null
+++ b/matlab/@pricklypear/x.m
@@ -0,0 +1,7 @@
+function output = x(obj)
+%X Generate a simulation timeline (used as x axis for plotting)
+%   x(PP)
+
+    output = (0:obj.dt:obj.tstop)';
+
+end
diff --git a/matlab/flagParser.m b/matlab/flagParser.m
new file mode 100644
index 0000000..84993f2
--- /dev/null
+++ b/matlab/flagParser.m
@@ -0,0 +1,31 @@
+function [S args] = flagParser(args, flags)
+%FLAGPARSER Parse varargin for flags (BEFORE INPUTPARSER)
+%   [S args] = flagParser(args, flags) parses the cell array args,
+%   checking for any flags provided in the cell array flags, then returns
+%   the args without the flags so they can be processed in inputParser and
+%   the struct S, containing all the booleans specifying the absence or
+%   presence of the flags in the args.
+%
+%   Please provide the flags without the leading dash.
+    
+    % Convert input to cell array if needed
+    if ischar(flags)
+        flags = {flags};
+    end
+    
+    % tempFlags has a dash prepended to all flags
+    tempFlags = cellfun(@(x)(['-' x]), flags, 'UniformOutput', false);
+    
+    % tempArgs contains only strings
+    tempArgs = cellfun(@num2str, args, 'UniformOutput', false);
+
+    % Check for every potential flag
+    [hasFlags idxFlag] = ismember(tempFlags, tempArgs);
+    
+    % Make the output struct S
+    S = cell2struct(num2cell(hasFlags), flags, 2);
+    
+    % Remove the flags from the args (compatible with inputParser)
+    args(idxFlag(hasFlags)) = [];
+
+end
diff --git a/matlab/pp_dir.template.txt b/matlab/pp_dir.template.txt
new file mode 100644
index 0000000..86180d9
--- /dev/null
+++ b/matlab/pp_dir.template.txt
@@ -0,0 +1 @@
+C:\matlab\pricklypear\
\ No newline at end of file