Algorithm for product of vectors

Consider

uu=[
 [ ((1, 2), (3, -2)), ((1, 3), (3, -3))],
 [ ((1, -1), (2, 1))],
 [((1, 1), (2, -1)), ((1, -1), (2, 1))],
 [ ((2, -2), (4, 2)), ((2, -3), (4, 3))],
 [((3, 1), (4, -1)),  ((3, 4), (4, -4))],
 [((3, 1), (4, -1)), ((3, 2), (4, -2)), ((3, 3), (4, -3))]
]

and the following algorithm

function sum_exp(uu)
    sum_dict = Dict{Int, Int}()
    for pairs in uu
        for pair in pairs
            first_elem = pair[1]
            second_elem = pair[2]
            if haskey(sum_dict, first_elem)
                sum_dict[first_elem] += second_elem
            else
                sum_dict[first_elem] = second_elem
            end
        end
    end
    return sum_dict
end

function cartesian_product
(uu)
    res = []
    g = Iterators.product(uu...)
    for gi in g
        if all(values(sum_exp(gi)) .== 0)
            push!(res, gi)
        end
    end
    return res
end

This code works fine but it is slow and take a lot of memory.
I suspect the root cause lies in the execution of g = Iterators.product(uu...) prior to calculating sum_exp(gi)) .== 0 within the cartesian_product function

Code Functionality:

  • Each tuple within uu represents a polynomial, with ((1, 2), (3, -2)) corresponding to x1^2 * x3^(-2) .
  • The code aims to multiply each tuple in uu[1] by all tuples within uu[i] for i > 1 .
  • The sum_exp function returns terms with zero exponents.

what is the range for the integers?, in particular, i’m interested if the range for the values of the first elements are positive (and if possible, less than 100)

edit: nvm i saw the polynomial meaning

without tacking the Iterators.product problem or the res = []::Vector{Any} container, this implementation seems to be a little bit faster

function sum_exp2!(res,uu)
    for pairs in uu
        for pair in pairs
            p1,p2 = pair
            len = length(res)
            if p1 > len
                resize!(res,p1) #resizes vector to new size
                for i in (len + 1):p1
                    res[i] = 0
                end
            end
            res[p1] += p2
        end
    end
    return res
end

function cartesian_product2(uu)
    sum_exp_cache = Int[]
    res = []
    g = Iterators.product(uu...)
    for gi in g
        sum_exp_cache .= 0
        sum_exp2!(sum_exp_cache,gi)
        if all(iszero,sum_exp_cache)
            push!(res, gi)
        end
    end
    return res
end

2 Likes

Making them positive can them fast?
first element are variables, so positive and second elements can be make positive by adding the maximum value to all of them. the range is less than 100

using SpliApplyCombine
filterview(  x->all(==(0),x),collect(groupsum(first,last ,flatten(gi)) for gi in  Iterators.product(uu...)))
@btime Iterators.filter(x->all(==(0),groupsum(first,last ,flatten(x))), Iterators.product(uu...))
  302.024 ns (4 allocations: 256 bytes)
Base.Iterators.Filter{var"#31#32", Base.Iterators.ProductIterator{NTuple{6, Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}}}}(var"#31#32"(), Base.Iterators.ProductIterator{NTuple{6, Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}}}(([((1, 2), (3, -2)), ((1, 3), (3, -3))], [((1, -1), (2, 1))], [((1, 1), (2, -1)), ((1, -1), (2, 1))], [((2, -2), (4, 2)), ((2, -3), (4, 3))], [((3, 1), (4, -1)), ((3, 4), (4, -4))], [((3, 1), (4, -1)), ((3, 2), (4, -2)), ((3, 3), (4, -3))])))

if you need to collect

julia> collect(Iterators.filter(x->all(==(0),groupsum(first,last ,flatten(x))), Iterators.product(uu...)))
1-element Vector{NTuple{6, Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}}:
 (((1, 2), (3, -2)), ((1, -1), (2, 1)), ((1, -1), (2, 1)), ((2, -2), (4, 2)), ((3, 1), (4, -1)), ((3, 1), (4, -1)))

This give error

MethodError: no method matching count(::Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, ::Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, ::Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, ::Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, ::Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, ::Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}})

Closest candidates are:
  count(::Any, ::Union{Base.AbstractBroadcasted, AbstractArray}; dims, init)
   @ Base reducedim.jl:440
  count(::Any, ::Any; init)
   @ Base reduce.jl:1350
  count(::Union{Base.AbstractBroadcasted, AbstractArray}; dims, init)
   @ Base reducedim.jl:439
  ...


Stacktrace:
 [1] product(::Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, ::Vararg{Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}})
   @ Main ./In[78]:6
 [2] top-level scope
   @ In[85]:2

It’s not clear what you’re referring to

Thank you for your reply. I am using julia 1.9.1 and 1.9.3, but i got this error.
I try with the suggested name SplitApplyCombine but still getting the above error.

The following package names could not be resolved:
 * SpliApplyCombine (not found in project, manifest or registry)
   Suggestions: SplitApplyCombine

Stacktrace:
  [1] pkgerror(msg::String)
    @ Pkg.Types ~/Applications/Julia-1.9.1.app/Contents/Resources/julia/share/julia/stdlib/v1.9/Pkg/src/Types.jl:69

You are right.
My mistake, sorry,
The exact name is Split Apply Combine
You need to add the package.

import Pkg
Pkg.add("SplitApplyCombine")

julia> versioninfo()
Julia Version 1.9.0
Commit 8e63055292 (2023-05-07 11:25 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 8 × 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, tigerlake)
  Threads: 4 on 8 virtual cores
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 4

julia> 

julia> filterview(  x->all(==(0),x),collect(groupsum(first,last ,flatten(gi)) for gi in  Iterators.product(uu...)))
1-element view(::Array{Dictionaries.Dictionary{Int64, Int64}, 6}, CartesianIndex{6}[CartesianIndex(1, 1, 2, 1, 1, 1)]) with eltype Dictionaries.Dictionary{Int64, Int64}:
 {1 = 0, 3 = 0, 2 = 0, 4 = 0}

Regarding your output, it doesn’t tell which vector gives {1 = 0, 3 = 0, 2 = 0, 4 = 0}.
Also running the code gives error

MethodError: no method matching count(::Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, ::Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, ::Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, ::Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, ::Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, ::Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}})

Closest candidates are:
  count(::Any, ::Union{Base.AbstractBroadcasted, AbstractArray}; dims, init)
   @ Base reducedim.jl:440
  count(::Any, ::Any; init)
   @ Base reduce.jl:1350
  count(::Union{Base.AbstractBroadcasted, AbstractArray}; dims, init)
   @ Base reducedim.jl:439
  ...


Stacktrace:
 [1] product(::Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, ::Vararg{Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}})
   @ Main ./In[78]:6
 [2] top-level scope
   @ In[114]:3

An alternative method using properties of arithmetic and primes:

julia> uu=[
        [ ((1, 2), (3, -2)), ((1, 3), (3, -3))],
        [ ((1, -1), (2, 1))],
        [((1, 1), (2, -1)), ((1, -1), (2, 1))],
        [ ((2, -2), (4, 2)), ((2, -3), (4, 3))],
        [((3, 1), (4, -1)),  ((3, 4), (4, -4))],
        [((3, 1), (4, -1)), ((3, 2), (4, -2)), ((3, 3), (4, -3))]
       ];

julia> using .Iterators: flatten

julia> maximum(first, flatten(flatten(uu)))
4

julia> using Primes

julia> pnums = nextprimes(2,N)
4-element Vector{Int64}:
 2
 3
 5
 7

julia> uup = [[prod((pnums[b]//1)^e for (b,e) in pp) for pp in u] for u in uu]
6-element Vector{Vector{Rational{Int64}}}:
 [4//25, 8//125]
 [3//2]
 [2//3, 3//2]
 [49//9, 343//27]
 [5//7, 625//2401]
 [5//7, 25//49, 125//343]

julia> indices = CartesianIndices(Tuple(length.(uu)));

julia> qq = findall(t->prod(getindex.(uup, Tuple(t)))==1,indices)
1-element Vector{CartesianIndex{6}}:
 CartesianIndex(1, 1, 2, 1, 1, 1)

julia> q = Tuple(only(qq))
(1, 1, 2, 1, 1, 1)

julia> getindex.(uu,q)
6-element Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}:
 ((1, 2), (3, -2))
 ((1, -1), (2, 1))
 ((1, -1), (2, 1))
 ((2, -2), (4, 2))
 ((3, 1), (4, -1))
 ((3, 1), (4, -1))
2 Likes

the opposite of minimum

prime numbers are always beautiful

1 Like

To “extract” the information you could do this:

julia> uu=[
        [ ((1, 2), (3, -2)), ((1, 3), (3, -3))],
        [ ((1, -1), (2, 1))],
        [((1, 1), (2, -1)), ((1, -1), (2, 1))],
        [ ((2, -2), (4, 2)), ((2, -3), (4, 3))],
        [((3, 1), (4, -1)),  ((3, 4), (4, -4))],
        [((3, 1), (4, -1)), ((3, 2), (4, -2)), ((3, 3), (4, -3))] 
       ]
6-element Vector{Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}}:
 [((1, 2), (3, -2)), ((1, 3), (3, -3))]
 [((1, -1), (2, 1))]
 [((1, 1), (2, -1)), ((1, -1), (2, 1))]
 [((2, -2), (4, 2)), ((2, -3), (4, 3))]
 [((3, 1), (4, -1)), ((3, 4), (4, -4))]
 [((3, 1), (4, -1)), ((3, 2), (4, -2)), ((3, 3), (4, -3))]        

julia> using SplitApplyCombine
julia> idx=filterview(  x->all(==(0),x),collect(groupsum(first,last ,flatten(gi)) for gi in  Iterators.product(uu...)))
1-element view(::Array{Dictionaries.Dictionary{Int64, Int64}, 6}, CartesianIndex{6}[CartesianIndex(1, 1, 2, 1, 1, 1)]) with eltype Dictionaries.Dictionary{Int64, Int64}:
 {1 = 0, 3 = 0, 2 = 0, 4 = 0}

julia> [getindex.(uu,Tuple(t)) for t in idx.indices[1]]
1-element Vector{Vector{Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}}:
 [((1, 2), (3, -2)), ((1, -1), (2, 1)), ((1, -1), (2, 1)), ((2, -2), (4, 2)), ((3, 1), (4, -1)), ((3, 1), (4, -1))]

In any case it doesn’t seem so important to me since there is the modified version which directly provides the requested result.
Regarding the error, you should start a clean session of Julia and also post (your) code that produces the error
I take this opportunity to give a different modality that does not involve any dependencies

function cartesian_product(uu)
        res = []
        g = Iterators.product(uu...)
        for gi in g
            all(==(0),values(mergewith(+,Dict.(gi)...))) &&  push!(res, gi)
        end
        return res
    end

cartesian_product(uu)

#or 
julia> [  gi for gi in Iterators.product(uu...) if all(==(0),values(mergewith(+,Dict.(gi)...)))]
1-element Vector{NTuple{6, Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}}:
 (((1, 2), (3, -2)), ((1, -1), (2, 1)), ((1, -1), (2, 1)), ((2, -2), (4, 2)), ((3, 1), (4, -1)), ((3, 1), (4, -1)))
1 Like

Very beautifully done, I’am still wondering why it is working. (any link or paper?)
Thank you.

I’ll try to give an idea (not a formal proof) of the principles on which Dan’s proposal is based.
Instead of aggregating the 2-tuples with the same 1-element and sum of the values of the 2-element, a prime number (different for each 1-element) raised to the sum of the 2-elements is substituted.
Imposing that the product of p_k^{(v1k+v2k+v3k)}*p_h^{(v1h+v2h)}*p_j^{(v1j+v2j+v3j)} = 1 is equivalent to imposing that all exponents are equal to 0.

If the sum of the exponents is zero then each power is equal to 1 and the product is 1.
If the product is 1 then since p_k,p_h and p_j are prime then the exponents must be 0.

ps
It is not clear to me why the form (1/p)^(-e) is used instead of the equivalent form p^e.

2 Likes

(p//1)^e can be used. The (1//p) was just leftover from error when doing p^e, but the prime must be converted to a Rational first (fixed it in the solution).

Also note, that Int64 are can hold exponents up to a certain size before overflowing which limits this method. In the case more range is needed, (Int128(1)//p) can be used, or even BigInt.

Yes, certainly. I had verified that (p//1)^e remained of the rational type. :smile:
For this reason I asked myself whether there was another reason for the form (1//p)^(-e).
As for the need to try to make this solution “robust”, I wouldn’t pose the problem so much as, despite being surprisingly original, I would hardly use it as a solution to the problem since in addition to the initial transformation it requires an inverse transformation at the end to obtain the solution vector ( in this case of only one element) and not the vector of indices as expressly requested by the OP.

actually, (p//1)^e is the same as (1//p)^(-e), so the resulting fractions are identical. Even if they were not, only the indices of the ==1 elements are needed, the indices and ==1 do not depend on an inverse transformation.
As for speed, the 1//p rational transformation, is basically a no-op, as it simply stores the numerator and denominator. For the usual scalar numbers division is an expensive computation.

1 Like

You are right, I chose it because it gave a different idea and I think I can improve it. Even though @longemen3000 code is the fastest one. using @time.