% clear all
close all
tall = tic;

% a1_init=2.4; % use for MoC order 3
% b1_init=2.58; % use for MoC order 3

a1_init=2.4; % use for MoC order 4
b1_init=2.58; % use for MoC order 4


ra1= (b1_init-a1_init).*rand(40,1) + a1_init;

lam1=ra1;
N=3;
order=3;
numBr = 2^order - 1;

rimax = [ 0.296 0.318 0.337 0.358 0.384 0.414 0.445 0.484 0.539 0.608 ...
    0.692 0.793 0.913 1.052 1.203 1.374];

rinitial=[];

for k=1:order
    rinitial=[rinitial;rimax(k)*ones(2.^(order-k),1)];
end
% % % 
for jj=1:N
pert(:,jj)=0.01*rand(1,length(rinitial))';
rinitial1(:,jj)=rinitial+pert(:,jj);
end

psi1=zeros(N,length(lam1));

for j=1:N

parfor i=1:length(lam1)

[~,~,~,~,~,~, summ(i)]=coup(order,lam1(i),rinitial1(:,j));

end

     psi1(j,:) = summ;
%      parsave(sprintf('output%d.mat', j),psi(j,:));
end
lam_dm=repmat(lam1',N,1);
histogram2(psi1,lam_dm,40,'XBinLimits',[0,1],'YBinLimits',[2.4 2.58],'FaceColor','g'); view(2)
%%
% A = [0];
% B=[0];
% A1=[0];
% B1=[0];
% 
% f=zeros(N,length(lam1));
% for j=1:N
% for i=1:length(lam1)
% if p_bot1(j,i)>25
% f(j,i)=1;
% elseif p_bot1(j,i) >=10 && p_bot1(j,i) <=25
%  f(j,i)=0.5;
% else f(j,i)=0;
% end
% end
% end
%lam_moc=repmat(lam1,100,1);
% % lam_dm=repmat(lam',N,1);
% for ii = 1:N
% idx= find(f(ii,:)==1);
% idx1=find(f(ii,:)==0.5);
% % b1=psi1(ii,idx)';
% % a1=lam_moc(ii,idx)';
% b=psi1(ii,idx)';
% b11=psi1(ii,idx1)';
% a=lam_moc(ii,idx)';
% a11=lam_moc(ii,idx1)';
% %      b = ones(randi([1,10],1,1), 1); % here you put your random-sized vector
%      if (size(A(:,1)) == 1) 
% %          A1=b11;
% %          B1=a11;
%          A = b;
%          B=a;
%      else  
%          if (length(b) <  length(A(:,1))) 
% 
% %              A1 = [A1 [b1; zeros(length(A1(:,1)) - length(b1),1)]];
% %              B1= [B1 [a1; zeros(length(B1(:,1)) - length(a1),1)]];
%              A = [A [b; zeros(length(A(:,1)) - length(b),1)]];
%              B= [B [a; zeros(length(B(:,1)) - length(a),1)]];
%          end
%          if (length(b) == length(A(:,1)))
% %                 A1 = [A1 b1];
% %                 B1= [B1 a1];
%              A = [A b];
%              B = [B a];
%          end
%          if (length(b) >  length(A(:,1)))
% 
% %                A1 = [[A1; zeros(length(b1) - length(A1(:,1)), length(A1(1,:)))] b1];
% %                B1 = [[B1; zeros(length(a1) - length(B1(:,1)), length(B1(1,:)))] a1];
%              A = [[A; zeros(length(b) - length(A(:,1)), length(A(1,:)))] b];
%              B = [[B; zeros(length(a) - length(B(:,1)), length(B(1,:)))] a];
%          end
%      end
% end
%      if (size(A1(:,1)) == 1) 
%          A1=b11;
%          B1=a11;
%          %A = b;
%          %B=a;
%      else  
%          if (length(b11) <  length(A1(:,1))) 
% 
% %              A1 = [A1 [b1; zeros(length(A1(:,1)) - length(b1),1)]];
% %              B1= [B1 [a1; zeros(length(B1(:,1)) - length(a1),1)]];
%              A1 = [A1 [b11; zeros(length(A1(:,1)) - length(b11),1)]];% creates matrix of size (length(b11) 2) 
%              B1= [B1 [a11; zeros(length(B1(:,1)) - length(a11),1)]]; %creates matrix of size (length(a11) N) 
%          end
%          if (length(b11) == length(A1(:,1)))
% %                 A1 = [A1 b1];
% %                 B1= [B1 a1];
%              A1 = [A1 b11];
%              B1 = [B1 a11];
%          end
%          if (length(b11) >  length(A1(:,1)))
% 
% %                A1 = [[A1; zeros(length(b1) - length(A1(:,1)), length(A1(1,:)))] b1];
% %                B1 = [[B1; zeros(length(a1) - length(B1(:,1)), length(B1(1,:)))] a1];
%              A1 = [[A1; zeros(length(b11) - length(A1(:,1)), length(A1(1,:)))] b11];
%              B1 = [[B1; zeros(length(a11) - length(B1(:,1)), length(B1(1,:)))] a11];
%          end
%      end
% end
% % 
% % % use [5 7.5] for DM order 6, use [8 15] for DM order 4
% h=histogram2(psi,lam_dm,40,'XBinLimits',[0,1],'YBinLimits',[5 7.5]);%,'Normalization','probability');
% h.FaceColor='Flat';
% % h.EdgeColor='black';
% view(2)
% hold on
% h1=histogram2(A',B',40,'XBinLimits',[0,1],'YBinLimits',[5 7.5],'FaceColor','r');
% % h1.FaceColor='Flat';
% view(2)
% hold on
% h2=histogram2(A1',B1',40,'XBinLimits',[0,1],'YBinLimits',[6 14],'FaceColor','g');
% % h2.FaceColor='Flat';
% view(2)
%%
 toc(tall)
% profile viewer
% 
% for i=1:11
% hold on
% % xlim([0 0.2])
% ylim([0 100])
% subplot(3,4,i)
% histogram(psi(:,i),20)
% title(['DM \lambda =', num2str(lam1(i))])
% end

% for i=1:11
% hold on
% ylim([0 100])
% subplot(3,4,i)
% histogram(psi1(:,i),11)
% title(['MoC \lambda =', num2str(lam1(i))])
% end
% 
% numBr = 2^order - 1;
% 
%     % Define the edges of the graph first
%     node1 = zeros(1,numBr);
%     node2 = zeros(1,numBr);
%     
% %     e.g. for an order 4: node1 = [1 2 2 3 3 4 4 5 5 6 6 7 7 8 8];
%     node1(1) = 1;
%     for i = 1:(numBr-1)/2
%         node1(2*i:2*i+1) = [i+1 i+1];
%     end
%    
%     % e.g. for an order 4: node2 = [2 3 4 5 6 7 8 9 10 11 12 13 14 15 16];
%     for i = 1:numBr
%         node2(i) = i+1;
%     end
% %     
% % %%    %    
 %
% axis tight manual
% ax = gca;
% ax.NextPlot = 'replaceChildren';
% F(length(lam)) = struct('cdata',[],'colormap',[]); 
% v = VideoWriter('DM.avi','Motion JPEG AVI');
%  v.Quality = 100;
%  v.FrameRate = 2;
%  open(v);
% for j=1:length(lam)
% %   v = VideoWriter('DM.avi');   
% % open(v); 
% 
% 
% lam(j)
% labels=num2cell(air(:,j))';
% labels_flip=fliplr(labels);
% A=fliplr(air(j,:));
% % Fl=fliplr(flow1_mean(:,j)'./(-50))';
% G = graph(node1,node2,A);
% figure(3);
% % u = uicontrol('Style','slider','Position',[10 50 20 340],...
% %     'Min',1,'Max',loops,'Value',1);
% h=plot(G,'LineWidth',5);
% 
% labelnode(h,[node1 node2],'')
% % nodes=[0 1 2 2 1 5 5];
% % treeplot(nodes);
% % c = get(gca, 'Children'); % get handles to children
% % grab X and Y coords from the second child (the first one is axes)
% 
% h.EdgeCData=A;
% 
% colormap jet
% 
% ind = zeros(1,numBr);
%     for i = 1:numBr
%         ind(i) = i;
%     end
%    
% labeledge(h,ind,G.Edges.Weight(ind));
% % u.Value = j;
% F(j) = getframe(gcf);
% % size(F(j).cdata);
% % writeVideo(v,F(j));
% % x = get(h, 'XData');
% % y = get(h, 'YData');
% % c = colorbar;
% %     c.Label.String = 'r';
% %     caxis([0,0.32]);
% % text(x, y, labels_flip, 'VerticalAlignment','bottom', ...
% %                          'HorizontalAlignment','right')
% % highlight(h,GF,'EdgeColor',[0.9 0.3 0.1],'NodeColor',[0.9 0.3 0.1])
% 
% 
% pause
% 
% end