""" DoubleClassify.py command_options list_of_data_records

For each of the records, first estimate a classification (apnea,
borderline, or normal) for the entire record, then classify (normal or
apnea) each minute in the record.

EG: python DoubleClassify.py --Amodel=data/mod_A --BCmodel=data/mod_C
           --Lmodel=data/mod_L 2.2 --Mmodel=data/mod_M 1.15
           --Hmodel=data/mod_H 1.05 --Annotations=data/summary_of_training
           --Dpath=data/Apnea --Results=data/reportTrain data/scoreTrain
           x01 x02 x03 x04 x05 x06 ...  x34 x35

Copyright (c) 2005, 2008 Andrew Fraser
This file is part of HMM_DS_Code.

HMM_DS_Code is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.

HMM_DS_Code is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License along with
this program; if not, write to the Free Software Foundation, Inc.,
59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

"""
import numpy, optparse, ApOb, pickle, sys
OP = optparse.OptionParser(
usage=
"""Usage: %prog [options] records""",
description=
"""Train an HMM with an observation model from ApOb.py.  Use the training data specified by "records".""")
OP.add_option('--Amodel',            help='filename')
OP.add_option('--BCmodel',           help='filename')
OP.add_option('--Lmodel', nargs=3,   help='filename Fudge Pow')
OP.add_option('--Mmodel', nargs=3,   help='filename Fudge Pow')
OP.add_option('--Hmodel', nargs=3,   help='filename Fudge Pow')
OP.add_option('--Annotations',
              help='File of expert annotations for all records')
OP.add_option('--Results',  nargs=2, help='filenames for results')
OP.add_option('--Dpath',             help='Path to .lphr and .resp files')
OP.add_option('--Single', action='store_true',
              help='Only do the first pass classification')
opt,records = OP.parse_args()

SamPerMin = 10

Amod  = pickle.load(open(opt.Amodel, 'r'))
BCmod = pickle.load(open(opt.BCmodel, 'r'))
Lmod  = ApOb.read_mod(*opt.Lmodel) # *opt.Lmodel is File_name, Fudge, Power
Mmod  = ApOb.read_mod(*opt.Mmodel)
Hmod  = ApOb.read_mod(*opt.Hmodel)
for mod in [Amod,BCmod,Lmod,Mmod,Hmod]:
    mod.Py_calc = mod.Py_wo_class

routines = [ApOb.read_lphr, ApOb.read_resp]
paths = [opt.Dpath+'/%s.lphr', opt.Dpath+'/%s.resp']
args = [4, None] # The AR order of the HR models is 4
#Ys,Yall = ApOb.read_records(routines,paths,args,records)
#print 'len(Ys)=%d'%len(Ys), 'Yall[0].shape=',Yall[0].shape
if opt.Results:
    old_stdout = sys.stdout
    sys.stdout = open(opt.Results[0],'w')
for record in records:
    HR = ApOb.read_lphr(opt.Dpath+'/%s.lphr',record,4)
    # The AR order of the HR models is 4
    Resp = ApOb.read_resp(opt.Dpath+'/%s.resp',record)
    T = min(len(HR),len(Resp))
    Y = [HR[:T],Resp[:T]]
    # Calculate 70% cumulative peak height and average variation
    lp = ApOb.read_data(opt.Dpath+'/%s.lphr'%record)[2,:T] # LP Heart Rate
    peaks = []
    L1 = 0
    W = 5
    for t in xrange(T):
        L1 += abs(lp[t])
        # Look at peaks
        if t < W or t+W >= T:
            continue
        if lp[t] < 0 or max(lp[t-W:t+W]) > lp[t]:
            continue
        peaks.append(lp[t])
    peaks.sort()
    L1 = L1/T
    R = peaks[int(.74*len(peaks))]/L1
    # Calculate the log likelihood ratio
    Amod.Py_calc(Y)
    A = Amod.forward()
    BCmod.Py_calc(Y)
    BC = BCmod.forward()
    llr = (A - BC)/T

    stat = R + .5*llr          # Was 0.5
    if stat < 2.39:            # Was 2.39
        Name = opt.Lmodel[0]
        model = Lmod
    elif stat > 2.55:          # Was 2.55
        Name = opt.Hmodel[0]
        model = Hmod
    else:
        Name = opt.Mmodel[0]
        model = Mmod
    print '%5s # use %-10s stat= %6.3f llr= %6.3f R= %6.3f L1=% 6.3f'%(
        record,Name,stat,llr,R,L1)
    if opt.Single:  # If you want only the first pass classification
        continue
    # Now classify the minutes
    assert model.__class__ is ApOb.Both_HMM
    N,AR1 = model.HR.A.shape
    assert N is model.N
    HR = ApOb.read_lphr(opt.Dpath+'/%s.lphr',record,AR1-1)
    Y = [HR[:T],Resp[:T]]
    Cseq = 2*model.class_decode(Y)-1 # +/- 1
    tot = 0
    for n in xrange(len(Cseq)):
        tot += Cseq[n]            # Add up Cseq over each minute
        if n%(SamPerMin*60) == 0: #If on a hour boundary
            if not n == 0:
                print string
            print n/(SamPerMin*60), " ", #Print number of hours done
            string = ''
        if n%SamPerMin == SamPerMin-1:   #If on a minute boundary
            if tot > 0:
                string += 'A'
            else:
                string += 'N'
            tot = 0
    print string
print >> sys.stderr, "pass one finished"
if not opt.Single:
    if opt.Results:
        sys.stdout = old_stdout
    print ApOb.score(opt.Annotations,opt.Results[0],records,verbose=True)


#Local Variables:
#mode:python
#End: