114 |
16 Jun 06 |
enell |
1 |
/* |
114 |
16 Jun 06 |
enell |
$Id$ |
114 |
16 Jun 06 |
enell |
3 |
|
114 |
16 Jun 06 |
enell |
Copyright (C) 2006 Johan Enell |
114 |
16 Jun 06 |
enell |
5 |
|
114 |
16 Jun 06 |
enell |
This file is part of BASE - BioArray Software Environment. |
114 |
16 Jun 06 |
enell |
Available at http://base.thep.lu.se/ |
114 |
16 Jun 06 |
enell |
8 |
|
114 |
16 Jun 06 |
enell |
BASE is free software; you can redistribute it and/or modify it |
114 |
16 Jun 06 |
enell |
under the terms of the GNU General Public License as published by |
114 |
16 Jun 06 |
enell |
the Free Software Foundation; either version 2 of the License, or |
114 |
16 Jun 06 |
enell |
(at your option) any later version. |
114 |
16 Jun 06 |
enell |
13 |
|
114 |
16 Jun 06 |
enell |
BASE is distributed in the hope that it will be useful, but |
114 |
16 Jun 06 |
enell |
WITHOUT ANY WARRANTY; without even the implied warranty of |
114 |
16 Jun 06 |
enell |
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
114 |
16 Jun 06 |
enell |
General Public License for more details. |
114 |
16 Jun 06 |
enell |
18 |
|
114 |
16 Jun 06 |
enell |
You should have received a copy of the GNU General Public License |
114 |
16 Jun 06 |
enell |
along with this program; if not, write to the Free Software |
114 |
16 Jun 06 |
enell |
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA |
114 |
16 Jun 06 |
enell |
02111-1307, USA. |
114 |
16 Jun 06 |
enell |
23 |
*/ |
114 |
16 Jun 06 |
enell |
24 |
package lowess; |
114 |
16 Jun 06 |
enell |
25 |
|
114 |
16 Jun 06 |
enell |
26 |
import java.util.ArrayList; |
114 |
16 Jun 06 |
enell |
27 |
import java.util.Arrays; |
114 |
16 Jun 06 |
enell |
28 |
import java.util.Collections; |
114 |
16 Jun 06 |
enell |
29 |
import java.util.Comparator; |
114 |
16 Jun 06 |
enell |
30 |
import java.util.Iterator; |
114 |
16 Jun 06 |
enell |
31 |
import java.util.List; |
114 |
16 Jun 06 |
enell |
32 |
|
114 |
16 Jun 06 |
enell |
33 |
import basefile.BASEFileException; |
114 |
16 Jun 06 |
enell |
34 |
import basefile.BASEFileSpotSection; |
114 |
16 Jun 06 |
enell |
35 |
|
130 |
09 Aug 06 |
enell |
36 |
/** |
130 |
09 Aug 06 |
enell |
* This is a normalizer using the lowess algorithm. This implementation has been |
130 |
09 Aug 06 |
enell |
* enhanced with the functionality of spike and pin-based normalization. |
130 |
09 Aug 06 |
enell |
39 |
* |
130 |
09 Aug 06 |
enell |
* @author Johan Enell, johan.enell@med.lu.se, Dept Oncology, Lund University |
130 |
09 Aug 06 |
enell |
41 |
*/ |
114 |
16 Jun 06 |
enell |
42 |
public class Normalizer |
114 |
16 Jun 06 |
enell |
43 |
{ |
114 |
16 Jun 06 |
enell |
44 |
private final double ywRangeFactor = 6; |
114 |
16 Jun 06 |
enell |
45 |
|
130 |
09 Aug 06 |
enell |
46 |
/** |
130 |
09 Aug 06 |
enell |
* The smoother span. This gives the proportion of spots which influence the |
130 |
09 Aug 06 |
enell |
* smooth at each value. Larger values give more smoothness. |
130 |
09 Aug 06 |
enell |
49 |
*/ |
114 |
16 Jun 06 |
enell |
50 |
private final Float fitFraction; |
114 |
16 Jun 06 |
enell |
51 |
|
130 |
09 Aug 06 |
enell |
52 |
/** |
130 |
09 Aug 06 |
enell |
* The number of robustifying iterations which should be performed. Using |
130 |
09 Aug 06 |
enell |
* smaller values of iter will make lowess run faster. |
130 |
09 Aug 06 |
enell |
55 |
*/ |
130 |
09 Aug 06 |
enell |
56 |
private final Integer iter; |
114 |
16 Jun 06 |
enell |
57 |
|
130 |
09 Aug 06 |
enell |
58 |
/** |
130 |
09 Aug 06 |
enell |
* Values of x which lie within delta of each other are replaced by a single |
130 |
09 Aug 06 |
enell |
* value in the output from lowess. Defaults to 1/100th of the range of x. |
130 |
09 Aug 06 |
enell |
61 |
*/ |
114 |
16 Jun 06 |
enell |
62 |
private final Float delta; |
114 |
16 Jun 06 |
enell |
63 |
|
130 |
09 Aug 06 |
enell |
64 |
/** |
130 |
09 Aug 06 |
enell |
* When using spike normalization this list of reporters is used to |
130 |
09 Aug 06 |
enell |
* include/exclude spots for the lowess computation. |
130 |
09 Aug 06 |
enell |
67 |
*/ |
130 |
09 Aug 06 |
enell |
68 |
private final List<String> reporters; |
130 |
09 Aug 06 |
enell |
69 |
|
130 |
09 Aug 06 |
enell |
70 |
/** |
130 |
09 Aug 06 |
enell |
* When using spike normalization the mode tells the algorithm if the |
130 |
09 Aug 06 |
enell |
* reporters in the reporterlist should be used to include or exclude spots |
130 |
09 Aug 06 |
enell |
* in the lowess computation. |
130 |
09 Aug 06 |
enell |
74 |
*/ |
114 |
16 Jun 06 |
enell |
75 |
private final ReporterMode mode; |
114 |
16 Jun 06 |
enell |
76 |
|
130 |
09 Aug 06 |
enell |
77 |
/** |
130 |
09 Aug 06 |
enell |
* Creates a normalizer with the specified parameters. |
130 |
09 Aug 06 |
enell |
79 |
* |
130 |
09 Aug 06 |
enell |
* @param fitFraction |
130 |
09 Aug 06 |
enell |
* the smoother span. This gives the proportion of spots which |
130 |
09 Aug 06 |
enell |
* influence the smooth at each value. Larger values give more |
130 |
09 Aug 06 |
enell |
* smoothness. |
130 |
09 Aug 06 |
enell |
* @param delta |
130 |
09 Aug 06 |
enell |
* values of x which lie within delta of each other are replaced by a |
130 |
09 Aug 06 |
enell |
* single value in the output from lowess. Defaults to 1/100th of the |
130 |
09 Aug 06 |
enell |
* range of x. |
130 |
09 Aug 06 |
enell |
* @param iter |
130 |
09 Aug 06 |
enell |
* the number of robustifying iterations which should be performed. |
130 |
09 Aug 06 |
enell |
* Using smaller values of iter will make lowess run faster. |
130 |
09 Aug 06 |
enell |
* @param mode |
130 |
09 Aug 06 |
enell |
* when using spike normalization the mode tells the algorithm if the |
130 |
09 Aug 06 |
enell |
* reporters in the reporterlist should be used to include or exclude |
130 |
09 Aug 06 |
enell |
* spots in the lowess computation. |
130 |
09 Aug 06 |
enell |
* @param reporters |
130 |
09 Aug 06 |
enell |
* when using spike normalization this list of reporters is used to |
130 |
09 Aug 06 |
enell |
* include/exclude spots for the lowess computation. |
130 |
09 Aug 06 |
enell |
98 |
*/ |
114 |
16 Jun 06 |
enell |
99 |
public Normalizer(Float fitFraction, Float delta, Integer nIter, String mode, String[] reporters) |
114 |
16 Jun 06 |
enell |
100 |
{ |
114 |
16 Jun 06 |
enell |
101 |
this.fitFraction = fitFraction; |
114 |
16 Jun 06 |
enell |
102 |
this.delta = delta; |
130 |
09 Aug 06 |
enell |
103 |
this.iter = nIter; |
114 |
16 Jun 06 |
enell |
104 |
this.mode = ReporterMode.get(mode); |
114 |
16 Jun 06 |
enell |
105 |
this.reporters = Arrays.asList(reporters); |
114 |
16 Jun 06 |
enell |
106 |
} |
114 |
16 Jun 06 |
enell |
107 |
|
130 |
09 Aug 06 |
enell |
108 |
/** |
130 |
09 Aug 06 |
enell |
* This method will normalize the spot section. It will start by sorting the |
130 |
09 Aug 06 |
enell |
* spots by group en each group by A. Then it will normalize each group. If |
130 |
09 Aug 06 |
enell |
* spike normalization is choosen the reporters specified in the reportelist |
130 |
09 Aug 06 |
enell |
* is either included or excluded from the lowess computation. Each spot |
130 |
09 Aug 06 |
enell |
* will be fitted to the closest spot in the lowess smooth curve. |
130 |
09 Aug 06 |
enell |
114 |
* |
130 |
09 Aug 06 |
enell |
* @param bfss |
130 |
09 Aug 06 |
enell |
* the spot section to be normalized |
130 |
09 Aug 06 |
enell |
* @throws BASEFileException |
130 |
09 Aug 06 |
enell |
118 |
*/ |
114 |
16 Jun 06 |
enell |
119 |
public void normalize(BASEFileSpotSection<Reporter, Spot> bfss) |
114 |
16 Jun 06 |
enell |
120 |
throws BASEFileException |
114 |
16 Jun 06 |
enell |
121 |
{ |
114 |
16 Jun 06 |
enell |
122 |
bfss.sortSpot(new Comparator<Spot>() |
114 |
16 Jun 06 |
enell |
123 |
{ |
114 |
16 Jun 06 |
enell |
124 |
public int compare(Spot s1, Spot s2) |
114 |
16 Jun 06 |
enell |
125 |
{ |
114 |
16 Jun 06 |
enell |
126 |
if (s1.getBlock() < s2.getBlock()) return -1; |
114 |
16 Jun 06 |
enell |
127 |
else if (s1.getBlock() == s2.getBlock()) |
114 |
16 Jun 06 |
enell |
128 |
{ |
114 |
16 Jun 06 |
enell |
129 |
if (s1.getA() < s2.getA()) return -1; |
114 |
16 Jun 06 |
enell |
130 |
else if (s1.getA() == s2.getA()) |
114 |
16 Jun 06 |
enell |
131 |
{ |
114 |
16 Jun 06 |
enell |
132 |
return 0; |
114 |
16 Jun 06 |
enell |
133 |
} |
114 |
16 Jun 06 |
enell |
134 |
else return 1; |
114 |
16 Jun 06 |
enell |
135 |
} |
114 |
16 Jun 06 |
enell |
136 |
else return 1; |
114 |
16 Jun 06 |
enell |
137 |
} |
114 |
16 Jun 06 |
enell |
138 |
}); |
114 |
16 Jun 06 |
enell |
139 |
|
114 |
16 Jun 06 |
enell |
140 |
Iterator<List<Spot>> groupIterator = new GroupIterator(bfss); |
114 |
16 Jun 06 |
enell |
141 |
while (groupIterator.hasNext()) |
114 |
16 Jun 06 |
enell |
142 |
{ |
114 |
16 Jun 06 |
enell |
143 |
List<Spot> group = groupIterator.next(); |
114 |
16 Jun 06 |
enell |
144 |
if (group.size() == 0) |
114 |
16 Jun 06 |
enell |
145 |
{ |
114 |
16 Jun 06 |
enell |
146 |
throw new BASEFileException("One of the groups is to small. Increas the groupsize."); |
114 |
16 Jun 06 |
enell |
147 |
} |
114 |
16 Jun 06 |
enell |
148 |
List<MASpot> fit = null; |
114 |
16 Jun 06 |
enell |
149 |
if (mode == ReporterMode.all) |
114 |
16 Jun 06 |
enell |
150 |
{ |
114 |
16 Jun 06 |
enell |
151 |
fit = lowess(group); |
114 |
16 Jun 06 |
enell |
152 |
} |
114 |
16 Jun 06 |
enell |
153 |
else |
114 |
16 Jun 06 |
enell |
154 |
{ |
114 |
16 Jun 06 |
enell |
155 |
List<Spot> normList = new ArrayList<Spot>(); |
114 |
16 Jun 06 |
enell |
156 |
if (mode == ReporterMode.include) |
114 |
16 Jun 06 |
enell |
157 |
{ |
114 |
16 Jun 06 |
enell |
158 |
for (Spot s : group) |
114 |
16 Jun 06 |
enell |
159 |
{ |
114 |
16 Jun 06 |
enell |
160 |
if (reporters.contains(s.getReporter().getReporterId())) |
114 |
16 Jun 06 |
enell |
161 |
{ |
114 |
16 Jun 06 |
enell |
162 |
normList.add(s); |
114 |
16 Jun 06 |
enell |
163 |
} |
114 |
16 Jun 06 |
enell |
164 |
} |
114 |
16 Jun 06 |
enell |
165 |
} |
114 |
16 Jun 06 |
enell |
166 |
else if (mode == ReporterMode.exclude) |
114 |
16 Jun 06 |
enell |
167 |
{ |
114 |
16 Jun 06 |
enell |
168 |
for (Spot s : group) |
114 |
16 Jun 06 |
enell |
169 |
{ |
114 |
16 Jun 06 |
enell |
170 |
if (!reporters.contains(s.getReporter().getReporterId())) |
114 |
16 Jun 06 |
enell |
171 |
{ |
114 |
16 Jun 06 |
enell |
172 |
normList.add(s); |
114 |
16 Jun 06 |
enell |
173 |
} |
114 |
16 Jun 06 |
enell |
174 |
} |
114 |
16 Jun 06 |
enell |
175 |
} |
114 |
16 Jun 06 |
enell |
176 |
if (normList.size() == 0) |
114 |
16 Jun 06 |
enell |
177 |
{ |
130 |
09 Aug 06 |
enell |
178 |
throw new BASEFileException( |
130 |
09 Aug 06 |
enell |
179 |
"The number of reporters to normalize in is to small. Change the include/exclude parameter."); |
114 |
16 Jun 06 |
enell |
180 |
} |
114 |
16 Jun 06 |
enell |
181 |
fit = lowess(normList); |
114 |
16 Jun 06 |
enell |
182 |
} |
114 |
16 Jun 06 |
enell |
183 |
|
114 |
16 Jun 06 |
enell |
// Fit spot to closest fit point |
114 |
16 Jun 06 |
enell |
185 |
int j = 0; |
114 |
16 Jun 06 |
enell |
186 |
MASpot mas = fit.get(0); |
114 |
16 Jun 06 |
enell |
187 |
for (Spot s : group) |
114 |
16 Jun 06 |
enell |
188 |
{ |
114 |
16 Jun 06 |
enell |
189 |
if (j < fit.size() - 1) |
114 |
16 Jun 06 |
enell |
190 |
{ |
114 |
16 Jun 06 |
enell |
191 |
MASpot nextMas = fit.get(j + 1); |
114 |
16 Jun 06 |
enell |
192 |
if (Math.abs(s.getA() - nextMas.getA()) <= Math.abs(s.getA() - mas.getA())) |
114 |
16 Jun 06 |
enell |
193 |
{ |
114 |
16 Jun 06 |
enell |
194 |
j++; |
114 |
16 Jun 06 |
enell |
195 |
mas = fit.get(j); |
114 |
16 Jun 06 |
enell |
196 |
} |
114 |
16 Jun 06 |
enell |
197 |
} |
114 |
16 Jun 06 |
enell |
198 |
s.setM(s.getM() - mas.getM()); |
114 |
16 Jun 06 |
enell |
199 |
} |
114 |
16 Jun 06 |
enell |
200 |
} |
114 |
16 Jun 06 |
enell |
201 |
} |
114 |
16 Jun 06 |
enell |
202 |
|
114 |
16 Jun 06 |
enell |
203 |
/** |
130 |
09 Aug 06 |
enell |
* Computes a lowess smooth curve from a list of {@link Spot Spots}. The list must be sorted on A. |
114 |
16 Jun 06 |
enell |
205 |
* |
114 |
16 Jun 06 |
enell |
* @param spots |
130 |
09 Aug 06 |
enell |
* @return |
114 |
16 Jun 06 |
enell |
208 |
*/ |
114 |
16 Jun 06 |
enell |
209 |
private List<MASpot> lowess(List<Spot> spots) |
114 |
16 Jun 06 |
enell |
210 |
{ |
114 |
16 Jun 06 |
enell |
211 |
Float[] smoothCurve = new Float[spots.size()]; |
114 |
16 Jun 06 |
enell |
212 |
int windowSize = Math.min(spots.size(), (int) (spots.size() * fitFraction + 0.5)); |
114 |
16 Jun 06 |
enell |
213 |
|
114 |
16 Jun 06 |
enell |
214 |
List<Double> wFit = new ArrayList<Double>(); |
114 |
16 Jun 06 |
enell |
215 |
wFit.addAll(Collections.nCopies(spots.size(), 1D)); |
130 |
09 Aug 06 |
enell |
216 |
for (int iteration = 0; iteration < iter; iteration++) |
114 |
16 Jun 06 |
enell |
217 |
{ |
114 |
16 Jun 06 |
enell |
218 |
int windowStart = 0; |
114 |
16 Jun 06 |
enell |
219 |
int i = 0; |
114 |
16 Jun 06 |
enell |
220 |
int prevI = -1; |
114 |
16 Jun 06 |
enell |
221 |
while (prevI < spots.size() - 1) |
114 |
16 Jun 06 |
enell |
222 |
{ |
114 |
16 Jun 06 |
enell |
// center window around the i:th value in x |
130 |
09 Aug 06 |
enell |
// while distance from windowStart to i is greater then distance |
130 |
09 Aug 06 |
enell |
// from i to windowEnd: move window |
114 |
16 Jun 06 |
enell |
226 |
while (windowStart + windowSize < spots.size() && spots.get(i).getA() - spots.get(windowStart).getA() > spots.get( |
114 |
16 Jun 06 |
enell |
227 |
windowStart + windowSize).getA() - spots.get(i).getA()) |
114 |
16 Jun 06 |
enell |
228 |
{ |
114 |
16 Jun 06 |
enell |
229 |
windowStart++; |
114 |
16 Jun 06 |
enell |
230 |
} |
114 |
16 Jun 06 |
enell |
231 |
|
114 |
16 Jun 06 |
enell |
232 |
List<Spot> window = spots.subList(windowStart, windowStart + windowSize); |
114 |
16 Jun 06 |
enell |
233 |
List<Double> wFitWindow = wFit.subList(windowStart, windowStart + windowSize); |
114 |
16 Jun 06 |
enell |
234 |
|
114 |
16 Jun 06 |
enell |
235 |
List<Double> w = calculateWeights(window, spots.get(i).getA(), wFitWindow); |
114 |
16 Jun 06 |
enell |
236 |
double[] km = weightedLeastSquaresRegression(window, w); |
114 |
16 Jun 06 |
enell |
237 |
smoothCurve[i] = (float) (km[0] * spots.get(i).getA() + km[1]); |
114 |
16 Jun 06 |
enell |
238 |
|
114 |
16 Jun 06 |
enell |
// Interpolate skipped points due to delta |
114 |
16 Jun 06 |
enell |
240 |
if (prevI + 1 < i) |
114 |
16 Jun 06 |
enell |
241 |
{ |
114 |
16 Jun 06 |
enell |
242 |
double d = spots.get(i).getA() - spots.get(prevI).getA(); |
114 |
16 Jun 06 |
enell |
243 |
if (d == 0) |
114 |
16 Jun 06 |
enell |
244 |
{ |
114 |
16 Jun 06 |
enell |
245 |
for (int j = prevI + 1; j < i; j++) |
114 |
16 Jun 06 |
enell |
246 |
{ |
114 |
16 Jun 06 |
enell |
247 |
smoothCurve[j] = smoothCurve[i]; |
114 |
16 Jun 06 |
enell |
248 |
} |
114 |
16 Jun 06 |
enell |
249 |
} |
114 |
16 Jun 06 |
enell |
250 |
else |
114 |
16 Jun 06 |
enell |
251 |
{ |
114 |
16 Jun 06 |
enell |
252 |
for (int j = prevI + 1; j < i; j++) |
114 |
16 Jun 06 |
enell |
253 |
{ |
114 |
16 Jun 06 |
enell |
254 |
double a = (spots.get(j).getA() - spots.get(prevI).getA()) / d; |
114 |
16 Jun 06 |
enell |
255 |
smoothCurve[j] = (float) (a * smoothCurve[i] + (1D - a) * smoothCurve[prevI]); |
114 |
16 Jun 06 |
enell |
256 |
} |
114 |
16 Jun 06 |
enell |
257 |
} |
114 |
16 Jun 06 |
enell |
258 |
} |
114 |
16 Jun 06 |
enell |
259 |
|
130 |
09 Aug 06 |
enell |
// increase i, next x value must be at least delta greater then |
130 |
09 Aug 06 |
enell |
// current |
114 |
16 Jun 06 |
enell |
262 |
prevI = i; |
114 |
16 Jun 06 |
enell |
263 |
double cut = spots.get(i).getA() + delta; |
114 |
16 Jun 06 |
enell |
264 |
while (i < spots.size() - 1 && spots.get(i).getA() <= cut) |
114 |
16 Jun 06 |
enell |
265 |
{ |
114 |
16 Jun 06 |
enell |
266 |
i++; |
114 |
16 Jun 06 |
enell |
267 |
} |
114 |
16 Jun 06 |
enell |
268 |
} |
114 |
16 Jun 06 |
enell |
269 |
|
114 |
16 Jun 06 |
enell |
270 |
double invYWRange = 1D / (ywRangeFactor * medianCorrection(spots, smoothCurve)); |
114 |
16 Jun 06 |
enell |
271 |
for (int j = 0; j < wFit.size(); j++) |
114 |
16 Jun 06 |
enell |
272 |
{ |
114 |
16 Jun 06 |
enell |
273 |
double w = Math.abs((smoothCurve[j] - spots.get(j).getM()) * invYWRange); |
114 |
16 Jun 06 |
enell |
274 |
wFit.set(j, w < 1 ? Math.pow(Math.pow(1 - w, 2), 2) : 0); |
114 |
16 Jun 06 |
enell |
275 |
} |
114 |
16 Jun 06 |
enell |
276 |
} |
114 |
16 Jun 06 |
enell |
277 |
|
114 |
16 Jun 06 |
enell |
278 |
ArrayList<MASpot> ret = new ArrayList<MASpot>(smoothCurve.length); |
114 |
16 Jun 06 |
enell |
279 |
for (int i = 0; i < spots.size(); i++) |
114 |
16 Jun 06 |
enell |
280 |
{ |
114 |
16 Jun 06 |
enell |
281 |
ret.add(new MASpot(smoothCurve[i], spots.get(i).getA())); |
114 |
16 Jun 06 |
enell |
282 |
} |
114 |
16 Jun 06 |
enell |
283 |
return ret; |
114 |
16 Jun 06 |
enell |
284 |
} |
114 |
16 Jun 06 |
enell |
285 |
|
114 |
16 Jun 06 |
enell |
286 |
private double medianCorrection(List<Spot> spots, Float[] smoothCurve) |
114 |
16 Jun 06 |
enell |
287 |
{ |
114 |
16 Jun 06 |
enell |
288 |
List<Float> temp = new ArrayList<Float>(spots.size()); |
114 |
16 Jun 06 |
enell |
289 |
for (int i = 0; i < spots.size(); i++) |
114 |
16 Jun 06 |
enell |
290 |
{ |
114 |
16 Jun 06 |
enell |
291 |
temp.add(Math.abs(spots.get(i).getM() - smoothCurve[i])); |
114 |
16 Jun 06 |
enell |
292 |
} |
114 |
16 Jun 06 |
enell |
293 |
Collections.sort(temp); |
114 |
16 Jun 06 |
enell |
294 |
if (temp.size() % 2 == 0) |
114 |
16 Jun 06 |
enell |
295 |
{ |
114 |
16 Jun 06 |
enell |
296 |
double v1 = temp.get((temp.size() / 2) - 1); |
114 |
16 Jun 06 |
enell |
297 |
double v2 = temp.get((temp.size() / 2)); |
114 |
16 Jun 06 |
enell |
298 |
return 0.5 * v1 * v2; |
114 |
16 Jun 06 |
enell |
299 |
} |
114 |
16 Jun 06 |
enell |
300 |
else |
114 |
16 Jun 06 |
enell |
301 |
{ |
114 |
16 Jun 06 |
enell |
302 |
return temp.get((temp.size() - 1) / 2); |
114 |
16 Jun 06 |
enell |
303 |
} |
114 |
16 Jun 06 |
enell |
304 |
} |
114 |
16 Jun 06 |
enell |
305 |
|
114 |
16 Jun 06 |
enell |
306 |
private double[] weightedLeastSquaresRegression(List<Spot> spots, List<Double> w) |
114 |
16 Jun 06 |
enell |
307 |
{ |
114 |
16 Jun 06 |
enell |
308 |
double k; |
114 |
16 Jun 06 |
enell |
309 |
double m; |
114 |
16 Jun 06 |
enell |
310 |
double sumX = 0; |
114 |
16 Jun 06 |
enell |
311 |
double sumY = 0; |
114 |
16 Jun 06 |
enell |
312 |
double sumXX = 0; |
114 |
16 Jun 06 |
enell |
313 |
double sumXY = 0; |
114 |
16 Jun 06 |
enell |
314 |
double sumW = 0; |
114 |
16 Jun 06 |
enell |
315 |
for (int j = 0; j < spots.size(); j++) |
114 |
16 Jun 06 |
enell |
316 |
{ |
114 |
16 Jun 06 |
enell |
317 |
double localW = w.get(j); |
114 |
16 Jun 06 |
enell |
318 |
sumX += spots.get(j).getA() * localW; |
114 |
16 Jun 06 |
enell |
319 |
sumY += spots.get(j).getM() * localW; |
114 |
16 Jun 06 |
enell |
320 |
sumXX += spots.get(j).getA() * spots.get(j).getA() * localW; |
114 |
16 Jun 06 |
enell |
321 |
sumXY += spots.get(j).getA() * spots.get(j).getM() * localW; |
114 |
16 Jun 06 |
enell |
322 |
sumW += localW; |
114 |
16 Jun 06 |
enell |
323 |
} |
114 |
16 Jun 06 |
enell |
324 |
if (sumW <= 0) |
114 |
16 Jun 06 |
enell |
325 |
{ |
114 |
16 Jun 06 |
enell |
326 |
System.out.println("dont panic"); |
114 |
16 Jun 06 |
enell |
327 |
} |
114 |
16 Jun 06 |
enell |
328 |
double denom = sumW * sumXX - sumX * sumX; |
114 |
16 Jun 06 |
enell |
329 |
if (denom != 0.) |
114 |
16 Jun 06 |
enell |
330 |
{ |
114 |
16 Jun 06 |
enell |
331 |
k = (sumW * sumXY - sumX * sumY) / denom; |
114 |
16 Jun 06 |
enell |
332 |
m = (sumY - k * sumX) / sumW; |
114 |
16 Jun 06 |
enell |
333 |
} |
114 |
16 Jun 06 |
enell |
334 |
else |
114 |
16 Jun 06 |
enell |
335 |
{ |
114 |
16 Jun 06 |
enell |
336 |
k = 0.; |
114 |
16 Jun 06 |
enell |
337 |
m = sumY / sumW; |
114 |
16 Jun 06 |
enell |
338 |
} |
114 |
16 Jun 06 |
enell |
339 |
|
114 |
16 Jun 06 |
enell |
340 |
return new double[] {k, m}; |
114 |
16 Jun 06 |
enell |
341 |
} |
114 |
16 Jun 06 |
enell |
342 |
|
114 |
16 Jun 06 |
enell |
343 |
private List<Double> calculateWeights(List<Spot> values, double x1, List<Double> wFit) |
114 |
16 Jun 06 |
enell |
344 |
{ |
114 |
16 Jun 06 |
enell |
345 |
List<Double> w = new ArrayList<Double>(); |
114 |
16 Jun 06 |
enell |
346 |
double invRadius = 0; |
114 |
16 Jun 06 |
enell |
347 |
for (Spot s : values) |
114 |
16 Jun 06 |
enell |
348 |
{ |
114 |
16 Jun 06 |
enell |
349 |
if (Math.abs(x1 - s.getA()) > invRadius) |
114 |
16 Jun 06 |
enell |
350 |
{ |
114 |
16 Jun 06 |
enell |
351 |
invRadius = Math.abs(x1 - s.getA()); |
114 |
16 Jun 06 |
enell |
352 |
} |
114 |
16 Jun 06 |
enell |
353 |
} |
114 |
16 Jun 06 |
enell |
354 |
invRadius = 1 / invRadius; |
114 |
16 Jun 06 |
enell |
355 |
for (int i = 0; i < values.size(); i++) |
114 |
16 Jun 06 |
enell |
356 |
{ |
114 |
16 Jun 06 |
enell |
357 |
Spot s = values.get(i); |
114 |
16 Jun 06 |
enell |
358 |
double distance = Math.abs(x1 - s.getA()) * invRadius; |
114 |
16 Jun 06 |
enell |
359 |
w.add((distance < 1 ? Math.pow(1D - Math.pow(Math.abs(distance), 3), 3) : 0) * wFit.get(i)); |
114 |
16 Jun 06 |
enell |
360 |
} |
114 |
16 Jun 06 |
enell |
361 |
return w; |
114 |
16 Jun 06 |
enell |
362 |
} |
114 |
16 Jun 06 |
enell |
363 |
|
114 |
16 Jun 06 |
enell |
364 |
private class GroupIterator implements Iterator<List<Spot>> |
114 |
16 Jun 06 |
enell |
365 |
{ |
114 |
16 Jun 06 |
enell |
366 |
BASEFileSpotSection<Reporter, Spot> bfss; |
114 |
16 Jun 06 |
enell |
367 |
|
114 |
16 Jun 06 |
enell |
368 |
private int start = 0; |
114 |
16 Jun 06 |
enell |
369 |
|
114 |
16 Jun 06 |
enell |
370 |
private int end = 0; |
114 |
16 Jun 06 |
enell |
371 |
|
114 |
16 Jun 06 |
enell |
372 |
private int currentGroup = -1; |
114 |
16 Jun 06 |
enell |
373 |
|
114 |
16 Jun 06 |
enell |
374 |
public GroupIterator(BASEFileSpotSection<Reporter, Spot> bfss) |
114 |
16 Jun 06 |
enell |
375 |
{ |
114 |
16 Jun 06 |
enell |
376 |
this.bfss = bfss; |
114 |
16 Jun 06 |
enell |
377 |
} |
114 |
16 Jun 06 |
enell |
378 |
|
114 |
16 Jun 06 |
enell |
379 |
public boolean hasNext() |
114 |
16 Jun 06 |
enell |
380 |
{ |
114 |
16 Jun 06 |
enell |
381 |
start = end; |
114 |
16 Jun 06 |
enell |
382 |
if (start >= bfss.getHeight()) |
114 |
16 Jun 06 |
enell |
383 |
{ |
114 |
16 Jun 06 |
enell |
384 |
return false; |
114 |
16 Jun 06 |
enell |
385 |
} |
114 |
16 Jun 06 |
enell |
386 |
currentGroup = bfss.getSpot(start).getBlock(); |
114 |
16 Jun 06 |
enell |
387 |
while (end < bfss.getSpotSize()) |
114 |
16 Jun 06 |
enell |
388 |
{ |
114 |
16 Jun 06 |
enell |
389 |
Spot s = bfss.getSpot(end); |
114 |
16 Jun 06 |
enell |
390 |
if (s.getBlock() != currentGroup) |
114 |
16 Jun 06 |
enell |
391 |
{ |
114 |
16 Jun 06 |
enell |
392 |
break; |
114 |
16 Jun 06 |
enell |
393 |
} |
114 |
16 Jun 06 |
enell |
394 |
end++; |
114 |
16 Jun 06 |
enell |
395 |
} |
114 |
16 Jun 06 |
enell |
396 |
return true; |
114 |
16 Jun 06 |
enell |
397 |
} |
114 |
16 Jun 06 |
enell |
398 |
|
114 |
16 Jun 06 |
enell |
399 |
public List<Spot> next() |
114 |
16 Jun 06 |
enell |
400 |
{ |
114 |
16 Jun 06 |
enell |
401 |
if (start >= bfss.getHeight()) |
114 |
16 Jun 06 |
enell |
402 |
{ |
114 |
16 Jun 06 |
enell |
403 |
return null; |
114 |
16 Jun 06 |
enell |
404 |
} |
114 |
16 Jun 06 |
enell |
405 |
return bfss.spotSubList(start, end); |
114 |
16 Jun 06 |
enell |
406 |
} |
114 |
16 Jun 06 |
enell |
407 |
|
114 |
16 Jun 06 |
enell |
408 |
public void remove() |
114 |
16 Jun 06 |
enell |
409 |
{ |
114 |
16 Jun 06 |
enell |
410 |
throw new UnsupportedOperationException(); |
114 |
16 Jun 06 |
enell |
411 |
} |
114 |
16 Jun 06 |
enell |
412 |
} |
114 |
16 Jun 06 |
enell |
413 |
|
114 |
16 Jun 06 |
enell |
414 |
private enum ReporterMode |
114 |
16 Jun 06 |
enell |
415 |
{ |
114 |
16 Jun 06 |
enell |
416 |
include, exclude, all; |
114 |
16 Jun 06 |
enell |
417 |
|
114 |
16 Jun 06 |
enell |
418 |
public static ReporterMode get(String mode) |
114 |
16 Jun 06 |
enell |
419 |
{ |
114 |
16 Jun 06 |
enell |
420 |
if (mode.equals("include")) |
114 |
16 Jun 06 |
enell |
421 |
{ |
114 |
16 Jun 06 |
enell |
422 |
return ReporterMode.include; |
114 |
16 Jun 06 |
enell |
423 |
} |
114 |
16 Jun 06 |
enell |
424 |
else if (mode.equals("exclude")) |
114 |
16 Jun 06 |
enell |
425 |
{ |
114 |
16 Jun 06 |
enell |
426 |
return ReporterMode.exclude; |
114 |
16 Jun 06 |
enell |
427 |
} |
114 |
16 Jun 06 |
enell |
428 |
else if (mode.equals("all")) |
114 |
16 Jun 06 |
enell |
429 |
{ |
114 |
16 Jun 06 |
enell |
430 |
return ReporterMode.all; |
114 |
16 Jun 06 |
enell |
431 |
} |
114 |
16 Jun 06 |
enell |
432 |
return null; |
114 |
16 Jun 06 |
enell |
433 |
} |
114 |
16 Jun 06 |
enell |
434 |
} |
114 |
16 Jun 06 |
enell |
435 |
} |