input: Write a Python function to optimize a model using Particle Swarm Optimization (PSO) with given hyperparameters.
  assert isinstance(optimize(model, X, y), type(model(**{'bias': 0, 'weights': [0]*X.shape[1]})))
  assert hasattr(optimize(model, X, y), 'predict')
output: def optimize(model, X, y, num_particles=30, max_iter=100, w=0.5, c1=1.5, c2=1.5):
      import numpy as np
  
      def fitness(particle):
          m = model(**particle)
          m.fit(X, y)
          predictions = m.predict(X)
          return np.mean((predictions - y) ** 2)
  
      # Initialize particles
      particles = [{'bias': np.random.randn(), 'weights': np.random.randn(X.shape[1])} for _ in range(num_particles)]
      velocities = [{'bias': np.random.randn(), 'weights': np.random.randn(X.shape[1])} for _ in range(num_particles)]
      personal_best = particles.copy()
      personal_best_scores = np.array([fitness(p) for p in particles])
      global_best = personal_best[np.argmin(personal_best_scores)]
  
      for iteration in range(max_iter):
          for i, particle in enumerate(particles):
              for key in particle:
                  velocities[i][key] = (w * velocities[i][key] 
                                        + c1 * np.random.rand() * (personal_best[i][key] - particle[key]) 
                                        + c2 * np.random.rand() * (global_best[key] - particle[key]))
                  particle[key] += velocities[i][key]
  
              current_score = fitness(particle)
              if current_score < personal_best_scores[i]:
                  personal_best[i] = particle.copy()
                  personal_best_scores[i] = current_score
          
          global_best = personal_best[np.argmin(personal_best_scores)]
  
      best_model = model(**global_best)
      best_model.fit(X, y)
      return best_model