How to define a function with multiple handle in a loop?

I have a situation where I need to use a function that can take different number of input (depnding upon the number of treatment considered).

The function has a general expression (Attached as Pic for readability)

Currently I am defining the function (for different number of treatment) seperately (as bellow)

%% Parameters
par.a=80/365;% %; %bitting frequiency (per day)
par.g=1/10;      %mosquito death rate (per day)
par.n=12;        %duration of sporogony in mopsquito
par.b=.5;        %transmission probability: mosquito to human
par.c=.23;       %transmission probability: human to mosquito
par.r=1/60;        %rate of blood stage infection clearance
par.omega=1/425  ;   %hypnozoites death rate
par.nu=8.5;   %Number of hypnozoites per infection
par.alpha=1/332;%1/332;

p_blood=.9; %treatment effect
p_rad=.9;    %treatment effect

% tau is the current time
% s1,s2,s3 are treatment time.

pA=@(tau) par.alpha*(exp(-par.r*tau)-exp(-(par.alpha+par.omega)*tau))/(par.alpha+par.omega-par.r); %Base function without treatment


pA_rad=@(tau,s1) (1-p_blood)*exp(-par.r*(tau-s1)).*pA(s1)+(1-p_rad)*(pA(tau)-exp(-par.r*(tau-s1)).*pA(s1)); %with one treatment
pA_rad2=@(tau,s1,s2) (1-p_blood)*exp(-par.r*(tau-s2)).*pA_rad(s2,s1)+(1-p_rad)^2*(pA(tau)-exp(-par.r*(tau-s2)).*pA(s2)); %with two treatment
pA_rad3=@(tau,s1,s2,s3) (1-p_blood)*exp(-par.r*(tau-s3)).*pA_rad2(s3,s1,s2)+(1-p_rad)^3*(pA(tau)-exp(-par.r*(tau-s3)).*pA(s3)); %with three

as I couldn't figure out how to define these in a single function. The main issue is that as the number of treatment increases, the function handle increases as well and the function with larger handles utilise the function with lower handles for different inputs (as is attached pic). This is making me crazy. Is there any way that I can define these in a loop or nested function? Any help is much appreciated.

Answer 1

[Jan]

Is there a need to use an anonymous function? Why not writing a standard function?

function y = pA(tau, s1, s2, s3)

a=80/365;% %; %bitting frequiency (per day)
g=1/10;      %mosquito death rate (per day)
n=12;        %duration of sporogony in mopsquito
b=.5;        %transmission probability: mosquito to human
c=.23;       %transmission probability: human to mosquito
r=1/60;        %rate of blood stage infection clearance
omega=1/425  ;   %hypnozoites death rate
nu=8.5;   %Number of hypnozoites per infection
alpha=1/332;%1/332;

p_blood=.9; %treatment effect
p_rad=.9;    %treatment effect

% tau is the current time
% s1,s2,s3 are treatment time.

y0 = alpha*(exp(-r*tau)-exp(-(alpha+omega)*tau))/(alpha+omega-par.r);

switch nargin
    case 1
        y = y0; %Base function without treatment
    case 2
        y = (1-p_blood)*exp(-r*(tau-s1)).*pA(s1)+(1-p_rad)*(y0-exp(-r*(tau-s1)).*pA(s1)); %with one treatment
    case 3
        y = (1-p_blood)*exp(-r*(tau-s2)).*pA(s2,s1)+(1-p_rad)^2*(y0-exp(-r*(tau-s2)).*pA(s2)); %with two treatment
    case 4
        y = (1-p_blood)*exp(-r*(tau-s3)).*pA(s3,s1,s2)+(1-p_rad)^3*(y0-exp(-r*(tau-s3)).*pA(s3)); %with three
    otherwise
        error('Too many inputs');
end
end

[Md Nurul Anwar]

Hi Jan, thanks for the idea. The reason I am looking for a general function as I will need to do a sensitivity analysis, I need to do at least 15/20 treatment which means 15/20 functions. It would be great if I could do that without standart method

For your cases, can I pass the par structure as argument? I tried in this way but looks like it's not recognising 'par'

function y = pA(tau,s1,par)
y0 = par.alpha*(exp(-par.gamma*tau)-exp(-(par.alpha+par.mu)*tau))/(par.alpha+par.mu-par.gamma);
switch nargin
    case 1
        error('need parameters')
    case 2
        y = y0; %Base function without treatment
    case 3
        y = (1-par.p_blood)*exp(-par.gamma*(tau-s1)).*pA(s1,par)+(1-par.p_rad)*(y0-exp(-par.gamma*(tau-s1)).*pA(s1,par)); %with one treatment
    case 4
        y = (1-par.p_blood)*exp(-par.gamma*(tau-s2)).*pA_all(s2,s1,par)+(1-par.p_rad)^2*(y0-exp(-par.gamma*(tau-s2)).*pA(s2,par)); %with two treatment
    case 5
        y = (1-par.p_blood)*exp(-par.gamma*(tau-s3)).*pA_all(s3,s1,s2,par)+(1-par.p_rad)^3*(y0-exp(-par.gamma*(tau-s3)).*pA(s3,par)); %with three
    otherwise
        error('Too many inputs');
end
end

I am getting this error

>> pA(100,par)
Not enough input arguments.

Error in pA (line 2)
y0 =
par.alpha*(exp(-par.gamma*tau)-exp(-(par.alpha+par.mu)*tau))/(par.alpha+par.mu-par.gamma); 

[Jan]

For the call pA(100, par) the 2nd argument is copied to the variable s1 for this function header:

function y = pA(tau,s1,par)

and the 3rd input is undefined.

Better:

function y = pA(tau, varargin)
if nargin < 2
   error('need parameters')
end

par = varargin{end};
y0 = par.alpha*(exp(-par.gamma*tau)-exp(-(par.alpha+par.mu)*tau))/(par.alpha+par.mu-par.gamma);
switch nargin
    case 2
        y = y0;
    case 3
        s1 = varargin{1};
        y = (1-par.p_blood)*exp(-par.gamma*(tau-s1)).*pA(s1,par)+(1-par.p_rad)*(y0-exp(-par.gamma*(tau-s1)).*pA(s1,par)); %with one treatment
    case 4
        s1 = varargin{1};
        s2 = varargin{2};
        y = (1-par.p_blood)*exp(-par.gamma*(tau-s2)).*pA(s2,s1,par)+(1-par.p_rad)^2*(y0-exp(-par.gamma*(tau-s2)).*pA(s2,par)); %with two treatment
    case 5
        s1 = varargin{1};
        s2 = varargin{2};
        s3 = varargin{3};
        y = (1-par.p_blood)*exp(-par.gamma*(tau-s3)).*pA(s3,s1,s2,par)+(1-par.p_rad)^3*(y0-exp(-par.gamma*(tau-s3)).*pA(s3,par)); %with three
    otherwise
        error('Too many inputs');
end
end

Cleaner:

function y = pA(tau, par, s1, s2, s3)
if nargin < 2
   error('need parameters')
end

par = varargin{end};
y0 = par.alpha*(exp(-par.gamma*tau)-exp(-(par.alpha+par.mu)*tau))/(par.alpha+par.mu-par.gamma);
switch nargin
    case 2
        y = y0;
    case 3
        y = (1-par.p_blood)*exp(-par.gamma*(tau-s1)).*pA(s1,par)+(1-par.p_rad)*(y0-exp(-par.gamma*(tau-s1)).*pA(s1,par)); %with one treatment
    ... same as above without using varargin.    

It looks strange, that the order of arguments is:

pA(s1, par)
pA(s2, s1, par)
pA(s3, s1, s2, par)

If you get 20 terms, using a loop would be smarter. What is the aczual formula?

[Md Nurul Anwar]

First of all, thank for the idea. The formula is attached. This oder of input is being used within (as marked by red)

I tried the cleaner vertion as above.

function y = pA(tau, par, s1, s2, s3)
if nargin < 2
   error('need parameters')
end

par = varargin{end};
y0 = par.alpha*(exp(-par.gamma*tau)-exp(-(par.alpha+par.mu)*tau))/(par.alpha+par.mu-par.gamma);
switch nargin
    case 2
        y = y0;
    case 3
        y = (1-par.p_blood)*exp(-par.gamma*(tau-s1)).*pA(s1,par)+(1-par.p_rad)*(y0-exp(-par.gamma*(tau-s1)).*pA(s1,par)); %with one treatment
    case 4
        y = (1-par.p_blood)*exp(-par.gamma*(tau-s2)).*pA(s2,s1,par)+(1-par.p_rad)^2*(y0-exp(-par.gamma*(tau-s2)).*pA(s2,par)); %with two treatment
    case 5
        y = (1-par.p_blood)*exp(-par.gamma*(tau-s3)).*pA(s3,s1,s2,par)+(1-par.p_rad)^3*(y0-exp(-par.gamma*(tau-s3)).*pA(s3,par)); %with three
    otherwise
        error('Too many inputs');
end
end 

But getting an error

>> pA(100,par)
Error: File: pA.m Line: 6 Column: 7
Invalid syntax for calling function 'varargin' on the path. Use a valid syntax or
explicitly initialize 'varargin' to make it a variable.

But the initial one (with using varargin) is working

[Jan]

If you provide par as 2nd input, there is no need to call varargin anymore. Omit line 6:

par = varargin{end};  % Not needed here

But then the other calls do not match anymore: pA(s2,s1,par) and pA(s3,s1,s2,par).

[Jan]

This is a recursive impelemtation of your formula:

function y = pA(t, p, varargin)
N = nargin - 2;
if N < 0
    error('need parameters')
end

if N == 0  % No s as inputs:
    y = p.alpha / (p.alpha + p.mu - p.gamma) * ...
        (exp(-p.gamma * t) - exp(-(p.alpha + p.mu) * t)) ;
else  % s1..sN are provided as inputs - recursive call:
    sN = varargin{N};
    y  = (1 - p.p_blood) * exp(-p.gamma * (t - sN)) * ...
         pA(sN, p, varargin{1:N-1});
end
end

I've abbreviated "tau" to "t" and "par" to "p", because it improves the readability.

A recursive function can be converted to a loop in every case. But at least this implementation is much cleaner than a pile of nested anonymous functions.

During editing I struggled with the naming conventions: "par.alpha + par.omega - par.r". In the formula the "omega" is a mu and "r" is a gamma. The "t" is not a greek "tau".

I do no see, that your anonymous function for s1 matchs the formula:

pA_1 = @(tau,s1) (1 - p.p_blood) * exp(-p.gamma * (tau-s1)) * pA_0(s1) + ...
        (1-p.p_rad)*(pA_0(tau) - exp(-p.gamma*(tau-s1)).*pA_0(s1)); %with one treatment

While the first line is equivalent to

the 2nd line is not clear to me. p_rad does not appear in the formula.

[Md Nurul Anwar]

Appology for the confusion. I used pA_rad for the subscript 'r'. The complete formula is below ( as in code).

The 2nd line as you mentioned is the y0 (here pA :no treatment)

Thanks a lot for these idea so far.

[Jan]

@Md Nurul Anwar: Please care for posting the correct formula directly in the question the next time.

function y = pA(t, p, varargin)
N = nargin - 2;
if N == 0  % No s as inputs:
    y = p.alpha / (p.alpha + p.mu - p.gamma) * ...
        (exp(-p.gamma * t) - exp(-(p.alpha + p.mu) * t)) ;
elseif N > 0  % s1..sN are provided as inputs - recursive call:
    sN = varargin{N};
    y  = (1 - p.p_blood) * exp(-p.gamma * (t - sN)) * ...
         pA(sN, p, varargin{1:N-1}) + ...
         (1 - p.p_rad)^N * (pA(t, p) - exp(-p.gamma * (t - sN)) * pA(sN, p));
else
    error('need parameters')
end
end

Please test this, if it produces the correct answers. It works with an arbitrary number of inputs s1, ...sN limited by the recursion depth of 500. If this is a problem or the function is time-critical, convert it to a loop.

An option would be to provide s as vector instead of a list of scalars.

[Md Nurul Anwar]

It works! Thanks a lot Jan.