function main2
neurogenesis(0);
end % main
function objfun = neurogenesis(run_num)
enrich = [1 2];
% enrich = [3 4];
% parameters
% output
TEXT_OUT = 1;
PLOT_OUT = 1; ANIME_OUT = 1;
PLOT2D_OUT = 0; ANIME2D_OUT = 0;
% stimulus
sim = 40;
odor_names = 'limonene(+)_ster limonene(_)_ster propylpropionate_es3 ethylbutyrate_es3 isopropylbenzene_ModuleC1 cyclohexanone_SG18 acetone methylacetate_SG19 cycloheptanelow_cycloalk propanol_simp_2500 isoamylbutyrate_est1 butyricacid_aci1 hexanal_ald1 ethylbenzene_HC';
choose = [enrich enrich 5:14];
% network
non_lin = 0;
conn = 4;
CS = 0.002;
% survival
ts = 0.1;
gamma = 5/ts;
th = 0.2;
rm = 0;
rg = 0;
% stepping
cont_density = 0;
exp_time = 5000; step = round(exp_time/100); dt = 5;
% testing
testing = 1;
test_Sstr = 1;
test_choose = [1 2]; Ns_test = 2;
% tracking vs T
tracking = 2; % 0: none, 1: in S, 2: in test S
track_pairs = [1 2];
% ...
S0 = 1; Sstr = 1; MC_per_Glom = 1;
TYPE = 1; % 1: pearson 2: L2
prob_conn = 0; % only for cont_density = 0
Na = 1;
perm_ratio = 0; % only for cont_density = 0
minv = 0; maxv = 1;
% setup stimulus
% setup stimulus
if TEXT_OUT == 1
fprintf('\n--- Start ---\n - initializing odor\n');
end
[Sall, coord, metric, name] = gara(14, MC_per_Glom, odor_names, sim);
Nc = size(Sall,1);
Ns = length(choose);
Senrich = S0+Sstr*(Sall(:,choose));
S = Senrich;
S(:,1:4) = 0;
S_name = cell(1,size(S,2));
for i = 1:size(S,2)
S_name{i} = name{choose(i)};
end
S_test = zeros(Nc,Ns_test);
for i = 1:Ns_test
S_test(:,i) = Sall(:,test_choose(i));
end
S_test = S0 + test_Sstr*S_test;
S_test_name = cell(1,size(S_test,2));
for i = 1:size(S_test,2)
S_test_name{i} = '';
end
if TEXT_OUT == 1
fprintf(' - Nc = %d\n', Nc);
fprintf(' corr = %f\n', mean_excluNaN(uptri_1d(corr(S,TYPE))));
end
if PLOT_OUT == 1;
setup_Pplot(S,corr(S,TYPE),corr(S',TYPE),rg,1);
if testing == 1
setup_Pplot(S_test,corr(S_test,TYPE), corr(S_test,TYPE),rg,1001);
end
drawnow;
end
if PLOT2D_OUT == 1;
setup_Pplot2D(S,coord,101,S_name);
if testing == 1
setup_Pplot2D(S_test,coord,1101,S_test_name);
end
drawnow;
end
% setup network
if TEXT_OUT == 1
fprintf(' - initializing network\n');
end
if cont_density == 1
option = odeset('Stats','off','RelTol',1e-3,'AbsTol',1e-8);
Isize = nchoosek(Nc,conn);
perm = nchoosek(1:Nc,conn);
C = zeros(Nc,Isize);
for i = 1:Isize
for c = 1:conn
C(perm(i,c),i) = 1;
end
end
N = rand(Isize,1);
Wmg = C*diag(N); Wgm = C';
Iage = ones(1,Isize);
Imark = zeros(1,Isize);
else
Isize = 2*Nc; N = ones(Isize,1);
Wmg = zeros(Nc, Isize); Wgm = zeros(Isize, Nc);
Iage = -ones(1,Isize);
Imark = zeros(1,Isize);
end
[P, I] = cal_activity(non_lin,CS,Wmg,Wgm,S,S,rm,rg);
time_axis = 0:step:exp_time;
N_t = NaN*ones(length(time_axis),Isize);
Pcorr_t = NaN*ones(1,length(time_axis));
Tcorr_t = NaN*ones(1,length(time_axis));
Pangle_t = NaN*ones(1,length(time_axis));
Tangle_t = NaN*ones(1,length(time_axis));
Pfoc_t = NaN*ones(1,length(time_axis));
CV_t = NaN*ones(1,length(time_axis));
CVid_t = NaN*ones(Ns,length(time_axis));
F_t = NaN*ones(Ns,Ns,length(time_axis));
if cont_density == 1
N_t(1,:) = N;
end
if PLOT_OUT == 1;
HP1d = setup_Pplot(P,corr(P,TYPE),corr(P',TYPE),rg,2,Wmg,Wgm,time_axis,N_t);
if testing == 1
[P_test, I_test] = cal_activity(non_lin,CS,Wmg,Wgm,S_test,S_test,rm,rg);
HPT1d = setup_Pplot(P_test,corr(P_test,TYPE),corr(P_test',TYPE),rg,1002);
end
[HI, Iaxis] = setup_Iplot(cont_density,time_axis,I,corr(I(Iage>=0,:),TYPE),Iage,Wmg,N_t(1,:),3);
VST = cell(1,2);
VST{1} = Pcorr_t; VST{2} = Tcorr_t;
line_style = cell(1,2);
line_style{1} = '-b'; line_style{2} = '-r';
line_name = cell(1,2);
line_name{1} = Pcorr_t; line_name{2} = Tcorr_t;
Hinfo = setup_Infoplot(time_axis,VST,line_style,line_name,corr(S,TYPE),corr(P,TYPE),4);
drawnow;
end
if PLOT2D_OUT == 1;
setup_Pplot2D(P,coord,102,S_name);
if testing == 1
[P_test, I_test] = cal_activity(non_lin,CS,Wmg,Wgm,S_test,S_test,rm,rg);
setup_Pplot2D(P_test,coord,1102,S_test_name);
end
drawnow;
end
% step
if TEXT_OUT == 1
fprintf(' - running\n');
end
Pcorr_t(1) = mean_excluNaN(uptri_1d(corr(P,TYPE)));
Pangle_t(1) = mean_excluNaN(uptri_1d(corr_angle(corr(S,TYPE),corr(P,TYPE))));
Pfoc_t(1) = mean_excluNaN(focality(P,metric));
CV_t(1) = std(mean(P,2))/mean(mean(P,2));
CVid_t(:,1) = std(P)./mean(P);
F_t(:,:,1) = corr(P);
if tracking == 1
Tcorr_t(1) = mean_excluNaN(cal_track_corr(track_pairs,P));
Tangle_t(1) = mean_excluNaN(corr_angle(cal_track_corr(track_pairs,S),cal_track_corr(track_pairs,P)));
elseif tracking == 2
[P_test, I_test] = cal_activity(non_lin,CS,Wmg,Wgm,S_test,S_test,rm,rg);
Tcorr_t(1) = mean_excluNaN(cal_track_corr(track_pairs,P_test));
Tangle_t(1) = mean_excluNaN(corr_angle(cal_track_corr(track_pairs,S_test),cal_track_corr(track_pairs,P_test)));
end
if TEXT_OUT == 1
fprintf(' corr = %f\n', Pcorr_t(1));
end
for i = 1:round(exp_time/step)
if i>round(exp_time/step)/2
S = Senrich;
setup_Pplot(S,corr(S,TYPE),corr(S',TYPE),rg,1);
end
if TEXT_OUT == 1
fprintf(' - run num = %d, time = %f\n',run_num,i*step);
end
if cont_density == 1
[ignore,N] = ode23(@RHS,[0 step],N_t(i,:),option);
N_t(i+1,:) = N(end,:);
Wmg = C*diag(N(end,:)); Wgm = C';
else
for j = 1:round(step/dt)
add_cell;
[P, I] = cal_activity(non_lin,CS,Wmg,Wgm,S,P,rm,rg);
remove_cell;
end
end
[P, I] = cal_activity(non_lin,CS,Wmg,Wgm,S,P,rm,rg);
Pcorr_t(i+1) = mean_excluNaN(uptri_1d(corr(P,TYPE)));
Pangle_t(i+1) = mean_excluNaN(uptri_1d(corr_angle(corr(S,TYPE),corr(P,TYPE))));
Pfoc_t(i+1) = mean_excluNaN(focality(P,metric));
CV_t(i+1) = std(mean(P,2))/mean(mean(P,2));
CVid_t(:,i+1) = std(P)./mean(P);
F_t(:,:,i+1) = corr(P);
if tracking == 1
Tcorr_t(i+1) = mean_excluNaN(cal_track_corr(track_pairs,P));
Tangle_t(i+1) = mean_excluNaN(corr_angle(cal_track_corr(track_pairs,S),cal_track_corr(track_pairs,P)));
elseif tracking == 2
[P_test, I_test] = cal_activity(non_lin,CS,Wmg,Wgm,S_test,P_test,rm,rg);
Tcorr_t(i+1) = mean_excluNaN(cal_track_corr(track_pairs,P_test));
Tangle_t(i+1) = mean_excluNaN(corr_angle(cal_track_corr(track_pairs,S_test),cal_track_corr(track_pairs,P_test)));
end
if TEXT_OUT == 1
fprintf(' corr = %f\n', Pcorr_t(i+1));
end
if ANIME_OUT == 1;
update_Pplot(P,corr(P,TYPE),corr(P',TYPE),rg,HP1d,Wmg,Wgm,N_t);
if testing == 1
[P_test, I_test] = cal_activity(non_lin,CS,Wmg,Wgm,S_test,P_test,rm,rg);
update_Pplot(P_test,corr(P_test,TYPE),corr(P_test',TYPE),rg,HPT1d);
end
update_Iplot(cont_density,time_axis,I,corr(I(Iage>=0,:),TYPE),Iage,Wmg,N_t(i,:),HI);
VST = cell(1,2);
VST{1} = Pcorr_t; VST{2} = Tcorr_t;
update_Infoplot(VST,corr(P,TYPE),Hinfo);
drawnow;
end
if ANIME2D_OUT == 1;
setup_Pplot2D(P,coord,102,S_name);
if testing == 1
[P_test, I_test] = cal_activity(non_lin,CS,Wmg,Wgm,S_test,S_test,rm,rg);
setup_Pplot2D(P_test,coord,1102,S_test_name);
end
drawnow;
end
end
% end
if PLOT_OUT == 1;
update_Pplot(P,corr(P,TYPE),corr(P',TYPE),rg,HP1d,Wmg,Wgm,N_t);
if testing == 1
[P_test, I_test] = cal_activity(non_lin,CS,Wmg,Wgm,S_test,P_test,rm,rg);
update_Pplot(P_test,corr(P_test,TYPE),corr(P_test',TYPE),rg,HPT1d);
end
update_Iplot(cont_density,time_axis,I,corr(I(Iage>=0,:),TYPE),Iage,Wmg,N_t(i,:),HI);
VST = cell(1,2);
VST{1} = Pcorr_t; VST{2} = Tcorr_t;
update_Infoplot(VST,corr(P,TYPE),Hinfo);
drawnow;
end
if PLOT2D_OUT == 1;
setup_Pplot2D(P,coord,102,S_name);
if testing == 1
[P_test, I_test] = cal_activity(non_lin,CS,Wmg,Wgm,S_test,S_test,rm,rg);
setup_Pplot2D(P_test,coord,1102,S_test_name);
end
drawnow;
end
objfun = return_val;
if TEXT_OUT == 1
fprintf('--- End ---\n');
end
% nested function definition
function dN_ = RHS(ignore,N1_)
N_ = N1_;
Wmg = C*diag(N_); Wgm = C';
[ignore,G_] = cal_activity(non_lin,CS,Wmg,Wgm,S,P,rm,rg);
P_ = survival(G_);
P_ = P_ + 1e-10*randn(size(P_));
B_ = Na*Nc/Isize/conn;
dN_ = B_ + log(P_).*N_;
end % RHS
function prob_ = survival(G_)
Ca_ = sum(rec(G_,th,0),2);
prob_ = (tanh((Ca_-ts)*gamma)+1)*(maxv-minv)/2+minv;
end % survival
function val = remove_cell
prob_ = survival(I);
surv_ = floor(prob_ + rand(size(prob_)));
IX_ = find((surv_==0)&(Iage'>=0));
Iage(IX_) = -1;
Imark(IX_) = 0;
Wmg(:,IX_) = 0;
Wgm(IX_,:) = 0;
val = length(IX_);
end % remove_cell
function val = add_cell
IX_ = find(Iage>=0);
Iage(IX_) = Iage(IX_)+dt;
IX_ = find(Iage<0);
if length(IX_) < round(dt*Na*Nc)
add_space(round(dt*Na*Nc)-length(IX_));
IX_ = find(Iage<0);
end
if prob_conn == 0
temp_ = [ones(conn,1); zeros(Nc-conn,1)];
for i_ = 1:round(dt*Na*Nc)
Iage(IX_(i_)) = 0;
temp_ = temp_(randperm(Nc));
Wmg(:, IX_(i_)) = temp_;
IX2_ = randperm(length(temp_));
temp1_ = temp_(IX2_(1:round(perm_ratio*length(temp_))));
temp2_ = temp_(IX2_(1+round(perm_ratio*length(temp_)):end));
Wgm(IX_(i_),IX2_) = [temp1_(randperm(length(temp1_))); temp2_]';
end
else
conn_prob_ = conn/Nc;
for i_ = 1:round(dt*Na*Nc)
Iage(IX_(i_)) = 0;
if (marking == 1) && (i*step >= marking_t(1)) && (i*step < marking_t(2))
Imark(IX_(i_)) = 1;
end
temp_ = floor(rand(Nc,1)+conn_prob_);
Wmg(:, IX_(i_)) = temp_;
IX2_ = randperm(length(temp_));
temp1_ = temp_(IX2_(1:round(perm_ratio*length(temp_))));
temp2_ = temp_(IX2_(1+round(perm_ratio*length(temp_)):end));
Wgm(IX_(i_),IX2_) = [temp1_(randperm(length(temp1_))); temp2_]';
end
end
val = round(dt*Na*Nc);
end % add_cell
function add_space(short)
previous_size = Isize;
tempI = I;
tempIage = Iage;
tempImark = Imark;
tempWmg = Wmg;
tempWgm = Wgm;
Isize = Isize + 5*short;
I = zeros(Isize, Ns);
Iage = -ones(1, Isize);
Imark = zeros(1, Isize);
Wmg = zeros(Nc, Isize);
Wgm = zeros(Isize, Nc);
I(1:previous_size, :) = tempI;
Iage(1, 1:previous_size) = tempIage;
Imark(1, 1:previous_size) = tempImark;
Wmg(:, 1:previous_size) = tempWmg;
Wgm(1:previous_size,:) = tempWgm;
if testing == 1
tempItest = I_test;
I_test = zeros(Isize, size(I_test,2));
I_test(1:previous_size, :) = tempItest;
end
if ANIME_OUT==1
set(Iaxis, 'YLim', [1 Isize]);
end
end % add_space
function val = cal_track_corr(track_pairs_,S_)
corr_ = zeros(1,size(track_pairs_,1));
for i_ = 1:size(track_pairs_,1)
temp_ = corr(S_(:,track_pairs_(i_,:)),TYPE);
corr_(i_) = temp_(1,2);
end
val = corr_;
end % cal_track_corr
function val = return_val
val = 0;
end % return_val
end % neurogenesis
