format short;
clc;
clear;
close all;

%%
%Problem Parameters & objective
Gamma = linspace(0,1,50);
returns = @(x) means * x';
risks = @(x) x*Q*x';
nVar=31;
VarSize=[1 nVar];
VarMin = zeros(VarSize);
VarMax = ones(VarSize);
CostFunction = @(x,gam) gam * risks(x) - (1-gam) * returns(x);
%%
%readdata
dat = readmatrix("mean-var.txt"); %change with mean-var2.txt, etc.
dat2 = readmatrix("rho.txt"); %change with rho2.txt, etc.
dat3 = readmatrix("portef1.txt"); %change with portef2.txt, etc.
umeans = dat3(:,1);
urisk = dat3(:,2);
means = dat(:,1)';
sd = dat(:,2)';
rho = dat2(:,3);
Q = zeros([nVar nVar]);
for i = 1:nVar
    for j = i:nVar
        Q(i,j) = rho(round((i-1)*(2*nVar+2-i)/2)+j-i+1)*sd(i)*sd(j);
        Q(j,i) = Q(i,j);
    end
end

%%
%BBO parameters
maxgen= 50;
MaxIt=1500*nVar;
nPop=100;
nKeep=5;
nNew=nPop-nKeep;
lambda=linspace(1,0,nPop+1);
mu=1-lambda;
temp = 0;
pMax = 0.05;
for i = 1:nPop
    temp = temp + prod(lambda(1:i))/prod(mu(2:i+1));
end
prob0 = 1/(temp+1);
for i = 1:nPop
    prob(i) = prob0 * prod(lambda(1:i))/prod(mu(2:i+1));
end
probMax = max(prob);
for i = 1:nPop
    pMutation(i) = pMax * (1-prob(i)/probMax);
end
mu=linspace(1,1/nPop,nPop);
lambda = 1-mu;
%%
%main loop
habitat.Position=[];
habitat.Cost=[];
tic
for iter = 1:maxgen
CostF =  @(x)CostFunction(x,Gamma(iter));
pop=repmat(habitat,nPop,1);
for i=1:nPop
    pop(i).Asset = randperm(nVar,10);
    pop(i).Position = zeros(VarSize);
    pop(i).Position(pop(i).Asset) = rand([1 10]);
    [pop(i).Position, pop(i).Cost,~] = modify(pop(i),CostF);
end
[~, SortOrder]=sort([pop.Cost]);
pop=pop(SortOrder);
for it=1:MaxIt
    EP=lambda(2:10);
    EP(i)=0;
    EP=EP/sum(EP);
    j=RouletteWheelSelection(EP);
    newpop.Position = zeros(VarSize);
    j = j+1;
    for k = 1:nVar
    if pop(1).Position(k)>0 || pop(j).Position(k)>0
    newpop.Position(k) = 0.9 * pop(1).Position(k) + 0.1 * pop(j).Position(k);
    end
    if rand < pMutation(j)
        newpop.Position(k) = rand;
    end
    end
    [~,idx3] = find(newpop.Position>0);
    newpop.Asset = idx3;
    idx = pop(1).Asset;
    idx2 = pop(j).Asset;
    idx4 = setdiff(union(idx,idx2),idx3);
    numR = length(idx3);
    if numR > 10
        idx1 = randperm(numR,numR-10);
        newpop.Position(idx3(idx1)) = 0;
        newpop.Asset = setdiff(newpop.Asset, idx3(idx1));
    elseif numR < 10
        numA = length(idx4);
        idx1 = randperm(numA, 10-numR);
        newpop.Position(idx4(idx1)) = 0;
        newpop.Asset = union(newpop.Asset,idx4(idx1));
    end
    [newpop.Position, newpop.Cost,~] = modify(newpop, CostF);
    pop(100) = newpop;
    [~, SortOrder]=sort([pop.Cost]);
    pop=pop(SortOrder);
end
[~,~,w] = modify(pop(1),CostF);
optx(iter) = risks(w);
opty(iter) = returns(w);
disp(iter)
end
runtime=toc
%%
metric
plotcef