diff --git a/p3/data.go b/p3/data.go index d211f1d..2d9f037 100644 --- a/p3/data.go +++ b/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 +} diff --git a/p3/lrls/recursive.go b/p3/lrls/recursive.go index 314ea48..c20395d 100644 --- a/p3/lrls/recursive.go +++ b/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 } diff --git a/p3/run.go b/p3/run.go index 0d0c88e..2f82d77 100644 --- a/p3/run.go +++ b/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