Crystal version of Monte Carlo simulation on when two networks or more networks fail together and for how many minutes (how robust does OOB network need to be, to have reasonable probability that both networks are not down at the same time)

netmonte.cr

module NetMonte

STANDARD_DEVIATION = 2.5
MINUTES = 365*24*60

class Simulate
  def initialize(years, *networks)
    years.times do |year|
      nets = init_networks networks
      roulette year, nets
    end
  end

  private def roulette(year, nets)
    MINUTES.times do |minute|
      nets.each do |net|
        net.run minute
      end
      if nets.all? { |net| net.down }
        puts "outage at year #{year}, minute #{minute}"
      end
    end
  end

  private def init_networks(networks)
    nets = [] of Network
    networks.each do |net|
      outage_count  = RandomGaussian.new net[:outage_count],  STANDARD_DEVIATION
      outage_length = RandomGaussian.new net[:outage_length], STANDARD_DEVIATION
      nets << Network.new(outage_count.rand, outage_length)
    end
    nets
  end

end

class Network
  getter :down
  def initialize(outage_count : Float64, outage_length)
    @down       = false
    @outage_end = 0
    @outage_probability = (1 / ((outage_count/365)/1440)).to_i64
    @outage_length = outage_length
  end
  def run(minute)
    if @down && (minute > @outage_end)
      @down = false
    end
    if rand(@outage_probability) == 42
      @down = true
      @outage_end = minute + @outage_length.rand
    end
  end
end

#stole from SE
class RandomGaussian
  def initialize(mean = 0.0, sd = 1.0, rng = ->(){ Random.rand })
    @mean, @sd, @rng = mean, sd, rng
    @compute_next_pair = false
    @g0 = 0.0 # wat
    @g1 = 0.0 # wat
  end

  def rand : Float64
    if (@compute_next_pair = !@compute_next_pair)
      # Compute a pair of random values with normal distribution.
      # See http://en.wikipedia.org/wiki/Box-Muller_transform
      theta = 2 * Math::PI * @rng.call
      scale = @sd * Math.sqrt(-2 * Math.log(1 - @rng.call))
      @g1 = @mean + (scale * Math.sin(theta))
      @g0 = @mean + (scale * Math.cos(theta))
      @g0.abs()
    else
      @g1.abs()
    end
  end
end

end

begin
  netA = { outage_count: 3.0, outage_length: 180.0 }
  netO = { outage_count: 5.0, outage_length:  60.0 }
  NetMonte::Simulate.new(10000, netA, netO)
end