def test_SolveUsingSpline_3D(self):
        xKnotPointsShouldBe = numpy.array([0.607, 0.607, 0.607, 3.017, 3.017, 3.017])
        yKnotPointsShouldBe = numpy.array([1.984, 1.984, 1.984, 3.153, 3.153, 3.153])
        coefficientsShouldBe = numpy.array(
            [2.33418963, 1.80079612, 5.07902936, 0.54445029, 1.04110843, 2.14180324, 0.26992805, 0.39148852, 0.8177307]
        )
        testEvaluationShouldBe = numpy.array([0.76020577997])
        model = pyeq2.Models_3D.Spline.Spline(inSmoothingFactor=1.0, inXOrder=2, inYOrder=2)
        pyeq2.dataConvertorService().ConvertAndSortColumnarASCII(DataForUnitTests.asciiDataInColumns_3D, model, False)
        fittedParameters = pyeq2.solverService().SolveUsingSpline(model)

        # example of later using the saved spline knot points and coefficients
        unFittedSpline = scipy.interpolate.fitpack2.SmoothBivariateSpline(
            model.dataCache.allDataCacheDictionary["X"],
            model.dataCache.allDataCacheDictionary["Y"],
            model.dataCache.allDataCacheDictionary["DependentData"],
            s=model.smoothingFactor,
            kx=model.xOrder,
            ky=model.yOrder,
        )
        unFittedSpline.tck = fittedParameters
        testEvaluation = unFittedSpline.ev(2.5, 2.5)

        self.assertTrue(numpy.allclose(testEvaluation, testEvaluationShouldBe, rtol=1.0e-10, atol=1.0e-300))
        self.assertTrue(numpy.equal(fittedParameters[0], xKnotPointsShouldBe).all())
        self.assertTrue(numpy.equal(fittedParameters[1], yKnotPointsShouldBe).all())
        self.assertTrue(numpy.allclose(fittedParameters[2], coefficientsShouldBe, rtol=1.0e-06, atol=1.0e-300))

 def test_manual_bounds(self, cuda=False):
     device = torch.device("cuda") if cuda else torch.device("cpu")
     for dtype in (torch.float, torch.double):
         # get a test module
         train_x = torch.tensor([[1.0, 2.0, 3.0]], device=device, dtype=dtype)
         train_y = torch.tensor([4.0], device=device, dtype=dtype)
         likelihood = GaussianLikelihood()
         model = ExactGP(train_x, train_y, likelihood)
         model.covar_module = RBFKernel(ard_num_dims=3)
         model.mean_module = ConstantMean()
         model.to(device=device, dtype=dtype)
         mll = ExactMarginalLogLikelihood(likelihood, model)
         # test the basic case
         x, pdict, bounds = module_to_array(
             module=mll, bounds={"model.covar_module.raw_lengthscale": (0.1, None)}
         )
         self.assertTrue(np.array_equal(x, np.zeros(5)))
         expected_sizes = {
             "likelihood.noise_covar.raw_noise": torch.Size([1]),
             "model.covar_module.raw_lengthscale": torch.Size([1, 3]),
             "model.mean_module.constant": torch.Size([1]),
         }
         self.assertEqual(set(pdict.keys()), set(expected_sizes.keys()))
         for pname, val in pdict.items():
             self.assertEqual(val.dtype, dtype)
             self.assertEqual(val.shape, expected_sizes[pname])
             self.assertEqual(val.device.type, device.type)
         lower_exp = np.full_like(x, 0.1)
         for p in ("likelihood.noise_covar.raw_noise", "model.mean_module.constant"):
             lower_exp[_get_index(pdict, p)] = -np.inf
         self.assertTrue(np.equal(bounds[0], lower_exp).all())
         self.assertTrue(np.equal(bounds[1], np.full_like(x, np.inf)).all())

 def _get_ind_under_point(self, event):
     'get the index of the vertex under point if within epsilon tolerance'
     try:
         x, y = zip(*self._poly.xy)
         
         # display coords
         xt, yt = self._poly.get_transform().numerix_x_y(x, y)
         d = np.sqrt((xt-event.x)**2 + (yt-event.y)**2)
         indseq = np.nonzero(np.equal(d, np.amin(d)))
         ind = indseq[0]
     
         if d[ind]>=self._epsilon:
             ind = None
     
         return ind
     except:
         # display coords
         xy = np.asarray(self._poly.xy)
         xyt = self._poly.get_transform().transform(xy)
         xt, yt = xyt[:, 0], xyt[:, 1]
         d = np.sqrt((xt-event.x)**2 + (yt-event.y)**2)
         indseq = np.nonzero(np.equal(d, np.amin(d)))[0]
         ind = indseq[0]
         
         if d[ind]>=self._epsilon:
             ind = None
         
         return ind

 def applyMorphologicalCleaning(self, image):
 	"""
 	Applies a variety of morphological operations to improve the detection
 	of worms in the image.
 	Takes 0.030 s on MUSSORGSKY for a typical frame region
 	Takes 0.030 s in MATLAB too
 	"""
     # start with worm == 1
     image = image.copy()
     segmentation.clear_border(image)  # remove objects at edge (worm == 1)
     # fix defects in the thresholding by closing with a worm-width disk
     # worm == 1
     wormSE = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,
                                        (self.wormDiskRadius+1,
                                        	self.wormDiskRadius+1))
     imcl = cv2.morphologyEx(np.uint8(image), cv2.MORPH_CLOSE, wormSE)
     imcl = np.equal(imcl, 1)
     # fix defects by filling holes
     imholes = ndimage.binary_fill_holes(imcl)
     imcl = np.logical_or(imholes, imcl)
     # fix barely touching regions
     # majority with worm pixels == 1 (median filter same?)
     imcl = nf.median_filter(imcl, footprint=[[1, 1, 1],
                                              [1, 0, 1],
                                              [1, 1, 1]])
     # diag with worm pixels == 0
     imcl = np.logical_not(bwdiagfill(np.logical_not(imcl)))
     # open with worm pixels == 1
     openSE = cv2.getStructuringElement(cv2.MORPH_RECT, (1, 1))
     imcl = cv2.morphologyEx(np.uint8(imcl), cv2.MORPH_OPEN, openSE)
     return np.equal(imcl, 1)

def merge_img_array(arr_1, nband_1, arr_2, nband_2):
    """
    Merge images from difference arrays of 2D images. (ex. one has 7, the other has 2, both of them has 500 temporal images, then 
    we merge to a list of 500 temporal images with 9 bands.) The variable arr_1 is the main one.
    """
    m = arr_1.shape[0]
    n = arr_1.shape[1]
    l1 = arr_1.shape[2]
    l2 = arr_2.shape[2]
    merged_arr = np.empty((m, n, l1+l2))
    
    
    for i in range(int(l1/nband_1)):
        position_1 = i*nband_1 + i*nband_2
        position_2 = position_1 + nband_1
        merged_arr[:, :, position_1 : position_1+nband_1] = arr_1[:, :, i*nband_1 : i*nband_1+nband_1]
        merged_arr[:, :, position_2 : position_2+nband_2] = arr_2[:, :, i*nband_2 : i*nband_2+nband_2]
        
    # Check before done
    clause_1 = np.all(np.equal(merged_arr[:, :, 0:nband_1], arr_1[:, :, 0:nband_1]))
    clause_2 = np.all(np.equal(merged_arr[:, :, -nband_2:], arr_2[:, :, -nband_2:]))
    if (clause_1 == True) and (clause_2 == True):
        #print('done')
        pass
    else:
        raise ValueError('A very specific bad thing happened. Maybe the number of given band is wrong?')
    return merged_arr

def rgb_to_hsv( r,g,b ):

    maxc = numpy.maximum(r,numpy.maximum(g,b))
    minc = numpy.minimum(r,numpy.minimum(g,b))

    v = maxc

    minc_eq_maxc = numpy.equal(minc,maxc)

    # compute the difference, but reset zeros to ones to avoid divide by zeros later.
    ones = numpy.ones((r.shape[0],r.shape[1]))
    maxc_minus_minc = numpy.choose( minc_eq_maxc, (maxc-minc,ones) )

    s = (maxc-minc) / numpy.maximum(ones,maxc)
    rc = (maxc-r) / maxc_minus_minc
    gc = (maxc-g) / maxc_minus_minc
    bc = (maxc-b) / maxc_minus_minc

    maxc_is_r = numpy.equal(maxc,r)
    maxc_is_g = numpy.equal(maxc,g)
    maxc_is_b = numpy.equal(maxc,b)

    h = numpy.zeros((r.shape[0],r.shape[1]))
    h = numpy.choose( maxc_is_b, (h,4.0+gc-rc) )
    h = numpy.choose( maxc_is_g, (h,2.0+rc-bc) )
    h = numpy.choose( maxc_is_r, (h,bc-gc) )

    h = numpy.mod(h/6.0,1.0)

    hsv = numpy.asarray([h,s,v])
    
    return hsv

    def error_type_voxelWise(self, reference):
        '''Returns a dictionary with tagged error type
        for each regions from both reference and this label.'''
        
        self.data = np.where(self.data>0, 1, 0)
        reference.data = np.where(reference.data>0, 1, 0)
        
        interSet = self.intersection(reference).get(1, dict()).get(1, set())
        interArray = np.asarray(list(interSet))
        stats = {'reference': dict(), 'self': dict()}

        sIdx = np.where(self.data>0)
        rIdx = np.where(reference.data>0)
        
        sIdxArray = np.transpose(np.asarray(sIdx))
        rIdxArray = np.transpose(np.asarray(rIdx))
            
        for s in sIdxArray:
            if interArray.any():
                if np.equal(s, interArray).all(axis=1).any():
                    stats['self'][str(s)] = {'type': 'TP', 'regions': 1}
                else:
                    stats['self'][str(s)] = {'type': 'FP'}
            else:
                stats['self'][str(s)] = {'type': 'FP'}

        for r in rIdxArray:
            if interArray.any():
                if not np.equal(r, interArray).all(axis=1).any():
                    stats['reference'][str(r)] = {'type': 'FN'}
            else:
                stats['reference'][str(r)] = {'type': 'FN'}

        return stats

def long_test_2():
    # Set up the market
    bids = np.arange(1, 2, .10)
    offers = np.arange(1.1, 2.1, .10)
    bid_vols = 2000000 * np.ones(10)
    offer_vols = 1000000 * np.ones(10)
    mean_spread = 0.1
    mean_range = 0.0
    # Set up signals for single long trade
    signals  = [1,0,0,0,0,0,0,0,0,-1]
    # TEST 2
    # Test enter and exit with one buy and sell - no position carry, trade only on signal = False
    # Test for opportunistic profit
    # NB - NO CARRY will average closing prices which again distorts pnl in monotonic markets like the test sets so closing pnl will be exaggerated here
    (pnl, position_deltas, position_running, closing_position, closing_pnl, ignored_signals, m2m_pnl) = simulate2.execute_aggressive(ts, bids, offers, bid_vols, offer_vols, signals, currency_pair, signal_window_time=1, min_window_signals=1, min_profit_prct=0.0001, carry_position = False, default_trade_size = 1, max_position=5, fill_function=None, cut_long = -(mean_spread+mean_range)*2, cut_short= -(mean_spread+mean_range)*2, usd_transaction_cost= 0, trade_file='/tmp/testoutshit', take_profit_pct=0.0001)
    # Profit here is convoluted since the price array is artificial and changes very quickly so when averaging the price to close the position the true price is distorted - i.e. first 4 offers average to 1.30
    assert(round(sum(pnl),2) == 0.35)

    pos_deltas_test = [ 1.,  0., -1.,  0.,  0.,  0.,  0.,  0.,  0., 1.]
    pos_deltas_out = np.equal(position_deltas, pos_deltas_test)
    assert(pos_deltas_out.all() == True)
    assert(sum(position_deltas) + closing_position == 0)

    pos_run_test  = [ 1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]
    pos_run_out = np.equal(position_running, pos_run_test)
    assert(pos_run_out.all() == True)
    

    assert(round(closing_pnl,2) == 0.25)
    assert(closing_position == -1.)

def make_features_in_layers(board):
  feature_layers = [] # These are effectively 'colours' for the CNN
  #print(board)
  
  #print("Board mask")
  mask     = np.greater( board[:, :], 0 )*1.
  feature_layers.append( mask.astype('float32') )
  
  # This works out whether each cell is the same as the cell 'above it'
  for shift_down in [1,2,3,4,5,]:
    #print("\n'DOWN' by %d:" % (shift_down,))
    sameness = np.zeros_like(board, dtype='float32')
    
    sameness[:,:-shift_down] = np.equal( board[:, :-shift_down], board[:, shift_down:] )*1.
    #print(sameness)

    feature_layers.append( sameness )
  
  # This works out whether each cell is the same as the cell in to columns 'to the left of it'
  for shift_right in [1,2,3,]:
    #print("\n'RIGHT' by %d:" % (shift_right,))
    sameness = np.zeros_like(board, dtype='float32')
    
    sameness[:-shift_right,:] = np.equal( board[:-shift_right, :], board[shift_right:, :] )*1.
    #print(sameness)

    feature_layers.append( sameness )
  
  stacked = np.dstack( feature_layers )
  return np.rollaxis( stacked, 2, 0 )

def mix_test1():
    signals  = [1,1,1,0,-1,0,-1,0,1,-1]
    mean_spread = 1
    mean_range = 1
    # Reset the market 
    bids = np.arange(1, 2, .10)
    offers = np.arange(1.1, 2.1, .10)
    bid_vols = 2000000 * np.ones(10)
    offer_vols = 1000000 * np.ones(10)

    # TEST 7
    # No take profit - trade only on signals 
    # No cuts - set spread and range high 
    # DO carry position - elimintates averaging problem of closing trde - pnl more transparent
    # No min profit - trade on all signals - not just profitable ones
    # This should trade on ALL signals regardless of profitability and never cut 
    (pnl, position_deltas, position_running, closing_position, closing_pnl, ignored_signals, m2m_pnl) = simulate2.execute_aggressive(ts, bids, offers, bid_vols, offer_vols, signals, currency_pair, signal_window_time=1, min_window_signals=1, min_profit_prct=None, carry_position = True, default_trade_size = 1, max_position=5, fill_function=None, cut_long = -(mean_spread+mean_range)*2, cut_short= -(mean_spread+mean_range)*2, usd_transaction_cost= 0, trade_file='/tmp/testoutshit', take_profit_pct=None)

    assert(round(sum(pnl), 2) == .15)

    pos_deltas_test = [ 1.,  1.,  1.,  0., -3.,  0., -1.,  0.,  1., -1.] 
    pos_deltas_out = np.equal(position_deltas, pos_deltas_test)
    assert(pos_deltas_out.all() == True)

    # End with position since carry_position is True here
    pos_run_test  = [ 1.,  2.,  3.,  3.,  0.,  0., -1.,  -1.,  0., -1.]
    pos_run_out = np.equal(position_running, pos_run_test)
    assert(pos_run_out.all() == True)

    assert(closing_position == -1.)
    assert(closing_pnl == 0)

def mix_test4():
    signals  = [1,1,1,0,-1,0,-1,0,1,-1]
    # Reset the market 
    bids = np.arange(1, 2, .10)
    offers = np.arange(1.1, 2.1, .10)
    bid_vols = 2000000 * np.ones(10)
    offer_vols = 1000000 * np.ones(10)
    mean_spread = 0
    mean_range = 0.0001
    # TEST 10
    # Allow take profit - fully exit profitable position
    # Allow cuts 
    mean_spread = 0
    mean_range = 0.0001
    # DO NOT carry position - close avg price of long/short against avg of past market
    # No min profit - trade on all signals - not just profitable ones
    (pnl, position_deltas, position_running, closing_position, closing_pnl, ignored_signals, m2m_pnl) = simulate2.execute_aggressive(ts, bids, offers, bid_vols, offer_vols, signals, currency_pair, signal_window_time=1, min_window_signals=1, min_profit_prct=None, carry_position = True, default_trade_size = 1, max_position=5, fill_function=None, cut_long = -(mean_spread+mean_range)*2, cut_short= -(mean_spread+mean_range)*2, usd_transaction_cost= 0, trade_file='/tmp/testoutshit', take_profit_pct=0.0001)

    assert(round(sum(pnl), 2) == 0.02)

    pos_deltas_test = [ 1.,  1.,  1., -3., -1.,  1., -1.,  1.,  1., -1.] 
    pos_deltas_out = np.equal(position_deltas, pos_deltas_test)
    assert(pos_deltas_out.all() == True)

    pos_run_test = [ 1.,  2.,  3.,  0., -1.,  0., -1.,  0.,  1.,  0.] 
    pos_run_out = np.equal(position_running, pos_run_test)
    assert(pos_run_out.all() == True)

    assert(closing_position == 0.)
    assert(closing_pnl == 0.)

def short_test3():
    # Set up the market
    # Cut level is offer[6]
    offers = [ 2.1,  2. ,  1.9,  1.8,  1.7,  1.6,  2.5 ,  1.4,  1.3,  1.2]
    bids = [ 2. ,  1.9,  1.8,  1.7,  1.6,  1.5,  1.4,  1.3,  1.2,  1.1]
    bid_vols = 2000000 * np.ones(10)
    offer_vols = 1000000 * np.ones(10)
    mean_spread = 0.1
    mean_range = 0.0
    # Start with short and end with long
    signals = [-1,0,0,0,0,0,0,0,0,1]
    # TEST 6
    # Test short trade with cutoff
    # NB - NO CARRY will average closing prices which again distorts pnl in monotonic markets like the test sets so closing pnl will be exaggerated here
    mean_spread = 0.0
    mean_range = 0.0001
    (pnl, position_deltas, position_running, closing_position, closing_pnl, ignored_signals, m2m_pnl) = simulate2.execute_aggressive(ts, bids, offers, bid_vols, offer_vols, signals, currency_pair, signal_window_time=1, min_window_signals=1, min_profit_prct=0.0001, carry_position = False, default_trade_size = 1, max_position=5, fill_function=None, cut_long = -(mean_spread+mean_range)*2, cut_short= -(mean_spread+mean_range)*2, usd_transaction_cost= 0, trade_file='/tmp/testoutshit', take_profit_pct=None)
    assert(round(sum(pnl), 2) == -0.55)

    pos_deltas_test = [-1.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  -1.]
    pos_deltas_out = np.equal(position_deltas, pos_deltas_test)
    assert(pos_deltas_out.all() == True)

    pos_run_test  = [ -1.,  -1.,  -1.,  -1.,  -1.,  -1.,  0.,  0.,  0.,  0.]
    pos_run_out = np.equal(position_running, pos_run_test)
    assert(pos_run_out.all() == True)

    assert(round(closing_pnl,2) == -0.05)
    assert(closing_position == 1.)

def short_test2():
    # Set up the market
    bids = np.arange(2, 1, -.10)
    offers = np.arange(2.1, 1.1, -.10)
    bid_vols = 2000000 * np.ones(10)
    offer_vols = 1000000 * np.ones(10)
    mean_spread = 0.1
    mean_range = 0.0
    # Start with short and end with long
    signals = [-1,0,0,0,0,0,0,0,0,1]
    # TEST 5
    # Test short trade first - no carry position, trade only on signal = False
    # Short and take profit - then trade last frame long and close position with short (no carry)
    # NB - NO CARRY will average closing prices which again distorts pnl in monotonic markets like the test sets so closing pnl will be exaggerated here
    (pnl, position_deltas, position_running, closing_position, closing_pnl, ignored_signals, m2m_pnl) = simulate2.execute_aggressive(ts, bids, offers, bid_vols, offer_vols, signals, currency_pair, signal_window_time=1, min_window_signals=1, min_profit_prct=0.0001, carry_position = False, default_trade_size = 1, max_position=5, fill_function=None, cut_long = -(mean_spread+mean_range)*2, cut_short= -(mean_spread+mean_range)*2, usd_transaction_cost= 0, trade_file='/tmp/testoutshit', take_profit_pct=0.0001)
    assert(round(sum(pnl), 2) == 0.35)

    pos_deltas_test = [ -1.,  0., 1.,  0.,  0.,  0.,  0.,  0.,  0., -1.]
    pos_deltas_out = np.equal(position_deltas, pos_deltas_test)
    assert(pos_deltas_out.all() == True)

    pos_run_test  = [ -1.,  -1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]
    pos_run_out = np.equal(position_running, pos_run_test)
    assert(pos_run_out.all() == True)

    assert(round(closing_pnl,2) == 0.25)
    assert(closing_position == 1.)

def short_test1():
    # Set up the market
    bids = np.arange(2, 1, -.10)
    offers = np.arange(2.1, 1.1, -.10)
    bid_vols = 2000000 * np.ones(10)
    offer_vols = 1000000 * np.ones(10)
    mean_spread = 0.1
    mean_range = 0.0
    # Start with short and end with long
    signals = [-1,0,0,0,0,0,0,0,0,1]

    # TEST 4
    # Test short trade first - no carry position, trade only on signal = True
    (pnl, position_deltas, position_running, closing_position, closing_pnl, ignored_signals, m2m_pnl) = simulate2.execute_aggressive(ts, bids, offers, bid_vols, offer_vols, signals, currency_pair, signal_window_time=1, min_window_signals=1, min_profit_prct=0.0001, carry_position = False, default_trade_size = 1, max_position=5, fill_function=None, cut_long = -(mean_spread+mean_range)*2, cut_short= -(mean_spread+mean_range)*2, usd_transaction_cost= 0, trade_file='/tmp/testoutshit', take_profit_pct=None)
    assert(round(sum(pnl), 2) == 0.8)

    pos_deltas_test = [ -1.,  0., 0.,  0.,  0.,  0.,  0.,  0.,  0., 1.]
    pos_deltas_out = np.equal(position_deltas, pos_deltas_test)
    assert(pos_deltas_out.all() == True)

    pos_run_test  = [ -1.,  -1.,  -1.,  -1.,  -1.,  -1.,  -1.,  -1.,  -1.,  0.]
    pos_run_out = np.equal(position_running, pos_run_test)
    assert(pos_run_out.all() == True)

    assert(round(closing_pnl,2) == 0.0)
    assert(round(closing_position,2) == 0.0)

 def test_support(self):
     self.assertIsInstance(self.f0.support, Domain)
     self.assertIsInstance(self.f1.support, Domain)
     self.assertIsInstance(self.f2.support, Domain)
     self.assertEqual(self.f0.support.size, 0)
     self.assertTrue(np.equal(self.f1.support,[-1,1]).all())
     self.assertTrue(np.equal(self.f2.support,[-1,2]).all())

	def fit(self,X,y):
		#TODO: check X,y
		
		self.classes = np.unique(y)
		#calculate class prior probabilities: P(y=ck)
		if self.class_prior == None:
                        class_num = len(self.classes)
			if not self.fit_prior:
				self.class_prior = [1.0/class_num for _ in range(class_num)]  #uniform prior
			else:
				self.class_prior = []
				sample_num = float(len(y))
				for c in self.classes:
					c_num = np.sum(np.equal(y,c))
					self.class_prior.append(
                                                (c_num+self.alpha)/(sample_num+class_num*self.alpha))
		
		#calculate Conditional Probability: P( xj | y=ck )
		self.conditional_prob = {}  # like { c0:{ x0:{ value0:0.2, value1:0.8 }, x1:{} }, c1:{...} }
		for c in self.classes:
			self.conditional_prob[c] = {}
			for i in range(len(X[0])):  #for each feature
                                feature = X[np.equal(y,c)][:,i]
				self.conditional_prob[c][i] = self._calculate_feature_prob(feature)
		return self

def calc_stats_class(image, mask = None, index = None):
    ''' funcion auxiliar para el calculo estadistico de cada categoria.'''
    import spectral
    from numpy import zeros, transpose,compress, indices, reshape, not_equal,mean,std
    from spectral.io import typecode
    
    typechar = typecode(image)
    (nrows, ncols, B) = image.shape
    (nr,nc) = mask.shape
    mask_i = numpy.equal(mask, index)
    mask_array = mask.reshape(nr*nc)
    mask_index = numpy.equal(mask_array, index)
    nSamples = sum(mask_index.ravel())
    
    
    inds = transpose(indices((nrows, ncols)), (1, 2, 0))
    inds = reshape(inds, (nrows * ncols, 2))
    inds = compress(not_equal(mask_i.ravel(), 0), inds, 0).astype('h')

    vector=zeros((inds.shape[0], B), 'd')
 
    for i in range(inds.shape[0]):
        x = image[inds[i][0], inds[i][1]]
        vector[i] = x     
    media = mean(vector,axis=0)
    desv = std(vector, axis=0)
    return(media,desv)

def find_intersection(x, tr_bounds, lb, ub):
    """Find intersection of trust-region bounds and initial bounds.

    Returns
    -------
    lb_total, ub_total : ndarray with shape of x
        Lower and upper bounds of the intersection region.
    orig_l, orig_u : ndarray of bool with shape of x
        True means that an original bound is taken as a corresponding bound
        in the intersection region.
    tr_l, tr_u : ndarray of bool with shape of x
        True means that a trust-region bound is taken as a corresponding bound
        in the intersection region.
    """
    lb_centered = lb - x
    ub_centered = ub - x

    lb_total = np.maximum(lb_centered, -tr_bounds)
    ub_total = np.minimum(ub_centered, tr_bounds)

    orig_l = np.equal(lb_total, lb_centered)
    orig_u = np.equal(ub_total, ub_centered)

    tr_l = np.equal(lb_total, -tr_bounds)
    tr_u = np.equal(ub_total, tr_bounds)

    return lb_total, ub_total, orig_l, orig_u, tr_l, tr_u

def issorted(pot, variables=[]):
    """Check whether the variables in the Potential is sorted.

    Parameters :
        pot : Potential :
            The target potential.

        variables : sequence[n_variables, ] or np.ndarray[n_variables, ], optional, default : None :
            The sorted sequence of variables to be compared with.

    Returns :
        issorted : boolean :
            True for pot is sorted, False otherwise.

    Raises :
        None

    Notes :
        If variables is not None, compare the variables of the pot with the
        given variables.
    """
    if variables is None:
        variables = np.array([])
    else:
        variables = np.array(variables)
    originVar = pot.variables
    if np.size(variables) == 0:
        return np.equal(np.sort(originVar), originVar)
    else:
        return np.equal(originVar, variables)

def colon(r1, inc, r2):
    """
      Matlab's colon operator, althought it doesn't although inc is required

    """

    s = np.sign(inc)

    if s == 0:
        return_value = np.zeros(1)
    elif s == 1:
        n = ((r2 - r1) + 2 * np.spacing(r2 - r1)) // inc
        return_value = np.linspace(r1, r1 + inc * n, n + 1)
    else:  # s == -1:
        # NOTE: I think this is slightly off as we start on the wrong end
        # r1 should be exact, not r2
        n = ((r1 - r2) + 2 * np.spacing(r1 - r2)) // np.abs(inc)
        temp = np.linspace(r2, r2 + np.abs(inc) * n, n + 1)
        return_value = temp[::-1]

    # If the start and steps are whole numbers, we should cast as int
    if(np.equal(np.mod(r1, 1), 0) and
       np.equal(np.mod(s, 1), 0) and
       np.equal(np.mod(r2, 1), 0)):
        return return_value.astype(int)
    else:
        return return_value