text:
# Comparing Python Script Performance Across Platforms

## Wrangling cProfile Output

Organizing the results from Python’s profiling module helps compare them across disparate platforms, especially when the code path is different between them.

I ran into a case where a small script was taking much longer on a particular operating system than others. Using Python’s built-in cProfile module, I was easily able to create a table of the functions that were taking the most time:

>python -O -mcProfile -scumtime load-all.py
330729 function calls (324916 primitive calls) in 0.680 seconds

Ordered by: cumulative time

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
236/1    0.020    0.000    0.681    0.681 {built-in method builtins.exec}
1    0.000    0.000    0.681    0.681 load-all.py:1(<module>)
278/3    0.002    0.000    0.457    0.152 <frozen importlib._bootstrap>:966(_find_and_load)
278/3    0.001    0.000    0.457    0.152 <frozen importlib._bootstrap>:936(_find_and_load_unlocked)
238/3    0.001    0.000    0.456    0.152 <frozen importlib._bootstrap>:651(_load_unlocked)
...

However, the relative times spent in various functions were not the same across the platforms I was testing. When the distribution of cumulative time in functions is different, the comparison becomes complicated and a bit of reorganization helps to identify problem areas.

## Gathering Data: Profiling the Script on Multiple Platforms

The over-all goal of this work was to improve the beginning to end run time of a script. Therefore, import time and all set-up and tear-down work was important and to be measured. Luckily, this made it extremely easy to profile right from the command line with a simple:

> python -O -mcProfile -oload-data.linux.prof load-data.py

So the output the cProfile run on the script load-data.py was output to a file called load-data.linux.prof. Great! I ran a similar command on Windows to get a file I called load-data.win.prof. Noticed I used the -O option when running this to avoid any if __debug__ or assert statements since these are really only useful for debugging and should not be relied upon in production code. In fact, I’m going to go ahead and claim that normal operation of any Python script should be done with the -O option, dropping it only when debugging, altering or writing code.

One thing we’ll have to do is to strip the leading directories on all function identifiers on the system where the profiling was done. This is because the pstats.Stats.strip_dirs() method does not currently handle windows paths when analyzed on a Linux machine. This can be done using the interpreter fairly easily:

>>> import pstats

and the profiling data should be overwritten without the leading paths in the filenames.

## Initial Comparison

The first pass in comparing the timings measured by cProfile is easy: just list them side by side. Using the pstats Python module, I used this script:

import os, pstats, sys

class ProfileResults(pstats.Stats):
def __init__(self, filename):
super().__init__(filename)
self.name = os.path.splitext(os.path.basename(filename))[0]
self = self.strip_dirs()

results = [ProfileResults(f) for f in sys.argv[1:]]
results[0].sort_stats('cumtime').print_stats(20)

to print the timings sorted by cumulative time spent in each function. If this script seems overly complicated, don’t worry about it quite yet - recall we are setting up for the comparison between two profiling results. Anyways, here are my results for Linux:

load-all.linux.prof

527040 function calls (521169 primitive calls) in 0.812 seconds

Ordered by: cumulative time
List reduced from 1954 to 20 due to restriction <20>

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
236/1    0.018    0.000    0.815    0.815 {built-in method builtins.exec}
1    0.000    0.000    0.815    0.815 load-all.py:1(<module>)
278/3    0.002    0.000    0.427    0.142 <frozen importlib._bootstrap>:966(_find_and_load)
278/3    0.001    0.000    0.427    0.142 <frozen importlib._bootstrap>:936(_find_and_load_unlocked)
238/3    0.001    0.000    0.426    0.142 <frozen importlib._bootstrap>:651(_load_unlocked)
193/3    0.001    0.000    0.426    0.142 <frozen importlib._bootstrap_external>:672(exec_module)
339/3    0.000    0.000    0.425    0.142 <frozen importlib._bootstrap>:211(_call_with_frames_removed)
896/40    0.001    0.000    0.418    0.010 <frozen importlib._bootstrap>:997(_handle_fromlist)
96/7    0.000    0.000    0.418    0.060 {built-in method builtins.__import__}
50/20    0.000    0.000    0.371    0.019 contextlib.py:49(inner)
1    0.000    0.000    0.257    0.257 session.py:207(connect)
1    0.000    0.000    0.257    0.257 tecutil_connector.py:124(connect)
1    0.000    0.000    0.240    0.240 tecutil_client.py:2(<module>)
17    0.001    0.000    0.228    0.013 __init__.py:1(<module>)
150    0.002    0.000    0.221    0.001 preamble.py:20(_fn)
150    0.204    0.001    0.216    0.001 tecutil_client.py:362(sndrcv)
1    0.000    0.000    0.210    0.210 message_pb2.py:4(<module>)
1    0.000    0.000    0.194    0.194 __init__.py:33(<module>)

and for Windows:

load-all.win.prof

290648 function calls (283683 primitive calls) in 1.947 seconds

Ordered by: cumulative time
List reduced from 2044 to 20 due to restriction <20>

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
245/1    0.003    0.000    1.948    1.948 {built-in method builtins.exec}
1    0.000    0.000    1.948    1.948 load-all.py:1(<module>)
50/20    0.000    0.000    1.332    0.067 contextlib.py:71(inner)
130    0.001    0.000    1.330    0.010 preamble.py:20(_fn)
130    1.295    0.010    1.328    0.010 tecutil_client.py:362(sndrcv)
306/3    0.003    0.000    0.591    0.197 <frozen importlib._bootstrap>:978(_find_and_load)
306/3    0.001    0.000    0.591    0.197 <frozen importlib._bootstrap>:948(_find_and_load_unlocked)
244/3    0.001    0.000    0.587    0.196 <frozen importlib._bootstrap>:663(_load_unlocked)
202/3    0.001    0.000    0.587    0.196 <frozen importlib._bootstrap_external>:722(exec_module)
380/3    0.000    0.000    0.586    0.195 <frozen importlib._bootstrap>:211(_call_with_frames_removed)
832/45    0.001    0.000    0.573    0.013 <frozen importlib._bootstrap>:1009(_handle_fromlist)
117/7    0.001    0.000    0.572    0.082 {built-in method builtins.__import__}
1    0.000    0.000    0.378    0.378 session.py:207(connect)
1    0.000    0.000    0.378    0.378 tecutil_connector.py:124(connect)
1    0.000    0.000    0.355    0.355 tecutil_client.py:2(<module>)
10    0.000    0.000    0.345    0.034 layout.py:332(active_frame)
16    0.000    0.000    0.339    0.021 __init__.py:1(<module>)
1    0.000    0.000    0.290    0.290 message_pb2.py:4(<module>)

This script calls into native libraries quite a bit and so the tottime is mostly time spent in compiled C++ code. If we can identify this as the culprit, then we could fire up valgrind or similar to start profiling that. But since we are really interested in the beginning-to-end run time of the script, we sort first on cumulative time as shown above.

At this point, we can take this information and start working on optimizing the longest running methods and work our way down. But let’s say we’ve done quite a bit of this already and we now want to understand the differences when running on these two platforms. One reason for wanting to do this would be if a script runs many times slower on a specific system and we want to understand why.

## Reorganization for a Better Understanding

The first thing to notice is that the script takes about twice as long on Windows as it does on Linux. The computer hardware was similar but not exactly so and we should try to normalize the times across the platforms tested. The only normalization factor that makes sense at this point in the analysis is the total time as shown in the results above: 0.812 s on Linux and 1.947 s on Windows. Adding a normalization method to our script effectively transforms all times to fractions of total execution time:

import os, pstats, sys

class ProfileResults(pstats.Stats):
def __init__(self, filename):
super().__init__(filename)
self.name = os.path.splitext(os.path.basename(filename))[0]
self = self.strip_dirs()
self.normalize_times()

def normalize_times(self):
total_tt = self.total_tt
oldstats = self.stats
self.stats = {}
for func, (cc, nc, tt, ct, oldcallers) in oldstats.items():
tt /= total_tt
ct /= total_tt
callers = {}
for func2, (cc2, nc2, tt2, ct2) in oldcallers.items():
tt2 /= total_tt
ct2 /= total_tt
callers[func2] = (cc2, nc2, tt2, ct2)
self.stats[func] = (cc, nc, tt, ct, callers)

results = [ProfileResults(f) for f in sys.argv[1:]]
results[0].sort_stats('cumtime').print_stats(20)

And Now our results look like this for Linux:

load-all.linux.prof

527040 function calls (521169 primitive calls) in 0.812 seconds

Ordered by: cumulative time
List reduced from 1954 to 20 due to restriction <20>

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
236/1    0.022    0.000    1.003    1.003 {built-in method builtins.exec}
1    0.000    0.000    1.003    1.003 load-all.py:1(<module>)
278/3    0.002    0.000    0.525    0.175 <frozen importlib._bootstrap>:966(_find_and_load)
278/3    0.002    0.000    0.525    0.175 <frozen importlib._bootstrap>:936(_find_and_load_unlocked)
238/3    0.002    0.000    0.524    0.175 <frozen importlib._bootstrap>:651(_load_unlocked)
193/3    0.001    0.000    0.524    0.175 <frozen importlib._bootstrap_external>:672(exec_module)
339/3    0.000    0.000    0.524    0.175 <frozen importlib._bootstrap>:211(_call_with_frames_removed)
896/40    0.001    0.000    0.515    0.013 <frozen importlib._bootstrap>:997(_handle_fromlist)
96/7    0.000    0.000    0.515    0.074 {built-in method builtins.__import__}
50/20    0.001    0.000    0.457    0.023 contextlib.py:49(inner)
1    0.000    0.000    0.316    0.316 session.py:207(connect)
1    0.000    0.000    0.316    0.316 tecutil_connector.py:124(connect)
1    0.000    0.000    0.296    0.296 tecutil_client.py:2(<module>)
17    0.001    0.000    0.281    0.017 __init__.py:1(<module>)
150    0.002    0.000    0.272    0.002 preamble.py:20(_fn)
150    0.251    0.002    0.266    0.002 tecutil_client.py:362(sndrcv)
1    0.000    0.000    0.259    0.259 message_pb2.py:4(<module>)
1    0.000    0.000    0.239    0.239 __init__.py:33(<module>)

and for Windows:

load-all.win.prof

290648 function calls (283683 primitive calls) in 1.947 seconds

Ordered by: cumulative time
List reduced from 2044 to 20 due to restriction <20>

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
245/1    0.002    0.000    1.000    1.000 {built-in method builtins.exec}
1    0.000    0.000    1.000    1.000 load-all.py:1(<module>)
50/20    0.000    0.000    0.684    0.034 contextlib.py:71(inner)
130    0.000    0.000    0.683    0.005 preamble.py:20(_fn)
130    0.665    0.005    0.682    0.005 tecutil_client.py:362(sndrcv)
306/3    0.001    0.000    0.303    0.101 <frozen importlib._bootstrap>:978(_find_and_load)
306/3    0.001    0.000    0.303    0.101 <frozen importlib._bootstrap>:948(_find_and_load_unlocked)
244/3    0.001    0.000    0.302    0.101 <frozen importlib._bootstrap>:663(_load_unlocked)
202/3    0.001    0.000    0.302    0.101 <frozen importlib._bootstrap_external>:722(exec_module)
380/3    0.000    0.000    0.301    0.100 <frozen importlib._bootstrap>:211(_call_with_frames_removed)
832/45    0.000    0.000    0.294    0.007 <frozen importlib._bootstrap>:1009(_handle_fromlist)
117/7    0.000    0.000    0.294    0.042 {built-in method builtins.__import__}
1    0.000    0.000    0.194    0.194 session.py:207(connect)
1    0.000    0.000    0.194    0.194 tecutil_connector.py:124(connect)
1    0.000    0.000    0.182    0.182 tecutil_client.py:2(<module>)
10    0.000    0.000    0.177    0.018 layout.py:332(active_frame)
16    0.000    0.000    0.174    0.011 __init__.py:1(<module>)
1    0.000    0.000    0.149    0.149 message_pb2.py:4(<module>)

From here, we can finally start to compare directly the results of the profiling and first off, we notice that importing Python modules seems to be taking more time on Linux – recall that from now on, when we say “time spent,” it is always in fractions of total execution time. Read another way, we can say that the operations in and around the tecutil_client.py:362(sndrcv) method are taking longer on Windows than on Linux.

After some thought, I realized what I wanted to know was the relative difference in times for the various functions. So, focusing on cumulative time, we now print out the times for each platform along with the difference between them, sorting on the absolute magnitude of this difference. Here is the new comparison script:

import math, os, pstats, sys

class ProfileResults(pstats.Stats):
def __init__(self, filename):
super().__init__(filename)
self.name = os.path.splitext(os.path.basename(filename))[0]
self = self.strip_dirs()
self.normalize_times()

def normalize_times(self):
total_tt = self.total_tt
oldstats = self.stats
self.stats = {}
for func, (cc, nc, tt, ct, oldcallers) in oldstats.items():
tt /= total_tt
ct /= total_tt
callers = {}
for func2, (cc2, nc2, tt2, ct2) in oldcallers.items():
tt2 /= total_tt
ct2 /= total_tt
callers[func2] = (cc2, nc2, tt2, ct2)
self.stats[func] = (cc, nc, tt, ct, callers)

def compare_cumtime(self, other):
res = []
for func, (_, _, _, ct0, _) in self.stats.items():
try:
ct1 = other.stats[func][3]
except KeyError:
ct1 = math.nan
res.append(('{}:{}({})'.format(*func), ct0, ct1))
return res

results = [ProfileResults(f) for f in sys.argv[1:]]
ctcomp = results[0].compare_cumtime(results[1])
ctcomp = sorted(ctcomp, key=lambda x: -abs(x[2] - x[1]))

print('''cumulative time:
{}  {}   diff   function'''.format(*[r.name for r in results]))
for func, ct0, ct1 in ctcomp[:20]:
print(' {0: ^{1}.3f}  {2: ^{3}.3f}  {4: >+6.3f}  {5}'.format(
ct0, len(results[0].name), ct1, len(results[1].name),
ct1 - ct0, func))

and here are the results:

cumulative time:
0.012          0.103      +0.091  tecutil_rpc.py:1637(DataSetIsAvailableForFrame)
0.010          0.088      +0.078  tecutil_rpc.py:2431(FrameGetActiveID)
0.012          0.088      +0.077  tecutil_rpc.py:3898(PageGetUniqueID)
0.005          0.044      +0.040  tecutil_rpc.py:2443(FrameGetDataSetUniqueIDByFrameID)
0.005          0.045      +0.040  tecutil_rpc.py:5142(StringListAlloc)
0.005          0.044      +0.039  tecutil_rpc.py:2439(FrameGetCount)
0.006          0.044      +0.038  tecutil_rpc.py:5167(StringListDealloc)
0.006          0.044      +0.038  tecutil_rpc.py:5152(StringListAppendString)
0.096          0.074      -0.022  tecutil_rpc.py:3845(NewLayout)
0.001          0.000      -0.001  tecutil_rpc.py:2(<module>)
0.000          0.000      -0.000  tecutil_rpc.py:21(TecUtilRPC)
0.006           nan         +nan  tecutil_rpc.py:5177(StringListGetCount)
0.005           nan         +nan  tecutil_rpc.py:5182(StringListGetRawStringPtr)
0.030          0.009      -0.021  pyparsing.py:2738(__init__)
0.033          0.014      -0.019  pyparsing.py:61(<module>)
0.014          0.006      -0.009  pyparsing.py:5353(pyparsing_common)
0.000          0.006      +0.006  reflection.py:46(<module>)
0.007          0.003      -0.004  pyparsing.py:4210(__init__)
0.007          0.002      -0.004  pyparsing.py:3702(leaveWhitespace)

Reading from top, down we see that, on Windows, the function DataSetIsAvailableForFrame is taking up almost an order of magnitude more time than on Linux along with several StringList related functions. Also, it looks like StringListGetCount is not called at all on Windows - this may be OK, but it would be best to understand why it’s only called on Linux. It may mean that the code makes some assumption that is good on one platform but not on the other, and the point is not these specific functions but rather that we now have a direction to go in the investigation of script performance.

## Generalization of Profile Comparisons and Conclusion

It may not always be possible to test the exact code against multiple platforms, especially if you do not have direct access to the platform you are debugging. However, if you can run similar code and you still want to compare the profiling results, you can make one more change to the script above to ensure a functions match between datasets, and that is: remove the line number in the function key. To do this, I set the line numbers to zero and then call pstats.Stats.strip_dirs() one more time. The final script looks like this:

import math, os, pstats, sys

class ProfileResults(pstats.Stats):
def __init__(self, filename):
super().__init__(filename)
self.name = os.path.splitext(os.path.basename(filename))[0]
self = self.strip_dirs()
self.normalize_times()
self = self.strip_dirs()

def normalize_times(self):
total_tt = self.total_tt
oldstats = self.stats
self.stats = {}
for func, (cc, nc, tt, ct, oldcallers) in oldstats.items():
func = (func[0], 0, func[2])
tt /= total_tt
ct /= total_tt
callers = {}
for func2, (cc2, nc2, tt2, ct2) in oldcallers.items():
func2 = (func2[0], 0, func2[2])
tt2 /= total_tt
ct2 /= total_tt
callers[func2] = (cc2, nc2, tt2, ct2)
self.stats[func] = (cc, nc, tt, ct, callers)

def compare_cumtime(self, other):
res = []
for func, (_, _, _, ct0, _) in self.stats.items():
try:
ct1 = other.stats[func][3]
except KeyError:
ct1 = math.nan
res.append(('{0}({2})'.format(*func), ct0, ct1))
return res

results = [ProfileResults(f) for f in sys.argv[1:]]
ctcomp = results[0].compare_cumtime(results[1])
ctcomp = sorted(ctcomp, key=lambda x: -abs(x[2] - x[1]))

print('''cumulative time:
{}  {}   diff   function'''.format(*[r.name for r in results]))
for func, ct0, ct1 in ctcomp[:20]:
print(' {0: ^{1}.3f}  {2: ^{3}.3f}  {4: >+6.3f}  {5}'.format(
ct0, len(results[0].name), ct1, len(results[1].name),
ct1 - ct0, func))

and the results from the data above now looks like this:

cumulative time:
0.012          0.103      +0.091  tecutil_rpc.py(DataSetIsAvailableForFrame)
0.010          0.088      +0.078  tecutil_rpc.py(FrameGetActiveID)
0.012          0.088      +0.077  tecutil_rpc.py(PageGetUniqueID)
0.005          0.044      +0.040  tecutil_rpc.py(FrameGetDataSetUniqueIDByFrameID)
0.005          0.045      +0.040  tecutil_rpc.py(StringListAlloc)
0.005          0.044      +0.039  tecutil_rpc.py(FrameGetCount)
0.006          0.044      +0.038  tecutil_rpc.py(StringListDealloc)
0.006          0.044      +0.038  tecutil_rpc.py(StringListAppendString)
0.096          0.074      -0.022  tecutil_rpc.py(NewLayout)
0.001          0.000      -0.001  tecutil_rpc.py(<module>)
0.000          0.000      -0.000  tecutil_rpc.py(TecUtilRPC)
0.006           nan         +nan  tecutil_rpc.py(StringListGetCount)
0.005           nan         +nan  tecutil_rpc.py(StringListGetRawStringPtr)
0.062          0.043      -0.019  requirements.py(<module>)
0.033          0.014      -0.019  pyparsing.py(<module>)
0.018          0.001      -0.017  version.py(__init__)
0.018          0.002      -0.016  version.py(parse)
0.022          0.006      -0.016  specifiers.py(<module>)
0.019          0.005      -0.014  specifiers.py(Specifier)
0.014          0.006      -0.009  pyparsing.py(pyparsing_common)

So while the first actions taken to address performance of any code is to go after the slowest or longest-running sections, it may also help to evaluate and compare what these are across several of the target platforms.


text:
# What are QML algorithms using less than 8 qubits and provide a quantum advantage?

So this is more of a soft question.

I've been trying to find some quantum machine learning algorithms can both be run with less than 8 qubits and provide a quantum advantage to classical machine learning algorithms. However, I have not been able to find any.

Could the readers of this please tell me the names of some algorithms that satisfy the two qualities that you know of, if there are any?

Thank you.

Note that 8 qubits is not enough to show any quantum advantage. For such number of qubits, you have only $$2^8=256$$ basis states. Multiplying matrices and vectors of that size can be carried out in reasonable time on a classical computer.

format: plain


text:
# $7^x-3^y=4$ in the integers

I have to find all solutions to $$7^x-3^y=4$$ in the integers, I've already proven that $$x$$ and $$y$$ have the same parity and that they cannot be even. But I'm stuck in the case when $$x$$ and $$y$$ are odd. Could someone show me how to solve for all solutions? I know there is at least one, $$(x,y)=(1,1)$$, but I don't know if there's any other solution

• I remember this as an Indian National Math Olympiad problem – N.S.JOHN Oct 20 '18 at 20:08
• I found out that $y$ must be of the form $6k+1$. Don't know if it will be of any help – Jakobian Oct 20 '18 at 20:24

## 2 Answers

If $$y\leq 1$$, the only solution is $$(1,1)$$, which you found. If $$y\geq 2$$, then consider the equation $$\bmod 9$$. We see

$$7^x\equiv 4\bmod 9$$

We have

$$7^0\equiv 1,\ 7^1\equiv 7,\ 7^2\equiv 4,\ 7^3\equiv 1,$$

so $$7^x\equiv 4\bmod 9$$ iff $$x\equiv 2\bmod 3$$. As a result, we have $$x\equiv 5\bmod 6$$ (as you have already shown $$x$$ is odd). So, as

$$7^{12}\equiv 1\bmod 13,$$

we have

$$7^x\equiv 7^5\mathrm{\ or\ }7^{11}\equiv \pm 2\bmod 13$$

since $$x\equiv 5\bmod 6\implies x\equiv 5\mathrm{\ or\ }11\bmod 12$$. This means that

$$3^y=7^x-4\equiv (\pm 2)-4\in\{7,11\}\bmod 13.$$

However

$$3^0\equiv 1,3^1\equiv 3,3^2\equiv 9,3^3\equiv 1\bmod 13,$$

so we get a contradiction.

• Looks good to me, +1 – Jakobian Oct 20 '18 at 21:38
• I didn't quite get why $7^x\equiv7^5$ and why $3^y\in\lbrace7,11\rbrace\mod13$ – Bruno Andrades Oct 20 '18 at 23:27
• @BrunoAndrades I've added a bit more explanation to my post. – Carl Schildkraut Oct 20 '18 at 23:43

We first show that if x and y have a common divisor like c, then a number such as $$a^x-b^y$$ has a factor like $$a-b$$(it does not mean if x and y have no common divisor $$a^x-b^y$$ definitely has no common divisor like $$a-b$$):

We know that:

$$(a^x-1, a^y-1)=a^c-1$$

$$(b^x-1, b^y-1)=b^c-1$$

Now take one say $$a^x-1$$ which has a factor like $$a^c-1$$ and $$b^y-1$$ which has a factor like $$b^c-1$$ and subtract them:

$$gcd(a^x-b^y)≡a^c-b^c=(a-b)(a^{c-1}+a^{c-2}b+ . . .)$$

So we may write:

$$a^x-b^y=k(a-b)$$

It is clear that for a certain set of values for a, b,x and y there exist only one k to satisfy this relation. Hence there are infinitely many such unique relation for certain values of involved parameters. That means equation $$7^x-3^y=4$$ in which $$a=7,. b=3,.$$ and $$k=1$$ are certain values has only one solution $$x=y=1$$.


text:
# LaplaceTransform vs Integrate with inconsistent result

I am trying to evaluate the following integral, which happens to be a laplace transform

$$G(x,x';z)=\frac{1}{i\hslash}\int_0^{\infty}dt \exp\left(\frac{izt}{\hslash}\right)\sqrt{\frac{m}{2\pi i \hslash t}}\exp\left(-\frac{m(x-x')^2}{2i\hslash t}\right),$$ with $z=\frac{k^2\hslash^2}{2m}+i\epsilon$, where $\epsilon$ is positive, then taking the limit $\epsilon\rightarrow 0$.

I'm going to list my attempts and results. The answer to this integral should be $$\frac{m}{ik\hslash^2}\exp{i k Abs(x)}.$$

Attempt 1:

2/(I \[HBar]) Integrate[E^(I ((z) u^2)/\[HBar]) Sqrt[m/(2\[Pi] I \[HBar])] E^(I (x^2 m)/(2\[HBar] u^2))/.z->((k \[HBar])^2/(2m)),{u,0,\[Infinity]},Assumptions->{k>0,\[HBar]>0,m>0,x\[Element]Reals,x!=0,z\[Element]Reals}]

Here I made the variable tranformation $t\rightarrow u^2$

results:

If I set $z$ as real I get the correct answer.

If I set $z=\frac{k^2\hslash^2}{2m}+i\epsilon$

2/(I \[HBar]) Integrate[E^(I ((z) u^2)/\[HBar]) Sqrt[m/(2\[Pi] I \[HBar])] E^(I (x^2 m)/(2\[HBar] u^2))/.z->((k \[HBar])^2/(2m)+I \[Epsilon]),{u,0,\[Infinity]},Assumptions->{k>0,\[HBar]>0,m>0,x\[Element]Reals,x!= 0,\[Epsilon]>0}]

I get a cossine instead of the exponential. TrigToExp doesn't solve it, the answer is indeed incorrect.

Attempt 2(Laplace Transform then Full Simplify):

1/(I \[HBar]) LaplaceTransform[HeavisideTheta[t] Sqrt[m/(2\[Pi] I \[HBar] t)]Exp[(I x^2 m)/(2\[HBar] t)],t,-I(e+I \[Epsilon])/\[HBar]]

FullSimplify[%/.{e->(k^2 \[HBar]^2)/(2m),\[Epsilon]-> 0},{k>0,\[HBar]>0,m>0,x\[Element]Reals,x!= 0,z\[Element]Reals}]

Here I get the right result.

I wouldn't be worried if Mathematica couldn't evaluate the integrals above, but it seems to be doing so and giving the wrong output. What am I doing wrong?

• Note that setting z reals goes againt the method, z should definetly be complex. But when I set z real is when I get the right answer, which is very weird. – Pedro Portugal Jun 9 '17 at 19:01


text:
School has started up, and I’m in a class called Financial Computing. I thought it might be interesting to share some of my assignments and explain what I learn along the way. Most of the posts won’t be describing the Stochastic Calculus involved in each assignment, but will instead focus on the details of the implementation in R. I don’t claim to have the best way for any of these assignments, but perhaps you can learn something!

### What are we doing today?

In this post, I’ll work through simulating paths of a stock that follows the log normal distribution used in the Black Scholes model. Importantly, I’ll explain my thought process as I tried to optimize the implementation from loops to vectors.

As an added bonus, at the very bottom is some extra content on a basic function that I have created to replicate the concept of broadcasting from Python. Someone could (should?) probably create an entire package out of this idea.

### The Model

Under the Black Scholes model, the stock price, $S_t$, at time t follows a Geometric Brownian Motion, which means that it satisfies the Stochastic Differential Equation:

$dS_t = r S_t dt + \sigma S_t dW_t$

Where:

• $r =$ Drift - The average log return per time unit
• $\sigma =$ Volatility - How much variance is in that drift
• $W_t =$ Brownian Motion - Random noise from a normal distribution with mean 0 and variance t

Interestingly, we actually have the solution to this equation (one of the few we have analytical solutions for):

$S_t = S_0 \times e^{(r - \frac{1}{2} \sigma^2) t + \sigma W_t}$

More generally, this can be written as a formula providing us with the recursive equation:

$S_{t_i} = S_{t_{i-1}} \times e^{(r - \frac{1}{2} \sigma^2) ({t_i - t_{i-1}}) + \sigma (W_{t_i} - W_{t_{i-1}}) }$

If you want to know how to get the solution, this is a pretty good explanation, but be prepared to learn about Ito’s Lemma.

### Really, dude?

Okay, that’s a lot without any explanation, and I get that. But the point of this post is more to explain how to simulate paths of $S_t$. So how do we do that?

• We will start from time 0 with an initial stock price, then we will generate the next stock price from that using the recursive formula, and so on.
• The only random piece is the brownian motion increment (dW), which we will generate at each time point using draws from a normal distribution.
• The time increment will be a constant value (dt) to keep things simple.

I was given some starting parameters:

# Set a seed too so I can reproduce the exact numbers
set.seed(123)

# Parameters
r       <- 0.028
sigma   <- 0.255
time_T  <- 0.5
dt      <- 1/12
t_total <- time_T / dt
s_0     <- 100
N       <- 10000

Where time_T and dt mean each simulation path will go from time 0 to time_T by increments of dt. Dividing time_T / dt gets us the total number of time steps required. N is the number of paths to simulate.

### First attempt

We all know loops are to be avoided when you can in R, and that you should instead vectorize the operations. At first, I thought this wasn’t going to be possible, as this is a recursive type of formula where the next value relies on the previous one. With that in mind, I created the following implementation.

First off, set up a matrix to fill with the 10000 simulations (one per row), each one having 6 time steps (7 columns total including the initial stock price).

# Create a 10000x7 matrix of NA's to fill in
# Each row is a simulation
s <- matrix(NA_real_, nrow = N, ncol = t_total+1)

# The first column is just the initial price
s[,1] <- s_0

head(s)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]  100   NA   NA   NA   NA   NA   NA
## [2,]  100   NA   NA   NA   NA   NA   NA
## [3,]  100   NA   NA   NA   NA   NA   NA
## [4,]  100   NA   NA   NA   NA   NA   NA
## [5,]  100   NA   NA   NA   NA   NA   NA
## [6,]  100   NA   NA   NA   NA   NA   NA

So far so good, now let’s create a matrix for dW, our brownian motion increments. A very important fact is that these are all independent of each other, so the generation of them is straightforward. Each increment:

$W_{t_i} - W_{t_{i-1}}$

can be drawn from a normal distribution with mean 0 and variance $t_i - t_{i-1}$ (which is what I have defined as dt because it is a constant).

# ~N(0, dt)
# To fill in 10000 simulations, and move forward 6 steps, we need a 10000x6 matrix
dW <- rnorm(N * t_total, mean = 0, sd = sqrt(dt))
dW <- matrix(dW, N, t_total)
head(dW)
##             [,1]        [,2]        [,3]        [,4]        [,5]
## [1,] -0.16179538  0.68436942 -0.24141807 -0.05588936  0.13929365
## [2,] -0.06644652 -0.04815446 -0.06367394  0.07452070  0.20824004
## [3,]  0.44996033  0.26759071 -0.60723241 -0.15539745 -0.14658701
## [4,]  0.02035402 -0.16401128 -0.48145457 -0.34036612 -0.01867981
## [5,]  0.03732215  0.06497791 -0.31695458  0.25999452  0.37589000
## [6,]  0.49509662  0.32677618 -0.48082343 -0.00469084 -0.06311503
##             [,6]
## [1,]  0.07493710
## [2,]  0.26486253
## [3,] -0.20851535
## [4,] -0.23333150
## [5,] -0.04080481
## [6,]  0.65154366

Based on this setup, I thought I would need a loop. The algorithm would step through the 10000 simulations all at once, but would have to loop through the 6 time steps one at a time, because each time step depended on the previous one. So, following the formula (below again for reference) I did this:

$S_{t_i} = S_{t_{i-1}} \times e^{(r - \frac{1}{2} \sigma^2) ({t_i - t_{i-1}}) + \sigma (W_{t_i} - W_{t_{i-1}}) }$

for(i in 1:(t_total)) {
s[,i+1] <- s[,i] * exp((r - 1/2 * sigma^2) * dt + sigma * dW[,i])
}

head(s)
##      [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
## [1,]  100  95.92209 114.16839 107.31161 105.75330 109.53595 111.60722
## [2,]  100  98.28292  97.04695  95.44803  97.24258 102.50728 109.62854
## [3,]  100 112.11600 119.98823 102.73710  98.70849  95.05115  90.09528
## [4,]  100 100.48258  96.33055  85.16909  78.05933  77.65918  73.14576
## [5,]  100 100.91830 102.56581  94.56668 101.01083 111.13033 109.93864
## [6,]  100 113.41388 123.22299 108.96314 108.79196 107.01479 126.30943

And that does work! But can we avoid the loop? YES WE CAN!

### MATH

To avoid the loop, we are going to manipulate a simple case, and then apply it generally. One key element for this vectorization is that brownian motion increments are independent. Let’s think about what actually happens from time 0->1 and from 1->2.

$S_1 = S_0 \times e^{(r - \frac{1}{2} \sigma^2) (t_1 - t_0) + \sigma (W_1 - W_0)}$

$S_2 = S_1 \times e^{(r - \frac{1}{2} \sigma^2) (t_2 - t_1) + \sigma (W_2 - W_1)}$

If we plug the equation for S_1 into the equation for S_2

$S_2 = (S_0 \times e^{(r - \frac{1}{2} \sigma^2) (t_1 - t_0) + \sigma (W_1 - W_0)}) \times e^{(r - \frac{1}{2} \sigma^2) (t_2 - t_1) + \sigma (W_2 - W_1)}$

And then combine exponents…

$S_2 = (S_0 \times e^{(r - \frac{1}{2} \sigma^2) (t_1 - t_0 + t_2 - t_1) + \sigma (W_1 - W_0 + W_2 - W_1)})$

Notice that some of the t and W terms cancel:

$S_2 = (S_0 \times e^{(r - \frac{1}{2} \sigma^2) (t_2 - t_0) + \sigma (W_2 - W_0)})$

And by definition t_0 and W_0 are 0:

$S_2 = S_0 \times e^{(r - \frac{1}{2} \sigma^2) t_2 + \sigma W_2}$

This is actually the form that was proposed as the solution to the geometric brownian motion stochastic differential equation:

$S_t = S_0 \times e^{(r - \frac{1}{2} \sigma^2) t + \sigma W_t}$

It looks like we can actually generate S_2 without needing to know S_1 at all. Notice that the exponent now contains t_2 and W_2. t_2 is known beforehand, but W_2 seems like it would rely on W_1 in a way that has to be recursively calculated. Actually, if we think of W_2 as a sum of brownian motion increments (I told you this would help):

$W_2 = (W_2 - W_1) + (W_1 - W_0) = dW_2 + dW_1$

then W_2 is just the cumulative sum of the increments, and, by definition, each increment is independent of the previous increment so we can generate them all before hand (we already did this when we created the dW matrix).

# Rowwise cumulative sum of dW generates W1, W2, W3, ... for each simulation
W  <- plyr::aaply(dW, 1, cumsum)
head(W)
##
## X1            1          2          3           4          5          6
##   1 -0.16179538  0.5225740  0.2811560  0.22526661  0.3645603  0.4394974
##   2 -0.06644652 -0.1146010 -0.1782749 -0.10375422  0.1044858  0.3693483
##   3  0.44996033  0.7175510  0.1103186 -0.04507883 -0.1916658 -0.4001812
##   4  0.02035402 -0.1436573 -0.6251118 -0.96547795 -0.9841578 -1.2174893
##   5  0.03732215  0.1023001 -0.2146545  0.04534000  0.4212300  0.3804252
##   6  0.49509662  0.8218728  0.3410494  0.33635853  0.2732435  0.9247872

Unlike the recursive formula from before where dt was used, the time that we are currently at, t, is used instead so we will need that as well.

time_steps <- matrix(seq(from = dt, to = time_T, by = dt),
nrow = N, ncol = t_total, byrow = TRUE)
head(time_steps)
##            [,1]      [,2] [,3]      [,4]      [,5] [,6]
## [1,] 0.08333333 0.1666667 0.25 0.3333333 0.4166667  0.5
## [2,] 0.08333333 0.1666667 0.25 0.3333333 0.4166667  0.5
## [3,] 0.08333333 0.1666667 0.25 0.3333333 0.4166667  0.5
## [4,] 0.08333333 0.1666667 0.25 0.3333333 0.4166667  0.5
## [5,] 0.08333333 0.1666667 0.25 0.3333333 0.4166667  0.5
## [6,] 0.08333333 0.1666667 0.25 0.3333333 0.4166667  0.5

Now it’s a vectorized one-liner to calculate the stock price at each time!

# Stock price simulation
s_t <- s_0 * exp((r - 1/2 * sigma^2) * time_steps + sigma * W)

# Add the original stock price onto the front
s_t <- cbind(s_0, s_t)

# Add 0 as the column name for initial value (it's important I promise)
colnames(s_t)[1] <- "0"

head(s_t)
##     0         1         2         3         4         5         6
## 1 100  95.92209 114.16839 107.31161 105.75330 109.53595 111.60722
## 2 100  98.28292  97.04695  95.44803  97.24258 102.50728 109.62854
## 3 100 112.11600 119.98823 102.73710  98.70849  95.05115  90.09528
## 4 100 100.48258  96.33055  85.16909  78.05933  77.65918  73.14576
## 5 100 100.91830 102.56581  94.56668 101.01083 111.13033 109.93864
## 6 100 113.41388 123.22299 108.96314 108.79196 107.01479 126.30943

Just as a sanity check, this should have produced the same results as the for loop

# ignore the dimname attributes
all.equal(s, s_t, check.attributes = FALSE)
## [1] TRUE

### Now what?

There are a number of interesting things we could do with these results. One is to calculate the fair price of a European Option on this stock. I think I’ll save that for the next post.

Something else we might do is visualize the distribution of $S_T$, the stock price at the terminal (final) time. Because the stock price is modeled as an exponential of a normal random variable (W_t), the stock price itself has a log-normal distribution. For practicality, this means that it is right tailed and can’t drop below 0 (good properties of a stock).

library(tidyverse)

# let's just take a moment and admire the fact that I can put LaTeX in ggplots
library(latex2exp)

tibble::tibble(s_T = s_t[,6]) %>%
ggplot(mapping = aes(x=s_T)) +
geom_histogram(bins = 500) +
labs(x = TeX('S_T'), y = NULL, title = TeX('Log-normal Distribution of S_T') )

We could also look at the 6-step path of 100 of our simulations.

as_tibble(s_t) %>%
rownames_to_column(var = "iteration") %>%
gather(time_step, stock_price, -iteration) %>%
mutate(time_step = as.numeric(time_step),
iteration = as.factor(iteration)) %>%
filter(iteration %in% 1:100) %>%

ggplot(aes(x = time_step, y = stock_price, group = iteration)) +
geom_line() +
labs(x = "Time", y = "Stock Price", title = "100 Simulated Paths")

### Conclusion + Extra Content

If I haven’t bored you to tears yet, allow me to thank you for sticking around this long. I think these posts are useful because they force me to try and understand a concept a bit more than if I was just reading it from a book.

As you may have noted in the post above, I had to create a large matrix time_steps for R to perform the matrix addition I wanted correctly. I thought it would have been simple. Ideally I could create a 1x6 matrix of times, and add it to a 10000x6 matrix of the brownian motions and have the times matrix broadcasted to each row of the brownian motion matrix, adding element by element to each row. This works in Python and Matlab, but R has a mind of it’s own.

# First try with a vector
x1 <- c(1,2)
x1
## [1] 1 2
x2 <- matrix(c(1,2,3,4), nrow = 2, ncol = 2, byrow = TRUE)
x2
##      [,1] [,2]
## [1,]    1    2
## [2,]    3    4
# I want to add c(1,2) to each row of x2
x1 + x2
##      [,1] [,2]
## [1,]    2    3
## [2,]    5    6

That’s not right, it’s adding 1 to the first column of x2, then 2 to the second column of x2. To get what I want I could butcher it like this:

t(x1+t(x2))
##      [,1] [,2]
## [1,]    2    4
## [2,]    4    6

If x1 was a matrix instead of a vector, then it gives a non-conformable array error.

x1 <- matrix(x1, ncol = 2)
x1
##      [,1] [,2]
## [1,]    1    2
x1 + x2
## Error in x1 + x2: non-conformable arrays

So R is super strict here. That’s fine and all, but with other languages able to do this, and with it being such a natural way of thinking about this type of addition, I decided to roll my own function that allows me to add matrices together that meet certain conditions by broadcasting one of them over the other.

# Broadcast addition
# One of the special % % operators
%+% <- function(e1, e2) {

stopifnot(is.matrix(e1))
stopifnot(is.matrix(e2))

# e1 - e2 & 1 has more rows & equal cols
if(nrow(e1) >= nrow(e2) & ncol(e1) == ncol(e2)) {
case <- "1"

# e1 - e2 & 2 has more rows & equal cols
} else if(nrow(e1) < nrow(e2) & ncol(e1) == ncol(e2)) {
case <- "2"

# e1 - e2 & 1 has more cols & equal rows
} else if(ncol(e1) >= ncol(e2) & nrow(e1) == nrow(e2)) {
case <- "3"

# e1 - e2 & 2 has more cols & equal rows
} else if(ncol(e1) < ncol(e2) & nrow(e1) == nrow(e2)) {
case <- "4"

# Fail
} else {
stop("Incorrect dims")
}

switch(case,
"1" = t(apply(e1, 1, function(x) {x + e2})),
"2" = t(apply(e2, 1, function(x) {x + e1})),
"3" = t(apply(e1, 2, function(x) {x + e2})),
"4" = t(apply(e2, 2, function(x) {x + e1})))
}

Let’s see what this thing can do!

x1 %+% x2
##      [,1] [,2]
## [1,]    2    4
## [2,]    4    6

Nice! That’s what I want. One thing to note is that order of operations don’t work quite as you’d expect because of the precedence of the special %+% operator in relation to + and *, so you have to be really explicit.

# This tries to do addition first
x1 * 2 %+% x2
## Error: is.matrix(e1) is not TRUE
# Explicit parenthesis
(x1 * 2) %+% x2
##      [,1] [,2]
## [1,]    3    6
## [2,]    5    8

Armed with the ability to broadcast addition, let’s redo the last step of the stock price simulation.

# Instead of a massive matrix, just create a 1x6
time_steps <- matrix(seq(from = dt, to = time_T, by = dt), nrow = 1, ncol = t_total)
time_steps
##            [,1]      [,2] [,3]      [,4]      [,5] [,6]
## [1,] 0.08333333 0.1666667 0.25 0.3333333 0.4166667  0.5
# Remember that W is 10000x6
head(W)
##
## X1            1          2          3           4          5          6
##   1 -0.16179538  0.5225740  0.2811560  0.22526661  0.3645603  0.4394974
##   2 -0.06644652 -0.1146010 -0.1782749 -0.10375422  0.1044858  0.3693483
##   3  0.44996033  0.7175510  0.1103186 -0.04507883 -0.1916658 -0.4001812
##   4  0.02035402 -0.1436573 -0.6251118 -0.96547795 -0.9841578 -1.2174893
##   5  0.03732215  0.1023001 -0.2146545  0.04534000  0.4212300  0.3804252
##   6  0.49509662  0.8218728  0.3410494  0.33635853  0.2732435  0.9247872
# Add using broadcasted addition, making sure to be careful about parenthesis!
s_t <- s_0 * exp(((r - 1/2 * sigma^2) * time_steps) %+% (sigma * W))

s_t <- cbind(s_0, s_t, deparse.level = 0)
head(s_t)
##   [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
## 1  100  95.92209 114.16839 107.31161 105.75330 109.53595 111.60722
## 2  100  98.28292  97.04695  95.44803  97.24258 102.50728 109.62854
## 3  100 112.11600 119.98823 102.73710  98.70849  95.05115  90.09528
## 4  100 100.48258  96.33055  85.16909  78.05933  77.65918  73.14576
## 5  100 100.91830 102.56581  94.56668 101.01083 111.13033 109.93864
## 6  100 113.41388 123.22299 108.96314 108.79196 107.01479 126.30943

So much better! I have used this a few times in the past month or so. Credit to Alex Hayes for teaching me a bit about why broadcasting is awesome. I created the base for %+% in response to his comments here.


text:
# Review/rate my new graph class

I have written a PHP class called "graph". It is a class that performs RESTful-like commands to a MySQL database. I have posted the GitHub repo.

Here is the code as well:

config.php

<?php
//standard database connection data

//these obviously need to be changed....

define("DBHOST","localhost");
define("MAIN_DB","main");

?>

database.class.php

<?php
require_once('config.php');  //set your database settings in this

class database{

public static function dbConnect($db){$link = mysql_connect(DBHOST, DBUSERNAME, DBPASSWORD);

if (!$link) {$response = array(
"object"        => "database",
"instrcut"      => "dbConnect",
"status"        => "fail",
"message"       =>  mysql_error()

);

}

if($db == ''){ //sets the$db var to the constant MAIN_DB if not otherwise set.
$db = MAIN_DB;$db_select = mysql_select_db($db); }else{$db_select = mysql_select_db($db); } if(!$db_select){
$response = array( "object" => "database", "instrcut" => "dbConnect", "status" => "fail", "message" => mysql_error() ); }else{$response = array(
"object"        => "database",
"instrcut"      => "dbConnect",
"status"        => "pass",
"message"       =>  "The database '$db' was connected to successfully." ); } return$response;

}

public static function conQuer($sql,$db){

$query = self::dbConnect($db);
$result = mysql_query($sql);

if($query&&$result){
$response = array( "object" => "database", "instrcut" => "conQuer", "status" => "pass", "message" => "The connection was established and the query was run successfully.", "data" =>$result
);
}else{
$response = array( "object" => "database", "instrcut" => "conQuer", "status" => "fail", "message" => "The conQuer method failed because there was a bad query." ); } return$response;
}

//end database class
}

?>

Finally, the big one:

graph.class.php

 <?php
//PHP-GRAPH class

require_once('database.class.php');

class graph{
//used to retrieve existing data

public static function get($data){$jsonData = json_encode($data); //decode json data$data = json_decode($jsonData); //////////////////////////////$db = $data->db;$table = $data->table;$field = $data->field;$value = $data->value;$request = array($field,$value);
//////////////////////////////

if($db&&$table&&$field&&$value){
////////////GET ALL THE AVAILABLE FIELDS FROM SELECTED TABLE////////////

//query to get the available fields
$result = database::conQuer("SHOW COLUMNS FROM$table",$db);$result = $result['data']; if (!$result) {

//return the failure response
$response = array("obj"=>"graph","instruct"=>"get","status"=>"fail","msg"=>mysql_error()); } if (mysql_num_rows($result) > 0) {
while ($row = mysql_fetch_assoc($result)) {

$field_array[] =$row['Field'];

}
}
/////////////////////////////////////////////////////////////////////////////

////////////GET ALL THE DATA FROM THE SELECTED FIELDS ////////////

$getQuery = database::conQuer("SELECT * FROM$table WHERE $field = '$value'",$db);$getQuery = $getQuery['data'];$queryCount = mysql_numrows($getQuery); while($row = mysql_fetch_array($getQuery)){ foreach($field_array as $fa){$value = $row["$fa"]; //get the value of the field
$field =$fa; //get the actual field name

//store in an array to be easily accessable
if($queryCount > 1){$getInfo[$field][] =$value;
}else{
$getInfo[$field] = $value; } } } //need to check if any data was actually recieved if(empty($getInfo)){
$response = array("obj"=>"$table","instruct"=>"get","status"=>"pass","msg"=>"Graph ran successfully, but no data was retrieved. Check the graph coordinates and try again.","data"=>"empty");
}else{
$response = array("obj"=>"$table","instruct"=>"get","request"=>$request,"data"=>$getInfo,"status"=>"pass","msg"=>"Data was successfully retrieved.");
}

}else{
$response = array("obj"=>"$table","instruct"=>"get","status"=>"fail","msg"=>"Graph couldn't get any data because all required parameters weren't passed.");
}

//if successful, this will return a JSON array, with the final value of "data".

return json_encode($response); } public static function post($data){
//used to create new data

//required parameters
$db =$data['db'];
$table =$data['table'];
$fieldData =$data['field'];

if($db&&$table){
if($fieldData){ //need to create the new record here$field_id = $table."_id"; //apart of the naming standards setup that the system must use to work...$createNewRecord = database::conQuer("INSERT INTO $db.$table ($field_id) VALUES (NULL);",$db);

if($createNewRecord){$newFieldId = mysql_insert_id(); //gets the lasat id of the query

foreach($fieldData as$value){
$field = array_search("$value", $fieldData);$insert = database::conQuer("UPDATE $db.$table SET $field = '$value' WHERE $table.$field_id =$newFieldId LIMIT 1 ;",$db);
}

$response = array("obj"=>"$table","instruct"=>"post","data"=>$fieldData,"status"=>"pass","msg"=>"Graph posted the data successfully."); }else{$response = array("obj"=>"$table","instruct"=>"post","status"=>"fail","msg"=>"Graph didn't run because there was a query error when creating a new value."); } }else{$response = array("obj"=>"$table","instruct"=>"post","status"=>"fail","msg"=>"Graph didn't run because no data values were passed."); } }else{$response = array("obj"=>"$table","instruct"=>"post","status"=>"fail","msg"=>"Graph didn't run because the required parameters were not passed."); } return json_encode($response);
}

public static function put($data){ //used to update existing data$db = $data['db'];$table = $data['table'];$field = $data['field'];$value = $data['value'];$ref = $data['ref'];$refval = $data['refval']; if($db&&$table&&$field&&$value&&$ref&&$refval){$put = database::conQuer("UPDATE $db.$table SET $field = '$value' WHERE $table.$ref = '$refval';",$db);
if($put){$putData = array($field,$value,$ref,$refval);

$response = array("obj"=>"$table","instruct"=>"put","data"=>$putData,"status"=>"pass","msg"=>"Graph updated the data successfully."); }else{$response = array("obj"=>"$table","instruct"=>"put","status"=>"fail","msg"=>"Graph didn't put any data in because the query failed."); } }else{$response = array("obj"=>"$table","instruct"=>"put","status"=>"fail","msg"=>"Graph didn't put any data in because the required parameters were not passed."); } return json_encode($response);

}

public static function delete($data){ //used to remove existing data //used to update existing$db = $data['db'];$table = $data['table'];$ref = $data['ref'];$refval = $data['refval']; if($db&&$table&&$ref&&$refval){$delete = database::conQuer("DELETE FROM $table WHERE $table.$ref = '$refval'",$db); if($delete){
$deleteData = array($ref,$refval);$response = array("obj"=>"$table","instruct"=>"delete","data"=>$deleteData,"status"=>"pass","msg"=>"Graph deleted the data successfully.");
}else{
$response = array("obj"=>"$table","instruct"=>"delete","status"=>"fail","msg"=>"Graph didn't delete any data in because the query failed.");
}
}else{
$response = array("obj"=>"$table","instruct"=>"delete","status"=>"fail","msg"=>"Graph didn't delete any data in because the required parameters were not passed.");
}

return json_encode($response); } //end graph class } ?>  I would really appreciate the community's view on the class. What do you guys think about it? What can I do better? Critiques are welcome. I have used the class in several sites that I've built and it really seemed to help me. I hope that maybe it could help someone else. Please post any comments, suggestions, etc. • Looks clean but not sure about mysql_*. Users on newer versions of PHP are going to get deprecated warnings everywhere, I'd switch to PDO. – elclanrs Oct 9 '13 at 1:45 • Please include the code you'd like to have reviewed. The link could die at any time, rendering the question useless. – Jamal Oct 9 '13 at 1:46 • Yes. That is a good idea. I will update that... – Daniel Sikes Oct 9 '13 at 1:46 ## 1 Answer A few points: 1. You should avoid mysql_query. It's deprecated. Use PDO or mysqli. • Read up on SQL injection attacks and use parametrized queries. 2. Your indentation seems to be all over the place (not sure if that's a result of copy-n-paste, haven't checked the github sources) • Use more spaces - they are cheap (e.g. }else{ should be at least } else {); 3. Comments should mainly describe reasons for non-obvious decision. Most of your comments seem to just describe what the code below is doing which you can see by looking at the code. Example: //decode json data$data = json_decode($jsonData);  Those comments add no value. Furthermore reiterating what the code does in the comments will usually lead to comments conflicting with code once a few rounds of refactoring have taken place. 4. Why is this called Graph? Graph means usually one of two things to me: Either a data structure with nodes and edges connecting the nodes or the visualization of some data. I can see neither of those here really. 5. You repeat a lot of string literals in your response creation. This should be refactored into a method. Makes refactoring later easier and avoids typos. Maybe even create a class encapsulating the response so you can do r = response::build(...); r.Add(key, value); r.toJson(); or something along those lines. • Most prominently you use instrcut in your database class but instruct in your graph class. I'd actually use command (instructions in computer terms are usually associated with more low level operations but YMMV) 6. Consider dealing with error cases first to reduce nesting. This makes the code easier to read when there are multiple levels. E.g.: if (!$db || !$table) { return response::build(...).toJson(); } if (!$fieldData)
{
return response::build(...).toJson();
}

// handle the main case

return response::build().toJson();

In regards to the response refactoring:

class response
{
private data = array();

public static function build($object,$command, $status,$message)
{
$res = new response();$res->data['object'] = $object; ... return$res;
}

public function toJson()
{
return json_encode(\$this->data);
}
}

Disclaimer: Above code might not be 100% correct but should get the idea across. It's been a long time since I worked with PHP.

• Thanks for the help Chris Wue! To answer your question about why I named it graph. Originally I visualized the setup like this : The database represented a x value, the table represented a y value, and the field represented a z value. In my mind, I was in a way, building a graph of a database. It could be called anything really... I just left it what I had originally called it. You are the 2nd person to mention the switch to PDO. I will be adding that, and then uploading it again. Also, I copied and pasted the code so thats why the indenting is crazy. It looked fine in my editor. – Daniel Sikes Oct 9 '13 at 2:43
• On point # 5. Could you please be more specific? What should a response look like? What do you mean it should be a method? Show me what you are talking about. I want to make this class better, for sure. So that others can use it. Do you like the idea of being able to perform rest commands to a database? Thanks, Dan – Daniel Sikes Oct 9 '13 at 3:12
• @DanielSikes: Fair enough. I updated my answer – ChrisWue Oct 9 '13 at 3:37


text:
# Field norm

In mathematics, the (field) norm is a particular mapping defined in field theory, which maps elements of a larger field into a subfield.

## Formal definition

Let K be a field and L a finite extension (and hence an algebraic extension) of K.

The field L is then a finite dimensional vector space over K.

Multiplication by α, an element of L,

${\displaystyle m_{\alpha }\colon L\to L}$
${\displaystyle m_{\alpha }(x)=\alpha x}$,

is a K-linear transformation of this vector space into itself.

The norm, NL/K(α), is defined as the determinant of this linear transformation.[1]

If L/K is a Galois extension, one may compute the norm of α ∈ L as the product of all the Galois conjugates of α:

${\displaystyle \operatorname {N} _{L/K}(\alpha )=\prod _{\sigma \in \operatorname {Gal} (L/K)}\sigma (\alpha ),}$

where Gal(L/K) denotes the Galois group of L/K.[2] (Note that there may be a repetition in the terms of the product)

For a general field extension L/K, and nonzero α in L,

let σ1(α), ..., σn(α) be the roots of the minimal polynomial of α over K (roots listed with multiplicity and lying in some extension field of L); then

${\displaystyle \operatorname {N} _{L/K}(\alpha )=\left(\prod _{j=1}^{n}\sigma _{j}(\alpha )\right)^{[L:K(\alpha )]}}$.

If L/K is separable, then each root appears only once in the product (though the exponent, the degree [L:K(α)], may still be greater than 1).

## Examples

One of the basic examples of norms comes from quadratic field extensions ${\displaystyle \mathbb {Q} ({\sqrt {a}})/\mathbb {Q} }$ where ${\displaystyle a}$ is a square-free integer.

Then, the multiplication map by ${\displaystyle {\sqrt {a}}}$ on an element ${\displaystyle x+y\cdot {\sqrt {a}}}$ is

${\displaystyle {\sqrt {a}}\cdot (x+y\cdot {\sqrt {a}})=y\cdot a+x\cdot {\sqrt {a}}.}$

The element ${\displaystyle x+y\cdot {\sqrt {a}}}$ can be represented by the vector

${\displaystyle {\begin{bmatrix}x\\y\end{bmatrix}},}$

since there is a direct sum decomposition ${\displaystyle \mathbb {Q} ({\sqrt {a}})=\mathbb {Q} \oplus \mathbb {Q} \cdot {\sqrt {a}}}$ as a ${\displaystyle \mathbb {Q} }$-vector space.

The matrix of ${\displaystyle m_{\sqrt {a}}}$ is then

${\displaystyle m_{\sqrt {a}}={\begin{bmatrix}0&a\\1&0\end{bmatrix}}}$

and the norm is ${\displaystyle N_{\mathbb {Q} ({\sqrt {a}})/\mathbb {Q} }=-a}$, since it is the determinant of this matrix.

In this example the norm was the square of the usual Euclidean distance norm in ${\displaystyle \mathbb {C} }$.

In general, the field norm is very different from the usual distance norm.

We will illustrate that with an example where the field norm can be negative.

#### Norm of Q(√2)

Consider the number field ${\displaystyle K=\mathbb {Q} ({\sqrt {2}})}$.

The Galois group of ${\displaystyle K}$ over ${\displaystyle \mathbb {Q} }$ has order ${\displaystyle d=2}$ and is generated by the element which sends ${\displaystyle {\sqrt {2}}}$ to ${\displaystyle -{\sqrt {2}}}$.

So the norm of ${\displaystyle 1+{\sqrt {2}}}$ is:

${\displaystyle (1+{\sqrt {2}})(1-{\sqrt {2}})=-1.}$

The field norm can also be obtained without the Galois group.

Fix a ${\displaystyle \mathbb {Q} }$-basis of ${\displaystyle \mathbb {Q} ({\sqrt {2}})}$, say:

${\displaystyle \{1,{\sqrt {2}}\}}$.

Then multiplication by the number ${\displaystyle 1+{\sqrt {2}}}$ sends

1 to ${\displaystyle 1+{\sqrt {2}}}$ and
${\displaystyle {\sqrt {2}}}$ to ${\displaystyle 2+{\sqrt {2}}}$.

So the determinant of "multiplying by ${\displaystyle 1+{\sqrt {2}}}$" is the determinant of the matrix which sends the vector

${\displaystyle {\begin{bmatrix}1\\0\end{bmatrix}}}$ (corresponding to the first basis element, i.e., 1) to ${\displaystyle {\begin{bmatrix}1\\1\end{bmatrix}}}$,
${\displaystyle {\begin{bmatrix}0\\1\end{bmatrix}}}$ (corresponding to the second basis element, i.e., ${\displaystyle {\sqrt {2}}}$) to ${\displaystyle {\begin{bmatrix}2\\1\end{bmatrix}}}$,

viz.:

${\displaystyle {\begin{bmatrix}1&2\\1&1\end{bmatrix}}.}$

The determinant of this matrix is −1.

### K-th root field extensions

Another easy class of examples comes from field extensions of the form ${\displaystyle \mathbb {Q} ({\sqrt[{p}]{a}})/\mathbb {Q} }$ where the prime factorization of ${\displaystyle a\in \mathbb {Q} }$ contains no ${\displaystyle p}$-th powers.

The multiplication map by ${\displaystyle {\sqrt[{p}]{a}}}$ of an element is

{\displaystyle {\begin{aligned}m_{\sqrt[{p}]{a}}(x)&={\sqrt[{p}]{a}}\cdot (a_{1}+a_{2}{\sqrt[{p}]{a}}+a_{3}{\sqrt[{p}]{a^{2}}}+\cdots +a_{p-1}{\sqrt[{p}]{a^{p-1}}})\\&=a_{1}{\sqrt[{p}]{a}}+a_{2}{\sqrt[{p}]{a^{2}}}+a_{3}{\sqrt[{p}]{a^{3}}}+\cdots +a_{p-1}a\end{aligned}}}

giving the matrix

${\displaystyle {\begin{bmatrix}0&0&\cdots &0&a\\1&0&\cdots &0&0\\0&1&\cdots &0&0\\\vdots &\vdots &\ddots &\vdots &\vdots \\0&0&\cdots &1&0\end{bmatrix}}}$

The determinant gives the norm

${\displaystyle N_{\mathbb {Q} ({\sqrt[{p}]{a}})/\mathbb {Q} }({\sqrt[{p}]{a}})=(-1)^{p-1}a=a.}$

### Complex numbers over the reals

The field norm from the complex numbers to the real numbers sends

x + iy

to

x2 + y2,

because the Galois group of ${\displaystyle \mathbb {C} }$ over ${\displaystyle \mathbb {R} }$ has two elements,

• the identity element and
• complex conjugation,

and taking the product yields (x + iy)(xiy) = x2 + y2.

### Finite fields

Let L = GF(qn) be a finite extension of a finite field K = GF(q).

Since L/K is a Galois extension, if α is in L, then the norm of α is the product of all the Galois conjugates of α, i.e.[3]

${\displaystyle \operatorname {N} _{L/K}(\alpha )=\alpha \cdot \alpha ^{q}\cdot \alpha ^{q^{2}}\cdots \alpha ^{q^{n-1}}=\alpha ^{(q^{n}-1)/(q-1)}.}$

In this setting we have the additional properties,[4]

• ${\displaystyle \forall \alpha \in L,\quad \operatorname {N} _{L/K}(\alpha ^{q})=\operatorname {N} _{L/K}(\alpha )}$
• ${\displaystyle \forall a\in K,\quad \operatorname {N} _{L/K}(a)=a^{n}.}$

## Properties of the norm

Several properties of the norm function hold for any finite extension.[5][6]

### Group homomorphism

The norm NL/K : L* → K* is a group homomorphism from the multiplicative group of L to the multiplicative group of K, that is

${\displaystyle \operatorname {N} _{L/K}(\alpha \beta )=\operatorname {N} _{L/K}(\alpha )\operatorname {N} _{L/K}(\beta ){\text{ for all }}\alpha ,\beta \in L^{*}.}$

Furthermore, if a in K:

${\displaystyle \operatorname {N} _{L/K}(a\alpha )=a^{[L:K]}\operatorname {N} _{L/K}(\alpha ){\text{ for all }}\alpha \in L.}$

If aK then ${\displaystyle \operatorname {N} _{L/K}(a)=a^{[L:K]}.}$

### Composition with field extensions

Additionally, the norm behaves well in towers of fields:

if M is a finite extension of L, then the norm from M to K is just the composition of the norm from M to L with the norm from L to K, i.e.

${\displaystyle \operatorname {N} _{M/K}=\operatorname {N} _{L/K}\circ \operatorname {N} _{M/L}.}$

### Reduction of the norm

The norm of an element in an arbitrary field extension can be reduced to an easier computation if the degree of the field extension is already known. This is

${\displaystyle N_{L/K}(\alpha )=N_{K(\alpha )/K}(\alpha )^{[L:K(\alpha )]}}$[6]

For example, for ${\displaystyle \alpha ={\sqrt {2}}}$ in the field extension ${\displaystyle L=\mathbb {Q} ({\sqrt {2}},\zeta _{3}),K=\mathbb {Q} }$, the norm of ${\displaystyle \alpha }$ is

{\displaystyle {\begin{aligned}N_{\mathbb {Q} ({\sqrt {2}},\zeta _{3})/\mathbb {Q} }({\sqrt {2}})&=N_{\mathbb {Q} ({\sqrt {2}})/\mathbb {Q} }({\sqrt {2}})^{[\mathbb {Q} ({\sqrt {2}},\zeta _{3}):\mathbb {Q} ({\sqrt {2}})]}\\&=(-2)^{2}\\&=4\end{aligned}}}

since the degree of the field extension ${\displaystyle L/K(\alpha )}$ is ${\displaystyle 2}$.

### Detection of units

An element ${\displaystyle \alpha \in {\mathcal {O}}_{K}}$ is a unit if and only if ${\displaystyle N_{K/\mathbb {Q} }(\alpha )=\pm 1}$.

For instance

${\displaystyle N_{\mathbb {Q} (\zeta _{3})/\mathbb {Q} }(\zeta _{3})=1}$

where

${\displaystyle \zeta _{3}^{3}=1}$.

Then any number field ${\displaystyle {\mathcal {O}}_{K}}$ containing ${\displaystyle \zeta _{3}}$ has it as a unit.

## Further properties

The norm of an algebraic integer is again an integer, because it is equal (up to sign) to the constant term of the characteristic polynomial.

In algebraic number theory one defines also norms for ideals. This is done in such a way that if I is a nonzero ideal of OK, the ring of integers of the number field K, N(I) is the number of residue classes in ${\displaystyle O_{K}/I}$ – i.e. the cardinality of this finite ring. Hence this ideal norm is always a positive integer.

When I is a principal ideal αOK then N(I) is equal to the absolute value of the norm to Q of α, for α an algebraic integer.


text:
# A Brief Survey of Generative Models

A high-level summary of various generative models including Variational Autoencoders (VAE), Generative Adverserial Networks (GAN), and their notable extentions and generalizations, such as f-GAN, Adversarial Variational Bayes (AVB), Wasserstein GAN, Wasserstein Auto-Encoder (WAE), Cramer GAN and etc

## Overview

This page is a high-level summary of various generative models with little explanations. Models to cover are as follows:

Variational Autoencoders (VAE)

• Kingma, Diederik P., and Max Welling. "Auto-encoding variational bayes." arXiv preprint arXiv:1312.6114 (2013).

Extention to VAE to use non-Gaussian encoders

• Mescheder, Lars, Sebastian Nowozin, and Andreas Geiger. "Adversarial variational bayes: Unifying variational autoencoders and generative adversarial networks." arXiv preprint arXiv:1701.04722 (2017).

• Goodfellow, Ian, et al. "Generative adversarial nets." Advances in neural information processing systems. 2014.

Generalized divergence minimization GAN ($f$-GAN)

• Nowozin, Sebastian, Botond Cseke, and Ryota Tomioka. "f-gan: Training generative neural samplers using variational divergence minimization." Advances in Neural Information Processing Systems. 2016.

Wasserstein GAN (WGAN)

• Arjovsky, Martin, Soumith Chintala, and Léon Bottou. "Wasserstein gan." arXiv preprint arXiv:1701.07875 (2017).

• Makhzani, Alireza, et al. "Adversarial autoencoders." arXiv preprint arXiv:1511.05644 (2015).

Wasserstein Auto-Encoder (WAE)

• Tolstikhin, Ilya, et al. "Wasserstein Auto-Encoders." arXiv preprint arXiv:1711.01558 (2017).

Cramer GAN

• Bellemare, Marc G., et al. "The Cramer Distance as a Solution to Biased Wasserstein Gradients." arXiv preprint arXiv:1705.10743 (2017).

## VAE

### Model setup:

• Recognition model: $q_\phi(z \vert x) = \mathcal N(\mu=h_1(x), \sigma^2 \mathbf I=h_2(x)\mathbf I)$
• Assumed fixed prior: $p(z) = \mathcal N(0,\mathbf I)$
• Generation model: $p_\theta(x \vert z) = \mathcal N(\mu=g_1(z), \sigma^2 \mathbf I=g_2(z)\mathbf I)$
• Implied (but intractable) posterior: $p_\theta(z\vert x)$

### Key equations:

$$$$\begin{split} \log p_\theta (x^i) = D_{KL}(q_\phi(z\vert x^i)\| p_\theta(z\vert x^i) + \mathcal L(\theta, \phi, x^i) \end{split}$$$$ $$$$\begin{split} \mathcal L(\theta, \phi, x^i) &= \mathbb{E}_{z\sim q_\phi(z\vert x^i)}[\log p_\theta(x^i,z) - \log q_\phi(z\vert x^i)]\\\\ &= \mathbb{E}_{z\sim q_\phi(z\vert x^i)}[\log p_\theta(x^i,z)] + H[q_\phi(z\vert x^i)]\\\\ &= \mathbb{E}_{z\sim q_\phi(z\vert x^i)}[\log p_\theta(x^i \vert z)] - D_{KL}[q_\phi(z\vert x^i) \| p(z)]\\\\ \end{split} \label{vae_elbo}$$$$

### Optimization objective:

$$\hat\theta, \hat\phi = \argmax_{\theta, \phi} \sum_i \mathcal L(\theta, \phi, x^i)$$

Difficulties in calculating $\mathcal L(\theta, \phi, x^i)$:

• Due to the generality of $q$ and $p$ (typically a neural network), the expectation in $\ref{vae_elbo}$ does not have an analytical form. So we need to resort to Monte Carlo estimation.
• Furthermore, direct sampling $z$ according to $q$ poses difficulty in taking derivative against parameters $\phi$ that parameterizes the distribution $q$.

Solution: Reparameterization Trick

Find smooth and invertible transformation $z=g_\phi(\epsilon)$ such that with $\epsilon$ drawn from a fixed (non-parameterized) distribution $p(\epsilon)$ we have $z \sim q(z; \phi)$, so

$$\mathbb{E}_{z\sim q(z;\phi)}[f(z)] = \mathbb{E}_{\epsilon\sim p(\epsilon)}[f(g_\phi(\epsilon))]$$

For the Normal distribution used here ($q_\phi(z\vert x)$), it is convenient to use location-scale transformation, $z=\mu+\sigma * \epsilon$ with $\epsilon \sim \mathcal N(0,\mathbf I)$.

$$$$\widetilde{\mathcal{L}}(\theta, \phi, x^i) = \frac{1}{L} \sum_{l=1}^L \left( \log p_\theta(x^i \vert z^{i,l})] - D_{KL}[q_\phi(z^{i,l}\vert x^i) \| p(z^{i,l}) \right)$$$$ $$z^{i,l} = \mu_{x^i} + \sigma_{x^i} * \epsilon ^{i,l} ~~\text{and}~~ \epsilon^l \sim \mathcal N(0,\mathbf I)$$

For total $N$ data points with mini batch size $M$:

$$$$\begin{split} {\mathcal L}(\theta, \phi; X) = \sum_{i=1}^N \mathcal L(\theta, \phi, x^i) \approx \widetilde {\mathcal L^M}(\theta, \phi; X) = \frac{N}{M} \sum_{i=1}^M \widetilde {\mathcal L}(\theta, \phi, x^i) \end{split}$$$$

For sufficiently large batch size $M$, the inner loop sample size $L$ can be set to 1. Due to stochastic mini batch gradient descent and stochastic expectation estimation, this is also called doubly stochastic estimation.

### Using non-Gaussian encoders

Todo: discuss AVB paper

### Gumble trick for discrete latent variables

Ref for this section:

1. Gumble max trick https://hips.seas.harvard.edu/blog/2013/04/06/the-gumbel-max-trick-for-discrete-distributions/
2. Balog, Matej, et al. "Lost Relatives of the Gumbel Trick." arXiv preprint arXiv:1706.04161 (2017).
3. Jang, Eric, Shixiang Gu, and Ben Poole. "Categorical reparameterization with gumbel-softmax." arXiv preprint arXiv:1611.01144 (2016).

Gumble distribution:

## $f$-GAN and GAN

### Prelude on $f$-divergence and its variational lower bound

The f-divergence family

$$$$D_f = \int_{\mathcal X} q(x) ~ f\left( \frac{p(x)}{q(x)}\right) dx \label{fdiv}$$$$

where the generator function $f: \mathbb{R}_{+} \rightarrow \mathbb{R}$ is a convex, lower-semicontinuous function satisfying $f(1) = 0$.

Every convex, lower-semicontinuous function has a convex conjugate function $f^c$, also known as Fenchel conjugate. This function is defined as

$$$$f^c(t) = \underset {u \in \text{dom}_f}{\text{sup}} \{ut - f(u)\}$$$$

Function $f^c$ is again convex and lower-semicontinuous and the pair $(f,f^c)$ is dual to each other, i.e. $\left(f^{c}\right)^c=f$. So we can represent $f$ as

$$$$f(t) = \underset {t \in \text{dom}_{f^c}}{\text{sup}} \{tu - f^c(t)\}$$$$

With this we can establish a lower bound for estimating the f-divergence in general

$$$$\begin{split} D_f(P \| Q) &= \int_{\mathcal X} q(x) \underset {t \in \text{dom}_{f^c}}{\text{sup}} \{t \frac{p(x)}{q(x)} - f^c(t)\} dx \\\\ & \ge \underset {T \in {\mathcal T}} {\text{sup}} \int_{\mathcal X} \left( p(x)T(x) - q(x)f^c(T(x)) \right) dx \\\\ & = \underset {T \in {\mathcal T}} {\text{sup}} \left( \mathbb{E}_{x\sim P}[T(x)] - \mathbb{E}_{x\sim Q}[f^c(T(x))] \right) \end{split} \label{f_lowerbound}$$$$

where $\mathcal T$ is an arbitrary class of functions $T: \mathcal X \rightarrow \mathbb R$. The inequality is due to Jessen's inequality and constraints imposed by $\mathcal T$.

The bond is tight for

$$$$T^*(x) = f' \left(\frac{p(x)}{q(x)} \right)$$$$

Suppose our goal is to come up with a distribution $Q$ (model) that is close to $P$ (the data distribution) and the similarity score (loss) is measured by $D_f(P | Q)$. However the direct calculation of $\ref{fdiv}$ is intractable, such as the case where the functional form of $P$ is unknown and $Q$ is a complex model parameterized by a neural network.

To be specific:

• Evaluating $q(x)$ at any $x$ is easy, but integrating it is hard due to lack of easy functional form.
• For $p(x)$, we do not know how to evaluate it at any $x$
• Sampling from both $P$ and $Q$ are easy. Because drawing from data set approximates $x \sim P$ and we can make the model $Q$ take random vectors as input which are easy to produce.

Fortunately, we can sample from both of them easily. In this case, $\ref{f_lowerbound}$ offers a way to estimate the lower bound of the divergence. We would need to maximize this lower bound by changing $T$ so that it is close to the true divergence, then minimize it over $Q$. This is formally stated as follows.

$$$$F(\theta, \omega) = \mathbb{E}_{x\sim P}[T_\omega(x)] + \mathbb{E}_{x\sim Q_\theta}[-f^c(T_\omega(x))]$$$$ $$$$\hat \theta = \argmin_\theta \max_\omega F(\theta, \omega)$$$$

To ensure that the output of $T_\omega$ respects the domain of ${f^c}$, we define $T_\omega(x) = g_f(V_\omega(x))$, where $V_\omega: \mathcal X \rightarrow \mathbb R$ without any range constraints on the output and $g_f: \mathbb R \rightarrow \text{dom}_{f^c}$ is an output activation function specific to the $f$-divergence used with suitable output ranges.

### GAN

For the original GAN, with a divergence target similar to Jensen-Shannon

$$$$F(\theta, \omega) = \mathbb{E}_{x\sim P}[\log D_\omega(x)] + \mathbb{E}_{x\sim Q_\theta}[\log(1-D_\omega(x))]$$$$

with $D_\omega = 1/(1+e^{-V_\omega(x)})$
which corresponds to the following

$$g_f(\nu)= \log(1/(1+e^{-\nu}))$$ $$T_\omega(x) = \log (D_\omega(x)) = g_f(V_\omega(x))$$ $$f^c(t) = -\log (1-\exp(t))$$ $$\log (1-D_\omega(x)) = -f^c(T_\omega(x))$$

### Practical considerations in adversarial training

Todo: log trick, DCGAN heuristics
Name $D_f(P\vert Q)$ Generator $f(u)$ $T^*(x)$
Forward KL $\int p(x) \log \frac{p(x)}{q(x)} dx$ $u\log u$ $1 +\log \frac{q(x)}{p(x)}$
Reverse KL $\int q(x) \log \frac{p(x)}{q(x)} dx$ $-\log u$ $- \frac{q(x)}{p(x)}$
Jensen-Shannon $\frac{1}{2} \int p(x) \log \frac{2p(x)}{p(x)+q(x)} + q(x) \log \frac{2q(x)}{p(x)+q(x)} dx$ $u\log u - (u+1) \log \frac{u+1}{2}$ $\log \frac{2p(x)}{p(x)+q(x)}$
GAN $\int p(x) \log \frac{2p(x)}{p(x)+q(x)} + q(x) \log \frac{2q(x)}{p(x)+q(x)} dx -\log(4)$ $u\log u - (u+1) \log (u+1)$ $\log \frac{p(x)}{p(x)+q(x)}$
Name Conjugate $f^c(t)$ $\text{dom}_{f^c}$ Output activation $g_f$ $f'(1)$
Forward KL $\exp(t-1)$ $\mathbb R$ $\nu$ $1$
Reverse KL $-1-\log(-t)$ $\mathbb R_{-}$ $-\exp(\nu)$ $-1$
Jensen-Shannon $-\log(2-\exp(t))$ $t < \log(2)$ $\log(2) - \log(1+\exp(-\nu))$ $0$
GAN $-\log(1-\exp(t))$ $\mathbb R_{-}$ $- \log(1+\exp(-\nu))$ $-\log(2)$

## WGAN and WAE

### Optimal transport (OT)

Kantorovich formulated the optimization target in optimal transport problems as follows

$$$$W_c(P_X, P_G) = \underset{\Gamma \in \mathcal P(x \sim P_X, y \sim P_Y)}{\text{inf}} \mathbb{E}_{x,y \sim \Gamma}[c(x,y)]$$$$

where $\mathcal P(X\sim P_X, Y\sim P_Y)$ is a set of all join distributions of $(X,Y)$ with marginals $P_X$ and $P_Y$.

### Wasserstein distance

When $c(x,y) = | x-y | ^p$ for $p \ge 1$, $W_c^{1/p}$ is called p-Wasserstein distance.

$$$$W_p(P_X, P_G) = \underset{\Gamma \in \mathcal P(x \sim P_X, y \sim P_Y)}{\text{inf}} \mathbb{E}_{x,y \sim \Gamma}[\|x - y\|^p]$$$$

The optimization problem is highly intractable in general, due to the constraint. However when $p=1$, Kantorovich-Rubinstein duality holds:

$$$$W_1(P_X, P_G) = \underset{f \in \text{\{1-Lipschitz\}}}{\text{sup}} \mathbb{E}_{x\sim P_X}[f(x)] - \mathbb{E}_{y\sim P_Y}[f(y)]$$$$

The family of divergences from $f$-divergence only consider the relative probability (the ratio between two probability density functions) and do not measure the closeness of the underlying outcomes. With disjoint support or overlapping support but intersections that yield zero measure, the divergence between a target distribution and a $\theta$-parameterized distribution might not be continuous with respect to $\theta$. Wasserstein distance on the other hand does take into account the underlying topology of the outcomes and is continuous and differentiable almost everywhere with respect to $\theta$ and thus almost always provide useful gradient for optimization.

### Wasserstein GAN (WGAN)

Following the dual form of $W_1$, we can form a generative-adversarial model for a data distribution $P_D$ and model $Q_\theta$ with auxiliary function $f$ that is 1-Lipschitz continuous.

$$$$\hat \theta = \argmin_\theta \underset{f \in \text{\{1-Lipschitz\}}}{\text{sup}} \mathbb{E}_{x\sim P_D}[f(x)] - \mathbb{E}_{x\sim Q_\theta}[f(x)]$$$$

Practical considerations for WGAN

Todo: Gradient clipping with K-Lipschitz constraint on $f$; Soft gradient penalty (WGAN-GP)

### Wasserstein Auto-encoder (WAE)

Rather than working with the dual or Wasserstein distance, which only holds for $W_1$, we can also work with the primal form directly. As shown in Tolstikhin, Ilya, et al. “Wasserstein Auto-Encoders.” the following holds when we have a deterministic decoder mapping latent variable $Z$ to $Y$ through $y=G(z)$:

$$$$\begin{split} W_c(P_X, P_G) = W_c^\dagger(P_X, P_G) &= \underset{P \in \mathcal P(x \sim P_X, z \sim P_Z)}{\text{inf}} \mathbb{E}_{x,y \sim P}[c(x, G(z)]\\\\ &= \underset{Q: Q_Z = P_Z}{\text{inf}} \mathbb{E}_{x\sim P_X} \mathbb{E}_{z \sim Q(Z\vert X)}[c(x, G(z)] \end{split}$$$$

The constraint put on $Q(Z\vert X)$ is that its marginal needs to equal to $P(Z)$. To have a feasible optimization problem we relax this constraint with the following constraint-free optimization target with a penalty that assess the closeness between $Q(Z)$ and $P(Z)$ via any reasonable divergence. This new objective is named penalized optimal transport (POT).

$$$$D_{POT/WAE}(P_X, P_G) := \underset{Q \in \mathcal Q}{\text{inf}} \mathbb{E}_{x\sim P_X} \mathbb{E}_{z \sim Q(Z\vert X)}[c(x, G(z)] + \lambda \cdot D_{Z} (Q_Z, P_Z)$$$$

If the divergence between $P_Z$ and $Q_Z$ is intractable to directly calculate, we could use generative-adversarial training to approximate it (see $f$-GAN).

Note: if decoder is probabilistic instead of deterministic, we would only have $W_c(P_X, P_G) \le W_c^\dagger(P_X, P_G)$, so we are minimizing an upper bound of the true OT cost.
Thought: the original paper used JS divergence for $D_Z$, how about we use Wasserstein distance for $D_Z$.
Todo: discuss connections to AAE paper


text:
# $\mathbb{Z}$ under finite-complement topology [CSIR-2015]

Consider the set $$\mathbb{Z}$$ of integers with the topology $$\tau$$ in which a subset is closed iff it is empty, or $$\mathbb{Z}$$, or finite. Which of the following is true?

1. $$\tau$$ is the subspace topology induced from the usual topology on $$\mathbb{R}$$.

2. $$\mathbb{Z}$$ is compact in the topology $$\tau$$.

3. $$\mathbb{Z}$$ is hausdorff in the topology $$\tau$$.

4. Every infinite subset of $$\mathbb{Z}$$ is dense in the topology $$\tau$$.

I think 1 is true since any closed set in $$\tau$$ is $$\mathbb{Z}\cap [a,b]$$(or a finite union of these sets of this type) a,b integers(for finite sets) or the entire real line when the closed set is entire $$\mathbb{Z}$$ .

Since any non-trivial open set is basically $$\mathbb{Z}$$-{a finite set of integers}.

3 is false because given two integers $$a,b$$ if i try to look for two disjoint open sets $$U_a$$,$$U_b$$ they will always have a infinite number of points in their intersection.

4 is true since if i take any open set in $$\tau$$(which is of the above type) it must intersect a infinite subset of $$\mathbb{Z}$$.

I am not able to say anything about compactness. Further are all of the above reasoning correct?

• Closed set in $\mathbb{R}$ is not necessarily a finite union of $[a,b]$. Jun 23 '15 at 15:52
• How was the exam overall ?How much did you score?@Prayagdeep Jun 23 '15 at 16:16
• @learnmore i guess i will end up with a tip over 100.Is this the right place to continue this discussion...i have some more doubts on the questions, is there another place where we can discuss..? Jun 23 '15 at 17:28
• A few years late, but in case anyone else is wondering: MSE has a chat room, which would be the best place for discussing contests after they've completed. Jul 5 '17 at 5:49

Since when $\mathbb Z$ is given subspace topology it becomes discrete i.e every one point sets are both closed and open which is not true with $\tau$.
$\tau$ is compact since if you consider the open cover $\{U_{\alpha}:\alpha \in I\}$ then $\mathbb Z\setminus U_{\alpha}$ is finite and hence the remaining points can be covered with finite number of elements from $I$
(1) is false. The subspace topology on $\Bbb Z$ induced by the standard topology on $\Bbb R$ is the discrete topology. So, $\{1\}$ is open in the subspace topolobgy, but not in the finite complement topology.


text:
# Thread: HSC 2017 MX2 Integration Marathon (archive)

1. ## Re: HSC 2017 MX2 Integration Marathon

$\noindent a) \int \frac{\text{d}x}{x \sqrt{x^{10} \pm x^5 - 1}} \\\\ b) \int \frac{\text{d}x}{x \sqrt{x^{10} \pm x^5 + 1}}$
$a) \quad Let\quad I=\int { \frac { dx }{ x\sqrt { { x }^{ 10 }\pm { x }^{ 5 }+1 } } } \\ Let\quad u={ x }^{ 5 }\\ du=5{ x }^{ 4 }dx\\ I=\frac { 1 }{ 5 } \int { \frac { du }{ u\sqrt { u^{ 2 }\pm u-1 } } } \\ \quad =\frac { 1 }{ 5 } \int { \frac { du }{ u\sqrt { { (u\pm \frac { 1 }{ 2 } ) }^{ 2 }-\frac { 5 }{ 4 } } } }\quad$
$Let\quad I=\quad =\frac { 1 }{ 5 } \int { \frac { du }{ u\sqrt { { (u\pm \frac { 1 }{ 2 } ) }^{ 2 }+\frac { 3 }{ 4 } } } } \\ Let\quad w=u\pm \frac { 1 }{ 2 } \\ dw=du\\ I=\frac { 1 }{ 5 } \int { \frac { dw }{ (w\mp \frac { 1 }{ 2 } )\sqrt { { w }^{ 2 }+\frac { 3 }{ 4 } } } } \\ =\frac { 1 }{ 5 } (ln\left| w\mp \frac { 1 }{ 2 } \right| -ln\left| 2\sqrt { { w }^{ 2 }+\frac { 3 }{ 4 } } \pm w+\frac { 3 }{ 2 } \right| )+C\quad \\ =\frac { 1 }{ 5 } (ln\left| u \right| -ln\left| 2\sqrt { { (u\pm \frac { 1 }{ 2 } ) }^{ 2 }+\frac { 3 }{ 4 } } \pm u\mp \frac { 1 }{ 2 } +\frac { 3 }{ 2 } \right| )+C\\ =\frac { 1 }{ 5 } (ln\left| { x }^{ 5 } \right| -ln\left| 2\sqrt { { ({ x }^{ 5 }\pm \frac { 1 }{ 2 } ) }^{ 2 }+\frac { 3 }{ 4 } } \pm { x }^{ 5 }\mp \frac { 1 }{ 2 } +\frac { 3 }{ 2 } \right| )+C\\ =\frac { 1 }{ 5 } (ln\left| { x }^{ 5 } \right| -ln\left| 2\sqrt { { x }^{ 10 }\pm { x }^{ 5 }+1 } \pm { x }^{ 5 }\mp \frac { 1 }{ 2 } +\frac { 3 }{ 2 } \right| )+C\\ \\$
$b)\quad Let\quad I=\int { \frac { dx }{ x\sqrt { { x }^{ 10 }\pm { x }^{ 5 }-1 } } } \\ Let\quad u={ x }^{ 5 }\\ du=5{ x }^{ 4 }dx\\ I=\frac { 1 }{ 5 } \int { \frac { du }{ u\sqrt { u^{ 2 }\pm u-1 } } } \\ \quad =\frac { 1 }{ 5 } \int { \frac { du }{ u\sqrt { { (u\pm \frac { 1 }{ 2 } ) }^{ 2 }-\frac { 5 }{ 4 } } } } \\ Let\quad w=u\pm \frac { 1 }{ 2 } \\ dw=du\\ I=\frac { 1 }{ 5 } \int { \frac { dw }{ (w\mp \frac { 1 }{ 2 } )\sqrt { { w }^{ 2 }-\frac { 5 }{ 4 } } } } \\ I=\pm \frac { 1 }{ 5 } \tan ^{ -1 }{ (\frac { \frac { w }{ 2 } \mp \frac { 5 }{ 4 } }{ \sqrt { { w }^{ 2 }-\frac { 5 }{ 4 } } } )+C } \\ I=\pm \frac { 1 }{ 5 } \tan ^{ -1 }{ (\frac { \frac { u }{ 2 } \pm \frac { 1 }{ 4 } \mp \frac { 5 }{ 4 } }{ \sqrt { { u }^{ 2 }\pm u-1 } } )+C } \\ I=\pm \frac { 1 }{ 5 } \tan ^{ -1 }{ (\frac { u\pm 2 }{ 2\sqrt { { u }^{ 2 }\pm u-1 } } )+C } \\ I=\pm \frac { 1 }{ 5 } \tan ^{ -1 }{ (\frac { { x }^{ 5 }\mp 2 }{ 2\sqrt { { x }^{ 10 }\pm { x }^{ 5 }-1 } } )+C }$
though in reality these questions are trivial and can easily be done by inspection

2. ## Re: HSC 2017 MX2 Integration Marathon

If b>a>0 evaluate:

$\int_{\tfrac{1}{b}}^{\tfrac{1}{a}} \dfrac{\tan^{-1}{ax} - \cot^{-1}{bx}}{x} \, \text{d}x$

3. ## Re: HSC 2017 MX2 Integration Marathon

If b>a>0 evaluate:

$\int_{\tfrac{1}{b}}^{\tfrac{1}{a}} \dfrac{\tan^{-1}{ax} - \cot^{-1}{bx}}{x} \, \text{d}x$
Trivial by inspection

$\noindent I=\int_{\tfrac{1}{b}}^{\tfrac{1}{a}} \dfrac{\tan^{-1}{ax} - \cot^{-1}{bx}}{x} \, \text{d}x \\\\ Let x=\dfrac{1}{y} \\\\ \text{d}x = \dfrac{-\text{d}y}{y^2} \\\\ I=\int _b^a \dfrac{\tan^{-1}{\tfrac{a}{y}}-\cot^{-1}{\tfrac{b}{y}}}{\tfrac{1}{y}} \dfrac{-\text{d}y}{y^2} \\\\ I =\int_a^b \dfrac {\tan^{-1}{\tfrac{a}{y}}- \cot^{-1}{\tfrac{b}{y}}}{y} \text{d}y \\\\ I =\int_a^b \dfrac{\cot^{-1}{\tfrac{y}{a}} - \tan^{-1}{\tfrac{y}{b}}}{y} \text{d}y \\\\ Let y=abx \\\\I =\int_{\tfrac{1}{b}}^{\tfrac{1}{a}} \dfrac{\cot^{-1}{bx} - \tan^{-1}{ax}}{abx} \text{d}x \\\\ \quad =-ab\int _{ \tfrac{1}{b}}^{ \tfrac{1}{a}} \frac{\tan^{-1}{ax}-\cot^{-1}{bx}}{x} \text{d}x \\\\ \quad I=-abI\\\\ \therefore \quad I=0 since ab > 0 > -1$
Note the potential tautology if ab=-1 cannot occur as b>a>0.

4. ## Re: HSC 4U Integration Marathon 2017

$\noindent \int \frac{\sqrt{\cos{2x}}}{\sin{x}} \text{d}x = \int \frac{\sqrt{2\cos^2{x}-1}}{1-\cos^2{x}} \sin{x} \text{d}x \\\\ \sqrt{2} \cos{x} = \sec{\theta} \Rightarrow -\sqrt{2} \sin{x} \text{d}x = \sec{\theta} \tan{\theta} \text{d}\theta \\\\ \sqrt{2}\int \frac{\sqrt{(\sqrt{2} \cos{x})^2 -1}}{(2\cos{x})^2 -2} (-\sqrt{2} \sin{x} \text{d}x) \\\\ = \sqrt{2} \int \frac{\sec{\theta} \tan^2{\theta}}{\sec^2{\theta} -2} \text{d}\theta = \sqrt{2} \int \frac{\sin^2{\theta} \cos{\theta} \, \, \text{d}(\theta)}{(1-\sin^2{\theta}) (2\sin^2{\theta}-1)} \\\\ \sin^2{\theta} \equiv (2\sin^2{\theta}-1) + (1-\sin^2{\theta}) \\\\ \Rightarrow \frac{\sin^2{\theta}}{(1-\sin^2{\theta})(2\sin^2{\theta}-1)} \equiv \sec^2{\theta} + \frac{1}{2\sin^2{\theta}-1} \\\\ \sqrt{2}\int \sec{\theta} \, \text{d}\theta - \sqrt{2} \int \frac{\text{d}(\sin{\theta})}{1-2\sin^2{\theta}}$

$\noindent = \sqrt{2} \log{(\sec{\theta} + \tan{\theta})} + \frac{1}{2} \log{\left(\frac{1-\sqrt{2}\sin{\theta}}{1+\sqrt{2}\sin{\theta}} \right)} + \mathcal{C} \\\\ \sec{\theta} = \sqrt{2}\cos{x},\, \tan{\theta} = \sqrt{\cos{2x}}, \sin{\theta} = \frac{\sqrt{\cos{2x}}}{\sqrt{2} \cos{x}} \\\\ \sqrt{2} \log(\sqrt{2}\cos{x} + \sqrt{\cos{2x}}) + \frac{1}{2} \log{\left( \frac{\cos{x} - \sqrt{\cos{2x}}}{\cos{x} + \sqrt{\cos{2x}}} \right)} + \mathcal{C}$

Simplifying the derivative in Mathematica confirms the result.
Hey, sorry for being so bad at this subject, but how do you know to use that substitution? Thinking process? What are some key things to look out for when attempting a question?

5. ## Re: HSC 4U Integration Marathon 2017

Originally Posted by EngineeringHelp
Hey, sorry for being so bad at this subject, but how do you know to use that substitution? Thinking process? What are some key things to look out for when attempting a question?
Eliminate irrational functions using whatever you have at your disposal.

The only three ways of doing so are:

Classical Trigonometric Substitution (High School level)

Hyperbolic Trigonometric Substitution (University level)

Euler Substitution (Extracurricular)

6. ## Re: HSC 4U Integration Marathon 2017

Eliminate irrational functions using whatever you have at your disposal.

The only three ways of doing so are:

Classical Trigonometric Substitution (High School level)

Hyperbolic Trigonometric Substitution (University level)

Euler Substitution (Extracurricular)
Ooh, so for the HSC all i should worry about is a Trig Sub? As it wouldn't be right to be asked questions requiring the other two, right?

7. ## Re: HSC 4U Integration Marathon 2017

Also, does anyone have any 'not so complicated' reduction-formula questions? If so, please post. Im assuming Paradoxica will be supplying solutions xD <3

8. ## Re: HSC 4U Integration Marathon 2017

Originally Posted by EngineeringHelp
Also, does anyone have any 'not so complicated' reduction-formula questions? If so, please post. Im assuming Paradoxica will be supplying solutions xD <3
http://prntscr.com/fnxmsv

9. ## Re: HSC 4U Integration Marathon 2017

Originally Posted by pikachu975
$\int_0^1 \frac{x^n}{\sqrt{x+1}} \, \text{d}x$
...ew

10. ## Re: HSC 4U Integration Marathon 2017

Originally Posted by pikachu975
Lmao use Mathematical Induction to prove? Stuff In

11. ## Re: HSC 2017 MX2 Integration Marathon

\begin{align*}I_n&=\int_0^1\left(\frac{x^n+x^{n-1}-x^{n-1}}{\sqrt{x+1}}\right)dx\\ &= \int_0^1 \frac{x^{n-1}(x+1)}{\sqrt{x+1}}-I_{n-1}\\ &\stackrel{IBP}{=}\frac{x^n}n \sqrt{x+1} \bigg|_0^1 - \frac{1}{2n}\int_0^1\frac{x^n}{\sqrt{x+1}}-I_{n-1}\\ \implies \left(1+\frac{1}{2n}\right)I_n&=\frac{\sqrt2}{n} - I_{n-1}\\ (2n+1)I_n&=2\sqrt2 - 2nI_{n-1} \end{align*}

12. ## Re: HSC 2017 MX2 Integration Marathon

Try this one.
$\cot 3x\sec^2(ln \sqrt{sin 3x})$

13. ## Re: HSC 2017 MX2 Integration Marathon

Originally Posted by stupid_girl
Try this one.
$\cot 3x\sec^2(ln \sqrt{sin 3x})$
d(lnsqrt(sin3x)) = 3/2 * cot3x * dx
2/3 * d(lnsqrt(sin3x)) = cot3x dx

Therefore integrand becomes 2/3 * sec^2 (lnsqrt(sin3x)) * d(lnsqrt(sin3x))
Integrate:
= 2/3 * tan(lnsqrt(sin3x)) + C

14. ## Re: HSC 2017 MX2 Integration Marathon

Originally Posted by pikachu975
d(lnsqrt(sin3x)) = 3/2 * cot3x * dx
2/3 * d(lnsqrt(sin3x)) = cot3x dx

Therefore integrand becomes 2/3 * sec^2 (lnsqrt(sin3x)) * d(lnsqrt(sin3x))
Integrate:
= 2/3 * tan(lnsqrt(sin3x)) + C

15. ## Re: HSC 2017 MX2 Integration Marathon

long time no integrals
here's a couple of good ones to try out!

$\noindent \int \frac{x^2-2}{\left ( x^4+5x^2+4 \right )\tan^{-1}\left ( x+\frac{2}{x} \right )} dx \\ \\ \\ \int e^x\left ( \frac{1-x}{1+x^2} \right ) ^2 dx \\ \\ \\ \int \frac{\sin x - \cos x}{\sqrt{\sin(2x)}} dx \\ \\ \\\int \left ( \frac{1}{\ln x}+\ln \left ( \ln x \right ) \right ) dx$

16. ## Re: HSC 2017 MX2 Integration Marathon

Originally Posted by jathu123
long time no integrals
here's a couple of good ones to try out!

$\noindent \int \frac{x^2-2}{\left ( x^4+5x^2+4 \right )\tan^{-1}\left ( x+\frac{2}{x} \right )} dx \\ \\ \\ \int e^x\left ( \frac{1-x}{1+x^2} \right ) ^2 dx \\ \\ \\ \int \frac{\sin x - \cos x}{\sqrt{\sin(2x)}} dx \\ \\ \\\int \left ( \frac{1}{\ln x}+\ln \left ( \ln x \right ) \right ) dx$
$\noindent 1) By inspection \frac{\left(x^2-2\right)}{\left(x^4+5x^2+4\right)} is the derivative of tan^{-1}\left(x+\frac{2}{x}\right) hence the integral becomes: \\ \\ \quad \frac{d\left(tan^{-1}\left(x+\frac{2}{x}\right)\right)}{tan^{-1}\left(x+\frac{2}{x}\right)} which equals ln\left|tan^{-1}\left(x+\frac{2}{x}\right)\right|+C$

$\noindent 4) let \quad u = lnx , \quad e^u du = dx \\ Integral becomes \int \:\left(\frac{1}{u}+lnu\right)e^udu and applying reverse product rule you get \quad xln(lnx) + C.$

17. ## Re: HSC 2017 MX2 Integration Marathon

Originally Posted by jathu123
long time no integrals
here's a couple of good ones to try out!

$\noindent \int \frac{x^2-2}{\left ( x^4+5x^2+4 \right )\tan^{-1}\left ( x+\frac{2}{x} \right )} dx \\ \\ \\ \int e^x\left ( \frac{1-x}{1+x^2} \right ) ^2 dx \\ \\ \\ \int \frac{\sin x - \cos x}{\sqrt{\sin(2x)}} dx \\ \\ \\\int \left ( \frac{1}{\ln x}+\ln \left ( \ln x \right ) \right ) dx$
$\noindent The second one can also be found using the reverse product rule. Begin by noting that$

\begin{align*}\left (\frac{1 - x}{1 + x^2} \right )^2 &= \frac{(1 - x)^2}{(1 + x^2)^2} = \frac{1 + x^2}{(1 + x^2)^2} - \frac{2x}{(1 + x^2)^2} = \frac{1}{1 + x^2} - \frac{2x}{(1 + x^2)^2} = f(x) + f'(x),\end{align*}

$\noindent where f(x) = \frac{1}{1 + x^2}.$

$\noindent So the integral becomes$

\begin{align*}\int e^x \left (\frac{1 - x}{1 + x^2} \right )^2 \, dx &= \int e^x (f(x) + f'(x)) \, dx = e^x f(x) + C = \frac{e^x}{1 + x^2} + C.\end{align*}

18. ## Re: HSC 2017 MX2 Integration Marathon

Originally Posted by jathu123
long time no integrals
here's a couple of good ones to try out!

$\noindent \int \frac{x^2-2}{\left ( x^4+5x^2+4 \right )\tan^{-1}\left ( x+\frac{2}{x} \right )} dx \\ \\ \\ \int e^x\left ( \frac{1-x}{1+x^2} \right ) ^2 dx \\ \\ \\ \int \frac{\sin x - \cos x}{\sqrt{\sin(2x)}} dx \\ \\ \\\int \left ( \frac{1}{\ln x}+\ln \left ( \ln x \right ) \right ) dx$
$\noindent There is quite a nice way to do the third one. We begin by noting that$

\begin{align*}\sin 2x &= 2 \sin x \cos x = 1 + 2 \sin x \cos x - 1 = (\cos^2 x + \sin^2 x) + 2\sin x \cos x - 1 = (\cos x + \sin x)^2 - 1.\end{align*}

$\noindent So, for the integral we have$

\begin{align*}\int \frac{\sin x - \cos x}{\sqrt{\sin 2x}} \, dx &= \int \frac{\sin x - \cos x}{\sqrt{(\cos x + \sin x)^2 - 1}} \, dx.\end{align*}

$\noindent Letting u = \cos x + \sin x, du = -(\sin x - \cos x) \, dx so$

\begin{align*}\int \frac{\sin x - \cos x}{\sqrt{\sin 2x}} \, dx &= -\int \frac{du}{\sqrt{u^2 - 1}}\\ &=-\ln |u + \sqrt{u^2 - 1}| + C\\ &= -\ln |\cos x + \sin x + \sqrt{\sin 2x}| + C.\end{align*}

19. ## Re: HSC 2017 MX2 Integration Marathon

Integrate 1/[(sinx)sqrt(sin2x)]

20. ## Re: HSC 2017 MX2 Integration Marathon

The new 2018 marathon has been created.

dan

Page 11 of 11 First ... 91011

There are currently 1 users browsing this thread. (0 members and 1 guests)

#### Posting Permissions

• You may not post new threads
• You may not post replies
• You may not post attachments
• You may not edit your posts
•


text:
# Events - April 2019

Thursday, April 25, 2019 — 3:00 PM EDT

## Geometry & Topology Seminar

**Rescheduled from April 19. Note new day/time/room**

Shubham Dwivedi, Department of Pure Mathematics, University of Waterloo

"A gradient flow of isometric $G_2$ structures"

Wednesday, April 17, 2019 — 3:30 PM EDT

## Computability Learning Seminar

Michael Deveau, Department of Pure Mathematics, University of Waterloo

"Non-triviality of the Bounded Jump Hierarchy - Part 3"

We finish the proof that for all $B$, there is some set $A$ with $B <_{bT} A <_{bT} B^b$. We determine $\hat{h}_A$ and $\hat{h}_C$ as promised, taking into account the restrictions imposed by the construction, and prove that our choices are correct.

MC 5413

Monday, April 15, 2019 — 11:30 AM EDT

## Model Theory Seminar

Nigel Pynn-Coates, University of Illinois at Urbana-Champaign

"Asymptotic valued differential fields and differential-henselianity"

Wednesday, April 10, 2019 — 3:30 PM EDT

## Analysis Seminar

Clifford Bearden, University of Texas at Tyler

"A module version of the weak expectation property"

Wednesday, April 10, 2019 — 3:30 PM EDT

## Computability Learning Seminar

Michael Deveau, Department of Pure Mathematics, University of Waterloo

"Non-triviality of the Bounded Jump Hierarchy - Part 2"

We continue with the proof that for all $B$, there is some set $A$ with $B <_{bT} A <_{bT} B^b$. We start the main construction itself, check that the requirements eventually settle and are met, and we investigate the correct $\hat{h}_A$ and $\hat{h}_C$ as promised.

MC 5413

Friday, April 5, 2019 — 2:30 PM EDT

## Geometry & Topology Seminar

Spiro Karigiannis, Department of Pure Mathematics, University of Waterloo

"A curious system of second order nonlinear PDEs for $\mathrm{U}(m)$-structures on manifolds"

Thursday, April 4, 2019 — 12:00 PM EDT

Nick Manor, Unviersity of Waterloo

"Exact Groups"

There are two notions of exactness for groups, but they are not known to be equivalent outside the discrete case. I will show how we can extend this equivalence to include two larger classes of groups.

MC 5403

Wednesday, April 3, 2019 — 4:30 PM EDT

## Pure Math and C&O Joint Colloquium

Ben Moore, University of Waterloo

"A proof of the Hell-Neˇsetˇril Dichotomy via Siggers Polymorphism"

Wednesday, April 3, 2019 — 3:30 PM EDT

## Analysis Seminar

Gholamreza Aghamollaei, Shahid Bahonar University of Kerman, Iran

"The numerical range of quaternion matrices with respect to involutions"

Wednesday, April 3, 2019 — 3:30 PM EDT

## Computability Learning Seminar

Michael Deveau, Department of Pure Mathematics, University of Waterloo

“Non-triviality of the Bounded Jump Hierarchy - Part 1”

Monday, April 1, 2019 — 4:00 PM EDT

## Pure Math Colloquium

Jan Minac, University of Western Ontario

"The 13th mysterious room of a palace of absolute Galois groups"

### May 2019

S M T W T F S
28
29
30
1
3
4
5
6
7
8
11
12
18
19
20
22
25
26
27
29
1
1. 2019 (176)
1. November (10)
2. October (19)
3. September (13)
4. August (7)
5. July (12)
6. June (18)
7. May (22)
8. April (11)
9. March (25)
10. February (17)
11. January (22)
2. 2018 (219)
1. December (2)
2. November (32)
3. October (27)
4. September (26)
5. August (4)
6. July (9)
7. June (13)
8. May (17)
9. April (13)
10. March (28)
11. February (27)
12. January (21)
3. 2017 (281)
4. 2016 (335)
5. 2015 (209)
6. 2014 (235)
7. 2013 (251)
8. 2012 (135)


text:
# Glossary

## Newick format

Newick format offers a way of representing binary trees in text form by using parentheses and commas. Here's an example of a straightforward rooted tree.

According to this tree, taxa $a$ and $b$ appear to be closely related, as do $c$ and $d$. More formally, $a$ and $b$ are said to be neighbors since they are adjacent to the same internal node of the tree. We label this internal node using the syntax $(a, b)$, then delete $a$ and $b$ and iterate the process. In the end, the above tree can be represented using Newick format as $((a,b),(c,d))$.

The above rooted tree, which is more complex, can be represented by the following steps:

• consolidate $a$ and $c$ into $(a, c)$;
• consolidate $(a, c)$ and $b$ into $((a, c), b)$;
• consolidate $d$ and $e$ into $(d, e)$;
• consolidate $f$ and $g$ into $(f, g)$;
• consolidate $(d, e)$ and $(f, g)$ into $((d, e), (f, g))$;
• consolidate the final two leaf nodes, $((a, c), b)$ and $((d, e), (f, g))$, into the tree's final Newick representation $(((a, c), b), ((d, e), (f, g)))$.

For unrooted trees, we may arbitrarily construct a root on any one of the tree's edges. root picked arbitrarily at one of the internal edges. As a result, we will wind up with different possible Newick formats. For example, the tree above tree can be represented as $(((a,b),c),(d,e))$ or $((a,b),(c,(d,e)))$; try to find another representation of it.