docs: def stencil(**kwargs):
      """
      Applying genotype calls to multi-way alignment incidence matrix
  
      :param alnfile: alignment incidence file (h5),
      :param gtypefile: genotype calls by GBRS (tsv),
      :param grpfile: gene ID to isoform ID mapping info (tsv)
      :return: genotyped version of alignment incidence file (h5)
      """
      alnfile = kwargs.get('alnfile')
      gtypefile = kwargs.get('gtypefile')
      grpfile = kwargs.get('grpfile')
      if grpfile is None:
          grpfile2chk = os.path.join(DATA_DIR, 'ref.gene2transcripts.tsv')
          if os.path.exists(grpfile2chk):
              grpfile = grpfile2chk
          else:
              print >> sys.stderr, '[gbrs::stencil] A group file is *not* given. Genotype will be stenciled as is.'
  
      # Load alignment incidence matrix ('alnfile' is assumed to be in multiway transcriptome)
      alnmat = emase.AlignmentPropertyMatrix(h5file=alnfile, grpfile=grpfile)
  
      # Load genotype calls
      hid = dict(zip(alnmat.hname, np.arange(alnmat.num_haplotypes)))
      gid = dict(zip(alnmat.gname, np.arange(len(alnmat.gname))))
      gtmask = np.zeros((alnmat.num_haplotypes, alnmat.num_loci))
      gtcall_g = dict.fromkeys(alnmat.gname)
      with open(gtypefile) as fh:
          if grpfile is not None:
              gtcall_t = dict.fromkeys(alnmat.lname)
              for curline in dropwhile(is_comment, fh):
                  item = curline.rstrip().split("\t")
                  g, gt = item[:2]
                  gtcall_g[g] = gt
                  hid2set = np.array([hid[c] for c in gt])
                  tid2set = np.array(alnmat.groups[gid[g]])
                  gtmask[np.meshgrid(hid2set, tid2set)] = 1.0
                  for t in tid2set:
                      gtcall_t[alnmat.lname[t]] = gt
          else:
              for curline in dropwhile(is_comment, fh):
                  item = curline.rstrip().split("\t")
                  g, gt = item[:2]
                  gtcall_g[g] = gt
                  hid2set = np.array([hid[c] for c in gt])
                  gtmask[np.meshgrid(hid2set, gid[g])] = 1.0
  
      alnmat.multiply(gtmask, axis=2)
      for h in xrange(alnmat.num_haplotypes):
          alnmat.data[h].eliminate_zeros()
      outfile = kwargs.get('outfile')
      if outfile is None:
          outfile = 'gbrs.stenciled.' + os.path.basename(alnfile)
      alnmat.save(h5file=outfile)