Algorithm for product of vectors

blociss · February 10, 2024, 10:03pm

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.

longemen3000 · February 10, 2024, 11:57pm

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

longemen3000 · February 11, 2024, 1:02am

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

blociss · February 11, 2024, 9:45am

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

rocco_sprmnt21 · February 11, 2024, 11:52am

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)))

blociss · February 11, 2024, 12:26pm

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

rocco_sprmnt21 · February 11, 2024, 2:22pm

It’s not clear what you’re referring to

blociss · February 11, 2024, 3:47pm

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

rocco_sprmnt21 · February 11, 2024, 3:49pm

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}

blociss · February 11, 2024, 4:01pm

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

Dan · February 11, 2024, 6:19pm

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))

rocco_sprmnt21 · February 11, 2024, 9:58pm

the opposite of minimum

rocco_sprmnt21 · February 11, 2024, 10:10pm

prime numbers are always beautiful

rocco_sprmnt21 · February 12, 2024, 7:56am

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)))

blociss · February 12, 2024, 1:32pm

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

rocco_sprmnt21 · February 12, 2024, 2:53pm

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.

Dan · February 12, 2024, 4:56pm

(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.

rocco_sprmnt21 · February 12, 2024, 5:21pm

Yes, certainly. I had verified that (p//1)^e remained of the rational type.
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.

Dan · February 12, 2024, 5:43pm

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.

blociss · February 12, 2024, 10:06pm

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.