Tag: simulation

optimization – Why does a kalman filter work better in a simulation than in real life?

I am trying to use a classical
Kalman filter
for getting indoor location. I am using geolocation API of HTML to get
latitude, longitude, position accuracy and speed. I am also using the
smartphone’s accelerometer to get acceleration. The Kalman filter
combines this sensor data. I tried making a simulation in python with
gaussian error comparable to actual measurement. The GPS measurement
has accuracy of around 20 meters.

In simulation using Kalman filter I was getting consistent accuracy of
1-2 metres for 10000 iterations. When I tried applying it in real life
the location just stays at the same place. There is almost no movement
when I move from the original location.

My question is how can I make the Kalman filter respond to the movement?

It should at least follow the direction even if absolute location may not
be accurate initially.

Here is the python snippet:

  • vx, vy are velocity,
  • ax, ay are acceleration, and
  • kx,ky are the output position of the Kalman filter.
  • numpy.eye(n) is just $I_n$, the identity matrix.
class KalmanFilterLinear:
    def __init__(self, _A, _B, _H, _x, _P, _Q, _R):
        self.A = _A  # State transition matrix.
        self.B = _B  # Control matrix.
        self.H = _H  # Observation matrix.
        self.current_state_estimate = _x  # Initial state estimate.
        self.current_prob_estimate = _P  # Initial covariance estimate.
        self.Q = _Q  # Estimated error in process.
        self.R = _R  # Estimated error in measurements.

    def GetCurrentState(self):
        return self.current_state_estimate

    def Step(self, control_vector, measurement_vector, accuracymatrix):
        # ---------------------------Prediction step-----------------------------
        predicted_state_estimate = (
            self.A * self.current_state_estimate + self.B * control_vector
        )
        predicted_prob_estimate = (
            self.A * self.current_prob_estimate
        ) * numpy.transpose(self.A) + self.Q
        # --------------------------Observation step-----------------------------
        innovation = measurement_vector - self.H * predicted_state_estimate
        innovation_covariance = (
            self.H * predicted_prob_estimate * numpy.transpose(self.H) + accuracymatrix
        )
        # -----------------------------Update step-------------------------------
        kalman_gain = (
            predicted_prob_estimate
            * numpy.transpose(self.H)
            * numpy.linalg.inv(innovation_covariance)
        )
        self.current_state_estimate = (
            predicted_state_estimate + kalman_gain * innovation
        )
        # We need the size of the matrix so we can make an identity matrix.
        size = self.current_prob_estimate.shape(0)
        # eye(n) = nxn identity matrix.
        self.current_prob_estimate = (
            numpy.eye(size) - kalman_gain * self.H
        ) * predicted_prob_estimate
        # -------Creating the Kalman filter--------------------------
        timeslice = 0.0043
        state_transition = numpy.matrix(
            ((1, timeslice, 0, 0), (0, 1, 0, 0), (0, 0, 1, timeslice), (0, 0, 0, 1))
        )
        control_matrix = numpy.matrix(
            ((1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1))
        )
        observation_matrix = numpy.eye(4)
        initial_state = numpy.matrix(((20.0), (0.0), (20.0), (0.0)))
        initial_probability = numpy.eye(4)
        process_covariance = numpy.matrix(
            (
                (timeslice ** 4 / 4, 0, timeslice ** 3 / 2, 0),
                (timeslice ** 3 / 2, 0, timeslice ** 2, 0),
                (0, timeslice ** 4 / 4, 0, timeslice ** 3 / 2),
                (0, timeslice ** 3 / 2, 0, timeslice ** 2),
            )
        )
        measurement_covariance = numpy.eye(4) * numpy.array((20, 0.2, 20, 0.2))


kx = 0.0
ky = 0.0

kf = KalmanFilterLinear(
    state_transition,
    control_matrix,
    observation_matrix,
    initial_state,
    initial_probability,
    process_covariance,
    measurement_covariance,
)

# --------updating the kalman filter----------------------

control_vector = numpy.matrix(
    (
        (0.5 * ax * timeslice * timeslice),
        (ax * timeslice),
        (0.5 * ay * timeslice * timeslice),
        (ay * timeslice),
    )
)
accuracymatrix = numpy.eye(4) * numpy.array((accuracy, accuracy, accuracy, accuracy))

kf.Step(control_vector, numpy.matrix(((x), (vx), (y), (vy))), accuracymatrix)

kx = kf.GetCurrentState()(0, 0)
ky = kf.GetCurrentState()(2, 0)

python – Efficient implementation of comparison operators for cards in a poker simulation

Just as a learning exercise, I am implementing a poker game simulator. I am purposefully not searching for ready made tutorials showing how to do this, relying mostly on the python docs – that way I hope to learn in a more memorable, albeit slower, way. I’m looking for advice on how to code the comparison between cards, so that e.g. Five beats Four.

My current implementation works, but I feel it is a little hacky and very inefficient in one particular detail and there must be a better way. In the Card class init, I am setting an attribute order which is just a list of the members of the CardValue enum, in the order they are declared.

While I could just use enum.auto(), or manually set the values as integers and hence simply use those for comparisons (still with eq etc), that would lose me the ability to make use of the enum values to form the human readable card value names.

Is there a better way to go about this than the below code? I feel like the comparison operators should be in the enum itself, but so far have not worked out how to get that to work with my meaningful values for the members of the enum.

from enum import Enum
from functools import total_ordering

class CardValue(Enum):
    TWO = "2"
    THREE = "3"
    ...
    ...
    KING = "King"
    ACE = "Ace"    

@total_ordering
class Card:
    def __init__(self, value, suit):
        self.value = value
        self.suit = suit
        self.order = (en for en in CardValue) ### This seems wrong and very inefficient!

    def __eq__(self, other):
        if self.__class__ is other.__class__:
            return self.order.index(self.value) == self.order.index(other.value)

    def __lt__(self, other):
        if self.__class__ is other.__class__:
            return self.order.index(self.value) < self.order.index(other.value)

simulation – Simulating a custom card game in Fate Core

It’s actually a thing with a history worth learning from! But before I start in on it, let me open with a warning: don’t do The Poker Episode. You know the one. We’re hunting down the ten secret keys to Dark Stobolous’s planet-cracker all across the galaxy, this one crashed into the Chance Domes of Rota X, and in the opening shot it turns out that First Lieuteneral Twilliam has been sharping cards all over the breakrooms but after the closing credits we will never mention it again?

Fate Core assumes that the things that define PCs will stay relevant over the course of a campaign. If you take characters made without the knowledge that this card game exists or that gambling is a thing and throw them into a scenario where they have to engage with it, odds are none of their stunts will be relevant and none of their aspects will apply, and all you’ll have done is make the game more boring for a little while.

Instead I’ll be assuming not that, where people were aware of the cultural significance of gambling in the setting and made characters to that effect and somebody may even have written “Free-Spirited Card Sharp” in the high concept box or something. What’s that going to look like?

A Skill Or Not A Skill

Back before Fate Core was even a thing, Spirit of the Century came out and billed itself as a game of pulp adventure, and the pulps were lousy with glitzy casinos and backroom games and fortunes turning on a single card. So they tried to drop in a Gambling skill.

This is a setup with an entire extra slant on the skill pyramid and about three times the stunts you get in Fate Core, mind, and they still had trouble making Gambling really stand up on its own as a thing you’d spotlight on a character. A followup book, Strange Tales of the Century, can’t muster even a single new Gambling stunt, with all its expansion content.

Probably “not a skill” then. Fate Core suggested converting characters who were trained in skills that just dropped out, like Gambling or Leadership, to use aspects and write stunts that reflected the ways you’d use the remaining core skills to lead people or make bets. Contacts to find a likely gambler, Deceive to bluff, Empathy to read people, Notice to look out for tells, Lore to remember obscure rules, Provoke to force mistakes, Burglary to sleight some hands. (They provided examples in this wiki.)

You can use the idea of “aspects granting permission” here – a Free-Spirited Card Sharp is going to be able to use the existing core skills to gamble, casually, in a way other people can’t. Or if you’re in the Chance Domes of Rota X, where Everyone Bets Everything, anyone can gamble.

If you think you can make it work as a skill, more power to you; this advice is probably going to apply either way, but I wanted to get that out there.

The Casual Gambler

Spirit of the Century had this to say about setting up scenes to use its Gambling skill:

In practice, Gambling requires striking a balance between cool scenes and boring play. The moment that matters in a gambling scene is the last one, when everything is on the line and the last card gets turned over. (…)

Gambling scenes should usually be picked up in medias res. Give a quick rundown of who’s at the table, making sure to include their body language, before picking up with the gambler character and the fall of the dice.

What’s gambling for, though? In the casual case, getting stuff, same as Resources, but with a different spread of costs and benefits – there are things people will and won’t gamble, and instead of temporarily being short on cash as a cost to failure, you can find yourself with a shady debt to a shady character.

You can use this to introduce various aspects of the game, the concept of side bets, et cetera, before whatever big gambling scenes you have planned. Even fairly high-stakes gambles for big pricey things can be resolved in a single roll – the drama’s going to come in as Fate Points trickle into play and bonk off the various character and scene aspects in order to actually put up the necessary numbers.

The Big Gambling Scene

When you want to make gambling into the spotlight of the session, rather than a single dramatic story beat, the first thing you need to ask yourself is “what’s the point?” I mean, maybe just “win the game” is the point? Planet-cracker key and all? But other things can be the point, too. Maybe the game is a distraction for a casino heist, and staying in for as long as possible is the goal, but winning would be pretty sweet too. Maybe the game is the point, but you can’t win with Dark Stobolous’s Doom Generator broadcasting misfortune waves and you need to hold on as long as possible while your team takes it out.

The second thing you need to ask yourself is “what’s everyone else doing”? Besides First Lieuteneral Twilliam, I mean, who is through to the finals. This is going to follow from the first question. I really don’t think side bets for side bets’ sake are going to cut it, at this point.

Gambling, The Challenge

If gambling is one part of a multistep plan – say, Twilliam’s out there on the casino floor being fancy while everyone else tries to track down and gain entrance to the high rollers’ suite to bet on him for the planet-cracker – then use the rules for resolving a challenge. Give everyone something to do and Twilliam’s part will be one gambling-related roll, probably Deceive but any of those earlier mentioned stand-in skills will likely do. Just keep in mind the basic rules of challenges – repetition is highly discouraged.

Gambling, The Contest

You can intersperse rounds of a contest with cutaways to anybody else doing anything else, so if you want to know if Twilliam actually wins anything while keeping up the pace elsewhere, you might consider that. Of the gambling-related skills mentioned earlier, Deceive, Provoke, and Empathy will likely make the relevant roll in the contest, and Notice (vs. Deceive), Burglary (vs. Notice), and Lore (vs. an escalating obstacle as relevant loopholes get more obscure) can be used to jockey for advantage.

Gambling, The Conflict

Or if you want things to work like a full-on contest of wills, you can just run a social conflict, straight-up, but expand the scopes of allowable attacks and defenses to reflect the fact that you can also stress out the opposition by being good at gambling.

You might want to run it in a wider scope that encompasses multiple successive hands of cards or even entire games in a single roll on Twilliam’s part, giving people the opportunity to support them with Create an Advantage and such in between the big moments on the casino floor.

opengl – Box2D simulation doesn’t work

I previously used Box2D and it always worked fine until recently I decided to test how it would work in my custom 2D game engine, I just wanted to test the physics updates without any GUI interaction, as you can see in the below code I just try to print the plain position values of the dynamic body and it just doesn’t move. All it does is print the initialisation position I set in the initialiser, afaik all I’m trying to do it print the values in a simple loop that runs more than 60 times per frame, the box2d code doesn’t interact with the rendering in anyway. IDK what’s wrong with box 2d to run fine in a simple loop. I’m really confused why the simulation isn’t happening. Let me know if you need more info about anything.

Freefall.h

 #include <fireworks/fireworks.h>
 #include <box2d/box2d.h>

using namespace fireworks;

class FreeFall : public Fireworks
{
private:
    Window*         m_Window;
    Layer*          defaultLayer;

    b2Vec2          m_Gravity;
    const double    m_PhysicsTimeStep = 1.0f / 60.0f;
    unsigned int    m_VelocityIterations;
    unsigned int    m_PositionIterations;
public:
    b2World* world;

    b2BodyDef       groundBodyDef;
    b2Body*         groundBody;
    b2PolygonShape  groundShape;
    b2FixtureDef    groundFixtureDef;

    b2BodyDef       dynBoxBodyDef;
    b2Body*         dynBoxBody;
    b2PolygonShape  dynBoxShape;
    b2FixtureDef    dynBoxFixtureDef;

    Sprite*         ground;
    Sprite*         dynBox;
public:
    FreeFall()
        : m_Gravity(b2Vec2(0.0f, -29.81f)), m_VelocityIterations(6), m_PositionIterations(2)
    {
        world = new b2World(m_Gravity);
        // Static ground body
        groundBodyDef.position.Set(0.0f, -10.0f);
        groundBody = world->CreateBody(&groundBodyDef);
        groundShape.SetAsBox(20.0f, 4.0f);
        groundFixtureDef.shape = &groundShape;
        groundFixtureDef.density = 1.0f;
        groundFixtureDef.friction = 0.3f;
        groundBody->CreateFixture(&groundFixtureDef);

        // Dynamic simulation box
        dynBoxBodyDef.type = b2_dynamicBody;
        dynBoxBodyDef.position.Set(-1.0f, 4.0f);
        dynBoxBody = world->CreateBody(&dynBoxBodyDef);
        dynBoxShape.SetAsBox(2.0f, 2.0f);
        dynBoxFixtureDef.shape = &dynBoxShape;
        groundFixtureDef.density = 1.5f;
        dynBoxFixtureDef.friction = 0.25f;
        dynBoxBody->CreateFixture(&dynBoxFixtureDef);
    }

    ~FreeFall()
    {
        delete defaultLayer;
        delete world;
    }

    // Runs once per initialisation
    void init() override
    {
        m_Window = createWindow("Freefall physics sim", 800, 600);
        glClearColor(0.8, 0.8f, 0.2f, 1.0f);

       

    }
    // Runs once per second
    void tick() override { }

    // Runs 60 times per second
    void update() override { }

    // Runs as fast as possible
    void render() override
    {
        //Physics Update
        world->Step(m_PhysicsTimeStep, m_VelocityIterations, m_PositionIterations);
        b2Vec2 dynPos = dynBoxBody->GetPosition();
        std::cout << "dynnamic Box2d Box position X : " << dynPos.x << " and Y is : " << dynPos.y << std::endl;
    }
};

MainGame.cpp

#include "physics-sims/Freefall.h"

int main()
{
    FreeFall game;
    game.start();
    return 0;
}

This is the output I get :
dynnamic Box2d Box position X : -1 and Y is : 4 for as long as the loop runs, this is soo infuriating, IDK what’s breaking what.

performance – Project Euler #645 — speed up Monte-Carlo simulation in Python

I am trying to solve Q645. While the logic used for my code seems to be appropriate, the code itself is way too slow for the large number required in this question. May I ask for suggestion(s) to improve my code’s performance?

The question is as in the link: https://projecteuler.net/problem=645

My Python code is as follow:

def Exp(D):
    day_list = (0)*D
    num_emperor = 0
    while all((d == 1 for d in day_list)) == False:
        #the birthday of the emperors are independent and uniformly distributed throughout the D days of the year
        bday = np.random.randint(0,D)
        day_list(bday) = 1
        num_emperor+=1
        #indices of d in day_list where d == 0
        zero_ind = (i for i,v in enumerate(day_list) if v == 0)
        for ind in zero_ind:
            try:
                if day_list(ind-1) and day_list(ind+1) == 1:
                    day_list(ind) = 1
            except IndexError:
                if ind == 0:
                    if day_list(-1) and day_list(1) == 1:
                        day_list(0) = 1
                elif ind == len(day_list)-1:
                    if day_list(len(day_list)-2) and day_list(0) == 1:
                        day_list(len(day_list)-1) = 1
    return num_emperor

def my_mean(values):
    n = 0
    summ = 0.0
    for value in values:
        summ += value
        n += 1
    return summ/n

def monte_carlo(iters, D):
    iter = 0
    n_emperor = 0
    while iter < iters:
        n_emperor = Exp(D)
        yield n_emperor
        iter += 1

avg_n_emperor = my_mean(monte_carlo(iters,D))
print(avg_n_emperor)

And my logic is as follow:

For the day_list inside the Exp(D) function, where D is the number of days in a year, zeros mean no holiday, and ones mean holiday. Initially the day_list is all zeros since there is no holiday to begin with.

The rules of defining a random day (d) as a holiday is as follow:

  1. At the beginning of the reign of the current Emperor, his birthday is declared a holiday from that year onwards.

  2. If both the day before and after a day d are holidays, then d also becomes a holiday.

I then subsequently implement the rules stated for the question, to gradually add holidays (ones) into the day_list. After num_emperor number of emperors, all the days (d) in day_list will become 1, i.e. all days will become holiday. This is the point to quit the while_loop in Exp(D) function and count the number of emperors required. To get the average number of emperors required for all the days to become holidays (avg_n_emperor), I then apply the monte-carlo method.

For my current code, the time takes is as follow:

avg_n_emperor = my_mean(monte_carlo(iters=100000,D=5)) #6-7 seconds

avg_n_emperor = my_mean(monte_carlo(iters=1000000,D=5)) #about 62 seconds

in which the time takes increase approx. linearly with the iters.

However,

avg_n_emperor = my_mean(monte_carlo(iters=1000,D=365)) #about 68 seconds

already takes about 68 seconds, and the question is asking for D=10000. Not to mention that the iters required for the answer to be accurate within 4 digits after the decimal points (as required by the question) would be much larger than 1000000 too…

Any suggestion(s) to speed up my code would be appreciated! 🙂

performance – N-body simulation in javascript

I’ve worked on an n-body simulation using javascript in the past few years and recently I decided to update it to use a barnes-hut tree to speed up performance. It currently handles about 2,000 bodies at roughly 30fps, which I believe is quite underwhelming despite being double the performance it had before the tree structure was implemented.

The following code is contained in the PhysicsWorld implementation. I’d like some insight on ways I could improve the performance of this class, and ideas on how I can take responsibility out of it, as I believe it’s doing too much currently but I’m not yet entirely sure how to move code to other classes.

Note that most of the time spent here is wasted on PhysicsWorld.traverseTree() so it’s a good idea to focus effort there.

Thanks in advance.

class PhysicsWorld {
    constructor (debugRenderer) {
        this.G = 0.01
        this.iterations = 1
        this.bodies = ()
        this.intersections = ()
        this.garbage = ()
        this.barnesHutTree = null
        this.theta = 0.5
        this.computationsPerIteration = 0
        this.debugRenderer = debugRenderer
        this.adaptiveDomainBoundary = new AABB(new Point(0, 0), new Size(innerWidth, innerHeight))
        this.enforceSquareNodes = false
        this.useAdaptiveDomainBoundary = true
    }
    createBarnesHutTree () {
        const position = new Point(0, 0),
            size = new Size(window.innerWidth, window.innerHeight),
            boundary = new AABB(position, size)
        this.barnesHutTree = new BarnesHutTree(this.adaptiveDomainBoundary, this.debugRenderer)

        for (let body of this.bodies) {
            this.barnesHutTree.insert(body)
        }
    }
    applyVelocityToPosition () {
        let index = 0

        let position = new Point(Infinity, Infinity),
            size = new Size(-Infinity, -Infinity)

        for (const body of this.bodies) {

            body.userData.bodiesArrayIndex = index

            body.pastPosition.x = body.position.x
            body.pastPosition.y = body.position.y
            body.position.x += body.velocity.dx / this.iterations
            body.position.y += body.velocity.dy / this.iterations

            // Store maximum and minimum body positions to use as adaptive domain boundary:
            position.x = Math.min(body.position.x, position.x)
            position.y = Math.min(body.position.y, position.y)
            size.width = Math.max(body.position.x, size.width)
            size.height = Math.max(body.position.y, size.height)

            ++index
        }

        // Use stored body positions for adaptive domain boundary:
        const margin = 1 // fixes floating point errors from preventing body from being added to tree
        position.x -= margin
        position.y -= margin
        size.width = size.width - position.x + margin / 2
        size.height = size.height - position.y + margin / 2
        if (this.enforceSquareNodes) {
            size.width = Math.max(size.width, size.height)
            size.height = Math.max(size.width, size.height)
        }
        if (this.useAdaptiveDomainBoundary)
            this.adaptiveDomainBoundary = new AABB(position, size)
    }
    traverseTree () {
        this.computationsPerIteration = 0
        for (const body of this.bodies) {
            this.barnesHutTree.forEachNode(node => {
                const distanceX = node.centerOfMass.x - body.position.x,
                    distanceY = node.centerOfMass.y - body.position.y,
                    distanceSquare = distanceX * distanceX + distanceY * distanceY,
                    averageNodeSideLength = Math.max(node.boundary.size.width, node.boundary.size.height),
                    averageNodeSideLengthSquare = averageNodeSideLength * averageNodeSideLength,
                    thetaSquare = this.theta * this.theta,
                    isNodeFarEnoughToApproximateAsSingleBody = averageNodeSideLengthSquare / distanceSquare < thetaSquare

                ++this.computationsPerIteration

                if (isNodeFarEnoughToApproximateAsSingleBody || node.isEndNode) {

                    if (node.isEndNode && node.isPopulated) {
                        if (node.body === body) return false

                        const radiiSum = node.body.shape.radius + body.shape.radius,
                            radiiSumSquare = radiiSum * radiiSum

                        if (radiiSumSquare > distanceSquare) {
                            if (body.collidable && node.body.collidable) 
                                this.markAsIntersection(body, node.body)

                            return false
                        }
                    }

                    const distance = Math.sqrt(distanceSquare),
                        force = this.G * ((body.mass * node.mass) / distanceSquare),
                        forceByIteration = force / this.iterations

                    body.velocity.dx += (forceByIteration / body.mass) * distanceX / distance
                    body.velocity.dy += (forceByIteration / body.mass) * distanceY / distance

                    return false
                } else return true
            })
        }
    }
    markAsIntersection (bodyA, bodyB) {
        const isAlreadyIntersectingWithAnotherBody = {
                bodyA: bodyA.contact !== null,
                bodyB: bodyB.contact !== null
            },
            numberOfBodiesWithPreviousIntersections = 
                isAlreadyIntersectingWithAnotherBody.bodyA + 
                isAlreadyIntersectingWithAnotherBody.bodyB

        switch (numberOfBodiesWithPreviousIntersections) {
            case 0:
                const intersection = (bodyA, bodyB),
                    index = this.intersections.length

                bodyA.contact = index
                bodyB.contact = index

                this.intersections.push(intersection)
                break
            case 1:
                if (isAlreadyIntersectingWithAnotherBody.bodyA) {
                    bodyB.contact = bodyA.contact
                    this.intersections(bodyA.contact).push(bodyB)
                } else {
                    bodyA.contact = bodyB.contact
                    this.intersections(bodyB.contact).push(bodyA)
                }
                break
            case 2:
                if (bodyA.contact !== bodyB.contact) { // no need to do anything if they're already in the same intersection
                    const oldIndex = bodyB.contact
                    for (const bodyC of this.intersections(bodyB.contact)) {
                        this.intersections(bodyA.contact).push(bodyC)
                        bodyC.contact = bodyA.contact
                    }
                    this.intersections(oldIndex).length = 0 // can't remove array element otherwise all indexes after it would get messed up
                }
                break
        }
    }
    mergeIntersectingBodies () {
        for (const intersection of this.intersections) {
            if (intersection.length > 0) {
                let bodyA = intersection(0)
                for (let i = 1; i < intersection.length; ++i) {
                    let bodyB = intersection(i)

                    if (bodyB.mass > bodyA.mass) {
                        const buffer = bodyA
                        bodyA = bodyB
                        bodyB = buffer
                    }

                    bodyA.position.x = (bodyA.position.x * bodyA.mass + bodyB.position.x * bodyB.mass) / (bodyA.mass + bodyB.mass)
                    bodyA.position.y = (bodyA.position.y * bodyA.mass + bodyB.position.y * bodyB.mass) / (bodyA.mass + bodyB.mass)

                    bodyA.velocity.dx = (bodyA.velocity.dx * bodyA.mass + bodyB.velocity.dx * bodyB.mass) / (bodyA.mass + bodyB.mass)
                    bodyA.velocity.dy = (bodyA.velocity.dy * bodyA.mass + bodyB.velocity.dy * bodyB.mass) / (bodyA.mass + bodyB.mass)

                    bodyA.shape.volume += bodyB.shape.volume / bodyA.density

                    this.garbage.push(bodyB.userData.bodiesArrayIndex)

                    bodyB.contact = null
                }
                bodyA.contact = null
            }
        }
        this.intersections.length = 0
    }
    collectGarbage () {
        let numberOfDeletedItems = 0

        const sortAlgorithmDescendingOrder = (a, b) => b - a
        this.garbage.sort(sortAlgorithmDescendingOrder) // prevents from having to shift the index of subsequent items

        for (const index of this.garbage) {
            const body = this.bodies(index)

            this.bodies.splice(index, 1)
            body.destroy()

            ++numberOfDeletedItems
        }
        this.garbage.length = 0
    }
    step () {
        for (let i = 0; i < this.iterations; ++i) {
            
            this.applyVelocityToPosition()

            this.createBarnesHutTree()
            this.traverseTree()

            this.mergeIntersectingBodies()
            this.collectGarbage()
        }
    }
}

If you want to take a look on the entire project please check out its github repository.

plotting – Erro in Monte Carlo simulation. Write: file is not a string, stream or lists of strings and streams

I’ve been trying to code a kinetic monte carlo simulation for a system of consecutive chemical reactions using a fixed (repaired) code that the user flinty gives me in the last questions I post:
Strange error not found in a monte carlo code, Infinite expression 1/0 encountered.
In that post, the original code for one chemical reaction fail because of the use of Import and FormatType function, and it was fixed using Reap and Sow.
Now I have a code that works perfectly using again Import and FormatType but it fails when I try to implement Reap and Sow on it.

The aim is to obtain the concentration of every chemical compound involved in a system of 2 chemical reactions.

The Reason why I don’t want to use Import and FormatType is because when I try to add mor chemical reactions it fails, I mean, If you change something of that code, it gets broken also.

I don’t undesrtand whats wrong if the Reap and Sow structure is similar to the code for one chemical reaction

Here the Code that works using Reap and Sow for One chemical Reaction, Dimerization, 2M–>D

Mo = 5; (*mol/L*) (*Initial concentration of M, Monomer*)
DDo = 0; (*mol/L*) (*Initial concentration of D, Dimer*)
R = 1.98 ;(*J/(mol K), Ideal Gas constant*)
Tk = 383; (*K, Tmperature of the reaction in kelvin*)
kdim = 2.52*10^4*E^(-22347/(R*Tk)); (*L/(mol s), kinetic constant*)
tf = 30000*60*60; (*s, final time in seconds, 30'000 hours of reaction*)

M = Mo/(1 + 2*kdim*t*Mo); (*Analytical solution for the ODE system, for M*)

DD = (2*kdim*t* Mo^2 + DDo + 2*kdim*t* Mo*DDo)/(1 + 2*kdim*t*Mo); (*Analytical solution for the ODE system, for D*)

p8 = Plot({M, DD}, {t, 0, tf}, PlotLabel -> "Concentration vs Time", 
  AxesLabel -> {"Time (s)", "Concentration mol/L"}, PlotStyle -> {Red, Green}, PlotLegends -> {"(M)", "(D)"});

(*Kinetic Monte Carlo loop starts here*)

Mo = 5;
R = 1.98;
Tk = 383;
kdim = 2.52*10^4*Exp(-22347/(R*Tk));

tf = 30000*60*60

M = Mo/(1 + 2*kdim*t*Mo);

p6 = Plot(M, {t, 0, tf}, PlotStyle -> Red, 
   AxesLabel -> {"time (s)", "Concentration (mol/L)"}, 
   PlotLegends -> {"(M)"}, PlotLabel -> "Concentration of M");

NAV = 2*10^3; (*Avogadro number times system volume*)

XM = Mo*NAV; (*Number of M molecules*)

kdimMC = 2*(kdim/NAV); (*Kinetic Monte carlo constant*)

tb = 0.0;

M = Mo;

Niter = 0;

toc = Timing(
   result = 
    Reap(While(tb <= tf, Sow({tb, M}); Niter = Niter + 1; 
       a0 = kdimMC*XM*(XM - 1)/2; XM = XM - 2; M = XM/NAV; 
       tb = tb - Log(RandomReal())/a0))((2, 1)));

Export("Result3.dat", result);

Print("iterations for M:", Niter, "ntiming for M", First(toc))

p9 = ListPlot(result, DataRange -> {0, tf}, PlotStyle -> Blue, 
   PlotLegends -> {"kMC (M)"});

Mo = 5;
DDo = 0;
R = 1.98;
Tk = 383;

kdim = 2.52*10^4*Exp(-22347/(R*Tk));

tf = 30000*60*60

M = Mo/(1 + 2*kdim*t*Mo);

DD = (2*kdim*t*Mo^2 + DDo + 2*kdim*t*Mo*DDo)/(1 + 2*kdim*t*Mo);

p7 = Plot(DD, {t, 0, tf}, PlotStyle -> Green, 
   AxesLabel -> {"time (s)", "Concentration (mol/L)"}, 
   PlotLegends -> {"(D)"}, PlotLabel -> "Concentration of D");

NAV = 2*10^3; (*Avogadro number times system volume*)

XM = Mo*NAV; (*Initial number of molecules of M*)

XDD = DDo*NAV/2; (*It's divide by 2 because there must be half of the number of molecules of D regarding the number of M molecules, due to D is a Dimer of M*)

kdimMC = 2*(kdim/NAV); (*kinetic monte carlo constant*)

tb = 0.0;

M = Mo;

DD = DDo;

Niter = 0;

toc = Timing(result = Reap(While(tb <= tf, Sow({tb, DD});
       Niter = Niter + 1;
       a0 = kdimMC*XM*(XM - 1)/2;
       XM = XM - 2;
       XDD = XDD + 1;
       M = XM/NAV;
       DD = Mo - M;(*by stoichiometry, we can rewrite D concentration like this, just a mass balance*)
       tb = tb - Log(RandomReal())/a0))((2, 1)));

Export("Result3.dat", result);

Print("iterations for D:", Niter, "ntiming for D", First(toc))

p10 = ListPlot(result, DataRange -> {0, tf}, PlotStyle -> Orange, 
   PlotLegends -> {"kMC (D)"});

Show(p8, p9, p10, PlotRange -> {{0, tf}, {0, 5}})

Here the code that works using Import and FormatType for the consecutive reactions A–>B; B–>C

A0 = 1.00;(*mol/L, initial concentration of A*)
B0 = 0.00;(*mol/L, initial concentration of B*)
C0 = 0.00;(*mol/L, Initial concentration of C*)

k1 = 2*10^(-4);(* s^-1, kinetic constant for the reaction A-->B*)
k2 = 1*10^(-4);(* s^-1, kinetic constant for the reaction B-->C*)

tf = 3600*5;(*s, final time*)

ta = {t, 0, tf}; (*Interval of time*)

A = A0*Exp(-k1*t); (*Analytical solution for A*)

B = B0*Exp(-k2*t) + (A0*k1/(k2 - k1))*(Exp(-k1*t) - Exp(-k2*t)); (*Analytical solution for B*)

Cs = A0 + B0 + C0 - A - B; (*Analytical solution for C*)

p1 = Plot(A, ta ,  PlotStyle -> {Dashed, Red});
p2 = Plot(B, ta , PlotStyle -> {Dashed, Blue});
p3 = Plot(Cs, ta , PlotStyle -> {Dashed, Gray});

NAV = 2*10^3; (*Avogadro number times volume system*)

XA = A0*NAV; (*Number of molecules of A*)
XB = B0*NAV; (*Number of molecules of B*)
XC = C0*NAV; (*Number of molecules of C*)

k1MC = k1; (*Kinetic monte carlo constant*)
k2MC = k2; (*Kinetic Monte Carlo constant*)

Maxr = 2;  (*Maximum number of reactions*)

kmax = Ceiling(Log2(Maxr)); (*Number of Steps, or movements, needed to occur one of the two possible reactions*)

VR = Table(0.0, {i, 1, 1 kmax + 1}, {j, 1, 2^kmax}); (*Vector of Reactions*)

tb = 0.0;
Niter = 0;

file = OpenWrite("Result.dat", FormatType -> OutputForm)
toc = Timing(
   While(tb < tf,
    Niter = Niter + 1;
    
    
    VR((kmax + 1, 1)) = k1MC*XA; (*Reaction 1*)
    
    VR((kmax + 1, 2)) = k2MC*XB; (*Reaction 2*)

    For(i = 1, i <= kmax + 1, i++,
     For(j = 1, j <= 2^(kmax - i), j++,
       VR((kmax + 1 - i, j)) = 
         VR((kmax + 2 - i, 2*j)) + VR((kmax + 2 - i, 2*j - 1));
       ););
    
    a0 = VR((1, 1)); (*Total reaction rate*)
    ra0 = RandomReal()*a0;
    mv = 1;
    
    For(i = 1, i <= kmax, i++,
     If(ra0 <= VR((i + 1, 2*mv - 1)), mv = 2*mv - 1, 
      ra0 = ra0 - VR((i + 1, 2*mv - 1)); mv = 2*mv));
    
    If (mv == 1, XA = XA - 1; XB = XB + 1;);
    If (mv == 2, XC = XC + 1; XB = XB - 1;);

    A = XA/NAV;
    B = XB/NAV;
    Cs = XC/NAV;

    Write(file, tb, " ", A, " ", B, " ", Cs);
    tb = tb - Log(RandomReal())/a0));

Close(file);
Print(Niter);
Result = Import("Result.dat");

tb = Result((All, 1));
A = Result((All, 2));
B = Result((All, 3));
Cs = Result((All, 4));

p4 = ListPlot(Transpose@{tb, A}, PlotStyle -> {Dotted, Red});
p5 = ListPlot(Transpose@{tb, B}, PlotStyle -> {Dotted, Blue});
p6 = ListPlot(Transpose@{tb, Cs}, PlotStyle -> {Dotted, Gray});
Show(p1, p2, p3, p4, p5, p6, PlotRange -> {0, 1})

Here the Broken Code for the system A–>B; B–>C using Reap and Sow

A0 = 1.00; (*mol/L*)
B0 = 0.00;(*mol/L*)
C0 = 0.00;(*mol/L*)

k1 = 2*10^-4; (*s^-1*)
k2 = 1*10^-4; (*s^-1*)

tf = 3600*5; (*s*)
ta = {t, 0, tf};
A = A0*Exp(-k1*t);
B = B0*E^(-k2*t) + (A0*k1)/(k2 - k1)*(E^(-k1*t) - E^(-k2*t));
Cs = A0 + B0 + C0 - A - B;

p1 = Plot(A, ta, PlotStyle -> {Dashed, Red});
p2 = Plot(B, ta, PlotStyle -> {Dashed, Blue});
p3 = Plot(Cs, ta, PlotStyle -> {Dashed, Gray});

NAV = 2*10^3; (* (molecules L / mol) Avogadro number times Volume*)

XA = A0*NAV;  (* Number of A molecules *)
XB = B0*NAV;   (* Number of B molecules *)
XC = C0*NAV;   (* Number of C momlecules *)

k1MC = k1; (*Kinetic monte carlo reaction constants*)
k2MC = k2;

Maxr = 2;  (* Number of reactions *)

kmax = Ceiling(
   Log2(Maxr));   (*Number of steps (or movements) to determine which 
chemical reactions proceeds, according to the binary tree*)

VR = Table(0.0, {i, 1, kmax + 1}, {j, 1, 2^kmax});

tb = 0.0;

Niter = 0;

toc = Timing(result = Reap(
      While(tb < tf, Sow({tb, , A, B, Cs});
       Niter = Niter + 1;
       
       VR((kmax + 1, 1)) = k1MC*XA; (*R1, reaction 1, 
       VR is Vector of reactions*)
       
       VR((kmax + 1, 2)) = k2MC*XB; (*R2, reaction 2*)
       
       For(i = 1, i <= kmax + 1, i++,
        For(j = 1, j <= 2^(kmax - i), j++,
          
          VR((kmax + 1 - i, j)) = 
            VR((kmax + 2 - i, 2*j)) + VR((kmax + 2 - i, 2*j - 1));
          ););
       a0 = VR((1, 1)); (*Total vector of reactionr, VR*)
       ra0 = RandomReal()*a0;
       mv = 1;
       
       For(i = 1, i <= kmax, i++,
        If(ra0 <= VR((i + 1, 2*mv - 1)), mv = 2*mv - 1,
         ra0 = ra0 - VR((i + 1, 2*mv - 1)); mv = 2*mv));
       
       If(mv == 1, XA = XA - 1; XB = XB + 1;);
       If(mv == 2, XC = XC + 1; XB = XB - 1;);
       
       A = XA/NAV; (*Concentration of A*)
       B = XB/NAV;
       Cs = XC/NAV;
       
       Write(file, tb, " ", A, " ", B, " ", Cs);
       
       tb = tb - Log(RandomReal())/a0))((2, 1)));

Export("Result3.dat", result);
Print("iterations:", Niter, "ntiming", First(toc))

p4 = ListPlot(result, DataRange -> {0, tf});

Close(file);
Print(Niter);

Show(p1, p2, p3, p4, PlotRange -> {0, 1})

OutPut of the Broken Code using Reap and Sow

Game of life GUI simulation using java.swing

I am fairly new to GUI programming and I am interested in improving my code from every aspect. Performance, security, readability, conciseness, and the look-and-feel aspects are all important to me.

Here is the non-GUI part of the code:

import java.util.Random;

public class Universe
{
    private int generation;
    private int alive;
    private boolean()() currentGeneration;
    private boolean()() nextGeneration;
    private Random random;

public Universe(int height, int width, int seed, String pattern)
{
   this.currentGeneration = new boolean(height)(width);
   if (pattern.equalsIgnoreCase("Random"))
   {
       random = new Random(seed);
       for (int i = 0; i < height; i++) 
       {
           for (int j = 0; j < width; j++) 
           {
               currentGeneration(i)(j) = random.nextBoolean();
           }
       }
   }
   else if (pattern.equalsIgnoreCase("glider"))
   {
       getGlider(currentGeneration);
   }
   //to add more cases here
    nextGeneration = generateNextGeneration(currentGeneration);
    generation = 1;
    alive = calculateAlive(currentGeneration);
}

//Getters and instance methods

int getGeneration()
{
    return generation;
}
int getAlive()
{
    return alive;
}
boolean()() getCurrentGeneration()
{
    return currentGeneration;
}
boolean()() getNextGeneration()
{
    return nextGeneration;
}
void moveToNextState()
{
    boolean()() temp = generateNextGeneration(nextGeneration);
    currentGeneration = nextGeneration;
    nextGeneration = temp;
    alive = calculateAlive(currentGeneration);
    generation++;
}
void reset(int h, int w, int seed)
{
    this.currentGeneration = new boolean(h)(w);
    random = new Random(seed);
    for (int i = 0; i < h; i++)
    {
        for (int j = 0; j < w; j++)
        {
            this.currentGeneration(i)(j) = random.nextBoolean();
        }
    }
    nextGeneration = generateNextGeneration(currentGeneration);
    generation = 1;
    alive = calculateAlive(currentGeneration);
}

//Utility methods

static int calculateNeighbours(boolean()() grid, int row, int column)
{
    int neighbours = 0, r, c;
    int N = grid.length;
    int M = grid(0).length;

    for (int p = -1; p <= 1; p++)
    {
        for (int m = -1; m <= 1; m++)
        {
            r = row + p;
            c = column + m;

            if (r < 0)
                r = N - 1;
            if (r > N - 1)
                r = 0;
            if (c < 0)
                c = M - 1;
            if (c > M - 1)
                c = 0;

            if (grid(r)(c) && (p != 0 || m != 0))
                neighbours++;
        }
    }
    return neighbours;
}

static int calculateAlive(boolean()() grid)
{
    int alive = 0;
    for (int i = 0; i < grid.length; i++)
    {
        for (int j = 0; j < grid(0).length; j++)
        {
            if (grid(i)(j))
                alive++;
        }
    }
    return alive;
}

static boolean()() generateNextGeneration(boolean()() currentGeneration)
{
    int N = currentGeneration.length;
    int M = currentGeneration(0).length;
    boolean()() nextGeneration = new boolean(N)(M);
    int neighbours;
    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < M; j++)
        {
            neighbours = calculateNeighbours(currentGeneration, i, j);

            if (neighbours == 3 || (currentGeneration(i)(j) && neighbours == 2))
                nextGeneration(i)(j) = true;
            else
                nextGeneration(i)(j) = false;
        }
    }
    return nextGeneration;
}

static boolean()() generateNthGeneration(boolean()() currentGeneration, int X)
{
    if (X == 0)
        return currentGeneration;
    else
        return generateNthGeneration(generateNextGeneration(currentGeneration), X - 1);
}
static void printGeneration(boolean()() generation)
{
    for (int i = 0; i < generation.length; i++)
    {
        for (int j = 0; j < generation(0).length; j++)
            System.out.print(generation(i)(j)? "O" : " ");
        System.out.println();
    }

}

static void getGlider(boolean currentGeneration()())
{
    for(int i = 0; i < 60; i++)
    {
        for (int j =0; j < 90; j++)
            currentGeneration(i)(j) = false;
    }
    currentGeneration(1)(3) = true;
    currentGeneration(2)(3) = true;
    currentGeneration(3)(3) = true;
    currentGeneration(2)(1) = true;
    currentGeneration(3)(2) = true;
}

}

And here is the GUI part of the code:

import javax.swing.*;
import java.awt.*;
import java.awt.event.ActionListener;

class Cells extends JPanel
{
    boolean()() grid;
    int h, w;
    Cells(boolean()() grid)
    {
       this.grid = grid;
        h = grid.length;
        w = grid(0).length;
    }
    {
        setBounds(50, 20, 961, 620);
        setBackground(Color.DARK_GRAY);
    }
    @Override
    public void paintComponent(Graphics g)
    {
        super.paintComponent(g);
        Graphics2D g2 = (Graphics2D) g;
        //g2.setColor(Color.BLUE);

        for (int x = 0; x < w * 10; x+=10)
        {
            for (int y = 0; y < h * 10; y+=10)
            {
                if (grid(y/10)(x/10))
                {
                    g2.setColor(Color.BLUE);
                    g2.fillRect(x, y, 10, 10);
                }
                else
                {
                    g2.setColor(Color.gray);
                    g2.drawRect(x, y, 10, 10);
                }
            }
        }
    }
}
public class GameOfLife extends JFrame
{
    public final int height = 60;
    public final int width = 90;

    Universe universe = new Universe(height, width, (int) Math.random(), "Random");

    Cells cells = new Cells(universe.getCurrentGeneration());

    JLabel generationLabel = new JLabel("Generation#" + universe.getGeneration());
    JLabel aliveLabel = new JLabel("Alive: " + universe.getAlive());

    JButton resetButton, speedUpButton, slowDownButton;
    JToggleButton playToggleButton;

    String() items = {"random", "Glider", "Gun", "Spaceship", "Beacon", "Pulsar"}; //to be added
    JComboBox patterns = new JComboBox(items); //to be added

    ActionListener repaint = e ->
    {
        universe.moveToNextState();
        generationLabel.setText("Generation #" + universe.getGeneration());
        aliveLabel.setText("Alive: " + universe.getAlive());
        cells.grid = universe.getCurrentGeneration();
        repaint();
        setVisible(true);
    };

    int speed = 100;
    Timer timer = new Timer(speed, repaint);

    public GameOfLife()
    {
        setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        setSize(1000, 700);
        setResizable(false);
        setLocationRelativeTo(null);
        setLayout(null);
        setBackground(Color.darkGray);
        getContentPane().setBackground(Color.darkGray);

        generationLabel.setName("GenerationLabel");
        aliveLabel.setName("AliveLabel");
        resetButton = new JButton("Reset");
        resetButton.setName("ResetButton");
        playToggleButton = new JToggleButton("Pause");
        playToggleButton.setName("PlayToggleButton");
        speedUpButton = new JButton("Speed+");
        slowDownButton = new JButton("Speed-");

        add(cells);
        addLabels();
        addButtons();
        addFunctionality();

        timer.start();

        setVisible(true);
    }

    void addLabels()
    {
        JPanel labels = new JPanel()
        {
            {
                setBounds(50, 636, 200, 40);
                setLayout(new BoxLayout(this, BoxLayout.Y_AXIS));
                setBackground(Color.DARK_GRAY);
                generationLabel.setForeground(Color.LIGHT_GRAY);
                aliveLabel.setForeground(Color.LIGHT_GRAY);
                add(generationLabel);
                add(aliveLabel);
            }
        };
        add(labels);
    }
    void addButtons()
    {
        JPanel buttons = new JPanel()
        {
            {
                setBounds(250, 636, 500, 40);
                setLayout(new BoxLayout(this, BoxLayout.X_AXIS));
                setBackground(Color.DARK_GRAY);
                resetButton.setForeground(Color.darkGray);
                playToggleButton.setForeground(Color.darkGray);
                speedUpButton.setForeground(Color.darkGray);
                slowDownButton.setForeground(Color.darkGray);

                resetButton.setBackground(Color.LIGHT_GRAY);
                playToggleButton.setBackground(Color.LIGHT_GRAY);
                speedUpButton.setBackground(Color.LIGHT_GRAY);
                slowDownButton.setBackground(Color.LIGHT_GRAY);

                add(resetButton);
                add(playToggleButton);
                add(speedUpButton);
                add(slowDownButton);
            }
        };
        add(buttons);
    }
    void addFunctionality()
    {
        playToggleButton.addActionListener(e ->
        {
             if (playToggleButton.getText().equals("Play") && !timer.isRunning())
             {
                 timer.restart();
                 playToggleButton.setText("Pause");
             }
             else if (playToggleButton.getText().equals("Pause") && timer.isRunning())
             {
                 timer.stop();
                 playToggleButton.setText("Play");
             }
        });
        speedUpButton.addActionListener(e ->
        {
            if (speed == 0)
            {}
            else
                timer.setDelay(speed -= 50);
        });
        slowDownButton.addActionListener(e -> timer.setDelay(speed += 50));
        resetButton.addActionListener(e -> universe.reset(height, width, (int) Math.random()));
    }
    public static void main(String() args)
    {
        new GameOfLife();
    }
}

Any comment, even if it was on one small issue with my code will be appreciated.

open source – Travelling Salesman Problem simulation

I apologise if I’m wrong to ask this here but I couldn’t think of a better forum.

I am learning about the travelling sales man algorithm and wish to implement a demonstration.

My initial thought was to develop it in a webpage but as I am natively a PHP developer, I feel it might not be the best solution to use javascript/jQuery to animate characters from point to point through a graph.

I wonder if there is an open-source 2D game engine I can use and write PHP/javascript logic for my algorithm implementation?

simulation – modeling DNS server CPU utilization on simulator framework

I’m creating my simulation framework, the DNS system. I have two related questions about the name server.

  1. I would like to know how much time the Authoritative Name server spends time to resolve one packet when the packet arrives at the server?.i

  2. I’m modeling the resource utilization at the DNS server in my simulator, I need to specify many parameters like maximum/avg/minimum CPU usage per one process. Those parameters need to be gotten from the real server. However, sadly, I don’t know where is I can get from? It would be nice if you guide me on how to get them.

Because I’m simulating the Authoritative Name server in SimPy, but it seems like I need some specific values from the real server to define to my simulator.