File: find.m.R

package info (click to toggle)
r-cran-blockmodeling 1.1.5-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 912 kB
  • sloc: ansic: 2,024; f90: 952; sh: 13; makefile: 5
file content (140 lines) | stat: -rw-r--r-- 7,255 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
#' @encoding UTF-8
#' @title Computing the threshold
#' 
#' @description
#' The functions compute the maximum  value of \code{m/cut} where a certain  block is still classified as \code{alt.blocks} and not "null".
#' The difference between \code{find.m} and \code{find.m2} it that \code{find.m} uses an optimization  approach and is faster and more precise
#' than \code{find.m2}. However, \code{find.m} only supports regular ("reg") and complete ("com") as \code{alt.blocks}, while \code{find.m2} supports all block types.
#' Also, \code{find.m} does not always work, especially if \code{cormet} is not "none".
#' 
# #' @usage find.m(M, clu, alt.blocks = "reg", diag = !is.list(clu),
# #' cormet = "none", half = TRUE, FUN = "max")
# #' find.m2(M, clu, alt.blocks = "reg", neval = 100, half = TRUE,
# #' ms = NULL, ...)
# #' find.cut(M, clu, alt.blocks = "reg", cuts = "all", ...)
#'
#' @param M A matrix representing the (usually valued) network. For now, only one-relational networks are supported.
#' The network can have one or more modes (different  kinds of units with no ties among themselves.
#' If the network is not two-mode, the matrix must be square.
#' @param clu A partition. Each unique value represents one cluster.
#' If the network  is one-mode, then this should be a vector, else a list of vectors, one for each mode.
#' @param alt.blocks Only one of allowed blocktypes, as alternative to the null block:\cr
#' "com" - complete block\cr
#' "rdo", "cdo" - row and column-dominant blocks (binary, valued, and implicit approach only)\cr
#' "reg" - (f-)regular block\cr
#' "rre", "cre" - row and column-(f-)regular blocks\cr
#' "rfn", "cfn" - row and column-dominant blocks (binary, valued, and implicit approach only)\cr
#' "den" - density block (binary approach only)\cr
#' "avg" - average block (valued approach only).
#' @param diag (default = \code{TRUE}) Should the special status  of diagonal be acknowledged.
#' @param cormet Which method should be used to correct for different maximum  error contributions\cr
#' "none" - no correction\cr
#' "censor" - censor values larger than \code{M}\cr
#' "correct" -  so that the maximum  possible error contribution of the cell is the same regardless of a condition (either that something  must be 0 or at least \code{M}).
#' @param FUN (default = "max") Function f used in row-f-regular, column-f-regular, and f-regular blocks.
#' @param cuts The cuts, which should be evaluated. If \code{cuts="all"} (default), all unique values are evaluated.
#' @param neval A number of different \code{m} values to be evaluated.
#' @param half Should the returned value of m be one half of the value where the inconsistencies are the same.
#' @param ms The values of m where the function should be evaluated.
#' @param \dots Other parameters to \code{critFunC}.
#'
#' @return A matrix of maximal \code{m/cut} values.
#'
#' @references
#' Doreian, P., Batagelj, V. & Ferligoj, A. \enc{Anuška}{Anuska} (2005). Generalized blockmodeling, (Structural analysis in the social sciences, 25). Cambridge [etc.]: Cambridge University Press.
#'
#' \enc{Žiberna, A.}{Ziberna, A.} (2007). Generalized Blockmodeling of Valued Networks. Social Networks, 29(1), 105-126. doi: 10.1016/j.socnet.2006.04.002
#' 
#' \enc{Žiberna, A.}{Ziberna, A.} (2008). Direct and indirect approaches to blockmodeling of valued networks in terms of regular equivalence. Journal of Mathematical Sociology, 32(1), 57-84. doi: 10.1080/00222500701790207
#' 
#' \enc{Žiberna, A.}{Ziberna, A.} (2014). Blockmodeling of multilevel networks. Social Networks, 39(1), 46-61. doi: 10.1016/j.socnet.2014.04.002
#'
#' @author \enc{Aleš Žiberna}{Ales Ziberna}
#' 
#' @seealso \code{\link{critFunC}} and maybe also \code{\link{optParC}}, \code{\link{plotMat}}
#' 
#' @keywords cluster
#' @importFrom stats optimize
#' @export

"find.m" <-
function(
	M,	#matrix of a network
	clu,	#partition
	alt.blocks="reg", #alternative block to null block (for now only "reg" is supported)
	diag=!is.list(clu) ,#allow diagonal blocks
	cormet="none", #should we correct for diferent maxismum error contributins
			# "censor" - censor values larger than m
			# "correct" -  so that the maxsimum possible error contribution of the cell is the same regardles of a condition (either that somthing must be o or at least m)
	half = TRUE,	# should the returned value of m be one half of the value where the incosnistencies are the same, otherwise, the m is restricted to max(M)
	FUN="max"
){
	mx<-max(M)*(1+ half)
	mn<-min(M)
  diag=diag
	if(is.list(clu)){
		k<-sapply(clu,function(x)length(unique(x)))
		clu<-lapply(clu,function(x)as.integer(factor(x)))
		if(length(k)>2) {
			for(i in 2:length(clu)){
				clu[[i]]<-clu[[i]] + max(clu[[i-1]])
	  		}
	  		clu<-unlist(clu)
	  		clu<-list(clu,clu)
	  	}
	} else {
		clu<-as.integer(factor(clu))
		clu<-list(clu,clu)
		k<-sapply(clu,function(x)length(unique(x)))
	}
	
	m<-matrix(NA,nrow=k[1],ncol=k[2])
	err<-list(
		reg=function(B,m,FUN){
			nr<-dim(B)[1]	#numer of rows
			nc<-dim(B)[2]	#numer of colums
			sr<-apply(B,1,FUN);er<-m-sr[sr<m]
			sc<-apply(B,2,FUN);ec<-m-sc[sc<m]
			return(sum(er)*nc+ sum(ec)*nr - sum(pmin(rep(er,times=length(ec)),rep(ec,each=length(er)))))	#regular block error
		},
		com=function(B,m,FUN){sumpos(m-B)}
	)
	
	errd<-list(
		com=function(B,m,FUN){sumpos(m-B) + min(0,sum(diag(B))-sumpos(m-diag(B)))},
		reg=function(B,m,FUN){
			nr<-dim(B)[1]	#numer of rows
			nc<-dim(B)[2]	#numer of colums
			sr<-apply(B,1,FUN);er<-m-sr[sr<m]
			sc<-apply(B,2,FUN);ec<-m-sc[sc<m]
			return(sum(er)*nc+ sum(ec)*nr - sum(pmin(rep(er,times=length(ec)),rep(ec,each=length(er))))) #regular block error
		}
	)
	errd.null<-function(B,m){sum(B) - min(0,sum(diag(B))-sumpos(m-diag(B)))}
	
	for(i in 1:k[1]){
		for(j in 1:k[2]){
			B<-M[clu[[1]]==i,clu[[2]]==j, drop=FALSE]
			if(ss(B)==0) m[i,j]<-min(2*B[1,1],mx) else if(i==j&&diag&&sum(dim(B))>1){
				if(errd.null(B,m=mx)>=errd[[alt.blocks]](B,mx,FUN)*ifelse(cormet=="correct",(mx - 0)/(mx - mn),1)){
					m[i,j]<-mx
				}else{ 
					m[i,j]<-optimize(f=function(m,B,alt.blocks,FUN,cormet,mx,mn){corf<-ifelse(cormet=="correct", (mx - 0)/(m - mn),1); if(cormet=="censor") B[B>m]<-m;(errd.null(B,m)-errd[[alt.blocks]](B,m,FUN)*corf)^2},lower=ifelse(cormet=="censor",mn,0),upper=mx,B=B,FUN=FUN,alt.blocks=alt.blocks,cormet=cormet,mx=mx,mn=mn)$minimum
					if(cormet=="correct" && errd.null(B)<err[[alt.blocks]](B,m[i,j],FUN)*(mx - 0)/(m[i,j] - mn)) m[m<=min(B)]<-0				
			}
				
			}else{
				if(sum(B)>=err[[alt.blocks]](B,mx,FUN)*ifelse(cormet=="correct",(mx - 0)/(mx - mn),1)){
					m[i,j]<-mx 
				}else{
					m[i,j]<-optimize(f=function(m,B,alt.blocks,FUN,cormet,mx,mn){corf<-ifelse(cormet=="correct", (mx - 0)/(m - mn),1); if(cormet=="censor") B[B>m]<-m;(sum(B)-err[[alt.blocks]](B,m,FUN)*corf)^2},lower=ifelse(cormet=="censor",mn,0),upper=mx,B=B,FUN=FUN,alt.blocks=alt.blocks,cormet=cormet,mx=mx,mn=mn)$minimum
					if(cormet=="correct" && sum(B)<err[[alt.blocks]](B,m[i,j],FUN)*(mx - 0)/(m[i,j] - mn)) m[m<=min(B)]<-0				
				}
			}
		}
	}
	if(cormet=="censor") m[m<min(M[M>0])]<-0
	if(half) m<-m/2
	return(m)
}