text:
Definitions

# Simulated annealing

Simulated annealing (SA) is a generic probabilistic meta-algorithm for the global optimization problem, namely locating a good approximation to the global optimum of a given function in a large search space. It is often used when the search space is discrete (e.g., all tours that visit a given set of cities). For certain problems, simulated annealing may be more effective than exhaustive enumeration — provided that the goal is merely to find an acceptably good solution in a fixed amount of time, rather than the best possible solution.

The name and inspiration come from annealing in metallurgy, a technique involving heating and controlled cooling of a material to increase the size of its crystals and reduce their defects. The heat causes the atoms to become unstuck from their initial positions (a local minimum of the internal energy) and wander randomly through states of higher energy; the slow cooling gives them more chances of finding configurations with lower internal energy than the initial one.

By analogy with this physical process, each step of the SA algorithm replaces the current solution by a random "nearby" solution, chosen with a probability that depends on the difference between the corresponding function values and on a global parameter T (called the temperature), that is gradually decreased during the process. The dependency is such that the current solution changes almost randomly when T is large, but increasingly "downhill" as T goes to zero. The allowance for "uphill" moves saves the method from becoming stuck at local minima—which are the bane of greedier methods.

The method was independently described by S. Kirkpatrick, C. D. Gelatt and M. P. Vecchi in 1983, and by V. Černý in 1985. The method is an adaptation of the Metropolis-Hastings algorithm, a Monte Carlo method to generate sample states of a thermodynamic system, invented by N. Metropolis et al in 1953.

## Overview

In the simulated annealing (SA) method, each point s of the search space is analogous to a state of some physical system, and the function E(s) to be minimized is analogous to the internal energy of the system in that state. The goal is to bring the system, from an arbitrary initial state, to a state with the minimum possible energy.

### The basic iteration

At each step, the SA heuristic considers some neighbour s' of the current state s, and probabilistically decides between moving the system to state s' or staying in state s. The probabilities are chosen so that the system ultimately tends to move to states of lower energy. Typically this step is repeated until the system reaches a state that is good enough for the application, or until a given computation budget has been exhausted.

### The neighbours of a state

The neighbours of each state (the candidate moves) are specified by the user, usually in an application-specific way. For example, in the traveling salesman problem, each state is typically defined as a particular tour (a permutation of the cities to be visited); and one could define the neighbours of a tour as those tours that can be obtained from it by exchanging any pair of consecutive cities.

### Acceptance probabilities

The probability of making the transition from the current state $s$ to a candidate new state $s\text{'}$ is specified by an acceptance probability function $P\left(e, e\text{'}, T\right)$, that depends on the energies $e = E\left(s\right)$ and $e\text{'} = E\left(s\text{'}\right)$ of the two states, and on a global time-varying parameter $T$ called the temperature.

One essential requirement for the probability function $P$ is that it must be nonzero when $e\text{'} > e$, meaning that the system may move to the new state even when it is worse (has a higher energy) than the current one. It is this feature that prevents the method from becoming stuck in a local minimum—a state that is worse than the global minimum, yet better than any of its neighbours.

On the other hand, when $T$ goes to zero, the probability $P\left(e, e\text{'}, T\right)$ must tend to zero if $e\text{'} > e$, and to a positive value if $e\text{'} < e$. That way, for sufficiently small values of $T$, the system will increasingly favor moves that go "downhill" (to lower energy values), and avoid those that go "uphill". In particular, when $T$ becomes 0, the procedure will reduce to the greedy algorithm—which makes the move only if it goes downhill.

In the original description of SA, the probability $P\left(e, e\text{'}, T\right)$ was defined as 1 when $e\text{'} < e$ — i.e., the procedure always moved downhill when it found a way to do so, irrespective of the temperature. Many descriptions and implementations of SA still take this condition as part of the method's definition. However, this condition is not essential for the method to work, and one may argue that it is both counterproductive and contrary to its spirit.

The $P$ function is usually chosen so that the probability of accepting a move decreases when the difference $e\text{'}-e$ increases—that is, small uphill moves are more likely than large ones. However, this requirement is not strictly necessary, provided that the above requirements are met.

Given these properties, the evolution of the state $s$ depends crucially on the temperature $T$. Roughly speaking, the evolution of $s$ is sensitive to coarser energy variations when $T$ is large, and to finer variations when $T$ is small.

### The annealing schedule

Another essential feature of the SA method is that the temperature is gradually reduced as the simulation proceeds. Initially, $T$ is set to a high value (or infinity), and it is decreased at each step according to some annealing schedule—which may be specified by the user, but must end with $T=0$ towards the end of the allotted time budget. In this way, the system is expected to wander initially towards a broad region of the search space containing good solutions, ignoring small features of the energy function; then drift towards low-energy regions that become narrower and narrower; and finally move downhill according to the steepest descent heuristic.

 Example illustrating the effect of cooling schedule on the performance of simulated annealing. The problem is to rearrange the pixels of an image so as to minimize a certain potential energy function, which causes similar colours to attract at short range and repel at slightly larger distance. The elementary moves swap two adjacent pixels. The images were obtained with fast cooling schedule (left) and slow cooling schedule (right), producing results similar to amorphous and crystalline solids, respectively.

It can be shown that for any given finite problem, the probability that the simulated annealing algorithm terminates with the global optimal solution approaches 1 as the annealing schedule is extended. This theoretical result, however, is not particularly helpful, since the time required to ensure a significant probability of success will usually exceed the time required for a complete search of the solution space.

## Pseudocode

The following pseudocode implements the simulated annealing heuristic, as described above, starting from state s0 and continuing to a maximum of kmax steps or until a state with energy emax or less is found. The call neighbour(s) should generate a randomly chosen neighbour of a given state s; the call random() should return a random value in the range $\left[0, 1\right)$. The annealing schedule is defined by the call temp(r), which should yield the temperature to use, given the fraction r of the time budget that has been expended so far.

``` ```

``s := s0; e := E(s)                           // Initial state, energy.``
``` sb := s; eb := e // Initial "best" solution k := 0 // Energy evaluation count. while k < kmax and e > emax // While time remains & not good enough: sn := neighbour(s) // Pick some neighbour. en := E(sn) // Compute its energy. if en < eb then // Is this a new best? sb := sn; eb := en // Yes, save it. if P(e, en, temp(k/kmax)) > random() then // Should we move to it? s := sn; e := en // Yes, change state. k := k + 1 // One more evaluation done return sb // Return the best solution found. ```

Actually, the "pure" SA algorithm does not keep track of the best solution found so far: it does not use the variables sb and eb, it lacks the first if inside the loop, and, at the end, it returns the current state s instead of sb. While saving the best state is a standard optimization, that can be used in any metaheuristic, it breaks the analogy with physical annealing — since a physical system can "store" a single state only.

Saving the best state is not necessarily an improvement, since one may have to specify a smaller kmax in order to compensate for the higher cost per iteration. However, the step sb := sn happens only on a small fraction of the moves. Therefore, the optimization is usually worthwhile, even when state-copying is an expensive operation.

## Selecting the parameters

In order to apply the SA method to a specific problem, one must specify the following parameters: the state space, the energy (goal) function E(), the candidate generator procedure neighbour(), the acceptance probability function P(), and the annealing schedule temp(). These choices can have a significant impact on the method's effectiveness. Unfortunately, there are no choices of these parameters that will be good for all problems, and there is no general way to find the best choices for a given problem. The following sections give some general guidelines.

### Diameter of the search graph

Simulated annealing may be modeled as a random walk on a search graph, whose vertices are all possible states, and whose edges are the candidate moves. An essential requirement for the neighbour() function is that it must provide a sufficiently short path on this graph from the initial state to any state which may be the global optimum. (In other words, the diameter of the search graph must be small.) In the traveling salesman example above, for instance, the search space for $n=20$ cities has $n! = 2 432 902 008 176 640 000$ (2.5 quintillion) states; yet the neighbour generator function that swaps two consecutive cities can get from any state (tour) to any other state in $n\left(n-1\right)/2 = 190$ steps.

### Transition probabilities

For each edge $\left(s,s\text{'}\right)$ of the search graph, one defines a transition probability, which is the probability that the SA algorithm will move to state $s\text{'}$ when its current state is $s$. This probability depends on the current temperature as specified by temp(), by the order in which the candidate moves are generated by the neighbour() function, and by the acceptance probability function P(). (Note that the transition probability is not simply $P\left(e, e\text{'}, T\right)$, because the candidates are tested serially.)

### Acceptance probabilities

The specification of neighbour(), P(), and temp() is partially redundant. In practice, it's common to use the same acceptance function P() for many problems, and adjust the other two functions according to the specific problem.

In the formulation of the method by Kirkpatrick et al., the acceptance probability $P\left(e, e\text{'}, T\right)$ was defined as 1 if $e\text{'} < e$, and $exp\left(\left(e-e\text{'}\right)/T\right)$ otherwise. This formula corresponds to the Metropolis-Hastings algorithm, in the case where the proposal distribution of Metropolis-Hastings is symmetric. However, this acceptance probability is often used for simulated annealing even when the neighbour() function, which is analogous to the proposal distribution in Metropolis-Hastings, is not symmetric, or not probabilistic at all. Individual transitions of the simulated annealing algorithm do not correspond to the short-term evolution of a physical system, but rather the long-term distribution over states of the algorithm at a particular temperature corresponds to the probability distribution over states of a physical system at a particular temperature.

### Efficient candidate generation

When choosing the candidate generator neighbour(), one must consider that after a few iterations of the SA algorithm, the current state is expected to have much lower energy than a random state. Therefore, as a general rule, one should skew the generator towards candidate moves where the energy of the destination state $s\text{'}$ is likely to be similar to that of the current state. This heuristic (which is the main principle of the Metropolis-Hastings algorithm) tends to exclude "very good" candidate moves as well as "very bad" ones; however, the latter are much more common than the former, so the heuristic is generally quite effective.

In the traveling salesman problem above, for example, swapping two consecutive cities in a low-energy tour is expected to have a modest effect on its energy (length); whereas swapping two arbitrary cities is far more likely to increase its length than to decrease it. Thus, the consecutive-swap neighbour generator is expected to perform better than the arbitrary-swap one, even though the latter could provide a somewhat shorter path to the optimum (with $n-1$ swaps, instead of $n\left(n-1\right)/2$).

A more precise statement of the heuristic is that one should try first candidate states $s\text{'}$ for which $P\left(E\left(s\right), E\left(s\text{'}\right), T\right)$ is large. For the "standard" acceptance function $P$ above, it means that $E\left(s\text{'}\right) - E\left(s\right)$ is on the order of $T$ or less. Thus, in the traveling salesman example above, one could use a neighbour() function that swaps two random cities, where the probability of choosing a city pair vanishes as their distance increases beyond $T$.

### Barrier avoidance

When choosing the candidate generator neighbour() one must also try to reduce the number of "deep" local minima — states (or sets of connected states) that have much lower energy than all its neighbouring states. Such "closed catchment basins" of the energy function may trap the SA algorithm with high probability (roughly proportional to the number of states in the basin) and for a very long time (roughly exponential on the energy difference between the surrounding state and the bottom of the basin).

As a rule, it is impossible to design a candidate generator that will satisfy this goal and also prioritize candidates with similar energy. On the other hand, one can often vastly improve the efficiency of SA by relatively simple changes to the generator. In the traveling salesman problem, for instance, it is not hard to exhibit two tours $A$, $B$, with nearly equal lengths, such that (0) $A$ is optimal, (1) every sequence of city-pair swaps that converts $A$ to $B$ goes through tours that are much longer than both, and (2) $A$ can be transformed into $B$ by flipping (reversing the order of) a set of consecutive cities. In this example, $A$ and $B$ lie in different "deep basins" if the generator performs only random pair-swaps; but they will be in the same basin if the generator performs random segment-flips.

### Cooling schedule

The physical analogy that is used to justify SA assumes that the cooling rate is low enough for the probability distribution of the current state to be near thermodynamic equilibrium at all times. Unfortunately, the relaxation time—the time one must wait for the equilibrium to be restored after a change in temperature—strongly depends on the "topography" of the energy function and on the current temperature. In the SA algorithm, the relaxation time also depends on the candidate generator, in a very complicated way. Note that all these parameters are usually provided as black box functions to the SA algorithm.

Therefore, in practice the ideal cooling rate cannot be determined beforehand, and should be empirically adjusted for each problem. The variant of SA known as thermodynamic simulated annealing tries to avoid this problem by dispensing with the cooling schedule, and instead automatically adjusting the temperature at each step based on the energy difference between the two states, according to the laws of thermodynamics.

## Restarts

Sometimes it is better to move back to a solution that was significantly better rather than always moving from the current state. This is called restarting. To do this we set `s` and `e` to `sb` and `eb` and perhaps restart the annealing schedule. The decision to restart could be based on a fixed number of steps, or based on the current energy being too high from the best energy so far.

## Related methods

• Quantum annealing uses "quantum fluctuations" instead of thermal fluctuations get through high but thin barriers in the target function.
• Stochastic tunneling attempts to overcome the increasing difficulty simulated annealing runs have in escaping from local minima as the temperature decreases, by 'tunneling' through barriers.
• Tabu search normally moves to neighbouring states of lower energy, but will take uphill moves when it finds itself stuck in a local minimum; and avoids cycles by keeping a "taboo list" of solutions already seen.
• Stochastic gradient descent runs many greedy searches from random initial locations.
• Genetic algorithms maintain a pool of solutions rather than just one. New candidate solutions are generated not only by "mutation" (as in SA), but also by "combination" of two solutions from the pool. Probabilistic criteria, similar to those used in SA, are used to select the candidates for mutation or combination, and for discarding excess solutions from the pool.
• Ant colony optimization (ACO) uses many ants (or agents) to traverse the solution space and find locally productive areas.
• The cross-entropy method (CE) generates candidates solutions via a parameterized probability distribution. The parameters are updated via cross-entropy minimization, so as to generate better samples in the next iteration.
• Harmony search mimics musicians in improvisation process where each musician plays a note for finding a best harmony all together.
• Stochastic optimization is an umbrella set of methods that includes simulated annealing and numerous other approaches.

## References

• S. Kirkpatrick and C. D. Gelatt and M. P. Vecchi, Optimization by Simulated Annealing, Science, Vol 220, Number 4598, pages 671-680, 1983. http://www.cs.virginia.edu/cs432/documents/sa-1983.pdf or http://citeseer.ist.psu.edu/kirkpatrick83optimization.html .
• V. Cerny, A thermodynamical approach to the travelling salesman problem: an efficient simulation algorithm. Journal of Optimization Theory and Applications, 45:41-51, 1985
• N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller. "Equations of State Calculations by Fast Computing Machines". Journal of Chemical Physics, 21(6):1087-1092, 1953.
• A. Das and B. K. Chakrabarti (Eds.), Quantum Annealing and Related Optimization Methods, Lecture Note in Physics, Vol. 679, Springer, Heidelberg (2005)
• E. Weinberger, Correlated and Uncorrelated Fitness Landscapes and How to Tell the Difference, Biological Cybernetics, 63, No. 5, 325-336 (1990).
• J. De Vicente, J. Lanchares, R. Hermida, "Placement by Thermodynamic Simulated Annealing”, Physics Letters A,Vol. 317, Issue 5-6, pp.415-423, 2003.


text:
# Regulatory Reforms - Can we make bankers more like pilots?

The banking industry today faces many challenges, and the need for strong internal governance and control has never been more important. Controls, rules and limits within a bank – and in particular, the interactions of those controls, rules and limits – do more than just limit risk; they also create incentives.

The following speech by Wayne Byres discusses about why incentives matter. But in diverse and complex organisations such as today’s banks, ensuring that incentives are working in the right direction is easier said than done. Effective regulation and supervision is about helping to ensure those incentives are appropriately aligned, for the benefit not just of the bank but also for society as a whole.

[gview file="http://www.bis.org/speeches/sp121024.pdf" save="1"]

# R Programming Bundle: 25% OFF

Get our R Programming - Data Science for Finance Bundle for just $29$39.
Get it now for just \$29

format: plain


text:
# teaching machines

## CS 268: Lecture 8 – Position

Dear students:

Prior to 1998, CSS was mostly a vehicle for factoring the coloring, sizing, and alignment attributes out of HTML into a reusable stylesheet. But in 1998, CSS gained a new property: position. This opened the door to thinking of CSS as a tool for layout. Eventually we got robust layout systems like Flexbox and Grid, but position is still useful to achieve certain effects. We begin exploring those effects today.

But first, let’s complete an exercise!

### Activity

If the children of a flex container have all been shrunk to their utmost and there still is not enough room for them, they can be made to wrap along the cross axis using the flex-wrap property. Valid values include nowrap, wrap, and wrap-reverse. We’ll see these at work in Flexercise.

With exactly one neighbor, claim your task on Crowdsource. (Everyone will be assigned task 1.) Then recreate the following page:

Apply styling to the following starter code to match the page above:

<!DOCTYPE html>

<title>Gallery</title>

Notice how the content is centered. We have used Flexbox to achieve vertical and horizontal centering. To achieve horizontal centering only, we have another option:

• Set the width of the content to be centered.
• Set the horizontal margin of the content to be centered to auto.

Notice also the shadowing, which was achieved with the box-shadow property.

Submit your solution on Crowdsource.

### Position

Normally, our boxes of content are positioned within the flow of the document. An inline element appears to the right of its predecessor box if there’s room, or on the next line if there’s not. A block element appears below its predecessor box and above its successor box. What is true of a page where flow strictly follows this pattern of succession? It’s flat and rigid.

More interesting layouts can be achieved if we allow elements to break out of this flow. We can do this with the position property in CSS. By default, our elements are all have static positions. There are three other widely supported options: relative, absolute, and fixed. Each lets us anchor an element to some box on our page. Where they differ is in what that anchor box is. Once the anchor is identified, we position and size an element using offsets from that anchor.

Let’s examine each of these position values in turn.

#### Relative

An element with a relative position is anchored against the box it would have if it had a static position. But the element’s box can be shaped relative to the static box with one or more of the four edge properties: left, right, top, and bottom. The box that is shaped appears to be the margin-box, though I’ve only deduced this from the browsers I have handy.

For example, we can shift an element down from its normal position with this style:

.shifted {
position: relative;
top: 15px;
}

Note that the shifted content does cause the neighboring content to reflow. This can lead to stacking.

#### Fixed

An element with a fixed position is anchored against the viewport of the browser. This kind of positioning is often used to make what is called a heads-up display in games, which are UI elements like mini-maps and health bars that are laid over the screen and don’t move with the rest of the content.

For example, we can make a message appear at the bottom of the screen alerting a user that a site uses cookies, like this:

Unlike relative positioning, the element does not leave a shadow, statically-positioned box consuming space in the flow.

#### Absolute

An element with an absolute position is anchored against the box of its nearest non-static ancestor. This kind of positioning is often used to make child elements appear in a certain relationship to their parents.

Unlike relative positioning, the element does not leave a shadow, statically-positioned box consuming space in the flow.

We’re going to devote our next lecture to effects we can achieve with absolute positions.

#### Overconstraining

If all three of width, left, and height are defined, then right is ignored when the flow is left-to-right. If all three of height, top, and bottom are defined, then bottom is ignored.

For relative positions, if both top and bottom are set, bottom is ignored, regardless of whether or not height is defined. But I don’t know why.

#### z-index

These three non-static positions can leading to situations were elements stack on top of one another. The ordering of stacked elements is determined by the z-index property. Elements with a larger z-index appear on top of elements with a smaller z-index. If two elements have the same z-index, whichever appeared later in the document is on top.

### TODO

Here’s your TODO list for next time:

• Play CSS Diner to learn some new CSS selectors that we haven’t discussed.

See you next time.

Sincerely,

P.S. It’s time for a haiku! This one’s more personal. It’s an homage to one of sons, who uses a relative positioning scheme.

A chair’s under him
I wouldn’t say he’s sitting
There’s too much offset

P.P.S. Here’s the code we wrote together in class…

#### gallery.html

<!DOCTYPE html>

<title>Gallery</title>

#### bernie.html

<!DOCTYPE html>

<title>...</title>


text:
# Exact algorithm for NAE-3SAT

The NAE-3SAT problem is to determine whether a given 3CNF formula has a satisfying assignment that gives each clause at least one false (and at least one true) literal. The problem is NP-complete. One can reduce 3SAT to it pretty easily so that the number of variables becomes roughly the number of clauses in the original instance, and so because of sparsification, under ETH you won't be able to get a subexponential time algorithm. (Also, it's known that NAE-SAT for unbounded length clauses requires $2^n$ time under Strong ETH.)

My question is, what is the best exponential time running time in terms of the number of variables?

There is a trivial reduction to 3SAT that does not increase the number of variables $n$ so that the fastest 3SAT algorithm has the same running time on NAE-3SAT instances (the best known is by Hertli, $1.308^n$). Is there a faster known algorithm for NAE-3SAT than the one for 3SAT? (Or is there a reduction that shows that the two problems are equivalent with respect to exact algorithm running times?)

What about monotone NAE-3SAT? Here there are no negated variables. It's known that this problem is also NP-complete. What's the fastest algorithm for it?

• I suppose the following paper is your source of the statement "NAE-SAT for unbounded length clauses requires $2^n$ time under Strong ETH": Cygan et al. On problems as hard as CNF-SAT. Please note down any other highly relevant papers. A starting point could be of immense help to seach in the literature. Mar 28 at 5:47

S. Cliff Liu [1] reports a deterministic algorithm for NAE-$$k$$-SAT that is faster than the best deterministic algorithm for $$k$$-SAT on all $$k\geq 3$$.

(Quote from Abstract)

Our main results include:
• A deterministic algorithm for NAE-$$k$$-SAT that is faster than the best deterministic algorithm for $$k$$-SAT on all $$k\geq 3$$. Previously, no NAE-$$k$$-SAT algorithm is known to be faster than $$k$$-SAT algorithms. For $$k = 3$$, we achieve an upper bound of $$1.326^n$$. The corresponding bound for 3-SAT is $$1.328^n$$.
...
Our finding sheds new light on the following question: Is NAE-k-SAT easier than k-SAT? The answer might be affirmative at least on solving the problems exactly and deterministically.

See also [2] for improvements of run-time of $$k$$-SAT.

## References

[1] Liu, S. Cliff. "The Curse and Blessing of Not-All-Equal in k-Satisfiability." arXiv preprint arXiv:1809.04312 (2018).

[2] Liu, S. Cliff. "Chain, generalization of covering code, and deterministic algorithm for k-sat." arXiv preprint arXiv:1804.07901 (2018).


text:
# cf.DimensionCoordinate.period¶

DimensionCoordinate.period(*value)[source]

Set the period for cyclic coordinates.

Parameters: value: data-like or None, optional The period. The absolute value is used. {+data-like-scalar} TODO out: cf.Data or None The period prior to the change, or the current period if no value was specified. In either case, None is returned if the period had not been set previously.

Examples:

>>> print(c.period())
None
>>> c.Units
<Units: degrees_east>
>>> print(c.period(360))
None
>>> c.period()
<CF Data(): 360.0 'degrees_east'>
>>> import math


text:
# Stokes Theorem cone oriented downwards

1. Jul 28, 2012

### gtfitzpatrick

1. The problem statement, all variables and given/known data

Verify stokes theorem where F(xyz) = -yi+xj-2k and s is the cone $z^2 = x^2 + y^2$ , 0≤ Z ≤ 4 oriented downwards

2. Relevant equations

$\oint_{c} F.dr = \int\int_{s} (curlF).dS$

3. The attempt at a solution

Firstly the image of the widest part of the cone on the xy plane is the circle ofradius 4
we parametize this circle using the parameters x=4cost , y=4sint , z=0 and also dx=-4sintdt , dy=4costdt

so we get
$\oint_{c} F.dr = \int^{2\pi}_{0} -ydx + xdy -2dz$
$\oint_{c} F.dr = \int^{2\pi}_{0} (-4sint)(-4sint) + (4cost)(4cost)dt$
$= 16 \int^{2\pi}_{0} (sin^2 t) + (cos^2 t)dt$
$= 16 \int^{2\pi}_{0} dt = 16(2\pi- 0) = 32\pi$

next curlF = 2k

so $= \int\int_{s} (curlF).dS = \int\int_{D} (2).dA$
=$2\int^{4}_{0}\int^{2\pi}_{0} rdtdr = 2\int^{4}_{0} 2\pi r dr$
=$2(\pi r^2)^{4}_{0} = 2(16\pi) =32\pi$

ok so i got the same answer both ways which verifies Stoke but the thing im not sure about is the question says the cone is oriented downward so should i have reversed the limits as in 0≥r≥-4 and because it is downward orientated should i have put a minus infront of each side? but then the would have canceled each other out anyways?

2. Jul 28, 2012

### LCKurtz

First mistake. The bounding circle of the cone is not z = 0, it is z = 4, up on the cone.

Nitpick, should still be $\int_C$ for the second integral since it isn't expressed in terms of dt yet.

More than a nitpick here. At this point you must think about orientation. If the normal is downward and you use the right hand rule, what direction do you go around the circle? And the limits must reflect that.

That's a vector $d\vec S$ and you must choose it to be oriented correctly. Is it? How do you calculate $d\vec S$ and check its orientation?

No, not that way. See above.

3. Jul 29, 2012

### gtfitzpatrick

Hi LCkurtz,
If the normal is pointed downward and i use the righthand rule i am going clockwise around the circle. So my limits reflecting this should be 2∏ → 0? but then if i use these limits my answer becomes $-32\pi$?? but it should be just $32\pi$ ???

4. Jul 29, 2012

### LCKurtz

You didn't respond about the orientation of $d\vec S\$that I highlighted in red. If you had, you might not have had to ask that last question. In your equation $$\int\int_{s} (curlF).dS = \int\int_{D} (2).dA$$you haven't made clear how you got the orientation of $d\vec S$ or how you got from $d\vec S$ to $dA$. You have a dot in front of dA but you can't dot scalars. How do you know both sides don't give you $-32\pi$? (Disclaimer: I haven't yet checked the calculation details since there may be no problem once you understand proper orientation.)

5. Jul 29, 2012

### gtfitzpatrick

Sorry the dot shouldnt be there.
Using the right hand rule my limits should be -c or 2∏ to 0. dS is calculated by the normal or $\sqrt{(\frac{dx}{dt})^2 + (\frac{dy}{dt})^2 + (\frac{dz}{dt})^2 } dA$ where x=rcost,y=rsint,z=r

Thanks for the help LCkurtz

6. Jul 29, 2012

### LCKurtz

Yes, the limits for the line integral should be $-C$ or $2\pi$ to $0$. But we are talking about the surface integral now. And yes, $d\vec S$ is the normal surface area vector. You must calculate it and make sure it agrees with the orientation. That square root formula you gave is not a vector and is not the correct formula.

Either look in your text for the correct formula for $d\vec S$ or read my post #13 here:

7. Jul 30, 2012

### gtfitzpatrick

Thanks that clears up alot for me.
Thanks a million for all your help so far but can i ask another question?
if i was dealing with a hemisphere that was orientated in the postive y or x axis would $\hat{n} dS = \left\langle -\frac{ \partial f}{\partial x}, 1 ,-\frac{ \partial f}{\partial z}, \right\rangle dz dx$ say if it was oriented in the positive y direction? and set my limits for z and x accordingly, well i'd have to paramaterize it but is my thinking right?

8. Jul 30, 2012

### gtfitzpatrick

Or if i use $d\vec S =\hat n dS =\pm\frac{\vec R_u\times\vec R_v}{|\vec R_u\times\vec R_v|}|\vec R_u\times \vec R_v|dudv=\pm{\vec R_u\times\vec R_v}dudv$ it doesnt matter which axis it is orientated towards. i just need to choose my ± accordingly?

9. Jul 30, 2012

### LCKurtz

Before changing the subject in this thread, how about posting how you did the surface integral, showing how you got from $d\vec S$ to $dA$. Did you get $-32\pi$ like you got in the line integral? Then I will know if you really understand it.

10. Jul 30, 2012

### LCKurtz

You don't have to learn any such formulas. If $y = f(x,z)$ your parameterization is$$\vec R(x,z) =\langle x, f(x,z), z)\rangle$$So $\vec R_x = \langle 1, f_x,0\rangle$ and $\vec R_z = \langle 0, f_z, 1\rangle$. The cross product is$$\vec R_x \times \vec R_z = \langle f_x, -1, f_z\rangle$$
Notice that this is opposite the vector you have, but, yes, you would think about the sign of the $y$ component when determining the orientation, and choose the $\pm$ sign accordingly when calculating $d\vec S$.

11. Aug 1, 2012

### gtfitzpatrick

Hi LCkurtz,I got from dS⃗ to dA

dS = $\sqrt{(\frac{dz}{dx})^2 + (\frac{dz}{dy})^2 + 1 } dxdy$

x=rcost,y=rsint and sub these values in i get rdrdt or rdA?

thats how i did it.
but im still abit confused, i need to look at again.

12. Aug 1, 2012

### LCKurtz

So you have a parameterization$$\vec R(r,t) = \langle r\cos t, r\sin t, r\rangle$$

That argument still doesn't address the orientation because what you are calling dS is not a vector.

I have highlighted in red some of your earlier observations. You know curl F and you have a parameterization. Why don't you actually try using the formula I gave you in that post #13 reference I gave you in post #6 above? Now your parameters are r and t instead of u and v or x and y. And don't forget to think about the orientation choice of $\pm$ when you do it.

Last edited: Aug 1, 2012
13. Aug 2, 2012

### gtfitzpatrick

right so $\vec R(r,t) = \langle r\cos t, r\sin t, r\rangle$

$\vec R_{r} = \langle cos t, sin t, 1\rangle$ and $\vec R_{t} = \langle -rsin t, cos t, 0\rangle$

$\vec R_r \times \vec R_t = \langle -rcost, -rsint, r\rangle$ = $=\hat n dS$ im going to take this as + positive because although my cone is orientated downwards im still going around it it anticlock wise. like wise my limits will be from 0 to ∏
so then i get $\int^{4}_{0}\int^{2\pi}_{0} \langle 0,0,2\rangle .\langle -rcost, -rsint, r\rangle dtdr$
= $\int^{4}_{0}\int^{2\pi}_{0} 2r dtdr$ =32$\pi$

Last edited: Aug 2, 2012
14. Aug 2, 2012

### LCKurtz

No. You know that's wrong because the line integral gives $-32\pi$. Forget the anticlockwise and the line integral. You are doing the surface integral. The vector $\pm \langle -rcost, -rsint, r\rangle$ must agree with the orientation direction. The $z$ component is $r$ which is positive, hence upward. So you must choose the minus sign and use $-\langle -rcost, -rsint, r\rangle = \langle rcost, rsint, -r\rangle$ to agree with the given downward orientation. That will give you $-32\pi$, in agreement with the circuit integral. If you look back on how we have done the flux surface integral, you will see it really isn't that difficult. I hope you will study what we have done and learn from it.

15. Aug 2, 2012

### gtfitzpatrick

Thank you for your patience! i think i see now. I was thinking the normal was off the disk at the base of the cone in the positive direction but it should have been orientated in the same direction as the cone. I'll do a few more question. hopefully things will go a bit smoother.

Actually i have one more question if it was a cylinder instead of a cone that was open at both ends i would have to get the line integrals of both ends and add add them together to equal the surface integral. right?

16. Aug 2, 2012

### LCKurtz

That isn't a surface bounded by a closed curve. But imagine cutting it lengthwise to unroll it, but don't actually unroll it, just separate it very slightly. Now it is a surface bounded by a closed curve. If the cylinder was oriented outward you can use the right hand rule to determine the orientation of the boundary. Notice that the top and bottom circles are oriented in opposite directions and the slit is oriented once in each direction so can be ignored. So you can work such problems but you need to be careful with the orientation of things. But that shouldn't bother you much since you are now expert at orientation.


text:
Q. 2.7

A compression member constructed from steel pipe has an outside diameter of 90 mm and a cross sectional area of 1580 mm². What axial force P will cause the outside diameter to increase by 0.0094 mm? Given E = 200 GPa and 1/m = 0.3.

Verified Solution

First, we shall calculate longitudinal strain

$\varepsilon_{\text {long }}=\frac{\sigma}{E}$

But,      $\text { stress }=\sigma=\frac{P}{1580}$

$\begin{gathered} \varepsilon_{\text {long }}=\frac{P}{1580 \times 200 \times 10^3}=3.16 \times 10^{-9} P \\ \frac{1}{m}=\frac{\varepsilon_{\text {lat }}}{\varepsilon_{\text {long }}}, \quad \therefore \quad \varepsilon_{\text {lat }}=\frac{\text { change in diameter }}{\text { original diameter }}=\frac{0.0094}{90}=1.04 \times 10^{-4} \end{gathered}$

∴        $0.3=\frac{1.04 \times 10^{-4}}{3.16 \times 10^{-9} P}, \quad \therefore \quad P=110014.8 N =110.01 kN$


text:
# Changeset 3415 for src

Ignore:
Timestamp:
Jan 14, 2014, 3:09:15 PM (7 years ago)
Message:
• changes for proceedings of TACAS 2014
Location:
src/ASM/TACAS2013-policy
Files:
3 edited

Unmodified
Removed
• ## src/ASM/TACAS2013-policy/algorithm.tex

 r3393 return a smaller or equal solution. Experiments on the gcc 2.3.3 test suite of C programs have shown that on average, about 25 percent of jumps are encoded as short or absolute jumps by the algorithm. As not all instructions are jumps, this does not make for a large reduction in size, but it can make for a reduction in execution time: if jumps Experimenting with our algorithm on the test suite of C programs included with gcc 2.3.3 has shown that on average, about 25 percent of jumps are encoded as short or absolute jumps by the algorithm. As not all instructions are jumps, this does not make for a large reduction in size, but it can make for a reduction in execution time: if jumps are executed multiple times, for example in loops, the fact that short jumps take less cycles to execute than long jumps can have great effect. As for complexity, there are at most $2n$ iterations, with $n$ the number of As for complexity, there are at most $2n$ iterations, where $n$ is the number of branch instructions. Practical tests within the CerCo project on small to medium pieces of code have shown that in almost all cases, a fixed point is \If {$instr$ is a backward jump to $j$} \State $length \gets \mathrm{jump\_size}(pc,sigma_1(labels(j)))$ \Comment compute jump distance \ElsIf {$instr$ is a forward jump to $j$} \State $length \gets \mathrm{jump\_size}(pc,old\_sigma_1(labels(j))+added)$ \State $old\_length \gets \mathrm{old\_sigma_1}(ppc)$ \State $new\_length \gets \mathrm{max}(old\_length, length)$ \Comment length must never decrease \State $old\_size \gets \mathrm{old\_sigma_2}(ppc)$ \State $new\_size \gets \mathrm{instruction\_size}(instr,new\_length)$ \Comment compute size in bytes \State $new\_added \gets added+(new\_size-old\_size)$ \Comment keep track of total added bytes \State $new\_sigma \gets old\_sigma$ \State $new\_sigma_1(ppc+1) \gets pc+new\_size$ \State $new\_sigma_2(ppc) \gets new\_length$ \\ \State $new\_sigma_2(ppc) \gets new\_length$ \Comment update $\sigma$ \\ \Return $\langle new\_added, pc+new\_size, new\_sigma \rangle$ \EndFunction \item $instr$, the instruction currently under consideration; \item $ppc$, the pseudo-address of $instr$; \item $acc$, the fold accumulator, which contains $pc$ (the highest memory address reached so far), $added$ (the number of bytes added to the program size with respect to the previous iteration), and of course $sigma$, the \item $acc$, the fold accumulator, which contains $added$ (the number of bytes added to the program size with respect to the previous iteration), $pc$ (the highest memory address reached so far), and of course $sigma$, the $\sigma$ function under construction. \end{itemize}
• ## src/ASM/TACAS2013-policy/main.tex

 r3361 \usepackage{amsmath} \usepackage[british]{babel} \usepackage{caption} \usepackage{hyperref} \usepackage[utf8]{inputenc}
• ## src/ASM/TACAS2013-policy/proof.tex

 r3393 &&&&& (pc + \mathtt{instruction\_size}\ sigma\ ppc\ instruction < 2^{16}\ \vee\notag\\ &&&&& (\forall ppc'.ppc' < |instr\_list| \rightarrow ppc < ppc' \rightarrow \notag\\ &&&&& \mathbf{let}\ instruction' \equiv \mathtt{fetch\_pseudo\_instruction}\ instr\_list\ ppc'\ ppc\_ok'\ \mathbf{in} \notag\\ &&&&& \mathbf{let}\ instruction' \equiv \mathtt{fetch\_pseudo\_instruction}\ instr\_list\ ppc'\ \mathbf{in} \notag\\ &&&&&\ \mathtt{instruction\_size}\ sigma\ ppc'\ instruction' = 0)\ \wedge \notag\\ &&&&& pc + \mathtt{instruction\_size}\ sigma\ ppc\ instruction = 2^{16})
Note: See TracChangeset for help on using the changeset viewer.


text:
# Differential amplifier with current source load

Status
Not open for further replies.

#### sutapanaki

For example Razavi, "Design of analog CMOS Integrated Circuits" 2nd ed. page 80, eq.(3.116). Or Carusone, David Johns Ken Martin "Analog Integrated Circuit Design" 2nd ed. page 126, eq.(3.20)

I'm attaching here a short derivation of the source resistance for the common gate stage.

#### LvW

Thank you for the literature information and your derivation.
However, one short comment seems to be in order: For further signal processing, the output (drain node) of the most right transistor must be connected to the 2nd stage. This results, of course, in a drastic change of the impedance effective at this node and - hence - of the input impedance at the source node.

#### sutapanaki

Yes, what you're saying is correct, especially at high frequencies since the load from the next stage is mostly capacitive in most of the cases. The derivation in my previous post is for DC and low frequencies. At high frequencies, when the load capacitance starts to show up, the source impedance will converge to 1/gm and then at even higher frequencies there will be also a pole at the source.

Status
Not open for further replies.


text:
A duck, pursued by a fox, escapes to the center of a perfectly circular pond. The fox cannot swim, and the duck cannot take flight from the water. The fox is $k$ times faster than the duck. Assuming the fox pursues an optimal strategy, how can the duck reach the edge of the pond and fly away without being eaten?
First, find a strategy by which the duck can escape if $k=4$. Then, determine the smallest value of $k$ such that it is impossible for the duck to escape.


text:
Author's Note:  Let it be known that the material in this article is intended for recreational use only.  Attempting to record bands without prior approval is considered nasty, and could get you into a load of trouble.  You've been warned.

I'm sure that anybody fortunate enough to suffer through prime time television in the 70's will remember the show What's Happening!  The most memorable episode that comes to mind is the one where Rerun  sneaks a tape recorder into a Doobie Brothers concert.  Half way through their masterful performance, the recorder falls out of Rerun's pocket, the Doobs stop playing and then the moral of the story kicks in.  I don't really recall what the moral was, but I do remember that it was rather nasty and put recording shows in a bad light.  Much the same way that photographers take pictures at shows to chronicle an event, people who record shows, for the most part, do it for their own personal recreation.

THE RECORDER

Much like purchasing a 4 track, there's a multitude of options.

Sony WMD3  This model looks almost identical to a real Walkman.   And like its lower models, it has only a two head design, and the microphone doesn't really have a wide range.  Like all the other models that will be mentioned, the D3 has a recording level control knob.  Some fancier recorders have separate knobs for the left and right channel, but this just always seemed like a bit too much work for me.  Regarding the D3, I never had much luck with them, especially when it was so damned easy to accidentally bump the record switch off.  However, for many, this suits their needs just fine.  It costs around $300, and like the other machines mentioned, can be found at finer electronic stores. However, for my own peace of mind, I've purchased all of my recorders from J&R Music World in Maspeth, New York. [call 1-800-221-8180 to be sent a free catalog] I do this in large part due to the fact that they've got a very good return policy, a massive selection from inventory, and you can purchase an additional warranty to cover electronic foul-ups. Sony WMD6 If there is an Econoline van of the portable tape recording world, this would have to be it. It's sturdy, reliable and will last forever. The 3 head design is terrific and, providing that you're not a perfectionist like most, it will suit almost every person's needs. It has select switches for tape type and Dolby, line in/out ports and an attenuation switch. Apart from the Dolby, none of these features are on the wimpier D3. The attenuation switch is of particular importance. If the show is loud (e.g. requires earplugs), the -20dB attenuation helps a lot. For quieter, more acoustic-y shows (sans plugs) the 0dB switch will be the one to go with. The D6 is used more widely than any other recording device on the market. If Whitney, soundman for such bands as the Jesus Lizard, Blues Explosion and Jawbox, uses it to record their shows, that should give you a good idea of how well these machines travel and perform. Sony WMD7, WMD8 and all other digital recording devices Why get a D7? Because it's one better than a D6! Hoooo! But seriously, I've never been able to figure out the need to record live shows with digital equipment, so my personal knowledge on the subject is rather limited. I do know that with the ever increasing bootleg CD market, DATs are used almost exclusively as master sources. In addition, with regular analog tapes, quality can be lost in two to three generations from the master. With DAT, there's no loss. Big deal. A D7 (or any digital recorder for that matter) will cost almost double of a D6 or similar analog recorder. And if you make copies for fellow traders, you must also remember that you'll need a second DAT machine. Also, you have routine maintenance and additional equipment which could run your total set up cost well over$1,500!  For just a hobby, this seems a bit outrageous, so I've never considered entering the digital...


text:
[Tex/LaTex] Convert LaTeX to SVG online

conversiononlinesvg

What websites allow users to quickly and simply convert mathematical equations into SVG without having to install or run any software?

1. Type: $$E = mc^2$$


text:
Question

# The sides of a triangle are in A.P. and the greatest angle exceeds the least by $$90^0$$, prove that the sides are proportional to $$\sqrt7+1, \sqrt 7$$ and $$\sqrt7-1$$.

Solution

## Let A be the greatest and C the least angle of $$\triangle ABC.$$ It is given that a, b, c are in A.P. so that a+c= 2b                ............(1)Also A-C =$$90^0$$From (1), Sin A + sinC = 2 sin B= 2 sin (A+C)   ........(2)$$\therefore 2 sin\dfrac{A+C}{2} cos\dfrac{A-C}{2}=4 sin\dfrac{A+C}{2}cos \dfrac{A+C}{2}$$or $$cos\dfrac{A -C}{2}=2 cos \dfrac{A+C}{2}=2 sin{B}{2}$$$$\therefore 2 sin(B/2)= cos(90^0/2)=cos45^0 =1/\sqrt2$$,         by(2)Hence $$sin\dfrac{B}{2}=\dfrac{1}{2\sqrt 2}$$and $$cos\dfrac{B}{2}=\sqrt{(1-\dfrac{1}{8})}=\dfrac{\sqrt7}{2\sqrt(2)}$$$$\therefore sinB = 2 sin\dfrac{B}{2}cos\dfrac{B}{2}$$= $$2.\dfrac{1}{2\sqrt{(2)}}.\dfrac{\sqrt7}{2\sqrt{(2)}}=\dfrac{\sqrt7}{4}.$$$$\therefore sinA + sinC = 2sinB=\sqrt7/2$$            .....(3)  Also $$sin A - sin C = 2 cos [(A-C)/2]. sin[(A-C)/2] = 2 sin (B/2) sin45^0 = 2. \dfrac{1}{2\sqrt{(2)}}.\dfrac{1}{\sqrt{(2)}} =\dfrac{1}{2}$$ .....(4)Addding and subtracting (3) and (4), we get $$2sinA= (\sqrt 7+1)/2$$i.e.$$sin A= (\sqrt 7+1)/4$$and $$2 sin C= (\sqrt 7-1)/2$$i.e.$$sinC =(\sqrt 7-1)/4$$Hence a:b:c = sinA: sinB: sinC=$$\sqrt 7+ 1: \sqrt7: \sqrt7 -1$$Maths

Suggest Corrections

0

Similar questions
View More

People also searched for
View More


text:
# Signal Source Separation Using W-Net Architecture

This example shows how to separate two mixed signal sources using a deep learning network. Source separation is a common and complex signal processing problem that finds use in audio, vibration analysis, and biomedical applications. It consists of separating the signal components of a signal mixture when only the mixture is available.

An important source separation problem consists of discerning fetal and maternal electrocardiogram (ECG) signals present in noninvasive measurements taken on the abdominal area of a pregnant patient. This is an important problem because, if solved correctly, it can allow physicians to monitor the fetal ECG with minimum risk. Fetal cardiac monitoring and assessment during pregnancy are used for the early detection of fetal cardiac conditions.

This example uses simulated noninvasive abdominal ECG measurements on pregnant patients to illustrate how to solve the difficult problem of separating the fetal ECG and maternal ECG signals using a deep network. The source separation deep learning architecture used in this example is not limited to ECG signals and can be used in many other applications.

### FECGSYN Data Set

This example uses the FECGSYN PhysioNet data set [1], [2], which contains simulated adult and noninvasive fetal ECG signals. The data is generated using the FECGSYN simulator [3]. The simulator represents maternal and fetal hearts as punctual dipoles with different magnitudes and spatial positions. It obtains fetal–maternal mixtures by treating each abdominal signal and noise component as an individual source whose signal is propagated onto the observational points (electrodes). This database is able to provide separate waveform files for each signal source, making it ideal to test a source separation deep learning model.

The FECGSYN consists of simulated ECG signals corresponding to ten different subjects. For each subject, simulations produced a fetal ECG (`fECG`), a maternal ECG (`mECG`), and two noise sources, all sampled at a rate of 250 Hz for five minutes. The original data set repeats simulations five times for five different SNR levels, for 34 ECG channels or "electrodes", and for five different measurement scenarios or cases. In this example we use a subset of the data set and consider all ten subjects, a single channel (channel 19 from the original data set), four SNR levels (3, 6, 9, and 12 dB), and three different measurement cases labeled C0, C1, and C3. As mentioned before, the simulation was repeated over five iterations for each combination of subject, SNR value, and measurement case, yielding a total of 10 subjects × 4 SNRs × 3 cases × 5 iterations = 600 files. There are three different measurement cases:

• Case 0 (C0) — Baseline ECG signals

• Case 1 (C1) — Fetal movement + C0

• Case 3 (C3) — Signals with varying maternal and fetal heart rates + Noise from uterine contractions

The data set contains one MAT-file for each combination of subject, SNR level, iteration, and measurement case. The filenames use the format Ij_Ck.mat, where `j` is the iteration number (1 to 5) and `k` is the measurement case identifier (0, 1, 3). Each MAT-file contains these variables:

• `mECG` — Maternal ECG signal

• `fECG` — Fetal ECG signal

• `mECG_QRS` — QRS peak locations for the maternal ECG signal as annotated by an expert system

• `fECG_QRS` — QRS peak locations for the fetal ECG signal as annotated by an expert system

• `noise1` — First noise source

• `noise2` — Second noise source

All signals have been bandpass filtered into the frequency range from 5 Hz to 90 Hz.

The abdominal ECG signal (`aECG`) for each file is computed as the following mixture:

`$\mathrm{aECG}=\mathrm{mECG}+\mathrm{fECG}+\mathrm{noise1}+\mathrm{noise2}$`

The `mECG_QRS` and `fECG_QRS` variables contain QRS peak locations of the maternal and fetal ECG signals and can be used to validate the efficacy of a source separation algorithm to identify correct heartbeat locations in time.

This example uses the data from the first nine subjects to train a deep network and the data from the tenth subject to test the network performance. The training data size is about 1.15 GB, and the training of the deep learning network takes a few hours even when run on a GPU. If you want to skip downloading the training data and the training process, set the `trainNetworkFlag` flag to `false`. If the flag is set to `false`, the example downloads a pretrained network that can be used to perform source separation on the test data. The example always downloads the test data corresponding to subject 10.

`trainNetworkFlag = true;`

Download the train and test data sets using the `downloadSupportFile` function. The data will be unzipped to the `tempdir` directory. If you want the data at a different location, change `trainingDatasetFolder` and `testDatasetFolder` to the desired locations.

```if trainNetworkFlag % Download training data set trainingDatasetZipFile = matlab.internal.examples.downloadSupportFile('SPT','data/fetal-ecg-source-separation-trainingData.zip'); trainingDatasetFolder = fullfile(tempdir,'fetal-ecg-source-separation-trainingData'); if ~exist(trainingDatasetFolder,'dir') unzip(trainingDatasetZipFile,trainingDatasetFolder); end end % Download test data set testDatasetZipFile = matlab.internal.examples.downloadSupportFile('SPT','data/fetal-ecg-source-separation-testData.zip'); testDatasetFolder = fullfile(tempdir,'fetal-ecg-source-separation-testData'); if ~exist(testDatasetFolder,'dir') unzip(testDatasetZipFile,testDatasetFolder); end```

Create a signal datastore to access the files in the test data set. Specify the names of the variables that you want the datastore to read from each file.

```testDS = signalDatastore(testDatasetFolder,IncludeSubfolders=true, ... SignalVariableNames=["mECG" "fECG" "noise1" "noise2" "mECG_QRS" "fECG_QRS"]);```

Plot the first 2048 samples of the ECG signals for case C1 and SNR of 3 dB. Overlay the annotated QRS peaks for each signal.

```idx = contains(testDS.Files,fullfile("snr03dB","I1_C1.mat")); sds3dBC1 = subset(testDS,idx); data = preview(sds3dBC1); [mECG,fECG,noise1,noise2,mECG_QRS,fECG_QRS] = data{:}; % Abdominal ECG mixture aECG = mECG(1:2048) + fECG(1:2048) + noise1(1:2048) + noise2(1:2048); figure subplot(3,1,1) plot(aECG) xline(fECG_QRS(1:50),":",Color="#77AC30") xline(mECG_QRS(1:50),":",Color="#D95319") axis([0 2048 -0.6 1]) title("aECG (red = mECG QRS peaks, green = fECG QRS peaks)") subplot(3,1,2) plot(fECG(1:2048)) xline(fECG_QRS(1:50),":",Color="#77AC30") axis([0 2048 -0.6 1]) title("fECG") subplot(3,1,3) plot(mECG(1:2048)) xline(mECG_QRS(1:50),":",Color="#D95319") axis([0 2048 -0.6 1]) title("mECG")```

Notice the large difference in scale between the `mECG` and `fECG` signals.

### Prepare Training Data

This example uses the data from the first nine subjects to train a deep network and the data from the tenth subject to test the network performance. To train the network, each signal is broken into segments of 1024 samples for a total of 73 segments per signal. Set up a signal datastore that reads the ECG signals and noise realizations for subjects 1 to 9. Transform the datastore to obtain `aECG`, `mECG`, and `fECG` signal segments of length 1024 samples. Each read to `trainDS` returns 73 segments of length 1024 `aECG`, `mECG`, and `mECG` signals formatted using `CBT` (channel-batch-time) dimensions. In addition to segmenting the signals into 1024-sample segments, the transform function, `getECGSegments`, also normalizes each segment using the `rescale` function to bring the signal levels to between –1 and 1. The rescaled segments are then centered using their median value.

```segmentLength = 1024; if trainNetworkFlag trainDS = signalDatastore(trainingDatasetFolder,IncludeSubfolders=true,SignalVariableNames=["mECG" "fECG" "noise1" "noise2"]); trainDS = transform(trainDS,@(d,f)getECGSegments(d,segmentLength)); end```

To speed up training, read all the training data into memory so that the signal segmentation and normalization happens only once. If you have a Parallel Computing Toolbox™ license, use the `UseParallel` parameter so that the read operations are done in parallel. Create an array datastore to iterate through the training signal segments.

```if trainNetworkFlag trainData = readall(trainDS,UseParallel=true); trainDS = arrayDatastore(trainData,OutputType="same"); end```

### W-Net Architecture for Source Separation

This example uses a so-called W-Net architecture to perform source separation [4]. W-Net consists of two U-Net autoencoders [5] that have been modified to operate on 1-D signal inputs. A U-Net autoencoder is a deep network that encodes signal features reducing its size at each step and then decodes the features to recreate the original input signal. You can think of the encoder branch of the autoencoder as a feature extraction branch. The main idea of the W-Net architecture is to have one auto encoder to reproduce an `fECG` signal (`fECG` autoencoder) and another to reproduce an `mECG` signal (`mECG` autoencoder) when the input to the autoencoders is set to an `aECG` mixture. The connection between the two autoencoders happens in the encoding branches. You subtract the features obtained by the `mECG` autoencoder from the features obtained by the `fECG` autoencoder, effectively achieving separation of the `mECG` component from the `aECG` input and yielding the desired separated `fECG` signal. This figure shows the architecture in detail.

Following reference [4], for the ECG source separation problem at hand set the filter size of the 1D convolutional layers to 4 for the `fECG` side and 35 for the `mECG` side. The number of filters used at the input 1D convolutional layers, `N` in the figure above, is set to `16`. The input size, `P` in the figure above, has already been described as 1024. Create the W-Net network architecture using the createWNet function.

```if trainNetworkFlag filterSize_fECG = 4; filterSize_mECG = 35; numFilters_fECG = 16; numFilters_mECG = 16; lgraph = createWNet(segmentLength,filterSize_fECG,numFilters_fECG,filterSize_mECG,numFilters_mECG); wNet = dlnetwork(lgraph); end```

### Training Loop

You need a training loop to train the W-Net model because you need to define a loss that combines the losses of the `fECG` and `mECG` branches of the network. The `modelLoss` function computes the training loss as the weighted sum of the mean absolute deviation between the actual and predicted ECG signals:

`$\mathrm{loss}=\mathrm{fECGWeight}×\mathrm{mean}\left(|\mathrm{fECGactual}-\mathrm{fECGpredicted}|\right)+\mathrm{mECGWeight}×\mathrm{mean}\left(|\mathrm{mECGactual}-\mathrm{mECGpredicted}|\right);$`

Set `fECGWeight` to a value greater than `mECGWeight` to reflect the fact that the primary signals of interest are the fetal ECGs.

Use an Adam optimizer to update the network learnable parameters and specify an initial learn rate, a decay factor, the number of epochs, and the mini-batch size. The `minibatchqueue` outputs `miniBatchSize` batches of `aECG`, `mECG`, and `fECG` signal segments.

Due to the large size of the data set, the training process may take several hours. If your machine has a GPU and Parallel Computing Toolbox™, set the `useGPUflag` flag to `true` to speed up the training process.

```useGPUflag = true; if trainNetworkFlag NumEpochs = 100; miniBatchSize = 512; learnRate = 0.0005; decay = 0.25; mECGWeight = 0.25; fECGWeight = 0.75; mbqTrain = minibatchqueue(trainDS, 3, ... MiniBatchSize=miniBatchSize,... MiniBatchFormat={'CBT','CBT','CBT'}, ... MiniBatchFcn=@processMB, ... DispatchInBackground=true); if useGPUflag mbqTrain.OutputEnvironment = "gpu"; end % Initialize some training loop variables trailingAvg = []; trailingAvgSq = []; iteration = 0; lossByIteration = 0; minLoss = Inf; % Loop over epochs and store the lowest loss network, reshuffle the % mini-batch queue at each epoch for epoch = 1:NumEpochs reset(mbqTrain) shuffle(mbqTrain) % Loop over mini-batches while hasdata(mbqTrain) iteration = iteration + 1; % Get the next mini-batch [aECGbatch,mECGbatch,fECGbatch] = next(mbqTrain); % Evaluate the model gradients and loss [loss,gradients,state] = dlfeval(@modelLoss,wNet,aECGbatch, ... mECGbatch,fECGbatch,mECGWeight,fECGWeight); lossByIteration(iteration) = loss; % Update the network state wNet.State = state; % Update the network parameters using an Adam optimizer [wNet,trailingAvg,trailingAvgSq] = adamupdate(wNet,gradients, ... trailingAvg,trailingAvgSq,iteration,learnRate,decay); end if loss < minLoss minLoss = loss; bestModel = wNet; % Uncomment the line below to save the best model so far %save Model.mat wNet end end wNet = bestModel; % Plot the loss by iteration figure plot(1:iteration,mag2db(lossByIteration)) grid on title("Training Loss by Iteration") xlabel("Iteration") ylabel("Loss (dB)") axis tight end```

Load a pretrained model if `trainNetworkFlag` is `false`. The model file will be unzipped to the `tempdir` directory. If you want the model at a different location, change `modelFolder` to the desired value.

```if ~trainNetworkFlag % Download the pre-trained network modelZipFile = matlab.internal.examples.downloadSupportFile('SPT','data/fetal-ecg-source-separation-model.zip'); modelFolder = fullfile(tempdir,'fetal-ecg-source-separation-model'); if ~exist(modelFolder,'dir') unzip(modelZipFile,modelFolder); end modelFile = fullfile(modelFolder,'fetal-ecg-source-separation-model','Model.mat'); load(modelFile) end```

### Test Model

To test the trained network, use the previously created test datastore, `testDS`, that points to data from subject 10. This datastore reads the ECG data and the QRS peak location annotations so they can be used to validate the predicted `mECG` and `fECG` signals. As was done for the training datastore, transform the test datastore to get segmented and normalized `aECG`, `mECG`, and `fECG` signals.

`testDS = transform(testDS,@(d,f)getECGSegments(d,segmentLength));`

Call the predict method of the trained network to get separated `mECG` and `fECG` signals from an `aECG` input. Take for example iteration 3 of case C1 with and SNR of 9 dB. Estimate the fetal and maternal waveforms for that case as follows.

```idx = contains(string(testDS.UnderlyingDatastores{1}.Files),fullfile("snr09dB","I3_C1.mat")); ds = subset(testDS,idx); data = read(ds); mECG_QRS = data(:,4); fECG_QRS = data(:,5); [aECGbatch,mECGbatch,fECGbatch] = processMB(data(:,1),data(:,2),data(:,3)); % Move the aECGbatch into a dlarray and call the predict method of the % trained network to estimate the source signals dlaECG = dlarray(aECGbatch,"CBT"); [dlpred_fECG,dlpred_mECG] = predict(wNet,dlaECG); pred_fECG = squeeze(extractdata(dlpred_fECG))'; pred_mECG = squeeze(extractdata(dlpred_mECG))'; pred_fECG = pred_fECG(:); pred_mECG = pred_mECG(:);```

Plot a few samples of the predicted waveforms. Overlay the annotated true QRS peaks using dotted lines.

```figure subplot(2,1,1) plot(pred_fECG(1:2048)) xline([fECG_QRS{1}; fECG_QRS{2}],":k") title("Predicted fECG") axis([1 2048 -1.5 1]) subplot(2,1,2) plot(pred_mECG(1:2048)) xline([mECG_QRS{1}; mECG_QRS{2}],":k") title("Predicted mECG") axis([1 2048 -1 2])```

Plot predicted ECG signals for the case of high (12 dB) and low (3 dB) SNRs, for iteration 4 measurements, and for all three measurement cases using the `plotPredictedECGs` function. `N` and `M` can be set to plot segments `N` to `N`+`M` for the case at hand.

```% Plot segment 4 for each case N = 4; M = 1; plotPredictedECGs(wNet,testDS,"12","C0","I4",N,M)```

`plotPredictedECGs(wNet,testDS,"12","C1","I4",N,M)`

`plotPredictedECGs(wNet,testDS,"12","C3","I4",N,M)`

`plotPredictedECGs(wNet,testDS,"03","C0","I4",N,M)`

`plotPredictedECGs(wNet,testDS,"03","C1","I4",N,M)`

`plotPredictedECGs(wNet,testDS,"03","C3","I4",N,M)`

`fECG` signals have faster heart rates than `mECGs` so we show fewer `fECG` points just for better visualization. The dotted lines on the plots correspond to annotated ground truth QRS peak locations. Proper location of QRS peaks is as important as the estimation of the overall signal shape. QRS peak locations allow estimation of heart rate and conditions like arrhythmia. Proper peak location should be considered when evaluating the performance of the source separation procedure.

Recall that the main purpose of this network is to extract fetal ECG signals, which are the most difficult to obtain from the mixture. In the W-Net architecture the primary target is the one estimated by the left U-Net branch, which corresponds to the first output of the network built in this example using the `modelLoss` function. Overall, the network does a very good job in estimating QRS peak locations and waveform shapes for both high- and low-SNR cases and different measurement conditions.

There are extreme measurement cases where the combination of noise, fetal movements, and heart rate variations are too severe for the network. For example, plot the ECG estimates for an SNR of 6 dB, measurement case C3, and iteration 5. In this case, the network fails to predict an acceptable `fECG` waveform.

`plotPredictedECGs(wNet,testDS,"06","C3","I5",N,M)`

Plot the mean absolute deviation of the estimated `fECG` and `mECG` signals for all measurements of subject 10 using the `computeErrorsForAllCases` function.

`computeErrorsForAllCases(wNet,testDS)`

The errors do not decrease monotonically with SNR because of the variability of all the different combinations of noise, fetal movement, and heart-rate irregularities.

### Conclusion

This example implements a W-Net architecture suitable for source separation of a mixture of two signals. The example analyzes the performance of the network using synthetic signal mixtures comprised of fetal and maternal ECG waveforms. The example shows that, in most scenarios, the network does a good job separating ECG signals and estimating correct waveform shapes and QRS peak locations.

### References

[1] Goldberger, Ary L., Luis A. N. Amaral, Leon Glass, Jeffrey M. Hausdorff, Plamen Ch. Ivanov, Roger G. Mark, Joseph E. Mietus, George B. Moody, Chung-Kang Peng, and H. Eugene Stanley. “PhysioBank, PhysioToolkit, and PhysioNet.” Circulation 101, no. 23 (June 13, 2000): e215–20. https://doi.org/10.1161/01.CIR.101.23.e215.

[2] F. Andreotti, J. Behar, and G. D. Clifford. Fetal ECG Synthetic Database v1.0.0 (physionet.org), April 29, 2016, Version 1.0.0.

[3] F. Andreotti, J. Behar, S. Zaunseder, J. Oster, and G. D. Clifford. "An Open-Source Framework for Stress-Testing Non-Invasive Foetal ECG Extraction Algorithms." Physiological Measurement, Volume 37, Number 5, 2016.

[4] K. J. Lee and B. Lee, "End-to-End Deep Learning Architecture for Separating Maternal and Fetal ECGs Using W-Net," IEEE Access, Volume 10, pp. 39782-39788, 2022.

[5] O. Ronneberger, P. Fischer, and T. Brox. "U-Net: Convolutional Networks for Biomedical Image Segmentation", MICCAI, May 18, 2015.

### Appendix: Helper Functions

The functions listed in this section are only for use in this example. They may change or be removed in a future release.

#### getECGSegments

This function creates `aECG` mixtures from `mECG`, `fECG`, and noise signals. The function breaks the ECG signals into segments of length `segmentLength`. Each segment is normalized and reshaped to a `CBT` format with `C` and `B` equal to 1. When the input to the function contains QRS peak locations, the function breaks the locations according to the start and end index of each segment.

```function outputCell = getECGSegments(cellInput,segmentLength) mECG = cellInput{1}; fECG = cellInput{2}; noise1 = cellInput{3}; noise2 = cellInput{4}; aECG = mECG + fECG + noise1 + noise2; % Segment the data and keep indices so that we can also segment the QRS % peak locations [idxs,~] = buffer(1:size(mECG,1),segmentLength); mECG = single(mECG(idxs)'); fECG = single(fECG(idxs)'); aECG = single(aECG(idxs)'); % Normalize for idx = 1:size(mECG,1) mECG(idx,:) = rescale(mECG(idx,:),-1,1); mECG(idx,:) = mECG(idx,:) - median(mECG(idx,:)); fECG(idx,:) = rescale(fECG(idx,:),-1,1); fECG(idx,:) = fECG(idx,:) - median(fECG(idx,:)); aECG(idx,:) = rescale(aECG(idx,:),-1,1); aECG(idx,:) = aECG(idx,:) - median(aECG(idx,:)); end numRows = size(mECG,1); % CBT format C=1 B=numRows T=segmentLength mECG = reshape(mECG,1,numRows,[]); fECG = reshape(fECG,1,numRows,[]); aECG = reshape(aECG,1,numRows,[]); % Create cell array with individual elements --> CBT format C=1 B=1 T=segmentLength mECGCell = mat2cell(mECG,1,ones(numRows,1),segmentLength)'; fECGCell = mat2cell(fECG,1,ones(numRows,1),segmentLength)'; aECGCell = mat2cell(aECG,1,ones(numRows,1),segmentLength)'; outputCell = [aECGCell mECGCell fECGCell]; if numel(cellInput) == 6 mECG_QRSTmp = cellInput{5}; fECG_QRSTmp = cellInput{6}; segmentLimits = [idxs(1,:)' idxs(end,:)']; numSegments = size(segmentLimits,1); mECG_QRS = cell(numSegments,1); fECG_QRS = cell(numSegments,1); for idx = 1:numSegments mECG_QRS{idx} = mECG_QRSTmp(mECG_QRSTmp >= segmentLimits(idx,1) & ... mECG_QRSTmp <= segmentLimits(idx,2)); fECG_QRS{idx} = fECG_QRSTmp(fECG_QRSTmp >= segmentLimits(idx,1) & ... fECG_QRSTmp <= segmentLimits(idx,2)); end outputCell = [outputCell mECG_QRS fECG_QRS]; end end```

#### processMB

This function converts cell array inputs, containing ECG segments, to mini-batches with `CBT` format.

```function [aECGbatch,mECGbatch,fECGbatch] = processMB(aECGCell,mECGCell,fECGCell) aECGbatch = cat(2,aECGCell{:}); mECGbatch = cat(2,mECGCell{:}); fECGbatch = cat(2,fECGCell{:}); end ```

#### modelLoss

This function feeds an `aECG` input to the network and computes the gradient and resulting loss.

```function [loss,grads,state] = modelLoss(net,aECG,mECG,fECG,mECGWeight,fECGWeight) [fECGpred,mECGpred,state] = net.forward(aECG); loss = stripdims(fECGWeight*mean(abs(fECG-fECGpred),"all") + ... mECGWeight*mean(abs(mECG-mECGpred),"all")); grads = dlgradient(loss,net.Learnables); loss = double(gather(extractdata(loss))); end```

#### plotPredictedECGs

This function plots actual and predicted ECG signals for a specified measurement case, iteration, and SNR value. The function plots segments `N` to `N`+`M`.

```function plotPredictedECGs(wNet,testDS,SNRstr,caseStr,iterStr,N,M) % testDS is datastore pointing to test data % SNRstr can be "12", "09", "06", "03" % iterStr can be "I1", "I2", "I3", "I4", "I5" % caseStr can be "C0", "C1", "C3" dataIdx = N:N+M; % Get a datastore with the requested case idx = contains(string(testDS.UnderlyingDatastores{1}.Files), ... fullfile("snr"+SNRstr+"dB",iterStr+"_"+caseStr+".mat")); ds = subset(testDS,idx); data = read(ds); data = data(dataIdx,:); mECG_QRS = data(:,4); fECG_QRS = data(:,5); [aECGbatch,mECGbatch,fECGbatch] = processMB(data(:,1),data(:,2),data(:,3)); % Move the aECGbatch into a dlarray and call the predict method of the % trained network to estimate the source signals dlaECG = dlarray(aECGbatch,"CBT"); [dlpred_fECG,dlpred_mECG] = predict(wNet,dlaECG); pred_fECG = extractdata(dlpred_fECG); pred_mECG = extractdata(dlpred_mECG); % Plot the results aECG = squeeze(aECGbatch)'; mECG = squeeze(mECGbatch)'; fECG = squeeze(fECGbatch)'; pred_fECG = squeeze(pred_fECG)'; pred_mECG = squeeze(pred_mECG)'; aECG = aECG(:); mECG = mECG(:); fECG = fECG(:); pred_mECG = pred_mECG(:); pred_fECG = pred_fECG(:); mECG_QRS = cat(1,mECG_QRS{:}); fECG_QRS = cat(1,fECG_QRS{:}); mECG_QRS = mECG_QRS - ((N-1)*1024-1) - 1; fECG_QRS = fECG_QRS - ((N-1)*1024-1) - 1; titleStr = "SNR = "+SNRstr+" dB, Case = "+caseStr+", Iteration "+iterStr; figure subplot(3,2,[1 2]) plot(aECG) title("aECG mixture, "+titleStr) minECG = min(aECG); maxECG = max(aECG); axis([1 length(aECG) minECG-abs(minECG*0.35) maxECG+maxECG*0.35]) minfECG = min(fECG); maxfECG = max(fECG); minPredfECG = gather(min(pred_fECG)); maxPredfECG = gather(max(pred_fECG)); minECG = min(minfECG,minPredfECG); maxECG = max(maxfECG,maxPredfECG); subplot(3,2,3) plot(fECG) xline(fECG_QRS,":k") title("fECG target") axis([1 floor(length(fECG)/2) minECG-abs(minECG*0.35) maxECG+maxECG*0.35]) subplot(3,2,4) plot(pred_fECG) xline(fECG_QRS,":k") title("fECG predicted") axis([1 floor(length(pred_fECG)/2) minECG-abs(minECG*0.35) maxECG+maxECG*0.35]) minmECG = min(mECG); maxmECG = max(mECG); minPredmECG = gather(min(pred_mECG)); maxPredmECG = gather(max(pred_mECG)); minECG = min(minmECG,minPredmECG); maxECG = max(maxmECG,maxPredmECG); subplot(3,2,5) plot(mECG) xline(mECG_QRS,":k") title("mECG target") axis([1 length(mECG) minECG-abs(minECG*0.35) maxECG+maxECG*0.35]) subplot(3,2,6) plot(pred_mECG) xline(mECG_QRS,":k") title("mECG predicted") axis([1 length(pred_mECG) minECG-abs(minECG*0.35) maxECG+maxECG*0.35]) end```

#### computeErrorsForAllCases

This function computes the mean absolute error between actual `fECG` and `mECG` signals and predicted ones for all SNR values, measurement cases, and iterations of subject 10.

```function computeErrorsForAllCases(wNet,testDS) % testDS is datastore pointing to test data of subject 10 SNRVect = ["03" "06" "09" "12"]; caseVect = ["C0" "C1" "C3"]; for SNRidx = 1:numel(SNRVect) SNRstr = SNRVect(SNRidx); for caseIdx = 1:numel(caseVect) caseStr = caseVect(caseIdx); idx = contains(string(testDS.UnderlyingDatastores{1}.Files),fullfile("snr"+SNRstr+"dB","I1_"+caseStr+".mat")); idx = idx | contains(string(testDS.UnderlyingDatastores{1}.Files),fullfile("snr"+SNRstr+"dB","I2_"+caseStr+".mat")); idx = idx | contains(string(testDS.UnderlyingDatastores{1}.Files),fullfile("snr"+SNRstr+"dB","I3_"+caseStr+".mat")); idx = idx | contains(string(testDS.UnderlyingDatastores{1}.Files),fullfile("snr"+SNRstr+"dB","I4_"+caseStr+".mat")); idx = idx | contains(string(testDS.UnderlyingDatastores{1}.Files),fullfile("snr"+SNRstr+"dB","I5_"+caseStr+".mat")); ds = subset(testDS,idx); data = readall(ds); [aECGbatch,mECGbatch,fECGbatch] = processMB(data(:,1),data(:,2),data(:,3)); dlaECG = dlarray(aECGbatch,"CBT"); [dlpred_fECG,dlpred_mECG] = predict(wNet,dlaECG); pred_mECG = extractdata(dlpred_mECG); pred_fECG = extractdata(dlpred_fECG); errMtx(caseIdx,SNRidx) = 0.5*mean(abs(mECGbatch - pred_mECG),'all') + 0.5*mean(abs(fECGbatch - pred_fECG),'all'); errMtxFecg(caseIdx,SNRidx) = mean(abs(fECGbatch - pred_fECG),'all'); errMtxMecg(caseIdx,SNRidx) = mean(abs(mECGbatch - pred_mECG),'all'); end end figure subplot(2,1,1) plot([3 6 9 12],errMtxFecg'); title("fECG mean absolute errors") xlabel("SNR") ylabel("MAE") legend("C0","C1","C3") grid on axis tight subplot(2,1,2) plot([3 6 9 12],errMtxMecg'); title("mECG mean absolute errors") xlabel("SNR") ylabel("MAE") legend("C0","C1","C3") grid on axis tight end```

#### createWNet

This function implements a W-Net architecture and returns a layer graph.

```function lgraph = createWNet(inputSize,filterSizeLeft,numFiltersLeft,filterSizeRight,numFiltersRight) lgraph = layerGraph; inputLayer = sequenceInputLayer(1,MinLength=inputSize,Name="inputMixture"); lgraph = addLayers(lgraph,inputLayer); % Define left and right U-Net branches % Layer name conventions - left means it belongs to left U-Net % - ds means down sample, us means upsample branch, % bridge is the final row in the autoencoder % - i_j means ith row, jth layer % Add left branch U-Net lgraph = createUNet(lgraph,filterSizeLeft,numFiltersLeft,"left"); lgraph = connectLayers(lgraph,'inputMixture','conv1d_left_ds_1_1'); % Add right branch U-Net lgraph = createUNet(lgraph,filterSizeRight,numFiltersRight,"right"); lgraph = connectLayers(lgraph,'inputMixture','conv1d_right_ds_1_1'); % Connect right U-Net encoder branch to subtraction layers lgraph = connectLayers(lgraph,"avgpool1d_right_1_to_2","subtraction_2/in2"); lgraph = connectLayers(lgraph,"avgpool1d_right_2_to_3","subtraction_3/in2"); lgraph = connectLayers(lgraph,"avgpool1d_right_3_to_4","subtraction_4/in2"); lgraph = connectLayers(lgraph,"avgpool1d_right_4_to_5","subtraction_5/in2"); end```

#### createUNet

This function implements the left and right U-Net branches needed to build a W-Net architecture.

```function lgraph = createUNet(lgraph,filterSize,numFilters,branchStr) % branchStr can be "left" or "right" % % Layer name conventions - left means it belongs to left U-Net % - ds means down sample, us means upsample branch % - i_j means ith row, jth layer numFiltScale = 1 + double(branchStr == "right"); if branchStr == "left" branchStrOutput = "outputLayer_left_targetSignal"; else branchStrOutput = "outputLayer_right_secondarySignal"; end unet = [ % Row 1 encoder branch convolution1dLayer(filterSize, numFilters, Padding="same", Name="conv1d_"+branchStr+"_ds_1_1") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_ds_1_1") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_ds_1_1") convolution1dLayer(filterSize, numFilters, Padding="same", Name="conv1d_"+branchStr+"_ds_1_2") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_ds_1_2") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_ds_1_2") convolution1dLayer(filterSize, numFilters, Padding="same",Name="conv1d_"+branchStr+"_ds_1_3") batchNormalizationLayer("Name","batchnorm_"+branchStr+"_ds_1_3") leakyReluLayer(0.01,"Name","leakyrelu_"+branchStr+"_ds_1_3") averagePooling1dLayer(2, Padding="same", Stride=2, Name="avgpool1d_"+branchStr+"_1_to_2") ]; % Row 2 encoder branch if branchStr == "left" unet = [unet functionLayer(@minus,NumInputs=2,Formattable=true,Acceleratable=true,Name="subtraction_2"); tanhLayer(Name="tanh_2") ]; end unet = [unet convolution1dLayer(filterSize, numFilters*2, Padding="same", Name="conv1d_"+branchStr+"_ds_2_1") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_ds_2_1") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_ds_2_1") convolution1dLayer(filterSize, numFilters*2, Padding="same", Name="conv1d_"+branchStr+"_ds_2_2") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_ds_2_2") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_ds_2_2") averagePooling1dLayer(2, Padding="same", Stride=2, Name="avgpool1d_"+branchStr+"_2_to_3") ]; % Row 3 encoder branch if branchStr == "left" unet = [unet functionLayer(@minus,NumInputs=2,Formattable=true,Acceleratable=true,Name="subtraction_3"); tanhLayer(Name="tanh_3")]; end unet = [unet convolution1dLayer(filterSize, numFilters*4, Padding="same", Name="conv1d_"+branchStr+"_ds_3_1") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_ds_3_1") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_ds_3_1") convolution1dLayer(filterSize, numFilters*4, Padding="same", Name="conv1d_"+branchStr+"_ds_3_2") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_ds_3_2") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_ds_3_2") averagePooling1dLayer(2, Padding="same", Stride=2, Name="avgpool1d_"+branchStr+"_3_to_4")]; % Row 4 encoder branch if branchStr == "left" unet = [unet functionLayer(@minus,NumInputs=2,Formattable=true,Acceleratable=true,Name="subtraction_4"); tanhLayer(Name="tanh_4") ]; end unet = [unet convolution1dLayer(filterSize, numFilters*8, Padding="same", Name="conv1d_"+branchStr+"_ds_4_1") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_ds_4_1") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_ds_4_1") convolution1dLayer(filterSize, numFilters*8, Padding="same", Name="conv1d_"+branchStr+"_ds_4_2") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_ds_4_2") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_ds_4_2") averagePooling1dLayer(2, Padding="same", Stride=2, Name="avgpool1d_"+branchStr+"_4_to_5") ]; % Row 5 encoder branch if branchStr == "left" unet = [unet functionLayer(@minus,NumInputs=2,Formattable=true,Acceleratable=true,Name="subtraction_5"); tanhLayer(Name="tanh_5") ]; end unet = [unet convolution1dLayer(filterSize, numFilters*16, Padding="same", Name="conv1d_"+branchStr+"_ds_5_1") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_ds_5_1") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_ds_5_1") convolution1dLayer(filterSize, numFilters*16, Padding="same", Name="conv1d_"+branchStr+"_ds_5_2") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_ds_5_2") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_ds_5_2") averagePooling1dLayer(2, Padding="same", Stride=2, Name="avgpool1d_"+branchStr+"_5_to_6") % Row 6 - bridge convolution1dLayer(filterSize, numFilters*16*numFiltScale, Padding="same", Name="conv1d_"+branchStr+"_bridge_6_1") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_bridge_6_1") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_bridge_6_1") convolution1dLayer(filterSize, numFilters*16*numFiltScale, Padding="same", Name="conv1d_"+branchStr+"_bridge_6_2") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_bridge_6_2") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_bridge_6_2") transposedConv1dLayer(filterSize, numFilters*16*numFiltScale, Stride=2, Cropping="same", Name="transconv1d_"+branchStr+"_us_6_to_5") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_us_6_to_5") % Row 5 decoder branch concatenationLayer(1, 2, Name="concat_"+branchStr+"_5") convolution1dLayer(filterSize, numFilters*16, Padding="same", Name="conv1d_"+branchStr+"_us_5_1") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_us_5_1") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_us_5_1") convolution1dLayer(filterSize, numFilters*16, Padding="same", Name="conv1d_"+branchStr+"_us_5_2") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_us_5_2") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_us_5_2") transposedConv1dLayer(filterSize, numFilters*16, Stride=2, Cropping="same", Name="transconv1d_"+branchStr+"_us_5_to_4") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_us_5_to_4") % Row 4 decoder branch concatenationLayer(1, 2, Name="concat_"+branchStr+"_4") convolution1dLayer(filterSize, numFilters*8, Padding="same", Name="conv1d_"+branchStr+"_us_4_1") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_us_4_1") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_us_4_1") convolution1dLayer(filterSize, numFilters*8, Padding="same", Name="conv1d_"+branchStr+"_us_4_2") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_us_4_2") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_us_4_2") transposedConv1dLayer(filterSize, numFilters*8, Stride=2, Cropping="same", Name="transconv1d_"+branchStr+"_us_4_to_3") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_us_4_to_3") % Row 3 decoder branch concatenationLayer(1, 2, Name="concat_"+branchStr+"_3") convolution1dLayer(filterSize, numFilters*4, Padding="same", Name="conv1d_"+branchStr+"_us_3_1") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_us_3_1") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_us_3_1") convolution1dLayer(filterSize, numFilters*4, Padding="same", Name="conv1d_"+branchStr+"_us_3_2") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_us_3_2") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_us_3_2") transposedConv1dLayer(filterSize, numFilters*4, Stride=2, Cropping="same", Name="transconv1d_"+branchStr+"_us_3_to_2") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_us_3_to_2") % Row 2 decoder branch concatenationLayer(1, 2, Name="concat_"+branchStr+"_2") convolution1dLayer(filterSize, numFilters*2, Padding="same", Name="conv1d_"+branchStr+"_us_2_1") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_us_2_1") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_us_2_1") convolution1dLayer(filterSize, numFilters*2, Padding="same", Name="conv1d_"+branchStr+"_us_2_2") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_us_2_2") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_us_2_2") transposedConv1dLayer(filterSize, numFilters*2, Stride=2, Cropping="same", Name="transconv1d_"+branchStr+"_us_2_to_1") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_us_2_to_1") % Row 1 decoder branch concatenationLayer(1, 2, Name="concat_"+branchStr+"_1") convolution1dLayer(filterSize, numFilters, Padding="same", Name="conv1d_"+branchStr+"_us_1_1") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_us_1_1") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_us_1_1") convolution1dLayer(filterSize, numFilters, Padding="same", Name="conv1d_"+branchStr+"_us_1_2") batchNormalizationLayer(Name="batchnorm_"+branchStr+"_us_1_2") leakyReluLayer(0.01,Name="leakyrelu_"+branchStr+"_us_1_2") convolution1dLayer(filterSize, numFilters, Padding="same",Name="conv1d_"+branchStr+"_us_1_3") batchNormalizationLayer("Name","batchnorm_"+branchStr+"_us_1_3") leakyReluLayer(0.01,"Name","leakyrelu_"+branchStr+"_us_1_3") convolution1dLayer(filterSize, 1, Padding="same",Name=branchStrOutput) ]; lgraph = addLayers(lgraph,unet); lgraph = connectLayers(lgraph,"leakyrelu_"+branchStr+"_ds_5_2","concat_"+branchStr+"_5/in2"); lgraph = connectLayers(lgraph,"leakyrelu_"+branchStr+"_ds_4_2","concat_"+branchStr+"_4/in2"); lgraph = connectLayers(lgraph,"leakyrelu_"+branchStr+"_ds_3_2","concat_"+branchStr+"_3/in2"); lgraph = connectLayers(lgraph,"leakyrelu_"+branchStr+"_ds_2_2","concat_"+branchStr+"_2/in2"); lgraph = connectLayers(lgraph,"leakyrelu_"+branchStr+"_ds_1_3","concat_"+branchStr+"_1/in2"); end```


text:
Proving Poisson's Equation

ozone
Sorry guys I have a lot of trouble with proofs and could use your help with this one.

This is right out of Griffiths (problem 2.29) I have a solution manual, but I'd like to try to get a nudge in the right direction before I turn to it.

Poisson's eq:

$\nabla ^2 V = - \frac{ \rho }{\epsilon }$

Using an identity out of Griffiths we have $\nabla ^2\left(\frac{1}{r}\right) = -4\pi \delta ^3(r)$

finally we know that $V(r) = \frac{1}{4\pi \epsilon }\int \frac{p(r')}{r} \, dr'$

distributing the gradient into our function will yield
$V(r) = -\frac{1}{\epsilon }\int p(r')\delta ^3(r)dr'$

But I do not know what to do with the dirac delta in the final function we have arrived at.

Answers and Replies

Staff Emeritus
Science Advisor
Homework Helper
Education Advisor
Sorry guys I have a lot of trouble with proofs and could use your help with this one.

This is right out of Griffiths (problem 2.29) I have a solution manual, but I'd like to try to get a nudge in the right direction before I turn to it.

Poisson's eq:

$\nabla ^2 V = - \frac{ \rho }{\epsilon }$

Using an identity out of Griffiths we have $\nabla ^2\left(\frac{1}{r}\right) = -4\pi \delta ^3(r)$
The ##r## on the lefthand side is the distance between two points whereas the ##r## in the argument of the delta function should be the vector ##\vec{r}##.

finally we know that $V(r) = \frac{1}{4\pi \epsilon }\int \frac{p(r')}{r} \, dr'$
The potential V isn't a function of only radial distance. It's a function of all three coordinates. In my edition of Griffiths, he uses the notation V(P) to indicate it's the potential at point P. He doesn't write V(r), which implies something much different. The symbols ##p## and ##\rho## look similar, but they aren't the same symbol. Be consistent. You should also define exactly what you mean by ##r## inside the integral. How is it related to the variables of integration? Finally, ##dr'## should actually be the volume element ##d\tau##.

distributing the gradient into our function will yield
$V(r) = -\frac{1}{\epsilon }\int p(r')\delta ^3(r)dr'$

But I do not know what to do with the dirac delta in the final function we have arrived at.
The lefthand side should not be the potential V.

Electromagnetism is a fairly math-intensive course, so you should really make an effort to be precise with your notation. The sloppiness that you could get away with in earlier courses will really start to hinder you now. Clean it up so you can avoid making dumb errors and confusing yourself.


text:
# What is the period of f(t)=sin( t /2 )+ cos( (13t)/24 ) ?

Apr 25, 2016

$52 \pi$

#### Explanation:

The period of both sin kt and cos kt is $\frac{2 \pi}{k}$.

So, separately, the periods of the two terms in f(t) are $4 \pi \mathmr{and} \left(\frac{48}{13}\right) \pi$.

For the sum, the compounded period is given by $L \left(4 \pi\right) = M \left(\left(\frac{48}{13}\right) \pi\right)$, making the common value as the least integer multiple of $\pi$.

L=13 and M=1. The common value = $52 \pi$;

Check: $f \left(t + 52 \pi\right) = \sin \left(\left(\frac{1}{2}\right) \left(t + 52 \pi\right)\right) + \cos \left(\left(\frac{24}{13}\right) \left(t + 52 \pi\right)\right)$
$= \sin \left(26 \pi + \frac{t}{2}\right) + \cos \left(96 \pi + \left(\frac{24}{13}\right) t\right)$
$= \sin \left(\frac{t}{2}\right) + \cos \left(\frac{24}{13} t\right) = f \left(t\right)$..


text:
# What is the best collision detection method or algorithm for a 2d platformer?

I'm working of a 2d platformer. It's a preety simple game, but the ground is flat and it looks boring. I did the collision detection for it and it's all working when the boxes are axis aligned, but I want to add slopes so it's more interesting. I researched a lot and I found the separating axis theorem but I read that it is not good for platformers because it's not side specific. So, my question is: what is the best way to detect collision with curves and slopes? I use javascript if you want to provide some code. Thanks in advance!

For horizontal boxes AABB is great, for slopes, you need to add a twist to it.

The height of a slope that's positioned between x1 and x2 at a specific point can be calculated the following way (x is the point we want to get the height at, h1 is the height at x1 and h2 is the hight at x2):

height := (x - x1) / (x2 - x1) * (h2 - h1) + h1

Now, first you need to check if the player collides with the aabb of the slope (this is the broad phase), if it does, then get the height of it at the player's position.

Next, check if the player's height is less than the result, if it is and he's fallong, then set his y coordinate to the height of the slope and if he has one, set his vertical velocity to 0.

This doesn't work with the underside of the slopes, if you want that too, you need to do the same for that too.

• Made an edit to the above equation. Since the algorithm doesn't follow traditional PEMDAS, I added parentheses around division operation as this is required to happen first before multiplication operation. So it should look like this: height := ((x - x1) / (x2 - x1)) * (h2 - h1) + h1 – SaviorXTanren Jan 27 '17 at 3:43
• @Savior Uhm, no. Divisions and multiplications are cummilatibve, it doesn'f matter in which order I do them. Get your maths right. – Bálint Jan 27 '17 at 6:31
• No need to get touchy dude, but your statement isn't fully correct. I also think you mean "commutative" instead of "cumulative", but that's only true for multiplication, not division. If you leave out the parenthesis, some compilers will help apply them for you, but not all are guaranteed to. Additionally, putting the parenthesis on the division or the multiplication produce different results.: (noParen = 5.5; divParen = 5.5; mulParen = 4.041). Although some programming languages might be able to run PEMDAS on an equation, not all will. It's better to call out what you want done. – SaviorXTanren Jan 27 '17 at 19:16
• @SaviorXTanren What I wanted to say was that it doesn't matter in which order you apply the stuff here. (x - x1) / (x2 - x1) * (h2 - h1) always equals (x - x1) / ((x2 - x1) * (h2 - h1)) and ((x - x1) / (x2 - x1)) * (h2 - h1). I was a bit tired at the time I wrote that comment, so I didn't express myself good enough – Bálint Jan 27 '17 at 19:34
• Understandable, but as I just showed above, that's not the case. So you can't put the parenthesis in either spot and get the same result. If you type out your first equation directly into some C# code, it gives you the "right" answer. Even if you move the multiplication part to in front of the division, it still gives you the right answer. But that's also a computer doing the right thing. In the event that someone sees this and decides to breakup the equation into multiple parts, you want to make sure they are aware that depending on how they do it, they can get a wrong answer. – SaviorXTanren Jan 27 '17 at 19:37