import pandas as pd
from pandas import DataFrame
from psychopy import data, core, gui, misc
import numpy as np
import seaborn as sns
#from ggplot import *
import scipy.stats as stats

import statsmodels.formula.api as smf
import statsmodels.api as sm
from __future__ import division

from pivottablejs import pivot_ui
%pylab inline

#plotting params
from matplotlib import rcParams
rcParams['font.family'] = 'ubuntu'
sns.set(style="whitegrid", color_codes=True)
sns.set_context("talk")

"""
Explanation: Analysis Gestalt appreciation & numerosity estimation experiment
Researchers: Elise Berbiers (MA student), Rebecca Chamberlain, Johan Wagemans, Sander Van de Cruys

Replication & extension of experiment 1 in: 

Topolinski, S., Erle, T. M., & Reber, R. (2015). Necker’s smile: Immediate affective consequences of early perceptual processes. Cognition, 140, 1–13.

Notes: 

All errorbars are 95% CIs
Only a 100 milliseconds exposure duration was used, instead of the two exposure durations (25 ms and 100 ms) in Topolinski, Erle, & Reber (2015). Note that Topolinski, Erle, & Reber found the effect for both exposure durations.
30 coherent "Gestalt" images, 30 scrambled versions, 23 inverted Gestalt images (7 images don't have a singular canonical orientation).

On the construction of the scrambled stimuli Topolinski & Strack write: 

"First, we developed and tested a set of pictorial stimuli that were useful for intuitive judgments because they were so degraded that they could only rarely be visually recognized (Bower et al., 1990). We used 30 black-and-white drawings of everyday objects randomly chosen from the inventory by Snodgrass and Vanderwart (1980), with the only constraint being that depicted objects were visually not too simple (e.g., a circle). Following Volz and von Cramon (2006), these stimuli were visually degraded by a filter that masked the black picture on the white background by increasing the white pixels by 75%. These pictures were the object condition (Volz & von Cramon, 2006) since they depicted visually degraded real objects. Then, these pictures were divided into nine equal rectangles (3x3); and these rectangles were randomly rotated within the picture (Volz & von Cramon, 2006; cf., Bower et al., 1990; Wippich, 1994). Thus, these pictures contained the same pixel information as in the object condition and even contained local collinearities (Volz & von Cramon, 2006), but the picture as a whole depicted a physically impossible and thus meaningless object. These pictures were used in the nonobject condition."



Interestingly, there is no effect of condition on density, but there is one on hull. This is probably due to the stimulus generation: the fact that the scrambled are made by scrambling 'tiles' of the source (which may keep the nearest dots together, and density measures the average distance to the nearest dots).


I agree we should address the zygomaticus data. Because the participants liked the inverted condition as much as the gestalt I think this is fairly easily addressed - we would predict the same pattern of results for inverted vs. scrambled as gestalt vs. scrambled (i.e. more positive affect). The authors state in the paper, 'We interpret this differential finding in the way that a success in early Gestalt integration immediately evokes an increase in positive affectivity’ (p.6). I think we could argue the same for the inverted condition, but it is low level grouping/clustering rather than higher-level gestalt that is producing the positive affect. Compact/dense figures are more easily processed, hence liked?


I also have some doubts on the way I analyzed it. Basically, I've done a simple anova of condition on ratings (as in Topolinski), where scramble is significant. Then, to look at the low-level measures, I computed mean rating for each image, and added variables for the extra measures of this image. In this case, an anova of condition on ratings is not significant.  If I add hull as predictor, then scramble becomes significant (0.045) and hull is too. If I add density (but not hull), there's no scramble effect, but there is a density effect. Finally, including all three leaves nothing significant.


At first glance I think the first approach (adding mean hull/density etc. for the mean image ratings) is more appropriate, but I agree we should check this with the lab. However, I’m not sure if adding them all in to the same ANCOVA is the right way to go. Aren’t there dependencies between some the measures (i.e. don’t we run into collinearity problems?). Should it then be three separate ANCOVAs with the image statistics as covariates?


If I now use the full dataset to add the low-level measures, so I just add the columns of density/hull to every existing trial, then scramble stays significant (note that it already was significant in this dataset), and additionally scramble*hull interaction and main effect of density become significant. So it seems these low-level measures are not sufficient/appropriate to account for the scramble-effect.


Overall I think there is a good story here but it’s worth presenting to the lab and checking:

How to analyse the effect of the image statistics
How to present the numerosity effects so they are in line with the story
End of explanation
"""


# get data file names
files = gui.fileOpenDlg("../data")

dfs = []

for filename in files:
    
    #print(filename)
    df = pd.read_table(filename, sep=",") #only Appreciation trials 85
   
    dfs.append(df)
    
    
df = pd.concat(dfs, ignore_index=True)    

len(df)

dfApp = df[pd.notnull(df.images)]       
dfNum = df[pd.notnull(df.imgA)]

print('total len:', len(dfApp))
print('total len:', len(dfNum))

"""
Explanation: Importing data
End of explanation
"""


dfApp = dfApp.rename(columns={'rating.response': 'rating'})
dfApp = dfApp.rename(columns={'rating.rt': 'rt'})

#add var for img
dfApp.loc[:,"img"]= dfApp.images.str.extract("(\d+)", expand=False).astype(int)


print(dfApp.columns)
dfApp.to_csv("dataLiking.csv")


dfApp = pd.read_csv("dataLiking.csv")

"""
Explanation: Appreciation experiment
Relevant vars:

condition
rating.response
rating.rt
End of explanation
"""


print(len(dfApp['participant'].unique()))
len(dfApp)

dfApp.groupby(["condition"]).rating.describe()

dfApp.groupby(["condition"]).rt.describe()

sns.distplot(dfApp.groupby(['participant']).rating.mean(), bins=20);

"""
Explanation: Descriptive stats
End of explanation
"""


#sns.violinplot(x="condition", y="rating.response", data=dfApp);
sns.pointplot(x="condition", y="rating", unit="participant", data=dfApp);


sns.violinplot(x="condition", y="rating", unit="participant", data=dfApp);

# GLM test

model = smf.glm(formula="rating ~ condition", data=dfApp)

results = model.fit()
print(results.summary())

"""
Explanation: Does gestaltness influence appreciation?
End of explanation
"""


# compute diff score to correlate with PNS score

print(len(dfApp['participant'].unique()))

def diffScore(df):
    
    gestaltm =  df[df.condition=='gestalt'].rating.mean()
    scrambledm = df[df.condition=='scrambled'].rating.mean()
    
    diff= gestaltm - scrambledm
    
    #df['id'].iloc[0]
    dfout = pd.DataFrame(data=[(gestaltm, scrambledm, diff)],\
                         columns=['gestaltm', 'scrambledm', 'diff'])
    
    return dfout

dfdiff = dfApp.groupby('participant').apply(diffScore)
dfdiff = dfdiff.reset_index()





# add PNS scores
dfPNS = pd.read_table("ScoringPNS.csv", sep=",") 
dfPNS = dfPNS.iloc[4:,:]
dfPNS["participant"] = pd.to_numeric(dfPNS["participant"])
dfmerged= pd.merge(dfdiff, dfPNS, how='outer', on='participant')
dfmerged.head()

"""
Explanation: Interim conclusion:

scrambling effect is significant, in expected direction (less liked).
Inversion does not seem to matter, meaning it probably is not an (implicitly processed) "familiar Gestalt" effect.
Note: no item effect included. Add it? (but then: incomplete data, bc inverted)

Correlation between individual difference scores (appreciation of Gestalt minus scrambled patterns) and PNS scores
End of explanation
"""


sns.jointplot(x="Totaalscore", y="diff", data=dfmerged, kind="reg");
dfmerged["diff"].mean()

"""
Explanation: Does higher PNS (personal need for structure) imply higher relative appreciation of Gestalt vs scrambled patterns?
End of explanation
"""


dfNum = dfNum.rename(columns={'numResp.corr': 'acc','numResp.rt': 'rt', 'numResp.Alocation': 'Alocation' })
#add var for img
dfNum.loc[:,"img"]= dfNum.imgA.str.extract("(\d+)", expand=False).astype(int)

dfNum.to_csv('dfNum.csv', sep='\t')

dfNum = pd.read_csv("dfNum.csv", sep='\t')

print dfNum.columns
print dfNum.dtypes

"""
Explanation: Interim conclusion

Nonsign correlation
Trend in direction opposite to expected: higher PNS -> less liking of Gestalt vs scrambled.
This outcome is what would be expected if the "Gestalt-effect" is not actually caused by coherent Gestalt, but by more basic image feature not related to PNS?
main effect on rating?

Numerosity experiment

2-IFC: which has the most dots (first or second image)
Note that all comparison images have exactly the same numbers of dots in reality
Hence any deviation from .50 is caused by the difference in configuration

Relevant vars:

acc: percentage choice for first image
rt
condition
End of explanation
"""


sns.pointplot(x="condition", y="acc",unit="participant", data=dfNum);
#sns.axlabel("Condition", "Percentage higher")

#sns.stripplot(x="condition", y="acc",unit="participant", data=dfNum, jitter=True);

model = smf.glm(formula="acc ~ condition", data=dfNum,  family=sm.families.Binomial())

results = model.fit()
print(results.summary())

dfimnum = dfNum.groupby(['imgA','imgB','condition'])['acc'].mean()
dfimnum= dfimnum.reset_index()
dfimnum.head()

"""
Explanation: Does Gestaltness (and inversion) influence numerosity estimation?
End of explanation
"""


# prerequisites

import os, re
import matplotlib.pyplot as plt
from skimage.data import data_dir
from skimage.util import img_as_ubyte
from skimage.morphology import skeletonize, convex_hull_image, disk
from skimage import io
from skimage.measure import label, regionprops
from skimage.filters.rank import entropy
from skimage.feature import peak_local_max



def plot_comparison(original, filtered, filter_name):

    fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(8, 4), sharex=True,
                                   sharey=True)
    ax1.imshow(original, cmap=plt.cm.gray)
    ax1.set_title('original')
    ax1.axis('off')
    ax1.set_adjustable('box-forced')
    ax2.imshow(filtered, cmap=plt.cm.gray)
    ax2.set_title(filter_name)
    ax2.axis('off')
    ax2.set_adjustable('box-forced')


"""
Explanation: Notes

The numerosity exp might be nice to add, because it is consistent with the liking results. Basically, people noticed the equal number of dots for the comparison Gestalt vs inverted (it's around 50%, where it theoretically always should be because the compared images always contained an equal amount of dots). In the other comparisons you see they (substantially) overestimate the dots in the scrambled condition (or underestimate the Gestalt/inverted condition). 

Checking for low-level structural measures of the test images
Using scikit-image:

Stéfan van der Walt, Johannes L. Schönberger, Juan Nunez-Iglesias, François Boulogne, Joshua D. Warner, Neil Yager, Emmanuelle Gouillart, Tony Yu and the scikit-image contributors. scikit-image: Image processing in Python. PeerJ 2:e453 (2014) http://dx.doi.org/10.7717/peerj.453
End of explanation
"""


from itertools import repeat

def distance(p1,p2):
    """Euclidean distance between two points."""
    x1,y1 = p1
    x2,y2 = p2
    return np.hypot(x2 - x1, y2 - y1)


# get data file names
files = gui.fileOpenDlg("../ExpCode/images/)")#, allowed="bmp files (*.bmp)|*.bmp")

d = []

for filename in files:
    
    #print(filename)
    stim = io.imread(filename, as_grey=True)
    stim = img_as_ubyte(io.imread(filename, as_grey=True))
    #fig, ax = plt.subplots()
    #ax.imshow(stim, cmap=plt.cm.gray)
    
    #compute convex hull for this img
    
    chull = convex_hull_image(stim == 0)
    #plot_comparison(stim, chull, 'convex hull')
    label_img = label(chull)
    regions_img = regionprops(label_img)
    region = regions_img[0]

    hull= float(region.convex_area)/(chull.shape[0]*chull.shape[1])
    #print "percentage pixels of convex hull image: ", hull
    
    #compute density for this img
    
    stim = numpy.invert(stim) # peaks should be white
    # extract peaks aka dots
    coords = peak_local_max(stim, min_distance=1) 
    #print coords
    #print len(coords)
    #print len(stim[stim>.5])
    
    #compute density as mean of distances to 3 closest neighbors
    
    density = 0
    for p1 in coords:
        
        dists = [distance(*pair) for pair in zip(repeat(p1),coords)]
        sorted_dists = np.sort(dists)[1:6] # take only distances to 3 closest neighbors 
        av_dist = np.mean(sorted_dists)
        density += av_dist
        #print "dists: ", len(dists)
        #print "sorted: ", sorted_dists
        #print "average: ", av_dist 

    density /= len(coords)  
    
    d.append({'images': filename.split("GestaltAppreciation/ExpCode/")[1] ,'hull': hull, 'density': 1/density})

#print(d)                 
stims = pd.DataFrame(d)


# make dataset where row is image

dfim = dfApp.groupby(['images','condition'])['rating'].mean()
dfim= dfim.reset_index()
dfim.head()

dfmerged= pd.merge(dfim, stims, how='outer', on='images')
dfmerged.head()

dfmerged.to_csv("dataLikingWithLowlevel.csv")

dfApp= pd.merge(dfApp, stims, how='outer', on='images')

# GLM test
model = smf.glm(formula="rating ~ condition * hull * density", data=dfApp)
results = model.fit()
print(results.summary())

"""
Explanation: Compute measures: convex hull area and density
Convex hull area

Convex hull docs: putting an elastic band around all elements to form the smallest possible convex polygon
Intuitively the area depends on how much the points are spread out.
One possible measure of clustering/compactness

Density

For each dot: mean of distances to 3 closest neighbors, averaged over all dots
End of explanation
"""


dfnummerged= pd.merge(dfimnum, stims, how='left', left_on='imgA', right_on='images')
dfnummerged = dfnummerged.rename(columns={'hull': 'Ahull','density': 'Adensity', 'images': 'imageA' })
dfnummerged= pd.merge(dfnummerged, stims, how='left', left_on='imgB', right_on='images')
dfnummerged = dfnummerged.rename(columns={'hull': 'Bhull','density': 'Bdensity', 'images': 'imageB' })
dfnummerged["hullDiff"]= dfnummerged.Ahull-dfnummerged.Bhull
dfnummerged["densDiff"]= dfnummerged.Adensity-dfnummerged.Bdensity
dfnummerged["img"]= dfnummerged['imgA'].str.extract("(\d+)", expand=False).astype(int)

dfnummerged.tail()      

lldf = dfnummerged[dfnummerged.condition == "Gestalt-scrambled"]
lldf = lldf[['img','hullDiff', 'densDiff']]
lldf.head()


model = smf.glm(formula="acc ~ hullDiff*condition", data=dfnummerged, family=sm.families.Binomial())

results = model.fit()
print(results.summary())

"""
Explanation: Do low-level measures influence on numerosity judgments?
End of explanation
"""



g = sns.distplot(dfmerged[dfmerged.condition=="inverted"].hull, bins=15, label="inverted");
sns.distplot(dfmerged[dfmerged.condition=="scrambled"].hull, bins=15, label="scrambled");
sns.distplot(dfmerged[dfmerged.condition=="gestalt"].hull, bins=15, label="Gestalt");
plt.legend();

g.set(xlabel='Convex hull area', ylabel='Frequency');


"""
Explanation: How is convex hull area distributed? (over images)
End of explanation
"""


g = sns.pointplot(x="condition", y="hull", data=dfmerged, linestyles=["-"]);
g.set(xlabel='Condition', ylabel='Average convex hull area');

# GLM test

model = smf.glm(formula="hull ~ condition", data=dfmerged)
#model = smf.ols(formula="rating ~ hull", data=dfmerged)

results = model.fit()
print(results.summary())

"""
Explanation: Does convex hull area differ across conditions?
End of explanation
"""


sns.set_style("whitegrid",{"axes.grid": False})
sns.set_context("talk")

g = sns.lmplot(x="hull", y="rating", hue="condition",\
           data=dfmerged, x_estimator=np.mean, palette="Set1");
g.set(xlabel='Convex hull area', ylabel='Rating');
sns.set_style("whitegrid",{"axes.grid": False})

# GLM test

model = smf.glm(formula="rating ~ condition * hull", data=dfmerged)
#model = smf.ols(formula="rating ~ hull", data=dfmerged)

results = model.fit()
print(results.summary())

"""
Explanation: Interim conclusion

As expected, convex hull area is significantly higher in the scrambled condition (compared to Gestalt or inverted)
Inversion does not matter (obviously). It is not exactly the same as Gestalt because some images were not inverted (see notes at top)

Are ratings dependent on convex hull area, in the different conditions?
End of explanation
"""



sns.distplot(dfmerged[dfmerged.condition=="inverted"].density, bins=15, label="inverted");
sns.distplot(dfmerged[dfmerged.condition=="scrambled"].density, bins=15, label="scrambled");
sns.distplot(dfmerged[dfmerged.condition=="gestalt"].density, bins=15, label="Gestalt");
plt.legend();
g.set(xlabel='Density', ylabel='Frequency');


"""
Explanation: Interim conclusion:

If we add hull area to the model (and interaction with condition) for DV rating:
Hull significantly influences ratings
Interestingly the main scrambling effect also becomes significant here
Interaction scramblingxhull is also (marginally) significant. As can be seen in the graph, this is because within the scrambling condition, hull does not seem to influence ratings much. This is only the case for the inverted and the Gestalt condition, specifically, both show a (similar) negative correlation.
Together this seems to indicate that ratings are determined by low-level structural differences rather than high-level Gestalt or coherence.

How is density distributed? (over images)
End of explanation
"""


sns.pointplot(x="condition", y="density", data=dfmerged, linestyles=["-"]);
g.set(xlabel='Condition', ylabel='Average density');

# GLM test

model = smf.glm(formula="density ~ condition", data=dfmerged)
#model = smf.ols(formula="rating ~ hull", data=dfmerged)

results = model.fit()
print(results.summary())

"""
Explanation: Does density differ across conditions?
End of explanation
"""


sns.set_style("whitegrid",{"axes.grid": False})
sns.set_context("talk")

g= sns.lmplot(x="density", y="rating",\
              hue="condition", data=dfmerged, x_estimator=np.mean, palette="Set1");
g.set(xlabel='Density', ylabel='Rating');


# GLM test

model = smf.glm(formula="rating ~ condition * density", data=dfmerged)
#model = smf.ols(formula="rating ~ hull", data=dfmerged)

results = model.fit()
print(results.summary())

# GLM test

model = smf.glm(formula="rating ~ condition * hull * density", data=dfmerged)
#model = smf.ols(formula="rating ~ hull", data=dfmerged)

results = model.fit()
print(results.summary())

"""
Explanation: Interim conclusion

Density differs in the three conditions.
Inverted imgs should have same density but as Gestalt but because some images were not inverted this is not the case (see notes at top)

Are ratings dependent on density, in the different conditions?
End of explanation
"""


sns.jointplot(x="density", y="hull", data=dfmerged, kind="reg")

"""
Explanation: Interim conclusion:

If we add density to the model (and interaction with condition) for DV rating:
Density significantly influences ratings (higher density -> higher rating)
No other significant relations.

Is there a correlation between density and convex hull area? (Yes, a negative one, as expected)
End of explanation
"""


ndf = dfNum.pivot_table(index=['participant','img'], columns=['condition'], values='acc')
ndf = ndf.reset_index()
adf = dfApp.pivot_table(index=['participant','img'], columns=['condition'], values='rating')
adf = adf.reset_index()
adf["diffGestalt"] = adf.gestalt - adf.scrambled
adf["diffInv"] = adf.inverted - adf.scrambled
andf = pd.merge(ndf, adf, how='outer', on=['participant', 'img'])
andf = pd.merge(andf, lldf, how = 'left', on = 'img')
andf.head()

andf.to_csv("data_per_img_participant.csv")

andf = pd.read_csv("data_per_img_participant.csv")

stats.ttest_1samp(andf[andf.diffInv.notnull()].diffInv, 0, axis=0)
stats.ttest_1samp(andf.diffGestalt, 0, axis=0)

sns.jointplot(x="hullDiff", y="diffGestalt", data=andf, kind="reg");

sns.lmplot(x="hullDiff", y="diffGestalt", hue="participant", units="participant", 
           y_jitter=0.5, x_jitter=0.02, data=andf, ci=False, fit_reg=False, legend=False, scatter_kws={'alpha':0.2}, line_kws={'alpha':0.2});

sns.lmplot(x="hullDiff", y="diffGestalt", hue="participant",\
           scatter=False, legend=False, ci=False, data=andf, line_kws={'alpha':0.2});

"""
Explanation: Does difference in hull/density determines difference in liking on a per-image basis?
Make combined dataset:
End of explanation
"""


partHigh = andf.groupby('participant').filter(lambda x: x['diffGestalt'].mean() >= .3)
partLow = andf.groupby('participant').filter(lambda x: x['diffGestalt'].mean() < .1)
print(len(andf))
print(len(partHigh))
len(partHigh.participant.unique())

sns.lmplot(x="hullDiff", y="diffGestalt", hue="participant",\
           scatter=False, legend=False, ci=False, data=partHigh, line_kws={'alpha':0.3}, scatter_kws={'alpha':0.3});


model = smf.glm(formula="diffGestalt ~ hullDiff", data=andf)

results = model.fit()
print(results.summary())

"""
Explanation: Are participants with greatest effect also those with greatest effect of hulldiff/densdiff?
Select participants with greatest average effect (over images):
End of explanation
"""


sns.pointplot(x="Gestalt-scrambled", y="diffGestalt", unit="participant", data=andf);

sns.pointplot(x="Inverted-scrambled", y="diffInv", unit="participant", data=andf);

stats.ttest_ind(andf.diffGestalt[andf["Gestalt-scrambled"]==0], 
                andf.diffGestalt[andf["Gestalt-scrambled"]==1])

model = smf.glm(formula="diffGestalt ~ hullDiff", data=andf)

results = model.fit()
print(results.summary())

"""
Explanation: Does numerosity predict liking?
For those that overestimated the number of dots in the scrambled condition, did they like the scrambled images less?
End of explanation
"""


