text:
# February 2008 Archives

## Sifakas, Ringtaileds, and Sportives, Aye-Aye!

Dr. Anne Yoder and her group at the Duke Lemur Center have produced lemur family tree tree. Greg Laden has the goods.

Note: Prof. Steve Steve was not involved in the study.

## Rendering Equations with ImageMagick’s Convert

Last October I described my algorithm for rendering LaTeX equations for use on this blog: “Rendering Equations in Movable Type”. For kicks, here are instructions on how to manually run the algorithm. The only thing missing is the middle alignment part of the algorithm.

First create a file, called eqn.tex using the following skeleton. Insert your LaTeX equation code where it says to.

\nonstopmode
\documentclass{article}
\usepackage[T1]{fontenc}
\usepackage{amsmath,amsfonts,amssymb,wasysym,latexsym,marvosym,txfonts}
\usepackage[pdftex]{color}
\pagestyle{empty}
\begin{document}
\fontsize{12}{24}
\selectfont
\color{white}
\pagecolor{black}
$% % EQUATION GOES HERE %$
\end{document}

Next process the file with pdflatex eqn.tex, which will render your equation as eqn.pdf. Now we will use the convert command line tool from the ImageMagick library to turn eqn.pdf into eqn.png.

convert $$-density DENSITY eqn.pdf -trim +repage$$ \
$$+clone -fuzz 100% -fill FG -opaque black$$ \
+swap -compose copy-opacity -composite \
$$+clone -fuzz 100% -fill white -opaque BG +matte$$ \
+swap -compose over -composite eqn.png

Here DENSITY, FG, and BG are user tunable variables. DENSITY tells convert how many pixels per inch to use when rendering the pdf to a png. This determines how big of a png you have for a given rendered equation. Note that a density of 300 is print quality, 96 is windows monitor standard, and 72 is the mac standard. FG is the color of the equation text in the final image, and BG is the color of the background. Setting BG to none will make the background transparent.

Now for an explanation of the options:

$$-density DENSITY eqn.pdf -trim +repage$$ \ loads the pdf into memore, removes excess background, and saves it into the image stack position 0.

$$+clone -fuzz 100% -fill FG -opaque black$$ \ copies the image at position 0, fills it with the FG color, and saves it at position 1. We clone and fill to ensure that image 1 has the same dimensions as the trimmed image in position 0.

+swap -compose copy-opacity -composite \ uses the black-and-white values in position 0 to determine how transparent pixels are in 1, clears the stack, and saves the result to pos 0. We actually swap positions 0 and 1 to put them in the right order for the -composite operator.

$$+clone -fuzz 100% -fill white -opaque BG +matte$$ \ copies image 0 again and fills it with the BG color, discarding all transparency information.

+swap -compose over -composite eqn.png overlays image 0 onto image 1 and saves result to eqn.png.

And now you have an equation rendered as a png that can be included on any webpage.

## Estimating Local Ancestry

Last month, one of the fellow postdocs in our lab went back home to Denmark and began professoring again. So it’s nice to find out that one of his collaborators in Denmark is a blogger and has a really cool post about Estimating local ancestry. He even connects the post to some research that was done in our lab, which makes it double plus one good.

## Hiding Email Addresses From Robots with Javascript

| 1 Comment | 1 TrackBack

Spammers use automated programs to search the internet for email addresses to harvest. There are several techniques out there to hide email addresses from spam “robots”: Address Munging. Although, all these techniques can be beaten—some rather easily—you probably still want to use them because the average harvester is not designed to handle them. Sites with exposed email addresses outnumber those that protect theirs, and as long as there are easy nuts, the spammers are not going to invest the resources and time to crack the hard nuts. (Think Darwin’s Finches.)

The simplest technique is to add characters to the email address that most humans will recognize and remove or replace to form the valid address:

jason [at] theargonauts [dot] com
NOSPAMjason@theargonauts.com

This is a very effective technique; however, it requires user intervention, and the most popular techniques can still be interpreted with little effort by a robot. If user intervention is not a concern, then one nearly flawless technique is to simply display the email in an image. In my opinion this should only be used as a last resort because it has several drawbacks: no click to email ability, difficulty in copying the email address, increased download size, etc.

Transparent munging techniques can be used to allow users to click and use email address directly. The simplest technique encodes the email address using html entities:

jason@theargonauts.com

Obviously, this is not that hard to crack because everyone knows the entity-to-character translation rules ahead of time. However, we can use javascript to make things difficult to understand but easy to use. The logic here is to encrypt email addresses in webpages and use javascript to force browsers to decrypt them. The advantage here is that to all humans visiting your site, the email address behaves like a normal mailto: link. However, robots won’t understand it because they don’t implement javascript. Using a simple rot13 cypher would produces the following address:

wnfba@gurnetbanhgf.pbz

From this, javascript can be used to tell the browser that the above email address is really jason@theargonauts.com. Although generic harvesters will probably never implement the complexities of javascript, they can be programed to decode the email addresses without javascript once they learn the algorithms. However, as long as each site implements its own encode and decode algorithms, it is uneconomical for robots to specialize and hard to be comprehensive.

I use the following javascript framework to decode email addresses on my blogs:

1. Find all links with a mailto: specification and for each:
2. Match an encoded address or return
4. Split the decoded string into the real address and an optional url text
6. Set the inner html text if it was specified.

Using the tricks of the jQuery library and the assumption that email addresses are encoded as hexadecimal numbers (characters 0-9 and a-f), we have the following frame work. You will just need to supply your own decryption routine.

$(function(){$('a[href^=mailto:]').each(function(i) {
var xr = this.href.match(/^mailto:([0-9a-f]+)$/); if(!xr) return; var x = xr[1]; // encoded address var s = ''; // decoded address /* Here be decoding of 'x' into 's'. */ xr = s.split(/\|/); // optionally split real address from url text this.href = "mailto:" + xr[0]; if(xr.length > 1)$(this).html(xr[1]);
});
});

Here is an example of an encoded email address used on this page that is decoded by the script above:

  Email Me

Very cryptic and yet the page still validates—unlike other javascript email hiding implementations—while the user hardly ever notices that anything is up. Of course, if you decide to use this technique you will need a way to encode your email addresses in the first place. I wrote movabletype plugin to produce my encoded email addresses. Another option would be to setup a form on your website that would process addresses and return the proper encoding.

## Green Openlab

Bora and I met today to pack and mail copies of the Open Laboratory 2007 to most of our contributors. We sent the international ones off today, and Bora will mail the rest on Monday.

What was cool about this batch was that Lulu filled the bulk order using a different printer than the previous batches. The result was green covers! Yes, when ordered in bulk, the Openlab 2007 comes with a proper green cover. When ordered in small batches, it comes with a blue-gray cover. Unfortunately, we were not able to send green books to everyone, we had to use up our stock of blue-gray books.first. Some of the authors who ordered multiple copies will actually get both colors, just because.

I’m slowly revising the 2007 version to create a corrected, second edition, which I hope will be available during the summer. I’m hoping to fix the cover color issues at that time. We may even try a different publisher/printer since we still retain the rights to the book.

I am also working with Bora and Jennifer to extend the layout and styles we developed for this edition to next year’s.

## External Linking with XHTML and jQuery

| 1 Comment

Many of you may be familiar with the behavior on this blog and the Panda’s Thumb that links to external websites open a new window or tab. In old fashioned html, this would be accomplished by adding a “target” attribute to your link:

Scitus

However, with the advent of XHTML, content was separated from style and behavior, and the target attribute was officially retired; although, it is still recognized by standard browsers. Instead with XHTML, the proper solution is to use the “rel” tag to show how the linked page is connected to the current page:

Scitus

Perhaps you are using a browser that understands this but are probably not; therefore, we use a line of javascript to tell the browser what to do with those links:

$(function(){$('a[href][rel*=external]').each(function(i){this.target = "_blank";});});

The above code is very cryptic because it relies on shortcuts provided by the jQuery library. With prototype you can do use something equally cryptic like this:

function externalLinks() {
('a[href][rel~=external]').each( function(value, index) {value.target = "_blank";});
}

Sure, it’s more complicated than the original, but it wouldn’t be progress if unless they made things harder to do.

## WΔZ: Evolutionary Equations on the Big Screen

According to Rotten Tomatoes, the movie WΔZ starts showing today in the UK. The movie is a psychological thriller/horror movie and has been compared to Se7en. What makes this movie interesting is the fact that the screenplay was inspired by Price’s Equation:

Price’s Equation is a broader version of Fisher’s Fundamental Theorem of Natural Selection. It describes how the change in trait with phenotypes is related to the phenotypes’ fitnesses, . Note that the genetics of the trait (mutation, ploidy, etc.) is contained in the second term. See Wikipedia for more details.

Now according to the Rotten Tomatoes exclusive on WΔZ:

The script comes from City of Vice scribe Clive Bradley, who claims to have come up with the movie’s premise after flicking through a book on Darwinism. “It featured a mathematical equation—W Delta Z—formulated by American population geneticist George R. Price,” he explains. “It supposedly shows that there’s no real altruism in nature; no such thing as selflessness. Price was so upset by his findings that he ended up giving away all his possessions to the poor and, eventually homeless himself, committed suicide with a pair of nail scissors in a filthy London squat.”

The study of the evolution of altruism goes beyond the description above, and I hope moviegoers won’t be seduced by this fictional account of evolutionary theory. (I’m waiting to see what demagoguery that AiG, DI, and the Expelled frauds come up with about this movie.) Now, it is true that according to Price’s Equation, altruistic behavior that benefits a species at the cost of individual fitness is selected against. (Note that a deleterious phenotype can still exist in a population through mutation-selection balance or genetic drift.) However, if the altruism only benefits certain members of the species (e.g. relatives), then altruism can be selected for.

This is represented by Hamilton’s rule: . This describes under what conditions an altruistic allele will invade a population. is the cost of the allele to the “actor”, is the relatedness of the receiver to the actor, and is the benefit that the receiver receives by the actor being altruistic. The consequence of Hamilton’s rule is that selfish genes can still be altruistic. There is a lot of interesting literature about the evolution of altruism, including how punishment can reinforce altruism. I recommend Sean Rice’s Evolutionary Theory, Chapter 10, as a good starting point.

So if anyone in the UK goes to see this movie this weekend, please send us an overview/review.

## jQuery

I’ve migrated the javascripts on this site from the Prototype library to jQuery. Let me know if anything breaks.

## Compression Enabled

Last night I modified the server config to add compression to html, text, css, and javascript files. Let me know if you see any issues. Last time we tried this on PT, it broke IE’s refresh control.

## Jack Tice, Tragic Pioneer at Iowa State

Most people who follow this blog probably know Iowa State as the soon to be former employer of “intelligent-design” activist Guillermo Gonzalez. However, it actually has an interesting progressive history. George Washington Carver was both the first African-American student and the first African American faculty member there. Today, ESPN ran a story about Jack Tice, the namesake of the football stadium and the first African-American athlete at ISU: “ISU only I-A school to honor African-American in stadium name.”

When Euseph Messiah arrived at Iowa State, he barely noticed a weather-beaten statue outside of his team’s home stadium, much less cared whom the facility was named after.

Jack Trice’s figure can be found outside the stadium that bears his name, despite a fleeting college career that lasted only two games for the Cyclones.

That might not seem like much. But for those who have learned about Trice’s tragic story of sacrifice, it remains one of the most compelling in the history of college football.

In an age when most stadiums are named after a megabuck donor or a corporate sponsor, the Cyclones’ home facility honors Trice, who died from injuries on Oct. 8, 1923, two days after he was trampled in a football game against Minnesota.

It is a very interesting story, and well worth a read.

## Happy Darwin Day!

Today is Darwin Day. Get out and learn some biology!

(HT: PZ)

## The Market

Tiffany started her first full-time lawyer job this week. Yay, now we can pay off our credit cards and get some new stuff—50-inch plasma screen!—for the house.

I’m in the job market as well; I applied for eight faculty positions last fall. I got my first official rejection email today, which doesn’t bother me. So far, I’ve had one interview, and it’s looking like that is the only one that I’m going to have this year. I came away from the interview thinking that I’d like to work there more than any of the other places that I had applied to. I haven’t heard either way about that position, but I really hope that I get an offer because I really think it is a perfect position with a lot of opportunities for both me and Tiffany.

| 1 Comment | 1 TrackBack

The Scientist has a new article on the history of the scientific literature and controversy of the HOTHEAD mutation in Arabidopsis: “Mendel Upended? How the behavior of an Arabidopsis gene could overturn the classical laws of genetics.” Long time readers of this blog will remember that when the original paper came out I blogged my thoughts on the study and eventually turned my blog thoughts into a scientific coauthorship. (See Blog About Hothead and Get an Easy Paper for more background.)

Comai and I are mentioned in The Scientist’s article, although not by name:

In November, 2005, researchers at the University of Washington and the University of Georgia proposed that, rather than a novel genetic mechanism causing the reversion, the HTH gene product might be a metabolic enzyme, which when mutated produces toxic and mutagenic compounds that act on the DNA. In particular, they suggested that the HTH protein may be in a class of proteins that, given a deleterious mutation, revert precisely to the wild-type from mutagenic effects.

Other authors with more fantastical explanations get named and cited, but our reasonable example is tossed aside. (Of course, the most reasonable example, pollen contamination, gets a lot of discussion.) I guess one consolation prise is that if you search for “Susan Lolle” my blog comes up first.

What strikes me in all the discussions about HOTHEAD is that most of the geneticists quoted appear to have not studied population genetics beyond classical, Mendelian work. The way they frame the discussion makes me think that they’ve never studied well known phenomena like meiotic drive. (I’m currently participating in a seminar class, which is reading Burt and Trivers’s Genes in Conflict and meiotic drive gets first billing. It is very refreshing to be participating in popgen discussions again.)

Most organisms that people care about are diploid; they contain two copies of each gene. During reproduction, each copy of their genes has a 50% chance of being passed to their offspring. However, some genes show meiotic drive, which occurs when one copy gets passed to offspring a majority of the time. Consider two, unlinked genes—Aa and Bb—where A drives at 90% and B is normal. Therefore, an AaBb individual will produce four different gametes, which have the following frequencies: AB = 45%, Ab = 45%, aB = 5%, and ab = 5%. Now an AABb will have two different gametes: AB = 50% and Ab = 50%. An aaBb individual is similar.

Unaffected by other forces, a driving allele will go to fixation because it has the advantage over other alleles. However, meiotic drive often works by physically killing gametes that contain the other allele. The simple model of how this works is that the A allele is actually a pair of genes: one produces a time-delayed “poison” before gamete formation, and the other provides an “antidote” after gamete formation. Because the “poison” is in all gametes, but the antidote is only in gametes that contain the A allele, the a-gametes are weakened and thus the A-gametes get more than 50% of the fertilizations. Because nearly half of their gametes are removed, individuals with driving alleles can suffer reduced fertility, which can make it difficult for driving genes to go to fixation. Instead a polymorphic equilibrium is reached and the driving alleles persist in the population as a “chronic infection.”

If you think that meiotic drive is too exotic to be common, most organisms that have been studied extensively have shown its signs, the classic examples being the t-haplotypes in mice and sd inversions in Drosophila.

## Super Bowl XL2: Carl is Unpissed.

I was really rooting for someone to finally go 19-0, but as the last of three boys, seeing Eli win a super bowl was a nice consolation prize for me. Of course Carl is excited that the Giants won.

## Carl is Ready for Sunday.

It’s the Pats versus the Giants in the Super Bowl this weekend. Carl is psyched for another New York versus Boston championship game:

Personally, I dislike both cities, but since I had Brady on my victorious fantasy football team this year—worst to first, baby—I’m pulling for the Pats to finally make everyone stop talking about the 1972 Miami Dolphins.

What’s even better than Carl’s weekly video diary is the fact that Aqua Teen Hunger Force is now airing new episodes.—Casting John Kruk as one of three beautiful sirens was hilarious.—To add the cherry to my joy sunday, [adult swim] just released ATHF:Volume Five on DVD.

Yay!


text:
# B Kinetic Energy of two positive charges flying apart

1. May 8, 2016

### UMath1

Usually when we calculate the kinetic energy acquired by a positive charge as it moves in the direction of the electric field produced by another positive charge, we assume that the source charge remains stationary and so the loss in potential energy for the moving charge is equal to its gain in kinetic energy.

But what if both charges are allowed to move? Is the original potential energy shared between the two charges? Would you use half of the original potential energy for each of the two charges?

2. May 8, 2016

### Simon Bridge

In the lab frame:
You can work out the potential energy between the initial and final positions - this will tell you the total change in kinetic energy.
The symmetry dictates how the kinetic energy is shared.
You can also work the problem by brute force using Lagrangian mechanics.
A lot depends on what, exactly, you are trying to calculate.

3. May 8, 2016

### Khashishi

It's not hard with two charges. If they are identical charges then you can just accelerate both of them by the same amount. If not, you can convert it into a 1 body problem using the reduced mass.

4. May 8, 2016

### Simon Bridge

For N particles, positions $\vec r_i$ Newton's second law for the jth particle is $$k\sum_{i\neq j}^N \frac{q_iq_j}{|\vec r_j-\vec r_i|^3}(\vec r_j-\vec r_i) = m_j\frac{d^2}{dt^2}(\vec r_j-\vec r_i)$$
... the result is a system of N differential equations which you can solve for each $r_i(t)$ by the usual methods.

For the case of two identical particles, $q_1=q_2=q$ and $m_1=m_2=m$ this simplifies to:
$$\frac{kq^2}{|\vec r_2-\vec r_1|^3}(\vec r_2-\vec r_1) = m\frac{d^2}{dt^2}(\vec r_2-\vec r_1)\\ \frac{kq^2}{|\vec r_1-\vec r_2|^3}(\vec r_1-\vec r_2) = m\frac{d^2}{dt^2}(\vec r_1-\vec r_2)\\ \vec r_1(0)=\vec r_{01}, \vec r_2(0)=\vec r_{02}, \frac{d}{dt} \vec r_1(0) = \frac{d}{dt} \vec r_2(0) = 0$$

It is usually easier to change coordinates so that one axis (usually z) lies along $\vec r_2-\vec r_1$ with the origin half way between the particles.
This way the particles start at $\pm z_0$ given by $z_0 = \frac{1}{2}|\vec r_1-\vec r_2|$.

The equations become:
$$\frac{kq^2}{(z_2-z_1)^2}\hat k = m\frac{d^2z_2}{dt^2}\hat k\\ \frac{kq^2}{(z_2-z_1)^2}\hat k = -m\frac{d^2z_1}{dt^2}\hat k\\ z_1(0)=-z_0, z_2(0)=z_0, \dot z_1=\dot z_2 = 0$$ ... something like that.

But if the particles start out stationary at $\pm z_0$ at some time, and they are later at $\pm z_1$ so that $z_1>z_0$, then, ceteris paribus, we can work out their velocities by conservation of energy:
$$v^2 = \frac{kq^2}{m}\left(\frac{1}{z_0}-\frac{1}{z_1}\right)\\ \vec v_2=v\hat k = -\vec v_1$$

Last edited: May 8, 2016

format: plain


text:
# Magic informationally complete POVMs with permutations

Eigenstates of permutation gates are either stabilizer states (for gates in the Pauli group) or magic states, thus allowing universal quantum computation [M. Planat and Rukhsan-Ul-Haq, Preprint 1701.06443]. We show in this paper that a subset of such magic states, when acting on the generalized Pauli group, define (asymmetric) informationally complete POVMs. Such IC-POVMs, investigated in dimensions $2$ to $12$, exhibit simple finite geometries in their projector products and, for dimensions $4$ and $8$ and $9$, relate to two-qubit, three-qubit and two-qutrit contextuality.
Submitted 10 Apr 2017 to Quantum Physics [quant-ph]
Published 11 Apr 2017
Updated 20 Apr 2017
Author comments: 11 pages, 6 figures, 1 table, a few complements and corrections
Msc class: 81P50, 81P68, 81P13, 81P45, 20B05
http://arxiv.org/abs/1704.02749
http://arxiv.org/pdf/1704.02749.pdf

Blake Stacey Apr 12 2017 00:34 UTC

Clarification request: Are all the IC-POVMs in this paper minimal? That is, does the number of elements in each POVM equal the square of the dimension? If so, I am confused about the quoted value of the inner product between projectors for the equiangular IC-POVM in dimension 5.

Planat Apr 12 2017 13:58 UTC

Yes, the IC-POVMs under consideration are minimal. The IC-POVM in dimension 5 is equiangular but is also not a SIC. In particular the trace product relation of a SIC is not satisfied. For the equiangular IC-POVM in dimension 7, we have a similar result.

Blake Stacey in reply to Planat Apr 12 2017 18:20 UTC

This is why I am confused (it is probably just a reading comprehension error on my part): If the POVM is IC, it must have at least $d^2$ elements. If it is a minimal IC-POVM, it must have exactly $d^2$ elements. But if it is minimal, IC and equiangular, then the angle is fixed by the requirement that the elements sum to the identity. Suppose that the trace of $\Pi_i \Pi_j$ is $\alpha$ whenever $i \neq j$. Summing this over all $j$ yields $1 + (d^2-1)\alpha$. But the projectors $\Pi_i$ themselves must sum to $dI$, so the value of $\alpha$ is fixed to $1/(d+1)$.

Planat Apr 13 2017 07:14 UTC

The trace of pairwise product of (distinct) projectors is not constant. For example, with the state $(0,1,-1,-1,1)$, one gets an equiangular IC-POVM in which the trace is trivalued: it is either $1/16$, or $(7 \pm 3\sqrt{5})/32$. For the state (0,1,i,-i,-1), there are five values of the trace.
We should explicit this observation in the next version of the paper.

Robin Blume-Kohout in reply to Planat Apr 13 2017 15:40 UTC (1 points)

This appears to be an odd and nonstandard definition of "equiangular", unless I'm missing something? Most references I'm aware of, including [Wikipedia][1] and [Renes et al 2004][2] agree that "a set of lines is called equiangular if every pair of lines makes the same angle". For unit vectors (rays), that means they have the same inner product (or absolute value thereof).

In the case you're describing, it sounds like the angles are not all equal, but chosen from a small discrete set, as (e.g.) in [this recent paper][3], or some of the older references therein. MUBs and stabilizer states are examples of k-angle tight frames, albeit not minimal ones.

[1]: https://en.wikipedia.org/wiki/Equiangular_lines
[2]: https://arxiv.org/abs/quant-ph/0310075
[3]: https://arxiv.org/abs/1605.09429

Planat Apr 13 2017 16:49 UTC

To define the complex angle, we used the (cyclotomic) field norm to the power one over the degree of the field, as stated in the introduction. It recovers the particular case of angles for SICs. In this sense "equiangular" means that all pairs of distinct lines make the same angle.

Marcus Appleby Apr 13 2017 18:54 UTC (2 points)

I am a little confused by this. As I use the term, lines are equiangular if and only if the "trace of pairwise product of (distinct) projectors is constant". You seem to be using the word in a different sense. It might be helpful if you were to explain exactly what is that sense.

Marcus Appleby Apr 13 2017 19:16 UTC (2 points)

I worded that badly, since you clearly have explained the sense in which you are using the word. I am wondering, however, how your definition relates to the usual one. Is it a generalization? Or just plain different? For instance, would a SIC be equiangular relative to your definition (using SIC field norm instead of cyclotomic field norm)?

Steve Flammia Apr 13 2017 19:16 UTC (2 points)

Just to clarify Michel's earlier remark, the field norm for the cyclotomics defines the norm in which these vectors are equiangular, and then they will generally **not** be equiangular in the standard norm based on the Hilbert-Schmidt inner product. In the example that he quotes,
$$\|(7\pm 3 \sqrt{5})/32\|_{\mathbb{Q}(\sqrt{5})} = \left[\frac{(7\pm 3 \sqrt{5})}{32} \frac{(7\mp 3 \sqrt{5})}{32}\right]^{1/\deg(\mathbb{Q}[\sqrt{5})]} = \frac{1}{16}.$$
It might be helpful in v2 of the paper if these vectors are called "generalized equiangular" or "equiangular with respect to the field norm", as this will help avoid confusion.

Blake Stacey in reply to Steve Flammia Apr 14 2017 00:35 UTC

I agree with Steve Flammia's comment. The field norm is a nice generalization of the standard norm. (I haven't yet thought about whether there might be a physics motivation for it, rather than a purely mathematical one, but that's not important right now.) To avoid confusion, some phrase like "equiangular with respect to the field norm" or "field-norm equiangular" should be used.

Robin Blume-Kohout Apr 14 2017 03:03 UTC

Okay, I see the resolution to my confusion now (and admit that I was confused). Thanks to Michel, Marcus, Blake, and Steve!

Since I don't know the first thing about cyclotomic field norms... can anybody explain the utility of this norm, for this problem? I mean, just to be extreme, I could define a trivial norm that is 1 for all vectors except $\vec{0}$, and then all rank-1 POVMs would be equiangular. I'm not by any means saying that this is what's done here! My point is that there exist norms for which equi-angularity is less interesting than others. The Hilbert-Schmidt norm is very relevant for quantum states in Hilbert spaces, because it's what appears in Born's rule. What can I do with this field norm that makes it interesting and relevant?

P.S. @Marcus, if I'm understanding this correctly, then whenever two pairs have equal Hilbert-Schmidt norm, they will have equal field norm (but different H-S norms can correspond to equal field norms). So SICs should still be equiangular in field norm. Unless I'm misunderstanding again!

Planat Apr 14 2017 08:11 UTC (1 points)

First of all, thanks to all for helping to clarify some hidden points of our paper.
As you can see, the field norm generalizes the standard Hilbert-Schmidt norm.
It works for SIC [e.g. d=2, d=3 (the Hesse) and d=8 (the Hoggar)].

The first non-trivial case is with d=4 when one needs to extend the rational field
by a 12th root of unity, i.e. n=GCD(d,r)=GCD(4,3)=12, that is r=3 for defining the appropriate
fiducial state and d=4 to allow the action of the two-qubit Pauli group on it.
Then one needs the field norm in the so defined cyclotomic extension to normalize the vectors of the
resulting IC-POVM. This IC is dichotomic in angles and traces of paiwise products.

Incidently, such a 4-dimensional IC is related to the Mermin square through the traces of triple products.


text:
# Rewrite function as Taylor series equal to natural logarithm of some value

How do I rewrite $$\sum_{n=1}^\infty\frac{\frac{1}{3}^n}{n}$$ as a Taylor series of a function for some x in order to find the sum? I also know that the sum is the natural log of three-halves.

I have tried squaring the Taylor series of the natural logarithm of 1+x in order to get rid of the alternating part of the series, but this does not equal the correct sum.

That sum is equal to $$\sum_{n=1}^\infty\frac{x^n}n$$, when $$x=\frac13$$. But$$x\in(-1,1)\implies\sum_{n=1}^\infty\frac{x^n}n=-\log(1-x)$$and so the sum of your series is$$-\log\left(\frac23\right)=\log\left(\frac32\right).$$


text:
# Solving $\lim\limits_{x \to 0^+} \frac{\ln[\cos(x)]}{x}$

I had a test today which involved solving the following limit:

$$\lim\limits_{x \to 0^+} \frac{\ln[\cos(x)]}{x}$$

I didn't figure out how to solve. After the test, I asked a couple of classmates and they told me that it was supposed to be solved by first transporting the multiplication by $\frac{1}{x}$ inside the logarithm as an exponent and then replacing $\cos(x)$ with $\sqrt{1-\sin^2(x)}$, giving the following expression:

$$\lim\limits_{x \to 0^+} \ln[(1-\sin^2(x))^\frac{1}{2x}]$$

However, I don't know how to proceed from here. It's not like it matters a lot at this point, but I'd still like to know how I'm supposed to solve this. I'm not sure if it's applicable here, but we haven't learned L'Hôpital's rule.

-

Note that $\ln(\cos(0))=0$.

So we can write our limit as $$\lim_{x\to 0^+} \frac{\ln(\cos(x))-\ln(\cos(0))}{(x-0)}.$$

Note that the above expression is almost the usual expression for the derivative of $\ln(\cos(x))$ at $x=0$. (If necessary, go back and look up the definition of the derivative of $f(x)$ at $x=a$). The only difference is the use of a one-sided limit. But what about if the two-sided limit existed?

If we are very lucky and the derivative of $\ln(\cos(x))$ at $0$ exists, the value of that derivative at $x=0$ will be our answer.

So differentiate $\ln(\cos(x))$ in the usual way. Everything works out nicely, the derivative is $0$.

Added: I expect there is no issue in finding the derivative, but here are the details. Using the Chain Rule, we get $$-(\sin(x))\frac{1}{\cos(x)}.$$ At $x=0$ this is $0$.

So the $0^+$ turns out to be unnecessary, plain old $0$ will do. The manipulations suggested by classmates are not needed, everything follows from the definition of derivative, if we know a couple of differentiation rules.

Comment: This is not really how I would do it, if I needed to know the answer. The "natural" approach is to use the power series expansions of $\cos(x)$ and $\ln(1+u)$. But since you mentioned that you had not yet done L'Hospital's Rule, I assumed that you would not yet have been exposed to power series.

But the power series approach is very much worth knowing. The power series expansion of $\cos x$ is $$1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!} +\cdots.$$

So very informally, if $x$ is near $0$, $\cos x$ is about $1-x^2/2$.

The power series expansion of $\ln(1+u)$ is $$u-\frac{u^2}{2}+\frac{u^3}{3} -\frac{u^4}{4}+\cdots.$$ (This expansion is only valid when $-1 \lt u \le 1$.)

So when $x$ is near $0$, $\ln(\cos(x))$ is about $-x^2/2$. Divide by $x$. We get $-x/2$, which approaches $0$ as $x$ approaches $0$.

-
Forgive me if I'm missing something, but what does the derivative have to do with this? –  Javier Badia Jul 16 '11 at 1:50
@Javier: what's the definition of the derivative? –  lentic catachresis Jul 16 '11 at 1:53
@Bruno: So what this is doing is sort of taking the expression and trying to transform it into something that will fit the definition of the derivative and then work from there? –  Javier Badia Jul 16 '11 at 1:56
@Javier: Exactly! You try to see your limit as the derivative of a function on a given point. Once you got it that way, you calculate the derivative of the recognized function using standard methods, evaluate it on the point, and voilà. –  lentic catachresis Jul 16 '11 at 1:59
@Javier Badia: There is seldom a "best" answer. As I pointed out, in real life I would use series, to get a handle on the error term. However, exams are not real life. If this was on a calculus exam or problem set, after derivatives but before series, the derivative approach is undoubtedly what you were expected to use. If the problem was medium early in an analysis course, something like what Aryabhata did would be the expected kind of approach. –  André Nicolas Jul 16 '11 at 14:38

We can also do this without the knowledge of derivatives (which I presume you haven't yet gotten to).

We know that $\displaystyle \lim_{x \to 0} (1-x)^{1/x} = e^{-1}$. Thus for $x$ sufficiently close to $\displaystyle 0$

we have that $\displaystyle (1-x)^{1/x} \gt e^{-2}$ and thus taking logs, we get

$$\frac{\ln (1-x)}{x} \ge -2$$ and so

$$\ln(1-x) \ge -2x$$ for $\displaystyle x$ sufficiently close to $\displaystyle 0$.

Thus

$\displaystyle \ln (\cos x) = \frac{1}{2} \ln (1 - \sin^2 x) \ge -\sin^2 x$

Thus for $\displaystyle x$ sufficiently close to $\displaystyle 0$ we have that

$$0 \ge \frac{\ln (\cos x)}{x} \ge \frac{-\sin^2 x}{x}$$

Since we know that $\displaystyle \lim_{x \to 0} \frac{\sin x}{x} = 1$, by the sandwich theorem we have that the limit you seek is $\displaystyle 0$, as $\displaystyle \frac{\sin^2 x}{x} \to 0$ as $\displaystyle x \to 0$.

-
This seems to make sense, but why did you write $\frac{\ln(1-\cos x)}{x}$ instead of $\frac{\ln(\cos x)}{x}$ in the last step, and how do you know that it's less than or equal to 0? –  Javier Badia Jul 16 '11 at 2:39
@Javier: That is a typo :) And $\ln t \le 0$ for $0 \lt t \lt 1$, –  Aryabhata Jul 16 '11 at 2:41
Of course, should have realized that. Thanks. –  Javier Badia Jul 16 '11 at 2:45

HINT $\$ Many limits can be simply calculated by recognizing them as instances of first derivatives and then calculating the derivatively rotely using known derivative rules. You can find a handful of examples in my prior posts starting here.

-
this is a good point to know (+1). –  Chris's sis Jun 20 '12 at 18:06

Let \begin{align*} \ell = \lim_{x \to 0+} (\cos{x})^{\frac{1}{x}} &= \lim_{x \to 0+} (1+ (\cos{x}-1))^{\frac{1}{x}} \\&= \lim_{x \to 0+}(1+(\cos{x}-1))^{\frac{1}{\cos{x}-1} \times \frac{\cos{x}-1}{x}} \\ &= e^{\scriptsize{\displaystyle\lim_{x \to 0+} \frac{\cos{x}-1}{x}}} \\ &= e^{0} =1 \qquad\qquad\Bigl[ \because \small\lim_{x \to 0} \frac{\cos{x}-1}{x} = \lim_{x \to 0}\frac{-2\:\sin^{2}\frac{x}{2}}{x} = -\:\lim_{x \to 0} \frac{\sin\frac{x}{2}}{\frac{x}{2}} \cdot \sin\frac{x}{2} \Bigr] \end{align*}

Hence $\displaystyle \lim_{x \to 0+} \ln(\ell) = \ln(1) =0$.

-

$$\lim_{x \to 0} \dfrac{\ln(\cos(x))}{x} = \lim_{x \to 0} \dfrac{\ln(1-2\sin^2(x/2))}{x} = \lim_{x \to 0} \dfrac{\ln(1-2\sin^2(x/2))}{2 \sin^2(x/2)} \dfrac{\sin^2(x/2)}{x/2} = \lim_{x \to 0} \dfrac{\ln(1-2\sin^2(x/2))}{2 \sin^2(x/2)} \lim_{x \to 0} \dfrac{\sin(x/2)}{x/2} \lim_{x \to 0} \sin(x/2) = -1 \times 1 \times 0 = 0$$

-


text:
# Definition

Suppose $$X \sim \chi^2_m$$ is a Chi-squared random variate on $$m$$ degrees of freedom. Then the distribution of $Y = \sqrt{\frac{X}{m}}$ is the Kay distribution on $$m$$ degrees of freedom, written as $Y \sim K_m.$ Its density is $f(y) = \left\{\begin{array}{lcl} \frac{m^{\frac{m}{2}} y ^{m-1} e^{-\frac{1}{2} m y^2}} {2^{\frac{m}{2}-1} \Gamma(\frac{m}{2})} &~~~& \text{for} ~~ 0 \le y < \infty \\ &&\\ 0 && \text{otherwise}. \end{array} \right.$

The $$K_m$$ density has some very attractive features over the $$\chi^2_m$$ density:

• $$K_m$$ has a much more symmetric density than had the $$\chi^2_m$$, for any $$m$$;
• like a $$\chi^2_m$$ density a $$K_m$$ density, as $$m \rightarrow \infty$$ $$K_m$$ becomes more symmetric and nearly normally distributed but does both faster than a $$\chi^2_m$$;
• as $$m \uparrow$$, the $$K_m$$ density concentrates around the value $$y=1$$, rather than heading off to $$\infty$$ like the $$\chi^2_m$$;

These values were calculated using the dkay(...) density function. For example, dkay(1.0, df=10) = 1.7546737.

## Normal theory relations

Perhaps the most obvious relation between a normal random variate and a $$K_m$$ is that if $$Z \sim N(0,1)$$, then $$|Z|\sim K_1$$, the half-normal.

More important in applications is that distribution of the estimator of the sample standard deviation is proportional to a $$K_m$$. To be precise, if $$Y_1, \ldots, Y_n$$ are independent and identically distributed as $$N(\mu, \sigma^2)$$ random variates, with realizations $$y_1, \ldots, y_n$$ and the usual estimates $$\widehat{\mu} = \sum y_i /n$$ and $$\widehat{\sigma} = \sqrt{\sum (y_i - \widehat{\mu})^2/(n-1)}$$, then the corresponding estimators $$\widetilde{\mu}$$ and $$\widetilde{\sigma}$$ are distributed as $\widetilde{\mu} \sim N(\mu, \frac{\sigma^2}{n}) ~~~~~\text{and} ~~~~~ \frac{\widetilde{\sigma}}{\sigma} \sim K_{n-1}.$ The latter shows that $$K_m$$ is used for inference (e.g. tests and confidence intervals) about $$\sigma$$.

This is handy because the $$K_m$$ quantiles vary much less than do those of $$\chi^2_m$$. For example, condider the following table of the cumulative distribution.

df p=0.05 p=0.5 p=0.95
1 0.0627068 0.6744898 1.959964
2 0.2264802 0.8325546 1.730818
3 0.3424648 0.8880642 1.613973
4 0.4215220 0.9160641 1.540108
5 0.4786390 0.9328944 1.487985
6 0.5220764 0.9441152 1.448654
7 0.5564364 0.9521263 1.417601
8 0.5844481 0.9581311 1.392269
9 0.6078297 0.9627987 1.371090
10 0.6277180 0.9665308 1.353035
15 0.6957463 0.9777136 1.290886
20 0.7365735 0.9832962 1.253205
25 0.7644974 0.9866425 1.227232
30 0.7851255 0.9888719 1.207932
35 0.8011601 0.9904636 1.192858
40 0.8140839 0.9916570 1.180662

Unlike the $$\chi^2_m$$ distribution, the quantiles in this table stabilize, allowing $$1 \pm 0.20$$ being not a bad rule of thumb for a $$90\%$$ probability of the ratio $$\widetilde{\sigma}/\sigma$$.

These values were calculated using the qkay(...) quantile function. For example, qkay(0.05, df=5) = 0.478639. These would be used to construct interval estimates for $$\sigma$$.

To get observed significance levels, the cumulative distribution function pkay(...) would be used. For example, SL = 1- pkay(1.4, df=10) = 1 - 0.9667287 = 0.0332713.

### The Student t distribution

For the standard normal theory, the Student $$t_m$$ distribution can be defined as follows. If $$Z \sim N(0,1)$$ and $$Y \sim K_m$$ is distributed independently of $$Z$$, then the ratio $T=\frac{Z}{Y} = \frac{N(0,1)}{K_m} = t_m$ which is fairly easy to remember.

For the estimators from the above model $\frac{\widetilde{\mu} - \mu} {\widetilde{\sigma} / \sqrt{n}} = \frac{ \frac{\widetilde{\mu}-\mu} {\sigma/\sqrt{n}} } {\frac{\widetilde{\sigma}} {\sigma} } = \frac{N(0,1)}{K_{n-1}} = t_{n-1}$ is used to construct interval estimates and tests for the value of the parameter $$\mu$$.

# The functions

As with every other distribution in R four functions are provided for the $$K_m$$ distribution. These are

• dkay(x, df=m, ...) which evalutes the density of $$K_m$$ at $$x$$,
• pkay(x, df=m, ...) which evalutes the distribution of $$K_m$$ at $$x$$,
• qkay(p, df=m, ...) which evalutes the quantile of $$K_m$$ at the proportion $$p$$,
• rkay(n, df=m, ...) which generates $$n$$ pseudo-random realizations from $$K_m$$.

The parameters in the ellipsis include a non-centrality parameter. All functions rely on the corresponding $$\chi^2_m$$ functions in base R.

We briefly illustrate each below.

## The density dkay(x, df, ...)

x <- seq(0,2,0.01)
plot(x, dkay(x, df=10), type="l", col="steelblue",
main="Density", xlab="x", ylab="f(x)")
abline(v=1.0, lty=2, col="grey")

## The cumulative distribution function pkay(x, df, ...)

x <- seq(0,2,0.01)
plot(x, pkay(x, df=10), type="l", col="steelblue",
main="Distribution", xlab="x", ylab="F(x)")
abline(v=1.0, lty=2, col="grey")

## The quantile function qkay(p, df, ...)

x <- seq(0,2,0.01)
p <- pkay(x, df=10)
plot(p, qkay(p, df=10), type="l", col="steelblue",
main="Quantile function", xlab="p", ylab="Q(p)")
abline(h=1.0, lty=2, col="grey")

## Pseudo-random realizations rkay(n, df, ...)

x <- rkay(1000, df=10)
hist(x, col="steelblue",
main="Pseudo-random numbers", xlab="x")
abline(v=1.0, lty=2, col="grey")


text:
# How to calculate volume

## How to calculate volume

Height × width × depth = volume

If the height, width and depth are measured in cm, the answer will be cm³.

If the height, width and depth are measured in m, the answer will be m³


text:
# Forward simulation the 1D Magnetotelluric (MT) problem¶

In the previous notebook, we walked through how to discretize and solve the 1D Magnetotelluric (MT) problem using a finite difference approach. In this notebook, we will use the numerical simulation to simulate MT data and explore concepts including

• Mesh Design: we will compare our numerical result to an analytic and see what happens when cell sizes are too large, or the domain doesn't extend far enough
• Non-uniqueness: prior to inverting geophysical data, it is important to set expectations on what we hope to recover from those data. Forward modelling is a powerful tool for getting a handle on this.
In [1]:
import numpy as np
import scipy.sparse as sp
import matplotlib.pyplot as plt
from SimPEG import Mesh, Utils, Solver
from scipy.constants import mu_0, epsilon_0
import matplotlib
matplotlib.rcParams["font.size"] = 13
%matplotlib inline

In the previous notebook, we generated the function simulateMT in the file MTforward.py, we will import that and use it in this notebook.

In [2]:
from MTforward import simulateMT

## Set up Model Parameters¶

To define an MT simulation, we will set up a conductivity model a frequency range over which we wish to simulate. We will use:

• a background resistivity of $\rho = 100 \Omega m$ ($10^{-2} Sm$)
• 25 frequencies between $10^{-3}$ Hz and 100 Hz
In [3]:
rho_half = 100.  # Resistivity of the halfspace in Ohm-m
sigma_half = 1./rho_half  # Conductivity is the inverse of conductivity

In [4]:
frequency = np.logspace(-3, 2, 25)  # frequencies at which to simulate the MT problem

## Skin depth and Mesh Design¶

When setting up a mesh, we want to make sure our fine cells are fine enough to capture the behaviour at the highest frequencies and that the domain extends far enough so that the fields have sufficiently decayed by the time they reach the boundary. To gauge this, we will examine the skin depth at the highest and lowest frequencies. Skin depth ($\delta$) is the distance at which the amplitude of an EM wave propagating through a homogeneous medium will have decayed by a factor of $1/e$

Skin Depth $$\delta = \frac{500}{\sqrt{\sigma f}}$$

In [5]:
def skin_depth(sigma, f):
"""
Depth at which the fields propagating through a homogeneous
medium have decayed by a factor of 1/e for a given
frequency, f and conductivity, sigma
"""
return 500./np.sqrt(sigma * f)

In [6]:
skin_depth_min = skin_depth(sigma_half, frequency.max())
skin_depth_max = skin_depth(sigma_half, frequency.min())

print("The minimum skin depth is {:1.2f}m".format(skin_depth_min))
print("The maximum skin depth is {:1.2e}m".format(skin_depth_max))

The minimum skin depth is 500.00m
The maximum skin depth is 1.58e+05m

### Mesh parameters¶

We start by choosing a mesh with parameters as follows:

• smallest cell size: $$\Delta z_{\rm core} = \delta_{\rm min} / 4$$

• domain extent: $$z_{\rm max} = 2 \delta_{\rm max}$$

For this example, we will be exploring model variations within the top 5km, so we will extend the core region (the region of the having uniform cells with width $\Delta z_{\rm core}$) of the mesh to a depth of 5km.

In [7]:
cs = skin_depth_min / 4.
core_extent = 5000.
domain_extent = 2 * skin_depth_max

print("The smallest cell size is {:1.2f}m".format(cs))
print("The core region of the mesh extends {:1.2e}m".format(core_extent))
print("The mesh should extend at least {:1.2e}m".format(domain_extent))

The smallest cell size is 125.00m
The core region of the mesh extends 5.00e+03m
The mesh should extend at least 3.16e+05m

We will use a tensor mesh, which means we can use non-uniform cells for the padding, that is, expanding the width of the cells with depth. We can get away with this because EM fields and fluxes are diffusive and high frequencies will be attenuated as they move through the conductive earth (eg. in seismic, you wouldn't want to do this as the seismic response is dominated by wave propagaion).

We expand by a factor of 1.3 until we are beyond the desired domain extent. Here, we write a small while loop to figure out how many padding cells we should use.

In [8]:
npad = 1  # start with 1 cell
padding_fact = 1.3  # the amount by which we will expand each cell of the padding

"""
given a number of padding cells, this computes how far
"""

# keep adding padding until we are beyond the desired extent
while padding_z < domain_extent:

print(
"{:1.0f} padding cells extends {:1.2e}m > {:1.2e}m "
"(2 skin depths)".format(
)
)

25 padding cells extends 3.82e+05m > 3.16e+05m (2 skin depths)

#### Construct the Tensor mesh¶

Now that we have defined all of the mesh parameters, we use the Mesh Class in SimPEG to construct a mesh. This will define all of the geometries, provide the necessary differential operators, and some handy functions like plotting.

In [9]:
ncz = np.ceil(core_extent / cs)  # number of cells in the core domain
hz = [(cs, npad, -1.3), (cs, ncz)]  # define how to construct the cell widths
mesh = Mesh.TensorMesh([hz], x0='N')  # construct a 1D Tensor Mesh

print(
"There are {:1.0f} cells in the mesh. The mest extends {:1.2e}m".format(
ncz, mesh.hx.sum()
)
)

There are 40 cells in the mesh. The mest extends 3.87e+05m

In [10]:
# plot the mesh
fig, ax = plt.subplots(1,1, figsize=(8, 3))
mesh.plotGrid(centers=True, faces=True, ax=ax)
ax.legend(["centers", "faces"])
ax.grid(which="both", linewidth=0.5)
ax.invert_xaxis()  # so that the surface is on our left hand side
ax.set_xlabel('z (m)')

Out[10]:
Text(0.5,0,'z (m)')

## Set up a model¶

In the previous notebook, we showed that for a half-space, we expect the apparent resistivity computed from impedance data to be the same as the half-space resistivity, and the phase to by $45^\circ$ across the entire frequency range.

What happens when we include a conductive target in our model?

Here we will use the following model parameters

• Resistivity of the background: $\rho_{\text{halfspace}} = 100 \Omega m$ (which is the same as $\sigma_{\text{halfspace}} = 10^{-2} S/m$)
• Resistivity of the target: $\rho_{\text{target}} = 10 \Omega m$ (which is the same as $\sigma_{\text{halfspace}} = 10^{-1} S/m$)
• Thickness of the target layer: 1000m
• Depth of the target layer: 2000m
In [11]:
rho_target = 10.  # resistivity in Ohm-m
depth = 2000.  # depth to the top of the target in m
thickness = 1000.  # thickness of the target in m

In [12]:
# put the model on the mesh
sigma = 1./rho_half * np.ones(mesh.nC)

# find the indices of the layer
layer_inds = (
(mesh.vectorCCx<=-depth) &
(mesh.vectorCCx>-(depth+thickness))
)
sigma[layer_inds] = 1./rho_target

In [13]:
# plot the model
fig, ax = plt.subplots(1, 1, figsize=(8, 3))

# trickery to plot from node to node rather than at cell centers
z = np.repeat(mesh.vectorNx[1:-1], 2, axis=0)
z = np.r_[mesh.vectorNx[0], z, mesh.vectorNx[-1]]
sigma_plt = np.repeat(sigma, 2, axis=0)

ax.semilogy(z, sigma_plt,"C0", lw=2)
ax.grid(which="both", linewidth=0.5)
ax.set_xlim([-5000., 0.])
ax.set_ylim([5e-3, 1])
ax.invert_xaxis() # plot the surface on the left

ax.set_xlabel("Elevation (m)", fontsize=14)
ax.set_ylabel("Conductivity (S/m)", fontsize=14)

Out[13]:
Text(0,0.5,'Conductivity (S/m)')

## Compare with an analytic¶

Prior to drawing conclusions from our simulations, we want to first make sure they are correct!

So here, we compare the numerical results computed with our method simulateMT with an analytic for the MT response over a layered earth. A more complete description, including the derivation, is available on EM GeoSci

In [14]:
from SimPEG.EM.Analytics import MT_LayeredEarth

In [15]:
# the analytic takes the frequencies, layer thicknesses and layer conductivities
sigma_layers = np.r_[
1./rho_half,
1./rho_target,
1./rho_half]
h = np.r_[depth, thickness]

app_res_ana, app_phase_ana = MT_LayeredEarth(
frequency, h, sigma_layers, 'Res-Phase'
)

In [16]:
# # Uncomment to read the documentation on MT_LayeredEarth
# MT_LayeredEarth??

In [17]:
# numerically compute the response
app_res, app_phase = simulateMT(mesh, sigma, frequency)

In [18]:
def plot_with_analytic(frequency, app_res, app_phase, app_res_ana, app_phase_ana):
# Plot and compare the results
fig, ax = plt.subplots(2, 1, figsize=(8, 3*2))

# apparent resistivity
ax[0].loglog(frequency, app_res, label='Numeric')
ax[0].loglog(frequency, app_res_ana, 'k.', label='Analytic')
ax[0].set_ylabel("$\\rho_a \ (\Omega m)$", fontsize=14)
ax[0].set_ylim([2e1, 3e2])

# phase
ax[1].semilogx(frequency, app_phase, label='Numeric')
ax[1].semilogx(frequency, app_phase_ana, 'k.', label='Analytic')
ax[1].set_ylabel("$\phi \ (^{\circ})$", fontsize=13)
ax[1].set_ylim([0., 90.])

for a in ax:
a.grid(True, which='both', linewidth=0.3)
a.set_xlim(frequency.max(), frequency.min())
a.set_xlabel("Frequency (Hz)", fontsize=14)
a.legend(fontsize=11)

plt.tight_layout()
plt.show()

In [19]:
plot_with_analytic(frequency, app_res, app_phase, app_res_ana, app_phase_ana)

Across the entire frequency band, the apparent resistivity and phase agree with the analytic. What happens if

1. we use a coarse cell size?
2. we reduce how far the padding extends?

To answer those questions, we have a small widget.

In [20]:
import ipywidgets

In [21]:
def mesh_design_app(domain_extent, cs):
ncz = np.ceil(core_extent / cs)  # number of cells in the core domain

# keep adding padding until we are beyond the desired extent
while padding_z < domain_extent*1e3:

hz = [(cs, npad, -1.3), (cs, ncz)]  # define how to construct the cell widths
mesh = Mesh.TensorMesh([hz], x0='N')  # construct a 1D Tensor Mesh

# put the model on the mesh
sigma = 1./rho_half * np.ones(mesh.nC)

# find the indices of the layer
layer_inds = (
(mesh.vectorCCx<=-depth) &
(mesh.vectorCCx>-(depth+thickness))
)
sigma[layer_inds] = 1./rho_target

app_res, app_phase = simulateMT(mesh, sigma, frequency)

plot_with_analytic(frequency, app_res, app_phase, app_res_ana, app_phase_ana)

In [22]:
print(
"It is helpful to remember that \n"
"the minimum skin depth is {:1.0f}m \n"
"and the maximum skin depth is {:1.0f}km".format(
skin_depth_min, skin_depth_max*1e-3
)
)

It is helpful to remember that
the minimum skin depth is 500m
and the maximum skin depth is 158km

• What happens if you make the domain smaller? Which frequencies are impacted first?
• What happens if you use a larger cell size? (Pay attention to the phase at high frequencies)
• notice that there will be two regions impacted: high frequencies, and the frequencies impacted by the conductive layer
In [23]:
skin_depth_min

Out[23]:
500.0
In [24]:
ipywidgets.interact(
mesh_design_app,
domain_extent=ipywidgets.FloatSlider(min=0, max=400, step=25, value=300, description="domain extent (km)"),
cs=ipywidgets.FloatSlider(min=25, max=1000, step=25, value=125, description="min cell size (m)")
)

Out[24]:
<function __main__.mesh_design_app>

### Considerations for the inversion¶

Prior to considering an inversion, first making sure that a suitable mesh is designed is an important step - otherwise, you can introduce artifacts and structures simply because the mesh doesn't capture the physics! Comparing against analytics that model a simplified, but similar, scenario can be a powerful tool for mesh design.

## Non-uniqueness¶

The second concept we will explore with forward modelling is non-uniqueness: when models are different, but the data are similar.

In the Magnetotelluric problem, a classic example that demonstrates non-uniqueness of the data is the non-uniquness associated with the conductivity-thickness product of a thin layer. If we start with a layer that has a conductivity of $\sigma$, halve its thickness and double its conductivity, the resulting data will be similar... (but don't just take our word for it!)

In [25]:
# compute models with equivalent conductivity-thickness product
conductivity_thickness = thickness * 1./rho_target
layer_thicknesses = (np.arange(5)+1.)*300.
print("we consider layers having thicknesses (m) {}".format(layer_thicknesses))

we consider layers having thicknesses (m) [ 300.  600.  900. 1200. 1500.]

In [26]:
sig_targets = conductivity_thickness / layer_thicknesses
print(
"To preserve the conductivity thickness product,\n"
"these layer have conductivities (S/m) of {}".format(sig_targets)
)
print("This is the same as resistivities (Ohm m) of {}".format(1./sig_targets))

To preserve the conductivity thickness product,
these layer have conductivities (S/m) of [0.33333333 0.16666667 0.11111111 0.08333333 0.06666667]
This is the same as resistivities (Ohm m) of [ 3.  6.  9. 12. 15.]

In [27]:
# construct conductivity models on the mesh
sigma_list = []

# Compute the midpoint and make sure the layers we will compare with are all centered at the same depth
depth_mid = depth + thickness / 2.

for i in range(len(layer_thicknesses)):
layer_inds = (
(mesh.vectorCCx<-depth_mid+layer_thicknesses[i]*0.5) &
(mesh.vectorCCx>-depth_mid-layer_thicknesses[i]*0.5)
)
sigma_temp = np.ones(mesh.nC) * 1./rho_half
sigma_temp[layer_inds] = sig_targets[i]
sigma_list.append(sigma_temp)

In [28]:
# compute the apparent resistivity and phase for each model
app_res_list = []
phase_list = []
for i in range(len(layer_thicknesses)):
app_res, phase = simulateMT(mesh, sigma_list[i], frequency)
app_res_list.append(app_res)
phase_list.append(phase)

In [29]:
def plot_apparent_resistivities_and_phases(
frequency, sigma_list, layer_thicknesses, app_res_list, phase_list, app_res_uncert=0., phase_uncert=0.
):
"""
Plot the apparent resistivity and phase for multiple models

frequency: numpy array of frequencies
sigma_list: list of conductivity models to plot
layer_thicknesses: thickness of the ~equivalent layers
app_res_list: list of apparent resistivities computed from the conductivity models
phase_list: list of phases computed fromt the conductivity models
app_res_uncert: uncertainty (in percentage) of the apparent resistivity
phase_uncert: uncertainty (in degrees) of the phase
"""
fig, ax = plt.subplots(3, 1, figsize=(8, 3*3), dpi=400)

for i in range(len(sigma_list)):
# trickery to plot from node to node rather than at cell centers
z = np.repeat(mesh.vectorNx[1:-1], 2, axis=0)
z = np.r_[mesh.vectorNx[0], z, mesh.vectorNx[-1]]
sigma_plt = np.repeat(sigma_list[i], 2, axis=0)

# noise models, apparent resistivity is a percent, phase in units of degrees
app_res_noise = (app_res_uncert/100.)*np.abs(app_res_list[i])*np.random.randn(len(app_res_list[i]))
phase_noise = phase_uncert*np.random.randn(len(app_res_list[i]))

ax[0].plot(z, sigma_plt, 'C{}'.format(i), lw=2, label="{:1.0f} m".format(layer_thicknesses[i]))
ax[1].loglog(frequency, app_res_list[i]+app_res_noise, 'C{}-'.format(i), lw=1, label="{:1.0f} m".format(layer_thicknesses[i]))
ax[2].semilogx(frequency, phase_list[i]+phase_noise, 'C{}-'.format(i), lw=1, label="{:1.0f} m".format(layer_thicknesses[i]))
# conductivity plot
ax[0].set_ylabel("Conductivity (S/m)", fontsize=13)
ax[0].set_xlabel("Depth (m)", fontsize=14)
ax[0].set_xlim([-5000., 0.])
ax[0].set_ylim([0, 0.4])
ax[0].legend(fontsize=10, bbox_to_anchor=(1.18, 1.01))

# apparent resistivity plot
ax[1].set_ylabel("Apparent Resistivity $(\Omega m)$", fontsize=13)
ax[1].set_ylim([2e1, 2e2])

# phase plot
ax[2].set_ylabel("Phase $(^{\circ})$", fontsize=13)
ax[2].set_ylim([0., 90.])

# apparent resistivity and phase plots
for a in ax[1:]:
a.set_xlabel("Frequency (Hz)", fontsize=13)
a.set_xlim([frequency.min(), frequency.max()])

for a in ax:
a.grid(True, which='both', linewidth=0.3)
a.invert_xaxis()

plt.tight_layout()
plt.show()

In [30]:
plot_apparent_resistivities_and_phases(frequency, sigma_list, layer_thicknesses, app_res_list, phase_list)

### What happens if we add noise??¶

In [31]:
ipywidgets.interact(
lambda app_res_uncert, phase_uncert: plot_apparent_resistivities_and_phases(
frequency, sigma_list, layer_thicknesses,
app_res_list, phase_list, app_res_uncert, phase_uncert
),
app_res_uncert = ipywidgets.FloatSlider(
min=0., max=20., step=1, value=0., description="$\delta \\rho_a (\%)$"
),
phase_uncert = ipywidgets.FloatSlider(
min=0., max=10., step=0.5, value=0., description="$\delta \phi (^\circ)$"
)
)

Out[31]:
<function __main__.<lambda>>

When setting expectations for what we can hope to get out of an inversion, it is important to keep non-uniquness in mind... Without additional information (eg. additional geophysical data, well logs, etc) that constrain the physical properties and geometries of our geologic setting, we cannot necessarily expect to nail down both the conductivity and geometry of a given geologic unit.


text:
# non-split extension and Schur multiplier

Let $G$ be a central extension of the group $K$ by the group $H$. If we know that this extension is non-split, is it true that the order of $K$ must divide the Schur multiplier of the group $H$?

-
All the studying groups are finite. –  Tina Jul 4 '13 at 14:54
I suspect that you meant "does the order of $K$ divide the Schur multiplier of $K$"; here $K$ is the quotient, $H$ is the subgroup, the way you have described the setup. The answer is "No". Once you have found one such $G$, take $G\times A$ for any abelian group $A$. This is a non-split central extension of $K$ by $H\times A$. –  Alex B. Jul 4 '13 at 16:28
Dear Alex, thank you for your answer. I should mention that I mean $G/K=H$ and $G=K.H$, where $H$ is a non-abelian simple group. Does the order of $K$ divide the Schur multiplier of $H$? –  Tina Jul 4 '13 at 16:48
You are welcome! There was a typo in the first sentence of my comment, I meant "does the order of $H$ divide the Schur multiplier of $K$". –  Alex B. Jul 4 '13 at 16:50
Just a quick clarification: there are two standards for naming extensions, unfortunately. In one of them, "extension of $K$ by $H$" means that $K\triangleleft G$ and $G/K\cong H$. But in the other standard, "extension of $K$ by $H$" means that $H\triangleleft G$ and $G/H\cong K$ (hence Alex B.'s mistaking your meaning above). I would suggest always specifying which one you mean. –  Arturo Magidin Jul 4 '13 at 18:18

Just so we agree on the setup, you have an exact sequence $$1\rightarrow K\rightarrow G\rightarrow H\rightarrow 1,$$ where $K\leq Z(G)$, and you assume that the extension is non-split. You are asking whether it is true that $K$ divides the Schur multiplier of $H$.
The answer is "no", you can make $K$ arbitrarily large.
Indeed, having found one such extension, take $G\times A$ for an abelian group $A$. This is an extension of the form $$1\rightarrow K\times A\rightarrow G\times A\rightarrow H\rightarrow 1,$$ and clearly $K\times A\leq Z(G\times A)=Z(G)\times A$. Moreover, the extension still doesn't split, since the map $G\times A\rightarrow H$ factors through $G\rightarrow H$.
That's precisely why in the theory of Schur multipliers you have to demand in addition that $K\leq G'$, and the above example shows that you cannot drop this condition. The above example had $(G\times A)' = G'$. So we no longer had $K\times A\leq G'$, even if this was true for $K$ itself.
According to a paper by Dixon and Zalesski(Finite imprimitive linear groups of prime degree, J.Algebra (276)340-370), we have "Suppose that $H$ is a perfect group and consider $Z_p$ as a $ZH$-module with trivial action. Then $H^2(H,Z_p)\neq0$ that is, there exist non-split central extensions of $Z_p$ by $H$ if and only if $p$ divides the order of the Schur multiplier $M(H)$ of $H$." So in my question, if $H$ is a simple group and $K=Z_p$, then my conclusion is right. But, I am interested in the other cases for $Z_p$. Can I replace $Z_p$ with abelian groups or some other groups? –  Tina Jul 7 '13 at 9:59


text:
Exxamm.com-India’s most trusted exam preparation site for NDA | NEET | NCERT | ENGINEERING | DEFENCE | CBSE | MEDICAL | NTPC | CGL | RRB | SSC | BANKING | NCERT | OLYMPIAD
Q 1225812761.     A circular arc PQ of radius 1  subtends on angle of x radians, 0 < x < π as shown in the adjacent figure. The point R is the Intersection of the two tangent lines at points P and Q of the arc.O is the centre of the circle from which the arc is cut. Let us define the following functions as function of x:
S ( x ) = area of the sector OPQ

T ( x ) = area of the triangle PQR

U ( x ) = area of the shaded region

.The value of the limit  lim_(x→0^+) (U(x))/(T(x)) =

A

2/3

B

3/2

C

1

D

1/2

#### Access free resources including

• 100% free video lectures with detailed notes and examples
• Previous Year Papers
• Mock Tests
• Practices question categorized in topics and 4 levels with detailed solutions
• Syllabus & Pattern Analysis

• RRB-NTPC


text:
# What kind of document is suitable for a book if customer wants changes [closed]

I have a situation that I have to write a book from the material that my customer gave and she has not yet decided the layout of the book. What kind of document should I wrote that I can easily change for example

Paper size (A4 vs A5)

Font size and style

Margin sizes

Vertical space between lines

How many letters should the line on maximal length contain, such that the space between letters should be configurable?

-

## closed as off-topic by Martin Schröder, Tobi, Mico, Peter Jansson, ienisseiFeb 25 '14 at 16:53

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question does not fall within the scope of TeX, LaTeX or related typesetting systems as defined in the help center." – Martin Schröder, Tobi, Mico, Peter Jansson, ienissei
If this question can be reworded to fit the rules in the help center, please edit the question.

Unfortunately, there is no clear answer to this question, because this depends on the material to be shown in the book. – Christian Hupfer Feb 25 '14 at 15:02
Oh, well the material is mostly just text. Namely, I will write hand written letters to the book. It contain also a few fractions of the form $\frac{1}{3}$ and temperatures of the form $-15^\circ C$. – novicebookwriter Feb 25 '14 at 15:15
Well, I suggest scrbook, the design can be changed afterwards. – Christian Hupfer Feb 25 '14 at 15:18
@novicebookwriter, for temperaturs better use \SI{-15}{\celsius} with the siunitx package. – quinmars Feb 25 '14 at 15:20
Welcome to TeX.SX! You can have a look at our starter guide to familiarize yourself further with our format. – Martin Schröder Feb 25 '14 at 15:21


text:
## Machine Learning

A collection of articles on Machine Learning that teach you the concepts and how they’re implemented in practice.

Pairwise comparison of data-sets is very important. It allows us to compare two sets of data and decide whether: * one is better than the other, * one has more of some feature than the other, * the two sets are significantly different or not.

In the context of the weather data that you've been working with, we could test the following hypotheses: * Hypothesis 1: The mean temperature in China is greater than the mean temperature in Japan * Hypothesis 2: The mean humidity in Russia is greater than the mean humidity in Spain

There are many tests for rejecting the null hypothesis (that there is no difference between the two data-sets), which one you use depends on the properties of your data-set. Below is a brief list of some of the available tests:

Parametric Non-parametric
Paired t-test Wilcoxon signed-rank test
Unpaired t-test Mann-Whitney U test
Pearson correlation Spearman correlation
... ...

You can see they've been divided into two groups, parametric tests and non-parametric tests. Which group you use again depends on the properties of your data-set, in this case it's whether your data-sets follow a specific distribution. Typically, if your data-set doesn't follow a specific distribution, you want to use a non-parametric test. These tests don't make any assumptions about the distribution of your data.

Before we move onto running any of these tests - we need to determine a sufficient sample size. This is because we cannot make a reliable judgement if, for example, we only compare 4 temperature readings from Japan to 4 from China. The mean temperature of both countries may be the same - but was there enough data to make this claim?

## Sample Size Sufficiency

Let's get straight into it! We will be using numpy, matplotlib.pyplot, and pandas, so make sure to import these in the following cell. I've imported numpy for you. If you have access to it, I recommend also importing seaborn which automatically improves the quality of your plots. You may not have this in your environment so feel free to comment it out.

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

## Dataset

Now we want to load in our data-set from the CSV file. You can use the pd.read_csv() function. You will also need to add the following column names:

wt.columns = ['time', 'country', 'humidity', 'pressure', 'temp', 'max_temp', 'min_temp']

# load data into variable 'wt' here
wt = pd.read_csv("weather.csv", names=['time', 'country', 'humidity', 'pressure', 'temp', 'max_temp', 'min_temp'])

# output a sample of the data for inspection...
time country humidity pressure temp max_temp min_temp
0 15/03/14 NP 33.0 1020.00 297.15 297.15 297.15
1 15/03/14 VE 99.0 835.59 288.64 288.64 288.64
2 15/03/14 CN 12.0 613.65 280.79 280.79 280.79
3 15/03/14 TR 87.0 1021.00 283.42 284.15 283.15
4 15/03/14 DO 100.0 1020.65 295.14 295.14 295.14

## Sample size sufficiency

Let's check the sample size sufficiency for Hypothesis 1, for this we will need to pick out the temperatures for China and Japan, let's store them in variables named temp_JP and temp_CN respectively.

# variables holding the temperatures for Japan and China
temp_JP = wt[wt.country=='JP'].temp
temp_CN = wt[wt.country=='CN'].temp

So we have an idea of the number of samples for each data-set, let's print out the number of samples for inspection. There are man ways to do this, using .count(), .size, .shape, etc.

# print out the number of samples for inspection
print("number of temperature samples for Japan: {}".format(temp_JP.count()))
print("number of temperature samples for China: {}".format(temp_CN.count()))
number of temperature samples for Japan: 708
number of temperature samples for China: 282

Next, let's store the mean of temp_JP and temp_CN in variables named temp_JP_mean and temp_CN_mean respectively

# variables holding the mean temperature for Japan and China
temp_JP_mean = temp_JP.mean()
temp_CN_mean = temp_CN.mean()

# print out the means for inspection
print("mean temperature in Japan: {}".format(temp_JP_mean))
print("mean temperature in China: {}".format(temp_CN_mean))
mean temperature in Japan: 281.62871468926556
mean temperature in China: 288.6073404255319

It may be tempting at this point to state that the temperature in China is typically warmer than the temperature in Japan. But this could simply be by chance! We must be confident that we can make this claim with some statistical significance.

Now that we have the data-sets and their means, we can have a look at the relationship between sample sizes and the Standard Error of the Mean (SEM). The formula for this is:

$S{E}_{M}=\frac{s}{\sqrt{n}}$

We will implement this for all our possible sample sizes for Japan, and then China below.

SEM = []

for sample_size in range(3,len(temp_JP)):
sample = temp_JP[0:sample_size]
SEM.append(sample.std() / np.sqrt(sample_size))

plt.plot(SEM, color='magenta')
plt.plot(200,SEM[200], marker='o', color='cyan')

plt.ylabel("Standard Error of the Mean ($SE_M$)")
plt.xlabel("Sample size $(n)$")
plt.title("Relationship between $SE_M$ and $n$");
SEM = []

for sample_size in range(3,len(temp_CN)):
sample = temp_CN[0:sample_size]
SEM.append(sample.std() / np.sqrt(sample_size))

plt.plot(SEM, color='magenta')
plt.plot(200,SEM[200], marker='o', color='cyan')

plt.ylabel("Standard Error of the Mean ($SE_M$)")
plt.xlabel("Sample size $(n)$")
plt.title("Relationship between $SE_M$ and $n$");

In the two figures above we can see that we have selected sufficient sample sizes of 200 for Japan and China. This decision has been made visually based on the dispersion of the samples. There are many ways of doing this, but let's move onto using these samples to test our hypothesis!

## Normal Distribution Test

First, let's import scipy

from scipy import stats

Now let's check to see if our data-sets are normally distributed, so we know whether to use parametric or non-parametric tests to compare them. For this we will use the typical alpha, $\alpha =0.05$, as the cut-off for significance.

s, p = stats.normaltest(temp_JP) # p = p-value
if p < 0.05:
print('JP temperature data is not normal')
else:
print('JP temperature data is normal')
JP temperature data is not normal

Now let's check China.

s, p = stats.normaltest(temp_CN) # p = p-value
if p < 0.05:
print('CN temperature data is not normal')
else:
print('CN temperature data is normal')
CN temperature data is not normal

Let's plot the normalised histograms for the temperatures in Japan and China to look at this visually.

# plot histograms for temp_JP and temp_CN, make sure to set normed to True!
temp_JP.hist(color='magenta', density=True);
temp_CN.hist(color='cyan', alpha=0.5, density=True);

## Wilcoxon-signed Rank Test

Our data is supposedly not following a normal distribution, so we may decide to use a non-parametric test. There is much to consider when selecting a test to use, but we will use the Wilcoxon-signed rank test.

Using scipy.stats.wilcoxon, conduct a pairwise comparison of temperatures from Japan to temperatures from China to see if the null hypothesis is rejected.

s, p = stats.wilcoxon(temp_JP[0:200], temp_CN[0:200])

if p < 0.05:
print('null hypothesis rejected, significant difference between the data-sets')
else:
print('null hypothesis accepted, no significant difference between the data-sets')
null hypothesis rejected, significant difference between the data-sets

Now, let's use a quick if statement to resolve whether the temperature in our China sample is greater than our Japan sample. We can make this claim now that we have checked the significance of our results.

if(temp_CN[0:200].mean() > temp_JP[0:200].mean()):
print("Hypothesis 1 accepted, it is warmer in China")
else:
print("Hypothesis 1 rejected, it is warmer in Japan")
Hypothesis 1 accepted, it is warmer in China

Exercise

Copy this notebook and try testing Hypothesis 2, make sure your decision to use parametric or non-parametric testing is justified.

## Machine Learning

A collection of articles on Machine Learning that teach you the concepts and how they’re implemented in practice.


text:
# Stationary distribution of a DTMC that has recurrent and transient states

I need to construct a stationary distribution for $$\{X_n, n\in \mathcal{N}\}$$ with state space $$E=\{0,1,2,3,4,5\}$$ and transition matrix: $$\begin{bmatrix} 0.8 & 0.2 & 0&0&0&0\\ 0&0.1&0.4&0.5&0&0\\ 0&0.3&0.3&0.4&0&0\\ 0&0.6&0.3&0.1&0&0\\ 0&0&0&0&0.7&0.3\\ 0&0&0&0&0.5&0.5 \end{bmatrix}$$

In the previous steps of this problem I am asked to find the stationary distribution of all positive recurrent classes. For the positive recurrent class $$\{1,2,3\}$$ I was able to find the stationary distribution to be $$1/3$$ for each state, and for $$\{4,5\}$$ to be $$5/8$$ and $$3/8$$, respectively. However, I need to put these distributions together to find the stationary distribution for $$X_n$$ and I am not sure how to do so. This is just a practice problem as I study for my exam. Thank you for any help

States $$\{0,1,2,3\}$$ are isolated from $$\{4,5\}$$, so the stationary distribution isn't unique: any weighted combination of the stationary distributions for the subsystems will be stationary. The probabilities will evolve to a combination of the two, with weights equal to the total probability each subsystem started with. For instance, if the initial probabilities are $$(1/3,1/3,0,0,0,1/3),$$ the limiting probabilities are $$(0,2/9,2/9,2/9,5/24,3/24).$$
• Thank you and you were right about the limiting probabilities of $\{4,5\}$. I corrected my mistake. Dec 24, 2020 at 0:12