New Bonsai.ML.Pca.Torch Package for Batched and Online PCA

Bonsai.ML.sln

@@ -46,6 +46,10 @@ Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Bonsai.ML.Lds.Torch.Tests",
 EndProject
 Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Bonsai.ML.Lds.Torch.Design", "src\Bonsai.ML.Lds.Torch.Design\Bonsai.ML.Lds.Torch.Design.csproj", "{1F52DECD-1B2C-4F6C-996C-14C715283B80}"
 EndProject
+Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Bonsai.ML.Pca.Torch", "src\Bonsai.ML.Pca.Torch\Bonsai.ML.Pca.Torch.csproj", "{1A2DEED3-795E-4C28-9C5E-BA3D76B2A485}"
+EndProject
+Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Bonsai.ML.Pca.Torch.Tests", "tests\Bonsai.ML.Pca.Torch.Tests\Bonsai.ML.Pca.Torch.Tests.csproj", "{4ABCC6B2-024A-450F-85CB-2A9B2D2D2A10}"
+EndProject
 Global
 	GlobalSection(SolutionConfigurationPlatforms) = preSolution
 		Debug|Any CPU = Debug|Any CPU
@@ -120,6 +124,14 @@ Global
 		{1F52DECD-1B2C-4F6C-996C-14C715283B80}.Debug|Any CPU.Build.0 = Debug|Any CPU
 		{1F52DECD-1B2C-4F6C-996C-14C715283B80}.Release|Any CPU.ActiveCfg = Release|Any CPU
 		{1F52DECD-1B2C-4F6C-996C-14C715283B80}.Release|Any CPU.Build.0 = Release|Any CPU
+		{1A2DEED3-795E-4C28-9C5E-BA3D76B2A485}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{1A2DEED3-795E-4C28-9C5E-BA3D76B2A485}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{1A2DEED3-795E-4C28-9C5E-BA3D76B2A485}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{1A2DEED3-795E-4C28-9C5E-BA3D76B2A485}.Release|Any CPU.Build.0 = Release|Any CPU
+		{4ABCC6B2-024A-450F-85CB-2A9B2D2D2A10}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{4ABCC6B2-024A-450F-85CB-2A9B2D2D2A10}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{4ABCC6B2-024A-450F-85CB-2A9B2D2D2A10}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{4ABCC6B2-024A-450F-85CB-2A9B2D2D2A10}.Release|Any CPU.Build.0 = Release|Any CPU
 	EndGlobalSection
 	GlobalSection(SolutionProperties) = preSolution
 		HideSolutionNode = FALSE

src/Bonsai.ML.Pca.Torch/Bonsai.ML.Pca.Torch.csproj

@@ -0,0 +1,15 @@
+<Project Sdk="Microsoft.NET.Sdk">
+
+  <PropertyGroup>
+    <Description>Bonsai.ML.Pca.Torch Bonsai library.</Description>
+    <PackageTags>$(PackageTags) PCA Principal Component Analysis</PackageTags>
+    <TargetFrameworks>net472;netstandard2.0</TargetFrameworks>
+    <Nullable>enable</Nullable>
+  </PropertyGroup>
+
+  <ItemGroup>
+    <ProjectReference Include="..\Bonsai.ML\Bonsai.ML.csproj" />
+    <ProjectReference Include="..\Bonsai.ML.Torch\Bonsai.ML.Torch.csproj" />
+  </ItemGroup>
+
+</Project>

src/Bonsai.ML.Pca.Torch/CreatePca.cs

@@ -0,0 +1,117 @@
+using System;
+using System.ComponentModel;
+using System.Collections.Generic;
+using System.Linq;
+using System.Linq.Expressions;
+using System.Reactive.Linq;
+using System.Reflection;
+using System.Xml.Serialization;
+using Bonsai.Expressions;
+
+namespace Bonsai.ML.Pca.Torch;
+
+/// <summary>
+/// Represents an operator that creates a PCA model.
+/// </summary>
+[XmlInclude(typeof(Pca))]
+[XmlInclude(typeof(ProbabilisticPca))]
+[XmlInclude(typeof(OnlineProbabilisticPca))]
+[XmlInclude(typeof(OnlinePcaGha))]
+[WorkflowElementCategory(ElementCategory.Source)]
+[Description("Creates a PCA model.")]
+public class CreatePca : ZeroArgumentExpressionBuilder, INamedElement, ICustomTypeDescriptor
+{
+    string INamedElement.Name => $"CreatePca.{ModelType}";
+
+    /// <summary>
+    /// Gets or sets the PCA model to create.
+    /// </summary>
+    [Browsable(false)]
+    public PcaBaseModel Model { get; set; } = new Pca();
+
+    /// <summary>
+    /// Gets or sets the type of PCA model to create.
+    /// </summary>
+    [XmlIgnore]
+    [RefreshProperties(RefreshProperties.All)]
+    [Description("The type of PCA model to create.")]
+    [Category("Combinator")]
+    public PcaModelType ModelType
+    {
+        get => Model switch
+        {
+            ProbabilisticPca => PcaModelType.ProbabilisticPca,
+            OnlineProbabilisticPca => PcaModelType.OnlineProbabilisticPca,
+            OnlinePcaGha => PcaModelType.OnlinePcaGha,
+            _ => PcaModelType.Pca,
+        };
+        set
+        {
+            if (value == ModelType) return;
+            Model = value switch
+            {
+                PcaModelType.ProbabilisticPca => new ProbabilisticPca(),
+                PcaModelType.OnlineProbabilisticPca => new OnlineProbabilisticPca(),
+                PcaModelType.OnlinePcaGha => new OnlinePcaGha(),
+                _ => new Pca(),
+            };
+        }
+    }
+
+    /// <inheritdoc/>
+    public override Expression Build(IEnumerable<Expression> arguments)
+    {
+        var processMethod = typeof(CreatePca).GetMethod(
+            nameof(Process),
+            BindingFlags.Static | BindingFlags.NonPublic);
+
+        var modelType = Model.GetType();
+        var genericMethod = processMethod.MakeGenericMethod(modelType);
+
+        return Expression.Call(genericMethod, Expression.Constant(Model));
+    }
+
+    static IObservable<T> Process<T>(T model) where T : PcaBaseModel
+    {
+        return Observable.Using(
+            () => model,
+            m => Observable.Return(m)
+                .Concat(Observable.Never(m))
+                .Finally(() => m.Dispose())
+        );
+    }
+
+    PropertyDescriptorCollection ICustomTypeDescriptor.GetProperties(Attribute[]? attributes)
+    {
+        var properties = TypeDescriptor.GetProperties(this, attributes, true).Cast<PropertyDescriptor>();
+        var modelProperties = TypeDescriptor.GetProperties(Model, attributes)
+            .Cast<PropertyDescriptor>()
+            .Where(descriptor => !descriptor.IsReadOnly);
+        return new PropertyDescriptorCollection(properties.Concat(modelProperties).ToArray());
+    }
+
+    PropertyDescriptorCollection ICustomTypeDescriptor.GetProperties() =>
+        ((ICustomTypeDescriptor)this).GetProperties(Array.Empty<Attribute>());
+
+    object ICustomTypeDescriptor.GetPropertyOwner(PropertyDescriptor? pd) =>
+        pd == null || pd.ComponentType.IsInstanceOfType(Model) ? Model : this;
+
+    AttributeCollection ICustomTypeDescriptor.GetAttributes() => TypeDescriptor.GetAttributes(this, true);
+
+    string? ICustomTypeDescriptor.GetClassName() => TypeDescriptor.GetClassName(this, true);
+
+    TypeConverter ICustomTypeDescriptor.GetConverter() => TypeDescriptor.GetConverter(this, true);
+
+    EventDescriptor? ICustomTypeDescriptor.GetDefaultEvent() => TypeDescriptor.GetDefaultEvent(this, true);
+
+    PropertyDescriptor? ICustomTypeDescriptor.GetDefaultProperty() => TypeDescriptor.GetDefaultProperty(this, true);
+
+    object? ICustomTypeDescriptor.GetEditor(Type editorBaseType) => TypeDescriptor.GetEditor(this, editorBaseType, true);
+
+    EventDescriptorCollection ICustomTypeDescriptor.GetEvents(Attribute[]? attributes) =>
+        TypeDescriptor.GetEvents(this, attributes, true);
+
+    EventDescriptorCollection ICustomTypeDescriptor.GetEvents() => TypeDescriptor.GetEvents(this, true);
+
+    string? ICustomTypeDescriptor.GetComponentName() => TypeDescriptor.GetComponentName(this, true);
+}

src/Bonsai.ML.Pca.Torch/Fit.cs

@@ -0,0 +1,46 @@
+using System;
+using System.ComponentModel;
+using System.Reactive.Linq;
+using static TorchSharp.torch;
+
+namespace Bonsai.ML.Pca.Torch;
+
+/// <summary>
+/// Represents an operator that fits a PCA model to the input data.
+/// </summary>
+[Combinator]
+[Description("Fits a PCA model to the input data.")]
+[WorkflowElementCategory(ElementCategory.Sink)]
+public class Fit
+{
+    private void FitModel(PcaBaseModel model, Tensor data)
+    {
+        model.Fit(data);
+    }
+
+    /// <summary>
+    /// Fits a PCA model to the input data.
+    /// </summary>
+    /// <param name="source"></param>
+    /// <returns></returns>
+    public IObservable<Tuple<T, Tensor>> Process<T>(IObservable<Tuple<T, Tensor>> source) where T : PcaBaseModel
+    {
+        return source.Do((value) =>
+        {
+            FitModel(value.Item1, value.Item2);
+        });
+    }
+
+    /// <summary>
+    /// Fits a PCA model to the input data.
+    /// </summary>
+    /// <param name="source"></param>
+    /// <returns></returns>
+    public IObservable<Tuple<Tensor, T>> Process<T>(IObservable<Tuple<Tensor, T>> source) where T : PcaBaseModel
+    {
+        return source.Do((value) =>
+        {
+            FitModel(value.Item2, value.Item1);
+        });
+    }
+}

src/Bonsai.ML.Pca.Torch/FitAndTransform.cs

@@ -0,0 +1,48 @@
+using System;
+using System.ComponentModel;
+using System.Reactive.Linq;
+using static TorchSharp.torch;
+
+namespace Bonsai.ML.Pca.Torch;
+
+/// <summary>
+/// Represents an operator that fits a PCA model to the input data and transforms it.
+/// </summary>
+[Combinator]
+[Description("Fits a PCA model to the input data and transforms it.")]
+[WorkflowElementCategory(ElementCategory.Transform)]
+public class FitAndTransform
+{
+    private static Tensor FitModelAndTransformData(PcaBaseModel model, Tensor data)
+    {
+        return model.FitAndTransform(data);
+    }
+
+    /// <summary>
+    /// Fits a PCA model to the input data and transforms it.
+    /// </summary>
+    /// <param name="source"></param>
+    /// <returns></returns>
+    public IObservable<Tuple<T, Tensor>> Process<T>(IObservable<Tuple<T, Tensor>> source) where T : PcaBaseModel
+    {
+        return source.Select((value) =>
+        {
+            var transformed = FitModelAndTransformData(value.Item1, value.Item2);
+            return Tuple.Create(value.Item1, transformed);
+        });
+    }
+
+    /// <summary>
+    /// Fits a PCA model to the input data and transforms it.
+    /// </summary>
+    /// <param name="source"></param>
+    /// <returns></returns>
+    public IObservable<Tuple<Tensor, T>> Process<T>(IObservable<Tuple<Tensor, T>> source) where T : PcaBaseModel
+    {
+        return source.Select((value) =>
+        {
+            var transformed = FitModelAndTransformData(value.Item2, value.Item1);
+            return Tuple.Create(transformed, value.Item2);
+        });
+    }
+}

src/Bonsai.ML.Pca.Torch/OnlinePcaGha.cs

@@ -0,0 +1,106 @@
+using System.ComponentModel;
+using System.Xml.Serialization;
+using static TorchSharp.torch;
+
+namespace Bonsai.ML.Pca.Torch;
+
+/// <summary>
+/// A streaming/online PCA model based on the Generalized Hebbian Algorithm (GHA).
+/// </summary>
+public class OnlinePcaGha : PcaBaseModel
+{
+    /// <summary>
+    /// Gets the number of samples that have been used to fit the model.
+    /// </summary>
+    [XmlIgnore]
+    public int SampleCount { get; private set; } = 0;
+
+    /// <summary>
+    /// Gets the mean of the fitted data.
+    /// </summary>
+    [XmlIgnore]
+    public Tensor Mean { get; private set; } = empty(0);
+
+    /// <summary>
+    /// Gets or sets the learning rate.
+    /// </summary>
+    [Category("ModelParameters")]
+    [Description("The learning rate used when updating the components.")]
+    public double LearningRate { get; set; } = 0.1;
+
+    /// <summary>
+    /// Gets or sets the random number generator used for initializing the model.
+    /// </summary>
+    [XmlIgnore]
+    [Description("The random number generator used for initializing the model.")]
+    public Generator? Generator { get; set; }
+
+    /// <inheritdoc/>
+    public override void Fit(Tensor data)
+    {
+        base.Fit(data);
+
+        var numSamples = data.size(0);
+
+        using (no_grad())
+        using (NewDisposeScope())
+        {
+            // Initialize components randomly
+            if (Components.numel() == 0)
+                Components = randn([NumFeatures, NumComponents], dtype: Type, device: Device, generator: Generator);
+
+            if (Mean.numel() == 0)
+                Mean = data.mean([0], keepdim: true);
+            else
+            {
+                Mean *= (double)SampleCount / (SampleCount + numSamples);
+                Mean += data.mean([0], keepdim: true) * numSamples / (SampleCount + numSamples);
+            }
+
+            SampleCount += (int)numSamples;
+            var dataCentered = data - Mean;
+
+            var projection = dataCentered.matmul(Components);
+            var hebbianTerm = dataCentered.T.matmul(projection);
+            var crossTerm = projection.T.matmul(projection);
+            var upperTriangular = crossTerm.triu(0);
+            var correlation = Components.matmul(upperTriangular);
+            var componentsUpdate = (hebbianTerm - correlation) * (LearningRate / numSamples);
+            var weights = Components + componentsUpdate;
+            var norms = weights.norm(dim: 0, keepdim: true, p: 2).clamp_min(1e-12);
+
+            Components = linalg.qr(weights / norms, mode: linalg.QRMode.Reduced).Q.MoveToOuterDisposeScope();
+            Mean = Mean.MoveToOuterDisposeScope();
+        }
+    }
+
+    /// <inheritdoc/>
+    public override Tensor Transform(Tensor data)
+    {
+        base.Transform(data);
+        var dataCentered = data - Mean;
+        return dataCentered.matmul(Components);
+    }
+
+    /// <inheritdoc/>
+    public override Tensor Reconstruct(Tensor data)
+    {
+        base.Reconstruct(data);
+        return data.matmul(Components.T) + Mean;
+    }
+
+    /// <inheritdoc/>
+    public override Tensor FitAndTransform(Tensor data)
+    {
+        Fit(data);
+        return Transform(data);
+    }
+
+    /// <inheritdoc/>
+    public override void Dispose()
+    {
+        base.Dispose();
+        Mean = Utils.DisposeAndReset(Mean);
+        SampleCount = 0;
+    }
+}