text:
Energy levels (Ch. 13) Nov. 13, 2017

Reading for atomic energy level (Laser Fundamentals, William & Silvast)

• Schrodinger equation: Model particle behavior$\text{ }-\frac{{\text{ℏ}}^{2}}{2m}{\nabla }^{2}\psi \left(\stackrel{\to }{r},t\right)\text{\hspace{0.17em}}+\text{\hspace{0.17em}}V\left(\stackrel{\to }{r},t\right)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}-j\text{ℏ}\partial \frac{\psi \left(\stackrel{\to }{r},t\right)}{\partial t}$ where$\text{ }V\left(\stackrel{\to }{r},t\right)$ is potential set up by the coulomb force from nucleus,$\text{ }\psi$ complex wave function, m mass of electron.
• Time independent Schrodinger equation (TISE): For particle in stationary state (stay in orbit),$\text{ }\psi \left(\stackrel{\to }{r},t\right)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\psi \left(\stackrel{\to }{r}\right){e}^{jEt/\text{ℏ}}$ and the equation becomes similar to Helmholtz equation$\text{ }-\frac{{\text{ℏ}}^{2}}{2m}{\nabla }^{2}\psi \text{\hspace{0.17em}}+\text{\hspace{0.17em}}V\psi \text{\hspace{0.17em}}=\text{\hspace{0.17em}}E\psi$
where E is the eigenvalue of the equation corresponding to the energy level.
Electron probability over a volume dV =$|\psi {|}^{2}dV$ and$\text{ }\int |\psi {|}^{2}dV\text{\hspace{0.17em}}=\text{\hspace{0.17em}}1$
• Atomic levels for$\text{ }{H}_{2}$ : Hydrogen is the simplest atom. Its energy level$\text{ }{E}_{n}$ can be found with Bohr’s model by equating potential energy from Coulomb force to Kinetic energy of centrifugal force from orbiting electron$\text{ }{E}_{n}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{−}{M}_{r}\frac{{e}^{4}}{\left(4\pi {\epsilon }_{o}{\right)}^{2}2{\text{ℏ}}^{2}}\text{\hspace{0.17em}}\frac{1}{{n}^{2}}$ where reduced mass$\text{ }{M}_{r}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{m{M}_{p}}{m\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{M}_{p}}$ ,$\text{ }{M}_{p}$ is the mass of proton, n is principal quantum number. Also radius of the first Bohr orbit (n=1)$\text{ }{a}_{H}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{{\epsilon }_{o}{h}^{2}}{\pi m{e}^{2}}$
Hydrogen-like ions: Ions with one outer electron, e.g.$\text{ }H{e}^{+}$ ,$\text{ }L{i}^{2+}$ , $\text{ }{E}_{n}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{−}{M}_{r}\frac{{e}^{4}}{\left(4\pi {\epsilon }_{o}{\right)}^{2}2{\text{ℏ}}^{2}}\text{\hspace{0.17em}}{\left(\frac{Z}{n}\right)}^{2}$ where Z is the number of protons in nucleus.
• Detailed energy state: Solve TISE with separation of variable$\text{ }{\psi }_{nlm}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{R}_{nl}\left(r\right){\Theta }_{lm}\left(\theta \right){\Phi }_{m}\left(\phi \right)$ where$\text{ }l\text{\hspace{0.17em}}=\text{\hspace{0.17em}}0\left(s\right),\text{\hspace{0.17em}}1\left(p\right),\text{\hspace{0.17em}}...,n-1$ is the azimuthal quantum number and$\text{ }m\text{\hspace{0.17em}}=\text{\hspace{0.17em}}0,\text{\hspace{0.17em}}±1,±2,...,\text{\hspace{0.17em}}±l$ is the magnetic quantum number.

$\text{ }{R}_{nl}\left(r\right)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\rho }^{l}{e}^{-\rho /2}{L}_{n+l}^{2l+1}\left(\rho \right)$

where where$\text{ }{L}_{n+1}^{2l+1}\left(\rho \right)$ is the Laguerre polynomial with$\text{ }\rho \text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{2r}{{a}_{H}}$$\text{ }{\Theta }_{lm}\left(\theta \right)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{P}_{l}^{|m|}\left(x\right)\frac{\left(1\text{\hspace{0.17em}}-\text{\hspace{0.17em}}{x}^{2}{\right)}^{1/\left(2|m|\right)}}{2l!}\frac{{d}^{|m|+1}}{d{x}^{|m|+1}}\left({x}^{2}\text{\hspace{0.17em}}-\text{\hspace{0.17em}}1{\right)}^{l}$ where$\text{ }{P}_{l}^{|m|}\left(x\right)$ is the Legendre polynomials with$\text{ }x\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{cos}\theta$

$\text{ }{\Phi }_{m}\left(\phi \right)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{e}^{jm\phi }$

Spin quantum number$\text{ }s\text{\hspace{0.17em}}=\text{\hspace{0.17em}}±1/2$
Notice that we show outer unfilled subshell in energy state diagram as a shorthand notation e.g. excited C with$1{s}^{2}{2}^{2}2p3s$ is referred C (3s).
• Multi-electron atoms: Need to consider interaction of electrons including 1) kinetic energy (orbiting), 2) potential energy (Coulomb force), 3) mutual electrostatic energy of electrons, 4) spin-spin correlations and 5) spin-orbit energy.
1) and 2) are included in Bohr’s model
3) and 4) are accounted for by LS coupling which is commonly used for lasing transition.
5) is considered in J-J coupling (not common).
• LS coupling:$\text{ }{\text{}}^{2S+1}{L}_{J}$ where$\text{ }\stackrel{\to }{L}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\stackrel{\to }{L}}_{1}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\stackrel{\to }{L}}_{2}$ ,$\text{ }L\text{\hspace{0.17em}}=\text{\hspace{0.17em}}|\stackrel{\to }{L}|$ ,$\text{ }\stackrel{\to }{S}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\stackrel{\to }{S}}_{1}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\stackrel{\to }{S}}_{2}$ ,$\text{ }S=\text{\hspace{0.17em}}|\stackrel{\to }{S}|$ and$\text{ }J\text{\hspace{0.17em}}=\text{\hspace{0.17em}}|\stackrel{\to }{L}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{+}\text{\hspace{0.17em}}\stackrel{\to }{S}|$ , e.g. He and He Ne.
• Molecules: Each molecules has multiple atoms, i.e. more energy states for rotation and vibration.
Diatomic rotation ($\text{ }{E}_{r}$ with rotational quantum number$\text{ }r\text{\hspace{0.17em}}=\text{\hspace{0.17em}}0,\text{\hspace{0.17em}}1,\text{\hspace{0.17em}},2...$ ) and vibration ($\text{ }{E}_{v}$ with vibrational quantum number$\text{ }v\text{\hspace{0.17em}}=\text{\hspace{0.17em}}0,\text{\hspace{0.17em}}1,\text{\hspace{0.17em}},2...$ )
Triatomic molecule, e.g. $\text{ }C{O}_{2}$ with 3 types of vibration -- asymmetric stretch, symmetric stretch and bending.
Dye molecules are large and complex organic material$\text{ }\to$ more complicated manifolds with wide range of possible lasing freq.
• Solid-state: Atoms/molecules are closer and interact with lattice (a periodic structure). Low energy levels remain sharp while the outer levels smear and form bands as atoms getting closer.
Lower band is called valence band with energy ($\text{ }{E}_{1}$ or$\text{ }{E}_{v}$ ). Upper band is called conduction band with energy ($\text{ }{E}_{2}$ or$\text{ }{E}_{c}$ ). A forbidden band is in between with bandgap energy$\text{ }{E}_{g}\text{\hspace{0.17em}}={E}_{2}\text{\hspace{0.17em}}-\text{\hspace{0.17em}}{E}_{1}$ Three types of solid -- metal (partial occupied conduction band), semiconductor (empty conduction band at 0K) and insulator (like semiconductor but with large$\text{ }{E}_{g}$ ).
Solid-state laser materials -- doped dielectric media, a host medium (transparent dielectric) with dopant (transition-metal, lanthanide ions), e.g. Ruby$\text{ }C{r}^{3+}$ :$\text{ }A{l}_{2}{O}_{3}$ , Ti Sapphire, Neodymium:YAG$\text{ }N{d}^{3+}:{Y}_{3}A{l}_{5}{O}_{2}$ ($1.064\mu m$ ), Neodymium:glass ($1.053\mu m$ ).
Semiconductors -- Similar to other solids with conduction band, valence band and bandgap energy but have 2 types of charge carriers (electrons and holes).
Band is not flat (more like parabolic) in spatial freq domain owing to periodic structure of lattice.
Direct bandgap (e.g. GaAs) is more favorable as lasing media than indirect bandgap (e.g. Si).
Man-made structures -- building heterojunctions with semiconductors having matched lattice constant but different$\text{ }{E}_{g}$ quantum well for 1D electron confinement, quantum wire for 2D electron confinement and quantum dot for 3D electron confinement.
• Occupation of energy level: Thermal distribution of particle can be modeled by Boltzmann distribution.
Probability$\text{ }P\left({E}_{m}\right)$ of finding atom at energy level$\text{ }{E}_{m}$$\text{ }P\left({E}_{m}\right)\propto \text{exp}\left(-{E}_{m}/\left(kT\right)\right)$ where m = 1, 2, 3, ...
Notice that we require$\text{ }\underset{m}{\overset{\text{}}{\Sigma }}P\left({E}_{m}\right)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}1$ and is not an integration of continuous function.
For laser media,$\text{ }\frac{{N}_{2}}{{N}_{1}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{{g}_{2}}{{g}_{1}}\text{exp}\left(-\frac{{E}_{2}\text{\hspace{0.17em}}-\text{\hspace{0.17em}}{E}_{1}}{kT}\right)$ where$\text{ }{g}_{1,2}$ are degeneracies of energy levels 1 and 2.
Semiconductors follow Fermi-Dirac distribution$\text{ }f\left(E\right)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{1}{{e}^{\left(E\text{\hspace{0.17em}}-\text{\hspace{0.17em}}{E}_{f}\right)/kT}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}1}$ where$\text{ }{E}_{f}$ is the Fermi-level (half population point). If$\text{ }E\gg {E}_{f}$ , it can be approximated by Boltzmann distribution.


text:
# Isentropic Analysis

#### Objectives¶

First we need some grids of values to work with. We can do this by dowloading information from the latest run of the GFS available on Unidata's THREDDS data server. First we access the catalog for the half-degree GFS output, and look for the dataset called the "Best GFS Half Degree Forecast Time Series". This dataset combines multiple sets of model runs to yield a time series of output with the shortest forecast offset.

In [1]:
from siphon.catalog import TDSCatalog

cat = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/grib/'
'NCEP/GFS/Global_0p5deg/catalog.xml')
best = cat.datasets['Best GFS Half Degree Forecast Time Series']

Next, we set up access to request subsets of data from the model. This uses the NetCDF Subset Service (NCSS) to make requests from the GRIB collection and get results in netCDF format.

In [2]:
subset_access = best.subset()
query = subset_access.query()

Let's see what variables are available. Instead of just printing subset_access.variables, we can ask Python to only display variables that end with "isobaric", which is how the TDS denotes GRIB fields that are specified on isobaric levels.

In [3]:
sorted(v for v in subset_access.variables if v.endswith('isobaric'))

Out[3]:
['5-Wave_Geopotential_Height_isobaric',
'Absolute_vorticity_isobaric',
'Cloud_mixing_ratio_isobaric',
'Geopotential_height_isobaric',
'Graupel_snow_pellets_isobaric',
'Ice_water_mixing_ratio_isobaric',
'Icing_severity_isobaric',
'Ozone_Mixing_Ratio_isobaric',
'Rain_mixing_ratio_isobaric',
'Relative_humidity_isobaric',
'Snow_mixing_ratio_isobaric',
'Temperature_isobaric',
'Total_cloud_cover_isobaric',
'Vertical_velocity_geometric_isobaric',
'Vertical_velocity_pressure_isobaric',
'u-component_of_wind_isobaric',
'v-component_of_wind_isobaric']

Now we put together the "query"--the way we ask for data we want. We give ask for a wide box of data over the U.S. for the time step that's closest to now. We also request temperature, height, winds, and relative humidity. By asking for netCDF4 data, the result is compressed, so the download is smaller.

In [4]:
from datetime import datetime
query.time(datetime.utcnow())
query.variables('Temperature_isobaric', 'Geopotential_height_isobaric',
'u-component_of_wind_isobaric', 'v-component_of_wind_isobaric',
'Relative_humidity_isobaric')
query.lonlat_box(west=-130, east=-50, south=10, north=60)
query.accept('netcdf4')

Out[4]:
var=Relative_humidity_isobaric&var=Geopotential_height_isobaric&var=v-component_of_wind_isobaric&var=Temperature_isobaric&var=u-component_of_wind_isobaric&time=2020-01-17T15%3A13%3A24.543070&west=-130&east=-50&south=10&north=60&accept=netcdf4

Now all that's left is to actually make the request for data:

In [5]:
nc = subset_access.get_data(query)

Open the returned netCDF data using XArray:

In [6]:
from xarray.backends import NetCDF4DataStore
import xarray as xr
ds = xr.open_dataset(NetCDF4DataStore(nc))

### Isentropic Interpolation¶

Now let's take what we've downloaded, and use it to make an isentropic map. In this case, we're interpolating from one vertical coordinate, pressure, to another: potential temperature. MetPy has a function isentropic_interpolation that can do this for us. First, let's start with a few useful imports.

In [7]:
import metpy.calc as mpcalc
from metpy.units import units
import numpy as np

Let's parse out the metadata for the isobaric temperature and get the projection information. We also index with 0 to get the first, and only, time:

In [8]:
ds = ds.metpy.parse_cf()
temperature = ds['Temperature_isobaric'][0]

In [9]:
data_proj = temperature.metpy.cartopy_crs

Let's pull out the grids out into some shorter variable names.

In [10]:
lat = temperature.metpy.y
lon = temperature.metpy.x

# Need to adjust units on humidity because '%' causes problems
ds['Relative_humidity_isobaric'].attrs['units'] = 'percent'

rh = ds['Relative_humidity_isobaric'][0]
height = ds['Geopotential_height_isobaric'][0]
u = ds['u-component_of_wind_isobaric'][0]
v = ds['v-component_of_wind_isobaric'][0]

In [11]:
# Can have different vertical levels for wind and thermodynamic variables
# Find and select the common levels
press = temperature.metpy.vertical
common_levels = np.intersect1d(press, u.metpy.vertical)
temperature = temperature.metpy.sel(vertical=common_levels)
u = u.metpy.sel(vertical=common_levels)
v = v.metpy.sel(vertical=common_levels)

# Get common pressure levels as a data array
press = press.metpy.sel(vertical=common_levels)

Next, we perform the isentropic interpolation. At a minimum, this must be given one or more isentropic levels, the 3-D temperature field, and the pressure levels of the original field; it then returns the 3D array of pressure values (2D slices for each isentropic level). You can also pass addition fields which will be interpolated to these levels as well. Below, we interpolate the winds (and pressure) to the 320K isentropic level:

In [12]:
isen_level = np.array([320]) * units.kelvin
isen_press, isen_u, isen_v = mpcalc.isentropic_interpolation(isen_level, press,
temperature, u, v)

Let's plot the results and see what it looks like:

In [13]:
# Need to squeeze() out the size-1 dimension for the isentropic level
isen_press = isen_press.squeeze()
isen_u = isen_u.squeeze()
isen_v = isen_v.squeeze()

In [14]:
%matplotlib inline
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# Create a plot and basic map projection
fig = plt.figure(figsize=(14, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal(central_longitude=-100))

# Contour the pressure values for the isentropic level. We keep the handle
# for the contour so that we can have matplotlib label the contours
levels = np.arange(300, 1000, 25)
cntr = ax.contour(lon, lat, isen_press, transform=data_proj,
colors='black', levels=levels)
cntr.clabel(fmt='%d')

# Set up slices to subset the wind barbs--the slices below are the same as ::5
# We put these here so that it's easy to change and keep all of the ones below matched
# up.
lon_slice = slice(None, None, 5)
lat_slice = slice(None, None, 3)
ax.barbs(lon[lon_slice], lat[lat_slice],
isen_u[lat_slice, lon_slice].to('knots').magnitude,
isen_v[lat_slice, lon_slice].to('knots').magnitude,
transform=data_proj, zorder=2)

ax.set_extent((-120, -70, 25, 55), crs=data_proj)

/home/travis/miniconda/envs/unidata/lib/python3.7/site-packages/cartopy/io/__init__.py:260: DownloadWarning: Downloading: http://naciscdn.org/naturalearth/110m/physical/ne_110m_land.zip

EXERCISE: Let's add some moisture information to this plot. Feel free to choose a different isentropic level.
• Calculate the mixing ratio (using the appropriate function from mpcalc).
• Call isentropic_interpolation with mixing ratio. You should copy the one from above and add mixing ratio to the call so that it interpolates everything.
• contour (in green) or contourf your moisture information on the map alongside pressure.
You'll want to refer to the [MetPy API documentation](https://unidata.github.io/MetPy/latest/api/index.html) to see what calculation functions would help you.
In [15]:
# Needed to make numpy broadcasting work between 1D pressure and other 3D arrays
# Use .metpy.unit_array to get numpy array with units rather than xarray DataArray
pressure_for_calc = press.metpy.unit_array[:, None, None]

#
# YOUR CODE: Calculate mixing ratio using something from mpcalc
#

# Take the return and convert manually to units of 'dimenionless'
#mixing.ito('dimensionless')

#
# YOUR CODE: Interpolate all the data
#

# Squeeze the returned arrays
#isen_press = isen_press.squeeze()
#isen_mixing = isen_mixing.squeeze()
#isen_u = isen_u.squeeze()
#isen_v = isen_v.squeeze()

# Create Plot -- same as before
fig = plt.figure(figsize=(14, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal(central_longitude=-100))

levels = np.arange(300, 1000, 25)
cntr = ax.contour(lon, lat, isen_press, transform=data_proj,
colors='black', levels=levels)
cntr.clabel(fmt='%d')

lon_slice = slice(None, None, 8)
lat_slice = slice(None, None, 8)
ax.barbs(lon[lon_slice], lat[lat_slice],
isen_u[lat_slice, lon_slice].to('knots').magnitude,
isen_v[lat_slice, lon_slice].to('knots').magnitude,
transform=data_proj, zorder=2)

#
# YOUR CODE: Contour/Contourf the mixing ratio values
#

ax.set_extent((-120, -70, 25, 55), crs=data_proj)

SOLUTION
In [16]:
# %load solutions/mixing.py

# Cell content replaced by load magic replacement.
# Needed to make numpy broadcasting work between 1D pressure and other 3D arrays
# Use .metpy.unit_array to get numpy array with units rather than xarray DataArray
pressure_for_calc = press.metpy.unit_array[:, None, None]

# Calculate mixing ratio using something from mpcalc
mixing = mpcalc.mixing_ratio_from_relative_humidity(rh, temperature, pressure_for_calc)

# Take the return and convert manually to units of 'dimenionless'
mixing.ito('dimensionless')

# Interpolate all the data
isen_level = np.array([295]) * units.kelvin
ret = mpcalc.isentropic_interpolation(isen_level, press, temperature, mixing, u, v)
isen_press, isen_mixing, isen_u, isen_v = ret

# Squeeze the returned arrays
isen_press = isen_press.squeeze()
isen_mixing = isen_mixing.squeeze()
isen_u = isen_u.squeeze()
isen_v = isen_v.squeeze()

# Create Plot -- same as before
fig = plt.figure(figsize=(14, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal(central_longitude=-100))

levels = np.arange(300, 1000, 25)
cntr = ax.contour(lon, lat, isen_press, transform=data_proj,
colors='black', levels=levels)
cntr.clabel(fmt='%d')

lon_slice = slice(None, None, 8)
lat_slice = slice(None, None, 8)
ax.barbs(lon[lon_slice], lat[lat_slice],
isen_u[lat_slice, lon_slice].to('knots').magnitude,
isen_v[lat_slice, lon_slice].to('knots').magnitude,
transform=data_proj, zorder=2)

# Contourf the mixing ratio values
mixing_levels = [0.001, 0.002, 0.004, 0.006, 0.010, 0.012, 0.014, 0.016, 0.020]
ax.contourf(lon, lat, isen_mixing, transform=data_proj,
levels=mixing_levels, cmap='YlGn')

/home/travis/miniconda/envs/unidata/lib/python3.7/site-packages/metpy/interpolate/one_dimension.py:147: UserWarning: Interpolation point out of data bounds encountered

format: plain


text:
This post was also published on Future Processing - Technical Blog

I bet that all of us, developers, at least once started thinking about a domain model from designing data structures. This way of modelling applied to the complex domain almost always leads us to the solution that is not really effective in dealing with business cases. Moreover it is far from representing the real business domain - it represents a projection that is focused on our data store. The classical way of “producing” models for business domains treats events as a side effects of technical flows, not as a domain’s first class citizens. In the “event storming” domain modelling starts from events, which are assumed to be key points of business workflows. The rest of the domain is built to support events.

# TL;DR

1. Define events in the system
2. Define sources (commands, process managers) and actors
3. Group commands around aggregates
4. Implement…

Post-it, long piece of paper and people…

# Let’s play

Generally speaking “event storming” is about extracting domain from events which describe it. In most of the complex domains there are events that represent important changes in the state of our domain, simply begin with them. We’ll need the following things:

• Unlimited modelling space - sounds impossible, but in most cases very long piece of paper will be more than enough. White board is also OK, but it is hard to “serialize” it after the meeting.
• Post-its in various colours - one colour represents one class of objects.
• People who know the questions, and those who know the answers - in fact, in most cases, stakeholders will be assumed here, but the guy that is doing an actual work in the company may be useful too.
• A room without chairs - sitting = sleeping, it’s hard to go asleep during a meeting without chairs :). The key point here is to engage everyone into the work.

I’m assuming that there will be at least one person in the room who knows what event storming is, and this person will be the moderator. It’s a good practice for moderator to place the first event on the modelling space. Simply put one post-it with the event name, a chosen colour will represent events. The modelling space represented by a long piece of paper should have a direction which will represent the time, so we can arrange events chronologically.

But what event actually is? There is a “simple” definition:

Domain event is something that happened that domain experts care about.

Let business experts put events that are in their region of interest on the paper, and let them argue about how the business works (but don’t let them kill each other)! They should work out the clean description of the business flow. The value here is to clarify complex workflows by finding common points. If they are not approaching common point, they must have found the misconception in their business - it is an opportunity to improve. Or maybe they are simply talking about the same thing from completely different perspectives. If one expert is responsible for invoicing and another one for delivery - they can talk about the same thing, but in two separated subdomains, so that thing looks completely different for each of them.

When all events are in place we can go to the next phase.

# Reverse engineering

Starting from the events we want to get to the point in which we have all important structures clearly visible on the modelling space. When we will think about events that are happening all the time around us everything has a more or less obvious cause. Similarly, in the domain, nothing happens without cause. We need to find all causes for events.

Event can be triggered by a command, by time or even by another event - so basically every event has the source. Obviously, event may also have more than one source i.e.: AssessmentCompletedEvent may be triggered by time when Assessment’s end date expires or by user which manually completes it in the application. Our task is to discover sources of events. In this phase we need to put on the modelling space all commands, process managers and actors that are involved in generating events.

Command is a well-known pattern and even with a bunch of rules of thumb that comes from DDD there is no reason to define it again.

Process Manager is a simple creature that executes commands when event is raised, so if we need to execute command in reaction to event we need the Process Manager to do this work for us.

Actor is an external source of activity - users, external systems and even time (here are some controversies) may be treated as an actor.

The moment we have defined all sources for our events, we can group commands around aggregates changed by them. Basically that should be iterative process where we are deciding which command modifies which aggregate. From my perspective this might be the most difficult part of the modelling as designing good aggregate is always hard & tricky!

# Implement!

The process described above should generate something :) that theoretically should be easy to express as a code. We have domain aggregates, commands that operates on them and events that represents state changes in aggregates. Every artefact that we generate before can be easily represented by building block from DDD and match to CQRS/ES architecture. There is nothing better to do than implement this model.

# Is Event Storming good for my project?

From my point of view any domain can be modelled in this way.


text:
# Derivs heeelp

Printable View

• Dec 14th 2008, 08:58 AM
bballa273
Derivs heeelp
f(x)=[(x+2)/(x+1)](2x-7)

I keep getting [-2x^2+7x-3] / (x+1)
but that is not one of the answer choices.

Also, (x-2)/sqrt(x^2-1)

Find an equation for the tan line to the graph f at the point P(pi/4, f(pi/4)) when f(x)=9tan x

Let f(x) = (4x-1)^6(3x+3)^2 Find the x-intercept of the tan line to the graph of f at the pt where the graph of f crosses the y-axis
• Dec 14th 2008, 07:28 PM
Rapha
Hi

Quote:

Originally Posted by bballa273
f(x)=[(x+2)/(x+1)](2x-7)

I keep getting [-2x^2+7x-3] / (x+1)
but that is not one of the answer choices.

this is indeed really hard

First: Quotient rule

$[(x+2)/(x+1)]' = \frac{1*(x+1)-1*(x+2)}{(x+1)^2}=:v'$

----
(2x-7)' = 2 =: u'

Product rule (u*v)' = u'*v + v'*u

and so on...

But it is much easier, if you multiply (x+2)(2x-7) first. Then you only must use the quotient rule

Quote:

Originally Posted by bballa273
Also, (x-2)/sqrt(x^2-1)

quotient rule: $(\frac{u}{v})' = \frac{u'*v-v'*u}{v^2}$

here: u= x-2 => u'=1

$v = \sqrt{x^2-1} = (x^2-1)^{0.5}$

$v' = 0.5 * 2x * (x^2-1)^{-0.5}$

...

Quote:

Originally Posted by bballa273
Find an equation for the tan line to the graph f at the point P(pi/4, f(pi/4)) when f(x)=9tan x

Let f(x) = (4x-1)^6(3x+3)^2 Find the x-intercept of the tan line to the graph of f at the pt where the graph of f crosses the y-axis

I dont know. Maybe someone else does

Rapha
• Dec 14th 2008, 07:33 PM
mr fantastic
Quote:

Originally Posted by bballa273
[snip]
Find an equation for the tan line to the graph f at the point P(pi/4, f(pi/4)) when f(x)=9tan x

[snip]

You know two things about the line:

1. Point on line is $\left(\frac{\pi}{4}, \, 9 \right)$ (why?)

2. Gradient of line is $9 \sec^2 \left( \frac{\pi}{4} \right) = \, ....$ (why?)

So you should be able to get the equation of the line.
• Dec 14th 2008, 07:36 PM
mr fantastic
Quote:

Originally Posted by bballa273
[snip]
Let f(x) = (4x-1)^6(3x+3)^2 Find the x-intercept of the tan line to the graph of f at the pt where the graph of f crosses the y-axis

1. Find the y-intercept of y = f(x). This gives you a known point on the tangent line.

2. Evaluate the derivative of f(x) at the y-intercept. This gives you the gradient of the tangent line.

3. Use 1. and 2. to get the equation of the tangent line.

4. Use the equation in 3. to get the x-intercept of the tangent.


text:
Share on:

## Lewis Carroll's game of Doublets

Such a long time since my last post! Well, that's the working life for you.

Anyway, recently I was reading about Lewis Carroll - always one of my favourite people - and was reminded of his word game "Doublets" in which one word is turned into another by changing one letter at a time, each new word being English.

You can read Carroll's original description here. Note his last sentence:

"It is, perhaps, needless to state that it is de rigueur that the links should be English words, such as might be used in good society."

Carroll, it seemed, frowned on slang words, or "low" words - very much in keeping with his personality and with his social and professional positions. One of his examples was

which has an answer: WHEAT CHEAT CHEAP CHEEP CREEP CREED BREED BREAD.

Clearly we would want the length of the chain of words to be as short as possible; with a lower bound being the Hamming distance between the words: the number of places in which letters are different. This distance is 3 for WHEAT and BREAD and so the minimum changes will be 3. But in practice chains are longer, as the one above. English simply doesn't contain all possible words of 5 letters, and so we can't have, for example:

This form of word puzzle, so simple and addictive, has been resurrected many times, often under such names as "word ladders", or "laddergrams".

## Obtaining word lists

Every computer system will have a spell-check list on it; on a Linux system these are usually found under /usr/share/dict . My system, running Arch Linux has these lists:

 1$wc -l /usr/share/dict/* 2 123115 /usr/share/dict/american-english 3 127466 /usr/share/dict/british-english 4 189057 /usr/share/dict/catalan 5 54763 /usr/share/dict/cracklib-small 6 88328 /usr/share/dict/finnish 7 221377 /usr/share/dict/french 8 304736 /usr/share/dict/german 9 92034 /usr/share/dict/italian 10 76258 /usr/share/dict/ogerman 11 56329 /usr/share/dict/spanish 12 123115 /usr/share/dict/usa  These all come from the package words, which is a collection of spell-check dictionaries. Although the sizes of the English word lists may seem impressive, there are bigger lists available. One of the best is at SCOWL (Spell Checker Oriented Word Lists). You can download the compressed SCOWL file, and when uncompressed you'll find it contains a directory called final. In this directory, the largest files are those of the sort english-words.*, and here's how big they are:  1$ wc -l english-words.*
2   4373 english-words.10
3   7951 english-words.20
4  36103 english-words.35
5   6391 english-words.40
6  23796 english-words.50
7   6233 english-words.55
8  13438 english-words.60
9  33270 english-words.70
10 139209 english-words.80
11 219489 english-words.95
12 490253 total

These lists contain increasingly more abstruse and unusual words. Thus english-words.10 contains words that most adept speakers of English would know:

 1$shuf -n 10 english-words.10 2errors 3hints 4green 5connections 6still 7mouth 8category's 9pi 10won's 11varied  At the other end, english-words.95 consists of unusual, obscure words unlikely to be in a common vocabulary; many of them are seldom used, have very specific meanings, or are technical:  1$ shuf -n 10 english-words.95
2deutschemark's
3disingenious
4retanner
6shlimazl
7unpontifical
8nonrequirement
9peccancy's
10photozinco
11nonuniting

This is the list which contains some splendid biological terms: "bdelloid", which is a class of microscopic water animals called rotifers; "ctenizidae", a small class of spiders (but the list does not contain "ctenidae" the much larger class of wandering spiders, including the infamous Brazilian wandering spider); "cnidocyst" which forms the stinging mechanism in the cells of the tentacles of jellyfish.

Putting the lists 10, 25, 30 together makes a "small" list; adding in 40 and 50 gives a "medium" list; next include 55, 65, 70 for "large"; "80" for "huge", and 95 for "insane". The author of SCOWL claims that 60 is the right level for spell checking, while: "The 95 contains just about every English word in existence and then some. Many of the words at the 95 level will probably not be considered valid English words by most people." (This is not true. I've discovered some words not in any list. One such is "buckleys", as in "He hansn't got buckleys chance", sometimes spelled with a capital B, and meaning "He hasn't got a chance"; that is, no chance at all. "Buckley" is only in the list of proper names, but given this usage it should be at least in the Australian words list. Which it isn't.)

You will also see from above that some words contain apostrophes: these words will need to be weeded out. Here's how to make a list and clean it up:

1$cat english-words.10 english-words.25 english-words.30 > english-words-small.txt 2$ grep -v "'" english-words-small.txt > small.wds

This approach will yield five lists with the following numbers of words:

1\$ wc -l *.wds
2 234563 huge.wds
3 414365 insane.wds
4 107729 large.wds
5  59577 medium.wds
6  38013 small.wds
7 854247 total

These can be read into Julia. These lists are unsorted, but that won't be an issue for our use of them. But you can certainly include a sort along the way.

Another list is available at https://github.com/dwyl/english-words which claims to have "over 466k English words". However, this list is not as carefully curated as is SCOWL.

Finally note that our lists are not disjoint, as are the original lists. Each list includes its predecessor, so that insane.wds contains all of the words in all of the lists.

## Using graph theory

The computation of word ladders can easily be managed using the tools of graph theory, with vertices being the words, and two vertices being adjacent if their Hamming distance is 1. Then finding a word ladder is easily done by a shortest path.

There is a problem though, as Donald Knuth discovered when he launched the first computerization of this puzzle, of which an explanation is available in his 1994 book The Stanford GraphBase. This page, you'll notice, contains "the 5757 five-letter words of English". However, deciding what is and what isn't an English word can be tricky: American versus English spellings, dialect words, newly created words and so on. I touched on this in an earlier post.

Knuth also found that there were 671 words which were connected to no others; he called these "aloof words", with ALOOF, of course, being one of them.

## Using Julia

Although Python has its mature and powerful NetworkX package for graph theory and network analysis, Python is too slow for this application: we are looking at very large graphs of many thousands of vertices, and computing the edges is a non-trivial task. So our choice of language is Julia. Julia's graph theory packages are in a bit of state of flux, an old package Graphs.jl is unmaintained, as is LightGraphs.jl. However, this latter package is receiving a new lease of life with the unfortunately confusing name of Graphs.jl, and which is designed to be "functionally equivalent" to LightGraphs.jl.

This is the package I'll be using.

## Setting up the graph

We'll use the medium.wds dictionary since it's relatively small, and we'll look at five letter words. Using six-letter words or the larger list will then be a simple matter of changing a few parameters.

We start by open the list in Julia:

1Julia> f = open("medium.wds", "r")
3Julia> length(words)
4
5       59577

Now we can easily extract the five letter words, and set up the graph, first of all loading the Graphs package. We also need the StringDistances package to find the Hamming distance.

1Julia> using Graphs, StringDistances
2Julia> words5 = filter(x -> length(x)==5, words)
3Julia> w5 = length(words5)
4
5       4388
6
7Julia> G5 = Graph()

Now the edges:

1Julia> for i in 1:w5
2	   for j in i+1:w5
3	       if Hamming()(words5[i],words5[j]) == 1
5	       end
6	   end
7       end
8Julia>

Note that there is a Julia package MetaGraph.jl which allows you to add labels to edges. However, it's just as easy to use the vertex numbers as indices into the list of words.

We can't use the graph G5 directly, as it is not a connected graph (remember Donald Knuth's "aloof" words?) We'll do two things: find the aloof words, and choose from G5 the largest connected component. First the aloof words:

1Julia> aloofs = map(x->words5[x],findall(iszero, degree(G5)))
2Julia> show(aloofs)

I won't include this list, as it's too long - it contains 616 words. But if you do so, you'll see some surprises here: who'd have though that such innocuous words as "opera", "salad" or "wagon" were aloof? But they most certainly are, at least within this set of words.

And now for the connected component:

1Julia> CC = connected_components(G5)
2Julia> CL = map(length,CC)  # size of each component
3Julia? mx, idx = findmax(CL)
4Julia> C5,vmap = induced_subgraph(G5, CC[idx])

You will find that the maximum connected component C5 has 3315 vertices, and the value of idx above is 3. Here vmap is a list of length 3315 which is the set of indices into the original list. Or, if the original list of vertices consisted of the numbers 1,2, up to 4388, then vmap is a sublist of length 3315. And we can consider all those numbers as indices into the words5 list.

Now we can write a little function to produce word ladders; here using the Julia function a_star to find a shortest path:

1Julia> subwords = words5[vmap]
3	   i = findfirst(x -> words5[x] == w1, vmap)
4	   j = findfirst(x -> words5[x] == w2, vmap)
5	   P = a_star(C5,i,j)
6	   verts = append!([i],map(dst,P))
7	   subwords[verts]
8       end

And of course try it out:

1Julia> show(ladder("wheat","bread"))
2
3       ["wheat", "cheat", "cleat", "bleat", "bleak", "break", "bread"]

Note that a word ladder can only exist when both words have indices in the chosen largest connected component. For example:

1Julia> show(ladder("roses","scent"))
2
3       ["roses", "ruses", "rusts", "rests", "tests", "teats", "seats", "scats",
4	"scans", "scant", "scent"]

but

1Julia> show(ladder("chair","table"))

will produce an error, as neither of those words have indices in the largest connected component. In fact "chair" sits in a small connected component of 3 words, and "table" is in another connected component of 5 words.

We have seen that the length of a word ladder will often exceed the Hamming distance between the start and ending words. But what is the maximum length of such a ladder?

Here the choice of function is the eccentricities of a graph: for each vertex, find the shortest path to every other vertex. The length of the longest such path is the eccentricity of that vertex.

1Julia> eccs = eccentricity(C5)

This command will take a noticeable amount of time - only a few seconds, it's true, but it is far from instantaneous. The intense computation needed here is one of the reasons that I prefer Julia to Python for this experiment.

Now we can use the eccentricities to find the longest path. Since this is an undirected graph, there must be at least two vertices with equal largest eccentricities.

1Julia> mx = maximum(eccs)
2Julia> inds = findall(x -> x == mx, eccs)
4Julia> print(join(ld,",  "),"\nwhich has length ",length(ld))
5
6aloud,  cloud,  clout,  flout,  float,  bloat,  bleat,  cleat,  cleft,  clefs,  clews,
7slews,  sleds,  seeds,  feeds,  feels,  fuels,  furls,  curls,  curly,  curvy,  curve,
8carve,  calve,  valve,  value,  vague,  vogue,  rogue"
9which has length  29

Who'd have guessed?

## Results from other lists

With the 12478 five-letter words from the "large" list, there are 748 aloof words, and the longest ladder is

rayon, racon, bacon, baron, boron, boson, bosom, besom, besot, beset, beret, buret, curet, cures, curds, surds, suras, auras, arras, arias, arils, anils, anile, anole, anode, abode, abide, amide, amine, amino, amigo

which has length 31. We might in general expect a shorter "longest ladder" than from using smaller list: the "large" list has far more words, hence greater connectivity, which would lead in many cases to shorter paths.

The "huge" and "insane" lists have longest ladders of length 28 and respectively:

essay, assay, asway, alway, allay, alley, agley, aglet, ablet, abled, ailed, aired, sired, sered, seres, serrs, sears, scars, scary, snary, unary, unarm, inarm, inerm, inert, inept, inapt, unapt

and

entry, entsy, antsy, artsy, artly, aptly, apply, apple, ample, amole, amoke, smoke, smote, smite, suite, quite, quote, quott, quoit, qubit, oubit, orbit, orbic, urbic, ureic, ureid, ursid

For six-letter words, here are longest ladders in each of the lists:

Small, length 32:

steady, steamy, steams, steals, steels, steers, sheers, cheers, cheeks, checks, chicks, clicks, slicks, slices, spices, spites, smites, smiles, smiled, sailed, tailed, tabled, tabbed, dabbed, dubbed, rubbed, rubber, rubier, rubies, rabies, babies, babied

Medium, length 44:

trusty, trusts, crusts, crests, chests, cheats, cleats, bleats, bloats, floats, flouts, flours, floors, floods, bloods, broods, brooks, crooks, crocks, cracks, cranks, cranes, crates, crated, coated, boated, bolted, belted, belied, belies, bevies, levies, levees, levers, lovers, hovers, hovels, hotels, motels, models, modals, morals, morays, forays

Large, length 43:

uneasy, unease, urease, crease, creese, cheese, cheesy, cheeky, cheeks, creeks, breeks, breeds, breads, broads, broods, bloods, bloops, sloops, stoops, strops, strips, stripe, strive, shrive, shrine, serine, ferine, feline, reline, repine, rapine, raping, raring, haring, hiring, siring, spring, sprint, splint, spline, saline, salina, saliva

Huge, length 63:

aneath, sneath, smeath, smeeth, smeech, sleech, fleech, flench, clench, clunch, clutch, crutch, crotch, crouch, grouch, grough, trough, though, shough, slough, slouch, smouch, smooch, smooth, smoots, smolts, smalts, spalts, spaits, spains, spaing, spring, sprint, splint, spline, upline, uplink, unlink, unkink, unkind, unkend, unkent, unsent, unseat, unseal, unseel, unseen, unsewn, unsews, unmews, enmews, emmews, emmers, embers, embars, embark, imbark, impark, impart, import, impost, impose, impone

Insane, length 49:

ambach, ambash, ambush, embush, embusk, embulk, embull, emball, empall, emparl, empark, embark, embars, embers, emmers, emmews, enmews, endews, enders, eiders, ciders, coders, cooers, cooees, cooeed, coomed, cromed, cromes, crimes, crises, crisis, crisic, critic, iritic, iridic, imidic, amidic, amidin, amidon, amydon, amylon, amylin, amyrin, amarin, asarin, asaron, usaron, uzaron, uzarin

The length of time taken for the computations increases with size of the lists, and up to a certain extent, the length of the words. On my somewhat mature laptop (Lenovo X1 Carbon 3rd Generation), computing the eccentricities for six-letter words and the "insane" list took over 6 minutes.

## A previous experiment with Mathematica

You can see this (with some nice graph diagrams) at https://blog.wolfram.com/2012/01/11/the-longest-word-ladder-puzzle-ever/ from about 10 years ago. This experiment used Mathematica's English language list, about which I can find nothing other than it exists. However, in the comments, somebody shows that Mathematica 8 had 92518 words in its English dictionary. And this number is dwarfed by Finnish, with 728498 words. But that might just as much be an artifact of Finnish lexicography. Deciding whether or nor a word is English is very hard indeed, and any lexicographer will decide on whether two words need separate definitions, or whether one can be defined within the other, so to speak - that is, under the same headword. MOUSE, MOUSY, MOUSELIKE - do these require three definitions, or two, or just one?


text:
Please use this identifier to cite or link to this item: http://hdl.handle.net/10481/28416

 Title: AMIGA: very low environment galaxies in the local universe Authors: Verley, S.Combes, F.Verdes-Montenegro, L.Leon, S.Odewahn, S. C.Bergond, G.Espada, D.García, E.Lisenfeld, UteSabater, J.Sulentic, J. Issue Date: 2006 Abstract: The evolutionary history of galaxies is thought to be strongly conditioned by the environment. In order to quantify and set limits on the role of nurture one must identify and study an isolated sample of galaxies. But it is not enough to identify a small number of the "most isolated" galaxies. We begin with 950 galaxies from the Catalogue of Isolated Galaxies (Karachentseva 1973) and reevaluate isolation using an automated star-galaxy classification procedure on large digitised POSS-I fields. We define, compare and discuss various criteria to quantify the degree of isolation for these galaxies: Karachentseva's revised criterion, local surface density computations and an estimation of the external tidal force affecting each isolated galaxy. Comparison of multi-wavelength ISM properties, in particular the H$\alpha$ emission line, will allow us to separate the influence of the environment from the one due to the initial conditions at formation. Publisher: EDP Sciences Keywords: AMIGA (Analysis of the interstellar Medium of Isolated Galaxies)GalaxiesISM URI: http://hdl.handle.net/10481/28416 Rights : Creative Commons Attribution-NonCommercial-NoDerivs 3.0 License Citation: Verley, S.; et al. AMIGA: very low environment galaxies in the local universe. In: F. Casoli, et al. (eds.), SF2A 2005: Seamine de l'Astrophysique Francaise. Strasbourg (France), 27 june- 1 july 2005. pp. 237-238. [http://hdl.handle.net/10481/28416] Appears in Collections: DFTC - Comunicaciones congresos, conferencias,...

Files in This Item:

File Description SizeFormat


text:
# ABC is a right angled triangle in which ∠A = 90º and AB = AC. Find ∠B and ∠C.

Solution:

It is given that

$\mathrm{AB}=\mathrm{AC}$

$\Rightarrow \angle C=\angle B$ (Angles opposite to equal sides are also equal)

$\ln \Delta \mathrm{ABC}$

$\angle A+\angle B+\angle C=180^{\circ}$ (Angle sum property of a triangle)

$\Rightarrow 90^{\circ}+\angle B+\angle C=180^{\circ}$

$\Rightarrow 90^{\circ}+\angle B+\angle B=180^{\circ}$

$\Rightarrow 2 \angle B=90^{\circ}$

$\Rightarrow \angle B=45^{\circ}$

$\therefore \angle B=\angle C=45^{\circ}$


text:
Hi people

I have to work out the area and diagonals of a quadrilateral,
the side lengths are AB = 19.4, BC = 27.6, CD= 16.9, DA= 213 and
angle B = 65 degrees
can anybody help with the calculations/ answer as im completely stumped,
any advise would be greatly appreciated

thanks John

2. Originally Posted by ratherbsurfing
I have to work out the area and diagonals of a quadrilateral,
the side lengths are AB = 19.4, BC = 27.6, CD= 16.9, DA= 213 and
angle B = 65 degrees
Start by drawing the diagonal AC. You can then use the cosine rule to find the length AC and then the angles in the triangles ABC and ACD. The area of the quadrilateral is the sum of the areas of the two triangles (which you get from the formula " $\tfrac12ab\sin C$"). Finally, draw the other diagonal BD and get its length by using the cosine rule in the triangle BCD (because by then you will know the angle at C: it is the sum of the angles ACB and ACD).

3. Opalg

thats brilliant thanks Ill give it a try, think ill still struggle but here goes

thanks John

4. Please note that this is not an advanced geometry question thus it belongs in the pre-university forums, not the university forums.


text:
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

It is Sunday, it's raining and I have a few hours to spend before I am invited for lunch at my parents place. Hence, I thought I'd use the time to produce another post. It has been a while since the last post as I have been in Africa for about two months for yet another stint of fieldwork on the great Kilimanjaro mountain…

This time I want to introduce the strip() function which is part of the METVURST suite for plotting meteorological/climatological data (available from https://github.com/tim-salabim/metvurst)

This function is intended to facilitate two things:

1. to enable plotting of large (climatological) data sets at hourly resolution (the data need to be hourly observations) in a reasonably well defined space (meaning that you won't need pages and pages of paper to print the results).
2. to display the data in a way that interpretation is possible from hourly to decadal time scales.

The function is much like the calendarHeat() function from Revolution Analytics (see here) though, as mentioned above, it produces plots of hourly data. In fact, calendarHeat() would be a good alternative to use when you intend to plot daily values.
strip() is implemented in lattice (needs latticeExtra to be more precise) and essentially plots a levelplot of hour of day (y-axis) vs. day of year (x-axis).

In detail the function takes the following parameters:

• x (numeric): Object to be plotted (e.g. temperature).
• date (character or POSIXct): Date(time) of the observations. Format must be 'YYYY-MM-DD hh:mm(:ss)'
• fun: The function to be used for aggregation to hourly observations (if original is of higher frequency). Default is mean.
• cond (factor): Conditioning variable e.g. years.
• arrange (character): One of “wide” or “long”. Defaults to “long” which provides a layout that is easier to interpret, however, “wide” is better for use in presentation slides.
• colour (character): a vector of color names. Defaults to rev(brewer.pal(11, “Spectral”))
• …: Further arguments to be passed to levelplot (see ?lattice::levelplot for options). This can be quite handy to set the legend title using main = "YOUR TEXT HERE" (the legend is plotted on top of the actual plot).

Here's an example using hourly data from Fiji (Station is at Lautoka on the west coast). The data can be freely accessed through the Bureau of Meteorology in Australia (BOM).

## LOAD RColorBrewer FOR COLOR CREATION
library(RColorBrewer)

url <- "http://www.bom.gov.au/ntc/IDO70004/IDO70004_"

yr <- 1993:2012

## READ DATA FOR ALL YEARS FROM URL INTO LIST
fijilst <- lapply(seq(yr), function(i) {
read.csv(paste(url, yr[i], ".csv", sep = ""), na.strings = c(-9999, 999))
})

## TURN LIST INTO COMPLETE DATAFRAME AND CONVERT NA STRINGS TO NAs
fiji <- do.call("rbind", fijilst)
fiji[fiji == -9999.00] <- NA
fiji[fiji == -9999.0] <- NA
fiji[fiji == 999.0] <- NA

## CREATE POSIX DATETIME AND CONVERT UTC TO LOCAL FIJI TIME
dts <- as.POSIXct(strptime(fiji$Date...UTC.Time, format = "%d-%b-%Y %H:%M")) + 12 * 60 * 60 ## CREATE CONDITIONING VARIABLE (IN THIS CASE YEAR) year <- substr(as.character(dts), 1, 4)  Now, let's plot the temperatures ## SOURCE FUNCTION (ALSO AVAILABLE AT https://github.com/tim-salabim/metvurst) source("http://www.staff.uni-marburg.de/~appelhat/r_stuff/strip.R") ## PLOT STRIP FOR TEMPERATURE strip(x = fiji$Air.Temperature,
date = dts,
cond = year,
arrange = "long",
main = "Temperature")

##
##  Module   :  strip
##  Author   :  Tim Appelhans <[email protected]>, Thomas Nauss
##  Version  :  2012-01-06

We can clearly see both the diurnal and the seasonal signal.
Looking at other parameters such as wind direction and speed should also give an interesting insight into the seasonal and diurnal dynamics of the location.

I leave it up to you to explore this further…

I am off to a nice sunday roast now! yumyum

sessionInfo()

## R version 2.15.3 (2013-03-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
##
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
##  [7] LC_PAPER=C                 LC_NAME=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] grid      parallel  stats     graphics  grDevices utils     datasets
## [8] methods   base
##
## other attached packages:
##  [1] reshape_0.8.4       plyr_1.8            latticeExtra_0.6-19
##  [4] lattice_0.20-13     RColorBrewer_1.0-5  RWordPress_0.2-3
##  [7] rgdal_0.8-5         raster_2.0-41       sp_1.0-5
## [10] knitr_1.1
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.3   evaluate_0.4.3 formatR_0.7    markdown_0.5.4
## [5] RCurl_1.95-3   stringr_0.6.2  tools_2.15.3   XML_3.95-0
## [9] XMLRPC_0.2-5


text:
# Berry Curvature Engineering by Gating Two-Dimensional Antiferromagnets

@article{Du2019BerryCE,
title={Berry Curvature Engineering by Gating Two-Dimensional Antiferromagnets},
author={Shiqiao Du and Peizhe Tang and Jiaheng Li and Z. Lin and Y. Xu and W. Duan and A. Rubio},
journal={arXiv: Materials Science},
year={2019}
}
• Shiqiao Du, +4 authors A. Rubio
• Published 2019
• Physics, Materials Science
• arXiv: Materials Science
Recent advances in tuning electronic, magnetic, and topological properties of two-dimensional (2D) magnets have opened a new frontier in the study of quantum physics and promised exciting possibilities for future quantum technologies. In this study, we find that the dual-gate technology can well tune the electronic and topological properties of antiferromagnetic (AFM) even septuple-layer (SL) MnBi$_2$Te$_4$ thin films. Under an out-of-plane electric field that breaks $\mathcal{PT}$ symmetry… Expand
6 Citations


text:
For in GOD we live, and move, and have our being. - Acts 17:28

The Joy of a Teacher is the Success of his Students. - Samuel Dominic Chukwuemeka

# Solved Examples - Introductory Statistics

For ACT Students
The ACT is a timed exam...$60$ questions for $60$ minutes
This implies that you have to solve each question in one minute.
Some questions will typically take less than a minute a solve.
Some questions will typically take more than a minute to solve.
The goal is to maximize your time. You use the time saved on those questions you solved in less than a minute, to solve the questions that will take more than a minute.
So, you should try to solve each question correctly and timely.
So, it is not just solving a question correctly, but solving it correctly on time.
Please ensure you attempt all ACT questions.
There is no "negative" penalty for any wrong answer.

For JAMB and CMAT Students
Calculators are not allowed. So, the questions are solved in a way that does not require a calculator.

Attempt all questions.
Show all work.

(1.) The table below shows the voltage measurements from a home in the City of Truth or Consequences, New Mexico for $25$ days.

Day Voltage (volts)
$1$ $121.3$
$2$ $121.1$
$3$ $121.5$
$4$ $121.7$
$5$ $122.0$
$6$ $121.2$
$7$ $121.2$
$8$ $121.3$
$9$ $121.2$
$10$ $121.7$
$11$ $121.6$
$12$ $121.9$
$13$ $121.5$
$14$ $121.4$
$15$ $121.8$
$16$ $121.8$
$17$ $121.2$
$18$ $121.4$
$19$ $121.1$
$20$ $121.9$
$21$ $121.7$
$22$ $121.9$
$23$ $121.9$
$24$ $121.3$
$25$ $121.6$

$(a.)$ Draw a frequency distribution table for the data. Your table should have $5$ classes.
$(b.)$ Compute the statistical properties of the classes.

Maximum voltage value = $122.0$
Minimum voltage value = $121.1$
Range = $122.0 - 121.1 = 0.9$
Number of classes = $5$
Class Width = $\dfrac{0.9}{5} = 0.18 \approx 0.2$

Based on the note, Data Organization; If the class intervals (the $LCI$ and the $UCI$) are decimals rounded to one decimal place, then the difference between the lower class interval of a class and the upper class interval of the previous class is $0.1$

$LCI\:\:of\:\:2nd\:\:Class - UCI\:\:of\:\:1st\:\:Class = 0.1 \\[3ex] 1st\:\:Class = 121.1 - 121.2 \\[3ex] 2nd\:\:Class = 121.3 - 121.4 \\[3ex] UCB\:\:of\:\:1st\:\:Class = \dfrac{UCI\:\:of\:\:1st\:\:Class + LCI\:\:of\:\:2nd\:\:Class}{2} = \dfrac{121.2 + 121.3}{2} = \dfrac{242.5}{2} = 121.25 \\[5ex] LCB\:\:of\:\:1st\:\:Class = \dfrac{LCI\:\:of\:\:1st\:\:Class + UCI\:\:of\:\:previous\:\:Class}{2} \\[5ex]$ We do not have a previous class.
But, assuming we did; the $UCI$ of that class would be $121.2 - 0.2 = 121.0$

$\therefore LCB\:\:of\:\:1st\:\:Class = \dfrac{121.1 + 121.0}{2} = \dfrac{242.1}{2} = 121.05 \\[5ex] Class\:\:Boundary\:\:of\:\:1st\:\:Class = 121.05 - 121.25 \\[3ex]$ The Frequency Distribution Table is constructed as shown:

Voltage Intervals Tally Frequency, $F$ Class Midpoints Class Boundaries Relative Frequency, $RF$ Cumulative Frequency, $CF$
$121.1 - 121.2$ IIII I $6$ $\dfrac{121.1 + 121.2}{2} = 121.15$ $121.05 - 121.25$ $\dfrac{6}{25} = 0.24 = 24\%$ $6$
$121.3 - 121.4$ IIII $5$ $\dfrac{121.3 + 121.4}{2} = 121.35$ $121.25 - 121.45$ $\dfrac{5}{25} = \dfrac{1}{5} = 0.2 = 20\%$ $6 + 5 = 11$
$121.5 - 121.6$ IIII $4$ $\dfrac{121.5 + 121.6}{2} = 121.55$ $121.45 - 121.65$ $\dfrac{4}{25} = 0.16 = 16\%$ $11 + 4 = 15$
$121.7 - 121.8$ IIII $5$ $\dfrac{121.7 + 121.8}{2} = 121.75$ $121.65 - 121.85$ $\dfrac{5}{25} = \dfrac{1}{5} = 0.2 = 20\%$ $15 + 5 = 20$
$121.9 - 122.0$ IIII $5$ $\dfrac{121.9 + 122.0}{2} = 121.95$ $121.85 - 122.05$ $\dfrac{5}{25} = \dfrac{1}{5} = 0.2 = 20\%$ $20 + 5 = 25$
$\Sigma F = 25$ $\Sigma RF = 1 = 100\%$

(2.) The ages of the $2016$ United States presidential candidates from $4$ political parties are:

 $70$ $64$ $45$ $45$ $65$ $63$ $66$ $54$ $62$ $58$ $53$ $61$ $71$ $61$ $45$ $48$ $66$ $53$ $63$ $70$ $55$ $68$ $75$ $65$ $63$

Use $5$ classes to form a grouped data.

(3.) $120$ Nursing majors took a standardized test.
The scores are summarized in the Frequency Table as shown:

Scores Frequency Scores Cumulative Frequency
$160 - 179$ $17$ $Less\:\:than\:\:180$ $17$
$180 - 199$ $20$ $Less\:\:than\:\:200$ $37$
$200 - 219$ $19$ $Less\:\:than\:\:220$ $56$
$220 - 239$ $x$ $Less\:\:than\:\:240$ $70$
$240 - 259$ $17$ $Less\:\:than\:\:260$ $87$
$260 - 279$ $33$ $Less\:\:than\:\:280$ $y$

Calculate the values of $x$ and $y$

$56 + x = 70 \\[3ex] x = 70 - 56 \\[3ex] x = 14 \\[3ex] 87 + 33 = y \\[3ex] 120 = y \\[3ex] y = 120$

(4.) CSEC The cumulative frequency distribution of the volume of petrol needed to fill the tanks of $150$ different vehicles is shown below.

Volume (litres) Cumulative Frequency
$11 - 20$ $24$
$21 - 30$ $59$
$31 - 40$ $101$
$41 - 50$ $129$
$51 - 60$ $150$

(a.) For the class $21 - 30$, determine the
(i) lower class boundary
(ii) class width

(b.) How many vehicles were recorded in the class $31 - 40$?

(c.) A vehicle is chosen at random from the $150$ vehicles. What is the probability that the volume of petrol needed to fill its tank is more than $50.5$ litres? Leave your answer as a fraction.

(d.) Byron estimates the median amount of petrol to be $43.5$ liters. Explain why Byron's estimate is INCORRECT.

(e.) On the partially labelled grid below, construct a histogram to represent the distribution of the volume of petrol needed to fill the tanks of the $150$ vehicles.

(a.) For the class $21 - 30$,

$(i)\:\: LCB = \dfrac{21 + 20}{2} = \dfrac{41}{2} = 20.5 \\[5ex] (ii)\:\: CW = 21 - 11 = 10 \\[3ex]$ Let us re-draw the table to show the frequencies
Volume (litres) Frequency Cumulative Frequency, $CF$
$11 - 20$ $24$ $24$
$21 - 30$ $59 - 24 = 35$ $59$
$31 - 40$ $101 - 59 = 42$ $101$
$41 - 50$ $129 - 101 = 28$ $129$
$51 - 60$ $150 - 129 = 21$ $150$

(b.) The number of vehicles recorded in the $31 - 30$ class is the frequency of that class.
Therefore, $42$ vehicles were recorded in the class $31 - 40$

(c.) Prerequisite to answer (c.) is: Probability
The number of vehicles whose tank is more than $50.5$ litres is the frequency of the $51 - 60$ class
Let $E$ be the event of selecting a vehicle whose tank is more than $50.5$ litres
Let $S$ be the sample space - the total number of vehicles.

$n(E) = 21 \\[3ex] n(S) = 150 \\[3ex] P(E) = \dfrac{n(E)}{n(S)} \\[5ex] P(E) = \dfrac{21}{150} \\[5ex] P(E) = \dfrac{7}{50} \\[5ex]$ (d.) Prerequisite to answer (c.) is: Measures of Center
$\Sigma f = 150 \\[3ex] \dfrac{\Sigma f}{2} = \dfrac{150}{2} = 75 \\[5ex] Begin\:\:from\:\:the\:\:first\:\:class \\[3ex] Keep\:\:adding\:\:the\:\:frequencies\:\:till\:\:you\:\:get\:\:to\:\:10 \\[3ex] 24 + 35 = 59 \\[3ex] 59 + 42 = 101...stop \\[3ex] Or\:\:better still \\[3ex] CF\:\:of\:\:3rd\:\:class = 101...stop \\[3ex]$ Therefore, the class that contains the median is the $3rd$ class = $31 - 40$
Byron's estimate of $43.5$ is in the $4th$ class = $41 - 50$. That is incorrect.


text:
시간 제한 메모리 제한 제출 정답 맞은 사람 정답 비율
2 초 256 MB 6 1 1 25.000%

문제

In 2000, the contest activities sort-of culminated at the Czech Technical University. There were another three competitions again (for the last time, probably), including the third Central Europe Regional Contest in a row. The region then moved to Warsaw and FEL++ was temporarily suspended.

We will end our small historical excursion here and present you the fifth and last problem that originates from our rich archives. It appeared in the 2000 Central Europe Regional Contest.

Building the GSM network is a very expensive and complex task. Moreover, after the Base Transceiver Stations (BTS) are built and working, we need to perform many various measurements to determine the state of the network, and propose effective improvements to be made.

The ACM technicians have a special equipment for measuring the strength of electro-magnetic fields, the transceiver power and signal quality. This equipment is packed into a huge knapsack and the technician must move with it from one BTS to another. Unfortunately, the knapsack has not enough memory for storing all of the measured values. It has a small cache only, that can store values for several seconds. Then the values must be transmitted to the BTS by an infrared connection (IRDA). The IRDA needs direct visibility between the technician and the BTS.

Your task is to find the path between two neighboring BTSes such that at least one them is always visible. For simplicity, a town is modelled as a rectangular grid of P × Q square fields. Each field is exactly 1 meter wide. For each field, a non-negative integer Zi,j is given, representing the height of the terrain in that place, in meters. That means the town model is made of cubes, each of them being either solid or empty.

The technician is moving in steps (steps stands for Standard Technician’s Elementary Positional Shift). Each step is made between two neighboring square fields in North, South, West or East direction, it is not possible to move diagonally. The step between two fields A and B (step from A to B) is allowed only if the height of the terrain in the field B is not very different from the height in the field A. The technician can climb at most 1 meter up or descend at most 3 meters down in a single step.

At the end of each step, at least one of the two BTSes must be visible. However, there can be some point “in the middle of the step” where no BTS is visible (data are handled by the cache). The BTS is considered visible if there is a direct visibility between the unit cube just above the terrain on the BTSes coordinates and the cube just above the terrain on the square field where the technician is. Direct visibility between two cubes means that the line segment connecting the centers of the two cubes does not intersect any solid cube. However, the line can touch any number of solid cubes. In other words, consider both the BTS and the technician being points exactly half meter above the surface and in the center of the appropriate square field.

Note that the IRDA beam can also go between two cubes that touch each other by their edge, although there is no real space between them. It is because such a beam touches both of these two cubes but does not intersect them.

입력

There is a single positive integer T on the first line of input. It stands for the number of test cases to follow. The first line of each test case contains two integer numbers P and Q, separated by a single space, 1 ≤ P, Q ≤ 200. Then there are P lines, each containing Q integer numbers separated by a space. These numbers are Zi,j , where 1 ≤ i ≤ P, 1 ≤ j ≤ Q and 0 ≤ Zi,j ≤ 5000. After the terrain description, there are four numbers R1, C1, R2, C2 on the last line of each test case. These numbers represent positions of two BTSes, 1 ≤ R1, R2 ≤ P, 1 ≤ C1, C2 ≤ Q. The first coordinate determines the row of the town, the second coordinate determines the column.

출력

You are to find the shortest possible path meeting the above criteria. That is, all steps must be done between neighboring fields, the terrain must not elevate or descend too much, and at the end of each step, at least one BTS must be visible.

For each test case, print one line containing the sentence “The shortest path is M steps long.”, where M is the number of steps that must be made. If there is no such path, output the sentence “Mission impossible!”.

예제 입력

4
5 5
8 7 6 5 4
2 2 2 2 2
2 2 2 2 2
2 2 2 2 2
2 2 2 2 2
1 1 5 1
5 8
2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2
9 9 9 9 9 9 9 2
2 2 2 2 2 2 2 2
1 2 5 1
5 8
2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2
9 9 9 9 9 9 9 2
2 2 2 2 2 2 2 2
1 5 5 1
6 12
5 5 5 5 1 5 5 5 5 5 5 5
5 5 5 5 1 5 5 5 5 5 5 5
5 5 5 5 9 5 5 5 5 5 5 5
5 9 1 5 5 5 5 5 5 5 5 5
5 5 9 5 5 5 5 5 5 5 5 5
5 5 9 5 5 5 5 5 5 5 5 5
6 1 3 12

예제 출력

The shortest path is 10 steps long.
Mission impossible!
The shortest path is 14 steps long.
The shortest path is 18 steps long.


text:
27 questions linked to/from List of primes under a million
303 views

Shortest test for primality for byte [duplicate]

Write a function which gets as input a number (0 <= n <= 255) stored in 32-bit variable and returns non-zero integer if n ...
44k views

Is this number a prime?

Believe it or not, we do not yet have a code golf challenge for a simple primality test. While it may not be the most interesting challenge, particularly for "usual" languages, it can be nontrivial in ...
5k views

Return the Closest Prime Number

Challenge This is a simple one: Given a positive integer up to 1,000,000, return the closest prime number. If the number itself is prime, then you should return that number; if there are two primes ...
2k views

Prime counting function

Introduction The Prime Counting Function, also known as the Pi function $\pi(x)$, returns the amount of primes less than or equal to x. Challenge Your program will take an integer x which you ...
3k views

Sum of primes between given range

Write the shortest code for finding the sum of primes between a and b (inclusive). Input a ...
9k views

List Prime Numbers [duplicate]

Introduction Prime numbers are simple, right? Well, now you get your chance to find out! Challenge You must write a program or function that takes an input n and ...
3k views

Calculate the sum of the first n prime numbers

I'm surprised that this challenge isn't already here, as it's so obvious. (Or I'm surprised I couldn't find it and anybody will mark it as a duplicate.) Task Given a non-negative integer $n$, ...
32k views

List of first n prime numbers most efficiently and in shortest code [closed]

Rules are simple: First n primes (not primes below n), should be printed to standard output separated by newlines (primes should be generated within the code) primes cannot be generated by an inbuilt ...
2k views

Smallest n-digit prime containing only these digits

You will need to generate the smallest prime with n digits, and it will only contain digits specified in the list k. Examples: ...
4k views

Kolmogorov-mania

The Kolmogorov complexity of a string s is defined as the length of the shortest program P that outputs s. If the length of P is shorter than the length of s, then s is said to be compressible, ...
1k views

Variable Prime “Twins”

I am 2/3 twins with my brother, i.e. born on the same day of the same month but twelve years later. When I was 5, he was 17, both primes; the last pair of ages we can reasonably count on is [71, 83] ...
2k views

A curious prime fraction formula

Given a positive integer n output the integers a and b (forming reduced fraction a/b) such that: Where pk is the k th prime number (with p1 = 2). Examples: ...
1k views

Prime Divisor Table

Intro Something I've played around with in recreational mathematics has been construction of a divisor table to visually compare/contrast the prime divisors of a set of numbers. The set of input ...


text:
# Math Help - Finding the roots of this equation: x^2 + (a - 1/a)x-1 = 0

1. ## Finding the roots of this equation: x^2 + (a - 1/a)x-1 = 0

Any help would be greatly appreciated.

2. Originally Posted by -DQ-
Any help would be greatly appreciated.

The roots of this equation are

$x= \frac{-b\pm\sqrt{b^2-4ac}}{2a}$

where a = 1, b= a-1/a, c=-1

3. Correct, but it is not a good idea to use "a" to represent two different quantities!

Perhaps it would be better to write
"The roots to the equation $ix^2+ jx+ k= 0$ are $x= \frac{-j\pm\sqrt{j^2- 4ik}}{2i}$. For this problem, i= 1, j= (a-1)/a, k= -1".

(I am assuming that in the first post, "a- 1/a" meant "(a-1)/a", not "a- (1/a)".)

4. This quadratic factorises as $(x+a)(x-\tfrac1a)$. So the solutions of the equation are $x=-a$ and $x-\tfrac1a$.

(I am making the opposite assumption from HallsofIvy, namely that, "a- 1/a" meant "a- (1/a)", not "(a-1)/a".)


text:
# Minimal embedding of Hunchentoot¶

The following code provides a concise example how to embed Hunchentoot into a running FriCAS instance. The GET and POST method will be demonstrated by using urls in a web-browser and a short Python code respectively.

## The code¶

Paste the following code into a file minserver.lisp

(load "~/quicklisp/setup")
(ql:quickload :hunchentoot)

(in-package :boot)

;;; Config
(defparameter +port+ 4242)

;;; SPAD eval

(defun spad_eval (code)
(let ((*package* (find-package :boot))
(alg (boot::|parseAndEvalToString| code)))
(format nil "~{~A~%~}" alg)))

;;; WEB server
(hunchentoot:define-easy-handler (fricas-eval :uri "/eval") (code)
(setf (hunchentoot:content-type*) "text/plain")
(format nil "~A~%" (spad_eval code)))

(hunchentoot:start (make-instance 'hunchentoot:easy-acceptor :port +port+))

;;; add :address "localhost"  if you wish local access only!
;;; test example: http://localhost:4242/eval?code=D(x^n,x,6)

## Load server¶

We will do all steps manually for the sake of clarity. Start FriCAS and load the file minserver.lisp:

\$ fricas

FriCAS Computer Algebra System
Version: FriCAS 2016-08-28
Timestamp: Sam Sep 17 00:34:49 CEST 2016
---------------------------------------------------------------------------
Issue )copyright to view copyright notices.
Issue )summary for a summary of useful system commands.
Issue )quit to leave FriCAS and return to shell.
---------------------------------------------------------------------------

(1) -> )lisp (load "minserver")
To load "hunchentoot":
Load 1 ASDF system:
hunchentoot
; Loading "hunchentoot"
............
Value = T
(1) ->

## Open a web browser¶

Now open a web browser (some terminal based browsers like w3m, lynx or links, do not accept all urls as e.g. Firefox does, however, when using quotes most urls will work). Enter the following url:

http://localhost:4242/eval?code=D(x^n,x,6)

The result in the browser window should exactly look like as below:

  6      5      4       3       2         n - 6
(n  - 15n  + 85n  - 225n  + 274n  - 120n)x
Type: Expression(Integer)

## Explanation¶

The connection of the url example above with the lisp code is almost obvious:

(fricas-eval :uri "/eval") (code)

The easy handler expects a query ? and one variable code. Then this query will be evaluated and the result formatted and written to the client.

## GET and POST¶

Two commonly used methods for a request-response between a client and server are:

GET and POST.

GET - Requests data from a specified resource
POST - Submits data to be processed to a specified resource

### The GET Method¶

Note that the query string (name/value pairs) is sent in the URL of a GET request:

/test/demo_form.asp?name1=value1&name2=value2

Restrictions: The GET method adds the data to the URL; and the length of a

URL is limited (maximum URL length is 2048 characters and ASCII characters only are allowed.

Certain character have a special meaning in an URL, so url encoding must follow some rules (HTML URL Encoding Reference).

### The POST Method¶

Note that the query string (name/value pairs) is sent in the HTTP message body of a POST request:

POST /test/demo_form.asp HTTP/1.1
Host: w3schools.com
name1=value1&name2=value2

The restrictions of the GET method do not apply here, however, to get a full overview consult the link above. There also are other HTTP request possibilities.

## Python example for POST¶

import requests
url = 'http://localhost:4242/eval'
payload = {'code': 'D(x^n,x,8)'}

r = requests.post(url, data=payload)

print(r.text)

Live action:

C:\Users\nilqed>python
Python 2.7.10 (default, May 23 2015, 09:40:32) [MSC v.1500 32 bit (Intel)]
Type "help", "copyright", "credits" or "license" for more information.
>>> import requests
>>> url = 'http://localhost:4242/eval'
>>> payload = {'code': 'D(x^n,x,8)'}
>>> r = requests.post(url, data=payload)
>>> r
<Response [200]>
>>> print(r.text)

8      7       6        5        4         3         2          n - 8
(n  - 28n  + 322n  - 1960n  + 6769n  - 13132n  + 13068n  - 5040n)x

Type: Expression(Integer)

>>>

With the POST method we can use high payloads and there are many languages with HTTP support libraries.

## Conclusion¶

This skeleton example shows the general method and of course only uses the most basic features of Hunchentoot. We can define many more easy handlers and other - even more sophisticated - methods to access the internals of FriCAS by HTTP requests. The sky is the limit ;-)


text:
# Powershell search millions of files as fast as possible

I once asked a similar question but in C#. Now I have the same problem in powershell..

What is the fastest way, to search files newer than 15 minutes, in a file system with more than 1 million files?

Is there any faster way than using pipe?

Get-ChildItem -Path $path -Recurse | Select Name, PSIsContainer, Directory, LastWriteTime, Length| where {($_.LastWriteTime -gt (Get-Date).AddMinutes(-15))}

I already cut off some attributes to minimize the object size. It still takes ages.

First, you don't need to call Get-Date for every file. Just call it once at the beginning:

$t = (Get-Date).AddMinutes(-15) Get-ChildItem -Path$path -Recurse |
Select Name, PSIsContainer, Directory, LastWriteTime, Length |
where {($_.LastWriteTime -gt$t)}

That's saves about 10% (as measured by Measure-Command).

Secondly, you don't need to call Select-Object for each file either. Just change the processing order:

$t = (Get-Date).AddMinutes(-15) Get-ChildItem -Path$path -Recurse |
where {($_.LastWriteTime -gt$t)} |
Select Name, PSIsContainer, Directory, LastWriteTime, Length

Thirdly, try increasing the buffer size using the OutBuffer parameter:

$t = (Get-Date).AddMinutes(-15) Get-ChildItem -Path$path -Recurse -OutBuffer 1000 |
where {($_.LastWriteTime -gt$t)} |
Select Name, PSIsContainer, Directory, LastWriteTime, Length

I've used 1000, but you can experiment with the value.

Those three changes reduced the running time to under one half on my system.

• Thank you. This saves a lot of time. I know the question is about powershell. Nevertheless I found an even faster solution -> Batch. Batch does the whole thing 10 times faster. – greenhoorn Jan 22 '15 at 12:08
• @greenhoorn, do you mean batch files (cmd.exe)? – Dangph Jan 22 '15 at 12:21
• Yes, I do mean batch files. – greenhoorn Jan 22 '15 at 12:23
• @greenhoorn Do you want to post your solution? Using Batch files, that is. – Anthony Horne Feb 7 '17 at 16:56
• @AnthonyHorne This post is over 2 y old. I'm not in possession of the solution anymore. I'm sure you will find it with google like I did. – greenhoorn Feb 8 '17 at 7:34

On this website, I found an example that I prefer to most of the solutions that are easily found with Google search.

The key is in using foreach-object instead of pipelining where-object filters and sort operations.

$newdate = [DateTime]::MinValue$newfn = ""
$path = "." get-childitem$path | ForEach-Object {
if ($_.LastWriteTime -gt$newdate -and -not $_.PSIsContainer) {$newfn = $_.Name$newdate = $_.LastWriteTime } }$output = ""
if ($newfn -ne "") {$output += "nNewest: " + $newdate + " -- " +$newfn }
if ($output -eq "") {$output += "nFolder is empty." }
\$output + "n"
`


text:
# wx.lib.docview.DocPrintout¶

DocPrintout is a default wx.Printout that prints the first page of a document view.

## Class Hierarchy¶

Inheritance diagram for class DocPrintout:

## Known Superclasses¶

wx.Printout

## Methods Summary¶

 __init__ Constructor. GetPageInfo Indicates that the DocPrintout only has a single page. GetView Returns the DocPrintout’s view. HasPage Indicates that the DocPrintout only has a single page. OnPrintPage Prints the first page of the view.

## Class API¶

class DocPrintout(wx.Printout)

DocPrintout is a default wx.Printout that prints the first page of a document view.

### Methods¶

__init__(self, view, title="Printout")

Constructor.

GetPageInfo(self)

Indicates that the DocPrintout only has a single page.

GetView(self)

Returns the DocPrintout’s view.

HasPage(self, pageNum)

Indicates that the DocPrintout only has a single page.

OnPrintPage(self, page)

Prints the first page of the view.


text:
## Understanding the Hilbert curve

The Hilbert curve fills space with good properties for sorting N-dimensional data in a linear fashion. I present an IPython notebook with the complete code to follow the algorithm of C. Hamilton CS-2006-07.

The notebook can be downloaded on my GitHub account: see the notebook | download the notebook (a right-click may help).

Update on 12 june 2015: addition of compact Hilbert indices.

In [17]:
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import math
from mpl_toolkits.mplot3d import Axes3D
from __future__ import print_function, division

# Understanding the Hilbert curve¶

Author: Pierre de Buyl
This notebook is originally developed in a github repository and presented on my blog. The code is 3-clause BSD and content CC-BY.

In this notebook, I review the algorithm by C. H. Hamilton to compute Hilbert curves and indices. The Hilbert curve is a space-filling curve that has good data locality to map N-dimensional data onto a linear space. It is based on a basis curve that visits the $2^N$ vertices of the N-dimensional binary space and that is replicated recursively $M$ times to map $2^{N M}$ points. The curve is the series of points in N-dimensional space and the index of a point is its place on this curve, computing it is subtle. The algorithm is presented in Compact Hilbert Indices, Dalhousie University Technical Report CS-2006-07, July 2006 (here after [TR]).

My motivations (and yours, maybe) are:

• Understand the algorithm to implement it in a simulation code. The point is to sort particles spatially for better memory access.
• Provide a full Python-based implementation that can be used for reference and illustration purposes.
• Have a clear explanation for later reference.

The notebook is self-contained: it contains all the code for all the intermediate operations defined in [TR]. Hold on for the graphic illustrations, they will come as the algorithm is transformed into code.

## Bitwise operations¶

I review here the bitwise operations in Python, as they will be used throughout. Only integer numbers are used. We define a word length that corresponds to the dimension of the Hilbert curve.

The examples use N=3 so that each single integer in $[0:2^N-1]$ can be mapped to a point in 3 dimension that is one vertex of the unit cube. Bitwise operations act as transformations on the integers and equivalently as geometrical operations in 3 dimensions. The XOR operator represents a reflection with respect to one or several of the symmetry planes of the cube. The rotation operator represents a rotations of a third around the diagonal, i.e. the line joining (0,0,0) and (1,1,1).

The function bin_str helps to visualize bitwise operations by converting an integer to its three bits. The right-most bit represents the x-axis, the middle bit the y-axis and the left-most bit the z-axis. In Python, the AND &, OR | and XOR ^ operators are defined and also the bit shifts (>> to the right, or less significant bits and << to the left or more significant bits). The NOT operator ~ requires further masking: It changes the sign bit of the integer. We wish here to operate on a N-bit space and thus cancel any bit beyond the N less significant bits of the integer. The fact that we remain with N non-zero bits and that NOT is operating properly is checked below. We also need rotation operators (like the shift operators but that transfer the over- and under-flowing bits around the 3-bit space) and a function to extract the i-th bit from a number, bit_component.

In [2]:
N = 3

def bin_str(i):
"""Return a string representation of i with N bits."""
out = ''
for j in range(N-1,-1,-1):
if (i>>j) & 1 == 1:
out += '1'
else:
out += '0'
return out

def rotate_right(x, d):
"""Rotate x by d bits to the right."""
d = d % N
out = x >> d
for i in range(d):
bit = (x & 2**i)>>i
out |= bit << (N+i-d)
return out

def rotate_left(x, d):
"""Rotate x by d bits to the left."""
d = d % N
out = x << d
excess = out
out = out & (2**N-1)
for i in range(d):
bit = (x & 2**(N-1-d+1+i))>> (N-1-d+1+i)
out |= bit << i
return out

def bit_component(x, i):
"""Return i-th bit of x"""
return (x & 2**i) >> i
In [3]:
print('AND  ', bin_str(1), '&', bin_str(5), '=', bin_str(1 & 5))
print('OR   ', bin_str(1), '|', bin_str(5), '=', bin_str(1 | 5))
print('XOR  ', bin_str(1), '^', bin_str(5), '=', bin_str(1 ^ 5))
print('shift', bin_str(2), '>> 1  =', bin_str(2 >> 1))
print('shift', bin_str(2), '<< 1  =', bin_str(2 << 1))
print('rot  ', bin_str(1), '↻', bin_str(rotate_right(1, 1)), '↻',
bin_str(rotate_right(1,2)), '↻', bin_str(rotate_right(1,3)))
print('NOT  ', bin_str(1), '      =', bin_str(~1 & 2**N-1))

# verify that '~i & 2**N-1' performs the NOT operation in N-bit space
for i in range(2**N):
not_i = ~i & 2**N-1
assert not_i >=0
assert not_i < 2**N
assert i & not_i == 0
assert i | not_i == 2**N-1
AND   001 & 101 = 001
OR    001 | 101 = 101
XOR   001 ^ 101 = 100
shift 010 >> 1  = 001
shift 010 << 1  = 100
rot   001 ↻ 100 ↻ 010 ↻ 001
NOT   001       = 110

The technical reports describes with good details all the intermediate steps of the algorithm using binary operations and sequences. The gray code index of a number i, for instance, is given by gc(i) = i XOR (i ≫ 1).

To understand the building of the Hilbert curve, other quantities are needed. The graphical representation use arrows in each sub-hypercube of the curve to represent the path. The definitions of the gray curve index and of the arrows is given and we use these defintions to illustrate the Hilbert curve afterwards.

The interpretation of gc, e and f below is a binary representation in base N of the vertices of a hypercube (a square for N=2, a cube for N=3, etc). The bits represent the z, y and x coordinates.

In [4]:
# Build all edges of a square and of a cube

edges = [((0,),(1,))]

"""Extend the list of edges from an hypercube to the list of
edges of the hypercube of the next dimension."""

old_edges = list(edges)
old_points = set( [x[0] for x in old_edges] ) | set( [x[1] for x in old_edges] )
edges = [ ( (0,)+x[0], (0,)+x[1] ) for x in old_edges ]
edges.extend( [ ( (1,)+x[0], (1,)+x[1] ) for x in old_edges ] )
for e in old_points:
edges.append( ( (0,)+e, (1,)+e) )

return edges

cube_xyz = []
for single_edge in cube_edges:
cube_xyz.append(single_edge[0])
cube_xyz.append(single_edge[1])
# The np.inf value breaks the line plot into disconnected parts
cube_xyz.append((np.inf,)*3)
cube_xyz = np.array(cube_xyz)

def set_unit_cube(ax, side=1, set_view=(10,-67)):
"""Present the unit cube."""
ax.set_xlabel('x'); ax.set_ylabel('y'); ax.set_zlabel('z');
ax.set_xticks(range(side+1)); ax.set_yticks(range(side+1)); ax.set_zticks(range(side+1));
ax.set_xlim(0, side); ax.set_ylim(0,side); ax.set_zlim(0,side);
if set_view:
ax.view_init(*set_view)

fig = plt.figure(figsize=(6, 5))

ax.plot(*zip(*cube_xyz), color='b')

for i in range(2**N):
x = bit_component(i, 0)
y = bit_component(i, 1)
z = bit_component(i, 2)
ax.text(x, y+0.05, z+0.05, bin_str(i))

set_unit_cube(ax)
In [5]:
def gc(i):
"""Return the Gray code index of i."""
return i ^ (i >> 1)

def e(i):
"""Return the entry point of hypercube i."""
if i==0:
return 0
else:
return gc(2*int(math.floor((i-1)//2)))

def f(i):
"""Return the exit point of hypercube i."""
return e(2**N-1-i) ^ 2**(N-1)

def i_to_p(i):
"""Extract the 3d position from a 3-bit integer."""
return [bit_component(i,j) for j in (0,1,2)]
In [6]:
print('Gray code', map(lambda i: (gc(i), bin_str(gc(i))), range(2**N)),
'Entry point', map(lambda i: (e(i), bin_str(e(i))), range(2**N)),
'Exit point', map(lambda i: (f(i), bin_str(f(i))), range(2**N)),
sep='\n')
Gray code
[(0, '000'), (1, '001'), (3, '011'), (2, '010'), (6, '110'), (7, '111'), (5, '101'), (4, '100')]
Entry point
[(0, '000'), (0, '000'), (0, '000'), (3, '011'), (3, '011'), (6, '110'), (6, '110'), (5, '101')]
Exit point
[(1, '001'), (2, '010'), (2, '010'), (7, '111'), (7, '111'), (4, '100'), (4, '100'), (4, '100')]

## Visiting the Hilbert curve in the cube¶

The path of the Hilbert curve without recursion is encoded in the sequence gc(i), whose 3-bit values are interpreted a points in 3D (see above), and that visits the edges of the cube. These coordinates (scaled by 1/2) are stored in xyz.

The first figure displays the curve visiting the centers of the subcubes, one of the typical representations of the Hilbert curve. In [TR], the path of the Hilbert curve proceeds by entering subcube gc(i) at point e(i) and leaving it a point f(i). The entry and exit points e(i) and f(i) are stored in e_xyz and f_xyz and shown below. The last figure represents the full path as the sequences of arrows joining those points.

In [7]:
# Obtain the 3 d coordinates for gc(i), e(i) and f(i) to build the Hilbert curve
xyz = np.array(map(lambda i: i_to_p(gc(i)), range(2**N)))/2
e_xyz = np.array(map(lambda i: i_to_p(e(i)), range(2**N)))/2
f_xyz = np.array(map(lambda i: i_to_p(f(i)), range(2**N)))/2

fig = plt.figure(figsize=(14, 12))

ax1.plot(*zip(*(xyz+0.25)))
for i, (x,y,z) in enumerate((xyz+0.27)):
ax1.text(x, y, z, str(i))
ax1.set_title('The Hilbert curve, subcube centered')

edges = []
for e_point, f_point, origin in zip(e_xyz, f_xyz, xyz):
dir = (f_point-e_point)
# For the 3d quiver plot, the origin of the arrow is shifted
# which is why we use origin+e_point+dir instead of origin+e_point
edges.append(np.concatenate((origin+e_point+dir, dir)))

ax2.quiver(*zip(*edges), length=0.5)
ax2.set_title('Arrows to visit the Hilbert curve')

# A naive sequence of binary points in 3D
c = 0.25 + np.array(map(i_to_p, range(2**N)))/2

ax3.plot(*zip(*c), color='b')
for i, (x,y,z) in enumerate(c):
ax3.text(x, y, z, str(i))

ax3.set_title('A naive curve to visit the cube')

for ax in [ax1, ax2, ax3]:
ax.plot(*zip(*cube_xyz), color='k', zorder=1, alpha=0.5)
set_unit_cube(ax)

## Define the curve directions and transforms¶

To write the Hilbert curve algorithm, we still need to define the inter-subcube direction g and the intra-subcube direction d (the direction of the arrow that connect e and f).

In [8]:
def inverse_gc(g):
"""The inverse gray code."""
i = g
j = 1
while j<N:
i = i ^ (g >> j)
j = j + 1
return i

def g(i):
"""The direction between subcube i and the next one"""
return int(np.log2(gc(i)^gc(i+1)))

def d(i):
"""The direction of the arrow whithin a subcube."""
if i==0:
return 0
elif (i%2)==0:
return g(i-1) % N
else:
return g(i) % N

def T(e, d, b):
"""Transform b."""
out = b ^ e
return rotate_right(out, d+1)

def T_inv(e, d, b):
"""Inverse transform b."""
return T(rotate_right(e, d+1), N-d-2, b)

## Validating the relations¶

The algorithm for the Hilbert curve is based on a series of relations between e, f, d and g. We can check those easily now. In order, we verify that:

• Only the g(i)-th bit is changed between gc(i) and gc(i+1)
• The tranform T sends an entry point e(i) to zero and the exit point f(i) to $2^{N-1}$: $T_{e(i),d(i)}(e(i))=0$ and $T_{e(i),d(i)}(f(i))=f(i)$
• The entry point e(i) reflected along the direction d(i) is the exit point f(i): $e(i)$ ^ $d(i) = f(i)$
• The inverse transform $T^{-1}$ composed with the direct transform $T$ is the identity.
In [9]:
# only g(i)-th bit changes from gc(i) to gc(i+1)
for i in range(2**N-1):
assert gc(i) ^ (1 << g(i)) == gc(i+1)

# T sends e(i) to 0 and f(i) to 2**(N-1)
for i in range(2**N):
assert T(e(i), d(i), e(i))==0
assert T(e(i), d(i), f(i))==2**(N-1)

# e(i) reflected in direction d(i) is f(i)
for i in range(2**N):
assert e(i) ^ (1<<d(i)) == f(i)

# T_inv composed with T (and vice versa) is the identity operator
for i in range(2**N):
for b in range(2**N):
assert T(e(i), d(i), T_inv(e(i), d(i), b)) == b
assert T_inv(e(i), d(i), T(e(i), d(i), b)) == b

## Transform the curve for the subcubes¶

According to [TR], one can use the transform T to map a subcube into the larger cube. Using the inverse transform, it is thus possible to generate the recursed Hilbert curve by mapping the main curve into each subcube. The coordinates of this refined Hilbert curve lie within the subcubes and must be translated by the corresponding origins to fill space correctly.

We construct below the refined Hilbert curve by this manual procedure and display the individual subcubes. The coordinates are then used to show the directed path followed by the curve.

In [10]:
fig = plt.figure(figsize=(12, 18))

# store the Hilbert curve with no recursion
main_curve = [gc(i) for i in range(2**N)]
# list of the in-subcube curves
sub_curves = []
for i, h_idx in enumerate(main_curve):
# append points from the subcube (obtained with the inverse transform T_inv
points = np.array( [ i_to_p(T_inv(e(i), d(i), x)) for x in main_curve ] )
# Translate the points by the origin of the subcube
points += np.array(i_to_p(h_idx))*2
sub_curves.append( points )

# plot all subcurves separately
[ax.plot(*zip(*c)) for c in sub_curves]
ax.set_title('Individual Hilbert "subcubes"')
set_unit_cube(ax, 3,  set_view=False)

# prepare the data to draw arrows with the quiver command
cc = np.concatenate(sub_curves)
U = cc[1:,0]-cc[:-1,0]
X = cc[:-1,0]+U
V = cc[1:,1]-cc[:-1,1]
Y = cc[:-1,1]+V
W = cc[1:,2]-cc[:-1,2]
Z = cc[:-1,2]+W
ax.quiver(X,Y,Z,U,V,W,color=[matplotlib.cm.gnuplot(x) for x in np.linspace(0, 1,len(X)*3)], length=1)
ax.set_title('Connected Hilbert "subcubes"')
set_unit_cube(ax, 3, set_view=False)

## Implementation of the algorithms¶

Now that all the components of the algorithm have been reviewed, it is time for the actual code. In [TR], the algorithm is listed explicitly in pseudo-code.

• Algorithm 2, implemented in TR_algo2, transforms a set of integer coordinates into its Hilbert index.
• Algorithm 3, implemented in TR_algo3, transforms a Hilbert index into a set of coordinates.

Algorithm 2 works by computing the Hilbert index at subsequent refinement levels. Starting from the most significant bits of the coordinates p, it computes the subcube to which p belongs as a label l. l encodes in its N bits the subcube of interest as was show in the unit cube earlier. From l, it is easy to obtain the Hilbert index. The final Hilbert index is composed of this data packed in a single integer variable. As it proceeds in the subcube, the algorithm composes the transform T by updating the current entry point ve direction vd according to Lemma 2.13 of [TR].

Algorithm 3 works in a similar manner but extracts the relevant bits of the Hilbert index at each refinement level to construct the coordinates of the point.

In [11]:
M = 2
def TR_algo2(p):
"""Return the Hilbert index of point p"""
# h will contain the Hilbert index
h = 0
# ve and vd contain the entry point and dimension of the current subcube
# we choose here a main traversal direction N-2 (i.e. z for a cube) to match
# the illustrations
ve = 0
vd = 2
for i in range(M-1, -1, -1):
# the cell label is constructed in two steps
# 1. extract the relevant bits from p
l = [bit_component(px, i) for px in p]
# 2. construct a integer whose bits are given by l
l = sum( [lx*2**j for j, lx in enumerate(l)] )
# transform l into the current subcube
l = T(ve, vd, l)
# obtain the gray code ordering from the label l
w = inverse_gc(l)
# compose (see [TR] lemma 2.13) the transform of ve and vd
# with the data of the subcube
ve = ve ^ (rotate_left(e(w), vd+1))
vd = (vd + d(w) + 1) % N
# move the index to more significant bits and add current value
h = (h << N) | w
return h

def TR_algo3(h):
"""Return the coordinates for the Hilbert index h"""
ve = 0
vd = 2
p = [0]*N
for i in range(M-1, -1, -1):
w = [bit_component(h, i*N+ii) for ii in range(N)]
#print(i, w)
w = sum( [wx*2**j for j, wx in enumerate(w)] )
#print(i, w, gc(w))
l = gc(w)
l = T_inv(ve, vd, l)
for j in range(N):
p[j] += bit_component(l, j) << i
ve = ve ^ rotate_left(e(w), vd+1)
vd = (vd + d(w) + 1) % N
return p

## Displaying and checking the algorithm¶

A basic consistency check is to verify that the successive application of TR_algo3 and TR_algo2 returns the same number than the one given.

A visual inspection of the Hilbert curve on which we superimpose the Hilbert indices confirms the consistency of the algorithm.

In [12]:
# Verifying that the algorithm and its inverse agree
for h_idx in range(2**(N*M)):
assert TR_algo2(TR_algo3(h_idx)) == h_idx
In [13]:
fig = plt.figure(figsize=(12,9))

ax.quiver(X,Y,Z,U,V,W,color=[matplotlib.cm.gnuplot(x) for x in np.linspace(0, 1,len(X)*3)], length=1)
ax.set_title('Hilbert indices in 3D space')

for h in range(2**(N*M)):
x, y, z = TR_algo3(h)
ax.text(x, y, z+0.1, str(h))

set_unit_cube(ax, 3, set_view=(15,-55))

## Compact indices¶

One novelty in [TR] is the definition of a compact Hilbert index that enables efficient storage of indices for spaces that have unequal sides. Think of it as allowing non-cubic boxes in 3D.

The ideas underlying this compact index are:

• The ordering of the index with respect to the full Hilbert index should be conserved.
• Some bits will be meaningless to store as they always be equal to zero and as such they will be skipped. In practice, it means that depending on the refinement (or recursion) level, less than N bits may be stored per level.

The algorithms below are very similar to the previous ones, with the addition of masking to retain only meaningful bits in the index. After the implementations, we perform a check that the direct and inverse algorithms agree and illustrate the final result.

In [14]:
def gcr(i, mu, pi):
"""Compute the gray code rank of i given the mask mu.
Algorithm 4 in [TR]"""
r = 0
for k in range(N-1, -1, -1):
if bit_component(mu, k):
r = (r << 1) | bit_component(i, k)
return r

def gcr_inv(r, mu, pi):
"""Inverse of the gray code rank, given the mask mu and the pattern pi.
Algorithm 5 in [TR]"""
i = 0
g = 0
j = sum([bit_component(mu, k) for k in range(N)])-1
for k in range(N-1, -1, -1):
if bit_component(mu, k)==1:
i |= bit_component(r, j) << k
g |= ( (bit_component(i, k) + bit_component(i, k+1))%2 ) << k
j -= 1
else:
g |= bit_component(pi, k) << k
i |= ( (bit_component(g, k) + bit_component(i, k+1)) % 2) << k
return i

# definition of the size of the space. compact_M is of length N
compact_M = [3, 2, 2]

"""Extract the mask for iteration i of the algorithm.
Algorithm 6 in [TR]"""
mu = 0
for j in range(N-1, -1, -1):
mu = mu << 1
if compact_M[j] > i:
mu = mu | 1
return mu

def TR_algo7(p):
"""Compute the compact Hilbert index for point p"""
h = 0
ve = 0
vd = 2
m = max(compact_M)
for i in range(m-1, -1, -1):
mu_norm = sum([bit_component(mu, j) for j in range(N)])
mu = rotate_right(mu, vd+1)
pi = rotate_right(ve, vd+1) & ((~mu) & 2**N-1)
l = [bit_component(px, i) for px in p]
# 2. construct a integer whose bits are given by l
l = sum( [lx*2**j for j, lx in enumerate(l)] )
l = T(ve, vd, l)
w = inverse_gc(l)
r = gcr(w, mu, pi)
ve = ve ^ rotate_left(e(w), vd+1)
vd = (vd + d(w) + 1) % N
h = (h << mu_norm) | r
return h

def TR_algo8(h):
"""Compute the point with compact Hilbert index h"""
ve = 0
vd = 2
k = 0
p = [0,]*N
m = max(compact_M)
vM = sum(compact_M)
for i in range(m-1, -1, -1):
mu_norm = sum([bit_component(mu, j) for j in range(N)])
mu = rotate_right(mu, vd+1)
pi = rotate_right(ve, vd+1) & (~mu & 2**N-1)
r = [bit_component(h, vM - k - (j+1)) for j in range(mu_norm)][::-1]
r = sum( [rx*2**j for j, rx in enumerate(r)] )
k = k + mu_norm
w = gcr_inv(r, mu, pi)
l = gc(w)
l = T_inv(ve, vd, l)
for j in range(N):
p[j] |= bit_component(l, j) << i
ve = ve ^ (rotate_left(e(w), vd+1))
vd = (vd + d(w) + 1) % N
return p
In [15]:
# Verifying that the algorithm and its inverse agree

for h_idx in range(2**sum(compact_M)):
assert TR_algo7(TR_algo8(h_idx))==h_idx
In [16]:
fig = plt.figure(figsize=(14, 12))

test_data = [TR_algo8(i) for i in range(2**sum(compact_M))]
ax.plot(*zip(*test_data))

for i in range(2**sum(compact_M)):
x, y, z = TR_algo8(i)
ax.text(x, y, z, str(i))

ax.set_title('Compact Hilbert indices in 3D space')
ax.set_xlabel('x'); ax.set_ylabel('y'); ax.set_zlabel('z');
ax.set_xticks(range(2**compact_M[0])); ax.set_yticks(range(2**compact_M[1])); ax.set_zticks(range(2**compact_M[2]));
ax.set_xlim(0, 2**compact_M[0]-1); ax.set_ylim(0,2**compact_M[1]-1); ax.set_zlim(0,2**compact_M[2]-1);
ax.view_init(10, -74);
In [ ]: