""" DoubleClassify.py command_options list_of_data_files

EG: python DoubleClassify.py --Model=mod --Annotations=reference
--Dpath=data/Apnea --Power 0.5 2.65 0.1 --Fudge .8 1.61 .05 a06 b01 ...

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
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('--Model', help='filename')
OP.add_option('--Dpath', help='directory of *lphr and *resp files')
OP.add_option('--Annotations',
    help='File of expert annotations for all records')
OP.add_option("--Power", type="float", nargs=3,
    help='"From To Step" for range of power settings')
OP.add_option("--Fudge", type="float", nargs=3,
    help='"From To Step" for range of fudge settings')
OP.add_option('--Results', help='filename for results')
opt,records = OP.parse_args()

Y_dict = {}
Ann_dict = {}
out = open(opt.Results,'w')
Mod  = ApOb.read_mod(opt.Model,1.0,1.0) # Get N and AR1
assert Mod.__class__ is ApOb.Both_HMM
N,AR1 = Mod.HR.A.shape
assert N is Mod.N
for record in records:
    HR = ApOb.read_lphr(opt.Dpath+'/%s.lphr',record,AR1-1)
    Resp = ApOb.read_resp(opt.Dpath+'/%s.resp',record)
    T = min(len(HR),len(Resp))
    Y_dict[record] = [HR[:T],Resp[:T]]
    Ann_dict[record] = numpy.array(ApOb.fetch_ann(opt.Annotations,
                                                  record),numpy.bool)
SPM = ApOb.SamPerMin
count = numpy.ones(SPM,numpy.int32)
for fudge in numpy.arange(*opt.Fudge):
    Mod  = ApOb.read_mod(opt.Model,fudge,1.0)
    for power in numpy.arange(*opt.Power):
        Mod.Pow = power
        T_tot = 0
        wrong_tot = 0
        for record in records:
            Y = Y_dict[record]
            Ann = Ann_dict[record]
            Cseq = Mod.class_decode(Y) # ints: 0 for normal 1 for apnea
            T = len(Ann)
            Mseq = numpy.zeros(T,numpy.bool)
            for t in xrange(T):
                if (t+1)*SPM <= len(Cseq):
                    Mseq[t] = (numpy.dot(count,Cseq[t*SPM:(t+1)*SPM])
                           > SPM/2)
                else:
                    Mseq[t] = Mseq[t-1]
            try:
                wrong_tot += (Ann ^ Mseq).sum()
            except:
                print 'Ann.shape=',Ann.shape,'Mseq.shape=',Mseq.shape
                print 'len(Cseq)=',len(Cseq)
                print '(Ann ^ Mseq).sum()',(Ann ^ Mseq).sum()
                wrong_tot += (Ann ^ Mseq).sum()
            T_tot += T
        frac_right = 1.0 - float(wrong_tot)/float(T_tot)
        print >>out, '%5.3f  %5.3f  %6.4f'%(power,fudge,frac_right)

#Local Variables:
#mode:python
#End: