def adjust(p, method='bonferroni'):
    '''
    Usage:
        pvalues.adjust(p, method='bonferroni')
        returns p-values adjusted using one of several methods.
        p ... numeric vector of p-values (possibly with NAs)
        method ... correction method
            bonferroni: Use Bonferroni method to control the family-wise error rate strictly.
            BH: Use method of Benjamini and Hochberg to control the false discovery rate.
            fdr: same as 'BH'.
    '''
    try:
        p = np.array(p).astype(float)
        n = len(p)
    except ValueError:
        print 'Error: input p-values contain invalid string elements.'
        quit()
    # remove "n.a." values from input vector
    p0 = p
    not_na = ~np.isnan(p)
    p = p[not_na]
    lp = len(p)
    if lp <= 1:
        return p0
    if method == 'bonferroni':
        p0[not_na] = np.fmin(1., n*p)
    elif method == 'BH' or method == 'fdr':
        i = np.arange(lp+1)[:0:-1]
        o = np.argsort(p)[::-1]
        ro = np.argsort(o)
        p0[not_na] = np.fmin(1., np.minimum.accumulate(float(n) / i.astype(float) * p[o]))[ro]
 
    return p0

def compute_akt(pi3k_value, time_value, initial_values, mfs):
	"""Rules-
		If pi3k is high and time is high akt is high
		If pi3k is low or time is low then akt is low"""
	a1_1 = mfs[0][1][initial_values[0] == pi3k_value]
	a1_2 = mfs[2][1][initial_values[2] == time_value]

	if( a1_1.size == 0):
		a1_1 = mfs[0][1][ find_closest(initial_values[0], pi3k_value)]
	if( a1_2.size == 0):
		a1_2 = mfs[2][1][ find_closest(initial_values[2], time_value)]

	a1 = min(a1_1, a1_2)
	c1 = np.fmin( np.linspace(a1, a1, 100), mfs[1][1])

	a2_1 = mfs[0][0][initial_values[0] == pi3k_value]
	a2_2 = mfs[2][0][initial_values[2] == time_value]

	if( a2_1.size == 0):
		a2_1 = mfs[0][0][ find_closest(initial_values[0], pi3k_value)]
	if( a2_2.size == 0):
		a2_2 = mfs[2][0][ find_closest(initial_values[2], time_value)]

	a2 = max(a2_1, a2_2)
	c2 = np.fmin( np.linspace(a2, a2, 100), mfs[1][0])

	c_com = np.fmax(c1,c2)
	return fuzz.defuzz( initial_values[1], c_com, 'centroid')

def rule_base(b, f_mat):
    """
    Returns y values of output by clipping by an amount of output activations for output fuzzy subsets
    arguments:
    f_mat - rule_strength matrix
    b     - b[2] , y values of output fuzzy subsets

E / DEL_E|         NM      ||       NS        ||        Z         ||       PS        ||       PM
----------------------------------------------------------------------------------------------------------
    NM   | f_mat[0][0] NM  || f_mat[0][1] NM  ||  f_mat[0][2] NS  || f_mat[0][3] Z   || f_mat[0][4] PS
    NS   | f_mat[1][0] NM  || f_mat[1][1] NM  ||  f_mat[1][2] NS  || f_mat[1][3] PS  || f_mat[1][4] PM
    Z    | f_mat[2][0] NM  || f_mat[2][1] NS  ||  f_mat[2][2] Z   || f_mat[2][3] PS  || f_mat[2][4] PM
    PS   | f_mat[3][0] NM  || f_mat[3][1] NS  ||  f_mat[3][2] PS  || f_mat[3][3] PM  || f_mat[3][4] PM
    PM   | f_mat[4][0] NS  || f_mat[4][1] Z   ||  f_mat[4][2] PS  || f_mat[4][3] PM  || f_mat[4][4] PM

    """
    NM = max(f_mat[0][0], f_mat[0][1], f_mat[1][0], f_mat[1][1], f_mat[2][0], f_mat[3][0])
    b[2][0] = np.fmin(NM, b[2][0])
    NS = max(f_mat[0][2], f_mat[1][2], f_mat[2][1], f_mat[3][1], f_mat[4][0])
    b[2][1] = np.fmin(NS, b[2][1])
    Z = max(f_mat[0][3], f_mat[2][2], f_mat[4][1])
    b[2][2] = np.fmin(Z, b[2][2])
    PS = max(f_mat[0][4], f_mat[1][3], f_mat[2][3], f_mat[3][2], f_mat[4][2])
    b[2][3] = np.fmin(PS, b[2][3])
    PM = max(f_mat[1][4], f_mat[2][4], f_mat[3][4], f_mat[3][3], f_mat[4][3], f_mat[4][4])
    b[2][4] = np.fmin(PM, b[2][4])

    return b[2]

def fuzzy_min(x, a, y, b):
    """
    Finds minimum between fuzzy set a in universe x and fuzzy set b in
    universe y.

    Parameters
    ----------
    x : 1d array, length N
        Universe variable for fuzzy set a.
    a : 1d array, length N
        Fuzzy set for universe x.
    y : 1d array, length M
        Universe variable for fuzzy set b.
    b : 1d array, length M
        Fuzzy set for universe y.

    Returns
    -------
    z : 1d array
        Output variable.
    mfz : 1d array
        Fuzzy membership set for variable z.

    Note
    ----
    Uses Zadeh's Extension Principle from Ross, Fuzzy Logic w/Engineering
    Applications, (2010), pp.414, Eq. 12.17.

    """
    # a and x, and b and y, are formed into (MxN) matrices.  The former has
    # identical rows; the latter identical identical columns.
    n = len(b)
    aa = np.dot(np.atleast_2d(a).T, np.ones((1, n)))
    x = np.dot(np.atleast_2d(x).T, np.ones((1, n)))
    m = len(a)
    bb = np.dot(np.ones((m, 1)), np.atleast_2d(b))
    y = np.dot(np.ones((m, 1)), np.atleast_2d(y))

    # Take the element-wise minimum
    zz = np.fmin(x, y).ravel()
    zz_index = np.argsort(zz)
    zz = np.sort(zz)

    # Array min() operation
    c = np.fmin(aa, bb).ravel()
    c = c[zz_index]

    # Initialize loop
    z, mfz = np.empty(0), np.empty(0)
    idx = 0

    for i in range(len(c)):
        index = np.nonzero(zz == zz[idx])[0]
        z = np.hstack((z, zz[idx]))
        mfz = np.hstack((mfz, c[index].max()))
        if zz[idx] == zz.max():
            break
        idx = index.max() + 1

    return z, mfz

	def work(self, input_items, output_items):
		if self.mlimstage:
			for x in xrange(0, len(input_items)):
				try:
					if self.input_blocks[x] is self.null_src:
						continue
				except Exception as e:
					print 'x:%s len(input_items):%s len(input_blocks):%s' % (x, len(input_items), len(self.input_blocks))
					raise e
				ia = input_items[x]
				numpy.fmin(ia, self.limvals[x], ia)
				numpy.multiply(ia, self.mulvals[x], ia)
				input_items[x] = ia

		cur = None
		for x in xrange(0, len(input_items)):
			if self.input_blocks[x] is self.null_src:
				continue
			if cur is None:
				cur = input_items[x]
			else:
				numpy.add(cur, input_items[x], cur)
		bsize = len(output_items[0])
		if cur is None:
			output_items[0][:] = numpy.zeros(bsize, self.dtype)
		else:
			output_items[0][:] = cur
		return bsize

def periodic_3d_distance(x1, y1, z1, x2, y2, z2, Lbox):
    """
    Function computes the distance between two sets of coordinates 
    with the same number of points, accounting for PBCs. 

    Parameters 
    ------------
    x1, y1, z1 : array_like 
        Length-Npts arrays storing Cartesian coordinates 

    x2, y2, z2 : array_like 
        Length-Npts arrays storing Cartesian coordinates 

    Lbox : float 
        Box length defining the periodic boundary conditions 

    Returns 
    --------
    r : array_like 
        Length-Npts array storing the 3d distance between the input 
        points, accounting for box periodicity. 
    """
    dx = np.fabs(x1 - x2)
    dx = np.fmin(dx, Lbox - dx)
    dy = np.fabs(y1 - y2)
    dy = np.fmin(dy, Lbox - dy)
    dz = np.fabs(z1 - z2)
    dz = np.fmin(dz, Lbox - dz)
    return np.sqrt(dx*dx+dy*dy+dz*dz)

 def test_fmin(self):
     from numpy import fmin, array
     nnan, nan, inf, ninf = float('-nan'), float('nan'), float('inf'), float('-inf')
     a = array((complex(ninf, 10), complex(10, ninf),
                complex( inf, 10), complex(10,  inf),
                5+5j, 5-5j, -5+5j, -5-5j,
                0+5j, 0-5j, 5, -5,
                complex(nan, 0), complex(0, nan)), dtype = complex)
     b = [inf]*a.size
     res = [a[0 ], a[1 ], b[2 ], a[3 ],
            a[4 ], a[5 ], a[6 ], a[7 ],
            a[8 ], a[9 ], a[10], a[11],
            b[12], b[13]]
     assert (fmin(a, b) == res).all()
     b = [ninf]*a.size
     res = [b[0 ], b[1 ], b[2 ], b[3 ],
            b[4 ], b[5 ], b[6 ], b[7 ],
            b[8 ], b[9 ], b[10], b[11],
            b[12], b[13]]
     assert (fmin(a, b) == res).all()
     b = [0]*a.size
     res = [a[0 ], b[1 ], b[2 ], b[3 ],
            b[4 ], b[5 ], a[6 ], a[7 ],
            b[8 ], a[9 ], b[10], a[11],
            b[12], b[13]]
     assert (fmin(a, b) == res).all()

def clip_to_window(boxlist, window):
  """Clip bounding boxes to a window.

  This op clips input bounding boxes (represented by bounding box
  corners) to a window, optionally filtering out boxes that do not
  overlap at all with the window.

  Args:
    boxlist: BoxList holding M_in boxes
    window: a numpy array of shape [4] representing the
            [y_min, x_min, y_max, x_max] window to which the op
            should clip boxes.

  Returns:
    a BoxList holding M_out boxes where M_out <= M_in
  """
  y_min, x_min, y_max, x_max = np.array_split(boxlist.get(), 4, axis=1)
  win_y_min = window[0]
  win_x_min = window[1]
  win_y_max = window[2]
  win_x_max = window[3]
  y_min_clipped = np.fmax(np.fmin(y_min, win_y_max), win_y_min)
  y_max_clipped = np.fmax(np.fmin(y_max, win_y_max), win_y_min)
  x_min_clipped = np.fmax(np.fmin(x_min, win_x_max), win_x_min)
  x_max_clipped = np.fmax(np.fmin(x_max, win_x_max), win_x_min)
  clipped = np_box_list.BoxList(
      np.hstack([y_min_clipped, x_min_clipped, y_max_clipped, x_max_clipped]))
  clipped = _copy_extra_fields(clipped, boxlist)
  areas = area(clipped)
  nonzero_area_indices = np.reshape(np.nonzero(np.greater(areas, 0.0)),
                                    [-1]).astype(np.int32)
  return gather(clipped, nonzero_area_indices)

def fuzzy_min(x, A, y, B):
    """
    Finds minimum between fuzzy set A in universe x and fuzzy set B in
    universe y.

    Parameters
    ----------
    x : 1d array, length N
        Universe variable for fuzzy set A.
    A : 1d array, length N
        Fuzzy set for universe x.
    y : 1d array, length M
        Universe variable for fuzzy set B.
    B : 1d array, length M
        Fuzzy set for universe y.

    Returns
    -------
    z : 1d array
        Output variable.
    mfz : 1d array
        Fuzzy membership set for variable z.

    Note
    ----
    Uses Zadeh's Extension Principle from Ross, Fuzzy Logic w/Engineering
    Applications, (2010), pp.414, Eq. 12.17.

    """
    # A and x, and B and y, are formed into (MxN) matrices.  The former has
    # identical rows; the latter identical identical columns.
    N = len(B)
    AA = np.dot(np.atleast_2d(A).T, np.ones((1, N)))
    X = np.dot(np.atleast_2d(x).T, np.ones((1, N)))
    M = len(A)
    BB = np.dot(np.ones((M, 1)), np.atleast_2d(B))
    Y = np.dot(np.ones((M, 1)), np.atleast_2d(y))

    # Take the element-wise minimum
    Z = np.fmin(X, Y).ravel()
    Z_index = np.argsort(Z)
    Z = np.sort(Z)

    # Array min() operation
    C = np.fmin(AA, BB).ravel()
    C = C[Z_index]

    # Initialize loop
    z, mfz = np.empty(0), np.empty(0)
    idx = 0

    for i in range(len(C)):
        index = np.nonzero(Z == Z[idx])[0]
        z = np.hstack((z, Z[idx]))
        mfz = np.hstack((mfz, C[index].max()))
        if Z[idx] == Z.max():
            break
        idx = index.max() + 1

    return z, mfz

def compute_pi3k(egfr_value, erk_value, time_value, initial_values, mfs):
	"""Rules ---
	if egfr is high and erk is low and time is high then pi3k is high
	if egfr is low or erk is high or time is low then pi3k is low"""

	a1_1 = mfs[0][1][initial_values[0] == egfr_value] #egfr_high[egfr == egfr_value]
	a1_2 = mfs[1][0][ initial_values[1] == erk_value] #erk_low[erk == erk_value]
	a1_3 = mfs[3][1][ initial_values[3] == time_value]

	if( a1_1.size == 0):
		a1_1 = mfs[0][1][ find_closest(initial_values[0], egfr_value)]
	if( a1_2.size == 0):
		a1_2 = mfs[1][0][ find_closest(initial_values[1], erk_value)]
	if( a1_3.size == 0):
		a1_3 = mfs[3][1][ find_closest(initial_values[3], time_value)]

	a1 = min(a1_1 , a1_2, a1_3)
	c1 = np.fmin( np.linspace(a1, a1, 100), mfs[2][1])

	a2_1 = mfs[0][0][ initial_values[0] == egfr_value] #egfr_low[egfr == egfr_value]
	a2_2 = mfs[1][1][ initial_values[1] == erk_value] #erk_high[erk == erk_value]
	a2_3 = mfs[3][0][ initial_values[3] == time_value]

	if( a2_1.size == 0):
		a2_1 = mfs[0][0][ find_closest(initial_values[0], egfr_value)]
	if( a2_2.size == 0):
		a2_2 = mfs[1][1][ find_closest(initial_values[1], erk_value)]
	if( a2_3.size == 0):
		a2_3 = mfs[3][0][ find_closest(initial_values[3], time_value)]

	a2 = max(a2_1 , a2_2, a2_3)
	c2 = np.fmin( np.linspace(a2, a2, 100), mfs[2][0] )

	c_com = np.fmax(c1, c2)
	return fuzz.defuzz(initial_values[2], c_com, 'centroid')

def newton(B, b, g, x0, eps, kmax, s):
    k = 0
    x = np.copy(x0)
    x1 = np.copy(x0)
    err = eps + 1
    B1 = np.copy(B.toarray())
    while (k < kmax and err > eps):
        k = k + 1
        F = np.fmin(B.dot(x) - b, x - g)
        Fp = csr_matrix(B, copy=True)

        """
        #############
        F1 = np.fmin(np.dot(B1, x1) - b, x1 - g)
        Fp1 = np.eye(len(B1))
        i1 = np.where((np.dot(B1, x1) - b) <= (x1 - g))
        Fp1[i1, :] = B1[i1, :]
        x1 = x1 - np.linalg.solve(Fp1, F1)
        #############
        """

        i = np.where((B.dot(x) - b) > (x - g))
        for i_row in i[0][:]:
            csr_row_set_nz_to_val(Fp, i_row) # set row to identity
        x = x - sparse.linalg.spsolve(Fp, F)
        #x = x - np.linalg.solve(Fp, F)
        err = np.linalg.norm(np.fmin(B.dot(x) - b , x - g), np.inf)
        err1 = np.linalg.norm(np.fmin(np.dot(B1, x1) - b, x1 - g), np.inf)

        #txt = "Scheme (k=" + str(k) + ")"
        #plt.plot(s, x[992:1024], label=txt )
        #plt.show()
    return x

def compute_erk_change(raf_value, time_value, initial_values, mfs):
	"""Rules-
		If raf is high and time is high then positive_change_erk is high
		If raf is high1 and time is low then positive_change_erk is low
		If raf is low then positive_change_erk is low
		If raf is low and time is high then negative_change_erk is high
		If raf is low and time is low then negative_change_erk is low"""

	
	#Antecedent 1
	f = interp1d(initial_values[0][0], mfs[0][1])
	a1_1 = f(raf_value) #raf_high[raf == raf_value]
	f = interp1d(initial_values[2], mfs[2][1])
	a1_2 = f(time_value) #time_high[time == time_value]

	a1 = min(a1_1, a1_2)
	c1 = np.fmin( a1, mfs[1][3]) #mfs[1][3] is positive_change_erk_high

	#Antecedent 2
	f = interp1d(initial_values[0][0], mfs[0][6])
	a2_1 = f(raf_value)
	f = interp1d(initial_values[2], mfs[2][0]) #time_low[time == time_value]
	a2_2 = f(time_value)

	a2 = min(a2_1, a2_2)
	c2 = np.fmin( a2, mfs[1][2]) #mfs[1][2] is positive_change_raf_low

	c_com_positive = np.fmax(c1,c2)

	f = interp1d(initial_values[0][0], mfs[0][0])
	a3 = f(raf_value)
	c3 = np.fmin(a3, mfs[1][2])

	c_com_positive = np.fmax(c_com_positive, c3)
	pos_change = fuzz.defuzz( initial_values[1][1], c_com_positive, 'centroid') #initial_values[1][1] is positive_change_erk

	###Negative Change

	#Antecedent 3
	'''f = interp1d(initial_values[0][0], mfs[0][0])
	a3_1 = f(raf_value) #raf_low[raf == raf_value]
	a3_2 = a1_2 #time_high[time == time_value]

	a3 = min(a3_1,a3_2)
	c3 = np.fmin(a3, mfs[1][5]) #mfs[1][3] is negative_change_erk_high

	#Antecedent 4
	a4_1 = a3_1 #raf_low[raf == raf_value]
	a4_2 = a2_2 #time_low[time == time_value]

	a4 = min(a4_1, a4_2)
	c4 = np.fmin(a4, mfs[1][4]) #mfs[1][4] is negative_change_erk_low

	c_com_negative = np.fmax(c3, c4)
	neg_change = fuzz.defuzz(initial_values[1][2], c_com_negative, 'centroid') #initial_values[1][2] is negative_change_erk'''
	
	#print pos_change, neg_change
	#print pos_change
	return pos_change 

    def test_half_ufuncs(self):
        """Test the various ufuncs"""

        a = np.array([0, 1, 2, 4, 2], dtype=float16)
        b = np.array([-2, 5, 1, 4, 3], dtype=float16)
        c = np.array([0, -1, -np.inf, np.nan, 6], dtype=float16)

        assert_equal(np.add(a, b), [-2, 6, 3, 8, 5])
        assert_equal(np.subtract(a, b), [2, -4, 1, 0, -1])
        assert_equal(np.multiply(a, b), [0, 5, 2, 16, 6])
        assert_equal(np.divide(a, b), [0, 0.199951171875, 2, 1, 0.66650390625])

        assert_equal(np.equal(a, b), [False, False, False, True, False])
        assert_equal(np.not_equal(a, b), [True, True, True, False, True])
        assert_equal(np.less(a, b), [False, True, False, False, True])
        assert_equal(np.less_equal(a, b), [False, True, False, True, True])
        assert_equal(np.greater(a, b), [True, False, True, False, False])
        assert_equal(np.greater_equal(a, b), [True, False, True, True, False])
        assert_equal(np.logical_and(a, b), [False, True, True, True, True])
        assert_equal(np.logical_or(a, b), [True, True, True, True, True])
        assert_equal(np.logical_xor(a, b), [True, False, False, False, False])
        assert_equal(np.logical_not(a), [True, False, False, False, False])

        assert_equal(np.isnan(c), [False, False, False, True, False])
        assert_equal(np.isinf(c), [False, False, True, False, False])
        assert_equal(np.isfinite(c), [True, True, False, False, True])
        assert_equal(np.signbit(b), [True, False, False, False, False])

        assert_equal(np.copysign(b, a), [2, 5, 1, 4, 3])

        assert_equal(np.maximum(a, b), [0, 5, 2, 4, 3])
        x = np.maximum(b, c)
        assert_(np.isnan(x[3]))
        x[3] = 0
        assert_equal(x, [0, 5, 1, 0, 6])
        assert_equal(np.minimum(a, b), [-2, 1, 1, 4, 2])
        x = np.minimum(b, c)
        assert_(np.isnan(x[3]))
        x[3] = 0
        assert_equal(x, [-2, -1, -np.inf, 0, 3])
        assert_equal(np.fmax(a, b), [0, 5, 2, 4, 3])
        assert_equal(np.fmax(b, c), [0, 5, 1, 4, 6])
        assert_equal(np.fmin(a, b), [-2, 1, 1, 4, 2])
        assert_equal(np.fmin(b, c), [-2, -1, -np.inf, 4, 3])

        assert_equal(np.floor_divide(a, b), [0, 0, 2, 1, 0])
        assert_equal(np.remainder(a, b), [0, 1, 0, 0, 2])
        assert_equal(np.divmod(a, b), ([0, 0, 2, 1, 0], [0, 1, 0, 0, 2]))
        assert_equal(np.square(b), [4, 25, 1, 16, 9])
        assert_equal(np.reciprocal(b), [-0.5, 0.199951171875, 1, 0.25, 0.333251953125])
        assert_equal(np.ones_like(b), [1, 1, 1, 1, 1])
        assert_equal(np.conjugate(b), b)
        assert_equal(np.absolute(b), [2, 5, 1, 4, 3])
        assert_equal(np.negative(b), [2, -5, -1, -4, -3])
        assert_equal(np.positive(b), b)
        assert_equal(np.sign(b), [-1, 1, 1, 1, 1])
        assert_equal(np.modf(b), ([0, 0, 0, 0, 0], b))
        assert_equal(np.frexp(b), ([-0.5, 0.625, 0.5, 0.5, 0.75], [2, 3, 1, 3, 2]))
        assert_equal(np.ldexp(b, [0, 1, 2, 4, 2]), [-2, 10, 4, 64, 12])

    def coord_to_px(self, x, y, latlon=False, rounded=True, check_valid=True):
        """ Convert x,y coordinates into pixel coordinates of raster.

        x,y may be either in native coordinate system of raster or lat/lon.

        Parameters:
            x : float, x coordinate to convert.
            y : float, y coordinate to convert.
            latlon : boolean, default False. Set as True if bounds in lat/lon.
            rounded : if set to True, return the rounded pixel coordinates, otherwise return the float values
            check_valid : bool, if set to True, will check that all pixels are in the valid range. 

        Returns:
            (x_pixel,y_pixel)

        """

        # Convert coordinates to map system if provided in lat/lon and image
        # is projected (rather than geographic)
        if latlon == True and self.proj != None:
            x, y = self.proj(x, y)

        # Shift to the centre of the pixel
        x = np.array(x - self.xres / 2)
        y = np.array(y - self.yres / 2)

        g0, g1, g2, g3, g4, g5 = self.trans
        if g2 == 0:
            xPixel = (x - g0) / float(g1)
            yPixel = (y - g3 - xPixel * g4) / float(g5)
        else:
            xPixel = (y * g2 - x * g5 + g0 * g5 - g2 * g3) / float(g2 * g4 - g1 * g5)
            yPixel = (x - g0 - xPixel * g1) / float(g2)

        # Round if required
        if rounded == True:
            xPixel = np.round(xPixel)
            yPixel = np.round(yPixel)

        if check_valid == False:
            return xPixel, yPixel

        # Check that pixel location is not outside image dimensions
        nx = self.ds.RasterXSize
        ny = self.ds.RasterYSize

        xPixel_new = np.copy(xPixel)
        yPixel_new = np.copy(yPixel)
        xPixel_new = np.fmin(xPixel_new, nx)
        yPixel_new = np.fmin(yPixel_new, ny)
        xPixel_new = np.fmax(xPixel_new, 0)
        yPixel_new = np.fmax(yPixel_new, 0)

        if np.any(xPixel_new != xPixel) or np.any(yPixel_new != yPixel):
            print ("Warning : some points are out of domain for file")

        return xPixel_new, yPixel_new

def compute_akt_change(pi3k_value, time_value, initial_values, mfs):
	"""Rules-
		If pi3k is high and time is high then positive_change_akt is high
		If pi3k is high1 and time is low then positive_change_akt is low
		If pi3k is low then positive_change_pi3k is low
		If pi3k is low and time is high then negative_change_akt is high
		If pi3k is low and time is low then negative_change_akt is low"""

	###Positive Change
	#Antecedent 1
	f = interp1d(initial_values[0][0], mfs[0][1])	
	a1_1 = f(pi3k_value) #pi3k_high[pi3k == pi3k_value]
	f = interp1d(initial_values[2], mfs[2][1])
	a1_2 = f(time_value) #time_high[time == time_value]

	a1 = min(a1_1, a1_2)
	c1 = np.fmin(a1, mfs[1][3]) #positive_change_akt is high

	#Antecedent 2
	f = interp1d(initial_values[0][0], mfs[0][6])
	a2_1 = f(pi3k_value) #pi3k_high[pi3k == pi3k_value]
	f = interp1d(initial_values[2], mfs[2][0])
	a2_2 = f(time_value) #time_low[time == time_value]

	a2 = min(a2_1, a2_2)
	c2 = np.fmin( a2, mfs[1][2]) #positive_change_akt is low

	c_com_positive = np.fmax(c1,c2)

	f = interp1d(initial_values[0][0], mfs[0][0])
	a3 = f(pi3k_value)
	c3 = np.fmin(a3, mfs[1][2])
	c_com_positive = np.fmax(c_com_positive, c3)
	pos_change = fuzz.defuzz( initial_values[1][1], c_com_positive, 'centroid') #initial_values[1][1] is positive_change_akt

	###Negative Change

	#Antecedent 3
	'''f = interp1d(initial_values[0][0], mfs[0][0])
	a3_1 = f(pi3k_value) #pi3k_low[pi3k == pi3k_value]
	a3_2 = a1_2 #time_high[time == time_value]

	a3 = min(a3_1, a3_2)
	c3 = np.fmin(a3, mfs[1][5]) #mfs[1][5] is negative_change_akt_high

	#Antecedent 4
	a4_1 = a3_1 #pi3k_low[pi3k == pi3k_value]
	a4_2 = a2_2 #time_low[time == time_value]

	a4 = min(a4_1, a4_2)
	c4 = np.fmin(a4, mfs[1][4]) #mfs[1][4] is negative_change_akt_low

	c_com_negative = np.fmax(c3, c4)
	neg_change = fuzz.defuzz(initial_values[1][2], c_com_negative, 'centroid') #initial_values[1][2] is negative_change_akt'''

	return pos_change 

def calaculate_akt(pi3k_mfs, akt_mfs, akt, pi3k_index):
	a1 = pi3k_mfs[0][pi3k_index]
	c1 = np.fmin(a1, akt_mfs[0])
	a2 = pi3k_mfs[1][pi3k_index]
	c2 = np.fmin(a2, akt_mfs[1])
	c_com = np.fmax(c1, c2)
	try:
		akt_val = fuzz.defuzz(akt, c_com, 'centroid')
	except AssertionError as e:
		akt_val = 0
	return akt_val

def calculate_egfr(time_mfs, egfr_mfs, egfr, time_index):
	a1 = time_mfs[0][time_index]
	c1 = np.fmin(a1, egfr_mfs[0])
	a2 = time_mfs[1][time_index]
	c2 = np.fmin(a2, egfr_mfs[1])
	c_com = np.fmax(c1, c2)
	try:
		egfr_val = fuzz.defuzz(egfr, c_com, 'centroid')
	except AssertionError as e:
		egfr_val = 0
	return egfr_val

def calculate_erk(raf_mfs, erk_mfs, erk, raf_index):
	a1 = raf_mfs[0][raf_index]
	c1 = np.fmin(a1, erk_mfs[0])
	a2 = raf_mfs[1][raf_index]
	c2 = np.fmin(a2, erk_mfs[1])
	c_com = np.fmax(c1, c2)
	try:
		erk_val = fuzz.defuzz(erk, c_com, 'centroid')
	except AssertionError as e:
		erk_val = 0
	return erk_val

def test_contrast():
    a = np.r_[0, 0, 0, 0.3, 0.7, 1, 0.9, 0]
    z = contrast(a, 1.8)

    # Legacy slower code which should produce identical result
    p = 1.8
    m = 0.5
    ymin = np.fmin(a, m)
    ymax = np.fmax(a, m)
    w = np.arange(len(a))
    wmax = w[ymax > m]
    wmin = w[ymax <= m]
    ymin = 2 ** (p - 1) * ymin ** p
    ymax = 1 - 2 ** (p - 1) * (1 - ymax) ** p
    ymin[wmax] = 0
    ymax[wmin] = 0

    assert_allclose(z, ymin + ymax)

    # Legacy slower code which should produce identical result
    p = 0.5
    m = 0.5
    z = contrast(a, 0.5)
    ymin = np.fmin(a, m)
    ymax = np.fmax(a, m)
    w = np.arange(len(a))
    wmax = w[ymax > m]
    wmin = w[ymax <= m]
    ymin = 2 ** (p - 1) * ymin ** p
    ymax = 1 - 2 ** (p - 1) * (1 - ymax) ** p
    ymin[wmax] = 0
    ymax[wmin] = 0

    assert_allclose(z, ymin + ymax)

    # Legacy slower code which should produce identical result
    p = 2.
    m = 0.5
    z = contrast(a, 2.)
    ymin = np.fmin(a, m)
    ymax = np.fmax(a, m)
    w = np.arange(len(a))
    wmax = w[ymax > m]
    wmin = w[ymax <= m]
    ymin = 2 ** (p - 1) * ymin ** p
    ymax = 1 - 2 ** (p - 1) * (1 - ymax) ** p
    ymin[wmax] = 0
    ymax[wmin] = 0

    assert_allclose(z, ymin + ymax)

def rule_actividad(value, graficar=False):

    # recibe los datos para poder seguir el proceso difuso para encontrar la pertenencia
    # retorna el valor al que pertenece en la clase segun la hora
    mf_actividad = generar_actividad(False)
    mf_calorico = generar_calorico(False)

    # se usa para encontrar los grados de pertenencia del valor, borrificacion
    actividad_nivel_rest = fuzz.interp_membership(mf_actividad['intensity'], mf_actividad['rest'], value)
    actividad_nivel_std = fuzz.interp_membership(mf_actividad['intensity'], mf_actividad['active'], value)
    actividad_nivel_work = fuzz.interp_membership(mf_actividad['intensity'], mf_actividad['workout'], value)

    # regla: si rest -> low
    rest_activation = np.fmin(actividad_nivel_rest, mf_calorico['low'])
    # regla: si active -> std
    active_activation = np.fmin(actividad_nivel_std, mf_calorico['standard'])
    # regla: si workout -> high
    workout_activation = np.fmin(actividad_nivel_work, mf_calorico['high'])

    # deborrificacion
    agregado = np.fmax(rest_activation, np.fmax(active_activation, workout_activation))
    caloric = fuzz.defuzz(mf_calorico['caloric'], agregado, 'centroid')


    # # graficar
    # if graficar:
    #     select_caloric = fuzz.interp_membership(mf_calorico['caloric'], agregado, caloric)
    #     caloric0 = np.zeros_like(mf_calorico['caloric'])
    #     fig, ax0 = plt.subplots(figsize=(8, 3))
    #
    #     ax0.plot(mf_calorico['caloric'], mf_calorico['low'], 'b', linewidth=0.5, linestyle='--', label='Low')
    #     ax0.plot(mf_calorico['caloric'], mf_calorico['standard'], 'g', linewidth=0.5, linestyle='--', label='Standard')
    #     ax0.plot(mf_calorico['caloric'], mf_calorico['high'], 'r', linewidth=0.5, linestyle='--', label='High')
    #     ax0.fill_between(mf_calorico['caloric'], caloric0, agregado, facecolor='Orange', alpha=0.7)
    #     ax0.plot([caloric, caloric], [0, select_caloric], 'k', linewidth=1.5, alpha=0.9)
    #     ax0.set_title('Dish Classification and Result (line)')
    #     ax0.legend()
    #     # Turn off top/right axes
    #     for ax in (ax0,):
    #         ax.spines['top'].set_visible(False)
    #         ax.spines['right'].set_visible(False)
    #         ax.get_xaxis().tick_bottom()
    #         ax.get_yaxis().tick_left()
    #
    #     plt.tight_layout()
    #     #fig.savefig('graphs/result_calorico.png', bbox_inches='tight')
    # 
    return caloric, "caloric"