def normal_test_case():
    '''
    Runs a test case with simulated data from a normal distribution.
    '''
    obs, fa, dur = [], [], []
    for n in range(15):
        d, f, o = make_test_data(
            5, split=min(plt.rand()*50+120, 170),
            intercept=plt.rand()*50 + 225,
            slope1=1 + plt.randn()/0.75, slope2=plt.randn()/.75)
        obs.append(o+n)
        fa.append(f)
        dur.append(d)
        plt.plot(f, d, 'o', alpha=0.1)

    dur, fa, obs = (np.hstack(dur)[:, np.newaxis],
                    np.hstack(fa)[:, np.newaxis],
                    np.hstack(obs)[:, np.newaxis])

    dur_mean = dur.mean()
    dur_std = dur.std()
    dur = (dur-dur_mean)/dur_std

    m = normal_model(dur, fa, obs)
    trace = sample_model(m, 5000)
    predict(trace, 5, 2500, {'mean': dur_mean, 'std': dur_std})
    plt.figure()
    traceplot(trace, 2, 2500)
    return dur, fa, obs, (dur_mean, dur_std), trace

def gamma_test_case():
    '''
    Runs a test case with simulated data from a normal distribution.
    '''
    obs, fa, dur = [], [], []
    delta_angle = np.arange(180)
    for n in range(15):
        mode = piecewise_predictor(
            delta_angle,
            100 + plt.randn()*20,
            250 + plt.randn()*20,
            1 + plt.randn()/2.0,
            -1 + plt.randn()/2.0)
        a, b = np_gamma_params(mode, 10)
        for _ in range(10):
            d = gamma.rvs(a=a, scale=1.0/b)
            fa.append(delta_angle)
            dur.append(d)
            obs.append(d*0+n)
    dur, fa, obs = np.concatenate(dur), np.concatenate(fa), np.concatenate(obs)
    m = gamma_model(dur, fa, obs.astype(int))
    trace = sample_model(m, 5000)
    predict(trace, 5, 2500 )
    plt.figure()
    traceplot(trace, 2, 2500)
    return dur, fa, obs, trace

 def __init__(self, numstates, numclasses):
     self.numstates = numstates
     self.numclasses = numclasses
     self.numparams = self.numclasses * self.numstates\
                      + self.numstates\
                      + self.numstates**2
     self.params = zeros(self.numparams, dtype=float)
     self.logEmissionProbs = \
                     self.params[:self.numclasses * self.numstates].reshape(
                                            self.numstates, self.numclasses)
     self.logInitProbs = self.params[self.numclasses * self.numstates:
                                     self.numclasses * self.numstates+
                                     self.numstates]
     self.logTransitionProbs = self.params[-self.numstates**2:].reshape(
                                             self.numstates, self.numstates)
     self.logEmissionProbs[:] = \
                 ones((self.numstates, self.numclasses),dtype=numpy.double)\
                                 /numpy.double(self.numclasses) +\
                         randn(self.numstates, self.numclasses)*0.001
     self.logEmissionProbs /= self.logEmissionProbs.sum(1)[:,newaxis]
     self.logEmissionProbs = log(self.logEmissionProbs)
     self.logInitProbs[:] = ones(self.numstates, dtype=float) \
                                                     / self.numstates+\
                         randn(self.numstates)*0.001
     self.logInitProbs /= self.logInitProbs.sum()
     self.logInitProbs[:] = log(self.logInitProbs)
     self.logTransitionProbs[:] = ones((self.numstates, self.numstates), 
                                              dtype=float)/self.numstates+\
                                  randn(self.numstates,self.numstates)*0.001
     self.logTransitionProbs /= self.logTransitionProbs.sum(1)[:,newaxis]
     self.logTransitionProbs[:] = log(self.logTransitionProbs)

def plot_F_and_pi(F, pi, causes, title=""):
    N, T, J = F.shape
    pl.figure(figsize=(0.5 * T, 1.0 * J))

    left = 2.0 / (T + 5.0)
    right = 1 - 0.05 / T
    bottom = 2.0 / (T + 5.0)
    top = 1 - 0.05 / T

    xmax = F.max()

    dj = (top - bottom) / J
    dt = (right - left) / T

    ax = {}
    for jj, j in enumerate(sorted(range(J), key=lambda j: pi[:, :, j].mean())):
        for t in range(T):
            pl.axes([left + t * dt, bottom + jj * dj, dt, dj])
            pl.plot(pl.randn(N), F[:, t, j], "b.", alpha=0.5, zorder=-100)
            pl.plot(0.5 * pl.randn(N), pi[:N, t, j], "g.", alpha=0.5, zorder=100)

            # pi[:,t,j].sort()
            # below = pi[:, t, j].mean() - pi[:,t,j][.025*N]
            # above = pi[:,t,j][.975*N] - pi[:, t, j].mean()
            # pl.errorbar([0], pi[:, t, j].mean(), [[below], [above]],
            #            fmt='gs', ms=10, mew=1, mec='white', linewidth=3, capsize=10,
            #            zorder=100)
            pl.text(
                -2.75,
                xmax * 0.9,
                "%.0f\n%.0f\n%.0f"
                % (
                    100 * F[:, t, j].mean(),
                    100 * pi[:, t, j].mean(),
                    100 * pi[:, t, j].mean() - 100 * F[:, t, j].mean(),
                ),
                va="top",
                ha="left",
            )
            pl.xticks([])
            pl.yticks([0.25, 0.5, 0.75])
            if jj == 0:
                pl.xlabel("%d\n%.0f" % (t + 1980, 100 * F[:, t, :].sum() / N))

            if t > 0:
                pl.yticks([])
            else:
                pl.ylabel(causes[j])

            pl.axis([-3, 3, 0, xmax])
    if title:
        pl.figtext(0.01, 0.99, title, va="top", ha="left")

def example_histogram_3():
    # first create a single histogram
    mu, sigma = 200, 25
    x = mu + sigma * plb.randn(10000)
        
    plb.figure(6)
    
    # finally: make a multiple-histogram of data-sets with different length
    x0 = mu + sigma * plb.randn(10000)
    x1 = mu + sigma * plb.randn(7000)
    x2 = mu + sigma * plb.randn(3000)
    
    n, bins, patches = plb.hist([x0, x1, x2], 10, histtype='bar')
    plb.show()

    def simulate(self, f_u, x0, tf):
        """
        Simulate the system.

        Parameters
        ----------
        f_u: The input function  f_u(t, x, i)
        x0: The initial state.
        tf: The final time.

        Return
        ------
        data : A StateSpaceDataArray object.

        """
        #pylint: disable=too-many-locals, no-member
        x0 = pl.matrix(x0)
        assert x0.shape[1] == 1
        t = 0
        x = x0
        dt = self.dt
        data = StateSpaceDataList([], [], [], [])
        i = 0
        n_x = self.A.shape[0]
        n_y = self.C.shape[0]
        assert pl.matrix(f_u(0, x0, 0)).shape[1] == 1
        assert pl.matrix(f_u(0, x0, 0)).shape[0] == n_y

        # take square root of noise cov to prepare for noise sim
        if pl.norm(self.Q) > 0:
            sqrtQ = scipy.linalg.sqrtm(self.Q)
        else:
            sqrtQ = self.Q

        if pl.norm(self.R) > 0:
            sqrtR = scipy.linalg.sqrtm(self.R)
        else:
            sqrtR = self.R

        # main simulation loop
        while t + dt < tf:
            u = f_u(t, x, i)
            v = sqrtR.dot(pl.randn(n_y, 1))
            y = self.measurement(x, u, v)
            data.append(t, x, y, u)
            w = sqrtQ.dot(pl.randn(n_x, 1))
            x = self.dynamics(x, u, w)
            t += dt
            i += 1
        return data.to_StateSpaceDataArray()

def genExampleAxisAligned():
    # first choose a label
    y = ''
    feats = {}
    if pylab.rand() < 0.5:   # negative example
        # from Nor([-1,0], 1)
        y = '-1'
        feats['x'] = pylab.randn() - 1
        feats['y'] = pylab.randn()
    else:
        # from Nor([+1,0], 1)
        y = '1'
        feats['x'] = pylab.randn() + 1
        feats['y'] = pylab.randn()
    return (y, feats)

def matplotlib_plot(fname):
    Xs, Ys = 6, 6
    fig = pylab.figure(figsize=(Xs, Ys), frameon=False)
    ax = fig.add_axes([0.0, 0.0, 1.0, 1.0], frameon=False)
    ax.set_xticks([])
    ax.set_yticks([])    
    N = 1000
    x, y = pylab.randn(N), pylab.randn(N)
    color = pylab.randn(N)
    size = abs(400*pylab.randn(N))
    p = pylab.scatter(x, y, c=color, s=size, alpha=0.75)
    pylab.xlim(-2.0, 2.0)
    pylab.ylim(-2.0, 2.0)
    pylab.savefig(fname, dpi=200)
    return Xs, Ys

    def estimate_vol_vol(self):

        # dummy time series: must be replaced by a reader from an external source
        rng = pd.date_range(start = '2010-01-01', end = '2011-01-01', freq='D')
        tmp_0 = 100000. + 10000.*pl.randn(len(rng))
        ts_volume = pd.Series(tmp_0, index=rng)
        tmp_1 = 0.02 + 0.002*pl.randn(len(rng))
        ts_volatility = pd.Series(tmp_1, index=rng)

        # estimation of the daily volume and of the daily volatility as a flat average of the previuos n_days_mav
        n_days_mav = 10
        period_start = pd.to_datetime('2010-03-01') + pd.DateOffset(days=-(n_days_mav+1))
        period_end = pd.to_datetime('2010-03-01') + pd.DateOffset(days=-1)
        self.volume_est = ts_volume[period_start:period_end].mean()
        self.volatility_est = ts_volatility[period_start:period_end].mean()

def test_speriodogram_2d():
    data = randn(1024,2)
    speriodogram(data)


    data = np.array([marple_data, marple_data]).reshape(64,2)
    speriodogram(data)

 def __init__(self,numin,numclasses):
     self.numin  = numin
     self.numclasses = numclasses
     self.params = 0.01 * randn(self.numin*self.numclasses+self.numclasses)
     self.scorefunc = logreg_score(self.numin,self.numclasses,self.params)
     self.scorefuncs = [scorefunc]
     Contrastive.__init__(self,normalizeacrosscliques=False)

 def __init__(self,numstates, numdims):
   self.numstates = numstates
   self.numdims = numdims
   self.numparams = self.numdims * self.numstates\
       + self.numdims**2 * self.numstates\
       + self.numstates\
       + self.numstates**2
   self.params = zeros(self.numparams, dtype=float)
   self.means = self.params[:self.numdims*self.numstates].reshape(
         self.numdims,self.numstates)
   self.covs = self.params[self.numdims * self.numstates:
                 self.numdims*self.numstates+
                 self.numdims**2 * self.numstates].reshape(
                     self.numdims, self.numdims, self.numstates)
   self.logInitProbs = self.params[self.numdims * self.numstates +
                         self.numdims**2 * self.numstates:
                         self.numdims * self.numstates +
                         self.numdims**2 * self.numstates+
                         self.numstates]
   self.logTransitionProbs = self.params[-self.numstates**2:].reshape(
                               self.numstates, self.numstates)
   self.means[:] = 0.1 * randn(self.numdims, self.numstates)
   for k in range(self.numstates):
     self.covs[:,:,k] = eye(self.numdims) * 0.1;
   self.logInitProbs[:] = log(ones(self.numstates, dtype=float) /
                             self.numstates)
   self.logTransitionProbs[:] = log(ones((self.numstates,
                                         self.numstates), dtype=float) /
                                   self.numstates)

 def __init__(self,numdims):
     self.numdims = numdims
     self.params = 0.001 * randn(numdims*2)
     self.wp = self.params[:numdims]
     self.wm = self.params[numdims:]
     self.b = 0.0
     self.t = 1.0

def synthbeats(duration, meanhr=60, stdhr=1, samplingfreq=250, sinfreq=None):
    #Minimaly based on the parameters from:
    #http://physionet.cps.unizar.es/physiotools/ecgsyn/Matlab/ecgsyn.m
    #If frequ exist it will be used to generate a sin instead of using rand
    #Inputs: duration in seconds
    #Returns: signal, peaks

    t = np.arange(duration * samplingfreq) / float(samplingfreq)
    signal = np.zeros(len(t))

    print(len(t))
    print(len(signal))

    if sinfreq == None:

        npeaks = 1.2 * (duration * meanhr / 60)
        # add 20% more beats for some cummulative error
        hr = pl.randn(npeaks) * stdhr + meanhr
        peaks = pl.cumsum(60. / hr) * samplingfreq
        peaks = peaks.astype('int')
        peaks = peaks[peaks < t[-1] * samplingfreq]

    else:
        hr = meanhr + sin(2 * pi * t * sinfreq) * float(stdhr)
        index = int(60. / hr[0] * samplingfreq)
        peaks = []
        while index < len(t):
            peaks += [index]
            index += int(60. / hr[index] * samplingfreq)

    signal[peaks] = 1.0

    return t, signal, peaks

def generate_fe(out_fname='data.csv'):
    """ Generate random data based on a fixed effects model

    This function generates data for all countries in all regions, based on the model::

        Y_r,c,t = beta * X_r,c,t

        beta = [10., -.5, .1, .1, -.1, 0., 0., 0., 0., 0.]

        X_r,c,t[0] = 1
        X_r,c,t[1] = t - 1990.
        X_r,c,t[k] ~ N(0, 1) for k >= 2
    """
    c4 = countries_by_region()

    a = 20.
    beta = [10., -.5, .1, .1, -.1, 0., 0., 0., 0., 0.]
    data = col_names()
    for t in time_range:
        for r in c4:
            for c in c4[r]:
                x = [1] + [t-1990.] + list(pl.randn(8))
                y = float(pl.dot(beta, x))
                se = 0.
                data.append([r, c, t, a, y, se] + list(x))

    write(data, out_fname)

def test_scatter_hist():
    import pylab

    X = pylab.randn(1000)
    Y = pylab.randn(1000)
    scatter_hist(X, Y)

    df = pd.DataFrame({"X": X, "Y": Y, "size": X, "color": Y})
    scatter_hist(df, hist_position="left")

    df = pd.DataFrame({"X": X})
    try:
        scatter_hist(df, hist_position="left")
        assert False
    except:
        assert True

def _smooth_demo():
    from numpy import linspace, sin, ones
    from pylab import subplot, plot, hold, axis, legend, title, show, randn

    t = linspace(-4, 4, 100)
    x = sin(t)
    xn = x + randn(len(t)) * 0.1
    y = smooth(x)

    ws = 31

    subplot(211)
    plot(ones(ws))

    windows = ["flat", "hanning", "hamming", "bartlett", "blackman"]

    hold(True)
    for w in windows[1:]:
        eval("plot(" + w + "(ws) )")

    axis([0, 30, 0, 1.1])

    legend(windows)
    title("The smoothing windows")
    subplot(212)
    plot(x)
    plot(xn)
    for w in windows:
        plot(smooth(xn, 10, w))
    l = ["original signal", "signal with noise"]
    l.extend(windows)

    legend(l)
    title("Smoothing a noisy signal")
    show()

def plot_axes(fig):
    # create some data to use for the plot
    dt = 0.001
    t = arange(0.0, 10.0, dt)
    r = exp(-t[:1000]/0.05)               # impulse response
    x = randn(len(t))
    s = convolve(x,r,mode=2)[:len(x)]*dt  # colored noise

    # the main axes is subplot(111) by default
    axes = fig.gca()
    axes.plot(t, s)
    axes.set_xlim((0, 1))
    axes.set_ylim((1.1*min(s), 2*max(s)))
    axes.set_xlabel('time (s)')
    axes.set_ylabel('current (nA)')
    axes.set_title('Gaussian colored noise')

    # this is an inset axes over the main axes
    a = fig.add_axes([.65, .6, .2, .2], axisbg='y')
    n, bins, patches = a.hist(s, 400, normed=1)
    a.set_title('Probability')
    a.set_xticks([])
    a.set_yticks([])

    # this is another inset axes over the main axes
    a = fig.add_axes([.2, .6, .2, .2], axisbg='y')
    a.plot(t[:len(r)], r)
    a.set_title('Impulse response')
    a.set_xlim((0, 0.2))
    a.set_xticks([])
    a.set_yticks([])

def groupConfidenceWeight(AllData):
  """
  weights answers by confidence for different groups
  """
  # everybodygetup
  subjects = range(len(AllData[1]['correct']))
  
  distribution = np.array(py.zeros([20,len(subjects)]))
  for i in subjects:
    newdist = getIndConf(AllData, i)
    print(len(newdist))
    distribution[:,i] = newdist
  m,n = py.shape(distribution)
  for i in xrange(m):
    for j in xrange(n):
      distribution[i,j] = distribution[i,j] + py.randn(1)*0.05
  
  print(distribution)
  
  fig = py.figure()
  ax20 = fig.add_subplot(111)
  for i in xrange(n):
    ax20.hist(distribution[:,i], bins=20, color='c', alpha=0.2,
              edgecolor='none')
  ax20.set_title('Weighted answers')
  ax20.set_xlabel('Distribution')
  ax20.set_ylabel('Counts')

 def _pvoc(self, X_hat, Phi_hat=None, R=None):
     """
     ::
       a phase vocoder - time-stretch
       inputs:
         X_hat - estimate of signal magnitude
         [Phi_hat] - estimate of signal phase
         [R] - resynthesis hop ratio
       output:
         updates self.X_hat with modified complex spectrum
     """
     N = self.nfft
     W = self.wfft
     H = self.nhop
     R = 1.0 if R is None else R
     dphi = (2*P.pi * H * P.arange(N/2+1)) / N
     print "Phase Vocoder Resynthesis...", N, W, H, R
     A = P.angle(self.STFT) if Phi_hat is None else Phi_hat
     phs = A[:,0]
     self.X_hat = []
     n_cols = X_hat.shape[1]
     t = 0
     while P.floor(t) < n_cols:
         tf = t - P.floor(t)            
         idx = P.arange(2)+int(P.floor(t))
         idx[1] = n_cols-1 if t >= n_cols-1 else idx[1]
         Xh = X_hat[:,idx]
         Xh = (1-tf)*Xh[:,0] + tf*Xh[:,1]
         self.X_hat.append(Xh*P.exp( 1j * phs))
         U = A[:,idx[1]] - A[:,idx[0]] - dphi
         U = U - P.np.round(U/(2*P.pi))*2*P.pi
         phs += (U + dphi)
         t += P.randn()*P.sqrt(PVOC_VAR*R) + R # 10% variance
     self.X_hat = P.np.array(self.X_hat).T