def plot_delta():     
    beta = 0.99
    N = 1000
    u = lambda c: sp.sqrt(c)
    W = sp.linspace(0,1,N)
    X, Y = sp.meshgrid(W,W)
    Wdiff = sp.transpose(X-Y)
    index = Wdiff <0
    Wdiff[index] = 0
    util_grid = u(Wdiff)
    util_grid[index] = -10**10
    
    Vprime = sp.zeros((N,1))
    delta = sp.ones(1)
    tol = 10**-9
    it = 0
    max_iter = 500
    
    while (delta[-1] >= tol) and (it < max_iter):
        V = Vprime
        it += 1;
        print(it)
        val = util_grid + beta*sp.transpose(V)
        Vprime = sp.amax(val, axis = 1)
        Vprime = Vprime.reshape((N,1))
        delta = sp.append(delta,sp.dot(sp.transpose(Vprime - V),Vprime-V))
        
    plt.figure()
    plt.plot(delta[1:])
    plt.ylabel(r'$\delta_k$')
    plt.xlabel('iteration')
    plt.savefig('convergence.pdf')

    def func(self, X, V):
        k = self.C.TFdata.k
        v1 = self.C.TFdata.v1
        w1 = self.C.TFdata.w1
        
        if k >=0:
            J_coords = self.F.sysfunc.J_coords
            w = sqrt(k)
        
            q = v1 - (1j/w)*matrixmultiply(self.F.sysfunc.J_coords,v1)
            p = w1 + (1j/w)*matrixmultiply(transpose(self.F.sysfunc.J_coords),w1)
            
            p /= linalg.norm(p)
            q /= linalg.norm(q)

            p = reshape(p,(p.shape[0],))
            q = reshape(q,(q.shape[0],))
            
            direc = conjugate(1/matrixmultiply(transpose(conjugate(p)),q))
            p = direc*p

            l1 = firstlyapunov(X, self.F.sysfunc, w, J_coords=J_coords, p=p, q=q)
            
            return array([l1])
        else:
            return array([1])

def sqcover(A,n):
    edge = sp.sqrt(A) # the length of an edge
    d = edge/n # the distance between two adjacent points
    r = d/2 # the "radius of "
    end = edge - r # end point
    base = sp.linspace(r, end, n)
    first_line = sp.transpose(sp.vstack((base, r*sp.ones(n))))
    increment = sp.transpose(sp.vstack((sp.zeros(n), d*sp.ones(n))))
    pts = first_line
    y_diff = increment
    for i in range(n-1):
        pts = sp.vstack((pts, first_line + y_diff))
        y_diff = y_diff + increment
    
    # Color matter
    colors = []
    for p in pts:
        cval = n*p[0] + p[1] # the x-coord has a higher weight
        cval = colormap.Spectral(cval/((n+1)*end)) # normalize by the max value that cval can take.
        colors.append(cval)

    colors = sp.array(colors)

    cover = (pts, r, colors)
    return cover

def trueFeatureStats(T, R, fMap, discountFactor, stateProp=1, MAT_LIMIT=1e8):
    """ Gather the statistics needed for LSTD,
    assuming infinite data (true probabilities).
    Option: if stateProp is  < 1, then only a proportion of all 
    states will be seen as starting state for transitions """
    dim = len(fMap)
    numStates = len(T)
    statMatrix = zeros((dim, dim))
    statResidual = zeros(dim)
    ss = range(numStates)
    repVersion = False
    
    if stateProp < 1:
        ss = random.sample(ss, int(numStates * stateProp))
    elif dim * numStates**2 < MAT_LIMIT:
        repVersion = True
    
    # two variants, depending on how large we can afford our matrices to become.        
    if repVersion:    
        tmp1 = tile(fMap, (numStates,1,1))
        tmp2 = transpose(tmp1, (2,1,0))
        tmp3 = tmp2 - discountFactor * tmp1            
        tmp4 = tile(T, (dim,1,1))
        tmp4 *= transpose(tmp1, (1,2,0))
        statMatrix = tensordot(tmp3, tmp4, axes=[[0,2], [1,2]]).T
        statResidual = dot(R, dot(fMap, T).T)
    else:
        for sto in ss:
            tmp = fMap - discountFactor * repmat(fMap[:, sto], numStates, 1).T
            tmp2 = fMap * repmat(T[:, sto], dim, 1)
            statMatrix += dot(tmp2, tmp.T)             
            statResidual += R[sto] * dot(fMap, T[:, sto])
    return statMatrix, statResidual

def Problem3Real():
    beta = 0.9
    N = 1000
    u = lambda c: sp.sqrt(c)
    W = sp.linspace(0,1,N)
    X, Y = sp.meshgrid(W,W)
    Wdiff = sp.transpose(X-Y)
    index = Wdiff <0
    Wdiff[index] = 0
    util_grid = u(Wdiff)
    util_grid[index] = -10**10
    
    Vprime = sp.zeros((N,1))
    psi = sp.zeros((N,1))
    delta = 1.0
    tol = 10**-9
    it = 0
    max_iter = 500
    
    while (delta >= tol) and (it < max_iter):
        V = Vprime
        it += 1;
        #print(it)
        val = util_grid + beta*sp.transpose(V)
        Vprime = sp.amax(val, axis = 1)
        Vprime = Vprime.reshape((N,1))
        psi_ind = sp.argmax(val,axis = 1)
        psi    = W[psi_ind]
        delta = sp.dot(sp.transpose(Vprime - V),Vprime-V)
    
    return psi

 def updatedata(self, A):
     if self.update:
         if self.corr:
             self.data.B = self.data.w*(linalg.norm(A,1)/linalg.norm(self.data.w,1))
             self.data.C = self.data.v*(linalg.norm(A,Inf)/linalg.norm(self.data.v,1))
         else:
             # Note: Problem when singular vectors switch smallest singular value (See NewLorenz).
             #       To overcome this, I have implemented a 1e-8 random nudge.
             try:
                 ALU = linalg.lu_factor(A)
                 BC = linalg.lu_solve(ALU, c_[linalg.lu_solve(ALU, self.data.B + 1e-8*self.data.Brand), \
                                      self.data.C + 1e-8*self.data.Crand], trans=1)
                 C = linalg.lu_solve(ALU, BC[:,-1*self.data.q:])
                 B = BC[:,0:self.data.p]
             except:
                 if self.C.verbosity >= 1:
                     print 'Warning: Problem updating border vectors.  Using svd...'
                 U, S, Vh = linalg.svd(A)
                 B = U[:,-1*self.data.p:]
                 C = num_transpose(Vh)[:,-1*self.data.q:]
         
             bmult = cmult = 1
             if matrixmultiply(transpose(self.data.B), B) < 0:
                 bmult = -1
             if matrixmultiply(transpose(self.data.C), C) < 0:
                 cmult = -1
             self.data.B = bmult*B*(linalg.norm(A,1)/linalg.norm(B))
             self.data.C = cmult*C*(linalg.norm(A,Inf)/linalg.norm(C))

def mlr(x,y,order):
    """Multiple linear regression fit of the columns of matrix x 
    (dependent variables) to constituent vector y (independent variables)
    
    order -     order of a smoothing polynomial, which can be included 
                in the set of independent variables. If order is
                not specified, no background will be included.
    b -         fit coeffs
    f -         fit result (m x 1 column vector)
    r -         residual   (m x 1 column vector)
    """
    
    if order > 0:
        s=scipy.ones((len(y),1))
        for j in range(order):
            s=scipy.concatenate((s,(scipy.arange(0,1+(1.0/(len(y)-1)),1.0/(len(y)-1))**j)[:,nA]),1)
        X=scipy.concatenate((x, s),1)
    else:
        X = x
    
    #calc fit b=fit coefficients
    b = scipy.dot(scipy.dot(scipy.linalg.pinv(scipy.dot(scipy.transpose(X),X)),scipy.transpose(X)),y)
    f = scipy.dot(X,b)
    r = y - f

    return b,f,r

def update_image(original_im, ci_red, ci_green, ci_blue):

    # diagnostics = dict()

    original_im = scipy.transpose(original_im)
    # diagnostics['original_im'] = original_im
    # diagnostics['ci_red'] = ci_red
    # diagnostics['ci_green'] = ci_green
    # diagnostics['ci_blue'] = ci_blue

    new_r = scipy.multiply(original_im[0], original_im[0] > ci_red)

    new_g = scipy.multiply(original_im[1], original_im[1] > ci_green)

    new_b = scipy.multiply(original_im[2], original_im[2] > ci_blue)

    new_im = (new_r, new_g, new_b)

    new_im = scipy.transpose(new_im)
    # diagnostics['new_im'] = new_im

    # with open('/Users/lages/Documents/sauceda/pictures_processed/diagnostics'
    #           '.p', 'wb') as f:
    #     pickle.dump(diagnostics, f)

    return new_im

def plot_disc_policy():
    #First compute policy function...==========================================
    N = 500
    w = sp.linspace(0,100,N)
    w = w.reshape(N,1)
    u = lambda c: sp.sqrt(c)
    util_vec = u(w)
    alpha = 0.5
    alpha_util = u(alpha*w)
    alpha_util_grid = sp.repeat(alpha_util,N,1)
    
    m = 20
    v = 200
    f = discretelognorm(w,m,v)
    
    VEprime = sp.zeros((N,1))
    VUprime    = sp.zeros((N,N))
    EVUprime = sp.zeros((N,1))
    psiprime = sp.ones((N,1))
    gamma = 0.1
    beta = 0.9
    
    m = 15
    tol = 10**-9
    delta = 1+tol
    it = 0
    while (delta >= tol):
        it += 1
        
        psi = psiprime.copy()
        arg1 = sp.repeat(sp.transpose(VEprime),N,0)
        arg2 = sp.repeat(EVUprime,N,1)
        arg = sp.array([arg2,arg1])
        psiprime = sp.argmax(arg,axis = 0) 
        
        for j in sp.arange(0,m):
            VE = VEprime.copy()
            VU = VUprime.copy()
            EVU = EVUprime.copy()
            VEprime = util_vec + beta*((1-gamma)*VE + gamma*EVU)
            arg1 = sp.repeat(sp.transpose(VE),N,0)*psiprime
            arg2 = sp.repeat(EVU,N,1)*(1-psiprime)
            arg = arg1+arg2
            VUprime = alpha_util_grid + beta*arg
            EVUprime = sp.dot(VUprime,f)  
    
        
    
        delta = sp.linalg.norm(psiprime -psi) 

    wr_ind = sp.argmax(sp.diff(psiprime), axis = 1)
    wr = w[wr_ind]
    print w[250],wr[250]
        
    #Then plot=================================================================
    plt.plot(w,psiprime[250,:]) 
    plt.ylim([-.5,1.5])      
    plt.xlabel(r'$w\prime$')
    plt.yticks([0,1])
    plt.savefig('disc_policy.pdf')

def __interpolateBetweenBinaryObjects(obj1, obj2, slices):
    """
    Takes two binary objects and puts slices slices in-between them, each of which
    contains a smooth binary transition between the objects.
    @note private inner function
    """
    if not obj1.shape == obj2.shape:
        raise AttributeError(
            "The two supplied objects have to be of the same shape, not {} and {}.".format(obj1.shape, obj2.shape)
        )

    # constant
    offset = 0.5  # must be a value smaller than the minimal distance possible
    temporal_dimension = 3

    # get all voxel position
    obj1_voxel = scipy.nonzero(obj1)
    obj2_voxel = scipy.nonzero(obj2)

    # get smallest pairwise distances between all object voxels
    distances = cdist(scipy.transpose(obj1_voxel), scipy.transpose(obj2_voxel))

    # keep for each True voxel of obj1 only the smallest distance to a True voxel in obj2
    min_distances = distances.min(1)

    # test if all seems to work
    if len(min_distances) != len(obj1_voxel[0]):
        raise Exception("Invalid number of minimal distances received.")

    # replace True voxels in obj1 with their respective distances to the True voxels in obj2
    thr_obj = obj1.copy()
    thr_obj = thr_obj.astype(scipy.float_)
    thr_obj[obj1_voxel] = min_distances
    thr_obj[obj1_voxel] += offset  # previous steps distances include zeros, therefore this is required

    # compute the step size for each slice that is added
    maximum = min_distances.max()
    step = maximum / float(slices + 1)
    threshold = maximum

    # control step: see if thr_obj really corresponds to obj1
    if not scipy.all(thr_obj.astype(scipy.bool_) == obj1.astype(scipy.bool_)):
        raise Exception("First created object does not correspond to obj1.")

    # assemble return volume
    return_volume = [thr_obj.astype(scipy.bool_)]  # corresponds to obj1
    for _ in range(slices):
        threshold -= step
        # remove all value higher than the threshold
        thr_obj[thr_obj > threshold] = 0
        # add binary volume to list /makes a copy)
        return_volume.append(thr_obj.astype(scipy.bool_))

    # add last slice (corresponds to es obj2 slice)
    thr_obj[thr_obj > offset] = 0
    return_volume.append(thr_obj.astype(scipy.bool_))

    # return binary scipy array
    return scipy.rollaxis(scipy.asarray(return_volume, dtype=scipy.bool_), 0, temporal_dimension + 1)

def baseline2(myarray):
    """Subtract average of the first and last bin from each bin
    """
    size = myarray.shape
    take_array = scipy.transpose(
        scipy.resize(scipy.transpose((myarray[:, 0] + myarray[:, size[1] - 1]) / 2), (size[1], size[0]))
    )
    return myarray - take_array

def calculate_ld(snps):
        #filter non binary snps
        snps_t = scipy.transpose(snps)
        snps_stand = scipy.transpose((snps_t - scipy.mean(snps, 1)) / scipy.std(snps, 1))
        r2_values =scipy.dot(snps_stand, scipy.transpose(snps_stand))
        r2_values *= (1.0 / snps.shape[1])
        r2_values **= 2
        return r2_values

def mls(p):

    data = ascii.read("datos.dat")  
    
    x=data["col1"]
    y=data["col2"]
    z=data["col3"]
    
    sig=30*(10**(-6))
    
    Y=sy.subtract(y,1)
    
    A=[]

    

    for m in x:
        
        pol=[]
        for i in range(p+1):
            pol.append(m**i)
        if m <0.4 or m>0.7:
            pol.append(0)
            A.append(pol) 
      
        else:
            pol.append(-1)
            A.append(pol)
            
    
    
    theta= sy.dot( sy.linalg.inv(sy.dot( sy.transpose(A),A )) , sy.dot(sy.transpose(A),Y) )

    modelo=[]
    
    for i in x:
        
        poli=1
        for s in range(p+1):
            poli+=(theta[s]*(i**s))    
        
        e=sy.random.normal(0,sig)
        if i <0.4 or i>0.7:
            modelo.append(poli)
      
        else:
            modelo.append(poli - theta[len(theta)-1])   
        
    
    chi2=0

    for h in range(len(x)):
        
        chi2+= ((y[h]-modelo[h]) / (sig) ) **2       
        
    return modelo, theta , len(x) ,sig ,chi2

def trim(image):  # 255 - white
    tr_image = transpose(image)
    start = 0
    while sum(tr_image[start]) == 255 * len(tr_image[start]):
        # condition on i is not needed, because of the balance between white and black
        start += 1
    finish = len(tr_image) - 1
    while sum(tr_image[finish]) == 255 * len(tr_image[finish]):
        finish -= 1
    return transpose(tr_image[start: finish + 1])

def Problem6Real():
    N = 500
    w = sp.linspace(0,100,N)
    w = w.reshape(N,1)
    u = lambda c: sp.sqrt(c)
    util_vec = u(w)
    alpha = 0.5
    alpha_util = u(alpha*w)
    alpha_util_grid = sp.repeat(alpha_util,N,1)
    
    m = 20
    v = 200
    f = discretelognorm(w,m,v)
    
    VEprime = sp.zeros((N,1))
    VUprime    = sp.zeros((N,N))
    EVUprime = sp.zeros((N,1))
    psiprime = sp.ones((N,1))
    gamma = 0.1
    beta = 0.9
    
    m = 15
    tol = 10**-9
    delta = 1+tol
    it = 0
    while (delta >= tol):
        it += 1
        
        psi = psiprime.copy()
        arg1 = sp.repeat(sp.transpose(VEprime),N,0)
        arg2 = sp.repeat(EVUprime,N,1)
        arg = sp.array([arg2,arg1])
        psiprime = sp.argmax(arg,axis = 0)    
        
        for j in sp.arange(0,m):
            VE = VEprime.copy()
            VU = VUprime.copy()
            EVU = EVUprime.copy()
            VEprime = util_vec + beta*((1-gamma)*VE + gamma*EVU)
            arg1 = sp.repeat(sp.transpose(VE),N,0)*psiprime
            arg2 = sp.repeat(EVU,N,1)*(1-psiprime)
            arg = arg1+arg2
            VUprime = alpha_util_grid + beta*arg
            EVUprime = sp.dot(VUprime,f)  
    
        
    
        delta = sp.linalg.norm(psiprime -psi)
        #print(delta)    
        
    wr_ind = sp.argmax(sp.diff(psiprime), axis = 1)
    wr = w[wr_ind]
    plt.plot(w,wr)
    plt.show()
    return wr

def shift_back_front(x,phi,neq,dx):
	n = len(x)/neq
	u = scipy.reshape(x,(neq,n))
	u_left = scipy.transpose([u[:,0]])
	u_right = scipy.transpose([u[:,-1]])
	u_ext = scipy.c_[u_left,u,u_right]
	dudx = 1./(2*dx)*(u_ext[:,2:]-u_ext[:,:-2])
	#u_ext = scipy.c_[u_left,u_left,u,u_right,u_right]
	#dudx = 1./(12*dx)*(-u_ext[:,4:]+8*u_ext[:,3:-1]-8*u_ext[:,1:-3]+u_ext[:,:-4])
	x_shifted=scipy.reshape(u+phi*dudx,(neq*n,))
	return x_shifted

def param(y):
    Y=sy.subtract(y,1)    
    
    A=[]
    
    for m in x:
        A.append([1,m,m**2,m**3,m**4,m**5]) 
    
    
    theta= sy.dot( sy.linalg.inv(sy.dot( sy.transpose(A),A )) , sy.dot(sy.transpose(A),Y) )
    return theta

 def shplot3d(self):
     fig = plt.figure()
     space = axs.Axes3D(fig)
     # plot the wall 
     points = self.xyzxyz()
     top = sp.transpose(points[0])
     space.scatter(top[0],top[1],top[2])
     # plot the pan bottom
     bottom = sp.transpose(points[1])
     space.scatter(bottom[0],bottom[1],bottom[2])
     plt.show()

    def __locate_newton(self, X, C):
        """x[0:self.dim] = (x,alpha)
           x[self.dim] = beta
           x[self.dim+1:2*self.dim] = p
        """
        J_coords = C.CorrFunc.jac(X[0:C.dim], C.coords)
        J_params = C.CorrFunc.jac(X[0:C.dim], C.params)

        return r_[C.CorrFunc(X[0:C.dim]) + X[C.dim]*X[C.dim+1:], \
                  matrixmultiply(transpose(J_coords),X[C.dim+1:]), \
                  matrixmultiply(transpose(X[C.dim+1:]),J_params), \
                  matrixmultiply(transpose(X[C.dim+1:]),X[C.dim+1:]) - 1]

def param(y):
    #Y=sy.subtract(y,1)
    Y=y    
    
    A=[]
    
    for m in x:
        A.append([1,m,m**2,m**3,m**4,m**5,m**6,m**7,m**8,m**9,m**10,m**11,m**12,m**13,m**14,m**15]) 
    
    
    theta= sy.dot( sy.linalg.inv(sy.dot( sy.transpose(A),A )) , sy.dot(sy.transpose(A),Y) )
    return theta