src/plugins/core/net/sf/basedb/plugins/LowessNormalization.java

Code
Comments
Other
Rev Date Author Line
2036 21 Feb 06 enell 1 /*
2036 21 Feb 06 enell 2   $Id$
2036 21 Feb 06 enell 3   
4889 06 Apr 09 nicklas 4   Copyright (C) 2006 Johan Enell, Jari Häkkinen, Nicklas Nordborg, Gregory Vincic
4889 06 Apr 09 nicklas 5   Copyright (C) 2007, 2008 Jari Häkkinen, Nicklas Nordborg
2036 21 Feb 06 enell 6   
2304 22 May 06 jari 7   This file is part of BASE - BioArray Software Environment.
2304 22 May 06 jari 8   Available at http://base.thep.lu.se/
2036 21 Feb 06 enell 9
2036 21 Feb 06 enell 10   BASE is free software; you can redistribute it and/or
2036 21 Feb 06 enell 11   modify it under the terms of the GNU General Public License
4480 05 Sep 08 jari 12   as published by the Free Software Foundation; either version 3
2036 21 Feb 06 enell 13   of the License, or (at your option) any later version.
2036 21 Feb 06 enell 14
2036 21 Feb 06 enell 15   BASE is distributed in the hope that it will be useful,
2036 21 Feb 06 enell 16   but WITHOUT ANY WARRANTY; without even the implied warranty of
2036 21 Feb 06 enell 17   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
2036 21 Feb 06 enell 18   GNU General Public License for more details.
2036 21 Feb 06 enell 19
2036 21 Feb 06 enell 20   You should have received a copy of the GNU General Public License
4513 11 Sep 08 jari 21   along with BASE. If not, see <http://www.gnu.org/licenses/>.
2036 21 Feb 06 enell 22 */
3410 30 May 07 jari 23
2036 21 Feb 06 enell 24 package net.sf.basedb.plugins;
2036 21 Feb 06 enell 25
2036 21 Feb 06 enell 26 import net.sf.basedb.core.BaseException;
2036 21 Feb 06 enell 27 import net.sf.basedb.core.BioAssay;
2036 21 Feb 06 enell 28 import net.sf.basedb.core.BioAssaySet;
2036 21 Feb 06 enell 29 import net.sf.basedb.core.DbControl;
2036 21 Feb 06 enell 30 import net.sf.basedb.core.DynamicResultIterator;
2036 21 Feb 06 enell 31 import net.sf.basedb.core.DynamicSpotQuery;
2036 21 Feb 06 enell 32 import net.sf.basedb.core.FloatParameterType;
4068 17 Dec 07 nicklas 33 import net.sf.basedb.core.Include;
2036 21 Feb 06 enell 34 import net.sf.basedb.core.IntegerParameterType;
6372 06 Dec 13 nicklas 35 import net.sf.basedb.core.IntensityTransform;
2036 21 Feb 06 enell 36 import net.sf.basedb.core.Item;
4068 17 Dec 07 nicklas 37 import net.sf.basedb.core.ItemParameterType;
4068 17 Dec 07 nicklas 38 import net.sf.basedb.core.ItemQuery;
2036 21 Feb 06 enell 39 import net.sf.basedb.core.Job;
2722 11 Oct 06 nicklas 40 import net.sf.basedb.core.Permission;
2036 21 Feb 06 enell 41 import net.sf.basedb.core.PluginParameter;
2036 21 Feb 06 enell 42 import net.sf.basedb.core.ProgressReporter;
4068 17 Dec 07 nicklas 43 import net.sf.basedb.core.ReporterList;
2036 21 Feb 06 enell 44 import net.sf.basedb.core.RequestInformation;
2036 21 Feb 06 enell 45 import net.sf.basedb.core.SpotBatcher;
2036 21 Feb 06 enell 46 import net.sf.basedb.core.Transformation;
3877 25 Oct 07 nicklas 47 import net.sf.basedb.core.Type;
2036 21 Feb 06 enell 48 import net.sf.basedb.core.VirtualColumn;
3615 31 Jul 07 nicklas 49 import net.sf.basedb.core.plugin.AbstractAnalysisPlugin;
2036 21 Feb 06 enell 50 import net.sf.basedb.core.plugin.GuiContext;
2036 21 Feb 06 enell 51 import net.sf.basedb.core.plugin.InteractivePlugin;
2722 11 Oct 06 nicklas 52 import net.sf.basedb.core.plugin.Permissions;
2036 21 Feb 06 enell 53 import net.sf.basedb.core.plugin.Request;
2036 21 Feb 06 enell 54 import net.sf.basedb.core.plugin.Response;
2036 21 Feb 06 enell 55 import net.sf.basedb.core.query.Dynamic;
3877 25 Oct 07 nicklas 56 import net.sf.basedb.core.query.Expression;
2036 21 Feb 06 enell 57 import net.sf.basedb.core.query.Expressions;
4068 17 Dec 07 nicklas 58 import net.sf.basedb.core.query.Hql;
4068 17 Dec 07 nicklas 59 import net.sf.basedb.core.query.JoinType;
2036 21 Feb 06 enell 60 import net.sf.basedb.core.query.Orders;
2036 21 Feb 06 enell 61 import net.sf.basedb.core.query.Restriction;
2036 21 Feb 06 enell 62 import net.sf.basedb.core.query.Restrictions;
2036 21 Feb 06 enell 63 import net.sf.basedb.core.query.Selects;
2036 21 Feb 06 enell 64 import net.sf.basedb.core.query.SqlResult;
4068 17 Dec 07 nicklas 65 import net.sf.basedb.core.query.WhenStatement;
5405 10 Sep 10 nicklas 66 import net.sf.basedb.core.signal.EnhancedThreadSignalHandler;
5405 10 Sep 10 nicklas 67 import net.sf.basedb.core.signal.Signal;
4811 13 Mar 09 nicklas 68 import net.sf.basedb.core.signal.SignalException;
4078 14 Jan 08 nicklas 69 import net.sf.basedb.core.signal.SignalHandler;
4078 14 Jan 08 nicklas 70 import net.sf.basedb.core.signal.SignalTarget;
4118 01 Feb 08 nicklas 71 import net.sf.basedb.core.signal.ThreadSignalHandler;
3615 31 Jul 07 nicklas 72 import net.sf.basedb.util.Values;
6372 06 Dec 13 nicklas 73 import net.sf.basedb.util.ma.MACalculator;
2036 21 Feb 06 enell 74
3877 25 Oct 07 nicklas 75 import java.sql.SQLException;
2036 21 Feb 06 enell 76 import java.util.ArrayList;
2036 21 Feb 06 enell 77 import java.util.Arrays;
2722 11 Oct 06 nicklas 78 import java.util.Collection;
2036 21 Feb 06 enell 79 import java.util.Collections;
2722 11 Oct 06 nicklas 80 import java.util.EnumSet;
2722 11 Oct 06 nicklas 81 import java.util.HashSet;
2036 21 Feb 06 enell 82 import java.util.List;
2036 21 Feb 06 enell 83 import java.util.Set;
4813 16 Mar 09 nicklas 84 import java.util.concurrent.Callable;
4813 16 Mar 09 nicklas 85 import java.util.concurrent.CompletionService;
4813 16 Mar 09 nicklas 86 import java.util.concurrent.ExecutionException;
4813 16 Mar 09 nicklas 87 import java.util.concurrent.ExecutorCompletionService;
4813 16 Mar 09 nicklas 88 import java.util.concurrent.ExecutorService;
4813 16 Mar 09 nicklas 89 import java.util.concurrent.Executors;
4813 16 Mar 09 nicklas 90 import java.util.concurrent.Future;
4813 16 Mar 09 nicklas 91 import java.util.concurrent.ThreadFactory;
2036 21 Feb 06 enell 92
3410 30 May 07 jari 93 /**
3877 25 Oct 07 nicklas 94    @author enell, Nicklas
3410 30 May 07 jari 95    @version 2.0
3410 30 May 07 jari 96    @base.modified $Date$
3410 30 May 07 jari 97 */
3410 30 May 07 jari 98
2036 21 Feb 06 enell 99 public class LowessNormalization
3615 31 Jul 07 nicklas 100   extends AbstractAnalysisPlugin
4078 14 Jan 08 nicklas 101   implements InteractivePlugin, SignalTarget
2036 21 Feb 06 enell 102 {
2036 21 Feb 06 enell 103
2036 21 Feb 06 enell 104   
2722 11 Oct 06 nicklas 105   private static final Set<Permissions> permissions = new HashSet<Permissions>();
2722 11 Oct 06 nicklas 106   
2036 21 Feb 06 enell 107   private static final PluginParameter<Float> fitFractionParameter = new PluginParameter<Float>
2036 21 Feb 06 enell 108   (
2036 21 Feb 06 enell 109     "fitFraction",
2036 21 Feb 06 enell 110     "Window size (fraction of points)",
2036 21 Feb 06 enell 111     "Smoothing parameter. The higher it is the smoother the funtion will be",
2062 23 Feb 06 gregory 112     new FloatParameterType(0F, 1F, 0.33F, true)
2036 21 Feb 06 enell 113   );
2036 21 Feb 06 enell 114   
2036 21 Feb 06 enell 115   private static final PluginParameter<Float> deltaParameter = new PluginParameter<Float>
2036 21 Feb 06 enell 116   (
2036 21 Feb 06 enell 117     "delta",
4371 03 Jul 08 nicklas 118     "Minimum log10(intensity) step",
2036 21 Feb 06 enell 119     "",
4371 03 Jul 08 nicklas 120     new FloatParameterType(0F, null, 0.043F, true)
2036 21 Feb 06 enell 121   );
2036 21 Feb 06 enell 122   
2036 21 Feb 06 enell 123   private static final PluginParameter<Integer> iterParameter = new PluginParameter<Integer>
2036 21 Feb 06 enell 124   (
2036 21 Feb 06 enell 125     "iter",
2036 21 Feb 06 enell 126     "Iterations",
2036 21 Feb 06 enell 127     "",
2036 21 Feb 06 enell 128     new IntegerParameterType(1, null, 4, true)
2036 21 Feb 06 enell 129   );
2036 21 Feb 06 enell 130
2132 30 Mar 06 enell 131   private static final PluginParameter<Integer> blockGroupParameter = new PluginParameter<Integer>
2132 30 Mar 06 enell 132   (
2132 30 Mar 06 enell 133     "blockGroup",
2132 30 Mar 06 enell 134     "Blockgroup size",
4061 14 Dec 07 nicklas 135     "The number of arraydesign blocks in each group.\n0 = use all blocks as a single group.",
4061 14 Dec 07 nicklas 136     new IntegerParameterType(0, null, 0, true)
2132 30 Mar 06 enell 137   );
3615 31 Jul 07 nicklas 138   
3615 31 Jul 07 nicklas 139   /**
3615 31 Jul 07 nicklas 140     Section for lowess options
3615 31 Jul 07 nicklas 141   */
3615 31 Jul 07 nicklas 142   protected static final PluginParameter<String> lowessSection = new PluginParameter<String>(
3615 31 Jul 07 nicklas 143     "lowessSection",
3615 31 Jul 07 nicklas 144     "Lowess options",
3615 31 Jul 07 nicklas 145     "", 
3615 31 Jul 07 nicklas 146     null
3615 31 Jul 07 nicklas 147   );
2036 21 Feb 06 enell 148
2036 21 Feb 06 enell 149   private static double ywRangeFactor = 6;
2036 21 Feb 06 enell 150   
2036 21 Feb 06 enell 151   private RequestInformation configureJob;
2036 21 Feb 06 enell 152
5405 10 Sep 10 nicklas 153   private EnhancedThreadSignalHandler signalHandler;
4078 14 Jan 08 nicklas 154   
2036 21 Feb 06 enell 155   /*
2036 21 Feb 06 enell 156     From the Plugin interface
2036 21 Feb 06 enell 157     -------------------------------------------
2036 21 Feb 06 enell 158   */
6127 14 Sep 12 nicklas 159   @Override
2512 10 Aug 06 nicklas 160   public boolean supportsConfigurations()
2512 10 Aug 06 nicklas 161   {
2512 10 Aug 06 nicklas 162     return false;
2512 10 Aug 06 nicklas 163   }
6127 14 Sep 12 nicklas 164   @Override
2512 10 Aug 06 nicklas 165   public boolean requiresConfiguration()
2512 10 Aug 06 nicklas 166   {
2512 10 Aug 06 nicklas 167     return false;
2512 10 Aug 06 nicklas 168   }
2722 11 Oct 06 nicklas 169   /**
2722 11 Oct 06 nicklas 170      Request use access to Experiment:s and read access to Job:s.
3410 30 May 07 jari 171     @since 2.0.2
2722 11 Oct 06 nicklas 172   */
6127 14 Sep 12 nicklas 173   @Override
2722 11 Oct 06 nicklas 174   public Collection<Permissions> getPermissions()
2722 11 Oct 06 nicklas 175   {
2722 11 Oct 06 nicklas 176     if (permissions.size() == 0)
2722 11 Oct 06 nicklas 177     {
2722 11 Oct 06 nicklas 178       permissions.add(new Permissions(Item.EXPERIMENT, null, EnumSet.of(Permission.USE)));
2722 11 Oct 06 nicklas 179       permissions.add(new Permissions(Item.JOB, null, EnumSet.of(Permission.READ)));
4068 17 Dec 07 nicklas 180       permissions.add(new Permissions(Item.REPORTERLIST, null, EnumSet.of(Permission.READ)));
2722 11 Oct 06 nicklas 181     }
2722 11 Oct 06 nicklas 182     return permissions;
2722 11 Oct 06 nicklas 183   }
2036 21 Feb 06 enell 184
6127 14 Sep 12 nicklas 185   @Override
2036 21 Feb 06 enell 186   public void run(Request request, Response response, ProgressReporter progress)
2036 21 Feb 06 enell 187   {
2036 21 Feb 06 enell 188     String command = request.getCommand();
2036 21 Feb 06 enell 189     if (command.equals(Request.COMMAND_EXECUTE))
2036 21 Feb 06 enell 190     {
5405 10 Sep 10 nicklas 191       if (signalHandler != null) signalHandler.setWorkerThread();
3615 31 Jul 07 nicklas 192       DbControl dc = null;
2036 21 Feb 06 enell 193       try
2036 21 Feb 06 enell 194       {
3615 31 Jul 07 nicklas 195         dc = sc.newDbControl();
3615 31 Jul 07 nicklas 196         BioAssaySet source = getSourceBioAssaySet(dc);
3615 31 Jul 07 nicklas 197         String childName = Values.getString((String)job.getValue(CHILD_NAME), source.getName());
3615 31 Jul 07 nicklas 198         String childDescription = (String)job.getValue(CHILD_DESCRIPTION);
3615 31 Jul 07 nicklas 199         int blockGroupSize = (Integer) job.getValue(blockGroupParameter.getName());
2036 21 Feb 06 enell 200         float delta = (Float) job.getValue(deltaParameter.getName());
2036 21 Feb 06 enell 201         float fitFraction = (Float) job.getValue(fitFractionParameter.getName());
3877 25 Oct 07 nicklas 202         int iterations = (Integer) job.getValue(iterParameter.getName());
4068 17 Dec 07 nicklas 203         ReporterList excludeReporters = (ReporterList)job.getValue("excludeReporters");
3877 25 Oct 07 nicklas 204         Job thisJob = getCurrentJob(dc);
2133 30 Mar 06 enell 205         
4068 17 Dec 07 nicklas 206         BioAssaySet child = normalize(dc, source, thisJob, fitFraction, delta, iterations, 
4068 17 Dec 07 nicklas 207           blockGroupSize, excludeReporters, progress);
3877 25 Oct 07 nicklas 208         child.setName(childName);
3877 25 Oct 07 nicklas 209         child.setDescription(childDescription);
2040 22 Feb 06 nicklas 210         dc.commit();
3877 25 Oct 07 nicklas 211         int normalizedSpots = child.getNumSpots();
2285 18 May 06 nicklas 212         if (progress != null) progress.display(100, normalizedSpots + " spots normalized\n");
3615 31 Jul 07 nicklas 213         response.setDone(normalizedSpots + " spots normalized, " + (source.getNumSpots() - normalizedSpots) + " spots removed");
2036 21 Feb 06 enell 214       }
5405 10 Sep 10 nicklas 215       catch (SignalException ex)
5405 10 Sep 10 nicklas 216       {
5405 10 Sep 10 nicklas 217         if (signalHandler != null && signalHandler.hasReceived(Signal.SHUTDOWN))
5405 10 Sep 10 nicklas 218         {
5405 10 Sep 10 nicklas 219           response.setContinue(command);
5405 10 Sep 10 nicklas 220         }
5405 10 Sep 10 nicklas 221         else
5405 10 Sep 10 nicklas 222         {
5405 10 Sep 10 nicklas 223           response.setError("Aborted by user", Arrays.asList(ex));
5405 10 Sep 10 nicklas 224         }
5405 10 Sep 10 nicklas 225       }
2036 21 Feb 06 enell 226       catch (Throwable ex)
2036 21 Feb 06 enell 227       {
2036 21 Feb 06 enell 228         response.setError(ex.getMessage(), Arrays.asList(ex));
2036 21 Feb 06 enell 229       }
2036 21 Feb 06 enell 230       finally
2036 21 Feb 06 enell 231       {
2040 22 Feb 06 nicklas 232         if (dc != null) dc.close();
2036 21 Feb 06 enell 233       }
2036 21 Feb 06 enell 234     }
2036 21 Feb 06 enell 235     else
2036 21 Feb 06 enell 236     {
2036 21 Feb 06 enell 237       response.setError("Unknown command: " + command, null);
2036 21 Feb 06 enell 238     }
2036 21 Feb 06 enell 239   }
2036 21 Feb 06 enell 240   // -------------------------------------------
2036 21 Feb 06 enell 241
2036 21 Feb 06 enell 242   /*
2036 21 Feb 06 enell 243     From the InteractivePlugin interface
2036 21 Feb 06 enell 244     -------------------------------------------
2036 21 Feb 06 enell 245   */
6127 14 Sep 12 nicklas 246   @Override
2195 26 Apr 06 nicklas 247   public String isInContext(GuiContext context, Object item)
2036 21 Feb 06 enell 248   {
3615 31 Jul 07 nicklas 249     String message = super.isInContext(context, item);
3615 31 Jul 07 nicklas 250     if (message == null)
2036 21 Feb 06 enell 251     {
2195 26 Apr 06 nicklas 252       BioAssaySet bas = (BioAssaySet)item;
2195 26 Apr 06 nicklas 253       int channels = bas.getRawDataType().getChannels();
2195 26 Apr 06 nicklas 254       if (channels != 2)
2195 26 Apr 06 nicklas 255       {
2195 26 Apr 06 nicklas 256         message = "Lowess normalization requires 2-channel data, not " + channels + "-channel.";
2195 26 Apr 06 nicklas 257       }
4337 17 Jun 08 nicklas 258       else if (bas.getNumSpots() <= 0)
2993 01 Dec 06 nicklas 259       {
4337 17 Jun 08 nicklas 260         message = "Lowess normalization requires spot data in the database";
4337 17 Jun 08 nicklas 261       }
4337 17 Jun 08 nicklas 262       else if (bas.getMaxRawMappingsForSpot() != 1)
4337 17 Jun 08 nicklas 263       {
2993 01 Dec 06 nicklas 264         message = "Lowess normalization can't be done on merged data.";
2993 01 Dec 06 nicklas 265       }
2195 26 Apr 06 nicklas 266     }
2195 26 Apr 06 nicklas 267     return message;
2036 21 Feb 06 enell 268   }
2036 21 Feb 06 enell 269
6127 14 Sep 12 nicklas 270   @Override
2036 21 Feb 06 enell 271   public RequestInformation getRequestInformation(GuiContext context, String command) throws BaseException
2036 21 Feb 06 enell 272   {
2036 21 Feb 06 enell 273     RequestInformation requestInformation = null;
2512 10 Aug 06 nicklas 274     if (command.equals(Request.COMMAND_CONFIGURE_JOB))
2036 21 Feb 06 enell 275     {
3615 31 Jul 07 nicklas 276       requestInformation = getConfigureJobParameters();
2036 21 Feb 06 enell 277     }
2036 21 Feb 06 enell 278     return requestInformation;
2036 21 Feb 06 enell 279   }
2036 21 Feb 06 enell 280
6127 14 Sep 12 nicklas 281   @Override
2036 21 Feb 06 enell 282   public void configure(GuiContext context, Request request, Response response)
2036 21 Feb 06 enell 283   {
2036 21 Feb 06 enell 284     String command = request.getCommand();
2036 21 Feb 06 enell 285     try
2036 21 Feb 06 enell 286     {
2512 10 Aug 06 nicklas 287       if (command.equals(Request.COMMAND_CONFIGURE_JOB))
2036 21 Feb 06 enell 288       {
3615 31 Jul 07 nicklas 289         // Validate job parameters
3615 31 Jul 07 nicklas 290         RequestInformation ri = getConfigureJobParameters();
3615 31 Jul 07 nicklas 291         List<Throwable> errors = validateRequestParameters(ri.getParameters(), request);
2036 21 Feb 06 enell 292         if (errors != null)
2036 21 Feb 06 enell 293         {
3615 31 Jul 07 nicklas 294           response.setError(errors.size() + " invalid parameter(s) were found in the request", errors);
2036 21 Feb 06 enell 295           return;
2036 21 Feb 06 enell 296         }
3615 31 Jul 07 nicklas 297
3615 31 Jul 07 nicklas 298         // Source bioassay set
3615 31 Jul 07 nicklas 299         storeValue(job, request, ri.getParameter(SOURCE_BIOASSAYSET));
3615 31 Jul 07 nicklas 300                 
3615 31 Jul 07 nicklas 301         // Child name, description and transformation
3615 31 Jul 07 nicklas 302         storeValue(job, request, ri.getParameter(CHILD_NAME));
3615 31 Jul 07 nicklas 303         storeValue(job, request, ri.getParameter(CHILD_DESCRIPTION));
3615 31 Jul 07 nicklas 304
3615 31 Jul 07 nicklas 305         // Lowess options
2036 21 Feb 06 enell 306         storeValue(job, request, fitFractionParameter);
2036 21 Feb 06 enell 307         storeValue(job, request, deltaParameter);
2036 21 Feb 06 enell 308         storeValue(job, request, iterParameter);
2132 30 Mar 06 enell 309         storeValue(job, request, blockGroupParameter);
4068 17 Dec 07 nicklas 310         storeValue(job, request, ri.getParameter("excludeReporters"));
5481 08 Nov 10 nicklas 311         
5481 08 Nov 10 nicklas 312         BioAssaySet bioAssaySet = (BioAssaySet)job.getValue(SOURCE_BIOASSAYSET);
5481 08 Nov 10 nicklas 313         response.setSuggestedJobName("Lowess normalization on bioassay set '" + 
5481 08 Nov 10 nicklas 314             bioAssaySet.getName() + "'");
5481 08 Nov 10 nicklas 315         
2036 21 Feb 06 enell 316         response.setDone("Job configuration complete", Job.ExecutionTime.SHORT);
2036 21 Feb 06 enell 317       }
2036 21 Feb 06 enell 318     }
2036 21 Feb 06 enell 319     catch (Throwable ex)
2036 21 Feb 06 enell 320     {
2036 21 Feb 06 enell 321       response.setError(ex.getMessage(), Arrays.asList(ex));
2036 21 Feb 06 enell 322     }
2036 21 Feb 06 enell 323   }
2036 21 Feb 06 enell 324   // -------------------------------------------
2036 21 Feb 06 enell 325
4078 14 Jan 08 nicklas 326   /*
4078 14 Jan 08 nicklas 327     From the SignalTarget interface
4078 14 Jan 08 nicklas 328     -------------------------------------------
4078 14 Jan 08 nicklas 329   */
6127 14 Sep 12 nicklas 330   @Override
4078 14 Jan 08 nicklas 331   public SignalHandler getSignalHandler()
4078 14 Jan 08 nicklas 332   {
5405 10 Sep 10 nicklas 333     if (signalHandler == null) 
5405 10 Sep 10 nicklas 334     {
5405 10 Sep 10 nicklas 335       signalHandler = new EnhancedThreadSignalHandler(Signal.ABORT, Signal.SHUTDOWN);
5405 10 Sep 10 nicklas 336     }
4078 14 Jan 08 nicklas 337     return signalHandler;
4078 14 Jan 08 nicklas 338   }
4078 14 Jan 08 nicklas 339   // -------------------------------------------
4078 14 Jan 08 nicklas 340
3877 25 Oct 07 nicklas 341   /**
3877 25 Oct 07 nicklas 342     Normalise the source bioassay set using LOWESS normalization.
4068 17 Dec 07 nicklas 343     All reporters are included.
4068 17 Dec 07 nicklas 344     @since 2.5
4068 17 Dec 07 nicklas 345     @see #normalize(DbControl, BioAssaySet, Job, float, float, int, int, ReporterList, ProgressReporter)
4068 17 Dec 07 nicklas 346   */
4068 17 Dec 07 nicklas 347   public BioAssaySet normalize(DbControl dc, BioAssaySet source, Job job, float fitFraction, float delta, int iterations, int blockGroupSize, ProgressReporter progress)
4068 17 Dec 07 nicklas 348   {
4068 17 Dec 07 nicklas 349     return normalize(dc, source, job, fitFraction, delta, iterations, blockGroupSize, null, progress);
4068 17 Dec 07 nicklas 350   }
4068 17 Dec 07 nicklas 351   
4068 17 Dec 07 nicklas 352   
4068 17 Dec 07 nicklas 353   /**
4068 17 Dec 07 nicklas 354     Normalise the source bioassay set using LOWESS normalization.
3877 25 Oct 07 nicklas 355     @param dc The DbControl to use for database access
3877 25 Oct 07 nicklas 356     @param source The source bioassay set that is going to be normalized
3877 25 Oct 07 nicklas 357     @param job The job that is doing the normalization, or null
4061 14 Dec 07 nicklas 358     @param blockGroupSize The number of blocks to group and normalise at the same time,
6898 12 May 15 nicklas 359       if &lt;=0, all blocks on a bioassay are grouped to a single dataset
4068 17 Dec 07 nicklas 360     @param excludeReporters A reporter list containing reporters that should not be
4068 17 Dec 07 nicklas 361       used in when normalizing, or null to use all spots
4034 05 Dec 07 martin 362      @param progress Progress reporter used by the caller to keep track of the progress. Null is allowed   
3877 25 Oct 07 nicklas 363     @return The normalized bioassayset
4069 17 Dec 07 nicklas 364     @since 2.6
3877 25 Oct 07 nicklas 365   */
4068 17 Dec 07 nicklas 366   public BioAssaySet normalize(DbControl dc, BioAssaySet source, Job job, float fitFraction, float delta, int iterations, int blockGroupSize, ReporterList excludeReporters, ProgressReporter progress)
3877 25 Oct 07 nicklas 367   {
3877 25 Oct 07 nicklas 368     if (progress != null) progress.display(0, "Preparing to normalize...");
3877 25 Oct 07 nicklas 369     
3877 25 Oct 07 nicklas 370     // Create Transformation
3877 25 Oct 07 nicklas 371     Transformation t = source.newTransformation(job);
5612 18 Apr 11 nicklas 372     t.setName("Normalization: Lowess");
3877 25 Oct 07 nicklas 373     dc.saveItem(t);
3877 25 Oct 07 nicklas 374     
3877 25 Oct 07 nicklas 375     // Create the normalized bioassay set
3877 25 Oct 07 nicklas 376     BioAssaySet child = t.newProduct(null, "new", true);
6372 06 Dec 13 nicklas 377     IntensityTransform transform = source.getIntensityTransform();
6372 06 Dec 13 nicklas 378     child.setIntensityTransform(transform);
3877 25 Oct 07 nicklas 379     dc.saveItem(child);
3877 25 Oct 07 nicklas 380     
6372 06 Dec 13 nicklas 381     MACalculator mac = transform.getMACalculator();
6372 06 Dec 13 nicklas 382     
3877 25 Oct 07 nicklas 383     // Batcher for inserting normalized data
3877 25 Oct 07 nicklas 384     SpotBatcher batcher = child.getSpotBatcher();
3877 25 Oct 07 nicklas 385
3877 25 Oct 07 nicklas 386     // Expressions used to get data
3877 25 Oct 07 nicklas 387     Expression block = Dynamic.rawData("block");
6372 06 Dec 13 nicklas 388     Expression ch1 = Dynamic.column(VirtualColumn.channelRaw(1));
6372 06 Dec 13 nicklas 389     Expression ch2 = Dynamic.column(VirtualColumn.channelRaw(2));
3877 25 Oct 07 nicklas 390
3877 25 Oct 07 nicklas 391     // Create restriction: column = :bioAssayColumn
3877 25 Oct 07 nicklas 392     Restriction bioAssayRestriction = Restrictions.eq(
3877 25 Oct 07 nicklas 393         Dynamic.column(VirtualColumn.COLUMN), 
3877 25 Oct 07 nicklas 394         Expressions.parameter("bioAssayColumn")
3877 25 Oct 07 nicklas 395       );
3877 25 Oct 07 nicklas 396       
4813 16 Mar 09 nicklas 397     // Build query for loading spots
3877 25 Oct 07 nicklas 398     DynamicSpotQuery query = source.getSpotData();
6372 06 Dec 13 nicklas 399     // Create restriction: ch1 > 0 and ch2 > 0
6372 06 Dec 13 nicklas 400     if (!transform.isTransformed())
6372 06 Dec 13 nicklas 401     {
6372 06 Dec 13 nicklas 402       Restriction intensityRestriction = Restrictions.and(
6372 06 Dec 13 nicklas 403           Restrictions.gt(ch1, Expressions.aFloat(0)), 
6372 06 Dec 13 nicklas 404           Restrictions.gt(ch2, Expressions.aFloat(0))
6372 06 Dec 13 nicklas 405         );
6372 06 Dec 13 nicklas 406       query.restrictPermanent(intensityRestriction);
6372 06 Dec 13 nicklas 407     }
3877 25 Oct 07 nicklas 408     
4068 17 Dec 07 nicklas 409     // If some reporters should be excluded, create an expression that:
4068 17 Dec 07 nicklas 410     // = 0 if spot's reporter is not present in list, = 1 if it is
4068 17 Dec 07 nicklas 411     Expression isPresent = Expressions.integer(0);
4068 17 Dec 07 nicklas 412     if (excludeReporters != null)
4068 17 Dec 07 nicklas 413     {
4068 17 Dec 07 nicklas 414       query.joinReporterList(excludeReporters, JoinType.LEFT);
4068 17 Dec 07 nicklas 415       isPresent = Expressions.caseWhen(isPresent, 
4068 17 Dec 07 nicklas 416         new WhenStatement(Dynamic.isPartOf(excludeReporters), Expressions.integer(1)));
4068 17 Dec 07 nicklas 417     }
4068 17 Dec 07 nicklas 418     
4068 17 Dec 07 nicklas 419     // Create query to retrieve spot data: COLUMN, POSITION, ch1, ch2, block, present
3877 25 Oct 07 nicklas 420     // We use a parameter to restrict the query to return data for one bioassay at a time
4118 01 Feb 08 nicklas 421     // The 'present' value is 0 to not exclude and 1 to exclude a spot
3877 25 Oct 07 nicklas 422     query.select(Dynamic.select(VirtualColumn.POSITION));
3877 25 Oct 07 nicklas 423     query.select(Selects.expression(ch1, "ch1"));
3877 25 Oct 07 nicklas 424     query.select(Selects.expression(ch2, "ch2"));
3877 25 Oct 07 nicklas 425     query.select(Selects.expression(block, "block"));
4068 17 Dec 07 nicklas 426     query.select(Selects.expression(isPresent, "present"));
3877 25 Oct 07 nicklas 427     query.restrict(bioAssayRestriction);
3877 25 Oct 07 nicklas 428     query.order(Orders.asc(block));
4358 01 Jul 08 nicklas 429     
3877 25 Oct 07 nicklas 430     // Normalize one bioassay at a time
3877 25 Oct 07 nicklas 431     List<BioAssay> assays = source.getBioAssays().list(dc);
4811 13 Mar 09 nicklas 432     int numAssays = assays.size();
4813 16 Mar 09 nicklas 433     int nextAssayIndex = 0;
4813 16 Mar 09 nicklas 434     int maxThreads = Math.min(4, Runtime.getRuntime().availableProcessors());
4813 16 Mar 09 nicklas 435     int maxQueueSize = 2 * maxThreads;
4813 16 Mar 09 nicklas 436     int numDone = 0;
4813 16 Mar 09 nicklas 437     int numWorking = 0;
4811 13 Mar 09 nicklas 438     
4813 16 Mar 09 nicklas 439     // Set up execution service for thread management
4813 16 Mar 09 nicklas 440     final ThreadGroup workers = new ThreadGroup("Lowess-threads");
4813 16 Mar 09 nicklas 441     ThreadFactory threadFactory = new ThreadFactory()
3877 25 Oct 07 nicklas 442     {
4813 16 Mar 09 nicklas 443       private int i=0;
4813 16 Mar 09 nicklas 444       @Override
4813 16 Mar 09 nicklas 445       public Thread newThread(Runnable r) 
3877 25 Oct 07 nicklas 446       {
4813 16 Mar 09 nicklas 447         ++i;
4813 16 Mar 09 nicklas 448         return new Thread(workers, r, "Lowess-worker-"+i);
4811 13 Mar 09 nicklas 449       }
4813 16 Mar 09 nicklas 450     };
4813 16 Mar 09 nicklas 451     final ExecutorService service = Executors.newFixedThreadPool(maxThreads, threadFactory);
4813 16 Mar 09 nicklas 452     final CompletionService<List<SpotData>> cService = 
4813 16 Mar 09 nicklas 453       new ExecutorCompletionService<List<SpotData>>(service);
4813 16 Mar 09 nicklas 454     try
4813 16 Mar 09 nicklas 455     {
4813 16 Mar 09 nicklas 456       // Loop until all bioassays has been processed
4813 16 Mar 09 nicklas 457       while (numDone < numAssays)
4811 13 Mar 09 nicklas 458       {
5405 10 Sep 10 nicklas 459         ThreadSignalHandler.checkInterrupted();
4813 16 Mar 09 nicklas 460
4813 16 Mar 09 nicklas 461         // Get data for the next bioassay unless the working queue
4813 16 Mar 09 nicklas 462         // is full
4813 16 Mar 09 nicklas 463         if (nextAssayIndex < numAssays && numWorking < maxQueueSize)
3877 25 Oct 07 nicklas 464         {
6372 06 Dec 13 nicklas 465           List<SpotData> data = readData(dc, query, assays.get(nextAssayIndex), mac);
4813 16 Mar 09 nicklas 466           nextAssayIndex++;
4813 16 Mar 09 nicklas 467           if (data.size() > 0) 
4813 16 Mar 09 nicklas 468           {
6372 06 Dec 13 nicklas 469             CallableWorker worker = new CallableWorker(data, mac, fitFraction, delta, iterations, blockGroupSize);
4813 16 Mar 09 nicklas 470             cService.submit(worker);
4813 16 Mar 09 nicklas 471             numWorking++;
4813 16 Mar 09 nicklas 472           }
4813 16 Mar 09 nicklas 473           else
4813 16 Mar 09 nicklas 474           {
4813 16 Mar 09 nicklas 475             numDone++;
4813 16 Mar 09 nicklas 476           }
4813 16 Mar 09 nicklas 477         }
4813 16 Mar 09 nicklas 478         
4813 16 Mar 09 nicklas 479         // If the working queue is full we use take() to wait until a result is
4813 16 Mar 09 nicklas 480         // ready, otherwise we use poll() which let's us continue if no results
4813 16 Mar 09 nicklas 481         // are completed
4813 16 Mar 09 nicklas 482         Future<List<SpotData>> result = numWorking >= maxQueueSize ? 
4813 16 Mar 09 nicklas 483           cService.take() : cService.poll();
4813 16 Mar 09 nicklas 484         if (result != null)
4813 16 Mar 09 nicklas 485         {
4813 16 Mar 09 nicklas 486           numDone++;
4813 16 Mar 09 nicklas 487           numWorking--;
4813 16 Mar 09 nicklas 488           List<SpotData> data = result.get();
4811 13 Mar 09 nicklas 489           for (SpotData spot : data)
3877 25 Oct 07 nicklas 490           {
4811 13 Mar 09 nicklas 491             batcher.insert(spot.column, spot.position, spot.newCh1, spot.newCh2);
3877 25 Oct 07 nicklas 492           }
4813 16 Mar 09 nicklas 493           batcher.flush();
4811 13 Mar 09 nicklas 494           if (progress != null)
3877 25 Oct 07 nicklas 495           {
4813 16 Mar 09 nicklas 496             int percent = (100 * numDone) / numAssays;
4813 16 Mar 09 nicklas 497             progress.display(percent, numDone + " of " + numAssays + 
4813 16 Mar 09 nicklas 498               " assays done.");
3877 25 Oct 07 nicklas 499           }
3877 25 Oct 07 nicklas 500         }
3877 25 Oct 07 nicklas 501       }
3877 25 Oct 07 nicklas 502       batcher.close();
3877 25 Oct 07 nicklas 503     }
3877 25 Oct 07 nicklas 504     catch (SQLException e)
3877 25 Oct 07 nicklas 505     {
3877 25 Oct 07 nicklas 506       throw new BaseException(e);
3877 25 Oct 07 nicklas 507     }
4813 16 Mar 09 nicklas 508     catch (ExecutionException ex)
4813 16 Mar 09 nicklas 509     {
4813 16 Mar 09 nicklas 510       throw new BaseException(ex.getCause());
4813 16 Mar 09 nicklas 511     }
4811 13 Mar 09 nicklas 512     catch (InterruptedException ex)
4811 13 Mar 09 nicklas 513     {
5405 10 Sep 10 nicklas 514       if (signalHandler != null)
5405 10 Sep 10 nicklas 515       {
5405 10 Sep 10 nicklas 516         signalHandler.checkForSignals();
5405 10 Sep 10 nicklas 517       }
5405 10 Sep 10 nicklas 518       throw new SignalException("Aborted by user", ex);
4811 13 Mar 09 nicklas 519     }
4811 13 Mar 09 nicklas 520     finally
4811 13 Mar 09 nicklas 521     {
4813 16 Mar 09 nicklas 522       service.shutdownNow();
4811 13 Mar 09 nicklas 523     }
3877 25 Oct 07 nicklas 524     return child;
3877 25 Oct 07 nicklas 525   }
3877 25 Oct 07 nicklas 526   
6372 06 Dec 13 nicklas 527   private List<SpotData> readData(DbControl dc, DynamicSpotQuery query, BioAssay assay, MACalculator mac)
4811 13 Mar 09 nicklas 528     throws SQLException
4811 13 Mar 09 nicklas 529   {
4811 13 Mar 09 nicklas 530     // Prepare list for holding data
4811 13 Mar 09 nicklas 531     int assaySpots = assay.getNumSpots();
4811 13 Mar 09 nicklas 532     List<SpotData> data = new ArrayList<SpotData>(assaySpots);
4811 13 Mar 09 nicklas 533     
4811 13 Mar 09 nicklas 534     // Load spot data for this bioassay
4811 13 Mar 09 nicklas 535     short bioassayColumn = assay.getDataCubeColumnNo();
4811 13 Mar 09 nicklas 536     query.setParameter("bioAssayColumn", (int)bioassayColumn, Type.INT);
4811 13 Mar 09 nicklas 537     
4811 13 Mar 09 nicklas 538     DynamicResultIterator it = query.iterate(dc);
4811 13 Mar 09 nicklas 539     int positionIndex = it.getIndex(VirtualColumn.POSITION.getName());
4811 13 Mar 09 nicklas 540     int ch1Index = it.getIndex("ch1");
4811 13 Mar 09 nicklas 541     int ch2Index = it.getIndex("ch2");
4811 13 Mar 09 nicklas 542     int blockIndex = it.getIndex("block");
4811 13 Mar 09 nicklas 543     int presentIndex = it.getIndex("present");
4811 13 Mar 09 nicklas 544     
4811 13 Mar 09 nicklas 545     // Copy bioassay data to SpotData objects
4811 13 Mar 09 nicklas 546     while (it.hasNext())
4811 13 Mar 09 nicklas 547     {
5405 10 Sep 10 nicklas 548       ThreadSignalHandler.checkInterrupted();
4811 13 Mar 09 nicklas 549       SqlResult r = it.next();
4811 13 Mar 09 nicklas 550       SpotData spot = new SpotData(bioassayColumn, r.getInt(positionIndex),
4811 13 Mar 09 nicklas 551         r.getFloat(ch1Index), r.getFloat(ch2Index), r.getInt(blockIndex), 
6372 06 Dec 13 nicklas 552         r.getInt(presentIndex) == 1, mac);
4811 13 Mar 09 nicklas 553       data.add(spot);
4811 13 Mar 09 nicklas 554     }
4811 13 Mar 09 nicklas 555     it.close();
4811 13 Mar 09 nicklas 556     return data;
4811 13 Mar 09 nicklas 557   }
4811 13 Mar 09 nicklas 558
3615 31 Jul 07 nicklas 559   private RequestInformation getConfigureJobParameters()
3615 31 Jul 07 nicklas 560   {
3615 31 Jul 07 nicklas 561     DbControl dc = null;
3615 31 Jul 07 nicklas 562     try
3615 31 Jul 07 nicklas 563     {
3615 31 Jul 07 nicklas 564       if (configureJob == null)
3615 31 Jul 07 nicklas 565       {
3615 31 Jul 07 nicklas 566         dc = sc.newDbControl();
3615 31 Jul 07 nicklas 567         BioAssaySet bas = getCurrentBioAssaySet(dc);
3615 31 Jul 07 nicklas 568         List<PluginParameter<?>> parameters = new ArrayList<PluginParameter<?>>(3);
3615 31 Jul 07 nicklas 569         
3615 31 Jul 07 nicklas 570         // Source and child bioassay set parameters
3615 31 Jul 07 nicklas 571         parameters.add(getSourceBioAssaySetParameter(null, null));
3615 31 Jul 07 nicklas 572         parameters.add(getChildNameParameter(null, null, bas.getName()));
3615 31 Jul 07 nicklas 573         parameters.add(getChildDescriptionParameter(null, null, null));        
3615 31 Jul 07 nicklas 574         
3615 31 Jul 07 nicklas 575         // Lowess options
3615 31 Jul 07 nicklas 576         parameters.add(lowessSection);
3615 31 Jul 07 nicklas 577         parameters.add(fitFractionParameter);
3615 31 Jul 07 nicklas 578         parameters.add(deltaParameter);
3615 31 Jul 07 nicklas 579         parameters.add(iterParameter);
3615 31 Jul 07 nicklas 580         parameters.add(blockGroupParameter);
4068 17 Dec 07 nicklas 581         PluginParameter<ReporterList> excludeReporters = 
4068 17 Dec 07 nicklas 582           getExcludeReportersParameters(dc);
4068 17 Dec 07 nicklas 583         if (excludeReporters != null) parameters.add(excludeReporters);
3615 31 Jul 07 nicklas 584   
3615 31 Jul 07 nicklas 585         configureJob = new RequestInformation
3615 31 Jul 07 nicklas 586         (
3615 31 Jul 07 nicklas 587           Request.COMMAND_CONFIGURE_JOB,
3615 31 Jul 07 nicklas 588           "Configure job",
3615 31 Jul 07 nicklas 589           "Set window size, minimum step and number of iterations",
3615 31 Jul 07 nicklas 590           parameters
3615 31 Jul 07 nicklas 591         );
3615 31 Jul 07 nicklas 592       }
3615 31 Jul 07 nicklas 593     }
3615 31 Jul 07 nicklas 594     finally
3615 31 Jul 07 nicklas 595     {
3615 31 Jul 07 nicklas 596       if (dc != null) dc.close();
3615 31 Jul 07 nicklas 597     }
3615 31 Jul 07 nicklas 598     return configureJob;
3877 25 Oct 07 nicklas 599   }
2036 21 Feb 06 enell 600
4068 17 Dec 07 nicklas 601   private PluginParameter<ReporterList> getExcludeReportersParameters(DbControl dc)
4068 17 Dec 07 nicklas 602   {
4068 17 Dec 07 nicklas 603     // Load reporter types and initialise the reporterTypeParameter
4068 17 Dec 07 nicklas 604     ItemQuery<ReporterList> query = ReporterList.getQuery();
4068 17 Dec 07 nicklas 605     query.order(Orders.asc(Hql.property("name")));
4068 17 Dec 07 nicklas 606     query.include(Include.MINE, Include.SHARED, Include.IN_PROJECT);
4068 17 Dec 07 nicklas 607     List<ReporterList> reporterLists = new ArrayList<ReporterList>(query.list(dc));
4068 17 Dec 07 nicklas 608     if (reporterLists.size() == 0) return null;
4068 17 Dec 07 nicklas 609     
4068 17 Dec 07 nicklas 610     ItemParameterType<ReporterList> reporterListType = 
4068 17 Dec 07 nicklas 611       new ItemParameterType<ReporterList>(ReporterList.class, null, false, 1, reporterLists);
4068 17 Dec 07 nicklas 612     PluginParameter<ReporterList> reporterListParameter = new PluginParameter<ReporterList>(
4068 17 Dec 07 nicklas 613       "excludeReporters",
4068 17 Dec 07 nicklas 614       "Exclude reporters",
4068 17 Dec 07 nicklas 615       "Spots with reporters in the selected list will not be used in the lowess calculations. " +
4068 17 Dec 07 nicklas 616       "If no list is selected all spots are used.",
4068 17 Dec 07 nicklas 617       reporterListType);
4068 17 Dec 07 nicklas 618     return reporterListParameter;
4068 17 Dec 07 nicklas 619   }
4068 17 Dec 07 nicklas 620   
4811 13 Mar 09 nicklas 621   private static List<Double> lowess(List<SpotData> data, double f, int iter, double delta)
2036 21 Feb 06 enell 622   {
3877 25 Oct 07 nicklas 623     int dataSize = data.size();
3877 25 Oct 07 nicklas 624     Double[] smoothCurve = new Double[dataSize];
3877 25 Oct 07 nicklas 625     int windowSize = Math.min(dataSize, (int) (dataSize * f + 0.5));
2036 21 Feb 06 enell 626
2036 21 Feb 06 enell 627     List<Double> wFit = new ArrayList<Double>();
3877 25 Oct 07 nicklas 628     wFit.addAll(Collections.nCopies(dataSize, 1D));
3877 25 Oct 07 nicklas 629     for (int iteration = 0; iteration < iter; ++iteration)
2036 21 Feb 06 enell 630     {
4811 13 Mar 09 nicklas 631       if (Thread.interrupted()) return null;
2036 21 Feb 06 enell 632       int windowStart = 0;
2036 21 Feb 06 enell 633       int i = 0;
2036 21 Feb 06 enell 634       int prevI = -1;
3877 25 Oct 07 nicklas 635       while (prevI < dataSize - 1)
2036 21 Feb 06 enell 636       {
3877 25 Oct 07 nicklas 637         double Ai = data.get(i).A;
3877 25 Oct 07 nicklas 638         // center window around the i:th value in A
2036 21 Feb 06 enell 639         // while distance from windowStart to i is greater then distance from i to windowEnd: move window
3877 25 Oct 07 nicklas 640         while ((windowStart + windowSize < dataSize) && 
4371 03 Jul 08 nicklas 641             (Ai - data.get(windowStart).A > data.get(windowStart + windowSize-1).A - Ai))
2036 21 Feb 06 enell 642         {
2036 21 Feb 06 enell 643           windowStart++;
2036 21 Feb 06 enell 644         }
2061 23 Feb 06 enell 645         
3877 25 Oct 07 nicklas 646         List<SpotData> window = data.subList(windowStart, windowStart + windowSize);
2061 23 Feb 06 enell 647         List<Double> wFitWindow = wFit.subList(windowStart, windowStart + windowSize);
2061 23 Feb 06 enell 648         
3877 25 Oct 07 nicklas 649         List<Double> w = calculateWeights(window, Ai, wFitWindow);
3877 25 Oct 07 nicklas 650         double[] km = weightedLeastSquaresRegression(window, w);
3877 25 Oct 07 nicklas 651         smoothCurve[i] = km[0] * Ai + km[1];
2036 21 Feb 06 enell 652
2036 21 Feb 06 enell 653         // Interpolate skipped points due to delta
2036 21 Feb 06 enell 654         if (prevI + 1 < i)
2036 21 Feb 06 enell 655         {
3877 25 Oct 07 nicklas 656           double d = Ai - data.get(prevI).A;
2036 21 Feb 06 enell 657           if (d == 0)
2036 21 Feb 06 enell 658           {
2036 21 Feb 06 enell 659             for (int j = prevI + 1; j < i; j++)
2036 21 Feb 06 enell 660             {
2036 21 Feb 06 enell 661               smoothCurve[j] = smoothCurve[i];
2036 21 Feb 06 enell 662             }
2036 21 Feb 06 enell 663           }
2036 21 Feb 06 enell 664           else
2036 21 Feb 06 enell 665           {
2036 21 Feb 06 enell 666             for (int j = prevI + 1; j < i; j++)
2036 21 Feb 06 enell 667             {
3877 25 Oct 07 nicklas 668               double t = (data.get(j).A - data.get(prevI).A) / d;
3877 25 Oct 07 nicklas 669               smoothCurve[j] = t * smoothCurve[i] + (1D - t) * smoothCurve[prevI];
2036 21 Feb 06 enell 670             }
2036 21 Feb 06 enell 671           }
2036 21 Feb 06 enell 672         }
2036 21 Feb 06 enell 673
2036 21 Feb 06 enell 674         // increase i, next x value must be at least delta greater then current
2036 21 Feb 06 enell 675         prevI = i;
3877 25 Oct 07 nicklas 676         double cut = Ai + delta;
3877 25 Oct 07 nicklas 677         while (i < dataSize-1 && data.get(i).A <= cut)
2036 21 Feb 06 enell 678         {
2036 21 Feb 06 enell 679           i++;
2036 21 Feb 06 enell 680         }
2036 21 Feb 06 enell 681       }
2036 21 Feb 06 enell 682       
3877 25 Oct 07 nicklas 683       double invYWRange = 1D/(ywRangeFactor  * medianCorrection(data, smoothCurve));
2036 21 Feb 06 enell 684       for (int j = 0; j < wFit.size(); j++)
2036 21 Feb 06 enell 685       {
3877 25 Oct 07 nicklas 686         double w = Math.abs((smoothCurve[j] - data.get(j).M) * invYWRange);
4371 03 Jul 08 nicklas 687         wFit.set(j, w < 1 ? Math.pow(1D - w*w, 2) : 0);
2036 21 Feb 06 enell 688       }
2036 21 Feb 06 enell 689     }
2036 21 Feb 06 enell 690     return Arrays.asList(smoothCurve);
2036 21 Feb 06 enell 691   }
2036 21 Feb 06 enell 692
3877 25 Oct 07 nicklas 693   private static double medianCorrection(List<SpotData> data, Double[] smoothCurve)
2036 21 Feb 06 enell 694   {
3877 25 Oct 07 nicklas 695     List<Double> temp = new ArrayList<Double>(data.size());
3877 25 Oct 07 nicklas 696     for (int i = 0; i < data.size(); i++)
2036 21 Feb 06 enell 697     {
3877 25 Oct 07 nicklas 698       temp.add(Math.abs(data.get(i).M - smoothCurve[i]));
2036 21 Feb 06 enell 699     }
2036 21 Feb 06 enell 700     Collections.sort(temp);
2036 21 Feb 06 enell 701     if (temp.size() % 2 == 0)
2036 21 Feb 06 enell 702     {
2036 21 Feb 06 enell 703       double v1 = temp.get((temp.size() / 2) - 1);
2036 21 Feb 06 enell 704       double v2 = temp.get((temp.size() / 2));
4371 03 Jul 08 nicklas 705       return 0.5 * (v1 + v2);
2036 21 Feb 06 enell 706     }
4371 03 Jul 08 nicklas 707     return temp.get(temp.size() / 2);
2036 21 Feb 06 enell 708   }
2036 21 Feb 06 enell 709
2036 21 Feb 06 enell 710   
3877 25 Oct 07 nicklas 711   private static double[] weightedLeastSquaresRegression(List<SpotData> data, List<Double> w)
2036 21 Feb 06 enell 712   {
2036 21 Feb 06 enell 713     double k;
2036 21 Feb 06 enell 714     double m;
3877 25 Oct 07 nicklas 715     double sumA = 0;
3877 25 Oct 07 nicklas 716     double sumM = 0;
3877 25 Oct 07 nicklas 717     double sumAA = 0;
3877 25 Oct 07 nicklas 718     double sumAM = 0;
2036 21 Feb 06 enell 719     double sumW = 0;
3877 25 Oct 07 nicklas 720     for (int j = 0; j < data.size(); j++)
2036 21 Feb 06 enell 721     {
4068 17 Dec 07 nicklas 722       SpotData spot = data.get(j);
4068 17 Dec 07 nicklas 723       if (!spot.exclude)
4068 17 Dec 07 nicklas 724       {
4068 17 Dec 07 nicklas 725         double localW = w.get(j);
4068 17 Dec 07 nicklas 726         double localA = spot.A;
4068 17 Dec 07 nicklas 727         double localM = spot.M;
4068 17 Dec 07 nicklas 728         sumA += localA * localW;
4068 17 Dec 07 nicklas 729         sumM += localM * localW;
4068 17 Dec 07 nicklas 730         sumAA += localA * localA * localW;
4068 17 Dec 07 nicklas 731         sumAM += localA * localM * localW;
4068 17 Dec 07 nicklas 732         sumW += localW;
4068 17 Dec 07 nicklas 733       }
2036 21 Feb 06 enell 734     }
2036 21 Feb 06 enell 735     if (sumW <= 0)
2036 21 Feb 06 enell 736     {
4062 14 Dec 07 nicklas 737       throw new BaseException("Sum of weigths in line_fit is not positive. " +
4062 14 Dec 07 nicklas 738           "This may happen if there are too few data points in a block group. Please " +
4062 14 Dec 07 nicklas 739           "increase the 'block group' parameter, or set it to 0 to include all blocks " +
4062 14 Dec 07 nicklas 740           "in a single group");
2036 21 Feb 06 enell 741     }
3877 25 Oct 07 nicklas 742     double denom = sumW * sumAA - sumA * sumA;
2036 21 Feb 06 enell 743     if (denom != 0.)
2036 21 Feb 06 enell 744     {
3877 25 Oct 07 nicklas 745       k = (sumW * sumAM - sumA * sumM) / denom;
3877 25 Oct 07 nicklas 746       m = (sumM - k * sumA) / sumW;
2036 21 Feb 06 enell 747     }
2036 21 Feb 06 enell 748     else
2036 21 Feb 06 enell 749     {
2036 21 Feb 06 enell 750       k = 0.;
3877 25 Oct 07 nicklas 751       m = sumM / sumW;
2036 21 Feb 06 enell 752     }
2036 21 Feb 06 enell 753
2036 21 Feb 06 enell 754     return new double[] {k, m};
2036 21 Feb 06 enell 755   }
2036 21 Feb 06 enell 756
3877 25 Oct 07 nicklas 757   private static List<Double> calculateWeights(List<SpotData> data, double A1, List<Double> wFit)
2036 21 Feb 06 enell 758   {
3877 25 Oct 07 nicklas 759     List<Double> w = new ArrayList<Double>(data.size());
2036 21 Feb 06 enell 760     double invRadius = 0;
3877 25 Oct 07 nicklas 761     for (SpotData spot : data)
2036 21 Feb 06 enell 762     {
3877 25 Oct 07 nicklas 763       double abs = Math.abs(A1 - spot.A);
2561 22 Aug 06 nicklas 764       if (abs > invRadius) invRadius = abs;
2036 21 Feb 06 enell 765     }
2036 21 Feb 06 enell 766     invRadius = 1 / invRadius;
3877 25 Oct 07 nicklas 767     for (int i = 0; i < data.size(); i++)
2036 21 Feb 06 enell 768     {
4068 17 Dec 07 nicklas 769       if (!data.get(i).exclude)
4068 17 Dec 07 nicklas 770       {
4068 17 Dec 07 nicklas 771         double A2 = data.get(i).A;
4068 17 Dec 07 nicklas 772         double distance = Math.abs(A1 - A2) * invRadius;
4068 17 Dec 07 nicklas 773         w.add((distance < 1 ? Math.pow(1D - Math.pow(Math.abs(distance), 3), 3) : 0) * wFit.get(i));
4068 17 Dec 07 nicklas 774       }
4068 17 Dec 07 nicklas 775       else
4068 17 Dec 07 nicklas 776       {
4068 17 Dec 07 nicklas 777         w.add(0.0);
4068 17 Dec 07 nicklas 778       }
2036 21 Feb 06 enell 779     }
2036 21 Feb 06 enell 780     return w;
2036 21 Feb 06 enell 781   }
3877 25 Oct 07 nicklas 782
3877 25 Oct 07 nicklas 783   private static class SpotData
3887 29 Oct 07 nicklas 784     implements Comparable<SpotData>
3877 25 Oct 07 nicklas 785   {
4811 13 Mar 09 nicklas 786     final short column;
3877 25 Oct 07 nicklas 787     final int position;
3877 25 Oct 07 nicklas 788     final float ch1;
3877 25 Oct 07 nicklas 789     final float ch2;
3877 25 Oct 07 nicklas 790     final double M;
3877 25 Oct 07 nicklas 791     final double A;
3877 25 Oct 07 nicklas 792     final int block;
4068 17 Dec 07 nicklas 793     final boolean exclude;
4811 13 Mar 09 nicklas 794     float newCh1;
4811 13 Mar 09 nicklas 795     float newCh2;
3879 26 Oct 07 nicklas 796     static final double INV_LN2 = 1 / Math.log(2);
3877 25 Oct 07 nicklas 797
6372 06 Dec 13 nicklas 798     public SpotData(short column, int position, float ch1, float ch2, int block, boolean exclude, MACalculator mac)
3877 25 Oct 07 nicklas 799     {
4811 13 Mar 09 nicklas 800       this.column = column;
3877 25 Oct 07 nicklas 801       this.position = position;
3877 25 Oct 07 nicklas 802       this.ch1 = ch1;
3877 25 Oct 07 nicklas 803       this.ch2 = ch2;
3877 25 Oct 07 nicklas 804       this.block = block;
4068 17 Dec 07 nicklas 805       this.exclude = exclude;
6372 06 Dec 13 nicklas 806       double[] MA = mac.MA(ch1, ch2);
6372 06 Dec 13 nicklas 807       this.M = MA[0];
6372 06 Dec 13 nicklas 808       this.A = MA[1];
3877 25 Oct 07 nicklas 809     }
3887 29 Oct 07 nicklas 810
3887 29 Oct 07 nicklas 811     /**
3887 29 Oct 07 nicklas 812       Sort by A in ascending order.
3887 29 Oct 07 nicklas 813     */
6127 14 Sep 12 nicklas 814     @Override
3887 29 Oct 07 nicklas 815     public int compareTo(SpotData o)
3887 29 Oct 07 nicklas 816     {
3887 29 Oct 07 nicklas 817       return this.A < o.A ? -1 : 1;
3887 29 Oct 07 nicklas 818     }
3887 29 Oct 07 nicklas 819     
3877 25 Oct 07 nicklas 820   }
4813 16 Mar 09 nicklas 821   
5384 13 Aug 10 nicklas 822   static class CallableWorker
4813 16 Mar 09 nicklas 823     implements Callable<List<SpotData>>
4811 13 Mar 09 nicklas 824   {
4813 16 Mar 09 nicklas 825   
4813 16 Mar 09 nicklas 826     private final List<SpotData> in;
6372 06 Dec 13 nicklas 827     private final MACalculator mac;
4811 13 Mar 09 nicklas 828     private final float fitFraction;
4811 13 Mar 09 nicklas 829     private final float delta;
4811 13 Mar 09 nicklas 830     private final int iterations;
4811 13 Mar 09 nicklas 831     private int blockGroupSize;
4811 13 Mar 09 nicklas 832     
6372 06 Dec 13 nicklas 833     CallableWorker(List<SpotData> in, MACalculator mac,
4811 13 Mar 09 nicklas 834       float fitFraction, float delta, int iterations, int blockGroupSize)
4811 13 Mar 09 nicklas 835     {
4811 13 Mar 09 nicklas 836       this.in = in;
6372 06 Dec 13 nicklas 837       this.mac = mac;
4811 13 Mar 09 nicklas 838       this.fitFraction = fitFraction;
4811 13 Mar 09 nicklas 839       this.delta = delta;
4811 13 Mar 09 nicklas 840       this.iterations = iterations;
4811 13 Mar 09 nicklas 841       this.blockGroupSize = blockGroupSize;
4811 13 Mar 09 nicklas 842     }
4813 16 Mar 09 nicklas 843   
4811 13 Mar 09 nicklas 844     @Override
4813 16 Mar 09 nicklas 845     public List<SpotData> call()
4811 13 Mar 09 nicklas 846     {
6372 06 Dec 13 nicklas 847       normalize(in, mac);
4813 16 Mar 09 nicklas 848       return in;
4811 13 Mar 09 nicklas 849     }
4811 13 Mar 09 nicklas 850     
4811 13 Mar 09 nicklas 851     /**
4811 13 Mar 09 nicklas 852       Normalize the spot data.
4811 13 Mar 09 nicklas 853       @return TRUE if all was ok, FALSE if worker has been interrupted and should
4811 13 Mar 09 nicklas 854         stop executing
4811 13 Mar 09 nicklas 855     */
6372 06 Dec 13 nicklas 856     private boolean normalize(List<SpotData> data, MACalculator mac)
4811 13 Mar 09 nicklas 857     {
4811 13 Mar 09 nicklas 858       int dataSize = data.size();
4811 13 Mar 09 nicklas 859       
4811 13 Mar 09 nicklas 860       // Get range of block numbers - NOTE! query must return spots sorted in block order
4811 13 Mar 09 nicklas 861       int minBlock = data.get(0).block;
4811 13 Mar 09 nicklas 862       int maxBlock = data.get(data.size()-1).block;
4811 13 Mar 09 nicklas 863       // If 0, we should use all blocks as a single data set
4811 13 Mar 09 nicklas 864       if (blockGroupSize <= 0) blockGroupSize = maxBlock - minBlock + 1;
4811 13 Mar 09 nicklas 865       
4811 13 Mar 09 nicklas 866       int fromIndex = 0;
4811 13 Mar 09 nicklas 867       int toIndex = 0;
4811 13 Mar 09 nicklas 868       int fromBlock = minBlock;
4811 13 Mar 09 nicklas 869       // Normalize each block range independently: fromBlock + blockGroupSize --> toBlock
4811 13 Mar 09 nicklas 870       while (fromBlock <= maxBlock)
4811 13 Mar 09 nicklas 871       {
4811 13 Mar 09 nicklas 872         if (Thread.interrupted()) return false;
4811 13 Mar 09 nicklas 873         
4811 13 Mar 09 nicklas 874         // Find start and end index for current block range
4811 13 Mar 09 nicklas 875         int toBlock = fromBlock + blockGroupSize - 1;
4811 13 Mar 09 nicklas 876         if (toBlock > maxBlock) toBlock = maxBlock;
4811 13 Mar 09 nicklas 877         fromIndex = toIndex;
5689 11 Aug 11 nicklas 878         // Data is sorted by block. Find index of last spot with: block <= toBlock
4811 13 Mar 09 nicklas 879         // spot given by toIndex should not be included
4811 13 Mar 09 nicklas 880         while (toIndex < dataSize && data.get(toIndex).block <= toBlock)
4811 13 Mar 09 nicklas 881         {
4811 13 Mar 09 nicklas 882           ++toIndex;
4811 13 Mar 09 nicklas 883         }
4811 13 Mar 09 nicklas 884         
4811 13 Mar 09 nicklas 885         if (toIndex > fromIndex)
4811 13 Mar 09 nicklas 886         {
4811 13 Mar 09 nicklas 887           List <SpotData> toNormalize = data.subList(fromIndex, toIndex);
4811 13 Mar 09 nicklas 888           Collections.sort(toNormalize);
4811 13 Mar 09 nicklas 889           List<Double> smoothCurve = lowess(toNormalize, fitFraction, iterations, delta);
4811 13 Mar 09 nicklas 890           if (smoothCurve == null) return false;
4811 13 Mar 09 nicklas 891           for (int j = 0; j < smoothCurve.size(); ++j)
4811 13 Mar 09 nicklas 892           {
4811 13 Mar 09 nicklas 893             if (Thread.interrupted()) return false;
4811 13 Mar 09 nicklas 894             SpotData spot = toNormalize.get(j);
6372 06 Dec 13 nicklas 895             double factor = smoothCurve.get(j) * 0.5;
6372 06 Dec 13 nicklas 896             double[] newCh = mac.correct(spot.ch1, spot.ch2, factor);
6372 06 Dec 13 nicklas 897             spot.newCh1 = (float)(newCh[0]);
6372 06 Dec 13 nicklas 898             spot.newCh2 = (float)(newCh[1]);
4811 13 Mar 09 nicklas 899           }
4811 13 Mar 09 nicklas 900         }
4811 13 Mar 09 nicklas 901         fromBlock = toBlock + 1;
4811 13 Mar 09 nicklas 902       }
4811 13 Mar 09 nicklas 903       return true;
4811 13 Mar 09 nicklas 904     }
4811 13 Mar 09 nicklas 905   }
4813 16 Mar 09 nicklas 906     
2036 21 Feb 06 enell 907 }

Generated on Sat Oct  7 02:40:44 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)