Firstly, to solve your problem in full would be a rather large post so I will focus on the methodology that you should follow.
As others have said, you would use some kind of filtering method to solve your problem; particle filtering, Kalman filtering, they are all the same thing but different variants that have different strengths. You could also use Markov Model methods that operate in discrete space rather than continuous space.
Firstly let's define the task: 3D (X, Y positions and yaw angle) localisation using an a priori map. This is different from SLAM (Simultaneous Localisation and Mapping) because you already have a map.
That being said, you could probably use an existing slam algorithm with a high-belief prior set by the map you have.
This leads to the definition of the states you are interested in (called the state vector):
$$\mathbf{x}_{k} = \big[ \mathbf{p}^\top, \mathbf{v}^\top, \psi \big]^\top$$
Where $\mathbf{x_k}$ is the state vector at discrete time-step $k$, $\mathbf{p}$ and $\mathbf{v}$ are the position and velocity, and $\psi$ is the yaw angle of the robot (all measured with respect to, and expressed in a global frame $\{G\}$). Velocity is added there because it is necessary for the prediction stage of the filter.
Second let's define the observations you can make (I exclude sensor measurements unnecessary for solving this 3DOF problem):
- Range measurements of the relative position to the wall ahead of the robot: $y_u(\mathbf{x_k})$
- Inertial measurements of the motion of the body:
- Acceleration (2 axes, x and y): $\mathbf{a}_k$
- Angular rotation rate (1 axis, z) $\mathbf{\omega}_k$
- Magnetic field strength (1 axis, z) $m_k$
- All measurements are made in the body frame $\{B\}$ which is attached to the IMU (note there is a transform between the ultrasound and body which I am not addressing, but should be addressed in an actual solution).
Now the issue is how do you use a filtering methodology to localise your robot?
Regardless of your choice of filter you will need a model that describes how the robot evolves through time (the process model).
Then you need a measurement model for each of your sensors, which describes what measurement is expected given an arbitrary state $\mathbf{x}$ of the robot.
For such a problem, your measurement models would describe (1) the ultrasound range, (2) the accelerometer, (3) the rate gyro, (4) the magnetometer.
Because of the high update rates of the IMU it is commonly used in the process model directly so that the update stage of the filter is not used too frequently (it is more computationally expensive than the predict stage). I'll refer you to equations (3) - (8) in this paper for this process model (in this application there is no need for quaternions so I'd advise to use Euler angle representations of the orientation instead).
Additionally, you must model the bias states in the IMU, which then adds another two states to your state vector:
$$\mathbf{x}_{k} = \big[ \mathbf{p}^\top, \mathbf{v}^\top, \psi, \mathbf{b}_a^\top, b_\omega \big]^\top$$
Where $\mathbf{b}_a$ and $b_\omega$ are the biases in the accelerometers and gyroscope, respectively.
The measurement model for the magnetometer can be found in equation 5.10 in this thesis.
Let's simplify your measurement model for your ultrasound, and pretend it is a laser range finder that simply returns the distance of the object (wall) it is pointing to. I assume you can do the calculation to convert the pulse times to distance. Your model is then:
$$y_u(\mathbf{x}_k) = |\mathbf{p}_{w,B}^B|_2 = |R^B_G(\mathbf{p}_{w,G}^G - \mathbf{p}_{B,G}^G)|_2$$
Where $ \mathbf{p}_{w,B}^B $ is the position of the wall that is being pointed to, with respect to (superscript) and measured from (subscript) the $\{B\}$ frame. $R^B_G$ is the rotation matrix that rotates a vector in the $\{G\}$ frame to the $\{B\}$ frame, and $\mathbf{p}_{B,G}^G$ is the position of the $\{B\}$ frame with respect to and expressed in the $\{G\}$ frame (which is the robot location represented by $p$ in the state vector).
You'll have to think about how to parameterize and calculate $\mathbf{p}_{w,G}^G$ (which is a function of the wall geometry, $p$ and $\psi$), I would think it would work with a look-up table or a solution to a set of line intersection equations representing the range beam and the walls. Either way it is discontinuous and nonlinear so this supports the use of a particle filter over a Kalman-based filter. You could possibly have a different measurement model for each wall if you really want to use a Kalman Filter, but that introduces a new data-association problem which is non-trivial.
I hope this plus NBCKLY's PF psuedocode gets you on your way.
A good resource on IMU models and noise.
Maybe of interest is this use of highly non-linear maps in PF's.