#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 .
*/
#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 oldColumns;
private ICollection 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();
}
public FilterSavitzkyGolayCommand(DataSet dataSet, string columnGroupName, int[] affectedColumns)
: base(dataSet, columnGroupName, affectedColumns) {
oldColumns = new Dictionary();
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(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 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().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;
}
///
/// 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
///
/// number of samples to the left
/// number of samples to the right
/// order of derivative (smoothing=0)
/// order of the polynomial to fit
/// resulting coefficients for convolution, in correct order (t-nl, ... t-1, t+0, t+1, ... t+nr)
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;
}
}
}
}