#StackBounty: #algorithm #math #discrete-mathematics Counting integer tuples with specific property

Bounty: 50

Given two numbers M_max,M_min count all the r-tuples of prime numbers (M_1,…,M_r) such that:

  1. M_min < M_i < M_max
  2. every M_i-1 is divisible by at least one square of an odd prime
  3. lcm(M_1-1,…,M_r-1)=(M_1-1)*(M_2-1)…(M_r-1)/2^(r-1)

Attempt so far: I have a lookup table of the primes M_i such that M_i-1 is not odd-squarefree as requested by 2. Also condition 3. as @btilly pointed out is equivalent to the numbers n_i=(M_i-1)/2 being coprime which leads us to:

Rephrased Question: How many r-tuples among the k numbers M_1,M_2,…,M_k are coprime?

Remark: In the general case the problem is equivalent to counting the number of matchings in a bipartite graph defined by the prime divisors p_j connected to the corresponding M_i it divides but I am not sure if something about the particular structure helps simplify it.

Main Question: Can we at least do something when r is small?

We can certainly search all (n binom r) tuples in that case giving us O(n^r) but can we do better?


Get this bounty!!!

#StackBounty: #javascript #math #graph #3d #2d JavaScript 3D Chart Library for Math/Engineering. Projection of 3D figures to 2D

Bounty: 100

I’m looking for a JavaScript library for plotting 3D charts with figures like Cylinder, Box, Torus, Sphere etc for mathematical/engineering purposes.

The library I’m looking for should have ability to mark some important sections of the graph like function name, intersections, crossing points etc.

There is no need for any 3D shading to show a shape of a 3D figure. Instead the library should draw dashed lines in case they are covered by the front of a 3D figure.

I found a good 2D library for drawing geometry JSX Graph that is well suited for the purpose but unlucky it lacks drawing 2D projection of 3D figures.

Below some examples from JSX Graph to picture some of the features that are required.

enter image description here

Adding auxiliary lines as dashed, adding labels to graphs, marking an important joint point

enter image description here

Adding points of lines’ cross, marking an intersection area.

In short the library should be able to make 3D graphs like these you can find in books or a technical documentation.
Do you know any library like JSX graph but with a support of projecting 3D figures to the 2D?


Get this bounty!!!

#StackBounty: #python #numpy #math #pca #voxel Why my PCA is not invariant to rotation and axis swap?

Bounty: 50

I have a voxel (np.array) with size 3x3x3, filled with some values, this setup is essential for me. I want to have rotation-invariant representation of it. For this case, I decided to try PCA representation which is believed to be invariant to orthogonal transformations. another

For simplicity, I took some axes swap, but in case I’m mistaken there can be np.rot90.

I have interpereted my 3d voxels as a set of weighted 3d cube point vectors which I incorrectly called "basis", total 27 (so that is some set of 3d point in space, represented by the vectors, obtained from cube points, scaled by voxel values).

import numpy as np

voxel1 = np.random.normal(size=(3,3,3))
voxel2 =  np.transpose(voxel1, (1,0,2)) #np.rot90(voxel1) #


basis = []
for i in range(3):
    for j in range(3):
        for k in range(3):
            basis.append([i+1, j+1, k+1]) # avoid 0
basis = np.array(basis)


voxel1 = voxel1.reshape((27,1))
voxel2 = voxel2.reshape((27,1))

voxel1 = voxel1*basis # weighted basis vectors
voxel2 = voxel2*basis
print(voxel1.shape)
(27, 3)

Then I did PCA to those 27 3-dimensional vectors:

def pca(x):
    center = np.mean(x, 0)
    x = x - center

    cov = np.cov(x.T) / x.shape[0]

    e_values, e_vectors = np.linalg.eig(cov)

    order = np.argsort(e_values)

    v = e_vectors[:, order].transpose()

    return x.dot(v)

vp1 = pca(voxel1)
vp2 = pca(voxel2)

But the results in vp1 and vp2 are different. Perhaps, I have a mistake (though I beleive this is the right formula), and the proper code must be

x.dot(v.T)

But in this case the results are very strange. The upper and bottom blocks of the transofrmed data are the same up to the sign:

>>> np.abs(np.abs(vp1)-np.abs(vp2)) > 0.01
array([[False, False, False],
       [False, False, False],
       [False, False, False],
       [ True,  True,  True],
       [ True,  True,  True],
       [ True,  True,  True],
       [ True,  True,  True],
       [ True,  True,  True],
       [ True, False,  True],
       [ True,  True,  True],
       [ True,  True,  True],
       [ True,  True,  True],
       [False, False, False],
       [False, False, False],
       [False, False, False],
       [ True,  True,  True],
       [ True,  True,  True],
       [ True,  True,  True],
       [ True,  True,  True],
       [ True,  True,  True],
       [ True, False,  True],
       [ True,  True,  True],
       [ True,  True,  True],
       [ True,  True,  True],
       [False, False, False],
       [False, False, False],
       [False, False, False]])

What I’m doing wrong?

What I want to do is to find some invariant representation of my weighted voxel, something like positioning according to the axes of inertia or principal axes. I would really appreciate if someone helps me.

UPD: Found the question similar to mine, but code is unavailable

EDIT2: Found the code InertiaRotate and managed to monkey-do the following:

import numpy as np

# https://github.com/smparker/orient-molecule/blob/master/orient.py

voxel1 = np.random.normal(size=(3,3,3))
voxel2 =  np.transpose(voxel1, (1,0,2))

voxel1 = voxel1.reshape((27,))
voxel2 = voxel2.reshape((27,))


basis = []
for i in range(3):
    for j in range(3):
        for k in range(3):
            basis.append([i+1, j+1, k+1]) # avoid 0
basis = np.array(basis)
basis = basis - np.mean(basis, axis=0)



def rotate_func(data, mass):

    #mass = [ masses[n.lower()] for n in geom.names ]

    inertial_tensor = -np.einsum("ax,a,ay->xy", data, mass, data)
    # negate sign to reverse the sorting of the tensor
    eig, axes = np.linalg.eigh(-inertial_tensor)
    axes = axes.T

    # adjust sign of axes so third moment moment is positive new in X, and Y axes
    testcoords = np.dot(data, axes.T) # a little wasteful, but fine for now
    thirdmoment = np.einsum("ax,a->x", testcoords**3, mass)

    for i in range(2):
        if thirdmoment[i] < 1.0e-6:
            axes[i,:] *= -1.0

    # rotation matrix must have determinant of 1
    if np.linalg.det(axes) < 0.0:
        axes[2,:] *= -1.0

    return axes

axes1 = rotate_func(basis, voxel1)
v1 = np.dot(basis, axes1.T)
axes2 = rotate_func(basis, voxel2)
v2 = np.dot(basis, axes2.T)


print(v1)
print(v2)

It seems to use basis (coordinates) and mass separately. The results are quite similar to my problem above: some parts of the transformed data match up to the sign, I believe those are some cube sides

print(np.abs(np.abs(v1)-np.abs(v2)) > 0.01)

[[False False False]
 [False False False]
 [False False False]
 [ True  True  True]
 [ True  True  True]
 [ True  True  True]
 [ True  True  True]
 [False False False]
 [ True  True  True]
 [ True  True  True]
 [ True  True  True]
 [ True  True  True]
 [False False False]
 [False False False]
 [False False False]
 [ True  True  True]
 [ True  True  True]
 [ True  True  True]
 [ True  True  True]
 [False False False]
 [ True  True  True]
 [ True  True  True]
 [ True  True  True]
 [ True  True  True]
 [False False False]
 [False False False]
 [False False False]]

Looking for some explanation. This code is designed for molecules, and must work…

UPD: Tried to choose 3 vectors as a new basis from those 24 – the one with biggest norm, the one with the smallest and their cross product. Combined them into the matrix V, then used the formula V^(-1)*X to transform coordinates, and got the same problem – the resulting sets of vectors are not equal for rotated voxels.


Get this bounty!!!

#StackBounty: #neural-network #deep-learning #backpropagation #cost-function #math Chain function in backpropagation

Bounty: 100

I’m reading a Neural Networks Tutorial. In order to answer my question you might have to take a brief look at it.

I understand everything until something they declare as "chain function":

enter image description here

It’s located under 4.5 Backpropagation in depth.

I know that the chain rule says that the derivative of a composite function equals to the product of the inner function derivative and the outer function derivative, but still don’t understand how they got to this equation.

Update:

Please read the article and try to only use the terms and signs that are declared there. It’s confusing and tough enough anyway, so I really don’t need extra-knowledge now. Just use the same "language" as in the article and explain how they arrived to that equation.


Get this bounty!!!

#StackBounty: #code-golf #math #abstract-algebra Finding Radical Ideals

Bounty: 50

Sandbox

Finding radicals of arbitrary ideals in rings is hard! However, it’s a doable task in polynomial rings, since they are Noetherian (all ideals are finitely generated)

Background

A quick introduction to radical ideals

A ring is a set, paired with two binary operations + and * such that 0 is a additive identity, 1 is a multiplicative identity, addition is commutative, every element has an additive inverse, and multiplication distributes over addition. For a more proper introduction, see Wikipedia. In this challenge, we’ll be dealing with rings of polynomials over the complex numbers (for example, C[x,y,z]), so multiplication will also be commutative.

An ideal is a non-empty set closed under addition that contains additive inverses for all elements, and such that for any r in the ring and i in the ideal, r*i is in the ideal.

Given some polynomials p1,…,pn, the ideal generated by the polynomials is the set of all polynomials p such that there exists polynomials q1,…,qn such that p = p1q1+p2q2+…+pnqn.

In other words, the ideal generated by a set of polynomials is the smallest ideal containing the set.

Given ideal I of R, the radical ideal of I, or root I, √I, is the ideal of elements r in R such that there is some natural n, r^n is in I.

challenge:

Given an ideal generated by polynomials in any number of variables in the complex numbers, you need to return a set of polynomials that generate the radical ideal.

Important Notes:

  • While the polynomials all will belong to the space of complex coefficients, I’ll only give you inputs of polynomials with integer coefficients.

  • Variables will always be single letters.

  • The output must be algebraically independent (there is no non-zero linear combination of the polynomials that sums to 0). The input will be algebraically independent.

  • Obviously, the output is non-unique.

Hints

A less algebraic way of thinking about radical ideals is: given that the input polynomials all evaluate to 0, what is the set of all polynomials we know that evaluate to 0. That set is the radical ideal.

Another nice characterization of the radical ideal is the intersection of all prime ideals containing the original generators.

There are no good known ways to find radicals in infinite fields, like the complex numbers. However, there are several algorithms that have reasonable average case times (many of which are listed here)

Singular has all of these algorithms implemented, if you want to experiment before building your own.

Of course, this is code golf, so inefficient solutions are permitted.

Input/Output

I/O should be pretty flexible. Possible forms of input/output:

Ideal: (x^2y-z^3, 2xy^2-2w^3, xy-zw)

I/O 0: "x^2y-z^3, 2xy^2-2w^3, xy-zw"

I/O 1: "x2y-z3,2xy2-2w3,xy-zw"

I/O 2: ["x2y-z3","2xy2-2w3","xy-zw"]

I/O 3: [[[1,2,1,0,0],[-1,0,0,3,0]],[[2,1,2,0,0],[-2,0,0,0,3]],[[1,1,1,0,0],[-1,0,0,1,1]]]

If there’s some other more convenient I/O, let me know.

Test Cases

As I said earlier, these are not the only valid outputs for these inputs.

x^3 → x

2x^4y^2z^3 → xyz

xy-z^2, x-y → x-y, x^2-z^2

x^2, xz^2-x, y-z → x, y-z

y^2+z, yz → y, z

x^2+y^2-1, w^2+z^2-1, xw+yz, xz-yw-1 → y+w, x-z, y^2+z^2-1

xy, (x-y)z → xy, yz, xz

xy, x-yz → x, yz

x^2y-z^3, 2xy^2-2w^3, xy-zw → z^2-xw, yz-w^2, xy-zw

Explained Examples

Upon request in the comments, I’ll work a few examples out:

Example 1: 2x^4y^2z^3 → xyz

Note that (xyz)^4 = x^4y^4z^4, while 0.5*y^2z*(2x^4y^2z^3)=x^4y^4z^4, so this tells us that xyz is in the radical of 2x^4y^2z^3. Now, consider some polynomial p with a term that doesn’t have x as a factor. Then, p^n will also have such a term, so it can’t be in the ideal generated by 2x^4y^2z^3. So, every polynomial in the radical ideal has x as a divisor, and by symmetry, has y and z as divisors, so has xyz as a divisor.

Example 2: x^2, xz^2-x, y-z → x, y-z

Note that x^2 is in the ideal, so x is in the radical ideal. So, x*(z^2-1)=xz^2-x will be in the radical ideal since it is a multiple of x. Clearly, x and y-z are independent, and since they are terms of degree 1, they generate the radical ideal.


Get this bounty!!!

#StackBounty: #string #math #combinatorics #fastest-code Average number of strings with Levenshtein distance up to 4

Bounty: 50

This is a version of this question which should not have such a straightforward solution and so should be more of an interesting coding challenge. It seems, for example, very likely there is no easy to find closed form solution, even though we have only increased the bound by one from the previous version. Having said that, you never know…

The Levenshtein distance between two strings is the minimum number of single character insertions, deletions, or substitutions to convert one string into the other one. Given a binary string $S$ of length $n$, we are a interested in the number of different strings of length $n$ which have distance at most $4$ from $S$.

For example, if $S = 0000$ there are four strings with Levenshtein distance exactly $3$ from $S$, six with distance exactly $2$, four with distance exactly $1$ and exactly one with distance $0$. This makes a total of $15$ distinct strings with distance at most $3$ from the string $0000$. The only string with distance exactly $4$ is $1111$.

For this task the input is a value of $n geq 4$. Your code must output the average number of binary strings of length $n$ which have Levenshtein distance at most $4$ from a uniform and randomly sampled string $S$. Your answer can be output in any standard way you choose but it must be exact.

Examples

  • n = 4. Average $16$.
  • n = 5. Average 31 $frac{11}{16}$.
  • n = 6. Average 61 $frac{21}{32}$.
  • n = 7. Average 116 $frac{7}{8}$.
  • n = 8. Average 214 $frac{43}{128}$.
  • n = 9. Average 378 $frac{49}{246}$.
  • n = 10. Average 640 $frac{301}{512}$.
  • n = 11. Average 1042 $frac{1}{16}$.
  • n = 12. Average 1631 $frac{1345}{2048}$.
  • n = 13. Average 2466 $frac{3909}{4096}$.
  • n = 14. Average 3614 $frac{563}{8192}$

Score

Your score is the highest value of $n$ you can reach in less than a minute on my linux box. If two answers get the same value then the first posted (or amended) wins. The timing is for each value separately.

My CPU is an Intel(R) Xeon(R) CPU X5460.


Get this bounty!!!

#StackBounty: #math #pseudocode Evenly distributing points on the surface of hyperspheres in higher dimensions

Bounty: 250

I am interested in evenly distributing N points on the surface of spheres in dimensions 3 and higher.

To be more specific:

  • Given a number of points N and number of dimensions D (where D > 1, N > 1)
  • The distance of every point to the origin must be 1
  • The minimum distance between any two points should be as large as possible
  • The distance of each point to it’s closest neighbour doesn’t necessarily have to be the same for every point (indeed it’s not possible for it to be the same unless the number of points forms the vertices of a platonic solid or if N <= D).

I am not interested in:

  • Creating a uniform random distribution on a hypersphere, because I want the minimum distance between any two points to be as large as possible instead of being randomly distributed.
  • Particle repulsion simulation type methods, because they are hard to implement and take an extremely long time to run for large N (Ideally the method should be deterministic and in O(n)).

One method that satisfies these criteria is called the fibonacci lattice, but I have only been able to find code implementations for that in 2d and 3d.

The method behind the fibonacci lattice (also known as the fibonacci spiral) is to generate a 1d line that spirals around the surface of the sphere such that the surface area covered by the line is roughly the same at every turn. You can then drop N points equally distributed on the spiral and they will roughly be evenly distributed on the surface of the sphere.

In this answer there is a python implementation for 3 dimensions that generates the following:

enter image description here

I wanted to know whether the fibonacci spiral could be extended to dimensions higher than 3 and posted a question on the maths stack exchange. To my surprise I received two amazing answers which as far as I can tell (because I don’t fully understand the maths shown) shows it’s indeed possible to extend this method to N dimensions.

Unfortunately I don’t understand enough of the maths shown to be able to turn either answer into (pseudo)code. I am an experienced computer programmer, but my maths background only goes so far.

I will copy in what I believe to be the most important part of one of the answers below (unfortunately SO doesn’t support mathjax so I had to copy as an image)

enter image description here

Difficulties presented by the above that I struggle with:

  • How to resolve the inverse function used for ψn?
  • The example given is for d = 3. How do I generate formulae for arbitrary d?

Would anyone here who understands the maths involved be able to make progress towards a pseudo code implementation of either answer? I understand a full implementation may be quite difficult so I’d be happy with a part implementation that leads me far enough to be able to complete the rest myself.

To make it easier, I’ve already coded a function that takes spherical coordinates in N dimensions and turns them into cartesian coordinates, so the implementation can output either one as I can easily convert.

Additionally I see that one answer uses the next prime number for each additional dimension. I can easily code a function that outputs each successive prime, so you can assume that’s already implemented.

Failing an implementation of the fibonacci lattice in N dimensions. I’d be happy to accept a different method that satisfies the above constraints.


Get this bounty!!!