[orientation] Use SumLog in EvaluateQ()

base.go

@@ -115,21 +115,29 @@ func Round(input float64) float64 {
 	return math.Floor(input + 0.5)
 }
 
+// SumLog returns the kernel of sum of log likelihood
+func SumLog(a []int) float64 {
+	sum := 0.0
+	for _, val := range a {
+		sum += math.Log(float64(val))
+	}
+	return sum
+}
+
 // HmeanInt returns the harmonic mean
 // That is:  n / (1/x1 + 1/x2 + ... + 1/xn)
-func HmeanInt(a []int, amin, amax int) int {
+func HmeanInt(a []int, amin, amax int) float64 {
 	size := len(a)
 	sum := 0.0
-	for i := 0; i < size; i++ {
-		val := a[i]
+	for _, val := range a {
 		if val > amax {
 			val = amax
 		} else if val < amin {
 			val = amin
 		}
 		sum += 1.0 / float64(val)
 	}
-	return int(Round(float64(size) / sum))
+	return float64(size) / sum
 }
 
 // GoldenArray is given list of ints, we aggregate similar values so that it becomes an

clm.go

@@ -69,7 +69,7 @@ type OrientedPair struct {
 type Contact struct {
 	strandedness int
 	nlinks       int
-	meanDist     int
+	meanDist     float64
 }
 
 // TigF stores the index to activeTigs and size of the tig
@@ -192,7 +192,8 @@ func (r *CLM) ParseClm() {
 
 		// Store all these info in contacts
 		gdists := GoldenArray(line.links)
-		meanDist := HmeanInt(line.links, GRLB, GRUB)
+		// meanDist := Hmean(line.links, GRLB, GRUB)
+		meanDist := SumLog(line.links)
 		strandedness := 1
 		if line.ao != line.bo {
 			strandedness = -1

orientation.go

@@ -11,6 +11,7 @@ package allhic
 
 import (
 	"fmt"
+	"math"
 
 	"github.com/gonum/matrix/mat64"
 )
@@ -155,9 +156,10 @@ func (r *CLM) Q() [][]GArray {
 // plus the actual link distances. Maximize Sum(1 / distance) for all links.
 // For performance consideration, we actually use a histogram to approximate
 // all link distances. See goldenArray() for details.
-func (r *CLM) EvaluateQ() (score float64) {
+func (r *CLM) EvaluateQ() float64 {
 	tour := r.Tour
 	Q := r.Q()
+	score := 0.0
 
 	size := tour.Len()
 	cumsize := make([]int, size)
@@ -181,10 +183,11 @@ func (r *CLM) EvaluateQ() (score float64) {
 				break
 			}
 			for k := 0; k < BB; k++ {
-				score += float64(Q[a][b][k]) / float64(GR[k]+dist)
+				// score += float64(Q[a][b][k]) / float64(GR[k]+dist)
+				score -= float64(Q[a][b][k]) * math.Log(float64(GR[k]+dist))
 			}
 			// fmt.Println(r.Tigs[a], r.Tigs[b], Q[a][b], score)
 		}
 	}
-	return
+	return score
 }