Nothing
R/copstressMin.R
In cops: Cluster Optimized Proximity Scaling
#' Fitting a COPS-C Model (COPS Variant 1).
#'
#' Minimizing Copstress to obtain a clustered ratio, interval or ordinal PS configuration with given explicit power transformations theta. The function allows mix-and-match of explicit (via theta) and implicit (via type) transformations by setting the kappa, lambda, nu (or theta) and type arguments.
#'
#' This is an extremely flexible approach to least squares proximity scaling: It supports ratio power stress; ratio, interval and ordinal r stress and ratio, interval and ordinal MDS with or without a COPS penalty. Famous special cases of these models that can be fitted are multiscale MDS if kappa->0 and delta=log(delta), Alscal MDS (sstress) with lambda=kappa=2, sammon type mapping with weightmat=delta and nu=-1, elastic scaling with weightmat=delta and nu=-2. Due to mix-and-match this function also allows to fit models that have not yet been published, such as for example an "elastic scaling ordinal s-stress with cops penalty".
#'
#' If one wants to fit these models without the cops penalty, we recommend to use \code{\link[smacofx]{powerStressMin}} (for ratio and interval MDS with any power transformation for weights, dissimilarities and distances) or \code{\link[smacofx]{rStressMin}} (for ratio, interval and ordinal MDS with power transformations for distances and weights) as these use majorization.
#'
#' @rdname copstressMin
#'
#' @param delta numeric matrix or dist object of a matrix of proximities
#' @param kappa power transformation for fitted distances
#' @param lambda power transformation for proximities (only used if type="ratio" or "interval")
#' @param nu power transformation for weights
#' @param theta the theta vector of powers; the first is kappa (for the fitted distances if it exists), the second lambda (for the observed proximities if it exist and type="ratio" or "interval"), the third is nu (for the weights if it exists). If less than three elements are is given as argument, it will be recycled. Defaults to 1 1 1. Will override any kappa, lambda, nu parameters if they are given and do not match.
#' @param type what type of MDS to fit. Currently one of "ratio", "interval" or "ordinal". Default is "ratio".
#' @param ties the handling of ties for ordinal (nonmetric) MDS. Possible are "primary" (default), "secondary" or "tertiary".
#' @param weightmat (optional) a matrix of nonnegative weights; defaults to 1 for all off diagonals
#' @param ndim number of dimensions of the target space
#' @param init (optional) initial configuration
#' @param stressweight weight to be used for the fit measure; defaults to 0.975
#' @param cordweight weight to be used for the cordillera; defaults to 0.025
#' @param q the norm of the cordillera; defaults to 1
#' @param minpts the minimum points to make up a cluster in OPTICS, see \code{\link[dbscan]{optics}} where it is called \code{minPts}; defaults to ndim+1.
#' @param epsilon the epsilon parameter of OPTICS, the neighbourhood that is checked, see \code{\link[dbscan]{optics}}; defaults to 10 (which is plenty for the explicit normalization we use). Note this means we do not expect any noise objects per default. This number will rarely be exceeded if we standardize the configuration as is the default in cops. However if no standardization is applied or there is a procrustes adjustment to a configuration with variance of 10 or more on any of the axes, it can have the effect of being too small. In that case just set a much higher epsilon.
#' @param dmax The winsorization limit of reachability distances in the OPTICS Cordillera. If supplied, it should be either a numeric value that matches 'max(rang)' or 'NULL'; if 'NULL' it is found as 1.5 times (for kappa >1) or 1 times (for kappa <=1) the maximum reachbility value of the power torgerson model with the same lambda. If 'dmax' and 'rang' are supplied and 'dmax' is not 'max(rang)', a warning is given and 'rang' takes precedence.
#' @param rang range of the reachabilities to be considered. If missing it is found from the initial configuration by taking 0 as the lower boundary and dmax (see above) as upper boundary. See also \code{\link[cordillera]{cordillera}}
#' @param optimmethod What optimizer to use? Choose one string of 'Newuoa' (\code{\link[minqa]{newuoa}}), 'NelderMead' (see \code{\link[stats]{optim}}), 'hjk' (Hooke-Jeeves algorithm from \code{\link[dfoptim]{hjk}}), 'solnl' (from \code{\link[NlcOptim]{solnl}}), 'solnp' (from \code{\link[Rsolnp]{solnp}}), 'subplex' (from \code{\link[subplex]{subplex}}), 'SANN' (simulated annealing, \code{\link[stats]{optim}}), 'BFGS' (see \code{\link[stats]{optim}}), 'snomadr' (from \code{\link[crs]{snomadr}}), 'genoud' (from \code{\link[rgenoud]{genoud}}), 'gensa' (from \code{\link[GenSA]{GenSA}}), 'cmaes' (from \code{\link[cmaes]{cma_es}}) and 'direct' (from \code{\link[nloptr]{direct}}). See the linked functions for details on these solvers. There are also combinations that proved to work well good, like 'hjk-Newuoa', 'hjk-BFGS', 'BFGS-hjk', 'Newuoa-hjk', 'direct-Newuoa' and 'direct-BFGS'. Usually everything with 'hjk', 'BFGS', 'Newuoa', 'subplex' and 'solnl' in it work rather well in an acceptable time frame (depending on the smoothness of copstress). Default is 'hjk-Newuoa'.
#' @param verbose numeric value hat prints information on the fitting process; >2 is very verbose
#' @param normed should the Cordillera be normed; defaults to TRUE.
#' @param scale Scale the configuration (in MDS stress is invariant up to a scaling factor). One of "none" (so no extra scaling of the configuration but normalized to sum delta^2=1), "sd" (configuration divided by the highest standard deviation of any the columns), "proc" (procrustes adjustment to the initial fit) and "rmsq" (configuration divided by the maximum root mean square of the columns). Default is "sd" which often gives a nicer spread on the axes. Note that the scaled configuration is returned as $conf and the unscaled as $usconf, so manual calculation of the OC should be done with $conf.
#' @param accuracy numerical accuracy, defaults to 1e-7.
#' @param itmax maximum number of iterations. Defaults to 10000. For the two-step algorithms if itmax is exceeded by the first solver, the second algorithm is run for at least 0.1*itmax (so overall itmax may be exceeded by a factor of 1.1).
#' @param stresstype which stress to use in the copstress. Defaults to stress-1. If anything else is set, explicitly normed stress which is (stress-1)^2 is used. Using stress-1 puts more weight on MDS fit.
#' @param principal If ‘TRUE’, principal axis transformation is applied to the final configuration.
#' @param ... additional arguments to be passed to the optimization procedure
#'
#' @return A copsc object (inheriting from smacofP). A list with the components
#' \itemize{
#' \item delta: the original untransformed dissimilarities
#' \item tdelta: the explicitly transformed dissimilarities
#' \item dhat: the explicitly transformed dissimilarities (dhats), optimally scaled and normalized (which are approximated by the fit)
#' \item confdist: Configuration distances, the transformed fitted distances
#' \item conf: the configuration (normed) and scaled as specified in scale.
#' \item usconf: the unscaled configuration (normed to sum delta^2=1). Scaling applied to usconf gives conf.
#' \item parameters, par, pars : the theta vector of powers tranformations (kappa, lambda, nu)
#' \item niter: number of iterations of the optimizer.
#' \item stress: the square root of explicitly normalized stress (calculated for confo).
#' \item spp: stress per point
#' \item ndim: number of dimensions
#' \item model: Fitted model name
#' \item call: the call
#' \item nobj: the number of objects
#' \item type, loss, losstype: stresstype
#' \item stress.m: The stress used for copstress. If stresstype="stress-1" this is like $stress else it is stress^2
#' \item copstress: the copstress loss value
#' \item resmat: the matrix of residuals
#' \item weightmat: the matrix of untransformed weights
#' \item tweightmat: the transformed weighting matrix (here weightmat^nu)
#' \item OC: the (normed) OPTICS Cordillera object (calculated for scaled conf)
#' \item OCv: the (normed) OPTICS Cordillera value alone (calculated for scaled conf)
#' \item optim: the object returned from the optimization procedure
#' \item stressweight, cordweight: the weights of the stress and OC respectively (v_1 and v_2)
#' \item optimmethod: The solver used
#' \item type: the type of MDS fitted
#'}
#'
#'
#'
#' @details
#' Some optimizers (including the default hjk-Newuoa) will print a warning if itmax is (too) small or if there was no convergence. Consider increasing itmax then.
#'
#' For some solvers theresometimes may be an error [NA/NaN/Inf in foreign function call (arg 3)] stemming from smacof::transform(). This happens when the algorithm places two object at exactly the same place so their fitted distance is 0. This is good from an OPTICS Cordillera point of view (as it is more clustered) which is why some solvers like to pick that up, but it can lead to an issue in the optimal scaling in smacof. This can usually be mitigated when specifying the model by either using less cordweight, less itmax, less accuracy or combining the two offending objects into one (so include them as a combined row in the distance matrix).
#'
#' We might eventually switch to newuoa in nloptr.
#'
#' @examples
#' dis<-as.matrix(smacof::kinshipdelta)
#'
#' set.seed(1)
#' ## Copstress with equal weight to stress and cordillera
#' res1<-copstressMin(dis,stressweight=0.5,cordweight=0.5,
#' itmax=100) #use higher itmax about 10000
#' res1
#' summary(res1)
#' plot(res1) #super clustered
#'
#' ##Alias name
#' res1<-copsc(dis,stressweight=0.5,
#' cordweight=0.5,itmax=100)
#'
#'
#' ## Elastic scaling ordinal s-stress with cops penalty
#' res1<-copsc(dis,type="ordinal",kappa=2,nu=-2,weightmat=dis,
#' stressweight=0.5, cordweight=0.5,itmax=100)
#'
#'
#' @import cordillera
#' @importFrom utils tail
#' @importFrom stats dist as.dist optim sd
#' @importFrom dfoptim hjk
#' @importFrom NlcOptim solnl
#' @importFrom Rsolnp solnp
#' @importFrom subplex subplex
#' @importFrom crs snomadr
#' @importFrom cmaes cma_es
#' @importFrom rgenoud genoud
#' @importFrom GenSA GenSA
#' @importFrom nloptr direct
#' @importFrom minqa newuoa
#' @importFrom smacof transPrep transform
#' @importFrom smacofx spp
#'
#'
#' @keywords clustering multivariate
#' @export
copstressMin <- function (delta, kappa=1, lambda=1, nu=1, theta=c(kappa,lambda,nu), type=c("ratio","interval","ordinal"), ties="primary", weightmat=1-diag(nrow(delta)), ndim = 2, init=NULL, stressweight=0.975,cordweight=0.025,q=1,minpts=ndim+1,epsilon=max(10,max(delta)),dmax=NULL,rang,optimmethod=c("NelderMead","Newuoa","BFGS","SANN","hjk","solnl","solnp","subplex","snomadr","hjk-Newuoa","hjk-BFGS","BFGS-hjk","Newuoa-hjk","cmaes","direct","direct-Newuoa","direct-BFGS","genoud","gensa"),verbose=0,scale=c("sd","rmsq","proc","none"),normed=TRUE, accuracy = 1e-7, itmax = 10000, stresstype=c("stress-1","stress"),principal=FALSE,...)
{
if(inherits(delta,"dist") || is.data.frame(delta)) delta <- as.matrix(delta)
if(!isSymmetric(delta)) stop("Delta is not symmetric.\n")
if(inherits(weightmat,"dist") || is.data.frame(weightmat)) weightmat <- as.matrix(weightmat)
if(!isSymmetric(weightmat)) stop("weightmat is not symmetric.\n")
## -- Setup for MDS type
if(missing(type)) type <- "ratio"
type <- match.arg(type, c("ratio", "interval", "ordinal"),several.ok = FALSE)
trans <- type
typo <- type
if (trans=="ratio"){
trans <- "none"
}
#TODO: if we want other optimal scalings as well
else if (trans=="ordinal" & ties=="primary"){
trans <- "ordinalp"
typo <- "ordinal (primary)"
} else if(trans=="ordinal" & ties=="secondary"){
trans <- "ordinals"
typo <- "ordinal (secondary)"
} else if(trans=="ordinal" & ties=="tertiary"){
trans <- "ordinalt"
typo <- "ordinal (tertiary)"
#} else if(trans=="spline"){
# trans <- "mspline"
}
if(type =="ordinal") theta <- c(kappa,1,nu) #We dont allow powers for dissimilarities in nonmetric MDS
kappa <- theta[1]
lambda <- theta[2]
nu <- theta[3]
#plot <- FALSE
n <- nrow(delta)
if (ndim > (n - 1)) stop("Maximum number of dimensions is n-1!")
if(is.null(rownames(delta))) rownames(delta) <- 1:n
labos <- rownames(delta)
if(verbose>0) cat(paste("Minimizing",type,"copstress with kappa=",kappa,"lambda=",lambda,"nu=",nu,".\n"))
if(missing(optimmethod)) optimmethod <- "hjk-Newuoa"
if(missing(scale)) scale <- "sd"
if(missing(stresstype)) stresstype <- "stress-1"
##-- Prepare for dissimilarity scaling
r <- kappa/2
deltaorig <- delta
delta <- delta^lambda
weightmato <- weightmat
weightmat <- weightmat^nu
weightmat[!is.finite(weightmat)] <- 0 #new
deltaold <- delta
disobj <- smacof::transPrep(as.dist(delta), trans = trans, spline.intKnots = 2, spline.degree = 2)#spline.intKnots = spline.intKnots, spline.degree = spline.degree) #FIXME: only works with dist() style object
## Add an intercept to the spline base transformation
#if (trans == "mspline") disobj$base <- cbind(rep(1, nrow(disobj$base)), disobj$base)
delta <- delta / enorm (delta, weightmat) #normalize to sum=1
## --- starting rang if not given
if(missing(rang))
#perhaps put this into the optimization function?
{
if(is.null(dmax))
{
if(is.null(init))
{
if(verbose>1) cat ("Fitting configuration for finding rang argument. \n")
initsol <- smacof::torgerson(delta,p=ndim)
init0 <- initsol
} else init0 <- init
init0 <- init0/enorm(init0)
# if(scale=="std") init0 <- scale(init0) #standardizes config before cordillera
if(scale=="none") init0 <- init0
if(scale=="sd") #scales config to sd=1 for most spread dimension before cordillera
{
init0 <- init0/max(apply(init0,2,stats::sd))
}
if(scale=="rmsq") #scales config to rmsq=1 for most spread dimension before cordillera
{
testso <- base::scale(init0,center=FALSE)
init0 <- init0/max(attr(testso,"scaled:scale"))
}
if(scale=="proc") #scales config by procrusting to init
{
if(is.null(init)) init <- initsol
procr <- smacof::Procrustes(init,init0)
init0 <- procr$Yhat
}
crp <- cordillera::cordillera(init0,q=q,minpts=minpts,epsilon=epsilon,scale=FALSE)$reachplot
cin <- max(crp)
dmax <- ifelse(kappa>1,1.5*cin,1.1*cin)
}
rang <- c(0,dmax) #This sets the range to (0,dmax)
if(verbose>1) cat("rang is",rang,".\n")
}
if(isTRUE(max(rang)!=dmax)) warning("Note: The supplied dmax and rang do not match. I took supplied rang as rang.\n")
## --- starting values
if(is.null(init))
{
if(exists("init0")) init <- init0 else init <- smacof::torgerson(delta,p=ndim)
#was: cops::powerStressFast(delta,kappa=kappa,lambda=lambda,nu=nu,ndim=ndim)$conf
}
xstart <- xold <- init
xold <- xold/enorm(xold)
dhat2 <- NA
#labs <- row.names(delta)
copsf <- function(x,delta,disobj,r,n,ndim,weightmat,stressweight,cordweight,q,minpts,epsilon,rang,scale,normed,init,...)
{
#init is used here only for Procrustes
if(!is.matrix(x)) x <- matrix(x,ncol=ndim)
delta <- delta/enorm(delta,weightmat)
x <- x/enorm(x)
dnew <- sqdist (x)
e <- as.dist(mkPower(dnew,r)) #was bug prior to 1.12-1
dhat2 <<- smacof::transform(e, disobj, w = as.dist(weightmat), normq = 0.5) ##I use <<- to change this also in the parent environment because the copsf function should only return a scalar for the optimizers but I the last dhats2 and also delta later on
dhatt <- dhat2$res
dhatd <- structure(dhatt, Size = n, call = quote(as.dist.default(m=b)), class = "dist", Diag = FALSE, Upper = FALSE)
delta <<- as.matrix(dhatd) ##I use <<- to save this also in the parent environment
rnew <- sum (weightmat * delta * mkPower (dnew, r))
nnew <- sum (weightmat * mkPower (dnew, 2*r))
anew <- rnew / nnew
stressi <- 1 - 2 * anew * rnew + (anew ^ 2) * nnew
if(stresstype=="stress-1") stressi <- sqrt(stressi)
if(scale=="none") x <- x
#if(scale=="std") x <- base::scale(x) #standardizes config before cordillera
if(scale=="sd") #scales config to sd=1 for most spread dimension before cordillera
{
x <- x/max(apply(x,2,stats::sd))
}
if(scale=="rmsq") #scales config to rmsq=1 for most spread dimension before cordillera
{
testso <- base::scale(x,center=FALSE)
x <- x/max(attr(testso,"scaled:scale"))
}
if(scale=="proc") #scales config by procrusting to init
{
procr <- smacof::Procrustes(init,x)
x <- procr$Yhat
}
corrd <- cordillera::cordillera(x,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=FALSE,...)
struc <- corrd$raw
if(normed) {
struc <- corrd$normed
}
ic <- stressweight*stressi - cordweight*struc
if(verbose>1) cat("copstress =",ic,"mdsloss =",stressi,"OC =",struc,"minpts=",minpts,"kappa =",kappa,"lambda =",lambda,"nu=",nu,"\n")
ic
}
if(verbose>1) cat("Starting Minimization with",optimmethod,":\n")
if(optimmethod=="Newuoa") {
suppressWarnings(optimized <- minqa::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmax,rhoend=accuracy,iprint=verbose-2),...))
itel <- itel+optimized$feval
ovalue <-optimized$fval
#optimized <- nloptr::newuoa(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmax,rhoend=accuracy,iprint=verbose-2),...)
#xnew <- matrix(optimized$par,ncol=ndim)
#itel <- optimized$iter
#ovalue <-optimized$value
}
if(optimmethod=="direct") {
xold <- as.vector(xold)
optimized <- nloptr::direct(function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),lower=rep(5*min(xold),length(xold)),upper=rep(5*max(xold),length(xold)),nl.info=isTRUE(verbose>1),control=list(maxeval=itmax,xtol_rel=accuracy),...)
# optimized <- nloptr::direct(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,#n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=),#lower=rep(5*min(xold),length(xold)),upper=rep(5*max(xold),length(xold)), nl.info=isTRUE(verbose>1),control=list(maxeval=itmax,xtol_rel=accuracy),...)
xnew <- matrix(optimized$par,ncol=ndim)
itel <- optimized$iter
ovalue <-optimized$value
}
if(optimmethod=="direct-Newuoa") {
xold <- as.vector(xold)
optimized1 <- nloptr::direct(function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),lower=rep(5*min(xold),length(xold)),upper=rep(5*max(xold),length(xold)),nl.info=isTRUE(verbose>1),control=list(maxeval=itmax,xtol_rel=accuracy),...)
xnew <- optimized1$par
itel1 <- optimized1$iter
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
suppressWarnings(optimized <- minqa::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmaxreduced,rhoend=accuracy,iprint=verbose-2),...))
itel <- itel1+optimized$feval
ovalue <-optimized$fval
#optimized <- nloptr::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),nl.info=isTRUE(verbose>2),control=list(maxeval=itmaxreduced,xtol_rel=accuracy),...)
#xnew <- matrix(optimized$par,ncol=ndim)
#itel <- itel1+optimized$iter
#ovalue <-optimized$value
}
if(optimmethod=="direct-BFGS") {
xold <- as.vector(xold)
optimized1 <- nloptr::direct(function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),lower=rep(5*min(xold),length(xold)),upper=rep(5*max(xold),length(xold)),nl.info=isTRUE(verbose>1),control=list(maxeval=itmax,xtol_rel=accuracy),...)
xnew <- matrix(optimized1$par,ncol=ndim)
itel1 <- optimized1$iter
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
optimized <- optim(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),method="BFGS",control=list(maxit=itmaxreduced,trace=0,reltol=accuracy),...)
xnew <- optimized$par
itel <- itel1+optimized$counts[[1]]
ovalue <-optimized$val
}
if(optimmethod=="genoud") {
optimized <- rgenoud::genoud(function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),starting.values=as.vector(xold),nvars=length(xold),print.level=verbose,...)
xnew <- matrix(optimized$par,ncol=ndim)
itel <- optimized$generations
ovalue <-optimized$value
}
if(optimmethod=="gensa") {
optimized <- GenSA::GenSA(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),lower=rep(5*min(xold),length(xold)),upper=rep(5*max(xold),length(xold)),control=list(max.call=itmax,verbose=isTRUE(verbose>1),smooth=FALSE),...)
xnew <- matrix(optimized$par,ncol=ndim)
itel <- optimized$counts
ovalue <-optimized$value
}
if(optimmethod=="cmaes") {
xold <- as.vector(xold)
optimized <- cmaes::cma_es(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxit=itmax),...)
xnew <- matrix(optimized$par,ncol=ndim)
itel <- optimized$counts[1]
ovalue <-optimized$value
}
if(optimmethod=="cmaes-Newuoa") {
xold <- as.vector(xold)
optimized1 <- cmaes::cma_es(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxit=itmax),...)
xnew <- matrix(optimized1$par,ncol=ndim)
itel1 <- optimized1$counts[1]
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
suppressWarnings(optimized <- minqa::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmaxreduced,rhoend=accuracy,iprint=verbose-2),...))
xnew <- matrix(optimized$par,ncol=ndim)
itel <- itel1+optimized$feval
ovalue <-optimized$fval
#optimized <- nloptr::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),nl.info=isTRUE(verbose>2),control=list(maxeval=itmaxreduced,xtol_rel=accuracy),...)
#xnew <- matrix(optimized$par,ncol=ndim)
#itel <- itel1+optimized$iter
#ovalue <-optimized$value
}
if(optimmethod=="Newuoa-cmaes") {
suppressWarnings(optimized1 <- minqa::newuoa(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmax,rhoend=accuracy,iprint=verbose-2),...))
#optimized1 <- nloptr::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),nl.info=isTRUE(verbose>2),control=list(maxeval=itmax,xtol_rel=accuracy),...)
xnew <- optimized1$par
itel1 <- optimized1$feval
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
optimized <- cmaes::cma_es(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxit=itmaxreduced),...)
xnew <- matrix(optimized$par,ncol=ndim)
itel <- optimized$counts[1]+itel1
ovalue <-optimized$value
}
if(optimmethod=="NelderMead") {
optimized <- optim(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxit=itmax,trace=0,reltol=accuracy),...)
xnew <- optimized$par
itel <- optimized$counts[[1]]
ovalue <-optimized$val
}
if(optimmethod=="BFGS") {
optimized <- optim(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),method="BFGS",control=list(maxit=itmax,trace=0,reltol=accuracy),...)
xnew <- optimized$par
itel <- optimized$counts[[1]]
ovalue <-optimized$val
}
if(optimmethod=="SANN") {
optimized <- optim(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),method="SANN",control=list(maxit=itmax,trace=0,reltol=accuracy),...)
xnew <- optimized$par
itel <- optimized$counts[[1]]
ovalue <-optimized$val
}
if(optimmethod=="hjk") {
optimized <- dfoptim::hjk(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfeval=itmax),...)
xnew <- optimized$par
itel <- optimized$feval
ovalue <-optimized$value
}
if(optimmethod=="hjk-Newuoa") { #twostep1
#cat("Before HJK Iterations:", itmax,"\n")
optimized1 <- dfoptim::hjk(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfeval=itmax),...)
xnew <- optimized1$par
itel1 <- optimized1$feval
#cat("OpVal after HJK:",optimized1$value,"\n")
#cat("After HJK Iterations:", itel1,"\n")
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
#cat("Before Newuoa Iterations:", itmaxreduced,"\n")
suppressWarnings(optimized <- minqa::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmaxreduced,rhoend=accuracy,iprint=verbose-2),...))
xnew <- matrix(optimized$par,ncol=ndim)
itel <- itel1+optimized$feval
ovalue <-optimized$fval
#optimized <- nloptr::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),nl.info=isTRUE(verbose>2),control=list(maxeval=itmaxreduced,xtol_rel=accuracy),...)
##BUG: Was init= before in all Newuoa (no argument). Changed in 1.3-6
#xnew <- matrix(optimized$par,ncol=ndim) #
##test function for return value fr <- function(x) { sum((x - matrix(c(8,5,3,4,6,7),ncol=2))^2) } if returned is matrix(c(8,5,3,4,6,7),ncol=2) it is cool
#cat("Newuoa Iterations:", optimized$iter,"\n")
#cat("Newuoa Iterations:", optimized$feval,"\n")
#itel <- itel1+optimized$iter
#cat("After Newuoa Iterations:", itel,"\n")
#ovalue <-optimized$value
#cat("OpVal after all:",ovalue,"\n")
}
if(optimmethod=="hjk-BFGS") { #twostep2
optimized1 <- dfoptim::hjk(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfeval=itmax,trace=0),...)
xnew <- optimized1$par
itel1 <- optimized1$feval
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
optimized <- optim(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),method="BFGS",control=list(maxit=itmaxreduced,trace=0,reltol=accuracy),...)
xnew <- optimized$par
itel <- itel1+optimized$counts[[1]]
ovalue <-optimized$val
}
if(optimmethod=="BFGS-hjk") { #twostep6
optimized1 <- optim(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),method="BFGS",control=list(maxit=itmax,trace=0),...)
xnew <- optimized1$par
itel1 <- optimized1$counts[[1]]
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
optimized <- dfoptim::hjk(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfeval=itmaxreduced,trace=0,tol=accuracy),...)
xnew <- optimized$par
itel <- itel1+optimized$feval
ovalue <-optimized$val
}
if(optimmethod=="BFGS-Newuoa") { #twostep5
optimized1 <- optim(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),method="BFGS",control=list(maxfeval=itmax,trace=0),...)
itel1 <- optimized1$counts[[1]]
xnew <- optimized1$par
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
suppressWarnings(optimized <- minqa::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmaxreduced,rhoend=accuracy,iprint=verbose-2),...))
xnew <- matrix(optimized$par,ncol=ndim)
itel <- itel1+optimized$feval
ovalue <-optimized$fval
#optimized <- nloptr::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),nl.info=isTRUE(verbose>2),control=list(maxeval=itmaxreduced,xtol_rel=accuracy),...)
#xnew <- matrix(optimized$par,ncol=ndim) #
#itel <- itel1+optimized$iter
#ovalue <-optimized$value
}
if(optimmethod=="hjk-solnl") {#twostep3
optimized1 <- dfoptim::hjk(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfeval=itmax,trace=0),...)
xnew <- optimized1$par
itel1 <- optimized1$feval
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
optimized <- NlcOptim::solnl(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),maxnFun=itmaxreduced,tolFun=accuracy,...)
xnew <- optimized$par
itel <- itel1+optimized$counts[[1]]
ovalue <-optimized$fn
}
if(optimmethod=="hjk-subplex") {#twostep4
optimized1 <- dfoptim::hjk(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfeval=itmax,trace=0),...)
xnew <- optimized1$par
itel1 <- optimized1$feval
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
optimized <- subplex::subplex(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxit=itmaxreduced,reltol=accuracy),...)
xnew <- matrix(optimized$par,ncol=ndim)
itel <- itel1+optimized$count
ovalue <-optimized$value
}
# if(optimmethod=="dfsane") {
# optimized <- BB::dfsane(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxit=itmax,trace=0,tol=accuracy),...)
# xnew <- optimized$par
# itel <- optimized$feval
# ovalue <-optimized$residual
# }
if(optimmethod=="solnl") {
optimized <- NlcOptim::solnl(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),maxnFun=itmax,tolFun=accuracy,...)
xnew <- optimized$par
itel <- optimized$counts[[1]]
ovalue <-optimized$fn
}
## if(optimmethod=="isres") {
## optimized <- nloptr::isres(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),maxeval=itmax,trace=verbose-2,xtol_rel=accuracy,...)
## xnew <- optimized$par
## itel <- optimized$iter
## ovalue <-optimized$value
## }
if(optimmethod=="solnp") {
optimized <- Rsolnp::solnp(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(outer.iter=itmax,trace=0,tol=accuracy),...)
xnew <- matrix(optimized$pars,ncol=ndim)
itel <- optimized$nfuneval
ovalue <-utils::tail(optimized$values,1)
}
if(optimmethod=="subplex") {
optimized <- subplex::subplex(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxit=itmax,reltol=accuracy),...)
xnew <- matrix(optimized$par,ncol=ndim)
itel <- optimized$count
ovalue <-optimized$value
}
if(optimmethod=="snomadr") {
copsf2 <- function(x,params)
{
copsf(x,delta=params[[1]],disobj=params[[2]],r=params[[3]],n=params[[4]],ndim=params[[5]],weightmat=params[[6]],stressweight=params[[7]],cordweight=params[[8]],q=params[[9]],minpts=params[[10]],epsilon=params[[11]],rang=params[[12]],scale=params[[13]],normed=params[[14]],init=params[[15]])
}
params <- list(delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init)
optimized <- crs::snomadr(copsf2,params=params,,n=dim(xold)[1],x0=xold,print.output=isTRUE(verbose-2>0),...)
xnew <- optimized$solution
itel <- optimized$iterations
ovalue <-optimized$objective
}
# if(optimmethod=="snomadr-Newuoa") {
# optimized1 <- crs::snomadr(function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),n=dim(xold)[1],x0=xold,print.output=isTRUE(verbose-2>0),...)
# xnew <- optimized1$solution
# itel <- optimized1$iterations
# optimized <- minqa::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmax-itel,rhoend=accuracy,iprint=verbose-2),...)
# xnew <- matrix(optimized$par,ncol=ndim)
# itel <- itel+optimized$feval
# ovalue <-optimized$fval
# }
# if(optimmethod=="snomadr-BFGS") {
# optimized1 <- crs::snomadr(function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=w#eightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),n=dim(xold)[1],x0=xold,print.output=isTRUE(verbose-2>0),...)
# xnew <- optimized1$solution
# itel <- optimized1$iterations
# optimized <- optim(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),method="BFGS",control=list(maxit=itmax-itel,trace=0,reltol=accuracy),...)
# xnew <- optimized$par
# itel <- itel+optimized$counts[[1]]
# ovalue <-optimized$val
# }
xnew <- xnew/enorm(xnew)
rownames(delta) <- labos
dnew <- sqdist (xnew)
rnew <- sum (weightmat * delta * mkPower (dnew, r))
nnew <- sum (weightmat * mkPower (dnew, 2*r))
anew <- rnew / nnew
stress.m <- snew <- 1 - 2 * anew * rnew + (anew ^ 2) * nnew
stress <- sqrt(stress.m)
if(stresstype=="stress-1") {
stress.m <- sqrt(stress.m)
}
attr(xnew,"dimnames")[[1]] <- rownames(delta)
attr(xnew,"dimnames")[[2]] <- paste("D",1:ndim,sep="")
#doutm <- (2*sqrt(sqdist(xnew)))^kappa #fitted powered euclidean distance but times two
#doutm <- as.matrix(dist(xnew)^kappa)
doutm <- mkPower(sqdist(xnew),r)
#deltam <- delta
#deltaorigm <- deltaorig
#deltaoldm <- deltaold
delta <- stats::as.dist(delta)
deltaorig <- stats::as.dist(deltaorig)
deltaold <- stats::as.dist(deltaold)
#doute <- doutm/enorm(doutm)
#doute <- stats::as.dist(doute)
dout <- stats::as.dist(doutm)
#resmat <- as.matrix((delta - doute)^2)
#spp <- colMeans(resmat)
#weightmatm <-weightmat
weightmat <- stats::as.dist(weightmat)
weightmato <- stats::as.dist(weightmato)
spoint <- smacofx::spp(delta, dout, weightmat)
resmat<-spoint$resmat
rss <- sum(spoint$resmat[lower.tri(spoint$resmat)])
spp <- spoint$spp
#stressen <- sum(weightmat*(delta-doute)^2)#raw stress on the normalized proximities and normalized distances
#if(scale=="std") xnews <- base::scale(xnew) #standardizes config before cordillera
if(scale=="sd") #scales config to sd=1 for most spread dimension before cordillera
{
xnews <- xnew/max(apply(xnew,2,stats::sd))
}
if(scale=="rmsq") #scales config to rmsq=1 for most spread dimension before cordillera
{
testso <- base::scale(xnew,center=FALSE)
xnews <- xnew/max(attr(testso,"scaled:scale"))
}
if(scale=="proc") #scales config by procrusting to init
{
procr <- smacof::Procrustes(init,xnew)
xnews <- procr$Yhat
}
if(scale=="none") xnews <- xnew #no standardisation
if (principal) {
xnew_svd <- svd(xnew)
xnew <- xnew %*% xnew_svd$v
xnews_svd <- svd(xnews)
xnews <- xnews %*% xnews_svd$v
}
if(verbose>0) cat("*** Stress:",stress.m,"; Stress-1:",stress,"; from optimization: ",ovalue,"\n")
#like smacofP
out <- list(delta=deltaorig, tdelta=deltaold, dhat=delta, confdist=dout, iord=dhat2$iord.prim, conf = xnews, stress=stress, spp=spp, ndim=ndim, weightmat=weightmato, resmat=resmat, rss=rss, init=xstart, model="COPS-C", niter = itel, nobj = n, type = type, call=match.call(), stress.m=snew, alpha = anew, sigma = snew, pars=c(kappa=kappa,lambda=lambda,nu=nu),parameters=c(kappa=kappa,lambda=lambda,nu=nu),theta=c(kappa=kappa,lambda=lambda,nu=nu),tweightmat=weightmat)
#extra slots for class copsc
out$gamma <- NA
out$ties <- ties
out$usconf <- xnew
#out$loss <- "copstress"
out$OC <- cordillera::cordillera(out$conf,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=FALSE)
out$OCv <- ifelse(normed,out$OC$normed,out$OC$raw)
out$copstress <- ovalue
out$optim <- optimized
out$stressweight <- stressweight
out$cordweight <- cordweight
out$optimethod <- optimmethod
#out$typo <- typo
#out$nobj <- dim(out$conf)[1]
class(out) <- c("copsc","cops","smacofP","smacofB","smacof")
out
}
#' @rdname copstressMin
#' @export
copsc <- copstressMin
#' @rdname copstressMin
#' @export
copStressMin <- copstressMin
Try the cops package in your browser
Any scripts or data that you put into this service are public.
cops documentation built on Sept. 30, 2024, 9:22 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#' Fitting a COPS-C Model (COPS Variant 1).
#'
#' Minimizing Copstress to obtain a clustered ratio, interval or ordinal PS configuration with given explicit power transformations theta. The function allows mix-and-match of explicit (via theta) and implicit (via type) transformations by setting the kappa, lambda, nu (or theta) and type arguments.
#'
#' This is an extremely flexible approach to least squares proximity scaling: It supports ratio power stress; ratio, interval and ordinal r stress and ratio, interval and ordinal MDS with or without a COPS penalty. Famous special cases of these models that can be fitted are multiscale MDS if kappa->0 and delta=log(delta), Alscal MDS (sstress) with lambda=kappa=2, sammon type mapping with weightmat=delta and nu=-1, elastic scaling with weightmat=delta and nu=-2. Due to mix-and-match this function also allows to fit models that have not yet been published, such as for example an "elastic scaling ordinal s-stress with cops penalty".
#'
#' If one wants to fit these models without the cops penalty, we recommend to use \code{\link[smacofx]{powerStressMin}} (for ratio and interval MDS with any power transformation for weights, dissimilarities and distances) or \code{\link[smacofx]{rStressMin}} (for ratio, interval and ordinal MDS with power transformations for distances and weights) as these use majorization.
#'
#' @rdname copstressMin
#'
#' @param delta numeric matrix or dist object of a matrix of proximities
#' @param kappa power transformation for fitted distances
#' @param lambda power transformation for proximities (only used if type="ratio" or "interval")
#' @param nu power transformation for weights
#' @param theta the theta vector of powers; the first is kappa (for the fitted distances if it exists), the second lambda (for the observed proximities if it exist and type="ratio" or "interval"), the third is nu (for the weights if it exists). If less than three elements are is given as argument, it will be recycled. Defaults to 1 1 1. Will override any kappa, lambda, nu parameters if they are given and do not match.
#' @param type what type of MDS to fit. Currently one of "ratio", "interval" or "ordinal". Default is "ratio".
#' @param ties the handling of ties for ordinal (nonmetric) MDS. Possible are "primary" (default), "secondary" or "tertiary".
#' @param weightmat (optional) a matrix of nonnegative weights; defaults to 1 for all off diagonals
#' @param ndim number of dimensions of the target space
#' @param init (optional) initial configuration
#' @param stressweight weight to be used for the fit measure; defaults to 0.975
#' @param cordweight weight to be used for the cordillera; defaults to 0.025
#' @param q the norm of the cordillera; defaults to 1
#' @param minpts the minimum points to make up a cluster in OPTICS, see \code{\link[dbscan]{optics}} where it is called \code{minPts}; defaults to ndim+1.
#' @param epsilon the epsilon parameter of OPTICS, the neighbourhood that is checked, see \code{\link[dbscan]{optics}}; defaults to 10 (which is plenty for the explicit normalization we use). Note this means we do not expect any noise objects per default. This number will rarely be exceeded if we standardize the configuration as is the default in cops. However if no standardization is applied or there is a procrustes adjustment to a configuration with variance of 10 or more on any of the axes, it can have the effect of being too small. In that case just set a much higher epsilon.
#' @param dmax The winsorization limit of reachability distances in the OPTICS Cordillera. If supplied, it should be either a numeric value that matches 'max(rang)' or 'NULL'; if 'NULL' it is found as 1.5 times (for kappa >1) or 1 times (for kappa <=1) the maximum reachbility value of the power torgerson model with the same lambda. If 'dmax' and 'rang' are supplied and 'dmax' is not 'max(rang)', a warning is given and 'rang' takes precedence.
#' @param rang range of the reachabilities to be considered. If missing it is found from the initial configuration by taking 0 as the lower boundary and dmax (see above) as upper boundary. See also \code{\link[cordillera]{cordillera}}
#' @param optimmethod What optimizer to use? Choose one string of 'Newuoa' (\code{\link[minqa]{newuoa}}), 'NelderMead' (see \code{\link[stats]{optim}}), 'hjk' (Hooke-Jeeves algorithm from \code{\link[dfoptim]{hjk}}), 'solnl' (from \code{\link[NlcOptim]{solnl}}), 'solnp' (from \code{\link[Rsolnp]{solnp}}), 'subplex' (from \code{\link[subplex]{subplex}}), 'SANN' (simulated annealing, \code{\link[stats]{optim}}), 'BFGS' (see \code{\link[stats]{optim}}), 'snomadr' (from \code{\link[crs]{snomadr}}), 'genoud' (from \code{\link[rgenoud]{genoud}}), 'gensa' (from \code{\link[GenSA]{GenSA}}), 'cmaes' (from \code{\link[cmaes]{cma_es}}) and 'direct' (from \code{\link[nloptr]{direct}}). See the linked functions for details on these solvers. There are also combinations that proved to work well good, like 'hjk-Newuoa', 'hjk-BFGS', 'BFGS-hjk', 'Newuoa-hjk', 'direct-Newuoa' and 'direct-BFGS'. Usually everything with 'hjk', 'BFGS', 'Newuoa', 'subplex' and 'solnl' in it work rather well in an acceptable time frame (depending on the smoothness of copstress). Default is 'hjk-Newuoa'.
#' @param verbose numeric value hat prints information on the fitting process; >2 is very verbose
#' @param normed should the Cordillera be normed; defaults to TRUE.
#' @param scale Scale the configuration (in MDS stress is invariant up to a scaling factor). One of "none" (so no extra scaling of the configuration but normalized to sum delta^2=1), "sd" (configuration divided by the highest standard deviation of any the columns), "proc" (procrustes adjustment to the initial fit) and "rmsq" (configuration divided by the maximum root mean square of the columns). Default is "sd" which often gives a nicer spread on the axes. Note that the scaled configuration is returned as $conf and the unscaled as $usconf, so manual calculation of the OC should be done with $conf.
#' @param accuracy numerical accuracy, defaults to 1e-7.
#' @param itmax maximum number of iterations. Defaults to 10000. For the two-step algorithms if itmax is exceeded by the first solver, the second algorithm is run for at least 0.1*itmax (so overall itmax may be exceeded by a factor of 1.1).
#' @param stresstype which stress to use in the copstress. Defaults to stress-1. If anything else is set, explicitly normed stress which is (stress-1)^2 is used. Using stress-1 puts more weight on MDS fit.
#' @param principal If ‘TRUE’, principal axis transformation is applied to the final configuration.
#' @param ... additional arguments to be passed to the optimization procedure
#'
#' @return A copsc object (inheriting from smacofP). A list with the components
#' \itemize{
#' \item delta: the original untransformed dissimilarities
#' \item tdelta: the explicitly transformed dissimilarities
#' \item dhat: the explicitly transformed dissimilarities (dhats), optimally scaled and normalized (which are approximated by the fit)
#' \item confdist: Configuration distances, the transformed fitted distances
#' \item conf: the configuration (normed) and scaled as specified in scale.
#' \item usconf: the unscaled configuration (normed to sum delta^2=1). Scaling applied to usconf gives conf.
#' \item parameters, par, pars : the theta vector of powers tranformations (kappa, lambda, nu)
#' \item niter: number of iterations of the optimizer.
#' \item stress: the square root of explicitly normalized stress (calculated for confo).
#' \item spp: stress per point
#' \item ndim: number of dimensions
#' \item model: Fitted model name
#' \item call: the call
#' \item nobj: the number of objects
#' \item type, loss, losstype: stresstype
#' \item stress.m: The stress used for copstress. If stresstype="stress-1" this is like $stress else it is stress^2
#' \item copstress: the copstress loss value
#' \item resmat: the matrix of residuals
#' \item weightmat: the matrix of untransformed weights
#' \item tweightmat: the transformed weighting matrix (here weightmat^nu)
#' \item OC: the (normed) OPTICS Cordillera object (calculated for scaled conf)
#' \item OCv: the (normed) OPTICS Cordillera value alone (calculated for scaled conf)
#' \item optim: the object returned from the optimization procedure
#' \item stressweight, cordweight: the weights of the stress and OC respectively (v_1 and v_2)
#' \item optimmethod: The solver used
#' \item type: the type of MDS fitted
#'}
#'
#'
#'
#' @details
#' Some optimizers (including the default hjk-Newuoa) will print a warning if itmax is (too) small or if there was no convergence. Consider increasing itmax then.
#'
#' For some solvers theresometimes may be an error [NA/NaN/Inf in foreign function call (arg 3)] stemming from smacof::transform(). This happens when the algorithm places two object at exactly the same place so their fitted distance is 0. This is good from an OPTICS Cordillera point of view (as it is more clustered) which is why some solvers like to pick that up, but it can lead to an issue in the optimal scaling in smacof. This can usually be mitigated when specifying the model by either using less cordweight, less itmax, less accuracy or combining the two offending objects into one (so include them as a combined row in the distance matrix).
#'
#' We might eventually switch to newuoa in nloptr.
#'
#' @examples
#' dis<-as.matrix(smacof::kinshipdelta)
#'
#' set.seed(1)
#' ## Copstress with equal weight to stress and cordillera
#' res1<-copstressMin(dis,stressweight=0.5,cordweight=0.5,
#' itmax=100) #use higher itmax about 10000
#' res1
#' summary(res1)
#' plot(res1) #super clustered
#'
#' ##Alias name
#' res1<-copsc(dis,stressweight=0.5,
#' cordweight=0.5,itmax=100)
#'
#'
#' ## Elastic scaling ordinal s-stress with cops penalty
#' res1<-copsc(dis,type="ordinal",kappa=2,nu=-2,weightmat=dis,
#' stressweight=0.5, cordweight=0.5,itmax=100)
#'
#'
#' @import cordillera
#' @importFrom utils tail
#' @importFrom stats dist as.dist optim sd
#' @importFrom dfoptim hjk
#' @importFrom NlcOptim solnl
#' @importFrom Rsolnp solnp
#' @importFrom subplex subplex
#' @importFrom crs snomadr
#' @importFrom cmaes cma_es
#' @importFrom rgenoud genoud
#' @importFrom GenSA GenSA
#' @importFrom nloptr direct
#' @importFrom minqa newuoa
#' @importFrom smacof transPrep transform
#' @importFrom smacofx spp
#'
#'
#' @keywords clustering multivariate
#' @export
copstressMin <- function (delta, kappa=1, lambda=1, nu=1, theta=c(kappa,lambda,nu), type=c("ratio","interval","ordinal"), ties="primary", weightmat=1-diag(nrow(delta)), ndim = 2, init=NULL, stressweight=0.975,cordweight=0.025,q=1,minpts=ndim+1,epsilon=max(10,max(delta)),dmax=NULL,rang,optimmethod=c("NelderMead","Newuoa","BFGS","SANN","hjk","solnl","solnp","subplex","snomadr","hjk-Newuoa","hjk-BFGS","BFGS-hjk","Newuoa-hjk","cmaes","direct","direct-Newuoa","direct-BFGS","genoud","gensa"),verbose=0,scale=c("sd","rmsq","proc","none"),normed=TRUE, accuracy = 1e-7, itmax = 10000, stresstype=c("stress-1","stress"),principal=FALSE,...)
{
if(inherits(delta,"dist") || is.data.frame(delta)) delta <- as.matrix(delta)
if(!isSymmetric(delta)) stop("Delta is not symmetric.\n")
if(inherits(weightmat,"dist") || is.data.frame(weightmat)) weightmat <- as.matrix(weightmat)
if(!isSymmetric(weightmat)) stop("weightmat is not symmetric.\n")
## -- Setup for MDS type
if(missing(type)) type <- "ratio"
type <- match.arg(type, c("ratio", "interval", "ordinal"),several.ok = FALSE)
trans <- type
typo <- type
if (trans=="ratio"){
trans <- "none"
}
#TODO: if we want other optimal scalings as well
else if (trans=="ordinal" & ties=="primary"){
trans <- "ordinalp"
typo <- "ordinal (primary)"
} else if(trans=="ordinal" & ties=="secondary"){
trans <- "ordinals"
typo <- "ordinal (secondary)"
} else if(trans=="ordinal" & ties=="tertiary"){
trans <- "ordinalt"
typo <- "ordinal (tertiary)"
#} else if(trans=="spline"){
# trans <- "mspline"
}
if(type =="ordinal") theta <- c(kappa,1,nu) #We dont allow powers for dissimilarities in nonmetric MDS
kappa <- theta[1]
lambda <- theta[2]
nu <- theta[3]
#plot <- FALSE
n <- nrow(delta)
if (ndim > (n - 1)) stop("Maximum number of dimensions is n-1!")
if(is.null(rownames(delta))) rownames(delta) <- 1:n
labos <- rownames(delta)
if(verbose>0) cat(paste("Minimizing",type,"copstress with kappa=",kappa,"lambda=",lambda,"nu=",nu,".\n"))
if(missing(optimmethod)) optimmethod <- "hjk-Newuoa"
if(missing(scale)) scale <- "sd"
if(missing(stresstype)) stresstype <- "stress-1"
##-- Prepare for dissimilarity scaling
r <- kappa/2
deltaorig <- delta
delta <- delta^lambda
weightmato <- weightmat
weightmat <- weightmat^nu
weightmat[!is.finite(weightmat)] <- 0 #new
deltaold <- delta
disobj <- smacof::transPrep(as.dist(delta), trans = trans, spline.intKnots = 2, spline.degree = 2)#spline.intKnots = spline.intKnots, spline.degree = spline.degree) #FIXME: only works with dist() style object
## Add an intercept to the spline base transformation
#if (trans == "mspline") disobj$base <- cbind(rep(1, nrow(disobj$base)), disobj$base)
delta <- delta / enorm (delta, weightmat) #normalize to sum=1
## --- starting rang if not given
if(missing(rang))
#perhaps put this into the optimization function?
{
if(is.null(dmax))
{
if(is.null(init))
{
if(verbose>1) cat ("Fitting configuration for finding rang argument. \n")
initsol <- smacof::torgerson(delta,p=ndim)
init0 <- initsol
} else init0 <- init
init0 <- init0/enorm(init0)
# if(scale=="std") init0 <- scale(init0) #standardizes config before cordillera
if(scale=="none") init0 <- init0
if(scale=="sd") #scales config to sd=1 for most spread dimension before cordillera
{
init0 <- init0/max(apply(init0,2,stats::sd))
}
if(scale=="rmsq") #scales config to rmsq=1 for most spread dimension before cordillera
{
testso <- base::scale(init0,center=FALSE)
init0 <- init0/max(attr(testso,"scaled:scale"))
}
if(scale=="proc") #scales config by procrusting to init
{
if(is.null(init)) init <- initsol
procr <- smacof::Procrustes(init,init0)
init0 <- procr$Yhat
}
crp <- cordillera::cordillera(init0,q=q,minpts=minpts,epsilon=epsilon,scale=FALSE)$reachplot
cin <- max(crp)
dmax <- ifelse(kappa>1,1.5*cin,1.1*cin)
}
rang <- c(0,dmax) #This sets the range to (0,dmax)
if(verbose>1) cat("rang is",rang,".\n")
}
if(isTRUE(max(rang)!=dmax)) warning("Note: The supplied dmax and rang do not match. I took supplied rang as rang.\n")
## --- starting values
if(is.null(init))
{
if(exists("init0")) init <- init0 else init <- smacof::torgerson(delta,p=ndim)
#was: cops::powerStressFast(delta,kappa=kappa,lambda=lambda,nu=nu,ndim=ndim)$conf
}
xstart <- xold <- init
xold <- xold/enorm(xold)
dhat2 <- NA
#labs <- row.names(delta)
copsf <- function(x,delta,disobj,r,n,ndim,weightmat,stressweight,cordweight,q,minpts,epsilon,rang,scale,normed,init,...)
{
#init is used here only for Procrustes
if(!is.matrix(x)) x <- matrix(x,ncol=ndim)
delta <- delta/enorm(delta,weightmat)
x <- x/enorm(x)
dnew <- sqdist (x)
e <- as.dist(mkPower(dnew,r)) #was bug prior to 1.12-1
dhat2 <<- smacof::transform(e, disobj, w = as.dist(weightmat), normq = 0.5) ##I use <<- to change this also in the parent environment because the copsf function should only return a scalar for the optimizers but I the last dhats2 and also delta later on
dhatt <- dhat2$res
dhatd <- structure(dhatt, Size = n, call = quote(as.dist.default(m=b)), class = "dist", Diag = FALSE, Upper = FALSE)
delta <<- as.matrix(dhatd) ##I use <<- to save this also in the parent environment
rnew <- sum (weightmat * delta * mkPower (dnew, r))
nnew <- sum (weightmat * mkPower (dnew, 2*r))
anew <- rnew / nnew
stressi <- 1 - 2 * anew * rnew + (anew ^ 2) * nnew
if(stresstype=="stress-1") stressi <- sqrt(stressi)
if(scale=="none") x <- x
#if(scale=="std") x <- base::scale(x) #standardizes config before cordillera
if(scale=="sd") #scales config to sd=1 for most spread dimension before cordillera
{
x <- x/max(apply(x,2,stats::sd))
}
if(scale=="rmsq") #scales config to rmsq=1 for most spread dimension before cordillera
{
testso <- base::scale(x,center=FALSE)
x <- x/max(attr(testso,"scaled:scale"))
}
if(scale=="proc") #scales config by procrusting to init
{
procr <- smacof::Procrustes(init,x)
x <- procr$Yhat
}
corrd <- cordillera::cordillera(x,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=FALSE,...)
struc <- corrd$raw
if(normed) {
struc <- corrd$normed
}
ic <- stressweight*stressi - cordweight*struc
if(verbose>1) cat("copstress =",ic,"mdsloss =",stressi,"OC =",struc,"minpts=",minpts,"kappa =",kappa,"lambda =",lambda,"nu=",nu,"\n")
ic
}
if(verbose>1) cat("Starting Minimization with",optimmethod,":\n")
if(optimmethod=="Newuoa") {
suppressWarnings(optimized <- minqa::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmax,rhoend=accuracy,iprint=verbose-2),...))
itel <- itel+optimized$feval
ovalue <-optimized$fval
#optimized <- nloptr::newuoa(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmax,rhoend=accuracy,iprint=verbose-2),...)
#xnew <- matrix(optimized$par,ncol=ndim)
#itel <- optimized$iter
#ovalue <-optimized$value
}
if(optimmethod=="direct") {
xold <- as.vector(xold)
optimized <- nloptr::direct(function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),lower=rep(5*min(xold),length(xold)),upper=rep(5*max(xold),length(xold)),nl.info=isTRUE(verbose>1),control=list(maxeval=itmax,xtol_rel=accuracy),...)
# optimized <- nloptr::direct(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,#n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=),#lower=rep(5*min(xold),length(xold)),upper=rep(5*max(xold),length(xold)), nl.info=isTRUE(verbose>1),control=list(maxeval=itmax,xtol_rel=accuracy),...)
xnew <- matrix(optimized$par,ncol=ndim)
itel <- optimized$iter
ovalue <-optimized$value
}
if(optimmethod=="direct-Newuoa") {
xold <- as.vector(xold)
optimized1 <- nloptr::direct(function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),lower=rep(5*min(xold),length(xold)),upper=rep(5*max(xold),length(xold)),nl.info=isTRUE(verbose>1),control=list(maxeval=itmax,xtol_rel=accuracy),...)
xnew <- optimized1$par
itel1 <- optimized1$iter
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
suppressWarnings(optimized <- minqa::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmaxreduced,rhoend=accuracy,iprint=verbose-2),...))
itel <- itel1+optimized$feval
ovalue <-optimized$fval
#optimized <- nloptr::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),nl.info=isTRUE(verbose>2),control=list(maxeval=itmaxreduced,xtol_rel=accuracy),...)
#xnew <- matrix(optimized$par,ncol=ndim)
#itel <- itel1+optimized$iter
#ovalue <-optimized$value
}
if(optimmethod=="direct-BFGS") {
xold <- as.vector(xold)
optimized1 <- nloptr::direct(function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),lower=rep(5*min(xold),length(xold)),upper=rep(5*max(xold),length(xold)),nl.info=isTRUE(verbose>1),control=list(maxeval=itmax,xtol_rel=accuracy),...)
xnew <- matrix(optimized1$par,ncol=ndim)
itel1 <- optimized1$iter
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
optimized <- optim(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),method="BFGS",control=list(maxit=itmaxreduced,trace=0,reltol=accuracy),...)
xnew <- optimized$par
itel <- itel1+optimized$counts[[1]]
ovalue <-optimized$val
}
if(optimmethod=="genoud") {
optimized <- rgenoud::genoud(function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),starting.values=as.vector(xold),nvars=length(xold),print.level=verbose,...)
xnew <- matrix(optimized$par,ncol=ndim)
itel <- optimized$generations
ovalue <-optimized$value
}
if(optimmethod=="gensa") {
optimized <- GenSA::GenSA(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),lower=rep(5*min(xold),length(xold)),upper=rep(5*max(xold),length(xold)),control=list(max.call=itmax,verbose=isTRUE(verbose>1),smooth=FALSE),...)
xnew <- matrix(optimized$par,ncol=ndim)
itel <- optimized$counts
ovalue <-optimized$value
}
if(optimmethod=="cmaes") {
xold <- as.vector(xold)
optimized <- cmaes::cma_es(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxit=itmax),...)
xnew <- matrix(optimized$par,ncol=ndim)
itel <- optimized$counts[1]
ovalue <-optimized$value
}
if(optimmethod=="cmaes-Newuoa") {
xold <- as.vector(xold)
optimized1 <- cmaes::cma_es(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxit=itmax),...)
xnew <- matrix(optimized1$par,ncol=ndim)
itel1 <- optimized1$counts[1]
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
suppressWarnings(optimized <- minqa::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmaxreduced,rhoend=accuracy,iprint=verbose-2),...))
xnew <- matrix(optimized$par,ncol=ndim)
itel <- itel1+optimized$feval
ovalue <-optimized$fval
#optimized <- nloptr::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),nl.info=isTRUE(verbose>2),control=list(maxeval=itmaxreduced,xtol_rel=accuracy),...)
#xnew <- matrix(optimized$par,ncol=ndim)
#itel <- itel1+optimized$iter
#ovalue <-optimized$value
}
if(optimmethod=="Newuoa-cmaes") {
suppressWarnings(optimized1 <- minqa::newuoa(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmax,rhoend=accuracy,iprint=verbose-2),...))
#optimized1 <- nloptr::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),nl.info=isTRUE(verbose>2),control=list(maxeval=itmax,xtol_rel=accuracy),...)
xnew <- optimized1$par
itel1 <- optimized1$feval
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
optimized <- cmaes::cma_es(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxit=itmaxreduced),...)
xnew <- matrix(optimized$par,ncol=ndim)
itel <- optimized$counts[1]+itel1
ovalue <-optimized$value
}
if(optimmethod=="NelderMead") {
optimized <- optim(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxit=itmax,trace=0,reltol=accuracy),...)
xnew <- optimized$par
itel <- optimized$counts[[1]]
ovalue <-optimized$val
}
if(optimmethod=="BFGS") {
optimized <- optim(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),method="BFGS",control=list(maxit=itmax,trace=0,reltol=accuracy),...)
xnew <- optimized$par
itel <- optimized$counts[[1]]
ovalue <-optimized$val
}
if(optimmethod=="SANN") {
optimized <- optim(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),method="SANN",control=list(maxit=itmax,trace=0,reltol=accuracy),...)
xnew <- optimized$par
itel <- optimized$counts[[1]]
ovalue <-optimized$val
}
if(optimmethod=="hjk") {
optimized <- dfoptim::hjk(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfeval=itmax),...)
xnew <- optimized$par
itel <- optimized$feval
ovalue <-optimized$value
}
if(optimmethod=="hjk-Newuoa") { #twostep1
#cat("Before HJK Iterations:", itmax,"\n")
optimized1 <- dfoptim::hjk(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfeval=itmax),...)
xnew <- optimized1$par
itel1 <- optimized1$feval
#cat("OpVal after HJK:",optimized1$value,"\n")
#cat("After HJK Iterations:", itel1,"\n")
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
#cat("Before Newuoa Iterations:", itmaxreduced,"\n")
suppressWarnings(optimized <- minqa::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmaxreduced,rhoend=accuracy,iprint=verbose-2),...))
xnew <- matrix(optimized$par,ncol=ndim)
itel <- itel1+optimized$feval
ovalue <-optimized$fval
#optimized <- nloptr::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),nl.info=isTRUE(verbose>2),control=list(maxeval=itmaxreduced,xtol_rel=accuracy),...)
##BUG: Was init= before in all Newuoa (no argument). Changed in 1.3-6
#xnew <- matrix(optimized$par,ncol=ndim) #
##test function for return value fr <- function(x) { sum((x - matrix(c(8,5,3,4,6,7),ncol=2))^2) } if returned is matrix(c(8,5,3,4,6,7),ncol=2) it is cool
#cat("Newuoa Iterations:", optimized$iter,"\n")
#cat("Newuoa Iterations:", optimized$feval,"\n")
#itel <- itel1+optimized$iter
#cat("After Newuoa Iterations:", itel,"\n")
#ovalue <-optimized$value
#cat("OpVal after all:",ovalue,"\n")
}
if(optimmethod=="hjk-BFGS") { #twostep2
optimized1 <- dfoptim::hjk(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfeval=itmax,trace=0),...)
xnew <- optimized1$par
itel1 <- optimized1$feval
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
optimized <- optim(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),method="BFGS",control=list(maxit=itmaxreduced,trace=0,reltol=accuracy),...)
xnew <- optimized$par
itel <- itel1+optimized$counts[[1]]
ovalue <-optimized$val
}
if(optimmethod=="BFGS-hjk") { #twostep6
optimized1 <- optim(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),method="BFGS",control=list(maxit=itmax,trace=0),...)
xnew <- optimized1$par
itel1 <- optimized1$counts[[1]]
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
optimized <- dfoptim::hjk(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfeval=itmaxreduced,trace=0,tol=accuracy),...)
xnew <- optimized$par
itel <- itel1+optimized$feval
ovalue <-optimized$val
}
if(optimmethod=="BFGS-Newuoa") { #twostep5
optimized1 <- optim(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),method="BFGS",control=list(maxfeval=itmax,trace=0),...)
itel1 <- optimized1$counts[[1]]
xnew <- optimized1$par
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
suppressWarnings(optimized <- minqa::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmaxreduced,rhoend=accuracy,iprint=verbose-2),...))
xnew <- matrix(optimized$par,ncol=ndim)
itel <- itel1+optimized$feval
ovalue <-optimized$fval
#optimized <- nloptr::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),nl.info=isTRUE(verbose>2),control=list(maxeval=itmaxreduced,xtol_rel=accuracy),...)
#xnew <- matrix(optimized$par,ncol=ndim) #
#itel <- itel1+optimized$iter
#ovalue <-optimized$value
}
if(optimmethod=="hjk-solnl") {#twostep3
optimized1 <- dfoptim::hjk(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfeval=itmax,trace=0),...)
xnew <- optimized1$par
itel1 <- optimized1$feval
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
optimized <- NlcOptim::solnl(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),maxnFun=itmaxreduced,tolFun=accuracy,...)
xnew <- optimized$par
itel <- itel1+optimized$counts[[1]]
ovalue <-optimized$fn
}
if(optimmethod=="hjk-subplex") {#twostep4
optimized1 <- dfoptim::hjk(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfeval=itmax,trace=0),...)
xnew <- optimized1$par
itel1 <- optimized1$feval
itmaxreduced <- itmax-itel1
if(itel1>itmax) itmaxreduced <- 0.1*itmax
optimized <- subplex::subplex(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxit=itmaxreduced,reltol=accuracy),...)
xnew <- matrix(optimized$par,ncol=ndim)
itel <- itel1+optimized$count
ovalue <-optimized$value
}
# if(optimmethod=="dfsane") {
# optimized <- BB::dfsane(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxit=itmax,trace=0,tol=accuracy),...)
# xnew <- optimized$par
# itel <- optimized$feval
# ovalue <-optimized$residual
# }
if(optimmethod=="solnl") {
optimized <- NlcOptim::solnl(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),maxnFun=itmax,tolFun=accuracy,...)
xnew <- optimized$par
itel <- optimized$counts[[1]]
ovalue <-optimized$fn
}
## if(optimmethod=="isres") {
## optimized <- nloptr::isres(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),maxeval=itmax,trace=verbose-2,xtol_rel=accuracy,...)
## xnew <- optimized$par
## itel <- optimized$iter
## ovalue <-optimized$value
## }
if(optimmethod=="solnp") {
optimized <- Rsolnp::solnp(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(outer.iter=itmax,trace=0,tol=accuracy),...)
xnew <- matrix(optimized$pars,ncol=ndim)
itel <- optimized$nfuneval
ovalue <-utils::tail(optimized$values,1)
}
if(optimmethod=="subplex") {
optimized <- subplex::subplex(xold,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxit=itmax,reltol=accuracy),...)
xnew <- matrix(optimized$par,ncol=ndim)
itel <- optimized$count
ovalue <-optimized$value
}
if(optimmethod=="snomadr") {
copsf2 <- function(x,params)
{
copsf(x,delta=params[[1]],disobj=params[[2]],r=params[[3]],n=params[[4]],ndim=params[[5]],weightmat=params[[6]],stressweight=params[[7]],cordweight=params[[8]],q=params[[9]],minpts=params[[10]],epsilon=params[[11]],rang=params[[12]],scale=params[[13]],normed=params[[14]],init=params[[15]])
}
params <- list(delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init)
optimized <- crs::snomadr(copsf2,params=params,,n=dim(xold)[1],x0=xold,print.output=isTRUE(verbose-2>0),...)
xnew <- optimized$solution
itel <- optimized$iterations
ovalue <-optimized$objective
}
# if(optimmethod=="snomadr-Newuoa") {
# optimized1 <- crs::snomadr(function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),n=dim(xold)[1],x0=xold,print.output=isTRUE(verbose-2>0),...)
# xnew <- optimized1$solution
# itel <- optimized1$iterations
# optimized <- minqa::newuoa(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),control=list(maxfun=itmax-itel,rhoend=accuracy,iprint=verbose-2),...)
# xnew <- matrix(optimized$par,ncol=ndim)
# itel <- itel+optimized$feval
# ovalue <-optimized$fval
# }
# if(optimmethod=="snomadr-BFGS") {
# optimized1 <- crs::snomadr(function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=w#eightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),n=dim(xold)[1],x0=xold,print.output=isTRUE(verbose-2>0),...)
# xnew <- optimized1$solution
# itel <- optimized1$iterations
# optimized <- optim(xnew,function(par) copsf(par,delta=delta,disobj=disobj,r=r,n=n,ndim=ndim,weightmat=weightmat,stressweight=stressweight,cordweight=cordweight,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=scale,normed=normed,init=init),method="BFGS",control=list(maxit=itmax-itel,trace=0,reltol=accuracy),...)
# xnew <- optimized$par
# itel <- itel+optimized$counts[[1]]
# ovalue <-optimized$val
# }
xnew <- xnew/enorm(xnew)
rownames(delta) <- labos
dnew <- sqdist (xnew)
rnew <- sum (weightmat * delta * mkPower (dnew, r))
nnew <- sum (weightmat * mkPower (dnew, 2*r))
anew <- rnew / nnew
stress.m <- snew <- 1 - 2 * anew * rnew + (anew ^ 2) * nnew
stress <- sqrt(stress.m)
if(stresstype=="stress-1") {
stress.m <- sqrt(stress.m)
}
attr(xnew,"dimnames")[[1]] <- rownames(delta)
attr(xnew,"dimnames")[[2]] <- paste("D",1:ndim,sep="")
#doutm <- (2*sqrt(sqdist(xnew)))^kappa #fitted powered euclidean distance but times two
#doutm <- as.matrix(dist(xnew)^kappa)
doutm <- mkPower(sqdist(xnew),r)
#deltam <- delta
#deltaorigm <- deltaorig
#deltaoldm <- deltaold
delta <- stats::as.dist(delta)
deltaorig <- stats::as.dist(deltaorig)
deltaold <- stats::as.dist(deltaold)
#doute <- doutm/enorm(doutm)
#doute <- stats::as.dist(doute)
dout <- stats::as.dist(doutm)
#resmat <- as.matrix((delta - doute)^2)
#spp <- colMeans(resmat)
#weightmatm <-weightmat
weightmat <- stats::as.dist(weightmat)
weightmato <- stats::as.dist(weightmato)
spoint <- smacofx::spp(delta, dout, weightmat)
resmat<-spoint$resmat
rss <- sum(spoint$resmat[lower.tri(spoint$resmat)])
spp <- spoint$spp
#stressen <- sum(weightmat*(delta-doute)^2)#raw stress on the normalized proximities and normalized distances
#if(scale=="std") xnews <- base::scale(xnew) #standardizes config before cordillera
if(scale=="sd") #scales config to sd=1 for most spread dimension before cordillera
{
xnews <- xnew/max(apply(xnew,2,stats::sd))
}
if(scale=="rmsq") #scales config to rmsq=1 for most spread dimension before cordillera
{
testso <- base::scale(xnew,center=FALSE)
xnews <- xnew/max(attr(testso,"scaled:scale"))
}
if(scale=="proc") #scales config by procrusting to init
{
procr <- smacof::Procrustes(init,xnew)
xnews <- procr$Yhat
}
if(scale=="none") xnews <- xnew #no standardisation
if (principal) {
xnew_svd <- svd(xnew)
xnew <- xnew %*% xnew_svd$v
xnews_svd <- svd(xnews)
xnews <- xnews %*% xnews_svd$v
}
if(verbose>0) cat("*** Stress:",stress.m,"; Stress-1:",stress,"; from optimization: ",ovalue,"\n")
#like smacofP
out <- list(delta=deltaorig, tdelta=deltaold, dhat=delta, confdist=dout, iord=dhat2$iord.prim, conf = xnews, stress=stress, spp=spp, ndim=ndim, weightmat=weightmato, resmat=resmat, rss=rss, init=xstart, model="COPS-C", niter = itel, nobj = n, type = type, call=match.call(), stress.m=snew, alpha = anew, sigma = snew, pars=c(kappa=kappa,lambda=lambda,nu=nu),parameters=c(kappa=kappa,lambda=lambda,nu=nu),theta=c(kappa=kappa,lambda=lambda,nu=nu),tweightmat=weightmat)
#extra slots for class copsc
out$gamma <- NA
out$ties <- ties
out$usconf <- xnew
#out$loss <- "copstress"
out$OC <- cordillera::cordillera(out$conf,q=q,minpts=minpts,epsilon=epsilon,rang=rang,scale=FALSE)
out$OCv <- ifelse(normed,out$OC$normed,out$OC$raw)
out$copstress <- ovalue
out$optim <- optimized
out$stressweight <- stressweight
out$cordweight <- cordweight
out$optimethod <- optimmethod
#out$typo <- typo
#out$nobj <- dim(out$conf)[1]
class(out) <- c("copsc","cops","smacofP","smacofB","smacof")
out
}
#' @rdname copstressMin
#' @export
copsc <- copstressMin
#' @rdname copstressMin
#' @export
copStressMin <- copstressMin
Try the cops package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.