def blitz_inline(arr):
    """Prints the given 3D array by using blitz converters which
    provides a numpy-like syntax for accessing the numpy data.

    Notice the following:
      1. '\\n' to escape generating a newline in the C++ code.
      2. rows, cols = Narr[0], Narr[1].
      3. Array access using arr(i, j, k).

    """

    code = """
    int rows = Narr[0];
    int cols = Narr[1];
    int depth = Narr[2];
    for (int i=0; i < rows; i++)
    {
        for (int j=0; j < cols; j++)
        {
            printf("img[%3d][%3d]=", i, j);
            for (int k=0; k< depth; ++k)
            {
                printf(" %3d", arr(i, j, k));
            }
            printf("\\n");
        }
    }
    """

    weave.inline(code, ['arr'], type_converters=converters.blitz)

def pure_inline(arr):
    """Prints the given 3D array by accessing the raw numpy data and
    without using blitz converters.

    Notice the following:
      1. '\\n' to escape generating a newline in the C++ code.
      2. rows, cols = Narr[0], Narr[1].
      3. Array access using arr[(i*cols + j)*depth + k].

    """

    code = """
    int rows = Narr[0];
    int cols = Narr[1];
    int depth = Narr[2];
    for (int i=0; i < rows; i++)
    {
        for (int j=0; j < cols; j++)
        {
            printf("img[%3d][%3d]=", i, j);
            for (int k=0; k< depth; ++k)
            {
                printf(" %3d", arr[(i*cols + j)*depth + k]);
            }
            printf("\\n");
        }
    }
    """

    weave.inline(code, ['arr'])

def shift_sum(v1, shifts, bins):
    real_type = real_same_precision_as(v1)
    shifts = numpy.array(shifts, dtype=real_type)
    
    bins = numpy.array(bins, dtype=numpy.uint32)
    blen = len(bins) - 1
    v1 = numpy.array(v1.data, copy=False)
    slen = len(v1)

    if v1.dtype.name == 'complex64':
        code = point_chisq_code_single
    else:
        code = point_chisq_code_double
    
    n = int(len(shifts))
    
    # Create some output memory
    chisq =  numpy.zeros(n, dtype=real_type)
    
    inline(code, ['v1', 'n', 'chisq', 'slen', 'shifts', 'bins', 'blen'],
                    extra_compile_args=[WEAVE_FLAGS] + omp_flags,
                    libraries=omp_libs
          )
          
    return  chisq

def logistic_map(x0, r, T):
    """
    Returns a time series of length T using the logistic map
    x_(n+1) = r*x_n(1-x_n) at parameter r and using the initial condition x0.

    INPUT: x0 - Initial condition, 0 <= x0 <= 1
            r - Bifurcation parameter, 0 <= r <= 4
            T - length of the desired time series
    """
    #  Initialize the time series array
    timeSeries = np.empty(T)

    r = float(r)

    code = r"""
    int i;
    double xn;

    // Set initial condition
    timeSeries(0) = x0;

    for (i = 1; i < T; i++) {
        xn = timeSeries(i-1);
        timeSeries(i) = r * xn * (1 - xn);
    }
    """
    args = ['x0', 'r', 'T', 'timeSeries']
    weave.inline(code, arg_names=args, type_converters=weave.converters.blitz,
                 compiler='gcc', extra_compile_args=['-O3'])

    return timeSeries

def solve(n,m,t0,t1,dt,nu,f,verbose):
    """
    takes in the arguments needed to calculate the heat diffusion.
    includes verbose mode. 
    """
    if verbose:
        print "Calculation from {0} to {1} with dt={2}".format(t0,t1,dt)
    if not isinstance(f,np.ndarray):
        f1=f
        f=np.zeros((m,n))
        f.fill(f1)

    u=np.zeros((m,n))
    u_new=np.zeros((m,n))
    t=float(t0)
    expr = """
    int i;
    int j;
    while(t<t1){
    for(i=0;i<100;i++){
    for(j=0;j<50;j++){
    u_new(i,j)=u(i,j) + dt*(nu*u(i-1,j) + nu*u(i, j-1) - 4*nu*u(i, j) + nu*u(i,j+1) + nu*u(i+1, j) + f(i, j));
    }
  }
    u=u_new;
    t=t+dt;

  }"""
    if verbose:
        print "Uses weaves inline function to run the program in c for performance";
    weave.inline(expr, ['n','m','u','u_new','nu','dt','t','t1','f'], type_converters=weave.converters.blitz, compiler='gcc')
    if verbose:
        print "Returns the numpy array with the calculations"
    return u_new

def recent_moving_average(x, axis = 0):
    """
    Fast computation of recent moving average, where

        frac = 1/sqrt(t)
        a[t] = (1-frac)*a[t-1] + frac*x[t]
    """

    import weave  # ONLY WORKS IN PYTHON 2.X !!!
    if x.ndim!=2:
        y = recent_moving_average(x.reshape(x.shape[0], x.size//x.shape[0]), axis=0)
        return y.reshape(x.shape)

    assert x.ndim == 2 and axis == 0, 'Only implemented for a special case!'
    result = np.zeros(x.shape)
    code = """
    int n_samples = Nx[0];
    int n_dim = Nx[1];
    for (int i=0; i<n_dim; i++)
        result[i] = x[i];
    int ix=n_dim;
    for (int t=1; t<n_samples; t++){
        float frac = 1./sqrt(t+1);
        for (int i=0; i<n_dim; i++){
            result[ix] = (1-frac)*result[ix-n_dim] + frac*x[ix];
        }
        ix += 1;
    }
    """
    weave.inline(code, ['x', 'result'], compiler = 'gcc')
    return result

def heat_equation(t0, t1, dt, n, m, u, f, nu):

	# Used as a placeholder
	u_new = np.zeros((n, m))

	code = """

		int stop = t1/dt;

		// Timestep loop
		for(t0; t0 < stop; t0++)
		{
			int i, j;
			for(i = 1; i < n - 1; i++)
			{
				for(j = 1; j < m - 1; j++)
				{
					U_NEW2(i,j) = U2(i, j) + dt*(nu*U2(i-1,j) + nu*U2(i,j-1) - 4*nu*U2(i,j) + nu*U2(i,j+1) + nu*U2(i+1,j) + F2(i,j));
				}
			}

			// Update U2 with the values from this run.
			// This can probably be done with  pointer swap also
			for(i = 1; i < n - 1; i++)
			{
				for(j = 1; j < m - 1; j++)
				{
					U2(i, j) = U_NEW2(i, j);
				}
			}
		}
		"""

	inline(code, ['t0', 't1', 'dt', 'n', 'm', 'u', 'f', 'nu', 'u_new'])
	return u

def decode(llr):
    N = llr.size//2
    x = (llr[:N*2].reshape(-1,2,1)*output_map_soft).sum(1)
    msg = np.empty(N, np.uint8)
    weave.inline("""
    const int M = 128;
    int64_t cost[M*2], scores[M] = {/* zero-initialized */};
    uint8_t bt[N][M];
    for (int k=0; k<N; k++) {
        for (int i=0; i<M; i++) {
            cost[2*i+0] = scores[((i<<1) & 127) | 0] + x(k, i);
            cost[2*i+1] = scores[((i<<1) & 127) | 1] + x(k, i);
        }
        for (int i=0; i<M; i++) {
            int a = cost[2*i+0];
            int b = cost[2*i+1];
            bt[k][i] = (a<b) ? 1 : 0;
            scores[i] = (a<b) ? b : a;
        }
    }
    int i = (scores[0] < scores[1]) ? 1 : 0;
    for (int k=N-1; k>=0; k--) {
        int j = bt[k][i];
        msg(k) = i >> 6;
        i = ((i<<1)&127) + j;
    }
    """, ['N','x','msg'], type_converters=weave.converters.blitz)
    return msg

def aggregate(group_idx, a, func='sum', size=None, fill_value=0, order='C',
             dtype=None, axis=None, **kwargs):
    func = get_func(func, aliasing, optimized_funcs)
    if not isstr(func):
        raise NotImplementedError("generic functions not supported, in the weave implementation of aggregate")

    # Preparations for optimized processing
    group_idx, a, flat_size, ndim_idx, size = input_validation(group_idx, a,
                                                               size=size,
                                                               order=order,
                                                               axis=axis)
    dtype = check_dtype(dtype, func, a, len(group_idx))
    check_fill_value(fill_value, dtype)
    nans = func.startswith('nan')

    if nans:
        flat_size += 1

    if func in ('sum', 'any', 'len', 'anynan', 'nansum', 'nanlen'):
        ret = np.zeros(flat_size, dtype=dtype)
    elif func in ('prod', 'all', 'allnan', 'nanprod'):
        ret = np.ones(flat_size, dtype=dtype)
    else:
        ret = np.full(flat_size, fill_value, dtype=dtype)

    # In case we should get some ugly fortran arrays, convert them
    inline_vars = dict(group_idx=np.ascontiguousarray(group_idx), a=np.ascontiguousarray(a),
                       ret=ret, fill_value=fill_value)
    # TODO: Have this fixed by proper raveling
    if func in ('std', 'var', 'nanstd', 'nanvar'):
        counter = np.zeros_like(ret, dtype=int)
        inline_vars['means'] = np.zeros_like(ret)
        inline_vars['ddof'] = kwargs.pop('ddof', 0)
    elif func in ('mean', 'nanmean'):
        counter = np.zeros_like(ret, dtype=int)
    else:
        # Using inverse logic, marking anyting touched with zero for later removal
        counter = np.ones_like(ret, dtype=bool)
    inline_vars['counter'] = counter

    if np.isscalar(a):
        func += 'scalar'
        inline_vars['a'] = a
    inline(c_funcs[func], inline_vars.keys(), local_dict=inline_vars, define_macros=c_macros, extra_compile_args=c_args)

    # Postprocessing
    if func in ('sum', 'any', 'anynan', 'nansum') and fill_value != 0:
        ret[counter] = fill_value
    elif func in ('prod', 'all', 'allnan', 'nanprod') and fill_value != 1:
        ret[counter] = fill_value

    if nans:
        # Restore the shifted return array
        ret = ret[1:]

    # Deal with ndimensional indexing
    if ndim_idx > 1:
        ret = ret.reshape(size, order=order)
    return ret

def Ramp_list1(result, start, end):
    code = """
           const int size = result.len();
           const double step = (end-start)/(size-1);
           for (int i = 0; i < size; i++)
               result[i] = start + step*i;
           """
    weave.inline(code, ["result","start", "end"], verbose=2)

def inner_inline_real(self, other):
    x = _np.array(self._data, copy=False)
    y = _np.array(other, copy=False)
    total = _np.array([0.], dtype=float64)
    N = len(self)
    inline(inner_code, ['x', 'y', 'total', 'N'], libraries=omp_libs,
           extra_compile_args=code_flags)
    return total[0]

    def embed_time_series_array(self, time_series_array, dimension, delay):
        """
        Return a :index:`delay embedding` of all time series.

        .. note::
           Only works for scalar time series!

        **Example:**

        >>> ts = Surrogates.SmallTestData().original_data
        >>> Surrogates.SmallTestData().embed_time_series_array(
        ...     time_series_array=ts, dimension=3, delay=2)[0,:6,:]
        array([[ 0.        ,  0.61464833,  1.14988147],
               [ 0.31244015,  0.89680225,  1.3660254 ],
               [ 0.61464833,  1.14988147,  1.53884177],
               [ 0.89680225,  1.3660254 ,  1.6636525 ],
               [ 1.14988147,  1.53884177,  1.73766672],
               [ 1.3660254 ,  1.6636525 ,  1.76007351]])

        :type time_series_array: 2D array [index, time]
        :arg time_series_array: The time series array to be normalized.
        :arg int dimension: The embedding dimension.
        :arg int delay: The embedding delay.
        :rtype: 3D array [index, time, dimension]
        :return: the embedded time series.
        """
        if self.silence_level <= 1:
            print "Embedding all time series in dimension", dimension, \
                  "and with lag", delay, "..."
        (N, n_time) = time_series_array.shape

        embedding = np.empty((N, n_time - (dimension - 1)*delay, dimension))

        code = r"""
        int i, j, k, max_delay, len_embedded, index;

        //  Calculate the maximum delay
        max_delay = (dimension - 1)*delay;
        //  Calculate the length of the embedded time series
        len_embedded = n_time - max_delay;

        for (i = 0; i < N; i++) {
            for (j = 0; j < dimension; j++) {
                index = j*delay;
                for (k = 0; k < len_embedded; k++) {
                    embedding(i,k,j) = time_series_array(i,index);
                    index++;
                }
            }
        }
        """
        args = ['N', 'n_time', 'dimension', 'delay', 'time_series_array',
                'embedding']
        weave.inline(code, arg_names=args,
                     type_converters=weave.converters.blitz, compiler='gcc',
                     extra_compile_args=['-O3'])

        return embedding

def batch_correlate_execute(self, y):
    num_vectors = self.num_vectors # pylint:disable=unused-variable
    size = self.size # pylint:disable=unused-variable
    x = numpy.array(self.x.data, copy=False) # pylint:disable=unused-variable
    z = numpy.array(self.z.data, copy=False) # pylint:disable=unused-variable
    y = numpy.array(y.data, copy=False)        
    inline(batch_correlator_code, ['x', 'y', 'z', 'size', 'num_vectors'],
                extra_compile_args=[WEAVE_FLAGS] + omp_flags,
                libraries=omp_libs)

def Ramp_numeric1(result,start,end):
    code = """
           const int size = Nresult[0];
           const double step = (end-start)/(size-1);
           double val = start;
           for (int i = 0; i < size; i++)
               *result++ = start + step*i;
           """
    weave.inline(code,['result','start','end'],compiler='gcc')

def correlate_simd(ht, st, qt):
    htilde = _np.array(ht.data, copy = False).view(dtype = float32)
    stilde = _np.array(st.data, copy = False).view(dtype = float32) # pylint:disable=unused-variable
    qtilde = _np.array(qt.data, copy = False).view(dtype = float32) # pylint:disable=unused-variable
    arrlen = len(htilde) # pylint:disable=unused-variable
    inline(corr_simd_code, ['htilde', 'stilde', 'qtilde', 'arrlen'],
           extra_compile_args = [WEAVE_FLAGS],
           #extra_compile_args = ['-mno-avx -mno-sse2 -mno-sse3 -mno-ssse3 -mno-sse4 -mno-sse4.1 -mno-sse4.2 -mno-sse4a -O2 -w'],
           #extra_compile_args = ['-msse3 -O3 -w'],
           support_code = corr_support, auto_downcast = 1)

def step_indices(group_idx):
    """ Get the edges of areas within group_idx, which are filled 
        with the same value
    """
    ilen = step_count(group_idx) + 1
    indices = np.empty(ilen, int)
    indices[0] = 0
    indices[-1] = group_idx.size
    inline(c_step_indices, ['group_idx', 'indices'], define_macros=c_macros, extra_compile_args=c_args)
    return indices

def fast_second_phase(invec, indices, N1, N2):
    """
    This is the second phase of the FFT decomposition that actually performs
    the pruning. It is an explicit calculation for the subset of points. Note
    that there seem to be some numerical accumulation issues at various values
    of N1 and N2.

    Parameters
    ----------
    invec :
        The result of the first phase FFT
    indices : array of ints
        The index locations to calculate the FFT
    N1 : int
        The length of the second phase "FFT"
    N2 : int
        The length of the first phase FFT

    Returns
    -------
    out : array of floats
    """
    invec = numpy.array(invec.data, copy=False)
    NI = len(indices) # pylint:disable=unused-variable
    N1=int(N1)
    N2=int(N2)
    out = numpy.zeros(len(indices), dtype=numpy.complex64)

    # Note, the next step if this needs to be faster is to invert the loops
    code = """
        float pi = 3.14159265359;
        for(int i=0; i<NI; i++){
            float sp, cp;
            std::complex<double> val= (0, 0);

            unsigned int k = indices[i];
            int N = N1*N2;
            float k2 = k % N2;

            float phase_inc = 2 * pi * float(k) / float(N);
            sincosf(phase_inc, &sp, &cp);
            std::complex<float> twiddle_inc = std::complex<float>(cp, sp);
            std::complex<float> twiddle = std::complex<float>(1, 0);

            for (float n1=0; n1<N1; n1+=1){
                val += twiddle * invec[int(k2 + N2*n1)];
                twiddle *= twiddle_inc;
            }
            out[i] = val;
        }
    """
    weave.inline(code, ['N1', 'N2', 'NI', 'indices', 'out', 'invec'],
                       extra_compile_args=[WEAVE_FLAGS])
    return out

def correlate_parallel(ht, st, qt):
    htilde = _np.array(ht.data, copy = False)
    stilde = _np.array(st.data, copy = False) # pylint:disable=unused-variable
    qtilde = _np.array(qt.data, copy = False) # pylint:disable=unused-variable
    arrlen = len(htilde) # pylint:disable=unused-variable
    segsize = default_segsize # pylint:disable=unused-variable
    inline(corr_parallel_code, ['htilde', 'stilde', 'qtilde', 'arrlen', 'segsize'],
           extra_compile_args = [WEAVE_FLAGS] + omp_flags, libraries = omp_libs,
           #extra_compile_args = ['-mno-avx -mno-sse2 -mno-sse3 -mno-ssse3 -mno-sse4 -mno-sse4.1 -mno-sse4.2 -mno-sse4a -O2 -w'],
           #extra_compile_args = ['-msse3 -O3 -w'],
           support_code = corr_support, auto_downcast = 1)

def in_place_mult(num, mat):
    """In-place multiplication of a matrix by a scalar.
    """
    nrow, ncol = mat.shape
    code = """
for(int i=0;i<nrow;++i)
    for(int j=0;j<ncol;++j)
	mat(i,j) *= num;

"""
    inline(code, ["num", "mat", "nrow", "ncol"], type_converters=converters.blitz)

def encode(y):
    output = np.empty(y.size*2, np.uint8)
    weave.inline("""
    int sh = 0, N = y.extent(blitz::firstDim);
    for (int i=0; i<N; i++) {
        sh = (sh>>1) ^ ((int)y(i) << 6);
        output(2*i+0) = output_map(0,sh);
        output(2*i+1) = output_map(1,sh);
    }
    """, ['y','output','output_map'], type_converters=weave.converters.blitz)
    return output