implicit_ODE_With_Jac.py 文件源码

def implicit_black_box(propensities,V,X,w,h,deter_vector,stoc_positions, positions, valid,jac):

    # Adjustment for systems reaching steady state
    """
    temp = derivative_G(propensities,V,X,w,deter_vector,stoc_positions,positions,valid,jac)
    #pdb.set_trace()
    valid_adjust_pos = np.where(np.sum(np.abs(temp),axis=0) < 1e-10,True,False)

    valid_adjust = valid[:,:]
    valid_adjust[valid_adjust_pos,:] = False

    print(" Reached Steady State %d"%(np.sum(valid_adjust_pos)))
    """
    from scipy.integrate import ode
    #pdb.set_trace()
    deter_ode = ode(f).set_integrator('vode',method='bdf', with_jacobian=False)
    deter_ode.set_initial_value(X[deter_vector,:].flatten(), 0)
    #deter_ode.set_f_params([propensities,V,X,deter_vector,stoc_positions, positions,valid_adjust,w,jac])
    deter_ode.set_f_params([propensities,V,X,deter_vector,stoc_positions, positions,valid,w,jac])

    #pdb.set_trace()
    while deter_ode.successful() and deter_ode.t < h:
        deter_ode.integrate(h)
    #print("Black Box: \n"+ str(deter_ode.y))

    #print("iterator : \n:"+str(next_X[deter_vector,:]))
    X[deter_vector,:] = deter_ode.y.reshape(( np.sum(deter_vector),X.shape[1]))


    # Another adjust to compensate for non negative
    X = np.where(X < 0.0,0.0,X)

    return X