# CARPool - Clusters - First year project

# Research Meetings
## 12/04/2023
### Paper things
1 - $\chi^2$

I have removed the degree of freedom in the denominator for this, the trend is similar, but now, of course there are higher values.

Adding the TNG errorbars to the covariance
$$
\Sigma = C_{\rm CARPoolGP} + C_{\rm TNG}
$$
Where
$$
C_{ii, \rm TNG} = \sigma_i^2/N_i\delta_{ii}
$$
and $\sigma_i^2$ is the variance in the TNG bin, and $N_i$ is the number of halos in that bin.


We find that the columns remain linearly dependent now, but the eigenvalues are still very small after 8 points blowing up the determinant.

**Parameter dependencies**
*It is a bit unnerving to me that in some cases there is a strong discrepancy between Trho and P. The most prominent example are the P slope of the low IMF case, which isn't reflected in T and rho, and the big dip in rho for the low ASN1 case, which isn't reflected in P. Any idea how come? I'm thinking either due to imperfections of the emulations or some funny result of the averaging of complex distributions of these quantities. We recently discussed how to define the mean pressure in the shell; can you please remind me how the averaging is done for T & rho? Minor point: maybe better to exclude r>Rvir from the plots (do you physically trust that regime?)?*


Ok, I have changed these a little bit... First, I was training CARPoolGP and finding that the errors in the emulations were very high. I think this was because the values of the data were very small (eg 10^{-22} in cgs or 10^6K). When I train on the Log of the values I get results with smaller error bars, and seem more trustworthy. Especially look at IMF slope, which now has a density profile matching that of the Pressure as we would expect. Also interesting is the temperature of the ASN1, which seems to have had a reduction in the derivative, matching more what we find in the actual ASN1 simulations. I have also checked that the $P \sim nk_BT$. It seems like I am always a factor of 2 off in this calculation. But the trends seem to match quite well.


*Is the lower metallicity of the high ASN1 case low enough to bring the winds into the more efficient regime?*
Yes, somewhat. So for the highest ASN1, we are above the turning point in the Vogelsberger plot, but not quite at the fully efficient mode



*You write that the higher density as well as the lower temperature in the higher ASN1 case are because the winds pushed more gas into the CGM. This sounds perfectly reasonable, but perhaps it's also that the lower BHs of this case have not cause as much heating, and have not pushed gas completely out of the halo as much as in the lower ASN1 case? If you agree that's a possibility, I think it's worth presenting both (not necessarily figuring out in detail which dominates).*

I do agree, and I think in light of these new emulations + the simulation results this is what we are seeing.

*Is the emulation done in linear or log space? Am I reading the plots correctly that the averaging here (for the thick lines) is in linear space?*
The emulation is now done in log space. The simulations mean was calculated as the mean in linear space, but plotted in log space

*"The inner halo is the largest variance region and could maybe be why we see some convergence between low and high ASN in emulations" -- I'm not following the argument here.*
I am trying to say that the scatter in the temperature profiles is highest at the inner 10% of the halos, so maybe the emulator is also highest variance here (it is)

*The errors on the thick line in particular for T must be quite large given the scatter. I'd say it's hard to say if it is consistent with the emulation or not without adding error bars.*
This is what i was trying to say (thank you!)

*Is n_e on these plots the same quantity as rho in the emulation plots? Or is the former taking into account the ionization fraction, which the latter doesn't? Also (and depending on the answer to the previous question), let's discuss the treatment of star-forming gas in terms of its electron fraction.*
I am computing this as a volume weighted average... Which I dont think is right now that I think about it.
$$
n_{e, i} = \dfrac{\sum_k\dfrac{X_{e,k}X_H\rho_k}{m_p} dV_k}{V_i}
$$
where $X_e$ is the internal energy, $X_H$ is the hydrogen fraction (0.76). I am currently ignoring all star forming cells in all calculations.


This ones trippy




*Could the stellar fraction play a role in the discrepancy between the higher rho but similar f_b with higher ASN2? Also, it could be interesting to use the simulations to examine f_b.
*


## 11/30/2023
### Paper things
* Lets talk about Daisuke's comments, what do you think of the updated figures, and the added (very brief appendix)
* Still some minor edits to do in the first few sections
* Still need to make conclusion
* Ana marias results...
* Bens section
* Interpretations - I think this maybe helps?

I am getting confused again with IMF
### Chi2
First I compute chi2 with 15 points in the Y-M relation:
- We find that there are imaginary eignevalues in the covariance matrix which is due to roundoff errors. So this matrix has linearly dependent columns.

## 11/21/2023

### Paper things!
* Structure, discussion vs results, rename discussion?
* Individual comments
* What do we call predictive uncertainty?
* "The benefits of both the base samples, which can explore individual sample variances in a wide volume of parameter locations and the surrogate samples, which can provide an accurate estimation of the sample variance at a small number of parameter locations, are connected to one another through matrix $M_{ij}$."
* Histogram plot
* 
* " A comparison of this type then requires some level of smoothness along individual parameter directions and should be considered a very stringent test on CARPoolGP."
* $h^{-1}$ in Fisher
* Symmetrizing the errors in Fisher
* Marginalizing
* 
* 

### IMFslope confusion
Metal sheilding? Why is the temperature increasing? Also the baryon fraction seems weird? for redder IMFslope, we expect higher metalicity (which we find). and with higher metalicity we expect smaller mass loading, so

Look at how black hole parameters effect temperature. (increasing black home masses increase temperatures)

## 11/15/2023 - Ben and Shy
This just wasnt possible before!

Why do some of these parameter be better constrained than others??
Check out will coulton and Bens paper on noisy derivatives. Could be that the noise on derivatives is making the constraints better

Upper and lower bound of fisher matrices? Look into sonme software for immplementing ths?

How to estimate fisher matrices from simulations

Can also look at marginalizing over the rest of the parameters. look into this! just do all 28 paramers.

Look at the eigenvectors of the covariance matrix, might have small eigenvalue - These might tell us which points to removee.

find the largest number of points that I can with believeable results, do we get the same results as witht he three points.

Include the $\sigma_q^2$ in the chi^2

$C_{emu} + I\sigma(M)/\sqrt{N}$

Include this, start from 3 and include more low mass buns, see how the chi2 changes.

Convinced that theres a 2 sigma error (from shy) - statistical fluctuation

## 11/8/2023
### Constraints
the last time we saw that delta log y is 0.15, so the $(\Delta \log Y)^2$ should be of order $0.01$, but we forgot to divide by 2! so $\Delta \log Y =0.075$ and $(\Delta \log Y)^2 \sim 0.001$. This is what I am finding in my covariance matrices. I need to remake this plot with shared x and shared y... For some reason its not cooperating with me.


### fbar on ASN
Current TNG emulation


Error



correlation matrix


Just using the diagonals, the $\chi^2\sim 17$ (or $\chi^2_{dof}\sim 1.2$)

I cant figure out how to incorporate the diagonals because the determinant of the covariance matrix is so small. This makes the inverse $C^{-1}$ massive.

### Comparing fbar with ASN runs


### Ana marias plots
What does this meaaaaaaaan? is this just that the cosmic variance was not the main cause of the small importance on high masses?

### Interpretations (see the windows)

| Parameter | BH | metalicity | fbar | Y |
| -------- | -------- | -------- | --- | --- |
| ASN1 high | low | low | high | high |
| ASN2 high | low | low | NaN | high |
| IMF high | high | high | low | low |

## 11/3/2023
Fisher results from updated covariance matrix. I feel pretty good abojtu this. The derivatives are the same as when computed from finite differencing as well, so i am confident in the derivatives. Heres the full process.

I compute the derivatives with jax (or finite differencing)
$$
\dfrac{\partial \log Y}{\partial f_i}
$$

Where
$$
f_i = \begin{cases}
\dfrac{\theta_i - l_i}{l_i-u_i} + 1, \ \ \ \ \ \ \ \ \ \ \ \ \ \ &LogFlag_i=0 \\
\dfrac{\log(\theta_i) - \log(l_i)}{\log(l_i)-\log(u_i)} + 1, \ \ \ &LogFlag_i=1
\end{cases}
$$

Where l and u are lower and upper bounds, and log flag is determined from the parameter information.
so
$$
\dfrac{\partial f_i}{\partial \theta_i} = \begin{cases}
\dfrac{1}{l_i-u_i} &LogFlag_i=0\\
\dfrac{1}{\theta_i(\log(l_i)-\log(u_i))} \ \ \ &LogFlag_i=1
\end{cases}
$$

So from this with the chain rule
$$
\dfrac{\partial \log Y}{\partial \log\theta_i} =\dfrac{\partial \log Y}{\partial f_i}\dfrac{\partial f_i}{\partial \log\theta_i} = \theta_i\dfrac{\partial \log Y}{\partial f_i}\dfrac{\partial f_i}{\partial \theta_i} = \dfrac{\partial \log Y}{\partial f_i}\begin{cases}
\dfrac{\theta}{l_i-u_i} &LogFlag_i=0\\
\dfrac{1}{(\log(l_i)-\log(u_i))} \ \ \ &LogFlag_i=1
\end{cases}
$$

Then I compute the inverse Covariance matrix with diagonal elements from Shivam for three mass bins (13-13.5, 13.5-14, 14-14.5)
$$
C_{ii}^{-1} = \left(\dfrac{0.5\Delta Y}{Y}\right)^{-2}
$$
Where Log Y is the value of Y computed at the center of the bin from CARPoolGP. Then I compute the Fisher which has a diagonal covariance with three values
$$
F_{ab} = \sum_{i=1}^{3} \dfrac{\partial \log Y_{i}}{\partial \log\theta_a}C^{-1}_{ii}\dfrac{\partial \log Y_{i}}{\partial \log\theta_b}
$$

For the supernova and the AGN parameters I add a prior with diagonal elements $$
\sigma_{\log(\theta_i)}^2 = 1
$$
and then invert the matrix. I use this to compute to generate the contours below



### IMF slope
I have been trying too understand the trends better. For IMFSlope we see an increase in the SZ signal when we decrease IMF slope.




### Profiles?
Idk, I thought maybe looking at actual density profiles, temperature profiles, Ysz profiles might help better understand this? train CARPoolGP on a bunch of bins and see how different regions change.


## 10/27/2023



## 10/20/2023
Parameter consistency

Somehow the variance reduction gets worse with batch??

Comparing predictions sampled uniformly across grid with sample at the fiducial.

## 10/17/2023


## 10/9/2023

I messed up a little. The change to selecting the correct halos did have an influence

with slightly different parameter dependencies


I started looking at ffisher forecasting using the covariance from wadekar:

But I cant figure out how to get to the 1sigma plot that he had

## 9/21/2023 - Chat with Ben and Shy

Fbar results: this is converged

Adding the CV set halos

For Y:

The correlation matrix of the TNG predictions, shows correlation between first few data points
With deep kernel:


Which has different apparent parameter dependencies!

## Todo list:
* Compare all 3 batches so far in error reduction/emulations in general
* Show unncertainty on rainbow plots with errorbars on two most extreme cases (lowest vs highest)
* Compare the results of f_bar with halos fromt he 1P-set
* Put CPGP on rusty/popeye
* run at larger learning rate to compare
* run for 10^6 iterations and then look at the parameters at the end of this
* Make a rainbow plot to show how adding simulations effects the emulator for f_bar
* Look at noise as a function of Mass
* Compare the VW plot to the good TNG plot... Why do they look different?
* Look at changing individual parameters, predicting something and then seeing how important that parameter is.
* Look at changing M, Increasing this would likely degrade the results!
* Write Simulations and results section
* Reach out to nick and emily for covariance matrices etc.

## 9/5/2023
### F_{bar}
Last time we talked about seeing what these values were converging too. It looks like the loss function is not converging, and the emulator predictions arent converging either, I get something like this.


But I got and found that by switching optimizers to ADAM, and using an exponent in the scales (what i was doing originally), I converge on a solution much quicker. I actually do converge in these scenarios, and these are the results.

Which dont look great/look quite biased

Then i figured that if the CV set was at fiducial, it should match TNG, so I computed the fbar for this set and found this.

So still looks biased (though there are box size effects). I then figured, well what if I add in these data points as well? meaning I have ~640 halos from the Zooms and ~200 halos at the fiducial, with no correlations between the two.

Which is clearly better. So Im not sure what to think. We have a biased result from the zooms,

The parameter dependencies look very similar though

## 8/30/2023
### $Y_{SZ}$

Y_500
This is what we would expect from the Pop + paper. I find a slope of 1.68 for the high mass (compared to the 1.71 found in the paper)
Now the emulator gives the following result, which matches the data quite well.

Then, looking at the ratio for each parameter -

Y_200 we can do the same thing


### $f_bar$

## 8/16/2023

For the case with $V=W$
$V\neq W$

## 8/10/2023 - Meeting with Ben

**The good**
Ok, We have some results from the end of batch 1 and 2, and into batch 3!

This is our predictions from batch 1->2 and our results
So we dont perfectly match but are quite close.

Here are the predictions on Y_sz currently and what we expect in ~10 days

As well as the other parameter dependencies that we are finding:

Or, as a function of mass, we can see how each parameter changes Y

**The Bad**
Things seem to be working universally well for the Y_SZ which is great, but for other quantities there are some convergence issues in the GP.
For example:

Which is a good prediction, but with pretty high variance.

But I am having trouble getting the parameters to converge:

**The ugly**

How do we include multiple quantities into the covariance matrix? We take the distance wrt parameters in the space for each of the kernels. but to do this for the quantities, we would have to have some notion of distance in the quantity parameter space. Im not sure what the best way to do that is?

We have used the paramter $\Delta p$ from your notes to denote the differences in functional form between surrogates and base simulations, but

## 8/9/2023
Y results

predictions for stage 2.5


and these are the parameter fits with the data (blue), and surrogates (green)

I think there is something wrong here though with the data. I am only finding ~300 points in my data set. When I looked into it, it looked like there were alot of files where there are no star particles or BH particles.

## 8/2/2023
Results from Batch 2!

We find that sigma reduction is ~1.5x


I am finding similar results fro other quantities as well

Temperature

Central Entropy

We can also compare our predictions to the reality. By this i mean... Take the data from the first batch used to predict the next set of points. Then train a network and compute predictive variance. Use that trained network to predict the what the quantity should be at the next set of points. Then retrain the network and compute predictive variance. finally, use the actual data we have from the second batch to train a network and compute predictive variance.


## 7/26/2023 - Ben and Shy
### Quick recap
CARPool regression has been shown to work on a 1D case to reduce variance on estimates, Active learning approach has been effective at performing a better sampling of the parameter space when applied to a quantity of interest. So there are three things which we should discuss.

1. **What we expect** from using active learning with this zoom-in suite
2. **What we are finding so far** from the first batch of simulations that has finished
3. **What statistics best represent the improvement** and show our ability to sample the space. It sounded like Ben had some thoughts on this.

### What we expect
We can perform active learning given a toy model to try and predict what we expect from CARPool. The model itsself is trained on the 2048 sobol sequence where I use the largest halo in each simulations to compute the Y_sz and build an emulator to be the true "model." We do the active learning in 4 steps with the following layout:
1. Sobol sequence over all 29 parameters
2. 128 active learning HR sims and their 128 surrogates
3. 128 active learning HR sims and their 128 surrogate
4. 128 active learning HR sims and their 128 surrogate

At each step, we test the emulator on ~1000 points on a sobol sequence in parameter space that were evaluated using our model We find the error on these points and show them below. So here we are comparing a bunch of things.

In red we are comparing how the active learning reduces variance wrt iteration number. It looks like it is quickly converging after ~3 iterations.

We also compare active learning to vanilla carpool regression (in grey) and see that while the peak of the distribution is near the same, vanilla carpool has an extended tail to higher variances.

Finally we compare active learnings final step (dark red), CARPool standard (grey), and just a standard no carpool process (black).


But this was all using $Y_{SZ}$ as our quantity to actively learn. So the next question is them can we still benefit from our active learning on a different quantity like central entropy. To do this we train a CARPool network that has the parameter locations determined by the $Y_{SZ}$ training, but quantities that are evaluated using a model for entropy which was trained on the SB28 sets largest halos. The question then is: *do we find a stronger variance reduction in entropy using parameter locations trained on $Y_{SZ}$ than we would if we used a standard sobol sequence, or with no carpool at all*

The results are as follows: If we train on $Y_{SZ}$ and test on entropy, we do slightly worse than if we had just trained on entropy in terms of variance reduction. We also do slightly better than if we never used CARPool. It does appear that the standard CARPool performs very well for S though, and out performs the "trained on Y" case.

### What we find
#### Emulating on SB28 set
We have so far completed ~250 simulations. ~125 are HR simulations and the remaining are surrogates which have been distributed to 64 islands in parameter space. From these simulations we have computed $Y_{SZ}$ and have tried to emulate the value for halos that are in the SB28 set. When we look at the residuals of the predictions with the truth we find that the distribution is biased negative with a positive skew.


Much of the negative biasing is coming from the low mass halos that are in the SB28 set which our network cant predict well. But still we are underpredicting the value of the smallest halos? even in the other mass regimes we are under predicting.


One potential problem here is that we are using a zoom-in with modes from a 200 Mpc^3 box to predict on a 25 Mpc box. This may have some effect. So we further test the emulation by doing a jackknife sampling on the zoom-ins that we have, where we train the emulator with one zoom left out and predict that zoom. We do this for each of the HR zooms in the set and find that the mean goes to zero, but the tail is still extended and skewed positive.

#### Next active learning step
We use the data that we recieve from the first 256 simulations and perform an active learning step to predict where to place the next 128 simulations and their surrogates. The results are interesting, in that there are a few parameters that strongly favor their edges, while many other simulations are similar to uniformly drawn. We do find a significant favoring for low mass halos.

While this seems strange, especially for the variable wind spec momentum, we find that the network is predicting a much lower variance when these data points are added to the covariance matrix. We find a significant reduction in all predictions. There are still moderately high errors and some NaN predictions in the lowest masses where the network is extrapolating below the range of its training data.

**What statistics best represent the improvement**
I had heard that Ben has recently begun developing a statistic that describes the sampling in a high dimensional parameter space, I want to discuss this and see if it is applicable to this work!

**Next Steps**
Paper time!

## 7/19/2023



## 7/12/2023

### Timing problem
Rusty
Timestep issues - 95, 2, *207, 187, 239*

Popeye
Timestep issues - 238, 10, 242, 102, 93, 110, 169, *153, 148, 190, 4, 51, 16, 162*

Theta_R

### Voronoi mesh Problem
138, *211*

### Timings of current sims

### Timing tests on fiducial

### CARPool
I have cleaned up the code for this and now its as easy as...



Some variance reduction from active learning

How about now using a new sobol sequence that the system has never seen?

Now looking at its prediction on S


So what parameters is it actually choosing?

Which surrogates is it choosing?

Compared to a regular sobol sequence

Still todo - Train this/find the parameters that would come from entropy instead of Y. I am currently running the script. It takes ~10 hours to do the paramter-space search, but should be done soon.
#### debug
216, 168, 254, 138, 158, 250, 212, 19, 10, 93, 102
138 - has the voronoi problem. Just deleted ics, trying again.

27 - Restart with 2 nodes
2 - timestep problem, rerun with 2x timestep
95 - also timestep problem
59 - ???
## 7/7/2023
**Funky business - Popeye edition**
High z zooms -> this is because we restarted (for me 46) simulations which had bad lagrange regions.

zoom84 - Hydro canceled almost instantly with this error ```TERMINATE: ******!!!!!****** Code termination on task=6, function byten_hash(), file src/restart.c, line 2116: read error```

zoom 126 - There is no Hydro, it looks like it went straight to DMO. This is the error 
zoom 212, 250, 158, 168, 216 - Similar there is no Hydro, but ther is this massive table in the outfile followed byu this error message before cancelling. 

zoom 165, 197, 177 - Accidentally left stop files in these runs, so they just stopped.

zooms 93, 19 - It looks like something happend when I restarted the run?

zoom, 77, 103, 60, 44, 138 - Metallicity problem

zoom 10 - BH accretion problem? 

In general -  Is this right?

## 7/5/2023
all at fiducial
* 10^13, 10^14, 10^14.5 m cluster on different nodes, but make sure theyre the same. A few clusters.
* build a grid for this! cascadelake is going to be similar to icelake.
* 10^14 2 rome nodes and 10^14.5 2 rome nodes

WRITING! get on it.

Whats going on with zoom_111???

TERMINATE: ******!!!!!****** Code termination on task=3, function do_stellar_evolution(), file src/GFM/stellar_evolution_main.c, line 421: GFM_STELLAR_EVOLUTION: negative metal mass element returned: ID=
1086204892 P[iPart].Mass=0.00657737 sed->metal_mass_released[k_elem]=-1.59999e-06 initial_mass=0.00657737 k_elem=0 age_of_star_in_Gyr=0 dtime_in_Gyr=0.00324526 metallicity=0.179821 initial_metals=0.427
717 0.392462 0.0258001 0.00424771 0.0920791 0.0330225 0.00883146 0.00674583 0.00584146 0.0032527 metal_mass_released=-1.59999e-06 0.000293713 4.3554e-05 3.81243e-06 6.14696e-05 2.05627e-05 4.53125e-06 4.
69117e-06 3.56355e-06 2.14573e-06

Sims that need metallicity restart:
27
75
2
82
15
95
22
59
31

166
143
223

High Z! 
155 - 2e15 cluster! 
251 - 6e14 cluster
171 - 6e14 cluster
235 - 2e15 cluster

Stange path in 111? 
This is a fluke, I set a path wrong for zoom_111, but caught it and stopped the run and reset it. The output is showing this when I stopped the run. No other zooms have this, and the run itsself is not effected. (I searched "production" in the files as the strange path said "zoom_75_off_production")


## 6/28/2023
### Progress
**Timings**
for a $10^{14}M_\odot$ halo
![](https://hackmd.io/_uploads/ByKSyjKOh.png)

for a $10^{13}M_\odot$ halo
![](https://hackmd.io/_uploads/SJp_1sFdn.png)

But it seems that the halos that I zoomed in on are differnt?

They are the same parent boxes, the masses are identical between the zoomed halos, but the halos are different. 
![](https://hackmd.io/_uploads/H14YBituh.png)
![](https://hackmd.io/_uploads/SJN3rjKu2.png)


**Temperature**
So theres still this extended region of hot gas,
![](https://hackmd.io/_uploads/SkS4dsFOn.png)
But when we look at the dark matter
![](https://hackmd.io/_uploads/ryv0aiKdh.png)
outside of the red circle, the contamination rises to above 10%
So maybe this is an artifact? also, the parameters for this run were kind of nonsense... Like i had omega m=0.19, and messed up the omega_lambda. idk, I feel like this may just be a bad case of the whoopsies. 

for the 13 solar mass halo, we have contamination that looks like this
![](https://hackmd.io/_uploads/rJQ7HnKOn.png)

The gas seems reasonable
![](https://hackmd.io/_uploads/rkf5rnKu2.png)


following the darkmatter
![](https://hackmd.io/_uploads/Hkx3HhFO3.png)

and the 13M on icelake
![](https://hackmd.io/_uploads/Skv1PnFO2.png)
![](https://hackmd.io/_uploads/SJLjL2Kuh.png)
![](https://hackmd.io/_uploads/HJd28nt_n.png)

**Simulation Scripts**
I wrote a code to divy up the tasks to different core types depending on the mass. I did this to write the first 128 simulations and 128 surrogates. we can adjust the mass thresholds for the different types of cores as we see fit. I also put the largest masses in a bin for 2nodes, but Im not sure how best to implement this in an automated way. any thoughts?

I am also running a test right now for a 1e14 Msun halo to see how well the scaling is on two nodes

**CARPool**
I normalized the parameters before doing training on the 29 set and am finding pretty good results. 

This is active learning matching what we will be following in our sims (4 batches)
![](https://hackmd.io/_uploads/rJro8aKOn.png)

The variance is really reduced with just one instance of active learning and then decreases still after that. 
![](https://hackmd.io/_uploads/SJUKU6t_h.png)

Now training on the log of interior entropy I am finding that there is a moderate benefit, but that this is not nearly as pronounced as the Y scenario. 
![](https://hackmd.io/_uploads/SkSmWAtdn.png)

## 6/23/2023
### Progress

**Simulation Scripts** - The code now has all 28 parameters plus a mass parameter, the omega_lambda bug was patched up, and the linking to code in the directories now links and then runs the code from the link. I am running a test right now with all 29 parameters at the fiducial to see if this works. 

**CARPool** - Found a bunch of bugs in the 28 parameter version. Fixed them up and find that we can get a pretty good fit with reduced variance using the active learning method. 
![](https://hackmd.io/_uploads/HJpUUbmd2.png)
![](https://hackmd.io/_uploads/SkpqLWmdn.gif)

Theres still some issues here though. It is extremely sensitive to the initial guess of kernel parameters. So there was a bunch of tuning. I think this could be problematic? I was thinking, that maybe this could be because of the different parameter ranges. I think it would be good to regularize the parameters with some invertable function so that they are all on the same scale. I saw on tinygp website, dan foreman-mackey does this in an example with a small neural network, maybe this could be done? 

I have also started writing the code seriously and generally and have much of it built. I think with this it will be easier to debug. 

I havent had time to debug and look more into the cross predictions, but thats on the list for today! 


## 6/21/2023
We're getting close to the day where we press play... I think. 
### Progress
**Simulation scripts** - On popeye, Rusty, can be used for Icelake, skylake, cascadelake, and rome. Iam in the process of running tests, they are mostly all still running. Though theres a problem on popeye for the hydro runs right now, I am getting an error that the parameter IMFSlope was not defined. 

Theres a couple other things too. Time limits - The script stopped running after 1 week. Is there any way to by pass that? The. hydro sim went for about 6 of the days and died on snap32 on rome (~18000 core hours...) This also seems like alot of core hours. I think there is still something happening though where high res or gas particles are far from the center. This is the temperature map from the simulation. 
![](https://hackmd.io/_uploads/Sk9Qswgd3.png)
Which looks pretty extended. The extent was : HEXTENT: [0.10855746, 0.13159636, 0.11618054] Which i thought was reasonable
![](https://hackmd.io/_uploads/Bk_fnvxu2.png)
I think the time makes sense given the rpevious runs, 
![](https://hackmd.io/_uploads/B1pv2vldh.png)
But just wanted to make sense it was sensible. 

**Contamination** - After chatting a bit with Daisuke and Daniel at the conference, they were interested in the level of contamination out to certain radii to see what science uses could be worked with right away. 

I made this plot for the test halo I simulated (so this is at snap - 32). I use all of the DM particles here, not just the ones indexed from the subfind catalog. So as we go out to high radius we expect low res particles to take over. The splashback radius I read was about 0.9-1.5 R200 typically, so I put the band there, and the shock radius, I saw somewhere was ~3 R200, but was not super sure about this. It looks like after this the contamination goes beyondy 10%. This halo is about 0.9 Mpc, so that gives us ~3 Mpc before too much contamination. 
![](https://hackmd.io/_uploads/r19p3Ped2.png)

I also wasnt sure if this was the best way to look at it and thought maybe density plots would be better? 
![](https://hackmd.io/_uploads/rJuRpwxd3.png)
So here Im showing the density of just the high res DM particles, gas, low res, and then I put all three into an array and compute the density with all of them together to find the total. 


**CARPool** - New active learning method. I have implemented this following what we have discussed, where I go through an array of test points and add them one by one to the covariance matrix. I then use this new covariance matrix to predict at a range of test points, and sum up the variance over all of the test points. The added points that have produce the smallest sum of variance are used to generate new simulations. 

Heres an example on the 1P set. I am using just 3 active learning steps to mimic what we will likely be doing (256 simulations per step)
![](https://hackmd.io/_uploads/HJgtO_edn.gif)

I have also tried to test on the 28 parameter set and am finding not such great results. Something weird is happening that I haent seen in other cases. Im not sure why the variance in some places is 0 but in others is not. Its a bit confusing. But I think this all points to me going through and making this code really general, easy to debug, and just be generally usable. 
![](https://hackmd.io/_uploads/H1CmB_l_2.gif)

Finally, I have looked at the active learning (on 6 param) trained on Y and applied to T
![](https://hackmd.io/_uploads/Bk-nQtlO2.png)

It seems to do OK, though there is not as much of a benefit as we got from the Y clearly. 



### Chat with Ben

Ok there are a few things I want to discuss today,**N-dimensional tests, the case for the 10% islands, Active Learning in a toy model, Active Learning with the 1P set**

#### N-Dimensional tests
We want to test how well CARPoolGP extrapolates to Ndimesnions, and can it be used on something like the 1P set. 
To do this we 
1. Find the $Y_{SZ}$ from the largest halo in each of the LH simulations. 
2. Train a simple gaussian process to learn this function $Y_{SZ}(\vec{\theta}$). This is now our Model which we want to predict. 
3. We now draw data from this model and add noise (which is ~ the scatter in the Y-M relation)
4. Our goal is to recover the model given our draws
5. We do this with CARPoolGP (correlations between noise in islands)
6. And we do this with a vanilla GP (no correlations)


![](https://hackmd.io/_uploads/HJcfp40r2.png)

When I do this, I get very similar results that they both perform pretty well
![](https://hackmd.io/_uploads/B1d4seBrh.png)

I then iterate on this. I look at how well carpool does as a function of number of islands, and consider 20 experiments at each island. I look at the ratio of the $\chi^2$ for each to find when CARPool does best. The solid black line here is mean, but it clearly looks like the median would be lower. We ultimately get a factor of 2 times better in $\chi^2$ when using optimal number of islands

![](https://hackmd.io/_uploads/Hy21AlHSn.png)

In this plot, the horizontal line seperates where CARPool does better (below) and where No CARPool is better (above). The x-axis is the number of islands used. Each of the thin grey lines is a different random seed, and the black line is the average over random seeds. The "cup" then shows that using too few or too many islands not optimal! For a fixed size of 1024 simulations, where 512 are used as High Res simulations, and 512 are as surrogates, it is optimal to split the 512 simulations into ~50 islands. So there are ~10 simulations per island. This might change for different clustering algorightms as well. 

#### Active learning in a toy model
Now we we want to ask the quesion - can we improve our ability to emulate everywhere if we simulate at more "important" points in the space. 

![](https://hackmd.io/_uploads/BJ8ixbBB2.png)
Just to orient our selves. The red are sims with unique initial conditions, the blue are surrogates that live on islands with matched ICs to the nearest red points. This is a fit with only 24 simulations. But, say we had access to 240 simulations, and we knew this was the initial fit. How best should we place those remaining simulations?

Heres what we do: 
1. Generate a grid of points in the span of the parameter space, randomly drawn, lets say 100 points. 
2. Compute the variance at these points using the current kernel trained on available data. 
3. Find the N points with the highest variance
4. Find the Nearest surrogate islands to each of these points. 
5. Simulate both new data at each of the N points and each of the associated islands with correlated noise. 
6. Re train the network

![](https://hackmd.io/_uploads/S1V_-WrBn.png)
That peak around $\theta=9$ has high variance so the automation decides to add 8 points around there which in turn bolsters the island. Already things look better, We can do this again and again, untill we reach our desired 240 simulations. Heres a fun animation showing it. (I think the seed is different but otherwise the same concept)

![](https://hackmd.io/_uploads/BkILsGaBh.gif)

#### Active learning on the 1P set. 
Now we combine both of these thoughts. Lets say we can run a total of 1024 simulations. And we can do so in batches of 64 simulations. We take the model from before to be the truth for $Y_{SZ}(\vec{\theta})$. Then we generate 64 simulations on a random sobol sequence. We also create 50 islands with surrogates (10% of the number of surrogates we will end up with). These islands are done using an aglomorative clustering method. We create and train the GP on these 128 data points (64 high res and 64 surrogates). Then we compute the variance at points across this space using a sobol sequence, find the highest variance points, add 32 new points, add 32 new surrogates, and iterate. Heres how we end up doing. 

![](https://hackmd.io/_uploads/HykPsfarh.gif)


### 5/19/2023
#### Recap
Last time we talked about the differences in the subfind used on SIMBA vs TNG. This was what resulted in the results where SIMBA had apparently very cool halos. We also talked a bit about the Ndimensional test results on CARPool where we found that CARPool was marginally better, but there were some issues with the addition of noise in this setup, and it did not include mass as a training parameter. 

#### SIMBA ASTRID TNG Comparisons
I have updated the code to work on SIMBA. Here I get the ids for the largest halo, and then remove the particles that have a duplicate and are more than 4R200 away from the center of the halo. This is still a bit sketchy, as it only works for the largest halo. The rest of the halos seem to be cut off at high radius. That is ok though because I am just investigating the largest halo. 

![](https://hackmd.io/_uploads/Hk6f_eBSn.png)


Now looking at different values, The mass weighted temperature:
![](https://hackmd.io/_uploads/HkcIOlSHn.png)
Y_sz
![](https://hackmd.io/_uploads/Hy5ideHS2.png)

For extended SIMBA: I remove the first 6 parameters to keep the axes fixed and still see changes. 
![](https://hackmd.io/_uploads/B18QFgrB2.png)
![](https://hackmd.io/_uploads/HJKDtxHrn.png)

For extended TNG
![](https://hackmd.io/_uploads/rJL19lHBh.png)
![](https://hackmd.io/_uploads/H1l0YlBHh.png)

It seems like SIMBA parameters have a larger influence in the properties of the halo, but this is not by a massive amoung. 

#### N-D tests
I have fixed up the ND tests a bit. I have
* Added mass as a parameter
* Changed the noise in the model to be $\sigma^2=0.3$ from looking at the scatter in the Y-M plot
* Used an agglomerative clustering method instead of my homebrewed nearest neighbors

When I do this, I get very similar results that they both perform pretty well
![](https://hackmd.io/_uploads/B1d4seBrh.png)

I then iterate on this. I look at how well carpool does as a function of number of islands, and consider 20 experiments at each island. I look at the ratio of the $\chi^2$ for each to find when CARPool does best. 

![](https://hackmd.io/_uploads/Hy21AlHSn.png)

In this plot, the horizontal line seperates where CARPool does better (below) and where No CARPool is better (above). The x-axis is the number of islands used. Each of the thin grey lines is a different random seed, and the black line is the average over random seeds. The "cup" then shows that using too few or too many islands not optimal! For a fixed size of 1024 simulations, where 512 are used as High Res simulations, and 512 are as surrogates, it is optimal to split the 512 simulations into ~50 islands. So there are ~10 simulations per island. This might change for different clustering algorightms as well. 

#### Adaptive learning 1D
Ok, this is really cool... Adaptive learning. I made a toy model which using carpool gives me the following fit. 
![](https://hackmd.io/_uploads/BJ8ixbBB2.png)
Just to orient our selves. The red are sims with unique initial conditions, the blue are surrogates that live on islands with matched ICs to the nearest red points. This is a fit with only 24 simulations. But, say we had access to 240 simulations, and we knew this was the initial fit. How best should we place those remaining simulations?

Well lets assume that the regions with the highest variance need the most work. So we predict a bunch of points across this parameter space and find the 8 (or 10 or 100 etc) points with the highest variance. We would then want to run simulations there, and, run surrogate simulations at the nearest island. 
![](https://hackmd.io/_uploads/S1V_-WrBn.png)
That peak around $\theta=9$ has high variance so the automation decides to add 8 points around there which in turn bolsters the island. Already things look better, We can do this again and again, untill we reach our desired 240 simulations. Heres a fun animation showing it. (I think the seed is different but otherwise the same concept)

![](https://hackmd.io/_uploads/Hk5M1-rS2.gif)






### 5/12/2023
Its been so long! I have some updates this week, finally.
#### SIMBA ASTRID TNG Comparisons
I have compared the 1P set across simba astrid and TNG, though I am convinced there is something wrong here. I will show the comparison plots and then why I think they are wrong. 

I extract the largest halo for each of the 1P sets and compute the Temperature and $Y_{SZ}$. I do this for TNG, Astrid and SIMBA. I only look at the original 1P set here:


![](https://hackmd.io/_uploads/B1cZVToVn.png)
![](https://hackmd.io/_uploads/rkLfNasV3.png)

I would expect SIMBA to have a higher temperature consistintly when compared to TNG
I went digging into the halos a bit and found that there is something wrong with the way I am extracting halos I think. If I plot the TNG halo for the first subfind group, I find: 

![](https://hackmd.io/_uploads/SyrqEToV3.png)

But when I do the same thing for SIMBA I get:

![](https://hackmd.io/_uploads/rkyCVpsV3.png)

So something is off. The orange point is the group center, and the green point is an R_200 away. 

I can redo this with only including points that are withing 2Rvir or 3Rvir, but is there an easier way to get the particles associated with the halo?

There is a similar problem when I look at astrid:
![](https://hackmd.io/_uploads/H1pHrTsN2.png)
The halo seems much more extended than the radius implies. Any thoughts here?

#### CARPool with LH model
I applied CARPool to a model where I trained on the largest halo in the LH set. So for each LH run, I extracted the largest halos temperature and Y_Sz. Ill show T results. I estimate the variance on the measurement by looking at the variance among the 10 largest halos at ~1e13. The results of using a regular GP on the set I find something that looks close to the 1P set but biased cooler. 
![](https://hackmd.io/_uploads/rkoUoaj4h.png)


Then I run carpool, where I have one fine grain grid of 512 points and coars sequence of 256 points. Each sample has a surrogate

I compare this to a standard GP with 1024 points

I test them both on 512 new parameter points and find residuals are pretty consistent
![](https://hackmd.io/_uploads/SJZoT6o42.png)

Similar errors too
![](https://hackmd.io/_uploads/SJ2-ATsN2.png)

When I check on the "1P modelled" set i find pretty consistent results as well. 
![](https://hackmd.io/_uploads/S127Rpi4n.png)

For $Y_{SZ}$ its a little different, with CARPool doing better
![](https://hackmd.io/_uploads/SkSF0Tj43.png)


### 4/21/2023
### Updates:
Not too much to report, I have been looking into the N-D tests for CARPool. Heres what Im finding. 

This is the data i test on. This comes from computing $f(\theta)$ on the parameters from the 1P set. 
![](https://i.imgur.com/qTblOzZ.png)

I use noise that correspondes to the variance among temperatures at the fiducial which is ~0.5 here. 

I then train both carpool (512 points 512 surrogates) and a no regular gaussian process (1024 points).

![](https://i.imgur.com/3h7wxrh.png)
Find marginal benefit. Computing the score of the regression
$$
R^2 = 1- u/v
$$
where $u=\sum(G(\theta, x) - data)^2$ and $v=\sum(data - \langle data \rangle)^2$

$R_{CP} = 0.995$
$R_{GP} = 0.984$

Now applying this to a wide range of test points:
![](https://i.imgur.com/OUdhDo0.png)
$R_{CP} = 0.972$
$R_{GP} = 0.969$
![](https://hackmd.io/_uploads/HJcfp40r2.png)

## 4/18/2023
### Last week recap
Last week we had two meetings, one with ben and one without. 

Here we discussed our project with ben, testing carpool, checking profiles before and after excising the SFR region, running simba vs TNG, running a test with multiple resolutions. 

**Little Ben**: We want to have him write a density fit to $10^{12}M_\odot$ halos, so we should get him to read a bit on this. 

**Testing CARPool**: ND tests, active learning tests, and how to place surrogates. 
We also talked about placing surrogates with big ben and the different ways to split the space. We thought of a few ways to make islands of surrogates. We could either have a coarser sobol sequence whicch we find the nearest neighbors for in parameter space, but this might lead to uneven islands, or we could do an eigenvector split where we divide the space along principle axes to the desired number of cuts. We should check both but not dwell too much on this. It is more important to get the ND test and the active learning example set up!

**SFR particles/AGN?**: We have some characteristic hump in our plots of halo pressure as a function of mass at $10^{12}$. We think this is when the AGN turns on and causes some of this scatter. but we want to clear this up and look at pressure profiles before and after excising SFR particles
 
**SIMBA vs TNG**: SIMBA's AGN is more efficient, so we may want to use SIMBA instead of TNG galaxy formation. We should compare the T, Y, S, P etc between the two. 

**Resolution tests**: We want to run a bunch of halos at a few points in parameter space to see if there is a strong correlation in resolution. 


### This weeks progress
#### Resolution test
I have four parameter points that I get from a sobol sequence
![](https://i.imgur.com/fGJ4jOP.png)
And at each point I sample four times
![](https://i.imgur.com/JzNVrZf.png)
I do this at two different resolutions, where I am only changing the highest resolution level (so one is at $512^3$, one is at $1024^3$)
![](https://i.imgur.com/wpWNe6A.png)

I think this next part is maybe not good... I also used mass as a parameter in the sobol sequence: 
![](https://i.imgur.com/0EYh85v.png)
But these masses are associated with the parameters, so the first om0, s8, agn etc combination is run with the different seeds and different resolutions but always for a  $3.7335\times 10^{14}M_\odot$ halo... Maybe I should have kept all the halos the same mass, or done just a totally random selection across all halos, or matched them for each parameter point?

We will see i guess if this is not good.

Looking at one of the zooms:
![](https://i.imgur.com/K2zeXcl.png)

some are still running, so not quite done yet. Heres what I have so far.
![](https://i.imgur.com/88vklKO.png)

#### ND-carpool tests:
I generate polynomial fits to the mean temperature over mean temperatures of 20 halos at $\sim 10^{12} M_\odot$ in the 1P set. 

![](https://i.imgur.com/rv95sPW.png)

I then use each of these fits as a function to generate a model:
$$
f(\vec{\theta}) = f(\theta_1)f(\theta_2)..f(\theta_N) + \epsilon
$$
This is the model we want to fit with CARPoolGP. I generate test data using a sobol sequence across the parameter points, I make 512 points for the sobol sequence using the model above, and add noise to the final data. 

I then generate a coarse grid sobol sequence of 128 points with a different initial seed to act as surrogate island locations. I find the nearest island for each high resolution point and create a simulation at the island with matched noise realization. I run a gaussian process with scale parameters for each of the dimensions and find the best fit parameters for my covariance matrix. 

We want to test this against a process that has more data with no correlations. So I generate a new sobol sequence with 1024 points (more than double). I run a gaussian process on this data set as well.

I then create a new realization of points to test on. I test both processes on 512 new points. and compare the histogram of residuals: 

![](https://i.imgur.com/8hJ7zgh.png)


Also comparing the variances on the predictions, the CARPool method seems to be much much smaller
![](https://i.imgur.com/ep2Kkzi.png)

Looking at $\chi^2$ there is ~4x decrease when using carpool. And this is with fewer samples!
![](https://i.imgur.com/nWI0e2O.png)

#### Pressure differences
![](https://i.imgur.com/Pehtr2q.png)

## 3/31/2023
### Last week recap
Last week I showed the newest zoom, which was a 5e13 mass halo, but had an error in the feedback parameters. I also had found that there was an error in the way I was dealing with the edgees of the periodic box. We talked about the new set of camels simulations which would be run with a larger box and therefore contain the 1e13 mass halos, so for our project, going to higher halo mass is important. 

### This weeks progress

#### Newest Zoom
I fixed up the halo function to wrap around the box, and then I ran another test for a more massive halo, 6e14. I ran this with the fiducial cosmology and values. 

![](https://i.imgur.com/shv0VMm.png)
![](https://i.imgur.com/WMyFjnP.png)

Looking at the values of the Y (same trend for P, S, and T)
![](https://i.imgur.com/eqcX7y0.png)

Which has a trend close to $M\propto Y%{3/5}$
![](https://i.imgur.com/k12ZQcV.png)


So I think things are working here?

#### Finding parameters
I added the fix suggested to remove the influence of star forming particles in the calculation of the profiles. This reduced alot of the discrepancy between the pressure scatter as a function of mass. But theres still a pretty significant spread (3 orders of mag). Is this ok? I guess for what were looking at, only the highest mass's matter, so maybe we just limit our attention to this?
![](https://i.imgur.com/qTbsXBp.png)

I have a new set of plots with these changes

#### Carpool
I applied the new kernels to the 1P set, but got mostly the same results, that there is not much effect. I still think this is because theres just not enough data points. So we need a way to try this out on a larger data set. 


## 3/31/2023
### Last week recap
We found that there was some bugs in my zooms. I had been running my simulations for the parent with a zoom-like configuration, I was not using tracers, had some inconsistencies between the double precision, and had no compression on the output. But there are a bunch of minor things that we have been missing which I wanted to discuss today

### This weeks progress

#### New Zoom
I fixed up the configs etc and ran a new zoom, I did it for a slightly larger mass
HYDRO run time rome 128, mass 5e13, 119960.472548 (33 hrs)
DMO run time rome 128, mass 5e13, 29789.639496 (8 hours)
![](https://i.imgur.com/e8pLp9Q.png)

there are 25 contaminant points, 34 for the DMO run
![](https://i.imgur.com/jZQ3fUC.png)

One of the parameters AAGN2 - ```RadioFeedbackReiorientationFactor``` I am setting to 0.7 or something very small when compared to values that I have seen used (20). But the LH set is uniformly really small. Should I be multiplying 20 by this value?

Are there any other tests we should do before running the mini test? 

4 points in parameter space - 10 sims at each point with different ICs. 10 matching sims at each point with lower resolution. 

Remaining todo here:
- [ ] Fix the halo function to allow for surrogates. 
- [x] Fix the halo function to prevent halos on the edges of the box. 
- [ ] Choose parameters
- [ ] Run mini test, run bigger test, run total

#### Finding parameters
I computed for matching halos across the 1P-28 set the T, P, Y, S. So now we can look at which parameters have the largest effect on the halos, (notebook mean_quantities2)... These dont include the zoom ins yet though... I need to do these as well, but we have old plots which kind of show which are the most important [old plot](#11/2/2022)
![](https://i.imgur.com/FqQpkZi.png) 

(a bit different but still maybe helpful)

Parameter with little influence
 ![](https://i.imgur.com/WtAU6VZ.png)
 
 
lets open the notebook and just look at the parameters and which have the most effect. 
<!-- ![](https://i.imgur.com/GSkYMZS.png)
![](https://i.imgur.com/k5YtkZs.png)
![](https://i.imgur.com/x7Wo7J2.png)
![](https://i.imgur.com/JY4u8eV.png)

 Parameters that I have found with seemingly significant influence on P. 
![](https://i.imgur.com/Siv6tQm.png)
![](https://i.imgur.com/FJ6CA1c.png)
![](https://i.imgur.com/bbM14P9.png)
![](https://i.imgur.com/ORSzChK.png) -->


#### CARPool
I implemented the new kernel into the 1P set, am not finding good results... I think again this is just that we have too few data points to be doing this regression on. I think I am going to try to use the fits to generate a new function and do ND tests now though, Im not sure. I almost think that we wont know if it works untill we run the whole simulations suite