\(\(\(
\)\( (*This\ file\ contains\ all\ the\ programs\ needed\ to\ implement\ the\ algorithms\ for\ "Identifying Community Structures from Network Data via Maximum Likelihood Methods", \ Copic, \ 
      Jackson, \ and\ Kirman, \ 
      December\ 2006*) \)\(
\)\( (*In\ order\ to\ run\ these\ programs, \ one\ needs\ to\ have\ installed\ Mathematica, \ 
      version\ 5.2\ or\ later . \ 
          To\ run\ the\ programs\ select\ the\ Kernel\ from\ the\ menu, \ 
      and\ click\ on\ the\ Evaluation\  -> \ 
        Evaluate\ notebook . \ 
            The\ path\ to\ the\ data\ files\ needs\ to\ be\ correctly\ specified\ for\ the\ programs\ to\ \(\(work\)\(.\)\)*) \)\(
\)\( (*All\ commands\ are\ documented . \ 
          The\ exact\ greedy\ method\ to\ search\ all\ neighborhoods\ of\ a\ given\ node\ in\ order\ to\ determine\ the\ closest\ partition\ to\ the\ given\ partition, \ function\ "closest", \ 
      is\ included\ only\ for\ illustration, \ 
      as\ it\ can\ only\ be\ used\ on\ graphs\ of\ size\ up\ to\ approximately\ 15\ nodes\ due\ to\ its\ computational\ demands\ \(\((the\ problem\ is\ very\ computationally\ demanding\ to\ solve\ exactly)\)\(.\)\)*) \)\(\ \)\(
\)\( (*\ 
      Please\ email\ jernej . copic@yale . 
          edu\ if\ you\ have\ any\ questions\ or\ \(\(problems\)\(.\)\)*) \)\)\)

neighbors[s_List, x_] := 
      Flatten[Position[s[[x]], 1]]

     (*s is a 0 - 1 matrix of connections,  
      x is a node (row),  
      output is the vector of nodes to which node x is connected*)

Clear[components]
\(components[s_List] := 
      Module[{g = s, checked = {}, notchecked = {}, \ nodes, rest, comp, l, 
          x}, 
comp = {}; 
l = 
          Length[s]; 
nodes = 
          Table[i, {i, 1, l}]; 
rest = 
          nodes; 
While[Not[rest == {}], 
          x = rest[\([1]\)]; 
		checked = {}; 
		notchecked = {x}; \ 
		While[
            Not[notchecked == {}], 
			x = 
              notchecked[\([1]\)]; 
			checked = 
              Union[checked, {x}]; 
			notchecked = 
			Union[notchecked, 
                Complement[neighbors[g, x], 
                  checked]]; 
			notchecked = 
              Complement[
                notchecked, {x}]
		]; 
		checked = Union[checked, {x}]; 
		comp = 
            Union[comp, {checked}]; 
		rest = 
            Complement[rest, checked]]; 
Return[
          comp]]\)

Clear[mat]

    mat[s_List, matr_List] := 
      Part[Map[Part[#, s] &, matr], s]

     (*mat takes a list s and a matrix matr and returns matr_s^
          s *) 

     (*that is it returns rows and columns specified by s*)

isComplete[matr_List] := 
        matr == 
          Table[1, {i, 1, Length[matr]}, {i, 1, 
              Length[matr]}];

complete[s_List] := 
      Total[Table[1, {i, 1, Length[s]}, {i, 1, Length[s]}] - s, 
        2]

<< LinearAlgebra`MatrixManipulation`
diag[a_List, b_List] := 
      BlockMatrix[{{a, 
            Table[0, {i, Length[a]}, {j, Length[b]}]}, {Table[
              0, {i, Length[b]}, {j, Length[a]}], 
            b}}]

uncomplete[s_List, n_List] := 
      Total[Abs[
            mat[Join[n, Complement[Table[i, {i, 1, Length[s]}], n]], s] - 
              diag[mat[n, s], 
                mat[Complement[Table[i, {i, 1, Length[s]}], n], s]]], 2]/
        2

maxdegs[s_List] := 
      Flatten[Position[(Plus @@ #) & /@ s, 
          Max[(Plus @@ #) & /@ s]]]

elts[a_List, b_List] := Extract[a, {#} & /@ b]

\(\(\(
\)\(\(Clear[naslednik]\)\(
\)
    \(naslednik[s_List, 
        x_List] := 
		Module[{r, pat, g = s, cur = x, k, 
          l}, 
		If[
          x == Reverse[Take[s, \(-Length[x]\)]], 
          r = Take[s, 
              Length[x] + 
                1], 
		pat = {}; 
		l = Length[g]; 
		If[Position[cur, l] == {}, 
            k = 0, 
					k = \(Position[cur, 
                  l]\)[\([1, 1]\)]]; 
		If[k > 0, 
            While[k <= Length[x] && 
                g[\([l]\)] == 
                  cur[\([k]\)]\ , 
				pat = 
                Union[pat, {g[\([l]\)]}]; 
				l = 
                l - 1; 
				k = 
                k + 1], ]; 
		k = 
            k - 1; 
		If[
            k == Length[x], 
			pom = 
              Complement[g, 
                Take[x, 
                  k - Length[pat] - 
                    1]]; 
			i = \(Position[pom, 
                    x[\([k - Length[pat]]\)]]\)[\([1, 1]\)] + 
                1; 
			r = 
              Join[Take[x, 
                  k - Length[pat] - 
                    1], {pom[\([i]\)]}, 
					Take[
                  Complement[pom, {pom[\([i]\)]}], 
                  Length[pat]]], 
			pom = 
              Complement[g, 
                Drop[x, \(-1\)]]; 
			i = \(Position[
                    pom, Last[x]]\)[\([1, 1]\)] + 
                1; 
			If[
              i > Length[
                  pom], 
				pat = {Last[
                    pom]}; 
				t = 
                Length[x] - 1; 
				While[
                Last[pat] < x[\([t]\)], pat = Union[pat, {x[\([t]\)]}]; 
                t = t - 1]; 
				pom = 
                Union[pom, 
                  pat, {x[\([t]\)]}]; 
				i = \(Position[pom, x[\([t]\)]]\)[\([1, 1]\)] + 1; 
				r = 
                Join[Take[x, 
                    t - 1], {pom[\([i]\)]}, 
				Take[
                    Complement[pom, {pom[\([i]\)]}], 
                    Length[x] - t]], 
				r = 
                Join[Drop[
                    x, \(-1\)], {pom[\([i]\)]}]]
		]
		]; 
		\ \ \ Return[
          r]
	]\)\(
\)
    \( (*\ s\ is\ a\ set\ of\ integers, \ 
      x\ is\ an\ ordered\ subset\ of\ s, \ 
      output\ is\ the\ next\ ordered\ subset\ of\ s, \ after\ x, \ 
      in\ the\ lexicographic\ ordering*) \)\)\)\)

\(closest[s_List, 
        nodes_List] := 
Module[{g, comp, n = nodes, 
          ni = Table[j, {j, 1, Length[s]}], m, r, neig, \ st, cur, curmin, 
          nex, t, pat, pati, i}, 
comp = 
          components[s]; 
If[
          Length[comp] > 1, 
cur = 
            Transpose[\(closest[mat[#, s], elts[n, #]] &\) /@ 
                comp]; 
r = {Flatten[cur[\([1]\)], 1], 
              Plus @@ cur[\([2]\)]},  (*If\ not - 
              from\ here*) 
g = s; 
If[
            isComplete[g], 
            r = {{n}, 0}, 
curmin = {{n}, 
                complete[g]}; 
m = 
              maxdegs[s]; 
For[i = 1, 
              i <= Length[m], \(i++\), 
neig = 
                Union[Complement[
                    neighbors[g, 
                      m[\([i]\)]], {m[\([i]\)]}]]; 
st = {First[neig]}; 
pat = 
                Extract[
                  n, \({#} &\) /@ 
                    Append[st, m[\([i]\)]]]; 
pati = 
                Extract[
                  ni, \({#} &\) /@ 
                    Append[st, m[\([i]\)]]]; 
nex = 
                Complement[ni, pati]; 
cur = 
                closest[mat[nex, s], 
                  Extract[n, \({#} &\) /@ nex]]; 
t = 
                complete[mat[pati, s]] + uncomplete[s, pati] + 
                  cur[\([2]\)]; 
If[t < curmin[\([2]\)], 
                curmin = {Append[cur[\([1]\)], pat], 
                    t}, ]; 
While[
                Not[st == Reverse[neig]], 
st = 
                  naslednik[neig, st]; \ 
pat = 
                  Extract[
                    n, \({#} &\) /@ 
                      Append[st, m[\([i]\)]]]; 
pati = 
                  Extract[
                    ni, \({#} &\) /@ 
                      Append[st, m[\([i]\)]]]; 
nex = 
                  Complement[ni, pati]; 
If[
                  complete[mat[pati, s]] + 
                      uncomplete[s, pati] ≥ 
                    curmin[\([2]\)], , 
cur = 
                    closest[mat[nex, s], 
                      Extract[n, \({#} &\) /@ nex]]; 
t = 
                    complete[mat[pati, s]] + uncomplete[s, pati] + 
                      cur[\([2]\)]; 
If[
                    t < curmin[\([2]\)], 
                    curmin = {Append[cur[\([1]\)], pat], 
                        
                        t}, ];
] (**) 
] (*While*) 
];  (*For*) 
r = 
              curmin;
] (*If*) 
];  (*If*) 
Return[r]]	\)

closestN[s_List, 
        nodes_List] := 
Module[{g, comp, n = nodes, 
          ni = Table[j, {j, 1, Length[s]}], m, r, neig,  st, cur, curmin, 
          nex, t, pat, pati, i, max = Length[s], 
          l}, 
comp = components[s]; 
If[
          Length[comp] > 1, 
cur = 
            Transpose[closestN[mat[#, s], elts[n, #]] & /@ 
                comp]; 
r = {Flatten[cur[[1]], 1], 
              Plus @@ cur[[2]]},  (*If not - 
              from here*) 
g = s; 
If[
            isComplete[g], 
            r = {{n}, 0}, 
curmin = {{n}, 
                complete[g]}; 
m = 
              maxdegs[s]; 
i = 
              Random[Integer, {1, Length[m]}]; 
neig = 
              Union[Complement[
                  neighbors[g, 
                    m[[i]]], {m[[i]]}]]; 
st = {First[neig]}; 
pat = 
              Extract[n, {#} & /@ 
                  Append[st, m[[i]]]]; 
pati = 
              Extract[ni, {#} & /@ 
                  Append[st, m[[i]]]]; 
nex = 
              Complement[ni, pati]; 
cur = 
              closestN[mat[nex, s], 
                Extract[n, {#} & /@ nex]]; 
t = 
              complete[mat[pati, s]] + uncomplete[s, pati] + 
                cur[[2]]; 
If[t < curmin[[2]], 
              curmin = {Append[cur[[1]], pat], 
                  t}, ]; 
For[l = 1, 
              l ≤ max, l++, 
st = 
                Union[Table[
                    Random[Integer, {1, Length[neig]}], {j, 1, 
                      Length[neig]}]];  
pat = 
                Extract[
                  n, {#} & /@ 
                    Append[st, m[[i]]]]; 
pati = 
                Extract[
                  ni, {#} & /@ 
                    Append[st, m[[i]]]]; 
nex = 
                Complement[ni, pati]; 
If[
                complete[mat[pati, s]] + uncomplete[s, pati] ≥ 
                  curmin[[2]], , 
cur = 
                  closestN[mat[nex, s], 
                    Extract[n, {#} & /@ nex]]; 
t = 
                  complete[mat[pati, s]] + uncomplete[s, pati] + 
                    cur[[2]]; 
If[t < curmin[[2]], 
                  curmin = {Append[cur[[1]], pat], 
                      t}, ]];  (**) 
];  (*For l*) 
r = 
              curmin;
] (*If*) 
];  (*If*) 
Return[r]]

closestNb[s_List, 
        nodes_List] := 
Module[{g, comp, n = nodes, 
          ni = Table[j, {j, 1, Length[s]}], m, r, neig,  st, cur, curmin, 
          str, nex, t, pat, pati, i, max = Length[s], maxd = 3, l, 
          ll}, 
comp = components[s]; 
If[
          Length[comp] > 1, 
cur = 
            Transpose[closestN[mat[#, s], elts[n, #]] & /@ 
                comp]; 
r = {Flatten[cur[[1]], 1], 
              Plus @@ cur[[2]]},  (*If not - 
              from here*) 
g = s; 
If[
            isComplete[g], 
            r = {{n}, 0}, 
curmin = {{n}, 
                complete[g]}; 
m = 
              maxdegs[s]; 
i = 
              Random[Integer, {1, Length[m]}]; 
neig = 
              Union[Complement[
                  neighbors[g, 
                    m[[i]]], {m[[i]]}]]; 
st = {First[neig]}; 
pat = 
              Extract[n, {#} & /@ 
                  Append[st, m[[i]]]]; 
pati = 
              Extract[ni, {#} & /@ 
                  Append[st, m[[i]]]]; 
nex = 
              Complement[ni, pati]; 
cur = 
              closestN[mat[nex, s], 
                Extract[n, {#} & /@ nex]]; 
t = 
              complete[mat[pati, s]] + uncomplete[s, pati] + 
                cur[[2]]; 
If[t < curmin[[2]], 
              curmin = {Append[cur[[1]], pat], 
                  t}, ]; 
For[l = 1, 
              l ≤ max, l++, 
st = 
                Union[Table[
                    Random[Integer, {1, Length[neig]}], {j, 1, 
                      Length[neig]}]];  
 (*If[
                    Random[] ≥   .1, 
For[
                      ll = 1, 
                      ll ≤ 
                        maxd,  ll++,   
] (*For*) , 
]; (*If*) *) 
While[Random[] >  .8, 
                str = neighbors[g, 
                    Take[st, Random[Integer,  {1, Length[st]}]]];  
                st = Union[st, str]]; 
pat = 
                Extract[
                  n, {#} & /@ 
                    Append[st, m[[i]]]]; 
pati = 
                Extract[
                  ni, {#} & /@ 
                    Append[st, m[[i]]]]; 
nex = 
                Complement[ni, pati]; 
If[
                complete[mat[pati, s]] + uncomplete[s, pati] ≥ 
                  curmin[[2]], , 
cur = 
                  closestN[mat[nex, s], 
                    Extract[n, {#} & /@ nex]]; 
t = 
                  complete[mat[pati, s]] + uncomplete[s, pati] + 
                    cur[[2]]; 
If[t < curmin[[2]], 
                  curmin = {Append[cur[[1]], pat], 
                      t}, ]];  (**) 
];  (*For l*) 
r = 
              curmin;
] (*If*) 
];  (*If*) 
Return[r]]

findNclosest[s_List, max_, st_] := 
      Module[{cur, t, curmin, 
          n = Table[j, {j, 1, Length[s]}]}, 
t = 
          complete[s]; 
curmin = {{n}, 
            t}; 
For[f = st, 
          f ≤ max + st, f++, 
SeedRandom[
            f]; 
cur = 
            closestN[s, n]; 
If[
            cur[[2]] < curmin[[2]], 
            curmin = cur, ]
]; 
Return[
          curmin]
]

findNclosestb[s_List, max_, st_] := 
      Module[{cur, t, curmin, 
          n = Table[j, {j, 1, Length[s]}]}, 
t = 
          complete[s]; 
curmin = {{n}, 
            t}; 
For[f = st, 
          f ≤ max + st, f++, 
SeedRandom[
            f]; 
cur = 
            closestNb[s, n]; 
If[
            cur[[2]] < curmin[[2]], 
            curmin = cur, ]
]; 
Return[
          curmin]
]
 (*s is a 0 - 1 matrix of connections,  
      max is the number of iterations,  
      st is the seedrandom number,  
      output is the closestN partition,  
      after running the closestN[] max times with different seeds*)

Clear[consolidation]
consolidation[net_List, par_List] := 
      Table[Plus @@ Flatten[net[[par[[i]], par[[j]]]]], {i, 1, 
          Length[par]}, {j, 1, Length[par]}]

Clear[blockdiag]
blockdiag[{a_List, b_List}] := 
      BlockMatrix[{{a, 
            Table[0, {i, Length[a]}, {j, Length[b]}]}, {Table[
              0, {i, Length[b]}, {j, Length[a]}], 
            b}}]
blockdiag[a_List] := 
      blockdiag[{a[[1]], blockdiag[Rest[a]]}]
blockdiag[{a_List}] := a

Clear[Novo]
Novo[par_List, matr_List] := 
      blockdiag[mat[#, matr] & /@ par]

pies[par_List, net_List, siz_List] := 
        N[(Plus @@ Flatten[Novo[par, net]])/(Plus @@ 
                Flatten[Novo[par, siz]])];

      
poes[par_List, net_List, siz_List] := 
        N[((Plus @@ Flatten[Novo[{Join @@ par}, net]]) - (Plus @@ 
                    Flatten[Novo[par, net]]))/((Plus @@ 
                    Flatten[Novo[{Join @@ par}, siz]]) - (Plus @@ 
                    Flatten[
                      Novo[par, siz]]))];

bestParOnce[s_List, g_List, pin_, 
          pout_] := 
Module[{m, i, j, l = Length[s], 
            r}, 
m = 
            Table[If[i == j, 1, 0], {i, 1, l}, {j, 1, 
                l}]; 
For[i = 1, 
            i < l, i++, 
For[j = i + 1, 
              j ≤ 
                l, j++, 
If[(g[[i, j]] + 
                          g[[j, i]]) 
                      Log[pin] + (s[[i, j]] + s[[j, i]] - 
                          g[[i, j]] - g[[j, i]]) 
                      Log[1 - 
                          pin] ≥ 
(g[[i, 
                              j]] + g[[j, i]]) 
                      Log[pout] + (s[[i, j]] + s[[j, i]] - 
                          g[[i, j]] - g[[j, i]]) 
                      Log[1 - pout], 
m[[i, j]] = 1;  
                m[[j, i]] = 
                  1, ]
]
]; 
r = findNclosest[m, 100, 
                0][[1]]; 
Return[
            r];
];
 (*s is a matrix of maximal intensities,  g is a matrix of connection intensities,  
      pin and pout are in -  and out -  
        citation probabilities .  Output is the approx .  
            ml partition,  
      obtained by applying findNclosest.*) 

bestParIt[s_List, g_List, pin_, 
        pout_] := 
Module[{pie, poe, l = Length[s], cons, 
          cong, ro, rn, rc, conl, lik}, 
rn = 
          bestParOnce[s, g, pin, pout]; 
conl = 
          l; 
rc = rn; 
ro = 
          rn; 
While[
          Length[rn] < conl, 
cons = 
            consolidation[s, rn]; 
cong = 
            consolidation[g, rn]; 
rn = 
            bestParOnce[cons, cong, pin, pout]; 
rc = 
            Map[Flatten, Extract[rc, #] & /@ Flatten[#, 1] & /@ 
                Map[{#} &, rn, 2]]; 
If[
            lrte[ro, rc, g, s, pin, pout] > 0, rc = ro; conl = Length[rn], 
            ro = rc; 
            conl = Length[
                cons]];
]; 
pie = 
          pies[rc, g, s]; 
poe = 
          poes[rc, g, s]; 
lik = 
          loglik[rc, g, s, pin, pout]; 
Return[{rc, lik, 
            pin,  pout, pie, 
            poe}]
]

kin[pin_, po_] = 
        Log[((pin (1 - po))/(po (1 - pin)))];
 (*a constant used in ml - method*) 
ko[pin_, po_] = 
        Log[((1 - pin)/(1 - po))];
  
Clear[lrte, par1, par2]
lrte[par1_List, par2_List, net_List, siz_List, pin_,  
        po_] := 
((Plus @@ 
                  Flatten[Novo[par1, net]]) - (Plus @@ 
                  Flatten[Novo[par2, net]])) kin[pin, po] + 
        ko[pin, po] ((Plus @@ Flatten[Novo[par1, siz]]) - (Plus @@ 
                  Flatten[Novo[par2, siz]]))
  

loglik[par_List, net_List,  siz_List,  pin_, 
        po_] := (Plus @@ Flatten[Novo[par, net]]) 
          Log[(pin/(1 - pin))] + (Plus @@ Flatten[Novo[par, siz]]) 
          Log[1 - pin] + (Plus @@ Flatten[net] - 
              Plus @@ Flatten[Novo[par, net]]) 
          Log[(po/(1 - po))] + (Plus @@ Flatten[siz] - 
              Plus @@ Flatten[Novo[par, siz]]) 
          Log[1 - po]

mlPartitionGrid[net_List,  siz_List,  pin_List,  pout_List,  
        deltain_, 
        deltaout_] := 
Module[{a, curm, i, 
          j}, 
a = 
          Table[bestParIt[siz, net, pi, po], {pi, pin[[1]], pin[[2]], 
              deltain}, {po, pout[[1]], pout[[2]], 
              deltaout}]; 
For[curm = a[[1, 1]]; i = 1; 
          j = 1, i ≤ Length[a]  &&  
            j ≤ 
              Length[a[[1]]], i++; j++, 
If[
            a[[i, j, 2]] > curm[[2]], 
            curm = a[[i, 
                  j]], ]
]; 
Return[{curm, 
            a}];
]

     (*net is a matrix of intensities,  
      siz is a matrix of maximal intensities,  
      pin is the interval for in - rates,  
      pout is the interval for out rates,  
      delta is the grid size,  output is the partition,  
      and all the estimates of the maximal partition on the grid *)

pseudoML[s_List, g_List, pin_, pout_] := 
      Module[{i, j, l = Length[s], r = 0}, 
For[i = 1, 
          i ≤  l, i++, 
For[j = 1, 
            j ≤ 
              l, j++, 
If[(g[[i, j]] + 
                          g[[j, i]]) 
                      Log[pin] + (s[[i, j]] + s[[j, i]] - 
                          g[[i, j]] - g[[j, i]]) 
                      Log[1 - 
                          pin] ≥ 
(g[[i, 
                              j]] + g[[j, i]]) 
                      Log[pout] + (s[[i, j]] + s[[j, i]] - 
                          g[[i, j]] - g[[j, i]]) 
                      Log[1 - pout]  ||  i == j, 
              r += (g[[i, j]]) 
                    Log[pin] + (s[[i, j]] + s[[j, i]] - 
                        g[[i, j]] - g[[j, i]]) 
                    Log[1 - pin], 
r += (g[[i, j]] + 
                        g[[j, i]]) 
                    Log[pout] + (s[[i, j]] + s[[j, i]] - 
                        g[[i, j]] - g[[j, i]]) 
                    Log[1 - 
                        pout]
] (*If*) 
]
]; 
Return[
          r];]

pseudoMLNew[s_List, g_List, pin_, pout_] := 
      Module[{i, j, l = Length[s], r = 0}, 
For[i = 1, 
          i ≤  l, i++, 
For[j = 1, 
            j ≤ 
              l, j++, 
If[(g[[i, j]] + 
                          g[[j, i]]) 
                      Log[pin] + (s[[i, j]] + s[[j, i]] - 
                          g[[i, j]] - g[[j, i]]) 
                      Log[1 - 
                          pin] ≥ 
(g[[i, 
                              j]] + g[[j, i]]) 
                      Log[pout] + (s[[i, j]] + s[[j, i]] - 
                          g[[i, j]] - g[[j, i]]) 
                      Log[1 - pout]  ||  i == j, 
              r += (g[[i, j]]) 
                    Log[pin] + (s[[i, j]] - g[[i, j]]) 
                    Log[1 - pin], 
r += (g[[i, j]]) 
                    Log[pout] + (s[[i, j]] - g[[i, j]]) 
                    Log[1 - 
                        pout]
] (*If*) 
]
]; 
Return[
          r];]

Clear[blockdiag]
blockdiag[{a_List, b_List}] := 
      BlockMatrix[{{a, 
            Table[0, {i, Length[a]}, {j, Length[b]}]}, {Table[
              0, {i, Length[b]}, {j, Length[a]}], 
            b}}]
blockdiag[a_List] := 
      blockdiag[{a[[1]], blockdiag[Rest[a]]}]
blockdiag[{a_List}] := a

Clear[Novo]
Novo[par_List, matr_List] := 
      blockdiag[mat[#, matr] & /@ par]

pseudoMLMatrix[s_List, g_List, pin_, pout_] := 
      Module[{i, j, l = Length[s], 
          r = Table[
              0, {n, 1, Length[s]}, {m, 1, 
                Length[s]}]}, 
For[i = 1, 
          i ≤  l, i++, 
For[j = 1, 
            j ≤ 
              l, j++, 
If[(g[[i, j]] + 
                          g[[j, i]]) 
                      Log[pin] + (s[[i, j]] + s[[j, i]] - 
                          g[[i, j]] - g[[j, i]]) 
                      Log[1 - 
                          pin] ≥ 
(g[[i, 
                              j]] + g[[j, i]]) 
                      Log[pout] + (s[[i, j]] + s[[j, i]] - 
                          g[[i, j]] - g[[j, i]]) 
                      Log[1 - pout]  ||  i == j, 
              r[[i, j]] = 
                1, 
] (*If*) 
]
]; 
Return[
          r];]

koko = Import["all.csv", "CSV"];

sizes = 
        Import["Workbook3.txt", "Table"];
 (*This command imports the data on maximal possible numbers of citations between each pair of journals;  
      the data are stored in the matrix sizes .  Here,  
      the data are in  . 
          txt format .  IMPORTANT :  when using this matrix,  
      it has to be multiplied by 9,  
      as the data in the matrix are based on sizes of journals for 1  year,  
      while the data on intesities are for 3 years .  
          This is just a scaling factor which does not affect the partition but just the estimated probabilities and likelihoods .  
          Again,  
      remember to correctly set the directory with SetDirectory[]*)

pb = {{2}, {6}, {13}, {23}, {33}, {35}, {42}, {22, 
            34}, {32}, {15}, {16}, {30}, {36}, {41}, {14}, {8, 12, 18}, {24, 
            26}, {28, 39, 9, 27}, {37}, {21, 29}, {11, 19, 40, 1, 31}, {4, 3, 
            20}, {25, 5, 38}, {17, 10, 7}};

\(\(journalNames = {"American Economic Review (AER)", "American Journal of Agricultural Economics (AJAE)", "Brooking Papers of Economic Activity (BPEA)", "Cambridge Journal of Economics (CJE)", "Canadian Journal of Economics (CJE)", "Economic Geography (EG)", "Econometrica (E)", "Economic History Review (EHR)", "Economic Inquiry (EI)", "Economic Journal (EJ)", "European Economic Review (EER)", "Exploration of Economic History (EEH)", "Health Economics (HE)", "International Economic Review (IER)", "Journal of Comparitive Economics (JCE)", "Journal of Development Economics (JDE)", "Journal of Econometrics (JE)", "Journal of Economic History (JEH)", "Journal of Economic Literature (JEL)", "Journal of Economic Perspectives (JEP)", "Journal of Economic Theory (JET)", "Journal of Environmental Economics and Management (JEEM)", "Journal of Health Economics (JHE)", "Journal of Human Resources (JHR)", "Journal of International Economics (JIE)", "Journal of Labor Economics (JLE)", "Journal of Law and Economics (JLE2)", "Journal of Law, Economics and Organization (JLEO)", "Journal of Mathematical Economics (JME)", "Journal of Monetary Economics (JME)", "Journal of Political Economy (JPE)", "Journal of Public Economics (JPE2)", "Journal of Urban Economics (JUE)", "Land Economics (LAE)", "National Tax Journal (NTJ)", "Oxford Bulletin of Economics and Statistics (OBES)", "Oxford Economic Papers (OEP)", "Quarterly Journal of Economics (QJE)", "Rand (Bell) Journal of Economics (RJE)", "Review of Economic Studies (RES)", "Review of Economics and Statistics (RES2)", "World Development (WD)"};\)\)

Export["partitionPB.tex", 
      TableForm[journalNames[[#]] & /@ 
          Reverse[Sort[Sort /@ pb]]], "TeX"]

heuristic = {{1, 3, 5, 11, 19, 20, 30, 31, 38}, {2}, {4, 10, 
            37}, {6}, {7, 14, 17, 21, 29, 39, 40}, {8, 12, 18}, {9, 27, 
            28}, {13}, {15}, {16}, {22, 34}, {23}, {24, 
            26}, {25}, {32}, {33}, {35}, {36}, {41}, {42}};
 (*This partition is included for a benchmark comparison,  
      as it is not obtained by the methods described here,  
      instead it was obtained by a heuristic and extensive search (searching all the nodes of the decision tree for the subtrees with less than 15 nodes,  and performing some consolidations) .  
          This partition attains higher likelihoods,  however,  
      the heuristic algorithm is not nicely implementable.*)