2

I have the following code. The beggining is quite long, but only serves to generate data. The problem happens with a few lines at the end.

##### Import packages
import numpy as np
import scipy.linalg as la

##### Initial conditions
N = 5
lamda = 7
mu = 2
a = 0.5
r = - np.log(a).copy()
St_Sp = np.arange(- N, N + 1)
Card = St_Sp.shape[0]

##### Define infintesimal generator
def LL(x, y):
    if x == N or x == - N: re = 0
    elif x - y == - 1: re = lamda
    elif x - y == 1: re = mu
    elif x - y == 0: re = - (mu + lamda)
    else: re = 0
    return re

def L(x):
    return - LL(x, x)

##### Define function Phi
def Phi(x): return max(x, 0)
Phi = np.vectorize(Phi)

##### Define vector b
b = Phi(St_Sp).copy()

##### Define function Psi
def Psi(x): return L(x) / (L(x) + r)
Psi = np.vectorize(Psi)

##### Generate a Boolean vector whose all elements are False
d = np.array([0] * Card).astype(bool)

##### Define matrix A
A = np.zeros((Card, Card))
for i in range(Card):
    for j in range(Card):
        if (i != j) & (L(St_Sp[i]) != 0):
            A[i, j] = LL(St_Sp[i], St_Sp[j]) / L(St_Sp[i])
        elif (i != j) & (L(St_Sp[i]) == 0):
            A[i, j] = 0
        elif (i == j) & (Psi(St_Sp[i]) != 0):
            A[i, j] = - 1 / Psi(St_Sp[i])
        else: A[i, j] = 1

##### Row names of A
rows = np.arange(0, Card)

##### Define matrix B
B = np.zeros((Card, Card))
for i in range(Card):
    for j in range(Card):
        if i != j:
            B[i, j] = LL(St_Sp[i], St_Sp[j])
        else: B[i, j] = LL(St_Sp[i], St_Sp[j]) - r

##### Generate I_0
I = [np.array([1] * Card).astype(bool), d.copy()]
Z = b.copy()       

index0 = np.matmul(B, Z) <= 0
index1 = ~ index0

##### Generate I_1
I = [index0, index1]
Z = b.copy()

if np.sum(I[1]) > 0:
    order = np.concatenate((rows[I[1]], rows[~ I[1]]))
    A1 = A[np.ix_(rows[I[1]], order)]    
    A2 = la.lu(A1)[2]
    p = np.atleast_2d(A1).shape[0]    
    B1 = A2[:, range(p)]          
    B2 = - np.matmul(A2[:, p:], Z[I[0]])

    print('Before being assigned new values, Z is \n', Z)

    print('\n The index I[1] of elements of Z to be change \n', I[1])

    M = la.solve_triangular(B1, B2,  lower = False)

    print('\n The values to be assigned to Z[I[1]] is \n', M)

    Z[I[1]] = M

    print('\n After being assigned new values, Z is \n', Z)

with result

Before being assigned new values, Z is 
 [0 0 0 0 0 0 1 2 3 4 5]

 The index I[1] of elements of Z to be change 
 [False False False False False  True  True  True  True  True False]

 The values to be assigned to Z[I[1]] is 
 [2.08686055 2.88974949 3.40529229 3.88978577 4.41338306]

 After being assigned new values, Z is 
 [0 0 0 0 0 2 2 3 3 4 5]

It's very weird to me that the command Z[I[1]] = M does not assign new values from M to the postion of Z indexed by I[1]. Could you please elaborate on why this problem arises and how to resolve it?

2
  • I posted an answer, but for the next time, please write a minimal reproducible example and look up on how to ask a good question. Thanks :) Commented Jun 17, 2020 at 17:14
  • @Dorian I very would like to do so, but the data generated is pertaining to how is its generated. Honestly, It's extremely hard to reproduce the problem without quoting the entire code. Commented Jun 17, 2020 at 17:19

1 Answer 1

5

The datatype of your array Z is int, to the values are typecasted by python automatically, resulting in the interger values of int([2.08686055 2.88974949 3.40529229 3.88978577 4.41338306]) = [2 2 3 3 4 5].

If you want to change that behavour, you just need to add a line to change the type of your original array:

Z = Z.astype(float)
Sign up to request clarification or add additional context in comments.

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.