In [1]:
import numpy as np
import numpy.ma as ma
import scipy.stats as stat
import random

In [ ]:
SM2xy = np.genfromtxt('SM2xy_RED_nPEG_37C_72pH_S1.csv', delimiter=",")

In [ ]:
SM2xy[10, 100:200]

In [ ]:
x10 = ma.log(SM2xy[10, :]) #take the log of the MSDS at a single timpoint

In [ ]:
x= np.mean(x10) #take the mean of the logs (geometric mean)
print(x)

In [ ]:
y= stat.sem(x10) #find the standard error around the mean of the logs (exponentiate to get confidence intervals)

In [ ]:
random.random()

In [ ]:
#an example, pretending to find weighted averages over 11 videos
SEM = np.zeros([1, 10])
theta = np.zeros([1, 10])
w = np.zeros([1, 10])

#create some fake averages and SEMs for demonstration purposes
for i in range(0, 11):
    theta[0, i - 1] = x + random.random()
    SEM[0, i - 1] = y + random.random()

wi = np.sum(1./(SEM*SEM))

for i in range(0, 11):
    w[0, i -1] = (1/(SEM[0, i-1]*SEM[0, i-1]))/wi

In [ ]:
otheta = np.sum(w*theta)

In [ ]:
oSEM = np.sqrt(1/(np.sum(SEM*SEM)))

In [ ]:
oSEM

In [ ]:
#where real calculations begin, not limited to one row.

#performing averages over an entire dataset of MSDs at all timepoints
geoMSD = np.mean(ma.log(SM2xy), axis=1)

In [ ]:
geoMSD

In [ ]:
gSEM = stat.sem(ma.log(SM2xy), axis=1)

In [ ]:
frames = 651

SEM = np.zeros([frames, 10])
theta = np.zeros([frames, 10])
w = np.zeros([frames, 10])

for i in range(0, 11):
    theta[:, i-1] = geoMSD + np.random.rand(frames)*0.01
    SEM[:, i-1] = gSEM + np.random.rand(frames)*0.0001

wi = np.sum((1./(SEM*SEM)), axis=1)

for i in range(0, 11):
    w[:, i-1] = (1/(SEM[:, i-1]*SEM[:, i-1]))/wi

In [ ]:
otheta = np.sum(w*theta, axis=1)
oSEM = np.sqrt(1/(np.sum(1/(SEM*SEM), axis=1)))

In [ ]:
oetheta = np.exp(otheta)

In [ ]:
loetheta = np.exp(otheta - 1.96*oSEM)
hoetheta = np.exp(otheta + 1.96*oSEM)

In [ ]:
loetheta

In [ ]:


In [15]:
parameters = {}
parameters["channels"] = ["RED"]
parameters["surface functionalities"] = ["nPEG", "PEG"]
parameters["slices"] = ["S1", "S2", "S3"]
parameters["videos"] = [1, 2, 3, 4, 5]

channels = parameters["channels"]
surface_functionalities = parameters["surface functionalities"]
slices = parameters["slices"]
videos = parameters["videos"]

geoM2xy = {}
gSEM = {}
SM1x = {}
SM1y = {}
SM2xy = {}
npar = {}

In [16]:
DIR = "./10mM/redo/"
cond = '37C_pH72'

for channel in channels:
    for surface_functionality in surface_functionalities:
        slice_counter = 0
        for slic in slices:
            for video in videos:
                sample_name = "{}_{}_{}_{}_{}".format(channel, surface_functionality, cond, slic, video)                
                #SM2xy[sample_name] = np.genfromtxt('SM2xy_{}.csv'.format(sample_name, delimiter=","))
                
                #npar[sample_name] = SM2xy[sample_name].shape
                geoM2xy[sample_name] = np.genfromtxt(DIR + 'geoM2xy_{}.csv'.format(sample_name, delimiter=","))
                gSEM[sample_name] = np.genfromtxt(DIR + 'gSEM_{}.csv'.format(sample_name, delimiter=","))

In [17]:
geo_slices = {}
gSEM_slices = {}
w_slices = {}
wo_slices = {}

In [18]:
#Calculate the precision weights over videos
#Alternately, can weight by the number of particles in each video
nvids = 5
nslices = 4
frames = 651

for channel in channels:
    for surface_functionality in surface_functionalities:
        slice_counter = 0
        for slic in slices:
            video_counter = 0
            w_holder = np.zeros((nvids, frames))
            sample_name = "{}_{}_{}_{}".format(channel, surface_functionality, cond, slic)
            for key in geoM2xy:
                if sample_name in key:
                    w_holder[video_counter, :] = 1/(gSEM[key]*gSEM[key])
                    video_counter = video_counter + 1
            wo_slices[sample_name] = np.sum(w_holder, axis=0)
            slice_counter = slice_counter + 1


C:\Users\koolk\Miniconda3\lib\site-packages\ipykernel\__main__.py:16: RuntimeWarning: divide by zero encountered in true_divide

In [19]:
#Calculate the weights SEMs and means over videos
#Remember to use alternate if not wanting to use precision weights at this level.
for channel in channels:
    for surface_functionality in surface_functionalities:
        slice_counter = 0
        for slic in slices:
            geo_holder = np.zeros((nvids, frames))
            gSEM_holder = np.zeros((nvids, frames))
            w_holder = np.zeros((nvids, frames))
            video_counter = 0
            sample_name = "{}_{}_{}_{}".format(channel, surface_functionality, cond, slic)
            for key in geoM2xy:
                if sample_name in key:
                    w_holder[video_counter, :] = (1/(gSEM[key]*gSEM[key]))/wo_slices[sample_name]
                    geo_holder[video_counter, :] = w_holder[video_counter, :] * geoM2xy[key]
                    gSEM_holder[video_counter, :] = (1/(gSEM[key]*gSEM[key]))
                    video_counter = video_counter + 1
            geo_slices[sample_name] = np.sum(geo_holder, axis=0)
            gSEM_slices[sample_name] = np.sqrt((1/np.sum(gSEM_holder, axis=0)))
            slice_counter = slice_counter + 1


C:\Users\koolk\Miniconda3\lib\site-packages\ipykernel\__main__.py:14: RuntimeWarning: divide by zero encountered in true_divide
C:\Users\koolk\Miniconda3\lib\site-packages\ipykernel\__main__.py:14: RuntimeWarning: invalid value encountered in true_divide
C:\Users\koolk\Miniconda3\lib\site-packages\ipykernel\__main__.py:16: RuntimeWarning: divide by zero encountered in true_divide

In [20]:
geo = {}
gS = {}
w_slices = {}
wo_slices = {}

In [21]:
#Calculate the precision weights over slices
for channel in channels:
    counter = 0
    for surface_functionality in surface_functionalities:
        w_holder = np.zeros((nslices, frames))
        slice_counter = 0
        sample_name = "{}_{}_{}".format(channel, surface_functionality, cond)
        for key in geo_slices:
            if sample_name in key:
                w_holder[slice_counter, :] = 1/(gSEM_slices[key]*gSEM_slices[key])
                slice_counter = slice_counter + 1
        wo_slices[sample_name] = np.sum(w_holder, axis=0)
        counter = counter + 1


C:\Users\koolk\Miniconda3\lib\site-packages\ipykernel\__main__.py:10: RuntimeWarning: divide by zero encountered in true_divide

In [22]:
#Calculate the weights SEMs and means over slices
for channel in channels:
    counter = 0
    for surface_functionality in surface_functionalities:
        geo_holder = np.zeros((nslices, frames))
        gSEM_holder = np.zeros((nslices, frames))
        w_holder = np.zeros((nslices, frames))
        slice_counter = 0
        sample_name = "{}_{}_{}".format(channel, surface_functionality, cond)
        for key in geo_slices:
            if sample_name in key:
                w_holder[slice_counter, :] = (1/(gSEM_slices[key]*gSEM_slices[key]))/wo_slices[sample_name]
                geo_holder[slice_counter, :] = w_holder[slice_counter, :] * geo_slices[key]
                gSEM_holder[slice_counter, :] = (1/(gSEM_slices[key]*gSEM_slices[key]))
                slice_counter = slice_counter + 1
        geo[sample_name] = np.sum(geo_holder, axis=0)
        gS[sample_name] = np.sqrt((1/np.sum(gSEM_holder, axis=0)))
        counter = counter + 1


C:\Users\koolk\Miniconda3\lib\site-packages\ipykernel\__main__.py:12: RuntimeWarning: divide by zero encountered in true_divide
C:\Users\koolk\Miniconda3\lib\site-packages\ipykernel\__main__.py:12: RuntimeWarning: invalid value encountered in true_divide
C:\Users\koolk\Miniconda3\lib\site-packages\ipykernel\__main__.py:14: RuntimeWarning: divide by zero encountered in true_divide

In [26]:
np.exp(geo['RED_PEG_37C_pH72'])


Out[26]:
array([        nan,  0.13670857,  0.19321679,  0.2444675 ,  0.28632928,
        0.32309273,  0.35631178,  0.38686723,  0.41561198,  0.44259056,
        0.46831452,  0.49278114,  0.51637739,  0.53924429,  0.56128842,
        0.58302006,  0.60385827,  0.62442879,  0.64419151,  0.66317858,
        0.68146361,  0.70040877,  0.71865672,  0.73649938,  0.75394626,
        0.7713587 ,  0.78822224,  0.80421202,  0.82104115,  0.83678071,
        0.85358275,  0.86902333,  0.88431794,  0.89930005,  0.91472886,
        0.92847019,  0.9434108 ,  0.95776271,  0.97162306,  0.9841274 ,
        0.99782002,  1.01248894,  1.02612868,  1.03963354,  1.05440659,
        1.06688431,  1.07871468,  1.09140542,  1.10490419,  1.11758375,
        1.12921682,  1.14017415,  1.15191053,  1.16436164,  1.17459904,
        1.18548434,  1.19562437,  1.20749267,  1.21831019,  1.22884059,
        1.23857241,  1.24966727,  1.25990564,  1.27173169,  1.28084883,
        1.28960888,  1.29992384,  1.30973131,  1.31722999,  1.32864683,
        1.33543382,  1.34462269,  1.35256304,  1.36218063,  1.37009912,
        1.37887477,  1.38673603,  1.39495443,  1.39973281,  1.40753539,
        1.41638563,  1.42698317,  1.43331034,  1.44210258,  1.45135655,
        1.45918112,  1.46569584,  1.46868695,  1.47614179,  1.48187382,
        1.49025777,  1.49868802,  1.50568632,  1.51175486,  1.51825041,
        1.52408895,  1.53081092,  1.53643984,  1.54212806,  1.54808375,
        1.55343809,  1.55835667,  1.56189448,  1.56523192,  1.57079086,
        1.57436885,  1.58148782,  1.58636625,  1.5858371 ,  1.59181646,
        1.59352512,  1.59792195,  1.60275209,  1.60672967,  1.60678843,
        1.60939049,  1.61083873,  1.61770021,  1.62054035,  1.62518401,
        1.62599859,  1.62858037,  1.63208439,  1.63151563,  1.63193159,
        1.63684923,  1.63839474,  1.64000255,  1.63811443,  1.64353697,
        1.6425965 ,  1.64481529,  1.6408682 ,  1.64129664,  1.64003765,
        1.6440806 ,  1.64383788,  1.64129237,  1.642678  ,  1.6423004 ,
        1.64161355,  1.64106262,  1.64077691,  1.64262023,  1.64165227,
        1.64587172,  1.64682636,  1.63937126,  1.64255216,  1.63772312,
        1.64097208,  1.64008769,  1.63707222,  1.62915868,  1.63172395,
        1.62816198,  1.62524693,  1.62419166,  1.62168551,  1.62230701,
        1.61437608,  1.61065677,  1.60358086,  1.60592694,  1.59806196,
        1.59757608,  1.59281098,  1.58432687,  1.58057024,  1.57617109,
        1.56425625,  1.55782891,  1.55809409,  1.54818665,  1.54137356,
        1.53781771,  1.52327975,  1.51910078,  1.51640504,  1.51177171,
        1.50445609,  1.49846929,  1.4938189 ,  1.48900958,  1.48368963,
        1.48110811,  1.47698412,  1.47633139,  1.46129703,  1.44936435,
        1.44642656,  1.44689899,  1.4358171 ,  1.42745985,  1.41739549,
        1.41737929,  1.40921906,  1.3965389 ,  1.38507094,  1.37372769,
        1.37047617,  1.36566732,  1.36046584,  1.34890452,  1.33608414,
        1.33120465,  1.32192363,  1.31285755,  1.32105018,  1.3089774 ,
        1.29509276,  1.28480119,  1.27582283,  1.26630988,  1.25959634,
        1.24797662,  1.24360877,  1.24495925,  1.23463063,  1.23192054,
        1.22115238,  1.21363213,  1.19798886,  1.18891238,  1.18048115,
        1.17533502,  1.16597028,  1.16585512,  1.15192212,  1.13811564,
        1.12823712,  1.11736587,  1.10490742,  1.10393128,  1.0947495 ,
        1.08748349,  1.08325151,  1.06635262,  1.05770622,  1.05196743,
        1.0452512 ,  1.03531401,  1.02726311,  1.02239972,  1.01815377,
        1.01424923,  0.99124849,  0.98237769,  0.96519971,  0.96014199,
        0.94057124,  0.94208805,  0.92817064,  0.92955209,  0.90713151,
        0.91209058,  0.90139438,  0.89508904,  0.89622415,  0.87282731,
        0.86800964,  0.85308997,  0.84646711,  0.84785591,  0.84496597,
        0.83840946,  0.82320815,  0.82013873,  0.80914126,  0.79567523,
        0.77903837,  0.78182309,  0.75919962,  0.75830077,  0.7468711 ,
        0.73699138,  0.73029254,  0.72998439,  0.71440989,  0.70075438,
        0.70709335,  0.69576093,  0.67344167,  0.67051203,  0.65744525,
        0.6572111 ,  0.66171674,  0.64107865,  0.64205485,  0.62792506,
        0.62413666,  0.60893319,  0.61529163,  0.59988716,  0.5907089 ,
        0.57952831,  0.57826067,  0.57157903,  0.57047163,  0.56217544,
        0.54595113,  0.54371161,  0.54008441,  0.53485626,  0.51831994,
        0.51976325,  0.51493067,  0.50334024,  0.49635081,  0.4921304 ,
        0.48108684,  0.48298511,  0.47122721,  0.47217111,  0.46372013,
        0.45773019,  0.44480953,  0.44510708,  0.43435613,  0.42574079,
        0.4264631 ,  0.41612462,  0.42129178,  0.42694167,  0.41303316,
        0.41954474,  0.40848485,  0.40373389,  0.39829441,  0.40268585,
        0.39711345,  0.39324085,  0.38304264,  0.38607708,  0.36928146,
        0.36164261,  0.35718964,  0.3558992 ,  0.35693844,  0.35335277,
        0.34760298,  0.34251258,  0.34092634,  0.33510096,  0.33258103,
        0.32973265,  0.32180458,  0.32438427,  0.31577121,  0.31328199,
        0.31542571,  0.30510102,  0.29692796,  0.29375057,  0.29408902,
        0.2843268 ,  0.282707  ,  0.28331609,  0.28118975,  0.28895289,
        0.28966227,  0.28427088,  0.28233683,  0.27929427,  0.27574279,
        0.27236342,  0.27191776,  0.26209909,  0.26859924,  0.25790991,
        0.25411939,  0.25005041,  0.25117471,  0.24843479,  0.25276848,
        0.2442414 ,  0.23512936,  0.2380273 ,  0.2283052 ,  0.22442217,
        0.22906747,  0.22451806,  0.2183545 ,  0.21713141,  0.21117267,
        0.21030044,  0.20301453,  0.20646321,  0.2057766 ,  0.19999596,
        0.19590724,  0.19380491,  0.18920863,  0.19209118,  0.19041206,
        0.18680251,  0.18616957,  0.17727696,  0.17774812,  0.18391369,
        0.18097748,  0.17355315,  0.17352287,  0.16465162,  0.1619001 ,
        0.15698626,  0.15908503,  0.15764792,  0.15523537,  0.14889705,
        0.14218063,  0.14283155,  0.14059172,  0.14463271,  0.14088465,
        0.13921936,  0.13511071,  0.13870755,  0.13531037,  0.13783559,
        0.13933063,  0.13373896,  0.13030483,  0.12922844,  0.12281846,
        0.12086341,  0.11659694,  0.11812624,  0.11348518,  0.11244042,
        0.11063644,  0.11323213,  0.10628357,  0.10364105,  0.10292453,
        0.10365815,  0.10076164,  0.09761793,  0.10128539,  0.10112561,
        0.09731407,  0.10057195,  0.09581001,  0.09455015,  0.09659836,
        0.09728287,  0.09712462,  0.09423086,  0.09532476,  0.09120938,
        0.0886034 ,  0.0880975 ,  0.09008021,  0.08646673,  0.08390762,
        0.08092832,  0.08157348,  0.08227306,  0.07794486,  0.07525438,
        0.07416916,  0.07299046,  0.07531351,  0.07317646,  0.06837008,
        0.06898513,  0.06733015,  0.06573817,  0.06570333,  0.06320009,
        0.06568534,  0.06276663,  0.0639829 ,  0.06267979,  0.06222481,
        0.0592191 ,  0.06100559,  0.05989935,  0.05964053,  0.06140024,
        0.05844104,  0.05732757,  0.05686968,  0.05544604,  0.05492326,
        0.056406  ,  0.05702772,  0.05532687,  0.05538327,  0.0524982 ,
        0.05398283,  0.05363615,  0.05228956,  0.05483727,  0.05292606,
        0.05110442,  0.05128705,  0.05056023,  0.04986873,  0.0509998 ,
        0.0454241 ,  0.04610519,  0.04436916,  0.04509698,  0.04744068,
        0.04728795,  0.04660616,  0.04793482,  0.04623789,  0.04749428,
        0.0445956 ,  0.04532226,  0.04363784,  0.04458074,  0.04299269,
        0.04142843,  0.04010611,  0.04098117,  0.04113286,  0.039689  ,
        0.03920625,  0.03805561,  0.03798202,  0.03773647,  0.03755681,
        0.03500887,  0.03470736,  0.03663853,  0.03383548,  0.03381696,
        0.03432824,  0.03402318,  0.03294947,  0.03220661,  0.03292106,
        0.03229529,  0.03355765,  0.03164112,  0.03087207,  0.03107044,
        0.02947993,  0.03136507,  0.03082258,  0.02959129,  0.03070677,
        0.02975417,  0.02933892,  0.02871981,  0.02937343,  0.02948654,
        0.03061156,  0.02939322,  0.02999776,  0.03012016,  0.03016618,
        0.029304  ,  0.02947728,  0.02809732,  0.02796792,  0.02707366,
        0.02710942,  0.02649647,  0.02539036,  0.02589228,  0.02627013,
        0.0265808 ,  0.02480714,  0.02623182,  0.02457979,  0.02522502,
        0.02386945,  0.02478408,  0.02435383,  0.02451832,  0.02455294,
        0.02511304,  0.02424592,  0.02395895,  0.02463111,  0.02376692,
        0.02400483,  0.02502662,  0.02338679,  0.02360558,  0.02197604,
        0.02212892,  0.02297961,  0.02193878,  0.02216552,  0.02159748,
        0.02088463,  0.020747  ,  0.02042062,  0.01987623,  0.01954617,
        0.01873012,  0.01858156,  0.01848652,  0.0191478 ,  0.01845855,
        0.01839952,  0.01930515,  0.01872719,  0.01832982,  0.01963752,
        0.01896712,  0.01850039,  0.01915982,  0.01819473,  0.01737948,
        0.01825058,  0.01829619,  0.01803805,  0.01941819,  0.01879071,
        0.01781747,  0.01697377,  0.01795109,  0.01691593,  0.01676375,
        0.01671739,  0.01682227,  0.01877459,  0.01675345,  0.01578177,
        0.01704496,  0.01628361,  0.0176589 ,  0.01717665,  0.01728755,
        0.01633907,  0.01728263,  0.01619983,  0.01520867,  0.01608906,
        0.01610706,  0.01547745,  0.01697538,  0.01600066,  0.01587225,
        0.01427545,  0.01446711,  0.01440449,  0.01413024,  0.01371949,
        0.01413985,  0.01286495,  0.01355754,  0.01230611,  0.01190423,
        0.00890605])

In [25]:
for key in geo:
    np.savetxt(DIR+'geoM2xy_{}.csv'.format(key), geo[key], delimiter=',')
    np.savetxt(DIR+'gSEM_{}.csv'.format(key), gS[key], delimiter=',')

In [ ]: