Optimizing Values to Equalize Data Sets in C++

[zbREy60M]

I have an interesting problem for you guys, it involves using C++ to optimize a few constants in a function.

I have 4 sets of pertinent data. I have stored this data in 4 separate arrays. The first array is called DMODELArray[120], the second array is DArray[120], the third array is ZArray[120], and finally the last array is DMArray[120].

Each set of data corresponds to an individual pulsar. I need to get the data in DMODELArray to match the data in DArray.

The math is set up in my code as so:

for (int i = 0; i < 120; i++) //loop to the Dmodel distance, z coordinate, and the ne for each pulsar; in order of pulsar name #
	{
                double A, B, C;
		ZArray[i] = DArray[i]*SINLATAArray[i];
		double ne = (A + B*exp(-ZArray[i]/C));
		DMODELArray[i] = (DMArray[i]/ne);
		cout << DMODELArray[i] << "   " << DArray[i] << "\n";
	}

This means that the DMODEL value needs to equal the actual D value. The difference is a double called "ne". ne is the electron density of the galaxy, while the DM is the dispersion measure. I need to take the dispersion measure (as seen in the 4th line of the loop) and modify it by the electron density. The problem is that I need to come up with a good formula for ne. The above equation is what I think the formula is going to look like, because we can assume that the electron density probably gets smaller as you get farther from the galactic plane.

My problem is that I need some way to write code to optimize the A, B, and C values. I need some sort of C++ algorithm to go through and find the appropriate values for these so that DMODELArray[i] matches DArray[i] as closely as possible.

How does one go about coding something like that?

Thank you very much

[Gaminic]

Since you don't really seem to know what kind of problem you're optimizing, it's difficult for us to weigh in with specifics.

Easiest way: simulation. Write your code, wrap it in a function, then for-loop the crap out of each variable. Print out each combination and its result, then pick the best one.

[zbREy60M]

Here is some of the data I am using:


The first 10 values DMArray are so:

4.333		
11.38	
2.38		
73.779
34.797
26.833
57.1420
2.64476	
39.903		
14.495

The first 10 values in DArray are so:

.28
1.03
.21
2.6
2.3
1.0
1.0
.156
.4
1.18

The first 10 values in ZArray are so:

-0.236435
-0.966765
-0.204467
-0.183174
-0.0846803
-0.0212911
0.0141365
-0.104302
-0.224194
0.155039

I am trying to get DMODEL to equal D as closely as possible.
DMODEL is just DM divided by ne, which is a function of Z.
I'm trying to find the best value of ne(function of Z) so that DMODEL and D match.

Does this help?

[mik2718]

Define closely, this means defining an error function.

A common one would be sum(DMODEL-D)^2 the smaller this is the better your match.

(It may be the case that you don't like this error function, maybe you want max{|DMODEL-D|} or something.
This will affect how hard it will be to find the optimum)

So now you have a function E(A,B,C) that you want to minimize. You could try a canned minimization routine, some can be found in the Gnu Scientific Library (gsl). There may be others available. This is a whole numerical area.

For the case where E(A,B,C) = sum(DMODEL-D)^2 its possible that the optimum can be found analytically - find dE/dA, dE/dB and dE/dC and calculate when they are zero (I haven't checked to see how easy this is)

If you want to roll your own optimiser here is an idea.

A,B,C space is a 3-dimensional space and we want to find a point in it where E is minimised. If we take two points in each if the A,B,C axes then this defines a cube in the space. Divide the cube into two in each direction giving 27 points and calculate E at each point.

case 1 : the middle point of the cube has the lowest E value. In this case we want to shrink the cube so as to find a more accurate optimum point.

case 2 : The E values along one side of the cube are smaller than the others. In this case we want to shift the cube until we get into a case 1 situation.

repeat until you have zeroed in on the optimum A,B,C

[Gaminic]

Don't think it's possible to solve analytically. The objective function is highly non-linear.