decodingSupervised.py

#  Copyright (c) 2019 MindAffect B.V. 
#  Author: Jason Farquhar <jason@mindaffect.nl>
# This file is part of pymindaffectBCI <https://github.com/mindaffect/pymindaffectBCI>.
#
# pymindaffectBCI 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 3 of the License, or
# (at your option) any later version.
#
# pymindaffectBCI 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 pymindaffectBCI.  If not, see <http://www.gnu.org/licenses/>

import numpy as np
from mindaffectBCI.decoder.normalizeOutputScores import normalizeOutputScores
from mindaffectBCI.decoder.zscore2Ptgt_softmax import zscore2Ptgt_softmax
#@function
def decodingSupervised(Fy, softmaxscale=3.5, marginalizemodels=True, 
                       marginalizedecis=False,
                       prior=None,
                       nocontrolamplitude=None,
                       tiebreaking_noise=1e-3, **kwargs):
  """    true-target estimator and error-probility estimator for each trial

   Args:
      Fy (nModel,nTrl,nEp,nY): The output scores.
      softmaxscale (float, optional): the scale length to pass to zscore2Ptgt_softmax.py, this is the inverse *noise* sigma. Defaults to 3.5.
      badFyThresh (int, optional): threshold for detction of bad Fy entry, in std-dev. Defaults to 4.
      centFy (bool, optional): bool, do we center Fy before computing  *important*  (true). Defaults to True.
      detrendFy (bool, optional): [description]. Defaults to True.
      nEpochCorrection (int, optional): int, number of epochs to use a base for correction of number epochs in the sum.
               such that: zscore = Fy / ( sigma * ( 1+1/(nEpoch/nEpochErrorCorrection)) ),
               basically, nEpoch < nEpochCorrection have highly increased Perr, so unlikely to be selected. Defaults to 100.
      minDecisLen (int, optional): int, number of epochs to use as base for distribution of time-based decision points.
                         i.e. decisions at, [1,2,4,...2^n]*exptDistDecis
                    OR: minDecisLen<0 => decision point every abs(minDeicsLen) epochs. Defaults to 0.
      maxDecisLen (int, optional): maximum number of epochs for a decision. Defaults to 0.
      bwdAccumulate (bool, optional): accumulate data backwards from last epoch gathered. Defaults to False.
      marginalizemodels (bool, optional): [description]. Defaults to True.
      marginalizedecis (bool, optional): [description]. Defaults to False.
      prior ([type], optional): [description]. Defaults to None.
      nocontrolamplitude ([type], optional): [description]. Defaults to None.
      priorsigma (tuple, optional): (sigma,N) prior estimate of sigma2 and number pseudo-points. Defaults to (-1,120).
      tiebreaking_noise ([type], optional): [description]. Defaults to 1e-3.

   Raises:
      NotImplementedError: [description]

   Returns:
      Yest (nTrl,nDecis): the most likely / minimum error output for each decision point
      Perr (nTrl,nDecis): the probability that this selection is an ERROR for each decision point
      Ptgt (nTrl,nDecis,nY): the probability each target is the true target for each decision point
      decisMdl, 
      decisEp

    Copyright (c) MindAffect B.V. 2018
  """   
  if Fy is None:
      return -1, 1, None, None, None
  
  #print("decodingSup args={}".format(kwargs))

  # get the info on which outputs are zero in each trial
  validTgt = np.any(Fy != 0, axis=-2) # valid if non-zero for any epoch..# (nModel,nTrl,nY)  
  # normalize the raw scores for each model to have nice distribution
  ssFy,varsFy,decisIdx,nEp,nY=normalizeOutputScores(Fy, validTgt=validTgt, **kwargs)

  if nocontrolamplitude is not None:
    raise NotImplementedError('no-control signal not yet implemented correctly')
    # # add a no-control pseudo-output to simulate not looking
    # #Median rather than mean?
    # mussFy=np.sum(ssFy,0)/np.maximum(2,np.sum(validTgt,0,keepdims=True)) # mean score for each trial/decisPt [ 1 x nDecis x nTrl x nMdl ]
    # # add to the scores & validTgt info
    # ssFy=np.concatenate((ssFy,nocontrolamplitude+mussFy),0)
    # validtgtTrl=np.concatenate((validtgtTrl,np.ones([1,size(validtgtTrl,2),size(validtgtTrl,3),size(validtgtTrl,4)])))

  # compute the target probabilities over output for each model+trial
  # use the softmax approach to get Ptgt for all outputs
  Ptgt = zscore2Ptgt_softmax(ssFy,
                              softmaxscale,
                              validTgt=validTgt,
                              marginalizemodels=marginalizemodels,
                              marginalizedecis=marginalizedecis,
                              prior=prior) # (nM,nTrl,nDecis,nY)
  # extract the predicted output and it's probability of being the target
  Ptgt2d = Ptgt.reshape((np.prod(Ptgt.shape[:-1]), Ptgt.shape[-1])) # make 2d-copy
  # add tie-breaking noise
  if tiebreaking_noise > 0:
      Ptgt2d = Ptgt2d + np.random.standard_normal(Ptgt2d.shape)*tiebreaking_noise

  Yestidx = np.argmax(Ptgt2d, -1) # max over outputs, i.e. models, decisPts, etc..
  Ptgt_max = Ptgt2d[np.arange(Ptgt2d.shape[0]), Yestidx] # value at max, indexing trick to find..
  Yestidx = Yestidx.reshape(Ptgt.shape[:-1]) # -> (nM,nTrl,nY)
  Ptgt_max = Ptgt_max.reshape(Ptgt.shape[:-1])# -> (nM,nTrl,nY)

  # pick a model for this maxY
  decisMdl = Yestidx if ssFy.ndim>3 else 1
  decisEp = 1   
  Yest = Yestidx
  # p(maxz != tgt ) = 1 - p(maxz==tgt)
  Perr = 1 - Ptgt_max
  return Yest, Perr, Ptgt, decisMdl, decisEp


def decodingSupervised_streamed(Fy,softmaxscale=3,nocontrolamplitude=None,
                                marginalizemodels=False):
    '''
    true-target estimator and error-probility estimator for each trial
    Inputs:
      Fy   - (nModel,nTrl,nEp,nY) [#Y x #Epoch x #Trials x nModels]
      softmaxscale - the scale length to pass to zscore2Ptgt_softmax.py
      badFyThresh - threshold for detction of bad Fy entry, in std-dev
      centFy   - bool, do we center Fy before computing  *important*  (true)
      nEpochCorrection - int, number of epochs to use a base for correction of number epochs in the sum.
               such that: zscore = Fy / ( sigma * ( 1+1/(nEpoch/nEpochErrorCorrection)) ),
               basically, nEpoch < nEpochCorrection have highly increased Perr, so unlikely to be selected
      minDecisLen - int, number of epochs to use as base for distribution of time-based decision points.
                         i.e. decisions at, [1,2,4,...2^n]*exptDistDecis
                    OR: minDecisLen<0 => decision point every abs(minDeicsLen) epochs
      bwdAccumulate - [bool], accumulate data backwards from last epoch gathered
      maxDecisLen   - maximum number of epochs for a decision
    Outputs:
      Yest - (nTrl,nDecis) [ nDecis x nTrl ] the most likely / minimum error output for each decision point
      Perr - (nTrl,nDecis) [ nDecix x nTrl ] the probability that this selection is an ERROR for each decision point
      Ptgt - (nTrl,nDecis,nY) [ nY x nDecis x nTrl ] the probability each target is the true target for each decision point

    Copyright (c) MindAffect B.V. 2018
    '''
    if Fy is None:
        return -1, 1, None, None, None
    
    # get the info on which outputs are zero in each trial
    validTgt = np.any(Fy != 0, axis=-2) # valid if non-zero for any epoch..# (nModel,nTrl,nY) [ nY x nTrl x nModel ]   
    # normalize the raw scores for each model to have nice distribution
    stdFy,varsFy,N,nEp,nY=normalizeOutputScores_streamed(Fy, validTgt=validTgt,
                                                         badFyThresh=badFyThresh, centFy=centFy,
                                                         minDecisLen=minDecisLen, filtLen=filtLen)



    # compute the target probabilities over output for each block+model+trial
    # use the softmax approach to get Ptgt for all outputs
    Ptgt = zscore2Ptgt_softmax(ssFy,
                               softmaxscale,
                               validTgt=validTgt,
                               marginalizemodels=False) # (nM,nTrl,nDecis,nY) [ nY x nDecis x nTrl]
    # extract the predicted output and it's probability of being the target
    Ptgt2d = Ptgt.reshape((np.prod(Ptgt.shape[:-1]), Ptgt.shape[-1])) # make 2d-copy
    Yestidx = np.argmax(Ptgt2d, -1) # max over outputs
    Ptgt_max = Ptgt2d[np.arange(Ptgt2d.shape[0]), Yestidx] # value at max, indexing trick to find..
    Yestidx = Yestidx.reshape(Ptgt.shape[:-1]) # -> (nM,nTrl,nY)
    Ptgt_max = Ptgt_max.reshape(Ptgt.shape[:-1])# -> (nM,nTrl,nY)

    # pick a model for this maxY
    decisMdl = Yestidx if ssFy.ndim>3 else 1
    decisEp = 1   
    Yest = Yestidx
    # p(maxz != tgt ) = 1 - p(maxz==tgt)
    Perr = 1 - Ptgt_max
    return Yest, Perr, Ptgt, decisMdl, decisEp



#@function
def testcase():
    from normalizeOutputScores import mktestFy
    Fy,nEp=mktestFy(startupNoisefrac=0,trlenfrac=0,sigstr=.4)
    from decodingSupervised import decodingSupervised
    print("Fy={}".format(Fy.shape))
    Yest,Perr,Ptgt,decismdl,decisEp=decodingSupervised(Fy)
    print("Fy={}".format(Fy[:,[0],:,:].shape))
    Yest,maxp,ssFy,decismdl,decisEp=decodingSupervised(Fy[:,[0],:,:])
    np.mean(Yest == 0,1) 
    from decodingCurveSupervised import decodingCurveSupervised
    decodingCurveSupervised(Fy);

if __name__=="__main__":
    testcase()