% there are several ways for creating these arrays
% 1 - squared brackets
A = [1 3 5; 3 6 3; 9 8 5]
% 2 - functions for array creation
A = rand(3)
A = magic(3) 

b = [3; 4: 8] % 3 x 1 - column vector!

A * b(1)
A * b % matrix multiplication
A.*b % element-wise multiplication

B = [b b b] % horizontal concatenation

A = A.*B

A(2,1)
A(3,:)
A(2:3,2:3)

t = linspace(0,2*pi,100); % use linspace to get n equidistant numbers
tnew = sin(t);
tc = asin(tnew);
t == tc % you could use `all(t==tc)` to check whether all are equal
notequalpos = find(t~=tc);

% use the documentation to find the functions needed
max(tnew)
min(tnew)
abs(tnew)

sqrt(t)
log10(t)

str = 'Sven'

upper(str)
lower(str)

cstr = {str, 'EMBL'}

%% import the data using autogenerated code
% assign the columns to data vectors
[time, R, L, RL] = importdata('RLBinding.csv',1, 17)

%% plot using the autogenerated code from the plot tools
createfigure(time, R, L, RL) % the exact function signature 
% depends on the generated code and your sequence of steps

%% get a list of csv files in the current folder
files = dir('*.csv');
% calculate the best number of rows and columns for 
% the subplot
nfiles = length(files);
nrows = round(sqrt(nfiles));
ncols = ceil(sqrt(nfiles));

%% loop through all files

for ii = 1:nfiles
    %% import the data using autogenerated code
    % assign the columns to data vectors
    [time, R, L, RL] = importdata(files(ii).name,1, 17)

    %% plot using the autogenerated code from the plot tools
    subplot(nrows,ncols,ii), createfigure(time, R, L, RL) 
    % the exact function signature 
    % depends on the generated code and your sequence of steps

end

% the ODE model
drdt = - kf*r*l;
dldt = - kf*r*l;

drldt = + kf*r*l;

function rhs = oderhs(t,x)
% wrap a function header around the previous code block and 
% save as oderhs.m

kf = 1;
% extract the input vector into the state variables
r = x(1);
l = x(2);
rl = x(3);

% ODE system
drdt = - kf*r*l;
dldt = - kf*r*l;
drldt = + kf*r*l;

% any MATLAB ODE solver expects a column vector as output argument
rhs = [drdt; dldt; drldt];

function [t,x] = odeeuler(t0, tf, x0, h)
% solves an ODE system using the Euler method
% if you want it to be a general function that could be 
% used for any ODE system, specify the ODE rhs function
% as an input argument:
% function [t,x] = odeeuler(@oderhs, t0, tf, x0, h)

% create the vector of time points
% at which to solve
t = t0:h:tf;
nsteps = length(t);
% preallocate the solution array
x = zeros(nsteps,3);
x(1,:) = x0;

% loop from 2 to the number of wanted solution points
for ii = 2:nsteps
    x(ii,:) = x(ii-1,:) + h*oderhs(t(ii),x(ii-1,:));
end

x0 = [10 8 0];
[t,x] = odeeuler(0, 20, x0, 0.1);
% extract into column vectors
Rs = x(:,1);
Ls = x(:,2);
RLs = x(:,3);

% we could be using the createfigure function we created earlier,
% or just call plot()
plot(t,Rs,'r-d', t,Ls,'g-o', t,RLs,'b-x')

function obj = objfunction(parameter)

[time, R, L, RL] = importdata(files(ii).name,1, 17);
% put the dataset into a single matrix so that we simply
% do a matrix subtraction
data = [R L RL];
ndata = length(time);

% simulate the model using the parameter 
% instead of our own ODE solver, we'll use one
% of the in-built solvers and define the model rhs as 
% a nested function to be able to pass in the parameter
% (matlab ODE solvers can only use functions that accept
% exactly two arguments)

[t, x] = ode15s(@odefunc,time, data(1,:));
res = data-x;
% calculate the least-squares term
% when using lsqnonlin in the calling function,
% simply return the column-vectorized residual 
% (here: res)
obj = 1/(ndata-1)*sum(sqrt((res(:)).^2));

% nested ode function
function rhs = odefunc(t,x)
    % the 'parameter' variable in the parent function
    % is visible in this nested function (see the 
    % section on scope)
    kf = parameter;
    % extract the input vector into the state variables
    r = x(1);
    l = x(2);
    rl = x(3);

    % ODE system
    drdt = - kf*r*l;
    dldt = - kf*r*l;
    drldt = kf*r*l;
end

% when using nested functions the main function needs
% to be closed with 'end'
end

function bestparameter = fitparameter()

[time, R, L, RL] = importdata(files(ii).name,1, 17);
data = [R L RL]; % this dataset (and time) is/are
% visible in objfunction

% when we use lsqnonlin (available in the optimization
% toolbox, we need to return only a vector of the 
% residuals) 
% the third and fourth argument are the lower and upper
% bounds on the argument that should be minimized
[bestparameter] = lsqnonlin(@objfunction, p0, 0, Inf);

function res = objfunction(parameter)

    ndata = length(time);
    [t, x] = ode15s(@odefunc,time, data(1,:));
    res = data-x;
    % convert the ndata x 3 matrix to a column vector
    % using the : operator
    res = res(:);

    function rhs = odefunc(t,x)
        kf = parameter;
        % extract the input vector into the state variables
        r = x(1);
        l = x(2);
        rl = x(3);

        % ODE system
        drdt = - kf*r*l;
        dldt = - kf*r*l;
        drldt = kf*r*l;
    end
end
end

files = dir('*.tif');

for ii = 1:length(files)
    iminfo = imfinfo(files(ii).name);
    % allocate the array imstack
    imstack = zeros(iminfo(1).Height, iminfo(1).Width,length(iminfo);

    % read in the tif into imstack using imread
    for jj = 1:length(iminfo)
        imstack(:,:,jj) = imread(files(ii).name,'Index',jj);
    end

    % extract the fluorescent images into another image stack
    imstackf = imstack(:,:,1:end/2);

    % maximum intensity projection
    imm = max(imstackf,[],3);

    % convert it to (0,1) intensity range
    img = mat2gray(imm);

    % median filtering to remove noise
    img = medfilt2(img);

    % sharpen the image
    img = imsharpen(img);

    % calculate the threshold
    threshold = graythresh(img);

    % segment the image using the threshold
    imbw = im2bw(img, threshold);

    % identify the objects
    imobjs = bwlabel(imbw);

    % extracting properties of the objects
    improps = regionprops(imobjs,'all');
    ncells = length(improps);

    % the number of cells is the number of identified objects
    sprintf('This image contains %d cells', ncells)
end