1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022011 Heuristic and Evolutionary Algorithms Laboratory (HEAL)


4  *


5  * This file is part of HeuristicLab.


6  *


7  * HeuristicLab is free software: you can redistribute it and/or modify


8  * it under the terms of the GNU General Public License as published by


9  * the Free Software Foundation, either version 3 of the License, or


10  * (at your option) any later version.


11  *


12  * HeuristicLab is distributed in the hope that it will be useful,


13  * but WITHOUT ANY WARRANTY; without even the implied warranty of


14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


15  * GNU General Public License for more details.


16  *


17  * You should have received a copy of the GNU General Public License


18  * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.


19  */


20  #endregion


21 


22  using System;


23  using System.Collections.Generic;


24  using System.Linq;


25  using System.Text;


26  using System.Xml;


27  using HeuristicLab.DataImporter.Data;


28  using HeuristicLab.DataImporter.Data.CommandBase;


29  using HeuristicLab.DataImporter.Data.Model;


30  using HeuristicLab.DataImporter.Command.View;


31  using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;


32 


33  namespace HeuristicLab.DataImporter.Command {


34  [StorableClass]


35  [ViewableCommandInfoAttribute("Savitzky Golay", 1, ColumnGroupState.DoubleColumnSelected, "Change Values", Position = 2,


36  OptionsView = typeof(FilterSavitzkyGolayCommandView))]


37  public class FilterSavitzkyGolayCommand : ColumnGroupCommandWithAffectedColumnsBase {


38  private Dictionary<int, DoubleColumn> oldColumns;


39  private ICollection<int> oldSortedColumnIndices;


40 


41  public int WindowLeft { get; set; }


42 


43  public int WindowRight { get; set; }


44 


45  public int Order { get; set; }


46 


47  public int OrderOfDerivative { get; set; }


48 


49  private FilterSavitzkyGolayCommand()


50  : base(null, string.Empty, null) {


51  oldColumns = new Dictionary<int, DoubleColumn>();


52  }


53 


54  public FilterSavitzkyGolayCommand(DataSet dataSet, string columnGroupName, int[] affectedColumns)


55  : base(dataSet, columnGroupName, affectedColumns) {


56  oldColumns = new Dictionary<int, DoubleColumn>();


57  this.WindowLeft = 16;


58  this.WindowRight = 16;


59  this.Order = 2;


60  this.OrderOfDerivative = 0;


61  }


62 


63  public override string Description {


64  get { return "Filter Column  Savitzky Golay"; }


65  }


66 


67  public override void Execute() {


68  base.Execute();


69  if (this.WindowLeft >= WindowRight)


70  throw new CommandExecutionException("Window size for filter commands must no be smaller than 1.", this);


71  foreach (int col in AffectedColumns) {


72  DoubleColumn doubleCol = ColumnGroup.GetColumn(col) as DoubleColumn;


73  if (doubleCol == null) throw new CommandExecutionException("Filtering is only supported for double columns.", this);


74  if (doubleCol.ContainsNullValues) throw new CommandExecutionException("Filtering is not supported for double columns with missing values.", this);


75  }


76  DoubleColumn column;


77  oldSortedColumnIndices = new List<int>(ColumnGroup.SortedColumnIndexes);


78  foreach (int col in AffectedColumns) {


79  if (ColumnGroup.GetColumn(col) is DoubleColumn) {


80  column = (DoubleColumn)ColumnGroup.GetColumn(col);


81  oldColumns.Add(col, column);


82  ColumnGroup.ReplaceColumn(col, CalcNewColumn(column, this.WindowLeft, this.WindowRight, this.OrderOfDerivative, this.Order));


83  }


84  }


85  ColumnGroup.SortedColumnIndexes = oldSortedColumnIndices;


86  ColumnGroup.FireChanged();


87  ColumnGroup = null;


88  }


89 


90  public override void UndoExecute() {


91  base.UndoExecute();


92  foreach (KeyValuePair<int, DoubleColumn> pair in oldColumns)


93  ColumnGroup.ReplaceColumn(pair.Key, pair.Value);


94 


95  ColumnGroup.SortedColumnIndexes = oldSortedColumnIndices;


96  oldSortedColumnIndices = null;


97  oldColumns.Clear();


98  ColumnGroup.FireChanged();


99  ColumnGroup = null;


100  }


101 


102  private DoubleColumn CalcNewColumn(DoubleColumn oldColumn, int windowLeft, int windowRight, int derivativeOrder, int order) {


103  double[] c;


104  SavitzkyGolay(Math.Abs(windowLeft), Math.Abs(windowRight), derivativeOrder, order, out c);


105 


106  DoubleColumn newCol = (DoubleColumn)oldColumn.CreateCopyOfColumnWithoutValues();


107  double[] data = oldColumn.ValuesEnumerable.Cast<double>().ToArray();


108  double[] ans = new double[oldColumn.TotalValuesCount];


109  alglib.convr1d(data, data.Length, c, c.Length, out ans);


110  for (int i = Math.Abs(windowRight); i < ans.Length  Math.Abs(windowLeft); i++) {


111  newCol.AddValue(ans[i]);


112  }


113  newCol.SortOrder = oldColumn.SortOrder;


114  return newCol;


115  }


116 


117  /// <summary>


118  /// Calculates coefficients for SavitzkyGolay filtering. (Numerical Recipes, page 769), one important change is that the coefficients are returned in normal order instead of wraparound order


119  /// </summary>


120  /// <param name="nl">number of samples to the left</param>


121  /// <param name="nr">number of samples to the right</param>


122  /// <param name="ld">order of derivative (smoothing=0)</param>


123  /// <param name="order">order of the polynomial to fit</param>


124  /// <param name="c">resulting coefficients for convolution, in correct order (tnl, ... t1, t+0, t+1, ... t+nr)</param>


125  private void SavitzkyGolay(int nl, int nr, int ld, int order, out double[] c) {


126  int np = nl + nr + 1;


127 


128  int j, k, imj, ipj, kk, mm;


129  double fac = 0;


130  double sum = 0;


131  if (nl < 0  nr < 0  ld > order  nl + nr < order) throw new ArgumentException();


132 


133  double[,] a = new double[order + 1, order + 1];


134  double[] b = new double[order + 1];


135  c = new double[np];


136 


137  for (ipj = 0; ipj <= (order << 1); ipj++) {


138  sum = (ipj > 0 ? 0.0 : 1.0);


139  for (k = 1; k <= nr; k++) sum += Math.Pow((double)k, (double)ipj);


140  for (k = 1; k <= nl; k++) sum += Math.Pow((double)k, (double)ipj);


141  mm = Math.Min(ipj, 2 * order  ipj);


142  for (imj = mm; imj <= mm; imj += 2)


143  a[(ipj + imj) / 2, (ipj  imj) / 2] = sum;


144  }


145  for (j = 0; j < order + 1; j++) b[j] = 0;


146  b[ld] = 1.0;


147  alglib.densesolverreport rep;


148  int info;


149  double[] x = new double[b.Length];


150  alglib.rmatrixsolve(a, b.Length, b, out info, out rep, out x);


151 


152  for (kk = 0; kk < np; kk++) c[kk] = 0.0;


153  for (k = nl; k <= nr; k++) {


154  sum = x[0];


155  fac = 1.0;


156  for (mm = 1; mm <= order; mm++) sum += x[mm] * (fac *= k);


157  kk = k + nl;


158  c[kk] = sum;


159  }


160  }


161  }


162  }

