|
| 1 | +package golsv |
| 2 | + |
| 3 | +import ( |
| 4 | + "log" |
| 5 | + "math" |
| 6 | + "time" |
| 7 | +) |
| 8 | + |
| 9 | +// xxx optimization 8/1 |
| 10 | +// |
| 11 | +// - precompute matrix of row operations, store in sparse form, load and convert to dense. |
| 12 | +// - pass in Dense matrix to aligner |
| 13 | +// |
| 14 | +// - use optimized Dense * Sparse multiplication |
| 15 | +// |
| 16 | +// - (maybe) add special SubMultiplyRight to Dense, works same as |
| 17 | +// the optimized MultiplyRight but only uses part of each column |
| 18 | +// slice. |
| 19 | +// |
| 20 | +// - then in the aligner use SubMultiplyRight to obtain the tail |
| 21 | +// of vector w, and only continue to compute the rest of w if |
| 22 | +// that tail is nonzero (meaning the vector is not in the image |
| 23 | +// the cumulative B to that point). |
| 24 | +// |
| 25 | + |
| 26 | +type Aligner struct { |
| 27 | + B1smith *Sparse |
| 28 | + P *DenseBinaryMatrix |
| 29 | + B1colops []Operation |
| 30 | + Z1 *Sparse |
| 31 | + U BinaryMatrix // output matrix |
| 32 | + minWeight int |
| 33 | +} |
| 34 | + |
| 35 | +func NewAligner(B1smith *Sparse, P *DenseBinaryMatrix, B1colops []Operation, |
| 36 | + Z1 *Sparse) *Aligner { |
| 37 | + return &Aligner{ |
| 38 | + B1smith: B1smith, |
| 39 | + P: P, |
| 40 | + B1colops: B1colops, |
| 41 | + Z1: Z1, |
| 42 | + U: NewSparseBinaryMatrix(B1smith.NumRows(), 0), |
| 43 | + minWeight: math.MaxInt, |
| 44 | + } |
| 45 | +} |
| 46 | + |
| 47 | +// Align computes a matrix U such that the columns of U, together with |
| 48 | +// the columns of Bsmith, form a basis for Z. |
| 49 | +func (a *Aligner) Align() BinaryMatrix { |
| 50 | + BIsSmith, Brank := a.B1smith.IsSmithNormalForm() |
| 51 | + if !BIsSmith { |
| 52 | + log.Fatal("matrix B1smith is not in Smith normal form") |
| 53 | + } |
| 54 | + log.Printf("B1smith = %v, Brank = %d", a.B1smith, Brank) |
| 55 | + |
| 56 | + // nb. unexpectedly B1colops keeps coming out empty. xxx why? for |
| 57 | + // now assume it is true. |
| 58 | + if len(a.B1colops) > 0 { |
| 59 | + panic("B1colops not empty") |
| 60 | + } |
| 61 | + // Algorithm: |
| 62 | + // |
| 63 | + // let B = B1smith |
| 64 | + // let Z = Z1 |
| 65 | + // let d = rank of B |
| 66 | + // |
| 67 | + // for each column vector v of Z: |
| 68 | + // - compute w = P(v) (P is the row operations matrix) |
| 69 | + // - let k = max index of support of w |
| 70 | + // - if k > d, then |
| 71 | + // - v is not in the span of B |
| 72 | + // - append column w to B |
| 73 | + // - reduce new column of B; |
| 74 | + // this involves row swapping to create the pivot in |
| 75 | + // position d and adding columns to clear column d |
| 76 | + // - apply the same row operations to P |
| 77 | + // - increment d |
| 78 | + // - append column v to U |
| 79 | + // - else |
| 80 | + // - discard v since it is in the span of B |
| 81 | + // |
| 82 | + // in our main example, this should produce 19 column vectors in U, |
| 83 | + // since we already know that dim H_1 = 19. |
| 84 | + |
| 85 | + log.Printf("beginning alignment") |
| 86 | + Z := a.Z1 |
| 87 | + rows := Z.NumRows() |
| 88 | + |
| 89 | + statInterval := 1 |
| 90 | + timeStart := time.Now() |
| 91 | + timeLastReport := timeStart |
| 92 | + doneLastReport := 0 |
| 93 | + |
| 94 | + for j := 0; j < Z.NumColumns(); j++ { |
| 95 | + v := Z.Submatrix(0, rows, j, j+1) |
| 96 | + w := a.P.MultiplyRight(v) |
| 97 | + k := w.(*DenseBinaryMatrix).MaxColumnSupport(0) |
| 98 | + if k >= Brank { |
| 99 | + a.handleIndependentVector(v, w, Brank, k) |
| 100 | + Brank++ |
| 101 | + } |
| 102 | + doneLastReport++ |
| 103 | + if doneLastReport >= statInterval { |
| 104 | + now := time.Now() |
| 105 | + reportElapsed := now.Sub(timeLastReport) |
| 106 | + totalElapsed := now.Sub(timeStart) |
| 107 | + cRate := float64(doneLastReport) / reportElapsed.Seconds() |
| 108 | + tRate := float64(j) / totalElapsed.Seconds() |
| 109 | + log.Printf("align: processed %d/%d (%.2f%%) cols; crate=%1.0f trate=%1.0f found=%d", |
| 110 | + j, Z.NumColumns(), 100.0 * float64(j) / float64(Z.NumColumns()), |
| 111 | + cRate, tRate, a.U.NumColumns()) |
| 112 | + timeLastReport = now |
| 113 | + doneLastReport = 0 |
| 114 | + if reportElapsed.Seconds() < 1 { |
| 115 | + statInterval = statInterval * 2 |
| 116 | + } else if reportElapsed.Seconds() > 10 { |
| 117 | + statInterval = 1 + statInterval / 2 |
| 118 | + } |
| 119 | + } |
| 120 | + } |
| 121 | + log.Printf("done; min weight is %d", a.minWeight) |
| 122 | + return a.U |
| 123 | +} |
| 124 | + |
| 125 | +func (a *Aligner) handleIndependentVector(v BinaryMatrix, w BinaryMatrix, Brank int, k int) { |
| 126 | + a.U.(*Sparse).AppendColumn(v) |
| 127 | + log.Printf("found independent vector; found=%d total", a.U.NumColumns()) |
| 128 | + // sneak peak at systole |
| 129 | + weight := v.(*Sparse).ColumnWeight(0) |
| 130 | + if weight < a.minWeight { |
| 131 | + a.minWeight = weight |
| 132 | + log.Printf("new min weight %d", weight) |
| 133 | + } |
| 134 | + a.B1smith.AppendColumn(w) |
| 135 | + a.makePivot(Brank, k) |
| 136 | + log.Printf("clearing column") |
| 137 | + a.clearColumn(Brank) |
| 138 | + log.Printf("done clearing column") |
| 139 | +} |
| 140 | + |
| 141 | +func (a *Aligner) makePivot(Brank int, k int) { |
| 142 | + // we're increasing the rank of B by 1, so we want a pivot in |
| 143 | + // position Brank. |
| 144 | + if k > Brank { |
| 145 | + a.B1smith.SwapRows(Brank, k) |
| 146 | + a.P.SwapRows(Brank, k) |
| 147 | + } |
| 148 | +} |
| 149 | + |
| 150 | +// xxx needs test! |
| 151 | +func (a *Aligner) clearColumn(Brank int) { |
| 152 | + col := Brank |
| 153 | + row := Brank + 1 |
| 154 | + ops := 0 |
| 155 | + reportInterval := 1000 |
| 156 | + for { |
| 157 | + row = a.B1smith.ScanDown(row, col) |
| 158 | + if row == -1 { |
| 159 | + break |
| 160 | + } |
| 161 | + a.B1smith.Set(row, col, 0) |
| 162 | + a.P.AddRow(col, row) |
| 163 | + row++ |
| 164 | + ops++ |
| 165 | + if ops % reportInterval == 0 { |
| 166 | + log.Printf("did %d row ops", ops) |
| 167 | + } |
| 168 | + } |
| 169 | + log.Printf("did %d row ops", ops) |
| 170 | +} |
0 commit comments