R/CircE.BFGS.R

Defines functions CircE.BFGS

Documented in CircE.BFGS

CircE.BFGS <-
function(R,
                        v.names,
                        
                        m,
                        N,r=1,
                        equal.com=FALSE,
                        equal.ang=FALSE, 
                        mcsc="unconstrained",     
                        start.values="IFA",
                        ci.level=0.95,
                        factr=1e9,pgtol=0,lmm=NULL,
                  iterlim=250, upper=NULL,lower=NULL,print.level=1,file=NULL,title="Circumplex Estimation",try.refit.BFGS=FALSE){
      if(!is.null(file))  sink(file,append=FALSE,split=TRUE)        
  

if(is.null(N)) stop('No default value for argument N. Specify sample size.')
if(is.null(m)) stop('No default value for argument m. Specify the number of free parameters in the Fourier correlation function (m>=1).')

if(m<=0) stop('The number of free parameters in the Fourier correlation function must be m>=1')

if (!is.null(lower)&equal.ang==TRUE) stop('You are trying to impose bounds on equally spaced polar angles! Set equal.ang=FALSE')
if (!is.null(upper)&equal.ang==TRUE) stop('You are trying to impose bounds on equally spaced polar angles! Set equal.ang=FALSE')
if(!is.null(upper) & !is.null(lower)){
for (i in 1:length(upper)){
if(upper[i]<lower[i]) stop ('lower bound greater than corresponding upper bound.')
} 
}


      p=dim(R)[1]
      is.triang <- function(R) {
      is.matrix(R) && (nrow(R) == ncol(R)) && 
            (all(0 == R[upper.tri(R)])) || (all(0 == R[lower.tri(R)]))
        }    
      is.symm <- function(R) {
        is.matrix(R) && (nrow(R) == ncol(R)) && all(R == t(R))
        }
    if (is.triang(R)) R <- R + t(R) - diag(diag(R))
    if (!is.symm(R)) stop('R MUST BE A SQUARE TRIANGULAR OR SYMMETRIC MATRIX !')
          
      k=3
      K=pi/180
cat("Date:",date(),"\n")
cat("Data:",title,"\n")
cat("Model:")
if(equal.com==TRUE & equal.ang==TRUE) {cat("Constrained model: equal spacing and equal radius \n")}
if(equal.com==TRUE & equal.ang==FALSE) {cat("Equal radius \n")}
if(equal.com==FALSE & equal.ang==TRUE) {cat("Equal spacing \n")}
if(equal.com==FALSE & equal.ang==FALSE){cat("Unconstrained model \n")}
cat("Reference variable at 0 degree:",v.names[r],"\n")
cat("\n")

ifa<-function(rr,mm) {
	if (length(which(eigen(rr)$values < 0)) != 0) {
            cat("WARNING!", "\n")
            cat("INPUT COVARIANCE/CORRELATION MATRIX IS NOT POSITIVE DEFINITE.", "\n")
            cat("STARTING VALUES CANNOT BE COMPUTED USING 'IFA': SET start.values='PFA'", "\n")
            stop("Make sure the listwise, not pairwise, missing data treatment has been selected in computing the input matrix\n")
             }

rinv <- solve(rr) 
sm2i <- diag(rinv)
smrt <- sqrt(sm2i)
dsmrt <- diag(smrt)
rsr <- dsmrt %*% rr %*% dsmrt
reig <- eigen(rsr)
vlamd <- reig$va
vlamdm <- vlamd[1:mm]
qqm <- as.matrix(reig$ve[, 1:mm])
theta <- mean(vlamd[(mm + 1):nrow(qqm)])
dg <- sqrt(vlamdm - theta)
fac <- diag(1/smrt) %*% qqm %*% diag(dg)
list(vlamd = vlamd, theta = theta,
fac = fac)
  }

pfa<-function (rr, mm =3, min.err = 0.001, max.iter = iterlim, 
    symmetric = TRUE, warnings = TRUE) 
{

        if (!is.matrix(rr)) {
            rr <- as.matrix(rr)
        }
        sds <- sqrt(diag(rr))
        rr <- rr/(sds %o% sds)
    
    
    
    rr.mat <- rr

    orig <- diag(rr)
    comm <- sum(diag(rr.mat))
    err <- comm
    i <- 1
    comm.list <- list()
    while (err > min.err) {
        eigens <- eigen(rr.mat, symmetric = symmetric)
        if (mm > 1) {
            loadings <- eigens$vectors[, 1:mm] %*% diag(sqrt(eigens$values[1:mm]))
        }
        else {
            loadings <- eigens$vectors[, 1] * sqrt(eigens$values[1])
        }
        model <- loadings %*% t(loadings)
        new <- diag(model)
        comm1 <- sum(new)
        diag(rr.mat) <- new
        err <- abs(comm - comm1)
        if (is.na(err)) {
            warning("imaginary eigen value condition encountered in factor.pa, exiting")
            break
        }
        comm <- comm1
        comm.list[[i]] <- comm1
        i <- i + 1
        if (i > max.iter) {
            if (warnings) {
                message("maximum iteration exceeded")
            }
            err <- 0
        }
    }
    if (!is.double(loadings)) {
        warning("the matrix has produced imaginary results -- proceed with caution")
        loadings <- matrix(as.double(loadings), ncol = mm)
    }
    if (mm > 1) {
        maxabs <- apply(apply(loadings, 2, abs), 2, which.max)
        sign.max <- vector(mode = "numeric", length = mm)
        for (i in 1:mm) {
            sign.max[i] <- sign(loadings[maxabs[i], i])
        }
        loadings <- loadings %*% diag(sign.max)
    }
    else {
        mini <- min(loadings)
        maxi <- max(loadings)
        if (abs(mini) > maxi) {
            loadings <- -loadings
        }
        loadings <- as.matrix(loadings)
    }
    loadings[loadings == 0] <- 10^-15
    model <- loadings %*% t(loadings)


list(fac=loadings)
}

        
start.valuesA=function(R,k,start.value=start.values){
        
        one=matrix(1,p,1)
        Lambda=if(start.value=="IFA")ifa(R,k)$fac else if(start.value=="PFA")pfa(rr=R, mm =k, min.err = 0.001, max.iter = iterlim)$fac

        Diagzeta=diag(diag(Lambda%*%t(Lambda)))
        Dzeta=sqrt(Diagzeta)
        
        uniq=diag(1,p,p)-(Lambda%*%t(Lambda))
        Dpsi=diag(diag(uniq))
        Dv=solve(Dzeta)%*%Dpsi
        
        Lambdax=solve(Dzeta)%*%Lambda
        Lambdaxhat=1/p*(t(Lambdax)%*%one)
        C=1/p*t(Lambdax-one%*%t(Lambdaxhat))%*%(Lambdax-one%*%t(Lambdaxhat))
        U=eigen(C)$vectors
        Lambdatilde=Lambdax%*%U
        
        if(equal.ang==FALSE){
                
        L..=Lambdatilde[,1:2]
        L..[,]=0
        for(i in 1:p){
                
                L..[i,]=Lambdatilde[i,1:2]/sqrt(Lambdatilde[i,1]^2+Lambdatilde[i,2]^2)
                
                }
   
   r=r
   
   sin.angles=rep(0,p)
   for(i in 1:p){
        
        sin.angles[i]=L..[i,2]*L..[r,1]-L..[i,1]*L..[r,2]
        
        }

   
   polar.angles=rep(0,p)
   for(i in 1:p){
        
        if(sin.angles[i]>=0){
        polar.angles[i]=L..[i,1]*L..[r,1]+L..[i,2]*L..[r,2]
        polar.angles[i]=acos(polar.angles[i])
        }
        else if(sin.angles[i]<0){
        polar.angles[i]=L..[i,1]*L..[r,1]+L..[i,2]*L..[r,2]
        polar.angles[i]=2*pi-acos(polar.angles[i])
        }
        }
  polar.angles=ifelse(polar.angles=="NaN",0,polar.angles) 
  polar.angles=polar.angles/K}

     if(equal.ang==TRUE){
           polar.angles=rep(0,p)
   for(i in 1:p){
        polar.angles[i]=(i-1)*(2*pi)/p
        }
        
polar.angles=polar.angles/K}
        

if(mcsc=="unconstrained"){
betas=matrix(0,1,m+1);
betas[,1]=(1/p*sum(Lambdatilde[,3]))^2
betas[,2]=1-betas[,1]
if(m>k){betas[,3:m]=0} ;if(m==k){ betas[,3]=0}
betas=betas/betas[,2]
                         }


if(mcsc=="-1"){
betas=matrix(0,1,m+1);
betas[,c(seq(1,(m+1),by=2))]=0
betas[,2]=1-betas[,1]
if(m>k){betas[,3:m]=0} ;if(m==k){ betas[,3]=0}
betas=betas/betas[,2] }



if(mcsc=="0"){
betas=matrix(0,1,m+1);
betas[,c(seq(1,(m+1)/2,by=2))]=1/(length(c(seq(1,(m+1)/2,by=2)))*2)
betas[,2]=1-betas[,1]
if(m>k){betas[,3:m]=0} ;if(m==k){ betas[,3]=0}
betas=betas/betas[,2] }


z=sqrt(diag(Lambda%*%t(Lambda)))
uniq=diag(1,p,p)-(diag(z^2))
        Dpsi=diag(diag(uniq))
        Dv=solve(Dzeta)%*%Dpsi
v=diag(Dv)
if(equal.com==TRUE){v=mean(v)} else v=v
if(equal.ang==FALSE)par=c(polar.angles[-c(1)]*K,c(betas[-c(2)]),v,z)
 if(equal.ang==TRUE) par=c(c(betas[-c(2)]),v,z)
 attributes(par)=list(parA=par,polar.angles=polar.angles,betas=betas)
 }
 
 parA=start.valuesA(R,k)$parA
 betas=start.valuesA(R,k)$betas
 
 
 ASK<-
function (arg) 
{
         value <- readline(paste("Continue? (YES/NO)", deparse(substitute(arg)), 
            ": "))
        if (value == "NO") 
           stop("STOP CURRENT COMPUTATION.",call.=FALSE)
     
}
if(length(parA)>1/2*(p*(p+1))){
	cat("ERROR: NUMBER OF PARAMETERS,",length(parA),", GREATER THAN NUMBER OF NONDUPLICATED ELEMENTS OF INPUT MATRIX,",1/2*(p*(p+1)),"\n* THE MODEL IS UNDERIDENTIFIED: Consider simplifying the model by reducing 'm' or by using equality constraints ('equal.com=TRUE' and/or 'equal.ang=TRUE'); \n* THE HESSIAN MATRIX MAY NOT BE POSITIVE DEFINITE:  THE STANDARD ERRORS OF THE MODEL PARAMETER ESTIMATES MAY NOT BE CALCULATED; \n* THE PROGRAM  MAY GENERATE NON-SENSICAL ESTIMATES OR DISPLAY VERY LARGE STANDARD ERRORS; \n* DEGREE OF FREEDOM < 0: SEVERAL FIT INDEXES CANNOT BE CALCULATED;...")
	ASK()}
if(length(parA)==1/2*(p*(p+1))){
	cat("ERROR: NUMBER OF PARAMETERS,",length(parA),", EQUAL TO THE NUMBER OF NONDUPLICATED ELEMENTS OF INPUT MATRIX,",1/2*(p*(p+1)),"\n* THE MODEL IS JUST IDENTIFIED (SATURATED): Consider simplifying the model by reducing 'm' or by using equality constraints ('equal.com=TRUE' and/or 'equal.ang=TRUE'); \n* THE STANDARD ERRORS OF THE MODEL PARAMETER ESTIMATES MAY NOT BE TRUSTWORTHY; \n* DEGREE OF FREEDOM = 0: SEVERAL FIT INDEXES CANNOT BE CALCULATED;...")
	ASK()}
 
 
append<-
function (x, values, after = length(x)) 
{
    lengx <- length(x)
    if (after <= 0) 
        c(values, x)
    else if (after >= lengx) 
        c(x, values)
    else c(x[1:after], values, x[(after + 1):lengx])
}


objective1max<- 
function(par){
        ang1=par[1:(p-1)]
        ang=append(ang1,0,r-1)  
      if(mcsc=="unconstrained"){if(m==1){b=  c(par[((p-1)+1)],1)} else b=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]); b=b/sum(b)}
      if(mcsc=="-1" ){if(m==1){b=  c(0,1)} else if(m>=3){b=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);b[c(seq(1,(m+1),by=2))]=0; b=b/sum(b)} else {b=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);b[c(seq(1,(m+1),by=2))]=0; b=b/sum(b)}}
      if(mcsc=="0" ){ if(m==1){b=  c(0.5,0.5)} else if(m>=3){b=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);b[1]=sum(b[seq(2,m+1,by=2)])-sum(b[seq(3,m+1,by=2)]); b=b/sum(b)} else {b=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);b[1]=sum(b[seq(2,m+1,by=2)])-b[3]; b=b/sum(b)}}


        
        v=  par[((p-1)+(m)+1)]
                                  z=  par[((p-1)+(m)+1+1):((p-1)+(m)+1+p)]
  
 K=pi/180
 M=matrix(c(0),p,p,byrow=TRUE)
  for(i in 1:p){
        for(j in 1:p){
        M[i,j]=c(b[-c(1)])%*%cos(c(1:m)*(ang[j]-ang[i]))
        }}
 Pc=M+matrix(b[1],p,p)
 Dv=diag(v,p)
 Dz=diag(z); f=-1*(sum(diag(R%*%solve(Dz%*%(Pc+Dv)%*%Dz)))+log(det(Dz%*%(Pc+Dv)%*%Dz))-log(det(R))-p);
 S1=Dz%*%(Pc+Dv)%*%Dz;S2=diag(1/sqrt(diag(S1)))%*%Dz;S=S2%*%(Pc+Dv)%*%S2
 attributes(f)=list(f=f,S=S,Cs=Dz%*%(Pc+Dv)%*%Dz,Pc=Pc)
 f}
 
objective1gr<-
function(par){
        ang1=par[1:(p-1)]
        ang=append(ang1,0,r-1)  
      if(mcsc=="unconstrained"){if(m==1){alpha=  c(par[((p-1)+1)],1)} else alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]); b=alpha/sum(alpha)}
      if(mcsc=="-1" ){if(m==1){b=alpha=  c(0,1)} else if(m>=3){alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)} else {alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)}}
      if(mcsc=="0" ){ if(m==1){b=alpha=  c(0.5,0.5)} else if(m>=3){alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]); b=alpha/sum(alpha)} else {alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-alpha[3]; b=alpha/sum(alpha)}}          
        v=  par[((p-1)+(m)+1)]
                                  z=  par[((p-1)+(m)+1+1):((p-1)+(m)+1+p)]
 K=pi/180
 M=matrix(c(0),p,p,byrow=TRUE)
  for(i in 1:p){
        for(j in 1:p){
        M[i,j]=c(b[-c(1)])%*%cos(c(1:m)*(ang[j]-ang[i]))
        }}
 Pc=M+matrix(b[1],p,p)
 Dv=diag(v,p)
 Dz=diag(z); 
 S=Dz%*%(Pc+Dv)%*%Dz
 
 hessian=matrix(0,length(par),length(par))

 gradientz=rep(0,p);Dpz=matrix(0,p*p,p)
for(i in 1:p){
        J=matrix(0,p,p)
        J[i,i]=1;
        dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
        gradientz[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpz[,i]=c(dp)
        }
gradientv=rep(0,p);Dpv=matrix(0,p*p,1)
        J=diag(1,p,p) 
        dp=J%*%(Dz)^2
        gradientv=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpv=c(dp)
       
 if((mcsc=="unconstrained")){       
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
for(i in 2:(m+1)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(    (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))  -1*(alpha[1]/sum(alpha)^2+(alpha[-c(1)]/sum(alpha))^2%*%cos(c(1:m)*(ang[j]-ang[z])))       )*Dz[z,z]*Dz[j,j] 
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
        
        
        for(i in 1:1){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)-alpha[i]/sum(alpha)^2 -1*sum( (alpha[-c(1)]/sum(alpha)^2)*cos(c(1:m)*(ang[j]-ang[z]) ) ) )*Dz[z,z]*Dz[j,j]
        }}
dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        
        Dpalph[,i]=c(dp)
        }  
}

 if((mcsc=="-1")){      
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(   (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))  -1*(alpha[1]/sum(alpha)^2+(alpha[-c(1)]/sum(alpha))^2%*%cos(c(1:m)*(ang[j]-ang[z])))      )*Dz[z,z]*Dz[j,j]  
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
           }
   
        
                 }
if((mcsc=="-1" & m<=2)){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
}
if((mcsc=="0")){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha) 
        
        - (sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]))*1/sum(alpha) 
        
        -(sum(alpha[-c(1,i)]*1/sum(alpha)*cos(c(seq(1,m,by=1)[-c(i-1)])*(ang[j]-ang[z]))))  
        
        +(1/sum(alpha)-alpha[i]*1/sum(alpha))*cos(c(i-1)*(ang[j]-ang[z]))   )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        } 
}
if(  m>=2 ){ 
for(i in seq(3,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)*cos(c(i-1)*(ang[j]-ang[z])) -(1/sum(alpha))  )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
}
  

if((mcsc=="0" & m==1)){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
} 
}

gradientang=rep(0,p);Dpang=matrix(0,p*p,length(ang))
for(i in 1:p){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
                
                if(z==i) s=1
                if(z!=i) s=0
                if(j==i) sj=1
                if(j!=i) sj=0
        M[z,j]= (s*c(1:m)%*%(b[-c(1)]*sin(c(1:m)*(ang[j]-ang[i]))) - sj*c(1:m)%*%(b[-c(1)]*sin(c(1:m)*(ang[i]-ang[z]))))*Dz[z,z]*Dz[j,j]
        }}
        dp=M;
        Dpang[,i]=c(dp)
        
        gradientang[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        
        
        };Dpang=Dpang[,-c(1)];
        gradient=c(gradientang[-c(r)], if((mcsc=="unconstrained") | (mcsc=="-1" & m>=3) | (mcsc=="0" & m>=2)){       
gradientalph[-c(2)]} else if((mcsc=="-1" & m==1) | (mcsc=="0" & m==1)){0} else if((mcsc=="-1" & m==2)){c(0,0)},gradientv,gradientz)
        

 
 gradient}

objective1hess<-
function(par){
        ang1=par[1:(p-1)]
        ang=append(ang1,0,r-1)  
      if(mcsc=="unconstrained"){if(m==1){alpha=  c(par[((p-1)+1)],1)} else alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]); b=alpha/sum(alpha)}
      if(mcsc=="-1" ){if(m==1){b=alpha=  c(0,1)} else if(m>=3){alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)} else {alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)}}
      if(mcsc=="0" ){if(m==1){b=alpha=c(0.5,0.5)} else if(m>=3){alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]); b=alpha/sum(alpha)} else {alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-alpha[3]; b=alpha/sum(alpha)}}

        v=  par[((p-1)+(m)+1)]
                                  z=  par[((p-1)+(m)+1+1):((p-1)+(m)+1+p)]
         K=pi/180
 M=matrix(c(0),p,p,byrow=TRUE)
  for(i in 1:p){
        for(j in 1:p){
        M[i,j]=c(b[-c(1)])%*%cos(c(1:m)*(ang[j]-ang[i]))
        }}
 Pc=M+matrix(b[1],p,p)
 Dv=diag(v,p)
 Dz=diag(z); 

 S=Dz%*%(Pc+Dv)%*%Dz;S2=diag(1/sqrt(diag(S)))%*%Dz;S1=S2%*%(Pc+Dv)%*%S2
 hessian=matrix(0,length(par),length(par))

 gradientz=rep(0,p);Dpz=matrix(0,p*p,p)
for(i in 1:p){
        J=matrix(0,p,p)
        J[i,i]=1;
        dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
        gradientz[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpz[,i]=c(dp)
        }
gradientv=rep(0,p);Dpv=matrix(0,p*p,1)
        J=diag(1,p,p) 
        dp=J%*%(Dz)^2
        gradientv=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpv=c(dp)


 if((mcsc=="unconstrained")){       
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
for(i in 2:(m+1)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(   (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))      )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
        
        
        for(i in 1:1){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)-alpha[i]/sum(alpha)^2 -1*sum( (alpha[-c(1)]/sum(alpha)^2)*cos(c(1:m)*(ang[j]-ang[z]) ) ) )*Dz[z,z]*Dz[j,j]
        }}
dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        
        Dpalph[,i]=c(dp)
        };if(mcsc=="unconstrained"){Dpalph=Dpalph[,-c(2)]} ;  
}
 if((mcsc=="-1")){      
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(   (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))  -1*(alpha[1]/sum(alpha)^2+(alpha[-c(1)]/sum(alpha))^2%*%cos(c(1:m)*(ang[j]-ang[z])))     )*Dz[z,z]*Dz[j,j]  
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
           Dpalph=Dpalph[,-c(2)]}
   
    
if(m<=2){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
Dpalph=Dpalph[,-c(2)]   }
                 
                  }


if((mcsc=="0")){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha) 
        
        - (sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]))*1/sum(alpha) 
        
        -(sum(alpha[-c(1,i)]*1/sum(alpha)*cos(c(seq(1,m,by=1)[-c(i-1)])*(ang[j]-ang[z]))))  
        
        +(1/sum(alpha)-alpha[i]*1/sum(alpha))*cos(c(i-1)*(ang[j]-ang[z]))   )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        } 
}
if(  m>=2 ){
for(i in seq(3,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)*cos(c(i-1)*(ang[j]-ang[z])) -(1/sum(alpha))  )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
}
  Dpalph=Dpalph[,-c(2)]

}
if( (mcsc=="0" & m==1)){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(b))
 Dpalph=Dpalph[,-c(2)]
}


gradientang=rep(0,p);Dpang=matrix(0,p*p,length(ang))
for(i in 1:p){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
                
                if(z==i) s=1
                if(z!=i) s=0
                if(j==i) sj=1
                if(j!=i) sj=0
        M[z,j]= (s*c(1:m)%*%(b[-c(1)]*sin(c(1:m)*(ang[j]-ang[i]))) - sj*c(1:m)%*%(b[-c(1)]*sin(c(1:m)*(ang[i]-ang[z]))))*Dz[z,z]*Dz[j,j]
        }}
        dp=M;
        Dpang[,i]=c(dp)
        
        gradientang[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        
        
        };Dpang=Dpang[,-c(1)];
        gradient=c(gradientang[-c(r)], if((mcsc=="unconstrained") | (mcsc=="-1" & m>=3) | (mcsc=="0" & m>=2)){       
gradientalph[-c(2)]}  else if((mcsc=="-1" & m==1) | (mcsc=="0" & m==1)){0} else if((mcsc=="-1" & m==2)){c(0,0)},gradientv,gradientz)
        if((mcsc=="unconstrained") | (mcsc=="-1" & m>=1) | (mcsc=="0" & m>=1)){
        DerivAnal=data.frame(Dpang,Dpalph,Dpv,Dpz);
        for(i in 1:length(par)){
        for(j in 1:length(par)){
     
                      hessian[i,j]=-1*sum(diag(  solve(S)%*%matrix(DerivAnal[,i],p,p)%*%solve(S)%*%matrix(DerivAnal[,j],p,p) ))

  }}} 
 
 
  hessian}




objective2max<-
function(par){
        ang1=par[1:(p-1)]
        ang=append(ang1,0,r-1)  
      if(mcsc=="unconstrained"){if(m==1){b=  c(par[((p-1)+1)],1)} else b=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);b=b/sum(b)}
      if(mcsc=="-1" ){if(m==1){b=  c(0,1)} else if(m>=3){b=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);b[c(seq(1,(m+1),by=2))]=0; b=b/sum(b)} else {b=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);b[c(seq(1,(m+1),by=2))]=0; b=b/sum(b)}}
      if(mcsc=="0" ){ if(m==1){b=  c(0.5,0.5)} else if(m>=3){b=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);b[1]=sum(b[seq(2,m+1,by=2)])-sum(b[seq(3,m+1,by=2)]); b=b/sum(b)} else {b=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);b[1]=sum(b[seq(2,m+1,by=2)])-b[3]; b=b/sum(b)}}

       
        v=  par[((p-1)+(m)+1):((p-1)+(m)+p)]
                                  z=  par[((p-1)+(m)+p+1):((p-1)+(m)+p+p)]
 K=pi/180
 M=matrix(c(0),p,p,byrow=TRUE)
  for(i in 1:p){
        for(j in 1:p){
        M[i,j]=c(b[-c(1)])%*%cos(c(1:m)*(ang[j]-ang[i]))
        }}
 Pc=M+matrix(b[1],p,p)
 Dv=diag(v,p)
 Dz=diag(z); f=-1*(sum(diag(R%*%solve(Dz%*%(Pc+Dv)%*%Dz)))+log(det(Dz%*%(Pc+Dv)%*%Dz))-log(det(R))-p);
 S1=Dz%*%(Pc+Dv)%*%Dz;S2=diag(1/sqrt(diag(S1)))%*%Dz;S=S2%*%(Pc+Dv)%*%S2
 attributes(f)=list(f=f,S=S,Cs=Dz%*%(Pc+Dv)%*%Dz,Pc=Pc)
 f}
 
objective2gr<-
function(par){
        ang1=par[1:(p-1)]
        ang=append(ang1,0,r-1)  
        if(mcsc=="unconstrained"){if(m==1){alpha=  c(par[((p-1)+1)],1)} else alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);b=alpha/sum(alpha)}
      if(mcsc=="-1" ){if(m==1){b=alpha=  c(0,1)} else if(m>=3){alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)} else {alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)}}
      if(mcsc=="0" ){ if(m==1){b=alpha= c(0.5,0.5)} else if(m>=3){alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]); b=alpha/sum(alpha)} else {alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-alpha[3]; b=alpha/sum(alpha)}}          

        v=  par[((p-1)+(m)+1):((p-1)+(m)+p)]
                                  z=  par[((p-1)+(m)+p+1):((p-1)+(m)+p+p)]
 K=pi/180
 M=matrix(c(0),p,p,byrow=TRUE)
  for(i in 1:p){
        for(j in 1:p){
        M[i,j]=c(b[-c(1)])%*%cos(c(1:m)*(ang[j]-ang[i]))
        }}
 Pc=M+matrix(b[1],p,p)
 Dv=diag(v,p)
 Dz=diag(z); 
 S=Dz%*%(Pc+Dv)%*%Dz
 
 hessian=matrix(0,length(par),length(par))

 gradientz=rep(0,p);Dpz=matrix(0,p*p,p)
for(i in 1:p){
        J=matrix(0,p,p)
        J[i,i]=1;
        dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
        gradientz[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpz[,i]=c(dp)
        }
gradientv=rep(0,p);Dpv=matrix(0,p*p,p)
for(i in 1:p){
        J=matrix(0,p,p)
        J[i,i]=1; 
        dp=J%*%(Dz)^2
        gradientv[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpv[,i]=c(dp)
        }

 if((mcsc=="unconstrained")){       
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
for(i in 2:(m+1)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(   (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))  -1*(alpha[1]/sum(alpha)^2+(alpha[-c(1)]/sum(alpha))^2%*%cos(c(1:m)*(ang[j]-ang[z])))     )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
        
        
        for(i in 1:1){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)-alpha[i]/sum(alpha)^2 -1*sum( (alpha[-c(1)]/sum(alpha)^2)*cos(c(1:m)*(ang[j]-ang[z]) ) ) )*Dz[z,z]*Dz[j,j]
        }}
dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        
        Dpalph[,i]=c(dp)
        }
}

 if((mcsc=="-1")){      
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(    (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))  -1*(alpha[1]/sum(alpha)^2+(alpha[-c(1)]/sum(alpha))^2%*%cos(c(1:m)*(ang[j]-ang[z])))    )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
           }
   
        
                 }
   
        
                 
if((mcsc=="-1" & m<=2)){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
}
if((mcsc=="0")){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha) 
        
        - (sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]))*1/sum(alpha) 
        
        -(sum(alpha[-c(1,i)]*1/sum(alpha)*cos(c(seq(1,m,by=1)[-c(i-1)])*(ang[j]-ang[z]))))  
        
        +(1/sum(alpha)-alpha[i]*1/sum(alpha))*cos(c(i-1)*(ang[j]-ang[z]))   )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        } 
}
if(  m>=2 ){ 
for(i in seq(3,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)*cos(c(i-1)*(ang[j]-ang[z])) -(1/sum(alpha))  )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
}
 
}
if((mcsc=="0" & m==1)){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
 
}
gradientang=rep(0,p);Dpang=matrix(0,p*p,length(ang))
for(i in 1:p){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
                
                if(z==i) s=1
                if(z!=i) s=0
                if(j==i) sj=1
                if(j!=i) sj=0
        M[z,j]= (s*c(1:m)%*%(b[-c(1)]*sin(c(1:m)*(ang[j]-ang[i]))) - sj*c(1:m)%*%(b[-c(1)]*sin(c(1:m)*(ang[i]-ang[z]))))*Dz[z,z]*Dz[j,j]
        }}
        dp=M;
        Dpang[,i]=c(dp)
        
        gradientang[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        
        
        };Dpang=Dpang[,-c(1)];
        gradient=c(gradientang[-c(r)], if((mcsc=="unconstrained") | (mcsc=="-1" & m>=3) | (mcsc=="0" & m>=2)){       
gradientalph[-c(2)]} else if((mcsc=="-1" & m==1) | (mcsc=="0" & m==1)){0} else if((mcsc=="-1" & m==2)){c(0,0)},gradientv,gradientz)
        
 
 
 gradient}

objective2hess<-
function(par){
        ang1=par[1:(p-1)]
        ang=append(ang1,0,r-1)  
        if(mcsc=="unconstrained"){if(m==1){alpha=  c(par[((p-1)+1)],1)} else alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);b=alpha/sum(alpha)}
     if(mcsc=="-1" ){if(m==1){b=alpha=  c(0,1)} else if(m>=3){alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)} else {alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)}}
      if(mcsc=="0" ){if(m==1){b=alpha=  c(0.5,0.5)} else if(m>=3){alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]); b=alpha/sum(alpha)} else {alpha=  c(par[((p-1)+1)],1,par[((p-1)+2):((p-1)+(m))]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-alpha[3]; b=alpha/sum(alpha)}}

        v=  par[((p-1)+(m)+1):((p-1)+(m)+p)]
                                  z=  par[((p-1)+(m)+p+1):((p-1)+(m)+p+p)]
 K=pi/180
 M=matrix(c(0),p,p,byrow=TRUE)
  for(i in 1:p){
        for(j in 1:p){
        M[i,j]=c(b[-c(1)])%*%cos(c(1:m)*(ang[j]-ang[i]))
        }}
 Pc=M+matrix(b[1],p,p)
 Dv=diag(v,p)
 Dz=diag(z); 

 S=Dz%*%(Pc+Dv)%*%Dz;S2=diag(1/sqrt(diag(S)))%*%Dz;S1=S2%*%(Pc+Dv)%*%S2
 hessian=matrix(0,length(par),length(par))

 gradientz=rep(0,p);Dpz=matrix(0,p*p,p)
for(i in 1:p){
        J=matrix(0,p,p)
        J[i,i]=1;
        dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
        gradientz[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpz[,i]=c(dp)
        }
gradientv=rep(0,p);Dpv=matrix(0,p*p,p)
for(i in 1:p){
        J=matrix(0,p,p)
        J[i,i]=1; 
        dp=J%*%(Dz)^2
        gradientv[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpv[,i]=c(dp)
        }


 if((mcsc=="unconstrained")){       
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
for(i in 2:(m+1)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(   (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))  -1*(alpha[1]/sum(alpha)^2+(alpha[-c(1)]/sum(alpha))^2%*%cos(c(1:m)*(ang[j]-ang[z])))     )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
        
        
        for(i in 1:1){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)-alpha[i]/sum(alpha)^2 -1*sum( (alpha[-c(1)]/sum(alpha)^2)*cos(c(1:m)*(ang[j]-ang[z]) ) ) )*Dz[z,z]*Dz[j,j]
        }}
dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        
        Dpalph[,i]=c(dp)
        };Dpalph=Dpalph[,-c(2)] 
}
 if((mcsc=="-1")){      
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(   (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))  -1*(alpha[1]/sum(alpha)^2+(alpha[-c(1)]/sum(alpha))^2%*%cos(c(1:m)*(ang[j]-ang[z])))     )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
           Dpalph=Dpalph[,-c(2)]}
   
    
if(m<=2){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
Dpalph=Dpalph[,-c(2)]   }
                 
                  }


if((mcsc=="0")){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha) 
        
        - (sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]))*1/sum(alpha) 
        
        -(sum(alpha[-c(1,i)]*1/sum(alpha)*cos(c(seq(1,m,by=1)[-c(i-1)])*(ang[j]-ang[z]))))  
        
        +(1/sum(alpha)-alpha[i]*1/sum(alpha))*cos(c(i-1)*(ang[j]-ang[z]))   )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        } 
}
if(  m>=2 ){
for(i in seq(3,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)*cos(c(i-1)*(ang[j]-ang[z])) -(1/sum(alpha))  )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
}
  Dpalph=Dpalph[,-c(2)]

}
if( (mcsc=="0" & m==1)){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(b))
 Dpalph=Dpalph[,-c(2)]
}


gradientang=rep(0,p);Dpang=matrix(0,p*p,length(ang))
for(i in 1:p){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
                
                if(z==i) s=1
                if(z!=i) s=0
                if(j==i) sj=1
                if(j!=i) sj=0
        M[z,j]= (s*c(1:m)%*%(b[-c(1)]*sin(c(1:m)*(ang[j]-ang[i]))) - sj*c(1:m)%*%(b[-c(1)]*sin(c(1:m)*(ang[i]-ang[z]))))*Dz[z,z]*Dz[j,j]
        }}
        dp=M;
        Dpang[,i]=c(dp)
        
        gradientang[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        
        
        };Dpang=Dpang[,-c(1)];
        gradient=c(gradientang[-c(r)], if((mcsc=="unconstrained") | (mcsc=="-1" & m>=3) | (mcsc=="0" & m>=2)){       
gradientalph[-c(2)]}  else if((mcsc=="-1" & m==1) | (mcsc=="0" & m==1)){0} else if((mcsc=="-1" & m==2)){c(0,0)},gradientv,gradientz)
        if((mcsc=="unconstrained") | (mcsc=="-1" & m>=1) | (mcsc=="0" & m>=1)){
        DerivAnal=data.frame(Dpang,Dpalph,Dpv,Dpz);
        for(i in 1:length(par)){
        for(j in 1:length(par)){
     
                      hessian[i,j]=-1*sum(diag(  solve(S)%*%matrix(DerivAnal[,i],p,p)%*%solve(S)%*%matrix(DerivAnal[,j],p,p) ))

  }}} 
 
 hessian}




objective3max<-
function(par){
    ang=c(start.valuesA(R,k)$polar.angles)*K
      if(mcsc=="unconstrained"){if(m==1){alpha=c(par[(1)],1)
} else alpha=  c(par[(1)],1,par[(2):((m))]);b=alpha/sum(alpha)}
     if(mcsc=="-1" ){if(m==1){b=  c(0,1)} else if(m>=3){alpha=  c(par[(1)],1,par[(2):(m)]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)} else {alpha=  c(par[(1)],1,par[(2):(m)]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)}}
      if(mcsc=="0" ){if(m==1){b=  c(0.5,0.5)} else if(m>=3){alpha=  c(par[(1)],1,par[(2):(m)]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]); b=alpha/sum(alpha)} else {alpha=  c(par[(1)],1,par[(2):(m)]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-alpha[3]; b=alpha/sum(alpha)}}

        


         v=  par[((m)+1)]
                                   z=  par[((m)+1+1):((m)+1+p)]

K=pi/180
M=matrix(c(0),p,p,byrow=TRUE)
  for(i in 1:p){
        for(j in 1:p){
        M[i,j]=c(b[-c(1)])%*%cos(c(1:m)*(ang[j]-ang[i]))
        }}
 Pc=M+matrix(b[1],p,p)
 Dv=diag(v,p)
 Dz=diag(z); f=-1*(sum(diag(R%*%solve(Dz%*%(Pc+Dv)%*%Dz)))+log(det(Dz%*%(Pc+Dv)%*%Dz))-log(det(R))-p) ;
 S1=Dz%*%(Pc+Dv)%*%Dz;S2=diag(1/sqrt(diag(S1)))%*%Dz;S=S2%*%(Pc+Dv)%*%S2
 attributes(f)=list(f=f,S=S,Cs=Dz%*%(Pc+Dv)%*%Dz,Pc=Pc)

 f}
 
objective3gr<-
function(par){
    ang=c(start.valuesA(R,k)$polar.angles)*K
      if(mcsc=="unconstrained"){if(m==1){alpha=c(par[(1)],1)
} else alpha=  c(par[(1)],1,par[(2):((m))]);b=alpha/sum(alpha)}
     if(mcsc=="-1" ){if(m==1){b=alpha=  c(0,1)} else if(m>=3){alpha=  c(par[(1)],1,par[(2):(m)]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)} else {alpha=  c(par[(1)],1,par[(2):(m)]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)}}
      if(mcsc=="0" ){if(m==1){b=alpha=  c(0.5,0.5)} else if(m>=3){alpha=  c(par[(1)],1,par[(2):(m)]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]); b=alpha/sum(alpha)} else {alpha=  c(par[(1)],1,par[(2):(m)]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-alpha[3]; b=alpha/sum(alpha)}}

         v=  par[((m)+1)]
                                   z=  par[((m)+1+1):((m)+1+p)]

K=pi/180

M=matrix(c(0),p,p,byrow=TRUE)
  for(i in 1:p){
        for(j in 1:p){
        M[i,j]=c(b[-c(1)])%*%cos(c(1:m)*(ang[j]-ang[i]))
        }}
 Pc=M+matrix(b[1],p,p)
 Dv=diag(v,p)
 Dz=diag(z); 
 S=Dz%*%(Pc+Dv)%*%Dz


 hessian=matrix(0,length(par),length(par))

 gradientz=rep(0,p);Dpz=matrix(0,p*p,p)
for(i in 1:p){
        J=matrix(0,p,p)
        J[i,i]=1;
        dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
        gradientz[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpz[,i]=c(dp)
        }
gradientv=rep(0,p);Dpv=matrix(0,p*p,1)
        J=diag(1,p,p) 
        dp=J%*%(Dz)^2
        gradientv=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpv=c(dp)

 if((mcsc=="unconstrained")){       
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
for(i in 2:(m+1)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(    (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))  -1*(alpha[1]/sum(alpha)^2+(alpha[-c(1)]/sum(alpha))^2%*%cos(c(1:m)*(ang[j]-ang[z])))       )*Dz[z,z]*Dz[j,j] 
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
        
        
        for(i in 1:1){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)-alpha[i]/sum(alpha)^2 -1*sum( (alpha[-c(1)]/sum(alpha)^2)*cos(c(1:m)*(ang[j]-ang[z]) ) ) )*Dz[z,z]*Dz[j,j]
        }}
dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        
        Dpalph[,i]=c(dp)
        }  
}

 if((mcsc=="-1")){      
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(   (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))  -1*(alpha[1]/sum(alpha)^2+(alpha[-c(1)]/sum(alpha))^2%*%cos(c(1:m)*(ang[j]-ang[z])))      )*Dz[z,z]*Dz[j,j]  
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
           }
   
if(( m<=2)){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
}
        
                 }


if((mcsc=="0")){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha) 
        
        - (sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]))*1/sum(alpha) 
        
        -(sum(alpha[-c(1,i)]*1/sum(alpha)*cos(c(seq(1,m,by=1)[-c(i-1)])*(ang[j]-ang[z]))))  
        
        +(1/sum(alpha)-alpha[i]*1/sum(alpha))*cos(c(i-1)*(ang[j]-ang[z]))   )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        } 
}
if(  m>=2 ){ 
for(i in seq(3,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)*cos(c(i-1)*(ang[j]-ang[z])) -(1/sum(alpha))  )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
}
  
}
if((mcsc=="0" & m==1)){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
}

gradient=c( if((mcsc=="unconstrained") | (mcsc=="-1" & m>=3) | (mcsc=="0" & m>=2)){       
gradientalph[-c(2)]} else if((mcsc=="-1" & m==1) | (mcsc=="0" & m==1)){0} else if((mcsc=="-1" & m==2)){c(0,0)},gradientv,gradientz)        

 
 gradient}

objective3hess<-
function(par){
    ang=c(start.valuesA(R,k)$polar.angles)*K
      if(mcsc=="unconstrained"){if(m==1){alpha=c(par[(1)],1)
} else alpha=  c(par[(1)],1,par[(2):((m))]);b=alpha/sum(alpha)}
     if(mcsc=="-1" ){if(m==1){b=alpha=  c(0,1)} else if(m>=3){alpha=  c(par[(1)],1,par[(2):(m)]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)} else {alpha=  c(par[(1)],1,par[(2):(m)]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)}}
      if(mcsc=="0" ){if(m==1){b=alpha=  c(0.5,0.5)} else if(m>=3){alpha=  c(par[(1)],1,par[(2):(m)]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]); b=alpha/sum(alpha)} else {alpha=  c(par[(1)],1,par[(2):(m)]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-alpha[3]; b=alpha/sum(alpha)}}

         v=  par[((m)+1)]
                                   z=  par[((m)+1+1):((m)+1+p)]

K=pi/180
M=matrix(c(0),p,p,byrow=TRUE)
  for(i in 1:p){
        for(j in 1:p){
        M[i,j]=c(b[-c(1)])%*%cos(c(1:m)*(ang[j]-ang[i]))
        }}
 Pc=M+matrix(b[1],p,p)
 Dv=diag(v,p)
 Dz=diag(z); 

 S=Dz%*%(Pc+Dv)%*%Dz;S2=diag(1/sqrt(diag(S)))%*%Dz;S1=S2%*%(Pc+Dv)%*%S2
 

 hessian=matrix(0,length(par),length(par))

 gradientz=rep(0,p);Dpz=matrix(0,p*p,p)
for(i in 1:p){
        J=matrix(0,p,p)
        J[i,i]=1;
        dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
        gradientz[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpz[,i]=c(dp)
        }
gradientv=rep(0,p);Dpv=matrix(0,p*p,1)
        J=diag(1,p,p) 
        dp=J%*%(Dz)^2
        gradientv=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpv=c(dp)

 if((mcsc=="unconstrained")){       
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
for(i in 2:(m+1)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(    (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))  -1*(alpha[1]/sum(alpha)^2+(alpha[-c(1)]/sum(alpha))^2%*%cos(c(1:m)*(ang[j]-ang[z])))       )*Dz[z,z]*Dz[j,j] 
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
        
        
        for(i in 1:1){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)-alpha[i]/sum(alpha)^2 -1*sum( (alpha[-c(1)]/sum(alpha)^2)*cos(c(1:m)*(ang[j]-ang[z]) ) ) )*Dz[z,z]*Dz[j,j]
        }}
dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        
        Dpalph[,i]=c(dp)
        };if(mcsc=="unconstrained"){Dpalph=Dpalph[,-c(2)]} ;  
}

 if((mcsc=="-1")){      
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(   (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))  -1*(alpha[1]/sum(alpha)^2+(alpha[-c(1)]/sum(alpha))^2%*%cos(c(1:m)*(ang[j]-ang[z])))      )*Dz[z,z]*Dz[j,j]  
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
           };Dpalph=Dpalph[,-c(2)]
   
        
                 }
if((mcsc=="-1" & m<=2)){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
Dpalph=Dpalph[,-c(2)]}

if((mcsc=="0")){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha) 
        
        - (sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]))*1/sum(alpha) 
        
        -(sum(alpha[-c(1,i)]*1/sum(alpha)*cos(c(seq(1,m,by=1)[-c(i-1)])*(ang[j]-ang[z]))))  
        
        +(1/sum(alpha)-alpha[i]*1/sum(alpha))*cos(c(i-1)*(ang[j]-ang[z]))   )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        } 
}
if(  m>=2 ){ 
for(i in seq(3,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)*cos(c(i-1)*(ang[j]-ang[z])) -(1/sum(alpha))  )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
}
  Dpalph=Dpalph[,-c(2)]
}

if((mcsc=="0" & m==1)){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
 Dpalph=Dpalph[,-c(2)]
}

        gradient=c( if((mcsc=="unconstrained") | (mcsc=="-1" & m>=3) | (mcsc=="0" & m>=2)){       
gradientalph[-c(2)]}  else if((mcsc=="-1" & m==1) | (mcsc=="0" & m==1)){0} else if((mcsc=="-1" & m==2)){c(0,0)},gradientv,gradientz)
        if((mcsc=="unconstrained") | (mcsc=="-1" & m>=1) | (mcsc=="0" & m>=1)){
        DerivAnal=data.frame(Dpalph,Dpv,Dpz);
        for(i in 1:length(par)){
        for(j in 1:length(par)){
     
                      hessian[i,j]=-1*sum(diag(  solve(S)%*%matrix(DerivAnal[,i],p,p)%*%solve(S)%*%matrix(DerivAnal[,j],p,p) ))
 
  }}} 
 
 hessian}
 



objective4max=
function(par){
    ang=c(start.valuesA(R,k)$polar.angles)*K
      if(mcsc=="unconstrained"){if(m==1){alpha=  c(par[(1)],1)} else alpha=  c(par[(1)],1,par[(2):((m))]);b=alpha/sum(alpha)}
      if(mcsc=="-1" ){if(m==1){b=  c(0,1)} else if(m>=3){alpha=  c(par[(1)],1,par[(2):(m)]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)} else {alpha=  c(par[(1)],1,par[(2):(m)]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)}}
      if(mcsc=="0" ){if(m==1){b=c(0.5,0.5)} else if(m>=3){alpha=  c(par[(1)],1,par[(2):(m)]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]); b=alpha/sum(alpha)} else {alpha=  c(par[(1)],1,par[(2):(m)]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-alpha[3]; b=alpha/sum(alpha)}}

         v=  par[((m)+1):((m)+p)]
                                   z=  par[((m)+p+1):((m)+p+p)]

 K=pi/180
 M=matrix(c(0),p,p,byrow=TRUE)
  for(i in 1:p){
        for(j in 1:p){
        M[i,j]=c(b[-c(1)])%*%cos(c(1:m)*(ang[j]-ang[i]))
        }}
 Pc=M+matrix(b[1],p,p)
 Dv=diag(v,p)
 Dz=diag(z); f=-1*(sum(diag(R%*%solve(Dz%*%(Pc+Dv)%*%Dz)))+log(det(Dz%*%(Pc+Dv)%*%Dz))-log(det(R))-p) ;
 S1=Dz%*%(Pc+Dv)%*%Dz;S2=diag(1/sqrt(diag(S1)))%*%Dz;S=S2%*%(Pc+Dv)%*%S2
 attributes(f)=list(f=f,S=S,Cs=Dz%*%(Pc+Dv)%*%Dz,Pc=Pc)
 f}
 
 objective4gr<-
function(par){
          ang=c(start.valuesA(R,k)$polar.angles)*K
      if(mcsc=="unconstrained"){if(m==1){alpha=  c(par[(1)],1)} else alpha=  c(par[(1)],1,par[(2):((m))]);b=alpha/sum(alpha)}
     if(mcsc=="-1" ){if(m==1){b=alpha=  c(0,1)} else if(m>=3){alpha=  c(par[(1)],1,par[(2):(m)]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)} else {alpha=  c(par[(1)],1,par[(2):(m)]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)}}
      if(mcsc=="0" ){if(m==1){b=alpha=  c(0.5,0.5)} else if(m>=3){alpha=  c(par[(1)],1,par[(2):(m)]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]); b=alpha/sum(alpha)} else {alpha=  c(par[(1)],1,par[(2):(m)]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-alpha[3]; b=alpha/sum(alpha)}}
         v=  par[((m)+1):((m)+p)]
                                   z=  par[((m)+p+1):((m)+p+p)]

 K=pi/180
 M=matrix(c(0),p,p,byrow=TRUE)
  for(i in 1:p){
        for(j in 1:p){
        M[i,j]=c(b[-c(1)])%*%cos(c(1:m)*(ang[j]-ang[i]))
        }}
 Pc=M+matrix(b[1],p,p)
 Dv=diag(v,p)
 Dz=diag(z); 
 S=Dz%*%(Pc+Dv)%*%Dz
 

 hessian=matrix(0,length(par),length(par))

 gradientz=rep(0,p);Dpz=matrix(0,p*p,p)
for(i in 1:p){
        J=matrix(0,p,p)
        J[i,i]=1;
        dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
        gradientz[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpz[,i]=c(dp)
        }
gradientv=rep(0,p);Dpv=matrix(0,p*p,p)
for(i in 1:p){
        J=matrix(0,p,p)
        J[i,i]=1; 
        dp=J%*%(Dz)^2
        gradientv[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpv[,i]=c(dp)
        }

 if((mcsc=="unconstrained")){       
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
for(i in 2:(m+1)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(   (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))  -1*(alpha[1]/sum(alpha)^2+(alpha[-c(1)]/sum(alpha))^2%*%cos(c(1:m)*(ang[j]-ang[z])))     )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
        
        
        for(i in 1:1){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)-alpha[i]/sum(alpha)^2 -1*sum( (alpha[-c(1)]/sum(alpha)^2)*cos(c(1:m)*(ang[j]-ang[z]) ) ) )*Dz[z,z]*Dz[j,j]
        }}
dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        
        Dpalph[,i]=c(dp)
        }  
}

 if((mcsc=="-1")){      
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(    (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))  -1*(alpha[1]/sum(alpha)^2+(alpha[-c(1)]/sum(alpha))^2%*%cos(c(1:m)*(ang[j]-ang[z])))    )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
           }
   
        
                 }
   
        
                 
if((mcsc=="-1" & m<=2)){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
}
if((mcsc=="0")){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha) 
        
        - (sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]))*1/sum(alpha) 
        
        -(sum(alpha[-c(1,i)]*1/sum(alpha)*cos(c(seq(1,m,by=1)[-c(i-1)])*(ang[j]-ang[z]))))  
        
        +(1/sum(alpha)-alpha[i]*1/sum(alpha))*cos(c(i-1)*(ang[j]-ang[z]))   )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        } 
}
if(  m>=2 ){ 
for(i in seq(3,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)*cos(c(i-1)*(ang[j]-ang[z])) -(1/sum(alpha))  )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
}
  
}
if((mcsc=="0" & m==1)){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
 
}
gradient=c( if((mcsc=="unconstrained") | (mcsc=="-1" & m>=3) | (mcsc=="0" & m>=2)){       
gradientalph[-c(2)]} else if((mcsc=="-1" & m==1) | (mcsc=="0" & m==1)){0} else if((mcsc=="-1" & m==2)){c(0,0)},gradientv,gradientz)        

 
 gradient}
 
objective4hess<-
function(par){
          ang=c(start.valuesA(R,k)$polar.angles)*K
      if(mcsc=="unconstrained"){if(m==1){alpha=  c(par[(1)],1)} else alpha=  c(par[(1)],1,par[(2):((m))]);b=alpha/sum(alpha)}
       if(mcsc=="-1" ){if(m==1){b=alpha=  c(0,1)} else if(m>=3){alpha=  c(par[(1)],1,par[(2):(m)]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)} else {alpha=  c(par[(1)],1,par[(2):(m)]);alpha[c(seq(1,(m+1),by=2))]=0; b=alpha/sum(alpha)}}
      if(mcsc=="0" ){if(m==1){b=alpha=  c(0.5,0.5)} else if(m>=3){alpha=  c(par[(1)],1,par[(2):(m)]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]); b=alpha/sum(alpha)} else {alpha=  c(par[(1)],1,par[(2):(m)]);alpha[1]=sum(alpha[seq(2,m+1,by=2)])-alpha[3]; b=alpha/sum(alpha)}}
         v=  par[((m)+1):((m)+p)]
                                   z=  par[((m)+p+1):((m)+p+p)]


 K=pi/180
 M=matrix(c(0),p,p,byrow=TRUE)
  for(i in 1:p){
        for(j in 1:p){
        M[i,j]=c(b[-c(1)])%*%cos(c(1:m)*(ang[j]-ang[i]))
        }}
 Pc=M+matrix(b[1],p,p)
 Dv=diag(v,p)
 Dz=diag(z); 
 
 S=Dz%*%(Pc+Dv)%*%Dz;S2=diag(1/sqrt(diag(S)))%*%Dz;S1=S2%*%(Pc+Dv)%*%S2
 

 hessian=matrix(0,length(par),length(par))

 gradientz=rep(0,p);Dpz=matrix(0,p*p,p)
for(i in 1:p){
        J=matrix(0,p,p)
        J[i,i]=1;
        dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
        gradientz[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpz[,i]=c(dp)
        }
gradientv=rep(0,p);Dpv=matrix(0,p*p,p)
for(i in 1:p){
        J=matrix(0,p,p)
        J[i,i]=1; 
        dp=J%*%(Dz)^2
        gradientv[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpv[,i]=c(dp)
        }

 if((mcsc=="unconstrained")){       
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
for(i in 2:(m+1)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(   (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))  -1*(alpha[1]/sum(alpha)^2+(alpha[-c(1)]/sum(alpha))^2%*%cos(c(1:m)*(ang[j]-ang[z])))     )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
        
        
        for(i in 1:1){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)-alpha[i]/sum(alpha)^2 -1*sum( (alpha[-c(1)]/sum(alpha)^2)*cos(c(1:m)*(ang[j]-ang[z]) ) ) )*Dz[z,z]*Dz[j,j]
        }}
dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        
        Dpalph[,i]=c(dp)
        };if(mcsc=="unconstrained"){Dpalph=Dpalph[,-c(2)]} ;  
}

 if((mcsc=="-1")){      
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(    (1/sum(alpha)-alpha[i]/sum(alpha)^2)*cos((i-1)*(ang[j]-ang[z]))  -1*(alpha[1]/sum(alpha)^2+(alpha[-c(1)]/sum(alpha))^2%*%cos(c(1:m)*(ang[j]-ang[z])))    )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
           };Dpalph=Dpalph[,-c(2)]
   
        
                 }
   
        
                 
if((mcsc=="-1" & m<=2)){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
Dpalph=Dpalph[,-c(2)]}

if((mcsc=="0")){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
if(  m>=3 ){ 
for(i in seq(4,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha) 
        
        - (sum(alpha[seq(2,m+1,by=2)])-sum(alpha[seq(3,m+1,by=2)]))*1/sum(alpha) 
        
        -(sum(alpha[-c(1,i)]*1/sum(alpha)*cos(c(seq(1,m,by=1)[-c(i-1)])*(ang[j]-ang[z]))))  
        
        +(1/sum(alpha)-alpha[i]*1/sum(alpha))*cos(c(i-1)*(ang[j]-ang[z]))   )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        } 
}
if(  m>=2 ){ 
for(i in seq(3,m+1,by=2)){
        M=matrix(c(0),p,p,byrow=TRUE)
  for(z in 1:p){
        for(j in 1:p){
        M[z,j]=(  1/sum(alpha)*cos(c(i-1)*(ang[j]-ang[z])) -(1/sum(alpha))  )*Dz[z,z]*Dz[j,j]
        }}
        dp=M
        gradientalph[i]=1*sum(diag(  solve(S)%*%(R-S)%*%solve(S)%*%dp  ))
        Dpalph[,i]=c(dp)
        }
}
if((mcsc=="0" & m==1)){ 
gradientalph=rep(0,(m+1));Dpalph=matrix(0,p*p,length(alpha))
}

  Dpalph=Dpalph[,-c(2)]
}

gradient=c( if((mcsc=="unconstrained") | (mcsc=="-1" & m>=3) | (mcsc=="0" & m>=2)){       
gradientalph[-c(2)]} else if((mcsc=="-1" & m==1) | (mcsc=="0" & m==1)){0} else if((mcsc=="-1" & m==2)){c(0,0)},gradientv,gradientz)        
        if((mcsc=="unconstrained") | (mcsc=="-1" & m>=1) | (mcsc=="0" & m>=1)){
        DerivAnal=data.frame(Dpalph,Dpv,Dpz);
        for(i in 1:length(par)){
        for(j in 1:length(par)){
     
                      hessian[i,j]=-1*sum(diag(  solve(S)%*%matrix(DerivAnal[,i],p,p)%*%solve(S)%*%matrix(DerivAnal[,j],p,p) ))

  }}}
 
 hessian}
 
 
 
 numericGradient <- function(f, t0, eps=1e-6, ...) {
   NPar <- length(t0)
   NVal <- length(f0 <- f(t0, ...))
   grad <- matrix(NA, NVal, NPar)
   row.names(grad) <- names(f0)
   colnames(grad) <- names(t0)
   for(i in 1:NPar) {
      t2 <- t1 <- t0
      t1[i] <- t0[i] - eps/2
      t2[i] <- t0[i] + eps/2
      grad[,i] <- (f(t2, ...) - f(t1, ...))/eps
   }
   return(grad)
}
 



maxL.BFGS.B <- function(fn, grad=NULL,
                    start,up,low,...) {
                    factr
                    pgtol
                    print.level
                    iterlim
                    
   message <- function(c) {
      switch(as.character(c),
               "0" = "successful convergence",
               "10" = "degeneracy in Nelder-Mead simplex"
               )
   };
if(equal.ang==FALSE ){
	if(m==1){par.names=c(v.names[-c(1)],paste(rep("a",m),c(0)),if(equal.com==TRUE) rep("v",1) else  paste(rep("v",p),v.names),paste(rep("z",p),v.names))}
 else par.names=c(v.names[-c(1)],paste(rep("a",m),c(0,2:m)),if(equal.com==TRUE) rep("v",1) else  paste(rep("v",p),v.names),paste(rep("z",p),v.names))}
   
if(equal.ang==TRUE ){
	if(m==1){par.names=c(paste(rep("a",m),c(0)),if(equal.com==TRUE) rep("v",1) else  paste(rep("v",p),v.names),paste(rep("z",p),v.names))} else par.names=c(paste(rep("a",m),c(0,2:m)),if(equal.com==TRUE) rep("v",1) else  paste(rep("v",p),v.names),paste(rep("z",p),v.names))}



   nParam <- length(start)
   func <- function(theta, ...) {
      sum(fn(theta, ...))
   }
   gradient <- function(theta, ...) {
      if(!is.null(grad)) {
         g <- grad(theta, ...)
         if(!is.null(dim(g))) {
            if(ncol(g) > 1) {
               return(colSums(g))
            }
         } else {
            return(g)
         }
      }
      g <- numericGradient(fn, theta, ...)
      if(!is.null(dim(g))) {
         return(colSums(g))
      } else {
         return(g)
      }
   }
   G1 <- gradient(start, ...)
   if(any(is.na(G1))) {
      stop("Na in the initial gradient")
   }
   if(length(G1) != nParam) {
      stop( "length of gradient (", length(G1),
         ") not equal to the no. of parameters (", nParam, ")" )
   }
   
   
   if( print.level == 1 ) {
   	  cat( "    -------------------------------\n")
      cat( "          Initial parameters:      \n")
      cat( "    -------------------------------\n")
      a <- cbind(start, G1, up,low)
      dimnames(a) <- list(par.names, c("parameter", "initial gradient",
                                          "upper","lower"))
      print(a)
      cat( "                               \n")
      cat( "                               \n")
      cat( "   Constrained (L-BFGS-B) Optimization\n")

         }

   type <- "L-BFGS-B maximisation"
   control <- list(trace=print.level,
                    REPORT=1,
                    fnscale=-1,
                    abstol=1e6,
                    maxit=iterlim,
                    factr=factr,
                    pgtol=pgtol,
                    lmm=if(is.null(lmm)){length(parA)} else lmm)
   a <- optim(start, func, gr=gradient, control=control, method="L-BFGS-B",
      hessian=FALSE,upper=up,lower=low, ...)


Active.Bounds.Up<-which(a$par==up)
Active.Bounds.Low<-which(a$par==low)
{if(length(Active.Bounds.Up)!=0 | length(Active.Bounds.Low)!=0){
Active.Bounds<-c(if(length(Active.Bounds.Up)!=0)Active.Bounds.Up,if(length(Active.Bounds.Low)!=0)Active.Bounds.Low)
Active.Bounds.names<-par.names[Active.Bounds]}
else Active.Bounds.names<-NULL}



   result <- list(
                   maximum=a$value,
                   estimate=a$par,
                   gradient=gradient(a$par),
                   hessian=a$hessian,
                   code=a$convergence,
                   message=paste(message(a$convergence), a$message),
                   last.step=NULL,
                   iterations=a$counts,
                   type=type,
                   Active.Bounds.names=Active.Bounds.names)
   class(result) <- "maximisation"
   invisible(result)
}


 if(is.vector(upper)) {up=c(upper[-c(1)]*K,rep(Inf,length(parA[-c(1:p-1)])))}
   if(is.null(upper)) {up=rep(Inf,length(parA))}
 if(is.vector(lower)){low=c(lower[-c(1)]*K,rep(0,length(parA[-c(1:p-1)])))}
   if(is.null(lower) & equal.ang==FALSE) {low=c(rep(-Inf,length(parA[c(1:p-1)])),rep(0,length(parA[-c(1:p-1)])))}
   if(is.null(lower) & equal.ang==TRUE) {low=c(rep(0,length(parA)))}

timing=system.time(res<-try(maxL.BFGS.B(if(equal.com==TRUE & equal.ang==FALSE)objective1max  else if(equal.com==FALSE & equal.ang==FALSE)objective2max else if(equal.com==TRUE & equal.ang==TRUE) objective3max else objective4max,grad=if(equal.com==TRUE & equal.ang==FALSE)objective1gr  else if(equal.com==FALSE & equal.ang==FALSE)objective2gr else if(equal.com==TRUE & equal.ang==TRUE) objective3gr else objective4gr,start=parA,up,low)))

if(is.character(res)==TRUE & try.refit.BFGS==TRUE){

cat("","\n");
cat("Fail to converge. Re-estimate removing default box constraints on z,v and a parameters... \n");
cat("","\n");

 if(is.vector(upper)) {up=c(upper[-c(1)]*K,rep(Inf,length(parA[-c(1:p-1)])))}
   if(is.null(upper)) {up=rep(Inf,length(parA))}
 if(is.vector(lower)){low=c(lower[-c(1)]*K,rep(-Inf,length(parA[-c(1:p-1)])))}
   if(is.null(lower) & equal.ang==FALSE) {low=c(rep(-Inf,length(parA)))}
   if(is.null(lower) & equal.ang==TRUE) {low=c(rep(-Inf,length(parA)))}



timing=system.time(res<-maxL.BFGS.B(if(equal.com==TRUE & equal.ang==FALSE)objective1max  else if(equal.com==FALSE & equal.ang==FALSE)objective2max else if(equal.com==TRUE & equal.ang==TRUE) objective3max else objective4max,grad=if(equal.com==TRUE & equal.ang==FALSE)objective1gr  else if(equal.com==FALSE & equal.ang==FALSE)objective2gr else if(equal.com==TRUE & equal.ang==TRUE) objective3gr else objective4gr,start=parA,up,low))
}

if(is.character(res)==TRUE){
cat("","\n");
cat("CONVERGENCE FAILURE: REDUCE m AND/OR lmm VALUES (default: lmm =",length(parA),"; suggested 3=<lmm<=20)... \n");
cat("","\n");
}

if(res$maximum>0){cat("Fail to converge, discrepancy function value is negative. Try to reduce m \n"); break}
if(print.level == 1) {
	  cat("\n")
      cat("Final gradient value:\n")
      print(res$gradient)
   }
param=res$est 
if(equal.ang==FALSE){res$est[1:(p-1)]=res$est[1:(p-1)]/K

for(i in 1:(p-1)){
if((res$est[i])<0)
res$est[i]=res$est[i]+360