text:
## Saturday, January 4, 2020

### Small Blending Problem in PuLP

In [1] a small blending problem is proposed.

There are different raw steel materials that have to be mixed (blended) into a final material that has certain specifications. In this case the specifications are limits on the elements Carbon (C), Copper (Cu) and Manganese (Mn). We assume things blend linearly.

Blending is a traditional linear programming application, and models are found in many text books.

The problem is small so let's try PuLP here.

I'll try to write a bit about indexing using strings (instead of integers), and compare PuLP with CVXPY and GAMS. As the model is small, I'll also add some data manipulation (using data frames) and some simple reporting (also using data frames). I use string indexing in Pulp (with GAMS this is standard). Of course CVXPY is very different: it is matrix based. So there things are position dependent. The idea it to extract data from the data frame in such a way that positions are predictable.

#### Problem data

The data for the problem is as follows:

Demand: 5000 Kg

Specification of final material:

Element      %Minimum %Max
Carbon       2         3
Copper       0.4       0.6
Manganese    1.2       1.65

Raw material inventory:

Alloy          C%   Cu%   Mn%     Stocks kg Price € / kg
Iron alloy     2.50 0.00  1.30    4000      1.20
Iron alloy     3.00 0.00  0.80    3000      1.50
Iron alloy     0.00 0.30  0.00    6000      0.90
Copper alloy   0.00 90.00 0.00    5000      1.30
Copper alloy   0.00 96.00 4.00    2000      1.45
Aluminum alloy 0.00 0.40  1.20    3000      1.20
Aluminum alloy 0.00 0.60  0.00   2,500      1.00

#### Mathematical Model

The basic model is:

Blending Model
\begin{align}\min & \sum_i \color{darkblue}{\mathit{Cost}}_i\cdot \color{darkred}{\mathit{Use}}_i\\ & \color{darkblue}{\mathit{Min}}_j \le \frac{\displaystyle\sum_i \color{darkblue}{\mathit{Element}}_{i,j}\cdot \color{darkred}{\mathit{Use}}_i}{\displaystyle \sum_i \color{darkred}{\mathit{Use}}_i} \le \color{darkblue}{\mathit{Max}}_j \\ & \sum_i \color{darkred}{\mathit{Use}}_i = \color{darkblue}{\mathit{Demand}} \\ & 0 \le \color{darkred}{\mathit{Use}}_i \le \color{darkblue}{\mathit{Available}}_i \end{align}

Here we use: $\begin{cases} i & \text{Types of raw material in stock}\\ j & \text{Element with limits}\\ {\mathit{Cost}}_i & \text{Unit cost of raw material} \\ {\mathit{Min}}_j, {\mathit{Max}}_j & \text{Limits on element content in final product} \\ {\mathit{Element}}_{i,j} & \text{Content of elements in raw material}\\ {\mathit{Demand}} & \text{Demand for final product}\\ {\mathit{Available}}_i & \text{Availability of raw material} \\ {\mathit{Use}}_i & \text{Decision variable: how much raw material to use} \end{cases}$ The blending constraint is nonlinear: we divide by the total weight of the final product to calculate the percentages. We can linearize this fraction in two ways:

1. multiply all sides by $$\sum_i \mathit{Use}_i$$. This leads to $\mathit{Min}_j \cdot\sum_i \mathit{Use}_i \le \sum_i \mathit{Element}_{i,j}\cdot\mathit{Use}_i \le \mathit{Max}_j \sum_i \mathit{Use}_i$ This first reformulation is especially useful when the total final product is not constant. Note that this formulation is sometimes difficult to recognize when rewritten as something like: $\begin{cases} \displaystyle \sum_i (\mathit{Element}_{i,j}-\mathit{Min}_j )\mathit{Use}_i \ge 0 \>\>\forall j \\ \displaystyle \sum_i (\mathit{Max}_j - \mathit{Element}_{i,j})\mathit{Use}_i \ge 0 \>\>\forall j \end{cases}$
2. In our case we know that $$\sum_i \mathit{Use}_i$$ is constant: it is always equal to $$\mathit{Demand}$$. I.e. $\mathit{Min}_j \le \frac{1}{\mathit{Demand}} \sum_i \mathit{Element}_{i,j}\cdot\mathit{Use}_i \le \mathit{Max}_j$

We often need to split sandwich equations into two simple inequalities. That often leads to duplicate expressions: $\begin{cases} \displaystyle \frac{1}{\mathit{Demand}}\sum_i \mathit{Element}_{i,j}\cdot\mathit{Use}_i \ge \mathit{Min}_j \>\>\forall j \\ \displaystyle \frac{1}{\mathit{Demand}}\sum_i \mathit{Element}_{i,j} \cdot\mathit{Use}_i \le \mathit{Max}_j \>\>\forall j \end{cases}$ For small problems this is not an issue. For larger problems, I prefer to introduce extra variables that prevents these duplicate expressions. We end up with the following linear programming formulation:

Linear Programming Formulation
\begin{align}\min & \sum_i \color{darkblue}{\mathit{Cost}}_i \cdot\color{darkred}{\mathit{Use}}_i\\ & \color{darkblue}{\mathit{Demand}} \cdot \color{darkred}{\mathit{Content}}_j = \sum_i \color{darkblue}{\mathit{Element}}_{i,j} \cdot\color{darkred}{\mathit{Use}}_i \\ & \sum_i \color{darkred}{\mathit{Use}}_i = \color{darkblue}{\mathit{Demand}} \\ & \color{darkred}{\mathit{Use}}_i \in [0, \color{darkblue}{\mathit{Available}}_i]\\ & \color{darkred}{\mathit{Content}}_j \in [\color{darkblue}{\mathit{Min}}_j,\color{darkblue}{\mathit{Max}}_j] \end{align}

Even for small, almost trivial models, it makes sense to first develop a mathematical model. Especially, if you are not very experienced in developing linear programming models. Starting with a pen and a piece of paper is sometimes better than immediately start coding.

#### Implementation in Python/Pulp

An implementation using PuLP can look like:

from io import StringIO
import pandas as pd
import pulp as lp

# for inputting tabular data below
def table(s):

#------------------------------------------------------------------
# data
#------------------------------------------------------------------

demand = 5000

requirements = table("""
ID  Element      Min   Max
C   Carbon       2     3
Cu  Copper       0.4   0.6
Mn  Manganese    1.2   1.65
""")

supplyData = table("""
ID  Alloy             C       Cu     Mn     Stock   Price
A   "Iron alloy"      2.50    0.00   1.30   4000    1.20
B   "Iron alloy"      3.00    0.00   0.80   3000    1.50
C   "Iron alloy"      0.00    0.30   0.00   6000    0.90
D   "Copper alloy"    0.00   90.00   0.00   5000    1.30
E   "Copper alloy"    0.00   96.00   4.00   2000    1.45
F   "Aluminum alloy"  0.00    0.40   1.20   3000    1.20
G   "Aluminum alloy"  0.00    0.60   0.00   2500    1.00
""")

print("----- Data-------")
print(requirements)
print(supplyData)

#------------------------------------------------------------------
# derived data
#------------------------------------------------------------------

# our sets are stockItems ["A","B",..] and elements ["C","Cu",...]
Items = supplyData.index
Elements = requirements.index

print("----- Indices-------")
print(Items)
print(Elements)

#------------------------------------------------------------------
# LP Model
#------------------------------------------------------------------

use = lp.LpVariable.dicts("Use",Items,0,None,cat='Continuous')
content = lp.LpVariable.dicts("Content",Elements,0,None,cat='Continuous')

model = lp.LpProblem("Steel", lp.LpMinimize)

# objective : minimize cost
model += lp.lpSum([use[i]*supplyData.loc[i,'Price'] for i in Items ])

# upper bounds wrt availability
for i in Items:
model += use[i] <= supplyData.loc[i,'Stock']

# final content of elements and their bounds
for j in Elements:
model += demand*content[j] == lp.lpSum([use[i]*supplyData.loc[i,j] for i in Items])
model += content[j] >= requirements.loc[j,'Min']
model += content[j] <= requirements.loc[j,'Max']

# meet demand
model += lp.lpSum([use[i] for i in Items]) == demand

# for debugging
#print(model)

#------------------------------------------------------------------
# Solve and reporting
#------------------------------------------------------------------

model.solve()

print("----- Model Results-------")
print("Status:", lp.LpStatus[model.status])
print("Objective:",lp.value(model.objective))

# collect results
L = []
for i in Items:
L.append(['use',i,0.0,use[i].varValue,supplyData.loc[i,'Stock']])
for j in Elements:
L.append(['content',j,requirements.loc[j,'Min'],content[j].varValue,requirements.loc[j,'Max']])
results = pd.DataFrame(L,columns=['Variable','Index','Lower','Value','Upper'])
print(results)

Notes:

• We input the basic data as data frames. Data frames are a standard way to handle tabular data. Data frames are originally from the R statistical software system.
• Usually read_csv is for CSV files. Here we use it to read from a string. Blanks are used as separator to make the table more readable for humans.
• For each data frame we added an index column. This index will allow us to select a row from the data frame. Note that the index is a string. In general using strings as index is safer than using an index number. We see much earlier that things are wrong when making a mistake like using $$j$$ (element) instead of $$i$$ (raw material).
• Python Pandas allows duplicate indices. We can check for this using duplicated() function.
• Because we access the data by name, it would not matter if the rows or columns are in a different position. This is more like a database table, where we assume no particular ordering.
• We also use a data frame for reporting. Data frames are printed in a nicer way than Python arrays, and they can be exported to CVS files or spreadsheet with one function call.
• The variables are also indexed by names. This is accomplished by lp.LpVariable.dicts(). This is safer than using a standard array of variables.
• AFAIK, PuLP can only handle a single scalar bound in the LpVariable statement (e.g. all lower bounds for a variables are zero). This means: we have to specify a number of bounds as explicit singleton constraints or use a var.bounds statement.

The results look like:

----- Data-------
Element  Min   Max
ID
C      Carbon  2.0  3.00
Cu     Copper  0.4  0.60
Mn  Manganese  1.2  1.65
Alloy    C    Cu   Mn  Stock  Price
ID
A       Iron alloy  2.5   0.0  1.3   4000   1.20
B       Iron alloy  3.0   0.0  0.8   3000   1.50
C       Iron alloy  0.0   0.3  0.0   6000   0.90
D     Copper alloy  0.0  90.0  0.0   5000   1.30
E     Copper alloy  0.0  96.0  4.0   2000   1.45
F   Aluminum alloy  0.0   0.4  1.2   3000   1.20
G   Aluminum alloy  0.0   0.6  0.0   2500   1.00
----- Indices-------
Index(['A', 'B', 'C', 'D', 'E', 'F', 'G'], dtype='object', name='ID')
Index(['C', 'Cu', 'Mn'], dtype='object', name='ID')
----- Model Results-------
Status: Optimal
Objective: 5887.57427835
Variable Index  Lower        Value    Upper
0      use     A    0.0  4000.000000  4000.00
1      use     B    0.0     0.000000  3000.00
2      use     C    0.0   397.763020  6000.00
3      use     D    0.0     0.000000  5000.00
4      use     E    0.0    27.612723  2000.00
5      use     F    0.0   574.624260  3000.00
6      use     G    0.0     0.000000  2500.00
7  content     C    2.0     2.000000     3.00
8  content    Cu    0.4     0.600000     0.60
9  content    Mn    1.2     1.200000     1.65

#### Safety

We will get an error if we misspell things. E.g. if we use in the second table Mgn instead of Mn, we will see:

KeyError                                  Traceback (most recent call last)
/usr/local/lib/python3.6/dist-packages/pandas/core/indexes/base.py in get_loc(self, key, method, tolerance)
2896             try:
-> 2897                 return self._engine.get_loc(key)
2898             except KeyError:

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc()

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

KeyError: 'Mn'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
9 frames
<ipython-input-3-d87ed8f34980> in <module>()
59 # final content of elements and their bounds
60 for j in Elements:
---> 61   model += demand*content[j] == lp.lpSum([use[i]*supplyData.loc[i,j] for i in Items])
62   model += content[j] >= requirements.loc[j,'Min']
63   model += content[j] <= requirements.loc[j,'Max']

<ipython-input-3-d87ed8f34980> in <listcomp>(.0)
59 # final content of elements and their bounds
60 for j in Elements:
---> 61   model += demand*content[j] == lp.lpSum([use[i]*supplyData.loc[i,j] for i in Items])
62   model += content[j] >= requirements.loc[j,'Min']
63   model += content[j] <= requirements.loc[j,'Max']

/usr/local/lib/python3.6/dist-packages/pandas/core/indexing.py in __getitem__(self, key)
1416                 except (KeyError, IndexError, AttributeError):
1417                     pass
-> 1418             return self._getitem_tuple(key)
1419         else:
1420             # we by definition only have the 0th axis

/usr/local/lib/python3.6/dist-packages/pandas/core/indexing.py in _getitem_tuple(self, tup)
803     def _getitem_tuple(self, tup):
804         try:
--> 805             return self._getitem_lowerdim(tup)
806         except IndexingError:
807             pass

/usr/local/lib/python3.6/dist-packages/pandas/core/indexing.py in _getitem_lowerdim(self, tup)
959                     return section
960                 # This is an elided recursive call to iloc/loc/etc'
--> 961                 return getattr(section, self.name)[new_key]
962
963         raise IndexingError("not applicable")

/usr/local/lib/python3.6/dist-packages/pandas/core/indexing.py in __getitem__(self, key)
1422
1423             maybe_callable = com.apply_if_callable(key, self.obj)
-> 1424             return self._getitem_axis(maybe_callable, axis=axis)
1425
1426     def _is_scalar_access(self, key: Tuple):

/usr/local/lib/python3.6/dist-packages/pandas/core/indexing.py in _getitem_axis(self, key, axis)
1848         # fall thru to straight lookup
1849         self._validate_key(key, axis)
-> 1850         return self._get_label(key, axis=axis)
1851
1852

/usr/local/lib/python3.6/dist-packages/pandas/core/indexing.py in _get_label(self, label, axis)
154             # but will fail when the index is not present
155             # see GH5667
--> 156             return self.obj._xs(label, axis=axis)
157         elif isinstance(label, tuple) and isinstance(label[axis], slice):
158             raise IndexingError("no slices here, handle elsewhere")

/usr/local/lib/python3.6/dist-packages/pandas/core/generic.py in xs(self, key, axis, level, drop_level)
3735             loc, new_index = self.index.get_loc_level(key, drop_level=drop_level)
3736         else:
-> 3737             loc = self.index.get_loc(key)
3738
3739             if isinstance(loc, np.ndarray):

/usr/local/lib/python3.6/dist-packages/pandas/core/indexes/base.py in get_loc(self, key, method, tolerance)
2897                 return self._engine.get_loc(key)
2898             except KeyError:
-> 2899                 return self._engine.get_loc(self._maybe_cast_indexer(key))
2900         indexer = self.get_indexer([key], method=method, tolerance=tolerance)
2901         if indexer.ndim > 1 or indexer.size > 1:

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc()

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

KeyError: 'Mn'

Pulp is not giving a very well-formed error message here (not sure if Pulp can actually do that -- Python is issuing this before Pulp sees what is happening). But at least we are alerted (rather heavy-handedly) there is something wrong when we use Mn.  Careful inspection of the stack frame shows we have a problem in the constraint

model += demand*content[j] == lp.lpSum([use[i]*supplyData.loc[i,j] for i in Items])

This error is actually generating an exception inside an exception handler!

Of course a much better and simpler error message would be: "supplyData.loc["A","Mn"]: Column "Mn" not found in data frame supplyData." IMHO programmers do not pay enough attention to provide meaningful error messages.

#### Solver log

The default solver is CBC via a DLL call. I don't think it is possible to see the solver log using this setup. I prefer to see the solver log, just to make sure there are no surprises. For this I used:

model.solve(lp.COIN_CMD(msg=1))

This will call the CBC executable (via an MPS file) and it will show the CBC log:

Welcome to the CBC MILP Solver
Version: 2.9
Build Date: Jan  6 2019

command line - cbc.exe fbbd73baa2494109b8e990ce26eb79b6-pulp.mps branch printingOptions all solution fbbd73baa2494109b8e990ce26eb79b6-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 22 COLUMNS
At line 64 RHS
At line 82 BOUNDS
At line 83 ENDATA
Problem MODEL has 17 rows, 10 columns and 34 elements
Coin0008I MODEL read with 0 errors
Presolve 4 (-13) rows, 7 (-3) columns and 18 (-16) elements
0  Obj 479.87991 Primal inf 10199.217 (4)
3  Obj 5887.5743
Optimal - objective value 5887.5743
After Postsolve, objective 5887.5743, infeasibilities - dual 1275.5592 (2), primal 0 (0)
Presolved model was optimal, full model needs cleaning up
Optimal - objective value 5887.5743
Optimal objective 5887.574275 - 3 iterations time 0.012, Presolve 0.00
Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.03   (Wallclock seconds):       0.03

The solver log shows that the presolver removes 13 of the 17 rows. This high reduction rate is related to the singleton constraints. When we look at the model we generated 3+3+7=13 bound constraints. The presolver is getting rid of these and makes proper bounds of them.

Note that: msg=1 may not always work when running as a Jupyter notebook.

#### Debugging

For debugging PuLP models, I recommend:

1. print(model). Printing the model shows how PuLP interpreted the constraints.
2. Writing an LP file: model.writeLP("steel.lp").

The output of print(model) is:

Steel:
MINIMIZE
1.2*Use_A + 1.5*Use_B + 0.9*Use_C + 1.3*Use_D + 1.45*Use_E + 1.2*Use_F + 1.0*Use_G + 0.0
SUBJECT TO
_C1: Use_A <= 4000

_C2: Use_B <= 3000

_C3: Use_C <= 6000

_C4: Use_D <= 5000

_C5: Use_E <= 2000

_C6: Use_F <= 3000

_C7: Use_G <= 2500

_C8: 5000 Content_C - 2.5 Use_A - 3 Use_B = 0

_C9: Content_C >= 2

_C10: Content_C <= 3

_C11: 5000 Content_Cu - 0.3 Use_C - 90 Use_D - 96 Use_E - 0.4 Use_F
- 0.6 Use_G = 0

_C12: Content_Cu >= 0.4

_C13: Content_Cu <= 0.6

_C14: 5000 Content_Mn - 1.3 Use_A - 0.8 Use_B - 4 Use_E - 1.2 Use_F = 0

_C15: Content_Mn >= 1.2

_C16: Content_Mn <= 1.65

_C17: Use_A + Use_B + Use_C + Use_D + Use_E + Use_F + Use_G = 5000

VARIABLES
Content_C Continuous
Content_Cu Continuous
Content_Mn Continuous
Use_A Continuous
Use_B Continuous
Use_C Continuous
Use_D Continuous
Use_E Continuous
Use_F Continuous
Use_G Continuous

The LP file looks like:

\* Steel *\
Minimize
OBJ: 1.2 Use_A + 1.5 Use_B + 0.9 Use_C + 1.3 Use_D + 1.45 Use_E + 1.2 Use_F
+ Use_G
Subject To
_C1: Use_A <= 4000
_C10: Content_C <= 3
_C11: 5000 Content_Cu - 0.3 Use_C - 90 Use_D - 96 Use_E - 0.4 Use_F
- 0.6 Use_G = 0
_C12: Content_Cu >= 0.4
_C13: Content_Cu <= 0.6
_C14: 5000 Content_Mn - 1.3 Use_A - 0.8 Use_B - 4 Use_E - 1.2 Use_F = 0
_C15: Content_Mn >= 1.2
_C16: Content_Mn <= 1.65
_C17: Use_A + Use_B + Use_C + Use_D + Use_E + Use_F + Use_G = 5000
_C2: Use_B <= 3000
_C3: Use_C <= 6000
_C4: Use_D <= 5000
_C5: Use_E <= 2000
_C6: Use_F <= 3000
_C7: Use_G <= 2500
_C8: 5000 Content_C - 2.5 Use_A - 3 Use_B = 0
_C9: Content_C >= 2
End

The information is basically the same, but the ordering of the rows is a bit different.

#### Comparison to CVXPY

We can try to model and solve the same problem using CVXPY. CVXPY is matrix oriented, so very different than PuLP. Here is my attempt:

from io import StringIO
import pandas as pd
import numpy as np
import cvxpy as cp

# for inputting tabular data below
def table(s):

#------------------------------------------------------------------
# data
#------------------------------------------------------------------

demand = 5000

requirements = table("""
ID  Element      Min   Max
C   Carbon       2     3
Cu  Copper       0.4   0.6
Mn  Manganese    1.2   1.65
""")

supplyData = table("""
ID  Alloy             C       Cu     Mn     Stock   Price
A   "Iron alloy"      2.50    0.00   1.30   4000    1.20
B   "Iron alloy"      3.00    0.00   0.80   3000    1.50
C   "Iron alloy"      0.00    0.30   0.00   6000    0.90
D   "Copper alloy"    0.00   90.00   0.00   5000    1.30
E   "Copper alloy"    0.00   96.00   4.00   2000    1.45
F   "Aluminum alloy"  0.00    0.40   1.20   3000    1.20
G   "Aluminum alloy"  0.00    0.60   0.00   2500    1.00
""")

print("----- Data-------")
print(requirements)
print(supplyData)

#------------------------------------------------------------------
# derived data
#------------------------------------------------------------------

# our sets are stockItems ["A","B",..] and elements ["C","Cu",...]
Items = supplyData.index
Elements = requirements.index

# extract arrays (make sure order is identical)
Min = requirements.loc[Elements,"Min"]
Max = requirements.loc[Elements,"Max"]
Cost = supplyData.loc[Items,"Price"]
Avail = supplyData.loc[Items,"Stock"]
Element = supplyData.loc[Items,Elements]

# counts
NumItems = np.shape(Items)[0]
NumElements  = np.shape(Elements)[0]

# reshape into proper Numpy column vectors to make cvxpy happy
Min = np.reshape(Min.to_numpy(),(NumElements,1))
Max = np.reshape(Max.to_numpy(),(NumElements,1))
Cost = np.reshape(Cost.to_numpy(),(NumItems,1))
Avail = np.reshape(Avail.to_numpy(),(NumItems,1))
Element = Element.to_numpy()

#------------------------------------------------------------------
# LP Model
#------------------------------------------------------------------

use = cp.Variable((NumItems,1),"Use",nonneg=True)
content = cp.Variable((NumElements,1),"Content",nonneg=True)

model = cp.Problem(cp.Minimize(Cost.T @ use),
[cp.sum(use) == demand,
cp.multiply(demand,content) == Element.T @ use,
content >= Min,
content <= Max,
use <= Avail
])

#------------------------------------------------------------------
# Solve and reporting
#------------------------------------------------------------------

model.solve(solver=cp.ECOS,verbose=True)

print("----- Model Results-------")
print("status:",model.status)
print("objective:",model.value)
results = pd.DataFrame({'variable':'use',
'index': Items,
'lower':0,
'level':use.value.flatten(),
'upper':Avail.flatten()
})
results = results.append(pd.DataFrame({'variable':'content',
'index': Elements,
'lower':Min.flatten(),
'level':content.value.flatten(),
'upper':Max.flatten()
}))
print(results)

Notes:

• I did my best to make sure that the ordering of rows and columns in the data frames is not significant.
• We convert the information in the data frames to standard NumPy arrays for the benefit of CVXPY. (A column in a dataframe is a pandas series).
• If we don't do proper shaping of the arrays, we may see error messages like: ValueError: Cannot broadcast dimensions  (3,) (3, 1)
• The model is compact, but we needed to put more effort in data extraction. In optimization, it is not at all unusual that data stuff takes more effort than the model equations.

The results look like:

----- Data-------
Element  Min   Max
ID
C      Carbon  2.0  3.00
Cu     Copper  0.4  0.60
Mn  Manganese  1.2  1.65
Alloy    C    Cu   Mn  Stock  Price
ID
A       Iron alloy  2.5   0.0  1.3   4000   1.20
B       Iron alloy  3.0   0.0  0.8   3000   1.50
C       Iron alloy  0.0   0.3  0.0   6000   0.90
D     Copper alloy  0.0  90.0  0.0   5000   1.30
E     Copper alloy  0.0  96.0  4.0   2000   1.45
F   Aluminum alloy  0.0   0.4  1.2   3000   1.20
G   Aluminum alloy  0.0   0.6  0.0   2500   1.00

ECOS 2.0.7 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS

It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
0  +6.293e+03  -4.362e+04  +9e+04  1e-01  8e-02  1e+00  4e+03    ---    ---    1  1  - |  -  -
1  +5.462e+03  -6.110e+04  +7e+04  2e-01  5e-02  2e+03  3e+03  0.5361  8e-01   0  0  0 |  0  0
2  +5.497e+03  +1.418e+03  +7e+03  1e-02  4e-03  5e+02  3e+02  0.9313  3e-02   0  0  0 |  0  0
3  +4.981e+03  +3.654e+03  +2e+03  4e-03  1e-03  2e+02  1e+02  0.6947  5e-02   0  0  0 |  0  0
4  +5.687e+03  +3.974e+03  +2e+03  7e-03  9e-04  3e+02  9e+01  0.4022  7e-01   0  0  0 |  0  0
5  +5.653e+03  +5.326e+03  +5e+02  1e-03  2e-04  2e+01  2e+01  0.9127  1e-01   0  0  0 |  0  0
6  +5.692e+03  +5.535e+03  +2e+02  5e-04  8e-05  1e+01  1e+01  0.5874  1e-01   0  0  0 |  0  0
7  +5.791e+03  +5.642e+03  +2e+02  6e-04  5e-05  2e+01  7e+00  0.7361  5e-01   0  0  0 |  0  0
8  +5.843e+03  +5.798e+03  +6e+01  2e-04  2e-05  5e+00  2e+00  0.9890  4e-01   0  0  0 |  0  0
9  +5.886e+03  +5.883e+03  +4e+00  1e-05  1e-06  3e-01  2e-01  0.9454  1e-02   0  0  0 |  0  0
10  +5.888e+03  +5.888e+03  +5e-02  2e-07  3e-08  4e-03  2e-03  0.9890  2e-03   0  0  0 |  0  0
11  +5.888e+03  +5.888e+03  +6e-04  2e-09  5e-10  5e-05  2e-05  0.9890  1e-04   1  0  0 |  0  0
12  +5.888e+03  +5.888e+03  +6e-06  2e-11  6e-12  5e-07  3e-07  0.9890  1e-04   1  0  0 |  0  0

OPTIMAL (within feastol=1.9e-11, reltol=1.1e-09, abstol=6.3e-06).
Runtime: 0.000566 seconds.

----- Model Results-------
status: optimal
objective: 5887.574272281105
variable index  lower         level    upper
0      use     A    0.0  4.000000e+03  4000.00
1      use     B    0.0  1.283254e-06  3000.00
2      use     C    0.0  3.977630e+02  6000.00
3      use     D    0.0  5.135476e-07  5000.00
4      use     E    0.0  2.761272e+01  2000.00
5      use     F    0.0  5.746243e+02  3000.00
6      use     G    0.0  5.163966e-06  2500.00
0  content     C    2.0  2.000000e+00     3.00
1  content    Cu    0.4  5.999999e-01     0.60
2  content    Mn    1.2  1.200000e+00     1.65

The results are not rounded. That often means that the solution of an interior point algorithm looks a bit ugly. In essence this is the same solution as we found with PuLP/CBC.

#### Comparison to GAMS

The GAMS model is closer to PuLP:

*---------------------------------------------------
* data
*---------------------------------------------------

set
i 'items from inventory' /A*G/
j 'elements' /C,Cu,Mn/
;

table requirements(j,*)
Min   Max
C     2     3
Cu    0.4   0.6
Mn    1.2   1.65
;

table supplyData(i,*)
C      Cu     Mn     Stock   Price
A    2.5            1.3    4000    1.20
B    3              0.8    3000    1.50
C            0.3           6000    0.90
D           90             5000    1.30
E           96      4      2000    1.45
F            0.40   1.20   3000    1.20
G            0.60          2500    1.00
;

scalar demand /5000/;

*---------------------------------------------------
* LP Model
*---------------------------------------------------

positive variables
use(i)       'usage of raw material'
content(j)   'characteristics of final product'
;

* bounds
use.up(i) = supplyData(i,'Stock');
content.lo(j) = requirements(j,'Min');
content.up(j) = requirements(j,'Max');

variable totalCost 'objective variable';

equations
obj  'objective'
calcContent(j) 'calculate contents of final product'
meetDemand 'total demand must be met'
;

obj.. totalCost =e= sum(i, use(i)*supplyData(i,'Price'));

calcContent(j).. demand*content(j) =e= sum(i, use(i)*supplyData(i,j));

meetDemand.. sum(i, use(i)) =e= demand;

model blending /all/;

*---------------------------------------------------
* Solve and reporting
*---------------------------------------------------

solve blending minimizing totalCost using lp;

parameter results(*,*,*);
results('use',i,'lower') = 0;
results('use',i,'level') = use.l(i);
results('use',i,'upper') = supplyData(i,'Stock');
results('content',j,'lower') = requirements(j,'Min');
results('content',j,'level') = content.l(j);
results('content',j,'upper') = requirements(j,'Max');
display results;

Notes:

• In GAMS we start with the sets, followed by the data. With external data we can extract sets from data.
• The rows and columns in the tables are not order dependent: we can change the order without changing the meaning of the model.
• We can specify bounds directly: no need for singleton constraints.
• The statement results('use',i,'lower') = 0; is not really needed as the default value is zero. (To be more precise: all data is stored sparse, where "does not exist" is the same as zero). However, this statement forces the set element 'lower' to be before 'level' and 'upper'. Another way to enforce this, is to introduce a dummy set: set dummy /lower,level,upper/. We don't have to use the set. Just the declaration would instill an element ordering.
• I used * in the data tables. This is conveniently disabling domain checking, but it is also dangerous. In this case, if we change in the supplyData table "Mn" to "Mng", we will not see an error. It would be safer to use full domain checking, which requires to add a few sets.
• The safe version of the data part in GAMS would be:
set
i    'items from inventory' /A*G/
d    'demand limits' /Min,Max/
s    'supply data' /C,Cu,Mn,Stock,Price/
j(s) 'elements' /C,Cu,Mn/
;

table requirements(j,d)
Min   Max
C     2     3
Cu    0.4   0.6
Mn    1.2   1.65
;

table supplyData(i,s)
C      Cu     Mng    Stock   Price
A    2.5            1.3    4000    1.20
B    3              0.8    3000    1.50
C            0.3           6000    0.90
D           90             5000    1.30
E           96      4      2000    1.45
F            0.40   1.20   3000    1.20
G            0.60          2500    1.00
;

Here the typo "Mng" will immediately give a good error message:
  20  table supplyData(i,s)
21         C      Cu     Mng    Stock   Price
****                         \$170
**** 170  Domain violation for element
For production models it is best to use domain checking throughout.
• IMHO, this GAMS model looks a bit more streamlined than the Python models.

The output looks like:

----     70 PARAMETER results

lower       level       upper

use    .A                 4000.000    4000.000
use    .B                             3000.000
use    .C                  397.763    6000.000
use    .D                             5000.000
use    .E                   27.613    2000.000
use    .F                  574.624    3000.000
use    .G                             2500.000
content.C        2.000       2.000       3.000
content.Cu       0.400       0.600       0.600
content.Mn       1.200       1.200       1.650

#### References

1. Translating a LP from Excel to Python, https://stackoverflow.com/questions/59579342/translating-a-lp-from-excel-to-python-pulp

#### 1 comment:

1. > Min = np.reshape(Min.to_numpy(),(NumElements,1))

You can do this via Min = Min.values.reshape(-1,1)


text:
# Congruence with variable modulus

I have been working on some problems and one of them has been particularly challenging. The problem is as follows.

Find a non-trivial (meaning more than 1 digit) positive integer a that satisfies:

$a/10 + 3121(a \mod 10) = 0 \mod a$

Here the division operator is meant to denote integer division so algebraically this is equivalent to:

$(\lfloor a/10\rfloor + 3121(a - 10*\lfloor a/10\rfloor)) - a*\lfloor\lfloor a/10\rfloor + 3121(a - 10*\lfloor a/10\rfloor))/a)\rfloor = 0$

I believe there should be some straightforward way to calculate this value. Using a brute force guess and check I found that $a = 101$ satisfied this congruence but I wanted to know if there was a way to analytically work this one out?

• Is the sort of topic covered in a normal number theory book? – frogeyedpeas Jan 14 '13 at 23:01
• If so, what book should I read to learn how to solve this type of problem – frogeyedpeas Jan 14 '13 at 23:01
• Nobody? :( (*thanks Jon for the formatting) – frogeyedpeas Jan 15 '13 at 1:12
• I think your question, as it is now, is pretty confusing: what does "$\,3121(a\pmod {10})\,$" mean? . Then what you say your original expression is equivalent too is waaaay more confusing (at least for me). In fact, I've no idea what is this related to and how did this beast come into life... – DonAntonio Jan 15 '13 at 3:25
• well the 3121 is multiplied to a mod 10 – frogeyedpeas Jan 15 '13 at 3:38

What about $\,a=3\,$? According to what you wrote, and I'm not saying I fully understand how, where, when...of it, we get:

$$\left[\frac{3}{10}\right]+3,121\cdot 3=0+3\cdot 3,121=0\pmod 3$$

And the above , of course, works the same for any positive integer $\,a<10\,$ ...

Added: If you want $\,a\,$ with more than one digit take $\,a=103\,$ :

$$\left[\frac{103}{10}\right]+3,121\cdot (103\pmod{10})=10+3\cdot 3,121=9,373=103\cdot 91$$

• yea, those were what i meant by trivial answers (1 digit) – frogeyedpeas Jan 16 '13 at 0:04

\begin{align} \left \lfloor \dfrac{a}{10} \right \rfloor + 3121 \left (a - 10 \left \lfloor \dfrac{a}{10} \right \rfloor \right ) &= 0 \pmod a \\ \left \lfloor \dfrac{a}{10} \right \rfloor - 31210 \left \lfloor \dfrac{a}{10} \right \rfloor &= 0 \pmod a \\ 31209 \left \lfloor \dfrac{a}{10} \right \rfloor &= 0 \pmod a \\ a \; &\bigg| \; 31209 \left \lfloor \dfrac{a}{10} \right \rfloor \end{align}

Let $a = 10x + y,$ where $1 \le y \le 9,$ then $(10x + y) \; \bigg| \; 31209x$

The divisors of $31209,$ sorted by their units digit, are

   1    101
3    103  303  10403
309  31209

Let $g = \gcd \left(a, \left \lfloor \dfrac{a}{10} \right \rfloor \right) = \gcd(10x+y, x) = \gcd(x, y).$ It follows that the natural numbers $\dfrac{10x+y}{g}$ and $\dfrac x g$ are relatively prime to each other.

\begin{align} (10x + y) \; &\bigg| \; 31209x \cr \left (\dfrac{10x+y}{g}\right ) \; &\bigg| \; 31209 \dfrac x g \cr \left (\dfrac{10x+y}{g}\right ) \; &\bigg| \; 31209 \cr \left( 10\dfrac x g + \dfrac y g \right) \; &\bigg| \; 31209 \cr \end{align}

Define $\Delta(u)$ to be the set of those divisors of 31209 which have a units digit of $u$.

Note that only $\Delta(1)$, $\Delta(3)$, and $\Delta(9)$ are the only non-empty sets.

So, if $g$ is a common divisor of $a = 10x + y$ and $\left \lfloor \dfrac{a}{10} \right \rfloor = x$, then $\dfrac a g = 10\dfrac x g + \dfrac y g$ is a divisor of $31209.$ It follows that the solutions to $a \; \bigg| \; 31209 \left \lfloor \dfrac{a}{10} \right \rfloor$ that correspond to $g$ are

$$a \in g \; \Delta \left( \dfrac y g \right)$$

We can run through the $10$ possible values of y pretty quicky now.

## $y = 0$

$10 \mid 31209 \implies a \in \varnothing$

## $y = 1$

$g = 1$: $\; (10x + 1) \mid 31209 \implies a \in \Delta(1) = \{1, 101\}$

## $y = 2$

$g = 1$: $\; (10x + 2) \mid 31209 \implies a \in \Delta(2) = \varnothing$

$g = 2$: $\; (10(x/2) + 1) \mid 31209 \implies a \in 2\Delta(1) = \{2, 202 \}$

## $y = 3$

$g = 1$: $\; (10x + 3) \mid 31209 \implies a \in \Delta(3) = \{3, 103, 303, 10403 \}$

$g = 3$: $\; (10(x/3) + 1) \mid 31209 \implies a \in 3\Delta(1) = \{3, 303 \}$

## $y = 4$

$g = 1$: $\; (10x + 4) \mid 31209 \implies a \in \Delta(4) = \varnothing$

$g = 2$: $\; (10(x/2) + 2) \mid 31209 \implies a \in 2\Delta(2) = \varnothing$

$g = 4$: $\; (10(x/4) + 1) \mid 31209 \implies a \in 4\Delta(1) = \{4, 404 \}$

## $y = 5$

$g = 1$: $\; (10x + 5) \mid 31209 \implies a \in \Delta(5) = \varnothing$

$g = 5$: $\; (10(x/5) + 1) \mid 31209 \implies a \in 5\Delta(1) = \{5, 505 \}$

## $y = 6$

$g = 1$: $\; (10x + 6) \mid 31209 \implies a \in \Delta(6) = \varnothing$

$g = 2$: $\; (10(x/2) + 3) \mid 31209 \implies a \in 2\Delta(3) = \{ 6, 206, 606, 20806 \}$

$g = 3$: $\; (10(x/3) + 2) \mid 31209 \implies a \in 4\Delta(2) = \varnothing$

$g = 6$: $\; (10(x/6) + 1) \mid 31209 \implies a \in 6\Delta(1) = \{ 6, 606 \}$

## $y = 7$

$g = 1$: $\; (10x + 7) \mid 31209 \implies a \in \Delta(7) = \varnothing$

$g = 7$: $\; (10(x/7) + 1) \mid 31209 \implies a \in 7\Delta(1) = \{7, 707 \}$

## $y = 8$

$g = 1$: $\; (10x + 8) \mid 31209 \implies a \in \Delta(8) = \varnothing$

$g = 2$: $\; (10(x/2) + 4) \mid 31209 \implies a \in 2\Delta(4) = \varnothing$

$g = 4$: $\; (10(x/4) + 2) \mid 31209 \implies a \in 4\Delta(2) = \varnothing$

$g = 8$: $\; (10(x/8) + 1) \mid 31209 \implies a \in 8\Delta(1) = \{8, 808 \}$

## $y = 9$

$g = 1$: $\; (10x + 9) \mid 31209 \implies a \in \Delta(9) = \{309, 31209 \}$

$g = 3$: $\; (10(x/3) + 3) \mid 31209 \implies a \in 3\Delta(3) = \{9, 309, 909, 31209 \}$

$g = 9$: $\; (10(x/9) + 1) \mid 31209 \implies a \in 9\Delta(1) = \{9, 909 \}$

## Solution Set

    1    2      3     4     5     6     7     8     9
103               206               309
101   202   303   404   505   606   707   808   909
10403             20806             31209

format: plain


text:
# Extracting expression data in r from GEO and got a S4 object

I want to open the expression data GSE9838 dataset from GEO using GEOquery R package.

gset <- getGEO(filename="GSE9838_family.soft.gz")

The result of this code said it's S4 object of class GSE.

How can I get the expression data of GSE9838?

• Welcome to the site. Have you looked for the tutorial of GEOquery? You can see them with vignette(package="GEOquery") . Also have you looked around on this site and in support.bioconductor.org for how have other people done this?
– llrs
Mar 28 '19 at 12:42
• Possible duplicate of Extracting expression data from GSE dataset downloaded from GEO
– llrs
Mar 28 '19 at 12:43
• Not sure if that method works in this case
– Oka
Mar 28 '19 at 19:15
• @Oka Why shouldn't it work? Have Qiran tested it? At least mentioning that it didn't work it would help the people trying to answer the question (and future readers)
– llrs
Mar 29 '19 at 9:38
• @llrs I´ve checked if it would, and it didn´t (at least on the first trial). The method worked fine with GSE used in the example, but didn´t seem to work in this case for some reason..
– Oka
Mar 29 '19 at 13:00

Well, if you use a SOFT formatted files:

#You can read your file
gds <- getGEO(filename="GSE9838_family.soft.gz")

# and check it
class(gds)
show(gds)

#then you usually would try to use  an exprs() accessor to retrieve the expression matrix:
exprs (gds[[1]])

# but in this case you will get an error revealing that it is a GSE datastructure
# eg. is composite structure having GSM and GPL object lists:
names(GSMList(gds))
names(GPLList(gds))

# now GPLList is probably not interesting in this case
# but you can look at the GSM-objects from GSMList:
first <- GSMList(gds)[[1]]
first
class(first)
Meta(first)

# at this point if you read the results you´ll notice that
#    title
#    [1] "X_intercept"

# and
#   An object of class "GEODataTable"
#   ****** Column Descriptions ******
#   [1] Column      Description
#   <0 rows> (or 0-length row.names)
#   ****** Data Table ******
#   data frame with 0 columns and 0 rows

So here you would probably wonder if the data you want is these objects. If you go back go GEO page for GSE9838, you can find that "Raw data provided as supplementary file. Processed data is available on Series record". So one may want to check series record instead:

# you can query the GSE9838 like this:
data = getGEO("GSE9838")

# then you can check what you´ve got
names(data)
#[1] "GSE9838_series_matrix.txt.gz"
class(data$GSE9838_series_matrix.txt.gz) # and check what´s inside: exprs(data$GSE9838_series_matrix.txt.gz)
# but apparently it just shows one header row.

So it wasn´t there either. And if you check the size of the file, you´ll find that it is just 4.8 kb. My educated guess would be that the expression values are not in the dataset itself but probably in some of the supplementary files?

filePaths = getGEOSuppFiles("GSE9838")


text:
# Find the 6

Find the total cost of 10 computers at $2100 each and 7 boxes of diskettes at$12 each

Correct result:

c =  21084

#### Solution:

We would be pleased if you find an error in the word problem, spelling mistakes, or inaccuracies and send it to us. Thank you!

## Next similar math problems:

• Startup
Jaxon’s start up business makes a profit of $450 during the first month. However, the company records a profit of -$60 per month for the next four months and a profit of $125 for the final month. What is the total profit for the first six months of Jaxon’ • Theater The theater has in each row with 19 seats. Ticket to the first 10 rows is for 30 USD. In the next rows is for 15 USD. The theater was completely sold out. Revenue was 12255 USD. How many rows are in the theater? • Loan If you take a bank loan$ 10000 and we want to repay after the year, we have to pay the total amount $10320/ What is the annual interest rate on this loan? • An investor An investor purchased 75 shares of stock at$42.80 per share, and 6 months later sold the shares for $47.15 per share. The$20 stockbroker fee will not be charged for this sale. Calculate the total profit or loss of the transaction.
• Telco company
The upstairs communications company offers customers a special long distance calling rate that includes a $0.10 per minute charge. Which of the following represents this fee scheduale where m represents the number of minutes and c is the overall cost of t • Exchange rates If the Canadian dollar appreciated by C$0.005 relative to the US dollar, what would be the new value of the Canadian dollar per US dollar? Assume the current exchange rate was US$1 = C$0.907.
• Balance of account
Theo had a balance of -$4 in his savings account. After making a deposit, he has$25 in his account. What is the overall change to his account?
• Saving 9
An amount if $2000 is invested at an interest of 5% per month. if$ 200 is added at the beginning of each successive month but no withdrawals. Give an expression for the value accumulated after n months. After how many months will the amount have accumul
• The Stolen Money
A man walks into a store and steals a $100 bill. 5 minutes later, he returns to the store and buys stuff worth$70. He pays with the bill that he had stolen, so the owner of the store returns him $30. How many dollars did the store owner lose? • Milk bill Mrs Tara buys 2 liters of milk daily. If 1 liter of milk cost$0.27. What will be her milk bill for 30 days.
• Nine books
Nine books are to be bought by a student. Art books cost $6.00 each and biology books cost$6.50 each . If the total amount spent was $56.00, how many of each book was bought? • Conference 2 9 people attended a conference on behalf of their company. The conference fee was £520 per person and the company paid a total of £856 in travel costs for the 9 people. How much did the conference cost the company altogether? • College 2 College student is moving into a dormitory. The student rent a truck for$19.95 plus $0.99 per mile. Before returning the truck the student fills the tank with gasoline, which cost$65.32. Total cost $144.67. Using a linear equation, explain the process t • Andy had Andy had$5.95 in his pocket. His father gave him another $3.55. How much does Andy have altogether now? • Tennis balls Can of tennis balls contains 3 balls per can and cost$7 how much will it cost for 36 tennis balls?
• Stacks
Annie has a total of \$ 702. The money must be divided into stacks so that each buyer has the same amounth of dollars. How many options she have?
• Juice box
In the box is 0.3 liters of juice. How many liters of juice contains 3 these boxes?


text: 5 added 5 characters in body edited Feb 26 '14 at 9:15 Rodriguez 15622 silver badges77 bronze badges Yes. Though the energy will not be unbounded, but bounded from above, if my calculation is correct. For real scalar field under $$(+---)$$ metric, besides the negative classical kinetic energy for the Lagrangian $$\mathcal{L}=-\frac{1}{2} \partial^{\mu} \phi \partial_{\mu} \phi - \frac{1}{2} m^2 \phi^2 \tag{1}$$, the classical equation of motion will be $$(\square - m^2 )\phi=0 . \tag{2}$$ For plane wave $$\phi ~\sim e^{ipx}$$, it gives $$p^2+m^2 = (p^0)^2 - \mathbf{p}^2+m^2=0$$ which is inconsistent with relativistic energy momentum relation. I am not sure if it is necessary to quantize it. Though the energy-momentum-relation argument will not work for the Dirac field, we can quantize it to see the energy will be negative definite. $$\mathcal{L} = \bar{\psi}( -i \gamma^{\mu} \partial_{\mu} - m ) \psi \tag{3}$$ The classical equation of motion is $$(i \gamma^{\mu} \partial_{\mu} +m) \psi=0 \tag{4}$$ To preserve all properties of $$u(p)$$$$u(\mathbf{p})$$ and $$v(p)$$$$v(\mathbf{p})$$, we define $$\psi =: u(\mathbf{p}) e^{ipx}$$ $$\psi =: v(\mathbf{p}) e^{-ipx}$$ Thus we can replace the $$u(\mathbf{p})$$ as $$v(\mathbf{p})$$ and $$v(\mathbf{p})$$ as $$u(\mathbf{p})$$ in the expansion of $$\psi$$ and $$\bar{\psi}$$. By $$\pi = -i \bar{\psi} \gamma^0$$ then $$H= \int d^3 x \bar{\psi} ( i \gamma^i \partial_i \psi + m ) \psi$$ Plug in expansionexpansions of spinors in the Schrodinger picture $$\psi = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( a_{\mathbf{p}}^s v^s (\mathbf{p}) e^{-i\mathbf{p} \cdot \mathbf{x} } + b_{\mathbf{p}}^{s\dagger} u^s(\mathbf{p}) e^{i \mathbf{p} \cdot \mathbf{x} } \right)$$ $$\bar{\psi} = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( b_{\mathbf{p}}^s \bar{u}^s (\mathbf{p}) e^{-i\mathbf{p} \cdot \mathbf{x}} + a_{\mathbf{p}}^{s\dagger} \bar{v}^s(\mathbf{p}) e^{i\mathbf{p} \cdot \mathbf{x}} \right)$$ we have $$H = \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( - \gamma^i p_i +m) u^s(\mathbf{p}) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( \gamma^i p_i +m) v^s(\mathbf{p})$$$$H = \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( - \gamma^i p_i +m) u^s(\mathbf{p}) + a_{\mathbf{p}}^{s'\dagger} a_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( \gamma^i p_i +m) v^s(\mathbf{p})$$ $$= \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( \gamma^0 p_0 ) u^s(\mathbf{p}) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( - \gamma^0 p_0 ) v^s(\mathbf{p})$$$$= \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( \gamma^0 p_0 ) u^s(\mathbf{p}) + a_{\mathbf{p}}^{s'\dagger} a_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( - \gamma^0 p_0 ) v^s(\mathbf{p})$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} E_{\mathbf{p}} ( b_{\mathbf{p}}^{s} b_{\mathbf{p}}^{s\dagger} - a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} )$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} - E_{\mathbf{p}} (b_{\mathbf{p}}^{s\dagger} b_{\mathbf{p}}^{s} + a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} ) - \infty$$ Changing anticommutator into commutator will make the spectrum unbounded. Yes. Though the energy will not be unbounded, but bounded from above, if my calculation is correct. For real scalar field under $$(+---)$$ metric, besides the negative classical kinetic energy for the Lagrangian $$\mathcal{L}=-\frac{1}{2} \partial^{\mu} \phi \partial_{\mu} \phi - \frac{1}{2} m^2 \phi^2 \tag{1}$$, the classical equation of motion will be $$(\square - m^2 )\phi=0 . \tag{2}$$ For plane wave $$\phi ~\sim e^{ipx}$$, it gives $$p^2+m^2 = (p^0)^2 - \mathbf{p}^2+m^2=0$$ which is inconsistent with relativistic energy momentum relation. I am not sure if it is necessary to quantize it. Though the energy-momentum-relation argument will not work for the Dirac field, we can quantize it to see the energy will be negative definite. $$\mathcal{L} = \bar{\psi}( -i \gamma^{\mu} \partial_{\mu} - m ) \psi \tag{3}$$ The classical equation of motion is $$(i \gamma^{\mu} \partial_{\mu} +m) \psi=0 \tag{4}$$ To preserve all properties of $$u(p)$$ and $$v(p)$$, we define $$\psi =: u(\mathbf{p}) e^{ipx}$$ $$\psi =: v(\mathbf{p}) e^{-ipx}$$ Thus we can replace the $$u(\mathbf{p})$$ as $$v(\mathbf{p})$$ and $$v(\mathbf{p})$$ as $$u(\mathbf{p})$$ in the expansion of $$\psi$$ and $$\bar{\psi}$$. By $$\pi = -i \bar{\psi} \gamma^0$$ then $$H= \int d^3 x \bar{\psi} ( i \gamma^i \partial_i \psi + m ) \psi$$ Plug in expansion of spinors in Schrodinger picture $$\psi = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( a_{\mathbf{p}}^s v^s (\mathbf{p}) e^{-i\mathbf{p} \cdot \mathbf{x} } + b_{\mathbf{p}}^{s\dagger} u^s(\mathbf{p}) e^{i \mathbf{p} \cdot \mathbf{x} } \right)$$ $$\bar{\psi} = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( b_{\mathbf{p}}^s \bar{u}^s (\mathbf{p}) e^{-i\mathbf{p} \cdot \mathbf{x}} + a_{\mathbf{p}}^{s\dagger} \bar{v}^s(\mathbf{p}) e^{i\mathbf{p} \cdot \mathbf{x}} \right)$$ we have $$H = \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( - \gamma^i p_i +m) u^s(\mathbf{p}) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( \gamma^i p_i +m) v^s(\mathbf{p})$$ $$= \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( \gamma^0 p_0 ) u^s(\mathbf{p}) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( - \gamma^0 p_0 ) v^s(\mathbf{p})$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} E_{\mathbf{p}} ( b_{\mathbf{p}}^{s} b_{\mathbf{p}}^{s\dagger} - a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} )$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} - E_{\mathbf{p}} (b_{\mathbf{p}}^{s\dagger} b_{\mathbf{p}}^{s} + a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} ) - \infty$$ Changing anticommutator into commutator will make the spectrum unbounded. Yes. Though the energy will not be unbounded, but bounded from above, if my calculation is correct. For real scalar field under $$(+---)$$ metric, besides the negative classical kinetic energy for the Lagrangian $$\mathcal{L}=-\frac{1}{2} \partial^{\mu} \phi \partial_{\mu} \phi - \frac{1}{2} m^2 \phi^2 \tag{1}$$, the classical equation of motion will be $$(\square - m^2 )\phi=0 . \tag{2}$$ For plane wave $$\phi ~\sim e^{ipx}$$, it gives $$p^2+m^2 = (p^0)^2 - \mathbf{p}^2+m^2=0$$ which is inconsistent with relativistic energy momentum relation. I am not sure if it is necessary to quantize it. Though the energy-momentum-relation argument will not work for the Dirac field, we can quantize it to see the energy will be negative definite. $$\mathcal{L} = \bar{\psi}( -i \gamma^{\mu} \partial_{\mu} - m ) \psi \tag{3}$$ The classical equation of motion is $$(i \gamma^{\mu} \partial_{\mu} +m) \psi=0 \tag{4}$$ To preserve all properties of $$u(\mathbf{p})$$ and $$v(\mathbf{p})$$, we define $$\psi =: u(\mathbf{p}) e^{ipx}$$ $$\psi =: v(\mathbf{p}) e^{-ipx}$$ Thus we can replace the $$u(\mathbf{p})$$ as $$v(\mathbf{p})$$ and $$v(\mathbf{p})$$ as $$u(\mathbf{p})$$ in the expansion of $$\psi$$ and $$\bar{\psi}$$. By $$\pi = -i \bar{\psi} \gamma^0$$ then $$H= \int d^3 x \bar{\psi} ( i \gamma^i \partial_i \psi + m ) \psi$$ Plug in expansions of spinors in the Schrodinger picture $$\psi = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( a_{\mathbf{p}}^s v^s (\mathbf{p}) e^{-i\mathbf{p} \cdot \mathbf{x} } + b_{\mathbf{p}}^{s\dagger} u^s(\mathbf{p}) e^{i \mathbf{p} \cdot \mathbf{x} } \right)$$ $$\bar{\psi} = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( b_{\mathbf{p}}^s \bar{u}^s (\mathbf{p}) e^{-i\mathbf{p} \cdot \mathbf{x}} + a_{\mathbf{p}}^{s\dagger} \bar{v}^s(\mathbf{p}) e^{i\mathbf{p} \cdot \mathbf{x}} \right)$$ we have $$H = \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( - \gamma^i p_i +m) u^s(\mathbf{p}) + a_{\mathbf{p}}^{s'\dagger} a_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( \gamma^i p_i +m) v^s(\mathbf{p})$$ $$= \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( \gamma^0 p_0 ) u^s(\mathbf{p}) + a_{\mathbf{p}}^{s'\dagger} a_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( - \gamma^0 p_0 ) v^s(\mathbf{p})$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} E_{\mathbf{p}} ( b_{\mathbf{p}}^{s} b_{\mathbf{p}}^{s\dagger} - a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} )$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} - E_{\mathbf{p}} (b_{\mathbf{p}}^{s\dagger} b_{\mathbf{p}}^{s} + a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} ) - \infty$$ Changing anticommutator into commutator will make the spectrum unbounded. 4 revise quantization procedure in Schrodinger picture edited Feb 26 '14 at 9:08 Rodriguez 15622 silver badges77 bronze badges Yes. Though the energy will not be unbounded, but bounded from above, if my calculation is correct. For real scalar field under $$(+---)$$ metric, besides the negative classical kinetic energy for the Lagrangian $$\mathcal{L}=-\frac{1}{2} \partial^{\mu} \phi \partial_{\mu} \phi - \frac{1}{2} m^2 \phi^2 \tag{1}$$, the classical equation of motion will be $$(\square - m^2 )\phi=0 . \tag{2}$$ For plane wave $$\phi ~\sim e^{ipx}$$, it gives $$p^2+m^2 = (p^0)^2 - \mathbf{p}^2+m^2=0$$ which is inconsistent with relativistic energy momentum relation. I am not sure if it is necessary to quantize it. Though the energy-momentum-relation argument will not work for the Dirac field, we can quantize it to see the energy will be negative definite. $$\mathcal{L} = \bar{\psi}( -i \gamma^{\mu} \partial_{\mu} - m ) \psi \tag{3}$$ The classical equation of motion is $$(i \gamma^{\mu} \partial_{\mu} +m) \psi=0 \tag{4}$$ To preserve all properties of $$u(p)$$ and $$v(p)$$, we define $$\psi =: u(p) e^{ipx}$$$$\psi =: u(\mathbf{p}) e^{ipx}$$ $$\psi =: v(p) e^{-ipx}$$$$\psi =: v(\mathbf{p}) e^{-ipx}$$ Thus we can replace the $$u(p)$$$$u(\mathbf{p})$$ as $$v(p)$$$$v(\mathbf{p})$$ and $$v(p)$$$$v(\mathbf{p})$$ as $$u(p)$$$$u(\mathbf{p})$$ in the expansion of $$\psi$$ and $$\bar{\psi}$$. By $$\pi = -i \bar{\psi} \gamma^0$$ then $$H= \int d^3 x \bar{\psi} ( i \gamma^i \partial_i \psi + m ) \psi$$ Plug in $$\psi = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( a_{\mathbf{p}}^s v^s (\mathbf{p}) e^{-ipx} + b_{\mathbf{p}}^{s\dagger} u^s(\mathbf{p}) e^{ipx} \right)$$expansion of spinors in Schrodinger picture $$\bar{\psi} = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( b_{\mathbf{p}}^s \bar{u}^s (\mathbf{p}) e^{-ipx} + a_{\mathbf{p}}^{s\dagger} \bar{v}^s(\mathbf{p}) e^{ipx} \right)$$$$\psi = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( a_{\mathbf{p}}^s v^s (\mathbf{p}) e^{-i\mathbf{p} \cdot \mathbf{x} } + b_{\mathbf{p}}^{s\dagger} u^s(\mathbf{p}) e^{i \mathbf{p} \cdot \mathbf{x} } \right)$$ $$\bar{\psi} = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( b_{\mathbf{p}}^s \bar{u}^s (\mathbf{p}) e^{-i\mathbf{p} \cdot \mathbf{x}} + a_{\mathbf{p}}^{s\dagger} \bar{v}^s(\mathbf{p}) e^{i\mathbf{p} \cdot \mathbf{x}} \right)$$ we have $$H = \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( - \gamma^i p_i +m) u^s(p) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( \gamma^i p_i +m) v^s(p)$$$$H = \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( - \gamma^i p_i +m) u^s(\mathbf{p}) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( \gamma^i p_i +m) v^s(\mathbf{p})$$ $$= \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( \gamma^0 p_0 ) u^s(p) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( - \gamma^0 p_0 ) v^s(p)$$$$= \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( \gamma^0 p_0 ) u^s(\mathbf{p}) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( - \gamma^0 p_0 ) v^s(\mathbf{p})$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} E_{\mathbf{p}} ( b_{\mathbf{p}}^{s} b_{\mathbf{p}}^{s\dagger} - a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} )$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} - E_{\mathbf{p}} (b_{\mathbf{p}}^{s\dagger} b_{\mathbf{p}}^{s} + a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} ) - \infty$$ Changing anticommutator into commutator will make the spectrum unbounded. Yes. Though the energy will not be unbounded, but bounded from above, if my calculation is correct. For real scalar field under $$(+---)$$ metric, besides the negative classical kinetic energy for the Lagrangian $$\mathcal{L}=-\frac{1}{2} \partial^{\mu} \phi \partial_{\mu} \phi - \frac{1}{2} m^2 \phi^2 \tag{1}$$, the classical equation of motion will be $$(\square - m^2 )\phi=0 . \tag{2}$$ For plane wave $$\phi ~\sim e^{ipx}$$, it gives $$p^2+m^2 = (p^0)^2 - \mathbf{p}^2+m^2=0$$ which is inconsistent with relativistic energy momentum relation. I am not sure if it is necessary to quantize it. Though the energy-momentum-relation argument will not work for the Dirac field, we can quantize it to see the energy will be negative definite. $$\mathcal{L} = \bar{\psi}( -i \gamma^{\mu} \partial_{\mu} - m ) \psi \tag{3}$$ The classical equation of motion is $$(i \gamma^{\mu} \partial_{\mu} +m) \psi=0 \tag{4}$$ To preserve all properties of $$u(p)$$ and $$v(p)$$, we define $$\psi =: u(p) e^{ipx}$$ $$\psi =: v(p) e^{-ipx}$$ Thus we can replace the $$u(p)$$ as $$v(p)$$ and $$v(p)$$ as $$u(p)$$ in the expansion of $$\psi$$ and $$\bar{\psi}$$. By $$\pi = -i \bar{\psi} \gamma^0$$ then $$H= \int d^3 x \bar{\psi} ( i \gamma^i \partial_i \psi + m ) \psi$$ Plug in $$\psi = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( a_{\mathbf{p}}^s v^s (\mathbf{p}) e^{-ipx} + b_{\mathbf{p}}^{s\dagger} u^s(\mathbf{p}) e^{ipx} \right)$$ $$\bar{\psi} = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( b_{\mathbf{p}}^s \bar{u}^s (\mathbf{p}) e^{-ipx} + a_{\mathbf{p}}^{s\dagger} \bar{v}^s(\mathbf{p}) e^{ipx} \right)$$ we have $$H = \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( - \gamma^i p_i +m) u^s(p) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( \gamma^i p_i +m) v^s(p)$$ $$= \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( \gamma^0 p_0 ) u^s(p) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( - \gamma^0 p_0 ) v^s(p)$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} E_{\mathbf{p}} ( b_{\mathbf{p}}^{s} b_{\mathbf{p}}^{s\dagger} - a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} )$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} - E_{\mathbf{p}} (b_{\mathbf{p}}^{s\dagger} b_{\mathbf{p}}^{s} + a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} ) - \infty$$ Changing anticommutator into commutator will make the spectrum unbounded. Yes. Though the energy will not be unbounded, but bounded from above, if my calculation is correct. For real scalar field under $$(+---)$$ metric, besides the negative classical kinetic energy for the Lagrangian $$\mathcal{L}=-\frac{1}{2} \partial^{\mu} \phi \partial_{\mu} \phi - \frac{1}{2} m^2 \phi^2 \tag{1}$$, the classical equation of motion will be $$(\square - m^2 )\phi=0 . \tag{2}$$ For plane wave $$\phi ~\sim e^{ipx}$$, it gives $$p^2+m^2 = (p^0)^2 - \mathbf{p}^2+m^2=0$$ which is inconsistent with relativistic energy momentum relation. I am not sure if it is necessary to quantize it. Though the energy-momentum-relation argument will not work for the Dirac field, we can quantize it to see the energy will be negative definite. $$\mathcal{L} = \bar{\psi}( -i \gamma^{\mu} \partial_{\mu} - m ) \psi \tag{3}$$ The classical equation of motion is $$(i \gamma^{\mu} \partial_{\mu} +m) \psi=0 \tag{4}$$ To preserve all properties of $$u(p)$$ and $$v(p)$$, we define $$\psi =: u(\mathbf{p}) e^{ipx}$$ $$\psi =: v(\mathbf{p}) e^{-ipx}$$ Thus we can replace the $$u(\mathbf{p})$$ as $$v(\mathbf{p})$$ and $$v(\mathbf{p})$$ as $$u(\mathbf{p})$$ in the expansion of $$\psi$$ and $$\bar{\psi}$$. By $$\pi = -i \bar{\psi} \gamma^0$$ then $$H= \int d^3 x \bar{\psi} ( i \gamma^i \partial_i \psi + m ) \psi$$ Plug in expansion of spinors in Schrodinger picture $$\psi = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( a_{\mathbf{p}}^s v^s (\mathbf{p}) e^{-i\mathbf{p} \cdot \mathbf{x} } + b_{\mathbf{p}}^{s\dagger} u^s(\mathbf{p}) e^{i \mathbf{p} \cdot \mathbf{x} } \right)$$ $$\bar{\psi} = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( b_{\mathbf{p}}^s \bar{u}^s (\mathbf{p}) e^{-i\mathbf{p} \cdot \mathbf{x}} + a_{\mathbf{p}}^{s\dagger} \bar{v}^s(\mathbf{p}) e^{i\mathbf{p} \cdot \mathbf{x}} \right)$$ we have $$H = \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( - \gamma^i p_i +m) u^s(\mathbf{p}) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( \gamma^i p_i +m) v^s(\mathbf{p})$$ $$= \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( \gamma^0 p_0 ) u^s(\mathbf{p}) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( - \gamma^0 p_0 ) v^s(\mathbf{p})$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} E_{\mathbf{p}} ( b_{\mathbf{p}}^{s} b_{\mathbf{p}}^{s\dagger} - a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} )$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} - E_{\mathbf{p}} (b_{\mathbf{p}}^{s\dagger} b_{\mathbf{p}}^{s} + a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} ) - \infty$$ Changing anticommutator into commutator will make the spectrum unbounded. 3 added 115 characters in body edited Feb 26 '14 at 8:27 Rodriguez 15622 silver badges77 bronze badges Yes. Though the energy will not be unbounded, but bounded from above, if my calculation is correct. For real scalar field under $$(+---)$$ metric, besides the negative classical kinetic energy for the Lagrangian $$\mathcal{L}=-\frac{1}{2} \partial^{\mu} \phi \partial_{\mu} \phi - \frac{1}{2} m^2 \phi^2 \tag{1}$$, the classical equation of motion will be $$(\square - m^2 )\phi=0 . \tag{2}$$ For plane wave $$\phi ~\sim e^{ipx}$$, it gives $$p^2+m^2 = (p^0)^2 - \mathbf{p}^2+m^2=0$$ which is inconsistent with relativistic energy momentum relation. I am not sure if it is necessary to quantize it. Though the energy-momentum-relation argument will not work for the Dirac field, we can quantize it to see the energy will be negative definite. $$\mathcal{L} = \bar{\psi}( -i \gamma^{\mu} \partial_{\mu} - m ) \psi \tag{3}$$ The classical equation of motion is $$(i \gamma^{\mu} \partial_{\mu} +m) \psi=0 \tag{4}$$ To preserve all properties of $$u(p)$$ and $$v(p)$$, we define This implies$$\psi =: u(p) e^{ipx}$$ $$\psi =: v(p) e^{-ipx}$$ Thus we can replace the $$u(p)$$ as $$v(p)$$ and $$v(p)$$ as $$u(p)$$ in the expansion of $$\psi$$ and $$\bar{\psi}$$. By $$\pi = -i \bar{\psi} \gamma^0$$ then $$H= \int d^3 x \bar{\psi} ( i \gamma^i \partial_i \psi + m ) \psi$$ Plug in $$\psi = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( a_{\mathbf{p}}^s v^s (\mathbf{p}) e^{-ipx} + b_{\mathbf{p}}^{s\dagger} u^s(\mathbf{p}) e^{ipx} \right)$$ $$\bar{\psi} = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( b_{\mathbf{p}}^s \bar{u}^s (\mathbf{p}) e^{-ipx} + a_{\mathbf{p}}^{s\dagger} \bar{v}^s(\mathbf{p}) e^{ipx} \right)$$ we have $$H = \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( - \gamma^i p_i +m) u^s(p) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( \gamma^i p_i +m) v^s(p)$$ $$= \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( \gamma^0 p_0 ) u^s(p) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( - \gamma^0 p_0 ) v^s(p)$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} E_{\mathbf{p}} ( b_{\mathbf{p}}^{s} b_{\mathbf{p}}^{s\dagger} - a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} )$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} - E_{\mathbf{p}} (b_{\mathbf{p}}^{s\dagger} b_{\mathbf{p}}^{s} + a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} ) - \infty$$ Changing anticommutator into commutator will make the spectrum unbounded. Yes. Though the energy will not be unbounded, but bounded from above, if my calculation is correct. For real scalar field under $$(+---)$$ metric, besides the negative classical kinetic energy for the Lagrangian $$\mathcal{L}=-\frac{1}{2} \partial^{\mu} \phi \partial_{\mu} \phi - \frac{1}{2} m^2 \phi^2 \tag{1}$$, the classical equation of motion will be $$(\square - m^2 )\phi=0 . \tag{2}$$ For plane wave $$\phi ~\sim e^{ipx}$$, it gives $$p^2+m^2 = (p^0)^2 - \mathbf{p}^2+m^2=0$$ which is inconsistent with relativistic energy momentum relation. I am not sure if it is necessary to quantize it. Though the energy-momentum-relation argument will not work for the Dirac field, we can quantize it to see the energy will be negative definite. $$\mathcal{L} = \bar{\psi}( -i \gamma^{\mu} \partial_{\mu} - m ) \psi \tag{3}$$ The classical equation of motion is $$(i \gamma^{\mu} \partial_{\mu} +m) \psi=0 \tag{4}$$ This implies we can replace the $$u(p)$$ as $$v(p)$$ and $$v(p)$$ as $$u(p)$$ in the expansion of $$\psi$$ and $$\bar{\psi}$$. By $$\pi = -i \bar{\psi} \gamma^0$$ then $$H= \int d^3 x \bar{\psi} ( i \gamma^i \partial_i \psi + m ) \psi$$ Plug in $$\psi = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( a_{\mathbf{p}}^s v^s (\mathbf{p}) e^{-ipx} + b_{\mathbf{p}}^{s\dagger} u^s(\mathbf{p}) e^{ipx} \right)$$ $$\bar{\psi} = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( b_{\mathbf{p}}^s \bar{u}^s (\mathbf{p}) e^{-ipx} + a_{\mathbf{p}}^{s\dagger} \bar{v}^s(\mathbf{p}) e^{ipx} \right)$$ we have $$H = \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( - \gamma^i p_i +m) u^s(p) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( \gamma^i p_i +m) v^s(p)$$ $$= \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( \gamma^0 p_0 ) u^s(p) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( - \gamma^0 p_0 ) v^s(p)$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} E_{\mathbf{p}} ( b_{\mathbf{p}}^{s} b_{\mathbf{p}}^{s\dagger} - a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} )$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} - E_{\mathbf{p}} (b_{\mathbf{p}}^{s\dagger} b_{\mathbf{p}}^{s} + a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} ) - \infty$$ Changing anticommutator into commutator will make the spectrum unbounded. Yes. Though the energy will not be unbounded, but bounded from above, if my calculation is correct. For real scalar field under $$(+---)$$ metric, besides the negative classical kinetic energy for the Lagrangian $$\mathcal{L}=-\frac{1}{2} \partial^{\mu} \phi \partial_{\mu} \phi - \frac{1}{2} m^2 \phi^2 \tag{1}$$, the classical equation of motion will be $$(\square - m^2 )\phi=0 . \tag{2}$$ For plane wave $$\phi ~\sim e^{ipx}$$, it gives $$p^2+m^2 = (p^0)^2 - \mathbf{p}^2+m^2=0$$ which is inconsistent with relativistic energy momentum relation. I am not sure if it is necessary to quantize it. Though the energy-momentum-relation argument will not work for the Dirac field, we can quantize it to see the energy will be negative definite. $$\mathcal{L} = \bar{\psi}( -i \gamma^{\mu} \partial_{\mu} - m ) \psi \tag{3}$$ The classical equation of motion is $$(i \gamma^{\mu} \partial_{\mu} +m) \psi=0 \tag{4}$$ To preserve all properties of $$u(p)$$ and $$v(p)$$, we define $$\psi =: u(p) e^{ipx}$$ $$\psi =: v(p) e^{-ipx}$$ Thus we can replace the $$u(p)$$ as $$v(p)$$ and $$v(p)$$ as $$u(p)$$ in the expansion of $$\psi$$ and $$\bar{\psi}$$. By $$\pi = -i \bar{\psi} \gamma^0$$ then $$H= \int d^3 x \bar{\psi} ( i \gamma^i \partial_i \psi + m ) \psi$$ Plug in $$\psi = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( a_{\mathbf{p}}^s v^s (\mathbf{p}) e^{-ipx} + b_{\mathbf{p}}^{s\dagger} u^s(\mathbf{p}) e^{ipx} \right)$$ $$\bar{\psi} = \int \frac{ d^3 p }{ (2\pi)^3} \frac{1}{ \sqrt{2 E_{\mathbf{p}}}} \sum_s \left( b_{\mathbf{p}}^s \bar{u}^s (\mathbf{p}) e^{-ipx} + a_{\mathbf{p}}^{s\dagger} \bar{v}^s(\mathbf{p}) e^{ipx} \right)$$ we have $$H = \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( - \gamma^i p_i +m) u^s(p) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( \gamma^i p_i +m) v^s(p)$$ $$= \sum_{ss'} \int \frac{d^3p}{ (2\pi)^3 2E_{\mathbf{p}} } b_{\mathbf{p}}^{s'} b_{\mathbf{p}}^{s\dagger} \bar{u}^{s'}(\mathbf{p}) ( \gamma^0 p_0 ) u^s(p) + a_{\mathbf{p}}^{s'\dagger} s_{\mathbf{p}}^{s} \bar{v}^{s'}(\mathbf{p}) ( - \gamma^0 p_0 ) v^s(p)$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} E_{\mathbf{p}} ( b_{\mathbf{p}}^{s} b_{\mathbf{p}}^{s\dagger} - a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} )$$ $$= \sum_s \int \frac{ d^3p}{ (2\pi)^3} - E_{\mathbf{p}} (b_{\mathbf{p}}^{s\dagger} b_{\mathbf{p}}^{s} + a_{\mathbf{p}}^{s\dagger} a_{\mathbf{p}}^{s} ) - \infty$$ Changing anticommutator into commutator will make the spectrum unbounded. 2 added 34 characters in body edited Feb 26 '14 at 4:31 Rodriguez 15622 silver badges77 bronze badges 1 answered Feb 26 '14 at 3:16 Rodriguez 15622 silver badges77 bronze badges


text:
# Latex Bibliography Style Apacite Package

Apacite provides a BibTeX style and a LaTeX package which are designed to match the requirements of the American Psychological Association's style for citations. The package follows the 6th edition of the APA manual, and is designed to work with the apa6 class. A test document is provided. The package is compatible with chapterbib and (to some extent) with hyperref (for limits of compatibility, see the documentation). The package also includes a means of generating an author index for a document.

Download the contents of this package in one zip archive (798.5k).

more

## Suggestions

Maybe you are interested in the following packages as well.

more

At the beginning I used

and I really like the author-year notation. The result looks fine, but it's English and this time my seminar thesis is written in german. Within the text and the bibliography I need to change the conjunction "and" into "und", and substitute some more vocab like "pages, volume" etc.

There were some suggestions to this issue and the most promising is:

I only need two implementations:

• a cite-version "Smith et al. (2000)" in running text
• and a citep-version "(Smith et al., 2000)" in parathesis

At the moment the apacite-solution alternative citet doesnt work for me, so all my \cite-refs are missing (questions marks within my text).

Another small but important issue are two options for the \citep

• \citep[zitiert nach:][S. 6]{sedlacek2011economics}

Doesn't work anymore. On the other hand a minimal example works very well for me:

My main-file blg. report doesn't list any mistakes. I still try to figure out the MWE/mainfile difference ..

Thank you

biblatexbibliographiesnatbibapa-style


text:
# Exponential distribution

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

2010 Mathematics Subject Classification: Primary: 60E99 [MSN][ZBL]

A continuous distribution of a random variable $X$ defined by the density

$$\tag{1 } p( x) = \left \{ \begin{array}{ll} \lambda e ^ {- \lambda x } , & x \geq 0 , \\ 0 , & x < 0. \\ \end{array} \right .$$

The density $p( x)$ is dependent on the positive scale parameter $\lambda$. The formula for the moments is ${\mathsf E} X ^ {n} = n!/ \lambda ^ {n}$, and, in particular, the expectation equals ${\mathsf E} X = 1/ \lambda$ and the variance equals ${\mathsf D} X = 1/ \lambda ^ {2}$; the characteristic function is $( 1- it/ \lambda ) ^ {-} 1$.

The exponential distribution belongs to the family of gamma-distributions (cf. Gamma-distribution) which are defined by the densities

$$p( x) = \frac{\lambda ^ \alpha x ^ {\alpha - 1 } }{\Gamma ( \alpha ) } e ^ {- \lambda x } ,\ \ x \geq 0,\ \alpha > 0 ;$$

the $n$- fold convolution of the density (1) is equal to the gamma-density with the same parameter $\lambda$ and with $\alpha = n$.

The exponential distribution is the unique distribution having the property of no after-effect: For any $x > 0$, $y > 0$ one has

$$\tag{2 } {\mathsf P} \{ X > x + y \mid X > y \} = \ {\mathsf P} \{ X > x \} ,$$

where ${\mathsf P} \{ X > x + y \mid X > y \}$ is the conditional probability of the event $X > x + y$ subject to the condition $X > y$. Property (2) is also called the lack-of-memory property.

In a homogeneous Poisson process, the distances between successive events have an exponential distribution. Conversely, a renewal process with exponential lifetime (1) is a Poisson process. An exponential distribution often arises as a limit process on the superposition or extension of renewal processes, as well as in high-level intersection problems in various random-path schemes, in critical branching processes, etc.

The above features explain why the exponential distribution is widely used in calculating various systems in queueing theory and reliability theory. One assumes that the lifetimes of the devices are independent random variables with exponential distributions, and then the property (2) enables one to examine a queueing system by means of finite or denumerable Markov chains with continuous time. Similarly, one uses Markov chains in reliability theory, where the fault-free operating times of the individual devices can often be taken as independent and as having exponential distributions.

#### References

 [F] W. Feller, "An introduction to probability theory and its applications" , 2 , Wiley (1971)


text:
# F(x) as function or value of the function?

Tags:
1. Dec 31, 2014

### Mr Davis 97

I have a question regarding notational conventions. From previous research, I have determined that, for the general function ƒ(x), x is the input, ƒ is the relation or "law of correspondence," and ƒ(x) is the output, the result of applying the rule ƒ to the variable x. However, almost ubiquitously in other texts and websites, they refer to the entire function as ƒ(x), e.g., "let's examine the function ƒ(x) = x2 + 1." My question, then, is what is right? Why do some authors say that ƒ is the function, while others say that ƒ(x) is the function? This distinction might seem rather pedantic, but I am nonetheless still bothered by it.

2. Dec 31, 2014

### jbunniii

"Let's examine the function $f(x) = x^2 + 1$" is a slightly sloppy (but much shorter) way of saying "let's examine the function $f:\mathbb{R} \to \mathbb{R}$ defined by $f(x) = x^2 + 1$". The shorter statement doesn't generally cause any ambiguity as long as the domain and codomain are understood. It would be tedious if we always insisted on pedantic language such as "the function $f:\mathbb{R} \to\mathbb{R}$ defined by $f(x) = e^x$" instead of simply "the function $e^x$".

Of course, there may be cases where we need to distinguish between the function $f$ and its value $f(x)$ at some particular $x$, and we should write more carefully in those cases.

3. Dec 31, 2014

### Joshua L

$f(x)$ is a notation that to implies "function $f$ is dependent of or written in terms of the variable $x$." The notation makes it somewhat easier to see how the variable plays into determining each value of the function.

I'm going to use an example. It will involve physics because I like physics.

Let's look at the equation of the position of a massive particle moving at a constant acceleration. $$\vec r_f = \vec r_i + \vec v_i \Delta t + \frac{1}{2} \vec a \Delta t^2$$
Since the initial position, the initial velocity, the initial time measurement, and the constant acceleration are "fixed" variables, unlike the variable $t$, then it can be rewritten as $$\vec r (t) = \vec r_i + \vec v_i (t-t_i) + \frac{1}{2} \vec a (t-t_i)^2$$
This shows that the particle's position $\vec r$ is dependent on time $t$.

Authors have the ability to switch the two notations you've listed above because they can. It implies the same thing.

4. Dec 31, 2014

### pwsnafu

A function is a triple $(A,B,\Gamma)$ such that
1. A is a set called the domain,
2. B is a set called the co-domain,
3. $\Gamma \subset A \times B$ is a set called the graph,
4. for all $x \in A$ there exists $(x, y) \in \Gamma$ ("every function is defined on its domain"),
5. if $(x_1, y) \in \Gamma$ and $(x_2, y) \in \Gamma$ then $x_1 = x_2$ ("the vertical line test").
Notations such as $f(x)$ are attempts at defining functions without writing the whole thing down. They are just shortcuts. Different books will have different conventions.

Now as to the question concerning $f$ vs $f(x)$. When I was much younger I came across an abstract algebra textbook focusing on ring theory dating back to the early 1930s. It used the $xf$ notation: that is $x$ is a point in the domain and $f$ is the map (homomorphism). The clear advantage is composition. $xfg$ is apply f then g. Today we have to write $g\circ f$. This is great for algebra, horrible for applied analysis. What is $2f$? Twice the functions? Or the function evaluated at 2? To top it off, as matrices became more standard, we have abandoned $xf$ notation for $f(x)$, but the split between "f" and "x" stayed in pure mathematics circles. So as a consequence pure mathematicians consider the function to be "f" and the valued to be "f(x)".

But I can't find any old applied mathematics papers that doesn't use "f(x)" is the function. For an applied mathematician what a function does at each point is of prime importance. It's my belief that applied mathematics has used "f(x)" before pure mathematics has used "f" as the function.

5. Dec 31, 2014

### Mr Davis 97

That makes sense. I've seen it written in some definitions for a function, that the domain and range must be declared in order to completely specify a function. But there still is another distinction which I am (pedantically) concerned about. Is there a way to specify what a function $f(~)$ does without having to place a dummy variable, $x$, in its argument? Is it even possible to talk about what the majority of functions do to inputs without having to apply the function to some free variable, like $x$?

Thank you, this really helps. I did sort of realize that applied texts, such as physics textbooks, often always define a function simply as $f(x)$, while pure mathematics texts used the convention of $f$ being the function while $f(x)$, the value at $x$.

6. Dec 31, 2014

### pwsnafu

Some functions sure. The identity function from $X\to X$ is the function with the property such that for any $f:X\to X$ we have $f \circ Id = Id \circ f$.
Remember category theory can defined things like injectivity without any reference to the variable.

Majority? Depends on what information you want from it, right?

Edit for grammar

Last edited: Dec 31, 2014
7. Dec 31, 2014

### Mr Davis 97

Majority as in functions you learn in high school, e.g., polynomial, transcendental, and combinations of the two. To talk about what the function, say, $f(x) = sin(x) + x^2 - 1$ does without reference to any input, what could we do? How could be talk about the $f$ function without reference to any dummy input? Do we only usually define functions like $f(x)$, with reference to the dummy input, because that's easier than saying $f(~): \text{take the sin, add the square, subtract 1}$? (Please excuse if this is a dumb question; I'm just trying to figure out why we define functions based on the output, $f(x)$, rather than just $f(~)$).

8. Dec 31, 2014

### Stephen Tashi

As you said, it usually is simplest to define a function using a dummy variable. If you think about computer programming a function is given a dummy argument so that the procedure of computing the function can be defined. After a function has been defined, it can be referred to as "the function f" or "the function f()". That notation is actually used.

The question of "what is a variable" is not simple to answer. Within a certain "scope" a variable is treated like a symbol that represents a specific number. In computer programming, the "scope" of a variable is clear. If you see an "x" within the scope of a function, you know it isn't the "x" referred to in the code for a different function. In mathematical writing, many authors are careless. For example, you may be told that f(x) = 3x + 2 and g(y) = sin(y) + y^2 are "functions". But then your are told something like y = 2x, enforcing a relation between y and x,. Many discussions in physics that speak of "functions" are actually discussion of "relations" or discussion of several different functions, all denoted as the same function. The use of specific dummy variables facilitates such ambiguous but sometimes convenient writing.


text:
# Minimum Cuts

Let $$G$$ be an arbitrary graph, and let $$s$$ and $$t$$ be two vertices of $$G$$. An $$(s, t)$$-cut (or more formally, an $$(s,t)$$-edge-cut) is a subset $$X$$ of the edges of $$G$$ that intersects every path from $$s$$ to $$t$$ in $$G$$. A minimum $$(s,t)$$-cut is an $$(s,t)$$-cut of minimum size, or of minimum total weight if the edges of $$G$$ are weighted. (In this context, edge weights are normally called capacities.)

The fastest method to compute minimum $$(s, t)$$-cuts in arbitrary graphs is to compute a maximum $$(s, t)$$-flow and then exploit the classical proof of the maxflow-mincut theorem. In undirected planar graphs, however, this dependency is reversed; the fastest method to compute maximum flows actually computes minimum cuts first.

## Duality: Shortest essential cycle

Let $$\Sigma$$ be an undirected planar map, each of whose edges $$e$$ has a non-negative capacity $$c(e)$$, and let $$s$$ and $$t$$ be vertices of $$\Sigma$$. The first step in our fast minimum-cut algorithm is to view the problem through the lens of duality. It is helpful to think of the source vertex $$s$$ and the target vertex $$t$$ as punctures or obstacles — points that are missing from the plane — and similarly to think of the corresponding faces $$s^*$$ and $$t^*$$ as holes in the dual map $$\Sigma^*$$. In other words, we should think of the dual map $$\Sigma^*$$ as a decomposition of the annulus $$A = \mathbb{R}^2 \setminus (s^*\cup t^*)$$ rather than a map on the plane or a disk. Without loss of generality, assume that $$t^*$$ is the outer face of $$\Sigma^*$$.

A simple cycle $$\gamma$$ in $$\Sigma^*$$ is essential if it satisfies any of the following equivalent conditions:

• $$\gamma$$ separates $$s^*$$ from $$t^*$$.
• $$\gamma$$ has winding number $$\pm 1$$ around $$s^*$$.
• $$\gamma$$ is homotopic in $$A$$ to the boundary of $$s^*$$.
• $$\gamma$$ is homotopic in $$A$$ to the boundary of $$t^*$$.

Each edge $$e^*$$ in the dual map $$\Sigma^*$$ has a cost or length $$c^*(e^*)$$ equal to the capacity of the corresponding primal edge $$e$$. Whitney’s duality between simple cycles in $$\Sigma$$ and and minimal cuts (bonds) in $$\Sigma^*$$ immediately implies the following:

Lemma:
A subset $$X$$ of edges is a minimum $$(s,t)$$-cut in $$\Sigma$$ if and only if the corresponding set $$X^*$$ of dual edges is a minimum-cost essential cycle in $$\Sigma^*$$.

## Crossing at most once

Now let $$\pi$$ be a shortest path in $$\Sigma^*$$ from any vertex of $$s^*$$ to any vertex of $$t^*$$. We can measure the winding number of any directed cycle $$\gamma$$ in $$\Sigma^*$$ by counting the number of times $$\gamma$$ crosses $$\pi$$ in either direction. We have to define “crossing” carefully here, because $$\gamma$$ and $$\pi$$ could share edges.

Suppose $$\pi = p_0\mathord\to p_1\mathord\to \cdots \mathord\to p_k$$, where $$p_0$$ lies on $$s^*$$ and $$p_k$$ lies on $$t^*$$. To simplify the following definition, we add two “ghost” darts $$p_{-1}\mathord\to p_0$$ and $$p_k\mathord\to p_{k+1}$$, where $$p_{-1}$$ lies inside $$s^*$$ and $$p_{k+1}$$ lies inside $$t^*$$. We say that a dart $$q\mathord\to p_i$$ enters $$\pi$$ from the left (resp. from the right) if the darts $$p_{i-1}\mathord\to p_i$$, $$q\mathord\to p_i$$, and $$p_{i+1}\mathord\to p_i$$ are ordered clockwise (resp. counterclockwise) around $$p_i$$. Symmetrically, we say that a dart leaves $$\pi$$ to the left (resp. to the right) if its reversal enters $$\pi$$ from the left (resp. from the right). The same dart can leave $$\pi$$ to the right and enter $$\pi$$ to the left, or vice versa.

A positive crossing between $$\pi$$ and $$\gamma$$ is a subpath of $$\gamma$$ that starts with a dart entering $$\pi$$ from the right and ends with a dart leaving $$\pi$$ to the left, and a negative crossing is defined similarly. Intuitively, for purposes of defining crossings, we are shifting the path $$\pi$$ very slightly to the left, so that it intersects the edges of $$\Sigma^*$$ transversely. It follows that the winding number $$\textsf{wind}(\gamma, s^*)$$ is the number of darts in $$\gamma$$ that leave $$\pi$$ to the left, minus the number of darts in $$\gamma$$ that enter $$\pi$$ from the left.

Crossing Lemma [Itai and Shiloach (1979)]:
The shortest essential cycle in $$\Sigma^*$$ crosses $$\pi$$ exactly once.
Proof:
We follow the same intuition that we used for shortest homotopic paths in the plane. Let $$\gamma$$ be any essential cycle in $$\Sigma^*$$, directed so that $$\textsf{wind}(\gamma, s^*) = 1$$, that crosses $$\pi$$ more than once. Then somewhere $$\gamma$$ must have be a negative crossing followed immediately by a positive crossing. It follows that $$\gamma$$ has a subpath $$p_i\mathord\to q\mathord\to\cdots\mathord\to r\mathord\to p_j$$, where $$p_i\mathord\to q$$ leaves $$\pi$$ to the right, $$r\mathord\to p_j$$ enters $$\pi$$ from the right. Let $$\gamma’$$ be the cycle obtained from $$\gamma$$ by replacing this subpath with the subpath of $$\pi$$ from $$p_i$$ to $$p_j$$. Because $$\pi$$ is a shortest path, $$\gamma’$$ must be shorter than $$\gamma$$. We conclude that $$\gamma’$$ is not the shortest essential cycle in $$\Sigma^*$$. $$\qquad\square$$

## Slicing along a path

Now let $$\Delta := \Sigma^* \backslash\!\!\backslash \pi$$ denote the planar map obtained by slicing the annular map $$\Sigma^*$$ along path $$\pi$$. The slicing operation replaces $$\pi$$ with two copies $$\pi^+$$ and $$\pi^-$$. Then for every vertex $$p_i$$ of $$\pi$$, all edges incident to $$p_i$$ on the left are redirected to $$p_i^+$$, and all edges incident to $$p_i$$ on the left are redirected to $$p_i^-$$. The channel between two two paths $$\pi^+$$ and $$\pi^-$$ joins $$s^*$$ and $$t^*$$ into a single outer face. Thus, we should think of $$\Delta$$ as being embedded on a disk. Every face of $$\Sigma^*$$ except $$s^*$$ and $$t^*$$ appears as a face in $$\Delta$$.

For any index $$i$$, let $$\sigma_i$$ denote the shortest path in $$\Delta$$ from $$p_i^+$$ to $$p_i^-$$. The shortest essential cycle $$\gamma$$ in $$\Sigma^*$$ appears in $$\Delta$$ as one of the shortest paths $$\sigma_i$$. Thus, to compute the minimum $$(s,t)$$-cut in our original planar map $$\Sigma$$, it suffices to compute the length of every shortest path $$\sigma_i$$ in $$\Delta$$.

## Algorithms

The simplest way to compute these $$k$$ shortest-path distances is to run Dijkstra’s algorithm at each vertex $$p_i^+$$. Assuming $$\pi$$ has $$k$$ edges, so there are $$k+1$$ terminal pairs $$p_i^\pm$$, this algorithm runs in $$O(kn\log n)$$ time, which is $$O(n^2\log n)$$ time in the worst case. We can reduce the running time to $$O(kn\log\log n)$$ using the faster shortest-path algorithm we described in the previous lecture note, or even to $$O(kn)$$ using a linear-time shortest-path algorithm.

Alternatively, we can compute all $$k$$ of these shortest paths using a multiple-source shortest-paths algorithm. The parametric MSSP algorithms of Klein and Cabello et al. both require $$O(n\log n)$$ time, which is faster in the worst case than $$O(kn)$$, but slower when $$k$$ is small. In particular, even when $$k=2$$, that algorithm could require $$\Omega(n)$$ pivots. The more recent recursive contraction algorithm runs in $$O(n\log n\log k)$$ time if we use Dijkstra’s algorithm to compute shortest paths, or in $$O(n\log k)$$ time if we use a linear-time shortest-path algorithm.

An older and simpler divide-and-conquer algorithm, proposed by John Reif in 1983, exactly matches the recursive contraction MSSP algorithm. Reif’s algorithm computes the median shortest path $$\sigma_m$$, where $$m = \lfloor k/2 \rfloor$$, and then recurses in each component of the sliced map $$\Delta \backslash\!\!\backslash \sigma_m$$. One of these components contains the terminal pairs $$p_0^\pm, p_1^\pm, \dots, p^\pm_m$$; the other contains the terminal pairs $$p_m^\pm, p_{m+1}^\pm, \dots, p_k^\pm$$.

Reif’s algorithm falls back on Dijkstra’s algorithm in two base cases. First, if $$k\le 2$$, we can obviously compute each of the $$k$$ shortest paths directly. A more subtle base case happens when the “floor” and “ceiling” paths collide. Let $$\alpha$$ denote the boundary path from $$p_0^+$$ to $$p_0^-$$, and let $$\beta$$ denote the boundary path from $$p_k^+$$ to $$p_k^-$$. If $$\alpha$$ and $$\beta$$ share a vertex $$x$$, then for every index $$i$$ we have $$\textsl{dist}(p_i^+, p_i^-) = \textsl{dist}(p_i^+, x) + \textsl{dist}(x, p_i^-)$$; thus, instead of recursing, we can compute all $$k$$ shortest-path distances by computing a single shortest-path tree rooted at $$x$$. (This second base case is not necessary for the correctness of Reif’s algorithm, but it is necessary for efficiency.)

Let $$T(n,k)$$ denote the running time of Reif’s algorithm, where $$k+1$$ is the number of terminal pairs and $$n$$ is the total number of vertices in the map $$\Delta$$. This function obeys the recurrence $T(n,k) = T(n_1, \lfloor k/2 \lfloor) + T(n_1, \lceil k/2 \rceil) + O(n\log n).$ where $$n_1$$ and $$n_2$$ are the number of vertices in the two components of $$\Delta \backslash\!\!\backslash \sigma_m$$. The second base case ensures that each vertex and edge of $$\Delta$$ appears in at most $$O(1)$$ subproblems at any level of the recursion tree.1 Thus, the total work at any level of recursion is $$O(n\log n)$$. The recursion tree has depth at most $$O(\log k)$$, so the overall algorithm runs in $$O(n\log n\log k)$$ time.

If we use a linear-time shortest-path algorithm instead of Dijkstra, the running time improves to $$O(n\log k)$$. This improvement was first described by Greg Frederickson in 1987, as one of the earliest applications of $$r$$-divisions.

## Faster Shortest Paths with Negative Lengths

[[[This is still really sketchy!]]]

In the previous lecture note, we described a 2009 algorithm by Klein, Mozes, and Weimann to compute shortest paths in directed planar maps, possibly with negative dart lengths, in $$O(n\log^2 n)$$ time. In 2010, Shay Mozes and Christian Wulff-Nilsen improved this algorithm even further by using a good $$r$$-division at each level of recursion (with $$r \approx n/\log n$$) instead of just one separator cycle; their improved algorithm runs in $$O(n\log^2n /\log\log n)$$.

Recall the high-level description of the Klein-Mozes-Weimann algorithm:

1. Split the map $$\Sigma$$ into two smaller maps $$A$$ and $$B$$ along a single cycle separator $$S$$.
2. Recursively compute $$\textsf{dist}_A(r, A)$$ and $$\textsf{dist}_B(r, B)$$
3. Compute $$\textsf{dist}_A(S,S)$$ and $$\textsf{dist}_B(S,S)$$ using MSSP in $$O(n\log n)$$ time.
4. Compute $$\textsf{dist}_\Sigma(r,S)$$ using FR-Bellman-Ford in $$O(n\alpha(n))$$ time.
5. Compute $$\textsf{dist}_\Sigma(r,\Sigma)$$ using reweighting and $$r$$-divisions in $$O(n\log\log n)$$ time.
6. Compute $$\textsf{dist}_\Sigma(s,\Sigma)$$ using reweighting and $$r$$-divisions in $$O(n\log\log n)$$ time.

Steps 1, 4, and 5 of Mozes and Wulff-Nilsen’s faster algorithm are identical. In step 2, the only minor different is that we need to run MSSP around the boundary of each hole of each piece, instead of just once around the sole cycle separator. The total time for each piece is $$O(r\log r)$$, so the total time for this step across the entire $$r$$-division is $$O(n/r)\cdot O(r\log r) = O(n\log r)$$.

The only major change to the algorithm is in step 3; we need to modify FR-Bellman-Ford to deal with holes in each piece of the $$r$$-division.

In step 2, we compute all boundary-to-boundary distances within each piece of the $$r$$-division using MSSP. Specifically, within each piece $$P$$, we run MSSP around each of the $$O(1)$$ boundary cycles of $$P$$. The total time for each piece is $$O(r\log r)$$, so the total time for this step across the entire $$r$$-division is $$O(n/r)\cdot O(r\log r) = O(n\log r)$$.

Modifying step 3 is more interesting. Recall that a good $$r$$-division $$\mathcal{R}$$ breaks $$\Sigma$$ into $$O(n/r)$$ pieces, each with $$O(r)$$ vertices, $$O(\sqrt{r})$$ boundary vertices, and $$O(1)$$ holes; thus, overall, $$\mathcal{R}$$ has $$O(n/\sqrt{r})$$ boundary vertices, organized into $$O(n/r)$$ cycles. We construct a complete directed multigraph $$\hat{R}$$ over the boundary vertices of the $$r$$-division, with an edge $$u{\to}v$$ with weight $$\textsf{dist}_P(u,v)$$ for every piece $$P$$ containing both $$u$$ and $$v$$.

Bellman-Ford computes shortest-path distances in $$\hat{R}$$ in $$O(n/\sqrt{r})$$ phases; in each phase, we find and relax all tense edges. Our modification of FR-Bellman-Ford breaks each relaxation phase into subphases:

for each piece P:
for each boundary cycle A of P:
for each boundary cycle B of P:
relax every edge within P from A to B

We already described how to relax every edge in $$P$$ from a boundary cycle to itself in $$O(k\alpha(k))$$ time (or $$O(k)$$ expected time). It remains to show how to relax all edges from one boundary cycle $$A$$ to a different boundary cycle $$B$$.

Without loss of generality, let’s assume that $$A$$ and $$B$$ are the only two boundary cycles in $$P$$; we’ll treat any other boundary cycles as faces of $$P$$. Thus, $$P$$ is a map of an annulus, that is, a disk with one hole. Let $$a_1, a_2, \dots, a_k$$ and $$b_1, b_2, \dots, b_l$$ denote the sequences of vertices of $$A$$ and $$B$$, respectively.

Let $$\pi$$ be the shortest path in $$P$$ from $$a_1$$ to $$b_1$$, and let $$D = P \mathbin{\backslash\!\!\backslash} \pi$$.

Lemma: For all indices $$i$$ and $$j$$, the shortest paths $$\textsf{path}_P(a_i,b_j)$$ and $$\textsf{path}_P(a_1,b_1)$$ cross at most once.

Let $$D’$$ be a duplicate of $$D$$, and let $$DD$$ denote the map obtained by gluing the “left” copy of $$\pi$$ in $$D$$ to the “right” copy of $$\pi$$ in $$D’$$.

Lemma: For all indices $$i$$ and $$j$$, we have $\textsf{dist}_P(a_i,b_j) = \min\big\{ \textsf{dist}_{DD}(a_i,b_j),~ \textsf{dist}_{DD}(a’_i,b_j),~ \textsf{dist}_{DD}(a_i,b’_j) \big\}$.

The matrix of distances from the “top edge” of $$DD$$ to the “bottom edge” of $$DD$$ is Monge! So we can relax all the edges from $$A$$ to $$B$$ in $$O(k+l)$$ time with one invocation of SMAWK.

Suppose $$P$$ has $$h = O(1)$$ holes, with $$k_1, k_2, \dots, k_h$$ boundary vertices; recall that $$k = \sum_i k_i = O(\sqrt{r})$$. Then the total time to relax all boundary-to-boundary edges in $$P$$ is $\sum_{i=1}^h O(k_i\alpha(k_i)) + \sum_{i=1}^h \sum_{j=1}^h O(k_i + k_j) ~=~ O(k\alpha(k) + kh) ~=~ O(\sqrt{r}\,\alpha(r)).$ It follows that the time to relax all edges of $$\mathcal{R}$$ is $$O((n/\sqrt{r})\alpha(r)) = o(n)$$

## Faster Minimum Cuts: FR-Dijkstra

Frederickson held the record for fastest planar minimum-cut algorithm for almost two and a half decades; the record was finally broken in 2011 by two independent pairs of researchers, who ultimately published their result jointly: Giuseppe Italiano, Yahav Nussbaum, Piotr Sankowski, and Christian Wulff-Nilsen. Their $$O(n\log\log n)$$-time algorithm relies on an improvement to Dijkstra’s algorithm in dense distance graphs, proposed by Jittat Fakcharoenphol and Satish Rao in 2001, and now usually called FR-Dijkstra.

Recall from the previous lecture on shortest paths that we can compute a dense distance graph for a nice $$r$$-division in $$O(n\log r)$$ time. The dense distance graph has $$O(n/\sqrt{r})$$ vertices—the boundary vertices of the pieces of the $$r$$-division—and $$O(n)$$ edges. So Dijkstra’s algorithm with a Fibonacci heap runs in $$O(E + V\log V) = O(n + (n/\sqrt{r})\log n)$$ time.

FR-Dijkstra removes the $$O(n)$$ term from this running time. Specifically, within each piece of the $$r$$-division, the algorithm exploits the Monge structure in the boundary-to-boundary distances to avoid looking at every pair of boundary vertices. This is the same high-level strategy that we previously used with FR-Bellman-Ford, but with one significant difference: We do not know the relevant Monge arrays in advance. Instead, portions of each Monge array are revealed each time the Dijkstra wavefront touches the corresponding piece of the $$r$$-division.

I’ll discuss FR-Dijkstra in detail, along with the faster planar minimum-cut algorithm, in the next lecture.

## References

1. Greg N. Frederickson. Fast algorithms for shortest paths in planar graphs with applications. SIAM J. Comput. 16(8):1004–1004, 1987.

2. Alon Itai and Yossi Shiloach. Maximum flow in planar networks. SIAM J. Comput. 8:135–150, 1979.

3. Giuseppe F. Italiano, Yahav Nussbaum, Piotr Sankowski, and Christian Wulff-Nilsen. Improved algorithms for min cut and max flow in undirected planar graphs. Proc. 43rd Ann. ACM Symp. Theory Comput., 313–322, 2011.

4. Shay Mozes and Christian Wulff-Nilsen. Shortest paths in planar graphs with real lengths in $$O(n \log^2 n/\log\log n)$$ time. Proc. 18th Ann. Europ. Symp. Algorithms, 206–217, 2010. Lecture Notes Comput. Sci. 6347, Springer-Verlag. arXiv:0911.4963.

5. John Reif. Minimum $$s$$-$$t$$ cut of a planar undirected network in $$O(n\log^2 n)$$ time. SIAM J. Comput. 12(1):71–81, 1983.

6. Hassler Whitney. Non-separable and planar graphs. Trans. Amer. Math. Soc. 34:339–362, 1932.

7. Hassler Whitney. Planar graphs. Fund. Math. 21:73–84, 1933.

## Aptly Named Sir Not

• Maximum cuts (or minimum cuts with negative capacities) [Hadlock 1975]
• Deriving maximum flows from minimum cuts [Hassin Johnson 1985]
• $$O(n)$$ time for unweighted graphs [Weihe 1997, Eisenstat Klein 2013, Balzotti Franciosa 2022]
• Global minimum cuts (dual to shortest weighted cycle) [Łącki Sankowski 2011]
• Minimum cuts in directed planar graphs, via double cover [Janiga Koubek 1992 but fixed, Erickson 2011, Fox 2013]
• Faster directed planar minimum cuts [Mozes Nikolaev Nussbaum Weimann SODA 2018]

1. Without the second base case, it is possible for a constant fraction of the vertices to appear in a constant fraction of the recursive subproblems, leading to a running time of $$O(kn\log n)$$.↩