DaMAT.ttObject

   1# -*- coding: utf-8 -*-
   2"""
   3This is a python class for tensors in TT-format.
   4Through this object you can compute TT-decomposition of multidimensional arrays in
   5one shot using [TTSVD algorithm](https://epubs.siam.org/doi/epdf/10.1137/090752286)
   6or incrementally using [TT-ICE algorithm](https://arxiv.org/abs/2211.12487).
   7
   8Furthermore, this object allows exporting and importing TT-cores using native
   9format (`.ttc`) and intermediary formats (`.txt`).
  10"""
  11
  12import warnings
  13import time
  14
  15import numpy as np
  16
  17from logging import warning
  18from .utils import deltaSVD, ttsvd, coreContraction, mode_n_unfolding
  19from pickle import dump, load
  20
  21
  22class ttObject:
  23    """
  24    Python object for tensors in Tensor-Train format.
  25
  26    This object computes the TT-decomposition of multidimensional arrays using `TTSVD`_,
  27    `TT-ICE`_, and `TT-ICE*`_.
  28    It furthermore contains an inferior method, ITTD, for benchmarking purposes.
  29
  30    This object handles various operations in TT-format such as dot product and norm.
  31    Also provides supporting operations for uncompressed multidimensional arrays such as
  32    reshaping and transpose.
  33    Once the tensor train approximation for a multidimensional array is computed, you can
  34    compute projection of appropriate arrays onto the TT-cores, reconstruct projected
  35    or compressed tensors and compute projection/compression accuracy.
  36
  37    You can also save/load tensors as `.ttc` or `.txt` files.
  38
  39    Attributes
  40    ----------
  41    inputType: type
  42        Type of the input data. This determines how the object is initialized.
  43    originalData:
  44        Original multidimensional input array. Generally will not be stored
  45        after computing
  46        an initial set of TT-cores.
  47    method: :obj:`str`
  48        Method of computing the initial set of TT-cores. Currently only accepts
  49        `'ttsvd'` as input
  50    ttEpsilon: :obj:`float`
  51        Desired relative error upper bound.
  52    ttCores: :obj:`list` of :obj:`numpy.array`
  53        Cores of the TT-decomposition. Stored as a list of numpy arrays.
  54    ttRanks: :obj:`list` of :obj:`int`
  55        TT-ranks of the decomposition. Stored as a list of integers.
  56    nCores: :obj:`int`
  57        Number of TT-cores in the decomposition.
  58    nElements: :obj:`int`
  59        Number of entries present in the current `ttObject`.
  60    originalShape: :obj:`tuple` or :obj:`list`
  61        Original shape of the multidimensional array.
  62    reshapedShape: :obj:`tuple` or :obj:`list`
  63        Shape of the multidimensional array after reshaping. Note that
  64        this attribute is only meaningful before computing a TT-decomposition
  65    indexOrder: :obj:`list` of :obj:`int`
  66        Keeps track of original indices in case of transposing the input array.
  67    .. _TTSVD:
  68        https://epubs.siam.org/doi/epdf/10.1137/090752286
  69        _TT-ICE:
  70        https://arxiv.org/abs/2211.12487
  71        _TT-ICE*:
  72        https://arxiv.org/abs/2211.12487
  73
  74    """
  75
  76    def __init__(
  77        self,
  78        data,
  79        epsilon: float = None,
  80        keepData: bool = False,
  81        samplesAlongLastDimension: bool = True,
  82        method: str = "ttsvd",
  83    ) -> None:
  84        """
  85        Initializes the ttObject.
  86
  87        Parameters
  88        ----------
  89        data: :obj:`numpy.array` or :obj:`list`
  90            Main input to the ttObject. It can either be a multidimensional `numpy array`
  91            or `list of numpy arrays`.
  92            If list of numpy arrays are presented as input, the object will interpret it
  93            as the TT-cores of an existing decomposition.
  94        epsilon: :obj:`float`, optional
  95            The relative error upper bound desired for approximation. Optional for cases
  96            when `data` has type `list`.
  97        keepData: :obj:`bool`, optional
  98            Optional boolean variable to determine if the original array will be kept
  99            after compression. Set to `False` by default.
 100        samplesAlongLastDimension: :obj:`bool`, optional
 101            Boolean variable to ensure if the samples are stacked along the last
 102            dimension. Assumed to be `True` since it is one of the assumptions in TT-ICE.
 103        method: :obj:`str`, optional
 104            Determines the computation method for tensor train decomposition of
 105            the multidimensional array presented as `data`. Set to `'ttsvd'` by default.
 106
 107            Currently the package only has support for ttsvd, additional support such as
 108            `ttcross` might be included in the future.
 109
 110        """
 111        # self.ttRanks=ranks
 112        self.ttCores = None
 113        self.nCores = None
 114        self.nElements = None
 115        self.inputType = type(data)
 116        self.method = method
 117        self.keepOriginal = keepData
 118        self.originalData = data
 119        self.samplesAlongLastDimension = samplesAlongLastDimension
 120        if self.inputType is np.memmap:
 121            self.inputType = np.ndarray
 122
 123        if self.inputType == np.ndarray:
 124            self.ttEpsilon = epsilon
 125            self.originalShape = list(data.shape)
 126            self.reshapedShape = self.originalShape
 127            self.indexOrder = [idx for idx in range(len(self.originalShape))]
 128        elif self.inputType == list:
 129            self.nCores = len(data)
 130            self.ttCores = data
 131            self.reshapedShape = [core.shape[1] for core in self.ttCores]
 132        else:
 133            raise TypeError("Unknown input type!")
 134
 135    @property
 136    def coreOccupancy(self) -> None:
 137        """
 138        :obj:`list` of :obj:`float`: A metric showing the *relative rank* of each
 139        TT-core. This metric is used for a heuristic enhancement tool in `TT-ICE*`
 140        algorithm
 141        """
 142        try:
 143            return [
 144                core.shape[-1] / np.prod(core.shape[:-1]) for core in self.ttCores[:-1]
 145            ]
 146        except ValueError:
 147            warnings.warn(
 148                "No TT cores exist, maybe forgot calling object.ttDecomp?", Warning
 149            )
 150            return None
 151
 152    @property
 153    def compressionRatio(self) -> float:
 154        """
 155        :obj:`float`: A metric showing how much compression with respect to the
 156        original multidimensional array is achieved.
 157        """
 158        originalNumEl = 1
 159        compressedNumEl = 0
 160        for core in self.ttCores:
 161            originalNumEl *= core.shape[1]
 162            compressedNumEl += np.prod(core.shape)
 163        return originalNumEl / compressedNumEl
 164
 165    def changeShape(self, newShape: tuple or list) -> None:
 166        """
 167        Function to change shape of input tensors and keeping track of the reshaping.
 168        Reshapes `originalData` and saves the final shape in `reshapedShape`
 169
 170        Note
 171        ----
 172        A simple `numpy.reshape` would be sufficient for this purpose but in order to keep
 173        track of the shape changes the `reshapedShape` attribute also needs to be updated
 174        accordingly.
 175
 176        Parameters
 177        ----------
 178        newShape:obj:`tuple` or `list`
 179            New shape of the tensor
 180
 181        Raises
 182        ------
 183        warning
 184            If an attempt is done to modify the shape after computing a TT-decomposition.
 185            This is important since it will cause incompatibilities in other functions
 186            regarding the shape of the uncompressed tensor.
 187
 188        """
 189        if self.ttCores is not None:
 190            warning(
 191                "Warning! You are reshaping the original data after computing a\
 192                TT-decomposition! We will proceed without reshaping self.originalData!!"
 193            )
 194            return None
 195        self.reshapedShape = newShape
 196        self.originalData = np.reshape(self.originalData, self.reshapedShape)
 197        self.reshapedShape = list(self.originalData.shape)
 198        # Changed reshapedShape to a list for the trick in ttICEstar
 199        if self.samplesAlongLastDimension:
 200            self.singleDataShape = self.reshapedShape[:-1]
 201            # This line assumes we keep the last index as the samples index and don't
 202            # interfere with its shape
 203
 204    def computeTranspose(self, newOrder: list) -> None:
 205        """
 206        Transposes the axes of `originalData`.
 207
 208        Similar to `changeShape`, a simple
 209        `numpy.transpose` would be sufficient for this purpose but in order to
 210        keep track of the transposition order `indexOrder` attribute also needs
 211        to be updated accordingly.
 212
 213        Parameters
 214        ----------
 215        newOrder:obj:`list`
 216            New order of the axes.
 217
 218        Raises
 219        ------
 220        ValueError
 221            When the number of transposition axes are not equal to the number of
 222            dimensions of `originalData`.
 223        """
 224        assert self.inputType == np.ndarray and self.ttCores is None
 225        if len(newOrder) != len(self.indexOrder):
 226            raise ValueError(
 227                "size of newOrder and self.indexOrder does not match. \
 228                    Maybe you forgot reshaping?"
 229            )
 230        self.indexOrder = [self.indexOrder[idx] for idx in newOrder]
 231        self.originalData = self.originalData.transpose(newOrder)
 232        self.reshapedShape = list(self.originalData.shape)
 233
 234    def indexMap(self) -> None:
 235        """
 236        Function to map the original indices to the reshaped indices.
 237        Currently not implemented.
 238        """
 239        self.A = 2
 240
 241    def primeReshaping(self) -> None:
 242        """
 243        Function to reshape the first d dimensions to prime factors.
 244        Currently not implemented
 245        """
 246        self.A = 2
 247
 248    def saveData(
 249        self, fileName: str, directory="./", justCores=True, outputType="ttc"
 250    ) -> None:
 251        """
 252        Writes the computed TT-cores to an external file.
 253
 254        Parameters
 255        ----------
 256        fileName:obj:`str`
 257
 258        directory:obj:`str`
 259            Location to save files with respect to the present working directory.
 260        justCores:obj:`bool`
 261            Boolean variable to determine if `originalData` will be discarded
 262            or not while saving.
 263        outputType:obj:`str`
 264            Type of the output file. `ttc` for pickled `ttObject`, `txt` for
 265            individual text files for each TT-core.
 266        """
 267        if justCores:
 268            if outputType == "ttc":
 269                with open(directory + fileName + ".ttc", "wb") as saveFile:
 270                    temp = ttObject(self.ttCores)
 271                    for attribute in vars(self):
 272                        if attribute != "originalData":
 273                            setattr(temp, attribute, eval(f"self.{attribute}"))
 274                    dump(temp, saveFile)
 275            elif outputType == "txt":
 276                for coreIdx, core in enumerate(self.ttCores):
 277                    np.savetxt(
 278                        directory + f"{fileName}_{coreIdx}.txt",
 279                        core.reshape(-1, core.shape[-1]),
 280                        header=f"{core.shape[0]} {core.shape[1]} {core.shape[2]}",
 281                        delimiter=" ",
 282                    )
 283            elif outputType == "npy":
 284                for coreIdx, core in enumerate(self.ttCores):
 285                    np.save(
 286                        directory + f"{fileName}_{coreIdx}.npy",
 287                        core,
 288                    )
 289            else:
 290                raise ValueError(f"Output type {outputType} is not supported!")
 291        else:
 292            if outputType == "txt":
 293                raise ValueError(
 294                    ".txt type outputs are only supported for justCores=True!!"
 295                )
 296            if self.method == "ttsvd":
 297                with open(directory + fileName + ".ttc", "wb") as saveFile:
 298                    dump(self, saveFile)
 299            else:
 300                raise ValueError("Unknown Method!")
 301
 302    @staticmethod
 303    def loadData(fileName: str, numCores=None) -> "ttObject":
 304        """
 305        Loads data from a `.ttc`, `.npy`, or `.txt` file
 306
 307
 308        Static method to load TT-cores into a ttObject object.
 309        Note
 310        ----
 311        If data is stored in {coreFile}_{coreIdx}.txt format,
 312        the input fileName should just be coreFile.txt
 313
 314        Parameters
 315        ----------
 316        fileName:obj:`str`
 317            Name of the file that will be loaded.
 318        numCores:obj:`int`
 319            Number of cores that the resulting `ttObject` will have
 320            (Only required when input data format is `.txt`)
 321        """
 322        fileExt = fileName.split(".")[-1]
 323        if fileExt == "ttc":
 324            with open(fileName, "rb") as f:
 325                dataSetObject = load(f)
 326            return dataSetObject
 327        elif fileExt == "txt":
 328            if numCores is None:
 329                raise ValueError("Number of cores are not defined!!")
 330            fileBody = fileName.split(".")[0]
 331            coreList = []
 332            for coreIdx in range(numCores):
 333                with open(f"{fileBody}_{coreIdx}.{fileExt}"):
 334                    coreShape = f.readline()[2:-1]
 335                    coreShape = [int(item) for item in coreShape.split(" ")]
 336                coreList.append(
 337                    np.loadtxt(f"{fileBody}_{coreIdx}.{fileExt}").reshape[coreShape]
 338                )
 339            return coreList
 340        elif fileExt == "npy":
 341            if numCores is None:
 342                raise ValueError("Number of cores are not defined!!")
 343            fileBody = fileName.split(".")[0]
 344            coreList = []
 345            for coreIdx in range(numCores):
 346                coreList.append(
 347                    np.load(f"{fileBody}_{coreIdx}.{fileExt}")
 348                )
 349            return ttObject(coreList)
 350        else:
 351            raise ValueError(f"{fileExt} files are not supported!")
 352
 353    @staticmethod
 354    def ttDot(tt1, tt2) -> float():
 355        """
 356        Computes the dot product of two tensors in TT format.
 357
 358        Parameters
 359        ----------
 360        tt1:obj:`ttObject`
 361            Tensor in TT-format
 362        tt2:obj:`ttObject`
 363            Tensor in TT-format
 364
 365        Returns
 366        -------
 367        prod:obj:`float`
 368            Dot product of `tt1` and `tt2`.
 369
 370        Raises
 371        ------
 372        AttributeError
 373            When either one of the input parameters are not `ttObject`s.
 374
 375        """
 376        if not isinstance(tt1, ttObject) or not isinstance(tt2, ttObject):
 377            if isinstance(tt1, list):
 378                tt1 = ttObject(tt1)
 379            if isinstance(tt2, list):
 380                tt2 = ttObject(tt2)
 381            if not isinstance(tt1, ttObject) or not isinstance(tt2, ttObject):
 382                raise AttributeError("One of the passed objects is not in TT-format!")
 383        prod = np.kron(tt1.ttCores[0][:, 0, :], tt2.ttCores[0][:, 0, :])
 384        for i1 in range(1, tt1.ttCores[0].shape[1]):
 385            prod += np.kron(tt1.ttCores[0][:, i1, :], tt2.ttCores[0][:, i1, :])
 386        for coreIdx in range(1, len(tt1.ttCores)):
 387            p = []
 388            for ik in range(tt1.ttCores[coreIdx].shape[1]):
 389                p.append(
 390                    prod
 391                    @ (
 392                        np.kron(
 393                            tt1.ttCores[coreIdx][:, ik, :], tt2.ttCores[coreIdx][:, ik, :]
 394                        )
 395                    )
 396                )
 397            prod = np.sum(p, axis=0)
 398        return prod.item()
 399
 400    @staticmethod
 401    def ttNorm(tt1) -> float:
 402        """
 403        Computes norm of a tensor in TT-format.
 404
 405        Parameters
 406        ----------
 407        tt1:obj:`ttObject`
 408            Tensor in TT-format
 409
 410        Returns
 411        -------
 412        norm:obj:`float`
 413            Norm of `tt1`.
 414
 415        Raises
 416        ------
 417        AttributeError
 418            When the input parameter is not a `ttObject`.
 419        """
 420        if not isinstance(tt1, ttObject):
 421            if isinstance(tt1, list):
 422                tt1 = ttObject(tt1)
 423            else:
 424                raise AttributeError("Passed object is not in TT-format!")
 425        try:
 426            norm = ttObject.ttDot(tt1, tt1)
 427            norm = np.sqrt(norm)
 428        except MemoryError:
 429            norm = np.linalg.norm(tt1.ttCores[-1])
 430        return norm
 431
 432    def projectTensor(self, newData: np.array, upTo=None) -> np.array:
 433        """
 434        Projects tensors onto basis spanned by TT-cores.
 435
 436        Given a tensor with appropriate dimensions, this function leverages
 437        the fact that TT-cores obtained through `TTSVD`_ and `TT-ICE`_ are
 438        column-orthonormal in the mode-2 unfolding.
 439
 440        Note
 441        ----
 442        This function will not yield correct results if the TT-cores are
 443        not comprised of orthogonal vectors.
 444
 445        Parameters
 446        ----------
 447        newData:obj:`np.aray`
 448            Tensor to be projected
 449        upTo:obj:`int`, optional
 450            Index that the projection will be terminated. If an integer is
 451            passed as this parameter, `newTensor` will be projected up to
 452            (not including) the core that has index `upTo`. Assumes 1-based
 453            indexing.
 454
 455         .. _TTSVD:
 456            https://epubs.siam.org/doi/epdf/10.1137/090752286
 457            _TT-ICE:
 458            https://arxiv.org/abs/2211.12487
 459            _TT-ICE*:
 460            https://arxiv.org/abs/2211.12487
 461        """
 462        for coreIdx, core in enumerate(self.ttCores):
 463            if (coreIdx == len(self.ttCores) - 1) or coreIdx == upTo:
 464                break
 465            newData = (
 466                core.reshape(np.prod(core.shape[:2]), -1).transpose()
 467            ) @ newData.reshape(
 468                self.ttRanks[coreIdx] * self.ttCores[coreIdx].shape[1], -1
 469            )
 470        return newData
 471
 472    def reconstruct(self, projectedData, upTo=None):
 473        """
 474        Reconstructs tensors using TT-cores.
 475
 476        Assumes that `projectedData` is a slice from the last
 477        TT-core.
 478        While reconstructing any projected tensor from  `projectTensor`,
 479        this function leverages the fact that TT-cores obtained through
 480        `TTSVD`_ and `TT-ICE`_ are column-orthonormal in the mode-2
 481        unfolding.
 482
 483        Note
 484        ----
 485        This function might not yield correct results for projected tensors
 486        if the TT-cores are not comprised of orthogonal vectors.
 487
 488        Parameters
 489        ----------
 490        projectedData:obj:`np.aray`
 491            A tensor slice (or alternatively an array)
 492        upTo:obj:`int`, optional
 493            Index that the reconstruction will be terminated. If an integer is
 494            passed as this parameter, `projectedData` will be projected up to
 495            (not including) the core that has index `upTo`. Assumes 1-based
 496            indexing.
 497        .. _TTSVD:
 498            https://epubs.siam.org/doi/epdf/10.1137/090752286
 499            _TT-ICE:
 500            https://arxiv.org/abs/2211.12487
 501
 502        """
 503        if upTo is None:
 504            upTo = len(self.ttCores) - 1  # Write the core index in 1-indexed form!!
 505        for core in self.ttCores[:upTo][::-1]:
 506            projectedData = np.tensordot(core, projectedData, axes=(-1, 0))
 507        return projectedData
 508
 509    def updateRanks(self) -> None:
 510        """
 511        Updates the ranks of the `ttObject` after incremental updates.
 512        """
 513        self.ttRanks = [1]
 514        for core in self.ttCores:
 515            self.ttRanks.append(core.shape[-1])
 516        return None
 517
 518    def computeRelError(self, newTensor: np.array, useExact=True) -> np.array:
 519        """
 520        Computes relative error by projecting data onto TT-cores.
 521
 522        This function computes the error induced by projecting `data` onto the TT-cores
 523        of `ttObject`. The last index of `newTensor` is assumed to be the number of
 524        individual observations.
 525
 526        Note
 527        ----
 528        - In order to compute the projection onto the TT-cores, the dimensions of `data`
 529        should match that of `ttObject`.
 530        - If a single observation will be passed as `newTensor`, an additional
 531        index/dimension should be introduced either through reshaping or [:,None]
 532
 533        Parameters
 534        ----------
 535        newTensor:obj:`np.array`
 536            Tensor for which the projection error is computed
 537        useExact:`bool`
 538            Boolean flag determining whether a projection error will be
 539            returned for the entire batch or each observation in the batch.
 540
 541        Returns
 542        -------
 543        relError:obj:`np.array`
 544            Array of relative errors
 545        """
 546        elementwiseNorm = np.linalg.norm(newTensor, axis=0)
 547        for _ in range(len(newTensor.shape) - 2):
 548            elementwiseNorm = np.linalg.norm(elementwiseNorm, axis=0)
 549        projectedData = self.projectTensor(newTensor)
 550        reconstructedData = self.reconstruct(projectedData).reshape(newTensor.shape)
 551        difference = newTensor - reconstructedData
 552        differenceNorm = np.linalg.norm(difference, axis=0)
 553        for _ in range(len(difference.shape) - 2):
 554            differenceNorm = np.linalg.norm(differenceNorm, axis=0)
 555        
 556        if useExact:
 557            relError = np.linalg.norm(differenceNorm)/np.linalg.norm(elementwiseNorm)
 558        else:
 559            relError = differenceNorm / elementwiseNorm
 560        return relError
 561
 562    def computeRecError(self, data: np.array, start=None, finish=None, useExact=True) -> None:
 563        """
 564        Function to compute relative error by reconstructing data from slices
 565        of TT-cores.
 566
 567        Parameters
 568        ----------
 569        data:obj:`np.array`
 570            Tensor for which the reconstruction error is computed
 571        start:`int`
 572            Index for the start of the batch in the compressed last core
 573        finish:`int`
 574            Index for the end of the batch in the compressed last core 
 575        useExact:`bool`
 576            Boolean flag determining whether a reconstruction error will be
 577            returned for the entire batch or each observation in the batch.
 578
 579        Returns
 580        -------
 581        recError:obj:`float` or `np.array`
 582            Reconstruction error of the data from cores. Returns an array
 583            of reconstruction errors for each observation in the batch if
 584            `useExact` is `False`.
 585        """
 586        if start is None:
 587            start = self.ttCores[-1].shape[1]-1
 588        if finish is None:
 589            finish = start+1
 590        rec=self.reconstruct(self.ttCores[-1][:,start:finish,:]).reshape(data.shape)
 591        elementwiseNorm=np.linalg.norm(data,axis=0)
 592        for idx in range(len(data.shape)-2):
 593            elementwiseNorm=np.linalg.norm(elementwiseNorm,axis=0)
 594        difference=data-rec
 595        differenceNorm=np.linalg.norm(difference,axis=0)
 596        for idx in range(len(difference.shape)-2):
 597            differenceNorm=np.linalg.norm(differenceNorm,axis=0)
 598
 599        if useExact:
 600            recError=np.linalg.norm(differenceNorm)/np.linalg.norm(elementwiseNorm)
 601        else:
 602            recError=differenceNorm/elementwiseNorm
 603
 604        return recError
 605
 606    def ttDecomp(self, norm=None, dtype=np.float32) -> "ttObject.ttCores":
 607        """
 608        Computes TT-decomposition of a multidimensional array using `TTSVD`_ algorithm.
 609
 610        Currently only supports `ttsvd` as method. In the future additional formats may
 611        be covered.
 612
 613        Parameters
 614        ----------
 615        norm:obj:`float`, optional
 616            Norm of the tensor to be compressed
 617        dtype:obj:`type`, optional
 618            Desired data type for the compression. Intended to allow lower precision
 619            if needed.
 620
 621        Raises
 622        ------
 623        ValueError
 624            When `method` is not one of the admissible methods.
 625
 626
 627        The following attributes are modified as a result of this function:
 628        -------
 629        - `ttObject.ttCores`
 630        - `ttObject.ttRanks`
 631        - `ttObject.compressionRatio`
 632
 633        .. _TTSVD:
 634            https://epubs.siam.org/doi/epdf/10.1137/090752286
 635        """
 636        if norm is None:
 637            norm = np.linalg.norm(self.originalData)
 638        if self.method == "ttsvd":
 639            startTime = time.time()
 640            self.ttRanks, self.ttCores = ttsvd(
 641                self.originalData, norm, self.ttEpsilon, dtype=dtype
 642            )
 643            self.compressionTime = time.time() - startTime
 644            self.nCores = len(self.ttCores)
 645            self.nElements = 0
 646            for cores in self.ttCores:
 647                self.nElements += np.prod(cores.shape)
 648            if not self.keepOriginal:
 649                self.originalData = None
 650            return None
 651        else:
 652            raise ValueError("Method unknown. Please select a valid method!")
 653
 654    def ttICE(self, newTensor, epsilon=None, tenNorm=None) -> None:
 655        """
 656        `TT-ICE`_ algorithmn without any heuristic upgrades.
 657
 658        Given a set of TT-cores, this function provides incremental updates
 659        to the TT-cores to approximate `newTensor` within a relative error
 660        defined in `epsilon`
 661
 662        Note
 663        ----
 664        This algorithm/function relies on the fact that TT-cores are columnwise
 665        orthonormal in the mode-2 unfolding.
 666
 667        Parameters
 668        ----------
 669        newTensor:obj:`np.array`
 670            New/streamed tensor that will be used to expand the orthonormal bases defined
 671            in TT-cores
 672        epsilon:obj:`float`, optional
 673            Relative error upper bound for approximating `newTensor` after incremental
 674            updates. If not defined, `ttObject.ttEpsilon` is used.
 675        tenNorm:obj:`float`, optional
 676            Norm of `newTensor`
 677
 678        Notes
 679        -------
 680        **The following attributes are modified as a result of this function:**
 681        - `ttObject.ttCores`
 682        - `ttObject.ttRanks`
 683        - `ttObject.compressionRatio`
 684        .. _TT-ICE:
 685            https://arxiv.org/abs/2211.12487
 686        """
 687        if tenNorm is None:
 688            tenNorm = np.linalg.norm(newTensor)
 689        if epsilon is None:
 690            epsilon = self.ttEpsilon
 691        newTensor = newTensor.reshape(list(self.reshapedShape[:-1]) + [-1])[None, :]
 692        newTensorSize = len(newTensor.shape) - 1
 693        newTensor = newTensor.reshape(self.reshapedShape[0], -1)
 694        Ui = self.ttCores[0].reshape(self.reshapedShape[0], -1)
 695        Ri = newTensor - Ui @ (Ui.T @ newTensor)
 696        for coreIdx in range(0, len(self.ttCores) - 2):
 697            URi, _, _ = deltaSVD(Ri, tenNorm, newTensorSize, epsilon)
 698            self.ttCores[coreIdx] = np.hstack((Ui, URi)).reshape(
 699                self.ttCores[coreIdx].shape[0], self.reshapedShape[coreIdx], -1
 700            )
 701            self.ttCores[coreIdx + 1] = np.concatenate(
 702                (
 703                    self.ttCores[coreIdx + 1],
 704                    np.zeros(
 705                        (
 706                            URi.shape[-1],
 707                            self.reshapedShape[coreIdx + 1],
 708                            self.ttRanks[coreIdx + 2],
 709                        )
 710                    ),
 711                ),
 712                axis=0,
 713            )
 714            Ui = self.ttCores[coreIdx].reshape(
 715                self.ttCores[coreIdx].shape[0] * self.reshapedShape[coreIdx], -1
 716            )
 717            newTensor = (Ui.T @ newTensor).reshape(
 718                np.prod(self.ttCores[coreIdx + 1].shape[:-1]), -1
 719            )
 720            Ui = self.ttCores[coreIdx + 1].reshape(
 721                self.ttCores[coreIdx].shape[-1] * self.reshapedShape[coreIdx + 1], -1
 722            )
 723            Ri = newTensor - Ui @ (Ui.T @ newTensor)
 724        coreIdx = len(self.ttCores) - 2
 725        URi, _, _ = deltaSVD(Ri, tenNorm, newTensorSize, epsilon)
 726        self.ttCores[coreIdx] = np.hstack((Ui, URi))
 727        self.ttCores[coreIdx + 1] = np.concatenate(
 728            (
 729                self.ttCores[coreIdx + 1],
 730                np.zeros(
 731                    (
 732                        URi.shape[-1],
 733                        self.ttCores[coreIdx + 1].shape[1],
 734                        self.ttRanks[coreIdx + 2],
 735                    )
 736                ),
 737            ),
 738            axis=0,
 739        )
 740        newTensor = self.ttCores[coreIdx].T @ newTensor
 741        self.ttCores[coreIdx] = self.ttCores[coreIdx].reshape(
 742            self.ttCores[coreIdx - 1].shape[-1], self.reshapedShape[coreIdx], -1
 743        )
 744        coreIdx += 1
 745        Ui = self.ttCores[coreIdx].reshape(self.ttCores[coreIdx].shape[0], -1)
 746        self.ttCores[coreIdx] = np.hstack((Ui, newTensor)).reshape(
 747            self.ttCores[coreIdx].shape[0], -1, 1
 748        )
 749        self.updateRanks()
 750        return None
 751
 752    def ttICEstar(
 753        self,
 754        newTensor: np.array,
 755        epsilon: float = None,
 756        tenNorm: float = None,
 757        elementwiseNorm: np.array = None,
 758        elementwiseEpsilon: np.array = None,
 759        heuristicsToUse: list = ["skip", "subselect", "occupancy"],
 760        occupancyThreshold: float = 0.8,
 761        simpleEpsilonUpdate: bool = False,
 762        surrogateThreshold=True,
 763        sampleEpsilon=False,
 764        samplePercent=0.1
 765    ) -> None:
 766        """
 767        `TT-ICE*`_ algorithmn with heuristic performance upgrades.
 768
 769        Given a set of TT-cores, this function provides incremental updates
 770        to the TT-cores to approximate `newTensor` within a relative error
 771        defined in `epsilon`.
 772
 773        Note
 774        ----
 775        This algorithm/function relies on the fact that TT-cores are columnwise
 776        orthonormal in the mode-2 unfolding.
 777
 778        Parameters
 779        ----------
 780        newTensor:obj:`np.array`
 781            New/streamed tensor that will be used to expand the orthonormal bases defined
 782            in TT-cores
 783        epsilon:obj:`float`, optional
 784            Relative error upper bound for approximating `newTensor` after incremental
 785            updates. If not defined, `ttObject.ttEpsilon` is used.
 786        tenNorm:obj:`float`, optional
 787            Norm of `newTensor`.
 788        elementwiseNorm:obj:`np.array`, optional
 789            Individual norms of the observations in `newTensor`.
 790        elementwiseEpsilon:obj:`np.array`, optional
 791            Individual relative projection errors of the observations in `newTensor`.
 792        heuristicsToUse:obj:`list`, optional
 793            List of heuristics to use while updating TT-cores. Currently only accepts
 794            `'skip'`, `'subselect'`, and `'occupancy'`.
 795        occupancyThreshold:obj:`float`, optional
 796            Threshold determining whether to skip updating a single core or not. Not used
 797            if `'occupancy'` is not in `heuristicsToUse`
 798        simpleEpsilonUpdate:obj:`bool`, optional
 799            Uses the simple epsilon update equation. *Warning*: this relies on the
 800            assumption that all observations in `newTensor` have similar norms.
 801
 802        Notes
 803        -------
 804        **The following attributes are modified as a result of this function:**
 805        - `ttObject.ttCores`
 806        - `ttObject.ttRanks`
 807        - `ttObject.compressionRatio`
 808        .. _TT-ICE*:
 809            https://arxiv.org/abs/2211.12487
 810        """
 811        if epsilon is None:
 812            epsilon = self.ttEpsilon
 813        if ("subselect" in heuristicsToUse) and (newTensor.shape[-1] == 1):
 814            warning(
 815                "The streamed tensor has only 1 observation in it. \
 816                    Subselect heuristic will not be useful!!"
 817            )
 818        newTensor = newTensor.reshape(list(self.reshapedShape[:-1]) + [-1])[None, :]
 819        updEpsilon = epsilon
 820        newTensorSize = len(newTensor.shape) - 1
 821
 822        if elementwiseEpsilon is None:
 823            elementwiseEpsilon = self.computeRelError(newTensor)
 824        if "skip" in heuristicsToUse:
 825            if surrogateThreshold:
 826                if sampleEpsilon:
 827                    sampleSize=int(newTensor.shape[-1]*samplePercent)
 828                    skipCriterion=np.mean(np.random.choice(elementwiseEpsilon,size=sampleSize,replace=False))
 829                else:
 830                    skipCriterion=np.mean(elementwiseEpsilon)
 831                if  skipCriterion<= epsilon:
 832                    newTensor = self.projectTensor(newTensor)
 833                    self.ttCores[-1] = np.hstack(
 834                        (self.ttCores[-1].reshape(self.ttRanks[-2], -1), newTensor)
 835                    ).reshape(self.ttRanks[-2], -1, 1)
 836                    return None
 837            else:
 838                errorNorm=elementwiseEpsilon*elementwiseNorm
 839                errorNorm=np.sqrt(np.sum(errorNorm**2))/np.linalg.norm(elementwiseNorm)
 840                if errorNorm <= epsilon:
 841                    newTensor = self.projectTensor(newTensor)
 842                    self.ttCores[-1] = np.hstack(
 843                        (self.ttCores[-1].reshape(self.ttRanks[-2], -1), newTensor)
 844                    ).reshape(self.ttRanks[-2], -1, 1)
 845                    return None
 846        if tenNorm is None and elementwiseNorm is None:
 847            tenNorm = np.linalg.norm(newTensor)
 848        elif tenNorm is None:
 849            tenNorm = np.linalg.norm(elementwiseNorm)
 850
 851        select = [True] * newTensor.shape[-1]
 852        discard = [False] * newTensor.shape[-1]
 853        if "subselect" in heuristicsToUse:
 854            select = elementwiseEpsilon > epsilon
 855            discard = elementwiseEpsilon <= epsilon
 856            if simpleEpsilonUpdate:
 857                updEpsilon = (
 858                    epsilon * newTensor.shape[-1]
 859                    - np.mean(elementwiseEpsilon[discard]) * discard.sum()
 860                ) / (select.sum())
 861            else:
 862                if elementwiseNorm is None:
 863                    elementwiseNorm = np.linalg.norm(newTensor, axis=0)
 864                    for _ in range(len(self.ttCores) - 1):
 865                        elementwiseNorm = np.linalg.norm(elementwiseNorm, axis=0)
 866                allowedError = (self.ttEpsilon * np.linalg.norm(elementwiseNorm)) ** 2
 867                discardedError = np.sum(
 868                    (elementwiseEpsilon[discard] * elementwiseNorm[discard]) ** 2
 869                )
 870                updEpsilon = np.sqrt(
 871                    (allowedError - discardedError)
 872                    / (np.linalg.norm(elementwiseNorm[select]) ** 2)
 873                )
 874        self.reshapedShape[-1] = np.array(
 875            select
 876        ).sum()  # a little trick for ease of coding
 877
 878        indexString = "["
 879        for _ in range(len(self.reshapedShape)):
 880            # this heuristic assumes that the last dimension is for observations
 881            indexString += ":,"
 882        selectString = indexString + "select]"
 883        selected = eval("newTensor" + selectString)
 884
 885        selected = selected.reshape(list(self.reshapedShape[:-1]) + [-1])[None, :]
 886        for coreIdx in range(0, len(self.ttCores) - 1):
 887            selected = selected.reshape(np.prod(self.ttCores[coreIdx].shape[:-1]), -1)
 888            if ("occupancy" in heuristicsToUse) and (
 889                self.coreOccupancy[coreIdx] >= occupancyThreshold
 890            ):
 891                pass
 892            else:
 893                Ui = self.ttCores[coreIdx].reshape(
 894                    self.ttCores[coreIdx].shape[0] * self.reshapedShape[coreIdx], -1
 895                )
 896                Ri = selected - Ui @ (Ui.T @ selected)
 897                if (elementwiseNorm is None) or ("subselect" not in heuristicsToUse):
 898                    URi, _, _ = deltaSVD(
 899                        Ri, np.linalg.norm(selected), newTensorSize, updEpsilon
 900                    )
 901                else:
 902                    URi, _, _ = deltaSVD(
 903                        Ri,
 904                        np.linalg.norm(elementwiseNorm[select]),
 905                        newTensorSize,
 906                        updEpsilon,
 907                    )
 908                self.ttCores[coreIdx] = np.hstack((Ui, URi))
 909                self.ttCores[coreIdx + 1] = np.concatenate(
 910                    (
 911                        self.ttCores[coreIdx + 1],
 912                        np.zeros(
 913                            (
 914                                URi.shape[-1],
 915                                self.ttCores[coreIdx + 1].shape[1],
 916                                self.ttRanks[coreIdx + 2],
 917                            )
 918                        ),
 919                    ),
 920                    axis=0,
 921                )
 922
 923            self.ttCores[coreIdx] = self.ttCores[coreIdx].reshape(
 924                np.prod(self.ttCores[coreIdx].shape[:-1]), -1
 925            )
 926            # #project onto existing core and reshape for next core
 927            selected = (self.ttCores[coreIdx].T @ selected).reshape(
 928                self.ttCores[coreIdx + 1].shape[0] * self.reshapedShape[coreIdx + 1], -1
 929            )
 930            # fold back the previous core
 931            self.ttCores[coreIdx] = self.ttCores[coreIdx].reshape(
 932                -1, self.reshapedShape[coreIdx], self.ttCores[coreIdx].shape[-1]
 933            )
 934            # self.ttCores[coreIdx]=self.ttCores[coreIdx].reshape(self.ttCores[coreIdx].shape[0],self.reshapedShape[coreIdx],-1)
 935            # #fold back the previous core
 936        self.updateRanks()
 937        coreIdx += 1
 938        # coreIdx=len(self.ttCores), i.e working on the last core
 939        self.ttCores[coreIdx] = self.ttCores[coreIdx].reshape(
 940            self.ttCores[coreIdx].shape[0], -1
 941        )
 942        self.ttCores[coreIdx] = np.hstack(
 943            (self.ttCores[coreIdx], self.projectTensor(newTensor))
 944        ).reshape(self.ttCores[coreIdx].shape[0], -1, 1)
 945        return None
 946
 947    def ttRound(self, norm, epsilon=0) -> None:
 948        """
 949        Reorthogonalizes and recompresses TT-cores as per `Oseledets 2011`_.
 950
 951        This is the TT-rounding algorithm described in `TTSVD`_ paper. This
 952        function first reorthogonalizes the cores from last to first using
 953        QR decomposition and then recompresses the reorthogonalized TT-cores
 954        from first to last using `TTSVD`_ algorithm.
 955
 956        Parameters
 957        ----------
 958        norm:obj:`float`
 959            Norm of the TT-approximation. If there is no information on the
 960            uncompressed tensor, `ttNorm` function can be used.
 961        epsilon:obj:`float`, optional
 962            Relative error upper bound for the recompression step. Note that
 963            this relative error upper bound is not with respect to the original
 964            tensor but to the `compressed` tensor. Set to 0 by default.
 965
 966        Notes
 967        -----
 968        **The following attributes are modified as a result of this function:**
 969        - `ttObject.ttCores`
 970        - `ttObject.ttRanks`
 971        - `ttObject.compressionRatio`
 972
 973        .. _Oseledets 2011:
 974            https://epubs.siam.org/doi/epdf/10.1137/090752286
 975            TTSVD:
 976            https://epubs.siam.org/doi/epdf/10.1137/090752286
 977        """
 978        d = [core.shape[1] for core in self.ttCores]
 979        for coreIdx in np.arange(len(self.ttCores))[::-1]:
 980            currCore = self.ttCores[coreIdx]
 981            # Compute mode 1 unfolding of the tt-core
 982            currCore = currCore.reshape(currCore.shape[0], -1)
 983            # Using the transpose results in a row orthogonal Q matrix !!!
 984            q, r = np.linalg.qr(currCore.T)
 985            q, r = q.T, r.T
 986            self.ttCores[coreIdx] = q
 987            # Contract previous tt-core and R matrix
 988            self.ttCores[coreIdx - 1] = np.tensordot(
 989                self.ttCores[coreIdx - 1], r, axes=(-1, 0)
 990            )
 991            if coreIdx == 1:
 992                break
 993                # pass #TODO: write the case for the last core here.
 994        # Compression of the orthogonalized representation #
 995        ranks = [1]
 996        for coreIdx in range(len(self.ttCores) - 1):
 997            core = self.ttCores[coreIdx]
 998            core = core.reshape(np.prod(core.shape[:-1]), -1)
 999            self.ttCores[coreIdx], sigma, v = deltaSVD(
1000                core, norm, len(self.ttCores), epsilon
1001            )
1002            ranks.append(self.ttCores[coreIdx].shape[-1])
1003            # fold matrices into tt-cores
1004            self.ttCores[coreIdx] = self.ttCores[coreIdx].reshape(
1005                ranks[coreIdx], d[coreIdx], -1
1006            )
1007            self.ttCores[coreIdx + 1] = (
1008                (np.diag(sigma) @ v) @ self.ttCores[coreIdx + 1]
1009            ).reshape(ranks[-1] * d[coreIdx + 1], -1)
1010        self.ttCores[-1] = self.ttCores[-1].reshape(-1, d[-1], 1)
1011        self.ttRanks = [1]
1012        for core in self.ttCores:
1013            self.ttRanks.append(core.shape[-1])
1014        return None
1015
1016    @staticmethod
1017    def ittd(tt1, tt2, rounding=True, epsilon=0) -> "ttObject":
1018        """
1019        Alternative incremental tensor decomposition algorithm used to benchmark.
1020
1021        Note
1022        ----
1023        This method is empirically shown to be inferior to `TT-ICE`_ and
1024        `TT-ICE*`_.
1025
1026        Parameters
1027        ----------
1028        tt1:obj:`ttObject`
1029            Tensor in TT-format.
1030        tt2:obj:`ttObject`
1031            Tensor in TT-format.
1032        rounding:obj:`bool`, optional
1033            Determines if `ttRounding` will be done to the incremented TT-core.
1034            Note that for unbounded streams, rounding is essential. Otherwise
1035            TT-ranks will grow unboundedly.
1036        epsilon:obj:`float`, optional
1037            Relative error upper bound for the recompression step. Note that
1038            this relative error upper bound is not with respect to the original
1039            tensor but to the `compressed` tensor. Set to 0 by default. Not used
1040            if `rounding` is set to `False`.
1041
1042        Notes
1043        -----
1044        **The following attributes are modified as a result of this function:**
1045        - `ttObject.ttCores`
1046        - `ttObject.ttRanks`
1047        - `ttObject.compressionRatio`
1048
1049        .. _TT-ICE:
1050            https://arxiv.org/abs/2211.12487
1051            _TT-ICE*:
1052            https://arxiv.org/abs/2211.12487
1053        """
1054        if not isinstance(tt1, ttObject) or not isinstance(tt2, ttObject):
1055            if isinstance(tt1, list):
1056                tt1 = ttObject(tt1)
1057            if isinstance(tt2, list):
1058                tt2 = ttObject(tt2)
1059            if not isinstance(tt1, ttObject) or not isinstance(tt2, ttObject):
1060                raise AttributeError("One of the passed objects is not in TT-format!")
1061        tt1.ttCores[-1] = tt1.ttCores[-1].squeeze(-1)
1062        tt2.ttCores[-1] = tt2.ttCores[-1].squeeze(-1)
1063        rank1, numObs1 = tt1.ttCores[-1].shape
1064        rank2, numObs2 = tt2.ttCores[-1].shape
1065        tt1.ttCores[-1] = np.concatenate(
1066            (tt1.ttCores[-1], np.zeros((rank1, numObs2))), axis=1
1067        )
1068        tt2.ttCores[-1] = np.concatenate(
1069            (np.zeros((rank2, numObs1)), tt2.ttCores[-1]), axis=1
1070        )
1071        # Addition in TT-format #
1072        for coreIdx in range(len(tt1.ttCores)):
1073            if coreIdx == 0:
1074                tt1.ttCores[coreIdx] = np.concatenate(
1075                    (tt1.ttCores[coreIdx].squeeze(0), tt2.ttCores[coreIdx].squeeze(0)),
1076                    axis=1,
1077                )[None, :, :]
1078            elif coreIdx == len(tt1.ttCores) - 1:
1079                tt1.ttCores[coreIdx] = np.concatenate(
1080                    (tt1.ttCores[coreIdx], tt2.ttCores[coreIdx]), axis=0
1081                )[:, :, None]
1082            else:
1083                s11, s12, s13 = tt1.ttCores[coreIdx].shape
1084                s21, s22, s23 = tt2.ttCores[coreIdx].shape
1085                tt1.ttCores[coreIdx] = np.concatenate(
1086                    (tt1.ttCores[coreIdx], np.zeros((s11, s12, s23))), axis=2
1087                )
1088                tt2.ttCores[coreIdx] = np.concatenate(
1089                    (np.zeros((s21, s22, s13)), tt2.ttCores[coreIdx]), axis=2
1090                )
1091                tt1.ttCores[coreIdx] = np.concatenate(
1092                    (tt1.ttCores[coreIdx], tt2.ttCores[coreIdx]), axis=0
1093                )
1094        tt1.ttRanks = [1]
1095        for core in tt1.ttCores:
1096            tt1.ttRanks.append(core.shape[-1])
1097        if rounding:
1098            tenNorm = ttObject.ttNorm(tt1)
1099            tt1.ttRound(tenNorm, epsilon=epsilon)
1100        return tt1
1101
1102    @staticmethod
1103    def ttfoa6d(data,ttRanks,ttCores=None,S=None,forgettingFactor=None):
1104        nDims=len(data.shape)
1105        dims=list(data.shape)
1106        coreSwap=False
1107        if forgettingFactor is None:
1108            forgettingFactor=0.7
1109
1110        if len(ttRanks) != nDims-1:
1111            raise ValueError("Number of dimensions do not match with number of TT-ranks")
1112        
1113        if ttCores is None:
1114            # Initialize TT-cores randomly if no core is given.
1115            ttCores=[]
1116            ttCores.append(np.random.randn(1,dims[0],ttRanks[0]))
1117            for dimIdx in range(1,nDims-1):
1118                ttCores.append(np.random.randn(ttRanks[dimIdx-1],dims[dimIdx],ttRanks[dimIdx]))
1119            ttCores.append(np.empty((ttRanks[-1],1)))
1120            coreSwap=True
1121        # else:
1122        ttCores_old=ttCores.copy()
1123        if S is None:
1124            # Initialize S matrices as identity matrices if no S is given from previous step.
1125            S=[]
1126            S.append(100*np.eye(ttRanks[0]))
1127            for dimIdx in range(1,nDims-1):
1128                S.append(100*np.eye(ttRanks[dimIdx-1]*ttRanks[dimIdx]))
1129
1130        # Update the last core G6
1131        
1132        Ht_1=coreContraction(ttCores_old[:-1])
1133        # Ht_1=Ht_1.reshape(-1,ttRanks[-1])
1134        Ht_1=mode_n_unfolding(Ht_1,5).T #this might be 5
1135        g6=np.linalg.lstsq(Ht_1,data.reshape(-1,order='F'),rcond=None)[0].reshape(-1,1,order='F')
1136        if coreSwap:
1137            ttCores[-1]=g6
1138        else:
1139            ttCores[-1]=np.hstack((ttCores[-1],g6))
1140        recData=(Ht_1@g6).reshape(dims,order='F')
1141        delta=data-recData
1142
1143        # Update G1
1144        deltaunfold=mode_n_unfolding(delta,0)
1145        Bk=coreContraction(ttCores_old[1:-1]+[g6])
1146        Wk=mode_n_unfolding(Bk,0)
1147        S[0]=forgettingFactor*S[0]+Wk@Wk.T
1148        # Vk=np.linalg.lstsq(S[0],Wk,rcond=None)[0]
1149        Vk=np.linalg.pinv(S[0])@Wk
1150        ttCores[0]=ttCores_old[0]+deltaunfold@Vk.T
1151
1152        # Update G2
1153        deltaunfold=mode_n_unfolding(delta,1)
1154        Ak=mode_n_unfolding(ttCores_old[0],2)#.T
1155        Bk=coreContraction(ttCores_old[2:-1]+[g6])
1156        Bk=mode_n_unfolding(Bk,0)
1157        Wk=np.kron(Bk,Ak)
1158        S[1]=forgettingFactor*S[1]+Wk@Wk.T
1159        # Vk=np.linalg.lstsq(S[1],Wk,rcond=None)[0]
1160        Vk=np.linalg.pinv(S[1])@Wk
1161        ttCores[1]=mode_n_unfolding(ttCores_old[1],1)
1162        ttCores[1]=ttCores[1]+deltaunfold@Vk.T
1163        ttCores[1]=ttCores[1].reshape(dims[1],ttRanks[0],ttRanks[1],order='F').transpose(1,0,2)
1164
1165        # Update G3
1166        deltaunfold=mode_n_unfolding(delta,2)
1167        Ak=mode_n_unfolding(coreContraction(ttCores_old[:2]),2)
1168        Bk=coreContraction(ttCores_old[3:-1]+[g6]) #insert a mode 0 unfolding here if necessary
1169        Bk=mode_n_unfolding(Bk,0)
1170        Wk=np.kron(Bk,Ak)
1171        S[2]=forgettingFactor*S[2]+Wk@Wk.T
1172        # Vk=np.linalg.lstsq(S[2],Wk,rcond=None)[0]
1173        Vk=np.linalg.pinv(S[2])@Wk
1174        ttCores[2]=mode_n_unfolding(ttCores_old[2],1)
1175        ttCores[2]=ttCores[2]+deltaunfold@Vk.T
1176        ttCores[2]=ttCores[2].reshape(dims[2],ttRanks[1],ttRanks[2],order='F').transpose(1,0,2)
1177
1178        # Update G4
1179        deltaunfold=mode_n_unfolding(delta,3)
1180        Ak=mode_n_unfolding(coreContraction(ttCores_old[:3]),3)
1181        Bk=coreContraction(ttCores_old[4:-1]+[g6]) #insert a mode 0 unfolding here if necessary
1182        Wk=np.kron(Bk,Ak)
1183        S[3]=forgettingFactor=S[3]+Wk@Wk.T
1184        # Vk=np.linalg.lstsq(S[3],Wk,rcond=None)[0]
1185        Vk=np.linalg.pinv(S[3])@Wk
1186        ttCores[3]=mode_n_unfolding(ttCores_old[3],1)
1187        ttCores[3]=ttCores[3]+deltaunfold@Vk.T
1188        ttCores[3]=ttCores[3].reshape(dims[3],ttRanks[2],ttRanks[3],order='F').transpose(1,0,2)
1189
1190        # Update G5
1191        deltaunfold=mode_n_unfolding(delta,4)
1192        Ak=mode_n_unfolding(coreContraction(ttCores_old[:4]),4)
1193        # Bk=coreContraction(g6) #insert a mode 0 unfolding here if necessary
1194        Wk=np.kron(g6,Ak)
1195        S[4]=forgettingFactor=S[4]+Wk@Wk.T
1196        # Vk=np.linalg.lstsq(S[4],Wk,rcond=None)[0]
1197        Vk=np.linalg.pinv(S[4])@Wk
1198        ttCores[4]=mode_n_unfolding(ttCores_old[4],1)
1199        ttCores[4]=ttCores[4]+deltaunfold@Vk.T
1200        ttCores[4]=ttCores[4].reshape(dims[4],ttRanks[3],ttRanks[4],order='F').transpose(1,0,2)
1201
1202        return S,ttCores
1203
1204    @staticmethod
1205    def ttfoa4d(data,ttRanks,ttCores=None,S=None,forgettingFactor=None):
1206        nDims=len(data.shape)
1207        dims=list(data.shape)
1208        coreSwap=False
1209        if forgettingFactor is None:
1210            forgettingFactor=0.7
1211
1212        if len(ttRanks) != nDims-1:
1213            raise ValueError("Number of dimensions do not match with number of TT-ranks")
1214        
1215        if ttCores is None:
1216            # Initialize TT-cores randomly if no core is given.
1217            ttCores=[]
1218            ttCores.append(np.random.randn(1,dims[0],ttRanks[0]))
1219            for dimIdx in range(1,nDims-1):
1220                ttCores.append(np.random.randn(ttRanks[dimIdx-1],dims[dimIdx],ttRanks[dimIdx]))
1221            ttCores.append(np.empty((ttRanks[-1],1)))
1222            coreSwap=True
1223        # else:
1224        ttCores_old=ttCores.copy()
1225        if S is None:
1226            # Initialize S matrices as identity matrices if no S is given from previous step.
1227            S=[]
1228            S.append(100*np.eye(ttRanks[0]))
1229            for dimIdx in range(1,nDims-1):
1230                S.append(100*np.eye(ttRanks[dimIdx-1]*ttRanks[dimIdx]))
1231
1232        Ht_1=coreContraction(ttCores_old[:-1])
1233        # Ht_1=Ht_1.reshape(-1,ttRanks[-1])
1234        Ht_1=mode_n_unfolding(Ht_1,3).T #this might be 4
1235        g4=np.linalg.lstsq(Ht_1,data.reshape(-1,order='F'),rcond=None)[0].reshape(-1,1,order='F')
1236        if coreSwap:
1237            ttCores[-1]=g4
1238        else:
1239            ttCores[-1]=np.hstack((ttCores[-1],g4))
1240        recData=(Ht_1@g4).reshape(dims,order='F')
1241        delta=data-recData
1242
1243        # Update G1
1244        deltaunfold=mode_n_unfolding(delta,0)
1245        # Bk=coreContraction(ttCores_old[1:-1]+[g4])
1246        Bk=np.tensordot(np.tensordot(ttCores_old[1],ttCores_old[2],axes=[-1,0]),g4,axes=[-1,0])
1247        Wk=mode_n_unfolding(Bk,0)
1248        S[0]=forgettingFactor*S[0]+Wk@Wk.T
1249        # Vk=np.linalg.lstsq(S[0],Wk,rcond=None)[0]
1250        try:
1251            Vk=np.linalg.solve(S[0],Wk)
1252        except np.linalg.LinAlgError:
1253            print('Singular S[0] matrix, using Moore-Penrose inverse')
1254            Vk=np.linalg.pinv(S[0])@Wk
1255        ttCores[0]=ttCores_old[0]+deltaunfold@Vk.T
1256
1257        # Update G2
1258        deltaunfold=mode_n_unfolding(delta,1)
1259        Ak=mode_n_unfolding(ttCores_old[0],2)#.T
1260        Bk=coreContraction(ttCores_old[2:-1]+[g4])
1261        Bk=mode_n_unfolding(Bk,0)
1262        Wk=np.kron(Bk,Ak)
1263        S[1]=forgettingFactor*S[1]+Wk@Wk.T
1264        # Vk=np.linalg.lstsq(S[1],Wk,rcond=None)[0]
1265        try:
1266            Vk=np.linalg.solve(S[1],Wk)
1267        except np.linalg.LinAlgError:
1268            print('Singular S[1] matrix, using Moore-Penrose inverse')
1269            Vk=np.linalg.pinv(S[1])@Wk
1270        ttCores[1]=mode_n_unfolding(ttCores_old[1],1)
1271        ttCores[1]=ttCores[1]+deltaunfold@Vk.T
1272        ttCores[1]=ttCores[1].reshape(dims[1],ttRanks[0],ttRanks[1],order='F').transpose(1,0,2)
1273
1274        # Update G3
1275        deltaunfold=mode_n_unfolding(delta,2)
1276        Ak=mode_n_unfolding(coreContraction(ttCores_old[:2]),2)
1277        # Bk=coreContraction(ttCores_old[3:-1]+[g4]) #insert a mode 0 unfolding here if necessary
1278        # Bk=mode_n_unfolding(Bk,0)
1279        Wk=np.kron(g4,Ak)
1280        S[2]=forgettingFactor*S[2]+Wk@Wk.T
1281        # Vk=np.linalg.lstsq(S[2],Wk,rcond=None)[0]
1282        try:
1283            Vk=np.linalg.solve(S[2],Wk)
1284        except np.linalg.LinAlgError:
1285            print('Singular S[2] matrix, using Moore-Penrose inverse')
1286            Vk=np.linalg.pinv(S[2])@Wk
1287        ttCores[2]=mode_n_unfolding(ttCores_old[2],1)
1288        ttCores[2]=ttCores[2]+deltaunfold@Vk.T
1289        ttCores[2]=ttCores[2].reshape(dims[2],ttRanks[1],ttRanks[2],order='F').transpose(1,0,2)
1290
1291        return S,ttCores