Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Be lens calc follow up #8

Merged
merged 21 commits into from
Nov 2, 2020

Conversation

cristinasewell
Copy link
Contributor

The lens_set file has been changed from using numpy to using json to save and read the data from the file.
The plan_set function has been added
A new lens_transmission function has been added too

Motivation and Context

How Has This Been Tested?

It is hard to test this, tested against the old function and seems to work the same, but more testing needs to be done....

Where Has This Been Documented?

Docstring

logger.info('\n %s', resstring)


def lens_transmission(r, fwhm, n=1, energy=None, id="Be", lens_thicknes=None):
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the old code we have ID='IF1' - is the ID here the material? that's what is seems like but IF1 is not a valid material for the current get_att_len function
https://github.com/pcdshub/old-hutch-configs/blob/18c7769d66a647ab558ad7ebc34d019fc98fbf48/cxi/blutil/calc.py#L1189

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Small horror story here: the old hutch python actually used a modified version of periodictable. So there very well could have been an added entry named IF1.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you think we'll need to support it?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I hope not. I'll quickly go looking for this illicit version to see what it does.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, a grep for IF1 is turning up empty

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ohh noo 😭

Copy link
Contributor Author

@cristinasewell cristinasewell Oct 21, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, I am going to add the aliases here then for the chemical formulas that they have. But in order to make it work with xrayDB I have to provide the density too, so I'll have to bring in the Density as well...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

added this in - 1fcd6fb - but I don't like it, should we have those densities and material in a separate file?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Keeping them in code is fine for now, in my opinion. Then yaml or JSON down the line.



# with good_sets = [False, True, True, False, False, True, False, False, False]
def test_plan_set():
Copy link
Contributor Author

@cristinasewell cristinasewell Oct 21, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is hard to make the plan_set generate some True for good_sets to be able to test the results.... I've been testing it with dummy good_sets against the old function and I've been getting the similar results, but obviously I have to find some good parameters to test this

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems good enough to me. The values provided to that function don't seem to change the control flow, meaning testing arbitrary-ish values against the old implementation is likely sufficient.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with Ken, this is good

lens_thicknes = lens_thicknes or APEX_DISTANCE
id_material = check_id(id_material)
# TODO: what should I do with the energy here
# energy = getE(energy=energy, correctEv=True)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this function requires a PV, should we have that in the calculations or is this up to to user to pass in the correct energy when using the classes?
https://github.com/pcdshub/old-hutch-configs/blob/18c7769d66a647ab558ad7ebc34d019fc98fbf48/cxi/blutil/calc.py#L10-L40

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There shouldn't be anything that requires a PV in this codebase - values only.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The previous codebase checked a PV if you didn't provide an energy value. I'm not sure exactly how/where to implement this behavior... perhaps not until we get to pcdsdevices and you need to move to a specific beam size?

Generally when you're planning an experiment the current x-ray energy does not match what it will be during your beam time.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ZLLentz what's the issue in requiring the user to be explicit about the energy as an input to these functions? Or rather, why do you think that functionality from the old codebase need to be included here?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When you call these functions you should call them with the energy you expect to have during the beam time. There may be some cases during a beam time when you want to call a function with "current beam energy". These are mostly localized around "I want to move my calculation motor without repeatedly updating the beam energy parameter"

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see - perhaps our LCLS object should have an easily toggleable subscription that calls the new equivalent of setE for setting the default energy when none is supplied by kwarg?

@klauer when you get a chance can you please elaborate on this^?

Will it look similar to this: (in the LCSL object)

photon_ev_hxr = Cpt(EpicsSignalRO, 'SIOC:SYS0:ML00:AO627', kind='hinted',
                        doc='Photon eV HXR [eV]')

    @photon_ev_hxr.sub_value
    def _energy_changed(self, *args, value, ** kwargs):
        self.set_energy(value)
        print(f'----VALUE CHANGED -------{value}')

    def set_energy(self, value):
        global _photon_energy
        _photon_energy = value

    def get_energy(self):
        return _photon_energy

And then in the LensStack object I would get this value from the LCLS object?

lcls = LCLS()
lcls.get_energy()

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was thinking something along the lines of:

photon_ev_hxr.subscribe(_set_current_energy)

def get_energy():
   if user_specified_energy is not None:
       return user_specified_energy
    return _current_energy

def other_calculation_function(, ..., energy=None):
    if energy is None:
        energy = get_energy()

It's a lot of customization of one parameter - but it's apparently an important one to be flexible with here.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm still a bit confused here - i'm so sorry.
Here is what i understand:
The photon_ev_hxr.subscribe(_set_current_energy) is in the beam_stats LCLS object. Here is where I will have access to the current_energy - but I don't need to know here about the user_specified_energy - or do i? - i'm not doing any calculations here.
I would need to know about the user specified energy in the LensStack and in the calculations...

I need to think about this a bit more, i might be overcomplicating things here

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The user_specified_energy is the user-provided energy that acts as an override for the current energy provided by the LCLS subscription callback.

Maybe I'm missing something here, too - not sure where the confusion is coming in. (@ZLLentz?)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like (3) - let's get this rolling and in the hands of those that can use it 👍

Sounds good, at least if people start using it we'll get some feedback (maybe) from them and hopefully some of the questions/concerns will be answered...

>>> list_of_sets = [[3, 0.0001, 1, 0.0002],
[1, 0.0001, 1, 0.0003, 1, 0.0005],
[2, 0.0001, 1, 0.0005]]
>>> set_lens_set_to_file(list_of_sets, ../path/to/lens_set)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
>>> set_lens_set_to_file(list_of_sets, ../path/to/lens_set)
>>> set_lens_set_to_file(list_of_sets, '../path/to/lens_set')

if set_number_top_to_bot not in range(1, sets.shape[0]):
with open(filename) as lens_file:
sets = json.loads(lens_file.read())
if set_number_top_to_bot not in range(1, len(sets)):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This and other errors in the function should probably be bumped up to exceptions (ValueError) as they'd break the calling function by returning None

logger.info('\n %s', resstring)


def lens_transmission(r, fwhm, n=1, energy=None, id="Be", lens_thicknes=None):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Keeping them in code is fine for now, in my opinion. Then yaml or JSON down the line.



# with good_sets = [False, True, True, False, False, True, False, False, False]
def test_plan_set():
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems good enough to me. The values provided to that function don't seem to change the control flow, meaning testing arbitrary-ish values against the old implementation is likely sufficient.

)
thickness[i, j] = (x[i] ** 2 + y[j] ** 2) / r + n * lens_thicknes
d = density.get(id_material)
att_length = get_att_len(energy, id_material, d)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need to revise this one, I'm using the get_att_len here (and tested it with get_att_len in the old code too), but after doing some research I realized that this is not the same as attenuationLength function from https://github.com/pcdshub/old-hutch-configs/blob/18c7769d66a647ab558ad7ebc34d019fc98fbf48/xpp/prod/utilitiesCalc.py#L362
which brings with it a few more other functions.... that require periodictable...

pcdscalc/be_lens_calcs.py Outdated Show resolved Hide resolved

Parameters
----------
energy : number
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Most of the input parameters here need to have units associated with them. It's probably good to check other functions as well so that we don't get confused later.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, i added to some of these functions in 279f107 but not to all of them (because I am not 100 % sure)
I am almost sure that the radius that we're using from the lens_set and in most of the functions are in meters, so i can add that too but I am not 100 % convinced.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok we need to figure these out then, if the wrong units are used then the calcs will not work as expected

Copy link
Contributor Author

@cristinasewell cristinasewell Oct 28, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah I agree - i'm going to go back and try to figure out all of them. Based on the old curren_set file we have these:

[(3, 0.0001, 1, 0.0002),
 (1, 0.0001, 1, 0.0003, 1, 0.0005),
 (2, 0.0001, 1, 0.0005)]

So it make sense that they would be in meters, because we have the LENS_RADII = [50e-6, 100e-6, 200e-6, 300e-6, 500e-6, 1000e-6, 1500e-6, 2000e-6, 3000e-6] which are in um
I need to double check here though ...

@cristinasewell cristinasewell marked this pull request as ready for review October 27, 2020 17:34
pcdscalc/be_lens_calcs.py Outdated Show resolved Hide resolved


def plan_set(energy, z_offset, z_range, beam_size_unfocused, size_horizontal,
size_vertical=None, exclude=[], max_tot_number_of_lenses=25,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It doesn't matter here, but sometimes you can run into trouble if you pass a mutable object like a list as the default argument to a function. This is because if you modify this list in the function call, the modified list will be used as the default value in the next call of the function, rather than the original list.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, thank you for pointing this out. I can make it None and inside the function do something like this:
exclude = exclude or []?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, that is the typical way to handle this, using an immutable sentinel like None

@@ -42,7 +45,7 @@ def configure_defaults(fwhm_unfocused=None, disk_thickness=None,
Parameters
-----------
fwhm_unfocused : float, optional
Full width at half maximum unfocused
Full width at half maximum unfocused.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion to match the unit in the other docstrings

Suggested change
Full width at half maximum unfocused.
Full width at half maximum unfocused in meters.

size_fwhm : float
distance : float
n_max : int, optional
max_each : int, optional
lens_radii : list, optional
fwhm_unfocused : float, optional
This is about 400 microns at XPP. Defaults to 500e-6.
This is about 400 microns at XPP. Defaults to 500e-6. (in meters)
eff_rad0 : float, optional
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this docstring is missing some descriptions too, though they are the same parameters as other functions in the file

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not sure here what eff_rad0 is - thinking that could be Effective Radius of curvature? but not sure
any ideas?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also I am not sure of these guys' units:

# Disk Thickness
DISK_THICKNESS = 1.0e-3
# Apex of the lens
APEX_DISTANCE = 30e-6

again i am assuming that they have to be meters - will come back to this one in a sec

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

1mm and 30µm sound like reasonable values for those, so I'd guess meters is correct.

@ZLLentz
Copy link
Member

ZLLentz commented Oct 28, 2020

I want to say that this is really close

@cristinasewell
Copy link
Contributor Author

cristinasewell commented Oct 28, 2020

I still don't feel comfortable for the Lens Transmission in the plan_set function, but I think I will know better when Matt gets back to me. I think we could be using the calc_trans_lens_set or calc_trans_for_single_lens - but we can adjust this later when we have more info..?

@ZLLentz
Copy link
Member

ZLLentz commented Oct 28, 2020

Yeah, we have a few options:

  1. Be patient for a response from the experts
  2. Merge, but don't tag, and fix discrepancies in a follow-up PR
  3. Merge, tag, and start including the functions in hutch-python/pcdsdevices as necessary, with the knowledge we may need to revise the base calculations later.

I'm leaning toward merging- otherwise this continues to get more difficult to review as it grows in scope.

@klauer
Copy link
Contributor

klauer commented Oct 28, 2020

I like (3) - let's get this rolling and in the hands of those that can use it 👍

@cristinasewell
Copy link
Contributor Author

I thought I added a comment earlier today here but I don't see it - I don't know what happen to it.
Anyway, I agree with you about option 3 - when people will start using these calculations we'll probably get feedback and then we can adjust things if needed.


import numpy as np
import xraydb as xdb

logger = logging.getLogger(__name__)
# global variable to be help with testing the plan_set function
_plan_set_test_res = None
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this was a good solution and doesn't need a revision. Another way to do this is to have testable helper functions with whatever return values we want, and use the helpers in the functions without good return values.

err_msg = ('Provided invalid path for lens set file: %s',
lens_file_path)
logger.error(err_msg)
raise FileNotFoundError(err_msg)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The error message here is not substituted correctly, it comes out as a tuple

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ohh i see... interesting, thanks for catching this.

Before I commit things though, would something like this be acceptable:

 err_msg = f'Provided invalid path for lens set file: {lens_file_path}'
        logger.error(err_msg)
        raise FileNotFoundError(err_msg)

I know we don't like using the f-Strings for logging.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, personally I've given up on doing the python logging "properly" with the %s substitutions because I find myself wanting to re-use the strings... I'm fine with f-strings and with .format

My understanding is that when you do logger.error('Bad thing: %s', bad_thing) and it happens many times you can filter on Bad thing: %s to catch all logging messages like this... but unless the log message happens a lot I think this is worth overlooking for formatting clarity/consistency.

@klauer might have an opinion too

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

An opinion? Of course!

Use:

  • logging.error('%s', value) if the string isn't going to be reused (generally speaking, I think this is the majority of the cases)
  • f'{value}' for pretty much everything else.

.format() most of the time ends up not being as nice as an f-string, as it requires you to visually match up the arguments instead of seeing how they're going to be interpolated directly in the string. If you find yourself justifying the use of .format() over an f-string because you need to do some calculation/cleanup before formatting, I'd recommend just making a temporary variable and using an f-string with that.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you, this should be fixed in ef27fed

beam_size_unfocused : float
Radial size of x-ray beam before focusing in meters
size_horiz
size_vert
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These two are missing parameter info and their names don't match the signature

Copy link
Contributor Author

@cristinasewell cristinasewell Nov 2, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I honestly don't know what these two fields are... I was thinking that they could be the Field of View horizontal, Field of View Vertical - which are just numbers, in meters or millimeter

But looking at the code, the np.max is taking size_horizontal as the a - array like, and the size_vertical as the axis - so this is confusing (unless they meant to use np.maximum?? - but don't think so)

  if size_vertical is not None:
        size_calc_1 = np.max(size_horizontal, size_vertical)
    else:
        size_calc_1 = size_horizontal

That's why I had it empty for now until I get more feedback from the scientists.... unless you guys know what they are supposed to be?

All i can say for now i guess is that the way we use it, we can conclude that:

size_horizontal : array_like
size_vertical : int or tuple of ints

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fair enough! Thanks for looking into the types

@ZLLentz
Copy link
Member

ZLLentz commented Nov 2, 2020

I've submitted my last few nitpicks- I'd be happy to merge/tag

@ZLLentz
Copy link
Member

ZLLentz commented Nov 2, 2020

Ok let's put this to rest

@ZLLentz ZLLentz merged commit 35c023d into pcdshub:master Nov 2, 2020
@klauer klauer mentioned this pull request Nov 3, 2020
2 tasks
@cristinasewell cristinasewell mentioned this pull request Dec 4, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants