HMSBEAGLE  1.0.0
BeagleCPUImpl.h
1 /*
2  * BeagleCPUImpl.h
3  * BEAGLE
4  *
5  * Copyright 2009 Phylogenetic Likelihood Working Group
6  *
7  * This file is part of BEAGLE.
8  *
9  * BEAGLE is free software: you can redistribute it and/or modify
10  * it under the terms of the GNU Lesser General Public License as
11  * published by the Free Software Foundation, either version 3 of
12  * the License, or (at your option) any later version.
13  *
14  * BEAGLE is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with BEAGLE. If not, see
21  * <http://www.gnu.org/licenses/>.
22  *
23  * @author Andrew Rambaut
24  * @author Marc Suchard
25  * @author Daniel Ayres
26  */
27 
28 #ifndef __BeagleCPUImpl__
29 #define __BeagleCPUImpl__
30 
31 #ifdef HAVE_CONFIG_H
32 #include "libhmsbeagle/config.h"
33 #endif
34 
35 #include "libhmsbeagle/BeagleImpl.h"
36 #include "libhmsbeagle/CPU/Precision.h"
37 #include "libhmsbeagle/CPU/EigenDecomposition.h"
38 
39 #include <vector>
40 
41 #define BEAGLE_CPU_GENERIC REALTYPE, T_PAD, P_PAD
42 #define BEAGLE_CPU_TEMPLATE template <typename REALTYPE, int T_PAD, int P_PAD>
43 
44 #define BEAGLE_CPU_FACTORY_GENERIC REALTYPE
45 #define BEAGLE_CPU_FACTORY_TEMPLATE template <typename REALTYPE>
46 
47 
48 #define T_PAD_DEFAULT 1 // Pad transition matrix rows with an extra 1.0 for ambiguous characters
49 #define P_PAD_DEFAULT 0 // No partials padding necessary for non-SSE implementations
50 
51 
52 namespace beagle {
53 namespace cpu {
54 
55 BEAGLE_CPU_TEMPLATE
56 class BeagleCPUImpl : public BeagleImpl {
57 
58 protected:
59  int kBufferCount;
60 
61  int kTipCount;
62 
69  int kPartialsPaddedStateCount;
70  int kEigenDecompCount;
72  int kScaleBufferCount;
73 
74  int kPartialsSize;
76 
78 
79  long kFlags;
80 
81  REALTYPE realtypeMin;
82  int scalingExponentThreshhold;
83 
85 
86  double* gCategoryRates; // Kept in double-precision until multiplication by edgelength
87  double* gPatternWeights;
88 
89  REALTYPE** gCategoryWeights;
90  REALTYPE** gStateFrequencies;
91 
92  //@ the size of these pointers are known at alloc-time, so the partials and
93  // tipStates field should be switched to vectors of vectors (to make
94  // memory management less error prone
95  REALTYPE** gPartials;
96  int** gTipStates;
97  REALTYPE** gScaleBuffers;
98 
99  signed short** gAutoScaleBuffers;
100 
101  int* gActiveScalingFactors;
102 
103  // There will be kMatrixCount transitionMatrices.
104  // Each kStateCount x (kStateCount+1) matrix that is flattened
105  // into a single array
106  REALTYPE** gTransitionMatrices;
107 
108  REALTYPE* integrationTmp;
109  REALTYPE* firstDerivTmp;
110  REALTYPE* secondDerivTmp;
111 
112  REALTYPE* outLogLikelihoodsTmp;
113  REALTYPE* outFirstDerivativesTmp;
114  REALTYPE* outSecondDerivativesTmp;
115 
116  REALTYPE* ones;
117  REALTYPE* zeros;
118 
119 public:
120  virtual ~BeagleCPUImpl();
121 
122  // creation of instance
123  int createInstance(int tipCount,
124  int partialsBufferCount,
125  int compactBufferCount,
126  int stateCount,
127  int patternCount,
128  int eigenDecompositionCount,
129  int matrixCount,
130  int categoryCount,
131  int scaleBufferCount,
132  int resourceNumber,
133  long preferenceFlags,
134  long requirementFlags);
135 
136  // initialization of instance, returnInfo can be null
137  int getInstanceDetails(BeagleInstanceDetails* returnInfo);
138 
139  // set the states for a given tip
140  //
141  // tipIndex the index of the tip
142  // inStates the array of states: 0 to stateCount - 1, missing = stateCount
143  int setTipStates(int tipIndex,
144  const int* inStates);
145 
146  // set the partials for a given tip
147  //
148  // tipIndex the index of the tip
149  // inPartials the array of partials, stateCount x patternCount
150  int setTipPartials(int tipIndex,
151  const double* inPartials);
152 
153 
154  int setPartials(int bufferIndex,
155  const double* inPartials);
156 
157  int getPartials(int bufferIndex,
158  int scaleBuffer,
159  double* outPartials);
160 
161  // sets the Eigen decomposition for a given matrix
162  //
163  // matrixIndex the matrix index to update
164  // eigenVectors an array containing the Eigen Vectors
165  // inverseEigenVectors an array containing the inverse Eigen Vectors
166  // eigenValues an array containing the Eigen Values
167  int setEigenDecomposition(int eigenIndex,
168  const double* inEigenVectors,
169  const double* inInverseEigenVectors,
170  const double* inEigenValues);
171 
172  int setStateFrequencies(int stateFrequenciesIndex,
173  const double* inStateFrequencies);
174 
175  int setCategoryWeights(int categoryWeightsIndex,
176  const double* inCategoryWeights);
177 
178  int setPatternWeights(const double* inPatternWeights);
179 
180  // set the vector of category rates
181  //
182  // categoryRates an array containing categoryCount rate scalers
183  int setCategoryRates(const double* inCategoryRates);
184 
185  int setTransitionMatrix(int matrixIndex,
186  const double* inMatrix,
187  double paddedValue);
188 
189  int setTransitionMatrices(const int* matrixIndices,
190  const double* inMatrices,
191  const double* paddedValues,
192  int count);
193 
194  int getTransitionMatrix(int matrixIndex,
195  double* outMatrix);
196 
197  // calculate a transition probability matrices for a given list of node. This will
198  // calculate for all categories (and all matrices if more than one is being used).
199  //
200  // nodeIndices an array of node indices that require transition probability matrices
201  // edgeLengths an array of expected lengths in substitutions per site
202  // count the number of elements in the above arrays
203  int updateTransitionMatrices(int eigenIndex,
204  const int* probabilityIndices,
205  const int* firstDerivativeIndices,
206  const int* secondDerivativeIndices,
207  const double* edgeLengths,
208  int count);
209 
210  // calculate or queue for calculation partials using an array of operations
211  //
212  // operations an array of triplets of indices: the two source partials and the destination
213  // dependencies an array of indices specify which operations are dependent on which (optional)
214  // count the number of operations
215  // rescale indicate if partials should be rescaled during peeling
216  int updatePartials(const int* operations,
217  int operationCount,
218  int cumulativeScalingIndex);
219 
220  // Block until all calculations that write to the specified partials have completed.
221  //
222  // This function is optional and only has to be called by clients that "recycle" partials.
223  //
224  // If used, this function must be called after an updatePartials call and must refer to
225  // indices of "destinationPartials" that were used in a previous updatePartials
226  // call. The library will block until those partials have been calculated.
227  //
228  // destinationPartials - List of the indices of destinationPartials that must be calculated
229  // before the function returns
230  // destinationPartialsCount - Number of destinationPartials (input)
231  //
232  // return error code
233  int waitForPartials(const int* destinationPartials,
234  int destinationPartialsCount);
235 
236 
237  int accumulateScaleFactors(const int* scalingIndices,
238  int count,
239  int cumulativeScalingIndex);
240 
241  int removeScaleFactors(const int* scalingIndices,
242  int count,
243  int cumulativeScalingIndex);
244 
245  int resetScaleFactors(int cumulativeScalingIndex);
246 
247  int copyScaleFactors(int destScalingIndex,
248  int srcScalingIndex);
249 
250  // calculate the site log likelihoods at a particular node
251  //
252  // rootNodeIndex the index of the root
253  // outLogLikelihoods an array into which the site log likelihoods will be put
254  int calculateRootLogLikelihoods(const int* bufferIndices,
255  const int* categoryWeightsIndices,
256  const int* stateFrequenciesIndices,
257  const int* cumulativeScaleIndices,
258  int count,
259  double* outSumLogLikelihood);
260 
261  // possible nulls: firstDerivativeIndices, secondDerivativeIndices,
262  // outFirstDerivatives, outSecondDerivatives
263  int calculateEdgeLogLikelihoods(const int* parentBufferIndices,
264  const int* childBufferIndices,
265  const int* probabilityIndices,
266  const int* firstDerivativeIndices,
267  const int* secondDerivativeIndices,
268  const int* categoryWeightsIndices,
269  const int* stateFrequenciesIndices,
270  const int* cumulativeScaleIndices,
271  int count,
272  double* outSumLogLikelihood,
273  double* outSumFirstDerivative,
274  double* outSumSecondDerivative);
275 
276  int getSiteLogLikelihoods(double* outLogLikelihoods);
277 
278  int getSiteDerivatives(double* outFirstDerivatives,
279  double* outSecondDerivatives);
280 
281  int block(void);
282 
283  virtual const char* getName();
284 
285  virtual const long getFlags();
286 
287 protected:
288  virtual void calcStatesStates(REALTYPE* destP,
289  const int* states1,
290  const REALTYPE* matrices1,
291  const int* states2,
292  const REALTYPE* matrices2);
293 
294 
295  virtual void calcStatesPartials(REALTYPE* destP,
296  const int* states1,
297  const REALTYPE* matrices1,
298  const REALTYPE* partials2,
299  const REALTYPE* matrices2);
300 
301  virtual void calcPartialsPartials(REALTYPE* destP,
302  const REALTYPE* partials1,
303  const REALTYPE* matrices1,
304  const REALTYPE* partials2,
305  const REALTYPE* matrices2);
306 
307  virtual int calcRootLogLikelihoods(const int bufferIndex,
308  const int categoryWeightsIndex,
309  const int stateFrequenciesIndex,
310  const int scaleBufferIndex,
311  double* outSumLogLikelihood);
312 
313  virtual int calcRootLogLikelihoodsMulti(const int* bufferIndices,
314  const int* categoryWeightsIndices,
315  const int* stateFrequenciesIndices,
316  const int* scaleBufferIndices,
317  int count,
318  double* outSumLogLikelihood);
319 
320  virtual int calcEdgeLogLikelihoods(const int parentBufferIndex,
321  const int childBufferIndex,
322  const int probabilityIndex,
323  const int categoryWeightsIndex,
324  const int stateFrequenciesIndex,
325  const int scalingFactorsIndex,
326  double* outSumLogLikelihood);
327 
328  virtual int calcEdgeLogLikelihoodsMulti(const int* parentBufferIndices,
329  const int* childBufferIndices,
330  const int* probabilityIndices,
331  const int* categoryWeightsIndices,
332  const int* stateFrequenciesIndices,
333  const int* scalingFactorsIndices,
334  int count,
335  double* outSumLogLikelihood);
336 
337  virtual int calcEdgeLogLikelihoodsFirstDeriv(const int parentBufferIndex,
338  const int childBufferIndex,
339  const int probabilityIndex,
340  const int firstDerivativeIndex,
341  const int categoryWeightsIndex,
342  const int stateFrequenciesIndex,
343  const int scalingFactorsIndex,
344  double* outSumLogLikelihood,
345  double* outSumFirstDerivative);
346 
347  virtual int calcEdgeLogLikelihoodsSecondDeriv(const int parentBufferIndex,
348  const int childBufferIndex,
349  const int probabilityIndex,
350  const int firstDerivativeIndex,
351  const int secondDerivativeIndex,
352  const int categoryWeightsIndex,
353  const int stateFrequenciesIndex,
354  const int scalingFactorsIndex,
355  double* outSumLogLikelihood,
356  double* outSumFirstDerivative,
357  double* outSumSecondDerivative);
358 
359  virtual void calcStatesStatesFixedScaling(REALTYPE *destP,
360  const int *child0States,
361  const REALTYPE *child0TransMat,
362  const int *child1States,
363  const REALTYPE *child1TransMat,
364  const REALTYPE *scaleFactors);
365 
366  virtual void calcStatesPartialsFixedScaling(REALTYPE *destP,
367  const int *child0States,
368  const REALTYPE *child0TransMat,
369  const REALTYPE *child1Partials,
370  const REALTYPE *child1TransMat,
371  const REALTYPE *scaleFactors);
372 
373  virtual void calcPartialsPartialsFixedScaling(REALTYPE *destP,
374  const REALTYPE *child0States,
375  const REALTYPE *child0TransMat,
376  const REALTYPE *child1Partials,
377  const REALTYPE *child1TransMat,
378  const REALTYPE *scaleFactors);
379 
380  virtual void calcPartialsPartialsAutoScaling(REALTYPE* destP,
381  const REALTYPE* partials1,
382  const REALTYPE* matrices1,
383  const REALTYPE* partials2,
384  const REALTYPE* matrices2,
385  int* activateScaling);
386 
387  virtual void rescalePartials(REALTYPE *destP,
388  REALTYPE *scaleFactors,
389  REALTYPE *cumulativeScaleFactors,
390  const int fillWithOnes);
391 
392  virtual void autoRescalePartials(REALTYPE *destP,
393  signed short *scaleFactors);
394 
395  virtual int getPaddedPatternsModulus();
396 
397  void* mallocAligned(size_t size);
398 
399 };
400 
401 BEAGLE_CPU_FACTORY_TEMPLATE
403 public:
404  virtual BeagleImpl* createImpl(int tipCount,
405  int partialsBufferCount,
406  int compactBufferCount,
407  int stateCount,
408  int patternCount,
409  int eigenBufferCount,
410  int matrixBufferCount,
411  int categoryCount,
412  int scaleBufferCount,
413  int resourceNumber,
414  long preferenceFlags,
415  long requirementFlags,
416  int* errorCode);
417 
418  virtual const char* getName();
419  virtual const long getFlags();
420 };
421 
422 //typedef BeagleCPUImplGeneral<double> BeagleCPUImpl;
423 
424 } // namespace cpu
425 } // namespace beagle
426 
427 // now that the interface is defined, include the implementation of template functions
428 #include "libhmsbeagle/CPU/BeagleCPUImpl.hpp"
429 
430 #endif // __BeagleCPUImpl__