function ncp(testproblem)
% This function optimizes the reliability of the network, subject to
% constraints on the budget. 
%
% Inputs: a     - probability vector of purchasing links   
%         N1    - sample size used in estimating  reliability at end
%         c     - link costs                                   
%         Cmax  - total budget                             
%         rho   - rarity parameter
%         p     - reliability of links                     
%         alpha - smoothing parameter
%         N     - sample size                              
%         beta  - stopping constant
%         K     - number of samples used in estimation
%
% Outputs: a    - probability vector from final iteration    
%          time - CPU time
%          rhat - estimated reliability of a system          
%          RE   - relative error

for iter=1:1
iter

if nargin < 1
    %model parameters
    n=6; % number of nodes
    m=(n-1)*n/2; % number of links
    rand('state',1234);
    c=uniformrand(0.2*ones(1,m),0.3*ones(1,m),1); % link costs
    Cmax=1; % total budget
    c(n-1) = Cmax+0.3;   %make the direct link  too expensive
    p=uniformrand(0.6*ones(1,m),0.7*ones(1,m),1); % link reliabilities    
    a=0.5*ones(1,m); % initial probability vector
    N=300; N1=1000; 
    rho=0.1; % rarity parameter
    alpha=0.7; % smoothing parameter
    beta=0.02; %stopping parameter
    K=50; % number of samples used to estimate reliability
else
    load(testproblem);
end

rand('state',1234567 + 100*iter);
m=size(a,2); n=(1+sqrt(1+8*m))/2; idx=ceil((1-rho)*N);
P=vect2mat(p);
t=0; %interation counter
table=[]; 
prob=zeros(N,1); 
pn=zeros(N,1);
fid=fopen('result.txt','aw'); % refer to result.txt for results.
fprintf(fid,'n: %d  N: %d  N1: %d  K: %d  Cmax: %3.2f  rho: %1.2f  alpha: %1.2f  beta: %1.2f\n',n,N,N1,K,Cmax,rho,alpha,beta);
%fprintf(fid,'%1.4f & %1.4f & %1.4f & %1.4f & %1.4f & %1.4f & %1.4f & %1.4f & %1.4f & %1.4f\\\\\n',c);
%fprintf(fid,'%1.4f & %1.4f & %1.4f & %1.4f & %1.4f & %1.4f & %1.4f & %1.4f & %1.4f & %1.4f\\\\\n',p);

tic;

%%%%% Main CE loop
while max(min(a,1-a))> beta
    t=t+1;
    pv=uniperm(N,a,c,Cmax); % generate a random sample using uniform permutation
    prob=estrel1(pv,p,K);
    sortSX=sortrows([pv prob],m+1); % sort from smallest to largest0
    gamhat=sortSX(idx,m+1); % compute (1-rho) sample quantile
    table=[table; gamhat];
    x=sortSX(idx:N,1:m);
    a1=mean(x); % update probability vector a
    a=alpha.*a1+(1-alpha).*a; % smoothed updating of a 
    disp([t gamhat a]) % display current results
end
time=toc

x=round(a); pn=sum(2.^(m-find(x)));
network(x);

pv2=uniperm(N1,a,c,Cmax); prob1=zeros(N1,1);
prob1=estrel1(pv2,p,K);
rhat=mean(prob1); % Estimate reliability of a network 
RE=sqrt((1-rhat)/(N1*rhat)); % compute relative error

disp(['rhat:  ',num2str(rhat)])
disp(['RE: ',num2str(RE)])
fprintf('95 perc. CI : (%g,%g)\n', rhat - rhat*RE*1.96, rhat+rhat*RE*1.96)
fprintf(fid,'T: %d  pn: %-16.0f  CPU time: %2.3f  rhat: %2.5f  RE: %2.5f  lower bound: %2.5f  upper bound: %2.5f\n',t,pn,time,rhat,RE,rhat-rhat*RE*1.96,rhat+rhat*RE*1.96);
fprintf(fid,'\n');
fclose(fid);

beep
end
