Misunderstanding Multiscattering

Today we will build a state-of-the-art multiscattering GGX BRDF for physically-based rendering. 

Already done you say? Yes, you're right, by people who understand maths and physics, that's cheating. We instead will try to do the same trying to learn as little as possible...

If you want to do this the right way instead, these are some excellent articles:

• Stephen Hill - "A Multifaceted Exploration"
• Danny Chan "Material advances in Call of Duty: WWII"
• McAuley "A journey through implementing Multiscattering BRDFs and Area Lights"

Ok, ready? 

Let's start from your vanilla contemporary GGX with all the correlated-Smith fixings. The following sequence of images is generated with Disney's BRDF explorer environment lighting mode. They show a metallic GGX withy varying roughness, from 0.3 to 0.9, using the common parametrization alpha = roughness squared.

Note: please always indicate the parametrization you are using in your articles and presentations. There are so many variants out there. Here I use one of the most common ones, even if, for what is worth probably the best idea is alpha = roughess^4, which is more perceptually linear and also manages to be a great approximation for the one-over-square-root Blinn-Phong Gloss to GGX mapping many engines prefer.

We can notice how at high roughness the material looks noticeably darker. This is annoying, and smart people have worked hard to find and fix this issue. 

Let's compare the image above with Kulla/Conty's Multiscattering GGX BRDF.

Fixed! 

Now, we could try to understand how this solution work, the math is already there. And you probably should, it's both interesting and useful to understand other related problems. But, let's try not to...

First of all, we should ask why is this a problem? How do we know it is a problem? And also, if it is a problem, how big of a problem it is?

Let's start at the top. Why do we like "PBR"? It's not because we like physics, at least, I don't... And if we did like physics, this would be a small error to focus on, considering that we are surely committing worse sins in our end-to-end image pipeline...

Computer graphics is not predictive rendering. We don't try to simulate physics to make accurate simulations, that's not a goal (and the people who care about physics are an order of magnitude ahead of us). We make pretty pictures.

At a given point, we noticed that using some physics to make pretty pictures allowed for easier workflows, decoupling materials from lighting, reducing the number of hacks and parameters, allowing to use libraries of real materials and so on. 
We thought artists would be happier and be able to work more efficiently, while still being able to create a lot of different scenes. It took a while to move over all our workflows, but it was the right call.

So, when we look for right and wrong, we have to start with art and artists. If there is a need we can identify there, then we can look at physics to see if the answers are there. In this case, we can say that yes, indeed we have a (small) problem. 
It would be nice if our material parameters were orthogonal, and having things go darker as we change roughness is not ideal.

Let's move to step two then. Can physics help? Is this darkening physically correct, or not? How do we test? Enter the furnace! Let's put our "GGX coated" metallic object in a uniformly lit environment and see what happens:

Light hits our surface, hits our microfacets. The microfacets reflect back some light, and refract some other, according to Fresnel. If we're assuming a metal people told us that the refracted light, that goes "inside" the surface, is quickly absorbed and never comes out (becomes heat). 

It's reasonable that even in a furnace our metallic object has a color, as some of the energy will be absorbed. But what if we take our microfacets and make them always reflect all the light, no absorption? 
This means we need to set our f0 to 1 (remember, Fresnel is what controls how much light the microfacets scatter), let's try that and see what happens:

The object is still not white. Something's wrong! Ok, the inquisitive reader might still say - how do you know it's wrong. Perhaps certain directions scatter more light and others less, so we still see some shading even if all the light eventually comes out... 

If we think of BRDF and lobes it's not easy to get correct intuition. Let's instead think of how a light path, starting from the camera, looks like. It might hit some of the microfacets, one or many times, bounce around then eventually escape and connect to the furnace environment which will always be a emitting a given contant energy.

How much of that energy reaches the camera? All of it! Because regardless of which and how many microfacets we hit, all the energy is reflected and none is absorbed.

Now we have an intuition of why the image above should be white, it isn't, so we have a problem in our math...

If you studied BRDFs you know that in microfacet model we have a masking-shadowing visibility function that models which microfacets are occluded by others. What we don't typically model though is the fact that these occlusions are themselves microfacets, so the light should bounce around and eventually get out, not just be discarded. 
This is what the multiscattering models model and fix, and indeed if we did put in a furnace Kulla and Conty's GGX that we have shown before, it would generate a boring and correct white image for a fully reflective material, regardless of its roughness.

Kulla's model is not trivial though, and in many cases is likely not worth using to fix such a relatively minor problem. So. Can we be more ignorant? What if we knew how much light our BRDF gives out in a furnace, for a given roughness and viewing angle (fixing then again f0 = 1)? Could we just take that value and normalize the BRDF with it? 

Spoiler alert, we can. And not only we can but it's also trivial to do, because we already have this "furnace" value in most modern engines, in the look-up tables used for the popular split-sum image based lighting approximation. 

The split-sum table boils down the "BRDF in a furnace" (also known as directional albedo or directional-hemispherical reflectance) to a scale and bias (add) factor to be applied to the Fresnel f0 value.
In our case, we want to normalize considering f0=1 so all we do is to scale our BRDF lobe by one over bias(roughness,ndotv) + scale(roughness,ndotv). This is the result:

Indeed, we're getting some energy back at high roughness, and if we tested this in a furnace it would come up white, correctly, at f0 = 1. But it's also different from Kulla's, in particular, color is not quite as saturated in the rough materials. How comes? Well again, if you studied this problem already (cheater!) you know the answer.

Proper multiscattering adds saturation because as light hits more and more microfacets before escaping the surface, we pick up more color (raising a color to a power results in a more saturated color). But so how can this be physically wrong, but still correct in our test? 

Well, it's simple, really. By not simulating this extra saturation we are still energy conserving, but we changed the meaning of our BRDF parameters. The f0 "meaning" in our "ignorant" multiscattering BRDF is not the same as Kulla's, it results in a different albedo, but the BRDF itself is still energy conserving, it's just a different parametrization. 
Most importantly, I'd argue it's a better parametrization! Then again, remember our objectives. We don't do physics for physics' sake, we do it to help our production.

We wanted to make our parameters more orthogonal so that artists don't need to artificially "brighten" our BRDF at high roughness. If we went with the "more correct" solution (about that, this recent post by Narkowicz is a great read) we would add a different dependency, now roughness instead of darkening our materials it makes them more saturated, which would kind of defeat the purpose. 
You can imagine some scenarios where this might be desirable, but I'd say it's almost always wrong for our use-cases.

If we wanted to simulate the added saturation, there are a few easy ways. Again, going for the most "ignorant" (simple) we can just scale the BRDF by 1+f0*(1/(bias(roughness,ndotv) + scale(roughness,ndotv)) - 1), resulting in the following:

Now we are close to Kulla's solution. This is still not entirely correct if you wonder why Kulla's approximation is better, as it might not show in images, it's because ours doesn't respect reciprocity.

This might be important for Sony Imageworks, as certain offline path-tracing light transport algorithms do require it, but it's fairly irrelevant for us.

Ok, so now that we have found an approximation we're happy with we can (and should) go a step further and see what we are really doing. Yes, we have a formula, but it depends on some look-up tables. 
This isn't a big issue as we need these tables around for image-based lighting anyways, but it would be an extra texture fetch and in general it's always important to double-check our math, so let's visualize the 1/(bias(roughness,ndotv) + scale(roughness,ndotv)) function we're using:

It looks remarkably simple! Indeed, it's so simple it has a trivial approximation, you don't even need any sophisticated tool to find it so I'll just show it: 1 + 2*alpha*alpha * ndotv. This fits very well:

Approximation compared to the correct normalization factor (gray surface)

You can see that there is a bit of error in the furnace test. We could improve by doing a proper polynomial fit - yes, it turns out "2" and "1" in the formula above are not exactly the best constants. 
Actually defining "best" would be a problem all on its own, because doing a simple mean-square minimization on the normalization function doesn't really make a lot of sense (we should care about end visuals, perceptive measures, which angles matter more and so on), and I think we already spent too much time for such a small fix.
Moreover, the actual rendered images are really hard to distinguish from the table-based solution.

What is interesting is to see what could we achieve if we wanted to be even simpler, dropping the dependency from ndotv. It turns out in this case 1 + alpha*alpha does the trick decently as well.
This means that we're just applying a multiplicative factor to our BRDF to brighten it at high roughness, which just makes a lot of sense.

This of course adds more error and it starts to show as the BRDF shape changes, we get more energy at grazing angles on rough surfaces, but it might just be good enough depending on your needs:

Under-exposed to highlight that with the simpler approximation sometimes we lose light, sometimes we add.

What makes "10x" engineers. A complete hypothesis.

This happened on Twitter recently:

It's a long thread and I censored the author's name because it doesn't matter and I don't want to add to the hate that already naturally happens on any social media these days. 

To be honest, I'm still not entirely sure it's not just a parody, but I think it isn't. Regardless, it seemed a good excuse for a blog post.

Let's for a moment forget about how good or bad the term is. I don't particularly love discussions around words. I've never used "10x" in my career and I don't find it particularly great (also because it's linked to a certain start-up culture that I don't particularly enjoy), but I don't really want to open that can of worms.

What's wrong with a description of a "10x" engineer like the one above? 

I hope it's obvious but I'll spell it out. It tries to describe an engineer not in terms of the work they do, the actual output and results, but in some kind of mystical terms.

Like you can infer people's skills by looking at how they dress, what they eat, what time they come to the office and so on.

It doesn't help that the attributes this person uses as signals for a "10x engineer" are all kinda on the introverted side of the spectrum, to be charitable, as introversion per-se is a test of some kind of smartness. It's close to the old Hollywood trope of the computer guy being a weird, overweight, sexually frustrated white guy.

Ok, now that we closed the chapter on the critique of random people's tweets, let's get to something more interesting... 

Do "10x" engineers even exist? What are they? And what makes them?

Let's start from that last one. 

What's a 10x engineer? Typically we think of "10x" engineers as people that are much more productive than the average when working with code. Code "wizards", "ninjas", "rockstars" or whatever else cringe-worthy moniker teenagers use.

In my experience, 10x engineers do exist. Even controlling for seniority, knowledge, and skills, productivity is not uniform among people. I hope this is not controversial, one can know something, be even experienced in doing a given thing with a good track record, and yet not be as effective at doing it as others.

Should you hire only "10x" people? Definitely! Sort-of. In a way... We all look for excellent people, of course, and being able to distinguish good from great among a given seniority level is certainly important.

That said, there are lots of ways to be excellent beyond coding. In a decently-sized team other aspects might even be more important, mentoring, coordination, project management and so on. As things grow code usually tends towards being an implementation detail, so to speak, secondary to product and people concerns.

Even if we just look at coding, there are lots of kinds of engineers, people who are great at handling huge, foreign code-bases, people who are great at fixing things, people who are great at creating new things, people who are great architects and so on...

Lots of things in which you can be "10x" - but still, the concept of productivity being separate from skill generally holds.

These multipliers are also of the hardest to assess in interviews because again we're saying it doesn't correlate with simply what a person has on the CV or their ability to answer technical questions. Correctly characterizing where this productivity comes from is thus of utmost importance.

Why are some people more effective than others, given the same skillset? Where does that effectiveness come from? Is it their choice of editors? Their typing speed? Some sort of flow-related supernatural focus? I think not.

First Hypothesis: Output = Skill * Effort * Allocation

Skill is knowledge and experience, what we usually mostly correlate with seniority levels. As an analogy, I would say this is the value of the chips you have, at a gambling table.

Effort, given the same workday, is mostly focus. It's a time management skill, the ability to execute your tasks in good-sized chunks. In the gambling analogy, this would be the number of chips you have.

Focus is partially environment-related, but what we don't say often enough is that a lot of it is a skill. A lot of it also relates to how much a person likes doing a given thing, how fit he is for the job at hand. Effective managers that try to hire and allocate people to do the things they are passionate about, can thus help to get the most from the effort multiplier.

The last aspect, what I called allocation. This is how you spend the chips you have. Second hypothesis: correct allocation is what "makes" a 10x engineer.

In other words, we all have a number of chips to spend each day on our tasks. And controlling for seniority levels more or less effectively capture this number, there isn't that much variance in that.

The part that has a lot of variance is the allocation. Not in what to do, as in prioritizing this bugfix to that feature (even if such skills are also fundamental and have very high variance we capture them well in job categories, think technical director versus principal engineer for example) but how we do things.

Do I use a scripting language, should I implement things in C, or maybe I should learn that fancy new language everyone's talking about? Do I rely on a library or write from scratch? Do I need to understand the overall architecture of this software? Do I need to understand the specifics of the functions I'm calling? When it's appropriate to be sloppy? Should I jump into prototyping or I need to learn about the state of the art first? Should I go deep or wide?

We always have a limited amount of resources, ability to keep things in our brain, of doing work. And software design has a lot of different dimensions, abstraction versus specificity, generalization versus integration, high-level versus low-level concerns and so on. 
The ability to navigate this design space and selecting the right tools for the job, both in terms of concrete artifacts (code, libraries, languages, IDEs) and of abstract methodologies, makes a huge difference. One thing is to know about things, the other is to be able to critically evaluate tradeoffs and allocate (your) resources well.

Third Hypothesis. The reason why "10x" skills look mystical is that we don't have a solid theory for allocation choices.

When we have a design space that lacks a solid theoretical framework to navigate it, all success looks random, unteachable, and mystical. This is why "10x" engineers are both rare and sometimes described in horrible ways like it was done in the twitter thread above.

Too much of programming is still an art, eventually some people "get it" after lots of exercise, but we don't really know how to replicate that success.

We accept that somehow of the huge talent pool of software engineers, some will somehow find a way to be productive at some things, and some will be less successful.

---

Update: I was made aware of this, which is a much more concise way of vehiculating the same message:

It's great that I'm not alone in this...

The Value of Pixels (presentation slides)

Presented at the bay area game tech meetup, hosted at Roblox offices.
If you want to be notified of future meetups, join here.

Download the slides here!

Seeing the whole Physically-Based picture.

Subtitle: Building our rendering on solidly shaky grounds.

Physically-Based Rendering has won. There is no question about it, after an initial period of reluctance, even artists have been converted and I don't think you can find many rendering systems nowadays, either offline or real-time, that hasn't embraced PBR. And PBR proved itself to even be able to adapt to multiple art styles, outside strict adherence to photorealism.

But, really, how much physics is there in our PBR renderers? Let's have a look.

- Optics and Photometry.

Starting from the top, we have to define our physical framework. Physics are models, made to "fit" reality in order to make predictions. Different models are appropriate for different contexts and problems. For rendering, we work with a framework called "geometrical optics".

In G.O. light is composed of multiple frequencies which are assumed to be independent. Light travels in straight lines (in homogeneous media). It changes direction at changes of media (changes of IOR), where it can be absorbed, reflected or refracted. It travels instantaneously and it follows the path of least time (Fermat's principle). 

Is this a good framework? It's already making a lot of assumptions, and we know it cannot model all light behavior even when it comes to things that are easily visible: diffraction, interference, fluorescence, phosphorescence. But we say that these phenomena are not that common in everyday materials, and we might be right.

That's not all though, even before we start rendering our first triangle, we make more assumptions. First, we define a color space, usually a trichromatic one, because of the visual system metamerism. Fine, but we know that's not correct for rendering. We know spectral rendering has in even sometimes dramatically different results, but we trust our artists can tune lighting, materials, and post-processing in the right way (even if the two things shouldn't be related) to generate nice images even if we restrict ourselves to RGB. Or at least, we hope.

- Scattering

Next, we have to define what happens when the light "hits" something (an IOR discontinuity). Well, who knows, light is really hard! Some electrons... resonate? Get polarized? Please let it not be something to do with quantum stuff... Anyhow, eventually they scatter some energy back... waves? particles? There is some interference at around the atomic level. Who knows, luckily, we have another framework that comes to rescue: microfacet theory.

Surfaces are made of microfacets, like a microscopic landscape, light rays hit, bounce around and eventually come out. If we integrate the behavior of said microfacets over a small area, we can compute a scattering probability (BRDF) from the distribution of the microfacets themselves and a lot of math and voila', rendering happens.

Over a small area? How small, by the way? Well, Naty Hoffman and Eric Heitz say around the order of magnitude of the projected area of a pixel. I say, around the order of magnitude of a light wavelength, and then the projected area thing is antialiasing applied "after". So probably it's the pixel thing that's right.

What are these microfacets made of? Ideal reflectors obeying only the Fresnel law for how much light is reflected and how much refracted. The refracted part gets into the material (for dielectrics, that somehow allow this behavior), scatters some more and eventually comes out. If it comes out still "near enough" we call that "diffuse" reflection.

Otherwise, we call that subsurface scattering. But how does the light scatter inside the material? It hits particles. Microflakes? But microfacet based diffuse models (e.g. Oren-Nayar) simply swap the facets from ideal reflectors to ideal diffusers (Lambert)...

Regardless. We know all these things! We have blog posts, Siggraph talks, and books. Physics... And this still is well in that "geometrical optics" framework. Rays of light hit things. So much so that we can create raytracers to brute-force these microscopic interactions and create our own BRDFs!

But, it is still reasonable to use geometrical optics for these interactions? They seem to be quite... small. Maybe diffraction now matters? It turns out, it does, and it's a well-known thing (if you read the papers from the sixties... Beckmann-Spizzichino), but we sweep it under the rug.

And well, we can't really derive the equations from the microfacets, that integral is itself hard, so the BRDFs that we use introduce, typically, a bunch of other assumptions. For example, they might disregard the fact that light can bounce around multiple times before "coming out".

But who cares, nice pictures can be generated with this theory, and that's what matters. Moreover, someone did try to fit the resulting equations to real-world materials, right? The MERL database? I wonder how much error there is in that. Or how much it samples "well" real-world materials. Or how perceptual is the error metric used in estimating the error... Better to not think too much.

- Fiat Lux!

Are we done now? Far from it! In practice, we cannot just use the BRDF and brute-force light rays, not for real-time rendering, we're not Arnold. We need to compute a few more integrals!

We need to integrate over the light source, and over the surface area that is "seen" by the pixel, we're considering (pixel footprint). And that is incredibly hard, so hard we don't even try before having introduced a bunch more assumptions and approximations.

First of all, when we talk about pixel footprint, we really mean that we consider some statistics of the surface normals. We don't consider the fact that, for example, the "view rays" themselves change (and the light ones too), or that the surface normals don't really exist as an entity separate from actual surface geometry (which would cause shadowing and all other fun things). We assume these effects to be small.

Then, when we talk about light, we mostly mean simple geometric shapes that emit light from their surface. For example, a sphere. At each point, the light is emitted equally in all directions, and most often, it's also emitted with the same intensity over the surface.

And even then it's not enough to compute everything in closed-form. In fact, the more complex the light is, typically, the more approximated the BRDF will become. And then we'll fit these approximated BRDFs to the "real" one, and sum everything up. And sprinkle some of that pixel footprint thing on top somehow as well, but really that's done once on the "real" BRDF, even if we never actually use that!

So we have an approximation for very small lights, and maybe a good one for spheres, one for lines and capsules with some more handwaving, even more for polygonal lights, especially if textured and lastly one for far, "environment" light... We have approximations for "diffuse" and for "specular", for each of these. And maybe for static versus dynamic objects? A lot of math and different approximations.

We compare them and make sure that more-or-less the material looks the same under different kinds of light, and call it a day... The most ambitious of us might even export a scene to a path-tracer and compute some sort of ground-truth, to visually make sure things are at least reasonable...

- We're done, right?

So... we get our final image. In all its Physically Based, 60fps, HDR glory! Spectacular.

Year after year people come up with better equations, tighter approximations, and we make shinier pixels as a result.

Is that all? Of course not! We are just getting started! 

In practice, materials are not just one surface... They can have layers! And they are never optically uniform! They sparkle! They are anisotropic, they have scratches. Really, look around, look at most things. Most things are sparkly and anisotropic, due to the way they are fabricated.

And nothing is a surface, really. It's mostly volume and particles. Even... air! So we need fog and volumetric models. But that's not just about the light that scatters in the air back to our virtual cameras, we should also consider how this scattering affects lighting at surfaces. Our rays of light are not that straight anymore! Participating volumes make our light sources more "diffuse". Bigger. All of them, also things like environment lighting! And... that should affect shadows too right?

And now that we think about shadows... all this complexity and unknowns are still only for what we call "direct" lighting! What about global illumination? What about the million other hacks and assumptions that we rely upon to render each or our frames?

- Conclusions

So. How much physics is there in a frame, really? And more importantly, what's the point of all this? Should we be ashamed of not knowing physics that well? Should we do physics more? Less?

I don't know. I personally do not know physics well and I'm not too ashamed. A lot of what we've been doing is reasonable, really. We went with GGX because its "tail" helps. All the lighting improvements served our products. All the assumptions, individually, looked reasonable.

But, there is a value I think in looking at our math and our approximations holistically, now that we are getting so good at photorealism.

Perhaps there is not too much value for example in going off the deep end of complexity when we think of BRDFs, if we can't then integrate them with complex lighting, or, in order to do so, we have to approximate them again.

Similarly, the features we focus on should be evaluated. Is it more important to have non-uniform emission in our light sources, or a different "tail" in GGX/T-R? Anisotropic surfaces or sparkles? Spectral sampling? Thin-film? Non-lambertian diffuse? Of which kind? Accurate energy conservation or multi-bounce in microfacets?

Is it better to use the best possible approximation for a given integral, even if we end up with many different ones, or should we just use a bunch of spherical gaussians, or LTCs and such, but keep the same representation everywhere? And in general, is most of our error still in the materials, or in the lights? This is very hard to tell from just looking at artist-made pictures, because artists will compensate any way they can!

But even more importantly - How much can we keep relying on simplifying assumptions in order to make our math work?

I suspect to answer these questions we'll need more data. Acquire from the real world. Brute force solutions. Then look at the data and understand what matters, what matters perceptually, what errors are we committing end-to-end, and what we should approximate better and how...

And we should not assume that because somewhere we have a bit of physics, we are doing things correctly. We are, after all, a field that forgot for decades basic things like color spaces and gamma.

NumPy by Example

This originally was in my Scientific Python 101 article, I've split it now as it was a long article and sometimes I need just to have a look at this code as a reminder of how things work.
If you're interested in a similar "by example" introduction for Mathematica, check this other one out.

Execute each section in an IPython cell (copy and paste, then shift-enter)

# Set-up the environment:
# First, we'll have to import NumPy and SciPy packages.

# IPython has "magics" (macros) to help common tasks like this:

%pylab

# This will avoid scientific notation when printing, see also %precision

np.set_printoptions(suppress = True) 

# Optional: change the way division works:
#from __future__ import division 
# 1/2 = 0.5 instead of integer division... problem is that then 2/2 would be float as well if we enable this, so, beware.

# In IPython notebook, cell execution is shift-enter

# In IPython-QT enter evaluates and control-enter is used for multiline input

# Tab is used for completion, shift-tab after a function name shows its help

# Note that IPython will display the output of the LAST -expression- in a cell

# but variable assignments in Python are NOT expressions, so they will

# suppress output, unlike Matlab or Mathematica. You will need to add print

# statements after the assignments in the following examples to see the output

data = [[1,2],[3,4],[5,6]]

# if you want to suppress output in the rare cases it's generated, use ;

data; # try without the ; to see the difference

# There is a magic for help as well

%quickref

# Also, you can use ? and ?? for details on a symbol

get_ipython().magic(u'pinfo %pylab')

# In addition to the Python built-in functions help() and dir(symbol)

# Other important magics are: %reset, %reset_selective, %whos

# %prun (profile), %pdb (debug), %run, %edit

1*3 # evaluate this in a separate cell...

_*3 # ...now _ refers to the output of last evaluated cell

# again, not that useful because you can't refer to the previous expression

# evaluated inside a cell, just to the previous evaluation of an entire cell

# A numpy array can be created from a homogeneous list-like object

arr = np.array(data)

# Or an uninitialized one can be created by specifying its shape

arr = np.ndarray((3,3))

# There are also many other "constructors"

arr = np.identity(5)

arr = np.zeros((4,5))

arr = np.ones((4,))

arr = np.linspace(2,3) # see also arange and logspace

arr = np.random.random((4,4))

# Arrays can also be created from functions

arr = np.fromfunction((lambda x: x*x), shape = (10,))

# ...or by parsing a string

arr = np.fromstring('1, 2', dtype=int, sep=',')

# Arrays are assigned by reference

arr2 = arr

# To create a copy, use copy

arr2 = arr.copy()

# An array shape is a descriptor, it can be changed with no copy

arr = np.zeros((4,4))

arr = arr.reshape((2,8)) # returns the same data in a new view

# numpy also supports matrices, which are arrays contrained to be 2d

mat = np.asmatrix(arr)

# Not all operations avoid copies, flatten creates a copy

arr.flatten()

# While ravel doesn't

arr = np.zeros((4,4))

arr.ravel()

# By default numpy arrays are created in C order (row-major), but 

# you can arrange them in fortran order as well on creation,

# or make a fortran-order copy of an existing array

arr = np.asfortranarray(arr)

# Data is always contiguously packed in memory

# Arrays can be indexed as with python lists/tuples

arr = np.zeros((4,4))

arr[0][0]

# Or with a multidimensional index

arr[0,0]

# Negative indices start from the end

arr[-1,-1] # same as arr[3,3] for this array

# Or, like matlab, via splicing of a range: start:end

arr = arange(10)

arr[1:3] # elements 1,2

arr[:3] # 0,1,2

arr[5:] # 5,6,7,8,9

arr[5:-1] # 5,6,7,8

arr[0:4:2] # step 2: 0,3

arr = arr.reshape(2,5)

arr[0,:] # first row

arr[:,0] # first column

# Also, splicing works on a list of indices (see also choose)

arr = arr.reshape(10)

arr[[1,3,5]]

# and with a numpy array of bools (see also where)

arr=np.array([1,2,3])

arr2=np.array([0,3,2])

arr[arr > arr2]

# flat returns an 1D-iterator

arr = arr.reshape(2,5)

arr.flat[3] # same as arr.reshape(arr.size)[3]

# Core operations on arrays are "ufunc"tions, element-wise

# vectorized operations

arr = arange(0,5)

arr2 = arange(5,10)

arr + arr2

# Operations on arrays of different shapes follow "broadcasting rules"

arr = np.array([[0,1],[2,3]]) # shape (2,2)

arr2 = np.array([1,1]) # shape (1,)

# If we do arr+arr2 arr2.ndim

# an input can be used if shape in all dimensions sizes, or if 

# the dimensions that don't match have size 1 in its shape

arr + arr2 # arr2 added to each row of arr!

# Broadcasting also works for assignment:

arr[...] = arr2 # [[1,1],[1,1]] note the [...] to access all contents

arr2 = arr # without [...] we just say arr2 refers to the same object as arr

arr[1,:] = 0 # This now is [[1,1],[0,0]]

# flat can be used with broadcasting too

arr.flat = 3 # [[3,3],[3,3]]

arr.flat[[1,3]] = 2 # [[3,2],[2,3]]

# broadcast "previews" broadcasting results

np.broadcast(np.array([1,2]), np.array([[1,2],[3,4]])).shape

# It's possible to manually add ones in the shape using newaxis

# See also: expand_dims

print(arr[np.newaxis(),:,:].shape) # (1,2,2)

print(arr[:,np.newaxis(),:].shape) # (2,1,2)

# There are many ways to generate list of indices as well

arr = arange(5) # 0,1,2,3,4

arr[np.nonzero(arr)] += 2 # 0,3,4,5,6

arr = np.identity(3)

arr[np.diag_indices(3)] = 0 # diagonal elements, same as np.diag(arg.shape[0])

arr[np.tril_indices(3)] = 1 # lower triangle elements

arr[np.unravel_index(5,(3,3))] = 2 # returns an index given the flattened index and a shape

# Iteration over arrays can be done with for loops and indices

# Oterating over single elements is of course slower than native

# numpy operators. Prefer vector operations with splicing and masking

# Cython, Numba, Weave or Numexpr can be used when performance matters.

arr = np.arange(10)

for idx in range(arr.size):

    print idx

# For multidimensiona arrays there are indices and iterators

arr = np.identity(3)

for idx in np.ndindex(arr.shape):

    print arr[idx]

for idx, val in np.ndenumerate(arr):

    print idx, val # assigning to val won't change arr

for val in arr.flat:

    print val # assigning to val won't change arr

for val in np.nditer(arr):

    print val # same as before

for val in np.nditer(arr, op_flags=['readwrite']):

    val[...] += 1 # this changes arr

# Vector and Matrix multiplication are done with dot

arr = np.array([1,2,3])

arr2 = np.identity(3)

np.dot(arr2, arr)

# Or by using the matrix object, note that 1d vectors

# are interpreted as rows when "seen" by a matrix object

np.asmatrix(arr2) * np.asmatrix(arr).transpose()

# Comparisons are also element-wise and generate masks

arr2 = np.array([2,0,0])

print (arr2 > arr)

# Branching can be done checking predicates for any true or all true

if (arr2 > arr).any():

    print "at least one greater"

# Mapping a function over an array should NOT be done w/comprehensions

arr = np.arange(5)

[x*2 for x in arr] # this will return a list, not an array

# Instead use apply_along_axis, axis 0 is rows

np.apply_along_axis((lambda x: x*2), 0, arr)
# apply_along_axis is equivalent to a python loop, for simple expressions like the above, it's much slower than broadcasting

# It's also possible to vectorize python functions, but they won't execute faster

def test(x):

    return x*2

testV = np.vectorize(test)

testV(arr)

# Scipy adds a wealth of numerical analysis functions, it's simple

# so I won't write about it.

# Matplotlib (replicates Matlab plotting) is worth having a quick look.

# IPython supports matplotlib integration and can display plots inline

%matplotlib inline

# Unfortunately if you chose the inline backend you will lose the ability

# of interacting with the plots (zooming, panning...)

# Use a separate backend like %matplotlib qt if interaction is needed

# Simple plots are simple

test = np.linspace(0, 2*pi)

plt.plot(test, np.sin(test)) # %pylab imports matplotib too

plt.show() # IPython will show automatically a plot in a cell, show() starts a new plot

plt.plot(test, np.cos(test)) # a separate plot

plt.plot(test, test) # this will be part of the cos plot

#plt.show()

# Multiple plots can also be done directly with a single plot statement

test = np.arange(0., 5., 0.1)

# Notes that ** is exponentiation. 

# Styles in strings: red squares and blue triangles

plt.plot(test, test**2, 'rs', test, test**3, 'b^')

plt.show()

# It's also possible to do multiple plots in a grid with subplot

plt.subplot(211)

plt.plot(test, np.sin(test))

plt.subplot(212)

plt.plot(test, np.cos(test), 'r--')

#plt.show()

# Matplotlib plots use a hierarchy of objects you can edit to

# craft the final image. There are also global objects that are

# used if you don't specify any. This is the same as before

# using separate objects instead of the global ones

fig = plt.figure()

axes1 = fig.add_subplot(211)

plt.plot(test, np.sin(test), axes = axes1)

axes2 = fig.add_subplot(212)

plt.plot(test, np.cos(test), 'r--', axes = axes2)

#fig.show()

# All components that do rendering are called artists

# Figure and Axes are container artists

# Line2D, Rectangle, Text, AxesImage and so on are primitive artists

# Top-level commands like plot generate primitives to create a graphic

# Matplotlib is extensible via toolkits. Toolkits are very important.

# For example mplot3d is a toolkit that enables 3d drawing:

from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')

points = np.random.random((100,3))

# note that scatter wants separate x,y,z arrays

ax.scatter(points[:,0], points[:,1], points[:,2])

#fig.show()

# Matplotlib is quite extensive (and ugly) and can do much more

# It can do animations, but it's very painful, by changing the

# data artists use after having generated them. Not advised.

# What I end up doing is to generate a sequence on PNG from

# matplotlib and play them in sequence. In general MPL is very slow.

# If you want to just plot functions, instead of discrete data,

# sympy plots are a usable alternative to manual sampling:

from sympy import symbols as sym

from sympy import plotting as splt

from sympy import sin

x = sym('x')

splt.plot(sin(x), (x, 0, 2*pi))