00001 #pragma ident "$Id: EquationSystem.cpp 2386 2010-04-14 10:03:34Z yanweignss $"
00002
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "EquationSystem.hpp"
00032 #include <iterator>
00033
00034 namespace gpstk
00035 {
00036
00037
00038
00039
00040 WhiteNoiseModel EquationSystem::whiteNoiseModel;
00041
00042
00043
00044
00045
00046
00047
00048
00049 EquationSystem& EquationSystem::addEquation( const Equation& equation )
00050 {
00051
00052
00053 equationDescriptionList.push_back(equation);
00054
00055
00056 isPrepared = false;
00057
00058 return (*this);
00059
00060 }
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073 EquationSystem& EquationSystem::removeEquation( const Variable& indterm )
00074 {
00075
00076
00077 std::list<Equation> backupList;
00078
00079
00080 for( std::list<Equation>::const_iterator itEq =
00081 equationDescriptionList.begin();
00082 itEq != equationDescriptionList.end();
00083 ++itEq )
00084 {
00085
00086
00087 if ( (*itEq).getIndependentTerm() != indterm )
00088 {
00089 backupList.push_back(*itEq);
00090 }
00091
00092 }
00093
00094
00095 clearEquations();
00096
00097
00098 for( std::list<Equation>::const_iterator itEq = backupList.begin();
00099 itEq != backupList.end();
00100 ++itEq )
00101 {
00102 addEquation(*itEq);
00103 }
00104
00105
00106 isPrepared = false;
00107
00108 return (*this);
00109
00110 }
00111
00112
00113
00114
00115 EquationSystem& EquationSystem::clearEquations()
00116 {
00117
00118 equationDescriptionList.clear();
00119
00120 isPrepared = false;
00121
00122 return (*this);
00123
00124 }
00125
00126
00127
00128
00129
00130
00131
00132
00133 EquationSystem& EquationSystem::Prepare( gnssRinex& gData )
00134 {
00135
00136
00137 gnssDataMap myGDSMap;
00138
00139
00140 myGDSMap.addGnssRinex( gData );
00141
00142
00143 return (Prepare(myGDSMap));
00144
00145 }
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 EquationSystem& EquationSystem::Prepare( gnssDataMap& gdsMap )
00156 {
00157
00158
00159 oldUnknowns = currentUnknowns;
00160
00161
00162 varUnknowns = currentUnknowns;
00163
00164
00165 currentUnknowns = prepareCurrentUnknownsAndEquations(gdsMap);
00166
00167
00168 allUnknowns.clear();
00169 for(VariableSet::const_iterator it = currentUnknowns.begin();
00170 it != currentUnknowns.end();
00171 it++)
00172 { allUnknowns.push_back(*it); }
00173
00174 currentUnknowns.clear();
00175 rejectUnknowns.clear();
00176 for(std::list<Variable>::const_iterator it = allUnknowns.begin();
00177 it != allUnknowns.end();
00178 it++)
00179 {
00180 if((*it).getTypeIndexed())
00181 {
00182 currentUnknowns.insert(*it);
00183 }
00184 else
00185 {
00186 rejectUnknowns.insert(*it);
00187 }
00188 }
00189
00190
00191
00192 varUnknowns.insert( currentUnknowns.begin(), currentUnknowns.end() );
00193
00194
00195 getPhiQ(gdsMap);
00196
00197
00198 getPrefit(gdsMap);
00199
00200
00201 getGeometryWeights(gdsMap);
00202
00203
00204 isPrepared = true;
00205
00206 return (*this);
00207
00208 }
00209
00210
00211
00212
00213 void EquationSystem::prepareCurrentSourceSat( gnssDataMap& gdsMap )
00214 {
00215
00216
00217 currentSatSet.clear();
00218 currentSourceSet.clear();
00219
00220
00221 currentSatSet = gdsMap.getSatIDSet();
00222
00223
00224 currentSourceSet = gdsMap.getSourceIDSet();
00225
00226
00227 return;
00228
00229 }
00230
00231
00232
00233
00234 VariableSet EquationSystem::prepareCurrentUnknownsAndEquations(
00235 gnssDataMap& gdsMap )
00236 {
00237
00238
00239 currentEquationsList.clear();
00240
00241
00242 VariableSet currentUnkSet;
00243
00244
00245
00246 prepareCurrentSourceSat( gdsMap );
00247
00248
00249
00250 for( std::list<Equation>::const_iterator itEq =
00251 equationDescriptionList.begin();
00252 itEq != equationDescriptionList.end();
00253 ++itEq )
00254 {
00255
00256
00257 SourceIDSet equSourceSet;
00258
00259
00260 if ( (*itEq).getEquationSource() == Variable::allSources )
00261 {
00262 equSourceSet = currentSourceSet;
00263 }
00264 else
00265 {
00266
00267
00268 if ( (*itEq).getEquationSource() == Variable::someSources )
00269 {
00270
00271
00272
00273 SourceIDSet tempSourceSet( (*itEq).getSourceSet() );
00274
00275
00276
00277 std::insert_iterator< SourceIDSet >
00278 itOut( equSourceSet, equSourceSet.begin() );
00279
00280
00281 set_intersection( tempSourceSet.begin(), tempSourceSet.end(),
00282 currentSourceSet.begin(), currentSourceSet.end(),
00283 itOut );
00284
00285 }
00286 else
00287 {
00288
00289
00290 equSourceSet.insert( (*itEq).getEquationSource() );
00291 }
00292
00293 }
00294
00295
00296 SatIDSet equSatSet = (*itEq).getSatSet();
00297
00298
00299
00300
00301
00302
00303
00304 for( SourceIDSet::const_iterator itSource = equSourceSet.begin();
00305 itSource != equSourceSet.end();
00306 ++itSource )
00307 {
00308
00309
00310 SatIDSet visibleSatSet;
00311
00312
00313 for( gnssDataMap::const_iterator it = gdsMap.begin();
00314 it != gdsMap.end();
00315 ++it )
00316 {
00317
00318
00319 sourceDataMap::const_iterator sdmIter(
00320 (*it).second.find( (*itSource) ) );
00321
00322
00323 if( sdmIter != (*it).second.end() )
00324 {
00325
00326
00327 for( satTypeValueMap::const_iterator stvmIter =
00328 (*sdmIter).second.begin();
00329 stvmIter != (*sdmIter).second.end();
00330 stvmIter++ )
00331 {
00332
00333 if((equSatSet.size() > 0) &&
00334 (equSatSet.find((*stvmIter).first) == equSatSet.end()))
00335 {
00336 continue;
00337 }
00338
00339
00340 visibleSatSet.insert( (*stvmIter).first );
00341
00342 }
00343
00344 }
00345
00346 }
00347
00348
00349
00350
00351
00352 Equation tempEquation( (*itEq) );
00353
00354
00355 tempEquation.clear();
00356
00357
00358 tempEquation.header.equationSource = (*itSource);
00359
00360
00361 for( VariableSet::const_iterator itVar = (*itEq).body.begin();
00362 itVar != (*itEq).body.end();
00363 ++itVar )
00364 {
00365
00366
00367 Variable var( (*itVar) );
00368
00369
00370
00371
00372 if( var.getSourceIndexed() )
00373 {
00374 var.setSource( (*itSource) );
00375 }
00376
00377
00378
00379
00380 tempEquation.addVariable(var);
00381
00382
00383
00384 if( !var.getSatIndexed() )
00385 {
00386
00387
00388 currentUnkSet.insert(var);
00389
00390 }
00391 else
00392 {
00393
00394
00395
00396 for( SatIDSet::const_iterator itSat = visibleSatSet.begin();
00397 itSat != visibleSatSet.end();
00398 ++itSat )
00399 {
00400
00401
00402 var.setSatellite( (*itSat) );
00403
00404
00405
00406 currentUnkSet.insert(var);
00407 }
00408
00409 }
00410
00411 }
00412
00413
00414
00415
00416
00417
00418 for( SatIDSet::const_iterator itSat = visibleSatSet.begin();
00419 itSat != visibleSatSet.end();
00420 ++itSat )
00421 {
00422 tempEquation.header.equationSat = (*itSat);
00423
00424
00425 currentEquationsList.push_back( tempEquation );
00426 }
00427
00428 }
00429
00430 }
00431
00432
00433
00434
00435 int eqListSize( currentEquationsList.size() );
00436 for( int i = 0; i < eqListSize; ++i )
00437 {
00438
00439
00440 Equation tempEqu( currentEquationsList.front() );
00441
00442
00443 currentEquationsList.pop_front();
00444
00445
00446 VariableSet varSet( tempEqu.body );
00447
00448
00449 tempEqu.clear();
00450
00451
00452
00453 for( VariableSet::iterator itVar = varSet.begin();
00454 itVar != varSet.end();
00455 ++itVar )
00456 {
00457
00458
00459 if( !(*itVar).getSatIndexed() )
00460 {
00461
00462 tempEqu.addVariable( (*itVar) );
00463 }
00464 else
00465 {
00466
00467
00468 Variable var( (*itVar) );
00469 var.setSatellite( tempEqu.header.equationSat );
00470
00471 tempEqu.addVariable( var );
00472 }
00473
00474 }
00475
00476
00477 currentEquationsList.push_back( tempEqu );
00478
00479 }
00480
00481
00482
00483 return currentUnkSet;
00484
00485 }
00486
00487
00488
00489
00490 void EquationSystem::getPhiQ( const gnssDataMap& gdsMap )
00491 {
00492
00493 const int numVar( varUnknowns.size() );
00494
00495
00496 phiMatrix.resize( numVar, numVar, 0.0);
00497 qMatrix.resize( numVar, numVar, 0.0);
00498
00499
00500 int i(0);
00501
00502
00503 for( VariableSet::const_iterator itVar = varUnknowns.begin();
00504 itVar != varUnknowns.end();
00505 ++itVar )
00506 {
00507
00508
00509 if( currentUnknowns.find( (*itVar) ) != currentUnknowns.end() )
00510 {
00511
00512
00513 gnssRinex gRin( gdsMap.getGnssRinex( (*itVar).getSource() ) );
00514
00515
00516 (*itVar).getModel()->Prepare( (*itVar).getType(),
00517 (*itVar).getSatellite(),
00518 gRin );
00519
00520
00521 if( oldUnknowns.find( (*itVar) ) != oldUnknowns.end() )
00522 {
00523
00524 phiMatrix(i,i) = (*itVar).getModel()->getPhi();
00525 qMatrix(i,i) = (*itVar).getModel()->getQ();
00526 }
00527 else
00528 {
00529
00530
00531 phiMatrix(i,i) = 0.0;
00532 qMatrix(i,i) = (*itVar).getInitialVariance();
00533 }
00534
00535 }
00536 else
00537 {
00538
00539
00540 phiMatrix(i,i) = whiteNoiseModel.getPhi();
00541 qMatrix(i,i) = whiteNoiseModel.getQ();
00542 }
00543
00544
00545 ++i;
00546 }
00547
00548
00549 return;
00550
00551 }
00552
00553
00554
00555
00556 void EquationSystem::getPrefit( gnssDataMap& gdsMap )
00557 {
00558
00559
00560 std::vector<double> tempPrefit;
00561
00562
00563 for( std::list<Equation>::const_iterator itEq =
00564 currentEquationsList.begin();
00565 itEq != currentEquationsList.end();
00566 ++itEq )
00567 {
00568
00569
00570 tempPrefit.push_back( gdsMap.getValue( (*itEq).header.equationSource,
00571 (*itEq).header.equationSat,
00572 (*itEq).header.indTerm.getType() ) );
00573
00574
00575 }
00576
00577
00578 measVector = tempPrefit;
00579
00580 return;
00581
00582 }
00583
00584
00585
00586
00587 void EquationSystem::getGeometryWeights( gnssDataMap& gdsMap )
00588 {
00589
00590
00591 hMatrix.resize( measVector.size(), varUnknowns.size(), 0.0);
00592 rMatrix.resize( measVector.size(), measVector.size(), 0.0);
00593
00594
00595 gnssDataMap gds2( gdsMap.frontEpoch() );
00596
00597
00598 int row(0);
00599 for( std::list<Equation>::const_iterator itRow =
00600 currentEquationsList.begin();
00601 itRow != currentEquationsList.end();
00602 ++itRow )
00603 {
00604
00605
00606 SourceID source( (*itRow).header.equationSource );
00607 SatID sat( (*itRow).header.equationSat );
00608
00609
00610
00611 TypeIDSet typeSet;
00612
00613
00614 bool found( false );
00615
00616
00617 for( gnssDataMap::const_iterator itGDS = gds2.begin();
00618 itGDS != gds2.end() && !found;
00619 ++itGDS )
00620 {
00621
00622 sourceDataMap::const_iterator itSDM = (*itGDS).second.find(source);
00623 if( itSDM != (*itGDS).second.end() )
00624 {
00625
00626 typeSet = (*itSDM).second.getTypeID();
00627 found = true;
00628 }
00629 }
00630
00631
00632
00633
00634
00635 if( typeSet.find(TypeID::weight) != typeSet.end() )
00636 {
00637
00638 rMatrix(row,row) = (*itRow).header.constWeight
00639 * gds2.getValue(source, sat, TypeID::weight);
00640 }
00641 else
00642 {
00643
00644 rMatrix(row,row) = (*itRow).header.constWeight;
00645 }
00646
00647
00648 int col(0);
00649 for( VariableSet::const_iterator itCol = varUnknowns.begin();
00650 itCol != varUnknowns.end();
00651 ++itCol )
00652 {
00653
00654
00655
00656 if( (*itRow).body.find( (*itCol) ) != (*itRow).body.end() &&
00657 currentUnknowns.find( (*itCol) ) != currentUnknowns.end() )
00658 {
00659
00660
00661
00662 if( (*itCol).isDefaultForced() )
00663 {
00664
00665 hMatrix(row,col) = (*itCol).getDefaultCoefficient();
00666 }
00667 else
00668 {
00669
00670
00671
00672 TypeID type( (*itCol).getType() );
00673
00674
00675 if( typeSet.find(type) != typeSet.end() )
00676 {
00677
00678 hMatrix(row,col) = gds2.getValue(source, sat, type);
00679 }
00680 else
00681 {
00682
00683
00684 hMatrix(row,col) = (*itCol).getDefaultCoefficient();
00685 }
00686
00687 }
00688
00689 }
00690
00691
00692 ++col;
00693
00694 }
00695
00696
00697 for( VariableSet::const_iterator itCol = (*itRow).body.begin();
00698 itCol != (*itRow).body.end();
00699 ++itCol )
00700 {
00701
00702 VariableSet::const_iterator itr = rejectUnknowns.find( (*itCol) );
00703 if( itr == rejectUnknowns.end() || (*itr).getTypeIndexed()) continue;
00704
00705 Variable var(*itr);
00706
00707 col = 0;
00708 for( VariableSet::const_iterator it = varUnknowns.begin(); it != varUnknowns.end(); it++)
00709 {
00710 if(((*itCol).getType() == (*it).getType()) &&
00711 ((*itCol).getModel() == (*it).getModel()) &&
00712 ((*itCol).getSourceIndexed() == (*it).getSourceIndexed())&&
00713 ((*itCol).getSatIndexed() == (*it).getSatIndexed()) &&
00714 ((*itCol).getSource() == (*it).getSource()) &&
00715 ((*itCol).getSatellite() == (*it).getSatellite())
00716 )
00717 {
00718 break;
00719 }
00720
00721 col++;
00722 }
00723
00724
00725
00726
00727 if( (*itCol).isDefaultForced() )
00728 {
00729
00730 hMatrix(row,col) = (*itCol).getDefaultCoefficient();
00731 }
00732 else
00733 {
00734
00735
00736
00737 TypeID type( (*itCol).getType() );
00738
00739
00740 if( typeSet.find(type) != typeSet.end() )
00741 {
00742
00743 hMatrix(row,col) = gds2.getValue(source, sat, type);
00744 }
00745 else
00746 {
00747
00748
00749 hMatrix(row,col) = (*itCol).getDefaultCoefficient();
00750 }
00751
00752 }
00753
00754 }
00755
00756
00757 ++row;
00758
00759 }
00760
00761
00762 return;
00763
00764 }
00765
00766
00767
00768
00769
00770
00771
00772
00773 int EquationSystem::getTotalNumVariables() const
00774 throw(InvalidEquationSystem)
00775 {
00776
00777
00778 if (!isPrepared)
00779 {
00780 GPSTK_THROW(InvalidEquationSystem("EquationSystem is not prepared"));
00781 }
00782
00783 return varUnknowns.size();
00784
00785 }
00786
00787
00788
00789
00790
00791
00792
00793
00794 VariableSet EquationSystem::getVarUnknowns() const
00795 throw(InvalidEquationSystem)
00796 {
00797
00798
00799 if (!isPrepared)
00800 {
00801 GPSTK_THROW(InvalidEquationSystem("EquationSystem is not prepared"));
00802 }
00803
00804 return varUnknowns;
00805
00806 }
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816 int EquationSystem::getCurrentNumVariables() const
00817 throw(InvalidEquationSystem)
00818 {
00819
00820
00821 if (!isPrepared)
00822 {
00823 GPSTK_THROW(InvalidEquationSystem("EquationSystem is not prepared"));
00824 }
00825
00826 return currentUnknowns.size();
00827
00828 }
00829
00830
00831
00832
00833
00834
00835
00836
00837 VariableSet EquationSystem::getCurrentUnknowns() const
00838 throw(InvalidEquationSystem)
00839 {
00840
00841
00842 if (!isPrepared)
00843 {
00844 GPSTK_THROW(InvalidEquationSystem("EquationSystem is not prepared"));
00845 }
00846
00847 return currentUnknowns;
00848
00849 }
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859 int EquationSystem::getCurrentNumSources() const
00860 throw(InvalidEquationSystem)
00861 {
00862
00863
00864 if (!isPrepared)
00865 {
00866 GPSTK_THROW(InvalidEquationSystem("EquationSystem is not prepared"));
00867 }
00868
00869 return currentSourceSet.size();
00870
00871 }
00872
00873
00874
00875
00876
00877
00878
00879
00880 SourceIDSet EquationSystem::getCurrentSources() const
00881 throw(InvalidEquationSystem)
00882 {
00883
00884
00885 if (!isPrepared)
00886 {
00887 GPSTK_THROW(InvalidEquationSystem("EquationSystem is not prepared"));
00888 }
00889
00890 return currentSourceSet;
00891
00892 }
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902 int EquationSystem::getCurrentNumSats() const
00903 throw(InvalidEquationSystem)
00904 {
00905
00906
00907 if (!isPrepared)
00908 {
00909 GPSTK_THROW(InvalidEquationSystem("EquationSystem is not prepared"));
00910 }
00911
00912 return currentSatSet.size();
00913
00914 }
00915
00916
00917
00918
00919
00920
00921
00922
00923 SatIDSet EquationSystem::getCurrentSats() const
00924 throw(InvalidEquationSystem)
00925 {
00926
00927
00928 if (!isPrepared)
00929 {
00930 GPSTK_THROW(InvalidEquationSystem("EquationSystem is not prepared"));
00931 }
00932
00933 return currentSatSet;
00934
00935 }
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945 Vector<double> EquationSystem::getPrefitsVector() const
00946 throw(InvalidEquationSystem)
00947 {
00948
00949
00950 if (!isPrepared)
00951 {
00952 GPSTK_THROW(InvalidEquationSystem("EquationSystem is not prepared"));
00953 }
00954
00955 return measVector;
00956
00957 }
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967 Matrix<double> EquationSystem::getGeometryMatrix() const
00968 throw(InvalidEquationSystem)
00969 {
00970
00971
00972 if (!isPrepared)
00973 {
00974 GPSTK_THROW(InvalidEquationSystem("EquationSystem is not prepared"));
00975 }
00976
00977 return hMatrix;
00978
00979 }
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989 Matrix<double> EquationSystem::getWeightsMatrix() const
00990 throw(InvalidEquationSystem)
00991 {
00992
00993
00994 if (!isPrepared)
00995 {
00996 GPSTK_THROW(InvalidEquationSystem("EquationSystem is not prepared"));
00997 }
00998
00999 return rMatrix;
01000
01001 }
01002
01003
01004
01005
01006
01007
01008
01009
01010 Matrix<double> EquationSystem::getPhiMatrix() const
01011 throw(InvalidEquationSystem)
01012 {
01013
01014
01015 if (!isPrepared)
01016 {
01017 GPSTK_THROW(InvalidEquationSystem("EquationSystem is not prepared"));
01018 }
01019
01020 return phiMatrix;
01021
01022 }
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032 Matrix<double> EquationSystem::getQMatrix() const
01033 throw(InvalidEquationSystem)
01034 {
01035
01036
01037 if (!isPrepared)
01038 {
01039 GPSTK_THROW(InvalidEquationSystem("EquationSystem is not prepared"));
01040 }
01041
01042 return qMatrix;
01043
01044 }
01045
01046
01047
01048 }