p3: implement Recursive Least Squares estimation

p3/data.go

@@ -4,7 +4,11 @@ import (
 	"encoding/csv"
 	"errors"
 	"io"
+	"io/fs"
+	"io/ioutil"
+	"log"
 	"os"
+	"regexp"
 	"strconv"
 )
 
@@ -80,12 +84,15 @@ func readFile(s *string) ([][]float64, error) {
 	return data, nil
 }
 
+// nolint: gocognit
 func saveStuff(
 	meanU, meanY,
 	varianceU, varianceY,
 	cov float64,
 	autocorrelationU, autocorrelationY,
-	mutCorrelationUY, mutCorrelationYU []float64,
+	mutCorrelationUY, mutCorrelationYU,
+	errVals []float64,
+	theta, thetaT [][]float64,
 ) error {
 	fFnames := map[string]float64{
 		"mean_u":        meanU,
@@ -99,6 +106,11 @@ func saveStuff(
 		"autocorrelation_y":     autocorrelationY,
 		"mutual_correlation_uy": mutCorrelationUY,
 		"mutual_correlation_yu": mutCorrelationYU,
+		"estimate_error":        errVals,
+	}
+	smFnames := map[string][][]float64{
+		"theta":  theta,
+		"thetaT": theta,
 	}
 	prefix := "data/"
 	suffix := ".txt"
@@ -137,6 +149,68 @@ func saveStuff(
 		}
 	}
 
+	for k, v := range smFnames {
+		if k == "thetaT" {
+			for i := range thetaT {
+				th := theta[i]
+
+				f, err := os.Create(prefix + k + "_p" + strconv.Itoa(i) + suffix)
+				if err != nil {
+					return err
+				}
+
+				defer f.Close()
+
+				w := csv.NewWriter(f)
+
+				s := recordsFromFloat64Slice(th)
+
+				err = w.WriteAll(s)
+				if err != nil {
+					return err
+				}
+			}
+		}
+
+		f, err := os.Create(prefix + k + suffix)
+		if err != nil {
+			return err
+		}
+
+		w := csv.NewWriter(f)
+
+		s := recordsFromSliceFloat64Slice(v)
+
+		err = w.WriteAll(s)
+		if err != nil {
+			return err
+		}
+
+		// since we need the file soon, close it now.
+		err = f.Close()
+		if err != nil {
+			return err
+		}
+
+		{
+			reread, err := ioutil.ReadFile(prefix + k + suffix)
+			if err != nil {
+				if errors.Is(err, fs.ErrNotExist) {
+					log.Println("haha")
+				}
+
+				return err
+			}
+
+			// remove double quotes.
+			cleared := regexp.MustCompile(`("?;"|";"?)|"`).ReplaceAll(reread, []byte("${1}"))
+
+			if err = ioutil.WriteFile(prefix+k+suffix, cleared, 0o600); err != nil {
+				return err
+			}
+		}
+	}
+
 	return nil
 }
 
@@ -151,3 +225,24 @@ func recordsFromFloat64Slice(f []float64) [][]string {
 
 	return s
 }
+
+func recordsFromSliceFloat64Slice(f [][]float64) [][]string {
+	s := make([][]string, 0, len(f))
+
+	for i := range f {
+		row := ""
+
+		for j := range f[i] {
+			tmp := strconv.FormatFloat(f[i][j], 'f', 64, 64)
+			if j == 0 {
+				row += tmp
+			} else {
+				row += ", " + tmp
+			}
+		}
+
+		s = append(s, []string{row})
+	}
+
+	return s
+}

p3/lrls/recursive.go

@@ -1,5 +1,183 @@
 package lrls
 
+import (
+	"gonum.org/v1/gonum/mat"
+)
+
 func recursive(u, y []float64) ([][]float64, []float64, error) {
-	return nil, nil, nil
+	dl := len(u)
+	if len(y) < dl {
+		dl = len(y)
+	}
+
+	// number of params.
+	p := 4
+	e := make([]float64, dl)
+
+	covMs, err := covMatrices(dl, p)
+	if err != nil {
+		return nil, nil, err
+	}
+
+	theta, err := newTheta(dl, p)
+	if err != nil {
+		return nil, nil, err
+	}
+
+	phi, err := newPhi(dl, p)
+	if err != nil {
+		return nil, nil, err
+	}
+
+	// the guts.
+	// in matlab:
+	// phi(:, k) = [-y(k-1) -y(k-2) u(k-1) u(k-2)]';
+	// e(k) = y(k) - phi(:, k)' * thm(:, k-1);
+	// C(:,:,k) = C(:,:,k-1) - (C(:,:,k-1)*phi(:,k)*phi(:,k)'*C(:,:,k-1))/(1+phi(:,k)'*C(:,:,k-1)*phi(:,k));
+	// thm(:,k) = thm(:,k-1) + C(:,:,k)*phi(:,k)*e(k);
+	for k := 2; k <= dl-1; k++ {
+		// phi[k] = []float64{-1.0 * y[k-1], -1 * y[k-2], -1 * u[k-1], -1 * u[k-2]}
+		// phi[k] = []float64{y[k-1], y[k-2], u[k-1], u[k-2]}
+		phi[k] = []float64{-1 * (y[k-1]), -1 * (y[k-2]), u[k-1], u[k-2]}
+		phiK := phi[k]
+		th := theta[k-1]
+
+		// phiT[k].theta[k-1].
+		phiThetaDotProd := phiThetaDotProd(p, phiK, th)
+
+		// get current error.
+		currentError := y[k] - phiThetaDotProd
+
+		// save current error.
+		e[k] = currentError
+
+		// next CM:
+		// cm(k) = cm(k-1) - [cm(k-1).phi(k).phiT(k).cm(k-1) / 1 + phiT(k).cm(k-1).phi(k)]
+		// next theta:
+		// theta(k) = theta(k-1) + cm(k).phi(k).e(k)
+
+		cmPrev := &mat.Dense{}
+		cmPrev.CloneFrom(covMs[k-1])
+
+		cmK := nextCM(p, cmPrev, phiK)
+
+		covMs[k] = cmK
+
+		tmp := &mat.Dense{}
+		phiVec := mat.NewVecDense(p, phiK)
+
+		tmp.Mul(cmK, phiVec)
+
+		// tmp.Scale(currentError, cmK) // -> this errs with dimens mismatch.
+		for i := range tmp.RawMatrix().Data {
+			tmp.RawMatrix().Data[i] *= currentError
+		}
+
+		tmpV := mat.NewVecDense(p, tmp.RawMatrix().Data)
+
+		thVec := mat.NewVecDense(p, th)
+		nuTh := mat.NewVecDense(p, nil)
+
+		nuTh.AddVec(thVec, tmpV)
+
+		theta[k] = nuTh.RawVector().Data
+	}
+
+	return theta, e, nil
+}
+
+func phiThetaDotProd(p int, phi, th []float64) float64 {
+	var res float64
+
+	for i := 0; i < p; i++ {
+		res += phi[i] * th[i]
+	}
+
+	return res
+}
+
+// nextCM calculates the value of the next Covariance Matrix.
+// cm(k) = cm(k-1) - [ cm(k-1).phi(k).phiT(k).cm(k-1) / 1 + phiT(k).cm(k-1).phi(k) ]
+// matlab:
+// C(k) = C(k-1) - (C(:,:,k-1)*phi(:,k)*phi(:,k)'*C(:,:,k-1))/(1+phi(:,k)'*C(:,:,k-1)*phi(:,k)).
+func nextCM(p int, cmPrev *mat.Dense, phi []float64) *mat.Dense {
+	phiVec := mat.NewVecDense(p, phi)
+	phiVecTransposed := phiVec.TVec()
+	cmTmp := &mat.Dense{}
+
+	// r, c := cmPrev.Dims()
+	// log.Printf("cmPrev dims: %d, %d", r, c)
+	// r, c = phiVec.Dims()
+	// log.Printf("phiVec dims: %d, %d", r, c)
+
+	cmTmp.Mul(phiVec, phiVecTransposed)
+
+	cmTmp.Mul(cmTmp, cmPrev)
+
+	cmTmp.Mul(cmPrev, cmTmp)
+
+	// -----
+
+	cmTmp2 := &mat.Dense{}
+
+	cmTmp2.Mul(phiVecTransposed, cmPrev)
+
+	vecTmp := mat.NewVecDense(p, cmTmp2.RawMatrix().Data)
+
+	denominator := 1 + doVector(
+		vecTmp.RawVector().Data,
+		phiVec.RawVector().Data,
+	)
+
+	for i := range cmTmp.RawMatrix().Data {
+		cmTmp.RawMatrix().Data[i] /= denominator
+	}
+
+	cm := &mat.Dense{}
+
+	cm.Sub(cmPrev, cmTmp)
+	cm = subMat(cmPrev, cmTmp)
+
+	return cm
+}
+
+// doVector returns a scalar product of column vector `cv` and regular vector
+// `rv`.
+func doVector(cv, rv []float64) float64 {
+	var res float64
+
+	for i := range cv {
+		res += cv[i] * rv[i]
+	}
+
+	return res
+}
+
+func subMat(m1, m2 *mat.Dense) *mat.Dense {
+	r, c := m1.Dims()
+	m := mat.NewDense(r, c, nil)
+
+	for i := range m1.RawMatrix().Data {
+		m.RawMatrix().Data[i] = m1.RawMatrix().Data[i] - m2.RawMatrix().Data[i]
+	}
+
+	return m
+}
+
+func TransposeTheta(theta [][]float64) [][]float64 {
+	thLen := len(theta)
+	elemLen := len(theta[0])
+	nuTheta := make([][]float64, elemLen)
+
+	for i := range theta[0] {
+		p := make([]float64, thLen)
+
+		for j := range theta {
+			p[j] = theta[j][i]
+		}
+
+		nuTheta[i] = p
+	}
+
+	return nuTheta
 }

p3/run.go

@@ -6,6 +6,7 @@ import (
 	"os/exec"
 
 	"git.dotya.ml/wanderer/ak9im/p2/stats"
+	"git.dotya.ml/wanderer/ak9im/p3/lrls"
 )
 
 var (
@@ -46,6 +47,13 @@ func run() error {
 
 		cov := stats.Covariance(data[u], data[y])
 
+		theta, errVals, err := lrls.Estimate(1, data[u], data[y])
+		if err != nil {
+			return err
+		}
+
+		thetaT := lrls.TransposeTheta(theta)
+
 		log.Printf("len(data): %d", len(data[u]))
 
 		log.Printf("means - u: %v, y: %v", meanU, meanY)
@@ -66,6 +74,8 @@ func run() error {
 			cov,
 			autocorrelationU, autocorrelationY,
 			mutCorrelationUY, mutCorrelationYU,
+			errVals,
+			theta, thetaT,
 		)
 		if err != nil {
 			return err