text:
# Understanding the differences: A direct comparison between Cython and C++ [closed]

Introduction:

The following two code samples are a direct comparison of performance between Cython and C++. I pieced them together with the help of others-- I barely understand how they work. Their speed in my computer is: 1.225s vs 3.989s (C++ being the faster). The purpose of this comparison is to decide which language I will be using for a series of small but computationally intensive scripts (~200-500 lines).

Question:

As a novice Python programmer dealing with data types, memory (de)allocation, and other low level processes for the very first time, what should I keep in mind when writing Cython modules? Is there anything in the Cython code that can be improved?

Cython script:

import cython
import numpy as np

@cython.boundscheck(False)
@cython.wraparound(False)

cdef class MathObject:
cdef double a

def __init__(self, a):
self.a = a

# chaotic function (simple math that can't be simplified out by the compiler).
cpdef double logMap(self):
self.a = 4.0 * self.a * (1.0 - self.a)
return 0

def get_value(self):
return self.a

def run():
cdef int i, ii, iii
cdef int one_thousand = 1000
cdef double x = 0.2  # logMap has complex behavior at this value

cdef MathObject objInstance

objArray = np.empty(one_thousand, dtype=MathObject)
for i in range(len(objArray)):
objArray[i] = MathObject(x)

cdef MathObject[:] objArrayiew = objArray

for i in range(one_thousand):
for ii in range(one_thousand):
for iii in range(one_thousand):
objInstance = objArrayiew[iii]
objInstance.logMap()
print objArray[i].get_value()

return None

C++ script:

#include <iostream>
using namespace std;

class MathObject {
public:

double a;
MathObject(double b) {
a=b;
};

// simple chaotic function
double logMap()
{
a = 4.0*a*(1.0-a);
} // logMap
};

int main(int argc, char *argv[])
{
int ii,jj,kk;

// create an array of 1000 objects
int N = 1000;
MathObject *objArray [N];  // 1000 pointers to objects
for (ii=0; ii<N; ++ii)  {
objArray[ii] = new MathObject(0.2);  // construct new objects
} // for ii

for (kk=0; kk<1000; ++kk)  {
for (jj=0; jj<1000; ++jj)  {
for (ii=0; ii<N; ++ii)  {
objArray[ii] -> logMap();
} // for ii
} // for jj
} // for kk

for (ii=0; ii<N; ++ii)  {
cout << objArray[ii] -> a << "\n";
} // for ii
}

(Note: I know that this question is a little vague, but most existing tutorials for Cython are written for an audience of C++ programmers. Beginner level information is astoundingly helpful!)

(Second Note: A good explanation of part of the Cython script can be found in the answer to this question.)

-

## closed as off-topic by jonrsharpe, Emily L., syb0rg, Edward, palacsintJun 21 '14 at 22:40

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

• "Questions must involve real code that you own or maintain. Questions seeking an explanation of someone else's code are off-topic. Pseudocode, hypothetical code, or stub code should be replaced by a concrete example." – jonrsharpe, Emily L., syb0rg, Edward, palacsint
If this question can be reworded to fit the rules in the help center, please edit the question.

Translating python to C++ (just a direct line for line translation).

import cython                                       #include <iostream>
import numpy as np                                  #include <vector>
#include "Range.h"
@cython.boundscheck(False)
@cython.wraparound(False)

### MathObject

cdef class MathObject:                              struct MathObject {
cdef double a                                       double a;

def __init__(self, a):                              MathObject(double a)
self.a = a                                          : a(a)
{}
cpdef double logMap(self):                          double logMap() {
self.a = 4.0 * self.a * (1.0 - self.a)              a = 4.0 * a * (1.0 - a);
return 0                                            return 0;
}
def get_value(self):                                double get_value() {
return self.a                                       return a;
}
};

### run

def run():                                          void run() {
cdef int i, ii, iii
cdef int one_thousand = 1000                        int one_thousand = 1000;
cdef double x = 0.2                                 double x = 0.2;

cdef MathObject objInstance

objArray = np.empty(one_thousand, dtype=MathObject) std::vector<MathObject> objArray(one_thousand, x);
for i in range(len(objArray)):
objArray[i] = MathObject(x)

cdef MathObject[:] objArrayiew = objArray

for i in range(one_thousand):                       for(int i: Range(one_thousand)) {
for ii in range(one_thousand):                      for(int ii: Range(one_thousand)) {
for iii in range(one_thousand):                     for(int iii: Range(one_thousand)) {
objInstance = objArrayiew[iii]                      MathObject& objInstance = objArray[iii];
objInstance.logMap()                                objInstance.logMap();
}
}
print objArray[i].get_value()                       std::cout << objArray[i].get_value() << "\n";
}

return None
}

### main

if __name__ == "__main__":                          int main() {
run()                                               run();
}

In C++ there is no built in Range type. But it is easy to write:

class Range
{
int beginValue;
int endValue;

struct RangeIterator
{
int value;
RangeIterator(int v) : value(v) {}
RangeIterator& operator++()     {++value;return *this;}
int            operator*()      {return value;}
bool           operator!=(RangeIterator const& rhs) const {return value != rhs.value;}
};

public:
typedef RangeIterator           iterator;
Range(int e)
: Range(0, e)
{}
Range(int b, int e)
: beginValue(b), endValue(e)
{}

iterator begin() { return RangeIterator(beginValue);}
iterator end()   { return RangeIterator(endValue);}
};

A couple of notes on the direct translation that I would change:

If you are not going to return a meaningful value (and we notice from usage that the return value is ignored).

double logMap() {
a = 4.0 * a * (1.0 - a);
return 0;
}

Then I would make this a void function so you know its not supposed to return anything.

void logMap() {
a = 4.0 * a * (1.0 - a);
}

I don't like geters they break encapsulation:

double get_value() {
return a;
}

You usually find that if you are getting values from an object doing stuff with it there is usally a member function (or member of the public API) that can do that. In this case you are only using the value to print. So why not just create a function that knows how to print the value.

friend std::ostream& operator<<(std::ostream& stream, MathObject const& data)
{
return stream << data.a;
}

If values don't change then mark them as const.

int    const one_thousand = 1000;
double const x            = 0.2;

-

C++ script:

• The data member a in MathObject should be private. This is necessary because classes are supposed to hide their data members and maintain their state from within the class.

• Instead of assigning variables in a constructor:

MathObject(double b) {
a=b;
};

use an initializer list:

MathObject(double a) : a(a) {}

This will also allow you to use the same name for the argument.

• Consider splitting main() into separate functions, rather than having it do everything. You could have a function for constructing the arrays, looping through them, and displaying them.

• Do not use a dynamically-allocated array when you can use an std::vector. These will do the memory management for you, which you've failed to do yourself (you use new to allocate memory but never free it with delete).

std::vector<MathObject> objArray(N, 0.2);

You also don't need to use pointers to objects; just use objects themselves.

• Loop counter variables should be initialized within the loop statements, not outside. At first, I thought those variables were for something else entirely.

• The nested for loops should be indented inside of each other, not stacked. It's more readable that way as it clearly shows that they're nested.

• Remove the "for" comments after each closing curly brace. They're not needed if the loop bodies are small and properly indented if nested.

-

Ignoring the Cython part of things for the moment, the C++ code you've posted is open to quite a bit of improvement, with respect to both speed and readability.

using namespace std;

This is generally considered a rather poor habit that it's best to break quickly (or, by strong preference, never form as a habit in the first place).

class MathObject {
public:

double a;
MathObject(double b) {
a=b;
};

You usually want to use a member initializer list instead of an assignment inside the body of the constructor when you can (and you can here). That gives code that looks more like this: MathObject(double b) : a(b) {}

// simple chaotic function
double logMap()
{
a = 4.0*a*(1.0-a);
} // logMap

Two points. First, this is defined as returning a double, but doesn't actually return a value at all. If you try to invoke it like double x = foo.logMap();, you get undefined behavior.

Second, this is the only function in the class other than the constructor. Since you're basically creating an object, and just using it to invoke this one function, it may make sense to consider implementing this as an overload of operator() so you can use this as a function object:

void operator()() {
a = 4.0*a*(1.0-a);
}

As a final point, the comment noting what a closing brace is closing can be helpful when the opening and closing brace are a long ways apart, but strike me as pretty pointless for a one-line function. About the only time I normally use them is on the closing brace of a namespace (particularly problematic because I don't normally indent all the code inside a namespace).

int main(int argc, char *argv[])

You don't seem to be using argc or argv in this case, so you can just as well leave them out, so the function header would be something like: int main()

int ii,jj,kk;

These don't strike me as the best possible variable names. Given that they're just loop indexes i, j and k are probably enough. If you want longer names, use the extra length to make the names more descriptive. Right now, you have the worst of both worlds: extra length and typing with zero gain in readability. As shown in the code below, a for loop index should normally be defined in the for loop in any case.

// create an array of 1000 objects
int N = 1000;
MathObject *objArray [N];  // 1000 pointers to objects

Two points here. First of all, using the built-in array type in C++ is another of those things that you want to avoid as a rule. You're almost always better off with an alternative, most often std::vector.

Second, avoid using pointers in C++. A lot of teachers and older books can make it seem as if nearly all C++ code needs to use pointers nearly everywhere, but this is actually false. In well-written C++, you'll see very few (if any) pointers at all. Likewise the concomitant use of new--most code shouldn't have to use new (directly) at all.

for (ii=0; ii<N; ++ii)  {
objArray[ii] = new MathObject(0.2);  // construct new objects
} // for ii

In this case, there's no real reason to use an array of pointers at all. I'd just use a vector of objects: std::vector<MathObject> objArray(N, 0.2);

for (kk=0; kk<1000; ++kk)  {
for (jj=0; jj<1000; ++jj)  {
for (ii=0; ii<N; ++ii)  {
objArray[ii] -> logMap();
} // for ii
} // for jj
} // for kk

With the change above from array of pointers to vector of objects, as well as adding the overload of operator() to the class, the syntax of the code in the loop changes a little, leaving you with something like this:

for (int i=0; i<1000; i++)
objArray[i]();

Likewise, where you print out the results:

for (ii=0; ii<N; ++ii)  {
cout << objArray[ii] -> a << "\n";
} // for ii

This can be changed to use a range-based for loop, much more like Python's, giving something like this:

for (MathObject &&m : objArray)
std::cout << m.a << "\n";

One more point: in this case, you're not really gaining much by using the class at all. Your entire code seems to work out to something like this:

double a = 0.2;

for (int k=0; k<1000; k++)
for (int j=0; j<1000; j++)
a = 4.0*a*(1.0-a);

std::cout << a < "\n";

Although you created 1000 separate objects, you seem to be starting them all from the same initial value and invoking each precisely the same number of times. The function may be chaotic, but the result is still deterministic, so you're basically just computing the same value 1000 times over.

Perhaps you intended to so something different, but at least as I read things, that seems to be equivalent to what you're actually doing. If I had to guess, it would be that what you intended was something more like this:

std::vector<double> a;

a.push_back(0.2);

for (int i=0; i<1000; i++)
a.push_back(a.back() * 4.0 * (1.0 - a.back());

This will compute and save the values from 1000 consecutive computations. Oh, although I forgot to mention it previously, also note the added white-space inside the computation. IMO, it improves readability somewhat.

-
+1 for the "don't need a class here." Even if he had 1000 different values, they can be stored in an array or vector of doubles, without all the OOP overhead, which I don't think adds anything here. – Jaime Jun 21 '14 at 16:24
@Jaime: What OOP overhead? Are you talking about source or object code. If you are talking source fine. If you are talking object there is none. The size of MathObject is the same as double. – Loki Astari Jun 21 '14 at 16:53
You can't remove the nested loops like that. Unfortunately the call to logMap has side affects and mutates the object. So you need to simulate 1000 * 1000 calls to each object in the array of 1000 elements. – Loki Astari Jun 21 '14 at 17:22
@LokiAstari: You may well be right--it was a guess and could easily be wrong. Maybe he intended to capture the last 1000 iterations, or maybe he intended to capture every thousandth iteration. As far as "OOP overhead" goes: calling a member function typically means loading the object's this pointer, then accessing the data via that pointer. Although it doesn't use extra data, it can still slow the function call. – Jerry Coffin Jun 21 '14 at 17:27
This was a simple case to test performance. At first, I wasn't sure if Cython would be able to handle this scenario at all! Function calls to a collection of classes is something I will be using a lot. It felt like a good idea for a test. Thank you for the C++ knowledge! – Paul Hill Jun 21 '14 at 17:49


text:
## Ramsey-like cardinals. II.(English)Zbl 1222.03055

Summary: This paper continues the study of the Ramsey-like large cardinals introduced in Part I [V. Gitman, ibid. 76, No. 2, 519–540 (2011; Zbl 1222.03054)] and in [I. Sharpe and P. D. Welch, Ann. Pure Appl. Logic 162, No. 11, 863–902 (2011; Zbl 1270.03071)]. Ramsey-like cardinals are defined by generalizing the characterization of Ramsey cardinals via the existence of elementary embeddings. Ultrafilters derived from such embeddings are fully iterable and so it is natural to ask about large cardinal notions asserting the existence of ultrafilters allowing only $$\alpha$$-many iterations for some countable ordinal $$\alpha$$. Here we study such $$\alpha$$-iterable cardinals. We show that the $$\alpha$$-iterable cardinals form a strict hierarchy for $$\alpha \leq \omega _{1}$$, that they are downward absolute to $$L$$ for $$\alpha <\omega _{1}^{L}$$, and that the consistency strength of Schindler’s remarkable cardinals is strictly between 1-iterable and 2-iterable cardinals.
We show that the strongly Ramsey and super Ramsey cardinals from Part I are downward absolute to the core model $$K$$. Finally, we use a forcing argument from a strongly Ramsey cardinal to separate the notions of Ramsey and virtually Ramsey cardinals. These were introduced in [Sharpe and Welch, loc. cit.] as an upper bound on the consistency strength of the Intermediate Chang’s Conjecture.

### MSC:

 3e+55 Large cardinals

### Citations:

Zbl 1222.03054; Zbl 1270.03071
Full Text:

### References:

 [1] Axiomatic set theory (Proceedings of the Symposium on Pure Mathematics, vol. XIII, part II, University of California, Los Angeles, California, 1967) pp 33– (1974) [2] DOI: 10.1016/0168-0072(90)90028-Z · Zbl 0719.03025 [3] DOI: 10.1016/0168-0072(89)90030-4 · Zbl 0714.03048 [4] Coding the universe 47 (1982) · Zbl 0468.03031 [5] Inner models and large cardinals 5 (2002) [6] Annals of Pure and Applied Logic [7] Ramsey-like cardinals 76 pp 519– (2011) · Zbl 1222.03054 [8] DOI: 10.1002/malq.200410002 · Zbl 1058.03052 [9] Ramsey cardinals and constructibility 44 pp 260– (1979) [10] Annals of Pure and Applied Logic 1 pp 179– (1970) [11] The higher infinite (2003) [12] Saturated ideals 43 pp 65– (1978) [13] Chains of end elementary extensions of models of set theory 63 pp 1116– (1998) · Zbl 0915.03034
This reference list is based on information provided by the publisher or from digital mathematics libraries. Its items are heuristically matched to zbMATH identifiers and may contain data conversion errors. It attempts to reflect the references listed in the original paper as accurately as possible without claiming the completeness or perfect precision of the matching.

format: plain


text:
# Percentage - Excercise 1

## Questions

1.   Entry fee in an exhibition was Rs. 1. Later, this was reduced by 25% which increased the sale by 20%. The percentage increase in the number of visitors is :

• 20.00%
• 40.00%
• 60.00%
• 80.00%

Let the total original sale be Rs. 100.
Then, original number of visitors = 100.
New number of visitors = ${\frac{120}{0.75}}$ = 160.
Increase %= 60%

2.   Helpers are needed to prepare for the fete. Each helper can make either 2 large cakes or 35 small cakes per hour. The kitchen is available for 3 hours and 20 large cakes and 700 small cakes are needed. How many helpers are required?

• 10
• 15
• 20
• 25

20 large cakes will require the equivalent of 10 helpers working for one hour.
700 small cakes will require the equivalent of 20 helpers working for one hour.
This means if only one hour were available we would need 30 helpers.
But since three hours are available we can use 10 helpers.

3.   The current birth rate per thousand is 32, whereas corresponding death rate is 11 per thousand. The net growth rate in terms of population increase in percent is given by :

• 1.20%
• 2.10%
• 3.10%
• 1.30%

Net growth on= ${1000=(32-11)\frac{}{}}$= 21
. Net growth on =${\left(100=\begin{array}{c}\frac{21}{1000}\times100\end{array}\right)}$= 2.1%.

4.   The price of a cycle is reduced by 25 per cent. The new price is reduced by a further 20 per cent. The two reductions together are equal to a single reduction of

• 45.00%
• 40.00%
• 35.00%
• 32.50%

Let the original price of the cycle be 100.
After the first reduction the price will be 75.
This new price is then reduced by ${\frac{20}{100}=0.8\times75}$
60 represents a reduction of 40 percent on the original.

5.   If the price of sugar rises from Rs. 6 per kg to Rs. 7.50 per kg, a person, to have no increase in his expenditure on sugar, will have to reduce his consumption of sugar by

• 15
• 20
• 25
• 30

Let original consumption = 100 kg and new consumption = x kg,
So, $(100\times6)=(x\times7.50)$
=x =80 kg
Reduction in consumption = 20%

6.   An agent, gets a commission of 2.5% on the sales of cloth. If on a certain day, he gets Rs. 12.50 as commission, the cloth sold through him on that day is worth

• 300
• 500
• 700
• 900

Let the total sale be Rs. x
Then 2.5%. of x =12.50
=(25/10$\times$1/100$\times$x) = 125/10
x = 500

7.   A dress on sale in a shop is marked at $D. During the discount sale its price is reduced by 15%. Staff are allowed a further 10% reduction on the discounted price. If a staff member buys the dress what will she have to pay in terms of D ? • 0.75D • 0.76D • 0.765D • None of the above The Correct answer is :0.765D If the price is reduced by 15 %, then the new price will be 0.85D nIf this new price is further reduced by 10% The discounted price will be=${(0.9\times0.85D)\frac{}{}}$= 0.765D 8. What percent decrease in salaries would exactly cancel out the 20 percent increase ? • (16 + 2/3) % • (15 + 2/3)% • (14 +1/3)% • (13 + 4/5)% The Correct answer is :(16 + 2/3) % Let the orginal salary = Rs. 100. New's salary = Rs.120. Decrease on 120 = 20. Decrease on 100 =(20/120$\times$100) = 16+2/3 9. The difference of two numbers is 20% of the larger number, if the smaller number is 20, then the larger number is : • 15 • 25 • 35 • 45 The Correct answer is :25 Let the large number be x. nThen${(x-20)=\frac{20x}{100}}{x-\frac{x}{5}}$= 20 x = 25. 10. To a sugar solution of 3 litres containing 40% sugar, one litre of water is added. The percentage of sugar in the new solution is : • 20 • 30 • 40 • 50 The Correct answer is :30 Quantity of sugar =${\left(\begin{array}{c}\frac{40}{100}\times3\end{array}\right)}$kg = 1.2 kg. New percentage =${\left(\begin{array}{c}\frac{1.2}{4}\times100\end{array}\right)}$= 30%. 11. In an examination, 34% of the students failed in Mathematics and 42% failed in English. If 20% of the students failed in both the subjects, then the percentage of students who passed in both the subjects was : • 34.00% • 44.00% • 54.00% • 64.00% The Correct answer is :44.00% Failed in mathematics, n(A) = 34 Failed in English, n(B) = 42 n(A∪B)=n(A)+n(B)−n(A∩B)=34+42−20=56 Failed in either or both subjects are 56 Percentage passed = (100−56)%=44% 12. How many litres of pure acid are there in 8 litres of a 20% solution ? • 1.4 • 1.5 • 1.6 • 1.7 The Correct answer is :1.6 Quantity of pure acid =${\left(\begin{array}{c}\frac{20}{100}\times8\end{array}\right)}$= 1.6 litres. 13. Ramesh started a software business by investing Rs. 50,000. After six months, Rajesh joined him with a capital of Rs. 75,000. After 3 years, they earned a profit of Rs. 27,000. What was Rajesh's share in the profit? • 10000 • 15000 • 20000 • 25000 The Correct answer is :15000 Ramesh : Rajesh =(50000 * 36):(25000 * 30)=4 : 5. Therefore, Rajesh's share Rs=(27000 * 5/9) =Rs.15,000. 14. In an election only two candidates contested 20% of the voters did not vote and 120 votes were declared as invalid. The winner got 200 votes more than his opponent thus he secured 41% votes of the total voters on the voter list. Percentage votes of the defeated candidate out of the total votes casted is: • 47.50% • 41.00% • 38.00% • 45.00% The Correct answer is :45.00% Let there be x voters and k votes goes to loser then 0.8x - 120 = k + (k + 200) k+200 = 0.41x k = 1440 and (k + 200) =1640 Therefore${(1440)\times(3200)\times(100)\frac{}{}}\$ = 45%

15.   600 students took the test on physics and chemistry. 35% students failed in physics and 45% students failed in chemistry and 40% of those who passed in chemistry also passed in physics, then how many students failed in both:

• 162
• 138
• 60
• None of these
The Correct answer is :None of these

Physics =>Failed= 35% Passed => 65%
Chemistry=>Failed=45%=> Passed=55%
Passed in both = 22% of total student
Percentage of students who are passed in any of the physics or chemistry or both =(65+55) - 22 = 98%
>So, the percentage of students who are failed in both = 2%
Therefore, total failed (in both subjects) students = 12


text:
# Citing this Document

This document can be cited as follows:

# Getting Started

This markdown document provides code and detailed examples associated with the Computer Applications and Quantitative Methods in Archaeology Workshop hosted by Matt Peeples and Tom Brughmans on March 13th, 2017 in Atlanta at Georgia State University. This document assumes you have a basic knowledge of network science terminology and at least some basic familiarity with R and R Studio. If you are new to R, we recommend you first run through the brief Code School examples focused on R available here (http://tryr.codeschool.com/). The examples below can be run by simply copying and pasting the code below but modifications will require a bit of additional experience. For those new to network science in general, we recommend you start by reading the Brughmans (2014) review article (see recommended bibliography for more details).

# Getting Our Data into R

For the purposes of this workshop we will be using one real archaeological dataset as an example (and some simulated data). These are data from research Peeples has conducted and can be used to demonstrate some of the most common methods for archaeological network analysis. These data come from the Southwest Social Networks Project (http://www.southwestsocialnetworks.net) and represent counts of decorated ceramic wares for one 50 year interval (A.D. 1050-1100) from a series of 223 settlements with public architectural features located across the northern U.S. Southwest (Chacoan Great Houses and Great Kivas for those who like such things). This dataset is divided into two .csv files, the first with the ceramic counts by site and a second with the attributes of those sites, including their locations. To ensure the security of site locational information as required by the state data repositories charged with maintaining these data, we have randomly relocated each settlement between 7-10 kilometers from their actual locations. The code below imports these data from the .csv files into objects we can use in R.

Right click and choose “save as” to download the ceramic data and attribute data.

Before running the code below, however, we need to ensure that our R session is set to the correct working directory (the location where you placed the .csv files for this workshop). To do that, go to the menu bar at the top and click Session > Set Working Directory > Choose Directory and navigate to the place on your hard drive where these files reside (alternatively you can hit Ctrl + Shift + H and then navigate to the appropriate directory). Once you have done that, you will be able to execute the code below by simply copying the text and then pasting it into the R console.

Let’s start by importing our ceramic and attribute data.

# Import chaco ceramic data into an object named chaco. row.names=1 sets the
# first column values to the row names for each item.
chaco <- read.csv(file = "AD1050cer.csv", row.names = 1)
# Import attribute data into an object named chaco.attr
chaco.attr <- read.csv("AD1050attr.csv", row.names = 1)

# Installing and initializing the necessary libraries

For the purposes of this workshop, we will rely on a few pre-existing R packages. In order to use these packages in a new installation of R and R-studio, we first need to install them. Note that you will only need to do this once on each new installation of R/R-Studio. To install packages, you can click on the “Packages” tab in the window in the bottom right of R studio, then click the “Install” button at the top and type the names of the packages separated by commas. Alternatively you can install packages from the console by simply typing “install.packages(‘nameofpackagehere’)” without the quotation marks.

For the purposes of this workshop we will rely on the statnet and tnet packages which also install several other required packages like sna, network, and igraph. Copy the line below into the console or follow the instructions to install using the packages tab:

install.packages(‘statnet’,‘tnet’)

Once you have installed these packages we use the library console command to initialize our packages.

library(statnet)
library(tnet)

# Defining custom functions for our analyses

In addition to the pre-existing R packages, there are also a few additional functions that we will need to define manually for several of our analyses below. Specifically, we want to be able to take our ceramic frequency data and convert that into a symmetric similarity or distance matrix that we will use to define and weight our networks. For the examples below, we will rely on three methods for defining these matrices: 1) Co-presence of types, 2) Brainerd-Robinson similarity, and 3) $$\chi^{2}$$ distances. These are only a few among a wide variety of options.

#### Co-presence

The first simple measure we will use for defining networks here is the number of categories that are co-present at pairs of sites. We can calculate this measure through simple matrix algebra. First, we create a binary incidence matrix of ceramic counts by site by simply dichotomizing our count data (a matrix with 1 for all categories that are present and 0 elsewhere). We can define that incidence matrix as $$A$$ and the transpose (switching rows and columns) of that matrix as $$A^{T}$$ we can find the number of categories that overlap between sites $$P$$ as: $P=A * A^{T}$

The result of this procedure is a symmetric matrix with the number of rows and columns determined by the number of nodes comparing each site to every other site with each cell represents counts of co-occurrence between pairs of sites. The diagonal of this matrix represents the total number of categories present for the site denoted by that row. In many cases, we may not want to count rare occurrences as “present” and may instead want to create a threshold absolute value or proportion that a given category must exceed to be counted as “present” for creating our matrix. In the example we will use here today a category must represent at least 10% of a row to be considered “present” in this calculation. The code below could be easily modified to use a different cutoff and the choice of cutoff is an important decision with substantive impacts on our networks.

The chunk of code below defines a procedure for calculating simple co-presence for all categories representing more than 10% of a given row. We then create a new object called “chacoP” which represents this matrix of co-occurrence from site to site.

co.p <- function(x, thresh = 0.1) {
# create matrix of proportions from chaco
temp <- prop.table(as.matrix(x), 1)
# define anything with greater than or equal to 0.1 as present (1)
temp[temp >= thresh] <- 1
# define all other cells as absent (0)
temp[temp < 1] <- 0
# matrix algebraic calculation to find co-occurence (%*% indicates matrix
# multiplication)
out <- temp %*% t(temp)
return(out)
}

# run the function
chacoP <- co.p(chaco)

#### Brainerd-Robinson Similarity

The next metric we will use here is a rescaled version of the Brainerd-Robinson (BR) similarity metric. This BR measure is commonly used in archaeology including in a number of recent (and not so recent) network studies. This measure represents the total similarity in proportional representation of categories and is defined as:

$S = {\frac{2-\sum_{k} \left|x_{k} - y_{k}\right|} {2}}$

where, for all categories $$k$$, $$x$$ is the proportion of $$k$$ in the first assemblage and $$y$$ is the proportion of $$k$$ in the second. This provides a scale of similarity from 0-1 where 1 is perfect similarity and 0 indicates no similarity. The chunk below defines the code for calculating this modified BR similarity measure. Note that this code could certainly be reduced and streamlined in many ways but I have chosen to use the “expanded” version for the sake of clarity for those new to R. By default, this code rounds the resulting similarity metric to three digits past the decimal for ease of viewing the output. This chunk ends by running this function for our sample dataset defining a new object “chacoBR” with the resulting similarity matrix.

sim.mat <- function(x) {
# get names of sites
names <- row.names(x)
x <- na.omit(x)  # remove any rows with missing data
x <- prop.table(as.matrix(x), 1)  # convert to row proportions
rd <- dim(x)[1]
# create an empty symmetric matrix of 0s
results <- matrix(0, rd, rd)
# the following dreaded double for loop goes through every cell in the
# output data table and calculates the BR value as descried above
for (s1 in 1:rd) {
for (s2 in 1:rd) {
x1Temp <- as.numeric(x[s1, ])
x2Temp <- as.numeric(x[s2, ])
results[s1, s2] <- 2 - (sum(abs(x1Temp - x2Temp)))
}
}
row.names(results) <- names  # assign row names to output
colnames(results) <- names  # assign column names to output
results <- results/2  # rescale results between 0 and 1
results <- round(results, 3)  # round results
return(results)
}  # return the final output table

# Run the function defined above on our sample data
chacoBR <- sim.mat(chaco)

#### $$\chi^{2}$$ Distance

The next measure we will use is the $$\chi^{2}$$ distance metric which is the basis of correspondence analysis and related methods commonly used for frequency seriation in archaeology (note that this should probably really be called the $$\chi$$ distance since the typical form we use is not squared, but the name persists this way in the literature so that’s what we use here). This measure is defined as:

$\chi_{jk} = \sqrt{\sum \frac 1{c_{j}} ({x_{j}-y_{j})^{2}}}$

where $$c_j$$ denotes the $$j_{th}$$ element of the average row profile (the proportional abundance of $$j$$ across all rows) and $$x$$ and $$y$$ represent row profiles for the two sites under comparison. This metric therefore takes raw abundance (rather than simply proportional representation) into account when defining distance between sites. The definition of this metric is such that rare categories play a greater role in defining distances among sites than common categories (as in correspondence analysis). This measure has a minimum value of 0 and no theoretical upper limit.

The code for calculating $$\chi^{2}$$ distances is defined in the chunk below and a new object called “chacoX” is created using this measure. It is sometimes preferable to rescale this measure so that it is bounded between 0 and 1. We create a second object called “chacoX01” which represents rescaled distances by simply dividing the matrix by the maximum observed value (there are many other ways to convert a distance to a similarity measure but this will be fine for our current purposes).

chi.dist <- function(x) {
rowprof <- x/apply(x, 1, sum)  # calculates the profile for every row
avgprof <- apply(x, 2, sum)/sum(x)  # calculates the average profile
# creates a distance object of $\chi^{2}$ distances
chid <- dist(as.matrix(rowprof) %*% diag(1/sqrt(avgprof)))
# return the reults
return(as.matrix(chid))
}

# Run the script and then create the rescaled 0-1 version
chacoX <- chi.dist(chaco)
chacoX01 <- chacoX/max(chacoX)

# Creating and plotting network objects from our similarity/distance matrices

Now that we have defined our three measures of similarity or distance, the next step is to convert these into network objects that our R packages will be able to work with. We can do this by either creating binary networks (where ties are either present or absent) or weighted networks (which in many cases are simply the raw similarity/distance matrices we calculated above). We will provide examples of both approaches, starting with simple binary networks. There are many ways to define networks from matrices like those we generated above and our examples below should not been seen as an exhaustive set of approaches.

#### Creating binary network objects

First, we will create binary networks using our ceramic ware co-occurrence matrix. Sites that share a single co-occurrence will be defined as connected here. We then use the “network” function to create a network object using the argument “directed=F” to let the function know that our data are from a symmetric matrix and that ties always extend in both directions between pairs of linked nodes. See help(network) for more details on options here.

After we create this new network object we plot it first using the default graph layout (Fruchterman-Reingold - We’ll discuss what this is in the workshop) and then based on the geographic location of our nodes. We’re not going to spend a lot of time exploring network plotting procedures in this workshop, but for those of you who may be interested in making far prettier graphs than those used here for demonstration purposes, I would recommend exploring the new “ggraph” package. (https://github.com/thomasp85/ggraph)

# create network object from co-occurrence
Pnet <- network(chacoP, directed = F)
# Now let's add names for our nodes based on the row names of our original
# matrix
Pnet %v% "vertex.names" <- row.names(chacoP)
# look at the results
Pnet
##  Network attributes:
##   vertices = 223
##   directed = FALSE
##   hyper = FALSE
##   loops = FALSE
##   multiple = FALSE
##   bipartite = FALSE
##   total edges= 16042
##     missing edges= 0
##     non-missing edges= 16042
##
##  Vertex attribute names:
##     vertex.names
##
##  Edge attribute names not shown
# set up for plotting two plots, side by side by setting the pot to 1 row, 2
# columns
par(mfrow = c(1, 2))
# plot network using default layout
plot(Pnet, edge.col = "gray", edge.lwd = 0.25, vertex.cex = 0.5, main = "co-presence network")
# plot network using geographic coordinates
plot(Pnet, edge.col = "gray", edge.lwd = 0.25, vertex.cex = 0.5, coord = chaco.attr[,
9:10])

par(mfrow = c(1, 1))  # return to single plotting mode

The next chunk of code will produce a network based on our BR similarity matrix. In this example, we define ties as present between pairs of sites when they share more than 75% commonality (BR > 0.75) in terms of the proportions of ceramic wares recovered from both sites in a dyad.

In the code below, the event2dichot function (from the statnet/network package) takes our matrix and divides it into 1s and 0s based on the cut off we choose. Here we’re using and ‘absolute’ cut off meaning we’re assigning a specific value to use as the cut off (0.75) and defining all similarity values higher than that cutoff as 1 and all other dyads as 0. We then send the output of this function to the network function just as before. After examining our new network we then plot it both using a graph layout and geographic locations.

# Define our binary network object from BR similarity
BRnet <- network(event2dichot(chacoBR, method = "absolute", thresh = 0.75),
directed = F)
# Now let's add names for our nodes based on the row names of our original
# matrix
BRnet %v% "vertex.names" <- row.names(chacoBR)
# look at the results.
BRnet
##  Network attributes:
##   vertices = 223
##   directed = FALSE
##   hyper = FALSE
##   loops = FALSE
##   multiple = FALSE
##   bipartite = FALSE
##   total edges= 7245
##     missing edges= 0
##     non-missing edges= 7245
##
##  Vertex attribute names:
##     vertex.names
##
##  Edge attribute names not shown
par(mfrow = c(1, 2))  # set up for plotting two plots, side by side
# plot network using default layout
plot(BRnet, edge.col = "gray", edge.lwd = 0.25, vertex.cex = 0.5, main = "BR network")
# plot network using geographic coordinates
plot(BRnet, edge.col = "gray", edge.lwd = 0.25, vertex.cex = 0.5, coord = chaco.attr[,
9:10])

par(mfrow = c(1, 1))  # return to single plotting mode

In the next chunk of code we will use the $$\chi^2$$ distances to create binary networks. This time, we will not use an absolute value to define ties as present, but instead will define those similarities (1-distances) greater than 80 percent of all similarities as present. We will then once again plot just as above.

# Note we use 1 minus chacoX01 here so to convert a distance to a similarity
Xnet <- network(event2dichot(1 - chacoX01, method = "quantile", thresh = 0.8),
directed = F)
# Once again add vertext names
Xnet %v% "vertex.names" <- row.names(chacoX01)
# look at the results
Xnet
##  Network attributes:
##   vertices = 223
##   directed = FALSE
##   hyper = FALSE
##   loops = FALSE
##   multiple = FALSE
##   bipartite = FALSE
##   total edges= 4861
##     missing edges= 0
##     non-missing edges= 4861
##
##  Vertex attribute names:
##     vertex.names
##
##  Edge attribute names not shown
par(mfrow = c(1, 2))  # set up for plotting two plots, side by side
# plot network using default layout
plot(Xnet, edge.col = "gray", edge.lwd = 0.25, vertex.cex = 0.5, main = "Chi-squared network")
# plot network using geographic coordinates
plot(Xnet, edge.col = "gray", edge.lwd = 0.25, vertex.cex = 0.5, coord = chaco.attr[,
9:10])

par(mfrow = c(1, 1))  # return to single plotting mode

#### Creating and plotting weighted network objects

It is also possible to use R to create weighted networks where individual edges are valued. We have found that this works reasonably well with networks of co-presence or something similar (counts of mentions in texts or monuments for example) but this does not perform well when applied to similarity or distance matrices (because every possible link has a value, however small, so the network gets unwieldy very fast). In the latter case, we have found it is better to just work directly with the underlying similarity/distance matrix.

Creating a weighted network object in R is easy and only requires a slight modification from the procedure above. In the chunk of code below, we will simply add the arguments “ignore.eval=F” and “names.eval=‘weight’” to let the network function know we would like weights to be retained and we would like that attribute called ‘weight’. We will apply this to the matrix of co-presence defined above and then plot the result showing the weights of individual ties. Although it is difficult to tell with a network this size, the lines defining the edges are scaled based on their edge weights. We further color code the nodes based on the “CSN_macro_group” (sub-region) variable in our attribute file. Although we only show this approach for the co-presence network, these same arguments work with the other matrices defined above.

# create weighted network object from co-occurrence matrix by adding the
# ignore.eval=F argument
Pnet2 <- network(chacoP, directed = F, ignore.eval = F, names.eval = "weight")
Pnet2 %v% "vertex.names" <- row.names(chacoP)
Pnet2
##  Network attributes:
##   vertices = 223
##   directed = FALSE
##   hyper = FALSE
##   loops = FALSE
##   multiple = FALSE
##   bipartite = FALSE
##   total edges= 16042
##     missing edges= 0
##     non-missing edges= 16042
##
##  Vertex attribute names:
##     vertex.names
##
##  Edge attribute names not shown
par(mfrow = c(1, 2))  # set up for plotting two plots, side by side
# create object to add color by subregion
reg.col <- as.factor(chaco.attr\$CSN_macro_group)
# plot weighted network using default layout
plot(Pnet2, edge.col = "gray", edge.lwd = "weight", vertex.cex = 0.5, vertex.col = reg.col,
main = "co-presence weighted network")
# plot weighted network using geographic coordinates
plot(Pnet2, edge.col = "gray", edge.lwd = "weight", vertex.cex = 0.5, coord = chaco.attr[,
9:10], vertex.col = reg.col)

par(mfrow = c(1, 1))  # return to single plotting mode

Now for the sake of visualization, we’ll create and plot a smaller random weighted network to show how edge weights are depicted. There are many more plotting options within R base graphics and ggplot/ggraph and we direct you to the help documents for those packages for more details.

set.seed(45645)  # set random seed for repeatability
num_nodes <- 40  # set the number of nodes
num_edges <- 80  # set the number of edges
# create a blank 'edgelist' for this network
edgelist <- matrix("", nrow = num_edges, ncol = 2)
# populate the edgelist with random node number assignments
for (i in 1:num_edges) edgelist[i, ] <- sample(x = seq(1, num_nodes), size = 2,
replace = F)
# convert to network object (empty)
net2 <- network.initialize(num_nodes)
# define ties based on edglist as 1
net2[as.matrix(edgelist)] <- 1
# set edge weights using a random uniform draw
edge_weights <- round(runif(num_edges, min = 1, max = 5))
# apply edge weights to network object
set.edge.value(net2, "weight", edge_weights)
# create a scale for color coding edges from gray to darkblue
edge.cols <- colorRampPalette(c("gray", "darkblue"))(5)
# plot network and add legend (get.edge.value uses the values of weight to
# color code and define line thickness
plot(net2, usearrows = F, edge.lwd = get.edge.value(net2, "weight"), edge.col = edge.cols[get.edge.value(net2,
"weight")], main = "random weighted network")
legend("bottomright", legend = seq(1, 5), lwd = seq(1, 5), col = edge.cols)

# Calculating graph-level and node-level metrics

One of the most common kinds of analysis for archaeological and other networks is the calculation of measures of node/edge centrality and graph centralization. There are many different measures in the literature each appropriate for different kinds of research questions and data formats (see Borgatti and Everett 2006). We will not cover the interpretation of these network metrics in depth in this workshop, but we will be using several common measures as our means of assessing the impact of missing data and other kinds of uncertainty on our interpretations of archaeological networks.

In this section, we briefly describe the network metrics we will use and then define a function to easily calculate multiple measures simultaneously on multiple networks. The primary measures we will use are the binary and weighted versions of: 1) degree centrality, 2) betweenness centrality, and 3) eigenvector centrality. We direct readers to Peeples and Roberts (2013) and especially the online supplemental materials for that article for more details on the calculation of these metrics.

Degree centrality for a node is defined as the total number of direct connections in which that node is involved. In weighted networks, this is simply the total weight of all connections for that node (minus 1 to remove self-loops). Betweenness centrality is defined as the number of shortest paths between pairs of nodes in a network involving the target node divided by the total number of shortest paths in the network as a whole. For binary networks, shortest paths are defined as the smallest number of direct ties that must be crossed to get from one specific node to another. Calculating betweenness centrality for weighted networks is a bit more complicated. The method we use here come from Opsahl and others (2010). This approach defines shortest paths as the “path of least resistance” between pairs of target nodes. In other words, the stronger the path (the higher the edge weights) the “shorter” it is according to this algorithm. Refer to Peeples and Roberts (2013) for more details. Eigenvector centrality is a measure of a node’s importance in a network defined in relation to the first eigenvector of the adjacency matrix of nodes for both binary and weighted networks. For both binary and weighted networks, a node’s eigenvector centrality score is proportional to the summed scores of other nodes to which it is connected. In other words, eigenvector centrality for a node will increase if a node is either connected to lots of other nodes, or if a node is connected to highly central nodes. We rescale this measure as is commonly seen in the network science literature so that the sum of squared scores is equal to the total number of nodes in that network.

The chunk of code below contains two functions: 1) the first calculates the three measures of centrality described above for binary networks and 2) the second calculates these measures for weighted networks (or similarity matrices). The results are by default returned in a matrix with a column for each measure and a row for each node. Note that the statnet/sna package uses some of the same names for functions as tnet/igraph. Thus, we must specify which package we mean to use in the code below (e.g., sna::degree means the degree calculation method used in the sna package which is initialized through statnet. If we wished to use the igraph version we would instead use igraph::degree).

# Calculate centrality scores for binary networks
net.stats <- function(y) {
# calculate degree centrality
dg <- as.matrix(sna::degree(y, gmode = "graph"))
# calculate and scale eigenvector centrality
eg <- as.matrix(sna::evcent(y))
eg <- sqrt((eg^2) * length(eg))
# calculate betweenness centrality
bw <- sna::betweenness(y, gmode = "graph")
# combine centrality scores into matrix
output <- cbind(dg, eg, bw)
rownames(output) <- rownames(as.matrix(y))
colnames(output) <- c("dg", "eg", "bw")
return(output)
}  # return results of this function

# Calculate centrality scores for weighted networks (similarity matrices)
net.stats.wt <- function(y) {
# calculate weighted degree as the sum of weights - 1
dg.wt <- as.matrix(rowSums(y) - 1)
# calculate weighted eigenvector centrality and rescale
eg.wt <- as.matrix(sna::evcent(y))
eg.wt <- sqrt((eg.wt^2) * length(eg.wt))
# calculate weighted betweenness from the tnet package (we use the
# suppressWarnings package to avoid notifications)
bw.wt <- suppressWarnings(betweenness_w(y, directed = F))[, 2]
output <- cbind(dg.wt, eg.wt, bw.wt)
rownames(output) <- rownames(as.matrix(y))
colnames(output) <- c("dg.wt", "eg.wt", "bw.wt")
return(output)
}  # return results of this function

Now let’s calculate these measures for our networks defined above and look at the first couple of rows for the final example based on the $$\chi^2$$ network.

# net stats for binary co-presence network
co.p.stats <- net.stats(Pnet)
# net stats for binary BR similarity network
BR.stats <- net.stats(BRnet)
# net stats for binary X^2 similarity network (1-distance)
X.stats <- net.stats(Xnet)
head(X.stats)
##                dg          eg        bw
## Aztec           3 0.017140077  89.73043
## Pueblo Pintado 92 1.638913576  70.11709
## Morris 39      37 0.011199709 947.53069
## Holmes Group   37 0.011199709 947.53069
## Jaquez         13 0.066034620 925.91999
## Fort Wingate    5 0.002894922 284.80973

And now let’s calculate the weighted versions.

# net stats for weighted co-presence network
co.pw.stats <- net.stats.wt(chacoP)
# net stats for weighted BR similarity network
BR.stats.w <- net.stats.wt(chacoBR)
# net stats for X^2 similarity (1-distance)
X.stats.w <- net.stats.wt(1 - chacoX01)
head(X.stats.w)
##                   dg.wt     eg.wt bw.wt
## Aztec          180.9051 1.0567527     0
## Pueblo Pintado 186.9819 1.0961256     0
## Morris 39      179.5033 1.0495716     0
## Holmes Group   179.4261 1.0490867     0
## Jaquez         182.8987 1.0688679     0
## Fort Wingate   167.2580 0.9757783     0

It is also often informative to evaluate graph-level measures of centralization. In the code below, we calculate centralization measures associated with the node-level centrality measures described above. These measures are essentially a measure of how central all nodes in our network are in relation to a theoretical maximum. We can calculate most measures using existing functions but we need to create a custom function for calculating weighted betweenness centralization. Again Peeples and Roberts (2013) provide more details.

## calculate centralization measures for binary network
cent.bin <- function(net) {
output <- matrix(0, 1, 3)
colnames(output) <- c("degree", "between", "eigen")
# calculate binary degree centralization
output[1, 1] <- centralization(net, sna::degree, normalize = T)
# calculate binary eigenvector centralization
output[1, 2] <- centralization(net, sna::evcent, normalize = T)
# calculate binary betweenness centralization
output[1, 3] <- centralization(net, sna::betweenness, normalize = T)
return(output)
}

# define function for calculating weighted betweenness centralization
bw.cent <- function(x) {
Cstar <- max(x)
Csum <- (Cstar - x)
num <- 2 * (sum(Csum))
den <- ((length(x) - 1)^2) * (length(x) - 2)
out <- num/den
return(out)
}  # output result of this function

## calculate centralization measures for weighted network
cent.wt <- function(sim) {
output <- matrix(0, 1, 3)
colnames(output) <- c("degree.wt", "eigen.wt", "between.wt")
# calculate degree centralization
output[1, 1] <- centralization(sim, sna::degree, normalize = T)
# calculate eigenvector centralization
output[1, 2] <- centralization(sim, sna::evcent, normalize = T)
# calculate betweenness centralization
output[1, 3] <- bw.cent(sim)
return(output)
}

Now let’s see examples of a couple of results using these functions focusing on the BR measure for now.

# BR binary net centralization
cent.bin(BRnet)
##         degree    between      eigen
## [1,] 0.2500917 0.05083638 0.04374436
# BR similarity centralization
cent.wt(chacoBR)
##      degree.wt   eigen.wt between.wt
## [1,] 0.1610866 0.03173933 4.6861e-10

# Assessing the potential impact of missing nodes

Now (finally) we’re ready to get into the meat of this workshop; assessing the impact of various forms of uncertainty on archaeological networks. The first example we will work with is the potential impact of missing nodes. In almost any archaeological network study, the networks we create are incomplete (i.e., we know that we are missing nodes for various reasons; site destruction, lack of survey coverage, looting, etc.). How might the fact that our networks are samples of a larger and typically unobtainable “complete network” influence our interpretations of network structure and node position? In this section, we are going to follow the lead of a study previously published in the journal Social Networks focused on this same issue (Costenbader and Valente 2003; see also Borgatti et al. 2006).

The overarching goal of the work by Costenbader and Valente (2003) was to take an existing network sampled from a (perhaps unknowable) population and determine how robust different structural and positional aspects of that network are to sub-sampling. This procedure, along with a careful consideration of the nature and quality of our sample, allows us to estimate the degree to which certain aspects of our sampled network do or do not approximate the same aspects of the original “complete” network from which it was drawn. Essentially, the idea is that if a particular metric is robust to sub-sampling in the data we have, we may have greater faith that a particular measure approximates the value in the complete network. Importantly, when we use the same underlying data to produce different networks as we have here (binary and weighted networks, networks of co-presence, BR similarity, and $$\chi^2$$ distance) we can further determine the degree to which different network definitions are vulnerable to different kinds of sampling issues. Due to this, we argue that it is good practice to explore a range of network definitions.

The basic approach advocated by Costenbader and Valente is to take a given network and then sub-sample it at intervals of 10% (producing some large number of replicates for each sampling fraction) and then calculate a set of network centrality and centralization metrics for every replicate of at each sampling fraction. We then calculate rank-order correlations (Spearman’s $$\rho$$) between the metrics for the original network sample and every sub-sample. This allows us to evaluate the average correlation and error at each sampling fraction.

In order to do this, we first provide a modifiable chunk of code that conducts this procedure for a single centrality (or centralization measure) and then plots the results as a boxplot. You may get some warning messages when you run this code if your random sub-samples end up creating a situation where some nodes connected to nothing (but the output will still be correct). For the sake of time in this workshop, we are only making 25 random replicates, but good practice suggests using many more (set the script running and go to lunch has always been my strategy).

#### Single measure example

nsim <- 25  # set number of replicates

# define function for testing betweenness centrality stability across
# sub-samples
resamp.test <- function(x) {
# the next line can be replaced with any centrality measure you'd like
bw.g <- sna::betweenness(x, gmode = "graph")
mat <- as.matrix(x)
dim.x <- dim(mat)[1]
out.mat <- matrix(NA, nsim, 9)  # create matrix for output and then name
colnames(out.mat) <- c("S90", "S80", "S70", "S60", "S50", "S40", "S30",
"S20", "S10")
# this double loop goes through each sampling fraction and each random
# replicate to cacluate centrality statistics and runs a Spearman's rho
# correlation between the resulting centrality values and the original
# sample
for (j in 1:9) {
for (i in 1:nsim) {
# this sampling procedure samples without replacement from a sequence from 1
# to the total number of nodes, the size of the sample being determined by
# 10-j/10. For example if j=1, 10-1/10 = 0.9 or a 90% sub-sample.
sub.samp <- sample(seq(1, dim.x), size = round(dim.x * ((10 - j)/10),
0), replace = F)
# calculate the betweenness statistic for the matrix reduced to only include
# the sub-sampled rows/columns
temp.stats <- sna::betweenness(mat[sub.samp, sub.samp], gmode = "graph")
# calcuate spearman's rho for replicate
out.mat[i, j] <- suppressWarnings(cor(temp.stats, bw.g[sub.samp],
method = "spearman"))
}
}
return(out.mat)
}  # return the result

# display results as boxplot by sampling fraction
boxplot(resamp.test(BRnet), ylim = c(0, 1), main = "BR - betweenness", xlab = "sampling fraction",
ylab = "Spearmans rho")

#### Combined example

The next chunk of code provides a single function that automates all of the steps described above using the centrality measures previously defined for both binary and weighted networks.

# The following function does the same calculation as above but is set up to
# work with the output of net.stats and net.stats.wt
nsim <- 25
samp.frac <- c("S90", "S80", "S70", "S60", "S50", "S40", "S30", "S20", "S10")

cv.resamp.bin <- function(x) {
# calculate all network stats for the original network
stats.g <- net.stats(x)
mat <- as.matrix(x)
dim.x <- dim(mat)[1]  # count number of rows (nodes)
# define empty matrices for output
dg.mat <- matrix(NA, nsim, 9)
ev.mat <- matrix(NA, nsim, 9)
bw.mat <- matrix(NA, nsim, 9)
# add column names based on sampling fraction
colnames(dg.mat) <- samp.frac
colnames(ev.mat) <- samp.frac
colnames(bw.mat) <- samp.frac

# this double loop goes through each sampling fraction and each random
# replicate to cacluate centrality statistics and runs a Spearman's rho
# correlation between the resulting centrality values and the original
# sample
for (j in 1:9) {
for (i in 1:nsim) {
sub.samp <- sample(seq(1, dim.x), size = round(dim.x * ((10 - j)/10),
0), replace = F)
temp.stats <- net.stats(mat[sub.samp, sub.samp])
dg.mat[i, j] <- suppressWarnings(cor(temp.stats[, 1], stats.g[sub.samp,
1], method = "spearman"))
ev.mat[i, j] <- suppressWarnings(cor(temp.stats[, 2], stats.g[sub.samp,
2], method = "spearman"))
bw.mat[i, j] <- suppressWarnings(cor(temp.stats[, 3], stats.g[sub.samp,
3], method = "spearman"))
}
}
out.list <- list()  # create list for output and populate it
out.list[[1]] <- dg.mat
out.list[[2]] <- ev.mat
out.list[[3]] <- bw.mat
return(out.list)
}  # return the resulting list

# And now for the version for weighted networks (similarity matrices)
cv.resamp.wt <- function(x) {
# calculate network stats for the original network
stats.g <- net.stats.wt(x)
mat <- as.matrix(x)
dim.x <- dim(mat)[1]
# create empty matrices for output
dg.mat <- matrix(NA, nsim, 9)
ev.mat <- matrix(NA, nsim, 9)
bw.mat <- matrix(NA, nsim, 9)
# add column names by sampling fraction
colnames(dg.mat) <- samp.frac
colnames(ev.mat) <- samp.frac
colnames(bw.mat) <- samp.frac

# this double loop goes through each sampling fraction and each random
# replicate to cacluate centrality statistics and runs a Spearman's rho
# correlation between the resulting centrality values and the original
# sample
for (j in 1:9) {
for (i in 1:nsim) {
sub.samp <- sample(seq(1, dim.x), size = round(dim.x * ((10 - j)/10),
0), replace = F)
temp.stats <- net.stats.wt(mat[sub.samp, sub.samp])
dg.mat[i, j] <- suppressWarnings(cor(temp.stats[, 1], stats.g[sub.samp,
1], method = "spearman"))
ev.mat[i, j] <- suppressWarnings(cor(temp.stats[, 2], stats.g[sub.samp,
2], method = "spearman"))
bw.mat[i, j] <- suppressWarnings(cor(temp.stats[, 3], stats.g[sub.samp,
3], method = "spearman"))
}
}
out.list <- list()  # create list for output
out.list[[1]] <- dg.mat
out.list[[2]] <- ev.mat
out.list[[3]] <- bw.mat
return(out.list)
}  # return the results

Now let’s run this for each of our binary and weighted networks defined above.

cop.rs <- cv.resamp.bin(Pnet)
BR.rs <- cv.resamp.bin(BRnet)
X.rs <- cv.resamp.bin(Xnet)

copw.rs <- cv.resamp.wt(chacoP)
BRw.rs <- cv.resamp.wt(chacoBR)
Xw.rs <- cv.resamp.wt(1 - chacoX01)

#### Plotting the results and assessing the impact

Once we’ve calculated the statistics for all of these sub-samples, it is then possible to plot the results to see how each network definition and each measure behaves in relation to sub-sampling. Let’s do that first for the three binary networks.

par(mfrow = c(3, 3))  # set up for 3 by 3 plotting
# plot boxplots by sampling fraction for each measure and each network type
boxplot(cop.rs[[1]], ylim = c(0, 1), main = "co-presence - degree", xlab = "sampling fraction",
ylab = "Spearmans rho")
boxplot(cop.rs[[2]], ylim = c(0, 1), main = "co-presence - eigenvector", xlab = "sampling fraction")
boxplot(cop.rs[[3]], ylim = c(0, 1), main = "co-presence - betweenness", xlab = "sampling fraction")
boxplot(BR.rs[[1]], ylim = c(0, 1), main = "BR - degree", xlab = "sampling fraction",
ylab = "Spearmans rho")
boxplot(BR.rs[[2]], ylim = c(0, 1), main = "BR - eigenvector", xlab = "sampling fraction")
boxplot(BR.rs[[3]], ylim = c(0, 1), main = "BR - betweenness", xlab = "sampling fraction")
boxplot(X.rs[[1]], ylim = c(0, 1), main = "Chi squared - degree", xlab = "sampling fraction",
ylab = "Spearmans rho")
boxplot(X.rs[[2]], ylim = c(0, 1), main = "Chi squared - eigenvector", xlab = "sampling fraction")
boxplot(X.rs[[3]], ylim = c(0, 1), main = "Chi squared - betweenness", xlab = "sampling fraction")


text:
Timezone: »

Hölder Bounds for Sensitivity Analysis in Causal Reasoning
Serge Assaad · Shuxi Zeng · Henry Pfister · Fan Li · Lawrence Carin

We examine interval estimation of the effect of a treatment T on an outcome Y given the existence of an unobserved confounder U. Using Hölder's inequality, we derive a set of bounds on the confounding bias |E[Y|T=t]-E[Y|do(T=t)]| based on the degree of unmeasured confounding (i.e., the strength of the connection U->T, and the strength of U->Y). These bounds are tight either when U⊥T or U⊥Y | T (when there is no unobserved confounding). We focus on a special case of this bound depending on the total variation distance between the distributions p(U) and p(U|T=t), as well as the maximum (over all possible values of U) deviation of the conditional expected outcome E[Y|U=u,T=t] from the average expected outcome E[Y|T=t]. We discuss possible calibration strategies for this bound to get interval estimates for treatment effects, and experimentally validate the bound using synthetic and semi-synthetic datasets.


text:
# Determine the equation of the plane in the form z=f(xy) that passes through the points ...

## Question:

Determine the equation of the plane in the form {eq}z=f(xy) {/eq} that passes through the points{eq}P(5,10,-9), \;Q(-1,-4,-10) {/eq} and {eq}R(4,-4,10) {/eq}.

## Planes:

To find the equation of a plane we need a normal vector of the plane and a point on the plane. We use the cross product of two vectors parallel to the plane to find the plane's normal vector.

Given three points {eq}P( 5, \, 10 , \, -9 ), \; Q( -1, \, -4 , \, -10 ), \; R( 4, \, -4 , \, 10 ) {/eq}, we find two vectors connecting two pairs...

Become a Study.com member to unlock this answer! Create your account


text:
## 過去の記録

#### 統計数学セミナー

13:00-14:10   数理科学研究科棟(駒場) 006号室
FMSP共催.参加をご希望される方は鎌谷 (阪大基礎工); kamatani at sigmath.es.osaka-u.ac.jpまでご連絡ください.
Ernst August Frhr. v. Hammerstein 氏 (Albert-Ludwigs-Universität Freiburg)
Laplace and Fourier based valuation methods in exponential Levy models (JAPANESE)
[ 講演概要 ]
A fundamental problem in mathematical finance is the explicit computation of expectations which arise as prices of derivatives. Closed formulas that can easily be evaluated are typically only available in models driven by a Brownian motion. If one considers more sophisticated jump-type Levy processes as drivers, the problem quickly becomes rather nontrivial and complicated. Starting with the paper of Carr and Madan (1999) and the PhD thesis of Raible (2000), Laplace and Fourier based methods have been used to derive option pricing formulas that can be evaluated very efficiently numerically. In this talk we review the initial idea of Raible (2000), show how it can be generalized and discuss under which precise mathematical assumptions the Laplace and Fourier approach work. We then give several examples of specific options and Levy models to which the general framework can be applied to. In the last part, we present some formulas for pricing options on the supremum and infimum of the asset price process that use the Wiener-Hopf factorization.

FMSP Lectures
http://faculty.ms.u-tokyo.ac.jp/~fmsp/jpn/conferences/fmsp.html
[ 参考URL ]
http://www.sigmath.es.osaka-u.ac.jp/~kamatani/statseminar/2012/13.html

http://faculty.ms.u-tokyo.ac.jp/~fmsp/jpn/conferences/fmsp.html

#### GCOEセミナー

16:00-17:00   数理科学研究科棟(駒場) 270号室
The obstacle problem for a shallow membrane-Justification and stability (ENGLISH)
[ 講演概要 ]
This lecture is twofold.
In the first part we recall the difference between the three-dimensional unilateral contact problem (the so-called Signorini problem) and its two-dimensional limit (the obstacle problem) in the case of an elastic shell as it was considered in [1] and [2].
In the second part we consider a simplified set of equations which describe the equilibrium equations of a shallow membrane (as justified in [3]) in contact with a plane obstacle and we study the stability of the contact zone with respect to small changes of the applied force, which amounts to studying the variation of the boundary of this contact zone. This kind of stability was first established in the scalar case for the Laplacian operator [4] then for the biharmonic operator [5]. The interest of the vectorial case considered here is due to the coupling effects between the in-plane and the transverse components of the displacement field in the framework of linearized Marguerre-von K´arm´an shell model.
This is a joint work with Alain L´eger, CNRS, Laboratoire de M´ecanique et d’Acoustique, 13402, Marseille, France.

### 2013年01月26日(土)

#### 調和解析駒場セミナー

13:30-18:00   数理科学研究科棟(駒場) 128号室
このセミナーは,月に1度程度,不定期に開催されます.

(東京大学) 13:30-15:00
On Triebel-Lizorkin spaces on Stratified Lie groups
(ENGLISH)
[ 講演概要 ]
We introduce the notion of Triebel-Lizorkin spaces
$\\dot{F}^{s}_{p,q}(G)$ on a stratified Lie group $G$
in terms of a Littlewood-Paley-type decomposition
with respect to a sub-Laplacian $\\mathscr{L}$ of $G$,
for $s \\in \\mathbb{R}$, $0 We show that the scale of these spaces is actually independent of the precise choice of the sub-Laplacian and the Littlewood-Paley-type decomposition. As we shall see, many properties of the classical Triebel-Lizorkin spaces on$\\mathbb{R}^{n}$, e.g., lifting property, embeddings and dual property, can be extended to the setting of stratified Lie groups without too much effort. We then study the boundedness of convolution operators on these spaces and finally, we obtain a Hormander type spectral multipliers theorem. 小野寺 有紹 (九州大学) 15:30-17:00 Profiles of solutions to an integral system related to the weighted Hardy-Littlewood-Sobolev inequality (JAPANESE) [ 講演概要 ] We study the Euler-Lagrange system for a variational problem associated with the weighted Hardy-Littlewood-Sobolev inequality of Stein and Weiss. We show that all the nonnegative solutions to the system are radially symmetric and have particular profiles around the origin and the infinity. This work extends previous results obtained by other authors to the general case. ### 2013年01月25日(金) #### 談話会・数理科学講演会 16:30-17:30 数理科学研究科棟(駒場) 002号室 旧記録は、上記セミナーURLにあります。 お茶&Coffee&お菓子: 16:00～16:30 (コモンルーム)。 中尾充宏 氏 (佐世保工業高等専門学校(校長)・九州大学名誉教授) 偏微分方程式の解に対する数値的検証の現状と動向 (JAPANESE) [ 講演概要 ] 本講演では、非線形偏微分方程式の解に対する 数値的検証法(精度保証付き数値計算法)について、現在までの研究の状況を概観し、今後の研究動向を考察する。特に、講演者らによってこれまでに得られた楕円型方程式を中心とする結果を総括し、最近の発展方程式に関する拡張についてもその現状と展望を述べる。さらに近年の力学系研究グループによる偏微分方程式に対する結果にも言及しその動向を探る。 ### 2013年01月24日(木) #### GCOEセミナー 14:00-15:00 数理科学研究科棟(駒場) 056号室 Leevan Ling 氏 (Hong Kong Baptist University) Global radial basis functions method and some adaptive techniques (ENGLISH) [ 講演概要 ] It is now commonly agreed that the global radial basis functions method is an attractive approach for approximating smooth functions. This superiority does not come free; one must find ways to circumvent the associated problem of ill-conditioning and the high computational cost for solving dense matrix systems. In this talk, we will overview different variants of adaptive methods for selecting proper trial subspaces so that the instability caused by inappropriately shaped parameters were minimized. #### GCOEセミナー 15:00-16:00 数理科学研究科棟(駒場) 056号室 Christian Clason 氏 (Graz University) Parameter identification problems with non-Gaussian noise (ENGLISH) [ 講演概要 ] For inverse problems subject to non-Gaussian (such as impulsive or uniform) noise, other data fitting terms than the standard L^2 norm are statistically appropriate and more robust. However, these formulations typically lead to non-differentiable problems which are challenging to solve numerically. This talk presents an approach that combines an iterative smoothing procedure with a semismooth Newton method, which can be applied to parameter identification problems for partial differential equations. The efficiency of this approach is illustrated for the inverse potential problem. ### 2013年01月23日(水) #### 諸分野のための数学研究会 10:30-11:30 数理科学研究科棟(駒場) 056号室 北海道大学のHPには、第1回(2005年6月22日)～第22回(2009年2月18日)までの情報が掲載されております。 Chun Liu 氏 (Pennsylvania State University) Ionic Fluids and Their Transport: From Kinetic Descriptions to Continuum Models (ENGLISH) [ 講演概要 ] The problems of distribution and transport of charged particles and ionic fluids are ubiquitous in our daily life and crucial in physical and biological applications. The multiscale and multiphysics nature of these problems provides major challenges and motivations. In this talk, I will present the derivation of the continuum models such as Poisson-Nerest-Planck (PNP) system from the kinetic formulations. Our focus is on the multiple species solutions and the corresponding boundary conditions for bounded containers. #### 作用素環セミナー 16:30-18:00 数理科学研究科棟(駒場) 118号室 David Evans 氏 (Cardiff University) Exotic subfactors and conformal field theories (ENGLISH) #### GCOEセミナー 16:30-18:00 数理科学研究科棟(駒場) 118号室 David Evans 氏 (Cardiff University) Exotic subfactors and conformal field theories (ENGLISH) [ 参考URL ] https://www.ms.u-tokyo.ac.jp/~yasuyuki/tokyo-seminar.htm #### GCOEセミナー 17:15-18:15 数理科学研究科棟(駒場) 370号室 Volker Schulz 氏 (Trier University) Shape and topology optimization in application (ENGLISH) [ 講演概要 ] Shape and topology optimization currently is of high interest for applications but also from a theoretical point of view. Recently, new developments in the shape calculus and in a related calculus for topology have enabled successful solutions of challenging optimization problems. This talk specifically reports on parameter free shape optimization in aerodynamics, thermoelastics and acoustics. Furthermore, novel results for the elastic topology optimization of the interior of wings are presented. We will try to give insight into the challenges in this field as well as the numerical solution approaches. ### 2013年01月22日(火) #### Lie群論・表現論セミナー 16:30-18:00 数理科学研究科棟(駒場) 126号室 Simon Goodwin 氏 (Birmingham University) Representation theory of finite W-algebras (ENGLISH) [ 講演概要 ] There has been a great deal of recent research interest in finite W-algebras motivated by important connection with primitive ideals of universal enveloping algebras and applications in mathematical physics. There have been significant breakthroughs in the rerpesentation theory of finite W-algebras due to the research of a variety of mathematicians. In this talk, we will give an overview of the representation theory of finite W-algebras focussing on W-algebras associated to classical Lie algebras (joint with J. Brown) and W-algebras associated to general linear Lie superalgebras (joint with J. Brown and J. Brundan). #### トポロジー火曜セミナー 16:30-18:00 数理科学研究科棟(駒場) 056号室 Tea: 16:00 - 16:30 コモンルーム Jarek Kedra 氏 (University of Aberdeen) On the autonomous metric of the area preserving diffeomorphism of the two dimensional disc. (ENGLISH) [ 講演概要 ] Let D be the open unit disc in the Euclidean plane and let G:=Diff(D, area) be the group of smooth compactly supported area-preserving diffeomorphisms of D. A diffeomorphism is called autonomous if it is the time one map of the flow of a time independent vector field. Every diffeomorphism in G is a composition of a number of autonomous diffeomorphisms. The least amount of such diffeomorphisms defines a norm on G. In the talk I will investigate geometric properties of such a norm. In particular I will construct a bi-Lipschitz embedding of the free abelian group of arbitrary rank to G. I will also show that the space of homogeneous quasi-morphisms vanishing on all autonomous diffeomorphisms in G is infinite dimensional. This is a joint work with Michael Brandenbursky. ### 2013年01月21日(月) #### トポロジー火曜セミナー 16:30-17:30 数理科学研究科棟(駒場) 002号室 日時と場所がいつもとは違います。 加藤 直樹 氏 (東京大学大学院数理科学研究科) ベキ零リー環を横断構造に持つリー葉層構造について (JAPANESE) [ 講演概要 ] リー$\\mathfrak{g}$-葉層構造が与えられたとき,それに付随した構造リー環と呼 ばれる$\\mathfrak{g}$の部分リー環$\\mathfrak{h}$が定まる.本講演ではその逆, すなわち,リー環$\\mathfrak{g}$とその部分リー環$\\mathfrak{h}$が与えられた とき構造リー環が$\\mathfrak{h}$となるリー$\\mathfrak{g}$-葉層構造が存在す るかという問題についてて,$\\mathfrak{g}$がベキ零リー環の場合に説明をする. #### トポロジー火曜セミナー 17:30-18:30 数理科学研究科棟(駒場) 002号室 日時と場所がいつもとは違います。 石田 智彦 氏 (東京大学大学院数理科学研究科) Quasi-morphisms on the group of area-preserving diffeomorphisms of the 2-disk (JAPANESE) [ 講演概要 ] Gambaudo and Ghys constructed linearly independent countably many quasi- morphisms on the group of area-preserving diffeomorphisms of the 2-disk from quasi-morphisms on braid groups. In this talk, we will explain that their construction is injective as a homomorphism between vector spaces of quasi-morphisms. If time permits, we introduce an application by Brandenbursky and K\\c{e} dra. #### 複素解析幾何セミナー 10:30-12:00 数理科学研究科棟(駒場) 126号室 雨宮 卓史 氏 (東大数理) コンパクト複素多様体への有理型写像の値分布について (JAPANESE) [ 講演概要 ] In a late paper of J. Noguchi and J. Winkelmann they showed the condition of being Kähler or non-Kähler of the image space to make a difference in the value distribution theory of meromorphic mappings into compact complex manifolds. In the present talk, we will discuss the order of meromorphic mappings to a Hopf surface which is more general than dealt with by Noguchi-Winkelmann, and an Inoue surface (they are non-Kähler surfaces). For a general Hopf surface$S$, we prove that there exists a differentiably non-degenerate holomorphic mapping$f:\mathbf{C}^2 \to S$whose order satisfies$\rho_{f}\leq 1$. For an Inoue surface$S'$, we prove that every non-constant meromorphic mapping$f:\mathbf{C}^n \to S'$is holomorphic and its order satisfies$\rho_{f}\geq 2$. ### 2013年01月16日(水) #### 幾何コロキウム 10:30-12:00 数理科学研究科棟(駒場) 122号室 開始時間と開催場所などは変更されることがあるので, セミナーごとにご確認ください. 田中祐二 氏 (京都大学) A construction of Spin(7)-instantons (JAPANESE) [ 講演概要 ] Spin(7)-instantons are elliptic gauge fields on 8-dimensional Spin(7)-manifolds, which minimize the Yang-Mills action. Analytic properties of Spin(7)-instantons have been studied by Gang Tian and others, but little was known about the existence of examples of Spin(7)-instantons other than an Oxford Ph.D thesis by Christopher Lewis in 1998. There are two known constructions of compact Spin(7)-manifolds both obtained by Dominic Joyce. The first one begins with a torus orbifold of a special kind with non-isolated singularities. The Spin(7)-manifold is obtained by resolving the singularities with the aid of algebraic geometry techniques. The second one begins with a Calabi-Yau four-orbifold with isolated singular points of a special kind and an anti-holomorphic involution fixing only the singular points. The Spin(7)-manifold is obtained by gluing ALE Spin(7)-manifolds with anti-holomorphic involutions fixing only the origins to each singular point. Christopher Lewis studied the problem of constructing Spin(7)-instantons on Spin(7)-manifolds coming from Joyce's first construction. This talk describes a general construction of Spin(7)-instantons on examples of compact Spin(7)-manifolds coming from Joyce's second construction. Starting with certain Hermitian-Einstein connections on the Calabi-Yau four-orbifold and on ALE Spin(7)-manifolds, we glue them together simultaneously with the underlying pieces to make a Spin(7)-instanton on the compact Spin(7)-manifold by Joyce. #### 代数学コロキウム 18:00-19:00 数理科学研究科棟(駒場) 002号室 大久保俊 氏 (東京大学数理科学研究科) On differential modules associated to de Rham representations in the imperfect residue field case (ENGLISH) [ 講演概要 ] Let K be a CDVF of mixed characteristic (0,p) and G the absolute Galois group of K. When the residue field of K is perfect, Laurent Berger constructed a p-adic differential equation N_dR(V) for any de Rham representation V of G. In this talk, we will generalize his construction when the residue field of K is not perfect. We also explain some ramification properties of our N_dR, which are due to Adriano Marmora in the perfect residue field case. (本講演は「東京北京パリ数論幾何セミナー」として、インターネットによる東大数理, Morningside Center of MathematicsとIHESの双方向同時中継で行います.) #### 講演会 10:00-11:00 数理科学研究科棟(駒場) 123号室 R\'emi Boutonnet 氏 (ENS Lyon)$W^*$-superrigidity of mixing Gaussian actions of rigid groups (ENGLISH) #### 講演会 11:30-12:30 数理科学研究科棟(駒場) 123号室 Tim de Laat 氏 (University of Copenhagen) The Approximation Property for Lie groups (ENGLISH) #### 講演会 14:40-15:40 数理科学研究科棟(駒場) 118号室 Arnaud Brothier 氏 (KU Leuven) Unique Cartan decomposition for II$_1$factors arising from cross section equivalence relations (ENGLISH) #### 講演会 15:55-16:55 数理科学研究科棟(駒場) 118号室 Michael Hartglass 氏 (UC Berkeley) Rigid$C^*$tensor categories of bimodules over interpolated free group factors (ENGLISH) #### 講演会 17:10-18:10 数理科学研究科棟(駒場) 118号室 James Tener 氏 (UC Berkeley) Manifestly unitary conformal field theory (ENGLISH) #### GCOEセミナー 14:40-15:40 数理科学研究科棟(駒場) 118号室 このセミナーはMiniworkshop on Operator Algebras IIの一環として行われます Arnaud Brothier 氏 (KU Leuven) Unique Cartan decomposition for II$_1\$ factors arising from cross section equivalence relations (ENGLISH)
[ 参考URL ]
https://www.ms.u-tokyo.ac.jp/~yasuyuki/mini2013-2.htm

#### GCOEセミナー

15:55-16:55   数理科学研究科棟(駒場) 118号室
このセミナーはMiniworkshop on Operator Algebras IIの一環として行われます
Michael Hartglass 氏 (UC Berkeley)
TBA (ENGLISH)
[ 参考URL ]
https://www.ms.u-tokyo.ac.jp/~yasuyuki/mini2013-2.htm

#### GCOEセミナー

16:00-17:00   数理科学研究科棟(駒場) 370号室
Oleg Emanouilov 氏 (Colorado State University)
3-D Calderon's Problem with partial Dirichlet-to Neumann map (ENGLISH)
[ 講演概要 ]
We present new results for the uniqueness of recovery of a potential in three dimensional Calderon's problem with partial Dirichlet-to-Neumann map.
The proof is based on complex geometric optics solutions and the Radon transform.


text: Abstractan overview on the triangular arbitrage, followed by an example presents the operation and its practical nuances, including order issued. its most relevant aspects are depicted in pseudo code. we move one to obtain the profit conditions, for the given example. conditions are generalized, a framework for comparing multiple opportunities is presented, and a trading algorithm is proposed based on the obtained conditions.


text:
# Source code for spynnaker8.models.connectors.distance_dependent_prob_connector

# Copyright (c) 2017-2021 The University of Manchester
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
from spynnaker.pyNN.models.neural_projections.connectors import (
DistanceDependentProbabilityConnector as
_BaseClass)
from spynnaker.pyNN.utilities.utility_calls import moved_in_v6

[docs]class DistanceDependentProbabilityConnector(_BaseClass):
""" Make connections using a distribution which varies with distance.

.. deprecated:: 6.0
Use
:py:class:spynnaker.pyNN.models.neural_projections.connectors.DistanceDependentProbabilityConnector
"""
__slots__ = []

def __init__(
self, d_expression, allow_self_connections=True, safe=True,
verbose=False, n_connections=None, rng=None, callback=None):
"""
:param str d_expression:
the right-hand side of a valid python expression for
probability, involving d, e.g. "exp(-abs(d))", or "d<3",
that can be parsed by :py:func:eval, that computes the distance
dependent distribution
:param bool allow_self_connections:
if the connector is used to connect a
Population to itself, this flag determines whether a neuron is
allowed to connect to itself, or only to other neurons in the
Population.
:param bool safe: if True, check that weights and delays have valid
values. If False, this check is skipped.
:param bool verbose:
Whether to output extra information about the connectivity to a
CSV file
:param int n_connections:
The number of efferent synaptic connections per neuron.
:param ~pyNN.random.NumpyRNG rng: random number generator
:param callable callback:
if given, a callable that display a progress bar on the terminal.

.. note::
Not supported by sPyNNaker.
"""
# pylint: disable=too-many-arguments
moved_in_v6("spynnaker8.models.connectors"
".DistanceDependentProbabilityConnector",
"spynnaker.pyNN.models.neural_projections.connectors"
".DistanceDependentProbabilityConnector")
_BaseClass.__init__(
self, d_expression=d_expression,
allow_self_connections=allow_self_connections,
safe=safe, verbose=verbose, n_connections=n_connections)


text:
# Show that the statement " If $x$ is real number such that $x^3+4x=0$, then $x$ is 0." is true using contradiction method.
" If $x$ is real number such that $x^3+4x=0$, then $x$ is 0."
Let $p :$ If $x$ is real number such that $x^3+4x=0. and$q:x=0$. To show that the above statement is true using contradiction method, Let us assume that$p$is true. and$q$is not true.$i.e.,x$is a real number such that$x^3+4x=0$and$x\neq 0\Rightarrow\:\:x(x^2+4)=0\Rightarrow\:$either$x=0$or$x^2+4=0$But our assumption is$x\neq 0$and also$x^2+4$cannot be$0$since$x^2 \geq 0\:\: \forall x \in R\thereforex^3+4x$cannot be equal to zero which is contradiction.$\thereforex$should be equal to$0$.$\therefore\$ The given statement is true.


text:
# Installing¶

The package works both on Linux and on Windows. You can find the last stable version at PyPI and install it with the usual Python installing mechanism:

python setup.py install

If you have pip, you can go directly (use -U for update or reinstall):

pip install audiolazy

pip install -U .

To install from a path that has the setup.py file and the package data uncompressed previously.

For the bleeding-edge version, you can install directly from the github repository (requires git for cloning):

pip install -U git+git://github.com/danilobellini/audiolazy.git

For older versions, you can install from the PyPI link or directly from the github repository, based on the repository tags. For example, to install the version 0.04 (requires git for cloning):

pip install -U git+git://github.com/danilobellini/audiolazy.git@v0.04

The package doesn’t have any strong dependency for its core besides the Python itself (versions 2.7, 3.2 or newer) as well as its standard library, but you might need:

• PyAudio: needed for playing and recording audio (AudioIO class);
• NumPy: needed for doing some maths, such as finding the LSFs from a filter or roots from a polynomial;
• MatPlotLib: needed for all default plotting, like in LinearFilter.plot method and several examples;
• SciPy (testing and examples only): used as an oracle for LTI filter testing and for the Butterworth filter example;
• Sympy (testing only): used for testing linear filters with time-varying matrices of symbolic coeffs where the Stream samples are these matrices;
• tox for testing all at once, or pytest, pytest-cov and pytest-timeout for testing in a single environment (testing only): runs test suite and shows code coverage status;
• wxPython (example only): used by one example with FM synthesis in an interactive GUI;
• Tkinter (example only): needed for the pitch follower based on the zero-crossing rate example GUI;
• Music21 (example only): there’s one example that gets the Bach chorals from that package corpora for synthesizing and playing;
• Sphinx (documentation only): it can create the software documentation in several different file formats.

Beside examples and tests, only the filter plotting with plot and zplot methods needs MatPlotLib. Also, the routines that needs NumPy up to now are:

• Root finding with zeros and poles properties (filter classes) or with roots property (Poly class);
• Some Linear Predictive Coding (lpc) strategies: nautocor, autocor and covar;
• Line Spectral Frequencies lsf and lsf_stable functions.


text:
Project Togo began in February 1998 when a Toronto-based ringing operation was dismantled and one of its participants informed the police that he was previously employed by a Montreal businessman who was also active in the resale of stolen vehicles. This initial tip was corroborated soon after by a thief who had been arrested while driving a stolen vehicle. By December 1998, the Togo investigation was under way. It spanned into February 1999 and 20 cars that were destined for France, Ghana, and local buyers in southern Quebec were retrieved. NETWORK 1-mode matrix 33 x 33 person by person. Undirected ties. Ties are communication exchanges between criminals. Data comes from police wiretapping.

covert_47

igraph object

Source

Available from Manchester (https://sites.google.com/site/ucinetsoftware/datasets/covert-networks), reconstructed from Morselli’s book, Inside Criminal Networks http://www.springer.com/social+sciences/criminology/book/978-0-387-09525-7 book page 187, Appendix

References

Morselli, C., 2009. Inside criminal networks. New York: Springer.


text:
Detect if all required glyphs are available

In XeTeX, a message is printed to the log file if a glyph is not available with the current font. The input

\documentclass{scrartcl}

\begin{document}
α
\end{document}

leads to the following entry in the log:

Missing character: There is no α in font cmr10!

Is there a way to check for a given input if all glyphs required to render it are available, without relying on the log file? Apparently, there is no Command to flush log file, and (in a scenario where I pipe into a running XeLaTeX process) I'd like to know before XeTeX is terminated.

• No, it's not possible; typesetting a character happens at a very deep stage in processing, when there's no user's control any more. I guess something can be done with LuaTeX, but not with XeTeX, that only provides \iffontchar\font<charcode> to check if the <charcode> corresponds to an existing glyph in the current font. – egreg Mar 15 '14 at 12:28
• @egreg: Thanks. Can XeTeX be told to treat this condition as an error, so that the process would abort when -halt-on-error is given on the command line? – krlmlr Mar 15 '14 at 13:15
• No, this is not possible; but you'd have to look at the log file anyway, wouldn't you? – egreg Mar 15 '14 at 14:17
• @egreg: It seems that I can only reliably parse the log file after XeTeX has terminated. Remember that I intend to keep XeTeX open and feed it via its stdin. See github.com/yihui/tikzDevice/issues/5 for some more detail, I'd be glad for any input. – krlmlr Mar 15 '14 at 14:37
• I think you're doomed, sorry. :( – egreg Mar 15 '14 at 14:41

Longer post

This task was a big issue for me many years ago when I was typesetting and generating preview books of fonts where I needed to know all kind information about the glyphs before actual typesetting. I divided the answer into three small parts each with compilable TeX file.

Part 1: xelatex and its tool

There is a way in xelatex, it is mentioned indirectly on page 8 in the reference manual. We can use the \XeTeXcharglyph command and check if the output is 0 (or .notdef if checking glyph name). Once any glyph is not defined the zero value is assigned to the character slot.

After that step we can use \ifnum...\else...\fi statement to decide what to typeset, if anything at all. In the example I am testing four letters, ? A \ a, two are defined, the remaining two are not. I have downloaded and installed RodgauApesInitials from Manfred Klein's collection, http://moorstation.org/typoasis/designers/klein04/deco/rodgau_apes.htm

We run xelatex mal-kanji1.tex and terminal says:

The glyph slot of 003F(hex) is: 0 (.notdef)
The glyph slot of 0041(hex) is: 3 (A)
The glyph slot of 005C(hex) is: 0 (.notdef)
The glyph slot of 0061(hex) is: 29 (a)

This is the code:

% run: xelatex mal-kanji1.tex
\documentclass[a4paper]{article}
\pagestyle{empty}
\usepackage{xltxtra}
\usepackage{fontspec}
\usepackage{pgffor}
\begin{document}
%\font\malfont=RodGauApes.ttf
\font\klein="RodGauApes Initials"
\font\cmr=cmr10%
\message{^^J}% One \n to the terminal...
% U+0041, A, dec 65
\def\checkme#1{%
\klein\message{The glyph slot of #1(hex) is: \the\XeTeXcharglyph"#1 \space
(\XeTeXglyphname\klein\XeTeXcharglyph"#1)^^J}%
%! Cannot use \XeTeXcharglyph with cmr10; not a native platform font.
\ifnum\XeTeXcharglyph"#1=0%
\cmr Undefined glyph!%
\else
\klein\char"#1%
\fi% End of \ifnum...
}% End of \checkme...
\foreach \malglyph in {003F,0041,005C,0061}% ? A \ a
{\checkme\malglyph\cmr\par%
}% End of \foreach...
\end{document}

I needed and I used this strategy when I typesetted a preview book of kanji. I was using many different font faces and I needed to know exactly how many glyphs were going to be typeset. If a glyph or more were missing it was up to me to decide if to typeset that rectangle, or not. I should say that it is normal that a font is missing some CJKV characters, there are so many of them. I enclose a snapshot from the book, it is a preview of a single Japanese character/kanji.

The big disadvantage is that we can use this strategy only with the xelatex format and we can check only newer fonts (TTF, OTF), it is not working with PFB. That's a serious limitation in work. Let me show you other approach.

Part 2: testfont.tex

One of tools in the TeX distribution is testfont.tex file. If we run, e.g.:

pdftex testfont.tex

We are asked about the name of the font, let's use dmjhira and hit Enter. There is a list of options when we use \help and hit Enter. The common option is to write \table\bye followed by Enter. We are getting the testfont.pdf file and this is a preview of it.

A problem is that we cannot customize it and we cannot enter TTF and OTF fonts without installing and setting up the fonts. It is worth checking another approach.

Part 3: Measuring a TeX box (Hiragana and Katakana as test cases)

The core of TeX is in measuring boxes. It is the approach we are going to use and test. We will test a single glyph, but we could easily measure almost anything. After defining a new box (\newbox), we virtually fulfill the box (\setbox) and we can measure width, height and depth of the box (\wd, \ht and \dp).

In the example below we typeset Hiragana, a Japanese syllabary, but we may spot several defects. There is an opening space, lines don't break ideally, spaces are shrinking and stretching, we have no idea how many glyphs were typesetted. But it is a starting point for further explorations.

Worth mentioning is that we can run any major latex engine and if we use e.g. fontspec package, we can test glyphs in the TTF, OTF and some dfont files when running lualatex or xelatex. I enclose the TeX code and a preview of the result.

% run: *latex mal-kanji2.tex
\documentclass[a4paper]{article}
\rightskip=6cm% And TeXie, :-), please narrow the text width somehow...
\begin{document}
\font\hira=dmjhira at 2ex% Setting up a new font face (Hiragana)...
\newcount\malcounter% Setting up a new counter...
\malcounter=-1% The initialization of the counter
\loop% The core of this example...
\advance\malcounter by 1% Move on to the next glyph...
{\hira\char\malcounter} % Show me the glyph, add space...
%\discretionary{}{}{}% Allow a page break after the glyph... (the first alternative)
\allowbreak% (the second alternative)
\ifnum\malcounter<255\repeat% Run all 256 characters in the font...
\end{document}

We typeset a glyph only when glyph width is a nonzero value. There are some exceptions in ornament fonts, but it is considered bad fontography. We add a counter (\newcount) taking care of a number of typesetted glyphs, we limit spaces to fixed dimension (\fontdimen 3, 4 and 7), we allow line breaks, simply said we are having full control over the output now.

I enclose an example and a preview of the PDF file with Katakana, a Japanase syllabary.

% run: *latex mal-kanji3.tex
\documentclass[a4paper]{article}
\pagestyle{empty}
\parindent=0pt
\begin{document}
\newbox\emptybox
\setbox\emptybox=\hbox{}
\newbox\malbox
\font\kata=dmjkata% A new font face to be used (Katakana)...
\newcount\counter
\counter=-1% The counter of glyphs...
\fontdimen3\font=0pt \fontdimen4\font=0pt \fontdimen7\font=0pt
% Eliminate the stretch in spaces, or, we could use \makebox...
\loop% Process all the glyphs...
\advance\counter by 1% Go to the next glyph...
\setbox\malbox=\hbox{\kata\char\counter}% Measure the glyph...
\ifnum\wd\malbox=\wd\emptybox\relax\else% Is width 0pt? Height is not tested (\ht, \dp)...
\texttt{\ifnum\counter<10 0\fi% Add zero in front of number <10...
\the\counter.\copy\malbox\ }% Show me the glyph...
\discretionary{}{}{}%Allow the line break, or, \allowbreak...
\fi% End of \ifnum condition...
\ifnum\counter<255\repeat% Show me all the glyphs in range 0-255...
\end{document}

Conclusions

There are some limitations of these methods. We cannot test if two glyphs are not just rotated, flipped or scaled, we are also not listing glyphs outside 0-255 region. In case you need something like that, use FontForge, it can be done there as it has own scripting language and a support for Python scripts.

In case you are wondering what was that splash in the very first example in the last letter under the letter: Oh, yes, that's an ape from the great Manfred Klein legendary font creator himself! I enclose a screenshot from FontForge and a preview of his font with 26+26 letters in it.

As an extension of Malipivo's excellent answer, I wanted to mention that you can also use \iffontchar as shown in the following example. Note that LuaLaTeX gives different results from XeLaTeX here -- it looks like XeLaTeX doesn't like glyphs in the very upper ranges of Unicode here (above U+F0000 perhaps).

The code should be fairly evident, but the idea is that you call

\fontrange {xits-math.otf} {1} {"FFFFF} { \texttt{\small #2}:\quad #3\par}

where the first argument is the fontspec font, the second is the starting glyph point, the second is the ending glyph point, and the fourth argument defines how the result is typeset (#2 is the hex code of the glyph and #3 is the character corresponding to the glyph itself).

\documentclass[a4paper]{article}
\usepackage{fontspec,multicol}
\usepackage[margin=2cm]{geometry}
\begin{document}

\ExplSyntaxOn
\DeclareDocumentCommand \fontrange { m m m +m }
{
\group_begin:
\cs_set:Npn \__typeset:nnn ##1 ##2 ##3 {#4}
\fontspec{#1}
\int_step_inline:nnnn {#2} {1} {#3}
{
\iffontchar\font ##1\relax
\__typeset:nnn {##1} { \int_to_Hex:n {##1} } { \char ##1 \relax }
\fi
}
\group_end:
}
\ExplSyntaxOff

\begin{multicols}{5}
\fontrange{xits-math.otf}{1}{"FFFFF}

• I'm not sure this fits into fontspec; it seems more like fonttable. – egreg May 27 '14 at 6:35
• @egreg I agree to some extent — but all I'd be providing in fontspec would be something like that \fontrange command and it'd be up to whomever to do something useful with it. Making it flexible would probably be difficult, though; so perhaps not. – Will Robertson May 28 '14 at 2:38


text:
# How can I solve this probability question?

An experiment consisting of testing calculators one after the other with no replacement until either 2 defective calculators are found or 4 are tested. Find the probability of:

a)the Event E consisting of all outcomes where 1 defective calculator is tested.

I tried analyzing every possible outcome, by looking at every possible combination of defective calculators and successful ones, but I couldn't come up with the right solution. How would I solve the problem?

-
You haven't actually stated a problem. –  Gerry Myerson Dec 9 '12 at 4:55
@GerryMyerson Oops, sorry. I have edited the question accordingly. –  asf Dec 9 '12 at 4:58
You have to know the probability of a calculator being defective. –  Gerry Myerson Dec 9 '12 at 5:00

Let 0 means perfect, 1 means faulty.

The possible outcomes are

$0000,0001$;$0010,0011$ ; $011$; $0100,0101$; $11$;$101$;$1001,1000$;

Among them $0001,0010,0100,1000$ have exactly 1 defective calculator.

So, the required probability will be $\frac4{11}$

If we continue up to $4$ tests the number possible cases are $2^4=16$

Exactly 3 $1$s can occur in $\binom 43=4$ ways.

Exactly 4 $1$s can occur in $\binom 44=1$ way.

So, the number of available cases are $=16-(4+1)=11$

Among $4$ tests, exactly 1 defective calculator can occur in $\binom 41=4$ ways.

-
Thank You! Just one question: is there a way to do this without listing out all possibilities? A generalized formula for the problem maybe? Because if say the the problem was scaled significantly, listing out all possibilities wouldn't be an option. –  asf Dec 9 '12 at 5:32
@asf, please let me know your idea with the edited answer. –  lab bhattacharjee Dec 9 '12 at 5:51
It's perfect, thank you again! –  asf Dec 9 '12 at 5:54

Write D for defective, O for OK. The possible outcomes are DD, DOD, DOOD, ODD, ODOD, OODD, OODO, OOOD, OOOO. Presumably what you want to do is figure out the probability of each of these 9 possible outcomes, and add up the ones that belong to your event E.

-
The answer books says otherwise. It doesn't give a solution but the answer says: 4/11. You also missed several possible outcomes. E.g. DDO –  asf Dec 9 '12 at 5:01
Or am I stupidly misunderstanding your answer? –  asf Dec 9 '12 at 5:07
You stop when you get 2 defectives, right? So, DDO can't happen --- you stop after the DD. But I did leave out DOOO and ODOO. –  Gerry Myerson Dec 9 '12 at 5:28
Ahh I see. The book says there are 11 total outcomes, not 9. Is it wrong? My teacher confirmed the book is correct. –  asf Dec 9 '12 at 5:30