In [2]:

    

%matplotlib inline
from matplotlib import pylab as pl
from scipy.optimize import curve_fit
import numpy as np

data= [i.strip().split() for i in open('a2_4-C_2-OS-CT-CT.hist')][2:]
x, y= [float(i[0]) for i in data], [float(i[1]) for i in data]



def gauss(x,mu,sigma,A):
    return A*np.exp(-(x-mu)**2/2/sigma**2)

def hund_modal(x,*params):
    result = []
    for i in range(-180,180,10):
        print params
        print params[i]
        result += gauss(x, i, params[i][0], params[i][1])
    return result
expected = [(0.2,2) for i in range(36)]
print expected
params,cov=curve_fit(hund_modal,x,y,expected)
sigma=np.sqrt(np.diag(cov))
print params
pl.plot(x,y)
pl.plot(x,hund_modal(x,*params),color='red',lw=3,label='model')
pl.legend()
print(params,'\n',sigma)


    

    




[(0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2), (0.2, 2)]
(0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0, 0.20000000000000001, 2.0)






    




---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-2-04aa0b8023cd> in <module>()
     21 expected = [(0.2,2) for i in range(36)]
     22 print expected
---> 23 params,cov=curve_fit(hund_modal,x,y,expected)
     24 sigma=np.sqrt(np.diag(cov))
     25 print params

/usr/local/lib/python2.7/dist-packages/scipy/optimize/minpack.pyc in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
    674         # Remove full_output from kwargs, otherwise we're passing it in twice.
    675         return_full = kwargs.pop('full_output', False)
--> 676         res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
    677         popt, pcov, infodict, errmsg, ier = res
    678         cost = np.sum(infodict['fvec'] ** 2)

/usr/local/lib/python2.7/dist-packages/scipy/optimize/minpack.pyc in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    375     if not isinstance(args, tuple):
    376         args = (args,)
--> 377     shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
    378     m = shape[0]
    379     if n > m:

/usr/local/lib/python2.7/dist-packages/scipy/optimize/minpack.pyc in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
     24 def _check_func(checker, argname, thefunc, x0, args, numinputs,
     25                 output_shape=None):
---> 26     res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
     27     if (output_shape is not None) and (shape(res) != output_shape):
     28         if (output_shape[0] != 1):

/usr/local/lib/python2.7/dist-packages/scipy/optimize/minpack.pyc in func_wrapped(params)
    453     if weights is None:
    454         def func_wrapped(params):
--> 455             return func(xdata, *params) - ydata
    456     else:
    457         def func_wrapped(params):

<ipython-input-2-04aa0b8023cd> in hund_modal(x, *params)
     16     for i in range(-180,180,10):
     17         print params
---> 18         print params[i]
     19         result += gauss(x, i, params[i][0], params[i][1])
     20     return result

IndexError: tuple index out of range

In [3]:

    

%%bash
head a2_4-C_2-OS-CT-CT.hist


    

    




# histogram values for [ a2_4-C_2-OS-CT-CT ]

-180.0	45372.0
-170.0	115639.0
-160.0	127833.0
-150.0	90430.0
-140.0	51008.0
-130.0	27227.0
-120.0	24347.0
-110.0	41734.0

In [38]:

    

from scipy.interpolate import InterpolatedUnivariateSpline
import numpy as np 
import matplotlib.pyplot as plt

from scipy.interpolate import interp1d

def normalized(a, axis=-1, order=2):
    l2 = np.atleast_1d(np.linalg.norm(a, order, axis))
    l2[l2==0] = 1
    return a / np.expand_dims(l2, axis)

pi = np.pi
data= [i.strip().split() for i in open('a2_4-C_2-OS-CT-CT.hist')][2:]
x, y= [float(i[0]) for i in data], [float(i[1]) for i in data]
x2 = x + [i+360 for i in x]
y2 = y + y
x = x2
y = y2


f2 = InterpolatedUnivariateSpline(x, y)
#Get dervative
der = []
for i in range(len(y)):
    h = 1e-4
    der.append( ( f2(x[i]+h)-f2(x[i]-h) )/(2*h) )
der = np.array(normalized(der)[0])
print y
yn = normalized(y)[0]
print x
xn = list(np.linspace(0,360,100))

fyn = interp1d(x, yn, kind='cubic') 
ynn = fyn(xn)
print yn
print len(der)
print len(x), len(xn)
plt.plot(xn, ynn , 'g', x, [0 for i in x],'b')#-np.sin(x)   x, der , 'r',
mu = []
for i in range(len(der)-1):
    if (der[i] > 0 and der[i+1] < 0):
        print der[i], der[i+1]
        if der[i] + np.abs(der[i+1]) < 0.025:
            continue
        mu.append(np.average([x[i],x[i+1]]))
        pl.vlines(np.average([x[i],x[i+1]]),0,0.1)
        
print mu
mu = [i for i in mu if i <=360 and i >= 0]
print mu
def gauss(x,sigma,A):
    mu = i
    return A*np.exp(-(x-mu)**2/2/sigma**2)

for i in mu:
    print i
    
    location = [k for k in zip(list(xn), list(ynn)) if k[0] > i-30 and k[0] < i+30]
    xn2 = np.array([k[0] for k in location])
    ynn2 = np.array([k[1] for k in location])
    params,cov=curve_fit(gauss,xn2,ynn2)
    pl.plot(xn,gauss(xn,*params),'r',lw=3,label='model')
plt.show()


    

    




[45372.0, 115639.0, 127833.0, 90430.0, 51008.0, 27227.0, 24347.0, 41734.0, 78048.0, 128000.0, 134983.0, 55630.0, 9907.0, 1034.0, 88.0, 9.0, 1.0, 0.0, 10003.0, 7.0, 38.0, 448.0, 3775.0, 15886.0, 24921.0, 12281.0, 2351.0, 251.0, 34.0, 11.0, 1.0, 0.0, 2.0, 49.0, 748.0, 8005.0, 45372.0, 115639.0, 127833.0, 90430.0, 51008.0, 27227.0, 24347.0, 41734.0, 78048.0, 128000.0, 134983.0, 55630.0, 9907.0, 1034.0, 88.0, 9.0, 1.0, 0.0, 10003.0, 7.0, 38.0, 448.0, 3775.0, 15886.0, 24921.0, 12281.0, 2351.0, 251.0, 34.0, 11.0, 1.0, 0.0, 2.0, 49.0, 748.0, 8005.0]
[-180.0, -170.0, -160.0, -150.0, -140.0, -130.0, -120.0, -110.0, -100.0, -90.0, -80.0, -70.0, -60.0, -50.0, -40.0, -30.0, -20.0, -10.0, 0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0, 120.0, 130.0, 140.0, 150.0, 160.0, 170.0, 180.0, 190.0, 200.0, 210.0, 220.0, 230.0, 240.0, 250.0, 260.0, 270.0, 280.0, 290.0, 300.0, 310.0, 320.0, 330.0, 340.0, 350.0, 360.0, 370.0, 380.0, 390.0, 400.0, 410.0, 420.0, 430.0, 440.0, 450.0, 460.0, 470.0, 480.0, 490.0, 500.0, 510.0, 520.0, 530.0]
[  1.06513910e-01   2.71470555e-01   3.00096813e-01   2.12290683e-01
   1.19744810e-01   6.39172666e-02   5.71562673e-02   9.79734529e-02
   1.83223081e-01   3.00488857e-01   3.16881933e-01   1.30595275e-01
   2.32573680e-02   2.42738655e-03   2.06586089e-04   2.11281228e-05
   2.34756920e-06   0.00000000e+00   2.34827347e-02   1.64329844e-05
   8.92076295e-05   1.05171100e-03   8.86207372e-03   3.72934843e-02
   5.85037720e-02   2.88304973e-02   5.51913519e-03   5.89239869e-04
   7.98173527e-05   2.58232612e-05   2.34756920e-06   0.00000000e+00
   4.69513840e-06   1.15030891e-04   1.75598176e-03   1.87922914e-02
   1.06513910e-01   2.71470555e-01   3.00096813e-01   2.12290683e-01
   1.19744810e-01   6.39172666e-02   5.71562673e-02   9.79734529e-02
   1.83223081e-01   3.00488857e-01   3.16881933e-01   1.30595275e-01
   2.32573680e-02   2.42738655e-03   2.06586089e-04   2.11281228e-05
   2.34756920e-06   0.00000000e+00   2.34827347e-02   1.64329844e-05
   8.92076295e-05   1.05171100e-03   8.86207372e-03   3.72934843e-02
   5.85037720e-02   2.88304973e-02   5.51913519e-03   5.89239869e-04
   7.98173527e-05   2.58232612e-05   2.34756920e-06   0.00000000e+00
   4.69513840e-06   1.15030891e-04   1.75598176e-03   1.87922914e-02]
72
72 100
0.180847112435 -0.0783847994574
0.182222386727 -0.196650902407
0.00250830488024 -0.00933496101084
2.67114356379e-05 -0.0347311437731
0.062767495087 -0.0111881763922
1.28560293503e-05 -4.8510744145e-06
0.220552337346 -0.0890237824048
0.182223441713 -0.196651185046
0.00250830527 -0.00933496111291
2.67114317111e-05 -0.0347311437753
0.0627674926005 -0.0111881671761
0.000341878967254 -0.000652789887686
[-165.0, -85.0, 5.0, 55.0, 195.0, 275.0, 365.0, 415.0]
[5.0, 55.0, 195.0, 275.0]
5.0
55.0
195.0
275.0






    

In [ ]:

    




    

In [ ]: