Magnetic reflectivity py #131
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Empty file.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,265 @@ | ||
| import numba | ||
| import math | ||
|
|
||
| Z_EPS = 1e-6 | ||
|
|
||
| ALIGN_MAGNETIC_SIG = 'i4(f8[:], f8[:], f8[:], f8[:], f8[:], f8[:], f8[:], f8[:], f8[:,:])' | ||
|
|
||
|
|
||
| # @numba.njit(ALIGN_MAGNETIC_SIG, parallel=False, cache=True) | ||
| @numba.njit(cache=True) | ||
|
There was a problem hiding this comment. Why no function signature? If it is because the inputs are not of the correct time and we are relying on conversion during call, then I would prefer to trigger an error here and force the caller to produce the correct types. There was a problem hiding this comment. for some of these there was explicit polymorphism encoded in the wrapper, and it seemed like a reasonable policy in numba to address that by not specifying a signature. I did some testing and was not able to see any performance difference between signature/no signature. For this function I suspect the usage of the variables is pretty unambiguous to numba. |
||
| def align_magnetic(d, sigma, rho, irho, dM, sigmaM, rhoM, thetaM, output_flat): | ||
| # ignoring thickness d on the first and last layers | ||
| # ignoring interface width sigma on the last layer | ||
| nlayers = len(d) | ||
| nlayersM = len(dM) | ||
| noutput = nlayers + nlayersM | ||
|
|
||
| # making sure there are at least two layers | ||
| if nlayers < 2 or nlayersM < 2: | ||
| raise ValueError("only works with more than one layer") | ||
|
|
||
| output = output_flat | ||
| magnetic = 0 # current magnetic layer index | ||
| nuclear = 0 # current nuclear layer index | ||
| z = 0.0 # current interface depth | ||
| next_z = 0.0 # next nuclear interface | ||
| next_zM = 0.0 # next magnetic interface | ||
| # int active = 3# active interfaces, active&0x1 for nuclear, active&0x2 for magnetic | ||
|
|
||
| k = 0 # current output layer index | ||
| while True: | ||
| # repeat over all nuclear/magnetic layers | ||
| if (k == noutput): | ||
| return -1 # exceeds capacity of output | ||
|
|
||
| # printf("%d: %d %d %g %g %g\n", k, nuclear, magnetic, z, next_z, next_zM) | ||
| # printf("%g %g %g %g\n", rho[nuclear], irho[nuclear], rhoM[magnetic], thetaM[magnetic]) | ||
| # printf("%g %g %g %g\n", d[nuclear], sigma[nuclear], dM[magnetic], sigmaM[magnetic]) | ||
|
|
||
| # Set the scattering strength using the current parameters | ||
| output[k][2] = rho[nuclear] | ||
| output[k][3] = irho[nuclear] | ||
| output[k][4] = rhoM[magnetic] | ||
| output[k][5] = thetaM[magnetic] | ||
|
|
||
| # Check if we are at the last layer for both nuclear and magnetic | ||
| # If so set thickness and interface width to zero. We are doing a | ||
| # center of the loop exit in order to make sure that the final layer | ||
| # is added. | ||
|
|
||
| if (magnetic == nlayersM-1 and nuclear == nlayers-1): | ||
| output[k][0] = 0. | ||
| output[k][1] = 0. | ||
| k += 1 | ||
| break | ||
|
|
||
| # Determine if we are adding the nuclear or the magnetic interface next, | ||
| # or possibly both. The order of the conditions is important. | ||
| # | ||
| # Note: the final value for next_z/next_zM is not defined. Rather than | ||
| # checking if we are on the last layer we simply add the value of the | ||
| # last thickness to z, which may be 0, nan, inf, or anything else. This | ||
| # doesn't affect the algorithm since we don't look at next_z when we are | ||
| # on the final nuclear layer or next_zM when we are on the final magnetic | ||
| # layer. | ||
| # | ||
| # Note: averaging nearly aligned interfaces can lead to negative thickness | ||
| # Consider nuc = [1-a, 0, 1] and mag = [1+a, 1, 1] for 2a < Z_EPS. | ||
| # On the first step we set next_z to 1-a, next_zM to 1+a and z to the | ||
| # average of 1-a and 1+a, which is 1. On the second step next_z is | ||
| # still 1-a, so the thickness next_z - z = -a. Since a is tiny we can just | ||
| # pretend that -a == zero by setting thickness to fmax(next_z - z, 0.0). | ||
|
|
||
| if (nuclear == nlayers-1): | ||
| # No more nuclear layers... play out the remaining magnetic layers. | ||
| output[k][0] = max(next_zM - z, 0.0) | ||
| output[k][1] = sigmaM[magnetic] | ||
| magnetic += 1 | ||
| next_zM += dM[magnetic] | ||
| elif (magnetic == nlayersM-1): | ||
| # No more magnetic layers... play out the remaining nuclear layers. | ||
| output[k][0] = max(next_z - z, 0.0) | ||
| output[k][1] = sigma[nuclear] | ||
| nuclear += 1 | ||
| next_z += d[nuclear] | ||
| elif (math.fabs(next_z - next_zM) < Z_EPS and math.fabs(sigma[nuclear]-sigmaM[magnetic]) < Z_EPS): | ||
| # Matching nuclear/magnetic boundary, with almost identical interfaces. | ||
| # Increment both nuclear and magnetic layers. | ||
| output[k][0] = max(0.5*(next_z + next_zM) - z, 0.0) | ||
| output[k][1] = 0.5*(sigma[nuclear] + sigmaM[magnetic]) | ||
| nuclear += 1 | ||
| next_z += d[nuclear] | ||
| magnetic += 1 | ||
| next_zM += dM[magnetic] | ||
| elif (next_zM < next_z): | ||
| # Magnetic boundary comes before nuclear boundary, so increment magnetic. | ||
| output[k][0] = max(next_zM - z, 0.0) | ||
| output[k][1] = sigmaM[magnetic] | ||
| magnetic += 1 | ||
| next_zM += dM[magnetic] | ||
| else: | ||
| # Nuclear boundary comes before magnetic boundary | ||
| # OR nuclear and magnetic boundaries match but interfaces are different. | ||
| # so increment nuclear. | ||
| output[k][0] = max(next_z - z, 0.0) | ||
| output[k][1] = sigma[nuclear] | ||
| nuclear += 1 | ||
| next_z += d[nuclear] | ||
|
|
||
| z += output[k][0] | ||
| k += 1 | ||
|
|
||
| return k | ||
|
|
||
|
|
||
| CONTRACT_MAG_SIG = 'i4(f8[:], f8[:], f8[:], f8[:], f8[:], f8[:], f8)' | ||
|
|
||
|
|
||
| # @numba.njit(CONTRACT_MAG_SIG, parallel=False, cache=True) | ||
| @numba.njit(cache=True) | ||
| def contract_mag(d, sigma, rho, irho, rhoM, thetaM, dA): | ||
|
There was a problem hiding this comment. I believe this code is broken and the magnetic folk always run with dA=0. I've added a ticket. |
||
| n = len(d) | ||
| i = newi = 1 # /* Skip the substrate */ | ||
| while (i < n): | ||
|
|
||
| # /* Get ready for the next layer */ | ||
| # /* Accumulation of the first row happens in the inner loop */ | ||
| dz = weighted_dz = 0 | ||
| rhoarea = irhoarea = rhoMarea = thetaMarea = 0.0 | ||
| rholo = rhohi = rho[i] | ||
| irholo = irhohi = irho[i] | ||
| maglo = maghi = rhoM[i] * math.cos(thetaM[i] * math.pi / 180.0) | ||
|
|
||
| # /* Accumulate slices into layer */ | ||
| while True: | ||
| # /* Accumulate next slice */ | ||
| dz += d[i] | ||
| rhoarea += d[i] * rho[i] | ||
| irhoarea += d[i] * irho[i] | ||
|
|
||
| # /* Weight the magnetic signal by the in -plane contribution | ||
| # * when accumulating rhoM and thetaM. */ | ||
| weight = math.cos(thetaM[i]*math.pi/180.) | ||
| mag = rhoM[i]*weight | ||
| rhoMarea += d[i]*rhoM[i]*weight | ||
| thetaMarea += d[i]*thetaM[i]*weight | ||
| weighted_dz += d[i]*weight | ||
|
|
||
| # /* If no more slices or sigma != 0, break immediately */ | ||
| i += 1 | ||
| if (i == n or sigma[i-1] != 0.): | ||
| break | ||
|
|
||
| # /* If next slice exceeds limit then break */ | ||
| if (rho[i] < rholo): | ||
| rholo = rho[i] | ||
| if (rho[i] > rhohi): | ||
| rhohi = rho[i] | ||
| if ((rhohi-rholo)*(dz+d[i]) > dA): | ||
| break | ||
|
|
||
| if (irho[i] < irholo): | ||
| irholo = irho[i] | ||
| if (irho[i] > irhohi): | ||
| irhohi = irho[i] | ||
| if ((irhohi-irholo)*(dz+d[i]) > dA): | ||
| break | ||
|
|
||
| if (mag < maglo): | ||
| maglo = mag | ||
| if (mag > maghi): | ||
| maghi = mag | ||
| if ((maghi-maglo)*(dz+d[i]) > dA): | ||
| break | ||
|
|
||
| # /* Save the layer */ | ||
| assert(newi < n) | ||
| d[newi] = dz | ||
| if (i == n): | ||
| # /* Last layer uses surface values */ | ||
| rho[newi] = rho[n-1] | ||
| irho[newi] = irho[n-1] | ||
| rhoM[newi] = rhoM[n-1] | ||
| thetaM[newi] = thetaM[n-1] | ||
| # /* No interface for final layer */ | ||
| else: | ||
| # /* Middle layers uses average values */ | ||
| rho[newi] = rhoarea / dz | ||
| irho[newi] = irhoarea / dz | ||
| rhoM[newi] = rhoMarea / weighted_dz | ||
| thetaM[newi] = thetaMarea / weighted_dz | ||
| sigma[newi] = sigma[i-1] | ||
| # /* First layer uses substrate values */ | ||
| newi += 1 | ||
|
|
||
| return newi | ||
|
|
||
|
|
||
| CONTRACT_BY_AREA_SIG = 'i4(f8[:], f8[:], f8[:], f8[:], f8)' | ||
|
|
||
|
|
||
| # @numba.njit(CONTRACT_BY_AREA_SIG, parallel=False, cache=True) | ||
| @numba.njit(cache=True) | ||
| def contract_by_area(d, sigma, rho, irho, dA): | ||
| n = len(d) | ||
| i = newi = 1 # /* Skip the substrate */ | ||
| while (i < n): | ||
|
|
||
| # /* Get ready for the next layer */ | ||
| # /* Accumulation of the first row happens in the inner loop */ | ||
| dz = rhoarea = irhoarea = 0.0 | ||
| rholo = rhohi = rho[i] | ||
| irholo = irhohi = irho[i] | ||
|
|
||
| # /* Accumulate slices into layer */ | ||
| while True: | ||
| # /* Accumulate next slice */ | ||
| dz += d[i] | ||
| rhoarea += d[i]*rho[i] | ||
| irhoarea += d[i]*irho[i] | ||
|
|
||
| # /* If no more slices or sigma != 0, break immediately */ | ||
| i += 1 | ||
| if (i == n or sigma[i-1] != 0.): | ||
| break | ||
|
|
||
| # /* If next slice won't fit, break */ | ||
| if (rho[i] < rholo): | ||
| rholo = rho[i] | ||
| if (rho[i] > rhohi): | ||
| rhohi = rho[i] | ||
| if ((rhohi-rholo)*(dz+d[i]) > dA): | ||
| break | ||
|
|
||
| if (irho[i] < irholo): | ||
| irholo = irho[i] | ||
| if (irho[i] > irhohi): | ||
| irhohi = irho[i] | ||
| if ((irhohi-irholo)*(dz+d[i]) > dA): | ||
| break | ||
|
|
||
| # /* dz is only going to be zero if there is a forced break due to | ||
| # * sigma, or if we are accumulating a substrate. In either case, | ||
| # * we want to accumulate the zero length layer | ||
| # */ | ||
| # /* if (dz == 0) continue; */ | ||
|
|
||
| # /* Save the layer */ | ||
| assert(newi < n) | ||
| d[newi] = dz | ||
| if (i == n): | ||
| # /* printf("contract: adding final sld at %d\n",newi); */ | ||
| # /* Last layer uses surface values */ | ||
| rho[newi] = rho[n-1] | ||
| irho[newi] = irho[n-1] | ||
| # /* No interface for final layer */ | ||
| else: | ||
| # /* Middle layers uses average values */ | ||
| rho[newi] = rhoarea / dz | ||
| irho[newi] = irhoarea / dz | ||
| sigma[newi] = sigma[i-1] | ||
| # /* First layer uses substrate values */ | ||
| newi += 1 | ||
|
|
||
| return newi | ||
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is 3.8 now required? Numpy is sitting at 3.7 as the oldest version, so ideally we would match that (though not important).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think I was using 3.8 as the target for the dataclasses branch, but I don't remember why. I don't think it's important - 3.7 is fine.