/*
 Custom CPAP/Oximetry Calculations Header
 Copyright (c)2011 Mark Watkins <jedimark@users.sourceforge.net>
 License: GPL
*/

#include <QMutex>
#include <QFile>
#include <QDataStream>
#include <QTextStream>
#include <cmath>
#include <algorithm>

#include "calcs.h"
#include "profiles.h"

extern double round(double number);

bool SearchApnea(Session *session, qint64 time, qint64 dist)
{
    if (session->SearchEvent(CPAP_Obstructive,time,dist)) return true;
    if (session->SearchEvent(CPAP_Apnea,time,dist)) return true;
    if (session->SearchEvent(CPAP_ClearAirway,time,dist)) return true;
    if (session->SearchEvent(CPAP_Hypopnea,time,dist)) return true;

    return false;
}

// Sort BreathPeak by peak index
bool operator<(const BreathPeak & p1, const BreathPeak & p2)
{
    return p1.start < p2.start;
}

//! \brief Filters input to output with a percentile filter with supplied width.
//! \param samples Number of samples
//! \param width number of surrounding samples to consider
//! \param percentile fractional percentage, between 0 and 1
void percentileFilter(EventDataType * input, EventDataType * output, int samples, int width, EventDataType percentile)
{
    if (samples<=0)
        return;

    if (percentile>1)
        percentile=1;

    QVector<EventDataType> buf(width);
    int s,e;
    int z1=width/2;
    int z2=z1+(width % 2);
    int nm1=samples-1;
    //int j;

    // Scan through all of input
    for (int k=0;k < samples;k++) {

        s=k-z1;
        e=k+z2;

        // Cap bounds
        if (s < 0) s=0;
        if (e > nm1) e=nm1;

        //
        int j=0;
        for (int i=s; i < e; i++) {
            buf[j++]=input[i];
        }
        j--;

        EventDataType val=j * percentile;
        EventDataType fl=floor(val);

        // If even percentile, or already max value..
        if ((val==fl) || (j>=width-1)) {
            nth_element(buf.begin(),buf.begin()+j,buf.begin()+width-1);
            val=buf[j];
        } else {
            // Percentile lies between two points, interpolate.
            double v1,v2;
            nth_element(buf.begin(),buf.begin()+j,buf.begin()+width-1);
            v1=buf[j];
            nth_element(buf.begin(),buf.begin()+j+1,buf.begin()+width-1);
            v2=buf[j+1];

            val = v1 + (v2-v1)*(val-fl);
        }
        output[k]=val;
    }
}

void xpassFilter(EventDataType * input, EventDataType * output, int samples, EventDataType weight)
{
    // prime the first value
    output[0]=input[0];

    for (int i=1;i<samples;i++) {
        output[i]=weight*input[i] + (1.0-weight)*output[i-1];
    }
    //output[samples-1]=input[samples-1];
}

FlowParser::FlowParser()
{
    m_session=NULL;
    m_flow=NULL;
    m_filtered=NULL;
    m_gain=1;
    m_samples=0;
    m_startsUpper=true;

    // Allocate filter chain buffers..
    m_filtered=(EventDataType *) malloc(max_filter_buf_size);
}
FlowParser::~FlowParser()
{
    free (m_filtered);
//    for (int i=0;i<num_filter_buffers;i++) {
//        free(m_buffers[i]);
//    }
}
void FlowParser::clearFilters()
{
    m_filters.clear();
}

EventDataType * FlowParser::applyFilters(EventDataType * data, int samples)
{
    EventDataType *in=NULL,*out=NULL;
    if (m_filters.size()==0) {
        //qDebug() << "Trying to apply empty filter list in FlowParser..";
        return NULL;
    }
    for (int i=0;i<num_filter_buffers;i++) {
        m_buffers[i]=(EventDataType *) malloc(max_filter_buf_size);
    }

    int numfilt=m_filters.size();

    for (int i=0;i<numfilt;i++) {
        if (i==0) {
            in=data;
            out=m_buffers[0];
            if (in==out) {
                //qDebug() << "Error: If you need to use internal m_buffers as initial input, use the second one. No filters were applied";
                return NULL;
            }
        } else {
            in=m_buffers[(i+1) % num_filter_buffers];
            out=m_buffers[i % num_filter_buffers];
        }

        // If final link in chain, pass it back out to input data
        if (i==numfilt-1) {
            out=data;
        }

        Filter & filter=m_filters[i];
        if (filter.type==FilterNone) {
            // Just copy it..
            memcpy(in,out,samples*sizeof(EventDataType));
        } else if (filter.type==FilterPercentile) {
            percentileFilter(in,out,samples,filter.param1,filter.param2);
        } else if (filter.type==FilterXPass) {
            xpassFilter(in,out,samples,filter.param1);
        }
    }

    for (int i=0;i<num_filter_buffers;i++) {
        free(m_buffers[i]);
    }
    return out;
}
void FlowParser::openFlow(Session * session, EventList * flow)
{
    if (!flow) {
        qDebug() << "called FlowParser::processFlow() with a null EventList!";
        return;
    }
    m_session=session;
    m_flow=flow;
    m_gain=flow->gain();
    m_rate=flow->rate();
    m_samples=flow->count();
    EventStoreType *inraw=flow->rawData();

    // Make sure we won't overflow internal buffers
    if (m_samples > max_filter_buf_size) {
        qDebug() << "Error: Sample size exceeds max_filter_buf_size in FlowParser::openFlow().. Capping!!!";
        m_samples=max_filter_buf_size;
    }

    // Begin with the second internal buffer
    EventDataType * buf=m_filtered;
    // Apply gain to waveform
    EventStoreType *eptr=inraw+m_samples;

    // Convert from store type to floats..
    for (; inraw < eptr; inraw++) {
        *buf++ = EventDataType(*inraw) * m_gain;
    }

    // Apply the rest of the filters chain
    buf=applyFilters(m_filtered, m_samples);

    calcPeaks(m_filtered, m_samples);
}

// Calculates breath upper & lower peaks for a chunk of EventList data
void FlowParser::calcPeaks(EventDataType * input, int samples)
{
    if (samples<=0)
        return;

    EventDataType min=0,max=0, c, lastc=0;

    EventDataType zeroline=0;

    double rate=m_flow->rate();

    double flowstart=m_flow->first();
    double lasttime,time;

    double peakmax=flowstart, peakmin=flowstart;

    time=lasttime=flowstart;
    breaths.clear();

    // Estimate storage space needed using typical average breaths per minute.
    m_minutes=double(m_flow->last() - m_flow->first()) / 60000.0;
    const double avgbpm=20;
    int guestimate=m_minutes*avgbpm;
    breaths.reserve(guestimate);

    // Prime min & max, and see which side of the zero line we are starting from.
    c=input[0];
    min=max=c;
    lastc=c;
    m_startsUpper=(c >= zeroline);

    qint32 start=0,middle=0;//,end=0;


    int sps=1000/m_rate;
    int len=0, lastk=0; //lastlen=0

    qint64 sttime=time;//, ettime=time;

    // For each samples, find the peak upper and lower breath components
    for (int k=0; k < samples; k++) {
        c=input[k];

        if (c >= zeroline) {

            // Did we just cross the zero line going up?
            if (lastc < zeroline) {
                // This helps filter out dirty breaths..
                len=k-start;
                if ((max>3) && ((max-min) > 8) && (len>sps) && (middle > start))  {

                    // peak detection may not be needed..
                    breaths.push_back(BreathPeak(min, max, start, middle, k)); //, peakmin, peakmax));

                    // Set max for start of the upper breath cycle
                    max=c;
                    peakmax=time;


                    // Starting point of next breath cycle
                    start=k;
                    sttime=time;
                }
            } else if (c > max) {
                // Update upper breath peak
                max=c;
                peakmax=time;
            }
        }

        if (c < zeroline) {
            // Did we just cross the zero line going down?

            if (lastc >= zeroline) {
                // Set min for start of the lower breath cycle
                min=c;
                peakmin=time;
                middle=k;

            } else if (c < min) {
                // Update lower breath peak
                min=c;
                peakmin=time;
            }

        }
        lasttime=time;
        time+=rate;
        lastc=c;
        lastk=k;
    }
}


//! \brief Calculate Respiratory Rate, TidalVolume, Minute Ventilation, Ti & Te..
// These are grouped together because, a) it's faster, and b) some of these calculations rely on others.
void FlowParser::calc(bool calcResp, bool calcTv, bool calcTi, bool calcTe, bool calcMv)
{
    if (!m_session) {
        return;
    }

    // Don't even bother if only a few breaths in this chunk
    const int lowthresh=4;
    int nm=breaths.size();
    if (nm < lowthresh)
        return;

    const qint64 minute=60000;


    double start=m_flow->first();
    double time=start;

    int bs,be,bm;
    double st,et,mt;

    /////////////////////////////////////////////////////////////////////////////////
    // Respiratory Rate setup
    /////////////////////////////////////////////////////////////////////////////////
    EventDataType minrr=0,maxrr=0;
    EventList * RR=NULL;
    quint32 * rr_tptr=NULL;
    EventStoreType * rr_dptr=NULL;
    if (calcResp) {
        RR=m_session->AddEventList(CPAP_RespRate,EVL_Event);
        minrr=RR->Min(), maxrr=RR->Max();
        RR->setGain(0.2);
        RR->setFirst(time+minute);
        RR->getData().resize(nm);
        RR->getTime().resize(nm);
        rr_tptr=RR->rawTime();
        rr_dptr=RR->rawData();
    }

    int rr_count=0;

    double len, st2,et2,adj, stmin, b, rr=0;
    double len2;
    int idx;



    /////////////////////////////////////////////////////////////////////////////////
    // Inspiratory / Expiratory Time setup
    /////////////////////////////////////////////////////////////////////////////////
    double lastte2=0,lastti2=0,lastte=0, lastti=0, te, ti, ti1, te1,c;
    EventList * Te=NULL, * Ti=NULL;
    if (calcTi) {
        Ti=m_session->AddEventList(CPAP_Ti,EVL_Event);
        Ti->setGain(0.02);
    }
    if (calcTe) {
        Te=m_session->AddEventList(CPAP_Te,EVL_Event);
        Te->setGain(0.02);
    }


    /////////////////////////////////////////////////////////////////////////////////
    // Tidal Volume setup
    /////////////////////////////////////////////////////////////////////////////////
    EventList * TV=NULL;
    EventDataType mintv=0, maxtv=0, tv=0;
    double val1, val2;
    quint32 * tv_tptr=NULL;
    EventStoreType * tv_dptr=NULL;
    int tv_count=0;
    if (calcTv) {
        TV=m_session->AddEventList(CPAP_TidalVolume,EVL_Event);
        mintv=TV->Min(), maxtv=TV->Max();
        TV->setGain(20);
        TV->setFirst(start);
        TV->getData().resize(nm);
        TV->getTime().resize(nm);
        tv_tptr=TV->rawTime();
        tv_dptr=TV->rawData();
    }
    /////////////////////////////////////////////////////////////////////////////////
    // Minute Ventilation setup
    /////////////////////////////////////////////////////////////////////////////////
    EventList * MV=NULL;
    EventDataType mv;
    if (calcMv) {
        MV=m_session->AddEventList(CPAP_MinuteVent,EVL_Event);
        MV->setGain(0.125);
    }

    EventDataType sps=(1000.0/m_rate); // Samples Per Second


    qint32 timeval=0; // Time relative to start

    for (idx=0;idx<nm;idx++) {
        bs=breaths[idx].start;
        bm=breaths[idx].middle;
        be=breaths[idx].end;

        // Calculate start, middle and end time of this breath
        st=start+bs * m_rate;
        mt=start+bm * m_rate;

        timeval=be * m_rate;

        et=start+timeval;



        /////////////////////////////////////////////////////////////////////
        // Calculate Inspiratory Time (Ti) for this breath
        /////////////////////////////////////////////////////////////////////
        if (calcTi) {
            ti=((mt-st)/1000.0)*50.0;
            ti1=(lastti2+lastti+ti)/3.0;
            Ti->AddEvent(mt,ti1);
            lastti2=lastti;
            lastti=ti;
        }
        /////////////////////////////////////////////////////////////////////
        // Calculate Expiratory Time (Te) for this breath
        /////////////////////////////////////////////////////////////////////
        if (calcTe) {
            te=((et-mt)/1000.0)*50.0;
            // Average last three values..
            te1=(lastte2+lastte+te)/3.0;
            Te->AddEvent(mt,te1);

            lastte2=lastte;
            lastte=te;
        }
        /////////////////////////////////////////////////////////////////////
        // Calculate TidalVolume
        /////////////////////////////////////////////////////////////////////
        if (calcTv) {
            val1=0, val2=0;
            // Scan the upper breath
            for (int j=bs;j<bm;j++)  {
                // convert flow to ml/s to L/min and divide by samples per second
                c=double(qAbs(m_filtered[j])) * 1000.0 / 60.0 / sps;
                val2+=c;
                //val2+=c*c; // for RMS
            }
            // calculate root mean square
            //double n=bm-bs;
            //double q=(1/n)*val2;
            //double x=sqrt(q)*2;
            //val2=x;
            tv=val2;

            if (tv < mintv) mintv=tv;
            if (tv > maxtv) maxtv=tv;
            *tv_tptr++ = timeval;
            *tv_dptr++ = tv / 20.0;
            tv_count++;

        }

        /////////////////////////////////////////////////////////////////////
        // Respiratory Rate Calculations
        /////////////////////////////////////////////////////////////////////
        if (calcResp) {
            stmin=et-minute;
            if (stmin < start)
                stmin=start;
            len=et-stmin;
            rr=0;
            len2=0;
            //if ( len >= minute) {
                //et2=et;

                // Step back through last minute and count breaths
                for (int i=idx;i>=0;i--) {
                    st2=start + double(breaths[i].start) * m_rate;
                    et2=start + double(breaths[i].end) * m_rate;
                    if (et2 < stmin)
                        break;

                    len=et2-st2;
                    if (st2 < stmin) {
                        // Partial breath
                        st2=stmin;
                        adj=et2 - st2;
                        b=(1.0 / len) * adj;
                        len2+=adj;
                    } else {
                        b=1;
                        len2+=len;
                    }

                    rr+=b;
                }
                if (len2 < minute) {
                    rr*=minute/len2;
                }
                // Calculate min & max
                if (rr < minrr)
                    minrr=rr;
                if (rr > maxrr)
                    maxrr=rr;

                // Add manually.. (much quicker)
                *rr_tptr++ = timeval;

                // Use the same gains as ResMed..

                *rr_dptr++ = rr * 5.0;
                rr_count++;

                //rr->AddEvent(et,br * 50.0);
            //}
        }
        if (calcMv && calcResp && calcTv) {
            mv=(tv/1000.0) * rr;
            MV->AddEvent(et,mv * 8.0);
        }
    }

    /////////////////////////////////////////////////////////////////////
    // Respiratory Rate post filtering
    /////////////////////////////////////////////////////////////////////

    if (calcResp) {
        RR->setMin(minrr);
        RR->setMax(maxrr);
        RR->setFirst(start);
        RR->setLast(et);
        RR->setCount(rr_count);
    }
    /////////////////////////////////////////////////////////////////////
    // Tidal Volume post filtering
    /////////////////////////////////////////////////////////////////////

    if (calcTv) {
        TV->setMin(mintv);
        TV->setMax(maxtv);
        TV->setFirst(start);
        TV->setLast(et);
        TV->setCount(tv_count);
    }
}

void FlowParser::flagEvents()
{
    if (!PROFILE.cpap->userEventFlagging())  return;

    int numbreaths=breaths.size();
    if (numbreaths<5) return;

    EventDataType val,mx,mn;
    QVector<EventDataType> br;

    QVector<qint32> bstart;
    QVector<qint32> bend;
    //QVector<EventDataType> bvalue;

    bstart.reserve(numbreaths*2);
    bend.reserve(numbreaths*2);
    //bvalue.reserve(numbreaths*2);
    br.reserve(numbreaths*2);

    double start=m_flow->first();
   // double sps=1000.0/m_rate;
    double st,et, dur; //mt
    qint64 len;

    bool allowDuplicates=PROFILE.cpap->userEventDuplicates();

    for (int i=0;i<numbreaths;i++) {
        mx=breaths[i].max;
        mn=breaths[i].min;
        br.push_back(qAbs(mx));
        br.push_back(qAbs(mn));
    }
    //EventList * uf2=m_session->AddEventList(CPAP_UserFlag2,EVL_Event);
    //EventList * uf3=m_session->AddEventList(CPAP_UserFlag3,EVL_Event);

    const EventDataType perc=0.6;
    int idx=float(br.size())*perc;
    nth_element(br.begin(),br.begin()+idx,br.end()-1);

    EventDataType peak=br[idx];//*(br.begin()+idx);

    EventDataType cutoffval=peak * (PROFILE.cpap->userFlowRestriction()/100.0);

    int bs,bm,be, bs1, bm1, be1;
    for (int i=0;i<numbreaths;i++) {
        bs=breaths[i].start;
        bm=breaths[i].middle;
        be=breaths[i].end;

        mx=breaths[i].max;
        mn=breaths[i].min;
        val=mx - mn;

//        if (qAbs(mx) > cutoffval) {
            bs1=bs;
            for (;bs1<be;bs1++) {
                if (qAbs(m_filtered[bs1]) > cutoffval) {
                    break;
                }
            }

            bm1=bm;
            for (;bm1>bs;bm1--) {
                if (qAbs(m_filtered[bm1]) > cutoffval) {
                    break;
                }
            }
            if (bm1>=bs1) {
                bstart.push_back(bs1);
                bend.push_back(bm1);
            }
  //      }
    //    if (qAbs(mn) > cutoffval) {
            bm1=bm;
            for (;bm1<be;bm1++) {
                if (qAbs(m_filtered[bm1]) > cutoffval) {
                    break;
                }
            }
            be1=be;
            for (;be1>bm;be1--) {
                if (qAbs(m_filtered[be1]) > cutoffval) {
                    break;
                }
            }
            if (be1>=bm1) {
                bstart.push_back(bm1);
                bend.push_back(be1);
            }
//        }
         st=start + bs1 * m_rate;
         et=start + be1  * m_rate;
//         uf2->AddEvent(st,0);
//         uf3->AddEvent(et,0);

    }


    EventDataType duration=PROFILE.cpap->userEventDuration();
    //double lastst=start, lastet=start;
    //EventDataType v;
    int bsize=bstart.size();
    EventList * uf1=NULL;

    for (int i=0;i<bsize-1;i++) {
        bs=bend[i];
        be=bstart[i+1];
        st=start + bs * m_rate;
        et=start + be  * m_rate;

        len=et-st;
        dur=len/1000.0;
        if (dur>=duration) {
            if (allowDuplicates || !SearchApnea(m_session,et-len/2,15000)) {
                if (!uf1) {
                    uf1=m_session->AddEventList(CPAP_UserFlag1,EVL_Event);
                }
                uf1->AddEvent(et-len/2,dur);
            }
        }
    }
}

void calcRespRate(Session *session, FlowParser * flowparser)
{
    if (session->machine()->GetType()!=MT_CPAP) return;

//    if (session->machine()->GetClass()!=STR_MACH_PRS1) return;

    if (!session->eventlist.contains(CPAP_FlowRate)) {
        //qDebug() << "calcRespRate called without FlowRate waveform available";
        return; //need flow waveform
    }

    bool trashfp;
    if (!flowparser) {
        flowparser=new FlowParser();
        trashfp=true;
        //qDebug() << "calcRespRate called without valid FlowParser object.. using a slow throw-away!";
        //return;
    } else {
        trashfp=false;
    }

    bool calcResp=!session->eventlist.contains(CPAP_RespRate);
    bool calcTv=!session->eventlist.contains(CPAP_TidalVolume);
    bool calcTi=!session->eventlist.contains(CPAP_Ti);
    bool calcTe=!session->eventlist.contains(CPAP_Te);
    bool calcMv=!session->eventlist.contains(CPAP_MinuteVent);


    int z=(calcResp ? 1 : 0) + (calcTv ? 1 : 0) + (calcMv ? 1 : 0);

    // If any of these three missing, remove all, and switch all on
    if (z>0 && z<3) {
        if (!calcResp && !calcTv && !calcMv)
        calcTv=calcMv=calcResp=true;

        QVector<EventList *> & list=session->eventlist[CPAP_RespRate];
        for (int i=0;i<list.size();i++) {
            delete list[i];
        }
        session->eventlist[CPAP_RespRate].clear();

        QVector<EventList *> & list2=session->eventlist[CPAP_TidalVolume];
        for (int i=0;i<list2.size();i++) {
            delete list2[i];
        }
        session->eventlist[CPAP_TidalVolume].clear();

        QVector<EventList *> & list3=session->eventlist[CPAP_MinuteVent];
        for (int i=0;i<list3.size();i++) {
            delete list3[i];
        }
        session->eventlist[CPAP_MinuteVent].clear();
    }

    flowparser->clearFilters();

    // No filters works rather well with the new peak detection algorithm..
    // Although the output could use filtering.

    //flowparser->addFilter(FilterPercentile,7,0.5);
    //flowparser->addFilter(FilterPercentile,5,0.5);
    //flowparser->addFilter(FilterXPass,0.5);
    EventList *flow;
    for (int ws=0; ws < session->eventlist[CPAP_FlowRate].size(); ws++) {
        flow=session->eventlist[CPAP_FlowRate][ws];
        if (flow->count() > 20) {
            flowparser->openFlow(session, flow);
            flowparser->calc(calcResp, calcTv, calcTi ,calcTe, calcMv);
            flowparser->flagEvents();
        }
    }
    if (trashfp) {
        delete flowparser;
    }
}


EventDataType calcAHI(Session *session,qint64 start, qint64 end)
{
    bool rdi=PROFILE.general->calculateRDI();

    double hours,ahi,cnt;
    if (start<0) {
        // much faster..
        hours=session->hours();
        cnt=session->count(CPAP_Obstructive)
                +session->count(CPAP_Hypopnea)
                +session->count(CPAP_ClearAirway)
                +session->count(CPAP_Apnea);

        if (rdi) {
            cnt+=session->count(CPAP_RERA);
        }

        ahi=cnt/hours;
    } else {
        hours=double(end-start)/3600000L;
        if (hours==0) return 0;
        cnt=session->rangeCount(CPAP_Obstructive,start,end)
                +session->rangeCount(CPAP_Hypopnea,start,end)
                +session->rangeCount(CPAP_ClearAirway,start,end)
                +session->rangeCount(CPAP_Apnea,start,end);

        if (rdi) {
            cnt+=session->rangeCount(CPAP_RERA,start,end);
        }

        ahi=cnt/hours;
    }
    return ahi;
}

int calcAHIGraph(Session *session)
{
    bool calcrdi=session->machine()->GetClass()=="PRS1";
    //PROFILE.general->calculateRDI()


    const qint64 window_step=30000; // 30 second windows
    double window_size=p_profile->cpap->AHIWindow();
    qint64 window_size_ms=window_size*60000L;

    bool zeroreset=p_profile->cpap->AHIReset();

    if (session->machine()->GetType()!=MT_CPAP) return 0;

    bool hasahi=session->eventlist.contains(CPAP_AHI);
    bool hasrdi=session->eventlist.contains(CPAP_RDI);
    if (hasahi && hasrdi)
        return 0; // abort if already there

    if (!(!hasahi && !hasrdi)) {
        session->destroyEvent(CPAP_AHI);
        session->destroyEvent(CPAP_RDI);
    }

    if (!session->channelExists(CPAP_Obstructive) &&
            !session->channelExists(CPAP_Hypopnea) &&
            !session->channelExists(CPAP_Apnea) &&
            !session->channelExists(CPAP_ClearAirway) &&
            !session->channelExists(CPAP_RERA)
            ) return 0;

    qint64 first=session->first(),
           last=session->last(),
           f;

    EventList *AHI=new EventList(EVL_Event);
    AHI->setGain(0.02);
    session->eventlist[CPAP_AHI].push_back(AHI);

    EventList *RDI=NULL;

    if (calcrdi) {
        RDI=new EventList(EVL_Event);
        RDI->setGain(0.02);
        session->eventlist[CPAP_RDI].push_back(RDI);
    }

    EventDataType ahi,rdi;

    qint64 ti=first,lastti=first;

    double avgahi=0;
    double avgrdi=0;
    int cnt=0;

    double events;
    double hours=(window_size/60.0);
    if (zeroreset) {
        // I personally don't see the point of resetting each hour.
        do {
            // For each window, in 30 second increments
            for (qint64 t=ti;t < ti+window_size_ms; t+=window_step) {
                if (t > last)
                    break;
                events=session->rangeCount(CPAP_Obstructive,ti,t)
                        +session->rangeCount(CPAP_Hypopnea,ti,t)
                        +session->rangeCount(CPAP_ClearAirway,ti,t)
                        +session->rangeCount(CPAP_Apnea,ti,t);

                ahi = events / hours;

                AHI->AddEvent(t,ahi * 50);
                avgahi+=ahi;

                if (calcrdi) {
                    events+=session->rangeCount(CPAP_RERA,ti,t);
                    rdi=events / hours;
                    RDI->AddEvent(t,rdi * 50);
                    avgrdi+=rdi;
                }

                cnt++;
            }
            lastti=ti;
            ti+=window_size_ms;
        } while (ti<last);

    } else {
        for (ti=first;ti<last;ti+=window_step) {
            f=ti-window_size_ms;
            //hours=window_size; //double(ti-f)/3600000L;

            events=session->rangeCount(CPAP_Obstructive,f,ti)
                    +session->rangeCount(CPAP_Hypopnea,f,ti)
                    +session->rangeCount(CPAP_ClearAirway,f,ti)
                    +session->rangeCount(CPAP_Apnea,f,ti);

            ahi=events/hours;
            avgahi+=ahi;
            AHI->AddEvent(ti,ahi * 50);

            if (calcrdi) {
                events+=session->rangeCount(CPAP_RERA,f,ti);
                rdi=events/hours;
                RDI->AddEvent(ti,rdi * 50);
                avgrdi+=rdi;
            }

            cnt++;
            lastti=ti;
            ti+=window_step;
        }
    }
    AHI->AddEvent(lastti,0);
    if (calcrdi)
        RDI->AddEvent(lastti,0);

    if (!cnt) {
            avgahi=0;
            avgrdi=0;
    } else {
        avgahi/=double(cnt);
        avgrdi/=double(cnt);
    }
    cnt++;
    session->setAvg(CPAP_AHI,avgahi);
    if (calcrdi)
        session->setAvg(CPAP_RDI,avgrdi);
    return cnt;
}
struct TimeValue {
    TimeValue() {
        time=0;
        value=0;
    }
    TimeValue(qint64 t, EventStoreType v) {
        time=t;
        value=v;
    }
    TimeValue(const TimeValue & copy) {
        time=copy.time;
        value=copy.value;
    }
    qint64 time;
    EventStoreType value;
};

struct zMaskProfile {
public:
    zMaskProfile(MaskType type, QString name);
    ~zMaskProfile();

    void reset() { pressureleaks.clear(); }
    void scanLeaks(Session * session);
    void scanPressure(Session * session);
    void updatePressureMin();
    void updateProfile(Session * session);

    void load(Profile * profile);
    void save();

    QMap<EventStoreType, EventDataType> pressuremax;
    QMap<EventStoreType, EventDataType> pressuremin;
    QMap<EventStoreType, EventDataType> pressurecount;
    QMap<EventStoreType, EventDataType> pressuretotal;
    QMap<EventStoreType, EventDataType> pressuremean;
    QMap<EventStoreType, EventDataType> pressurestddev;

    QVector<TimeValue> Pressure;

    EventDataType calcLeak(EventStoreType pressure);

protected:
    static const quint32 version=0;
    void scanLeakList(EventList * leak);
    void scanPressureList(EventList * el);

    MaskType    m_type;
    QString     m_name;
    Profile     * m_profile;
    QString     m_filename;

    QMap<EventStoreType, QMap<EventStoreType,quint32> > pressureleaks;
    EventDataType maxP,minP,maxL,minL,m_factor;
};


bool operator<(const TimeValue &p1, const TimeValue & p2) {
    return p1.time < p2.time;
}

zMaskProfile::zMaskProfile(MaskType type, QString name)
    :m_type(type), m_name(name)
{
}
zMaskProfile::~zMaskProfile()
{
    save();
}

void zMaskProfile::load(Profile * profile)
{
    m_profile=profile;
    m_filename=m_profile->Get("{"+STR_GEN_DataFolder+"}/MaskProfile.mp");
    QFile f(m_filename);
    if (!f.open(QFile::ReadOnly))
        return;

    QDataStream in(&f);
    in.setVersion(QDataStream::Qt_4_6);
    in.setByteOrder(QDataStream::LittleEndian);

    quint32 m,v;
    in >> m;
    in >> v;
    if (m!=magic) {
        qDebug() << "Magic wrong in zMaskProfile::load";
    }

    in >> pressureleaks;
    f.close();
}
void zMaskProfile::save()
{
    QFile f(m_filename);
    if (!f.open(QFile::WriteOnly))
        return;

    QDataStream out(&f);
    out.setVersion(QDataStream::Qt_4_6);
    out.setByteOrder(QDataStream::LittleEndian);

    out << (quint32)magic;
    out << (quint32)version;

    out << pressureleaks;
    f.close();
}

void zMaskProfile::scanPressureList(EventList * el)
{
    qint64 start=el->first();
    int count=el->count();
    EventStoreType * dptr=el->rawData();
    EventStoreType * eptr=dptr+count;
    quint32 * tptr=el->rawTime();
    qint64 time;
    EventStoreType pressure;

    for (;dptr < eptr; dptr++) {
        time=start+ *tptr++;
        pressure=*dptr;
        Pressure.push_back(TimeValue(time,pressure));
    }
}
void zMaskProfile::scanPressure(Session * session)
{
    Pressure.clear();

    int prescnt=0;
    //EventStoreType pressure;
    if (session->eventlist.contains(CPAP_Pressure)) {
        prescnt=session->count(CPAP_Pressure);
        Pressure.reserve(prescnt);
        for (int j=0;j<session->eventlist[CPAP_Pressure].size();j++) {
            QVector<EventList *> & el=session->eventlist[CPAP_Pressure];
            for (int e=0;e<el.size();e++) {
                scanPressureList(el[e]);
            }
        }
    } else if (session->eventlist.contains(CPAP_IPAP)) {
        prescnt=session->count(CPAP_IPAP);
        Pressure.reserve(prescnt);
        for (int j=0;j<session->eventlist[CPAP_IPAP].size();j++) {
            QVector<EventList *> & el=session->eventlist[CPAP_IPAP];
            for (int e=0;e<el.size();e++) {
                scanPressureList(el[e]);
            }
        }
    }
    qSort(Pressure);
}
void zMaskProfile::scanLeakList(EventList * el)
{
    qint64 start=el->first();
    int count=el->count();
    EventStoreType * dptr=el->rawData();
    EventStoreType * eptr=dptr+count;
    quint32 * tptr=el->rawTime();
    //EventDataType gain=el->gain();

    EventStoreType pressure,leak;

    //EventDataType fleak;
    QMap<EventStoreType, EventDataType>::iterator pmin;
    qint64 ti;
    bool found;
    for (;dptr<eptr;dptr++) {
        leak=*dptr;
        ti=start + *tptr++;

        found=false;
        pressure=Pressure[0].value;
        if (Pressure.size()>1) {
            for (int i=0;i<Pressure.size()-1;i++) {
                const TimeValue & p1=Pressure[i];
                const TimeValue & p2=Pressure[i+1];

                if ((p2.time > ti) && (p1.time <= ti)) {
                    pressure=p1.value;
                    found=true;
                    break;
                } else if (p2.time==ti) {
                    pressure=p2.value;
                    found=true;
                    break;
                }
            }
        } else {
            found=true;
        }
        if (found) {
            pressureleaks[pressure][leak]++;
//            pmin=pressuremin.find(pressure);
//            fleak=EventDataType(leak) * gain;
//            if (pmin==pressuremin.end()) {
//                pressuremin[pressure]=fleak;
//                pressurecount[pressure]=pressureleaks[pressure][leak];
//            } else {
//                if (pmin.value() > fleak) {
//                    pmin.value()=fleak;
//                    pressurecount[pressure]=pressureleaks[pressure][leak];
//                }
//            }
        } else {
            //int i=5;
        }
    }


}
void zMaskProfile::scanLeaks(Session * session)
{
    QVector<EventList *> & elv=session->eventlist[CPAP_LeakTotal];
    for (int i=0;i<elv.size();i++) {
        scanLeakList(elv[i]);
    }
}
void zMaskProfile::updatePressureMin()
{
    QMap<EventStoreType, QMap<EventStoreType,quint32> >::iterator it;

    EventStoreType pressure;
    //EventDataType perc=0.1;
    //EventDataType tmp,l1,l2;
    //int idx1, idx2,
    int lks;
    double SN;
    double percentile=0.40;

    double p=100.0*percentile;
    double nth,nthi;

    int sum1,sum2,w1,w2,N,k;

    for (it=pressureleaks.begin();it!=pressureleaks.end();it++) {
        pressure=it.key();
        QMap<EventStoreType,quint32> & leakmap = it.value();
        lks=leakmap.size();
        SN=0;

        // First sum total counts of all leaks
        for (QMap<EventStoreType,quint32>::iterator lit=leakmap.begin();lit!=leakmap.end();lit++) {
            SN+=lit.value();
        }


        nth=double(SN)*percentile; // index of the position in the unweighted set would be
        nthi=floor(nth);

        sum1=0,sum2=0;
        w1=0,w2=0;

        EventDataType v1=0,v2;

        N=leakmap.size();
        k=0;

        bool found=false;
        double total=0;
        for (QMap<EventStoreType,quint32>::iterator lit=leakmap.begin();lit!=leakmap.end();lit++,k++) {
            total+=lit.value();
        }
        pressuretotal[pressure]=total;
        for (QMap<EventStoreType,quint32>::iterator lit=leakmap.begin();lit!=leakmap.end();lit++,k++) {
        //for (k=0;k < N;k++) {
            v1=lit.key();
            w1=lit.value();
            sum1+=w1;

            if (sum1 > nthi) {
                pressuremin[pressure]=v1;
                pressurecount[pressure]=w1;
                found=true;
                break;
            }
            if (sum1 == nthi) {
                break; // boundary condition
            }
        }
        if (found) continue;
        if (k>=N) {
            pressuremin[pressure]=v1;
            pressurecount[pressure]=w1;
            continue;
        }

        v2=(leakmap.begin()+(k+1)).key();
        w2=(leakmap.begin()+(k+1)).value();
        sum2=sum1+w2;
        // value lies between v1 and v2

        double px=100.0/double(SN); // Percentile represented by one full value

        // calculate percentile ranks
        double p1=px * (double(sum1)-(double(w1)/2.0));
        double p2=px * (double(sum2)-(double(w2)/2.0));

        // calculate linear interpolation
        double v=v1 + ((p-p1)/(p2-p1)) * (v2-v1);
        pressuremin[pressure]=v;
        pressurecount[pressure]=w1;

    }
}

EventDataType zMaskProfile::calcLeak(EventStoreType pressure)
{
    if (maxP==minP)
        return pressuremin[pressure];

    EventDataType leak=(pressure-minP) * (m_factor) + minL;
    return leak;

}

void zMaskProfile::updateProfile(Session * session)
{
    scanPressure(session);
    scanLeaks(session);
    updatePressureMin();


    if (pressuremin.size()<=1) {
        maxP=minP=0;
        maxL=minL=0;
        m_factor=0;
        return;
    }
    EventDataType p,l,tmp,mean,sum;
    minP=250,maxP=0;
    minL=1000, maxL=0;

    long cnt=0, vcnt;
    int n;

    EventDataType maxcnt, maxval, lastval, lastcnt;

    for (QMap<EventStoreType, QMap<EventStoreType,quint32> >::iterator it=pressureleaks.begin();it!=pressureleaks.end();it++)  {
        p=it.key();
        l=pressuremin[p];
        QMap<EventStoreType,quint32> & leakval=it.value();
        cnt=0;
        vcnt=-1;
        n=leakval.size();
        sum=0;

        maxcnt=0, maxval=0, lastval=0, lastcnt=0;
        for (QMap<EventStoreType,quint32>::iterator lit=leakval.begin();lit!=leakval.end();lit++) {
            cnt+=lit.value();
            if (lit.value() > maxcnt) {
                lastcnt=maxcnt;
                maxcnt=lit.value();
                lastval=maxval;
                maxval=lit.key();

            }
            sum+=lit.key() * lit.value();
            if (lit.key()==(EventStoreType)l) {
                vcnt=lit.value();
            }
        }
        pressuremean[p]=mean=sum / EventDataType(cnt);
        if (lastval > 0) {
            maxval=qMax(maxval,lastval); //((maxval*maxcnt)+(lastval*lastcnt)) / (maxcnt+lastcnt);
        }
        pressuremax[p]=lastval;
        sum=0;
        for (QMap<EventStoreType,quint32>::iterator lit=leakval.begin();lit!=leakval.end();lit++) {
            tmp=lit.key()-mean;
            sum+=tmp * tmp;
        }
        pressurestddev[p]=tmp=sqrt(sum / EventDataType(n));

        if (vcnt>=0) {
            tmp=1.0 / EventDataType(cnt) * EventDataType(vcnt);

        }
    }
    QMap<EventStoreType, EventDataType> pressureval;
    QMap<EventStoreType, EventDataType> pressureval2;
    EventDataType max=0,tmp2,tmp3;
    for (QMap<EventStoreType, EventDataType>::iterator it=pressuretotal.begin();it!=pressuretotal.end();it++) {
        if (max < it.value()) max=it.value();
    }
    for (QMap<EventStoreType, EventDataType>::iterator it=pressuretotal.begin();it!=pressuretotal.end();it++) {
        p=it.key();
        tmp=pressurecount[p];
        tmp2=it.value();

        tmp3=(tmp / tmp2) * (tmp2 / max);
        if (tmp3 > 0.05) {
            tmp=pressuremean[p];
            if (tmp>0) {
                pressureval[p]=tmp; // / tmp2;
                if (p < minP) minP = p;
                if (p > maxP) maxP = p;
                if (tmp < minL) minL = tmp;
                if (tmp > maxL) maxL = tmp;
            }
        }
    }

    if ((maxP==minP)  || (minL==maxL) || (minP==250) || (minL==1000)) {
        // crappy data set
        maxP=minP=0;
        maxL=minL=0;
        m_factor=0;
        return;
    }
    m_factor = (maxL - minL) / (maxP - minP);

//    for (QMap<EventStoreType, EventDataType>::iterator it=pressuremin.begin();it!=pressuremin.end()-1;it++) {
//        p1=it.key();
//        p2=(it+1).key();
//        l1=it.value();
//        l2=(it+1).value();
//        if (l2 > l1) {
//            factor=(l2 - l1) / (p2 - p1);
//            sum+=factor;
//            cnt++;
//        }
//    }

//    m_factor=sum/double(cnt);

//    int i=0;
//    if (i==1) {
//        QFile f("/home/mark/leaks.csv");
//        f.open(QFile::WriteOnly);
//        QTextStream out(&f);
//        EventDataType p,l,c;
//        QString fmt;
//        for (QMap<EventStoreType, EventDataType>::iterator it=pressuremin.begin();it!=pressuremin.end();it++) {
//            p=EventDataType(it.key()/10.0);
//            l=it.value();
//            fmt=QString("%1,%2\n").arg(p,0,'f',1).arg(l);
//            out << fmt;
//        }

    // cruft
//        for (QMap<EventStoreType, QMap<EventStoreType,quint32> >::iterator it=pressureleaks.begin();it!=pressureleaks.end();it++)  {
//            QMap<EventStoreType,quint32> & leakval=it.value();
//            for (QMap<EventStoreType,quint32>::iterator lit=leakval.begin();lit!=leakval.end();lit++) {
//                l=lit.key();
//                c=lit.value();
//                fmt=QString("%1,%2,%3\n").arg(p,0,'f',2).arg(l).arg(c);
//                out << fmt;
//            }
//        }
//        f.close();
//    }
}

QMutex zMaskmutex;
zMaskProfile mmaskProfile(Mask_NasalPillows,"ResMed Swift FX");
bool mmaskFirst=true;

int calcLeaks(Session *session)
{

    if (session->machine()->GetType()!=MT_CPAP) return 0;
    if (session->eventlist.contains(CPAP_Leak)) return 0; // abort if already there
    if (!session->eventlist.contains(CPAP_LeakTotal)) return 0; // can't calculate without this..

    zMaskmutex.lock();
    zMaskProfile * maskProfile=&mmaskProfile;
    if (mmaskFirst) {
        mmaskFirst=false;
        maskProfile->load(p_profile);
    }
//    if (!maskProfile) {
//        maskProfile=new zMaskProfile(Mask_NasalPillows,"ResMed Swift FX");
//    }
    maskProfile->reset();
    maskProfile->updateProfile(session);

    EventList *leak=session->AddEventList(CPAP_Leak,EVL_Event,1);

    for (int i=0;i<session->eventlist[CPAP_LeakTotal].size();i++) {
        EventList & el=*session->eventlist[CPAP_LeakTotal][i];
        EventDataType gain=el.gain(),tmp,val;
        int count = el.count();
        EventStoreType *dptr=el.rawData();
        EventStoreType *eptr=dptr+count;
        quint32 * tptr=el.rawTime();
        qint64 start=el.first(),ti;
        EventStoreType pressure;
        tptr=el.rawTime();
        start=el.first();

        bool found;
        for (; dptr < eptr; dptr++) {
            tmp=EventDataType(*dptr) * gain;
            ti=start+ *tptr++;

            found=false;
            pressure=maskProfile->Pressure[0].value;
            for (int i=0;i<maskProfile->Pressure.size()-1;i++) {
                const TimeValue & p1=maskProfile->Pressure[i];
                const TimeValue & p2=maskProfile->Pressure[i+1];

                if ((p2.time > ti) && (p1.time <= ti)) {
                    pressure=p1.value;
                    found=true;
                    break;
                } else if (p2.time==ti) {
                    pressure=p2.value;
                    found=true;
                    break;
                }
            }
            if (found) {
                val=tmp-maskProfile->calcLeak(pressure);

                if (val < 0) {
                    val=0;
                }

                leak->AddEvent(ti,val);
            }
        }
    }
    zMaskmutex.unlock();
    return leak->count();
}


int calcPulseChange(Session *session)
{
    if (session->eventlist.contains(OXI_PulseChange)) return 0;

    QHash<ChannelID,QVector<EventList *> >::iterator it=session->eventlist.find(OXI_Pulse);
    if (it==session->eventlist.end()) return 0;

    EventDataType val,val2,change,tmp;
    qint64 time,time2;
    qint64 window=PROFILE.oxi->pulseChangeDuration();
    window*=1000;

    change=PROFILE.oxi->pulseChangeBPM();

    EventList *pc=new EventList(EVL_Event,1,0,0,0,0,true);
    pc->setFirst(session->first(OXI_Pulse));
    qint64 lastt;
    EventDataType lv=0;
    int li=0;

    int max;
    for (int e=0;e<it.value().size();e++) {
        EventList & el=*(it.value()[e]);

        for (unsigned i=0;i<el.count();i++) {
            val=el.data(i);
            time=el.time(i);


            lastt=0;
            lv=change;
            max=0;

            for (unsigned j=i+1;j<el.count();j++) { // scan ahead in the window
                time2=el.time(j);
                if (time2 > time+window) break;
                val2=el.data(j);
                tmp=qAbs(val2-val);
                if (tmp > lv) {
                    lastt=time2;
                    if (tmp>max) max=tmp;
                    //lv=tmp;
                    li=j;
                }
            }
            if (lastt>0) {
                qint64 len=(lastt-time)/1000.0;
                pc->AddEvent(lastt,len,tmp);
                i=li;

            }
        }
    }
    if (pc->count()==0) {
        delete pc;
        return 0;
    }
    session->eventlist[OXI_PulseChange].push_back(pc);
    session->setMin(OXI_PulseChange,pc->Min());
    session->setMax(OXI_PulseChange,pc->Max());
    session->setCount(OXI_PulseChange,pc->count());
    session->setFirst(OXI_PulseChange,pc->first());
    session->setLast(OXI_PulseChange,pc->last());
    return pc->count();
}


int calcSPO2Drop(Session *session)
{
    if (session->eventlist.contains(OXI_SPO2Drop)) return 0;

    QHash<ChannelID,QVector<EventList *> >::iterator it=session->eventlist.find(OXI_SPO2);
    if (it==session->eventlist.end()) return 0;

    EventDataType val,val2,change,tmp;
    qint64 time,time2;
    qint64 window=PROFILE.oxi->spO2DropDuration();
    window*=1000;
    change=PROFILE.oxi->spO2DropPercentage();

    EventList *pc=new EventList(EVL_Event,1,0,0,0,0,true);
    qint64 lastt;
    EventDataType lv=0;
    int li=0;

    // Fix me.. Time scale varies.
    //const unsigned ringsize=30;
    //EventDataType ring[ringsize]={0};
    //qint64 rtime[ringsize]={0};
    //int rp=0;
    int min;
    int cnt=0;
    tmp=0;

    qint64 start=0;

    // Calculate median baseline
    QList<EventDataType> med;
    for (int e=0;e<it.value().size();e++) {
        EventList & el=*(it.value()[e]);
        for (unsigned i=0;i<el.count();i++) {
            val=el.data(i);
            time=el.time(i);
            if (val>0) med.push_back(val);
            if (!start) start=time;
            if (time > start+3600000) break; // just look at the first hour
            tmp+=val;
            cnt++;
        }
    }
    EventDataType baseline=0;
    if (med.size()>0) {
        qSort(med);

        int midx=float(med.size())*0.90;
        if (midx>med.size()-1) midx=med.size()-1;
        if (midx<0) midx=0;
        baseline=med[midx];
    }
    session->settings[OXI_SPO2Drop]=baseline;
    //EventDataType baseline=round(tmp/EventDataType(cnt));
    EventDataType current;
    qDebug() << "Calculated baseline" << baseline;

    for (int e=0;e<it.value().size();e++) {
        EventList & el=*(it.value()[e]);
        for (unsigned i=0;i<el.count();i++) {
            current=el.data(i);
            if (!current) continue;
            time=el.time(i);
            /*ring[rp]=val;
            rtime[rp]=time;
            rp++;
            rp=rp % ringsize;
            if (i<ringsize)  {
                for (unsigned j=i;j<ringsize;j++) {
                    ring[j]=val;
                    rtime[j]=0;
                }
            }
            tmp=0;
            cnt=0;
            for (unsigned j=0;j<ringsize;j++) {
                if (rtime[j] > time-300000) {  // only look at recent entries..
                    tmp+=ring[j];
                    cnt++;
                }
            }
            if (!cnt) {
                unsigned j=abs((rp-1) % ringsize);
                tmp=(ring[j]+val)/2;
            } else tmp/=EventDataType(cnt); */

            val=baseline;
            lastt=0;
            lv=val;


            min=val;
            for (unsigned j=i;j<el.count();j++) { // scan ahead in the window
                time2=el.time(j);
                //if (time2 > time+window) break;
                val2=el.data(j);

                if (val2 > baseline-change) break;
                lastt=time2;
                li=j+1;
            }
            if (lastt>0) {
                qint64 len=(lastt-time);
                if (len>=window) {
                    pc->AddEvent(lastt,len/1000,val-min);
                    i=li;
                }
            }
        }
    }
    if (pc->count()==0) {
        delete pc;
        return 0;
    }
    session->eventlist[OXI_SPO2Drop].push_back(pc);
    session->setMin(OXI_SPO2Drop,pc->Min());
    session->setMax(OXI_SPO2Drop,pc->Max());
    session->setCount(OXI_SPO2Drop,pc->count());
    session->setFirst(OXI_SPO2Drop,pc->first());
    session->setLast(OXI_SPO2Drop,pc->last());
    return pc->count();
}