sentient ai

# Define parameters
N = 100  # number of particles
D = 2  # number of dimensions
steps = 100  # number of time steps
step_size = 0.1  # fixed distance particles move at each step
radius = 1.0  # radius of the circle particles start in

# Generate initial positions
theta = 2 * np.pi * np.random.rand(N)
r = radius * np.sqrt(np.random.rand(N))
x = r * np.cos(theta)
y = r * np.sin(theta)
pos = np.column_stack((x, y))

# Generate random directions for each particle
angles = 2 * np.pi * np.random.rand(N)
directions = np.column_stack((np.cos(angles), np.sin(angles)))

# Generate random walk trajectories
traj = np.zeros((N, D, steps + 1))
traj[:, :, 0] = pos
for i in range(steps):
    pos += step_size * directions
    traj[:, :, i + 1] = pos

# Plot trajectories
fig, ax = plt.subplots()
ax.set_xlim(-radius, radius)
ax.set_ylim(-radius, radius)
for i in range(N):
    ax.plot(traj[i, 0, :], traj[i, 1, :])
plt.show()

# Bing chat [1 prompt]

import numpy as np
import matplotlib.pyplot as plt

N = 100  # number of particles
T = 100  # number of timesteps
R = 5  # radius of circle

# initialize positions
theta = np.random.uniform(0, 2 * np.pi, N)
r = R * np.sqrt(np.random.uniform(0, 1, N))
x = r * np.cos(theta)
y = r * np.sin(theta)
positions = np.stack([x, y], axis=1)

# plot initial positions
plt.scatter(x, y)
plt.xlim(-R, R)
plt.ylim(-R, R)

# random walk
for t in range(T):
    # generate random directions
    theta = np.random.uniform(0, 2 * np.pi, N)
    # update positions
    positions += np.stack([np.cos(theta), np.sin(theta)], axis=1)
    # plot positions
    plt.scatter(positions[:, 0], positions[:, 1])
    plt.xlim(-R - np.sqrt(T), R + np.sqrt(T))
    plt.ylim(-R - np.sqrt(T), R + np.sqrt(T))

plt.show()

N = 1000
Nt = 100

phis = np.random.uniform(0, 2 * np.pi, (N, Nt))
rs = np.sqrt(np.random.random(N))
thetas = np.random.uniform(0, 2 * np.pi, N)

f = lambda proj: np.tile(rs * proj(thetas), (Nt, 1)).T + np.cumsum(proj(phis), axis=1)
xs = f(np.cos)
ys = f(np.sin)

plt.plot(xs.T, ys.T, "k", alpha=0.1, lw=0.5)
plt.gca().set_aspect("equal")

def perlin_noise_2d(shape, res):
    def surflet(grid_x, grid_y, x, y):
        dist_x, dist_y = abs(x - grid_x), abs(y - grid_y)
        poly_x = 1 - 6 * dist_x**5 + 15 * dist_x**4 - 10 * dist_x**3
        poly_y = 1 - 6 * dist_y**5 + 15 * dist_y**4 - 10 * dist_y**3
        grad = grid[(grid_x % shape[0]), (grid_y % shape[1])]
        dot_prod = dist_x * grad[0] + dist_y * grad[1]
        return dot_prod * poly_x * poly_y

    def interp(x, func):
        return lambda idx: func(idx) * (1 - x % 1) + func(idx + 1) * (x % 1)

    # Generate a grid of vectors
    grid = np.random.normal(size=(shape[0] + 1, shape[1] + 1, 2))
    # Generate a coordinate grid of shape res
    coords = np.mgrid[0 : shape[0] : res[0] * 1j, 0 : shape[1] : res[1] * 1j].astype(
        np.int32
    )
    # Generate an empty noise image
    noise = np.zeros(shape)

    for i in range(res[0]):
        for j in range(res[1]):
            # Calculate the coordinates of the 4 surrounding grid points
            x_min, x_max = coords[0, i, j], coords[0, i + 1, j]
            y_min, y_max = coords[1, i, j], coords[1, i, j + 1]
            # Calculate the contributions from the 4 surrounding surflets
            s = (
                surflet(x_min, y_min, *coords[:, i, j])
                + surflet(x_max, y_min, *coords[:, i + 1, j])
                + surflet(x_min, y_max, *coords[:, i, j + 1])
                + surflet(x_max, y_max, *coords[:, i + 1, j + 1])
            )
            # Interpolate between those values based on the location within the square
            noise[x_min:x_max, y_min:y_max] += interp(
                coords[0, i, j] / res[0], interp(coords[1, i, j] / res[1], s)
            )

    return noise


# Generate the noise image
noise = np.zeros((512, 512))
for i in range(4):
    noise += perlin_noise_2d((512, 512), res=(16 * 2**i, 16 * 2**i)) / 2**i

# Plot the noise image
plt.imshow(noise, cmap="gray")
plt.show()

def perlin(x, y, seed=0):
    # permutation table
    np.random.seed(seed)
    p = np.arange(256, dtype=int)
    np.random.shuffle(p)
    p = np.stack([p, p]).flatten()
    # coordinates of the top-left
    xi = x.astype(int)
    yi = y.astype(int)
    # internal coordinates
    xf = x - xi
    yf = y - yi
    # fade factors
    u = fade(xf)
    v = fade(yf)
    # noise components
    n00 = gradient(p[p[xi] + yi], xf, yf)
    n01 = gradient(p[p[xi] + yi + 1], xf, yf - 1)
    n11 = gradient(p[p[xi + 1] + yi + 1], xf - 1, yf - 1)
    n10 = gradient(p[p[xi + 1] + yi], xf - 1, yf)
    # combine noises
    x1 = lerp(n00, n10, u)
    x2 = lerp(n01, n11, u)  # FIX1: I was using n10 instead of n01
    return lerp(x1, x2, v)  # FIX2: I also had to reverse x1 and x2 here


def lerp(a, b, x):
    "linear interpolation"
    return a + x * (b - a)


def fade(t):
    "6t^5 - 15t^4 + 10t^3"
    return 6 * t**5 - 15 * t**4 + 10 * t**3


def gradient(h, x, y):
    "grad converts h to the right gradient vector and return the dot product with (x,y)"
    vectors = np.array([[0, 1], [0, -1], [1, 0], [-1, 0]])
    g = vectors[h % 4]
    return g[:, :, 0] * x + g[:, :, 1] * y


lin = np.linspace(0, 5, 512, endpoint=False)
x, y = np.meshgrid(lin, lin)  # FIX3: I had to swap x and y here

plt.imshow(perlin(x, y, seed=2), origin="upper")
plt.show()

def generate_perlin_noise_2d(shape, res):
    def f(t):
        return 6 * t**5 - 15 * t**4 + 10 * t**3

    delta = (res[0] / shape[0], res[1] / shape[1])
    d = (shape[0] // res[0], shape[1] // res[1])
    grid = np.mgrid[0 : res[0] : delta[0], 0 : res[1] : delta[1]].transpose(1, 2, 0) % 1
    # Gradients
    angles = 2 * np.pi * np.random.rand(res[0] + 1, res[1] + 1)
    gradients = np.dstack((np.cos(angles), np.sin(angles)))
    g00 = gradients[0:-1, 0:-1].repeat(d[0], 0).repeat(d[1], 1)
    g10 = gradients[1:, 0:-1].repeat(d[0], 0).repeat(d[1], 1)
    g01 = gradients[0:-1, 1:].repeat(d[0], 0).repeat(d[1], 1)
    g11 = gradients[1:, 1:].repeat(d[0], 0).repeat(d[1], 1)
    # Ramps
    n00 = np.sum(grid * g00, 2)
    n10 = np.sum(np.dstack((grid[:, :, 0] - 1, grid[:, :, 1])) * g10, 2)
    n01 = np.sum(np.dstack((grid[:, :, 0], grid[:, :, 1] - 1)) * g01, 2)
    n11 = np.sum(np.dstack((grid[:, :, 0] - 1, grid[:, :, 1] - 1)) * g11, 2)
    # Interpolation
    t = f(grid)
    n0 = n00 * (1 - t[:, :, 0]) + t[:, :, 0] * n10
    n1 = n01 * (1 - t[:, :, 0]) + t[:, :, 0] * n11
    return np.sqrt(2) * ((1 - t[:, :, 1]) * n0 + t[:, :, 1] * n1)


noise = generate_perlin_noise_2d((512, 512), (16, 16))

plt.imshow(noise, cmap="viridis")

#include "diffusion_solver.h"
#include <algorithm>
#include <cmath>
#include <iostream>

void plot(const std::vector<double> &y) {
  std::cout << "y_max = " << *std::max_element(y.begin(), y.end()) << std::endl;
  const auto y_max = 1.0;
  const auto ny = 20;
  const auto nx = 50;
  std::vector<std::string> plot_str;
  for (std::size_t i{0}; i < ny; ++i) {
    plot_str.push_back(std::string(nx, ' '));
  }
  for (std::size_t i{0}; i < y.size(); ++i) {
    const auto x_i = static_cast<int>(nx * i / y.size());
    const auto y_i = static_cast<int>(ny * y[i] / y_max);
    if (x_i < nx && y_i < ny && x_i >= 0 && y_i >= 0) {
      plot_str[ny - y_i - 1][x_i] = '.';
    }
  }
  for (std::size_t i{0}; i < ny; ++i) {
    std::cout << plot_str[i] << std::endl;
  }
}

int main(int argc, char **argv) {
  // Define the problem parameters
  const int num_points = 100;
  const double domain_length = 1.0;
  const double delta_x = domain_length / (num_points - 1);
  const double delta_t = 0.00005;

  std::cout << "CFL: " << 2.0 * delta_t / (delta_x * delta_x) << std::endl;

  // Create the diffusion solver
  DiffusionSolver solver(num_points, delta_x, delta_t);

  // Set the initial condition
  solver.set_initial_condition(
      [](double x) { return std::pow(x * (1 - x) / 0.25, 5); });

  // Set the boundary conditions
  solver.set_boundary_conditions([](double x) { return 0; },
                                 [](double x) { return 0; });

  // Advance the solution in time
  const int num_time_steps = atoi(argv[1]);
  for (int i = 0; i < num_time_steps; i++) {
    solver.advance();
  }

  // Print the final solution
  const auto &solution = solver.get_solution();
  plot(solution);

  return 0;
}

#include <functional>
#include <vector>

class DiffusionSolver {
public:
  DiffusionSolver(int num_points, double delta_x, double delta_t);
  void set_initial_condition(std::function<double(double)> func);
  void set_boundary_conditions(std::function<double(double)> left,
                               std::function<double(double)> right);
  void advance();
  const std::vector<double> &get_solution() const;

private:
  std::vector<double> m_solution;
  std::vector<double> m_previous;
  double m_delta_x;
  double m_delta_t;
  double m_coefficient;
  std::function<double(double)> m_left_bc;
  std::function<double(double)> m_right_bc;
};

#include "diffusion_solver.h"

DiffusionSolver::DiffusionSolver(int num_points, double delta_x, double delta_t)
    : m_solution(num_points), m_previous(num_points), m_delta_x(delta_x),
      m_delta_t(delta_t) {
  m_coefficient = delta_t / (delta_x * delta_x);
}

void DiffusionSolver::set_initial_condition(
    std::function<double(double)> func) {
  for (int i = 0; i < m_solution.size(); i++) {
    m_solution[i] = func(i * m_delta_x);
    m_previous[i] = m_solution[i];
  }
}

void DiffusionSolver::set_boundary_conditions(
    std::function<double(double)> left, std::function<double(double)> right) {
  m_left_bc = left;
  m_right_bc = right;
}

void DiffusionSolver::advance() {
  // Apply boundary conditions
  m_solution[0] = m_left_bc(0);
  m_solution[m_solution.size() - 1] =
      m_right_bc((m_solution.size() - 1) * m_delta_x);

  // Compute the next time step
  for (int i = 1; i < m_solution.size() - 1; i++) {
    m_solution[i] =
        m_previous[i] + m_coefficient * (m_previous[i - 1] - 2 * m_previous[i] +
                                         m_previous[i + 1]);
  }

  // Update the previous time step
  std::swap(m_solution, m_previous);
}

const std::vector<double> &DiffusionSolver::get_solution() const {
  return m_solution;
}

cmake_minimum_required(VERSION 3.10)
project(diffusion_solver)

set(CMAKE_CXX_STANDARD 11)

add_executable(diffusion_solver main.cpp diffusion_solver.cpp)

#ifndef DIFFUSIONSOLVER_H
#define DIFFUSIONSOLVER_H

#include <functional>
#include <vector>

class DiffusionSolver {
public:
  DiffusionSolver(int n_points, double dx, double dt);
  void set_initial_condition(std::function<double(double)> f);
  void set_boundary_conditions(std::function<double(double)> f);
  void advance();
  std::vector<double> get_solution() const;

private:
  int n_points;
  double dx;
  double dt;
  std::vector<double> u;
};

#endif

#include "diffusion_solver.h"

DiffusionSolver::DiffusionSolver(int n_points, double dx, double dt)
    : n_points(n_points), dx(dx), dt(dt), u(n_points) {}

void DiffusionSolver::set_initial_condition(std::function<double(double)> f) {
  for (int i = 0; i < n_points; ++i) {
    u[i] = f(i * dx);
  }
}

void DiffusionSolver::set_boundary_conditions(std::function<double(double)> f) {
  u[0] = f(0);
  u[n_points - 1] = f((n_points - 1) * dx);
}

void DiffusionSolver::advance() {
  std::vector<double> u_new(n_points);
  for (int i = 1; i < n_points - 1; ++i) {
    u_new[i] = u[i] + dt / (dx * dx) * (u[i + 1] - 2 * u[i] + u[i - 1]);
  }
  u = u_new;
}

std::vector<double> DiffusionSolver::get_solution() const { return u; }

cmake_minimum_required(VERSION 3.10)
project(diffusion_solver)

set(CMAKE_CXX_STANDARD 11)

add_executable(diffusion_solver main.cpp diffusion_solver.cpp)