def intercept(self, y, u):
     if self.aspherics is not None:
         return Interface.intercept(self, y, u) # expensive iterative
     # replace the newton-raphson with the analytic solution
     c, k = self.curvature, self.conic
     if c == 0:
         return -y[:, 2]/u[:, 2] # flat
     if not k:
         uy = (u*y).sum(1)
         uu = 1.
         yy = np.square(y).sum(1)
     else:
         k = np.array([(1, 1, 1 + k)])
         uy = (u*y*k).sum(1)
         uu = (np.square(u)*k).sum(1)
         yy = (np.square(y)*k).sum(1)
     d = c*uy - u[:, 2]
     e = c*uu
     f = c*yy - 2*y[:, 2]
     g = np.sqrt(np.square(d) - e*f)
     if self.alternate_intersection:
         g *= -1
     #g *= np.sign(u[:, 2])
     s = -(d + g)/e
     return s

def predictions(weather_turnstile):


    features_df = pandas.DataFrame({'Hour': weather_turnstile['Hour'], 
                                    'rain': weather_turnstile['rain'],
                                    'meantempi': weather_turnstile['meantempi'],
                                    'meanwindspdi': weather_turnstile['meanwindspdi'],
                                    'precipi': weather_turnstile['precipi'],
                                    'HourSquared': np.square(weather_turnstile['Hour']),
                                    'meantempiSquared': np.square(weather_turnstile['meantempi']),
                                    'precipiSquared': np.square(weather_turnstile['precipi'])})
    label = weather_turnstile['ENTRIESn_hourly']

    # Adds y-intercept to model
    features_df = sm.add_constant(features_df)

    # add dummy variables of turnstile units to features
    dummy_units = pandas.get_dummies(weather_turnstile['UNIT'], prefix='unit')
    features_df = features_df.join(dummy_units)
    model = sm.OLS(label,features_df)

    results = model.fit()

    prediction = results.predict(features_df)
    return prediction

    def compute_distances_no_loops(self, X):
        """
        Compute the distance between each test point in X and each training point
        in self.X_train using no explicit loops.

        Input / Output: Same as compute_distances_two_loops
        """
        num_test = X.shape[0]
        num_train = self.X_train.shape[0]
        dists = np.zeros((num_test, num_train))
        #########################################################################
        # TODO:                                                                 #
        # Compute the l2 distance between all test points and all training      #
        # points without using any explicit loops, and store the result in      #
        # dists.                                                                #
        #                                                                       #
        # You should implement this function using only basic array operations; #
        # in particular you should not use functions from scipy.                #
        #                                                                       #
        # HINT: Try to formulate the l2 distance using matrix multiplication    #
        #       and two broadcast sums.                                         #
        #########################################################################
        X_train_square = np.square(self.X_train).sum(axis=1)
        test_square = np.square(X).sum(axis=1)
        M = -2 * np.dot(X, self.X_train.T)
        dists = np.sqrt(M + X_train_square + np.matrix(test_square).T)
        #########################################################################
        #                         END OF YOUR CODE                              #
        #########################################################################
        return dists

    def test_meanSquaredDisplacement(self):
        from getMeanSquareDisplacement import getMeanSquareDisplacement
        numTradingDays = 4*getNumTradingDaysPerYear()
        growthRate     = 0.5
        logVolatility  = 2.0**-8
        numCols = 15000
        
        p = makeFakeDailyPrices(growthRate,logVolatility,numTradingDays,numCols)
        
        
        Ex21 = np.mean(getMeanSquareDisplacement(np.log(p),10),1)
        t = np.arange(10)
        volTerm = t*np.square(logVolatility)
        driftTerm = np.square(t*np.log(1.0+growthRate)/getNumTradingDaysPerYear())
        Ex20 = volTerm + driftTerm

        error = np.sqrt(np.mean(np.square((Ex21[1:] - Ex20[1:])/Ex20[1:])))

        self.assertLess(error,0.002)

        if 0:    
            import matplotlib.pyplot as plt
            print 'MSE =',np.round(100*error,3),'%'
            plt.plot(t,Ex21,'ro ')
            plt.plot(t,Ex20,'g:')
            plt.show()

    def time_std(self):
        if hasattr(self, '_time_std'):
            return self._time_std
        if self.savedir is not None:
            try:
                with open(join(self.savedir, 'time_std.pkl'),
                          'rb') as f:
                    time_std = pickle.load(f)
            except IOError:
                pass
            else:
                # Same protocol as the averages. Make sure the
                # std is a single 4D (zyxc) array and if not just
                # re-calculate the time std.
                if isinstance(time_std, np.ndarray):
                    self._time_std = time_std
                    return self._time_std

        sums = np.zeros(self.frame_shape)
        sums_squares = np.zeros(self.frame_shape)
        counts = np.zeros(self.frame_shape)
        for frame in it.chain.from_iterable(self):
            sums += np.nan_to_num(frame)
            sums_squares += np.square(np.nan_to_num(frame))
            counts[np.isfinite(frame)] += 1
        means = old_div(sums, counts)
        mean_of_squares = old_div(sums_squares, counts)
        std = np.sqrt(mean_of_squares-np.square(means))
        if self.savedir is not None and not self._read_only:
            with open(join(self.savedir, 'time_std.pkl'), 'wb') as f:
                pickle.dump(std, f, pickle.HIGHEST_PROTOCOL)
        self._time_std = std
        return self._time_std

    def setup(self):
        self.add_parameter(FloatParameter('audio-brightness', 1.0))
        self.add_parameter(FloatParameter('audio-stripe-width', 100.0))
        self.add_parameter(FloatParameter('audio-speed', 0.0))
        self.add_parameter(FloatParameter('speed', 0.01))
        self.add_parameter(FloatParameter('angle-speed', 0.1))
        self.add_parameter(FloatParameter('stripe-width', 20))
        self.add_parameter(FloatParameter('center-orbit-distance', 200))
        self.add_parameter(FloatParameter('center-orbit-speed', 0.1))
        self.add_parameter(FloatParameter('hue-step', 0.1))
        self.add_parameter(IntParameter('posterization', 8))
        self.add_parameter(StringParameter('color-gradient', "[(0,0,1), (0,0,1), (0,1,1), (0,1,1), (0,0,1)]"))
        self.add_parameter(FloatParameter('stripe-x-center', 0.5))
        self.add_parameter(FloatParameter('stripe-y-center', 0.5))
        self.hue_inner = random.random() + 100
        self._center_rotation = random.random()
        self.stripe_angle = random.random()

        cx, cy = self.scene().center_point()
        self.locations = self.scene().get_all_pixel_locations()
        x,y = self.locations.T
        x -= cx
        y -= cy
        self.pixel_distances = np.sqrt(np.square(x) + np.square(y))
        self.pixel_angles = (math.pi + np.arctan2(y, x)) / (2 * math.pi)
        self.pixel_distances /= max(self.pixel_distances)

        super(StripeGradient, self).setup()

def outlierCleaner(predictions, ages, net_worths):
    """
        clean away the 10% of points that have the largest
        residual errors (difference between the prediction
        and the actual net worth)

        return a list of tuples named cleaned_data where 
        each tuple is of the form (age, net_worth, error)
    """
    
    cleaned_data = []
    
    # calculate threshold of 90%    
    residual_errors = net_worths - predictions 
    residual_errors_square = np.square(residual_errors)
    residual_errors_square.sort(axis = 0)
#    print residual_errors_square
    
    percentile_90_index = int(len(residual_errors_square) * .9)
    percentile_90_threshold = residual_errors_square[percentile_90_index - 1][0]
#    print "threshold", percentile_90_threshold
    
    cleaned_data_all = zip(ages[:, 0].tolist(), net_worths[:, 0].tolist(), residual_errors[:, 0].tolist())
    
#    count = 0
    
    for e in cleaned_data_all: 
        (age, net_worth, error) = e
        if np.square(error) <= percentile_90_threshold: 
#            print error, percentile_90_threshold
            cleaned_data.append(e)
#            count += 1
#    print count
    return cleaned_data

  def test_infer(self):
    kmeans = self.kmeans
    kmeans.fit(input_fn=self.input_fn(), relative_tolerance=1e-4)
    clusters = kmeans.clusters()

    # Make a small test set
    num_points = 10
    points, true_assignments, true_offsets = make_random_points(clusters,
                                                                num_points)
    # Test predict
    assignments = kmeans.predict(input_fn=self.input_fn(
        batch_size=num_points, points=points))
    self.assertAllEqual(assignments, true_assignments)

    # Test score
    score = kmeans.score(
        input_fn=lambda: (constant_op.constant(points), None), steps=1)
    self.assertNear(score, np.sum(true_offsets), 0.01 * score)

    # Test transform
    transform = kmeans.transform(
        input_fn=lambda: (constant_op.constant(points), None))
    true_transform = np.maximum(
        0,
        np.sum(np.square(points), axis=1, keepdims=True) - 2 * np.dot(
            points, np.transpose(clusters)) +
        np.transpose(np.sum(np.square(clusters), axis=1, keepdims=True)))
    self.assertAllClose(transform, true_transform, rtol=0.05, atol=10)

def sphDist(ra1, dec1, ra2, dec2):
    """Calculate distance on the surface of a unit sphere.

    Input and Output are in radians.

    Notes
    -----
    Uses the Haversine formula to preserve accuracy at small angles.

    Law of cosines approach doesn't work well for the typically very small
    differences that we're looking at here.
    """
    # Haversine
    dra = ra1-ra2
    ddec = dec1-dec2
    a = np.square(np.sin(ddec/2)) + \
        np.cos(dec1)*np.cos(dec2)*np.square(np.sin(dra/2))
    dist = 2 * np.arcsin(np.sqrt(a))

    # This is what the law of cosines would look like
#    dist = np.arccos(np.sin(dec1)*np.sin(dec2) + np.cos(dec1)*np.cos(dec2)*np.cos(ra1 - ra2))

    # Could use afwCoord.angularSeparation()
    #  but (a) that hasn't been made accessible through the Python interface
    #  and (b) I'm not sure that it would be faster than the numpy interface.
    #    dist = afwCoord.angularSeparation(ra1-ra2, dec1-dec2, np.cos(dec1), np.cos(dec2))

    return dist

    def fitness(self, recordings):
        """
        Calculates the sum squared difference between each spike in the
        signal and the closest spike in the reference spike train, plus the
        vice-versa case

        `analysis` -- The analysis object containing all recordings and
                      analysis of them [analysis.AnalysedRecordings]
        """
        spikes = recordings.get_analysed_signal().spikes()
        inner = spikes[numpy.where(
            (spikes >= (self.time_start + self.time_buffer)) &
            (spikes <= (self.time_stop - self.time_buffer)))]
        # If no spikes were generated create a dummy spike that is guaranteed
        # to be further away from a reference spike than any within the time
        # window
        if len(spikes) == 0:
            spike_t = self.time_stop + self.time_start
            spikes = neo.SpikeTrain([spike_t], spike_t, units=spike_t.units)
        fitness = 0.0
        for spike in inner:
            fitness += float(numpy.square(self.ref_spikes - spike).min())
        for ref_spike in self.ref_inner:
            fitness += float(numpy.square(spikes - ref_spike).min())
        return fitness

 def K(self, X, X2=None,alpha=None,variance=None):
     """
     Computes the covariance matrix cov(X[i,:],X2[j,:]).
     
     Args:
         X: Matrix where each row is a point.
         X2: Matrix where each row is a point.
         alpha: It's the scaled alpha.
         Variance: Sigma hyperparameter.
         
     """
     if alpha is None:
         alpha=self.alpha
     if variance is None:
         variance=self.variance
     
     if X2 is None:
         X=X*alpha/self.scaleAlpha
         Xsq=np.sum(np.square(X), 1)
         r=-2.*np.dot(X, X.T) + (Xsq[:, None] + Xsq[None, :])
         r = np.clip(r, 0, np.inf)
         return variance*np.exp(-0.5*r)
     else:
         X=X*alpha/self.scaleAlpha
         X2=X2*alpha/self.scaleAlpha
         r=-2.*np.dot(X, X2.T) + (np.sum(np.square(X), 1)[:, None] + np.sum(np.square(X2), 1)[None, :])
         r = np.clip(r, 0, np.inf)
         return variance*np.exp(-0.5*r)

 def hit(self, ray):
     # assume sphere at origin, so translate ray:
     raypoint = ray.point - self.point
     p0 = raypoint[0]
     p1 = raypoint[1]
     p2 = raypoint[2]
     v0 = ray.vector[0]
     v1 = ray.vector[1]
     v2 = ray.vector[2]
     a = ((N.square(v0))/(N.square(self.A))) + ((N.square(v1))/(N.square(self.B))) + ((N.square(v2))/(N.square(self.C)))
     b = ((2*p0*v0)/(N.square(self.A))) + ((2*p1*v1)/(N.square(self.B))) + ((2*p2*v2)/(N.square(self.C)))
     c = ((N.square(p0))/(N.square(self.A))) + ((N.square(p1))/(N.square(self.B))) + ((N.square(p2))/(N.square(self.C))) - 1
     disc = b*b - 4*a*c
     if disc > 0.0:
         t = (-b-N.sqrt(disc))/(2*a)
         if t > EPSILON:
             p = ray.pointAt(t)
             n = normalize(self.normalAt(p))
             return (t, p, n, self)
         t = (-b+N.sqrt(disc))/(2*a)
         if t > EPSILON:
             p = ray.pointAt(t)
             n = normalize(self.normalAt(p))
             return (t, p, n, self)
     return (None, None, None, None)

   def analyse(self, a):
      global motion_detected, motion_timestamp, motion_array, motion_array_mask
      # calcuate length of motion vectors of mpeg macro blocks
      a = np.sqrt(
          np.square(a['x'].astype(np.float)) +
          np.square(a['y'].astype(np.float))
          ).clip(0, 255).astype(np.uint8)
      a = a * motion_array_mask
      # If there're more than 'sensitivity' vectors with a magnitude greater
      # than 'threshold', then say we've detected motion
      th = ((a > motion_threshold).sum() > motion_sensitivity)
      now = time.time()
      # motion logic, trigger on motion and stop after 2 seconds of inactivity
      if th:
         motion_timestamp = now

      if motion_detected:
          if (now - motion_timestamp) >= video_postseconds:
               motion_detected = False
      else:
        if th:
             motion_detected = True
        if debug:
                idx = a > motion_threshold
                a[idx] = 255
                motion_array = a

def energy(x, y, z):
    ex = np.sqrt(np.sum(np.square(np.subtract(x,mean(x)))))
    ey = np.sqrt(np.sum(np.square(np.subtract(y,mean(y)))))
    ez = np.sqrt(np.sum(np.square(np.subtract(z,mean(z)))))
    
    e = (1/(3 * len(x))) * (ex + ey + ez)
    return e

	def test_quarticspike(self):
		rr = np.square(self.X) + np.square(self.Y)
		r = np.sqrt(rr)
		res = blowup.quartic_spike(r)
		npt.assert_allclose(res[0,0],0.)
		npt.assert_allclose(res[0,self.N//2], 0.)
		npt.assert_allclose(res[self.N//2, self.N//2],1.)

    def __iter__(self):
        
        MAX_X,MAX_Y = self.dimensions
        MIN_V, MAX_V = self.velocity
        
        wt_min = 0.
        
        if self.init_stationary:

            x, y, x_waypoint, y_waypoint, velocity, wt = \
                init_random_waypoint(self.nr_nodes, MAX_X, MAX_Y, MIN_V, MAX_V, wt_min, 
                             (self.wt_max if self.wt_max is not None else 0.))

        else:

            NODES = np.arange(self.nr_nodes)
            print NODES
            x = U(0, MAX_X, NODES)
            y = U(0, MAX_Y, NODES)
            x_waypoint = U(0, MAX_X, NODES)
            y_waypoint = U(0, MAX_Y, NODES)
            wt = np.zeros(self.nr_nodes)
            velocity = U(MIN_V, MAX_V, NODES)

        theta = np.arctan2(y_waypoint - y, x_waypoint - x)
        costheta = np.cos(theta)
        sintheta = np.sin(theta)
        
        while True:
            # update node position
            x += velocity * costheta
            y += velocity * sintheta
            # calculate distance to waypoint
            d = np.sqrt(np.square(y_waypoint-y) + np.square(x_waypoint-x))
            # update info for arrived nodes
            arrived = np.where(np.logical_and(d<=velocity, wt<=0.))[0]
            
            # step back for nodes that surpassed waypoint
            x[arrived] = x_waypoint[arrived]
            y[arrived] = y_waypoint[arrived]
            
            if self.wt_max:
                velocity[arrived] = 0.
                wt[arrived] = U(0, self.wt_max, arrived)
                # update info for paused nodes
                wt[np.where(velocity==0.)[0]] -= 1.
                # update info for moving nodes
                arrived = np.where(np.logical_and(velocity==0., wt<0.))[0]
            
            if arrived.size > 0:
                x_waypoint[arrived] = U(0, MAX_X, arrived)
                y_waypoint[arrived] = U(0, MAX_Y, arrived)
                velocity[arrived] = U(MIN_V, MAX_V, arrived)
                theta[arrived] = np.arctan2(y_waypoint[arrived] - y[arrived], x_waypoint[arrived] - x[arrived])
                costheta[arrived] = np.cos(theta[arrived])
                sintheta[arrived] = np.sin(theta[arrived])
            
            self.velocity = velocity
            self.wt = wt
            yield np.dstack((x,y))[0]

def genEmpCov_kernel(sigma, width, sample_set, knownMean = True):
    timesteps = sample_set.__len__()
#    print timesteps
    mean_tile = 0
    K_sum = 0
    if knownMean != True:
        for j in range(int(max(0,timesteps-width)),timesteps):            
            K =  np.exp(-np.square(timesteps-j-1)/sigma)
            samplesPerStep = sample_set[j].shape[1]
            mean_tile = mean_tile + K* sample_set[j]
            K_sum = K_sum + K
            
        mean_tile =  np.sum(mean_tile, axis = 1)/samplesPerStep
        mean_tile = np.tile(mean_tile, (samplesPerStep,1)).T
    K_sum = 0
    S = 0
#    print 'timesteps and width is %d, %d, %d'%(timesteps,width, max(0,timesteps- width))
   
    for j in range(int(max(0,timesteps-width)),timesteps):
        K = np.exp(-np.square(timesteps-j-1)/sigma)
#        print 'j = ',j, 'K = ', K
        samplesPerStep = sample_set[j].shape[1]
        S = S + K*np.dot(sample_set[j]- mean_tile, (sample_set[j] -  mean_tile).T)/samplesPerStep
        K_sum = K_sum + K
    S = S/K_sum
    return S

def pick_triplets_impl(q_in, q_out):
  more = True
  while more:
      deq = q_in.get()
      if deq is None:
        more = False
      else:
        embeddings, emb_start_idx, nrof_images, alpha = deq
        print('running', emb_start_idx, nrof_images, os.getpid())
        for j in xrange(1,nrof_images):
            a_idx = emb_start_idx + j - 1
            neg_dists_sqr = np.sum(np.square(embeddings[a_idx] - embeddings), 1)
            for pair in xrange(j, nrof_images): # For every possible positive pair.
                p_idx = emb_start_idx + pair
                pos_dist_sqr = np.sum(np.square(embeddings[a_idx]-embeddings[p_idx]))
                neg_dists_sqr[emb_start_idx:emb_start_idx+nrof_images] = np.NaN
                all_neg = np.where(np.logical_and(neg_dists_sqr-pos_dist_sqr<alpha, pos_dist_sqr<neg_dists_sqr))[0]  # FaceNet selection
                #all_neg = np.where(neg_dists_sqr-pos_dist_sqr<alpha)[0] # VGG Face selecction
                nrof_random_negs = all_neg.shape[0]
                if nrof_random_negs>0:
                    rnd_idx = np.random.randint(nrof_random_negs)
                    n_idx = all_neg[rnd_idx]
                    #triplets.append( (a_idx, p_idx, n_idx) )
                    q_out.put( (a_idx, p_idx, n_idx) )
        #emb_start_idx += nrof_images
  print('exit',os.getpid())

    def update(self):
        rho = self.opt_config.rho
        epsilon = self.opt_config.epsilon
        lr = self.opt_config.lr
        clip = self.opt_config.clip

        all_norm = 0.0
        for param_name in self.apollo_net.active_param_names():
            param = self.apollo_net.params[param_name]
            grad = param.diff
            all_norm += np.sum(np.square(grad))
        all_norm = np.sqrt(all_norm)

        for param_name in self.apollo_net.active_param_names():
            param = self.apollo_net.params[param_name]
            grad = param.diff

            if all_norm > clip:
                grad = clip * grad / all_norm

            if param_name in self.sq_grads:
                self.sq_grads[param_name] = (1 - rho) * np.square(grad) + rho * self.sq_grads[param_name]
                rms_update = np.sqrt(self.sq_updates[param_name] + epsilon)
                rms_grad = np.sqrt(self.sq_grads[param_name] + epsilon)
                update = -rms_update / rms_grad * grad

                self.sq_updates[param_name] = (1 - rho) * np.square(update) + rho * self.sq_updates[param_name]
            else:
                self.sq_grads[param_name] = (1 - rho) * np.square(grad)
                update = np.sqrt(epsilon) / np.sqrt(epsilon + self.sq_grads[param_name]) * grad
                self.sq_updates[param_name] = (1 - rho) * np.square(update)

            param.data[...] += lr * update
            param.diff[...] = 0

def ratio_err(top,bottom,top_low,top_high,bottom_low,bottom_high):
    #uses simple propagation of errors (partial derivatives)
    #note it returns errorbars, not interval

    #-make sure input is numpy arrays-
    top = np.array(top)
    top_low = np.array(top_low)
    top_high = np.array(top_high)
    bottom = np.array(bottom)
    bottom_low = np.array(bottom_low)
    bottom_high = np.array(bottom_high)

    #-calculate errorbars-
    top_errlow = np.subtract(top,top_low)
    top_errhigh = np.subtract(top_high,top)
    bottom_errlow = np.subtract(bottom,bottom_low)
    bottom_errhigh = np.subtract(bottom_high,bottom)

    #-calculate ratio_low-
    ratio_low  = np.sqrt( np.square(np.divide(top_errlow,bottom)) + np.square( np.multiply(np.divide(top,np.square(bottom)),bottom_errlow)) )
    #-calculate ratio_high-
    ratio_high = np.sqrt( np.square(np.divide(top_errhigh,bottom)) + np.square( np.multiply(np.divide(top,np.square(bottom)),bottom_errhigh)) )
#    ratio_high = ((top_errhigh/bottom)**2.0 + (top/(bottom**2.0))*bottom_errhigh)**2.0)**0.5

    # return two vectors, err_low and err_high
    return ratio_low,ratio_high