text:
# SYDE 556/750: Simulating Neurobiological Systems¶

## Temporal Representation¶

• In the previous lectures, we covered representation of a vector space by a population of neurons
• We converted a vector $x$ into neural activity values $a_i$, and then converted those back to the original value $\hat{x}$.
• What happens if $x$ changes over time?
• I.e. if we have $x(t)$, how do we get $a_i(t)$?
• Seems pretty easy:
• $a=G[\alpha e \cdot x + J^{bias}]$, so...
• $a(t)=G[\alpha e \cdot x(t) + J^{bias}]$
• where $G[J(t)]= {1 \over {\tau_{ref}-\tau_{RC}ln(1-{1 \over J(t)})}}$
In [1]:
%pylab inline

//anaconda/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

Populating the interactive namespace from numpy and matplotlib

In [2]:
import nengo
import numpy

n = nengo.neurons.LIFRate()

x = numpy.zeros(100)
for i in range(10):
x[i*10:(i+1)*10]=numpy.random.uniform(-1,1)

figure()
title('$x$ value changing over time')
plot(x)
ylabel('$x$')
ylim(-1,1)
xlabel('Time (ms)')

figure()
title('Firing rate changing over time')
plot(n.rates(x, gain=1.5, bias=1))
ylabel('Firing rate (Hz)')
xlabel('Time (ms)')
show()

• Unfortunately, that formula for $G[J(t)]$ is an average firing rate over a long period of time.
• What does it mean to be firing at 10Hz for 10ms?
• We need to model what's actually happening to that neuron over time

## Leaky Integrate-and-Fire neurons¶

• Pretty much the standard computationally efficient neuron model
• Simple
• Produces spikes
• Limiting case of other more complex neuron models
• All parameters map onto known properties of real neurons

• Bilipid cell membrane acts as a capacitor

• Flow of ions through the cell membrane ion channels is a resistor

• This is the "leak" current
• Resistance and capacitance together give $\tau_{RC}$
• $\tau_{RC}=R C$
• When voltage passes some threshold, neuron emits a "spike"
• Very fast, sharp response that's pretty much the same each time
• Spike at time $t_n$ modelled as $\delta(t-t_n)$
• Neuron recovers over a short period of time
• Refractory time $\tau_{ref}$
• During this time the voltage is reset to some resting level

• As an RC circuit with a nonlinearity:

• Typical values

• $\tau_{RC}$: 0.02 seconds
• $\tau_{ref}$: 0.002 seconds
• Reset voltage (a.k.a. resting potential): -70mV
• Firing threshold (a.k.a. threshold potential): -55mV
• NOTE: we normalize this so that reset is at 0 and firing is at 1
• This does not affect anything about the behaviour of the model
• For lots more info, see Unit 2 in this online course

### The LIF RC circuit¶

• Let's ignore the spiking elements for now, so here's the circuit

### Useful physical laws¶

• Kirchhoff's Current Law

• Currents at a point have to add up to zeros
• Current coming into the neuron has two paths
• $J_{M} = J_R+J_C$
• Ohm's Law

• Formula for current through a resistor
• $J_R = {V \over R}$
• Capacitor

• $C = {Q \over V}$
• $Q = C V$
• Change in these quantities over time gives
• ${{dQ} \over {dt}} = C {{dV} \over {dt}}$
• Rate of change of charge is current
• $J_C = C {{dV} \over {dt}}$
• Applying these laws to RC circuits:
• $J_{M} = {V \over R} + C {{dV} \over {dt}}$
• Then solve for change in voltage with input

\begin{align} {{dV} \over {dt}} &= {1 \over C}(J_{M} - {V \over R}) \\ &= {1 \over {RC}}[R J_{M} - {V}] \\ &= {1 \over {\tau_{RC}}} [R J_{M} - {V}] \end{align}

### LIF RC Circuit Behaviour¶

• If we solve the general ODE, we get a convolution:

\begin{align} V(t) = \frac{R}{\tau_{RC}} \int_0^t e^{-(t-t')/\tau_{RC}} J_M(t')\; dt' \end{align}

• I.e. the voltage right now (at $t$) depends on all past current input, $J_M(t')$, where each input is weighted by a function that exponentially decays as it gets further away (in the past) from the current time (a kind of memory).
• If $J_M$ is constant:

\begin{align} V(t) &= \frac{R J_M}{\tau^{RC}} e^{-t/\tau_{RC}}\int_0^t e^{t'/\tau_{RC}} \; dt' \\ &= J_M R (1-e^{-t/\tau_{RC}}) \end{align}

• How long does it take to get to threshold?
• How long between spikes?
• Firing rate?
• How long does it take to get to threshold?

• $V_{th} = J_MR(1-e^{-t_{th}/\tau_{RC}})$
• $t_{th} = -\tau_{RC} ln (1-{V_{th} \over {J_MR}})$
• How long between spikes?

• $t_{th} + \tau_{ref}$
• Firing rate

• $a = {1 \over {t_{th} + \tau_{ref}}}$
• $a = {1 \over {\tau_{ref}-\tau_{RC}ln(1-{V_{th} \over {J_MR}})}}$
• Simplify by assuming $V_{th}=1$ and $R=1$

• $a = {1 \over {\tau_{ref}-\tau_{RC}ln(1-{1 \over J_M})}}$
• This is the rate approximation of the neuron
• But we can't assume a constant $J_M$ in a real brain!

• It will definitely change over time (since $x(t)$ also changes over time)
• So we need to do a computer simulation of the differential equation explicitly
• ${{dV} \over {dt}} = {1 \over {\tau_{RC}}} [R J_{input}(t) - {V(t)}]$

## Spikes¶

• When simulating the LIF neuron, we need to monitor threshold, as that when the 'spike' occurs
• A LIF neuron treats all spikes the same, and just 'pastes' them in
• As opposed to simulating the nonlinear differential equations that describe the actual voltage changes
In [59]:
import nengo
from nengo.utils.matplotlib import rasterplot

model = nengo.Network(label='Spiking Neurons')

with model:
stimulus = nengo.Node(0)
ens = nengo.Ensemble(1, dimensions=1,
encoders = [[1]],
intercepts = [-.5],
max_rates= [100])

nengo.Connection(stimulus, ens)

spikes_p = nengo.Probe(ens.neurons, 'spikes')
voltage_p = nengo.Probe(ens.neurons, 'voltage')
stim_p = nengo.Probe(stimulus)

sim = nengo.Simulator(model)
sim.run(.5)

t = sim.trange()
figure(figsize=(12, 8))
ax = gca()
ax.plot(t, sim.data[voltage_p],'g')
rasterplot(t, sim.data[spikes_p], ax=ax)
ylim((0,2))
ylabel("Voltage")
xlabel("Time");

Building finished in 0:00:01.
Simulating finished in 0:00:01.

In [60]:
import numpy as np
from nengo.utils.functions import piecewise

x = dict(zip(np.linspace(0,.1,10), np.random.uniform(-1,1,10)))

with model:
stimulus.output = piecewise(x)

sim = nengo.Simulator(model)
sim.run(.1)

t = sim.trange()

figure()
title('$x$ value changing over time')
plot(t, sim.data[stim_p])
ylabel('$x$')
ylim(-1,1)
xlabel('Time (s)')

figure()
title('Firing rate changing over time')
plot(ens.neuron_type.rates(sim.data[stim_p], gain=sim.data[ens].gain, bias=sim.data[ens].bias))
ylabel('Firing rate (Hz)')
xlabel('Time (ms)')

figure(figsize=(6, 3))
ax = gca()
rasterplot(t, sim.data[spikes_p], ax=ax)
ax.plot(t, sim.data[voltage_p],'g')#, drawstyle='steps-post') #if want to see exact timing
ylim((0,2))
title('Voltage and spikes over time')
ylabel("Voltage")
xlabel("Time (s)");

Building finished in 0:00:01.
Simulating finished in 0:00:01.

• Let's try a sine wave input
In [61]:
#have to run previous cells

with model:
stimulus.output = lambda t: np.sin(10*t)

sim = nengo.Simulator(model)
sim.run(.6)

nengo.utils.matplotlib.plot_tuning_curves(ens, sim)

t = sim.trange()
figure(figsize=(10,4))
plot(t, sim.data[stim_p], label='x')
ax = gca()
ax.plot(t, sim.data[voltage_p],'g', label='v')
ylim((-1,2))
ylabel('$x$')
xlabel('time (s)')
legend(loc='lower left');

rasterplot(t, sim.data[spikes_p], ax=ax.twinx(), use_eventplot=True)
title('Voltage and spikes over time')
ylabel("Voltage");
xlabel("Time");
ylabel('$neurons$')
ylim((-1,2))

Building finished in 0:00:01.
Simulating finished in 0:00:01.

Out[61]:
(-1, 2)
• So this is how we're going to model neurons over time
• Weaknesses of this approach
• These are "point" neurons (i.e. we're assuming the body of the neuron is just one infinitessimal point, rather than having lots and lots of compartments, each with its own RC circuit)
• Assumes R is a constant
• Assumes $V_{th}$ is a constant
• Assumes no adaptation (it's harder for a neuron to fire after it has recently fired)
• Doesn't consider nonlinearities at the input
• Each of these could be added
• Adding these increases computational cost
• Not going to worry about it now
• But we should make sure everything we do would also work with more detailed neuron models

### The time vs. rate-coding debate¶

• There is a standard ongoing debate in the neuroscience literature about whether neurons use a "rate" code or a "timing" code
• Rate code
• The important thing is the firing rate over some window of time (~100ms)
• Time code
• The precise pattern of spike generation carries information
• The NEF gets rid of this distinction
• Some people call what we do rate coding, some call it temporal coding
• The important thing is how to do the decoding usefully ... in short:
• 'Slow' decoders will seem like a rate code
• 'Fast' decoders will seem like a time code

## Temporal Decoding¶

• So that's what happens when we allow $x(t)$ to vary over time

• We keep $J(t)=\alpha e \cdot x(t) + J^{bias}$
• Feed that into the neuron model $G[J(t)]$
• We get out a sequence of spikes
• How do decode this?

• How do we get an estimate of $x$ given that sequence of spikes?
• Can we stick with linear decoding?

• Let's consider a case with 2 neurons
• For simplicity, use the same $\alpha$ and $J^{bias}$, but use $e_1=1$ and $e_2=-1$
In [62]:
#assumes has done %pylab inline
import numpy as np
import nengo
from nengo.utils.matplotlib import rasterplot

model = nengo.Network(label='Two Neurons')

with model:
stim = nengo.Node(lambda t: np.sin(10*t))
ens = nengo.Ensemble(2, dimensions=1,
encoders = [[1],[-1]],
intercepts = [-.5, -.5],
max_rates= [100, 100])
nengo.Connection(stim, ens)

stim_p = nengo.Probe(stim)
spikes_p = nengo.Probe(ens.neurons, 'spikes')

sim = nengo.Simulator(model)
sim.run(.6)

nengo.utils.matplotlib.plot_tuning_curves(ens, sim)

t = sim.trange()
figure(figsize=(12, 6))
plot(t, sim.data[stim_p],'g')
ax = gca()
ylabel("Output")
xlabel("Time");
rasterplot(t, sim.data[spikes_p], ax=ax.twinx(), use_eventplot=True)
ylabel("Neuron");

Building finished in 0:00:01.
Simulating finished in 0:00:01.

• Can we still decode $x$ using a sum of activities?

• $\hat{x}(t)=\sum_i{a_i(t) d_i}$
• Let's use exactly the same technique we did before to find $d$

• $d = \Gamma^{-1} \Upsilon$
• $\Upsilon_i = \sum_x a_i x dx$
• $\Gamma_{ij} = \sum_x a_i a_j dx$
In [63]:
import numpy as np
import nengo

model = nengo.Network(label='Decoding Neurons')

with model:
stim = nengo.Node(lambda t: np.sin(10*t))
ens = nengo.Ensemble(2, dimensions=1,
encoders = [[1],[-1]],
intercepts = [-.5, -.5],
max_rates = [100, 100])
temp = nengo.Ensemble(10, dimensions=1)

nengo.Connection(stim, ens)
connection = nengo.Connection(ens, temp) #This is just to generate the decoders

stim_p = nengo.Probe(stim)
spikes_p = nengo.Probe(ens.neurons, 'spikes')

sim = nengo.Simulator(model)
sim.run(.6)

t = sim.trange()

x = sim.data[stim_p][:,0]

A = sim.data[spikes_p]

gamma=np.dot(A.T,A)
upsilon=np.dot(A.T,x)
d = np.dot(np.linalg.pinv(gamma),upsilon)

xhat = np.dot(A, d)

figure(figsize=(8,4))
plot(t, x, label='x')
plot(t, xhat)
ylabel('$x$')
#ylim(-1,1)
xlabel('Time');

Building finished in 0:00:01.
Simulating finished in 0:00:01.

• What if we add more neurons?
In [65]:
import numpy as np
import nengo
from nengo.utils.matplotlib import rasterplot
from nengo.dists import Uniform

model = nengo.Network(label='Decoding Neurons')

N = 50

with model:
stim = nengo.Node(lambda t: np.sin(10*t))
ens = nengo.Ensemble(N, dimensions=1,
max_rates=Uniform(100,200))
temp = nengo.Ensemble(10, dimensions=1)

nengo.Connection(stim, ens)
connection = nengo.Connection(ens, temp) #This is just to generate the decoders

stim_p = nengo.Probe(stim)
spikes_p = nengo.Probe(ens.neurons, 'spikes')

sim = nengo.Simulator(model)
sim.run(.6)

x = sim.data[stim_p][:,0]

A = sim.data[spikes_p]

gamma=np.dot(A.T,A)
upsilon=np.dot(A.T,x)
d = np.dot(np.linalg.pinv(gamma),upsilon)

xhat = np.dot(A, d)

t = sim.trange()
figure(figsize=(12, 6))
ax = gca()
plot(t, sim.data[stim_p],'b')
plot(t, xhat,'g')
ylabel("Output")
xlabel("Time");
rasterplot(t, sim.data[spikes_p], ax=ax.twinx(), use_eventplot=True, color='k')
ylabel("Neuron");

Building finished in 0:00:01.
Simulating finished in 0:00:01.

• What's happening here?
• At a given instant in time $t$, only a very small number of neurons are firing
• If no neurons are firing at that instant, $\hat{x}=0$
• There are usually only one or two neurons firing at any given time
• As $dt$ gets smaller, the problem gets even worse
• What can we do?

## Temporal filtering¶

• We want a spike to contribute to other points in time
• So we want to convert a train of spikes into some continuous function
• Convolution:

• $\hat{x}(t)=\sum_i{a_i(t) * h(t) d_i}$
• What should we use for $h(t)$?

• Maybe a gaussian?
• $h(t) = e^{-t^2 / {2\sigma^2}}$
In [66]:
import numpy as np

dt = 0.001
sigma = 0.007
t_h = np.arange(200)*dt-0.1
h = np.exp(-t_h**2/(2*sigma**2))
h = h/norm(h,1)

figure()
plot(t_h, h)
xlabel('t')
ylabel('h(t)');

In [67]:
import nengo
from nengo.utils.matplotlib import rasterplot

model = nengo.Network(label='Decoding Neurons')

with model:
stim = nengo.Node(lambda t: np.sin(10*t))
ens = nengo.Ensemble(2, dimensions=1,
encoders = [[1],[-1]],
intercepts = [-.5, -.5],
max_rates = [100, 100])
temp = nengo.Ensemble(10, dimensions=1)

nengo.Connection(stim, ens)
connection = nengo.Connection(ens, temp) #This is just to generate the decoders

stim_p = nengo.Probe(stim)
spikes_p = nengo.Probe(ens.neurons, 'spikes')

sim = nengo.Simulator(model)
sim.run(.6)

t = sim.trange()
x = sim.data[stim_p]
d = sim.data[connection].weights.T

fspikes1 = np.convolve(sim.data[spikes_p][:,0], h, mode='same')
fspikes2 = np.convolve(sim.data[spikes_p][:,1], h, mode='same')

A = np.array([fspikes1, fspikes2]).T

xhat = np.dot(A, d)

figure(figsize=(8, 4))
ax = gca()
plot(t, sim.data[stim_p],'g')
ylabel("Output")
xlabel("Time");
rasterplot(t, sim.data[spikes_p], ax=ax.twinx(), use_eventplot=True, color='k')
ylabel("Neuron")

figure(figsize=(8,4))
plot(t, x, label='x')
plot(t, fspikes1*d[0])
plot(t, fspikes2*d[1])
ylabel('$x$')
xlabel('time (s)')

figure(figsize=(8,4))
plot(t, x, label='x')
plot(t, xhat, label='x')
ylabel('$x$')
xlabel('time (s)');

Building finished in 0:00:01.
Simulating finished in 0:00:01.

• Better, but how can we improve this?
• Different shapes for $h(t)$?
• Different $\sigma$?
• Is there a general solution?

## Optimal filters¶

• Let's try to find the optimal filter for a particular special case
• Two neurons
• Same $\alpha$ and $J^{bias}$ but different $e$ (-1 and +1)

$$\hat{x}(t)=\sum_i^N{a_i(t) * h(t) d_i}$$

• Since there are two neurons and they are opposites of each other, we can assume that $d_1=-d_2$ (i.e. they both contribute equally, but oppositely)

• $\hat{x}(t)= \sum_i{a_i(t) * h(t) d_i}$
• $\hat{x}(t)= a_1(t)*h(t)d_1 + a_2(t)*h(t)d_2$
• $\hat{x}(t) = (a_1(t)-a_2(t))*h(t) d_1$
• $\hat{x}(t) = (a_1(t)-a_2(t))*h(t)$ (since we can roll the constant into $h(t)$
• Let's call $a_1(t)-a_2(t)=r(t)$, for response

• In this case, it's a difference of spike trains

$\hat{x}(t)=r(t) * h(t)$

### Random inputs¶

• Can't just optimize for decoding a sine wave

• Needs to work for any random input
• What do we mean by a random input?

In [145]:
import numpy as np
x = np.random.uniform(-1,1,500)
plot(np.arange(500)*0.001, x);
xlim(0,.5)

Out[145]:
(0, 0.5)
• Maybe that's a bit too random
• Many signals don't change that quickly
• How do we get a random signal that doesn't change so quickly?
• I.e. a signal with a limit on its maximum frequency?
In [155]:
from nengo.processes import WhiteSignal

dt = 0.001
max_freq = 10

model = nengo.Network('White Noise')
with model:
sig = nengo.Node(output=WhiteSignal(period=1, high=max_freq, rms=0.5))

sig_p = nengo.Probe(sig)

sim = nengo.Simulator(model)
sim.run(1.0)

t = sim.trange()
x = sim.data[sig_p]

figure(figsize=(8, 4))
ax = gca()
plot(t, sim.data[sig_p])
ylabel("Output")
xlabel("Time");

Building finished in 0:00:01.
Simulating finished in 0:00:01.

• What is this doing?
• Generate random Fourier coefficients and convert back to time domain
• Randomly generate coefficients from a Gaussian distribution
• For $\omega$ values outside of some range, set to zero
In [156]:
import numpy as np

T=1.0
dt=0.001
signal = sim.data[sig_p][:,0]

sig_freqs = np.fft.fftshift(np.fft.fft(signal))

t = np.arange(int(T/dt))*dt
freq = np.arange(int(T/dt))/T - (T/dt)/2

figure()
plot(freq, np.abs(sig_freqs))
xlim(-500, 500)
xlabel('freq (Hz)')
ylabel('$|X(\omega)|$')

figure()
plot(freq, np.abs(sig_freqs))
xlim(-50, 50)
xlabel('freq (Hz)')
ylabel('$|X(\omega)|$')

figure()
plot(t, signal)
xlabel('time (s)')
ylabel('$x(t)$');

• Things to remember
• $X(\omega)$ values are complex numbers
• When we plot $X(\omega)$ we actually plot $|X(\omega)|$
• Since $x(t)$ should be real (i.e. no imaginary component), we expect $X(\omega)$ to be constrained
• This constraint is that $X(\omega)=X(-\omega)^*$
• (i.e., the complex conjugate: the real parts are equal, and the imaginary parts switch sign)

# Finding an optimal decoder¶

• $\hat{x}(t)=r(t) * h(t)$
• $r(t)$ is the response of the two neurons combined together
• $r(t)=a_1(t)-a_2(t)$
• $h(t)$ is the filter we are trying to find
• How do we find $h(t)$?
• Convolution is annoying; let's get rid of it
• One nice thing about convolution is that it turns into multiplication when you do a Fourier transform

$\hat{X}(\omega) = R(\omega)H(\omega)$

• We want to find $H(\omega)$ that minimizes the error

• $E = (X(\omega)-\hat{X}(\omega))^2$ (added up over all time/frequency)
• $E = (X(\omega)-R(\omega)H(\omega))^2$ (added up over all time/frequency)
• But, we know that our signal $x(t)$ can be written in the frequency domain as a sum of sine waves (that's the Fourier transform)

• We've already computed those for our signal $X$; that's how we made it in the first place
• So, we can write our error in terms of different frequency values $\omega$

• $E(\omega) = (X(\omega)-R(\omega)H(\omega))^2$

• Now we can take the derivative and set it equal to zero to do the minimization

• $H(\omega)= {{X(\omega)R^*(\omega)} \over {|R(\omega)|^2}}$

• Note the complex conjugate that appears there
• So now we can find $H(\omega)$ given the Fourier transform of our signal $X(\omega)$ and the Fourier transform of the spiking response $R(\omega)$

### Notes for comparison to the textbook¶

• Here we're using the convention that lower-case letters are for time-domain functions and upper-case letters are used for Fourier transfomr
• $x(t) \rightarrow X(\omega)$
• $r(t) \rightarrow R(\omega)$
• $h(t) \rightarrow H(\omega)$
• The textbook uses $A(\omega)$ instead of $X(\omega)$
• $A$ for amplitude
• This has caused confusion with $A$ as the matrix for activities of neurons, so I've switched it here
In [182]:
import nengo
from nengo.processes import WhiteSignal

dt = 0.001
max_freq = 10
T = 1.0

model = nengo.Network('White Noise')
with model:
in_stim = nengo.Node(output=WhiteSignal(T, high=max_freq, rms=0.5))
ens = nengo.Ensemble(2, dimensions=1,
encoders = [[1],[-1]],
intercepts = [-.3, -.3],
max_rates= [100, 100])
temp = nengo.Ensemble(1, dimensions=1)

nengo.Connection(in_stim, ens)
connection = nengo.Connection(ens, temp) #This is just to generate the decoders

stim_p = nengo.Probe(in_stim)
spikes_p = nengo.Probe(ens.neurons, 'spikes')
dec_p = nengo.Probe(ens)

sim = nengo.Simulator(model)
sim.run(T)

r_plot = (sim.data[spikes_p][:,0]-sim.data[spikes_p][:,1])*dt
r = sim.data[dec_p][:,0] #properly scaled r
signal = sim.data[stim_p][:,0]

t = sim.trange()
figure(figsize=(8, 6))
plot(t, sim.data[stim_p],'g', linewidth=2)
plot(t, r_plot,'r')
ylabel("Output")
xlabel("Time")
axis('tight');

Building finished in 0:00:01.
Simulating finished in 0:00:01.

In [183]:
freq = np.arange(int(T/dt)) - (T/dt)/2
omega = freq*2*numpy.pi

X = np.fft.fftshift(np.fft.fft(signal))
R = np.fft.fftshift(np.fft.fft(r))

figure()
plot(omega, np.abs(X))
xlabel('$\omega$')
ylabel('$|X(\omega)|$')
xlim(-200,200)

figure()
plot(omega, np.abs(R))
xlabel('$\omega$')
ylabel('$|R(\omega)|$');
#xlim(-200,200);

• Now we can find our optimal filter

$H(\omega)= {{X(\omega)R^*(\omega)} \over {|R(\omega)|^2}}$

In [184]:
Rmag = R*R.conjugate() #=np.abs(R)**2
#Rmag[(Rmag.shape[0]-1)/2 + 1] = (Rmag.shape[0]-1)/2 #set zero DC term
H = X*R.conjugate() / Rmag
h = np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(H))).real

figure()
plot(omega, np.abs(X*R.conjugate()))
xlim(-200,200)

figure()
plot(omega, np.abs(Rmag))

figure()
plot(omega, np.abs(H))
xlabel('$\omega$')
ylabel('$|H(\omega)|$')
xlim(-200,200)

figure()
plot(freq*dt,h)
xlabel('t')
ylabel('h(t)')
xlim(-.5,.5);

Does this new filter work?

In [185]:
fspikes1 = np.convolve(sim.data[spikes_p][:,0], h, mode='same')
fspikes2 = np.convolve(sim.data[spikes_p][:,1], h, mode='same')

A = np.array([fspikes1, fspikes2]).T
d = sim.data[connection].weights.T

xhat = np.dot(A, d)[:,0]

figure()
plot(t, signal)
plot(t, xhat)
show()
print 'RMSE:', np.average((signal-xhat)**2)

RMSE: 0.00726875227603

### Testing how well it generalizes¶

• If it's a good decoder, it should work for other inputs too
In [186]:
sim = nengo.Simulator(model)
sim.run(T)
sig = sim.data[stim_p][:,0]

fspikes1 = np.convolve(sim.data[spikes_p][:,0], h, mode='same')
fspikes2 = np.convolve(sim.data[spikes_p][:,1], h, mode='same')

A = np.array([fspikes1, fspikes2]).T
d = sim.data[connection].weights.T #not needed

xhat = np.dot(A, d)[:,0]

t = sim.trange()
figure()
plot(t, sig)
plot(t, xhat)
show()
print 'RMSE:', np.average((sig-xhat)**2)

Building finished in 0:00:01.
Simulating finished in 0:00:01.

RMSE: 0.0261412538446

• It's a few times worse on a different signal
• Not bad, but could use some improvement

### Improving the optimal decoder¶

• We're only training on a pretty small set of data
• So let's increase $T$
• This helps, but can get computationally expensive to generate more and more data
• There's also a bit of a shortcut
• We've already created a pretty long signal ($T=1.0$ seconds)
• We know the filter should probably decay towards zero
• A spike shouldn't affect values far away from when it happens
• So let's take our one existing signal and chop it up into lots of little signals
In [187]:
sig2 = np.zeros_like(sig)
r2 = np.zeros_like(r)

sig2[400:600] = sig[100:300]
r2[400:600] = r[100:300]

figure()
plot(t, sig)
plot(t, r)
vlines(0.1, -10, 10, linewidth=3)
vlines(0.3, -10, 10, linewidth=3)

figure()
plot(t, sig2)
plot(t, r2);

• Add up a whole bunch of these with different slices
• But, this introduces some weird edge effects at the beginning and ending of the window
• So we can do some sort of smoothing of the data (as described in Appendix B.3)
• We use a Gaussian because its Fourier transform is local in frequency (unlike a square)
In [188]:
sigma = 0.06 #This should be ~3x longer than the neuron memory
window = np.exp(-(t-0.5)**2/(sigma**2))

sig2w = np.roll(sig, 300)*window
r2w = np.roll(r, 300)*window

figure()
plot(t, sig)
plot(t, r)
vlines(0.1, -10, 10, linewidth=3)
vlines(0.3, -10, 10, linewidth=3)

figure()
plot(t, sig2w)
plot(t, r2w)
show()

• The last graph can be thought of as one trial, which we perform at every $\Delta t$ and then average over
• So we're going to take our original data ($x(t)$ and $r(t)$) and multiply it by another function (the Gaussian) in the time domain
• Consequently, in the frequency domain, this multiplication becomes a...

$$H(\omega)= {{(X(\omega)R^*(\omega)) * W(\omega)} \over {|R(\omega)|^2*W(\omega)}}$$

• So we're windowing and averaging in time because convolution in frequency (it's kind of unusual to see convolution in frequency), which gives an efficient method for estimation.
In [190]:
# original code ignoring DC
H = (X*R.conjugate()) / (R*R.conjugate())

# new code
sigma_t = 0.025
W2 = np.exp(-omega**2*sigma_t**2)
H = (np.convolve(X*R.conjugate(),W2, 'same')) / (np.convolve(R*R.conjugate(),W2, 'same'))

How well does this work?

In [211]:
in_stim.output = WhiteSignal(T, high=2, rms=0.5)

sim = nengo.Simulator(model)
sim.run(T)
sig1 = sim.data[stim_p][:,0]

r = sim.data[dec_p][:,0] #properly scaled r

X = np.fft.fftshift(np.fft.fft(sig1))
R = np.fft.fftshift(np.fft.fft(r))

sigma_t = 0.06
W2 = np.exp(-omega**2*sigma_t**2)
Rmag = R*R.conjugate()
Rmag[(Rmag.shape[0]-1)/2 + 1] = .1 #remove zero DC term
H = (np.convolve(X*R.conjugate(),W2, 'same')) / (np.convolve(Rmag,W2, 'same'))
h = np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(H))).real

Building finished in 0:00:01.
Simulating finished in 0:00:01.

• What does this new filter look like?
• How does it compare without windowing?
In [212]:
figure()
plot(omega, np.abs(H))
xlabel('$\omega$')
ylabel('$|H(\omega)|$')
xlim(-500, 500)

figure()
plot(t-T/2, h)
xlabel('t')
ylabel('h(t)')

figure()
plot(t-T/2, h)
xlabel('t')
ylabel('h(t)')
xlim(-0.2, 0.2)

Out[212]:
(-0.2, 0.2)
In [213]:
spikes1 = sim.data[spikes_p][:,0]
fspikes1_1 = np.convolve(sim.data[spikes_p][:,0], h, mode='same')
fspikes2_1 = np.convolve(sim.data[spikes_p][:,1], h, mode='same')
spikes_1 = (sim.data[spikes_p][:,0]-sim.data[spikes_p][:,1])*sim.dt

sim = nengo.Simulator(model)
sim.run(T)
sig2 = sim.data[stim_p][:,0]

fspikes1_2 = np.convolve(sim.data[spikes_p][:,0], h, mode='same')
fspikes2_2 = np.convolve(sim.data[spikes_p][:,1], h, mode='same')
spikes_2 = (sim.data[spikes_p][:,0]-sim.data[spikes_p][:,1])*sim.dt

A1 = numpy.array([fspikes1_1, fspikes2_1]).T
A2 = numpy.array([fspikes1_2, fspikes2_2]).T

d = sim.data[connection].weights.T

x1hat = np.dot(A1, d)[:,0]
x2hat = np.dot(A2, d)[:,0]

t = sim.trange()
figure()
plot(t, sig1)
plot(t, x1hat)
plot(t, spikes_1, 'k.')
ylim(-1.5,1.5)
show()
print 'RMSE:', np.average((sig1-x1hat)**2)

figure()
plot(t, sig2)
plot(t, x2hat)
plot(t, spikes_2, 'k.')
ylim(-1.5,1.5)
print 'RMSE:', np.average((sig2-x2hat)**2)

Building finished in 0:00:01.
Simulating finished in 0:00:01.

RMSE: 0.0111085632116
RMSE: 0.0122735844775

• It's a bit better on the original signal
• And is 'less worse' when tested on a different signal
• What about higher frequency inputs? Or lower frequency?
In [208]:
in_stim.output = WhiteSignal(T, high=10, rms=0.5)
sim = nengo.Simulator(model)
sim.run(T)
sig3 = sim.data[stim_p][:,0]

fspikes1_3 = np.convolve(sim.data[spikes_p][:,0], h, mode='same')
fspikes2_3 = np.convolve(sim.data[spikes_p][:,1], h, mode='same')
spikes_3 = (sim.data[spikes_p][:,0]-sim.data[spikes_p][:,1])*sim.dt

A3 = np.array([fspikes1_3, fspikes2_3]).T

d = sim.data[connection].weights.T

x3hat = np.dot(A3, d)[:,0]

t = sim.trange()
figure()
plot(t, sig3)
plot(t, x3hat)
plot(t, spikes_3, 'k.')
ylim(-1.5,1.5)
print 'RMSE:', np.average((sig3-x3hat)**2)

Building finished in 0:00:01.
Simulating finished in 0:00:01.
RMSE: 0.0156375697164

• The decoder works better if it is made for the right range of frequencies
• What does this tell us about rate vs timing codes?
• Doesn't make much sense without thinking about the signal properties
• These methods can account for both kinds of codes and smoothly vary between

## Optimal filters reconsidered¶

• So we have now found the optimal filter to apply to the spike data to decode out information
• Note a limitation: we've only done this for the case of 2 neurons!
• But we could use the same filter for many more neurons
• Recall what this windowed filter looks like
In [131]:
figure()
plot(t-T/2, h)
xlabel('t')
ylabel('h(t)')
xlim(-0.2, 0.2)

Out[131]:
(-0.2, 0.2)


text:
With Phoenix 1.3 well behind us now, Elixir developers understand how to approach most problems “The Elixir Way”:

• Create a well-named context module for the problem, e.g. “Payment”
• Put functions in it, which don’t expose how data is persisted
• Profit

This is great, but it can be taken even further. A Payment module is a great example, in fact. Many developers1 would scan the landscape for the best payment gateway out there2, find its companion library in Hex, and build their whole Payment module around its API. They would use the gateway library’s structs throughout their application, and name fields in the database using gateway-specific terms.

This gets the job done, but it can be shortsighted. Eventually, the stakeholders always want to do crazy things like:

• Switch payment gateways (“Authorize.net gave us a better deal!”)
• Route some payments through gateway A and others through gateway B

This kind of change is very common whenever the third party service you’re integrating is generic and commoditized, i.e. offered by more than one vendor. For example:

• Image manipulation
• Payment processing
• SMS delivery
• Email delivery
• Content-delivery networks
• Analytics
• Performance & error monitoring

If your app isn’t designed to anticipate this kind of change, you’ll be caught flatfooted whenever you need to change vendors for any of these services.

In Elixir, the best way to insulate your app from vendor changes is to isolate your vendor-specific code into Adapter modules.

Before you get to that though, you’d first define app-specific structs for your problem space. For example, you might create a MyApp.Payment.Charge struct to represent credit card charges.

You’d use these structs exclusively throughout your app, never depending on vendor-specific structs or types.

%MyApp.Payment.Charge{
gateway_id: "[id of charge on gateway]",
amount_cents: 1_000,
currency: "usd",
# ...
}

Next, you’d identify which features you need from your vendor, and codify those into an Adapter behaviour. In this case, Gateway seems to be a good name.

defmodule MyApp.Payment.Gateway do
alias MyApp.Payment.{Charge, GatewayError}

@callback create_charge(Charge.t) :: {:ok, Charge.t} | {:error, GatewayError.t}
end

You’d then implement that behaviour for your chosen vendor.

defmodule MyApp.Payment.StripeGateway do
@behaviour MyApp.Payment.Gateway

@impl true
def create_charge(charge) do
# ...
end
end

Inside your context, when you need to make a call to your vendor, you fetch an Adapter module from configuration and call it:

defmodule MyApp.Payment do
def create_charge(attrs \\ %{}) do
%MyApp.Payment.Charge{}
|> struct(attrs)
|> gateway().create_charge()
end

defp gateway do
:my_app
|> Application.get_env(__MODULE__, [])
|> Keyword.fetch!(:gateway)
end
end

Finally, you’d configure your context to use the desired Adapter.

config :my_app, MyApp.Payment,
gateway: MyApp.Payment.StripeGateway

Congratulations, you now have an adapter-based design!

## Benefits

Now that your vendor is isolated behind a behaviour, you can much more easily support multiple vendors. Beyond that, you also get a major benefit in your tests.

You can now easily implement a Test implementation of your Adapter which you only use in the :test environment. This allows you to thoroughly test all the possible error states, without having to make live API calls. This is particularly nice when your vendor doesn’t provide a good “test” version of their API.

## Examples

In conclusion, here are some stand-out libraries on Hex which implement the pattern if you want to see it in action:

• Ecto: Adapter-based relational database library
• Ueberauth: OAuth library, see the Strategy behaviour

format: plain


text:
12
Q:

# 476 ** 0 is divisible by both 3 and 11.The non zero digits in the hundred's and ten's places are respectively:

 A) 6 and 2 B) 8 and 2 C) 6 and 5 D) 8 and 5

Explanation:

Let the number  be 476ab0

476ab0 is divisible by 3

=> 4 + 7 + 6 + a + b + 0 is divisible by 3

=> 17 + a + b is divisible by 3 ------------------------(i)

476ab0 is divisible by 11

[(4 + 6 + b) -(7 + a + 0)] is 0 or divisible by 11

=> [3 + (b - a)] is 0 or divisible by 11  --------------(ii)

Substitute the values of a and b with the values given in the choices and select the values which satisfies both Equation 1 and Equation 2.

if a=6 and b=2,

17 + a + b = 17 + 6 + 2 = 25 which is not divisible by 3 --- Does not meet equation(i).Hence this is not the answer

if a=8 and b=2,

17 + a + b = 17 + 8 + 2 = 27 which is divisible by 3 --- Meet equation(i)

[3 + (b - a)] = [3 + (2 - 8)] = -3 which is neither 0 nor divisible by 11---Does not meet equation(ii).Hence this is not the answer

if a=6 and b=5,

17 + a + b = 17 + 6 + 5 = 28 which is not divisible by 3 --- Does not meet equation (i) .Hence this is not the answer

if a=8 and b=5,

17 + a + b = 17 + 8 + 5 = 30 which is divisible by 3 --- Meet equation 1

[3 + (b - a)] = [3 + (5 - 8)] = 0 ---Meet equation 2

Since these values satisfies both equation 1 and equation 2, this is the answer

Q:

Which expression is equivalent to ${\mathbf{i}}^{\mathbf{233}}$?

 A) i B) 1 C) -i D) -1

Explanation:

We know that,

0 99
Q:

The product of a number and its multiplicative inverse is

 A) 1 B) 0 C) -1 D) Infinity

Explanation:

The mutiplicative inverse of a number is nothing but a reciprocal of a number.

Now, the product of a number and its multiplicative inverse is always equal to 1.

For example :

Let the number be 15

Multiplicative inverse of 15 = 1/15

The product of a number and its multiplicative inverse is = 15 x 1/15 = 1.

0 99
Q:

25 divided by 7

 A) 3.571 B) 35.71 C) 0.351 D) 0.0357

Explanation:

Here we have 25 divided by 7.

25 will not go directly in 7

Hence, we get the result in decimals.

25/7 = 3.571.

2 117
Q:

Which of the following is not a prime number?

 A) 5 B) 11 C) 21 D) 37

Explanation:

A prime number is a whole number greater than 1 whose only factors are 1 and itself.

Factors of 5 are 1, 5

Factors of 11 are 1, 11

Factors of 21 are 1, 3, 7, 21

Factors of 37 are 1, 37.

Hence, according to the definition of a prime number, 21 is not a prime number as it has more than two factors.

1 98
Q:

What are the Multiples of 4?

The Multiples of 4 are 4, 8, 12, 16, 20, 24, 28, 32, 36, 40 upto 40. Mutiples of 4 means which can be divided by 4 leaving remainder '0'.

Common Multiples of 4 & 6 are 12, 24, 36, 48, 60 upto 60.

615
Q:

The sum of three consecutive odd numbers is 93. What is the middle number?

 A) 31 B) 33 C) 29 D) 27

Explanation:

Let the three consecutive odd numbers be x, x+2, x+4

Then,

x + x + 2 + x + 4 = 93

=> 3x + 6 = 93

=> 3x = 87

=> x = 29 => 29, 31, 33 are three consecutive odd numbers.

Therefore, the middle number is 31.

7 1262
Q:

How many numbers up to 300 and 700 are divisible by 2, 3 and 7 both together ?

 A) 7 B) 9 C) 11 D) 5

Explanation:

LCM of 2, 3, 7 is 42.

=> (700 – 300)/42 = 9  22/42 => 9 Numbers.

15 2700
Q:

In a triangle ABC, angle 'a' is the smallest one & A+7B = 112. Find the range of value 'c' ?

 A) 138 to 164 B) 125 to 147 C) 148 to 158 D) 152 to 164


text:
# zbMATH — the first resource for mathematics

##### Examples
 Geometry Search for the term Geometry in any field. Queries are case-independent. Funct* Wildcard queries are specified by * (e.g. functions, functorial, etc.). Otherwise the search is exact. "Topological group" Phrases (multi-words) should be set in "straight quotation marks". au: Bourbaki & ti: Algebra Search for author and title. The and-operator & is default and can be omitted. Chebyshev | Tschebyscheff The or-operator | allows to search for Chebyshev or Tschebyscheff. "Quasi* map*" py: 1989 The resulting documents have publication year 1989. so: Eur* J* Mat* Soc* cc: 14 Search for publications in a particular source with a Mathematics Subject Classification code (cc) in 14. "Partial diff* eq*" ! elliptic The not-operator ! eliminates all results containing the word elliptic. dt: b & au: Hilbert The document type is set to books; alternatively: j for journal articles, a for book articles. py: 2000-2015 cc: (94A | 11T) Number ranges are accepted. Terms can be grouped within (parentheses). la: chinese Find documents in a given language. ISO 639-1 language codes can also be used.

##### Operators
 a & b logic and a | b logic or !ab logic not abc* right wildcard "ab c" phrase (ab c) parentheses
##### Fields
 any anywhere an internal document identifier au author, editor ai internal author identifier ti title la language so source ab review, abstract py publication year rv reviewer cc MSC code ut uncontrolled term dt document type (j: journal article; b: book; a: book article)
Decay of weak solutions to the 2D dissipative quasi-geostrophic equation. (English) Zbl 1194.76040
Summary: We address the decay of the norm of weak solutions to the 2D dissipative quasi-geostrophic equation. When the initial data ${\theta }_{0}$ is in ${L}^{2}$ only, we prove that the ${L}^{2}$ norm tends to zero but with no uniform rate, that is, there are solutions with arbitrarily slow decay. For ${\theta }_{0}$ in ${L}^{p}\cap {L}^{2}$, with $1\le p<2$, we are able to obtain a uniform decay rate in ${L}^{2}$. We also prove that when the ${L}^{\frac{2}{2\alpha -1}}$ norm of ${\theta }_{0}$ is small enough, the ${L}^{q}$ norms, for $q>\frac{2}{2\alpha -1}$, have uniform decay rates. This result allows us to prove decay for the ${L}^{q}$ norms, for $q\ge \frac{2}{2\alpha -1}$, when ${\theta }_{0}$ is in ${L}^{2}\cap {L}^{\frac{2}{2\alpha -1}}$.

##### MSC:
 76D05 Navier-Stokes equations (fluid dynamics) 35B40 Asymptotic behavior of solutions of PDE 35Q35 PDEs in connection with fluid mechanics 76U05 Rotating fluids 86A05 Hydrology, hydrography, oceanography 86A10 Meteorology and atmospheric physics
##### References:
 [1] Berselli L.C. (2002). Vanishing viscosity limit and long-time behaviour for 2D quasi-geostrophic equations. Indiana Univ. Math. J. 51: 905–930 · Zbl 1044.35055 · doi:10.1512/iumj.2002.51.2075 [2] Cafarelli L., Kohn H. and Nirenberg L. (1982). Partial regularity of suitable weak solutions of the Navier-Stokes equations. Comm. Pure Appl. Math. 35: 771–831 · Zbl 0509.35067 · doi:10.1002/cpa.3160350604 [3] Carpio A. (1996). Large-time behaviour in incompressible Navier-Stokes equations. SIAM J. Math. Anal 27: 449–475 · Zbl 0845.76019 · doi:10.1137/S0036141093256782 [4] Carrillo J. and Ferreira L.C.F. (2007). Self-similar solutions and large time asymptotics for the dissipative quasi-geostrophic equations. Monatshefte für Math 15: 111–142 · Zbl 1126.35048 · doi:10.1007/s00605-007-0447-7 [5] Carrillo J. and Ferreira L.C.F. (2006). Convergence towards self-similar asymptotic behaviour for the dissipative quasi-geostrophic equations. Banach Center Publ. 74: 95–115 · doi:10.4064/bc74-0-5 [6] Carrillo, J., Ferreira, L.C.F.: Asymptotic behaviour for the subcritical dissipative quasi-geostrophic equations. Preprint UAB, 2006 [7] Chae D. (2003). The quasi-geostrophic equation in the Triebel-Lizorkin spaces. Nonlinearity 16(2): 479–495 · Zbl 1029.35006 · doi:10.1088/0951-7715/16/2/307 [8] Chae D. and Lee J. (2003). Global well-posedness in the super-critical dissipative quasi-geostrophic equations. Commun. Math. Phys. 233(2): 297–311 [9] Constantin, P., Cordoba, D., Wu, J.: On the critical dissipative quasi-geostrophic equation. Indiana Univ. Math. J. 50, Special Issue, 97–107 (2001) [10] Constantin P., Majda A. and Tabak E. (1994). Formation of strong fronts in the 2-D quasigeostrophic thermal active scalar. Nonlinearity 7: 1495–1533 · Zbl 0809.35057 · doi:10.1088/0951-7715/7/6/001 [11] Constantin P. and Wu J. (1999). Behaviour of solutions of 2D Quasi-geostrophic equations. SIAM J. Math. Anal. 30: 937–948 (electronic) · Zbl 0957.76093 · doi:10.1137/S0036141098337333 [12] Córdoba A. and Córdoba D. (2004). A maximum principle applied to quasi-geostrophic equations. Comm. Math. Phys. 249: 511–528 [13] Heywood, J.: Open problems in the theory of Navier-Stokes equations of viscous incompressible flow. The Navier-Stokes equations (Oberwolfach, 1988), Lecture Notes in Math. 1431, Berlin: Springer, 1990. pp. 1–22 [14] Ju N. (2004). Existence and uniqueness of the solution to the dissipative 2D quasi-geostrophic equations in the Sobolev space. Commun. Math. Phys. 251: 365–376 · Zbl 1106.35061 · doi:10.1007/s00220-004-1062-2 [15] Ju N. (2005). The maximum principle and the global attractor for the 2D dissipative quasi-geostrophic equation. Commun. Math. Phys. 255: 161–182 · Zbl 1088.37049 · doi:10.1007/s00220-004-1256-7 [16] Ju N. (2005). On the two dimensional quasi-geostrophic equations. Indiana Univ. Math. J. 54: 897–926 · Zbl 1185.35189 · doi:10.1512/iumj.2005.54.2518 [17] Kato T. (1984). Strong L p solutions of the Navier-Stokes equation in ${ℝ}^{m}$ , with applications to weak solutions. Math. Z. 187: 471–480 · Zbl 0545.35073 · doi:10.1007/BF01174182 [18] Kato T. and Fujita H. (1962). On the nonstationary Navier-Stokes system. Rend. Sem. Mat. Univ. Padova 32: 243–260 [19] Ogawa T., Rajopadhye S. and Schonbek M. (1997). Energy decay for a weak solution of the Navier-Stokes equation with slowly varying external forces. J. Funct. Anal. 144: 325–358 · Zbl 0873.35064 · doi:10.1006/jfan.1996.3011 [20] Pedloskym J. (1987). Geophysical Fluid Dynamics. Springer Verlag, New York [21] Resnick, S.: Dynamical problems in non-linear advective partial differential equarions. Ph. D. Thesis, University of Chicago, 1995 [22] Schonbek M. (1985). L 2 decay for weak solutions of the Navier-Stokes equations. Arch. Rat. Mech. Anal. 88: 209–222 · Zbl 0602.76031 · doi:10.1007/BF00752111 [23] Schonbek M. (1986). Large time behaviour of solutions to the Navier-Stokes equations. Comm. Partial Diff. Eqs. 11: 733–763 · Zbl 0607.35071 · doi:10.1080/03605308608820443 [24] Schonbek M. and Schonbek T. (2003). Asymptotic behavior to dissipative quasi-geostrophic flows. SIAM J. Math. Anal. 35: 357–375 (electronic) · Zbl 1126.76386 · doi:10.1137/S0036141002409362 [25] Schonbek M. and Schonbek T. (2005). Moments and lower bounds in the far-field of solutions to quasi-geostrophic flows. Disc. Contin. Dyn. Syst. 13: 1277–1304 · Zbl 1091.35070 · doi:10.3934/dcds.2005.13.1277 [26] Serrin, J.: The initial value problem for the Navier-Stokes equations. In: Nonlinear problems (Proc. Sympos., Madison, Wis.), Madison, WI: Oniv. ofvisc. Press, pp. 69–98, 1963 [27] Wu, J.: Dissipative quasi-geostrophic equations with L p data. Electron. J. Diff. Eq. (2001), No. 56, 13 pp. (electronic) [28] Wu J. (2002). The 2D dissipative quasi-geostrophic equation. Appl. Math. Letters 15: 925–930 · Zbl 1016.35060 · doi:10.1016/S0893-9659(02)00065-4 [29] Wu J. (2002). The quasi-geostrophic equation and its two regularizations. Comm. Partial Diff. Eqs. 27(5-6): 1161–1181 · Zbl 1012.35067 · doi:10.1081/PDE-120004898 [30] Wu J. (2004). Global solutions of the 2D dissipative quasi-geostrophic equation in Besov spaces. SIAM J. Math. Anal. 36(3): 1014–1030. (electronic) · Zbl 1083.76064 · doi:10.1137/S0036141003435576 [31] Wu J. (2005). Solutions of the 2D quasi-geostrophic equation in Hölder spaces. Nonlinear Anal. 62(4): 579–594 · Zbl 1116.35348 · doi:10.1016/j.na.2005.03.053 [32] Wu J. (2005). The two-dimensional quasi-geostrophic equation with critical or supercritical dissipation. Nonlinearity 18(1): 139–154 · Zbl 1067.35002 · doi:10.1088/0951-7715/18/1/008 [33] Zhang L. (1995). Sharp rates of decay of solutions to 2-dimensional Navier-Stokes equations. Comm. Partial Diff. Eqs. 20: 119–127 · Zbl 0823.35145 · doi:10.1080/03605309508821089


text:
# Extract posteriors and compute odds ratio in multinomial regression

Hi! I’am fitting a multinomial logistic regression with a within-subject factor (mimicry) as predictor. My response variable is a 3 level factor (perc) and I would like to see if there is an effect on each perc probability as a function of mimicry.
This is my dataset:

library(tidyverse)
library(tidybayes)
library(brms)

perc <- c("happy", "mixed", "neutral")
mimicry <- c("blocked", "free")
subject <- 1:10
trial <- 1:10

dat <- expand_grid(subject, trial, perc, mimicry)

dat\$resp <- sample(perc, nrow(dat), replace = T)

Then I fit the model:

fit_ex <- brm(resp ~ mimicry + (1|subject),
data = dat,
family = categorical(link = "logit"))

summary(fit_exp)

Now my model parameters are not directly what I’m looking for because I would like to have the odds ratio between let’s say happy (from perc) in the mimicry|blocked condition and the mimicry|free condition. So using the tidybayes package I have extracted the posteriors and computed the odds ratio as \frac{\frac{p(happy|free)}{1-p(happy|free)}}{\frac{p(happy|blocked)}{1-p(happy|blocked)}}

new_dat <- expand_grid(perc, mimicry)

post <- add_fitted_draws(new_dat, fit_ex, re_formula = NA)

post %>%
group_by(mimicry, .category, .draw) %>%
summarise(.value = mean(.value)) %>%
ungroup() %>%
pivot_wider(names_from = mimicry, values_from = .value) %>%
mutate(odds_ratio = (blocked/(1-blocked))/(free/(1-free))) %>%
group_by(.category) %>%
median_qi(odds_ratio)

Is that a correct approach using the posterior distributions?
Thanks

Hi there,

I initially thought something was off with your code but I was thinking of multinomial rather than categorical models (they can be confusing). But for categorical I think you are correct. The key is to compute your derived variable of interest for each draw and then summarise at the end, which is what you did.

I think the only thing I’d change is that your new_dat doesn’t need perc. You can just do

new_dat <- tibble(mimicry=mimicry)
post <- add_fitted_draws(new_dat, fit_ex, re_formula = NA)

This way there’s also no need to call summarise . Depending on your ultimate goal you may also choose to include the uncertainty from the subject effects with something like

new_dat <- expand_grid(mimicry=mimicry, subject=1000) # any new subject level
post <- add_fitted_draws(new_dat, fit_ex, re_formul = NULL, allow_new_levels=TRUE)

Hugo

EDIT: PS, tagging the post with ‘brms’ may be a good idea next time :)

2 Likes

Thanks Hugo! Yes I have realized that multinomial and categorical are a distinct type of stuff. Thanks also for the code improvement suggestion! I would like to ask you a clarification about this part

new_dat <- expand_grid(mimicry=mimicry, subject=1000) # any new subject level

What’s the role of the subject inside the newdata object?

If you wanted to ‘predict’, in a sense, or see how things may look for a new subject you can add a new ‘level’ to the grouping variable and then use one of the sample_new_levels options to account for this uncertainty . In your simulated data it hardly makes a difference but in real data it will likely widen your 95% HDI.

A nice discussion here: https://github.com/paul-buerkner/brms/issues/82#issuecomment-231440994

And a recent update:


text:
Goto Chapter: Top 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 Bib Ind

### 9 Bistellar flips

#### 9.1 Theory

Since two combinatorial manifolds are already considered distinct to each other as soon as they are not combinatorially isomorphic, a topological PL-manifold is represented by a whole class of combinatorial manifolds. Thus, a frequent question when working with combinatorial manifolds is whether two such objects are PL-homeomorphic or not. One possibility to approach this problem, i. e. to find combinatorially distinct members of the class of a PL-manifold, is a heuristic algorithm using the concept of bistellar moves.

(Bistellar moves [Pac87])
Let M be a combinatorial d-manifold (d-pseudomanifold), γ = ⟨ v_0 , ... , v_k ⟩ a k-face and δ = ⟨ w_0 , ... , w_d-k ⟩ a (d-k+1)-tuple of vertices of M that does not span a (d-k)-face in M, 0 ≤ k ≤ d, such that { v_0 , ..., v_k } ∩ { w_0 , ... , w_d-k } = ∅ and { v_0 , ... , v_k, w_0 , ... w_k-d } spans exactly d-k+1 facets. Then the operation

κ_(γ,δ) ( M ) = M ∖ (γ ⋆ ∂ δ) ∪ (∂ γ ⋆ δ)

is called a bistellar (d-k)-move.

In other words: If there exists a bouquet D ⊂ M of d-k+1 facets on a subset of vertices W ⊂ V of order d+2 with a common k-face γ and the complement δ of the vertices of γ in W does not span a (d-k)-face in M we can remove D and replace it by a bouquet of k+1 facets E ⊂ M with vertex set W with a common face spanned by δ. By construction ∂ D = ∂ E and the altered complex is again a combinatorial d-manifold (d-pseudomanifold). See Fig. 11 for a bistellar 1-move of a 2-dimensional complex, see Fig. 12 for all bistellar moves in dimension 3.

A bistellar 0-move is a stellar subdivision, i. e. the subdivision of a facet δ into d+1 new facets by introducing a new vertex at the center of δ (cf. Fig. 12 on the left). In particular, the vertex set of a combinatorial manifold (pseudomanifold) is not invariant under bistellar moves. For any bistellar (d-k)-move κ_(γ,δ) we have an inverse bistellar k-move κ^-1_(γ,δ) = κ_(δ,γ) such that κ_(δ,γ) ( κ_(γ,δ) (M)) = M. If for two combinatorial manifolds M and N there exist a sequence of bistellar moves that transforms one into the other, M and N are called bistellarly equivalent. So far bistellar moves are local operations on combinatorial manifolds that change its combinatorial type. However, the strength of the concept in combinatorial topology is a consequence of the following

(Bistellar moves [Pac87])
Two combinatorial manifolds (pseudomanifolds) M and N are PL homeomorphic if and only if they are bistellarly equivalent.

Unfortunately Pachners theorem does not guarantee that the search for a connecting sequence of bistellar moves between M and N terminates. Hence, using bistellar moves, we can not prove that M and N are not PL-homeomorphic. However, there is a very effective simulated annealing approach that is able to give a positive answer in a lot of cases. The heuristic was first implemented by Bjoerner and Lutz in [BL00]. The functions presented in this chapter are based on this code which can be used for several tasks:

In many cases the heuristic reduces a given triangulation but does not reach a minimal triangulation after a reasonable amount of flips. Thus, we usually can not expect the algorithm to terminate. However, in some cases the program normally stops after a small number of flips:

Technical note: Since bistellar flips do not respect the combinatorial properties of a complex, no attention to the original vertex labels is payed, i. e. the flipped complex will be relabeled whenever its vertex labels become different from the standard labeling (for example after every reverse 0-move).

#### 9.2 Functions for bistellar flips

##### 9.2-1 SCBistellarOptions
 ‣ SCBistellarOptions ( global variable )

Record of global variables to adjust output an behavior of bistellar moves in SCIntFunc.SCChooseMove (9.2-4) and SCReduceComplexEx (9.2-14) respectively.

1. BaseRelaxation: determines the length of the relaxation period. Default: 3

2. BaseHeating: determines the length of the heating period. Default: 4

3. Relaxation: value of the current relaxation period. Default: 0

4. Heating: value of the current heating period. Default: 0

5. MaxRounds: maximal over all number of bistellar flips that will be performed. Default: 500000

6. MaxInterval: maximal number of bistellar flips that will be performed without a change of the f-vector of the moved complex. Default: 100000

7. Mode: flip mode, 0=reducing, 1=comparing, 2=reduce as sub-complex, 3=randomize. Default: 0

8. WriteLevel: 0=no output, 1=storing of every vertex minimal complex to user library, 2=e-mail notification. Default: 1

9. MailNotifyIntervall: (minimum) number of seconds between two e-mail notifications. Default: 24 ⋅ 60 ⋅ 60 (one day)

10. MaxIntervalIsManifold: maximal number of bistellar flips that will be performed without a change of the f-vector of a vertex link while trying to prove that the complex is a combinatorial manifold. Default: 5000

11. MaxIntervalRandomize := 50: number of flips performed to create a randomized sphere. Default: 50

 gap> SCBistellarOptions.BaseRelaxation;
3
gap> SCBistellarOptions.BaseHeating;
4
gap> SCBistellarOptions.Relaxation;
0
gap> SCBistellarOptions.Heating;
0
gap> SCBistellarOptions.MaxRounds;
500000
gap> SCBistellarOptions.MaxInterval;
100000
gap> SCBistellarOptions.Mode;
0
gap> SCBistellarOptions.WriteLevel;
0
gap> SCBistellarOptions.MailNotifyInterval;
86400
gap> SCBistellarOptions.MaxIntervalIsManifold;
5000
gap> SCBistellarOptions.MaxIntervalRandomize;
50

##### 9.2-2 SCEquivalent
 ‣ SCEquivalent( complex1, complex2 ) ( method )

Returns: true or false upon success, fail or a list of type [ fail, SCSimplicialComplex, Integer, facet list] otherwise.

Checks if the simplicial complex complex1 (which has to fulfill the weak pseudomanifold property with empty boundary) can be reduced to the simplicial complex complex2 via bistellar moves, i. e. if complex1 and complex2 are PL-homeomorphic. Note that in general the problem is undecidable. In this case fail is returned.

It is recommended to use a minimal triangulation complex2 for the check if possible.

Internally calls SCReduceComplexEx (9.2-14) (complex1,complex2,1,SCIntFunc.SCChooseMove);

 gap> SCBistellarOptions.WriteLevel:=0;; # do not save complexes to disk
gap> obj:=SC([[1,2],[2,3],[3,4],[4,5],[5,6],[6,1]]);; # hexagon
gap> refObj:=SCBdSimplex(2);; # triangle as a (minimal) reference object
gap> SCEquivalent(obj,refObj);
#I  SCReduceComplexEx: complexes are bistellarly equivalent.
true

##### 9.2-3 SCExamineComplexBistellar
 ‣ SCExamineComplexBistellar( complex ) ( method )

Returns: simplicial complex passed as argument with additional properties upon success, fail otherwise.

Computes the face lattice, the f-vector, the AS-determinant, the dimension and the maximal vertex label of complex.

 gap> obj:=SC([[1,2],[2,3],[3,4],[4,5],[5,6],[6,1]]);
[SimplicialComplex

Properties known: Dim, FacetsEx, Name, Vertices.

Name="unnamed complex 7"
Dim=1

/SimplicialComplex]
gap> SCExamineComplexBistellar(obj);
[SimplicialComplex

Properties known: AltshulerSteinberg, BoundaryEx, Dim, FVector,
FacetsEx, HasBoundary, IsPseudoManifold, IsPure,
Name, NumFaces[], SkelExs[], Vertices.

Name="unnamed complex 7"
Dim=1
AltshulerSteinberg=0
FVector=[ 6, 6 ]
HasBoundary=false
IsPseudoManifold=true
IsPure=true

/SimplicialComplex]

##### 9.2-4 SCIntFunc.SCChooseMove
 ‣ SCIntFunc.SCChooseMove( dim, moves ) ( function )

Returns: a bistellar move, i. e. a pair of lists upon success, fail otherwise.

Since the problem of finding a bistellar flip sequence that reduces a simplicial complex is undecidable, we have to use an heuristic approach to choose the next move.

The implemented strategy SCIntFunc.SCChooseMove first tries to directly remove vertices, edges, i-faces in increasing dimension etc. If this is not possible it inserts high dimensional faces in decreasing co-dimension. To do this in an efficient way a number of parameters have to be adjusted, namely SCBistellarOptions.BaseHeating and SCBistellarOptions.BaseRelaxation. See SCBistellarOptions (9.2-1) for further options.

If this strategy does not work for you, just implement a customized strategy and pass it to SCReduceComplexEx (9.2-14).

See SCRMoves (9.2-10) for further information.

##### 9.2-5 SCIsKStackedSphere
 ‣ SCIsKStackedSphere( complex, k ) ( method )

Returns: a list upon success, fail otherwise.

Checks, whether the given simplicial complex complex that must be a PL d-sphere is a k-stacked sphere with 1≤ k≤ ⌊fracd+22⌋ using a randomized algorithm based on bistellar moves (see [Eff11b], [Eff11a]). Note that it is not checked whether complex is a PL sphere -- if not, the algorithm will not succeed. Returns a list upon success: the first entry is a boolean, where true means that the complex is k-stacked and false means that the complex cannot be k-stacked. A value of -1 means that the question could not be decided. The second argument contains a simplicial complex that, in case of success, contains the trigangulated (d+1)-ball B with ∂ B=S and operatornameskel_d-k(B)=operatornameskel_d-k(S), where S denotes the simplicial complex passed in complex.

Internally calls SCReduceComplexEx (9.2-14).

 gap> SCLib.SearchByName("S^4~S^1");
[ [ 463, "S^4~S^1 (VT)" ], [ 1474, "S^4~S^1 (VT)" ], [ 1475, "S^4~S^1 (VT)" ],
[ 2477, "S^4~S^1 (VT)" ], [ 4395, "S^4~S^1 (VT)" ],
[ 4396, "S^4~S^1 (VT)" ], [ 4397, "S^4~S^1 (VT)" ],
[ 4398, "S^4~S^1 (VT)" ], [ 4399, "S^4~S^1 (VT)" ],
[ 4402, "S^4~S^1 (VT)" ], [ 4403, "S^4~S^1 (VT)" ],
[ 4404, "S^4~S^1 (VT)" ] ]
gap> c:=SCLib.Load(last[1][1]);;
gap> l:=SCLink(c,1);
[SimplicialComplex

Properties known: Dim, FacetsEx, Name, Vertices.

Name="lk([ 1 ]) in S^4~S^1 (VT)"
Dim=4

/SimplicialComplex]
gap> SCIsKStackedSphere(l,1);
#I  SCIsKStackedSphere: checking if complex is a 1-stacked sphere...
#I  SCIsKStackedSphere: try 1/1
#I  SCIsKStackedSphere: complex is a 1-stacked sphere.
[ true, [SimplicialComplex

Properties known: Dim, FacetsEx, Name, Vertices.

Name="Filled 1-stacked sphere (lk([ 1 ]) in S^4~S^1 (VT))"
Dim=5

/SimplicialComplex] ]

##### 9.2-6 SCBistellarIsManifold
 ‣ SCBistellarIsManifold( complex ) ( method )

Returns: true or false upon success, fail otherwise.

Tries to prove that a closed simplicial d-pseudomanifold is a combinatorial manifold by reducing its vertex links to the boundary of the d-simplex.

false is returned if it can be proven that there exists a vertex link which is not PL-homeomorphic to the standard PL-sphere, true is returned if all vertex links are bistellarly equivalent to the boundary of the simplex, fail is returned if the algorithm does not terminate after the number of rounds indicated by SCBistellarOptions.MaxIntervallIsManifold.

Internally calls SCReduceComplexEx (9.2-14) (link,SCEmpty(),0,SCIntFunc.SCChooseMove); for every link of complex. Note that false is returned in case of a bounded manifold.

See SCIsManifoldEx (12.1-18) and SCIsManifold (12.1-17) for alternative methods for manifold verification.

 gap> c:=SCBdCrossPolytope(3);;
gap> SCBistellarIsManifold(c);
true

##### 9.2-7 SCIsMovableComplex
 ‣ SCIsMovableComplex( complex ) ( method )

Returns: true or false upon success, fail otherwise.

Checks if a simplicial complex complex can be modified by bistellar moves, i. e. if it is a pure simplicial complex which fulfills the weak pseudomanifold property with empty boundary.

 gap> c:=SCBdCrossPolytope(3);;
gap> SCIsMovableComplex(c);
true

Complex with non-empty boundary:

 gap> c:=SC([[1,2],[2,3],[3,4],[3,1]]);;
gap> SCIsMovableComplex(c);
false

##### 9.2-8 SCMove
 ‣ SCMove( c, move ) ( method )

Returns: simplicial complex of type SCSimplicialComplex upon success, fail otherwise.

Applies the bistellar move move to a simplicial complex c. move is given as a (r+1)-tuple together with a (d+1-r)-tuple if d is the dimension of c and if move is a r-move. See SCRMoves (9.2-10) for detailed information about bistellar r-moves.

Note: move and c should be given in standard labeling to ensure a correct result.

 gap> obj:=SC([[1,2],[2,3],[3,4],[4,1]]);
[SimplicialComplex

Properties known: Dim, FacetsEx, Name, Vertices.

Name="unnamed complex 5"
Dim=1

/SimplicialComplex]
gap> moves:=SCMoves(obj);
[ [ [ [ 1, 2 ], [  ] ], [ [ 1, 4 ], [  ] ], [ [ 2, 3 ], [  ] ],
[ [ 3, 4 ], [  ] ] ],
[ [ [ 1 ], [ 2, 4 ] ], [ [ 2 ], [ 1, 3 ] ], [ [ 3 ], [ 2, 4 ] ],
[ [ 4 ], [ 1, 3 ] ] ] ]
gap> obj:=SCMove(obj,last[2][1]);
[SimplicialComplex

Properties known: Dim, FacetsEx, Name, NumFaces[], SkelExs[],
Vertices.

Name="unnamed complex 6"
Dim=1

/SimplicialComplex]

##### 9.2-9 SCMoves
 ‣ SCMoves( complex ) ( method )

Returns: a list of list of pairs of lists upon success, fail otherwise.

See SCRMoves (9.2-10) for further information.

 gap> c:=SCBdCrossPolytope(3);;
gap> moves:=SCMoves(c);
[ [ [ [ 1, 3, 5 ], [  ] ], [ [ 1, 3, 6 ], [  ] ], [ [ 1, 4, 5 ], [  ] ],
[ [ 1, 4, 6 ], [  ] ], [ [ 2, 3, 5 ], [  ] ], [ [ 2, 3, 6 ], [  ] ],
[ [ 2, 4, 5 ], [  ] ], [ [ 2, 4, 6 ], [  ] ] ],
[ [ [ 1, 3 ], [ 5, 6 ] ], [ [ 1, 4 ], [ 5, 6 ] ], [ [ 1, 5 ], [ 3, 4 ] ],
[ [ 1, 6 ], [ 3, 4 ] ], [ [ 2, 3 ], [ 5, 6 ] ], [ [ 2, 4 ], [ 5, 6 ] ],
[ [ 2, 5 ], [ 3, 4 ] ], [ [ 2, 6 ], [ 3, 4 ] ], [ [ 3, 5 ], [ 1, 2 ] ],
[ [ 3, 6 ], [ 1, 2 ] ], [ [ 4, 5 ], [ 1, 2 ] ], [ [ 4, 6 ], [ 1, 2 ] ] ]
, [  ] ]

##### 9.2-10 SCRMoves
 ‣ SCRMoves( complex, r ) ( method )

Returns: a list of pairs of the form [ list, list ], fail otherwise.

A bistellar r-move of a d-dimensional combinatorial manifold complex is a r-face m_1 together with a d-r-tuple m_2 where m_1 is a common face of exactly (d+1-r) facets and m_2 is not a face of complex.

The r-move removes all facets containing m_1 and replaces them by the (r+1) faces obtained by uniting m_2 with any subset of m_1 of order r.

The resulting complex is PL-homeomorphic to complex.

 gap> c:=SCBdCrossPolytope(3);;
gap> moves:=SCRMoves(c,1);
[ [ [ 1, 3 ], [ 5, 6 ] ], [ [ 1, 4 ], [ 5, 6 ] ], [ [ 1, 5 ], [ 3, 4 ] ],
[ [ 1, 6 ], [ 3, 4 ] ], [ [ 2, 3 ], [ 5, 6 ] ], [ [ 2, 4 ], [ 5, 6 ] ],
[ [ 2, 5 ], [ 3, 4 ] ], [ [ 2, 6 ], [ 3, 4 ] ], [ [ 3, 5 ], [ 1, 2 ] ],
[ [ 3, 6 ], [ 1, 2 ] ], [ [ 4, 5 ], [ 1, 2 ] ], [ [ 4, 6 ], [ 1, 2 ] ] ]

##### 9.2-11 SCRandomize
 ‣ SCRandomize( complex[[, rounds][, allowedmoves]] ) ( function )

Returns: a simplicial complex upon success, fail otherwise.

Randomizes the given simplicial complex complex via bistellar moves chosen at random. By passing the optional array allowedmoves, which has to be a dense array of integer values of length SCDim(complex), certain moves can be allowed or forbidden in the proccess. An entry allowedmoves[i]=1 allows (i-1)-moves and an entry allowedmoves[i]=0 forbids (i-1)-moves in the randomization process.

With optional positive integer argument rounds, the amount of randomization can be controlled. The higher the value of rounds, the more bistellar moves will be randomly performed on complex. Note that the argument rounds overrides the global setting SCBistellarOptions.MaxIntervalRandomize (this value is used, if rounds is not specified). Internally calls SCReduceComplexEx (9.2-14).

 gap> c:=SCRandomize(SCBdSimplex(4));
[SimplicialComplex

Properties known: Dim, FacetsEx, Name, Vertices.

Name="Randomized S^3_5"
Dim=3

/SimplicialComplex]
gap> c.F;
[ 16, 65, 98, 49 ]

##### 9.2-12 SCReduceAsSubcomplex
 ‣ SCReduceAsSubcomplex( complex1, complex2 ) ( method )

Returns: SCBistellarOptions.WriteLevel=0: a triple of the form [ boolean, simplicial complex, rounds performed ] upon termination of the algorithm.

SCBistellarOptions.WriteLevel=1: A library of simplicial complexes with a number of complexes from the reducing process and (upon termination) a triple of the form [ boolean, simplicial complex, rounds performed ].

SCBistellarOptions.WriteLevel=2: A mail in case a smaller version of complex1 was found, a library of simplicial complexes with a number of complexes from the reducing process and (upon termination) a triple of the form [ boolean, simplicial complex, rounds performed ].

Returns fail upon an error.

Reduces a simplicial complex complex1 (satisfying the weak pseudomanifold property with empty boundary) as a sub-complex of the simplicial complex complex2.

Main application: Reduce a sub-complex of the cross polytope without introducing diagonals.

Internally calls SCReduceComplexEx (9.2-14) (complex1,complex2,2,SCIntFunc.SCChooseMove);

 gap> c:=SCFromFacets([[1,3],[3,5],[4,5],[4,1]]);;
gap> SCBistellarOptions.WriteLevel:=0;; # do not save complexes
gap> SCReduceAsSubcomplex(c,SCBdCrossPolytope(3));
[ true, [SimplicialComplex

Properties known: Dim, FacetsEx, Name, Vertices.

Name="unnamed complex 36"
Dim=1

/SimplicialComplex], 1 ]

##### 9.2-13 SCReduceComplex
 ‣ SCReduceComplex( complex ) ( method )

Returns: SCBistellarOptions.WriteLevel=0: a triple of the form [ boolean, simplicial complex, rounds performed ] upon termination of the algorithm.

SCBistellarOptions.WriteLevel=1: A library of simplicial complexes with a number of complexes from the reducing process and (upon termination) a triple of the form [ boolean, simplicial complex, rounds performed ].

SCBistellarOptions.WriteLevel=2: A mail in case a smaller version of complex1 was found, a library of simplicial complexes with a number of complexes from the reducing process and (upon termination) a triple of the form [ boolean, simplicial complex, rounds performed ].

Returns fail upon an error..

Reduces a pure simplicial complex complex satisfying the weak pseudomanifold property via bistellar moves. Internally calls SCReduceComplexEx (9.2-14) (complex,SCEmpty(),0,SCIntFunc.SCChooseMove);

 gap> obj:=SC([[1,2],[2,3],[3,4],[4,5],[5,6],[6,1]]);; # hexagon
gap> SCBistellarOptions.WriteLevel:=0;; # do not save complexes
gap> tmp := SCReduceComplex(obj);
[ true, [SimplicialComplex

Properties known: Dim, FacetsEx, Name, Vertices.

Name="unnamed complex 27"
Dim=1

/SimplicialComplex], 3 ]

##### 9.2-14 SCReduceComplexEx
 ‣ SCReduceComplexEx( complex, refComplex, mode, choosemove ) ( function )

Returns: SCBistellarOptions.WriteLevel=0: a triple of the form [ boolean, simplicial complex, rounds ] upon termination of the algorithm.

SCBistellarOptions.WriteLevel=1: A library of simplicial complexes with a number of complexes from the reducing process and (upon termination) a triple of the form [ boolean, simplicial complex, rounds ].

SCBistellarOptions.WriteLevel=2: A mail in case a smaller version of complex1 was found, a library of simplicial complexes with a number of complexes from the reducing process and (upon termination) a triple of the form [ boolean, simplicial complex, rounds ].

Returns fail upon an error.

Reduces a pure simplicial complex complex satisfying the weak pseudomanifold property via bistellar moves mode = 0, compares it to the simplicial complex refComplex (mode = 1) or reduces it as a sub-complex of refComplex (mode = 2).

choosemove is a function containing a flip strategy, see also SCIntFunc.SCChooseMove (9.2-4).

The currently smallest complex is stored to the variable minComplex, the currently smallest f-vector to minF. Note that in general the algorithm will not stop until the maximum number of rounds is reached. You can adjust the maximum number of rounds via the property SCBistellarOptions (9.2-1). The number of rounds performed is returned in the third entry of the triple returned by this function.

This function is called by

1. SCReduceComplex (9.2-13),

2. SCEquivalent (9.2-2),

3. SCReduceAsSubcomplex (9.2-12),

4. SCBistellarIsManifold (9.2-6).

5. SCRandomize (9.2-11).

Please see SCMailIsPending (15.2-3) for further information about the email notification system in case SCBistellarOptions.WriteLevel is set to 2.

 gap> c:=SCBdCrossPolytope(4);;
gap> SCBistellarOptions.WriteLevel:=0;; # do not save complexes
gap> SCReduceComplexEx(c,SCEmpty(),0,SCIntFunc.SCChooseMove);
[ true, [SimplicialComplex

Properties known: Dim, FacetsEx, Name, Vertices.

Name="unnamed complex 13"
Dim=3

/SimplicialComplex], 7 ]
gap> SCReduceComplexEx(c,SCEmpty(),0,SCIntFunc.SCChooseMove);
[ true, [SimplicialComplex

Properties known: Dim, FacetsEx, Name, Vertices.

Name="unnamed complex 18"
Dim=3

/SimplicialComplex], 9 ]
gap> SCMailSetAddress("johndoe@somehost");
true
gap> SCMailIsEnabled();
true
gap> SCReduceComplexEx(c,SCEmpty(),0,SCIntFunc.SCChooseMove);
[ true, [SimplicialComplex

Properties known: Dim, FacetsEx, Name, Vertices.

Name="unnamed complex 23"
Dim=3

/SimplicialComplex], 7 ]

Content of sent mail:

Greetings master,

this is simpcomp 2.1.9 running on comp01.maths.fancytown.edu

I have been working hard for 0 seconds and have a message for you, see below.

#### START MESSAGE ####

SCReduceComplex:

Computed locally minimal complex after 7 rounds:

[SimplicialComplex

Properties known: Boundary, Chi, Date, Dim, F, Faces, Facets, G, H,
HasBoundary, Homology, IsConnected, IsManifold, IsPM, Name, SCVertices,
Vertices.

Name="ReducedComplex_5_vertices_7"
Dim=3
Chi=0
F=[ 5, 10, 10, 5 ]
G=[ 0, 0 ]
H=[ 1, 1, 1, 1 ]
HasBoundary=false
Homology=[ [ 0, [ ] ], [ 0, [ ] ], [ 0, [ ] ], [ 1, [ ] ] ]
IsConnected=true
IsPM=true

/SimplicialComplex]

##### END MESSAGE #####

That's all, I hope this is good news! Have a nice day.

##### 9.2-15 SCReduceComplexFast
 ‣ SCReduceComplexFast( complex ) ( function )

Returns: a simplicial complex upon success, fail otherwise.

Same as SCReduceComplex (9.2-13), but calls an external binary provided with the simpcomp package.

Goto Chapter: Top 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 Bib Ind

generated by GAPDoc2HTML


text:
# When should normalization never be used?

Lately, there have been numerous questions about normalization

What are some of the situations where you never ever ever should normalize your data, and what are the alternatives?

• What do you precisely mean by normalization? – user88 Jul 27 '10 at 20:13

Whether one can normalize a non-normal data set depends on the application. For example, data normalization is required for many statistical tests (i.e. calculating a z-score, t-score, etc.) Some tests are more prone to failure when normalizing non-normal data, while some are more resistant ("robust" tests).

One less-robust statistic is the mean, which is sensitive to outliers (i.e. non-normal data). Alternatively, the median is less sensitive to outliers (and therefore more robust).

A great example of non-normal data when many statistics fail is bi-modally distributed data. Because of this, it's always good practice to visualize your data as a frequency distribution (or even better, test for normality!)

Of course one should never try to blindly normalize data if the data does not follow a (single) normal distribution.

For example one might want to rescale observables $X$ to all be normal with $(X-\mu)/\sigma$, but this can only work if the data is normal and if both $\mu$ and $\sigma$ are the same for all data points (e.g. $\sigma$ doesn't depend on $\mu$ in a particular $X$ range).

I thought this was too obvious, until I saw this question!

When you normalise data, make sure you always have access to the raw data after normalisation. Of course, you could break this rule if you have a good reason, e.g. storage.


text:
m
• D

F Nous contacter

0

# Vidéo de la semaine

## Post-edited

Enriques-Kodaira classification Kodaira dimension quartics complete intersections double plane

# Dernières Vidéos

## Post-edited Bayesian inference and mathematical imaging - Part 3: probability and convex optimisation

This course presents an overview of modern Bayesian strategies for solving imaging inverse problems. We will start by introducing the Bayesian statistical decision theory framework underpinning ...

## Post-edited On the existence of algebraic approximations of compact Kähler manifolds

Let $X$ be a compact Kähler manifold. The so-called Kodaira problem asks whether $X$ has arbitrarily small deformations to some projective varieties. While Kodaira proved that such deformations ...

## Post-edited Intermittent weak solutions of the 3D Navier-Stokes equations

I will discuss recent developments concerning the non-uniqueness of distributional solutions to the Navier-Stokes equation.

## Multi angle Topics on $K3$ surfaces - Lecture 5: Finite automorphism groups

Aim of the lecture is to give an introduction to $K3$ surfaces, that are special algebraic surfaces with an extremely rich geometry. The most easy example of such a surface is the Fermat quartic in ...

## Multi angle Topics on $K3$ surfaces - Lecture 3: Basic properties of $K3$ surfaces

Aim of the lecture is to give an introduction to $K3$ surfaces, that are special algebraic surfaces with an extremely rich geometry. The most easy example of such a surface is the Fermat quartic in ...

## Multi angle Topics on $K3$ surfaces - Lecture 6: Classification

Aim of the lecture is to give an introduction to $K3$ surfaces, that are special algebraic surfaces with an extremely rich geometry. The most easy example of such a surface is the Fermat quartic in ...

## Multi angle Topics on $K3$ surfaces - Lecture 4: Nèron-Severi group and automorphisms

Aim of the lecture is to give an introduction to $K3$ surfaces, that are special algebraic surfaces with an extremely rich geometry. The most easy example of such a surface is the Fermat quartic in ...
Z


text:
# Does a random sequence of vectors span a Hilbert space?

Let $$\mathcal{H}$$ be a separable Hilbert space. Let $$v$$ be a random variable taking values in $$\mathcal{H}$$ such that $$P(v \perp h) < 1$$ for all $$h \in \mathcal{H}.$$ Suppose we sample an infinite sequence $$v_1, v_2, \ldots.$$ Is it the case that, almost surely, the closed span of $$v_1, v_2, \ldots$$ is all of $$\mathcal{H}?$$

• and if they are iid the probability of being in a closed hyperplane $(h)^\perp$ is $P(v_k\perp h, k=1,2,\dots)=0$ – Pietro Majer Apr 21 '19 at 17:07
• @Pietro Majer: this is the probability that the vectors all lie in a given closed hyperplane. – Anthony Quas Apr 21 '19 at 17:10
• Whoever voted to close this, I'm pretty sure you don't understand the question. This is subtle and interesting. – Anthony Quas Apr 21 '19 at 17:18
• Maybe I'm missing something obvious, but is it even clear that the event you are interested in is measurable? – Jochen Glueck Apr 21 '19 at 17:21
• @JochenGlueck: Yes it's measurable: Let $(y_n)$ be a dense sequence in $\mathcal H$. Then the event is: for all $m>0$, for all $n>0$, there exist $k>0$ and rational $t_1,\ldots,t_k$ such that $\|t_1v_1+\ldots+t_kv_k-y_n\|<1/m$. – Anthony Quas Apr 21 '19 at 19:18

(This may turn out to be a simplified version of J. E. Pascoe's answer).

The support of (the distribution of) $$v$$, that we denote by $$\operatorname{supp} v$$, is the set of vectors $$h \in \mathcal{H}$$ such that $$P(v \in B(h, \varepsilon)) > 0$$ for every $$\varepsilon > 0$$. We list some properties of this set.

1. The set $$\operatorname{supp} v$$ is the complement of the union of all open sets $$B$$ such that $$P(v \in B) = 0$$. Thus, the support is a closed set.

2. Since $$\mathcal{H}$$ is a separable metric space, it has a countable topological base $$\mathcal{B}$$, and $$\operatorname{supp} v$$ is the complement of the union of all $$B \in \mathcal{B}$$ such that $$P(v \in B) = 0$$. By countable additivity, it follows that $$P(v \in \operatorname{supp} v) = 1$$ (the support is a set of full measure).

3. With probability one, the closure of the random set $$V = \{v_1, v_2, \ldots\}$$ contains $$\operatorname{supp} v$$. Indeed, let $$\{h_1, h_2, \ldots\}$$ be a countable, dense subset of $$\operatorname{supp} v$$. For every $$i, n = 1, 2, \ldots$$ we have $$P(v \in B(h_i, \tfrac{1}{n})) > 0$$, and thus, by Borel–Cantelli, $$P(V \cap B(h_i, \tfrac{1}{n}) = \varnothing) = 0$$. It follows that $$h_i \in \overline{V}$$ for every $$i = 1, 2, \ldots$$, and consequently $$\operatorname{supp} v \subseteq \overline{V}$$.

4. For every $$h \in \mathcal{H}$$, we have $$P(h \perp v) < 1$$, and therefore $$h$$ is not orthogonal to $$\operatorname{supp} v$$. It follows that the closed span of $$\operatorname{supp} v$$ is $$\mathcal{H}$$.

It remains to note that the closed span of $$V$$ is the same as the closed span of the closure of $$V$$, which with probability one contains the closed span of $$\operatorname{supp} v$$, which we have shown to be equal to $$\mathcal{H}$$.

(Item 1 is valid for any topological space; items 2 and 3 work in an arbitrary separable metric space.)

• It would perhaps be even clearer to say that for each fixed vector $h\in\textrm{supp}\; v$, we have $P(h\in \overline{V})=1$. (This shows that the assumption that $H$ is separable is used, and how you avoid the pitfall from Pietro's comment above.) – Christian Remling Apr 22 '19 at 19:53
• @ChristianRemling: I edited my answer to emphasize where separability is essential. Let me know if anything is wrong. Thanks! – Mateusz Kwaśnicki Apr 23 '19 at 21:39
• Looks good to me. I didn't mean to suggest that such an extensive edit was needed, of course, just wanted to draw attention to where your answer differs from Pietro's comment. – Christian Remling Apr 23 '19 at 23:13

Another Try

We say a $$\mathcal{H}$$-valued random variable $$h$$ is a random vector if $$P(h \perp g)<1$$ for all $$g\in \mathcal{H}.$$

If $$h_1, h_2, \ldots$$ is a sequence independent identically distributed of random vectors, then, almost surely, the closed span of the $$h_i$$ is equal to $$\mathcal{H}.$$

First we will need a lemma.

Lemma 1 Let $$h$$ be a random vector. There is a countable subset $$A$$ of $$\mathcal{H}$$ such that the closed span of the elements of $$A$$ is equal to $$\mathcal{H}$$ and for every point $$a\in A,$$ $$P(h\in U)>0$$ for any neighborhood $$U$$ of $$a.$$

Proof For any subset $$A$$ such that for every point $$a\in A,$$ $$P(h\in U)>0$$ for any neighborhood $$U$$ of $$a,$$ and the closed span of the elements of $$A$$ is not equal to $$\mathcal{H},$$ we will show that we can grow $$A$$ by a single element which is not in closed span of the elements of $$A.$$ We can only do this a countable number of times because the Hilbert space dimension of $$\mathcal{H}$$ is countable. (Otherwise, via Gram-Schmidt, we could construct an uncountable orthonormal set by transfinite induction.)

Choose $$g$$ such that $$g \perp a$$ for all $$a\in A.$$ Now, $$P(h \perp g)<1.$$ So there must be a point $$b$$ such that $$P(h\in U) >0$$ for every neighborhood of $$b$$ and $$b$$ is not perpendicular to $$g,$$ therefore, $$b$$ is not in the span of the elements of $$A.$$ QED

Suppose $$h_1, h_2, \ldots$$ is a sequence independent identically distributed of random vectors. Let $$A$$ be as in Lemma 1. Index $$A$$ a a sequence $$a_n.$$ Let $$B_{m,n}$$ be a ball of radius $$1/m$$ centered at $$a_n$$ Almost surely, the sequence $$h_i$$ must visit $$B_{m,n}$$ infinitely often, as $$P(h_i\in B_{m,n})>0$$. Therefore $$A$$ is a subset of the closure of the values of the sequence. (We have essentially the fact that a random function $$f:\mathbb{N}\rightarrow \mathbb{N}^2$$ is surjective with infinite multiplicity.)

• The point is that you "keep going" by transfinite induction. The process must stop at some countable ordinal before $\omega_1$ as the space is countable dimensional. – J. E. Pascoe Apr 21 '19 at 19:59
• That is for each $\alpha < \omega_1$ there would be $A_\alpha$ such that if $\alpha < \beta,$ $A_\alpha^\perp \cap A_\beta \neq \{0\}.$ – J. E. Pascoe Apr 21 '19 at 20:06
• "Let $A$ be as in Lemma 1" ... "we may take $A$ to consist of only isolated points, as in a Polish space a countable closed set is equal to the closure of its isolated points." But Lemma 1 does not guarantee that $A$ is closed. – Iosif Pinelis Apr 21 '19 at 20:07
• That does seem to be a gap @IosifPinelis . Ideas for closing it? – J. E. Pascoe Apr 21 '19 at 20:10
• Could you elaborate a bit further on why the set $A$ in Lemma 1 is countable? Of course, there exist linearly independent sets in $H$ that are uncountable -- just choose your favourite Hamel basis of $H$. – Jochen Glueck Apr 21 '19 at 20:18


text:
# JEE Novice - (31)

Geometry Level 4

A straight line given by $$3x+2y=6$$ intersects $$x$$ axis and $$y$$ axis at $$A,B$$ respectively.If the distance between the incentre and orthocentre of $$\Delta OAB$$ can be expressed in the form $$\dfrac{a\sqrt{b}}{c+\sqrt{d}}$$ , find $$a+b+c+d$$ where $$b,d$$ are square free.

×


text:
# hermite_int_full_domain: Outputs integral of the orthonormal Hermite functions on the... In hermiter: Efficient Sequential and Batch Estimation of Univariate and Bivariate Probability Density Functions and Cumulative Distribution Functions along with Quantiles (Univariate) and Nonparametric Correlation (Bivariate)

 hermite_int_full_domain R Documentation

## Outputs integral of the orthonormal Hermite functions on the full domain

### Description

The method calculates \int_{-∞}^{∞} h_k(t) dt for k=0,…,N.

### Usage

hermite_int_full_domain(N)

### Arguments

 N An integer number.

### Value

A numeric matrix with N+1 rows and 1 columns.

### Author(s)

Michael Stephanou <michael.stephanou@gmail.com>

hermiter documentation built on Nov. 16, 2022, 1:08 a.m.


text:
# How to rate the similarity of aspect ratio of rectangles

Given two rectangles, and a maximum difference score, tell me whether the aspect ratio of these rectangles are considered similar or not?

For example a $$2\times4$$ rectangle [ratio $$0.5$$] and a $$2\times3$$ [ratio $$0.66$$] are considered similar, meaning that the difference of $$0.16$$ is still considered similar.

Also a $$4\times2$$ rectangle [ratio $$2$$] and a $$3\times2$$ rectangle [ratio $$1.5$$] are both considered similar and the difference is $$0.5$$ which is much larger than $$0.16$$ and yet is considered similar.

Another example of a $$3\times2.9$$ rectangle and a $$2.8\times3.1$$ rectangle are also considered similar.

How can one give one formula for calculating such a thing?

Further clarification

Let us represent a completely identical ratio with a score of $$0$$, meaning that there is no change between them

Let's represent a completely different ratio with a score of $$1$$

Let's say that the largest rectangle can be up to $$1000$$ width or up to $$1000$$ height

Example of the same ratio: $$2\times3$$ rectangle and $$4\times6$$ rectangle and the like

An example of a completely different ratio: a $$1\times1000$$ rectangle and a $$1000\times1$$ rectangle

There is no mathematical definition for the similarity of aspect ratios. You might consider the ratio of the aspect ratios as a measure of aspect ratio similarity. The closer to $$1$$ the more similar.


text:
# Simple-iteration method

A method for approximately solving a system of linear algebraic equations $Ax = b$ that can be transformed to the form $x = Bx + c$ and whose solution is looked for as the limit of a sequence $x ^ {k+} 1 = B x ^ {k} + c$, $k = 0 , 1 \dots$ where $x ^ {0}$ is an initial approximation. In order that the simple-iteration method converges for any initial approximation $x ^ {0}$ it is necessary and sufficient that all eigenvalues of $B$ are less than one in modulus; it is sufficient that some norm of $B$ is less than one. If in some norm, compatible with the norm of a vector $x$, $B$ satisfies $\| B \| \leq \rho < 1$, then the simple-iteration method converges at the rate of a geometric series and the estimate

$$\| x ^ {m} - x \| \leq \rho ^ {m} \| x ^ {0} - x \|$$

holds for its error.

In the case of a cubic, octahedral or spherical vector norm, the condition $\| B \| \leq \rho$ is fulfilled if

1) $\sum _ {j=} 1 ^ {n} | b _ {ij} | \leq \rho$, $i = 1 \dots n$;

2) $\sum _ {i=} 1 ^ {n} | b _ {ij} | \leq \rho$, $j = 1 \dots n$;

3) $\sum _ {i , j = 1 } ^ {n} b _ {ij} ^ {2} \leq \rho ^ {2}$.

The simplest version of the method corresponds to the case when one takes $I - A$, where $I$ is the identity matrix, for $B$. If all diagonal entries of $A$ are non-zero, then, choosing $b = D ^ {-} 1 ( D - A )$ and $c = D ^ {-} 1 b$, where $D$ is the diagonal matrix with as diagonal entries those of $A$, one obtains the Jacobi method or the method of simultaneous displacement.

A particular case of the simple-iteration method is the method with $B = I - \tau A$ and $c = \tau b$, where $\tau$ is an iteration parameter, chosen from the condition that the norm of $I - \tau A$ is minimal with respect to $\tau$. If $\gamma _ {1}$ and $\gamma _ {2}$ are the minimal and maximal eigenvalues of a symmetric positive-definite matrix $A$ and $\tau = 2 / ( \gamma _ {1} + \gamma _ {2} )$, then one has for the matrix $B$ in the spherical norm the estimate $\| B \| \leq \rho$, with $\rho = ( \gamma _ {2} - \gamma _ {1} ) / ( \gamma _ {2} + \gamma _ {1} ) < 1$.

For a system of non-linear algebraic equations

$$\phi _ {i} ( x) = 0 ,\ 1 \leq i \leq n ,\ x = ( x _ {1} \dots x _ {n} ) ,$$

the simple-iteration method has the form

$$x _ {i} ^ {k+} 1 = x _ {i} ^ {k} - \tau \phi _ {i} ( x ^ {k} ) ,\ 1 \leq i \leq n ,\ k \geq 0 .$$

The problem of choosing the iteration parameter $\tau$ is solved in dependence on the differentiability properties of the $\phi _ {i}$. Often it is subjected to the requirement that the method converges locally in a neighbourhood of a solution.

#### References

 [1] D.K. Faddeev, V.N. Faddeeva, "Computational methods of linear algebra" , Freeman (1963) (Translated from Russian) [2] I.S. Berezin, N.P. Zhidkov, "Computing methods" , Pergamon (1973) (Translated from Russian) [3] J.M. Ortega, W.C. Rheinboldt, "Iterative solution of non-linear equations in several variables" , Acad. Press (1970) [4] A.A. Samarskii, E.S. Nikolaev, "Numerical methods for grid equations" , 1–2 , Birkhäuser (1989) (Translated from Russian)
How to Cite This Entry:
Simple-iteration method. Encyclopedia of Mathematics. URL: http://encyclopediaofmath.org/index.php?title=Simple-iteration_method&oldid=48702
This article was adapted from an original article by E.S. Nikolaev (originator), which appeared in Encyclopedia of Mathematics - ISBN 1402006098. See original article


text:
# IEEE Transactions on Nuclear Science

## Issue 4  Part 1 • Aug. 2012

This issue contains several parts.Go to:  Part 2  | Part 3

## Filter Results

Displaying Results 1 - 25 of 83
• ### [Front cover]

Publication Year: 2012, Page(s): C1
| PDF (76 KB)
• ### IEEE Transactions on Nuclear Science publication information

Publication Year: 2012, Page(s): C2
| PDF (138 KB)

Publication Year: 2012, Page(s):689 - 693
| PDF (197 KB)
• ### Comments by the Chairmen

Publication Year: 2012, Page(s):694 - 695
| PDF (238 KB) | HTML
• ### Comments by the Editors

Publication Year: 2012, Page(s): 696
| PDF (22 KB) | HTML
• ### List of reviewers

Publication Year: 2012, Page(s):697 - 698
| PDF (83 KB)
• ### Comments by the Awards Chairmen

Publication Year: 2012, Page(s):699 - 700
| PDF (36 KB)
• ### Modeling the Effects of Hydrogen on the Mechanisms of Dose Rate Sensitivity

Publication Year: 2012, Page(s):701 - 706
Cited by:  Papers (9)
| | PDF (712 KB) | HTML

The effects of hydrogen on dose-rate sensitivity are simulated using a one-dimensional (1-D) model that incorporates the physical mechanisms contributing to dose-rate effects in the metal-oxide-semiconductor (MOS) system of gated lateral pnp (GLPNP) bipolar transistors. Calculations show that molecular hydrogen cracking at positively charged defects may be a key reaction relating hydrogen and dose... View full abstract»

• ### Influence of 60-MeV Proton-Irradiation on Standard and Strained n- and p-Channel MuGFETs

Publication Year: 2012, Page(s):707 - 713
Cited by:  Papers (5)
| | PDF (967 KB) | HTML

In this work the proton irradiation influence on Multiple Gate MOSFETs (MuGFETs) performance is investigated. This analysis was performed through basic and analog parameters considering four different splits (unstrained, uniaxial, biaxial, uniaxial+biaxial). Although the influence of radiation is more pronounced for p-channel devices, in pMuGFETs devices, the radiation promotes a higher immunity t... View full abstract»

• ### Geant4 Analysis of n-Si Nuclear Reactions From Different Sources of Neutrons and Its Implication on Soft-Error Rate

Publication Year: 2012, Page(s):714 - 722
Cited by:  Papers (14)
| | PDF (1879 KB) | HTML

This work examines nuclear events resulting from the interaction of atmospheric neutrons at ground level and different atmospheric-like sources with a silicon layer. Using extensive Geant4 simulations and in-depth data analysis, this study provides a detailed comparison between several facilities and natural environment in terms of nuclear processes, secondary ion production and fragment energy di... View full abstract»

• ### Modeling the Non-Uniform Distribution of Radiation-Induced Interface Traps

Publication Year: 2012, Page(s):723 - 727
Cited by:  Papers (7)  |  Patents (1)
| | PDF (1066 KB) | HTML

The incorporation of total ionizing dose (TID) effects into surface-potential-based compact models requires calculating the dependence of surface potential $({psi}_{s})$ on radiation-induced defect densities. This dependence is obtained through the introduction of the oxide trapped charge View full abstract»

• ### Effect of Carrier Transport in Oxides Surrounding Active Devices on SEU in 45 nm SOI SRAM

Publication Year: 2012, Page(s):728 - 734
Cited by:  Papers (1)
| | PDF (1080 KB) | HTML

Three-dimensional NanoTCAD simulations of a 45 nm SOI SRAM cell indicate that the contributions of charge generation and transport in the surrounding oxides have a potentially meaningful effect on the SEU response in nanotechnologies. View full abstract»

• ### Comparison of Charge Pumping and $1/f$ Noise in Irradiated Ge pMOSFETs

Publication Year: 2012, Page(s):735 - 741
Cited by:  Papers (1)
| | PDF (1194 KB) | HTML

Irradiated Ge pMOSFETs have been characterized via charge pumping $(I_{rm CP})$ and $1/f$ noise. The noise increases much more with irradiation than does $I_{rm CP}$ for devices with eight Si monolayer... View full abstract»

• ### Spatial, LET and Range Dependence of Enhanced Charge Collection by Single Ion Strike in 4H-SiC MESFETs

Publication Year: 2012, Page(s):742 - 748
Cited by:  Papers (2)
| | PDF (908 KB) | HTML

The spatial dependence of transient currents induced in 4H-SiC MESFETs by various ions with the same Linear Energy Transfer (LET) but different projected ranges is evaluated. The largest signal is observed when an ion strikes the Gate-to-Drain (G-D) recess. In addition, enhanced transient currents are detected when an ion strikes the gate and the drain. The charge enhancement mechanism is clarifie... View full abstract»

• ### Rediscovery of Single-Event Gate Rupture Mechanism in Power MOSFETs

Publication Year: 2012, Page(s):749 - 754
Cited by:  Papers (1)
| | PDF (528 KB) | HTML

The catastrophic failure mode caused by single-event gate rupture phenomenon observed in power MOSFETs still remains a critical issue for devices to be used in space radiation environments. Detailed analyses of the devices damaged by the phenomenon suggested a new possible mechanism. A preliminary model for the mechanism was also proposed. View full abstract»

• ### Radiation-Induced Oxide Charge in Low- and High-H$_{2}$ Environments

Publication Year: 2012, Page(s):755 - 759
Cited by:  Papers (4)
| | PDF (593 KB) | HTML

Electronic structure calculations and irradiation measurements are used to obtain insight into oxide trapped charge mechanisms in varying hydrogen ambients. Quantitative agreement between measured and simulated oxide and interface-trap charge densities is obtained over a wide range of H$_{2}$ concentrations by implementing first-pri... View full abstract»

• ### Influence of the Manufacturing Process on the Radiation Sensitivity of Fluorine-Doped Silica-Based Optical Fibers

Publication Year: 2012, Page(s):760 - 766
Cited by:  Papers (7)
| | PDF (979 KB) | HTML

In this work, we analyze the origins of the observed differences between the radiation sensitivities of fluorine-doped optical fibers made with different fabrication processes. We used several experimental techniques, coupling in situ radiation-induced absorption measurements with post mortem confocal microscopy luminescence measurements. Our data showed that the silica intrinsic defects are gener... View full abstract»

• ### Comparative Analysis of MIS Capacitance Structures With High-k Dielectrics Under Gamma, $^{16}$ O and p Radiation

Publication Year: 2012, Page(s):767 - 772
| | PDF (1042 KB) | HTML

MIS capacitance structures, with Hafnium Oxide, Alumina and nanolaminate as dielectrics were studied under gamma photons $^{60}{hbox{Co}}$, 25 MeV oxygen ions and 10 MeV protons radiation using capacitance-voltage (C-V) characterization. View full abstract»

• ### 3-D Numerical Simulation of Bipolar Amplification in Junctionless Double-Gate MOSFETs Under Heavy-Ion Irradiation

Publication Year: 2012, Page(s):773 - 780
Cited by:  Papers (14)
| | PDF (1080 KB) | HTML

The bipolar amplification and charge collection in Junctionless Double-Gate MOSFETs (JL-DGFET) submitted to heavy-ion irradiation are investigated. The transient response of JL-DGFET is compared to that of conventional devices operating in inversion-mode (IM-DGFET). We show that the bipolar amplification is higher in junctionless devices than in conventional inversion-mode devices mainly due to th... View full abstract»

• ### Accelerated Oxidation of Silicon Due to X-ray Irradiation

Publication Year: 2012, Page(s):781 - 785
Cited by:  Papers (7)
| | PDF (571 KB) | HTML

Enhanced rates of oxide growth have been observed on silicon upon exposure to 10-keV X-ray irradiation. Oxide thicknesses were determined using spectroscopic ellipsometry on irradiated and control samples, and confirmed via X-ray photoelectron spectroscopy. The oxidation rate varied with the radiation total dose and dose rate. The increased oxidation rate is attributed to the generation of ozone, ... View full abstract»

• ### Effect of Ion Energy on Power MOSFET's Oxide Reliability

Publication Year: 2012, Page(s):786 - 791
Cited by:  Papers (6)
| | PDF (626 KB) | HTML

Heavy ion-induced Power MOSFET's reliability degradation has been studied. Irradiations were realized at low and high energy with the same electronic stopping power at the surface of the die. For both energies, a decrease of the charge to breakdown has been observed after irradiation. However, an enhanced degradation of the oxide layer reliability is observed at low energy. The lifetime reduction ... View full abstract»

• ### Gamma Radiation Effects in Vertically Aligned Carbon Nanotubes

Publication Year: 2012, Page(s):792 - 796
Cited by:  Papers (5)
| | PDF (1750 KB) | HTML

This paper describes an experimental study of gamma radiation effects in low-density arrays of vertically aligned carbon nanotubes. These arrays are characterized by excellent anti-reflective and absorbing properties for wavelengths from UV to IR, which makes them an interesting option for stray light control in optical space applications. Gamma irradiation equivalent to an estimated surface lifet... View full abstract»

• ### Contribution to SER Prediction: A New Metric Based on RC Transient Simulations

Publication Year: 2012, Page(s):797 - 802
Cited by:  Papers (1)
| | PDF (921 KB) | HTML

This work focuses on speeding up simulation time of SEU systematic detection in a 90 nm SRAM cell. Simulations were run in order to validate a simplified approach based on the injection of a noise source current at the sensitive node of an analytical RC circuit. Moreover, a new SEU reliability metric, mandatory for reliability studies, is introduced. It is based on based on transient I–V si... View full abstract»

• ### Single-Event-Hardened CMOS Operational Amplifier Design

Publication Year: 2012, Page(s):803 - 810
Cited by:  Papers (9)
| | PDF (1473 KB) | HTML

Novel RHBD techniques are described that utilize charge sharing to mitigate single-event voltage transients in a folded cascode operational amplifier. These techniques are analyzed using a new layout aware single-event model with the Cadence Spectre circuit simulator. The techniques are applied in each of the three stages of an operational amplifier (op amp). The sensitive node active charge cance... View full abstract»

• ### Constrained Placement Methodology for Reducing SER Under Single-Event-Induced Charge Sharing Effects

Publication Year: 2012, Page(s):811 - 817
Cited by:  Papers (16)
| | PDF (648 KB) | HTML

This paper presents a methodology to reduce the impact of double faults in a circuit by constraining the placement of its standard cells. A fault-injection emulation platform is used to analyze the single-event-induced charge sharing effect in every pair of nodes. Based on the sensitivity of each pair, guidelines are set in a commercial standard cell placement by using constraints. Results show th... View full abstract»

## Aims & Scope

IEEE Transactions on Nuclear Science focuses on all aspects of the theory and applications of nuclear science and engineering, including instrumentation for the detection and measurement of ionizing radiation; particle accelerators and their controls; nuclear medicine and its application; effects of radiation on materials, components, and systems; reactor instrumentation and controls; and measurement of radiation in space.

Full Aims & Scope

## Meet Our Editors

Editor-in-Chief
Paul Dressendorfer
11509 Paseo del Oso NE
Albuquerque, NM  87111  USA


text:
#### Evaluation on effects of the thermal aging and irradiation of Austenitis stainless steel welds in PWRs = 가동 원전 오스테나이트계 스테인리스강 용접부의 열취화 및 조사취화 평가

Cited 0 time in Cited 0 time in
• Hit : 74
Austenitic stainless steel welds are widely used in pressurizer surge-line, weld-overlay cladding and reactor vessel internals of pressurized water reactors (PWRs). Around 5 - 20 vol. % of residual $\delta-ferrite$ formed during solidification process of welding procedure have some advantages on high strength, corrosion resistance and weldability. However, it was known that presence of $\delta-ferrite$ could be susceptible to thermal aging and neutron irradiation due to microstructural evolution. After long-term service at operating temperature and neutron irradiation condition, Cr-rich ($\alpha'$) phase and secondary precipitate like G-phase formed in the $\delta-ferrite$ deteriorate the integrity of austenitic stainless steel weld especially for loss of fracture toughness because of combined effect on thermal aging and neutron irradiation. In order to simulate long-term exposure at 320 ℃ and irradiation dose , the E308 weld was thermally aged at accelerated aging temperature, such as 400 ℃, for 20,000 h prior to irradiation. Proton irradiation at a dose rate of 1.1 × $10^{-5} dpa/s$ was conducted in MIBL at 360 ℃ to 10 dpa. After thermal aging at 400 ℃ for 20,000 h, phase separation into Cr-rich ($\alpha'$) and Fe-rich ($\alpha$) phase induced by spinodal decomposition was observed in the $\delta-ferrite$ while any secondary precipitate have not beed found. After proton irradiation up to 10 dpa, the degree of spinodal decomposition was enhanced and especially secondary G-phase precipitate was formed. Consequently, due to synergism between thermal aging and irradiation, strength of $\delta-ferrite$ was further increased than that of only thermally aged. The result of nano-mechanical tests have shown that irradiation can further accelerate the degradation of mechanical property of $\delta-ferrite$ coming from huge amounts of G-phase formation at pressurized water reactor environment.
Jang, Changheuiresearcher장창희researcher
Description
한국과학기술원 :원자력및양자공학과,
Publisher
한국과학기술원
Issue Date
2019
Identifier
325007
Language
eng
Description

학위논문(석사) - 한국과학기술원 : 원자력및양자공학과, 2019.2,[v, 59 p. :]

Keywords

austenitic stainless steel welds▼athermal aging▼aneutron irradidation▼aJ-R curve▼airradiation embrittlement; 오스테나이트계 스테인리스강 용접부▼a열취화▼aJ-R 시험▼a조사취화

URI
http://hdl.handle.net/10203/266576
http://library.kaist.ac.kr/search/detail/view.do?bibCtrlNo=843336&flag=dissertation
Appears in Collection
NE-Theses_Master(석사논문)
Files in This Item
There are no files associated with this item.


text:
+0

0
89
2

In a single-file queue of n people with distinct heights, define a blocker to be someone who is either taller than the person standing immediately behind them, or the last person in the queue. For example, suppose that Ashanti has height a, Blaine has height b, Charlie has height c, Dakota has height d, and Elia has height e, and that a $$n \ge 1$$ and k ≥ 0, let Q(n,k) be the number of ways that n people can queue up such that there are exactly k blockers.

(a) Show that $$Q(3,2)= 2 \cdot Q(2,2)+ 2 \cdot Q(2,1).$$

(b) Show that for n ≥ 2 and $$k \ge 1, Q(n,k)=k \cdot Q(n-1,k)+(n-k+1) \cdot Q(n-1,k-1).$$(You can assume that Q(1,1)=1, and that Q(n,0)=0 for all n.)

______________

Help would be appreciated; this site doesn't have any complete answers, and I'm not sure how to proceed. Perhaps for part (a) we could bash the numbers; any ideas?

Thank you in advance to those helping!

Apr 26, 2022

#1
-1

This is trivial if you use recursion.

Apr 26, 2022
#2
0

Hi, thanks for your response; however, this isn't supposed to be done using recursion. Would you care to show your method/explain it further?

Guest Apr 26, 2022