def negSamplingCostAndGradient(predicted, target, outputVectors,
dataset, K = 10):
    klist = []
    for k in range(K):
        randomId = dataset.sampleTokenIdx()
        while randomId == target:
            randomId = dataset.sampleTokenIdx()
        klist.append(randomId)

    u0 = outputVectors[target]
    uks = -outputVectors[klist]
    vc = predicted
    U = np.vstack((u0, uks))

    dot = U.dot(vc)
    sigmoid_value = sigmoid(dot)
    cost = -np.sum(np.log(sigmoid_value))

    # print "cost is %f" % (cost, )

    gradPred = np.zeros(predicted.shape)
    grad = np.zeros(outputVectors.shape)

    temp = sigmoid(u0.dot(vc)) - 1
    intermediate = (sigmoid_value-1).reshape(-1, 1) * np.vstack((u0, -uks))
    gradPred += intermediate[0] - np.sum(intermediate[1:,], axis=0)
    grad[target] += temp * vc

    counter_dictionary = Counter(klist)
    unique_ks = list(counter_dictionary.keys())
    frequency_count = np.array(list(counter_dictionary.values()))
    grad[unique_ks] += (sigmoid(-outputVectors[unique_ks].dot(vc)) -1).reshape(-1,1) * -vc
    grad[unique_ks] *= frequency_count.reshape(-1, 1)

    return cost, gradPred, grad

def negSamplingCostAndGradient(predicted, target, outputVectors, dataset, 
    K=10):
    """ Negative sampling cost function for word2vec models """

    # Implement the cost and gradients for one predicted word vector  
    # and one target word vector as a building block for word2vec     
    # models, using the negative sampling technique. K is the sample  
    # size. You might want to use dataset.sampleTokenIdx() to sample  
    # a random word index. 
    # 
    # Note: See test_word2vec below for dataset's initialization.
    #                                       
    # Input/Output Specifications: same as softmaxCostAndGradient     
    # We will not provide starter code for this function, but feel    
    # free to reference the code you previously wrote for this        
    # assignment!
    
    ### YOUR CODE HERE
    grad = np.zeros(outputVectors.shape)
    gradPred = np.zeros(predicted.shape)

#    indices = [target]
#    for k in xrange(K):
#        newidx = dataset.sampleTokenIdx()
#        while newidx == target:
#            newidx = dataset.sampleTokenIdx()
#        indices += [newidx]
#
#    labels = np.array([1] + [-1 for k in xrange(K)])
#    vecs = outputVectors[indices,:]
#
#    t = sigmoid(vecs.dot(predicted) * labels)
#    cost = -np.sum(np.log(t))
#
#    delta = labels * (t - 1)
#    gradPred = delta.reshape((1,K+1)).dot(vecs).flatten()
#    gradtemp = delta.reshape((K+1,1)).dot(predicted.reshape(
#        (1,predicted.shape[0])))
#    for k in xrange(K+1):
#        grad[indices[k]] += gradtemp[k,:]

    t = sigmoid(predicted.dot(outputVectors[target,:]))
    cost = -np.log(t)
    delta = t - 1

    gradPred += delta * outputVectors[target, :]
    grad[target, :] += delta * predicted

    for k in xrange(K):
        idx = dataset.sampleTokenIdx()
    
        t = sigmoid(-predicted.dot(outputVectors[idx,:]))
        cost += -np.log(t)
        delta = 1 - t
    
        gradPred += delta * outputVectors[idx, :]
        grad[idx, :] += delta * predicted
    ### END YOUR CODE
    
    return cost, gradPred, grad

def negSamplingCostAndGradient(predicted, target, outputVectors, dataset,
    K=10):
    """ Negative sampling cost function for word2vec models """

    # Implement the cost and gradients for one predicted word vector
    # and one target word vector as a building block for word2vec
    # models, using the negative sampling technique. K is the sample
    # size. You might want to use dataset.sampleTokenIdx() to sample
    # a random word index.
    #
    # Note: See test_word2vec below for dataset's initialization.
    #
    # Input/Output Specifications: same as softmaxCostAndGradient
    # We will not provide starter code for this function, but feel
    # free to reference the code you previously wrote for this
    # assignment!

    ### YOUR CODE HERE
    grad = np.zeros(outputVectors.shape)

    s = sigmoid(np.dot(outputVectors[target,:], predicted))
    cost = -np.log(s)
    gradPred = - sigmoid_grad(s)/s*outputVectors[target,:]
    grad[target,:] = - sigmoid_grad(s)/s*predicted

    for k in range(K):
        i = dataset.sampleTokenIdx()
        s = sigmoid( - np.dot(outputVectors[i,:], predicted))
        cost -= np.log(s)
        gradPred += sigmoid_grad(s)/s*outputVectors[i,:]
        grad[i,:] += sigmoid_grad(s)/s*predicted

    ### END YOUR CODE

    return cost, gradPred, grad

def negSamplingCostAndGradient(predicted, target, outputVectors, dataset, 
    K=10):
    """ Negative sampling cost function for word2vec models """

    # Implement the cost and gradients for one predicted word vector  
    # and one target word vector as a building block for word2vec     
    # models, using the negative sampling technique. K is the sample  
    # size.
    #                                       
    # Input/Output Specifications: same as softmaxCostAndGradient     

    u_o = outputVectors[target,:]
    sigmoid_o = sigmoid(np.dot(u_o, predicted))
    cost = - np.log(sigmoid_o)
    gradPred = -u_o*(1-sigmoid_o)
    grad = np.zeros_like(outputVectors)
    grad[target,:] = - predicted*(1-sigmoid_o)
    
    for _ in range(K):
        k = dataset.sampleTokenIdx()
        while k == target:
            k = dataset.sampleTokenIdx()
        sigmoid_k = sigmoid(-np.dot(outputVectors[k,:],predicted))
        cost += - np.log(sigmoid_k)
        gradPred += outputVectors[k,:] * (1-sigmoid_k)
        grad[k,:] += predicted * (1-sigmoid_k)
    
    return cost, gradPred, grad

def negSamplingCostAndGradient(predicted, target, outputVectors, dataset, 
    K=10):
    """ Negative sampling cost function for word2vec models """

    # Implement the cost and gradients for one predicted word vector  
    # and one target word vector as a building block for word2vec     
    # models, using the negative sampling technique. K is the sample  
    # size. You might want to use dataset.sampleTokenIdx() to sample  
    # a random word index. 
    # 
    # Note: See test_word2vec below for dataset's initialization.
    #                                       
    # Input/Output Specifications: same as softmaxCostAndGradient     
    # We will not provide starter code for this function, but feel    
    # free to reference the code you previously wrote for this        
    # assignment!
    
    ### YOUR CODE HERE
    negativeSamples = [dataset.sampleTokenIdx() for i in range(K)]
    sigmoidTargetPred = sigmoid(outputVectors[target,:].transpose().dot(predicted))
    cost = -np.log(sigmoidTargetPred)   
    gradPred = (sigmoidTargetPred - 1.0)*outputVectors[target,:]
    grad = np.zeros(outputVectors.shape)
    grad[target,:] = predicted * (sigmoidTargetPred - 1.0)

    for sample in negativeSamples:
        sigmoidSamplePredicted = sigmoid(-outputVectors[sample,:].transpose().dot(predicted))
        cost -= np.log(sigmoidSamplePredicted)
        gradPred += (1.0 - sigmoidSamplePredicted)*outputVectors[sample,:]
        grad[sample,:] += (1.0 - sigmoidSamplePredicted)*predicted.transpose()
    ### END YOUR CODE
    
    return cost, gradPred, grad

def negSamplingCostAndGradient(predicted, target, outputVectors, dataset, 
    K=10):
    """ Negative sampling cost function for word2vec models """

    # Implement the cost and gradients for one predicted word vector  
    # and one target word vector as a building block for word2vec     
    # models, using the negative sampling technique. K is the sample  
    # size. You might want to use dataset.sampleTokenIdx() to sample  
    # a random word index. 
    # 
    # Note: See test_word2vec below for dataset's initialization.
    #                                       
    # Input/Output Specifications: same as softmaxCostAndGradient     
    # We will not provide starter code for this function, but feel    
    # free to reference the code you previously wrote for this        
    # assignment!
    
    ### YOUR CODE HERE
    '''
    Keep track of dims:
    
    D - dim of word vector
    V - number of words
    
    predicted     :  (D, )
    target        :  integer
    outputVectors :  (V, D)
    
    cost          :  float
    gradPred      :  (D, )
    grad          :  (V, D)
    '''
    predicted = predicted.reshape(-1, 1)                          # (D ,1)
    
    sampledIndices = [dataset.sampleTokenIdx() for i in xrange(K)]
    sampledVectors = outputVectors[sampledIndices, :]             # (K, D)

    outputVec = outputVectors[target, :]                          # (D, )

    prob_out = sigmoid(outputVec.dot(predicted))                  # float
    probs_negative = sigmoid(-sampledVectors.dot(predicted))      # (K, 1)
    
    cost = - np.log(prob_out) - np.sum(np.log(probs_negative))
    
    gradPred = (prob_out - 1) * outputVec - np.sum((probs_negative-1) * sampledVectors, axis=0)    # (D, )

    grad = np.zeros_like(outputVectors)                           # (V, D)
    grad[target, :] = (prob_out - 1) * predicted.reshape(-1)

    # Note that sampledIndices may have repeated indices, we may loop over all K samples.
    # And target should not be appeared in sampledIndices, but it's ok in gradient check,
    # because we use += to update each output vector, and no grads will be missed.
    for i in xrange(K):
        grad[sampledIndices[i], :] += (1 - probs_negative[i]) * predicted.reshape(-1)

    ### END YOUR CODE
    
    return cost, gradPred, grad

def test_sigmoid_gradient(dim_1, dim_2):
    a1    = np.random.normal(loc=0., scale=20., size=(dim_1,dim_2))
    shift = np.random.uniform(low=1e-9, high=1e-5, size=(dim_1,dim_2))
    ap = a1 + shift
    am = a1 - shift

    dsigmoid = (sigmoid(ap) - sigmoid(am)) / (2*shift)
    assert np.abs(np.max(dsigmoid - sigmoid_grad(sigmoid(a1)))) <= 1e-7
    assert np.abs(np.min(dsigmoid - sigmoid_grad(sigmoid(a1)))) <= 1e-7

def test_sigmoid_permutation_axis1(dim_1):
    a1          = np.random.normal(size=(1,dim_1))
    s1          = sigmoid(a1)

    permutation = np.random.permutation(dim_1)
    inverse_permutation = np.argsort(permutation)

    s1_perm     = sigmoid(a1.ravel()[permutation])
    assert rel_error(s1_perm.ravel()[inverse_permutation], s1) <= 1e-8

def test_sigmoid_permutation_axis0(dim_1, execution_number):
    """ sigmoid needs to be applied element-wise;"""
    a1          = np.random.normal(size=(dim_1,1))
    s1          = sigmoid(a1)

    permutation = np.random.permutation(dim_1)
    inverse_permutation = np.argsort(permutation)

    s1_perm     = sigmoid(a1[permutation])
    assert rel_error(s1_perm[inverse_permutation], s1) <= 1e-8

def test_sigmoid_shape(dim):
    testing_shape = []
    for y in range(0,dim):
        testing_shape.append(np.random.randint(3,8))
    shape = tuple(testing_shape)
    #z = np.random.randn(*testing_shape)
    x = np.random.standard_normal(shape)
    y = np.copy(x)
    assert x.shape == sigmoid(y).shape
    assert x.shape == sigmoid_grad(sigmoid(y)).shape

def sigmoid_forward(x):
    """
    Computes the forward pass for a sigmoid activation.

    Inputs:
    - x: Input data, numpy array of arbitary shape;

    Returns a tuple (out, cache)
    - out: output of the same shape as x
    - cache: identical to out; required for backpropagation
    """
    return sigmoid(x), sigmoid(x)

def negSamplingCostAndGradient(predicted, target, outputVectors, dataset, 
    K=10):
    """ Negative sampling cost function for word2vec models """

    # Implement the cost and gradients for one predicted word vector  
    # and one target word vector as a building block for word2vec     
    # models, using the negative sampling technique. K is the sample  
    # size. You might want to use dataset.sampleTokenIdx() to sample  
    # a random word index. 
    # 
    # Note: See test_word2vec below for dataset's initialization.
    #                                       
    # Input/Output Specifications: same as softmaxCostAndGradient     
    # We will not provide starter code for this function, but feel    
    # free to reference the code you previously wrote for this        
    # assignment!
    
    ### YOUR CODE HERE
    v_c = predicted
    u_o = outputVectors[target]

    suv = sigmoid(u_o.dot(v_c))
    pos = np.log(suv) # positive sample

    # sample w/ dataset.sampleTokenIdx() method iteratively
    n = []  # indexes of negative samples
    while len(n) < K:
        x = dataset.sampleTokenIdx()
        if x != target:
            n.append(x)
    neg_samples = outputVectors[n]


    skv = sigmoid((neg_samples.dot(v_c)))
    neg = np.sum(np.log(1-skv))

    cost = -pos - neg

    # neg_samples: K x d, skv: 1 x K
    gradPred = -(1 - suv) * u_o + (neg_samples.T * (skv)).sum(axis=1)
    
    grad = np.zeros(outputVectors.shape)
    grad[target] += -(1-suv) * v_c
    negGrad = np.outer(skv, v_c)
    # sum grads together when they have been sampled with replacement
    for i,x in enumerate(n):
        grad[x] += negGrad[i]
    ### END YOUR CODE

    return cost, gradPred, grad

def test_sigmoidgrad():
    """ Original sigmoid gradient test defined in q2_sigmoid.py; """
    x = np.array([[1, 2], [-1, -2]])
    f = sigmoid(x)
    g = sigmoid_grad(f)
    assert rel_error(g, np.array([[0.19661193, 0.10499359],
        [0.19661193, 0.10499359]])) <= 1e-7

def forward_backward_prop(data, labels, params, dimensions):
    """ 
    Forward and backward propagation for a two-layer sigmoidal network 
    
    Compute the forward propagation and for the cross entropy cost,
    and backward propagation for the gradients for all parameters.
    """

    ### Unpack network parameters (do not modify)
    ofs = 0
    Dx, H, Dy = (dimensions[0], dimensions[1], dimensions[2])

    W1 = np.reshape(params[ofs:ofs+ Dx * H], (Dx, H))
    ofs += Dx * H
    b1 = np.reshape(params[ofs:ofs + H], (1, H))
    ofs += H
    W2 = np.reshape(params[ofs:ofs + H * Dy], (H, Dy))
    ofs += H * Dy
    b2 = np.reshape(params[ofs:ofs + Dy], (1, Dy))

    ### YOUR CODE HERE: forward propagation
    # data: N x Dx, W1: Dx x H, b: 1 x H 
    a = data.dot(W1) + b1
    h = sigmoid(a)
    # h: N x H, W2: H x Dy, b2: 1 x Dy
    t = h.dot(W2) + b2
    y_hat = softmax(t)
    # y_hat: N x Dy, labels: N x Dy (as int)
    probs = labels * y_hat
    cost = np.sum(-np.log(probs.sum(axis=1)))
    ### END YOUR CODE
    
    ### YOUR CODE HERE: backward propagation
    # obtain the softmax gradient
    dJdt = (y_hat - labels) # N x Dy

    # b2 grad is sum along each index of the Dy vectors
    gradb2 = np.sum(dJdt, 0) 

    # h: N x H, dJdt: N x Dy
    gradW2 = h.T.dot(dJdt) # H x Dy

    # dJdt: N x Dy, W2: H x Dy
    dJdh = dJdt.dot(W2.T)
    # h: N x H
    dhda = sigmoid_grad(h)

    # data: N x Dx, dhda: N x H, DJdh: N x H
    gradW1 = data.T.dot(dhda * dJdh)
    
    # dhda: N x H, DJdh: N x H
    gradb1 = np.sum(dhda * dJdh, 0)
    ### END YOUR CODE
    
    ### Stack gradients (do not modify)
    grad = np.concatenate((gradW1.flatten(), gradb1.flatten(), 
        gradW2.flatten(), gradb2.flatten()))

    return cost, grad

def forward_backward_prop(data, labels, params, dimensions):
    """ 
    Forward and backward propagation for a two-layer sigmoidal network 
    
    Compute the forward propagation and for the cross entropy cost,
    and backward propagation for the gradients for all parameters.
    """

    ### Unpack network parameters (do not modify)
    ofs = 0
    Dx, H, Dy = (dimensions[0], dimensions[1], dimensions[2])

    W1 = np.reshape(params[ofs:ofs+ Dx * H], (Dx, H))
    ofs += Dx * H
    b1 = np.reshape(params[ofs:ofs + H], (1, H))
    ofs += H
    W2 = np.reshape(params[ofs:ofs + H * Dy], (H, Dy))
    ofs += H * Dy
    b2 = np.reshape(params[ofs:ofs + Dy], (1, Dy))

    ### YOUR CODE HERE: forward propagation
    #z1 = data.dot(W1) + b1
    #hidden = sigmoid(z1)
    #z2 = hidden.dot(W2) + b2
    #print 'z2.shape: ', z2.shape
    #prediction = softmax(z2)
    ### END YOUR CODE
    
    hidden = sigmoid(data.dot(W1) + b1)
    prediction = softmax(hidden.dot(W2) + b2)
    cost = -np.sum(np.log(prediction) * labels)

    
    ### YOUR CODE HERE: backward propagation
    #print 'NN: ', Dx, H, Dy
    #print 'b1.shape: ', b1.shape
    #print 'prediction.shape: ', prediction.shape
    #print 'labels.shape : ', labels.shape
    #print 'W2.shape: ', W2.shape
    #print 'hidden.shape: ', hidden.shape
    #print 'hidden.T.shape: ', hidden.T.shape
    #print 'delta.shape: ', delta.shape
    #print 'W1.shape: ', W1.shape
    #print 'data.shape: ', data.shape
    #gradW2 = delta * hidden
    #print 'sigmoid_grad(hidden).shape: ', sigmoid_grad(hidden).shape
    delta = prediction - labels
    gradW2 = hidden.T.dot(delta)
    gradb2 = np.sum(delta, axis = 0)
    hidden_delta = delta.dot(W2.T) * sigmoid_grad(hidden)
    gradW1 = data.T.dot(hidden_delta)
    gradb1 = np.sum(hidden_delta, axis = 0)
    ### END YOUR CODE
    
    ### Stack gradients (do not modify)
    grad = np.concatenate((gradW1.flatten(), gradb1.flatten(), 
        gradW2.flatten(), gradb2.flatten()))
    
    return cost, grad

def negSamplingCostAndGradient(predicted, target, outputVectors, dataset, 
    K=10):
    """ Negative sampling cost function for word2vec models """

    # Implement the cost and gradients for one predicted word vector  
    # and one target word vector as a building block for word2vec     
    # models, using the negative sampling technique. K is the sample  
    # size. You might want to use dataset.sampleTokenIdx() to sample  
    # a random word index. 
    # 
    # Note: See test_word2vec below for dataset's initialization.
    #                                       
    # Input/Output Specifications: same as softmaxCostAndGradient     
    # We will not provide starter code for this function, but feel    
    # free to reference the code you previously wrote for this        
    # assignment!
    
    ### YOUR CODE HERE
    K_set = []
    while (len(K_set) < K):
        candidateIndex = dataset.sampleTokenIdx()
        if (candidateIndex != target):
            K_set += [candidateIndex]

    cost = 0
    grad = np.zeros(outputVectors.shape)
    gradPred = np.zeros(predicted.shape)
    for k in K_set:
        score = sigmoid(- outputVectors[k].dot(predicted))
        cost += - np.log(score)
        grad[k,:] += - (score - 1) * predicted 
        gradPred += - (score - 1) * outputVectors[k]

    score = sigmoid(outputVectors[target].dot(predicted))
    cost += - np.log(score)
    grad[target,:] = (score - 1) * predicted
    gradPred += (score - 1) * outputVectors[target]

    # ugly fix ...
    gradPred = gradPred[np.newaxis, :]
        

    #raise NotImplementedError
    ### END YOUR CODE
    
    return cost, gradPred, grad

def test_sigmoid(dim_1, dim_2):
    a1       = np.random.normal(loc=0., scale=20., size=(dim_1,dim_2))
    a1_copy  = a1.copy()

    s_a1     = sigmoid(a1)
    s_sol_a1 = sigmoid_sol(a1_copy)

    assert rel_error(sigmoid_grad(s_a1), sigmoid_grad_sol(s_sol_a1)) <= 1e-10

def forward_backward_prop(data, labels, params, dimensions):
    """
    Forward and backward propagation for a two-layer sigmoidal network

    Compute the forward propagation and for the cross entropy cost,
    and backward propagation for the gradients for all parameters.
    """

    ### Unpack network parameters (do not modify)
    ofs = 0
    Dx, H, Dy = (dimensions[0], dimensions[1], dimensions[2])

    W1 = np.reshape(params[ofs:ofs + Dx * H], (Dx, H))
    ofs += Dx * H
    b1 = np.reshape(params[ofs:ofs + H], (1, H))
    ofs += H
    W2 = np.reshape(params[ofs:ofs + H * Dy], (H, Dy))
    ofs += H * Dy
    b2 = np.reshape(params[ofs:ofs + Dy], (1, Dy))

    ### YOUR CODE HERE: forward propagation
    # data : N * Dx
    # W1   : Dx * H
    # b1   : 1 * H
    # W2   : H * Dy
    # b2   : 1 * Dy
    N = data.shape[0]

    z1 = data.dot(W1) + b1
    a1 = sigmoid(z1)  # N * H
    z2 = a1.dot(W2) + b2
    a2 = softmax(z2)  # N * Dy

    cost = np.sum(-np.log(a2[labels == 1])) / N

    ### END YOUR CODE

    ### YOUR CODE HERE: backward propagation
    delta_score = a2 - labels  # 1 * Dy
    delta_score /= N

    gradW2 = np.dot(a1.T, delta_score)  # H * 1 * 1 * Dy = H * Dy
    gradb2 = np.sum(delta_score, axis=0)

    grad_h = np.dot(delta_score, W2.T)  # 1 * Dy * Dy * H = 1 * H
    grad_h = sigmoid_grad(a1) * grad_h

    gradW1 = np.dot(data.T, grad_h)
    gradb1 = np.sum(grad_h, axis=0)

    ### END YOUR CODE

    ### Stack gradients (do not modify)
    grad = np.concatenate((gradW1.flatten(), gradb1.flatten(),
                           gradW2.flatten(), gradb2.flatten()))

    return cost, grad

def negSamplingCostAndGradient(predicted, target, outputVectors, dataset,
                               K=10):
    """ Negative sampling cost function for word2vec models

    Implement the cost and gradients for one predicted word vector
    and one target word vector as a building block for word2vec
    models, using the negative sampling technique. K is the sample
    size.

    Note: See test_word2vec below for dataset's initialization.

    Arguments/Return Specifications: same as softmaxCostAndGradient
    """

    # Sampling of indices is done for you. Do not modify this if you
    # wish to match the autograder and receive points!
    indices = [target]
    indices.extend(getNegativeSamples(target, dataset, K))

    ### YOUR CODE HERE
    #raise NotImplementedError
    U = outputVectors
    u_o = U[target]
    v_c = predicted 
    N = U.shape[0]
    
    sig_o_c = sigmoid(np.dot(u_o, v_c))
    cost = -np.log(sig_o_c)
    gradPred = (sig_o_c - 1)*u_o
    grad = np.zeros((N, len(v_c)))
    grad[target] += (sig_o_c - 1)*v_c

    for k in indices:
        u_k = U[k]
        sig_k_c = sigmoid(-np.dot(u_k, v_c))
        cost += -np.log(sig_k_c)
        gradPred += -(sig_k_c - 1)*u_k
        grad[k] += -(sig_k_c - 1)*v_c

    ### END YOUR CODE

    return cost, gradPred, grad

def negSamplingCostAndGradient(predicted, target, outputVectors, dataset,
    K=10):
    """ Negative sampling cost function for word2vec models """

    # Implement the cost and gradients for one predicted word vector
    # and one target word vector as a building block for word2vec
    # models, using the negative sampling technique. K is the sample
    # size. You might want to use dataset.sampleTokenIdx() to sample
    # a random word index.
    #
    # Note: See test_word2vec below for dataset's initialization.
    #
    # Input/Output Specifications: same as softmaxCostAndGradient
    # We will not provide starter code for this function, but feel
    # free to reference the code you previously wrote for this
    # assignment!

    ### YOUR CODE HERE
    W,D = outputVectors.shape

    UK = np.zeros((K+1, D))
    indices = [target]
    for i in xrange(K):
        k = dataset.sampleTokenIdx()
        while k == target:
            k = dataset.sampleTokenIdx()
        indices.append(k)
    for i,ix in enumerate(indices):
        UK[i] = outputVectors[ix]

    u_o = outputVectors[target] # (D,)
    cost = - np.log(sigmoid(np.dot(u_o, predicted))) - np.sum(np.log(sigmoid(-np.dot(UK[1:], predicted))))
    gradPred = (sigmoid(np.dot(u_o,predicted))-1) * u_o + np.dot(UK[1:].T,sigmoid(np.dot(UK[1:], predicted))) # dJ/dV_c, (D,)

    y = np.zeros(K+1); y[0] = 1.0 #
    grad = np.zeros(outputVectors.shape)
    gradK = np.outer(sigmoid(np.dot(UK, predicted)) - y, predicted)
    for i,ix in enumerate(indices):
        grad[ix] += gradK[i]
    ### END YOUR CODE

    return cost, gradPred, grad