function [x fval exitflag output lambda g h]=findMLGRAN41(X0,N,N1,T,Iyy,GX,GU,GY,GYinv,TR,LDET, ...
    DvGdRHO1,DvGdC,DvGdG,DvGdT1,D,Dinv,In2Knn,En,coslambda,cos2lambda,cos3lambda,cos4lambda);


In=sparse(eye(N));
% options = optimset('LargeScale','off','GradObj','off','display','iter','TolX',0.00000000001);
options = optimset('LargeScale','on','GradObj','on','Hessian','on','display','iter','TolCon',0.0,'TolFun',0.000000000000000000001);




%A=[eye(10); -eye(10)];
%b=[20 20 20 10 10 10 .99 .99 .99 .99 0 0 0 0 0 0 .99 .99 .99 .99]' ;
lb=[-10 -10 -10 -10    0.001 0.0001 0.0001 0.0001  -.98 -.98 -.98 -.98 -.98 -.98 -.98 -.98 -.98 -.98]';
ub=[ 10  10  10  10   50     50      50      50    .98  .98  .98  .98  .98  .98 .98  .98  .98  .98]';

[x fval exitflag output lambda g h]=fmincon(@loglik,X0,[],[],[],[],lb,ub,[],options);
%[x fval exitflag]=fmincon(@loglik,X0,A,b,[],[],[],[],[],options);
%[x fval exitflag]=fmincon(@loglik,X0,A,b,options);



    function [f,g,h] = loglik(x)
        


            Chat=x(1:N)';
            GAMhat=(x(N1:2*N)');
            RHO1hat=x(2*N+1);
            RHO2hat=x(2*N+2);
            RHOS1hat=x(2*N+2+1:3*N+2);
            RHOS2hat=x(3*N+2+1:4*N+2);
            PHI=[1 -RHO1hat -RHO2hat];
            Vf=1;
            Vi=zeros(N,1);
            for i=1:N,
                PHIS(i,:)=[1 -RHOS1hat(i) -RHOS2hat(i)];
                Vi(i)=GAMhat(i);
            end





            GXdenom=coj(PHI,PHI);
            GXdenom2=multh(GXdenom,GXdenom);
            GU1denom=coj(PHIS(1,:),PHIS(1,:));
            GU2denom=coj(PHIS(2,:),PHIS(2,:));
            GU3denom=coj(PHIS(3,:),PHIS(3,:));
            GU4denom=coj(PHIS(4,:),PHIS(4,:));
            h1=[-2*RHO1hat(1)    -2*(RHO2hat(1)-1)];
            h2=[-2*RHO2hat(1)    -2*RHO1hat(1)      2];
            k11=[-2*RHOS1hat(1) -2*(RHOS2hat(1)-1)];
            k12=[-2*RHOS1hat(2) -2*(RHOS2hat(2)-1)];
            k13=[-2*RHOS1hat(3) -2*(RHOS2hat(3)-1)];
            k14=[-2*RHOS1hat(4) -2*(RHOS2hat(4)-1)];
            k21=[-2*RHOS2hat(1)    -2*RHOS1hat(1)      2];
            k22=[-2*RHOS2hat(2)    -2*RHOS1hat(2)      2];
            k23=[-2*RHOS2hat(3)    -2*RHOS1hat(3)      2];
            k24=[-2*RHOS2hat(4)    -2*RHOS1hat(4)      2];
            CkC=kron(Chat,Chat);
            CC=Chat*Chat';

            GX=Vf./(GXdenom(1)+GXdenom(2)*coslambda+GXdenom(3)*cos2lambda)';
            GUU=repmat(Vi,1,T)./[GU1denom(1)+GU1denom(2)*coslambda+GU1denom(3)*cos2lambda; ...
                                 GU2denom(1)+GU2denom(2)*coslambda+GU2denom(3)*cos2lambda; ...
                                 GU3denom(1)+GU3denom(2)*coslambda+GU3denom(3)*cos2lambda; ...
                                 GU4denom(1)+GU4denom(2)*coslambda+GU4denom(3)*cos2lambda];
            GU=reshape(En'*GUU,N,N,T);


              GY=reshape((repmat(CC,T,1).*kron(GX,ones(N,N)))',N,N,T)+GU;


              
              for i=1:T
                  GYinv(:,:,i)=inv(GY(:,:,i));
                  TR(i)=trace(GYinv(:,:,i)*Iyy(:,:,i));
                  LDET(i)=log(det(GY(:,:,i)));
              end
            f=-0.5*(sum(TR)+sum(LDET))-N*T/2*log(2*pi);

            % Gradient
            INFMAT=zeros(4*N+2,4*N+2);
            DFDP=zeros(4*N+2,1);
        

            for i=1:T
        
                M=kron(GYinv(:,:,i),GYinv(:,:,i).');
 
                GXdenom2l=GXdenom2(1)+GXdenom2(2)*coslambda(i)+GXdenom2(3)*cos2lambda(i)+GXdenom2(4)*cos3lambda(i) ...
                    +GXdenom2(5)*cos4lambda(i);
    %             GXdenom2lalt=(GX(i)/Vf)^2;
    %              GXdenom2l-GXdenom2lalt
    %              pause
                h1l=h1(1)+h1(2)*coslambda(i);
                h2l=h2(1)+h2(2)*coslambda(i)+h2(3)*cos2lambda(i);
                DvGdRHO1(:,i)=CkC*((Vf*h1l./GXdenom2l));
                DvGdRHO2(:,i)=CkC*((Vf*h2l./GXdenom2l));
                DvGdC(:,:,i)=(In2Knn)*kron(eye(N),Chat*GX(i));
                DvGdG(:,:,i)=En'*diag((diag(GU(:,:,i))./Vi));
                k11l=k11(1)+k11(2)*coslambda(i);
                k12l=k12(1)+k12(2)*coslambda(i);
                k13l=k13(1)+k13(2)*coslambda(i);
                k14l=k14(1)+k14(2)*coslambda(i);
                k1l=[k11l k12l k13l k14l]';
                k21l=k21(1)+k21(2)*coslambda(i)+k21(3)*cos2lambda(i);
                k22l=k22(1)+k22(2)*coslambda(i)+k22(3)*cos2lambda(i);
                k23l=k23(1)+k23(2)*coslambda(i)+k23(3)*cos2lambda(i);
                k24l=k24(1)+k24(2)*coslambda(i)+k24(3)*cos2lambda(i);
                k2l=[k21l k22l k23l k24l]';
 

                DvGdT1(:,:,i)=En'*diag((k1l.*diag(GU(:,:,i).^2)./Vi));
                DvGdT2(:,:,i)=En'*diag((k2l.*diag(GU(:,:,i).^2)./Vi));

                DGDP=[DvGdC(:,:,i) DvGdG(:,:,i)  DvGdRHO1(:,i) DvGdRHO2(:,i) DvGdT1(:,:,i) DvGdT2(:,:,i)];
                INFMAT=INFMAT+DGDP.'*M*conj(DGDP);
            DFDP=DFDP+DGDP.'*vec(GYinv(:,:,i).'*Iyy(:,:,i).'*GYinv(:,:,i).'-GYinv(:,:,i).');
            end
            h=INFMAT/2;
            g=-DFDP/2;
            f=-f;



        end



end


