We’re rewarding the question askers & reputations are being recalculated! Read more.

Hot answers tagged

48

Python (as of 2.6 and 3.0) now searches in the ~/.local directory for local installs, which do not require administrative privileges to install, so you just need to point your installer to that directory. If you have already downloaded the package foo and would like to install it manually, type: cd path/to/foo python setup.py install --user If you are ...


40

I'm going to break up my answer into three parts. Profiling, speeding up the python code via c, and speeding up python via python. It is my view that Python has some of the best tools for looking at what your code's performance is then drilling down to the actual bottle necks. Speeding up code without profiling is about like trying to kill a deer with an ...


37

I work in a lab that does global optimization of mixed-integer and non-convex problems. My experience with open source optimization solvers has been that the better ones are typically written in a compiled language, and they fare poorly compared to commercial optimization packages. If you can formulate your problem as an explicit system of equations and ...


32

fmincon(), as you mentioned, employs several strategies that are well-known in nonlinear optimization that attempt to find a local minimum without much regard for whether the global optimum has been found. If you're okay with this, then I think you have phrased the question correctly (nonlinear optimization). The best package I'm aware of for general ...


32

First, if your undergraduates are like ours and had no prior introduction to computers, expect to spend some time teaching them how to use basic stuff like using a proper editor (i.e., not MS Word), the command line, etc. I think the answer somewhat depends on where you set the focus of your course (or what you are required to teach). For example: How ...


31

This is indeed called catastrophic cancellation. In fact, this particular case is very easy: rewrite the function using the equivalent, numerically stable expression $$ \frac{t}{1+\sqrt{1-t^2}}. $$ Since you probably need a reference, this is discussed in most numerical methods textbooks in relation to the formula for solving quadratic equations (that ...


29

Ease of learning Python and Fortran are both relatively easy-to-learn languages. It's probably easier to find good Python learning materials than good Fortran learning materials because Python is used more widely, and Fortran is currently considered a "specialty" language for numerical computing. I believe the transition from Python to Fortran would be ...


23

You have to subclass the rv_continuous class in scipy.stats import scipy.stats as st class my_pdf(st.rv_continuous): def _pdf(self,x): return 3*x**2 # Normalized over its range, in this case [0,1] my_cv = my_pdf(a=0, b=1, name='my_pdf') now my_cv is a continuous random variable with the given PDF and range [0,1] Note that in this example ...


23

I am not super familiar with f2py internals, but I am very familiar with wrapping Fortran. F2py just automates some or all of the things below. You first need to export to C using the iso_c_binding module, as described for example here: http://fortran90.org/src/best-practices.html#interfacing-with-c Disclaimer: I am the main author of the fortran90.org ...


22

In 2014, I would've said Python. In 2017, I wholeheartedly believe that the language to teach undergraduates is Julia. Teaching is always about a tradeoff. On one hand, you want to choose something that is simple enough that it is easy to grasp. But secondly, you want to teach something that has staying power, i.e. something that can grow with you. The ...


19

You can use the Python builtin ctypes module as described on fortran90.org. It is pretty straight forward and doesn't require any external dependencies. Also, the ndpointer arg type helper is very handy.


19

I will try to answer your question considering that you are asking for Python specifically. I will describe my own method of tackling a simulation problem. Strategies for faster simulations are given in this description. First, I prototype new simulations in Python. Of course, I try to make use of NumPy and SciPy as much as I can. Whereas NumPy provides a ...


19

Let me try and break down your requirements: Maintainability Reading/writing text data Strong interfaces/capability for LU factorizations Sparse linear solvers Performance and scalability to large data From this list, I would consider the following languages: C, C++, Fortran, Python, MATLAB, Java Julia is a promising new language, but the community is ...


19

Joblib does what you want. The basic usage pattern is: from joblib import Parallel, delayed def myfun(arg): do_stuff return result results = Parallel(n_jobs=-1, verbose=verbosity_level, backend="threading")( map(delayed(myfun), arg_instances)) where arg_instances is list of values for which myfun is computed in parallel. The main ...


18

There are two issues that you are likely to be encountering. Ill-conditioning First, the problem is ill-conditioned, but if you only provide a residual, Newton-Krylov is throwing away half your significant digits by finite differencing the residual to get the action of the Jacobian: $$ J[x] y \approx \frac{F(x+\epsilon y) - F(x)}{\epsilon}$$ If you ...


17

Here is the Numba solution. On my machine the Numba version is >1000x faster than the python version without the decorator (for a 200x200 matrix, 'k' and 200-length vector 'a'). You can also use the @autojit decorator which adds about 10 microseconds per call so that the same code will work with multiple types. from numba import jit, autojit @jit('f8[:]...


15

You should check out sympy.stats. It provides an interface to deal with random variables. The following example provides a random variable X defined on the unit interval with density 2x In [1]: from sympy.stats import * In [2]: x = Symbol('x') In [3]: X = ContinuousRV(x, 2*x, Interval(0, 1)) In [4]: P(X>.5) Out[4]: 0.750000000000000 In [5]: Var(X) # ...


15

There are a number of issues in your question. Do not use Gaussian Elimination (LU factorization) to calculate the numerical rank of a matrix. LU factorization is unreliable for this purpose in floating-point arithmetic. Instead, use a rank-revealing QR decomposition (such as xGEQPX or xGEPQY in LAPACK, where x is C, D, S, or Z, though those routines are ...


15

A difficulty with any of these types of questions is that the answer is highly community-dependent. To answer some of your questions in haphazard order: MATLAB is used a lot both in academia and in industry. One of the reasons it's used quite a bit in industry is because it is taught in academia. I know for a fact that MATLAB is used at Lincoln Laboratory ...


15

The question has two very different subquestions. I will address the first one only. Matlab's version runs on average 24 times faster than my python equivalent! The second one is subjective. I would say that letting know the user that there is some problem with the integral is a good thing and this SciPy behavior outperforms the Matlab`s one to keep it ...


14

Update: see the new GEKKO package that we just released. APM Python is a free optimization toolbox that has interfaces to APOPT, BPOPT, IPOPT, and other solvers. It provides first (Jacobian) and second (Hessian) information to the solvers and provides an optional web-interface to view results. The APM Python client is installed with pip: pip install ...


14

For the first part of my question, I found this very useful comparison for performance of different linear interpolation methods using python libraries: http://nbviewer.ipython.org/github/pierre-haessig/stodynprog/blob/master/stodynprog/linear_interp_benchmark.ipynb Below is list of methods collected so far. Standart interpolation, structured grid: http:/...


14

Here is R1, as computed in MATLAB: 1.0e+07 * -7.382605957465515 -9.599867106092937 -2.830412177259742 -0.000000000002830 -0.000000000002830 -1.230434326244253 -1.599977851015490 -0.471735362876624 -0.000000000000472 -0.000000000000472 3.691302978732758 4.799933553046468 1.415206088629871 0.000000000001415 0.000000000001415 -5....


14

First, see Mark L. Stone's answers, which is completely correct. Second, realize that this is the reason why people told you to use relative errors in your numerical analysis class. :) Third, the real question here is why the results do not coincide exactly, since both languages call some BLAS library functions for their computations. There are several very ...


13

Algorithms for rank-1 updates of the SVD (also called incremental SVD) do exist, but I haven't been able to find a LAPACK-like implementation anywhere. The one I've seen mentioned repeatedly is that of Brand (2003). Judging from this website, it seems as though Brand's algorithm is relatively simple to implement using existing LAPACK and BLAS routines as ...


13

To the best of my knowledge, Numpy does not support independent streams. Indeed, getting independent streams from the Mersenne Twister (Pythons RNG) is notoriously difficult although it can be done. Consider using the RandomGen package. It is fully compatible with Numpy, and provides you with the PCG64 generator, supporting up to $2^{63}$ independent ...


12

For logging that allows full reproducibility, I highly recommend the Sumatra python package. It nicely links the version control commit number, machine state, and output files to each program run and has a django web interface to interact with the database of run info. The python API makes it very easy to include logging in my scripts.


12

I will address only the comparison of C to C++. While it is true that anything written in C can be ported to C++ with a few syntactic touch-ups, the communities have different values. The C library community, more than almost any other, values binary stability. Binary stability is critical for low-level libraries to avoid inflicting constant pain on those ...


11

If you need higher-order derivatives, you won't get good results using equidistant data points. If you can sample your function at arbitrary nodes, I would recommend using Chebyshev points, i.e. $$x_k = \cos\left( \pi \frac{k}{n} \right), \quad k=0\dots n$$ for a polynomial of degree $n$. You can evaluate the polynomial stably using Barycentric ...


11

This is a widely-held concern in the scientific programming community, and I would consider the performance uncertainty to be one of the major "myths" in computational science. As @fcruz discusses, petsc4py is a wrapper to the PETSc libraries, not a reimplementation of PETSc in Python. Therefore, you can expect any performance penalties to come from ...


Only top voted, non community-wiki answers of a minimum length are eligible