CVXPY in Python

August 9, 2024

CVXPY in Python

Here’s a basic example of using $CVXPY$, a $Python$ library for convex optimization.

This example solves a simple linear programming problem:

Problem

Minimize the function $( c^T x )$ subject to the constraint $( Ax \leq b )$ and $( x \geq 0 )$, where:

$( c )$ is a vector of coefficients for the objective function.
$( A )$ is a matrix of coefficients for the inequality constraints.
$( b )$ is a vector representing the upper bounds for the constraints.

Sample Code

import cvxpy as cp
import numpy as np

# Define the data for the problem
c = np.array([1, 2])        # Coefficients for the objective function
A = np.array([[1, 1], [1, -1], [-1, 2]])  # Coefficients for the constraints
b = np.array([2, 1, 2])     # Right-hand side values for the constraints

# Define the optimization variables
x = cp.Variable(2)

# Define the objective function: minimize c^T x
objective = cp.Minimize(c @ x)

# Define the constraints: Ax <= b and x >= 0
constraints = [A @ x <= b, x >= 0]

# Formulate the problem
problem = cp.Problem(objective, constraints)

# Solve the problem
problem.solve()

# Print the results
print("Status:", problem.status)
print("Optimal value of the objective function:", problem.value)
print("Optimal values of the variables x:", x.value)

Explanation

Objective Function:
c @ x is the dot product of the vector c with the variable vector x. We aim to minimize this value.
Constraints:
A @ x <= b represents the inequality constraints, and x >= 0 ensures that the variables are non-negative.
Optimization:
problem.solve() solves the optimization problem, and the optimal solution is stored in x.value.

Output

When you run this code, it will output the status of the optimization (e.g., “optimal”), the optimal value of the objective function, and the optimal values of the decision variables $( x )$.

1
2
3

Status: optimal
Optimal value of the objective function: 3.698490338406144e-10
Optimal values of the variables x: [3.59660015e-10 5.09450927e-12]

This example is a basic introduction, but $CVXPY$ can handle more complex problems, including quadratic programming, mixed-integer programming, and other types of convex optimization problems.

Prophet in Python

August 8, 2024

Prophet in Python

$Prophet$ is a popular library developed by $Facebook$ for time series forecasting.

It’s particularly effective for data that has strong seasonal effects and multiple seasonality with daily observations.

Here’s a basic example of how to use $Prophet$ in $Python$:

Step-by-Step Example:

Install Prophet (if you haven’t already):
1
pip install prophet

Import Required Libraries:

1
2
3

import pandas as pd
from prophet import Prophet
import matplotlib.pyplot as plt

Load Your Data:
For this example, let’s create a simple time series dataset.

# Create a simple time series dataset
dates = pd.date_range(start='2022-01-01', periods=365)
data = pd.DataFrame({
    'ds': dates,
    'y': 100 + (dates.dayofyear - 183) ** 2 / 100 + np.random.randn(365) * 5
})

Initialize and Fit the Prophet Model:

# Initialize the Prophet model
model = Prophet()

# Fit the model to the data
model.fit(data)

Make Predictions:
You can make future predictions using the model by specifying the number of days into the future you want to forecast.

# Create a dataframe for future predictions
future = model.make_future_dataframe(periods=30)  # Predict 30 days into the future

# Predict future values
forecast = model.predict(future)

Visualize the Forecast:
$Prophet$ has a built-in plot function to visualize the forecasted data.
1
2
3
# Plot the forecast
model.plot(forecast)
plt.show()
Plot Components:
You can also plot the components (trend, weekly seasonality, yearly seasonality) of the forecast.
1
2
3
# Plot the forecast components
model.plot_components(forecast)
plt.show()

Full Code Example:

import pandas as pd
import numpy as np
from prophet import Prophet
import matplotlib.pyplot as plt

# Step 1: Create a simple time series dataset
dates = pd.date_range(start='2022-01-01', periods=365)
data = pd.DataFrame({
    'ds': dates,
    'y': 100 + (dates.dayofyear - 183) ** 2 / 100 + np.random.randn(365) * 5
})

# Step 2: Initialize the Prophet model
model = Prophet()

# Step 3: Fit the model to the data
model.fit(data)

# Step 4: Create a dataframe for future predictions
future = model.make_future_dataframe(periods=30)  # Predict 30 days into the future

# Step 5: Predict future values
forecast = model.predict(future)

# Step 6: Plot the forecast
model.plot(forecast)
plt.show()

# Step 7: Plot the forecast components
model.plot_components(forecast)
plt.show()

Explanation

ds: The column containing the dates.
y: The column containing the values to be forecasted.
fit: The method to train the model with your time series data.
make_future_dataframe: Prepares a dataframe to hold future predictions.
predict: Generates predictions for the given dates.
plot: Visualizes the forecast along with the observed data.
plot_components: Breaks down the forecast into its components (e.g., trend, weekly seasonality).

Result

Running this code will generate a plot of the time series data with the forecasted values and their uncertainty intervals, as well as a breakdown of the forecast components.

DEAP in Python

August 7, 2024

The $DEAP (Distributed Evolutionary Algorithms in Python)$ library is a powerful and flexible library used for implementing evolutionary algorithms.

It is widely used for optimization tasks, particularly in the fields of genetic algorithms, genetic programming, and other evolutionary strategies.

Steps to Use DEAP Library in Python

Installation
You can install $DEAP$ using pip:
1
pip install deap
Basic Structure of a $DEAP$ Program
A typical $DEAP$ program involves the following steps:
- Importing $DEAP$ modules
- Defining the problem (fitness function)
- Creating individuals and population
- Defining the evolutionary operators (selection, mutation, crossover)
- Setting up the algorithm flow
- Running the algorithm

Example: Simple Genetic Algorithm

Below is an example of a simple genetic algorithm using $DEAP$ to maximize the function f(x) = x^2, where x is an integer.

import random
from deap import base, creator, tools, algorithms

# Step 1: Define the problem (Fitness Function)
creator.create("FitnessMax", base.Fitness, weights=(1.0,))
creator.create("Individual", list, fitness=creator.FitnessMax)

# Step 2: Define the individual and population
toolbox = base.Toolbox()
toolbox.register("attr_int", random.randint, 0, 100)
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_int, 10)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)

# Step 3: Define the fitness function
def evalOneMax(individual):
   return sum(individual),  # Summing all values in the individual

toolbox.register("evaluate", evalOneMax)
toolbox.register("mate", tools.cxTwoPoint)
toolbox.register("mutate", tools.mutFlipBit, indpb=0.05)
toolbox.register("select", tools.selTournament, tournsize=3)

# Step 4: Set up the evolutionary algorithm
def main():
   population = toolbox.population(n=300)
   NGEN = 40
   CXPB = 0.7
   MUTPB = 0.2

   # Evaluate the entire population
   fitnesses = list(map(toolbox.evaluate, population))
   for ind, fit in zip(population, fitnesses):
       ind.fitness.values = fit

   # Step 5: Run the algorithm
   for gen in range(NGEN):
       offspring = toolbox.select(population, len(population))
       offspring = list(map(toolbox.clone, offspring))

       for child1, child2 in zip(offspring[::2], offspring[1::2]):
           if random.random() < CXPB:
               toolbox.mate(child1, child2)
               del child1.fitness.values
               del child2.fitness.values

       for mutant in offspring:
           if random.random() < MUTPB:
               toolbox.mutate(mutant)
               del mutant.fitness.values

       invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
       fitnesses = map(toolbox.evaluate, invalid_ind)
       for ind, fit in zip(invalid_ind, fitnesses):
           ind.fitness.values = fit

       population[:] = offspring

       fits = [ind.fitness.values[0] for ind in population]

       length = len(population)
       mean = sum(fits) / length
       sum2 = sum(x*x for x in fits)
       std = abs(sum2 / length - mean**2)**0.5

       print(f"Gen {gen}: Min {min(fits)} Max {max(fits)} Avg {mean} Std {std}")

   best_ind = tools.selBest(population, 1)[0]
   print("Best individual is:", best_ind, best_ind.fitness.values)

if __name__ == "__main__":
   main()

Explanation

Problem Definition: The fitness function (evalOneMax) is defined to maximize x^2.
Individual and Population: An individual is a list of integers, and the population is a list of individuals.
Evolutionary Operators:
- mate: Uses two-point crossover.
- mutate: Uses bit flip mutation with a small probability.
- select: Uses tournament selection.
Algorithm Flow: The algorithm runs for a specified number of generations (NGEN), performing selection, crossover, mutation, and evaluation in each generation.
Output: The algorithm prints the minimum, maximum, average, and standard deviation of fitness values at each generation and finally outputs the best individual found.

Output

Gen 0: Min 309.0 Max 787.0 Avg 562.7166666666667 Std 85.22716540060568
Gen 1: Min 358.0 Max 832.0 Avg 627.3066666666666 Std 76.52445767349342
Gen 2: Min 500.0 Max 862.0 Avg 680.6666666666666 Std 69.2917182801974
Gen 3: Min 521.0 Max 868.0 Avg 741.4666666666667 Std 60.18418027651102
Gen 4: Min 596.0 Max 924.0 Avg 781.9566666666667 Std 53.113226120136446
Gen 5: Min 624.0 Max 929.0 Avg 820.27 Std 50.70007001967588
Gen 6: Min 567.0 Max 941.0 Avg 849.4633333333334 Std 57.31010954059929
Gen 7: Min 681.0 Max 948.0 Avg 883.84 Std 42.49165878930357
Gen 8: Min 645.0 Max 971.0 Avg 904.7033333333334 Std 44.194969422119485
Gen 9: Min 744.0 Max 971.0 Avg 920.1233333333333 Std 34.01541006988307
Gen 10: Min 663.0 Max 976.0 Avg 930.0766666666667 Std 35.58638488086127
Gen 11: Min 671.0 Max 978.0 Avg 936.6533333333333 Std 38.082539597856474
Gen 12: Min 745.0 Max 981.0 Avg 947.5033333333333 Std 31.812628345080054
Gen 13: Min 854.0 Max 983.0 Avg 957.0966666666667 Std 25.950606715236244
Gen 14: Min 764.0 Max 987.0 Avg 960.02 Std 33.497755148667416
Gen 15: Min 691.0 Max 987.0 Avg 967.7833333333333 Std 32.976401494941406
Gen 16: Min 587.0 Max 988.0 Avg 966.6833333333333 Std 43.08553185880072
Gen 17: Min 784.0 Max 988.0 Avg 973.6033333333334 Std 27.84276067889405
Gen 18: Min 692.0 Max 988.0 Avg 974.3433333333334 Std 33.028151460366374
Gen 19: Min 690.0 Max 988.0 Avg 973.64 Std 39.178187809035656
Gen 20: Min 785.0 Max 988.0 Avg 980.99 Std 27.707578385703428
Gen 21: Min 789.0 Max 988.0 Avg 980.0 Std 31.1358314486693
Gen 22: Min 690.0 Max 988.0 Avg 978.8 Std 32.96391966984678
Gen 23: Min 794.0 Max 988.0 Avg 981.1233333333333 Std 26.293499619149838
Gen 24: Min 691.0 Max 988.0 Avg 981.0566666666666 Std 33.116463008937174
Gen 25: Min 788.0 Max 988.0 Avg 979.43 Std 31.13858967048373
Gen 26: Min 694.0 Max 988.0 Avg 981.4366666666666 Std 30.476209118300122
Gen 27: Min 788.0 Max 988.0 Avg 975.47 Std 38.43005464476852
Gen 28: Min 689.0 Max 988.0 Avg 977.1333333333333 Std 40.12283915953909
Gen 29: Min 789.0 Max 988.0 Avg 979.78 Std 32.672694001772946
Gen 30: Min 688.0 Max 988.0 Avg 978.4633333333334 Std 35.322164744281395
Gen 31: Min 697.0 Max 988.0 Avg 977.1433333333333 Std 35.65991384672285
Gen 32: Min 691.0 Max 988.0 Avg 977.4666666666667 Std 39.729530019313195
Gen 33: Min 689.0 Max 988.0 Avg 977.0733333333334 Std 35.092752654769896
Gen 34: Min 694.0 Max 988.0 Avg 979.7866666666666 Std 33.611919843347344
Gen 35: Min 689.0 Max 988.0 Avg 981.41 Std 29.468094045369302
Gen 36: Min 788.0 Max 988.0 Avg 979.11 Std 34.46947296763569
Gen 37: Min 694.0 Max 988.0 Avg 980.1366666666667 Std 31.0949833395823
Gen 38: Min 788.0 Max 988.0 Avg 981.07 Std 30.00364033468625
Gen 39: Min 691.0 Max 988.0 Avg 981.0266666666666 Std 32.25077707935549
Best individual is: [98, 100, 100, 97, 94, 99, 100, 100, 100, 100] (988.0,)

The output represents the results of an evolutionary algorithm running for $40$ generations.

Here’s a breakdown of what each part of the output means:

Generations (Gen 0, Gen 1, etc.):
- Each “Gen” line corresponds to a generation in the evolutionary process. The algorithm starts at Gen 0 and evolves through to Gen $39$.
Min, Max, Avg, Std:
- Min: The minimum fitness value in the population for that generation.
- Max: The maximum fitness value in the population for that generation.
- Avg: The average fitness value across the entire population for that generation.
- Std: The standard deviation of fitness values in the population, indicating how spread out the fitness values are.
Trends over Generations:
- Min: The minimum fitness value generally increases over time, suggesting that the worst-performing individuals are improving.
- Max: The maximum fitness value also increases, indicating that the best individuals are becoming fitter.
- Avg: The average fitness steadily rises, showing that the overall population is improving.
- Std: The standard deviation decreases over time, meaning the population is becoming more homogenous, with individuals having similar fitness values.
Final Generation (Gen 39):
- Min 691.0, Max 988.0, Avg 981.0266666666666, Std 32.25077707935549: By the final generation, the population has a high average fitness, with the best individual achieving a fitness of $988.0$.
Best Individual:
- The best individual found during the evolutionary process is [98, 100, 100, 97, 94, 99, 100, 100, 100, 100] with a fitness value of 988.0. This individual likely represents an optimal or near-optimal solution to the problem.

In summary, the algorithm successfully evolves a population over $40$ generations, improving the fitness of individuals, and converging towards an optimal solution, as indicated by the high average and maximum fitness values in the later generations.

Customizing DEAP

Custom Fitness Functions: You can create custom fitness functions based on the problem at hand.
Custom Selection, Crossover, and Mutation: $DEAP$ allows you to define and use custom operators.
Multi-Objective Optimization: $DEAP$ also supports multi-objective optimization, where the fitness function can return multiple values (e.g., minimizing both time and cost).

This example demonstrates how to set up and run a basic genetic algorithm using $DEAP$.

You can extend it by adjusting the problem, operators, and parameters to suit different optimization tasks.

Polars in Python

August 6, 2024

Polars in Python

$Polars$ is a powerful DataFrame library in Python, known for its speed and efficiency, especially with large datasets.

Here’s a useful sample demonstrating some key features of $Polars$, including reading data, data manipulation, and aggregations.

We’ll use a sample dataset to illustrate these functionalities.

Sample Data

Let’s assume we have a CSV file named data.csv with the following content:

id,name,age,salary,department
1,Alice,30,70000,HR
2,Bob,40,80000,Engineering
3,Charlie,25,65000,Marketing
4,Diana,35,72000,HR
5,Edward,50,90000,Engineering
6,Frank,45,85000,Marketing

Code Example

import polars as pl

# Read CSV file into a Polars DataFrame
df = pl.read_csv('data.csv')

# Display the DataFrame
print("Original DataFrame:")
print(df)

# Select specific columns
selected_df = df.select(['name', 'age', 'salary'])
print("\nSelected Columns (name, age, salary):")
print(selected_df)

# Filter rows where age > 30
filtered_df = df.filter(pl.col('age') > 30)
print("\nFiltered DataFrame (age > 30):")
print(filtered_df)

# Add a new column 'bonus' which is 10% of the salary
df = df.with_column((pl.col('salary') * 0.1).alias('bonus'))
print("\nDataFrame with Bonus Column:")
print(df)

# Group by 'department' and calculate the average salary and total bonus
grouped_df = df.groupby('department').agg([
    pl.col('salary').mean().alias('avg_salary'),
    pl.col('bonus').sum().alias('total_bonus')
])
print("\nGrouped by Department with Aggregations:")
print(grouped_df)

# Sort the DataFrame by 'salary' in descending order
sorted_df = df.sort('salary', reverse=True)
print("\nSorted DataFrame by Salary (Descending):")
print(sorted_df)

Explanation

Reading CSV:
- pl.read_csv('data.csv') reads the CSV file into a $Polars$ DataFrame.
Displaying the DataFrame:
- print(df) shows the original DataFrame.
Selecting Specific Columns:
- df.select(['name', 'age', 'salary']) selects only the name, age, and salary columns.
Filtering Rows:
- df.filter(pl.col('age') > 30) filters the DataFrame to include only rows where the age is greater than $30$.
Adding a New Column:
- df.with_column((pl.col('salary') * 0.1).alias('bonus')) adds a new column bonus which is $10$% of the salary.
Grouping and Aggregating:
- df.groupby('department').agg([pl.col('salary').mean().alias('avg_salary'), pl.col('bonus').sum().alias('total_bonus')]) groups the DataFrame by department and calculates the average salary and total bonus for each department.
Sorting:
- df.sort('salary', reverse=True) sorts the DataFrame by salary in descending order.

Output

The output will look something like this:

Original DataFrame:
shape: (6, 5)
┌─────┬─────────┬─────┬────────┬─────────────┐
│ id  │ name    │ age │ salary │ department  │
│ --- │ ---     │ --- │ ---    │ ---         │
│ i64 │ str     │ i64 │ i64    │ str         │
├─────┼─────────┼─────┼────────┼─────────────┤
│ 1   │ Alice   │ 30  │ 70000  │ HR          │
│ 2   │ Bob     │ 40  │ 80000  │ Engineering │
│ 3   │ Charlie │ 25  │ 65000  │ Marketing   │
│ 4   │ Diana   │ 35  │ 72000  │ HR          │
│ 5   │ Edward  │ 50  │ 90000  │ Engineering │
│ 6   │ Frank   │ 45  │ 85000  │ Marketing   │
└─────┴─────────┴─────┴────────┴─────────────┘

Selected Columns (name, age, salary):
shape: (6, 3)
┌─────────┬─────┬────────┐
│ name    │ age │ salary │
│ ---     │ --- │ ---    │
│ str     │ i64 │ i64    │
├─────────┼─────┼────────┤
│ Alice   │ 30  │ 70000  │
│ Bob     │ 40  │ 80000  │
│ Charlie │ 25  │ 65000  │
│ Diana   │ 35  │ 72000  │
│ Edward  │ 50  │ 90000  │
│ Frank   │ 45  │ 85000  │
└─────────┴─────┴────────┘

Filtered DataFrame (age > 30):
shape: (4, 5)
┌─────┬─────────┬─────┬────────┬─────────────┐
│ id  │ name    │ age │ salary │ department  │
│ --- │ ---     │ --- │ ---    │ ---         │
│ i64 │ str     │ i64 │ i64    │ str         │
├─────┼─────────┼─────┼────────┼─────────────┤
│ 2   │ Bob     │ 40  │ 80000  │ Engineering │
│ 4   │ Diana   │ 35  │ 72000  │ HR          │
│ 5   │ Edward  │ 50  │ 90000  │ Engineering │
│ 6   │ Frank   │ 45  │ 85000  │ Marketing   │
└─────┴─────────┴─────┴────────┴─────────────┘

DataFrame with Bonus Column:
shape: (6, 6)
┌─────┬─────────┬─────┬────────┬─────────────┬───────┐
│ id  │ name    │ age │ salary │ department  │ bonus │
│ --- │ ---     │ --- │ ---    │ ---         │ ---   │
│ i64 │ str     │ i64 │ i64    │ str         │ f64   │
├─────┼─────────┼─────┼────────┼─────────────┼───────┤
│ 1   │ Alice   │ 30  │ 70000  │ HR          │ 7000  │
│ 2   │ Bob     │ 40  │ 80000  │ Engineering │ 8000  │
│ 3   │ Charlie │ 25  │ 65000  │ Marketing   │ 6500  │
│ 4   │ Diana   │ 35  │ 72000  │ HR          │ 7200  │
│ 5   │ Edward  │ 50  │ 90000  │ Engineering │ 9000  │
│ 6   │ Frank   │ 45  │ 85000  │ Marketing   │ 8500  │
└─────┴─────────┴─────┴────────┴─────────────┴───────┘

Grouped by Department with Aggregations:
shape: (3, 3)
┌─────────────┬────────────┬────────────┐
│ department  │ avg_salary │ total_bonus│
│ ---         │ ---        │ ---        │
│ str         │ f64        │ f64        │
├─────────────┼────────────┼────────────┤
│ HR          │ 71000      │ 14200      │
│ Engineering │ 85000      │ 17000      │
│ Marketing   │ 75000      │ 15000      │
└─────────────┴────────────┴────────────┘

Sorted DataFrame by Salary (Descending):
shape: (6, 6)
┌─────┬─────────┬─────┬────────┬─────────────┬───────┐
│ id  │ name    │ age │ salary │ department  │ bonus │
│ --- │ ---     │ --- │ ---    │ ---         │ ---   │
│ i64 │ str     │ i64 │ i64    │ str         │ f64   │
├─────┼─────────┼─────┼────────┼─────────────┼───────┤
│ 5   │ Edward  │ 50  │ 90000  │ Engineering │ 9000  │
│ 6   │ Frank   │ 45  │ 85000  │ Marketing   │ 8500  │
│ 2   │ Bob     │ 40  │ 80000  │ Engineering │ 8000  │
│ 4   │ Diana   │ 35  │ 72000  │ HR          │ 7200  │
│ 

1   │ Alice   │ 30  │ 70000  │ HR          │ 7000  │
│ 3   │ Charlie │ 25  │ 65000  │ Marketing   │ 6500  │
└─────┴─────────┴─────┴────────┴─────────────┴───────┘

This example covers some fundamental operations with the $Polars$ library, making it a useful starting point for data manipulation and analysis tasks.

SciPy in Python

August 5, 2024

SciPy in Python

Here are a few useful SciPy sample codes showcasing different aspects of the library, including optimization, integration, interpolation, and signal processing.

1. Optimization: Minimizing a Function

This example demonstrates how to minimize a mathematical function using scipy.optimize.

import numpy as np
from scipy.optimize import minimize

# Define a function to minimize
def objective_function(x):
    return x**2 + 10*np.sin(x)

# Initial guess
x0 = 2.0

# Perform the minimization
result = minimize(objective_function, x0, method='BFGS')

print("Minimum value of the function:", result.fun)
print("Value of x that minimizes the function:", result.x)

[Output]

Minimum value of the function: 8.31558557947746
Value of x that minimizes the function: [3.83746713]

2. Integration: Calculating the Area Under a Curve

This example shows how to calculate the definite integral of a function using scipy.integrate.quad.

from scipy.integrate import quad

# Define a function to integrate
def integrand(x):
    return np.exp(-x**2)

# Compute the integral from 0 to infinity
result, error = quad(integrand, 0, np.inf)

print("Integral result:", result)

[Output]

Integral result: 0.8862269254527579

3. Interpolation: Interpolating Data Points

This example demonstrates how to interpolate data points using scipy.interpolate.interp1d.

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

# Sample data points
x = np.linspace(0, 10, 10)
y = np.sin(x)

# Create an interpolation function
f = interp1d(x, y, kind='cubic')

# Generate new x values and interpolate
x_new = np.linspace(0, 10, 100)
y_new = f(x_new)

# Plot the results
plt.plot(x, y, 'o', label='Data points')
plt.plot(x_new, y_new, '-', label='Cubic interpolation')
plt.legend()
plt.show()

[Output]

4. Signal Processing: Filtering a Signal

This example shows how to apply a low-pass filter to a noisy signal using scipy.signal.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, lfilter

# Generate a noisy signal
np.random.seed(0)
t = np.linspace(0, 1, 500, endpoint=False)
signal = np.sin(2 * np.pi * 7 * t) + 0.5 * np.random.randn(500)

# Design a low-pass filter
def butter_lowpass(cutoff, fs, order=5):
    nyquist = 0.5 * fs
    normal_cutoff = cutoff / nyquist
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

# Apply the filter to the signal
def lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y

# Parameters
cutoff = 3.5  # Desired cutoff frequency of the filter, Hz
fs = 50.0     # Sample rate, Hz

# Filter the signal
filtered_signal = lowpass_filter(signal, cutoff, fs)

# Plot the results
plt.plot(t, signal, label='Noisy signal')
plt.plot(t, filtered_signal, label='Filtered signal', linewidth=2)
plt.legend()
plt.show()

[Output]

These examples provide a glimpse of the powerful tools available in SciPy for scientific computing, including optimization, integration, interpolation, and signal processing.

PuLP in Python

August 4, 2024

PuLP in Python

PuLP is a Python library used for linear programming (LP) and mixed-integer linear programming (MILP).

It allows you to define and solve optimization problems where you want to minimize or maximize a linear objective function subject to linear constraints.

Installation

First, you need to install PuLP. You can do this using pip:

1	pip install pulp

Basic Usage of PuLP

Here’s a step-by-step guide to solving a simple linear programming problem with PuLP.

Example Problem:

Suppose you are a factory manager and you want to determine the optimal number of products $A$ and $B$ to produce to maximize your profit.
Each product requires a certain amount of resources (time, materials) and yields a certain profit.

Objective: Maximize profit.
Constraints:
- You have $100$ units of material.
- You have $80$ hours of labor.
- Product $A$ requires $4$ units of material and $2$ hours of labor.
- Product $B$ requires $3$ units of material and $5$ hours of labor.
Profit:
- Product $A$ gives $20 profit per unit.
- Product $B$ gives $25 profit per unit.

Step-by-Step Implementation

import pulp

# Step 1: Define the problem
# Create a maximization problem
problem = pulp.LpProblem("Maximize_Profit", pulp.LpMaximize)

# Step 2: Define the decision variables
# Let x be the number of product A to produce, y be the number of product B to produce
x = pulp.LpVariable('x', lowBound=0, cat='Continuous')
y = pulp.LpVariable('y', lowBound=0, cat='Continuous')

# Step 3: Define the objective function
# Objective function: Maximize 20x + 25y
problem += 20*x + 25*y, "Profit"

# Step 4: Define the constraints
# Constraint 1: 4x + 3y <= 100 (Material constraint)
problem += 4*x + 3*y <= 100, "Material"

# Constraint 2: 2x + 5y <= 80 (Labor constraint)
problem += 2*x + 5*y <= 80, "Labor"

# Step 5: Solve the problem
problem.solve()

# Step 6: Display the results
print("Status:", pulp.LpStatus[problem.status])
print("Optimal number of product A to produce:", pulp.value(x))
print("Optimal number of product B to produce:", pulp.value(y))
print("Maximum Profit:", pulp.value(problem.objective))

Explanation of the Code:

Define the Problem: We create a linear programming problem with the objective to maximize profit.
Decision Variables: We define the decision variables x and y for the number of products A and B to produce.
These variables are continuous and non-negative.
Objective Function: We define the objective function to maximize, which is the total profit: 20*x + 25*y.
Constraints: We add constraints based on the available resources:
- Material constraint: 4*x + 3*y <= 100
- Labor constraint: 2*x + 5*y <= 80
Solve the Problem: We solve the linear programming problem using PuLP’s solve method.
Display the Results: We print the optimal values of x and y, and the maximum profit.

Interpretation of Output:

Output

Status: Optimal
Optimal number of product A to produce: 18.571429
Optimal number of product B to produce: 8.5714286
Maximum Profit: 585.714295

Interpretation

The optimal solution is to produce $18.5$ units of product A and $8.57$ units of product B, which yields a maximum profit of $585.71.

This is a basic example of using PuLP in Python.

The library is powerful and can handle more complex constraints, variables, and objectives, including mixed-integer programming.

Altair

August 3, 2024

Altair

Here’s a advanced and useful sample code using Altair, which demonstrates how to create a layered chart with interactivity, such as tooltips and selection, along with custom encoding:

import altair as alt
from vega_datasets import data

# Load dataset
source = data.cars()

# Define a brush selection
brush = alt.selection(type='interval')

# Base chart with a scatter plot
base = alt.Chart(source).mark_point().encode(
    x='Horsepower:Q',
    y='Miles_per_Gallon:Q',
    color=alt.condition(brush, 'Origin:N', alt.value('lightgray')),
    tooltip=['Name:N', 'Origin:N', 'Horsepower:Q', 'Miles_per_Gallon:Q']
).properties(
    width=600,
    height=400
).add_selection(
    brush
)

# Layer a bar chart on top that shows the distribution of 'Origin' for selected points
bars = alt.Chart(source).mark_bar().encode(
    x='count():Q',
    y='Origin:N',
    color='Origin:N'
).transform_filter(
    brush
)

# Combine the scatter plot and the bar chart
chart = base & bars
chart

Result

Explanation

Brush Selection:
An interactive selection that allows users to drag over the scatter plot to select a region.
Scatter Plot:
A basic plot where points are colored by their origin, and the color changes based on the selection.
Bar Chart:
Shows the distribution of car origins, updated based on the selection made in the scatter plot.
Interactivity:
Tooltips provide detailed information when hovering over the points, and the chart updates dynamically based on the user’s selection.

This example combines interactivity with advanced encoding and layout, making it useful for exploratory data analysis.

マルチレイヤーグラフ Altair

August 2, 2024

Altairを使用してPythonで複雑なグラフを作成する例を紹介します。

以下のコードでは、サンプルデータを作成し、複数のデータ系列を組み合わせた複雑な可視化を行います。

この例では、散布図と折れ線グラフを組み合わせたマルチレイヤーグラフを描きます。

import altair as alt
import pandas as pd
import numpy as np

# サンプルデータの作成
np.random.seed(0)
x = np.linspace(0, 100, 100)
y1 = np.sin(x / 10) + np.random.normal(0, 0.1, 100)
y2 = np.cos(x / 10) + np.random.normal(0, 0.1, 100)
category = np.random.choice(['A', 'B', 'C'], 100)

df = pd.DataFrame({'x': x, 'y1': y1, 'y2': y2, 'category': category})

# 散布図の作成
scatter = alt.Chart(df).mark_point().encode(
    x='x',
    y='y1',
    color='category',
    tooltip=['x', 'y1', 'category']
).properties(
    title='Scatter Plot with Line Chart Overlay'
)

# 折れ線グラフの作成
line = alt.Chart(df).mark_line(color='red').encode(
    x='x',
    y='y2'
)

# 散布図と折れ線グラフのオーバーレイ
combined_chart = scatter + line

# グラフの表示
combined_chart.show()

コードの説明

データの作成:
numpyを使用して、$100$点のデータを生成しています。
y1はノイズの入ったサイン波、y2はノイズの入ったコサイン波です。
また、categoryはランダムに割り当てたカテゴリです。
散布図の作成:
alt.Chart(df).mark_point()で散布図を作成し、x, y1, categoryでデータをエンコードしています。
折れ線グラフの作成:
alt.Chart(df).mark_line(color='red')で折れ線グラフを作成し、x, y2でデータをエンコードしています。
グラフのオーバーレイ:
scatter + lineで、散布図と折れ線グラフをオーバーレイして複合グラフを作成しています。
グラフの表示: combined_chart.show()でグラフを表示します。

[実行結果]

Altairを使用することで、簡単に複雑な可視化を作成できます。

このコードをGoogle Colabなどの環境で実行すると、インタラクティブなグラフが表示されます。

グラフ作成 Altair

August 1, 2024

グラフ作成 Altair

Altairを使った簡単なグラフ作成のサンプルコードを紹介します。

Altairは宣言型で直感的にグラフを描けるため、データの視覚化が非常に簡単です。

以下のサンプルでは、Altairを使って基本的な散布図を作成します。

サンプルコード: Altairでの散布図作成

!pip install altair_viewer
import altair as alt
import pandas as pd

# サンプルデータの作成
data = pd.DataFrame({
    'x': [1, 2, 3, 4, 5],
    'y': [2, 3, 7, 6, 9],
    'category': ['A', 'A', 'B', 'B', 'C']
})

# Altairで散布図を作成
chart = alt.Chart(data).mark_circle(size=60).encode(
    x='x',
    y='y',
    color='category',  # カテゴリーに基づく色分け
    tooltip=['x', 'y', 'category']  # ホバー時に表示されるデータ
).properties(
    title='Simple Scatter Plot with Altair'
)

# グラフを表示
chart.show()

コードの解説

ライブラリのインポート:
- import altair as alt: Altairライブラリをインポートします。
- import pandas as pd: サンプルデータを作成するためにPandasをインポートします。
データの作成:
- Pandasのデータフレームとしてデータを作成します。
  このデータにはx軸、y軸の値と、各点のcategory情報が含まれています。
Altairでの散布図作成:
- alt.Chart(data): データフレームをAltairのチャートオブジェクトに変換します。
- mark_circle(size=60): 円形のマークを使用して、各データポイントをプロットします。
- encode: プロットの設定を行います。
  - x='x': x軸の値を設定。
  - y='y': y軸の値を設定。
  - color='category': カテゴリーごとに色を分けます。
  - tooltip=['x', 'y', 'category']: ホバー時に表示するデータの情報を設定します。
- properties: グラフのタイトルやその他のプロパティを設定します。
グラフの表示:
- chart.show(): グラフを表示します（Jupyter Notebookを使っている場合はchartとするだけで表示されます）。

実行結果

このコードを実行すると、カテゴリーごとに色分けされた散布図が表示されます。

各データポイントにマウスをホバーすると、ツールチップに$x$, $y$の値とカテゴリーが表示されます。

Altairは、簡単なコードでインタラクティブなグラフを作成できる強力なツールです。

さらに複雑なグラフも、同様の構文で簡単に作成できます。

最適化問題 SciPy

July 31, 2024

最適化問題 SciPy

SciPyを用いて複雑な最適化問題を解き、その結果をグラフ化する例を示します。

以下では、非線形の目的関数と制約条件を持つ最適化問題を設定します。

問題設定

目的関数:
$$
f(x) = (x_0 - 2)^2 + (x_1 - 3)^2 + \cos(x_0 \cdot x_1)
$$

制約条件:

非線形等式制約:
$$
g_1(x) = x_0^2 + x_1^2 - 4 = 0
$$
非線形不等式制約:
$$
g_2(x) = x_0 \cdot x_1 - 1 \geq 0
$$

初期推定値として$ ( x_0 = 1 ), ( x_1 = 1 ) $を使用します。

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# 目的関数
def objective_function(x):
    return (x[0] - 2)**2 + (x[1] - 3)**2 + np.cos(x[0] * x[1])

# 制約条件
def constraint1(x):
    return x[0]**2 + x[1]**2 - 4

def constraint2(x):
    return x[0] * x[1] - 1

constraints = [{'type': 'eq', 'fun': constraint1},
               {'type': 'ineq', 'fun': constraint2}]

initial_guess = [1, 1]
result = minimize(objective_function, initial_guess, constraints=constraints)

print("最適解:", result.x)
print("最適値:", result.fun)

# 結果のプロット
x0_range = np.linspace(-3, 3, 400)
x1_range = np.linspace(-3, 3, 400)
X0, X1 = np.meshgrid(x0_range, x1_range)
Z = objective_function([X0, X1])

plt.figure(figsize=(10, 6))
contour = plt.contour(X0, X1, Z, levels=50, cmap='viridis')
plt.colorbar(contour)
plt.plot(result.x[0], result.x[1], 'ro', label='Optimal solution')
plt.title('Objective Function Contour and Optimal Solution')
plt.xlabel('$x_0$')
plt.ylabel('$x_1$')
plt.legend()
plt.grid()
plt.show()

このコードでは、SciPyのminimize関数を使って非線形の目的関数と制約条件を持つ最適化問題を解き、その結果をグラフ化しています。

結果として、目的関数の等高線図と最適解をプロットしています。

[実行結果]

ソースコード解説

このソースコードは、非線形制約付きの最適化問題をSciPyを使って解き、その結果を等高線プロットとして視覚化しています。

以下に、コードの各部分を詳しく説明します。

1. インポート

まず、必要なライブラリをインポートします。

1
2
3

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

numpy は数値計算をサポートするライブラリです。
scipy.optimize の minimize 関数は最適化問題を解くために使用します。
matplotlib.pyplot はプロットを作成するためのライブラリです。

2. 目的関数の定義

最適化の対象となる関数を定義します。

1 2	def objective_function(x): return (x[0] - 2)2 + (x[1] - 3)2 + np.cos(x[0] * x[1])

この関数は、2変数$ ( x[0] ) $と$ ( x[1] ) $を持ち、それらの二乗誤差と余弦関数の和として定義されます。

3. 制約条件の定義

2つの制約条件を定義します。

def constraint1(x):
    return x[0]**2 + x[1]**2 - 4

def constraint2(x):
    return x[0] * x[1] - 1

constraint1 は等式制約です。
これは$ ( x[0]^2 + x[1]^2 = 4 ) $という円の方程式に対応します。
constraint2 は不等式制約です。
これは$ ( x[0] \cdot x[1] \geq 1 ) $という条件です。

4. 制約条件の設定

制約条件をリストとして設定します。

1 2	constraints = [{'type': 'eq', 'fun': constraint1}, {'type': 'ineq', 'fun': constraint2}]

minimize 関数に渡すために、制約条件を辞書形式で定義します。

5. 初期推定値の設定

最適化アルゴリズムの初期値を設定します。

1	initial_guess = [1, 1]

最適化の開始点を$ ([1, 1]) $とします。

6. 最適化の実行

minimize 関数を用いて最適化を実行します。

1	result = minimize(objective_function, initial_guess, constraints=constraints)

objective_function を最小化し、与えられた制約条件を満たす解を求めます。

7. 結果の表示

最適解と最適値を出力します。

1 2	print("最適解:", result.x) print("最適値:", result.fun)

8. 結果のプロット

最適化結果を等高線プロットとして視覚化します。

x0_range = np.linspace(-3, 3, 400)
x1_range = np.linspace(-3, 3, 400)
X0, X1 = np.meshgrid(x0_range, x1_range)
Z = objective_function([X0, X1])

等高線プロットのためのグリッドを作成し、目的関数の値を計算します。

9. プロットの作成

等高線プロットを作成し、最適解をプロットします。

plt.figure(figsize=(10, 6))
contour = plt.contour(X0, X1, Z, levels=50, cmap='viridis')
plt.colorbar(contour)
plt.plot(result.x[0], result.x[1], 'ro', label='Optimal solution')
plt.title('Objective Function Contour and Optimal Solution')
plt.xlabel('$x_0$')
plt.ylabel('$x_1$')
plt.legend()
plt.grid()
plt.show()

等高線プロットを描画し、カラーバーを追加します。
最適解を赤い点としてプロットし、ラベルを付けます。
タイトルやラベルを設定し、プロットを表示します。

このコード全体の流れは、非線形最適化問題を定義し、それをSciPyの minimize 関数で解き、その結果を視覚的に確認するための等高線プロットを作成するものです。

結果解説

[実行結果]

このグラフは、非線形最適化問題の解法の結果を示しています。

目的関数の等高線図（Contour plot）と最適解の位置がプロットされています。

問題設定

目的関数:
$$
f(x) = (x_0 - 2)^2 + (x_1 - 3)^2 + \cos(x_0 \cdot x_1)
$$

制約条件:

非線形等式制約:
$$
g_1(x) = x_0^2 + x_1^2 - 4 = 0
$$
非線形不等式制約:
$$
g_2(x) = x_0 \cdot x_1 - 1 \geq 0
$$

初期推定値:
$$
x_0 = 1, x_1 = 1
$$

解の結果

最適解:
$$
x = [1.21825985, 1.5861409]
$$

最適値:
$$
f(x) = 2.265404399189813
$$

グラフの解説

等高線図（Contour plot）:
- $x$軸と$y$軸は、それぞれ変数$ (x_0) $と$ (x_1) $を表しています。
- カラーマップは、目的関数 $ ( f(x) ) $の値を表しています。
  色が濃いほど、目的関数の値が高いことを示します。
- 等高線は、同じ目的関数の値を持つ点を繋いだ線です。
  等高線が密な部分は勾配が急なことを示し、疎な部分は勾配が緩やかなことを示します。
最適解の位置:
- 赤い点は、最適解を示しています。
  この点は、与えられた制約条件を満たしつつ、目的関数の値を最小にする点です。
- グラフ上で赤い点の位置は、近くの等高線と一致しています。
  この点が最適解であることを示しています。

解説

この最適化問題では、2つの非線形制約条件の下で、非線形な目的関数を最小化する最適な解を求めています。
最適解は、与えられた初期推定値から始まり、制約条件を満たす範囲内で目的関数の値を最小化する点として見つかりました。
グラフ上では、この最適解の位置が目的関数の等高線と共に示されており、目的関数の形状や最適解の位置関係が視覚的に確認できます。

これにより、問題の最適化プロセスと結果を直感的に理解することができます。