text:
# Imaginary Numbers and DFT

I'm not a math guy per se, but i am trying to understand the DFT. I get to the point where imaginary numbers are used with Euler's formula. What I don't understand is why we need an imaginary plane or "i" to begin with since we can plot any point on a real plane using sin and cos ? I am using this site and I get to the Euler's formula page before I crash. I am a Musician/Programmer trying to broaden my knowledge.

https://jackschaedler.github.io/circles-sines-signals/

I think that in the sense you are asking, you don't in fact "need" the complex plane. You can just use the real plane with pairs of coordinates. That's choosing to write $$(\cos \theta , \sin \theta)$$ instead of $$e^{i \theta} .$$ But if you do eschew the complex number formalism you miss out on some compact notation and some algebra that helps you understand how the DFT works. In particular, the identity $$e^{i (\sigma + \theta)} = e^{i \sigma} e^{i \theta}$$ is much more suggestive than the addition formulas for $$\sin$$ and $$\cos$$.

It's true that in the end you can use the real Fourier series or the complex one and they are equal. I can think of one argument (which I think is very important) why the complex Fourier series is important:

Consider some function $$f$$ for which you want to find the real Fourier series. To do that you have to compute all the coefficients $$a_i, b_i$$ in the series $$f(x) = a_0 + \Sigma_{i = 1}^{\infty}a_icos(ix) + b_isin(ix)$$: The coefficients can be calculated as integrals, for the $$a_k$$s: $$a_k = 2/T\int_{-T/2}^{T/2}f(t)cos(\frac{kt2\pi}{T})dt$$ where $$T$$ denotes the fundamental period of the function $$f$$. Sometimes these integrals are hard to calculate but it is easier to calculate the integrals which are needed for the complex Fourier series coefficients (using the residue theorem of complex analysis for example).

Of course you can compute the real coefficients once you found the complex coefficients. But the thing is that for the sake of computation the complex coefficients are often very useful.

Since you are a programmer: In order to understand libraries which implement the DFT you will certainly have to learn about the complex Fourier series.

Since you are a musician: The real Fourier series is more interesting here because it's coefficients are easily interpreted. The sound of a particular instrument is determined by its waveform. Its waveform is a periodic function for which you can find the real Fourier series. In this representation the coefficient tell you which overtones are the most relevant for this instrument. For example: You are playing an A($$440$$Hz) on your instrument, hence your fundamental period is $$1/440s$$. Now assume that all the $$a_k$$'s are zero and that $$b_2$$ is significantly larger than other coefficients in the real Fourier series representation for your waveform. This means that the period $$1/880s$$ is important meaning that the overtone at $$880$$Hz is strongly represented in your instrument.

You don't need to use the complex Fourier series for many applications. The Wikipedia article Discrete cosine transorm describes some of them as follows:

DCTs are important to numerous applications in science and engineering, from lossy compression of audio (e.g. MP3) and images (e.g. JPEG) (where small high-frequency components can be discarded), to spectral methods for the numerical solution of partial differential equations

The article goes on to state that among the variants of the transform are

the discrete sine transform (DST), which is equivalent to a DFT of real and odd functions

Later on in more detail it states

The obvious distinction between a DCT and a DFT is that the former uses only cosine functions, while the latter uses both cosines and sines (in the form of complex exponentials).

You can read the article for more technical details and pointers to references for more information. The summary is that the complex Fourier series is elegant, but for many real world applications, there are clear advantages to using DCT. This does not negate the great importance of the DFT for theoretical purposes which is based on the use of Euler's formula to model simple harmonic motion. In other words, simple harmonic motion is closely related to circular motion using the important identity $$\, \sin(t)^2 + \cos(t)^2 = 1\,$$ and that $$\, \frac{d}{dt}\sin(t) = \cos(t), \: \frac{d}{dt} \cos(t) = -\sin(t).\,$$ The fact that these are naturally modeled by complex numbers is one of the most striking results in all of mathematics.


text:
I've implemented the Nelder-Mead optimization algorithm in C++. I found this document to be a very good explanation of how the algorithm works, but I'll try my best to explain it (with an example) here so it's clear for everyone.

The problem I'm trying to optimize is very simple : I'm trying to find the best a,b,c coefficients for a quadratic function to fit points on a quadratic function of form f(x) = 1.5*x^2 + 3.2*x + 7. So, of course, the optimal parameters I'm looking for are a=1.5, b=3.2, c=7.

In order to compute the error between the coefficients that I have and the expected coefficients, I use the following function :

#include "Eigen/Dense"
static double quadraticEquationError(Array<double, 1, 3> parameters)
{
double searchedParamters[] = {1.5, 3.2, 7};

auto x = Array<double, 100, 1>().setLinSpaced(100, 0, 90);
auto expected = searchedParamters[0] * x*x + searchedParamters[1] * x + searchedParamters[2];
auto result = parameters(0, 0) * x * x + parameters(0, 1) * x + parameters(0, 2);

auto diff = result.matrix() - expected.matrix();

return diff.norm();
}

It's very simple, but it's also a way for me to check if my algorithm works. The goal of this error function is to return 0, where I would have found my optimal parameters.

The Nelder-Mead algorithm takes as input an error function and some initial parameters (in my case, a=0, b=0, c=0). Then, it generates a simplex of N+1 dimensions (here, N = 3 since we are looking for 3 parameters) where each vertex of the simplex represents a new parameters set. After that, depending on the values of the worst, second worst and best sets, it replaces the worst point by a new value.

There are five operations that can be used : reflection, expansion, inside contraction, outside contraction and shrinking. The following figure shows the result of these operations for a 2D simplex (of 3 points).

In order to figure out which operation to use, we can look at the following diagram :

Progressively, the error will go down until a convergence criterion is reached. In my case, I decided to stop when the best error is under a certain value.

So, I hope this is pretty clear as to how the algorithm works, here's my code :

main.cpp

#include "library.h"
#include <iostream>
#include "Eigen/Dense"

using namespace Optimization;

static double quadraticEquationError(Array<double, 1, 3> parameters)
{
double searchedParamters[] = {1.5, 3.2, 7};

auto x = Array<double, 100, 1>().setLinSpaced(100, 0, 90);
auto expected = searchedParamters[0] * x*x + searchedParamters[1] * x + searchedParamters[2];
auto result = parameters(0, 0) * x * x + parameters(0, 1) * x + parameters(0, 2);

auto diff = result.matrix() - expected.matrix();

return diff.norm();
}

int main() {
Array<double, 1, 3> initialParameters;
initialParameters.row(0) << 0.0, 0.0, 0.0;

return 0;
}

#include "Eigen/Dense"
#include <vector>
#include <limits>
#include <cstring>
#include <iostream>

using namespace Eigen;

namespace Optimization
{
template<int nDims>
public:
typedef Array<double, 1, nDims> FunctionParameters;
typedef std::function<double(FunctionParameters)> ErrorFunction;

double minError, double initialEdgeLength,
double shrinkCoeff = 1, double contractionCoeff = 0.5,
double reflectionCoeff = 1, double expansionCoeff = 1)
: errorFunction(errorFunction), minError(minError),
shrinkCoeff(shrinkCoeff), contractionCoeff(contractionCoeff),
reflectionCoeff(reflectionCoeff), expansionCoeff(expansionCoeff),
worstValueId(-1), secondWorstValueId(-1), bestValueId(-1)
{
this->errors = std::vector(nDims + 1, std::numeric_limits<double>::max());

const double b = initialEdgeLength / (nDims * SQRT2) * (sqrt(nDims + 1) - 1);
const double a = initialEdgeLength / SQRT2;

this->values = initial.replicate(nDims + 1, 1);

for (int i = 0; i < nDims; i++)
{
FunctionParameters simplexRow;
simplexRow.setConstant(b);
simplexRow(0, i) = a;
simplexRow += initial;

this->values.row(i+1) = simplexRow;
}
}

void optimize()
{
for (int i = 0; i < nDims+1; i++)
{
this->errors.at(i) = this->errorFunction(this->values.row(i));
}

this->invalidateIdsCache();

while (this->errors.at(this->bestValueId) > this->minError)
{
step();
auto bestError = this->errorFunction(this->best());
auto worstError = this->errorFunction(this->worst());

std::cout << "Best error " << std::to_string(bestError) << " with : " << this->best();
std::cout << " Worst error " << std::to_string(worstError) << " with : " << this->worst();
std::cout << '\n';
}
}

void step()
{
auto meanWithoutWorst = this->getMeanWithoutWorst();

auto reflectionOfWorst = this->getReflectionOfWorst(meanWithoutWorst);
auto reflectionError = this->errorFunction(reflectionOfWorst);

FunctionParameters newValue = reflectionOfWorst;

bool shrink = false;

if (reflectionError < this->errors.at(this->bestValueId))
{
auto expansionValue = this->expansion(meanWithoutWorst, reflectionOfWorst);
double expansionError = this->errorFunction(expansionValue);

if (expansionError < this->errors.at(this->bestValueId))
{
newValue = expansionValue;
}
}
else if (reflectionError > this->errors.at(this->worstValueId))
{
newValue = this->insideContraction(meanWithoutWorst);

if (newError > this->errors.at(this->worstValueId)) { shrink = true; }
}
else if (reflectionError > this->errors.at(this->secondWorstValueId))
{
newValue = this->outsideContraction(meanWithoutWorst);

if (newError > reflectionError) { shrink = true; }
}
else
{
newValue = reflectionOfWorst;
}

if (shrink)
{
this->shrink();
this->invalidateIdsCache();
return;
}

this->values.row(this->worstValueId) = newValue;
this->invalidateIdsCache();
}

inline FunctionParameters worst() { return this->values.row(this->worstValueId); }
inline FunctionParameters best() { return this->values.row(this->bestValueId); }

private:

void shrink()
{
auto bestVertex = this->values.row(this->bestValueId);

for (int i = 0; i < nDims + 1; i++)
{
if (i == this->bestValueId) { continue; }

this->values.row(i) = bestVertex + this->shrinkCoeff * (this->values.row(i) - bestVertex);
this->errors.at(i) = this->errorFunction(this->values.row(i));
}
}

inline FunctionParameters expansion(FunctionParameters meanWithoutWorst, FunctionParameters reflection)
{
return reflection + this->expansionCoeff * (reflection - meanWithoutWorst);
}

inline FunctionParameters insideContraction(FunctionParameters meanWithoutWorst)
{
return meanWithoutWorst - this->contractionCoeff * (meanWithoutWorst - this->worst());
}

inline FunctionParameters outsideContraction(FunctionParameters meanWithoutWorst)
{
return meanWithoutWorst + this->contractionCoeff * (meanWithoutWorst - this->worst());
}

FunctionParameters getReflectionOfWorst(FunctionParameters meanWithoutWorst)
{
return meanWithoutWorst + this->reflectionCoeff * (meanWithoutWorst - this->worst());
}

FunctionParameters getMeanWithoutWorst()
{
FunctionParameters mean(0);
for (int i = 0; i < nDims + 1; i++)
{
if (i == this->worstValueId) { continue; }

mean += this->values.row(i);
}

// Not divided by nDims+1 because there's one ignored value.
mean /= nDims;

return mean;
}

void invalidateIdsCache() {
double worstError = std::numeric_limits<double>::min();
int worstId = -1;
double secondWorstError = std::numeric_limits<double>::max();
int secondWorstId = -1;
double bestError = std::numeric_limits<double>::max();
int bestId = -1;

for (int i = 0; i < nDims + 1; i++)
{
auto error = this->errors.at(i);

if (error > worstError)
{
secondWorstError = worstError;
secondWorstId = worstId;
worstError = error;
worstId = i;
}
else if (error > secondWorstError)
{
secondWorstError = error;
secondWorstId = i;
}

if (error < bestError)
{
bestError = error;
bestId = i;
}
}

// If we deal with a problem in 1D, it won't be set.
if (secondWorstId == -1)
{
secondWorstId = worstId;
}

this->bestValueId = bestId;
this->worstValueId = worstId;
this->secondWorstValueId = secondWorstId;
}

ErrorFunction errorFunction;
Array<double, nDims + 1, nDims> values;
std::vector<double> errors;

int worstValueId;
int secondWorstValueId;
int bestValueId;

double minError;

double shrinkCoeff;
double expansionCoeff;
double contractionCoeff;
double reflectionCoeff;

const double SQRT2 = sqrt(2);
};
}

Now, the thing is I don't know much about C++, but I know a lot about C# so I think my code structure is fine. I tried to use templates and typedefs to make code more readable/generic, but I'd like my C++ skills to improve a lot since I'm starting a new job soon that will require me to write C++ (my employer knows I'm not super strong in C++ but I'd like to get up to speed).

I'm looking for parts of my code that aren't using the full potential of C++, although I'm open to any part of the code reviewed.

# Usage

## Documentation

I'll consider the text contents of the post as documentation and it looks good.

## Consumers

It is important to have a straightforward consumption process. Usually it is having a CMakeLists script that searches for required dependencies and sets up required (!) compilation flags. The stuff like optimization levels and warnings can go into toolchain file.

## Usage flow

At the moment it looks like a toy program. What is the usecase for the current implementation? I would expect something like the following:

template <typename T, std::size_t Dimensions>
using coefficients_t = Array<T, 1, Dimensions>;

template <typename T, std::size_t Dimensions>
struct optimization_state_t {
coefficients_t<T, Dimensions> coefficients;
std::function<double(coefficients_t<T, Dimensions>)> error_function;
/*errors and what other internal state is needed */
};

template <typename T, std::size_t Dimensions>
optimization_state_t make_starting_state(coefficients_t<T, Dimensions> start = {}, std::function<double(coefficients_t<T, Dimensions>)> error_function = default_error_function);

template <typename T, std::size_t Dimensions>
void perform_step(optimization_state_t<T, Dimensions>& state);

template <typename T, std::size_t Dimensions>
optimization_state_t<T, Dimensions> optimize_until(coefficients_t<T, Dimensions> coefficients);

This way, my flow will be

decide on initial parameters -> optimize one step and see if the convergence speed is good -> try again or abort -> feed results somewhere else

or

decide on initial parameters -> optimize until some condition -> feed results somewhere else

I believe interfaces should support desired usage flow instead of just providing bits and pieces that together form a crooked road.

Note that error function is part of the optimization state, as the error values will lose meaning in incompatible error functions. As Scott Meyers said, interfaces should be "easy to use correctly and hard to use incorrectly".

# Details

## Proper types

I'm not an experienced user of Eigen, I believe Vector3d is a better fit due to the conversion into matrix in the error function and generally solving a system of equations. There is also a library called Blaze, which from what I remember makes better use of intel MKL (needless to say that on AMD CPUs performance will be subpar).

## Use standard constants

There are well defined constants since C++20. Computing it from scratch using library features might change the value from release to release.

## Consider the cost of accepting by value

To accept by value, a new instance of the type has to be made. If the type does dynamic memory allocation (new, malloc, etc) and there is no real need to copy/move/construct a new object, then it is better to avoid paying the performance cost of additional memory allocation by accepting by reference. On the other hand, if Array<double, 1, 3> does not do dynamic memory allocation, it is better to accept it by value since it is small and will not cause stackoverflow

## Do not use checking version of subscripts

Subscript operators that do bounds checking prevent SIMD instruction generation for compilers. It is better to just do plain assert and call non-checking version.

## Do not import names into current scope without the need

using namespace something is usually a bad idea inside a header. Aside from obvious name clashes, there is ADL. The bottom line is that if Eigen decides to add similar function to yours and the call is not qualified with a namespace, two erroneous outcomes will happen. First and the best one is that ambigious call error will be issued because both functions are at equal rank in overload resolution. The worst case is that the compiler will silently select unintended function because the programmer didn't know about this (and most people don't know about ADL).

# Style

## Usage of this

this is usually frowned upon unless there are some clever mixin tricks need to be done. I do not really know a strong reason why it is frowned upon, but arguments thrown around are that it makes it easy to confuse argument with member variable, that this can be a bit weird with lambdas and so on.

## Naming

The standard uses snake_case for almost everything except template parameter names which are Pascal case.

There is also Google style guide, MISRA and the others. I believe it is important to just pick a well-known one and stick with it, because clang-format configs are already written for them.

 double searchedParamters[] = {1.5, 3.2, 7};

You don't change the values within the function, so it ought to be const and there's no reason to create a local (stack-based) variable and copy the three elements into it on every invocation. Make this constexpr.

Now, the thing is I don't know much about C++, but I know a lot about C#

So, I can guess that you'll be completely overlooking the existence of const, and will probably not consider that C++ has value semantics.

static double quadraticEquationError(Array<double, 1, 3> parameters)

I don't know what Array is (your code posted is incomplete), but I thing you should not be passing it by value.

typedef std::function<double(FunctionParameters)> ErrorFunction;
Use using ErrorFunction = std::function<double(FunctionParameters)>; rather than the old typedef syntax.

        NelderMead(ErrorFunction errorFunction, FunctionParameters initial,
double minError, double initialEdgeLength,
double shrinkCoeff = 1, double contractionCoeff = 0.5,
double reflectionCoeff = 1, double expansionCoeff = 1)
: errorFunction(errorFunction), minError(minError),
shrinkCoeff(shrinkCoeff), contractionCoeff(contractionCoeff),
reflectionCoeff(reflectionCoeff), expansionCoeff(expansionCoeff),
worstValueId(-1), secondWorstValueId(-1), bestValueId(-1)

The final three members that have constant initializers should use inline direct initializers in the class, and then you don't need to name them here.
Uniform Initialization is to use curly braces, but parens still work.

this->errors = std::vector(nDims + 1, std::numeric_limits<double>::max()); etc.
Don't use this-> to access members.

this->errors.at(i) = this->errorFunction(this->values.row(i));
really, don't do that!

The class member: const double SQRT2 = sqrt(2);
should (at least) be static, as it's always the same so there is no reason for each instance to have its own. Having a const member will disable the autogenerated assignment operator.

inline FunctionParameters expansion(FunctionParameters meanWithoutWorst, FunctionParameters reflection)
Don't use the inline keyword for functions defined inside the class. I notice you only use that for some of them, anyway. Such functions are automatically inline. (inline means they can be in the .h file, BTW. Nothing more.)

Your sample main calls optimize but doesn't read out any results. As a synopsis of how the class is to be used, that's not very useful since that essentially doesn't do anything.

• Array is a class from the Eigen library. And how should I access members other than using this? I'd like to understand why I shouldn't do X! :) Jan 14 at 15:56
• @IEatBagels Members are in scope within member functions. Just delete the characters this-> and use the member name directly. Isn't that the same in C#? E.g. errors.at(i) = errorFunction(values.row(i)); Jan 14 at 16:16
• When I compiled, Array was resolved as Eigen::Array due to the using namespace. But agreed, it's better to be explicit. Jan 15 at 10:07

# Naming things

FunctionParameters is a bad name. There are many functions and many parameters. A much better name is Coordinates, as the type represents the coordinates of a vertex. You could even consider naming it Vertex, but consider that not every coordinate is an actual vertex (like meanWithoutWorst).

Use verbs for function names, and prefer nouns for variable names. Also be consistent with naming. For example, instead of expansion(), write getExpansion(), so it's a verb and matches getMeanWithoutWorst().

Also avoid having both variables and functions with the same name, like shrink. I would recommend keeping shrink() as the function name, but renaming the variable to doShrink. It's hard to follow the noun/verb split here.

# Rethink how you are templating your code

Instead of templating your algorithm on the number of dimensions, and then forcing coordinates of vertices to be Array<double, 1, Ndims>, consider that the Nelder-Mead algorithm only needs some basic vector algebra and doesn't really care about the type of things, apart from the dimension. You could template it on the vertex coordinate type instead, and have it derive the dimensions of the coordinates from that:

template <typename Coordinates>
static constexpr std::size_t nDims = Coordinates::ColsAtCompileTime + 1;
...

This still assumes Coordinates is a type that has a ColsAtCompileTime member value, but with a bit of work you could make this more generic. You can now use arrays of float or even std::complex for the vertex coordinates.

Also consider that you might not want to constrain the return type of the error function. Consider making the type of the error function a template parameter as well. You can use concepts to restrict it to only function types that take Coordinates as arguments:

template <typename Coordinates, typename ErrorFunction>
requires std::invocable<ErrorFunction, Coordinates>
static constexpr std::size_t nVertices = Coordinates::ColsAtCompileTime + 1;
using Error = std::invoke_result_t<ErrorFunction, Coordinates>;
...
public:
...
private:
ErrorFunction errorFunction;
std::array<Coordinate, nVertices> values;
std::array<Error, nVertices> errors;
...
};

With this, class template argument deduction (CTAD) works and you no longer have to specify any template arguments when creating an instance of the NelderMead class. It also shows you how to avoid hardcoding double in your code.

# Consider keeping values sorted by error

Instead of keeping track of bestId, worstId and secondWorstId, and introducing indirection in a lot of code, it might be better to keep the coordinates of the vertices sorted by their error. This is rather easy to do if you create a new struct to hold both a vertex's coordinates and error, and use that instead of having separate arrays/vectors for values and errors:

struct Vertex {
Coordinates coordinates;
ErrorType error;
};

std::array<Vertex, nVertices> vertices;

Sorting it on the error is especially easy with C++20's std::ranges::sort():

std::ranges::sort(vertices, {}, &Vertex::error);

This way, the best vertex is always vertices[0], second best vertices[1], and so on. Basically, you should sort whenever you called invalidateIdsCache(). Sorting is $$\O(N \log N)\$$, whereas your invalidateIdsCache() is $$\O(N)\$$, so when your coordinates have a large number of dimensions, it might not be worth it, but at least up to a few tens of dimensions it should be fine.

# Potential infinite loop

The Nelder-Mead algorithm doesn't guarantee that the error of the best vertex of the simplex will decrease until zero. Like many other optimization algorithms, it can end up in a local minimum, which might be higher than the minError you provide it. Consider adding a maxSteps parameter to limit the number of steps done, or some other heuristic to detect that the error is not improving significantly anymore.

• Correct me if I'm wrong, but isn't the use of a member-pointer quite inefficient? Jan 15 at 2:28
• @Deduplicator I'm assuming it will be inlined? Jan 15 at 7:25
• True, if it is all inlined, and considering it is the only time it is called with that signature in the program, it is virtually guaranteed, it will be fine. Jan 15 at 7:42
• @G.Sliepen all of the new range-based algorithms take this "projection" argument, which saves code. Since it's now part of the standard, I would think the compiler's optimizers would be sure to handle these, at least in the specific form that the standard library has them. Jan 18 at 14:51
• @JDługosz Yes. In fact, if you'd pass a lambda or any other invocable then without inlining it would be very inefficient as well. Jan 18 at 20:04

main.cpp includes "library.h" but doesn't seem to use anything from it (just as well, as you've not presented it for review). It also includes <iostream> and doesn't use it (though see Don't mix I/O with calculation below).

optimization/NelderMead.h includes <cstring> and doesn't appear to use it. It does however need <string> instead, for std::to_string() - though that is unnecessary, as shown below. It's also missing <cmath>, for std::sqrt().

# Don't using namespace in headers

using namespace Eigen;

That causes every translation unit that includes the header to have its global namespace filled with Eigen identifiers. That's harmful, and in the worst cases can cause bugs to appear when a newer version of the library is used. So don't import namespaces at global scope.

There's also using namespace Optimization; in the implementation file. That's less harmful, though I still recommend against it. Move the import inside main(), or (better) drop it altogether. After all, it adds more text and cognitive load than it saves!

# Choose template type wisely

We have

template<int nDims>

What happens if nDims is negative? Consider using an unsigned type for the dimension count.

# Write initialiser lists in the order of execution

Object members are initialised in the order they are declared in the class. It helps readers if you use the same order in the constructor's initialiser list; this is particularly important when members depend on earlier members or the initialisers have side-effects.

We should be initialising all members (including errors and values); we can initialise members that don't depend on constructor arguments in their definitions.

That gives us:

    std::vector<double> errors{std::vector(nDims + 1, std::numeric_limits<double>::max())};
int worstValueId = -1;
int secondWorstValueId = -1;
int bestValueId = -1;

    NelderMead(ErrorFunction errorFunction, FunctionParameters initial,
double minError, double initialEdgeLength,
double shrinkCoeff = 1, double contractionCoeff = 0.5,
double reflectionCoeff = 1, double expansionCoeff = 1)
: errorFunction{errorFunction},
values{initial.replicate(nDims + 1, 1)},
minError{minError},
shrinkCoeff{shrinkCoeff},
expansionCoeff{expansionCoeff},
contractionCoeff{contractionCoeff},
reflectionCoeff{reflectionCoeff}
{

# We don't need √2 in every object

We have a constant SQRT2 which doesn't need to be an instance member. Make it a static const and change its name so that it doesn't look like a preprocessor macro (we conventionally use all-caps for macro definitions, to identify hard-to-debug text substitution).

Given that it's only used in the constructor, we can reduce its scope significantly:

    static const auto root2 = std::sqrt(2);
auto const a = initialEdgeLength / root2;
auto const b = initialEdgeLength / (nDims * root2) * (std::sqrt(nDims + 1) - 1);

At that point, it's probably better to do away with the named constant entirely, and write std::sqrt(2) in both those places - any decent compiler will use constant values there, and it's just as clear to the reader.

# Don't scatter this-> all over the code

In C++, when an object member is referenced, it means the one belonging to this, unless shadowed by a local variable (just don't do that). So cluttered statements such as

            this->errors.at(i) = this->errorFunction(this->values.row(i));

can be much simpler and therefore clearer:

            errors.at(i) = errorFunction(values.row(i));

There's no need to hand-code this loop:

        for (int i = 0; i < nDims+1; i++)
{
this->errors.at(i) = this->errorFunction(this->values.row(i));
}

Instead, include <algorithm> and <ranges> (C++20) and use the standard-library transform():

        std::ranges::transform(values.rowwise(), errors.begin(), errorFunction);

# Don't mix I/O with calculation

NelderMead::optimize() runs the algorithm and prints the result. This makes it inflexible, because calling code can't do anything different with the result. Separate the actions, so that the function either returns all the results (perhaps as an aggregate type) or sets appropriate member variables so that the optimum can be accessed via the object.

Removing printing from this function then removes the need to include <iostream> from the header (and likely transfers it to main.cpp, if we want to print the results from there).

While looking at that output operation, I see you're unnecessarily constructing strings for types that have << operator defined. There's no need to do that:

            std::cout << "Best error " << bestError << " with : " << best()
<< " Worst error " << worstError << " with : " << worst()
<< '\n';

# Use const appropriately

worst() and best() don't need to change the state of the object, so should be declared const:

    FunctionParameters worst() const { return values.row(worstValueId); }
FunctionParameters best() const { return values.row(bestValueId); }

I removed inline from these declarations, because functions declared in-class are automatically inline (i.e. can be defined in multiple translation units).

The same applies to a number of other functions.

# Algorithms again

Instead of the logic in getMeanWithoutWorst() to sum all but one value, it is simpler to add all the values using std::accumulate (in <numeric>), then subtract the one we want to ignore:

    FunctionParameters getMeanWithoutWorst() const
{
auto const& rows = values.rowwise();
auto const total = std::accumulate(rows.begin(), rows.end(), FunctionParameters{})
- values.row(worstValueId);

}

A C++23 standard library should get std::ranges::accumulate() which would simplify further:

    FunctionParameters getMeanWithoutWorst() const
{
auto const total = std::ranges::accumulate(values.rowwise(), FunctionParameters{})
- values.row(worstValueId);

searchedParamters should be searchedParameters.
• std::sqrt appears to not be constexpr as of the C++20 standard. Jan 19 at 19:52
• Your use of "inlinable" might give the impression that such functions are nominated to be inlined. We should teach clearly that inline means "it can go in a header", and the historical meaning is no longer applicable. The compiler inlines whatever it wants, whether declared inline or not. Jan 19 at 19:56
• Oops, fell foul of my library's std::sqrt() going further than standard guarantees. I've improved the wording around inline, too. Jan 19 at 20:04
• FWIW, godbolt.org/z/jrn1h783c writing static const... inside a function, gcc will store it at compile time and multiplying uses a single instruction. msvc will also store it at compile time, but instead of just multipying by root2 in memory, it shuffles registers, loads the value from memory into another register, then multiplies. clang I'm sad to say does compute the value at compile time, but keeps the full machinery to initialize the variable on the first call. Jan 20 at 16:39

format: plain


text:
Version 1.5
$\newcommand{\lyxlock}{}$
Subsection 7.1.1: Fixed Displacement Up Section 7.1: Boundary Conditions Subsection 7.1.3: Fixed Shell Rotation

### 7.1.2 Fixed Shell Displacement

Certain shell element formulations used in FEBio have two sets of displacement degrees of freedom. One for the front-face of the shell, and one for the back face. The usual displacement degrees of freedom are used for the front-face of the shell, and the shell displacement degrees of freedom are used for the back-face of the shell. This boundary condition allows you to fix the back-face shell degrees of freedom.
Subsection 7.1.1: Fixed Displacement Up Section 7.1: Boundary Conditions Subsection 7.1.3: Fixed Shell Rotation


text:
# Is there a name for this algorithm for 2d shape collision?

A while ago I was considering my possibilities for implementing collisions between 2D shapes, and after much googling (which mostly turned up SAT), I came up with this:

let A, B = two shapes
let V = Speed of A relative to B

for a:point in vertices of A
for b:segment in segments of B
if segment(a, a+V) intersects b
return intersection point
return false

Now, I'm far from the smartest person out there, so I am 99.8% sure someone must have come up with that before me and, what's more important, since I didn't find anything describing this method, I assume it's a bad one.

Can anybody identify this algorithm and/or point me towards any literature on it explaining if, when and why it should be used?

I've failed to see any significant problems at first sight; all the complex calculations can be hoisted outside of the loops and then it's just matrix multiplications within both loops, and I couldn't find any algorithm that makes do without any nested loops over two sets of elements (be they vertices, segments or normals)

Feel free to trash-talk my algo all you like, I never expected it to be the second greatest discovery right after general relativity or anything ;)

PS: Sorry for the poor title, but I can hardly explain the entire algo in one line :)

• This algorithm will return false if A is completely within B or vice versa, by the way. – Philipp Jul 9 at 16:55
• @Philipp of course, thanks for pointing that out, but for B to get inside A it must have collided at some point, so for me that's acceptable. It also allows for a shape to be "trapped" within another, like a bottle, for example. – DarkWiiPlayer Jul 9 at 16:56
• You might like to search for "swept" or "continuous" collision detection. – Jay Jul 14 at 6:33


text:
# In the context of psychosomatic medicine, Flanders Dunbar and Franz Alexander maintained that conflicts produce anxiety, which becomes unconscious and takes a physiological toll on the body via the__________--

###### Question:

In the context of psychosomatic medicine, Flanders Dunbar and Franz Alexander maintained that conflicts produce anxiety, which becomes unconscious and takes a physiological toll on the body via the__________--

### The following qualities refer to a good ___________ ; golden brown crust, uniform attractive edge, blistery pebbly surface. cookie cake pastry

The following qualities refer to a good ___________ ; golden brown crust, uniform attractive edge, blistery pebbly surface. cookie cake pastry...

### Identify (y-3)^2=8(x+1)

Identify (y-3)^2=8(x+1)...

### Which sentence is written correctly? Es su pastel Es sus pastel Son sus pasteles

Which sentence is written correctly? Es su pastel Es sus pastel Son sus pasteles...

### 10.Check the 2 boxes that DO NOT FIT with the other examples of how Spanish exploration and colonization impacted the Americas. -Population Growth -Spread of Christianity -Spread of Disease -Increased Power For Native Tribes -Spanish Cultural Influences -European Technology -Introduced New Plants and Animals -Mestizo Culture Created -Ending of Slavery in the Americas

10.Check the 2 boxes that DO NOT FIT with the other examples of how Spanish exploration and colonization impacted the Americas. -Population Growth -Spread of Christianity -Spread of Disease -Increased Power For Native Tribes -Spanish Cultural Influences -European Technology -Introduced New Plants an...

### In this assignment, you will implement an online banking system. Users can sign-up with the system, log in to the system, change their password, and delete their account. They can also update their bank account balance and transfer money to another userâ€™s bank account. Youâ€™ll implement functions related to File I/O and dictionaries. The first two functions require you to import files and create dictionaries. User information will be imported from the "users.txt" file and account information

In this assignment, you will implement an online banking system. Users can sign-up with the system, log in to the system, change their password, and delete their account. They can also update their bank account balance and transfer money to another userâ€™s bank account. Youâ€™ll implement...

### Larry filled a milk carton completely with water and tightly secured the cap. He then placed it in the freezer. When he looked at the container two days later, he saw that the frozen water had expanded, breaking the container. What type of weathering is modeled by Larry’s experiment? A. chemical weathering B. biological chemical weathering C. physical weathering D. biological physical weathering

Larry filled a milk carton completely with water and tightly secured the cap. He then placed it in the freezer. When he looked at the container two days later, he saw that the frozen water had expanded, breaking the container. What type of weathering is modeled by Larry’s experiment? A. chemical w...

### What attracted European trading companies to India?

What attracted European trading companies to India?...

### _____ position refers to the baby's position in the uterus that causes the buttocks to be the first part to emerge from the vagina.

_____ position refers to the baby's position in the uterus that causes the buttocks to be the first part to emerge from the vagina....

### Need help please i dont get it

Need help please i dont get it...

### To make sure your reader remembers the point of your essay, which of the following is most important? A. Choosing an interesting subject B. The introduction C. The order in which you present your ideas D. The ending

To make sure your reader remembers the point of your essay, which of the following is most important? A. Choosing an interesting subject B. The introduction C. The order in which you present your ideas D. The ending...


text:
In [1]:

import numpy as np
import fst
import matplotlib.pyplot as P
%matplotlib inline

# Intro to HMMs

Hidden Markov Models (HMMs) are a well known mathematical construct for modeling sequences of events using probabilistic reasoning. Its advantage is that it's very intuitive and not too hard to implement. You can probably find much more information on places like Wikipedia and dozens of other sites you can google (e.g. this one is one of my favorites).

For speech recognition, the go to paper is an old one by Rabiner which is written in a very simple and easy to understand way with examples specifically tuned to its application in speech recognition. If this is your first time hearing about it, you have to go and read that paper at once. This notebook will simply go through some of the major points written there and show their implementation in Python.

## HMM algorithms

The Rabiner paper mentions 3 problems that are solved using HMMs:

1. Find the joint probability of a sequence of states given an HMM model - solved using the Forward algorithm
2. Finding the most likely sequence of hidden states given the probability above - solved using the Viterbi algorithm
3. Fine-tuning the parameters of the model by observing the likelihood of each state transition at each point in time - the Baum-Welch algorithm based around the Forward-Backward mechanic

As a reminder, the HMM is defined by the following parameters:

• observation probability - classically is the output of a Gaussian Mixture Model (GMM), but we want to use the output of the ANN (e.g. MLP) here instead. For consistency with the Rabiner paper we will instead use a discrete model here, which means that observations will simply be some descreete events and the model will gives us a simple value of how probable this event is given a particular state (from a lookup table)
• hidden states - in our case, these are the units we try to recognize: e.g. phonemes, but also possibly things like words, which are a combination of phonemes (models can be combined in several layers - more on that later)
• transitions - as the name suggests, transitions from one state to another - these can contain weights that modify the probability (i.e. it is more likely to transition between certain phonemes than others)
• priors - initial probailities of starting a sequence in a certain state (this is often skipped whithout much problems)

To keep things readable, we will use the same notation as in Rabiner's paper:

• for observations we will use the letter b
• for transitions the letter a
• for priors the word pi

## Example model

Here we will generate some random values to demonstrate the algorithms. Our HMM wil have 3 states and we will feed 5 observations into it. Each observation is a descrete event of on of 4 possible values. The model will also connect all the states to all the others. We are allowed to start and finish at any state.

#### Constants

These are the constants from above using notation from the Rabiner paper:

In [2]:

N=3
M=4
T=5

#### Observation sequence

We will generate a random sequence of observations:

In [3]:

O=np.random.randint(0,M,T)
print O

[0 1 2 1 0]

The observation matrix will convert these events into probablities for individual states. The matrix is a 3x4 matrix, where the first dimension is the state and second is the observation id:

In [4]:

b=np.random.random((N,M))
b/=np.array([b.sum(axis=0)])
print b

[[ 0.53375578  0.10345127  0.33152859  0.66017725]
[ 0.26331038  0.65987904  0.40656356  0.11079022]
[ 0.20293384  0.23666969  0.26190785  0.22903253]]

Note that probabilites for individual events have to sum to 1, i.e. each event is modeled by a probability distribution over the states (not the other way around):

In [5]:

print b.sum(axis=0)

[ 1.  1.  1.  1.]

Given thes probabilities, we can easily estimate the observation probabilities in time:

In [6]:

ob=np.zeros((T,N))
for t in range(T):
for i in range(N):
ob[t,i]=b[i,O[t]]
print ob

[[ 0.53375578  0.26331038  0.20293384]
[ 0.10345127  0.65987904  0.23666969]
[ 0.33152859  0.40656356  0.26190785]
[ 0.10345127  0.65987904  0.23666969]
[ 0.53375578  0.26331038  0.20293384]]

Or in simpler, Numpy terms:

In [7]:

ob=b.T[O]
print ob

[[ 0.53375578  0.26331038  0.20293384]
[ 0.10345127  0.65987904  0.23666969]
[ 0.33152859  0.40656356  0.26190785]
[ 0.10345127  0.65987904  0.23666969]
[ 0.53375578  0.26331038  0.20293384]]

These also sum to 1 at each time step:

In [8]:

print ob.sum(axis=1)

[ 1.  1.  1.  1.  1.]

Please note that in our case, we will reaplace this lookup table with a proper model, e.g. an MLP or something better. This simple model will let us demonstrate the Baum-Welch algorithm, which is normally not used when talking about hybrid ANN-HMM models.

#### State transition matrix

The state transition matrix is 3x3 and its rows also have to sum to 1 (i.e. the sum of probabilities of going from one state to all the others has to be equal one):

In [9]:

a=np.random.random((N,N))
a/=np.array([a.sum(axis=-1)]).T
print a
print a.sum(axis=1)

[[ 0.44056368  0.08175695  0.47767937]
[ 0.4364632   0.40292622  0.16061058]
[ 0.00281623  0.88340507  0.1137787 ]]
[ 1.  1.  1.]

#### Priors

We will also generate the priors. You can use random values or set everything to equally likely. The code below should work either way:

In [10]:

pi=np.random.random(N)
pi/=pi.sum()
print pi

[ 0.13060479  0.574861    0.29453421]

#### Final model

To make things easier to read below, we will merge all the values above into a single tuple:

In [11]:

hmm=(O,pi,a,b,N,M,T)

Here is some short code to save the model to a file:

In [12]:

import pickle

with open('../data/hmm.pkl','w') as f:
pickle.dump(hmm,f,pickle.HIGHEST_PROTOCOL)

### Forward algorithm

Implementing the algorithm directly using Rabiner's paper would look like this. Note that Rabiner iterates his vectors/matrices starting from 1 (which is a standard in mathematics), but Python (as most computer programming languages) iterates from 0, so there is a little offset everywhere:

In [13]:

def forward(O,pi,a,b,N,M,T):
fwd=np.zeros((T,N))

#initiazliation
for i in range(N):
fwd[0,i]=pi[i]*b[i,O[0]]

#induction
for t in range(T-1):
for j in range(N):
s=0
for i in range(N):
s+=fwd[t,i]*a[i,j]
fwd[t+1,j]=s*b[j,O[t+1]]

return fwd

print forward(*hmm)

[[ 0.06971106  0.15136687  0.05977096]
[ 0.01002924  0.07884961  0.01524421]
[ 0.01288864  0.01872524  0.00502583]
[ 0.00143438  0.00860381  0.0023042 ]
[ 0.00234515  0.00147968  0.00047268]]

Numpy allows us to implement this a bit more efficiently:

In [14]:

def forward(O,pi,a,b,N,M,T):
fwd=np.zeros((T,N))

#initialization
fwd[0]=pi*b[:,O[0]]

#induction:
for t in range(T-1):
fwd[t+1]=np.dot(fwd[t],a)*b[:,O[t+1]]

return fwd

print forward(*hmm)

[[ 0.06971106  0.15136687  0.05977096]
[ 0.01002924  0.07884961  0.01524421]
[ 0.01288864  0.01872524  0.00502583]
[ 0.00143438  0.00860381  0.0023042 ]
[ 0.00234515  0.00147968  0.00047268]]

In [15]:

def forward(O,pi,a,b,N,M,T):
fwd=np.zeros((T,N))

#initialization
fwd[0]=pi*b[:,O[0]]

#induction:
for t in range(T-1):
fwd[t+1]=np.dot(fwd[t],a)*b[:,O[t+1]]

return fwd

print forward(*hmm)

[[ 0.06971106  0.15136687  0.05977096]
[ 0.01002924  0.07884961  0.01524421]
[ 0.01288864  0.01872524  0.00502583]
[ 0.00143438  0.00860381  0.0023042 ]
[ 0.00234515  0.00147968  0.00047268]]

This allows us to determine the joint probability of the sequence given the model (simply a sum of all the forward probabilities at the last time step):

In [16]:

def full_prob(fwd):
return fwd[-1].sum()

print full_prob(forward(*hmm))

0.00429750109789

### Viterbi algorithm

Note that usually when we classifiy using MLPs directly (without HMMs), we just pick the most probable class at each timestep. This gives us the following sequnce:

In [17]:

np.argmax(ob,axis=1)

Out[17]:

array([0, 1, 1, 1, 0])

The sequence above doesn't take any context into account, however. The HMM model tells us that it may be more likely to choose certrain transitions than others. The Viterbi algorithm lets us find the most likely sequence of states (classes), while taking the parameters of the HMM into account:

In [18]:

def viterbi(O,pi,a,b,N,M,T):
d=np.zeros((T,N))
ph=np.zeros((T,N),dtype=np.int)

#initialization
for i in range(N):
d[0,i]=pi[i]*b[i,O[0]]
ph[0,i]=0

#recursion
for t in range(1,T):
for j in range(N):
m=np.zeros(N)
for i in range(N):
m[i]=d[t-1,i]*a[i,j]
ph[t,j]=m.argmax()
d[t,j]=m.max()*b[j,O[t]]

#termination
m=np.zeros(N)
for i in range(N):
m[i]=d[T-1,i]
Pv=m.max()

#path back-tracking
Q=np.zeros(T,dtype=np.int)
Q[T-1]=m.argmax()
for t in reversed(range(T-1)):
Q[t]=ph[t+1,Q[t+1]]

return Q

print viterbi(*hmm)

[1 1 1 1 0]

And a slightly faster Numpy implementation:

In [19]:

def viterbi(O,pi,a,b,N,M,T):
d=np.zeros((T,N))
ph=np.zeros((T,N),dtype=np.int)

#initialization
d[0]=pi*b[:,O[0]]
ph[0]=0

#recursion
for t in range(1,T):
m=d[t-1]*a.T
ph[t]=m.argmax(axis=1)
d[t]=m[np.arange(N),ph[t]]*b[:,O[t]]

#termination
Q=np.zeros(T,dtype=np.int)
Q[T-1]=np.argmax(d[T-1])
Pv=d[T-1,Q[T-1]]

#path back-tracking
for t in reversed(range(T-1)):
Q[t]=ph[t+1,Q[t+1]]

return Q

print viterbi(*hmm)

[1 1 1 1 0]

It should also be clear that if we set both the priors and transitions to be equiprobable (all the same value) the best observation and Viterbi paths will be the same.

### Forward-backward routine

The forward algorithm is already implemented above, so here, we will show the backward here:

In [20]:

def backward(O,pi,a,b,N,M,T):
bk=np.zeros((T,N))
#initialization
for i in range(N):
bk[T-1,i]=1
#induction:
for t in reversed(range(T-1)):
for i in range(N):
s=0
for j in range(N):
s+=a[i,j]*b[j,O[t+1]]*bk[t+1,j]
bk[t,i]=s
return bk

print backward(*hmm)

[[ 0.00983429  0.01427416  0.02428114]
[ 0.04165409  0.03925479  0.05146331]
[ 0.06524467  0.12455948  0.22368039]
[ 0.35361817  0.37165273  0.25720245]
[ 1.          1.          1.        ]]

And in Numpy:

In [21]:

def backward(O,pi,a,b,N,M,T):
bk=np.zeros((T,N))
#initialization:
bk[T-1]=1
#induction:
for t in reversed(range(T-1)):
bk[t]=np.dot(bk[t+1]*b[:,O[t+1]],a.T)
return bk

print backward(*hmm)

[[ 0.00983429  0.01427416  0.02428114]
[ 0.04165409  0.03925479  0.05146331]
[ 0.06524467  0.12455948  0.22368039]
[ 0.35361817  0.37165273  0.25720245]
[ 1.          1.          1.        ]]

When we have both Forward and Backward probabilities computed, we can use them together to estimate the probability of each state, at each point in time, by simply multuplying those two probablilities together. This works, because is the probability of a sequence from time 0 to time $t$ and backward is the probabilty from time $t$ to the end $T$. Adding them together gives a complete probability of passing that state in time. Addtionally, we normalize this by the joint probability $P$ to make sure everything sums up nicely (i.e. each timestep sums to 1):

In [22]:

def gamma(fwd,bk,fp):
gm=np.zeros((T,N))
for t in range(T):
for i in range(N):
gm[t,i]=(fwd[t,i]*bk[t,i])/fp
return gm

print gamma(forward(*hmm),backward(*hmm),full_prob(forward(*hmm)))

[[ 0.15952497  0.50276549  0.33770954]
[ 0.09720975  0.72023827  0.18255198]
[ 0.19567537  0.54273542  0.26158921]
[ 0.11802761  0.74406748  0.13790491]
[ 0.5457006   0.34431105  0.10998835]]

Numpy one-liner:

In [23]:

def gamma(fwd,bk,fp):
return (fwd*bk)/fp

print gamma(forward(*hmm),backward(*hmm),full_prob(forward(*hmm)))

[[ 0.15952497  0.50276549  0.33770954]
[ 0.09720975  0.72023827  0.18255198]
[ 0.19567537  0.54273542  0.26158921]
[ 0.11802761  0.74406748  0.13790491]
[ 0.5457006   0.34431105  0.10998835]]

We can compare the three probailities using a plot. You can notice on the plot (and it should be obvious from the equations) that forward probabilities decay exponentially in time and similarly (but in the opposite direction) with the backward probabilities. Adding them together, evens everything out and the full path denotes clearly what the most likely state at each step is.

In [24]:

fig,ax=P.subplots(1,3,figsize=(15,3),sharey=True)

ax[0].pcolormesh(forward(*hmm).T,cmap=P.cm.gray)
ax[0].set_yticks(np.arange(0,N)+.5)
ax[0].set_yticklabels(np.arange(0,N))
ax[0].set_title('Forward')

ax[1].pcolormesh(backward(*hmm).T,cmap=P.cm.gray)
ax[1].set_title('Backward')

g=gamma(forward(*hmm),backward(*hmm),full_prob(forward(*hmm)))
ax[2].pcolormesh(g.T,cmap=P.cm.gray)
ax[2].set_title('Full')

Out[24]:

<matplotlib.text.Text at 0x7fb43f5a47d0>

Please note there are many sequences we can extract now and they are all different and NOT equivalent in any way, so you cannot use them interchangebly.

• we can get the best observation at each timestep, which doesn't take the transitions or priors of the HMM into account
• we can take the best forward or backward probability at each time step, which takes only the history/future into account
• we can take the best full-probability at each time step, but this isn't actually the most likely sequence - it is just the most likely state, i.e. there are many possible sequences and that is the average probability of all the sequences (likely or not) that go through that state
• finally, we can get simply the most likely sequence using the Viterbi algorithm

The problem with full probability argmax not being equivalent to the Viterbi sequence can also be explained by the fact that the two states in the full prob. sequence don't have to even be a part of the same sequence.

In [25]:

print 'Best observation sequence: {}'.format(np.argmax(ob,axis=1))
print 'Best forward prob. sequence: {}'.format(np.argmax(forward(*hmm),axis=1))
print 'Best backward prob. sequence: {}'.format(np.argmax(backward(*hmm),axis=1))
g=gamma(forward(*hmm),backward(*hmm),full_prob(forward(*hmm)))
print 'Best full prob. sequence: {}'.format(np.argmax(g,axis=1))
print 'Best Viterbi sequence: {}'.format(viterbi(*hmm))

Best observation sequence: [0 1 1 1 0]
Best forward prob. sequence: [1 1 1 1 0]
Best backward prob. sequence: [2 2 2 1 0]
Best full prob. sequence: [1 1 1 1 0]
Best Viterbi sequence: [1 1 1 1 0]

### Baum-Welch

Now we can also determine the sequence probablity, which would help us with fine-tuning the state transition likelihood. You can read about the Baum-Welch algorithm in Rabiner's paper, but this is something that is not really used much in hybrid ANN-HMM models. We are mentioning it here more for sake of curiousity:

The transition probability computed for each timestep gives us a 3D matrix of size 4x3x3, since there are $N^2$ transitions at each timestep and only $T-1$ transitions time steps (i.e. number of transitions between 5 consecutive timesteps). The matrix is computed thusly:

In [26]:

def xi(fwd,bk,fp,O,pi,a,b,N,M,T):
ret=np.zeros((T-1,N,N))
for t in range(T-1):
for i in range(N):
for j in range(N):
ret[t,i,j]=(fwd[t,i]*a[i,j]*b[j,O[t+1]]*bk[t+1,j])/fp
return ret

print xi(forward(*hmm),backward(*hmm),full_prob(forward(*hmm)),*hmm)

[[[  3.07955404e-02   3.43532117e-02   9.43762150e-02]
[  6.62454248e-02   3.67618460e-01   6.89016034e-02]
[  1.68785883e-04   3.18266595e-01   1.92741635e-02]]

[[  2.22395996e-02   9.66233633e-03   6.53078151e-02]
[  1.73219688e-01   3.74381469e-01   1.72637109e-01]
[  2.16084372e-04   1.58691610e-01   2.36442879e-02]]

[[  4.83359429e-02   6.01336344e-02   8.72057951e-02]
[  6.95711863e-02   4.30564850e-01   4.25993795e-02]
[  1.20484393e-04   2.53368996e-01   8.09973173e-03]]

[[  7.84874698e-02   7.18524727e-03   3.23548965e-02]
[  4.66407165e-01   2.12406845e-01   6.52534701e-02]
[  8.05964047e-04   1.24718959e-01   1.23799829e-02]]]

Numpy version:

In [27]:

def xi(fwd,bk,fp,O,pi,a,b,N,M,T):
return fwd[:-1].reshape((T-1,N,1))*a.reshape((1,N,N))*b[:,O[1:]].T.reshape((T-1,1,N))*bk[1:].reshape((T-1,1,N))/fp

print xi(forward(*hmm),backward(*hmm),full_prob(forward(*hmm)),*hmm)

[[[  3.07955404e-02   3.43532117e-02   9.43762150e-02]
[  6.62454248e-02   3.67618460e-01   6.89016034e-02]
[  1.68785883e-04   3.18266595e-01   1.92741635e-02]]

[[  2.22395996e-02   9.66233633e-03   6.53078151e-02]
[  1.73219688e-01   3.74381469e-01   1.72637109e-01]
[  2.16084372e-04   1.58691610e-01   2.36442879e-02]]

[[  4.83359429e-02   6.01336344e-02   8.72057951e-02]
[  6.95711863e-02   4.30564850e-01   4.25993795e-02]
[  1.20484393e-04   2.53368996e-01   8.09973173e-03]]

[[  7.84874698e-02   7.18524727e-03   3.23548965e-02]
[  4.66407165e-01   2.12406845e-01   6.52534701e-02]
[  8.05964047e-04   1.24718959e-01   1.23799829e-02]]]

The expected prior value can be estimated by simply looking at the expected value $\gamma$ at the first point in time:

In [28]:

def exp_pi(gamma):
return gamma[0]

print exp_pi(gamma(forward(*hmm),backward(*hmm),full_prob(forward(*hmm))))

[ 0.15952497  0.50276549  0.33770954]

The expected transistions is the sum (for different timestep) of transition probabilities $\xi$ normalized by the corresponding emitting state probabilities $\gamma$:

In [29]:

def exp_a(gamma,xi,N):
e_a=np.zeros((N,N))
for i in range(N):
for j in range(N):
e_a[i,j]=np.sum(xi[:,i,j])/np.sum(gamma[:-1,i])

return e_a

fw=forward(*hmm)
bk=backward(*hmm)
fp=full_prob(fw)
g=gamma(fw,bk,fp)
x=xi(fw,bk,fp,*hmm)

ea=exp_a(g,x,N)

print exp_a(g,x,N)

[[ 0.3152992   0.19517369  0.48952711]
[ 0.30896542  0.55182403  0.13921055]
[ 0.00142573  0.92964492  0.06892936]]

Numpy version:

In [30]:

def exp_a(gamma,xi,N):
return xi[:].sum(axis=0)/gamma[:-1].sum(axis=0).reshape(N,1)

fw=forward(*hmm)
bk=backward(*hmm)
fp=full_prob(fw)
g=gamma(fw,bk,fp)
x=xi(fw,bk,fp,*hmm)

print exp_a(g,x,N)

[[ 0.3152992   0.19517369  0.48952711]
[ 0.30896542  0.55182403  0.13921055]
[ 0.00142573  0.92964492  0.06892936]]

Finally, the expected observation likelihood is sum of probabilities of given states at times where the specific observation occured normalized by the sum of all the probilities in time for the given state:

In [31]:

def exp_b(gamma,O,N,M):
e_b=np.zeros((N,M))
for j in range(N):
for k in range(M):
e_b[j,k]=np.sum(gamma[O==k,j])/np.sum(gamma[:,j])
return e_b

fw=forward(*hmm)
bk=backward(*hmm)
fp=full_prob(fw)
g=gamma(fw,bk,fp)

print exp_b(g,O,N,M)

[[ 0.63184425  0.19284112  0.17531463  0.        ]
[ 0.29679103  0.51305023  0.19015874  0.        ]
[ 0.43476621  0.31120054  0.25403325  0.        ]]

Numpy version (couldn't figure out how to get rid of the k loop):

In [32]:

def exp_b(gamma,O,N,M):
return np.array(map(lambda k: np.sum(gamma[O==k],axis=0)/np.sum(gamma,axis=0), np.arange(M))).T

fw=forward(*hmm)
bk=backward(*hmm)
fp=full_prob(fw)
g=gamma(fw,bk,fp)

print exp_b(g,O,N,M)

[[ 0.63184425  0.19284112  0.17531463  0.        ]
[ 0.29679103  0.51305023  0.19015874  0.        ]
[ 0.43476621  0.31120054  0.25403325  0.        ]]

The Baum-Welch algorithm now works by repeadatly updating the model parameters in this two step EM (first Expectation, then Modification) procedure. The goal of the algorithm is to maximize the $P(O|\Lambda)$, where $O$ is the observation sequence and $\Lambda$ is the model. We will also print out the mean square of the difference of model parameters vs their expected value:

In [49]:

print 'Initial probability: {}'.format(full_prob(forward(*hmm)))

hmm_new=hmm
for i in range(15):
fw=forward(*hmm_new)
bk=backward(*hmm_new)
fp=full_prob(fw)
g=gamma(fw,bk,fp)
x=xi(fw,bk,fp,*hmm_new)

pi_new=exp_pi(g)
a_new=exp_a(g,x,N)
b_new=exp_b(g,O,N,M)

err=np.concatenate(((pi_new-hmm_new[1]).ravel(),(a_new-hmm_new[2]).ravel(),(b_new-hmm_new[3]).ravel()))

hmm_new=(O,pi_new,a_new,b_new,N,M,T)

print 'Update #{} probability: {} -- mean error: {}'.format(i+1,full_prob(forward(*hmm_new)),np.mean(err**2))

Initial probability: 0.00429750109789
Update #1 probability: 0.00509988830851 -- mean error: 0.0312950732496
Update #2 probability: 0.00619259912726 -- mean error: 0.00250484514145
Update #3 probability: 0.00862175382808 -- mean error: 0.00354749899731
Update #4 probability: 0.0136389139542 -- mean error: 0.00476954018424
Update #5 probability: 0.0189215113683 -- mean error: 0.00249158166715
Update #6 probability: 0.0234132960784 -- mean error: 0.00156689403255
Update #7 probability: 0.037457169758 -- mean error: 0.00485260629198
Update #8 probability: 0.109155856481 -- mean error: 0.00896292597035
Update #9 probability: 0.21859410717 -- mean error: 0.00389884346838
Update #10 probability: 0.238835338999 -- mean error: 0.000297478434192
Update #11 probability: 0.244890426118 -- mean error: 7.89287339611e-05
Update #12 probability: 0.247751721958 -- mean error: 1.83136278322e-05
Update #13 probability: 0.249030029832 -- mean error: 3.72290018409e-06
Update #14 probability: 0.24958516681 -- mean error: 7.08709470271e-07
Update #15 probability: 0.249823249994 -- mean error: 1.30928780085e-07

Now, this doesn't do anything for optimizing the Viterbi sequence recognized by the model, but there is another way HMMs are trained for speech recognition. Usually, each phonetical unit is modelled by a simple 3-5 state left-to-right model. When we present the recognizer a sentence, we actually combine all these phonetic models into one long chain. Then we try to maximize the $P(O|\Lambda_{\rho})$ for each phoneme (indexed by $\rho$) individually, using the above E-M reestimation procedure.

With hybird ANN-HMM systems, we generally don't use EM training. ANNs are usually trained discriminatevly. In that case, we need a frame-wise alignment to train the model. We can either iteratively re-align the data (using Viterbi) and train the system, or we can use some other cost function (e.g. CTC), skipping the HMM part entirely.

## HMM topologies for speech recognition

To actually perform speech recognition using HMMs, we need to expand the definition of the model considerably. First of all, this model will allow us only to recognize a simple sequence of events.

Say we create an "acoustic model" (instead of the observation martix) that will recognize phonemes based on some acoustic observations (e.g. STFT filterbank outputs through time). We can easily do this using a simple MLP. Such a model can use HMMs to cconvert the raw MLP outputs into a sequence of phonemes. If we wish to recognize words, we need to create another model to convert that sequence of phonemes into their word representation.

Since a sequence of words is the same kind of sequence as the sequence of phonemes (albeit a bit more sparse), we can model it simply as another layer of HMMs on top of the phoneme one. This layer will not allow any sequence of phonemes, but only the ones that form actual words from a lexicon. We can also add one more layer that will define how words can interact with eachother in a sentence - the so-called language model.

All these layers can get very tricky. In practice, we may decide to acutally recreate this layered model in a single-layered complex graph containg all the weights and contraints of both the phonetic, lexical and language layers together. This is named, by some reaserchers, an "Integrated Network" or IN for short. This turns out to be not enough in order to fully realize our solution, as the hidden states are still phonemes and not words. We can therefore add special non-emitting states (states that don't move the time-line) to denote word bnoundaries - times where a word has to be output.

A slightly more efficient solution would be to instead add an additional field to each state, apart from the phoneme that we are trying to consume, a word that we are trying to output. For most states, this auxilary value would be blank, but every once in a while we would have a word that needs to be output.

Thus we have successfully redefined our classical HMM solution that relied on Finite-State Automata (FSA) to a much more convinient solution that relies of Finite-State Transducers (FST).

# WFST

Weighted Finite State Transducers are a type of automaton that has a sequence of states with transitions between them, just like an HMM. Unlike an HMM, a transducer has two types of symbols: an input and an output.

HMM problems are usually presented as a graph with hidden states, denoting our hidden value that we are trying to deduce, and a weight related to the transition between these states. A different and pretty much equivalent method to represent this same model is to rewrite the hidden symbols on the transitions, together with the weights. This is also how a transducer is normally depicted, with a sequence of transitions between states, where each transition has two symbols (input and output) and a weight.

All the details about WFST and their application to speech recognition can be found in this paper by Mohri, et al. One of the prime examples of WFST for speech recognition is the Kaldi project.

### Transducer example

An actual transducer is an automaton that converts symbols from one set of values to another. Take this transducer that converts symbols in a following manner:

• convert a to b
• convert b to c
• convert c to a

By using the excellent PyFst library, we can visualize this transducer like this:

In [34]:

f=fst.Transducer()
f[0].final=True
f

Out[34]:

FST

0

0

0->0

a:b

0->0

b:c

0->0

c:a

Now let's make another transducer, but this will have the same input and output symbols and they will form a sequence.

Note that a transducer with identical input and output symbols is simply an automaton (FSA), or otherwise known as an acceptor.

We will use a helpful linear_chain method to generate such an acceptor:

In [35]:

i=fst.linear_chain(['a','a','c','b'],syms=f.isyms)
i

Out[35]:

FST

0

0

1

1

0->1

a

2

2

1->2

a

3

3

2->3

c

4

4

3->4

b

### Composition

The most useful method used in a transducer is composition. This operation takes two transducers, e.g. one converting from symbol space A into B, and second converting from B into C and creates a third transducer that converts from A into C that is like a combination of the original ones.

We can use this method (denoted by the >> operator) to process the input chain transducer using our token converter above. The project_output method simply copies the output symbols to the input in order to discard the input symbols:

In [36]:

o=i>>f
o.project_output()
o

Out[36]:

FST

0

0

1

1

0->1

b

2

2

1->2

b

3

3

2->3

a

4

4

3->4

c

It pays off to think about this for a while until it's clear what is actually going on here. For more help, I recommend reading some of the docs at the OpenFST website.

We can also process several sequences at once. Here we create a transducer that is a union of 3 sub-transducers. We will also expand our definition of the new symbols. We have to remove epsilon, because (apart from being uneccessary) jupyter notebook has issues with printing unicode characters in SVG:

In [37]:

f=fst.Transducer()
f[0].final=True

i1=fst.linear_chain(['a','b','c','d','e'],syms=f.isyms)
i2=fst.linear_chain(['a','b','x','y','z'],syms=f.isyms)
i3=fst.linear_chain(['g','h','c','d','e'],syms=f.isyms)
i=i1.union(i2)
i=i.union(i3)
i.remove_epsilon()
i

Out[37]:

FST

0

0

6

6

0->6

a

11

11

0->11

g

1

1

0->1

a

7

7

6->7

b

12

12

11->12

h

2

2

1->2

b

3

3

2->3

c

4

4

3->4

d

5

5

4->5

e

8

8

7->8

x

9

9

8->9

y

10

10

9->10

z

13

13

12->13

c

14

14

13->14

d

15

15

14->15

e

Now, this graph can be made more efficient by using determinization and minimization. Note that the graph below is equivalent to the one above, but has less states and transitions. This is especially useful with big graphs:

In [38]:

i=i.determinize()
i.minimize()
i

Out[38]:

FST

0

0

1

1

0->1

a

2

2

0->2

g

3

3

1->3

b

4

4

2->4

h

6

6

3->6

c

5

5

3->5

x

4->6

c

8

8

6->8

d

7

7

5->7

y

9

9

7->9

z

8->9

e

And we can still produce an output, like before:

In [39]:

o=i>>f
o.project_output()
o

Out[39]:

FST

0

0

1

1

0->1

b

2

2

0->2

e

3

3

1->3

c

4

4

2->4

f

5

5

3->5

a

6

6

3->6

x

4->5

a

7

7

5->7

f

8

8

6->8

y

9

9

7->9

g

8->9

z

This method is also useful to combine a chain of operations into one complete transducer, for example to create a speech recognition pipeline as mentioned in the chapter about HMMs above. We can create a separate transducer for modeling phonemes, one for converting phonemes into words and another one for modeling sequences of words and combine all three into one large transducer converting context-indepentend phonemes into context-depndent words.

## HMM example using WFST

Let's reitrate our HMM example above using WFSTs. First lets convert the observation matrix into an input sequence:

In [40]:

def obs_mat(O,b,N,T):
f=fst.StdAcceptor()

for t in range(T):
for j in range(N):

f[T].final=True

return f

obs_fsa=obs_mat(O,b,N,T)

A thing to note is that weights in FSTs are usually represented in a different way than in HMMs. We generally define weights in a specific semiring: in this case the tropical semiring.

This semiring defines that the combination of two weigts is implemented using the sum. Probabilities are usually added togehter by multiplying them with each other. If we convert them into the log domain, the multiplication changes into a sum. Taking a logarithm of a value between 0 and 1 (i.e. probability) caonstraints the result to values between $-\inf$ to 0. The value of 0 is the best, and $-\inf$ is the worst. The tropical semiring works in the opposite direction, so we also have to take the negative of the logarithm to redefine values from 0 to $\inf$ and make the smaller values better than higher:

In [41]:

obs_fsa

Out[41]:

FST

0

0

1

1

0->1

0/0.627817

0->1

1/1.33442

0->1

2/1.59488

2

2

1->2

0/2.26865

1->2

1/0.415699

1->2

2/1.44109

3

3

2->3

0/1.10404

2->3

1/0.900015

2->3

2/1.33976

4

4

3->4

0/2.26865

3->4

1/0.415699

3->4

2/1.44109

5

5

4->5

0/0.627817

4->5

1/1.33442

4->5

2/1.59488

This FSA simply chains several sequences in parallel. Each state represent a time step and the transitions in between are the individual observations. It's shape is what made it known as the sausage.

Now, we can also incorporate the priors and transitions in another FST that looks like this:

In [42]:

def trans_mat(pi,a,N):
f=fst.StdAcceptor()

for i in range(N):

for i in range(N):
for j in range(N):

for i in range(N):
f[i+1].final=True

return f

tr_fsa=trans_mat(pi,a,N)

In [43]:

tr_fsa

Out[43]:

FST

0

0

1

1

0->1

0/2.03558

2

2

0->2

1/0.553627

3

3

0->3

2/1.22236

1->1

0/0.8197

1->2

1/2.504

1->3

2/0.738816

2->1

0/0.829051

2->2

1/0.909002

2->3

2/1.82877

3->1

0/5.87236

3->2

1/0.123971

3->3

2/2.1735

Now lets combine the input sausage with the above transitions:

In [44]:

out_fsa=obs_fsa>>tr_fsa
out_fsa

Out[44]:

FST

0

0

1

1

0->1

0/2.6634

2

2

0->2

1/1.88805

3

3

0->3

2/2.81724

4

4

1->4

0/3.08835

5

5

1->5

1/2.9197

6

6

1->6

2/2.17991

2->4

0/3.09771

2->5

1/1.3247

2->6

2/3.26986

3->4

0/8.14101

3->5

1/0.53967

3->6

2/3.61459

7

7

4->7

0/1.92374

8

8

4->8

1/3.40402

9

9

4->9

2/2.07858

5->7

0/1.93309

5->8

1/1.80902

5->9

2/3.16854

6->7

0/6.9764

6->8

1/1.02399

6->9

2/3.51326

10

10

7->10

0/3.08835

11

11

7->11

1/2.9197

12

12

7->12

2/2.17991

8->10

0/3.09771

8->11

1/1.3247

8->12

2/3.26986

9->10

0/8.14101

9->11

1/0.53967

9->12

2/3.61459

13

13

10->13

0/1.44752

14

14

10->14

1/3.83843

15

15

10->15

2/2.33369

11->13

0/1.45687

11->14

1/2.24342

11->15

2/3.42365

12->13

0/6.50017

12->14

1/1.45839

12->15

2/3.76838

What we get is a graph representing all the possible transitions among between all possible states with their corresponging weights. If we use the shortest_path method on this graph, we will get the same path as by using Viterbi:

In [45]:

sh_p=out_fsa.shortest_path()
sh_p

Out[45]:

FST

5

5

4

4

5->4

1/1.88805

3

3

4->3

1/1.3247

0

0

1

1

1->0

0/1.45687

2

2

2->1

1/1.3247

3->2

1/1.80902

In [46]:

shp_l=[]
for p in sh_p.paths():
for s in p:
shp_l.append(sh_p.osyms.find(s.olabel))
print shp_l

[u'1', u'1', u'1', u'1', u'0']

In [47]:

print viterbi(*hmm)

[1 1 1 1 0]

In [ ]:


text:
# Deriving the Time Complexity of Ryser's Algorithm for Evaluating the Permanent of a Matrix

Ryser has shown that the permanent of an $$n \times n$$ matrix $$A=(a_{ij})$$ can be expressed as

\begin{align} Perm(A) = (-1)^n \sum_{s \subset [n]} (-1)^{|s|} \prod_{i=1}^n \sum_{j \in s} a_{ij}, \end{align}

where $$[n]=\{1,2,\dots,n\}$$. This algorithm runs in $$\mathcal{O}(2^{n-1}n^2)$$ time. I've been trying to derive this, but can't quite get the result. Here's my work so far.

The outer sum is over all non-empty subsets of $$[n]$$, of which there are $$2^n-1$$. We recall that the number of subsets of size $$r$$ is $${n \choose r}=\frac{n!}{r!(n-r)!}$$.

For each set $$s$$ in the outer sum we compute $$\sum_{j \in s} a_{ij}$$, which uses $$|s|-1$$ additions. Next, this sum sees the product $$\prod_{i=1}^n$$, which takes $$n-1$$ multiplications. This happens for each non-empty subset of $$[n]$$, so there are the following number of total additions:

\begin{align} \sum_{s \subset [n]} \left(|s|-1\right) &= \sum_{k=1}^n (k-1) {n \choose k} \\ &= \sum_{k=1}^n k {n \choose k} - \sum_{k=1}^n {n \choose k} \\ &= n \sum_{k=1}^n {n-1 \choose k-1} - \left(2^n - 1\right) \\ &= n \left( 2^{n-1}-1 \right) - \left(2^n - 1\right) \\ &=n2^{n-1} - 2^n -n +1. \end{align}

The total number of multiplications is $$(n-1)(2^n-1)$$.

There is also the $$(-1)^{|s|}$$, which takes another $$2^n-1$$ operations in total.

This gives us a grand total of $$n2^n + n 2^{n-1} - 2^n -2n+1$$, which is not correct. It looks like I'm off by a factor of $$n$$ on the second term here. Where am I going wrong?

For a given $$s$$, you do $$(|s|-1)n$$ additions, not $$|s|-1$$ additions: you have to compute the inner sum for each $$i\in [n]$$. Thus, you are missing a factor of $$n$$.


text:
# Notes for JEE Main & Advanced Mathematics

## Recent Notes

Matrices of Rotation of Axes
Geometrical Transformations
Cayley-Hamilton Theorem
Consistency of a System of Linear Equation $\mathbf{AX=B,}$ where $\mathbf{A}$ is a square matrix
Homogeneous and Non-homogeneous Systems of Linear Equations
Echelon Form of a Matrix
Rank of Matrix
Inverse of a Matrix
Special Types of Matrices
Transpose of a Matrix
Positive Integral Powers of a Matrix
Multiplication of Matrices
Scalar Multiplication of Matrices
Trace of a Matrix
Types of Matrices
Equality of Matrices
Order of a Matrix
Definition
Some Special Determinants
Application of Determinants in Solving a System of Linear Equations
Differentiation and Integration of Determinants
Product of Two Determinants
Minors and Cofactors
Properties of Determinants
Definition
Lagrange's Mean Value Theorem
Rolle's Theorem
Greatest and Least Values of a Function Defined on an Interval $[a,\,\,b]$
Properties of Maxima and Minima
Higher Order Derivative Test
Sufficient Criteria for Extreme Values (1st Derivative Test)
Necessary Condition for Extreme Values
Definition
Definition
Length of Perpendicular from Origin to the Tangent
Length of Intercept Made on Axes by The Tangent
Length of Tangent, Normal, Subtangent, Subnormal
Angle of Intersection of Two Curves
Equation of the Tangent and Normal
Slope of the Tangent and Normal
Velocity and Acceleration in Rectilinear Motion
Derivative as the Rate of Change
Deduction of Euler?s Theorem
Euler's Theorem on Homogeneous Functions
Higher Partial Derivatives
Definition
Leibnitz?s Theorem
Differentiation of Integral Function
${{n}^{th}}$ Derivative Using Partial Fractions
Successive Differentiation or Higher Order Derivatives
Differentiation of a Function with Respect to Another Function
Methods of Differentiation
Theorems for Differentiation
Some Standard Differentiation
Introduction
Relation Between A.P., G.P. and H.P.
Properties of Arithmetic, Geometric, Harmonic Means Between Two Given Numbers
Special Series
Method of Difference
Method for Finding Sum
Sum of A.G.P.
${{n}^{th}}$term of A.G.P.
Definition
Properties of H.P.
Harmonic Mean
General Term of an H.P.
Definition
Properties of G.P.
Geometric Mean
Sum of infinite terms of a G.P.
Sum of first 'n' terms of a G.P.
Selection of Terms in a G.P.
General term of a G.P.
Definition
Properties of A.P.
Arithmetic Mean
Sum of n terms of an A.P.
Selection of Terms in an A.P.
General Term of an A.P.
Definition
Some Important Results
Some Properties Related to Circle
Some Terminology Related to Height and Distance
Solutions of Triangles
Regular Polygon
Ex-central Triangle
Pedal Triangle
Circle Connected With Triangle
Half Angle Formulae
Area of Triangle
Napier's Analogy (Law of Tangents)
Theorem of The Medians: (Apollonius Theorem)
Projection Formulae
The Law of Cosines or Cosine Rule
Relation Between Sides and Angles
Periodic Functions
Important Points to be Taken in Case of While Solving Trigonometrical Equations
Method for Finding Principal Value
General Solution of The Form $\mathbf{acos}\,\mathbf{\theta +bsin}\,\mathbf{\theta =c}$
Trigonometrical Equations With Their General Solution
Trigonometrical Equations and Inequations Definition
Differentiability of a Function at a Point
Discontinuous Function
Continuity From Left and Right
Continuity of a Function at a Point
Introduction
Methods of Evaluation of Limits
Fundamental Theorems on Limits
Limit of a Function
Inverse Function
Composite Function
Periodic Function
Even and Odd Function
Conditional Trigonometrical Identities
Kinds of function
Algebra of Functions
Domain, Co-domain and Range of Function
Definition of Function
Intervals
Maximum and Minimum Value of a $\mathbf{cos}\,\,\mathbf{\theta }\,\,\mathbf{+}\,\mathbf{b}\,\,\mathbf{sin}\,\,\mathbf{\theta }$
Some Important Definitions
Trigonometric Ratio of Sub-multiple of an Angle
Trigonometric Ratio of Multiple of an Angle
Formulae to Rransform The Product Into Sum or Difference
Trigonometrical Ratios In Terms of Each Other
Trigonometrical Ratios for Some Special Angles
Trigonometrical Ratios for Various Angles
Trigonometrical Ratios of Allied Angles
Trigonometrical Ratios or Functions
Domain and Range of a Trigonometrical Function
Angle of Intersection of Two Spheres
Intersection of Straight Line and a Sphere
Condition of Tangency of a Plane to a Sphere
Relation Between an arc and an Angle
Relation Between Three Systems of Measurement of an Angle
System of Measurement of Angles
Section of a Sphere by a Plane
Equation in Sphere in Various Forms
General Equation of Sphere
Projection of a Line on a Plane
Angle Between Line and Plane
Intersection Point of a Line and Plane
Introduction
Equation of Plane Through a Given Line
Coplanar Lines
Image of a Point in a Plane
Equation of Planes Bisecting Angle Between Two Given Planes
Locus
Angle Between Two Planes
Foot of Perpendicular from a Point $\mathbf{A(\alpha ,}\,\,\mathbf{\beta ,}\,\,\mathbf{\gamma )}$to a given plane $\mathbf{ax+by+cz+d=0}$.
Transformation of Axes
Equation of Plane Passing Through the Given Point
Area of Some Geometrical Figures
Definition of Plane and Its Equations
Shortest Distance Between two Straight Lines
Foot of Perpendicular from a Point to The Line
Intersection of Two Lines
Changing Unsymmetrical form to Symmetrical Form
Some Points of a Triangle
Straight Line in Space
Angle Between Two Lines
Projection
Direction Cosines and Direction Ratios
Triangle and Tetrahedron
Section Formulae
Section Formula
Distance Formula
Properties of Some Geometrical Figures
Co-ordinates of a Point in Space
Distance Formula
Polar Co-ordinates
Cartesian Co-ordinates of a Point
Introduction
Binomial Distribution
Total Probability and Baye's Rule
Conditional Probability


text:
# Flickering dashed circle in rotated Manipulate

Bug introduced in 12.0 or earlier and persisting through 12.0.0 or later

I'm drawing some objects (lines, circles etc) in a Manipulate. I am rotating the Manipulate upon showing all the graphics. I have found for a dashed circle any rotation angle other than 0 results in the circle flickering / not being rendered. This does not happen for a solid circle.

I have made a minimal simplified demo below. If you scale the radius of the circles, with a rotation angle of 0 things are fine, dashed and solid line circles are visible. However, for other rotation angles, the dashed, but not solid, circle flickers in and out of visibility when the radius is changed. When static (radius or rotation angle not being changed with sliders), the dashed circle vanishes.

Clear["Global*"]

Manipulate[
(* Make a line and circle *)

thePlaneCircle = Graphics[{Red, Dashed, Thick, Circle[{0, -5}, r]}];
thePlaneCircleSolid = Graphics[{Red, Thick, Circle[{0, 5}, r]}];
theLine = Graphics[Line[{{0, -50}, {0, 50}}]];

(* Rotate and show *)
Rotate[
Show[thePlaneCircle, thePlaneCircleSolid, theLine]
, rot Degree]

, {{r, 5, "Radius"}, 1, 10, Appearance -> "Labeled"}
, {{rot, 0, "Rotation Angle"}, 0, 360, Appearance -> "Labeled"}

]

I'm aware that I could avoid this if I simply did not use Rotate around Show, and exchanged my x and y coordinates, but for the purposes of my more complicated original code, this would be a pain.

However, I am puzzled as to why only the dashed circle flickers in this simplified example, regardless.

I'm reasonably new to Mathematica, so apologies if I have made a glaring silly error.

I'm using Mathematica version 12 on macOS Mojave 10.14.4.

Thanks

Update: I've submitted this a report of a possible bug. In the report pointing to this thread and showing Carls minimal example. Thanks for your help. Will report back if I get a response.

Update: Has been confirmed as a bug by Wolfram. No specific workaround available.

• Looks like a bug in M12 to me. I would report it to support. A simpler example: Rotate[Graphics[{Red, Dashed, Thick, Circle[{0, -5}, 1]}], 74. Degree] Apr 26 '19 at 18:41
• Try adding TrackedSymbols :> {r, rot}. I can't test it myself because I don't have V12 Apr 26 '19 at 21:19
• Thanks @Carl Woll. I'll submit a bug report. Apr 27 '19 at 15:22

I don't know if this will be useful to you, but on V11.3, I get better behavior with this code:

Manipulate[
thePlaneCircle = {Red, Dashed, Thick, Circle[{0, -5}, r]};
thePlaneCircleSolid = {Red, Thick, Circle[{0, 5}, r]};
theLine = Line[{{0, -50}, {0, 50}}];
Rotate[
Graphics[{thePlaneCircle, thePlaneCircleSolid, theLine}],
rot Degree],
{{r, 5, "Radius"}, 1, 10, Appearance -> "Labeled"},
{{rot, 0, "Rotation Angle"}, 0, 360, 1, Appearance -> "Labeled"},
TrackedSymbols :> {r, rot},
Alignment -> Center,
ContentSize -> {475, 475}]
`

• Thanks. Unfortunately this doesn't seem to work with v12. Same if I fix the rotation angle with a static variable (it is simply 90 Degrees; I added it to the manipulate to see what angles were effected). Looks like a bug in Mathematica I guess. I'll see if there are any further comments and submit a bug report if there are not. Apr 27 '19 at 15:19
• Bug report submitted and post updated. I assume I don't "accept" any answers on the thread? Yours is the closest to a solution per se. Not sure of the rules / etiquette. Apr 27 '19 at 15:36


text:
Write the Equation of a Linear Function Part II

Learning Outcomes

• Write the equation of a linear function from a graph
• Write the equation for a linear function from slope given two points

Now that we have written equations for linear functions using both slope-intercept form and point-slope form, we can choose which method to use based on the information we are given. That information may be provided in the form of a graph, a point and a slope, two points, and so on. Look at the graph of the function f below.

We are not given the slope of the line, but we can choose any two points on the line to find the slope. Choose $(0, 7)$ and $(4, 4)$. We can use these points to calculate the slope.

$\begin{array}{l} m & =\dfrac{{y}_{2}-{y}_{1}}{{x}_{2}-{x}_{1}}\hfill \\ & =\dfrac{4 - 7}{4 - 0}\hfill \\ & =-\dfrac{3}{4} \end{array}$

Now we can substitute the slope and the coordinates of one of the points into point-slope form.

$\begin{array}{l}y-{y}_{1}=m\left(x-{x}_{1}\right)\hfill \\ y - 4=-\dfrac{3}{4}\left(x - 4\right)\hfill \end{array}$

If we want to rewrite the equation in slope-intercept form, we would do the following:

$\begin{array}{l}y - 4=-\dfrac{3}{4}\left(x - 4\right)\hfill \\ y - 4=-\dfrac{3}{4}x+3\hfill \\ \text{ }y=-\dfrac{3}{4}x+7\hfill \end{array}$

If we wanted to find the equation of the line in slope-intercept form without first using point-slope form, we could have recognized that the line crosses the y-axis when the output value is $7$. Therefore, b = $7$. We now have the initial value b and the slope so we can substitute m and b into slope-intercept form of a line.

The function is $f\left(x\right)=-\dfrac{3}{4}x+7$, and the linear equation would be $y=-\dfrac{3}{4}x+7$.

How To: Given the graph of a linear function, write an equation to represent the function

1. Identify two points on the line.
2. Use the two points to calculate the slope.
3. Determine where the line crosses the y-axis to identify the y-intercept by visual inspection.
4. Substitute the slope and y-intercept into slope-intercept form of a linear equation.

Example

Write an equation for a linear function f given its graph shown below.

In the following video, we show an example of how to write the equation of a line given its graph.

Example

If f is a linear function, with $f\left(3\right)=-2$ , and $f\left(8\right)=1$, find an equation for the function in slope-intercept form.

In the video we show another example of how to write a linear function given two points written with function notation.

Summary

• Given a graph, you can write the equation of the line it represents.
• This method only works for graphs that have points that are easy to verify by visual inspection.


text:
Question-and-Answer Resource for the Building Energy Modeling Community
Get s tarted with the Help page

# Revision history [back]

### Controlling Setpoint with python pluging

Hello everyone,

I want to get the value of setpoint in an EnergyPlus model with the python plugin. I tried this code but it returns -1:

    self.T1_handle = self.api.exchange.get_variable_handle(state, "Zone Thermostat Cooling Setpoint Temperature", "Perimeter_ZN_1")
T1 = self.api.exchange.get_variable_value(state, self.T1_handle)

Based on this post, the number of inputs is fine I think. When I change the variable_names to others such as "Zone Mean Air Temperature", and "People Occupant Count" It works. But it does not work for "Zone Thermostat Cooling Setpoint Temperature".

Thanks!

### Controlling Setpoint with python pluging

Hello everyone,

I want to get the value of setpoint in an EnergyPlus model with the python plugin. I tried this code but it returns -1:

    self.T1_handle = self.api.exchange.get_variable_handle(state, "Zone Thermostat Cooling Setpoint Temperature", "Perimeter_ZN_1")
T1 = self.api.exchange.get_variable_value(state, self.T1_handle)

Based on this post, the number of inputs is fine I think. When I change the variable_names to others such as "Zone Mean Air Temperature", and "People Occupant Count" It works. But it does not work for "Zone Thermostat Cooling Setpoint Temperature".Temperature". I found the variable_names from .eso file.

Thanks!

### Controlling Setpoint with python pluging

Hello everyone,

I want to get the value of setpoint in an EnergyPlus model with the python plugin. I tried this code but it returns -1:

    self.T1_handle = self.api.exchange.get_variable_handle(state, "Zone Thermostat Cooling Setpoint Temperature", "Perimeter_ZN_1")
T1 = self.api.exchange.get_variable_value(state, self.T1_handle)

Based on this post, the number of inputs is fine I think. When I change the variable_names to others such as "Zone Mean Air Temperature", and "People Occupant Count" It works. But it does not work for "Zone Thermostat Cooling Setpoint Temperature". I found the variable_names from .eso file.

Thanks!

### Controlling getting Setpoint value with python pluging

Hello everyone,

I want to get the value of setpoint in an EnergyPlus model with the python plugin. I tried this code but it returns -1:

    self.T1_handle = self.api.exchange.get_variable_handle(state, "Zone Thermostat Cooling Setpoint Temperature", "Perimeter_ZN_1")
T1 = self.api.exchange.get_variable_value(state, self.T1_handle)

Based on this post, the number of inputs is fine I think. When I change the variable_names to others such as "Zone Mean Air Temperature", and "People Occupant Count" It works. But it does not work for "Zone Thermostat Cooling Setpoint Temperature". I found the variable_names from .eso file.

Thanks!

 5 retagged shorowit 7405 ●63 ●8 https://www.nrel.gov/b...

### getting Setpoint value with python pluging

Hello everyone,

I want to get the value of setpoint in an EnergyPlus model with the python plugin. I tried this code but it returns -1:

    self.T1_handle = self.api.exchange.get_variable_handle(state, "Zone Thermostat Cooling Setpoint Temperature", "Perimeter_ZN_1")
T1 = self.api.exchange.get_variable_value(state, self.T1_handle)

Based on this post, the number of inputs is fine I think. When I change the variable_names to others such as "Zone Mean Air Temperature", and "People Occupant Count" It works. But it does not work for "Zone Thermostat Cooling Setpoint Temperature". I found the variable_names from .eso file.

Thanks!

### getting Setpoint value with python pluging

Hello everyone,

I want to get the value of setpoint in an EnergyPlus model with the python plugin. I tried this code but it returns -1:

    self.T1_handle = self.api.exchange.get_variable_handle(state, "Zone Thermostat Cooling Setpoint Temperature", "Perimeter_ZN_1")
T1 = self.api.exchange.get_variable_value(state, self.T1_handle)

Based on this post, the number of inputs is fine I think. When I change the variable_names to others such as "Zone Mean Air Temperature", and "People Occupant Count" It works. But it does not work for "Zone Thermostat Cooling Setpoint Temperature". I found the variable_names from .eso file.

Thanks!


text:
DiffBind - dba.count - Reads overlapping multiple peaks question
1
0
Entering edit mode
@julienbryois-11773
Last seen 5.7 years ago

Dear Rory,

I am currently using DiffBind to detect differential chromatin accessibility using ATAC-seq data. We are testing a few different ways of doing this analysis and I have a question on the behaviour of dba.count for one of the options.

This option is to create non-overlapping sliding windows of 300bp at the location of the consensus peak set. The issue is that we have relatively long reads (150bp) so many reads are expected to overlap with 2 adjacent windows.

My question is: how is dba.count handling reads overlapping with multiple peaks? Are they counted for each peak? Discarded because they map to multiple peaks? Or maybe the fraction of read length mapping to each peak is counted (e.g: 80% peak 1 and 20% peak 2)?

I unfortunately could not find information on this in the manual or online. I apologise if I missed this information somewhere.

Best regards,

Julien

diffbind • 812 views
0
Entering edit mode
Gord Brown ▴ 650
@gord-brown-5664
Last seen 22 months ago
United Kingdom

Hi,

There are two distinct bits of counting code.  The original code (which will probably be removed when we get to it) counts reads if they overlap the interval by more than half the length of the read (or the interval if it's shorter than the read).

The other code uses "summarizeOverlaps" from the GenomicAlignments package and can be configured by the user, as per the "dba.count" manual page.  (Look for the "bUseSummarizeOverlaps" parameter to dba.count.)

Cheers,

- Gord


text:
# Can we perform matching on post-treatment variables?

I want to match followers of some seed accounts with some random users on Twitter based on observed covariates using the coarsened exact matching method. My goal is to test whether following those seed Twitter accounts would change the outcome or not.

As for observed covariates, I can get, for example, # of tweets and # of followers. But these are post-treatment values, because the Twitter API gives me the current values of these variables. So can I use these post-treatment variables in my matching procedure, or I should only match on pre-treatment variables (e.g. # of tweets before the treatment, # of followers before the treatment, etc.)

I would suggest to match only on pre-treatment variables and not on post-treatment variables. If we match on post-treatment variables it is extremely plausible we induce selection bias in our sample and thus we get bogus insights. Selection bias can appear when one of the variables $C$ we use in matching is related to the treatment $A$ as well as the outcome of interest $Y$.

More formally and quoting Hernan & Robins directly: "selection bias can be defined as the bias resulting from conditioning on the common effect of two variables, one of which is either the treatment or associated with the treatment, and the other is either the outcome or associated with the outcome". Chapt. 8 "Selection Bias" from their upcoming book "Causal Inference". Therefore conditioning on post-treatment variables is a clear case of conditioning on common effects.

Notice that selection bias can be present not only because we have certain eligibility criteria (as those imposed by matching on a post-treatment variable $C$) but also through loss to follow-up; i.e. the drop-out rates (stopping the use of Twitter) might be different between the groups examined. The small paper on "Loss to follow-up" by Dettori offers a very brief, succinct introduction on this.

• So, then what should I do with those post-treatment variables? How should I control for them? Say I count # of tweets before and after treatment. Then matched users according to # of tweets before treatment. Then how should I account for # of tweets after treatment? – msmazh Mar 24 '18 at 19:24
• We do not directly do anything with the post-treatment variables. That being said, they are useful, we can use these post-treatment variables as additional descriptors in the study or maybe use them to propose new hypotheses. We should be very careful not to use the post-treatment variables to immediately test new hypotheses though because that would amount to post-hoc theorising which is another problem in itself. – usεr11852 Mar 25 '18 at 1:26

As others have pointed out, adjusting for a post-treatment covariate can introduce bias in the estimated treatment effect. However, I would like to add that there are corner-case scenarios where adjusting for post-treatment covariate is acceptable. The basic idea is that if the covariate is unaffected of the treatment, the adjusting for it won't introduce bias. By adjusting for it, the precision of the estimated average treatment effect can then be improved. Furthermore, if you have confounding in the form of an unobserved pre-treatment covariate, then adjusting for a post-treatment covariate related to it can reduce the bias of your estimate under certain conditions. The details can be found in Rosenbaum (1984).

That said, I would be very reluctant to adjust for anything that is post-treatment, as the cases where such adjustment is beneficial are likely rare. As a general principle, avoid doing it. And if you plan on doing it, I strongly suggest reading the provided reference first.

References:

Rosenbaum, P. R. (1984). The consquences of adjustment for a concomitant variable that has been affected by the treatment. Journal of the Royal Statistical Society. Series A (General). 147(5): 656-666.

Generally, your intuition is correct. Matching on post-treatment variables, i.e. variables that are determined after the treatment, carries the risk of biasing your estimates due to the bad control problem. This problem results from the fact that, loosely speaking, your control variables are influenced by the treatment and could be regarded as outcome variables themselves. For details on the bad control problem, see, for example, Angrist & Pischke (Mostly Harmless Econometrics, 2008).

• Selection bias is caused by conditioning on a variable $C$ that is a common effect of the treatment $A$ or a cause of the treatment $A$ and of the outcome $Y$ or a cause of the outcome $Y$. In the quoted sentence I refer to post-treatment variables (as the question of the OP does) which are a clear case of potentially common effect. Please check the book Causal Inference by Hernan & Robins, Chapt. 8. "Selection Bias" for more details. – usεr11852 Mar 30 '18 at 23:43
• I added a direct quotation from H&R in my answer that I hope clarifies the point you raised. – usεr11852 Mar 31 '18 at 0:10


text:
What functions, if any, can NOT be approximated by a Taylor series? [duplicate]

I am taking some courses online (machine learning) and the instructor noted that a Taylor series can be used to approximate any function. A student left a comment that a Taylor series cannot approximate ANY function, but noted that the instructor probably stated it unequivocally because it can certainly approximate any function we would need for the class.

I am curious, what functions cannot be approximated by a Taylor series?

Two things lead me to believe that any function CAN be approximated by a Taylor series:

1. In context of regression, it seems I can always continue to add terms and adjust the coefficients on those terms (or solve for them to reduce the error) to continually get a better approximation of points in my data set.
2. The ambiguity of the word approximate suggests that we don't have to meet a strict criteria to say that a Taylor series approximates a function.

With those things in mind, I can also see a point to the opposite view. For example, a function with large discontinuities such as

$f(x) = \begin{cases}sin(x)/x^2 & \text{if$x < 0$} \\ 2x & \text{if$x > 100$}\end{cases}$

seems like it wouldn't be well approximated by a Taylor series. But then again, without some defenition of approximate, I don't see how we can say either way.

marked as duplicate by Xander Henderson, Namaste calculus StackExchange.ready(function() { if (StackExchange.options.isMobile) return; $('.dupe-hammer-message-hover:not(.hover-bound)').each(function() { var$hover = $(this).addClass('hover-bound'),$msg = $hover.siblings('.dupe-hammer-message');$hover.hover( function() { $hover.showInfoMessage('', { messageElement:$msg.clone().show(), transient: false, position: { my: 'bottom left', at: 'top center', offsetTop: -7 }, dismissable: false, relativeToBody: true }); }, function() { StackExchange.helpers.removeMessages(); } ); }); }); Jul 8 '18 at 14:17

• Taylor series approximations only "work" in a neighborhood of some point---they only converge on some open disk centered at the point where the expansion is based. Assuming that a function has enough derivatives, the approximation can be made as good as possible on that disk, and will blow up away from that disk. – Xander Henderson Jul 8 '18 at 3:00
• See e.g. MSE 2125155. You may also wish to read up on analytic functions (you want functions that are not analytic...). Wikilink: en.wikipedia.org/wiki/Analytic_function – Benjamin Dickman Jul 8 '18 at 3:01
• Sorry I did search the site before posting by I think my wording was poor. I will refer to the other post. – Hanzy Jul 8 '18 at 3:02
• I guess I took differentiability for granted. Maybe that was the point of contention by the other student, but it wasn’t stated specifically. – Hanzy Jul 8 '18 at 3:08
• @Hanzy In the case of a better search you may also want to see the "related posts " on the right side of the web page or else you can search your questions using approach0.xyz , hope this will help u search better questions again. – BAYMAX Jul 8 '18 at 3:30

Notice that to have a taylor series, our function must be smooth (e.g. infinitely continuously differentiable). Hence, a counter example will be any function that isn't smooth. Indeed, a simple counter example is \begin{align} f(x) = |x| \end{align} We cannot obtain a taylor series around $x=0$ because $f'(0)$ is not continuous.

For more sophiscated counter examples where $f$ is smooth, https://en.wikipedia.org/wiki/Non-analytic_smooth_function.

But a curious fact is that for a complex valued function $g$, then $g$ is smooth if and only if it's analytic(e.g. has a taylor series).

Here's a function that is smooth (infinitely differentiable) but fails to have a Taylor series at $x=0$:

Consider $f(x)=e^{-\frac{1}{x^2}}$ for $x \in \mathbb{R}-0$ and $f(0)=0$.

Show that $f^{(n)}=0$ and therefore the Taylor series for this function, if it existed, must be $0$. While this function is clearly not constant in any neighborhood of $0$.

In general, analytic functions have Taylor series. This is pretty much the definition of an analytic function as far as I'm aware.