Examples¶

Guess initial responsibility¶

The first step is to make an initial guess of responsibility. sati.GuessInitRsp helps this process.

# importing packages
import numpy as np
import sati

# Image to be analyzed needs to be given as numpy.ndarray.
image = np.loadtxt('yourdata.csv')

# Suppose the image has 5 terraces of different heights.
# When difference of values between adjacent pixels is below
# a threshold, the pixels are regared in the same terrace.
# The unit of threshold is that of the height of the image (e.g., nm).
initrsp = sati.GuessInitRsp(image, n=5, threshold=5e-3)

# Show a guess of initial responsibility.
initrsp.show()

If it works, each terrace is assigned to a class. Since the seeds (starting points) of terrace search are randomly chosen, results are different every time. Once you get a nice initial responsibility, you can reproduce it from the seeds used for the guess that is saved in the attribute seeds.

Note

If there are two (or more) terraces that are supposed to be the same height in an image, all of them should be assigned to the same class, or only one of them should be assigned to a class and the other(s) should be left unassigned.

Subtract a linear plane¶

Using the initial responsibility above, analysis with a linear plane and Cauchy distribution can be done as follows.

# Polynomial plane of degree 1
poly = sati.planes.Poly()

# Cauchy distribution
dist = sati.distributions.Cauchy()

# Create an model and optimize it.
# initrsp.guess is the initial responsibility above.
m = sati.Model(image, rsp=initrsp.guess, dist=dist, poly=poly)
m.optimize()

Optimized parameters are given in the attributes. In the above example,

m.subtracted a plane-subtracted image
m.rsp estimated responsibility, i.e., soft clustering of terraces
m.dist.loc estimated location parameters of the distribution, i.e., heights of terraces
m.dist.scale estimated scale parameters of the distribution.

See API documentation for more details.

Subtract decays¶

Logarithmic / exponential decays are included by setting the decay.

m = sati.Model(image, rsp=initrsp.guess, poly=sati.planes.Poly(),
               decay=sati.planes.Decay(tau=500, coef=-2e-11),
               dist=sati.distributions.Cauchy())
m.optimize()

The unit of tau and coef are pixel and the unit of the height (e.g., nm), respectively. Depending on scan direction of an image, you need to set the orgdrct arguments to specify chronological order of pixels. See API documentation for more details.

Note

When kind=='exp', the sign of coef is opposite to the z-direction to which the decay goes. For example, when the z-direction is positive (the tip height gets larger as time goes on), coef is negative.

Use a previous result as an initial value¶

For example, when you subtract a polynomial surface of degree > 1, you may want to use coefficients obtained for a model with a lower degree as initial values to prevent the optimization from being trapped in a small local maximum. Here is such an example.

# Linear plane
m = sati.Model(image, rsp=initrsp.guess, poly=sati.planes.Poly(),
               dist=sati.distributions.Cauchy())
m.optimize()

# Use the result as an initial guess and optimize again.
m.poly = sati.planes.Poly(degree=2, coef=m.poly.coef)
m.optimize()

Estimate a unit height of steps¶

Estimating a unit height of steps is achieved by adopting von Mises-Fisher distribution as a prior distribution for the location parameter of a model distribution.

# Optimize a model to get a set of initial values for the next step
m = sati.Model(image, rsp=initrsp.guess, poly=sati.planes.Poly(),
               dist=sati.distributions.Cauchy())
m.optimize()

# Set von Mises-Fisher distribution as a prior probability.
m.prior = sati.distributions.VonMises(scale=2., kappa=[.1]*rsp.shape[0]))
m.optimize()

Note

len(kappa) must be the same as the number of terraces.

The scale parameter takes an initial guess of the lattice constant. The estimated lattice constant is given in the scale attribute (m.prior.scale in the above example).

Save interim result¶

You can save a result so that you can reuse it as a set of initial parameters for successive calculations. This is useful, for example, when you change parameters at the end of a series of long calculations.

For example, when you estimate an unit height of steps, you would repeat calculations with different kappa. If calculations before applying a prior distribution take long time, you may want to save the result to reuse it.

m = sati.Model(image, rsp=initrsp.guess, poly=sati.planes.Poly(),
               dist=sati.distributions.Cauchy())
m.optimize()

# Save the result for subsequent calculations
m.pickle('tmp.pickle')

# kappa = 0.01
m.prior = sati.distributions.VonMises(scale=2., kappa=[.01]*rsp.shape[0]))
m.optimize()

Suppose you decide to calculate with larger kappa after seeing the result of kappa=0.01. You can directly calculate the final step as follows.

m = sati.Model.unpickle('tmp.pickle')
m.prior = sati.distributions.VonMises(scale=2., kappa=[.1]*rsp.shape[0]))
m.optimize()

m = sati.Model.unpickle('tmp.pickle')
m.prior = sati.distributions.VonMises(scale=2., kappa=[1.]*rsp.shape[0]))
m.optimize()

If you reuse a result without saving it, you can also use copy.deepcopy of the copy module.