#include "SoloReadFeature.h"
#include "binarySearch2.h"
#include <math.h>

bool inputFeatureUmi(fstream *strIn, int32 featureType, uint32 &feature, uint32 &umi, array<vector<uint64>,2> &sjAll)
{
    if (!(*strIn >> umi)) //end of file
        return false;
    
    if (featureType==0) {//gene
        *strIn >> feature;
    } else if (featureType==1) {//sj
        uint32 sj[2];
        *strIn >> sj[0] >> sj[1];
        feature=(uint32) binarySearch2(sj[0],sj[1],sjAll[0].data(),sjAll[1].data(),sjAll[0].size());
    };
    
    return true;       
};

void SoloReadFeature::inputRecords(uint32 **cbP, uint32 *cbReadCountExactTotal) 
{
    {//load exact matches
        strU_0->flush();
        strU_0->seekg(0,ios::beg);
        uint32 cb, feature, umi;
        while (inputFeatureUmi(strU_0, featureType, feature, umi, P.sjAll)) {
            *strU_0 >> cb;
            if (feature != (uint32)(-1)){
                cbP[cb][0]=feature;
                cbP[cb][1]=umi;
                cbP[cb]+=2;
                stats.V[stats.nExactMatch]++;
            };
        };
    };
    
    {//1 match
        strU_1->flush();
        strU_1->seekg(0,ios::beg);
        uint32 cb, feature, umi;
        while (inputFeatureUmi(strU_1,featureType, feature, umi, P.sjAll)) {
            *strU_1 >> cb;
            if (cbReadCountExactTotal[cb]>0) {
                if (feature != (uint32)(-1)){
                    cbP[cb][0]=feature;
                    cbP[cb][1]=umi;
                    cbP[cb]+=2;
                };
            } else {
                stats.V[stats.nNoExactMatch]++;
            };
        };
    };
     
    {//2 matches
        strU_2->flush();
        strU_2->seekg(0,ios::beg);
        uint32 cb=0, feature, umi, ncb;
        while (inputFeatureUmi(strU_2,featureType, feature, umi, P.sjAll)) {
            if (feature == (uint32) (-1)) {
                strU_2->ignore((uint32) (-1),'\n');//ignore until the end of the line
                continue; //nothing to record
            };
            *strU_2 >> ncb;
            float ptot=0.0,pmax=0.0;
            for (uint32 ii=0; ii<ncb; ii++) {
                uint32 cbin;
                char  qin;
                float pin;
                *strU_2 >> cbin >> qin;
                if (cbReadCountExactTotal[cbin]>0) {//otherwise this cbin does not work
                    qin -= pSolo.QSbase;
                    qin = qin < pSolo.QSmax ? qin : pSolo.QSmax;
                    pin=cbReadCountExactTotal[cbin]*pow(10.0,-qin/10.0);
                    ptot+=pin;
                    if (pin>pmax) {
                        cb=cbin;
                        pmax=pin;
                    };
                };
            };
            if (ptot>0.0 && pmax>=pSolo.cbMinP*ptot) {
                cbP[cb][0]=feature;
                cbP[cb][1]=umi;
                cbP[cb]+=2;
            } else {
                stats.V[stats.nTooMany]++;
            };        
        };
    };
};