Error: unrecognized function or variable "bi2de"

50 views (last 30 days)
Abdullahi
Abdullahi on 27 Aug 2024
Commented: Abdullahi on 27 Aug 2024
Hi there.. have been trying to run this code to fit a model but it keeps producing an error.. Please someone should help
% Data fitting for the covid-19 outbreak in Nigeria, with reference to
% Lagos State Using the Genetic Algorithm (GA).
function covid19_GA
clear all
clc
global a x1 beta_c ep_g c_g sigma nu theta psi d_o d_d gamma_i alpha ...
delta ep_s ep_i ep_a best E0 A0 I0
ep_g = 0; c_g = 0; delta = 0; ep = 0; ep_s = ep; ep_a=ep; ep_i=ep;
sigma = 1/5.2; nu = 0.5;
gamma_i = 1/15;
gamma_a = 0.13978; gamma_o = 0.13978;
d_o = 0.015; d_d = 0.015;
alpha = 0.5;
%==========================================================================
% Lagos data, starting from March 16, 2020 till May 3, 2020.
day = [1 3 4 6 7 8 9 10 11 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 ...
43 44 45 46 47 48];
Totalcases = [1 5 9 19 22 25 29 32 44 59 68 81 82 91 98 109 109 120 120 130 145 158 166 174 176 189 214 232 251 283 ...
306 376 376 430 504 582 657 689 731 764 844 931 976 1006 1068];
Activecases = [1 5 8 19 22 25 29 32 43 58 67 75 76 85 81 86 86 91 87 96 111 117 118 119 116 123 139 140 154 182 200 ...
266 266 309 383 460 525 548 584 602 682 718 756 753 791];
%==========================================================================
pop_size = 100; chrom_len = 100; pm = 0.01; pc = 0.5; max_gen = 50;
%suggested run: [P,best] = ga(100,100,0.01,0.5,200)
[P,best] = gaCovid(pop_size,chrom_len,pm,pc,max_gen);
[beta_c,theta,psi, E0, A0, I0] = gene_to_values(best);
%Control Reproduction Number
R_c = beta_c*(1-ep_g*c_g)*(1-delta)*(1-ep_s)*((alpha*nu*(1-ep_a)/...
(theta+gamma_a)) + ((1-nu)*(1-ep_i)/(psi+d_o+gamma_o)))
[t,y] = ode15s(@covidmodel1, [0,52],[14368332,E0,A0,I0,1,0,0]);
plot(t,y(:,7),'-b',day,Totalcases,'rs')
xlabel('Time (days)'),ylabel('Cumulative number of reported cases')
legend('Model','Lagos COVID-19 Data')
% plot(t,y(:,5),'-b',day,Activecases,'rs')
% subplot(211), plot(t,y(:,7),'-b',day,Totalcases,'rs')
% subplot(212), plot(t,y(:,5),'-b',day,Activecases,'rs')
function dx = covidmodel1(t,x)
dx = [0;0;0;0;0;0;0];
N = x(1)+x(2)+x(3)+x(4)+x(5)+x(6);
beta = beta_c*(1-ep_g*c_g)*(1-delta)*(1-ep_s);
dx(1) = -beta*x(1)*(alpha*(1-ep_a)*x(3)+(1-ep_i)*x(4))/(N-x(5));
dx(2) = beta*x(1)*(alpha*(1-ep_a)*x(3)+(1-ep_i)*x(4))/(N-x(5)) - sigma*x(2);
dx(3) = nu*sigma*x(2) - theta*x(3) - gamma_a*x(3);
dx(4) = (1-nu)*sigma*x(2) - psi*x(4) - d_o*x(4) - gamma_o*x(4);
dx(5) = theta*x(3) + psi*x(4) - gamma_i*x(5) -d_d*x(5); % Active cases, detected (asymptomatic and symptomatic)
dx(6) = gamma_i*x(5) + gamma_a*x(3) + gamma_o*x(4);
dx(7) = theta*x(3) + psi*x(4); % Cumulative number of new cases reported (asymptomatic and symptomatic)
end
function [P,best] = gaCovid(pop_size,chrom_len,pm,pc,max_gen)
P = initialize(pop_size,chrom_len);
fit = maxones_fitness(P);
gen = 1;
while gen<=max_gen && max(fit) < chrom_len
P = tournament_selection(P,fit,2);
P = two_point_crossover(P,pc);
P = point_mutation(P,pm);
fit = fitness(P);
max_fit(gen) = max(fit);
mean_fit(gen) = mean(fit);
disp_info(gen,max_fit(gen));
gen = gen+1;
end
disp(sprintf('Generation: %d',gen));
disp(sprintf('Best fitness: %d\n',max_fit(end)));
[m,ind] = max(fit);
best = P(ind,:);
end
function [P] = initialize(pop_size,chrom_length)
P = round(rand(pop_size,chrom_length));
end
function [fit] = maxones_fitness(P)
for i=1:size(P,1)
fit(i) = length(find(P(i,:)));
end
end
function [P_new] = tournament_selection(P,fit,tourn_size)
for i=1:size(P,1)
t = ceil(rand(1,tourn_size)*size(P,1));
[max_fit,winner] = max(fit(t));
P_new(i,:) = P(t(winner),:);
end
end
function [P_new] = two_point_crossover(P,pc)
mating_list = randperm(size(P,1));
P_new = [];
while ~isempty(mating_list)
pair = mating_list(1:2);
mating_list(1:2) = [];
if rand<pc
crossover_points = ceil(rand(1,2)*(size(P,2)));
point1 = min(crossover_points);
point2 = max(crossover_points);
individual1 = P(pair(1),:);
individual2 = P(pair(2),:);
individual1(point1:point2) = P(pair(2),point1:point2);
individual2(point1:point2) = P(pair(1),point1:point2);
P_new = [P_new;individual1;individual2];
else
P_new = [P_new;P(pair,:)];
end
end
end
function [P_new] = point_mutation(P,pm)
r = rand(size(P));
mutation_list = find(r<pm);
P_new = P;
P_new(mutation_list(find(P(mutation_list)==1))) = 0;
P_new(mutation_list(find(P(mutation_list)==0))) = 1;
end
function [] = disp_info(gen,fit)
if mod(gen,10)==0
disp(sprintf('Generation: %d',gen));
disp(sprintf('Best fitness: %d\n',fit));
end
end
function [P_new] = roulette_selection(P,fit)
fit = (fit - min(fit)).^2;
fit = cumsum(fit);
fit = fit/max(fit);
P_new = [];
for i=1:size(P,1)
f = find(fit>rand);
P_new = [P_ne';P(f(1),:)];
end
end
function fit = fitness(P)
fit = zeros(size(P,1),1);
for i=1:size(P,1)
[beta_c, theta, psi, E0, A0, I0] = gene_to_values(P(i,:));
[t,y] = ode15s(@covidmodel1, [1,day(end)],[14368332,E0,A0,I0,1,0,0]);
error = squared_error(day,Totalcases,t,y(:,7));
%error = squared_error(day,Activecases,t,y(:,5));
% error = squared_error(day,Totalcases,Activecases,t,y(:,7),y(:,5)); % 2 dataset
fit(i) = error^-1;
end
end
function [beta_c, theta, psi, E0, A0, I0] = gene_to_values(gene)
beta_c_min_exponent = -7;%-7;
beta_c_max_exponent = 1;
theta_min_exponent = -3;
theta_max_exponent = 1;
psi_min_exponent = -3;%-7;
psi_max_exponent = 1;
E0_min_exponent = 0;
E0_max_exponent = 2;
A0_min_exponent = 0;
A0_max_exponent = 2;
I0_min_exponent = 0;
I0_max_exponent = 2;
gene_len = size(gene, 2)/2;
beta_c_decimal = bi2de(gene(1:gene_len));
theta_decimal = bi2de(gene((gene_len+1):end));
psi_decimal = bi2de(gene((gene_len+1):end));
E0_decimal = bi2de(gene((gene_len+1):end));
A0_decimal = bi2de(gene((gene_len+1):end));
I0_decimal = bi2de(gene((gene_len+1):end));
beta_c_exponent = beta_c_min_exponent +...
((beta_c_max_exponent-beta_c_min_exponent)*...
beta_c_decimal/2^gene_len);
theta_exponent = theta_min_exponent +...
((theta_max_exponent-theta_min_exponent)*...
theta_decimal/2^gene_len);
psi_exponent = psi_min_exponent +...
((psi_max_exponent-psi_min_exponent)*...
psi_decimal/2^gene_len);
E0_exponent = E0_min_exponent +...
((E0_max_exponent-E0_min_exponent)*...
E0_decimal/2^gene_len);
A0_exponent = A0_min_exponent +...
((A0_max_exponent-A0_min_exponent)*...
A0_decimal/2^gene_len);
I0_exponent = I0_min_exponent +...
((I0_max_exponent-I0_min_exponent)*...
I0_decimal/2^gene_len);
beta_c = 10^beta_c_exponent;
theta = 10^theta_exponent;
psi = 10^psi_exponent;
E0 = 10^E0_exponent;
A0 = 10^A0_exponent;
I0 = 10^I0_exponent;
end
%=========================================================================
% Fitting 1 dataset
function error = squared_error(data_time,data_points,fn_time,fn_points)
fn_values = interp1(fn_time,fn_points,data_time);
error = sum((fn_values - data_points).^2);
end
%==========================================================================
% % Fitting 2 datasets
% function error = squared_error(data_time,data_points1,data_points2,fn_time,fn_points1,fn_points2)
%
% fn_values1 = interp1(fn_time,fn_points1,data_time); fn_values2 = interp1(fn_time,fn_points2,data_time);
%
% error = sum((fn_values1 - data_points1).^2) + sum((fn_values2 - data_points2).^2);
% end
%=========================================================================
end
  3 Comments
Show 1 older comment
Torsten
Torsten on 27 Aug 2024
Edited: Torsten on 27 Aug 2024
Works in R2024a (see above).
But this looks strange:
P_new = [];
for i=1:size(P,1)
f = find(fit>rand);
P_new = [P_ne';P(f(1),:)];
end
Shouldn't it be
P_new = [P_new;P(f(1),:)];
?
Mathieu NOE
Mathieu NOE on 27 Aug 2024
bi2de : looks like a binary to decimal conversion function
have you tried to use the bin2dec function instead ? : Convert text representation of binary integer to double value - MATLAB bin2dec - MathWorks France

Sign in to comment.

Sign in to answer this question.

Accepted Answer

Voss
Voss on 27 Aug 2024
Ran in:
bi2de is a function in the Communications Toolbox; based on the error message you received, I guess that toolbox is not installed on your system or is not licensed. You can use bin2dec instead, which is a function in base MATLAB.
However, bi2de takes numeric vector input, whereas bin2dec expects character vector input, so you have to convert the numeric vector of 0s and 1s to a character vector of '0's and '1's. You can make your own function called bi2de (so you don't have to change your existing code) which does the numeric to char conversion and then calls bin2dec.
That function might look like this:
function out = bi2de(in)
out = bin2dec(char(in+'0'));
end
Including that function in a function file on your path or inside your covid19_GA function (as I've done here) allows covid19_GA to run:
covid19_GA
Generation: 10 Best fitness: 2.739967e-05 Generation: 20 Best fitness: 2.788751e-05 Generation: 30 Best fitness: 2.744057e-05 Generation: 40 Best fitness: 2.801428e-05 Generation: 50 Best fitness: 2.744704e-05 Generation: 51 Best fitness: 2.744704e-05
R_c = 1.8435
% Data fitting for the covid-19 outbreak in Nigeria, with reference to
% Lagos State Using the Genetic Algorithm (GA).
function covid19_GA
% clear all
clc
global a x1 beta_c ep_g c_g sigma nu theta psi d_o d_d gamma_i alpha ...
delta ep_s ep_i ep_a best E0 A0 I0
ep_g = 0; c_g = 0; delta = 0; ep = 0; ep_s = ep; ep_a=ep; ep_i=ep;
sigma = 1/5.2; nu = 0.5;
gamma_i = 1/15;
gamma_a = 0.13978; gamma_o = 0.13978;
d_o = 0.015; d_d = 0.015;
alpha = 0.5;
%==========================================================================
% Lagos data, starting from March 16, 2020 till May 3, 2020.
day = [1 3 4 6 7 8 9 10 11 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 ...
43 44 45 46 47 48];
Totalcases = [1 5 9 19 22 25 29 32 44 59 68 81 82 91 98 109 109 120 120 130 145 158 166 174 176 189 214 232 251 283 ...
306 376 376 430 504 582 657 689 731 764 844 931 976 1006 1068];
Activecases = [1 5 8 19 22 25 29 32 43 58 67 75 76 85 81 86 86 91 87 96 111 117 118 119 116 123 139 140 154 182 200 ...
266 266 309 383 460 525 548 584 602 682 718 756 753 791];
%==========================================================================
pop_size = 100; chrom_len = 100; pm = 0.01; pc = 0.5; max_gen = 50;
%suggested run: [P,best] = ga(100,100,0.01,0.5,200)
[P,best] = gaCovid(pop_size,chrom_len,pm,pc,max_gen);
[beta_c,theta,psi, E0, A0, I0] = gene_to_values(best);
%Control Reproduction Number
R_c = beta_c*(1-ep_g*c_g)*(1-delta)*(1-ep_s)*((alpha*nu*(1-ep_a)/...
(theta+gamma_a)) + ((1-nu)*(1-ep_i)/(psi+d_o+gamma_o)))
[t,y] = ode15s(@covidmodel1, [0,52],[14368332,E0,A0,I0,1,0,0]);
plot(t,y(:,7),'-b',day,Totalcases,'rs')
xlabel('Time (days)'),ylabel('Cumulative number of reported cases')
legend('Model','Lagos COVID-19 Data')
% plot(t,y(:,5),'-b',day,Activecases,'rs')
% subplot(211), plot(t,y(:,7),'-b',day,Totalcases,'rs')
% subplot(212), plot(t,y(:,5),'-b',day,Activecases,'rs')
function out = bi2de(in)
out = bin2dec(char(in+'0'));
end
function dx = covidmodel1(t,x)
dx = [0;0;0;0;0;0;0];
N = x(1)+x(2)+x(3)+x(4)+x(5)+x(6);
beta = beta_c*(1-ep_g*c_g)*(1-delta)*(1-ep_s);
dx(1) = -beta*x(1)*(alpha*(1-ep_a)*x(3)+(1-ep_i)*x(4))/(N-x(5));
dx(2) = beta*x(1)*(alpha*(1-ep_a)*x(3)+(1-ep_i)*x(4))/(N-x(5)) - sigma*x(2);
dx(3) = nu*sigma*x(2) - theta*x(3) - gamma_a*x(3);
dx(4) = (1-nu)*sigma*x(2) - psi*x(4) - d_o*x(4) - gamma_o*x(4);
dx(5) = theta*x(3) + psi*x(4) - gamma_i*x(5) -d_d*x(5); % Active cases, detected (asymptomatic and symptomatic)
dx(6) = gamma_i*x(5) + gamma_a*x(3) + gamma_o*x(4);
dx(7) = theta*x(3) + psi*x(4); % Cumulative number of new cases reported (asymptomatic and symptomatic)
end
function [P,best] = gaCovid(pop_size,chrom_len,pm,pc,max_gen)
P = initialize(pop_size,chrom_len);
fit = maxones_fitness(P);
gen = 1;
while gen<=max_gen && max(fit) < chrom_len
P = tournament_selection(P,fit,2);
P = two_point_crossover(P,pc);
P = point_mutation(P,pm);
fit = fitness(P);
max_fit(gen) = max(fit);
mean_fit(gen) = mean(fit);
disp_info(gen,max_fit(gen));
gen = gen+1;
end
disp(sprintf('Generation: %d',gen));
disp(sprintf('Best fitness: %d\n',max_fit(end)));
[m,ind] = max(fit);
best = P(ind,:);
end
function [P] = initialize(pop_size,chrom_length)
P = round(rand(pop_size,chrom_length));
end
function [fit] = maxones_fitness(P)
for i=1:size(P,1)
fit(i) = length(find(P(i,:)));
end
end
function [P_new] = tournament_selection(P,fit,tourn_size)
for i=1:size(P,1)
t = ceil(rand(1,tourn_size)*size(P,1));
[max_fit,winner] = max(fit(t));
P_new(i,:) = P(t(winner),:);
end
end
function [P_new] = two_point_crossover(P,pc)
mating_list = randperm(size(P,1));
P_new = [];
while ~isempty(mating_list)
pair = mating_list(1:2);
mating_list(1:2) = [];
if rand<pc
crossover_points = ceil(rand(1,2)*(size(P,2)));
point1 = min(crossover_points);
point2 = max(crossover_points);
individual1 = P(pair(1),:);
individual2 = P(pair(2),:);
individual1(point1:point2) = P(pair(2),point1:point2);
individual2(point1:point2) = P(pair(1),point1:point2);
P_new = [P_new;individual1;individual2];
else
P_new = [P_new;P(pair,:)];
end
end
end
function [P_new] = point_mutation(P,pm)
r = rand(size(P));
mutation_list = find(r<pm);
P_new = P;
P_new(mutation_list(find(P(mutation_list)==1))) = 0;
P_new(mutation_list(find(P(mutation_list)==0))) = 1;
end
function [] = disp_info(gen,fit)
if mod(gen,10)==0
disp(sprintf('Generation: %d',gen));
disp(sprintf('Best fitness: %d\n',fit));
end
end
function [P_new] = roulette_selection(P,fit)
fit = (fit - min(fit)).^2;
fit = cumsum(fit);
fit = fit/max(fit);
P_new = [];
for i=1:size(P,1)
f = find(fit>rand);
P_new = [P_ne';P(f(1),:)];
end
end
function fit = fitness(P)
fit = zeros(size(P,1),1);
for i=1:size(P,1)
[beta_c, theta, psi, E0, A0, I0] = gene_to_values(P(i,:));
[t,y] = ode15s(@covidmodel1, [1,day(end)],[14368332,E0,A0,I0,1,0,0]);
error = squared_error(day,Totalcases,t,y(:,7));
%error = squared_error(day,Activecases,t,y(:,5));
% error = squared_error(day,Totalcases,Activecases,t,y(:,7),y(:,5)); % 2 dataset
fit(i) = error^-1;
end
end
function [beta_c, theta, psi, E0, A0, I0] = gene_to_values(gene)
beta_c_min_exponent = -7;%-7;
beta_c_max_exponent = 1;
theta_min_exponent = -3;
theta_max_exponent = 1;
psi_min_exponent = -3;%-7;
psi_max_exponent = 1;
E0_min_exponent = 0;
E0_max_exponent = 2;
A0_min_exponent = 0;
A0_max_exponent = 2;
I0_min_exponent = 0;
I0_max_exponent = 2;
gene_len = size(gene, 2)/2;
beta_c_decimal = bi2de(gene(1:gene_len));
theta_decimal = bi2de(gene((gene_len+1):end));
psi_decimal = bi2de(gene((gene_len+1):end));
E0_decimal = bi2de(gene((gene_len+1):end));
A0_decimal = bi2de(gene((gene_len+1):end));
I0_decimal = bi2de(gene((gene_len+1):end));
beta_c_exponent = beta_c_min_exponent +...
((beta_c_max_exponent-beta_c_min_exponent)*...
beta_c_decimal/2^gene_len);
theta_exponent = theta_min_exponent +...
((theta_max_exponent-theta_min_exponent)*...
theta_decimal/2^gene_len);
psi_exponent = psi_min_exponent +...
((psi_max_exponent-psi_min_exponent)*...
psi_decimal/2^gene_len);
E0_exponent = E0_min_exponent +...
((E0_max_exponent-E0_min_exponent)*...
E0_decimal/2^gene_len);
A0_exponent = A0_min_exponent +...
((A0_max_exponent-A0_min_exponent)*...
A0_decimal/2^gene_len);
I0_exponent = I0_min_exponent +...
((I0_max_exponent-I0_min_exponent)*...
I0_decimal/2^gene_len);
beta_c = 10^beta_c_exponent;
theta = 10^theta_exponent;
psi = 10^psi_exponent;
E0 = 10^E0_exponent;
A0 = 10^A0_exponent;
I0 = 10^I0_exponent;
end
%=========================================================================
% Fitting 1 dataset
function error = squared_error(data_time,data_points,fn_time,fn_points)
fn_values = interp1(fn_time,fn_points,data_time);
error = sum((fn_values - data_points).^2);
end
%==========================================================================
% % Fitting 2 datasets
% function error = squared_error(data_time,data_points1,data_points2,fn_time,fn_points1,fn_points2)
%
% fn_values1 = interp1(fn_time,fn_points1,data_time); fn_values2 = interp1(fn_time,fn_points2,data_time);
%
% error = sum((fn_values1 - data_points1).^2) + sum((fn_values2 - data_points2).^2);
% end
%=========================================================================
end

More Answers (1)

Mathieu NOE
Mathieu NOE on 27 Aug 2024
Thanks to this post it was easy to (re)create the missing / unknown bi2de function (see end of the code) : Most efficient way to convert from binary to dec without bin2dec or bi2de - MATLAB Answers - MATLAB Central (mathworks.com)
Credits to Walter Robertson !!
this is a good and fast alternative to the native bin2dec matlab function I mentionned above (but you could also use it)
updated code
% Data fitting for the covid-19 outbreak in Nigeria, with reference to
% Lagos State Using the Genetic Algorithm (GA).
function covid19_GA
clear all
clc
global a x1 beta_c ep_g c_g sigma nu theta psi d_o d_d gamma_i alpha ...
delta ep_s ep_i ep_a best E0 A0 I0
ep_g = 0; c_g = 0; delta = 0; ep = 0; ep_s = ep; ep_a=ep; ep_i=ep;
sigma = 1/5.2; nu = 0.5;
gamma_i = 1/15;
gamma_a = 0.13978; gamma_o = 0.13978;
d_o = 0.015; d_d = 0.015;
alpha = 0.5;
%==========================================================================
% Lagos data, starting from March 16, 2020 till May 3, 2020.
day = [1 3 4 6 7 8 9 10 11 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 ...
43 44 45 46 47 48];
Totalcases = [1 5 9 19 22 25 29 32 44 59 68 81 82 91 98 109 109 120 120 130 145 158 166 174 176 189 214 232 251 283 ...
306 376 376 430 504 582 657 689 731 764 844 931 976 1006 1068];
Activecases = [1 5 8 19 22 25 29 32 43 58 67 75 76 85 81 86 86 91 87 96 111 117 118 119 116 123 139 140 154 182 200 ...
266 266 309 383 460 525 548 584 602 682 718 756 753 791];
%==========================================================================
pop_size = 100; chrom_len = 100; pm = 0.01; pc = 0.5; max_gen = 50;
%suggested run: [P,best] = ga(100,100,0.01,0.5,200)
[P,best] = gaCovid(pop_size,chrom_len,pm,pc,max_gen);
[beta_c,theta,psi, E0, A0, I0] = gene_to_values(best);
%Control Reproduction Number
R_c = beta_c*(1-ep_g*c_g)*(1-delta)*(1-ep_s)*((alpha*nu*(1-ep_a)/...
(theta+gamma_a)) + ((1-nu)*(1-ep_i)/(psi+d_o+gamma_o)))
[t,y] = ode15s(@covidmodel1, [0,52],[14368332,E0,A0,I0,1,0,0]);
plot(t,y(:,7),'-b',day,Totalcases,'rs')
xlabel('Time (days)'),ylabel('Cumulative number of reported cases')
legend('Model','Lagos COVID-19 Data')
% plot(t,y(:,5),'-b',day,Activecases,'rs')
% subplot(211), plot(t,y(:,7),'-b',day,Totalcases,'rs')
% subplot(212), plot(t,y(:,5),'-b',day,Activecases,'rs')
function dx = covidmodel1(t,x)
dx = [0;0;0;0;0;0;0];
N = x(1)+x(2)+x(3)+x(4)+x(5)+x(6);
beta = beta_c*(1-ep_g*c_g)*(1-delta)*(1-ep_s);
dx(1) = -beta*x(1)*(alpha*(1-ep_a)*x(3)+(1-ep_i)*x(4))/(N-x(5));
dx(2) = beta*x(1)*(alpha*(1-ep_a)*x(3)+(1-ep_i)*x(4))/(N-x(5)) - sigma*x(2);
dx(3) = nu*sigma*x(2) - theta*x(3) - gamma_a*x(3);
dx(4) = (1-nu)*sigma*x(2) - psi*x(4) - d_o*x(4) - gamma_o*x(4);
dx(5) = theta*x(3) + psi*x(4) - gamma_i*x(5) -d_d*x(5); % Active cases, detected (asymptomatic and symptomatic)
dx(6) = gamma_i*x(5) + gamma_a*x(3) + gamma_o*x(4);
dx(7) = theta*x(3) + psi*x(4); % Cumulative number of new cases reported (asymptomatic and symptomatic)
end
function [P,best] = gaCovid(pop_size,chrom_len,pm,pc,max_gen)
P = initialize(pop_size,chrom_len);
fit = maxones_fitness(P);
gen = 1;
while gen<=max_gen && max(fit) < chrom_len
P = tournament_selection(P,fit,2);
P = two_point_crossover(P,pc);
P = point_mutation(P,pm);
fit = fitness(P);
max_fit(gen) = max(fit);
mean_fit(gen) = mean(fit);
disp_info(gen,max_fit(gen));
gen = gen+1;
end
disp(sprintf('Generation: %d',gen));
disp(sprintf('Best fitness: %d\n',max_fit(end)));
[m,ind] = max(fit);
best = P(ind,:);
end
function [P] = initialize(pop_size,chrom_length)
P = round(rand(pop_size,chrom_length));
end
function [fit] = maxones_fitness(P)
for i=1:size(P,1)
fit(i) = length(find(P(i,:)));
end
end
function [P_new] = tournament_selection(P,fit,tourn_size)
for i=1:size(P,1)
t = ceil(rand(1,tourn_size)*size(P,1));
[max_fit,winner] = max(fit(t));
P_new(i,:) = P(t(winner),:);
end
end
function [P_new] = two_point_crossover(P,pc)
mating_list = randperm(size(P,1));
P_new = [];
while ~isempty(mating_list)
pair = mating_list(1:2);
mating_list(1:2) = [];
if rand<pc
crossover_points = ceil(rand(1,2)*(size(P,2)));
point1 = min(crossover_points);
point2 = max(crossover_points);
individual1 = P(pair(1),:);
individual2 = P(pair(2),:);
individual1(point1:point2) = P(pair(2),point1:point2);
individual2(point1:point2) = P(pair(1),point1:point2);
P_new = [P_new;individual1;individual2];
else
P_new = [P_new;P(pair,:)];
end
end
end
function [P_new] = point_mutation(P,pm)
r = rand(size(P));
mutation_list = find(r<pm);
P_new = P;
P_new(mutation_list(find(P(mutation_list)==1))) = 0;
P_new(mutation_list(find(P(mutation_list)==0))) = 1;
end
function [] = disp_info(gen,fit)
if mod(gen,10)==0
disp(sprintf('Generation: %d',gen));
disp(sprintf('Best fitness: %d\n',fit));
end
end
function [P_new] = roulette_selection(P,fit)
fit = (fit - min(fit)).^2;
fit = cumsum(fit);
fit = fit/max(fit);
P_new = [];
for i=1:size(P,1)
f = find(fit>rand);
P_new = [P_ne';P(f(1),:)];
end
end
function fit = fitness(P)
fit = zeros(size(P,1),1);
for i=1:size(P,1)
[beta_c, theta, psi, E0, A0, I0] = gene_to_values(P(i,:));
[t,y] = ode15s(@covidmodel1, [1,day(end)],[14368332,E0,A0,I0,1,0,0]);
error = squared_error(day,Totalcases,t,y(:,7));
%error = squared_error(day,Activecases,t,y(:,5));
% error = squared_error(day,Totalcases,Activecases,t,y(:,7),y(:,5)); % 2 dataset
fit(i) = error^-1;
end
end
function [beta_c, theta, psi, E0, A0, I0] = gene_to_values(gene)
beta_c_min_exponent = -7;%-7;
beta_c_max_exponent = 1;
theta_min_exponent = -3;
theta_max_exponent = 1;
psi_min_exponent = -3;%-7;
psi_max_exponent = 1;
E0_min_exponent = 0;
E0_max_exponent = 2;
A0_min_exponent = 0;
A0_max_exponent = 2;
I0_min_exponent = 0;
I0_max_exponent = 2;
gene_len = size(gene, 2)/2;
beta_c_decimal = bi2de(gene(1:gene_len));
theta_decimal = bi2de(gene((gene_len+1):end));
psi_decimal = bi2de(gene((gene_len+1):end));
E0_decimal = bi2de(gene((gene_len+1):end));
A0_decimal = bi2de(gene((gene_len+1):end));
I0_decimal = bi2de(gene((gene_len+1):end));
beta_c_exponent = beta_c_min_exponent +...
((beta_c_max_exponent-beta_c_min_exponent)*...
beta_c_decimal/2^gene_len);
theta_exponent = theta_min_exponent +...
((theta_max_exponent-theta_min_exponent)*...
theta_decimal/2^gene_len);
psi_exponent = psi_min_exponent +...
((psi_max_exponent-psi_min_exponent)*...
psi_decimal/2^gene_len);
E0_exponent = E0_min_exponent +...
((E0_max_exponent-E0_min_exponent)*...
E0_decimal/2^gene_len);
A0_exponent = A0_min_exponent +...
((A0_max_exponent-A0_min_exponent)*...
A0_decimal/2^gene_len);
I0_exponent = I0_min_exponent +...
((I0_max_exponent-I0_min_exponent)*...
I0_decimal/2^gene_len);
beta_c = 10^beta_c_exponent;
theta = 10^theta_exponent;
psi = 10^psi_exponent;
E0 = 10^E0_exponent;
A0 = 10^A0_exponent;
I0 = 10^I0_exponent;
end
%=========================================================================
% Fitting 1 dataset
function error = squared_error(data_time,data_points,fn_time,fn_points)
fn_values = interp1(fn_time,fn_points,data_time);
error = sum((fn_values - data_points).^2);
end
%==========================================================================
% % Fitting 2 datasets
% function error = squared_error(data_time,data_points1,data_points2,fn_time,fn_points1,fn_points2)
%
% fn_values1 = interp1(fn_time,fn_points1,data_time); fn_values2 = interp1(fn_time,fn_points2,data_time);
%
% error = sum((fn_values1 - data_points1).^2) + sum((fn_values2 - data_points2).^2);
% end
%=========================================================================
function out = bi2de(in)
% binary to decimal conversion
% see FEX : https://fr.mathworks.com/matlabcentral/answers/470466-most-efficient-way-to-convert-from-binary-to-dec-without-bin2dec-or-bi2de
out = sum(in.*(2.^(size(in,2)-1:-1:0)),2);
end
end
  2 Comments
Mathieu NOE
Mathieu NOE on 27 Aug 2024
Hmm Houston we have a problem
while it seems that our bi2de function gives similar results (tested on one example only I admit) as the regular bin2dec , the overall result in this code gives strange results :
so , a quick comparison with the other option below (with bin2dec) shows that the latter seems to give the appropriate results :
- some further investigation about the other post - have to check if that is indeed linked to the limitations of 52 bits as highlighted by John D'Errico : With only the minor caveat that if you have more than 52 bits, it will fail, since the result will be in the form of a double.
so please use this code below instead of my first answer above !
% Data fitting for the covid-19 outbreak in Nigeria, with reference to
% Lagos State Using the Genetic Algorithm (GA).
function covid19_GA2
clear all
clc
global x1 beta_c ep_g c_g sigma nu theta psi d_o d_d gamma_i alpha ...
delta ep_s ep_i ep_a best E0 A0 I0
ep_g = 0; c_g = 0; delta = 0; ep = 0; ep_s = ep; ep_a=ep; ep_i=ep;
sigma = 1/5.2; nu = 0.5;
gamma_i = 1/15;
gamma_a = 0.13978; gamma_o = 0.13978;
d_o = 0.015; d_d = 0.015;
alpha = 0.5;
%==========================================================================
% Lagos data, starting from March 16, 2020 till May 3, 2020.
day = [1 3 4 6 7 8 9 10 11 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 ...
43 44 45 46 47 48];
Totalcases = [1 5 9 19 22 25 29 32 44 59 68 81 82 91 98 109 109 120 120 130 145 158 166 174 176 189 214 232 251 283 ...
306 376 376 430 504 582 657 689 731 764 844 931 976 1006 1068];
Activecases = [1 5 8 19 22 25 29 32 43 58 67 75 76 85 81 86 86 91 87 96 111 117 118 119 116 123 139 140 154 182 200 ...
266 266 309 383 460 525 548 584 602 682 718 756 753 791];
%==========================================================================
pop_size = 100; chrom_len = 100; pm = 0.01; pc = 0.5; max_gen = 50;
%suggested run: [P,best] = ga(100,100,0.01,0.5,200)
[P,best] = gaCovid(pop_size,chrom_len,pm,pc,max_gen);
[beta_c,theta,psi, E0, A0, I0] = gene_to_values(best);
%Control Reproduction Number
R_c = beta_c*(1-ep_g*c_g)*(1-delta)*(1-ep_s)*((alpha*nu*(1-ep_a)/...
(theta+gamma_a)) + ((1-nu)*(1-ep_i)/(psi+d_o+gamma_o)))
[t,y] = ode15s(@covidmodel1, [0,52],[14368332,E0,A0,I0,1,0,0]);
plot(t,y(:,7),'-b',day,Totalcases,'rs')
xlabel('Time (days)'),ylabel('Cumulative number of reported cases')
legend('Model','Lagos COVID-19 Data')
% plot(t,y(:,5),'-b',day,Activecases,'rs')
% subplot(211), plot(t,y(:,7),'-b',day,Totalcases,'rs')
% subplot(212), plot(t,y(:,5),'-b',day,Activecases,'rs')
function dx = covidmodel1(t,x)
dx = [0;0;0;0;0;0;0];
N = x(1)+x(2)+x(3)+x(4)+x(5)+x(6);
beta = beta_c*(1-ep_g*c_g)*(1-delta)*(1-ep_s);
dx(1) = -beta*x(1)*(alpha*(1-ep_a)*x(3)+(1-ep_i)*x(4))/(N-x(5));
dx(2) = beta*x(1)*(alpha*(1-ep_a)*x(3)+(1-ep_i)*x(4))/(N-x(5)) - sigma*x(2);
dx(3) = nu*sigma*x(2) - theta*x(3) - gamma_a*x(3);
dx(4) = (1-nu)*sigma*x(2) - psi*x(4) - d_o*x(4) - gamma_o*x(4);
dx(5) = theta*x(3) + psi*x(4) - gamma_i*x(5) -d_d*x(5); % Active cases, detected (asymptomatic and symptomatic)
dx(6) = gamma_i*x(5) + gamma_a*x(3) + gamma_o*x(4);
dx(7) = theta*x(3) + psi*x(4); % Cumulative number of new cases reported (asymptomatic and symptomatic)
end
function [P,best] = gaCovid(pop_size,chrom_len,pm,pc,max_gen)
P = initialize(pop_size,chrom_len);
fit = maxones_fitness(P);
gen = 1;
while gen<=max_gen && max(fit) < chrom_len
P = tournament_selection(P,fit,2);
P = two_point_crossover(P,pc);
P = point_mutation(P,pm);
fit = fitness(P);
max_fit(gen) = max(fit);
mean_fit(gen) = mean(fit);
disp_info(gen,max_fit(gen));
gen = gen+1;
end
disp(sprintf('Generation: %d',gen));
disp(sprintf('Best fitness: %d\n',max_fit(end)));
[m,ind] = max(fit);
best = P(ind,:);
end
function [P] = initialize(pop_size,chrom_length)
P = round(rand(pop_size,chrom_length));
end
function [fit] = maxones_fitness(P)
for i=1:size(P,1)
fit(i) = length(find(P(i,:)));
end
end
function [P_new] = tournament_selection(P,fit,tourn_size)
for i=1:size(P,1)
t = ceil(rand(1,tourn_size)*size(P,1));
[max_fit,winner] = max(fit(t));
P_new(i,:) = P(t(winner),:);
end
end
function [P_new] = two_point_crossover(P,pc)
mating_list = randperm(size(P,1));
P_new = [];
while ~isempty(mating_list)
pair = mating_list(1:2);
mating_list(1:2) = [];
if rand<pc
crossover_points = ceil(rand(1,2)*(size(P,2)));
point1 = min(crossover_points);
point2 = max(crossover_points);
individual1 = P(pair(1),:);
individual2 = P(pair(2),:);
individual1(point1:point2) = P(pair(2),point1:point2);
individual2(point1:point2) = P(pair(1),point1:point2);
P_new = [P_new;individual1;individual2];
else
P_new = [P_new;P(pair,:)];
end
end
end
function [P_new] = point_mutation(P,pm)
r = rand(size(P));
mutation_list = find(r<pm);
P_new = P;
P_new(mutation_list(find(P(mutation_list)==1))) = 0;
P_new(mutation_list(find(P(mutation_list)==0))) = 1;
end
function [] = disp_info(gen,fit)
if mod(gen,10)==0
disp(sprintf('Generation: %d',gen));
disp(sprintf('Best fitness: %d\n',fit));
end
end
function [P_new] = roulette_selection(P,fit)
fit = (fit - min(fit)).^2;
fit = cumsum(fit);
fit = fit/max(fit);
P_new = [];
for i=1:size(P,1)
f = find(fit>rand);
P_new = [P_ne';P(f(1),:)];
end
end
function fit = fitness(P)
fit = zeros(size(P,1),1);
for i=1:size(P,1)
[beta_c, theta, psi, E0, A0, I0] = gene_to_values(P(i,:));
[t,y] = ode15s(@covidmodel1, [1,day(end)],[14368332,E0,A0,I0,1,0,0]);
error = squared_error(day,Totalcases,t,y(:,7));
%error = squared_error(day,Activecases,t,y(:,5));
% error = squared_error(day,Totalcases,Activecases,t,y(:,7),y(:,5)); % 2 dataset
fit(i) = error^-1;
end
end
function [beta_c, theta, psi, E0, A0, I0] = gene_to_values(gene)
beta_c_min_exponent = -7;%-7;
beta_c_max_exponent = 1;
theta_min_exponent = -3;
theta_max_exponent = 1;
psi_min_exponent = -3;%-7;
psi_max_exponent = 1;
E0_min_exponent = 0;
E0_max_exponent = 2;
A0_min_exponent = 0;
A0_max_exponent = 2;
I0_min_exponent = 0;
I0_max_exponent = 2;
gene_len = size(gene, 2)/2;
beta_c_decimal = bin2dec(join(string(gene(1:gene_len))));
theta_decimal = bin2dec(join(string(gene((gene_len+1):end))));
psi_decimal = bin2dec(join(string(gene((gene_len+1):end))));
E0_decimal = bin2dec(join(string(gene((gene_len+1):end))));
A0_decimal = bin2dec(join(string(gene((gene_len+1):end))));
I0_decimal = bin2dec(join(string(gene((gene_len+1):end))));
beta_c_exponent = beta_c_min_exponent +...
((beta_c_max_exponent-beta_c_min_exponent)*...
beta_c_decimal/2^gene_len);
theta_exponent = theta_min_exponent +...
((theta_max_exponent-theta_min_exponent)*...
theta_decimal/2^gene_len);
psi_exponent = psi_min_exponent +...
((psi_max_exponent-psi_min_exponent)*...
psi_decimal/2^gene_len);
E0_exponent = E0_min_exponent +...
((E0_max_exponent-E0_min_exponent)*...
E0_decimal/2^gene_len);
A0_exponent = A0_min_exponent +...
((A0_max_exponent-A0_min_exponent)*...
A0_decimal/2^gene_len);
I0_exponent = I0_min_exponent +...
((I0_max_exponent-I0_min_exponent)*...
I0_decimal/2^gene_len);
beta_c = 10^beta_c_exponent;
theta = 10^theta_exponent;
psi = 10^psi_exponent;
E0 = 10^E0_exponent;
A0 = 10^A0_exponent;
I0 = 10^I0_exponent;
end
%=========================================================================
% Fitting 1 dataset
function error = squared_error(data_time,data_points,fn_time,fn_points)
fn_values = interp1(fn_time,fn_points,data_time);
error = sum((fn_values - data_points).^2);
end
%==========================================================================
% % Fitting 2 datasets
% function error = squared_error(data_time,data_points1,data_points2,fn_time,fn_points1,fn_points2)
%
% fn_values1 = interp1(fn_time,fn_points1,data_time); fn_values2 = interp1(fn_time,fn_points2,data_time);
%
% error = sum((fn_values1 - data_points1).^2) + sum((fn_values2 - data_points2).^2);
% end
%=========================================================================
end

Sign in to comment.

Categories

Find more on Biological and Health Sciences in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!