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