← Back to branch summary

~ubuntu-branches/debian/jessie/ugene/jessie

« back to all changes in this revision

Viewing changes to src/plugins_3rdparty/hmm3/src/tests/uhmmer3SearchTests.cpp

  • Committer: Package Import Robot
  • Author(s): Steffen Moeller
  • Date: 2011-11-02 13:29:07 UTC
  • mfrom: (1.2.1) (3.1.11 natty)
  • Revision ID: package-import@ubuntu.com-20111102132907-o34gwnt0uj5g6hen
Tags: 1.9.8+repack-1
* First release to Debian
  - added README.Debian
  - increased policy version to 3.9.2
  - added URLs for version control system
* Added debug package.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/plugins_3rdparty/hmm3/src/tests/uhmmer3SearchTests.cpp
 
1
/**
 
2
 * UGENE - Integrated Bioinformatics Tools.
 
3
 * Copyright (C) 2008-2011 UniPro <ugene@unipro.ru>
 
4
 * http://ugene.unipro.ru
 
5
 *
 
6
 * This program is free software; you can redistribute it and/or
 
7
 * modify it under the terms of the GNU General Public License
 
8
 * as published by the Free Software Foundation; either version 2
 
9
 * of the License, or (at your option) any later version.
 
10
 *
 
11
 * This program is distributed in the hope that it will be useful,
 
12
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 
14
 * GNU General Public License for more details.
 
15
 *
 
16
 * You should have received a copy of the GNU General Public License
 
17
 * along with this program; if not, write to the Free Software
 
18
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 
19
 * MA 02110-1301, USA.
 
20
 */
 
21
 
 
22
#include "uhmmer3SearchTests.h"
 
23
#include <gobject/uHMMObject.h>
 
24
 
 
25
#include <U2Core/DocumentModel.h>
 
26
#include <U2Core/IOAdapter.h>
 
27
#include <U2Core/AppContext.h>
 
28
#include <U2Core/DNASequenceObject.h>
 
29
#include <U2Core/TextUtils.h>
 
30
 
 
31
#include <QtCore/QList>
 
32
#include <memory>
 
33
 
 
34
namespace U2 {
 
35
 
 
36
/*******************************
 
37
* GTest_GeneralUHMM3Search
 
38
********************************/
 
39
 
 
40
const QString GTest_UHMM3Search::SEQ_DOC_CTX_NAME_TAG        = "seqDoc";
 
41
const QString GTest_UHMM3Search::HMM_FILENAME_TAG            = "hmm";
 
42
const QString GTest_UHMM3Search::HMMSEARCH_TASK_CTX_NAME_TAG = "taskCtxName";
 
43
const QString GTest_UHMM3Search::ALGORITHM_TYPE_OPTION_TAG   = "algo";
 
44
const QString GTest_UHMM3Search::SW_CHUNK_SIZE_OPTION_TAG    = "chunk";
 
45
 
 
46
const QString GTest_UHMM3Search::SEQ_E_OPTION_TAG            = "seqE";
 
47
const QString GTest_UHMM3Search::SEQ_T_OPTION_TAG            = "seqT";
 
48
const QString GTest_UHMM3Search::Z_OPTION_TAG                = "z";
 
49
const QString GTest_UHMM3Search::DOM_E_OPTION_TAG            = "domE";
 
50
const QString GTest_UHMM3Search::DOM_T_OPTION_TAG            = "domT";
 
51
const QString GTest_UHMM3Search::DOM_Z_OPTION_TAG            = "domZ";
 
52
const QString GTest_UHMM3Search::USE_BIT_CUTOFFS_OPTION_TAG  = "ubc";
 
53
const QString GTest_UHMM3Search::INC_SEQ_E_OPTION_TAG        = "incE";
 
54
const QString GTest_UHMM3Search::INC_SEQ_T_OPTION_TAG        = "incT";
 
55
const QString GTest_UHMM3Search::INC_DOM_E_OPTION_TAG        = "incdomE";
 
56
const QString GTest_UHMM3Search::INC_DOM_T_OPTION_TAG        = "incdomT";
 
57
const QString GTest_UHMM3Search::MAX_OPTION_TAG              = "max";
 
58
const QString GTest_UHMM3Search::F1_OPTION_TAG               = "f1";
 
59
const QString GTest_UHMM3Search::F2_OPTION_TAG               = "f2";
 
60
const QString GTest_UHMM3Search::F3_OPTION_TAG               = "f3";
 
61
const QString GTest_UHMM3Search::NOBIAS_OPTION_TAG           = "nobias";
 
62
const QString GTest_UHMM3Search::NONULL2_OPTION_TAG          = "nonull2";
 
63
const QString GTest_UHMM3Search::SEED_OPTION_TAG             = "seed";
 
64
const QString GTest_UHMM3Search::REMOTE_MACHINE_VAR          = "MACHINE";
 
65
 
 
66
static void setDoubleOption( double& num, const QDomElement& el, const QString& optionName, TaskStateInfo& si ) {
 
67
    if( si.hasError() ) {
 
68
        return;
 
69
    }
 
70
    QString numStr = el.attribute( optionName );
 
71
    if( numStr.isEmpty() ) {
 
72
        return;
 
73
    }
 
74
    
 
75
    bool ok = false;
 
76
    double ret = numStr.toDouble( &ok );
 
77
    if( !ok ) {
 
78
        si.setError( QString( "cannot_parse_double_number_from %1. Option: %2" ).arg( numStr ).arg( optionName ) );
 
79
        return;
 
80
    }
 
81
    num = ret;
 
82
}
 
83
 
 
84
static void setUseBitCutoffsOption( int& ret, const QDomElement& el, const QString& opName, TaskStateInfo& si ) {
 
85
    if( si.hasError() ) {
 
86
        return;
 
87
    }
 
88
 
 
89
    QString str = el.attribute( opName ).toLower();
 
90
 
 
91
    if( "ga" == str ) {
 
92
        ret = p7H_GA;
 
93
    } else if( "nc" == str ) {
 
94
        ret = p7H_NC;
 
95
    } else if( "tc" == str ) {
 
96
        ret = p7H_TC;
 
97
    } else if( !str.isEmpty() ) {
 
98
        si.setError( QString( "unrecognized_value_in %1 option" ).arg( opName ) );
 
99
    }
 
100
}
 
101
 
 
102
static void setBooleanOption( int& ret, const QDomElement& el, const QString& opName, TaskStateInfo& si ) {
 
103
    if( si.hasError() ) {
 
104
        return;
 
105
    }
 
106
    QString str = el.attribute( opName ).toLower();
 
107
 
 
108
    if( !str.isEmpty() &&  "n" != str && "no" != str ) {
 
109
        ret = TRUE;
 
110
    } else {
 
111
        ret = FALSE;
 
112
    }
 
113
}
 
114
 
 
115
static void setIntegerOption( int& num, const QDomElement& el, const QString& optionName, TaskStateInfo& si ) {
 
116
    if( si.hasError() ) {
 
117
        return;
 
118
    }
 
119
    QString numStr = el.attribute( optionName );
 
120
    if( numStr.isEmpty() ) {
 
121
        return;
 
122
    }
 
123
 
 
124
    bool ok = false;
 
125
    int ret = numStr.toInt( &ok );
 
126
    if( !ok ) {
 
127
        si.setError( QString( "cannot_parse_integer_number_from %1. Option: %2" ).arg( numStr ).arg( optionName ) );
 
128
        return;
 
129
    }
 
130
    num = ret;
 
131
}
 
132
 
 
133
void GTest_UHMM3Search::setSearchTaskSettings( UHMM3SearchSettings& settings, const QDomElement& el, TaskStateInfo& si ) {
 
134
    setDoubleOption( settings.e,       el, GTest_UHMM3Search::SEQ_E_OPTION_TAG,     si );
 
135
    setDoubleOption( settings.t,       el, GTest_UHMM3Search::SEQ_T_OPTION_TAG,     si );
 
136
    setDoubleOption( settings.z,       el, GTest_UHMM3Search::Z_OPTION_TAG,         si );
 
137
    setDoubleOption( settings.f1,      el, GTest_UHMM3Search::F1_OPTION_TAG,        si );
 
138
    setDoubleOption( settings.f2,      el, GTest_UHMM3Search::F2_OPTION_TAG,        si );
 
139
    setDoubleOption( settings.f3,      el, GTest_UHMM3Search::F3_OPTION_TAG,        si );
 
140
    setDoubleOption( settings.domE,    el, GTest_UHMM3Search::DOM_E_OPTION_TAG,     si );
 
141
    setDoubleOption( settings.domT,    el, GTest_UHMM3Search::DOM_T_OPTION_TAG,     si );
 
142
    setDoubleOption( settings.domZ,    el, GTest_UHMM3Search::DOM_Z_OPTION_TAG,     si );
 
143
    setDoubleOption( settings.incE,    el, GTest_UHMM3Search::INC_SEQ_E_OPTION_TAG, si );
 
144
    setDoubleOption( settings.incT,    el, GTest_UHMM3Search::INC_SEQ_T_OPTION_TAG, si );
 
145
    setDoubleOption( settings.incDomE, el, GTest_UHMM3Search::INC_DOM_E_OPTION_TAG, si );
 
146
    setDoubleOption( settings.incDomT, el, GTest_UHMM3Search::INC_DOM_T_OPTION_TAG, si );
 
147
 
 
148
    setBooleanOption( settings.doMax,        el, GTest_UHMM3Search::MAX_OPTION_TAG,     si );
 
149
    setBooleanOption( settings.noBiasFilter, el, GTest_UHMM3Search::NOBIAS_OPTION_TAG,  si );
 
150
    setBooleanOption( settings.noNull2,      el, GTest_UHMM3Search::NONULL2_OPTION_TAG, si );
 
151
 
 
152
    setIntegerOption( settings.seed, el, GTest_UHMM3Search::SEED_OPTION_TAG, si );
 
153
 
 
154
    setUseBitCutoffsOption( settings.useBitCutoffs, el, GTest_UHMM3Search::USE_BIT_CUTOFFS_OPTION_TAG, si );
 
155
}
 
156
 
 
157
static void setSearchAlgoType( GTest_UHMM3SearchAlgoType& alType, const QString& s ) {
 
158
    QString str = s.toLower();
 
159
 
 
160
    if( "general" == str ) {
 
161
        alType = GENERAL_SEARCH;
 
162
    } else if( "sw" == str ) {
 
163
        alType = SEQUENCE_WALKER_SEARCH;
 
164
    } else {
 
165
        alType = UNKNOWN_SEARCH;
 
166
    }
 
167
}
 
168
 
 
169
static P7_HMM * takeHmmFromDoc( Document * doc ) {
 
170
    assert( NULL != doc );
 
171
    QList< GObject* > objs = doc->getObjects();
 
172
    assert( 1 == objs.size() );
 
173
    UHMMObject * hmmObj = qobject_cast< UHMMObject* >( objs.at( 0 ) );
 
174
    if( NULL == hmmObj ) {
 
175
        return NULL;
 
176
    }
 
177
    return hmmObj->takeHMM();
 
178
}
 
179
 
 
180
void GTest_UHMM3Search::init( XMLTestFormat *tf, const QDomElement& el ) {
 
181
    Q_UNUSED( tf );
 
182
 
 
183
    hmmFilename = el.attribute( HMM_FILENAME_TAG );
 
184
    searchTaskCtxName = el.attribute( HMMSEARCH_TASK_CTX_NAME_TAG );
 
185
 
 
186
    searchTaskToCtx     = NULL;
 
187
    hmm = NULL;
 
188
    
 
189
    seqDocCtxName = el.attribute( SEQ_DOC_CTX_NAME_TAG );
 
190
    setSearchAlgoType( algo, el.attribute( ALGORITHM_TYPE_OPTION_TAG ) );
 
191
    setSearchTaskSettings( settings.inner, el, stateInfo );
 
192
    
 
193
    swChunk = UHMM3SWSearchTask::DEFAULT_CHUNK_SIZE;
 
194
    QString chunkStr = el.attribute(SW_CHUNK_SIZE_OPTION_TAG);
 
195
    if( !chunkStr.isEmpty() ) {
 
196
        bool ok = false;
 
197
        int candidate = chunkStr.toInt(&ok);
 
198
        if(ok && candidate > 0) {
 
199
            swChunk = candidate;
 
200
        }
 
201
    }
 
202
    
 
203
    cleanuped = false;
 
204
    ctxAdded = false;
 
205
    
 
206
    machinePath = env->getVar( REMOTE_MACHINE_VAR );
 
207
    if( !machinePath.isEmpty() ) {
 
208
        algo = SEQUENCE_WALKER_SEARCH;
 
209
    }
 
210
}
 
211
 
 
212
void GTest_UHMM3Search::setAndCheckArgs() {
 
213
    assert( !stateInfo.hasError() );
 
214
    if( hmmFilename.isEmpty() ) {
 
215
        stateInfo.setError( "hmm_filename_is_empty" );
 
216
        return;
 
217
    }
 
218
    hmmFilename = env->getVar( "COMMON_DATA_DIR" ) + "/" + hmmFilename;
 
219
 
 
220
    if( searchTaskCtxName.isEmpty() ) {
 
221
        stateInfo.setError( "task_ctx_name_is_empty" );
 
222
        return;
 
223
    }
 
224
 
 
225
    if( seqDocCtxName.isEmpty() ) {
 
226
        stateInfo.setError( "sequence_document_ctx_name_is_empty" );
 
227
        return;
 
228
    }
 
229
 
 
230
    if( UNKNOWN_SEARCH == algo ) {
 
231
        stateInfo.setError( "unknown_algorithm_type" );
 
232
        return;
 
233
    }
 
234
 
 
235
    Document* seqDoc = getContext<Document>( this, seqDocCtxName );
 
236
    if( NULL == seqDoc ) {
 
237
        stateInfo.setError( QString( "context %1 not found" ).arg( seqDocCtxName ) );
 
238
        return;
 
239
    }
 
240
    QList< GObject* > objsList = seqDoc->findGObjectByType( GObjectTypes::SEQUENCE );
 
241
    if( objsList.isEmpty() ) {
 
242
        stateInfo.setError( "no_dna_sequence_objects_in_document" );
 
243
        return;
 
244
    }
 
245
    DNASequenceObject* seqObj = qobject_cast< DNASequenceObject* >( objsList.first() );
 
246
    if( NULL == seqObj ) {
 
247
        stateInfo.setError( "cannot_cast_to_dna_object" );
 
248
        return;
 
249
    }
 
250
    sequence = seqObj->getDNASequence();
 
251
    if( !sequence.length() ) {
 
252
        stateInfo.setError( "empty_sequence_given" );
 
253
        return;
 
254
    }
 
255
    
 
256
    if( !machinePath.isEmpty() ) {
 
257
        machinePath = env->getVar( "COMMON_DATA_DIR" ) + "/" + machinePath;
 
258
    }
 
259
}
 
260
 
 
261
void GTest_UHMM3Search::prepare() {
 
262
    assert( !hasError() );
 
263
    setAndCheckArgs();
 
264
    if( hasError() ) {
 
265
        return;
 
266
    }
 
267
 
 
268
    switch( algo ) {
 
269
    case GENERAL_SEARCH:
 
270
        searchTaskToCtx = new UHMM3LoadProfileAndSearchTask(settings, hmmFilename, sequence.seq);
 
271
        addSubTask( searchTaskToCtx );
 
272
        break;
 
273
    case SEQUENCE_WALKER_SEARCH:
 
274
        if( machinePath.isEmpty() ) { /* search task on local machine */
 
275
            searchTaskToCtx = new UHMM3SWSearchTask( hmmFilename, sequence, settings, swChunk );
 
276
            addSubTask( searchTaskToCtx );
 
277
        } else { /* search on remote machine */
 
278
            addSubTask( LoadDocumentTask::getDefaultLoadDocTask( hmmFilename ) );
 
279
        }
 
280
        break;
 
281
    default:
 
282
        assert( 0 && "undefined_algorithm_type" );
 
283
    }
 
284
}
 
285
 
 
286
QList< Task* > GTest_UHMM3Search::onSubTaskFinished( Task * sub ) {
 
287
    assert( NULL != sub );
 
288
    QList< Task* > res;
 
289
    LoadDocumentTask * loadHmmTask = qobject_cast<LoadDocumentTask*>( sub );
 
290
    if( NULL == loadHmmTask ) {
 
291
        return res;
 
292
    }
 
293
    if( loadHmmTask->hasError() ) {
 
294
        setError( loadHmmTask->getError() );
 
295
        return res;
 
296
    }
 
297
    
 
298
    hmm = takeHmmFromDoc( loadHmmTask->getDocument() );
 
299
    assert( NULL != hmm );
 
300
    return res;
 
301
}
 
302
 
 
303
Task::ReportResult GTest_UHMM3Search::report() {
 
304
    if( stateInfo.hasError() ) {
 
305
        return ReportResult_Finished;
 
306
    }
 
307
    assert( NULL != searchTaskToCtx );
 
308
    if( !searchTaskToCtx->hasError() && !searchTaskToCtx->isCanceled() ) {
 
309
        addContext( searchTaskCtxName, searchTaskToCtx );
 
310
        ctxAdded = true;
 
311
    }
 
312
    return ReportResult_Finished;
 
313
}
 
314
 
 
315
void GTest_UHMM3Search::cleanup() {
 
316
    if( cleanuped ) {
 
317
        return;
 
318
    }
 
319
    if( ctxAdded ) {
 
320
        removeContext( searchTaskCtxName );
 
321
    }
 
322
    if( NULL != hmm ) {
 
323
        p7_hmm_Destroy( hmm );
 
324
    }
 
325
    cleanuped = true;
 
326
}
 
327
 
 
328
GTest_UHMM3Search::~GTest_UHMM3Search() {
 
329
    if( !cleanuped ) {
 
330
        cleanup();
 
331
    }
 
332
}
 
333
 
 
334
/**************************
 
335
* GTest_GeneralUHMM3SearchCompare
 
336
**************************/
 
337
 
 
338
const int   BUF_SZ      = 2048;
 
339
const char  TERM_SYM    = '\0';
 
340
 
 
341
static void readLine( IOAdapter* io, QByteArray& to, QStringList* tokens = NULL ) {
 
342
    assert( NULL != io );
 
343
    to.clear();
 
344
    QByteArray buf( BUF_SZ, TERM_SYM );
 
345
    bool there = false;
 
346
    int bytes = 0;
 
347
    while( !there ) {
 
348
        int ret = io->readUntil( buf.data(), BUF_SZ, TextUtils::LINE_BREAKS, IOAdapter::Term_Include, &there );
 
349
        if( 0 > ret ) {
 
350
            throw QString( "read_error_occurred" );
 
351
        }
 
352
        if( 0 == ret ) {
 
353
            break;
 
354
        }
 
355
        to.append( QByteArray( buf.data(), ret ) );
 
356
        bytes += ret;
 
357
    }
 
358
    to = to.trimmed();
 
359
    if( 0 == bytes ) {
 
360
        throw QString( "unexpected_end_of_file_found" );
 
361
    }
 
362
 
 
363
    if( NULL != tokens ) {
 
364
        *tokens = QString( to ).split( QRegExp( "\\s+" ), QString::SkipEmptyParts );
 
365
    }
 
366
}
 
367
 
 
368
static QByteArray getNextToken( QStringList& tokens ) {
 
369
    if( tokens.isEmpty() ) {
 
370
        throw QString( "unexpected_end_of_line:token_is_missing" );
 
371
    }
 
372
    return tokens.takeFirst().toAscii();
 
373
}
 
374
 
 
375
static double getDouble( const QByteArray& numStr ) {
 
376
    bool ok = false;
 
377
    double ret = numStr.toDouble( &ok );
 
378
    if( ok ) {
 
379
        return ret;
 
380
    }
 
381
    throw QString( GTest_UHMM3SearchCompare::tr( "cannot_parse_double_number_from_string:%1" ).arg( QString( numStr ) ) );
 
382
}
 
383
 
 
384
static float getFloat( const QByteArray& numStr ) {
 
385
    return (float)getDouble( numStr );
 
386
}
 
387
 
 
388
static bool getSignificance( const QByteArray& str ) {
 
389
    if( "!" == str ) {
 
390
        return true;
 
391
    } else if( "?" == str ) {
 
392
        return false;
 
393
    }
 
394
    throw QString( GTest_UHMM3SearchCompare::tr( "cannot_parse_significance:%1" ).arg( QString( str ) ) );
 
395
}
 
396
 
 
397
static UHMM3SearchSeqDomainResult getDomainRes( QStringList& tokens ) {
 
398
    UHMM3SearchSeqDomainResult res;
 
399
 
 
400
    getNextToken( tokens );
 
401
    res.isSignificant = getSignificance( getNextToken( tokens ) );
 
402
    res.score   = getFloat( getNextToken( tokens ) );
 
403
    res.bias    = getFloat( getNextToken( tokens ) );
 
404
    res.cval    = getDouble( getNextToken( tokens ) );
 
405
    res.ival    = getDouble( getNextToken( tokens ) );
 
406
 
 
407
    int hmmFrom = (int)getFloat( getNextToken( tokens ) );
 
408
    int hmmTo   = (int)getFloat( getNextToken( tokens ) );
 
409
    res.queryRegion = U2Region( hmmFrom, hmmTo - hmmFrom );
 
410
    getNextToken( tokens );
 
411
 
 
412
    int aliFrom = (int)getFloat( getNextToken( tokens ) );
 
413
    int aliTo   = (int)getFloat( getNextToken( tokens ) );
 
414
    res.seqRegion = U2Region( aliFrom - 1, aliTo - aliFrom + 1 );
 
415
    getNextToken( tokens );
 
416
 
 
417
    int envFrom = (int)getFloat( getNextToken( tokens ) );
 
418
    int envTo   = (int)getFloat( getNextToken( tokens ) );
 
419
    res.envRegion = U2Region( envFrom, envTo - envFrom );
 
420
    getNextToken( tokens );
 
421
 
 
422
    res.acc = getDouble( getNextToken( tokens ) );
 
423
    return res;
 
424
}
 
425
 
 
426
const double COMPARE_PERCENT_BORDER = 0.01; // 1 percent
 
427
 
 
428
template<class T>
 
429
static bool compareNumbers( T f1, T f2 ) {
 
430
    bool ret = false;
 
431
    if( 0 == f1 ) {
 
432
        ret = 0 == f2 ? true : f2 < COMPARE_PERCENT_BORDER;
 
433
    } else if( 0 == f2 ) {
 
434
        ret = f1 < COMPARE_PERCENT_BORDER;
 
435
    } else {
 
436
        ret = ( qAbs( f1 - f2 ) ) / f1 < COMPARE_PERCENT_BORDER;
 
437
    }
 
438
 
 
439
    if( !ret ) {
 
440
        qDebug() << "!!! compare numbers mismatch: " << f1 << " and " << f2 << " !!!\n";
 
441
    }
 
442
    
 
443
    return ret;
 
444
}
 
445
 
 
446
void GTest_UHMM3SearchCompare::generalCompareResults( const UHMM3SearchResult& myRes, const UHMM3SearchResult& trueRes, TaskStateInfo& ti ) {
 
447
    const UHMM3SearchCompleteSeqResult& myFull = myRes.fullSeqResult;
 
448
    const UHMM3SearchCompleteSeqResult& trueFull = trueRes.fullSeqResult;
 
449
    
 
450
    if( myFull.isReported != trueFull.isReported ) {
 
451
        ti.setError( QString( "reported_flag_not_matched: %1 and %2" ).arg( myFull.isReported ).arg( trueFull.isReported ) );
 
452
        return;
 
453
    }
 
454
    
 
455
    if( myFull.isReported ) {
 
456
        if( !compareNumbers<float>( myFull.bias, trueFull.bias )   ) { 
 
457
            ti.setError( QString( "full_seq_bias_not_matched: %1 and %2" ).arg( myFull.bias ).arg( trueFull.bias ) );  return; 
 
458
        }
 
459
        if( !compareNumbers<double>( myFull.eval, trueFull.eval )  ) {
 
460
            ti.setError( QString( "full_seq_eval_not_matched: %1 and %2" ).arg( myFull.eval ).arg( trueFull.eval ) );  return; 
 
461
        }
 
462
        if( !compareNumbers<float>( myFull.score, trueFull.score ) ) {
 
463
            ti.setError( QString( "full_seq_score_not_matched: %1 and %2" ).arg( myFull.score ).arg( trueFull.score ) ); return; 
 
464
        }
 
465
        if( !compareNumbers<float>( myFull.expectedDomainsNum, trueFull.expectedDomainsNum ) ) {
 
466
            ti.setError( QString( "full_seq_exp_not_matched: %1 and %2" ).arg( myFull.expectedDomainsNum ).arg( trueFull.expectedDomainsNum ) ); 
 
467
            return; 
 
468
        }
 
469
        if( myFull.reportedDomainsNum != trueFull.reportedDomainsNum ) { 
 
470
            ti.setError( QString( "full_seq_n_not_matched: %1 and %2" ).arg( myFull.reportedDomainsNum ).arg( trueFull.reportedDomainsNum ) );
 
471
            return; 
 
472
        }
 
473
    }
 
474
    
 
475
    const QList< UHMM3SearchSeqDomainResult >& myDoms = myRes.domainResList;
 
476
    const QList< UHMM3SearchSeqDomainResult >& trueDoms = trueRes.domainResList;
 
477
    if( myDoms.size() != trueDoms.size() ) {
 
478
        ti.setError( QString( "domain_res_number_not_matched: %1 and %2" ).arg( myDoms.size() ).arg( trueDoms.size() ) );
 
479
        return;
 
480
    }
 
481
    for( int i = 0; i < myDoms.size(); ++i ) {
 
482
        UHMM3SearchSeqDomainResult myCurDom = myDoms.at( i );
 
483
        UHMM3SearchSeqDomainResult trueCurDom = trueDoms.at( i );
 
484
        if( !compareNumbers<double>( myCurDom.acc, trueCurDom.acc ) )   { 
 
485
            ti.setError( QString( "dom_acc_not_matched: %1 and %2" ).arg( myCurDom.acc ).arg( trueCurDom.acc ) );   return; 
 
486
        }
 
487
        if( !compareNumbers<float>( myCurDom.bias, trueCurDom.bias ) )  { 
 
488
            ti.setError( QString( "dom_bias_not_matched: %1 and %2" ).arg( myCurDom.bias ).arg( trueCurDom.bias ) );  return; 
 
489
        }
 
490
        if( !compareNumbers<double>( myCurDom.cval, trueCurDom.cval ) )  { 
 
491
            ti.setError( QString( "dom_cval_not_matched: %1 and %2" ).arg( myCurDom.cval ).arg( trueCurDom.cval ) );  return; 
 
492
        }
 
493
        if( !compareNumbers<double>( myCurDom.ival, trueCurDom.ival ) )  { 
 
494
            ti.setError( QString( "dom_ival_not_matched: %1 and %2" ).arg( myCurDom.ival ).arg( trueCurDom.ival ) );  return; 
 
495
        }
 
496
        if( !compareNumbers<float>( myCurDom.score, trueCurDom.score ) ) { 
 
497
            ti.setError( QString( "dom_score_not_matched: %1 and %2" ).arg( myCurDom.score ).arg( trueCurDom.score ) ); return; 
 
498
        }
 
499
        if( myCurDom.envRegion != trueCurDom.envRegion ) { 
 
500
            ti.setError( QString( "dom_env_region_not_matched: %1---%2 and %3---%4" ).
 
501
                arg( myCurDom.envRegion.startPos ).arg( myCurDom.envRegion.length ).arg( trueCurDom.envRegion.startPos ).
 
502
                arg( trueCurDom.envRegion.length ) ); return;
 
503
        }
 
504
        if( myCurDom.queryRegion != trueCurDom.queryRegion ) { 
 
505
            ti.setError( QString( "dom_hmm_region_not_matched: %1---%2 and %3---%4" ).
 
506
                arg( myCurDom.queryRegion.startPos ).arg( myCurDom.queryRegion.length ).arg( trueCurDom.queryRegion.startPos ).
 
507
                arg( trueCurDom.queryRegion.length ) ); return;
 
508
        }
 
509
        if( myCurDom.seqRegion != trueCurDom.seqRegion ) {
 
510
            ti.setError( QString( "dom_seq_region_not_matched: %1---%2 and %3---%4" ).
 
511
                arg( myCurDom.seqRegion.startPos ).arg( myCurDom.seqRegion.length ).arg( trueCurDom.seqRegion.startPos ).
 
512
                arg( trueCurDom.seqRegion.length ) ); return;
 
513
        }
 
514
        if( myCurDom.isSignificant != trueCurDom.isSignificant ) { 
 
515
            ti.setError( QString( "dom_sign_not_matched: %1 and %2" ).arg( myCurDom.isSignificant ).arg( trueCurDom.isSignificant ) );
 
516
            return; 
 
517
        }
 
518
    }
 
519
}
 
520
 
 
521
static QList<int>
 
522
findEqualDomain(const QList<UHMM3SWSearchTaskDomainResult>& res, const UHMM3SearchSeqDomainResult & dres, bool compareSeqRegion) {
 
523
    QList<int> diff;
 
524
    for(int i = 0; i < res.size(); ++i) {
 
525
        UHMM3SearchSeqDomainResult dom = res.at(i).generalResult;
 
526
        int count = 0;
 
527
        if( !compareNumbers<double>( dom.acc, dres.acc ) )    { count++; }
 
528
        if( !compareNumbers<float>( dom.bias, dres.bias ) )   { count++; }
 
529
        if( !compareNumbers<double>( dom.cval, dres.cval ) )  { count++; }
 
530
        if( !compareNumbers<double>( dom.ival, dres.ival ) )  { count++; }
 
531
        if( !compareNumbers<float>( dom.score, dres.score ) ) { count++; }
 
532
        if( dom.queryRegion != dres.queryRegion ) { count++; }
 
533
        if( compareSeqRegion && dom.seqRegion != dres.seqRegion ) { count++; }
 
534
        if( compareSeqRegion && dom.envRegion != dres.envRegion ) { count++; }
 
535
        if( dom.isSignificant != dres.isSignificant ) { count++; }
 
536
        diff << count;
 
537
    }
 
538
    return diff;
 
539
}
 
540
 
 
541
static QString seqDomainResult2String(const UHMM3SearchSeqDomainResult & r) {
 
542
    return QString("score=%1, eval=%2, bias=%3, acc=%4, query=%5 seq=%6").arg(r.score).arg(r.ival).arg(r.bias).arg(r.acc).
 
543
        arg(QString("%1..%2").arg(r.queryRegion.startPos).arg(r.queryRegion.endPos())).
 
544
        arg(QString("%1..%2").arg(r.seqRegion.startPos).arg(r.seqRegion.endPos()));
 
545
}
 
546
 
 
547
/* we compare here that every domain of trueResult is included in myResult */
 
548
void
 
549
GTest_UHMM3SearchCompare::swCompareResults( const QList<UHMM3SWSearchTaskDomainResult>& myR, const UHMM3SearchResult& trueR, 
 
550
                                            TaskStateInfo& ti, bool compareSeqRegion ) {
 
551
    int sz = trueR.domainResList.size();
 
552
    int i = 0;
 
553
    for( i = 0; i < sz; ++i ) {
 
554
        const UHMM3SearchSeqDomainResult & trueDom = trueR.domainResList.at(i);
 
555
        if(trueDom.score < 2) {
 
556
            continue;
 
557
        }
 
558
        QList<int> diff = findEqualDomain(myR, trueDom, compareSeqRegion);
 
559
        if(!diff.contains(0)) {
 
560
            int minPos = 0;
 
561
            int min = 1000000;
 
562
            for(int j = 0; j < myR.size(); ++j) {
 
563
                float d = qAbs(myR.at(j).generalResult.score - trueR.domainResList.at(i).score);
 
564
                if( d < min ) {
 
565
                    min = d;
 
566
                    minPos = j;
 
567
                }
 
568
            }
 
569
            if(!myR.isEmpty()) {
 
570
                ti.setError( QString( "Cannot find result #%1: %2, most close result: %3" ).
 
571
                    arg(i).
 
572
                    arg(seqDomainResult2String(trueR.domainResList.at(i))).
 
573
                    arg(seqDomainResult2String(myR.at(minPos).generalResult)));
 
574
            } else {
 
575
                ti.setError( QString( "Cannot find result #%1: %2" ).
 
576
                    arg(i).arg(seqDomainResult2String(trueR.domainResList.at(i))));
 
577
            }
 
578
            return;
 
579
        }
 
580
    }
 
581
}
 
582
 
 
583
const QString GTest_UHMM3SearchCompare::SEARCH_TASK_CTX_NAME_TAG = "searchTask";
 
584
const QString GTest_UHMM3SearchCompare::TRUE_OUT_FILE_TAG        = "trueOut";
 
585
 
 
586
void GTest_UHMM3SearchCompare::init( XMLTestFormat *tf, const QDomElement& el ) {
 
587
    Q_UNUSED( tf );
 
588
 
 
589
    searchTaskCtxName = el.attribute( SEARCH_TASK_CTX_NAME_TAG );
 
590
    trueOutFilename = el.attribute( TRUE_OUT_FILE_TAG );
 
591
}
 
592
 
 
593
void GTest_UHMM3SearchCompare::setAndCheckArgs() {
 
594
    assert( !hasError() );
 
595
 
 
596
    if( searchTaskCtxName.isEmpty() ) {
 
597
        stateInfo.setError( "search_task_ctx_name_is_empty" );
 
598
        return;
 
599
    }
 
600
 
 
601
    if( trueOutFilename.isEmpty() ) {
 
602
        stateInfo.setError( "true_out_filename_is_empty" );
 
603
        return;
 
604
    }
 
605
    trueOutFilename = env->getVar( "COMMON_DATA_DIR" ) + "/" + trueOutFilename;
 
606
 
 
607
    Task* searchTask = getContext<Task>( this, searchTaskCtxName );
 
608
    if( NULL == searchTask ) {
 
609
        stateInfo.setError( tr( "cannot_find_search_task_in_context" ) );
 
610
        return;
 
611
    }
 
612
    
 
613
    generalTask = qobject_cast< UHMM3LoadProfileAndSearchTask* >( searchTask );
 
614
    swTask      = qobject_cast< UHMM3SWSearchTask*  >( searchTask );
 
615
    
 
616
    if( NULL != generalTask ) {
 
617
        algo = GENERAL_SEARCH;
 
618
    } else if (NULL != swTask) {
 
619
        algo = SEQUENCE_WALKER_SEARCH;
 
620
    } else {
 
621
        assert( 0 && "cannot_cast_task_to_search_task" );
 
622
    }
 
623
}
 
624
 
 
625
bool GTest_UHMM3SearchCompare::searchResultLessThan(const UHMM3SearchSeqDomainResult & r1, const UHMM3SearchSeqDomainResult & r2) {
 
626
    if( r1.score == r2.score ) {
 
627
        if(r1.seqRegion == r2.seqRegion) {
 
628
            return &r1 < &r2;
 
629
        }
 
630
        return r1.seqRegion < r2.seqRegion;
 
631
    }
 
632
    return r1.score > r2.score;
 
633
}
 
634
 
 
635
UHMM3SearchResult GTest_UHMM3SearchCompare::getOriginalSearchResult( const QString & filename ) {
 
636
    assert( !filename.isEmpty() );
 
637
    
 
638
    IOAdapterFactory* iof = AppContext::getIOAdapterRegistry()->getIOAdapterFactoryById( BaseIOAdapters::url2io( filename ) );
 
639
    std::auto_ptr< IOAdapter > io( iof->createIOAdapter() );
 
640
    if( NULL == io.get() ) {
 
641
        throw QString( "cannot_create_io_adapter_for_'%1'_file" ).arg( filename );
 
642
    }
 
643
    if( !io->open( filename, IOAdapterMode_Read ) ) {
 
644
        throw QString( "cannot_open_'%1'_file" ).arg( filename );
 
645
    }
 
646
    
 
647
    UHMM3SearchResult res;
 
648
    QByteArray buf;
 
649
    QStringList tokens;
 
650
    bool wasHeader = false;
 
651
    bool wasFullSeqResult = false;
 
652
    readLine( io.get(), buf ); /* the first line. starts with # search or # phmmer */
 
653
    do {
 
654
        readLine( io.get(), buf );
 
655
        if( buf.isEmpty() ) { /* but no error - empty line here */
 
656
            continue;
 
657
        }
 
658
        if( buf.startsWith( "# HMMER 3" ) ) {
 
659
            wasHeader = true;
 
660
            continue;
 
661
        }
 
662
        if( buf.startsWith( "Scores for complete sequences" ) ) {
 
663
            if( !wasHeader ) {
 
664
                throw QString( "hmmer_output_header_is_missing" );
 
665
            }
 
666
            UHMM3SearchCompleteSeqResult& fullSeqRes = res.fullSeqResult;
 
667
            readLine( io.get(), buf );
 
668
            readLine( io.get(), buf );
 
669
            readLine( io.get(), buf );
 
670
            readLine( io.get(), buf, &tokens );
 
671
            if( buf.startsWith( "[No hits detected" ) ) {
 
672
                fullSeqRes.isReported = false;
 
673
                break;
 
674
            } else {
 
675
                fullSeqRes.eval     = getDouble( getNextToken( tokens ) );
 
676
                fullSeqRes.score    = getFloat( getNextToken( tokens ) );
 
677
                fullSeqRes.bias     = getFloat( getNextToken( tokens ) );
 
678
                /* skip best domain res. we will check it later */
 
679
                getNextToken( tokens );getNextToken( tokens );getNextToken( tokens );
 
680
                fullSeqRes.expectedDomainsNum = getFloat( getNextToken( tokens ) );
 
681
                fullSeqRes.reportedDomainsNum = (int)getFloat( getNextToken( tokens ) );
 
682
                fullSeqRes.isReported = true;
 
683
                wasFullSeqResult = true;
 
684
            }
 
685
            continue;
 
686
        }
 
687
        if( buf.startsWith( "Domain annotation for each sequence" ) ) {
 
688
            if( !wasFullSeqResult ) {
 
689
                throw QString( "full_seq_result_is_missing" );
 
690
            }
 
691
            readLine( io.get(), buf );
 
692
            readLine( io.get(), buf );
 
693
            readLine( io.get(), buf );
 
694
            QList< UHMM3SearchSeqDomainResult >& domainResList = res.domainResList;
 
695
            assert( domainResList.isEmpty() );
 
696
 
 
697
            int nDomains = res.fullSeqResult.reportedDomainsNum;
 
698
            int i = 0;
 
699
            for( i = 0; i < nDomains; ++i ) {
 
700
                readLine( io.get(), buf, &tokens );
 
701
                domainResList << getDomainRes( tokens );
 
702
            }
 
703
            break;
 
704
        }
 
705
    } while ( 1 );
 
706
    return res;
 
707
}
 
708
 
 
709
Task::ReportResult GTest_UHMM3SearchCompare::report() {
 
710
    assert( !hasError() );
 
711
    setAndCheckArgs();
 
712
    if( hasError() ) {
 
713
        return ReportResult_Finished;
 
714
    }
 
715
    
 
716
    UHMM3SearchResult trueRes;
 
717
    try {
 
718
        trueRes = getOriginalSearchResult( trueOutFilename );
 
719
    } catch( const QString& ex ) {
 
720
        stateInfo.setError( ex );
 
721
    } catch(...) {
 
722
        stateInfo.setError( "undefined_error_occurred" );
 
723
    }
 
724
 
 
725
    if( hasError() ) {
 
726
        return ReportResult_Finished;
 
727
    }
 
728
 
 
729
    switch( algo ) {
 
730
    case GENERAL_SEARCH:
 
731
        assert( NULL != generalTask );
 
732
        generalCompareResults( generalTask->getResult(), trueRes, stateInfo );
 
733
        break;
 
734
    case SEQUENCE_WALKER_SEARCH:
 
735
        {
 
736
            QList<UHMM3SWSearchTaskDomainResult> result;
 
737
            if( NULL != swTask ) {
 
738
                result = swTask->getResults();
 
739
            } else {
 
740
                assert( false );
 
741
            }
 
742
            qSort(trueRes.domainResList.begin(), trueRes.domainResList.end(), searchResultLessThan);
 
743
            swCompareResults( result, trueRes, stateInfo );
 
744
        }
 
745
        break;
 
746
    default:
 
747
        assert( 0 && "unknown_algo_type" );
 
748
    }
 
749
    
 
750
    return ReportResult_Finished;
 
751
}
 
752
 
 
753
} // U2