#! /usr/bin/env python # def mono_unrank_grlex ( m, rank ): #*****************************************************************************80 # ## MONO_UNRANK_GRLEX computes the monomial of given grlex rank. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 25 October 2015 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, the spatial dimension. # 1 <= M. # # Input, integer RANK, the rank of the monomial. # # Output, integer X[M], the monomial. # import numpy as np from i4_choose import i4_choose from sys import exit # # Ensure that 1 <= M. # if ( m < 1 ): print ( '' ) print ( 'MONO_UNRANK_GRLEX - Fatal error!' ) print ( ' M < 1' ) exit ( 'MONO_UNRANK_GRLEX - Fatal error!' ) # # Ensure that 1 <= RANK. # if ( rank < 1 ): print ( '' ) print ( 'MONO_UNRANK_GRLEX - Fatal error!' ) print ( ' RANK < 1' ) exit ( 'MONO_UNRANK_GRLEX - Fatal error!' ) x = np.zeros ( m, dtype = np.int32 ) # # Special case M = 1. # if ( m == 1 ): x[0] = rank - 1 return x # # Determine the appropriate value of NM. # Do this by adding up the number of compositions of sum 0, 1, 2, # ..., without exceeding RANK. Moreover, RANK - this sum essentially # gives you the rank of the composition within the set of compositions # of sum NM. And that's the number you need in order to do the # unranking. # rank1 = 1 nm = -1 while ( True ): nm = nm + 1 r = i4_choose ( nm + m - 1, nm ) if ( rank < rank1 + r ): break rank1 = rank1 + r rank2 = rank - rank1 # # Convert to KSUBSET format. # Apology: an unranking algorithm was available for KSUBSETS, # but not immediately for compositions. One day we will come back # and simplify all this. # ks = m - 1 ns = nm + m - 1 xs = np.zeros ( ks, dtype = np.int32 ) nksub = i4_choose ( ns, ks ) j = 1 for i in range ( 1, ks + 1 ): r = i4_choose ( ns - j, ks - i ) while ( r <= rank2 and 0 < r ): rank2 = rank2 - r j = j + 1 r = i4_choose ( ns - j, ks - i ) xs[i-1] = j j = j + 1 # # Convert from KSUBSET format to COMP format. # x[0] = xs[0] - 1 for i in range ( 2, m ): x[i-1] = xs[i-1] - xs[i-2] - 1 x[m-1] = ns - xs[ks-1] return x def mono_unrank_grlex_test ( ): #******************************************************************************/ # ## MONO_UNRANK_GRLEX_TEST tests MONO_UNRANK_GRLEX. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 25 October 2014 # # Author: # # John Burkardt # import numpy as np import platform from i4_uniform_ab import i4_uniform_ab from mono_upto_enum import mono_upto_enum from mono_upto_next_grlex import mono_upto_next_grlex m = 3 print ( '' ) print ( 'MONO_UNRANK_GRLEX' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' MONO_UNRANK_GRLEX is given a rank, and returns the corresponding' ) print ( ' monomial in the sequence of all monomials in M dimensions' ) print ( ' in grlex order.' ) print ( '' ) print ( ' For reference, print a monomial sequence with ranks.' ) n = 4 rank_max = mono_upto_enum ( m, n ) print ( '' ) print ( ' Let M = %d' % ( m ) ) print ( ' N = %d' % ( n ) ) print ( '' ) x = np.zeros ( m, dtype = np.int32 ) i = 1 while ( True ): print ( ' %2d ' % ( i ) ), for j in range ( 0, m ): print ( '%2d' % ( x[j] ) ), print ( '' ) if ( x[0] == n ): break mono_upto_next_grlex ( m, n, x ) i = i + 1 print ( '' ) print ( ' Now choose random ranks between 1 and %d' % ( rank_max ) ) print ( '' ) seed = 123456789 test_num = 5 for test in range ( 0, test_num ): rank, seed = i4_uniform_ab ( 1, rank_max, seed ) x = mono_unrank_grlex ( m, rank ) print ( ' %2d ' % ( rank ) ), for j in range ( 0, m ): print ( '%2d' % ( x[j] ) ), print ( '' ) # # Terminate. # print ( '' ) print ( 'MONO_UNRANK_GRLEX_TEST' ) print ( ' Normal end of execution.' ) return if ( __name__ == '__main__' ): from timestamp import timestamp timestamp ( ) mono_unrank_grlex_test ( ) timestamp ( )