def dateRange(first, last):

    # Type check, float --> int
    if isinstance(first[0], float):
        temp = np.zeros_like(first, dtype='int')
        for i in xrange(temp.size):
            temp[i] = first[i]
        first = tuple(temp)

    if isinstance(last[0], float):
        temp = np.zeros_like(last, dtype='int')
        for i in xrange(temp.size):
            temp[i] = last[i]
        last = tuple(temp)

    # Initialize date dictionary
    dateList = {}

    # Populate dictionary
    first = dt.datetime(*first[:6])
    last = dt.datetime(*last[:6])
    n = (last + dt.timedelta(days=1) - first).days
    dateList['year'] = np.array([(first + dt.timedelta(days=i)).year for i in xrange(n)])
    dateList['month'] = np.array([(first + dt.timedelta(days=i)).month for i in xrange(n)])
    dateList['day'] = np.array([(first + dt.timedelta(days=i)).day for i in xrange(n)])

    return dateList

def viterbi_decode(score, transition_params):
  """Decode the highest scoring sequence of tags outside of TensorFlow.

  This should only be used at test time.

  Args:
    score: A [seq_len, num_tags] matrix of unary potentials.
    transition_params: A [num_tags, num_tags] matrix of binary potentials.

  Returns:
    viterbi: A [seq_len] list of integers containing the highest scoring tag
        indicies.
    viterbi_score: A float containing the score for the Viterbi sequence.
  """
  trellis = np.zeros_like(score)
  backpointers = np.zeros_like(score, dtype=np.int32)
  trellis[0] = score[0]

  for t in range(1, score.shape[0]):
    v = np.expand_dims(trellis[t - 1], 1) + transition_params
    trellis[t] = score[t] + np.max(v, 0)
    backpointers[t] = np.argmax(v, 0)

  viterbi = [np.argmax(trellis[-1])]
  for bp in reversed(backpointers[1:]):
    viterbi.append(bp[viterbi[-1]])
  viterbi.reverse()

  viterbi_score = np.max(trellis[-1])
  return viterbi, viterbi_score

def reg(psf_model, parms):
    """
    Regularization and derivative.
    """
    eps = parms.eps
    if (eps is None):
        return np.zeros_like(psf_model)

    psf_shape = psf_model.shape
    d = np.zeros_like(psf_model)
    r = np.zeros_like(psf_model)
    for i in range(psf_shape[0]):
        for j in range(psf_shape[1]): 
            if i > 0:
                r[i, j] += (psf_model[i, j] - psf_model[i - 1, j]) ** 2.
                d[i, j] += 2. * (psf_model[i, j] - psf_model[i - 1, j]) 
            if j > 0:
                r[i, j] += (psf_model[i, j] - psf_model[i, j - 1]) ** 2.
                d[i, j] += 2. * (psf_model[i, j] - psf_model[i, j - 1]) 
            if i < psf_shape[0] - 1:
                r[i, j] += (psf_model[i, j] - psf_model[i + 1, j]) ** 2.
                d[i, j] += 2. * (psf_model[i, j] - psf_model[i + 1, j]) 
            if j < psf_shape[1] - 1:
                r[i, j] += (psf_model[i, j] - psf_model[i, j + 1]) ** 2.
                d[i, j] += 2. * (psf_model[i, j] - psf_model[i, j + 1]) 
    r *= eps
    d *= eps
    return r, d

    def finalize(self):
        """Calculates the flux, inverse variance and resolution for this spectrum.

        Uses the accumulated data from all += operations so far but does not prevent
        further accumulation.  This is the expensive step in coaddition so we make
        it something that you have to call explicitly.  If you forget to do this,
        the flux,ivar,resolution attributes will be None.

        If the coadded resolution matrix is not invertible, a warning message is
        printed and the returned flux vector is zero (but ivar and resolution are
        still valid).
        """
        # Convert to a dense matrix if necessary.
        if scipy.sparse.issparse(self.Cinv):
            self.Cinv = self.Cinv.todense()
        # What pixels are we using?
        mask = (np.diag(self.Cinv) > 0)
        keep = np.arange(len(self.Cinv_f))[mask]
        keep_t = keep[:,np.newaxis]
        # Initialize the results to zero.
        self.flux = np.zeros_like(self.Cinv_f)
        self.ivar = np.zeros_like(self.Cinv_f)
        R = np.zeros_like(self.Cinv)
        # Calculate the deconvolved flux,ivar and resolution for ivar > 0 pixels.
        self.ivar[mask],R[keep_t,keep] = decorrelate(self.Cinv[keep_t,keep])
        try:
            R_it = scipy.linalg.inv(R[keep_t,keep].T)
            self.flux[mask] = R_it.dot(self.Cinv_f[mask])/self.ivar[mask]
        except np.linalg.linalg.LinAlgError:
            self.log.warning('resolution matrix is singular so no coadded fluxes available.')
        # Convert R from a dense matrix to a sparse one.
        self.resolution = desispec.resolution.Resolution(R)

def totalvalue(cash_ini,orderform,valueform):
    
    trades = pd.read_csv(orderform,header=None,sep=',')
    trades = trades.dropna(axis = 1, how='all')
    trades.columns = ['Year','Month','Day','Symbol','Order','Share']
    dateall = []
    for i in np.arange(len(trades.Year)):
        dateall.append(dt.datetime(trades['Year'][i],trades['Month'][i],trades['Day'][i],16))
    dateall = pd.to_datetime(dateall)
    trades=trades.drop(['Year','Month','Day'],axis=1)
    trades['Date']=dateall
    trades.set_index('Date',inplace=True)
    
    ls_symbols = []
    for symbol in trades.Symbol:
        if symbol not in ls_symbols:
            ls_symbols.append(symbol)
            
    startdate = dateall[0]
    enddate = dateall[-1]
    dt_timeofday = dt.timedelta(hours=16)
    ldt_timestamps = du.getNYSEdays(startdate,enddate+dt_timeofday,dt_timeofday)
    ls_keys = 'close'
    c_dataobj = da.DataAccess('Yahoo')
    price = c_dataobj.get_data(ldt_timestamps, ls_symbols, ls_keys)
    orders = price*np.NaN
    orders = orders.fillna(0)
    for i in np.arange(len(trades.index)):
        ind = trades.index[i]
        if trades.ix[i,'Order']=='Buy':
            orders.loc[ind,trades.ix[i,'Symbol']]+=trades.ix[i,'Share']
        else:
            orders.loc[ind,trades.ix[i,'Symbol']]+=-trades.ix[i,'Share']
    #    keys = ['price','orders']
    #    trading_table = pd.concat([ldf_data,orders],keys=keys,axis=1)
    cash = np.zeros(np.size(price[ls_symbols[0]]),dtype=np.float)
    cash[0] = cash_ini
    # updating the cash value
    for i in np.arange(len(orders.index)):
        if i == 0: 
            cash[i] = cash[i] - pd.Series.sum(price.ix[i,:]*orders.ix[i,:])
        else:
            cash[i] = cash[i-1] - pd.Series.sum(price.ix[i,:]*orders.ix[i,:])
    # updating ownership
    ownership = orders*np.NaN
    for i in np.arange(len(orders.index)):
        ownership.ix[i,:]=orders.ix[:i+1,:].sum(axis=0) 
        
    # updating total portofolio value
    value = np.zeros_like(cash)
    for i in np.arange(len(ownership.index)):
        value[i] = pd.Series.sum(price.ix[i,:]*ownership.ix[i,:]) 
    keys = ['price','orders','ownership']
    trading_table = pd.concat([price,orders,ownership],keys = keys, axis=1)
    trading_table[('value','CASH')]=cash
    trading_table[('value','STOCK')]=value
    total = np.zeros_like(cash)
    total = cash + value
    trading_table[('value','TOTAL')]=total
    trading_table[('value','TOTAL')].to_csv(valueform)

    def fit_deriv(cls, x, amplitude, x_0, width):
        """One dimensional Box model derivative with respect to parameters"""

        d_amplitude = cls.evaluate(x, 1, x_0, width)
        d_x_0 = np.zeros_like(x)
        d_width = np.zeros_like(x)
        return [d_amplitude, d_x_0, d_width]

def dataVtk_3dMatrix(points,bounds,vectors):
    """
    Function that turns a vtk output formated data to 3d field matrix data
    from [(x1,y1,z1),...,(xn,yn,zn)]
    to [[[[x1,y1,z1],[...],[x3,y1,z1]],[[x1,y2,z1],[...],[...]],[[x1,y3,z1],[...],[...]]]
    ,[[[x1,y1,z2],[...],[...]],[...],[...]] , [.........]]    
    -points => list of the coordinates of the poitns where the data is located.
    -bounds => bounds of the data.(Xmin,Xmax,Ymin,Ymax,Zmin,Zmax)
    -vectors => vector data of the field at the 'points'
    """
    #asign variables
    (xmin,xmax,ymin,ymax,zmin,zmax) = bounds
        
    #generate the output arrays
    grid3d = N.mgrid[zmin:zmax+1, ymin:ymax+1, xmin:xmax+1]
    pnts3d = N.zeros_like(grid3d[0],dtype= N.ndarray)
    vect3d = N.zeros_like(grid3d[0],dtype= N.ndarray)
    
    #loop and rearange
    for i in range(len(points)):
        x_t = points[i][0]
        y_t = points[i][1]
        z_t = points[i][2]
        pnts3d[z_t+zmax][y_t+ymax][x_t+xmax] = points[i]
        vect3d[z_t+zmax][y_t+ymax][x_t+xmax] = vectors[i]
        
    return {'points':pnts3d,'vectors':vect3d}

    def __init__(self, network, **kwargs):
        # due to the way that theano handles updates, we cannot update a
        # parameter twice during the same function call. so, instead of handling
        # everything in the updates for self.f_learn(...), we split the
        # parameter updates into two function calls. the first "prepares" the
        # parameters for the gradient computation by moving the entire model one
        # step according to the current velocity. then the second computes the
        # gradient at that new model position and performs the usual velocity
        # and parameter updates.

        self.params = network.params(**kwargs)
        self.momentum = kwargs.get('momentum', 0.5)

        # set up space for temporary variables used during learning.
        self._steps = []
        self._velocities = []
        for param in self.params:
            v = param.get_value()
            n = param.name
            self._steps.append(theano.shared(np.zeros_like(v), name=n + '_step'))
            self._velocities.append(theano.shared(np.zeros_like(v), name=n + '_vel'))

        # step 1. move to the position in parameter space where we want to
        # compute our gradient.
        prepare = []
        for param, step, velocity in zip(self.params, self._steps, self._velocities):
            prepare.append((step, self.momentum * velocity))
            prepare.append((param, param + step))

        logging.info('compiling NAG adjustment function')
        self.f_prepare = theano.function([], [], updates=prepare)

        super(NAG, self).__init__(network, **kwargs)

def get_jk_coulomb(mol, dm, hermi=1, coulomb_allow='SSSS',
                   opt_llll=None, opt_ssll=None, opt_ssss=None):
    if coulomb_allow.upper() == 'LLLL':
        logger.info(mol, 'Coulomb integral: (LL|LL)')
        j1, k1 = _call_veff_llll(mol, dm, hermi, opt_llll)
        n2c = j1.shape[1]
        vj = numpy.zeros_like(dm)
        vk = numpy.zeros_like(dm)
        vj[...,:n2c,:n2c] = j1
        vk[...,:n2c,:n2c] = k1
    elif coulomb_allow.upper() == 'SSLL' \
      or coulomb_allow.upper() == 'LLSS':
        logger.info(mol, 'Coulomb integral: (LL|LL) + (SS|LL)')
        vj, vk = _call_veff_ssll(mol, dm, hermi, opt_ssll)
        j1, k1 = _call_veff_llll(mol, dm, hermi, opt_llll)
        n2c = j1.shape[1]
        vj[...,:n2c,:n2c] += j1
        vk[...,:n2c,:n2c] += k1
    else: # coulomb_allow == 'SSSS'
        logger.info(mol, 'Coulomb integral: (LL|LL) + (SS|LL) + (SS|SS)')
        vj, vk = _call_veff_ssll(mol, dm, hermi, opt_ssll)
        j1, k1 = _call_veff_llll(mol, dm, hermi, opt_llll)
        n2c = j1.shape[1]
        vj[...,:n2c,:n2c] += j1
        vk[...,:n2c,:n2c] += k1
        j1, k1 = _call_veff_ssss(mol, dm, hermi, opt_ssss)
        vj[...,n2c:,n2c:] += j1
        vk[...,n2c:,n2c:] += k1
    return vj, vk

def sor(A, b):
    sol = []
    
    n = len(A)
    D = np.zeros_like(A)
    L = np.zeros_like(A)
    
    for i in range(0,n):
        D[i][i] = A[i][i];
        
    for i in range(0,n):
        for j in range(0,i):
            L[i][j] = -A[i][j];
    
    omega = omegafind(A,D)   
    Q = D/omega -L
    Tj = np.linalg.inv(Q).dot(Q-A)
    c = np.linalg.inv(Q).dot(b)
    x = np.zeros_like(b)
    

    for itr in range(ITERATION_LIMIT):
        x=Tj.dot(x) + c;

    sol = x
    
    return list(sol)

 def filter_frames(self, data):
     data = data[0]
     lp = gaussian_filter(data, 100)
     hp = data - lp # poormans background subtraction
     hp -= np.min(hp)
     sh = hp.shape
     print "here"
     hp = hp.astype('uint32')
     hp = flex.int(hp)
     print "here now"
     
     mask = flex.bool(np.ones_like(hp).astype('bool'))
     print "here now"
     result1 = flex.bool(np.zeros_like(hp).astype('bool'))
     spots = np.zeros_like(hp).astype('bool')
     print "here now"
     
     for i in range(3, self.parameters['spotsize'], 5):
         print "here now"
         algorithm = DispersionThreshold(sh, (i, i), 1, 1, 0, -1)
         print "here now"
         print type(hp), type(mask), type(result1)
         thing = algorithm(hp, mask, result1)
         print "here now"
         spots = spots + result1.as_numpy_array()
     return [data, spots*data]

def initialize_adam(parameters):
    '''
    初始化v和s,他们都是字典类型的向量,都包含了以下字段
    -key:'dW1','db1',...'dWL','dbL'
    -values:与对应的梯度/参数相同维度的值为零的numpy矩阵
    
    :param parameters: -包含了以下参数的字典变量
            parameters['W'+str(l)] = W1
            parameters['b'+str(l)] = bl
    :return: 
    v - 包含梯度的指数加权平均值,字段如下:
        v['dW'+str(l)] = ...
        v['db'+str(l)] = ...
    s - 包含平方梯度的指数加权平均值,字段如下:
        s['dW'+str(l)] = ...
        s['db'+str(l)] = ...
    '''

    L = len(parameters)//2
    v= {}
    s = {}

    for l in range(L):
        v['dW'+str(l+1)] = np.zeros_like(parameters['W'+str(l+1)])
        v['db'+str(l+1)] = np.zeros_like(parameters['b'+str(l+1)])

        s['dW'+str(l+1)] = np.zeros_like(parameters['W'+str(l+1)])
        s['db'+str(l+1)] = np.zeros_like(parameters['b'+str(l+1)])

    return(v,s)

def compute_normals(im_pos, n_offset=3):
    """
    Converts an XYZ image to a Normal Image
    --Input--
    im_pos : ndarray (NxMx3)
        Image with x/y/z values for each pixel
    n_offset : int
        Smoothness factor for calculating the gradient
    --Output--
    normals : ndarray (NxMx3)
        Image with normal vectors for each pixel
    """
    gradients_x = np.zeros_like(im_pos)
    gradients_y = np.zeros_like(im_pos)
    for i in range(3):
        gradients_x[:, :, i], gradients_y[:, :, i] = np.gradient(im_pos[:, :, i], n_offset)

    gradients_x /= np.sqrt(np.sum(gradients_x**2, -1))[:, :, None]
    gradients_y /= np.sqrt(np.sum(gradients_y**2, -1))[:, :, None]

    normals = np.cross(gradients_x.reshape([-1, 3]),
                       gradients_y.reshape([-1, 3])).reshape(im_pos.shape)
    normals /= np.sqrt(np.sum(normals**2, -1))[:, :, None]
    normals = np.nan_to_num(normals)

    return normals

    def calculate_feed_tonnage_constrain(self,schedule,opening = None,closing = None):
        if opening is None:
            opening = np.zeros(self.ndp,dtype=np.int)
        else:
            assert(schedule.shape[0] == opening.shape[0])
            
        if closing is None:
            closing = np.zeros(self.ndp,dtype=np.int)
            closing[:] = self.nperiods - 1
        else:
            assert(schedule.shape[0] == closing.shape[0])


        production_period = self.calculate_feed_tonnage(schedule,opening,closing)

        #calculate the deviation from feed targets
        #logger.debug("minimum_feed_production=%f",self.minimum_feed_production)
        #logger.debug("maximum_feed_production=%f",self.maximum_feed_production)

        minp = np.zeros_like(production_period)
        indices = np.where(production_period < self.minimum_feed_production)[0]
        if len(indices) > 0:
            minp[indices] = self.minimum_feed_production - production_period[indices]

        maxp = np.zeros_like(production_period)
        indices = np.where(production_period > self.maximum_feed_production)[0]
        if len(indices) > 0:
            maxp[indices] = production_period[indices] - self.maximum_feed_production

            
        return tuple(maxp) + tuple(minp)

	def predictionToPosition(self,pi, dim = 64):
		pirescale = np.expand_dims(pi, axis=1)
		pirescale = np.append(pirescale, np.zeros_like(pirescale), axis=1)
		positions = np.zeros_like(pirescale)
		positions[:,0] = pirescale[:,0] // dim
		positions[:,1] = pirescale[:,0] % dim
		return positions

def fdtd(input_grid, steps):
    grid = input_grid.copy()
    old_grid = np.zeros_like(input_grid)
    previous_grid = np.zeros_like(input_grid)

    l_x = grid.shape[0]
    l_y = grid.shape[1]

    for i in range(steps):
        np.copyto(previous_grid, old_grid)
        np.copyto(old_grid, grid)

        for x in range(l_x):
            for y in range(l_y):
                grid[x,y] = 0.0
                if 0 < x+1 < l_x:
                    grid[x,y] += old_grid[x+1,y]
                if 0 < x-1 < l_x:
                    grid[x,y] += old_grid[x-1,y]
                if 0 < y+1 < l_y:
                    grid[x,y] += old_grid[x,y+1]
                if 0 < y-1 < l_y:
                    grid[x,y] += old_grid[x,y-1]

                grid[x,y] /= 2.0
                grid[x,y] -= previous_grid[x,y]

    return grid

    def _transform_dense(self, X):
        non_zero = (X != 0.0)
        X_nz = X[non_zero]

        X_step = np.zeros_like(X)
        X_step[non_zero] = np.sqrt(X_nz * self.sample_interval_)

        X_new = [X_step]

        log_step_nz = self.sample_interval_ * np.log(X_nz)
        step_nz = 2 * X_nz * self.sample_interval_

        for j in range(1, self.sample_steps):
            factor_nz = np.sqrt(step_nz /
                                np.cosh(np.pi * j * self.sample_interval_))

            X_step = np.zeros_like(X)
            X_step[non_zero] = factor_nz * np.cos(j * log_step_nz)
            X_new.append(X_step)

            X_step = np.zeros_like(X)
            X_step[non_zero] = factor_nz * np.sin(j * log_step_nz)
            X_new.append(X_step)

        return np.hstack(X_new)

def curvesSimilar(t1, y1, t2, y2, tol):
    """
    This function returns True if the two given curves are similar enough within tol. Otherwise returns False.

    t1: time/domain of standard curve we assume to be correct
    y1: values of standard curve, usually either temperature in (K) or log of a mol fraction
    t2: time/domain of test curve
    y2: values of test curve, usually either temperature in (K) or log of a mol fraction

    The test curve is first synchronized to the standard curve using geatNearestTime function. We then calculate the value of
    abs((y1-y2')/y1), giving us a normalized difference for every point. If the average value of these differences is less
    than tol, we say the curves are similar.

    We choose this criteria because it is compatible with step functions we expect to see in ignition systems.
    """
    # Make synchornized version of t2,y2 called t2sync,y2sync.
    t2sync=numpy.zeros_like(t1)
    y2sync=numpy.zeros_like(t1)
    for i, timepoint1 in enumerate(t1):
        time_index = findNearest(t2, timepoint1)
        t2sync[i]=t2[time_index]
        y2sync[i]=y2[time_index]

    # Get R^2 value equivalent:
    normalizedError=abs((y1-y2sync)/y1)
    normalizedError=sum(normalizedError)/len(y1)

    if normalizedError > tol:
        return False
    else:
        return True

def num_sdss_rand_both_catalogs(hemi, grid):
    d_rand = load_sdss_rand_both_catalogs(hemi)
    n_noisy = grid.num_from_radecz(d_rand['ra'],d_rand['dec'], d_rand['z'])
    # get the distance-from-observer 3d grid
    d_obs = grid.distance_from_observer()
    d_obs_max = d_obs[n_noisy>0].max()
    d_obs_1d = np.linspace(0, d_obs_max+1., 100)
    n_1d = np.zeros_like(d_obs_1d)
    delta_d_obs = d_obs_1d[1]-d_obs_1d[0]
    for i,this_d_obs in enumerate(d_obs_1d):
        wh=np.where((np.abs(d_obs-this_d_obs)<(0.5*delta_d_obs)) & (n_noisy>0))
        if len(wh[0])==0: continue
        n_1d[i] = np.median(n_noisy[wh])
    # now interpolate n_1d onto 3d grid
    from scipy import interpolate
    f = interpolate.interp1d(d_obs_1d, n_1d)
    n_median = np.zeros_like(n_noisy)
    wh_ok_interp = np.where((d_obs>np.min(d_obs_1d))&(d_obs<np.max(d_obs_1d)))    
    n_median[wh_ok_interp] = f(d_obs[wh_ok_interp])

    weight = np.zeros_like(n_median)
    weight[n_noisy>12]=1.
    n_median *= weight
    #pl.figure(1); pl.clf(); pl.imshow(n_noisy[:,:,128], vmin=0,vmax=n_median.max()); pl.colorbar()
    #pl.figure(2); pl.clf(); pl.imshow(n_median[:,:,128], vmin=0,vmax=n_median.max()); pl.colorbar()
    #ipdb.set_trace()
    return n_median, weight

 def pop(self, index=-1):
     if not isinstance(index, int):
         msg = "{0} indices must be integers, not {1}"
         raise TypeError(msg.format(self.__class__.__name__,
                                    index.__class__.__name__))
     if index < 0:
         index += self._len
     if index < 0:
         raise IndexError
     if self.storage == "numpy":
         ret_data = datacopy(self._data[index])
         ret_lengths = None
         if self._arity > 1:
             ret_lengths = _get_lenarray_empty(ret_data.shape)
         ret = self._element(
                 _mode="from_numpy",
                 data_store=ret_data,
                 len_data_store=ret_lengths,
             )
         self._data[index:self._len-1] = self._data[index+1:self._len]
         try:
             self._data[self._len-1] = np.zeros_like(self._data[self._len-1])
         except ValueError: # numpy bug
             for field in self._data.dtype.fields:
                 self._data[self._len-1][field] = np.zeros_like(self._data[self._len-1][field])
         self._del_child_lengths(index)
     elif self._elementary:
         ret = self._data[:self._len][index]
         self._data.__delitem__(index)
     else:
         ret = self._children[:self._len][index].copy()
         self._children.__delitem__(index)
         self._data.__delitem__(index)
     self._len -= 1
     return ret