Nothing
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
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.