move from pystan to cmdstanpy

Author: ihrke

docs/modeling.ipynb

@@ -1776,7 +1776,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.6"
+   "version": "3.8.3"
   },
   "rise": {
    "enable_chalkboard": true,

docs/symp_talk_uit2019.ipynb

@@ -8,6 +8,8 @@
 import scipy.signal as signal
 import scipy
 import math
+import os
+import cmdstanpy
 
 import scipy.interpolate
 from scipy.interpolate import interp1d, splrep, splev
@@ -290,7 +292,11 @@ def prominence_to_lambda(w, lam_min=1, lam_max=100):
     # load or compile model
     vprint(10, "Compiling Stan model")
 
-    sm = StanModel_cache(stan_code_baseline_model_asym_laplac)
+    fname="stan/baseline_model_asym_laplac.stan"
+    fpath=os.path.join(os.path.split(__file__)[0], fname)
+    sm = cmdstanpy.CmdStanModel(stan_file=fpath)
+    sm.compile()
+    #sm = StanModel_cache(stan_code_baseline_model_asym_laplac)
 
     ## put the data for the model together
     data={
@@ -307,8 +313,8 @@ def prominence_to_lambda(w, lam_min=1, lam_max=100):
 
     ## variational optimization
     vprint(10, "Optimizing Stan model")
-    opt=sm.vb(data=data)
-    vbc=opt["mean_pars"]
+    opt=sm.variational(data=data)
+    vbc=opt.stan_variable("coef")
     meansigvb=np.dot(B, vbc)
     vprint(10, "Done optimizing Stan model")
 
@@ -350,8 +356,8 @@ def prominence_to_lambda(w, lam_min=1, lam_max=100):
 
     ##  variational optimization
     vprint(10, "2nd Optimizing Stan model")
-    opt=sm.vb(data=data2)
-    vbc2=opt["mean_pars"]
+    opt=sm.variational(data=data2)
+    vbc2=opt.stan_variable("coef")
     meansigvb2=np.dot(B2, vbc2)  
     vprint(10, "Done 2nd Optimizing Stan model")
 
@@ -475,7 +481,7 @@ def baseline_pupil_model(tx,sy,event_onsets, fs=1000, lp1=2, lp2=0.2):
     event_onsets_ix=np.argmin(np.abs(np.tile(event_onsets, (sy.size,1)).T-tx), axis=1)
 
     # set up a single regressor
-    x1=np.zeros(sy.size, dtype=np.float)
+    x1=np.zeros(sy.size, dtype=float)
     x1[event_onsets_ix]=1
     kernel=pupil_kernel(4, fs=fs)
     x1=np.convolve(x1, kernel, mode="full")[0:x1.size]

pypillometry/baseline.py

@@ -11,8 +11,6 @@
 def nprange(ar):
     return (ar.min(),ar.max())
 
-
-import pystan
 import pickle
 from hashlib import md5
 
@@ -68,25 +66,6 @@ def trans_logistic_vec(x, a, b, inverse=False):
 
     return x
 
-def StanModel_cache(model_code, model_name=None, **kwargs):
-    """Use just as you would `stan`"""
-    code_hash = md5(model_code.encode('ascii')).hexdigest()
-    if model_name is None:
-        cache_fn = 'cached-model-{}.pkl'.format(code_hash)
-    else:
-        cache_fn = 'cached-{}-{}.pkl'.format(model_name, code_hash)
-    try:
-        sm = pickle.load(open(cache_fn, 'rb'))
-    except:
-        sm = pystan.StanModel(model_code=model_code)
-        with open(cache_fn, 'wb') as f:
-            pickle.dump(sm, f)
-    else:
-        print("Using cached StanModel")
-    return sm
-
-
-
 def plot_pupil_ipy(tx, sy, event_onsets=None, overlays=None, overlay_labels=None, 
                    blinks=None, interpolated=None,
                    figsize=(16,8), xlab="ms", nsteps=100):

pypillometry/convenience.py

@@ -101,7 +101,7 @@ def generate_pupil_data(event_onsets, fs=1000, pad=5000, baseline_lowpass=0.2,
 
     ### real events regressor
     ## scaling
-    event_ix=(np.array(event_onsets)/1000.*fs).astype(np.int)
+    event_ix=(np.array(event_onsets)/1000.*fs).astype(int)
     #a, b = (myclip_a - my_mean) / my_std, (myclip_b - my_mean) / my_std
     delta_weights=stats.truncnorm.rvs(-1/response_fluct_sd,np.inf, loc=1, scale=response_fluct_sd, size=event_ix.size)
     x1=np.zeros_like(sy)

pypillometry/fakedata.py

@@ -139,7 +139,7 @@ def detect_blinks_zero(sy, min_duration, blink_val=0):
     -------
     np.array (nblinks x 2) containing the indices of the start/end of the blinks
     """
-    x=np.r_[0, np.diff((sy==blink_val).astype(np.int))]
+    x=np.r_[0, np.diff((sy==blink_val).astype(int))]
     starts=np.where(x==1)[0]
     ends=np.where(x==-1)[0]-1
     if sy[0]==blink_val: ## first value missing?
@@ -206,7 +206,7 @@ def blink_onsets_mahot(sy, blinks, smooth_winsize, vel_onset, vel_offset, margin
         onset=max(0, onset-margin[0]) # avoid overflow to the left
 
         # find start of "reversal period" and move forward until it drops back
-        offset_ix=np.argmin(np.abs(((offsets-endl<0)*np.iinfo(np.int).max)+(offsets-endl)))
+        offset_ix=np.argmin(np.abs(((offsets-endl<0)*np.iinfo(int).max)+(offsets-endl)))
         while(offset_ix<(len(offsets)-1) and offsets[offset_ix+1]-1==offsets[offset_ix]):
             offset_ix+=1        
         offset=offsets[offset_ix]

pypillometry/preproc.py

@@ -102,7 +102,7 @@ def pupil_kernel(duration=4000, fs=1000, npar=10.1, tmax=930.0):
         sampled version of h(t) over the interval [0,`duration`] with sampling rate `fs`
     """
     n=int(duration/1000.*fs)
-    t = np.linspace(0,duration, n, dtype = np.float)
+    t = np.linspace(0,duration, n, dtype = float)
     h=pupil_kernel_t(t,npar,tmax)
     #h = t**(npar) * np.exp(-npar*t / tmax)   #Erlang gamma function Hoek & Levelt (1993)
     #hmax=np.exp(-npar)*tmax**npar ## theoretical maximum
@@ -159,7 +159,7 @@ def pupil_build_design_matrix(tx,event_onsets,fs,npar,tmax,max_duration="estimat
     h=pupil_kernel(duration=max_duration, fs=fs, npar=npar, tmax=tmax) ## pupil kernel
 
     # event-onsets for each event
-    x1 = np.zeros((event_onsets.size, tx.size), dtype=np.float) # onsets
+    x1 = np.zeros((event_onsets.size, tx.size), dtype=float) # onsets
 
     # event-onsets as indices of the txd array
     evon_ix=np.argmin(np.abs(np.tile(event_onsets, (tx.size,1)).T-tx), axis=1)
@@ -174,7 +174,7 @@ def pupil_build_design_matrix(tx,event_onsets,fs,npar,tmax,max_duration="estimat
     # h=pupil_kernel(duration=max_duration, fs=fs, npar=npar, tmax=tmax) ## pupil kernel
     # 
     # # event-onsets for each event
-    # x1 = np.zeros((event_onsets.size, tx.size), dtype=np.float) # onsets
+    # x1 = np.zeros((event_onsets.size, tx.size), dtype=float) # onsets
     # 
     # # event-onsets as indices of the txd array
     # evon_ix=np.argmin(np.abs(np.tile(event_onsets, (tx.size,1)).T-tx), axis=1)

pypillometry/pupil.py

@@ -165,15 +165,15 @@ def __init__(self,
             sometimes, when the eyetracker loses signal, no entry in the EDF is made; 
             when this option is True, such entries will be made and the signal set to 0 there
         """
-        self.sy=np.array(pupil, dtype=np.float)
+        self.sy=np.array(pupil, dtype=float)
         if sampling_rate is None and time is None:
             raise ValueError("you have to specify either sampling_rate or time-vector (or both)")
 
         if time is None:
             maxT=len(self)/sampling_rate*1000.
             self.tx=np.linspace(0,maxT, num=len(self))
         else:
-            self.tx=np.array(time, dtype=np.float)
+            self.tx=np.array(time, dtype=float)
 
         if sampling_rate is None:
             self.fs=np.round(1000./np.median(np.diff(self.tx)))
@@ -216,9 +216,9 @@ def __init__(self,
                 self.sy=nsy
 
         if event_onsets is None:
-            self.event_onsets=np.array([], dtype=np.float)
+            self.event_onsets=np.array([], dtype=float)
         else:
-            self.event_onsets=np.array(event_onsets, dtype=np.float)
+            self.event_onsets=np.array(event_onsets, dtype=float)
 
         # check whether onsets are in range
         for onset in self.event_onsets:
@@ -253,12 +253,12 @@ def __init__(self,
         self.response_estimated=False
 
         ## initialize blinks
-        self.blinks=np.empty((0,2), dtype=np.int)
-        self.blink_mask=np.zeros(len(self), dtype=np.int)
+        self.blinks=np.empty((0,2), dtype=int)
+        self.blink_mask=np.zeros(len(self), dtype=int)
 
         ## interpolated mask
-        self.interpolated_mask=np.zeros(len(self), dtype=np.int)
-        self.missing=np.zeros(len(self), dtype=np.int)
+        self.interpolated_mask=np.zeros(len(self), dtype=int)
+        self.missing=np.zeros(len(self), dtype=int)
         self.missing[self.sy==0]=1
 
         self.original=None
@@ -364,8 +364,8 @@ def sub_slice(self, start: float=-np.inf, end: float=np.inf, units: str="sec"):
         slic.event_onsets=slic.event_onsets[keepev]
         slic.event_labels=slic.event_labels[keepev]
         ## just remove all detected blinks (need to rerun `detect_blinks()`)
-        slic.blinks=np.empty((0,2), dtype=np.int)
-        slic.blink_mask=np.zeros(len(slic), dtype=np.int)
+        slic.blinks=np.empty((0,2), dtype=int)
+        slic.blink_mask=np.zeros(len(slic), dtype=int)
         return slic
 
     def summary(self) -> dict:
@@ -596,8 +596,8 @@ def _plot(self, plot_range, overlays, overlay_labels, units, interactive, highli
         overlays=(ov[startix:endix] for ov in overlays)
 
         if interactive:
-            blinks=np.empty((0,2), dtype=np.int)
-            interpolated=np.empty((0,2), dtype=np.int)
+            blinks=np.empty((0,2), dtype=int)
+            interpolated=np.empty((0,2), dtype=int)
             if highlight_blinks:
                 blinks=[]
                 for sblink,eblink in self.blinks:
@@ -966,7 +966,7 @@ def blinks_detect(self, min_duration:float=20, blink_val:float=0,
         blinks=helper_merge_blinks(blinks_vel, blinks_zero)
         obj.blinks=np.array([[on,off] for (on,off) in blinks if off-on>=min_duration_ix])
 
-        obj.blink_mask=np.zeros(self.sy.size, dtype=np.int)
+        obj.blink_mask=np.zeros(self.sy.size, dtype=int)
 
         for start,end in obj.blinks:
             obj.blink_mask[start:end]=1

pypillometry/pupildata.py

@@ -0,0 +1,60 @@
+// Stan model for pupil-baseline estimation
+//
+functions{
+    // asymmetric laplace function with the mu, sigma, tau parametrization
+    real my_skew_double_exponential_lpdf(real y, real mu, real sigma, real tau) {
+      return log(tau) + log1m(tau)
+        - log(sigma)
+        - 2 * ((y < mu) ? (1 - tau) * (mu - y) : tau * (y - mu)) / sigma;
+    }
+    
+    // zero-centered asymmetric laplace function with the mu, lambda, kappa parametrization
+    real skew_double_exponential2_lpdf(real y, real lam, real kappa) {
+      return log(lam) - log(kappa+1/kappa)
+        + ((y<0) ? (lam/kappa) : (-lam*kappa))*(y);
+    }
+}
+data{ 
+  int<lower=1> n; // number of timepoints in the signal
+  vector[n] sy;   // the pupil signal
+    
+  int<lower=1> ncol; // number of basis functions (columns in B)
+  matrix[n,ncol] B;  // spline basis functions
+    
+  int<lower=1> npeaks;          // number of lower peaks in the signal
+  int<lower=1> peakix[npeaks];  // index of the lower peaks in sy
+  vector<lower=0>[npeaks] lam_prominences; // lambda-converted prominence values
+    
+  real<lower=0> lam_sig;    // lambda for the signal where there is no peak
+  real<lower=0,upper=1> pa; // proportion of allowed distribution below 0
+}
+
+transformed data{
+    vector[n] lam;       // lambda at each timepoint
+    real<lower=0> kappa; // calculated kappa from pa
+    kappa=sqrt(pa)/sqrt(1-pa);
+    
+    lam=rep_vector(lam_sig, n); 
+    for(i in 1:npeaks){
+        lam[peakix[i]]=lam_prominences[i];
+    }
+}
+parameters {
+    vector[ncol] coef; // coefficients for the basis-functions
+}
+
+transformed parameters{
+    
+}
+
+model {
+    {
+    vector[n] d;
+        
+    coef ~ normal(0,5);
+    d=sy-(B*coef); // center at estimated baseline
+    for( i in 1:n ){
+        d[i] ~ skew_double_exponential2(lam[i], kappa);
+    }
+    }
+}

pypillometry/stan/baseline_model_asym_laplac.stan

@@ -1,6 +1,6 @@
 numpy
 scipy
-pystan
+cmdstanpy
 matplotlib
 pandas
 typing

requirements.txt

@@ -18,12 +18,19 @@
     long_description=long_description,
     long_description_content_type="text/markdown",
     url="https://github.com/ihrke/pypillometry",
-    packages=setuptools.find_packages(),
+    packages=setuptools.find_namespace_packages(),#where="pypillometry", exclude="tests"),
     classifiers=[
         "Programming Language :: Python :: 3",
         "License :: OSI Approved :: MIT License",
         "Operating System :: OS Independent",
     ],
+    #include_package_data=True,
+    package_dir={"pypillometry": "pypillometry"},
+    package_data={
+        "pypillometry": [],
+        "pypillometry.stan": ['*.stan'],
+        "pypillometry.tests": []
+        },
     install_requires=requirements,
-    python_requires='>=3.6',
+    python_requires='>=3.10',
 )