%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Baldwin effect by JV Stone, January 2010.
% Copyright: This software is licensed under the Creative Commons License
% (CC-GNU GPL version 2.0 or later). 
%     Contact: j.v.stone at shef.ac.uk
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copied from main_evol2.m 10/3/05.
% Implements the Baldwin effect as described in Hinton's paper.
% G. E. Hinton and S. J. Nowlan. How learning can guide evolution. Complex
% Systems, 1:495--502, 1987.

% code [1 2 3 ]==[1 0 ?]
clear all;
% jrand_set_seed(9);
ngene 		= 20; 	% num genes
norg		= 1000; % num organisms
proplearn 	= 0.5; 	% proportion of genes unspecified (ie state = ?)
pop = ones(ngene,norg); % one org in each col.
nsetgene = proplearn*ngene;

ngen = 250; % num generations.
ntrial = 1000; ntrial0=ntrial; % num learning trials.
nmating	= norg; % num matings

% Init state of each organism
pop(1:nsetgene/2,:)=2;
pop(nsetgene+1:ngene,:)=3; % -1 => ? (not set)
f(1);imagesc(pop); title('Proportion 0 1 and ?'); 

%%%%%%
% Use this
%pop = jshuffle_each_col(pop);
% or
% function a=jshuffle_each_col(a);
a=pop;
[r c]=size(a);
for i=1:c
	newr=randperm(r);
	col = a(:,i);
	a(:,i)=col(newr,:);
end;
pop=a;
%%%%%%%

f(2);imagesc(pop);title('Population, one org per column');

stats=zeros(ngen,3);

for g=1:ngen
	orgscores = zeros(norg,1); % how close to solution is org?
	orgprobs = zeros(norg,1);	% prob of mating for each org.

	for i=1:norg		
		org = pop(:,i);
		org0=org;
		% check if org has enuf ones to do task in principle ...
		if sum(org==1)<proplearn*ngene ntrial=1;else ntrial=ntrial0; end;
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%		
		% Learn
		maxscore=0;
		ind = find(org==3); % indices of switchablegenes
		nsw = length(ind);
		for t=1:ntrial
			% Set switchable genes to 1 or 2.
			switchablegenes = rand(nsw,1)>0.5;
			ind2=(switchablegenes==0);
			switchablegenes(ind2)=2;
			org(ind) = switchablegenes;
			% [org0 org] pr;
			% Score genome
			score = sum(org==1);
			if score>maxscore
				maxscore = score;
				% bestorg = org;
			end;
			if score==ngene break; end;
		end; % learning trials
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		% Record score and best org.
		orgscores(i) = maxscore;
		orgprobs(i) = 1 + (ngene-1)*(ntrial-t)/ntrial;
		% orgprobs(i) = (ntrial-t)/ntrial;
		% [t orgprobs(i)] pr;
		% [orgscores(i) orgprobs(i)]
		% NO - keep pop same! pop(:,i) = bestorg;
	end; % org
	% Assign prob to each org
	% fprintf('%e\n',orgprobs);
	k=sum(orgprobs); 
	if k==0 
		orgprobs=ones(size(orgprobs))/norg; 
	else
		orgprobs = orgprobs/k;
	end;
	orgcumprobs = cumsum(orgprobs);
	% Make each org have up to 
	f(1); hist(orgscores);  title('scores');
	f(2); plot(orgcumprobs); title('orgprobs'); drawnow;%pr;
	
	% Mating
	newpop=pop;
	for m=1:nmating
		% choose parent orgs
		p=rand; b=(orgcumprobs<p);	orgid1=sum(b)+1;
		p=rand; b=(orgcumprobs<p);	orgid2=sum(b)+1;
		org1=pop(:,orgid1);
		org2=pop(:,orgid2);
		% Make child
		mask=rand(ngene,1)>0.5;
		child1 = org1.*mask;
		child2 = org2.*~mask;
		child=child1+child2;
		% [child1 child2 child] pr;
		newpop(:,m)=child;
	end;
	pop=newpop;

	% imagesc(pop);pr;
	% Count num of each allele type.
	p=pop(:);
	k=norg*ngene;
	n1=sum(p==1)/k;
	n0=sum(p==2)/k;
	nm1 = sum(p==3)/k;
	stats(g,1)=nm1;	stats(g,2)=n0;	stats(g,3)=n1;
	f(3);plot(stats(1:g,:));legend('?','0','1');title('Gene frequency');drawnow;
end; % gen

stats(1:g,:)

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%