source: branches/HeuristicLab.DataImporter/HeuristicLab.DataImporter.Command/ChangeValues/FilterSavitzkyGolayCommand.cs @ 7267

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2012 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
 *
 * This file is part of HeuristicLab.
 *
 * HeuristicLab is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * HeuristicLab is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
 */
#endregion

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Xml;
using HeuristicLab.DataImporter.Data;
using HeuristicLab.DataImporter.Data.CommandBase;
using HeuristicLab.DataImporter.Data.Model;
using HeuristicLab.DataImporter.Command.View;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;

namespace HeuristicLab.DataImporter.Command {
  [StorableClass]
  [ViewableCommandInfoAttribute("Savitzky Golay", 1, ColumnGroupState.DoubleColumnSelected, "Change Values", Position = 2,
    OptionsView = typeof(FilterSavitzkyGolayCommandView))]
  public class FilterSavitzkyGolayCommand : ColumnGroupCommandWithAffectedColumnsBase {
    private Dictionary<int, DoubleColumn> oldColumns;
    private ICollection<int> oldSortedColumnIndices;

    public int WindowLeft { get; set; }

    public int WindowRight { get; set; }

    public int Order { get; set; }

    public int OrderOfDerivative { get; set; }

    private FilterSavitzkyGolayCommand()
      : base(null, string.Empty, null) {
      oldColumns = new Dictionary<int, DoubleColumn>();
    }

    public FilterSavitzkyGolayCommand(DataSet dataSet, string columnGroupName, int[] affectedColumns)
      : base(dataSet, columnGroupName, affectedColumns) {
      oldColumns = new Dictionary<int, DoubleColumn>();
      this.WindowLeft = -16;
      this.WindowRight = 16;
      this.Order = 2;
      this.OrderOfDerivative = 0;
    }

    public override string Description {
      get { return "Filter Column - Savitzky Golay"; }
    }

    public override void Execute() {
      base.Execute();
      if (this.WindowLeft >= WindowRight)
        throw new CommandExecutionException("Window size for filter commands must no be smaller than 1.", this);
      foreach (int col in AffectedColumns) {
        DoubleColumn doubleCol = ColumnGroup.GetColumn(col) as DoubleColumn;
        if (doubleCol == null) throw new CommandExecutionException("Filtering is only supported for double columns.", this);
      }
      DoubleColumn column;
      oldSortedColumnIndices = new List<int>(ColumnGroup.SortedColumnIndexes);
      foreach (int col in AffectedColumns) {
        if (ColumnGroup.GetColumn(col) is DoubleColumn) {
          column = (DoubleColumn)ColumnGroup.GetColumn(col);
          oldColumns.Add(col, column);
          ColumnGroup.ReplaceColumn(col, CalcNewColumn(column, this.WindowLeft, this.WindowRight, this.OrderOfDerivative, this.Order));
        }
      }
      ColumnGroup.SortedColumnIndexes = oldSortedColumnIndices;
      ColumnGroup.FireChanged();
      ColumnGroup = null;
    }

    public override void UndoExecute() {
      base.UndoExecute();
      foreach (KeyValuePair<int, DoubleColumn> pair in oldColumns)
        ColumnGroup.ReplaceColumn(pair.Key, pair.Value);

      ColumnGroup.SortedColumnIndexes = oldSortedColumnIndices;
      oldSortedColumnIndices = null;
      oldColumns.Clear();
      ColumnGroup.FireChanged();
      ColumnGroup = null;
    }

    private DoubleColumn CalcNewColumn(DoubleColumn oldColumn, int windowLeft, int windowRight, int derivativeOrder, int order) {
      double[] c;
      SavitzkyGolay(Math.Abs(windowLeft), Math.Abs(windowRight), derivativeOrder, order, out c);

      DoubleColumn newCol = (DoubleColumn)oldColumn.CreateCopyOfColumnWithoutValues();
      double[] data = oldColumn.ValuesEnumerable.Cast<double?>().Select(x => x.HasValue ? x.Value : double.NaN).ToArray();
      double[] ans = new double[oldColumn.TotalValuesCount];
      alglib.convr1d(data, data.Length, c, c.Length, out ans);
      for (int i = Math.Abs(windowRight); i < ans.Length - Math.Abs(windowLeft); i++) {
        if (!double.IsNaN(ans[i]))
          newCol.AddValue(ans[i]);
        else
          newCol.AddValueOrNull(null);
      }
      newCol.SortOrder = oldColumn.SortOrder;
      return newCol;
    }

    /// <summary>
    /// Calculates coefficients for Savitzky-Golay filtering. (Numerical Recipes, page 769), one important change is that the coefficients are returned in normal order instead of wraparound order
    /// </summary>
    /// <param name="nl">number of samples to the left</param>
    /// <param name="nr">number of samples to the right</param>
    /// <param name="ld">order of derivative (smoothing=0)</param>
    /// <param name="order">order of the polynomial to fit</param>
    /// <param name="c">resulting coefficients for convolution, in correct order (t-nl, ... t-1, t+0, t+1, ... t+nr)</param>
    private void SavitzkyGolay(int nl, int nr, int ld, int order, out double[] c) {
      int np = nl + nr + 1;

      int j, k, imj, ipj, kk, mm;
      double fac = 0;
      double sum = 0;
      if (nl < 0 || nr < 0 || ld > order || nl + nr < order) throw new ArgumentException();

      double[,] a = new double[order + 1, order + 1];
      double[] b = new double[order + 1];
      c = new double[np];

      for (ipj = 0; ipj <= (order << 1); ipj++) {
        sum = (ipj > 0 ? 0.0 : 1.0);
        for (k = 1; k <= nr; k++) sum += Math.Pow((double)k, (double)ipj);
        for (k = 1; k <= nl; k++) sum += Math.Pow((double)-k, (double)ipj);
        mm = Math.Min(ipj, 2 * order - ipj);
        for (imj = -mm; imj <= mm; imj += 2)
          a[(ipj + imj) / 2, (ipj - imj) / 2] = sum;
      }
      for (j = 0; j < order + 1; j++) b[j] = 0;
      b[ld] = 1.0;
      alglib.densesolverreport rep;
      int info;
      double[] x = new double[b.Length];
      alglib.rmatrixsolve(a, b.Length, b, out info, out rep, out x);

      for (kk = 0; kk < np; kk++) c[kk] = 0.0;
      for (k = -nl; k <= nr; k++) {
        sum = x[0];
        fac = 1.0;
        for (mm = 1; mm <= order; mm++) sum += x[mm] * (fac *= k);
        kk = k + nl;
        c[kk] = sum;
      }
    }
  }
}