-
Notifications
You must be signed in to change notification settings - Fork 3
/
main.c
128 lines (109 loc) · 3.68 KB
/
main.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
/*=================================================================*/
/*= main.c =*/
/*= main file for treekin =*/
/*= --------------------------------------------------------- =*/
/*= Last changed Time-stamp: <2006-11-27 13:09:09 mtw> =*/
/*= $Id: main.c,v 1.24 2006/11/27 13:47:57 mtw Exp $ =*/
/*= --------------------------------------------------------- =*/
/*= (c) Michael Thomas Wolfinger =*/
/*= mtw@tbi.univie.ac.at =*/
/*= treekin =*/
/*=================================================================*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <getopt.h>
#include <errno.h>
#include <math.h>
#include <time.h>
#include "calc.h"
#include "mxccm.h"
#include "barparser.h"
#include "globals.h"
int
main (int argc, char **argv)
{
clock_t clck1 = clock();
BarData *Data=NULL;
/* U - matrix (Q+I)^T, where Q is infetisimal generator (^T - transposed)
S - eigenvectors of U
p0 - distribution in the begining
p8 - stable (end) distribution
R - tmp rates matrix (transposed) */
double *U=NULL, *S=NULL, *p0=NULL, *p8=NULL, *R=NULL;
int dim;
parse_commandline(argc, argv);
switch (opt.method) {
case 'F': dim = ParseInfile(opt.INFILE, &R); break;
case 'L': dim = ParseLocfile(opt.INFILE, &R); break;
case 'I':
dim = ParseBarfile (opt.INFILE, &Data);
ParseRatesFile(&R, dim);
break;
case 'A': dim = ParseBarfile (opt.INFILE, &Data); break;
}
MxInit (dim);
U = MxBar2Matrix(Data, R);
MxGetSpace(&p8);
MxStartVec(&p0);
fprintf(stderr, "Time to initialize: %.2f secs.\n", (clock() - clck1)/(double)CLOCKS_PER_SEC);
clck1 = clock();
if(opt.method == 'F')
MxEqDistrFULL (E, p8);
else {
//MxEqDistr (Data, &p8); /* FIX THIS */
//MxEqDistrFromLocalBalance((TESTING?uU:U), &p8);
MxEqDistrFromLinSys(U, &p8);
}
clck1 = clock();
// first passage times computation
if(opt.fpt) {
// output file?
FILE *fpt = stderr;
if (opt.fpt_file!=NULL) {
fpt = fopen(opt.fpt_file, "w");
if (fpt==NULL) {
fprintf(stderr, "Cannot open file \"%s\" for FPT output, using \"stderr\" instead\n", opt.fpt_file);
fpt = stderr;
}
}
// all or only one state?
if ((opt.fpt_num == -1) || opt.absrb) { // all ftp's
MxFPT(U, p8, fpt);
} else { // only to state opt.ftp_num
double *res = MxFPTOneState(U, opt.fpt_num-1); // starting from 0
if (res != NULL) {
fprintf(fpt, "First passage times to state number %d:\n", opt.fpt_num);
int k;
for (k = 0; k < dim; k++) fprintf(fpt, "%10.7f ", res[k]);
fprintf(fpt,"\n---\n");
free(res);
}
}
fclose(fpt);
fprintf(stderr, "Time to compute fpt: %.2f secs.\n", (clock() - clck1)/(double)CLOCKS_PER_SEC);
clck1 = clock();
}
// iteration
if(opt.matexp) MxExponent(p0,p8,U);
else {
if(opt.rrecover) /* read pre-calculated eigenv{alues,ectors} */
MxRecover(&S, p8);
else {
clck1 = clock();
MxDiagonalize (U, &S, p8);
fprintf(stderr, "Time to diagonalize: %.2f secs.\n", (clock() - clck1)/(double)CLOCKS_PER_SEC);
clck1 = clock();
}
MxIterate (p0, p8, S);
fprintf(stderr, "Time to iterate: %.2f secs.\n", (clock() - clck1)/(double)CLOCKS_PER_SEC);
free(S);
free(p8);
}
MxMemoryCleanUp();
if (opt.pini != NULL) free(opt.pini);
free(U);
if(Data != NULL) free(Data);
return 0;
}
/* End of file */