NCL: nxscxxdiscretematrix.cpp Source File

00001 //  Copyright (C) 2008 Mark Holder
00002 //
00003 //  This file is part of NCL (Nexus Class Library) version 2.1
00004 //
00005 //  NCL is free software; you can redistribute it and/or modify
00006 //  it under the terms of the GNU General Public License as published by
00007 //  the Free Software Foundation; either version 2 of the License, or
00008 //  (at your option) any later version.
00009 //
00010 //  NCL is distributed in the hope that it will be useful,
00011 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00013 //  GNU General Public License for more details.
00014 //
00015 //  You should have received a copy of the GNU General Public License
00016 //  along with NCL; if not, write to the Free Software Foundation, Inc.,
00017 //  59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00018 //
00019 
00020 #include "ncl/nxscxxdiscretematrix.h"
00021 #include "ncl/nxsutilcopy.h"
00022 #include <cassert>
00023 using std::string;
00024 using std::vector;
00025 using std::cout;
00026 using std::endl;
00027 
00028 
00029 NxsCXXDiscreteMatrix::NxsCXXDiscreteMatrix(const NxsCharactersBlock & cb, bool gapsToMissing)
00030     {
00031     Initialize(&cb, gapsToMissing);
00032     }
00033 
00034 void NxsCXXDiscreteMatrix::Initialize(const NxsCharactersBlock * cb, bool gapsToMissing)
00035 {
00036     this->nativeCMatrix.stateList = 0L;
00037     this->nativeCMatrix.stateListPos = 0L;
00038     this->nativeCMatrix.matrix = 0L;
00039     this->nativeCMatrix.symbolsList = 0L;
00040     this->nativeCMatrix.nStates = 0;
00041     this->nativeCMatrix.nChar = 0;
00042     this->nativeCMatrix.nTax = 0L;
00043     this->nativeCMatrix.nObservedStateSets = 0;
00044     this->nativeCMatrix.datatype = NxsAltGeneric_Datatype;
00045     this->symbolsStringAlias.clear();
00046     this->matrixAlias.Initialize(0, 0);
00047     this->stateListAlias.clear();
00048     this->stateListPosAlias.clear();
00049     this->intWts.clear();
00050     this->dblWts.clear();
00051     this->activeExSet.clear();
00052     if (cb == NULL)
00053         return;
00054     std::vector<const NxsDiscreteDatatypeMapper *> mappers = cb->GetAllDatatypeMappers();
00055     if (mappers.size() > 1)
00056         throw NxsException("too many mappers");
00057     if (mappers.empty() || mappers[0] == NULL)
00058         throw NxsException("no mappers");
00059 
00060     const NxsDiscreteDatatypeMapper & mapper = *mappers[0];
00061     const NxsDiscreteStateMatrix & rawMatrix = cb->GetRawDiscreteMatrixRef();
00062 
00063     NxsCharactersBlock::DataTypesEnum inDatatype = mapper.GetDatatype();
00064     if (inDatatype < LowestNxsCDatatype || inDatatype > HighestNxsCDatatype)
00065         throw NxsException("Datatype cannot be converted to NxsCDiscreteMatrix");
00066     this->nativeCMatrix.datatype = NxsAltDatatypes(inDatatype);
00067     this->nativeCMatrix.nStates = mapper.GetNumStates();
00068     const std::string fundamentalSymbols = mapper.GetSymbols();
00069     const std::string fundamentalSymbolsPlusGaps = mapper.GetSymbolsWithGapChar();
00070     const bool hadGaps = !(fundamentalSymbols == fundamentalSymbolsPlusGaps);
00071 
00072     this->symbolsStringAlias = fundamentalSymbols;
00073     char missingSym = cb->GetMissingSymbol();
00074     NxsCDiscreteState_t newMissingStateCode = this->nativeCMatrix.nStates;
00075     NCL_ASSERT((int)NXS_MISSING_CODE < 0);
00076     NCL_ASSERT((int)NXS_GAP_STATE_CODE < 0);
00077     NxsDiscreteStateCell sclOffsetV;
00078     if (hadGaps)
00079         sclOffsetV = std::min((NxsDiscreteStateCell)NXS_GAP_STATE_CODE, (NxsDiscreteStateCell)NXS_MISSING_CODE);
00080     else
00081         sclOffsetV = NXS_MISSING_CODE;
00082     const NxsDiscreteStateCell sclOffset(sclOffsetV);
00083 
00084     const NxsDiscreteStateCell negSCLOffset = -sclOffset;
00085     const unsigned nMapperStateCodes = mapper.GetNumStateCodes();
00086     const unsigned recodeVecLen = nMapperStateCodes;
00087     const unsigned nMapperPosStateCodes = nMapperStateCodes + sclOffset;
00088     std::vector<NxsCDiscreteState_t> recodeVec(recodeVecLen + negSCLOffset, -2);
00089     NxsCDiscreteState_t * recodeArr = &recodeVec[negSCLOffset];
00090 
00091     if (fundamentalSymbols.length() < this->nativeCMatrix.nStates)
00092         throw NxsException("Fundamental states missing from the symbols string");
00093     const unsigned nfun_sym = (const unsigned)fundamentalSymbols.length();
00094     for (NxsCDiscreteState_t i = 0; i < (NxsCDiscreteState_t)this->nativeCMatrix.nStates; ++i)
00095         {
00096         if (i < (NxsCDiscreteState_t)nfun_sym && (NxsCDiscreteState_t)fundamentalSymbols[i] == '\0' && mapper.PositionInSymbols(fundamentalSymbols[i]) != (NxsDiscreteStateCell) i)
00097             {
00098             NCL_ASSERT(i >= (NxsCDiscreteState_t)nfun_sym || fundamentalSymbols[i] == '\0' || mapper.PositionInSymbols(fundamentalSymbols[i]) == (NxsDiscreteStateCell) i);
00099             }
00100 #       if defined (NDEBUG)
00101             const std::set<NxsDiscreteStateCell>     & ss =  mapper.GetStateSetForCode(i);
00102             NCL_ASSERT(ss.size() == 1);
00103             NCL_ASSERT(*ss.begin() == i);
00104 #       endif
00105         stateListAlias.push_back(1);
00106         stateListAlias.push_back(i);
00107         stateListPosAlias.push_back((unsigned) 2*i);
00108         recodeArr[i] = i;
00109         }
00110 
00111     //NXS_INVALID_STATE_CODE
00112 
00113     if (hadGaps)
00114         recodeArr[NXS_GAP_STATE_CODE] = ((hadGaps && gapsToMissing)? newMissingStateCode : -1);
00115 
00116     if (missingSym == '\0')
00117         missingSym = (hadGaps ? mapper.GetGapSymbol() : '?');
00118     else
00119         {
00120         NCL_ASSERT(NXS_MISSING_CODE == mapper.GetStateCodeStored(missingSym));
00121         }
00122     recodeArr[NXS_MISSING_CODE] = newMissingStateCode;
00123     this->symbolsStringAlias.append(1, missingSym);
00124     const unsigned nCodesInMissing  = this->nativeCMatrix.nStates + (gapsToMissing ?  0 : 1);
00125     stateListPosAlias.push_back(2*this->nativeCMatrix.nStates);
00126     stateListAlias.push_back(nCodesInMissing);
00127     if (!gapsToMissing)
00128         stateListAlias.push_back(-1);
00129     for (NxsCDiscreteState_t i = 0; i < (NxsCDiscreteState_t)this->nativeCMatrix.nStates; ++i)
00130         stateListAlias.push_back(i);
00131 
00132     NxsCDiscreteState_t nextStateCode = newMissingStateCode + 1;
00133     for (NxsDiscreteStateCell i = (NxsDiscreteStateCell)this->nativeCMatrix.nStates; i < (NxsDiscreteStateCell) nMapperPosStateCodes; ++i)
00134         {
00135         const std::set<NxsDiscreteStateCell>     &ss = mapper.GetStateSetForCode( i);
00136         const unsigned ns = (const unsigned)ss.size();
00137         const bool mapToMissing  = (!mapper.IsPolymorphic(i) && (nCodesInMissing + 1 == ns || nCodesInMissing == ns));
00138         if (mapToMissing)
00139             recodeArr[i] = newMissingStateCode;
00140         else
00141             {
00142             recodeArr[i] = nextStateCode++;
00143             stateListPosAlias.push_back((unsigned)stateListAlias.size());
00144             stateListAlias.push_back(ns);
00145             for (std::set<NxsDiscreteStateCell>::const_iterator sIt = ss.begin(); sIt != ss.end(); ++sIt)
00146                 stateListAlias.push_back((NxsCDiscreteState_t) *sIt);
00147             std::string stateName = mapper.StateCodeToNexusString(i);
00148             if (stateName.length() != 1)
00149                 this->symbolsStringAlias.append(1, ' ');
00150             else
00151                 this->symbolsStringAlias.append(1, stateName[0]);
00152             }
00153         }
00154     NCL_ASSERT(stateListPosAlias.size() == (unsigned)nextStateCode);
00155     NCL_ASSERT(symbolsStringAlias.size() == (unsigned)nextStateCode);
00156     this->nativeCMatrix.nObservedStateSets = nextStateCode;
00157 
00158     this->nativeCMatrix.nTax = (unsigned)rawMatrix.size();
00159     this->nativeCMatrix.nChar = (this->nativeCMatrix.nTax == 0 ? 0 : (unsigned)rawMatrix[0].size());
00160     this->matrixAlias.Initialize(this->nativeCMatrix.nTax, this->nativeCMatrix.nChar);
00161     nativeCMatrix.matrix = matrixAlias.GetAlias();
00162     const unsigned nt = this->nativeCMatrix.nTax;
00163     const unsigned nc = this->nativeCMatrix.nChar;
00164     for (unsigned r = 0; r < nt; ++r)
00165         {
00166         NxsCDiscreteStateSet     * recodedRow = nativeCMatrix.matrix[r];
00167         const std::vector<NxsDiscreteStateCell> & rawRowVec = rawMatrix[r];
00168         if (rawRowVec.empty())
00169             {
00170             NxsCDiscreteState_t recodedMissing = recodeArr[NXS_MISSING_CODE];
00171             for (unsigned c = 0; c < nc; ++c)
00172                 *recodedRow++ = recodedMissing;
00173             }
00174         else
00175             {
00176             NCL_ASSERT(rawRowVec.size() == nc);
00177             const NxsDiscreteStateCell * rawRow = &rawRowVec[0];
00178             for (unsigned c = 0; c < nc; ++c)
00179                 {
00180                 const NxsDiscreteStateCell rawC = *rawRow++;
00181                 if ((unsigned)(rawC +  negSCLOffset) >= recodeVecLen)
00182                     {
00183                     NCL_ASSERT((unsigned)(rawC +  negSCLOffset) < recodeVecLen);
00184                     }
00185                 NCL_ASSERT(rawC >= sclOffset);
00186                 const NxsCDiscreteState_t recodedC = recodeArr[rawC];
00187                 NCL_ASSERT(recodedC > -2);
00188                 NCL_ASSERT(recodedC < nextStateCode);
00189                 *recodedRow++ = recodedC;
00190                 }
00191             }
00192         }
00193     nativeCMatrix.symbolsList = symbolsStringAlias.c_str();
00194     nativeCMatrix.stateListPos = &stateListPosAlias[0];
00195     nativeCMatrix.stateList = &stateListAlias[0];
00196 
00197     intWts.clear();
00198     dblWts.clear();
00199     const NxsTransformationManager &tm = cb->GetNxsTransformationManagerRef();
00200     intWts = tm.GetDefaultIntWeights();
00201     if (intWts.empty())
00202         dblWts = tm.GetDefaultDoubleWeights();
00203     activeExSet = cb->GetExcludedIndexSet();
00204 }
00205 
00210 NxsCXXDiscreteMatrix::NxsCXXDiscreteMatrix(const NxsCDiscreteMatrix & mat)
00211     :nativeCMatrix(mat),//aliases pointers, but we'll fix this below
00212     symbolsStringAlias(mat.symbolsList),
00213     matrixAlias(mat.nTax, mat.nChar),
00214     stateListPosAlias(mat.stateListPos, (mat.stateListPos + mat.nObservedStateSets))
00215     {
00216     nativeCMatrix.symbolsList = symbolsStringAlias.c_str();
00217     nativeCMatrix.stateListPos = &stateListPosAlias[0];
00218     if (mat.nObservedStateSets > 0)
00219         {
00220         const unsigned lastStateIndex = nativeCMatrix.stateListPos[nativeCMatrix.nObservedStateSets - 1];
00221         const unsigned lenAmbigList = lastStateIndex + mat.stateList[lastStateIndex] + 1;
00222         //  cout << "lenAmbigList = "<< lenAmbigList <<endl;
00223         stateListAlias.reserve(lenAmbigList);
00224         ncl_copy(mat.stateList, (mat.stateList + lenAmbigList), std::back_inserter(stateListAlias));
00225         }
00226     nativeCMatrix.stateList = &stateListAlias[0];
00227     nativeCMatrix.matrix = matrixAlias.GetAlias();
00228 
00229     // cout << "Matrix in NxsCXXDiscreteMatrix ctor:" << mat.nTax << ' '<< mat.nChar<< endl;
00230     for (unsigned i = 0; i < mat.nTax; ++i)
00231         {
00232         if (mat.nChar > 0)
00233             ncl_copy(mat.matrix[i], mat.matrix[i] + mat.nChar, nativeCMatrix.matrix[i]);
00234         }
00235 
00236     }
00237