Approximating a distribution with an exponential family using a custom loss function

Exponential Family is a family of (normalized) functions which is log-linear. In other words, if we take logarithm of each point, we will get a vector space. Various statistical distribution families like Normal, Gamma, Beta are examples of exponential families. 
Here I look exponential families of the form 1/(1+Exp[polynomial of degree n]), and give a procedure to fit such exponential family to given distribution using a custom loss function.

The function

g[x_]=0.5 1/(1+Exp[-100 x])+0.25+0.2Sin[20 x];
points=Table[i,{i,-1,1,0.01}];
Show[Plot[{0.5,g[x]},{x,-1,1}],ListPlot[Thread[{points,g[#]&/@points}]]]

Losses

Log - loss truncated to avoid division by 0 errors

loglosst=-#1Log[((Min[#2,0.999])/(1-Min[#2,0.999]))]-Log[1-Min[#2,.999]]&;

Plot[loglosst[0.6,x],{x,0,1},PlotLabel->"Log loss"]

sqloss=(((#1-#2))^2)&;
Plot[sqloss[0.6,x],{x,0,1},PlotLabel->"Squared error"]

Fitter

myfit[loss_,dims_,points_,g_]:=Module[{params,polys,f,timing,sol,pt},
params=Thread[(Subscript[p, #])&/@Table[i,{i,0,dims-1}]];
polys=LegendreP[#,x]*(2#+1)&/@Table[i,{i,0,dims-1}]; (* orthonormal polynomials*)
f=1/(1+Exp[-Plus@@Thread[polys*params]]);
realloss=Plus@@(loss[g[#],f/.{x->#}]&/@points) ;
p0=Table[0,{i,0,dims-1}];
p0[[2]]=10; (* set starting point to the optimal classifier *)
timing=Timing[pt=FindMinimum[realloss,Thread[{params,p0}]]];
sol=pt[[2]];
{timing[[1]],f/.sol,params/.sol}
];

Demo

logLossFit=myfit[loglosst,10,points,g][[2]];
sqLossFit=myfit[sqloss,10,points,g][[2]];

Plot[{logLossFit,sqLossFit,.5}//Evaluate,{x,-1,1},PlotLegend->{"logloss","squared error","threshold"},PlotLabel->"10 parameter exponential family, squared loss minimizer overshoots"]

Rank-1 decomposition to improve GraphPlot

Consider applying Thorp shuffle to 4 cards.This can be viewed as a random walk on a graph with 4! nodes which we can visualize using GraphPlot. The graph is highly regular, however, GraphPlot doesn' t detect this regularity, and the result looks quite jumbled. We can fix this by using rank-1 decomposition (SVD) to suggest a natural grouping of nodes. Group nodes together if they are source nodes in the same rank-1 component of the adjacency matrix, and a nice octahedral structure emerges.

swap[l_,p_Integer]:=Module[{ll=l},ll[[{2 p-1,2 p}]]=ll[[{2 p,2 p-1}]];ll];
swap[l_,p_List]:=Fold[swap[#1,#2]&,l,Position[p,True]//Flatten]
thorp[l_]:=Module[{k,ll},k=Length[l]/2;
ll=Riffle[l[[;;k]],l[[k+1;;]]];
FromDigits[swap[ll,#]]&/@Tuples[{False,True},k]]
thorpEdges=Thread[FromDigits[#]->thorp[#]]&/@Permutations[Range[4]]//Flatten;
graph1=GraphPlot[thorpEdges,DirectedEdges->True,PlotLabel->"original"];
edges2mat[edges_,size_]:=Module[{perms=Permutations[Range[size]]},Array[If[MemberQ[edges,FromDigits[perms[[#1]]]->FromDigits[perms[[#2]]]],1,0]&,{Length[perms],Length[perms]}]];
thorpMat=edges2mat[thorpEdges,4];
{u,d,v}={#1,Diagonal[#2],#3}&@@(SingularValueDecomposition[thorpMat]);
rank1mats=Table[d[[i]]*Outer[Times,u[[All,i]],v[[All,i]]],{i,1,6}];
outNodes[mat_]:=MemberQ[#,1]&/@mat//Position[#,True]&//Flatten;
partitions=outNodes[#]&/@rank1mats;
mergeGraph[adjMat_,partitions_]:=Module[{l=Length[partitions]},Table[Unitize[Total[adjMat[[partitions[[i]],partitions[[j]]]],2]],{i,1,l},{j,1,l}]];
graph2=mergeGraph[thorpMat,partitions]//GraphPlot[#,DirectedEdges->True]&;
p2=Permutations[partitions][[66]];
rules[i_]:=Thread[p2[[i]]->Table[{Cos[i 2 Pi/6],Sin[i 2 Pi/6]}+2/20 {Cos[j 2 Pi/4],Sin[j 2 Pi/4]},{j,1,4}]];
allrules=rules[#]&/@Range[6]//Flatten;
graph3=GraphPlot[thorpMat,DirectedEdges->True,VertexCoordinateRules->allrules,PlotLabel->"Improved vertex positioning"];
GraphicsColumn[{graph1,graph3}]

Closed form for Logistic Regression

The code below constructs and solves maximum likelihood equations for saturated Logistic regression model over k binary variables.

solvelog[k_]:=Module[{design,params,eq,eqs,odd,odds,ts,subst,gs,sol1,tsol,sol2},
design=KroneckerProduct@@Table[{{1,0},{1,1}},{i,1,k}]; (* form design matrix corresponding to conjuction features *)
params=(Subscript[θ, #])&/@Range[(2^k)];
(* equation for i'th feature *)
eq[i_]:=Sum[Transpose[design][[i]][[x]](((Exp[design[[x]].params])/(Exp[design[[x]].params]+1))-g[x]),{x,1,(2^k)}]==0;
eqs=eq[#]&/@Range[1,(2^k)];
(* predicted odds for given x *)
odd[x_]:=((Exp[design[[x]].params])/(Exp[design[[x]].params]+1));
odds=odd[#]&/@Range[(2^k)];
ts=(Subscript[t, #])&/@Range[(2^k)];
subst=Thread[odds->ts];
gs=g[#]&/@Range[(2^k)] ;(* these are empirical probabilities of y=1 for each x *)
sol1=Solve[eqs/.subst,rightvars] ;(* solution for odds in terms of empirical odds, trivial for saturated models *)
tsol=MapThread[Solve,{Thread[odds==ts],params}]//Flatten ;
sol2=tsol[[All,2]]//.tsol;
sol2/.sol1
];
solvelog[2]//Simplify

Visualizing Fisher' s Linear Discriminant

Fisher' s linear discriminant attempts to find an optimal separating hyperplane by maximizing the ratio of between - class scatter to within - class - scatter. The dynamic section shows separating boundary + class outlines

(* Draws points of a particular class *)
drawpoints[p_,class_]:=Module[{pp,pt},
pp=Disk[#,.03]&/@p;
Graphics[{If[class==1,Red,Blue],Thick}~Join~pp]
];
(* Converts contour of a quadratic form to an ellipse (ContourPlot too slow) *)
quad2ellipse[mat_,center_,threshold_, class_]:=Module[{evec,evals,c},
evec=Last@Eigenvectors[mat];
evals=Eigenvalues[mat];
c=Circle[center,{Sqrt[threshold/First@evals],Sqrt[threshold/Last@evals]}];
Graphics[{If[class==1,Red,Blue],Rotate[c,-ArcTan@(Divide@@evec)]}]
]
(* Draws Fisher discriminant for a given set of points *)
fishers[points_]:=Module[{means,d,scatter,Sb,Sw,normal,midpoint,w0,classContours,classBoundary,gpoints},
means=((Total[#])/n)&/@points;
d=means[[1]]-means[[2]];
(* Compute the scatter matrices *)
Sb=Outer[Times,d,d];
scatter[xs_,m_]:=Total[Table[Outer[Times,xs[[i]]-m,xs[[i]]-m],{i,1,n}]];
Sw=Total[scatter[points[[#]],means[[#]]]&/@Range[2]];
(* Solve the generalized eigenvalue problem *)
normal=Eigenvectors[Inverse[Sw].Sb][[1]];
(* Offset is mid-point between class means, projected onto the discriminant vector *)
midpoint=(means[[2]]+means[[1]])/2;
w0=midpoint.normal;
classContours:=quad2ellipse[Inverse[Sw],means[[#]],(midpoint-means[[1]]).Inverse[Sw].(midpoint-means[[1]]), #]&/@Range[2];
gpoints:=drawpoints[points[[#]],#]&/@Range[2];
classBoundary:=ContourPlot[normal.{x1,x2}==w0,{x1,-3,3},{x2,-3,3},ContourStyle->{Thick}];
Show[gpoints~Join~classContours~Append~classBoundary,PlotRange->{{-1.1,1.1},{-1.1,1.1}}]
];
n=5;
(* initial random points *)
ipoints=RandomReal[{-1,0},{2*n,2}];
ipoints[[n+1;;2n]]=ipoints[[n+1;;2n]]+1;
Manipulate[fishers[{pppoints[[1;;n]],pppoints[[n+1;;2n]]}],{{pppoints,ipoints},Locator},SaveDefinitions->True]

Computing Effective Resistance between 2 points in a network

I programmed this in order to access a Russian physics forum which requires you to find resistance between end - nodes of a Wheatstone bridge in order to gain access. In the process of researching this issue I found 7 other ways to compute this resistance which I' d like to implement sometime

mat={{0,1,3,0},{1,0,6,2},{3,6,0,3},{0,2,3,0}};
edges={1->2,1->3,3->2,2->4,3->4};
Off[Part::"pspec"];
edgeMap=edges/.{(a_->b_)->((a->b)->{a->b,mat[[a,b]]})};
edgeMapWeight=edges/.{(a_->b_)->((a->b)->mat[[a,b]])};
GraphPlot[edges/.edgeMap,VertexLabeling->True,PlotLabel->"Resistor Network"]

lapl[m_]:=DiagonalMatrix[Array[1&,Length[m]].m]-m;
resmat[m_]:=Module[{n=Length[m]},
g=PseudoInverse[lapl[m/.{x_?Positive->1/x}]];
Table[g[[i,i]]+g[[j,j]]-g[[i,j]]-g[[j,i]],{i,1,n},{j,1,n}]
];
resmat[mat]//MatrixForm

MaxEnt distribution with a skewed sufficient statistic

If we maximize entropy subject to moment constraints we can get Normal, Exponential, Gamma and a few other distributions, depending on a particular form of constraints. If our constraints are on first, second and third moments, then there' s no closed form for the maxent distribution. Below I find the maxent solution numerically. I used NullSpace to turn this into lower-dimensional problem, LinearSolve to find a starting solution satisfying positivity constraints, and FindMaximum for the final solution

getdist[n_,var_,skew_]:=Module[{},
h[x_]=-Piecewise[{{x Log[x],x>0}},0];
SetAttributes[h,Listable];
ent[x_]:=Total[h[x]];
vals=Table[i,{i,1/n,1,1/n}];
vals=Reverse[-vals]~Join~vals;
cons={vals^0,vals^1,vals^2,vals^3}//N;
ns=Orthogonalize[NullSpace[cons]]//Transpose;
vars=Table[(Subscript[x, i]),{i,1,Dimensions[ns]//Last}];
sol=LinearSolve[cons,{1,0,var,skew}];
realsol=LinearProgramming[ConstantArray[1,Length[vars]],ns,-sol,-Infinity]//Chop;
distr=(ns.vars+sol);
poscons=And@@(#>=0&/@distr);
obj=ent[distr];
distsol=distr/.(FindMaximum[{obj,poscons},Thread[{vars,realsol}]]//Last)
]

dd=getdist[10,.4,.23];
ListPlot[dd,Filling->Axis,PlotRange->{{-1.1,1.1},{0,.3}},DataRange->{-1,1},PlotStyle->PointSize[Medium]]

Density Plot of the magnitudes in Farey Sequence (Subscript[F, 12])

Actually of the logarithm of the magnitudes. Self-similar structure emerges

Clear[prod];
farey=Block[{
A={{1,1},{0,1}},
B={{1,0},{1,1}},
n=12,
prods},
prods=Reap[
prod[{}]=IdentityMatrix[2];
prod[l_]:=(Sow[prod[l]=If[l[[1]]==1,A,B].prod[Rest[l]]]);
MatrixForm[prod[#]]&/@Tuples[{0,1},n];
][[2,1]];
Divide@@Plus@@Transpose[#]&/@prods
];
Block[{lf=Sort[Log[farey],LessEqual]},
ListPlot[Thread[{MovingAverage[lf,2],1/Differences[lf]}],PlotStyle->{PointSize[.001]}]
]

Generating a random doubly stochastic matrix

Method 1, iterative proportional fitting

This method normalizes rows/columns in alternation. It converges with probability one and produces a doubly stochastic matrix that is closest to the random starting matrix in KL - divergence

n=3;m=RandomReal[{1},{n,n}];
e=Array[1&,n];
result=NestWhile[normalize,m,Total[Abs[e.#-e]]+Total[Abs[#.e-e]]>10^-10&];
MatrixForm/@({{e}.result,result,result.Transpose[{e}]})

Method 2, MCMC (From paper by Diaconis, "Algebraic Algorithms For Sampling ...")

This method constructs a Markov chain which randomly perturbs the matrix at each step, without destroying double - stochasticity. It is initialized with a random permutation matrix.

randCorner[n_]:=Module[{corn1,corn2},
{corn1,corn2}=RandomInteger[{1,n},{2,2}];
If[Or@@MapThread[Equal,{corn1,corn2}],
randCorner[n],
{corn1,corn2}
]
];
makeStep[n_]:=Module[{s,bounds,r,i,j,k,l,mat},
Function[{m},
mat=m;
(* pick 2 corners *)
{{i,j},{k,l}}=randCorner[n];
(* find range of acceptable modification *)
bounds={Max[-m[[i,j]],-m[[k,l]],m[[i,l]]-1,m[[k,j]]-1],
Min[1-m[[i,j]],1-m[[k,l]],m[[i,l]],m[[k,j]]]};
If[bounds[[2]]>bounds[[1]],
r=RandomReal[bounds];
mat[[{i,k},{j,l}]]+={{r,-r},{-r,r}};mat,
mat
]
]
];
(* gives random permutation matrix *)
randPerm[n_]:=IdentityMatrix[n][[RandomSample[Range[n]]]];
(* run for 100 steps *)
n=3;
e=Array[1&,n];
f=makeStep[n];
result=Nest[f,randPerm[n],100];
MatrixForm/@({{e}.result,result,result.Transpose[{e}]})

Visualize Gibbs sampling of a Hard - Core Lattice Gas model

Hard - core model pretends each molecule can occupy only one of a set of fixed sites on a lattice, and neighbouring sites can not be simultaneously occupied. Code below implements Gibbs sampling to draw random configurations from this model. Visualization shows current configuration, and plot of number of sites occupied in previous iterations of MCMC run

(* takes makes 0-1 matrix symmetric *)
symmetrize[mat_]:=Unitize[mat+Transpose[mat]];
(* Convert list of edges to adjacency matrix *)
graph2mat[edges_]:=Module[{a,nodes},
(* custom flatten to allow edges be lists *)
nodes=List@@@edges//Flatten//Union;
Array[If[MemberQ[edges,nodes[[#1]]->nodes[[#2]]],1,0]&,{Length[nodes],Length[nodes]}]
]
comboPlot[g_,vals_, k_]:=Module[{window=50},
m=Max[k-window,1];
g1=ListPlot[Total[vals[[m;;k]],{2}],Ticks->{{},{{1,1},{12,12}}},PlotRange->{{m,m+window},{0,13}},DataRange->{m,k},AxesOrigin->{m,0},Axes->{False,True}];
g2=twoColorGraph[g,vals[[k]]];
Graphics[{Inset[g1,{-.3,0},Center,{1,1}],Inset[g2,{0,-1.2},Center,{2,2}]},ImageSize->{400,500},ImageMargins->None,PlotRangePadding->None]
];
twoColorGraph[g_,seen_]:=Module[{radius=.22},
fv=Function[{p,v},{EdgeForm[Thick],If[seen[[v]]==1,Pink,White],Disk[p,radius]}];
GraphPlot[g,VertexLabeling->True,VertexRenderingFunction->fv]
];
makeChain[h_,l_]:=Module[{connectLevel,connectLevels,intraconnect,interconnect},
connectLevel[start_,height_]:=x[#]->x[#+1]&/@Range[start,start+height-2];
connectLevels[start_,height_]:=x[#]->x[#+height]&/@Range[start,start+height-1];
intraconnect=connectLevel[1+h(#-1),h]&/@Range[l];
interconnect=connectLevels[1+h(#-1),h]&/@Range[l-1];
Flatten[intraconnect~Join~interconnect,{1,2}]
(*intraconnect~Join~interconnect*)
];
makeMCMC[g_,l_]:=Module[{n=Length[g]},
Function[{x},
p=RandomInteger[{1,n},{1}]//First;
Table[If[i==p,
Boole[Total@x[[Position[g[[i]],1]//Flatten]]==0]Boole[RandomReal[]<l/(1+l)],
x[[i]]],
{i,n}]
]
];
n=6;g=makeChain[n,n]//graph2mat//symmetrize;
f=makeMCMC[g,4];
x0=ConstantArray[0.,Length[g]];
vals=NestList[f,x0,500];
pics=Table[comboPlot[g,vals,Floor[k]],{k,1,Length[vals],1}];
ListAnimate[pics,20]

Bethe approximation for Ising spin-glass

The code below computes a sequence of Bethe - Peierls - like approximations for computing mean magnetization for each site in a planar Ising spin - glass model. The idea for the approximation scheme came from Alan Sokal' s "The Repulsive Lattice Gas ... " paper and is based on limited memory self-avoiding walks. For the simple graph below, memory-5 walks recover magnetization values exactly. The graph plots estimated magnetization values for 5 nodes, last number compares approximate result with exact, obtained by explicitly summering over all configurations

(* For empty path, children are all vertices *)
children[p_]:=If[Length[p]==0,Range[n],Select[Position[g[[All,Last[p]]],1]//Flatten,FreeQ[p,#]&]];
(* True if given path ends with a leaf *)
isLeaf[p_]:=Length[p]>1&&children[p[[-2;;]]]==={};
(* Whether given path ends in a dead-end (leaf or cycle) *)
deadend[p_]:=Length[children[p]]==0;

(* All self-avoiding walks up to length 2 *)
makeSAWs[k_]:=Module[{},
makeSAWsr[prefix_,j_]:=Module[{},
If[j≥1,
AppendTo[saws,prefix];
makeSAWsr[prefix~Join~{#},j-1]&/@children[prefix]
];
];
saws={};
makeSAWsr[{#},k]&/@Range[n];
saws//Sort
];

(* gets all labelings over n nodes with i'th node set to p *)
labelings[n_,i_,p_]:=Select[Tuples[{1,-1},n],#[[i]]==p&];

(* Converts a list of edges to adjacency matrix *)
graph2mat[edges_]:=Module[{a,nodes},
(* custom flatten to allow edges be lists *)
nodes=List@@@edges//Flatten//Union;
Array[If[MemberQ[edges,nodes[[#1]]->nodes[[#2]]],1,0]&,{Length[nodes],Length[nodes]}]
]

(* Converts a list of edges to adjacency matrix, allows Lists as nodes *)
graph2mat2[edges_]:=Module[{a,nodes},
(* custom flatten to allow edges be lists *)
nodes=List@@@edges//Flatten[#,{{1,2},{3}}]&//Union//Sort;
Array[If[MemberQ[edges,nodes[[#1]]->nodes[[#2]]],1,0]&,{Length[nodes],Length[nodes]}]
]
(* makes 0-1 matrix symmetric *)
symmetrize[mat_]:=Unitize[mat+Transpose[mat]];
tail[p_,l_]:=Take[p,-Min[l,Length[p]]]
tail2[p_]:=Take[p,-Min[2,Length[p]]];
(* gets possible continuations of given memory-1 walk *)
childrenL[p_]:=Select[Position[g[[All,Last[p]]],1]//Flatten,FreeQ[tail2[p],#]&];
(* Returns position of first occurence of p//Last, or -1 if no such *)
loopPos[p_]:=(
pos=Position[p[[;;-2]],p//Last]//Flatten;
If[pos=={},-1,pos//First]
)
(* returns loop type for given path (0, no loop, 1 opening, -1 closing *)
loopType[p_]:=Module[{pp},
pp=loopPos[p];
answer=
Which[
pp==-1,0,
(* closing, right order *)
p[[pp+1]]<p[[-2]],-1,
(* opening, wrong order *)
p[[pp+1]]>p[[-2]],1
];
(*Print["called for",p,answer];*)
answer
];


(* returns True if given node is a leaf *)
isLeaf[p_]:=Length[p]>1&&children[p[[-2;;]]]==={}

(* Makes a thick chain of edges *)
x=.;
makeChain[h_,l_]:=Module[{connectLevel,connectLevels,intraconnect,interconnect},
connectLevel[start_,height_]:=x[#]->x[#+1]&/@Range[start,start+height-2];
connectLevels[start_,height_]:=x[#]->x[#+height]&/@Range[start,start+height-1];
intraconnect=connectLevel[1+h(#-1),h]&/@Range[l];
interconnect=connectLevels[1+h(#-1),h]&/@Range[l-1];
Flatten[intraconnect~Join~interconnect,{1,2}]
(*intraconnect~Join~interconnect*)
];

(* almost regular tree, first level has d-1 children, the rest have d children each *)
makeTree[degree_,depth_]:=Module[{},
Ceiling[(#)/degree]->#&/@Range[2,degree^depth]
];

(* Procedures for working with dictionaries *)
keys[dict_]:=Map[#[[1]]/.{HoldPattern[dict[x_]]->x}&,DownValues[dict]]//ReleaseHold;
values[dict_]:=Map[#[[-1]]&,DownValues[dict]];

(* For empty path, children are all vertices *)
children[p_,g_]:=Module[{n=Length[g]},If[Length[p]==0,Range[n],Select[Position[g[[All,Last[p]]],1]//Flatten,FreeQ[p,#]&]]];
(* gets possible continuations of given memory-1 walk *)
childrenL[p_,g_]:=Select[Position[g[[All,Last[p]]],1]//Flatten,FreeQ[tail2[p],#]&];
(* True if given path ends with a leaf *)
isLeaf[p_,g_]:=Length[p]>1&&children[p[[-2;;]],g]==={};
(* Whether given path ends in a dead-end (leaf or cycle) *)
deadend[p_,g_]:=Length[children[p,g]]==0;

(* All self-avoiding walks up to length 2 *)
makeSAWs[k_,g_]:=Module[{makeSAWsr,saws,n=Length[g]},
makeSAWsr[prefix_,j_]:=Module[{},
If[j≥1,
AppendTo[saws,prefix];
makeSAWsr[prefix~Join~{#},j-1]&/@children[prefix,g]
];
];
saws={};
makeSAWsr[{#},k]&/@Range[n];
saws//Sort
];

(* Example graphs *)
gs={makeChain[1,4],makeTree[2,3],{1->2,2->3,3->1},{1->2,1->3,2->3,3->4,4->5,4->6,5->6},{1->2,2->3,1->3,2->4,2->5,4->5},{1->2,1->3,2->3,2->4,3->4},{1->2,2->3,3->4,4->5,1->3,2->4,3->5},makeChain[2,4]};

(* Takes matrix of interactions strengths, magnetic fields, outputs marginals *)
exactMarginals[J_,h_]:=Module[{yy,vars,H,labelings,potm,n=Length[J]},
yy[i_]:=ToExpression["y"<>ToString[i]];
vars=Array[yy[#]&,n];
H=Function@@{vars,Expand[vars.J.vars/2]+vars.h};
(* gets all labelings over n nodes with i'th node set to p *)
labelings[n_,i_,p_]:=Select[Tuples[{1,-1},n],#[[i]]==p&];

potm[i_,p_]:=Total[Exp[H@@#&/@labelings[n,i,p]]];
Table[1/2Log[potm[i,1]/potm[i,-1]],{i,n}]
];
trueM[edges_]:=Module[{g,h},
g=edges//graph2mat//symmetrize;
h=Array[Boole[#!=1]&,{Length[g]}];
exactMarginals[-g,h]//N
];
(* Makes bottom-up SAW equations for loops of given girth *)
isingSABP[J_,h_,mem_,iters_]:=Module[{makeFactor,f,bm,update,val,g=Unitize[J]},
(* constructs edge factor *)
makeFactor[p_,child_]:=Module[{path},
path=p~Join~{child};
Switch[loopType[path],
(* Regular message *)
0,m[tail[path,mem]],
(* Loop opening edge *)
1,J[[i,child]],
(* Loop closing edge *)
-1,-J[[i,child]]
]
];

(* Equation for convolution step *)
f[z_]:=Function[{x},ArcTanh[Tanh[z]Tanh[x]]];

(* Equation for message p *)
bm[p_]:=Module[{},
ch=childrenL[p,g];
contract=If[Length[p]>1,f[J[[p[[-2]],p[[-1]]]]],#&];
i=p//Last;
contract[Total[makeFactor[p,#]&/@ch]+h[[i]]]
];

vars=makeSAWs[mem,g];
(* Create fixed point equations *)
Clear[m];
eqs=Hold@@{bm/@vars};
(* Initialize all messages to 0 *)
(m[#]=0)&/@vars;
(* Single step of update iteration *)
val=ConstantArray[0,{Length[vars]}];
update:=(
val=ReleaseHold[eqs];
(m[vars[[#]]]=val[[#]])&/@Range[Length[vars]];
);
Do[update,{iters}];
val
];
approxM[edges_,mem_]:=Module[{g,h,n},
g=edges//graph2mat//symmetrize//N;
n=Length[g];
h=Array[Boole[#!=1]&,{Length[g]}];
isingSABP[-g,h,mem,50][[;;n]]//N
];
approximate=Table[approxM[gs[[7]],k],{k,2,6}];
exact=trueM[gs[[7]]];
GraphPlot[gs[[7]],VertexLabeling->True]
ArrayPlot[approximate]
(Last[approximate]-exact)//Total

SetDirectory[NotebookDirectory[]];
Export["wri-appl.html",SelectedNotebook[],"HTML"]