Skip to content
Snippets Groups Projects
Commit 9ae06572 authored by lykos98's avatar lykos98
Browse files

added working kstar density

parent 8507be87
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
coherency_tests.ipynb
+ 33
− 23
View file @ 9ae06572
%% Cell type:code id:d1f5f17d-9949-4707-b553-e03c15bc59ee tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
import dadac
```
%% Cell type:code id:aeb8a19b-be6d-48d8-86e0-479b4d6d99cf tags:
``` python
x = np.fromfile("bb/ordered_data.npy", dtype = np.float64)
x = x.reshape((x.shape[0]//5, 5))
x.shape
```
%% Output
(1842842, 5)
(921421, 5)
%% Cell type:code id:550c24bc-062c-4e9a-83ab-5fb51aa44dc4 tags:
``` python
data = dadac.Data(x)
```
%% Output
You are running in a notebook maybe the timing output will break, but everything should be fine
%% Cell type:code id:83e48b9f-1739-4b4c-aeea-87fb732c02e2 tags:
``` python
data.compute_distances(299)
data.compute_id_2NN()
#data.id = 4
data.id = 4
data.compute_density_kstarNN()
```
%% Output
Building the KDtree v2:
Total time: 1.741s
Total time: 1.607s
knn search:
Total time: 18.676s
Total time: 8.044s
ID estimation:
ID value: 3.920865
Total time: 0.391s
ID value: 3.862658
Total time: 0.180s
Density and k* estimation:
Total time: 5.015s
Total time: 2.085s
%% Cell type:code id:314b516a-02ec-412a-bc7e-1de8f790e429 tags:
``` python
den_gt = data.log_den
den_comp = np.fromfile("bb/ordered_density.npy", np.float64)
print(data.id)
```
%% Output
3.920865231328582
4
%% Cell type:code id:92225324-e798-4388-a022-b6bb8906768c tags:
``` python
np.average(np.abs(den_comp - den_gt))
```
%% Output
0.0001499206205206319
8.791685787264807e-08
%% Cell type:code id:67a8bfdf-f1c1-421a-ab06-0ed9dc57b1b5 tags:
``` python
den_gt
```
%% Output
array([-3.64725383, -4.53602916, -4.30573383, ..., 7.18148125,
-6.25416517, -3.54944958])
array([-3.06640128, -3.86980622, -3.66751667, ..., -0.64454652,
-3.35926525, -7.39476774])
%% Cell type:code id:a7093d5f-d40b-43d9-a4d0-1f70638d9bbd tags:
``` python
den_comp
```
%% Output
array([-3.64714384, -4.53592253, -4.30561972, ..., 7.18179655,
-6.25407934, -3.54934192])
array([-3.06640124, -3.86980629, -3.66751671, ..., -0.64454651,
-3.35926533, -7.39476776])
%% Cell type:code id:f02c64d5-18b8-4b88-a756-e00fa9295e64 tags:
``` python
a = np.argmax(np.abs(den_comp - den_gt))
print(den_comp[a], den_gt[a])
print(den_comp[a], den_gt[a], "max diff", np.abs(den_comp[a] - den_gt[a]))
a
```
%% Output
-0.5825421810150146 -1.2760780561214355
-8.463785171508789 -8.463785648321064 max diff 4.768122749965187e-07
1789174
674003
%% Cell type:code id:10dcc6ee-af0f-412d-8be9-d7500be6cd08 tags:
%% Cell type:code id:cd441b5f-f7d8-476c-9b8f-2c1a7f92380f tags:
``` python
np.abs(den_comp - den_gt)[a]
data.kstar[a]
```
%% Output
0.6935358751064209
14
%% Cell type:code id:19b1c9e6-bf00-44f9-822f-f13412263ce5 tags:
%% Cell type:code id:1338cf31-e49d-47da-bbb4-1225db9aefd7 tags:
``` python
data.kstar[a]
plt.plot(np.abs(den_comp - den_gt), ".-")
```
%% Output
130
[<matplotlib.lines.Line2D at 0x7485f7bad5d0>]
%% Cell type:code id:cd441b5f-f7d8-476c-9b8f-2c1a7f92380f tags:
%% Cell type:code id:a309b4e5-0ae4-4b55-a971-6fb05320c49b tags:
``` python
```
......
src/common/common.c
+ 4
− 0
View file @ 9ae06572
......@@ -20,6 +20,7 @@ void get_context(global_context_t* ctx)
ctx -> n_halo_points_send = NULL;
ctx -> halo_datapoints = NULL;
ctx -> local_datapoints = NULL;
ctx -> __recv_heap_buffers = NULL;
}
void free_context(global_context_t* ctx)
......@@ -39,14 +40,17 @@ void free_context(global_context_t* ctx)
{
for(int i = 0; i < ctx -> world_size; ++i)
{
/*
for(int j = 0; j < ctx -> n_halo_points_recv[i]; ++j)
{
FREE_NOT_NULL(ctx -> halo_datapoints[i][j].ngbh.data);
}
*/
FREE_NOT_NULL(ctx -> halo_datapoints[i]);
}
}
FREE_NOT_NULL(ctx -> halo_datapoints);
FREE_NOT_NULL(ctx -> __recv_heap_buffers);
if(ctx -> idx_halo_points_recv)
{
......
src/common/common.h
+ 4
− 0
View file @ 9ae06572
......@@ -29,9 +29,12 @@ typedef struct datapoint_info_t {
#define float_t double
#endif
#define MY_TRUE 1
#define MY_FALSE 0
#define CHECK_ALLOCATION(x) if(!x){printf("[!!!] %d rank encountered failed allocation at line %s ", ctx -> mpi_rank, __LINE__ ); exit(1);};
#define DB_PRINT(...) printf(__VA_ARGS__)
#ifdef NDEBUG
#undef DB_PRINT(...)
......@@ -126,6 +129,7 @@ struct global_context_t
int* rank_n_points;
char processor_mame[MPI_MAX_PROCESSOR_NAME];
MPI_Comm mpi_communicator;
heap_node* __recv_heap_buffers;
};
struct pointset_t
......
src/tree/tree.c
+ 14
− 5
View file @ 9ae06572
......@@ -2366,12 +2366,20 @@ void find_foreign_nodes(global_context_t* ctx, datapoint_info_t* dp, datapoint_i
{
for(int j = 0; j < n_heap_to_recv[i]; ++j)
{
/*
foreign_dp[i][j].array_idx = array_indexes_to_request[i][j];
init_heap(&(foreign_dp[i][j].ngbh));
allocate_heap(&(foreign_dp[i][j].ngbh), k);
foreign_dp[i][j].ngbh.N = k;
foreign_dp[i][j].ngbh.count = k;
memcpy(foreign_dp[i][j].ngbh.data, heap_buffer_to_recv + k * (j + rdispls[i]), k * sizeof(heap_node));
*/
foreign_dp[i][j].array_idx = array_indexes_to_request[i][j];
//init_heap(&(foreign_dp[i][j].ngbh));
foreign_dp[i][j].ngbh.N = k;
foreign_dp[i][j].ngbh.count = k;
foreign_dp[i][j].ngbh.data = heap_buffer_to_recv + k * (j + rdispls[i]);
if(foreign_dp[i][j].ngbh.data[0].array_idx != array_indexes_to_request[i][j])
{
......@@ -2402,7 +2410,8 @@ void find_foreign_nodes(global_context_t* ctx, datapoint_info_t* dp, datapoint_i
free(count_to_request);
free(capacities);
free(heap_buffer_to_recv);
/* free(heap_buffer_to_recv); this needs to be preserved*/
ctx -> __recv_heap_buffers = heap_buffer_to_recv;
free(heap_buffer_to_send);
free(idx_buffer_to_send);
free(idx_buffer_to_recv);
......@@ -2701,7 +2710,7 @@ void simulate_master_read_and_scatter(int dims, size_t n, global_context_t *ctx)
//ctx -> n_points = 48*5*2000;
ctx->n_points = ctx->n_points / ctx->dims;
// ctx->n_points = (ctx->n_points * 0.5) / 10;
ctx->n_points = (ctx->n_points * 0.1) / 10;
// ctx -> n_points = ctx -> world_size * 1000;
//ctx -> n_points = 10000000 * ctx -> world_size;
......@@ -2817,7 +2826,7 @@ void simulate_master_read_and_scatter(int dims, size_t n, global_context_t *ctx)
elapsed_time = TIME_STOP;
//id = 3.920865231328582;
//id = 4.008350298212649;
//id = 4.;
id = 4.;
LOG_WRITE("ID estimate", elapsed_time)
MPI_DB_PRINT("ID %lf \n",id);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment