import math


def calcFrag(values):
    """Given a set of parsed values that were generated by macs2 predictd,
    see get_size fn, this fn calculates the estimated fragment size and 
    the sd.
    
    **IMPORTANT: this fn is a python translation of the R code in 
    chilin2/modules/macs2_fragment/qc.py -- stat_frag_std

    RETURNS: (estimated frag, sd)
    """
    #calculate the estimated frag size: xmax
    ymax = max(values['ycorr'])
    i_ymax = values['ycorr'].index(ymax)
    xmax = values['xcorr'][i_ymax]
    #print(ymax, xmax)

    #find expected
    p_expect=sum([x* p/100.0 for (x,p) in zip(values['x'],values['positive'])])
    m_expect=sum([x* m/100.0 for (x,m) in zip(values['x'],values['minus'])])
    #print(p_expect, m_expect)

    #calc sd
    p_sd = math.sqrt(sum([((x - p_expect)**2)* p/100.0 \
                              for (x,p) in zip(values['x'],values['positive'])]))
    m_sd = math.sqrt(sum([((x - m_expect)**2)* m/100.0 \
                              for (x,m) in zip(values['x'],values['minus'])]))
    #print(p_sd, m_sd)

    #FINAL avg std error
    avg_sd = (p_sd + m_sd) /2.0

    return (xmax, avg_sd)