def empirical_ci(y,alpha=0.05):
	"""Computes an empirical (alpha/2,1-alpha/2) confidence interval for the distributional data in x.

	Parameters:
	------------
	x : numpy array, required
		set of data to produce empirical upper/lower bounds for

	alpha : float, optional
		sets desired CI range

	Returns:
	------------
	lb, ub : floats, lower and upper bounds for x

	"""
	ytilde = sort(y)
	xl = (alpha/2)*len(y)
	xu = (1.0 - alpha/2)*len(y)
	l1 = int(floor(xl))
	l2 = int(ceil(xl))
	u1 = int(floor(xu))
	u2 = int(ceil(xu))
	lb = interp(xl,[l1,l2],[ytilde[l1],ytilde[l2]])
	ub = interp(xu,[u1,u2],[ytilde[u1],ytilde[u2]])
	return lb,ub

def fit(dataset):

    # f = gzip.open('../../../datasets/Mnist/mnist.pkl.gz', 'rb')
    # train_set, _, _ = pickle.load(f)
    # f.close()
    #
    # _, labels = train_set

    # features = pickle.load(open('../../../datasets/Mnist/convnet_train_features.p', 'rb'))
    samples = pickle.load(open('../../../datasets/Mnist/bag_train_features.p', 'rb'))

    features = []
    labels = []
    for sample in samples:
        features.append(sample['features'])
        labels.append(sample['label'])
    features = np.array(features)
    labels = np.array(labels)

    model = Oasis(n_iter=100000, do_psd=True, psd_every=3,
                  save_path="/tmp/gwtaylor/oasis_test").fit(features, labels,
                                                            verbose=True)
    W = model._weights.view()

    W.shape = (np.int(np.sqrt(W.shape[0])), np.int(np.sqrt(W.shape[0])))

    # pickle.dump(W, open('../convnet/oasis_weights.p', 'wb'))
    pickle.dump(W, open('../bag/oasis_weights.p', 'wb'))

    def _TO_DELETE_initialize_drifters(self, driftersPerOceanModel):
        """
        Initialize drifters and attach them for each particle.
        """
        self.driftersPerOceanModel = np.int32(driftersPerOceanModel)
        
        # Define mid-points for the different drifters 
        # Decompose the domain, so that we spread the drifters as much as possible
        sub_domains_y = np.int(np.round(np.sqrt(self.driftersPerOceanModel)))
        sub_domains_x = np.int(np.ceil(1.0*self.driftersPerOceanModel/sub_domains_y))
        self.midPoints = np.empty((driftersPerOceanModel, 2))
        for sub_y in range(sub_domains_y):
            for sub_x in range(sub_domains_x):
                drifter_id = sub_y*sub_domains_x + sub_x
                if drifter_id >= self.driftersPerOceanModel:
                    break
                self.midPoints[drifter_id, 0]  = (sub_x + 0.5)*self.nx*self.dx/sub_domains_x
                self.midPoints[drifter_id, 1]  = (sub_y + 0.5)*self.ny*self.dy/sub_domains_y
              
        # Loop over particles, sample drifters, and attach them
        for i in range(self.numParticles+1):
            drifters = GPUDrifterCollection.GPUDrifterCollection(self.gpu_ctx, self.driftersPerOceanModel,
                                                 observation_variance=self.observation_variance,
                                                 boundaryConditions=self.boundaryConditions,
                                                 domain_size_x=self.nx*self.dx, domain_size_y=self.ny*self.dy)

            initPos = np.empty((self.driftersPerOceanModel, 2))
            for d in range(self.driftersPerOceanModel):
                initPos[d,:] = np.random.multivariate_normal(self.midPoints[d,:], self.initialization_cov_drifters)
            drifters.setDrifterPositions(initPos)
            self.particles[i].attachDrifters(drifters)

    def _calcNeighbours(self):
        """ Calculates the neighbours and creates the bitstrings """
        self._logger.info("Creating neighbourhood bitstrings")  

        bitstrings = []

        # apply BS func on all c alphas in peptide chains
        cas = self._df.loc[(self._df['an'] == 'CA') & (self._df['peptideChain'] == True) & (self._df['surface'] == True), ['chain', 'resi', 'resn', 'inscode']]
        for tpl in cas.itertuples():
            idx, chain, resi, resn, inscode = tpl

            # other c alphas in same chain, including insertion mutants (i.e. same resi but different inscode)
            others = self._df.loc[(self._df['an']=='CA') & (-(self._df['resi']==resi) | -(self._df['inscode']==inscode)) & (self._df['chain']==chain)].index
            first = second = None            
            # copy distance values from matrix and sort
            distances = self._distMatrix.loc[idx, others].copy()
            distances.sort()
            # first two entries
            minEntries = self._df.loc[distances.iloc[:2].index, 'resn']
            first, second = minEntries

            bs = np.NaN
            # set the bits            
            if resn in PPIBitstrings.AADICT:
                bs = np.zeros(60, dtype=np.int)
                bs[PPIBitstrings.AADICT[resn][0]] = np.int(1)
                for k,v in enumerate([first, second]):
                    if v in PPIBitstrings.AADICT:
                        bs[(k+1)*20 + PPIBitstrings.AADICT[v][0]] = np.int(1)   
            bitstrings.append(bs)
        
        # create series
        s = pd.Series(bitstrings, index=cas.index, name='bitstring', dtype=np.object)
        # join to our df
        self._df = self._df.join(pd.DataFrame(s))       

def qd_read_picks_from_hyp_file(filename):
    f = open(filename, "r")
    lines = f.readlines()
    f.close()

    for iline in range(len(lines)):
        line = lines[iline]
        words = line.split()
        if words[0] == "PHASE":
            iline_phase = iline
            break

    phases = {}
    for line in lines[iline + 1 :]:
        words = line.split()
        try:
            if words[4] == "P":
                station = words[0]
                year = np.int(words[6][0:4])
                month = np.int(words[6][4:6])
                day = np.int(words[6][6:8])
                hour = np.int(words[7][0:2])
                minute = np.int(words[7][2:4])
                seconds = np.float(words[8])
                ptime = utcdatetime.UTCDateTime(year, month, day, hour, minute, seconds)
                phases[station] = ptime
        except IndexError:
            pass

    return phases

    def testInt(self):
        num = np.int(2562010)
        self.assertEqual(np.int(ujson.decode(ujson.encode(num))), num)

        num = np.int8(127)
        self.assertEqual(np.int8(ujson.decode(ujson.encode(num))), num)

        num = np.int16(2562010)
        self.assertEqual(np.int16(ujson.decode(ujson.encode(num))), num)

        num = np.int32(2562010)
        self.assertEqual(np.int32(ujson.decode(ujson.encode(num))), num)

        num = np.int64(2562010)
        self.assertEqual(np.int64(ujson.decode(ujson.encode(num))), num)

        num = np.uint8(255)
        self.assertEqual(np.uint8(ujson.decode(ujson.encode(num))), num)

        num = np.uint16(2562010)
        self.assertEqual(np.uint16(ujson.decode(ujson.encode(num))), num)

        num = np.uint32(2562010)
        self.assertEqual(np.uint32(ujson.decode(ujson.encode(num))), num)

        num = np.uint64(2562010)
        self.assertEqual(np.uint64(ujson.decode(ujson.encode(num))), num)

def dispims_color(M, border=0, bordercolor=[0.0, 0.0, 0.0], savePath=None, *imshow_args, **imshow_keyargs):
    """ Display an array of rgb images. 

    The input array is assumed to have the shape numimages x numpixelsY x numpixelsX x 3
    """
    bordercolor = numpy.array(bordercolor)[None, None, :]
    numimages = len(M)
    M = M.copy()
    for i in range(M.shape[0]):
        M[i] -= M[i].flatten().min()
        M[i] /= M[i].flatten().max()
    height, width, three = M[0].shape
    assert three == 3
    
    n0 = numpy.int(numpy.ceil(numpy.sqrt(numimages)))
    n1 = numpy.int(numpy.ceil(numpy.sqrt(numimages)))
    im = numpy.array(bordercolor)*numpy.ones(
                             ((height+border)*n1+border,(width+border)*n0+border, 1),dtype='<f8')
    for i in range(n0):
        for j in range(n1):
            if i*n1+j < numimages:
                im[j*(height+border)+border:(j+1)*(height+border)+border,
                   i*(width+border)+border:(i+1)*(width+border)+border,:] = numpy.concatenate((
                  numpy.concatenate((M[i*n1+j,:,:,:],
                         bordercolor*numpy.ones((height,border,3),dtype=float)), 1),
                  bordercolor*numpy.ones((border,width+border,3),dtype=float)
                  ), 0)
    imshow_keyargs["interpolation"]="nearest"
    pylab.imshow(im, *imshow_args, **imshow_keyargs)
    
    if savePath == None:
        pylab.show()
    else:
        pylab.savefig(savePath)

    def _subplot_dims(self):
        """Determine size of subplot grid, given possible
        constraints on nrows, ncols"""
        nrows = self.subplot_opts.get('nrows', None)
        ncols = self.subplot_opts.get('ncols', None)

        #if 2 keys provided, rows and cols are fixed
        if len(self._keys) == 2:
            nr = len(self._key_index[0])
            nc = len(self._key_index[1])

            if ((nrows is not None and nrows != nr) or
                (ncols is not None and ncols != nc)):
                raise ValueError("Two keys specified: (nrows, ncols) must be "
                                 "(%i, %i) " % (nr, nc))
            return nr, nc

        sz = len(self._key_index[0])

        #if 1 key provided, just need nrows * ncols >= nfacets
        if nrows is None:
            if ncols is None:
                nrows = max(1, np.int(np.sqrt(sz)))
            else:
                nrows = np.int(np.ceil(1. * sz / ncols))
        if ncols is None:
            ncols = np.int(np.ceil(1. * sz / nrows))
        if nrows * ncols < sz:
            raise ValueError("nrows (%i) and ncols (%i) not big enough "
                             "to plot %i facets" % (nrows, ncols, sz))
        return nrows, ncols

    def _grid(self, corners=False):
        """Create an xy grid of coordinates for heliographic array.

        Uses meshgrid. If corners is selected, this function will shift the array by half a pixel in both directions
        so that the corners of the normal array can be accessed easily.

        Args:
            corners (bool, optional): defaults to False, chooses whether to apply the corner calculation or not

        Returns:
            xg: 2D array containing the x-coordinates of each pixel
            yg: 2D array containing the y-coordinates of each pixel

        """
        # Retrieve integer dimensions and create arrays holding
        # x and y coordinates of each pixel
        x_dim = np.int(np.floor(self.im_raw.dimensions[0].value))
        y_dim = np.int(np.floor(self.im_raw.dimensions[1].value))

        if corners:
            x_row = (np.arange(0, x_dim + 1) - self.par['X0'] - 0.5) * self.par['xscale']
            y_row = (np.arange(0, y_dim + 1) - self.par['Y0'] - 0.5) * self.par['yscale']
            xg, yg = mnp.meshgrid(x_row, y_row)
            rg = mnp.sqrt(xg ** 2 + yg ** 2)
            self.Rg = rg
        else:
            x_row = (np.arange(0, x_dim) - self.par['X0']) * self.par['xscale']
            y_row = (np.arange(0, y_dim) - self.par['Y0']) * self.par['yscale']
            xg, yg = mnp.meshgrid(x_row, y_row)
            rg = mnp.sqrt(xg ** 2 + yg ** 2)
            self.xg = xg
            self.yg = yg
            self.rg = rg

        return xg, yg

def avhrr(scans_nb, scan_points,
          scan_angle=55.37, frequency=1 / 6.0, apply_offset=True):
    """Definition of the avhrr instrument.

    Source: NOAA KLM User's Guide, Appendix J
    http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/j/app-j.htm
    """
    # build the avhrr instrument (scan angles)
    avhrr_inst = np.vstack(((scan_points / 1023.5 - 1)
                            * np.deg2rad(-scan_angle),
                            np.zeros((len(scan_points),))))

    avhrr_inst = np.tile(
        avhrr_inst[:, np.newaxis, :], [1, np.int(scans_nb), 1])

    # building the corresponding times array
    # times = (np.tile(scan_points * 0.000025 + 0.0025415, [scans_nb, 1])
    #         + np.expand_dims(offset, 1))

    times = np.tile(scan_points * 0.000025, [np.int(scans_nb), 1])
    if apply_offset:
        offset = np.arange(np.int(scans_nb)) * frequency
        times += np.expand_dims(offset, 1)

    return ScanGeometry(avhrr_inst, times)

	def get_2state_gaussian_seq(lens,dims=2,means1=[2,2,2,2],means2=[5,5,5,5],vars1=[1,1,1,1],vars2=[1,1,1,1],anom_prob=1.0):
		
		seqs = co.matrix(0.0, (dims, lens))
		lbls = co.matrix(0, (1,lens))
		marker = 0

		# generate first state sequence
		for d in range(dims):
			seqs[d,:] = co.normal(1,lens)*vars1[d] + means1[d]

		prob = np.random.uniform()
		if prob<anom_prob:		
			# add second state blocks
			while (True):
				max_block_len = 0.6*lens
				min_block_len = 0.1*lens
				block_len = np.int(max_block_len*np.single(co.uniform(1))+3)
				block_start = np.int(lens*np.single(co.uniform(1)))

				if (block_len - (block_start+block_len-lens)-3>min_block_len):
					break

			block_len = min(block_len,block_len - (block_start+block_len-lens)-3)
			lbls[block_start:block_start+block_len-1] = 1
			marker = 1
			for d in range(dims):
				#print block_len
				seqs[d,block_start:block_start+block_len-1] = co.normal(1,block_len-1)*vars2[d] + means2[d]

		return (seqs, lbls, marker)

def olci(scans_nb, scan_points=None):
    """Definition of the OLCI instrument.

    Source: Sentinel-3 OLCI Coverage
    https://sentinel.esa.int/web/sentinel/user-guides/sentinel-3-olci/coverage
    """

    if scan_points is None:
        scan_len = 4000  # samples per scan
        scan_points = np.arange(4000)
    else:
        scan_len = len(scan_points)
    # scan_rate = 0.044  # single scan, seconds
    scan_angle_west = 46.5  # swath, degrees
    scan_angle_east = -22.1  # swath, degrees
    # sampling_interval = 18e-3  # single view, seconds
    # build the olci instrument scan line angles
    scanline_angles = np.linspace(np.deg2rad(scan_angle_west),
                                  np.deg2rad(scan_angle_east), scan_len)
    inst = np.vstack((scanline_angles, np.zeros(scan_len,)))

    inst = np.tile(inst[:, np.newaxis, :], [1, np.int(scans_nb), 1])

    # building the corresponding times array
    # times = (np.tile(scan_points * 0.000025 + 0.0025415, [scans_nb, 1])
    #         + np.expand_dims(offset, 1))

    times = np.tile(np.zeros_like(scanline_angles), [np.int(scans_nb), 1])
    # if apply_offset:
    #     offset = np.arange(np.int(scans_nb)) * frequency
    #     times += np.expand_dims(offset, 1)

    return ScanGeometry(inst, times)

def run():
    def predict_return(N,d,Pfade):
        anz_zurueck = 0
        for _ in range(Pfade):
            walks = rw(N,d)
            weite = np.zeros(d)
            for i in range(N):
                weite = weite + walks[i]
                if np.sum(np.abs(weite)) == 0:
                    anz_zurueck += 1
                    break
        print np.float(anz_zurueck)/np.float(Pfade) 
        return anz_zurueck/np.float(Pfade) 
    print '1D'
    d1 = [predict_return(np.int(np.exp(i)),1,200) for i in range(15)]
    #print d1
    print '2D'
    d2 = [predict_return(np.int(np.exp(i)),2,200) for i in range(15)]
    print '3D'
    d3 = [predict_return(np.int(np.exp(i)),3,200) for i in range(15)]
    plt.plot(range(15),d1,label = 'Prediction in 1D')
    plt.plot(range(15),d2,label = 'Prediction in 2D')
    plt.plot(range(15),d3,label = 'Prediction in 3D')
    plt.legend()
    plt.show()

def test_superpixel(aia171_test_map, aia171_test_map_with_mask):
    dimensions = (2, 2)*u.pix
    superpixel_map_sum = aia171_test_map.superpixel(dimensions)
    assert_quantity_allclose(superpixel_map_sum.dimensions[1], aia171_test_map.dimensions[1]/dimensions[1]*u.pix)
    assert_quantity_allclose(superpixel_map_sum.dimensions[0], aia171_test_map.dimensions[0]/dimensions[0]*u.pix)
    assert_quantity_allclose(superpixel_map_sum.data[0][0], (aia171_test_map.data[0][0] +
                                                             aia171_test_map.data[0][1] +
                                                             aia171_test_map.data[1][0] +
                                                             aia171_test_map.data[1][1]))

    superpixel_map_avg = aia171_test_map.superpixel(dimensions, func=np.mean)
    assert_quantity_allclose(superpixel_map_avg.dimensions[1], aia171_test_map.dimensions[1]/dimensions[1]*u.pix)
    assert_quantity_allclose(superpixel_map_avg.dimensions[0], aia171_test_map.dimensions[0]/dimensions[0]*u.pix)
    assert_quantity_allclose(superpixel_map_avg.data[0][0], (aia171_test_map.data[0][0] +
                                                             aia171_test_map.data[0][1] +
                                                             aia171_test_map.data[1][0] +
                                                             aia171_test_map.data[1][1])/4.0)

    # Test that the mask is respected
    superpixel_map_sum = aia171_test_map_with_mask.superpixel(dimensions)
    assert superpixel_map_sum.mask is not None
    assert_quantity_allclose(superpixel_map_sum.mask.shape[0],
                             aia171_test_map.dimensions[1]/dimensions[1])
    assert_quantity_allclose(superpixel_map_sum.mask.shape[1],
                             aia171_test_map.dimensions[0]/dimensions[0])

    # Test that the offset is respected
    superpixel_map_sum = aia171_test_map_with_mask.superpixel(dimensions, offset=(1, 1)*u.pix)
    assert_quantity_allclose(superpixel_map_sum.dimensions[1], aia171_test_map.dimensions[1]/dimensions[1]*u.pix - 1*u.pix)
    assert_quantity_allclose(superpixel_map_sum.dimensions[0], aia171_test_map.dimensions[0]/dimensions[0]*u.pix - 1*u.pix)

    dimensions = (7, 9)*u.pix
    superpixel_map_sum = aia171_test_map_with_mask.superpixel(dimensions, offset=(4, 4)*u.pix)
    assert_quantity_allclose(superpixel_map_sum.dimensions[0], np.int((aia171_test_map.dimensions[0]/dimensions[0]).value)*u.pix - 1*u.pix)
    assert_quantity_allclose(superpixel_map_sum.dimensions[1], np.int((aia171_test_map.dimensions[1]/dimensions[1]).value)*u.pix - 1*u.pix)

def _interpolate_2d(x, y, a1, b1, a2, b2, c):
  
  """
  parameters:
  (x,y) = where we want to estimate the function 
  a1 , b1 = lower bounds of the grid
  a2 , b2 = upper bounds of the grid
  c = spline coefficients
  """
  n1 = c.shape[0] - 3
  n2 = c.shape[1] - 3
  h1 = (b1 - a1)/n1
  h2 = (b2 - a2)/n2
  
  l1 = np.int((x - a1)/h1) + 1
  l2 = np.int((y - a2)/h2) + 1
  m1 = min(l1 + 3, n1 + 3)
  m2 = min(l2 + 3, n2 + 3)
  
  s = 0

  for i1 in xrange(l1, m1 + 1):
    u_x = u(x, i1, a1, h1)
    for i2 in xrange(l2, m2 + 1):
      u_y = u(y, i2, a2, h2)
      s += c[i1 - 1, i2 - 1] * u_x * u_y

  return s

def rand_jacobi_rotation(A):
    """Random Jacobi rotation of a sparse matrix.
    
    Parameters
    ----------
    A : spmatrix
        Input sparse matrix.
    
    Returns
    -------
    spmatrix
        Rotated sparse matrix.
    """
    if A.shape[0] != A.shape[1]:
        raise Exception("Input matrix must be square.")
    n = A.shape[0]
    angle = (2 * np.random.random() - 1) * np.pi
    a = 1.0 / np.sqrt(2) * np.exp(-1j * angle)
    b = 1.0 / np.sqrt(2) * np.exp(1j * angle)
    i = np.int(np.floor(np.random.random() * n))
    j = i
    while i == j:
        j = np.int(np.floor(np.random.random() * n))
    data = np.hstack(([a, -b, a, b], np.ones(n - 2, dtype=int)))
    diag = np.delete(np.arange(n), [i, j])
    rows = np.hstack(([i, i, j, j], diag))
    cols = np.hstack(([i, j, i, j], diag))
    R = sp.coo_matrix((data, (rows, cols)), shape=[n, n]).tocsr()
    A = R * A * R.conj().transpose()
    return A

def findValidFFTWDim( inputDims ):
    """
    Finds a valid dimension for which FFTW can optimize its calculations. The 
    return is a shape which is forced to be square, as this gives uniform pixel
    size in x-y in Fourier space.
    
    If you want a minimum padding size, call as findValidFFTWDim( image.shape + 128 ) 
    or similar.
    """
    dim = np.max( np.round( inputDims ) )
    maxPow2 = np.int( np.ceil( math.log( dim, 2 ) ) )
    maxPow3 = np.int( np.ceil( math.log( dim, 3 ) ) )
    maxPow5 = np.int( np.ceil( math.log( dim, 5 ) ) )
    maxPow7 = np.int( np.ceil( math.log( dim, 7 ) ) )   
    
    dimList = np.zeros( [(maxPow2+1)*(maxPow3+1)*(maxPow5+1)*(maxPow7+1)] )
    count = 0
    for I in np.arange(0,maxPow7+1):
        for J in np.arange(0,maxPow5+1):
            for K in np.arange(0,maxPow3+1):
                for L in np.arange(0,maxPow2+1):
                    dimList[count] = 2**L * 3**K * 5**J * 7**I
                    count += 1
    dimList = np.sort( np.unique( dimList ) )
    dimList = dimList[ np.argwhere(dimList < 2*dim)].squeeze()
    dimList = dimList.astype('int64')
    # Throw out odd image shapes, this just causes more problems with many 
    # functions
    dimList = dimList[ np.mod(dimList,2)==0 ]
    
    # Find first dim that equals or exceeds dim
    nextValidDim =  dimList[np.argwhere( dimList >= dim)[0,0]]
    return np.array( [nextValidDim, nextValidDim] )    

def chain2image(chaincode,start_pix):

    """
    Method to compute the pixel contour providing the chain code string
    and the starting pixel location [X,Y].
    Author: Xavier Bonnin (LESIA)
    """

    if (type(chaincode) != str):
        print "First input argument must be a string!"
        return None

    if (len(start_pix) != 2):
        print "Second input argument must be a 2-elements vector!"
        return None

    ardir = np.array([[-1,0],[-1,1],[0,1],[1,1],[1,0],[1,-1],[0,-1],[-1,-1]])
    ccdir = np.array([0,7,6,5,4,3,2,1])

    X=[start_pix[0]]
    Y=[start_pix[1]]
    for c in chaincode:
        if (abs(np.int8(c)) > 7):
            print "Wrong chain code format!"
            return None
        wc = np.where(np.int8(c) == np.int8(ccdir))[0]
        X.append(X[-1] + np.int(ardir[wc,0]))
        Y.append(Y[-1] + np.int(ardir[wc,1]))
    return X,Y

    def testIntMax(self):
        num = np.int(np.iinfo(np.int).max)
        self.assertEqual(np.int(ujson.decode(ujson.encode(num))), num)

        num = np.int8(np.iinfo(np.int8).max)
        self.assertEqual(np.int8(ujson.decode(ujson.encode(num))), num)

        num = np.int16(np.iinfo(np.int16).max)
        self.assertEqual(np.int16(ujson.decode(ujson.encode(num))), num)

        num = np.int32(np.iinfo(np.int32).max)
        self.assertEqual(np.int32(ujson.decode(ujson.encode(num))), num)

        num = np.uint8(np.iinfo(np.uint8).max)
        self.assertEqual(np.uint8(ujson.decode(ujson.encode(num))), num)

        num = np.uint16(np.iinfo(np.uint16).max)
        self.assertEqual(np.uint16(ujson.decode(ujson.encode(num))), num)

        num = np.uint32(np.iinfo(np.uint32).max)
        self.assertEqual(np.uint32(ujson.decode(ujson.encode(num))), num)

        if platform.architecture()[0] != '32bit':
            num = np.int64(np.iinfo(np.int64).max)
            self.assertEqual(np.int64(ujson.decode(ujson.encode(num))), num)

            # uint64 max will always overflow as it's encoded to signed
            num = np.uint64(np.iinfo(np.int64).max)
            self.assertEqual(np.uint64(ujson.decode(ujson.encode(num))), num)

def loss(x, method, a_vec=np.zeros(3)):
    rtn = 0
    if(sum(a_vec)==0):
        tmp = sorted(x)
        a_vec[0] = tmp[np.int(len(tmp)*0.5)]
        a_vec[1] = tmp[np.int(len(tmp)*0.75)]
        a_vec[2] = tmp[np.int(len(tmp)*0.85)]
    a = a_vec[0]
    b = a_vec[1]
    c = a_vec[2]
    if(sum(x<0)==0):
        rtn = np.zeros(len(x))
        
        rtn[x<=a] = x[x<=a]**2/2
        rtn[(x>a)*(x<=b)] = a*x[(x>a)*(x<=b)]-a*a/2
        rtn[(x>b)*(x<=c)] = a*(x[(x>b)*(x<=c)]-c)**2/(2*(b-c))+a*(b+c-a)/2
        rtn[x>c] = a*(b+c-a)/2
        
        '''
        tmp = x[x<=c]/c
        rtn[x<=c] = 1-(1-tmp**2)**3
        rtn[x>c] = 1
        '''
        #rtn = a**2*np.log(1+(x/a)**2)
        #rtn = a**2*(np.sqrt(1+(x/a)**2)-1)
    return(rtn)