Any way to add outer while loops programmatically?

New to Julia
1

I have never come across this need before. I think maybe with meta-programming it is possible, but I have never done that before.

I have a function which takes an integer input nknots. For every increase in nknots by one, I need to add one more outer while loop to the code than was there for the original value. (I also need to add one more line of code per nknot with the current version, but I think this could be handled by an interior for loop.) So far, I have been hard-coding a pseudo-method for every possible value of the integer behind if statements. However, I would like to make the code generic so that it can take any integer value larger than 2. The code looks very uniform with all the hard-coded numbers easy to turn into n’s, but I don’t know how to add the outer while loops in a loop.

I need help replacing the portion of the code below inside the if statement with something generic. If nknots==3, only the innermost while loop should be used. If nknots==4, only the inner two while loops should be used. If nknots==6, an outer while loop should be added resembling the others. Ad infinitum.

using ElasticArrays
# Computes all possible knot locations on a dataset with `il` number of points, for a minimum segment length of `ilsegmin` and `nknots` number of knots.
function allknots(il::Int,ilsegmin::Int,nknots::Int)

    # Initialize first knot combination
    knots=ElasticArray{Int}(undef,nknots,1) # Create data storage structure for knots
    j=1 # Possible knot location counter
    knots[1]=1
    for n=2:nknots-1
        knots[n]=knots[n-1]+ilsegmin-1
    end
    knots[nknots]=il
    
if nknots==5

    # Determine all remaining knot combinations
    while knots[nknots-3,j] < il-3*ilsegmin+3 # Move thrid from last knot incrementally
        while knots[nknots-2,j] < il-2*ilsegmin+2 # Move second to last knot incrementally
            while knots[nknots-1,j] < il-ilsegmin+1 # Move last knot incrementally
                j+=1
                append!(knots,knots[:,j-1])
                knots[nknots-1,j]=knots[nknots-1,j]+1
            end
            j+=1
            append!(knots,knots[:,j-1])
            knots[nknots-2,j]=knots[nknots-2,j]+1
            knots[nknots-1,j]=knots[nknots-2,j]+ilsegmin-1
        end
        j+=1
        append!(knots,knots[:,j-1])
        knots[nknots-3,j]=knots[nknots-3,j]+1
        knots[nknots-2,j]=knots[nknots-3,j]+ilsegmin-1
        knots[nknots-1,j]=knots[nknots-2,j]+ilsegmin-1
    end

end

return knots
end
allknots(12,3,5)

5×20 ElasticArray{Int64,2,1,Array{Int64,1}}:
  1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
  3   3   3   3   3   3   3   3   3   3   4   4   4   4   4   4   5   5   5   6
  5   5   5   5   6   6   6   7   7   8   6   6   6   7   7   8   7   7   8   8
  7   8   9  10   8   9  10   9  10  10   8   9  10   9  10  10   9  10  10  10
 12  12  12  12  12  12  12  12  12  12  12  12  12  12  12  12  12  12  12  12
2

The example is a bit complicated, so I’m not sure this would work, but have you thought about a function with a single while loop that is called recursively? Roughly:

function while_knots(nknots, j, k)
  if k == 1
    # the j+=1 and append code goes here
  else
    while_knots(nknots, j, k-1)
  end
end
4 Likes
3

@Nathan_Boyer using Base: @nexprs, @nloops, @nref may be useful to what you are trying to achieve. Check the documentation of each macro and you may be able to rewrite N-dimensional code more easily.

1 Like
4

I am not sure I fully understand the example, but if you can map to iterating over a CartesianIndices, which you generate with ntuple, I would go for that pattern. Otherwise, lispy recursion (see the implementation of CartesianIndices again), or the macros suggested by @juliohm.

1 Like
5

I updated the while loops so that they are all exactly the same now except for a new integer k. I also added a comment to each line to explain what is being done. I think @hendri54’s suggestion could work now, but I haven’t had time yet to take it that far. I will look into the other suggestions as well, but that will take longer since they are a bit over my head at the moment.

if nknots==5
    # Determine all remaining knot combinations
    k=1 # Initialize k
    while knots[nknots-k,j] < il-k*ilsegmin+k         # Move thrid from last knot incrementally
        while knots[nknots-k,j] < il-k*ilsegmin+k     # Move second to last knot incrementally
            while knots[nknots-k,j] < il-k*ilsegmin+k # Move last knot incrementally
                k=1                                   # Current knot
                j+=1                                  # Increment possible knot counter
                append!(knots,knots[:,j-1])           # New knot combination same as last knot combination except...
                knots[nknots-k,j]+=1                  #  increment current knot and
                for kk=k-1:-1:1                       #  reset downstream knots
                    knots[nknots-kk,j]=knots[nknots-(kk+1),j]+ilsegmin-1
                end
            end
            k=2
            j+=1
            append!(knots,knots[:,j-1])
            knots[nknots-k,j]+=1
            for kk=k-1:-1:1
                knots[nknots-kk,j]=knots[nknots-(kk+1),j]+ilsegmin-1
            end
        end
        k=3
        j+=1
        append!(knots,knots[:,j-1])
        knots[nknots-k,j]+=1
        for kk=k-1:-1:1
            knots[nknots-kk,j]=knots[nknots-(kk+1),j]+ilsegmin-1
        end
    end
end
6

Full code below working nicely now for any number of knots with @hendri54’s method:

using ElasticArrays
# Computes all possible knot locations on a dataset with `il` number of points, for a minimum segment length of `ilsegmin` and `nknots` number of knots.
function allknots(il::Int, ilsegmin::Int, nknots::Int)

    # Initialize first knot combination
    knots = ElasticArray{Int}(undef, nknots, 1) # Create data storage structure for knots
    j = 1 # Possible knot location counter
    knots[1] = 1
    for n = 2:nknots - 1
        knots[n] = knots[n - 1] + ilsegmin - 1
    end
    knots[nknots] = il

    # Ensure first knot combination is possible
    if (nknots < 2) || (knots[nknots] - knots[nknots - 1] < ilsegmin)
        println("No possible knots with the given parameters.")
        return knots = ElasticArray{Int}(undef, nknots, 0)
    end
    
    # Determine all remaining knot combinations
    function incrementknot(il::Int, ilsegmin::Int, nknots::Int, k::Int=nknots - 2)

        while knots[nknots - k, j] < il - k * ilsegmin + k # While current knot not yet at maximum location
            if k > 1
                incrementknot(il, ilsegmin, nknots, k - 1) # Call function recursively to obtain inner while loops
            elseif k < 1
                throw(ArgumentError("k must be a counting number"))
            end
            j += 1                                         # Increment possible knot counter
            append!(knots, knots[:, j - 1])                # Create new knot combination same as last knot combination except...
            knots[nknots - k, j] += 1                      #  increment current knot and
            for kk = k - 1:-1:1                            #  reset all downstream knots
                knots[nknots - kk, j] = knots[nknots - (kk + 1), j] + ilsegmin - 1
            end
        end
        return knots
    end

    incrementknot(il, ilsegmin, nknots)
    return knots
end

# Test
for n = 1:7
    display(allknots(12, 3, n))
end