In [1]:
import numpy as np
import matplotlib.pyplot as pp
In [2]:
N_g = int( round( np.random.uniform( low = 3, high = 7 ) ) )
print( N_g )

4
In [3]:
groups = []

r = 5
dr = r / 5

theta = 2 * np.pi / ( N_g - 1 )
dtheta = theta / 5

for i in range( N_g ):
    
    if i == 0:
        center = np.array( [ 0.0, 0.0 ] )
    else:
        curr_r = np.random.uniform( low = r - dr, high = r + dr )
        
        curr_theta = np.random.uniform(
            low  = ( i - 1 ) * theta - dtheta,
            high = ( i - 1 ) * theta + dtheta
        )
        
        center = curr_r *  np.array( [ np.cos( curr_theta ), np.sin( curr_theta ) ] )
    
    points = []
    for j in range( int( round( np.random.uniform( low = 5, high = 15 ) ) ) ):
        
        curr_r = np.random.uniform( low = 0, high = r / 2 )
        
        curr_theta = np.random.uniform(
            low  = 0,
            high = 2 * np.pi
        )
        
        p_center = center + curr_r *  np.array( [ np.cos( curr_theta ), np.sin( curr_theta ) ] )
        v = np.random.uniform( low = 2, high = 4 )
        h = np.random.uniform( low = 2, high = 4 )
        points.append( ( ( h, v ), p_center ) )
        
    groups.append( ( center, points ) )
In [4]:
pp.figure( figsize = ( 10, 10 ) )

colors = [ 'blue', 'red', 'black', 'cyan', 'green', 'orange', 'gray' ]

i = 0
for c in groups:
#     pp.scatter( *c[ 0 ], marker = "x", color = colors[ i ] )
    
    for p in c[ 1 ]:
        
        pp.gca().add_patch(
            pp.Rectangle(
                ( p[ 1 ][ 0 ], p[ 1 ][ 1 ] ),
                p[ 0 ][ 0 ] , p[ 0 ][ 1 ],
                color = colors[ i ],
                alpha = 0.1
            )
        )
#         pp.plot(
#             [ 
#                 p[ 1 ][ 0 ],
#                 p[ 1 ][ 0 ] + p[ 0 ][ 0 ],
#                 p[ 1 ][ 0 ] + p[ 0 ][ 0 ],
#                 p[ 1 ][ 0 ]                
#             ],
#             [ 
#                 p[ 1 ][ 1 ],
#                 p[ 1 ][ 1 ],
#                 p[ 1 ][ 1 ] + p[ 0 ][ 1 ],
#                 p[ 1 ][ 1 ] + p[ 0 ][ 1 ]                
#             ],
#             color = colors[ i ]
#         )
#         pp.scatter( *p[ 1 ], color = colors[ i ] )
        
    i = i + 1

pp.axis( ( -10, 10, -10, 10 ) )
    
pp.show()
In [5]:
def zero_trunc( x ):
    if x >= 0.0:
        return x
    else:
        return 0.0

def O_ij( x, h, i, j ):
    
    if x[ i ] >= x[ j ]:
        
        return ( zero_trunc( h[ j ] ** 2.0 - ( x[ i ] - x[ j ] ) ** 2.0 ) ** 2.0 ) / ( h[ j ] ** 4.0 )
    
    else:
        return ( zero_trunc( h[ i ] ** 2.0 - ( x[ i ] - x[ j ] ) ** 2.0 ) ** 2.0 ) / ( h[ i ] ** 4.0 )

def E_O( xy, h, v ):
    
    x = xy[ : : 2 ]
    y = xy[ 1 : : 2 ]
    n = int( xy.shape[ 0 ] / 2 )
    
    retval = 0.0
    
    for i in range( n ):
        
        curr_xi = x[ i ]
        curr_yi = y[ i ]
        
        curr_hi = h[ i ]
        curr_vi = v[ i ]
        
        for j in range( i + 1, n ):
            
            curr_xj = x[ j ]
            curr_yj = y[ j ]

            curr_hj = h[ j ]
            curr_vj = v[ j ]
            
            retval += O_ij( x, h, i, j ) * O_ij( y, v, i, j )
            
    return 2.0 * retval / ( n * ( n + 1 ) )
In [6]:
def gera_xyhv( groups ):
    
    xy = []
    h = []
    v = []
    
    for g in groups:
        
        for p in g[ 1 ]:
            
            h.append( p[ 0 ][ 0 ] )
            v.append( p[ 0 ][ 1 ] )
            xy.append( p[ 1 ][ 0 ] )
            xy.append( p[ 1 ][ 1 ] )
            
    return( np.array( xy ), np.array( h ), np.array( v ) )
In [7]:
( xy, h, v ) = gera_xyhv( groups )
In [8]:
print( E_O( xy, h, v ) )

0.08946818762626481

$$ \sum_{i = 0}^{n - 1}\sum_{j = i + 1}^{n - 1}1 = \sum_{i = 0}^{n - 1}n - 1 - ( i + 1 ) + 1 = \sum_{i = 0}^{n - 1}n - i - 1 $$

Seja $j = n - 1 - i$, então: $$ \sum_{i = 0}^{n - 1}\sum_{j = i + 1}^{n - 1}1 = \sum_{j = n - 1}^0 j = \sum_{j = 0}^{n - 1} j $$

In [9]:
def calcula_L( groups ):
    
    n_verts = 0
    for g in groups:
        for p in g[ 1 ]:
            n_verts += 1
            
    L = np.eye( n_verts )
    
    curr_idx = 0
    for g in groups:
        
        curr_n = len( g[ 1 ] )
        
        for i in range( curr_n ):
            for j in range( i + 1, curr_n ):
                    L[ curr_idx + i, curr_idx + j ] = -1.0 / ( curr_n - 1.0 )
                    L[ curr_idx + j, curr_idx + i ] = L[ curr_idx + i, curr_idx + j ]
        
        curr_idx += curr_n
        
    return L
    
L = calcula_L( groups )
In [10]:
delta_x = L @ xy[ : : 2 ]
delta_y = L @ xy[ 1 : : 2 ]
In [11]:
def E_N( xy, w ):
    
    retval = np.linalg.norm( L @ xy[ : : 2 ] - w * delta_x ) ** 2
    retval += np.linalg.norm( L @ xy[ 1 : : 2 ] - w * delta_y ) ** 2
    
    retval *= ( L.shape[ 0 ] ** 2 ) / \
              ( 2.0 * ( np.linalg.norm( delta_x ) ** 2 ) + np.linalg.norm( delta_y ) ** 2 )
        
    return retval
In [12]:
print( E_N( xy / 2.0, 0.5 ) )

0.0
In [13]:
import scipy.optimize
In [14]:
res = scipy.optimize.minimize( E_O, xy, ( h, v ) )
In [15]:
print( res.x )

[ -42.30722678   85.97017328    2.64308081  -61.81827828  -26.32355969
   25.20718784  -80.95467565  -16.77509433  -55.05285084  -11.01288943
   31.74369084  -14.19879609   14.07567587  -70.05793448  -38.30556668
    1.0031023    12.17971643    2.70089429    0.36339466 -100.43271758
   26.54648417  -45.7806287    63.9845504    44.27708055   41.23762484
   -2.97389535   63.53371268   50.81203328    6.70768413  -87.73217062
    0.76025928   56.2482034    45.39588912  -38.02015855  -51.41956357
   46.70444719   68.73070703   49.83845741   15.94543022    9.99950355
   42.4898151    57.25410262  -74.3451693    29.88713488  -59.84903865
    7.88965919    1.12891111  -14.38141334   11.03887298   68.07765425
   -9.34765233  -42.28738927  -83.15027197  -28.65493445   39.74474273
    2.2659919    45.01439826   45.50940431  -46.28602389  -24.16996902
  113.37234714   27.91901023  -32.34259225  -49.32260423   94.42799472
   30.584592    -42.90923333   -4.04346519  110.41113588  -28.54397655
  -90.23233522   41.43820173  -50.69651372  -23.75184769    3.78423378
   54.78412852  -52.97849242   10.018625     -9.85524112   55.39515771
   70.69621289  -79.70259997  -59.74042001  -78.35325845  -45.0011785
  -11.01334338]
In [16]:
pp.figure( figsize = ( 10, 10 ) )

colors = [ 'blue', 'red', 'black', 'cyan', 'green', 'orange', 'gray' ]

i = 0
j = 0
for c in groups:
    
    for p in c[ 1 ]:
        
        pp.gca().add_patch(
            pp.Rectangle(
                ( res.x[ 2 * j ], res.x[ 2 * j + 1 ] ),
                p[ 0 ][ 0 ] , p[ 0 ][ 1 ],
                color = colors[ i ],
                alpha = 1.0
            )
        )
        
        j = j + 1
        
    i = i + 1

pp.axis( ( -150, 150, -150, 150 ) )
    
pp.show()
In [17]:
print( E_O( res.x, h, v ) )

0.0
In [18]:
def obj_fun( z, alpha ):
    # z contém xy e w
    
    xy = z[ : -1 ]
    w = z[ -1 ]
    
    return alpha * E_N( xy, w ) + ( 1.0 - alpha ) * E_O( xy, h, v )
In [19]:
res = scipy.optimize.minimize( obj_fun, np.concatenate( ( xy, np.ones( ( 1, ) ) ) ), ( 0.9 ) )
In [20]:
pp.figure( figsize = ( 10, 10 ) )

colors = [ 'blue', 'red', 'black', 'cyan', 'green', 'orange', 'gray' ]

i = 0
j = 0
for c in groups:
    
    for p in c[ 1 ]:
        
        pp.gca().add_patch(
            pp.Rectangle(
                ( res.x[ 2 * j ], res.x[ 2 * j + 1 ] ),
                p[ 0 ][ 0 ] , p[ 0 ][ 1 ],
                color = colors[ i ],
                alpha = 1.0
            )
        )
        
        j = j + 1
        
    i = i + 1

pp.axis( ( -100, 100, -100, 100 ) )
    
pp.show()