from numpy import *
import inspect

class Quantity(object):
    """docstring for Quantity"""
    def __init__(self, value,std=None):
        
        if isinstance(value,str):
            
            parts=value.split('+-')
            if len(parts)==1:  # look for 456(23) notation
            
                parts=value.split('(')
                self.value=float(parts[0])
                
                parts[1]=parts[1].replace(')','').strip()
                
                parts[0]=parts[0].strip()
                for c in '123456789':
                    parts[0]=parts[0].replace(c,'0')
                    
                parts[0]=parts[0][:-len(parts[1])]+parts[1]
                
                self.std=float(parts[0])
                
            
            elif len(parts)==2:
                self.value=float(parts[0])
            
                if '%' in parts[1]:
                    percent=float(parts[1].replace('%',''))
                    self.std=abs(self.value*percent/100.0)
                else:            
                    self.std=float(parts[1])
        
        else:
            
            if std is None:
                std=0
                
            self.value=float(value)
            self.std=float(std)
            
        super(Quantity, self).__init__()
            
    def sample(self,N):
        
        r=random.randn(N)
        
        return (r*self.std)+self.value
            
    def __repr__(self):
        
        s=str(self.value)+" +- "+str(self.std)
        return s

def evaluate(equation_string,sample_size=10000):

    equation_string=''.join(equation_string.split())
    
    f=inspect.currentframe()
    out=inspect.getouterframes(f)
    parent_dict=out[1][0].f_locals
    vardict=parent_dict
    
    eval_dict={}
    for name in vardict:
        if isinstance(vardict[name],Quantity):
            eval_dict[name]=vardict[name].sample(sample_size)    
        else:
            eval_dict[name]=vardict[name]
    
    
    _val_=eval(equation_string,eval_dict)
    q=Quantity(_val_.mean(),_val_.std())
    
    return q
    
    
if __name__=="__main__":
    
    # examples taken from "Measurements and Uncertanties"
    # John Denker, http://www.av8n.com/physics/uncertainty.htm
    
    print "Example 1"
    # I do specify the values in several ways, to show the possibilities
    a=Quantity(2) # no uncertainty
    b=Quantity(0.4,0.04)  # 10 percent
    c=Quantity('0.4 +- 0.04')  # 10 percent
    d=Quantity('0.4 +- 10%')  # 10 percent
    e=b

    x=evaluate('(((a+b)+c)+d)+e')

    print "x = (((2 + 0.4) + 0.4) + 0.4) + 0.4 : ",x


    print "Example 2"

    Mg24_mass=Quantity("23.9850423(8)")
    Mg25_mass=Quantity("24.9858374(8)")
    Mg26_mass=Quantity("25.9825937(8)")

    Mg24_abundance=Quantity("78.99(4)")
    Mg25_abundance=Quantity("10.01(1)")
    Mg26_abundance=Quantity("11.01(3)")

    x=evaluate("""Mg26_mass*Mg26_abundance/100+
                  Mg25_mass*Mg25_abundance/100+
                  Mg24_mass*(100-Mg26_abundance-Mg25_abundance)/100""")

    print "Mg mass with sum rule: ",x


    x=evaluate("""Mg26_mass*Mg26_abundance/100+
                  Mg25_mass*Mg25_abundance/100+
                  Mg24_mass*Mg24_abundance/100""")

    print "Mg mass without sum rule: ",x