R/profitLikeModel.R

In ProFit: Fit Projected 2D Profiles to Galaxy Images

profitLikeModel=function(parm, Data, makeplots=FALSE, 
  whichcomponents=list(sersic="all",moffat="all",ferrer="all",pointsource="all"), rough=FALSE,
  cmap = rev(colorRampPalette(brewer.pal(9,'RdYlBu'))(100)), errcmap=cmap, plotchisq=FALSE, maxsigma=5,
  model=NULL, image=Data$image, sigma=Data$sigma, region=Data$region, like.func=Data$like.func,
  algo.func=Data$algo.func, verbose=Data$verbose) {

  priorsum=0
  if(is.null(model))
  {
    if(class(Data)!='profit.data'){stop("The Data must be of class profit.data, as generated by the profitSetupData function!")}
    finesample = 1L
    if(length(Data$finesample)>0) finesample = Data$finesample
    profitCheckIsPositiveInteger(finesample)
    
    remakeout=profitRemakeModellist(parm=parm, Data=Data)
    modellistnew=remakeout$modellist
    parm=remakeout$parm
    
    # Calculate priors with the new versus old modellist
    if(length(Data$priors)>0){
      priorsum=Data$priors(modellistnew,Data$modellist)
    }
    
    # This is strictly the PSF for convolution with extended sources
    # It should evaluate as false if only point sources are being fit, because then there's
    # no need to generate a PSF image or verify the PSF dimensions
    psfdim = dim(Data$psf)
    if(Data$fitpsf) {
      psf = NULL
    } else {
      psf = Data$psf
    }
    
    openclenv=Data$openclenv
    if(identical(openclenv,new("externalptr"))) openclenv = NULL
    if(Data$usecalcregion){
      model = profitMakeModel(modellist=modellistnew, magzero = Data$magzero, psf=psf, dim=dim(image), psfdim=psfdim,
        whichcomponents = whichcomponents, rough=rough, calcregion=Data$calcregion, docalcregion=Data$usecalcregion,
        magmu=Data$magmu,finesample=finesample, convopt=Data$convopt, openclenv=openclenv, omp_threads=Data$omp_threads,
        adjust_calcregion = FALSE)
    }else{
      model = profitMakeModel(modellist=modellistnew, magzero = Data$magzero, psf=psf, dim=dim(image), psfdim=psfdim,
        whichcomponents = whichcomponents, rough=rough,
        magmu=Data$magmu, finesample=finesample, convopt=Data$convopt, openclenv=openclenv, omp_threads=Data$omp_threads,
        adjust_calcregion = FALSE)
    }
  } else {
    stopifnot(is.list(model) && !is.null(model$z))
    stopifnot(identical(dim(image),dim(model$z)))
    modellistnew=list()
  }

  cutsig=sigma
  cutim=image
  cutmod=model$z
  if(!is.null(region)) {
    cutim=cutim[region]
    cutmod=cutmod[region]
    cutsig=cutsig[region] 
  }

  #Force like.func to be lower case:
  like.func=profitParseLikefunc(like.func)
  
  #Various allowed likelihoods:
  isnorm = like.func == "norm"
  ischisq = like.func == "chisq"
  ist = like.func == "t"
  isst = like.func == "st"
  fitst = isst && is.null(Data$skewtparm)
  if(isnorm || ischisq || ist || isst) {
    cutsig=(cutim-cutmod)/cutsig
  }
  if("chisq" %in% Data$mon.names) chisq = sum(cutsig^2)
  if(ist || fitst)
  {
    vardata = var(cutsig,na.rm = TRUE)
    dof=2*vardata/(vardata-1)
    #dof=interval(dof,0,Inf)
    dof=max(1, min(Inf, dof, na.rm = TRUE), na.rm = TRUE)
  }
  skewtparm = Data$skewtparm
  if(isnorm){
    LL=sum(dnorm(cutsig, log=TRUE))
  } else if(ischisq) {
    ndata = length(cutim)
    if(!exists("chisq")) chisq = sum(cutsig^2)
    LL=dchisq(chisq, ndata, log=TRUE)
  } else if(ist) {
    vardata = var(cutsig,na.rm = TRUE)
    dof=2*vardata/(vardata-1)
    #dof=interval(dof,0,Inf)
    dof=max(1, min(Inf, dof, na.rm = TRUE), na.rm = TRUE)
    LL=sum(dt(cutsig,dof,log=TRUE))
  } else if(isst) {
    dstlike <- function(parm,Data) { 
      LP = sum(sn::dst(Data$chi, dp=parm, log=TRUE))
      return(list(LP=LP,Dev=-2*LP,Monitor=numeric(0),yhat=1,parm=parm))
    }
    if(is.null(Data$skewtparm)) {
      skewtparm = c(mean=median(cutsig), dof=dof, alpha=0.1, omega=1)
      STData = list(mon.names=character(0),parm.names=names(skewtparm), N=length(cutsig),chi=cutsig)
      stfit = LaplaceApproximation(Model=dstlike, parm=skewtparm, Data=STData, Iterations=1e3,
        Method="HAR",sir=FALSE)
      skewtparm = stfit$Summary1[,"Mode"]
      stfit = LaplacesDemon(Model=dstlike, Initial.Values=skewtparm, Data=STData, Iterations=1e3,
        Status = 100, Algorithm = "CHARM", Specs = list(alpha.star=0.44), Thinning = 1,
        CheckDataMatrixRanks = FALSE)
      skewtparm = stfit$Posterior1
    } else {
      skewtparm = Data$skewtparm
    }
    if(is.matrix(skewtparm)) {
      best = FALSE
      if(best) {
        LL = -Inf
        for(i in 1:dim(skewtparm)[1]) {
          LLi = dstlike(skewtparm[i,], Data=list(chi=cutsig))$LP
          if(LLi > LL) LL = LLi
        }
      } else {
        LL = 0
        for(i in 1:dim(skewtparm)[1]) LL = LL + dstlike(skewtparm[i,],
          Data=list(chi=cutsig))$LP
        LL = LL/dim(skewtparm)[1]
      }
    } else {
      LL=dstlike(skewtparm, Data=list(chi=cutsig))$LP
    }
  }
  else if(like.func=="pois") {
    scale=max(abs(image)/abs(sigma)^2)
    if(scale<0.1 | scale>10){
      cutmod=cutmod*scale
      cutim=cutim*scale
    }
    LL=-sum(cutmod-cutim*log(cutmod))
  } else {
    stop(paste0("Error: unknown likelihood function: '",like.func,"'"))
  }
  
  if(makeplots){
    skylevel = 0
    if(!is.null(modellistnew$sky) && !is.null(modellistnew$sky$bg) &&
      is.numeric(modellistnew$sky$bg))skylevel = modellistnew$sky$bg
    profitMakePlots(image-skylevel,model$z-skylevel,region, sigma, cmap=cmap,
      errcmap=errcmap,plotchisq=plotchisq,maxsigma=maxsigma, skewtparm=skewtparm)
  }
  
  LP=as.numeric(LL+priorsum)
  if(Data$verbose) {
    toprint = parm
    if(isTRUE(Data$printparmdiff)) toprint = parm-Data$init
    toprint = c(toprint,LP)
    if(isTRUE(Data$printLPdiff) && !is.null(Data$initLP) && is.numeric(Data$initLP))
    {
      toprint[length(toprint)] = toprint[length(toprint)] - Data$initLP
    }
    print(toprint,digits = 5)
  }
  if(algo.func=='') {
    out = list(model=model,psf=psf)
    if(fitst) out$skewtparm = skewtparm
  }
  if(algo.func=='optim' | algo.func=='CMA'){out=LP}
  if(algo.func=='LA' | algo.func=='LD')
  {
    Monitor=c(LL=LL,LP=LP)
    if("time" %in% Data$mon.names){
      Monitor = c(Monitor,tend = proc.time()["elapsed"])
    }
    if("chisq" %in% Data$mon.names){
      Monitor[which(Data$mon.names=="chisq")] = chisq
      names(Monitor)[which(Data$mon.names=="chisq")]='chisq'
    }
    if("dof" %in% Data$mon.names){
      Monitor[which(Data$mon.names=="dof")] = dof
      names(Monitor)[which(Data$mon.names=="dof")]='dof'
    }
    out=list(LP=LP,Dev=-2*LL,Monitor=Monitor,yhat=1,parm=parm)
  }
  return=out
}