def row_func(i):
     pyrandom.seed(initseed + int(i))
     scirandom.seed(initseed + int(i))
     k = scirandom.binomial(N, p, 1)[0]
     cur_row[:] = 0.0
     cur_row[pyrandom.sample(myrange, k)] = weight
     return cur_row

    def test_gmres(self):
        # Ensure repeatability
        random.seed(0)
        
        #  For these small matrices, Householder and MGS GMRES should give the same result, 
        #  and for symmetric (but possibly indefinite) matrices CR and GMRES should give same result
        for maxiter in [1,2,3]:
            for case, symm_case in zip(self.cases, self.symm_cases):
                A = case['A'] 
                b = case['b']
                x0 = case['x0']
                A_symm = symm_case['A'] 
                b_symm = symm_case['b']
                x0_symm = symm_case['x0']
                
                # Test agreement between Householder and GMRES
                (x, flag) = gmres_householder(A,b,x0=x0,maxiter=min(A.shape[0],maxiter))
                (x2, flag2) = gmres_mgs(A,b,x0=x0,maxiter=min(A.shape[0],maxiter))
                assert_array_almost_equal(x/norm(x), x2/norm(x2), 
                        err_msg='Householder GMRES and MGS GMRES gave different results for small matrix')
                assert_equal(flag, flag2, 
                        err_msg='Householder GMRES and MGS GMRES returned different convergence flags for small matrix')

                # Test agreement between GMRES and CR
                if A_symm.shape[0] > 1:
                    residuals2 = []
                    (x2, flag2) = gmres_mgs(A_symm, b_symm, x0=x0_symm, maxiter=min(A.shape[0],maxiter),residuals=residuals2)
                    residuals3 = []
                    (x3, flag2) = cr(A_symm, b_symm, x0=x0_symm, maxiter=min(A.shape[0],maxiter),residuals=residuals3)
                    residuals2 = array(residuals2)
                    residuals3 = array(residuals3)
                    assert_array_almost_equal(residuals3/norm(residuals3), residuals2/norm(residuals2), 
                            err_msg='CR and GMRES yield different residual vectors')
                    assert_array_almost_equal(x2/norm(x2), x3/norm(x3), err_msg='CR and GMRES yield different answers')

  def reset(self, params, repetition):
    # if params['encoding'] == 'basic':
    #   self.inputEncoder = PassThroughEncoder()
    # elif params['encoding'] == 'distributed':
    #   self.outputEncoder = PassThroughEncoder()
    # else:
    #   raise Exception("Encoder not found")

    print params
    self.inputEncoder = PassThroughEncoder()
    self.outputEncoder = PassThroughEncoder()

    if params['dataset'] == 'nyc_taxi':
      self.dataset = NYCTaxiDataset()
    else:
      raise Exception("Dataset not found")

    self.testCounter = 0

    self.history = []
    self.resets = []
    self.currentSequence = self.dataset.generateSequence()

    random.seed(6)
    self.nDimInput = 3
    self.nDimOutput = 1
    self.net = buildNetwork(self.nDimInput, params['num_cells'], self.nDimOutput,
                       hiddenclass=LSTMLayer, bias=True, outputbias=True, recurrent=True)

  def reset(self, params, repetition):
    print params

    self.nDimInput = 3
    self.inputEncoder = PassThroughEncoder()

    if params['output_encoding'] == None:
      self.outputEncoder = PassThroughEncoder()
      self.nDimOutput = 1
    elif params['output_encoding'] == 'likelihood':
      self.outputEncoder = ScalarBucketEncoder()
      self.nDimOutput = self.outputEncoder.encoder.n

    if params['dataset'] == 'nyc_taxi' or params['dataset'] == 'nyc_taxi_perturb_baseline':
      self.dataset = NYCTaxiDataset(params['dataset'])
    else:
      raise Exception("Dataset not found")

    self.testCounter = 0
    self.resets = []
    self.iteration = 0

    # initialize LSTM network
    random.seed(6)
    if params['output_encoding'] == None:
      self.net = buildNetwork(self.nDimInput, params['num_cells'], self.nDimOutput,
                         hiddenclass=LSTMLayer, bias=True, outputbias=True, recurrent=True)
    elif params['output_encoding'] == 'likelihood':
      self.net = buildNetwork(self.nDimInput, params['num_cells'], self.nDimOutput,
                         hiddenclass=LSTMLayer, bias=True, outclass=SigmoidLayer, recurrent=True)

    (self.networkInput, self.targetPrediction, self.trueData) = \
      self.dataset.generateSequence(
      prediction_nstep=params['prediction_nstep'],
      output_encoding=params['output_encoding'])

 def setUp(self):
     self.numberOfSamples = 1000
     pkl_file = open(os.path.join(unittest_dir, "test_data", "kdeOrigin_xyz.pck"), "r")
     xp = cPickle.load(pkl_file)
     yp = cPickle.load(pkl_file)
     zp = cPickle.load(pkl_file)
     self.sampOrg = SamplingOrigin.SamplingOrigin(zp, xp, yp)
     random.seed(10)

def setRandomParameters(net,seed=None,randFunc=random.random):
    """
    Sets parameters to random values given by the function randFunc (by
    default, uniformly distributed on [0,1) ).
    """
    random.seed(seed)
    net.setOptimizables( randFunc(len(net.GetParameters())) )
    return net.GetParameters()

def GaussianRandomInitializer(gridShape, sigma=0.2, seed=None, slipSystem=None, slipPlanes=None, slipDirections=None, vacancy=None, smectic=None):

    oldgrid = copy.copy(gridShape)
   
    if len(gridShape) == 1:
	    gridShape = (128,)
    if len(gridShape) == 2:
	    gridShape = (128,128)
    if len(gridShape) == 3:
	    gridShape = (128,128,128)

    """ Returns a random initial set of fields of class type PlasticityState """
    if slipSystem=='gamma':
        state = SlipSystemState.SlipSystemState(gridShape,slipPlanes=slipPlanes,slipDirections=slipDirections)
    elif slipSystem=='betaP':
        state = SlipSystemBetaPState.SlipSystemState(gridShape,slipPlanes=slipPlanes,slipDirections=slipDirections)
    else:
        if vacancy is not None:
            state = VacancyState.VacancyState(gridShape,alpha=vacancy)
        elif smectic is not None:
            state = SmecticState.SmecticState(gridShape)
        else:
            state = PlasticityState.PlasticityState(gridShape)

    field = state.GetOrderParameterField()
    Ksq_prime = FourierSpaceTools.FourierSpaceTools(gridShape).kSq * (-sigma**2/4.)

    if seed is None:
        seed = 0
    n = 0
    random.seed(seed)

    Ksq = FourierSpaceTools.FourierSpaceTools(gridShape).kSq.numpy_array()

    for component in field.components:
        temp = random.normal(scale=gridShape[0],size=gridShape)
        ktemp = fft.rfftn(temp)*(sqrt(pi)*sigma)**len(gridShape)*exp(-Ksq*sigma**2/4.)
        field[component] = numpy.real(fft.irfftn(ktemp))
        #field[component] = GenerateGaussianRandomArray(gridShape, temp ,sigma)
        n += 1

    """
    t, s = LoadState("2dstate32.save", 0)
    for component in field.components:
        for j in range(0,32):
            field[component][:,:,j] = s.betaP[component].numpy_array()
    """

    ## To make seed consistent across grid sizes and convergence comparison
    gridShape = copy.copy(oldgrid)
    if gridShape[0] != 128:
        state = ResizeState(state,gridShape[0],Dim=len(gridShape))

    state = ReformatState(state)
    state.ktools = FourierSpaceTools.FourierSpaceTools(gridShape)
    
    return state 

    def test_minimal_residual(self):
        # Ensure repeatability
        random.seed(0)

        self.definite_cases.extend(self.spd_cases)

        for case in self.definite_cases:
            A = case['A']
            maxiter = case['maxiter']
            x0 = rand(A.shape[0],)
            b = zeros_like(x0)
            reduction_factor = case['reduction_factor']
            if A.dtype != complex:

                # This function should always decrease (assuming zero RHS)
                fvals = []

                def callback(x):
                    fvals.append(sqrt(dot(ravel(x),
                                 ravel(A*x.reshape(-1, 1)))))
                #
                (x, flag) = minimal_residual(A, b, x0=x0,
                                             tol=1e-16, maxiter=maxiter,
                                             callback=callback)
                actual_factor = (norm(ravel(b) - ravel(A*x.reshape(-1, 1))) /
                                 norm(ravel(b) - ravel(A*x0.reshape(-1, 1))))
                assert(actual_factor < reduction_factor)
                if A.dtype != complex:
                    for i in range(len(fvals)-1):
                        assert(fvals[i+1] <= fvals[i])

        # Test preconditioning
        A = pyamg.gallery.poisson((10, 10), format='csr')
        x0 = rand(A.shape[0], 1)
        b = zeros_like(x0)
        fvals = []

        def callback(x):
            fvals.append(sqrt(dot(ravel(x), ravel(A*x.reshape(-1, 1)))))
        #
        resvec = []
        sa = pyamg.smoothed_aggregation_solver(A)
        (x, flag) = minimal_residual(A, b, x0, tol=1e-8, maxiter=20,
                                     residuals=resvec, M=sa.aspreconditioner(),
                                     callback=callback)
        assert(resvec[-1] < 1e-8)
        for i in range(len(fvals)-1):
            assert(fvals[i+1] <= fvals[i])

    def test_steepest_descent(self):
        # Ensure repeatability
        random.seed(0)

        for case in self.spd_cases:
            A = case['A']
            b = case['b']
            x0 = case['x0']
            maxiter = case['maxiter']
            reduction_factor = case['reduction_factor']

            # This function should always decrease
            fvals = []

            def callback(x):
                fvals.append(0.5*dot(ravel(x), ravel(A*x.reshape(-1, 1))) -
                             dot(ravel(b), ravel(x)))

            (x, flag) = steepest_descent(A, b, x0=x0, tol=1e-16,
                                         maxiter=maxiter, callback=callback)
            actual_factor = (norm(ravel(b) - ravel(A*x.reshape(-1, 1))) /
                             norm(ravel(b) - ravel(A*x0.reshape(-1, 1))))
            assert(actual_factor < reduction_factor)

            if A.dtype != complex:
                for i in range(len(fvals)-1):
                    assert(fvals[i+1] <= fvals[i])

        # Test preconditioning
        A = pyamg.gallery.poisson((10, 10), format='csr')
        b = rand(A.shape[0], 1)
        x0 = rand(A.shape[0], 1)
        fvals = []

        def callback(x):
            fvals.append(0.5*dot(ravel(x), ravel(A*x.reshape(-1, 1))) -
                         dot(ravel(b), ravel(x)))

        resvec = []
        sa = pyamg.smoothed_aggregation_solver(A)
        (x, flag) = steepest_descent(A, b, x0, tol=1e-8, maxiter=20,
                                     residuals=resvec, M=sa.aspreconditioner(),
                                     callback=callback)
        assert(resvec[-1] < 1e-8)
        for i in range(len(fvals)-1):
            assert(fvals[i+1] <= fvals[i])

    def test_ishermitian(self):
        # make tests repeatable
        random.seed(0)
        casesT = []
        casesF = []
        # 1x1
        casesT.append(mat(rand(1, 1)))
        casesF.append(mat(1.0j*rand(1, 1)))
        # 2x2
        A = array([[1.0, 0.0], [2.0, 1.0]])
        Ai = 1.0j*A
        casesF.append(A)
        casesF.append(Ai)
        A = A + Ai
        casesF.append(A)
        casesT.append(A + A.conjugate().T)
        # 3x3
        A = mat(rand(3, 3))
        Ai = 1.0j*rand(3, 3)
        casesF.append(A)
        casesF.append(Ai)
        A = A + Ai
        casesF.append(A)
        casesT.append(A + A.H)

        for A in casesT:
            # dense arrays
            assert_equal(ishermitian(A, fast_check=False), True)
            assert_equal(ishermitian(A, fast_check=True), True)

            # csr arrays
            A = csr_matrix(A)
            assert_equal(ishermitian(A, fast_check=False), True)
            assert_equal(ishermitian(A, fast_check=True), True)

        for A in casesF:
            # dense arrays
            assert_equal(ishermitian(A, fast_check=False), False)
            assert_equal(ishermitian(A, fast_check=True), False)

            # csr arrays
            A = csr_matrix(A)
            assert_equal(ishermitian(A, fast_check=False), False)
            assert_equal(ishermitian(A, fast_check=True), False)

    def get_neural_net(input_number, output_number, NetworkType, HiddenLayerType, neurons_per_layer=[9], use_bias=False):
        random.seed(123)

        input = LinearLayer(input_number)
        output = SoftmaxLayer(output_number)

        neural_net = NetworkType()
        neural_net.addInputModule(input)
        neural_net.addOutputModule(output)

        if use_bias:
            bias = BiasUnit()
            neural_net.addModule(bias)

        prev = input
        for i in range(0, len(neurons_per_layer)):
            hidden = HiddenLayerType(neurons_per_layer[i])

            neural_net.addModule(hidden)
            neural_net.addConnection(FullConnection(prev, hidden))

            if use_bias:
                neural_net.addConnection(FullConnection(bias, hidden))

            prev = hidden

        neural_net.addConnection(FullConnection(prev, output))

        if use_bias:
            neural_net.addConnection(FullConnection(bias, output))

        neural_net.sortModules()

        fast_net = neural_net.convertToFastNetwork()

        if fast_net is not None:
            neural_net = fast_net
            print "Use fast C++ implementation"
        else:
            print "Use standard Python implementation"

        print "Create neural network with {} neurons ({} layers)".format(neurons_per_layer, len(neurons_per_layer))

        return neural_net

def SineWaveInitializer(gridShape, randomPhase=False, seed=None, slipSystem=False, slipPlanes=None, slipDirections=None):
    """
    Initialize a plasticity state by setting all its components in any dimension with a
    single period of a sine function.
    """
    if seed is not None:
        random.seed(seed)
    if slipSystem=='gamma':
        pass
#state = SlipSystemState.SlipSystemState(gridShape,slipPlanes=slipPlanes,slipDirections=slipDirections)
    elif slipSystem=='betaP':
        pass
#state = SlipSystemState.SlipSystemBetaPState(gridShape,slipPlanes=slipPlanes,slipDirections=slipDirections)
    else:
        state = PlasticityState.PlasticityState(gridShape)
    field = state.GetOrderParameterField()
    for component in field.components:
        field[component] = GenerateSineArray(gridShape, field.GridDimension(), randomPhase = randomPhase)
    return state 

  def train(self, params, verbose=False):

    if params['reset_every_training']:
      if verbose:
        print 'create lstm network'

      random.seed(6)
      if params['output_encoding'] == None:
        self.net = buildNetwork(self.nDimInput, params['num_cells'], self.nDimOutput,
                           hiddenclass=LSTMLayer, bias=True, outputbias=True, recurrent=True)
      elif params['output_encoding'] == 'likelihood':
        self.net = buildNetwork(self.nDimInput, params['num_cells'], self.nDimOutput,
                           hiddenclass=LSTMLayer, bias=True, outclass=SigmoidLayer, recurrent=True)

    self.net.reset()

    ds = SequentialDataSet(self.nDimInput, self.nDimOutput)
    networkInput = self.window(self.networkInput, params)
    targetPrediction = self.window(self.targetPrediction, params)

    # prepare a training data-set using the history
    for i in xrange(len(networkInput)):
      ds.addSample(self.inputEncoder.encode(networkInput[i]),
                   self.outputEncoder.encode(targetPrediction[i]))

    if params['num_epochs'] > 1:
      trainer = RPropMinusTrainer(self.net, dataset=ds, verbose=verbose)

      if verbose:
        print " train LSTM on ", len(ds), " records for ", params['num_epochs'], " epochs "

      if len(networkInput) > 1:
        trainer.trainEpochs(params['num_epochs'])

    else:
      self.trainer.setData(ds)
      self.trainer.train()

    # run through the training dataset to get the lstm network state right
    self.net.reset()
    for i in xrange(len(networkInput)):
      self.net.activate(ds.getSample(i)[0])

    def _set_params_and_init(self, mod, **kwargs):
        if not mod.converged: raise ValueError(" Error: This model has not converged, abort sampling ...")

        self.seed = 199
        self.verbose = False
        self.N = 1000
        if kwargs is not None:
            for key, value in kwargs.iteritems():
                setattr(self, key, value)

        random.seed(self.seed)

        if (mod.multi):
            #  N is recalculated in case of multi mass
            self.Nj = numpy.array(mod.Mj/mod.mj, dtype='int')
            self.N = sum(self.Nj)

            if min(self.Nj)==0: raise ValueError(" Error: One of the mass components has zero stars!")
        self.mod = mod
        self.ani = True if min(mod.raj)/mod.rt < 3 else False

    def test_condest(self):
        # make tests repeatable
        random.seed(0)

        # Should be exact for small matrices
        cases = []
        A = mat(array([2.14]))
        cases.append(A)
        A = mat(array([2.14j]))
        cases.append(A)
        A = mat(array([-1.2 + 2.14j]))
        cases.append(A)
        for i in range(1, 6):
            A = mat(rand(i, i))
            cases.append(A)
            cases.append(1.0j*A)
            A = mat(A + 1.0j*rand(i, i))
            cases.append(A)

        for A in cases:
            eigs = eigvals(A)
            exact = max(abs(eigs))/min(abs(eigs))
            c = condest(A)
            assert_almost_equal(exact, c)

def SmecticInitializer(gridShape, sigma=0.2, seed=None):
    if seed is None:
        seed = 0
    random.seed(seed)

    state = SmecticState.SmecticState(gridShape)
    field = state.GetOrderParameterField()

    Ksq = FourierSpaceTools.FourierSpaceTools(gridShape).kSq.numpy_array()

    for component in field.components:
        temp = random.normal(scale=gridShape[0],size=gridShape)
        ktemp = fft.rfftn(temp)*(sqrt(pi)*sigma)**len(gridShape)*exp(-Ksq*sigma**2/4.)
        field[component] = numpy.real(fft.irfftn(ktemp))

    ## To make seed consistent across grid sizes and convergence comparison
    gridShape = copy.copy(oldgrid)
    if gridShape[0] != 128:
        state = ResizeState(state,gridShape[0],Dim=len(gridShape))

    state = ReformatState(state)
    state.ktools = FourierSpaceTools.FourierSpaceTools(gridShape)
    
    return state 

    def test_cond(self):
        # make tests repeatable
        random.seed(0)

        # Should be exact
        cases = []
        A = mat(array([2.14]))
        cases.append(A)
        A = mat(array([2.14j]))
        cases.append(A)
        A = mat(array([-1.2 + 2.14j]))
        cases.append(A)
        for i in range(1, 6):
            A = mat(rand(i, i))
            cases.append(A)
            cases.append(1.0j*A)
            A = mat(A + 1.0j*rand(i, i))
            cases.append(A)

        for A in cases:
            U, Sigma, Vh = svd(A)
            exact = max(Sigma)/min(Sigma)
            c = cond(A)
            assert_almost_equal(exact, c)

def populate_geo_coord_polygon(polygon, number_of_points, seed=None,
                               exclude=None, randoms=None):
    """Populate a polygon with uniformly distributed points.
       Takes into account the curvature of the earth.

    Input:
       polygon - list of vertices of polygon, vertices in (lat, long)
       number_of_points - (optional) number of points
       seed - seed for random number generator (default=None)
       exclude - list of polygons (inside main polygon) from where points
       should be excluded

       randoms - array of numbers between 0 and 1 to be used as the
                 random numbers for calculating the latitude of the
                 points.  Obviously is this case the numbers are not
                 random.  This should just be used for writing tests.

    Output:
       points - list of points inside polygon

    Examples:
       populate_polygon( [[0,0], [-40,0], [-40,1], [0,1]], 6,
                         randoms = (0.0, 0.2, 0.4, 0.6, 0.8, 1.0))
       will return six not so randomly selected points inside the unit square
    """

    if seed is not None:
        random.seed(seed)

    if randoms is None:
        randoms = random.random(number_of_points)
    else:
        # Ensure randoms is a numpy array
        randoms = array(randoms)

    # print polygon
    points = []

    # Find outer extent of polygon
    max_x = min_x = polygon[0][1]
    max_y = min_y = polygon[0][0]
    for point in polygon[1:]:
        x = point[1]
        if x > max_x:
            max_x = x
        if x < min_x:
            min_x = x
        y = point[0]
        if y > max_y:
            max_y = y
        if y < min_y:
            min_y = y

    max_y_rad = radians(max_y)
    min_y_rad = radians(min_y)

    # Do this until we have a points array of points inside the polygon that
    # greater or equal the number of points requested (we'll slice at return)
    while len(points) < number_of_points:
        y_p_rad = arcsin(randoms * abs((sin(max_y_rad) - sin(min_y_rad))) +
                         sin(min_y_rad))
        y_p = degrees(y_p_rad)
        x_p = random.uniform(min_x, max_x, number_of_points)
        points_p = array([y_p, x_p]).T

        points_p = points_p[inside_polygon(points_p.tolist(), polygon)]
        if exclude is not None and len(points_p) > 0:
            for ex_poly in exclude:
                points_p = points_p[
                    outside_polygon(points_p.tolist(), ex_poly)]

        if len(points_p) > 0:
            points.extend(points_p.tolist())

        randoms = random.random(number_of_points)

    return points[:number_of_points]

import numpy as np
from scipy import random
from matplotlib import pyplot as plt

# Parameters
N = 200
A = 5
E = 1.0
L = 2.0
myseed = 42

# Initialize the random with seed
random.seed(myseed)

# Hidden Model
x = np.linspace(-L, L, N)
y1 = A * np.cos(0.5*np.pi*x/L)
y2 = E * random.normal(size=N)

# Show
plt.plot(x, y1+y2, 'bs', alpha=0.5, label="Medicion")
plt.plot(x, y1,    'k-', lw=2.0, label="Relacion determinista")
plt.xlim([-2.5,2.5])
plt.ylim([-5,10])
plt.xlabel("x []")
plt.ylabel("y []")
plt.legend(numpoints=1, loc="lower center")
#plt.savefig("images/dataN%d.png"%N)
plt.show()

# Shuffle the data

 def setUp(self):	
     random.seed(0) #make tests repeatable