# verfügbar unter:
# http://www.anicca-vijja.de/lg/R/skripte/optSchnitt.R

#########################
volldis <- function(distanz, both=FALSE)
{
# VOLLe DIStanzmatrix erzeugen aus Output von dist()
# mit Klassenattribut "class" = "dist"
 if(attr(distanz,"class") != "dist") stop("Keine mit dist() erzeugte Distanzmatrix")
 n <- attr(distanz,"Size")
 mat <- matrix(data=rep(0,n^2),nrow=n)
 mat[lower.tri(mat)] <- distanz
 if(both) mat <- mat + t(mat)
 return(mat)
}

#########################
# Protoypenberechnung = Maximum der Summe der Zeilen (Spalten)
# da 1 = connect, 0 = no connect
proto <- function(distanz)
 {
 return(apply(volldis(distanz, both=FALSE),1,sum))
 }

#########################
# Optimaler Schnitt durch eine Proximity-Matrix, hier Distanzmatrix
# cophcor() nach Oldenbürger (1981)
cophcor <- function(distanz, plotten = TRUE, dists = TRUE, verbose = FALSE)
{
# COPHenetic CORrelation = optimaler Schnitt durch eine Proximity-Matrix
# Quelle: Oldenbürger (1981)
# hier: Behandlung von mit dist() erzeugten Distanzmatrizen
# Es werden _keine_ allgemeinen Proximitymatrizen unterstützt
# braucht: volldis()
#
# (C) LG zuerst: 30.05.05, zuletzt: 24.10.05
# Lizenz: GPL http://www.gnu.org
#
# (1) Beispielaufruf zur Replierung, siehe set.seed():
#set.seed(1234)
#cophcor(dist(sample(1:100,10)))

 if(is.numeric(distanz) == FALSE) stop("Kein numerischer Inhalt")
 if(dists == TRUE) mat <- volldis(distanz, both=FALSE)

 n1 <- length(distanz)
 vek <- data.frame(distwert = sort(as.numeric(distanz)), corr = rep(0, n1),
                   ratio = rep(0,n1), sum.ne = rep(0,n1))
 NE <- numeric(n1)

 for(i in 1:(n1-1))
 {
  NE[i] <- 1
  vek$corr[i] <- cor(vek$distwert, NE)
  vek$sum.ne[i] <- sum(NE)
  vek$ratio[i] <- vek$sum.ne[i]/n1
  }

 if(dists) vek$corr <- abs(vek$corr)
 indize <- which(vek$corr == max(vek$corr))
 maxdw <- vek$distwert[indize]
 distanz[distanz <= maxdw] <- 1  # cutoff <= 1 ist ein progressives Kriterium
 distanz[distanz > maxdw] <- 0
 mat <- mat + t(volldis(distanz, both=FALSE))
 prot <- proto(distanz)
 maxi <- vek[indize,]
 rownames(maxi) <- "Maximum"
 if(plotten) plot(vek$distwert, vek$corr,
             main="Kophenetische Korrelation von\nDistanzmatrix und Null-Eins Matrix",
             ylim=c(0:1), xlab="Distanzwert", ylab="Korrelationskoeffizient (Pearson)")
 points(maxi$distwert, maxi$corr, pch=24, lty=2, lwd=1.2, col="red", bg="red")
 lines(vek$distwert, vek$corr, type="c", lty=4, col="blue")
 return(list("DistKorr"=vek,
             "MaxKorr"=cbind(indize,maxi),
             "NE"=mat,
             "ProtoV"=prot))
}