wip

lib/src/SlicedInverseRegression.cxx

@@ -38,12 +38,87 @@ SlicedInverseRegression::SlicedInverseRegression()
   // Nothing to do
 }
 
+SlicedInverseRegression::SlicedInverseRegression(const Sample & inputSample,
+                                                 const Sample & outputSample)
+  : PersistentObject()
+  , inputSample_(inputSample)
+  , outputSample_(outputSample)
+{
+  if (outputSample.getDimension() != 1)
+    throw InvalidArgumentException(HERE) << "Supervision variable must be of dimension 1";
+}
+
 /* Virtual constructor method */
 SlicedInverseRegression * SlicedInverseRegression::clone() const
 {
   return new SlicedInverseRegression(*this);
 }
 
+void SlicedInverseRegression::run()
+{
+  // const UnsignedInteger dimension = inputSample_.getDimension();
+  const UnsignedInteger size = inputSample_.getSize();
+
+  const Indices supervisionIndices = outputSample_.argsort();
+  Collection<Indices> list_chunk;
+  UnsignedInteger offset = 0;
+  Indices chunk_population;
+  for (UnsignedInteger i = 0; i < n_slices_; ++ i)
+  {
+    const UnsignedInteger localSize = std::min(size / n_slices_, size - offset);
+    Indices chunk(localSize);
+    std::copy(supervisionIndices.begin() + offset, supervisionIndices.begin() + offset + localSize, chunk.begin());
+    list_chunk.add(chunk);
+    chunk_population.add(localSize);
+    offset += localSize;
+  }
+  std::cout << "SlicedInverseRegression::run"<< list_chunk<< " pop="<<chunk_population<< std::endl;
+
+  const Point center(inputSample_.computeMean());
+  std::cout << "SlicedInverseRegression::run center="<<center<< std::endl;
+  Sample X_centered(inputSample_);
+  X_centered -= center;
+  CovarianceMatrix input_covariance(X_centered.computeCovariance());
+  std::cout << "SlicedInverseRegression::run input_covariance="<<input_covariance<< std::endl;
+  Matrix u;
+  Matrix vT;
+#if OPENTURNS_VERSION >= 102300
+  const Point singularValues = input_covariance.computeSVDInPlace(u, vT, false); // fullSVD
+#else
+  const Point singularValues = input_covariance.computeSVD(u, vT, false, false); // fullSVD, keepIntact
+#endif
+  Point s1(singularValues.getSize());
+  for(UnsignedInteger i = 0; i < s1.getSize(); ++ i)
+    s1[i] = 1.0 / singularValues[i];
+  Matrix us1(u);
+  for(UnsignedInteger j = 0; j < s1.getSize(); ++ j)
+    for(UnsignedInteger i = 0; i < s1.getSize(); ++ i)
+      us1(i, j) *= s1[j];
+  Matrix inverseCovariance(us1 * vT);
+  std::cout << "SlicedInverseRegression::run inverseCovariance="<<inverseCovariance<< std::endl;
+  const UnsignedInteger inputDimension = inputSample_.getDimension();
+  Matrix weighted_covariance(inputDimension, inputDimension);
+  for(UnsignedInteger j = 0; j < n_slices_; ++ j)
+  {
+    const Point meanSlice(inputSample_.select(list_chunk[j]).computeMean());
+    Matrix slice_moment(inputDimension, 1);
+    for(UnsignedInteger i = 0; i < inputDimension; ++ i)
+      slice_moment(i, 0) = meanSlice[i];
+    Matrix slice_covariance(slice_moment * slice_moment.transpose());
+    const Scalar w = chunk_population[j] * 1.0 / size;
+    weighted_covariance = weighted_covariance + slice_covariance * w;
+  }
+  std::cout << "SlicedInverseRegression::run weighted_covariance="<<weighted_covariance<< std::endl;
+  SquareMatrix eigen_vector;
+  SymmetricMatrix icwc((inverseCovariance * weighted_covariance).getImplementation());
+#if OPENTURNS_VERSION >= 102300
+  const Point eigen_value = icwc.computeEVInPlace(eigen_vector);
+#else
+  const Point eigen_value = icwc.computeEV(eigen_vector, false); // keepIntact
+#endif
+  std::cout << "SlicedInverseRegression::run eigen_value="<<eigen_value<< std::endl;
+}
+
 /* example of a func that return a point squared. */
 Point SlicedInverseRegression::square(Point& p) const
 {
@@ -67,13 +142,15 @@ String SlicedInverseRegression::__repr__() const
 /* Method save() stores the object through the StorageManager */
 void SlicedInverseRegression::save(Advocate & adv) const
 {
-  PersistentObject::save( adv );
+  PersistentObject::save(adv);
+  adv.saveAttribute( "n_slices_", n_slices_ );
 }
 
 /* Method load() reloads the object from the StorageManager */
 void SlicedInverseRegression::load(Advocate & adv)
 {
-  PersistentObject::load( adv );
+  PersistentObject::load(adv);
+  adv.loadAttribute( "n_slices_", n_slices_ );
 }
 
 

lib/src/otsliced/SlicedInverseRegression.hxx

@@ -23,7 +23,7 @@
 
 #include <openturns/PersistentObject.hxx>
 #include <openturns/StorageManager.hxx>
-#include <openturns/Point.hxx>
+#include <openturns/Sample.hxx>
 #include "otsliced/otslicedprivate.hxx"
 
 namespace OTSLICED
@@ -43,12 +43,17 @@ public:
   /** Default constructor */
   SlicedInverseRegression();
 
+  SlicedInverseRegression(const OT::Sample & inputSample,
+                          const OT::Sample & outputSample);
+
   /** Virtual constructor method */
   SlicedInverseRegression * clone() const override;
 
   /** example of a func that return a point squared. **/
   OT::Point square(OT::Point& p) const;
 
+  void run();
+
   /** String converter */
   OT::String __repr__() const override;
 
@@ -59,6 +64,11 @@ public:
   void load(OT::Advocate & adv) override;
 
 private:
+  OT::Sample inputSample_;
+  OT::Sample outputSample_;
+
+  OT::UnsignedInteger n_components_ = 0;
+  OT::UnsignedInteger n_slices_ = 10;
 
 }; /* class SlicedInverseRegression */
 

python/doc/examples/plot_example1.py

@@ -21,9 +21,9 @@
 
 # %%
 import openturns as ot
-import ottemplate
+import otsliced
 
-a = ottemplate.MyClass()
+a = otsliced.SlicedInverseRegression()
 p = ot.Point([2, 3])
 squared_p = a.square(p)
 print(squared_p)

python/test/t_SlicedInverseRegression_std.py

@@ -1,14 +1,19 @@
 #!/usr/bin/env python
 
-from __future__ import print_function
 import openturns as ot
 import otsliced
+import numpy as np
 
-a = otsliced.SlicedInverseRegression()
-print(a)
+np.random.seed(0)
 
-p = ot.Point([2, 3])
-print(p)
+# Make 2D dataset
+t = np.linspace(0,1,100).reshape(-1,1)
+X = np.random.normal(size=(100, 2))*0.1 + t @ np.array([[-1,2]])
+X = X - X.mean(0)
+y = (4*(X.T[0]+2*X.T[1])+2) + np.random.normal(size=100)*0.2
+
+print(ot.Sample(X), ot.Sample.BuildFromPoint(y))
+
+algo = otsliced.SlicedInverseRegression(X, ot.Sample.BuildFromPoint(y))
+algo.run()
 
-squared_p = a.square(p)
-print(squared_p)