# R/CircE.BFGS.R In CircE: Circumplex models Estimation

#### 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) {
}
else {
}
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
}
}
warning("the matrix has produced imaginary results -- proceed with caution")
}
if (mm > 1) {
sign.max <- vector(mode = "numeric", length = mm)
for (i in 1:mm) {
}
}
else {
if (abs(mini) > maxi) {
}
}

}

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

function (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;...")
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;...")

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))

for(i in 1:p){
J=matrix(0,p,p)
J[i,i]=1;
dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
Dpz[,i]=c(dp)
}
J=diag(1,p,p)
dp=J%*%(Dz)^2
Dpv=c(dp)

if((mcsc=="unconstrained")){
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
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

Dpalph[,i]=c(dp)
}
}

if((mcsc=="-1")){
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
Dpalph[,i]=c(dp)
}
}

}
if((mcsc=="-1" & m<=2)){
}
if((mcsc=="0")){
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
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
Dpalph[,i]=c(dp)
}
}

if((mcsc=="0" & m==1)){
}
}

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)

};Dpang=Dpang[,-c(1)];

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))

for(i in 1:p){
J=matrix(0,p,p)
J[i,i]=1;
dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
Dpz[,i]=c(dp)
}
J=diag(1,p,p)
dp=J%*%(Dz)^2
Dpv=c(dp)

if((mcsc=="unconstrained")){
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
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

Dpalph[,i]=c(dp)
};if(mcsc=="unconstrained"){Dpalph=Dpalph[,-c(2)]} ;
}
if((mcsc=="-1")){
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
Dpalph[,i]=c(dp)
}
Dpalph=Dpalph[,-c(2)]}

if(m<=2){
Dpalph=Dpalph[,-c(2)]   }

}

if((mcsc=="0")){
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
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
Dpalph[,i]=c(dp)
}
}
Dpalph=Dpalph[,-c(2)]

}
if( (mcsc=="0" & m==1)){
Dpalph=Dpalph[,-c(2)]
}

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)

};Dpang=Dpang[,-c(1)];
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))

for(i in 1:p){
J=matrix(0,p,p)
J[i,i]=1;
dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
Dpz[,i]=c(dp)
}
for(i in 1:p){
J=matrix(0,p,p)
J[i,i]=1;
dp=J%*%(Dz)^2
Dpv[,i]=c(dp)
}

if((mcsc=="unconstrained")){
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
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

Dpalph[,i]=c(dp)
}
}

if((mcsc=="-1")){
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
Dpalph[,i]=c(dp)
}
}

}

if((mcsc=="-1" & m<=2)){
}
if((mcsc=="0")){
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
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
Dpalph[,i]=c(dp)
}
}

}
if((mcsc=="0" & m==1)){

}
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)

};Dpang=Dpang[,-c(1)];

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))

for(i in 1:p){
J=matrix(0,p,p)
J[i,i]=1;
dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
Dpz[,i]=c(dp)
}
for(i in 1:p){
J=matrix(0,p,p)
J[i,i]=1;
dp=J%*%(Dz)^2
Dpv[,i]=c(dp)
}

if((mcsc=="unconstrained")){
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
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

Dpalph[,i]=c(dp)
};Dpalph=Dpalph[,-c(2)]
}
if((mcsc=="-1")){
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
Dpalph[,i]=c(dp)
}
Dpalph=Dpalph[,-c(2)]}

if(m<=2){
Dpalph=Dpalph[,-c(2)]   }

}

if((mcsc=="0")){
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
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
Dpalph[,i]=c(dp)
}
}
Dpalph=Dpalph[,-c(2)]

}
if( (mcsc=="0" & m==1)){
Dpalph=Dpalph[,-c(2)]
}

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)

};Dpang=Dpang[,-c(1)];
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))

for(i in 1:p){
J=matrix(0,p,p)
J[i,i]=1;
dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
Dpz[,i]=c(dp)
}
J=diag(1,p,p)
dp=J%*%(Dz)^2
Dpv=c(dp)

if((mcsc=="unconstrained")){
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
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

Dpalph[,i]=c(dp)
}
}

if((mcsc=="-1")){
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
Dpalph[,i]=c(dp)
}
}

if(( m<=2)){
}

}

if((mcsc=="0")){
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
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
Dpalph[,i]=c(dp)
}
}

}
if((mcsc=="0" & m==1)){
}

gradient=c( if((mcsc=="unconstrained") | (mcsc=="-1" & m>=3) | (mcsc=="0" & m>=2)){

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))

for(i in 1:p){
J=matrix(0,p,p)
J[i,i]=1;
dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
Dpz[,i]=c(dp)
}
J=diag(1,p,p)
dp=J%*%(Dz)^2
Dpv=c(dp)

if((mcsc=="unconstrained")){
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
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

Dpalph[,i]=c(dp)
};if(mcsc=="unconstrained"){Dpalph=Dpalph[,-c(2)]} ;
}

if((mcsc=="-1")){
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
Dpalph[,i]=c(dp)
}
};Dpalph=Dpalph[,-c(2)]

}
if((mcsc=="-1" & m<=2)){
Dpalph=Dpalph[,-c(2)]}

if((mcsc=="0")){
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
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
Dpalph[,i]=c(dp)
}
}
Dpalph=Dpalph[,-c(2)]
}

if((mcsc=="0" & m==1)){
Dpalph=Dpalph[,-c(2)]
}

gradient=c( if((mcsc=="unconstrained") | (mcsc=="-1" & m>=3) | (mcsc=="0" & m>=2)){
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))

for(i in 1:p){
J=matrix(0,p,p)
J[i,i]=1;
dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
Dpz[,i]=c(dp)
}
for(i in 1:p){
J=matrix(0,p,p)
J[i,i]=1;
dp=J%*%(Dz)^2
Dpv[,i]=c(dp)
}

if((mcsc=="unconstrained")){
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
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

Dpalph[,i]=c(dp)
}
}

if((mcsc=="-1")){
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
Dpalph[,i]=c(dp)
}
}

}

if((mcsc=="-1" & m<=2)){
}
if((mcsc=="0")){
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
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
Dpalph[,i]=c(dp)
}
}

}
if((mcsc=="0" & m==1)){

}
gradient=c( if((mcsc=="unconstrained") | (mcsc=="-1" & m>=3) | (mcsc=="0" & m>=2)){

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))

for(i in 1:p){
J=matrix(0,p,p)
J[i,i]=1;
dp=J%*%(Pc+Dv)%*%Dz+Dz%*%(Pc+Dv)%*%J
Dpz[,i]=c(dp)
}
for(i in 1:p){
J=matrix(0,p,p)
J[i,i]=1;
dp=J%*%(Dz)^2
Dpv[,i]=c(dp)
}

if((mcsc=="unconstrained")){
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
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

Dpalph[,i]=c(dp)
};if(mcsc=="unconstrained"){Dpalph=Dpalph[,-c(2)]} ;
}

if((mcsc=="-1")){
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
Dpalph[,i]=c(dp)
}
};Dpalph=Dpalph[,-c(2)]

}

if((mcsc=="-1" & m<=2)){
Dpalph=Dpalph[,-c(2)]}

if((mcsc=="0")){
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
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
Dpalph[,i]=c(dp)
}
}
if((mcsc=="0" & m==1)){
}

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

gradient=c( if((mcsc=="unconstrained") | (mcsc=="-1" & m>=3) | (mcsc=="0" & m>=2)){
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, ...))
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
}
}

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, ...))
}
if(!is.null(dim(g))) {
if(ncol(g) > 1) {
return(colSums(g))
}
} else {
return(g)
}
}
if(!is.null(dim(g))) {
return(colSums(g))
} else {
return(g)
}
}
if(any(is.na(G1))) {
}
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,
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")