Experiments benchmark.

[joaquinffernandez]

Models in increasing level of difficulty to process:
First Draft

Air Conditioners Population
Modelica model:

model airconds
  import math;
  constant Integer N = 4000;
  parameter Real CAP[N], RES[N], POT[N], THA = 32,pmax=0;
  Real th[N];
  discrete Real on[N];
  discrete Real tref[N];
  discrete Real nextSample[N];
  discrete Real noise[N];
  discrete Real nextTref[N];

  initial algorithm
    for i in 1:N loop
      th[i] := rand(4)+ 18;
      CAP[i] := rand(100)+ 550;
      RES[i] := rand(0.4)+ 1.8;
      POT[i] := rand(2)+ 13;
      pmax:=pmax+POT[i];
      nextSample[i] := 1;
      noise[i] := rand(2)-1;
      tref[i] := 20;
      nextTref[i] := 1000;
    end for;
    for i in 1:N loop
      if th[i] - tref[i]  - 0.5 > 0 then
        on[i] := 1;
      end if;
    end for;

  equation
    for i in 1:N loop  // Node 1
      der(th[i]) = (THA/RES[i]-POT[i]*on[i]-th[i]/RES[i]+noise[i]/RES[i])/CAP[i];
    end for;

  algorithm
    for i in 1:N loop // Node 2
      when th[i] - tref[i] + on[i] - 0.5 > 0 then
        on[i] := 1;
      elsewhen th[i] - tref[i] + on[i] - 0.5 < 0 then
        on[i] := 0;
      end when;
    end for;
    for i in 1:N loop // Node 3
      when time > nextTref[i] then
        if (nextTref[i] == 1000) then
          tref[i] := 20.5;
        else
          tref[i] := 20;
        end if;
        nextTref[i] := 2000;
      end when;
    end for;
    for i  in 1:N loop // Node 4
      when time > nextSample[i] then
        nextSample[i] := nextSample[i]+1;
        noise[i] := 2*abs(sin(i*time))-1;
      end when;
    end for;
end airconds;

Expected Partition:
This simple model can be splitted without any communication by distributing an equal number of elements of each node into a given partition. In fact, the distributive initial partition should be enough to partition this model. We should check that no iteration is done in this case.

Advection 1D
From this Modelica model:

model advection
  parameter Real alpha=0.5,mu=1000;
  constant Integer N = 20000;
  Real u[N];
  initial algorithm
  for i in 1:N/3 loop
    u[i]:=1;
  end for;

  equation
    der(u[1])=(-u[1]+1)*N-mu*u[1]*(u[1]-alpha)*(u[1]-1);
    for j in 2:N loop
      der(u[j])=(-u[j]+u[j-1])*N-mu*u[j]*(u[j]-alpha)*(u[j]-1);
  end for;
end advection;

We can also use this model as an example of how the SBG input is built. We have the following SBG Partitioner input JSON file:

Expected Partition:
For this simple experiment, best partition we can get is the one that given any two partitions only communicates the last element of the first partition with the first element of the second partition.

Advection 2D
Modelica Model:

model advection2D
  constant Integer N = 20;
  parameter Real ax;
  parameter Real ay;
  parameter Real r;
  parameter Real dx;
  parameter Real dy;
  Real u[N,N];

  initial algorithm

  ax := 1;
  ay := 1;
  r := 1;
  dx := 1;
  dy := 1;
  
  u[1,1] := 1; 

 equation
  
    // Origin equation
    der(u[1,1])=-ax*u[1,1]/dx - ay*u[1,1]/dy + r*(u[1,1]^2-u[1,1]^3);
    
    // Borders 
    for i in 2:N loop
      der(u[i,1])=-ax*u[i,1]/dx - ay*(u[i,1] - u[i-1,1])/dy + r*(u[i,1]^2-u[i,1]^3);
    end for;

    // Borders 
    for j in 2:N loop
      der(u[1,j])=-ax*(u[1,j] - u[1,j-1])/dx - ay*u[1,j]/dy + r*(u[1,j]^2-u[1,j]^3);
    end for;

    // Rest of the grid.
    for i in 2:N, j in 2:N loop
      der(u[i, j])=-ax*(u[i,j] - u[i,j-1])/dx - ay*(u[i,j] - u[i-1,j])/dy + r*(u[i,j]^2-u[i,j]^3);
	    end for;
end advection2D;

Expected Partition:

Air Conditioners Population with Controller
Modelica model:

model airconds_cont
  function getSection
    input Integer x;
    output Real y;
    external "C" y = getSection(x) annotation(
                                Include="#include \"sections.c\"");
  end getSection;

  import math;
  constant Integer N = 4000;
  constant Integer SECTIONS = 4;
  parameter Real CAP[N], RES[N], POT[N], THA = 32,pmax=0,Kp=1,Ki=1,tref=20;
  parameter Integer sections[N];
  parameter Integer sections_init[SECTIONS];
  parameter Integer sections_end[SECTIONS];
  Real th[N];
  Real ierr;
  discrete Real ptotals;
  discrete Real partTotal[SECTIONS];
  discrete Real sendPartTotal[SECTIONS];
  Real ptotal;
    
  discrete Real on[N];
  discrete Real dtref,pref(start=0.5);
  discrete Real nextSample(start=1);
  discrete Real partSample[SECTIONS](each start=1);
  discrete Real noise[N];
  discrete Real update[SECTIONS];
  discrete Real sampleNoise[N];
  discrete Real inputs[N];

  initial algorithm
    for i in 1:N loop
      th[i] := random(4)+ 18;
      CAP[i] := random(100)+ 550;
      RES[i] := random(0.4)+ 1.8;
      POT[i] := random(0.2)+ 13;
      pmax:=pmax+POT[i];
      noise[i] := 2*sin(i)^2;
      sampleNoise[i] := random(2);
      sections[i] := getSection(i);
    end for;
    for i in 1:SECTIONS loop
       sections_init[i] := (i-1)*N/SECTIONS;
       sections_end[i] := i*N/SECTIONS;
    end for;


    for i in 1:N loop
      if th[i] - tref   - 0.5> 0 then
        on[i] := 1;
        ptotal := ptotal + POT[i];
      end if;
    end for;

  equation
    for i in 1:N loop
      der(th[i]) = (THA/RES[i]-POT[i]*on[i]-th[i]/RES[i]+noise[i]/RES[i])/CAP[i];
    end for;
    der(ierr) = pref-ptotals/pmax;
   der(ptotal) = 0;

  algorithm
    for i in 1:SECTIONS loop
      when time > partSample[i] then
        partSample[i] := partSample[i]+1;
        update[i] := 1;
        sendPartTotal[i] := partTotal[i];
        //reinit(ptotal, ptotal + sendPartTotal[i]);        
        partTotal[i] := 0;
      end when  annotation(MMO_Event_Id="Ev_1", MMO_HH={(Ev_6[_d2],[1:SECTIONS, sections_init[i]:sections_end[i]])});
    end for;

    for i in 1:SECTIONS loop
      when update[i] > 0.5 then
        reinit(ptotal, ptotal + sendPartTotal[i]);
        update[i] := 0;
      end when;
    end for;

    when time > 1000 then
      pref := 0.4;
    end when;
    when time > 2000 then
      pref := 0.5;
    end when;

    when time > nextSample then
      nextSample := nextSample+1;
      ptotals := ptotal;
      dtref := Kp*(ptotals/pmax-pref)-Ki*ierr;
    end when;

    for i in 1:N/2 loop
      when th[i] - tref -dtref + on[i] - 0.5 > 0 then
        on[i] := 1;
        partTotal[1]:= partTotal[1] + POT[i];
      elsewhen th[i] - tref -dtref + on[i] - 0.5 < 0 then
        on[i] := 0;
        partTotal[1]:= partTotal[1] - POT[i];
    end when annotation(MMO_Event_Id="Ev_6",
                                                  MMO_LHS_DSC={partTotal[sections[i]]},
		   MMO_HH={Ev_1[sections[i]]});		
    end for;
    
    for i in N+1:N loop
      when th[i] - tref -dtref + on[i] - 0.5 > 0 then
        on[i] := 1;
        partTotal[2]:= partTotal[2] + POT[i];
      elsewhen th[i] - tref -dtref + on[i] - 0.5 < 0 then
        on[i] := 0;
        partTotal[2]:= partTotal[2] - POT[i];
    end when annotation(MMO_Event_Id="Ev_6",
                                                  MMO_LHS_DSC={partTotal[sections[i]]},
		   MMO_HH={Ev_1[sections[i]]});		
    end for;


    for i in 1:N loop
      when time > sampleNoise[i] then
        sampleNoise[i] := sampleNoise[i] + 1;
         inputs[i] := inputs[i] + 1;
        noise[i] := 2*sin(i*inputs[i])^2;
      end when;
    end for;
end airconds_cont;

Expected Partition:

Spiking Neural Network
Modelica model:

model spl_neurons_test
  function randomC
    input Integer x;
    output Real y;
    external "C" y = randomC(x) annotation(
                                Include="#include \"rneurons.c\"");
  end randomC;

  function randomN
    output Real y;
    external "C" y = randomN() annotation(
                                Include="#include \"rneurons.c\"");
  end randomN;

  function random100
    input Integer x;
    output Real y;
    external "C" y = random100(x) annotation(
                                Include="#include \"rneurons.c\"");
  end random100;

  import math;
  constant Integer N=300000, M=60000, NEURON_OFFSET=1000, NEURON_CONN=200;
  Real v[N](each start=-60),gex[N],ginh[N];
  parameter Real tau=20, vrest=-60, vthres=-50,Trefrac=5,Eex=0,Einh=-80,tauex=5,tauinh=10,dgex=0.4,dgin=1.6;
  discrete Real  active[N](each start=1),tfire[N](each start=10000);
  discrete Real  nextev[M](each start=random(1));
  parameter Real Period=10,lastTime=100;
  parameter Integer SNmap1[N];
  parameter Integer NNmap[200, N];
  parameter Real ex[N];
  discrete Real exReinit[N];
  discrete Real inhReinit[N];
  discrete Real nextGlobalTime[M];
  discrete Real inputs[M];
  initial algorithm
    for i in 1:M loop
      SNmap1[i]:=randomN();
    end for;

    for i in 1:200, j in 1:N loop
      NNmap[i, j] := randomC(i);
    end for;

    for i in 1:N loop
      ex[i]:=random100(i);
    end for;

  equation
    for i  in 1:N loop
      der(v[i])=active[i]*((vrest-v[i])+gex[i]*(Eex-v[i])+ginh[i]*(Einh-v[i]))/tau;
      der(gex[i])=-gex[i]/tauex;
      der(ginh[i])=-ginh[i]/tauinh;
      end for;

  algorithm
    for i in 1:N loop
      when v[i]*ex[i]>vthres then
        reinit(v[i],vrest);
        reinit(ginh[i],0);
        reinit(gex[i],0);
        for j in 1:NEURON_CONN loop
          exReinit[NEURON_OFFSET + j] := 1;
        end for;
        active[i]:=0;
        tfire[i]:=time;
      end when;
    end for;

    for i in 1:N loop
      when v[i]*(1-ex[i])>vthres then
        reinit(v[i],vrest);
        reinit(ginh[i],0);
        reinit(gex[i],0);
        for j in 1:NEURON_CONN loop
         inhReinit[NEURON_OFFSET+j] := 1;
        end for;
  
        active[i]:=0;
        tfire[i]:=time;
      end when;
    end for;

    for i in 1:N loop
      when exReinit[i] > 0.5 then
         reinit(gex[i],gex[i]+dgex);
         exReinit[i]:=0;
      end when;
    end for;

    for i in 1:N loop
      when inhReinit[i] > 0.5 then
         reinit(ginh[i],ginh[i]+dgin);
         inhReinit[i]:=0;
      end when;
    end for;

    for i in 1:N loop
      when time>tfire[i]+Trefrac then
        active[i]:=1;
      end when;
    end for;

    for i in 1:M loop
      when time>nextev[i] then
        exReinit[SNmap1[i]] := 1;
         inputs[i] := inputs[i] + 1;
        nextev[i]:=time+10*sin(i*inputs[i])^2;
      end when;
    end for;
end spl_neurons_test;

In this model we should also increase the number of connections and test the performance, i.e. the largest number of connections that support with the current implementation.

Expected Partition:

Spiking Neural Network with Random Connections
Modelica model:

model spl_neurons_test
  function randomC
    input Integer x;
    output Real y;
    external "C" y = randomC(x) annotation(
                                Include="#include \"rneurons.c\"");
  end randomC;

  function randomN
    output Real y;
    external "C" y = randomN() annotation(
                                Include="#include \"rneurons.c\"");
  end randomN;

  function random100
    input Integer x;
    output Real y;
    external "C" y = random100(x) annotation(
                                Include="#include \"rneurons.c\"");
  end random100;

  import math;
  constant Integer N=300000, M=60000, NEURON_CONN=200;
  Real v[N](each start=-60),gex[N],ginh[N];
  parameter Real tau=20, vrest=-60, vthres=-50,Trefrac=5,Eex=0,Einh=-80,tauex=5,tauinh=10,dgex=0.4,dgin=1.6;
  discrete Real  active[N](each start=1),tfire[N](each start=10000);
  discrete Real  nextev[M](each start=random(1));
  parameter Real Period=10,lastTime=100;
  parameter Integer SNmap1[N];
  parameter Integer NNmap[200, N];
  parameter Real ex[N];
  discrete Real exReinit[N];
  discrete Real inhReinit[N];
  discrete Real nextGlobalTime[M];
  discrete Real inputs[M];
  initial algorithm
    for i in 1:M loop
      SNmap1[i]:=randomN();
    end for;

    for i in 1:200, j in 1:N loop
      NNmap[i, j] := randomC(i);
    end for;

    for i in 1:N loop
      ex[i]:=random100(i);
    end for;

  equation
    for i  in 1:N loop
      der(v[i])=active[i]*((vrest-v[i])+gex[i]*(Eex-v[i])+ginh[i]*(Einh-v[i]))/tau;
      der(gex[i])=-gex[i]/tauex;
      der(ginh[i])=-ginh[i]/tauinh;
      end for;

  algorithm
    for i in 1:N loop
      when v[i]*ex[i]>vthres then
        reinit(v[i],vrest);
        reinit(ginh[i],0);
        reinit(gex[i],0);
        for j in 1:NEURON_CONN loop
          exReinit[NNmap[j, i]] := 1;
        end for;
        active[i]:=0;
        tfire[i]:=time;
      end when;
    end for;

    for i in 1:N loop
      when v[i]*(1-ex[i])>vthres then
        reinit(v[i],vrest);
        reinit(ginh[i],0);
        reinit(gex[i],0);
        for j in 1:NEURON_CONN loop
         inhReinit[NNmap[j, i]] := 1;
        end for;
  
        active[i]:=0;
        tfire[i]:=time;
      end when;
    end for;

    for i in 1:N loop
      when exReinit[i] > 0.5 then
         reinit(gex[i],gex[i]+dgex);
         exReinit[i]:=0;
      end when;
    end for;

    for i in 1:N loop
      when inhReinit[i] > 0.5 then
         reinit(ginh[i],ginh[i]+dgin);
         inhReinit[i]:=0;
      end when;
    end for;

    for i in 1:N loop
      when time>tfire[i]+Trefrac then
        active[i]:=1;
      end when;
    end for;

    for i in 1:M loop
      when time>nextev[i] then
        exReinit[SNmap1[i]] := 1;
         inputs[i] := inputs[i] + 1;
        nextev[i]:=time+10*sin(i*inputs[i])^2;
      end when;
    end for;
end spl_neurons_test;

This is just a stress test to check so that we can compare our results with the ones from other partitioners. We should use this experiment to explain the limitations of our proposed algorithm.

Expected Partition:
In this case there is no optimal partition.