source: branches/lib/GDE/MAFFT/mafft-7.055-with-extensions/core/dndpre.c

Last change on this file was 10371, checked in by aboeckma, 12 years ago

updated mafft version. Added extensions (no svn ignore, yet)
File size: 7.6 KB
Line 
1#include "mltaln.h"
2
3#define TEST 0
4
5static int treeout = 0;
6static int maxdist = 1;
7static int nadd = 0;
8
9#ifdef enablemultithread
10typedef struct _jobtable
11{
12    int i; 
13    int j; 
14} Jobtable;
15
16typedef struct _thread_arg
17{
18        int njob;
19        int thread_no;
20        float *selfscore;
21        double **mtx;
22        char **seq;
23        Jobtable *jobpospt;
24        pthread_mutex_t *mutex;
25} thread_arg_t;
26
27void *athread( void *arg )
28{
29        thread_arg_t *targ = (thread_arg_t *)arg;
30        int njob = targ->njob;
31        int thread_no = targ->thread_no;
32        float *selfscore = targ->selfscore;
33        double **mtx = targ->mtx;
34        char **seq = targ->seq;
35        Jobtable *jobpospt = targ->jobpospt;
36
37        int i, j;
38        float ssi, ssj, bunbo;
39        double mtxv;
40
41        if( njob == 1 ) return( NULL );
42       
43        while( 1 )
44        {
45                pthread_mutex_lock( targ->mutex );
46                j = jobpospt->j;
47                i = jobpospt->i;
48                j++;
49//              fprintf( stderr, "\n i=%d, j=%d before check\n", i, j );
50                if( j == njob )
51                {
52//                      fprintf( stderr, "\n j = %d, i = %d, njob = %d\n", j, i, njob );
53                        fprintf( stderr, "%4d/%4d (thread %4d), dndpre\r", i+1, njob, thread_no );
54                        i++;
55                        j = i + 1;
56                        if( i == njob-1 )
57                        {
58//                              fprintf( stderr, "\n i=%d, njob-1=%d\n", i, njob-1 );
59                                pthread_mutex_unlock( targ->mutex );
60                                return( NULL );
61                        }
62                }
63//              fprintf( stderr, "\n i=%d, j=%d after check\n", i, j );
64                jobpospt->j = j;
65                jobpospt->i = i;
66                pthread_mutex_unlock( targ->mutex );
67
68                ssi = selfscore[i];
69                ssj = selfscore[j];
70
71                bunbo = MIN( ssi, ssj );
72                if( bunbo == 0.0 )
73                        mtxv = maxdist;
74                else
75                        mtxv = maxdist * ( 1.0 - (double)naivepairscore11( seq[i], seq[j], penalty ) / bunbo );
76#if 1
77                if( mtxv > 9.0 || mtxv < 0.0 )
78                {
79                        fprintf( stderr, "Distance %d-%d is strange, %f.\n", i, j, mtxv );
80                        exit( 1 );
81                }
82#else // CHUUI!!!  2012/05/16
83                if( mtxv > 2.0 )
84                {
85                        mtxv = 2.0;
86                }
87                if( mtxv < 0.0 )
88                {
89                        fprintf( stderr, "Distance %d-%d is strange, %f.\n", i, j, mtxv );
90                        exit( 1 );
91                }
92#endif
93                mtx[i][j] = mtxv;
94        }
95}
96
97#endif
98
99void arguments( int argc, char *argv[] )
100{
101    int c;
102
103        nadd = 0;
104        nthread = 1;
105        alg = 'X';
106        fmodel = 0;
107        treeout = 0;
108        scoremtx = 1;
109        nblosum = 62;
110        dorp = NOTSPECIFIED;
111        inputfile = NULL;
112        ppenalty = NOTSPECIFIED; //?
113        ppenalty_ex = NOTSPECIFIED; //?
114        poffset = NOTSPECIFIED; //?
115        kimuraR = NOTSPECIFIED;
116        pamN = NOTSPECIFIED;
117
118    while( --argc > 0 && (*++argv)[0] == '-' )
119        {
120        while ( (c = *++argv[0]) )
121                {
122            switch( c )
123            {
124                                case 't':
125                                        treeout = '1';
126                                        break;
127                                case 'D':
128                                        dorp = 'd';
129                                        break;
130                                case 'a':
131                                        fmodel = 1;
132                                        break;
133                                case 'P':
134                                        dorp = 'p';
135                                        break;
136                                case 'K': // Hontou ha iranai. disttbfast.c, tbfast.c to awaserutame.
137                                        break;
138                                case 'I':
139                                        nadd = atoi( *++argv );
140                                        fprintf( stderr, "nadd = %d\n", nadd );
141                                        --argc;
142                                        goto nextoption;
143                                case 'f':
144                                        ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
145                                        --argc;
146                                        goto nextoption;
147                                case 'g':
148                                        ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
149                                        --argc;
150                                        goto nextoption;
151                                case 'h':
152                                        poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
153                                        --argc;
154                                        goto nextoption;
155                                case 'k':
156                                        kimuraR = atoi( *++argv );
157//                                      fprintf( stderr, "kimuraR = %d\n", kimuraR );
158                                        --argc;
159                                        goto nextoption;
160                                case 'b':
161                                        nblosum = atoi( *++argv );
162                                        scoremtx = 1;
163//                                      fprintf( stderr, "blosum %d\n", nblosum );
164                                        --argc;
165                                        goto nextoption;
166                                case 'j':
167                                        pamN = atoi( *++argv );
168                                        scoremtx = 0;
169                                        TMorJTT = JTT;
170                                        fprintf( stderr, "jtt %d\n", pamN );
171                                        --argc;
172                                        goto nextoption;
173                                case 'm':
174                                        pamN = atoi( *++argv );
175                                        scoremtx = 0;
176                                        TMorJTT = TM;
177                                        fprintf( stderr, "TM %d\n", pamN );
178                                        --argc;
179                                        goto nextoption;
180                                case 'i':
181                                        inputfile = *++argv;
182                                        fprintf( stderr, "inputfile = %s\n", inputfile );
183                                        --argc;
184                                        goto nextoption;
185                                case 'M':
186                                        maxdist = atoi( *++argv );
187                                        fprintf( stderr, "maxdist = %d\n", maxdist );
188                                        --argc; 
189                                        goto nextoption;
190                                case 'C':
191                                        nthread = atoi( *++argv );
192                                        fprintf( stderr, "nthread = %d\n", nthread );
193                                        --argc; 
194                                        goto nextoption;
195            }
196                }
197                nextoption:
198                        ;
199        }
200    if( argc != 0 ) 
201        {
202                fprintf( stderr, "options: Check source file !\n" );
203                exit( 1 );
204        }
205}
206
207int main( int argc, char **argv )
208{
209        int i, j, ilim;
210        char **seq;
211        static char **name;
212        static int nlen[M];
213        float *selfscore;
214        double **mtx;
215        double mtxv;
216        FILE *fp;
217        FILE *infp;
218        float ssi, ssj, bunbo;
219
220
221        arguments( argc, argv );
222#ifndef enablemultithread
223        nthread = 0;
224#endif
225
226        if( inputfile )
227        {
228                infp = fopen( inputfile, "r" );
229                if( !infp )
230                {
231                        fprintf( stderr, "Cannot open %s\n", inputfile );
232                        exit( 1 );
233                }
234        }
235        else
236                infp = stdin;
237
238#if 0
239        PreRead( stdin, &njob, &nlenmax );
240#else
241        getnumlen( infp );
242#endif
243        rewind( infp );
244
245        njob -= nadd; // atarashii hairetsu ha mushi
246
247        seq = AllocateCharMtx( njob, nlenmax+1 );
248        name = AllocateCharMtx( njob, B+1 );
249        mtx = AllocateDoubleMtx( njob, njob );
250        selfscore = AllocateFloatVec( njob );
251
252#if 0
253        FRead( stdin, name, nlen, seq );
254#else
255        readData_pointer( infp, name, nlen, seq );
256#endif
257        fclose( infp );
258
259        constants( njob, seq );
260
261#if 0
262        for( i=0; i<njob-1; i++ )
263        {
264                fprintf( stderr, "%4d/%4d\r", i+1, njob );
265                for( j=i+1; j<njob; j++ )
266                        mtx[i][j] = (double)substitution_hosei( seq[i], seq[j] );
267//                      fprintf( stderr, "i=%d,j=%d, l=%d &&&  %f\n", i, j, nlen[0], mtx[i][j] );
268        }
269#else // 061003
270        for( i=0; i<njob; i++ )
271        {
272                selfscore[i] = (float)naivepairscore11( seq[i], seq[i], penalty );
273        }
274#ifdef enablemultithread
275        if( nthread > 0 )
276        {
277                thread_arg_t *targ;
278                Jobtable jobpos;
279                pthread_t *handle;
280                pthread_mutex_t mutex;
281
282                jobpos.i = 0;
283                jobpos.j = 0;
284
285                targ = calloc( nthread, sizeof( thread_arg_t ) );
286                handle = calloc( nthread, sizeof( pthread_t ) );
287                pthread_mutex_init( &mutex, NULL );
288
289                for( i=0; i<nthread; i++ )
290                {
291                        targ[i].thread_no = i;
292                        targ[i].njob = njob;
293                        targ[i].selfscore = selfscore;
294                        targ[i].mtx = mtx;
295                        targ[i].seq = seq;
296                        targ[i].jobpospt = &jobpos;
297                        targ[i].mutex = &mutex;
298
299                        pthread_create( handle+i, NULL, athread, (void *)(targ+i) );
300                }
301
302                for( i=0; i<nthread; i++ )
303                {
304                        pthread_join( handle[i], NULL );
305                }
306                pthread_mutex_destroy( &mutex );
307        }
308        else
309#endif
310        {
311                ilim = njob-1;
312                for( i=0; i<ilim; i++ )
313                {
314                        ssi = selfscore[i];
315                        fprintf( stderr, "%4d/%4d\r", i+1, njob );
316                        for( j=i+1; j<njob; j++ )
317                        {
318                                ssj = selfscore[j];
319                                bunbo = MIN( ssi, ssj );
320                                if( bunbo == 0.0 )
321                                        mtxv = maxdist;
322                                else
323                                        mtxv = maxdist * ( 1.0 - (double)naivepairscore11( seq[i], seq[j], penalty ) / bunbo );
324//                                      mtxv = 1.0 - (double)naivepairscore11( seq[i], seq[j], penalty ) / MIN( selfscore[i], selfscore[j] );
325//                              fprintf( stderr, "i=%d,j=%d, l=%d### %f, score = %f, %f, %f\n", i, j, nlen[0], mtxv, naivepairscore11( seq[i], seq[j], penalty ), ssi, ssj  );
326
327#if 1
328                                if( mtxv > 9.0 || mtxv < 0.0 )
329                                {
330                                        fprintf( stderr, "Distance %d-%d is strange, %f.\n", i, j, mtxv );
331                                        exit( 1 );
332                                }
333#else // CHUUI!!!  2012/05/16
334                                if( mtxv > 2.0 )
335                                {
336                                        mtxv = 2.0;
337                                }
338                                if( mtxv < 0.0 )
339                                {
340                                        fprintf( stderr, "Distance %d-%d is strange, %f.\n", i, j, mtxv );
341                                        exit( 1 );
342                                }
343#endif
344                                mtx[i][j] = mtxv;
345                        }
346                }
347        }
348#endif
349       
350#if TEST
351        for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ ) 
352                fprintf( stdout, "i=%d, j=%d, mtx[][] = %f\n", i, j, mtx[i][j] );
353#endif
354
355        fp = fopen( "hat2", "w" );
356        WriteHat2_pointer( fp, njob, name, mtx );
357        fclose( fp );
358#if 0
359        if( treeout )
360        {
361                int ***topol;
362                double **len;
363
364                topol = AllocateIntCub( njob, 2, njob );
365                len = AllocateDoubleMtx( njob, njob );
366                veryfastsupg_double_outtree( njob, mtx, topol, len );
367        }
368#endif
369        SHOWVERSION;
370        exit( 0 );
371/*
372        res = system( ALNDIR "/spgsdl < hat2"  );
373        if( res ) exit( 1 );
374        else exit( 0 );
375*/
376}
Note: See TracBrowser for help on using the repository browser.