classdef WholeCellRecording
    properties
        %% Acquisition:
        fs

        %% Database:
        filename
        paradigms
        response_durations
        response_samples

        %% Stimulus:
        rates
        npulses

        %% Recording:
        times
        stimulus
        membrane_potential
        injected_current
        representative_response

        %% Cell:
        membrane_capacitance
        input_resistance
        resting_potential
        excitatory_reversal_potential 
        inhibitory_reversal_potential
        threshold_potential
        activation_potential
        steady_state_potential
        alpha_multiplier
        beta_multiplier
        spikes_per_stimulus
        reference_potential

        %% Analysis:
        membrane_current
        leakage_current
        alpha
        beta
        activation_current
        excitatory_conductance
        inhibitory_conductance
        excitatory_current
        inhibitory_current
        depolarizations
        hyperpolarizations
        excitation
        inhibition
        resultant_excitaiton
        resultant_inhibition
        estimated_membrane_potential
        estimated_time
        estimated_activation_current
        estimated_excitatory_current
        estimated_inhibitory_current
        estimated_leakage_current
        threshold_crossing
        estimated_membrane_current

        %% Filters
        lowpass
        notch

        %% Stats
        means
        maxima
        minima

        %% Plots
        estimations_figure
        stats_figure
    end

    methods(Access=public)
        %% Datastructure building methods
        function app = WholeCellRecording(filename, paradigms, response_durations, filter_parameters)
            if nargin > 0
                [m, n] = size(paradigms);
                app(m, n) = app;
                for i = 1: m
                    for j = 1: n
                        app(i, j).filename = filename;
                        app(i, j).paradigms = paradigms(i, j);
                        app(i, j).response_durations = response_durations(i, j);
                    end
                end
                app = app.read_data();
                app = app.read_parameters();
                app = app.build_filters(filter_parameters);
                app = app.build_analysis_arrays();
            end
        end

        function app = read_data(app)
            [m, n] = size(app);
            for i = 1: m
                for j = 1: n
                    data = readtable(app(i, j).filename, 'Sheet', app(i, j).paradigms, 'ReadVariableNames', 1, 'VariableNamingRule', 'preserve');
                    try
                        app(i, j).times = data.times;
                        app(i, j).fs = 1/(app(i, j).times(2, 1) - app(i, j).times(1, 1));
                        data.times = [];
                    catch
                        error(strcat('times array was not found in ', app(i, j).paradigms, '!'));
                    end
                    try
                        app(i, j).stimulus = data.stimulus;
                        data.stimulus = [];
                    catch
                        warning(strcat("stimulus array was not found in ", app(i, j).paradigms, "!"));
                        app(i, j).stimulus = [];
                    end
                    try
                        app(i, j).representative_response = data.representative .* 1e-3;
                        data.representative = [];
                    catch
                        warning(strcat("Representative array was not found in ", app(i, j).paradigms, "!"));
                        app(i, j).representative_response = [];
                    end
                    try
                        app(i, j).membrane_potential = table2array(data).*1e-3;
                    catch
                        error(strcat('membrane potentials was not found in ', app(i, j).paradigms, '!'));
                    end
                    app(i, j).response_samples = floor(app(i, j).fs * app(i, j).response_durations);
                    if app(i, j).response_samples > size(app(i, j).membrane_potential, 1)
                        app(i, j).response_samples = size(app(i, j).membrane_potential, 1);
                    end
                end
            end
        end

        function app = read_parameters(app)
            [m, n] = size(app);
            for i = 1: m
                for j = 1: n
                    parameters = readtable(app(i, j).filename, 'Sheet', strcat('parameters_', app(i, j).paradigms));
                    app(i, j).injected_current = parameters.Iinj' + zeros(size(app(i, j).membrane_potential));
                    app(i, j).membrane_capacitance = parameters.Cm' + zeros(size(app(i, j).membrane_potential));
                    app(i, j).input_resistance = parameters.Rin' + zeros(size(app(i, j).membrane_potential));
                    app(i, j).resting_potential = parameters.Er' + zeros(size(app(i, j).membrane_potential));
                    app(i, j).excitatory_reversal_potential = parameters.Ee' + zeros(size(app(i, j).membrane_potential));
                    app(i, j).inhibitory_reversal_potential = parameters.Ei' + zeros(size(app(i, j).membrane_potential));
                    app(i, j).threshold_potential = parameters.Et' + zeros(size(app(i, j).membrane_potential));
                    app(i, j).activation_potential = parameters.Eact' + zeros(size(app(i, j).membrane_potential));
                    app(i, j).steady_state_potential = parameters.Ess' + zeros(size(app(i, j).membrane_potential));
                    app(i, j).alpha_multiplier = parameters.xalpha' + zeros(size(app(i, j).membrane_potential));
                    app(i, j).beta_multiplier = parameters.xbeta' + zeros(size(app(i, j).membrane_potential));
                    app(i, j).spikes_per_stimulus = parameters.sps' + zeros(size(app(i, j).membrane_potential));
                    app(i, j).reference_potential = parameters.Eref' + zeros(size(app(i, j).membrane_potential));
                    app(i, j).rates = parameters.rate';
                    app(i, j).npulses = parameters.npulses';
                    app(i, j).threshold_crossing = 0;
                end
            end
        end

        function app = build_analysis_arrays(app)
            [m, n] = size(app);
            for i = 1: m
                for j = 1: n
                    app(i, j).membrane_current = zeros(size(app(i, j).membrane_potential));
                    app(i, j).leakage_current = zeros(size(app(i, j).membrane_potential));
                    app(i, j).alpha = zeros(size(app(i, j).membrane_potential));
                    app(i, j).beta = zeros(size(app(i, j).membrane_potential));
                    app(i, j).activation_current = zeros(size(app(i, j).membrane_potential));
                    app(i, j).excitatory_conductance = zeros(size(app(i, j).membrane_potential, 1), 1);
                    app(i, j).inhibitory_conductance = zeros(size(app(i, j).membrane_potential, 1), 1);
                    app(i, j).excitatory_current = zeros(size(app(i, j).membrane_potential));
                    app(i, j).inhibitory_current = zeros(size(app(i, j).membrane_potential));
                    app(i, j).depolarizations = zeros(size(app(i, j).membrane_potential));
                    app(i, j).hyperpolarizations = zeros(size(app(i, j).membrane_potential));
                    app(i, j).excitation = zeros(size(app(i, j).membrane_potential, 1), 1);
                    app(i, j).inhibition = zeros(size(app(i, j).membrane_potential, 1), 1);
                    app(i, j).resultant_excitaiton = zeros(size(app(i, j).membrane_potential, 1), 1);
                    app(i, j).resultant_inhibition = zeros(size(app(i, j).membrane_potential, 1), 1);
                    app(i, j).estimated_membrane_potential = app(i, j).resting_potential(:, 1);
                    app(i, j).estimated_time = zeros(size(app(i, j).membrane_potential, 1), 1);
                    app(i, j).estimated_activation_current = app(i, j).activation_current.*0;
                    app(i, j).estimated_excitatory_current = app(i, j).activation_current.*0;
                    app(i, j).estimated_inhibitory_current = app(i, j).activation_current.*0;
                    app(i, j).estimated_leakage_current = app(i, j).activation_current.*0;
                    app(i, j).estimated_membrane_current = app(i, j).membrane_potential.*0;
                end
            end
        end

        function app = adjust_membrane_potential_wrt_steady_state(app)
            [m, n] = size(app);
            for i = 1: m
                for j = 1: n
                    app(i, j).injected_current = (1./app(i, j).input_resistance).*(app(i, j).steady_state_potential - app(i, j).resting_potential);
                    app(i, j).membrane_potential = app(i, j).membrane_potential - app(i, j).membrane_potential(1, :) + app(i, j).steady_state_potential;
                end
            end
        end

    end
    %% Filter methods
    methods
        function app = build_notch(app, params)
            nyquist_frequency = app(1, 1).fs/2;
            normalized_stop_frequency = params.Fst/nyquist_frequency;
            normalized_stopband_width = params.BW/nyquist_frequency;
            [app(1, 1).notch.num, app(1, 1).notch.den] = iirnotch(normalized_stop_frequency, normalized_stopband_width);
        end
        
        function app = build_lowpass(app, params)
            app(1, 1).lowpass = designfilt("lowpassiir", ...
                "PassbandFrequency", params.Fp, ...
                "StopbandFrequency", params.Fst, ...
                "PassbandRipple", params.Ap, ...
                "StopbandAttenuation", params.Ast, ...
                "SampleRate", app(1, 1).fs);
        end

        function app = build_filters(app, params)
            app = app.build_lowpass(params.lowpass);
            app = app.build_notch(params.notch);
        end

        function app = filter_membrane_potential(app)
            [m, n] = size(app);
            for i = 1: m
                for j = 1: n
                    vm = app(i, j).membrane_potential - app(i, j).membrane_potential(1, :);
                    appender = zeros(size(app(i, j).membrane_potential));
                    vm = cat(1, appender+vm(1, :), vm, appender+vm(end, :));
                    vm = filtfilt(app(1, 1).lowpass, vm);
                    vm = filtfilt(app(1, 1).notch.num, app(1, 1).notch.den, vm);
                    app(i, j).membrane_potential = vm(size(app(i, j).membrane_potential, 1)+1: end-size(app(i, j).membrane_potential, 1), :) + app(i, j).membrane_potential(1, :);
                end
            end
        end

        function app = filter_membrane_current(app)
            [m, n] = size(app);
            lowpass_im = designfilt("lowpassiir", ...
                "PassbandFrequency", app(1, 1).lowpass.PassbandFrequency/5, ...
                "StopbandFrequency", app(1, 1).lowpass.StopbandFrequency/5, ...
                "PassbandRipple", app(1, 1).lowpass.PassbandRipple, ...
                "StopbandAttenuation", app(1, 1).lowpass.StopbandAttenuation, ...
                "SampleRate", app(1, 1).fs);
            for i = 1: m
                for j = 1: n
                    appender = zeros(size(app(i, j).membrane_current));
                    im = app(i, j).membrane_current - app(i, j).membrane_current(1, :);
                    im = cat(1, appender+im(1, :), im, appender+im(end, :));
                    im = filtfilt(lowpass_im, im);
                    im = im(size(appender, 1)+1: end-size(appender, 1), :);
                    im = im - im(1, :);
                    app(i, j).membrane_current = im; % + app(i, j).membrane_current(1, :);
                end
            end
        end

        function app = filter_activation_currents(app)
            [m, n] = size(app);
            lowpass_iact = designfilt("lowpassiir", ...
                "PassbandFrequency", app(1, 1).lowpass.PassbandFrequency/5, ...
                "StopbandFrequency", app(1, 1).lowpass.StopbandFrequency/5, ...
                "PassbandRipple", app(1, 1).lowpass.PassbandRipple, ...
                "StopbandAttenuation", app(1, 1).lowpass.StopbandAttenuation, ...
                "SampleRate", app(1, 1).fs);
            for i = 1: m
                for j = 1: n
                    appender = zeros(size(app(i, j).activation_current));
                    iact = app(i, j).activation_current - app(i, j).activation_current(1, :);
                    iact = cat(1, appender+iact(1, :), iact, appender+iact(end, :));
                    iact = filtfilt(lowpass_iact, iact);
%                     iact = filtfilt(app(1, 1).notch.num, app(1, 1).notch.den, iact);
                    iact = iact(size(appender, 1)+1: end-size(appender, 1), :);
                    iact = iact - iact(1, :);
                    app(i, j).activation_current = iact; % + app(i, j).membrane_current(1, :);
                    % app(i, j).activation_current(app(i, j).activation_current < 0) = 0;
                end
            end
        end

    end

    %% Computing current methods
    methods

        function app = compute_membrane_current(app)
            [m, n] = size(app);
            for i = 1: m
                for j = 1: n
                    vm = app(i, j).membrane_potential - app(i, j).membrane_potential(1, :);
                    appender = zeros(1, size(app(i, j).membrane_potential, 2));
                    dvmdt = diff(cat(1, appender+vm(1, :), vm, appender+vm(end, :)), 1);
                    dvmdt = dvmdt(2:end, :); 
                    im = app(i, j).membrane_capacitance.*dvmdt.*app(i, j).fs;
                    app(i, j).membrane_current = im - im(1, :);
                end
            end
        end

        function app = compute_leakage_current(app)
            [m, n] = size(app);
            for i = 1: m
                for j = 1: n
                    app(i, j).leakage_current = (1./app(i, j).input_resistance).*(app(i, j).membrane_potential - app(i, j).resting_potential);
                end
            end
        end

        function app = compute_active_conductances(app)
            [m, n] = size(app);
            for i = 1: m
                for j = 1: n
                    app(i, j).alpha = ((1./app(i, j).input_resistance)./(2.*(app(i, j).activation_potential-app(i, j).steady_state_potential)));
                    app(i, j).beta = app(i, j).alpha.*(app(i, j).threshold_potential - app(i, j).steady_state_potential);
                    app(i, j).alpha = app(i, j).alpha.*app(i, j).alpha_multiplier;
                    app(i, j).beta = app(i, j).beta.*app(i, j).beta_multiplier;
                end
            end
        end

        function app = compute_activation_currents(app)
            app = app.compute_active_conductances();
            [m, n] = size(app);
            for i = 1: 1: m
               for j = 1: 1: n
                  app(i, j).activation_current = (app(i, j).alpha.*(app(i, j).membrane_potential-app(i, j).threshold_potential).*(app(i, j).membrane_potential-app(i, j).resting_potential)) ...
                          + (app(i, j).beta.*(app(i, j).membrane_potential-app(i, j).resting_potential));
                  app(i, j).activation_current(app(i, j).membrane_potential > app(i, j).activation_potential) = 0;
                  app(i, j).activation_current(app(i, j).membrane_potential < app(i, j).steady_state_potential) = 0;
               end
            end
        end

        function app = compute_passive_conductances(app)
            [m, n] = size(app);
            for i = 1: m
                for j = 1: n
                    samples = size(app(i, j).membrane_potential, 1);
                    A = zeros(2, 2, samples);
                    B = zeros(2, 1, samples);
                    A(1, 1, :) = sum((app(i, j).membrane_potential - app(i, j).excitatory_reversal_potential).^2, 2);
                    A(1, 2, :) = sum((app(i, j).membrane_potential - app(i, j).excitatory_reversal_potential).*(app(i, j).membrane_potential - app(i, j).inhibitory_reversal_potential), 2);
                    A(2, 1, :) = A(1, 2, :);
                    A(2, 2, :) = sum((app(i, j).membrane_potential - app(i, j).inhibitory_reversal_potential).^2, 2);
                    C = app(i, j).membrane_current - app(i, j).injected_current - app(i, j).activation_current + app(i, j).leakage_current;
                    B(1, 1, :) = -sum(C.*(app(i, j).membrane_potential - app(i, j).excitatory_reversal_potential), 2);
                    B(2, 1, :) = -sum(C.*(app(i, j).membrane_potential - app(i, j).inhibitory_reversal_potential), 2);
                    G = pagemtimes(pageinv(A), B);
                    app(i, j).excitatory_conductance = reshape(G(1, 1, :), [samples, 1]);
                    app(i, j).inhibitory_conductance = reshape(G(2, 1, :), [samples, 1]);
                end
            end
        end

        function app = compute_passive_currents(app)
            app = app.compute_passive_conductances();
            [m, n] = size(app);
            for i = 1: m
                for j = 1: n
                    app(i, j).excitatory_current = -1.*app(i, j).excitatory_conductance.*(app(i, j).membrane_potential - app(i, j).excitatory_conductance);
                    app(i, j).inhibitory_current = -1.*app(i, j).inhibitory_conductance.*(app(i, j).membrane_potential - app(i, j).inhibitory_conductance);
                end
            end
        end
    end

    methods

        function app = compute_resultants(app)
            [m, n] = size(app);
            for i = 1: m
                for j = 1: n
                    app(i, j).excitation = (app(i, j).excitatory_conductance > 0).*app(i, j).excitatory_conductance;
                    app(i, j).inhibition = (app(i, j).inhibitory_conductance > 0).*app(i, j).inhibitory_conductance;
                    resultant = app(i, j).excitation - app(i, j).inhibition;
                    app(i, j).resultant_excitaiton = (resultant > 0).*resultant;
                    app(i, j).resultant_inhibition = (resultant < 0).*resultant;
                end
            end
        end

        function app = compute_polarizations(app)
            [m, n] = size(app);
            for i = 1: m
                for j = 1: n
                    polarization = app(i, j).membrane_potential - app(i, j).steady_state_potential;
                    app(i, j).depolarizations = (polarization > 0).*polarization;
                    app(i, j).hyperpolarizations = abs((polarization < 0).*polarization);
                end
            end
        end

        function app = compute_means(app, idx)
            [m, n] = size(app);
            for i = 1: m
                for j = 1: n
                    app(i, j).means.depolarization = mean(app(i, j).depolarizations(1:app(i, j).response_samples, idx), 1);
                    app(i, j).means.hyperpolarization = mean(app(i, j).hyperpolarizations(1:app(i, j).response_samples, idx), 1);
                    app(i, j).means.excitatory_conductance = mean(app(i, j).excitatory_conductance(1:app(i, j).response_samples), 1);
                    app(i, j).means.inhibitory_conductance = mean(app(i, j).inhibitory_conductance(1:app(i, j).response_samples), 1);
                    app(i, j).means.excitatory_current = mean(app(i, j).excitatory_current(1:app(i, j).response_samples, idx), 1);
                    app(i, j).means.inhibitory_current = mean(app(i, j).inhibitory_current(1:app(i, j).response_samples, idx), 1);
                    app(i, j).means.excitation = mean(app(i, j).excitation(1:app(i, j).response_samples), 1);
                    app(i, j).means.inhibition = mean(app(i, j).inhibition(1:app(i, j).response_samples), 1);
                    app(i, j).means.resultant_excitation = mean(app(i, j).resultant_excitaiton(1:app(i, j).response_samples), 1);
                    app(i, j).means.resultant_inhibition = mean(app(i, j).resultant_inhibition(1:app(i, j).response_samples), 1);
                    app(i, j).means.activation_current = mean(app(i, j).activation_current(1:app(i, j).response_samples, idx), 1);
                    app(i, j).means.membrane_current = mean(app(i, j).membrane_current(1:app(i, j).response_samples, idx), 1);
                    app(i, j).means.rates = app(i, j).rates(1, idx);
                    app(i, j).means.npulses = app(i, j).npulses(1, idx);
                    app(i, j).means.spikes_per_stimulus = app(i, j).spikes_per_stimulus(idx, 1);
                end
            end
        end

        function mean_stats = compute_mean_stats(app)
            means = [app.means];
            pps = [means.rates];
            pulses = [means.npulses];
            sps = [means.spikes_per_stimulus];
            depol = [means.depolarization];
            hyperpol = [means.hyperpolarization];
            ge = [means.excitation];
            gi = [means.inhibition];
            net_ge = [means.resultant_excitation];
            net_gi = [means.resultant_inhibition];
            Ie = [means.excitatory_current];
            Ii = [means.inhibitory_current];
            Iact = [means.activation_current];
            Im = [means.membrane_current];
            Eth_cross = [app.threshold_crossing];
            mean_stats = [pps; pulses; sps; depol; hyperpol; ge; gi; net_ge; net_gi; Ie; Ii; Iact; Im; Eth_cross];
        end

        function meta_stats = compute_meta_stats(app, save)
            if nargin < 2
                save = 0;
            end
            paradigms = [app.paradigms];
            row_names = ["rate"; "pulses"; "sps"; "depol"; "hyperpol"; "ge"; "gi"; "net_ge"; "net_gi"; "Ie"; "Ii"; "Iact"; "Im"; "EthCross"];
            computed_stats = app.compute_mean_stats();
            meta_stats = array2table(computed_stats, "RowNames", row_names, "VariableNames",paradigms);
            if save == 1
                writetable(meta_stats, app(1).filename, "Sheet", "stats", "WriteRowNames", true);
            end
        end
    end
    
    %% Call method for estimator
    methods
        function app = call(app)
            app = app.filter_membrane_potential();
            app = app.adjust_membrane_potential_wrt_steady_state();
            app = app.compute_membrane_current();
            app = app.filter_membrane_current();
            app = app.compute_leakage_current();
            app = app.compute_activation_currents();
            app = app.filter_activation_currents();
            app = app.compute_passive_currents();
            app = app.compute_resultants();
            app = app.compute_polarizations();
%             app = app.reverse_estimate(1, 0);
        end

        function app = reverse_estimate(app, excitation_multiplier, inhibition_multiplier)
            if nargin < 2
                excitation_multiplier = 1;
                inhibition_multiplier = 0;
            elseif nargin < 3
                    inhibition_multiplier = 0;
            end
            [m, n] = size(app);
            for i = 1: m
                for j = 1: n
                    cm = app(i, j).membrane_capacitance(1, 1);
                    iinj = app(i, j).injected_current(1, 1);
                    ee = app(i, j).excitatory_reversal_potential(1, 1);
                    ei = app(i, j).inhibitory_reversal_potential(1, 1);
                    rin = app(i, j).input_resistance(1, 1);
                    er = app(i, j).steady_state_potential(1, 1);
                    a = app(i, j).alpha(1, 1);
                    et = app(i, j).threshold_potential(1, 1);
                    b = app(i, j).beta(1, 1);
                    vm = app(i, j).steady_state_potential(1, 1);
                    ge = excitation_multiplier.*(app(i, j).excitatory_conductance - app(i, j).excitatory_conductance(1, 1));
                    gi = inhibition_multiplier.*(app(i, j).inhibitory_conductance - app(i, j).inhibitory_conductance(1, 1));
                    for k = 2: 1: size(app(i, j).estimated_membrane_potential, 1)
                        Ie = ge(k).*(vm - ee).*(ge(k) > 0);
                        Ii = gi(k).*(vm - ei).*(gi(k) > 0);
                        Iactive = (a.*(vm - er).*(vm-et) + b.*(vm - er));
                        Ileak = (1./rin).*(vm - er);
                        vm = vm + (1./app(i, j).fs)*(1./cm)*(iinj - Ie - Ii - Ileak + Iactive);
                        if vm > et
                            vm = et;
                        end
                        app(i, j).estimated_membrane_potential(k, 1) = vm;
                        app(i, j).estimated_activation_current(k, 1) = Iactive;
                        app(i, j).estimated_excitatory_current(k, 1) = Ie;
                        app(i, j).estimated_inhibitory_current(k, 1) = Ii;
                        app(i, j).estimated_leakage_current(k, 1) = Ileak;
                    end
                    app(i, j).estimated_time = app(i, j).times;
                    ppp = find(app(i, j).estimated_membrane_potential >= app(i, j).threshold_potential(1, 1), 1);
                    if ~isempty(ppp)
                        app(i, j).threshold_crossing = ppp/app(i, j).fs;
                    end
                    app(i, j).estimated_membrane_current = app(i, j).fs.*cm.*[zeros(1, size(app(i, j).estimated_membrane_potential, 2)); diff(app(i, j).estimated_membrane_potential, 1, 1)];
                end
            end
        end
    end

    %% Plot methods
    methods
        function app = build_plots(app)
            app(1, 1).estimations_figure = figure('Name', strcat('Estimations for ', app(1, 1).filename), 'Position', [1200 100 900 900]);
            app(1, 1).stats_figure = figure('Name', strcat('Means for ', app(1, 1).filename), 'Position', [1200 200 900 900]);
        end

        function app = plot_without_representative_response(app)
            nplots = 8;
            figure_handle = figure(app(1, 1).estimations_figure);
            [m, n] = size(app);
            for i = 1: m
                tiledlayout(nplots, n);
                ax = cell(nplots, n);
                for k = 1: nplots
                    for j = 1: n
                        ax{j, k} = nexttile;
                        switch k
                            case 1
                                plot(app(i, j).times, app(i, j).membrane_potential, app(i, j).estimated_time, app(i, j).estimated_membrane_potential, 'k');
                                hold on;
                                plot(app(i, j).times, app(i, j).threshold_potential(:, 1), '--k');
                                hold off;
                                ax{j, k}.Title.String = app(i, j).paradigms;
                                ax{j, k}.Subtitle.String = [strcat("Eact=", num2str(app(i, j).activation_potential(1, :))), strcat('Diff: O=', num2str(size(app(1, 1).lowpass.Coefficients, 1)-1), '; Band=(', num2str(app(1, 1).lowpass.PassbandFrequency), ', ', num2str(app(1, 1).lowpass.StopbandFrequency), ')')];
                                if j == 1
                                    ax{j, k}.YLabel.String = 'Vm (V)';
                                end
                            case 2
                                plot(app(i, j).times, app(i, j).activation_current, app(i, j).estimated_time, app(i, j).estimated_activation_current, 'k');
                                ax{j, k}.Subtitle.String = strcat('xalpha=', num2str(app(i, j).alpha_multiplier(1, :)), '; xbeta=', num2str(app(i, j).beta_multiplier(1, :)));
                                if j == 1
                                    ax{j, k}.YLabel.String = 'Iact (A)';
                                end
                            case 3
                                plot(app(i, j).times, app(i, j).leakage_current, app(i, j).estimated_time, app(i, j).estimated_leakage_current(:, 1), 'k');
                                ax{j, k}.Subtitle.String = strcat('Rin=', num2str(app(i, j).input_resistance(1, 1)));
                                if j == 1
                                    ax{j, k}.YLabel.String = 'Ileak (A)';
                                end
                            case 4
                                plot(app(i, j).times, app(i, j).membrane_current, app(i, j).estimated_time, app(i, j).estimated_membrane_current(:, 1), 'k');
                                ax{j, k}.Subtitle.String = strcat('Cm=', num2str(app(i, j).membrane_capacitance(1, 1)));
                                if j == 1
                                    ax{j, k}.YLabel.String = 'Im (A)';
                                end
                            case 5
                                plot(app(i, j).times, app(i, j).excitatory_conductance, 'r', app(i, j).times, app(i, j).inhibitory_conductance, 'b');
                                ax{j, k}.Subtitle.String = strcat('Ee=', num2str(app(i , j).excitatory_reversal_potential(1, 1)), '; Ei=', num2str(app(i, j).inhibitory_reversal_potential(1, 1)));
                                if j == 1
                                    ax{j, k}.YLabel.String = 'G (S)';
                                end
                            case 6
                                area(app(i, j).times, app(i, j).resultant_excitaiton, 'FaceColor', 'r');
                                hold on;
                                area(app(i, j).times, app(i, j).resultant_inhibition, 'FaceColor', 'b');
                                hold off;
                                if j == 1
                                   ax{j, k}.YLabel.String = 'Resultant G(S)';
                                end
                            case 7
                                plot(app(i, j).times, app(i, j).excitatory_current(:, 1), 'r', app(i, j).times, app(i, j).inhibitory_current(:, 1), 'b');
                                if j == 1
                                    ax{j, k}.YLabel.String = 'Isyn (A)';
                                end
                            case 8
                                if ~(isempty(app(i, j).stimulus))
                                    plot(app(i, j).times, app(i, j).stimulus, 'k');
                                end
                                if j == 1
                                   ax{j, k}.YLabel.String = 'Stimulus (V)';
                                end
                        end
                        linkaxes([ax{j, :}], 'x');
                    end
                    linkaxes([ax{:, k}], 'y');
                end
                ylim('tight');
                xlim('tight');
            end
            saveas(figure_handle, strcat(app(1, 1).filename, '_development_traces.fig'));
        end

        function app = plot_with_representative_response(app)
            nplots = 5;
            figure_handle = figure(app(1, 1).estimations_figure);
            [m, n] = size(app);
            for i = 1: m
                tiledlayout(nplots, n);
                ax = cell(nplots, n);
                for k = 1: nplots
                    for j = 1: n
                        ax{j, k} = nexttile;
                        switch k
                            case 1
                                if ~(isempty(app(i, j).representative_response))
                                    plot(app(i, j).times, app(i, j).representative_response, 'k', 'Color', [0, 0, 0]);
                                end
                                if j == 1
                                   ax{j, k}.YLabel.String = 'Rep. (V)';
                                end
                            case 2
                                plot(app(i, j).times, app(i, j).membrane_potential, 'Color', [128 128 128]./255);
                                ax{j, k}.Title.String = app(i, j).paradigms;
                                ax{j, k}.Subtitle.String = [strcat("Eact=", num2str(app(i, j).activation_potential(1, :))), strcat('Diff: O=', num2str(size(app(1, 1).lowpass.Coefficients, 1)-1), '; Band=(', num2str(app(1, 1).lowpass.PassbandFrequency), ', ', num2str(app(1, 1).lowpass.StopbandFrequency), ')')];
                                if j == 1
                                    ax{j, k}.YLabel.String = 'Vm (V)';
                                end
                            case 4
                                plot(app(i, j).times, app(i, j).excitatory_conductance, 'Color', [255 0 0]./255);
                                hold on;
                                plot(app(i, j).times, app(i, j).inhibitory_conductance, 'Color', [0 190 216]./255);
                                hold off;
                                ax{j, k}.Subtitle.String = strcat('Ee=', num2str(app(i , j).excitatory_reversal_potential(1, 1)), '; Ei=', num2str(app(i, j).inhibitory_reversal_potential(1, 1)));
                                if j == 1
                                    ax{j, k}.YLabel.String = 'G (S)';
                                end
                            case 3
                                area(app(i, j).times, app(i, j).resultant_excitaiton, 'FaceColor', [255 0 0]./255);
                                hold on;
                                area(app(i, j).times, app(i, j).resultant_inhibition, 'FaceColor', [0 190 216]./255);
                                hold off;
                                if j == 1
                                   ax{j, k}.YLabel.String = 'Resultant G(S)';
                                end
                            case 5
                                if ~(isempty(app(i, j).stimulus))
                                    plot(app(i, j).times, app(i, j).stimulus, 'k');
                                end
                                if j == 1
                                   ax{j, k}.YLabel.String = 'Stimulus (V)';
                                end
                        end
%                         linkaxes([ax{j, :}], 'x');
                    end
%                     linkaxes([ax{:, k}], 'y');
                end
                linkaxes([ax{:, :}], 'x');
                linkaxes([ax{:, 1}, ax{:, 2}], 'y');
                linkaxes([ax{:, 3}, ax{:, 4}], 'y');
                linkaxes([ax{:, 5}], 'y');
                xlim('tight');
                ylim('tight');
            end
            saveas(figure_handle, strcat(app(1, 1).filename, '_production_traces.fig'));
        end

        function app = plot_estimations(app, production)
            if nargin < 2
                production = 0;
            end
            if production == 1
                app.plot_with_representative_response();
            else
                app.plot_without_representative_response();
            end
            
        end

        function app = plot_meta_stats_production(app, meta_stats)
            x = categorical(meta_stats.Properties.VariableNames);
            meta_stats = rows2vars(meta_stats);
            nplots = 1;
            figure_handle = figure(app(1, 1).stats_figure);
            tiledlayout(nplots, 1);
            ax = cell(nplots, 1);
            ax{1, 1} = nexttile;
            bar(meta_stats.ge, "BarWidth", 0.5, "FaceColor", [1, 0, 0]);
            hold on;
            bar(-1.*meta_stats.gi, "BarWidth", 0.5, "FaceColor", [0 190 216]./255);
            bar(meta_stats.net_ge, "BarWidth", 0.25, "FaceColor", [1, 0, 0]);
            bar(meta_stats.net_gi, "BarWidth", 0.25, "FaceColor", [0 190 216]./255);
            ylim([-1*max(max(meta_stats.ge), max(meta_stats.gi)),max(max(meta_stats.ge), max(meta_stats.gi))]);
            hold off;
            ylabel("syn conductances (S)");
            ax{1, 1}.YAxis.MinorTick = 'on';
            xticklabels(x);
            yyaxis right;
            plot(meta_stats.sps, '-ok', 'MarkerFaceColor',[0, 0, 0], 'MarkerSize', 6);
            ylim([-1*(max(meta_stats.sps)+0.2), (max(meta_stats.sps)+0.2)]);
            ylabel("sps");
            ax{1, 1}.YAxis(2).MinorTick = 'on';
            saveas(figure_handle, strcat(app(1, 1).filename, '_production_stats.fig'));
        end

        function app = plot_meta_stats_analysis(app, meta_stats)
            x = categorical(meta_stats.Properties.VariableNames);
            meta_stats = rows2vars(meta_stats);
            nplots = 5;
            figure_handle = figure(app(1, 1).stats_figure);
            tiledlayout(nplots, 1);
            ax = cell(nplots, 1);
            for k = 1: nplots
                ax{k, 1} = nexttile;
                switch k
                    case 1
                        bar(meta_stats.depol, "BarWidth", 0.5, "FaceColor", [1, 0, 0]);
                        hold on;
                        bar(-1.*meta_stats.hyperpol, "BarWidth", 0.5, "FaceColor", [0, 0, 1]);
                        ylim([-1*max(max(meta_stats.depol), max(meta_stats.hyperpol)),max(max(meta_stats.depol), max(meta_stats.hyperpol))]);
                        hold off;
                        ylabel("polarization (V)")
                    case 2
                        bar(meta_stats.ge, "BarWidth", 0.5, "FaceColor", [1, 0, 0]);
                        hold on;
                        bar(-1.*meta_stats.gi, "BarWidth", 0.5, "FaceColor", [0, 0, 1]);
                        bar(meta_stats.net_ge, "BarWidth", 0.25, "FaceColor", [1, 0.5, 0]);
                        bar(meta_stats.net_gi, "BarWidth", 0.25, "FaceColor", [0, 0.5, 1]);
                        ylim([-1*max(max(meta_stats.ge), max(meta_stats.gi)),max(max(meta_stats.ge), max(meta_stats.gi))]);
                        hold off;
                        ylabel("syn conductances (S)")
                    case 3
                        bar(meta_stats.Ie, "BarWidth", 0.5, "FaceColor", [1, 0, 0]);
                        hold on;
                        bar(-1.*meta_stats.Ii, "BarWidth", 0.5, "FaceColor", [0, 0, 1]);
                        ylim([-1*max(max(meta_stats.Ie), max(meta_stats.Ii)),max(max(meta_stats.Ie), max(meta_stats.Ii))]);
                        hold off;
                        ylabel("syn currents (A)")
                    case 4
                        bar(meta_stats.Iact, "BarWidth", 0.5, "FaceColor", [1, 0, 0]);
                        max_abs_Iact = abs(max(meta_stats.Iact));
                        ylim([-1*max_abs_Iact,max_abs_Iact]);
                        ylabel("act current (A)")
                    case 5
                        bar(meta_stats.Im, "BarWidth", 0.5, "FaceColor", [1, 0, 0]);  
                        max_abs_Im = abs(max(meta_stats.Im));
                        ylim([-1*max_abs_Im,max_abs_Im]);
                        ylabel("mem currents (A)")
                end
                ax{k, 1}.YAxis.MinorTick = 'on';
                xticklabels(x);
                yyaxis right;
                plot(meta_stats.sps, '-ok', 'MarkerFaceColor',[0, 0, 0], 'MarkerSize', 6);
                ylim([-1*(max(meta_stats.sps)+0.2), (max(meta_stats.sps)+0.2)]);
                ax{k, 1}.YAxis(2).MinorTick = 'on';
                ylabel("sps");
            end
            saveas(figure_handle, strcat(app(1, 1).filename, '_dev_stats.fig'));
        end

        function app = plot_meta_stats(app, meta_stats, production)
            if nargin < 3
                production = 0;
            end
            if production == 1
                app.plot_meta_stats_production(meta_stats);
            else
                app.plot_meta_stats_analysis(meta_stats);
            end
            
        end

    end
end