Model fitting, getting fit-models as array

[Bernhard Schaffer]

Hi,

I'm currenlty trying to get a hold on concepts, possibilities and limitations of using Python & Hyperspy, without a specific goal in mind.

I'm also a total newcomer to both HyperSpy and Python, so please keep this in mind when some questions seem super trivial, I'm having a steep learning curve ahead of me...

Full disclosure: I'm currenlty also working for Gatan :c)

Now to my questions,

I've worked my way through various of the documenation and examples, and I'm currenlty dealing with the fitting routines.

What I have not been able to find out yet is, how I would get to a components fitted representation, i.e. the data values this particular component has (after the fit) over the data range.

I know how I can plot this, but how do I get the actual values of that plot? (As numpy array, I guess.)

My script so far.

import numpy as np
import hyperspy.api as hs
filename = "C:\\test.dm3"  # some STEM EELS SI
dataSI = hs.load(filename)
print('Source data:', dataSI )
dataSI.add_elements(('Ti', 'N', 'O'))
fitModel = dataSI.create_model()    # no LowLoss
print('Fit model:')
print(fitModel.components)
fitModel.smart_fit()
fitModel.quantify()
print('Fit model parameters:')
fitModel.print_current_values()
fitModel.plot()
print('----------\nAll Done\n----------')

At the moment I'm only interested in getting the "plot" of the various component fits, but I'm actually also interested in getting to the "maps" of multifits.

I've managed to perform the multifit and use the export_results() command to save the maps as RPL files (TIF didn't work for me, not quite sure why), but

how would I get to these map results as arrays within the Python script?

As an asside: Is there also a "raw" saving format, i.e. just a dump of the actual binary values? (In whatever ordering spatial-before-spectral or the the other way around.)

Or is this just not needed, because one could do that with the numpy arrays (once one gets them) ? If so, could you please provide a example code snippet?

I'm also struggling a bit with the documentation itself.

The link "Full HyperSpy API Documentation" sounds very promising, but somehow I seem to get only a list of many modul and package names and very little information what they contain or what the acording methods are.

Am I just using the documenation wrongly, or is it just not that complete?

best regards,

  Bernhard

[Francisco]

Hi Bernhard,

Nice to read you here.

If you just want to plot the components you could use m.plot(plot_components=True). If you would like to extract the data to plot it manually, simply call the components e.g. gaussian_data = gaussian_component(). To get the map results you can use the ``as_signal`` method of the component parameters e.g. gaussian_component.as_signal().

Best regards,

Francisco

[Bernhard Schaffer]

Hi Francisco,

thans for the quick reply. I'm trying to get a grip on all those concepts at once ( - just now figured that I should use the TAB for auto-complete and getting an idea of what I can do with what object ),

so at any time some code-snippets are welcome ;c)

From your message, I gather that I can/should do

import numpy as np
import hyperspy.api as hs
filename = "C:\\test.dm3"  # some STEM EELS SI
dataSI = hs.load(filename)
print('Source data:', dataSI )
dataSI.add_elements(('Ti', 'N', 'O'))
fitModel = dataSI.create_model()    # no LowLoss
print('Fit model:')
print(fitModel.components)
fitModel.smart_fit()
fitModel.quantify()
print('Fit model parameters:')
fitModel.print_current_values()

fitModel["PowerLaw"].r.as_signal

print('----------\nAll Done\n----------')

[Bernhard Schaffer]

Argh.

I accidently hit "enter" in the form and Google posted right away. Sorry, the email was not finished.
I still don't quite get how I would plot the "powerlaw" funcition over the fitted data range (and get the result as a numpy array).

Could you provide a little code snippet, please?

Hi Francisco,

thanks for the quick reply. I'm trying to get a grip on all those concepts at once ( - just now figured that I should use the TAB for auto-complete and getting an idea of what I can do with what object ),

so at any time some code-snippets are welcome ;c)

From your message, I gather that I can/should do

import numpy as np
import hyperspy.api as hs
filename = "C:\\test.dm3"  # some STEM EELS SI
dataSI = hs.load(filename)
print('Source data:', dataSI )
dataSI.add_elements(('Ti', 'N', 'O'))
fitModel = dataSI.create_model()    # no LowLoss
print('Fit model:')
print(fitModel.components)
fitModel.smart_fit()
fitModel.quantify()
print('Fit model parameters:')
fitModel.print_current_values()

fitModel["PowerLaw"].r.as_signal # ???

[Thomas Aarholt]

Hi Bernhard,

Quick note: You can access the github "chat" functionality here for faster conversation, although using the mail system is perfectly fine too:

https://gitter.im/hyperspy/hyperspy

The as_signal() method that Francisco gave you an example of provides you a map of that parameter's value across the sample. Imagine if you were fitting a Gaussian to the white line on a Core Loss L3 edge. Then core_loss_component.centre.as_signal() would provide you with a new hyperspy signal that showed how the centre of the component varied across the map. - but you would have to run m.multifit() first, in order to have fitted the component in all pixels on the sample.

I don't think that is the functionality you are looking for.

If I read your email correctly, you want to extract an array describing the the shape of the PowerLaw after fitting.

The `component.function(x)` method will do this for the current pixel you are viewing with the plot. It takes the signal energy-axis as input, so you have to give it that:

data_array = fitModel['Double Power Law'].function(fitModel.axis.axis)

Unfortunately, it does only give you this for the current pixel of your map - if you've fitted it in multiple areas, it might be nice to get all of these functions back at the same time.

I've added this functionality in my "Linear Fitting" branch, which I've submitted a pull request to hyperspy. It essentially adds a "multi=True" to get the array for all pixels on the sample.

However, you probably want a quicker workaround for now - the solution is to temporarily turn all other components off, and then extract only the component you want using fitModel.as_signal().

I'll send a code example in a minute.

- Thomas

--
You received this message because you are subscribed to the Google Groups "hyperspy-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hyperspy-user...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

Thomas M. Aarholt

Post Doc, University of Oslo 

 

Mobile: +47 92855828

Email: thomas...@gmail.com

[Bernhard Schaffer]

Hi Thomas,

Yes, the “function()” was what I was trying to get at, thanks.

I think the confusion was causes by me asking two things in a single post 8-], as I was also asking for the parameter-maps. Which I get with the “as_signal”, I guess.

 

However, will the as_signal returned values be a numpy array already, or do I need to call something of that?

I don’t get the “fitModel.axis.axis“ part yet, put I need to think/check through it in more detail myself first.

 

Thanks for the link to the chat, I’ll eventually jump on to that as well.

For now, I need my “time” thinking/testing things through and keeping a “paper trail” of information to review later, so email suits me a bit better at the moment.

 

Regards and thanks again (to both of you),

Bernhard

[Thomas Aarholt]

It turns out that since I last used it, the model.as_signal() method has been updated.

To extract the component array you want, you would first fit the component in the model, and then use the following code:

component_signal = fitModel.as_signal(component_list=['Power Law']) # returns a hyperspy signal

data = component_signal.data # returns a numpy array from the hyperspy signal

I'd like to also point out that smart_fit is very old and outdated, and you should just use fitModel.fit() and fitModel.multifit() instead.

multifit fits the model in every pixel rather than just the current one you have in m.plot(). Where did you find out about smart_fit?

So a typical code where you fit some data and wanted to extract just the background fitting function, would look like the following:

import numpy as np

import hyperspy.api as hs

filename = r"C:\test.dm3"  # The r on the front allows you to use windows paths

dataSI = hs.load(filename)

print('Source data:', dataSI )

dataSI.add_elements(('Ti', 'N', 'O'))

fitModel = dataSI.create_model()    # no LowLoss

print('Fit model:')

print(fitModel.components)

fitModel.fit()

fitModel.plot()

# Now take a look at the model before continuing, make sure looks ok

fitModel.assign_current_values_to_all() # This puts the result of the fit() as starting values for the multifit in all pixels

fitModel.multifit()

fitModel.plot()

# Now use arrow keys to wander around on the plot and inspect that it looks good

# For the function, we need to know the energy-axis (x-axis) of the eels signal

energy_axis = fitModel.axis.axis # this is a "convenience function" which I shouldn't have told you about :D

# You should just use energy_axis = fitModel.axes_manager.signal_axes[0].axis

background = fitModel['PowerLaw'].function(energy_axis) # this is a calibrated hyperspy signal

back_data = background.data # this is a numpy array

Best,

Thomas

[Bernhard Schaffer]

Hi Thomas,

Thanks a lot. I’ll digest that info with some trial & testing now and will eventually bug you with some more questions, if I may.

As for the smartfit(), I’m essentially piecing all information I get together, this particular from

 

http://hyperspy.org/hyperspy-doc/current/user_guide/eels.html#eels-curve-fitting

 

(It was also a bit unfortunate that I started with that section first and only then realized that a lot of what I need is in the “Model fitting” section before that.)

 

Bernhard

 

PS: In case if you’re wondering: I’m currently working (for fun & personal interest) on an old idea of mine, which I so far never had the chance trying: Calling out to Python scripts from DigitalMicrograph.

And at the same time, I want to compare some of the functionality and results between the programs. I’ve been tiptoeing around Python and HyperSpy for a while (-no longer doing scientific data analysis in application myself-), so I’m glad I finally got things started at least…

[Thomas Aarholt]

Sounds very cool! Good luck with that! Very nice to hear that you're interested in hyperspy :) 

We'll take a look at the docs and see if we can update things a bit! 

Best, 

Thomas 

[Francisco]

Thomas, why do you think that smart_fit is "very old and outdated"? It should be as useful as ever to fit EELS core loss data with fine structure.

[Thomas Aarholt]

Apologies, I am completely wrong in this regard. I thought smart_fit was a legacy fitting function for some reason.

When Bernhard mentioned it, I googled "smart fit hyperspy" (without underscores), and on the "model fitting" page there was no mention of smart fit. Therefore I assumed it was some outdated function. 

Reading the EELS docs, perhaps they could do with a better explanation of what the "smart_fit" actually does! 

I will definitely be taking a closer look at it! :) 

Den tir. 24. jul. 2018, 10:35 skrev Francisco <frande...@gmail.com>:

Thomas, why do you think that smart_fit is "very old and outdated"? It should be as useful as ever to fit EELS core loss data with fine structure.

[Bernhard Schaffer]

Hi,

I’ve made a bit more progress with the whole thing yesterday evening, but I stumbled about something I don’t know is either a bug or something I still don’t know about. Maybe you can enlighten me :c)

Following previous suggestions, I’ve cobbled together the following script (run on a 6x6x2048 EELS SI)

 

import numpy as np

import hyperspy.api as hs

filename = "C:\\EELS_1.0nA_CoreLoss_binned2.dm3" # Test-data 6x6x2048

dataSI = hs.load(filename)

print('Source data:', dataSI )

dataSI.add_elements(('Ti', 'N', 'O'))

fitModel = dataSI.create_model() # no LowLoss

print('Fit model:')

print(fitModel.components)

fitModel.smart_fit()

fitModel.quantify()

print('Fit model parameters:')

fitModel.print_current_values()

print('PERFORMING SI FIT. BE VERY PATIENT:')

fitModel.multifit(  show_progressbar=False )

 

# iterate over fitted parameters (info only)

for comp in fitModel:

    print( comp.trait_get('name').get('name') )

    for par in comp.free_parameters:

        print( "-",par.trait_get('name').get('name') )

 

print( 'Data type:', dataArray.dtype )

print( dataArray )

 

This gives me the “map” as float64 numpy array. However, something is fishy about the array.

Further processing of that array in a pipeline tested with other float64 arrays comes up with nonsensical numbers.

However, if I do a simple numpy-array copy and use the copied array, all works fine again.

 

print( 'Data type:', dataArray.dtype )

print( dataArray )

dataArray2 = dataArray.copy()

print( 'Data type:', dataArray2.dtype )

print( dataArray2 )

 

The two arrays seem “identical” on the surface of it, but I wouldn’t know how to compare them on deeper level or to understand what’s “wrong” with the original.

Any ideas?

The actual print-output (showing that there are no stupid NaNs etc in the array)

 

PERFORMING SI FIT. BE VERY PATIENT:

Data type: float64

[[-0.40869008  3.0840393   4.70018472  6.3182053   6.00980115  0.31944671]

[-0.36471608  2.78945432  2.48835507  7.07445396 11.31741117  2.40420006]

[ 0.08022761  3.60850085  2.83486668  6.99945442 12.28993755  2.70031988]

[ 0.40889855  4.07521977  2.36466091  8.88673975 12.83003037  2.4048956 ]

[-0.29364582  3.02107615  2.68336018  6.78201372 13.50250115  2.29713747]

[-0.33466722  3.1132138   5.0523999   6.68046332  6.53065061 -0.33731116]]

Data type: float64

[[-0.40869008  3.0840393   4.70018472  6.3182053   6.00980115  0.31944671]

[-0.36471608  2.78945432  2.48835507  7.07445396 11.31741117  2.40420006]

[ 0.08022761  3.60850085  2.83486668  6.99945442 12.28993755  2.70031988]

[ 0.40889855  4.07521977  2.36466091  8.88673975 12.83003037  2.4048956 ]

[-0.29364582  3.02107615  2.68336018  6.78201372 13.50250115  2.29713747]

[-0.33466722  3.1132138   5.0523999   6.68046332  6.53065061 -0.33731116]]

 

BTW, I did this initially with Numpy 1.14.0 but have then upgraded to 1.15.0 – still the same.

 

Best regards,

  Bernhard Schaffer

 

PS: If you have any other comments on my scripts, i.e. things that could be done much nicer than I did, please let me know! Still learning…

 

From: hypersp...@googlegroups.com <hypersp...@googlegroups.com> On Behalf Of Thomas Aarholt
Sent: 24 July 2018 09:51
To: hyperspy-users <hypersp...@googlegroups.com>
Subject: Re: [hyperspy-users] Re: Model fitting, getting fit-models as array

 

Apologies, I am completely wrong in this regard. I thought smart_fit was a legacy fitting function for some reason.

[Magnus Nord]

Hey Bernhard,

Nice to see you here :)

Could you include the EEL spectrum file as well, it makes it easier to debug/replicate the script results.

Also, in case you don't know about them, HyperSpy has several Jupyter Notebook demos here: https://github.com/hyperspy/hyperspy-demos , specifically the EELS one might be of interest to you: https://github.com/hyperspy/hyperspy-demos/blob/master/electron_microscopy/EELS/EELS_analysis.ipynb

Magnus

[Bernhard Schaffer]

Sure, it’s not a secret. Does it work as simple email attachment (now attached) or do I have to upload it somewhere?

Bernhard

[Francisco]

It is not easy to tell what the issue is because, although I missed something, we don't know how you get the faulty dataArray nor what the fault is. Could you provide more details?

I second Thomas on his gitter suggestion: this thread is getting very long and sharing code in Google groups is less than ideal. Bernhard, would you like to give gitter a try?

Francisco