1 | #region License Information
|
---|
2 | /* HeuristicLab
|
---|
3 | * Copyright (C) 2002-2013 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 HeuristicLab.DataImporter.Command.View;
|
---|
26 | using HeuristicLab.DataImporter.Data;
|
---|
27 | using HeuristicLab.DataImporter.Data.CommandBase;
|
---|
28 | using HeuristicLab.DataImporter.Data.Model;
|
---|
29 | using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
|
---|
30 |
|
---|
31 | namespace HeuristicLab.DataImporter.Command {
|
---|
32 | [StorableClass]
|
---|
33 | [ViewableCommandInfoAttribute("Savitzky Golay", 1, ColumnGroupState.DoubleColumnSelected, "Change Values", Position = 2,
|
---|
34 | OptionsView = typeof(FilterSavitzkyGolayCommandView))]
|
---|
35 | public class FilterSavitzkyGolayCommand : ColumnGroupCommandWithAffectedColumnsBase {
|
---|
36 | private Dictionary<int, DoubleColumn> oldColumns;
|
---|
37 | private ICollection<int> oldSortedColumnIndices;
|
---|
38 |
|
---|
39 | [Storable]
|
---|
40 | public int WindowLeft { get; set; }
|
---|
41 | [Storable]
|
---|
42 | public int WindowRight { get; set; }
|
---|
43 | [Storable]
|
---|
44 | public int Order { get; set; }
|
---|
45 | [Storable]
|
---|
46 | public int OrderOfDerivative { get; set; }
|
---|
47 |
|
---|
48 | [StorableConstructor]
|
---|
49 | protected FilterSavitzkyGolayCommand(bool deserializing)
|
---|
50 | : base(deserializing) {
|
---|
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 | }
|
---|
75 | DoubleColumn column;
|
---|
76 | oldSortedColumnIndices = new List<int>(ColumnGroup.SortedColumnIndexes);
|
---|
77 | foreach (int col in AffectedColumns) {
|
---|
78 | if (ColumnGroup.GetColumn(col) is DoubleColumn) {
|
---|
79 | column = (DoubleColumn)ColumnGroup.GetColumn(col);
|
---|
80 | oldColumns.Add(col, column);
|
---|
81 | ColumnGroup.ReplaceColumn(col, CalcNewColumn(column, this.WindowLeft, this.WindowRight, this.OrderOfDerivative, this.Order));
|
---|
82 | }
|
---|
83 | }
|
---|
84 | ColumnGroup.SortedColumnIndexes = oldSortedColumnIndices;
|
---|
85 | ColumnGroup.FireChanged();
|
---|
86 | ColumnGroup = null;
|
---|
87 | }
|
---|
88 |
|
---|
89 | public override void UndoExecute() {
|
---|
90 | base.UndoExecute();
|
---|
91 | foreach (KeyValuePair<int, DoubleColumn> pair in oldColumns)
|
---|
92 | ColumnGroup.ReplaceColumn(pair.Key, pair.Value);
|
---|
93 |
|
---|
94 | ColumnGroup.SortedColumnIndexes = oldSortedColumnIndices;
|
---|
95 | oldSortedColumnIndices = null;
|
---|
96 | oldColumns.Clear();
|
---|
97 | ColumnGroup.FireChanged();
|
---|
98 | ColumnGroup = null;
|
---|
99 | }
|
---|
100 |
|
---|
101 | private DoubleColumn CalcNewColumn(DoubleColumn oldColumn, int windowLeft, int windowRight, int derivativeOrder, int order) {
|
---|
102 | double[] c;
|
---|
103 | SavitzkyGolay(Math.Abs(windowLeft), Math.Abs(windowRight), derivativeOrder, order, out c);
|
---|
104 |
|
---|
105 | DoubleColumn newCol = (DoubleColumn)oldColumn.CreateCopyOfColumnWithoutValues();
|
---|
106 | double[] data = oldColumn.ValuesEnumerable.Cast<double?>().Select(x => x.HasValue ? x.Value : double.NaN).ToArray();
|
---|
107 | double[] ans = new double[oldColumn.TotalValuesCount];
|
---|
108 | alglib.convr1d(data, data.Length, c, c.Length, out ans);
|
---|
109 | for (int i = Math.Abs(windowRight); i < ans.Length - Math.Abs(windowLeft); i++) {
|
---|
110 | if (!double.IsNaN(ans[i]))
|
---|
111 | newCol.AddValue(ans[i]);
|
---|
112 | else
|
---|
113 | newCol.AddValueOrNull(null);
|
---|
114 | }
|
---|
115 | newCol.SortOrder = oldColumn.SortOrder;
|
---|
116 | return newCol;
|
---|
117 | }
|
---|
118 |
|
---|
119 | /// <summary>
|
---|
120 | /// 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
|
---|
121 | /// </summary>
|
---|
122 | /// <param name="nl">number of samples to the left</param>
|
---|
123 | /// <param name="nr">number of samples to the right</param>
|
---|
124 | /// <param name="ld">order of derivative (smoothing=0)</param>
|
---|
125 | /// <param name="order">order of the polynomial to fit</param>
|
---|
126 | /// <param name="c">resulting coefficients for convolution, in correct order (t-nl, ... t-1, t+0, t+1, ... t+nr)</param>
|
---|
127 | private void SavitzkyGolay(int nl, int nr, int ld, int order, out double[] c) {
|
---|
128 | int np = nl + nr + 1;
|
---|
129 |
|
---|
130 | int j, k, imj, ipj, kk, mm;
|
---|
131 | double fac = 0;
|
---|
132 | double sum = 0;
|
---|
133 | if (nl < 0 || nr < 0 || ld > order || nl + nr < order) throw new ArgumentException();
|
---|
134 |
|
---|
135 | double[,] a = new double[order + 1, order + 1];
|
---|
136 | double[] b = new double[order + 1];
|
---|
137 | c = new double[np];
|
---|
138 |
|
---|
139 | for (ipj = 0; ipj <= (order << 1); ipj++) {
|
---|
140 | sum = (ipj > 0 ? 0.0 : 1.0);
|
---|
141 | for (k = 1; k <= nr; k++) sum += Math.Pow((double)k, (double)ipj);
|
---|
142 | for (k = 1; k <= nl; k++) sum += Math.Pow((double)-k, (double)ipj);
|
---|
143 | mm = Math.Min(ipj, 2 * order - ipj);
|
---|
144 | for (imj = -mm; imj <= mm; imj += 2)
|
---|
145 | a[(ipj + imj) / 2, (ipj - imj) / 2] = sum;
|
---|
146 | }
|
---|
147 | for (j = 0; j < order + 1; j++) b[j] = 0;
|
---|
148 | b[ld] = 1.0;
|
---|
149 | alglib.densesolverreport rep;
|
---|
150 | int info;
|
---|
151 | double[] x = new double[b.Length];
|
---|
152 | alglib.rmatrixsolve(a, b.Length, b, out info, out rep, out x);
|
---|
153 |
|
---|
154 | for (kk = 0; kk < np; kk++) c[kk] = 0.0;
|
---|
155 | for (k = -nl; k <= nr; k++) {
|
---|
156 | sum = x[0];
|
---|
157 | fac = 1.0;
|
---|
158 | for (mm = 1; mm <= order; mm++) sum += x[mm] * (fac *= k);
|
---|
159 | kk = k + nl;
|
---|
160 | c[kk] = sum;
|
---|
161 | }
|
---|
162 | }
|
---|
163 | }
|
---|
164 | }
|
---|