← Back to branch summary

~ubuntu-branches/ubuntu/raring/libpal-java/raring

« back to all changes in this revision

Viewing changes to src/pal/substmodel/GTR.java

  • Committer: Package Import Robot
  • Author(s): Andreas Tille
  • Date: 2012-01-27 13:57:54 UTC
  • Revision ID: package-import@ubuntu.com-20120127135754-d63l3581f65fw9pk
Tags: upstream-1.5.1
Import upstream version 1.5.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/pal/substmodel/GTR.java
 
1
// GTR.java
 
2
//
 
3
// (c) 1999-2001 PAL Development Core Team
 
4
//
 
5
// This package may be distributed under the
 
6
// terms of the Lesser GNU General Public License (LGPL)
 
7
 
 
8
 
 
9
package pal.substmodel;
 
10
 
 
11
import pal.misc.*;
 
12
import pal.util.*;
 
13
import java.io.*;
 
14
 
 
15
 
 
16
/**
 
17
 * GTR (general time reversible) model of nucleotide evolution
 
18
 *
 
19
 * <i>Lanave, C., G. Preparata, C. Saccone, and G. Serio. 1984. A new method for calculating evolutionary substitution rates. J Mol Evol 20: 86-93. </i>
 
20
 *
 
21
 * @version $Id: GTR.java,v 1.13 2003/11/30 05:29:22 matt Exp $
 
22
 * <p>
 
23
 * <em>Parameters</em>
 
24
 * <ol>
 
25
 *  <li> A </li> <!-- 0 -->
 
26
 *  <li> B </li> <!-- 1 -->
 
27
 *  <li> C </li> <!-- 2 -->
 
28
 *  <li> D </li> <!-- 3 -->
 
29
 *  <li> E </li> <!-- 4 -->
 
30
 * </ol>
 
31
 * </p>
 
32
 *
 
33
 * @author Korbinian Strimmer
 
34
 */
 
35
public class GTR extends NucleotideModel implements Serializable, XMLConstants
 
36
{
 
37
 
 
38
        //
 
39
        // Private stuff
 
40
        //
 
41
 
 
42
        private boolean showSE;
 
43
        private double a, b, c, d, e;
 
44
        private double aSE, bSE, cSE, dSE, eSE;
 
45
 
 
46
        //
 
47
        // Serialization code
 
48
        //
 
49
        private static final long serialVersionUID= -8557884770092535699L;
 
50
 
 
51
        //serialver -classpath ./classes pal.substmodel.GTR
 
52
        private void writeObject(java.io.ObjectOutputStream out) throws java.io.IOException {
 
53
                out.writeByte(1); //Version number
 
54
                out.writeBoolean(showSE);
 
55
                out.writeDouble(a);
 
56
                out.writeDouble(b);
 
57
                out.writeDouble(c);
 
58
                out.writeDouble(d);
 
59
                out.writeDouble(e);
 
60
                out.writeDouble(aSE);
 
61
                out.writeDouble(bSE);
 
62
                out.writeDouble(cSE);
 
63
                out.writeDouble(dSE);
 
64
                out.writeDouble(eSE);
 
65
        }
 
66
 
 
67
        private void readObject(java.io.ObjectInputStream in) throws java.io.IOException, ClassNotFoundException{
 
68
                byte version = in.readByte();
 
69
                switch(version) {
 
70
                        default : {
 
71
                                showSE = in.readBoolean();
 
72
                                a = in.readDouble();
 
73
                                b = in.readDouble();
 
74
                                c = in.readDouble();
 
75
                                d = in.readDouble();
 
76
                                e = in.readDouble();
 
77
                                aSE = in.readDouble();
 
78
                                bSE = in.readDouble();
 
79
                                cSE = in.readDouble();
 
80
                                dSE = in.readDouble();
 
81
                                eSE = in.readDouble();
 
82
                                break;
 
83
                        }
 
84
                }
 
85
        }
 
86
 
 
87
 
 
88
        /**
 
89
         * constructor 1
 
90
         *
 
91
         * @param a entry in rate matrix
 
92
         * @param b entry in rate matrix
 
93
         * @param c entry in rate matrix
 
94
         * @param d entry in rate matrix
 
95
         * @param e entry in rate matrix
 
96
         * @param freq nucleotide frequencies
 
97
         */
 
98
        public GTR(double a, double b, double c, double d, double e, double[] freq)     {
 
99
                super(freq);
 
100
                this.a = a; this.b = b; this.c = c; this.d = d; this.e = e;
 
101
                setParameters(new double[] {a,b,c,d,e});
 
102
                showSE = false;
 
103
        }
 
104
 
 
105
        /**
 
106
         * constructor 2
 
107
         *
 
108
         * @param params parameter list
 
109
         * @param freq nucleotide frequencies
 
110
         */
 
111
        public GTR(double[] params, double[] freq)      {
 
112
                this(params[0], params[1], params[2],
 
113
                        params[3], params[4], freq);
 
114
        }
 
115
 
 
116
        public Object clone() {
 
117
                return new GTR(this);
 
118
        }
 
119
 
 
120
        private GTR(GTR gtr) {
 
121
                this(gtr.a, gtr.b, gtr.c, gtr.d, gtr.e, gtr.getEquilibriumFrequencies());
 
122
        }
 
123
 
 
124
        // Get numerical code describing the model type
 
125
        public int getModelID() {
 
126
                return 0;
 
127
        }
 
128
 
 
129
        // interface Report
 
130
 
 
131
        public void report(PrintWriter out)
 
132
        {
 
133
                out.println("Model of substitution: GTR (Lanave et al. 1984)");
 
134
 
 
135
                out.print("Parameter a: ");
 
136
                format.displayDecimal(out, a, 2);
 
137
                if (showSE)
 
138
                {
 
139
                        out.print("  (S.E. ");
 
140
                        format.displayDecimal(out, aSE, 2);
 
141
                        out.print(")");
 
142
                }
 
143
                out.println();
 
144
 
 
145
                out.print("Parameter b: ");
 
146
                format.displayDecimal(out, b, 2);
 
147
                if (showSE)
 
148
                {
 
149
                        out.print("  (S.E. ");
 
150
                        format.displayDecimal(out, bSE, 2);
 
151
                        out.print(")");
 
152
                }
 
153
                out.println();
 
154
 
 
155
                out.print("Parameter c: ");
 
156
                format.displayDecimal(out, c, 2);
 
157
                if (showSE)
 
158
                {
 
159
                        out.print("  (S.E. ");
 
160
                        format.displayDecimal(out, cSE, 2);
 
161
                        out.print(")");
 
162
                }
 
163
                out.println();
 
164
 
 
165
                out.print("Parameter d: ");
 
166
                format.displayDecimal(out, d, 2);
 
167
                if (showSE)
 
168
                {
 
169
                        out.print("  (S.E. ");
 
170
                        format.displayDecimal(out, dSE, 2);
 
171
                        out.print(")");
 
172
                }
 
173
                out.println();
 
174
 
 
175
                out.print("Parameter e: ");
 
176
                format.displayDecimal(out, e, 2);
 
177
                if (showSE)
 
178
                {
 
179
                        out.print("  (S.E. ");
 
180
                        format.displayDecimal(out, eSE, 2);
 
181
                        out.print(")");
 
182
                }
 
183
                out.println();
 
184
 
 
185
                out.println("                                   A  C  G  T");
 
186
                out.println("Corresponding rate matrix      ----------------");
 
187
                out.println("(shown without frequencies):     A    a  b  c");
 
188
                out.println("                                 C       d  e");
 
189
                out.println("                                 G          1");
 
190
 
 
191
                out.println();
 
192
                printFrequencies(out);
 
193
                printRatios(out);
 
194
        }
 
195
 
 
196
        // interface Parameterized
 
197
 
 
198
        public int getNumParameters()   {
 
199
                return 5;
 
200
        }
 
201
 
 
202
        public void setParameterSE(double paramSE, int n)       {
 
203
                switch(n)       {
 
204
                        case 0: aSE = paramSE; break;
 
205
                        case 1: bSE = paramSE; break;
 
206
                        case 2: cSE = paramSE; break;
 
207
                        case 3: dSE = paramSE; break;
 
208
                        case 4: eSE = paramSE; break;
 
209
 
 
210
                        default: throw new IllegalArgumentException();
 
211
                }
 
212
                showSE = true;
 
213
        }
 
214
 
 
215
        public double getLowerLimit(int n)      {
 
216
                return 0.0001;
 
217
        }
 
218
 
 
219
        public double getUpperLimit(int n)      {
 
220
                return 10000.0;
 
221
        }
 
222
 
 
223
        public double getDefaultValue(int n)    {
 
224
                return 1.0;
 
225
        }
 
226
 
 
227
        /**
 
228
         * @return the name of this rate matrix
 
229
         */
 
230
        public String getUniqueName() {
 
231
                return GTR;
 
232
        }
 
233
 
 
234
        public String getParameterName(int i) {
 
235
                switch (i) {
 
236
                        case 0: return A_TO_C;
 
237
                        case 1: return A_TO_G;
 
238
                        case 2: return A_TO_T;
 
239
                        case 3: return C_TO_G;
 
240
                        case 4: return C_TO_T;
 
241
                        default: return UNKNOWN;
 
242
                }
 
243
        }
 
244
 
 
245
        // Make REV model
 
246
        protected void rebuildRateMatrix(double[][] rate, double[] parameters){
 
247
                this.a = parameters[0];
 
248
                this.b = parameters[1];
 
249
                this.c = parameters[2];
 
250
                this.d = parameters[3];
 
251
                this.e = parameters[4];
 
252
                // Q matrix
 
253
                rate[0][1] = a; rate[0][2] = b; rate[0][3] = c;
 
254
                rate[1][2] = d; rate[1][3] = e;
 
255
                rate[2][3] = 1;
 
256
        }
 
257
}